text:
function gdfPlotLabels(tn::AbstractTimenum, args...;
                       shiftLabels::Float64 = 0.)
    ## create usual line plot data
    nams = names(tn)
    tn = narowrm(tn)
    vals = asArr(tn, Float64, NaN) # NA to NaN
    dats = dat2num(tn)
    df = composeDataFrame([dats vals], [:Idx; nams])
    stackedData = stack(df, nams)

    ## create label information
    nLines = length(nams)
    xVals = Float64[dats[end] + shiftLabels for ii=1:nLines]
    yVals = Float64[get(tn[end, ii], 1, 1) for ii=1:nLines]
    lineLabels = AssetMgmt.symbToStr(nams)

    ## put label info in single DataFrame
    dfLabels = DataFrame(x = xVals, y = yVals,
                         labels = lineLabels)

    Gadfly.plot(layer(stackedData, x="Idx",
                      y="value",
                      color="variable",
                      Gadfly.Geom.line
                      ),
                layer(dfLabels, x="x", y="y",
                      label="labels",
                      Gadfly.Geom.label(position=:dynamic,
                                        hide_overlaps=false)
                      ), args...)
end

##################################
## visualize colored line paths ##
##################################

function gdfGroupPlot(tn::AbstractTimenum, args...;
                      variableInfo::DataFrame = DataFrame(),
                      joinCol::Symbol = :AssetLabel,
                      colorCol::Symbol = :AssetClass,
                      labelCol::Symbol = :AssetLabel,
                      shiftLabels::Float64 = 0.)

    ## convert strings in key column to symbols for join
    if !(eltype(variableInfo[joinCol]) == Symbol)
        ## conversion required for join
        variableInfo = deepcopy(variableInfo)
        variableInfo[joinCol] = AssetMgmt.strToSymb(variableInfo[joinCol])
    end

    ## get variable names
    nams = names(tn)

    ## standard data manipulation:
    ## remove / replace NAs, convert dates
    tn = narowrm(tn)
    vals = asArr(tn, Float64, NaN) # NA to NaN
    dats = dat2num(tn)
    df = composeDataFrame([dats vals], [:Idx; nams])
    stackedData = stack(df, nams)

    ## create label information
    nLines = length(nams)
    xVals = Float64[dats[end] + shiftLabels for ii=1:nLines]
    yVals = Float64[get(tn[end, ii], 1, 1) for ii=1:nLines]
    colNames = nams
    dfLabelCoordinates = DataFrame(variable = colNames,
                                   x = xVals, y = yVals)

    # rename column name
    rename!(dfLabelCoordinates, :variable, joinCol)

    ## join labels from assetInfo table
    subInfo = variableInfo[:, [joinCol, labelCol]]
    dfLabelInfo = join(dfLabelCoordinates,
                       subInfo, on = joinCol)

    ## convert possible symbols in label column to strings
    if !(issubtype(eltype(dfLabelInfo[labelCol]), String))
        ## conversion required for join
        dfLabelInfo[labelCol] = AssetMgmt.symbToStr(dfLabelInfo[labelCol])
    end

    ## remove rows with empty label
    noLabel = Bool[dfLabelInfo[labelCol][ii] == "" for ii=1:size(dfLabelInfo, 1)]
    dfLabelInfo = dfLabelInfo[!noLabel, :]

    # rename column name
    rename!(stackedData, :variable, joinCol)

    # join asset class
    subInfo = variableInfo[:, [joinCol, colorCol]]
    xx = join(stackedData, subInfo, on = joinCol)

    p = plot(layer(xx, x="Idx", y="value",
                   group=joinCol, color=colorCol,
                   Geom.line),
             layer(dfLabelInfo, x="x", y="y",
                   label=labelCol,
                   Gadfly.Geom.label(position=:dynamic,
                                     hide_overlaps=false)),
             args...)
    return p
end

#################################
## visualize mu-sigma universe ##
#################################

function plotAssetMoments(moments::DataFrame, args...;
                          variableInfo::DataFrame = DataFrame(),
                          scalingFunc::Function = AssetMgmt.defaultMuSigmaScaling,
                          joinCol::Symbol = :AssetLabel,
                          colorCol::Symbol = :AssetClass,
                          labelCol::Symbol = :AssetLabel,
                          shiftLabels::Float64 = 0.)

    ## Visualize return moments.
    ## Return moments are scaled according to some given function.
    ## Also, points can be colored with respect to additional
    ## information.

    ## ## check for extended asset info
    ## extendedInfo = true
    ## if isempty(assetInfo)
    ##     extendedInfo = false
    ## end

    ## convert keys in moments table to strings
    if eltype(moments[joinCol]) == Symbol
        ## conversion required for join
        moments = deepcopy(moments)
        moments[joinCol] = AssetMgmt.symbToStr(moments[joinCol])
    end

    ## convert keys in asset info table to strings
    if eltype(variableInfo[joinCol]) == Symbol
        ## conversion required for join
        variableInfo = deepcopy(variableInfo)
        variableInfo[joinCol] = AssetMgmt.symbToStr(variableInfo[joinCol])
    end

    ## get number of asset
    nAss = size(moments, 1)

    ## scale moments
    momentsArr = convert(Array, moments[:, [:mu, :sigma]])
    scaledMus, scaledSigmas = scalingFunc(momentsArr[:, 1],
                                          momentsArr[:, 2])

    scaledMoments = DataFrame()
    scaledMoments[joinCol] = moments[joinCol]
    scaledMoments[:mu] = scaledMus
    scaledMoments[:sigma] = scaledSigmas

    ## join assetInfo table
    dfLabelInfo = join(scaledMoments,
                       variableInfo, on = joinCol)
    momentsTable = join(scaledMoments,
                        variableInfo, on = joinCol)

    ## convert possible symbols in label column to strings
    if !(issubtype(eltype(dfLabelInfo[labelCol]), String))
        ## conversion required for join
        dfLabelInfo[labelCol] = AssetMgmt.symbToStr(dfLabelInfo[labelCol])
    end

    ## remove rows with empty label
    noLabel = Bool[dfLabelInfo[labelCol][ii] == "" for ii=1:size(dfLabelInfo, 1)]
    dfLabelInfo = dfLabelInfo[!noLabel, :]

    p = plot(layer(momentsTable,
                   x="sigma", y="mu",
                   color=colorCol,
                   Geom.point),
             layer(dfLabelInfo, x="sigma", y="mu",
                   label=labelCol,
                   Gadfly.Geom.label(position=:dynamic,
                                     hide_overlaps=false)),
             Guide.title("μ-σ diagram"),
             Guide.xlabel("σ"), Guide.ylabel("μ"),
             args...)
    return p
end

function plotAssetMoments(mod::AssetMgmt.MuSigmaModel, args...;
                          scalingFunc::Function =
                          AssetMgmt.defaultMuSigmaScaling)
    ## get information
    moments = DataFrame()
    moments[:AssetLabel] = mod.names
    moments[:mu] = mod.mu
    moments[:sigma] = AssetMgmt.getVolas(mod)
    p = plotAssetMoments(moments, args...;
                         variableInfo = moments,
                         scalingFunc = scalingFunc,
                         colorCol = :AssetLabel,
                         labelCol = :AssetLabel)
    return p
end

## plot efficient frontier
##------------------------

function plotEff(mod::MuSigmaModel;
                 scalingFunc = AssetMgmt.defaultMuSigmaScaling,
                 nPoints = 50)

    ## extract and scale moments
    sigmas = AssetMgmt.getVolas(mod)
    scaledMus, scaledSigmas = scalingFunc(mod.mu, sigmas)

    ## calculate efficient frontier moments
    muEffUnscaled, sigmaEffUnscaled =
        getEffFrontier(mod, nPoints = nPoints)

    ## scale efficient frontier
    muEff, sigmaEff =
        scalingFunc(muEffUnscaled, sigmaEffUnscaled)

    p = plot(layer(x=sigmaEff, y=muEff, Geom.line),
             layer(x=scaledSigmas, y=scaledMus, Geom.point));

    return p
end

## plot asset and portfolio moments
##---------------------------------

function plotPfsAndUniverse(mod::MuSigmaModel,
                            wgts = Array{Float64, 2};
                            scalingFunc::Function = AssetMgmt.defaultMuSigmaScaling)
    ## portfolios should be given in rows

    ## extract and scale moments
    sigmas = AssetMgmt.getVolas(mod)
    scaledMus, scaledSigmas = scalingFunc(mod.mu, sigmas)

    ## get type names
    typNams = UTF8String["asset" for ii=1:length(sigmas)]

    ## get asset moments as DataFrame
    assetMomentsTable = DataFrame(mu = scaledMus,
                                  sigma = scaledSigmas,
                                  color = typNams)

    ## get portfolio moments
    muPfs = AssetMgmt.getPMean(wgts, mod.mu)
    sigmaPfs = sqrt(AssetMgmt.getPVar(wgts, mod.sigma))

    ## scale portfolio moments
    scaledMuPfs, scaledSigmaPfs =
        scalingFunc(muPfs, sigmaPfs)

    ## get type names
    typNams = UTF8String["portfolio" for ii=1:length(sigmaPfs)]

    ## get portfolio moments as DataFrame
    pfMomentsTable = DataFrame(mu = scaledMuPfs,
                               sigma = scaledSigmaPfs,
                               color = typNams)

    ## join tables
    momentsTable = [assetMomentsTable; pfMomentsTable]

    p = plot(momentsTable, x="sigma", y="mu",
             color="color",
             Geom.point, Guide.title("μ-σ diagram"),
             Guide.xlabel("σ"), Guide.ylabel("μ"));

    return p
end

import Gadfly.plot
function plot(tm::Timematr, settings...)
    ## plot Timematr object
    ## - multi-line plot
    ## - index must be dates convertible to strings
    ## - plot only some maximum number of randomly picked paths
    ## - plot some statistics

    dats = TimeData.dat2num(tm)
    (nObs, nAss) = size(tm)

    ## maximum number of paths plotted
    maxPaths = 10

    ## get values for paths
    if nAss > maxPaths
        ## get randomly drawn paths

        chosenAssets = randperm(nAss)[1:maxPaths]
        randPaths = core(tm)[:, chosenAssets]
        nams = names(tm)[chosenAssets]

        ## get statistics
        allVals = core(tm)
        averagePath = mean(allVals, 2)
        lowestEndingPath = indmin(allVals[end, :])
        lowPath = allVals[:, lowestEndingPath]
        highestEndingPath = indmax(allVals[end, :])
        highPath = allVals[:, highestEndingPath]

        vals = [averagePath lowPath highPath randPaths]
        valsDf = AssetMgmt.composeDataFrame(vals,
                                            [:mean, :low, :high,
                                             nams])
    else
        valsDf = copy(tm.vals)
    end

    nams = names(valsDf)
    valsDf[:Idx] = dats

    stackedData = stack(valsDf, nams)

    plot(stackedData, x="Idx", y="value",
         color="variable", Geom.line, settings...)
end

