text:
"""
Collection of utilities for working with Requests data
"""
module Requests

using Genie, Genie.Router, Genie.Input
using HTTP

export jsonpayload, rawpayload, filespayload, postpayload, getpayload, request, matchedroute, matchedchannel
export infilespayload, download, filename, payload, read

"""
    jsonpayload()

Processes an `application/json` `POST` request.
If it fails to successfully parse the `JSON` data it returns `nothing`. The original payload can still be accessed invoking `rawpayload()`
"""
@inline function jsonpayload()
  @params(Genie.PARAMS_JSON_PAYLOAD)
end

"""
    jsonpayload(::Type{T}) where {T}

Processes an `application/json` `POST` request attempting to convert the payload into a value of type `T`.
If it fails to successfully parse and convert the `JSON` data, it throws an exception. The original payload can still be accessed invoking `rawpayload()`
"""
@inline function jsonpayload(::Type{T})::T where {T}
  @params(Genie.PARAMS_JSON_PAYLOAD)::T
end

"""
    rawpayload() :: String

Returns the raw `POST` payload as a `String`.
"""
@inline function rawpayload() :: String
  @params(Genie.PARAMS_RAW_PAYLOAD)
end

"""
    filespayload() :: Dict{String,HttpFile}

Collection of form uploaded files.
"""
@inline function filespayload() :: Dict{String,HttpFile}
  @params(Genie.PARAMS_FILES)
end

"""
    filespayload(filename::Union{String,Symbol}) :: HttpFile

Returns the `HttpFile` uploaded through the `key` input name.
"""
@inline function filespayload(key::Union{String,Symbol}) :: HttpFile
  @params(Genie.PARAMS_FILES)[string(key)]
end

"""
    infilespayload(key::Union{String,Symbol}) :: Bool

Checks if the collection of uploaded files contains a file stored under the `key` name.
"""
@inline function infilespayload(key::Union{String,Symbol}) :: Bool
  haskey(filespayload(), string(key))
end

"""
    Base.download(file::HttpFile, name::String = file.name)

Saves uploaded `HttpFile` `file` to local storage under the `name` filename.
"""
@inline function Base.download(file::HttpFile, name::String = file.name) :: Int
  write(name, IOBuffer(file.data))
end

"""
    read(file::HttpFile)

Returns the content of `file` as string.
"""
@inline function Base.read(file::HttpFile, ::Type{String}) :: String
  file.data |> String
end

"""
    filename(file::HttpFile) :: String

Original filename of the uploaded `HttpFile` `file`.
"""
@inline function filename(file::HttpFile) :: String
  file.name
end

"""
    postpayload() :: Dict{Symbol,Any}

A dict representing the POST variables payload of the request (corresponding to a `form-data` request)
"""
@inline function postpayload() :: Dict{Symbol,Any}
  @params(Genie.PARAMS_POST_KEY)
end

"""
    postpayload(key::Symbol) :: Any

Returns the value of the POST variables `key`.
"""
@inline function postpayload(key::Symbol)
  postpayload()[key]
end

"""
    postpayload(key::Symbol, default::Any)

Returns the value of the POST variables `key` or the `default` value if `key` is not defined.
"""
@inline function postpayload(key::Symbol, default::Any)
  haskey(postpayload(), key) ? postpayload(key) : default
end

"""
    getpayload() :: Dict{Symbol,Any}

A dict representing the GET/query variables payload of the request (the part correspoding to `?foo=bar&baz=moo`)
"""
@inline function getpayload() :: Dict{Symbol,Any}
  @params(Genie.PARAMS_GET_KEY)
end

"""
    getpayload(key::Symbol) :: Any

The value of the GET/query variable `key`, as in `?key=value`
"""
@inline function getpayload(key::Symbol) :: Any
  getpayload()[key]
end

"""
    getpayload(key::Symbol, default::Any) :: Any

The value of the GET/query variable `key`, as in `?key=value`. If `key` is not defined, `default` is returned.
"""
@inline function getpayload(key::Symbol, default::Any) :: Any
  haskey(getpayload(), key) ? getpayload(key) : default
end

"""
    request() :: HTTP.Request

Returns the raw HTTP.Request object associated with the request.
"""
@inline function request() :: HTTP.Request
  @params(Genie.PARAMS_REQUEST_KEY)
end

"""
    payload() :: Any

Utility function for accessing the `@params` collection, which holds the request variables.
"""
@inline function payload() :: Any
  @params
end

"""
    payload(key::Symbol) :: Any

Utility function for accessing the `key` value within the `@params` collection of request variables.
"""
@inline function payload(key::Symbol) :: Any
  @params()[key]
end

"""
    payload(key::Symbol, default_value::T) :: Any

Utility function for accessing the `key` value within the `@params` collection of request variables.
If `key` is not defined, `default_value` is returned.
"""
@inline function payload(key::Symbol, default_value::T)::T where {T}
  haskey(@params(), key) ? @params()[key] : default_value
end

"""
    matchedroute() :: Route

Returns the `Route` object which was matched for the current request.
"""
@inline function matchedroute() :: Genie.Router.Route
  @params(Genie.PARAMS_ROUTE_KEY)
end

"""
    matchedchannel() :: Channel

Returns the `Channel` object which was matched for the current request.
"""
@inline function matchedchannel() :: Genie.Router.Channel
  @params(Genie.PARAMS_CHANNELS_KEY)
end

end


text:
let
        tot = 0
        prev = Dict()
        foundIt = false
        loopCount = 0
        while !foundIt
            loopCount += 1
            println(loopCount)
            open("data.txt") do file
                for line in eachline(file)
                    operator = line[1]
                    operand = parse(Int, line[2:end])
                    if operator == '+'
                        tot += operand
                    else
                        tot -= operand
                    end
                    if haskey(prev, tot)
                        foundIt = true
                        break;
                    end

                    prev[tot] = 1
                end
            end

        end
        println(tot)

end

ext: jl


text:
using DistributionVectors
using Documenter

# allow plots on headless server
# https://juliadocs.github.io/Documenter.jl/stable/man/syntax/#@example-block
ENV["GKSwstype"] = "100"

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

makedocs(;
    modules=[DistributionVectors],
    authors="Thomas Wutzler",
    repo="https://github.com/bgctw/DistributionVectors.jl/blob/{commit}{path}#{line}",
    sitename="DistributionVectors.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://bgctw.github.io/DistributionVectors.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Types" => "distributionvector.md",
        "Aggregation" => "sumdist.md",
        "Univariate" => [
            "Overview" => "univariate.md",
            "Normal" => "normal.md",
            "LogNormal" => "lognormal.md",
        ],
    ],
    doctest=false,
)

function tmpf()
    #https://discourse.julialang.org/t/deploy-documenter-docs-locally/5463/7
    #ENV["TRAVIS_REPO_SLUG"] = "github.com/bgctw/DistributionVectors.jl.git"
    ENV["TRAVIS_REPO_SLUG"] = "."
    ENV["TRAVIS_PULL_REQUEST"] = "false"
    ENV["TRAVIS_OS_NAME"] = "linux"
    ENV["TRAVIS_JULIA_VERSION"] = "1.5"
    ENV["TRAVIS_TAG"] = ""
    ENV["TRAVIS_BRANCH"] = "main"
    # ENV["DOCUMENTER_KEY"] = chomp(open(f->read(f, String), "docs/secrets/DOCUMENTER_KEY"));
    deploydocs(deps = nothing, make = nothing,
        repo = "github.com/bgctw/DistributionVectors.jl.git",
        target = "build",
        branch = "gh-pages",
        devbranch = ENV["TRAVIS_BRANCH"],
    )
    # does not work, because the build folder is not on github
end

deploydocs(;
    repo="github.com/bgctw/DistributionVectors.jl",
    devbranch="main",
)

function deploydocs_fromlocal()
    # pages directory must be setup to track gh-pages branch
    cd("../DistributionVectors_pages/")
    run(`git checkout gh-pages`)
    run(`git pull`)
    src = "../DistributionVectors/docs/build"
    dst = "dev"
    readdir(src)
    readdir(dst)
    cp("$dst/siteinfo.js","$src/siteinfo.js")
    cp(src, dst, force=true)
    run(`git add .`)
    run(`git commit -m"local doc generation"`)
    run(`git push`)
    cd("../DistributionVectors")
end


text:
using BoltzmannMachines
import BoltzmannMachinesPlots # only needed for plotting capabilities

# ==============================================================================
# Binary data: BernoulliRBM and BasicDBM
# ------------------------------------------------------------------------------

# Generate data and split it in training and test data set
nsamples = 100;
nvariables = 36;
using Random
Random.seed!(1);
x = barsandstripes(nsamples, nvariables);
x, xtest = splitdata(x, 0.3);
datadict = DataDict("Training data" => x, "Test data" => xtest);

# Fit a BernoulliRBM and collect monitoring data during training.
# For small models, the model's loglikelihood can be computed exactly.
monitor = Monitor();
rbm = fitrbm(x; nhidden = 12,
      epochs = 80, learningrate = 0.006,
      monitoring = (rbm, epoch) -> begin
            monitorexactloglikelihood!(monitor, rbm, epoch, datadict)
            monitorreconstructionerror!(monitor, rbm, epoch, datadict)
      end);

# Plot the collected data (requires Gadfly)
BoltzmannMachinesPlots.plotevaluation(monitor, monitorexactloglikelihood)
BoltzmannMachinesPlots.plotevaluation(monitor, monitorreconstructionerror)

# Simplified monitoring:
monitor, rbm = monitored_fitrbm(x; nhidden = 12,
      epochs = 80, learningrate = 0.006,
      monitoringdata = datadict,
      monitoring = [monitorexactloglikelihood!, monitorreconstructionerror!]);

Random.seed!(12);

# DBM-Fitting: Pretraining and Fine-Tuning combined in one function
dbm = fitdbm(x, nhiddens = [6;2], epochs = 20, learningrate = 0.05);

# .. with extensive monitoring
monitor = Monitor(); monitor1 = Monitor(); monitor2 = Monitor();
dbm = fitdbm(x, epochs = 20, learningrate = 0.05,
      monitoringdatapretraining = datadict,
      pretraining = [
            TrainLayer(nhidden = 6, monitoring = (rbm, epoch, datadict) ->
                  monitorreconstructionerror!(monitor1, rbm, epoch, datadict));
            TrainLayer(nhidden = 2, monitoring = (rbm, epoch, datadict) ->
                  monitorreconstructionerror!(monitor2, rbm, epoch, datadict))
            ],
      monitoring = (dbm, epoch) ->
            monitorexactloglikelihood!(monitor, dbm, epoch, datadict));

# Monitoring plots
BoltzmannMachinesPlots.plotevaluation(monitor1)
BoltzmannMachinesPlots.plotevaluation(monitor2)
BoltzmannMachinesPlots.plotevaluation(monitor)

# Evaluate final result with exact computation of likelihood
exactloglikelihood(dbm, xtest)

# Simplified monitoring
monitors, dbm = monitored_fitdbm(x; nhiddens = [6,2],
      epochs = 20, learningrate = 0.05,
      monitoringpretraining = monitorreconstructionerror!,
      monitoring = monitorexactloglikelihood!,
      monitoringdata = datadict);

# If the model has more parameters, in this case more hidden nodes,
# the exact calculation is not feasible any more:
# We need to calculate the likelihood using AIS.

Random.seed!(12);
monitor = Monitor();
rbm = fitrbm(x, nhidden = 36, epochs = 100,
      learningrates = [0.008*ones(80); 0.001*ones(20)],
      monitoring = (rbm, epoch) ->
            monitorloglikelihood!(monitor, rbm, epoch, datadict));

BoltzmannMachinesPlots.plotevaluation(monitor, monitorloglikelihood; sdrange = 3.0)

# Simplified monitoring:
monitor, rbm = monitored_fitrbm(x, nhidden = 36, epochs = 100,
      learningrates = [0.008*ones(80); 0.001*ones(20)],
      monitoring = monitorloglikelihood!,
      monitoringdata = datadict);

# In the DBM, the estimation of the log likelihood is much slower than in the
# RBM because it needs a AIS run for each sample.
# The lower bound of the log probability is faster to estimate:
# (But on the other hand, in some cases it might be too conservative to be
# useful as it is only a lower bound.)

Random.seed!(2);
monitor = Monitor();
dbm = stackrbms(x; nhiddens = [36;10;5], predbm = true, learningrate = 0.05);
traindbm!(dbm, x; epochs = 50, learningrate = 0.008,
      monitoring = (rbm, epoch) ->
            monitorlogproblowerbound!(monitor, rbm, epoch, datadict));

BoltzmannMachinesPlots.plotevaluation(monitor, monitorlogproblowerbound; sdrange = 3.0)

# Simplified monitoring:
monitor, dbm = monitored_traindbm!(dbm, x; epochs = 50, learningrate = 0.008,
      monitoring = monitorlogproblowerbound!,
      monitoringdata = datadict);

