text:
module IncrementalInference

# @info "Multithreaded  convolutions possible, Threads.nthreads()=$(Threads.nthreads()).  See `addFactor!(.;threadmodel=MultiThreaded)`."

using Distributed
using Requires
using Reexport

@reexport using Distributions
@reexport using KernelDensityEstimate
@reexport using ApproxManifoldProducts
# @reexport using Graphs
@reexport using LinearAlgebra

using Manifolds

export ℝ, AbstractManifold
export ProductRepr
# common groups -- preferred defaults at this time.
export TranslationGroup, RealCircleGroup
# common non-groups -- TODO still teething problems to sort out in IIF v0.25-v0.26.
export Euclidean, Circle

import NLsolve
import NLSolversBase
import Optim

using
  Dates,
  TimeZones,
  DistributedFactorGraphs,
  DelimitedFiles,
  Statistics,
  Random,
  StatsBase,
  BSON,
  FileIO,
  ProgressMeter,
  DocStringExtensions,
  FunctionalStateMachine,
  JSON2,
  Combinatorics,
  UUIDs,
  TensorCast

using StaticArrays

using ManifoldsBase

# for BayesTree
using MetaGraphs

using Logging

# bringing in BSD 3-clause ccolamd
include("ccolamd.jl")
using SuiteSparse.CHOLMOD: SuiteSparse_long # For CCOLAMD constraints.
using .Ccolamd

# likely overloads or not exported by the upstream packages
import Base: convert, ==
import Distributions: sample
import Random: rand, rand!
import KernelDensityEstimate: getBW
import KernelDensityEstimate: getPoints
import ApproxManifoldProducts: kde!, manikde!
import ApproxManifoldProducts: getBW
import ApproxManifoldProducts: mmd
import ApproxManifoldProducts: isPartial
import DistributedFactorGraphs: reconstFactorData
import DistributedFactorGraphs: addVariable!, addFactor!, ls, lsf, isInitialized
import DistributedFactorGraphs: compare, compareAllSpecial
import DistributedFactorGraphs: rebuildFactorMetadata!
import DistributedFactorGraphs: getDimension, getManifold, getPointType, getPointIdentity
import DistributedFactorGraphs: getPPE, getPPEDict
import DistributedFactorGraphs: getFactorOperationalMemoryType
import DistributedFactorGraphs: getPoint, getCoordinates
import DistributedFactorGraphs: getVariableType
import DistributedFactorGraphs: AbstractPointParametricEst, loadDFG
import DistributedFactorGraphs: getFactorType

# will be deprecated in IIF
import DistributedFactorGraphs: isSolvable

# must be moved to their own repos
const KDE = KernelDensityEstimate
const MB = ManifoldsBase
const AMP = ApproxManifoldProducts
const FSM = FunctionalStateMachine
const IIF = IncrementalInference

const InstanceType{T} = Union{Type{<:T},<:T}
const NothingUnion{T} = Union{Nothing, <:T}
const BeliefArray{T} = Union{<:AbstractMatrix{<:T}, <:Adjoint{<:T, AbstractMatrix{<:T}} } # TBD deprecate?

## =============================
# API Exports

# Package aliases
export KDE, AMP, DFG, FSM, IIF

# TODO temporary for initial version of on-manifold products
KDE.setForceEvalDirect!(true)

include("ExportAPI.jl")

## =============================
# Source code

# FIXME, move up to DFG
# abstract type AbstractManifoldMinimize <: AbstractRelative end

# regular
include("entities/SolverParams.jl")

include("entities/FactorOperationalMemory.jl")

include("Factors/GenericMarginal.jl")
# Special belief types for sampling as a distribution
include("entities/AliasScalarSampling.jl")
include("entities/OptionalDensities.jl")
include("entities/BeliefTypes.jl")
include("entities/FactorGradients.jl")

# Statistics helpers on manifolds
include("VariableStatistics.jl")

# factors needed for belief propagation on the tree
include("Factors/MsgPrior.jl")

include("entities/CliqueTypes.jl")
include("entities/JunctionTreeTypes.jl")

include("services/GraphInit.jl")
include("FactorGraph.jl")
include("services/BayesNet.jl")

# Serialization helpers
include("Serialization/entities/SerializingDistributions.jl")
include("Serialization/entities/AdditionalDensities.jl")
include("Serialization/services/SerializingDistributions.jl")
include("Serialization/services/SerializationMKD.jl")
include("Serialization/services/DispatchPackedConversions.jl")

include("FGOSUtils.jl")
include("CompareUtils.jl")
include("NeedsResolution.jl")

# tree and init related functions
include("SubGraphFunctions.jl")
include("JunctionTree.jl")
include("TreeMessageAccessors.jl")
include("TreeMessageUtils.jl")
include("TreeBasedInitialization.jl")

# included variables of IIF, easy to extend in user's context
include("Variables/DefaultVariables.jl")
include("Variables/Circular.jl")

# included factors, see RoME.jl for more examples
include("Factors/GenericFunctions.jl")
include("Factors/Mixture.jl")
include("Factors/DefaultPrior.jl")
include("Factors/LinearRelative.jl")
include("Factors/EuclidDistance.jl")
include("Factors/Circular.jl")
include("Factors/PartialPrior.jl")
include("Factors/PartialPriorPassThrough.jl")
include("DefaultNodeTypes.jl") # older file

# Refactoring in progress
include("services/CalcFactor.jl")
# gradient tools
include("services/FactorGradients.jl")
include("services/CliqueTypes.jl")

# solving graphs
include("SolverUtilities.jl")
include("NumericalCalculations.jl")
include("DeconvUtils.jl")
include("ExplicitDiscreteMarginalizations.jl")
# include("InferDimensionUtils.jl")
include("services/EvalFactor.jl")
include("services/ApproxConv.jl")

include("ConsolidateParametricRelatives.jl") # FIXME CONSOLIDATE

include("GraphProductOperations.jl")
include("SolveTree.jl")
include("TetherUtils.jl")
include("TreeDebugTools.jl")
include("CliqStateMachineUtils.jl")

#EXPERIMENTAL parametric
include("ParametricCSMFunctions.jl")
include("ParametricUtils.jl")

#X-stroke
include("CliqueStateMachine.jl")

include("CanonicalGraphExamples.jl")

include("AdditionalUtils.jl")
include("SolverAPI.jl")

# Symbolic tree analysis files.
include("AnalysisTools.jl")

include("ManifoldSampling.jl")

# deprecation legacy support
include("Deprecated.jl")

exportimg(pl) = error("Please do `using Gadfly` to allow image export.")
function __init__()
  @require InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" include("RequireInteractiveUtils.jl")
  @require Gadfly="c91e804a-d5a3-530f-b6f0-dfbca275c004" include("EmbeddedPlottingUtils.jl")
  @require DifferentialEquations="0c46a032-eb83-5123-abaf-570d42b7fbaa" include("ODE/DERelative.jl")
  @require Interpolations="a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" include("services/HeatmapSampler.jl")

  # combining neural networks natively into the non-Gaussian  factor graph object
  @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" begin
    include("Flux/FluxModelsDistribution.jl")
    include("Serialization/entities/FluxModelsSerialization.jl")
    include("Serialization/services/FluxModelsSerialization.jl") # uses BSON
  end
end

export setSerializationNamespace!, getSerializationModule, getSerializationModules

end


text:
module QMMM2

include("data.jl")

include("train.jl")

include("calculators.jl")

include("solve.jl")

include("sketches.jl")

# include("edip.jl")

end # module

ext: jl


text:
# 20200816
using Dates
using Printf

"""
是否回文
"""
function ishuiwen(n::Int64)
    sn::Int32 = 0 # 可以不用类型
    sn = n
    tn::Int32 = 0
    while sn != 0
        tn = tn * 10 + sn % 10
        sn = sn ÷ 10
    end
    if tn == n
        return true
    end
    return false
end

"""
回文10000000
"""
function huiwen(max::Int64=10000000)
    tx::Int32 = 0
    x::Int32 = 0
    for i = 0:max
        if ishuiwen(i) == true
            tx+=1
        end
    end
    return tx
end

function run()
    t1 = now()
    r = huiwen(10000000)
    t2 = now()
    tx = t2 - t1
    @printf("回文 结果:: %d == 时间:: %F \n", r, Dates.value(tx))
end

run()
# julia .\huiwen.jl; julia --version


text:
# adds a registered solver
function add(solver_name::AbstractString)
    @assert solver_name in SUPPORTED_SOLVERS string("The solver: ", solver_name, " is not supported")
    full_url = string(REMOTE_URL, solver_name, ".jl")
    Pkg.clone(full_url)
    Pkg.build(solver_name)
end

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
using Pkg
Pkg.add("InfoZIP")
Pkg.add("ExcelReaders")
Pkg.add("DataFrames")
Pkg.add("CSVFiles")
Pkg.add("CSV")
Pkg.add("StatsBase")
Pkg.add("Distributions")

using Mimi
using InfoZIP

function isa_url(x)
    return startswith(x, "https:")
end

#list of URLs of branches of packages to test
dependencies = [
    "https://github.com/fund-model/fund/archive/1768edf12aaaac3a41bbea081d5b51299121f993.zip",
    "https://github.com/anthofflab/mimi-rice-2010.jl/archive/2b5996b0a0c8be92290991f045c43af425c5a9c8.zip"
]

function run_dependency_tests(dependencies=dependencies)
    #list of failed tests to build as you go
    errors = []
    #make a temporary directory to run the tests in
    tmp_path = joinpath(@__DIR__,"tmp_testing/")
    mkdir(tmp_path)

    #loop through each dependent package
    for d in dependencies
        if  isa_url(d)
            zip_name = chomp(basename(d))
            zip_file_path = joinpath(tmp_path, zip_name)
            download(d, zip_file_path)
            InfoZIP.unzip(zip_file_path, tmp_path)
            rm(zip_file_path)
            #find the name of the unzipped package (this only works if the zip archive only has one directory, the package)
            package_name = readdir(tmp_path)[1]
            file_path = string(tmp_path, package_name)
        else
            package_name = basename(d)
            file_path = d
        end

        #first check for mimitests.jl, if not found default to runtests.jl
        if "mimitests.jl" in readdir(string(file_path, "/test/"))
            process = string(file_path, "/test/mimitests.jl")
        else
            process = string(file_path, "/test/runtests.jl")
        end

        #test the package
        try
            run(`$(Sys.BINDIR)/julia $process`)
        catch e
            append!(errors, [(package_name, e)])
        end
        #delete current package before testing next one (if it was a downloaded package)
        if isa_url(d)
            rm(joinpath(tmp_path, package_name), recursive=true)
        end
    end

    #remove the temporary directory
    rm(tmp_path, recursive=true)

    #report the errors that occurred
    num_errors = length(errors)
    error_message = "Failed tests: $num_errors"

    for (package_name, error) in errors
        error_message = string(error_message, "\n", "error in $package_name:", error)
    end

    if num_errors > 0
        error(error_message)
    else
        @info "All dependency tests passed."
    end
end


text:
"""
    time_evo_L <: AbstractProblem

Problem or finding the steady state of a ℒdagℒ matrix by computing
𝒞 = ∑|ρ(σ)|²|⟨⟨σ|ℒ |ρ⟩⟩|² using the sparse Liouvillian matrix.
"""
struct time_evo_L{B, SM} <: OpenTimeEvolutionProblem where {B<:Basis,
                                                 SM<:SparseMatrixCSC}
    HilbSpace::B            # 0
    L::SM
    ρss
end

time_evo_L(args...) = time_evo_L(STD_REAL_PREC, args...)
time_evo_L(T::Type{<:Number}, gl::GraphLindbladian) =
    time_evo_L(T, liouvillian(gl))
time_evo_L(T::Type{<:Number}, Liouv::SparseSuperOperator) =
    time_evo_L(first(Liouv.basis_l), sparse(transpose(Liouv.data)), 0.0)
time_evo_L(T::Type{<:Number}, Ham::DataOperator, cops::Vector) =
    time_evo_L(T, liouvillian(Ham, cops))

basis(prob::time_evo_L) = prob.HilbSpace