function plotPriceReturnSeries(pfRets::Timematr,
                               discRet::Timematr,
                               BMs...;
                               strName = :strategy)

    ## plot price evolution of assets and portfolios
    ##
    ## Inputs:
    ## 	pfRets	nObs x 1 Timematr of portfolio returns
    ## 	discRet 	nObs x nAss Timematr of discrete asset returns
    ## 	BMs		nBms tuple of benchmark returns
    ## 	strName	optional expression to label strategy in plots

    (nObs, nAss) = size(discRet)
    nBMs = length(BMs)

    ## get cumulated prices for individual assets
    cumPrices = cumprod(discRet .+ 1, 1)

    ############################
    ## asset price evolutions ##
    ############################

    pricePlot = AssetMgmt.plot(cumPrices, Guide.xlabel("time"),
                               Guide.ylabel("price"),
                               Guide.title("normalized asset prices"));

    ################################################
    ## portfolio price evolutions with benchmarks ##
    ################################################

    ## plot portfolio evolution with benchmarks
    bmPlot = AssetMgmt.plotWithBMs(pfRets, BMs...);

    #################################
    ## plot returns and volatility ##
    #################################

    returnPlot = AssetMgmt.plot(pfRets, Guide.xlabel("time"),
                                Guide.ylabel("return"),
                                Guide.title("portfolio returns"));

    ## get squared portfolio returns as proxy for volatility
    volaPlot = AssetMgmt.plot(pfRets.^2, Guide.xlabel("time"),
                              Guide.ylabel("squared returns"),
                              Guide.title("squared portfolio
returns"));

    ################################
    ## putting the parts together ##
    ################################

    pfEvolution = vstack(pricePlot, bmPlot, returnPlot, volaPlot)
    filename = string("dev_pics/pfEvolution_", strName, ".pdf")
    draw(PDF(filename, 10inch, 20inch), pfEvolution)
end

function plotTOverDiversification(rawInvs::Investments,
                                  fltInvs::Investments,
                                  discRet::Timematr;
                                  ## sectors::Dict;
                                  strName = :strategy)

    ## plot intended and realized turnover and diversification
    ## measures
    ##
    ## Inputs:
    ## 	rawInvs	nObs x nAss unfiltered Investments
    ## 	fltInvs	nObs x nAss turnover reducing Investments
    ## 	discRet	nObs x nAss Timematr of discrete asset returns
    ## 	sectors	nSect x 1 Dictionary of assets per sector
    ## 	strName	optional: expression with strategy name for plot
    ## 				labeling

    ##############
    ## turnover ##
    ##############

    ## get investment turnover
    tOver = AssetMgmt.turnover(fltInvs, discRet)
    intTOver = AssetMgmt.turnover(rawInvs, discRet)

    tOData = [tOver intTOver]

    tOPlot = AssetMgmt.plot(tOData, Guide.xlabel("time"),
               Guide.ylabel("turnover"),
               Guide.title("actual turnover vs intended turnover"));

    ##############################
    ## diversification measures ##
    ##############################

    ## get diversification measures
    divIndicators = AssetMgmt.diversification(fltInvs)

    ## plot diversification measures related to number of assets
    divNAssets = divIndicators[[:nAss, :nSignAss, :nShort]]
    divNAssetsPlot = AssetMgmt.plot(divNAssets, Guide.xlabel("time"),
                                    Guide.ylabel("number"),
                                    Guide.title("diversification measures"));

    ## plot diversification measures expressed as fractions
    divFractions = divIndicators[[:maxWgt, :Gini, :stdWgt]]
    divFractionsPlot = AssetMgmt.plot(divFractions, Guide.xlabel("time"),
               Guide.ylabel("value"),
               Guide.title("diversification measures"));

    divPlot = vstack(tOPlot, divNAssetsPlot, divFractionsPlot)
    filename = string("dev_pics/pfDiversif_", strName, ".pdf")
    draw(PDF(filename, 10inch, 12inch), divPlot)
end

function plotPfDistribution(pfRet::Timematr,
                            discRet::Timematr,
                            BMs...;
                            strName = :strategy)
    ## plot strategy performance indicators

    (nObs, nAss) = size(discRet)
    nBMs = length(BMs)

    #################################
    ## portfolio return statistics ##
    #################################

    retStats = AssetMgmt.returnStatistics(pfRet)

    distrPlots = Array(Any, (nBMs+1))

    distrPlots[1] = plot(pfRet.vals, x=string(strName),
         xintercept = [retStats.mu, retStats.VaR95, retStats.VaR99],
         Geom.histogram(minbincount=80),
         Geom.vline(color="orange"),
         Guide.xlabel("returns"),
         Guide.ylabel("frequency"),
         Guide.title("portfolio return distribution"));

    for ii=1:nBMs
        ## currBM = Timematr(BMs[ii].vals, idx(BMs[ii]))
        currBM = BMs[ii]
        title = string(string(names(currBM)[1]), " return distribution")
        retStats = AssetMgmt.returnStatistics(currBM)
        distrPlots[ii+1] = plot(currBM.vals, x=string(names(currBM.vals)[1]),
                             xintercept = [retStats.mu, retStats.VaR95, retStats.VaR99],
                             Geom.histogram(minbincount=80),
                             Geom.vline(color="orange"),
                             Guide.xlabel("returns"),
                             Guide.ylabel("frequency"),
                             Guide.title(title));
    end

    ###################
    ## mu-sigma-plot ##
    ###################

    retStats = AssetMgmt.returnStatistics(pfRet)

    ## get moments of assets
    assetMus = mean(core(discRet), 1)
    assetSigmas = std(core(discRet), 1)
    ## covs = cov(discRet)
    ## for ii=1:nAss
    ##     assetSigmas[names(discRet)[ii]] = sqrt(covs[ii, ii])
    ## end
    assetNames = Array(ASCIIString, nAss)
    for ii=1:nAss
        assetNames[ii] = string(names(discRet)[ii])
    end
    assetTyps = ["asset" for ii=1:nAss]

    ## get moments for benchmark strategies
    bmMus = Array(Float64, nBMs)
    bmSigmas = Array(Float64, nBMs)
    bmNames = Array(ASCIIString, nBMs)
    for ii=1:nBMs
        currBM = BMs[ii]
        bmMus[ii] = mean(core(currBM))
        bmSigmas[ii] = std(core(currBM))
        bmNames[ii] = string(names(currBM)[1])
    end
    bmTyps = ["Benchm" for ii=1:nBMs]

    ## get moments for portfolio return
    pfMu = retStats.mu
    pfSigma = retStats.sigma
    pfName = string(strName)
    pfTyp = "portf"

    moments = DataFrame()
    ## moments[:mus] = [assetMus[:], bmMus[:], pfMu]
    ## moments[:sigmas] = [assetSigmas[:], bmSigmas[:], pfSigma]
    ## moments[:names] = [assetNames[:], bmNames[:], pfName]
    ## moments[:typ] = [assetTyps[:], bmTyps[:], pfTyp]

    ## for debugging
    ## display(moments)

    ## new version should prioritize labels of portfolios over assets
    moments[:mus] = [bmMus[:], pfMu, assetMus[:]]
    moments[:sigmas] = [bmSigmas[:], pfSigma, assetSigmas[:]]
    moments[:names] = [bmNames[:], pfName, assetNames[:]]
    moments[:typ] = [bmTyps[:], pfTyp, assetTyps[:]]

    muSigmaPlot = plot(layer(moments,
                             x="sigmas", y="mus", color="typ",
                             Geom.point),
                       layer(moments, x="sigmas", y="mus",
                             label="names", Geom.label),
                       Guide.xlabel("sigma"),
                       Guide.ylabel("mu"),
                       Guide.title("mu-sigma combinations"))

    distrPlots = vstack(muSigmaPlot, distrPlots...)
    filename = string("dev_pics/pfDistrib_", strName, ".pdf")
    draw(PDF(filename, 10inch, 25inch), distrPlots)

end

function plotSectorAnalysis(fltInvs::Investments,
                            discRet::Timematr,
                            assetsInSector::Dict;
                            strName = :strategy)

        ## get mean sector return series
    sectorReturns = AssetMgmt.groupVars(discRet,
                                        assetsInSector,
                                        x -> mean(x, 2))

    ## get sector means and vars
    mus = mean(core(sectorReturns)*100, 1)
    sigmas = std(core(sectorReturns)*100, 1)

    ## get sector weights
    sectorWgts = mean(core(AssetMgmt.groupAssets(fltInvs, assetsInSector)),
                      1)
    absWgts = abs(sectorWgts)
    longAssets = 1*(sectorWgts .> 0)

    df = DataFrame()
    nams = [string(names(sectorReturns)[ii]) for ii=1:length(mus)]
    df[:names] = nams
    df[:mus] = mus[:]
    df[:sigmas] = sigmas[:]
    df[:absWgts] = absWgts[:]
    df[:weights] = sectorWgts[:]
    df[:long] = longAssets[:]

    muSigmaPlot = plot(layer(df, x="sigmas",
                             y="mus",
                             color="weights", Geom.point),
                       layer(df, x="sigmas", y="mus",
                             label="names", Geom.label),
                       Guide.xlabel("sigma"),
                       Guide.ylabel("mu"),
                       Guide.title("sector weights with sector moments")
                       )

    ## get weights for each sector
    sectorWgts = AssetMgmt.groupAssets(fltInvs, assetsInSector)

    title = string(string(strName), " portfolio sector weights")
    sectorWgtsPlot = AssetMgmt.plot(convert(Timematr, sectorWgts),
                                    Guide.xlabel("time"),
                                    Guide.ylabel("sector weight"),
                                    Guide.title(title))

    scaledSectorWgts = DataFrame()
    for eachSector in assetsInSector
        nAssPerSector = length(eachSector[2])
        sectorName = eachSector[1]
        scaledSectorWgts[sectorName] =
            (core(sectorWgts[sectorName]) ./ nAssPerSector)[:]
    end
    scaledSectorWgtsTm = Timematr(scaledSectorWgts, idx(fltInvs))

    title = string(string(strName), " average asset weight per sector")
    averageWgtsPlot = AssetMgmt.plot(scaledSectorWgtsTm,
                                     Guide.xlabel("time"),
                                     Guide.ylabel("average weight"),
                                     Guide.title(title))

    sectPlots = vstack(muSigmaPlot, sectorWgtsPlot, averageWgtsPlot)
    filename = string("dev_pics/pfSectors_", strName, ".pdf")
    draw(PDF(filename, 10inch, 20inch), sectPlots)

end


text:
module CommunalHelper

using ..Ahorn, Maple
using Cairo

# General
export hexToRGBA
# Dream Blocks
export renderDreamBlock
# Cassette Blocks
export renderCassetteBlock, cassetteColorNames, getCassetteColor
# Connected Blocks
export getExtensionRectangles, notAdjacent
# Depths
export depths
# mapdef hacks
export @mapdefdata, appendkwargs

"Celeste.Depths"
const depths = Dict{String,Integer}(
    "BGTerrain (10000)" => 10000,
    "BGDecals (9000)" => 9000,
    "BGParticles (8000)" => 8000,
    "Below (2000)" => 2000,
    "NPCs (1000)" => 1000,
    "Player (0)" => 0,
    "Dust (-50)" => -50,
    "Pickups (-100)" => -100,
    "Particles (-8000)" => -8000,
    "Above Particles (-8500)" => -8500,
    "Solids (-9000)" => -9000,
    "FGTerrain (-10000)" => -10000,
    "FGDecals (-10500)" => -10500,
    "DreamBlocks (-11000)" => -11000,
    "CrystalSpinners (-11500)" => -11500,
    "Chaser (-12500)" => -12500,
    "Fake Walls (-13000)" => -13000,
    "FGParticles (-50000)" => -50000,
)

"Celeste.SurfaceSound"
const surfaceSounds = Dict{String,Integer}(
    "Asphalt" => 1,
    "Car" => 2,
    "Dirt" => 3,
    "Snow" => 4,
    "Wood" => 5,
    "StoneBridge" => 6,
    "Girder" => 7,
    "Brick" => 8,
    "ZipMover" => 9,
    "DreamBlockInactive" => 11,
    "DreamBlockActive" => 12,
    "ResortWood" => 13,
    "ResortRoof" => 14,
    "ResortSinkingPlatforms" => 15,
    "ResortBasementTile" => 16,
    "ResortLinens" => 17,
    "ResortBoxes" => 18,
    "ResortBooks" => 19,
    "ClutterDoor" => 20,
    "ClutterSwitch" => 21,
    "ResortMagicButton" => 21,
    "ResortElevator" => 22,
    "CliffsideSnow" => 23,
    "CliffsideGrass" => 25,
    "CliffsideWhiteBlock" => 27,
    "Gondola" => 28,
    "AuroraGlass" => 32,
    "Grass" => 33,
    "CassetteBlock" => 35,
    "CoreIce" => 36,
    "CoreMoltenRock" => 37,
    "Glitch" => 40,
    "MoonCafe" => 42,
    "DreamClouds" => 43,
    "Moon" => 44,
)

"""
    hexToRGBA(hex)

Convert a hex color code to an RGBA tuple with value from 0.0-1.0

# Examples:
    hexToRGBA("ff00ff")
"""
hexToRGBA(hex) =
    tuple(Ahorn.argb32ToRGBATuple(parse(Int, hex; base=16))[1:3] ./ 255..., 1.0)

"""
    renderDreamBlock(ctx, x, y, width, height, data)

Draw a custom dreamblock based on the attributes present in `data`.

# Supported Attributes:
    `featherMode`::Bool
    `oneUse`::Bool
    `noCollide`::Bool
    `doubleRefill`::Bool OR `refillCount`::Int
"""
function renderDreamBlock(ctx::CairoContext, x::Number, y::Number, width::Number, height::Number, data::Dict{String, Any})
    save(ctx)

    set_antialias(ctx, 1)
    set_line_width(ctx, 1)

    fillColor =
        get(data, "featherMode", false) ? (0.31, 0.69, 1.0, 0.4) : (0.0, 0.0, 0.0, 0.4)
    #get(data, "noCollide", false)

    lineColor = (1.0, 1.0, 1.0, 1.0)
    if get(data, "doubleRefill", false)
        lineColor = (1.0, 0.43, 0.94, 1.0)
    else
        refillCount = Int(get(data, "refillCount", -1))
        if refillCount != -1
            lineColor = get(hairColors, refillCount + 1, (1.0, 0.43, 0.94, 1.0))
        end
    end

    if (get(data, "oneUse", false))
        set_dash(ctx, [0.6, 0.2])
    end

    Ahorn.drawRectangle(ctx, x, y, width, height, fillColor, lineColor)

    restore(ctx)
end

# Translated from Celeste.UserIO.GetSavePath
function getSavesDir()
    if Sys.islinux() || Sys.isfreebsd() || Sys.isopenbsd() || Sys.isnetbsd()
        envVar = get(ENV, "XDG_DATA_HOME", nothing)
        if !isnothing(envVar) && !isempty(envVar)
            return joinpath(envVar, "Celeste", "Saves")
        end
        envVar = get(ENV, "HOME", nothing)
        if !isnothing(envVar) && !isempty(envVar)
            return joinpath(envVar, ".local", "share", "Celeste", "Saves")
        end
    elseif Sys.isapple()
        envVar = get(ENV, "HOME", nothing)
        if !isnothing(envVar) && !isempty(envVar)
            return joinpath(envVar, "Library", "Application Support", "Celeste", "Saves")
        end
    elseif Sys.iswindows()
        return joinpath(Ahorn.config["celeste_dir"], "Saves")
    end
    throw(
        ErrorException(
            "Unsupported operating system. (how did you even get celeste working?)",
        ),
    )
end

function getPlayerHairColors()
    colors = [(0.27, 0.7, 1.0, 1.0), (0.67, 0.2, 0.2, 1.0), (1.0, 0.43, 0.94, 1.0)]
    try
        path = joinpath(getSavesDir(), "modsettings-MoreDasheline.celeste")
        if isfile(path)
            # can't hecking use YAML.jl because it can't read the hex values properly
            moreDasheline = Dict{String,String}(
                strip(key) => strip(value) for (key, value) in split.(readlines(path), ":")
            )

            colors = append!(colors, [
                    hexToRGBA(moreDasheline["ThreeDashColor"]),
                    hexToRGBA(moreDasheline["FourDashColor"]),
                    hexToRGBA(moreDasheline["FiveDashColor"]),
                ],
            )
            println("Using MoreDasheline hair colors for CommunalHelper DreamBlocks.")
            return colors
        end
    catch err
        @warn "Error loading MoreDasheline hair colors:\n$err"
    end

    println(
        "MoreDasheline hair colors not loaded, using default hair colors for CommunalHelper DreamBlocks.",
    )

    return colors
end

const hairColors = getPlayerHairColors()

const cassetteBlock = "objects/cassetteblock/solid"
const cassetteColors = Dict{Int,Ahorn.colorTupleType}(
    1 => (240, 73, 190, 255) ./ 255,
    2 => (252, 220, 58, 255) ./ 255,
    3 => (56, 224, 78, 255) ./ 255,
)
const defaultCassetteColor = (73, 170, 240, 255) ./ 255

const cassetteColorNames =
    Dict{String,Int}("Blue" => 0, "Rose" => 1, "Bright Sun" => 2, "Malachite" => 3)

getCassetteColor(index::Int) = get(cassetteColors, index, defaultCassetteColor)

function renderCassetteBlock(ctx::CairoContext, x, y, width, height, index, customColor=nothing)
    tileWidth = ceil(Int, width / 8)
    tileHeight = ceil(Int, height / 8)

    color = isnothing(customColor) ? get(cassetteColors, index, defaultCassetteColor) : customColor

    for i in 1:tileWidth, j in 1:tileHeight
        tx = (i == 1) ? 0 : ((i == tileWidth) ? 16 : 8)
        ty = (j == 1) ? 0 : ((j == tileHeight) ? 16 : 8)

        Ahorn.drawImage(ctx, cassetteBlock, x + (i - 1) * 8, y + (j - 1) * 8, tx, ty, 8, 8, tint=color)
    end
end

# Get Rectangles from SolidExtensions present in the room.
function getExtensionRectangles(room::Room)
    entities = filter(e -> e.name == "CommunalHelper/SolidExtension", room.entities)
    rects = []

    for e in entities
        push!(
            rects,
            Ahorn.Rectangle(
                Int(get(e.data, "x", 0)),
                Int(get(e.data, "y", 0)),
                Int(get(e.data, "width", 8)),
                Int(get(e.data, "height", 8)),
            ),
        )
    end

    return rects
end

"Check for collision with an array of rectangles at specified tile position"
function notAdjacent(x, y, ox, oy, rects)
    rect = Ahorn.Rectangle(x + ox + 4, y + oy + 4, 1, 1)

    for r in rects
        if Ahorn.checkCollision(r, rect)
            return false
        end
    end

    return true
end
notAdjacent(entity::Entity, ox, oy, rects) = notAdjacent(Ahorn.position(entity)..., ox, oy, rects)

function detectMod(mod)
    any(s -> occursin(mod, lowercase(s)), Ahorn.getCelesteModZips()) ||
        any(s -> occursin(mod, lowercase(s)), Ahorn.getCelesteModDirs())
end

# Don't worry about it
macro mapdefdata(entity::Symbol, id::String, type::Symbol, data::Symbol)
    expr = Base.eval(__module__, data)
    expr.args[1] = type
    return :(@mapdef $entity $id $expr) |> esc
end

function appendkwargs(func::Expr, kwargs::Expr...)
    return _appendkwargs!(deepcopy(func), deepcopy(kwargs)...)
end

function _appendkwargs!(func::Expr, kwargs::Expr...)
    for e in kwargs
        if e.head == :tuple
            _appendkwargs!(func, e.args...)
        else
            e.head = :kw
            push!(func.args, e)
        end
    end
    return func
end

end

ext: jl


text:
using Mamba

## Data
blocker = Dict{Symbol, Any}(
  :rt =>
    [3, 7, 5, 102, 28, 4, 98, 60, 25, 138, 64, 45, 9, 57, 25, 33, 28, 8, 6, 32,
     27, 22],
  :nt =>
    [38, 114, 69, 1533, 355, 59, 945, 632, 278, 1916, 873, 263, 291, 858, 154,
     207, 251, 151, 174, 209, 391, 680],
  :rc =>
    [3, 14, 11, 127, 27, 6, 152, 48, 37, 188, 52, 47, 16, 45, 31, 38, 12, 6, 3,
     40, 43, 39],
  :nc =>
    [39, 116, 93, 1520, 365, 52, 939, 471, 282, 1921, 583, 266, 293, 883, 147,
     213, 122, 154, 134, 218, 364, 674]
)
blocker[:N] = length(blocker[:rt])

## Model Specification
model = Model(

  rc = Stochastic(1,
    (mu, nc, N) ->
      begin
        pc = invlogit(mu)
        UnivariateDistribution[Binomial(nc[i], pc[i]) for i in 1:N]
      end,
    false
  ),

  rt = Stochastic(1,
    (mu, delta, nt, N) ->
      begin
        pt = invlogit(mu + delta)
        UnivariateDistribution[Binomial(nt[i], pt[i]) for i in 1:N]
      end,
    false
  ),

  mu = Stochastic(1,
    () -> Normal(0, 1000),
    false
  ),

  delta = Stochastic(1,
    (d, s2) -> Normal(d, sqrt(s2)),
    false
  ),

  delta_new = Stochastic(
    (d, s2) -> Normal(d, sqrt(s2))
  ),

  d = Stochastic(
    () -> Normal(0, 1000)
  ),

  s2 = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  )

)

## Initial Values
inits = [
  Dict(:rc => blocker[:rc], :rt => blocker[:rt], :d => 0, :delta_new => 0,
       :s2 => 1, :mu => zeros(blocker[:N]), :delta => zeros(blocker[:N])),
  Dict(:rc => blocker[:rc], :rt => blocker[:rt], :d => 2, :delta_new => 2,
       :s2 => 10, :mu => fill(2, blocker[:N]), :delta => fill(2, blocker[:N]))
]

## Sampling Scheme
scheme = [AMWG(:mu, 0.1),
          AMWG([:delta, :delta_new], 0.1),
          Slice([:d, :s2], 1.0)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, blocker, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)


text:
#!/usr/bin/env julia

import HTTP, JSON, OrderedCollections

vars = ["varno","codetype", "obstype", "sensor", "reporttype", "subtype", "obstype"]

for var in vars
  println("Creating $var.schema.json")
  href = "https://apps.ecmwf.int/odbgov/$var"
  resp = HTTP.get("$href.json")
  dict = JSON.parse(String(resp.body))

  findind(field) = findfirst(isequal(field), dict["fields"])

  data = sort(dict["data"], lt= (r1, r2) -> r1[findind("code")] < r2[findind("code")])

  code        = get.(data, findind("code"),nothing)
  description = get.(data, findind("description"),nothing)

  enum_titles = ["$(lpad(c,3)), $d" for (c,d) in zip(code, description)]

  r = OrderedCollections.OrderedDict()
  r[raw"$comment"] = "File created by create_odbschemas.jl"
  r["type"] = "integer"
  r["title"] = var
  r["enum"]  = code
  r["options"] = Dict("enum_titles"=>enum_titles)
  r["links"]   = [Dict("rel"=>"ECMWF documentation","href"=>href)]

  io = open("$var.schema.json","w")
  JSON.print(io, r, 4)
  close(io)
end


text:
for ijulia_behavior in [:warn, :clear, :append]

ProgressMeter.ijulia_behavior(ijulia_behavior)

println("Testing showvalues with a Dict (2 values)")
function testfunc1(n, dt, tsleep, desc, barlen)
    p = ProgressMeter.Progress(n, dt, desc, barlen)
    for i = 1:n
        sleep(tsleep)
        ProgressMeter.next!(p; showvalues = Dict(:i => i, :halfdone => (i >= n/2)))
    end
end
testfunc1(50, 1, 0.2, "progress  ", 70)

println("Testing showvalues with an Array of tuples (4 values)")
function testfunc2(n, dt, tsleep, desc, barlen)
    p = ProgressMeter.Progress(n, dt, desc, barlen)
    for i = 1:n
        sleep(tsleep)
        ProgressMeter.next!(p; showvalues = [(:i, i),
            (:constant, "foo"), (:isq, i^2), (:large, 2^i)])
    end
end
testfunc2(30, 1, 0.2, "progress  ", 60)

println("Testing showvalues when types of names differ (3 values)")
function testfunc3(n, dt, tsleep, desc, barlen)
    p = ProgressMeter.Progress(n, dt, desc, barlen)
    for i = 1:n
        sleep(tsleep)
        ProgressMeter.next!(p; showvalues = [(:i, i*10), ("constant", "foo"),
            ("foobar", round(i*tsleep, digits=4))])
    end
end
testfunc3(30, 1, 0.2, "progress  ", 70)

println("Testing progress with showing values when num values to print changes between iterations")
function testfunc4(n, dt, tsleep, desc, barlen)
    p = ProgressMeter.Progress(n, dt, desc, barlen)
    for i = 1:n
        sleep(tsleep)
        values = [(:i, i*10), ("constant", "foo"), ("foobar", round(i*tsleep, digits=4))]
        ProgressMeter.next!(p; showvalues = values[randn(3) .< 0.5])
    end
end
testfunc4(30, 1, 0.2, "opt steps  ", 70)

println("Testing showvalues with a different color (1 value)")
function testfunc5(n, dt, tsleep, desc, barlen)
    p = ProgressMeter.Progress(n, dt, desc, barlen)
    for i = 1:n
        sleep(tsleep)
        ProgressMeter.next!(p; showvalues = [(:large, 2^i)], valuecolor = :yellow)
    end
end
testfunc5(10, 1, 0.2, "progress  ", 40)

println("Testing showvalues with threshold-based progress")
prog = ProgressMeter.ProgressThresh(1e-5, "Minimizing:")
for val in 10 .^ range(2, stop=-6, length=20)
    ProgressMeter.update!(prog, val; showvalues = Dict(:margin => abs(val - 1e-5)))
    sleep(0.1)
end

println("Testing showvalues with online progress")
prog = ProgressMeter.ProgressUnknown("Entries read:")
for title in ["a", "b", "c", "d", "e"]
    ProgressMeter.next!(prog; showvalues = Dict(:title => title))
    sleep(0.5)
end
ProgressMeter.finish!(prog)

println("Testing showvalues with early cancel")
prog = ProgressMeter.Progress(100, 1, "progress: ", 70)
for i in 1:50
    ProgressMeter.update!(prog, i; showvalues = Dict(:left => 100 - i))
    sleep(0.1)
end
ProgressMeter.cancel(prog)

end # if


text:
"""

In this example, we compute the time correlation function of the
rotations of the side chain of a trp residue

"""

using M3GTools

println(" Loading Plots... ")
using Plots
ENV["GKSwstype"]="nul" # This supresses the need of a display while plotting

# Initialize simulation data and path to VMD executable
println(" Reading simulation data ... ")
mysim = M3GTools.init(psf="../simulation_files/structure.psf",
                  dcd="../simulation_files/structure.dcd",
                  vmd="vmd")

println(" Defining selections... ")

# Absorption vector, start and end (the center of mass of selected atoms is considered)

abs_start = M3GTools.select(mysim,"index 1548 1549")
abs_end = M3GTools.select(mysim,"index 1546")

# Emission vector, start and end
emi_start = M3GTools.select(mysim,"index 1548 1549")
emi_end = M3GTools.select(mysim,"index 1544")

# Optional: align coordinates to first frame, according to the alignment of
# the selection defined here. In this case, we remove the rotation of the protein.

CAs = M3GTools.select(mysim,"protein and name CA")

t, legendre, tcf = M3GTools.tcf(mysim,
                            abs_start,abs_end,emi_start,emi_end;
                            r0=0.4,
                            theta=34.24,
                            #maxdt=20.,
                            scaletime=0.001,
                            align=CAs)

# It is a good idea to save the data to a file:

file="tcf.dat"
println(" Writting data to file: ", file)
f = open(file,"w")
println(f," Time  2nd-Legendre TCF ")
for i in 1:length(t)
  println(f,t[i]," ",legendre[i]," ",tcf[i])
end
close(f)

println(" Plotting... ")
plot(t,legendre)
savefig("tcf.pdf")
println(" Created figure: tcf.pdf")


text:
function _io_print(print_fun, headline, msg, dat, file, more...)
    print_fun(headline,
        isempty(msg) ? "" : string("\n", msg),
        "\ndir: ", relpath(dirname(abspath(file))),
        "\nfile: ", basename(file),
        "\nsize: ", filesize(file), " bytes",
        "\ndata type: ", typeof(dat),
        more...,
        "\n"
    )
end

function _io_error_print(print_fun, err, file, more...)
    print_fun("ERROR",
        "\ndir: ", relpath(dirname(abspath(file))),
        "\nfile: ", basename(file),
        "\nerr: ", err_str(err),
        more...,
        "\n"
    )
end


text:
### Functions for loading/anntoting mutations/variants

using FASTX

### load mutations

function loadmutations(; file="hk1mutations.tsv", projdir=getprojectdir())
    filepath = joinpath(projdir, "data", file)
    df = CSV.read(filepath, DataFrame)

    df.Ref = [ismissing(s) ? dna"" : LongDNASeq(s) for s in df.Ref]
    df.Alt = [ismissing(s) ? dna"" : LongDNASeq(s) for s in df.Alt]

    df.PatientGroup = Int.(floor.(df.Patient))
    df[!, :ID] = string.(df.start, "_", df.Ref, "_", df.Alt)

    df
end

function parsechromloc(s)
    fields = split(replace(s, "," => ""), r"[:-]")
    fields[1], parse(Int, fields[2]):parse(Int, fields[3])
end

function loadregregion(;file="hk1_reg_region.fa",  projdir=getprojectdir())
    filepath = joinpath(projdir, "data", file)
    reader = open(FASTA.Reader, filepath)
    record = first(collect(reader))
    seq = sequence(record)
    chromloc = replace(first(split(description(record))), "range=" => "")
    chrom, loc = parsechromloc(chromloc)
    seq
    Interval(chrom, loc, '+', seq)
end

function annotateregion!(mutations, region)

    regc = leftposition(region):rightposition(region)
    mutreg = [s:e for (s, e) in zip(mutations.start, mutations.stop)]
    relreg = [mr .- leftposition(region) .+ 1 for mr in mutreg]
    mutations[!, :RelReg] = relreg
    mutations[!, :Contained] = length.(intersect.(Ref(regc), mutreg)) .== length.(mutreg)

    ### assert that ref alleles match sequence
    refseqs = [region.metadata[r] for r in mutations[mutations.Contained, :RelReg]]
    @assert refseqs == mutations.Ref[mutations.Contained]
    mutations

end

buildmutation(refreg, kind, ref, alt, contained) = (coords=refreg, kind=kind, ref=ref, alt=alt, contained=contained)

function buildaltseq(seq, mutation)
    if !mutation.contained
        return missing, missing, missing, missing
    end
    leftind = 1:(first(mutation.coords) - 1)
    refrightind = (last(mutation.coords) + 1):length(seq)
    if mutation.kind == "var"
        altind  = mutation.coords
        altrightind = refrightind
    elseif mutation.kind == "del"
        altind = first(mutation.coords) .+ (0:-1)
        altrightind = first(mutation.coords):(length(seq) - length(mutation.coords))
    else
        error("Mutation type: $mutation.kind not recognised")
    end

    altseq = seq[leftind]*mutation.alt*seq[refrightind]
    leftind, altind, refrightind, altrightind, altseq
end

function distinctmuts(mutations, region)
    distinctmutations = combine(groupby(@subset(mutations, :Contained), names(mutations, Not([:Patient, :Mutation, :Inheritance]))), nrow=>:count)

    muts = buildmutation.(distinctmutations.RelReg, distinctmutations.type, distinctmutations.Ref, distinctmutations.Alt, distinctmutations.Contained);

    refseq = region.metadata
    altseqindex = buildaltseq.(Ref(refseq), muts)
    distinctmutations.LeftInd = first.(altseqindex)
    distinctmutations.AltInd = getindex.(altseqindex, 2)
    distinctmutations.RefRightInd = getindex.(altseqindex, 3)
    distinctmutations.AltRightInd = getindex.(altseqindex, 4)

    altseqs = last.(altseqindex)

    distinctmutations, refseq, altseqs
end


text:
@testset "curl" begin

N = (3,4,5)
M = prod(N)
r = reshape(collect(1:3M), M, 3)'[:]  # index mapping from block matrix to narrowly banded matrix
Z = spzeros(M,M)

F = rand(ComplexF64, N..., 3)
G = similar(F)
g = zeros(ComplexF64, 3M)

@testset "create_curl and apply_curl!" begin
    for ci = CartesianIndices((false:true,false:true,false:true))
        # Construct Curl for a uniform grid and Bloch boundaries.
        isfwd = Vector{Bool}([ci.I...])
        ∂x = (nw = 1; create_∂(nw, isfwd[nw], [N...]))
        ∂y = (nw = 2; create_∂(nw, isfwd[nw], [N...]))
        ∂z = (nw = 3; create_∂(nw, isfwd[nw], [N...]))
        Curl_blk = [Z -∂z ∂y;
                    ∂z Z -∂x;
                    -∂y ∂x Z]

        Curl = create_curl(isfwd, [N...], order_cmpfirst=false)

        # Test the overall coefficients.
        @test size(Curl) == (3M,3M)
        @test all(any(Curl.≠0, dims=1))  # no zero columns
        @test all(any(Curl.≠0, dims=2))  # no zero rows
        @test all(sum(Curl, dims=2) .== 0)  # all row sums are zero, because Curl * ones(M) = 0

        @test Curl == Curl_blk

        # Construct Curl for a nonuniform grid and general boundaries.
        ∆l⁻¹ = rand.(N)  # isfwd = true (false) uses ∆l⁻¹ at dual (primal) locations
        isbloch = [true, false, false]
        e⁻ⁱᵏᴸ = rand(ComplexF64, 3)
        ∂x = (nw = 1; create_∂(nw, isfwd[nw], [N...], ∆l⁻¹[nw], isbloch[nw], e⁻ⁱᵏᴸ[nw]))
        ∂y = (nw = 2; create_∂(nw, isfwd[nw], [N...], ∆l⁻¹[nw], isbloch[nw], e⁻ⁱᵏᴸ[nw]))
        ∂z = (nw = 3; create_∂(nw, isfwd[nw], [N...], ∆l⁻¹[nw], isbloch[nw], e⁻ⁱᵏᴸ[nw]))
        Curl_blk = [Z -∂z ∂y;
                    ∂z Z -∂x;
                    -∂y ∂x Z]

        Curl = create_curl(isfwd, ∆l⁻¹, isbloch, e⁻ⁱᵏᴸ, order_cmpfirst=false)

        # Test Curl.
        @test Curl == Curl_blk

        # Test Cartesian-component-major ordering.
        Curl_cmpfirst = create_curl(isfwd, ∆l⁻¹, isbloch, e⁻ⁱᵏᴸ, order_cmpfirst=true)
        @test Curl_cmpfirst == Curl[r,r]

        # Test apply_curl!.
        f = F[:]
        mul!(g, Curl, f)
        apply_curl!(G, F, Val(:(=)), isfwd, ∆l⁻¹, isbloch, e⁻ⁱᵏᴸ)
        @test G[:] ≈ g

        # print("matrix: "); @btime mul!($g, $Curl, $f)
        # print("matrix-free: "); @btime apply_curl!($G, $F, Val(:(=)), $isfwd, $∆l⁻¹, $isbloch, $e⁻ⁱᵏᴸ)

        # Test 1×1 curls for 1D, and 1×2, 2×1 curls for 2D.
        nw_blk = [0 -3 2;
                  3 0 -1;
                  -2 1 0]

        cmp_shps = ([1], [2], [3], [1,2], [2,3], [3,1])
        for cmp_shp = cmp_shps
            if length(cmp_shp) == 1
                cmp_outs = setdiff(([1],[2],[3]), tuple(cmp_shp))
            else  # length(cmp_shp) == 2
                cmp_outs = (cmp_shp, setdiff([1,2,3], cmp_shp))
            end

            for cmp_out = cmp_outs
                if length(cmp_shp) == 1
                    cmp_in = [6 - only(cmp_shp) - only(cmp_out)]
                else  # length(cmp_shp) == 2
                    cmp_in = setdiff([1,2,3], cmp_out)
                end

                Ncmp = N[cmp_shp]
                Mcmp = prod(Ncmp)
                Zcmp = spzeros(Mcmp,Mcmp)

                isfwdcmp = isfwd[cmp_shp]
                ∆l⁻¹cmp = ∆l⁻¹[cmp_shp]
                isblochcmp = isbloch[cmp_shp]
                e⁻ⁱᵏᴸcmp = e⁻ⁱᵏᴸ[cmp_shp]

                Curl = create_curl(isfwdcmp, ∆l⁻¹cmp, isblochcmp, e⁻ⁱᵏᴸcmp;
                                   cmp_shp, cmp_out, cmp_in, order_cmpfirst=false)

                ∂_blks = [[(nw = nw_blk[nv,nu];
                            if iszero(nw)
                                ∂_blk = Z
                            else
                                sn = sign(nw)
                                nw = abs(nw)
                                ind_nw = findfirst(cmp_shp.==nw)
                                ∂_blk = sn * create_∂(ind_nw, isfwd[nw], [Ncmp...], ∆l⁻¹[nw], isbloch[nw], e⁻ⁱᵏᴸ[nw])
                            end;
                            ∂_blk) for nu = cmp_in] for nv = cmp_out]

                Curl_blk = cat(((row_blks->cat(row_blks..., dims=2)).(∂_blks))..., dims=1)

                @test Curl == Curl_blk

                Fcmp = rand(ComplexF64, Ncmp..., length(cmp_in))
                Gcmp = similar(Fcmp, Ncmp..., length(cmp_out))
                gcmp = zeros(ComplexF64, length(Gcmp))

                fcmp = Fcmp[:]
                mul!(gcmp, Curl, fcmp)
                apply_curl!(Gcmp, Fcmp, Val(:(=)), isfwdcmp, ∆l⁻¹cmp, isblochcmp, e⁻ⁱᵏᴸcmp;
                            cmp_shp, cmp_out, cmp_in)
                @test Gcmp[:] ≈ gcmp
            end
        end
    end
end  # @testset "create_curl and apply_curl!"

@testset "curl of curl" begin
    # Construct Cu and Cv for a uniform grid and Bloch boundaries.
    ∆ldual⁻¹ = ones.(N)  # inverse of ∆l evaluated at dual locations
    ∆lprim⁻¹ = ones.(N)  # inverse of ∆l evaluated at primal locations
    isbloch = [true, false, false]
    e⁻ⁱᵏᴸ = ones(3)

    Cu = create_curl([true,true,true], ∆ldual⁻¹, isbloch, e⁻ⁱᵏᴸ, order_cmpfirst=false)
    Cv = create_curl([false,false,false], ∆lprim⁻¹, isbloch, e⁻ⁱᵏᴸ, order_cmpfirst=false)

    # Test symmetry of each block.
    for i = 1:3
        for j = mod1.(i .+ [1,2], 3)
            -Cv[(i-1)*M+1:i*M,(j-1)*M+1:j*M]' == Cu[(i-1)*M+1:i*M,(j-1)*M+1:j*M]
        end
    end

    # Construct Cv * Cu.
    isbloch = fill(true, 3)
    Cu = create_curl([true,true,true], ∆ldual⁻¹, isbloch, e⁻ⁱᵏᴸ, order_cmpfirst=false)
    Cv = create_curl([false,false,false], ∆lprim⁻¹, isbloch, e⁻ⁱᵏᴸ, order_cmpfirst=false)
    A = Cv * Cu

    # Test curl of curl.
    @test all(diag(A) .== 4)  # all diagonal entries are 4
    @test all(sum(A.≠0, dims=2) .== 13)  # 13 nonzero entries per row
    @test A == A'  # Hermitian

    B = A - 4I
    @test all(abs.(B[B.≠0]).==1)  # all nonzero off-diagonal entries are ±1
end  # @testset "curl of curl"

@testset "curl of curl, mixed forward and backward" begin
    # Construct Cu and Cv for a uniform grid and Bloch boundaries.
    ∆ldual⁻¹ = ones.(N)
    ∆lprim⁻¹ = ones.(N)
    isbloch = [true, false, false]
    e⁻ⁱᵏᴸ = ones(3)

    isfwd = [true,false,true]
    Cu = create_curl(isfwd, ∆ldual⁻¹, isbloch, e⁻ⁱᵏᴸ, order_cmpfirst=false)
    Cv = create_curl(.!isfwd, ∆lprim⁻¹, isbloch, e⁻ⁱᵏᴸ, order_cmpfirst=false)

    # Test symmetry of each block.
    for i = 1:3
        for j = mod1.(i .+ [1,2], 3)
            -Cv[(i-1)*M+1:i*M,(j-1)*M+1:j*M]' == Cu[(i-1)*M+1:i*M,(j-1)*M+1:j*M]
        end
    end

    # Construct Cv * Cu.
    isbloch = fill(true, 3)
    Cu = create_curl(isfwd, ∆ldual⁻¹, isbloch, e⁻ⁱᵏᴸ, order_cmpfirst=false)
    Cv = create_curl(.!isfwd, ∆lprim⁻¹, isbloch, e⁻ⁱᵏᴸ, order_cmpfirst=false)
    A = Cv * Cu

    # Test curl of curl.
    @test all(diag(A) .== 4)  # all diagonal entries are 4
    @test all(sum(A.≠0, dims=2) .== 13)  # 13 nonzero entries per row
    @test A == A'  # Hermitian

    B = A - 4I
    @test all(abs.(B[B.≠0]).==1)  # all nonzero off-diagonal entries are ±1
end  # @testset "curl of curl, mixed forward and backward"

end  # @testset "curl"


text:
using SignalDecomposition
using Interpolations
using Plots

const _SDC = SignalDecomposition
const _INT = Interpolations

const freq = 50.0
const nw_id_default = 1

xfmr = Dict("voltage_harmonics" => [1,3],
            "current_harmonics" => [1,3], #this one has to match the data in the matpower file
            "N" => 50,
            "current_type" => :rectangular,
            "excitation_type" => :sigmoid,
            "inom" => 0.7,
            "ψmax" => 0.005,
            "voltage_type" => :rectangular,
            "dv" => [0.1,0.1],
            "vmin" => -[1.1,1.1],
            "vmax" => [1.1,1.1],
            "dθ" => [π/5,π/5],
            "θmin" => [0.0,0.0],
            "θmax" => [2π,2π])

excitation_flux_polar(V, θ, w, t) =
    sum(V[n] ./ w[n] .* sin.(w[n] .* t .+ θ[n]) for n in 1:length(V))
excitation_flux_rectangular(Vre, Vim, w, t) =
    excitation_flux_polar(hypot.(Vre,Vim), atan.(Vim,Vre), w, t)
excitation_current_sigmoid(inom, ψmax, ψ) =
    -inom .* log.(2.0 ./ (ψ ./ ψmax .+ 1.0) .- 1.0)

""
function sample_voltage_polar(voltage_harmonics, dv, vmin, vmax, dθ, θmin, θmax)
    S = reduce(vcat,[[vmin[ni]:dv[ni]:vmax[ni],θmin[ni]:dθ[ni]:θmax[ni]]
                      for (ni,nh) in enumerate(voltage_harmonics)])
    R = [1:length(s) for s in S]
    return S, R
end
function sample_voltage_rectangular(voltage_harmonics, dv, vmin, vmax)
    S = reduce(vcat,[[vmin[ni]:dv[ni]:vmax[ni],vmin[ni]:dv[ni]:vmax[ni]]
                      for (ni,nh) in enumerate(voltage_harmonics)])
    R = [1:length(s) for s in S]
    return S, R
end

xfmr_exc = xfmr

# HARMONICS
voltage_harmonics = xfmr_exc["voltage_harmonics"]
current_harmonics = xfmr_exc["current_harmonics"]

# TIME
N = xfmr_exc["N"]

# EXCITATION CURRENT
current_type = xfmr_exc["current_type"]

# EXCITATION FLUX
excitation_type = xfmr_exc["excitation_type"]
inom = xfmr_exc["inom"]
ψmax = xfmr_exc["ψmax"]

# EXCITATION VOLTAGE
voltage_type = xfmr_exc["voltage_type"]
dv   = xfmr_exc["dv"]
vmin = xfmr_exc["vmin"]
vmax = xfmr_exc["vmax"]
dθ   = xfmr_exc["dθ"]
θmin = xfmr_exc["θmin"]
θmax = xfmr_exc["θmax"]

# TIME
w = (2.0 * pi * freq) .* voltage_harmonics
t = 0.0:(1 / (freq * N * maximum(current_harmonics))):(5.0 / freq)

# DECOMPOSITION
fq = _SDC.Sinusoidal(freq .* current_harmonics)

# SAMPLE VOLTAGE
if voltage_type == :polar
    S, R = sample_voltage_polar(voltage_harmonics, dv, vmin, vmax, dθ, θmin, θmax)
elseif voltage_type == :rectangular
    S, R = sample_voltage_rectangular(voltage_harmonics, dv, vmin, vmax)
end
Ia = Dict(nh => zeros(R...) for nh in current_harmonics)
Ib = Dict(nh => zeros(R...) for nh in current_harmonics)

nr = [20,20,12,12]
sample = [S[ni][ns] for (ni,ns) in enumerate(nr)]

# if voltage_type == :polar
#     V, θ = sample[1:2:end], sample[2:2:end]
#     ψexc = excitation_flux_polar(V, θ, w, t)
# else voltage_type == :rectangular
Vre, Vim = sample[1:2:end], sample[2:2:end]
ψexc = excitation_flux_rectangular(Vre, Vim, w, t)
# end

# if excitation_type == :sigmoid
I_exc = excitation_current_sigmoid(inom, ψmax, ψexc)
# end

plot(t,ψexc)
plot(t,I_exc)

_SDC.decompose(t, I_exc, fq)

I, φ = fq.A[2:end], -fq.φ[2:end]


text:
function minibatch(x, y, batchsize; atype=Array{Float32}, xrows=784, yrows=10, xscale=255)
    xbatch(a)=convert(atype, reshape(a./xscale, xrows, div(length(a),xrows)))
    ybatch(a)=(a[a.==0]=10; convert(atype, sparse(convert(Vector{Int},a),1:length(a),one(eltype(a)),yrows,length(a))))
    xcols = div(length(x),xrows)
    xcols == length(y) || throw(DimensionMismatch())
    data = Any[]
    for i=1:batchsize:xcols-batchsize+1
        j=i+batchsize-1
        push!(data, (xbatch(x[1+(i-1)*xrows:j*xrows]), ybatch(y[i:j])))
    end
    return data
end

function loaddata()
    info("Loading MNIST...")
    gzload("train-images-idx3-ubyte.gz")[17:end],
    gzload("t10k-images-idx3-ubyte.gz")[17:end],
    gzload("train-labels-idx1-ubyte.gz")[9:end],
    gzload("t10k-labels-idx1-ubyte.gz")[9:end]
end

function gzload(file; path=Knet.dir("data",file), url="http://yann.lecun.com/exdb/mnist/$file")
    isfile(path) || download(url, path)
    f = gzopen(path)
    a = read(f)
    close(f)
    return(a)
end

max_stars_repo_name: ilkerkesen/adversarial

max_issues_repo_name: muhammetfaik/GAN

max_forks_repo_name: muhammetfaik/GAN


text:
function parse_line(line)
    id, x, y, w, h = parse.(Int, match(r"#(\d+) @ (\d+),(\d+): (\d+)x(\d+)", line).captures)
    (id, x .+ (1:w), y .+ (1:h))
end

input = parse_line.(readlines("03.input"))

cloth = zeros(Int, 1000, 1000)
for (id, xrange, yrange) in input
    cloth[xrange, yrange] .+= 1
end

# part 1: number of squares with >1 claim
@show sum(c > 1 for c in cloth)

# part 2: find only claim that has no intersection
# (re-use the "cloth" marked up for the first star to make this O(n))

@btime [id
        for (id, xrange, yrange)
        in input
        if all(isequal(1), view(cloth, xrange, yrange))]
#   452.265 μs (5281 allocations: 185.72 KiB)
# 1-element Array{Int64,1}:
#  724

@btime first(id
             for (id, xrange, yrange)
             in input
             if all(isequal(1), view(cloth, xrange, yrange)))
#   244.011 μs (2899 allocations: 101.86 KiB)
# 724


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".

function _co_aberration(jd::T, ra::T, dec::T, eps::T) where {T<:AbstractFloat}
    t = (jd - J2000) / JULIANCENTURY
    if isnan(eps)
        eps = true_obliquity(jd)
    end
    sunlong = sunpos(jd, radians=true)[3]
    e = @evalpoly t 0.016708634 -0.000042037 -0.0000001267
    pe = @evalpoly t 102.93735 1.71946 0.00046
    sd, cd = sincos(deg2rad(dec))
    ce = cos(eps)
    te = tan(eps)
    sp, cp = sincos(deg2rad(pe))
    ss, cs = sincos(sunlong)
    sa, ca = sincos(deg2rad(ra))
    t1 = (cs*ce*(te*cd - sa*sd) + ca*sd*ss)
    t2 = (cp*ce*(te*cd - sa*sd) + ca*sd*sp)
    d_ra = 20.49552*(e*(ca*cp*ce + sa*sp) - ca*cs*ce - sa*ss)/cd
    d_dec = 20.49552*(e*t2 - t1)
    return d_ra, d_dec
end

"""
    co_aberration(jd, ra, dec[, eps=NaN]) -> d_ra, d_dec

### Purpose ###

Calculate changes to right ascension and declination due to the effect
of annual aberration

### Explanation ###

With reference to Meeus, Chapter 23

### Arguments ###

* `jd`: julian date, scalar or vector
* `ra`: right ascension in degrees, scalar or vector
* `dec`: declination in degrees, scalar or vector
* `eps` (optional): true obliquity of the ecliptic (in radians). It will be
  calculated if no argument is specified.

### Output ###

The 2-tuple `(d_ra, d_dec)`:

* `d_ra`: correction to right ascension due to aberration, in arc seconds
* `d_dec`: correction to declination due to aberration, in arc seconds

### Example ###

Compute the change in RA and Dec of Theta Persei (RA = 2h46m,11.331s, Dec = 49d20',54.5'')
due to aberration on 2028 Nov 13.19 TD

```jldoctest
julia> using AstroLib

julia> jd = jdcnv(2028,11,13,4, 56)
2.4620887055555554e6

julia> co_aberration(jd,ten(2,46,11.331)*15,ten(49,20,54.54))
(30.04404628365077, 6.699400463119431)
```

d_ra = 30.04404628365103'' (≈ 2.003s)
d_dec = 6.699400463118504''

### Notes ###

Code of this function is based on IDL Astronomy User's Library.

The output d_ra is *not* multiplied by cos(dec), so that
apparent_ra = ra + d_ra/3600.

These formula are from Meeus, Chapters 23.  Accuracy is much better than 1
arcsecond. The maximum deviation due to annual aberration is 20.49'' and occurs when the
Earth's velocity is perpendicular to the direction of the star.

This function calls [`true_obliquity`](@ref) and [`sunpos`](@ref).
"""
co_aberration(jd::Real, ra::Real, dec::Real, eps::Real=NaN) =
    _co_aberration(promote(float(jd), float(ra), float(dec), float(eps))...)

function co_aberration(jd::AbstractVector{R}, ra::AbstractVector{R},
                       dec::AbstractVector{R}, eps::Real=NaN) where {R<:Real}
    @assert length(jd) == length(ra) == length(dec) "jd, ra and dec vectors should be of the same length"
    typejd = float(R)
    ra_out  = similar(ra,  typejd)
    dec_out = similar(dec, typejd)
    for i in eachindex(jd)
        ra_out[i], dec_out[i] = co_aberration(jd[i], ra[i], dec[i], eps)
    end
    return ra_out, dec_out
end


text:
# Used "Conjugate Bayesian analysis of the Gaussian distribution" by Murphy as
# a reference.  Note that there were some typos in that document so the code
# here may not correspond exactly.

struct NormalInverseWishart{T<:Real,V<:AbstractVector{T},M<:AbstractMatrix{T}} <: ContinuousUnivariateDistribution
    dim::Int
    zeromean::Bool
    mu::V
    kappa::T              # This scales precision (inverse covariance)
    Lamchol::Cholesky{T,M}  # Covariance matrix (well, sqrt of one)
    nu::T

    function NormalInverseWishart{T}(mu::AbstractVector{T}, kappa::T,
                                     Lamchol::Cholesky{T,M}, nu::T) where {T<:Real, M<:AbstractMatrix{T}}
        # Probably should put some error checking in here
        d = length(mu)
        zmean::Bool = true
        for i = 1:d
            if !iszero(mu[i])
                zmean = false
                break
            end
        end
        new{T,typeof(mu),M}(d, zmean, mu, T(kappa), Lamchol, T(nu))
    end
end

function NormalInverseWishart(mu::AbstractVector{U}, kappa::Real,
                                Lamchol::Cholesky{S}, nu::Real) where {S<:Real, U<:Real}
    T = promote_type(eltype(mu), typeof(kappa), typeof(nu), S)
    return NormalInverseWishart{T}(convert(AbstractVector{T}, mu),
                                   T(kappa),
                                   Cholesky{T}(Lamchol),
                                   T(nu))
end

function NormalInverseWishart(mu::AbstractVector{U}, kappa::Real,
                              Lambda::AbstractMatrix{S}, nu::Real) where {S<:Real, U<:Real}
    T = promote_type(eltype(mu), typeof(kappa), typeof(nu), S)
    return NormalInverseWishart{T}(convert(AbstractVector{T}, mu),
                                   T(kappa),
                                   Cholesky{T}(cholesky(Lambda)),
                                   T(nu))
end

function insupport(::Type{NormalInverseWishart}, x::AbstractVector{T}, Sig::AbstractMatrix{T}) where T<:Real
    return (all(isfinite, x) &&
           size(Sig, 1) == size(Sig, 2) &&
           isApproxSymmmetric(Sig) &&
           size(Sig, 1) == length(x) &&
           hasCholesky(Sig))
end

"""
    params(niw::NormalInverseWishart)

The parameters are
* μ::AbstractVector{T<:Real} the expected mean vector
* Λchol::Cholesky{T<:Real} the Cholesky decomposition of the scale matrix
* κ::T<:Real prior pseudocount for the mean
* ν::T<:Real prior pseudocount for the covariance
"""
params(niw::NormalInverseWishart) = (niw.mu, niw.Lamchol, niw.kappa, niw.nu)

pdf(niw::NormalInverseWishart, x::AbstractVector{T}, Sig::AbstractMatrix{T}) where T<:Real =
        exp(logpdf(niw, x, Sig))

function logpdf(niw::NormalInverseWishart, x::AbstractVector{T}, Sig::AbstractMatrix{T}) where T<:Real
    if !insupport(NormalInverseWishart, x, Sig)
        return -Inf
    else
        p = size(x, 1)

        nu = niw.nu
        kappa = niw.kappa
        mu = niw.mu
        Lamchol = niw.Lamchol
        hnu = 0.5 * nu
        hp = 0.5 * p

        # Normalization
        logp::T = hnu * logdet(Lamchol)
        logp -= hnu * p * log(2.)
        logp -= logmvgamma(p, hnu)
        logp -= hp * (log(2.0*pi) - log(kappa))

        # Inverse-Wishart
        logp -= (hnu + hp + 1.) * logdet(Sig)
        logp -= 0.5 * tr(Sig \ convert(AbstractMatrix, Lamchol))

        # Normal
        z = niw.zeromean ? x : x - mu
        logp -= 0.5 * kappa * invquad(PDMat(Sig), z)

        return logp

    end
end

function rand(niw::NormalInverseWishart)
    Sig = rand(InverseWishart(niw.nu, niw.Lamchol))
    mu = rand(MvNormal(niw.mu, Sig ./ niw.kappa))
    return (mu, Sig)
end


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

@doc raw"""PodSecurityContext holds pod-level security attributes and common container settings. Some fields are also present in container.securityContext.  Field values of container.securityContext take precedence over field values of PodSecurityContext.

    IoK8sApiCoreV1PodSecurityContext(;
        fsGroup=nothing,
        runAsGroup=nothing,
        runAsNonRoot=nothing,
        runAsUser=nothing,
        seLinuxOptions=nothing,
        supplementalGroups=nothing,
        sysctls=nothing,
        windowsOptions=nothing,
    )

    - fsGroup::Int64 : A special supplemental group that applies to all containers in a pod. Some volume types allow the Kubelet to change the ownership of that volume to be owned by the pod:  1. The owning GID will be the FSGroup 2. The setgid bit is set (new files created in the volume will be owned by FSGroup) 3. The permission bits are OR'd with rw-rw----  If unset, the Kubelet will not modify the ownership and permissions of any volume.
    - runAsGroup::Int64 : The GID to run the entrypoint of the container process. Uses runtime default if unset. May also be set in SecurityContext.  If set in both SecurityContext and PodSecurityContext, the value specified in SecurityContext takes precedence for that container.
    - runAsNonRoot::Bool : Indicates that the container must run as a non-root user. If true, the Kubelet will validate the image at runtime to ensure that it does not run as UID 0 (root) and fail to start the container if it does. If unset or false, no such validation will be performed. May also be set in SecurityContext.  If set in both SecurityContext and PodSecurityContext, the value specified in SecurityContext takes precedence.
    - runAsUser::Int64 : The UID to run the entrypoint of the container process. Defaults to user specified in image metadata if unspecified. May also be set in SecurityContext.  If set in both SecurityContext and PodSecurityContext, the value specified in SecurityContext takes precedence for that container.
    - seLinuxOptions::IoK8sApiCoreV1SELinuxOptions : The SELinux context to be applied to all containers. If unspecified, the container runtime will allocate a random SELinux context for each container.  May also be set in SecurityContext.  If set in both SecurityContext and PodSecurityContext, the value specified in SecurityContext takes precedence for that container.
    - supplementalGroups::Vector{Int64} : A list of groups applied to the first process run in each container, in addition to the container's primary GID.  If unspecified, no groups will be added to any container.
    - sysctls::Vector{IoK8sApiCoreV1Sysctl} : Sysctls hold a list of namespaced sysctls used for the pod. Pods with unsupported sysctls (by the container runtime) might fail to launch.
    - windowsOptions::IoK8sApiCoreV1WindowsSecurityContextOptions : The Windows specific settings applied to all containers. If unspecified, the options within a container's SecurityContext will be used. If set in both SecurityContext and PodSecurityContext, the value specified in SecurityContext takes precedence.
"""
mutable struct IoK8sApiCoreV1PodSecurityContext <: SwaggerModel
    fsGroup::Any # spec type: Union{ Nothing, Int64 } # spec name: fsGroup
    runAsGroup::Any # spec type: Union{ Nothing, Int64 } # spec name: runAsGroup
    runAsNonRoot::Any # spec type: Union{ Nothing, Bool } # spec name: runAsNonRoot
    runAsUser::Any # spec type: Union{ Nothing, Int64 } # spec name: runAsUser
    seLinuxOptions::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SELinuxOptions } # spec name: seLinuxOptions
    supplementalGroups::Any # spec type: Union{ Nothing, Vector{Int64} } # spec name: supplementalGroups
    sysctls::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1Sysctl} } # spec name: sysctls
    windowsOptions::Any # spec type: Union{ Nothing, IoK8sApiCoreV1WindowsSecurityContextOptions } # spec name: windowsOptions

    function IoK8sApiCoreV1PodSecurityContext(;fsGroup=nothing, runAsGroup=nothing, runAsNonRoot=nothing, runAsUser=nothing, seLinuxOptions=nothing, supplementalGroups=nothing, sysctls=nothing, windowsOptions=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("fsGroup"), fsGroup)
        setfield!(o, Symbol("fsGroup"), fsGroup)
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("runAsGroup"), runAsGroup)
        setfield!(o, Symbol("runAsGroup"), runAsGroup)
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("runAsNonRoot"), runAsNonRoot)
        setfield!(o, Symbol("runAsNonRoot"), runAsNonRoot)
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("runAsUser"), runAsUser)
        setfield!(o, Symbol("runAsUser"), runAsUser)
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("seLinuxOptions"), seLinuxOptions)
        setfield!(o, Symbol("seLinuxOptions"), seLinuxOptions)
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("supplementalGroups"), supplementalGroups)
        setfield!(o, Symbol("supplementalGroups"), supplementalGroups)
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("sysctls"), sysctls)
        setfield!(o, Symbol("sysctls"), sysctls)
        validate_property(IoK8sApiCoreV1PodSecurityContext, Symbol("windowsOptions"), windowsOptions)
        setfield!(o, Symbol("windowsOptions"), windowsOptions)
        o
    end
