text:
using Plots
#using CUDA
#using CUDA.CUSPARSE
using Printf

#CUDA.allowscalar(false)

function Q_DYNAMIC!(dψδ, ψδ, p, t)

    #@printf "\r\t%f" t/p.year_seconds

    reject_step = p.reject_step
    if reject_step[1]
        return
    end

    nn = p.nn
    δNp = p.δNp
    Lw = p.Lw
    Δτ = p.vars.Δτ
    u = p.vars.u
    ge = p.vars.ge
    vf = p.vars.vf
    M = p.ops.M̃
    K = p.ops.K
    H̃ = p.ops.H̃
    JI = p.ops.JI
    Crr = p.ops.Crr
    RS = p.RS
    metrics = p.metrics
    ops = p.ops
    η = metrics.η
    b = p.RS.b
    uf2 = p.vars.uf2
    t_begin = p.vars.t_prev[2]

    xf1 = metrics.facecoord[1][1]
    yf1 = metrics.facecoord[2][1]

    ψ  = @view ψδ[1:nn]
    δ =  @view ψδ[nn + 1 : 2nn]
    dψ = @view dψδ[1:nn]
    V = @view dψδ[nn + 1 : 2nn]

    mod_data!(δ, uf2, ge, K, vf, RS, metrics, Lw, t-t_begin)

    u[:] = M \ ge

    HI = ops.HI[1]
    G = ops.G[1]
    Γ = ops.Γ[1]
    L = ops.L[1]
    sJ = metrics.sJ[1]

    Δτ[:] = - (HI * G * u + Crr * Γ * (δ ./ 2 - L * u)) ./ sJ

    for n in 1:nn

        if n <= δNp

            ψn = ψ[n]
            bn = b[n]
            τn = Δτ[n]
            ηn = η[n]

            #Vn = (2 * RS.V0 * sinh(τn / (RS.σn * RS.a))) / (exp(ψn/RS.a))

            # rootfinding with radiation damping ηV
            VR = abs(τn / ηn)
            VL = -VR
            Vn = V[n]
            obj_rs(V) = rateandstateQ(V, ψn, RS.σn, τn, ηn, RS.a, RS.V0)
            (Vn, _, iter) = newtbndv(obj_rs, VL, VR, Vn; ftol = 1e-14,
            atolx = 1e-14, rtolx = 1e-14)

            if !isfinite(Vn)
                reject_step[1] = true
                return
            end

            V[n] = Vn

            if bn != 0
                dψ[n] = (bn * RS.V0 / RS.Dc) * (exp((RS.f0 - ψn) / bn) - abs(Vn) / RS.V0)
            else
                dψ[n] = 0
            end

            if !isfinite(dψ[n])
                dψ .= 0
                reject_step[1] = true
                return
            end
        else
            V[n] = RS.Vp
        end

    end

    nothing

end

function Q_STATIC_MMS!(p)

    nn = p.nn
    Δτ = p.Δτ
    u = p.u
    ge = p.ge
    vf = p.vf
    M = p.ops.M̃
    K = p.ops.K
    JH = p.ops.JH
    RS = p.RS
    MMS = p.MMS
    B_p = p.B_p
    metrics = p.metrics
    ops = p.ops
    η = metrics.η
    b = p.b
    δ = 0
    ys = p.year_seconds
    t = p.t_final

    static_data_mms!(δ, ge, K, JH, vf, MMS, B_p, RS, metrics, t)

    u[:] = M \ ge

end

function Q_DYNAMIC_MMS_NOROOT!(dψδ, ψδ, p, t)

    #@printf "\r\t%f" t/31556952

    reject_step = p.reject_step
    if reject_step[1]
        return
    end

    nn = p.nn
    Δτ = p.Δτ
    u = p.u
    ge = p.ge
    vf = p.vf
    M = p.ops.M̃
    K = p.ops.K
    H̃ = p.ops.H̃
    JH = p.ops.JH
    Crr = p.ops.Crr
    RS = p.RS
    MMS = p.MMS
    B_p = p.B_p
    metrics = p.metrics
    ops = p.ops
    η = metrics.η
    b = p.b
    f1_source = p.f1_source
    h = p.ops.hmin

    xf1 = metrics.facecoord[1][1]
    yf1 = metrics.facecoord[2][1]

    ψ  = @view ψδ[1:nn]
    δ =  @view ψδ[nn + 1 : 2nn]
    dψ = @view dψδ[1:nn]
    V = @view dψδ[nn + 1 : 2nn]

    mod_data_mms!(δ, ge, K, H̃, JH, vf, MMS, B_p, RS, metrics, t)

    u[:] = M \ ge

    HI = ops.HI[1]
    G = ops.G[1]
    Γ = ops.Γ[1]
    L = ops.L[1]
    sJ = metrics.sJ[1]

    Δτ[:] = - (HI * G * u + Crr * Γ * (δ ./ 2 - L*u)) ./ sJ

    for n in 1:nn

        ψn = ψ[n]
        bn = b[n]
        τn = Δτ[n]
        ηn = η[n]

        Vn = (2 * RS.V0 * sinh(τn / (RS.σn * RS.a))) / (exp(ψn/RS.a))

        if !isfinite(Vn)
            reject_step[1] = true
            return
        end

        V[n] = Vn

        if bn != 0
            dψ[n] = (bn * RS.V0 / RS.Dc) * (exp((RS.f0 - ψn) / bn) - abs(V[n]) / RS.V0)
            dψ[n] += S_rsh(xf1[n], yf1[n], t, bn, B_p, RS, MMS)[1]
        else
            dψ[n] = 0
        end

        if !isfinite(dψ[n])
            dψ .= 0
            reject_step[1] = true
            return
        end

    end

    nothing

end

function Q_DYNAMIC_MMS!(dψδ, ψδ, p, t)

    #@printf "%f\n" t
    @printf "\r\t%f" t/31556952

    reject_step = p.reject_step
    if reject_step[1]
        return
    end

    nn = p.nn
    Δτ = p.Δτ
    u = p.u
    ge = p.ge
    vf = p.vf
    M = p.ops.M̃
    K = p.ops.K
    H̃ = p.ops.H̃
    JH = p.ops.JH
    Crr = p.ops.Crr
    RS = p.RS
    MMS = p.MMS
    B_p = p.B_p
    metrics = p.metrics
    ops = p.ops
    η = metrics.η
    b = p.b
    f1_source = p.f1_source

    xf1 = metrics.facecoord[1][1]
    yf1 = metrics.facecoord[2][1]

    ψ  = @view ψδ[1:nn]
    δ =  @view ψδ[nn + 1 : 2nn]
    dψ = @view dψδ[1:nn]
    V = @view dψδ[nn + 1 : 2nn]

    mod_data_mms!(δ, ge, K, H̃, JH, vf, MMS, B_p, RS, metrics, t)

    u[:] = M \ ge

    HI = ops.HI[1]
    G = ops.G[1]
    Γ = ops.Γ[1]
    L = ops.L[1]
    sJ = metrics.sJ[1]

    Δτ[:] = - (HI * G * u + Crr * Γ * (δ ./ 2- L * u)) ./ sJ

    for n in 1:nn

        ψn = ψ[n]
        bn = b[n]
        τn = Δτ[n]
        ηn = η[n]

        VR = abs(τn / ηn)
        VL = -VR
        Vn = V[n]
        obj_rs(V) = rateandstateQ(V, ψn, RS.σn, τn, ηn, RS.a, RS.V0)
        (Vn, _, iter) = newtbndv(obj_rs, VL, VR, Vn; ftol = 1e-14,
                                 atolx = 1e-14, rtolx = 1e-14)

        if !isfinite(Vn)
            reject_step[1] = true
            return
        end

        V[n] = Vn

        if bn != 0
            dψ[n] = (bn * RS.V0 / RS.Dc) * (exp((RS.f0 - ψn) / bn) - abs(Vn) / RS.V0)
            dψ[n] += S_rsh(xf1[n], yf1[n], t, bn, B_p, RS, MMS)
        else
            dψ[n] = 0
        end

        if !isfinite(dψ[n])
            dψ .= 0
            reject_step[1] = true
            return
        end

    end

    nothing

end

# timestepping rejection function
function stepcheck(_, p, _)
    if p.reject_step[1]
        p.reject_step[1] = false
        return true
    end
    return false
end

function break_con(t_now, sim_seconds, cycle_flag, cycles, num_cycles)

    if cycle_flag == false
        return t_now < sim_seconds
    elseif cycle_flag == true
        return cycles <= num_cycles
    end

end

function PLOTFACE(ψδ,t,i)

    if isdefined(i,:fsallast)
        year_seconds = i.p.year_seconds
        yf1 = i.p.metrics.facecoord[2][1]
        xf1 = i.p.metrics.facecoord[1][1]
        MMS = i.p.MMS
        B_p = i.p.B_p
        nn = i.p.nn
        dψV = i.fsallast
        u = i.p.u
        L = i.p.ops.L[2]
        V = @view dψV[nn .+ (1:nn)]
        plot(V, yf1, legend=false, color =:blue, yflip=true)
        plot!(2*he_t(xf1, yf1, t, MMS), yf1, legend=false, color =:red, yflip=true)
        gui()
    end

    return false

end

# function for every accepted timstep with integrator stopping condition
function STOPFUN_Q(ψδ,t,i)

    if isdefined(i,:fsallast)

        dynamic_flag = i.p.dynamic_flag
        nn = i.p.nn
        δNp = i.p.δNp
        RS = i.p.RS
        τ = i.p.vars.Δτ[1:δNp]
        t_prev = i.p.vars.t_prev
        year_seconds = i.p.year_seconds
        u_prev = i.p.vars.u_prev
        u = i.p.vars.u
        fc = i.p.fc
        Lw = i.p.Lw
        io = i.p.io
        pf = i.p.io.pf
        η = i.p.metrics.η[1:δNp]
        sJ = i.p.metrics.sJ[1]
        cycles = i.p.cycles[1]
        HI = i.p.ops.HI[1]
        G = i.p.ops.G[1]

        dψV = i.fsallast
        ψ = @view ψδ[(1:δNp)]
        δ = @view ψδ[nn .+ (1:δNp)]
        V = @view dψV[nn .+ (1:δNp)]
        Vmax = maximum(abs.(V))

        if pf[1] % 40 == 0

            write_out_fault_data(io.remote_name,
                                 (i.p.ops.L[2] * u, zeros(nn)),
                                 0.0,
                                 t)

            write_out_fault_data(io.fault_name, (δ, V, τ .- η .* V, ψ), 0.0, t)

        end

        if pf[1] % 10 == 0

        write_out_stations(io.station_name,
                           io.stations,
                           fc,
                           (δ, V, τ .- η .* V, ψ),
                           maximum((η .* V) ./ (τ .- η .* V)),
                           t)

        end

        if io.slip_plot[1] != nothing && pf[1] % 120 == 0
            io.slip_plot[1] = plot!(io.slip_plot[1], δ, fc, linecolor=:blue, linewidth=.1)
            v_plot = plot(V, fc, legend = false, yflip=true, ylabel="Depth(Km)", xlabel="Slip rate (m/s)", color =:black)
            plot(io.slip_plot[1], v_plot, layout = (1,2))
            gui()
        end

        year_count = t/year_seconds

        if Vmax >= 1e-2  && dynamic_flag > 0
            return true
        end

        pf[1] += 1
        u_prev[:] = u
        t_prev[1] = t

    end

    return false

end

function MMS_WAVEPROP_CPU!(dq, q, p, t)

    #@printf "\r\t%f" t

    nn = p.nn
    fc = p.fc
    coord = p.coord
    R = p.R
    B_p = p.B_p
    MMS = p.MMS
    Λ = p.d_ops.Λ
    sJ = p.sJ
    Z̃f = p.d_ops.Z̃f
    L = p.d_ops.L
    H = p.d_ops.H
    P̃I = p.d_ops.P̃I
    JIHP = p.d_ops.JIHP
    CHAR_SOURCE = p.CHAR_SOURCE
    FORCE = p.FORCE

    u = q[1:nn^2]

    # compute all temporal derivatives
    dq[:] = Λ * q

    for i in 1:4
        fx = fc[1][i]
        fy = fc[2][i]
        S̃_c = sJ[i] .* CHAR_SOURCE(fx, fy, t, i, R[i], B_p, MMS)
        dq[2nn^2 + (i-1)*nn + 1 : 2nn^2 + i*nn] .+= S̃_c ./ (2*Z̃f[i])
        dq[nn^2 + 1:2nn^2] += L[i]' * H[i] * S̃_c ./ 2
    end
    dq[nn^2 + 1 : 2nn^2] = JIHP * dq[nn^2 + 1 : 2nn^2]
    dq[nn^2 + 1:2nn^2] += P̃I * FORCE(coord[1][:], coord[2][:], t, B_p, MMS)
end

function WAVEPROP!(dq, q, p, t)
    nn = p.nn
    Λ = p.d_ops.Λ
    JIHP = p.JIHP
    dq[:] = Λ * q
    dq[nn^2 + 1 : 2nn^2] = JIHP * dq[nn^2 + 1 : 2nn^2]
end

function MMS_FAULT_CPU!(dq, q, p, t)

    #@printf "\r\t%f" t

    nn = p.nn
    R = p.R
    fc = p.fc
    coord = p.coord
    B_p = p.B_p
    RS = p.RS
    b = p.b
    MMS = p.MMS
    Λ = p.d_ops.Λ
    sJ = p.sJ
    Z̃f = p.d_ops.Z̃f
    L = p.d_ops.L
    H = p.d_ops.H
    P̃I = p.d_ops.P̃I
    JIHP = p.d_ops.JIHP
    nCnΓ1 = p.d_ops.nCnΓ1
    HIGΓL1 = p.d_ops.HIGΓL1

    τ̃f = p.τ̃f

    u = @view q[1:nn^2]
    v = @view q[nn^2 + 1: 2nn^2]
    û1 = @view q[2nn^2 + 1 : 2nn^2 + nn]
    ψ = @view q[2nn^2 + 4nn + 1 : 2nn^2 + 5nn]
    dv = @view dq[nn^2 + 1: 2nn^2]
    dû1 = @view dq[2nn^2 + 1 : 2nn^2 + nn]
    dψ = @view dq[2nn^2 + 4nn + 1 : 2nn^2 + 5nn]

    # compute all temporal derivatives
    dq[:] = Λ * q

    # compute numerical traction on face 1
    #τ̃f .= sJ[1] .* τhe(fc[1][1], fc[2][1], t, 1, B_p, MMS)
    τ̃f[:] = HIGΓL1 * u + nCnΓ1 * û1
    #ψ .= ψe_d(fc[1][1], fc[2][1], t, B_p, RS, MMS)

    # Root find for RS friction
    for n in 1:nn

        vn = v[1 + nn * (n-1)]

        v̂_root(v̂) = rateandstateD(v̂,
                                  Z̃f[1][n],
                                  vn,
                                  sJ[1][n],
                                  ψ[n],
                                  RS.a,
                                  τ̃f[n],
                                  RS.σn,
                                  RS.V0)

        left = vn - τ̃f[n]/Z̃f[1][n]
        right = -left

        if left > right
            tmp = left
            left = right
            right = tmp
        end

        (v̂n, _, _) = newtbndv(v̂_root, left, right, vn; ftol = 1e-14,
                              atolx = 1e-14, rtolx = 1e-14)

        if isnan(v̂n)
            #println("Not bracketing root")
        end

        dû1[n] = v̂n
        #dû1[n] = he_t(fc[1][1][n], fc[2][1][n], t, MMS)
        dv[1 + (n - 1)*nn] +=  H[1][n, n] * (Z̃f[1][n] * v̂n)
        dψ[n] = (b[n] .* RS.V0 ./ RS.Dc) .* (exp.((RS.f0 .- ψ[n]) ./ b[n]) .- abs.(2 .* v̂n) ./ RS.V0)
    end

    # Non-fault Source injection
    for i in 2:4
        SOURCE = sJ[i] .* S_ch(fc[1][i], fc[2][i], t, i, R[i], B_p, MMS)
        dq[2nn^2 + (i-1)*nn + 1 : 2nn^2 + i*nn] .+= SOURCE ./ (2*Z̃f[i])
        dq[nn^2 + 1:2nn^2] .+= L[i]' * H[i] * SOURCE ./ 2
    end

    # psi source
    dψ .+= S_rsdh(fc[1][1], fc[2][1], b, t, B_p, RS, MMS)

    dq[nn^2 + 1:2nn^2] = JIHP * dq[nn^2 + 1:2nn^2]
    dq[nn^2 + 1:2nn^2] += P̃I * Forcing_hd(coord[1][:], coord[2][:], t, B_p, MMS)

