text:
SEED_Char(io::IO, BUF::SeisIOBuf, nb::UInt16) = replace(String(fastread(io, nb)),
                                            ['\r', '\0'] =>"")

function SEED_Unenc!(io::IO, S::GphysData, c::Int64, xi::Int64, nb::UInt16, nx::UInt16)
  buf = getfield(BUF, :buf)
  checkbuf_8!(buf, xi)
  x = getindex(getfield(S, :x), c)
  fast_readbytes!(io, buf, nb)
  T::Type = (if BUF.fmt == 0x01
      Int16
    elseif BUF.fmt == 0x03
      Int32
    elseif BUF.fmt == 0x04
      Float32
    else
      Float64
    end)
  xr = reinterpret(T, buf)
  if BUF.swap
    x[xi+1:xi+nx] .= bswap.(xr[1:nx])
  else
    copyto!(x, xi+1, xr, 1, nx)
  end
  setfield!(BUF, :k, Int64(nx))
  return nothing
end

function SEED_Geoscope!(io::IO, BUF::SeisIOBuf)
  mm = 0x0fff
  gm = BUF.fmt == 0x0d ? 0x7000 : 0xf000
  for i = 0x0001:BUF.n
    x = BUF.swap ? bswap(fastread(io, UInt16)) : fastread(io, UInt16)
    m = Int32(x & mm)
    g = Int32((x & gm) >> 12)
    ex = -1*g
    setindex!(BUF.x, ldexp(Float64(m-2048), ex), i)
  end
  BUF.k = BUF.n
  return nothing
end

function SEED_CDSN!(io::IO, BUF::SeisIOBuf)
  for i = 0x0001:BUF.n
    x = BUF.swap ? bswap(fastread(io, UInt16)) : fastread(io, UInt16)
    m = Int32(x & 0x3fff)
    g = Int32((x & 0xc000) >> 14)
    mult = 4^g * g==3 ? 2 : 1
    m -= 0x1fff
    setindex!(BUF.x, m*mult, i)
  end
  BUF.k = BUF.n
  return nothing
end

function SEED_SRO!(io::IO, BUF::SeisIOBuf)
  for i = 0x0001:BUF.n
    x = BUF.swap ? bswap(fastread(io, UInt16)) : fastread(io, UInt16)
    m = Int32(x & 0x0fff)
    g = Int32((x & 0xf000) >> 12)
    if m > 0x07ff
      m -= 0x1000
    end
    ex = -1*g + 10
    setindex!(BUF.x, ldexp(Float64(m), ex), i)
  end
  BUF.k = BUF.n
  return nothing
end

function SEED_DWWSSN!(io::IO, BUF::SeisIOBuf)
  for i = 0x0001:BUF.n
    x = signed(UInt32(BUF.swap ? bswap(fastread(io, UInt16)) : fastread(io, UInt16)))
    BUF.x[i] = x > 32767 ? x - 65536 : x
  end
  BUF.k = BUF.n
  return nothing
end

# Steim1 or Steim2
function SEED_Steim!(io::IO, BUF::SeisIOBuf, nb::UInt16)
  x = getfield(BUF, :x)
  buf = getfield(BUF, :buf)
  ff = getfield(BUF, :uint32_buf)
  nc = Int64(div(nb, 0x0040))
  ni = div(nb, 0x0004)
  fast_readbytes!(io, buf, nb)
  (ni > lastindex(ff)) && resize!(ff, ni)
  fillx_u32_be!(ff, buf, ni, 0)
  k = zero(Int64)       # number of values read
  x0 = zero(Float32)    # first data value
  xn = zero(Float32)    # last data value
  a = zero(UInt8)       # byte offset to first data in each UInt32
  b = zero(UInt8)       # number of reads in each UInt32
  c = zero(UInt8)       # length (in bits) of each read in each UInt32
  d = zero(UInt8)       # amount of right bitshift
  fq = zero(Float32)    # Float32 data placeholder
  m = zero(UInt8)       # counter to reads in each UInt32
  dnib = zero(UInt8)    # dnib; two-bit secondary encoding flag
  q = zero(Int32)       # Int32 data placeholder
  u = zero(UInt32)      # placeholder for bit-shifted UInt32
  ck = zero(UInt8)      # two-bit primary encoding flag
  z = zero(UInt32)      # UInt32 containing nibbles
  χ = zero(UInt32)      # packed UInt32
  r = zero(Int64)       # "row" index to "matrix" of UInt32s
  for i = 1:nc
    z = getindex(ff, 1+r)
    for j = 1:16
      χ = getindex(ff, j+r)
      ck = UInt8((z >> steim[j]) & 0x03)

      # Steim1 and Steim2 are the same here
      if ck == 0x01
        a = 0x00
        b = 0x08
        c = 0x04

      # Steim1 for ck > 0x01
      elseif BUF.fmt == 0x0a
        a = 0x00
        if ck == 0x02
          b = 0x10
          c = 0x02
        elseif ck == 0x03
          b = 0x20
          c = 0x01
        end

      # Steim2
      else
        dnib = UInt8(χ >> 0x0000001e)
        if ck == 0x02
          a = 0x02
          if dnib == 0x01
            b = 0x1e
            c = 0x01
          elseif dnib == 0x02
            b = 0x0f
            c = 0x02
          elseif dnib == 0x03
            b = 0x0a
            c = 0x03
          end
        elseif ck == 0x03
          if dnib == 0x00
            a = 0x02
            b = 0x06
            c = 0x05
          elseif dnib == 0x01
            a = 0x02
            b = 0x05
            c = 0x06
          else
            a = 0x04
            b = 0x04
            c = 0x07
          end
        end
      end
      if ck != 0x00
        u = χ << a
        m = zero(UInt8)
        d = 0x20 - b
        while m < c
          k = k + 1
          q = signed(u)
          q >>= d
          fq = Float32(q)
          setindex!(x, fq, k)
          m = m + 0x01
          u <<= b
        end
      end
      if i == 1
        if j == 2
          x0 = Float32(signed(χ))
        elseif j == 3
          xn = Float32(signed(χ))
        end
      end
    end
    r = r+16
  end

  # Cumsum by hand
  setindex!(x, x0, 1)
  xa = copy(x0)
  @inbounds for i1 = 2:k
    xa = xa + getindex(x, i1)
    setindex!(x, xa, i1)
  end

  # Check data values
  if isapprox(getindex(x, k), xn) == false
    println(stdout, string("RDMSEED: data integrity -- Steim-",
                            getfield(BUF, :fmt) - 0x09, " sequence #",
                            String(copy(getfield(BUF, :seq))),
                            " integrity check failed, last_data=",
                            getindex(getfield(BUF, :x), k),
                            ", should be xn=", xn))
  end
  setfield!(BUF, :k, k)
  return nothing
end


text:
module TestPrecompile
    using Memoize
    run = 0
    @memoize function forgetful(x)
        global run += 1
        return true
    end

    forgetful(1)
end # module

ext: jl


text:
using Memoize

infile = "./src/Day7/input.txt"
input = readlines(infile)

function parseinput(input)
    rules = Dict{String, Vector{Pair{String, Int}}}()
    splits = split.(input, " bags contain ")
    left, right = string.(first.(splits)) , string.(last.(splits))
    reg = Regex(" bag[s\\.,\\s]+")
    for i in 1:length(left)
        pieces = split(right[i], reg)
        edges = Vector{Pair{String, Int}}()
        for r in pieces
            if length(r) == 0
                continue
            end
            if occursin("no other", r)
                continue
            end
            val = parse(Int, match(r"\d", r).match)
            col = string(match(Regex(join(left, "|")), r).match)
            edge = (col => val)
            push!(edges, edge)
        end
        if length(edges) > 0
            rules[left[i]] = edges
        end
    end
    return rules
end

rules = parseinput(input)

function check_rules!(visited, rules, target)
    # find all bags which can contain the target
    containers = findall(x -> any(==(target), first.(x)), rules)
    # if no bags can contain the target, we're done, return the set found so far
    if length(containers) == 0
        println("Nope, no bags can contain $target")
        return visited
    end
    println("The following bags can contain $target:")
    println.(containers)
    # include the set of containers we've found in our list
    push!(visited, containers...)
    # now find *their* containers, and push into visited
    for p in containers
        check_rules!(visited, rules, p)
    end
    return visited
end

can_reach_shinygold = check_rules!(Set{String}(), rules, "shiny gold")

function check_contents(rules, target)
    if !haskey(rules, target)
        return 0
    end
    innards = rules[target]
    total = 0
    for (col, val) in innards
        #println("Adding $val bags of colour $col")
        total += val * (1 + check_contents(rules, col))
        #println(total)
    end
    return total
end

bag_count = check_contents(rules, "shiny gold")


text:
mutable struct PDEStruct
    codename::String; # Exasim
    version::String;  # Set version
    appname::String;  # application name
    platform::String; # CPU or GPU
    cpucompiler::String; # Path to CPU compiler
    mpicompiler::String; # Path to MPI compiler
    gpucompiler::String; # Path to GPU compiler
    mpirun::String;      # Path to MPI run command and MPI run options
    cpuflags::String;    # options for CPU compiler
    gpuflags::String;    # options for GGU compiler
    model::String;# used to indicate PDE model
    modelfile::String;# PDE model file name

    preprocessmode::IntP; # preprocessing mode
    mpiprocs::IntP; # number of MPI ranks
    nd::IntP; # physical dimension
    nc::IntP; # number of compoments of (u, q, p)
    ncu::IntP;# number of compoments of (u)
    ncq::IntP;# number of compoments of (q)
    ncp::IntP;# number of compoments of (p)
    nco::IntP;# number of compoments of (odg)
    nch::IntP;# number of compoments of (uhat)
    ncx::IntP;# number of compoments of (xdg)
    ncw::IntP;# number of compoments of (wdg)
    nce::IntP;# number of compoments of (output field)
    neb::IntP;# number of element blocks for parallel computation
    nfb::IntP;# number of face blocks for parallel computation
    elemtype::IntP; # type of elements
    nodetype::IntP; # type of nodes
    hybrid::IntP; # discretization method
    tdep::IntP; # flag for steady-state or time-dependent problem
    wave::IntP; # flag for wave problem
    linearproblem::IntP; # flag for linear problem
    debugmode::IntP; # flag for debug mode
    stgNmode::IntP; # number of synthetic turbulence generation modes

    porder::IntP; # polymnomial degree
    pgauss::IntP; # Gauss quadrature polynomial degree
    temporalscheme::IntP; # temporal scheme
    torder::IntP; # temporal order of accuracy
    nstage::IntP; # number of RK stages
    convStabMethod::IntP; # flag for convection stabilization method
    diffStabMethod::IntP; # flag for diffusion stabilization method
    rotatingFrame::IntP;  # flag for rotating frame
    viscosityModel::IntP; # flag for viscosity model
    SGSmodel::IntP; # flag for SGS model
    ALE::IntP; # flag for ALE
    AV::IntP; # flag for artificial viscosity
    nonlinearsolver::IntP; # flag for nonlinear solver (Newton default)
    linearsolver::IntP; # flag for linear solver (GMRES default)
    NLiter::IntP; # maximum number of nonlinear iterations
    linearsolveriter::IntP; # maximum number of linear iterations
    GMRESrestart::IntP; # number of GMRES restarts
    GMRESortho::IntP; # GMRES orthogonalization method
    preconditioner::IntP; # flag for preconditioner method
    precMatrixType::IntP; # flag for type of preconditioner matrix
    NLMatrixType::IntP;
    runmode::IntP; # flag for run mode
    tdfunc::IntP; # flag for time-dependent function associated with time-derivative
    source::IntP; # flag for source-term function

    matvecorder::IntP; # flag for order of accuracy for matrix-vector multiplication
    RBdim::IntP; # reduced-basis dimension
    saveSolFreq::IntP; # flag for how often the solution saved in binary files
    saveSolOpt::IntP; # option for how the solution be saved: 0 -> u only, 1 -> u and q
    timestepOffset::IntP; # for restarting the simulation from the saved solution
    saveSolBouFreq::IntP; # how often the solution be saved on a particular boundary
    ibs::IntP; # the boundary on which the solution be saved
    compudgavg::IntP; # flag if time-average solution is computed
    extFhat::IntP;
    extUhat::IntP;

    time::FloatP; # starting time (usually 0, however >0 if restarting from the saved solution)
    NLparam::FloatP;
    NLtol::FloatP; # nonlinear solver tolerance
    linearsolvertol::FloatP; # linear solver tolerance
    matvectol::FloatP; # matrix-vector product tolerance

    flag::Array{IntP,2};   # flag parameters
    problem::Array{IntP,2};# problem parameters
    boundaryconditions::Array{IntP,2};# a list of boundary condition numbers
    stgib::Array{IntP,2};  # synthetic turbulence

    dt::Array{FloatP,1};      # time steps
    tau::Array{FloatP,1}; # stabilization parameters
    factor::Array{FloatP,2};  # factors
    physicsparam::Array{FloatP,2}; # physical parameters
    solversparam::Array{FloatP,2}; # solvers parameters
    stgdata::Array{FloatP,2}; # synthetic turbulence
    stgparam::Array{FloatP,2}; # synthetic turbulence
    externalparam::Array{FloatP,2};    # external parameters
    uinf::Array{FloatP,2};

    gmsh::String;
    metis::String;
    paraview::String;
    visfilename::String;
    visscalars;
    visvectors;
    viselem;
    soltime;
    visdt::FloatP;
    PDEStruct() = new();
end