# Evaluate final result with AIS-estimated likelihood
loglikelihood(dbm, xtest)

# ==============================================================================
# Categorical data: Softmax0BernoulliRBM
# ------------------------------------------------------------------------------

# Some data set with values in the categories {0, 1, 2}
# Mixing variables with a different number of categories is also possible
# (see also example below for a MultimodalDBM with binary and categorical
# variables).
xcat = map(v -> max(v, 0),
      barsandstripes(100, 4) .+ barsandstripes(100, 4) .- barsandstripes(100, 4));

# Encode values for Softmax0BernoulliRBM:
# Each of the variables, which have three categories, are translated into
# two binary variables.
# This encoding is to similiar to the one-hot encoding with the deviation
# that a zero is encoded as all-zeros.
# (That way, the encoding of binary variables is not changed.)
x01 = oneornone_encode(xcat, 3) # 3 categories
rbm = fitrbm(x01, rbmtype = Softmax0BernoulliRBM, categories = 3,
      nhidden = 4, epochs = 20);

# The monitoring is omitted here for brevity:
# The same monitoring and training options as those for BernoulliRBMs
# are also available for Softmax0BernoulliRBMs.

# A DBM with a Softmax0BernoulliRBM in the first layer and binary hidden layers:
dbm = fitdbm(x01, pretraining = [
         TrainLayer(rbmtype = Softmax0BernoulliRBM, categories = 3, nhidden = 4);
         TrainLayer(nhidden = 4);
         TrainLayer(nhidden = 2)])

# Getting samples with values in the original sample space:
oneornone_decode(samples(dbm, 5), 3)

# ==============================================================================
# Mixing binary and categorical data: MultimodalDBM
# ------------------------------------------------------------------------------

# x1: binary data
x1 = barsandstripes(100, 4);
# x2: values with categories {0, 1, 2}
x2 = map(v -> max(v, 0), x1 .+ barsandstripes(100, 4) .- barsandstripes(100, 4));
x2 = oneornone_encode(x2, 3)
x = hcat(x1, x2);
x, xtest = splitdata(x, 0.3);
datadict = DataDict("Training data" => x, "Test data" => xtest);
monitor1 = Monitor(); monitor2 = Monitor(); monitor3 = Monitor(); monitor = Monitor();
dbm = fitdbm(x; epochspretraining = 50, epochs = 15,
      learningratepretraining = 0.05, learningrate = 0.1,
      monitoringdatapretraining = datadict,
      pretraining = [
            TrainPartitionedLayer([
               TrainLayer(nhidden = 5, nvisible = 4,
                     monitoring = (rbm, epoch, datadict) ->
                           monitorreconstructionerror!(monitor1, rbm, epoch, datadict));
               TrainLayer(nhidden = 7,
                     rbmtype = Softmax0BernoulliRBM, categories = 3,
                     monitoring = (rbm, epoch, datadict) ->
                           monitorreconstructionerror!(monitor2, rbm, epoch, datadict))
            ]),
            TrainLayer(nhidden = 6,
                  monitoring = (rbm, epoch, datadict) ->
                        monitorreconstructionerror!(monitor3, rbm, epoch, datadict))
      ],
      monitoring =
            (dbm, epoch) -> monitorlogproblowerbound!(monitor, dbm, epoch, datadict));

# The same with simplified monitoring:
Random.seed!(1)
monitors, dbm = monitored_fitdbm(x, epochspretraining = 50, epochs = 15,
      pretraining = [
            TrainPartitionedLayer([
                  TrainLayer(nhidden = 5, nvisible = 4);
                  TrainLayer(nhidden = 7,
                        rbmtype = Softmax0BernoulliRBM, categories = 2)]),
            TrainLayer(nhidden = 6)
      ],
      learningratepretraining = 0.05,
      learningrate = 0.1,
      monitoringdata = datadict,
      monitoringpretraining = monitorreconstructionerror!,
      monitoring = monitorlogproblowerbound!);

# Reconstructionerror of BernoulliRBM in first layer:
BoltzmannMachinesPlots.plotevaluation(monitors[1][1])
# Reconstructionerror of Softmax0BernoulliRBM:
BoltzmannMachinesPlots.plotevaluation(monitors[1][2])
# Reconstructionerror of second layer:
BoltzmannMachinesPlots.plotevaluation(monitors[2])
# Lower bound of likelihood during fine-tuning:
BoltzmannMachinesPlots.plotevaluation(monitors[3])

# It is also possible to combine binary and categorical variables in one RBM,
# as a Softmax0BernoulliRBM with variables having two categories is equivalent
# to a BernoulliRBM.
# Here we do not have a partition then. We can use one Softmax0BernoulliRBM for
# the visible layer and have to specify the number of categories for each of the
# variables separately:

categories = [fill(2, 4); fill(3, 4)];
dbm = fitdbm(x;
      pretraining = [
         TrainLayer(nhidden = 7, categories = categories,
               rbmtype = Softmax0BernoulliRBM);
         TrainLayer(nhidden = 6)
      ])

# ==============================================================================
# Real valued data: Intensities, GaussianBernoulliRBM or GaussianBernoulliRBM2
# ------------------------------------------------------------------------------

# Use "iris" dataset as example data for continuous data
using DelimitedFiles
x = convert(Matrix{Float64},
      readdlm(joinpath(dirname(pathof(BoltzmannMachines)), "..",
            "test/data/iris.csv"), ',', header = true)[1][:,1:4]);
Random.seed!(12);
x, xtest = splitdata(x, 0.3);
datadict = DataDict("Training data" => x, "Test data" => xtest);

# Using the intensities transformation and a standard BernoulliRBM
xintensities, xtrafo = intensities_encode(x)
rbm = fitrbm(xintensities)
# Transform samples back into the original sample space:
intensities_decode(samples(rbm, 5), xtrafo)

# Example using a GaussianBernoulliRBM
Random.seed!(1);
monitor = Monitor();
rbm = fitrbm(x, rbmtype = GaussianBernoulliRBM,
      nhidden = 3, epochs = 80, learningrate = 0.0005,
      monitoring = (rbm, epoch) ->
            monitorexactloglikelihood!(monitor, rbm, epoch, datadict));

BoltzmannMachinesPlots.plotevaluation(monitor, monitorexactloglikelihood)

# The alternative variant of an RBM with Gaussian visible nodes and binary hidden
# nodes, as described by Cho et. al. in "Improved learning of Gaussian-Bernoulli
# restricted Boltzmann machines", is also available as "GaussianBernoulliRBM2".
rbm = fitrbm(x, rbmtype = GaussianBernoulliRBM2,
      nhidden = 3, epochs = 70, learningrate = 0.001);

# The same monitoring options as for BernoulliRBMs are available
# for GaussianBernoulliRBMs and GaussianBernoulliRBM2s.

# ==============================================================================
# Data with binomially distributed values in {0,1,2}: Binomial2BernoulliRBM
# ------------------------------------------------------------------------------

# Generate data and split it in training and test data set
nsamples = 100;
nvariables = 25;
Random.seed!(12);
x = barsandstripes(nsamples, nvariables) + barsandstripes(nsamples, nvariables);
x, xtest = splitdata(x, 0.3);
datadict = DataDict("Training data" => x, "Test data" => xtest);

# Fit Binomial2BernoulliRBM.
monitor = Monitor();
rbm = fitrbm(x, rbmtype = Binomial2BernoulliRBM,
      nhidden = 10, epochs = 100, learningrates = [0.003*ones(60); 0.001*ones(40)],
      monitoring = (rbm, epoch) ->
            monitorexactloglikelihood!(monitor, rbm, epoch, datadict));

BoltzmannMachinesPlots.plotevaluation(monitor, monitorexactloglikelihood)

# ==============================================================================
# Examples for conditional sampling
# ------------------------------------------------------------------------------
nsamples = 100;
nvariables = 36;

x = barsandstripes(100, 16);
dbm = fitdbm(x);

# Conditional gibbs sampling of particles
particles = initparticles(dbm, 500)
variablestofixto1 = 1:3
particles[1][:, variablestofixto1] .= 1.0
gibbssamplecond!(particles, dbm, variablestofixto1)
particles[1]

# Short alternative, returns only the visible nodes
samples(dbm, 500, conditions = [i => 1.0 for i in 1:3])

# ==============================================================================
# Examples for cross-validation
# ------------------------------------------------------------------------------

nsamples = 500;
nvariables = 16;
x = barsandstripes(nsamples, nvariables);

# Determine the optimal number of training epochs for a RBM
monitor = crossvalidation(x,
      (x, datadict) ->
            begin
               monitor = BoltzmannMachines.Monitor()
               BoltzmannMachines.fitrbm(x, epochs = 50,
                     learningrate = 0.005,
                     monitoring = (rbm, epoch) ->
                           BoltzmannMachines.monitorexactloglikelihood!(
                                 monitor, rbm, epoch, datadict))
               monitor
            end);
BoltzmannMachinesPlots.crossvalidationcurve(monitor)

# Determine the optimal number of pretraining epochs for a DBM
# given the other parameters
using Distributed
@everywhere using BoltzmannMachines
@everywhere function my_pretraining_monitoring(
      x::Matrix{Float64}, datadict::DataDict, epoch::Int)

   monitor = Monitor()
   monitorlogproblowerbound!(monitor,
         stackrbms(x,
               nhiddens = [16; 12],
               learningrate = 0.001,
               predbm = true,
               epochs = epoch),
         epoch,
         datadict,
         parallelized = false)
end

monitor = crossvalidation(x, my_pretraining_monitoring, 10:10:200);
BoltzmannMachinesPlots.crossvalidationcurve(monitor, monitorlogproblowerbound)

# ==============================================================================
# Examples for custom optimization algorithm
# ------------------------------------------------------------------------------

x = barsandstripes(100, 16);

struct MyRegularizedOptimizer{R<:AbstractRBM} <: AbstractOptimizer{R}
   llopt::LoglikelihoodOptimizer{R}
end

function MyRegularizedOptimizer()
   MyRegularizedOptimizer(loglikelihoodoptimizer(learningrate = 0.01))
end

import BoltzmannMachines: initialized, computegradient!, updateparameters!

function initialized(myoptimizer::MyRegularizedOptimizer,
      rbm::R) where {R<: AbstractRBM}
   MyRegularizedOptimizer(initialized(myoptimizer.llopt, rbm))
end

function computegradient!(myoptimizer::MyRegularizedOptimizer,
      v, vmodel, h, hmodel, rbm)
   computegradient!(myoptimizer.llopt, v, vmodel, h, hmodel, rbm)
   # add L2 regularization term
   myoptimizer.llopt.gradient.weights .-= 0.1 .* rbm.weights
   nothing
end

function updateparameters!(rbm::R, myoptimizer::MyRegularizedOptimizer{R}
      ) where {R<: AbstractRBM}
   updateparameters!(rbm, myoptimizer.llopt)
end

# The optimizer can be used for fitting RBMs ...
rbm = fitrbm(x; optimizer = MyRegularizedOptimizer())

# and also for DBMs
dbm = fitdbm(x; optimizer = MyRegularizedOptimizer())


text:
#!/usr/bin/env julia
import JSON
using LightGraphs, MetaGraphs

using Catlab
import Catlab.Diagram.Graphviz
using SemanticFlowGraphs

""" Visualize concepts as "concept map" or "mind map" using GraphViz.
"""
function concept_map()::MetaDiGraph
  db = OntologyDB()
  load_concepts(db)

  # Create graph for concept map.
  g = MetaDiGraph()
  set_props!(g, Dict(
    :graph => Dict(
      :fontname => "helvetica",
    ),
    :node => Dict(
      :fontname => "helvetica",
      :shape => "box",
      :style => "rounded",
      :width => "0",
      :height => "0",
    ),
    :edge => Dict(
      :fontname => "helvetica",
      :arrowsize => "0.75",
    ),
  ))

  # Add node for object generators.
  for ob in generators(concepts(db), Monocl.Ob)
    name = first(ob)
    doc = concept_document(db, name)
    add_vertex!(g, Dict(
      :id => name,
      :label => doc["name"],
    ))
  end
  set_indexing_prop!(g, :id)

  # Add edges for subobject generators.
  for subob in generators(concepts(db), Monocl.SubOb)
    dom_v, codom_v = g[first(dom(subob)),:id], g[first(codom(subob)),:id]
    add_multi_edge!(g, dom_v, codom_v, Dict(
      :xlabel => "is",
    ))
  end

  # Add edges for morphism generators with simple domain and codomain.
  subobjects_graph = subobjects(concepts(db))
  primitive_types = concepts(db, ["array", "scalar"])
  for hom in generators(concepts(db), Monocl.Hom)
    if !(head(dom(hom)) == :generator && head(codom(hom)) == :generator)
      continue
    end

    # Handle primitive types in codomain specially: create new node for each
    # occurrence, rather than linking back to original node.
    # This heuristic significantly improves the layout of the graph.
    dom_v, codom_v = g[first(dom(hom)),:id], g[first(codom(hom)),:id]
    if any(is_subobject(subobjects_graph, codom(hom), ob) for ob in primitive_types)
      name = first(codom(hom))
      doc = concept_document(db, name)
      add_vertex!(g, Dict(
        :label => doc["name"],
      ))
      codom_v = nv(g)
    end

    name = first(hom)
    doc = concept_document(db, name)
    add_multi_edge!(g, dom_v, codom_v, Dict(
      :id => name,
      :xlabel => doc["name"]
    ))
  end
  return g