end

function FAULT_CPU!(dq, q, p, t)

    nn = p.nn
    RS = p.RS
    b = p.b
    Λ = p.d_ops.Λ
    sJ = p.sJ[1]
    Z̃f = p.d_ops.Z̃f[1]
    H = p.d_ops.H[1]
    JIHP = p.d_ops.JIHP
    nCnΓ1 = p.d_ops.nCnΓ1
    nBBCΓL1 = p.d_ops.nBBCΓL1
    τ̃f = p.τ̃f

    u = @view q[1:nn^2]
    v = @view q[nn^2 + 1 : 2nn^2]
    û1 = @view q[2nn^2 + 1 : 2nn^2 + nn]
    ψ = @view q[2nn^2 + 4nn + 1 : 2nn^2 + 5nn]
    dv = @view dq[nn^2 + 1 : 2nn^2]
    dû1 = @view dq[2nn^2 + 1 : 2nn^2 + nn]
    dψ = @view dq[2nn^2 + 4nn + 1 : 2nn^2 + 5nn]

    dq[:] = Λ * q
    # compute numerical traction on face 1
    τ̃f[:] = nBBCΓL1 * u + nCnΓ1 * û1

    # Root find for RS friction
    for n in 1:nn

        vn = v[1 + nn * (n-1)]

        v̂_root(v̂) = rateandstateD(v̂,
                                  Z̃f[n],
                                  vn,
                                  sJ[n],
                                  ψ[n],
                                  RS.a,
                                  τ̃f[n],
                                  RS.σn,
                                  RS.V0)

        left = vn - τ̃f[n]/Z̃f[n]
        right = -left

        if left > right
            tmp = left
            left = right
            right = tmp
        end

        (v̂n, _, _) = newtbndv(v̂_root, left, right, vn; ftol = 1e-12,
                              atolx = 1e-12, rtolx = 1e-12)

        if isnan(v̂n)
            #println("Not bracketing root")
        end
        dû1[n] = v̂n
        dv[1 + (n - 1)*nn] +=  H[n, n] * (Z̃f[n] * v̂n)
        dψ[n] = (b[n] .* RS.V0 ./ RS.Dc) .* (exp.((RS.f0 .- ψ[n]) ./ b[n]) .- abs.(2 .* v̂n) ./ RS.V0)
    end

    dv[:] = JIHP * dq[nn^2 + 1:2nn^2]

end

#=
function FAULT_GPU!(dq, q, p, t)

    #@printf "\r\t%f" t

    nn = p.nn
    δNp = p.δNp
    RS = p.RS
    b = p.b
    Λ = p.Λ
    sJ = p.sJ
    Z̃f1 = p.Z̃f1
    Z̃f2 = p.Z̃f2
    Z̃f3 = p.Z̃f3
    L2 = p.L2
    L3 = p.L3
    H = p.H
    JIHP = p.JIHP
    nCnΓ1 = p.nCnΓ1
    HIGΓL1 = p.HIGΓL1
    τ̃f = p.τ̃f

    source2 = p.source2
    source3 = p.source3

    threads = p.threads
    blocks = p.blocks

    u = @view q[1:nn^2]
    vf = @view q[nn^2 + 1: nn : 2nn^2]
    û1 = @view q[2nn^2 + 1 : 2nn^2 + nn]
    ψ = @view q[2nn^2 + 4nn + 1 : 2nn^2 + 5nn]
    dv = @view dq[nn^2 + 1 : 2nn^2]
    dû1 = @view dq[2nn^2 + 1 : 2nn^2 + nn]
    dψ = @view dq[2nn^2 + 4nn + 1 : 2nn^2 + 5nn]
    dvf = @view dq[nn^2 + 1: nn : 2nn^2]

    dq[:] = Λ * q

    # compute numerical traction on face 1
    τ̃f[:] = HIGΓL1 * u + nCnΓ1 * û1

    @cuda blocks=blocks threads=threads FAULT_PROBLEM!(dû1, dvf, vf, τ̃f, Z̃f1, H, sJ, ψ, dψ, b, RS)

    dq[2nn^2 + nn + 1 : 2nn^2 + 2nn] += source2 ./ (2*Z̃f2)
    dq[nn^2 + 1:2nn^2] += L2' * (H .* source2 ./ 2)

    dq[2nn^2 + 2nn + 1 : 2nn^2 + 3nn] += source3 ./ (2*Z̃f3)
    dq[nn^2 + 1:2nn^2] += L3' * (H .* source3 ./ 2)

    dv[:] = JIHP * dq[nn^2 + 1:2nn^2]

end
=#

function FAULT_PROBLEM!(dû1, dvf, vf, τ̃f, Z̃f, H, sJ, ψ, dψ, b, RS)

    a = RS[1]
    σn = RS[2]
    V0 = RS[3]
    Dc = RS[4]
    f0 = RS[5]
    Vp = RS[6]
    nn = RS[7]
    δNp = RS[8]

    n = blockDim().x * (blockIdx().x - 1) + threadIdx().x

    if n <= nn

        if n <= δNp

            vn = vf[n]
            Z̃n = Z̃f[n]
            sJn = sJ[n]
            ψn = ψ[n]
            τ̃n = τ̃f[n]
            bn = b[n]
            Hn = H[n]
            v̂nL = vn - τ̃n/Z̃n
            v̂nR = -v̂nL

            (fL, _) = rateandstateD_GPU(v̂nL, Z̃n, vn, sJn, ψn, a, τ̃n, σn, V0)
            (fR, _) = rateandstateD_GPU(v̂nR, Z̃n, vn, sJn, ψn, a, τ̃n, σn, V0)
            (f, df) = rateandstateD_GPU(vn, Z̃n, vn, sJn, ψn, a, τ̃n, σn, V0)

            dv̂nlr = v̂nR - v̂nL

            v̂n = vn

            count = 0
            for iter in 1:500

                dv̂n = -f / df
                v̂n  = v̂n + dv̂n

                if v̂n < v̂nL || v̂n > v̂nR || abs(dv̂n) / dv̂nlr < 0
                    v̂n = (v̂nR + v̂nL) / 2
                    dv̂n = (v̂nR - v̂nL) / 2
                end

                (f, df) = rateandstateD_GPU(v̂n, Z̃n, vn, sJn, ψn, a, τ̃n, σn, V0)

                if f * fL > 0
                    (fL, v̂nL) = (f, v̂n)
                else
                    (fR, v̂nR) = (f, v̂n)
                end
                dv̂nlr = v̂nR - v̂nL

                if abs(f) < 1e-12 && abs(dv̂n) < 1e-12 + 1e-12 * (abs(dv̂n) + abs(v̂n))
                    break
                end
                count += 1
            end

            dû1[n] = v̂n
            dvf[n] += Hn * (Z̃n * v̂n)
            dψ[n] = (bn .* V0 ./ Dc) .* (exp.((f0 .- ψn) ./ bn) .- abs.(2 .* v̂n) ./ V0)

        else

            Hn = H[n]
            τ̃n = τ̃f[n]

            dû1[n] = Vp/2
            dvf[n] += Hn * τ̃n

        end

    end

    nothing

end

# bracketed newton method
function newtbndv(func, xL, xR, x; ftol = 1e-6, maxiter = 500, minchange=0,
                  atolx = 1e-4, rtolx = 1e-4)

    (fL, _) = func(xL)
    (fR, _) = func(xR)
    if fL .* fR > 0
        return (typeof(x)(NaN), typeof(x)(NaN), -maxiter)
    end

    (f, df) = func(x)
    dxlr = xR - xL

    for iter = 1:maxiter
        dx = -f / df
        x  = x + dx

        if x < xL || x > xR || abs(dx) / dxlr < minchange
            x = (xR + xL) / 2
            dx = (xR - xL) / 2
        end

        (f, df) = func(x)

        if f * fL > 0
            (fL, xL) = (f, x)
        else
            (fR, xR) = (f, x)
        end
        dxlr = xR - xL

        if abs(f) < ftol && abs(dx) < atolx + rtolx * (abs(dx) + abs(x))
            return (x, f, iter)
        end
    end
    return (x, f, -maxiter)
end

# Dynamic roofinding problem on the fault
function rateandstateD(v̂, z̃, v, sJ, ψ, a, τ̃, σn, V0)
    Y = (1 / (2 * V0)) * exp(ψ / a)
    f = a * asinh(2v̂ * Y)
    dfdv̂  = a * (1 / sqrt(1 + (2v̂ * Y)^2)) * 2Y
    g = sJ * σn * f + τ̃ + z̃*(v̂ - v)
    dgdv̂ = z̃ + sJ * σn * dfdv̂

    return (g, dgdv̂)
end

# Dynamic roofinding problem on the fault
function rateandstateD_GPU(v̂, z̃, v, sJ, ψ, a, τ̃, σn, V0)
    Y = (1 / (2 * V0)) * exp(ψ / a)
    f = a * CUDA.asinh(2v̂ * Y)
    dfdv̂  = a * (1 / sqrt(1 + (2v̂ * Y)^2)) * 2Y
    g = sJ * σn * f + τ̃ + z̃*(v̂ - v)
    dgdv̂ = z̃ + sJ * σn * dfdv̂

    return (g, dgdv̂)
end

