text:
function auglag_kernel(n::Int, major_iter::Int, max_auglag::Int,
                       line_start::Int, scale::Float64,
                       mu_max::Float64,
                       u_curr::CuDeviceArray{Float64,1}, v_curr::CuDeviceArray{Float64,1},
                       l_curr::CuDeviceArray{Float64,1}, rho::CuDeviceArray{Float64,1},
                       wRIij::CuDeviceArray{Float64,1},
                       param::CuDeviceArray{Float64,2},
                       _YffR::CuDeviceArray{Float64,1}, _YffI::CuDeviceArray{Float64,1},
                       _YftR::CuDeviceArray{Float64,1}, _YftI::CuDeviceArray{Float64,1},
                       _YttR::CuDeviceArray{Float64,1}, _YttI::CuDeviceArray{Float64,1},
                       _YtfR::CuDeviceArray{Float64,1}, _YtfI::CuDeviceArray{Float64,1},
                       frBound::CuDeviceArray{Float64,1}, toBound::CuDeviceArray{Float64,1})
    tx = threadIdx().x
    ty = threadIdx().y
    I = blockIdx().x

    x = @cuDynamicSharedMem(Float64, n)
    xl = @cuDynamicSharedMem(Float64, n, n*sizeof(Float64))
    xu = @cuDynamicSharedMem(Float64, n, (2*n)*sizeof(Float64))

    pij_idx = line_start + 8*(I-1)

    if tx == 1 && ty == 1
        @inbounds for i=1:n
            xl[i] = -Inf
            xu[i] = Inf
        end

        @inbounds begin
            xl[5] = frBound[2*(I-1)+1]
            xu[5] = frBound[2*I]
            xl[6] = toBound[2*(I-1)+1]
            xu[6] = toBound[2*I]
            xl[9] = -2*pi
            xu[9] = 2*pi
            xl[10] = -2*pi
            xu[10] = 2*pi

            x[1] = u_curr[pij_idx]
            x[2] = u_curr[pij_idx+1]
            x[3] = u_curr[pij_idx+2]
            x[4] = u_curr[pij_idx+3]
            x[5] = min(xu[5], max(xl[5], u_curr[pij_idx+4]))
            x[6] = min(xu[6], max(xl[6], u_curr[pij_idx+5]))
            x[7] = wRIij[2*(I-1)+1]
            x[8] = wRIij[2*I]
            x[9] = min(xu[9], max(xl[9], u_curr[pij_idx+6]))
            x[10] = min(xu[10], max(xl[10], u_curr[pij_idx+7]))
        end
    end

    @inbounds begin
        YffR = _YffR[I]; YffI = _YffI[I]
        YftR = _YftR[I]; YftI = _YftI[I]
        YttR = _YttR[I]; YttI = _YttI[I]
        YtfR = _YtfR[I]; YtfI = _YtfI[I]

        param[1,I] = l_curr[pij_idx]
        param[2,I] = l_curr[pij_idx+1]
        param[3,I] = l_curr[pij_idx+2]
        param[4,I] = l_curr[pij_idx+3]
        param[5,I] = l_curr[pij_idx+4]
        param[6,I] = l_curr[pij_idx+5]
        param[7,I] = rho[pij_idx]
        param[8,I] = rho[pij_idx+1]
        param[9,I] = rho[pij_idx+2]
        param[10,I] = rho[pij_idx+3]
        param[11,I] = rho[pij_idx+4]
        param[12,I] = rho[pij_idx+5]
        param[13,I] = v_curr[pij_idx]
        param[14,I] = v_curr[pij_idx+1]
        param[15,I] = v_curr[pij_idx+2]
        param[16,I] = v_curr[pij_idx+3]
        param[17,I] = v_curr[pij_idx+4]
        param[18,I] = v_curr[pij_idx+5]

        param[25,I] = l_curr[pij_idx+6]
        param[26,I] = l_curr[pij_idx+7]
        param[27,I] = rho[pij_idx+6]
        param[28,I] = rho[pij_idx+7]
        param[29,I] = v_curr[pij_idx+6]
        param[30,I] = v_curr[pij_idx+7]
    end

    if major_iter == 1
        @inbounds param[24,I] = 10.0
        mu = 10.0
    else
        @inbounds mu = param[24,I]
    end

    CUDA.sync_threads()

    eta = 1 / CUDA.pow(mu, 0.1)
    omega = 1 / mu
    max_feval = 500
    max_minor = 200
    gtol = 1e-6

    it = 0
    terminate = false

    while !terminate
        it += 1

        # Solve the branch problem.
        status, minor_iter = tron_kernel(n, 0, max_feval, max_minor, gtol, scale, false, x, xl, xu,
                                         param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI)

        @inbounds begin
            # Check the termination condition.
            cviol1 = x[1] - (YffR*x[5] + YftR*x[7] + YftI*x[8])
            cviol2 = x[2] - (-YffI*x[5] - YftI*x[7] + YftR*x[8])
            cviol3 = x[3] - (YttR*x[6] + YtfR*x[7] - YtfI*x[8])
            cviol4 = x[4] - (-YttI*x[6] - YtfI*x[7] - YtfR*x[8])
            cviol5 = x[7]^2 + x[8]^2 - x[5]*x[6]
            cviol6 = x[9] - x[10] - atan(x[8], x[7])
        end

        cnorm = max(abs(cviol1), abs(cviol2), abs(cviol3), abs(cviol4), abs(cviol5), abs(cviol6))

        if cnorm <= eta
            if cnorm <= 1e-6
                terminate = true
            else
                if tx == 1 && ty == 1
                    @inbounds begin
                        param[19,I] += mu*cviol1
                        param[20,I] += mu*cviol2
                        param[21,I] += mu*cviol3
                        param[22,I] += mu*cviol4
                        param[23,I] += mu*cviol5
                        param[31,I] += mu*cviol6
                    end
                end

                eta = eta / CUDA.pow(mu, 0.9)
                omega  = omega / mu
            end
        else
            mu = min(mu_max, mu*10)
            eta = 1 / CUDA.pow(mu, 0.1)
            omega = 1 / mu
            @inbounds param[24,I] = mu
        end

        if it >= max_auglag
            terminate = true
        end

        CUDA.sync_threads()
    end

    @inbounds begin
        u_curr[pij_idx] = x[1]
        u_curr[pij_idx+1] = x[2]
        u_curr[pij_idx+2] = x[3]
        u_curr[pij_idx+3] = x[4]
        u_curr[pij_idx+4] = x[5]
        u_curr[pij_idx+5] = x[6]
        wRIij[2*(I-1)+1] = x[7]
        wRIij[2*I] = x[8]
        u_curr[pij_idx+6] = x[9]
        u_curr[pij_idx+7] = x[10]
        param[24,I] = mu
    end

    CUDA.sync_threads()

    return
end

function auglag_kernel_cpu(n::Int, nline::Int, major_iter::Int, max_auglag::Int,
                           line_start::Int, mu_max::Float64,
                           u_curr::Array{Float64}, v_curr::Array{Float64},
                           l_curr::Array{Float64}, rho::Array{Float64},
                           wRIij::Array{Float64},
                           param::Array{Float64},
                           YffR::Array{Float64}, YffI::Array{Float64},
                           YftR::Array{Float64}, YftI::Array{Float64},
                           YttR::Array{Float64}, YttI::Array{Float64},
                           YtfR::Array{Float64}, YtfI::Array{Float64},
                           frBound::Array{Float64}, toBound::Array{Float64})
    avg_auglag_it = 0
    avg_minor_it = 0

    x = zeros(n)
    xl = zeros(n)
    xu = zeros(n)

    @inbounds for I=1:nline
        @inbounds for i=1:n
            xl[i] = -Inf
            xu[i] = Inf
        end

        xl[5] = frBound[2*(I-1)+1]
        xu[5] = frBound[2*I]
        xl[6] = toBound[2*(I-1)+1]
        xu[6] = toBound[2*I]
        xl[9] = -2*pi
        xu[9] = 2*pi
        xl[10] = -2*pi
        xu[10] = 2*pi

        pij_idx = line_start + 8*(I-1)

        x[1] = u_curr[pij_idx]
        x[2] = u_curr[pij_idx+1]
        x[3] = u_curr[pij_idx+2]
        x[4] = u_curr[pij_idx+3]
        x[5] = min(xu[5], max(xl[5], u_curr[pij_idx+4]))
        x[6] = min(xu[6], max(xl[6], u_curr[pij_idx+5]))
        x[7] = wRIij[2*(I-1)+1]
        x[8] = wRIij[2*I]
        x[9] = min(xu[9], max(xl[9], u_curr[pij_idx+6]))
        x[10] = min(xu[10], max(xl[10], u_curr[pij_idx+7]))

        param[1,I] = l_curr[pij_idx]
        param[2,I] = l_curr[pij_idx+1]
        param[3,I] = l_curr[pij_idx+2]
        param[4,I] = l_curr[pij_idx+3]
        param[5,I] = l_curr[pij_idx+4]
        param[6,I] = l_curr[pij_idx+5]
        param[7,I] = rho[pij_idx]
        param[8,I] = rho[pij_idx+1]
        param[9,I] = rho[pij_idx+2]
        param[10,I] = rho[pij_idx+3]
        param[11,I] = rho[pij_idx+4]
        param[12,I] = rho[pij_idx+5]
        param[13,I] = v_curr[pij_idx]
        param[14,I] = v_curr[pij_idx+1]
        param[15,I] = v_curr[pij_idx+2]
        param[16,I] = v_curr[pij_idx+3]
        param[17,I] = v_curr[pij_idx+4]
        param[18,I] = v_curr[pij_idx+5]

        if major_iter == 1
            param[24,I] = 10.0
            mu = 10.0
        else
            mu = param[24,I]
        end

        param[25,I] = l_curr[pij_idx+6]
        param[26,I] = l_curr[pij_idx+7]
        param[27,I] = rho[pij_idx+6]
        param[28,I] = rho[pij_idx+7]
        param[29,I] = v_curr[pij_idx+6]
        param[30,I] = v_curr[pij_idx+7]

        function eval_f_cb(x)
            f= eval_f_kernel_cpu(I, x, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI)
            return f
        end

        function eval_g_cb(x, g)
            eval_grad_f_kernel_cpu(I, x, g, param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI)
            return
        end

        function eval_h_cb(x, mode, rows, cols, scale, lambda, values)
            eval_h_kernel_cpu(I, x, mode, scale, rows, cols, lambda, values,
                              param, YffR, YffI, YftR, YftI, YttR, YttI, YtfR, YtfI)
            return
        end

        eta = 1 / mu^0.1
        omega = 1 / mu
        max_feval = 500
        max_minor = 100
        gtol = 1e-6

        nele_hess = 33
        tron = ExaTron.createProblem(n, xl, xu, nele_hess, eval_f_cb, eval_g_cb, eval_h_cb;
                                     :tol => gtol, :matrix_type => :Dense, :max_minor => 200,
                                     :frtol => 1e-12)
        tron.x .= x
        it = 0
        avg_tron_minor = 0
        terminate = false

        while !terminate
            it += 1

            # Solve the branch problem.
            status = ExaTron.solveProblem(tron)
            x .= tron.x
            avg_tron_minor += tron.minor_iter

            # Check the termination condition.
            cviol1 = x[1] - (YffR[I]*x[5] + YftR[I]*x[7] + YftI[I]*x[8])
            cviol2 = x[2] - (-YffI[I]*x[5] - YftI[I]*x[7] + YftR[I]*x[8])
            cviol3 = x[3] - (YttR[I]*x[6] + YtfR[I]*x[7] - YtfI[I]*x[8])
            cviol4 = x[4] - (-YttI[I]*x[6] - YtfI[I]*x[7] - YtfR[I]*x[8])
            cviol5 = x[7]^2 + x[8]^2 - x[5]*x[6]
            cviol6 = x[9] - x[10] - atan(x[8], x[7])

            cnorm = max(abs(cviol1), abs(cviol2), abs(cviol3), abs(cviol4), abs(cviol5), abs(cviol6))

            if cnorm <= eta
                if cnorm <= 1e-6
                    terminate = true
                else
                    param[19,I] += mu*cviol1
                    param[20,I] += mu*cviol2
                    param[21,I] += mu*cviol3
                    param[22,I] += mu*cviol4
                    param[23,I] += mu*cviol5
                    param[31,I] += mu*cviol6

                    eta = eta / mu^0.9
                    omega  = omega / mu
                end
            else
                mu = min(mu_max, mu*10)
                eta = 1 / mu^0.1
                omega = 1 / mu
                param[24,I] = mu
            end

            if it >= max_auglag
                terminate = true
            end
        end

        avg_auglag_it += it
        avg_minor_it += (avg_tron_minor / it)
        u_curr[pij_idx] = x[1]
        u_curr[pij_idx+1] = x[2]
        u_curr[pij_idx+2] = x[3]
        u_curr[pij_idx+3] = x[4]
        u_curr[pij_idx+4] = x[5]
        u_curr[pij_idx+5] = x[6]
        wRIij[2*(I-1)+1] = x[7]
        wRIij[2*I] = x[8]
        u_curr[pij_idx+6] = x[9]
        u_curr[pij_idx+7] = x[10]
        param[24,I] = mu
    end

    return (avg_auglag_it / nline), (avg_minor_it / nline)
end


text:
using SeisIO, SeisIO.RandSeis, Test