end

function add_multi_edge!(g::MetaDiGraph, u::Int, v::Int, data::Dict=Dict())
  if has_edge(g, u, v)
    push!(get_prop(g, u, v, :edges), data)
  else
    add_edge!(g, u, v, :edges, [data])
  end
end

""" Compute all (explicit and implicit) subojects as directed graph.

TODO: This computer algebra belongs somewhere else, possibly in Catlab.
"""
function subobjects(concepts::Presentation)::MetaDiGraph
  # Encode (explicit) subobject relations as directed graph.
  g = MetaDiGraph()
  for ob in generators(concepts, Monocl.Ob)
    add_vertex!(g, :name, first(ob))
  end
  set_indexing_prop!(g, :name)
  for subob in generators(concepts, Monocl.SubOb)
    add_edge!(g, g[first(dom(subob)),:name], g[first(codom(subob)),:name])
  end

  # Compute reflexive, transitive closure to obtain all subobject relations.
  g_trans = MetaDiGraph(transitiveclosure(g.graph))
  for v in vertices(g)
    add_edge!(g_trans, v, v) # add self loops
    set_prop!(g_trans, v, :name, get_prop(g, v, :name))
  end
  set_indexing_prop!(g_trans, :name)
  return g_trans
end

function is_subobject(graph::MetaDiGraph,
                      ob1::Monocl.Ob{:generator}, ob2::Monocl.Ob{:generator})
  has_edge(graph, graph[first(ob1),:name], graph[first(ob2),:name])
end

function main(params::Dict)::Dict
  graphviz = Graphviz.to_graphviz(concept_map(); multigraph=true)
  return Dict(
    "graph" => sprint(Graphviz.pprint, graphviz),
    "prog" => "neato",
  )
end

if @__MODULE__() == Main
  const PARAMS = JSON.parse(isempty(ARGS) ? "{}" : ARGS[1])
  JSON.print(main(PARAMS))
end

max_stars_repo_name: sander/datascienceontology-backend

max_issues_repo_name: sander/datascienceontology-backend

max_forks_repo_name: sander/datascienceontology-backend


text:
"""set objective"""
function set_objective(sp, data::Dict)
    @stageobjective(sp, sum(values(sp.ext[:cost])))
end

"""add generators cost to cost dict"""
function add_gen_cost(sp, data::Dict)
    sp.ext[:cost][:gen_cost] = objective_function(sp)
end

"""add spillage cost to cost dict"""
function add_spill_cost(sp, data::Dict)
    sp.ext[:cost][:spill_cost] = sum(data["hydro"]["Hydrogenerators"][i]["spill_cost"]*sp[:spill][i] for i=1:data["hydro"]["nHyd"])
end

"""add deficit cost to cost dict"""
function add_deficit_cost(sp, data::Dict)
    sp.ext[:cost][:deficit_cost] = sum(sp[:deficit])*data["powersystem"]["baseMVA"]*data["powersystem"]["cost_deficit"]
end


text:
module MinLengthTest

using Hashids
using Test

const TESTDATA = [
    [0],
    [1, 2, 3],
    [10000000000, 0, 0, 0, 999999999999999],
    [827364827638428762424, 0x410498019283098083],
    [340282366920938463463374607431768211456]
]

function test_with_min_length(min_length::Int)
    config = Hashids.configure(min_length=min_length)
    for numbers in TESTDATA
        hashid = encode(config, numbers)
        @test length(hashid) ≥ min_length
        @test decode(config, hashid) == numbers
    end
end

@testset "Min Length" begin

@testset "0" begin
    test_with_min_length(0)
end  # of "0"

@testset "1" begin
    test_with_min_length(1)
end  # of "1"

@testset "10" begin
    test_with_min_length(10)
end  # of "10"

@testset "999" begin
    test_with_min_length(999)
end  # of "999"

@testset "1000" begin
    test_with_min_length(1000)
end  # of "1000"

end  # of "Min Length"

end  # module


text:
"""
---
A^{-1}_{22} = A^{22} - A^{21}(A^{11}^{-1})A^{12}
"""
function A_22inv(Ai, n1, nc)
    n2 = size(Ai, 1) - n1
    A22 = Ai[n1+1:end, n1+1:end]
    A12 = Ai[1:n1, n1+1:end]
    A11 = Ai[1:n1, 1:n1]

    v = A11 \ vec(A12[:, 1])
    A12'v
end


text:
# ---
# title: 151. Reverse Words in a String
# id: problem151
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/reverse-words-in-a-string/description/>
# hidden: true
# ---
#
# Given an input string `s`, reverse the order of the **words**.
#
# A **word** is defined as a sequence of non-space characters. The **words** in
# `s` will be separated by at least one space.
#
# Return _a string of the words in reverse order concatenated by a single
# space._
#
# **Note** that `s` may contain leading or trailing spaces or multiple spaces
# between two words. The returned string should only have a single space
# separating the words. Do not include any extra spaces.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "the sky is blue"
#     Output: "blue is sky the"
#
#
# **Example 2:**
#
#
#
#     Input: s = "  hello world  "
#     Output: "world hello"
#     Explanation: Your reversed string should not contain leading or trailing spaces.
#
#
# **Example 3:**
#
#
#
#     Input: s = "a good   example"
#     Output: "example good a"
#     Explanation: You need to reduce multiple spaces between two words to a single space in the reversed string.
#
#
# **Example 4:**
#
#
#
#     Input: s = "  Bob    Loves  Alice   "
#     Output: "Alice Loves Bob"
#
#
# **Example 5:**
#
#
#
#     Input: s = "Alice does not even like bob"
#     Output: "bob like even not does Alice"
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 104`
#   * `s` contains English letters (upper-case and lower-case), digits, and spaces `' '`.
#   * There is **at least one** word in `s`.
#
#
#
# **Follow up:**
#
#   * Could you solve it **in-place** with `O(1)` extra space?
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function handle_step(obj::FieldOrekitFixedStepHandler, arg0::FieldSpacecraftState, arg1::jboolean)
    return jcall(obj, "handleStep", void, (FieldSpacecraftState, jboolean), arg0, arg1)
end

function init(obj::FieldOrekitFixedStepHandler, arg0::FieldSpacecraftState, arg1::FieldAbsoluteDate, arg2::RealFieldElement)
    return jcall(obj, "init", void, (FieldSpacecraftState, FieldAbsoluteDate, RealFieldElement), arg0, arg1, arg2)
end


text:
# The format of simplegraph files is as follows:
# a one line header: <num_vertices>, <num_edges>, {"d" | "u"}, <name>
#   - num_vertices is an integer
#   - num_edges is an integer
#   - "d" for directed graph, "u" for undirected. Note that this
#       option does not perform any additional edge construction; it's
#       merely used to return the correct type of graph.
#   - name is a string
# header followed by a list of (comma-delimited) edges - src,dst.
# Multiple graphs may be present in one file.

function _lg_read_one_graph(f::IO, n_v::Integer, n_e::Integer, directed::Bool)
    if directed
        g = DiGraph(n_v)
    else
        g = Graph(n_v)
    end
    for i = 1:n_e
        line = chomp(readline(f))
        if length(line) > 0
            src_s, dst_s = split(line,r"\s*,\s*")
            src = parse(Int, src_s)
            dst = parse(Int, dst_s)
            add_edge!(g, src, dst)
        end
    end
    return g
end

function _lg_skip_one_graph(f::IO, n_e::Integer)
    for _ in 1:n_e
        readline(f)
    end
end

"""Returns a dictionary of (name=>graph) loaded from file `fn`."""
function loadlg_mult(io::IO)
    graphs = Dict{AbstractString, SimpleGraph}()
    while !eof(io)
        line = strip(chomp(readline(io)))
        if startswith(line,"#") || line == ""
            next
        else
            nvstr, nestr, dirundir, graphname = split(line, r"s*,s*", limit=4)
            n_v = parse(Int, nvstr)
            n_e = parse(Int, nestr)
            dirundir = strip(dirundir)
            graphname = strip(graphname)
            directed = !(dirundir == "u")

            g = _lg_read_one_graph(io, n_v, n_e, directed)
            graphs[graphname] = g
        end
    end
    return graphs
end

function loadlg(io::IO, gname::AbstractString)
    while !eof(io)
        line = strip(chomp(readline(io)))
        (startswith(line,"#") || line == "") && continue
        nvstr, nestr, dirundir, graphname = split(line, r"s*,s*", limit=4)
        n_v = parse(Int, nvstr)
        n_e = parse(Int, nestr)
        graphname = strip(graphname)
        if gname == graphname
            dirundir = strip(dirundir)
            directed = !(dirundir == "u")
            return _lg_read_one_graph(io, n_v, n_e, directed)
        else
            _lg_skip_one_graph(io, n_e)
        end
    end
    error("Graph $gname not found")
end

"""Writes a graph `g` with name `graphname` in a proprietary format
to the IO stream designated by `io`.

Returns 1 (number of graphs written).
"""
function savelg(io::IO, g::SimpleGraph, gname::AbstractString)
    # write header line
    dir = is_directed(g)? "d" : "u"
    line = join([nv(g), ne(g), dir, gname], ",")
    write(io, "$line\n")
    # write edges
    for e in edges(g)
        write(io, "$(src(e)),$(dst(e))\n")
    end
    return 1
end

"""Writes a dictionary of (name=>graph) to a file `fn`,
with default `GZip` compression.

Returns number of graphs written.
"""
function savelg_mult(io::IO, graphs::Dict)
    ng = 0
    for (gname, g) in graphs
        ng += savelg(io, g, gname)
    end
    return ng
end

# savelg(io::IO, g::SimpleGraph, n::AbstractString) =
#     savelg_mult(io, Dict(n=>g))

# write(g::Graph, fn::AbstractString; compress::Bool=true) = write(g, "graph", fn; compress=compress)
# write(g::DiGraph, fn::AbstractString; compress::Bool=true) = write(g, "digraph", fn; compress=compress)

filemap[:lg] = (loadlg, loadlg_mult, savelg, savelg_mult)


text:
function graphenefunction(z)
    f = 1e12
    c = 299792458.0
    μ₀ = 4π*1e-7
    ϵ₀ = 1/(μ₀*c^2)

    e = 1.602176565e-19
    kB = 1.3806488e-23
    hk = 1.05457168e-34
    vFe = 1e6
    muc = 0.05*e
    t = 0.135e-12
    T = 300
    ϵᵣ₁ = 1.0
    ϵᵣ₂ = 11.9

    ω = 2π*f
    k₀ = ω/c
    kᵣ₀ = -im*z*k₀

    Slo=-im*e^2*kB*T*log(2+2*cosh(muc/kB/T)) / (π*hk^2*(ω-im/t))

    a = -3*vFe^2*Slo/(4*(ω-im/t)^2)
    b = a/3

    Y1TM = ω*ϵᵣ₁*ϵ₀/sqrt(ϵᵣ₁*k₀^2 - kᵣ₀^2);
    Y2TM = ω*ϵᵣ₂*ϵ₀/sqrt(ϵᵣ₂*k₀^2 - kᵣ₀^2);
    YSTM = Slo + 1*a*kᵣ₀^2 + 1*b*kᵣ₀^2;

    # four Riemann sheets
    w = (Y1TM + Y2TM + YSTM)*(-Y1TM + Y2TM + YSTM)*(Y1TM - Y2TM + YSTM)*(-Y1TM - Y2TM + YSTM)
end

# Analysis parameters
xb = -100  # real part begin
xe = 400  # real part end
yb = -100  # imag part begin
ye = 400  # imag part end
r = 18  # initial mesh step

origcoords = rectangulardomain(complex(xb, yb), complex(xe, ye), r)

# matlab results from https://github.com/PioKow/GRPF for comparison
matlab_zroots = [-38.1777253145628 - 32.5295210454247im,
                 -32.1019622517269 - 27.4308619361753im,
                  32.1019622517269 + 27.4308619360714im,
                  38.17772531429 + 32.5295210455806im,
                  332.744888929695 + 282.243079954389im,
                  336.220287339074 + 285.191091013829im,
                  368.439467215558 + 312.522078059503im,
                  371.007570834263 + 314.700407676927im]

matlab_zpoles = [-2.30871731988513e-10 - 3.44963766202144im,
                 -2.65852297441317e-10 + 3.4496376622893im]

zroots, zpoles = grpf(graphenefunction, origcoords)

@test length(zroots) == 8
@test length(zpoles) == 2

@test approxmatch(zroots, matlab_zroots)
@test approxmatch(zpoles, matlab_zpoles)

pzroots, pzpoles, quadrants, phasediffs, tess, g2f = grpf(graphenefunction, origcoords, PlotData())

@test approxmatch(pzroots, matlab_zroots)
@test approxmatch(pzpoles, matlab_zpoles)

# Test with big origcoords
xb = big"-100"  # real part begin
xe = big"400"  # real part end
yb = big"-100"  # imag part begin
ye = big"400"  # imag part end
r = big"18"  # initial mesh step
tolerance = 1e-9

origcoords = rectangulardomain(complex(xb, yb), complex(xe, ye), r)

bzroots, bzpoles = grpf(graphenefunction, origcoords)

@test all(isa.(bzroots, Complex{BigFloat}))
@test all(isa.(bzpoles, Complex{BigFloat}))

@test approxmatch(bzroots, zroots)
@test approxmatch(bzpoles, zpoles)


text:
module LatentHeat

# If you want to add a new method here, feel free to do so.
# Remember to also export the function name in GeoParams.jl (in addition to here)

using Parameters, LaTeXStrings, Unitful
using ..Units
using GeoParams: AbstractMaterialParam
using ..MaterialParameters: MaterialParamsInfo
import Base.show, GeoParams.param_info

abstract type AbstractLatentHeat{T} <: AbstractMaterialParam end

export  compute_latent_heat,                  # calculation routines
        param_info,
        ConstantLatentHeat                  # constant

# Constant  -------------------------------------------------------
"""
    ConstantLatentHeat(Q_L=400kJ/kg)

Set a constant latent heat:
```math
    Q_L  = cst