function timestep_write!(q, f!, p, dt, (t0, t1), Δq = similar(q), Δq2 = similar(q))
    T = eltype(q)

    nn = p.nn
    δNp = p.δNp
    Lw = p.Lw
    fc = p.fc
    pf = p.io.pf
    τ̃f = p.τ̃f
    sJ = p.sJ
    Z̃f = p. Z̃f1
    io = p.io
    v̂ = p.v̂
    δ = p.δ
    v̂_cpu = p.v̂_cpu
    τ̂_cpu = p.τ̂_cpu
    τ̃_cpu = p.τ̃_cpu
    ψ_cpu = p.ψ_cpu
    η = p.η
    d_to_s = p.d_to_s
    HIG = p.HIG
    vf = @view q[nn^2 + 1: nn : 2nn^2]
    v = @view q[nn^2 + 1 : 2nn^2]
    u = @view q[1 : nn^2]
    ûf = @view q[2nn^2 + 1 : 2nn^2 + nn]
    ψ = @view q[2nn^2 + 4nn + 1 : 2nn^2 + 5*nn]

    RKA = [
        T(0),
        T(-567301805773 // 1357537059087),
        T(-2404267990393 // 2016746695238),
        T(-3550918686646 // 2091501179385),
        T(-1275806237668 // 842570457699),
    ]

    RKB = [
        T(1432997174477 // 9575080441755),
        T(5161836677717 // 13612068292357),
        T(1720146321549 // 2090206949498),
        T(3134564353537 // 4481467310338),
        T(2277821191437 // 14882151754819),
    ]

    RKC = [
        T(0),
        T(1432997174477 // 9575080441755),
        T(2526269341429 // 6820363962896),
        T(2006345519317 // 3224310063776),
        T(2802321613138 // 2924317926251),
    ]

    nstep = ceil(Int, (t1 - t0) / dt)
    dt = (t1 - t0) / nstep

    pf[1] = .1
    pf[2] = .01
    pf[3] = 1.0

    fill!(Δq, 0)
    fill!(Δq2, 0)
    for step in 1:nstep
        t = t0 + (step - 1) * dt
        for s in 1:length(RKA)
            f!(Δq2, q, p, t + RKC[s] * dt)
            v̂[:] = Δq2[2nn^2 + 1 : 2nn^2 + nn]
            Δq .+= Δq2
            q .+= (RKB[s] * dt) .* Δq
            Δq .*= RKA[s % length(RKA) + 1]
        end

        v̂_cpu[:] = Array(v̂)

        if step == ceil(Int, pf[1]/dt)

            δ[:] = Array(2ûf)
            τ̂_cpu[:] = Array(-τ̃f ./ sJ .- Z̃f .* (v̂ - vf) ./ sJ)
            ψ_cpu[:] = Array(ψ)

            if any(isnan, v̂_cpu)
                error("nan from dynamic rootfinder")
            end

            write_out_fault_data(io.fault_name,
                                 (δ[1:δNp], 2v̂_cpu[1:δNp], τ̂_cpu[1:δNp], ψ_cpu[1:δNp]),
                                 0.0,
                                 t)

            write_out_fault_data(io.remote_name,
                                 (Array(p.L2 * u), zeros(nn)),
                                 0.0,
                                 t)

            if io.vp == 1
                u_out = Array(u)
                v_out = Array(v)
                write_out_volume(io.volume_name, (u_out, v_out), 2v̂_cpu, nn, t)
            end

            pf[1] +=.1
        end

        if step == ceil(Int, pf[2]/dt)

            δ[:] = Array(2ûf)
            τ̂_cpu[:] = Array(-τ̃f ./ sJ .- Z̃f .* (v̂ - vf) ./ sJ)
            ψ_cpu[:] = Array(ψ)

            write_out_stations(io.station_name,
                               io.stations,
                               fc,
                               (δ, 2v̂_cpu, τ̂_cpu, ψ_cpu[1:δNp]),
                               maximum((η[1:δNp] .* 2v̂_cpu[1:δNp]) ./ τ̂_cpu[1:δNp]),
                               t)

        pf[2] += .01
        end

        if step == ceil(Int, pf[3]/dt)
           if io.slip_plot[1] != nothing
                io.slip_plot[1] = plot!(io.slip_plot[1], δ[1:δNp], fc[1:δNp], linecolor=:red, linewidth=.1)
                v_plot = plot(2v̂_cpu[1:δNp], fc[1:δNp], legend = false, yflip=true, ylabel="Depth(Km)", xlabel="Slip rate (m/s)", color =:black)
                plot(io.slip_plot[1], v_plot, layout = (1,2))
                gui()
            end
            pf[3] += 1.0
        end

        if maximum(v̂_cpu) < d_to_s/2
            return t
        end

    end

    nothing

end

function timestep!(q, f!, p, dt, (t0, t1), Δq = similar(q), Δq2 = similar(q))
    T = eltype(q)

    RKA = (
        T(0),
        T(-567301805773 // 1357537059087),
        T(-2404267990393 // 2016746695238),
        T(-3550918686646 // 2091501179385),
        T(-1275806237668 // 842570457699),
    )

    RKB = (
        T(1432997174477 // 9575080441755),
        T(5161836677717 // 13612068292357),
        T(1720146321549 // 2090206949498),
        T(3134564353537 // 4481467310338),
        T(2277821191437 // 14882151754819),
    )

    RKC = (
        T(0),
        T(1432997174477 // 9575080441755),
        T(2526269341429 // 6820363962896),
        T(2006345519317 // 3224310063776),
        T(2802321613138 // 2924317926251),
    )

    nstep = ceil(Int, (t1 - t0) / dt)
    dt = (t1 - t0) / nstep

    fill!(Δq, 0)
    fill!(Δq2, 0)
    for step in 1:nstep
        t = t0 + (step - 1) * dt
        for s in 1:length(RKA)
            f!(Δq2, q, p, t + RKC[s] * dt)
            Δq .+= Δq2
            q .+= RKB[s] * dt * Δq
            Δq .*= RKA[s % length(RKA) + 1]
        end
    end

    nothing

end

function euler!(q, f!, p, dt, t_span, Δq = similar(q))

    nstep = ceil(Int, (t_span[2] - t_span[1]) / dt)
    dt = (t_span[2] - t_span[1]) / nstep

    for step in 1:nstep
        t = t_span[1] + (step - 1) * dt
        fill!(Δq, 0)
        f!(Δq, q, p, t)
        q .+= dt * Δq
    end

end

# Quasi-Dynamic rootfinding problem on the fault
function rateandstateQ(V, ψ, σn, τn, ηn, a, V0)

    Y = (1 ./ (2 .* V0)) .* exp.(ψ ./ a)
    f = a .* asinh.(V .* Y)
    dfdV  = a .* (1 ./ sqrt.(1 + (V .* Y).^2)) .* Y
    g = σn .* f + ηn .* V - τn
    dgdV = σn .* dfdV + ηn
    (g, dgdV)
end


text:
mutable struct Player
    score::Int
    ante::Int
    wins::Int
    losses::Int
    strategy::Pair{String, Function}
end

randomchoicetostop(player, group) = rand(Bool)
variablerandtostop(player, group) = any(x -> x.score > player.score, group) ? rand() < 0.1 : rand(Bool)
overtwentystop(player, group) = player.ante > 20
over20unlesslosingstop(player, group) = player.ante > 20 && all(x -> x.score < 80, group)

const strategies = ("random choice to stop" => randomchoicetostop, "variable rand to stop" => variablerandtostop,
                  "roll to 20" => overtwentystop, "roll to 20 then if not losing stop" => over20unlesslosingstop)
const players = [Player(0, 0, 0, 0, s) for s in strategies]
const dice = collect(1:6)

function turn(player, verbose=false)
    playernum = findfirst(p -> p == player, players)
    scorewin() = for p in players if p == player p.wins += 1 else p.losses += 1 end; p.score = 0 end
    player.ante = 0
    while (r = rand(dice)) != 1
        player.ante += r
        verbose && println("Player $playernum rolls a $r.")
        if player.score + player.ante >= 100
            scorewin()
            verbose && println("Player $playernum wins.\n")
            return false
        elseif player.strategy[2](player, players)
            player.score += player.ante
            verbose && println("Player $playernum holds and has a new score of $(player.score).")
            return true
        end
    end
    verbose && println("Player $playernum rolls a 1, so turn is over.")
    true
end

function rungames(N)
    for i in 1:N
        verbose = (i == 3) ? true : false  # do verbose if it's game number 3
        curplayer = rand(collect(1:length(players)))
        while turn(players[curplayer], verbose)
            curplayer = curplayer >= length(players) ? 1 : curplayer + 1
        end
    end
    results = sort([(p.wins/(p.wins + p.losses), p.strategy[1]) for p in players], rev=true)
    println("             Strategy                % of wins (N = $N)")
    println("------------------------------------------------------------")
    for pair in results
        println(lpad(pair[2], 34), lpad(round(pair[1] * 100, digits=1), 18))
    end
end

rungames(1000000)

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Parameters

export hamiltonian, ϑ, ω, ωabs

@inline function u(t, q)
    q[4]
end

function hamiltonian(t, q, params)
    @unpack μ = params
    0.5 * u(t,q)^2 + μ*B(t,q)
end

dHdx₁(t,q,μ) = μ * dBdx₁(t,q)
dHdx₂(t,q,μ) = μ * dBdx₂(t,q)
dHdx₃(t,q,μ) = μ * dBdx₃(t,q)
dHdx₄(t,q,μ) = u(t,q)

function dH(t, q, dH, params)
    @unpack μ = params
    dH[1] = dHdx₁(t,q,μ)
    dH[2] = dHdx₂(t,q,μ)
    dH[3] = dHdx₃(t,q,μ)
    dH[4] = dHdx₄(t,q,μ)
    nothing
end

d²Hdx₁dx₁(t, q, μ) = μ * d²Bdx₁dx₁(t,q)
d²Hdx₁dx₂(t, q, μ) = μ * d²Bdx₁dx₂(t,q)
d²Hdx₁dx₃(t, q, μ) = μ * d²Bdx₁dx₃(t,q)
d²Hdx₁dx₄(t, q, μ) = zero(eltype(q))
d²Hdx₂dx₁(t, q, μ) = μ * d²Bdx₂dx₁(t,q)
d²Hdx₂dx₂(t, q, μ) = μ * d²Bdx₂dx₂(t,q)
d²Hdx₂dx₃(t, q, μ) = μ * d²Bdx₂dx₃(t,q)
d²Hdx₂dx₄(t, q, μ) = zero(eltype(q))
d²Hdx₃dx₁(t, q, μ) = μ * d²Bdx₃dx₁(t,q)
d²Hdx₃dx₂(t, q, μ) = μ * d²Bdx₃dx₂(t,q)
d²Hdx₃dx₃(t, q, μ) = μ * d²Bdx₃dx₃(t,q)
d²Hdx₃dx₄(t, q, μ) = zero(eltype(q))
d²Hdx₄dx₁(t, q, μ) = zero(eltype(q))
d²Hdx₄dx₂(t, q, μ) = zero(eltype(q))
d²Hdx₄dx₃(t, q, μ) = zero(eltype(q))
d²Hdx₄dx₄(t, q, μ) = one(eltype(q))

ϑ₁(t,q) = A₁(t,q) + u(t,q) * b₁(t,q)
ϑ₂(t,q) = A₂(t,q) + u(t,q) * b₂(t,q)
ϑ₃(t,q) = A₃(t,q) + u(t,q) * b₃(t,q)
ϑ₄(t,q) = zero(eltype(q))

function ϑ(t, q, ϑ, params)
    ϑ[1] = ϑ₁(t,q)
    ϑ[2] = ϑ₂(t,q)
    ϑ[3] = ϑ₃(t,q)
    ϑ[4] = ϑ₄(t,q)
    nothing
end

dϑ₁dx₁(t,q) = dA₁dx₁(t,q) + u(t,q) * db₁dx₁(t,q)
dϑ₁dx₂(t,q) = dA₁dx₂(t,q) + u(t,q) * db₁dx₂(t,q)
dϑ₁dx₃(t,q) = dA₁dx₃(t,q) + u(t,q) * db₁dx₃(t,q)
dϑ₁dx₄(t,q) = b₁(t,q)
dϑ₂dx₁(t,q) = dA₂dx₁(t,q) + u(t,q) * db₂dx₁(t,q)
dϑ₂dx₂(t,q) = dA₂dx₂(t,q) + u(t,q) * db₂dx₂(t,q)
dϑ₂dx₃(t,q) = dA₂dx₃(t,q) + u(t,q) * db₂dx₃(t,q)
dϑ₂dx₄(t,q) = b₂(t,q)
dϑ₃dx₁(t,q) = dA₃dx₁(t,q) + u(t,q) * db₃dx₁(t,q)
dϑ₃dx₂(t,q) = dA₃dx₂(t,q) + u(t,q) * db₃dx₂(t,q)
dϑ₃dx₃(t,q) = dA₃dx₃(t,q) + u(t,q) * db₃dx₃(t,q)
dϑ₃dx₄(t,q) = b₃(t,q)
dϑ₄dx₁(t,q) = zero(eltype(q))
dϑ₄dx₂(t,q) = zero(eltype(q))
dϑ₄dx₃(t,q) = zero(eltype(q))
dϑ₄dx₄(t,q) = zero(eltype(q))

d²ϑ₁dx₁dx₁(t,q) = d²A₁dx₁dx₁(t,q) + u(t,q) * d²b₁dx₁dx₁(t,q)
d²ϑ₁dx₁dx₂(t,q) = d²A₁dx₁dx₂(t,q) + u(t,q) * d²b₁dx₁dx₂(t,q)
d²ϑ₁dx₁dx₃(t,q) = d²A₁dx₁dx₃(t,q) + u(t,q) * d²b₁dx₁dx₃(t,q)
d²ϑ₁dx₁dx₄(t,q) = db₁dx₁(t,q)

d²ϑ₁dx₂dx₁(t,q) = d²A₁dx₂dx₁(t,q) + u(t,q) * d²b₁dx₂dx₁(t,q)
d²ϑ₁dx₂dx₂(t,q) = d²A₁dx₂dx₂(t,q) + u(t,q) * d²b₁dx₂dx₂(t,q)
d²ϑ₁dx₂dx₃(t,q) = d²A₁dx₂dx₃(t,q) + u(t,q) * d²b₁dx₂dx₃(t,q)
d²ϑ₁dx₂dx₄(t,q) = db₁dx₂(t,q)

d²ϑ₁dx₃dx₁(t,q) = d²A₁dx₃dx₁(t,q) + u(t,q) * d²b₁dx₃dx₁(t,q)
d²ϑ₁dx₃dx₂(t,q) = d²A₁dx₃dx₂(t,q) + u(t,q) * d²b₁dx₃dx₂(t,q)
d²ϑ₁dx₃dx₃(t,q) = d²A₁dx₃dx₃(t,q) + u(t,q) * d²b₁dx₃dx₃(t,q)
d²ϑ₁dx₃dx₄(t,q) = db₁dx₃(t,q)

d²ϑ₁dx₄dx₁(t,q) = db₁dx₁(t,q)
d²ϑ₁dx₄dx₂(t,q) = db₁dx₂(t,q)
d²ϑ₁dx₄dx₃(t,q) = db₁dx₃(t,q)
d²ϑ₁dx₄dx₄(t,q) = zero(eltype(q))

d²ϑ₂dx₁dx₁(t,q) = d²A₂dx₁dx₁(t,q) + u(t,q) * d²b₂dx₁dx₁(t,q)
d²ϑ₂dx₁dx₂(t,q) = d²A₂dx₁dx₂(t,q) + u(t,q) * d²b₂dx₁dx₂(t,q)
d²ϑ₂dx₁dx₃(t,q) = d²A₂dx₁dx₃(t,q) + u(t,q) * d²b₂dx₁dx₃(t,q)
d²ϑ₂dx₁dx₄(t,q) = db₂dx₁(t,q)

d²ϑ₂dx₂dx₁(t,q) = d²A₂dx₂dx₁(t,q) + u(t,q) * d²b₂dx₂dx₁(t,q)
d²ϑ₂dx₂dx₂(t,q) = d²A₂dx₂dx₂(t,q) + u(t,q) * d²b₂dx₂dx₂(t,q)
d²ϑ₂dx₂dx₃(t,q) = d²A₂dx₂dx₃(t,q) + u(t,q) * d²b₂dx₂dx₃(t,q)
d²ϑ₂dx₂dx₄(t,q) = db₂dx₂(t,q)

d²ϑ₂dx₃dx₁(t,q) = d²A₂dx₃dx₁(t,q) + u(t,q) * d²b₂dx₃dx₁(t,q)
d²ϑ₂dx₃dx₂(t,q) = d²A₂dx₃dx₂(t,q) + u(t,q) * d²b₂dx₃dx₂(t,q)
d²ϑ₂dx₃dx₃(t,q) = d²A₂dx₃dx₃(t,q) + u(t,q) * d²b₂dx₃dx₃(t,q)
d²ϑ₂dx₃dx₄(t,q) = db₂dx₃(t,q)

d²ϑ₂dx₄dx₁(t,q) = db₂dx₁(t,q)
d²ϑ₂dx₄dx₂(t,q) = db₂dx₂(t,q)
d²ϑ₂dx₄dx₃(t,q) = db₂dx₃(t,q)
d²ϑ₂dx₄dx₄(t,q) = zero(eltype(q))

d²ϑ₃dx₁dx₁(t,q) = d²A₃dx₁dx₁(t,q) + u(t,q) * d²b₃dx₁dx₁(t,q)
d²ϑ₃dx₁dx₂(t,q) = d²A₃dx₁dx₂(t,q) + u(t,q) * d²b₃dx₁dx₂(t,q)
d²ϑ₃dx₁dx₃(t,q) = d²A₃dx₁dx₃(t,q) + u(t,q) * d²b₃dx₁dx₃(t,q)
d²ϑ₃dx₁dx₄(t,q) = db₃dx₁(t,q)

d²ϑ₃dx₂dx₁(t,q) = d²A₃dx₂dx₁(t,q) + u(t,q) * d²b₃dx₂dx₁(t,q)
d²ϑ₃dx₂dx₂(t,q) = d²A₃dx₂dx₂(t,q) + u(t,q) * d²b₃dx₂dx₂(t,q)
d²ϑ₃dx₂dx₃(t,q) = d²A₃dx₂dx₃(t,q) + u(t,q) * d²b₃dx₂dx₃(t,q)
d²ϑ₃dx₂dx₄(t,q) = db₃dx₂(t,q)

d²ϑ₃dx₃dx₁(t,q) = d²A₃dx₃dx₁(t,q) + u(t,q) * d²b₃dx₃dx₁(t,q)
d²ϑ₃dx₃dx₂(t,q) = d²A₃dx₃dx₂(t,q) + u(t,q) * d²b₃dx₃dx₂(t,q)
d²ϑ₃dx₃dx₃(t,q) = d²A₃dx₃dx₃(t,q) + u(t,q) * d²b₃dx₃dx₃(t,q)
d²ϑ₃dx₃dx₄(t,q) = db₃dx₃(t,q)

d²ϑ₃dx₄dx₁(t,q) = db₃dx₁(t,q)
d²ϑ₃dx₄dx₂(t,q) = db₃dx₂(t,q)
d²ϑ₃dx₄dx₃(t,q) = db₃dx₃(t,q)
d²ϑ₃dx₄dx₄(t,q) = zero(eltype(q))

d²ϑ₄dx₁dx₁(t,q) = zero(eltype(q))
d²ϑ₄dx₁dx₂(t,q) = zero(eltype(q))
d²ϑ₄dx₁dx₃(t,q) = zero(eltype(q))
d²ϑ₄dx₁dx₄(t,q) = zero(eltype(q))

d²ϑ₄dx₂dx₁(t,q) = zero(eltype(q))
d²ϑ₄dx₂dx₂(t,q) = zero(eltype(q))
d²ϑ₄dx₂dx₃(t,q) = zero(eltype(q))
d²ϑ₄dx₂dx₄(t,q) = zero(eltype(q))

d²ϑ₄dx₃dx₁(t,q) = zero(eltype(q))
d²ϑ₄dx₃dx₂(t,q) = zero(eltype(q))
d²ϑ₄dx₃dx₃(t,q) = zero(eltype(q))
d²ϑ₄dx₃dx₄(t,q) = zero(eltype(q))

d²ϑ₄dx₄dx₁(t,q) = zero(eltype(q))
d²ϑ₄dx₄dx₂(t,q) = zero(eltype(q))
d²ϑ₄dx₄dx₃(t,q) = zero(eltype(q))
d²ϑ₄dx₄dx₄(t,q) = zero(eltype(q))

ω₁(t,q) = dϑ₃dx₂(t,q) - dϑ₂dx₃(t,q)
ω₂(t,q) = dϑ₁dx₃(t,q) - dϑ₃dx₁(t,q)
ω₃(t,q) = dϑ₂dx₁(t,q) - dϑ₁dx₂(t,q)

function ω(t, q, ω::Vector, params=NamedTuple())
    ω[1] = ω₁(t,q)
    ω[2] = ω₂(t,q)
    ω[3] = ω₃(t,q)
    nothing
end

function ω(t, q, Ω::Matrix, params=NamedTuple())
    Ω[1,1] = 0
    Ω[1,2] = dϑ₁dx₂(t,q) - dϑ₂dx₁(t,q)
    Ω[1,3] = dϑ₁dx₃(t,q) - dϑ₃dx₁(t,q)
    Ω[1,4] = dϑ₁dx₄(t,q) - dϑ₄dx₁(t,q)

    Ω[2,1] = dϑ₂dx₁(t,q) - dϑ₁dx₂(t,q)
    Ω[2,2] = 0
    Ω[2,3] = dϑ₂dx₃(t,q) - dϑ₃dx₂(t,q)
    Ω[2,4] = dϑ₂dx₄(t,q) - dϑ₄dx₂(t,q)

    Ω[3,1] = dϑ₃dx₁(t,q) - dϑ₁dx₃(t,q)
    Ω[3,2] = dϑ₃dx₂(t,q) - dϑ₂dx₃(t,q)
    Ω[3,3] = 0
    Ω[3,4] = dϑ₃dx₄(t,q) - dϑ₄dx₃(t,q)

    Ω[4,1] = dϑ₄dx₁(t,q) - dϑ₁dx₄(t,q)
    Ω[4,2] = dϑ₄dx₂(t,q) - dϑ₂dx₄(t,q)
    Ω[4,3] = dϑ₄dx₃(t,q) - dϑ₃dx₄(t,q)
    Ω[4,4] = 0

    nothing
end

function ωabs(t, q, params=NamedTuple())
    ω₁(t,q) * dϑ₁dx₄(t,q) + ω₂(t,q) * dϑ₂dx₄(t,q) + ω₃(t,q) * dϑ₃dx₄(t,q)
end

dω₁dx₁(t,q) = d²ϑ₃dx₂dx₁(t,q) - d²ϑ₂dx₃dx₁(t,q)
dω₁dx₂(t,q) = d²ϑ₃dx₂dx₂(t,q) - d²ϑ₂dx₃dx₂(t,q)
dω₁dx₃(t,q) = d²ϑ₃dx₂dx₃(t,q) - d²ϑ₂dx₃dx₃(t,q)
dω₁dx₄(t,q) = d²ϑ₃dx₂dx₄(t,q) - d²ϑ₂dx₃dx₄(t,q)
dω₂dx₁(t,q) = d²ϑ₁dx₃dx₁(t,q) - d²ϑ₃dx₁dx₁(t,q)
dω₂dx₂(t,q) = d²ϑ₁dx₃dx₂(t,q) - d²ϑ₃dx₁dx₂(t,q)
dω₂dx₃(t,q) = d²ϑ₁dx₃dx₃(t,q) - d²ϑ₃dx₁dx₃(t,q)
dω₂dx₄(t,q) = d²ϑ₁dx₃dx₄(t,q) - d²ϑ₃dx₁dx₄(t,q)
dω₃dx₁(t,q) = d²ϑ₂dx₁dx₁(t,q) - d²ϑ₁dx₂dx₁(t,q)
dω₃dx₂(t,q) = d²ϑ₂dx₁dx₂(t,q) - d²ϑ₁dx₂dx₂(t,q)
dω₃dx₃(t,q) = d²ϑ₂dx₁dx₃(t,q) - d²ϑ₁dx₂dx₃(t,q)
dω₃dx₄(t,q) = d²ϑ₂dx₁dx₄(t,q) - d²ϑ₁dx₂dx₄(t,q)

β₁(t,q,μ) = dHdx₄(t,q,μ) * ϑ₁(t,q)
β₂(t,q,μ) = dHdx₄(t,q,μ) * ϑ₂(t,q)
β₃(t,q,μ) = dHdx₄(t,q,μ) * ϑ₃(t,q)

function β(t, q, β, params)
    @unpack μ = params
    β[1] = β₁(t,q,μ)
    β[2] = β₂(t,q,μ)
    β[3] = β₃(t,q,μ)
    nothing
end

γ₁(t,q,μ) = ϑ₂(t,q) * dHdx₃(t,q,μ) - dHdx₂(t,q,μ) * ϑ₃(t,q)
γ₂(t,q,μ) = ϑ₃(t,q) * dHdx₁(t,q,μ) - dHdx₃(t,q,μ) * ϑ₁(t,q)
γ₃(t,q,μ) = ϑ₁(t,q) * dHdx₂(t,q,μ) - dHdx₁(t,q,μ) * ϑ₂(t,q)

function γ(t, q, γ, params)
    @unpack μ = params
    γ[1] = γ₁(t,q,μ)
    γ[2] = γ₂(t,q,μ)
    γ[3] = γ₃(t,q,μ)
    nothing
end

dβ₁dx₂(t,q,μ) = d²Hdx₄dx₂(t,q,μ) * ϑ₁(t,q) + dHdx₄(t,q,μ) * dϑ₁dx₂(t,q)
dβ₁dx₁(t,q,μ) = d²Hdx₄dx₁(t,q,μ) * ϑ₁(t,q) + dHdx₄(t,q,μ) * dϑ₁dx₁(t,q)
dβ₁dx₃(t,q,μ) = d²Hdx₄dx₃(t,q,μ) * ϑ₁(t,q) + dHdx₄(t,q,μ) * dϑ₁dx₃(t,q)
dβ₁dx₄(t,q,μ) = d²Hdx₄dx₄(t,q,μ) * ϑ₁(t,q) + dHdx₄(t,q,μ) * dϑ₁dx₄(t,q)
dβ₂dx₁(t,q,μ) = d²Hdx₄dx₁(t,q,μ) * ϑ₂(t,q) + dHdx₄(t,q,μ) * dϑ₂dx₁(t,q)
dβ₂dx₂(t,q,μ) = d²Hdx₄dx₂(t,q,μ) * ϑ₂(t,q) + dHdx₄(t,q,μ) * dϑ₂dx₂(t,q)
dβ₂dx₃(t,q,μ) = d²Hdx₄dx₃(t,q,μ) * ϑ₂(t,q) + dHdx₄(t,q,μ) * dϑ₂dx₃(t,q)
dβ₂dx₄(t,q,μ) = d²Hdx₄dx₄(t,q,μ) * ϑ₂(t,q) + dHdx₄(t,q,μ) * dϑ₂dx₄(t,q)
dβ₃dx₁(t,q,μ) = d²Hdx₄dx₁(t,q,μ) * ϑ₃(t,q) + dHdx₄(t,q,μ) * dϑ₃dx₁(t,q)
dβ₃dx₂(t,q,μ) = d²Hdx₄dx₂(t,q,μ) * ϑ₃(t,q) + dHdx₄(t,q,μ) * dϑ₃dx₂(t,q)
dβ₃dx₃(t,q,μ) = d²Hdx₄dx₃(t,q,μ) * ϑ₃(t,q) + dHdx₄(t,q,μ) * dϑ₃dx₃(t,q)
dβ₃dx₄(t,q,μ) = d²Hdx₄dx₄(t,q,μ) * ϑ₃(t,q) + dHdx₄(t,q,μ) * dϑ₃dx₄(t,q)

dγ₁dx₁(t,q,μ) = dϑ₂dx₁(t,q) * dHdx₃(t,q,μ) + ϑ₂(t,q) * d²Hdx₃dx₁(t,q,μ) - d²Hdx₂dx₁(t,q,μ) * ϑ₃(t,q) - dHdx₂(t,q,μ) * dϑ₃dx₁(t,q)
dγ₁dx₂(t,q,μ) = dϑ₂dx₂(t,q) * dHdx₃(t,q,μ) + ϑ₂(t,q) * d²Hdx₃dx₂(t,q,μ) - d²Hdx₂dx₂(t,q,μ) * ϑ₃(t,q) - dHdx₂(t,q,μ) * dϑ₃dx₂(t,q)
dγ₁dx₃(t,q,μ) = dϑ₂dx₃(t,q) * dHdx₃(t,q,μ) + ϑ₂(t,q) * d²Hdx₃dx₃(t,q,μ) - d²Hdx₂dx₃(t,q,μ) * ϑ₃(t,q) - dHdx₂(t,q,μ) * dϑ₃dx₃(t,q)
dγ₁dx₄(t,q,μ) = dϑ₂dx₄(t,q) * dHdx₃(t,q,μ) + ϑ₂(t,q) * d²Hdx₃dx₄(t,q,μ) - d²Hdx₂dx₄(t,q,μ) * ϑ₃(t,q) - dHdx₂(t,q,μ) * dϑ₃dx₄(t,q)
dγ₂dx₁(t,q,μ) = dϑ₃dx₁(t,q) * dHdx₁(t,q,μ) + ϑ₃(t,q) * d²Hdx₁dx₁(t,q,μ) - d²Hdx₃dx₁(t,q,μ) * ϑ₁(t,q) - dHdx₃(t,q,μ) * dϑ₁dx₁(t,q)
dγ₂dx₂(t,q,μ) = dϑ₃dx₂(t,q) * dHdx₁(t,q,μ) + ϑ₃(t,q) * d²Hdx₁dx₂(t,q,μ) - d²Hdx₃dx₂(t,q,μ) * ϑ₁(t,q) - dHdx₃(t,q,μ) * dϑ₁dx₂(t,q)
dγ₂dx₃(t,q,μ) = dϑ₃dx₃(t,q) * dHdx₁(t,q,μ) + ϑ₃(t,q) * d²Hdx₁dx₃(t,q,μ) - d²Hdx₃dx₃(t,q,μ) * ϑ₁(t,q) - dHdx₃(t,q,μ) * dϑ₁dx₃(t,q)
dγ₂dx₄(t,q,μ) = dϑ₃dx₄(t,q) * dHdx₁(t,q,μ) + ϑ₃(t,q) * d²Hdx₁dx₄(t,q,μ) - d²Hdx₃dx₄(t,q,μ) * ϑ₁(t,q) - dHdx₃(t,q,μ) * dϑ₁dx₄(t,q)
dγ₃dx₁(t,q,μ) = dϑ₁dx₁(t,q) * dHdx₂(t,q,μ) + ϑ₁(t,q) * d²Hdx₂dx₁(t,q,μ) - d²Hdx₁dx₁(t,q,μ) * ϑ₂(t,q) - dHdx₁(t,q,μ) * dϑ₂dx₁(t,q)
dγ₃dx₂(t,q,μ) = dϑ₁dx₂(t,q) * dHdx₂(t,q,μ) + ϑ₁(t,q) * d²Hdx₂dx₂(t,q,μ) - d²Hdx₁dx₂(t,q,μ) * ϑ₂(t,q) - dHdx₁(t,q,μ) * dϑ₂dx₂(t,q)
dγ₃dx₃(t,q,μ) = dϑ₁dx₃(t,q) * dHdx₂(t,q,μ) + ϑ₁(t,q) * d²Hdx₂dx₃(t,q,μ) - d²Hdx₁dx₃(t,q,μ) * ϑ₂(t,q) - dHdx₁(t,q,μ) * dϑ₂dx₃(t,q)
dγ₃dx₄(t,q,μ) = dϑ₁dx₄(t,q) * dHdx₂(t,q,μ) + ϑ₁(t,q) * d²Hdx₂dx₄(t,q,μ) - d²Hdx₁dx₄(t,q,μ) * ϑ₂(t,q) - dHdx₁(t,q,μ) * dϑ₂dx₄(t,q)

function v₁(t, q, v, params)
    @unpack μ = params
    v[1] = + dβ₃dx₂(t,q,μ)
    v[2] = - dβ₃dx₁(t,q,μ)
    v[3] = 0
    v[4] = 0
end

function v₂(t, q, v, params)
    @unpack μ = params
    v[1] = - dβ₂dx₃(t,q,μ)
    v[2] = 0
    v[3] = + dβ₂dx₁(t,q,μ)
    v[4] = 0
end

function v₃(t, q, v, params)
    @unpack μ = params
    v[1] = 0
    v[2] = + dβ₁dx₃(t,q,μ)
    v[3] = - dβ₁dx₂(t,q,μ)
    v[4] = 0
end

function v₄(t, q, v, params)
    @unpack μ = params
    v[1] = + dγ₁dx₄(t,q,μ)
    v[2] = 0
    v[3] = 0
    v[4] = - dγ₁dx₁(t,q,μ)
end

function v₅(t, q, v, params)
    @unpack μ = params
    v[1] = 0
    v[2] = + dγ₂dx₄(t,q,μ)
    v[3] = 0
    v[4] = - dγ₂dx₂(t,q,μ)
end

function v₆(t, q, v, params)
    @unpack μ = params
    v[1] = 0
    v[2] = 0
    v[3] = + dγ₃dx₄(t,q,μ)
    v[4] = - dγ₃dx₃(t,q,μ)
end

function v(t, q, v, params)
    @unpack μ = params
    v[1] =   dγ₁dx₄(t,q,μ) + dβ₃dx₂(t,q,μ) - dβ₂dx₃(t,q,μ)
    v[2] =   dγ₂dx₄(t,q,μ) + dβ₁dx₃(t,q,μ) - dβ₃dx₁(t,q,μ)
    v[3] =   dγ₃dx₄(t,q,μ) + dβ₂dx₁(t,q,μ) - dβ₁dx₂(t,q,μ)
    v[4] = - dγ₁dx₁(t,q,μ) - dγ₂dx₂(t,q,μ) - dγ₃dx₃(t,q,μ)
    nothing
end


text:
################################################################################

"sort hidden Markov model by amplitude"
function sortHMM!(self::HMM)

  ampVc = map(amplitude, self.dataM)
  sortedAmpVc = sort(ampVc)
  ampIxVc = map(sortedAmpVc) do y
   findall(x -> x == y, ampVc)
  end
  templ = Array{Array{T, 1} where T <: Number, 1}()

  for i in ampIxVc
    push!(templ, self.dataM[i[1]])
  end
  self.dataM = templ
end

################################################################################

max_stars_repo_name: grabherr/HiddenMarkovModelReaders

max_issues_repo_name: grabherr/HiddenMarkovModelReaders

max_forks_repo_name: grabherr/HiddenMarkovModelReaders


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
#
# contributed by David Campbell
# modified by Jarret Revels, Kristoffer Carlsson, Alex Arslan

mutable struct PushVector{T, A<:AbstractVector{T}} <: AbstractVector{T}
    v::A
    l::Int
end

PushVector{T}() where {T} = PushVector(Vector{T}(undef, 60), 0)

Base.IndexStyle(::Type{<:PushVector}) = IndexLinear()
Base.length(v::PushVector) = v.l
Base.size(v::PushVector) = (v.l,)
@inline function Base.getindex(v::PushVector, i)
    @boundscheck checkbounds(v, i)
    @inbounds v.v[i]
end

function Base.push!(v::PushVector, i)
    v.l += 1
    if v.l > length(v.v)
        resize!(v.v, v.l * 2)
    end
    v.v[v.l] = i
    return v
end

function Base.resize!(v::PushVector, l::Integer)
    # Only support shrinking for now, since that is all we need
    @assert l <= v.l
    v.l = l
end

const _revcompdata = Dict(
   'A'=> 'T', 'a'=> 'T',
   'C'=> 'G', 'c'=> 'G',
   'G'=> 'C', 'g'=> 'C',
   'T'=> 'A', 't'=> 'A',
   'U'=> 'A', 'u'=> 'A',
   'M'=> 'K', 'm'=> 'K',
   'R'=> 'Y', 'r'=> 'Y',
   'W'=> 'W', 'w'=> 'W',
   'S'=> 'S', 's'=> 'S',
   'Y'=> 'R', 'y'=> 'R',
   'K'=> 'M', 'k'=> 'M',
   'V'=> 'B', 'v'=> 'B',
   'H'=> 'D', 'h'=> 'D',
   'D'=> 'H', 'd'=> 'H',
   'B'=> 'V', 'b'=> 'V',
   'N'=> 'N', 'n'=> 'N',
)

const revcompdata = zeros(UInt8, 256)
for (k, v) in _revcompdata
    revcompdata[k%UInt8] = v%UInt8
end

function print_buff(outio, bb)
    b = resize!(bb.v, length(bb))
    l = length(b)
    length(b) == 0 && return

    br = reverse!(b)
    for i = 1:60:l
        if i+59 > l
            write(outio, @view br[i:end])
        else
            write(outio, @view br[i:i+59])
        end
        write(outio, '\n')
    end
end

function perf_revcomp(io=stdin, outio = stdout)
    buff = PushVector{UInt8}()
    while true
        line = codeunits(readline(io))
        if isempty(line)
            print_buff(outio, buff)
            return
        elseif first(line) == UInt8('>')
            print_buff(outio, buff)
            resize!(buff, 0)
            write(outio, line)
            write(outio, '\n')
        else
            l = length(line)
            @inbounds for c in line
                push!(buff, revcompdata[c%Int])
            end
        end
    end
end

perf_revcomp()


text:
function AbsoluteJacobiansMapper(arg0::JString, arg1::ParameterDriversList)
    return AbsoluteJacobiansMapper((JString, ParameterDriversList), arg0, arg1)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_additional_state_dimension(obj::JacobiansMapper)
    return jcall(obj, "getAdditionalStateDimension", jint, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_name(obj::AbstractJacobiansMapper)
    return jcall(obj, "getName", JString, ())
end

function get_parameters(obj::AbstractJacobiansMapper)
    return jcall(obj, "getParameters", jint, ())
end

function get_parameters_jacobian(obj::JacobiansMapper, arg0::SpacecraftState, arg1::Vector{Vector{jdouble}})
    return jcall(obj, "getParametersJacobian", void, (SpacecraftState, Vector{Vector{jdouble}}), arg0, arg1)
end

function get_state_jacobian(obj::JacobiansMapper, arg0::SpacecraftState, arg1::Vector{Vector{jdouble}})
    return jcall(obj, "getStateJacobian", void, (SpacecraftState, Vector{Vector{jdouble}}), arg0, arg1)
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function set_initial_jacobians(obj::JacobiansMapper, arg0::SpacecraftState, arg1::Vector{Vector{jdouble}}, arg2::Vector{Vector{jdouble}}, arg3::Vector{jdouble})
    return jcall(obj, "setInitialJacobians", void, (SpacecraftState, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{jdouble}), arg0, arg1, arg2, arg3)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
using GitHub: GitHub

function all_specified_statuses_passed(
    api::GitHub.GitHubAPI,
    registry::GitHub.Repo,
    pr::GitHub.PullRequest,
    sha::AbstractString,
    specified_status_contexts::AbstractVector{<:AbstractString};
    auth::GitHub.Authorization,
)
    # Keep track of the values of the specified statuses. If all of
    # these are switched over to true, the result is a success. Do not
    # care about the result of unspecified statuses.
    status_passed = Dict{String,Bool}(
        context => false for context in specified_status_contexts
    )
    combined_status = GitHub.status(api, registry, sha; auth=auth)
    all_statuses = combined_status.statuses
    for status in all_statuses
        context = status.context
        if haskey(status_passed, context)
            status_passed[context] = status.state == "success"
        end
    end
    return all(values(status_passed))
end

function all_specified_check_runs_passed(
    api::GitHub.GitHubAPI,
    registry::GitHub.Repo,
    pr::GitHub.PullRequest,
    sha::AbstractString,
    specified_checks::AbstractVector{<:AbstractString};
    auth::GitHub.Authorization,
)
    # Keep track of the results of the specified checks. If all of
    # these are switched over to true, the result is a success. Do not
    # care about the result of unspecified checks.
    check_passed = Dict{String,Bool}(context => false for context in specified_checks)
    endpoint = "/repos/$(registry.full_name)/commits/$(sha)/check-runs"
    check_runs = GitHub.gh_get_json(
        api,
        endpoint;
        auth=auth,
        headers=Dict("Accept" => "application/vnd.github.antiope-preview+json"),
    )
    for check_run in check_runs["check_runs"]
        name = check_run["name"]
        check_run_was_success =
            (check_run["status"] == "completed") && (check_run["conclusion"] == "success")
        if haskey(check_passed, name)
            check_passed[name] = check_run_was_success
        end
    end
    return all(values(check_passed))
end

pr_comment_is_blocking(c::GitHub.Comment) = !occursin("[noblock]", body(c))

function pr_has_no_blocking_comments(
    api::GitHub.GitHubAPI,
    registry::GitHub.Repo,
    pr::GitHub.PullRequest;
    auth::GitHub.Authorization,
)
    all_pr_comments = get_all_pull_request_comments(api, registry, pr; auth=auth)
    return !any(pr_comment_is_blocking.(all_pr_comments))
end

function pr_is_old_enough(
    pr_type::Symbol,
    pr_age::Dates.Period;
    pkg::AbstractString,
    new_package_waiting_period::Dates.Period,
    new_jll_package_waiting_period::Dates.Period,
    new_version_waiting_period::Dates.Period,
    new_jll_version_waiting_period::Dates.Period,
    pr_author,
    authorized_authors,
    authorized_authors_special_jll_exceptions,
)
    this_is_jll_package = is_jll_name(pkg)

    if this_is_jll_package
        if pr_author in authorized_authors_special_jll_exceptions
            this_pr_can_use_special_jll_exceptions = true
        else
            this_pr_can_use_special_jll_exceptions = false
        end
    else
        this_pr_can_use_special_jll_exceptions = false
    end

    if this_pr_can_use_special_jll_exceptions
        if pr_type == :NewPackage
            return pr_age > new_jll_package_waiting_period
        elseif pr_type == :NewVersion
            return pr_age > new_jll_version_waiting_period
        else
            throw(ArgumentError("pr_type must be either :NewPackage or :NewVersion"))
        end
    else
        if pr_type == :NewPackage
            return pr_age > new_package_waiting_period
        elseif pr_type == :NewVersion
            return pr_age > new_version_waiting_period
        else
            throw(ArgumentError("pr_type must be either :NewPackage or :NewVersion"))
        end
    end
end

function _get_all_pr_statuses(
    api::GitHub.GitHubAPI,
    repo::GitHub.Repo,
    pr::GitHub.PullRequest;
    auth::GitHub.Authorization,
)
    combined_status = GitHub.status(api, repo, pr.head.sha; auth=auth)
    all_statuses = combined_status.statuses
    return all_statuses
end

function _get_status_description(status::GitHub.Status)::String
    if hasproperty(status, :description)
        if status.description === nothing
            return ""
        else
            return status.description
        end
    else
        return ""
    end
end

function _postprocess_automerge_decision_status(status::GitHub.Status; whoami)
    @debug("status: ", status)
    @debug("status.creator: ", status.creator)
    new_package_passed_regex = r"New package. Approved. name=\"(\w*)\". sha=\"(\w*)\""
    new_version_passed_regex = r"New version. Approved. name=\"(\w*)\". sha=\"(\w*)\""
    status_description = _get_status_description(status)
    if status.state == "success" && occursin(new_package_passed_regex, status_description)
        m = match(new_package_passed_regex, status_description)
        passed_pkg_name = m[1]
        passed_pr_head_sha = m[2]
        return true, passed_pkg_name, passed_pr_head_sha, :NewPackage
    end
    if status.state == "success" && occursin(new_version_passed_regex, status_description)
        m = match(new_version_passed_regex, status_description)
        passed_pkg_name = m[1]
        passed_pr_head_sha = m[2]
        return true, passed_pkg_name, passed_pr_head_sha, :NewVersion
    end
    return false, "", "", :failing
end

function pr_has_passing_automerge_decision_status(
    api::GitHub.GitHubAPI,
    repo::GitHub.Repo,
    pr::GitHub.PullRequest;
    auth::GitHub.Authorization,
    whoami,
)
    all_statuses = _get_all_pr_statuses(api, repo, pr; auth=auth)
    for status in all_statuses
        if status.context == "automerge/decision"
            return _postprocess_automerge_decision_status(status; whoami=whoami)
        end
    end
    return false, "", "", :failing
end

function cron_or_api_build(
    api::GitHub.GitHubAPI,
    registry::GitHub.Repo;
    auth::GitHub.Authorization,
    authorized_authors::Vector{String},
    authorized_authors_special_jll_exceptions::Vector{String},
    merge_new_packages::Bool,
    merge_new_versions::Bool,
    new_package_waiting_period,
    new_jll_package_waiting_period,
    new_version_waiting_period,
    new_jll_version_waiting_period,
    whoami::String,
    all_statuses::AbstractVector{<:AbstractString},
    all_check_runs::AbstractVector{<:AbstractString},
    read_only::Bool,
)
    # first, get a list of ALL open pull requests on this repository
    # then, loop through each of them.
    all_currently_open_pull_requests = my_retry(
        () -> get_all_pull_requests(api, registry, "open"; auth=auth)
    )
    reverse!(all_currently_open_pull_requests)
    at_least_one_exception_was_thrown = false
    if isempty(all_currently_open_pull_requests)
        @info("There are no open pull requests.")
    else
        for pr in all_currently_open_pull_requests
            try
                my_retry() do
                    cron_or_api_build(
                        api,
                        pr,
                        registry;
                        auth=auth,
                        authorized_authors=authorized_authors,
                        authorized_authors_special_jll_exceptions=authorized_authors_special_jll_exceptions,
                        merge_new_packages=merge_new_packages,
                        merge_new_versions=merge_new_versions,
                        new_package_waiting_period=new_package_waiting_period,
                        new_jll_package_waiting_period=new_jll_package_waiting_period,
                        new_version_waiting_period=new_version_waiting_period,
                        new_jll_version_waiting_period=new_jll_version_waiting_period,
                        whoami=whoami,
                        all_statuses=all_statuses,
                        all_check_runs=all_check_runs,
                        read_only=read_only,
                    )
                end
            catch ex
                at_least_one_exception_was_thrown = true
                showerror(stderr, ex)
                Base.show_backtrace(stderr, catch_backtrace())
                println(stderr)
            end
        end
        if at_least_one_exception_was_thrown
            throw(
                AutoMergeCronJobError(
                    "At least one exception was thrown. Check the logs for details."
                ),
            )
        end
    end
    return nothing
end

function cron_or_api_build(
    api::GitHub.GitHubAPI,
    pr::GitHub.PullRequest,
    registry::GitHub.Repo;
    auth::GitHub.Authorization,
    authorized_authors::Vector{String},
    authorized_authors_special_jll_exceptions::Vector{String},
    merge_new_packages::Bool,
    merge_new_versions::Bool,
    new_package_waiting_period,
    new_jll_package_waiting_period,
    new_version_waiting_period,
    new_jll_version_waiting_period,
    whoami::String,
    all_statuses::AbstractVector{<:AbstractString},
    all_check_runs::AbstractVector{<:AbstractString},
    read_only::Bool,
)
    #       first, see if the author is an authorized author. if not, then skip.
    #       next, see if the title matches either the "New Version" regex or
    #               the "New Package regex". if it is not either a new
    #               package or a new version, skip.
    #       next, see if it is old enough. if it is not old enough, then skip.
    #       then, get the `automerge/decision` status and make sure it is passing
    #       then, get all of the pull request comments. if there is any comment that is
    #               (1) not by me, and (2) does not contain the text [noblock], then skip
    #       if all of the above criteria were met, then merge the pull request
    pr_number = number(pr)
    @info("Now examining pull request $(pr_number)")
    pr_author = author_login(pr)
    if pr_author ∉ vcat(authorized_authors, authorized_authors_special_jll_exceptions)
        @info(
            string(
                "Pull request: $(pr_number). ",
                "Decision: do not merge. ",
                "Reason: pull request author is not authorized to automerge.",
            ),
            pr_author,
            authorized_authors,
            authorized_authors_special_jll_exceptions
        )
        return nothing
    end

    if !(is_new_package(pr) || is_new_version(pr))
        @info(
            string(
                "Pull request: $(pr_number). ",
                "Decision: do not merge. ",
                "Reason: pull request is neither a new package nor a new version.",
            ),
            title(pr)
        )
        return nothing
    end

    if is_new_package(pr) # it is a new package
        pr_type = :NewPackage
        pkg, version = parse_pull_request_title(NewPackage(), pr)
    else # it is a new version
        pr_type = :NewVersion
        pkg, version = parse_pull_request_title(NewVersion(), pr)
    end
    pr_age = time_since_pr_creation(pr)
    this_pr_is_old_enough = pr_is_old_enough(
        pr_type,
        pr_age;
        pkg=pkg,
        new_package_waiting_period=new_package_waiting_period,
        new_jll_package_waiting_period=new_jll_package_waiting_period,
        new_version_waiting_period=new_version_waiting_period,
        new_jll_version_waiting_period=new_jll_version_waiting_period,
        pr_author=pr_author,
        authorized_authors=authorized_authors,
        authorized_authors_special_jll_exceptions=authorized_authors_special_jll_exceptions,
    )
    if !this_pr_is_old_enough
        @info(
            string(
                "Pull request: $(pr_number). ",
                "Type: $(pr_type). ",
                "Decision: do not merge. ",
                "Reason: mandatory waiting period has not elapsed.",
            ),
            pr_type,
            pr_age,
            _canonicalize_period(pr_age),
            pkg,
            is_jll_name(pkg),
            new_package_waiting_period,
            _canonicalize_period(new_package_waiting_period),
            new_jll_package_waiting_period,
            _canonicalize_period(new_jll_package_waiting_period),
            new_version_waiting_period,
            _canonicalize_period(new_version_waiting_period),
            new_jll_version_waiting_period,
            _canonicalize_period(new_jll_version_waiting_period),
            pr_author,
            authorized_authors,
            authorized_authors_special_jll_exceptions
        )
        return nothing
    end

    i_passed_this_pr, passed_pkg_name, passed_pr_head_sha, status_pr_type = pr_has_passing_automerge_decision_status(
        api, registry, pr; auth=auth, whoami=whoami
    )
    if !i_passed_this_pr
        @info(
            string(
                "Pull request: $(pr_number). ",
                "Type: $(pr_type). ",
                "Decision: do not merge. ",
                "Reason: automerge/decision status is not passing",
            ),
            whoami
        )
        return nothing
    end

    always_assert(pkg == passed_pkg_name)
    always_assert(pr.head.sha == passed_pr_head_sha)
    _statuses_good = all_specified_statuses_passed(
        api, registry, pr, passed_pr_head_sha, all_statuses; auth=auth
    )
    _checkruns_good = all_specified_check_runs_passed(
        api, registry, pr, passed_pr_head_sha, all_check_runs; auth=auth
    )
    if !(_statuses_good && _checkruns_good)
        @error(
            string(
                "Pull request: $(pr_number). ",
                "Type: $(pr_type). ",
                "Decision: do not merge. ",
                "Reason: ",
                "It is not the case that ",
                "all of the specified statuses and ",
                "check runs passed. ",
            )
        )
        return nothing
    end

    if !pr_has_no_blocking_comments(api, registry, pr; auth=auth)
        @info(
            string(
                "Pull request: $(pr_number). ",
                "Type: $(pr_type). ",
                "Decision: do not merge. ",
                "Reason: pull request has one or more blocking comments.",
            )
        )
        return nothing
    end

    if pr_type == :NewPackage # it is a new package
        always_assert(status_pr_type == :NewPackage)
        if merge_new_packages
            my_comment = comment_text_merge_now()
            @info(
                string(
                    "Pull request: $(pr_number). ",
                    "Type: $(pr_type). ",
                    "Decision: merge now.",
                )
            )
            if read_only
                @info "`read_only` mode on; skipping merge"
            else
                my_retry(() -> merge!(api, registry, pr, passed_pr_head_sha; auth=auth))
            end
        else
            @info(
                string(
                    "Pull request: $(pr_number). ",
                    "Type: $(pr_type). ",
                    "Decision: do not merge. ",
                    "Reason: ",
                    "This is a new package pull request. ",
                    "All of the criteria for automerging ",
                    "were met. ",
                    "However, merge_new_packages is false, ",
                    "so I will not merge. ",
                    "If merge_new_packages had been set to ",
                    "true, I would have merged this ",
                    "pull request right now.",
                )
            )
        end
    else # it is a new version
        always_assert(pr_type == :NewVersion)
        always_assert(status_pr_type == :NewVersion)
        if merge_new_versions
            my_comment = comment_text_merge_now()
            @info(
                string(
                    "Pull request: $(pr_number). ",
                    "Type: $(pr_type). ",
                    "Decision: merge now.",
                )
            )
            if read_only
                @info "`read_only` mode on; skipping merge"
            else
                my_retry(() -> merge!(api, registry, pr, passed_pr_head_sha; auth=auth))
            end
        else
            @info(
                string(
                    "Pull request: $(pr_number). ",
                    "Type: $(pr_type). ",
                    "Decision: do not merge. ",
                    "Reason: merge_new_versions is false",
                    "This is a new version pull request. ",
                    "All of the criteria for automerging ",
                    "were met. ",
                    "However, merge_new_versions is false, ",
                    "so I will not merge. ",
                    "If merge_new_versions had been set to ",
                    "true, I would have merged this ",
                    "pull request right now.",
                )
            )
        end
    end
    return nothing
end


text:
export StraightLineHomotopy

"""
    StraightLineHomotopy(G::System, F::System; gamma = 1.0)
    StraightLineHomotopy(G::AbstractSystem, F::AbstractSystem; gamma = 1.0)

Constructs the straight line homotopy ``H(x, t) = γ t G(x) + (1-t) F(x)`` where
``γ`` is `gamma`.
"""
struct StraightLineHomotopy{S<:AbstractSystem,T<:AbstractSystem} <: AbstractHomotopy
    start::S
    target::T
    γ::ComplexF64

    u::Vector{ComplexF64}
    ū::Vector{ComplexDF64}
    U::Matrix{ComplexF64}

    dv_start::LA.Transpose{ComplexF64,Matrix{ComplexF64}}
    dv_target::LA.Transpose{ComplexF64,Matrix{ComplexF64}}
end

function StraightLineHomotopy(start::ModelKit.System, target::ModelKit.System; kwargs...)
    StraightLineHomotopy(ModelKitSystem(start), ModelKitSystem(target); kwargs...)
end
function StraightLineHomotopy(
    start::AbstractSystem,
    target::AbstractSystem;
    γ = 1.0,
    gamma = γ,
)
    size(start) == size(target) ||
        throw(ArgumentError("Start and target do not have the same size, got $(size(start)) and $(size(target))"))

    m, n = size(start)
    u = zeros(ComplexF64, m)
    ū = zeros(ComplexDF64, m)
    U = zeros(ComplexF64, m, n)

    dv_start = LA.transpose(zeros(ComplexF64, 5, m))
    dv_target = LA.transpose(zeros(ComplexF64, 5, m))

    StraightLineHomotopy(start, target, ComplexF64(gamma), u, ū, U, dv_start, dv_target)
end

Base.size(H::StraightLineHomotopy) = size(H.start)

function Base.show(io::IO, mime::MIME"text/plain", H::StraightLineHomotopy)
    println(io, typeof(H), ":")
    println(io, "γ: ", H.γ)
    println(io, "\nG: ")
    show(io, mime, H.start)
    println(io, "\n\nF:")
    show(io, mime, H.target)
end

function evaluate!(u, H::StraightLineHomotopy, x::AbstractVector{T}, t) where {T}
    evaluate!(u, H.start, x)

    if T isa ComplexDF64 || T isa DoubleF64
        evaluate!(H.ū, H.target, x)
        ts, tt = H.γ .* t, 1 - t
        for i = 1:size(H, 1)
            @inbounds u[i] = ts * u[i] + tt * H.ū[i]
        end
    else
        evaluate!(H.u, H.target, x)
        ts, tt = H.γ .* t, 1 - t
        for i = 1:size(H, 1)
            @inbounds u[i] = ts * u[i] + tt * H.u[i]
        end
    end
    u
end

function evaluate_and_jacobian!(
    u,
    U,
    H::StraightLineHomotopy,
    x::AbstractVector{T},
    t,
) where {T}
    evaluate_and_jacobian!(u, U, H.start, x)
    evaluate_and_jacobian!(H.u, H.U, H.target, x)
    ts, tt = H.γ .* t, 1 - t
    for i = 1:size(H, 1)
        @inbounds u[i] = ts * u[i] + tt * H.u[i]
    end
    for j = 1:size(H, 2), i = 1:size(H, 1)
        @inbounds U[i, j] = ts * U[i, j] + tt * H.U[i, j]
    end
    nothing
end
function taylor!(u, ::Val{1}, H::StraightLineHomotopy, x, t)
    evaluate!(u, H.start, x)
    evaluate!(H.u, H.target, x)
    for i = 1:size(H, 1)
        @inbounds u[i] = H.γ * u[i] - H.u[i]
    end
    u
end
function taylor!(u, v::Val{K}, H::StraightLineHomotopy, tx::TaylorVector, t) where {K}
    taylor!(H.dv_start, v, H.start, tx)
    taylor!(H.dv_target, v, H.target, tx)

    for i = 1:size(H, 1)
        start = H.γ * (H.dv_start[i, K] + t * H.dv_start[i, K+1])
        target = (1 - t) * H.dv_target[i, K+1] - H.dv_target[i, K]
        u[i] = start + target
    end
    u
end


text:
struct ToSTDOUT end

struct IOConfig
    stdin::Union{String,Nothing}
    stdout::Union{String,Nothing}
    stderr::Union{String,Nothing,ToSTDOUT}
end

function ioconfig(input::AbstractDict{<:AbstractString})
    function pathfor(name)
        stream = get(input, name, nothing)
        stream === nothing && return nothing
        path = get(stream, "path", nothing)
        path == nothing || path isa AbstractString ||
            return ArgumentError("Keyword argument `$name.path` must be a `nothing` or `String`.")
        return path
    end

    function pathforstderr()
        path = pathfor("stderr")
        path === nothing || return path
        stream = get(input, "stderr", nothing)
        stream === nothing && return ToSTDOUT()
        to = get(stream, "to", "file")
        to === "file" && return path
        to in ("file", "stdout") ||
            return ArgumentError("""Keyword argument `stderr.to` must be `"file"` or `"stdout"`.""")
        @assert to === "stdout"
        return ToSTDOUT()
    end

    return IOConfig(
        (@return_exception pathfor("stdin")),
        (@return_exception pathfor("stdout")),
        (@return_exception pathforstderr()),
    )
end

function with_pipe(f)
    pipe = Base.link_pipe!(Pipe())
    reader = @async read(pipe, String)
    result = try
        f(pipe.in)
    finally
        close(pipe)
    end
    return (result, fetch(reader))
end

function with_file(f, path; kwargs...)
    result = open(path; kwargs...) do io
        fileno = fd(io)
        if ccall(:isatty, Cint, (Cint,), fileno) == 1
            tty = Base.TTY(Libc.dup(RawFD(fileno)))
            have_color = Base.have_color
            try
                @eval Base have_color = true
                f(tty)
            finally
                close(tty)
                @eval Base have_color = $have_color
            end
        else
            f(io)
        end
    end
    return (result, nothing)
end

with_stdout(f, ::Nothing) = with_pipe(io -> redirect_stdout(f, io))
with_stderr(f, ::Nothing) = with_pipe(io -> redirect_stderr(f, io))

with_stdout(f, path::AbstractString) =
    with_file(io -> redirect_stdout(f, io), path, write = true)
with_stderr(f, path::AbstractString) =
    with_file(io -> redirect_stderr(f, io), path, write = true)

with_stderr(f, ::ToSTDOUT) = redirect_stderr(f, stdout), nothing

with_stdin(f, ::Nothing) = f(), nothing  # use devnull?
with_stdin(f, path::AbstractString) =
    with_file(io -> redirect_stdin(f, io), path, read = true)

function with_stdio(f, ioconfig::IOConfig)
    (((result, _), err), out) = with_stdout(ioconfig.stdout) do
        with_stderr(ioconfig.stderr) do
            with_stdin(ioconfig.stdin) do
                f()
            end
        end
    end
    return (result, out, err)
end

function with_newlogger(f)
    if stderr isa Base.TTY
        logger = Logging.ConsoleLogger(stderr)
    else
        logger = Base.CoreLogging.SimpleLogger(stderr)
    end
    return Base.CoreLogging.with_logger(f, logger)
end


text:
# 1D linear convection problem

# Packages and inclusions
using ModelingToolkit,DiffEqOperators,DiffEqBase,LinearAlgebra,Test

# Tests

@testset "Test 00: Dt(u(t,x)) ~ -Dx(u(t,x))" begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ -Dx(u(t,x))
    bcs = [u(0,x) ~ (0.5/(0.2*sqrt(2.0*3.1415)))*exp(-(x  -0.75)^2/(2.0*0.2^2)),
           u(t,0) ~ 0.0,
           u(t,2) ~ 0.0]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,0.6),
               x ∈ IntervalDomain(0.0,2.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])

    # Method of lines discretization
    dx = 2/80
    order = 1
    discretization = MOLFiniteDifference(dx,order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    using OrdinaryDiffEq
    sol = solve(prob,Euler(),dt=.025,saveat=0.1)

    # Plot and save results
    # using Plots
    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,5]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,6]))
    # savefig("MOL_1D_Linear_Convection_Test00.png")

    # Test
    x_interval = domains[2].domain.lower+dx:dx:domains[2].domain.upper-dx
    u = @. (0.5/(0.2*sqrt(2.0*3.1415)))*exp(-(x_interval-(0.75+0.6))^2/(2.0*0.2^2))
    t_f = size(sol,3)

    @test sol[:,1,t_f] ≈ u atol = 0.1;

end

@testset "Test 01: Dt(u(t,x)) ~ -Dx(u(t,x)) + 0.01" begin

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ -Dx(u(t,x)) + 0.01
    bcs = [u(0,x) ~ (0.5/(0.2*sqrt(2.0*3.1415)))*exp(-(x  -0.75)^2/(2.0*0.2^2)),
           u(t,0) ~ 0.0,
           u(t,2) ~ 0.0]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,0.6),
               x ∈ IntervalDomain(0.0,2.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])

    # Method of lines discretization
    dx = 2/80
    order = 1
    discretization = MOLFiniteDifference(dx,order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    using OrdinaryDiffEq
    sol = solve(prob,Euler(),dt=.025,saveat=0.1)

    # Plot and save results
    # using Plots
    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,5]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,6]))
    # savefig("MOL_1D_Linear_Convection_Test01.png")

    # Test
    x_interval = domains[2].domain.lower+dx:dx:domains[2].domain.upper-dx
    u = @. (0.5/(0.2*sqrt(2.0*3.1415)))*exp(-(x_interval-(0.75+0.6))^2/(2.0*0.2^2))
    t_f = size(sol,3)
    @test sol[:,1,t_f] ≈ u atol = 0.1;

end

@testset "Test 02: Dt(u(t,x)) ~ -v*Dx(u(t,x))" begin
    # Parameters, variables, and derivatives
    @parameters t x v
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    v = 1.0

    # 1D PDE and boundary conditions
    eq  = Dt(u(t,x)) ~ -v*Dx(u(t,x))
    bcs = [u(0,x) ~ (0.5/(0.2*sqrt(2.0*3.1415)))*exp(-(x-0.75)^2/(2.0*0.2^2)),
           u(t,0) ~ 0.0,
           u(t,2) ~ 0.0]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,0.6),
               x ∈ IntervalDomain(0.0,2.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u])

    # Method of lines discretization
    dx = 2/80
    order = 1
    discretization = MOLFiniteDifference(dx,order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    using OrdinaryDiffEq
    sol = solve(prob,Euler(),dt=.025,saveat=0.1)

    # Plot and save results
    # using Plots
    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,5]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,6]))
    # savefig("MOL_1D_Linear_Convection_Test02.png")

    # Test
    x_interval = domains[2].domain.lower+dx:dx:domains[2].domain.upper-dx
    u = @. (0.5/(0.2*sqrt(2.0*3.1415)))*exp(-(x_interval-(0.75+v*0.6))^2/(2.0*0.2^2))
    t_f = size(sol,3)
    @test sol[:,1,t_f] ≈ u atol = 0.1;
end

@testset "Test 03: Dt(u(t,x)) ~ -Dx(v(t,x))*u(t,x)-v(t,x)*Dx(u(t,x)) with v(t,x)=1" begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables v(..) u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    # 1D PDE and boundary conditions
    eq  = [ Dt(u(t,x)) ~ -(Dx(v(t,x))*u(t,x)+v(t,x)*Dx(u(t,x))),
            v(t,x) ~ 1.0 ]
    bcs = [u(0,x) ~ (0.5/(0.2*sqrt(2.0*3.1415)))*exp(-(x-0.75)^2/(2.0*0.2^2)),
           u(t,0) ~ 0.0,
           u(t,2) ~ 0.0,
           v(0,x) ~ 1.0,
           v(t,0) ~ 1.0,
           v(t,2) ~ 1.0
           ]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,0.6),
               x ∈ IntervalDomain(0.0,2.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u,v])

    # Method of lines discretization
    dx = 2/80
    order = 1
    discretization = MOLFiniteDifference(dx,order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    using OrdinaryDiffEq
    sol = solve(prob,Euler(),dt=.025,saveat=0.1)

    #Plot and save results
    # using Plots
    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,5]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,6]))
    # savefig("MOL_1D_Linear_Convection_Test03.png")

    # Test
    x_interval = domains[2].domain.lower+dx:dx:domains[2].domain.upper-dx
    u = @. (0.5/(0.2*sqrt(2.0*3.1415)))*exp(-(x_interval-(0.75+1.0*0.6))^2/(2.0*0.2^2))
    t_f = size(sol,3)

    @test sol[:,1,t_f] ≈ u atol = 0.1;
end

@testset "Test 04: Dt(u(t,x)) ~ -Dx(v(t,x))*u(t,x)-v(t,x)*Dx(u(t,x)) with v(t,x)=0.999 + 0.001 * t * x " begin
    # Parameters, variables, and derivatives
    @parameters t x
    @variables v(..) u(..)
    Dt = Differential(t)
    Dx = Differential(x)

    # 1D PDE and boundary conditions
    eq  = [ Dt(u(t,x)) ~ -Dx(v(t,x))*u(t,x)-v(t,x)*Dx(u(t,x)),
            v(t,x) ~ 0.999 + 0.001 * t * x ]
    bcs = [u(0,x) ~ (0.5/(0.2*sqrt(2.0*3.1415)))*exp(-(x-0.75)^2/(2.0*0.2^2)),
           u(t,0) ~ 0.0,
           u(t,2) ~ 0.0,
           v(0,x) ~ 0.999,
           v(t,0) ~ 0.999,
           v(t,2) ~ 0.999 + 0.001 * t * 2.0]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,0.6),
               x ∈ IntervalDomain(0.0,2.0)]

    # PDE system
    pdesys = PDESystem(eq,bcs,domains,[t,x],[u,v])

    # Method of lines discretization
    dx = 2/80
    order = 1
    discretization = MOLFiniteDifference(dx,order)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    using OrdinaryDiffEq
    sol = solve(prob,Euler(),dt=.025,saveat=0.1)

    # Plot and save results
    # using Plots
    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,5]))
    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,6]))
    # savefig("MOL_1D_Linear_Convection_Test04.png")

    # Test
    x_interval = domains[2].domain.lower+dx:dx:domains[2].domain.upper-dx
    u = @. (0.5/(0.2*sqrt(2.0*3.1415)))*exp(-(x_interval-(0.75+1.0*0.6))^2/(2.0*0.2^2))
    t_f = size(sol,3)

    @test sol[:,1,t_f] ≈ u atol = 0.1;