function compute_Cloc(prob::time_evo_L, net::MatrixNet, σ, lnψ=net(σ), σp=deepcopy(σ))
    ℒ = prob.L
    set_index!(σp, index(σ))

    C_loc = zero(Complex{real(out_type(net))})
    # Iterate through all elements in row i_σ of the matrix computing
    # ⟨i_σ|ℒdagℒψ⟩ = Σ_{i_σp} ⟨i_σ|ℒdagℒ|i_σp⟩⟨i_σp|ψ⟩
    # NOTE: ℒdagℒ is CSC, but I would like a CSR matrix. Since it is hermitian I
    # can simply take the conjugate of the elements in the columns
    i_σ = index(σ)
    for row_id = ℒ.colptr[i_σ]:(ℒ.colptr[i_σ+1]-1)
      # Find nonzero elements s by doing <i_sp|ℒdagℒ|i_σ>
      i_σp = ℒ.rowval[row_id]
      # BackConvert to int
      set_index!(σp, i_σp)
      # Compute the log(ψ(σ)/ψ(σ')), by only computing differences.
      log_ratio = logψ(net, σp) - lnψ
      # Conj because I am taking the transpose... see the note above.
      C_loc  += (ℒ.nzval[row_id]) * exp(log_ratio)
    end
    return C_loc
end


text:
struct AllocError    <: Exception end
struct BoundsError   <: Exception end
struct Error         <: Exception end
struct RuntimeError  <: Exception end
struct LogicError    <: Exception end
struct ArgumentError <: Exception end

function ElError(code::Integer)
    if code == 0
        return nothing
    elseif code == 1
        throw(AllocError())
    elseif code == 2
        throw(BoundsError())
    elseif code == 3
        throw(ArgumentError())
    elseif code == 4
        throw(LogicError())
    elseif code == 5
        throw(RuntimeError())
    elseif code == -1
        # catchall error code
        return throw(Error())
    else
        throw(ArgumentError("Unknown Elemental Error Code: $code"))
    end
end

ErrorString(code::Integer) =
    String(ccall((:ElErrorString, libEl), Ptr{Cchar}, (Cint,), code))


text:
export SeisBandPass,
SeisDecimate,
SeisDelay,
SeisDiff,
SeisEnvelope,
SeisFKFilter,
SeisFXDecon,
SeisGain,
SeisKolmogoroff,
SeisMute,
SeisNMO,
SeisProcessFile,
SeisPWD,
SeisRadonFreqFor,
SeisRadonFreqInv,
SeisSemblance,
SeisSmooth1,
SeisSmooth2,
SeisSmoothGathers,
stacktraces,
SeisSincInterp1D
include("SeisBandPass.jl")
include("SeisDecimate.jl")
include("SeisDelay.jl")
include("SeisDiff.jl")
include("SeisEnvelope.jl")
include("SeisFKFilter.jl")
include("SeisFXDecon.jl")
include("SeisGain.jl")
include("SeisKolmogoroff.jl")
include("SeisMute.jl")
include("SeisNMO.jl")
include("SeisProcessFile.jl")
include("SeisPWD.jl")
include("SeisRadonFreqFor.jl")
include("SeisRadonFreqInv.jl")
include("SeisSemblance.jl")
include("SeisSmooth1.jl")
include("SeisSmooth2.jl")
include("SeisSmoothGathers.jl")
include("SeisStack.jl")
include("SeisSincInterp1D.jl")


text:
#!/usr/bin/env julia
@info "Installing Julia dependencies"

if VERSION >= v"0.7.0-DEV.2005"
    using Pkg
end

Pkg.activate(@__DIR__)
Pkg.instantiate()

@info "Installed Julia packages:"
Pkg.status()

max_stars_repo_name: cffbots/ESMValTool

max_issues_repo_name: cffbots/ESMValTool

max_forks_repo_name: cffbots/ESMValTool


