text:
abstract ColorScheme

Base.getindex(scheme::ColorScheme, i::Integer) = getColor(scheme, i)

export
    cgrad

cgrad() = default_gradient()

function cgrad(arg, values = nothing; alpha = nothing, scale = :identity)
    colors = ColorGradient(arg, alpha=alpha).colors
    values = if values != nothing
        values
    elseif scale in (:log, :log10)
        log10(linspace(1,10,30))
    elseif scale == :log2
        log2(linspace(1,2,30))
    elseif scale == :ln
        log(linspace(1,pi,30))
    elseif scale in (:exp, :exp10)
        (exp10(linspace(0,1,30)) - 1) / 9
    else
        linspace(0, 1, length(colors))
    end
    ColorGradient(colors, values)
end

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

getColor(scheme::ColorScheme) = getColor(scheme, 1)
getColorVector(scheme::ColorScheme) = [getColor(scheme)]

colorscheme(scheme::ColorScheme) = scheme
colorscheme(s::AbstractString; kw...) = colorscheme(Symbol(s); kw...)
colorscheme(s::Symbol; kw...) = haskey(_gradients, s) ? ColorGradient(s; kw...) : ColorWrapper(convertColor(s); kw...)
colorscheme{T<:Real}(s::Symbol, vals::AVec{T}; kw...) = ColorGradient(s, vals; kw...)
colorscheme(cs::AVec, vs::AVec; kw...) = ColorGradient(cs, vs; kw...)
colorscheme{T<:Colorant}(cs::AVec{T}; kw...) = ColorGradient(cs; kw...)
colorscheme(f::Function; kw...) = ColorFunction(f; kw...)
colorscheme(v::AVec; kw...) = ColorVector(v; kw...)
colorscheme(m::AMat; kw...) = size(m,1) == 1 ? map(c->colorscheme(c; kw...), m) : [colorscheme(m[:,i]; kw...) for i in 1:size(m,2)]'
colorscheme(c::Colorant; kw...) = ColorWrapper(c; kw...)

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

convertColor(c::AbstractString) = parse(Colorant, c)
convertColor(c::Symbol) = parse(Colorant, string(c))
convertColor(c::Colorant) = c
convertColor(cvec::AbstractVector) = map(convertColor, cvec)
convertColor(c::ColorScheme) = c
convertColor(v::Void) = RGBA(0,0,0,0)
convertColor(b::Bool) = b ? RGBA(0,0,0,1) : RGBA(0,0,0,0)

function convertColor(c, α::Real)
  c = convertColor(c)
  RGBA(RGB(getColor(c)), α)
end
convertColor(cs::AVec, α::Real) = map(c -> convertColor(c, α), cs)
convertColor(c, α::Void) = convertColor(c)

# backup... try to convert
getColor(c) = convertColor(c)

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

function darken(c, v=0.1)
    rgba = convert(RGBA, c)
    r = max(0, min(rgba.r - v, 1))
    g = max(0, min(rgba.g - v, 1))
    b = max(0, min(rgba.b - v, 1))
    RGBA(r,g,b,rgba.alpha)
end
function lighten(c, v=0.3)
    darken(c, -v)
end

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

const _rainbowColors = [colorant"purple", colorant"blue", colorant"green", colorant"orange", colorant"red"]
const _testColors = [colorant"darkblue", colorant"blueviolet",  colorant"darkcyan",colorant"green",
                     darken(colorant"yellow",0.3), colorant"orange", darken(colorant"red",0.2)]

const _gradients = KW(
    :blues        => [colorant"lightblue", colorant"darkblue"],
    :reds         => [colorant"lightpink", colorant"darkred"],
    :greens       => [colorant"lightgreen", colorant"darkgreen"],
    :redsblues    => [colorant"darkred", RGB(0.8,0.85,0.8), colorant"darkblue"],
    :bluesreds    => [colorant"darkblue", RGB(0.8,0.85,0.8), colorant"darkred"],
    :heat         => [colorant"lightyellow", colorant"orange", colorant"darkred"],
    :grays        => [RGB(.95,.95,.95),RGB(.05,.05,.05)],
    :rainbow      => _rainbowColors,
    :lightrainbow => map(lighten, _rainbowColors),
    :darkrainbow  => map(darken, _rainbowColors),
    :darktest     => _testColors,
    :lighttest    => map(c -> lighten(c, 0.3), _testColors),
  )

function register_gradient_colors{C<:Colorant}(name::Symbol, colors::AVec{C})
    _gradients[name] = colors
end

include("color_gradients.jl")

default_gradient() = ColorGradient(:inferno)

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

"Continuous gradient between values.  Wraps a list of bounding colors and the values they represent."
immutable ColorGradient <: ColorScheme
  colors::Vector
  values::Vector

  function ColorGradient{S<:Real}(cs::AVec, vals::AVec{S} = linspace(0, 1, length(cs)); alpha = nothing)
    if length(cs) == length(vals)
      return new(convertColor(cs,alpha), collect(vals))
    end

    # # otherwise interpolate evenly between the minval and maxval
    # minval, maxval = minimum(vals), maximum(vals)
    # vs = Float64[interpolate(minval, maxval, w) for w in linspace(0, 1, length(cs))]
    # new(convertColor(cs,alpha), vs)

    # interpolate the colors for each value
    vals = merge(linspace(0, 1, length(cs)), vals)
    grad = ColorGradient(cs)
    cs = [getColorZ(grad, z) for z in linspace(0, 1, length(vals))]
    new(convertColor(cs, alpha), vals)
  end
end

Base.getindex(cs::ColorGradient, i::Integer) = getColor(cs, i)
Base.getindex(cs::ColorGradient, z::Number) = getColorZ(cs, z)

# create a gradient from a symbol (blues, reds, etc) and vector of boundary values
function ColorGradient{T<:Real}(s::Symbol, vals::AVec{T} = 0:0; kw...)
  haskey(_gradients, s) || error("Invalid gradient symbol.  Choose from: ", sort(collect(keys(_gradients))))
  cs = _gradients[s]
  if vals == 0:0
    vals = linspace(0, 1, length(cs))
  end
  ColorGradient(cs, vals; kw...)
end

# function ColorGradient{T<:Real}(cs::AVec, vals::AVec{T} = linspace(0, 1, length(cs)); kw...)
#   ColorGradient(map(convertColor, cs), vals; kw...)
# end

function ColorGradient(grad::ColorGradient; alpha = nothing)
  ColorGradient(convertColor(grad.colors, alpha), grad.values)
end

# anything else just gets the default gradient
function ColorGradient(cw; alpha=nothing)
    ColorGradient(default_gradient(), alpha=alpha)
end

getColor(gradient::ColorGradient, idx::Int) = gradient.colors[mod1(idx, length(gradient.colors))]

function getColorZ(gradient::ColorGradient, z::Real)
  cs = gradient.colors
  vs = gradient.values
  n = length(cs)
  @assert n > 0 && n == length(vs)

  # can we just return the first color?
  if z <= vs[1] || n == 1
    return cs[1]
  end

  # find the bounding colors and interpolate
  for i in 2:n
    if z <= vs[i]
      return interpolate_rgb(cs[i-1], cs[i], (z - vs[i-1]) / (vs[i] - vs[i-1]))
    end
  end

  # if we get here, return the last color
  cs[end]
end

getColorVector(gradient::ColorGradient) = gradient.colors

# for 0.3
Colors.RGBA(c::Colorant) = RGBA(red(c), green(c), blue(c), alpha(c))
Colors.RGB(c::Colorant) = RGB(red(c), green(c), blue(c))

function interpolate_rgb(c1::Colorant, c2::Colorant, w::Real)
  rgb1 = RGBA(c1)
  rgb2 = RGBA(c2)
  r = interpolate(rgb1.r, rgb2.r, w)
  g = interpolate(rgb1.g, rgb2.g, w)
  b = interpolate(rgb1.b, rgb2.b, w)
  a = interpolate(rgb1.alpha, rgb2.alpha, w)
  RGBA(r, g, b, a)
end

function interpolate(v1::Real, v2::Real, w::Real)
  (1-w) * v1 + w * v2
end

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

"Wraps a function, taking an index and returning a Colorant"
immutable ColorFunction <: ColorScheme
  f::Function
end

getColor(scheme::ColorFunction, idx::Int) = scheme.f(idx)

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

"Wraps a function, taking an z-value and returning a Colorant"
immutable ColorZFunction <: ColorScheme
  f::Function
end

getColorZ(scheme::ColorZFunction, z::Real) = scheme.f(z)

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

"Wraps a vector of colors... may be vector of Symbol/String/Colorant"
immutable ColorVector <: ColorScheme
  v::Vector{Colorant}
  ColorVector(v::AVec; alpha = nothing) = new(convertColor(v,alpha))
end

getColor(scheme::ColorVector, idx::Int) = convertColor(scheme.v[mod1(idx, length(scheme.v))])
getColorVector(scheme::ColorVector) = scheme.v

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

"Wraps a single color"
immutable ColorWrapper <: ColorScheme
  c::RGBA
  ColorWrapper(c::Colorant; alpha = nothing) = new(convertColor(c, alpha))
end

ColorWrapper(s::Symbol; alpha = nothing) = ColorWrapper(convertColor(parse(Colorant, s), alpha))

getColor(scheme::ColorWrapper, idx::Int) = scheme.c
getColorZ(scheme::ColorWrapper, z::Real) = scheme.c
convertColor(c::ColorWrapper, α::Void) = c.c

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

isbackgrounddark(bgcolor::Color) = Lab(bgcolor).l < 0.5

# move closer to lighter/darker depending on background value
function adjustAway(val, bgval, vmin=0., vmax=100.)
  if bgval < 0.5 * (vmax+vmin)
    tmp = max(val, bgval)
    return 0.5 * (tmp + max(tmp, vmax))
  else
    tmp = min(val, bgval)
    return 0.5 * (tmp + min(tmp, vmin))
  end
end

# borrowed from http://stackoverflow.com/a/1855903:
lightnessLevel(c::Colorant) = 0.299 * red(c) + 0.587 * green(c) + 0.114 * blue(c)

isdark(c::Colorant) = lightnessLevel(c) < 0.5
islight(c::Colorant) = !isdark(c)

function convertHexToRGB(h::Unsigned)
  mask = 0x0000FF
  RGB([(x & mask) / 0xFF for x in  (h >> 16, h >> 8, h)]...)
end

# note: I found this list of hex values in a comment by Tatarize here: http://stackoverflow.com/a/12224359
const _masterColorList = [
    0xFFFFFF, 0x000000, 0x0000FF, 0x00FF00, 0xFF0000, 0x01FFFE, 0xFFA6FE, 0xFFDB66, 0x006401, 0x010067,
    0x95003A, 0x007DB5, 0xFF00F6, 0xFFEEE8, 0x774D00, 0x90FB92, 0x0076FF, 0xD5FF00, 0xFF937E, 0x6A826C,
    0xFF029D, 0xFE8900, 0x7A4782, 0x7E2DD2, 0x85A900, 0xFF0056, 0xA42400, 0x00AE7E, 0x683D3B, 0xBDC6FF,
    0x263400, 0xBDD393, 0x00B917, 0x9E008E, 0x001544, 0xC28C9F, 0xFF74A3, 0x01D0FF, 0x004754, 0xE56FFE,
    0x788231, 0x0E4CA1, 0x91D0CB, 0xBE9970, 0x968AE8, 0xBB8800, 0x43002C, 0xDEFF74, 0x00FFC6, 0xFFE502,
    0x620E00, 0x008F9C, 0x98FF52, 0x7544B1, 0xB500FF, 0x00FF78, 0xFF6E41, 0x005F39, 0x6B6882, 0x5FAD4E,
    0xA75740, 0xA5FFD2, 0xFFB167, 0x009BFF, 0xE85EBE
  ]
const _allColors = map(convertHexToRGB, _masterColorList)
const _darkColors = filter(isdark, _allColors)
const _lightColors = filter(islight, _allColors)
const _sortedColorsForDarkBackground = vcat(_lightColors, reverse(_darkColors[2:end]))
const _sortedColorsForLightBackground = vcat(_darkColors, reverse(_lightColors[2:end]))

const _defaultNumColors = 17

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

# Methods to automatically generate gradients for color selection based on
# background color and a short list of seed colors

# here are some magic constants that could be changed if you really want
const _lightness_darkbg = [80.0]
const _lightness_lightbg = [60.0]
const _lch_c_const = [60]

function adjust_lch(color, l, c)
    lch = convert(LCHab, color)
    convert(RGB, LCHab(l, c, lch.h))
end

function lightness_from_background(bgcolor)
  bglight = convert(LCHab, bgcolor).l
  bglight < 50.0 ? _lightness_darkbg[1] : _lightness_lightbg[1]
end

function gradient_from_list(cs)
    zvalues = Plots.get_zvalues(length(cs))
    indices = sortperm(zvalues)
    sorted_colors = map(RGBA, cs[indices])
    sorted_zvalues = zvalues[indices]
    ColorGradient(sorted_colors, sorted_zvalues)
end

function generate_colorgradient(bgcolor = colorant"white";
                               color_bases = color_bases=[colorant"steelblue",colorant"orangered"],
                               lightness = lightness_from_background(bgcolor),
                               chroma = _lch_c_const[1],
                               n = _defaultNumColors)
  seed_colors = vcat(bgcolor, map(c -> adjust_lch(c, lightness, chroma), color_bases))
  colors = distinguishable_colors(n,
      seed_colors,
      lchoices=Float64[lightness],
      cchoices=Float64[chroma],
      hchoices=linspace(0, 340, 20)
    )[2:end]
  gradient_from_list(colors)
end

function get_color_palette(palette, bgcolor::Union{Colorant,ColorWrapper}, numcolors::Integer)
  grad = if palette == :auto
    generate_colorgradient(bgcolor)
  else
    ColorGradient(palette)
  end
  zrng = get_zvalues(numcolors)
  RGBA[getColorZ(grad, z) for z in zrng]
end

function get_color_palette{C<:Colorant}(palette::Vector{C},
            bgcolor::Union{Colorant,ColorWrapper}, numcolors::Integer)
  palette
end

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

function getpctrange(n::Int)
    n > 0 || error()
    n == 1 && return zeros(1)
    zs = [0.0, 1.0]
    for i in 3:n
        sorted = sort(zs)
        diffs = diff(sorted)
        widestj = 0
        widest = 0.0
        for (j,d) in enumerate(diffs)
            if d > widest
                widest = d
                widestj = j
            end
        end
        push!(zs, sorted[widestj] + 0.5 * diffs[widestj])
    end
    zs
end

function get_zvalues(n::Int)
    offsets = getpctrange(ceil(Int,n/4)+1)/4
    offsets = vcat(offsets[1], offsets[3:end])
    zvalues = Float64[]
    for offset in offsets
        append!(zvalues, offset + [0.0, 0.5, 0.25, 0.75])
    end
    vcat(zvalues[1], 1.0, zvalues[2:n-1])