function initializepde(version)
    pde = PDEStruct();

    pde.cpucompiler = "g++";
    pde.mpicompiler = "mpicxx";
    pde.gpucompiler = "nvcc";
    pde.mpirun = "mpirun";
    pde.gmsh = "gmsh";
    pde.metis = "mpmetis";
    pde.paraview = "/Applications/ParaView-5.8.1.app/Contents/MacOS/paraview";

    pde.codename = "Exasim";
    pde.version = version;
    pde.appname = "app";
    pde.platform = "cpu";
    #pde.cpuflags = "-O2 -Wall -ldl -lm -lblas -llapack";
    pde.cpuflags = "-O2 -ldl -lm -lblas -llapack";
    pde.gpuflags = "-lcudart -lcublas";
    pde.modelfile = "";
    pde.preprocessmode = 1;
    pde.mpiprocs = 1;
    pde.nd = 1;
    pde.nc = 1;
    pde.ncu = 1;
    pde.ncq = 0;
    pde.ncp = 0;
    pde.nco = 0;
    pde.nch = 1;
    pde.ncx = 1;
    pde.ncw = 0;
    pde.nce = 0;
    pde.neb = 512*8;
    pde.nfb = 512*32;
    pde.elemtype = 1;
    pde.nodetype = 1;
    pde.model="ModelD";
    pde.hybrid = 0;
    pde.tdep = 0;
    pde.wave = 0;
    pde.linearproblem = 0;
    pde.debugmode = 0;
    pde.stgNmode = 0;
    pde.porder = 1;
    pde.pgauss = 2;
    pde.temporalscheme = 0;
    pde.torder = 1;
    pde.nstage = 1;
    pde.convStabMethod = 0;
    pde.diffStabMethod = 0;
    pde.rotatingFrame = 0;
    pde.viscosityModel = 0;
    pde.SGSmodel = 0;
    pde.ALE = 0;
    pde.AV = 0;
    pde.nonlinearsolver = 0;
    pde.linearsolver = 0;
    pde.NLiter = 20;
    pde.linearsolveriter = 200;
    pde.GMRESrestart = 25;
    pde.GMRESortho = 0;
    pde.preconditioner = 0;
    pde.precMatrixType = 0;
    pde.NLMatrixType = 0;
    pde.runmode = 0;
    pde.tdfunc = 1;
    pde.source = 1;
    pde.matvecorder = 1;
    pde.RBdim = 5;
    pde.saveSolFreq = 1;
    pde.saveSolOpt = 1;
    pde.timestepOffset = 0;
    pde.saveSolBouFreq = 0;
    pde.ibs = 0;
    pde.compudgavg = 0;
    pde.extFhat = 0;
    pde.extUhat = 0;

    pde.time = 0.0;
    pde.NLparam = 0.0;
    pde.NLtol = 1e-6;
    pde.linearsolvertol = 1e-3;
    pde.matvectol = 1e-3;

    pde.flag = [0 0];
    pde.problem = [0 0];
    pde.boundaryconditions = [0 0];
    pde.stgib = [0 0];

    pde.tau = [1.0]; # stabilization parameters
    pde.dt = [0.0];  # time steps
    pde.factor = [0.0 0.0];  # factors
    pde.physicsparam = [0.0 0.0]; # physical parameters
    pde.solversparam = [0.0 0.0]; # solvers parameters
    pde.stgdata = [0.0 0.0]; # synthetic turbulence
    pde.stgparam = [0.0 0.0]; # synthetic turbulence
    pde.externalparam = [0.0 0.0]; #
    pde.uinf = [0.0 0.0]; #

    pde.visfilename = "dataout/output";
    pde.visscalars = [];
    pde.visvectors = [];
    pde.viselem = [];
    pde.soltime = [1];
    pde.visdt = 1.0;
    return pde;
end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
module SDPNAL

using SparseArrays
using MATLAB

export sdpnalplus

# See solver_main_default/sdpnalplus.m
const ALLOWED_OPTIONS = [
    "tol",
    "sigma",
    "maxiter",
    "maxtime",
    "printlevel",
    "rescale",
    "stopoption",
    "AATsolve",
    "BBTsolve",
    "beta",
    "SSNprecond",
    "ADMmaxiter",
    "ADMtol"
]

_array(x::AbstractMatrix) = x
_array(x::Vector) = x
_array(x::Float64) = [x]

# TODO log in objective, OPTION, initial iterates X0, y0, Z0
# Solve the primal/dual pair
# min c'x,      max b'y
# s.t. Ax = b,   c - A'x ∈ K
#       x ∈ K
function sdpnalplus(
    blk::Matrix,
    At::Vector{<:Union{Matrix{Float64}, SparseMatrixCSC{Float64}}},
    C::Vector{<:Union{Matrix{Float64}, Vector{Float64}}},
    b::Vector{Float64},
    L = [], U = [],
    Bt = [], l = [], u = []; kws...)

    #C::Vector{<:Union{Vector{Float64}, SparseVector{Float64}}}, b::Vector{Float64})

    options = Dict{String, Any}(string(key) => value for (key, value) in kws)
    @assert all(i -> size(At[i], 2) == length(b), 1:length(At))
    @assert length(At) == size(blk, 1)
    #@assert all(i -> size(A[i], 1) == dim(blk[i, 1], blk[i, 2]), 1:length(A))
    #@assert all(i -> length(C[i], 1) == dim(blk[i, 1], blk[i, 2]), 1:length(A))
    # There are 10 output arguments so we use `10` below
    obj, X, s, y, Z1, Z2, y2, v, info, runhist = mxcall(:sdpnalplus, 10, blk, At, C, b, L, U, Bt, l, u, options)
    return obj, _array.(X), _array(s), _array(y), _array.(Z1), _array.(Z2), y2, v, info, runhist
end

include("MOI_wrapper.jl")

end # module


text:
## module Plots

## Deprecated...

## use nvd3.org javascript to make graphics
## pluses: easy to make
## minuses: not really composable. For example, no way to add scatter and plot, plot and plot,...
## so no plot! type functions

## Deprecate... using Plots is more flexible

export tickFormat
export plot, scatter

## Type to hold D3 commands
type D3Plot
    q
    receiver
    D3Plot() = new([""], "chart")
end

# Co-opt the REPL display
Base.display(::Base.REPL.REPLDisplay, ::MIME"text/plain", p::D3Plot) = display(_display, p)
function Base.writemime(io::IO, ::MIME"text/html", x::D3Plot)
    tpl = Mustache.template_from_file(Pkg.dir("DThree", "tpl", "d3.html"))
    Mustache.render(io, tpl, Dict(:script=>get(x), :title=>"Graphic"))
end

function browse(x::D3Plot; style="")
    f = tempname() * ".html"
    io = open(f, "w")
    writemime(io, "text/html", x)
    close(io)
    open_browser_window(f)
end

import Base.push!, Base.get, Base.*

push!(x::D3Plot, cmd::AbstractString) = push!(x.q, cmd)
push!(x::D3Plot, cmd::D3) = push!(x.q, cmd.render())
push!(x::D3Plot) = cmd -> push!(x, cmd)  ## d3... |> push!(q)
## p * cmd; not p = p*cmd
*(p::D3Plot, x::D3) = push!(p, x.render())
*(p::D3Plot, x::AbstractString) = push!(p, x)

function nv_addGraph(x::D3Plot)
    x.q = ["", "nv.addGraph(function() {", x.q[2:end]..., "nv.utils.windowResize(chart.update)", "return chart;});" ]
end

function get(x::D3Plot)
    if length(x.q) > 1
        out = join(x.q[2:end], "\n")
    else
        out = ""
    end
    out
end
clear(x::D3Plot) =  x.q=[""]

_recycle(xs, n) = [xs[mod1(i, length(xs))] for i in 1:n]

## return tickFormat based on variable type
"""

Customize format of ticks on axis.

Example `plot(..., tickformat = tickFormat(1.0))`
"""
tickFormat(x::Real) = D3().format(",.02f")
tickFormat(x::Integer) = D3().format(",.r")
tickFormat(x::DateTime) = asis("function(d) { return d3.time.format('%x')(new Date(d))}")

tickFormat{T}(x::Vector{T}) = tickFormat(x[1])
tickFormat{T}(x::DataArray{T}) = tickFormat(x[1])
tickFormat(x) = tickFormat(collect(x))

## do we convert time to seconds for display? Only if a time
iftime(x::Any) = x
iftime(x::DateTime) = (x - DateTime(1970,1,1,0,0,0)).value
iftime{T}(x::DataArray{T}) = iftime(x[1])

## Some utilities
function _labels(q::D3Plot;
                main::AbstractString="",
                x::AbstractString="",
                y::AbstractString="")

    d3 = D3()
    d3.receiver(q.receiver).__("xAxis").axisLabel(x)     |> push!(q)
    d3.receiver(q.receiver).__("yAxis").axisLabel(y)     |> push!(q)
end

## Set tickFormat for axis based on element type.
## use `nothing` to not customize the format
function _tickFormat(q::D3Plot, x=nothing, y=nothing)
    d3 = D3()
    x != nothing && d3.receiver(q.receiver).__("xAxis").tickFormat(x)     |> push!(q)
    y != nothing && d3.receiver(q.receiver).__("yAxis").tickFormat(y)     |> push!(q)
end

"show legend?"
legend(q::D3Plot, flag::Bool=true) = D3()._("showLegend",flag) |> push!(q)

"use interactive guideline"
guideline(q::D3Plot, flag::Bool=true) = D3()._("useInteractiveGuideline",flag) |> push!(q)

function _showDist(q::D3Plot, x::Bool=true, y::Bool=true)
    D3()._("showDistX",x)._("showDistY",y) |> push!(q)
end

function _margins(q::D3Plot, d::Dict)
    e = Dict()
    for key in [:left, :top, :right, :bottom]
        if (haskey(d, key) && d[key] != nothing) e[key] = d[key] end
    end

    length(e) > 0 && D3()._("margin", JSON.json(e)) |> push!(q)

end

datum(q::D3Plot, d::AsIs) = D3().data(d) |> push!(q)

function dict_from_dfrow(df::DataFrame)
    d = Dict{AbstractString, Any}()
    for nm in names(df)
        d[string(nm)] = df[1,nm]
    end
    d
end

## values from a data frame
values_from_df(df::DataFrame) = [dict_from_dfrow(df[i,:]) for i in 1:size(df)[1]]

##################################################

## mave a variable `var` of type `models.lineChart`. Other values pushed onto `var`
function line_chart(q::D3Plot, var="chart", receiver="nv")
    D3().var(var).receiver("nv").__("models")._("lineChart")  |> push!(q)
end

## create data for lineType
## basically push everything down to data frame stye which has
## d,:f to split into groups and d for just a single line.

## dataframe. No grouping variable
## legend and colors are scalar
function create_datum(::Type{Val{:lineType}},
                      d::DataFrame,
                      f::Void,
                      vars=[:x, :y];
                      legend=":x",
                      color="blue",
                      area::Bool=false
                      kwargs...)

    data = [Dict(:values=>values_from_df(d),
	         :key=>legend,
	         :color=>color,
                 :area=>area)
            ]

    data |> JSON.json |> asis
end

## dataframe with grouping by f
function create_datum(::Type{Val{:lineType}},
                      d::DataFrame,
                      f::Symbol,
                      vars=[:x, :y];
                      legend=Dict(),  # d[level] = "text"
                      colors=Dict(),  # d[level] = :color,
                      area = Dict()
                      kwargs...)
    ## who should group_by, but instead do this manually

    fs = d[f]
    levels = sort(unique(fs))

    if isa(legend, Vector)
        legs = copy(legend)
        legend = Dict()
        [legend[level] = legs[mod1(i, length(legs))] for (i, level) in enumerate(levels)]
    else
        for level in levels
            if !haskey(legend, level)  legend[level] = string(level)  end
        end
    end

    if isa(colors, Vector)
        cols = copy(colors)
        colors = Dict()
        [colors[level] = cols[mod1(i, length(cols))] for (i, level) in enumerate(levels)]
    else
        for (i, level) in enumerate(levels)
            if !haskey(colors, level) colors[level] = fifty_shades(i)  end
        end
    end

    if isa(area, Vector)
        ars = copy(area)
        area = Dict()
        [area[level] = ars[mod1(i, length(cols))] for (i, level) in enumerate(levels)]
    else
        for (i, level) in enumerate(levels)
            if !haskey(area, level) area[level] = false  end
        end
    end

    data = [Dict(:values=>values_from_df(d[fs .== level, vars]),
	         :key=>legend[level],
	         :color=>colors[level],
                 :area => area[level]
                 )
            for level in levels
            ]

    data |> JSON.json |> asis
end

function create_datum{T<:Real, S<:Real}(::Type{Val{:lineType}},
                                        xs::Vector{T},
                                        ys::Vector{S};
                                        kwargs...
                                        )

    d = DataFrame(x=xs, y=ys)
    create_datum(Val{:lineType},
                 d, nothing, [:x, :y];
                 kwargs...)

end

## lineTtype for multiple values specified with a factor
## legend and cols are dictionaries with keys as factors
function create_datum{T<:Real, S<:Real}(::Type{Val{:lineType}},
                                        xs::Vector{T}, ys::Vector{S},
                                        fs::Vector;
                                        kwargs...
                                        )

    d = DataFrame(x=xs, y=ys, f=fs)

    create_datum(Val{:lineType},
                 d, :f, [:x, :y];
                 kwargs...)

end

## create data from functions
## lineType should this be a parameter to dispatch on?
function create_datum(::Type{Val{:lineType}},
                      fs::Vector{Function},
                      from::Real, to::Real;
                      kwargs...
                      )

    m = length(fs)
    n = 151

    xs = linspace(from, to, n)
    d = DataFrame(x=repmat(xs, m))
    d[:y] = vcat([map(f, xs) for f in fs]...)
    d[:f] = repmat(map(string, fs), 1, n)'[:]

    create_datum(Val{:lineType},
                 d, :f, [:x, :y];
                 kwargs...)

