text:
## Load the packages.
using LoopVectorization
using Statistics
using BenchmarkTools
using PyPlot

## Create the refine function.
function refine_field_int!(hi, hi_tmp, lo, n_hi, n_lo)
    if n_hi ÷ n_lo == 2
        @tturbo for j in 0:size(lo, 2)-3
            for i in 0:size(lo, 1)-3
                i_lo = i + 2
                j_lo = j + 2
                i_hi = i*2 + 2
                j_hi = j*2 + 2
                for jj in 0:1, ii in 0:1
                    hi_tmp[i_hi+ii, j_hi+jj] = lo[i_lo, j_lo]
                end
            end
        end

        hi_tmp[1, :] = hi_tmp[end-1, :]
        hi_tmp[end, :] = hi_tmp[2, :]
        hi_tmp[:, 1] = hi_tmp[:, end-1]
        hi_tmp[:, end] = hi_tmp[:, 2]

        # hi[:, :] = hi_tmp[:, :]

        @tturbo for j in 2:size(hi, 2)-1
            for i in 2:size(hi, 1)-1
                hi[i, j] = ( 1//16*hi_tmp[i-1, j-1] + 1//8*hi_tmp[i, j-1] + 1//16*hi_tmp[i+1, j-1]
                           + 1// 8*hi_tmp[i-1, j  ] + 1//4*hi_tmp[i, j  ] + 1// 8*hi_tmp[i+1, j  ]
                           + 1//16*hi_tmp[i-1, j+1] + 1//8*hi_tmp[i, j+1] + 1//16*hi_tmp[i+1, j+1] )
            end
        end
    elseif n_hi ÷ n_lo == 3
        @tturbo for j in 0:size(lo, 2)-3
            for i in 0:size(lo, 1)-3
                i_lo = i + 2
                j_lo = j + 2
                i_hi = i*3 + 2
                j_hi = j*3 + 2
                for jj in 0:2, ii in 0:2
                    hi_tmp[i_hi+ii, j_hi+jj] = lo[i_lo, j_lo]
                end
            end
        end

        hi_tmp[1, :] = hi_tmp[end-1, :]
        hi_tmp[end, :] = hi_tmp[2, :]
        hi_tmp[:, 1] = hi_tmp[:, end-1]
        hi_tmp[:, end] = hi_tmp[:, 2]

        @tturbo for j in 2:size(hi, 2)-1
            for i in 2:size(hi, 1)-1
                hi[i, j] = ( 1//9*hi_tmp[i-1, j-1] + 1//9*hi_tmp[i, j-1] + 1//9*hi_tmp[i+1, j-1]
                           + 1//9*hi_tmp[i-1, j  ] + 1//9*hi_tmp[i, j  ] + 1//9*hi_tmp[i+1, j  ]
                           + 1//9*hi_tmp[i-1, j+1] + 1//9*hi_tmp[i, j+1] + 1//9*hi_tmp[i+1, j+1] )
            end
        end
    else
        throw(DomainError(n_hi/n_lo, "Refinement should be 2 or 3"))
    end
end

## Set up the grids.
n_hi = 12
n_lo = n_hi ÷ 3

a_lo = rand(n_lo, n_lo)
a_hi = zeros(n_hi, n_hi)

a_hi_gc = zeros(n_hi+2, n_hi+2)
a_lo_gc = zeros(n_lo+2, n_lo+2)

a_hi_gc[2:n_hi+1, 2:n_hi+1] = a_hi[:, :]
a_lo_gc[2:n_lo+1, 2:n_lo+1] = a_lo[:, :]

a_hi_gc[1, :] = a_hi_gc[end-1, :]
a_hi_gc[end, :] = a_hi_gc[2, :]
a_hi_gc[:, 1] = a_hi_gc[:, end-1]
a_hi_gc[:, end] = a_hi_gc[:, 2]

a_lo_gc[1, :] = a_lo_gc[end-1, :]
a_lo_gc[end, :] = a_lo_gc[2, :]
a_lo_gc[:, 1] = a_lo_gc[:, end-1]
a_lo_gc[:, end] = a_lo_gc[:, 2]

dx_hi = 1/n_hi; dx_lo = 1/n_lo
x_hi = dx_hi/2:dx_hi:1 |> collect
x_lo = dx_lo/2:dx_lo:1 |> collect

## Compute.
a_hi_tmp = zeros(size(a_hi_gc))
@btime refine_field_int!(a_hi_gc, a_hi_tmp, a_lo_gc, n_hi, n_lo)
a_hi_int = a_hi_gc[2:end-1, 2:end-1]

println(mean(a_lo) ≈ mean(a_hi_int))

## Plot the output.
figure()
plot(x_hi, a_hi_int[1, :], "C0-o")
plot(x_hi, a_hi_int[2, :], "C1-+")
if n_hi ÷ n_lo == 3
    plot(x_hi, a_hi_int[3, :], "C2-^")
end
plot(x_lo, a_lo[1, :], "k:")
display(gcf())

xh_hi = 0:dx_hi:1 |> collect
xh_lo = 0:dx_lo:1 |> collect

figure(figsize=(10, 5))
subplot(121)
pcolormesh(xh_lo, xh_lo, a_lo, vmin=0, vmax=1)
subplot(122)
pcolormesh(xh_hi, xh_hi, a_hi_int, vmin=0, vmax=1)
tight_layout()
display(gcf())
show()


text:
module Aerostructural

## Package imports
#==========================================================================================#

using StaticArrays
using LinearAlgebra
using ComponentArrays
using TimerOutputs
using SparseArrays
using SplitApplyCombine

# Structures
import ..Beams: tube_stiffness_matrix, AbstractBeam, Beam, Tube

# Conversions
import ..Laplace: freestream_to_cartesian

# Panelling
import ..PanelGeometry: make_panels, normal_vector

import ..AircraftGeometry: WingMesh

# VLM Aerodynamics
import ..VortexLattice: velocity, induced_velocity, induced_trailing_velocity, Horseshoe, horseshoe_normal, collocation_point, bound_leg_center, influence_coefficient, influence_matrix, boundary_condition, geometry_to_wind_axes, bound_leg_vector, kutta_joukowsky, surface_forces

import ..AeroMDAO: solve_linear, solve_nonlinear, solve_nonlinear!

struct AerostructWing{S,T}
    aerodynamics :: WingMesh{S}
    structures   :: Beam{T}
end

struct AerostructSystem{T}
    vortices :: Array{Horseshoe{T}}

end

## Aerodynamic analysis
#==========================================================================================#

include("aerodynamics.jl")

## Structural analysis
#==========================================================================================#

include("structures.jl")

## Load-displacement transfer mechanisms
#==========================================================================================#

include("transfers.jl")

## Weights, engine, and fuel loads
#==========================================================================================#

# ???

## Coupled residual systems
#==========================================================================================#

include("residuals.jl")

end

ext: jl


text:
"""
Basic representation of the state of the age structured model:
u[wa_index,age_group,disease state]

This is row major unpacked so
first 1,...,n_wa entries are 0-4 year old susceptibles in wider areas 1,...,n_wa
then n_wa+1,...,2n_wa are 5-9 year old susceptibles in wider areas 1,...,n_wa
.
.
.
then n_wa*n_s + 1, ...,  n_wa*n_s + n_wa entries are 0-4 year old exposed in wider areas 1,...,n_wa
.
.
.

States:
1 -> S
2 -> E
3 -> P
4 -> A(symptomatic)
5 -> M(ild)
6 -> (se)V(ere)
7 -> R(ecovered)
8 -> H(ospitalised)
9 -> Cumulative A
10 -> Cumulative M
11 -> Cumulative V
12 -> Cum. V->H

Events for each wider area and age group:

1-> Transmission
2-> Incubation into Prodromal
3 -> P->A
3-> P->M
4-> P->V
....
5-> Hospitalised/treated recover
6-> Diseased recover D->R
7-> Asymptomatics recover A->R
8->Hospitalised/treated->death

"""

dc_age = zeros(Int64,n_wa*n_a*n_s,n_ta*n*n_a)
# asymp_indices = zeros(Bool,n_wa,n_a,n_s)
# asymp_indices[:,:,3] .= true;
# f_asymp_indices = findall(asymp_indices[:])
# diseased_indices = zeros(Bool,n_wa,n_a,n_s)
# diseased_indices[:,:,4] .= true;
# f_diseased_indices = findall(diseased_indices[:])
# asymp_indices = 0;#free memory
# diseased_indices = 0;

function calculate_infection_rates!(u,p::CoVParameters_AS,t)
    I_A = @view u[:,:,3]
    I_D = @view u[:,:,4]
    mul!(p.Î,p.T,p.ϵ*I_A .+ p.ϵ_D*I_D)#Local infecteds **if** everyone moved around
    p.Î[:,immobile_age_indices] .= p.ϵ*I_A[:,immobile_age_indices] .+ p.ϵ_D*I_D[:,immobile_age_indices]#This corrects for immobility
    mul!(p.λ_loc,p.β*p.c_t(t).*(p.Î ./p.N̂),p.M)#Local force of infection due to age-mixing --- M is in to (row), from (col) format
    mul!(p.λ,p.T',p.λ_loc)#this accounts for mobile susceptibles contracting away from home
    p.λ[:,immobile_age_indices] .= p.λ_loc[:,immobile_age_indices]#This corrects for immobility
    p.λ[ind_mombasa_as,:] .+= p.ext_inf_rate*import_rate_mom(t,p.into_mom,p.global_prev)
    p.λ[ind_nairobi_as,:] .+= p.ext_inf_rate*import_rate_nai(t,p.into_nai,p.global_prev)
    return nothing
end

function rates(out,u,p::CoVParameters_AS,t)
    @unpack λ,γ,σ,δ,τ,μ₁,χ,rel_detection_rate,clear_quarantine = p
    calculate_infection_rates!(u,p,t)
    for k = 1:length(out)
        i,a,eventtype = Tuple(index_as_events[k])
        if eventtype ==1
            out[k] = χ[a]*λ[i,a]*u[i,a,1] #Transmission
        end
        if eventtype ==2
            out[k] = (1-(δ*rel_detection_rate[a]))*σ*u[i,a,2] # E->A
        end
        if eventtype ==3
            out[k] = δ*rel_detection_rate[a]*σ*u[i,a,2] # E->D
        end
        if eventtype ==4
            out[k] = τ*u[i,a,4] # D->H
        end
        if eventtype ==5
            out[k] = clear_quarantine*u[i,a,5] # H->R
        end
        if eventtype ==6
            out[k] = γ*u[i,a,4] # D->R
        end
        if eventtype ==7
            out[k] = γ*u[i,a,3] # A->R
        end
        if eventtype ==8
            out[k] = μ₁*u[i,a,5] # H->death
        end
    end
end

function change_matrix(dc)
    d1,d2 = size(dc)
    for k = 1:d2
        i,a,eventtype = Tuple(index_as_events[k])
        if eventtype ==1 #Transmission
            ind_S = linear_as[i,a,1]
            ind_E = linear_as[i,a,2]
            dc[ind_S,k] = -1
            dc[ind_E,k] = 1
        end
        if eventtype ==2 # E->A
            ind_E = linear_as[i,a,2]
            ind_A = linear_as[i,a,3]
            ind_cumA = linear_as[i,a,7]
            dc[ind_E,k] = -1
            dc[ind_A,k] = 1
            dc[ind_cumA,k] = 1
        end
        if eventtype ==3 # E->D
            ind_E = linear_as[i,a,2]
            ind_D = linear_as[i,a,4]
            ind_cumD = linear_as[i,a,8]
            dc[ind_E,k] = -1
            dc[ind_D,k] = 1
            dc[ind_cumD,k] = 1
        end
        if eventtype ==4# D->H
            ind_D = linear_as[i,a,4]
            ind_H = linear_as[i,a,5]
            dc[ind_D,k] = -1
            dc[ind_H,k] = 1
        end
        if eventtype ==5 # H->R
            ind_H = linear_as[i,a,5]
            ind_R = linear_as[i,a,6]
            dc[ind_H,k] = -1
            dc[ind_R,k] = 1
        end
        if eventtype ==6 # D->R
            ind_D = linear_as[i,a,4]
            ind_R = linear_as[i,a,6]
            dc[ind_D,k] = -1
            dc[ind_R,k] = 1
        end
        if eventtype ==7 # A->R
            ind_A = linear_as[i,a,3]
            ind_R = linear_as[i,a,6]
            dc[ind_A,k] = -1
            dc[ind_R,k] = 1
        end
        if eventtype ==8 # H->death
            ind_H = linear_as[i,a,5]
            ind_cumdead = linear_as[i,a,9]
            dc[ind_H,k] = -1
            dc[ind_cumdead,k] = 1
        end
    end
end

#This caps the Poisson processes at causing no more transitions than the group being effected
function max_change(out,u,p::CoVParameters_AS)
    for i = 1:n_wa,a = 1:n_a
        ind_trans = linear_as_events[i,a,1]
        ind_EA = linear_as_events[i,a,2]
        ind_ED = linear_as_events[i,a,3]
        ind_DH = linear_as_events[i,a,4]
        ind_HR = linear_as_events[i,a,5]
        ind_DR = linear_as_events[i,a,6]
        ind_AR = linear_as_events[i,a,7]
        ind_Hdeath = linear_as_events[i,a,8]

        out[ind_trans] = min(out[ind_trans],u[i,a,1])
        out[ind_HR] = min(out[ind_HR],u[i,a,5])
        out[ind_AR] = min(out[ind_AR],u[i,a,3])
        #spliting events using binomial sampling
        if out[ind_EA] + out[ind_ED] > u[i,a,2] #More incubations than actual rural E population
            out[ind_ED] = rand(Binomial(u[i,a,2],p.δ)) #Binomially distributed the rural incubations between Asympotomatic and symptomatic
            out[ind_EA] =  u[i,a,2] - out[ind_ED]
        end
        if out[ind_DR] + out[ind_DH] > u[i,a,4] # More end of infection events than diseased/symptomatic infecteds
            out[ind_DH] = rand(Binomial(u[i,a,4],p.τ/(p.τ + p.γ))) #Binomially distributed the end of infections between hospitalisation and recovery
            out[ind_DR] = u[i,a,4] - out[ind_DH]
        end
    end
end

function nonneg_tauleap(du,u,p::CoVParameters_AS,t)
    @unpack dc,dN,poi_rates,du_linear = p
    rates(poi_rates,u,p,t) #calculate rates of underlying Poisson processes
    PP_drivers(dN,poi_rates,p)#Generate Poisson rvs with rates scaled by time step dt
    max_change(dN,u,p)#Cap the size of the Poisson rvs to maintain non-negativity
    mul!(du_linear,dc,dN)#Calculates the effect on the state in the inplace du vector
    du .= reshape(du_linear,n_wa,n_a,n_s)
    du .+= u #Calculates how the state should change
end

function ode_model(du,u,p::CoVParameters_AS,t)
    @unpack λ,γ,σ,δ,τ,μ₁,χ,rel_detection_rate,clear_quarantine = p
    calculate_infection_rates!(u,p,t)
    for i = 1:n_wa,a = 1:n_a
        du[i,a,1] = (-1)*χ[a]*λ[i,a]*u[i,a,1]
        du[i,a,2] = χ[a]*λ[i,a]*u[i,a,1] - σ*u[i,a,2]
        du[i,a,3] = (1-(δ*rel_detection_rate[a]))*σ*u[i,a,2] - γ*u[i,a,3]
        du[i,a,4] = δ*rel_detection_rate[a]*σ*u[i,a,2] - γ*u[i,a,4] - τ*u[i,a,4]
        du[i,a,5] = τ*u[i,a,4] - clear_quarantine*u[i,a,5] - μ₁*u[i,a,5]
        du[i,a,6] = clear_quarantine*u[i,a,5] + γ*(u[i,a,4] + u[i,a,3])
        du[i,a,7] = (1-δ)*σ*u[i,a,2]
        du[i,a,8] = δ*σ*u[i,a,2]
        du[i,a,9] = μ₁*u[i,a,5]
     end
end


text:
module WANPerformanceAnalyzer

export Trace, Workload, rate, drate

include("traces.jl")

end

max_stars_repo_name: BenLauwens/WANPerformanceAnalyzer

max_issues_repo_name: BenLauwens/WANPerformanceAnalyzer

max_forks_repo_name: BenLauwens/WANPerformanceAnalyzer


text:
function _decamelize(cur_string::AbstractString)
  capital_letters = matchall(r"[A-Z]+", cur_string)
  capital_letters = map(x -> downcase(x), capital_letters)

  lowercase_phrases = split(cur_string, r"([A-Z]+)")

  cur_array = vcat(
    map(
      x -> join(x, "_"),
      [zip(lowercase_phrases, capital_letters)...]
    )...
  )

  cur_string = join(cur_array)

  if length(lowercase_phrases) > length(capital_letters)
    cur_string *= lowercase_phrases[end]
  end

  if length(lowercase_phrases) < length(capital_letters)
    error("Improperly handled string on decamelize: $cur_string")
  end

  ( startswith(cur_string, "_") ) && ( cur_string = cur_string[2:end] )
  ( startswith(cur_string, "-") ) && ( cur_string = cur_string[2:end] )

  cur_string
end


text:
export ispalindrome_string
export ispalindrome_integer
export ispalindrome
export reverse_sum
export reverse_sum_integer
export digitindex
export integer_to_string

"""
Returns whether a number is a palindrome or not, using String operations
to pull out digits.
"""
function ispalindrome_string(n::Integer)::Bool
    stringval = string(n)
    return stringval == reverse(stringval)
end

"""
Returns whether a number is a palindrome or not, using normal integer
operations (floor, mod) to pull out digits.
"""
function ispalindrome_integer(n::Integer)::Bool
    N = ndigits(n) - 1
    for d in N:-2:1
        pow = convert(typeof(n), 10)^d
        if fld(n, pow) != mod(n, 10)
            return false
        end
        n = fld(n - fld(n, pow) * pow, 10)
    end

    return true
end

"""
Returns whether a number is a palindrome or not, using Julia's built-in digits
function.
"""
function ispalindrome(n::Integer)::Bool
    ds = digits(n)
    if ds[1] == 0
        return false
    end
    numdigits = length(ds)
    for d in 1:fld(numdigits, 2)
        @inbounds if ds[d] != ds[numdigits-d+1]
            return false
        end
    end

    return true
end

"""
Returns the sum of the number and its reverse using Julia builtins.
"""
function reverse_sum(n::Integer)::Integer
    digs = digits(n)
    n_digits = length(digs)
    multipliers = convert(typeof(n), 10) .^ (n_digits-1:-1:0)
    return n + sum(multipliers .* digs)
end

"""
Returns the sum of the number and its reverse using mods and divs instead of
Julia built-ins.
"""
function reverse_sum_integer(n::Integer)::Integer
    n_digits = ndigits(n) - 1

    pow = convert(typeof(n), 10)^n_digits
    s = n
    for d in 1:n_digits
        s += (pow * mod(n, 10))
        n = div(n, 10)
        pow = div(pow, 10)
    end
    return s + n
end

"""
Returns the digit at index, e.g., digitindex(421513, 4) = 5.
"""
function digitindex(n::Integer, index::Integer)::Integer
    digs = digits(n)
    return digs[length(digs) - index + 1]
end

"""
Returns a plain English string version of a number.
"""
function integer_to_string(n::Integer)::String
    ones_words = Dict{Integer, String}( 1 => "one",
                                        2 => "two",
                                        3 => "three",
                                        4 => "four",
                                        5 => "five",
                                        6 => "six",
                                        7 => "seven",
                                        8 => "eight",
                                        9 => "nine")
    teen_words = Dict{Integer, String}(11 => "eleven",
                                       12 => "twelve",
                                       13 => "thirteen",
                                       14 => "fourteen",
                                       15 => "fifteen",
                                       16 => "sixteen",
                                       17 => "seventeen",
                                       18 => "eighteen",
                                       19 => "nineteen")
    tens_words = Dict{Integer, String}(10 => "ten",
                                       20 => "twenty",
                                       30 => "thirty",
                                       40 => "fourty",
                                       50 => "fifty",
                                       60 => "sixty",
                                       70 => "seventy",
                                       80 => "eighty",
                                       90 => "ninety")
    short_scale_words = Dict{Integer, String}(4 => "thousand",
                                              7 => "million",
                                              10 => "billion",
                                              13 => "trillion",
                                              16 => "quadrillion")
    hundred_word = "hundred"
    connector_word = "and"
    hyphen = "-"

    # Split each group of 3 digits in a large number into groups of 2 and 1
    # digits. So 234_423_032_001 will become [(2, 34), (4, 23), (0, 32), (0, 1)]
    # Rules for labeling are the same within each triplet of digits.
    desc = ""
    for ndig in 1:3:ndigits(n)
        hundreds_digit = mod(fld(n, 100), 10)
        tens_digits = mod(n, 100)
        triplet_desc = ""
        if haskey(ones_words, hundreds_digit)
            triplet_desc = triplet_desc * ones_words[hundreds_digit] *
                           " $hundred_word "
            if (tens_digits > 0)
                triplet_desc = triplet_desc * "$connector_word "
            end
        end
        if haskey(tens_words, tens_digits)
            # Ends in a multiple of ten
            triplet_desc = triplet_desc * tens_words[tens_digits] * " "
        elseif haskey(teen_words, tens_digits)
            # Ends in a teen
            triplet_desc = triplet_desc * teen_words[tens_digits] * " "
        else
            # Ends in something else
            ones_digit = mod(n, 10)
            tens_digit = tens_digits - ones_digit
            tens_suffix = ""
            if haskey(tens_words, tens_digit)
                if haskey(ones_words, ones_digit)
                    tens_suffix = tens_words[tens_digit] * "$hyphen"
                else
                    tens_suffix = tens_words[tens_digit] * " "
                end
            end
            if haskey(ones_words, ones_digit)
                tens_suffix = tens_suffix * ones_words[ones_digit] * " "
            end
            triplet_desc = triplet_desc * tens_suffix
        end
        if haskey(short_scale_words, ndig)
            triplet_desc = triplet_desc * short_scale_words[ndig] * " "
        end
        desc = triplet_desc * desc
        n = fld(n, 1000)
    end
    return desc
end

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
function solve()
    n = readline()
    ret = 0
    for i in 2:2:length(n)
        ret += if i < length(n)
            9 * 10^(i÷2-1)
        else
            h = i÷2
            f = parse(Int, n[1:h])
            s = parse(Int, n[h+1:end])
            x = s ≥ f ? f : f-1
            x - 10^(h-1) + 1
        end
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
@testset "CTQW model check" begin
  n = 5
  g = complete_graph(n)
  ctqw = CTQW(g)
  ctqw_qss = QWSearch(ctqw, [1,3])

  @test isapprox(sum(measure(ctqw_qss, initial_state(ctqw_qss))), 1.0)
  @test QuantumWalk.matrix(ctqw) == ctqw.matrix

  @testset "Default jumping rate" begin
    @test QuantumWalk.jumping_rate(ctqw) ≈ 1. /(n-1)
    @test QuantumWalk.jumping_rate(CTQW(Graph(n))) == 1.
    @test_throws AssertionError QuantumWalk.jumping_rate(CTQW(g, :laplacian))
  end

  @testset "CTQW conversion" begin
    ctqw = CTQW(smallgraph(:bull))
    qws = QWSearch(ctqw, [1], 0.5)

    function is_equal_mine(q1::QWSearch, q2::QWSearch)
      model(q1) == model(q2) &&
      parameters(q1)[:hamiltonian] ≈ parameters(q2)[:hamiltonian] &&
      penalty(q1) == penalty(q2) &&
      marked(q1) == marked(q2)
    end

    qss1 = QWSearch(qws, marked = [1])
    qss2 = QWSearch(qws, penalty = 0.5)
    qss3 = QWSearch(qws, marked = [2])
    qss4 = QWSearch(qws, penalty = 1.)
    qss5 = QWSearch(qws, marked = [2], penalty = 1.)
    qss6 = QWSearch(qws, marked = [2], penalty = 2.)

    qss1_ref = QWSearch(ctqw, [1], 0.5)
    qss2_ref = QWSearch(ctqw, [1], 0.5)
    qss3_ref = QWSearch(ctqw, [2], 0.5)
    qss4_ref = QWSearch(ctqw, [1], 1.)
    qss5_ref = QWSearch(ctqw, [2], 1.)

    @test is_equal_mine(qss1, qss1_ref)
    @test is_equal_mine(qss2, qss2_ref)
    @test is_equal_mine(qss3, qss3_ref)
    @test is_equal_mine(qss4, qss4_ref)
    @test is_equal_mine(qss5, qss5_ref)
    @test !is_equal_mine(qss6, qss5_ref)
  end
end

@testset "CTQWDense model check" begin
  n = 5
  g = complete_graph(n)
  ctqw = CTQWDense(g)
  ctqw_qss = QWSearch(ctqw, [1,3])

  @test isapprox(sum(measure(ctqw_qss, initial_state(ctqw_qss))), 1.0)
  @test QuantumWalk.matrix(ctqw) == ctqw.matrix

  @testset "Default jumping rate" begin
    @test QuantumWalk.jumping_rate(ctqw) ≈ 1. /(n-1)
    @test QuantumWalk.jumping_rate(CTQWDense(Graph(n))) == 1.
    @test_throws AssertionError QuantumWalk.jumping_rate(CTQWDense(g, :laplacian))
  end

  @testset "CTQWDense conversion" begin
    ctqw = CTQWDense(smallgraph(:bull))
    qws = QWSearch(ctqw, [1], 0.5)

    function is_equal_mine(q1::QWSearch, q2::QWSearch)
      model(q1) == model(q2) &&
      parameters(q1)[:hamiltonian] ≈ parameters(q2)[:hamiltonian] &&
      penalty(q1) == penalty(q2) &&
      marked(q1) == marked(q2)
    end

    qss1 = QWSearch(qws, marked = [1])
    qss2 = QWSearch(qws, penalty = 0.5)
    qss3 = QWSearch(qws, marked = [2])
    qss4 = QWSearch(qws, penalty = 1.)
    qss5 = QWSearch(qws, marked = [2], penalty = 1.)
    qss6 = QWSearch(qws, marked = [2], penalty = 2.)

    qss1_ref = QWSearch(ctqw, [1], 0.5)
    qss2_ref = QWSearch(ctqw, [1], 0.5)
    qss3_ref = QWSearch(ctqw, [2], 0.5)
    qss4_ref = QWSearch(ctqw, [1], 1.)
    qss5_ref = QWSearch(ctqw, [2], 1.)

    @test is_equal_mine(qss1, qss1_ref)
    @test is_equal_mine(qss2, qss2_ref)
    @test is_equal_mine(qss3, qss3_ref)
    @test is_equal_mine(qss4, qss4_ref)
    @test is_equal_mine(qss5, qss5_ref)
    @test !is_equal_mine(qss6, qss5_ref)
  end
end


text:
################################################################################
######      Cache holding the information generated along a sampling      ######
################################################################################
mutable struct MCMCSRLEvaluationCache{T,T2,TV,TVC,TM,TD,S} <: EvaluationSamplingCache
    Oave::TV#Vector{T}
    OOave::TM#Matrix{T}
    Eave::T
    EOave::TVC
    LLOave::TVC
    Zave::T2

    # Individual values to compute statistical correlators
    Evalues::Vector{T}

    # Caches to avoid allocating during computation
    LLO_i::TVC
    ∇lnψ::TD
    ∇lnψ2::TD
    σ::S
end

function MCMCSRLEvaluationCache(net::NeuralNetwork, prob)
    TC = Complex{real(out_type(net))}
    der_vec = grad_cache(net).tuple_all_weights

    Oave   = [zero(dvec) for dvec=der_vec]
    OOave  = [zeros(eltype(dvec), length(dvec), length(dvec)) for dvec=der_vec]
    EOave  = [zeros(TC, size(dvec)) for dvec=der_vec]
    LLOave = [zeros(TC, size(dvec)) for dvec=der_vec]
    LLO_i  = [zeros(TC, size(dvec)) for dvec=der_vec]

    cache = MCMCSRLEvaluationCache(Oave,
                                  OOave,
                                  zero(TC),
                                  EOave,
                                  LLOave,
                                  zero(real(TC)),
                                  Vector{TC}(),
                                  LLO_i,
                                  grad_cache(net),
                                  grad_cache(net),
                                  state(prob, net))
    zero!(cache)
end

SamplingCache(alg::SR, prob::LRhoSquaredProblem, net) = MCMCSRLEvaluationCache(net, prob)

function zero!(comp_vals::MCMCSRLEvaluationCache)
    comp_vals.Eave   = 0.0
    comp_vals.Zave = 0.0
    resize!(comp_vals.Evalues, 0)

    for i=1:length(comp_vals.Oave)
        comp_vals.Oave[i]   .= 0.0
        comp_vals.OOave[i]  .= 0.0
        comp_vals.EOave[i]  .= 0.0
        comp_vals.LLOave[i] .= 0.0
    end

    comp_vals
end

# Utility method utilised to accumulate results on a single variable
function add!(acc::MCMCSRLEvaluationCache, o::MCMCSRLEvaluationCache)
    acc.Eave   += o.Eave
    acc.Zave   += o.Zave
    append!(acc.Evalues, o.Evalues)

    for i=1:length(acc.Oave)
       acc.Oave[i]   .+= o.Oave[i]
       acc.OOave[i]  .+= o.OOave[i]
       acc.EOave[i]  .+= o.EOave[i]
       acc.LLOave[i] .+= o.LLOave[i]
   end
   acc
end

#TODO Things that allocate in this function are
# out.LVals ./= 1:length(vals.Evalues)
# S .= real.(OOave .- (conj(Oave) .* transpose(Oave)))
# F .= real.(EOave .- Eave .* conj(Oave))
# and they clearly should not. It's julia's fault.
##
function evaluation_post_sampling!(out::SREvaluation,
                                   vals::MCMCSRLEvaluationCache,
                                   sampler_steps = vals.Zave)
     Eave = vals.Eave   /= sampler_steps

     #TODO Here I am reallocating. Should think about how to fix it.
     out.L = Eave
     # Convergence history
     #out.LVals = cumsum(vals.Evalues)./(1:length(vals.Evalues))
     resize!(out.LVals, length(vals.Evalues))
     cumsum!(out.LVals, vals.Evalues)
     out.LVals ./= 1:length(vals.Evalues)

     for i=1:length(vals.Oave)
         Oave   = vals.Oave[i]   ./= sampler_steps
         OOave  = vals.OOave[i]  ./= sampler_steps
         EOave  = vals.EOave[i]  ./= sampler_steps
         LLOave = vals.LLOave[i] ./= sampler_steps

         S = out.S[i]; F = out.F[i];
         if eltype(S) <: Real
             S .= real.(OOave  .- (conj(Oave) .* transpose(Oave)))
             F .= real.(LLOave .- Eave .* conj(Oave))
         else
             S .=       OOave  .- (conj(Oave) .* transpose(Oave))
             F .=       LLOave .- Eave .* conj(Oave)
         end
     end
end


text:
@testset "228.summary-ranges.jl" begin
    @test max_profit_iii([1]) == 0
    @test max_profit_iii([7,6,4,3,1]) == 0
    @test max_profit_iii([91, 11, 6, 40, 14, 56, 30, 60, 80, 69, 19, 22, 35, 66, 98, 84, 79, 38, 20, 48, 48, 4, 50, 42, 95, 16, 27, 61, 52, 12, 84, 80, 64, 61, 88, 42, 39, 1, 62, 55, 31, 34, 71, 21, 14, 83, 57, 94, 10, 69, 56, 66, 32, 46, 29, 92, 70, 2, 21, 10, 57, 42, 74, 97, 51, 51, 52, 91, 63, 72, 89, 82, 93, 74, 97, 89, 69, 25, 51, 71, 9, 51, 37, 3, 53, 38, 13, 86, 30, 78, 12, 33, 96, 39, 87, 100, 77, 52, 11, 17, 95, 65, 1, 34, 88, 16, 33, 7, 84, 81, 47, 93, 59, 86, 39, 2, 87, 44, 48, 34, 33, 59, 57, 3, 74, 73, 66, 18, 38, 60, 47, 73, 21, 38, 90, 68, 28, 38, 32, 20, 61, 78, 28, 63, 62, 8, 62, 59, 78, 13, 30, 21, 27, 67, 17, 61, 39, 65, 61, 84, 7, 73, 52, 34, 18, 85, 62, 68, 75, 46, 45, 91, 90, 6, 8, 87, 83, 70, 6, 64, 42, 93, 1, 57, 37, 61, 68, 88, 32, 94, 94, 36, 42, 67, 1, 77, 47, 95, 64, 89]) == 193
    @test max_profit_iii([1,2,3,4,5]) == 4
end


text:
module C2_0
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    -l1_dot.*(3.53674108537383e-11*l1 - 2.16714051041859e-11*l2 - 1.38960161444059e-11*l3 + 3.04085216377704e-14) + l2_dot.*(2.16714051041859e-11*l1 - 2.75081888864459e-11*l2 + 5.69132167691212e-12*l3 + 2.23183844184968e-14) + l3_dot.*(1.38960161444059e-11*l1 + 5.69132167691212e-12*l2 - 1.96418867892903e-11*l3 + 8.44288384244248e-15)
end
end


text:
module GaussianLikelihoods

using ..GPs
using ..Kernels
using ..DataSets
using StaticArrays
using DocStringExtensions

import ..GPs: sethyperparameters!, gethyperparameters, hyperparametercount,
    loglikelihood, loglikelihood_dθ, loglikelihood_dθ!

struct GP{T, NX, NY, K <: Kernel{T}} <: AbstractGP{T}
    kernel::K # covariance function
    σn::MVector{1, T} # additive noise (MVector so it remains mutable)
    training::TrainingSet{NX, NY, T}
    K::UpperTriangular{T, Matrix{T}} # Covariance matrix
    L::UpperTriangular{T, Matrix{T}} # Cholesky factors of the covariance matrix
    α::Matrix{T} # α = K⁻¹y

    function GP(kernel::Kernel, σn::T, training::TrainingSet{NX, NY}) where {T <: Number, NX, NY}
        n = samplecount(training) # number of samples
        gp = new{T, NX, NY, typeof(kernel)}(kernel,
                                            MVector(σn),
                                            training,
                                            UpperTriangular(Matrix{T}(n, n)), # K
                                            UpperTriangular(Matrix{T}(n, n)), # L = chol(K)
                                            Matrix{T}(n, NY)) # α = K⁻¹y
        # Update with the specified hyperparameters/observations
        update!(gp)
        return gp
    end
end

export GP

function sethyperparameters!(gp::GP, θ)
    # Recursively set the hyperparameters
    gp.σn[1] = θ[1]
    sethyperparameters!(gp.kernel, θ[2:end])
    # Update the GP with the new hyperparameters
    update!(gp)
end

gethyperparameters(gp::GP) = [SVector(gp.σn); gethyperparameters(gp.kernel)]

hyperparametercount(gp::GP) = hyperparametercount(gp.kernel) + 1

function update!(gp::GP)
    # Prior covariance using the specified kernel
    covariance_matrix!(gp.K, gp.kernel, gp.training.x)
    # Add measurement noise at each point
    for i ∈ diagind(gp.K)
        gp.K[i] += gp.σn[1]^2
    end
    # Calculate the Cholesky decomposition
    copy!(gp.L, gp.K)
    cholfact!(Symmetric(gp.L)) # assumes L is upper triangular
    # Calculate quantities required for interpolation
    gp.α .= gp.training.ym'
    At_ldiv_B!(gp.L, gp.α) # L' \ y (intermediate)
    A_ldiv_B!(gp.L, gp.α) # α = L \ (L' \ y)
    return gp
end

#--- Posterior measures

import Base: mean, var, cov

"""
$(SIGNATURES)

Return the mean of the Gaussian Process at the specified points.
"""
function mean(gp::GP{T}, x::Vector{SVector{N, T}}) where {T, N}
    Kₛ = covariance_matrix(gp.kernel, x, gp.training.x)
    return At_mul_Bt(gp.α, Kₛ) # (Kₛ*α)ᵀ
end

mean(gp::GP{T}, x::Union{AbstractVector{T}, Matrix{T}}) where {T} = mean(gp, toSVector(x))

"""
$(SIGNATURES)

Return the variance of the Gaussian Process at the specified points.
"""
function var(gp::GP{T}, x::Vector{SVector{N, T}}) where {T, N}
    Kₛ = covariance_matrix(gp.kernel, gp.training.x, x)
    v = gp.L' \ Kₛ
    Kₛₛ = [covariance(gp.kernel, xᵢ, xᵢ) for xᵢ in x]
    return Kₛₛ - vec(sum(v.*v, 1)) + gp.σn[1]^2
end

var(gp::GP{T}, x::Union{AbstractVector{T}, Matrix{T}}) where {T} = var(gp, toSVector(x))

"""
$(SIGNATURES)

Return the covariance matrix of the Gaussian Process at the specified points.
"""
function cov(gp::GP{T}, x::Vector{SVector{N, T}}) where {T, N}
    Kₛ = covariance_matrix(gp.kernel, gp.training.x, x)
    v = gp.L' \ Kₛ
    Kₛₛ = covariance_matrix(gp.kernel, x)
    Kₚ = Kₛₛ - v'*v
    for i ∈ diagind(Kₚ)
        Kₚ[i] += gp.σn[1]^2
    end
    return Kₚ
end

cov(gp::GP{T}, x::Union{AbstractVector{T}, Matrix{T}}) where {T} = cov(gp, toSVector(x))

#--- Likelihoods

function loglikelihood(gp::GP{T}) where {T}
    # log likelihood of the GP itself (w.r.t. the hyperparameters)
    n = samplecount(gp.training)
    # log(det(K))
    logdetK = sum(log.(diag(gp.L)))
    # The log likelihood
    lik = T(0.5)*mean(sum(gp.training.ym'.*gp.α, 1)) + logdetK + T(0.5)*n*log(T(2)*π)
end

function loglikelihood_dθ!(dθ::AbstractVector{T}, gp::GP{T}) where {T}
    n = samplecount(gp.training)
    I = eye(T, n)
    At_ldiv_B!(gp.L, I) # I = gp.L' \ eye(n)
    A_ldiv_B!(gp.L, I) # I = gp.L \ (gp.L' \ eye(n))
    #inner = (gp.α*gp.α' - gp.L \ (gp.L' \ eye(n)))'
    Kdθ = @view dθ[2:end]
    for j ∈ 1:n
        for i ∈ 1:j
            mult = (mean(gp.α[i, :] .* gp.α[j, :]) - I[i, j]) * ifelse(i == j, T(0.5), T(1)) # use ifelse to avoid computing the symmetric counterparts
            Kdθ .-= covariance_dθ(gp.kernel, gp.training.x[i], gp.training.x[j]) .* mult
            dθ[1] -= ifelse(i == j, T(2) * gp.σn[1] * mult, zero(T))
        end
    end
    return dθ
end

loglikelihood_dθ(gp::GP{T}) where {T} = loglikelihood_dθ!(zeros(T, hyperparametercount(gp)), gp)

end


text:
using DiscreteEntropyEstimators
using Base.Test

@test entropy_ml(Frequencies([1, 1, 2])) == entropy(MLEE(), Frequencies([1, 1, 2]))
@test entropy_ml(Frequencies([1, 1, 2]), 2) == 1.5
@test entropy_ml(Frequencies([1, 1, 2]), e) == entropy_ml(Frequencies([1, 1, 2]))
@test entropy_ml(Frequencies([1, 1, 1, 1, 4, 4, 4]), 4) == 1.25
@test entropy_ml(Samples([1, 2, 3, 3])) == entropy_ml(Frequencies([1, 1, 2]))
@test entropy(MMBCEE(), Frequencies([1, 1, 2])) == entropy(MMBCEE(), Samples([1, 2, 3, 3]))


text:
#versioninfo()

#import Pkg; Pkg.add("Arpack")
#import Pkg; Pkg.add("ArgParse")
#import Pkg; Pkg.add("LinearMaps")

using LinearAlgebra
using Arpack
using SparseArrays
using ArgParse
using LinearMaps

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--N"
            help = "set Nsize (should be >=4)"
            arg_type = Int
            default = 8
        "--J1"
            help = "set J1"
            arg_type = Float64
            default = 1.0
        "--J2"
            help = "set J2"
            arg_type = Float64
            default = 0.0
    end
    return parse_args(s)
end

function make_hamiltonian(Jxx,Jzz,list_isite1,list_isite2,N,Nint,Nhilbert)
#    listki = ones(Int,(Nint+1)*Nhilbert)
#    loc = ones(Int,(Nint+1)*Nhilbert)
#    elemnt = zeros((Nint+1)*Nhilbert)
#    for k in 0:Nint
#        for i in 0:Nhilbert-1
#            listki[k*Nhilbert+i+1] = i+1
#        end
#    end
    function get_vec!(vecnew::AbstractVector,vec::AbstractVector)
        lngth = length(vec)
        length(vecnew) == lngth || throw(DimensionMismatch())
        for i in 0:Nhilbert-1
            vecnew[i+1] = 0.0
        end
        for k in 0:Nint-1
            isite1 = list_isite1[k+1]
            isite2 = list_isite2[k+1]
            is1 = 1<<isite1
            is2 = 1<<isite2
            is12 = is1 + is2
            wght = 2.0*Jxx[k+1]
            diag = Jzz[k+1]
            for i in 0:Nhilbert-1
                ibit = i & is12
#                loc[Nint*Nhilbert+i+1] = i+1
                if (ibit==0 || ibit==is12)
                    vecnew[i+1] += +diag*vec[i+1]
#                    elemnt[Nint*Nhilbert+i+1] += diag
                else
                    vecnew[i+1] += -diag*vec[i+1]
#                    elemnt[Nint*Nhilbert+i+1] -= diag
                    iexchg = xor(i,is12)
                    vecnew[i+1] += wght*vec[iexchg+1]
#                    elemnt[k*Nhilbert+i+1] = wght
#                    loc[k*Nhilbert+i+1] = iexchg+1
                end
            end
        end
        return vecnew
    end
    (vecnew,vec) -> get_vec!(vecnew,vec)
#    HamCSC = dropzeros(sparse(listki,loc,elemnt))
#    return HamCSC
end

function calc_zcorr(Nhilbert,Ncorr,list_corr_isite1,list_corr_isite2,psi)
    szz = zeros(Float64,Ncorr)
    for k in 0:Ncorr-1 # loop for all bonds for correlations
        isite1 = list_corr_isite1[k+1]
        isite2 = list_corr_isite2[k+1]
        is1 = 1<<isite1
        is2 = 1<<isite2
        is12 = is1 + is2
        corr = 0.0
        for i in 0:Nhilbert-1 # loop for all spin configurations
            ibit = i & is12
            if (ibit==0 || ibit==is12) # if (spin1,spin2) = (00) or (11): factor = +1
                factor = +1.0
            else # if (spin1,spin2) = (01) or (10): factor = -1
                factor = -1.0
            end
            corr += factor*abs(psi[i+1])^2
        end
        szz[k+1] = 0.25 * corr
        if (isite1==isite2)
            szz[k+1] = 0.25
        end
    end
    return szz
end

function calc_xcorr(Nhilbert,Ncorr,list_corr_isite1,list_corr_isite2,psi)
    sxx = zeros(Float64,Ncorr)
    for k in 0:Ncorr-1 # loop for all bonds for correlations
        isite1 = list_corr_isite1[k+1]
        isite2 = list_corr_isite2[k+1]
        is1 = 1<<isite1
        is2 = 1<<isite2
        is12 = is1 + is2
        corr = 0.0
        for i in 0:Nhilbert-1 # loop for all spin configurations with fixed Sz
            ibit = i & is12
            if (ibit==is1 || ibit==is2) # if (spin1,spin2) = (10) or (01)
                iexchg = xor(i,is12)
                corr += real(psi[i+1]*conj(psi[iexchg+1])) # psi[i]: real
            end
        end
        sxx[k+1] = 0.25 * corr
        if (isite1==isite2)
            sxx[k+1] = 0.25
        end
    end
    return sxx
end

function make_lattice(N,J1,J2)
#    Jxx = []
#    Jzz = []
#    list_isite1 = []
#    list_isite2 = []
    Jxx = Float64[]
    Jzz = Float64[]
    list_isite1 = Int[]
    list_isite2 = Int[]
    Nint = 0
    for i in 0:N-1
        site1 = i
        site2 = mod((i+1),N)
        site3 = mod((i+2),N)
#
        push!(list_isite1,site1)
        push!(list_isite2,site2)
        push!(Jxx,J1)
        push!(Jzz,J1)
        Nint += 1
#
        push!(list_isite1,site1)
        push!(list_isite2,site3)
        push!(Jxx,J2)
        push!(Jzz,J2)
        Nint += 1
    end
    return Jxx, Jzz, list_isite1, list_isite2, Nint
end

function main()
    parsed_args = parse_commandline()
#    println("Parsed args:")
#    for (arg,val) in parsed_args
#        println("  $arg  =>  $val")
#    end

#    N = 16 # should be N>=4
#    J1 = 1.0
#    J2 = 0.0

    N = parsed_args["N"]
    J1 = parsed_args["J1"]
    J2 = parsed_args["J2"]
    Nhilbert = 2^N

    println("J1=",J1)
    println("J2=",J2)
    println("N=",N)
    println("Nhilbert=",Nhilbert)

    println("")
    @time Jxx, Jzz, list_isite1, list_isite2, Nint = make_lattice(N,J1,J2)
    println(Jxx)
    println(Jzz)
    println(list_isite1)
    println(list_isite2)
    println("Nint=",Nint)

    println("")
#    @time Ham = make_hamiltonian(Jxx,Jzz,list_isite1,list_isite2,N,Nint,Nhilbert)
    @time get_vec_LM = make_hamiltonian(Jxx,Jzz,list_isite1,list_isite2,N,Nint,Nhilbert)
    @time Ham = LinearMap(get_vec_LM,Nhilbert;ismutating=true,issymmetric=true)
    #println(Ham)
    #HamElm = SparseArrays.sparse(Ham)
    #println(HamElm)

    #@time ene,vec = eigs(Ham, nev=5, which=:LM)
    @time ene,vec = eigs(Ham, nev=5, which=:SR)
    #println("# GS energy:",ene[1])
    println("# energy:",ene[1]," ",ene[2]," ",ene[3]," ",ene[4]," ",ene[5])
    #vec_sgn = sign(argmax(vec[:,1]))
    #println("# GS wave function:")
    #for i in 0:Nhilbert-1
    #    ii = list_1[i+1]
    #    binii = string(ii,base=2,pad=N)
    #    println(i," ",vec[i+1,1]*vec_sgn," ",binii)
    #end

    println("")
    Ncorr = N # number of total correlations
    list_corr_isite1 = [0 for k in 0:Ncorr-1] # site 1
    list_corr_isite2 = [k for k in 0:Ncorr-1] # site 2
    println(list_corr_isite1)
    println(list_corr_isite2)
    psi = vec[:,1] # choose the ground state
    @time szz = calc_zcorr(Nhilbert,Ncorr,list_corr_isite1,list_corr_isite2,psi)
    @time sxx = calc_xcorr(Nhilbert,Ncorr,list_corr_isite1,list_corr_isite2,psi)
    ss = szz+sxx+sxx
    stot2 = N*sum(ss)
    println("# szz:",szz)
    println("# sxx:",sxx)
    println("# ss:",ss)
    println("# stot(stot+1):",stot2)

end

main()


text:
# Quadratic terms (in objective) & constraints

"""
    add_qpterms!(model::Model, qr::Vector, qc::Vector, qv::Vector)
    add_qpterms!(model, H::SparseMatrixCSC{Float64})
    add_qpterms!(model, H::Matrix{Float64})

Add quadratic terms to the obejctive function
"""
function add_qpterms!(model::Model, qr::Vector{Cint}, qc::Vector{Cint}, qv::Vector{Float64})
    nnz = length(qr)
    (nnz == length(qc) == length(qv)) || error("Inconsistent argument dimensions.")
    if nnz > 0
        ret = @xprs_ccall(chgmqobj, Cint, (
            Ptr{Nothing},    # model
            Cint,         # nnz
            Ptr{Cint},    # qrow
            Ptr{Cint},    # qcol
            Ptr{Float64}  # qval
            ),
            model.ptr_model, nnz, qr .- Cint(1), qc .- Cint(1), qv)

        if ret != 0
            throw(XpressError(model))
        end
    end
    nothing
end
function add_qpterms!(model::Model, qr::Vector, qc::Vector, qv::Vector)
    add_qpterms!(model, ivec(qr), ivec(qc), fvec(qv))
end
function add_qpterms!(model, H::SparseMatrixCSC{Float64}) # H must be symmetric
    n = num_vars(model)
    (H.m == n && H.n == n) || error("H must be an n-by-n symmetric matrix.")

    nnz_h = nnz(H)
    qr = Array{Cint}(undef,  nnz_h)
    qc = Array{Cint}(undef,  nnz_h)
    qv = Array{Float64}(undef,  nnz_h)
    k = 0

    colptr::Vector{Int} = H.colptr
    nzval::Vector{Float64} = nonzeros(H)

    for i = 1 : n
        qi::Cint = convert(Cint, i)
        Hrows = rowvals(H)
        for j = colptr[i]:(colptr[i+1]-1)
            qj = convert(Cint, Hrows[j])

            if qi < qj
                k += 1
                qr[k] = qi
                qc[k] = qj
                qv[k] = nzval[j]#*2
            elseif qi == qj
                k += 1
                qr[k] = qi
                qc[k] = qj
                qv[k] = nzval[j]# * 0.5
            end
        end
    end

    add_qpterms!(model, qr[1:k], qc[1:k], qv[1:k])
end
function add_qpterms!(model, H::Matrix{Float64}) # H must be symmetric
    n = num_vars(model)
    size(H) == (n, n) || error("H must be an n-by-n symmetric matrix.")

    nmax = round(Int,n * (n + 1) / 2)
    qr = Array{Cint}(undef,  nmax)
    qc = Array{Cint}(undef,  nmax)
    qv = Array{Float64}(undef,  nmax)
    k::Int = 0

    for i = 1 : n
        qi = convert(Cint, i)

        v = H[i,i]
        if v != 0.
            k += 1
            qr[k] = qi
            qc[k] = qi
            qv[k] = v# * 0.5
        end

        for j = i+1 : n
            v = H[j, i]
            if v != 0.
                k += 1
                qr[k] = qi
                qc[k] = convert(Cint, j)
                qv[k] = v#*2
            end
        end
    end

    add_qpterms!(model, qr[1:k], qc[1:k], qv[1:k])
end

"""
    add_diag_qpterms!(model, H::Vector)
    add_diag_qpterms!(model, hv::Real)

Add quadratic diagonal terms to the objective
"""
function add_diag_qpterms!(model, H::Vector)  # H stores only the diagonal element
    n = num_vars(model)
    n == length(H) || error("Incompatible dimensions.")
    q = [convert(Cint,1):convert(Cint,n)]
    add_qpterms!(model, q, q, fvec(h))
end
function add_diag_qpterms!(model, hv::Real)  # all diagonal elements are H
    n = num_vars(model)
    q = [convert(Cint,1):convert(Cint,n)]
    add_qpterms!(model, q, q, fill(float64(hv), n))
end

"""
    delq!(model::Model)

delete all quadritic terms formthe objective
"""
function delq!(model::Model)

    n = num_vars(model)
    k = n*(n+1)÷2

    qr = zeros(Int, k)
    qc = zeros(Int, k)
    qv = zeros(k)

    cont = 1
    for i in 1:n
        for j in i:n
            qr[cont] = i
            qc[cont] = j
            cont += 1
        end
    end

    add_qpterms!(model, qr, qc, qv)
end

"""
    getq(model::Model)

Get quadratic terms form the objective function
"""
function getq(model::Model)
    U = getq_upper(model)
    Q = (U+U')/2
    return sparse(Q)
end
function getq_upper(model::Model)
    I,J,V = getq_triplets_upper(model::Model)
    n = num_vars(model)
    return sparse(I,J,V,n,n)
end
function getq_triplets_upper(model::Model)
    nnz = num_qnzs(model)
    n = num_vars(model)
    nels = Array{Cint}(undef,  1)
    nels[1] = nnz

    mstart = Array{Cint}(undef,  n+1)
    mclind = Array{Cint}(undef,  nnz)
    dobjval = Array{Float64}(undef,  nnz)

    ret = @xprs_ccall(getmqobj, Cint, (
        Ptr{Nothing},  # model
        Ptr{Cint},  # mstart
        Ptr{Cint},  # mclind
        Ptr{Float64},  # dobjval
        Cint, #size
        Ptr{Cint}, #nels
        Cint, #first
        Cint #last
        ),
        model.ptr_model,mstart,mclind,dobjval,nnz,nels,0,n-Cint(1))

    if ret != 0
        throw(XpressError(model))
    end
    triangle_nnz = convert(Int,nels[1])

    mstart[end] = triangle_nnz
    I = Array{Int}(undef,  triangle_nnz)
    J = Array{Int}(undef,  triangle_nnz)
    V = Array{Float64}(undef,  triangle_nnz)
    for i in 1:length(mstart)-1
        for j in (mstart[i]+1):mstart[i+1]
            I[j] = i
            J[j] = mclind[j]+1
            V[j] = dobjval[j]
        end
    end
    return I, J, V
end

"""
    add_qconstr!(model::Model, lind::Vector, lval::Vector, qr::Vector, qc::Vector,
    qv::Vector{Float64}, rel::GChars, rhs::Real)

Add quadratic contraint to the model
"""
function add_qconstr!(model::Model, lind::IVec, lval::FVec, qr::IVec, qc::IVec, qv::FVec, rel::Cchar, rhs::Float64)
    # in XPRESS quadratic matrices are added over existing linear constraints
    # ------------------------

    qnnz = length(qr)
    qnnz == length(qc) == length(qv) || error("Inconsistent argument dimensions.")

    lnnz = length(lind)
    lnnz == length(lval) || error("Inconsistent argument dimensions.")

    add_constr!(model, lind, lval, rel, rhs)
    m = num_constrs(model)

    if qnnz > 0
        ret = @xprs_ccall(addqmatrix, Cint, (
            Ptr{Nothing},    # model
            Cint,         # lin contraint o add matrix
            Cint,         # qnnz
            Ptr{Cint},    # qrow
            Ptr{Cint},    # qcol
            Ptr{Float64}  # qval
            ),
            model.ptr_model, m - Cint(1), qnnz, qr .- Cint(1), qc .- Cint(1), qv)

        if ret != 0
            throw(XpressError(model))
        end
    end
    nothing
end
function add_qconstr!(model::Model, lind::Vector, lval::Vector, qr::Vector, qc::Vector,
    qv::Vector{Float64}, rel::GChars, rhs::Real)

    add_qconstr!(model, ivec(lind), fvec(lval), ivec(qr), ivec(qc), fvec(qv), cchar(rel), Float64(rhs))
end

"""
    get_qrows(model::Model)

get indices of quadratic rows
"""
function get_qrows(model::Model)
# int XPRS_CC XPRSgetqrows(XPRSprob prob, int * qmn, int qcrows[]);

    qmn = num_qconstrs(model)

    if qmn > 0

        qcrows = Array{Cint}(undef,  qmn)

        ret = @xprs_ccall(getqrows, Cint, (
            Ptr{Nothing},    # model
            Ptr{Cint},    # qmn
            Ptr{Cint}     # qcrows
            ),
            model.ptr_model, C_NULL, qcrows)

        if ret != 0
            throw(XpressError(model))
        end

        # for i in 1:length(qcrows)
        #     qcrows += 1
        # end
        return qcrows .+ 1

    end

    return Cint[]

end

"""
    get_lrows(model::Model)

get indices of purely linear rows
"""
function get_lrows(model::Model)

    nrows = num_constrs(model)
    nlrows = num_linconstrs(model)
    nqrows = num_qconstrs(model)

    if nlrows == 0
        return Int[]
    end
    if nqrows == 0
        return collect(1:nlrows)
    end

    return setdiff(collect(1:nrows), get_qrows(model))
end

"""
    get_qrowmatrix(model::Model, row::Int)

Get matrix from quadratic row in CSC form
"""
function get_qrowmatrix(model::Model, row::Int)
    U = get_qrowmatrix_upper(model, row)
    Q = (U+U')/2
    return sparse(Q)
end
function get_qrowmatrix_upper(model::Model, row::Int)
    I,J,V = get_qrowmatrix_triplets_upper(model::Model, row::Int)
    n = num_vars(model)
    return sparse(I,J,V,n,n)
end
function get_qrowmatrix_triplets_upper(model::Model, row::Int)
#int XPRS_CC XPRSgetqrowqmatrixtriplets(XPRSprob prob, int irow, int *nqelem, int mqcol1[], int mqcol2[], double dqe[]);
    qrows = get_qrows(model)
    if row in qrows

        nqelem = Array{Cint}(undef, 1)
        ret = @xprs_ccall(getqrowqmatrixtriplets, Cint, (
            Ptr{Nothing},    # model
            Cint,
            Ptr{Cint},    # nqelem
            Ptr{Cint},    # mqcol1
            Ptr{Cint},    # mqcol2
            Ptr{Float64}  # dqe
            ),
            model.ptr_model, Cint(row-Cint(1)), nqelem, C_NULL, C_NULL, C_NULL)

        if ret != 0
            throw(XpressError(model))
        end

        mqcol1 = Array{Cint}(undef,  nqelem[1])
        mqcol2 = Array{Cint}(undef,  nqelem[1])
        dqe = Array{Float64}(undef,  nqelem[1])

        ret = @xprs_ccall(getqrowqmatrixtriplets, Cint, (
            Ptr{Nothing},    # model
            Cint,
            Ptr{Cint},    # nqelem
            Ptr{Cint},    # mqcol1
            Ptr{Cint},    # mqcol2
            Ptr{Float64}  # dqe
            ),
            model.ptr_model, Cint(row-Cint(1)), nqelem, mqcol1, mqcol2, dqe)

        if ret != 0
            throw(XpressError(model))
        end
        return mqcol1.+1, mqcol2.+1, dqe
    end
    return Cint[], Cint[], Float64[]
end


text:
using Pkg

Pkg.add("CSV")
Pkg.add("PDMats")
Pkg.add("Distributions")
Pkg.add("Random")
Pkg.add("StatsBase")
Pkg.add("Plots")
Pkg.add("Polynomials")
Pkg.add("MATLAB")
Pkg.add("DataFrames")
Pkg.add("SpecialFunctions")
Pkg.add("TimerOutputs")
Pkg.add("Parameters")

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
include("update.jl")

function randomwalk_generator(c, prob=0.5)
    x = 0
    put!(c, x)
    while true
        x = update(x, prob)
        put!(c, x)
    end
end

function main(nsteps, prob)
    c = Channel{Int}(32)
    @async randomwalk_generator(c, prob)
    for (i,x) in enumerate(Iterators.take(c, nsteps+1))
        println(i-1, " ", x)
    end
end
if abspath(PROGRAM_FILE) == @__FILE__
    nsteps = length(ARGS)>0 ? parse(Int, ARGS[1]) : 100
    prob = length(ARGS)>1 ? parse(Float64, ARGS[1]) : 0.5
    main(nsteps, prob)
end

max_stars_repo_name: JuliaTagBot/JuliaBook-Samples

max_issues_repo_name: JuliaTagBot/JuliaBook-Samples

max_forks_repo_name: JuliaTagBot/JuliaBook-Samples


text:
using StateSpaceReconstruction
using PerronFrobenius
using Test
using DelayEmbeddings
using StaticArrays
Vpts = [rand(4) for i = 1:50]
Dpts = Dataset([rand(4) for i = 1:50])
Spts = Dataset([rand(4) for i = 1:50])
Mpts = Dataset([rand(4) for i = 1:50])
pts = rand(50, 4)
Tpts = rand(4, 50)

@testset "RectangularInvariantMeasure, constructor" begin

	# On raw points
	@test typeof(rectangularinvariantmeasure(pts, RectangularBinning(3))) <: AbstractRectangularInvariantMeasure
	@test typeof(rectangularinvariantmeasure(pts, RectangularBinning(0.4))) <: AbstractRectangularInvariantMeasure
	@test typeof(rectangularinvariantmeasure(pts, RectangularBinning([3, 2, 2, 3]))) <: AbstractRectangularInvariantMeasure
	@test typeof(rectangularinvariantmeasure(pts, RectangularBinning([0.6, 0.5, 0.6, 0.5]))) <: AbstractRectangularInvariantMeasure

	# On embeddings
	@test typeof(rectangularinvariantmeasure(pts, RectangularBinning(3))) <: AbstractRectangularInvariantMeasure
	@test typeof(rectangularinvariantmeasure(pts, RectangularBinning(0.4))) <: AbstractRectangularInvariantMeasure
	@test typeof(rectangularinvariantmeasure(pts, RectangularBinning([3, 2, 2, 3]))) <: AbstractRectangularInvariantMeasure
	@test typeof(rectangularinvariantmeasure(pts, RectangularBinning([0.6, 0.5, 0.6, 0.5]))) <: AbstractRectangularInvariantMeasure
end

@testset "RectangularInvariantMeasure, different inputs" begin
	b = RectangularBinning(3)
	@test typeof(rectangularinvariantmeasure(Vpts, b)) <: AbstractRectangularInvariantMeasure
	@test typeof(rectangularinvariantmeasure(Dpts, b)) <: AbstractRectangularInvariantMeasure
	@test typeof(rectangularinvariantmeasure(Spts, b)) <: AbstractRectangularInvariantMeasure
	@test typeof(rectangularinvariantmeasure(Mpts, b)) <: AbstractRectangularInvariantMeasure
	@test typeof(rectangularinvariantmeasure(pts, b)) <: AbstractRectangularInvariantMeasure
	@test typeof(rectangularinvariantmeasure(Tpts, b)) <: AbstractRectangularInvariantMeasure
end


text:
@testset "Motzkin with $(factory.constructor)" for factory in sdp_factories
    @polyvar x y

    m = SOSModel(factory)

    p = x^4*y^2 + x^2*y^4 + 1 - 3*x^2*y^2

    @constraint m p >= 0

    JuMP.optimize!(m)

    @test JuMP.dual_status(m) == MOI.InfeasibilityCertificate

    M = SOSModel(factory)

    q = (x^2 + y^2) * p

    @constraint M q >= 0

    JuMP.optimize!(M)

    @test JuMP.primal_status(M) == MOI.FeasiblePoint
end


text:
using SeisReconstruction
using SeisProcessing
using LinearAlgebra
using Test

d = SeisLinearEvents()

#randomly decimate 50% of traces

mask = rand(1,size(d,2));
perc=50
mask[(LinearIndices(mask .< perc/100))[findall(mask .< perc/100)]] .= 0;
mask[(LinearIndices(mask .>= perc/100))[findall(mask .>= perc/100)] ] .= 1;

ddec = Array{Float64}(undef,size(d));

for it = 1 : size(d,1)
    ddec[it,:] = d[it:it,:].*mask;
end

dpocs = SeisPOCS(ddec,dt=0.004,fmax=100,Niter=100,p=1.5)

# test that quality factor is greater than 10 Decibels
quality_factor = 10*log10(norm(d[:],2)/norm(dpocs[:]-d[:],2))
println("Quality factor = ",quality_factor)
@test quality_factor > 5


text:
using DataFrames

function distortBinaryPrediction( y::BitArray{1}, distortionFactor::Float64 )
    res          = deepcopy( y );
    indices      = rand( 1:1:size( y, 1 ), round( Int64, distortionFactor * size( y, 1 ) ) );
    res[indices] = .!y[indices];
    return res;
end

function makeupPredictions()
  n    = 100000;
  y    = Bool.( rand( 0:1,n ) );
  yH1  = Float64.( distortBinaryPrediction( y, 0.20 ) );
  yH2  = Float64.( distortBinaryPrediction( y, 0.21 ) );
  yH3  = Float64.( distortBinaryPrediction( y, 0.22 ) );
  yH4  = Float64.( distortBinaryPrediction( y, 0.23 ) );
  yH5  = Float64.( distortBinaryPrediction( y, 0.24 ) );
  yH6  = Float64.( distortBinaryPrediction( y, 0.24 ) );
  yH7  = Float64.( distortBinaryPrediction( y, 0.26 ) );
  yH8  = Float64.( distortBinaryPrediction( y, 0.27 ) );
  yH9  = Float64.( distortBinaryPrediction( y, 0.28 ) );
  yH10 = Float64.( distortBinaryPrediction( y, 0.29 ) );
  yH11 = Float64.( distortBinaryPrediction( y, 0.30 ) );
  yH12 = Float64.( distortBinaryPrediction( y, 0.33 ) );
  yH13 = Float64.( distortBinaryPrediction( y, 0.34 ) );
  yH14 = Float64.( distortBinaryPrediction( y, 0.35 ) );
  yH15 = Float64.( distortBinaryPrediction( y, 0.36 ) );
  yH16 = Float64.( distortBinaryPrediction( y, 0.37 ) );
  y    = Float64.( y );
  limit           = round( Int64, 0.7 * size( y, 1 ) );
  predictions     = hcat( yH1, yH2, yH3, yH4, yH5, yH6, yH7, yH8, yH9, yH10, yH11, yH12, yH13, yH14, yH15, yH16 );
  predTraining    = predictions[1:limit,:];
  predEval        = predictions[(limit+1):end,:];
  predMatTraining = convert( Matrix, predTraining );
  predMatEval     = convert( Matrix, predEval );
  yTraining       = Float64.( y[1:limit] );
  yEval           = Float64.( y[(limit+1):end] );
  errMatTraining  = ( repeat(  Float64.( yTraining ),outer = [1,size(predMatTraining,2)] ) .- predMatTraining ).^2;
  return predMatTraining, predMatEval, yTraining, yEval, errMatTraining;
end


text:
using Test
using Random
using StaticArrays: SVector
Random.seed!(1234)

using ClimateMachine.VariableTemplates: @vars, varsize, Vars
using ClimateMachine.BalanceLaws
const BL = BalanceLaws
import ClimateMachine.BalanceLaws: vars_state, eq_tends, prognostic_vars

struct TestBL <: BalanceLaw end
struct X <: PrognosticVariable end
struct Y <: PrognosticVariable end
struct F1{PV} <: TendencyDef{Flux{FirstOrder}, PV} end
struct F2{PV} <: TendencyDef{Flux{SecondOrder}, PV} end
struct S{PV} <: TendencyDef{Source, PV} end

prognostic_vars(::TestBL) = (X(), Y())
eq_tends(::X, ::TestBL, ::Flux{FirstOrder}) = (F1{X}(),)
eq_tends(::Y, ::TestBL, ::Flux{FirstOrder}) = (F1{Y}(),)
eq_tends(::X, ::TestBL, ::Flux{SecondOrder}) = (F2{X}(),)
eq_tends(::Y, ::TestBL, ::Flux{SecondOrder}) = (F2{Y}(),)
eq_tends(::X, ::TestBL, ::Source) = (S{X}(),)
eq_tends(::Y, ::TestBL, ::Source) = (S{Y}(),)

@testset "BalanceLaws" begin
    bl = TestBL()
    @test prognostic_vars(bl) == (X(), Y())
    @test fluxes(bl, FirstOrder()) == (F1{X}(), F1{Y}())
    @test fluxes(bl, SecondOrder()) == (F2{X}(), F2{Y}())
    @test sources(bl) == (S{X}(), S{Y}())
    show_tendencies(bl)
    show_tendencies(bl; include_params = true)
end

vars_state(bl::TestBL, st::Prognostic, FT) = @vars begin
    ρ::FT
    ρu::SVector{3, FT}
end

vars_state(bl::TestBL, st::Auxiliary, FT) = @vars()

@testset "Prognostic-Primitive conversion (identity)" begin
    FT = Float64
    bl = TestBL()
    vs_prog = vars_state(bl, Prognostic(), FT)
    vs_prim = vars_state(bl, Primitive(), FT)
    vs_aux = vars_state(bl, Auxiliary(), FT)
    prim_arr = zeros(varsize(vs_prim))
    prog_arr = zeros(varsize(vs_prog))
    aux_arr = zeros(varsize(vs_aux))

    # Test prognostic_to_primitive! identity
    prog_arr .= rand(varsize(vs_prog))
    prog_0 = deepcopy(prog_arr)
    prim_arr .= 0
    BL.prognostic_to_primitive!(bl, prim_arr, prog_arr, aux_arr)
    BL.primitive_to_prognostic!(bl, prog_arr, prim_arr, aux_arr)
    @test all(prog_arr .≈ prog_0)

    # Test primitive_to_prognostic! identity
    prim_arr .= rand(varsize(vs_prim))
    prim_0 = deepcopy(prim_arr)
    prog_arr .= 0
    BL.primitive_to_prognostic!(bl, prog_arr, prim_arr, aux_arr)
    BL.prognostic_to_primitive!(bl, prim_arr, prog_arr, aux_arr)
    @test all(prim_arr .≈ prim_0)
end

max_stars_repo_name: jm-c/CLIMA

max_issues_repo_name: jm-c/CLIMA

max_forks_repo_name: jm-c/CLIMA


text:
# Model
include_model("biped")
model_fixed_time = model
model = free_time_model(no_slip_model(model))

# Visualize
# - Pkg.add any external deps from visualize.jl
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

# Horizon
T = 51
Tm = 26

# Time step
tf = 1.0
h = tf / (T - 1)
t = range(0.0, stop = tf, length = T)

# Configurations
# 1: x pos
# 2: z pos
# 3: torso angle (rel. to downward vertical)
# 4: thigh 1 angle (rel. to downward vertical)
# 5: calf 1 (rel. to thigh 1)
# 6: thigh 2 (rel. to downward vertical)
# 7: calf 2 (rel. to thigh 2)

# θ = pi / 10.0
# q1 = initial_configuration(model, θ) # generate initial config from θ
# k1 = kinematics_2(model, q1, body = :leg_1, mode = :ee)[1]
# k2 = kinematics_2(model, q1, body = :leg_2, mode = :ee)[1]
# d1 = abs(k1 - k2)
# qT = copy(q1)
# qT[1] += 2 * d1
# q1[3] -= pi / 30.0
# q1[4] += pi / 20.0
# q1[5] -= pi / 10.0
# q1, qM = loop_configurations(model, θ)
# qT[1] += 1.0

# θ = pi / 5.0
# q1 = zeros(model.nq)
# q1[3] = pi
# q1[4] = θ #+ pi / 20.0
# q1[5] = -2.0 * θ
# q1[6] = θ #- pi / 20.0
# q1[7] = -2.0 * θ
# q1[2] = model.l2 * cos(θ) + model.l3 * cos(θ)
# qT = copy(q1)
# qT[1] = 0.5
# # qT[2] = 1.0
# kinematics_2(model, q1, body = :leg_1, mode = :ee)[2]
# kinematics_2(model, q1, body = :leg_2, mode = :ee)[2]
function initial_configuration_1(model, θ_torso, θ_thigh_1, θ_leg_1)
    q1 = zeros(model.nq)
    q1[3] = θ_torso #pi - pi / 50.0
    q1[4] = θ_thigh_1 #pi / 7.5
    q1[5] = θ_leg_1 #- pi / 10.0
    z1 = model.l2 * cos(q1[4]) + model.l3 * cos(q1[4] + q1[5])

    q1[6] = - pi / 20.0
    q1[7] = -1.0 * acos((z1 - model.l4 * cos(q1[6])) / model.l5) - q1[6] #-pi / 20.0
    q1[2] = z1

    p1 = kinematics_2(model, q1, body = :leg_1, mode = :ee)
    p2 = kinematics_2(model, q1, body = :leg_2, mode = :ee)
    @show stride = abs(p1[1] - p2[1])

    q1[1] = -1.0 * p1[1]

    qM = copy(q1)
    qM[4] = q1[6]
    qM[5] = q1[7]
    qM[6] = q1[4]
    qM[7] = q1[5]
    qM[1] = abs(p2[1])

    pM_1 = kinematics_2(model, qM, body = :leg_1, mode = :ee)
    pM_2 = kinematics_2(model, qM, body = :leg_2, mode = :ee)

    qT = copy(q1)
    qT[1] = 2 * stride

    pT_1 = kinematics_2(model, qT, body = :leg_1, mode = :ee)
    pT_2 = kinematics_2(model, qT, body = :leg_2, mode = :ee)

    return q1, qM, qT
end

q1, qM, qT = initial_configuration_1(model, pi - pi / 50.0, pi / 7.5, -pi / 5.0)

visualize!(vis, model, [q1])

foot_2_1 = kinematics_2(model, q1, body = :leg_2, mode = :ee)
foot_2_M = kinematics_2(model, qM, body = :leg_2, mode = :ee)

foot_1_M = kinematics_2(model, qM, body = :leg_1, mode = :ee)
foot_1_T = kinematics_2(model, qT, body = :leg_1, mode = :ee)

zh = 0.05
foot_2_x = [range(foot_2_1[1], stop = foot_2_M[1], length = Tm)...,
    [foot_2_M[1] for t = 1:Tm-1]...]
foot_2_z = sqrt.((zh^2.0) * (1.0 .- ((foot_2_x).^2.0)
    ./ abs(foot_2_1[1])^2.0) .+ 1.0e-8)

foot_1_x = [[foot_1_M[1] for t = 1:Tm-1]...,
    range(foot_1_M[1], stop = foot_1_T[1], length = Tm)...]
foot_1_z = sqrt.((zh^2.0) * (1.0 .- ((foot_1_x .- abs(foot_1_T[1] / 2.0)).^2.0)
    ./ abs(foot_1_T[1] / 2.0)^2.0) .+ 1.0e-8)

using Plots
plot(foot_2_x, foot_2_z)
plot!(foot_1_x, foot_1_z, aspect_ratio = :equal)

plot(t, foot_2_x)
plot!(t, foot_1_x)

plot(t, foot_2_z)
plot!(t, foot_1_z)

# u1 = initial_torque(model_τ, qM, h)

# Bounds

# control
# u = (τ1..4, λ1..2, β1..4, ψ1..2, η1...4, s1)
# ul <= u <= uu
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= Inf
_uu[end] = 2.0 * h
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
_ul[end] = 0.01 * h
ul, uu = control_bounds(model, T, _ul, _uu)

foot_lift2 = zeros(model.nq)
foot_lift2[7] = -pi / 50.0
visualize!(vis, model, [q1 + 1.0 * foot_lift2])

foot_lift1 = zeros(model.nq)
foot_lift1[5] = -pi / 50.0
visualize!(vis, model, [qM + 1.0 * foot_lift1])

u1 = initial_torque(model_fixed_time, q1 + 1.0 * foot_lift2, h,
    tol_r = 1.0e-3, tol_d = 1.0e-3)[model.idx_u]
uM = initial_torque(model_fixed_time, qM + 1.0 * foot_lift1, h,
    tol_r = 1.0e-3, tol_d = 1.0e-3)[model.idx_u]

# # q1_mod = copy(q1)
# # q1_mod[3] = Inf
# _xl = -Inf * ones(model.n)
# _xu = Inf * ones(model.n)

xl, xu = state_bounds(model, T,
    x1 = [Inf * ones(model.nq); q1],
    xT = [Inf * ones(model.nq); qT])

xl[Tm][model.nq .+ (1:model.nq)] = copy(qM)
xu[Tm][model.nq .+ (1:model.nq)] = copy(qM)

for t = 1:T
    xl[t][3] = q1[3] - pi / 20.0
    xl[t][10] = q1[3] - pi / 20.0
    xu[t][3] = q1[3] #+ pi / 20.0
    xu[t][10] = q1[3] #+ pi / 20.0
end

# Objective
include_objective(["velocity", "nonlinear_stage"])

# q_ref = linear_interpolation(q1, qM, T)
q_ref = [linear_interpolation(q1, qM, Tm)...,
    linear_interpolation(qM, qT, Tm)[2:end]...]
x0 = configuration_to_state(q_ref)

vis = Visualizer()
render(vis)
visualize!(vis, model, q_ref, Δt = h)

# penalty on slack variable
obj_penalty = PenaltyObjective(1.0e5, model.m - 1)

# quadratic tracking objective
# Σ (x - xref)' Q (x - x_ref) + (u - u_ref)' R (u - u_ref)
obj_control = quadratic_time_tracking_objective(
    [Diagonal(1.0e-5 * ones(model.n)) for t = 1:T],
    [Diagonal([1.0e-3 * ones(model.nu)..., 0.0e-5 * ones(model.m - model.nu - 1)..., 0.0]) for t = 1:T-1],
    [x0[t] for t = 1:T],
    [zeros(model.m) for t = 1:T-1],
    1.0)

# quadratic velocity penalty
# Σ v' Q v
q_v = 100.0 * ones(model.nq)
# q_v[2] = 0.0
obj_velocity = velocity_objective(
    [Diagonal(q_v) for t = 1:T-1],
    model.nq,
    h = h,
    idx_angle = collect([3, 4, 5, 6, 7]))

# torso height
t_h = kinematics_1(model, q1, body = :torso, mode = :com)[2]
function l_stage_torso_h(x, u, t)
    return 1000.0 * (kinematics_1(model,
            view(x, 8:14), body = :torso, mode = :com)[2] - t_h)^2.0
end

l_terminal_torso_h(x) = 0.0
obj_th = nonlinear_stage_objective(l_stage_torso_h, l_terminal_torso_h)

# torso lateral
function l_stage_torso_lat(x, u, t)
    return (0.0 * (kinematics_1(model,
        view(x, 8:14), body = :torso, mode = :com)[1]
        - kinematics_1(model,
        view(x0[t], 8:14), body = :torso, mode = :com)[1])^2.0)
end
l_terminal_torso_lat(x) = 0.0
obj_tl = nonlinear_stage_objective(l_stage_torso_lat, l_terminal_torso_lat)

# foot 1 height
function l_stage_fh1(x, u, t)
    return (10000.0 * (kinematics_2(model,
        view(x, 1:7), body = :leg_1, mode = :ee)[2] - foot_1_z[t])^2.0
        + 10000.0 * (kinematics_2(model,
            view(x, 8:14), body = :leg_1, mode = :ee)[2] - foot_1_z[t])^2.0)
end
l_terminal_fh1(x) = 0.0
obj_fh1 = nonlinear_stage_objective(l_stage_fh1, l_terminal_fh1)

# foot 1 lateral
function l_stage_fl1(x, u, t)
    return (10.0 * (kinematics_2(model,
        view(x, 1:7), body = :leg_1, mode = :ee)[1] - foot_1_x[t])^2.0
        + 10.0 * (kinematics_2(model,
            view(x, 8:14), body = :leg_1, mode = :ee)[1] - foot_1_x[t])^2.0)
end
l_terminal_fl1(x) = 0.0
obj_fl1 = nonlinear_stage_objective(l_stage_fl1, l_terminal_fl1)

# foot 2 height
function l_stage_fh2(x, u, t)
    return (10000.0 * (kinematics_2(model,
        view(x, 1:7), body = :leg_2, mode = :ee)[2] - foot_2_z[t])^2.0
        + 10000.0 * (kinematics_2(model,
            view(x, 8:14), body = :leg_2, mode = :ee)[2] - foot_2_z[t])^2.0)
end
l_terminal_fh2(x) = 0.0
obj_fh2 = nonlinear_stage_objective(l_stage_fh2, l_terminal_fh2)

# foot 2 lateral
function l_stage_fl2(x, u, t)
    return (10.0 * (kinematics_2(model,
        view(x, 1:7), body = :leg_2, mode = :ee)[2] - foot_2_x[t])^2.0
        + 10.0 * (kinematics_2(model,
            view(x, 8:14), body = :leg_2, mode = :ee)[2] - foot_2_x[t])^2.0)
end
l_terminal_fl2(x) = 0.0
obj_fl2 = nonlinear_stage_objective(l_stage_fl2, l_terminal_fl2)

obj = MultiObjective([obj_penalty,
                      obj_control,
                      obj_velocity,
                      obj_th,
                      obj_tl,
                      obj_fh1,
                      obj_fh2,
                      obj_fl1,
                      obj_fl2])

# Constraints
include_constraints(["contact_no_slip", "loop", "free_time"])
con_loop = loop_constraints(model, collect([(2:7)...,(9:14)...]), 1, T)
con_contact = contact_no_slip_constraints(model, T)
con_free_time = free_time_constraints(T)
con = multiple_constraints([con_contact, con_free_time, con_loop])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               # h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# trajectory initialization
u0 = [[1.0e-4 * rand(model.m - 1); h] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

# Solve
optimize = true

if optimize
    include_snopt()

	@time z̄ , info = solve(prob, copy(z0),
		nlp = :SNOPT7,
		tol = 1.0e-3, c_tol = 1.0e-3, mapl = 5,
		time_limit = 60 * 1)
	@show check_slack(z̄, prob)
	x̄, ū = unpack(z̄, prob)
    tf, t, h̄ = get_time(ū)

	# projection
	# Q = [Diagonal(ones(model.n)) for t = 1:T]
	# R = [Diagonal(0.1 * ones(model.m)) for t = 1:T-1]
	# x_proj, u_proj = lqr_projection(model, x̄, ū, h̄[1], Q, R)
    #
	# @show tf
	# @show h̄[1]
	# @save joinpath(@__DIR__, "biped_gait_no_slip.jld2") x̄ ū h̄ x_proj u_proj
else
	# @load joinpath(@__DIR__, "biped_gait_no_slip.jld2") x̄ ū h̄ x_proj u_proj
end

# Visualize
# vis = Visualizer()
# render(vis)
visualize!(vis, model, state_to_configuration(x̄), Δt = h̄[1])

fh1 = [kinematics_2(model,
    state_to_configuration(x_proj)[t], body = :leg_1, mode = :ee)[2] for t = 1:T]
fh2 = [kinematics_2(model,
    state_to_configuration(x_proj)[t], body = :leg_2, mode = :ee)[2] for t = 1:T]

# plot(fh1, linetype = :steppost, label = "foot 1")
# plot!(fh2, linetype = :steppost, label = "foot 2")
#
# plot(hcat(ū...)[1:4, :]',
#     linetype = :steppost,
#     label = "",
#     color = :red,
#     width = 2.0)
#
# plot!(hcat(u_proj...)[1:4, :]',
#     linetype = :steppost,
#     label = "",
#     color = :black)
#
# plot(hcat(u_proj...)[5:6, :]',
#     linetype = :steppost,
#     label = "",
#     width = 2.0)
#
# plot(hcat(state_to_configuration(x̄)...)',
#     color = :red,
#     width = 2.0,
#     label = "")
#
# plot!(hcat(state_to_configuration(x_proj)...)',
#     color = :black,
#     label = "")


text:
function plot_geo_data(data_in, filename, settings; solution = Dict())
    io = open(filename, "w")
    println(io, string("<?xml version=",raw"\"","1.0",raw"\""," encoding=",raw"\"","UTF-8",raw"\"","?>")) #
    println(io, string("<kml xmlns=",raw"\"","http://earth.google.com/kml/2.1",raw"\"",">")) #
    println(io, string("<Document>"))
    draw_order = 1
    if haskey(solution, "solution")
        if haskey(solution["solution"],"nw")
            sol = solution["solution"]["nw"]["1"]
        else
            sol = solution["solution"]
        end
    end
    if haskey(data_in, "solution")
        if haskey(data_in, "nw")
            data = data_in["solution"]["nw"]["1"]
        else
            data = data_in["solution"]
        end
    else
        if haskey(data_in, "nw")
            data = data_in["nw"]["1"]
        else
            data = data_in
        end
    end
    if settings["add_nodes"] == true
        for (b, bus) in data["bus"]
            plot_bus(io, bus, b)
        end
    end
    for (b, branch) in data["branch"]
        plot_branch(io, branch, b, data; color_in = "blue")
    end
    if haskey(data, "ne_branch")
        for (b, branch) in data["ne_branch"]
            if haskey(settings, "plot_solution_only")
                if  sol["ne_branch"]["$b"]["built"] == 1
                    plot_branch(io, branch, b, data; color_in = "blue", name = "Candidate Line")
                end
            else
                plot_branch(io, branch, b, data; color_in = "green", name = "Candidate Line")
            end
        end
    end
    if haskey(data, "branchdc")
       for (bdc, branchdc) in data["branchdc"]
            plot_dc_branch(io, branchdc, bdc, data; color_in = "yellow")
       end
    end

    if haskey(data, "branchdc_ne")
        for (bdc, branchdc) in data["branchdc_ne"]
            if haskey(settings, "plot_solution_only")
                if  sol["branchdc_ne"]["$bdc"]["isbuilt"] == 1
                    plot_dc_branch(io, branchdc, bdc, data; color_in = "yellow", name = "Candidate DC Line")
                end
            else
                plot_dc_branch(io, branchdc, bdc, data; color_in = "red", name = "Candidate DC Line")
            end
        end
     end

     if haskey(data, "convdc")
        for (cdc, convdc) in data["convdc"]
             plot_dc_conv(io, convdc, cdc, data; color_in = "yellow")
        end
     end

     if haskey(data, "convdc_ne")
        for (cdc, convdc) in data["convdc_ne"]
            if haskey(settings, "plot_solution_only")
                if  sol["convdc_ne"]["$cdc"]["isbuilt"] == 1
                    plot_dc_conv(io, convdc, cdc, data; color_in = "yellow", name = "Candidate DC converter")
                end
            else
                plot_dc_conv(io, convdc, cdc, data; color_in = "red", name = "Candidate DC converter")
            end
        end
     end

     if haskey(data, "storage")
        for (s, storage) in data["storage"]
             plot_storage(io, storage, s, data; color_in = "yellow")
        end
     end

     if haskey(data, "ne_storage")
        for (s, storage) in data["ne_storage"]
            if haskey(settings, "plot_solution_only")
                if  sol["ne_storage"]["$s"]["isbuilt"] == 1
                    plot_storage(io, storage, s, data; color_in = "yellow", name = "Candidate storage")
                end
            else
                plot_storage(io, storage, s, data; color_in = "red", name = "Candidate storage")
            end
        end
     end

    println(io,  string("</Document>"))
    println(io,  string("</kml>"))
    close(io)
end

function plot_bus(io, bus, b)
    lat = bus["lat"]
    lon = bus["lon"]
    println(io, string("<Placemark> "));
    println(io, string("<name>Node","$b","</name> "));
    #println(io, string("<description>drawOrder=","$draw_order","</description>"));
    println(io, string("<ExtendedData> "));
    println(io, string("<SimpleData name=",raw"\"","Name",raw"\"",">Node 1</SimpleData>"));
    println(io, string("<SimpleData name=",raw"\"","Description",raw"\"","></SimpleData> "));
    println(io, string("<SimpleData name=",raw"\"","Latitude",raw"\"",">","$lat","</SimpleData>"));
    println(io, string("<SimpleData name=",raw"\"","Longitude",raw"\"",">","$lon","</SimpleData>"));
    println(io, string("<SimpleData name=",raw"\"","Icon",raw"\"","></SimpleData> "));
    println(io, string("</ExtendedData>"));
    println(io, string("<Point> "));
    println(io, string("<coordinates>","$lon",",","$lat",",","0","</coordinates>"));
    println(io, string("</Point> "));
    println(io, string("<Style id=",raw"\"","downArrowIcon",raw"\"",">"));
    println(io, string("<IconStyle> "));
    println(io, string("<Icon> "));
    println(io, string("<href>http://maps.google.com/mapfiles/kml/shapes/placemark_circle_highlight.png</href>"));
    println(io, string("</Icon> "));
    println(io, string("</IconStyle> "));
    println(io, string("</Style> "));
    println(io, string("</Placemark> "));
    return io
end

function plot_branch(io, branch, b, data; color_in = "blue", name = "Line")
    println(io, string("<Placemark> "));
    println(io, string("<name>",name,"$b","</name> "))
    println(io, string("<LineString>"))
    println(io, string("<tessellate>1</tessellate>"))
    println(io, string("<coordinates>"))
    fbus = branch["f_bus"]
    tbus = branch["t_bus"]
    fbus_lat = data["bus"]["$fbus"]["lat"]
    fbus_lon = data["bus"]["$fbus"]["lon"]
    tbus_lat = data["bus"]["$tbus"]["lat"]
    tbus_lon = data["bus"]["$tbus"]["lon"]
    println(io, string("$fbus_lon",",","$fbus_lat","0"))
    println(io, string("$tbus_lon",",","$tbus_lat","0"))
    println(io, string("</coordinates>"))
    println(io, string("</LineString>"))
    println(io, string("<Style>"))
    println(io, string("<LineStyle>"))
    if color_in == "green"
        color = "#FF14F000"
    else
        color = "#FFF00014"
    end
    println(io, string("<color>",color,"</color>"))
    #println(io, string("<description>drawOrder=","$draw_order","</description>"))
    println(io, string("<width>3</width>"))
    println(io, string("</LineStyle>"))
    println(io, string("</Style>"))
    println(io, string("</Placemark>"))
    return io
end

function plot_dc_branch(io, branch, b, data; color_in = "yellow", name = "DC Line")
    println(io, string("<Placemark> "));
    println(io, string("<name>",name,"$b","</name> "))
    println(io, string("<LineString>"))
    println(io, string("<tessellate>1</tessellate>"))
    println(io, string("<coordinates>"))
    fbusdc = branch["fbusdc"]
    tbusdc = branch["tbusdc"]
    if haskey(data, "convdc")
        for (c, conv) in data["convdc"]
            if conv["busdc_i"] == fbusdc
                global fbus = conv["busac_i"]
            end
            if conv["busdc_i"] == tbusdc
                global tbus = conv["busac_i"]
            end
        end
    end
    if haskey(data, "convdc_ne")
        for (c, conv) in data["convdc_ne"]
            if conv["busdc_i"] == fbusdc
                global fbus = conv["busac_i"]
            end
            if conv["busdc_i"] == tbusdc
                global tbus = conv["busac_i"]
            end
        end
    end
    fbus_lat = data["bus"]["$fbus"]["lat"]
    fbus_lon = data["bus"]["$fbus"]["lon"]
    tbus_lat = data["bus"]["$tbus"]["lat"]
    tbus_lon = data["bus"]["$tbus"]["lon"]
    println(io, string("$fbus_lon",",","$fbus_lat","0"))
    println(io, string("$tbus_lon",",","$tbus_lat","0"))
    println(io, string("</coordinates>"))
    println(io, string("</LineString>"))
    println(io, string("<Style>"))
    println(io, string("<LineStyle>"))
    if color_in == "red"
        color = "#FF1400FF"
    else
        color = "#FF14F0FF"
    end
    println(io, string("<color>",color,"</color>"))
    #println(io, string("<description>drawOrder=","$draw_order","</description>"))
    println(io, string("<width>3</width>"))
    println(io, string("</LineStyle>"))
    println(io, string("</Style>"))
    println(io, string("</Placemark>"))
    return io
end

function plot_dc_conv(io, conv, c, data; color_in = "yellow", name = "DC Converter")
    println(io, string("<Placemark> "));
    println(io, string("<name>",name,"$c","</name> "))
    println(io, string("<LineString>"))
    println(io, string("<tessellate>1</tessellate>"))
    println(io, string("<coordinates>"))
    bus = conv["busac_i"]
    bus_lat = data["bus"]["$bus"]["lat"]
    bus_lon = data["bus"]["$bus"]["lon"]
    bus_lon1 = bus_lon + 0.05
    bus_lat1 = bus_lat + 0.05
    println(io, string("$bus_lon1",",","$bus_lat1","0"))
    bus_lon1 = bus_lon + 0.05
    bus_lat1 = bus_lat - 0.05
    println(io, string("$bus_lon1",",","$bus_lat1","0"))
    bus_lon1 = bus_lon - 0.05
    bus_lat1 = bus_lat - 0.05
    println(io, string("$bus_lon1",",","$bus_lat1","0"))
    bus_lon1 = bus_lon - 0.05
    bus_lat1 = bus_lat + 0.05
    println(io, string("$bus_lon1",",","$bus_lat1","0"))
    bus_lon1 = bus_lon + 0.05
    bus_lat1 = bus_lat + 0.05
    println(io, string("$bus_lon1",",","$bus_lat1","0"))
    println(io, string("</coordinates>"))
    println(io, string("</LineString>"))
    println(io, string("<Style>"))
    println(io, string("<LineStyle>"))
    if color_in == "red"
        color = "#FF1400FF"
    else
        color = "#FF14F0FF"
    end
    println(io, string("<color>",color,"</color>"))
    #println(io, string("<description>drawOrder=","$draw_order","</description>"))
    println(io, string("<width>3</width>"))
    println(io, string("</LineStyle>"))
    println(io, string("</Style>"))
    println(io, string("</Placemark>"))
    return io
end

function plot_storage(io, storage, s, data; color_in = "yellow", name = "Storage")
    println(io, string("<Placemark> "));
    println(io, string("<name>",name,"$s","</name> "))
    println(io, string("<LineString>"))
    println(io, string("<tessellate>1</tessellate>"))
    println(io, string("<coordinates>"))
    bus = storage["storage_bus"]
    bus_lat = data["bus"]["$bus"]["lat"]
    bus_lon = data["bus"]["$bus"]["lon"]
    bus_lon1 = bus_lon + 0.05
    bus_lat1 = bus_lat + 0.05
    println(io, string("$bus_lon1",",","$bus_lat1","0"))
    bus_lon1 = bus_lon + 0.05
    bus_lat1 = bus_lat - 0.05
    println(io, string("$bus_lon1",",","$bus_lat1","0"))
    bus_lon1 = bus_lon - 0.025
    bus_lat1 = bus_lat - 0.025
    println(io, string("$bus_lon1",",","$bus_lat1","0"))
    bus_lon1 = bus_lon + 0.05
    bus_lat1 = bus_lat + 0.05
    # println(io, string("$bus_lon1",",","$bus_lat1","0"))
    # bus_lon1 = bus_lon + 0.03
    # bus_lat1 = bus_lat + 0.03
    println(io, string("$bus_lon1",",","$bus_lat1","0"))
    println(io, string("</coordinates>"))
    println(io, string("</LineString>"))
    println(io, string("<Style>"))
    println(io, string("<LineStyle>"))
    if color_in == "red"
        color = "#FF1400FF"
    else
        color = "#FF14F0FF"
    end
    println(io, string("<color>",color,"</color>"))
    #println(io, string("<description>drawOrder=","$draw_order","</description>"))
    println(io, string("<width>3</width>"))
    println(io, string("</LineStyle>"))
    println(io, string("</Style>"))
    println(io, string("</Placemark>"))
    return io
end


text:
module TestFunctionsTests

import ..Algopt # ensure we are using the correct Algopt

using Test

using Algopt.TestFunctions:
Rosenbrock

@testset "Algopt.TestFunctions.Rosenbrock" begin
    ros = Rosenbrock()
    @test 0 == ros.f(ros.argmin)
    @test [0, 0] == ros.∇f(ros.argmin)

    ros = Rosenbrock(a=100, b=16)
    @test 0 == ros.f(ros.argmin)
    @test [0, 0] == ros.∇f(ros.argmin)
end

end # module


text:
"""
    abstract type AbstractUInt16RNG <: AbstractRNG
Supertype for an `AbstractRNG` that is fed random UInt16s.
"""
abstract type AbstractUInt16RNG <: AbstractRNG end

"""
    getrand(r::AbstractUInt16RNG, T::Type, n::Integer)

Generate `n` random `T`s using UInt16s from `r`.
"""
function getrand(::AbstractUInt16RNG, ::Type{UInt16}, n::Integer)
    error("NotImplemented: getrand must be implemented for UInt16 as an AbstractUInt16RNG.")
end
getrand(r::AbstractUInt16RNG, ::Type{UInt32}, n::Integer) =
    Vector{UInt32}(getrand(r, UInt16, n)) .<< 16 .+ Vector{UInt32}(getrand(r, UInt16, n))
getrand(r::AbstractUInt16RNG, ::Type{UInt64}, n::Integer) =
    Vector{UInt64}(getrand(r, UInt32, n)) .<< 32 .+ Vector{UInt64}(getrand(r, UInt32, n))
getrand(r::AbstractUInt16RNG, ::Type{UInt128}, n::Integer) =
    Vector{UInt128}(getrand(r, UInt64, n)) .<< 64 .+ Vector{UInt128}(getrand(r, UInt64, n))

function getrand(r::AbstractUInt16RNG, ::Type{UInt8}, n::Integer)
    vals = Vector{UInt8}(undef, n)
    u16s = getrand(r, UInt16, ceil(Int, n / 2))
    for i in 1:length(u16s)-1
        vals[2i-1:2i] = digits(UInt8, u16s[i], base=256, pad=2)
    end
    if n % 2 == 0
        vals[end-1:end] = digits(UInt8, u16s[end], base=256, pad=2)
    else
        vals[end] = u16s[end] % UInt8
    end
    vals
end
function getrand(r::AbstractUInt16RNG, ::Type{Bool}, n::Integer)
    bools = Vector{Bool}(undef, n)
    uint16needed = ceil(Int, n / 16)
    uint16s = getrand(r, UInt16, uint16needed)
    for i in 1:(uint16needed-1)
        bools[16(i-1)+1:16i] = digits(Bool, uint16s[i], base=2, pad=16)
    end
    bools[16*(uint16needed-1)+1:n] = digits(Bool, uint16s[end], base=2, pad=16)[1:n%16]
    bools
end

function rand!(r::AbstractUInt16RNG, A::Array{Bool}, ::Random.SamplerType{Bool})
    A[:] = getrand(r, Bool, length(A))
    A
end

rand(r::AbstractUInt16RNG, ::Random.SamplerType{Bool}) = getrand(r, Bool, 1)[1]

const UInts = (UInt8, UInt16, UInt32, UInt64, UInt128)
const Ints = (Int8, Int16, Int32, Int64, Int128)

for T in UInts
    @eval function rand!(r::AbstractUInt16RNG, A::Array{$T}, ::Random.SamplerType{$T})
        A[:] = getrand(r, $T, length(A))
        A
    end
    @eval rand(r::AbstractUInt16RNG, ::Random.SamplerType{$T}) = getrand(r, $T, 1)[1]
end
for (T, U) in zip(Ints, UInts)
    @eval function rand!(r::AbstractUInt16RNG, A::Array{$T}, ::Random.SamplerType{$T})
        A[:] = reinterpret.($T, getrand(r, $U, length(A)))
        A
    end
    @eval rand(r::AbstractUInt16RNG, ::Random.SamplerType{$T}) = reinterpret.($T, getrand(r, $U, 1))[1]
end

Random.rng_native_52(::AbstractUInt16RNG) = Float64

const F64BITMASK1 = 0b0011111111111111111111111111111111111111111111111111111111111111
const F64BITMASK2 = 0b0011111111110000000000000000000000000000000000000000000000000000

const F32BITMASK1 = 0b00111111111111111111111111111111
const F32BITMASK2 = 0b00111111100000000000000000000000

function rand!(r::AbstractUInt16RNG, A::Array{Float64}, ::Random.SamplerTrivial{Random.CloseOpen01{Float64}, Float64})
    A[:] = reinterpret.(Float64, getrand(r, UInt64, length(A)) .& F64BITMASK1 .| F64BITMASK2) .- 1
    A
end
function rand(r::AbstractUInt16RNG, ::Random.CloseOpen12{Float64})
    reinterpret(Float64, getrand(r, UInt64, 1)[1] & F64BITMASK1 | F64BITMASK2)
end

function rand!(r::AbstractUInt16RNG, A::Array{Float32}, ::Random.SamplerTrivial{Random.CloseOpen01{Float32}, Float32})
    A[:] = reinterpret.(Float32, getrand(r, UInt32, length(A)) .& F32BITMASK1 .| F32BITMASK2) .- 1
    A
end

# Online and cached UInt16 sources

mutable struct OnlineUInt16RNG <: AbstractUInt16RNG
    label::String
    batchsize::Integer
    cache::Vector{UInt16}
    pointer::Integer
    fetchrand::Function
    prefillthreshold::Integer
    cacheprefill::Vector{UInt16}
    currentlyprefilling::Bool
    function OnlineUInt16RNG(label::String, fetchrand::Function, batchsize::Integer,
                             prefillthreshold::Integer; prefill=prefillthreshold > 0)
        rng = new(label, batchsize, UInt16[], 1, fetchrand,
                  prefillthreshold, UInt16[], false)
        if prefill
            Threads.@spawn prefillcache(rng)
        end
        rng
    end
end

const term256greentored = [119, 190, 220, 221, 215, 208, 202, 160, 124]
function show(io::IO, r::OnlineUInt16RNG)
    print(io, "$(r.label). ")
    unusedcached = length(r.cache) - r.pointer + 1
    color = if r.prefillthreshold == 0
        [:red, :yellow, :green][round(Int, 3*unusedcached/r.batchsize)]
    elseif unusedcached > r.prefillthreshold
        :green
    elseif unusedcached > 0.5r.prefillthreshold
        :yellow
    else
        :red
    end
    printstyled(io, unusedcached; color)
    if r.currentlyprefilling
        printstyled(io, '*'; color=:light_black)
    end
    printstyled(io, '/', r.batchsize; bold=true)
    print(io, " UInt16s cached.")
end

function getrand(r::OnlineUInt16RNG, ::Type{UInt16}, n::Integer)
    if r.batchsize ≤ 0
        r.fetchrand(n)
    else
        ncached = length(r.cache) - r.pointer + 1
        if ncached > n
            r.pointer += n
            if length(r.cacheprefill) == 0 &&
                length(r.cache) - r.pointer + 1 < r.prefillthreshold &&
                ! r.currentlyprefilling
                Threads.@spawn prefillcache(r)
            end
            r.cache[r.pointer-n:r.pointer-1]
        elseif ncached + length(r.cacheprefill) > n
            oldcache = r.cache[r.pointer:end]
            r.cache = r.cacheprefill
            r.cacheprefill = []
            r.pointer = 1
            [oldcache; getrand(r, UInt16, n - length(oldcache))]
        else
            rands = r.fetchrand(n - ncached + r.batchsize)
            res = [r.cache[r.pointer:end]; rands[r.batchsize+1:end]]
            r.cache = rands[1:r.batchsize]
            res
        end
    end
end

function prefillcache(r::OnlineUInt16RNG)
    r.currentlyprefilling = true
    if length(r.cache) == 0
        r.cache = r.fetchrand(r.batchsize)
    else
        r.cacheprefill = r.fetchrand(r.batchsize)
    end
    r.currentlyprefilling = false
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule seatd_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("seatd")
JLLWrappers.@generate_main_file("seatd", UUID("1322d9e5-2885-5d14-b40b-5d6c7dfe9ea1"))
end  # module seatd_jll


text:
"""
    ArrayIndices(iter)

Construct an `ArrayIndices <: AbstractIndices` from an arbitrary Julia iterable with unique
elements. Lookup uses naive iteration.

ArrayIndices make simple and flexible indices using naive search that may be optimal for
small collections. Larger collections are better handled by containers like `Indices`.
"""
struct ArrayIndices{I, Inds <: AbstractArray{I}} <: AbstractIndices{I}
    inds::Inds

    @inline function ArrayIndices{I, Inds}(a::Inds) where {I, Inds <: AbstractArray{I}}
        @boundscheck for i in LinearIndices(a)
            @inbounds x = a[i]
            for j in i+1:last(LinearIndices(a))
                if isequal(x, @inbounds(a[j]))
                    throw(IndexError("Indices are not unique"))
                end
            end
        end
        return new(a)
    end
end

@propagate_inbounds ArrayIndices() = ArrayIndices{Any}([])
@propagate_inbounds ArrayIndices{I}() where {I} = ArrayIndices{I, Vector{I}}(I[])
@propagate_inbounds ArrayIndices{I, Inds}() where {I, Inds} = ArrayIndices{I}(Inds())

@propagate_inbounds ArrayIndices(iter) = ArrayIndices(collect(iter))
@propagate_inbounds ArrayIndices{I}(iter) where {I} = ArrayIndices{I}(collect(I, iter))

@propagate_inbounds ArrayIndices(a::AbstractArray{I}) where {I} = ArrayIndices{I}(a)
@propagate_inbounds ArrayIndices{I}(a::AbstractArray{I}) where {I} = ArrayIndices{I, typeof(a)}(a)

Base.convert(::Type{AbstractIndices{I}}, inds::ArrayIndices) where {I} = convert(ArrayIndices{I}, inds)
Base.convert(::Type{ArrayIndices}, inds::AbstractIndices{I}) where {I} = convert(ArrayIndices{I}, inds)
Base.convert(::Type{ArrayIndices{I}}, inds::AbstractIndices) where {I} = convert(ArrayIndices{I, Vector{I}}, inds)
function Base.convert(::Type{ArrayIndices{I, Inds}}, inds::AbstractIndices) where {I, Inds <: AbstractArray{I}}
    a = convert(Inds, collect(I, inds))
    return @inbounds ArrayIndices{I, typeof(a)}(a)
end

Base.convert(::Type{ArrayIndices{I}}, inds::ArrayIndices{I}) where {I} = inds
Base.convert(::Type{ArrayIndices{I}}, inds::ArrayIndices{<:Any,Inds}) where {I, Inds <: AbstractArray{I}} = convert(ArrayIndices{I, Inds}, inds)
Base.convert(::Type{ArrayIndices{I, Inds}}, inds::ArrayIndices{I, Inds}) where {I, Inds <: AbstractArray{I}} = inds
function Base.convert(::Type{ArrayIndices{I, Inds}}, inds::ArrayIndices) where {I, Inds <: AbstractArray{I}}
    a = convert(Inds, parent(inds))
    return @inbounds ArrayIndices{I, Inds}(a)
end

Base.parent(inds::ArrayIndices) = getfield(inds, :inds)

# Basic interface
@propagate_inbounds function Base.iterate(i::ArrayIndices{I}, state...) where {I}
    iterate(parent(i), state...)
end

Base.in(i::I, inds::ArrayIndices{I}) where {I} = in(i, parent(inds))
Base.IteratorSize(::ArrayIndices) = Base.HasLength()
Base.length(inds::ArrayIndices) = length(parent(inds))

istokenizable(i::ArrayIndices) = true
tokentype(::ArrayIndices) = Int
@inline iteratetoken(inds::ArrayIndices, s...) = iterate(LinearIndices(parent(inds)), s...)
@inline function iteratetoken_reverse(inds::ArrayIndices)
    li = LinearIndices(parent(inds))
    if isempty(li)
        return nothing
    else
        t = last(li)
        return (t, t)
    end
end
@inline function iteratetoken_reverse(inds::ArrayIndices, t)
    li = LinearIndices(parent(inds))
    t -= 1
    if t < first(li)
        return nothing
    else
        return (t, t)
    end
end

@inline function gettoken(inds::ArrayIndices{I}, i::I) where {I}
    a = parent(inds)
    @inbounds for x in LinearIndices(a)
        if isequal(i, a[x])
            return (true, convert(Int, x))
        end
    end
    return (false, 0)
end
@propagate_inbounds gettokenvalue(inds::ArrayIndices, x::Int) = parent(inds)[x]

isinsertable(i::ArrayIndices) = true # Need an array trait here...
@inline function gettoken!(inds::ArrayIndices{I}, i::I, values = ()) where {I}
    a = parent(inds)
    @inbounds for x in LinearIndices(a)
        if isequal(i, a[x])
            return (true, convert(Int, x))
        end
    end
    push!(a, i)
    foreach(v -> resize!(v, length(v) + 1), values)
    return (false, last(LinearIndices(a)))
end

@inline function deletetoken!(inds::ArrayIndices, x::Int, values = ())
    deleteat!(parent(inds), x)
    foreach(v -> deleteat!(v, x), values)
    return inds
end

function Base.empty!(inds::ArrayIndices, values = ())
    empty!(parent(inds))
    foreach(empty!, values)
    return inds
end

empty_type(::Type{<:ArrayIndices}, ::Type{I}) where {I} = ArrayIndices{I, Vector{I}}

function Base.copy(inds::ArrayIndices, ::Type{I}) where {I}
    if I === eltype(inds)
        ArrayIndices{I}(copy(parent(inds)))
    else
        ArrayIndices{I}(convert(AbstractArray{I}, parent(inds)))
    end
end

function Base.filter!(pred, inds::ArrayIndices)
    filter!(pred, parent(inds))
    return inds
end

function randtoken(rng::Random.AbstractRNG, inds::ArrayIndices)
    return rand(rng, keys(parent(inds)))
end

function Base.sort!(inds::ArrayIndices; kwargs...)
    sort!(inds.inds; kwargs...)
    return inds
end


text:
@testset "Sonic.jl" begin

# import XL data
[~,dates,returns] = getXLData["hfma_matlab_data.xlsx','CTA Index"]

n = size(returns,1); # data points
y = fVAI[returns]; get VAMI values

figure; #create figure()
sDate = datenum(dates[1]); #set start date for x-axis()
eDate = datenum(dates[end]); #set end date for x-axis()
xData = range(sDate,eDate,length=n)

plot(xData,y); #plot line chart
xlabel!("Date"); #add x label
ylabel!("VAMI [$]"); #add y label

#box off
#ytick = get(gca,"YTick"); #format axes()
#set(gca,"YTickLabel', sprintf('.2f|",ytick))
#datetick("x','yyyy','keeplimits")
#annotation["textarrow',[.845,.845],[.7,.85],'String','$182.20"]; add arrow & text()

#row & text()

# & text()

#row & text()

end