end

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

make255(x) = round(Int, 255 * x)

function webcolor(c::Color)
    @sprintf("rgb(%d, %d, %d)", [make255(f(c)) for f in [red,green,blue]]...)
end
function webcolor(c::TransparentColor)
    @sprintf("rgba(%d, %d, %d, %1.3f)", [make255(f(c)) for f in [red,green,blue]]..., alpha(c))
end
webcolor(cs::ColorScheme) = webcolor(getColor(cs))
webcolor(c) = webcolor(convertColor(c))
webcolor(c, α) = webcolor(convertColor(getColor(c), α))

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

# converts a symbol or string into a colorant (Colors.RGB), and assigns a color automatically
function getSeriesRGBColor(c, sp::Subplot, n::Int)

  if c == :auto
    c = autopick(sp[:color_palette], n)
  end

  # c should now be a subtype of ColorScheme
  colorscheme(c)
end


text:
using Revise, Test, ReferenceTests

push!(LOAD_PATH, pwd()*"/main")
using eqws

tests = [
    "filter_bank.jl",
    "compute.jl"
]

@testset "eqws" begin
for t in tests
    @testset "$t" begin
        include(t)
    end
end
end

ext: jl

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
# Nonequivalent Groups : Linear methods
function ObservableStats(F::NEAT)
    x = F.rawX; v = F.rawV
    μx = mean(x); σx = std(x)
    μv = mean(v); σv = std(v)
    covxv = cov(x, v); corxv = cor(x, v)
    return μx, σx, μv, σv, covxv, corxv
end

function SummaryStats(F::NEAT)
    r(x) = round2(x; digits = 2)
    x = F.rawX; v = F.rawV
    μx = mean(x); σx = std(x); kx = kurtosis(x); sx = skewness(x);
    μv = mean(v); σv = std(v); kv = kurtosis(v); sv = skewness(v);
    df = DataFrame(
        test = ["X", "V"], mean = [μx, μv], sigma = [σx, σv], kurtosis = [kx, kv], skewness = [sx, sv],
        min = [minimum(x), minimum(v)], max = [maximum(x), maximum(v)], N = [length(x), length(v)]
    )
    println(df)
    return df;
end

mutable struct NEATEquateResult <: NEATEquateMethod
    method
    table
    synthetic
    estimates
    data
end

"""
    Tucker(X::NEAT, Y::NEAT; w₁ = length(X.rawX) / (length(X.rawX) + length(Y.rawX)), w₂ = 1.0 - w₁)

Tucker observed score method.
"""
function Tucker(X::NEAT, Y::NEAT; w₁ = length(X.rawX) / (length(X.rawX) + length(Y.rawX)), w₂ = 1.0 - w₁)
    W = w₁ + w₂
    w₁ = w₁ / W; w₂ = w₂ / W
    # test score
    x = X.rawX; xv = X.rawV
    y = Y.rawX; yv = Y.rawV
    # ObservableStats
    μx, σx, μxv, σxv, covxv, corxv = ObservableStats(X)
    μy, σy, μyv, σyv, covyv, coryv = ObservableStats(Y)
    # regression slope
    γ₁ = covxv / σxv^2
    γ₂ = covyv / σyv^2
    # synthetic mean and var
    μsX = μx - w₂*γ₁*(μxv-μyv)
    μsY = μy + w₁*γ₂*(μxv-μyv)
    σ²sX = σx^2 - w₂*γ₁^2*(σxv^2-σyv^2) + w₁*w₂*γ₁^2*(μxv-μyv)^2
    σ²sY = σy^2 + w₁*γ₂^2*(σxv^2-σyv^2) + w₁*w₂*γ₂^2*(μxv-μyv)^2
    # transformation
    slope = sqrt(σ²sY)/sqrt(σ²sX); intercept = μsY - slope*μsX
    lYx = @. X.tableX.scale * slope + intercept
    tbl = DataFrame(scaleX = X.tableX.scale, lYx = lYx)
    return NEATEquateResult(
        :Tucker,
        tbl,
        DataFrame(Group = [1, 2], μ = [μsX, μsY], σ = [sqrt(σ²sX), sqrt(σ²sY)], γ = [γ₁, γ₂], w = [w₁, w₂]),
        (slope = slope, intercept = intercept),
        (X = X, Y = Y)
    )
end
# Nonequivalent Groups : Levine under a classical congeneric model
"""
    LevineTrueScore(X::NEAT, Y::NEAT; common = :external)
Levine true score method applied to observed scores. This algorithm assumes Classical Congeneric model, which dosen't need to estimate true scores.
"""
function LevineTrueScore(X::NEAT, Y::NEAT; common = :external)
    # test score
    x = X.rawX; xv = X.rawV
    y = Y.rawX; yv = Y.rawV
    # ObservableStats
    μx, σx, μxv, σxv, covxv, corxv = ObservableStats(X)
    μy, σy, μyv, σyv, covyv, coryv = ObservableStats(Y)
    # regression slope
    if common == :internal
        γ₁ = σx^2 / covxv
        γ₂ = σy^2 / covyv
    elseif common == :external
        γ₁ = (σx^2 + covxv) / (σxv^2 + covxv)
        γ₂ = (σy^2 + covyv) / (σyv^2 + covyv)
    else
        println("`common` argument can recieve only :internal or :external")
    end
    # doesn't use the synthetic population.
    # transformation
    slope = γ₂/γ₁ ; intercept = μy + γ₂ * (μxv - μyv) - slope*μx
    lYx = @. X.tableX.scale * slope + intercept
    tbl = DataFrame(scaleX = X.tableX.scale, lYx = lYx)
    return NEATEquateResult(
        :LevineTrueScore,
        tbl,
        DataFrame(Group = [1, 2], μ = [μx, μy], σ = [σx, σy], γ = [γ₁, γ₂], w = [missing, missing]),
        (slope = slope, intercept = intercept),
        (X = X, Y = Y)
    )
end

"""

Levine observed score method. This algorithm assumes Classical Congeneric model, which dosen't need to estimate true scores.
"""
function LevineObservedScore(X::NEAT, Y::NEAT; common = :external, w₁ = length(X.rawX) / (length(X.rawX) + length(Y.rawX)), w₂ = 1.0 - w₁)
    W = w₁ + w₂
    w₁ = w₁ / W; w₂ = w₂ / W
    # test score
    x = X.rawX; xv = X.rawV
    y = Y.rawX; yv = Y.rawV
    # ObservableStats
    μx, σx, μxv, σxv, covxv, corxv = ObservableStats(X)
    μy, σy, μyv, σyv, covyv, coryv = ObservableStats(Y)
    # regression slope
    if common == :internal
        γ₁ = σx^2 / covxv
        γ₂ = σy^2 / covyv
    elseif common == :external
        γ₁ = (σx^2 + covxv) / (σxv^2 + covxv)
        γ₂ = (σy^2 + covyv) / (σyv^2 + covyv)
    else
        println("`common` argument can recieve only :internal or :external")
    end
    # synthetic mean and var
    μsX = μx - w₂*γ₁*(μxv-μyv)
    μsY = μy + w₁*γ₂*(μxv-μyv)
    σ²sX = σx^2 - w₂*γ₁^2*(σxv^2-σyv^2) + w₁*w₂*γ₁^2*(μxv-μyv)^2
    σ²sY = σy^2 + w₁*γ₂^2*(σxv^2-σyv^2) + w₁*w₂*γ₂^2*(μxv-μyv)^2
    # transformation
    slope = sqrt(σ²sY)/sqrt(σ²sX); intercept = μsY - slope*μsX
    lYx = @. X.tableX.scale * slope + intercept
    tbl = DataFrame(scaleX = X.tableX.scale, lYx = lYx)
    return NEATEquateResult(
        :LevineObservedScore,
        tbl,
        DataFrame(Group = [1, 2], μ = [μsX, μsY], σ = [sqrt(σ²sX), sqrt(σ²sY)], γ = [γ₁, γ₂], w = [w₁, w₂]),
        (slope = slope, intercept = intercept),
        (X = X, Y = Y)
    )
end

# Nonequivalent Groups : Chained linear Observed Score Equating
"""
    ChainedLinear(X::NEAT, Y::NEAT)

Equate 2 tests in a chain through the common test scores. Now supports only 2 form equating.

# Basic Algorithm

1. put X on the scale of V -call lV(x);
2. put V on the scale of Y - call lY(v);
3. obtain Y-equivalent as lY(x).

"""
function ChainedLinear(X::NEAT, Y::NEAT)
    # test score
    x = X.rawX; xv = X.rawV
    y = Y.rawX; yv = Y.rawV
    # ObservableStats
    μx, σx, μxv, σxv, covxv, corxv = ObservableStats(X)
    μy, σy, μyv, σyv, covyv, coryv = ObservableStats(Y)
    # regression slope
    γ₁ = σx / σxv
    γ₂ = σy / σyv
    # estimate
    slope = γ₂/γ₁
    intercept = μy + γ₂ *(μxv - μyv) - slope * μx
    lYx = @. X.tableX.scale * slope + intercept
    tbl =  DataFrame(scaleX = X.tableX.scale, lYx = lYx)
    return NEATEquateResult(
        :ChainedLinear,
        tbl,
        DataFrame(Group = [1,2], γ = [γ₁, γ₂]),
        (slope = slope, intercept = intercept),
        (X = X, Y = Y)
    )
end

"""
    ChainedMean(X::NEAT, Y::NEAT)
Chained mean equating, which results are identical for Levine's observed score and true score (not implemented) when w₁ = 1.0 (there is no option for the population weights in the mean equating).
"""
function ChainedMean(X::NEAT, Y::NEAT)
    # test score
    x = X.rawX; xv = X.rawV
    y = Y.rawX; yv = Y.rawV
    # ObservableStats
    μx, σx, μxv, σxv, covxv, corxv = ObservableStats(X)
    μy, σy, μyv, σyv, covyv, coryv = ObservableStats(Y)
    γ₂ = σy / σyv
    # estimate
    slope = 1.0
    intercept = μy + γ₂ *(μxv - μyv) - slope * μx
    lYx = @. X.tableX.scale * slope + intercept
    tbl =  DataFrame(scaleX = X.tableX.scale, lYx = lYx)
    return NEATEquateResult(
        :ChainedMean,
        tbl,
        DataFrame(Group = [1,2], γ = [γ₁, γ₂]),
        (slope = slope, intercept = intercept),
        (X = X, Y = Y)
    )