end

## parametric t -> (f, g)
function create_datum(::Type{Val{:lineType}},
                      f::Function, g::Function,
                      from::Real, to::Real;
                      kwargs...
                      )

    ts = linspace(from, to, 250)
    xs = map(f, ts)
    ys = map(g, ts)

    create_datum(Val{:lineType}, xs, ys;
                 kwargs...)
end

## Main function
## `args...`, `kwargs...` passed to `create_datum`
## other arguments for customizing graphic
function _line_chart(args...;
                     selector="#chart svg",
                     labels=["x", "y"],
                     tickformat = [tickFormat(1.0), tickFormat(1.0)],
                     margins=Dict(:left=>nothing, :top=>nothing, :right=>nothing, :bottom=>nothing),
                     raise::Bool=true,
                     kwargs...)

    ## produce data from fs, return value in asis
    ## legend and color values are embedded in the data
    data = create_datum(Val{:lineType}, args...; kwargs...)

    d3 = D3()
    q = D3Plot()

    ## http://nvd3.org/ghpages/line.html
    line_chart(q)
    _labels(q, x=labels[1], y=labels[2])
    _tickFormat(q, tickformat...)
    _margins(q, margins)

    d3.select(selector).datum(data).transition().duration(500).call(asis("chart")) |> push!(q)

    nv_addGraph(q)
    raise && scf()
    q
end

## plot interface

# plot(xs, ys)
# plot(xs, ys, fs)
# plot(fn, a, b)
# plot(fns, a, b)
# plot(d, f, labels=[:x, :y])

"""

Plot a dot-to-dot graph. Can have several interfaces:

`plot(xs, ys)`:    plot `(x,y)`, connect dots
`plot(xs, ys, fs)` plot `(x,y)` for each group defined by `fs`, connect dots
`plot(f, from, to)` plot `y =f(x)` over domain `[from, to]`.
`plot([fs], from to)` plot `y=f_i(x)` over domain `[from, to]` for each `f` in `fs`
`plot(d, nothing, vars=[:x, :y])` plot `xs` `ys` defined in data frame `d`
`plot(d, :f, vars=[:x, :y])` group data by `:f`, then for each group make line chart.

The multiple line chart functions take `legend` and `colors` as a dictionary or vector of suitable length.

The single line hcart functions, take `legend` and `colors` as scalars.

"""
function plot(xs::Union{Range, Vector}, ys::Union{Range, Vector};
              labels=["x", "y"],
              selector="#chart svg",
              kwargs...)
    _line_chart(collect(xs), collect(ys);
                labels=labels, selector=selector,
                kwargs...)
end

function plot(xs::Union{Range, Vector}, ys::Union{Range, Vector}, fs::Vector;
              labels=["x", "y"],
              selector="#chart svg",
              kwargs...
              )

    _line_chart(collect(xs), collect(ys), fs;
                labels=labels, selector=selector,
                kwargs...)
end

plot(f::Function, from::Real, to::Real; kwargs...) = plot([f], from, to; kwargs...)

function plot(fs::Vector, from::Real, to::Real;
              labels=["x", "y"],
              selector="#chart svg",
              kwargs...
              )

    _line_chart(fs, from, to;
                labels=labels,
                selector=selector,
                kwargs...)
end

## parametric plot
function plot(f::Function, g::Function, from::Real, to::Real;
              labels=["x", "y"],
              selector="#chart svg",
              kwargs...)

      _line_chart(d, f, vars;
                labels=labels,
                selector=selector,
                kwargs...)
end

function plot(d::DataFrame, f::Union{Void, Symbol}, vars=[:x, :y];
              labels=["x", "y"],
              selector="#chart svg",
              kwargs...)

    _line_chart(d, f, vars;
                labels=labels,
                selector=selector,
                kwargs...)
end

function plot(f::Function, g::Function, a::Real, b::Real;
              labels=["x", "y"],
              selector="#chart svg",
              kwargs...)
    _line_chart(f, g, a, b;
                labels=labels,
                selector=selector,
                kwargs...)
end

##################################################

function scatter_chart(q::D3Plot, var="chart", receiver="nv")
    D3().var(var).receiver("nv").__("models")._("scatterChart")  |> push!(q)
end

## scatter chart
## shapes = ['circle', 'cross', 'triangle-up', 'triangle-down', 'diamond', 'square'],

## dataframe. No grouping variable
## legend and colors are scalar
function create_datum(::Type{Val{:scatterType}},
                      d::DataFrame,
                      f::Void,
                      vars=[:x, :y];
                      legend=":x",
                      colors=[:blue],
                      shapes=[:circle]
                      kwargs...)

    data = [Dict(:values=>values_from_df(d),
	         :key=>legend,
	         :color=>color,
                 :area=>area)
            ]

    data |> JSON.json |> asis
end

## dataframe with grouping by f
function create_datum(::Type{Val{:lineType}},
                      d::DataFrame,
                      f::Symbol,
                      vars=[:x, :y];
                      legend=Dict(),  # d[level] = "text"
                      colors=Dict(),  # d[level] = :color,
                      area = Dict()
                      kwargs...)
    ## who should group_by, but instead do this manually

    fs = d[f]
    levels = sort(unique(fs))

    if isa(legend, Vector)
        legs = copy(legend)
        legend = Dict()
        [legend[level] = legs[mod1(i, length(legs))] for (i, level) in enumerate(levels)]
    else
        for level in levels
            if !haskey(legend, level)  legend[level] = string(level)  end
        end
    end

    if isa(colors, Vector)
        cols = copy(colors)
        colors = Dict()
        [colors[level] = cols[mod1(i, length(cols))] for (i, level) in enumerate(levels)]
    else
        for (i, level) in enumerate(levels)
            if !haskey(colors, level) colors[level] = fifty_shades(i)  end
        end
    end

    if isa(area, Vector)
        ars = copy(area)
        area = Dict()
        [area[level] = ars[mod1(i, length(cols))] for (i, level) in enumerate(levels)]
    else
        for (i, level) in enumerate(levels)
            if !haskey(area, level) area[level] = false  end
        end
    end

    data = [Dict(:values=>values_from_df(d[fs .== level, vars]),
	         :key=>legend[level],
	         :color=>colors[level],
                 :area => area[level]
                 )
            for level in levels
            ]

    data |> JSON.json |> asis
end

function create_datum{T<:Real,S<:Real}(::Type{Val{:scatterType}},
                                       xs::Vector{T},
                                       ys::Vector{S};
                                       legend="y",
                                       colors=["blue"],
                                       sizes=[1],
                                       shapes=["circle"]
                                       )

    ## recycle colors, sizes, shapes if not arrays
    colors = recycle(colors, length(xs))
    sizes  = recycle(sizes, length(xs))
    shapes = recycle(shapes, length(xs))

    d = DataFrame(x=xs, y=ys, color=colors, size=sizes, shape=shapes)
    data = [Dict(:values => values_from_df(d),
                 :key => legend)]

    data |> JSON.json |> asis
end

function create_datum(::Type{Val{:scatterType}},
                      d::DataFrame,
                      f::Void;
                      vars=[1,2], #[:x, :y], ...
                      color=["blue"],
                      sizes=[1],
                      shapes=["circle"],
                      kwargs...
                      )
    m = size(d)[1]

    colors = _recycle(color, m)
    sizes = _recycle(sizes, m)
    shapes = _recycle(shapes, m)

    df = DataFrame()
    df[:x] = d[vars[1]]
    df[:y] = d[vars[2]]
    df[:color] = colors
    df[:size] = sizes
    df[:shape] = shapes

    data = [Dict(
                :color  => color,
                :values => [values_from_df(df)]
                )]
    data |> JSON.json |> asis
end

function create_datum(::Type{Val{:scatterType}},
                      d::DataFrame,
                      f::Symbol;
                      vars=[1,2], #[:x, :y], ...
                      legend=AbstractString[],
                      colors=[],
                      sizes=[],
                      shapes=[]
                      )

    df = DataFrame()
    df[:x] = d[vars[1]]
    df[:y] = d[vars[2]]

    levels = sort(unique(d[f]))
    nlevels = length(levels)

    # group attributes
    colors = length(colors) == 0 ? map(fifty_shades, 1:nlevels)  : _recycle(colors, nlevels)
    sizes  = length(sizes)  == 0 ? [1 for _ in 1:nlevels]        : _recycle(sizes, nlevels)
    shapes = length(shapes) == 0 ? ["circle" for _ in 1:nlevels] : _recycle(shapes, nlevels)
    legend = length(legend) == 0 ? [string(_) for _ in levels]   : _recycle(legend, nlevels)

    data = [Dict(:values => [Dict(:x=>d[d[:f] .== level, vars[1]][j],
                                 :y=>d[d[:f] .== level, vars[2]][j],
                                 :color => colors[i],
                                 :shape => shapes[i],
                                 :size => sizes[i]) for j in 1:sum(d[:f] .== level)],
                 :key    => legend[i])
            for (i, level) in enumerate(levels)]

    data |> JSON.json |> asis

end

## DRY THIS UP: scatter and line differe just by call to a {val}
function _scatter_chart(args...;
               selector="#chart svg",
               labels=["x", "y"],
               tickformat = [tickFormat(1.0), tickFormat(1.0)],
               margins=Dict(:left=>nothing, :top=>nothing, :right=>nothing, :bottom=>nothing),
               kwargs...)

    ## produce data from fs, return value in asis
    ## legend and color values are embedded in the data
    data = create_datum(Val{:scatterType}, args...; kwargs...)

    d3 = D3()
    q = D3Plot()

    ## http://nvd3.org/ghpages/line.html
    scatter_chart(q)
    _labels(q, x=labels[1], y=labels[2])
    _tickFormat(q, tickformat...)
    _margins(q, margins)
    _showDist(q, true, true)

    d3.select(selector).datum(data).transition().duration(500).call(asis("chart")) |> push!(q)

    nv_addGraph(q)
    q
end

# scatter(xs, ys)
# scatter(xs, ys, fs)
# scatter(d, f, labels=[:x,:y])

function scatter(xs::Union{Range, Vector}, ys::Union{Range, Vector};
                 labels=["x", "y"],
                 selector="#chart svg",
                 legend="y",
                 colors="blue",
                 kwargs...
                 )

    _scatter_chart(collect(xs), collect(ys);
                   labels=labels, selector=selector,
                   legend=legend, colors=colors,
                   kwargs...)
end

## scatterChart.
## data frame with columns x, y and  grouping variable. Also can have columns shape, size
scatterChart(x::DataFrame; args...)  = scatterChart(x, nothing, args...)
function scatterChart(x::DataFrame, f::Union{Void, Symbol};
                      selector = "#chart svg", # in web template
                      labels=nothing,    #  axis
                      legend=nothing,    # when grouping
                      colors=nothing
                      )

    ## get data from x
    function df_to_data(x::Union{DataFrame, SubDataFrame};
                        key=nothing, color=nothing)
        function do_row(i)
            a = Dict(:x => iftime(x[i,1]), :y => x[i, 2]) # now add size and shape?
##            a["size"] = haskey(x[1,:], "size") ? x[1,"size"] : 0.5
##            if haskey(x[1,:], "shape")
##                shape = x[1,"shape"]
##                a[:shape] = isa(shape, Integer) ? shapes[1 + shape % length(shapes)] : shape
##            end
            a
        end

        key = (key == nothing ? "group" : key)
        color = (color == nothing ? "#ff00ff" : color)

       Dict(:key=>key,
        :color=> color,
        :values=>[do_row(i) for i in 1:nrow(x)]
       )
    end

    if isa(f, Void)
         data = [df_to_data(x)] |> JSON.json |> asis
    else
        gp = groupby(x, f)
        ## need colors, legend key, ...
        if isa(legend, Void)
            legend = [gp[i][1,f] for i in 1:length(gp)]
        end
        if isa(colors, Void)
            colors = map(fifty_shades, 1:length(gp))
        end

        data = [df_to_data(gp[i], key=legend[i], color=colors[i]) for i in 1:length(gp)] |> JSON.json |> asis
    end

    labels = (isa(labels, Void) ? names(x)[1:2] : labels)

    d3 = D3()
    q = D3Plot()

    ## http://nvd3.org/ghpages/scatter.html
    d3.var("chart").receiver("nv")._("models.scatterChart")._("showDistX", true)._("showDistY",true).color(D3().scale_category10().range()) |> push!(q)
    d3.receiver("chart.xAxis").axisLabel(labels[1]).tickFormat(tickFormat(x[1])) |> push!(q)
    d3.receiver("chart.yAxis").axisLabel(labels[2]).tickFormat(tickFormat(x[2])) |> push!(q)
    d3.select(selector).datum(data).transition().duration(500).call(asis("chart")) |> push!(q)

    nv_addGraph(q)
    q
end

function barChart{S <: Real, T <: AbstractString}(xs::Vector{S}, labels::Vector{T};
                                          selector = "#chart svg", # in web template
                                          colors=nothing, #  which are needed?
                                          legend=nothing)

    if isa(legend, Void) legend = "Nothing" end

    data = [Dict(:key =>legend,
             :values => [Dict(:value=>iftime(x), :label => label) for (x,label) in zip(xs, labels)]
                 )] |> JSON.json |> asis

    d3 = D3()
    q = D3Plot()

    ## http://nvd3.org/ghpages/discreteBar.html
    d3.var("chart").receiver("nv")._("models.discreteBarChart").
    _("x", asis("function(d) {return d.label}")).
    _("y", asis("function(d) {return d.value}"))[:staggerLabels](true)[:tooltips](false)[:showValues](true)  |> push!(q)

    d3.select(selector).datum(data).transition().duration(500).call(asis("chart"))  |> push!(q)

    nv_addGraph(q)
    q