"""
    check_tutorial_data(S::GphysData; prune::Bool=true)

Check tutorial data for possible errors. If `prune=true`, channels with data
errors are deleted.

Current checks:
* No data
* Data are all NaNs (common with bad filtering parameters)
* Data are all equal (very rare; theoretically possible but we've never seen it)
"""
function check_tutorial_data!(S::GphysData; prune::Bool=true)
  pull = falses(S.n)

  for i in 1:S.n
    T = eltype(S.x[i])

    # check for no data
    if length(S.x[i]) == 0
      @warn(string("Channel ", i, " [id = ", S.id[i], "] has no data."))
      pull[i] = true
      continue
    end

    # check for NaNs
    nn = length(findall(isnan.(S.x[i])))
    if nn > 0
      @warn(string("Channel ", i,
            " [id = ", S.id[i], ", data type = ", T,
            "]. Output contains ", nn, " NaNs."))
      pull[i] = true
    end

    # check for variation
    mx = maximum(abs.(diff(S.x[i])))
    if mx == zero(T)
      @warn(string("Channel ", i,
            " [id = ", S.id[i], ", data type = ", T,
            "]. Output is a constant."))
      pull[i] = true
    end
  end

  if prune
    inds = findall(pull.==true)
    if isempty(inds)
      @info("No data issues found.")
    else
      info_strings = join([string(i, ": ", S.id[i]) for i in inds], "\n")
      @info(string("Deleting channels: \n", info_strings))
      deleteat!(S, inds)
    end
  end
  return nothing
end

# run a test to ensure this is behaving correctly
function test_check_tutorial_data()
  @info(string("Testing check_tutorial_data(). Expect two Warnings and an Info string."))
  S = randSeisData(6, nx=32000, s=1.0)
  U = deepcopy(S)
  S.x[1][12:20] .= NaN
  fill!(S.x[3], 3*one(eltype(S.x[3])))
  S.x[6] = eltype(S.x[6])[]
  check_tutorial_data!(S)
  @test S.id == U.id[[2, 4, 5]]
  @test S.x == U.x[[2, 4, 5]]
  println("Test of check_tutorial_data! complete.")
  return nothing
end

ext: jl


text:
#This is new, need this

function calculate_rules_vector(data_dictionary::Dict{String,Any}, value_array)

	# Initialize -
	v = zeros(95)

	function estimate_default(data_dictionary)
		return data_dictionary["default_vmax_value"]
	end
	default() = estimate_default(data_dictionary)

	function calculate_fn_value(x, value_array)
		return value_array[x]
	end
	fn(x) = calculate_fn_value(x, value_array)

	v[1] = max(fn(6), fn(35))
	v[2] = fn(137)
	v[3] = max(max(fn(97), fn(76)), fn(57))
	v[4] = max(fn(5), fn(36))
	v[5] = max(fn(5), fn(36))
	v[6] = default()
	v[7] = fn(9)
	v[8] = min(min(fn(4), fn(16)), fn(17))
	v[9] = fn(87)
	v[10] = max(max(fn(7), fn(39)), fn(35))
	v[11] = default()
	v[12] = max(min(min(min(fn(112), fn(113)), fn(114)), min(min(min(min(fn(107), fn(108)), fn(109)), fn(110)), fn(111))), min(min(min(min(fn(112), fn(113)), fn(114)), min(min(min(min(fn(107), fn(108)), fn(109)), fn(110)), fn(111))), fn(115)))
	v[13] = default()
	v[14] = fn(137)
	v[15] = fn(11)
	v[16] = max(min(fn(31), fn(32)), min(fn(20), fn(21)))
	v[17] = max(fn(93), fn(105))
	v[18] = fn(88)
	v[19] = default()
	v[20] = default()
	v[21] = default()
	v[22] = default()
	v[23] = default()
	v[24] = default()
	v[25] = default()
	v[26] = default()
	v[27] = default()
	v[28] = default()
	v[29] = default()
	v[30] = default()
	v[31] = default()
	v[32] = default()
	v[33] = default()
	v[34] = default()
	v[35] = default()
	v[36] = default()
	v[37] = default()
	v[38] = default()
	v[39] = default()
	v[40] = max(max(fn(61), fn(51)), fn(90))
	v[41] = max(fn(119), fn(134))
	v[42] = max(fn(30), fn(85))
	v[43] = max(fn(30), fn(85))
	v[44] = min(min(min(fn(130), fn(131)), fn(132)), fn(133))
	v[45] = min(min(min(min(fn(54), fn(55)), fn(56)), fn(78)), fn(79))
	v[46] = max(max(fn(45), fn(129)), fn(44))
	v[47] = fn(104)
	v[48] = fn(58)
	v[49] = fn(52)
	v[50] = max(max(min(min(min(fn(80), fn(33)), fn(78)), fn(79)), min(min(min(min(fn(54), fn(55)), fn(56)), fn(78)), fn(79))), min(min(min(fn(80), fn(46)), fn(78)), fn(79)))
	v[51] = max(fn(116), fn(37))
	v[52] = min(min(fn(26), fn(25)), fn(24))
	v[53] = fn(50)
	v[54] = max(max(fn(53), fn(10)), fn(41))
	v[55] = min(fn(98), fn(99))
	v[56] = fn(127)
	v[57] = fn(60)
	v[58] = max(fn(27), fn(137))
	v[59] = fn(34)
	v[60] = fn(125)
	v[61] = max(fn(62), fn(38))
	v[62] = max(fn(124), fn(94))
	v[63] = fn(104)
	v[64] = fn(100)
	v[65] = fn(40)
	v[66] = fn(82)
	v[67] = min(min(min(min(min(min(min(min(min(min(min(min(fn(63), fn(64)), fn(65)), fn(66)), fn(67)), fn(68)), fn(69)), fn(70)), fn(71)), fn(72)), fn(73)), fn(74)), fn(75))
	v[68] = max(fn(123), min(fn(42), fn(43)))
	v[69] = max(fn(137), fn(8))
	v[70] = fn(137)
	v[71] = min(min(fn(2), fn(3)), fn(4))
	v[72] = max(fn(117), fn(49))
	v[73] = max(max(max(min(min(fn(28), fn(29)), fn(86)), min(fn(28), fn(29))), min(fn(28), fn(96))), min(fn(120), fn(121)))
	v[74] = fn(126)
	v[75] = fn(91)
	v[76] = fn(23)
	v[77] = max(max(fn(106), fn(136)), fn(22))
	v[78] = max(fn(95), fn(103))
	v[79] = fn(122)
	v[80] = fn(102)
	v[81] = fn(48)
	v[82] = max(fn(77), fn(81))
	v[83] = max(fn(59), fn(47))
	v[84] = default()
	v[85] = max(fn(101), fn(135))
	v[86] = max(fn(89), fn(128))
	v[87] = fn(104)
	v[88] = default()
	v[89] = min(min(min(fn(12), fn(13)), fn(14)), fn(15))
	v[90] = min(fn(18), fn(19))
	v[91] = max(fn(83), fn(1))
	v[92] = min(fn(42), fn(43))
	v[93] = max(fn(92), fn(84))
	v[94] = max(fn(92), fn(84))
	v[95] = fn(118)

	#v = v .== 1

	return v
end

function calculate_gene_rules_vector(data_dictionary::Dict{String,Any}, value_array::BitArray{1})

	# Initialize -
	v = zeros(95)

	function estimate_default(data_dictionary)
		return data_dictionary["default_vmax_value"]
	end
	default() = true

	function calculate_fn_value(x, value_array)
		return value_array[x]
	end
	fn(x) = calculate_fn_value(x, value_array)

	v[1] = max(fn(6), fn(35))
	v[2] = fn(137)
	v[3] = max(max(fn(97), fn(76)), fn(57))
	v[4] = max(fn(5), fn(36))
	v[5] = max(fn(5), fn(36))
	v[6] = default()
	v[7] = fn(9)
	v[8] = min(min(fn(4), fn(16)), fn(17))
	v[9] = fn(87)
	v[10] = max(max(fn(7), fn(39)), fn(35))
	v[11] = default()
	v[12] = max(min(min(min(fn(112), fn(113)), fn(114)), min(min(min(min(fn(107), fn(108)), fn(109)), fn(110)), fn(111))), min(min(min(min(fn(112), fn(113)), fn(114)), min(min(min(min(fn(107), fn(108)), fn(109)), fn(110)), fn(111))), fn(115)))
	v[13] = default()
	v[14] = fn(137)
	v[15] = fn(11)
	v[16] = max(min(fn(31), fn(32)), min(fn(20), fn(21)))
	v[17] = max(fn(93), fn(105))
	v[18] = fn(88)
	v[19] = default()
	v[20] = default()
	v[21] = default()
	v[22] = default()
	v[23] = default()
	v[24] = default()
	v[25] = default()
	v[26] = default()
	v[27] = default()
	v[28] = default()
	v[29] = default()
	v[30] = default()
	v[31] = default()
	v[32] = default()
	v[33] = default()
	v[34] = default()
	v[35] = default()
	v[36] = default()
	v[37] = default()
	v[38] = default()
	v[39] = default()
	v[40] = max(max(fn(61), fn(51)), fn(90))
	v[41] = max(fn(119), fn(134))
	v[42] = max(fn(30), fn(85))
	v[43] = max(fn(30), fn(85))
	v[44] = min(min(min(fn(130), fn(131)), fn(132)), fn(133))
	v[45] = min(min(min(min(fn(54), fn(55)), fn(56)), fn(78)), fn(79))
	v[46] = max(max(fn(45), fn(129)), fn(44))
	v[47] = fn(104)
	v[48] = fn(58)
	v[49] = fn(52)
	v[50] = max(max(min(min(min(fn(80), fn(33)), fn(78)), fn(79)), min(min(min(min(fn(54), fn(55)), fn(56)), fn(78)), fn(79))), min(min(min(fn(80), fn(46)), fn(78)), fn(79)))
	v[51] = max(fn(116), fn(37))
	v[52] = min(min(fn(26), fn(25)), fn(24))
	v[53] = fn(50)
	v[54] = max(max(fn(53), fn(10)), fn(41))
	v[55] = min(fn(98), fn(99))
	v[56] = fn(127)
	v[57] = fn(60)
	v[58] = max(fn(27), fn(137))
	v[59] = fn(34)
	v[60] = fn(125)
	v[61] = max(fn(62), fn(38))
	v[62] = max(fn(124), fn(94))
	v[63] = fn(104)
	v[64] = fn(100)
	v[65] = fn(40)
	v[66] = fn(82)
	v[67] = min(min(min(min(min(min(min(min(min(min(min(min(fn(63), fn(64)), fn(65)), fn(66)), fn(67)), fn(68)), fn(69)), fn(70)), fn(71)), fn(72)), fn(73)), fn(74)), fn(75))
	v[68] = max(fn(123), min(fn(42), fn(43)))
	v[69] = max(fn(137), fn(8))
	v[70] = fn(137)
	v[71] = min(min(fn(2), fn(3)), fn(4))
	v[72] = max(fn(117), fn(49))
	v[73] = max(max(max(min(min(fn(28), fn(29)), fn(86)), min(fn(28), fn(29))), min(fn(28), fn(96))), min(fn(120), fn(121)))
	v[74] = fn(126)
	v[75] = fn(91)
	v[76] = fn(23)
	v[77] = max(max(fn(106), fn(136)), fn(22))
	v[78] = max(fn(95), fn(103))
	v[79] = fn(122)
	v[80] = fn(102)
	v[81] = fn(48)
	v[82] = max(fn(77), fn(81))
	v[83] = max(fn(59), fn(47))
	v[84] = default()
	v[85] = max(fn(101), fn(135))
	v[86] = max(fn(89), fn(128))
	v[87] = fn(104)
	v[88] = default()
	v[89] = min(min(min(fn(12), fn(13)), fn(14)), fn(15))
	v[90] = min(fn(18), fn(19))
	v[91] = max(fn(83), fn(1))
	v[92] = min(fn(42), fn(43))
	v[93] = max(fn(92), fn(84))
	v[94] = max(fn(92), fn(84))
	v[95] = fn(118)

	#v = v .== 1

	return v
end