end # type IoK8sApiCoreV1PodSecurityContext

const _property_map_IoK8sApiCoreV1PodSecurityContext = Dict{Symbol,Symbol}(Symbol("fsGroup")=>Symbol("fsGroup"), Symbol("runAsGroup")=>Symbol("runAsGroup"), Symbol("runAsNonRoot")=>Symbol("runAsNonRoot"), Symbol("runAsUser")=>Symbol("runAsUser"), Symbol("seLinuxOptions")=>Symbol("seLinuxOptions"), Symbol("supplementalGroups")=>Symbol("supplementalGroups"), Symbol("sysctls")=>Symbol("sysctls"), Symbol("windowsOptions")=>Symbol("windowsOptions"))
const _property_types_IoK8sApiCoreV1PodSecurityContext = Dict{Symbol,String}(Symbol("fsGroup")=>"Int64", Symbol("runAsGroup")=>"Int64", Symbol("runAsNonRoot")=>"Bool", Symbol("runAsUser")=>"Int64", Symbol("seLinuxOptions")=>"IoK8sApiCoreV1SELinuxOptions", Symbol("supplementalGroups")=>"Vector{Int64}", Symbol("sysctls")=>"Vector{IoK8sApiCoreV1Sysctl}", Symbol("windowsOptions")=>"IoK8sApiCoreV1WindowsSecurityContextOptions")
Base.propertynames(::Type{ IoK8sApiCoreV1PodSecurityContext }) = collect(keys(_property_map_IoK8sApiCoreV1PodSecurityContext))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodSecurityContext }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PodSecurityContext[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodSecurityContext }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodSecurityContext[property_name]