end

## x =DataFrame(xvar=.., yvar=..., ..., f=factor)
function stackedAreaChart(x::DataFrame, f::Symbol;
                          selector="#chart svg",
                          labels = nothing
                          )

    gp = groupby(x, f)
    legend = [gp[i][1,f] for i in 1:length(gp)]

    data = [Dict(:key => legend[i],
             :values => [[iftime(x),y] for (x,y) in zip(gp[i][1], gp[i][2])]
            ) for i in 1:length(gp)] |> JSON.json |> asis

    d3 = D3()
    q = D3Plot()

    ## http://nvd3.org/ghpages/stackedArea.html
    d3.var("chart").receiver("nv")._("models.stackedAreaChart").
    _("x", asis("function(d) {return d[0]}")).
    _("y", asis("function(d) {return d[1]}"))[:clipEdge](true)  |> push!(q)

    d3.receiver("chart.xAxis")[:showMaxMin](false).tickFormat(tickFormat(gp[1][1][1])) |> push!(q)
    d3.receiver("chart.yAxis").tickFormat(tickFormat(gp[1][2][1])) |> push!(q)
    d3.select(selector).datum(data).transition().duration(500).call(asis("chart")) |> push!(q)

    nv_addGraph(q)
    q

end

## end (modlue)


text:
@testset "EarthEnv LandCover" begin
    using RasterDataSources: rastername, rasterurl, rasterpath

    @test rastername(EarthEnv{LandCover}, 2; discover=true) == "landcover_complete_2.tif"
    landcover_path = joinpath(ENV["RASTERDATASOURCES_PATH"], "EarthEnv", "consensus_landcover")
    @test rasterpath(EarthEnv{LandCover}) == landcover_path
    @test rasterpath(EarthEnv{LandCover}, 2; discover=true) == joinpath(landcover_path, "landcover_complete_2.tif")
    @test rasterurl(EarthEnv{LandCover}, 2; discover=true) ==
        URI(scheme="https", host="data.earthenv.org", path="/consensus_landcover/with_DISCover/consensus_full_class_2.tif")
    getraster(EarthEnv{LandCover}, 2)
    @test isfile(joinpath(landcover_path, "landcover_partial_2.tif"))
end


text:
using WindMixing
using Flux
using JLD2
using FileIO
using Plots

FILE_NAME = "parameter_optimisation_18sim_windcooling_windheating_5params_BFGS"
FILE_PATH = joinpath(pwd(), "training_output", "$(FILE_NAME).jld2")
@assert isfile(FILE_PATH)
OUTPUT_PATH = joinpath(pwd(), "extracted_training_output", "$(FILE_NAME)_extracted.jld2")
OUTPUT_PATH = joinpath("D:\\University Matters\\MIT\\CLiMA Project\\OceanParameterizations.jl\\training_output", "$(FILE_NAME)_extracted.jld2")

extract_parameters_modified_pacanowski_philander_optimisation(FILE_PATH, OUTPUT_PATH)

file = jldopen(OUTPUT_PATH)
@info file["parameters"]
losses = file["losses"]

ν₀_initial = 1f-4
ν₋_initial = 1f-1
ΔRi_initial = 1f-1
Riᶜ_initial = 0.25f0
Pr_initial = 1f0

scalings = 1 ./ [ν₀_initial, ν₋_initial, ΔRi_initial, Riᶜ_initial, Pr_initial]

@info file["parameters"] ./ scalings

close(file)

plot(1:length(losses), losses, yscale=:log10)


text:
"""
Exception type indicating that an API function
returned a non-success code.
"""
struct VulkanError <: Exception
    msg::String
    code
end

Base.showerror(io::IO, e::VulkanError) = print(io, e.code, ": ", e.msg)

"""
    @check vkCreateInstance(args...)

Assign the expression to a variable named `_return_code`. Then, if the value is not a success code, return a [`VulkanError`](@ref) holding the return code.

"""
macro check(expr)
    msg = string("failed to execute ", expr)
    esc(:(@check $expr $msg))
end

macro check(expr, msg)
    esc(quote
        _return_code = $expr
        _code_int = Int32(_return_code)
        if _code_int < 0
            return VulkanError($msg, Result(_code_int))
        end
        _return_code
    end)
end

"""
Evaluate the given expression, and reexecute it while the local variable `_return_code` equals to `VK_INCOMPLETE`.
"""
macro repeat_while_incomplete(expr)
    esc(quote
        $expr
        while _return_code == VK_INCOMPLETE
            $expr
        end
    end)
end


text:
module Tutorial
using Revise #useful when revising code frequently

#prerequisites:
# in the directory where you want to store projects, from the REPL type ] generate "Tutorial"
# Install the following packages by typing the following into the REPL
# ] add Distributions LinearAlgebra Random Revise CUDA GLM BenchmarkTools DataFrames

#load well established
using Distributions, LinearAlgebra, DataFrames, Random, CUDA
sleep(0.1)

#NOTE: Always wrap Julia code in functions

function basics()
  #hello world
  println("Hello world")

  #The info macro is specifically designed to echo info back to the user
  @info "Hello again"

  #Simple math and assignment
  x = 3 + 4

  #string interpolation
  @info "x=$x"

  #compound assignment
  x *= 2
  @info "now x=$x" #14

  #Many ways to create vectors
  v1 = ones(5)
  v2 = [1.0 for i ∈ 1:5]
  v3 = (_->1.0).(1:5) #the dot here means we are executing the function for each value in 1:5
                      #_->1.0 is an anonymous function that takes a value as input and returns 1.
  v4 = Vector{Float64}(undef, 5); v4 .= 1.0 #The dot here is used for element wise operations

  #Use assertion to verify something is true- it will thrown an error if it is not
  @assert v1 == v2 == v3 == v4
  @info "Since we got here, v1 == v2 == v3 == v4 = $(v1 == v2 == v3 == v4)"

  #conditionals are easy
  if iseven(x)
    @info "x is $x so it is even"
  elseif isodd(x)
    @info "x is $x so it is odd"
  else
    @assert false #this pattern is useful for debugging
  end

  #The compiler works best when each variable maintains its type
  #we can enforce this by explicitly declaring a variable's type
  local M::Matrix{Float64} #declare M is a matrix of Float64 values
  local v::Vector{Float64} #declare v is a vector of Float64 values

  #assign random values to M
  M = rand(10_000,1_000) #M is a 10^6x10 matrix of random values, drawn from the std uniform
  v = rand(Normal(), 10_000) #v is a 10^6 matrix of random values,drawn from the std normal

  #Matrix algebra is easy
  Mtvsum = M'*v |> sum #|> is a pipe operator that sends output to a single argument function
  @info "M'*v |> sum: $(Mtvsum)"

  #we can also do element wise operations.
  #The .* operator expands the dimensions of v, if possible, to allow for the element wise operatoin
  Movsum = M .* v |> sum
  @info "M .* v |> sum $(Movsum)"
  @assert Mtvsum ≈ Movsum #we can use the ≈ operator to check if the answers are within rounding error

  sqrmat(X) = X'*X #functions are easy to write and will adapt to the input type
  MtMsum = M |> sqrmat |> sum
  vtvsum = v |> sqrmat |> sum #this works even though v is a vector, not a matrix
  x2 = x |> sqrmat #it even works for scalars
  @info "MtMsum: $MtMsum, vtvsum: $vtvsum, x2=$x2"

  #Suppose we wanted to define it for strings
  try # use  atry block to catch exceptions
    sqrmat("7")
  catch err
    @warn "sqrtmat(\"7\") failed with error $err"
  end

  #we can define a method which handles trhis by specifying the type in the argument
  #note that we need to explicitly set the parsing type, so we set an optional type argument
  #this is a very high performance pattern for coding
  sqrmat(X::String, ::Type{T}=Float64) where T = parse(T, X) |> sqrmat
  floatstring2 = sqrmat("7.0")
  intstring2 = sqrmat("7", Int)
  @info "sqrmat(\"7.0\"): $floatstring2, sqrmat(\"7\", Int): $intstring2"

  #suppose we wanted to square 5 matrices
  #time the results with @time (note @btime from BenchmarkTools library does a better job)
  print("Time squaring 20 matrices")
  M20 = [rand(Float32, 10_000,1000) for i ∈ 1:20]
  M20² = (Mᵢ->similar(Mᵢ,(1000,1000))).(M20) #prellocation can improve performance
  @time begin #use begin/end to create scoping blocks
    for (m20², m20ᵢ) ∈ zip(M20², M20) #loops are fast!
      m20² .= m20ᵢ |> sqrmat
    end
  end

  #for faster performnace, we can tell Julia to multi-thread
  print("Time squaring 20 matrices (parallel)")
  M20²mt = M20² |> deepcopy#use deepcopy to recursively copy contents of objects
  @time begin #use begin/end to create scoping blocks
    Threads.@threads for i ∈ 1:20 #only simple indexing allowed in loops like this
      M20²mt[i] .= M20[i] |> sqrmat
    end
  end

  print("Time squaring 20 matrices (gpu)")
  dM20² = M20² |> deepcopy .|> CuMatrix{Float32}
  dM20 = CUDA.@sync M20 .|> CuMatrix{Float32}
  @time CUDA.@sync begin #use begin/end to create scoping blocks
    for i ∈ 1:20 #only simple indexing allowed in loops like this
      dM20²[i] .= dM20[i] |> sqrmat
    end
  end

  #The DataFrames package is essential for working with data
  df = DataFrame(x = [rand() > 0.5 ? rand() : missing for i ∈ 1:1000],
                 g = rand([1,2,3,4], 1000))
  #it basically functions like a combination of other data packages
  @info "mean(df.x): $(mean(df.x)), mean(df.x |> skipmissing): $(mean(df.x |> skipmissing))"

  #sum by group
  println(combine(groupby(df, :g), [:x,:g] .=> (c)->sum(skipmissing(c))))

  #https://dataframes.juliadata.org/stable/
  #see more advanced DataFrames features at
  #
  # https://github.com/bkamins/Julia-DataFrames-Tutorial (requires Jupiter notebooks)

  #NOT covered in this tutorial but worth learning
  # structs, constructors, parametric types - see Julia docs
  # Serialization, Dates (see Julia docs- standard library)
  #Many packages- see their respective githubs:
  # CUDA (gpu)
  # Zygote (analytical derivitvies)
  # GLM (least squares)
  # BenchmarkTools (benchmarking)
  # Optim (optimizaiton)
  # CSV (text file i/o)
  # Weave (markdown)
  #  Plots, Gadfly, VegaLite (these are various plotting packages, don't need to know all)

end

basics()
end #module end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Distributed, Distributions, Serialization
include("/srv/git/rys_nucleosomes/aws/aws_wrangler.jl")

sib_e_pth = "/bench/PhD/NGS_binaries/BMI/sib_e"
rys_e_pth = "/bench/PhD/NGS_binaries/BMI/rys_e"
combined_e_pth = "/bench/PhD/NGS_binaries/BMI/combined_e"

@info "Assembling worker pool..."

#DISTRIBUTED CLUSTERs CONSTANTS
remote_machine = "10.0.0.3"
no_local_processes = 2
no_remote_processes = 6

@info "Spawning local cluster workers..."
worker_pool=addprocs(no_local_processes, topology=:master_worker)
remote_pool=addprocs([(remote_machine,no_remote_processes)], tunnel=true, topology=:master_worker)

worker_pool=vcat(worker_pool, remote_pool)

#AWS PARAMS
@info "Setting up AWS wrangling..."

security_group_name="calc1"
security_group_desc="calculation group"
ami="ami-0bb0f5609b995d39e"
skeys="AWS"
instance_type="c5a.24xlarge"
zone,spot_price=get_cheapest_zone(instance_type)
no_instances=5
instance_workers=48
bid=spot_price+.1

@assert bid >= spot_price

@info "Wrangling AWS instances..."
aws_ips = spot_wrangle(no_instances, bid, security_group_name, security_group_desc, skeys, zone, ami, instance_type)
@info "Giving instances 150s to boot..."
sleep(150)

# aws_ips=["18.223.214.57","18.224.4.81", "3.17.13.150","3.21.113.229","18.222.211.44"]

@info "Spawning AWS cluster workers..."
for ip in aws_ips
    instance_pool=addprocs([(ip, instance_workers)], tunnel=true, topology=:master_worker, sshflags="-o StrictHostKeyChecking=no")
    global worker_pool=vcat(worker_pool, instance_pool)
end

@info "Loading worker libraries everywhere..."
@everywhere using BioMotifInference,Random
@everywhere Random.seed!(myid()*100000)

#JOB CONSTANTS
funcvec=full_perm_funcvec
models_to_permute=10000
func_limit=25
push!(funcvec, BioMotifInference.perm_src_fit_mix)
push!(funcvec, BioMotifInference.random_decorrelate)

min_clamps=fill(.01,length(funcvec))
min_clamps[2:3].=.1 #perm_src_fit_mix & permute_mix
min_clamps[8]=.1 #difference_merge
max_clamps=fill(.5,length(funcvec))
max_clamps[6:7].=.15 #ss and am
max_clamps[9]=.15 #sim merge

initial_weights= ones(length(funcvec))./length(funcvec)
# override_weights=fill(.034375,length(funcvec))
# override_weights[6:9].=.1;override_weights[13:14].=.1625
# override_time=20.