end


text:
module SpringCollab2020CustomSandwichLava

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/CustomSandwichLava" CustomSandwichLava(x::Integer, y::Integer,
    direction::String="CoreModeBased", speed::Number=20.0, sandwichGap::Number=160.0)

const directions = String["AlwaysUp", "AlwaysDown", "CoreModeBased"]

const placements = Ahorn.PlacementDict(
    "Custom Sandwich Lava (Spring Collab 2020)" => Ahorn.EntityPlacement(
        CustomSandwichLava
    )
)

Ahorn.editingOptions(entity::CustomSandwichLava) = Dict{String, Any}(
    "direction" => directions
)

function Ahorn.selection(entity::CustomSandwichLava)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CustomSandwichLava, room::Maple.Room)
    direction = get(entity.data, "direction", "CoreModeBased")

    if direction == "AlwaysUp"
        Ahorn.drawSprite(ctx, "ahorn/SpringCollab2020/lava_sandwich_up", 0, 0)
    elseif direction == "AlwaysDown"
        Ahorn.drawSprite(ctx, "ahorn/SpringCollab2020/lava_sandwich_down", 0, 0)
    else
        Ahorn.drawImage(ctx, Ahorn.Assets.lavaSanwitch, -12, -12)
    end
end

end


text:
@testset "733.flood-fill.jl" begin
    @test flood_fill!([[1, 1, 1] [1, 1, 0] [1, 0, 1]], 1, 1, 1) ==
        [[1, 1, 1] [1, 1, 0] [1, 0, 1]]
    @test flood_fill!([[1, 1, 1] [1, 1, 0] [1, 0, 1]], 1, 1, 2) ==
        [[2, 2, 2] [2, 2, 0] [2, 0, 1]]
    @test flood_fill!([[0, 0, 0] [0, 1, 0]], 1, 1, 2) == [[0, 0, 0] [0, 2, 0]]