function check_required(o::IoK8sApiCoreV1PodSecurityContext)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PodSecurityContext }, name::Symbol, val)
end


text:
using FastDMTransform
using Documenter

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

makedocs(;
    modules=[FastDMTransform],
    authors="Kiran Shila <me@kiranshila.com> and contributors",
    repo="https://github.com/kiranshila/FastDMTransform.jl/blob/{commit}{path}#{line}",
    sitename="FastDMTransform.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://kiranshila.github.io/FastDMTransform.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "API" => "api.md"
    ],
)

deploydocs(;
    repo="github.com/kiranshila/FastDMTransform.jl",
    devbranch="main",
)


text:
function SimpsonIntegrator()
    return SimpsonIntegrator(())
end

function SimpsonIntegrator(arg0::jdouble, arg1::jdouble, arg2::jint, arg3::jint)
    return SimpsonIntegrator((jdouble, jdouble, jint, jint), arg0, arg1, arg2, arg3)
end

function SimpsonIntegrator(arg0::jint, arg1::jint)
    return SimpsonIntegrator((jint, jint), arg0, arg1)
end


text:
struct STE <: AbstractLoss
    σ::Real
    constant::Float64

    function STE(;σ::T = 1/sqrt(2)) where T <: Real
        σ > 0 || throw(ArgumentError("σ in STE loss must be > 0"))
        new(σ, 1/σ^2)
    end