args=[Vector{Tuple{Symbol,Any}}() for i in 1:length(funcvec)]
args[end-1]=[(:weight_shift_freq,0.),(:length_change_freq,1.),(:length_perm_range,1:1)]
args[end]=[(:iterates,50),(:source_permute_freq,.3),(:mix_move_range,1:10)]

instruct = Permute_Instruct(funcvec, initial_weights, models_to_permute, func_limit;min_clmps=min_clamps, max_clmps=max_clamps, args=args)

display_rotation=[true,10,1,[[:tuning_disp,:lh_disp,:src_disp],[:conv_plot,:liwi_disp,:ens_disp]]]

@info "Beginning nested sampling for sib ensemble..."
sib_e=deserialize(sib_e_pth*"/ens")
sib_logZ = converge_ensemble!(sib_e, instruct, worker_pool, converge_criterion="compression", converge_factor=150., backup=(true,25), tuning_disp=true,lh_disp=true,src_disp=true, disp_rotate_inst=display_rotation)
sib_logZ_err=sqrt(sib_e.Hi[end]/length(sib_e.log_Li))
sib_e=[]; Base.GC.gc();

@info "Beginning nested sampling for rys ensemble..."
rys_e=deserialize(rys_e_pth*"/ens")
rys_logZ = converge_ensemble!(rys_e, instruct, worker_pool, converge_criterion="compression", converge_factor=150., backup=(true,25), tuning_disp=true,lh_disp=true,src_disp=true, disp_rotate_inst=display_rotation)
rys_logZ_err=sqrt(rys_e.Hi[end]/length(rys_e.log_Li))
rys_e=[]; Base.GC.gc();

@info "Beginning nested sampling for combined ensemble..."
combined_e=deserialize(combined_e_pth*"/ens")
combined_logZ = converge_ensemble!(combined_e, instruct, worker_pool, converge_criterion="compression", converge_factor=150., backup=(true,25), tuning_disp=true,lh_disp=true,src_disp=true,disp_rotate_inst=display_rotation)
combined_logZ_err=sqrt(combined_e.Hi[end]/length(combined_e.log_Li))
combined_e=[]; Base.GC.gc();

@info "log evidence ratio of joint sib+rys / combined:"
ratio=(sib_logZ+rys_logZ)-combined_logZ
ratio_err=sqrt(sib_logZ_err^2 + rys_logZ_err^2 + combined_logZ_err^2)
println("$ratio ± $ratio_err, $(ratio/ratio_err) standard deviation significance")

serialize("/bench/PhD/NGS_binaries/BMI/dif_pos_learner_report", "$ratio ± $ratio_err, $(ratio/ratio_err) standard deviation significance")

rmprocs(worker_pool)

@info "Job done!"


text:
#
# Categorical Sampler
#
# parameters:
# (1)-(n)  0 <= w_i <= 1
#

type CategoricalSampler <: DiscreteDistributionSampler
	numcategories::Int
	paramranges::Vector{Tuple{Float64,Float64}}
	distribution::Categorical
	function CategoricalSampler(numcategories::Int, params=Float64[])
		numcategories >= 1 || error("number of categories must be at least one")
		s = new(numcategories, fill((0.0,1.0), numcategories))
		# Note: could parameterise using number of parameters that is one less than the domain size
		# but there is not an obvious way to do this in a symmetrical way so as not bias any optimisation
		# of the parameters
		setparams!(s, isempty(params) ? fill(1.0, numcategories) : params)
		s
	end
end

function setparams!(s::CategoricalSampler, params)
	checkparamranges(s, params)
	totalweight = sum(params)
	if totalweight == 0.0
		weights = fill(1.0/s.numcategories, s.numcategories)
	else
		weights = params ./ totalweight
	end
	# TODO we assume that sum of params is close enough (even with rounding errors) to 1.0 to satisfy the Categorical constructor
	# so we do not need handle the rounding errors
	s.distribution = Categorical(weights)
end

getparams(s::CategoricalSampler) = copy(s.distribution.p)

function estimateparams!(s::CategoricalSampler, traces)
	samples = extractsamplesfromtraces(s, traces)
	if length(samples) >= 1
		s.distribution = fit(primarydistributiontype(s.distribution), samples)
		# if there is a sequence of values at the end of the categories with zero probabilities,
		# then the fit will reduce the number of categories which will cause get parameters to return too
		# few values - we fix this here by padding with zeros
		if length(s.distribution.p) < s.numcategories
			s.distribution = Categorical([copy(s.distribution.p); fill(0.0, s.numcategories - length(s.distribution.p))])
		end
	end
end


text:
module ConvolutionTools

import Base.conv

export conv, conv_valid, pad_const

"""
    pad_const(a, n[, c])

Pads any multidimensional array with `n[i]` constant values `c` on both sides.
Returns an array of size `size(a) + 2*n`.
"""
function pad_const(a::Array, n::Array{Int, 1}, c=0)
    shape = collect(size(a)) .+ 2 .* n
    b = zeros(eltype(a), shape...) .+ c
    block = [(n[i] + 1):(size(a)[i] + n[i]) for i in 1:length(n)]
    b[block...] = a
    return b
end

function conv{T<:Base.LinAlg.BlasFloat, N}(u::Array{T, N}, v::Array{T, N})
    su = [size(u)...]
    sv = [size(v)...]
    sp = su .+ sv .- 1 # for zero-padding
    # pad for efficient FFT
    sp2 = [n > 1024 ? nextprod([2,3,5], n) : nextpow2(n) for n in sp]
    upad = zeros(T, (sp2...))
    upad[[1:n for n in su]...] = u
    vpad = zeros(T, (sp2...))
    vpad[[1:n for n in sv]...] = v
    if T <: Real
        p = plan_rfft(upad)
    else
        p = plan_fft(upad)
    end
    # perform convolution
    y = p \ ((p * upad) .* (p * vpad))
    # undo FFT padding
    full = [1:sp[i] for i in 1:N]
    return y[full...]
end
conv{T<:Integer, N}(u::Array{T, N}, v::Array{T, N}) = round(Int, conv(float(u), float(v)))
conv{T<:Integer, S<:Base.LinAlg.BlasFloat, N}(u::Array{T, N}, v::Array{S, N}) = conv(float(u), v)
conv{T<:Integer, S<:Base.LinAlg.BlasFloat, N}(u::Array{S, N}, v::Array{T, N}) = conv(u, float(v))

function strip_invalid{T, S, N}(output::Array{T, N}, stencil::Array{S, N})
    shape_padded = [size(output)...]
    shape_stencil = [size(stencil)...]
    shape_valid = shape_padded .- 2 .* (shape_stencil .- 1)
    invalid = shape_stencil .- 1
    lower = invalid .+ 1
    upper = invalid .+ shape_valid
    center = [lower[i]:upper[i] for i in 1:N]
    return output[center...]
end

function conv_valid(data, stencil)
    full = conv(data, stencil)
    return strip_invalid(full, stencil)
end

end # module


text:
export Sphere

mutable struct Sphere{N,N²,D,T} <: Shape{N,N²,D,T}
    c::SVector{N,T}  # center of sphere
    r::T  # radius
    data::D  # auxiliary data
    # Sphere{N,N²,D,T}(c,r,data) where {N,N²,D,T<:Real} = new(c,r,data)  # suppress default outer constructor
end

Sphere(c::SVector{N,T}, r::T, data::D=nothing) where {N,D,T<:Real} = Sphere{N,N*N,D,T}(c, r, data)
Sphere(c::AbstractVector{<:Real}, r::Real, data=nothing) = (N = length(c); Sphere(SVector{N}(c), r, data))

Base.:(==)(s1::Sphere, s2::Sphere) = s1.c==s2.c && s1.r==s2.r && s1.data==s2.data
Base.isapprox(s1::Sphere, s2::Sphere) = s1.c≈s2.c && s1.r≈s2.r && s1.data==s2.data
Base.hash(s::Sphere, h::UInt) = hash(s.c, hash(s.r, hash(s.data, hash(:Sphere, h))))

Base.in(x::SVector{N,<:Real}, s::Sphere{N}) where {N} = sum(abs2, x - s.c) ≤ s.r^2

function surfpt_nearby(x::SVector{N,<:Real}, s::Sphere{N}) where {N}
    nout = x==s.c ? SVector(ntuple(k -> k==1 ? 1.0 : 0.0, Val(N))) :  # nout = e₁ for x == s.c
                    normalize(x-s.c)
    return s.c+s.r*nout, nout
end

translate(s::Sphere{N,N²,D}, ∆::SVector{N}) where {N,N²,D,T<:Real} = Sphere{N,N²,D,T}(s.c+∆, s.r, s.data)

bounds(s::Sphere) = (s.c.-s.r, s.c.+s.r)


text:
struct BinaryRelation{O, S, T}
    x::S
    y::T
end

BinaryRelation{O}(x::X, y::Y) where {O,X,Y} = BinaryRelation{O,X,Y}(x, y)

x = DiscreteVariable(:x, 1:5)
y = DiscreteVariable(:y, 2:6)

# BinaryRelation(==, x, y)
# BinaryRelation(==, x, 1)

# Use NodeVariables to break up binary relations like x + y ~ 3 and ternary such as x + y <= z
struct NodeVariable <: Var
    name
    op
    x
    y
    domain
    booleans
    varmap
end

Base.show(io::IO, var::NodeVariable) = print(io, "$(var.name) = $(var.op)($(var.x), $(var.y))")

find_domain(op, x, y) = sort(collect(Set([op(i, j) for i in domain(x) for j in domain(y)])))

domain(v::Var) = v.domain
booleans(v::Var) = v.booleans

domain(x::Real) = x:x

"All booleans inside an expression that must be added to the encoding"
recursive_booleans(v::NodeVariable) = v.booleans ∪ booleans(v.x) ∪ booleans(v.y)

function NodeVariable(op, x, y, name=gensym())
    domain = find_domain(op, x, y)
    booleans = [Variable(name, i) for i ∈ indices(domain)]
    varmap = Dict(i => v for (i, v) in zip(domain, booleans))
    NodeVariable(name, op, x, y, domain, booleans, varmap)
end

Base.:+(x::Var, y) = NodeVariable(+, x, y)
Base.:+(x, y::Var) = NodeVariable(+, x, y)

Base.:*(x::Var, y) = NodeVariable(*, x, y)
Base.:*(x, y::Var) = NodeVariable(*, x, y)

Base.:^(x::Var, y) = NodeVariable(^, x, y)
Base.:^(x, y::Var) = NodeVariable(^, x, y)

function clauses(var::NodeVariable)

    op = var.op
    x = var.x
    y = var.y

    clauses = exactly_one(var.booleans)

    # deal with 2 + X
    if x isa Real
        for j in domain(var.y)
            value = op(x, j)
            # y == j => var == op(x, j)
            push!(clauses, ¬(y.varmap[j]) ∨ var.varmap[value])

        end

    # deal with x + 2
    elseif y isa Real
        for i in domain(var.x)
            value = op(i, y)
            # x == i => var == op(i, y)
            push!(clauses, ¬(x.varmap[i]) ∨ var.varmap[value])

        end

    else

        # neither x nor y is real
        for i in domain(var.x)
            for j in domain(var.y)
                value = op(i, j)
                # x == i && y == j => var == op(i, j)
                push!(clauses, ¬(x.varmap[i]) ∨ ¬(y.varmap[j]) ∨ var.varmap[value])
            end
        end

    end

    return clauses

end

Base.:(==)(v::Var, w) = BinaryRelation{==}(v, w)
Base.:(<=)(v::Var, w) = BinaryRelation{<=}(v, w)

"Encode relation like x == 1"
function encode(rel::BinaryRelation{==, <:Var, <:Real})

    x = rel.x
    y = rel.y

    if y ∉ domain(x)
        error("$y is not in the domain of $x")
    end

    boolean = x.varmap[y]

    return [boolean]
end

"Encode relation like x != 1"
function encode(rel::BinaryRelation{!=, <:Var, <:Real})

    x = rel.x
    y = rel.y

    if y ∈ domain(x)
        boolean = ¬(x.varmap[y])
    end

    return [boolean]
end

"Encode relation like x == y"
function encode(rel::BinaryRelation{==, <:Var, <:Var})
    x = rel.x
    y = rel.y

    clauses = []

    for i in domain(x)
        if i in domain(y)
        # (x == i) => (y == i)
            push!(clauses, ¬(x.varmap[i]) ∨ (y.varmap[i]))

        else
            push!(clauses, ¬(x.varmap[i]))
        end
    end

    for i in domain(y)
        if i in domain(x)
        # (y == i) => (x == i)
            push!(clauses, ¬(y.varmap[i]) ∨ (x.varmap[i]))

        else
            push!(clauses, ¬(y.varmap[i]))
        end
    end

    return clauses
end

function encode(rel::BinaryRelation{!=, <:Var, <:Var})
    x = rel.x
    y = rel.y

    clauses = []

    for i in domain(x)
        if i in domain(y)
        # (x == i) => (y != i)
            push!(clauses, ¬(x.varmap[i]) ∨ ¬(y.varmap[i]))
        end
    end

    return clauses
end

function parse_expression(varmap, ex)
    op = operation(ex)
    args = arguments(ex)

    new_args = parse_expression.(Ref(varmap), args)

    return op, new_args
end

parse_expression(varmap, ex::Term) = varmap[ex]
parse_expression(varmap, ex::Sym) = varmap[ex]
parse_expression(varmap, ex::Num) = parse_expression(varmap, value(ex))
parse_expression(varmap, ex::Real) = ex

"Parse a symbolic expression into a relation"
function parse_relation(varmap, ex)
    op = operation(ex)
    args = arguments(ex)

    lhs = args[1]
    rhs = args[2]

    return BinaryRelation{op}(parse_expression(varmap, lhs), parse_expression(varmap, rhs))