function calculate_full_rules_vector(data_dictionary::Dict{String,Any}, value_array::BitArray{1})

	# Initialize -
	v = zeros(159)

	function estimate_default(data_dictionary)
		return 1.0
	end
	default() = true

	function calculate_fn_value(x, value_array)
		return value_array[x]
	end
	fn(x) = calculate_fn_value(x, value_array)

	v[1] = default()
	v[2] = ~(fn(158))
	v[3] = ~(fn(159))
	v[4] = max(~(fn(3)), fn(110))
	v[5] = max(~(fn(3)), fn(110))
	v[6] = default()
	v[7] = default()
	v[8] = default()
	v[9] = default()
	v[10] = fn(11)
	v[11] = ~(fn(169))
	v[12] = default()
	v[13] = default()
	v[14] = default()
	v[15] = default()
	v[16] = max(max(~(max(fn(152), fn(43))), fn(111)), fn(110))
	v[17] = max(max(~(max(fn(152), fn(43))), fn(111)), fn(110))
	v[18] = max(max(~(max(fn(152), fn(43))), fn(111)), fn(110))
	v[19] = max(max(~(max(fn(152), fn(43))), fn(111)), fn(110))
	v[20] = default()
	v[21] = default()
	v[22] = default()
	v[23] = default()
	v[24] = max(~(fn(43)), fn(152))
	v[25] = max(~(fn(43)), fn(152))
	v[26] = default()
	v[27] = default()
	v[28] = default()
	v[29] = default()
	v[30] = default()
	v[31] = default()
	v[32] = max(fn(152), min(fn(43), fn(111)))
	v[33] = max(fn(152), min(fn(43), fn(111)))
	v[34] = max(fn(152), min(fn(43), fn(111)))
	v[35] = default()
	v[36] = default()
	v[37] = max(~(fn(48)), ~(fn(2)))
	v[38] = default()
	v[39] = max(fn(163), ~(fn(160)))
	v[40] = max(max(~(fn(168)), ~(min(fn(168), fn(2)))), fn(110))
	v[41] = default()
	v[42] = default()
	v[43] = ~(fn(168))
	v[44] = default()
	v[45] = default()
	v[46] = default()
	v[47] = max(~(fn(163)), min(fn(167), ~(fn(111))))
	v[48] = ~(fn(163))
	v[49] = default()
	v[50] = default()
	v[51] = ~(fn(152))
	v[52] = ~(max(fn(152), fn(43)))
	v[53] = default()
	v[54] = ~(fn(2))
	v[55] = fn(2)
	v[56] = default()
	v[57] = ~(max(fn(67), fn(164)))
	v[58] = default()
	v[59] = default()
	v[60] = default()
	v[61] = max(fn(155), ~(fn(48)))
	v[62] = max(fn(155), ~(fn(48)))
	v[63] = max(fn(155), ~(fn(48)))
	v[64] = default()
	v[65] = default()
	v[66] = default()
	v[67] = fn(157)
	v[68] = default()
	v[69] = default()
	v[70] = default()
	v[71] = ~(max(fn(152), fn(43)))
	v[72] = ~(max(fn(152), fn(43)))
	v[73] = ~(max(fn(152), fn(43)))
	v[74] = ~(max(fn(152), fn(43)))
	v[75] = ~(max(fn(152), fn(43)))
	v[76] = ~(max(fn(152), fn(43)))
	v[77] = ~(max(fn(152), fn(43)))
	v[78] = ~(max(fn(152), fn(43)))
	v[79] = ~(max(fn(152), fn(43)))
	v[80] = ~(max(fn(152), fn(43)))
	v[81] = ~(max(fn(152), fn(43)))
	v[82] = ~(max(fn(152), fn(43)))
	v[83] = ~(max(fn(152), fn(43)))
	v[84] = default()
	v[85] = default()
	v[86] = default()
	v[87] = default()
	v[88] = default()
	v[89] = default()
	v[90] = default()
	v[91] = default()
	v[92] = default()
	v[93] = max(fn(152), min(fn(43), fn(111)))
	v[94] = default()
	v[95] = default()
	v[96] = default()
	v[97] = default()
	v[98] = default()
	v[99] = default()
	v[100] = default()
	v[101] = min(~(fn(152)), fn(103))
	v[102] = min(~(fn(152)), fn(103))
	v[103] = fn(160)
	v[104] = ~(fn(10))
	v[105] = max(fn(111), fn(43))
	v[106] = max(fn(111), fn(43))
	v[107] = ~(min(fn(156), fn(164)))
	v[108] = ~(min(fn(156), fn(164)))
	v[109] = ~(fn(152))
	v[110] = fn(171)
	v[111] = fn(154)
	v[112] = default()
	v[113] = default()
	v[114] = default()
	v[115] = min(min(fn(155), ~(fn(152))), fn(143))
	v[116] = ~(fn(152))
	v[117] = default()
	v[118] = default()
	v[119] = default()
	v[120] = default()
	v[121] = default()
	v[122] = default()
	v[123] = default()
	v[124] = default()
	v[125] = default()
	v[126] = default()
	v[127] = ~(fn(167))
	v[128] = fn(111)
	v[129] = default()
	v[130] = default()
	v[131] = default()
	v[132] = max(fn(152), fn(43))
	v[133] = max(fn(152), fn(43))
	v[134] = default()
	v[135] = default()
	v[136] = min(~(fn(138)), max(~(fn(152)), fn(2)))
	v[137] = min(~(fn(138)), max(~(fn(152)), fn(2)))
	v[138] = fn(39)
	v[139] = default()
	v[140] = default()
	v[141] = default()
	v[142] = max(max(fn(43), fn(111)), fn(143))
	v[143] = fn(144)
	v[144] = max(max(fn(170), fn(162)), fn(166))
	v[145] = max(fn(43), fn(143))
	v[146] = max(fn(43), fn(143))
	v[147] = max(fn(43), fn(143))
	v[148] = max(fn(43), fn(143))
	v[149] = default()
	v[150] = default()
	v[151] = default()
	v[152] = ~(fn(168))
	v[153] = default()
	v[154] = ~(fn(163))
	v[155] = ~(max(max(fn(163), fn(166)), fn(165)))
	v[156] = fn(157)
	v[157] = fn(127)
	v[158] = max(min(~(fn(172)), ~(max(max(fn(182), fn(181)), fn(178)))), fn(161))
	v[159] = min(~(max(fn(176), fn(175))), ~(max(max(max(max(fn(173), fn(179)), fn(177)), fn(174)), fn(180))))

	return v
end


text:
using CanEvolve
using Test

@testset "CanEvolve.jl" begin
	p, n = autoSimulation(log(2), 0.0, Nfinal = 1000, noise = "betabinomial", rho = 0.002, nsubclone_min = 1, nsubclone_max = 2, lower_cutoff = 0.1, upper_cutoff = 0.4)
	@test sum(p[3][1]) > 0.0
	@test sum(n[3][1]) == 0.0
end


text:
module BioMedQuery

# --------Database Utilities---------
module DBUtils
export init_mysql_database,
       assemble_vals,
       insert_row!,
       db_select,
       db_query,
       db_clean_string

include("DBUtils/DBUtils.jl")
end

#-----------Entrez----------------
module Entrez
export esearch,
       efetch,
       eparse,
       eparse_from_file,
       save_efetch_mysql,
       save_efetch_sqlite,
       citations_endnote,
       citations_bibtex,
       save_article_citations

#types
export PubMedArticle,
       MeshHeading,
       MeshHeadingList

include("Entrez/Entrez.jl")
include("Entrez/pubmed_article.jl")
include("Entrez/entrez_save.jl")
include("Entrez/citation_manager.jl")
end

#----------UMLS--------------------
module UMLS
export Credentials,
       get_tgt,
       search_umls,
       best_match_cui,
       get_semantic_type,
       populate_net_mysql,
       build_tree_dict
include("UMLS/UMLS.jl")
include("UMLS/semantic_network.jl")
end

#--------Clinical Trials------------
module CT
export search_ct
include("CT/CT.jl")
end
#-----------------------------------

#--------Processes------------
module Processes
export pubmed_search_and_save,
       map_mesh_to_umls_async!,
       map_mesh_to_umls!,
       umls_semantic_occurrences,
       filter_mesh_by_concept,
       export_citation
include("Processes/pubmed_search_and_save.jl")
include("Processes/pubmed_mesh_to_umls_map.jl")
include("Processes/pubmed_occurrance_filtering.jl")
include("Processes/pubmed_export_citations.jl")
end

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

end #BioMedQuery


text:
function nlp_cvx_109_012(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    model = Model(optimizer)

    @variable(model, x >= 0.00001, start = 0.1)
    @variable(model, y >= 0.00001, start = 0.1)

    @NLobjective(model, Max, log(x + y))
    @NLconstraint(model, (y - 2)^2 <= -x + 2)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, log(7 / 4 + 5 / 2), tol = objective_tol)
    return check_solution([x, y], [7 / 4, 5 / 2], tol = primal_tol)
end


text:
using Printf
function traverse(x)
    ceiling = trunc(Int, sqrt(x))
    for i in 3:2:ceiling
        if x % i == 0
            return false
        end
    end
    return true
end

function is_prime(x)
    if x == 2
        return true
    else
        if x % 2 == 0
            return false
        else
            return traverse(x)
        end
    end
end

print(is_prime(parse(Int, ARGS[1])))

max_stars_repo_name: KnugiHK/SingleScript

max_issues_repo_name: KnugiHK/SingleScript

max_forks_repo_name: KnugiHK/SingleScript


text:
## Convert julia mmd/md file into an Rmd file with julia code prepocessed
## TODO

##            code eval output
## .           o   o    o
## nocode      x   o    o
## noeval      o   x    x       (verbatim)
## noout       o   o    x
## Idea is we can mark our code fences as above and suppress the code, output or evaluation

## Helpers
rmd_code_chunk(buf, txt) =  println(buf, """

```{r eval=FALSE}
$txt

```
                                    """)

function rmd_output_chunk(buf, res)
    if VERSION >= v"0.5.0"
        tmpbuf = IOBuffer()
        context = IOContext(tmpbuf,multiline=true,limit=true)
        show(context, bestmime(res), res)
        out = String(take!(tmpbuf))
    else
        out = sprint(io -> show(io,  bestmime(res), res))
    end
    out = join(["## "*r for r in split(out, "\n")], "\n")

    println(buf, """

```{r eval=FALSE}
$out
```
""")
end

## for raw html, we need to strip out initial white space
function trim_white(tmp)
    tmp = split(tmp, "\n")
    tmp = join(map(lstrip, tmp), "\n")
end

## Main function to take a jmd file and turn into an Rmd file
function markdownToRmd(fname::AbstractString; TITLE="", kwargs...)

    dirnm, basenm = dirname(fname), basename(fname)
    newnm = replace(fname, r"[.].*" => ".Rmd")
    out = mdToRmd(fname; TITLE=TITLE, kwargs...)

    io = open(newnm, "w")
    write(io, out)
    close(io)
end

"""

Convert a markdown file into an Rmd file.

Processes `Question` types, graphical output of `Winston`, `Gadfly`,
and `PyPlot`. For the latter, one should call `gcf()` as the last
expressions in a cell block.

Markdown idiosyncracies:
* there is no `_underscore_` or `__double underscore__`
* for sections, best not to skip 2 or more lines, as they can get wrapped in `<p></p>` tags and not show
* Can use LaTeX markup

"""
function mdToRmd(fname::AbstractString; TITLE="", kwargs...)
    m = make_module()
    safeeval(m, Meta.parse("using LaTeXStrings, Plots; plotly()"))
    safeeval(m, Meta.parse("macro q_str(x)  \"`\$x`\" end"))

    process_block("using WeavePynb, LaTeXStrings", m)

    buf = IOBuffer()

    out = Markdown.parse_file(fname, flavor=Markdown.julia)

    for i in 1:length(out.content)
        ##println(out.content[i])
        if isa(out.content[i], Markdown.Code)
            ## Code Blocks are evaluated and their last value is added to the output
            ## this is different from IJulia, but similar.
            ## There are issues with Gadfly graphics (need the script...)
            ## and PyPlot, where we need an invocation to manage the figures

            txt = out.content[i].code
            lang = out.content[i].language

            ## we need to set
            ## nocode, noeval, noout
            langs = map(lstrip, split(lang, ","))

            docode, doeval, doout = true, true, true
            if "nocode" in langs
                docode = false
            end
            if "verbatim" in langs || "noeval" in langs
                doeval, doout = false, false
            end
            if "noout" in langs
                doout = false
            end

            ## language is used to pass in arguments
            result = nothing
            if doeval
                result = process_block(txt, m)
            end

            !docode && (txt = "")

            if isa(result, Outputonly)
                txt = ""
                result = result.x
            end

            ## special cases
            ## hsould use dispatch here, but we don't....
            if result == nothing
                "Do not show output, just input"
                txt = replace(txt, r"\nnothing$" => "") ## trim off trailing "nothing"
                docode && length(txt) > 0 && rmd_code_chunk(buf, txt)
            elseif isa(result, Invisible)
                "Do not show output or input"
                nothing
            elseif isa(result, HTMLoutput)
                "Do not execute input, show as is"
                txt = ""
                tmp = sprint(io -> show(io, "text/plain", result))
                tmp = trim_white(tmp)
                doout && println(buf, tmp) #show(buf, "text/plain", result)
            elseif isa(result, Verbatim)
                "Do not execute input, show as is"
                doout && show(buf, "text/plain", result)
            elseif isa(result, Bootstrap)
                "Show with Bootstrap formatting"
                doout && show(buf, "text/html", result)
            elseif isa(result, Question)
                "Show a question"
                tmp = sprint(io -> show(io, "text/html", result))
                tmp = trim_white(tmp)
                doout && println(buf, "\n", "----\n", "\n", tmp, "\n", "----\n", "\n")

            elseif isa(result, ImageFile)
                "show an image stored in a file, but embed"
                txt = ""
                img = gif_to_data(result.f, result.caption)
                doout && println(buf, trim_white(img))
            elseif isa(result, Plots.Plot)
                docode && length(txt) > 0 && rmd_code_chunk(buf, txt)

                if isa(result, Plots.Plot{Plots.PlotlyBackend})
                    Plots.prepare_output(result);
                    img= Plots.html_body(result)
                    img = trim_white(img)
                    doout && write(buf,  img)
                else
                    #                    img = stringmime("image/png", result.o)
                    imgfile = tempname() * ".png"
                    png(result, imgfile)
                    img = base64encode(readall(imgfile))
                    doout && println(buf, """""")
                end
            elseif string(typeof(result)) == "FramedPlot"
                length(txt) > 0 &&  rmd_code_chunk(buf, txt)
                img = stringmime("image/png", result)
                doout && println(buf, """""")
            elseif string(typeof(result)) == "Figure"
                "Must do gcf() for last line"
                if length(txt) > 0
                    txt1 = join(split(txt, "\n")[1:(end-1)], "\n") ## trim last line which is gcf()
                    docode &&   rmd_code_chunk(buf, txt)
                end
                img = stringmime("image/png", result)
                println(buf, """""")
            elseif string(typeof(result)) == "SymPy.Sym"
                length(txt) > 0 &&   rmd_code_chunk(buf, txt)

                rmd_output_chunk(buf, "... bug ...")
            else

                length(txt) > 0 &&   rmd_code_chunk(buf, txt)
                try
                    tmpbuf = IOBuffer()
                    if doout
                        rmd_output_chunk(tmpbuf, result)
                        println(buf, String(take!(tmpbuf)))
                    end
                catch e
                                        rethrow(e)
                    ## no output
                end
            end

        else
           if isa(out.content[i], Markdown.LaTeX)
               Markdown.latex(buf, out.content[i])
           else
               ## we drill down
               #               println(sprint(io -> Markdown.plain(io, out.content[i])))
               println(buf, "")
               Markdown.plain(buf,  out.content[i])
           end
        end

    end

    body = String(take!(buf)) #takebuf_string(buf)

    ## return string
    D = Dict()