end

@doc raw"""
    function kernel(loss::STE, X::AbstractMatrix)

Computes:

```math
K = \exp(\|X_i - X_j\|^2/(2\sigma^2)) \forall (i,j) \in {1,\ldots,n} \times {1,\ldots,n}.
```
"""
function kernel(loss::STE, X::Embedding)
    K = pairwise(SqEuclidean(), X.X, dims=2)
    c = -loss.constant / 2

    for ij in eachindex(K)
        @inbounds K[ij] = exp(c * K[ij])
    end
    return K
end

function gradient(loss::STE, triplets, X::AbstractMatrix)

    K = kernel(loss, X) # Triplet kernel values (in the STE loss)

    # We need to create an array to prevent race conditions
    # This is the best average solution for small and big Embeddings
    nthreads = Threads.nthreads()
    triplets_range = partition(length(triplets), nthreads)

    C = zeros(Float64, nthreads)
    ∇C = [zeros(Float64, size(X)) for _ in 1:nthreads]

    Threads.@threads for tid in 1:nthreads
        C[tid] = tgradient!(∇C[tid], loss, triplets, X, K, triplets_range[tid])
    end

    return sum(C), -sum(∇C)
end

function tgradient!(
    ∇C::Matrix{<:AbstractFloat},
    loss::STE,
    triplets::Triplets,
    X::AbstractMatrix,
    K::Matrix{<:AbstractFloat},
    triplets_range::UnitRange{Int64})

    C = 0.0

    for t in triplets_range
        @views @inbounds i, j, k = triplets[t][:i], triplets[t][:j], triplets[t][:k]

        @inbounds P = K[i,j] / (K[i,j] + K[i,k])
        C += -log(P)

        for d in 1:ndims(X)
            @inbounds ∂x_j = (1 - P) * (X[d,i] - X[d,j])
            @inbounds ∂x_k = (1 - P) * (X[d,i] - X[d,k])

            @inbounds ∇C[d,i] += - loss.constant * (∂x_j - ∂x_k)
            @inbounds ∇C[d,j] +=   loss.constant *  ∂x_j
            @inbounds ∇C[d,k] += - loss.constant *  ∂x_k
        end
    end
    return C