end


text:
include("imports.jl")

@testset ExtendedTestSet "XTree constructors" begin

    @test_nowarn XNode(:tag)
    @test_throws MethodError XNode(:tag, [])
    @test_nowarn XNode(:tag, XLeaf.(1:10))

end

@testset ExtendedTestSet "fold" begin
    x = XNode(:tag, XLeaf.(1:10))
    @test foldleaves(+, x, 0) == sum(1:10)
end

@testset ExtendedTestSet "cata" begin
    x = XNode(:tag, XLeaf.(1:10))
    x_ = cata(x) do node
        if node isa XLeaf
            return XLeaf(-node[])
        else
            return node
        end
    end
    @test x_ == XNode(:tag, XLeaf.(-1:-1:-10))
end

@testset ExtendedTestSet "catafold" begin
    x = XNode(:tag, XLeaf.(1:10))
    x_, n = catafold(x, 0) do node, state
        if node isa XLeaf
            return XLeaf(-node[]), state + 1
        else
            return node, state
        end
    end
    @test x_ == XNode(:tag, XLeaf.(-1:-1:-10))
    @test n == 10
end


text:
push!(LOAD_PATH, "Modules")
using Colors
using JuliRay

φ=MathConstants.φ

Base.@irrational ° 0.0174532925199432957692369076848861271344 (big(pi)/big(180))

