text:
import GPLinearODEMaker.powers_of_negative_one

"""
    rq_kernel(hyperparameters, δ, dorder; shift_ind=0)

Created by kernel_coder(). Requires 2 hyperparameters.
Likely created using rq_kernel_base() as an input.
Use with include("src/kernels/rq_kernel.jl").

# Arguments
- `hyperparameters::Vector`: The hyperparameter values. For this kernel, they should be `["α", "μ"]`
- `δ::Real`: The difference between the inputs (e.g. `t1 - t2`)
- `dorder::Vector{<:Integer}`: How many times to differentiate with respect to the inputs and the `hyperparameters` (e.g. `dorder=[0, 1, 0, 2]` would correspond to differentiating once w.r.t the second input and twice w.r.t `hyperparameters[2]`)
- `shift_ind::Integer=0`: If changed, the index of which hyperparameter is the `δ` shifting one
"""
function rq_kernel(
    hyperparameters::Vector{<:Real},
    δ::Real,
    dorder::Vector{<:Integer};
    shift_ind::Integer=0)

    include_shift = shift_ind!=0
    @assert length(hyperparameters)==2+Int(include_shift) "hyperparameters is the wrong length"
    dorder_len = 4 + Int(include_shift)
    @assert length(dorder)==dorder_len "dorder is the wrong length"
    @assert maximum(dorder) < 3 "No more than two derivatives for each time or hyperparameter can be calculated"

    dorder2 = dorder[2]
    dorder[2] += dorder[1]

    if include_shift
        dorder[2] += dorder[shift_ind+2]
        δ += hyperparameters[shift_ind]
        hyperparameters_view = view(hyperparameters, 1:3 .!= shift_ind)
        dorder_view = view(dorder, 2:dorder_len)
        dorder_view = view(dorder_view, 1:(dorder_len-1) .!= (shift_ind+1))
    else
        hyperparameters_view = hyperparameters
        dorder_view = view(dorder, 2:dorder_len)
    end

    α = hyperparameters_view[1]
    μ = hyperparameters_view[2]

    if dorder_view==[6, 0, 2]
        func = -90*μ*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 585*δ^2*μ^2*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + (-1935/4)*δ^4*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + (-465/4)*δ^6*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + (-39/4)*δ^8*μ^5*(-6 - α)*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-7 - α)/α^6 + (-1/4)*δ^10*μ^6*(-6 - α)*(-7 - α)*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-8 - α)/α^7
    end

    if dorder_view==[5, 0, 2]
        func = -90*δ*μ*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 165*δ^3*μ^2*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + (-255/4)*δ^5*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + (-15/2)*δ^7*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + (-1/4)*δ^9*μ^5*(-6 - α)*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-7 - α)/α^6
    end

    if dorder_view==[4, 0, 2]
        func = -6*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α - 42*δ^2*μ*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + (-123/4)*δ^4*μ^2*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + (-11/2)*δ^6*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + (-1/4)*δ^8*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5
    end

    if dorder_view==[3, 0, 2]
        func = -6*δ*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α - 12*δ^3*μ*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + (-15/4)*δ^5*μ^2*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + (-1/4)*δ^7*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4
    end

    if dorder_view==[2, 0, 2]
        func = -3*δ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + (-9/4)*δ^4*μ*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + (-1/4)*δ^6*μ^2*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3
    end

    if dorder_view==[1, 0, 2]
        func = -δ^3*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + (-1/4)*δ^5*μ*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2
    end

    if dorder_view==[0, 0, 2]
        func = (-1/4)*δ^4*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α
    end

    if dorder_view==[6, 1, 1]
        func = 45*μ^2*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 45*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 90*μ^2*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 - 45*μ^2*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 15*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-3 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + (375/2)*δ^2*μ^3*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + (375/2)*δ^2*μ^3*(-1 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 180*δ^2*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 555*δ^2*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 + (195/2)*δ^4*μ^4*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + (195/2)*δ^4*μ^4*(-1 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + (195/2)*δ^4*μ^4*(-1 - α)*(-2 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 75*δ^4*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 180*δ^2*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + (-15/2)*δ^2*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 45*δ^2*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-4 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + (735/2)*δ^4*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 + (27/2)*δ^6*μ^5*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + (27/2)*δ^6*μ^5*(-1 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + (27/2)*δ^6*μ^5*(-1 - α)*(-2 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + (27/2)*δ^6*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + 6*δ^6*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 - 75*δ^4*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + (-45/2)*δ^4*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 15*δ^4*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-5 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 60*δ^6*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^6 + (1/2)*δ^8*μ^6*(-6 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-7 - α)/α^6 + (1/2)*δ^8*μ^6*(-6 - α)*(-1 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-7 - α)/α^6 + (1/2)*δ^8*μ^6*(-6 - α)*(-1 - α)*(-2 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-7 - α)/α^6 + (1/2)*δ^8*μ^6*(-6 - α)*(-1 - α)*(-2 - α)*(-3 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-7 - α)/α^6 + (1/2)*δ^8*μ^6*(-6 - α)*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-7 - α)/α^6 - 6*δ^6*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*((-1/2)*δ^2*μ*(-6 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + (-15/2)*δ^6*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 - δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-6 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-6 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + (5/2)*δ^8*μ^6*(-6 - α)*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-7 - α)/α^7 + (-1/2)*δ^8*μ^6*(-6 - α)*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*((-1/2)*δ^2*μ*(-6 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-7 - α)/α^6
    end

    if dorder_view==[5, 1, 1]
        func = 45*δ*μ^2*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 45*δ*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 90*δ*μ^2*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 + (95/2)*δ^3*μ^3*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + (95/2)*δ^3*μ^3*(-1 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 40*δ^3*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 45*δ*μ^2*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 15*δ*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-3 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 135*δ^3*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 + 10*δ^5*μ^4*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 10*δ^5*μ^4*(-1 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 10*δ^5*μ^4*(-1 - α)*(-2 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 5*δ^5*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 40*δ^3*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + (-15/2)*δ^3*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 10*δ^3*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-4 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 35*δ^5*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 + (1/2)*δ^7*μ^5*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + (1/2)*δ^7*μ^5*(-1 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + (1/2)*δ^7*μ^5*(-1 - α)*(-2 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + (1/2)*δ^7*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 - 5*δ^5*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 5*δ^5*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - δ^5*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-5 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 2*δ^7*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^6 + (-1/2)*δ^7*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5
    end

    if dorder_view==[4, 1, 1]
        func = 6*μ*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + 6*μ*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^2 - 6*μ*(-1 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α - 3*μ^2*(-1 - α)*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-2 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + (39/2)*δ^2*μ^2*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 18*δ^2*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + (75/2)*δ^2*μ^2*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 + 7*δ^4*μ^3*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 7*δ^4*μ^3*(-1 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 4*δ^4*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 18*δ^2*μ^2*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + (-3/2)*δ^2*μ^2*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 6*δ^2*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-3 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 18*δ^4*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 + (1/2)*δ^6*μ^4*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + (1/2)*δ^6*μ^4*(-1 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + (1/2)*δ^6*μ^4*(-1 - α)*(-2 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 4*δ^4*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 3*δ^4*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - δ^4*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-4 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + (3/2)*δ^6*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 + (-1/2)*δ^6*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4
    end

    if dorder_view==[3, 1, 1]
        func = 6*δ*μ*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + 6*δ*μ*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^2 + (9/2)*δ^3*μ^2*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 3*δ^3*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 6*δ*μ*(-1 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α - 3*δ*μ^2*(-1 - α)*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-2 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + (15/2)*δ^3*μ^2*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 + (1/2)*δ^5*μ^3*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + (1/2)*δ^5*μ^3*(-1 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 3*δ^3*μ^2*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + (-3/2)*δ^3*μ^2*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - δ^3*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-3 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + δ^5*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 + (-1/2)*δ^5*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3
    end

    if dorder_view==[2, 1, 1]
        func = -((-1/2)*δ^2*μ*(-1 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-1 - α) - μ*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-1 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-1 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-1 - α) + 2*δ^2*μ*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + 2*δ^2*μ*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^2 + (1/2)*δ^4*μ^2*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 2*δ^2*μ*(-1 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + (-1/2)*δ^2*μ*(-1 - α)*((-1/2)*δ^2*μ*(-1 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α - δ^2*μ^2*(-1 - α)*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-2 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + (1/2)*δ^4*μ^2*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 + (-1/2)*δ^4*μ^2*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2
    end

    if dorder_view==[1, 1, 1]
        func = -δ*((-1/2)*δ^2*μ*(-1 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-1 - α) - δ*μ*((-1/2)*δ^2/(α*(1 + (1/2)*δ^2*μ/α)) + (-1/2)*δ^2*(-1 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) + (1/4)*δ^4*μ*(-1 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)^2))*(1 + (1/2)*δ^2*μ/α)^(-1 - α) + (-1/2)*δ^3*μ*(-1 - α)*((-1/2)*δ^2*μ*(-1 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α
    end

    if dorder_view==[0, 1, 1]
        func = (-1/2)*δ^2*(1 + (1/2)*δ^2*μ/α)^(-1 - α)*((1/2)*δ^2*μ/(α*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^2
    end

    if dorder_view==[6, 0, 1]
        func = -45*μ^2*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + (-375/2)*δ^2*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + (-195/2)*δ^4*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + (-27/2)*δ^6*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + (-1/2)*δ^8*μ^6*(-6 - α)*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-7 - α)/α^6
    end

    if dorder_view==[5, 0, 1]
        func = -45*δ*μ^2*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + (-95/2)*δ^3*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 10*δ^5*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + (-1/2)*δ^7*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5
    end

    if dorder_view==[4, 0, 1]
        func = -6*μ*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + (-39/2)*δ^2*μ^2*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 7*δ^4*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + (-1/2)*δ^6*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4
    end

    if dorder_view==[3, 0, 1]
        func = -6*δ*μ*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + (-9/2)*δ^3*μ^2*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + (-1/2)*δ^5*μ^3*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3
    end

    if dorder_view==[2, 0, 1]
        func = (-5/2)*δ^2*μ*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + (-1/2)*δ^4*μ^2*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - (1 + (1/2)*δ^2*μ/α)^(-1 - α)
    end

    if dorder_view==[1, 0, 1]
        func = -δ*(1 + (1/2)*δ^2*μ/α)^(-1 - α) + (-1/2)*δ^3*μ*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α
    end

    if dorder_view==[0, 0, 1]
        func = (-1/2)*δ^2*(1 + (1/2)*δ^2*μ/α)^(-1 - α)
    end

    if dorder_view==[6, 2, 0]
        func = -30*μ^3*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 60*μ^3*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 - 60*μ^3*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 - 90*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^4 + 30*μ^3*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 30*μ^3*(-1 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 90*δ^2*μ^4*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 90*δ^2*μ^4*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 90*δ^2*μ^4*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 60*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 - 15*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-3 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-3 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2))/α^2 - 15*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 270*δ^2*μ^4*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 - 270*δ^2*μ^4*(-1 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 - 270*δ^2*μ^4*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 - 30*δ^4*μ^5*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 30*δ^4*μ^5*(-2 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 30*δ^4*μ^5*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 30*δ^4*μ^5*(-1 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 30*δ^4*μ^5*(-1 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 30*δ^4*μ^5*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 540*δ^2*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^5 + 90*δ^2*μ^4*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 90*δ^2*μ^4*(-1 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 90*δ^2*μ^4*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 120*δ^4*μ^5*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 - 120*δ^4*μ^5*(-1 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 - 120*δ^4*μ^5*(-1 - α)*(-2 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 - 120*δ^4*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 - 2*δ^6*μ^6*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 - 2*δ^6*μ^6*(-2 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 - 2*δ^6*μ^6*(-2 - α)*(-3 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 - 2*δ^6*μ^6*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 - 2*δ^6*μ^6*(-1 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 - 2*δ^6*μ^6*(-1 - α)*(-3 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 - 2*δ^6*μ^6*(-1 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 - 2*δ^6*μ^6*(-1 - α)*(-2 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 - 2*δ^6*μ^6*(-1 - α)*(-2 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 - 2*δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + 270*δ^2*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 - 45*δ^2*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-4 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-4 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2))/α^3 - 45*δ^2*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 300*δ^4*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^6 + 30*δ^4*μ^5*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 30*δ^4*μ^5*(-1 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 30*δ^4*μ^5*(-1 - α)*(-2 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 30*δ^4*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 10*δ^6*μ^6*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^6 - 10*δ^6*μ^6*(-1 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^6 - 10*δ^6*μ^6*(-1 - α)*(-2 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^6 - 10*δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^6 - 10*δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^6 + 120*δ^4*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 - 15*δ^4*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-5 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-5 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2))/α^4 - 15*δ^4*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 30*δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^7 + 2*δ^6*μ^6*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*((-1/2)*δ^2*μ*(-6 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + 2*δ^6*μ^6*(-1 - α)*(-3 - α)*(-4 - α)*(-5 - α)*((-1/2)*δ^2*μ*(-6 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + 2*δ^6*μ^6*(-1 - α)*(-2 - α)*(-4 - α)*(-5 - α)*((-1/2)*δ^2*μ*(-6 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + 2*δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(-5 - α)*((-1/2)*δ^2*μ*(-6 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + 2*δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-6 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + 10*δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*((-1/2)*δ^2*μ*(-6 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^6 - δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-6 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-6 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2))/α^5 - δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*((-1/2)*δ^2*μ*(-6 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5
    end

    if dorder_view==[5, 2, 0]
        func = -30*δ*μ^3*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 60*δ*μ^3*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 - 60*δ*μ^3*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 - 20*δ^3*μ^4*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 20*δ^3*μ^4*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 20*δ^3*μ^4*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 90*δ*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^4 + 30*δ*μ^3*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 30*δ*μ^3*(-1 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 60*δ^3*μ^4*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 - 60*δ^3*μ^4*(-1 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 - 60*δ^3*μ^4*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 - 2*δ^5*μ^5*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 2*δ^5*μ^5*(-2 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 2*δ^5*μ^5*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 2*δ^5*μ^5*(-1 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 2*δ^5*μ^5*(-1 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 2*δ^5*μ^5*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 60*δ*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 - 15*δ*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-3 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-3 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2))/α^2 - 15*δ*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 120*δ^3*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^5 + 20*δ^3*μ^4*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 20*δ^3*μ^4*(-1 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 20*δ^3*μ^4*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 8*δ^5*μ^5*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 - 8*δ^5*μ^5*(-1 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 - 8*δ^5*μ^5*(-1 - α)*(-2 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 - 8*δ^5*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 + 60*δ^3*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 - 10*δ^3*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-4 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-4 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2))/α^3 - 10*δ^3*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 20*δ^5*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^6 + 2*δ^5*μ^5*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 2*δ^5*μ^5*(-1 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 2*δ^5*μ^5*(-1 - α)*(-2 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 2*δ^5*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 8*δ^5*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 - δ^5*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-5 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-5 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2))/α^4 - δ^5*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4
    end

    if dorder_view==[4, 2, 0]
        func = -6*μ^2*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^2 - 6*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^3 + 6*μ^2*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α - 12*δ^2*μ^3*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 6*μ^2*(-1 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^2 - 3*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-2 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-2 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2))/α - 3*μ^2*(-1 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α - 24*δ^2*μ^3*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 - 24*δ^2*μ^3*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 - 2*δ^4*μ^4*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 2*δ^4*μ^4*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 2*δ^4*μ^4*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 36*δ^2*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^4 + 12*δ^2*μ^3*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 12*δ^2*μ^3*(-1 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 6*δ^4*μ^4*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 - 6*δ^4*μ^4*(-1 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 - 6*δ^4*μ^4*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 + 24*δ^2*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 - 6*δ^2*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-3 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-3 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2))/α^2 - 6*δ^2*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 12*δ^4*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^5 + 2*δ^4*μ^4*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 2*δ^4*μ^4*(-1 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 2*δ^4*μ^4*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 6*δ^4*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 - δ^4*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-4 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-4 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2))/α^3 - δ^4*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3
    end

    if dorder_view==[3, 2, 0]
        func = -6*δ*μ^2*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^2 - 2*δ^3*μ^3*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 6*δ*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^3 + 6*δ*μ^2*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α - 4*δ^3*μ^3*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 - 4*δ^3*μ^3*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 + 6*δ*μ^2*(-1 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^2 - 3*δ*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-2 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-2 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2))/α - 3*δ*μ^2*(-1 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α - 6*δ^3*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^4 + 2*δ^3*μ^3*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 2*δ^3*μ^3*(-1 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 4*δ^3*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 - δ^3*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-3 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-3 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2))/α^2 - δ^3*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2
    end

    if dorder_view==[2, 2, 0]
        func = -μ*(1 + (1/2)*δ^2*μ/α)^(-1 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-1 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-1 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2)) - μ*((-1/2)*δ^2*μ*(-1 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-1 - α) - 2*δ^2*μ^2*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^2 - 2*δ^2*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^3 + 2*δ^2*μ^2*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + 2*δ^2*μ^2*(-1 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^2 - δ^2*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-2 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-2 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2))/α - δ^2*μ^2*(-1 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α
    end

    if dorder_view==[1, 2, 0]
        func = -δ*μ*(1 + (1/2)*δ^2*μ/α)^(-1 - α)*(δ^2*μ/(α^2*(1 + (1/2)*δ^2*μ/α)) + δ^2*μ*(-1 - α)/(α^3*(1 + (1/2)*δ^2*μ/α)) + (-1/4)*δ^4*μ^2*(-1 - α)/(α^4*(1 + (1/2)*δ^2*μ/α)^2)) - δ*μ*((-1/2)*δ^2*μ*(-1 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2*(1 + (1/2)*δ^2*μ/α)^(-1 - α)
    end

    if dorder_view==[0, 2, 0]
        func = (1 + (1/2)*δ^2*μ/α)^(-α)*((1/2)*δ^2*μ/(α*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))^2 + (1/4)*δ^4*μ^2*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^3
    end

    if dorder_view==[6, 1, 0]
        func = 15*μ^3*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 15*μ^3*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 30*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 - 15*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 45*δ^2*μ^4*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 45*δ^2*μ^4*(-1 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 45*δ^2*μ^4*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 135*δ^2*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 + 15*δ^4*μ^5*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 15*δ^4*μ^5*(-1 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 15*δ^4*μ^5*(-1 - α)*(-2 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 15*δ^4*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 45*δ^2*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 60*δ^4*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 + δ^6*μ^6*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + δ^6*μ^6*(-1 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + δ^6*μ^6*(-1 - α)*(-2 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 + δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5 - 15*δ^4*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + 5*δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^6 - δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*((-1/2)*δ^2*μ*(-6 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5
    end

    if dorder_view==[5, 1, 0]
        func = 15*δ*μ^3*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 15*δ*μ^3*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 30*δ*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 + 10*δ^3*μ^4*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 10*δ^3*μ^4*(-1 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 10*δ^3*μ^4*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 15*δ*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 30*δ^3*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 + δ^5*μ^5*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + δ^5*μ^5*(-1 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + δ^5*μ^5*(-1 - α)*(-2 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 + δ^5*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - 10*δ^3*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + 4*δ^5*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^5 - δ^5*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*((-1/2)*δ^2*μ*(-5 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4
    end

    if dorder_view==[4, 1, 0]
        func = 3*μ^2*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + 3*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^2 - 3*μ^2*(-1 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + 6*δ^2*μ^3*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 6*δ^2*μ^3*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 12*δ^2*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 + δ^4*μ^4*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + δ^4*μ^4*(-1 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 + δ^4*μ^4*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 6*δ^2*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + 3*δ^4*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^4 - δ^4*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*((-1/2)*δ^2*μ*(-4 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3
    end

    if dorder_view==[3, 1, 0]
        func = 3*δ*μ^2*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + 3*δ*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^2 + δ^3*μ^3*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 + δ^3*μ^3*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 3*δ*μ^2*(-1 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + 2*δ^3*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^3 - δ^3*μ^3*(-1 - α)*(-2 - α)*((-1/2)*δ^2*μ*(-3 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2
    end

    if dorder_view==[2, 1, 0]
        func = -μ*((-1/2)*δ^2*μ*(-1 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-1 - α) + δ^2*μ^2*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α + δ^2*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α^2 - δ^2*μ^2*(-1 - α)*((-1/2)*δ^2*μ*(-2 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α
    end

    if dorder_view==[1, 1, 0]
        func = -δ*μ*((-1/2)*δ^2*μ*(-1 - α)/(α^2*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))*(1 + (1/2)*δ^2*μ/α)^(-1 - α)
    end

    if dorder_view==[0, 1, 0]
        func = (1 + (1/2)*δ^2*μ/α)^(-α)*((1/2)*δ^2*μ/(α*(1 + (1/2)*δ^2*μ/α)) - log(1 + (1/2)*δ^2*μ/α))
    end

    if dorder_view==[6, 0, 0]
        func = -15*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 45*δ^2*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - 15*δ^4*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4 - δ^6*μ^6*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(-5 - α)*(1 + (1/2)*δ^2*μ/α)^(-6 - α)/α^5
    end

    if dorder_view==[5, 0, 0]
        func = -15*δ*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - 10*δ^3*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3 - δ^5*μ^5*(-1 - α)*(-2 - α)*(-3 - α)*(-4 - α)*(1 + (1/2)*δ^2*μ/α)^(-5 - α)/α^4
    end

    if dorder_view==[4, 0, 0]
        func = -3*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α - 6*δ^2*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2 - δ^4*μ^4*(-1 - α)*(-2 - α)*(-3 - α)*(1 + (1/2)*δ^2*μ/α)^(-4 - α)/α^3
    end

    if dorder_view==[3, 0, 0]
        func = -3*δ*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α - δ^3*μ^3*(-1 - α)*(-2 - α)*(1 + (1/2)*δ^2*μ/α)^(-3 - α)/α^2
    end

    if dorder_view==[2, 0, 0]
        func = -μ*(1 + (1/2)*δ^2*μ/α)^(-1 - α) - δ^2*μ^2*(-1 - α)*(1 + (1/2)*δ^2*μ/α)^(-2 - α)/α
    end

    if dorder_view==[1, 0, 0]
        func = -δ*μ*(1 + (1/2)*δ^2*μ/α)^(-1 - α)
    end

    if dorder_view==[0, 0, 0]
        func = (1 + (1/2)*δ^2*μ/α)^(-α)
    end

    dorder[2] = dorder2  # resetting dorder[2]
    return powers_of_negative_one(dorder2) * float(func)  # correcting for amount of t2 derivatives

end

return rq_kernel, 2  # the function handle and the number of kernel hyperparameters

max_stars_repo_name: christiangil/GPLinearODEMaker

max_issues_repo_name: christiangil/GPLinearODEMaker

max_forks_repo_name: christiangil/GPLinearODEMaker


text:
using Documenter, KNITRO

makedocs(
    modules = [KNITRO]
)

deploydocs(
    deps = Deps.pip("pygments", "mkdocs", "mkdocs-material", "python-markdown-math"),
    repo = "github.com/jump-dev/KNITRO.jl.git",
    julia = "0.4"
)

ext: jl


text:
mutable struct DESPOTBeliefUpdater{S,A,O} <: POMDPs.Updater
    pomdp::POMDP
    num_updates::Int64
    rng::AbstractRNG
    seed::UInt32
    rand_max::Int64
    belief_update_seed::UInt32
    particle_weight_threshold::Float64
    eff_particle_fraction::Float64

    #pre-allocated variables (TODO: add the rest at some point)
    n_particles::Int64
    next_state::S
    observation::O
    new_particle::DESPOTParticle{S}
    n_sampled::Int64
    obs_probability::Float64

    #default constructor
    function DESPOTBeliefUpdater{S,A,O}(pomdp::POMDP{S,A,O};
                                 seed::UInt32 = convert(UInt32, 42),
                                 rand_max::Int64 = 2147483647,
                                 n_particles = 500,
                                 particle_weight_threshold::Float64 = 1e-20,
                                 eff_particle_fraction::Float64 = 0.05,
                                 next_state::S = S(),
                                 observation::O = O(),
                                 rng::AbstractRNG=DESPOTDefaultRNG(convert(UInt32, seed ⊻ (n_particles+1)), rand_max)
                                ) where {S,A,O}
        this = new()
        this.pomdp = pomdp
        this.num_updates = 0
        this.belief_update_seed = seed ⊻ (n_particles + 1)
        this.rng = rng
        this.rand_max = rand_max
        this.particle_weight_threshold = particle_weight_threshold
        this.eff_particle_fraction = eff_particle_fraction
        this.n_particles = n_particles

        # init preallocated variables
        this.next_state = next_state
        this.observation = observation
        this.new_particle = DESPOTParticle{S}(this.next_state, 1, 1) #placeholder
        this.n_sampled = 0
        this.obs_probability = -1.0
        return this
    end
end

# Special create_belief version for DESPOTBeliefUpdater
create_belief{S,A,O}(bu::DESPOTBeliefUpdater{S,A,O}) =
    DESPOTBelief{S,A,O}(Array{DESPOTParticle{S}}(bu.n_particles), History{A,O}())

get_belief_update_seed(bu::DESPOTBeliefUpdater) = bu.seed ⊻ (bu.n_particles + 1)

reset_belief(bu::DESPOTBeliefUpdater) = bu.num_updates = 0

function initialize_belief{S,A,O}(bu::DESPOTBeliefUpdater{S,A,O},
                  state_distribution::ParticleDistribution{S})

    new_belief = create_belief(bu)
    n_particles = length(state_distribution.particles)

    # convert to DESPOTParticle type
    pool = Array{DESPOTParticle{S}}(n_particles)

    for i in 1:n_particles
        pool[i] = DESPOTParticle{S}(state_distribution.particles[i].state,
                                 i, # id
                                 state_distribution.particles[i].weight)
    end

    DESPOT.sample_particles!(new_belief.particles,
                             pool,
                             bu.n_particles,
                             bu.belief_update_seed,
                             bu.rand_max)

    #shuffle!(new_belief.particles) #TODO: uncomment if higher randomness is required
    return new_belief
end

#XXX this could probably be combined with the method above or made more efficient without the sample_particles!
function initialize_belief{S}(bu::DESPOTBeliefUpdater{S}, b)
    new_belief = create_belief(bu)

    pool = Array{DESPOTParticle{S}}(bu.n_particles)
    w = 1.0/bu.n_particles
    for i in 1:bu.n_particles
        pool[i] = DESPOTParticle{S}(rand(bu.rng, b), i, w)
    end

    DESPOT.sample_particles!(new_belief.particles,
                             pool,
                             bu.n_particles,
                             bu.belief_update_seed,
                             bu.rand_max)

    return new_belief
end

updater{S,A,O}(p::DESPOTPolicy{S,A,O}) = DESPOTBeliefUpdater{S,A,O}(p.pomdp,
                                                                    n_particles=p.solver.config.n_particles,
                                                                    next_state=p.solver.next_state,
                                                                    observation=p.solver.curr_obs,
                                                                    rng=p.solver.rng
                                                                   )

function normalize!{S}(particles::Vector{DESPOTParticle{S}})
    prob_sum = 0.0
    for p in particles
        prob_sum += p.weight
    end
    for p in particles
        p.weight /= prob_sum
    end
end

function update{S,A,O}(bu::DESPOTBeliefUpdater{S,A,O},
                current_belief::DESPOTBelief{S},
                action::A,
                obs::O)

    updated_belief::DESPOTBelief{S} = create_belief(bu)
    random_number = Array{Float64}(1)

    if bu.n_particles != length(current_belief.particles)
        err("belief size mismatch: belief updater - $(bu.n_particles) particles, belief - $(length(current_belief.particles))")
    end
    updated_belief.particles = []

    #reset RNG
    if isa(bu.rng, DESPOTDefaultRNG)
        bu.rng = DESPOTDefaultRNG(bu.belief_update_seed, bu.rand_max)
    end

    for p in current_belief.particles
        rand!(bu.rng, random_number)

        if isa(bu.rng, DESPOTDefaultRNG)
            rng = DESPOTRandomNumber(random_number[1])
        else
            rng = bu.rng
        end

        bu.next_state = POMDPs.generate_s(bu.pomdp, p.state, action, rng)

        #get observation distribution for (s,a,s') tuple
        od = POMDPs.observation(bu.pomdp, p.state, action, bu.next_state)

        bu.obs_probability = pdf(od, obs)

        if bu.obs_probability > 0.0
            bu.new_particle = DESPOTParticle(bu.next_state, p.id, p.weight * bu.obs_probability)
            push!(updated_belief.particles, bu.new_particle)
        end
    end

    normalize!(updated_belief.particles)

    if length(updated_belief.particles) == 0
        # No resulting state is consistent with the given observation, so create
        # states randomly until we have enough that are consistent.
        warn("Particle filter empty. Bootstrapping with random states")
        bu.n_sampled = 0
        resample_rng = DESPOTDefaultRNG(bu.belief_update_seed, bu.rand_max)
        particle_number::Int64 = 0
        while bu.n_sampled < bu.n_particles
            #TODO: see if this can be done better
            #Pick a random particle from the current belief state as the initial state
            rand!(resample_rng, random_number)
            particle_number = ceil(bu.n_particles * random_number[1])

            next_state = POMDPs.rand(resample_rng, states(bu.pomdp)) #generate a random state
            od = POMDPs.observation(bu.pomdp,
                               current_belief.particles[particle_number].state,
                               action,
                               next_state)
            bu.obs_probability = pdf(od, bu.observation)
            if bu.obs_probability > 0.0
                bu.n_sampled += 1
                bu.new_particle = DESPOTParticle(next_state, bu.obs_probability)
                push!(updated_belief.particles, bu.new_particle)
            end
        end
        normalize!(updated_belief.particles)
        return updated_belief.particles
    end

    # Remove all particles below the threshold weight
    viable_particle_indices = Array{Int64}(0)
    for i in 1:length(updated_belief.particles)
        if updated_belief.particles[i].weight >= bu.particle_weight_threshold
            push!(viable_particle_indices, i)
        end
    end
    updated_belief.particles = updated_belief.particles[viable_particle_indices]

    if length(updated_belief.particles) != 0
        normalize!(updated_belief.particles)
    end

    # Resample if we have < N particles or number of effective particles drops
    # below the threshold
    num_eff_particles = 0
    for p in updated_belief.particles
        num_eff_particles += p.weight^2
    end

    num_eff_particles = 1./num_eff_particles
    if (num_eff_particles < bu.n_particles * bu.eff_particle_fraction) ||
        (length(updated_belief.particles) < bu.n_particles)
        resampled_set = Array{DESPOTParticle{S}}(bu.n_particles)
        sample_particles!(resampled_set,
                          updated_belief.particles,
                          bu.n_particles,
                          bu.belief_update_seed,
                          bu.rand_max)
        updated_belief.particles = resampled_set
    end

    # Finally, update history
    add(updated_belief.history, action, obs)

    return updated_belief
end


text:
__precompile__(true)
module InterestRates

# 0.4 compat
if !isdefined(Core, :String)
    typealias String UTF8String
end

using BusinessDays
using Base.Dates

export
	DayCountConvention,
	CompoundingType,
	CurveMethod,
	AbstractIRCurve,
	ERF, ER, discountfactor, zero_rate, forward_rate,
	ERF_to_rate, discountfactor_to_rate,
	isnullcurve

include("types.jl")
include("nullcurve.jl")

############# DAYCOUNT #################

daycount(conv::BDays252, date_start::Date, date_end::Date) = Int(bdays(conv.hc, date_start, date_end))
daycount(::Actual360, date_start::Date, date_end::Date) = Int(date_end - date_start)
daycount(::Actual365, date_start::Date, date_end::Date) = Int(date_end - date_start)

advancedays(conv::BDays252, date_start::Date, daycount::Int) = advancebdays(conv.hc, date_start, daycount)
advancedays(::Actual360, date_start::Date, daycount::Int) = date_start + Day(daycount)
advancedays(::Actual365, date_start::Date, daycount::Int) = date_start + Day(daycount)

function advancedays(conv::DayCountConvention, date_start::Date, daycount_vec::Vector{Int})
	l = length(daycount_vec)
	result = Array(Date, l)
	for i in 1:l
		result[i] = advancedays(conv, date_start, daycount_vec[i])
	end
	return result
end

advancedays(curve::AbstractIRCurve, daycount) = advancedays(curve_get_daycount(curve), curve_get_date(curve), daycount)

daysperyear(::BDays252) = 252
daysperyear(::Actual360) = 360
daysperyear(::Actual365) = 365

yearfraction(conv::DayCountConvention, date_start::Date, date_end::Date) = daycount(conv, date_start, date_end) / daysperyear(conv)
yearfraction(curve::AbstractIRCurve, maturity::Date) = yearfraction(curve_get_daycount(curve), curve_get_date(curve), maturity)

######### COMPOUNDING TYPE ############

# Effective Rate Factor
_ERF(::ContinuousCompounding, r::Float64, t::Float64) = t == 0.0 ? 1.0 : exp(r*t)
_ERF(::SimpleCompounding, r::Float64, t::Float64) = t == 0.0 ? 1.0 : 1.0 + r*t
_ERF(::ExponentialCompounding, r::Float64, t::Float64) = t == 0.0 ? 1.0 : (1.0+r)^t
_ERF(ct::CompoundingType, dcc::DayCountConvention, r::Float64, date_start::Date, date_end::Date) = _ERF(ct, r, yearfraction(dcc, date_start, date_end))
ERF(curve::AbstractIRCurve, maturity::Date) = _ERF(curve_get_compounding(curve), curve_get_daycount(curve), zero_rate(curve, maturity), curve_get_date(curve), maturity)
ERF(curve::AbstractIRCurve, forward_date::Date, maturity::Date) = ERF(curve, maturity) / ERF(curve, forward_date)

_ERF_to_rate(::ContinuousCompounding, ERF::Float64, t::Float64) = log(ERF) / t
_ERF_to_rate(::SimpleCompounding, ERF::Float64, t::Float64) = (ERF-1.0) / t
_ERF_to_rate(::ExponentialCompounding, ERF::Float64, t::Float64) = ERF^(1.0/t) - 1.0
ERF_to_rate(curve::AbstractIRCurve, ERF::Float64, t::Float64) = _ERF_to_rate(curve_get_compounding(curve), ERF, t)

discountfactor_to_rate(c::CompoundingType, _discountfactor_::Float64, t::Float64) = _ERF_to_rate(c, 1.0 / _discountfactor_, t)

function discountfactor_to_rate(c::CompoundingType, _discountfactor_vec_::Vector{Float64}, t_vec::Vector{Float64})
	l = length(_discountfactor_vec_)

	if l != length(t_vec)
		error("_discountfactor_vec_ and t_vec must have the same length. ($l != $(length(t_vec)))")
	end

	result = Array(Float64, l)
	for i in 1:l
		result[i] = discountfactor_to_rate(c, _discountfactor_vec_[i], t_vec[i])
	end

	return result
end

# Effective Rate = [Effective Rate Factor] - 1
_ER(c::CompoundingType, r::Float64, t::Float64) = _ERF(c, r, t) - 1.0
_ER(ct::CompoundingType, dcc::DayCountConvention, r::Float64, date_start::Date, date_end::Date) = _ER(ct, r, yearfraction(dcc, date_start, date_end))
ER(curve::AbstractIRCurve, maturity::Date) = _ER(curve_get_compounding(curve), curve_get_daycount(curve), zero_rate(curve, maturity), curve_get_date(curve), maturity)
ER(curve::AbstractIRCurve, forward_date::Date, maturity::Date) =  ERF(curve, forward_date, maturity) - 1.0

# [Discount Factor] = 1 / [Effective Rate Factor]
_discountfactor(ct::CompoundingType, r::Float64, t::Float64) = 1.0 / _ERF(ct, r, t)
_discountfactor(ct::CompoundingType, dcc::DayCountConvention, r::Float64, date_start::Date, date_end::Date) = _discountfactor(ct, r, yearfraction(dcc, date_start, date_end))
discountfactor(curve::AbstractIRCurve, maturity::Date) = _discountfactor(curve_get_compounding(curve), curve_get_daycount(curve), zero_rate(curve, maturity), curve_get_date(curve), maturity)

# Vector function for ERF and discountfactor functions
for fun in (:ERF, :discountfactor)
	@eval begin
		function ($fun)(curve::AbstractIRCurve, maturity_vec::Vector{Date})
			l = length(maturity_vec)
			_zero_rate_vec_ = zero_rate(curve, maturity_vec)
			result = Array(Float64, l)
			for i = 1:l
				result[i] = $( Symbol("_" * string(fun)) )(curve_get_compounding(curve), curve_get_daycount(curve), _zero_rate_vec_[i], curve_get_date(curve), maturity_vec[i])
			end
			return result
		end
	end
end

forward_rate(curve::AbstractIRCurve, forward_date::Date, maturity::Date) = ERF_to_rate(curve, ERF(curve, forward_date, maturity), yearfraction(curve_get_daycount(curve), forward_date, maturity))

function days_to_maturity(curve::AbstractIRCurve, maturity::Date)
	const d = daycount(curve_get_daycount(curve), curve_get_date(curve), maturity)
	if d < 0
		error("Maturity date $(maturity) should be greater than curve observation date $(curve_get_date(curve))")
	end
	return d
end

# Let's use the curve's method to multiple-dispatch. Ugly methods _zero_rate are not exported.
zero_rate(curve::AbstractIRCurve, maturity::Date) = _zero_rate(curve_get_method(curve), curve, maturity)
zero_rate(curve::AbstractIRCurve, maturity_vec::Vector{Date}) = _zero_rate(curve_get_method(curve), curve, maturity_vec)

include("methods.jl") # implements various _zero_rate() methods for each CurveMethod

end # module InterestRates


text:
# We want to generate a geo file for curve surrounded by a rectangle
# defined by shift to the bounding box. If the box is tight some of the
# the curve's points will be on the boundary and we need to help gmsh
# to correctly draw the rectanngle

"""Distance between interior point p and box"""
function distance{D, T}(b::BoundingBox{D, T}, p::SVector{D, T})
    lx, hx = b.low[1], b.high[1]
    (lx <= p[1] <= hx) || throw("Point not inside")

    ly, hy = b.low[2], b.high[2]
    (ly <= p[2] <= hy) || throw("Point not inside")

    min(p[1] - lx, hx - p[1], p[2] - ly, hy - p[2])
end

"""Order points such that when visited one after other a loop is formed"""
function loop{T}(points::Vector{SVector{2, T}})
    # The idea here is to order according to angle - think polar coordinates
    x = map(first, points)
    x0 = mean(x)

    y = map(last, points)
    y0 = mean(y)

    angles = atan2.(x, y)
    loop = sortperm(angles)
    # Adding the last guy to convenienly makes lines with zip
    push!(loop, first(loop))
end

"""Generate a geo file for code contained in a padded bounding box"""
function gmsh_code{T}(c::Curve{2, T},
                      xpadding::Tuple{Number, Number},
                      ypadding::Tuple{Number, Number},
                      path::AbstractString)
    # Nongetative
    @assert first(xpadding) >= zero(T) && last(xpadding) >= zero(T)
    @assert first(ypadding) >= zero(T) && last(ypadding) >= zero(T)

    box = BoundingBox(c)
    low, high = box.low, box.high
    # Enclosing frame for the curve; shift to bbox
    frame_pts = [low - SVector(first(xpadding), first(ypadding)),
                 SVector(high[1], low[2]) + SVector(last(xpadding), -first(ypadding)),
                 box.high + SVector(last(xpadding), last(ypadding)),
                 SVector(low[1], high[2]) + SVector(-first(xpadding), last(ypadding))]

    # With points on a boundary we must think how to define the boundaing
    # rectangle in gmsh
    if !(all(xpadding .> 0) && all(ypadding .> 0))
        box = BoundingBox(frame_pts[1], frame_pts[3])
        # Now let's get all the points that are on the boundary
        d = p -> distance(box, p)
        distances = map(d, c.points)
        indices = sortperm(distances)
        # Visiting indices (ordered by distance)
        bdry_indices = Vector{Int}()
        for i in indices
            # Since we have ordering the first guys is not on bondary
            # is also the last one checked
            distances[i] > 0 && break
            # On boundary
            push!(bdry_indices, i)
        end

        bdry_pts = c.points[bdry_indices]
        # It might be that the frame points are already among the boundary ones
        # Then we wat to remove them
        for i in 1:length(frame_pts)
            last(frame_pts) ∈ bdry_pts && pop!(frame_pts)
        end

        # Using bdry and frame points we wamt to draw a loop
        loop_pts = Vector{SVector{2, T}}([bdry_pts..., frame_pts...])
        # Get the loop refered to by coordinates of loop_pts (local)
        loop_idx = loop(loop_pts)
        nbdry, nall = length(bdry_pts), length(c.points)
        # Translate to coords of curve points
        # If local refers to bdry point do a lookup for the global
        # otherwise the numbering is started after the curve points
        loop_idx = [(p <= nbdry) ? bdry_indices[p] : (nall + p - nbdry)
                    for p in loop_idx]
        # And the loop
        bdry_lines = collect(zip(loop_idx[1:end-1], loop_idx[2:end]))
    else
        n = length(c.points)
        bdry_lines = [(1+n, 2+n), (2+n, 3+n), (3+n, 4+n), (4+n, 1+n)]
    end
    # Done with prep work"""
    gmsh_code(c, frame_pts, bdry_lines, path)
end

"""Write the gmsh code"""
function gmsh_code{T}(c::Curve,
                      bdry_pts::Vector{SVector{2, T}},
                      bdry_lines::Vector{Tuple{Int, Int}},
                      path::AbstractString)

    f = open(path, "w")
    # Allow for different sizes on bdry and curve
    # NOTE: call gmsh -setnumber size VALUE0 -setnumber fsize VALUE1 ...
    write(f, "DefineConstant[size = {1}];\n")
    write(f, "DefineConstant[fsize = {1}];\n")

    # Write curve first
    for (index, p) in enumerate(c.points)
        write(f, "Point($(index)) = {$(p[1]), $(p[2]), 0, size};\n")
    end

    for (index, line) in enumerate(c.segments)
        write(f, "Line($(index)) = {$(line[1]), $(line[2])};\n")
    end
    # Mark crack as physical line 1
    crack = join(map(string, 1:length(c.segments)), ", ")
    write(f, "crack[] = {$(crack)};\n")
    write(f, "Physical Line(1) = {crack[]};\n")

    # Write bdry, points/lines start their enumeration after = consistent
    # with how bdry points are numbered in bdry lines
    n = length(c.points)
    for (i, p) in enumerate(bdry_pts)
        index = i + n
        write(f, "Point($(index)) = {$(p[1]), $(p[2]), 0, fsize};\n")
    end

    n = length(c.segments)
    for (i, line) in enumerate(bdry_lines)
        index = i + n
        write(f, "Line($(index)) = {$(line[1]), $(line[2])};\n")
    end

    loop = join(map(string, (n+1):(n+length(bdry_lines))), ", ")
    write(f, "Curve Loop(1) = {$(loop)};\n")
    write(f, "Plane Surface(1) = {1};\n")
    write(f, "Physical Surface(1) = {1};\n")
    # The magic
    write(f, "Line{crack[]} In Surface{1};\n")

    close(f)
end


text:
# examples UnitCells
chainunitcell = UnitCell{1}(1, [0.], +1.)
squareunitcell = UnitCell{2}(1, [(0, 0)], [(0, 1), (1, 0)])
triangularunitcell = UnitCell{2}(1, [(0, 0)], [(1, 0), (0.5, sin(pi/3))])
kagomestripunitcell = UnitCell{1}(3, [0., 0., 0.], +1.)
honeycombunitcell = UnitCell{2}(2, [(0, 0), (0, 1)], [(sqrt(3), 0), (sqrt(3)/2, 3/2)])

# example types
spinhalf = SpinType(2)

sz, sp, sm = spinoperators(1/2)

# example models
function j1j2model(lx::Int, j1::Float64, j2::Float64;
                   boundary::Symbol=:OBC,
                   symmetry::Type{<:AbstractCharge}=Trivial)

    heis = nbodyopexpansion(2,
                            0.25 * (ringexchangeoperator(2) - I(4)),
                            mode=:SYMBOL)

    heis1 = QModelInteraction{1, 2, Float64}(
        j1,
        (1, 1),
        ((0,), (1,)), heis)
    heis2 = QModelInteraction{1, 2, Float64}(
        j2,
        (1, 1),
        ((0,), (2,)), heis)
    terms = QModelInteraction[]
    j1 != 0 && push!(terms, heis1)
    j2 != 0 && push!(terms, heis2)

    UnitCellQModel(spinhalf,
                   QLattice(chainunitcell, lx, boundary),
                   symmetry,
                   terms)
end

# triangular
# ls is (ly, lx)
function triangularspinmodel(ls::Tuple{Int, Int},
                             j1::Float64,
                             j2::Float64,
                             j3::Float64,
                             k1::Float64,
                             k2::Float64,
                             k3::Float64;
                             boundary::Tuple{Symbol,Symbol}=(:PBC, :OBC),
                             symmetry::Type{<:AbstractCharge}=Trivial)

    heis = nbodyopexpansion(2,
                            0.25 * (ringexchangeoperator(2) - I(4)),
                            mode=:SYMBOL)
    ly, lx = ls

    heis1 = QModelInteraction{2, 2, Float64}(
        j1,
        (1, 1),
        ((0, 0), (0, 1)), heis)

    heis2 = QModelInteraction{2, 2, Float64}(
        j2,
        (1, 1),
        ((0, 0), (1, 0)), heis)

    heis3 = QModelInteraction{2, 2, Float64}(
        j3,
        (1, 1),
        ((0, 0), (1, -1)), heis)

    R4 = nbodyopexpansion(4,
                          ringexchangeoperator(4) - 0.25 * I(16),
                          mode=:SYMBOL)

    ring1 = QModelInteraction{2, 4, Float64}(
        k1,
        (1, 1, 1, 1),
        ((0, 0), (0, 1), (1, 1), (1, 0)),
        R4)

    ring2 = QModelInteraction{2, 4, Float64}(
        k2,
        (1, 1, 1, 1),
        ((0, 0), (0, 1), (1, 0), (1, -1)),
        R4)

    ring3 = QModelInteraction{2, 4, Float64}(
        k3,
        (1, 1, 1, 1),
        ((0, 0), (1, 0), (0, 1), (-1, 1)),
        R4)

    terms = QModelInteraction[]
    j1 != 0 && push!(terms, heis1)
    j2 != 0 && push!(terms, heis2)
    j3 != 0 && push!(terms, heis3)
    k1 != 0 && push!(terms, ring1)
    k2 != 0 && push!(terms, ring2)
    k3 != 0 && push!(terms, ring3)

    UnitCellQModel(spinhalf,
                   QLattice(triangularunitcell, (ly, lx), boundary),
                   symmetry,
                   terms)
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
using Test

@testset "spectral" begin
    RankedArray = MatrixNetworks.RankedArray

    @testset "RankedArray" begin
        n = 10
        v = rand(n)
        p = sortperm(v,rev=true)
        r = RankedArray(p)

        @test haskey(r,3) == true
        @test haskey(r,-1) == false
        @test haskey(r,0) == false
        @test haskey(r,n+1) == false

        @test_throws DimensionMismatch sweepcut(sparse(zeros(1,1)),v)
    end
    @testset "integers" begin
        n = 6
        x = collect(1.:Float64(n))
        A = sparse(1:n-1,2:n,1,n,n)
        A = A + A' + sparse(1.0I,n,n)
        profile = sweepcut(A,x)
        @test argmin(profile.conductance) == 3
        @test all(profile.cut .== 1)
    end
    @testset "floats" begin
        n = 6
        x = collect(1.:Float64(n))
        A = sparse(1:n-1,2:n,0.5,n,n)
        A = A + A' + sparse(1.5I,n,n)
        profile = sweepcut(A,x)
        @test argmin(profile.conductance) == 3
        @test all(profile.cut .== 0.5)
    end
    @testset "complement set" begin
        # Test an example where we need to take
        # the complement set
        n = 6
        A = sparse(1:n-1,2:n,1,n,n)
        A = A + A'
        x = [6,3,5,4,2,1]
        profile = sweepcut(A,x)
        @test length(bestset(profile)) == 2
    end
    @testset "spectral_cut" begin
        n = 50
        A = sparse(1:n-1,2:n,1,n,n)
        A = A + A'
        output = spectral_cut(A,true,true)

        (x,lam2) = fiedler_vector(sparse(zeros(1,1)))
        profile = sweepcut(sparse(zeros(1,1)),x)
        @test isempty(profile.conductance)

        @test_throws ArgumentError fiedler_vector(sparse(1.0I,2,2) + sparse([1],[2],1.,2,2))

        @test isempty(spectral_cut(sparse(zeros(5,5)),true,true).set)
        n = 50
        A = sparse(1:n-1,2:n,1.,n,n)
        A = A + A'
        @test length(spectral_cut(A).set) == 25

        n = 100
        A = sparse(ones(n,n))
        @test length(spectral_cut(A).set) == 50

        M = MatrixNetwork(A)
        @test length(spectral_cut(M).set) == 50

        @test_throws ArgumentError spectral_cut(sparse(1.0I,2,2) + spdiagm(1=>[1]), true, false)
        @test_throws ArgumentError spectral_cut(-sparse(1.0I,2,2), false, false)
    end
    dtol = 1.e-8 # default tolerance
    @testset "fiedler_vector" begin
        v,lam2 = fiedler_vector(sparse(zeros(1,1)))
        @test all(v .== [0.])
        # Newman's netscience graph in case anyone care ...
        A = sparse([2,3,4,5,16,44,113,131,250,259,1,3,1,2,1,5,13,14,15,16,44,45,46,47,61,126,127,128,146,152,153,154,164,165,166,176,177,249,250,274,313,314,323,324,330,371,373,374,1,4,15,16,44,45,46,47,176,177,199,201,202,204,231,235,236,237,238,249,250,254,298,313,314,373,374,7,8,6,8,190,191,192,193,6,7,26,62,63,64,65,137,189,342,343,344,10,11,12,9,11,12,67,68,69,9,10,12,9,10,11,4,14,15,16,17,18,19,20,274,4,13,15,16,4,5,13,14,16,45,46,47,176,177,278,279,334,366,367,368,1,4,5,13,14,15,45,46,47,153,154,176,177,249,250,313,314,323,324,371,373,13,18,29,58,172,201,258,261,365,13,17,58,172,201,258,261,365,13,20,208,13,19,22,23,24,33,109,220,221,232,233,268,287,288,21,23,24,21,22,24,50,51,52,54,55,220,227,228,21,22,23,79,140,220,229,232,233,268,26,27,28,8,25,27,28,40,95,104,105,106,107,108,124,125,155,197,198,231,234,251,273,295,296,297,306,315,316,317,25,26,28,25,26,27,17,31,32,33,34,35,36,51,76,219,30,32,33,34,30,31,33,34,35,36,51,76,216,217,218,219,307,369,370,21,30,31,32,34,35,36,109,219,220,221,30,31,32,33,30,32,33,36,219,30,32,33,35,38,37,304,305,40,26,39,173,174,175,239,240,282,326,42,43,241,242,243,244,41,43,52,170,241,242,243,244,275,276,277,364,41,42,275,276,277,1,4,5,66,4,5,15,16,46,47,176,177,323,324,4,5,15,16,45,47,176,177,4,5,15,16,45,46,176,177,49,48,67,74,181,182,183,226,23,51,52,53,23,30,32,50,52,76,95,100,160,169,170,307,308,322,337,23,42,50,51,76,100,116,117,169,170,266,267,291,364,50,170,23,55,23,54,57,100,110,194,195,56,17,18,59,60,58,60,58,59,304,357,358,359,4,164,165,166,8,63,64,65,8,62,64,8,62,63,8,62,264,265,299,300,301,302,345,346,44,100,101,102,110,111,10,49,68,69,70,71,72,73,74,75,80,81,121,122,123,169,248,285,363,10,67,69,80,81,10,67,68,70,71,75,80,81,285,67,69,71,72,73,75,119,150,214,303,328,329,348,349,350,351,378,379,67,69,70,72,73,75,67,70,71,73,119,303,67,70,71,72,150,49,67,67,69,70,71,285,30,32,51,52,169,170,78,79,77,79,24,77,78,67,68,69,81,67,68,69,80,121,122,123,83,84,85,86,87,88,89,82,84,85,86,87,88,89,262,263,82,83,85,86,87,88,89,82,83,84,86,87,88,89,106,260,262,263,82,83,84,85,87,88,89,106,260,262,263,82,83,84,85,86,88,89,82,83,84,85,86,87,89,106,260,262,263,82,83,84,85,86,87,88,91,92,90,92,130,131,132,133,134,188,90,91,130,131,132,133,134,188,94,95,96,97,98,99,93,95,96,97,98,99,26,51,93,94,96,97,98,99,178,179,180,286,308,337,338,339,362,93,94,95,97,98,99,141,93,94,95,96,98,99,178,362,93,94,95,96,97,99,93,94,95,96,97,98,51,52,56,66,101,102,103,110,111,145,194,195,66,100,102,103,311,66,100,101,159,280,360,361,100,101,26,105,106,107,108,26,104,106,107,26,85,86,88,104,105,107,185,260,26,104,105,106,108,283,284,372,26,104,107,125,155,156,158,21,33,56,66,100,111,194,195,66,100,110,113,114,115,1,112,114,115,131,135,136,189,259,312,352,353,354,355,356,112,113,115,131,135,312,112,113,114,172,347,52,117,318,319,320,321,52,116,119,120,209,70,72,118,120,214,303,348,349,350,351,118,119,67,81,122,123,67,81,121,123,169,248,67,81,121,122,26,125,26,108,124,234,255,256,4,127,128,129,327,330,4,126,128,129,327,330,4,126,127,129,327,330,340,341,126,127,128,91,92,131,132,133,134,1,91,92,113,114,130,132,133,134,135,136,259,91,92,130,131,133,134,188,91,92,130,131,132,134,91,92,130,131,132,133,113,114,131,189,352,353,354,355,356,113,131,8,138,342,343,137,140,24,139,206,207,229,230,96,143,144,145,142,144,145,142,143,145,100,142,143,144,194,4,147,148,146,148,146,147,150,151,270,293,70,73,149,151,270,271,293,294,149,150,4,4,16,154,4,16,153,26,108,156,157,158,108,155,155,108,155,102,51,161,162,163,322,160,162,160,161,163,160,162,4,61,165,166,4,61,164,166,4,61,164,165,168,169,170,167,169,170,205,51,52,67,76,122,167,168,170,205,248,289,290,291,292,42,51,52,53,76,167,168,169,266,267,336,364,172,17,18,115,171,325,347,40,174,175,40,173,175,40,173,174,282,4,5,15,16,45,46,47,177,4,5,15,16,45,46,47,176,95,97,179,180,95,178,95,178,49,182,183,49,181,183,49,181,182,226,185,186,187,309,310,106,184,186,187,309,310,184,185,187,309,310,184,185,186,91,92,132,8,113,135,7,191,192,193,7,190,7,190,193,7,190,192,56,100,110,145,195,196,56,100,110,194,194,26,26,251,5,200,201,202,203,204,258,199,201,202,5,17,18,199,200,202,203,204,245,252,253,254,258,298,5,199,200,201,203,204,258,199,201,202,5,199,201,202,245,298,168,169,140,207,140,206,375,376,377,19,118,211,212,213,214,215,210,212,213,214,215,210,211,213,214,215,222,223,224,225,210,211,212,214,215,70,119,210,211,212,213,215,303,348,349,350,351,210,211,212,213,214,32,217,218,32,216,218,32,216,217,30,32,33,35,21,23,24,33,221,21,33,220,212,223,224,225,212,222,224,225,212,222,223,225,212,222,223,224,49,183,23,228,23,227,24,140,140,5,26,232,233,234,235,236,237,238,239,240,246,257,297,21,24,231,233,268,21,24,231,232,268,26,125,231,5,231,5,231,237,238,239,240,245,246,247,257,5,231,236,238,5,231,236,237,40,231,236,240,246,257,326,40,231,236,239,246,257,326,41,42,242,243,244,41,42,241,243,244,41,42,241,242,244,41,42,241,242,243,201,204,236,246,247,298,231,236,239,240,245,247,257,236,245,246,67,122,169,4,5,16,1,4,5,16,313,314,26,198,201,253,201,252,5,201,125,256,125,255,231,236,239,240,246,17,18,199,201,202,1,113,131,85,86,88,106,17,18,83,85,86,88,263,83,85,86,88,262,65,265,65,264,299,300,301,302,52,170,267,52,170,266,21,24,232,233,270,271,272,149,150,269,271,272,293,294,150,269,270,294,269,270,26,4,13,42,43,276,277,42,43,275,277,42,43,275,276,15,279,15,278,102,281,360,361,280,40,175,107,284,107,283,67,69,75,95,21,288,21,287,169,290,169,289,291,292,52,169,290,169,290,149,150,270,150,270,271,26,296,26,295,26,231,5,201,204,245,65,265,300,301,302,65,265,299,301,65,265,299,300,65,265,299,70,72,119,214,328,329,348,349,350,351,378,379,38,60,305,357,358,359,38,304,26,32,51,51,95,337,184,185,186,310,184,185,186,309,101,113,114,4,5,16,250,314,4,5,16,250,313,26,316,317,26,315,317,26,315,316,116,319,320,321,116,318,320,321,116,318,319,321,116,318,319,320,51,160,4,16,45,324,4,16,45,323,172,40,239,240,126,127,128,70,303,329,70,303,328,4,126,127,128,332,333,334,335,331,333,334,335,331,332,334,335,15,331,332,333,335,366,367,368,331,332,333,334,170,51,95,308,95,339,95,338,128,341,128,340,8,137,343,8,137,342,8,65,346,65,345,115,172,70,119,214,303,349,350,351,70,119,214,303,348,350,351,70,119,214,303,348,349,351,70,119,214,303,348,349,350,113,135,353,354,355,356,113,135,352,354,355,356,113,135,352,353,355,356,113,135,352,353,354,356,113,135,352,353,354,355,60,304,358,359,60,304,357,359,60,304,357,358,102,280,361,102,280,360,95,97,67,42,52,170,17,18,15,334,367,368,15,334,366,368,15,334,366,367,32,370,32,369,4,16,107,4,5,16,374,4,5,373,207,376,377,207,375,377,207,375,376,70,303,379,70,303,378],
                   [1,1,1,1,1,1,1,1,1,1,2,2,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8,9,9,9,10,10,10,10,10,10,11,11,11,12,12,12,13,13,13,13,13,13,13,13,13,14,14,14,14,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,17,17,17,17,17,17,17,17,17,18,18,18,18,18,18,18,18,19,19,19,20,20,21,21,21,21,21,21,21,21,21,21,21,21,22,22,22,23,23,23,23,23,23,23,23,23,23,23,24,24,24,24,24,24,24,24,24,24,25,25,25,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,27,27,27,28,28,28,29,30,30,30,30,30,30,30,30,30,31,31,31,31,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,33,33,33,33,33,33,33,33,33,33,33,34,34,34,34,35,35,35,35,35,36,36,36,36,37,38,38,38,39,40,40,40,40,40,40,40,40,40,41,41,41,41,41,41,42,42,42,42,42,42,42,42,42,42,42,42,43,43,43,43,43,44,44,44,44,45,45,45,45,45,45,45,45,45,45,46,46,46,46,46,46,46,46,47,47,47,47,47,47,47,47,48,49,49,49,49,49,49,49,50,50,50,50,51,51,51,51,51,51,51,51,51,51,51,51,51,51,51,52,52,52,52,52,52,52,52,52,52,52,52,52,52,53,53,54,54,55,55,56,56,56,56,56,57,58,58,58,58,59,59,60,60,60,60,60,60,61,61,61,61,62,62,62,62,63,63,63,64,64,64,65,65,65,65,65,65,65,65,65,65,66,66,66,66,66,66,67,67,67,67,67,67,67,67,67,67,67,67,67,67,67,67,67,67,67,68,68,68,68,68,69,69,69,69,69,69,69,69,69,70,70,70,70,70,70,70,70,70,70,70,70,70,70,70,70,70,70,71,71,71,71,71,71,72,72,72,72,72,72,73,73,73,73,73,74,74,75,75,75,75,75,76,76,76,76,76,76,77,77,78,78,79,79,79,80,80,80,80,81,81,81,81,81,81,81,82,82,82,82,82,82,82,83,83,83,83,83,83,83,83,83,84,84,84,84,84,84,84,85,85,85,85,85,85,85,85,85,85,85,86,86,86,86,86,86,86,86,86,86,86,87,87,87,87,87,87,87,88,88,88,88,88,88,88,88,88,88,88,89,89,89,89,89,89,89,90,90,91,91,91,91,91,91,91,91,92,92,92,92,92,92,92,92,93,93,93,93,93,93,94,94,94,94,94,94,95,95,95,95,95,95,95,95,95,95,95,95,95,95,95,95,95,96,96,96,96,96,96,96,97,97,97,97,97,97,97,97,98,98,98,98,98,98,99,99,99,99,99,99,100,100,100,100,100,100,100,100,100,100,100,100,101,101,101,101,101,102,102,102,102,102,102,102,103,103,104,104,104,104,104,105,105,105,105,106,106,106,106,106,106,106,106,106,107,107,107,107,107,107,107,107,108,108,108,108,108,108,108,109,109,110,110,110,110,110,110,111,111,111,112,112,112,113,113,113,113,113,113,113,113,113,113,113,113,113,113,113,114,114,114,114,114,114,115,115,115,115,115,116,116,116,116,116,116,117,117,118,118,118,119,119,119,119,119,119,119,119,119,119,120,120,121,121,121,121,122,122,122,122,122,122,123,123,123,123,124,124,125,125,125,125,125,125,126,126,126,126,126,126,127,127,127,127,127,127,128,128,128,128,128,128,128,128,129,129,129,130,130,130,130,130,130,131,131,131,131,131,131,131,131,131,131,131,131,132,132,132,132,132,132,132,133,133,133,133,133,133,134,134,134,134,134,134,135,135,135,135,135,135,135,135,135,136,136,137,137,137,137,138,139,140,140,140,140,140,140,141,142,142,142,143,143,143,144,144,144,145,145,145,145,145,146,146,146,147,147,148,148,149,149,149,149,150,150,150,150,150,150,150,150,151,151,152,153,153,153,154,154,154,155,155,155,155,155,156,156,157,158,158,159,160,160,160,160,160,161,161,162,162,162,163,163,164,164,164,164,165,165,165,165,166,166,166,166,167,167,167,168,168,168,168,169,169,169,169,169,169,169,169,169,169,169,169,169,169,170,170,170,170,170,170,170,170,170,170,170,170,171,172,172,172,172,172,172,173,173,173,174,174,174,175,175,175,175,176,176,176,176,176,176,176,176,177,177,177,177,177,177,177,177,178,178,178,178,179,179,180,180,181,181,181,182,182,182,183,183,183,183,184,184,184,184,184,185,185,185,185,185,185,186,186,186,186,186,187,187,187,188,188,188,189,189,189,190,190,190,190,191,191,192,192,192,193,193,193,194,194,194,194,194,194,195,195,195,195,196,197,198,198,199,199,199,199,199,199,199,200,200,200,201,201,201,201,201,201,201,201,201,201,201,201,201,201,202,202,202,202,202,202,202,203,203,203,204,204,204,204,204,204,205,205,206,206,207,207,207,207,207,208,209,210,210,210,210,210,211,211,211,211,211,212,212,212,212,212,212,212,212,212,213,213,213,213,213,214,214,214,214,214,214,214,214,214,214,214,214,215,215,215,215,215,216,216,216,217,217,217,218,218,218,219,219,219,219,220,220,220,220,220,221,221,221,222,222,222,222,223,223,223,223,224,224,224,224,225,225,225,225,226,226,227,227,228,228,229,229,230,231,231,231,231,231,231,231,231,231,231,231,231,231,231,232,232,232,232,232,233,233,233,233,233,234,234,234,235,235,236,236,236,236,236,236,236,236,236,236,237,237,237,237,238,238,238,238,239,239,239,239,239,239,239,240,240,240,240,240,240,240,241,241,241,241,241,242,242,242,242,242,243,243,243,243,243,244,244,244,244,244,245,245,245,245,245,245,246,246,246,246,246,246,246,247,247,247,248,248,248,249,249,249,250,250,250,250,250,250,251,251,252,252,253,253,254,254,255,255,256,256,257,257,257,257,257,258,258,258,258,258,259,259,259,260,260,260,260,261,261,262,262,262,262,262,263,263,263,263,263,264,264,265,265,265,265,265,265,266,266,266,267,267,267,268,268,268,268,269,269,269,270,270,270,270,270,270,270,271,271,271,271,272,272,273,274,274,275,275,275,275,276,276,276,276,277,277,277,277,278,278,279,279,280,280,280,280,281,282,282,283,283,284,284,285,285,285,286,287,287,288,288,289,289,290,290,290,290,291,291,291,292,292,293,293,293,294,294,294,295,295,296,296,297,297,298,298,298,298,299,299,299,299,299,300,300,300,300,301,301,301,301,302,302,302,303,303,303,303,303,303,303,303,303,303,303,303,304,304,304,304,304,304,305,305,306,307,307,308,308,308,309,309,309,309,310,310,310,310,311,312,312,313,313,313,313,313,314,314,314,314,314,315,315,315,316,316,316,317,317,317,318,318,318,318,319,319,319,319,320,320,320,320,321,321,321,321,322,322,323,323,323,323,324,324,324,324,325,326,326,326,327,327,327,328,328,328,329,329,329,330,330,330,330,331,331,331,331,332,332,332,332,333,333,333,333,334,334,334,334,334,334,334,334,335,335,335,335,336,337,337,337,338,338,339,339,340,340,341,341,342,342,342,343,343,343,344,345,345,346,346,347,347,348,348,348,348,348,348,348,349,349,349,349,349,349,349,350,350,350,350,350,350,350,351,351,351,351,351,351,351,352,352,352,352,352,352,353,353,353,353,353,353,354,354,354,354,354,354,355,355,355,355,355,355,356,356,356,356,356,356,357,357,357,357,358,358,358,358,359,359,359,359,360,360,360,361,361,361,362,362,363,364,364,364,365,365,366,366,366,366,367,367,367,367,368,368,368,368,369,369,370,370,371,371,372,373,373,373,373,374,374,374,375,375,375,376,376,376,377,377,377,378,378,378,379,379,379],
                   1., 379, 379)
        v,lam2 = fiedler_vector(A)
        @test abs(lam2 - 3.026826881103339417e-03) <= dtol # result from LAPACK graph_eigs
        @test abs(v[1] - 1.878339361656544346e-02) <= dtol # result from Matlab
        # check inverse participation ratio score from graph_eigs
        @test abs.(sum(abs.(v).^4)/(sum(abs.(v).^2)^2) - 6.761596909820281540e-03) <= dtol

        n = 25
        A = sparse(1:n-1,2:n,1.,n,n)
        A = A + A'
        elam2 = 2*sin(pi/2*1/(n-1.))^2 # exact lam2
        v,lam2 = fiedler_vector(A)
        @show elam2
        @show lam2
        @test abs(lam2 - elam2) <= dtol

        n = 500
        A = sparse(1:n-1,2:n,1.,n,n)
        A = A + A'
        elam2 = 2*sin(pi/2*1/(n-1.))^2 # exact lam2
        v,lam2 = fiedler_vector(A;nev=6)
        @show elam2
        @show lam2
        @test abs(lam2 - elam2) <= dtol

        @test length(bestset(sweepcut(A,v))) == 250
    end
    @testset "Fiedler partition" begin
        # Make sure that the Fiedler partition will run properly on the largest
        # connected component.  The following network is the in-recip weighted motif
        # adjacency matrix for the Florida Bay food web.
        W = sparse([24,25,26,30,32,35,36,38,39,40,41,45,48,124,41,44,45,49,124,41,44,45,49,124,41,44,45,49,124,91,124,30,41,44,45,49,91,124,16,17,18,19,20,21,23,27,123,15,123,15,123,15,123,15,123,15,123,15,123,123,126,15,123,8,124,8,124,8,96,100,124,15,123,30,124,8,14,29,31,42,124,30,34,36,44,49,50,51,61,124,8,34,61,124,31,32,61,8,96,100,124,8,31,37,42,96,100,124,36,49,96,100,124,8,96,100,124,8,96,100,124,8,96,100,124,8,9,10,11,14,96,100,124,30,36,91,96,100,124,60,62,66,96,99,100,9,10,11,14,31,60,62,96,100,124,8,9,10,11,14,60,62,96,100,124,8,60,62,124,9,10,11,14,31,37,60,62,96,100,124,31,60,62,124,31,60,62,124,66,99,66,99,60,62,66,99,43,44,45,48,49,50,51,59,62,68,94,95,100,31,32,34,43,44,45,48,49,50,51,59,60,68,94,95,100,96,100,96,100,66,99,43,57,58,59,65,68,69,72,99,60,62,66,99,66,99,66,99,13,14,42,124,96,100,60,62,60,62,26,35,36,37,38,39,40,41,42,43,44,45,49,63,64,93,100,43,57,58,59,65,66,68,69,72,26,35,36,37,38,39,40,41,42,43,44,45,49,60,62,63,64,93,96,15,16,17,18,19,20,21,22,23,27,126,8,9,10,11,13,14,24,25,26,29,30,31,32,35,36,37,38,39,40,41,42,44,45,48,49,50,51,91,22,123],
        [8,8,8,8,8,8,8,8,8,8,8,8,8,8,9,9,9,9,9,10,10,10,10,10,11,11,11,11,11,13,13,14,14,14,14,14,14,14,15,15,15,15,15,15,15,15,15,16,16,17,17,18,18,19,19,20,20,21,21,22,22,23,23,24,24,25,25,26,26,26,26,27,27,29,29,30,30,30,30,30,30,31,31,31,31,31,31,31,31,31,32,32,32,32,34,34,34,35,35,35,35,36,36,36,36,36,36,36,37,37,37,37,37,38,38,38,38,39,39,39,39,40,40,40,40,41,41,41,41,41,41,41,41,42,42,42,42,42,42,43,43,43,43,43,43,44,44,44,44,44,44,44,44,44,44,45,45,45,45,45,45,45,45,45,45,48,48,48,48,49,49,49,49,49,49,49,49,49,49,49,50,50,50,50,51,51,51,51,57,57,58,58,59,59,59,59,60,60,60,60,60,60,60,60,60,60,60,60,60,61,61,61,62,62,62,62,62,62,62,62,62,62,62,62,62,63,63,64,64,65,65,66,66,66,66,66,66,66,66,66,68,68,68,68,69,69,72,72,91,91,91,91,93,93,94,94,95,95,96,96,96,96,96,96,96,96,96,96,96,96,96,96,96,96,96,99,99,99,99,99,99,99,99,99,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,123,123,123,123,123,123,123,123,123,123,123,124,124,124,124,124,124,124,124,124,124,124,124,124,124,124,124,124,124,124,124,124,124,124,124,124,124,124,124,126,126],
                   [1,1,1,1,1,1,1,1,1,1,1,1,1,13,1,1,1,1,4,1,1,1,1,4,1,1,1,1,4,1,1,1,1,1,1,1,1,6,1,1,1,1,1,1,1,1,8,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,5,1,1,1,1,1,1,1,1,6,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,4,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,5,1,1,1,1,1,3,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,5,1,1,1,1,1,1,1,1,1,5,1,1,1,1,1,1,1,1,1,1,1,1,1,1,6,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,12,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,12,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,8,1,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,16,1,1,1,1,1,8,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,16,8,1,1,1,1,1,1,1,1,1,1,13,4,4,4,1,6,1,1,1,1,5,6,1,1,4,2,1,1,1,5,3,5,5,1,6,1,1,3,1,1])

        sc = spectral_cut(W,true,true)
        @test length(sc.set) == 9
        @test length(intersect([57,58,59,65,66,68,69,72,99], sc.set)) == 9
    end
end


text:
mutable struct WallTimer
    starttime_ns::typeof(Base.time_ns())
    paused_elapsed_ns::typeof(Base.time_ns())
    WallTimer() = new(0,0)
end

function start!(timer::WallTimer)
    timer.starttime_ns = (Base.time_ns)()
    return nothing
end
started(timer::WallTimer) = (timer.starttime_ns ≠ 0)

""" Return seconds since timer was started or 0 if not yet started. """
function elapsed(timer::WallTimer)
    local elapsedtime_ns = (Base.time_ns)() - timer.starttime_ns
    return started(timer) * float(elapsedtime_ns) / 1000000000
end

function pause!(timer::WallTimer)
    timer.paused_elapsed_ns = (Base.time_ns)() - timer.starttime_ns
    return nothing
end
function unpause!(timer::WallTimer)
    timer.starttime_ns = (Base.time_ns)()
    timer.starttime_ns -= timer.paused_elapsed_ns;
    return nothing
end

t = WallTimer()
start!(t)
elapsed(t)
pause!(t)
unpause!(t)
elapsed(t)

# -----------------

""" Like WallTimer, but uses in-game time instead of wall time. """
mutable struct GameTimer
    elapsed_ns::typeof(Base.time_ns())
    started::Bool
    paused::Bool
    GameTimer() = new(0, false, false)
end

function update!(timer::GameTimer, dt)
    if timer.started && !timer.paused
        timer.elapsed_ns += round(dt * 1e9)  # secs to nanoseconds
    end
    return nothing
end

function start!(timer::GameTimer)
    timer.started = true
    return nothing
end
started(timer::GameTimer) = timer.started
elapsed(timer::GameTimer) = timer.elapsed_ns / 1e9

function pause!(timer::GameTimer)
    timer.paused = true
    return nothing
end
function unpause!(timer::GameTimer)
    timer.paused = false
    return nothing
end

t = GameTimer()
start!(t)
update!(t, 2.0)
elapsed(t)
pause!(t)
update!(t, 2.0)
unpause!(t)
elapsed(t)


text:
#####
##### Navier-Stokes and tracer advection equations
#####

""" Store previous value of the source term and calculate current source term. """
function calculate_interior_tendency_contributions!(G, arch, grid, coriolis, buoyancy, surface_waves, closure,
                                                    U, C, pHY′, K, F, clock)

    # Manually choose thread-block layout here as it's ~20% faster.
    # See: https://github.com/climate-machine/Oceananigans.jl/pull/308
    Nx, Ny, Nz = grid.Nx, grid.Ny, grid.Nz

    if Nx == 1
        Tx, Ty = 1, min(256, Ny)
        Bx, By, Bz = Tx, floor(Int, Ny/Ty), Nz
    elseif Ny == 1
        Tx, Ty = min(256, Nx), 1
        Bx, By, Bz = floor(Int, Nx/Tx), Ty, Nz
    else
        Tx, Ty = 16, 16
        Bx, By, Bz = floor(Int, Nx/Tx), floor(Int, Ny/Ty), Nz
    end

    @launch(device(arch), threads=(Tx, Ty), blocks=(Bx, By, Bz),
            calculate_Gu!(G.u, grid, coriolis, surface_waves, closure, U, C, K, F, pHY′, clock))

    @launch(device(arch), threads=(Tx, Ty), blocks=(Bx, By, Bz),
            calculate_Gv!(G.v, grid, coriolis, surface_waves, closure, U, C, K, F, pHY′, clock))

    @launch(device(arch), threads=(Tx, Ty), blocks=(Bx, By, Bz),
            calculate_Gw!(G.w, grid, coriolis, surface_waves, closure, U, C, K, F, clock))

    for tracer_index in 1:length(C)
        @inbounds Gc = G[tracer_index+3]
        @inbounds Fc = F[tracer_index+3]
        @inbounds  c = C[tracer_index]

        @launch(device(arch), threads=(Tx, Ty), blocks=(Bx, By, Bz),
                calculate_Gc!(Gc, grid, c, Val(tracer_index), closure, buoyancy, U, C, K, Fc, clock))
    end

    return nothing
end

"""
    calculate_velocity_tendencies_on_boundaries!(tendency_calculation_args...)

Calculate the velocity tendencies *on* east, north, and top boundaries, when the
x-, y-, or z- directions have a `Bounded` topology.
"""
function calculate_velocity_tendencies_on_boundaries!(tendency_calculation_args...)
    calculate_east_boundary_Gu!(tendency_calculation_args...)
    calculate_north_boundary_Gv!(tendency_calculation_args...)
    calculate_top_boundary_Gw!(tendency_calculation_args...)
    return nothing
end

# Fallbacks for non-bounded topologies in x-, y-, and z-directions
@inline calculate_east_boundary_Gu!(args...) = nothing
@inline calculate_north_boundary_Gv!(args...) = nothing
@inline calculate_top_boundary_Gw!(args...) = nothing

"""
    calculate_east_boundary_Gu!(G, arch, grid::AbstractGrid{FT, <:Bounded},
                                coriolis, buoyancy, surface_waves, closure,
                                U, C, pHY′, K, F, clock) where FT

Calculate `Gu` on east boundaries when the x-direction has `Bounded` topology.
"""
function calculate_east_boundary_Gu!(G, arch, grid::AbstractGrid{FT, <:Bounded},
                                     coriolis, buoyancy, surface_waves, closure,
                                     U, C, pHY′, K, F, clock) where FT

    @launch(device(arch), config=launch_config(grid, :yz),
            _calculate_east_boundary_Gu!(G.u, grid, coriolis, surface_waves,
                                         closure, U, C, K, F, pHY′, clock))

    return nothing
end

"""
    calculate_north_boundary_Gu!(G, arch, grid::AbstractGrid{FT, <:Bounded},
                                 coriolis, buoyancy, surface_waves, closure,
                                 U, C, pHY′, K, F, clock) where FT

Calculate `Gv` on north boundaries when the y-direction has `Bounded` topology.
"""
function calculate_north_boundary_Gv!(G, arch, grid::AbstractGrid{FT, TX, <:Bounded},
                                      coriolis, buoyancy, surface_waves, closure,
                                      U, C, pHY′, K, F, clock) where {FT, TX}

    @launch(device(arch), config=launch_config(grid, :xz),
            _calculate_north_boundary_Gv!(G.v, grid, coriolis, surface_waves,
                                          closure, U, C, K, F, pHY′, clock))

    return nothing
end

"""
    calculate_top_boundary_Gw!(G, arch, grid::AbstractGrid{FT, <:Bounded},
                               coriolis, buoyancy, surface_waves, closure,
                               U, C, pHY′, K, F, clock) where FT

Calculate `Gw` on top boundaries when the z-direction has `Bounded` topology.
"""
function calculate_top_boundary_Gw!(G, arch, grid::AbstractGrid{FT, TX, TY, <:Bounded},
                                    coriolis, buoyancy, surface_waves, closure,
                                    U, C, pHY′, K, F, clock) where {FT, TX, TY}

    @launch(device(arch), config=launch_config(grid, :xy),
            _calculate_top_boundary_Gw!(G.w, grid, coriolis, surface_waves,
                                        closure, U, C, K, F, clock))

    return nothing
end

#####
##### Tendency calculators for u-velocity, aka x-velocity
#####

""" Calculate the right-hand-side of the u-velocity equation. """
function calculate_Gu!(Gu, grid, coriolis, surface_waves, closure, U, C, K, F, pHY′, clock)
    @loop_xyz i j k grid begin
        @inbounds Gu[i, j, k] = u_velocity_tendency(i, j, k, grid, coriolis, surface_waves,
                                                    closure, U, C, K, F, pHY′, clock)
    end
    return nothing
end

""" Calculate the right-hand-side of the u-velocity equation on the east boundary. """
function _calculate_east_boundary_Gu!(Gu, grid, coriolis, surface_waves,
                                      closure, U, C, K, F, pHY′, clock)
    i = grid.Nx + 1
    @loop_yz j k grid begin
        @inbounds Gu[i, j, k] = u_velocity_tendency(i, j, k, grid, coriolis, surface_waves,
                                                    closure, U, C, K, F, pHY′, clock)
    end
    return nothing
end

#####
##### Tendency calculators for v-velocity
#####

""" Calculate the right-hand-side of the v-velocity equation. """
function calculate_Gv!(Gv, grid, coriolis, surface_waves, closure, U, C, K, F, pHY′, clock)
    @loop_xyz i j k grid begin
        @inbounds Gv[i, j, k] = v_velocity_tendency(i, j, k, grid, coriolis, surface_waves,
                                                    closure, U, C, K, F, pHY′, clock)
    end
    return nothing
end

""" Calculate the right-hand-side of the v-velocity equation on the north boundary. """
function _calculate_north_boundary_Gv!(Gv, grid, coriolis, surface_waves,
                                       closure, U, C, K, F, pHY′, clock)
    j = grid.Ny + 1
    @loop_xz i k grid begin
        @inbounds Gv[i, j, k] = v_velocity_tendency(i, j, k, grid, coriolis, surface_waves,
                                                    closure, U, C, K, F, pHY′, clock)
    end
    return nothing
end

#####
##### Tendency calculators for w-velocity
#####

""" Calculate the right-hand-side of the w-velocity equation. """
function calculate_Gw!(Gw, grid, coriolis, surface_waves, closure, U, C, K, F, clock)
    @loop_xyz i j k grid begin
        @inbounds Gw[i, j, k] = w_velocity_tendency(i, j, k, grid, coriolis, surface_waves,
                                                    closure, U, C, K, F, clock)
    end
    return nothing
end

""" Calculate the right-hand-side of the w-velocity equation. """
function _calculate_top_boundary_Gw!(Gw, grid, coriolis, surface_waves, closure, U, C, K, F, clock)
    k = grid.Nz + 1
    @loop_xy i j grid begin
        @inbounds Gw[i, j, k] = w_velocity_tendency(i, j, k, grid, coriolis, surface_waves,
                                                    closure, U, C, K, F, clock)
    end
    return nothing
end

#####
##### Tracer(s)
#####

""" Calculate the right-hand-side of the tracer advection-diffusion equation. """
function calculate_Gc!(Gc, grid, c, tracer_index, closure, buoyancy, U, C, K, Fc, clock)
    @loop_xyz i j k grid begin
        @inbounds Gc[i, j, k] = tracer_tendency(i, j, k, grid, c, tracer_index,
                                                closure, buoyancy, U, C, K, Fc, clock)
    end
    return nothing
end

#####
##### Boundary contributions to tendencies due to user-prescribed fluxes
#####

""" Apply boundary conditions by adding flux divergences to the right-hand-side. """
function calculate_boundary_tendency_contributions!(Gⁿ, arch, U, C, clock, state)

    # Velocity fields
    for i in 1:3
        apply_x_bcs!(Gⁿ[i], U[i], arch, clock, state)
        apply_y_bcs!(Gⁿ[i], U[i], arch, clock, state)
        apply_z_bcs!(Gⁿ[i], U[i], arch, clock, state)
    end

    # Tracer fields
    for i in 4:length(Gⁿ)
        apply_x_bcs!(Gⁿ[i], C[i-3], arch, clock, state)
        apply_y_bcs!(Gⁿ[i], C[i-3], arch, clock, state)
        apply_z_bcs!(Gⁿ[i], C[i-3], arch, clock, state)
    end

    return nothing
end

#####
##### Vertical integrals
#####

"""
Update the hydrostatic pressure perturbation pHY′. This is done by integrating
the `buoyancy_perturbation` downwards:

    `pHY′ = ∫ buoyancy_perturbation dz` from `z=0` down to `z=-Lz`
"""
function update_hydrostatic_pressure!(pHY′, grid, buoyancy, C)
    @loop_xy i j grid begin
        @inbounds pHY′[i, j, grid.Nz] = - ℑzᵃᵃᶠ(i, j, grid.Nz+1, grid, buoyancy_perturbation, buoyancy, C) * ΔzF(i, j, grid.Nz+1, grid)
        @unroll for k in grid.Nz-1 : -1 : 1
            @inbounds pHY′[i, j, k] =
                pHY′[i, j, k+1] - ℑzᵃᵃᶠ(i, j, k+1, grid, buoyancy_perturbation, buoyancy, C) * ΔzF(i, j, k+1, grid)
        end
    end
    return nothing
end

"""
Compute the vertical velocity w by integrating the continuity equation from the bottom upwards

    `w^{n+1} = -∫ [∂/∂x (u^{n+1}) + ∂/∂y (v^{n+1})] dz`
"""
function compute_w_from_continuity!(model)
    @launch(device(model.architecture), config=launch_config(model.grid, :xy),
            _compute_w_from_continuity!(datatuple(model.velocities), model.grid))
    return nothing
end

function _compute_w_from_continuity!(U, grid)
    @loop_xy i j grid begin
        # U.w[i, j, 1] = 0 is enforced via halo regions.
        @unroll for k in 2:grid.Nz
            @inbounds U.w[i, j, k] = U.w[i, j, k-1] - ΔzC(i, j, k, grid) * hdivᶜᶜᵃ(i, j, k-1, grid, U.u, U.v)
        end
    end
    return nothing
end


text:
# ---
# title: 847. Shortest Path Visiting All Nodes
# id: problem847
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Breadth-first Search
# link: <https://leetcode.com/problems/shortest-path-visiting-all-nodes/description/>
# hidden: true
# ---
#
# An undirected, connected graph of N nodes (labeled `0, 1, 2, ..., N-1`) is
# given as `graph`.
#
# `graph.length = N`, and `j != i` is in the list `graph[i]` exactly once, if
# and only if nodes `i` and `j` are connected.
#
# Return the length of the shortest path that visits every node. You may start
# and stop at any node, you may revisit nodes multiple times, and you may reuse
# edges.
#
#
#
# **Example 1:**
#
#
#
#     Input: [[1,2,3],[0],[0],[0]]
#     Output: 4
#     **Explanation** : One possible path is [1,0,2,0,3]
#
# **Example 2:**
#
#
#
#     Input: [[1],[0,2,4],[1,3,4],[2],[1,2]]
#     Output: 4
#     **Explanation** : One possible path is [0,1,4,2,3]
#
#
#
#
# **Note:**
#
#   1. `1 <= graph.length <= 12`
#   2. `0 <= graph[i].length < graph.length`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#the majority of the following code is inspired by
#the DataFrames.jl package (https://github.com/JuliaStats/DataFrames.jl)
#DataFrames.jl is licensed under the MIT License, a copy of which
#can be found in this package's repo in LICENSE.md.

#This file contains functions that are useful for crossjoining
#arbitrary arrays, and are thus useful for building a tensor product
#structure.

function rep(x::AbstractVector, times::Integer = 1, each::Integer = 1)
    res = similar(x, each * times * length(x))
    i = 1
    for jdx in 1:times
        for idx in 1:length(x)
            tmp = x[idx]
            for kdx in 1:each
                res[i] = tmp
                i += 1
            end
        end
    end
    return res
end

function rep(x::AbstractArray, times::Integer = 1, each::Integer = 1)
    res = similar(x, each * times * size(x,1), size(x,2))
    i = 1
    for jdx in 1:times
        for idx in 1:size(x,1)
            tmp = x[idx,:]
            for kdx in 1:each
                res[i,:] = tmp
                i += 1
            end
        end
    end
    return res
end

crossjoin(A::Array, B::Array) = hcat(rep(A, 1, size(B,1)), rep(B, size(A,1), 1))

crossjoin(A::Array...) = reduce(crossjoin, A)