end

function tgradient!(
    ∇C::Matrix{<:AbstractFloat},
    loss::STE,
    triplets::LabeledTriplets,
    X::AbstractMatrix,
    K::Matrix{<:AbstractFloat},
    triplets_range::UnitRange{Int64})

    C = 0.0

    for t in triplets_range
        @views @inbounds i, j, k = if triplets[t][:y] == -1
            triplets[t][:i], triplets[t][:j], triplets[t][:k]
        else
            triplets[t][:i], triplets[t][:k], triplets[t][:j]
        end

        @inbounds P = K[i,j] / (K[i,j] + K[i,k])
        C += -log(P)

        for d in 1:ndims(X)
            @inbounds ∂x_j = (1 - P) * (X[d,i] - X[d,j])
            @inbounds ∂x_k = (1 - P) * (X[d,i] - X[d,k])

            @inbounds ∇C[d,i] += - loss.constant * (∂x_j - ∂x_k)
            @inbounds ∇C[d,j] +=   loss.constant *  ∂x_j
            @inbounds ∇C[d,k] += - loss.constant *  ∂x_k
        end
    end
    return C
end


text:
using Cobweb
using Cobweb: h, Page, CSS
using Markdown

css = Dict(
    "html" => Dict(
        "font-family" => "Arial"
    )
)

page = h.html(
    h.head(
        h.meta(charset="UTF-8"),
        h.meta(name="viewport", content="width=device-width, initial-scale=1.0"),
        h.title("Cobweb.jl Docs"),
        CSS(css),
    ),
    h.body(
        h.h1("This page was built with ", h.code("Cobweb.jl"), "."),
        Markdown.parse("""
        Take a look at [`docs/make.jl`](https://github.com/joshday/Cobweb.jl/blob/main/docs/make.jl) inside the [`Cobweb.jl` repo](https://github.com/joshday/Cobweb.jl).
        """),
        h.button("Click Me for an alert!", onclick="buttonClicked()"),
        Cobweb.Javascript("const buttonClicked = () => alert('This button was clicked!')"),
    )
)

index_html = touch(joinpath(mkpath(joinpath(@__DIR__, "build")), "index.html"))

Cobweb.save(Page(page), index_html)


text:
"""
    density_2D( rho::Real, pixelSideLength::Real,
                Mass::Real=1.989e43, Length::Real=3.085678e21)

Computes the density in units of [g/cm^2].
"""
function density_2D(rho::Real, pixelSideLength::Real,
                    Mass::Real=1.989e43, Length::Real=3.085678e21)
    return rho * ( (Mass / Length^2 ) * pixelSideLength / Length )
end


text:
using JULES, Oceananigans
using Printf
using Plots
using NLsolve

Nx = Ny = 1
Nz = 32
L = 10e3

grid = RegularCartesianGrid(size=(Nx, Ny, Nz), halo=(2, 2, 2),
                            x=(0, L), y=(0, L), z=(0, L))
Δz = grid.Δz

pₛ = 100000
Tₐ = 293.15
g  = 9.80665

buoyancy = IdealGas()
Rᵈ, cₚ, γ = buoyancy.Rᵈ, buoyancy.cₚ, buoyancy.γ

####
#### Isothermal atmosphere
####

H = Rᵈ * Tₐ / g    # Scale height [m]
ρₛ = pₛ / (Rᵈ*Tₐ)  # Surface density [kg/m³]

p₀(x, y, z) = pₛ * exp(-z/H)
ρ₀(x, y, z) = ρₛ * exp(-z/H)

θ₀(x, y, z) = Tₐ * exp(z/H * Rᵈ/cₚ)
Θ₀(x, y, z) = ρ₀(x, y, z) * θ₀(x, y, z)

θ₀_prof = θ₀.(0, 0, grid.zC)
p₀_prof = p₀.(0, 0, grid.zC)
ρ₀_prof = ρ₀.(0, 0, grid.zC)

####
#### Objective function describing the ideal gas equation of state and hydrostatic balance
#### Set up a nonlinear system of equations F(x) = 0 and solve for pₖ and ρₖ.
####

# First Nz variables are pressure values pₖ, k = 1, 2, …, Nz
# Second Nz variables are density values ρₖ
function f!(F, x)
    @inbounds begin
        p₁, ρ₁, θ₁ = x[1], x[Nz+1], θ₀_prof[1]
        F[1] = p₁ - pₛ * (Rᵈ*ρ₁*θ₁ / pₛ)^γ
        F[Nz+1] = p₁ - pₛ + ρ₁*g*(Δz/2)

        for k in 2:Nz
            pₖ, pₖ₋₁, ρₖ, θₖ = x[k], x[k-1], x[Nz+k], θ₀_prof[k]
            F[k] = pₖ - pₛ * (Rᵈ*ρₖ*θₖ / pₛ)^γ
            F[Nz+k] = pₖ - pₖ₋₁ + ρₖ*g*Δz
        end
    end
end

guess = [p₀_prof; ρ₀_prof]
sol = nlsolve(f!, guess, show_trace=true)

####
#### Plot difference between initial condition and solution.
####

p∞_prof = sol.zero[1:Nz]
ρ∞_prof = sol.zero[Nz+1:2Nz]

p_plot = plot(p₀_prof, grid.zC, xlabel="pressure (Pa)", ylabel="z (m)", label="initial")
plot!(p_plot, p∞_prof, grid.zC, label="balanced")

ρ_plot = plot(ρ₀_prof, grid.zC, xlabel="rho (kg/m³)", ylabel="z (m)", label="initial")
plot!(ρ_plot, ρ∞_prof, grid.zC, label="balanced")

display(plot(p_plot, ρ_plot, show=true))

####
#### Set up compressible model
####

model = CompressibleModel(grid=grid, buoyancy=buoyancy, reference_pressure=pₛ,
                          thermodynamic_variable=ModifiedPotentialTemperature(),
                          tracers=(:Θᵐ,))

Θ∞_prof = ρ∞_prof .* θ₀_prof

set!(model.density, reshape(ρ∞_prof, (Nx, Ny, Nz)))
set!(model.tracers.Θᵐ, reshape(Θ∞_prof, (Nx, Ny, Nz)))

times = [model.clock.time]
ρw_ts = [model.momenta.W[1, 1, Int(Nz/2)]]

while model.clock.time < 50
    time_step!(model; Δt=0.5, Nt=10)

    @show model.clock.time
    Θ_prof = model.tracers.Θᵐ[1, 1, 1:Nz] .- Θ∞_prof
    W_prof = model.momenta.W[1, 1, 1:Nz+1]
    ρ_prof = model.density[1, 1, 1:Nz] .- ρ∞_prof

    Θ_plot = plot(Θ_prof, grid.zC, xlim=(-2.5, 2.5), xlabel="Theta_prime", ylabel="z (m)", label="")
    W_plot = plot(W_prof, grid.zF, xlim=(-1.2, 1.2), xlabel="rho*w", label="")
    ρ_plot = plot(ρ_prof, grid.zC, xlim=(-0.01, 0.01), xlabel="rho", label="")

    t_str = @sprintf("t = %d s", model.clock.time)
    display(plot(Θ_plot, W_plot, ρ_plot, title=["" "$t_str" ""], layout=(1, 3), show=true))

    push!(times, model.clock.time)
    push!(ρw_ts, model.momenta.W[1, 1, Int(Nz/2)])
end


text:
import OnlinePackage
import Documenter: doctest

doctest(OnlinePackage)