function POLYGON(vertices;r1=1.0,r2=0.05,r3=0.025)
    n=length(vertices)
    V=rgbColor(csgUnion((p->Sphere(p,r2)).(vertices)),RGB(0.1,0.1,0.1))
    E=rgbColor(csgUnion([Cylinder(vertices[i],vertices[mod(i,n)+1],r3) for i in 1:n]),RGB(0.1,0.1,0.1))
    F=rgbftColor(Polygon(vertices),RGB(0.5,1,1),FT(0.1,0.1))
    return csgUnion(V,E,F)
end

function Mirror(q,p1,p2,p3)
    e₃=NormalVector(p1,p2,p3)
    e₁=OrthogonalVector(e₃)
    e₂=cross(e₃,e₁)
    A=hcat(e₁,e₂,e₃)*transpose(hcat(e₁,e₂,-e₃))
    fixedpoint=(p1+p2+p3)/3
    b=fixedpoint-A*fixedpoint
    return A*q+b
end

function Put(F,i,O)
    p1=O
    n=length(F)
    p2,p3=F[mod(i-1,n)+1],F[mod(i,n)+1]
    return (p->Mirror(p,p1,p2,p3)).(F)
end

function f(x)
    ifelse(x>0,exp(-1/x),0)
end
function g(x)
    f(2x/√3)/(f(2x/√3)+f(2/√3-2x/√3))