end

# function parse_relation!(varmap, ex::Equation)

#     varmap[ex.lhs] = parse_expression(varmap, ex.rhs)

# end

is_variable(ex) = !istree(ex) || (operation(ex) == getindex)

function process(constraint)

    new_constraints = []

    constraint2 = value(constraint)

    op = operation(constraint2)
    args = arguments(constraint2)

    lhs = args[1]
    rhs = args[2]

    intermediates_generated = false

    if !is_variable(lhs)
        lhs = ReversePropagation.cse_equations(lhs)

        append!(new_constraints, lhs[1])
        lhs = lhs[2]

    end

    if !is_variable(rhs)
        rhs = ReversePropagation.cse_equations(rhs)

        append!(new_constraints, rhs[1])
        rhs = rhs[2]

    end

    # @show new_constraints
    # @show op, lhs, rhs
    push!(new_constraints, op(lhs, rhs))

    return new_constraints

end

process(constraints::Vector) = reduce(vcat, process(constraint) for constraint in constraints)

function parse_constraint!(domains, ex)

    expr = value(ex)
    op = operation(expr)

    new_constraints = []

    if op == ∈   # assumes right-hand side is an explicit set specifying the doomain

        args = arguments(expr)

        var = args[1]
        domain = args[2]

        domains[var] = domain

    else
        new_constraints = process(expr)

        println()
        # @show expr
        # @show new_constraints

    end

    return domains, new_constraints

end

function parse_constraints(constraints)
    # domains = Dict(var => -Inf..Inf for var in vars)

    additional_vars = []
    domains = Dict()
    all_new_constraints = []  # excluding domain constraints

    for constraint in constraints

        # binarize constraints:
        domains, new_constraints = parse_constraint!(domains, value(constraint))

        for statement in new_constraints
            if statement isa Assignment
                push!(additional_vars, statement.lhs)
            end
        end

        append!(all_new_constraints, new_constraints)
    end

    # append!(additional_vars, keys(domains))

    return domains, all_new_constraints, additional_vars
end

Base.isless(x::Sym, y::Sym) = isless(x.name, y.name)
Base.isless(x::Term, y::Term) = isless(string(x), string(y))
Base.isless(x::Term, y::Sym) = isless(string(x), string(y))
Base.isless(x::Sym, y::Term) = isless(string(x), string(y))

struct ConstraintSatisfactionProblem
    original_vars
    additional_vars   # from binarizing
    domains
    constraints
end

function ConstraintSatisfactionProblem(constraints)

    domains, new_constraints, additional_vars = parse_constraints(constraints)
    @show keys(domains)
    vars = sort(identity.(keys(domains)))
    additional_vars = sort(identity.(additional_vars))

    return ConstraintSatisfactionProblem(vars, additional_vars, domains, new_constraints)
end

struct DiscreteCSP
    original_vars
    additional_vars
    constraints
    varmap   # maps symbolic variables to the corresponding DiscreteVariable
end

DiscreteCSP(constraints) = DiscreteCSP(ConstraintSatisfactionProblem(constraints))

function DiscreteCSP(prob::ConstraintSatisfactionProblem)
    variables = []
    varmap = Dict()

    new_constraints = []

    for (var, domain) in prob.domains
        variable = DiscreteVariable(var, domain)
        push!(variables, variable)
        push!(varmap, var => variable)
    end

    for constraint in prob.constraints
        constraint = value(constraint)
        if constraint isa Assignment

            lhs = constraint.lhs

            op, new_args = parse_expression(varmap, constraint.rhs)   # makes a NodeVariable

            # @show op, new_args

            variable = NodeVariable(op, new_args[1], new_args[2], lhs)

            push!(variables, variable)
            push!(varmap, lhs => variable)

        else
            # @show constraint
            push!(new_constraints, parse_relation(varmap, constraint))
        end
    end

    original_vars = [varmap[var] for var in prob.original_vars]
    additional_vars = [varmap[var] for var in prob.additional_vars]

    return DiscreteCSP(original_vars, additional_vars, new_constraints, varmap)
end

function encode(prob::DiscreteCSP)
    all_variables = []
    all_clauses = []

    variables = Any[prob.original_vars; prob.additional_vars]

    domains = Dict(v.name => v.domain for v in variables)

    for var in variables
        append!(all_variables, var.booleans)
        append!(all_clauses, clauses(var))
    end

    # @show prob.constraints

    for constraint in prob.constraints
        # @show constraint
        append!(all_clauses, encode(constraint))
    end

    # @show all_variables
    # @show all_clauses

    return SymbolicSATProblem(identity.(all_variables), identity.(all_clauses))
end

function solve(prob::DiscreteCSP)
    symbolic_sat_problem = encode(prob)

    status, result_dict = solve(symbolic_sat_problem)

    (status == :unsat) && return status, missing

    # variables = Any[prob.original_vars; prob.additional_vars]
    variables = prob.original_vars

    return status, decode(prob, result_dict)
end

decode(prob::DiscreteCSP, result_dict) = Dict(v.name => decode(result_dict, v) for v in prob.original_vars)

function all_solutions(prob::DiscreteCSP)
    sat_problem = encode(prob)
    # variables = Any[prob.original_vars; prob.additional_vars]
    variables = prob.original_vars

    prob2 = encode(prob)
    prob3 = prob2.p   # SATProblem

    sat_solutions = all_solutions(prob3)

    if isempty(sat_solutions)
        return sat_solutions
    end

    solutions = []

    for solution in sat_solutions
        push!(solutions, decode(prob, decode(prob2, solution)))
    end

    return solutions

end


text:
struct Simulator{N,T<:AbstractFloat,W<:Wave{N,T},R<:AbstractResource,H<:Tuple}
    wave::W
    resource::R
    hooks::H
    duration::T
    maxiter::Int
end

function Simulator(wave::Wave{N,T}, hooks...; resource=CPU1(), duration=0.01) where {N,T}
    maxiter = ceil(Int, duration / wave.dt)
    Simulator{N,T,typeof(wave),typeof(resource),typeof(hooks)}(wave, resource, hooks, T(duration), maxiter)
end

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

function simulator_init(state::State, sim::Simulator)
    foreach(hook -> hook_init!(hook, state, sim), sim.hooks)
end

function update!(state::State, backend::Backend, sim::Simulator)
    backend_update!(state, backend, sim)
    state.previous, state.current = state.current, state.previous
    state.t    += sim.wave.dt
    state.iter += 1
    foreach(hook -> hook_update!(hook, state, sim), sim.hooks)
end

function simulate!(state::State, backend::Backend, sim::Simulator)
    simulator_init(state, sim)
    for _ in 1:sim.maxiter
        update!(state, backend, sim)
    end
    state
end

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

function simulate(f0, sim::Simulator{N,T}, domain::Domain{N}) where {N,T}
    backend = backend_init(sim.resource, domain, sim)
    state   = state_init(f0, backend, domain, sim)
    simulate!(state, backend, sim)
    backend_cleanup!(backend, state, sim)
    state
end

function simulate(sim::Simulator{N,T}, domain::Domain{N}) where {N,T}
    backend = backend_init(sim.resource, domain, sim)
    state   = state_init(backend, domain, sim)
    simulate!(state, backend, sim)
    backend_cleanup!(backend, state, sim)
    state
end

function simulate_gauss(sim::Simulator{N,T}, domain::Domain{N}) where {N,T}
    simulate(gauss(domain), sim, domain)
end


text:
using CSV, DataFrames, Plots, LinearAlgebra, Interact, StatsBase, Statistics

f(x) = typeof(x) <: Number ? Float32(x) : parse(Float32, x)

#Load data, clean data, recast data
kaggleData = CSV.read("/Users/javier/Desktop/SynData/synthetic-derivative/Dataset/kag_risk_factors_cervical_cancer.csv", DataFrame)
data = Array(kaggleData)
data[data .== "?"] .= -1
data = f.(data)

#Load data, clean data, recast Syn data
synData = CSV.read("/Users/javier/Desktop/SynData/synthetic-derivative/Dataset/synData-3.csv", DataFrame)
synData = CSV.read("/Users/javier/Desktop/SynData/synthetic-derivative/Dataset/synData-4.csv", DataFrame)
dataSyn = Array{Int}(synData)
# data[data .== "?"] .= -1
dataSyn = Array{Int}(f.(dataSyn[1:5000,:]))

#PC
function doPCA(arr)
    l = size(arr,1)
    # mean2=[arr[j,i]-sum(arr[:,i])/(l) for j=1:l,i=1:size(arr,2)]
    mean2=[arr[j,i]-mean(arr[:,i]) for j=1:l,i=1:size(arr,2)]

    #Cov matrix
    m=transpose(mean2)*mean2/l
    # plot(m, seriestype=:heatmap)

    #Eigenvalues and Eigenvectors
    d,v = eigen(m)
    dNorm=d/maximum(d)

    m, dNorm, v
end

mReal, dReal, vReal = doPCA(data)
mSyn, dSyn, vSyn = doPCA(dataSyn)
mean(abs2, mReal .- mSyn)
#Plot eigenvalues
plot(dReal, dSyn, ms=10, markershapes = [:circle], markerstrokewidth=0, lw=1, size=(500,500))
plot!(0:1,x->x, color=:red)

#Plot eigenvectors
@manipulate for i in 1:size(vReal,2)
    plot(abs.(vReal[:,i]), abs.(vSyn[:,i]), frame=:box, xrotation=0, ms=5, markershapes = [:circle], markerstrokewidth=0, lw=0, size=(500,500), tickfont = font(12, "Helvetica"), legendfont = font(12, "Helvetica"), label="$(dReal[i])", legend=(0.15,0.9), xlim=(-0.1,1), ylim=(-0.1,1), rightmargin = 5Plots.mm)
    plot!(0:1, x->x)
end

#Plot histogram
@manipulate for i in 1:size(dataSyn,2)
    lab = names(kaggleData)[i]
    plot(Int.(dataSyn[:,i]), seriestype=:barhist, frame=:box, width=0, opacity=0.5, normalize=true, title="$lab", label="Synth", tickfont = font(12, "Helvetica"), xlabel="Values", ylabel="Freq")
    plot!(data[:,i], seriestype=:barhist, width=0, opacity=0.5, normalize=true, label="Real")
end
typeof(dataSyn)
#=Plot eigenvector inner product as heatmap
    A = ∑|λᵢ⟩⟨cᵢ| and B = ∑|β⟩⟨c|
    Then, C = Aᵀ⋅B = ∑⟨λᵢ|βⱼ⟩|cᵢ⟩⟨cᵢ|
    Notice that ⟨λᵢ|βⱼ⟩ is the correlation between the eigenvector of the synthetic
    data covariance matrix (SDCM) and the real data covariance matrix (RDCM). Hence, a heatmap
    of C shows the correlation between all eigenvectors from the SDCM with the RDCM
    one-on-one.