text:
"Prints a benchmark results just like btime"
function print_benchmark(trialmin)
    trialtime = BenchmarkTools.time(trialmin)
    trialallocs = BenchmarkTools.allocs(trialmin)
    println(
        "  ",
        BenchmarkTools.prettytime(trialtime),
        " (",
        trialallocs,
        " allocation",
        trialallocs == 1 ? "" : "s",
        ": ",
        BenchmarkTools.prettymemory(BenchmarkTools.memory(trialmin)),
        ")",
    )
end

# Not all special functions have been implemented yet:
# https://github.com/JuliaMath/SpecialFunctions.jl/issues/19
# (MIT licensed)
#
# This expint implementation is copied from:
# https://github.com/mschauer/Bridge.jl/blob/a719e8a57aaa77fa4ca289c25978a7e080908826/src/expint.jl#L1
########################################################################
# Inlined, optimized code for the exponential integral E₁ in double precison
# by Steven G. Johnson (@stevengj), Code under MIT license.
# For more explanations, see course notes by Steven G. Johnson at
#     https://github.com/stevengj/18S096-iap17/blob/master/pset3/pset3-solutions.ipynb

# n coefficients of the Taylor series of E₁(z) + log(z), in type T:
function E₁_taylor_coefficients(::Type{T}, n::Integer) where T <: Number
    n < 0 && throw(ArgumentError("$n ≥ 0 is required"))
    n == 0 && return T[]
    n == 1 && return T[-eulergamma]
    # iteratively compute the terms in the series, starting with k=1
    term::T = 1
    terms = T[-eulergamma, term]
    for k in 2:n
        term = -term * (k - 1) / (k * k)
        push!(terms, term)
    end
    return terms
end

# inline the Taylor expansion for a given order n, in double precision
macro E₁_taylor64(z, n::Integer)
    c = E₁_taylor_coefficients(Float64, n)
    taylor = :(@evalpoly zz)
    append!(taylor.args, c)
    quote
        let zz = $(esc(z))
            $taylor - log(zz)
        end
    end
end

# for numeric-literal coefficients: simplify to a ratio of two polynomials:
# return (p,q): the polynomials p(x) / q(x) corresponding to E₁_cf(x, a...),
# but without the exp(-x) term
function E₁_cfpoly(n::Integer, ::Type{T}=BigInt) where T <: Real
    q = Polynomials.Polynomial(T[1])
    p = x = Polynomials.Polynomial(T[0,1])
    for i in n:-1:1
        p, q = x * p + (1 + i) * q, p # from cf = x + (1+i)/cf = x + (1+i)*q/p
        p, q = p + i * q, p     # from cf = 1 + i/cf = 1 + i*q/p
    end
    # do final 1/(x + inv(cf)) = 1/(x + q/p) = p/(x*p + q)
    return p, x * p + q
end

macro E₁_cf64(z, n::Integer)
    p, q = E₁_cfpoly(n, BigInt)
    num_expr =  :(@evalpoly zz)
    append!(num_expr.args, Float64.(Polynomials.coeffs(p)))
    den_expr = :(@evalpoly zz)
    append!(den_expr.args, Float64.(Polynomials.coeffs(q)))
    quote
        let zz = $(esc(z))
            exp(-zz) * $num_expr / $den_expr
        end
    end
end

# exponential integral function E₁(z)
function expint(z::Union{Float64,Complex{Float64}})
    x² = real(z)^2
    y² = imag(z)^2
    if real(z) > 0 && x² + 0.233 * y² ≥ 7.84 # use cf expansion, ≤ 30 terms
        if (x² ≥ 546121) & (real(z) > 0) # underflow
            return zero(z)
        elseif x² + 0.401 * y² ≥ 58.0 # ≤ 15 terms
            if x² + 0.649 * y² ≥ 540.0 # ≤ 8 terms
                x² + y² ≥ 4e4 && return @E₁_cf64 z 4
                return @E₁_cf64 z 8
            end
            return @E₁_cf64 z 15
        end
        return @E₁_cf64 z 30
    else # use Taylor expansion, ≤ 37 terms
        r² = x² + y²
        return r² ≤ 0.36 ? (r² ≤ 2.8e-3 ? (r² ≤ 2e-7 ? @E₁_taylor64(z,4) :
                                                        @E₁_taylor64(z,8)) :
                                                        @E₁_taylor64(z,15)) :
                                                        @E₁_taylor64(z,37)
    end
end
expint(z::Union{T,Complex{T},Rational{T},Complex{Rational{T}}}) where {T <: Integer} = expint(float(z))

######################################################################
# exponential integral Eₙ(z)

function expint(n::Integer, z)
    if n == 1
        return expint(z)
    elseif n < 1
        # backwards recurrence from E₀ = e⁻ᶻ/z
        zinv = inv(z)
        e⁻ᶻ = exp(-z)
        Eᵢ = zinv * e⁻ᶻ
        for i in 1:-n
            Eᵢ = zinv * (e⁻ᶻ + i * Eᵢ)
        end
        return Eᵢ
    else
        # forwards recurrence from E₁
        e⁻ᶻ = exp(-z)
        Eᵢ = expint(z)
        Eᵢ *= !isinf(Eᵢ)
        for i in 2:n
            Eᵢ = (e⁻ᶻ - z * Eᵢ) / (i - 1)
        end
        return Eᵢ
    end
end


text:
# ============================================================
# Coffee Creamer of Quadrature
# ============================================================

"""
    slerp(pt1, pt2, n)

Spherical linear interpolation
"""
function slerp(pt1::T, pt2::T, n::Integer) where {T<:AbstractArray{<:Real,1}}
    if norm(pt1 - pt2) < 1e-10 # same points
        return repeat(pt1, 1, n) # return n copies of that point
    end

    omega = acos(dot(pt1, pt2))
    t = collect(range(0, stop = 1, length = n))

    return (sin.((1 .- t) * omega) / sin.(omega))' .* pt1 +
           (sin.(t * omega) / sin.(omega))' .* pt2
end

"""
    unique(Points, Triangles)

Cleaner for all duplicate (non unique) entries of quadrature points and triangles

* @args Points : quadrature points
* @args Triangles : triangulation
* @return xyz & triangulation : new quadrature points and triangulation

"""
function unique(
    Points::AbstractMatrix{T1},
    Triangles::AbstractMatrix{T2},
) where {T1<:Real,T2<:Integer}

    nPoints = size(Points)[2]
    nTriangles = size(Triangles)[2]

    map = collect(1:nPoints)
    unique = ones(Int64, nPoints)

    """
    map is a vector, where each entry contains the mapping to the ID in the unique vector we want to have.
    We start with map(i) = i and then iterate through all the points
    if it turns out, that point 10 == point 4, we will set map(10) = 4 (and map(4)=4)
    if then also point 20 == point 10 == point 4, we will set map(20) == 4 as well.
    unique simply stores a boolean as for whether or not the value is unique or a duplicate.
    This does however keep the first occurrence as unique.
    Meaning for a vector
    [1 2 3 3 4 5 3]
    unique would be
    [1 1 1 0 1 1 0]
    """

    for i = 1:nPoints
        for j = i+1:nPoints
            dist = sqrt(
                (Points[1, i] - Points[1, j])^2 +
                (Points[2, i] - Points[2, j])^2 +
                (Points[3, i] - Points[3, j])^2,
            )
            if (dist < 1e-6) # equal points
                map[j] = min(map[j], i) # take the minimal ID
                unique[j] = 0 # this point is no longer unique
            end
        end
    end

    """
    We are still not done
    Consider again the vector from above
    [1 2 3 3 4 5 3]
    map was at first
    [0 1 2 3 4 5 6]
    and is now
    [0 1 2 2 4 5 2]
    meaning that we have to renumber certain entries because
    we want to get
    [0 1 2 2 3 4 2] , i.e. the "real count
    """

    uniqueCounter = 1 # the "real" count
    for i = 1:nPoints
        if unique[i] == 1
            map[i] = uniqueCounter
            uniqueCounter += 1
        else
            map[i] = map[map[i]]
        end
    end

    # now we go through the triangles and update the IDs
    triangulation = zeros(Int64, 3, nTriangles)
    for i = 1:nTriangles
        for j = 1:3
            idx = Triangles[j, i]
            triangulation[j, i] = min(idx, map[idx])
        end
    end

    # count unique points
    nQuadPoints = 0
    for i = 1:nPoints
        nQuadPoints += (unique[i] == 1)
    end

    # write unique points to XYZ
    xyz = zeros(Float64, 3, nQuadPoints)
    counter = 1
    for i = 1:nPoints
        if unique[i] == 1
            xyz[:, counter] = Points[:, i]
            counter += 1
        end
    end

    return xyz, triangulation

end

function area(
    A::T,
    B::T,
    C::T,
    geometry = :plane::Symbol,
) where {T<:AbstractVector{<:Real}}

    if geometry == :plane
        alpha = angle(B, A, C)
        lb = norm(B - A)
        la = norm(C - A)
        return 0.5 * sin(alpha) * lb * la
    elseif geometry == :sphere
        alpha = angle(B, A, C, :sphere)
        beta = angle(C, B, A, :sphere)
        gamma = angle(A, C, B, :sphere)
        #https://en.wikipedia.org/wiki/Spherical_trigonometry#Area_and_spherical_excess
        #Excess = max(0.0,alpha+beta+gamma-pi)
        #R = 1; # We always consider the uni sphere
        #return Excess*R*R;
        return alpha + beta + gamma - pi
    else
        throw("geometry has to be planar or sphere")
    end

end

"""
Args order (B,A,C) isn't a mistake

"""
function angle(
    B::T,
    A::T,
    C::T,
    geometry = :plane::Symbol,
) where {T<:AbstractVector{<:Real}}

    if geometry == :plane
        u, v = A - B, C - A
        return acos(dot(u, v) / norm(u, 2) / norm(v, 2))
    elseif geometry == :sphere
        c = distance(B, A, :sphere)
        b = distance(A, C, :sphere)
        a = distance(C, B, :sphere)
        if min(a, b, c) < 1e-10
            return 0
        end
        # https://en.wikipedia.org/wiki/Spherical_trigonometry#Cosine_rules_and_sine_rules
        tmp = (cos(a) - cos(b) * cos(c)) / (sin(b) * sin(c))
        # numerical error may lead that tmp is not inside [-1,1] but inside [-1-eps,1+eps]
        tmp = max(-1.0, min(1.0, tmp))

        angle = acos(tmp)
        return angle
    else
        throw("geometry has to be planar or sphere")
    end

end

function distance(
    v1::T,
    v2::T,
    geometry = :plane::Symbol,
) where {T<:AbstractVector{<:Real}}

    if geometry == :plane
        return norm(v1 - v2)
    elseif geometry == :sphere
        if sum(abs2, v1 - v2) < 1e-5
            return norm(v1 - v2)
        end
        #https://en.wikipedia.org/wiki/Great-circle_distance
        #return atan(norm(cross(v1, v2), 2) / dot(v1, v2))
        return acos(max(-1.0, dot(v1, v2)))
    else
        throw("geometry has to be planar or sphere")
    end

end

function muphi_xyz!(muphi::T, xyz::T) where {T<:AbstractMatrix{<:Real}}
    n = size(xyz, 1)
    for i = 1:n
        xyz[i, 1] = sqrt(1 - muphi[i, 1]^2) * cos(muphi[i, 2])
        xyz[i, 2] = sqrt(1 - muphi[i, 1]^2) * sin(muphi[i, 2])
        xyz[i, 3] = muphi[i, 1]
    end
end

function xyz_muphi!(xyz::T, muphi::T) where {T<:AbstractMatrix{<:Real}}
    n = size(xyz, 1)
    for i = 1:n
        muphi[i, 1] = xyz[i, 3]
        muphi[i, 2] = atan(xyz[i, 2], xyz[i, 1])
    end
end


text:
#! /usr/bin/julia

# Rosetta Code, Loops/Continue

for i in 1:10
    print(i)
    if i%5 == 0
        println()
        continue
    end
    print(", ")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
module mopkg

using Calculus

"minimum_nr() returns the minimum value of the given function, calculated by the Newton-Raphson method."

function minimum_nr(f, x, eps = 1e-5)
	fx = f(x)
	i = 0
	d = Calculus.derivative(f,x)
	dd = Calculus.second_derivative(f,x)
	h = d/dd
	while abs(h) >= eps
		d = Calculus.derivative(f,x)
		dd = Calculus.second_derivative(f,x)
		h = d/dd
		x = x - h
		fx = f(x)
		i = i + 1
		#println("x: ", x, "   f(x): ", fx, "   step: ", h, "   iteration: ", i)
	end
	return fx, x
end

export minimum_nr

end # module


text:
arg_dict = Dict{Any, Any}(
"physics" => "Euler",
"dimensions" => 2,
"use_DG" => true,
"operator_type" => "SBPOmega",
"run_type" => 1,
"real_time" => true,
"order" => 1,
"Flux_name" => "IRFlux",
"IC_name" => "ICPeriodicMMS",
"volume_integral_type" => 2,
"face_integral_type" => 2,
"Volume_flux_name" => "IRFlux",
"FaceElementIntegral_name" => "ESLW2FaceIntegral",
"vortex_x0" => 5.0,
"numBC" => 0,
"BC1_name" => "PeriodicMMSBC",
"CFL" => 0.25,
"t_max" => 10.0,
"SRCname" => "SRCPeriodicMMS",
"smb_name" => "SRCMESHES/square_curve.smb",
"dmg_name" => ".null",
"itermax" => 10000,
"output_freq" => 50,
"write_vis" => true,
"do_postproc" => true,
"exact_soln_func" => "ICPeriodicMMS",
"res_abstol" => 1e-8,
"res_reltol" => 1e-10,
"step_tol" => 1e-10,
"solve" => false,
)


text:
include("SortSparse.jl")
include("KoLesky.jl")
include("Cholesky.jl")
include("CovFuncs.jl")
include("Utils.jl")
using LinearAlgebra
using Base.Threads
using Distributions
using Distances
using Random
using Match
using JLD
using BlockArrays
using ArgParse
Random.seed!(123)

s = ArgParseSettings()
@add_arg_table s begin
  "maternOrder"
  help = "The order of the matern kernel to be used. Should be in (\"12\", \"32\", \"52\")"
  required = true
end
maternOrder = parse_args(s)["maternOrder"]