end
function Smooth(a,b,x)
    g((x-a)/(b-a))
end

## Tetrahedron

h=1/3 # 単位球面に内接させた際の内接球半径
n=3 # n角形
m=4 # m面体
M=60 # アニメーションの刻み数

r1=√(1-h^2) # n角形の外接円半径
for i ∈ 1:M
    θ=(1-Smooth(0,1,2(i-1)/M)+Smooth(1,2,2(i-1)/M))*π/2 # パラメータθは0→1→0の順で滑らかに変化
    R=1/cos(θ)
    O=[0,0,tan(θ)]
    N=O+[0,0,R]
    H=O-[0,0,√(R^2-r1^2)]

    F=[[[r1*cos(2π*i/n),r1*sin(2π*i/n),0.0]+H for i in 1:n] for i ∈ 1:m] # 最初のn角形の頂点
    F[2]=Put(F[1],1,O) # 以下はn角形を付け加える操作
    F[3]=Put(F[2],2,O)
    F[4]=Put(F[2],3,O)

    object=csgUnion([POLYGON(F[i],r2=0.03,r3=0.01) for i ∈ 1:m]) # 頂点座標をもとにオブジェクトを生成
    render(object,camera=LngLatCamera(lng=30°,lat=30°,pers=0.2,zoom=0.2,width=640,height=360),name="Np"*string(m),index=i)
end

max_stars_repo_name: hyrodium/JuliRay

max_issues_repo_name: hyrodium/JuliRay

max_forks_repo_name: hyrodium/JuliRay


text:
# lidHttpServer.jl
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

using Text
using Base.Test, Stage, Ollam, DataStructures
using HttpServer

filename = "./data/text.tsv"
if size(ARGS, 1) > 0
    filename = ARGS[1]
end

# ----------------------------------------------------------------------
# LID
# ----------------------------------------------------------------------

# train the model

train       = map(l -> split(chomp(l), '\t')[2], filelines(filename))
train_truth = map(l -> split(chomp(l), '\t')[1], filelines(filename))

bkgmodel, fextractor, model = tc_train(train, train_truth, lid_iterating_tokenizer, mincount = 2, cutoff = 1e10,
                                       trainer = (fvs, truth, init_model) -> train_mira(fvs, truth, init_model, iterations = 20, k = 2, C = 0.01, average = true),
                                       iteration_method = :eager)

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

# function to identify language from text using the model

function detectLang(text)
    fv      = fextractor(lid_iterating_tokenizer(text))
    scores  = score(model, fv)
    bidx, b = best(scores)
    return model.index_class[bidx]
end

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

# HTTP Server to identify language from PUT Request data

http = HttpHandler() do req::Request, res::Response
    if req.method != "PUT"
        return Response("{\"lang\": \"error - use PUT method\"}")
    end
    text = bytestring(req.data)
    # println(text)
    # println(detectLang(text))
    return Response(string("{\"lang\": \"", detectLang(text), "\"}"))
end
http.events["error"]  = (client, err) -> println(err)
http.events["listen"] = (saddr) -> println("Running on http://$saddr (Press CTRL+C to quit)")

server = Server(http)
run(server, host=IPv4(127,0,0,1), port=8000)


text:
module Search

import Moves
import Util

global node_count = 0

function staticEvaluation(board,mobility,colour)
    opp_colour = (colour == 'W') ? 'B' : 'W'
    score = count_ones(board[colour * 'P']) - count_ones(board[opp_colour * 'P'])
    score = score + (3 * (count_ones(board[colour * 'N']) + count_ones(board[colour * 'B']) - (count_ones(board[opp_colour * 'N']) + count_ones(board[opp_colour * 'B']))))
    score = score + (5 * (count_ones(board[colour * 'R']) - count_ones(board[opp_colour * 'R'])))
    score = score + (9 * (count_ones(board[colour * 'Q']) - count_ones(board[opp_colour * 'Q'])))
    score = score + (200 * (count_ones(board[colour * 'K']) - count_ones(board[opp_colour * 'K'])))
    score = score + (0.1 * mobility)
    under_attack = Moves.getAttackBoard(board,colour,false) & Util.getWhitePieces(board)
    opp_under_attack = Moves.getAttackBoard(board,opp_colour,false) & Util.getBlackPieces(board)
    #number of pieces in attack
    score = score + (0.2 * (count_ones(opp_under_attack) - count_ones(under_attack)))
    #find if pieces are protected that are under attack
   # score = score + (0.2 * (count_ones(under_attack & getAttackBoard(board,colour,false)) - count_ones(opp_under_attack & getAttackBoard(board,opp_colour,false))))
    return score
end

function search(depth,colour,history,board)
    global node_count
    w = @time negamax(-Inf,Inf,depth,colour,history,board,0)
    println("weight of move: ", w, ' ', best_move)
    println("nodes explored: ", node_count)
    node_count = 0
end

function negamax(alpha,beta,depth,colour,history,board, ply)
    opp_colour = (colour == 'W') ? 'B' : 'W'
    global node_count = node_count + 1
    #  println(depth,' ',alpha,' ',beta, ' ', moves)
    moves = Moves.getMoves(board,colour,history)
    if(depth == 0)
        #  println("staticEval: ",alpha, " ",beta, " ")
        x = staticEvaluation(board, length(moves),colour)
        # println("return staticEvaluation: ", x)
        return x
    end
 #   node_count = node_count + 1
    old_alpha = alpha
    bestmove_sofar = ""
    for move in moves
        copied = Moves.makeMoveUCI(move,board,colour)
        ply = ply + 1
        score = -negamax(-beta,-alpha,depth-1,opp_colour,history,copied,ply)
        ply = ply - 1
        # println("SCORE after -negamax: ", score, ' ',depth, " ",alpha, " ",beta, " ", move, " ", ply)

        if(score >= beta)
            #  println("returning beta: ", beta )
            return beta
        end
        if(score > alpha) #if better move
            alpha = score
            #   println("updating alpha, " , alpha)
            if(ply == 0)
                # println("updating best move : ", move)
                bestmove_sofar = move
            end
        end

    end
    if(old_alpha != alpha)
        global best_move = bestmove_sofar
    end
    if(length(moves) == 0)
        mate = Moves.isCheckMate(board,colour,history)
        if(mate === 1)
            return -100000
        elseif(mate === 2)
            return 0
        end
    end
    return alpha
end

end

max_stars_repo_name: grahamclyne/chessEngine

max_issues_repo_name: grahamclyne/chessEngine

max_forks_repo_name: grahamclyne/chessEngine


text:
module GadgetUnits

    using Unitful
    using UnitfulAstro
    using Cosmology
    using GadgetIO

    include("conversion_structs/with_units.jl")
    include("conversion_structs/without_units.jl")
    include("cosmo_utility/cosmo.jl")
    include("cosmo_utility/time.jl")
    include("cosmo_utility/radiation.jl")
    include("cosmo_utility/scale.jl")

    export GadgetPhysical,
           GadgetPhysicalUnits,
           cosmology,
           age, redshift,
           arcmin_to_kpc,
           mJy_to_W

end # module


text:
struct DefaultInit <: DiffEqBase.DAEInitializationAlgorithm end
struct NoInit <: DiffEqBase.DAEInitializationAlgorithm end

struct ShampineCollocationInit{T} <: DiffEqBase.DAEInitializationAlgorithm
	initdt::T
end
ShampineCollocationInit() = ShampineCollocationInit(nothing)

struct BrownFullBasicInit{T} <: DiffEqBase.DAEInitializationAlgorithm
	abstol::T
end
BrownFullBasicInit() = BrownFullBasicInit(1e-10)

## Notes

#=
differential_vars = [any(!iszero,x) for x in eachcol(M)]

A column should be zero for an algebraic variable, since that means that the
derivative term doesn't show up in any equations (i.e. is an algebraic variable).
The rows are not necessarily non-zero, for example a flux condition between two
differential variables. But if it's a condition that doesn't involve the algebraic
variable, then the system is not Index 1!

=#

## Expansion

function DiffEqBase.initialize_dae!(integrator::ODEIntegrator,
									initializealg = integrator.initializealg)
	_initialize_dae!(integrator, integrator.sol.prob,
					 initializealg,
					 Val(DiffEqBase.isinplace(integrator.sol.prob)))
end

## Default algorithms

function _initialize_dae!(integrator, prob::ODEProblem,
						 alg::DefaultInit, x::Val{true})
	_initialize_dae!(integrator, prob,
					BrownFullBasicInit(), x)
end

function _initialize_dae!(integrator, prob::ODEProblem,
						 alg::DefaultInit, x::Val{false})
	_initialize_dae!(integrator, prob,
					BrownFullBasicInit(), x)
end

function _initialize_dae!(integrator, prob::DAEProblem,
						 alg::DefaultInit, x::Val{false})
	if prob.differential_vars === nothing
		_initialize_dae!(integrator, prob,
						ShampineCollocationInit(), x)
	else
		_initialize_dae!(integrator, prob,
						BrownFullBasicInit(), x)
	end
end

function _initialize_dae!(integrator, prob::DAEProblem,
						 alg::DefaultInit, x::Val{true})
	if prob.differential_vars === nothing
		_initialize_dae!(integrator, prob,
						ShampineCollocationInit(), x)
	else
		_initialize_dae!(integrator, prob,
						BrownFullBasicInit(), x)
	end
end

## NoInit

function _initialize_dae!(integrator, prob::ODEProblem,
						 alg::NoInit, x::Val{true})
end

function _initialize_dae!(integrator, prob::ODEProblem,
						 alg::NoInit, x::Val{false})
end

function _initialize_dae!(integrator, prob::DAEProblem,
						 alg::NoInit, x::Val{false})
end

function _initialize_dae!(integrator, prob::DAEProblem,
						 alg::NoInit, x::Val{true})
end

## ShampineCollocationInit

#=
The method:

du = (u-u0)/h
Solve for `u`

=#