D[:Title] = TITLE
D[:day] = string(today())
    D[:style] = ""
    D[:body] = body
    for (k,v) in kwargs
        D[k] = v
    end
    Mustache.render(markdown_tpl, D)
end

"""
  Helper function to allow the md file to be a Mustache template.
  Expects a file `fname.jl` which defines a module `fname`.

  XXX Needs work XXX
"""
function mmd_to_rmd(fname::AbstractString; kwargs...)
    mmd_to_md(fname)

    bname = basename(fname)
    bname = replace(bname, r"\.mmd$" => "")
    markdownToRmd("$bname.md"; kwargs...)
end
export mmd_to_rmd

"""
Template for an HTML page with embedded questions.

Takes the following  arguments in template for Mustache

:style -- optional additional style values
:BRAND_HREF -- for upper left corner of nav bar
:BRAND_NAME
:body -- filled in

Uses Bootstrap styling and MathJaX for LaTeX markup.

"""
markdown_tpl = mt"""
{{{:body}}}
"""

# tufte = """
# title: "{{:Title}}"
# subtitle: "{{:SubTitle}}"
# author: "{{:Author}}"
# date: {{:day}}
# output:
#   tufte::tufte_html: default
#   tufte::tufte_handout:
#     citation_package: natbib
#     latex_engine: xelatex
#   tufte::tufte_book:
#     citation_package: natbib
#     latex_engine: xelatex
# {{#:bibfile}}bibliography: {{:bibfile}}.bib {{/:bibfile}}
# link-citations: yes
# """


text:
using Documenter, NeuralNetDiffEq

makedocs(
    sitename = "NeuralNetDiffEq.jl",
    authors="#",
    clean = true,
    doctest = false,
    modules = [NeuralNetDiffEq],

    format = Documenter.HTML(#analytics = "",
                             assets = ["assets/favicon.ico"],
                             canonical="#"),
    pages=[
        "Home" => "index.md",
        "Tutorials" => Any[
        ],
    ]
)

deploydocs(
   repo = "github.com/SciML/NeuralNetDiffEq.jl.git";
   push_preview = true
)


text:
# deprecated
function simulate(m::HeadTurn;
    onset::Real = 0,
    trial_duration::Real = 5,
    stochastic::Bool = false)

    @assert onset >= 0
    @assert trial_duration > 0

    # onset = m.onset
    duration = m.duration
    total_duration = trial_duration

    onset + duration <= total_duration || error("Movement extends beyond simulation event.")
    m.D ∈ ["left", "right"] || error("Direction not specified.")

    D = m.D == "left" ? -1 : 1
    # motion amplitude is noisy if stochastic simulation
    A = stochastic == 0 ? m.A : rand(Normal(m.A, 0.5))
    noise = m.sensor
    Δt = m.Δt
    f = 1/(duration) # frequency: single sinusiodal head turn

    timesteps = range(0, stop = total_duration, step = Δt)
    T = length(timesteps)

    κ1, κ2 = cupula_dynamics(Δt)
    # α = zeros(T)
    # ω = Array{Float64}(undef, T) # velocity
    # θ = Array{Float64}(undef, T) # position
    # c = Array{Float64}(undef, T) # cupula dynamics
    # y = Array{Float64}(undef, T) # sensory observations
    α = zeros(T)
    ω = zeros(T)
    θ = zeros(T)
    c = zeros(T)
    y = zeros(T)

    # start until motion onset
    j = 1
    t = round(timesteps[j]; digits = 3)

    while t < onset
        y[j] = ω[j] - c[j] + rand(noise)
        j += 1
        t = round(timesteps[j]; digits = 3)
    end

    for i ∈ (j+1):T
        t = round(timesteps[i]; digits = 3)
        # here, we add a small amount of noise to the acceleration if the
        # simulation is stochastic
        α[i] = t < onset + duration ? acceleration(D, A, f, t, onset, noisy = stochastic) : 0

        ω[i] = ω[i-1] + Δt * α[i-1]
        θ[i] = θ[i-1] + Δt * ω[i-1] + 1/2 * Δt^2 * α[i-1]
        c[i] = κ1 * c[i-1] + κ2 * ω[i]
        y[i] = ω[i] - c[i] + rand(noise)
    end
    out = (timesteps = collect(timesteps),
            y = y, α = α, ω = ω, θ = θ, c = c,
            Δt = Δt, onset = onset)
    return out
end

# deprecated. TODO: add deceleration
function simulate(m::ExternalDisturbance;
    onset::Real = 0, trial_duration::Real = 5)

    @assert onset >= 0
    @assert trial_duration > 0

    # onset = m.onset
    duration = m.duration
    total_duration = trial_duration

    onset + duration <= total_duration || error("Movement extends beyond simulation event.")
    m.D ∈ ["left", "right"] || error("Direction not specified.")

    D = m.D == "left" ? -1 : 1
    noise = m.sensor
    Δt = m.Δt
    f = 1/(duration) # frequency: single sinusiodal head turn

    timesteps = range(0, stop = total_duration, step = Δt)
    T = length(timesteps)

    κ1, κ2 = cupula_dynamics(Δt)
    # α = zeros(T)
    # ω = Array{Float64}(undef, T) # velocity
    # θ = Array{Float64}(undef, T) # position
    # c = Array{Float64}(undef, T) # cupula dynamics
    # y = Array{Float64}(undef, T) # sensory observations
    α = zeros(T)
    ω = zeros(T)
    θ = zeros(T)
    c = zeros(T)
    y = zeros(T)

    # start until motion onset
    j = 1
    t = round(timesteps[j]; digits = 3)

    while t < onset
        y[j] = ω[j] - c[j] + rand(noise)
        j += 1
        t = round(timesteps[j]; digits = 3)
    end

    for i ∈ (j+1):T
        t = round(timesteps[i]; digits = 3)
        α[i] = t < onset + duration ? D * m.A : 0

        ω[i] = ω[i-1] + Δt * α[i-1]
        θ[i] = θ[i-1] + Δt * ω[i-1] + 1/2 * Δt^2 * α[i-1]
        c[i] = κ1 * c[i-1] + κ2 * ω[i]
        y[i] = ω[i] - c[i] + rand(noise)
    end
    out = (timesteps = collect(timesteps),
            y = y, α = α, ω = ω, θ = θ, c = c,
            Δt = Δt, onset = onset)
    return out
end

# deprecated
function simulate(mᵤ::HeadTurn,
    mₑ::HeadTurn;
    onsetᵤ::Real = 0, onsetₑ::Real = 0,
    trial_duration::Real = 5,
    stochastic::Bool = false)

    @assert onsetᵤ >= 0
    @assert onsetₑ >= 0
    @assert trial_duration > 0

    # onset = m.onset
    duration1 = mᵤ.duration
    duration2 = mₑ.duration
    total_duration = trial_duration

    onsetᵤ + duration1 <= total_duration || error("Head turn beyond simulation event.")
    onsetₑ + duration2 <= total_duration || error("Disturbance extends beyond simulation event.")

    mᵤ.D ∈ ["left", "right"] || error("Direction not specified.")
    mₑ.D ∈ ["left", "right"] || error("Direction not specified.")

    D1 = mᵤ.D == "left" ? -1 : 1
    D2 = mₑ.D == "left" ? -1 : 1

    # motion amplitude is noisy if stochastic simulation
    A1 = stochastic == 0 ? mᵤ.A : rand(Normal(m1.A, 0.5))
    A2 = stochastic == 0 ? mₑ.A : rand(Normal(m2.A, 0.5))

    noise = mᵤ.sensor
    Δt = mᵤ.Δt
    f1 = 1/(duration1) # frequency: single sinusiodal head turn
    f2 = 1/(duration2)

    timesteps = range(0, stop = total_duration, step = Δt)
    T = length(timesteps)

    κ1, κ2 = cupula_dynamics(Δt)
    # α = zeros(T)
    # ω = Array{Float64}(undef, T) # velocity
    # θ = Array{Float64}(undef, T) # position
    # c = Array{Float64}(undef, T) # cupula dynamics
    # y = Array{Float64}(undef, T) # sensory observations
    α1 = zeros(T)
    α2 = zeros(T)
    α = zeros(T)
    ω1 = zeros(T)
    ω2 = zeros(T)
    ω = zeros(T)
    θ = zeros(T)
    c = zeros(T)
    y = zeros(T)

    for i ∈ 2:T

        t = round(timesteps[i]; digits = 3)

        α1[i] = (t > onsetᵤ) & (t < onsetᵤ + duration1) ? acceleration(D1, A1, f1, t, onsetᵤ, noisy = stochastic) : 0
        α2[i] = (t > onsetₑ) & (t < onsetₑ + duration2) ? acceleration(D2, A2, f2, t, onsetₑ, noisy = stochastic) : 0
        α[i] = α1[i] + α2[i]
        ω1[i] = ω1[i-1] + Δt * α1[i-1]
        ω2[i] = ω2[i-1] + Δt * α2[i-1]

        ω[i] =  ω1[i] + ω2[i]
        θ[i] = θ[i-1] + Δt * ω[i-1] + 1/2 * Δt^2 * α[i-1]
        c[i] = κ1 * c[i-1] + κ2 * ω[i]
        y[i] = ω[i] - c[i] + rand(noise)
    end

    out = (timesteps = collect(timesteps),
            y = y,
            αᵤ = α1, αₑ = α2,
            α = α,
            ωᵤ = ω1, ωₑ = ω2,
            ω = ω,
            θ = θ, c = c,
            Δt = Δt, onsetᵤ = onsetᵤ, onsetₑ = onsetₑ)
    return out

end


text:
module NNSparseCoding
using Random,Distributions
using LinearAlgebra
using Images, FileIO
using MultivariateStats

struct NNSC_optimizer{T}
  X::Matrix{T}
  AS::Matrix{T}
  A::Matrix{T}
  S::Matrix{T}
  lambda::T
end

function NNSC_optimizer(X::Matrix{<:Real},d::Integer,
    lambda::Real)
  nfull,nsampl=size(X)
  A,S = initialize_mats(X,d)
  project!(A)
  NNSC_optimizer(X,A*S,A,S,lambda)
end
function nnsc_objective(o::NNSC_optimizer)
  return 0.5 * mapreduce( (x,x̃) -> (x-x̃)^2,
    +,o.AS,o.X) + o.lambda * sum(o.S)
end

function nnsc_objective(X::M,A::M,S::M,lambda::Real) where M <: Matrix{<:Real}
  X̃ = A*S
  @. X̃ = (X̃ - X)^2
  return 0.5*sum(X̃) + lambda * sum(S)
end

function project!(A::Matrix{<:Real})
  # no negative values
  for i in eachindex(A)
    @inbounds begin
      if A[i] < 0.0
        A[i] = 0.0
  end end end
  # normalize columns
  for c in 1:size(A,2)
    @inbounds A[:,c] ./= norm(A[:,c])
  end
  return nothing
end
function projectA!(o::NNSC_optimizer)
  project!(o.A)
  return nothing
end

"""
    initialize_mats(n,b,k)
`n` : number of dimensions
`b` : number of basis vectors
`k` : number of samples
"""
function initialize_mats(n::I,b::I,k::I) where I<:Integer
  A = rand(distr_initA,(n,b))
  project!(A)
  B = rand(distr_initB,(b,k))
  return A,B
end
const distr_initA = Exponential(1.0)
const distr_initB = Exponential(1.0)

function initialize_mats(X,d)
  (nfull,ndat) = size(X)
  Xv = X[:]
  A = sqrt.(sample(Xv,(nfull,d))./ d)
  B = sqrt.(sample(Xv,(d,ndat)) ./ d)
  return A,B
end

function updateS!(S::M,A::M,X::M,lambda::Real) where M<:Matrix{<:Real}
  return S .*= (A' * X) ./ (A' * A * S .+ lambda)