text:
#-------------------------------------------------------------------
#* Model of a pump (simplified, used in distillation column model)
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#*--------------------------------------------------------------------
type pump
	pump()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of chemical components"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0,
				:PosY=>0.4727,
				:Protected=>true,
				:Symbol=>"_{in}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>1,
				:PosY=>0.1859,
				:Protected=>true,
				:Symbol=>"_{out}"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Increase",
				:Lower=>0,
				:DisplayUnit=>"kPa",
				:Symbol=>"P_{incr}"
			)),
			[
				:(Inlet.F = Outlet.F),
				:(Inlet.z = Outlet.z),
				:(Outlet.P = Inlet.P + Pincrease),
				:(Outlet.h = Inlet.h),
			],
			[
				"Molar Balance","Composittion","Pump head","Pump potency",
			],
			[:PP,:NComp,],
			[:Inlet,:Outlet,:Pincrease,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Inlet::stream
	Outlet::liquid_stream
	Pincrease::press_delta
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export pump
function setEquationFlow(in::pump)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
end
function atributes(in::pump,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Pump"
	fields[:Brief]="Model of a simplified pump, used in distillation column model."
	fields[:Info]="== ASSUMPTIONS ==
* Steady State;
* Only Liquid;
* Adiabatic;
* Isentropic.

== SPECIFY ==
* the inlet stream;
* the pump Pincrease.
"
	drive!(fields,_)
	return fields
end
pump(_::Dict{Symbol,Any})=begin
	newModel=pump()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
function usedbytes(v)
    n = 0
    for i in 1:length(v)
        n += sizeof(v[i])
    end
    nv = sizeof(v)
    if nv == n
        # we already accounted for the storage of the vectors v[i] and the
        # vector v because the elements of the vector v are stored in line
        return n
    else
        # the vectors v[i] are accounted for, but now we need to add the storage
        # for the vector v itself
        return n + nv
    end
end


text:
import .Distributions

Distributions.cdf(::Type{E}, args...; opts...) where {E<:AnyLog} = exp(E, Distributions.logcdf(args...; opts...))
Distributions.ccdf(::Type{E}, args...; opts...) where {E<:AnyLog} = exp(E, Distributions.logccdf(args...; opts...))
Distributions.pdf(::Type{E}, args...; opts...) where {E<:AnyLog} = exp(E, Distributions.logpdf(args...; opts...))


text:
# # Minimax estimation of linear functionals

# This tutorial seeks to illustrate Stein's one-dimensional subfamily heuristic for
# finding minimax optimal linear estimators.

# ## Preliminaries

# Let us load some required packages first
using MinimaxCalibratedEBayes
using MosekTools
using Setfield
using LaTeXStrings
using LaTeXTabulars

const MCEB = MinimaxCalibratedEBayes

# Some further preliminaries for the plotting setting:
using Plots
using PGFPlotsX
pgfplotsx()
def_size = (850,440)
third_col = RGB(68/255, 69/255, 145/255)
import Plots:pgfx_sanitize_string
pgfx_sanitize_string(s::AbstractString) = s

# # Recap of goal in this section
# Here we seek to demonstrate how we can solve optimization problems of the following form:
# Given a partition of the real line $(-\infty, \infty)$ into bins $I_{k}$, find
# the constant $Q_0$ and function $Q(\cdot)$ which optimize the following
# optimization problem:
#
# $$\min_{Q_0,Q} \left\{ \sup_{ G \in \mathcal{G}} \{ (Q_0 + E_{G}[Q(X_i)] - L(G))^2\} + \frac{1}{n}\int Q^2(x) \bar{f}(x) dx \right\}$$
# We solve this over all functions $Q(\cdot)$ that are piecewise constant on $I_k$. This
# corresponds to finding the linear estimator:
#
# $$\hat{L} = Q_0 + \frac{1}{n} \sum_{i=1}^n Q(X_i)$$
#
# that solves a worst-case bias, variance problem with respect
# to the following problem parameters.
# * The convex class of priors (effect size distributions) $\mathcal{G}$.
# * The linear functional$L(G)$ operating on $G \in \mathcal{G}$.
# * A pilot estimator for the marginal density of the observations $\bar{f}(\cdot)$.

# # Main plots of quasi-minimax linear estimators
# Here we reproduce Figures 3,4,5, of the paper and also
# show how the package can be used to derive quasi-minimax linear estimators.

# ## Marginal Density Target, Gaussian Mixture prior

# We start by defining the prior (effect size) distribution we will use throughout.
# This is used as the ground truth (and would typically not be available) in a real
# analysis.

prior_dbn = MixtureModel(Normal, [(-2,0.2), (2, 0.2)])

# Note that under this prior $G$, the true marginal density is the following:
marginalize(prior_dbn, StandardNormalSample(0.0))

# Next we define the discretization of the output domain by using the
# `DiscretizedStandardNormalSamples` type.
marginal_grid = -6:0.01:6
Zs_discr = DiscretizedStandardNormalSamples(marginal_grid);

# We populate the `var_proxy` field of `Zs_discr` based on the ground truth.
# That is, for $k$ and bin $I_k$ in the discretization, we compute
# $$\bar{\nu}(k) = \int_{I_k} \bar{f}(x) dx$$

Zs_discr_var = @set Zs_discr.var_proxy = pdf(marginalize(prior_dbn, Zs_discr));

# We define the linear functional we want to design a linear estimator for,
# namely the marginal density at 0.0, i.e. $f_{G}(0)$.
marginal_target = MarginalDensityTarget(StandardNormalSample(0.0))

# Furthermore we will do this based on 10000 samples.
n_marginal = 10_000

# Finally we define $\mathcal{G}$, the class of effect size distributions we consider:
gmix = GaussianMixturePriorClass(0.2, -3:0.01:3, Mosek.Optimizer, ["QUIET" =>true])

# Finally we are ready to compute the minimax estimator; the `MCEB.RMSE(n_marginal, 0.0)`
# option denotes that we want to optimize the mean-squared error trade-off based
# on the sample size `n_marginal`.

marginal_fit = SteinMinimaxEstimator(Zs_discr_var, gmix,
                                     marginal_target, MCEB.RMSE(n_marginal, 0.0));

# Statistical performance is often improved by adding constraints on the marginal
# density of the observations, i.e., when it is known that $\sup_x |f(x) - \bar{f}(x)| \leq c_m$
# for some $c_m$. This is possible to do in a data-driven way. Here we illustrate by using
# $c_m = 0.02$ around the ground truth density.

C∞ = 0.02
Zs_discr_nbhood = set_neighborhood(Zs_discr, prior_dbn; C∞ = C∞);

# Finally we recompute the minimax estimator under the additional localization constraint:
marginal_fit_nbhood = SteinMinimaxEstimator(Zs_discr_nbhood, gmix,
									 marginal_target, MCEB.RMSE(n_marginal, 0.0));

# We have computed two linear estimators for the marginal density target. A baseline to compare
# against is the Butucea-Comte estimator:

bc_marginal = MCEB.ButuceaComteEstimator(marginal_target;  n=n_marginal);

# Finally, we are ready to plot all the linear estimators, as well as the worst case
# densities. This is Figure 3 in the manuscript.

marginal_density_affine = steinminimaxplot(marginal_fit, marginal_fit_nbhood; size=def_size,
									                   ylim_relative_offset=1.4)
plot!(marginal_density_affine[1], marginal_grid, bc_marginal.(marginal_grid),
      linestyle=:dot, linecolor = third_col, label="Butucea-Comte")

#md plot(marginal_density_affine)

#jl savefig(marginal_density_affine, "marginal_density_affine.pdf")

# ## Prior Tail Probability, Gaussian Mixture prior
# Next we compute the quasi-minimax estimators in the same setting as above, with the following changes:
# The target we consider is instead the prior tail probability $\mathbb P_G[\mu_i \geq 0]$:

prior_cdf_target = PriorTailProbability(0.0)

# The sample size we consider is
n_prior_cdf = 200

# We compute the minimax estimators over the whole class $\mathcal{G}$, as well as the
# localized class $\mathcal{G}_m$.

prior_cdf_fit = SteinMinimaxEstimator(Zs_discr_var, gmix,
                                      prior_cdf_target, MCEB.RMSE(n_prior_cdf, 0.0));

prior_cdf_fit_nbhood = SteinMinimaxEstimator(Zs_discr_nbhood, gmix,
									  prior_cdf_target, MCEB.RMSE(n_prior_cdf, 0.0));

# Let us plot the results

prior_cdf_affine = steinminimaxplot(prior_cdf_fit, prior_cdf_fit_nbhood; size=def_size,
									  					 ylim_relative_offset=1.4)

#md plot(prior_cdf_affine)
#jl savefig(prior_cdf_affine, "prior_cdf_affine.pdf")

# ## Prior Density Target, Sobolev prior class

# As our next example we will keep the same sample size as above.
n_prior_density = 200
# However, the target will instead be the density at $0$, i.e., $L(G) = g(0)$:
prior_density_target = PriorDensityTarget(0.0)

# Furthermore, we will now take $\mathcal{G}$ to be the Sobolev class:
hmclass = MCEB.HermitePriorClass(qmax=90, sobolev_order=2, sobolev_bound=0.5,
                                 solver=Mosek.Optimizer, solver_params = ["QUIET" => true]);

# Finally, we will also use a different ground truth $G$
prior_dbn_normal = Normal(0, 2)

marginalize(prior_dbn_normal, StandardNormalSample(0.0))

# Let us compute again the marginal (discretized) distribution implies by prior_dbn_normal
Zs_discr_var_normal = @set Zs_discr.var_proxy = pdf(marginalize(prior_dbn_normal, Zs_discr));
Zs_discr_nbhood_normal = set_neighborhood(Zs_discr, prior_dbn_normal; C∞ = C∞);

# We now compute the minimax linear estimators (without and with localization):

prior_density_fit = SteinMinimaxEstimator(Zs_discr_var_normal, hmclass,
								          prior_density_target, MCEB.RMSE(n_prior_density,0.0));

prior_density_fit_nbhood = SteinMinimaxEstimator(Zs_discr_nbhood_normal, hmclass,
								   prior_density_target, MCEB.RMSE(n_prior_density,0.0));

# Let us also instantiate the Butucea-Comte estimator
bc_prior = MCEB.ButuceaComteEstimator(prior_density_target;  n=n_prior_density);

# Again we plot the results:

prior_density_affine = steinminimaxplot(prior_density_fit, prior_density_fit_nbhood; size=def_size,
									     ylim_relative_offset=1.3, ylim_panel_e= (0,0.33))

plot!(prior_density_affine[1], marginal_grid, bc_prior.(marginal_grid),
      linestyle=:dot, linecolor = third_col, label="Butucea-Comte")

#jl savefig(prior_density_affine, "prior_density_affine.pdf")


text:
module PermutationGroups

using AbstractAlgebra
using Markdown

import AbstractAlgebra: AbstractPermutationGroup, AbstractPerm, Group, GroupElem, Generic.Perm,
    mul!, @perm_str, degree, order, permtype, gens

export AbstractPermutationGroup, AbstractPerm, Group, GroupElem
export @perm_str, Perm, degree, order, permtype, gens

export AbstractOrbit, Orbit, Transversal, Schreier, StabilizerChain, PermGroup
export firstmoved, fixes, fixedpoints, lastmoved, nfixedpoints
export base, getinv, representative, schreier_sims, sgs, sift

include("types.jl")
include("utils.jl")
include("orbit.jl")
include("words.jl")
include("schreier.jl")
include("groups.jl")
include("stabchain.jl")
include("schreier-sims.jl")

end # module


text:
module VsRecorder

using Artifacts, Dates

using SimpleI18n
using VsRecorderBase
import Tesseract

include("io.jl")
export VsI18n
include("i18n.jl")
using .VsI18n

export PokemonBattle,
    PARSE_NONE, PARSE_MINIMAL, PARSE_OPPONENT, PARSE_BOTH_PLAYERS
include("config.jl")
include("types/base.jl")

export PokemonID, ItemID, MoveID, AbilityID, Nature,
    PokemonType,
    Gender, GENDER_NULL, GENDER_MALE, GENDER_FEMALE
include("types/identifiers.jl")
export Player, Stats, Pokemon, Team
include("types/pokemon.jl")

export poke_dex, move_dex,
    item_list, ability_list, poke_list, move_list,
    import_team, export_team
include("data/Data.jl")
using .Data

export BattleEnums, Battle, ParsedBattle
include("types/battle.jl")
include("types/context.jl")
include("types/result.jl")

include("scenes/Scenes.jl")

include("main.jl")

using Reexport: @reexport
@reexport using VsRecorderBase.API

include("boot.jl")

end # module


text:
include("test_Dirichlet.jl")
include("test_NormalInverseGamma.jl")


text:
################################################################################
# MIT License
#
# Copyright (c) 2021 Hajime Nakagami<nakagami@gmail.com>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
################################################################################

mutable struct Cursor <: DBInterface.Cursor
    stmt::Statement

    names::Vector{Symbol}
    types::Vector{Type}
    lookup::Dict{Symbol, Int}

    rows::Vector{Vector{Any}}
    current_rownumber::Int
end

struct Row <: Tables.AbstractRow
    cursor::Cursor
    rownumber::Int
end

getcursor(r::Row) = getfield(r, :cursor)
getrownumber(r::Row) = getfield(r, :rownumber)

Tables.columnnames(r::Row) = getcursor(r).names

Tables.getcolumn(r::Row, i::Int) = getcursor(r).rows[getrownumber(r)][i]
Tables.getcolumn(r::Row, nm::Symbol) = Tables.getcolumn(r, getcursor(r).lookup[nm])

Tables.isrowtable(::Type{<:Cursor}) = true
Tables.schema(c::Cursor) = Tables.Schema(c.names, c.types)

Base.eltype(c::Cursor) = Row
Base.length(c::Cursor) = length(c.rows)

function Base.iterate(cursor::Cursor, i=1)
    if length(cursor.rows) < i
        nothing
    else
        cursor.current_rownumber = i
        Row(cursor, i), i + 1
    end
end

"""
    DBInterface.close!(cursor)

Close a cursor.
"""
function DBInterface.close!(cur::Cursor)
    close!(cur.stmt)
end


text:
using FastGMM
using Documenter

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

makedocs(;
    modules=[FastGMM],
    authors="Andreas Hildebrandt <andreas.hildebrandt@uni-mainz.de> and contributors",
    repo="https://github.com/hildebrandtlab/FastGMM.jl/blob/{commit}{path}#{line}",
    sitename="FastGMM.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://hildebrandtlab.github.io/FastGMM.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/hildebrandtlab/FastGMM.jl",
    devbranch="main",
)


text:
using Dates

dir = mktempdir()

@test :xyz == file_format("/tmp/file.xyz")

txt = "ABC XYZ"
zip = joinpath(dir, "test.zip")

compress(zip, "file.txt") do io
    write(io, txt)
end

uncompress(zip) do name, io
    @test name == "file.txt"
    @test txt == read(io, String)
end

bad = joinpath(dir, "test.bad")

@test_throws ErrorException compress(bad, "file.txt") do io
end

@test_throws ErrorException uncompress(bad) do io
end

df = DataFrame(A = 1:4, B = ["M", "F", "F", "M"])

for format in df_formats()
    file = joinpath(dir, "test_df."*String(format))
    df_write(file, df; subformat=:csv)
    df2 = df_read(file)
    @test df == df2
end

for format in compression_formats()
    file = joinpath(dir, "test_df."*String(format))

    for subformat in df_formats()
        if subformat == :jld2 || subformat == :jld2c || subformat == :parquet || subformat == :jdf
            @test_throws ErrorException df_write(file, df; subformat=subformat)
            continue
        end

        df_write(file, df; subformat=subformat)
        df2 = df_read(file)
        @test df == df2
    end
end

@test_throws ErrorException df_write(bad, df)
@test_throws ErrorException df_read(bad)

# Test Issue #14
df = DataFrame()
df[!, :test] = [DateTime(2000,1,1,1,1,1), missing]
DataFrameTools.df_write(joinpath(dir, "test.jdf"), df)


text:
Module VREP

include("VREPConst/VREPConst.jl")
include("VREPFunc/VREPFunc.jl")

end


text:
## ------------------------------------------------------
mutable struct Container{T}
    v::Vector{T}
    jump_size::Int
    idx::Int

    function Container{T}(jump_size) where {T}
        @assert jump_size > 0
        new{T}(Vector{T}(undef, jump_size), jump_size, 1)
    end
    Container{T}() where {T} = Container{T}(10_000)
end

function Base.push!(cont::Container{T}, elms::T...) where {T}
    for elm in elms
        vlen = length(cont.v)
        if cont.idx <= vlen
            cont.v[cont.idx] = elm
            cont.idx += 1
        else
            resize!(cont.v, vlen + cont.jump_size)
            push!(cont, elm)
        end
    end
    return cont
end

Base.length(cont::Container) = (cont.idx - 1)
Base.firstindex(cont::Container) = firstindex(cont.v)
Base.lastindex(cont::Container) = (cont.idx - 1)
Base.isempty(cont::Container) = (cont.idx == 1)
reset!(cont::Container) = (cont.idx = firstindex(cont); cont)
Base.empty!(cont::Container) = (empty!(cont.v); reset!(cont))

concat!(conts::Container...) = vcat(vec!.(conts)...)

Base.vec(cont::Container) = view(cont.v, firstindex(cont):lastindex(cont))
Base.resize!(cont::Container) = resize!(cont.v, lastindex(cont))
vec!(cont::Container) = (resize!(cont); cont.v)
Base.iterate(cont::Container) = iterate(vec(cont))
Base.iterate(cont::Container, state) = iterate(vec(cont), state)


text:
mutable struct DDEIntegrator{algType,IIP,uType,tType,P,eigenType,absType,relType,
                             residType,tTypeNoUnits,tdirType,ksEltype,
                             SolType,F,CacheType,
                             IType,NType,O,tstopsType,
                             FSALType,EventErrorType} <: AbstractDDEIntegrator{algType,IIP,uType,tType}
    sol::SolType
    u::uType
    k::ksEltype
    t::tType
    dt::tType
    f::F
    p::P
    uprev::uType
    uprev2::uType
    tprev::tType
    prev_idx::Int
    prev2_idx::Int
    fixedpoint_abstol::absType
    fixedpoint_reltol::relType
    resid::residType # This would have to resize for resizing DDE to work
    fixedpoint_norm::NType
    max_fixedpoint_iters::Int
    saveat::tstopsType
    tracked_discontinuities::Vector{Discontinuity{tType}}
    alg::algType
    dtcache::tType
    dtchangeable::Bool
    dtpropose::tType
    tdir::tdirType
    eigen_est::eigenType
    EEst::tTypeNoUnits
    qold::tTypeNoUnits
    q11::tTypeNoUnits
    erracc::tTypeNoUnits
    dtacc::tType
    success_iter::Int
    iter::Int
    saveiter::Int
    saveiter_dense::Int
    cache::CacheType
    kshortsize::Int
    force_stepfail::Bool
    last_stepfail::Bool
    just_hit_tstop::Bool
    event_last_time::Int
    last_event_error::EventErrorType
    accept_step::Bool
    isout::Bool
    reeval_fsal::Bool
    u_modified::Bool
    opts::O
    integrator::IType
    fsalfirst::FSALType
    fsallast::FSALType

    # incomplete initialization without fsalfirst and fsallast
    function DDEIntegrator{algType,IIP,uType,tType,P,eigenType,absType,relType,
                           residType,tTypeNoUnits,
                           tdirType,ksEltype,SolType,F,CacheType,IType,
                           NType,O,tstopsType,FSALType,EventErrorType}(
                               sol,u,k,t,dt,f,p,uprev,uprev2,tprev,prev_idx,prev2_idx,
                               fixedpoint_abstol,fixedpoint_reltol,resid,fixedpoint_norm,
                               max_fixedpoint_iters,saveat,tracked_discontinuities,alg,
                               dtcache,dtchangeable,dtpropose,tdir,eigen_est,EEst,qold,
                               q11,erracc,dtacc,success_iter,iter,saveiter,saveiter_dense,
                               cache,kshortsize,force_stepfail,last_stepfail,
                               just_hit_tstop,event_last_time,last_event_error,
                               accept_step,isout,reeval_fsal,u_modified,opts,
                               integrator) where
        {algType,IIP,uType,tType,P,eigenType,absType,relType,residType,tTypeNoUnits,
        tdirType,ksEltype,
         SolType,F,CacheType,IType,NType,O,tstopsType,FSALType,EventErrorType}

        new{algType,IIP,uType,tType,P,eigenType,absType,relType,
                               residType,tTypeNoUnits,
                               tdirType,ksEltype,SolType,F,CacheType,IType,
                               NType,O,tstopsType,FSALType,EventErrorType}(
            sol,u,k,t,dt,f,p,uprev,uprev2,tprev,prev_idx,prev2_idx,fixedpoint_abstol,
            fixedpoint_reltol,resid,fixedpoint_norm,max_fixedpoint_iters,saveat,
            tracked_discontinuities,alg,dtcache,dtchangeable,dtpropose,tdir,eigen_est,
            EEst,qold,q11,erracc,dtacc,success_iter,iter,saveiter,saveiter_dense,cache,
            kshortsize,force_stepfail,last_stepfail,just_hit_tstop,accept_step,
            event_last_time,last_event_error,isout,reeval_fsal,u_modified,opts,integrator)
    end
end

function (integrator::DDEIntegrator)(t, deriv::Type=Val{0}; idxs=nothing)
    OrdinaryDiffEq.current_interpolant(t, integrator, idxs, deriv)
end

function (integrator::DDEIntegrator)(val::AbstractArray, t::Union{Number,AbstractArray},
                                     deriv::Type=Val{0}; idxs=nothing)
    OrdinaryDiffEq.current_interpolant!(val, t, integrator, idxs, deriv)
end


text:
# This file was generated, do not modify it. # hide
rand(Float16, Normal())

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
"""
Example for training on the maestro dataset for language modeling.
"""

using Flux
using Flux: onehotbatch, label_smoothing, update!
using Zygote

using Random
using FileIO
using Dates
using JLD2

using MusicTransformer
using Transformers
using Transformers.Basic
using Transformers.Datasets

Base.@kwdef struct HParams
    num_classes::Int=310
    lr::Float64=1e-6
    epochs::Int=10
    crop_len::Int=2048
    batchsize::Int=8
    shuffle::Bool=true
    α::Float64=0.1
    log_at_steps::Int=25
    checkpoint_at_steps::Int=1000
end

"""
    preprocess(data::Vector{Vector{Int}},
                crop_len::Int=2048,
                num_classes::Int=1,
                α::Float64=0.1)

Return X and Y batches of the data after preprocessing by taking random crops.
Label smoothing is applied to the targets.
"""
function preprocess(data::Vector{Vector{Int}},
                    crop_len::Int=2048,
                    num_classes::Int=310,
                    α::Float64=0.1)

    batchsize = length(data)
    batch = Matrix{Int}(undef, batchsize, crop_len+1)

    for i = 1:batchsize
        offset = length(data[i]) - crop_len
        crop_begin = rand(1:offset)
        crop_end = crop_begin + crop_len
        batch[i, :] = @view data[i][crop_begin:crop_end]
    end

    X = batch[:, 1:end-1]
    # Shift inputs right and one-hot encode them
    Y = onehotbatch(batch[:, 2:end], 1:num_classes)

    X, Float32.(label_smoothing(Y, α))
end

"""
    train!(model; hparams::HParams=HParams(), checkpoint_dir="/tmp/maestro_model")

Train the model with the given hyperparameters.
The model weights are saved after training.
"""
function train!(model; hparams::HParams=HParams(), checkpoint_dir="/tmp/maestro_model")
    !isdir(checkpoint_dir) && mkdir(checkpoint_dir)
    ps = params(model)
    opt = ADAM(hparams.lr)

    @info "Training..."
    for i = 1:hparams.epochs
        trainingstep = 1
        datareader = reader(Train, MaestroDataset(),
                        batchsize=hparams.batchsize, shuffle=hparams.shuffle)

        for data in datareader
            X, Y = preprocess(data, hparams.crop_len, hparams.num_classes, hparams.α)

            loss, back = Zygote.pullback(ps) do
                logits = model(X)
                logcrossentropy(Y, logsoftmax(logits))
            end
            grads = back(1f0)

            if trainingstep % hparams.log_at_steps == 0
                println("Epoch: $i   Step: $trainingstep   Loss: $loss")
            end

            if trainingstep % hparams.checkpoint_at_steps == 0
                save_path = joinpath(checkpoint_dir, "model_checkpoint_$(now()).jld2")
                JLD2.@save save_path weights=ps
            end

            update!(opt, ps, grads)
            trainingstep += 1
        end
    end

    save_path = joinpath(checkpoint_dir, "model_checkpoint_$(now()).jld2")
    JLD2.@save save_path weights=ps

    @info "Saved model weights at $save_path"
end

# Usage:
# model = BaselineMusicTransformer(512, 8, 2048, 6)
# hparams = HParams(epochs=15, crop_len=1028, log_at_steps=10)
# train!(model, hparams=hparams)


text:
function eosshow(io::IO, ::MIME"text/plain", Base.@nospecialize(model::EoSModel))
    print(io, typeof(model))
    length(model) == 1 && println(io, " with 1 component:")
    length(model) > 1 && println(io, " with ", length(model), " components:")
    for i in 1:length(model)
        print(io, " \"", model.components[i], "\"")
        println(io)
    end
    if hasfield(typeof(model),:params)
        paramnames = fieldnames(typeof(model.params))
        len_params = length(paramnames)
        !iszero(len_params) && print(io,"Contains parameters: ")
        firstloop = true
        for fieldname in paramnames
            firstloop == false && print(io, ", ")
            print(io, fieldname)
            firstloop = false
        end
    end
end

function eosshow(io::IO, Base.@nospecialize(model::EoSModel))
    print(io, typeof(model))
    firstloop = true
    print(io, "(")
    for i in 1:length(model.components)
        firstloop == false && print(io, ", ")
        print(io, "\"", model.components[i], "\"")
        firstloop = false
    end
    print(io, ")")
end

function gc_eosshow(io::IO, ::MIME"text/plain", Base.@nospecialize(model::EoSModel))
    print(io, typeof(model))
    length(model) == 1 && println(io, " with 1 component:")
    length(model) > 1 && println(io, " with ", length(model), " components:")
    for i in 1:length(model)
        print(io, " \"", model.components[i], "\": ")
        firstloop = true
        for k in 1:length(model.groups.groups[i])
            firstloop == false && print(io, ", ")
            print(io, "\"", model.groups.groups[i][k], "\" => ", model.groups.n_groups[i][k])
            firstloop = false
        end
        println(io)
    end
    print(io, "Contains parameters: ")
    firstloop = true
    for fieldname in fieldnames(typeof(model.params))
        firstloop == false && print(io, ", ")
        print(io, fieldname)
        firstloop = false
    end
end

export eosshow


text:
using DynamicHMC,
    LinearAlgebra,
    LogDensityProblems,
    MCMCDiagnosticTools,
    Optim,
    Pathfinder,
    Random,
    Statistics,
    StatsFuns,
    Test,
    TransformVariables

Random.seed!(1)

struct RegressionProblem{X,Y}
    x::X
    y::Y
end

function (prob::RegressionProblem)(θ)
    σ = θ.σ
    α = θ.α
    β = θ.β
    x = prob.x
    y = prob.y
    lp = normlogpdf(σ) + logtwo
    lp += normlogpdf(α)
    lp += sum(normlogpdf, β)
    y_hat = muladd(x, β, α)
    lp += sum(eachindex(y_hat, y)) do i
        return normlogpdf(y_hat[i], σ, y[i])
    end
    return lp
end

function mean_and_mcse(f, θs)
    zs = map(f, θs)
    ms = mean(zs)
    ses = map(mcse, eachrow(reduce(hcat, zs)))
    return ms, ses
end

function compare_estimates(f, xs1, xs2, α=0.05)
    nparams = length(first(xs1))
    α /= nparams  # bonferroni correction
    p = α / 2
    m1, s1 = mean_and_mcse(f, xs1)
    m2, s2 = mean_and_mcse(f, xs2)
    zs = @. (m1 - m2) / sqrt(s1^2 + s2^2)
    @test all(norminvcdf(p) .< zs .< norminvccdf(p))
end

@testset "DynamicHMC integration" begin
    A = Diagonal(rand(5))
    B = randn(5, 2)
    D = exp(Symmetric(randn(2, 2)))
    M⁻¹ = Pathfinder.WoodburyPDMat(A, B, D)

    @testset "Pathfinder.WoodburyLeftInvFactor" begin
        L = Pathfinder.WoodburyLeftInvFactor(M⁻¹)
        v = randn(5)
        V = randn(5, 2)
        @test L.A === M⁻¹
        @test size(L) === size(M⁻¹)
        Lmat = @inferred Matrix(L)
        @test Lmat isa Matrix
        @test Lmat * Lmat' ≈ inv(M⁻¹)
        @test L * v ≈ Lmat * v
        @test L * V ≈ Lmat * V
        @test AbstractMatrix{Float32}(L) isa Pathfinder.WoodburyLeftInvFactor{Float32}
        for i in axes(L, 1), j in axes(L, 2)
            @test L[i, j] ≈ Lmat[i, j]
        end
    end

    @testset "DynamicHMC.GaussianKineticEnergy" begin
        κ = DynamicHMC.GaussianKineticEnergy(M⁻¹)
        @test κ.M⁻¹ === M⁻¹
        @test κ.W === Pathfinder.WoodburyLeftInvFactor(M⁻¹)
        κdense = DynamicHMC.GaussianKineticEnergy(Symmetric(Matrix(M⁻¹)))
        κdense2 = DynamicHMC.GaussianKineticEnergy(Symmetric(Matrix(M⁻¹)), Matrix(κ.W))
        p = randn(5)
        q = randn(5)

        @test size(κ) == size(κdense)
        @test DynamicHMC.kinetic_energy(κ, p, q) ≈ DynamicHMC.kinetic_energy(κdense, p, q)
        @test DynamicHMC.calculate_p♯(κ, p, q) ≈ DynamicHMC.calculate_p♯(κdense, p, q)
        @test DynamicHMC.∇kinetic_energy(κ, p, q) ≈ DynamicHMC.∇kinetic_energy(κdense, p, q)
        @test DynamicHMC.rand_p(MersenneTwister(42), κ, q) ≈
            DynamicHMC.rand_p(MersenneTwister(42), κdense2, q)
    end

    @testset "DynamicHMC.mcmc_with_warmup" begin
        ndraws = 1_000
        x = 0:0.01:1
        y = sin.(x) .+ randn.() .* 0.2 .+ x
        X = [x x .^ 2 x .^ 3]
        prob = RegressionProblem(X, y)
        trans = as((σ=asℝ₊, α=asℝ, β=as(Array, size(X, 2))))
        P = TransformedLogDensity(trans, prob)
        ∇P = ADgradient(:ForwardDiff, P)
        rng = Random.GLOBAL_RNG

        result_hmc1 = mcmc_with_warmup(rng, ∇P, ndraws; reporter=NoProgressReport())

        logp(x) = LogDensityProblems.logdensity(P, x)
        ∇logp(x) = LogDensityProblems.logdensity_and_gradient(∇P, x)[2]
        result_pf = pathfinder(logp, ∇logp; dim=LogDensityProblems.dimension(P))

        @testset "Initial point" begin
            result_hmc2 = mcmc_with_warmup(
                rng,
                ∇P,
                ndraws;
                initialization=(; q=result_pf.draws[:, 1]),
                reporter=NoProgressReport(),
            )
            @test result_hmc2.κ.M⁻¹ isa Diagonal
            compare_estimates(identity, result_hmc2.chain, result_hmc1.chain)
        end

        @testset "Initial point and metric" begin
            result_hmc3 = mcmc_with_warmup(
                rng,
                ∇P,
                ndraws;
                initialization=(;
                    q=result_pf.draws[:, 1],
                    κ=GaussianKineticEnergy(result_pf.fit_distribution.Σ),
                ),
                warmup_stages=default_warmup_stages(; M=Symmetric),
                reporter=NoProgressReport(),
            )
            @test result_hmc3.κ.M⁻¹ isa Symmetric
            compare_estimates(identity, result_hmc3.chain, result_hmc1.chain)
        end

        @testset "Initial point and final metric" begin
            result_hmc4 = mcmc_with_warmup(
                rng,
                ∇P,
                ndraws;
                initialization=(;
                    q=result_pf.draws[:, 1],
                    κ=GaussianKineticEnergy(result_pf.fit_distribution.Σ),
                ),
                warmup_stages=default_warmup_stages(; middle_steps=0, doubling_stages=0),
                reporter=NoProgressReport(),
            )
            @test result_hmc4.κ.M⁻¹ === result_pf.fit_distribution.Σ
            compare_estimates(identity, result_hmc4.chain, result_hmc1.chain)
        end
    end
end


text:
print("Loading modules...")

print("Util...")
using Util
print("Params...")
using Params
print("DataFrames...")
using DataFrames
print("DataStructures...")
using DataStructures: DefaultDict
print("JLD...")
using HDF5, JLD
print("Gadfly...")
using Gadfly

println("DONE")

#=
Extracts and categorizes the data from saved jld files.

Returns a DataFrame containing all sufficient statistics
=#
function getdata()
    df = DataFrame()
    df[:goodratio] = Float64[]
    df[:pkt_max] = Int[]
    df[:n_agent] = Int[]
    df[:agent_type] = String[]
    df[:generated_packets] = Float64[]
    df[:tried_packets] = Float64[]
    df[:sent_packets] = Float64[]
    df[:buf_levels] = Matrix{Float64}[]
    df[:buf_overflow] = Float64[]
    df[:en] = Vector{Float64}[]
    df[:bits] = Vector{Float64}[]
    df[:prefix] = String[]
    df[:t] = Vector{Float64}[]
    df[:latency] = Matrix{Float64}[]
    df[:latencyhist] = Matrix{Dict{Int, Int}}[]
    df[:en_idle] = Float64[]
    df[:en_sw] = Float64[]
    df[:en_sense] = Float64[]
    # df[:en_sleep] = Float64[]
    df[:en_tx] = Float64[]
    df[:delta] = Float64[]
    df[:b_idle] = Float64[]
    df[:d_svd] = Int[]
    df[:trustQ] = Float64[]
    df[:lost_packets] = Float64[]
    df[:loss_rate] = Float64[]
    agent_types = ["CooperativeQ", "OptHighestSNR", "RandomChannel", "ContextQ"]
    for P = genparams()
        dir = joinpath("data/", P.prefix)
        if ! isdir(dir)
            continue
        end
        println("Processing ", P.prefix)
        # s=DataFrame(agent_type=ASCIIString[], buf=Float64[], en=Float64[], t=Int64[], bits=Float64[], ee=Float64[], cumee=Float64[], cumth=Float64[])
        for at in agent_types
            file = jldopen(joinpath(dir, "$at.jld"), "r")
            energy = vec(mean(read(file, "avg_energies"), [1]))
            bits = read(file, "avg_bits")
            buf_overflow = read(file, "buf_overflows")
            buf_levels = read(file, "buf_levels")
            generated_packets = mean(read(file, "generated_packets"))
            tried_packets = mean(read(file, "tried_packets"))
            sent_packets = mean(read(file, "sent_packets"))
            ee = mean(bits)/mean(energy)
            latency = read(file, "latencies")
            extra = jldopen(joinpath(dir, "$at-extra.jld"), "r")
            latencyhist = read(extra, "latencyhist")
            energies = jldopen(joinpath(dir, "$at-energies.jld"), "r")
            en_total = mean(energy)
            en_idle = mean(vec(mean(read(energies, "en_idle"), [1]))) / en_total
            en_sw = mean(vec(mean(read(energies, "en_sw"), [1]))) / en_total
            en_sense = mean(vec(mean(read(energies, "en_sense"), [1]))) / en_total
            # en_sleep = mean(vec(mean(read(energies, "en_sleep"), [1]))) / en_total
            en_tx = mean(vec(mean(read(energies, "en_tx"), [1]))) / en_total
            lostinchan = read(file, "lostinchan")
            close(energies)
            d = {
                (P.goodratio[1] / P.goodratio[2])
                P.pkt_max
                P.n_agent
                at
                generated_packets
                tried_packets
                sent_packets
                buf_levels
                buf_overflow
                energy
                vec(mean(bits, [1]))
                P.prefix
                [1:Params.t_total]
                latency
                latencyhist
                en_idle
                en_sw
                en_sense
                #en_sleep
                en_tx
                P.δ
                P.beta_idle
                P.d_svd
                P.trustQ
                lostinchan
                lostinchan / tried_packets
            }
            push!(df, d)
            b=vec(mean(bits, [1]))
            en=vec(mean(energy, [1]))
            close(file)
            avg_buf_overflow = mean(buf_overflow)
            avg_buf_levels = mean(buf_levels)
            println("$at\t$ee\t$avg_buf_levels\t$avg_buf_overflow")

        end
        prefix = P.prefix
        #= s[:t] .*= 10e-3
        theme = Theme(major_label_font_size=6mm, minor_label_font_size=4mm, key_title_font_size=6mm, key_label_font_size=4mm, line_width=1mm, default_point_size=0.8mm)
        yticks = {:th => Guide.yticks(ticks=[1000:1000:7000]), :ee => Guide.yticks(ticks=[3e6:1e6:6e6])}
        draw(PDF("plots/$prefix-ee-cumulative.pdf", 6inch, 4inch), plot(s[5:end,:], x="t", y="cumee", color="agent_type", Geom.smooth, yticks[:ee], Guide.xlabel("Time (s)"), Guide.ylabel("EE (b/J)"), Guide.colorkey("Algorithm"), theme))
        draw(PDF("plots/$prefix-ee-smoothed.pdf", 6inch, 4inch), plot(s, x="t", y="ee", color="agent_type", Geom.smooth, Guide.xlabel("Time (s)"), Guide.ylabel("EE (b/J)"), Guide.colorkey("Algorithm"), theme))
        draw(PDF("plots/$prefix-throughput.pdf", 6inch, 4inch), plot(s[1:100:end, :], x="t", y="bits", color="agent_type", Geom.smooth, yticks[:th], Guide.xlabel("Time (s)"), Guide.ylabel("Throughput (packets)"), Guide.colorkey("Algorithm"), theme))
        draw(PDF("plots/$prefix-buffer.pdf", 6inch, 4inch), plot(s[1:100:end, :], x="t", y="buf", color="agent_type", Geom.smooth, Guide.xlabel("Time (s)"), Guide.ylabel("Buffer Occupancy (packets)"), Guide.colorkey("Algorithm"), theme)) =#
    end
    df[:en_mean] = map(mean, df[:en])
    df[:en_std] = map(std, df[:en])
    df[:bits_mean] = map(mean, df[:bits])
    df[:bits_std] = map(std, df[:bits])
    ee = map(./, df[:bits], df[:en])
    df[:ee] = ee
    df[:ee_cum] = map(./, map(cumsum, df[:bits]), map(cumsum, df[:en]))
    df[:ee_mean] =  map(mean, ee)
    # retries per sent packet
    df[:retries] = (df[:tried_packets] ./ df[:sent_packets]) - 1
    df[:buf_mean] = map(mean, df[:buf_levels])
    df[:buf_max] = map(maximum, df[:buf_levels])
    df[:buf_std] = map(std, df[:buf_levels])
    df[:latency_mean] = map(mean, df[:latency])
    df[:latency_std] = map(std, df[:latency])
    df[:latency_max] = map(maximum, df[:latency])
    df
end

#=
Save sufficient statistics
=#
function savestats(df)
    stats = copy(df)
    for key in [:en, :bits, :t, :buf_levels, :latency, :latencyhist, :ee, :ee_cum]
        delete!(stats, key)
    end
    # save whole statistics
    writetable("suffstats.csv", stats)
end

function plotframes(df)
    theme = Theme(major_label_font_size=6mm, minor_label_font_size=4mm, key_title_font_size=6mm, key_label_font_size=4mm, line_width=1mm, default_point_size=0.8mm)
    yticks = {:th => Guide.yticks(ticks=[1000:1000:7000]), :ee => Guide.yticks(ticks=[3e6:1e6:6e6])}
    for data in groupby(df, :prefix)
        prefix = data[1, :prefix]
        println("Plotting $prefix")
        ∇ = DataFrame(agent_type=ASCIIString[], ee=Float64[], ee_cum=Float64[], bits=Float64[], t=Float64[], buf_levels=Float64[])
        for i=1:size(data)[1]
            dd = DataFrameRow(data, i)
            append!(∇, DataFrame(agent_type=fill(dd[:agent_type], length(dd[:ee])), ee=dd[:ee], ee_cum=dd[:ee_cum], bits=dd[:bits], t=dd[:t], buf_levels=vec(mean(dd[:buf_levels], 1))))
        end
        function sumlatencies(latencies)
            lat = DefaultDict(Int, Int, 0)
            for l in latencies
                for (k, v) in l
                    lat[k] += v
                end
            end
            highest = 0
            latarray = zeros(Float64, sum(map(i -> i[2], lat)))
            i = 1
            @inbounds for (k, v) in lat
                latarray[i:i+v-1] = k
                i += v
                highest = max(highest, v)
            end
            highest, latarray
        end
        highest = 0
        d = by(data, :agent_type) do df
            δ=DataFrame()
            δ[:agent_type] = ASCIIString[]
            δ[:latency] = Float64[]
            for i in size(df)[1]
                maxi, l = sumlatencies(df[i, :latencyhist])
                highest = max(highest, maxi)
                Δ = DataFrame(agent_type=fill(df[i, :agent_type], length(l)), latency=l)
                append!(δ, Δ)
            end
            δ
        end
        println("Plotting Latency histogram")
        # Plot latency histograms
        latencyhist = plot(d, x="latency", color="agent_type", Geom.histogram(position=:dodge, bincount=10), Guide.xlabel("Latency (time slots)"), Guide.colorkey("Algorithm"), theme, Guide.yticks(ticks=linspace(0, highest * 2, 5)))
        draw(PDF("plots/$prefix-latency-histogram.pdf", 10inch, 6inch), latencyhist)
        smooth = Geom.smooth(smoothing=0.2)
        println("Plotting other plots")
        draw(PDF("plots/$prefix-ee-cumulative.pdf", 6inch, 4inch), plot(∇, x="t", y="ee_cum", color="agent_type", smooth, yticks[:ee], Guide.xlabel("Time (s)"), Guide.ylabel("EE (b/J)"), Guide.colorkey("Algorithm"), theme))
        draw(PDF("plots/$prefix-ee-smoothed.pdf", 6inch, 4inch), plot(∇, x="t", y="ee", color="agent_type", smooth, Guide.xlabel("Time (s)"), Guide.ylabel("EE (b/J)"), Guide.colorkey("Algorithm"), theme))
        draw(PDF("plots/$prefix-throughput.pdf", 6inch, 4inch), plot(∇, x="t", y="bits", color="agent_type", smooth, yticks[:th], Guide.xlabel("Time (s)"), Guide.ylabel("Throughput (packets)"), Guide.colorkey("Algorithm"), theme))
        draw(PDF("plots/$prefix-buffer.pdf", 6inch, 4inch), plot(∇, x="t", y="buf_levels", color="agent_type", smooth, Guide.xlabel("Time (s)"), Guide.ylabel("Buffer Occupancy (packets)"), Guide.colorkey("Algorithm"), theme))
    end
end

df = getdata()
println("Saving statistics")
savestats(df)
# plotframes(df)


text:
Pkg.add("LLLplus")

#include("../src/MUMIMO.jl")
using MUMIMO
using Plots

mimoUplink(10000,4,[4],[4],[0],[12],12,[-5.0:5:20;])
savefig("benchmark/perfVsSNRqpsk4ant.png") # run from root directory
mimoUplink(10000,4,[8],[8],[0],[12],12,[-5.0:5:20;])
savefig("benchmark/perfVsSNRqpsk8ant.png") # run from root directory

max_stars_repo_name: christianpeel/MUMIMO

max_issues_repo_name: christianpeel/MUMIMO

max_forks_repo_name: christianpeel/MUMIMO


text:
module FluxUtils

using FluxRM

function __get_treedict!(dict, path...)
    next_dict = dict
    for node in path
        next_dict = get!(()-> Dict{String, Any}(), next_dict, node)
    end
    return next_dict
end

function __set_io_path(spec::FluxRM.JobSpec.Jobspec, iotype, name, path)
    system = spec.attributes.system
    if system.shell === nothing
        system.shell = Dict{String, Any}()
    end
    io = __get_treedict!(system.shell, "options", iotype, name)

    io["type"] = "file"
    io["path"] = path
end

function propagate_load_path!(env, project=Base.ACTIVE_PROJECT[])
    pathsep = Sys.iswindows() ? ";" : ":"
    if get(env, "JULIA_LOAD_PATH", nothing) === nothing
        env["JULIA_LOAD_PATH"] = join(LOAD_PATH, pathsep)
    end
    if get(env, "JULIA_DEPOT_PATH", nothing) === nothing
        env["JULIA_DEPOT_PATH"] = join(DEPOT_PATH, pathsep)
    end
    # Set the active project on workers using JULIA_PROJECT.
    # Users can opt-out of this by (i) passing `env = ...` or (ii) passing
    # `--project=...` as `exeflags` to addprocs(...).
    if project !== nothing && get(env, "JULIA_PROJECT", nothing) === nothing
        env["JULIA_PROJECT"] = project
    end
end

function juliaspec(args, dir; num_nodes=1, num_tasks_per_node=1, cores_per_task=4, project=dir)
    num_tasks = num_nodes*num_tasks_per_node
    cmd = `$(Base.julia_cmd()) $(args)`
    jobspec = FluxRM.JobSpec.from_command(cmd; num_nodes, num_tasks, cores_per_task)
    system = jobspec.attributes.system
    system.cwd = realpath(dir)
    env = Dict{String, String}()
    propagate_load_path!(env, realpath(project))
    system.environment = env
    __set_io_path(jobspec, "output", "stderr", "flux-{{id}}.err")
    __set_io_path(jobspec, "output", "stdout", "flux-{{id}}.out")
    @assert FluxRM.JobSpec.validate(jobspec, 1)
    jobspec
end

function nodes()
    rpc = fetch(FluxRM.RPC(Flux(), "resource.status", nodeid=0))
    R = first(rpc.R.execution.R_lite)

    hosts = FluxRM.IDSet(R.rank)
    Int(length(hosts))
end

struct Job
    workdir::String
    preload::String
    cmd::String
end

function Base.run(job::Job)
    args = String[]
    if job.preload != ""
        append!(args, ["-L", job.preload])
    end
    append!(args, ["-e", job.cmd])

    jobspec = juliaspec(args, job.workdir)
    submission = FluxRM.submit(Flux(), jobspec)
    return FluxRM.Job(submission)
end

end # module


text:
module Routing # Bukdu

using ..Bukdu: ApplicationController, Route, Conn
using ..Bukdu.Deps
using ..Bukdu.Naming
using ..Bukdu.System
using ..Bukdu.Plug

store = Dict{Symbol, Any}(
    :pipe => nothing,
    :routing_tables => Vector{Any}(),
    :routing_path => Dict{Tuple{Symbol,Symbol,Symbol},String}(),
)
routing_pipelines = Dict{Symbol, Vector{Function}}()

function handle_conn(req::Deps.Request, method::String)
    uri = Deps.HTTP.URIs.URI(req.target)
    segments = split(uri.path, '/'; keepempty=false)
    vals = [Val(Symbol(seg)) for seg in segments]
    verb = Val(Symbol(method))
    route(verb, vals...)
end

struct AbstractControllerError <: Exception
    msg
end

(::Type{Vector{Pair{String,Any}}})(p::Pair{String,T}...) where T = [p...]

route(args...) = Route(System.MissingController, System.not_found, Dict{Symbol,DataType}(), Vector{Pair{String,Any}}(), Vector{Function}())

# idea from HTTP.jl/src/Handlers.jl
function penetrate_segments(segments)
    vals = Expr[]
    path_params = Expr[]
    for seg in segments
        if startswith(seg, ':')
            param = seg[2:end]
            mangled = Symbol(param, :_)
            pair = :(Pair(String($param), String(first(typeof($mangled).parameters))))
            push!(path_params, pair)
            expr = :($mangled::Any)
        else
            expr = :(::Val{Symbol($seg)})
        end
        push!(vals, expr)
    end
    return (vals, path_params)
end

function add_route(verb, url::String, C::Type{<:ApplicationController}, action, param_types::Dict{Symbol,DataType}) # throw AbstractControllerError
    isabstracttype(C) && throw(AbstractControllerError(string("use concrete type")))
    pipe = store[:pipe]
    pipelines = get(routing_pipelines, pipe, [])
    segments = split(url, '/'; keepempty=false)
    (vals, path_params) = penetrate_segments(segments)
    method = Naming.verb_name(verb)
    @eval route(::Val{Symbol($method)}, $(vals...)) = Route($C, $action, $param_types, Vector{Pair{String,Any}}($(path_params...)), $pipelines)
    routing_tables = vcat(store[:routing_tables],
        Naming.verb_name(verb), url, nameof(C), nameof(action), (pipe isa Nothing ? "" : repr(pipe)))
    store[:routing_tables] = routing_tables
    store[:routing_path][Naming.routing_path_key(verb,C,action)] = url
end

function route_path(verb, C::Type{<:ApplicationController}, action)::Union{String,Nothing}
    get(store[:routing_path], Naming.routing_path_key(verb,C,action), nothing)
end

"""
    Routing.empty!()
"""
function empty!()
    store[:pipe] = nothing
    store[:routing_tables] = Vector{Any}()
    store[:routing_path] = Dict{Tuple{Symbol,Symbol,Symbol},String}()
    Base.empty!(routing_pipelines)
    Plug.ContentParsers.env[:decoders] = Plug.ContentParsers.default_content_decoders
    Plug.ContentParsers.env[:parsers] = Plug.ContentParsers.default_content_parsers
    ms = methods(route)
    for m in ms
        Base.delete_method(m)
    end
    @eval route(args...) = Route(System.MissingController, System.not_found, Dict{Symbol,DataType}(), Vector{Pair{String,Any}}(), Vector{Function}())
end

end # module Bukdu.Routing


text:
include("stochasticsimulation.jl")
include("quadrature.jl")
using Plots

# Demere
num = 100
for i in 1:10
    prob = demere(num)
    print(prob)
    print("\t")
end
print("\n")

# Demere
num = 100000
for i in 1:10
    prob1 = demere(num)
    print(prob1)
    print("\t")
end
print("\n")

# roulette example
prob2 = roulette(8, 23, 10000)
println(prob2)

# hit and miss monte carlo method
f(x) = exp(-x)
val = hitmiss(f)
println(val)

# crude monte carlo method
val = crudemc(f)
println(val)

# compare with standard method
N = 10 .^ (3:6)
fexact = 1.0 - exp(-1.0)
hmerr = []
mcerr = []
trerr = []
box = [0;1]
for i in N
    append!(hmerr,abs(hitmiss(f, box, i)-fexact))
    append!(mcerr,abs(crudemc(f, box, i)-fexact))
    append!(trerr,abs(ctr(f, box, i)-fexact))
end

# Plot the errors
pgfplotsx();
plot(N, hmerr, xaxis=:log, yaxis=:log)
plot!(N, mcerr, xaxis=:log, yaxis=:log)
plot!(N, trerr, xaxis=:log, yaxis=:log,linestyle = :dot,linealpha = 0.5, linewidth = 4, linecolor=:red)


text:
# Plotting

function layerPolygon(arr::Array{Point{T, 2}, 1}) where T
    xx = []
    yy = []
    for i in 1:length(arr)
        push!(xx, arr[i].coordinates[1])
        push!(yy, arr[i].coordinates[2])
    end
    return layer(x = xx, y = yy, Geom.polygon(preserve_order = true, fill = true))
end

function layer(arr::Array{Point{T, 2}, 1}) where T
    xx = []
    yy = []
    for i in 1:length(arr)
        push!(xx, arr[i].coordinates[1])
        push!(yy, arr[i].coordinates[2])
    end
    return layer(x = xx, y = yy)
end


text:
# Autogenerated wrapper script for OpenSSH_jll for aarch64-linux-musl
export scp, ssh, ssh_add, ssh_agent, ssh_keygen, ssh_keyscan

using Zlib_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("OpenSSH")
JLLWrappers.@declare_executable_product(scp)
JLLWrappers.@declare_executable_product(ssh)
JLLWrappers.@declare_executable_product(ssh_add)
JLLWrappers.@declare_executable_product(ssh_agent)
JLLWrappers.@declare_executable_product(ssh_keygen)
JLLWrappers.@declare_executable_product(ssh_keyscan)
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, OpenSSL_jll)
    JLLWrappers.@init_executable_product(
        scp,
        "bin/scp",
    )

    JLLWrappers.@init_executable_product(
        ssh,
        "bin/ssh",
    )

    JLLWrappers.@init_executable_product(
        ssh_add,
        "bin/ssh-add",
    )

    JLLWrappers.@init_executable_product(
        ssh_agent,
        "bin/ssh-agent",
    )

    JLLWrappers.@init_executable_product(
        ssh_keygen,
        "bin/ssh-keygen",
    )

    JLLWrappers.@init_executable_product(
        ssh_keyscan,
        "bin/ssh-keyscan",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using SimpleDiffEq, StaticArrays, OrdinaryDiffEq, Test

function loop(u, p, t)
    @inbounds begin
        σ = p[1]; ρ = p[2]; β = p[3]
        du1 = σ*(u[2]-u[1])
        du2 = u[1]*(ρ-u[3]) - u[2]
        du3 = u[1]*u[2] - β*u[3]
        return SVector{3}(du1, du2, du3)
    end
end
function liip(du, u, p, t)
    σ = p[1]; ρ = p[2]; β = p[3]
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
    return nothing
end

u0 = 10ones(3)
dt = 1e-2

odeoop = ODEProblem{false}(loop, SVector{3}(u0), (0.0, 100.0),  [10, 28, 8/3])
odeiip = ODEProblem{true}(liip, u0, (0.0, 100.0),  [10, 28, 8/3])

oop = init(odeoop,SimpleATsit5(),dt=dt)
step!(oop); step!(oop)
iip = init(odeiip,SimpleATsit5(),dt=dt)
step!(iip); step!(iip)

deoop = DiffEqBase.init(odeoop, Tsit5(); dt = dt)
step!(deoop); step!(deoop)
@test oop.u ≈ deoop.u atol=1e-5
@test oop.t ≈ deoop.t atol=1e-5

deiip = DiffEqBase.init(odeiip, Tsit5(); dt = dt)
step!(deiip); step!(deiip)
@test iip.u ≈ deiip.u atol=1e-5
@test iip.t ≈ deiip.t atol=1e-5

sol = solve(odeoop,SimpleATsit5(),dt=dt)

# Test keywords:
oop = init(odeoop,SimpleATsit5(),dt=dt, reltol = 1e-9, abstol = 1e-9)
step!(oop); step!(oop)
deoop = DiffEqBase.init(odeoop, Tsit5(); dt = dt, reltol=1e-9, abstol=1e-9)
step!(deoop); step!(deoop)

@test oop.u ≈ deoop.u atol=1e-5
@test oop.t ≈ deoop.t atol=1e-5

# Test reinit!
reinit!(oop, odeoop.u0; dt = dt)
reinit!(iip, odeiip.u0; dt = dt)
step!(oop); step!(oop)
step!(iip); step!(iip)

@test oop.u ≈ deoop.u atol=1e-5
@test oop.t ≈ deoop.t atol=1e-5
@test iip.u ≈ deiip.u atol=1e-5
@test iip.t ≈ deiip.t atol=1e-5

# Interpolation tests
uprev = copy(oop.u)
step!(oop)
@test uprev ≈ oop(oop.tprev) atol = 1e-12
@test oop(oop.t) ≈ oop.u atol = 1e-12

uprev = copy(iip.u)
step!(iip)
@test uprev ≈ iip(iip.tprev) atol = 1e-12
@test iip(iip.t) ≈ iip.u atol = 1e-12

# Interpolation tests comparing Tsit5 and SimpleATsit5
function f(du,u,p,t)
    du[1] = 2.0*u[1] + 3.0*u[2]
    du[2] = 4.0*u[1] + 5.0*u[2]
end
tmp = [1.0;1.0]
tspan = (0.0,1.0)
prob = ODEProblem(f,tmp,tspan)
integ1 = init(prob, SimpleATsit5(), abstol = 1e-6, reltol = 1e-6, save_everystep = false, dt = 0.1)
integ2 = init(prob, Tsit5(), abstol = 1e-6, reltol = 1e-6, save_everystep = false, dt = 0.1)
step!(integ2)
step!(integ1)
for i in 1:9
    x = i/10
    y = 1 - x
    @test integ1(x * integ2.t + y * integ2.tprev) ≈ integ2(x * integ2.t + y * integ2.tprev) atol=1e-7
end
step!(integ2)
step!(integ1)
for i in 1:9
    x = i/10
    y = 1 - x
    @test integ1(x * integ2.t + y * integ2.tprev) ≈ integ2(x * integ2.t + y * integ2.tprev) atol=1e-7
end
###################################################################################
# Internal norm test:
function moop(u, p, t)
    x = loop(u[:, 1], p, t)
    y = loop(u[:, 2], p, t)
    return hcat(x,y)
end
function miip(du, u, p, t)
    @views begin
        liip(du[:, 1], u[:, 1], p, t)
        liip(du[:, 2], u[:, 2], p, t)
    end
    return nothing
end

ran = rand(SVector{3})
odemoop = ODEProblem{false}(moop, SMatrix{3,2}(hcat(u0, ran)), (0.0, 100.0),  [10, 28, 8/3])
odemiip = ODEProblem{true}(miip, hcat(u0, ran), (0.0, 100.0),  [10, 28, 8/3])

using LinearAlgebra

oop = init(odemoop,SimpleATsit5(),dt=dt, internalnorm = (u, t) -> norm(u[:, 1]))
step!(oop); step!(oop)

iip = init(odemiip,SimpleATsit5(),dt=dt, internalnorm = (u, t) -> norm(u[:, 1]))
step!(iip); step!(iip)

@test oop.u ≈ iip.u atol=1e-5
@test oop.t ≈ iip.t atol=1e-5

###################################################################################
# VectorVector test:
function vvoop(du, u, p, t) # takes Vector{SVector}
    @inbounds for j in 1:2
        du[j] = loop(u[j], p, t)
    end
    return nothing
end
function vviip(du, u, p, t) # takes Vector{Vector}
    @inbounds for j in 1:2
        liip(du[j], u[j], p, t)
    end
    return nothing
end

ran = rand(3)
odevoop = ODEProblem{true}(vvoop, [SVector{3}(u0),  SVector{3}(ran)], (0.0, 100.0),  [10, 28, 8/3])
odeviip = ODEProblem{true}(vviip, [u0, ran], (0.0, 100.0),  [10, 28, 8/3])

viip = init(odeviip,SimpleATsit5(),dt=dt; internalnorm = (u, t) -> DiffEqBase.ODE_DEFAULT_NORM(u[1], t))
step!(viip); step!(viip)

iip = init(odeiip,SimpleATsit5(),dt=dt)
step!(iip); step!(iip)

@test iip.u ≈ viip.u[1] atol=1e-5

voop = init(odevoop,SimpleATsit5(),dt=dt,internalnorm = (u, t) -> DiffEqBase.ODE_DEFAULT_NORM(u[1], t))
step!(voop); step!(voop)

oop = init(odeoop,SimpleATsit5(),dt=dt)
step!(oop); step!(oop)

@test oop.u ≈ voop.u[1] atol=1e-5

# Final test that the states of both methods should be the same:
@test voop.u[2] ≈ viip.u[2] atol=1e-5

# viip = init(odeviip,SimpleATsit5(),dt=dt; internalnorm = u -> SimpleDiffEq.defaultnorm(u[1]))
# step!(viip); step!(viip)
# x = Float64[]; y = Float64[]
# for i in 1:1000
#     step!(viip)
#     push!(x, viip.u[1][1])
#     push!(y, viip.u[1][2])
# end
# using PyPlot
# plot(x,y)

###################################################################################
# Issue #46 check (step!(integ, dt, true))
# see https://github.com/SciML/SimpleDiffEq.jl/issues/46
using SimpleDiffEq, OrdinaryDiffEq, StaticArrays

@inbounds function duffing_rule(x, p, t)
    ω, f, d, β = p
    dx1 = x[2]
    dx2 = f*cos(ω*t) - β*x[1] - x[1]^3 - d * x[2]
    return SVector(dx1, dx2)
end
prob = ODEProblem(duffing_rule, SVector(0.1, 0.2), (0.0, 1e12), [0.3, 0.1, 0.2, -1])

T = 2π/0.3 # this is the period of the oscillator

integ2 = init(prob, SimpleATsit5(), reltol=1e-9)
step!(integ2, T*20, true)
v=zeros(200,2)
for k in 1:200
   step!(integ2, T, true)
   v[k,:] = integ2.u
end

@test length(unique(round.(v; digits = 3))) == 2


text:
#=
   Dice puzzle in Turing.jl

   From Berlin Bayesians:
   (https://app.slack.com/client/TFPMSKW3F/CFQHMRD6K)
   """
   What's likelier: Rolling at least one six in four throws of a single die,
   or rolling at least one double 6 in 24 throws of a pair of dice?
   """

   This is de Méré's classical dice puzzle (which induced the study of
   probability theory). See
   - https://mathworld.wolfram.com/deMeresProblem.html

   * 4 throws getting at least one 6_
     1-(5/6)**4  ~ 0.51774691358024691358

   * 24 throws getting at least one double 6
     1-(35/36)**24 ~ 0.49140387613090325958

   Puzzle 1:
   Distributions of variable p (num:0)
   1.00000 =>    5162  (0.516200)
   0.00000 =>    4838  (0.483800)

   Distributions of variable s (num:0)
   0.00000 =>    4838  (0.483800)
   1.00000 =>    3803  (0.380300)
   2.00000 =>    1183  (0.118300)
   3.00000 =>     173  (0.017300)
   4.00000 =>       3  (0.000300)

   Puzzle 2:
   Distributions of variable p (num:0)
   0.00000 =>    5014  (0.501400)
   1.00000 =>    4986  (0.498600)

   Distributions of variable s (num:0)
   0.00000 =>    5014  (0.501400)
   1.00000 =>    3571  (0.357100)
   2.00000 =>    1129  (0.112900)
   3.00000 =>     248  (0.024800)
   4.00000 =>      37  (0.003700)
   5.00000 =>       1  (0.000100)

   Cf ~/webppl/dice_puzzle.jl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model dice_puzzle1() = begin
    n = 4
    throws = tzeros(n)
    for t in 1:n
        throws[t] ~ DiscreteUniform(1,6)
    end
    s ~ Dirac(sum(throws[t] == 6 for t in 1:n))
    p ~ Dirac(s >= 1)
end

@model dice_puzzle2() = begin
    n = 24
    throws = tzeros(n,2)
    for t in 1:n
        throws[t,1] ~ DiscreteUniform(1,6)
        throws[t,2] ~ DiscreteUniform(1,6)
    end
    s ~ Dirac(sum(throws[t,1] == 6 && throws[t,2] == 6 for t in 1:n))
    p ~ Dirac(s >= 1)
end

function run_model(model)

    num_chns = 4

    # chns = sample(model(), Prior(), 10_000)
    # chns = sample(model(), MH(), 10_000)
    # chns = sample(model(), PG(15), 10_000)
    # chns = sample(model(), SMC(), 10_000)
    chns = sample(model(), IS(), 10_000)
    #

    # display(chns)
    # display(plot(chns))
    show_var_dist_pct(chns,:p)
    show_var_dist_pct(chns,:s)

end

println("Puzzle 1:")
@time run_model(dice_puzzle1)

println("\n\nPuzzle 2:")
@time run_model(dice_puzzle2)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
export AbstractFeatureConstructor, create_features, feature_size

"""
    AbstractFeatureCreator
An abstract feature creator, for feature transformation from the states.
"""
abstract type AbstractFeatureConstructor end

"""
    create_features
Actually create the features
"""
function create_features end

"""
    feature_size(fc::AbstractFeatureCreator)
Get size of feature vector the features assume exists.
"""
function feature_size(fc::AbstractFeatureConstructor)
    throw("Implement feature size for $(typof(fc))")
end


text:
import Oceananigans.TimeSteppers: calculate_tendencies!
import Oceananigans: tracer_tendency_kernel_function

using Oceananigans.Architectures: device_event
using Oceananigans: fields, prognostic_fields
using Oceananigans.Utils: work_layout

"""
    calculate_tendencies!(model::NonhydrostaticModel)

Calculate the interior and boundary contributions to tendency terms without the
contribution from non-hydrostatic pressure.
"""
function calculate_tendencies!(model::HydrostaticFreeSurfaceModel)

    # Calculate contributions to momentum and tracer tendencies from fluxes and volume terms in the
    # interior of the domain
    calculate_hydrostatic_free_surface_interior_tendency_contributions!(model)

    # Calculate contributions to momentum and tracer tendencies from user-prescribed fluxes across the
    # boundaries of the domain
    calculate_hydrostatic_boundary_tendency_contributions!(model.timestepper.Gⁿ,
                                                           model.grid,
                                                           model.architecture,
                                                           model.velocities,
                                                           model.free_surface,
                                                           model.tracers,
                                                           model.clock,
                                                           fields(model),
                                                           model.closure,
                                                           model.buoyancy)

    return nothing
end

function calculate_free_surface_tendency!(arch, grid, model, dependencies)

    Gη_event = launch!(arch, grid, :xy,
                       calculate_hydrostatic_free_surface_Gη!, model.timestepper.Gⁿ.η,
                       grid,
                       model.velocities,
                       model.free_surface,
                       model.tracers,
                       model.auxiliary_fields,
                       model.forcing,
                       model.clock;
                       dependencies = dependencies)

    return Gη_event
end

""" Calculate momentum tendencies if momentum is not prescribed. `velocities` argument eases dispatch on `PrescribedVelocityFields`."""
function calculate_hydrostatic_momentum_tendencies!(model, velocities; dependencies = device_event(model))

    arch = model.architecture
    grid = model.grid

    momentum_kernel_args = (grid,
                            model.advection.momentum,
                            model.coriolis,
                            model.closure,
                            velocities,
                            model.free_surface,
                            model.tracers,
                            model.buoyancy,
                            model.diffusivity_fields,
                            model.pressure.pHY′,
                            model.auxiliary_fields,
                            model.forcing,
                            model.clock)

    Gu_event = launch!(arch, grid, :xyz,
                       calculate_hydrostatic_free_surface_Gu!, model.timestepper.Gⁿ.u, momentum_kernel_args...;
                       dependencies = dependencies)

    Gv_event = launch!(arch, grid, :xyz,
                       calculate_hydrostatic_free_surface_Gv!, model.timestepper.Gⁿ.v, momentum_kernel_args...;
                       dependencies = dependencies)

    Gη_event = calculate_free_surface_tendency!(arch, grid, model, dependencies)

    events = [Gu_event, Gv_event, Gη_event]

    return events
end

using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVD, CATKEVDArray

# Fallback
@inline tracer_tendency_kernel_function(model::HydrostaticFreeSurfaceModel, closure, tracer_name) =
    hydrostatic_free_surface_tracer_tendency

@inline tracer_tendency_kernel_function(model::HydrostaticFreeSurfaceModel, closure::CATKEVD, ::Val{:e}) =
    hydrostatic_turbulent_kinetic_energy_tendency

function tracer_tendency_kernel_function(model::HydrostaticFreeSurfaceModel, closures::Tuple, ::Val{:e})
    if any(cl isa Union{CATKEVD, CATKEVDArray} for cl in closures)
        return hydrostatic_turbulent_kinetic_energy_tendency
    else
        return hydrostatic_free_surface_tracer_tendency
    end
end

top_tracer_boundary_conditions(grid, tracers) =
    NamedTuple(c => tracers[c].boundary_conditions.top for c in propertynames(tracers))

""" Store previous value of the source term and calculate current source term. """
function calculate_hydrostatic_free_surface_interior_tendency_contributions!(model)

    arch = model.architecture
    grid = model.grid

    barrier = device_event(model)

    events = calculate_hydrostatic_momentum_tendencies!(model, model.velocities; dependencies = barrier)

    top_tracer_bcs = top_tracer_boundary_conditions(grid, model.tracers)

    for (tracer_index, tracer_name) in enumerate(propertynames(model.tracers))
        @inbounds c_tendency = model.timestepper.Gⁿ[tracer_name]
        @inbounds c_advection = model.advection[tracer_name]
        @inbounds c_forcing = model.forcing[tracer_name]
        c_kernel_function = tracer_tendency_kernel_function(model, model.closure, Val(tracer_name))

        Gc_event = launch!(arch, grid, :xyz,
                           calculate_hydrostatic_free_surface_Gc!,
                           c_tendency,
                           c_kernel_function,
                           grid,
                           Val(tracer_index),
                           c_advection,
                           model.closure,
                           model.buoyancy,
                           model.velocities,
                           model.free_surface,
                           model.tracers,
                           top_tracer_bcs,
                           model.diffusivity_fields,
                           model.auxiliary_fields,
                           c_forcing,
                           model.clock;
                           dependencies = barrier)

        push!(events, Gc_event)
    end

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end

#####
##### Tendency calculators for u, v
#####

""" Calculate the right-hand-side of the u-velocity equation. """
@kernel function calculate_hydrostatic_free_surface_Gu!(Gu, grid, args...)
    i, j, k = @index(Global, NTuple)
    @inbounds Gu[i, j, k] = hydrostatic_free_surface_u_velocity_tendency(i, j, k, grid, args...)
end

""" Calculate the right-hand-side of the v-velocity equation. """
@kernel function calculate_hydrostatic_free_surface_Gv!(Gv, grid, args...)
    i, j, k = @index(Global, NTuple)
    @inbounds Gv[i, j, k] = hydrostatic_free_surface_v_velocity_tendency(i, j, k, grid, args...)
end

#####
##### Tendency calculators for tracers
#####

""" Calculate the right-hand-side of the tracer advection-diffusion equation. """
@kernel function calculate_hydrostatic_free_surface_Gc!(Gc, tendency_kernel_function, grid, args...)
    i, j, k = @index(Global, NTuple)
    @inbounds Gc[i, j, k] = tendency_kernel_function(i, j, k, grid, args...)
end

#####
##### Tendency calculators for an explicit free surface
#####

""" Calculate the right-hand-side of the free surface displacement (η) equation. """
@kernel function calculate_hydrostatic_free_surface_Gη!(Gη, grid, args...)
    i, j = @index(Global, NTuple)
    @inbounds Gη[i, j, 1] = free_surface_tendency(i, j, grid, args...)
end

#####
##### Boundary condributions to hydrostatic free surface model
#####

function apply_flux_bcs!(Gcⁿ, events, c, arch, barrier, args...)
    x_bcs_event = apply_x_bcs!(Gcⁿ, c, arch, barrier, args...)
    y_bcs_event = apply_y_bcs!(Gcⁿ, c, arch, barrier, args...)
    z_bcs_event = apply_z_bcs!(Gcⁿ, c, arch, barrier, args...)

    push!(events, x_bcs_event, y_bcs_event, z_bcs_event)

    return nothing
end

""" Apply boundary conditions by adding flux divergences to the right-hand-side. """
function calculate_hydrostatic_boundary_tendency_contributions!(Gⁿ, grid, arch, velocities, free_surface, tracers, args...)

    barrier = device_event(arch)

    events = []

    # Velocity fields
    for i in (:u, :v)
        apply_flux_bcs!(Gⁿ[i], events, velocities[i], arch, barrier, args...)
    end

    # Free surface
    apply_flux_bcs!(Gⁿ.η, events, displacement(free_surface), arch, barrier, args...)

    # Tracer fields
    for i in propertynames(tracers)
        apply_flux_bcs!(Gⁿ[i], events, tracers[i], arch, barrier, args...)
    end

    events = filter(e -> typeof(e) <: Event, events)

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end


text:
export gain, izero, pole

"""
### function gain(sys::StateSpace)

Compute the gain of a system in a State Space representation.

##### Args

* sys: State space.

##### Remarks

The idea was obtained from the routine TG04BX of Octave-Force project.

"""

function gain(sys::StateSpace)
    ( !issiso(sys) ) && error("gain() is only defined for SISO systems.")

    # To compute the gain, notice that:
    #
    #                                #zeros
    #                                 prod  (s-zi)
    #                  -1             i = 1
    #    G(s) = C.(sI-A) . B+D = K . --------------
    #                                #poles
    #                                 prod  (s-pi)
    #                                 i = 1
    #
    # Let s0 \in |R be a number that is neither a pole nor a zeros. Hence,
    #
    # #zeros                        #poles
    #  prod  (s0-zi) |= 0    and     prod  (s0-pi) != 0
    #  i = 1                         i = 1
    #
    # Finally, the system gain can be computed using:
    #
    #                        #poles
    #      _              _   prod  (s0-zi)
    #     |         -1     |  i = 1
    # K = |C.(s0.I-A) . B+D|.--------------
    #     |_              _| #zeros
    #                         prod  (s0-pi)
    #                         i = 1
    #

    # System poles.
    p = pole(sys)

    # System invariant zeros.
    z = izero(sys)

    # Find a real value for s0 that is neither a pole nor a zero.
    #
    # TODO: Improve this algorithm.
    #
    # This version was adapted from routine TG04BX of Octave-Forge project.
    s0 = 2*maximum(abs.([z;p]))

    # Compute:
    #
    #    #zeros
    #     prod  (s0-zi)
    #     i = 1
    #
    # Notice that this is a real number because s0 is real and if zi is complex
    # then its conjugate is also a zero, for all i.
    prod_zeros = prod(s0-z)

    # Compute:
    #
    #    #poles
    #     prod  (s0-pi)
    #     i = 1
    #
    # Notice that this is a real number because s0 is real and if pi is complex
    # then its conjugate is also a pole, for all i.
    prod_poles = prod(s0-p)

    # Compute:
    #              _              _
    #             |         -1     |
    #     G(s0) = |C.(s0.I-A) . B+D|
    #             |_              _|
    #
    A, B, C, D = sys.A, sys.B, sys.C, sys.D
    G_s0 = C*inv(s0*eye(A)-A)*B+D

    # Compute the gain.
    (G_s0*prod_poles/prod_zeros)[1]
end

"""
### function izero(sys::StateSpace)

Compute the invariant zeros of the system.

##### Args

* sys: Dynamic system in space state representation.

##### Returns

* The invariant zeros of the system.

"""

function izero(sys::StateSpace)
    A, B, C, D = sys.A, sys.B, sys.C, sys.D

    # Number of states.
    num_states = size(A,1)

    # Compute the invariant zeros using generalized eigenvalues.
    M = [A -B; C -D]

    N = zeros(M)
    N[1:num_states, 1:num_states] = eye(A)

    D,V = eig(M,N)

    # Remove the eigenvalues at +Inf and -Inf.
    D[ -Inf .< real(D) .< +Inf ]
end

"""
### function izero(G::TransferFunction)

Compute the invariant zeros of the system.

##### Args

* G: Transfer function of the system.

##### Returns

* The zeros of the transfer function.

"""

function izero(G::TransferFunction)
    # Compute the roots of the denominator.
    roots(G.num)
end

"""
### function pole(sys::StateSpace)

Compute the poles of the system.

##### Args

* sys: Dynamic system in state space representation.

##### Returns

* The poles of the system.

"""

function pole(sys::StateSpace)
    # Compute the eigenvalues of matrix A.
    eig(sys.A)[1]
end

"""
### function pole(G::TransferFunction)

Compute the poles of the system.

##### Args

* G: Transfer function of the system.

##### Returns

* The poles of the transfer function.

"""

function pole(G::TransferFunction)
    # Compute the roots of the denominator.
    roots(G.den)
end


text:
# sum of squares: (n*(n+1)*(2n+1))/6    https://trans4mind.com/personal_development/mathematics/series/sumNaturalSquares.htm

sumsSquare = n -> ((n*(n+1))/2)^2
sumOfSquares = n -> (n*(n+1)*(2n+1))/6

println(Int64(sumsSquare(100) - sumOfSquares(100)))

max_stars_repo_path: 006 - Sum Square Difference/006 - Sum_Square_Difference.jl

max_stars_repo_name: LanguageAdventures/Project-Euler-Julia

max_issues_repo_name: LanguageAdventures/Project-Euler-Julia

max_forks_repo_name: LanguageAdventures/Project-Euler-Julia


text:
function run_target(M,net_param,target_type)

N, p, g, b, dt, T = net_param;
N = Int(N);

t   = collect(dt:dt:T);

#---------- 1. periodic functions ----------#
if target_type == "periodic"
    u_periodic  = zeros(length(t),N);
    for j=1:N
        # A  = 0.5 + 1.5*rand();
        A  = 0.5 + 1.0*rand();
        T1 = 50. + 50.*rand();
        T2 = 10. + 40.*rand();
        t1 = T1*rand();
        t2 = T2*rand();

        u_periodic[:,j] = A*sin((t-t1)*(2*pi/T1)).*sin((t-t2)*(2*pi/T2));
    end

    return u_periodic
end

#---------- 2. rate network ----------#
if target_type == "ratemodel"
    u_rm   = zeros(length(t),N);
    g_rm = 5;
    M = g_rm*genw_sparse(N,0,1,p)/sqrt(N*p);
    b_rm = 1/4;

    # the mean of incoming synaptic connections to a neuron is zero
    for i = 1:N
      idx = abs(M[i,:]) .> 0;
      M[i,idx] = M[i,idx] - sum(M[i,idx])/sum(idx);
    end

    # random initial condition
    u_rm[1,:] = 0.2*rand(N)';

    function phi(x)
        x[x .< 0.] = 0.;
        y = sqrt(x)/pi;
        return y
    end

    for i = 1:length(t)-1
        u_rm[i+1,:] = (1-dt*b_rm)*u_rm[i,:] + (dt*b_rm*M*phi(u_rm[i,:]'))';
    end
    return u_rm
end

#---------- 3. Ornstein-Ulenbeck ----------#
if target_type == "ou"
    u_ou    = zeros(length(t),N);
    b_ou  = 1/20;
    mu    = 0.0;
    sig   = 0.3;

    for j = 1:N
        for i = 1:length(t)-1
            u_ou[i+1,j] = u_ou[i,j]+b_ou*(mu-u_ou[i,j])*dt + sig*sqrt(dt)*randn();
        end
    end
    return u_ou
end

end


text:
const QUANDL_URL = "https://www.quandl.com/api/v3/datasets"  # for querying quandl's servers

"""
Set up Quandl user account authorization. Run once passing your Quandl API key, and it will be saved for future use.

    quandl_auth(key::String=""; authfile::String=expanduser("~/quandl-auth"))::String

"""
function quandl_auth(key::String=""; authfile::String=expanduser("~/quandl-auth"))::String
    if key == ""
        if isfile(authfile)
            key = read(authfile, String)
        end
    else
        f = open(authfile, "w")
        write(f, key)
        close(f)
    end
    return key
end

"""
Download time series data from Quandl as a TS object.

```
quandl(code::String;
       from::String="",
       thru::String="",
       freq::Char='d',
       calc::String="none",
       sort::Char='a',
       rows::Int=0,
       auth::String=quandl_auth())::TS
```

# Arguments

- `code::String;`: quandl code to query
- `from::String=""`: start date from which the data should begin, expressed as a string of format yyyy-mm-dd
- `thru::String=""`: end date through which the data should continue, expressed as a string of format yyyy-mm-dd
- `freq::Char='d'`: frequency at which the data should occur (one of 'd', 'w', 'm', 'q', or 'a')
- `calc::String="none"`: the calculation passed through to the Quandl engine (one of "none", "diff", "rdiff", "cumul", or "normalize")
- `sort::Char='a'`: direction in which to sort results (one of 'a' or 'd')
- `rows::Int=0`: number of rows to return from query (default corresponds to all rows)
- `auth::String=quandl_auth()`: authorization token for the Quandl API dedicated to a user account

"""
function quandl(code::String;
                from::String="",
                thru::String="",
                freq::Char='d',
                calc::String="none",
                sort::Char='a',
                rows::Int=0,
                auth::String=quandl_auth())::TS
    # Check arguments =========================================================
    @assert from=="" || (from[5]=='-' && from[8]=='-') "Argument `from` has invlalid format."
    @assert thru=="" || (thru[5]=='-' && thru[8]=='-') "Argument `thru` has invlalid format."
    @assert freq in ['d', 'w', 'm', 'q', 'a'] "Invalid `freq` argument (must be in ['d', 'w', 'm', 'q', 'a'])."
    @assert calc in ["none","diff","rdiff","cumul","normalize"] "Invalid `calc` argument."
    @assert sort  == 'a' || sort == 'd' "Argument `sort` must be either \'a\' or \'d\'."
    if rows != 0 && (from != "" || thru != "")
        error("Cannot specify `rows` and date range (`from` or `thru`).")
    end
    # Format URL ===============================================================
    sort_arg = (sort=='a' ? "asc" : "des")
    freq_arg = (freq=='d' ? "daily" : (freq=='w' ? "weekly" : (freq=='m' ? "monthly" : (freq=='q' ? "quarterly" : (freq=='a' ? "annual" : "")))))
    if rows == 0
        fromstr = from == "" ? "" : "&start_date=$from"
        thrustr = thru == "" ? "" : "&end_date=$thru"
        url = "$QUANDL_URL/$code.csv?$(fromstr)$(thrustr)&order=$sort_arg&collapse=$freq_arg&transform=$calc&api_key=$auth"
    else
        url = "$QUANDL_URL/$code.csv?&rows=$rows&order=$sort_arg&collapse=$freq_arg&transform=$calc&api_key=$auth"
    end
    indata = csvresp(HTTP.get(url), sort=sort)
    return TS(indata[1], indata[2], indata[3][2:end])
end

"""
Download Quandl metadata for a database and dataset into a Julia Dict object.

    quandl_meta(database::String, dataset::String)::Dict{String,Any}

"""
function quandl_meta(database::String, dataset::String)::Dict{String,Any}
    resp = HTTP.get("$QUANDL_URL/$database/$dataset/metadata.json")
    @assert resp.status == 200 "Error downloading metadata from Quandl."
    return JSON.parse(String(resp.body))["dataset"]
end

"""
Search Quandl for data in a given database, `db`, or matching a given query, `qry`.

    quandl_search(;db::String="", qry::String="", perpage::Int=1, pagenum::Int=1)

"""
function quandl_search(;db::String="", qry::String="", perpage::Int=1, pagenum::Int=1)
    @assert db!="" || qry!="" "Must enter a database or a search query."
    dbstr = db   == "" ? "" : "database_code=$db&"
    qrystr = qry  == "" ? "" : "query=$(replace(qry, ' ' => '+'))&"
    resp = HTTP.get("$QUANDL_URL.json?$(dbstr)$(qrystr)per_page=$perpage&page=$pagenum")
    @assert resp.status == 200 "Error retrieving search results from Quandl"
    return JSON.parse(String(resp.body))
end