function _initialize_dae!(integrator, prob::ODEProblem,
						 alg::ShampineCollocationInit, ::Val{true})

	@unpack p, t, f = integrator
 	M = integrator.f.mass_matrix
	dtmax = integrator.opts.dtmax
	tmp = first(get_tmp_cache(integrator))
	u0 = integrator.u

	dt = t != 0 ? min(t/1000,dtmax) : dtmax # Haven't implemented norm reduction

 	nlequation! = function (out,u)
		update_coefficients!(M,u,p,t)
 		#M * (u-u0)/dt - f(u,p,t)
		@. tmp = (u - u0)/dt
		mul!(out,M,tmp)
		f(tmp,u,p,t)
		out .-= tmp
		nothing
 	end

	update_coefficients!(M,u0,p,t)
	algebraic_vars = [all(iszero,x) for x in eachcol(M)]
	algebraic_eqs  = [all(iszero,x) for x in eachrow(M)]
  (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return
	f(tmp,u0,p,t)
	tmp .= algebraic_eqs .* tmp

 	integrator.opts.internalnorm(tmp,t) <= integrator.opts.abstol && return

	integrator.u .= nlsolve(nlequation!, u0).zero
	recursivecopy!(integrator.uprev,integrator.u)
	if alg_extrapolates(integrator.alg)
		recursivecopy!(integrator.uprev2,integrator.uprev)
	end
  return nothing
end

function _initialize_dae!(integrator, prob::ODEProblem,
						 alg::ShampineCollocationInit, ::Val{false})

	@unpack p, t, f = integrator
	u0 = integrator.u
	M = integrator.f.mass_matrix
	dtmax = integrator.opts.dtmax

	dt = t != 0 ? min(t/1000,dtmax/10) : dtmax # Haven't implemented norm reduction

	update_coefficients!(M,u0,p,t)
	algebraic_vars = [all(iszero,x) for x in eachcol(M)]
	algebraic_eqs  = [all(iszero,x) for x in eachrow(M)]
  (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return
	du = f(u0,p,t)
	resid = du[algebraic_eqs]

  integrator.opts.internalnorm(resid,t) <= integrator.opts.abstol && return

	nlequation_oop = function (u)
		update_coefficients!(M,u,p,t)
		M * (u-u0)/dt - f(u,p,t)
	end

	nlequation! = (out,u) -> out .= nlequation_oop(u)

	integrator.u = nlsolve(nlequation!, u0).zero
	integrator.uprev = integrator.u
	if alg_extrapolates(integrator.alg)
		integrator.uprev2 = integrator.uprev
	end
  return
end

function _initialize_dae!(integrator, prob::DAEProblem,
						 alg::ShampineCollocationInit, ::Val{true})

	@unpack p, t, f = integrator
	u0 = integrator.u

	dtmax = integrator.opts.dtmax
	tmp = get_tmp_cache(integrator)[1]
	resid = get_tmp_cache(integrator)[2]

	dt = t != 0 ? min(t/1000,dtmax) : dtmax # Haven't implemented norm reduction

 	nlequation! = function (out,u)
 		#M * (u-u0)/dt - f(u,p,t)
		@. tmp = (u - u0)/dt
		f(out,tmp,u,p,t)
		nothing
 	end

	nlequation!(tmp,u0)
	f(resid,integrator.du,u0,p,t)
 	integrator.opts.internalnorm(resid,t) <= integrator.opts.abstol && return

	integrator.u .= nlsolve(nlequation!, u0).zero
	recursivecopy!(integrator.uprev,integrator.u)
	if alg_extrapolates(integrator.alg)
		recursivecopy!(integrator.uprev2,integrator.uprev)
	end
  return
end

function _initialize_dae!(integrator, prob::DAEProblem,
						 alg::ShampineCollocationInit, ::Val{false})

	@unpack p, t, f = integrator
	u0 = integrator.u
	dtmax = integrator.opts.dtmax

	dt = t != 0 ? min(t/1000,dtmax/10) : dtmax # Haven't implemented norm reduction

	nlequation_oop = function (u)
		f((u-u0)/dt,u,p,t)
	end

	nlequation! = (out,u) -> out .= nlequation_oop(u)

	resid = f(integrator.du,u0,p,t)
	integrator.opts.internalnorm(resid,t) <= integrator.opts.abstol && return

	integrator.u = nlsolve(nlequation!, u0).zero
	integrator.uprev = integrator.u
	if alg_extrapolates(integrator.alg)
		integrator.uprev2 = integrator.uprev
	end
end

## BrownFullBasic

#=
The method:

Keep differential variables constant
Solve for the algebraic variables

=#

function _initialize_dae!(integrator, prob::ODEProblem,
						 alg::BrownFullBasicInit, ::Val{true})
	@unpack p, t, f = integrator
	u = integrator.u
	M = integrator.f.mass_matrix
	update_coefficients!(M,u,p,t)
	algebraic_vars = [all(iszero,x) for x in eachcol(M)]
	algebraic_eqs  = [all(iszero,x) for x in eachrow(M)]
  (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return
	tmp = get_tmp_cache(integrator)[1]

	f(tmp,u,p,t)

	tmp .= algebraic_eqs .* tmp

	integrator.opts.internalnorm(tmp,t) <= alg.abstol && return
	alg_u = @view u[algebraic_vars]

	nlequation = (out, x) -> begin
		alg_u .= x
		f(tmp, u, p, t)
		out .= @view tmp[algebraic_eqs]
	end

	r = nlsolve(nlequation, u[algebraic_vars])
	alg_u .= r.zero

	recursivecopy!(integrator.uprev,integrator.u)
	if alg_extrapolates(integrator.alg)
		recursivecopy!(integrator.uprev2,integrator.uprev)
	end

	return
end

function _initialize_dae!(integrator, prob::ODEProblem,
						 alg::BrownFullBasicInit, ::Val{false})
	@unpack p, t, f = integrator

	u0 = integrator.u
	M = integrator.f.mass_matrix
	update_coefficients!(M,u0,p,t)
	algebraic_vars = [all(iszero,x) for x in eachcol(M)]
	algebraic_eqs  = [all(iszero,x) for x in eachrow(M)]
  (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return

	du = f(u0,p,t)
	resid = du[algebraic_eqs]

  integrator.opts.internalnorm(resid,t) <= alg.abstol && return

	if u0 isa Number
		# This doesn't fix static arrays!
		u = [u0]
	else
		u = u0
	end

	alg_u = @view u[algebraic_vars]

	nlequation = (out,x) -> begin
		alg_u .= x
		du = f(u,p,t)
		out .= @view du[algebraic_eqs]
	end
	r = nlsolve(nlequation, u0[algebraic_vars])
	alg_u .= r.zero

	if u0 isa Number
		# This doesn't fix static arrays!
		integrator.u = first(u)
	else
		integrator.u = u
	end

	integrator.uprev = integrator.u
	if alg_extrapolates(integrator.alg)
		integrator.uprev2 = integrator.uprev
	end

	return
end

function _initialize_dae!(integrator, prob::DAEProblem,
						 alg::BrownFullBasicInit, ::Val{true})
	@unpack p, t, f = integrator
	differential_vars = prob.differential_vars
	u = integrator.u
	du = integrator.du

	tmp = get_tmp_cache(integrator)[1]
	f(tmp, du, u, p, t)

	if integrator.opts.internalnorm(tmp,t) <= alg.abstol
		return
	elseif differential_vars === nothing
		error("differential_vars must be set for DAE initialization to occur. Either set consistent initial conditions, differential_vars, or use a different initialization algorithm.")
	end

	nlequation = (out, x) -> begin
		@. du = ifelse(differential_vars,x,du)
		@. u  = ifelse(differential_vars,u,x)
		f(out, du, u, p, t)
	end

	r = nlsolve(nlequation, ifelse.(differential_vars,du,u))

	@. du = ifelse(differential_vars,r.zero,du)
	@. u  = ifelse(differential_vars,u,r.zero)

	recursivecopy!(integrator.uprev,integrator.u)
	if alg_extrapolates(integrator.alg)
		recursivecopy!(integrator.uprev2,integrator.uprev)
	end

	return
end

function _initialize_dae!(integrator, prob::DAEProblem,
						  alg::BrownFullBasicInit, ::Val{false})
	@unpack p, t, f = integrator
	differential_vars = prob.differential_vars

	if integrator.opts.internalnorm(f(integrator.du, integrator.u, p, t),t) <= alg.abstol
		return
	elseif differential_vars === nothing
		error("differential_vars must be set for DAE initialization to occur. Either set consistent initial conditions, differential_vars, or use a different initialization algorithm.")
	end

	if integrator.u isa Number && integrator.du isa Number
		# This doesn't fix static arrays!
		u = [integrator.u]
		du = [integrator.du]
	else
		u = integrator.u
		du = integrator.du
	end

	nlequation = (out,x) -> begin
		@. du = ifelse(differential_vars,x,du)
		@. u  = ifelse(differential_vars,u,x)
		out .= f(du, u, p, t)
	end

	r = nlsolve(nlequation, ifelse.(differential_vars,du,u))

	@. du = ifelse(differential_vars,r.zero,du)
	@. u  = ifelse(differential_vars,u,r.zero)

	if integrator.u isa Number && integrator.du isa Number
		# This doesn't fix static arrays!
		integrator.u = first(u)
		integrator.du = first(du)
	else
		integrator.u = u
		integrator.du = du
	end

	integrator.uprev = integrator.u
	if alg_extrapolates(integrator.alg)
		integrator.uprev2 = integrator.uprev
	end

	return
end


text:
export sunrise, sunset

include("points.jl")

# The built-in deg2rad, rad2deg functions do
# weird things when z is very large
deg2radDP(z) = z * pi/180
rad2degDP(z) = z * 180/pi

isdaylight(dt::DateTime,p::Point,tz::TimeZone) = DateTime(sunrise(Date(dt),p,tz)) <= dt < DateTime(sunset(Date(dt),p,tz))

sunrise(d::Date,p::Point,tz::TimeZone) = calcSunriseSet(true,d,latitude(p),longitude(p),tz)
sunset(d::Date,p::Point,tz::TimeZone) = calcSunriseSet(false,d,latitude(p),longitude(p),tz)

function calcSunriseSet(rise,d::Date,latitude,longitude,tz)
    calcSunriseSet(rise,datetime2julian(DateTime(d)),latitude,longitude,tz)
end

function calcSunriseSet(rise,JD,latitude,longitude,tz)
    timeUTC = calcSunriseSetUTC(rise,JD,latitude,longitude)
    newtimeUTC = calcSunriseSetUTC(rise,JD+timeUTC/1440.0,latitude,longitude)
    dt = ZonedDateTime(julian2datetime(JD+newtimeUTC/1440.0),tz"UTC")
    astimezone(dt,tz)
end

function calcSunriseSetUTC(rise,JD,latitude,longitude)
    t = calcTimeJulianCent(JD)
    eqTime = calcEquationOfTime(t)
    solarDec = calcSunDeclination(t)
    hourAngle = calcHourAngleSunrise(latitude,solarDec)
    if !rise
        hourAngle *= -1
    end
    delta = longitude + rad2degDP(hourAngle)
    720  - (4.0 * delta) - eqTime
end

calcTimeJulianCent(jd) = (jd - 2451545.0)/36525.0

function calcEquationOfTime(t)
  epsilon = calcObliquityCorrection(t)
  l0 = calcGeomMeanLongSun(t)
  e = calcEccentricityEarthOrbit(t)
  m = calcGeomMeanAnomalySun(t)

  y = tan(deg2radDP(epsilon)/2.0);
  y *= y;

  sin2l0 = sin(2.0 * deg2radDP(l0))
  sinm   = sin(deg2radDP(m))
  cos2l0 = cos(2.0 * deg2radDP(l0))
  sin4l0 = sin(4.0 * deg2radDP(l0))
  sin2m  = sin(2.0 * deg2radDP(m))

  Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m
  rad2degDP(Etime)*4.0
end

function calcSunDeclination(t)
  e = calcObliquityCorrection(t)
  lambda = calcSunApparentLong(t)

  sint = sin(deg2radDP(e)) * sin(deg2radDP(lambda))
  rad2degDP(asin(sint))
end

function calcObliquityCorrection(t)
  e0 = calcMeanObliquityOfEcliptic(t)
  omega = 125.04 - 1934.136 * t
  e0 + 0.00256 * cosd(omega)
end

function calcGeomMeanLongSun(t)
  L0 = 280.46646 + t * (36000.76983 + t*(0.0003032))

  while L0 > 360.0
    L0 -= 360.0
  end

  while L0 < 0.0
    L0 += 360.0
  end

  L0
end

calcEccentricityEarthOrbit(t) = 0.016708634 - t * (0.000042037 + 0.0000001267 * t)

calcGeomMeanAnomalySun(t) = 357.52911 + t * (35999.05029 - 0.0001537 * t)

function calcSunApparentLong(t)
  o = calcSunTrueLong(t)
  omega = 125.04 - 1934.136 * t
  o - 0.00569 - 0.00478 * sind(omega)
end

function calcMeanObliquityOfEcliptic(t)
  seconds = 21.448 - t*(46.8150 + t*(0.00059 - t*(0.001813)))
  23.0 + (26.0 + (seconds/60.0))/60.0
end

function calcSunTrueLong(t)
  l0 = calcGeomMeanLongSun(t)
  c = calcSunEqOfCenter(t)
  l0 + c
end

function calcSunEqOfCenter(t)
  m = calcGeomMeanAnomalySun(t)
  mrad = deg2radDP(m)
  sinm = sin(mrad)
  sin2m = sin(2mrad)
  sin3m = sin(3mrad)
  C = sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) + sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289
end

function calcHourAngleSunrise(latitude,dec)
    latRad = deg2radDP(latitude)
    sdRad  = deg2radDP(dec)
    HAarg = (cos(deg2radDP(90.833))/(cos(latRad)*cos(sdRad))-tan(latRad) * tan(sdRad))
    HA = acos(HAarg)
end


text:
## PUT face database
## https://biometrics.cie.put.poznan.pl/index.php?option=com_content&view=article&id=4&Itemid=2&lang=en

immutable PutDataset <: FaceDataset
end

immutable PutFrontalDataset <: FaceDataset
end

function shape_to_image_path(shape_path)
    dir, filename = splitdir(shape_path)
    subjdir = basename(dir)
    basedir = dirname(dirname(dir))
    return joinpath(basedir, "images", "0" * subjdir[2:end],
                    replace(filename, ".yml", ".JPG"))
end

function load_images(::Type{PutFrontalDataset}, datadir::AbstractString;
                     downscale_times::Int=3)
    contourdir = joinpath(datadir, "contours")
    contourfiles = walkdir(contourdir)
    files = map(shape_to_image_path, contourfiles)
    function it()
        for file in files
            img = Images.load(file)
            # downscale if needed
            for i=1:downscale_times
                img = restrict(img)
            end
            arr = rgb2arr(img)
            produce(arr)
        end
    end
    return @task it()
end

function load_shape_put(path::AbstractString)
    # YAML parser couldn't parse it, fall back to a stupid manual method
    lines = split(open(readall, path), "\r\n")
    x_lines = filter(l -> startswith(l, "x"), map(strip, lines))
    xs = [parse(Float64, l[4:end]) for l in x_lines]
    y_lines = filter(l -> startswith(l, "y"), map(strip, lines))
    ys = [parse(Float64, l[4:end]) for l in y_lines]
    return hcat(ys, xs)  # ij coordinates
end

function load_shapes(::Type{PutFrontalDataset}, datadir::AbstractString;
                     downscale_times::Int=3)
    contourdir = joinpath(datadir, "contours")
    files = walkdir(contourdir)
    function it()
        for file in files
            shape = load_shape_put(file)
            downscaled = shape ./ (2^(downscale_times))
            produce(downscaled)
        end
    end
    return @task it()
end


text:
module TestPartitionLp
  using Base.Test
  using Brim

  A = [1 0; 0 1]
  S = [1 0; 0 1; 1 0; 0 1]

  M = partition_lp(A)

  # The modularity should be the diagonal matrix
  @test sum(abs(M.S .- S)) == 0

  # The original adjacency matrix should not change
  @test sum(abs(M.A .- A)) == 0

end