text:
using Dates
struct LSB
    lsb # используется для расчетов
    amp # только для справки, 0 = не задано
    div # только для справки, 0 = не задано
end

#если просто  Int, то потом выдает ошибку в parseChAttr
LSB(amp, div) = LSB(amp / div, amp,div)
LSB(lsb::Float64) = LSB(lsb, 0, 0)

tounits(lsb::LSB, pt) = pt * lsb.lsb # из точки в физ. величину
topoints(lsb::LSB, x) = x / lsb.lsb # из физ. величины в точки (дробные)
tointpoints(lsb::LSB, x) = floor(Int, x / lsb.lsb) # из физ. величины в точки (целые)

struct PhysicalInfo # или иметь только один вариант и не париться?
    lsb::LSB
    units::String # названия ед. измерения
    maxval::Float64 # ?? нужны ли ??
    minval::Float64 # ??
end

struct FeatureInfo
    mode::String
    featureset #:: ? names, masks, checks
end

abstract type DataSet{T} end

mutable struct InMemory_data{T} #мб RingBuffer???
    ind #индексы, откуда прочитано
    data:: Vector{T}
    maxlen::Int64
    function InMemory_data{T}(maxlen::Int) where T
        new{T}(0:0,Vector{T}(undef,0), maxlen) #, 0)
    end
end

mutable struct PhysicalDataSet{T} <: DataSet{T} # а и б, только режим разный
    type::String # “signal” / “param”
    lsb::PhysicalInfo
    name::String
    datatype

    data
    data_IM::InMemory_data{T}
end

mutable struct IndexDataSet{T} <: DataSet{T}
    TG::TimeGrid
    name::String
    datatype

    data
    data_IM::InMemory_data{T}
end

mutable struct IntervalDataSet{T} <: DataSet{T}
    offsetdata #::IndexDataSet
    TG::TimeGrid
    name::String
    datatype

    data
    data_IM::InMemory_data{T}
end

mutable struct FeatureDataSet{T} <: DataSet{T}
    mask::Dict
    name::String
    datatype::DataType

    data
    data_IM:: InMemory_data{T}
end

#поиск информации о данных по их имени в дереве
function findnode(tree::Dict, dataname::String) where T
    Nnodes = length(tree["nodes"])
    for k=1:Nnodes
        node = tree["nodes"][k]
        if dataname==node["name"]
            return dataTree = node
        else
            dataTree = findnode(tree["nodes"][k]["nodes"], dataname)
            if ~isempty(dataTree)
                return dataTree
            end
        end
    end
    return []
end
function findnode(tree::Array, dataname::String) where T
    Nnodes = length(tree)
    for k=1:Nnodes
        node = tree[k]
        if dataname==node["name"]
            return dataTree = node
        else
            dataTree = findnode(tree[k]["nodes"],dataname)
            if ~isempty(dataTree)
                return dataTree
            end
        end
    end
    return []
end
#если атрибуты уже были считаны
function parseAttr(attr::Dict)
    #преобразованное время начала
    if haskey(attr, "TimeStart")
        if isa(attr["TimeStart"],Union{Int64,Float64})
            TimeStart = valToDate(attr["TimeStart"])
        else
            TimeStart = attr["TimeStart"]
        end
    else
        TimeStart = valToDate(0)
    end
    if haskey(attr, "Freq")
        Fs = Float64(attr["Freq"])
    else
        Fs = 1
    end
    #собираем временную сетку
    TG = TimeGrid(TimeStart, Fs)
    if haskey(attr, "amp")
        lsb = LSB(attr["amp"],attr["div"])
    else
        lsb = LSB(1,1)
    end
    if haskey(attr, "unit")
        units = attr["unit"]
    else
        units = ""
    end

    #minmax по умолчанию пока зануляем
    PhI = PhysicalInfo(lsb,units,0,0)

    return TG,PhI
end

function valToDate(datenum::Union{Int64,Float64})
    #преобразуем время из дней в DateTime
    dateMs = trunc(Int,datenum*24*60*60*1000) #врем в миллисекундах
    datenum = dateMs + Dates.value(DateTime(1900,12,30,00,00,00))
    dateTime = Dates.epochms2datetime(datenum)
end

#конструктор лямбда - функций для чтения
function construct_dsreadInRange(baseIP::IPv4,port::Union{String,Int64},groupname::String,dsname::String)
    fc = (from,to,count) -> getStructData(baseIP,port,groupname,dsname,"0",from,to,count)
    return fc
end
function construct_dsread(baseIP::IPv4,port::Union{String,Int64},dsname::String)
    fc = (from,to,count) -> getData(baseIP,port,dsname,"0",from,to,count)
    return fc
end

#создаем новый датасет - точнее открывем его и читаем атрибуты поля
function ds_new(baseIP::IPv4,port::Union{String,Int64},groupname::String,dsname::String, attr::Dict,TG = [], PhInfo = [])
    #если ничего не известно про атрибуты
    #информация о канале
    if isempty(TG)
        TG, PhInfo = parseAttr(attr)
    end
    # @info d
    dstype = attr["dstype"] #signal/param/index/interval/feature
    if dstype == "feature"
        mask = getMask(parseType(attr))
        if isempty(mask)
            #если не было маски, то называем тип по имени данных
            mask[dsname] = AllMask(UInt32(1),UInt32(1))
        end
    else
        mask = Dict()
    end

    dstype = attr["dstype"]
    #просто чтение по номеру в массиве
    ds_reader= construct_dsread(baseIP,port,dsname)

    if dstype == "signal"
        ds_readerInRange =  construct_dsread(baseIP,port,dsname)
    else
        ds_readerInRange =  construct_dsreadInRange(baseIP,port,groupname,dsname)
    end

    DS_obj= ds_new(baseIP,port,dsname,attr,ds_reader,ds_readerInRange,TG,PhInfo,mask)
    return DS_obj
end

#создаем объект набора данных
function ds_new(baseIP::IPv4,port::Union{String,Int64},dsname::String,attr::Dict,ds_reader,ds_readerInRange,TG::TimeGrid,PhInfo::PhysicalInfo,mask::Dict)
    dstype = attr["dstype"] #signal/param/index/interval/feature
    datatype = attr["datatype"]
    # datatype = getDataType(baseIP,port,dsname)
    # @info dstype
    if dstype == "interval"
        if haskey(attr,"offsetdata")
            OD = attr["offsetdata"]
        else
            OD = nothing
        end
        DS_obj = IntervalDataSet{datatype}(OD, TG, dsname, datatype, (ds_reader, ds_readerInRange),InMemory_data{datatype}(100))
    elseif dstype == "feature"
        DS_obj = FeatureDataSet{datatype}(mask,dsname, datatype, (ds_reader, ds_readerInRange),InMemory_data{datatype}(100))
    elseif dstype == "index"
        DS_obj = IndexDataSet{datatype}(TG, dsname, datatype, (ds_reader, ds_readerInRange),InMemory_data{datatype}(100))
    elseif dstype == "signal" || dstype == "param"
        DS_obj = PhysicalDataSet{datatype}(dstype,PhInfo,dsname, datatype, (ds_reader, ds_readerInRange),InMemory_data{datatype}(100))
    end
    DS_obj
end


text:
"""Version of UCRL2 adapted to finite horizon
Confidence intervals taken from Ian Osband's TabulaRL
package"""
type UCRL2<:FiniteHorizonAgent
    policy::Matrix{UInt64}
    ns::Array{UInt64, 4}
    phat::Array{Float64, 4}
    n::Array{UInt64, 3}
    nstemp::Array{UInt64, 4}
    ntemp::Array{UInt64, 3}
    rewards::Array{Float64, 3}
    rewardstemp::Array{Float64, 3}
    Vopt::Matrix{Float64}
    δ::Float64
    rewards_known::Bool
    maxRet::Float64
    maxR::Float64
    immediate_updates::Bool
    update_necessary::Bool
end

function UCRL2(S, A, H, rewards, δ, immediate_updates=false)
    policy = zeros(UInt64, S, H)
    ns = zeros(UInt64, S, S, A, H)
    phat = zeros(S, S, A, H)
    phat[1, :, :, :] = 1.
    n = zeros(UInt64, S, A, H)
    nstemp = zeros(UInt64, S, S, A, H)
    ntemp = ones(UInt64, S, A, H)
    rewardstemp = zeros(S, A, H)
    maxRet = sum(maximum(reshape(rewards, S*A, H), 1))
    m = UCRL2(policy, ns, phat, n, nstemp, ntemp, rewards, rewardstemp, zeros(S, H+1), δ, true, maxRet, maximum(rewards[:]), immediate_updates, false)
    rand!(m.policy, 1:A) # Initialize with uniformly random policy
    m
end

function observe!(m::UCRL2, s, a, t, r, sn)
    m.nstemp[sn, s, a, t] += 1
    m.ntemp[s, a, t] +=  1
    m.rewardstemp[s, a, t] += r
    if m.immediate_updates || m.ntemp[s,a,t] >= max(1, m.n[s,a,t])
        m.update_necessary = true
    end
end
maxV(m::UCRL2) = m.maxRet
maxR(m::UCRL2) = m.maxR

function update_policy!(m::UCRL2, force=false)
    if !force && !m.update_necessary
        return
    end

    # First make all counts up-to-date
    m.n += m.ntemp
    m.ns += m.nstemp
    if !m.rewards_known
        m.rewards += m.rewardstemp
    end
    fill!(m.ntemp, 0)
    fill!(m.nstemp, 0)
    fill!(m.rewardstemp, 0.)
    m.update_necessary=false

    S = nS(m)
    H = horizon(m)
    A = nA(m)
    δ = m.δ
    Vmax = maxV(m)
    Rmax = maxR(m)
    Q = zeros(A)
    k = sum(m.n[:,:, 1]) # number of observed episodes
    # precomputations for confidence intervals
    lR = 4 * log(4 * S * H * A * (k+1) / δ)
    lC = 4 * S * log(2 * S * H * A * (k+1) / δ)

    m.phat .= m.ns ./ reshape(m.n, 1, S, A, H)
    # Extended value iteration
    for t=H:-1:1
        V = m.Vopt[:, t+1]
        ind = sortperm(V)
        curmaxV = min(V[ind[end]], Vmax)

        for s ∈ 1:S
            for a ∈ 1:A
                slack = √(lC / m.n[s, a, t])
                mass_left = 1.
                EV = 0.
                for i ∈ S:-1:1
                    ii = ind[i]
                    ph = m.phat[ii, s, a, t]
                    pmax = min(mass_left, ph + slack)
                    slack -= pmax - ph
                    mass_left -= pmax
                    EV += pmax * V[ii]
                    if mass_left < 4 * eps(typeof(mass_left))
                        break
                    end
                end

                r = m.rewards[s, a, t]
                if !m.rewards_known
                    min(m.rewards[s, a, t] / n[s, a, t] + Rmax * √(lR / m.n[s, a, t]), Rmax)
                end

                Q[a] = r + EV
            end
            bestA = indmax(Q)
            m.policy[s, t] = bestA
            m.Vopt[s, t] = Q[bestA]
        end
    end
end

ext: jl


text:
using PyPlot
using FileIO

mystep(r::AbstractRange) = step(r)
mystep(r::AbstractArray) = r[2] - r[1]

extendrange(r) = range( (first(r)-mystep(r)/2), last(r)+mystep(r)/2, length=length(r)+1)