NTrain = 5000
NTest = 100
NSamples = 10
N = NTrain + NTest
d = 2
l = 0.5
λ = 1.0
ρList = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]

meanRMSEPredFirst = zeros(length(ρList))
stdRMSEPredFirst = zeros(length(ρList))
meanRMSEPredLast = zeros(length(ρList))
stdRMSEPredLast = zeros(length(ρList))
meanRMSENoPred = zeros(length(ρList))
stdRMSENoPred = zeros(length(ρList))

# Creating uniform random data
xTrain = rand(d, NTrain)
xTest = 1.1 * ones(d,NTest) .+ 0.1 * rand(d, NTest)
x = hcat(xTrain, xTest)

covfunc12!(r) = matern12!(r,l)
covfunc32!(r) = matern32!(r,l)
covfunc52!(r) = matern52!(r,l)

@match maternOrder begin
  "12" => (covfunc! = covfunc12!)
  "32" => (covfunc! = covfunc32!)
  "52" => (covfunc! = covfunc52!)
end

for (ρInd, ρ)  in enumerate(ρList)

  global meanRMSEPredFirst
  global stdRMSEPredFirst
  global meanRMSEPredLast
  global stdRMSEPredLast
  global meanRMSENoPred
  global stdRMSENoPred

    for k = 1 : NSamples

    # =====================================
    # Now construct the proper covariance matrix
    # =====================================

    Θ = covfunc!(pairwise(Euclidean(), x, dims=2))
    chol_Θ = cholesky(Matrix(Θ))

    Θ = PseudoBlockArray(Θ, [NTrain, NTest], [NTrain, NTest])
    chol_Θ_train = cholesky(Θ[Block(1,1)])
    # =====================================
    # Starting with interpolation:
    # =====================================

    # Creating reference solution
    y = chol_Θ.L * randn(N)

    # Creating data
    y = PseudoBlockArray(y, [NTrain, NTest])

    # Computing the true posterior mean
    μ = Θ[Block(2,1)] * (chol_Θ_train \ y[Block(1)])

    # only needs to be done for one sample
    Σ = Θ[Block(2,2)] - Θ[Block(2,1)] * (chol_Θ_train \ Θ[Block(1,2)])

    # =====================================
    # Computing posterior mean for prediction variables last
    # =====================================

    # reordering points and forming skeletons
    colptr, rowval, P, revP, distances = sortSparseRev(xTrain, xTest, ρ, 1)
    xOrd = x[:, P]
    xOrdTest = xOrd[:, 1 : NTest]
    xOrdTrain = xOrd[:, (NTest + 1) : end]
    skeletons = construct_skeletons(colptr, rowval, distances, λ)
    μOrd = μ[P[1 : NTest] .- NTrain]
    ΘOrd = Θ[P, P]
    ΣOrd = Σ[P[1 : NTest] .- NTrain, P[1 : NTest] .- NTrain]
    yOrd = PseudoBlockArray(y[P], [NTest, NTrain])

    # Creating vecchia approximation
    L = assembleL(skeletons, xOrd, covfunc!)
    L22 = L[(NTest + 1) : end, (NTest + 1) : end]
    μPredLast = (L' \ ( L \ vcat(zero(y[Block(2)]), L22 * (L22' * (yOrd[Block(2)])))))[1:NTest]

    ΣPredLast = inv(Matrix(L[1 : NTest, 1 : NTest] * L[1 : NTest, 1 : NTest]'))

    ΘPredLast = inv(Matrix(L * L'))

    #added rescaling with std
    meanRMSEPredLast[ρInd] += sum((μOrd - μPredLast).^2 ./ diag(ΣOrd))

    # only needs to be done once
    stdRMSEPredLast[ρInd] = sqrt(mean((sqrt.(diag(ΣOrd)) - sqrt.(diag(ΣPredLast))).^2 / diag(ΣOrd)))

    # =====================================
    # Computing posterior mean for prediction variables first
    # =====================================

    # reordering points and forming skeletons
    colptr, rowval, P, revP, distances = sortSparseRev(xTrain, ρ, 1)
    xOrdTrain = xTrain[:, P]
    xOrd = hcat(xOrdTrain, xTest)
    skeletons = construct_skeletons(colptr, rowval, distances, λ)
    μOrd = μ
    ΘOrd = Θ[vcat(P, (NTrain + 1) : N), vcat(P, (NTrain + 1) : N)]
    ΣOrd = Σ
    yOrd = y[Block(1)][P]

    μPredFirst, ΣPredFirst = predict(skeletons, (NTrain+1): N, yOrd, xOrd, covfunc!)

    # added scaling by std
    meanRMSEPredFirst[ρInd] += sum((μOrd - μPredFirst).^2 ./ diag(ΣOrd))
    stdRMSEPredFirst[ρInd] = sqrt(mean((sqrt.(diag(ΣOrd)) - sqrt.(ΣPredFirst)).^2 ./ diag(ΣOrd)))

    # =====================================
    # Don't include prediction variables at all
    # =====================================

    # reordering points and forming skeletons
    colptr, rowval, P, revP, distances = sortSparseRev(xTrain, ρ, 1)
    xOrdTrain = xTrain[:, P]
    xOrd = hcat(xOrdTrain, xTest)
    skeletons = construct_skeletons(colptr, rowval, distances, λ)
    L = assembleL(skeletons, xOrdTrain, covfunc!)
    @show nnz(L) * 2 / NTrain / (NTrain + 1)
    μOrd = μ
    ΘOrd = Θ[vcat(P, (NTrain + 1) : N), vcat(P, (NTrain + 1) : N)]
    ΣOrd = Σ
    yOrd = y[Block(1)][P]

    TMat = pairwise(Euclidean(), xOrdTrain, xTest, dims=2)
    covfunc!(TMat)
    CMat = pairwise(Euclidean(), xTest, dims=2)
    covfunc!(CMat)

    ΣNoPred = CMat - TMat' * L * L' * TMat

    ΘNoPred = inv(Matrix(L * L'))
    ΘNoPred = vcat(ΘNoPred, TMat')
    ΘNoPred = hcat(ΘNoPred, vcat(TMat, CMat))

    stdRMSENoPred[ρInd] = sqrt.(mean((sqrt.(diag(ΣOrd)) - sqrt.(max.(diag(ΣNoPred), 0.0))).^2 ./ diag(ΣOrd)))

    μNoPred = vec(yOrd' * L * L' * TMat)
    meanRMSENoPred[ρInd] += sum((μOrd - μNoPred).^2 ./ diag(ΣOrd))
  end
end
meanRMSENoPred /= NSamples * NTest
meanRMSENoPred .= sqrt.(meanRMSENoPred)
meanRMSEPredFirst /= NSamples * NTest
meanRMSEPredFirst .= sqrt.(meanRMSEPredFirst)
meanRMSEPredLast /= NSamples * NTest
meanRMSEPredLast .= sqrt.(meanRMSEPredLast)

save("./out/jld/IncludePredFarVaryRho$maternOrder.jld",
     "meanRMSEPredFirst", meanRMSEPredFirst,
     "stdRMSEPredFirst", stdRMSEPredFirst,
     "meanRMSEPredLast", meanRMSEPredLast,
     "stdRMSEPredLast", stdRMSEPredLast,
     "meanRMSENoPred", meanRMSENoPred,
     "stdRMSENoPred", stdRMSENoPred,
     "xTrain", xTrain,
     "xTest", xTest,
     "ρList", ρList)


text:
__precompile__()
module Annealing
    using Base.Dates
    using Tsp.Model
    using Tsp.Debug

    function find_optimal(;
        grid = (:required, "grid"),
        s0 = (:required, "s0"),
        T0 = (:required, "temp0"),
        limit::Int = 1000,
        starttime = (:required, "starttime"),
        timelimit = -1,
        cooling_factor=0.997
    )

        check_args(grid, s0, T0, starttime)
        current::Route = s0
        best::Route = current
        best_round::Int = 0
        last_round::Int = 0
        T::Float64 = T0

        no_progress = 0
        cooldowns = 0
        reheats = 0

        init_reheat_threshold = 50000
        T0 += length(grid.cities)/1000 * 10
        cooling_factor -= length(grid.cities)/1000 * .003
        # if length(grid.cities) > 10000
        #     init_reheat_threshold = 100000
        # elseif length(grid.cities > 8000)

        reheat_threshold = init_reheat_threshold

        for step in 0:limit
            if no_progress >= reheat_threshold && reheats < 30
                reheat_threshold = init_reheat_threshold
                no_progress = 0
                current = best
                reheats += 1
                @info @printf("[%10d] [T=%8e] No progress - restarting from best known position.\n", step, T)

                T = T0/(reheats+5)
                if length(grid.cities) > 4000
                    cooling_factor -= .0001
                end
                # T=T0/(reheats*5)
                # T=temp(cooldowns - Int(round(cooldowns/2)), T0)
            end
            candidate, i, j  = random_neighbour(current, T, grid)

            diff = (candidate.score - s0.score) / s0.score * 100

            p::Float64 = probability(Float64(best.score), Float64(candidate.score), T)

            if p >= rand(0.0:0.0001:1.0)
                @info @printf("[%10d] [T=%8e] [p=%g] Move: (%5d, %5d). New score: %.8f (%+.4f%%)\n",
                        step, T, p, i, j, candidate.score, diff)
                current = candidate
                T = temp(T, cooling_factor)
                cooldowns += 1
            elseif candidate.score < current.score
                # @debug @printf("[%10d] [T=%8e] Switching to a better route: %.8f (%+.4f%%)\n",
                #     step, T, candidate.score, diff)
                current = candidate
                T = temp(T, cooling_factor)
                cooldowns += 1
                last_progress = step
                # reheat_threshold += 1000
            else
                no_progress += 1
            end
            if current.score < best.score
                reheat_threshold = init_reheat_threshold*5
                @log @printf("[%10d] [T=%8e] New best: %.8f (%+.4f%%)\n",
                  step, T, candidate.score, diff)
                best = candidate
                best_round = step
                no_progress = 0
            end

            done, elapsed = timeout(timelimit, starttime)
            if done
                @log @printf("[%10d] Timeout (%s)\n", step, elapsed)
                return best
            end

            if T <= 1.0e-6
                @log @printf("Temperature reached 0, stopping search in round %d\n", step)
                last_round = step
                break
            end

            if step > 0 && step % 50000 == 0
                @info @printf("[%10d] [T=%8e] [cd=%5d] [rh=%3d] last candidate: %.8f (%+.4f%%)\n",
                    step, T, cooldowns, reheats, candidate.score, diff)
            end
        end

        @log @printf("Solution found in round %d\n", best_round)
        return best
    end

    function timeout(timeLimit::Dates.Period, startTime::DateTime)
        elapsed = now()-startTime
        return elapsed > Millisecond(timeLimit) - Millisecond(1000), elapsed
    end

    timeout(timeLimit::Int, startTime::DateTime) = false, 0

    function random_neighbour(origin::Route, T::Float64, grid::Grid)::Tuple{Route,Int,Int}
        s1 = 1 + (rand(0:0.0001:0.01)) * origin.score
        s2 = 1 - (rand(0:0.0001:0.01)) * origin.score
        # if rand(0:0.0001:1.0) >= probability(s1, s2, T)
            _random_neighbour(origin)
        # else
            # _chunk_neighbour(origin, grid)
        # end
    end

    function _random_neighbour(origin::Route)::Tuple{Route,Int,Int}
        newroute = Route(origin)

        i = rand(2:length(origin)-1)
        j = rand(2:length(origin)-1)
        while j == i
            j = rand(2:length(origin)-1)
        end

        _reverse!(newroute, i, j)
        return newroute,i,j
    end

    function _chunk_neighbour(origin::Route, grid::Grid)::Tuple{Route,Int,Int}
        newroute = Route(origin)
        i = rand(2:length(origin)-1)
        src_city = newroute.cities[i]
        chunk = grid.chunks[src_city.chunk_x, src_city.chunk_y]
        if length(chunk.cities) < 2
            return _random_neighbour(origin)
        end

        dst_city = rand(chunk.cities)
        while dst_city.id == src_city.id || dst_city.id == 1
            dst_city = rand(chunk.cities)
        end

        j = findfirst(city->city.id == dst_city.id, newroute.cities)
        println("-------")
        println(grid.cities[i])
        println(grid.cities[j])
        println("-------")
        _swap!(newroute, i, j)
        return newroute, i,j
    end

    function _reverse!(route, i, j)
        idx = i < j ? (i:j) : (j:i)
        route.score -= distance(route.cities[first(idx)], route.cities[first(idx)-1]) + distance(route.cities[last(idx)], route.cities[last(idx)+1])
        route.cities[idx] = route.cities[reverse(idx)]
        route.score += distance(route.cities[first(idx)], route.cities[first(idx)-1]) + distance(route.cities[last(idx)], route.cities[last(idx)+1])
        route
    end

    function _swap!(route, i, j)
        route.score -= distance(route.cities[i], route.cities[i-1])
                    + distance(route.cities[i], route.cities[i+1])
                    + distance(route.cities[j], route.cities[j-1])
                    + distance(route.cities[j], route.cities[j+1])

        route.cities[i], route.cities[j] = route.cities[j], route.cities[i]

        route.score += distance(route.cities[i], route.cities[i-1])
                    + distance(route.cities[i], route.cities[i+1])
                    + distance(route.cities[j], route.cities[j-1])
                    + distance(route.cities[j], route.cities[j+1])
        route
    end

    """ get new temperature based on current value """
    temp(prev::Float64, rate::Float64)::Float64 = prev * rate

    """ get the temperature after n cooldowns """
    temp(n::Int, rate::Float64, T0::Float64) = n > 0 ? T0 * rate^n : T0

    # """
    #     Get the estimated temperature after given cooldown
    # """
    # temp(cooldown::Int, t0::Float64) = cooldown > 0 : t0 * rate^cooldown

    function probability(old::Float64, next::Float64, temp::Float64)
        delta =  next - old
        return 1.0/(1.0 + e ^ (delta / temp))
    end

    function check_args(args...)
        for arg in args
            if isa(arg, Tuple) && arg[1] == :required
                throw(ArgumentError("$(arg[2]) must be specified"))
            end
        end
    end

end