end
function updateS!(o::NNSC_optimizer)
  o.S .*= o.A'*o.X
  o.S ./= (o.A'*o.AS .+ o.lambda)
  updateAS!(o)
  return nothing
end
@inline function updateAS!(o::NNSC_optimizer)
  mul!(o.AS,o.A,o.S)
  return nothing
end

function updateA!(A::M,S::M,X::M,lambda::Real,
    stepsize::Real,
    obj_pre::Union{Nothing,<:Real}=nothing) where M<:Matrix{<:Real}
  obj_pre = something(obj_pre, nnsc_objective(X,A,S,lambda))
  Apost = A .- stepsize .* ((A*S-X)*S')
  project!(Apost)
  obj_post = nnsc_objective(X,Apost,S,lambda)
  if (obj_post < obj_pre || stepsize < 1E-9)
    copy!(A,Apost)
    return stepsize
  else
    return updateA!(A,S,X,lambda,0.5*stepsize,obj_pre)
  end
end

function updateA!(o::NNSC_optimizer,stepsize::Real,obj_pre::Real)
  Apre = copy(o.A)
  ASpre = copy(o.AS)
  o.A .-= stepsize .* ((o.AS-o.X)*o.S')
  projectA!(o)
  updateAS!(o)
  obj_post = nnsc_objective(o)
  delta_obj = obj_pre - obj_post
  if (delta_obj > 0.0 || stepsize < 1E-9)
    return (stepsize,obj_post,delta_obj)
  else
    copy!(o.A,Apre)
    copy!(o.AS,ASpre)
    return updateA!(o,0.5*stepsize,obj_pre)
  end
end

function nnsc_optimize(X::Matrix{<:Real},d::Integer,
  lambda::Real ; step_zero::Real=0.05, nstepmax::Integer=10_000,
    opt_threshold::Real=1E-4,verbose::Bool=false)
  o = NNSC_optimizer(X,d,lambda)
  stepsize=step_zero
  obj_temp = nnsc_objective(o)
  delta_obj = Inf
  nsteps=0
  while ( (nsteps < nstepmax) && (delta_obj > opt_threshold) )
    nsteps +=1
    (stepsize,obj_temp,delta_obj) = updateA!(o, stepsize,obj_temp)
    updateS!(o)
    if verbose && (nsteps == nstepmax)
      @info """NNSC completed:
      threshold NOT reached: $opt_threshold
      after : $nsteps (max) iterations
      objective : $obj_temp
      """
    elseif verbose && (delta_obj <= opt_threshold)
      @info """NNSC completed:
      threshold reached: $opt_threshold
      after : $nsteps iterations
      objective : $obj_temp
      """
    end
  end
  o.S[o.S .< 1E-7] .= 0.0
  return o
end

#################
# What follows is manipulaton of image patches and files

# regularization types
abstract type ReguType end
# thou shall not broadcast!
Base.Broadcast.broadcastable(reg::ReguType)=Ref(reg)

struct NoRegu <: ReguType end
struct StandardRegu <: ReguType end

function regularize_image(img,rt::Type{NoRegu})::Matrix{Float64}
    return convert(Array{Float64},Gray.(img))
end
function regularize_image(img,rt::Type{StandardRegu})::Matrix{Float64}
    im=regularize_image(img,NoRegu)
    _mu,_std = mean(im),std(im)
    @. im = (im -_mu)/_std
    return im
end

"""

"""
function read_natural_images(
    dir,
    sz,
    (regu::Type{<:ReguType}) = StandardRegu;
    rotate::Bool = true,
    verbose::Bool = false,
    format::String = ".jpg")::Vector{Matrix{Float64}}
  @assert isdir(dir) "directory $dir not found!"
  filenames = filter(f-> occursin(Regex("$format\\b"),f),  readdir(dir))
  filenames = joinpath.(dir,filenames)
  good_size(img, sz) = min(size(img)...) >= sz
  @assert !isempty(filenames) "no files with extension $format found in $dir"
  ret = Matrix{Any}[]
  _rot_fun = rotate ? rotation_functions : [identity]
  nfiles = length(filenames)
  for (k, file) in enumerate(filenames)
    verbose && println("reading file $k/$(nfiles)")
    img = load(file)
    if good_size(img, sz)
      for rot in _rot_fun
        push!(ret, rot(img))
      end
    else
      @warn "image $file has size smaller than $sz , ignored"
    end
  end
  return cut_patch.(regularize_image.(ret,regu),Ref(sz))
end

# rotate images exactly
function rot_90(img)
    r,c=size(img)
    out=similar(img,c,r)
    for cc in 1:c
        out[cc,:] = img[end:-1:1,cc]
    end
    out
end
function rot_m90(img)
    r,c=size(img)
    out=similar(img,c,r)
    for (cc,mcc) in zip(1:c,c:-1:1)
        out[cc,:] = img[:,mcc]
    end
    out
end
function rot_180(img)
    r,c=size(img)
    out=similar(img)
    for (cc,mcc) in zip(1:c,c:-1:1)
        out[:,cc] = img[r:-1:1,mcc]
    end
    out
end
const rotation_functions = [identity , rot_90, rot_m90, rot_180 ]

"""
        function cut_patch(img, sz)

cuts a random square section of size `sz` from the matrix `img`
"""
function cut_patch(img::AbstractArray{T}, sz) where T
  _rows, _cols = size(img)
  _crows,_ccols =_rows-sz+1, _cols-sz+1
  i,j= rand(1:_crows), rand(1:_ccols)
  return img[i:i+sz-1,j:j+sz-1]
end

"""
    sampling_tiles(n_samples::Integer, images::Vector{Matrix{Float64}},
        patch_size::Integer)
Generates `n_samples` patches of size `patch_size`,
"""
function get_random_patches(n_samples::Integer, images::Vector{Matrix{Float64}},
    patch_size::Integer)
  n_img=length(images)
  @assert !isempty(images) "images are missing!"
  #check how many patches per image
  n_full,n_rem = divrem(n_samples,n_img)
  n_imgsampl = fill(n_full,n_img)
  # add the rem, picking them randomly!
  if n_rem>0
    idx_rem=sample(1:n_img , n_rem ;replace=false)
    n_imgsampl[idx_rem] .+= 1
  end
  @assert  sum(n_imgsampl) == n_samples "wrong selection?"
  ret = Vector{Matrix{Float64}}(undef,n_samples)
  k=0
  for (im,ns) in zip(images,n_imgsampl) , s in 1:ns
    k+=1
    ret[k]=cut_patch(im,patch_size)
  end
  return ret
end

function regularize01!(mat)
    l,u = extrema(mat)
    @. mat = (mat-l)/(u-l)
end

function merge_mats(mats::Vector{Matrix{R}};
   zeromean::Bool=true) where R<:Real
  sz = size(mats[1],1)
  nmats = length(mats)
  ret = Matrix{R}(undef,sz*sz,nmats)
  for (k,mat) in enumerate(mats)
    ret[:,k] .= mat[:]
  end
  if zeromean
    broadcast!(-,ret,ret,mean(ret;dims=2))
  end
  return ret
end

# pretty much wrappers to MultivariateStats

function plus_minus_parts(M::Array{T}) where T<:Real
  plus,minus = zero(M),zero(M)
  for i in eachindex(M)
    @inbounds begin
    if M[i] > 0
      plus[i] = M[i]
    else
      minus[i] = -M[i]
  end end end
  return plus,minus
end

abstract type Preprocess end

struct PreprocessWhiten <: Preprocess end
function preprocess_plus_minus(X,::Type{PreprocessWhiten})
  M = fit(Whitening,X)
  @assert all(x-> isapprox(x,0.0 ; atol=1E-6),mean(M)) "The input should have zero mean!"
  return (M, plus_minus_parts(transform(M,X))...)
end
function revert_plus_minus(Mp,Mm,W::Whitening)
  return  W.W' \ (Mp - Mm)
end

struct PreprocessPCA <: Preprocess
  d::Integer
end

function preprocess_plus_minus(X,p::PreprocessPCA)
  M = fit(PCA,X;maxoutdim=p.d)
  @assert all(x-> isapprox(x,0.0 ; atol=1E-6),mean(M)) "The input should have zero mean!"
  return (M, plus_minus_parts(transform(M,X))...)
end
function revert_plus_minus(Mp,Mm,p::PCA)
  return reconstruct(p, Mp - Mm)
end

struct PreprocessNot <: Preprocess end
function preprocess_plus_minus(X,::Type{PreprocessNot})
  @assert all(x-> isapprox(x,0.0 ; atol=1E-6),mean(X;dims=2)) "The input should have zero mean!"
  return (PreprocessNot, plus_minus_parts(X)...)
end
function revert_plus_minus(Mp,Mm,::Type{PreprocessNot})
  return (Mp - Mm)
end

end #module


text:
# When an object like DimArray or DimStack is constructed,
# `format` is called on its dimensions check that the axes match,
# and fill in any `Auto-` fields.

"""
    format(dims, x) => Tuple{Vararg{<:Dimension,N}}

Format the passed-in dimension(s) `dims` to match the object `x`.

This means converting indexes of `Tuple` to `LinRange`, and running
`format`. Errors are also thrown if dims don't match the array dims or size.

If a [`LookupArray`](@ref) hasn't been specified, a lookup is chosen
based on the type and element type of the index.
"""
format(dims, A::AbstractArray) = format((dims,), A)
function format(dims::NamedTuple, A::AbstractArray)
    dims = map(keys(dims), values(dims)) do k, v
        rebuild(key2dim(k), v)
    end
    return format(dims, axes(A))
end
format(dims::Tuple{Vararg{<:Any,N}}, A::AbstractArray{<:Any,N}) where N =
    format(dims, axes(A))
@noinline format(dims::Tuple{Vararg{<:Any,M}}, A::AbstractArray{<:Any,N}) where {N,M} =
    throw(DimensionMismatch("Array A has $N axes, while the number of dims is $M: $(map(basetypeof, dims))"))
format(dims::Tuple, axes::Tuple) = map(_format, dims, axes)

_format(dimname::Symbol, axis::AbstractRange) = Dim{dimname}(NoLookup(axes(axis, 1)))
_format(::Type{D}, axis::AbstractRange) where D<:Dimension = D(NoLookup(axes(axis, 1)))
_format(dim::Dimension{Colon}, axis::AbstractRange) = rebuild(dim, NoLookup(axes(axis, 1)))
function _format(dim::Dimension, axis::AbstractRange)
    newlookup = format(val(dim), basetypeof(dim), axes(axis, 1))
    checkaxis(newlookup, axis)
    return rebuild(dim, newlookup)
end

format(val::AbstractArray, D::Type, axis::AbstractRange) = format(AutoLookup(), D, val, axis)
format(m::LookupArray, D::Type, axis::AbstractRange) = format(m, D, parent(m), axis)

# Format LookupArrays
# No more identification required for NoLookup
format(m::NoLookup, D::Type, index, axis::AbstractRange) = m
format(m::NoLookup, D::Type, index::AutoIndex, axis::AbstractRange) = NoLookup(axis)
# # AutoLookup
function format(m::AutoLookup, D::Type, index::AbstractArray{T}, axis::AbstractRange) where T
    # A mixed type index is Categorical
    m = if isconcretetype(T)
        Sampled(; order=order(m), span=span(m), sampling=sampling(m), metadata=metadata(m))
    else
        o = order(m) isa AutoOrder ? Unordered() : order(m)
        Categorical(; order=o, metadata=metadata(m))
    end
    format(m, D, index, axis)
end
function format(m::AutoLookup, D::Type, index::AbstractArray{<:CategoricalEltypes}, axis::AbstractRange)
    o = _format(order(m), D, index)
    return Categorical(index; order=o, metadata=metadata(m))
end
function format(m::Categorical, D::Type, index, axis::AbstractRange)
    i = _format(index, axis)
    o = _format(order(m), D, index)
    return rebuild(m; data=i, order=o)
end
# # Sampled
function format(m::AbstractSampled, D::Type, index, axis::AbstractRange)
    i = _format(index, axis)
    o = _format(order(m), D, index)
    sp = _format(span(m), D, index)
    sa = _format(sampling(m), sp, D, index)
    x = rebuild(m; data=i, order=o, span=sp, sampling=sa)
    return x
end
# # Transformed
format(m::Transformed, D::Type, index::AutoIndex, axis::AbstractRange) = rebuild(m; data=axis)
format(m::Transformed, D::Type, index, axis::AbstractRange) = m

# Index
_format(index::AbstractArray, axis::AbstractRange) = index
_format(index::AutoLookup, axis::AbstractRange) = axis
# Order
_format(order::Order, D::Type, index) = order
_format(order::AutoOrder, D::Type, index) = LookupArrays.orderof(index)
# Span
_format(span::AutoSpan, D::Type, index::Union{AbstractArray,Val}) =
    _format(Irregular(), D, index)
_format(span::AutoSpan, D::Type, index::AbstractRange) = Regular(step(index))
_format(span::Regular{AutoStep}, D::Type, index::Union{AbstractArray,Val}) = _arraynosteperror()
_format(span::Regular{AutoStep}, D::Type, index::AbstractRange) = Regular(step(index))
_format(span::Regular, D::Type, index::Union{AbstractArray,Val}) = span
function _format(span::Regular, D::Type, index::AbstractRange)
    step(span) isa Number && !(step(span) ≈ step(index)) && _steperror(index, span)
    return span
end
_format(span::Irregular{AutoBounds}, D, index) = Irregular(nothing, nothing)
_format(span::Irregular{<:Tuple}, D, index) = span
_format(span::Explicit, D, index) = span
# Sampling
_format(sampling::AutoSampling, span::Span, D::Type, index) = Points()
_format(sampling::AutoSampling, span::Explicit, D::Type, index) =
    Intervals(_format(locus(sampling), D, index))
_format(sampling::Points, span::Span, D::Type, index) = sampling
_format(sampling::Points, span::Explicit, D::Type, index) = _explicitpoints_error()
_format(sampling::Intervals, span::Span, D::Type, index) =
    rebuild(sampling, _format(locus(sampling), D, index))
# Locus
_format(locus::AutoLocus, D::Type, index) = Center()
# Time dimensions need to default to the Start() locus, as that is
# nearly always the _format and Center intervals are difficult to
# calculate with DateTime step values.
_format(locus::AutoLocus, D::Type{<:TimeDim}, index) = Start()
_format(locus::Locus, D::Type, index) = locus

_order(index) = first(index) <= last(index) ? ForwardOrdered() : ReverseOrdered()

@noinline _explicitpoints_error() =
    throw(ArgumentError("Cannot use Explicit span with Points sampling"))
@noinline _steperror(index, span) =
    throw(ArgumentError("lookup step $(step(span)) does not match index step $(step(index))"))
@noinline _arraynosteperror() =
    throw(ArgumentError("`Regular` must specify `step` size with an index other than `AbstractRange`"))

checkaxis(lookup::Transformed, axis) = nothing
function checkaxis(lookup, axis)
    if !(first(axes(lookup)) == axis)
        throw(DimensionMismatch(
            "axes of $(basetypeof(lookup)) of $(first(axes(lookup))) do not match array axis of $axis"
        ))
    end
end


text:
using HypergeometricFunctions
using SpecialFunctions

"""Kravchuk.jl contains the functions to evolve a set of N + 1 amplitude points
over the discrete and finite Kravchuk oscillator realization"""

# Kravchuk symmetric polynomial definition. This particular functions is using
# the resource given in HypergeometricFunctions.jl, the method 2, who avoids
# floating-points problems, but has a flaw for the definitions in the corners
# when n = 2j and m = j (energy and position respectively).
kravd(n, m, j) = HypergeometricFunctions._₂F₁general2(-n, -m, -2j, 2)

# Normalizing terms and binomial re-definition
A(n, j) = (-1)^n / 2^j
binom(a, j) = binomial(2j, a)

# Primitive Kravchuk symmetric function. Needs a fix for the flaw described above
kravp(n, m, j) = A(n, j) * sqrt(binom(n, j) * binom(m + j, j)) * kravd(n, m + j, j)

# Conditional ad-hoc solution for the problem in the _2F1 function
function krav(n, m, j)
    if n == 2j && m == j
        return kravp(0, -j, j)
    else
        return kravp(n, m, j)
    end
end

# Exponential counting the rotation angles of the fractional evolution
kerk(n, α) = exp(-1im * pi * α * n / 2)

# Evolution kernel, this is intended to be used as standalone
frk(n, m, α, j) = sum([krav(l, n, j) * kerk(l, α) * krav(l, m, j) for l in 0:2j])

# Evolution of initial amplitudes in the vector v. This functions is intended
# to be used when the evolution is required, just given v of size 2j + 1, and
# the numbers of steps p
function evol(v, p, j)
    θ = range(0, stop = 2, length = p)
    em = zeros(Float64, (2j + 1, p))
    for r in 1:p
        for s in -j:j
        em[s + 1 + j, r] = abs2(sum([frk(s, m - j - 1, θ[r], j) * v[m] for m in 1:(2j +1)]))
        end
    end
    return em
end


text:
"""
    call_tetgen(setup, surfaces)

Call Tetgen on surface geometry.
"""
function call_tetgen(surfaces, refinement = nothing)
    @unpack points, facets, facetmarkers, regions = surfaces

    nregion = size(regions, 2)

    input = RawTetGenIO{Cdouble}()
    input.pointlist = points

    facetlist!(input, facets)

    input.facetmarkerlist = facetmarkers

    input.regionlist = [
        regions
        collect(1:nregion)'
        fill(0.1, 1, nregion)
    ]

    if isnothing(refinement)
        tetgen = tetrahedralize(input, "pqA")
    else
        tetgen = tetrahedralize(input, "pqAa$(refinement)")
    end

    points = tetgen.pointlist
    facets = Int.(tetgen.trifacelist)
    facetmarkers = Int.(tetgen.trifacemarkerlist)
    elements = Int.(tetgen.tetrahedronlist)
    elementmarkers = Int.(tetgen.tetrahedronattributelist[:])

    (; points, facets, facetmarkers, elements, elementmarkers)
end


text:
"""
`lineplot(x, y; nargs...)` → `Plot`

Description
============

Draws a path through the given points on a new canvas.
It uses the first parameter `x`
(which should be a vector or range) to denote
the horizontal position of each point,
and the second parameter `y`
(which should also be a vector or range)
as their vertical position.
This means that the two vectors have to have the same length.

Usage
======

    lineplot(x, y; title = "", name = "", width = 40, height = 15, border = :solid, xlim = [0, 0], ylim = [0, 0], margin = 3, padding = 1, color = :auto, labels = true, canvas = BrailleCanvas, grid = true)

    lineplot(fun[, start, end]; title = "", name = "", width = 40, height = 15, border = :solid, xlim = [0, 0], ylim = [0, 0], margin = 3, padding = 1, color = :auto, labels = true, canvas = BrailleCanvas, grid = true)

Arguments
==========

- **`x`** : The horizontal dimension for each point. Can be a Real number or of type `TimeType`

- **`y`** : The vertical dimension for each point.

- **`fun`** : A function f: R -> R that should be drawn from `start` to `end`

- **`title`** : Text to display on the top of the plot.

- **`name`** : Annotation of the current drawing to displayed on the right

- **`width`** : Number of characters per row that should be used for plotting.

- **`height`** : Number of rows that should be used for plotting. Not applicable to `barplot`.

- **`border`** : The style of the bounding box of the plot. Supports `:solid`, `:bold`, `:dashed`, `:dotted`, `:ascii`, and `:none`.

- **`xlim`** : Plotting range for the x coordinate. `[0, 0]` stands for automatic.

- **`ylim`** : Plotting range for the y coordinate. `[0, 0]` stands for automatic.

- **`margin`** : Number of empty characters to the left of the whole plot.

- **`padding`** : Space of the left and right of the plot between the labels and the canvas.

- **`color`** : Color of the drawing. Can be any of `:blue`, `:red`, `:yellow`

- **`labels`** : Can be used to hide the labels by setting `labels=false`.

- **`canvas`** : The type of canvas that should be used for drawing.

- **`grid`** : Can be used to hide the gridlines at the origin

Returns
========

A plot object of type `Plot{T<:Canvas}`

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)
- Milktrader (Github: https://github.com/milktrader)

Examples
=========

    julia> lineplot([1, 2, 7], [9, -6, 8], title = "My Lineplot")

                        My Lineplot
         ┌────────────────────────────────────────┐
       9 │⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠⠊│
         │⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠁⠀⠀│
         │⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠊⠁⠀⠀⠀⠀│
         │⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠊⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠒⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠉⠉⠉⠉⡏⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⢉⠭⠋⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠁│
         │⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⡠⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⢀⠔⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠀⠀⠀⠸⡀⠀⠀⡠⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      -6 │⠀⠀⠀⠀⠀⠀⢇⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         └────────────────────────────────────────┘
         1                                        7

See also
=========

`Plot`, `scatter`, `stairs`, `BrailleCanvas`, `BlockCanvas`, `AsciiCanvas`
"""
function lineplot(
        x::AbstractVector{<:Number},
        y::AbstractVector{<:Number};
        color::Symbol = :auto,
        name::AbstractString = "",
        canvas::Type = BrailleCanvas,
        kw...)
    X = convert(Vector{Float64},x)
    Y = convert(Vector{Float64},y)
    new_plot = Plot(X, Y, canvas; kw...)
    color = color == :auto ? next_color!(new_plot) : color
    name == "" || annotate!(new_plot, :r, name, color)
    lines!(new_plot, X, Y, color)
end

function lineplot!(
        plot::Plot{<:Canvas},
        x::AbstractVector{<:Number},
        y::AbstractVector{<:Number};
        color::Symbol = :auto,
        name::AbstractString = "",
        kw...)
    X = convert(Vector{Float64},x)
    Y = convert(Vector{Float64},y)
    color = color == :auto ? next_color!(plot) : color
    name == "" || annotate!(plot, :r, name, color)
    lines!(plot, X, Y, color)
end

function lineplot!(
        plot::Plot{<:Canvas},
        intercept::Number,
        slope::Number;
        kw...)
    xmin = origin_x(plot.graphics)
    xmax = origin_x(plot.graphics) + width(plot.graphics)
    ymin = origin_y(plot.graphics)
    ymax = origin_y(plot.graphics) + height(plot.graphics)
    lineplot!(plot, [xmin, xmax], [intercept + xmin*slope, intercept + xmax*slope]; kw...)
end

function lineplot!(
        plot::Plot{<:Canvas},
        Y::Function,
        x::AbstractRange;
        kw...)
    X = collect(x)
    lineplot!(plot, Y, X; kw...)
end

function lineplot(
        Y::Function,
        x::AbstractRange;
        kw...)
    X = collect(x)
    lineplot(Y, X; kw...)
end

function lineplot!(
        plot::Plot{<:Canvas},
        Y::Function,
        X::AbstractVector{<:Number};
        name::AbstractString = "",
        kw...)
    y = convert(Vector{Float64}, [Y(i) for i in X])
    name = name == "" ? string(Y, "(x)") : name
    lineplot!(plot, X, y; name = name, kw...)
end

function lineplot(
        Y::Function,
        X::AbstractVector{<:Number};
        name::AbstractString = "",
        kw...)
    y = convert(Vector{Float64}, [Y(i) for i in X])
    name = name == "" ? string(Y, "(x)") : name
    new_plot = lineplot(X, y; name = name, kw...)
    xlabel!(new_plot, "x")
    ylabel!(new_plot, "f(x)")
end

function lineplot!(
        plot::Plot{<:Canvas},
        Y::Function,
        startx::Number,
        endx::Number;
        kw...)
    diff = abs(endx - startx)
    X = collect(startx:(diff/(3*ncols(plot.graphics))):endx)
    lineplot!(plot, Y, X; kw...)
end

function lineplot(Y::Function; kw...)
    lineplot(Y, -10, 10; kw...)
end

function lineplot!(
        plot::Plot{<:Canvas},
        Y::Function;
        kw...)
    lineplot!(plot, Y, origin_x(plot.graphics), origin_x(plot.graphics) + width(plot.graphics); kw...)
end

function lineplot(
        Y::Function,
        startx::Number,
        endx::Number;
        width::Int = 40,
        kw...)
    diff = abs(endx - startx)
    X = collect(startx:(diff/(3*width)):endx)
    lineplot(Y, X; width = width, kw...)
end

function lineplot(Y::AbstractVector{<:Function}; kw...)
    lineplot(Y, -10, 10; kw...)
end

function lineplot(
        Y::AbstractVector{<:Function},
        startx::Number,
        endx::Number;
        kw...)
    n = length(Y)
    @assert n > 0
    new_plot = lineplot(Y[1], startx, endx; kw...)
    for i = 2:n
        lineplot!(new_plot, Y[i], startx, endx; kw...)
    end
    new_plot
end

function lineplot(
        X::AbstractRange{<:Number},
        Y::AbstractRange{<:Number};
        kw...)
    lineplot(collect(X), collect(Y); kw...)
end

function lineplot(
        X::AbstractRange,
        Y::AbstractVector{<:Number};
        kw...)
    lineplot(collect(X), Y; kw...)
end

function lineplot(
        X::AbstractVector{<:Number},
        Y::AbstractRange;
        kw...)
    lineplot(X, collect(Y); kw...)
end

function lineplot(X::AbstractVector; kw...)
    lineplot(1:length(X), X; kw...)
end

function lineplot!(
        plot::Plot{<:Canvas},
        X::AbstractVector;
        kw...)
    lineplot!(plot, 1:length(X), X; kw...)
end

function lineplot!(
        plot::Plot{<:Canvas},
        X::AbstractRange{<:Number},
        Y::AbstractRange{<:Number};
        kw...)
    lineplot!(plot, collect(X), collect(Y); kw...)
end

function lineplot!(
        plot::Plot{<:Canvas},
        X::AbstractRange,
        Y::AbstractVector{<:Number};
        kw...)
    lineplot!(plot, collect(X), Y; kw...)
end

function lineplot!(
        plot::Plot{<:Canvas},
        X::AbstractVector{<:Number},
        Y::AbstractRange;
        kw...)
    lineplot!(plot, X, collect(Y); kw...)
end

function lineplot(
        X::AbstractVector{<:TimeType},
        Y::AbstractVector{<:Number};
        kw...)
    d = convert(Vector{Float64}, Dates.value.(X))
    new_plot = lineplot(d, Y; kw...)
    annotate!(new_plot, :bl, string(first(X)))
    annotate!(new_plot, :br, string(last(X)))
end


text:
export main

using ArgParse
using Printf

# julia -e "Pkg" ...

function main(args)

    s = ArgParseSettings("Calculate the CRC for the given files")

    @add_arg_table s begin
        "--list", "-l"
        help = "list available CRC algorithms"
        action = :store_true
        "--decimal", "-d"
        help = "show checksums as decimal values (default is hex)"
        action = :store_true
        "--crc", "-c"
        help = "name the CRC to use"
        default = "CRC_32"
        "--append", "-a"
        help = "combine the data from all files"
        action = :store_true
        "files"
        help = "the files to read (- for stdin)"
        nargs = '*'
    end

    parsed_args = parse_args(args, s)

    if parsed_args["list"]
        names = sort(collect(keys(ALL)), by=n->ALL[n])
        for name in names
            @printf("%s %s\n", name, ALL[name])
        end
    end

    name = Symbol(parsed_args["crc"])
    if !haskey(ALL, name)
        error("CRC $name is not defined")
    end

    append = parsed_args["append"]
    spec = ALL[name]
    c = crc(spec)
    if parsed_args["decimal"]
        fmt = x -> Int(x)
    else
        fmt = x -> "0x" * string(x, base=16, pad = 1 + div(spec.width-1, 4))
    end
    sum = 0
    for file in parsed_args["files"]
        if file == "-"
            sum = c(STDIN, append=append)
        else
            open(file, "r") do f
                sum = c(f, append=append)
            end
        end
        if !append
            println("$(fmt(sum)) $file")
        end
    end
    if append
        println("$(fmt(sum))")
    end

end


text:
using QUBOConstraints
using Documenter

DocMeta.setdocmeta!(QUBOConstraints, :DocTestSetup, :(using QUBOConstraints); recursive=true)

makedocs(;
    modules=[QUBOConstraints],
    authors="azzaare <jf@baffier.fr> and contributors",
    repo="https://github.com/JuliaConstraints/QUBOConstraints.jl/blob/{commit}{path}#{line}",
    sitename="QUBOConstraints.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaConstraints.github.io/QUBOConstraints.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaConstraints/QUBOConstraints.jl",
    devbranch="main",
)


text:
let
    roadway = get_test_roadway()
    trajdata = get_test_trajdata(roadway)

    Δt = 0.1
    rec = SceneRecord(5, Δt)
    @test capacity(rec) == 5
    @test nframes(rec) == 0
    @test !pastframe_inbounds(rec, 0)
    @test !pastframe_inbounds(rec, -1)
    @test !pastframe_inbounds(rec, 1)

    scene = get!(Scene(), trajdata, 1)
    update!(rec, scene)
    @test nframes(rec) == 1
    @test pastframe_inbounds(rec, 0)
    @test !pastframe_inbounds(rec, -1)
    @test !pastframe_inbounds(rec, 1)
    @test isapprox(get_elapsed_time(rec, 0), 0)
    @test rec[0][1].state == get_state(trajdata, 1, 1)
    @test rec[0][1].def == get_def(trajdata, 1)
    @test rec[0][2].state == get_state(trajdata, 2, 1)
    @test rec[0][2].def == get_def(trajdata, 2)
    show(IOBuffer(), rec)

    get!(scene, trajdata, 2)
    update!(rec, scene)
    @test nframes(rec) == 2
    @test pastframe_inbounds(rec, 0)
    @test pastframe_inbounds(rec, -1)
    @test !pastframe_inbounds(rec, 1)
    @test isapprox(get_elapsed_time(rec,  0),  0)
    @test isapprox(get_elapsed_time(rec, -1), Δt)
    @test isapprox(get_elapsed_time(rec, -1, 0), Δt)
    @test rec[0][1].state == get_state(trajdata, 1, 2)
    @test rec[0][1].def == get_def(trajdata, 1)
    @test rec[0][2].state == get_state(trajdata, 2, 2)
    @test rec[0][2].def == get_def(trajdata, 2)
    @test rec[-1][1].state == get_state(trajdata, 1, 1)
    @test rec[-1][1].def == get_def(trajdata, 1)
    @test rec[-1][2].state == get_state(trajdata, 2, 1)
    @test rec[-1][2].def == get_def(trajdata, 2)

    scene2 = get!(Scene(), rec)
    @test scene2[1].state == get_state(trajdata, 1, 2)
    @test scene2[1].def == get_def(trajdata, 1)
    @test scene2[2].state == get_state(trajdata, 2, 2)
    @test scene2[2].def == get_def(trajdata, 2)

    get!(scene2, rec, -1)
    @test scene2[1].state == get_state(trajdata, 1, 1)
    @test scene2[1].def == get_def(trajdata, 1)
    @test scene2[2].state == get_state(trajdata, 2, 1)
    @test scene2[2].def == get_def(trajdata, 2)

    empty!(rec)
    @test nframes(rec) == 0

    test_veh_state = VehicleState(VecSE2(7.0,7.0,2.0), roadway, 10.0)
    test_veh_def = VehicleDef(AgentClass.CAR, 5.0, 3.0)
    test_veh = Vehicle(test_veh_state, test_veh_def, 999)
    rec[-1][1] = test_veh
    @test rec[-1][1].state == test_veh_state
end


text:
"""
	Rocket
"""

abstract type Rocket{I, T} <: Model{I, T} end

struct RocketNominal{I, T} <: Rocket{I, T}
	n::Int
	m::Int
	d::Int

	g  # gravity
	m1 # mass
	l1 # length
	J  # inertia
end

struct RocketSlosh{I, T} <: Rocket{I, T}
	n::Int
	m::Int
	d::Int

	g  # gravity
	m1 # mass
	l1 # length
	J  # inertia

	m2 # mass of fuel
	l2 # length of fuel pendulum
	l3 # length from rocket COM to fuel pendulum
end

g = 9.81 # gravity
m1 = 1.0 # mass of rocket
l1 = 0.5 # length from COM to thruster
J = 1.0 / 12.0 * m1 * (2.0 * l1)^2.0 # inertia of rocket

m2 = 0.1 # mass of pendulum
l2 = 0.1 # length from COM to pendulum
l3 = 0.1 # length of pendulum

# models
m = 2

# nominal model
n_nominal = 6
model_nominal = RocketNominal{Midpoint, FixedTime}(n_nominal, m, n_nominal, g, m1, l1, J)

# slosh model
n_slosh = 8
model_slosh = RocketSlosh{Midpoint, FixedTime}(n_slosh, m, n_slosh, g, m1 - m2, l1, J, m2, l2, l3)

function k_thruster(model::Rocket, x)
	y, z, θ = x[1:3]
	py = y + model.l1 * sin(θ)
	pz = z - model.l1 * cos(θ)

	return [py;
		    pz]
end

function jacobian_thruster(model::RocketNominal, x)
	y, z, θ = x[1:3]

	return [1.0 0.0 model.l1 * cos(θ);
			0.0 1.0 model.l1 * sin(θ)]
end

function jacobian_thruster(model::RocketSlosh, x)
	y, z, θ = x[1:3]

	return [1.0 0.0 model.l1 * cos(θ) 0.0;
			0.0 1.0 model.l1 * sin(θ) 0.0]
end

function k_mass(model::RocketSlosh, x)
	yp = x[1] + l2 * sin(x[3]) + l3 * sin(x[4])
	zp = x[2] - l2 * cos(x[3]) - l3 * cos(x[4])

	return [yp;
			zp]
end

function jacobian_mass(model::RocketSlosh, x)
	return [1.0 0.0 l2 * cos(x[3]) l3 * cos(x[4]);
	        0.0 1.0 l2 * sin(x[3]) l3 * sin(x[4])]
end

function lagrangian(model::RocketNominal, q, q̇)
	return (0.5 * model.m1 * (q̇[1]^2 + q̇[2]^2.0) + 0.5 * model.J * q̇[3]^2.0
			- model.m1 * model.g * q[2])
end

function lagrangian(model::RocketSlosh, q, q̇)
	zp = k_mass(model, q)[2]
	vp = jacobian_mass(model, q) * q̇

	return (0.5 * model.m1 * (q̇[1]^2 + q̇[2]^2) + 0.5 * model.J * q̇[3]^2.0
			- model.m1 * model.g * q[2]
			+ 0.5 * model.m2 * vp' * vp
			- model.m2 * model.g * zp)
end

function dLdq(model::Rocket, q, q̇)
	Lq(x) = lagrangian(model, x, q̇)
	ForwardDiff.gradient(Lq, q)
end

function dLdq̇(model::Rocket, q, q̇)
	Lq̇(x) = lagrangian(model, q, x)
	ForwardDiff.gradient(Lq̇, q̇)
end

function f(model::Rocket, x, u, w)
	nq = convert(Int, floor(model.n / 2.0))
	q = x[1:nq]
	q̇ = x[nq .+ (1:nq)]
	tmp_q(z) = dLdq̇(model, z, q̇)
	tmp_q̇(z) = dLdq̇(model, q, z)
	[q̇;
	 ForwardDiff.jacobian(tmp_q̇,q̇) \ (-1.0 * ForwardDiff.jacobian(tmp_q, q) * q̇
	 	+ dLdq(model, q, q̇)
		+ jacobian_thruster(model, q)' * u[1:2])]
end

function state_output(model::RocketSlosh, x)
	x[collect([1, 2, 3, 5, 6, 7])]
end

function visualize_rocket!(vis, model::Rocket, x, x1, xT;
       Δt = 0.1, version = :slosh_nom, T_off = length(x))

	T = length(x)

	default_background!(vis)

	body = Cylinder(Point3f0(0.0, 0.0, -1.0 * model.l1),
		Point3f0(0.0, 0.0, 3.0 * model.l1),
		convert(Float32, 0.225))

	setobject!(vis["rocket"]["body"], body,
		MeshPhongMaterial(color = RGBA(1.0, 0.0, 0.0, 1.0)))

	obj_rocket = joinpath(pwd(), "models/rocket/space_x_booster.obj")
	mtl_rocket = joinpath(pwd(), "models/rocket/space_x_booster.mtl")
	ctm = ModifiedMeshFileObject(obj_rocket,mtl_rocket,scale=1.0)
	setobject!(vis["rocket"]["falcon9"], ctm)

	settransform!(vis["rocket"]["falcon9"],
		compose(Translation(([-xT[1] + 4.0; -4.0; xT[2] - 1.0])...),
		LinearMap(RotY(1.0 * xT[3]) * RotZ(pi) * RotX(pi / 2.0))))

	settransform!(vis["rocket"],
		compose(Translation(-1.0 * x1[1], 0.0, x1[2]),
		LinearMap(RotY(x1[3]))))

	setvisible!(vis["rocket"]["falcon9"], true)
	obj_platform = joinpath(pwd(), "models/rocket/space_x_platform.obj")
	mtl_platform = joinpath(pwd(), "models/rocket/space_x_platform.mtl")

	obj_platform = joinpath(pwd(), "models/rocket/space_x_platform.obj")
	mtl_platform = joinpath(pwd(), "models/rocket/space_x_platform.mtl")
	ctm_platform = ModifiedMeshFileObject(obj_platform,mtl_platform,scale=1.0)
	setobject!(vis["platform"],ctm_platform)
	settransform!(vis["platform"], compose(Translation(0.0,2.5,-0.85),LinearMap(RotZ(pi)*RotX(pi/2))))

	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    for t = 1:T
        MeshCat.atframe(anim,t) do
			if t >= T_off
				setvisible!(vis["rocket"]["body"], false)
			else
				setvisible!(vis["rocket"]["body"], true)
			end
			settransform!(vis["rocket"],
				compose(Translation(-1.0 * x[t][1], 0.0, x[t][2]),
				LinearMap(RotY(x[t][3]))))
        end
    end

	settransform!(vis["/Cameras/default"], compose(Translation(0.0, 20.0, -1.0),
		LinearMap(RotZ(pi / 2.0))))
	setvisible!(vis["/Grid"], false)

    MeshCat.setanimation!(vis, anim)
end

function visualize_rocket_ghost!(vis, model::Rocket, x)

	default_background!(vis)

	body = Cylinder(Point3f0(0.0, 0.0, -1.0 * model.l1),
		Point3f0(0.0, 0.0, 3.0 * model.l1),
		convert(Float32, 0.225))
	rkt_offset = [3.9,-6.35,0.2]
	y_shift = 2.35
	obj_rocket = joinpath(pwd(), "models/rocket/space_x_booster.obj")
	mtl_rocket = joinpath(pwd(), "models/rocket/space_x_booster.mtl")
	ctm = ModifiedMeshFileObject(obj_rocket,mtl_rocket,scale=1.0)

	for (t, _x) in enumerate(x)
		println("t= $t")
		setobject!(vis["rocket_$t"]["body"], body,
			MeshPhongMaterial(color = RGBA(1.0, 0.0, 0.0, 1.0)))
		setvisible!(vis["rocket_$t"]["body"], t == length(x) ? false : true)

		setobject!(vis["rocket_$t"]["falcon9"], ctm)
		settransform!(vis["rocket_$t"]["falcon9"],
			compose(Translation(([-x[end][1]- 0.1225; y_shift; x[end][2]-1.15] + rkt_offset)...),
			LinearMap(RotY(1.0 * x[end][3]) * RotZ(pi) * RotX(pi / 2.0))))
		setvisible!(vis["rocket_$t"]["falcon9"], true)

		settransform!(vis["rocket_$t"],
			compose(Translation(-1.0 * _x[1], 0.0, _x[2]),
			LinearMap(RotY(_x[3]))))
	end

	obj_platform = joinpath(pwd(), "models/rocket/space_x_platform.obj")
	mtl_platform = joinpath(pwd(), "models/rocket/space_x_platform.mtl")

	obj_platform = joinpath(pwd(), "models/rocket/space_x_platform.obj")
	mtl_platform = joinpath(pwd(), "models/rocket/space_x_platform.mtl")
	ctm_platform = ModifiedMeshFileObject(obj_platform,mtl_platform,scale=1.0)

	setobject!(vis["platform"],ctm_platform)
	settransform!(vis["platform"], compose(Translation(0.0,2.5,-0.85),LinearMap(RotZ(pi)*RotX(pi/2))))

	setvisible!(vis["/Grid"], false)

	settransform!(vis["/Cameras/default"], compose(Translation(0.0, 20.0, -1.0),
		LinearMap(RotZ(pi / 2.0))))
end

function visualize_simple!(vis, model::Rocket, x;
       Δt = 0.1)

	default_background!(vis)

	l1 = Cylinder(Point3f0(0.0, 0.0, -1.0 * model.l1),
		Point3f0(0.0, 0.0, 3.0 * model.l1),
		convert(Float32, 0.25))

    setobject!(vis["rocket"], l1,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

   	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    for t = 1:length(x)
        MeshCat.atframe(anim,t) do
			settransform!(vis["rocket"],
				compose(Translation(-1.0 * x[t][1], 0.0, x[t][2]),
				LinearMap(RotY(x[t][3]))))
        end
    end

    settransform!(vis["/Cameras/default"], compose(Translation(0.0, 25.0, -1.0),
		LinearMap(RotZ(pi / 2.0))))

    MeshCat.setanimation!(vis, anim)
end


text:
using Documenter, Revise

makedocs(
    modules = [Revise],
    clean = false,
    format = :html,
    sitename = "Revise.jl",
    authors = "Tim Holy",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "config.md",
        "limitations.md",
        "internals.md",
        "user_reference.md",
        "dev_reference.md",
    ],
    # # Use clean URLs, unless built as a "local" build
    # html_prettyurls = !("local" in ARGS),
#    html_canonical = "https://juliadocs.github.io/Revise.jl/stable/",
)

deploydocs(
    repo = "github.com/timholy/Revise.jl.git",
    target = "build",
    julia = "nightly",
    deps = nothing,
    make = nothing,
)


text:
# Super-hacked out constructor: DataVector[1, 2, NA]
function ref(::Type{DataVector}, vals...)
    # Get the most generic non-NA type
    toptype = _dv_most_generic_type(vals)

    # Allocate an empty DataVector
    lenvals = length(vals)
    res = DataArray(Array(toptype, lenvals), BitArray(lenvals))

    # Copy from vals into data and mask
    for i = 1:lenvals
        if isna(vals[i])
            res.data[i] = baseval(toptype)
            res.na[i] = true
        else
            res.data[i] = vals[i]
            res.na[i] = false
        end
    end

    return res
end

##############################################################################
##
## isna()
##
##############################################################################

isna(x::Any) = false

##############################################################################
##
## find()
##
##############################################################################

function find(dv::AbstractDataVector{Bool})
    n = length(dv)
    res = Array(Int, n)
    bound = 0
    for i in 1:length(dv)
        if !isna(dv[i]) && dv[i]
            bound += 1
            res[bound] = i
        end
    end
    return res[1:bound]
end

##############################################################################
##
## String representations and printing
##
## TODO: Inherit these from AbstractArray after implementing DataArray
##
##############################################################################

head(dv::AbstractDataVector) = dv[1:min(6, length(dv))]
tail(dv::AbstractDataVector) = dv[max(length(dv) - 6, 1):length(dv)]

##############################################################################
##
## Container operations
##
##############################################################################

# TODO: Fill in definitions for PooledDataVector's
# TODO: Macroize these definitions

function push!{T}(dv::DataVector{T}, v::NAtype)
    push!(dv.data, baseval(T))
    push!(dv.na, true)
    return v
end

function push!{S, T}(dv::DataVector{S}, v::T)
    push!(dv.data, v)
    push!(dv.na, false)
    return v
end

function pop!(dv::DataVector)
    d, m = pop!(dv.data), pop!(dv.na)
    if m
        return NA
    else
        return d
    end
end

function unshift!{T}(dv::DataVector{T}, v::NAtype)
    unshift!(dv.data, baseval(T))
    unshift!(dv.na, true)
    return v
end

function unshift!{S, T}(dv::DataVector{S}, v::T)
    unshift!(dv.data, v)
    unshift!(dv.na, false)
    return v
end

function shift!{T}(dv::DataVector{T})
    d, m = shift!(dv.data), shift(dv.na)
    if m
        return NA
    else
        return d
    end
end

function map(f::Function, dv::DataVector)
    n = length(dv)
    res = DataArray(Any, n)
    for i in 1:n
        res[i] = f(dv[i])
    end
    return res
end


text:
function run_BMA(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2)
    df = DataFrame()
    name = "BMA"
    df_file = "../../../../data/results/baseline/level$level-h$h-J$J/$name.csv"
    ## Forecast using Bayesian Model Averaging
    θ_vec = 0.05:0.01:0.99
    pre_conv = 1_000
    post_conv = 2_000
    θ, ŷ, RMSE, RMSE_cv = bma_oos_forc(J, h, h_lag, J_cv, Y, Y_lag, X_lag, θ_vec, pre_conv, post_conv)
    println("Prior parameter θ chosen for $(name): $θ")
    println("RMSE from Bayesian Model Averaging, h = $h: \nRMSE = $RMSE")
    println("RMSE_cv is $RMSE_cv")
    df[Symbol(name)] = ŷ
    CSV.write(df_file, df)
end

function run_BMA2(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2)
    df = DataFrame()
    name = "BMA2"
    df_file = "../../../../data/results/baseline/level$level-h$h-J$J/$name.csv"
    ## Forecast using Bayesian Model Averaging with interaction terms
    θ_vec = 0.05:0.01:0.99
    pre_conv = 1_000
    post_conv = 2_000
    X_lag = hcat(X_lag, X_lag2)
    θ, ŷ, RMSE, RMSE_cv = bma_oos_forc(J, h, h_lag, J_cv, Y, Y_lag, X_lag, θ_vec, pre_conv, post_conv)
    println("Prior parameter θ chosen for $(name) with interactions: $θ")
    println("RMSE from Bayesian Model Averaging with interaction terms, h = $h: \nRMSE = $RMSE")
    println("RMSE_cv is $RMSE_cv")
    df[Symbol(name)] = ŷ
    CSV.write(df_file, df)
end

function run_DFM2(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2)
    df = DataFrame()
    name = "DFM2"
    df_file = "../../../../data/results/baseline/level$level-h$h-J$J/$name.csv"
    #Forecast using a Dynamic Factor Model with interaction terms
    r1 = 1
    r2 = 1
    ŷ, RMSE = dfm_oos_forc(J, h, h_lag, Y, Y_lag, X_lag, X_lag2, r1, r2)
    println("Dynamic Factor with interaction terms r1 = $(r1), r2 = $(r2): \nRMSE = $RMSE")
    df[Symbol(name)] = ŷ
    CSV.write(df_file, df)
end

function run_LAS(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2)
    df = DataFrame()
    name = "LAS"
    df_file = "../../../../data/results/baseline/level$level-h$h-J$J/$name.csv"
    #Forecast using a LASSO model
    α = 1.0
    λ_vec = 0.1:0.1:5.
    β̂, λ, ŷ, RMSE, RMSE_cv = glmnet_oos_forc(J, h, h_lag, J_cv, Y, Y_lag, X_lag, λ_vec, α)
    println("Tuning Parameter Chosen for lasso: $λ")
    println("Percent of Slope Coefficients Set Equal to Zero by lasso: $(round(100 * (sum(β̂[2:end] .==  .0) / length(β̂[2:end])), 2))")
    println("RMSE from lasso Including all Variables, h = $h: \nRMSE = $RMSE")
    println("RMSE_cv is $RMSE_cv")
    df[Symbol(name)] = ŷ
    CSV.write(df_file, df)
end

function run_LAS2(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2)
    df = DataFrame()
    name = "LAS2"
    df_file = "../../../../data/results/baseline/level$level-h$h-J$J/$name.csv"
    #Forecast using a LASSO model with interaction terms
    α = 1.0
    λ_vec = 0.1:0.1:5.
    X_lag = hcat(X_lag, X_lag2)
    β̂, λ, ŷ, RMSE, RMSE_cv = glmnet_oos_forc(J, h, h_lag, J_cv, Y, Y_lag, X_lag, λ_vec, α)
    println("Tuning Parameter Chosen for lasso: $λ")
    println("Percent of Slope Coefficients Set Equal to Zero by lasso: $(round(100 * (sum(β̂[2:end] .==  .0) / length(β̂[2:end])), 2))")
    println("RMSE from lasso Including all Variables, h = $h: \nRMSE = $RMSE")
    println("RMSE_cv is $RMSE_cv")
    df[Symbol(name)] = ŷ
    CSV.write(df_file, df)
end

function run_RDG(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2)
    df = DataFrame()
    name = "RDG"
    df_file = "../../../../data/results/baseline/level$level-h$h-J$J/$name.csv"
    #Forecast using a Ridge model with interaction terms
    α = 0.0
    λ_vec = 10:0.1:20
    β̂, λ, ŷ, RMSE, RMSE_cv = glmnet_oos_forc(J, h, h_lag, J_cv, Y, Y_lag, X_lag, λ_vec, α)
    println("Tuning Parameter Chosen for ridge: $λ")
    println("RMSE from ridge including first order terms: \nRMSE = $RMSE")
    println("RMSE_cv is $RMSE_cv")
    df[Symbol(name)] = ŷ
    CSV.write(df_file, df)
end

function run_RDG2(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2)
    df = DataFrame()
    name = "RDG2"
    df_file = "../../../../data/results/baseline/level$level-h$h-J$J/$name.csv"
    #Forecast using a Ridge model with interaction terms
    α = 0.0
    λ_vec = 10:0.1:20
    X_lag = hcat(X_lag, X_lag2)
    β̂, λ, ŷ, RMSE, RMSE_cv = glmnet_oos_forc(J, h, h_lag, J_cv, Y, Y_lag, X_lag, λ_vec, α)
    println("Tuning Parameter Chosen for ridge: $λ")
    println("RMSE from ridge including second order terms: \nRMSE = $RMSE")
    println("RMSE_cv is $RMSE_cv")
    df[Symbol(name)] = ŷ
    CSV.write(df_file, df)
end

function run_RDF(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2)
    df = DataFrame()
    name = "RDF"
    df_file = "../../../../data/results/baseline/level$level-h$h-J$J/$name.csv"
    #Forecast using a Random Forrest Regression model
    m_vec = 1:1:25
    n_est = 50
    m, ŷ, RMSE, RMSE_cv = rfr_oos_forc(J , h, h_lag, J_cv, Y, Y_lag, X_lag, m_vec, n_est)
    println("Depth parameter chosen for Random Forest: $m")
    println("RMSE from Random Forest Including all Variables, h = $h: \nRMSE = $RMSE")
    df[Symbol(name)] = ŷ
    CSV.write(df_file, df)
end

function run_RDF2(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2)
    df = DataFrame()
    name = "RDF2"

    df_file = "../../../../data/results/baseline/level$level-h$h-J$J/$name.csv"
    #Forecast using a Random Forrest Regression model
    m_vec = 1:1:25
    n_est = 50
    X_lag = hcat(X_lag, X_lag2)
    m, ŷ, RMSE, RMSE_cv = rfr_oos_forc(J , h, h_lag, J_cv, Y, Y_lag, X_lag, m_vec, n_est)
    println("Depth parameter chosen for Random Forest with interaction terms: $m")
    println("RMSE from Random Forest Including all Variables, h = $h: \nRMSE = $RMSE")
    df[Symbol(name)] = ŷ
    CSV.write(df_file, df)
end

function run_TRD(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2)
    df = DataFrame()
    name = "TRD"
    df_file = "../../../../data/results/baseline/level$level-h$h-J$J/$name.csv"
    df[:REAL] = Y[J - (h_lag + h):end]

    #Forecast using the dynamic factor model
    r = 1
    ŷ, RMSE = dfm_oos_forc(J, h, h_lag, Y, Y_lag, X_lag, r)
    println("Dynamic Factor r = $r: \nRMSE = $RMSE")
    df[:DFM] = ŷ

    #Forecast using AR1 model #
    ŷ, RMSE = AR_oos_forc(J, h, h_lag, Y, Y_lag)
    println("AR(1) Model: \nRMSE = $RMSE")
    df[:ARM] =  ŷ

    #Forecast using OLS with Lagged Aggregate Variable
    ŷ, RMSE = OLS_oos_forc(J, h, h_lag, Y, Y_lag, X_lag)
    println("OLS Regression Including all Variables: \nRMSE = $RMSE")
    df[:OLS] = ŷ

    #Model Average out of sample estimation
    Q = 1
    ŷ, RMSE = modelavg_oos_forc(J, h, h_lag, Y, Y_lag, X_lag, Q)
    println("Equal Weights Model Averaging with Q = $Q: \nRMSE = $RMSE")
    df[:MAG] = ŷ

    #Model Average with interactions
    Q = 1
    X_lag = hcat(X_lag, X_lag2)
    ŷ, RMSE = modelavg_oos_forc(J, h, h_lag, Y, Y_lag, X_lag, Q)
    println("Equal Weights Model Averaging with interactions, Q = $Q: \nRMSE = $RMSE")
    df[:MAG2] = ŷ

    #Forecast using the random walk model
    ŷ, RMSE = rwm_oos_forc(J, h, h_lag, Y, Y_lag)
    println("Random Walk Model: \nRMSE = $RMSE")
    df[:RWM] = ŷ
    CSV.write(df_file, df)
end

max_stars_repo_name: jungangl/forecast-machine-learning

max_issues_repo_name: jungangl/forecast-machine-learning

max_forks_repo_name: jungangl/forecast-machine-learning


text:
using Plots

p = plot()
plot!(p, 1:5, rand(5))
plot!(p, 1:5, rand(5))
vline!(p, [1, 2, 3])