=#
plot(vSyn' * vSyn, seriestype=:heatmap, xticks=(0:35,reverse(names(kaggleData))), yticks=(0:35,reverse(names(kaggleData))), xrotation=80, bottommargin = 30Plots.mm,  size=(700,700))
plot(vReal' * vSyn, seriestype=:heatmap, xticks=(0:35,reverse(names(kaggleData))), yticks=(0:35,reverse(names(kaggleData))), xrotation=80, bottommargin = 30Plots.mm,  size=(700,700))
plot(abs.(vReal' * vSyn), seriestype=:heatmap, xticks=(0:35,reverse(names(kaggleData))), yticks=(0:35,reverse(names(kaggleData))), xrotation=80, bottommargin = 30Plots.mm,  size=(700,700))

#=
    From the previous, we now compute the eigenvectors of C. Notice that if C corresponded to a completely correlated system, i.e., ⟨λᵢ|βⱼ⟩ = δᵢⱼ, then C = Identity. In that case, the spectrum is degenerate and
    the eigenvalues all equal 1. In our case, notice the eigenvalues are complex, yet the absolute value
    are all close to 1.
=#
dCorr, vCorr = eigen(vReal' * vSyn)
dCorr, vCorr = eigen(vReal' * vReal .- 0.001*randn(36,36))
dCorr, vCorr = eigen(vReal' * randn(36,36))
dCorr, vCorr = eigen(vReal' * vSyn[sample(collect(1:36),36, replace=false),:])
dCorr, vCorr = eigen(zeros(36,36))

plot(abs.(dCorr))

@manipulate for i in 1:size(vReal,2)
    plot(abs.(vReal[:,i]), frame=:box, xrotation=80, ms=10, markershapes = [:circle], markerstrokewidth=0, lw=1, size=(500,500), tickfont = font(12, "Helvetica"), legendfont = font(12, "Helvetica"), label="$(dReal[i])", legend=(0.75,0.9), xticks=(0:35,names(kaggleData)), bottommargin = 40Plots.mm)
end

# using HTTP, JSON, BSON, Dates
#
# a = HTTP.request("GET", "https://data-live.flightradar24.com/zones/fcgi/feed.js")
#
# a.body
#
# aa = JSON.parse(String(a.body))

#TO EXPORT
PATH = "/Users/javier/Dropbox/Aplicaciones/Overleaf/SynthData/Figs/"

begin
    fig = plot([abs.(vReal[:,i]) for i in 31:36], [abs.(vSyn[:,i]) for i in 31:36], frame=:box, xrotation=0, ms=5, markershapes = [:circle], markerstrokewidth=0, lw=0, size=(500,500), tickfont = font(12, "Helvetica"), legendfont = font(12, "Helvetica"), label="", legend=(0.15,0.9), xlim=(-0.1,1), ylim=(-0.1,1), rightmargin = 5Plots.mm, xlabel="Real Data Covariance Matrix Eigenvectors", ylabel="Synth Data Covariance Matrix Eigenvectors", title="5 largest eigenvalues")
    fig = plot!(0:1, x->x, legend=:none, lw=2, c=:black)
    savefig(fig, PATH * "SR-DCM-EigenVectors.png")
end

begin
    fig = plot(dReal, dSyn, ms=10, markershapes = [:circle], markerstrokewidth=0, lw=1, size=(500,500), tickfont = font(12, "Helvetica"), frame=:box, legend=:topleft, label="Eigenvalues", legendfont = font(12, "Helvetica"), xlabel="Real Data Covariance Matrix", ylabel="Synth Data Covariance Matrix")
    fig = plot!(0:1,x->x, color=:black, lw=2, label="y=x")
    savefig(fig, PATH * "SR-DCM-EigenValues.png")
end

s(x) = split(x, " ")[1]
begin
    for i in 1:size(dataSyn,2)
        lab = names(kaggleData)[i]
        name = "$i-" * s.(names(kaggleData))[i]
        fig = plot(Int.(dataSyn[:,i]), seriestype=:barhist, frame=:box, width=0, opacity=0.5, normalize=true, title="$lab", label="Synth", tickfont = font(12, "Helvetica"), xlabel="Values", ylabel="Freq")
        fig = plot!(data[:,i], seriestype=:barhist, width=0, opacity=0.5, normalize=true, label="Real")
        savefig(fig, PATH * "Hist-$name.png")
    end
end

begin
    fig = plot(abs.(vReal' * vSyn), seriestype=:heatmap, xticks=(0:35,reverse(names(kaggleData))), yticks=(0:35,reverse(names(kaggleData))), xrotation=80, bottommargin = 30Plots.mm,  size=(700,700))
    savefig(fig, PATH * "SandR-DCMInnerProduct.png")
    fig = plot(vSyn' * vSyn, seriestype=:heatmap, xticks=(0:35,reverse(names(kaggleData))), yticks=(0:35,reverse(names(kaggleData))), xrotation=80, bottommargin = 30Plots.mm,  size=(700,700))
    savefig(fig, PATH * "SandS-DCMInnerProduct.png")
    fig = plot(abs.(vReal' * (vReal .- 0.4*randn(36,36))), seriestype=:heatmap, xticks=(0:35,reverse(names(kaggleData))), yticks=(0:35,reverse(names(kaggleData))), xrotation=80, bottommargin = 30Plots.mm,  size=(700,700))
    savefig(fig, PATH * "R-DCMInnerProductwNoise.png")
end

begin
    dCorr, vCorr = eigen(vReal' * vSyn)
    fig = plot(sort(abs.(dCorr), rev=true), markershapes = [:circle], markerstrokewidth=0, lw=1,
    frame=:box, legend=:none, xlabel="Eigenvalues", ylabel="Amplitude")
    savefig(fig, PATH * "eigenvaluesForSRDCMInnerproduct.png")
    dCorr, vCorr = eigen(vReal' * vReal .- 0.001*randn(36,36))
    fig = plot(sort(abs.(dCorr), rev=true), markershapes = [:circle], markerstrokewidth=0, lw=1,
    frame=:box, legend=:none, xlabel="Eigenvalues", ylabel="Amplitude")
    savefig(fig, PATH * "eigenvaluesForRDCMInnerproductwNoise.png")
    dCorr, vCorr = eigen(vReal' * vReal .- 0.000*randn(36,36))
    fig = plot(sort(abs.(dCorr), rev=true), markershapes = [:circle], markerstrokewidth=0, lw=1,
    frame=:box, legend=:none, xlabel="Eigenvalues", ylabel="Amplitude")
    savefig(fig, PATH * "eigenvaluesForRDCMInnerproduct.png")
    dCorr, vCorr = eigen(ones(36,36))
    fig = plot(sort(abs.(dCorr), rev=true), markershapes = [:circle], markerstrokewidth=0, lw=1,
    frame=:box, legend=:none, xlabel="Eigenvalues", ylabel="Amplitude")
    savefig(fig, PATH * "eigenvaluesForOnes.png")
    # dCorr, vCorr = eigen(vReal' * randn(36,36))
    # dCorr, vCorr = eigen(vReal' * vSyn[sample(collect(1:36),36, replace=false),:])

end
#################################
#PRIVACY
##################################
sample(1:10, 11, replace=false)
function replaceRow(ar1, row, idx)
    ar3 = copy(ar1)
    ar3[idx,:] = row
    return ar3
end

function replaceRows(ar1,ar2, nrows)
    samp1 = sample(1:size(ar1,1),nrows, replace=false)
    samp2 = sample(1:size(ar2,1),nrows, replace=false)
    # err = Vector{Float32}(undef, size(ar1,1))
    temp = copy(ar1)
    for i in 1:nrows
        temp = replaceRow(temp, ar2[samp2[i],:],samp1[i])
        # err[i] = mean(abs, ar1 .- temp)
    end
    temp
end

data1 = replaceRow(data, dataSyn[1,:],1)

err = [sum(abs, data .- replaceRows(data,dataSyn, j)) for i in 1:500, j in 0:100]
reshape(mean(err, dims=1),:)
plot(reshape(mean(err, dims=1),:), ribbon=reshape(std(err, dims=1),:))

mReal, dReal, vReal = doPCA(data)
mSyn, dSyn, vSyn = doPCA(replaceRows(data,dataSyn, 100))

@manipulate for nrows in 1:858
    mSyn, dSyn, vSyn = doPCA(replaceRows(data,dataSyn, nrows))
    plot(dReal, dSyn, ms=10, markershapes = [:circle], markerstrokewidth=0, lw=1, size=(500,500))
    plot!(0:1,x->x, color=:red)
end

max_stars_repo_name: jquetzalcoatl/SynthData

max_issues_repo_name: jquetzalcoatl/SynthData

max_forks_repo_name: jquetzalcoatl/SynthData


text:
# This takes various figures from the following paper
#
# Lippitsch, R., et al. 2003. Upper mantle structure beneath the Alpine orogen from high-resolution teleseismic tomography. J. Geophys. Res. 108, 2376. https://doi.org/10.1029/2002JB002016
#
# and transfers it into paraview format
#
# For convenience we created a zip file with all cross-sections and mapviews here:
#
#

using GeophysicalModelGenerator

# Process cross-sections of Figure 13. Note that we estimated some of the lon/lat locations

data_Fig13a         =   Screenshot_To_GeoData("Lippitsch_Fig13a.png",( 4.65,45.73, -400.0), (17.23, 43.80, 0.0))
Write_Paraview(data_Fig13a, "Lippitsch_Fig13a")

data_Fig13b         =   Screenshot_To_GeoData("Lippitsch_Fig13b.png",( 5.51,51.53, -400.0), (12.04, 43.68 , 0.0))
Write_Paraview(data_Fig13b, "Lippitsch_Fig13b")

data_Fig13c         =   Screenshot_To_GeoData("Lippitsch_Fig13c.png",(17.78,50.95, -400.0), (11.66, 43.68, 0.0))
Write_Paraview(data_Fig13c, "Lippitsch_Fig13c")

# Mapview images
Corner_LowerLeft    =   ( 3.5, 43.0 , -150.0)
Corner_UpperRight   =   (15.5, 50.0 , -150.0)
Corner_LowerRight   =   (15.5, 43.0 , -150.0)
Corner_UpperLeft    =   (3.5 , 50.0 , -150.0)
data_Fig13_map      =   Screenshot_To_GeoData("Fig13_mapview.png",Corner_LowerLeft, Corner_UpperRight, Corner_LowerRight=Corner_LowerRight,Corner_UpperLeft=Corner_UpperLeft)
Write_Paraview(data_Fig13_map, "Lippitsch_Fig13_mapview")

Depth                   =   -90.0;
Corner_LowerLeft        =   (Corner_LowerLeft[1],   Corner_LowerLeft[2],    Depth)
Corner_UpperRight       =   (Corner_UpperRight[1],  Corner_UpperRight[2],   Depth)
Corner_LowerRight       =   (Corner_LowerRight[1],  Corner_LowerRight[2],   Depth)
Corner_UpperLeft        =   (Corner_UpperLeft[1],   Corner_UpperLeft[2],    Depth)
data_Fig12_90km    =   Screenshot_To_GeoData("Fig12_90km.png",Corner_LowerLeft, Corner_UpperRight, Corner_LowerRight=Corner_LowerRight,Corner_UpperLeft=Corner_UpperLeft)
Write_Paraview(data_Fig12_90km, "Lippitsch_Fig12_90km")

Depth                   = -180.0;
Corner_LowerLeft        =   (Corner_LowerLeft[1],   Corner_LowerLeft[2],    Depth)
Corner_UpperRight       =   (Corner_UpperRight[1],  Corner_UpperRight[2],   Depth)
Corner_LowerRight       =   (Corner_LowerRight[1],  Corner_LowerRight[2],   Depth)
Corner_UpperLeft        =   (Corner_UpperLeft[1],   Corner_UpperLeft[2],    Depth)
data_Fig12_180km   =   Screenshot_To_GeoData("Fig12_180km.png",Corner_LowerLeft, Corner_UpperRight, Corner_LowerRight=Corner_LowerRight,Corner_UpperLeft=Corner_UpperLeft)
Write_Paraview(data_Fig12_180km, "Lippitsch_Fig12_180km")

Depth                   = -300.0;
Corner_LowerLeft        =   (Corner_LowerLeft[1],   Corner_LowerLeft[2],    Depth)
Corner_UpperRight       =   (Corner_UpperRight[1],  Corner_UpperRight[2],   Depth)
Corner_LowerRight       =   (Corner_LowerRight[1],  Corner_LowerRight[2],   Depth)
Corner_UpperLeft        =   (Corner_UpperLeft[1],   Corner_UpperLeft[2],    Depth)
data_Fig12_300km   =   Screenshot_To_GeoData("Fig12_300km.png",Corner_LowerLeft, Corner_UpperRight, Corner_LowerRight=Corner_LowerRight,Corner_UpperLeft=Corner_UpperLeft)
Write_Paraview(data_Fig12_300km, "Lippitsch_Fig12_300km")

Depth                   = -400.0;
Corner_LowerLeft        =   (Corner_LowerLeft[1],   Corner_LowerLeft[2],    Depth)
Corner_UpperRight       =   (Corner_UpperRight[1],  Corner_UpperRight[2],   Depth)
Corner_LowerRight       =   (Corner_LowerRight[1],  Corner_LowerRight[2],   Depth)
Corner_UpperLeft        =   (Corner_UpperLeft[1],   Corner_UpperLeft[2],    Depth)
data_Fig12_400km   =   Screenshot_To_GeoData("Fig12_400km.png",Corner_LowerLeft, Corner_UpperRight, Corner_LowerRight=Corner_LowerRight,Corner_UpperLeft=Corner_UpperLeft)
Write_Paraview(data_Fig12_400km, "Lippitsch_Fig12_400km")

# Example of how we can save this to a multibock *.VTM file (which allows you to open all files @ once in paraview)
vtmfile = vtk_multiblock("Lippitsch_CrossSections")
Write_Paraview(data_Fig12_90km, vtmfile)
Write_Paraview(data_Fig12_180km, vtmfile)
Write_Paraview(data_Fig12_300km, vtmfile)
Write_Paraview(data_Fig12_400km, vtmfile)
vtk_save(vtmfile)


text:
using PlotAxes

function PlotAxes.asplotable(x::AxisArray{<:SourceTracking},args...;component=1,kwds...)
  vals = cat((xi[:,:,component,:] for xi in x)...,dims=4)
  axs = AxisArrays.axes(x[1])
  start,stop = extrema(log.(ustrip.(freqs(x[1]))))
  logfreq = Axis{:logfreq}(range(start,stop=stop,length=length(freqs(x[1]))))
  y = AxisArray(vals,axs[1],logfreq,axs[4],AxisArrays.axes(x)...)
  asplotable(y,:time,:logfreq,:scale,:params,args...;kwds...)
end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: pi
using AWS.Compat
using AWS.UUIDs

"""
    DescribeDimensionKeys()

For a specific time period, retrieve the top N dimension keys for a metric.

# Required Parameters
- `EndTime`: The date and time specifying the end of the requested time series data. The value specified is exclusive - data points less than (but not equal to) EndTime will be returned. The value for EndTime must be later than the value for StartTime.
- `GroupBy`: A specification for how to aggregate the data points from a query result. You must specify a valid dimension group. Performance Insights will return all of the dimensions within that group, unless you provide the names of specific dimensions within that group. You can also request that Performance Insights return a limited number of values for a dimension.
- `Identifier`: An immutable, AWS Region-unique identifier for a data source. Performance Insights gathers metrics from this data source. To use an Amazon RDS instance as a data source, you specify its DbiResourceId value - for example: db-FAIHNTYBKTGAUSUZQYPDS2GW4A
- `Metric`: The name of a Performance Insights metric to be measured. Valid values for Metric are:    db.load.avg - a scaled representation of the number of active sessions for the database engine.    db.sampledload.avg - the raw number of active sessions for the database engine.
- `ServiceType`: The AWS service for which Performance Insights will return metrics. The only valid value for ServiceType is: RDS
- `StartTime`: The date and time specifying the beginning of the requested time series data. You can't specify a StartTime that's earlier than 7 days ago. The value specified is inclusive - data points equal to or greater than StartTime will be returned. The value for StartTime must be earlier than the value for EndTime.

# Optional Parameters
- `Filter`: One or more filters to apply in the request. Restrictions:   Any number of filters by the same dimension, as specified in the GroupBy or Partition parameters.   A single filter for any other dimension in this dimension group.
- `MaxResults`: The maximum number of items to return in the response. If more items exist than the specified MaxRecords value, a pagination token is included in the response so that the remaining results can be retrieved.
- `NextToken`: An optional pagination token provided by a previous request. If this parameter is specified, the response includes only records beyond the token, up to the value specified by MaxRecords.
- `PartitionBy`: For each dimension specified in GroupBy, specify a secondary dimension to further subdivide the partition keys in the response.
- `PeriodInSeconds`: The granularity, in seconds, of the data points returned from Performance Insights. A period can be as short as one second, or as long as one day (86400 seconds). Valid values are:    1 (one second)    60 (one minute)    300 (five minutes)    3600 (one hour)    86400 (twenty-four hours)   If you don't specify PeriodInSeconds, then Performance Insights will choose a value for you, with a goal of returning roughly 100-200 data points in the response.
"""
describe_dimension_keys(EndTime, GroupBy, Identifier, Metric, ServiceType, StartTime; aws_config::AWSConfig=global_aws_config()) = pi("DescribeDimensionKeys", Dict{String, Any}("EndTime"=>EndTime, "GroupBy"=>GroupBy, "Identifier"=>Identifier, "Metric"=>Metric, "ServiceType"=>ServiceType, "StartTime"=>StartTime); aws_config=aws_config)
describe_dimension_keys(EndTime, GroupBy, Identifier, Metric, ServiceType, StartTime, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = pi("DescribeDimensionKeys", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("EndTime"=>EndTime, "GroupBy"=>GroupBy, "Identifier"=>Identifier, "Metric"=>Metric, "ServiceType"=>ServiceType, "StartTime"=>StartTime), args)); aws_config=aws_config)

"""
    GetResourceMetrics()

Retrieve Performance Insights metrics for a set of data sources, over a time period. You can provide specific dimension groups and dimensions, and provide aggregation and filtering criteria for each group.

# Required Parameters
- `EndTime`: The date and time specifiying the end of the requested time series data. The value specified is exclusive - data points less than (but not equal to) EndTime will be returned. The value for EndTime must be later than the value for StartTime.
- `Identifier`: An immutable, AWS Region-unique identifier for a data source. Performance Insights gathers metrics from this data source. To use an Amazon RDS instance as a data source, you specify its DbiResourceId value - for example: db-FAIHNTYBKTGAUSUZQYPDS2GW4A
- `MetricQueries`: An array of one or more queries to perform. Each query must specify a Performance Insights metric, and can optionally specify aggregation and filtering criteria.
- `ServiceType`: The AWS service for which Performance Insights will return metrics. The only valid value for ServiceType is: RDS
- `StartTime`: The date and time specifying the beginning of the requested time series data. You can't specify a StartTime that's earlier than 7 days ago. The value specified is inclusive - data points equal to or greater than StartTime will be returned. The value for StartTime must be earlier than the value for EndTime.

# Optional Parameters
- `MaxResults`: The maximum number of items to return in the response. If more items exist than the specified MaxRecords value, a pagination token is included in the response so that the remaining results can be retrieved.
- `NextToken`: An optional pagination token provided by a previous request. If this parameter is specified, the response includes only records beyond the token, up to the value specified by MaxRecords.
- `PeriodInSeconds`: The granularity, in seconds, of the data points returned from Performance Insights. A period can be as short as one second, or as long as one day (86400 seconds). Valid values are:    1 (one second)    60 (one minute)    300 (five minutes)    3600 (one hour)    86400 (twenty-four hours)   If you don't specify PeriodInSeconds, then Performance Insights will choose a value for you, with a goal of returning roughly 100-200 data points in the response.
"""
get_resource_metrics(EndTime, Identifier, MetricQueries, ServiceType, StartTime; aws_config::AWSConfig=global_aws_config()) = pi("GetResourceMetrics", Dict{String, Any}("EndTime"=>EndTime, "Identifier"=>Identifier, "MetricQueries"=>MetricQueries, "ServiceType"=>ServiceType, "StartTime"=>StartTime); aws_config=aws_config)
get_resource_metrics(EndTime, Identifier, MetricQueries, ServiceType, StartTime, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = pi("GetResourceMetrics", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("EndTime"=>EndTime, "Identifier"=>Identifier, "MetricQueries"=>MetricQueries, "ServiceType"=>ServiceType, "StartTime"=>StartTime), args)); aws_config=aws_config)


text:
isempty(x::Union{AbstractDict, AbstractArray, AbstractString, Tuple, NamedTuple}) = length(x) == 0
isempty(::Number) = false
isempty(::Nothing) = true
isempty(x) = false

getconvert(::Type{String}) = string
number(::Type{T}, x::T) where {T} = x
number(T, x::String) = Parsers.parse(T, x)
number(T, x) = convert(T, x)
getconvert(::Type{T}) where {T <: Number} = x -> number(T, x)
getconvert(T) = x -> convert(T, x)

const BUF = IOBuffer()
write(obj; kwargs...) = (write(BUF, obj; kwargs...); return String(take!(BUF)))

function write(io::IO, obj::Dict{K, V}; kwargs...) where {K, V}
    Base.write(io, '{')
    isempty(obj) && @goto done
    len = length(obj)
    i = 1
    for (k, v) in obj
        write(io, string(k); kwargs...)
        Base.write(io, ':')
        write(io, v; kwargs...) # recursive
        i < len && Base.write(io, ',')
        i += 1
    end
    @label done
    Base.write(io, '}')
    return
end

function write(io::IO, obj::Union{AbstractArray,Tuple,AbstractSet}; kwargs...)
    # always written as single array
    Base.write(io, '[')
    len = length(obj)
    i = 1
    for x in obj
        write(io, x; kwargs...) # recursive
        i < len && Base.write(io, ',')
        i += 1
    end
    Base.write(io, ']')
    return
end

write(io::IO, obj::Function; kwargs...) = (Base.write(io, obj.str); return)
write(io::IO, obj::Number; kwargs...) = (Base.write(io, string(obj)); return)
write(io::IO, obj::AbstractFloat; kwargs...) = (Base.print(io, isfinite(obj) ? obj : "null"); return)
write(io::IO, obj::Date; dateformat=Dates.ISODateFormat, kwargs...) = (Base.write(io, "\"$(Dates.format(obj, dateformat))\""); return)
write(io::IO, obj::DateTime; write_datetimeformat=nothing, dateformat=Dates.ISODateTimeFormat, kwargs...) = (Base.write(io, "\"$(Dates.format(obj, something(write_datetimeformat, dateformat)))\""); return)
write(io::IO, obj::Nothing; kwargs...) = (Base.write(io, "null"); return)
write(io::IO, obj::Missing; kwargs...) = (Base.write(io, "null"); return)
write(io::IO, obj::Bool; kwargs...) = (Base.write(io, obj ? "true" : "false"); return)
write(io::IO, obj::Union{Char, Symbol, Enum, Type}; kwargs...) = write(io, string(obj); kwargs...)
write(io::IO, p::Pair; kwargs...) = write(io, Dict(Symbol(p.first)=>p.second); kwargs...)

# N = # of fields
function generate_write_body(N, inds, names, omitempties, converts)
    # @show N, inds, names
    body = Expr(:block, :(kwargs = JSON2.mergedefaultkwargs(obj; kwargs...)))
    push!(body.args, :(Base.write(io, '{')))
    push!(body.args, :(j = 1))

    vals = ((Symbol("val_$i") for i = 1:N)...,)
    foreach(1:N) do i
        push!(body.args, quote
            nm = $(names[i])
            if isdefined(obj, $(inds[i]))
                $(vals[i]) = getfield(obj, $(inds[i]))
                if !($(omitempties[i])) || !JSON2.isempty($(vals[i]))
                    j == 1 || Base.write(io, ',')
                    JSON2.write(io, string(nm); kwargs...)
                    Base.write(io, ':')
                    JSON2.write(io, $(converts[i])($(vals[i])); kwargs...)
                    j += 1
                end
            else
                if !$(omitempties[i])
                    j == 1 || Base.write(io, ',')
                    JSON2.write(io, string(nm); kwargs...)
                    Base.write(io, ':')
                    JSON2.write(io, nothing; kwargs...)
                    j += 1
                end
            end
        end)
    end
    push!(body.args, :(Base.write(io, '}'); return))
    # @show body
    return body
end

@generated function write(io::IO, obj::T; kwargs...) where T
    N = fieldcount(T)
    inds = Tuple(1:N)
    names = Tuple(string(fieldname(T, i)) for i in inds)
    omitempties = Tuple(false for i in inds)
    converts = Tuple(getconvert(fieldtype(T, i)) for i in inds)
    ex = generate_write_body(N, inds, names, omitempties, converts)
    return ex
end


text:
## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-an-episode

"""
    episodes_get_single(episode_id; market="US")

**Summary**: Get Spotify catalog information for a single episode identified by its unique Spotify ID.

# Arguments
- `episode_id` : The Spotify ID for the episode_id

# Optional keywords
- `market` : An ISO 3166-1 alpha-2 country code. Default is set to "US".

# Example
```julia-repl
julia> Spotify.episodes_get_single("512ojhOuo1ktJprKbVcKyQ")[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 19 entries:
  :audio_preview_url    => "https://p.scdn.co/mp3-preview/566fcc94708f39bcddc09e4ce84a8e5db8f07d4d"
  :description          => "En ny tysk bok granskar för första gången Tredje rikets drogberoende, från Führerns k…
  :duration_ms          => 1502795
  :explicit             => false
  :external_urls        => {…
```
"""
function episodes_get_single(episode_id; market="US")
    return spotify_request("episodes/$episode_id?market=$market")
end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-multiple-episodes

"""
    episodes_get_multiple(episode_ids; market="US")

**Summary**: Get Spotify catalog information for several episodes based on their Spotify IDs.

# Arguments
- `episode_ids` : A comma-separated list of the Spotify IDs for the episodes. Maximum: 50 IDs.

# Optional keywords
- `market` : An ISO 3166-1 alpha-2 country code. Default is set to "US".

# Example
```julia-repl
julia> Spotify.episodes_get_multiple("77o6BIVlYM3msb4MMIL1jH,0Q86acNRm6V9GYx55SXKwf")[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 1 entry:
  :episodes => JSON3.Object[{…
```
"""
function episodes_get_multiple(episode_ids; market="US")
    return spotify_request("episodes?ids=$episode_ids&market=$market")
end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-users-saved-episodes

"""
    episodes_get_saved(;limit::Int64=20, market::String="US", offset::Int64=0)

**Summary**: Get a list of the episodes saved in the current Spotify user's library.
             This API endpoint is in **beta** and could change without warning.

# Optional keywords
- `limit::Int64` : Maximum number of items to return, default is set to 20
- `market::String` : An ISO 3166-1 alpha-2 country code. If a country code is specified,
                     only episodes that are available in that market will be returned.
                     Default is set to "US".
- `offset::Int64` : Index of the first item to return, default is set to 0

# Example
```julia-repl
julia> Spotify.episodes_get_saved()[1]
[ Info: We try the request without checking if current grant includes scope user-library-read.
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 7 entries:
  :href     => "https://api.spotify.com/v1/me/episodes?offset=0&limit=20&market=US"
  :items    => Union{}[]
  :limit    => 20
"""
function episodes_get_saved(;limit::Int64=20, market::String="US", offset::Int64=0)

    return spotify_request("me/episodes?limit=$limit&market=$market&offset=$offset";
    scope = "user-library-read")

end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/check-users-saved-episodes

"""
    episodes_get_contains(ids)

**Summary**: Check if one or more episodes is already saved in the current Spotify user's 'Your Episodes' library.
             This API endpoint is in **beta** and could change without warning.

# Arguments
- `ids` : A comma-separated list of the Spotify episode IDs

# Example
```julia-repl
julia> Spotify.episodes_get_contains("77o6BIVlYM3msb4MMIL1jH,0Q86acNRm6V9GYx55SXKwf")[1]
[ Info: We try the request without checking if current grant includes scope user-library-read.
2-element JSON3.Array{Bool, Base.CodeUnits{UInt8, String}, Vector{UInt64}}:
 0
 0
```
"""
function episodes_get_contains(ids)
    return spotify_request("me/episodes/contains?ids=$ids"; scope = "user-library-read")
end


text:
abstract type FluidSolution end

mutable struct DiffusionSolution{T<:Union{Array{Float64,2},Nothing}} <: FluidSolution
    u::Array{Float64,2}

    u_analytic::T
    errors::T

    prob::DiffusionProblem
    retcode::Bool
end
function DiffusionSolution(prob::DiffusionProblem, retcode::Bool)
    u_analytic, errors = get_analytic(prob.model, prob.cache.u, prob.cache.t)
    return DiffusionSolution(prob.cache.u, u_analytic, errors, prob, retcode)
end

function get_analytic(model::DiffusionModel{T}, u, t) where T<:Function
    u_analytic = model.u_analytic(t)
    errors = abs.(u - u_analytic)
    return u_analytic, errors
end
get_analytic(model::DiffusionModel{Nothing}, u, t) = nothing, nothing
function solve(prob::DiffusionProblem)
    sol = solve!(prob)
    return sol
end

function solve!(prob::DiffusionProblem)
    while prob.cache.t < prob.model.parameters.tspan[2]
        timestep!(prob)
    end
    return DiffusionSolution(prob, true)
end
function timestep!(prob::DiffusionProblem)
    prob.cache.t += prob.model.Δt
    perform_step!(prob.model, prob.method, prob.cache)
    apply_BCs!(prob.cache)
    return nothing
end