end
# Nonequivalent Groups : Braun-Holland Linear Method
"""
    BraunHolland(X::NEAT, Y::NEAT; w₁ = length(X.rawX) / (length(X.rawX) + length(Y.rawX)), w₂ = 1.0 - w₁)

Conduct Braun-Holland Linear equating, which use means and standard deviations estimated by the synsetic population distributions.

# Arguments

- `X` frequency table generated by `freqtab(X, V)`, `V` is the test score of common item part
- `Y` same as above.
- `w₁`, `w₂` The synsetic population weights, that are constrained to be w₁ + w₂ = 1.
"""
function BraunHolland(X::NEAT, Y::NEAT; w₁ = length(X.rawX) / (length(X.rawX) + length(Y.rawX)), w₂ = 1.0 - w₁)
    # synthetic weight
    W = w₁ + w₂
    w₁ = w₁ / W; w₂ = w₂ / W
    # prior (the weights from common part)
    J = length(X.tableV.freq)
    h₁ = X.tableV.prob
    h₂ = Y.tableV.prob
    # Conditional distribution
    conditionalX = X.marginal ./ sum(X.marginal)
    conditionalY = Y.marginal ./ sum(Y.marginal)
    for v in 1:length(X.tableV.scale)
        conditionalX[:, v] ./= h₁[v]
    end
    for v in 1:length(Y.tableV.scale)
        conditionalY[:, v] ./= h₂[v]
    end
    # Marginal out
    f₂x = zeros(Float64, length(X.tableX.freq))
    g₁y = zeros(Float64, length(Y.tableX.freq))
    for j in 1:length(X.tableX.scale)
        f₂x[j] = conditionalX[j,:]'h₂
    end
    for j in 1:length(Y.tableX.scale)
        g₁y[j] = conditionalY[j,:]'h₁
    end
    # normalize
    f₂x = f₂x / sum(f₂x)
    g₁y = g₁y / sum(g₁y)
    # synthetic population
    fsx = @. w₁ * X.tableX.prob + w₂ * f₂x; fsx = fsx ./ sum(fsx)
    fsy = @. w₁ * g₁y + w₂ * Y.tableX.prob; fsy = fsy ./ sum(fsy)
    # @show f₂x, g₁y
    # synthetic pupulation parameter
    μsx = fsx' * X.tableX.scale
    σsx = sqrt(fsx' * (X.tableX.scale .- μsx) .^2)
    μsy = fsy' * Y.tableX.scale
    σsy = sqrt(fsy' * (Y.tableX.scale .- μsy) .^2)
    # external regression parameter
    slope = σsy / σsx; intercept = μsy - slope * μsx
    lYx = @. X.tableX.scale * slope + intercept
    tbl = DataFrame(scaleX = X.tableX.scale, lYx = lYx)
    return NEATEquateResult(
        :BraunHolland,
        tbl,
        DataFrame(Group = [1, 2], μ = [μsx, μsy], σ = [σsx, σsy], w = [w₁, w₂]),
        (slope = slope, intercept = intercept),
        (X = X, Y = Y)
    )
end

# Nonequivalent Goups : Frequency Estimation
"""
    FrequencyEstimation(X::NEAT, Y::NEAT; w₁ = length(X.rawX) / (length(X.rawX) + length(Y.rawX)), w₂ = 1.0 - w₁)

Conduct frequency estimation equipercentile equating, which assumes, for both form, the conditional distribution of total score given each common part score is the same in both populations.

"""
function FrequencyEstimation(X::NEAT, Y::NEAT; w₁ = length(X.rawX) / (length(X.rawX) + length(Y.rawX)), w₂ = 1.0 - w₁, case = :middle)
    # synthetic weight
    W = w₁ + w₂
    w₁ = w₁ / W; w₂ = w₂ / W
    # prior (the weights from common part)
    J = length(X.tableV.freq)
    h₁ = X.tableV.prob
    h₂ = Y.tableV.prob
    # Conditional distribution
    conditionalX = X.marginal ./ sum(X.marginal)
    conditionalY = Y.marginal ./ sum(Y.marginal)
    for v in 1:length(X.tableV.scale)
        conditionalX[:, v] ./= h₁[v]
    end
    for v in 1:length(Y.tableV.scale)
        conditionalY[:, v] ./= h₂[v]
    end
    # Marginal out
    f₂x = zeros(Float64, length(X.tableX.freq))
    g₁y = zeros(Float64, length(Y.tableX.freq))
    for j in 1:length(X.tableX.scale)
        f₂x[j] = conditionalX[j,:]'h₂
    end
    for j in 1:length(Y.tableX.scale)
        g₁y[j] = conditionalY[j,:]'h₁
    end
    # normalize
    f₂x = f₂x / sum(f₂x)
    g₁y = g₁y / sum(g₁y)
    # synthetic population
    fsx = @. w₁ * X.tableX.freq + w₂ * f₂x
    fsy = @. w₁ * g₁y + w₂ * Y.tableX.freq
    # Equipercentile Equating
    ftX = FreqTab(DataFrame(scale = X.tableX.scale, freq = fsx, cumprob = cumsum(fsx) ./ sum(fsx)),
                  X.rawX, X.intervalX, (; dummy = "dummy"))
    ftY = FreqTab(DataFrame(scale = Y.tableX.scale, freq = fsy, cumprob = cumsum(fsy) ./ sum(fsy)),
                  Y.rawX, Y.intervalX, (; dummy = "dummy"))
    tbl = Equipercentile(ftX, ftY; case = case)
    return NEATEquateResult(
        :FrequencyEstimation,
        tbl.table,
        DataFrame(Group = [1, 2], μ = [mean(fsx), mean(fsy)], σ = [std(fsx), std(fsy)], w = [w₁, w₂]),
        nothing,
        (X = X, Y = Y)
    )
end

# Nonequivalent Groups : Chained Equipercentile Method
"""
    ChainedEquipercentile(X::NEAT, Y::NEAT; case = :middle)

# Algorithm
Chained equipercentile equating is, in short, to conduct equipercentile equating consecutively.
First, find equipercentile relationship for score X to scores V on the first population.
This function is referred to as 𝑒V1(𝑥).
Second, fubd the equipercentile relationship for converting scores on the common items(V) to scores on the form Y based on examinees from the population 2.
Refer to the resulting function as 𝑒Y2(𝑣).
"""
function ChainedEquipercentile(X::NEAT, Y::NEAT; case = :lower)
    # Find equipercentile relationship for score V on the first population
    ftX = freqtab(X.rawX; interval = X.intervalX, scale = X.tableX.scale)
    ftXV = freqtab(X.rawV; interval = X.intervalV, scale = X.tableV.scale)
    ftYV = freqtab(Y.rawV; interval = Y.intervalV, scale = Y.tableV.scale)
    ftY = freqtab(Y.rawX; interval = Y.intervalX, scale = Y.tableX.scale)
    eV₁x = Equipercentile(ftX, ftXV; case = case)
    # Search percentile of score V on scale Y
    eYxu = zeros(Float64, length(eV₁x.table.scaleX))
    eYxl = zeros(Float64, length(eV₁x.table.scaleX))
    for (i,v) in enumerate(eV₁x.table.eYx)
        P = PRF(v, ftYV)
        eYxu[i] = PFu(P, ftY)
        eYxl[i] = PFl(P, ftY)
    end
    if case == :upper
        eYx = eYxu
    elseif case == :lower
        eYx = eYxl
    elseif case == :both
        eYx = string.(eYxu, "_", eYxl)
    elseif case == :middle
        eYx = (eYxu .+ eYxl) ./ 2.0
    end
    tbl = DataFrame(scaleX = eV₁x.table.scaleX, eYx = eYx)
    return NEATEquateResult(
        :ChainedEquipercentile,
        tbl,
        nothing,
        nothing,
        (X = X, Y = Y)
    )
end
#-----------------


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function derivassign(str::String, f, dstr)

n = length(f);
# for i=1:n
#     f[i] = sympy.simplify(f[i]);
# end
ts, fs = sympy.cse(f);

m = length(ts);
for i = 1:m
    str1 = sympy.ccode(ts[i][1]);
    str2 = sympy.ccode(ts[i][2]);
    str = str * "\t\tT " * str1 * " = " * str2 * ";\n";
end

if length(fs)==1
    for i = 1:n
        #str1 = "\t\tu[" * string(i-1) * "*ng+i" * "]";
        str1 = "\t\t" * dstr * "[" * string(i-1) * " + i*" * string(n) * "]";
        str2 = sympy.ccode(fs[1][i]);
        str = str * str1 * " = " * str2 * ";\n";
    end
elseif length(fs)>1
    for i = 1:n
        #str1 = "\t\tu[" * string(i-1) * "*ng+i" * "]";
        str1 = "\t\t"  * dstr * "[" * string(i-1) * " + i*" * string(n) * "]";
        str2 = sympy.ccode(fs[i]);
        str = str * str1 * " = " * str2 * ";\n";
    end
else
    error("Symbolic expressions are wrong");
end

return str;

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# ---
# title: 200. Number of Islands
# id: problem200
# author: Indigo
# date: 2021-01-20
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search, Union Find
# link: <https://leetcode.com/problems/number-of-islands/description/>
# hidden: true
# ---
#
# Given an `m x n` 2d `grid` map of `'1'`s (land) and `'0'`s (water), return
# _the number of islands_.
#
# An **island** is surrounded by water and is formed by connecting adjacent
# lands horizontally or vertically. You may assume all four edges of the grid
# are all surrounded by water.
#
#
#
# **Example 1:**
#
#
#
#     Input: grid = [
#       ["1","1","1","1","0"],
#       ["1","1","0","1","0"],
#       ["1","1","0","0","0"],
#       ["0","0","0","0","0"]
#     ]
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: grid = [
#       ["1","1","0","0","0"],
#       ["1","1","0","0","0"],
#       ["0","0","1","0","0"],
#       ["0","0","0","1","1"]
#     ]
#     Output: 3
#
#
#
#
# **Constraints:**
#
#   * `m == grid.length`
#   * `n == grid[i].length`
#   * `1 <= m, n <= 300`
#   * `grid[i][j]` is `'0'` or `'1'`.
#
#
## @lc code=start
using LeetCode

function lands_num(grid::Vector{Vector{String}})
    nr, nc = length(grid), length(grid[1])
    function dfs(grid::Vector{Vector{String}}, r, c)
        grid[r][c] = "0"
        (r + 1 <= nr) && grid[r+1][c] == "1" && dfs(grid, r+1, c)
        (r - 1 > 0) && grid[r-1][c] == "1" && dfs(grid, r-1, c)
        (c + 1 <= nc) && grid[r][c+1] == "1" && dfs(grid, r, c+1)
        (c - 1 > 0) && grid[r][c-1] == "1" && dfs(grid, r, c-1)
        nothing
    end
    res = 0
    for i in 1:nr
        for j in 1:nc
            if grid[i][j] == "1"
                res += 1
                dfs(grid, i, j)
            end
        end
    end
    res
end
## @lc code=end


text:
# Updated copy of the same code in Base, but with bugs fixed
using Core.Compiler: count_added_node!, add!, NewSSAValue, add_pending!,
    StmtRange, BasicBlock

Base.length(c::Core.Compiler.NewNodeStream) = Core.Compiler.length(c)
Base.setindex!(i::Instruction, args...) = Core.Compiler.setindex!(i, args...)
Core.Compiler.BasicBlock(x::UnitRange) =
    BasicBlock(StmtRange(first(x), last(x)))
Core.Compiler.BasicBlock(x::UnitRange, preds::Vector{Int}, succs::Vector{Int}) =
    BasicBlock(StmtRange(first(x), last(x)), preds, succs)
Base.size(x::Core.Compiler.UnitRange) = Core.Compiler.size(x)

function my_insert_node!(compact::IncrementalCompact, before, inst::NewInstruction, attach_after::Bool=false)
    @assert inst.effect_free_computed
    if isa(before, SSAValue)
        if before.id < compact.result_idx
            count_added_node!(compact, inst.stmt)
            line = something(inst.line, compact.result[before.id][:line])
            node = add!(compact.new_new_nodes, before.id, attach_after)
            node[:inst], node[:type], node[:line], node[:flag] = inst.stmt, inst.type, line, inst.flag
            return NewSSAValue(node.idx)
        else
            line = something(inst.line, compact.ir.stmts[before.id][:line])
            node = add_pending!(compact, before.id, attach_after)
            node[:inst], node[:type], node[:line], node[:flag] = inst.stmt, inst.type, line, inst.flag
            os = OldSSAValue(length(compact.ir.stmts) + length(compact.ir.new_nodes) + length(compact.pending_nodes))
            push!(compact.ssa_rename, os)
            push!(compact.used_ssas, 0)
            return os
        end
    elseif isa(before, OldSSAValue)
        pos = before.id
        if pos < compact.idx
            renamed = compact.ssa_rename[pos]
            count_added_node!(compact, inst.stmt)
            line = something(inst.line, compact.result[renamed.id][:line])
            node = add!(compact.new_new_nodes, renamed.id, attach_after)
            node[:inst], node[:type], node[:line], node[:flag] = inst.stmt, inst.type, line, inst.flag
            return NewSSAValue(node.idx)
        else
            if pos > length(compact.ir.stmts)
                #@assert attach_after
                info = compact.pending_nodes.info[pos - length(compact.ir.stmts) - length(compact.ir.new_nodes)]
                pos, attach_after = info.pos, info.attach_after
            end
            line = something(inst.line, compact.ir.stmts[pos][:line])
            node = add_pending!(compact, pos, attach_after)
            node[:inst], node[:type], node[:line], node[:flag] = inst.stmt, inst.type, line, inst.flag
            os = OldSSAValue(length(compact.ir.stmts) + length(compact.ir.new_nodes) + length(compact.pending_nodes))
            push!(compact.ssa_rename, os)
            push!(compact.used_ssas, 0)
            return os
        end
    elseif isa(before, NewSSAValue)
        before_entry = compact.new_new_nodes.info[before.id]
        line = something(inst.line, compact.new_new_nodes.stmts[before.id][:line])
        new_entry = add!(compact.new_new_nodes, before_entry.pos, attach_after)
        new_entry[:inst], new_entry[:type], new_entry[:line], new_entry[:flag] = inst.stmt, inst.type, line, inst.flag
        return NewSSAValue(new_entry.idx)
    else
        error("Unsupported")
    end
end


text:
@testset "Utils: push!, pushfirst!" begin
    for coarse_quantizer in [:naive, :hnsw]
        ivfadc = build_index_random_data(;coarse_quantizer=coarse_quantizer,
                                         index_type=UInt8)
        ol = length(ivfadc)
        nnv = 256 - nvectors

        # push!
        for i in 1:nnv
            push!(ivfadc, rand(nrows))
        end
        @test length(ivfadc) == ol + nnv
        @test_throws AssertionError push!(ivfadc, rand(nrows))  # index is full
        delete_from_index!(ivfadc, [1])
        @test_throws AssertionError push!(ivfadc, rand(nrows+1))  # wrong dimension

        # pushfirst!
        for i in 1:nnv-1
            delete_from_index!(ivfadc, [i])
        end
        for i in 1:nnv
            pushfirst!(ivfadc, rand(nrows))
        end
        @test length(ivfadc) == ol + nnv
        @test_throws AssertionError pushfirst!(ivfadc, rand(nrows))  # index is full
        delete_from_index!(ivfadc, [1])
        @test_throws AssertionError pushfirst!(ivfadc, rand(nrows+1))  # wrong dimension
    end
end

@testset "Utils: pop!, popfirst!" begin
    for coarse_quantizer in [:naive, :hnsw]
        ivfadc = build_index_random_data(;coarse_quantizer=coarse_quantizer,
                                         index_type=UInt8)
        ol = length(ivfadc)
        npops = 1
        # pop!
        for i in 1:npops
            v = pop!(ivfadc)
            @test v isa Vector{Float64}
            @test length(v) == size(ivfadc, 1)
        end
        @test length(ivfadc) == ol - npops

        # popfirst!
        ol = length(ivfadc)
        for i in 1:npops
            v=popfirst!(ivfadc)
            @test v isa Vector{Float64}
            @test length(v) == size(ivfadc, 1)
        end
        @test length(ivfadc) == ol - npops
    end
end

@testset "Utils: delete_from_index!" begin
    for coarse_quantizer in [:naive, :hnsw]
        ivfadc = build_index_random_data(coarse_quantizer=coarse_quantizer)
        ivfadc_copy = deepcopy(ivfadc)
        n = length(ivfadc)
        L1s, L1e = 1, 5
        L2s, L2e = 10, 30
        L3s, L3e = n-5, n
        first_list= collect(L1s : L1e)
        second_list = collect(L2s : L2e)
        third_list = collect(L3s : L3e)
        idxs_to_delete = vcat(first_list...,
                              second_list...,
                              third_list...)

        # Test that points are deleted
        delete_from_index!(ivfadc, idxs_to_delete)
        @test length(ivfadc) == n - length(idxs_to_delete)

        # Test that indexes are correctly updated
        mismatches = 0
        for (cl, ivlist) in enumerate(ivfadc_copy.inverse_index)
            cluster_indexes = ivlist.idxs
            cluster_indexes_del = ivfadc.inverse_index[cl].idxs  # after deletion
            found = intersect(cluster_indexes, idxs_to_delete .- 1)
            @test length(cluster_indexes) == length(cluster_indexes_del) + length(found)
            for (i, idx) in enumerate(cluster_indexes .+ 1)
                if idx > L1e && idx < L2s
                    shift = L1e - L1s + 1
                elseif idx > L2e && idx < L3s
                    shift = (L1e - L1s + 1) + (L2e - L2s + 1)
                else
                    shift = nothing
                end

                if shift != nothing
                    oldpos = i
                    newval = idx - shift - 1
                    newpos = findfirst(x->x==newval, cluster_indexes_del)
                    if ivfadc_copy.inverse_index[cl].codes[oldpos] !=
                            ivfadc.inverse_index[cl].codes[newpos]
                        mismatches += 1
                    end
                end
            end
        end
        @test mismatches == 0
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module NetworkSecurityGroupClient

using Random
using Dates
using HTTP
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_NetworkSecurityGroupsApi.jl")
include("api_SecurityRulesApi.jl")

# export models
export convert, NetworkSecurityGroupListResult
export convert, NetworkSecurityGroupPropertiesFormat
export convert, Resource
export convert, SecurityRuleListResult
export convert, SecurityRulePropertiesFormat
export convert, SubResource
export convert, NetworkInterface
export convert, NetworkSecurityGroup
export convert, SecurityRule
export convert, Subnet

# export operations
export convert, NetworkSecurityGroupsApi, SecurityRulesApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
using MovingWindow
using Test

@testset "MovingWindow.jl" begin
    # Write your tests here.
end


text:
####
# These are special routines to make operations involving +
# more efficient
####

const Add{Factors<:Tuple} = Applied{<:Any, typeof(+), Factors}

const AddArray{T,N,Factors<:Tuple} = ApplyArray{T,N,typeof(+), Factors}
const AddVector{T,Factors<:Tuple} = AddArray{T,1,Factors}
const AddMatrix{T,Factors<:Tuple} = AddArray{T,2,Factors}

AddArray(factors...) = ApplyArray(+, factors...)

"""
    Add(A1, A2, …, AN)

A lazy representation of `A1 + A2 + … + AN`; i.e., a shorthand for `applied(+, A1, A2, …, AN)`.
"""
Add(As...) = applied(+, As...)

for op in (:+, :-)
    @eval begin
        size(M::Applied{<:Any, typeof($op)}, p::Int) = size(M)[p]
        axes(M::Applied{<:Any, typeof($op)}, p::Int) = axes(M)[p]
        ndims(M::Applied{<:Any, typeof($op)}) = ndims(first(M.args))

        length(M::Applied{<:Any, typeof($op)}) = prod(size(M))
        size(M::Applied{<:Any, typeof($op)}) = length.(axes(M))
        axes(M::Applied{<:Any, typeof($op)}) = axes(first(M.args))

        eltype(M::Applied{<:Any, typeof($op)}) = promote_type(map(eltype,M.args)...)

        combine_mul_styles(::ApplyLayout{typeof($op)}) = IdentityMulStyle()
    end
end

getindex(M::Add, k::Integer) = sum(getindex.(M.args, k))
getindex(M::Add, k::Integer, j::Integer) = sum(getindex.(M.args, k, j))

getindex(M::Add, k::CartesianIndex{1}) = M[convert(Int, k)]
getindex(M::Add, kj::CartesianIndex{2}) = M[kj[1], kj[2]]

for MulAdd_ in [MatMulMatAdd, MatMulVecAdd]
    # `MulAdd{ApplyLayout{typeof(+)}}` cannot "win" against
    # `MatMulMatAdd` and `MatMulVecAdd` hence `@eval`:
    @eval begin
        function materialize!(M::$MulAdd_{ApplyLayout{typeof(+)}})
            α, A, B, β, C = M.α, M.A, M.B, M.β, M.C
            if C ≡ B
                B = copy(B)
            end
            _fill_lmul!(β, C)
            for a in arguments(A)
                C .= applied(+,applied(*,α, a,B), C)
            end
            C
        end
        function materialize!(M::$MulAdd_{ApplyLayout{typeof(-)}})
            α, A, B, β, C = M.α, M.A, M.B, M.β, M.C
            if C ≡ B
                B = copy(B)
            end
            _fill_lmul!(β, C)
            a1,a2 = arguments(A)
            C .= applied(+,applied(*,α, a1,B), C)
            C .= applied(+,applied(*,-α, a2,B), C)
            C
        end
    end
end

###
# views
####
_view_tuple(a, b::Tuple) = view(a, b...)
for op in (:+, :-)
    @eval begin
        sublayout(a::ApplyLayout{typeof($op)}, _) = a
        arguments(::ApplyLayout{typeof($op)}, a::SubArray) =
            _view_tuple.(arguments(parent(a)), Ref(parentindices(a)))
        call(::ApplyLayout{typeof($op)}, a::SubArray) = $op
    end
end


text:
using EzXML
@testset "Variables" begin
    X = parsexml("""
                 <param name="x" type="real" d1="1" d2="1" local="false" dynamics="any" controlled="true" />
                 """)
    vars = Dict()
    SpaceExParser.add_variable!(vars, root(X))
    @test collect(keys(vars)) == [:x]
    @test vars[:x]["local"] == false
    @test vars[:x]["controlled"] == true
    @test vars[:x]["dynamics"] == "any"
    @test vars[:x]["id"] == 1
end


text:
using MathProgComplex, DataStructures, OPFInstances

function main()

    problem_c, pt = import_from_dat(getinstancepath("Matpower", "QCQP", "case9"))
    problem = pb_cplx2real(problem_c)

    relax_ctx = set_relaxation(problem; hierarchykind=:Real,
                                        # symmetries=[PhaseInvariance],
                                        d = 1,
                                        params = Dict(:opt_outlev=>1))

    # println("\n--------------------------------------------------------")
    # println("problem = \n$problem")

    # println("\n--------------------------------------------------------")
    # println("relax_ctx = \n$relax_ctx")

    ########################################
    # Construction du sparsity pattern, extension chordale, cliques maximales.
    max_cliques = get_maxcliques(relax_ctx, problem)
    relax_ctx.relaxparams[:opt_nb_cliques] = length(max_cliques)

    # println("\n--------------------------------------------------------")
    # println("max cliques =\nmax_cliques")

    ########################################
    # Compute moment and localizing matrices parameters: order et variables
    momentmat_param, localizingmat_param = build_sparsity(relax_ctx, problem, max_cliques)

    ########################################
    # Build the moment relaxation problem
    momentrel = build_momentrelaxation(relax_ctx, problem, momentmat_param, localizingmat_param, max_cliques)

    # println("\n--------------------------------------------------------")
    # println("momentrel = $momentrel")

    ########################################
    # Convert to a primal SDP problem
    sosrel = build_SOSrelaxation(relax_ctx, momentrel)

    # println("\n--------------------------------------------------------")
    # println("sosrel = \n$sosrel")

    path = joinpath(pwd(), "Mosek_runs", "worksdp")
    mkpath(path)
    export_SDP(sosrel, path)

    sdp_instance = build_SDP_Instance_from_sdpfiles(path)

    primal = SortedDict{Tuple{String,String,String}, Float64}()
    dual = SortedDict{Tuple{String, String, String}, Float64}()

    primobj, dualobj = solve_mosek(sdp_instance::SDP_Problem, primal, dual, sol_info=relax_ctx.relaxparams)

end

main()


text:
using Documenter, ReinforcementLearning

makedocs(
    modules = [ReinforcementLearning],
    clean = false,
    format = :html,
#    assets = [
#        joinpath("assets", "favicon.ico"),
#        joinpath("assets", "style.css"),
#    ],
    sitename = "Reinforcement Learning",
    authors = "Christof Stocker",
    linkcheck = !("skiplinks" in ARGS),
    pages = Any[
        "Home" => "index.md",
        "Developer Documentation" => Any[
            joinpath("devdocs", "design.md"),
        ],
        hide("Indices" => "indices.md"),
        "LICENSE.md",
    ],
    html_prettyurls = !("local" in ARGS),
)

deploydocs(
    repo = "github.com/Evizero/ReinforcementLearning.jl.git",
    target = "build",
    julia = "0.6",
    deps = nothing,
    make = nothing,
)


text:
using Documenter, MLDatasets

## Commented out since gives warning
# DocMeta.setdocmeta!(MLDatasets, :DocTestSetup, :(using MLDatasets); recursive=true)

# Build documentation.
# ====================

makedocs(
    modules = [MLDatasets],
    doctest = true,
    clean = false,
    sitename = "MLDatasets.jl",
    format = Documenter.HTML(
        canonical = "https://juliadata.github.io/MLDatasets.jl/stable/",
        assets = ["assets/favicon.ico"],
        prettyurls = get(ENV, "CI", nothing) == "true",
        collapselevel=3,
    ),

    authors = "Hiroyuki Shindo, Christof Stocker, Carlo Lucibello",

    pages = Any[
        "Home" => "index.md",
        "Datasets" => Any[
            "Graphs" => "datasets/graphs.md",
            "Miscellaneous" => "datasets/misc.md",
            "Text" => "datasets/text.md",
            "Vision" => "datasets/vision.md",
        ],
        "Utils" => "utils.md",
        "Data Containers" => "containers/overview.md",
        "LICENSE.md",
    ],
    strict = true,
    checkdocs = :exports
)

deploydocs(repo = "github.com/JuliaML/MLDatasets.jl.git")


text:
using GDAXClient
using Base.Test
import GDAXClient: onWSMessage

user = GDAXUser("https://api.gdax.com", "wss://ws-feed.gdax.com", "", "", "")

@testset "Public" begin
    @test Requests.statuscode(getProducts(user)) == 200
    @test Requests.statuscode(getServerTime(user)) == 200
    @test Requests.statuscode(getProductOrderBook(user, "BTC-EUR")) == 200
    @test Requests.statuscode(getProductTicker(user, "BTC-EUR")) == 200
    @test Requests.statuscode(getProductTrades(user, "BTC-EUR")) == 200

    resp = getHistoricRates(user, "BTC-USD", DateTime(2018, 1, 4), DateTime(2018, 1, 4, 0, 2, 0), 60)
    Requests.json(resp)
    @test Requests.statuscode(resp) == 200
    @test length(Requests.json(resp)) == 2

    resp = get24HourStats(user, "BTC-USD")
    @test Requests.statuscode(resp) == 200
    @test !isempty(Requests.json(resp))

    @test Requests.statuscode(getCurrencies(user)) == 200
end

@testset "Orders" begin
    @test Requests.statuscode(listOrders(user)) == 401
    @test Requests.statuscode(placeOrder(user, "buy", "BTC-USD", 0.00010001, 12000.0)) == 401
    @test Requests.statuscode(cancelOrder(user, string(Base.Random.uuid4()))) == 401
end

mutable struct TestEventsHandler <: GDAXClient.AbstractGDAXMessageHandler
    message_counter::Int64
end

function onWSMessage(this::TestEventsHandler, msg::Dict{String, Any})::Void
    this.message_counter += 1
    return nothing
end

@testset "Websocket" begin
    subscription = Dict("type" => "subscribe",
                        "product_ids" => ["BTC-EUR", "BTC-USD"],
                        "channels" => ["heartbeat", "level2", "full"])

    client = GDAXWebSocketClient(user,
                            subscription,
                            TestEventsHandler(0))

    connect(client)
    subscribe(client)
    sleep(10)
    @test client.events_handler.message_counter > 0
    unsubscribe(client)
    sleep(10) #make sure to unsubscribe, though it might fail
    n = client.events_handler.message_counter
    sleep(15) #wait for new message that may arrive
    # test that there is no new messages
    @test n == client.events_handler.message_counter
end

using FIX
using DataStructures
import FIX: onFIXMessage

function onFIXMessage(this::TestEventsHandler, msg::OrderedDict{Int64, String})
    this.message_counter += 1
    return nothing
end

@testset "fix orders" begin
    gadx_keys = JSON.parsefile(joinpath(homedir(),".gdax_keys"))
    api_key = gadx_keys["Key"]
    api_secret = gadx_keys["Secret"]
    passphrase = gadx_keys["Passphrase"]

    reg_user = GDAXUser("https://api.gdax.com", "wss://ws-feed.gdax.com", api_key, api_secret, passphrase)

    handler = TestEventsHandler(0)

    client = GDAXClient.fixconnect(reg_user, handler)
    m, mstr = placeOrder(client, "buy", "BTC-EUR", 0.0001001, 10000.0)
    sleep(2)
    @test !isempty(FIX.getOpenOrders(client))
    cancelAll(client)
    sleep(1)
    @test isempty(FIX.getOpenOrders(client))
    logout(client)
    sleep(1)
    @test !isempty(client.m_messages.incoming.logout)
    close(client)
end

max_stars_repo_name: DatName/gdax-julia

max_issues_repo_name: DatName/gdax-julia

max_forks_repo_name: DatName/gdax-julia


text:
# #
# Reproduces the figures in
# the work of V. Duchêne and B. Mélinand
# on the quadratic pseudo-spectral method (WW2)
# #
export IntegrateWW2,Figure
using WaterWaves1D,FFTW,Plots,LinearAlgebra,ProgressMeter;
#using JLD #(uncomment if using @save)

#--- Integration

"""
	IntegrateWW2(;init,args)

Integrate in time the WW2 system with an initial data depending on the provided `init`
- if `init=1`, then surface deformation `η(t=0,x)=exp(-x^p)` and velocity `v(t=0,x)=0` (with `p` provided as an optional argument, by default `p=2`)
- if `init=2`, then surface deformation `η(t=0,x)=exp(-x^2)` and velocity `v(t=0,x)=exp.(-x^2)*(sin(x)+sin(K*x)/K^2)` (with `K` provided as an optional argument, by default `K=100`)

Other arguments are optional:
- `μ` the shallowness parameter (default is `1`),
- `ϵ` the nonlinearity parameter (default is `0.1`),
- `L` the half-length of the mesh (default is `20`),
- `N` the number of collocation points (default is `2^10`),
- `T` the final time of integration (default is `10`),
- `dt` the timestep (default is `0.001`),
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `1`, i.e. 2/3 rule, `0` means no dealiasing);
- `δ` the strength of the rectifier (default is `0.001`),
- `m` the order of the rectifier (as a regularizing operator, default is `-1`),
- `Ns` the number of stored computed times (default is all times).

Return `(problem,blowup_time,blowup,error_energy)` where
- `problem` contains all the information,
- `blowup_time` is the first time with NaN values (final time if there is none)
- `blowup` is a boolean indicating if NaN values occured
- `error_energy` is the relative energy preservation between first and final time
"""
function IntegrateWW2(;init=1,μ=1,ϵ=0.1,L=20,N=2^10,T=10,dt = 0.001,dealias=1,δ=0.001,m=-1,K=100,p=2,Ns=nothing)
	if isnothing(Ns)
		param = ( μ  = μ, ϵ  = ϵ,
				N  = N, L  = L,
	            T  = T, dt = dt )
	else
		param = ( μ  = μ, ϵ  = ϵ,
				N  = N, L  = L,
	            T  = T, dt = dt,
				Ns = Ns )
	end

	mesh=Mesh(param)
	if init == 1
		init = Init(x->exp.(-abs.(x).^p),x->zero(x))
	elseif init == 2
		init = Init(x->zero(x), x->exp.(-x.^2).*(sin.(x).+sin.(K*x)/K^2))
	else
		@error "argument init must be 1 or 2"
	end

	model = WWn(param;n = 2, δ = δ, m = m, dealias = dealias)
	problem = Problem(model, init, param)
	solve!( problem )

	(ηfin,vfin)   =  solution(problem)

	# check energy preservation
	x=mesh.x;k=mesh.k;
	Tmu = -1im*tanh.(sqrt(μ)*k)
	tmu = tanh.(sqrt(μ)*k)./k;tmu[1]=sqrt(μ);
	∂ₓ= 1im*k

	e(η,v) = 1/2* ( η.^2 .+ v.*ifft(tmu.*fft(v)) +ϵ*η.*(v.^2-ifft(Tmu.*fft(v)).^2) );
	e(η1,v1,η2,v2) = 1/2* ( (η1-η2).*(η1+η2) .+ (v1-v2).*ifft(tmu.*fft(v1+v2))
						+ϵ*(η1-η2).*(v1.^2-ifft(Tmu.*fft(v1)).^2)
						+ϵ*η2.*((v1-v2).*(v1+v2)-ifft(Tmu.*fft(v1-v2)).*ifft(Tmu.*fft(v1+v2))) );
	η0=init.η(x);v0=init.v(x)
	#error_energy_1 = abs(sum(e(ηfin,vfin)-e(η0,v0))/sum(e(η0,v0)))
	error_energy = abs(sum(e(ηfin,vfin,η0,v0))/sum(e(η0,v0)))
	@info "normalized preservation of the total energy: $error_energy\n"

	# compute blowup time (if any)
	blowup_index=0;blowup=false;
	for j in 1:problem.times.Ns
		blowup_index+=1
		if isnan(problem.data.U[j][1,1])
			blowup=true
			break
		end
	end
	blowup_time=	problem.times.ts[blowup_index]
	if blowup_time<T
		@info "The solution blowed at time t=$blowup_time (first NaN occurence)\n"
	end
	#display(plt)
	return problem,blowup_time,blowup,error_energy
end

#--- Figures

"""
	Figure(scenario;name,anim)

Several numerical experiments,
depending on the argument `scenario` (between 1 and 12)
corresponding to different figures in [DM]

- `scenario∈[1,7]`: spurious instability formation for smooth initial data.
    - `scenario=1`: no instability even without dealiasing and rectification, if small number of modes.
    - `scenario=2`: instabilities without dealiasing and rectification, if larger number of modes.
    - `scenario=3`: no instability with dealiasing and without rectification, if small number of modes.
    - `scenario=4`: instabilities with dealiasing and without rectification, if larger number of modes.
    - `scenario=5`: no instability in the presence of sufficiently regularizing rectifiers.
    - `scenario=6`: instabilities in the presence of insufficiently regularizing rectifiers.
    - `scenario=7`: mild instabilities in the presence of insufficiently strong rectifiers (δ small).
- `scenario∈[8,10]`: blowups for an initial data with a low-frequency and a high-frequency (with wavenumber K) component.
    - `scenario=8`: dependency of blowup times with respect to K (also `scenario=8.5` with a cut-off proportional to K)
    - `scenario=9`: dependency of blowup times with respect to ϵ.
    - `scenario=10`: an example of blowup with two different dealiasing.
- `scenario=11`: critical rectifier strength (δ) as a function of ϵ.
- `scenario≈12`: error of the rectified model (WW2) with respect to the water waves system, depending on δ and ϵ.
    - `scenario=12.1`: with initial data exp(-|x|)
    - `scenario=12.2`: with initial data exp(-|x|^2)
    - `scenario=12.3`: with initial data exp(-|x|^3)

Optional arguments are
- `compression`, an integer `m` used to reduce the size of figures by plotting only one every `m` points (when relevant),
- `name` (a string) used to save figures,
- `anim` generate animations (when relevant) if set to `true`.

Return relevant plots and problems depending on the situation.
"""
function Figure(scenario;compression=false,name=nothing,anim=false)

	if scenario == 1
		problem,=IntegrateWW2(init=1,μ=1,ϵ=0.1,L=20,N=2^9,T=10,dt = 0.001,dealias=0,δ=0,m=-1)
		plt=plot(problem;var=[:surface,:fourier],compression=compression,label="")
		if !isnothing(name)
			savefig(plt,"$name.pdf");savefig(plt,"$name.svg");
			if anim
				anim = @animate for t in LinRange(0,problem.times.tfin,101)
					plt = plot(problem;var=[:surface,:fourier],T=t,compression=compression)
					ylims!(plt[1],(-0.5,1.1))
				end
				gif(anim, "$name.gif", fps = 15)
			end
 		end
		return problem,plt

	elseif scenario == 2
		problem,=IntegrateWW2(init=1,μ=1,ϵ=0.1,L=20,N=2^11,T=1.5,dt = 0.001,dealias=0,δ=0,m=-1)
		plt=plot(problem,var=[:surface,:fourier],T=1.2;compression=compression,label="")
		if !isnothing(name)
			savefig(plt,"$name.pdf");savefig(plt,"$name.svg");
			if anim
				anim = @animate for t in LinRange(0,problem.times.tfin,101)
					plt = plot(problem;var=[:surface,:fourier],T=t,compression=compression)
					ylims!(plt[1],(-0.5,1.1))
				end
				gif(anim, "$name.gif", fps = 15)
			end
		end
		return problem,plt

	elseif scenario == 3
		problem,=IntegrateWW2(init=1,μ=1,ϵ=0.1,L=20,N=2^12,T=10,dt = 0.001,dealias=1,δ=0,m=-1)
		plt=plot(problem;var=[:surface,:fourier],compression=compression,label="")
		if !isnothing(name)
			savefig(plt,"$name.pdf");savefig(plt,"$name.svg");
			if anim
				anim = @animate for t in LinRange(0,problem.times.tfin,101)
					plt = plot(problem;var=[:surface,:fourier],T=t,compression=compression)
					ylims!(plt[1],(-0.5,1.1))
				end
				gif(anim, "$name.gif", fps = 15)
			end
		end
		return problem,plt

	elseif scenario == 4
		problem,=IntegrateWW2(init=1,μ=1,ϵ=0.1,L=20,N=2^14,T=1.5,dt = 0.001,dealias=1,δ=0,m=-1)
		plt=plot(problem,var=[:surface,:fourier],T=1.3;compression=compression,label="")
		if !isnothing(name)
			savefig(plt,"$name.pdf");savefig(plt,"$name.svg");
			if anim
				anim = @animate for t in LinRange(0,problem.times.tfin,101)
					plt = plot(problem;var=[:surface,:fourier],T=t,compression=compression)
					ylims!(plt[1],(-0.5,1.1))
				end
				gif(anim, "$name.gif", fps = 15)
			end
 		end
		return problem,plt

	elseif scenario == 5
		if anim Ns=nothing else Ns=1 end
		problem0,=IntegrateWW2(init=1,μ=1,ϵ=0.1,L=20,N=2^18,T=10,dt = 0.01,dealias=1,δ=0.01,m=-1/2,Ns=Ns)
		problem1,=IntegrateWW2(init=1,μ=1,ϵ=0.1,L=20,N=2^18,T=10,dt = 0.01,dealias=1,δ=0.01,m=-1,Ns=Ns)
		plt0=plot(problem0;var=[:surface,:fourier],compression=compression,label="")
		plt1=plot(problem1;var=[:surface,:fourier],compression=compression,label="")
		if !isnothing(name)
			savefig(plt0,"$name-r12.pdf");savefig(plt0,"$name-r12.svg");
			savefig(plt1,"$name-r1.pdf");savefig(plt1,"$name-r1.svg");
			if anim
				anim = @animate for t in LinRange(0,problem0.times.tfin,101)
					plt = plot([problem0,problem1];var=[:surface,:fourier],T=t,compression=compression,label=["m=-1/2" "m=-1"])
					ylims!(plt[1],(-0.5,1.1))
				end
				gif(anim, "$name.gif", fps = 15)
			end
		end
		return problem0,problem1,plt0,plt1

	elseif scenario == 6
		if anim Ns=nothing else Ns=10 end
		problem0,=IntegrateWW2(init=1,μ=1,ϵ=0.1,L=20,N=2^18,T=1,dt = 0.01,dealias=1,δ=0.01,m=-1/4,Ns=Ns)
		problem1,=IntegrateWW2(init=1,μ=1,ϵ=0.1,L=20,N=2^16,T=1,dt = 0.01,dealias=1,δ=0.01,m=-1/4,Ns=Ns)
		plt0=plot(problem0,var=[:surface,:fourier],T=0.6;compression=compression,label="")
		plt1=plot(problem1,var=[:surface,:fourier],T=0.6;compression=compression,label="")
		if !isnothing(name)
			savefig(plt0,"$name-N18.pdf");savefig(plt0,"$name-N18.svg");
			savefig(plt1,"$name-N16.pdf");savefig(plt1,"$name-N16.svg");
			if anim
				anim = @animate for t in LinRange(0,problem0.times.tfin,101)
					plt = plot([problem0,problem1];var=[:surface,:fourier],T=t,compression=compression,label=["N=2^18" "N=2^16"])
					ylims!(plt[1],(-0.5,1.1))
				end
				gif(anim, "$name.gif", fps = 15)
			end
		end
		return problem0,problem1,plt0,plt1

	elseif scenario in [7,7.1,7.2,7.3]
		if anim Ns=nothing else Ns=10 end
		p=round((scenario-7)*10);if p==0 p=2 end
		problem0,=IntegrateWW2(init=1,p=p,μ=1,ϵ=0.1,L=20,N=2^14,T=10,dt = 0.01,dealias=1,δ=0.01,m=-1,Ns=Ns)
		problem1,=IntegrateWW2(init=1,p=p,μ=1,ϵ=0.1,L=20,N=2^14,T=10,dt = 0.01,dealias=1,δ=0.002,m=-1,Ns=Ns)
		#problem2,=IntegrateWW2(init=1,p=p,μ=1,ϵ=0.1,L=20,N=2^14,T=2,dt = 0.01,dealias=1,δ=0.001,m=-1)
		plt0=plot(problem0;var=[:surface,:fourier],compression=compression,label="t=10")
		plot!(plt0,problem0,T=2;var=[:surface,:fourier],compression=compression,label="t=2")
		title!(plt0[1,1],"surface deformation")
		plt1=plot(problem1;var=[:surface,:fourier],compression=compression,label="t=10")
		plot!(plt1,problem1,T=2;var=[:surface,:fourier],compression=compression,label="t=2")
		title!(plt1[1,1],"surface deformation")
		if !isnothing(name)
			savefig(plt0,"$name-d01.pdf");savefig(plt0,"$name-d01.svg");
			savefig(plt1,"$name-d002.pdf");savefig(plt1,"$name-d002.svg");
			if anim
				anim = @animate for t in LinRange(0,problem0.times.tfin,101)
					plt = plot([problem0,problem1];var=[:surface,:fourier],T=t,compression=compression,label=["δ=0.01" "δ=0.002"])
					ylims!(plt[1],(-0.5,1.1))
				end
				gif(anim, "$name.gif", fps = 15)
			end
		end
		return problem0,problem1,plt0,plt1

	elseif scenario == 8
		blowups0=Real[];blowups1=Real[];
		K=100:20:800;iter=0;
		for k in K
			iter+=1;@info string("K=",k," (iteration ",iter,"/",length(K),")\n")
			~,blowup0=IntegrateWW2(init=2,K=k,μ=1,ϵ=0.15,L=20,N=2^14,T=10,dt = 0.001,dealias=1,δ=0,m=-1)
			~,blowup1=IntegrateWW2(init=2,K=k,μ=1,ϵ=0.2,L=20,N=2^14,T=10,dt = 0.001,dealias=1,δ=0,m=-1)
			push!(blowups0,blowup0);push!(blowups1,blowup1);
		end
		plt=scatter(K,[blowups0 blowups1],
				label=["ϵ=0.15" "ϵ=0.2"],
				marker=[:c :d],
				xlabel="K (in log scale)",
				ylabel="t (in log scale)",
				xscale=:log10,yscale=:log10)
		if !isnothing(name)
			savefig(plt,"$name.pdf");savefig(plt,"$name.svg");
			#@save(name,K,blowups0,blowups1)
		end
		return K,blowups0,blowups1,plt

	elseif scenario == 8.5
		K=100:20:800;
		blowups0=Real[];blowups1=Real[];iter=0;
		for k in K
			iter+=1;@info string("K=",k," (iteration ",iter,"/",length(K),")\n")
			N=2^14;d=N*π/k/20 # useful to define truncation at frequency K
			~,blowup0=IntegrateWW2(init=2,K=k,μ=1,ϵ=0.2,L=20,N=N,T=10,dt = 0.001,dealias=1,δ=0,m=-1)
			~,blowup1=IntegrateWW2(init=2,K=k,μ=1,ϵ=0.2,L=20,N=N,T=10,dt = 0.001,dealias=3/4*d-2,δ=0,m=-1)
			push!(blowups0,blowup0);push!(blowups1,blowup1);
		end
		plt=scatter(K,[blowups0 blowups1],
				label=["usual dealiasing" "adapted dealiasing"],
				marker=[:c :d],
				xlabel="K (in log scale)",
				ylabel="t (in log scale)",
				xscale=:log10,yscale=:log10)
		if !isnothing(name)
			savefig(plt,"$name.pdf");savefig(plt,"$name.svg");
			#@save(name,K,blowups0,blowups1)
		end
		return K,blowups0,blowups1,plt

	elseif scenario == 9
		Eps=10 .^(-1:0.02:0)
		blowups1=Real[];blowups2=Real[];blowups3=Real[];blowups4=Real[];iter=0;
		for eps in Eps
			iter+=1;@info string("ϵ=",eps," (iteration ",iter,"/",length(Eps),")\n")
			~,blowup1=IntegrateWW2(init=2,K=100,μ=1,ϵ=eps,L=20,N=2^14,T=10,dt = 0.001,dealias=1,δ=0,m=-1)
			~,blowup2=IntegrateWW2(init=2,K=200,μ=1,ϵ=eps,L=20,N=2^14,T=10,dt = 0.001,dealias=1,δ=0,m=-1)
			~,blowup3=IntegrateWW2(init=2,K=400,μ=1,ϵ=eps,L=20,N=2^14,T=10,dt = 0.001,dealias=1,δ=0,m=-1)
			~,blowup4=IntegrateWW2(init=2,K=800,μ=1,ϵ=eps,L=20,N=2^14,T=10,dt = 0.001,dealias=1,δ=0,m=-1)
			push!(blowups1,blowup1);push!(blowups2,blowup2);
			push!(blowups3,blowup3);push!(blowups4,blowup4);
		end
		plt=scatter(Eps,[blowups1 blowups2 blowups3 blowups4],
				label=["K=100" "K=200" "K=400" "K=800"],
				marker=[:c :d :ut :s],
				xlabel="ϵ (in log scale)",
				ylabel="t (in log scale)",
				xscale=:log10,yscale=:log10)
		if !isnothing(name)
			savefig(plt,"$name.pdf");savefig(plt,"$name.svg");
			#@save(name,Eps,blowups1,blowups2,blowups3,blowups4)
 		end
		return Eps,blowups1,blowups2,blowups3,blowups4,plt

	elseif scenario == 10
		K=400;N=2^14
		d=N*π/K/20 # useful to define truncation at frequency K
		problem0,=IntegrateWW2(init=2,K=K,μ=1,ϵ=0.2,L=20,N=N,T=2,dt = 0.001,dealias=1,δ=0,m=-1)
		problem1,=IntegrateWW2(init=2,K=K,μ=1,ϵ=0.2,L=20,N=N,T=2,dt = 0.001,dealias=3/4*d-2,δ=0,m=-1)

		plt0=plot(problem0, T=0.1;var=[:surface,:fourier],compression=compression,label="t=0.1")
		plot!(plt0,problem0,T=0.3;var=[:surface,:fourier],compression=compression,label="t=0.3")
		plot!(plt0,problem0,T=0.5;var=[:surface,:fourier],compression=compression,label="t=0.5")
		plot!(plt0[1,1],title="surface deformation")

		plt1=plot(problem1, T=0.1;var=[:surface,:fourier],compression=compression,label="t=0.1")
		plot!(plt1,problem1,T=0.3;var=[:surface,:fourier],compression=compression,label="t=0.3")
		plot!(plt1,problem1,T=0.5;var=[:surface,:fourier],compression=compression,label="t=0.5")
		plot!(plt1[1,1],title="surface deformation")
		if !isnothing(name)
			savefig(plt0,"$name-a.pdf");savefig(plt0,"$name-a.svg");
			savefig(plt1,"$name-b.pdf");savefig(plt1,"$name-b.svg");
			if anim
				anim = @animate for t in LinRange(0,problem0.times.tfin,101)
					plt = plot([problem0,problem1];var=[:surface,:fourier],T=t,compression=compression,label=["dealiasing" "low-pass filter"])
					ylims!(plt[1],(-0.5,1.1))
				end
				gif(anim, "$name.gif", fps = 15)

			end
		end
		return problem0,problem1,plt0,plt1

	elseif scenario == 11
		Eps=[0.01 0.0125 0.015 0.02 0.025 0.03 0.04 0.05 0.06 0.08 0.1 0.125 0.15 0.2 0.25 0.3 0.4 0.5 0.6 0.8]
		δc2=Real[];δc10=Real[];i=0;
		for eps in Eps
			i+=1
			#T=2
			δc_max=eps^2/2;δc_min=eps/500;
			for n=1:6
				@info string("iteration ",(i-1)*10+n,"/",10*length(Eps),"\n")
				@info string("critical ratio interval: ",(δc_min,δc_max))
				δc=sqrt(δc_max*δc_min)
				problem,blowuptime,blowup=IntegrateWW2(init=1,μ=1,ϵ=eps,L=20,N=2^20,T=2,dt = 0.005,dealias=1,δ=δc,m=-1,Ns=1)
				if blowup
					δc_min=δc
				else
					δc_max=δc
				end
			end
			push!(δc2,[δc_min,δc_max])
			#T=10
			δc_max=δc_min*1.5;
			for n=1:4
				@info string("iteration ",(i-1)*10+6+n,"/",10*length(Eps),"\n")
				@info string("critical ratio interval: ",(δc_min,δc_max))
				δc=sqrt(δc_max*δc_min)
				problem,blowuptime,blowup=IntegrateWW2(init=1,μ=1,ϵ=eps,L=20,N=2^20,T=10,dt = 0.005,dealias=1,δ=δc,m=-1,Ns=1)
				if blowup
					δc_min=δc
				else
					δc_max=δc
				end
			end
			push!(δc10,[δc_min,δc_max])
			#@save(name,Eps,δc2,δc10)
		end

		plt=plot()
		scatter!(plt,log10.(Eps[:]),[sqrt(δc10[i][2]*δc10[i][1]) for i in 1:length(Eps)],label="T=10",color=:1,marker=:c)
		scatter!(plt,log10.(Eps[:]),[sqrt(δc2[i][2]*δc2[i][1]) for i in 1:length(Eps)],label="T=2",color=:2,marker=:d)
		y2=OHLC[(δc2[i][1],δc2[i][1],δc2[i][2],δc2[i][2]) for i=1:length(Eps)]
		y10=OHLC[(δc10[i][1],δc10[i][1],δc10[i][2],δc10[i][2]) for i=1:length(Eps)]
		ohlc!(plt,log10.(Eps[:]),y10,label="",linecolor=:blue)
		ohlc!(plt,log10.(Eps[:]),y2,label="",linecolor=:orange)
		plot!(plt,log10.(Eps[:]),Eps[:].^2,label="δ=ϵ²",linecolor=:3)
		plot!(yscale=:log10,legend=:topleft,ylabel="δ (in log scale)",xlabel="ϵ (in log scale)")
		if !isnothing(name)
			savefig(plt,"$name.pdf");savefig(plt,"$name.svg");
			#@save(name,Eps,δc2,δc10)
		end
		return Eps,δc2,δc10,plt

	elseif scenario in [12,12.1,12.2,12.3]
		p=round((scenario-12)*10);if p==0 p=2 end

		Eps=[0.2,0.1,0.05]
		Delta=10 .^ (-2:0.1:0)
		Norms=Vector{Real}[];j=0;
		for eps in Eps
			norms=Real[];i=0;j+=1;
			for delta in Delta
				i+=1
				@info string("ϵ=",eps,", δ=",delta," (iteration ",(j-1)*length(Delta)+i,"/",length(Eps)*length(Delta),")")
				problem,=IntegrateWW2(init=1,p=p,μ=1,ϵ=eps,L=20,N=2^12,T=10,dt = 0.01,dealias=1,δ=delta,m=-1)
				param = (μ=1,ϵ=eps,L=20,N=2^12,T=10,dt = 0.01)
				model = WaterWaves(param; dealias = 1,method=1,maxiter=20)
				problem0 = Problem(model, problem.initial, param)
				solve!( problem0 )
				(η0,v0,x0)=solution(problem0)
				(η,v)=solution(problem,x=x0)
				push!(norms,norm(η-η0,Inf))
			end
			push!(Norms,norms)
		end

		plt=plot()
		markers=[:c :d :ut :s :r :lt :pent :hex :hep :oct :+ :x]
		for i in 1:length(Eps)
			scatter!(plt,Delta,Norms[i],label="ϵ=$Eps[i]",marker=markers[i])
		end
		plot!(plt,xlabel="δ",ylabel="error",legend=:topleft)
		if !isnothing(name)
			savefig(plt,"$name.pdf");savefig(plt,"$name.svg");
			plot!(plt,xscale=:log10,yscale=:log10)
			plot!(plt,xlabel="δ (in log scale)",ylabel="error (in log scale)")
			savefig(plt,"$name-log.pdf");savefig(plt,"$name-log.svg");
			#@save(name,Eps,Delta,Norms)
		end
		return Eps,Delta,Norms,plt

	else
		error("The first argument must be as specified in the documentation (basically between 1 and 12)")
	end

end
nothing


text:
"""
    isnull(c::Union{CXCursor,CLCursor}) -> Bool
Return true if cursor is null.
Wrapper for libclang's `clang_Cursor_isNull`.
"""
isnull(c::Union{CXCursor,CLCursor}) = clang_Cursor_isNull(c) != 0

"""
    getcursor() -> (NULL)CXCursor
Return the NULL CXCursor.
Wrapper for libclang's `clang_getNullCursor`.
"""
getcursor() = clang_getNullCursor()

# equality
Base.:(==)(c1::CXCursor, c2::CXCursor)::Bool = clang_equalCursors(c1, c2)
Base.hash(c::CXCursor) = clang_hashCursor(c)

"""
    kind(c::CXCursor) -> CXCursorKind
Return the kind of the given cursor.
Wrapper for libclang's `clang_getCursorKind`.
"""
kind(c::CXCursor) = clang_getCursorKind(c)

"""
    kind(c::CLCursor) -> CXCursorKind
Return the kind of the given cursor.
Note this method directly reads CXCursor's `kind` field, which won't invoke additional
`clang_getCursorKind` function calls.
"""
kind(c::CLCursor) = c.cursor.kind

"""
    isdecl(k::CXcursorKind) -> Bool
    isdecl(c::CLCursor) -> Bool
Return true if the given cursor kind represents a declaration.
Wrapper for libclang's `clang_isDeclaration`.
"""
isdecl(k::CXCursorKind)::Bool = clang_isDeclaration(k)
isdecl(c::CLCursor) = isdecl(kind(c))

"""
    isref(k::CXcursorKind) -> Bool
    isref(c::CLCursor) -> Bool
Return true if the given cursor kind represents a simple reference. Note that other kinds of
cursors (such as expressions) can also refer to other cursors. Use [`getref`](@ref) to determine
whether a particular cursor refers to another entity.
Wrapper for libclang's `clang_isReference`.
"""
isref(k::CXCursorKind)::Bool = clang_isReference(k)
isref(c::CLCursor) = isref(kind(c))

"""
    isexpr(k::CXcursorKind) -> Bool
    isexpr(c::CLCursor) -> Bool
Return true if the given cursor kind represents an expression.
Wrapper for libclang's `clang_isExpression`.
"""
isexpr(k::CXCursorKind)::Bool = clang_isExpression(k)
isexpr(c::CLCursor) = isexpr(kind(c))

"""
    isstmt(k::CXcursorKind) -> Bool
    isstmt(c::CLCursor) -> Bool
Return true if the given cursor kind represents a statement.
Wrapper for libclang's `clang_isStatement`.
"""
isstmt(k::CXCursorKind)::Bool = clang_isStatement(k)
isstmt(c::CLCursor) = isstat(kind(c))

"""
    isattr(k::CXcursorKind) -> Bool
    isattr(c::CLCursor) -> Bool
Return true if the given cursor kind represents an attribute.
Wrapper for libclang's `clang_isAttribute`.
"""
isattr(k::CXCursorKind)::Bool = clang_isAttribute(k)
isattr(c::CLCursor) = isattr(kind(c))

"""
    hasattr(c::Union{CXCursor,CLCursor}) -> Bool
Determine whether the given cursor has any attributes.
Wrapper for libclang's `clang_Cursor_hasAttrs`.
"""
hasattr(c::Union{CXCursor,CLCursor})::Bool = clang_Cursor_hasAttrs(c)

"""
    isvalid(k::CXcursorKind) -> Bool
    isvalid(c::CLCursor) -> Bool
Return true if the given cursor kind represents an valid cursor.
Wrapper for libclang's `clang_isInvalid`.
"""
Base.isvalid(k::CXCursorKind) = !(Bool(clang_isInvalid(k)))
Base.isvalid(c::CLCursor) = isvalid(kind(c))

"""
    is_translation_unit(k::CXcursorKind) -> Bool
    is_translation_unit(c::CLCursor) -> Bool
Return true if the given cursor kind represents a translation unit.
Wrapper for libclang's `clang_isTranslationUnit`.
"""
is_translation_unit(k::CXCursorKind)::Bool = clang_isTranslationUnit(k)
is_translation_unit(c::CLCursor) = is_translation_unit(kind(c))

"""
    ispreprocessing(k::CXcursorKind) -> Bool
    ispreprocessing(c::CLCursor) -> Bool
Return true if the given cursor kind represents a preprocessing element, such as a
preprocessor directive or macro instantiation.
Wrapper for libclang's `clang_isPreprocessing`.
"""
ispreprocessing(k::CXCursorKind)::Bool = clang_isPreprocessing(k)
ispreprocessing(c::CLCursor) = ispreprocessing(kind(c))

"""
    isunexposed(k::CXcursorKind) -> Bool
Return true if the given cursor kind represents a currently unexposed piece of the AST
(e.g., CXCursor_UnexposedStmt).
Wrapper for libclang's `clang_isUnexposed`.
"""
isunexposed(k::CXCursorKind)::Bool = clang_isUnexposed(k)
isunexposed(c::CLCursor) = isunexposed(kind(c))

"""
    linkage(c::Union{CXCursor,CLCursor}) -> CXLinkageKind
Return the linkage of the entity referred to by a given cursor.
Wrapper for libclang's `clang_getCursorLinkage`.
"""
linkage(c::Union{CXCursor,CLCursor})::CXLinkageKind = clang_getCursorLinkage(c)

## TODO:
# clang_getCursorVisibility
# clang_getCursorAvailability
# clang_getCursorPlatformAvailability
# clang_disposeCXPlatformAvailability
# clang_getCursorLanguage
# clang_getCursorTLSKind

"""
    get_translation_unit(c::Union{CXCursor,CLCursor}) -> CXTranslationUnit
Returns the translation unit that a cursor originated from.
Wrapper for libclang's `clang_Cursor_getTranslationUnit`.
"""
get_translation_unit(c::Union{CXCursor,CLCursor}) = clang_Cursor_getTranslationUnit(c)

## TODO: this is mainly for internal debugging?
# clang_createCXCursorSet
# clang_disposeCXCursorSet
# clang_CXCursorSet_contains
# clang_CXCursorSet_insert

"""
    get_semantic_parent(c::CXCursor) -> CXCursor
    get_semantic_parent(c::CLCursor) -> CLCursor
Return the semantic parent of the given cursor. Please checkout libclang's doc to know more.
Wrapper for libclang's `clang_getCursorSemanticParent`.
"""
get_semantic_parent(c::CXCursor) = clang_getCursorSemanticParent(c)
get_semantic_parent(c::CLCursor)::CLCursor = clang_getCursorSemanticParent(c)

"""
    get_lexical_parent(c::CXCursor) -> CXCursor
    get_lexical_parent(c::CLCursor) -> CLCursor
Return the lexical parent of the given cursor. Please checkout libclang's doc to know more.
Wrapper for libclang's `clang_getCursorLexicalParent`.
"""
get_lexical_parent(c::CXCursor) = clang_getCursorLexicalParent(c)
get_lexical_parent(c::CLCursor)::CLCursor = clang_getCursorLexicalParent(c)

## TODO:
# clang_getOverriddenCursors
# clang_disposeOverriddenCursors

"""
    get_included_file(c::Union{CXCursor,CLCursor}) -> CXFile
Return the file that is included by the given inclusion directive cursor.
Wrapper for libclang's `clang_getIncludedFile`.
"""
get_included_file(c::Union{CXCursor,CLCursor}) = clang_getIncludedFile(c)

"""
    location(c::Union{CXCursor,CLCursor}) -> CXSourceLocation
Return the physical location of the source constructor referenced by the given cursor.
Wrapper for libclang's `clang_getCursorLocation`.
"""
location(c::Union{CXCursor,CLCursor}) = clang_getCursorLocation(c)

"""
    extent(c::Union{CXCursor,CLCursor}) -> CXSourceRange
Return the physical extent of the source construct referenced by the given cursor.

The extent of a cursor starts with the file/line/column pointing at the first character
within the source construct that the cursor refers to and ends with the last character
within that source construct. For a declaration, the extent covers the declaration itself.
For a reference, the extent covers the location of the reference (e.g., where the referenced
entity was actually used).
Wrapper for libclang's `clang_getCursorExtent`.
"""
extent(c::Union{CXCursor,CLCursor}) = clang_getCursorExtent(c)

"""
    type(c::CXCursor) -> CXType
    type(c::CLCursor) -> CLType
Return the type of a CXCursor (if any). To get the cursor from a type, see [`typedecl`](@ref).
Wrapper for libclang's `clang_getCursorType`.
"""
type(c::CXCursor) = clang_getCursorType(c)
type(c::CLCursor)::CLType = clang_getCursorType(c)

"""
    underlying_type(c::CLTypedefDecl) -> CLType
Return the underlying type of a typedef declaration.
Wrapper for libclang's `clang_getTypedefDeclUnderlyingType`.
"""
underlying_type(c::CXCursor) = clang_getTypedefDeclUnderlyingType(c)
underlying_type(c::CLTypedefDecl)::CLType = clang_getTypedefDeclUnderlyingType(c)

"""
    integer_type(c::CLEnumDecl) -> CLType
Retrieve the integer type of an enum declaration.
Wrapper for libclang's `clang_getEnumDeclIntegerType`.
"""
integer_type(c::CXCursor) = clang_getEnumDeclIntegerType(c)
integer_type(c::CLEnumDecl)::CLType = clang_getEnumDeclIntegerType(c)

"""
    value(c::CLCursor) -> Int
Return the integer value of an enum constant declaration.
"""
function value(c::CLEnumConstantDecl)::Int
    typeKind = kind(type(c))
    if typeKind == CXType_Int || typeKind == CXType_Long || typeKind == CXType_LongLong
        return clang_getEnumConstantDeclValue(c)
    elseif typeKind == CXType_UInt || typeKind == CXType_ULong || typeKind == CXType_ULongLong
        return clang_getEnumConstantDeclUnsignedValue(c)
    end
    error("Unknown EnumConstantDecl type: ", typeKind, " cursor: ", kind(c))
end

"""
    bitwidth(c::CLFieldDecl) -> Int
Return the bit width of a bit field declaration as an integer.
Wrapper for libclang's `clang_getFieldDeclBitWidth`.
"""
bitwidth(c::CXCursor)::Int = clang_getFieldDeclBitWidth(c)
bitwidth(c::CLFieldDecl) = clang_getFieldDeclBitWidth(c)

"""
    argnum(c::Union{CXCursor,CLFunctionDecl,CLCXXMethod}) -> Int
Return the number of non-variadic arguments associated with a given cursor.
Wrapper for libclang's `clang_Cursor_getNumArguments`.
"""
argnum(c::Union{CXCursor,CLFunctionDecl,CLCXXMethod})::Int = clang_Cursor_getNumArguments(c)

"""
    argument(c::CXCursor, i::Integer) -> CXCursor
    argument(c::Union{CLFunctionDecl,CLCXXMethod}, i::Integer) -> CLCursor
Return the argument cursor of a function or method.
Wrapper for libclang's `clang_Cursor_getArgument`.
"""
argument(c::CXCursor, i::Integer) = clang_Cursor_getArgument(c, Unsigned(i))
argument(c::Union{CLFunctionDecl,CLCXXMethod}, i::Integer)::CLCursor = clang_Cursor_getArgument(c, Unsigned(i))

## TODO:
# clang_Cursor_getNumTemplateArguments
# clang_Cursor_getTemplateArgumentKind
# clang_Cursor_getTemplateArgumentType
# clang_Cursor_getTemplateArgumentValue
# clang_Cursor_getTemplateArgumentUnsignedValue

"""
    isfunctionlike(c::Union{CXCursor,CLCursor}) -> Bool
Determine whether a CXCursor that is a macro, is function like.
Wrapper for libclang's `clang_Cursor_isMacroFunctionLike`.
"""
isfunctionlike(c::Union{CXCursor,CLCursor})::Bool = clang_Cursor_isMacroFunctionLike(c)

"""
    isbuiltin(c::Union{CXCursor,CLCursor}) -> Bool
Determine whether a  CXCursor that is a macro, is a builtin one.
Wrapper for libclang's `clang_Cursor_isMacroBuiltin`.
"""
isbuiltin(c::Union{CXCursor,CLCursor})::Bool = clang_Cursor_isMacroBuiltin(c)

"""
    isinlined(c::Union{CXCursor,CLCursor}) -> Bool
Determine whether a CXCursor that is a function declaration, is an inline declaration.
Wrapper for libclang's `clang_Cursor_isFunctionInlined`.
"""
isinlined(c::Union{CXCursor,CLCursor})::Bool = clang_Cursor_isFunctionInlined(c)

"""
    result_type(c::CXCursor) -> CXType
    result_type(c::Union{CLFunctionDecl,CLCXXMethod}) -> CLType
Return the return type associated with a given cursor. This only returns a valid type if
the cursor refers to a function or method.
Wrapper for libclang's `clang_getCursorResultType`.
"""
result_type(c::CXCursor) = clang_getCursorResultType(c)
result_type(c::Union{CLFunctionDecl,CLCXXMethod})::CLType = clang_getCursorResultType(c)

## TODO:
# clang_getCursorExceptionSpecificationType
# clang_Cursor_getOffsetOfField
#

"""
    isanonymous(c::Union{CXCursor,CLCursor}) -> Bool
Return true if the given cursor represents an anonymous record declaration(C++).
Wrapper for libclang's `clang_Cursor_isAnonymous`.
"""
isanonymous(c::Union{CXCursor,CLCursor})::Bool = clang_Cursor_isAnonymous(c)

"""
    isbit(c::Union{CXCursor,CLCursor}) -> Bool
Return true if the cursor specifies a Record member that is a bitfield.
Wrapper for libclang's `clang_Cursor_isBitField`.
"""
isbit(c::Union{CXCursor,CLCursor}) = clang_Cursor_isBitField(c) != 0

"""
    isvirtual(c::Union{CXCursor,CLCursor}) -> Bool
Return true if the base class specified by the cursor with kind CX_CXXBaseSpecifier is virtual.
Wrapper for libclang's `clang_isVirtualBase`.
"""
isvirtual(c::Union{CXCursor,CLCursor})::Bool = clang_isVirtualBase(c)

## TODO:
# clang_getCXXAccessSpecifier
# clang_Cursor_getStorageClass
# clang_getNumOverloadedDecls
# clang_getOverloadedDecl
# clang_getIBOutletCollectionType
#

"""
    spelling(c::Union{CXCursor,CLCursor}) -> String
Return a name for the entity referenced by this cursor.
"""
function spelling(c::Union{CXCursor,CLCursor})
    cxstr = clang_getCursorSpelling(c)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

# clang_Cursor_getSpellingNameRange

"""
    name(c::Union{CXCursor,CLCursor}) -> String
Return the display name for the entity referenced by this cursor.
"""
function name(c::Union{CXCursor,CLCursor})
    cxstr = clang_getCursorDisplayName(c)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

"""
    getref(c::CXCursor) -> CXCursor
    getref(c::CLCursor) -> CLCursor
For a cursor that is a reference, retrieve a cursor representing the entity that it references.
Wrapper for libclang's `clang_getCursorReferenced`.
"""
getref(c::CXCursor) = clang_getCursorReferenced(c)
getref(c::CLCursor)::CLCursor = clang_getCursorReferenced(c)

"""
    getdef(c::CXCursor) -> CXCursor
    getdef(c::CLCursor) -> CLCursor
For a cursor that is either a reference to or a declaration of some entity, retrieve a cursor
that describes the definition of that entity.
Wrapper for libclang's `clang_getCursorDefinition`.
"""
getdef(c::CXCursor) = clang_getCursorDefinition(c)
getdef(c::CLCursor)::CLCursor = clang_getCursorDefinition(c)

"""
    isdef(c::Union{CXCursor,CLCursor}) -> Bool
Return true if the declaration pointed to by this cursor is also a definition of that entity.
Wrapper for libclang's `clang_isCursorDefinition`.
"""
isdef(c::Union{CXCursor,CLCursor})::Bool = clang_isCursorDefinition(c)

"""
    canonical(c::CXCursor) -> CXCursor
    canonical(c::CLCursor) -> CLCursor
Return the canonical cursor corresponding to the given cursor.
Wrapper for libclang's `clang_getCanonicalCursor`.
"""
canonical(c::CXCursor) = clang_getCanonicalCursor(c)
canonical(c::CLCursor)::CLCursor = clang_getCanonicalCursor(c)

## TODO:
# clang_Cursor_getObjCSelectorIndex
# clang_Cursor_isDynamicCall
# clang_Cursor_getReceiverType

"""
    isvariadic(c::Union{CXCursor,CLCursor}) -> Bool
Return true if the given cursor is a variadic function or method.
Wrapper for libclang's `clang_Cursor_isVariadic`.
"""
isvariadic(c::Union{CXCursor,CLCursor}) = clang_Cursor_isVariadic(c) != 0

# clang_Cursor_isExternalSymbol
# clang_Cursor_getCommentRange
# clang_Cursor_getRawCommentText
# clang_Cursor_getBriefCommentText

function spelling(k::CXCursorKind)
    cxstr = clang_getCursorKindSpelling(k)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

# helper
"""
    filename(c::Union{CXCursor,CLCursor}) -> String
Return the complete file and path name of the given file referenced by the input cursor.
"""
function filename(c::Union{CXCursor,CLCursor})
    file = Ref{CXFile}(C_NULL)
    line = Ref{Cuint}(0)
    column = Ref{Cuint}(0)
    offset = Ref{Cuint}(0)
    location = clang_getCursorLocation(c)
    clang_getExpansionLocation(location, file, line, column, offset)
    if file[] != C_NULL
        cxstr = clang_getFileName(file[])
        ptr = clang_getCString(cxstr)
        s = unsafe_string(ptr)
        clang_disposeString(cxstr)
        return s
    else
        return ""
    end
end

"""
    return_type(c::CLCursor, resolve::Bool=true) -> CXtype
Return the return type associated with a function/method cursor.
"""
function return_type(c::CLCursor, resolve::Bool=true)
    t = result_type(c)
    resolve && return resolve_type(t)
    return t
end

"""
    typedef_type(c::CLCursor) -> CXType
Return the underlying type of a typedef declaration.
"""
function typedef_type(c::CLCursor)
    t = underlying_type(c)
    kind(t) == CXType_Unexposed && return type(children(c)[1])
    return t
end

"""
    search(cursors::Vector{CLCursor}, ismatch::Function) -> Vector{CLCursor}
Return vector of CLCursors that match predicate. `ismatch` is a function that accepts a CLCursor argument.
"""
function search(cursors::Vector{T}, ismatch::Function) where {T<:CLCursor}
    matched = CLCursor[]
    for cursor in cursors
        ismatch(cursor) && push!(matched, cursor)
    end
    return matched
end
search(c::CLCursor, ismatch::Function) = search(children(c), ismatch)
search(c::CLCursor, s::String) = search(c, x->spelling(x)==s)
search(c::CLCursor, k::CXCursorKind) = search(c, x->kind(x)==k)

# visitor
function cu_children_visitor(cursor::CXCursor, parent::CXCursor, list)::Cuint
    push!(list, cursor)
    return CXChildVisit_Continue
end

"""
    children(cursor::CXCursor) -> Vector{CXCursor}
    children(cursor::CLCursor) -> Vector{CLCursor}
Return a child cursor vector of the given cursor.
"""
function children(cursor::CXCursor)
    # TODO: possible to use sizehint! here?
    list = CXCursor[]
    cu_visitor_cb = @cfunction(cu_children_visitor, Cuint, (CXCursor, CXCursor, Ref{Vector{CXCursor}}))
    GC.@preserve list ccall((:clang_visitChildren, LibClang.libclang), UInt32, (CXCursor, CXCursorVisitor, Any), cursor, cu_visitor_cb, list)
    return list
end
children(c::CLCursor)::Vector{CLCursor} = children(c.cursor)

"""
    function_args(cursor::CLCursor) -> Vector{CLCursor}
Return function arguments for a given cursor.
"""
function_args(cursor::CLCursor) = search(cursor, CXCursor_ParmDecl)

"""
    is_typedef_anon(current::CLCursor, next::CLCursor) -> Bool
Return true if the current cursor is an typedef anonymous struct/enum.
"""
function is_typedef_anon(current::CLCursor, next::CLCursor)
    !isempty(name(current)) && return false
    refback = children(next)
    if kind(next) == CXCursor_TypedefDecl &&
       !isempty(refback) && isempty(name(refback[1]))
        return true
    else
        return false
    end
end

"""
    is_forward_declaration(cursor::CLCursor) -> Bool
Return true if the current cursor is a forward declaration.
"""
function is_forward_declaration(cursor::CLCursor)
    definition = clang_getCursorDefinition(cursor)

    # If the definition is null, then there is no definition in this translation
    # unit, so this cursor must be a forward declaration.
    if clang_equalCursors(definition, clang_getNullCursor()) == 1
        return true
    end

    # If there is a definition, then the forward declaration and the definition
    # are in the same translation unit. This cursor is the forward declaration if
    # it is not the definition.
    return clang_equalCursors(cursor, definition) == 0
end


text:
module Gate
using ..Schrodinger
using SparseArrays

"""
    Gate.rotation(θ, n=(1,0,0))

Generate a qubit rotation operator giving a `θ` rad rotation about an axis defined by the vector ``\\vec{n}``. Note that `n` will be normalized, allowing for inputs like `(1,1,0)`. The rotation operator is defined as

```math
\\hat{R}_{\\vec{n}}(θ) = \\exp\\left(−iθ\\vec{n}⋅\\vec{σ}/2\\right) = \\cos\\frac{θ}{2}𝕀 - i\\sin\\frac{θ}{2}(n_x σ_x + n_y σ_y + n_z σ_z).
```
"""
function rotation(θ::Real,n::NTuple{3,Real}=(1,0,0))
    a = 1/√(n[1]^2+n[2]^2+n[3]^2)
    nx,ny,nz = a*n[1],a*n[2],a*n[3]
    c = cos(θ/2); s = sin(θ/2)
    R = Matrix{Complex{typeof(nz*s)}}(undef,2,2)
    R[1,1] = complex(c,-nz*s)
    R[2,1] = complex(ny*s,-nx*s)
    R[1,2] = complex(-ny*s,-nx*s)
    R[2,2] = complex(c,nz*s)
    return Operator(R,(2,))
end

"""
    Gate.H

The Hadamard gate:
```math
H = \\frac{1}{\\sqrt 2}
\\begin{pmatrix}
1 & 1 \\\\
1 & -1
\\end{pmatrix}
```

The Hadamard is a ``π`` rotation about the axis ``\\vec{n} = (1,0,1)``, plus a global ``i`` phase.
"""
const H = Operator([√0.5 √0.5; √0.5 -√0.5],(2,))

"""
    Gate.S

The phase gate:
```math
S =
\\begin{pmatrix}
1 & 0 \\\\
0 & i
\\end{pmatrix}
```

The phase gate is the square of the T gate: ``S = T^2``.
"""
const S = Operator(sparse([1 0; 0 1im]),(2,))

"""
    Gate.T

The T gate:
```math
T =
\\begin{pmatrix}
1 & 0 \\\\
0 & \\exp(iπ/4)
\\end{pmatrix}
```

The phase gate is the square of the T gate: ``S = T^2``.
"""
const T = Operator(sparse([1 0; 0 complex(√0.5,√0.5)]),(2,))

const cNOT = Operator(sparse([1 0 0 0;0 1 0 0;0 0 0 1;0 0 1 0]),(2,2))
const rcNOT = Operator(sparse([0 1 0 0;1 0 0 0;0 0 1 0;0 0 0 1]),(2,2))
const cZ = Operator(sparse([1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 -1]),(2,2))

end


text:
using ChainRulesCore
using LinearAlgebra
using KrylovKit

@Zygote.nograd StopFunction
@Zygote.nograd leftorth
@Zygote.nograd rightorth
@Zygote.nograd _initializect_square
@Zygote.nograd printstyled
@Zygote.nograd save
@Zygote.nograd load
@Zygote.nograd error

# patch since it's currently broken otherwise
function ChainRulesCore.rrule(::typeof(Base.typed_hvcat), ::Type{T}, rows::Tuple{Vararg{Int}}, xs::S...) where {T,S}
    y = Base.typed_hvcat(T, rows, xs...)
    function back(ȳ)
        return NO_FIELDS, NO_FIELDS, NO_FIELDS, permutedims(ȳ)...
    end
    return y, back
end

# improves performance compared to default implementation, also avoids errors
# with some complex arrays
function ChainRulesCore.rrule(::typeof(LinearAlgebra.norm), A::AbstractArray)
    n = norm(A)
    function back(Δ)
        return NO_FIELDS, Δ .* A ./ (n + eps(0f0)), NO_FIELDS
    end
    return n, back
end

"""
    ChainRulesCore.rrule(::typeof(leftenv), AL::AbstractArray{T}, M::AbstractArray{T}, FL::AbstractArray{T}; kwargs...) where {T}

```
           ┌──  AL ──┐
           │    │    │
dM   =  - FL ──   ── ξl
           │    │    │
           └──  AL ──┘

           ┌──     ──┐         ┌──  AL ──┐
           │    │    │         │    │    │
dAL  =  - FL ── M ── ξl   -   FL ── M ── ξl
           │    │    │         │    │    │
           └──  AL ──┘         └──     ──┘
```
"""
function ChainRulesCore.rrule(::typeof(leftenv), AL::AbstractArray{T}, M::AbstractArray{T}, FL::AbstractArray{T}; kwargs...) where {T}
    λl, FL = leftenv(AL, M, FL)
    # @show λl
    function back((dλ, dFL))
        ξl, info = linsolve(FR -> ein"ηpβ,βaα,csap,γsα -> ηcγ"(AL, FR, M, conj(AL)), permutedims(dFL, (3, 2, 1)), -λl, 1)
        # errL = ein"abc,cba ->"(FL, ξl)[]
        # abs(errL) > 1e-1 && throw("FL and ξl aren't orthometric. err = $(errL)")
        dAL = -ein"γcη,csap,γsα,βaα -> ηpβ"(FL, M, conj(AL), ξl) - ein"γcη,csap,ηpβ,βaα -> γsα"(FL, M, AL, ξl)
        dM = -ein"γcη,ηpβ,γsα,βaα -> csap"(FL, AL, conj(AL), ξl)
        return NO_FIELDS, dAL, dM, NO_FIELDS...
    end
    return (λl, FL), back
end

"""
    ChainRulesCore.rrule(::typeof(rightenv), AR::AbstractArray{T}, M::AbstractArray{T}, FR::AbstractArray{T}; kwargs...) where {T}

```
          ┌──  AR ──┐
          │    │    │
dM   =   ξr ──   ── FR
          │    │    │
          └──  AR ──┘

          ┌──     ──┐         ┌──  AR ──┐
          │    │    │         │    │    │
dAR  =   ξr ── M ── FR   +   ξr ── M ── FR
          │    │    │         │    │    │
          └──  AR ──┘         └──     ──┘
```
"""
function ChainRulesCore.rrule(::typeof(rightenv), AR::AbstractArray{T}, M::AbstractArray{T}, FR::AbstractArray{T}; kwargs...) where {T}
    λr, FR = rightenv(AR, M, FR)
    # @show λr
    function back((dλ, dFR))
        ξr, info = linsolve(FL -> ein"ηpβ,γcη,csap,γsα -> αaβ"(AR, FL, M, conj(AR)), permutedims(dFR, (3, 2, 1)), -λr, 1)
        # errR = ein"abc,cba ->"(ξr, FR)[]
        # abs(errR) > 1e-1 && throw("FR and ξr aren't orthometric. err = $(errR)")
        dAR = -ein"γcη,csap,γsα,βaα -> ηpβ"(ξr, M, conj(AR), FR) - ein"γcη,csap,ηpβ,βaα -> γsα"(ξr, M, AR, FR)
        dM = -ein"γcη,ηpβ,γsα,βaα -> csap"(ξr, AR, conj(AR), FR)
        return NO_FIELDS, dAR, dM, NO_FIELDS...
    end
    return (λr, FR), back
end

"""
    ChainRulesCore.rrule(::typeof(ACenv),AC::AbstractArray{T}, FL::AbstractArray{T}, M::AbstractArray{T}, FR::AbstractArray{T};

```
          ┌──  AC ──┐
          │    │    │
dM   =   FL ──   ── FR
          │    │    │
          └──  ξ  ──┘

          ┌──  AC ──┐
          │    │    │
dFL  =      ── M ── FR
          │    │    │
          └──  ξ  ──┘

          ┌──  AC ──┐
          │    │    │
dFR  =   FL ── M ──
          │    │    │
          └──  ξ  ──┘
```
"""
function ChainRulesCore.rrule(::typeof(ACenv),AC::AbstractArray{T}, FL::AbstractArray{T}, M::AbstractArray{T}, FR::AbstractArray{T};
    kwargs...) where {T}
    λAC, AC = ACenv(AC, FL, M, FR)
    # @show λAC
    function back((dλ, dAC))
        ξ, info = linsolve(AC -> ein"αaγ,αsβ,asbp,ηbβ -> γpη"(FL, AC, M, FR), dAC, -λAC, 1)
        # errAC = ein"abc,abc ->"(AC, ξ)[]
        # abs(errAC) > 1e-1 && throw("AC and ξ aren't orthometric. err = $(errAC)")
        dFL = -ein"ηpβ,βaα,csap,γsα -> γcη"(AC, FR, M, ξ)
        dM = -ein"γcη,ηpβ,γsα,βaα -> csap"(FL, AC, ξ, FR)
        dFR = -ein"ηpβ,γcη,csap,γsα -> βaα"(AC, FL, M, ξ)
        return NO_FIELDS, NO_FIELDS, dFL, dM, dFR
    end
    return (λAC, AC), back
end

"""
    ChainRulesCore.rrule(::typeof(ACenv),AC::AbstractArray{T}, FL::AbstractArray{T}, M::AbstractArray{T}, FR::AbstractArray{T};

```
          ┌──  C  ──┐
          │         │
dFL  =      ─────── FR
          │         │
          └──  ξ  ──┘

          ┌──  C  ──┐
          │         │
dFR  =   FL ───────
          │         │
          └──  ξ  ──┘
```
"""
function ChainRulesCore.rrule(::typeof(Cenv), C::AbstractArray{T}, FL::AbstractArray{T}, FR::AbstractArray{T}; kwargs...) where {T}
    λC, C = Cenv(C, FL, FR)
    # @show λC
    function back((dλ, dC))
        ξ, info = linsolve(C -> ein"αaγ,αβ,ηaβ -> γη"(FL, C, FR), dC, -λC, 1)
        # errC = ein"ab,ab ->"(C, ξ)[]
        # abs(errC) > 1e-1 && throw("C and ξ aren't orthometric. err = $(errC)")
        # @show info ein"ab,ab ->"(C,ξ)[] ein"γp,γp -> "(C,dC)[]
        dFL = -ein"ηβ,βaα,γα -> γaη"(C, FR, ξ)
        dFR = -ein"ηβ,γcη,γα -> βcα"(C, FL, ξ)
        return NO_FIELDS, NO_FIELDS, dFL, dFR
    end
    return (λC, C), back
end

# adjoint for QR factorization
# https://journals.aps.org/prx/abstract/10.1103/PhysRevX.9.031041 eq.(5)
function ChainRulesCore.rrule(::typeof(qrpos), A::AbstractArray{T,2}) where {T}
    Q, R = qrpos(A)
    function back((dQ, dR))
        M = R * dR' - dQ' * Q
        dA = (UpperTriangular(R + I * 1e-6) \ (dQ + Q * Symmetric(M, :L))' )'
        return NO_FIELDS, dA
    end
    return (Q, R), back
end

function ChainRulesCore.rrule(::typeof(lqpos), A::AbstractArray{T,2}) where {T}
    L, Q = lqpos(A)
    function back((dL, dQ))
        M = L' * dL - dQ * Q'
        dA = LowerTriangular(L + I * 1e-6)' \ (dQ + Symmetric(M, :L) * Q)
        return NO_FIELDS, dA
    end
    return (L, Q), back
end

"""
    ChainRulesCore.rrule(::typeof(bigleftenv), AL::AbstractArray{T}, M::AbstractArray{T}, FL::AbstractArray{T}; kwargs...) where {T}

```
       ┌── AL──┐     ┌── AL──┐
       │   │   │     │   │   │
       │──   ──│     │── M ──│
dM = - FL  │   ξl  - FL  │   ξl
       │── M ──│     │──   ──│
       │   │   │     │   │   │
       ┕── AL──┘     ┕── AL──┘

        ┌──   ──┐     ┌── AL──┐
        │   │   │     │   │   │
        │── M ──│     │── M ──│
dAL = - FL  │   ξl  - FL  │   ξl
        │── M ──│     │── M ──│
        │   │   │     │   │   │
        ┕── AL──┘     ┕──   ──┘

```
"""
function ChainRulesCore.rrule(::typeof(bigleftenv), AL::AbstractArray{T}, M::AbstractArray{T}, FL4::AbstractArray{T}; kwargs...) where {T}
    λl, FL4 = bigleftenv(AL, M, FL4; kwargs...)
    # @show λl
    function back((dλl, dFL4))
        ξl, info = linsolve(FR4 -> ein"fghi,def,ckge,bjhk,aji -> dcba"(FR4,AL,M,M,conj(AL)), dFL4, -λl, 1)
        # errL = ein"abc,cba ->"(FL4, ξl)[]
        # abs(errL) > 1e-1 && throw("FL and ξl aren't orthometric. err = $(errL)")
        dAL = -ein"dcba,ckge,bjhk,aji,fghi -> def"(FL4, M, M, AL, ξl) -ein"dcba,def,ckge,bjhk,fghi -> aji"(FL4, AL, M, M, ξl)
        dM = -ein"dcba,def,bjhk,aji,fghi -> ckge"(FL4, AL, M, AL, ξl) -ein"dcba,def,ckge,aji,fghi -> bjhk"(FL4, AL, M, AL, ξl)
        return NO_FIELDS, dAL, dM, NO_FIELDS...
    end
    return (λl, FL4), back
end

"""
    ChainRulesCore.rrule(::typeof(rightenv), AR::AbstractArray{T}, M::AbstractArray{T}, FR::AbstractArray{T}; kwargs...) where {T}

```
       ┌── AR──┐     ┌── AR──┐
       │   │   │     │   │   │
       │──   ──│     │── M ──│
dM = - ξr  │   FR  - ξr  │   FR
       │── M ──│     │──   ──│
       │   │   │     │   │   │
       ┕── AR──┘     ┕── AR──┘

        ┌──   ──┐     ┌── AR──┐
        │   │   │     │   │   │
        │── M ──│     │── M ──│
dAR = - ξr  │   FR  - ξr  │   FR
        │── M ──│     │── M ──│
        │   │   │     │   │   │
        ┕── AR──┘     ┕──   ──┘
```
"""
function ChainRulesCore.rrule(::typeof(bigrightenv), AR::AbstractArray{T}, M::AbstractArray{T}, FR4::AbstractArray{T}; kwargs...) where {T}
    λr, FR4 = bigrightenv(AR, M, FR4; kwargs...)
    # @show λr
    function back((dλ, dFR4))
        ξr, info = linsolve(FL4 -> ein"dcba,def,ckge,bjhk,aji -> fghi"(FL4,AR,M,M,conj(AR)), dFR4, -λr, 1)
        # errR = ein"abc,cba ->"(ξr, FR4)[]
        # abs(errR) > 1e-1 && throw("FR and ξr aren't orthometric. err = $(errR)")
        dAR = -ein"dcba,ckge,bjhk,aji,fghi -> def"(ξr, M, M, AR,FR4) -ein"dcba,def,ckge,bjhk,fghi -> aji"(ξr, AR, M, M, FR4)
        dM = -ein"dcba,def,bjhk,aji,fghi -> ckge"(ξr, AR, M, AR, FR4) -ein"dcba,def,ckge,aji,fghi -> bjhk"(ξr, AR, M, AR, FR4)
        return NO_FIELDS, dAR, dM, NO_FIELDS...
    end
    return (λr, FR4), back
end

@doc raw"
    num_grad(f, K::Real; [δ = 1e-5])

return the numerical gradient of `f` at `K` calculated with
`(f(K+δ/2) - f(K-δ/2))/δ`

# example

```jldoctest; setup = :(using TensorNetworkAD)
julia> TensorNetworkAD.num_grad(x -> x * x, 3) ≈ 6
true
```
"
num_grad(f, K::Real; δ::Real=1e-5) = (f(K + δ / 2) - f(K - δ / 2)) / δ

@doc raw"
    num_grad(f, K::AbstractArray; [δ = 1e-5])

return the numerical gradient of `f` for each element of `K`.

# example

```jldoctest; setup = :(using TensorNetworkAD, LinearAlgebra)
julia> TensorNetworkAD.num_grad(tr, rand(2,2)) ≈ I
true
```
"
function num_grad(f, a::AbstractArray; δ::Real=1e-5)
    map(CartesianIndices(a)) do i
        foo = x -> (ac = copy(a); ac[i] = x; f(ac))
        num_grad(foo, a[i], δ=δ)
    end
end


text:
using CSV
using DataFrames

function parse_committees(chamber)
    DF = DataFrame!(CSV.File("$(chamber)_committees.csv"))

    sub = combine(groupby(DF, :new_id)) do sdf
               sdf[argmax(sdf.new_id), :]
               end

    labels = sub[:, "party"]

    hedges = Vector{Vector{Int64}}()
    for sub ∈ groupby(DF, [:session, :committee])
        hedge = sub[!, :new_id]
        sort!(hedge)
        push!(hedges, hedge)
    end

    names = ["$f $l" for (f, l) in zip(sub[!, "first"], sub[!, "last"])]

    dir = string("$(chamber)-committees")
    if !isdir(dir)
        mkdir(dir)
    end

    # hyperedge data
    open("$dir/hyperedges-$dir.txt", "w") do f
        for hedge in hedges
            write(f, join(hedge, ','))
            write(f, "\n")
        end
    end

    # names of the integer labels
    open("$dir/label-names-$dir.txt", "w") do f
        write(f, "Democrat\n")
        write(f, "Republican\n")
    end

    # labels for each node
    open("$dir/node-labels-$dir.txt", "w") do f
        for party in labels
            write(f, "$(party)\n")
        end
    end

    # names of the nodes
    open("$dir/node-names-$dir.txt", "w") do f
        for name in names
            write(f, "$(name)\n")
        end
    end
end

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: nveldt/TestPkg

max_forks_repo_name: PhilChodrow/HypergraphModularity