function plot_heatmap_tracing_prob_vs_c(
  results,
  tracing_probs = 0:0.05:1,
  Cs=0:0.05:1;
  cmin=nothing,
  cmax=nothing,
  logscale=true,
  addcbar::Bool=true,
  polish::Bool=false)

  reduction = 1 .- Cs / 1.35  |> collect

  if nothing == cmax
    cmax = maximum(results)
  end

  if nothing == cmin
    cmin = minimum(results)
  end

  if logscale
    img = pcolor(
      extendrange(reduction),
      extendrange(tracing_probs),
      results,
      norm=matplotlib.colors.LogNorm(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
      clim(vmin=cmin)
  else
    img = pcolor(
      extendrange(reduction),
      extendrange(tracing_probs),
      results,
      norm=matplotlib.colors.Normalize(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
  end

  if addcbar
    cbar = colorbar()
  end

  if polish
    xlabel("f - Stopień redukcji kontaktów")
    ylabel("b - skuteczność śledzenia kontaktów")
  else
    xlabel("f - contact reduction rate")
    ylabel("b - contact tracing efficiency")
  end

  xticks(
    0:0.1:1,
    ["$(100*f)%" for f in 0:0.1:1],
    rotation=60
  )

  gca().invert_yaxis()
  gca().invert_xaxis()
  img
end

function plot_heatmap_mild_detection_vs_tracing_prob(
    results,
    mild_detections::AbstractVector{T} where T <: Real = 0:0.05:1,
    tracking_probs::AbstractVector{T} where T <: Real = 0:0.05:1;
    cmin=nothing,
    cmax=nothing,
    logscale=true,
    addcbar=true,
    polish::Bool=false
  )

  if nothing == cmax
    cmax = maximum(results)
  end

  if nothing == cmin
    cmin = minimum(results)
  end

  if logscale
    im = pcolor(
      extendrange(mild_detections),
      extendrange(tracking_probs),
      results,
      norm=matplotlib.colors.LogNorm(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
      clim(vmin=cmin)
  else
    im = pcolor(
      extendrange(mild_detections),
      extendrange(tracking_probs),
      results,
      norm=matplotlib.colors.Normalize(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
  end

  clim(vmin=cmin)
  if addcbar
    colorbar()
  end

  gca().invert_yaxis()
  gca().invert_xaxis()
  if polish
    xlabel("q' - prawdopodobieństwo wykrycia lekkich przypadków")
    ylabel("b - skuteczność śledzenia kontaktów")
  else
    xlabel("q' - mild case detection probability")
    ylabel("b - contact tracing efficiency")
  end
  im
end

function plot_heatmap_mild_detection_vs_c(
    results,
    mild_detection_probs::AbstractVector{T} where T<:Real = 0:0.05:1,
    Cs::AbstractVector{T} where T <:Real=0:0.05:1;
    cmin=nothing, cmax=nothing,
    logscale::Bool=true,
    addcbar::Bool=true,
    polish::Bool=true)
  reduction = 1 .- Cs / 1.35  |> collect

  if nothing == cmax
    cmax = maximum(results)
  end

  if nothing == cmin
    cmin = minimum(results)
  end

  if logscale
    im = pcolor(
      extendrange(reduction),
      extendrange(mild_detection_probs),
      results,
      norm=matplotlib.colors.LogNorm(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
      clim(vmin=cmin)
  else
    im = pcolor(
      extendrange(reduction),
      extendrange(mild_detection_probs),
      results,
      norm=matplotlib.colors.Normalize(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
  end
  if addcbar
    c = colorbar()
  end

  gca().invert_yaxis()
  gca().invert_xaxis()

  if polish
    xlabel("f - Stopień redukcji kontaktów")
    ylabel("q' - Skuteczność wykrywania lekkich przypadków")
  else
    xlabel("f - contact reduction rate")
    ylabel("q' - mild case detection probability")
  end

  xticks(
    0:0.1:1,
    ["$(Int(100*f))%" for f in 0:0.1:1],
    rotation=60)

  im
end

function plot_heatmap_c_vs_phone_tracing_usage(results, Cs=0:0.05:1, phone_tracking_usage = 0:0.05:1; tlim=0.0001, cmin=nothing, cmax=nothing, logscale=true, addcbar::Bool=true, polish::Bool=false)
  reduction = 1 .- Cs / 1.35

  if nothing == cmax
    cmax = maximum(results)
  end

  if nothing == cmin
    cmin = minimum(results)
  end

  if logscale
    im = pcolor(
      extendrange(reduction),
      extendrange(phone_tracking_usage),
      results,
      #norm=matplotlib.colors.SymLogNorm(linthresh=tlim, vmin=cmin, vmax=cmax),
      norm=matplotlib.colors.LogNorm(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
      clim(vmin=cmin)
  else
    im = pcolor(
      extendrange(reduction),
      extendrange(phone_tracking_usage),
      results,
      norm=matplotlib.colors.Normalize(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
  end

  if addcbar
    cbar = colorbar()
  end

  if polish
    xlabel("f - Stopień redukcji kontaktów")
    ylabel("u - część populacji używająca aplikacji \n do śledzenia kontaktów")
  else
    xlabel("f - contact reduction rate")
    ylabel("u - fraction of population \n using the contact tracing app")
  end

  xticks(
    0:0.1:1,
    ["$(Int(100*f))%" for f in 0:0.1:1],
    rotation=60)

  gca().invert_yaxis()
  gca().invert_xaxis()

  return im
end

function plot_heatmap_phone_tracking_usage_vs_tracking_prob(
  results,
  phone_tracking_usages = 0:0.05:1,
  tracking_probs=0:0.05:1;
  cmin=nothing,
  cmax=nothing,
  logscale=true,
  addcbar=true
)

  if nothing == cmax
    cmax = maximum(results)
  end

  if nothing == cmin
    cmin = minimum(results)
  end

  if logscale
    im = pcolor(
      extendrange(phone_tracking_usages),
      extendrange(tracking_probs),
      results,
      norm=matplotlib.colors.LogNorm(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
      clim(vmin=cmin)
  else
    im = pcolor(
      extendrange(phone_tracking_usages),
      extendrange(tracking_probs),
      results,
      norm=matplotlib.colors.Normalize(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
  end

  clim(vmin=cmin)

  if addcbar
    colorbar()
  end

  gca().invert_xaxis()
  gca().invert_yaxis()
  xlabel("u - część populacji używająca aplikacji śledzącej")
  ylabel("b - skuteczność śledzenia kontaktów")
  return im
end

function plot_initiators_vs_c(results, initiators = 0:2:50, Cs=0:0.05:1; cmin=nothing, cmax=nothing, logscale=true, addcbar::Bool=true, polish::Bool=false)
  #figure(figsize=(10,5))
  reduction = 1 .- Cs / 1.35  |> collect

  if nothing == cmax
    cmax = maximum(results)
  end

  if nothing == cmin
    cmin = minimum(results)
  end

  if logscale
    img = pcolor(
      extendrange(reduction),
      extendrange(initiators),
      results,
      norm=matplotlib.colors.LogNorm(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
      clim(vmin=cmin)
  else
    img = pcolor(
      extendrange(reduction),
      extendrange(initiators),
      results,
      norm=matplotlib.colors.Normalize(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
  end

  if addcbar
    cbar = colorbar()
  end

  if polish
    xlabel("f - stopień redukcji kontaktów")
    ylabel("n - liczba początkowych zarażeń")
  else
    xlabel("f - contact reduction rate")
    ylabel("n - number of initiators")
  end

  xticks(
    0:0.1:1,
    ["$(100*f)%" for f in 0:0.1:1],
    rotation=60
  )

  #gca().invert_yaxis()
  gca().invert_xaxis()
  img
end


text:
# Utilities to deal with fetching/serving actual Pkg resources

const uuid_re = raw"[0-9a-z]{8}-[0-9a-z]{4}-[0-9a-z]{4}-[0-9a-z]{4}-[0-9a-z]{12}(?-i)"
const hash_re = raw"[0-9a-f]{40}"
const meta_re     = Regex("^/meta\$")
const registry_re = Regex("^/registry/($uuid_re)/($hash_re)\$")
const resource_re = Regex("""
    ^/registries\$
  | ^/registry/$uuid_re/$hash_re\$
  | ^/package/$uuid_re/$hash_re\$
  | ^/artifact/$hash_re\$
""", "x")
const hash_part_re = Regex("/($hash_re)\$")

function get_registries(config, server::PkgStorageServer)
    regs = Dict{String, String}()
    response = HTTP.get("$(server.url)/registries")
    for line in eachline(IOBuffer(response.body))
        m = match(registry_re, line)
        if m !== nothing
            uuid, hash = m.captures
            regs[uuid] = hash
        else
            @error "invalid response" server=server.url resource="/registries" line=line Dates.now()
        end
    end
    return regs
end

"""
    write_atomic(f::Function, path::String)

Performs an atomic filesystem write by writing out to a file on the same
filesystem as the given `path`, then `move()`'ing the file to its eventual
destination.  Requires write access to the file and the containing folder.
Currently stages changes at "<path>.tmp.<randstring>".  If the return value
of `f()` is `false` or an exception is raised, the write will be aborted.
"""
function write_atomic(f::Function, path::String)
    temp_file = path * ".tmp." * randstring()
    try
        retval = open(temp_file, "w") do io
            f(temp_file, io)
        end
        if retval !== false
            mv(temp_file, path; force=true)
        end
        return retval
    catch e
        rm(temp_file; force=true)
        rethrow(e)
    end
end

# Current registry hashes.
const REGISTRY_HASHES = Dict{String,String}()

const last_registry_update = Ref{Float64}(0)

function update_registries(config)
    t = time()
    if t < last_registry_update[] + config.min_time_between_registry_updates
        return
    end

    changed = false
    # Collect current registry hashes from servers.
    for server in config.storage_servers
        for (uuid, hash) in get_registries(config, server)
            if get(REGISTRY_HASHES, uuid, "") != hash
                REGISTRY_HASHES[uuid] = hash
                changed = true
            end
        end
    end

    # Write new registry info to file.
    if changed
        write_atomic(joinpath(config.cache_dir, "registries")) do temp_file, io
            for uuid in sort!(collect(keys(REGISTRY_HASHES)))
                hash = REGISTRY_HASHES[uuid]
                println(io, "/registry/$uuid/$hash")
            end
            return true
        end
    end
    last_registry_update[] = t
    return
end

function fetch(config::Config, resource::AbstractString)
    resource == "/registries" && update_registries(config)
    path = config.cache_dir * resource
    isfile(path) && return path
    update_registries(config)
    servers = config.storage_servers
    isempty(servers) && throw(@error "fetch called with no servers" resource=resource)

    mkpath(dirname(path))
    for server in servers
        if download(config, server, resource, path)
            break
        end
    end
    success = isfile(path)
    success || @warn "download failed" resource=resource  Dates.now()
    return success ? path : nothing
end

function tarball_git_hash(tarball::String)
    local tree_hash
    mktempdir() do tmp_dir
        open(tarball) do io
            Tar.extract(GzipDecompressorStream(io), tmp_dir)
        end
        tree_hash = bytes2hex(Pkg.GitTools.tree_hash(tmp_dir))
        chmod(tmp_dir, 0o777, recursive = true)
    end
    return tree_hash
end

function download(config, server::StorageServer, resource::AbstractString,
                  path::AbstractString)
    @info "downloading resource" server=server resource=resource Dates.now()
    hash = let m = match(hash_part_re, resource)
        m !== nothing ? m.captures[1] : nothing
    end

    write_atomic(path) do temp_file, io
        if !get_resource_from_storage_server!(config, server, resource, io)
            return false
        end

        # If we're given a hash, then check tarball git hash
        if hash !== nothing
            tree_hash = tarball_git_hash(temp_file)
            # Raise warnings about resource hash mismatches
            if hash != tree_hash
                @warn "resource hash mismatch" server=server resource=resource hash=tree_hash Dates.now()
                return false
            end
        end

        return true
    end
end

function get_resource_from_storage_server!(config, server::PkgStorageServer,
                                           resource, io)
    response = HTTP.get(status_exception = false,
                        response_stream = io,
                        server.url * resource)

    # Raise warnings about bad HTTP response codes
    if response.status != 200
        @warn "response status $(response.status)" Dates.now()
        return false
    end

    return true
end

function serve_file(http::HTTP.Stream, path::String)
    HTTP.setheader(http, "Content-Length" => string(filesize(path)))
    # We assume that everything we send is gzip-compressed (since they're all tarballs)
    HTTP.setheader(http, "Content-Encoding" => "gzip")
    startwrite(http)

    # Open the path, write it out directly to the HTTP stream
    open(io -> write(http, read(io, String)), path)
end


text:
export examples, example

const Example = NamedTuple{(:name, :inv, :vars, :kwargs),Tuple{Symbol, Invariant, Vector{Symbol}, Any}}

_examples = Dict{Symbol,Example}()

example!(x::Example) = push!(_examples, x.name=>x)
examples() = values(_examples) |> collect
example(s) = _examples[s]

# Double
# x, y = 0, 1
# while true
#     x = 2x
#     y = 1/2 y + 1
# end
example!(Example((
    :double1,
    @invariant(x*y == 2x),
    [:x, :y],
    Dict()
)))
example!(Example((
    :double2,
    @invariant(x == 2y),
    [:x, :y],
    Dict()
)))

# Square
# a, b = 0, 0
# while true
#     a = a + 2b + 1
#     b = b + 1
# end
example!(Example((
    :square,
    @invariant(a == b^2),
    [:a, :b],
    Dict()
)))

example!(Example((
    :dblsquare,
    @invariant(a == 2*b^2),
    [:a, :b],
    Dict()
)))

# Sum1
# a, b, c = 0, 0, 1
# while true
#     a = a + 1
#     b = b + c
#     c = c + 2
# end
example!(Example((
    :sum1,
    @invariant(1 + 2a == c && 4b == (c-1)^2),
    [:a, :b, :c],
    Dict()
)))

# Sum2
# a, b, c, s = 0, 0, 1, 0
# while true
#     a = a + 1
#     b = b + c
#     c = c + 2
#     s = s + 2a + 1
# end
example!(Example((
    :sum2,
    @invariant(1 + 2a == c && b + c == 1 + s && c*(c+2) == 3 + 4s),
    [:s, :a, :b, :c],
    Dict()
)))

# eucliddiv
# r, q = x, 0
# while true
#     r = r - y
#     q = q + 1
# end
example!(Example((
    :eucliddiv,
    @invariant(x(0) == y(0)*q(n) + r(n)),
    [:r, :q, :x, :y],
    Dict(:params=>[:x, :y])
)))

# Integer Square Root - version 1
# k, j, m = 0, 1, 1
# while m<=n
#     k = k + 1
#     j = j + 2
#     m = m + j
# end
example!(Example((
    :intsqrt1,
    @invariant(j == 1 + 2k && (j+1)^2 == 4m),
    [:m, :k, :j],
    Dict()
)))

# Integer Square Root - version 2
# y, r = 1/2*a, 0
# while true
#     y = y - r
#     r = r + 1
# end
example!(Example((
    :intsqrt2,
    @invariant(a(0) + r == r^2 + 2y),
    [:y, :r, :a],
    Dict(:params=>[:a])
)))

# Integer Cubic Root
# x, r, s = a, 1, 13/4
# while true
#     x = x-s
#     s = s + 6r + 3
#     r = r + 1
# end
example!(Example((
    :intcbrt,
    @invariant(1/4 + 3r^2 == s && 1 + 4*a(0) + 6r^2 == 3r + 4r^3 + 4x),
    [:x, :s, :r, :a],
    Dict(:params=>[:a])
)))

# Consecutive Cubes
# n, x, y, z = 0, 0, 1, 6
# while true
#     x = x + y
#     y = y + z
#     z = z + 6
#     n = n + 1
# end
example!(Example((
    :cubes,
    @invariant(n^3 == x && 1 + 3n + 3n^2 == y && 6 + 6n == z),
    [:x, :y, :z, :n],
    Dict()
)))

# Petter 1
# x, y = 0, 0
# while true
#     x = x + y^1
#     y = y + 1
# end
example!(Example((
    :petter1,
    @invariant(y^2 == 2*x+y),
    [:x, :y],
    Dict()
)))

# Petter 2
# x, y = 0, 0
# while true
#     x = x + y^2
#     y = y + 1
# end
# 6 x == y * (-1 + 3 y - 2 y^2)

# Petter 3
# x, y = 0, 0
# while true
#     x = x + y^3
#     y = y + 1
# end
# 4 x - (-1 + y)^2 y^2

# Add (https://rise4fun.com/Dafny/Add)
# r = x
# n = y
# while n != 0
#     r = r + 1
#     n = n - 1
# end
example!(Example((
    :add1,
    @invariant(r(m) == x(0)+y(0)-n(m)),
    [:r, :n, :x, :y],
    Dict(:params=>[:x, :y])
)))

# r = 2x
# n = y
# while n != 0
#     r = r + 1
#     n = n - 1
# end
example!(Example((
    :add2,
    @invariant(r(m) == 2*x(0)+y(0)-n(m)),
    [:r, :n, :x, :y],
    Dict(:params=>[:x, :y])
)))

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

function Base.run(examples::Vector{Example}; solvers=[Yices,Z3], timeout=2, maxsol=1, shapes=[UnitUpperTriangular, UpperTriangular, FullSymbolic])
    progress = ProgressUnknown("Instances completed:")
    df = DataFrame(Solver = Type{<:NLSolver}[], Instance = Any[], Shape = MatrixShape[], Status = NLStatus[], Elapsed = Any[], Result = Union{RecSystem,Nothing}[])
    for ex in examples
        for shape in shapes
            strat = strategy_fixed(ex.inv, copy(ex.vars), shape, [length(ex.vars)+1]; constone=true, ex.kwargs...)
            for solver in solvers
                sols = solutions(strat; solver=solver, timeout=timeout, maxsol=maxsol)
                for res in sols
                    push!(df, (solver, ex.name, shape, res.status, res.elapsed, res.recsystem))
                    # push!(df, (solver, ex.name, shape, res.status, res.status == NLSat.timeout ? "-" : res.elapsed, res.recsystem))
                    next!(progress)
                end
            end
        end
    end
    finish!(progress)
    df
end

function wide(df)
    key = [:Instance, :Shape]
    tab = unstack(df, key, :Solver, :Elapsed)
    z3 = unstack(tab, :Instance, :Shape, :Z3Solver)
    yices = unstack(tab, :Instance, :Shape, :YicesSolver)
    rename!(z3,
        :UnitUpperTriangular => :Z3_UnitUpperTriangular,
        :UpperTriangular     => :Z3_UpperTriangular,
        :FullSymbolic        => :Z3_FullSymbolic
    )
    rename!(yices,
        :UnitUpperTriangular => :Yices_UnitUpperTriangular,
        :UpperTriangular     => :Yices_UpperTriangular,
        :FullSymbolic        => :Yices_FullSymbolic
    )
    join(yices, z3, on = :Instance, makeunique = true)
end


text:
using OptimNewton
using Test

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


text:
import Base.tanh
sigm(x::Float64, γ::Float64=1.) = 1.0/(1.0 + exp(-x*γ))
tanh(x::Real, γ::Real) = (exp(2x*γ)-1) / (exp(2x*γ)+1)
relu(x::Real) = max(zero(eltype(x)),x)

type Synapse
    source
    destination
    weight::Float64
    function Synapse(source, destination, weight::Float64)
        s = new(source, destination, weight);
        addSynapse(destination, s); s
    end
end

incoming(s::Synapse) = s.weight * s.source.output

function Base.show(io::IO, s::Synapse)
    @printf(io,"%d -> %3.5f -> %d\n", s.source.id, s.weight, s.destination.id)
end

type Neuron
    id::Int64
    synapses::Vector{Synapse}
    bias::Float64
    nType::Symbol       # [:INPUT, :OUTPUT, :HIDDEN]
    activation::Function  # [:sigm, :tanh, :relu]
    response::Float64   # default = 4.924273 (Stanley, p. 146)
    output::Float64     # for recurrent networks all neurons must have an "initial state"
    function Neuron(neurontype::Symbol, id::Int64, bias::Float64=0., activation::Symbol=:sigm, γ::Float64=1.)
        f = activation == :sigm?  x->sigm(x,γ):activation ==:tanh? x->tanh(x,γ):activation ==:relu? x->relu(x): x->x
        new(id,[], bias, neurontype, f, γ, 0.)
    end
end

addSynapse(n::Neuron, s::Synapse) = push!(n.synapses,s)   # adds the synapse to the destination neuron

function activate(n::Neuron)
#         "Activates the neuron"
        @assert n.nType != :INPUT
        return n.activation(updateActivation!(n) + n.bias)
end

function updateActivation!(n::Neuron)
    output = 0.
    for s in n.synapses
        output += incoming(s)
    end
    return output
end

# "Prints neuron's current state"
function currentOutput(n::Neuron)
    println("state: $(n.id) - output: $(n.output)")
    return n.output
end

function Base.show(io::IO, n::Neuron)
    @printf(io,"%d %s\n", n.id, n.nType)
end

type Network
    neurons::Vector{Neuron}
    synapses::Vector{Synapse}
    numInputs::Int64
    nntype::ChromoType
    Network(numInputs::Int64) = new([],[], numInputs, Recurrent())
    Network(neurons::Vector{Neuron}, synapses::Vector{Synapse}, numInputs::Int64, nntype::ChromoType) =
        new(neurons, synapses, numInputs, nntype)
end

addNeuron(network::Network, neuron::Neuron) = push!(network.neurons,neuron)

addSynapse(network::Network, synapse::Synapse) = push!(network.synapses,synapse)

flush!(network::Network) = for n in network.neurons  n.output = 0. end

function Base.show(io::IO, net::Network)
    @printf(io,"%s %d nodes and %d synapses\n",net.nntype, length(net.neurons), length(net.synapses))
end

function activate(::FeedForward, nn::Network, inputs::Vector)

    #=  Serial (asynchronous) network activation method. Mostly
    used  in classification tasks (supervised learning) in
    feedforward topologies. All neurons are updated (activated)
    one at a time following their order of importance, so if
    you're defining your own feedforward topology, make sure
    you got them in the right order of activation. =#

    @assert length(inputs) == nn.numInputs

    for i  = 1:nn.numInputs
        nn.neurons[i].output = inputs[i] # iterates over inputs
    end

    # activate all neurons in the network (except for the inputs)
    netOutput = zeros(0)
    for n in nn.neurons[nn.numInputs+1:end]
        n.output = activate(n)
        if n.nType == :OUTPUT  push!(netOutput,n.output) end
    end

    return netOutput

end

function activate(::Recurrent, nn::Network, inputs::Vector)

    #= Parallel (synchronous) network activation method. Mostly used
    for control and unsupervised learning (i.e., artificial life)
    in recurrent networks. All neurons are updated (activated)
    simultaneously. =#

    @assert length(inputs) == nn.numInputs

    # the current state is like a "photograph" taken at each time step
    # reresenting all neuron's state at that time (think of it as a clock)
    inputCnt = 0
    currentState = zeros(0)
    for n in nn.neurons
        if n.nType == :INPUT
            inputCnt += 1
            n.output = inputs[inputCnt]
        else
            push!(currentState,activate(n))
        end
    end

    netOutput = zeros(0)
    for i = nn.numInputs+1:length(nn.neurons)
        n = nn.neurons[i]
        n.output = currentState[i-nn.numInputs]
        if n.nType == :OUTPUT
            push!(netOutput,n.output)
        end
    end

    return netOutput
end

function createPhenotype(ch::Chromosome)
    # Receives a chromosome and returns its phenotype (a neural network)

    neurons = Neuron[]
    Idx2Neuron = Dict{Int64,Neuron}()
    if ch.node_gene_type == Recurrent()

        for ng in ch.node_genes
            n = Neuron(ng.ntype, ng.id, ng.bias, ng.activation, ng.response)
            push!(neurons, n)
            Idx2Neuron[ng.id] = n
        end

    else # FeedForward() type

        # first create inputs
        for ng in ch.node_genes[1:ch.inputCnt]
            n = Neuron(ng.ntype, ng.id, ng.bias, ng.activation, ng.response)
            push!(neurons,n)
            Idx2Neuron[ng.id] = n
        end

        # Add hidden nodes in the right order
        for id in ch.node_order
            ng = ch.node_genes[id]
            n = Neuron(ng.ntype, ng.id, ng.bias, ng.activation, ng.response)
            push!(neurons, n)
            Idx2Neuron[ng.id] = n
        end

        # finally the output
        for ng in ch.node_genes[ch.inputCnt+1:ch.inputCnt+ch.outputCnt]
            n = Neuron(ng.ntype, ng.id, ng.bias, ng.activation, ng.response)
            push!(neurons, n)
            Idx2Neuron[ng.id] = n
        end

    end

    @assert length(neurons) == length(ch.node_genes)

    synapses = Synapse[]
    for (k,cg) in ch.connection_genes
        if cg.enable push!(synapses, Synapse(Idx2Neuron[cg.inId], Idx2Neuron[cg.outId], cg.weight)) end
    end

    return Network(neurons, synapses, ch.inputCnt, ch.node_gene_type)

end

###################################
#         previous version        #
###################################

# import Base.+, Base.tanh

# sigm(x::Real) = 1.0 / (1.0 + exp(-x))
# sigm(x::Real, γ::Real) = 1.0 / (1.0 + exp(-γ*x))
# tanh(x::Real, γ::Real) = (exp(2x*γ)-1) / (exp(2x*γ)+1)
# relu(x::Real) = max(zero(x),x)

# abstract Gene

# type GeneConnection <: Gene
#     in::Gene
#     out::Gene
#     f::Function
#     w::Float64
#     expressed::Bool
#     inId::Int64
#     outId::Int64
#     id::Int64
#     function GeneConnection(inNode, inID, outID, id)
#        (gene = new(inNode,inNode,x->x,randn(),true,inID,outID,id);gene.f=()->gene.in.f()*gene.w;gene)
#     end
# end

# type BiasNode <: Gene
#     f::Function # always returns 1.0
#     BiasNode() = (g = new(y->y); g.f=()->1.; g)
# end

# type InputNode <: Gene
#     f::Function # always returns value of x
#     x::Float64
#     InputNode(value::Float64) = (g = new(y->y,value); g.f=()->g.x; g)
#     InputNode() = InputNode(0.)
# end

# +(x1::Float64,x2::GeneConnection) = x1 + x2.f() # needed for reduce(+,...)
# +(x1::GeneConnection,x2::GeneConnection) = x1.f() + x2.f()

# type GeneNode <: Gene
#     xs::Vector{GeneConnection}
#     bias::Float64
#     f::Function
#     g::Function # usually either a sigmoid, tanh, or relu function
#     expressed::Bool
#     id::Int64
#     function GeneNode(id::Int64, xs::Vector{GeneConnection}, bias::Float64)
#        (gene = new(xs,bias,x->x,x->relu(x),true,id); # default to rectified linear unit
#         gene.f=()-> gene.g(reduce(+, 0., gene.xs) + gene.bias);
#         map(x->x.out=gene,xs);
#         gene)
#     end
#     GeneNode(id::Int64, bias::Float64) = GeneNode(id, Array(NEAT.GeneConnection,0), bias)
#     GeneNode(id::Int64) = GeneNode(id, Array(NEAT.GeneConnection,0), 0.)
# end

# type Genome
#     numIns::Int64
#     numOuts::Int64
#     outputs::Vector{Gene}
#     inputs::Vector{Gene}
#     nodes::Vector{GeneNode}
#     connections::Vector{GeneConnection}
# end

# # # Initialize Network
# function initGenome(params::NeatParams)

# end

# function addNode(g::Genome, connIdx::Int64)

#     # get reference to the out node
#     c = g.connections[connIdx]
#     out = c.out

#     # remove from out node's forward pass function
# #     idx = findfirst(o->o.id==connIdx, out)

# #     c.expressed = false
#     # add to output functions
#     return out

# end

max_stars_repo_name: Andy-P/NEAT

max_issues_repo_name: Andy-P/NEAT

max_forks_repo_name: Andy-P/NEAT


text:
################################################################################
# Copyright (C) Marcus Greiff 2020
#
# @file example_attitude_sys_id.jl
# @author Marcus Greiff
# @date July 2020
#
# This scrips contains some data and methods used in the system ideantification
# of Patrik's home-made drone. The basic methods and ideas are quite genetral,
# but not documeneted here.
#
################################################################################
using LinearAlgebra
using LaTeXStrings

# In the first experiment, we measure the thrust of the rotors as a function of
# a normalized command sent to the motors (typically pwm duty cycle), here,
# pairs of (measured force, duty cycle) are given in F and D respectively, from
# which the a and b parameters of the polynomial map are using a LS-regression.

# Generate some synthetic data
Theta    = rand(2,1);
N        = 10^3;
D        = rand(N,1);
Phi      = [D D.^2];
F        =  Phi*Theta + 0.1*randn(N,1);

#length = 55.0;

D = hcat(collect(50:50:600))
W = hcat([
    6,
    19,
    37,
    60,
    87,
    119,
    148,
    180,
    211,
    242,
    271,
    305,
])

g = 9.81;
m = 560;

D = D./1000;
W = W./1000;

#f = ma

F = W.*g * 2 / 4.0;

Phi      = [D D.^2];
Cov      = inv(Phi'*Phi);
Thetahat = (Phi'*Phi)\(Phi'*F);

using Plots
dd = hcat(collect((0:N)/N))
scatter(D, F, marker=:circle, color=:green, markeralpha = 0.3, xaxis="Duty cycle", yaxis="Force", label="Measurements")
#plot!(dd, [dd dd.^2]*Theta, color=:red, linewidth=2, label="True function")
plot!(dd, [dd dd.^2]*Thetahat, color=:blue, linewidth=2, label="Estimated function")

#1.1885640614704445
#2.291161604353094

# In the first experiment, we estimate C using the above function and
# measurements of the rotational torque about the e3 axis as a function of the
# applied PWM duty cycle. Here, we whould expect a linear relationship between
# the rotor thrust and the generated torque. If this is not linear, then c
# should represent the linear approximation of this function at the single rotor
# thrust of mg/4.

tauz = 4*0.005*9.81*0.09; # Nm
fi   = 2*1.18701; # N
# tauz = c * fi where fi is the thrust of one of the rotors
c = tauz / fi;

# index 0 - front right - pos force, positive moment about z, negative moment x, negative moment y
# index 1 - back right  - pos force, negative moment about z, negative moment x, positive moment y
# index 2 - back left   - pos force, positive moment about z, positive moment x, positive moment y
# index 3 - front left  - pos force, negative moment about z, positive moment x, negative moment y
#
#     3       0
#       .   .
#         .
#       .   .
#     2       1
#
# long 15.5
# short 9.3

A = [+1.0, +1.0, +1.0, +1.0;
     -dy, -dy, +dy, +dy;
     -dx, +dx, +dx, -dx;
     + c, - c, + c, - c]

# Cylindrical cross (mass, length of one axis, radius of the cross)
MC = 0.284
LC = 0.24
RC = 0.01

ICX = MC * (RC^2 / 4 + LC^2 / 12) + MC*RC^2 / 2;
ICY = ICX;
ICZ = MC * (RC^2 / 4 + LC^2 / 12) + MC * (RC^2 / 4 + LC^2 / 12);

# Battery box (Mass, length (x), with (y), height (z), offset in z)
ME = 0.172
WE = 0.035
LE = 0.09
HE = 0.025
DE = 0.03

IEX = ME * (WE^2/12 + HE^2/12 + DE^2);
IEY = ME * (LE^2/12 + HE^2/12 + DE^2);
IEZ = ME * (LE^2/12 + HE^2/12);

IXX = ICX + IEX;
IYY = ICY + IEY;
IZZ = ICZ + IEZ;


text:
# ------------------------------------------------------------------------------------------
# # 循环
#
# 话题：
# 1. `while`循环
# 2. `for`循环
#
#
# ## `while` 循环
#
# `while` 循环语法：
#
# ```julia
# while <condition>
#     <loop body>
# end
# ```
#
# 举个例子，可以使用 `while` 来数数或者遍历一个数组。
# ------------------------------------------------------------------------------------------

n = 0
while n < 10
    n += 1
    println(n)
end
n

myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

i = 1
while i <= length(myfriends)
    friend = myfriends[i]
    println("Hi $friend, it's great to see you!")
    i += 1
end

# ------------------------------------------------------------------------------------------
# ## `for` 循环
#
# `for` 循环语法：
# ```julia
# for <var> in <loop iterable>
#     <loop body>
# end
# ```
#
# 上面两个例子都可以用 `for` 循环改写：
# ------------------------------------------------------------------------------------------

for n in 1:10
    println(n)
end

myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

for friend in myfriends
    println("Hi $friend, it's great to see you!")
end

# ------------------------------------------------------------------------------------------
# 让我们来试试用 `for` 循环创建加和表，表的每个元素为其行索引和列索引的和。
#
# 首先，来初始化一个全零的数组（矩阵）。
# ------------------------------------------------------------------------------------------

m, n = 5, 5
A = fill(0, (m, n))

for i in 1:m
    for j in 1:n
        A[i, j] = i + j
    end
end
A

# ------------------------------------------------------------------------------------------
# 趁机介绍一下像上面这种嵌套 `for` 循环的语法糖
# ------------------------------------------------------------------------------------------

B = fill(0, (m, n))

for i in 1:m, j in 1:n
    B[i, j] = i + j
end
B

# ------------------------------------------------------------------------------------------
# 创建同样的表但更具有“Julia”风格的方法是使用**数组解析式（array comprehension）**。
# ------------------------------------------------------------------------------------------

C = [i + j for i in 1:m, j in 1:n]

# ------------------------------------------------------------------------------------------
# ## 练习
#
# ### 4.1
# 循环遍历 1 到 100 的整数并打印它们的平方。
# ------------------------------------------------------------------------------------------

# 在此作答

# ------------------------------------------------------------------------------------------
# ### 4.2
# 给上一题的代码添一点以创建一个字典，使字典存着每个整数和它的平方的键值对，就像这样
#
# ```julia
# squares[10] == 100
# ```
# ------------------------------------------------------------------------------------------

# 在此作答

@assert squares[10] == 100
@assert squares[11] == 121

# ------------------------------------------------------------------------------------------
# ### 4.3
# 用数组解析式来创建一个数组 `squares_arr`，让它存着从1到100的整数的平方。
# ------------------------------------------------------------------------------------------

squares_arr = [i for i in 1:100]

@assert length(squares_arr) == 100
@assert sum(squares_arr) == 338350

# ------------------------------------------------------------------------------------------
# 请在完成练习运行 `@assert` 块以检测结果是否正确，若无报错即为正确。
# ------------------------------------------------------------------------------------------

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
module MCMC

export runGibbs

using DataFrames
using CategoricalArrays
using StatsModels
using MixedModels
using Distributions,LinearAlgebra #samplers
using StatsBase
using Printf

include("equations.jl")
include("runTime.jl")
include("samplers.jl")

runGibbs = function(formula,userHints,userData,userPedData,nChain,nBurn,nThin,blockThese,VCV,VCVM;map,genotypes...)
	idY,yVec,FE,RE,ME,regionSizes = equations.mme(formula,userHints,userData,userPedData,blockThese;paths2geno=genotypes)
        samplers.runSampler(idY,yVec,FE,RE,nChain,nBurn,nThin,VCV,VCVM,ME,map,regionSizes)
	println("region sizes: $regionSizes")
        return(yVec,FE,RE,ME)
end

end

max_stars_repo_name: animBreedTools/NextGP

max_issues_repo_name: animBreedTools/NextGP

max_forks_repo_name: animBreedTools/NextGP


text:
module TestSpawnAt
using Test
import ThreadPools: StaticPool
using ThreadPools

include("util.jl")

macro ifv1p4(expr)
    if VERSION >= v"1.4"
        thunk = esc(:(()->($expr)))
        quote
            $thunk()
        end
    end
end

@testset "@tspawnat" begin

    @testset "@normal operation" begin
        obj = TestObj(0)
        function fn!(obj)
            sleep(0.1)
            obj.data = Threads.threadid()
        end
        task = @tspawnat Threads.nthreads() fn!(obj)
        @test obj.data == 0
        wait(task)
        @test obj.data == Threads.nthreads()
    end

    @ifv1p4 begin
        @testset "interpolation" begin
            function foo(x)
                sleep(0.01)
                return x
            end

            x = 1
            expect_sum = 3
            t1 = @tspawnat max(1, Threads.nthreads()) foo($x)
            x += 1
            t2 = @tspawnat max(1, Threads.nthreads()-1) foo($x)

            test_sum = fetch(t1) + fetch(t2)
            @test expect_sum == test_sum
        end
    end

    @testset "@out of bounds" begin
        @test_throws AssertionError task = @tspawnat Threads.nthreads()+1 randn()
        @test_throws AssertionError task = @tspawnat 0 randn()
    end

    end
end # module


text:
using Plots

function plot_lanczos(model)
    plotly()
    i = 0;
    for set in model.p.C.sets
        if isa(set, COSMO.PsdConeTriangleLanczos)
            p1 = plot(hcat(set.subspace_dim_history,
                    set.n/2*ones(length(set.subspace_dim_history))),
                    title=string("Set #", i,":  Subspace sizes"),
                    label=["subspace size" "n/2"],
                    # line=:stem
                )
            p2 = plot(hcat(max.(set.residual_history, 1e-6),
                    max.(set.λ_rem_history, 1e-6)),
                    title=string("Set #", i,":  Residual history"),
                    yaxis=:log,
                    label=["total residual" "only from ignored subspace"],
                    # line=:stem
                )
            p3 = scatter(set.λ_rem_multiplications,
                    title=string("Set #", i,":  ARPACK multiplications"),
                    markersize = 2,
                )
            plot(p1, p2, p3)
            gui()
        end
        i = i + 1
	end
end


text:
using PyCall,DataFrames,RCall
unshift!(PyVector(pyimport("sys")["path"]), "")
@pyimport ncdata
len = length
filename = "edwards11.nc"

data = ncdata.get(filename)
specs = names!(DataFrame(data["spec"]),[Symbol(i)for i in data["sc"]])
rates = names!(DataFrame(data["rate"]),[Symbol(i)for i in data["rc"]])
rates = rates[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["rc"])]]
specs = specs[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["sc"])]]
rates = rates[7:len(rates)]
rates = rates[:,[Symbol(i) for i in filter(x->sum(rates[x])>0,names(rates))]]
joint = "\n"*(join(names(rates),"\n"))*"\n"

products= [split(i[4:len(i)],"+") for i in matchall(r"-->([A-z0-9+]*)",joint)]
reactants = [split(i[2:len(i)-1],"+") for i in matchall(r"\n([A-z0-9+]{1,60})([-->]{0,1})",joint)]

node_set = Set([ i[2:len(i)-1] for i in matchall(r"[\n+>-]{1}[A-z0-9]+[\n+>-]{1}",joint) ])
nodes = [i for i in node_set]
dict = Dict( [[i,x] for (x,i) in enumerate(nodes)])

flux = []
edges = []
for i in 1:length(reactants)
  rcol=[]
  #println(i)
  for j in reactants[i]
    coeff = match( r"([\d\s\.]*)(\D[\d\D]*)", j)
    dummy = specs[Symbol(coeff[2])]
    try
      push!(rcol,parse(Float32,coeff[1]))
    catch
      push!(rcol,1.0f0)
    end

    prod = 1
    for k in rcol
      prod*=k
    end

    push!(flux, prod*rates[i])
    # make edges array
    for l in products[i]
        push!(edges,[match( r"([\d\s\.]*)(\D[\d\D]*)", j)[2],match( r"([\d\s\.]*)(\D[\d\D]*)", l)[2],i])
    end

  end
end

normalise(x) = x/norm(x)

function sorteqn(x,delim)
    r,p = split(string(x),delim)
    r = join(sort(split(r,'+')),'+')
    p = join(sort(split(p,'+')),'+')
    return join([r,p],delim)
end

rateeqn = [Symbol(sorteqn(i,"-->")) for i in names(rates)]

grouprecorded = [dict[Symbol(sorteqn(i,"-->"))] for i in names(rates)]

################READ CATEGORIES##################
dict = Dict()
open("edgecat.json", "r") do f
    global dict
    dict=JSON.parse(readstring(f))  # parse and transform data
end

dict = Dict(Symbol( sorteqn( replace(key,"=","-->") , "-->" )) => value for (key, value) in dict)

reactiontypes= [try dict[i]; catch err; "missing" end  for i in rateeqn]

filter(i->ismatch(r".*NC7H16.*", i), [string(i) for i in keys(dict)])

smiles =DataFrame(readtable("smiles_mined.csv"))

t = 325

links = filter(i -> flux[i][t]>0 , 1:len(flux))
tflux = [log10(flux[i][t]) for i in links]
weight =  tflux - minimum(tflux)
weight = 1 - weight/maximum(weight)

#1+normalise(tflux)

newflux = [flux[i][t] for i in 1:len(flux)]
counter = 0
for i in links
  newflux[i] = weight[counter+=1]
end

edge = filter(i -> newflux[i[3]]>0 , edges)
source = [i[1] for i in edge]
target=[i[2] for i in edge]
weighted = [newflux[i[3]]+0.0001 for i in edge]
grouping  = [reactiontypes[i[3]] for i in edge]

JSONdata =  [(i[1],i[2],newflux[i[3]],reactiontypes[i[3]]) for i  in edge]

max_stars_repo_name: wolfiex/RopaJL

max_issues_repo_name: wolfiex/RopaJL

max_forks_repo_name: wolfiex/RopaJL


text:
# Autogenerated wrapper script for GoogleTest_jll for i686-linux-musl-cxx03
export libgtest, libgtest_main

JLLWrappers.@generate_wrapper_header("GoogleTest")
JLLWrappers.@declare_library_product(libgtest, "libgtest.so")
JLLWrappers.@declare_library_product(libgtest_main, "libgtest_main.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgtest,
        "lib/libgtest.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgtest_main,
        "lib/libgtest_main.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Inspired by Python's gettext, but in a Julian way.

"""
    Translation

This holds the contents of all localizations.
"""
struct Translation
    content::String # original English content
    translation::Dict{Symbol, Any}
end

Translation(content::String) = Translation(content, Dict{Symbol, Any}())

Base.getindex(text::Translation, locale::Symbol) = getindex(text.translation, locale)
Base.setindex!(text::Translation, trans, locale::Symbol) = setindex!(text.translation, trans, locale)
Base.show(io::IO, t::Translation) = Base.show(io, t.content)

# make gettext people happy
translate(text::Translation, locale::Symbol) = text.translation[locale]


text:
# ---
# title: 1671. Minimum Number of Removals to Make Mountain Array
# id: problem1671
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/minimum-number-of-removals-to-make-mountain-array/description/>
# hidden: true
# ---
#
# You may recall that an array `arr` is a **mountain array** if and only if:
#
#   * `arr.length >= 3`
#   * There exists some index `i` ( **0-indexed** ) with `0 < i < arr.length - 1` such that:
#     * `arr[0] < arr[1] < ... < arr[i - 1] < arr[i]`
#     * `arr[i] > arr[i + 1] > ... > arr[arr.length - 1]`
#
# Given an integer array `nums`, return _the **minimum** number of elements
# to remove to make _`nums __` ___a **mountain array**._
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,3,1]
#     Output: 0
#     Explanation: The array itself is a mountain array so we do not need to remove any elements.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [2,1,1,5,6,2,3,1]
#     Output: 3
#     Explanation: One solution is to remove the elements at indices 0, 1, and 5, making the array nums = [1,5,6,3,1].
#
#
# **Example 3:**
#
#
#
#     Input: nums = [4,3,2,1,1,2,3,1]
#     Output: 4
#
#
# **Example 4:**
#
#
#
#     Input: nums = [1,2,3,4,4,3,2,1]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `3 <= nums.length <= 1000`
#   * `1 <= nums[i] <= 109`
#   * It is guaranteed that you can make a mountain array out of `nums`.
#
#
## @lc code=start
using LeetCode

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


text:
using MacroTools

function clean_kws(args)
  map(args) do arg
    !isexpr(arg,:(=)) ? arg : Expr(:kw,arg.args...)
  end
end

function isdocstr_block(docstr)
  if isexpr(docstr,:block)
    filtered = filter(x -> !isexpr(x,:line),docstr.args)
    length(filtered) == 1 &&
    filtered[1] isa String
  end
end

function warn_doc(name)
  """

!!! warning "Do not call inside a package"

    Do not call @$(name) inside of a package or in tests. It should
    only be used in one-off scripts. If Weber$(name) is not currently installed
    it will be installed by this macro using `Pkg.add` which can lead to
    problems when called in packages or tests.

    If you want to include an extension without this behavior you can call
    @$(name)_safe which mimics @$(name) except that it will never call
    `Pkg.add`.

"""
end

"""
    @extension [Symbol] begin
      [docstring...]
    end

Registers a given Weber extension. This creates a macro called `@[Symbol]` which
imports `Weber[Symbol]` and calls `Weber[Symbol].InitExtension`, with the given
arguments. InitExtension should return either `nothing` or an extension
object.

The doc string is used to document the usage of the extension, and
should normally include a link to the website of a julia package for the
extension.
"""
macro extension(symbol,docstr)
  if (symbol isa Symbol) & isdocstr_block(docstr)
    mod = Symbol("Weber"string(symbol))
    args = gensym("args")
    safe_doc = "See documentation for @$(string(symbol))."
    expr = quote
      export $(Symbol("@"string(symbol))), $(Symbol("@"string(symbol)"_safe"))

      @doc $(filter(x -> !isexpr(x,:line),docstr.args)[1]*
             warn_doc(string(symbol))) ->
      macro $symbol($args...)
        init_call = :($$(Expr(:quote,mod)).InitExtension())
        append!(init_call.args,Weber.clean_kws($args))
        quote
          try
            $(Expr(:import,$(Expr(:quote,mod))))
          catch
            Pkg.add($$(string(mod)))
            $(Expr(:import,$(Expr(:quote,mod))))
          end
          $(esc(init_call))
        end
      end

      @doc $(safe_doc) ->
      macro $(Symbol(string(symbol)"_safe"))($args...)
        init_call = :($$(Expr(:quote,mod)).InitExtension())
        append!(init_call.args,Weber.clean_kws($args))
        quote
          $(Expr(:import,$(Expr(:quote,mod))))
          $(esc(init_call))
        end
      end
    end
    esc(expr)
  else
    error("Unexpected syntax: @extension expects a sybmol followed by a",
          " block containing a doc string.")
  end
end


text:
# Testing a pretrained normalizing flow for seismic denosing
# Authors: Ali Siahkoohi, alisk@gatech.edu
#          Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
#          Philipp Witte, pwitte3@gatech.edu
# Date: September 2020
"""
Example:
julia test_hint_denoising.jl --batchsize 256 --max_epoch 5
"""

using DrWatson
@quickactivate :FastApproximateInference

using InvertibleNetworks
using JLD
using JLD2
using Random
using Statistics
using ArgParse
using BSON
using Printf
using PyPlot
using Seaborn
using Flux: gpu, cpu
using Flux.Data: DataLoader

set_style("whitegrid")
rc("font", family="serif", size=13)
font_prop = matplotlib.font_manager.FontProperties(
    family="serif",
    style="normal",
    size=13
)
sfmt=matplotlib.ticker.ScalarFormatter(useMathText=true)
sfmt.set_powerlimits((0, 0))

# Random seed
# Random.seed!(19)

s = ArgParseSettings()
@add_arg_table! s begin
    "--max_epoch"
        help = "Maximum # of epochs"
        arg_type = Int
        default = 50
    "--epoch"
        help = "Epoch to load net params"
        arg_type = Int
        default = -1
    "--lr"
        help = "starting learning rate"
        arg_type = Float32
        default = 1f-3
    "--lr_step"
        help = "lr scheduler applied at evey lr_step epoch"
        arg_type = Int
        default = 1
    "--batchsize"
        help = "batch size"
        arg_type = Int
        default = 4
    "--n_hidden"
        help = "# hidden channels"
        arg_type = Int
        default = 64
    "--depth"
        help = "depth of the network"
        arg_type = Int
        default = 12
    "--sigma"
        help = "noise std"
        arg_type = Float32
        default = 2f-2
    "--sim_name"
        help = "simulation name"
        arg_type = String
        default = "seismic-denoising-256x256-0.3CS"
end
parsed_args = parse_args(s)

max_epoch = parsed_args["max_epoch"]
lr = parsed_args["lr"]
lr_step = parsed_args["lr_step"]
batchsize = parsed_args["batchsize"]
n_hidden = parsed_args["n_hidden"]
depth = parsed_args["depth"]
sigma = parsed_args["sigma"]
sim_name = parsed_args["sim_name"]

if parsed_args["epoch"] == -1
    parsed_args["epoch"] = parsed_args["max_epoch"]
end
epoch = parsed_args["epoch"]

# Define raw data directory
mkpath(datadir("training-data"))
data_path = datadir("training-data", "seismic_samples_256_by_256_num_10k.jld")
label_path = datadir("training-data", "seismic_samples_256_by_256_num_10k_labels.jld")

# Download the dataset into the data directory if it does not exist
if isfile(data_path) == false
    run(`wget https://www.dropbox.com/s/vapl62yhh8fxgwy/'
        'seismic_samples_256_by_256_num_10k.jld -q -O $data_path`)
end
if isfile(label_path) == false
    run(`wget https://www.dropbox.com/s/blxkh6tdszudlcq/'
        'seismic_samples_256_by_256_num_10k_labels.jld -q -O $label_path`)
end

# Load seismic images and create training and testing data
JLD.jldopen(data_path, "r") do file
    X_orig = read(file, "X")

    labels = JLD.jldopen(label_path, "r")["labels"][:]
    idx = findall(x -> x == 1, labels)

    X_orig = X_orig[:, :, :, idx]
    global nx, ny, nc, nsamples = size(X_orig)

    # Whiten data
    X_orig = whiten_input(X_orig)

    # Load observed data
    Y_orig = wload(
        datadir("testing-data", "noisy_seismic_samples_256_by_256_num_10k.jld2")
    )["Y_orig"]

    # Load split in training/testing
    train_idx = wload(
        datadir("testing-data", "noisy_seismic_samples_256_by_256_num_10k.jld2")
    )["train_idx"]

    # Split in training - testing
    global test_size = 5
    global test_idx = shuffle(setdiff(1:nsamples, train_idx))[1:5]

    # Dimensions after wavelet squeeze to increase no. of channels
    global nx = Int(nx/2)
    global ny = Int(ny/2)
    global n_in = Int(nc*4)

    # Apply wavelet squeeze (change dimensions to -> n1/2 x n2/2 x nc*4 x ntrain)
    global X_test = zeros(Float32, nx, ny, n_in, test_size)
    global Y_test = zeros(Float32, nx, ny, n_in, test_size)
    for (k, j) in enumerate(test_idx)
        X_test[:, :, :, k:k] = wavelet_squeeze(X_orig[:, :, :, j:j])
        Y_test[:, :, :, k:k] = wavelet_squeeze(Y_orig[:, :, :, j:j])
    end
end

# Create network
CH = NetworkConditionalHINT(
    n_in, n_hidden, depth, k1=3, k2=3, p1=1, p2=1
)

# Loading the experiment—only network weights and training loss
Params, fval = load_experiment(parsed_args)
put_params!(CH, convert(Array{Any,1}, Params))

Zx = randn(Float32, nx, ny, n_in, test_size)
Zy = randn(Float32, nx, ny, n_in, test_size)

# Precited model and data
X_, Y_ = CH.inverse(Zx, Zy)

# Now select single fixed sample from all Ys
X_fixed = X_test[:, :, :, 1:1]
Y_fixed = Y_test[:, :, :, 1:1]
Zy_fixed = CH.forward_Y(Y_fixed)

# Draw new Zx, while keeping Zy fixed
n_samples = 1000
X_post = zeros(Float32, nx, ny, n_in, n_samples)
CH = CH |> gpu

test_loader = DataLoader(
    (randn(Float32, nx, ny, n_in, n_samples), repeat(Zy_fixed, 1, 1, 1, n_samples)),
    batchsize=8, shuffle=false
)

for (itr, (X, Y)) in enumerate(test_loader)
    counter = (itr - 1)*size(X)[4] + 1

    X = X |> gpu
    Y = Y |> gpu

    X_post[:, :, :, counter:(counter+size(X)[4]-1)] = (CH.inverse(X, Y)[1] |> cpu)
    @printf("Sampling: [%d/%d]\r", counter, n_samples)
end

X_post = wavelet_unsqueeze(X_post)

# Some stats
X_post_mean = mean(X_post; dims=4)
X_post_std = std(X_post; dims=4)

# Unsqueeze all tensors
X_test = wavelet_unsqueeze(X_test)
Y_test = wavelet_unsqueeze(Y_test)

X_ = wavelet_unsqueeze(X_)
Y_ = wavelet_unsqueeze(Y_)
Zx = wavelet_unsqueeze(Zx)
Zy = wavelet_unsqueeze(Zy)

X_fixed = wavelet_unsqueeze(X_fixed)
Y_fixed = wavelet_unsqueeze(Y_fixed)
Zy_fixed = wavelet_unsqueeze(Zy_fixed)

save_dict = @strdict max_epoch epoch lr lr_step batchsize n_hidden depth sim_name
save_path = plotsdir(sim_name, savename(save_dict; digits=6))

# Training loss
fig = figure("training logs", figsize=(7, 2.5))
plot(fval, color="#d48955")
title("Negative log-likelihood")
ylabel("Training objective")
xlabel("Iterations")
safesave(joinpath(save_path, "log.png"), fig)
close(fig)

# Plot the true model
fig = figure("x", figsize=(5, 5))
imshow(
    X_fixed[:, :, 1, 1], vmin=-4.0, vmax=4.0, aspect=1, cmap="Greys", resample=true,
    interpolation="lanczos", filterrad=1
)
title(L"True model: $\mathbf{x} \sim \widehat{\pi}_x (\mathbf{x})$")
colorbar(fraction=0.047, pad=0.01, format=sfmt)
grid(false)
ax = plt.gca()
ax.axes.xaxis.set_visible(false)
ax.axes.yaxis.set_visible(false)
safesave(joinpath(save_path, "true_model" * string(test_idx[1]) * ".png"), fig)
close(fig)

# Plot the observed data
fig = figure("y", figsize=(5, 5))
imshow(
    Y_fixed[:, :, 1, 1], vmin=-4.0, vmax=4.0, aspect=1, cmap="Greys", resample=true,
    interpolation="lanczos", filterrad=1
)
title(L"Observed data: $\mathbf{y} \sim \widehat{\pi}_y (\mathbf{y})$")
colorbar(fraction=0.047, pad=0.01, format=sfmt)
grid(false)
ax = plt.gca()
ax.axes.xaxis.set_visible(false)
ax.axes.yaxis.set_visible(false)
safesave(joinpath(save_path, "observed_data" * string(test_idx[1]) * ".png"), fig)
close(fig)

# Plot the conditional mean estimate
fig = figure("x_cm", figsize=(5, 5))
imshow(
    X_post_mean[:, :, 1, 1], vmin=-4.0, vmax=4.0, aspect=1, cmap="Greys", resample=true,
    interpolation="lanczos", filterrad=1
)
title("Conditional mean")
colorbar(fraction=0.047, pad=0.01, format=sfmt)
grid(false)
ax = plt.gca()
ax.axes.xaxis.set_visible(false)
ax.axes.yaxis.set_visible(false)
safesave(joinpath(save_path, "conditional_mean" * string(test_idx[1]) * ".png"), fig)
close(fig)

# Plot the pointwise standard deviation
fig = figure("x_std", figsize=(5, 5))
imshow(
    X_post_std[:, :, 1, 1], vmin=.30,
    vmax=.42, aspect=1, cmap="OrRd", resample=true,
    interpolation="lanczos", filterrad=1
)
title("Pointwise STD")
colorbar(fraction=0.047, pad=0.01, format=sfmt)
grid(false)
ax = plt.gca()
ax.axes.xaxis.set_visible(false)
ax.axes.yaxis.set_visible(false)
safesave(joinpath(save_path, "pointwise_std" * string(test_idx[1]) * ".png"), fig)
close(fig)

fig = figure("err", figsize=(5, 5))
imshow(
    2f0.*(X_fixed[:, :, 1, 1] - X_post_mean[:, :, 1, 1])./(abs.(X_fixed[:, :, 1, 1])  + abs.(X_post_mean[:, :, 1, 1])), vmin=-2, vmax=2, aspect=1, cmap="twilight_shifted", resample=true,
    interpolation="lanczos", filterrad=1
)
title("Estimation error — accelerated")
colorbar(fraction=0.047, pad=0.01, format=sfmt)
grid(false)
ax = plt.gca()
ax.axes.xaxis.set_visible(false)
ax.axes.yaxis.set_visible(false)
safesave(joinpath(save_path, "conditional_mean_err.png"), fig)
close(fig)


text:
# ffts
function ChainRulesCore.frule((_, Δx, _), ::typeof(fft), x::AbstractArray, dims)
    y = fft(x, dims)
    Δy = fft(Δx, dims)
    return y, Δy
end
function ChainRulesCore.rrule(::typeof(fft), x::AbstractArray, dims)
    y = fft(x, dims)
    project_x = ChainRulesCore.ProjectTo(x)
    function fft_pullback(ȳ)
        x̄ = project_x(bfft(ChainRulesCore.unthunk(ȳ), dims))
        return ChainRulesCore.NoTangent(), x̄, ChainRulesCore.NoTangent()
    end
    return y, fft_pullback
end

function ChainRulesCore.frule((_, Δx, _), ::typeof(rfft), x::AbstractArray{<:Real}, dims)
    y = rfft(x, dims)
    Δy = rfft(Δx, dims)
    return y, Δy
end
function ChainRulesCore.rrule(::typeof(rfft), x::AbstractArray{<:Real}, dims)
    y = rfft(x, dims)

    # compute scaling factors
    halfdim = first(dims)
    d = size(x, halfdim)
    n = size(y, halfdim)
    scale = reshape(
        [i == 1 || (i == n && 2 * (i - 1) == d) ? 1 : 2 for i in 1:n],
        ntuple(i -> i == first(dims) ? n : 1, Val(ndims(x))),
    )

    project_x = ChainRulesCore.ProjectTo(x)
    function rfft_pullback(ȳ)
        x̄ = project_x(brfft(ChainRulesCore.unthunk(ȳ) ./ scale, d, dims))
        return ChainRulesCore.NoTangent(), x̄, ChainRulesCore.NoTangent()
    end
    return y, rfft_pullback
end

function ChainRulesCore.frule((_, Δx, _), ::typeof(ifft), x::AbstractArray, dims)
    y = ifft(x, dims)
    Δy = ifft(Δx, dims)
    return y, Δy
end
function ChainRulesCore.rrule(::typeof(ifft), x::AbstractArray, dims)
    y = ifft(x, dims)
    invN = normalization(y, dims)
    project_x = ChainRulesCore.ProjectTo(x)
    function ifft_pullback(ȳ)
        x̄ = project_x(invN .* fft(ChainRulesCore.unthunk(ȳ), dims))
        return ChainRulesCore.NoTangent(), x̄, ChainRulesCore.NoTangent()
    end
    return y, ifft_pullback
end

function ChainRulesCore.frule((_, Δx, _, _), ::typeof(irfft), x::AbstractArray, d::Int, dims)
    y = irfft(x, d, dims)
    Δy = irfft(Δx, d, dims)
    return y, Δy
end
function ChainRulesCore.rrule(::typeof(irfft), x::AbstractArray, d::Int, dims)
    y = irfft(x, d, dims)

    # compute scaling factors
    halfdim = first(dims)
    n = size(x, halfdim)
    invN = normalization(y, dims)
    twoinvN = 2 * invN
    scale = reshape(
        [i == 1 || (i == n && 2 * (i - 1) == d) ? invN : twoinvN for i in 1:n],
        ntuple(i -> i == first(dims) ? n : 1, Val(ndims(x))),
    )

    project_x = ChainRulesCore.ProjectTo(x)
    function irfft_pullback(ȳ)
        x̄ = project_x(scale .* rfft(real.(ChainRulesCore.unthunk(ȳ)), dims))
        return ChainRulesCore.NoTangent(), x̄, ChainRulesCore.NoTangent(), ChainRulesCore.NoTangent()
    end
    return y, irfft_pullback
end

function ChainRulesCore.frule((_, Δx, _), ::typeof(bfft), x::AbstractArray, dims)
    y = bfft(x, dims)
    Δy = bfft(Δx, dims)
    return y, Δy
end
function ChainRulesCore.rrule(::typeof(bfft), x::AbstractArray, dims)
    y = bfft(x, dims)
    project_x = ChainRulesCore.ProjectTo(x)
    function bfft_pullback(ȳ)
        x̄ = project_x(fft(ChainRulesCore.unthunk(ȳ), dims))
        return ChainRulesCore.NoTangent(), x̄, ChainRulesCore.NoTangent()
    end
    return y, bfft_pullback
end

function ChainRulesCore.frule((_, Δx, _, _), ::typeof(brfft), x::AbstractArray, d::Int, dims)
    y = brfft(x, d, dims)
    Δy = brfft(Δx, d, dims)
    return y, Δy
end
function ChainRulesCore.rrule(::typeof(brfft), x::AbstractArray, d::Int, dims)
    y = brfft(x, d, dims)

    # compute scaling factors
    halfdim = first(dims)
    n = size(x, halfdim)
    scale = reshape(
        [i == 1 || (i == n && 2 * (i - 1) == d) ? 1 : 2 for i in 1:n],
        ntuple(i -> i == first(dims) ? n : 1, Val(ndims(x))),
    )

    project_x = ChainRulesCore.ProjectTo(x)
    function brfft_pullback(ȳ)
        x̄ = project_x(scale .* rfft(real.(ChainRulesCore.unthunk(ȳ)), dims))
        return ChainRulesCore.NoTangent(), x̄, ChainRulesCore.NoTangent(), ChainRulesCore.NoTangent()
    end
    return y, brfft_pullback
end

# shift functions
function ChainRulesCore.frule((_, Δx, _), ::typeof(fftshift), x::AbstractArray, dims)
    y = fftshift(x, dims)
    Δy = fftshift(Δx, dims)
    return y, Δy
end
function ChainRulesCore.rrule(::typeof(fftshift), x::AbstractArray, dims)
    y = fftshift(x, dims)
    project_x = ChainRulesCore.ProjectTo(x)
    function fftshift_pullback(ȳ)
        x̄ = project_x(ifftshift(ChainRulesCore.unthunk(ȳ), dims))
        return ChainRulesCore.NoTangent(), x̄, ChainRulesCore.NoTangent()
    end
    return y, fftshift_pullback
end

function ChainRulesCore.frule((_, Δx, _), ::typeof(ifftshift), x::AbstractArray, dims)
    y = ifftshift(x, dims)
    Δy = ifftshift(Δx, dims)
    return y, Δy
end
function ChainRulesCore.rrule(::typeof(ifftshift), x::AbstractArray, dims)
    y = ifftshift(x, dims)
    project_x = ChainRulesCore.ProjectTo(x)
    function ifftshift_pullback(ȳ)
        x̄ = project_x(fftshift(ChainRulesCore.unthunk(ȳ), dims))
        return ChainRulesCore.NoTangent(), x̄, ChainRulesCore.NoTangent()
    end
    return y, ifftshift_pullback
end


text:
export
  ket,
  bra,
  ketbra,
  proj,
  res,
  unres,
  fourier_matrix

"""
    ket(index, size)

Return `index`-th base (column) vector in the `size`-dimensional vector space.
To be consistent with Julia indexing, `index` = 1, 2, ..., `size`.

# Examples

```jldoctest; setup = :(using QSWalk)
julia> ket(1, 2)
2-element SparseArrays.SparseVector{Int64,Int64} with 1 stored entry:
  [1]  =  1
```
"""
function ket(index::Int, size::Int)

  @argumentcheck size > 0 "vector size must be positive"
  @assert 1 <=  index <=  size "index must be positive and not bigger than vector size"

  ret = spzeros(Int, size)
  ret[index]  =  1
  ret
end

"""
    bra(index, size)

Return `index`-th base row vector in the `size`-dimensional vector space, with
`index` = 1, 2, ..., `size`.

# Examples

```jldoctest; setup = :(using QSWalk)
julia> bra(1, 2)
1×2 LinearAlgebra.Adjoint{Int64,SparseArrays.SparseVector{Int64,Int64}}:
 1  0
```
"""
function bra(index::Int, size::Int)
  ket(index, size)'
end

"""
    ketbra(irow, icol, size)

Return matrix acting on `size`-dimensional vector space, `irow`,
`icol` = 1, 2, ..., `size`. The matrix consists of single non-zero element
equal to one, located at position (`irow`, `icol`).

# Examples

```jldoctest; setup = :(using QSWalk)
julia> ketbra(1, 2, 3)
3×3 SparseArrays.SparseMatrixCSC{Int64,Int64} with 1 stored entry:
  [1, 2]  =  1

```
"""
function ketbra(irow::Int, icol::Int, size::Int)
  ket(irow, size)*bra(icol, size)
end

"""
    proj(index, size)

Return projector onto `index`-th base vector in `size`-dimensional vector space,
with `index` = 1, 2, ..., `size`. This is equivalent to `ketbra(index, index,
size)`.

```jldoctest; setup = :(using QSWalk)
julia> proj(1, 2)
2×2 SparseArrays.SparseMatrixCSC{Int64,Int64} with 1 stored entry:
  [1, 1]  =  1
```
"""
function proj(index::Int, size::Int)
  ketbra(index, index, size)
end

"""
    proj(vector)

Return projector onto the subspace spanned by vector `vector`.

# Examples
```jldoctest; setup = :(using QSWalk)
julia> v = 1/sqrt(2) * (ket(1, 3)+ket(3, 3))
3-element SparseArrays.SparseVector{Float64,Int64} with 2 stored entries:
  [1]  =  0.707107
  [3]  =  0.707107

julia> proj(v)
3×3 SparseArrays.SparseMatrixCSC{Float64,Int64} with 4 stored entries:
  [1, 1]  =  0.5
  [3, 1]  =  0.5
  [1, 3]  =  0.5
  [3, 3]  =  0.5
```
"""
function proj(vector::AbstractVector)
  vector*vector'
end

"""

    res(matrix)

Return vectorization of the `matrix` in the row order. This is equivalent to
`Base.vec(transpose(matrix))`.

# Examples

```jldoctest; setup = :(using QSWalk)
julia> M = reshape(1:9, (3, 3))'*1im
3×3 Array{Complex{Int64},2}:
 0+1im  0+2im  0+3im
 0+4im  0+5im  0+6im
 0+7im  0+8im  0+9im

julia> v = res(M)
9-element reshape(::LinearAlgebra.Transpose{Complex{Int64},Array{Complex{Int64},2}}, 9) with eltype Complex{Int64}:
 0 + 1im
 0 + 2im
 0 + 3im
 0 + 4im
 0 + 5im
 0 + 6im
 0 + 7im
 0 + 8im
 0 + 9im

julia> res(unres(v)) == v
true
```
"""
function res(matrix::AbstractMatrix)
  Base.vec(transpose(matrix))
end

"""

    unres(vector)

Return square matrix with elements from `vector`. The `vector` is expected to
have perfect square number of arguments.

# Examples
```jldoctest; setup = :(using QSWalk)
julia> v = collect(1:9)*im
9-element Array{Complex{Int64},1}:
 0 + 1im
 0 + 2im
 0 + 3im
 0 + 4im
 0 + 5im
 0 + 6im
 0 + 7im
 0 + 8im
 0 + 9im

julia> unres(v)
3×3 LinearAlgebra.Transpose{Complex{Int64},Array{Complex{Int64},2}}:
 0+1im  0+2im  0+3im
 0+4im  0+5im  0+6im
 0+7im  0+8im  0+9im

julia> res(unres(v)) == v
true
```
"""
function unres(vector::AbstractVector)
  dim = floor(Int64, sqrt(length(vector)))
  @argumentcheck dim*dim == length(vector) "Expected vector with perfect square number of elements."

  transpose(reshape(vector, (dim, dim)))
end

"""

    fourier_matrix(size)

Returns Fourier matrix of size `size`×`size`.

# Examples

```jldoctest; setup = :(using QSWalk)
julia> fourier_matrix(1)
1×1 SparseArrays.SparseMatrixCSC{Complex{Float64},Int64} with 1 stored entry:
  [1, 1]  =  1.0+0.0im

julia> fourier_matrix(2)
2×2 SparseArrays.SparseMatrixCSC{Complex{Float64},Int64} with 4 stored entries:
  [1, 1]  =  1.0+0.0im
  [2, 1]  =  1.0+0.0im
  [1, 2]  =  1.0+0.0im
  [2, 2]  =  -1.0+1.22465e-16im

```
"""
function fourier_matrix(size::Int)
  @argumentcheck size>0 "Size of the matrix needs to be positive"
  sparse([exp(2im*π*(i-1)*(j-1)/size) for i = 1:size, j = 1:size])
end


text:
include(normpath(joinpath(@__DIR__, "IOM", "src", "models", "EMOM", "EMOM.jl")))
using DataStructures
using NCDatasets
using MPI
using Formatting
using ArgParse, JSON

function runOneCmd(cmd)
    println(">> ", string(cmd))
    run(cmd)
end

function pleaseRun(cmd)
    if isa(cmd, Array)
        for i = 1:length(cmd)
            runOneCmd(cmd[i])
        end
    else
        runOneCmd(cmd)
    end
end

println("""
This program produces (1) init.jld and (2) Nz_bot.nc file given a config file.
Nz_bot.nc is generated through running make_Nz_bot.jl or make_Nz_bot_SOM.jl.

The NaNs of TEMP or SALT variables are assumed masked and thus topo file is
constructed.

This program also checks if the given initial profiles are consistent, meaning
their NaNs are located exactly at the same grids.
""")

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin

        "--init-profile-TEMP"
            help = "Initial ocean profile"
            arg_type = String
            required = true

        "--init-profile-SALT"
            help = "Initial ocean profile"
            arg_type = String
            required = true

        "--init-profile-HMXL"
            help = "Initial ocean profile"
            arg_type = String
            required = true

        "--HMXL-unit"
            help = "Unit of HMXL. By default I will try to convert it to meters."
            arg_type = String
            required = true

        "--config-file"
            help = "config TOML file"
            arg_type = String
            required = true

        "--SOM"
            help = "If set then data blow Nz_bot will be removed."
            action = :store_true

    end

    return parse_args(s)
end

parsed = DataStructures.OrderedDict(parse_args(ARGS, s))
parsed = parse_commandline()
JSON.print(parsed,4)

MPI.Init()

println("Processing data...")

using TOML
config = TOML.parsefile(parsed["config-file"])

domain_file = config["MODEL_CORE"]["domain_file"]

init_file  = config["MODEL_MISC"]["init_file"]
topo_file  = config["MODEL_CORE"]["topo_file"]

Nz = length(config["MODEL_CORE"]["z_w"]) - 1 # Layers used. Thickness ≈ 503m

Dataset(domain_file, "r") do ds
    global Nx = ds.dim["ni"]
    global Ny = ds.dim["nj"]
    global mask = ds["mask"][:]
end

Dataset(parsed["init-profile-TEMP"], "r") do ds
    global TEMP  = permutedims(nomissing(ds["TEMP"][:, :, 1:Nz, 1],  NaN), [3, 1, 2])
end

Dataset(parsed["init-profile-SALT"], "r") do ds
    global SALT  = permutedims(nomissing(ds["SALT"][:, :, 1:Nz, 1],  NaN), [3, 1, 2])
end

Dataset(parsed["init-profile-HMXL"], "r") do ds
    global HMXL  = nomissing(ds["HMXL"][:, :, 1],  NaN)
    if parsed["HMXL-unit"] == "m"
        println("It is already in meters.")
    elseif parsed["HMXL-unit"] == "cm"
        println("It is cm! Convert it...")
        HMXL ./= 100.0
    else
        thorw(ErrorException("Unknown unit: $(parsed["HMXL-unit"])"))
    end

    if any(HMXL .< 0)
        thorw(ErrorException("HMXL cannot be negative."))
    end
end

Nx, Ny = size(TEMP)[2:3]
Nz_bot = zeros(Int64, Nx, Ny)
mask_T = zeros(Float64, Nz, Nx, Ny)
mask_T[isfinite.(TEMP)] .= 1.0
valid_idx = mask_T .== 1.0

println("Making Nz_bot...")
if ! parsed["SOM"]
    println("This is not SOM.")
    Nz_bot .= sum(mask_T, dims=1)[1, :, :]
else
    println("This is SOM.")
    for i=1:Nx, j=1:Ny
        h = HMXL[i, j]
        if isnan(h)
            Nz_bot[i, j] = 0
        else
            Nz_bot[i, j] = findlast(z_w_top .>= - h)  # Use ">=" to avoid h = 0 that makes Nz = 0
        end
    end

end

# Test if Nz_bot is consistent with the mask
mask_from_Nz_bot = copy(Nz_bot)
mask_from_Nz_bot[mask_from_Nz_bot .!= 0] .= 1

println("Check if SOM's Nz_bot is consistent with the mask...")
if all(mask_from_Nz_bot .== mask)
    println("Yes they are consistent.")
else
    throw(ErrorException("Nz_bot and mask not consistent"))
end

println("Check if the NaNs of TEMP, SALT are both consistent")
if any(isfinite.(TEMP) .!= isfinite.(SALT))
    throw(ErrorException("Locations of NaNs of TEMP and SALT mismatch. Please check."))
end

println(format("Output file: {}.", topo_file))
Dataset(topo_file, "c") do ds

    defDim(ds, "Nx", Nx)
    defDim(ds, "Ny", Ny)
    defVar(ds, "Nz_bot", Nz_bot, ("Nx", "Ny", ), ; attrib = Dict(
        "long_name" => "z-grid idx of the deepest cell",
        "units"     => "none",
    ))
end

println("Create a model to save initial condition.")

ev = EMOM.Env(config["MODEL_CORE"])
mb = EMOM.ModelBlock(ev; init_core=false)

mb.fi.sv[:TEMP][valid_idx] .= TEMP[valid_idx]
mb.fi.sv[:SALT][valid_idx] .= SALT[valid_idx]
mb.fi.HMXL[1, :, :] .= HMXL

println(format("Output file: {}.", init_file))

EMOM.takeSnapshot(DateTimeNoLeap(1,1,1), mb, init_file)

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
"""
Created on Wed Apr 22 2020

@author: Yoann Pradat

Implementation of the alpha NMF algorithm under various constraints as described in Chapter 3 of
Andrzej CICHOCKI, et. al., Nonnegative Matrix and Tensor Factorizations. John Wiley & Sons, 2009.
"""

"""
    _normalize(X[, p])

Column-wise normalization.
"""
function _normalize(X::Matrix{T}; p::Integer=1) where T <: Real
    norms = sum(X .^ p, dims=1) .^ (1. / p)
    X .* repeat(norms .^ -1, size(X, 1), 1)
end

"""
    _update_nmf_α!(V, V_a, W, H, params)

Update the H and W matrices in that order using the update rules using algorithm 3.3 page 147
Andrzej CICHOCKI, et. al., Nonnegative Matrix and Tensor Factorizations. John Wiley & Sons, 2009.
"""
function _update_nmf_α!(V::Matrix{T}, V_a::Matrix{T}, W::Matrix{T}, H::Matrix{T}, params::CParams) where T <: Real
    #### compute regularization
    #### for J_H = ||H||₁, ∇J_H = J_{K,N} with J the matrix filled with ones
    #### for J_H = 0.5*||H||₂, ∇J_H = H
    regu_H = H -> params.α_H * (params.l₁ratio_H * ones(size(H)) + (1. - params.l₁ratio_H) * H) / length(H)
    regu_W = W -> params.α_W * (params.l₁ratio_W * ones(size(W)) + (1. - params.l₁ratio_W) * W) / length(W)

    #### orthonormalization term
    orth_H = H -> params.μ_H * (H - repeat(sum(H, dims=1), size(H)[1], 1))
    orth_W = W -> params.μ_W * (W - repeat(sum(W, dims=2), 1, size(W)[2]))

    if params.α == 0
        #### update H
        H .= (H .* exp.(params.ω * transpose(_normalize(W, p=1)) * log.(V ./ V_a))) .^ (1 + params.α_sparse_H)
        V_a .= W * H

        #### update W
        W .= (W .* exp.(params.ω * log.(V ./ V_a) * transpose(_normalize(H, p=1)))) .^ (1 + params.α_sparse_W)
        V_a .= W * H
    else
        #### update H
        H .= (H .* (transpose(W) * ((V ./ (V_a .+ eps(T))) .^ params.α) - regu_H(H) + orth_H(H)) .^ (params.ω / params.α)) .^
        (1 + params.α_sparse_H)
        V_a .= W * H

        #### update W (incl. normalize columns of W)
        W .= (W .* (((V ./ (V_a .+ eps(T))) .^ params.α) * transpose(H) - regu_W(W) + orth_W(W)) .^ (params.ω / params.α)) .^
        (1 + params.α_sparse_W)
        W .= _normalize(W, p=1)
        V_a .= W * H
    end
end

"""
    nmf_α(V, W_init, H_init, params)

Run the NMF α algorithm using the set of input parameters. For the parameters, see [`CParams`](@ref)

# References
- Andrzej CICHOCKI, et. al., Nonnegative Matrix and Tensor Factorizations. John Wiley & Sons, 2009.
"""
function nmf_α(V::Matrix{T}, W_init::Matrix{T}, H_init::Matrix{T}, params::CParams) where T <: Real
    W       = W_init[:,:]
    H       = H_init[:,:]
    V_a     = W * H

    #### metrics and stopping_crit
    metrics            = NMFMetrics(V, W_init, H_init, _α_cost, params)
    stopping_crit_met  = false

    while !stopping_crit_met
        #### update V_a, W, H
        _update_nmf_α!(V, V_a, W, H, params)

        #### update metrics
        _update_metrics!(metrics, V, W, H, params)

        #### check the stopping criterion
        stopping_crit_met = _check_stopping_crit(metrics, params)

        #### print intermediate if verbose
        if params.verbose && (metrics.cur_iter % (metrics.max_iter ÷ 10) == 0)
            @printf("iter %d/%d | cost: %.5g\n", metrics.cur_iter, metrics.max_iter, metrics.absolute_cost[metrics.cur_iter])
        end
    end

    #### convergence information
    if metrics.cur_iter == metrics.max_iter && params.stopping_crit != :none
        @warn "maximum number of iterations reached"
    end
    if metrics.cur_iter < metrics.max_iter
        @info "stopping criterion met at" metrics.cur_iter
    end

    NMFResults{T}(W, H, metrics, params)
end


text:
export SmoothingSpline, SmoothingSpline2

function SmoothingSpline(θ)
    F = [2.0 -1.0; 1.0 0.0]
    Q = fill(θ[1].^2,1,1)
    G = zeros(2,1)
    H = [1.0 0.0]
    R = fill(θ[2].^2,1,1)
    Γ = zeros(1,1)
    Θ = [1.0;0.0]
    μ = [0.0;0.0]
    Σ0 = [1.0 0.0; 0.0 1.0]
    LinearStateSpaceModel(F,Q,G,H,R,Γ,Θ,μ,Σ0)
end


text:
materialFilterEntry = Entry(
    primary_icon_name="edit-find-symbolic",
    secondary_icon_name="edit-clear-symbolic",
    placeholder_text="Search..."
)

@guarded signal_connect(materialFilterEntry, "changed") do widget
    text = Gtk.bytestring(GAccessor.text(widget))
    lcText = lowercase(text)

    searchRawNames = get(config, "search_entities_raw_name", true)
    searchAttributes = get(config, "search_entities_attributes", true)
    searchTooltips = get(config, "search_entities_tooltips", false)
    needsPlacements = searchRawNames || searchAttributes || searchTooltips

    # TODO - This relies on layers list having raw names
    targetLayerIndex = Ahorn.getSelected(Ahorn.layersList)
    targetLayer = targetLayerIndex !== nothing ? Ahorn.layersList.data[targetLayerIndex][1] : nothing

    if targetLayer !== nothing
        Ahorn.persistence["material_search_$(targetLayer)"] = text
    end

    cache = Dict{String, Any}()
    placements = nothing

    tooltipKey = ""

    if needsPlacements
        if targetLayer == "entities"
            cache = Ahorn.entityPlacementsCache
            placements = Ahorn.entityPlacements
            tooltipKey = "entities"

        elseif targetLayer == "triggers"
            cache = Ahorn.triggerPlacementsCache
            placements = Ahorn.triggerPlacements
            tooltipKey = "triggers"
        end
    end

    for row in 1:length(materialList.data)
        visible = false

        materialName = materialList.data[row][1]
        visible = occursin(lcText, lowercase(materialName))

        if !visible && placements !== nothing && haskey(placements, materialName)
            ep = placements[materialName]
            name = isa(ep.func, DataType) ? string(ep.func.parameters[1]) : nothing
            target = get(cache, materialName, nothing)

            if name === nothing && target === nothing
                target = Ahorn.getCachedPlacement!(cache, placements, materialName)
                name = target.name
            end

            if name !== nothing
                tooltips = get(Ahorn.langdata, ["placements", tooltipKey, name, "tooltips"])
                visible |= searchRawNames && occursin(lcText, lowercase(name))

                if !visible && searchAttributes
                    data = target !== nothing ? target.data : tooltips
                    for attr in keys(data)
                        visible |= occursin(lcText, lowercase(string(attr)))
                    end
                end

                if !visible && searchTooltips
                    for value in values(tooltips)
                        visible |= occursin(lcText, lowercase(value))
                    end
                end
            end
        end

        materialList.visible[row] = visible
    end

    # Don't select anything automatically, notify tool to do so instead
    Gtk.GLib.@sigatom filterContainer!(materialList, nothing)
    Ahorn.notifyMaterialsFiltered(text)
end

@guarded signal_connect(materialFilterEntry, "icon-press") do widget, n, event
    Gtk.GLib.@sigatom GAccessor.text(widget, "")
end

function handleFilterKeyPressed(widget::Gtk.GObject, event::Gtk.GdkEventKey)
    if widget == materialFilterEntry &&
        if event.keyval == Gtk.GdkKeySyms.Escape
            unfocusFilterEntry!()

            return true

        elseif event.keyval == Gtk.GdkKeySyms.Return
            unfocusFilterEntry!()
            Gtk.GLib.@sigatom GAccessor.text(widget, "")

            return true
        end
    end

    return false
end

function focusFilterEntry!(args...)
    GAccessor.focus(window, materialFilterEntry)
end

function unfocusFilterEntry!(args...)
    # Can't focus the canvas, focus something we know exists
    GAccessor.focus(window, scrollableWindowMaterialList)
end


text:
# Autogenerated wrapper script for Cubature_jll for i686-w64-mingw32
export libcubature

JLLWrappers.@generate_wrapper_header("Cubature")
JLLWrappers.@declare_library_product(libcubature, "libcubature.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcubature,
        "bin\\libcubature.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import pytest
from julia import Base

:q!


text:
@testset "4.4.1.3 (d cos)^m (a+b cot)^n" begin
    (a, b, x, ) = @variables a b x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(d*cos(e+f*x))^m*(a+b*cot(e+f*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(d*cos(e+f*x))^m*(a+b*cot(e+f*x))^n*with*a^2+b^2=0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cos(e+f*x)^m*(a+b*cot(e+f*x))^n*when*a^2+b^2=0=#

    #= ::Subsubsection:: =#
    #=n>0=#

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [cos(x)^4/(im + cot(x)), x, 5, -((im*x)/16) + 1/(32*(im - cot(x))^2) + im/(8*(im - cot(x))) - im/(24*(im + cot(x))^3) + 5/(32*(im + cot(x))^2) + (3*im)/(16*(im + cot(x)))]
    @test_int [cos(x)^3/(im + cot(x)), x, 9, (-(1/5))*cos(x)^5 - (1/3)*im*sin(x)^3 + (1/5)*im*sin(x)^5]
    @test_int [cos(x)^2/(im + cot(x)), x, 5, -((im*x)/8) + im/(8*(im - cot(x))) + 1/(8*(im + cot(x))^2) + im/(4*(im + cot(x)))]
    @test_int [cos(x)^1/(im + cot(x)), x, 8, (-(1/3))*cos(x)^3 - (1/3)*im*sin(x)^3]
    @test_int [sec(x)^1/(im + cot(x)), x, 8, (-im)*atanh(sin(x)) - cos(x) + im*sin(x)]
    @test_int [sec(x)^2/(im + cot(x)), x, 3, im*x - log(sin(x)) + log(tan(x)) - im*tan(x)]
    @test_int [sec(x)^3/(im + cot(x)), x, 8, (1/2)*im*atanh(sin(x)) + sec(x) - (1/2)*im*sec(x)*tan(x)]
    @test_int [sec(x)^4/(im + cot(x)), x, 4, tan(x)^2/2 - (1/3)*im*tan(x)^3]
    @test_int [sec(x)^5/(im + cot(x)), x, 9, (1/8)*im*atanh(sin(x)) + sec(x)^3/3 + (1/8)*im*sec(x)*tan(x) - (1/4)*im*sec(x)^3*tan(x)]
    @test_int [sec(x)^6/(im + cot(x)), x, 4, tan(x)^2/2 - (1/3)*im*tan(x)^3 + tan(x)^4/4 - (1/5)*im*tan(x)^5]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(d*cos(e+f*x))^m*(a+b*cot(e+f*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cos(e+f*x)^m*(a+b*cot(e+f*x))^n=#

    #= ::Subsubsection:: =#
    #=n>0=#

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [cos(x)^4/(a + b*cot(x)), x, 8, (a*(3*a^4 - 6*a^2*b^2 - b^4)*x)/(8*(a^2 + b^2)^3) - (a^4*b*log(b*cos(x) + a*sin(x)))/(a^2 + b^2)^3 + ((4*b*(2*a^2 + b^2) + a*(5*a^2 + b^2)*cot(x))*sin(x)^2)/(8*(a^2 + b^2)^2) - ((b + a*cot(x))*sin(x)^4)/(4*(a^2 + b^2))]
    @test_int [cos(x)^3/(a + b*cot(x)), x, 10, (a^3*b*atanh((a*cos(x) - b*sin(x))/sqrt(a^2 + b^2)))/(a^2 + b^2)^(5/2) - (a^2*b*cos(x))/(a^2 + b^2)^2 - (b*cos(x)^3)/(3*(a^2 + b^2)) - (a*b^2*sin(x))/(a^2 + b^2)^2 + (a*sin(x))/(a^2 + b^2) - (a*sin(x)^3)/(3*(a^2 + b^2))]
    @test_int [cos(x)^2/(a + b*cot(x)), x, 7, (a*(a^2 - b^2)*x)/(2*(a^2 + b^2)^2) - (a^2*b*log(b*cos(x) + a*sin(x)))/(a^2 + b^2)^2 + ((b + a*cot(x))*sin(x)^2)/(2*(a^2 + b^2))]
    @test_int [cos(x)^1/(a + b*cot(x)), x, 6, (a*b*atanh((a*cos(x) - b*sin(x))/sqrt(a^2 + b^2)))/(a^2 + b^2)^(3/2) - (b*cos(x))/(a^2 + b^2) + (a*sin(x))/(a^2 + b^2)]
    @test_int [sec(x)^1/(a + b*cot(x)), x, 6, atanh(sin(x))/a + (b*atanh((a*cos(x) - b*sin(x))/sqrt(a^2 + b^2)))/(a*sqrt(a^2 + b^2))]
    @test_int [sec(x)^2/(a + b*cot(x)), x, 3, -((b*log(a + b*cot(x)))/a^2) - (b*log(tan(x)))/a^2 + tan(x)/a]
    @test_int [sec(x)^3/(a + b*cot(x)), x, 9, atanh(sin(x))/(2*a) + (b^2*atanh(sin(x)))/a^3 + (b*sqrt(a^2 + b^2)*atanh((a*cos(x) - b*sin(x))/sqrt(a^2 + b^2)))/a^3 - (b*sec(x))/a^2 + (sec(x)*tan(x))/(2*a)]
    @test_int [sec(x)^4/(a + b*cot(x)), x, 3, -((b*(a^2 + b^2)*log(a + b*cot(x)))/a^4) - (b*(a^2 + b^2)*log(tan(x)))/a^4 + ((a^2 + b^2)*tan(x))/a^3 - (b*tan(x)^2)/(2*a^2) + tan(x)^3/(3*a)]

    #= Following*hangs*Mathematica*6 & 7: =#
    @test_int [sec(x)/(1 + 2*cot(x)), x, 6, (2*atanh((cos(x) - 2*sin(x))/sqrt(5)))/sqrt(5) + atanh(sin(x))]
    end

max_stars_repo_path: test/4 Trig functions/4.4 Cotangent/4.4.1.3 (d cos)^m (a+b cot)^n.jl


text:
struct Visual
    plottype::PlotFunc
    attributes::Dict{Symbol, Any}
end
Visual(plottype=Any; kwargs...) = Visual(plottype, Dict{Symbol, Any}(kwargs))

function (v::Visual)(e::Entry)
    plottype = AbstractPlotting.plottype(v.plottype, e.plottype)
    mappings = e.mappings
    attributes = merge(e.attributes, v.attributes)
    return Entry(plottype, mappings, attributes)
end

visual(plottype=Any; kwargs...) = Layer((Visual(plottype; kwargs...),))


text:
function Base.show(
            io::IO,
            ::MIME"image/svg+xml",
            f::Figure
        )::Nothing

    opath = savefig(f, tempname(); info=false)
    write(io, read(opath))
    flush(io)
    return
end


text:
using Test
using BEAST
using CompScienceMeshes
using SparseArrays

fn = joinpath(@__DIR__, "assets/sphere5.in")
m = readmesh(fn)

Id = BEAST.Identity()
X = BEAST.raviartthomas(m)

Z1 = assemble(Id, X, X, storage_policy=Val{:densestorage})
Z2 = assemble(Id, X, X, storage_policy=Val{:bandedstorage})
Z3 = assemble(Id, X, X, storage_policy=Val{:sparsedicts})

@test Z1 isa DenseMatrix
@test Z2 isa SparseMatrixCSC
@test Z3 isa SparseMatrixCSC

@test Z1 ≈ Z2 atol=1e-8
@test Z1 ≈ Z3 atol=1e-8


text:
struct Flip end

function (::Flip)(res, x, α, β)
  res .= (-α) .* x .+ β .* res
end

function test_callable()
  @testset ExtendedTestSet "Test callable" begin
    Mv = ones(2)
    op = LinearOperator(2, 2, true, true, Flip())
    @test op * ones(2) == -ones(2)
    @test op' * ones(2) == -ones(2)
    @test transpose(op) * ones(2) == -ones(2)
    v = ones(2)
    allocs = @allocated mul!(Mv, op, v)
    @test allocs == 0
    @test Mv == -ones(2)
  end
end

test_callable()


text:
module JuliaFMI

# Dependecies
using Libdl         # For using dlopen, dlclose and so on
using LightXML      # For parsing XML files

export simulateFMU

function simulateFMU(pathToFMU::String)
    main(pathToFMU)
end

end # module

max_stars_repo_name: cuihantao/JuliaFMI

max_issues_repo_name: cuihantao/JuliaFMI

max_forks_repo_name: cuihantao/JuliaFMI


text:
"""
Struct for Lewis Integration Method

		bsProcess=CarrMadanLewisMethod(σ::num1) where {num1 <: Number}

Where:\n
		σ	=	volatility of the process.
"""
mutable struct CarrMadanLewisMethod{num <: Number, num_1 <: Integer} <: AbstractFFTMethod
    A::num
    Npow::num_1
    function CarrMadanLewisMethod(A::num, N::num_1) where {num <: Number, num_1 <: Integer}
        @assert(A > 0.0, "A must be positive")
        @assert(N > 2, "N must be greater than 2")
        @assert(N < 24, "N will cause overflow")
        return new{num, num_1}(A, N)
    end
end

export CarrMadanLewisMethod;

"""
Documentation CarrMadanLewis Method
"""
function pricer(mcProcess::FinancialMonteCarlo.BaseProcess, StrikeVec::Array{U, 1}, zero_rate::FinancialMonteCarlo.AbstractZeroRateCurve, T::Number, method::CarrMadanLewisMethod, ::FinancialMonteCarlo.BaseMode = FinancialMonteCarlo.SerialMode()) where {U <: Number}
    Npow = method.Npow
    N = 2^Npow
    S0 = mcProcess.underlying.S0
    A = method.A
    d = FinancialMonteCarlo.integral(FinancialMonteCarlo.dividend(mcProcess), T) / T
    r = FinancialMonteCarlo.integral(zero_rate.r, T) / T
    cf = FinancialFFT.CharactheristicFunction(mcProcess, T)
    correction = (r - d) * T - FinancialFFT.CharactheristicExponent(-1im, mcProcess, T)
    CharFunc(v) = cf(v) * exp(v * 1im * correction)
    dx = A / N
    x = collect(0:(N-1)) * dx
    x[1] = 1e-312
    weights_ = @. 3 + (-1)^((0:(N-1)) + 1)
    @views weights_[1] = 1
    @views weights_[N] = 1

    dk = 2 * pi / A
    b = N * dk / 2
    ks = -b .+ dk * (0:(N-1))
    integrand = @. exp(-1im * b * x) * CharFunc(x - 0.5 * 1im) / (x^2 + 0.25) * weights_ * dx / 3
    fft!(integrand)
    integral_value = @. real_mod(integrand) / pi

    spline_cub = CubicSplineInterpolation(ks, integral_value)
    prices = @. S0 - sqrt(S0 * StrikeVec) * exp(-(r - d) * T) * spline_cub(log(StrikeVec / S0))
    return prices * exp(-d * T)
end


text:
#=
Copyright 2022 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

"""
    joinVessels(b :: Blood, vessels...)

Call either `solveConjunction` or `solveBifurcation` depending on the number of given
vessels.
"""
function joinVessels(b :: Blood, vessels...)
    if length(vessels) == 2
        solveConjunction(b, vessels[1], vessels[2])

    elseif length(vessels) == 3
        solveBifurcation(vessels[1], vessels[2], vessels[3])
    end
end

"""
    newtonRaphson(vessels :: Array{Vessel,1}, J, U, k, funW, funF)

Solve the non-linear junction system by means of Newton-Raphson method.
"""
function newtonRaphson(vessels :: Array{Vessel,1}, J, U, k, funW, funF)
    W = funW(U, k)
    F = funF(vessels, U, k, W)

    nr_toll_U = 1e-5
    nr_toll_F = 1e-5

    while true
        @fastmath dU = J\(-F)
        U_new = U .+ dU

        if any(isnan(dot(F,F)))
            e = "("
            for i = 1:length(vessels)
                e *= vessels[i].label
                e *= ", "
            end
            e = e[1:end-2]*")"
            error("\nNewton-Raphson doesn't converge at $e junction!")
        end

        u_ok = 0
        f_ok = 0
        for i in 1:length(dU)
            if abs(dU[i]) <= nr_toll_U || abs(F[i]) <= nr_toll_F
                u_ok += 1
                f_ok += 1
            end
        end

        if u_ok == length(dU) || f_ok == length(dU)
            return U_new

        else
            U = U_new
            W = funW(U, k)
            F = funF(vessels, U, k, W)
        end
    end
end


text:
using AdventOfCodeSolutions
using Test
import LinearAlgebra: norm

function input(puzzle::Puzzle{2019, 3, n}) where n
    io = openInput(puzzle)
    return read(io, String)
end

function parseInput(puzzle::Puzzle{2019, 3, n}, input) where n
    map(s -> split(s, ','), split(input, '\n', keepempty=false))
end

function path(turns)
    dirMap = Dict('R' => [1 0], 'L' => [-1 0], 'U' => [0 1], 'D' => [0 -1])
    path = [[0 0]]
    foreach(map(d -> (dirMap[d[1]], parse(Int, d[2:end])), turns)) do turn
        foreach(d -> push!(path, path[end] + turn[1]), 1:turn[2])
    end
    return path[2:end]
end

manhattan(p) = norm(p, 1)
function solve(p::Puzzle{2019, 3, 1}, input)
    turns = parseInput(p, input)

    return minimum(
        manhattan,
        intersect(path(turns[1]), path(turns[2]))
    )
end

function solve(p::Puzzle{2019, 3, 2}, input)
    turns = parseInput(p, input)

    p1 = path(turns[1])
    p2 = path(turns[2])
    intersections = intersect(p1, p2)
    minimum(intersections) do point
        i1 = minimum(findall(==(point), p1))
        i2 = minimum(findall(==(point), p2))

        return i1 + i2
    end
end

@testset "Day 3" begin
    TEST_INPUT_1 = """
R75,D30,R83,U83,L12,D49,R71,U7,L72
U62,R66,U55,R34,D71,R55,D58,R83
"""
    TEST_INPUT_2 = """
R98,U47,R26,D63,R33,U87,L62,D20,R33,U53,R51
U98,R91,D20,R16,D67,R40,U7,R15,U6,R7
"""
    @test(solve(Puzzle(2019, 3, 1), TEST_INPUT_1) == 159)
    @test(solve(Puzzle(2019, 3, 1), TEST_INPUT_2) == 135)
    @test(solve(Puzzle(2019, 3, 1), input(Puzzle(2019, 3, 1))) == 1084)
    @test(solve(Puzzle(2019, 3, 2), input(Puzzle(2019, 3, 2))) == 9240)
end


text:
module DynamicWalking2018

end


text:
#=
This file is part of Pseudospectra.jl.

Julia translation
copyright 2017 Ralph Smith

Portions derived from EigTool
Copyright (c) 2002-2014, The Chancellor, Masters and Scholars
of the University of Oxford, and the EigTool Developers. All rights reserved.

SPDX-License-Identifier: BSD-3-Clause
License-Filename: LICENSES/BSD-3-Clause_Eigtool
=#

"""
    numerical_range(A, nstep=20) -> Vector{Complex}

Compute points along the numerical range of a matrix.

Note: this solves an eigensystem for each point, so may be expensive.
"""
function numerical_range(A::AbstractMatrix, thmax=20)
    rayleighquotient(B,x) = dot(x, B * x) / dot(x,x)
    # m,n = size(A)
    T = eltype(A)
    CT = (T <: Real) ? Complex{T} : T
    z = zeros(CT,2*thmax+2)
    for i = 0:thmax
        # upstream allows for interactive cancellation

        # get Hermitian part of rotated A
        th = (i/thmax)*π
        Ath = exp(th*1im)*A
        H = (1/2)*(Ath + Ath')
        F = eigen(H)
        d,X = F.values,F.vectors

        # RQ's of A correspond to eigenvalues of H w/ extreme real parts
        k = sortperm(real(d))
        z[i+1] = rayleighquotient(A,X[:,k[1]])
        z[1+i+thmax] = rayleighquotient(A,X[:,k[end]])
    end
    z[end] = z[1] # close curve for plotting
    z
end

"""
    numrange!(ps_data,nstep=20)

compute the numerical range (a.k.a. field of values) of a dense matrix
then store it in the `ps_data`.

Note: this solves an eigensystem for each point, so may be expensive.
"""
function numrange!(ps_data::PSAStruct,thmax=20)
    ps_dict = ps_data.ps_dict
    # don't recompute if no change
    if !haskey(ps_dict,:fov) || (length(ps_dict[:fov]) != 2*(thmax+1))
        # Why does upstream enforce this constraint?
        if !haskey(ps_dict,:schur_mtx)
            throw(ArgumentError("only implemented for Schur-factored matrices"))
        end
        z = numerical_range(ps_dict[:schur_mtx],thmax)
        ps_dict[:fov] = z
    else
        z = ps_dict[:fov]
    end

    zoom = ps_data.zoom_list[ps_data.zoom_pos]
    ax = zoom.ax
    if mapreduce(w -> ((real(w) < ax[1]) | (real(w) > ax[2]) |
                       (imag(w) < ax[3]) | (imag(w) > ax[4])),
                 &, z, init=true)
        warn("The boundary of the numerical range is not visible on the "
             * "current axes; expand axis limits to see it. "
             * "A bounding box is $(extrema(real(z))) $(extrema(imag(z))).")
    end
end

"""
    numerical_abscissa(A)

Compute the numerical abscissa of a matrix `A`, `ω(A)`.

Uses `eigvals()`. `ω(A)` provides bounds and limiting behavior for
`norm(expm(t*A))`.
"""
function numerical_abscissa(A::AbstractMatrix)
    (1/2) * maximum(eigvals(A+A'))
end