```
where ``Q_L`` is the latent heat [``kJ/kg``].
"""
@with_kw_noshow struct ConstantLatentHeat{T,U} <: AbstractLatentHeat{T}
    Q_L::GeoUnit{T,U}         =   400kJ/kg                # Latent heat
end
ConstantLatentHeat(a...) = ConstantLatentHeat(convert.(GeoUnit,a)...)

function param_info(s::ConstantLatentHeat) # info about the struct
    return MaterialParamsInfo(Equation = L"Q_L = cst")
end

# Calculation routine
function compute_latent_heat(s::ConstantLatentHeat)
    @unpack_val Q_L   = s

    return Q_L
end

# Print info
function show(io::IO, g::ConstantLatentHeat)
    print(io, "Constant latent heat: Q_L=$(g.Q_L.val)")
end
#-------------------------------------------------------------------------

# Help info for the calculation routines
"""
    Ql = compute_latent_heat(s:<AbstractLatentHeat)

Returns the latent heat `Q_L`

"""
compute_latent_heat()

end


text:
using Documenter, MySample

makedocs(;
    modules=[MySample],
    format=Documenter.HTML(assets=String[]),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/terasakisatoshi/MySample.jl/blob/{commit}{path}#L{line}",
    sitename="MySample.jl",
    authors="SatoshiTerasaki <terasakisatoshi.math@gmail.com>",
)

deploydocs(;
    repo="github.com/terasakisatoshi/MySample.jl",
)


text:
var2(x) = mean(x.^2) .- mean(x)^2

"""
    waic(ll::AbstractArray{<:Real}; pointwise=false, log_lik="log_lik, kwargs...)

Compute the Widely Applicable Information Criterion (WAIC).

# Arguments
* `loglik::AbstractArray`   : A vector of posterior log likelihoods
* `pointwise::Bool`         : Compute WAIC pointwise, return a vector

# Returns
* `res::NamedTuple`: (WAIC=waics, lppd=lpd, penalty=pD, std_err=se) where

    WAIC                    : Sum of pointwise waic values (or pointwise vector)
    lppd                    : Log pointwise predictive density
    penalty                 : Penalty term ("overfitting penalty")
    std_err                 : Standard error of pointwise waic values

"""
function waic(ll::AbstractArray; pointwise=false)

    n_samples, n_obs = size(ll)
    pD = zeros(n_obs);

    lpd = reshape(logsumexp(ll .- log(n_samples); dims=1), n_obs);
    for i in 1:n_obs
        pD[i] = var2(ll[:,i])
    end

    waic_vec = (-2) .* ( lpd - pD );
    if pointwise == false
        waics = sum(waic_vec)
        lpd = sum(lpd)
        pD = sum(pD)
    else
        waics = waic_vec
    end

    local se
    try
        se = sqrt( n_obs*var2(waic_vec) )
    catch e
        println(e)
        se = nothing
    end

    (WAIC=waics, lppd=lpd, penalty=pD, std_err=se)
end


text:
using ArchGDAL
using StaticArrays

struct Geo
    origin::SVector{2,Float64}
    transform::SMatrix{2, 2, Float64}
end

"""
Find the origin and linear transform from pixel to x-y.
"""
function band_geo(band::ArchGDAL.AbstractRasterBand)
    geo = ArchGDAL.getgeotransform(band)
    Geo([geo[1], geo[4]], [geo[2] geo[3]; geo[5] geo[6]])
end

pixel_corner(geo::Geo, ij) = geo.transform * (ij .- 1) + geo.origin
pixel_center(geo::Geo, ij) = geo.transform * (ij .- 0.5) + geo.origin
# First column is upper-left (x y). Second column is lower-right (x, y).
raster_rect(geo::Geo, ullr) = pixel_corner(geo, hcat(ullr[:, 1], ullr[:, 2] .+ 1))
partial_pixel_containing(geo::Geo, xy) = (geo.transform \ (xy - origin)) .+ 1
offset_geo(geo::Geo, offset) = Geo(pixel_corner(geo, offset), geo.transform)

function ij_cover_rect(geo::Geo, rect; eps = 1e-7)
    ul = partial_pixel_containing(geo, rect[:, 1])
    ceil_dx = ceil.[ul] - ul
    ulij = Int32.(floor.(ul))
    for uli in 1:2
        if ceil_dx[uli] < eps
            ulij[uli] = ceil_dx[uli]
        end
    end
    lr = partial_pixel_containing(geo, rect[:, 2])
    floor_dx = floor.[lr] - lr
    lrij = Int32.(ceil.(lr))
    for lri in 1:2
        if floor_dx[lri] < eps
            lrij[lri] = floor_dx[lri]
        end
    end
    hcat(ulij, lrij)
end

function block_cover_rect(ul_lr, bounds, blocksize)
    # available bounds may be less than what covers this rectangle.
    ula = [  # first column is upper-left. second column is lower-right.
        max(ul_lr[1, 1], bounds[1][1])  min(ul_lr[1, 2], bounds[1][end])
        max(ul_lr[2, 1], bounds[2][1])  min(ul_lr[2, 2], bounds[2][end])
    ]
    # return as two UnitRanges, one for each axis.
    toblock(ij, blocksize) = (ij .- 1) ÷ blocksize .+ 1
    ul = toblock(ula[:, 1])
    lr = toblock(ula[:, 2])
    (ul[1]:lr[1], ul[2]:lr[2])
end

struct RasterBlock{T}
    A::Array{T}
    range::Tuple{UnitRange{Int64},UnitRange{Int64}}
    geo::Geo
    nodatavalue::T
end

"""
The processing function works one block at a time and returns whether
    it wrote to a block.
"""
function process_block(blocka::RasterBlock, blockb::RasterBlock)
    for pixela in iterblock(blocka)
        rect = pixel_rect(blocka, pixela)
        cnt = 0
        for pixel in inside_rect(blockb, rect)
            if pixel != blockb.nodatavalue
                cnt += 1
            end
        end
        rate = min(blocka.A[pixela] / cnt, 100.0)
        for pixel in inside_rect(blockb, rect)
            if pixel != blockb.nodatavalue
                blockb[pixel] = rate
            end
        end
    end
    false, true
end

function over_block_pairs(processor::Function, band, geo)
    blocksize = ArchGDAL.blocksize.(band)
    bounds = (bd -> [1:ArchGDAL.width(bd), 1:ArchGDAL.height(bd)]).(band)
    dtype = ArchGDAL.pixeltype.(band)
    # Separate these because they are different types.
    nodata1 = ArchGDAL.getnodatavalue(band1)
    nodata2 = ArchGDAL.getnodatavalue(band2)

    rect2 = raster_rect(geo[2], bounds[2])
    ul_lr1 = ij_cover_rect(geo[1], rect2)
    block_bounds1 = block_cover_rect(ul_lr1, bounds[1], blocksize[1])
    buffer1 = zeros(dtype[1], bounds[1])
    for j in block_bounds1[2]
        for i in block_bounds1[1]
            fill!(buffer1, nodata[1])
            readblock!(band[1], i - 1, j - 1, buffer1)
            block_geo1 = offset_geo(geo[1], blocksize[1] .* ([i, j] .- 1))

            block1_rect = raster_rect(block_geo1, axes(buffer1))
            ul_lr2 = ij_cover_rect(geo[2], block1_rect)
            block_bounds2 = block_cover_rect(ul_lr2, bounds[2], blocksize[2])
            xdim = block_bounds2[1][end] - block_bounds2[1][1]
            ydim = block_bounds2[2][end] - block_bounds2[2][1]
            buffer2 = zeros(dtype[2], xdim, ydim)
            smallbuffer = zeros(dtype[2], blocksize[2]...)
            for j2 in block_bounds2[2]
                for i2 in block_bounds2[1]
                    readblock!(band[2], i2 - 1, j2 - 1, smallbuffer)
                    buffer2[
                        ((i2 - 1) * blocksize[2][1] + 1):(i2 * blocksize[2][1]),
                        ((j2 - 1) * blocksize[2][2] + 1):(j2 * blocksize[2][2]),
                        ] = smallbuffer
                end
            end
            # These are missing the pixel ranges that are allowed.
            rb1 = RasterBlock{dtype[1]}(buffer1, block_geo1)
            rb2 = RasterBlock{dtype[2]}(buffer2, block_geo2)
            towrite = processor(rb1, rb2)
        end
    end
end

function adjust_hrsl_with_landscan(a::ArchGDAL.Dataset, b::ArchGDAL.Dataset)
    band = [ArchGDAL.getband(a, 1), ArchGDAL.getband(b, 1)]
    geo = band_geo.([a, b])
    over_block_pairs(process_block, band, geo)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -1.3634297333714107 <= q <= -0.20981772159177886)

                     add_NL_constraint(m, :(swish(-0.25086622252953283 + -0.4158912166699471*swish(0.18862560740304835 + -0.15744636994779526*swish(0.5575683410422836 + 0.3104874388021095*swish(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])) + 0.12330243724213741*swish(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))) + 0.6894873557418655*swish(0.5818927234969289 + 0.4060546195732435*swish(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])) + 0.22410021012329473*swish(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4])))) + 0.15556574286983027*swish(-0.6464546506721769 + -0.3052437329487665*swish(0.5575683410422836 + 0.3104874388021095*swish(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])) + 0.12330243724213741*swish(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))) + 0.4287563119683666*swish(0.5818927234969289 + 0.4060546195732435*swish(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])) + 0.22410021012329473*swish(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))) + swish(-0.9805026843034819 + -0.19393497198753096*swish(0.18862560740304835 + -0.15744636994779526*swish(0.5575683410422836 + 0.3104874388021095*swish(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])) + 0.12330243724213741*swish(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))) + 0.6894873557418655*swish(0.5818927234969289 + 0.4060546195732435*swish(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])) + 0.22410021012329473*swish(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4])))) + -0.01642706539433636*swish(-0.6464546506721769 + -0.3052437329487665*swish(0.5575683410422836 + 0.3104874388021095*swish(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])) + 0.12330243724213741*swish(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))) + 0.4287563119683666*swish(0.5818927234969289 + 0.4060546195732435*swish(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])) + 0.22410021012329473*swish(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Revise
using Distributed

N_PROCS=56
addprocs(N_PROCS)
@everywhere begin
    using POMDPs
    using GridInterpolations
    using POMDPModelTools
    using POMDPPolicies

    using Parameters
    using StaticArrays
    using DiscreteValueIteration

    using AutomotiveDrivingModels
    using AutoViz
    using AutomotiveSensors
    using AutomotivePOMDPs
    using PedestrianAvoidancePOMDP

    pomdp = SingleOCFPOMDP()
end

solver = ParallelValueIterationSolver(n_procs=N_PROCS, max_iterations=200, belres=1e-4, include_Q=true, verbose=true)

vi_policy = solve(solver, pomdp)
qmdp_policy = AlphaVectorPolicy(pomdp, vi_policy.qmat, vi_policy.action_map)

# save policy!
FileIO.save("policy.jld2", "policy", qmdp_policy)


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

JLLWrappers.@generate_main_file_header("Cairomm")
JLLWrappers.@generate_main_file("Cairomm", UUID("af74c99f-f0eb-54aa-aecc-a10e8fc65c17"))
end  # module Cairomm_jll


text:
# Kernels for unary array operations

fp = open("cuda1.cu","w")
#using Knet: unary_ops

include("gamma.jl")
print(fp,cuda1gammafamily())

function cuda1src(f, j=f, ex="$f(xi)"; seperate_impl=false, BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","$(f)_32"),("double","$(f)_64")]
            seperate_impl && (ex = "$F(xi)")
            print(s,
"""
__global__ void _$F(int n, $T *x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    $T xi = x[i];
    y[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  $DLLEXPORT void $F(int n, $T *x, $T *y) {
    if (n>0) _$F<<<$BLK,$THR>>>(n,x,y);
  }
}
""")
        end
    end
end

seperate_impl_ops = ["gamma_impl", "digamma_impl", "trigamma_impl"]
for a in unary_ops
    if !isa(a,Tuple); a=(a,); end
    print(fp,cuda1src(a...; seperate_impl=(a[1] in seperate_impl_ops)))
end

# Kernels used by setindex! and getindex: fill, xfill, xcopy:

function cuda1fill(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _fill_$F(int n, $T x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    y[i] = x;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  $DLLEXPORT void fill_$F(int n, $T x, $T *y) {
    if (n>0) _fill_$F<<<$BLK,$THR>>>(n,x,y);
  }
}
""")
        end
    end
end

print(fp,cuda1fill())

function cuda1xfill(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _xfill_$F(int nrows, int ncols, $T x, $T *y, int incy) {
  int row, col, yidx;
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = i % nrows;
    col = i / nrows;
    if (col >= ncols) break;
    yidx = row + col * incy;
    y[yidx] = x;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  $DLLEXPORT void xfill_$F(int nrows, int ncols, $T x, $T *y, int incy) {
    if (nrows>0 && ncols>0) _xfill_$F<<<$BLK,$THR>>>(nrows, ncols, x, y, incy);
  }
}
""")
        end
    end
end

print(fp,cuda1xfill())

function cuda1xcopy(; BLK=256, THR=256)
"""
__global__ void _xcopy(int nrows, int ncols, const char *x, int incx, char *y, int incy) {
  int row, col, xidx, yidx;
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = i % nrows;
    col = i / nrows;
    if (col >= ncols) break;
    xidx = row + col * incx;
    yidx = row + col * incy;
    y[yidx] = x[xidx];
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  $DLLEXPORT void xcopy(int nrows, int ncols, const void *x, int incx, void *y, int incy) {
    if (nrows>0 && ncols>0) _xcopy<<<$BLK,$THR>>>(nrows,ncols,(char*)x,incx,(char*)y,incy);
  }
}
"""
end

print(fp,cuda1xcopy())

function cuda1icat(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _icat_$F(int nrows, int ncols, $T **x, $T *y) {
  int row, col, yidx;
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = i % nrows;
    col = i / nrows;
    if (col >= ncols) break;
    yidx = row + col * nrows;
    y[yidx] = x[col][row];
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  $DLLEXPORT void icat_$F(int nrows, int ncols, $T **x, $T *y) {
    $T **xx;
    if (nrows>0 && ncols>0) {
      size_t s = ncols * sizeof($T *);
      cudaMalloc(&xx, s);
      cudaMemcpy(xx, x, s, cudaMemcpyHostToDevice);
      _icat_$F<<<$BLK,$THR>>>(nrows, ncols, xx, y);
      cudaFree(xx);
    }
  }
}
""")
        end
    end
end

print(fp,cuda1icat())

# This is for missing double atomicAdd()
print(fp,"""
static __inline__ __device__ float atomicAdd2(float *address, float val) {
  return atomicAdd(address, val);
}
#if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 600
static __inline__ __device__ double atomicAdd2(double *address, double val) {
  return atomicAdd(address, val);
}
#else
static __inline__ __device__ double atomicAdd2(double *address, double val) {
  unsigned long long int* address_as_ull = (unsigned long long int*)address;
  unsigned long long int old = *address_as_ull, assumed;
  if (val==0.0)
    return __longlong_as_double(old);
  do {
    assumed = old;
    old = atomicCAS(address_as_ull, assumed, __double_as_longlong(val +__longlong_as_double(assumed)));
  } while (assumed != old);
  return __longlong_as_double(old);
}
#endif
""")

function cuda1getcols(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _getcols_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T *y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % xrows;
    col = yidx / xrows;
    if (col >= ncols) break;
    xidx = row + (cols[col]-1) * xrows;
    y[yidx] = x[xidx];
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _setcols_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T *y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % xrows;
    col = yidx / xrows;
    if (col >= ncols) break;
    xidx = row + (cols[col]-1) * xrows;
    x[xidx] = y[yidx];
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _addcols_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T *y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % xrows;
    col = yidx / xrows;
    if (col >= ncols) break;
    xidx = row + (cols[col]-1) * xrows;
    atomicAdd2(&x[xidx], y[yidx]);
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _setcol1_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % xrows;
    col = yidx / xrows;
    if (col >= ncols) break;
    xidx = row + (cols[col]-1) * xrows;
    x[xidx] = y;
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _getrows_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T *y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % nrows;
    col = yidx / nrows;
    if (col >= xcols) break;
    xidx = rows[row] - 1 + col * xrows;
    y[yidx] = x[xidx];
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _setrows_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T *y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % nrows;
    col = yidx / nrows;
    if (col >= xcols) break;
    xidx = rows[row] - 1 + col * xrows;
    x[xidx] = y[yidx];
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _addrows_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T *y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % nrows;
    col = yidx / nrows;
    if (col >= xcols) break;
    xidx = rows[row] - 1 + col * xrows;
    atomicAdd2(&x[xidx], y[yidx]);
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _setrow1_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T y) {
  int row, col, xidx;
  int yidx = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = yidx % nrows;
    col = yidx / nrows;
    if (col >= xcols) break;
    xidx = rows[row] - 1 + col * xrows;
    x[xidx] = y;
    yidx += blockDim.x * gridDim.x;
  }
}
__global__ void _getents_$F(int n, int *ents, $T *x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    y[i] = x[ents[i]-1];
    i += blockDim.x * gridDim.x;
  }
}
__global__ void _setents_$F(int n, int *ents, $T *x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    x[ents[i]-1] = y[i];
    i += blockDim.x * gridDim.x;
  }
}
__global__ void _addents_$F(int n, int *ents, $T *x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    atomicAdd2(&x[ents[i]-1], y[i]);
    i += blockDim.x * gridDim.x;
  }
}
__global__ void _setent1_$F(int n, int *ents, $T *x, $T y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    x[ents[i]-1] = y;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
$DLLEXPORT void getcols_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T *y)
{ if (ncols>0 && xrows>0 && xcols>0) _getcols_$F<<<$BLK,$THR>>>(xrows,xcols,ncols,cols,x,y); }
$DLLEXPORT void setcols_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T *y)
{ if (ncols>0 && xrows>0 && xcols>0) _setcols_$F<<<$BLK,$THR>>>(xrows,xcols,ncols,cols,x,y); }
$DLLEXPORT void addcols_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T *y)
{ if (ncols>0 && xrows>0 && xcols>0) _addcols_$F<<<$BLK,$THR>>>(xrows,xcols,ncols,cols,x,y); }
$DLLEXPORT void setcol1_$F(int xrows, int xcols, int ncols, int *cols, $T *x, $T  y)
{ if (ncols>0 && xrows>0 && xcols>0) _setcol1_$F<<<$BLK,$THR>>>(xrows,xcols,ncols,cols,x,y); }
$DLLEXPORT void getrows_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T *y)
{ if (nrows>0 && xrows>0 && xcols>0) _getrows_$F<<<$BLK,$THR>>>(xrows,xcols,nrows,rows,x,y); }
$DLLEXPORT void setrows_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T *y)
{ if (nrows>0 && xrows>0 && xcols>0) _setrows_$F<<<$BLK,$THR>>>(xrows,xcols,nrows,rows,x,y); }
$DLLEXPORT void addrows_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T *y)
{ if (nrows>0 && xrows>0 && xcols>0) _addrows_$F<<<$BLK,$THR>>>(xrows,xcols,nrows,rows,x,y); }
$DLLEXPORT void setrow1_$F(int xrows, int xcols, int nrows, int *rows, $T *x, $T  y)
{ if (nrows>0 && xrows>0 && xcols>0) _setrow1_$F<<<$BLK,$THR>>>(xrows,xcols,nrows,rows,x,y); }
$DLLEXPORT void getents_$F(int n, int *ents, $T *x, $T *y)
{ if (n>0) _getents_$F<<<$BLK,$THR>>>(n,ents,x,y); }
$DLLEXPORT void setents_$F(int n, int *ents, $T *x, $T *y)
{ if (n>0) _setents_$F<<<$BLK,$THR>>>(n,ents,x,y); }
$DLLEXPORT void addents_$F(int n, int *ents, $T *x, $T *y)
{ if (n>0) _addents_$F<<<$BLK,$THR>>>(n,ents,x,y); }
$DLLEXPORT void setent1_$F(int n, int *ents, $T *x, $T  y)
{ if (n>0) _setent1_$F<<<$BLK,$THR>>>(n,ents,x,y); }
}
""")
        end
    end
end

print(fp,cuda1getcols())

# Dropout

function cuda1dropout(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _dropout_$F(int n, $T p, $T q, $T *x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    if (y[i] > p) {
      y[i] = x[i] * q;
    } else {
      y[i] = 0;
    }
    i += blockDim.x * gridDim.x;
  }
}
__global__ void _dropback_$F(int n, $T q, $T *y, $T *dy, $T *dx) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    if (y[i] == 0) {
        dx[i] = 0;
    } else {
        dx[i] = dy[i] * q;
    }
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  $DLLEXPORT void dropout_$F(int n, $T p, $T *x, $T *y) {
    if (n>0) _dropout_$F<<<$BLK,$THR>>>(n,p,1.0/(1.0-p),x,y);
  }
  $DLLEXPORT void dropback_$F(int n, $T p, $T *x, $T *y, $T *dy, $T *dx) {
    if (n>0) _dropback_$F<<<$BLK,$THR>>>(n,1.0/(1.0-p),y,dy,dx);
  }
}
""")
        end
    end
end

print(fp,cuda1dropout())

# This is still too slow compared to concat on cpu and copy to gpu
# Tested for 25 arrays of 200
function cuda1concat(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _concat_$F(int narrays, int *starts, int *lengths, $T **x, $T *y) {
  int array = blockIdx.x;
  int nelts = lengths[array];
  int offset = starts[array];
  for (int i = threadIdx.x; i < nelts; i += blockDim.x) {
    y[i+offset] = x[array][i];
  }
}
extern "C" {
  // julia is responsible for copying args to gpu
  $DLLEXPORT void concat_$F(int narrays, int *starts, int *lengths, $T **x, $T *y) {
    _concat_$F<<<narrays,$THR>>>(narrays, starts, lengths, x, y);
  }
}
""")
        end
    end
end

print(fp,cuda1concat())

close(fp)

# Here is the test script for cuda1concat:
# using Knet, BenchmarkTools
# using Knet: @knet8

# for S in (32,64); T = Symbol("Float$S"); F = "concat_$S"
# @eval function concat(A::KnetArray{$T}...)
#     nargs = length(A)
#     S = Array{Int32}(nargs)
#     L = Array{Int32}(nargs)
#     nelts = 0
#     @inbounds for i in 1:nargs
#         n = length(A[i])
#         S[i] = nelts
#         L[i] = n
#         nelts += n
#     end
#     S = KnetArray(S)
#     L = KnetArray(L)
#     X = KnetArray([map(pointer,A)...])
#     Y = KnetArray{$T}(nelts)
#     @knet8($F,(Cint,Ptr{Cint},Ptr{Cint},Ptr{Ptr{$T}},Ptr{$T}),nargs,S,L,X,Y)
#     return Y
# end
# end

# a = [ rand(Float32,200) for i=1:25 ]
# k = map(KnetArray,a)
# @show vcat(a...) == vcat(k...)
# @show vcat(a...) == concat(k...)
# @show @benchmark vcat(a...)
# @show @benchmark vcat(k...)
# @show @benchmark concat(k...)
# @show @benchmark KnetArray(vcat(a...))


text:
"""
perturbation.jl

General functions for perturbation
"""

function pairpotential_minimum(pert::Perturbation, u::Matrix{Pairpotential}, rmin::Matrix{Float64})::Matrix{Float64}
  N::Int = ncomp(pert)

  return [ u[i,j](rmin[i,j]) for i in 1:N, j in 1:N ]
end

function pairpotential_minimum(pert::Perturbation)
  u::Matrix{Pairpotential} = pairpotential(pert)
  rmin::Matrix{Float64} = pairpotential_minimizer(pert)

  return pairpotential_minimum(pert, u, rmin)
end

function hsdiameter_estimate(pert::Perturbation, u::Matrix{Pairpotential}, rmin::Matrix{Float64}, T::Float64)::Matrix{Float64}
  N::Int = ncomp(pert)
  rcore::Matrix{Float64} = coreradius(pert)

  ret = Matrix{Float64}(N,N)
  for i in 1:N, j in 1:N
    i > j && continue
    fopt(σ::Float64)::Float64 = abs(u[i,j](σ) - u[i,j](rmin[i,j]) - kB*T)
    opt = Optim.optimize(fopt, rcore[i,j], rmin[i,j])
    ret[i,j] = ret[j,i] = Optim.minimizer(opt)
  end

  return ret
end

function hsdiameter_estimate(pert::Perturbation, T::Float64)::Matrix{Float64}
  u::Matrix{Pairpotential} = pairpotential(pert)
  rmin::Matrix{Float64} = pairpotential_minimizer(pert)

  return hsdiameter_estimate(pert, u, rmin, T)
end


text:
# GPU memory management and pooling

using Printf
using TimerOutputs

## allocation statistics

mutable struct AllocStats
  # pool allocation requests
  pool_nalloc::Int
  pool_nfree::Int
  ## in bytes
  pool_alloc::Int

  # actual CUDA allocations
  actual_nalloc::Int
  actual_nfree::Int
  ## in bytes
  actual_alloc::Int
  actual_free::Int

  pool_time::Float64
  actual_time::Float64
end

const alloc_stats = AllocStats(0, 0, 0, 0, 0, 0, 0, 0, 0)

Base.copy(alloc_stats::AllocStats) =
  AllocStats((getfield(alloc_stats, field) for field in fieldnames(AllocStats))...)

## CUDA allocator

const alloc_to = TimerOutput()

alloc_timings() = (show(alloc_to; allocations=false, sortby=:name); println())

const usage = Ref(0)
const usage_limit = Ref{Union{Nothing,Int}}(nothing)

const allocated = Dict{CuPtr{Nothing},Mem.DeviceBuffer}()

function actual_alloc(bytes)::Union{Nothing,CuPtr{Nothing}}
  # check the memory allocation limit
  if usage_limit[] !== nothing
    if usage[] + bytes > usage_limit[]
      return nothing
    end
  end

  # try the actual allocation
  try
    alloc_stats.actual_time += Base.@elapsed begin
      @timeit alloc_to "alloc" buf = Mem.alloc(Mem.Device, bytes)
    end
    @assert sizeof(buf) == bytes
    alloc_stats.actual_nalloc += 1
    alloc_stats.actual_alloc += bytes
    usage[] += bytes
    ptr = convert(CuPtr{Nothing}, buf)
    @assert !haskey(allocated, ptr)
    allocated[ptr] = buf
    return ptr
  catch ex
    ex == CUDAdrv.ERROR_OUT_OF_MEMORY || rethrow()
  end

  return nothing
end

function actual_free(ptr::CuPtr{Nothing})
  buf = allocated[ptr]
  delete!(allocated, ptr)
  bytes = sizeof(buf)

  alloc_stats.actual_nfree += 1
  alloc_stats.actual_free += bytes
  usage[] -= bytes

  @timeit alloc_to "free"  begin
    alloc_stats.actual_time += Base.@elapsed Mem.free(buf)
  end

  return
end

## memory pool

const pool_to = TimerOutput()

macro pool_timeit(args...)
    TimerOutputs.timer_expr(__module__, false, :($CuArrays.pool_to), args...)
end

pool_timings() = (show(pool_to; allocations=false, sortby=:name); println())

# API:
# - init()
# - deinit()
# - alloc(sz)::CuPtr{Nothing}
# - free(::CuPtr{Nothing})
# - used_memory()
# - cached_memory()

include("memory/binned.jl")
include("memory/simple.jl")
include("memory/split.jl")
include("memory/dummy.jl")

const pool = Ref{Union{Nothing,Module}}(nothing)

const requested = Dict{CuPtr{Nothing},Int}()

@inline function alloc(sz)::CuPtr{Nothing}
  # 0-byte allocations shouldn't hit the pool
  sz == 0 && return CU_NULL

  @assert pool[] !== nothing "Cannot allocate before CuArrays has been initialized."
  alloc_stats.pool_time += Base.@elapsed begin
    @pool_timeit "pooled alloc" ptr = pool[].alloc(sz)
  end
  if ptr === nothing
    @error "Out of GPU memory trying to allocate $(Base.format_bytes(sz))"
    pool[].dump()
    memory_status()
    throw(OutOfMemoryError())
  end

  alloc_stats.pool_nalloc += 1
  alloc_stats.pool_alloc += sz
  @assert !haskey(requested, ptr)
  requested[ptr] = sz

  return ptr
end

@inline function free(ptr::CuPtr{Nothing})
  # 0-byte allocations shouldn't hit the pool
  ptr == CU_NULL && return

  @assert haskey(requested, ptr)
  sz = requested[ptr]
  delete!(requested, ptr)

  alloc_stats.pool_nfree += 1
  alloc_stats.pool_time += Base.@elapsed begin
    @pool_timeit "pooled free" pool[].free(ptr)
  end

  return
end

function memory_pool!(mod::Module=BinnedPool)
  if pool[] !== nothing
    pool[].deinit()
  end

  reset_timers!()

  pool[] = mod
  mod.init()

  return
end

pool_dump() = pool[].dump()

## utilities

macro allocated(ex)
    quote
        let
            local f
            function f()
                b0 = alloc_stats.pool_alloc
                $(esc(ex))
                alloc_stats.pool_alloc - b0
            end
            f()
        end
    end
end

"""
    @time ex

Run expression `ex` and report on execution time and GPU/CPU memory behavior. The GPU is
synchronized right before and after executing `ex` to exclude any external effects.

"""
macro time(ex)
    quote
        # @time might surround an application, so be sure to initialize CUDA before that
        # FIXME: this should be done in CUDAdrv (`synchronize(ctx=CuCurrentOrNewContext()`)
        #        but the CUDA initialization mechanics are part of CUDAnative.jl
        CUDAnative.maybe_initialize("@time")

        # coarse synchronization to exclude effects from previously-executed code
        CUDAdrv.synchronize()

        local gpu_mem_stats0 = copy(alloc_stats)
        local cpu_mem_stats0 = Base.gc_num()
        local cpu_time0 = time_ns()

        # fine-grained synchronization of the code under analysis
        local val = @sync $(esc(ex))

        local cpu_time1 = time_ns()
        local cpu_mem_stats1 = Base.gc_num()
        local gpu_mem_stats1 = copy(alloc_stats)

        local cpu_time = (cpu_time1 - cpu_time0) / 1e9
        local gpu_gc_time = gpu_mem_stats1.pool_time - gpu_mem_stats0.pool_time
        local gpu_alloc_count = gpu_mem_stats1.pool_nalloc - gpu_mem_stats0.pool_nalloc
        local gpu_lib_time = gpu_mem_stats1.actual_time - gpu_mem_stats0.actual_time
        local gpu_alloc_size = gpu_mem_stats1.pool_alloc - gpu_mem_stats0.pool_alloc
        local cpu_mem_stats = Base.GC_Diff(cpu_mem_stats1, cpu_mem_stats0)
        local cpu_gc_time = cpu_mem_stats.total_time / 1e9
        local cpu_alloc_count = Base.gc_alloc_count(cpu_mem_stats)
        local cpu_alloc_size = cpu_mem_stats.allocd

        Printf.@printf("%10.6f seconds", cpu_time)
        for (typ, gctime, libtime, bytes, allocs) in
            (("CPU", cpu_gc_time, 0, cpu_alloc_size, cpu_alloc_count),
             ("GPU", gpu_gc_time, gpu_lib_time, gpu_alloc_size, gpu_alloc_count))
          if bytes != 0 || allocs != 0
              allocs, ma = Base.prettyprint_getunits(allocs, length(Base._cnt_units), Int64(1000))
              if ma == 1
                  Printf.@printf(" (%d%s %s allocation%s: ", allocs, Base._cnt_units[ma], typ, allocs==1 ? "" : "s")
              else
                  Printf.@printf(" (%.2f%s %s allocations: ", allocs, Base._cnt_units[ma], typ)
              end
              print(Base.format_bytes(bytes))
              if gctime > 0
                  Printf.@printf(", %.2f%% gc time", 100*gctime/cpu_time)
                if libtime > 0
                    Printf.@printf(" of which %.2f%% spent allocating", 100*libtime/gctime)
                end
              end
              print(")")
          elseif gctime > 0
              Printf.@printf(", %.2f%% %s gc time", 100*gctime/cpu_time, typ)
          end
        end
        println()

        val
    end
end

function memory_status()
  free_bytes, total_bytes = CUDAdrv.Mem.info()
  used_bytes = total_bytes - free_bytes
  used_ratio = used_bytes / total_bytes

  @printf("Effective GPU memory usage: %.2f%% (%s/%s)\n",
          100*used_ratio, Base.format_bytes(used_bytes),
          Base.format_bytes(total_bytes))

  @printf("CuArrays GPU memory usage: %s", Base.format_bytes(usage[]))
  if usage_limit[] !== nothing
    @printf(" (capped at %s)", Base.format_bytes(usage_limit[]))
  end
  println()

  alloc_used_bytes = pool[].used_memory()
  alloc_cached_bytes = pool[].cached_memory()
  alloc_total_bytes = alloc_used_bytes + alloc_cached_bytes

  @printf("%s usage: %s (%s allocated, %s cached)\n", nameof(pool[]),
          Base.format_bytes(alloc_total_bytes), Base.format_bytes(alloc_used_bytes),
          Base.format_bytes(alloc_cached_bytes))

  requested_bytes = reduce(+, values(requested); init=0)

  @printf("%s efficiency: %.2f%% (%s requested, %s allocated)\n", nameof(pool[]),
          100*requested_bytes/usage[],
          Base.format_bytes(requested_bytes),
          Base.format_bytes(usage[]))

  # check if the memory usage as counted by the CUDA allocator wrapper
  # matches what is reported by the pool implementation
  discrepancy = usage[] - alloc_total_bytes
  if discrepancy != 0
    @debug "Discrepancy of $(Base.format_bytes(discrepancy)) between memory pool and allocator"
  end
end

## init

function __init_memory__()
  if haskey(ENV, "CUARRAYS_MEMORY_LIMIT")
    usage_limit[] = parse(Int, ENV["CUARRAYS_MEMORY_LIMIT"])
  end

  if haskey(ENV, "CUARRAYS_MEMORY_POOL")
    memory_pool!(
      if ENV["CUARRAYS_MEMORY_POOL"] == "binned"
        BinnedPool
      elseif ENV["CUARRAYS_MEMORY_POOL"] == "simple"
        SimplePool
      elseif ENV["CUARRAYS_MEMORY_POOL"] == "split"
        SplittingPool
      elseif ENV["CUARRAYS_MEMORY_POOL"] == "none"
        DummyPool
      else
        error("Invalid allocator selected")
      end)
  else
    memory_pool!()
  end

  # if the user hand-picked an allocator, be a little verbose
  if haskey(ENV, "CUARRAYS_MEMORY_POOL")
    atexit(()->begin
      Core.println("""
        CuArrays.jl $(nameof(pool[])) statistics:
         - $(alloc_stats.pool_nalloc) pool allocations: $(Base.format_bytes(alloc_stats.pool_alloc)) in $(round(alloc_stats.pool_time; digits=2))s
         - $(alloc_stats.actual_nalloc) CUDA allocations: $(Base.format_bytes(alloc_stats.actual_alloc)) in $(round(alloc_stats.actual_time; digits=2))s""")
    end)
  end
end

function reset_timers!()
  TimerOutputs.reset_timer!(alloc_to)
  TimerOutputs.reset_timer!(pool_to)
end


text:
# * Setup
module plotting
using PyPlot
using PyCall
using Compat

const COL = PyNULL()
const masked_array = PyNULL()
const pgf_backend = PyNULL()
const pyplot_collections = PyNULL()
const RegularPolyCollection = PyNULL()
rcParams = 0

function __init__()
    copy!(COL, pyimport_conda("matplotlib.colors", "matplotlib"))
    copy!(masked_array, pyimport_conda("numpy.ma", "numpy"))
    copy!(pgf_backend, pyimport_conda("matplotlib.backends.backend_pgf", "matplotlib"))
    copy!(pyplot_collections, pyimport_conda("matplotlib.collections", "matplotlib"))
    copy!(RegularPolyCollection, pyplot_collections.RegularPolyCollection)
    global rcParams
    rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
end

using UnicodeFun

import Jagot: meshgrid

# using GSL
using Printf

plot_style(style::String) = matplotlib.style.use(style)

# * Figure wrappers

function ax_common(;nox=false, noy=false)
    if nox
        no_tick_labels(:x)
        xlabel("")
    end
    if noy
        no_tick_labels(:y)
        ylabel("")
    end
end

function toggle_toolbar(visible)
    rcParams["toolbar"] = visible ? "toolbar2" : "None"
end

toggle_statusbar(visible, fig=gcf()) = hasproperty(fig.canvas, :window) && fig.canvas.window().statusBar().setVisible(visible)

function cfigure(fun::Function, figname;
                 clear=true,
                 nox=false, noy=false,
                 hide_bars=true, tight=true,
                 kwargs...)
    toggle_toolbar(!hide_bars)
    fig = figure(figname; kwargs...)
    clear && clf()
    fun()
    ax_common(nox=nox, noy=noy)
    toggle_statusbar(!hide_bars, fig)
    tight && tight_layout()
    fig
end

function csubplot(fun::Function, args...; projection=nothing, kwargs...)
    ax = subplot(args...; projection=projection)
    fun()
    ax_common(;kwargs...)
    ax
end

# * Colormaps

include("colormaps.jl")

function colorbar_hack(v::AbstractVector,cmap)
    Z = [0 0
         0 0]
    CS3 = contourf(Z, v, cmap=cmap)
    clf()
    CS3
end

# * Map plots

function filter_kwargs!(kwargs, sym, default=nothing)
    v = get(kwargs, sym, default)
    filter!(kv -> kv[1] != sym, kwargs)
    v
end

# Slightly ugly hack since `PyCall.size(o::PyObject)` returns the
# _length_ of `o`.

function shape(pyobj::PyCall.PyObject)
    :shape ∈ keys(pyobj) ||
        throw(ArgumentError("`:shape` not present among `PyObject`s `keys`."))
    pyobj.shape
end
shape(pyobj::PyCall.PyObject, i) = shape(pyobj)[i]
shape(a, args...) = size(a, args...)

function plot_map(args...; kwargs...)
    kwargs = Dict{Symbol,Any}(kwargs)
    set_default!(key, val) = (kwargs[key] = get(kwargs, key, val))
    kwargs[:cmap] = get_cmap(get(kwargs, :cmap, "viridis"))

    set_default!(:shading, "auto")

    plot_fun = if filter_kwargs!(kwargs, :contour, false)
        contourf
    else
        set_default!(:rasterized, true)
        pcolormesh
    end

    vmin = filter_kwargs!(kwargs, :vmin, nothing)
    vmax = filter_kwargs!(kwargs, :vmax, nothing)
    kwargs[:norm] = (get(kwargs, :norm, :lin) == :lin ?
                     COL.Normalize(vmin, vmax) :
                     COL.LogNorm(vmin, vmax))
    aw = filter_kwargs!(kwargs, :align_ticks, false)
    xtl = filter_kwargs!(kwargs, :xticklabels, nothing)
    ytl = filter_kwargs!(kwargs, :yticklabels, nothing)
    tickrotation = filter_kwargs!(kwargs, :tickrotation, 0)
    xtickrotation = filter_kwargs!(kwargs, :xtickrotation, tickrotation)
    ytickrotation = filter_kwargs!(kwargs, :ytickrotation, tickrotation)
    xticker,yticker = if aw != false
        x,y = if length(args) == 3
            x,y,z = args
            x,y
        elseif length(args) == 1
            z = args[1]
            m,n = shape(z)
            x = 1:n
            y = 1:m
            x,y
        else
            @warn "Don't know how to shift x/y axes in case of $(length(args)) `args`"
            nothing,nothing
        end
        xt = !isnothing(x) && aw != :y ? () -> xticks(x .- (x[2]-x[1])/2,
                                                     isnothing(xtl) ? x : xtl,
                                                     rotation=xtickrotation) : () -> ()
        yt = !isnothing(y) && aw != :x ? () -> yticks(y .- (y[2]-y[1])/2,
                                                     isnothing(ytl) ? x : xtl,
                                                     rotation=ytickrotation) : () -> ()
        xt,yt
    else
        () -> (), () -> ()
    end
    p = plot_fun(args...; kwargs...)
    margins(0,0)
    xticker()
    yticker()
    p
end

function plot_polar_map(r::AbstractVector, ϕ::AbstractVector, V::AbstractMatrix,
                        args...; kwargs...)
    nϕ = length(ϕ)
    R = repeat(r,1,nϕ)
    Φ = repeat(ϕ',length(r),1)
    X = R.*cos.(Φ)
    Y = R.*sin.(Φ)
    plot_map(X,Y,V, args...; kwargs...)
end

plot_polar_map(r::AbstractVector, V::AbstractMatrix, nϕ::Integer, args...; kwargs...) =
    plot_polar_map(r, range(0,stop=2π,length=nϕ), V, args...; kwargs...)

plot_polar_map(r::AbstractVector, v::AbstractVector, nϕ::Integer, args...; kwargs...) =
    plot_polar_map(r, repeat(v,1,nϕ), nϕ, args...; kwargs...)

# function spherical_harmonic_plot(fun::Function,
#                                  r::AbstractVector,
#                                  J::AbstractVector{Int},
#                                  v::AbstractMatrix,
#                                  nθ::Int,
#                                  args...;
#                                  wfn = false,
#                                  J_max = Inf,
#                                  m_min = 0, m_max = Inf,
#                                  kwargs...)
#     nr = length(r)
#     (nr,length(J)) == size(v) || (wfn && (nr,length(J)^2) == size(v)) || error("Dimension mismatch!")

#     θ = range(0,stop=2π,length=nθ)
#     cosθ = cos.(θ)

#     R = repmat(r,1,nθ)
#     Θ = repmat(θ',nr,1)

#     X = R.*cos.(Θ)
#     Y = R.*sin.(Θ)

#     V = zeros(eltype(v), nr, nθ)

#     leg_fun = (wfn ? sf_legendre_sphPlm : sf_legendre_Plm)

#     if size(v,2) == length(J)
#         for j in eachindex(J)
#             J[j] > J_max && break
#             PJ = map(x -> leg_fun(J[j], 0, x), cosθ)'
#             V += broadcast(*, v[:,j], PJ)
#         end
#     elseif !wfn
#         error("3d plot only implemented for wavefunctions")
#     else
#         j = 0
#         for ell in J
#             ell > J_max && break
#             for m in -ell:ell
#                 j += 1
#                 abs(m) > m_max && continue
#                 abs(m) < m_min && continue
#                 Pell_m = (m<0 ? (-1)^abs(m) : 1
#                           )*map(x -> leg_fun(ell, abs(m), x), cosθ)'
#                 V += broadcast(*, v[:,j], Pell_m)
#             end
#         end
#     end

#     plot_map(X, Y, fun.(V), args...; kwargs...)
# end
# spherical_harmonic_plot(r::AbstractVector,
#                         J::AbstractVector{Int},
#                         v::AbstractMatrix,
#                         nθ::Int,
#                         args...; kwargs...) = spherical_harmonic_plot(identity,
#                                                                       r::AbstractVector,
#                                                                       J::AbstractVector{Int},
#                                                                       v::AbstractMatrix,
#                                                                       nθ,
#                                                                       args...; kwargs...)
# * Matrix plots

function plot_matrix(A::AbstractMatrix{T}, args...; align_ticks=true, bias=0.5, kwargs...) where T
    vals = filter(!iszero, unique(A))
    if T <: Complex
        @warn "Plotting real part of complex-valued matrix"
        vals = real(vals)
    end
    vmin,vmax = isempty(vals) ? (1,1) : extrema(vals)
    if vmin == vmax
        δ = 0.1abs(vmin)
        vmin = vmin - 2bias*δ
        vmax = vmax + 2*(1-bias)*δ
    end
    aa = pycall(masked_array.masked_equal, Any, Matrix(A), 0)
    plot_map(aa, args...; align_ticks=align_ticks, vmin=vmin, vmax=vmax, kwargs...)
    gca().invert_yaxis()
    square_axs()
end

# # Inspired by/stolen from
# # https://github.com/tonysyu/mpltools/blob/master/mpltools/special/hinton.py

# @pydef mutable struct SquareCollection <: pyplot_collections.RegularPolyCollection
#     __init__(self; kwargs...) = begin
#         pyplot_collections.RegularPolyCollection.__init__(
#             self, 4, rotation=pi/4; kwargs...
#         )
#     end
#     get_transform(self) = begin
#         """Return transform scaling circle areas to data space."""
#         ax = self.axes
#         pts2pixels = 72.0 / ax.figure.dpi
#         scale_x = pts2pixels * ax.bbox.width / ax.viewLim.width
#         scale_y = pts2pixels * ax.bbox.height / ax.viewLim.height
#         matplotlib.transforms.Affine2D().scale(scale_x, scale_y)
#     end
# end

# function hinton_plot_matrix(a::AbstractMatrix;
#                             max_weight=nothing,
#                             pos_color=nothing,
#                             neg_color=nothing,
#                             bg_color="none",
#                             traditional=false)
#     if traditional
#         pos_color = "white"
#         neg_color = "black"
#         bg_color = "gray"
#     else
#         pos_color == nothing && (pos_color="C0")
#         neg_color == nothing && (neg_color="C1")
#     end

#     ax = gca()
#     ax.set_facecolor(bg_color)

#     height, width = size(a)
#     if max_weight == nothing
#         max_weight = 2^ceil(log2(maximum(abs, a)))
#     end

#     vals = clamp.(a/max_weight, -1, 1)
#     neg = vals .< 0
#     pos = vals .> 0

#     cols,rows = meshgrid(1:width,1:height)

#     for (sel,color) in zip([neg,pos], [neg_color,pos_color])
#         if any(sel)
#             xy = collect(zip(cols[sel], rows[sel]))
#             circle_areas = π/2 * abs.(vals[sel])
#             squares = SquareCollection(sizes=circle_areas,
#                                        offsets=xy, transOffset=ax.transData,
#                                        facecolor=color, edgecolor=color)
#             ax.add_collection(squares, autolim=true)
#         end
#     end
#     square_axs()
#     grid(false)
#     ax.set_xlim(0.5, width+0.5)
#     ax.set_ylim(0.5, height+0.5)
#     gca().invert_yaxis()
#     ax
# end
# hinton_plot_matrix(a;kwargs...) = hinton_plot_matrix(full(a);kwargs...)

# * Patches

function draw_patch(kind, args...; ax=gca(), kwargs...)
    p = getproperty(matplotlib.patches, kind)
    ax.add_patch(p(args...; kwargs...))
end

draw_ellipse((x,y), w, h; kwargs...) = draw_patch(:Ellipse, (x,y), w, h; kwargs...)
draw_circle((x,y), r; kwargs...) = draw_ellipse((x,y), r, r; kwargs...)
draw_arc((x,y), w, h, ϕ, θ₁, θ₂; kwargs...) =
    draw_patch(:Arc, (x,y), w, h, ϕ, θ₁, θ₂; kwargs...)
draw_arc((x,y), r, ϕ, θ₁, θ₂; kwargs...) =
    draw_arc((x,y), r, r, ϕ, θ₁, θ₂; kwargs...)

# * LaTeX/font setup

function set_pgf_to_pdf(preamble=[]; texsystem = "xelatex")
    matplotlib.rcdefaults()
    ion()
    matplotlib.backend_bases.register_backend("pdf", pgf_backend.FigureCanvasPgf)

    set_latex_serif()
    rc("pgf";
       texsystem = texsystem,
       preamble = preamble)
end

function set_font(; kwargs...)
    if :serif in keys(kwargs)
        kwargs[:family] = "serif"
    elseif :sans_serif in keys(kwargs)
        kwargs[:family] = "sans-serif"
    end

    kwargs = Dict(Symbol(replace(string(k), "_" => "-")) => v
                  for (k,v) in kwargs)

    rc("font"; kwargs...)
end

function set_times_new_roman(; kwargs...)
    set_font(serif = "Times New Roman"; kwargs...)
    rc("mathtext";
       rm = "serif",
       it = "serif:italic",
       bf = "serif:bold",
       fontset = "stix")
end

function set_latex_serif(; kwargs...)
    rc("text";
       usetex = true)
    rc("font";
       family = "serif")
    rc("pgf";
       rcfonts = false)
    rc("mathtext";
       rm = "serif",
       it = "serif:italic",
       bf = "serif:bold",
       fontset = "cm")
    set_font(; kwargs...)
end

function latex(o)
    io = IOBuffer()
    show(io, MIME"text/latex"(), o)
    String(take!(io))
end

# * Scientific notation

function latex_base10(v)
    if v==0
        return "0"
    end
    ex = floor(log10(abs(v)))
    v /= 10^ex
    if ex == 0
        @sprintf("%0.7f",v)
    elseif ex == 1
        @sprintf("%0.7f\\cdot10",v)
    else
        @sprintf("%0.7f\\cdot10^{%i}",v,ex)
    end
end
latex_base10(v::Vector) = map(latex_base10, v)

function base10(v)
    if v==0
        return "0"
    end
    ex = floor(log10(abs(v)))
    v /= 10^ex
    ex = floor(Int, ex)
    if v==1
        @sprintf("10%s",to_superscript(ex))
    else
        @sprintf("%0.3f×10%s",v,to_superscript(ex))
    end
end
base10(v::Vector) = map(base10, v)

# * Axes/tick(label)s

function axis_add_ticks(ticks, labels, ax = :x; kwargs...)
    a = gca()
    a2 = a[ax == :x ? :twiny : :twinx]()
    a2[ax == :x ? :set_xlim : :set_ylim](a[ax == :x ? :get_xlim : :get_ylim]())
    a2[ax == :x ? :set_xscale : :set_yscale](a[ax == :x ? :get_xscale : :get_yscale]())
    a2[ax == :x ? :set_xticks : :set_yticks](vec(collect(ticks)))
    a2[ax == :x ? :set_xticklabels : :set_yticklabels](vec(collect(labels)); kwargs...)
    sca(a)
end

function set_ticklabel_props(ax=:x; kwargs...)
    labels = gca()[ax == :x ? :get_xticklabels : :get_yticklabels]()
    setp(labels; kwargs...)
end

function π_frac_string(i; pi_sym = "\\pi")
    den_str = v -> denominator(v) != 1 ? "/$(denominator(v))" : ""
    if i == 0
        L"$0$"
    elseif numerator(i) == 1
        latexstring("\$$(pi_sym)$(den_str(i))\$")
    elseif numerator(i) == -1
        latexstring("\$-$(pi_sym)$(den_str(i))\$")
    else
        latexstring("\$$(numerator(i))$(pi_sym)$(den_str(i))\$")
    end
end

function π_labels(ax = :x; max_count = 10, divisor=4,
                  pi_sym = "\\pi", ca=gca())
    lims = ca[ax == :x ? :get_xlim : :get_ylim]()
    f = divisor
    mi,ma = map(l -> trunc(Int, l/(π/divisor)), lims)
    if ma-mi > max_count
        f = 1
        mi,ma = map(l -> trunc(Int, l/π), lims)
    end
    d = round(Int, max((ma-mi)/max_count,1))
    r = (mi:d:ma)//f
    ca[ax == :x ? :set_xticks : :set_yticks](collect(r)*π)

    tick_labels = map(r) do i
        π_frac_string(i, pi_sym = pi_sym)
    end
    ca[ax == :x ? :set_xticklabels : :set_yticklabels](tick_labels)
end

function frac_ticks(ts::AbstractVector{Rational{Int}}, axis = :x; sfrac = false)
    if axis == :x
        xticks(ts)
    else
        yticks(ts)
    end
    tls = map(ts) do t
        if t == 0
            L"0"
        elseif denominator(t)==1
            latexstring("$(numerator(t))")
        else
            if sfrac
                latexstring("\\sfrac{$(numerator(t))}{$(denominator(t))}")
            else
                latexstring("$(numerator(t))/$(denominator(t))")
            end
        end
    end
    gca()[axis == :x ? :set_xticklabels : :set_yticklabels](tls)
end

function sci_ticks(ax, min=0, max=0)
    ticklabel_format(style="sci", axis="$ax", scilimits=(min,max))
end

function colorbar_sci_ticks(cb, min=0, max=0)
    cb.formatter.set_scientific(true)
    cb.formatter.set_powerlimits((min, max))
    cb.update_ticks()
end

function square_axs()
    gca().set_aspect("equal")
    gca().autoscale(tight=true)
end

function axes_labels_opposite(axis = :y, ax = gca())
    if axis == :y
        ax.yaxis.tick_right()
        ax.yaxis.set_label_position("right")
    else
        ax.xaxis.tick_top()
        ax.xaxis.set_label_position("top")
    end
end

function no_tick_labels(axis = :x, ax = gca(), ticks = false)
    a,b,la,lb = if axis == :x
        :bottom,:top,:labelbottom,:labeltop
    else
        :left,:right,:labelleft,:labelright
    end
    ax.tick_params(; which="both",
                   Dict(a => ticks, b => ticks,
                        la => false, lb => false)...)
end

# ** Sqrt axes

mticker = matplotlib.ticker

@pydef mutable struct SqrtTransform <: matplotlib.transforms.Transform
    input_dims = 1
    output_dims = 1
    is_separable = true
    has_inverse = true
    __init__(self) = matplotlib.transforms.Transform.__init__(self)
    transform_non_affine(self, a) = begin
        v = similar(a)
        sel = a .>= 0
        v[sel] .= .√(a[sel])
        # v[.!sel] .= NaN
        v
    end
    inverted(self) = SquareTransform()
end

@pydef mutable struct SquareTransform <: matplotlib.transforms.Transform
    input_dims = 1
    output_dims = 1
    is_separable = true
    has_inverse = true
    __init__(self) = matplotlib.transforms.Transform.__init__(self)
    transform_non_affine(self, a) = a .^ 2
    inverted(self) = SqrtTransform()
end

@pydef mutable struct SqrtScale <: matplotlib.scale.ScaleBase
    name="sqrt"
    __init__(self, axis, args...; kwargs...) = begin
        matplotlib.scale.ScaleBase.__init__(axis)
    end
    get_transform(self) = SqrtTransform()
    set_default_locators_and_formatters(self, axis) = nothing
    limit_range_for_scale(self, vmin, vmax, minpos) = begin
        max(vmin, 0), max(vmax, 0)
    end
end

matplotlib.scale.register_scale(SqrtScale)

# * Misc

pyslice(args...) = pycall(pybuiltin("slice"), PyObject, args...)

function savefig_f(filename, args...; kwargs...)
    savefig(filename, args...; kwargs...)
    filename
end

reltext(x,y,string,args...;kwargs...) = text(x,y,string,transform=gca().transAxes,args...;kwargs...)

function disp()
    display(gcf())
    println()
end

GridSpec = matplotlib.gridspec.GridSpec

# https://stackoverflow.com/a/51207905/1079038
function next_color(ax=gca())
    pc = ax._get_lines.prop_cycler
    first(pc)["color"]
end

# * ICC support
include("save_pgf_with_icc.jl")

# * Exports

export plot_style,
    toggle_toolbar, toggle_statusbar,
    cfigure, csubplot,
    colormaps, colorbar_hack,
    plot_map, plot_polar_map, spherical_harmonic_plot, plot_matrix, hinton_plot_matrix,
    draw_patch, draw_ellipse, draw_circle, draw_arc,
    set_pgf_to_pdf, set_font, set_times_new_roman, set_latex_serif,
    latex, latex_base10, base10,
    axis_add_ticks, set_ticklabel_props, π_frac_string, π_labels, frac_ticks, sci_ticks, colorbar_sci_ticks,
    square_axs, axes_labels_opposite, no_tick_labels,
    pyslice, savefig_f, reltext, disp,
    GridSpec, next_color

end


text:
export mv_avg_2d, delta_h, envelope, mute, mute!, laplacian, offset_map

"""
    mv_avg_2d(x; k=5)

2D moving average with a square window of width k
"""
function mv_avg_2d(x::AbstractArray{T, 2}; k=5) where T
    out = 1f0 * x
    nx, ny = size(x)
    kl, kr = k÷2 + 1, k÷2
    for i=kl:nx-kr, j=kl:ny-kr
        out[i, j] = sum(x[i+k, j+l] for k=-kr:kr, l=-kr:kr)/k^2
    end
    out
end

"""
    envelope(x)

Envelope of a vector or a 2D matrix. The envelope over the first dimension is taken for a 2D matrix (see DSP `hilbert`)
"""
envelope(x::AbstractArray{T, 1}) where T = abs.(hilbert(x))

"""
    envelope(x)

Envelope of a vector or a 1D vector (see DSP `hilbert`)

"""
envelope(x::AbstractArray{T, 2}) where T = abs.(hilbert(x))

"""
    delta_h(ha, h, tol)

Compute the binary mask where `ha` is within `tol` of `h`.
"""
delta_h(ha::AbstractArray{T, 2}, h::Number, tol::Number) where T = Float32.(abs.(h .- ha) .<= tol)

"""
    mute!(shot, offsets;vp=1500, t0=1/10, dt=0.004)

In place direct wave muting of a shot record with water sound speed `vp`, time sampling `dt` and firing time `t0`.
"""
function mute!(shot::AbstractArray{Ts, 2}, offsets::Vector{To}; vp=1500, t0=1/10, dt=.004) where {To, Ts}
    length(offsets) == size(shot, 2) || throw(DimensionMismatch("Number of offsets has to match the number of traces"))
    inds = trunc.(Integer, (offsets ./ vp .+ t0) ./ dt)
    inds = min.(max.(1, inds), size(shot, 1))
    for (rx, i) = enumerate(inds)
        shot[1:i, rx] .= 0f0
    end
end

"""
    mute(shot, offsets;vp=1500, t0=1/10, dt=0.004)

Direct wave muting of a shot record with water sound speed `vp`, time sampling `dt` and firing time `t0`.
"""
function mute(shot::AbstractArray{Ts, 2}, offsets::Vector{To}; vp=1500, t0=1/10, dt=.004) where {To, Ts}
    out = Ts(1) .* shot
    mute!(out, offsets; vp=vp, t0=t0, dt=dt)
    out
end

"""
    laplacian(image; hx=1, hy=1)

2D laplacian of an image with grid spacings (hx, hy)
"""
function laplacian(image::AbstractArray{T, 2}; hx=1, hy=1) where T
    scale = 1/(hx*hy)
    out = 1 .* image
    @views begin
        out[2:end-1, 2:end-1] .= -4 .* image[2:end-1, 2:end-1]
        out[2:end-1, 2:end-1]  .+= image[1:end-2, 2:end-1] + image[2:end-1, 1:end-2]
        out[2:end-1, 2:end-1]  .+= image[3:end, 2:end-1] + image[2:end-1, 3:end]
    end
    return scale .* out
end

"""
    offset_map(rtm, rtmo; scale=0)

Return the regularized least-square division of rtm and rtmo. The regularization consists of the envelope and moving average
followed by the least-square division ['surface_gather'](@ref)
"""
function offset_map(rtm::AbstractArray{T, 2}, rtmo::AbstractArray{T, 2}; scale=0) where T
    rtmn = mv_avg_2d(envelope(rtm))
    rtmo = mv_avg_2d(envelope(rtmo))

    offset_map = rtmn .* rtmo ./ (rtmn .* rtmn .+ eps(Float32)) .- scale
    return offset_map
end


text:
## types for GASS
##

###########################
mutable struct  inputCfg
    file_parameters::String
    folder_results::String
    verbose::Bool
end

##########################
mutable struct  observation
    Array_Configuration_File::String
    Result_Folder::String
    Display_Verbose::Bool
    Observatory_Latitude::Float64
    Source_Declination::Float64
    Source_Hour_Angle::Float64
    Antenna_Number::Int
    Subarray_Number::Int
end

############################
mutable struct  subarrayParameters
        Pads_Per_Subarray::Vector{Int}
        Subarray_Name::Vector{String}
        Subrange::Array{UnitRange{Int64},1}
        Spatial_Resolution::Vector{Float64}
        Maximum_Recoverable_Scale::Vector{Float64}
        Elongation::Vector{Float64}
        Sidelobe_Level::Vector{Float64}
end

###############################
mutable struct  weight
    Weight_Subarray::Vector{Float64}
    Weight_Spatial_Resolution::Vector{Float64}
    Weight_Maximum_Recoverable_Scale::Vector{Float64}
    Weight_Elongation::Vector{Float64}
    Weight_Sidelobe_Levels::Vector{Float64}
end

############################
mutable struct GA
    Number_Iterations::Int
    Population_Size::Int
    Termination_Condition::Bool
    Threshold::Float64
    Mutation_Rate::Float64
    Tournament_Size::Int
    Number_Elitism::Int
end

############################
############################
mutable struct cfg
    arr::AbstractDataFrame
    obs::observation
    sub::subarrayParameters
    wei::weight
    ga::GA
end

############################
## size are in arcsecs
struct synthbeam
    bx::Float64       ## bx inarcsec
    by::Float64       ## by i arcsec
    ar::Float64       ## angular resolution in arcsec
    e::Float64        ## excentricity
    sidelobe::Float64 ## sidelobe levels in percentage
end

#############################
#############################
mutable struct population
    age::Int32                        ## age
    subarr::Array{Array{Int,1},2}     ## {Population,Subindices}
    fitness::Array{Float64,2}         ## fitness of each subarray
    score::Array{Float64,1}           ## global score of each set
    param::Array{Dict{String,Float64},2}    ## beam and mrs for each subarray
end

max_stars_repo_name: bosscha/GASS

max_issues_repo_name: bosscha/GASS

max_forks_repo_name: bosscha/GASS


text:
module SurfaceGenerator

    include("generator.jl")
    include("api.jl")

end

max_stars_repo_name: foldfelis/Replotter

max_issues_repo_name: foldfelis/Replotter

max_forks_repo_name: foldfelis/Replotter