text:
import StatsBase

using LinearAlgebra: norm, normalize, dot
using Statistics

const MINFACT = 0.1
const HITCORR = 0.9
const MINCORR = 0.25 # reasonable default
const MAXFUZZ = 1.8
const MINFUZZ = 1.0
const TARGETNGUESS = 200
const Guess{N} = NTuple{N,Tuple{Int64,Float64}} where N

"""
    overlay!(signal, chunks, positions)

Overlay a series of chunks on top of a signal at the
specified positions.
"""
function overlay!(signal, chunks, positions)
    for (i,c) in enumerate(chunks)
        p = positions[i]
        signal[p:(p+length(c)-1)] .+= c
    end
    signal
end

function alignments(signal, chunk, start_pos; tol=250, fuzziness=1.5,
                    minsig=MINFACT*norm(signal), mincorr=MINCORR)
    l = length(chunk)
    section = max(1, -tol+start_pos):min(length(signal), tol+l-1+start_pos)
    s = @view signal[section]
    # reject weak signals, important if signal only contains noise
    c = normalize(chunk)
    # return s,c
    corr = zeros(length(section)-l+1)
    lc = length(corr)
    lratio = l/length(signal)
    for i in eachindex(corr)
        sig = @view s[i:(i+l-1)]
        xc = dot(c,sig)
        corr[i] = xc < minsig*lratio ? 0.001 :  xc/norm(sig)
    end
    M = maximum(corr)
    if M < mincorr
        return Tuple{Int64,Float64}[]
    end
    # println(M)
    ps = section.start - 1
    let lc=lc, corr=corr, ps=ps, M=M
        Tuple{Int64,Float64}[ (i+ps, corr[i]) for i=eachindex(corr) if
                              # corr[i] > mincorr &&
                              (i==1 || corr[i] > corr[i-1]) &&
                              (i==lc || corr[i] > corr[i+1]) &&
                              M/corr[i] < fuzziness ]
    end
end

alignments(s::Spectrum, args...; kw...) = alignments(s[:],args...;kw...)
alignments(s::Spectrum, l::Spectrum, args...; kw...) =
    Array{Tuple{Int64,Float64},1}[ alignments(s[:], l[r], r.start, args...; kw...) for r in intrng_indices(l) ]

guesses(s::Spectrum, l::Spectrum; kw...) = vec(collect(Base.product(alignments(s, l; kw...)...)))
function guesses_adaptive(s::Spectrum, l::Spectrum;
                          ntarget=TARGETNGUESS, δ=ntarget/4, kw...)
    fuzz = (MAXFUZZ + MINFUZZ) / 2
    @binary_opt ( *(length.(alignments(s, l; fuzziness=fuzz, kw...))...) - ntarget ) fuzz MINFUZZ MAXFUZZ δ
    guesses(s, l; fuzziness=fuzz, kw...)
end

function projection(signal, chunk, start_pos)
    l = length(chunk)
    s = @view signal[start_pos:start_pos+l-1]
    dot(s, chunk)/norm(chunk)^2
end

positions(guess) = Int64[ g[1] for g in guess ]
fitness(guess) = Float64[ g[2] for g in guess ]

function fitness_weights(guess, λ=10.0)
    f = fitness(guess)
    F = maximum(f)
    exp.(λ/F*(f.-F))
end

function projection_weights(projs, fitness_weights=ones(length(projs)), η=2.0)
    # p = projections(s, l, positions(guess))
    fw = StatsBase.weights(fitness_weights)
    σ = std(projs, fw; corrected=false)/η
    m = mean(projs, fw)
    exp.(-(projs.-m).^2/2σ^2)
end

function projections(s::Spectrum, l::Spectrum, positions)
    sig = s[:]
    Float64[ projection(sig, r, p) for (p,r) in zip(positions, intrng_data(l)) ]
end

# function projection(s::Spectrum, l::Spectrum, positions, weights=ones(positions))
#     projs = projections(s, l, positions)
#     mean(projs, StatsBase.weights(weights))
# end

function projection(s::Spectrum, l::Spectrum, guess::NTuple{N,Tuple{Int64,Float64}}) where N
    projs = projections(s, l, positions(guess))
    if length(projs) == 1
        return first(projs)
    end
    fw = fitness_weights(guess)
    # projection(s, l, positions(guess), f(s, l, guess))
    mean(projs, StatsBase.weights(projection_weights(projs, fw)))
    # projection(s, l, positions(guess), Float64[f(ff) for ff in fitness(guess)])
end

function projection_score(s::Spectrum, l::Spectrum, guess::NTuple{N,Tuple{Int64,Float64}}) where N
    projs = projections(s, l, positions(guess))
    mean(projs)/std(projs, StatsBase.weights(fitness(guess)); corrected = false)
end

function fit_score(s::Spectrum, l::Spectrum, guess::NTuple{N,Tuple{Int64,Float64}}) where N
    if N == 0
        0.0
    else
        *((g[2] for g in guess)...)^(1/N)
    end
end

# overall score, not comparable between references
score(s, l, g) = fit_score(s, l, g) * projection_score(s, l, g)

synthesize(l::Spectrum, positions) =
    overlay!(zeros(length(l)), intrng_data(l), positions)
synthesize(l::Spectrum, guess::Guess{N}) where N =
    synthesize(l, positions(guess))

# function score(s::Spectrum, l::Spectrum, guess::AbstractArray{Tuple{Int,Float64}})
#     *((g[2] for g in guess)...)
# end

function aligned_signals(s::Spectrum, l::Spectrum; sloppiness=0, kw...)
    positions = alignments(s, l; kw...)
    matched = filter(p->!isempty(p), positions)
    if length(matched) == 0 ||  length(positions) - length(matched) > sloppiness
        return []
    end
    chunks = [l[r] for (i,r) in enumerate(intrng_indices(l)) if !isempty(positions[i])]
    (overlay!(zeros(length(s[:])), chunks, comb) for
        comb in Base.product(matched...))
end

struct DecompositionResult
    coefficients :: Vector{Float64}
    refnums :: Vector{Int}
    signal :: Vector{Float64}
    matrix :: Matrix{Float64}
end

function lsq_analyze(s::Spectrum, lib::AbstractArray{Spectrum}, found; kw...)
    ll = length(lib)
    gs = Array{Any}(undef, ll)
    fit_scores = Array{Array{Float64,1}}(undef, ll)
    scores = Array{Array{Float64,1}}(undef, ll)
    @Threads.threads for r=1:ll
        gs[r] = guesses_adaptive(s, lib[r]; kw...)
        fit_scores[r] = isempty(gs[r]) || r ∈ found ? [0.0] : fit_score.(Ref(s), Ref(lib[r]), gs[r])
        scores[r] = isempty(gs[r]) || r ∈ found ? [0.0] : score.(Ref(s), Ref(lib[r]), gs[r])
    end
    # find best guess per reference based on overall score
    bestinds = argmax.(scores)
    # find best reference based on fit_score
    max_score,max_ref = findmax(getindex.(fit_scores, bestinds))
    if max_score == 0.0
        return (s, 0, (), 0.0, Float64[])
    end
    # pss = score.(s, lib[maxind], gs[maxind])
    # _,maxguessind = findmax(pss)
    maxguess = gs[max_ref][bestinds[max_ref]]
    ss = deepcopy(s)
    l = synthesize(lib[max_ref], maxguess)
    p = projection(s, lib[max_ref], maxguess)
    ss[:] .-= p.*l
    ss, max_ref, maxguess, p, l
end

function lsq_analyze(s::Spectrum, lib::AbstractArray{Spectrum};
                     callback = refs -> println("Found: #$refs"), kw...)
    found = Int64[]
    coeffs = Float64[]
    vecs = Array{Float64,1}[]
    ss = copy(s[:])
    ms = MINFACT * norm(ss)
    while true
        s, m, g, p, v = lsq_analyze(s, lib, found; minsig=ms, kw...)
        if m == 0
            break
        end
        push!(found, m)
        push!(vecs, v)
        push!(coeffs, p)
        callback([m])
    end
    DecompositionResult(coeffs, found, ss, hcat(vecs...))
end

function decompose(d::DecompositionResult)
    if isempty(d.coefficients)
        return ([],zeros(length(d.signal)), d.signal)
    end
    recon = d.matrix * d.coefficients
    residue = d.signal .- recon
    components = d.coefficients'.*d.matrix
    (components, recon, residue)
end


text:
#=
This file is taken from: https://github.com/rikhuijzer/PlutoStaticHTML.jl/blob/main/src/PlutoStaticHTML.jl
in Oct 2021. The project has grown since with parallel builds/caching like features.
Kindly use PlutoStaticHTML.jl instead.
=#

using Pluto:
  Cell,
  CellOutput,
  Notebook,
  ServerSession,
  SessionActions,
  generate_html,
  load_notebook_nobackup,
  update_run!,
  notebook_to_js

"""
    IMAGEMIME
Union of MIME image types.
Based on Pluto.PlutoRunner.imagemimes.
"""
const IMAGEMIME = Union{
  MIME"image/svg+xml",
  MIME"image/png",
  MIME"image/jpg",
  MIME"image/jpeg",
  MIME"image/bmp",
  MIME"image/gif"
}

"""
    _escape_html(s::AbstractString)
Escape HTML.
Useful for showing HTML inside code blocks, see
https://github.com/rikhuijzer/PlutoStaticHTML.jl/issues/9.
"""
function _escape_html(s::AbstractString)
  s = replace(s, "
" => """
""")
  s = replace(s, '<' => "<")
  s = replace(s, '>' => ">")
  return s
end

function code_block(code; class = "language-julia")
  if code == ""
    return ""
  end
  code = _escape_html(code)
  return """
$code
"""
end

function output_block(s; class = "code-output")
  if s == ""
    return ""
  end
  return """
$s
"""
end

function _code2html(code::AbstractString, class, hide_md_code)
  if hide_md_code && startswith(code, "md\"")
    return ""
  end
  if contains(code, "# hideall")
    return ""
  end
  sep = '\n'
  lines = split(code, sep)
  filter!(!endswith("# hide"), lines)
  code = join(lines, sep)
  return code_block(code; class)
end

function _output2html(body, T::IMAGEMIME, class)
  encoded = base64encode(body)
  uri = "data:$T;base64,$encoded"
  return """"""
end

function _output2html(body, ::MIME"application/vnd.pluto.stacktracetobject", class)
  return error(body)
end

function _tr_wrap(elements::Vector)
  joined = join(elements, '\n')
  return "
\n$joined\n
"
end
_tr_wrap(::Array{String,0}) = "
\n...	\n
"

function _output2html(body::Dict{Symbol,Any}, ::MIME"application/vnd.pluto.table+object", class)
  rows = body[:rows]
  nms = body[:schema][:names]
  headers = _tr_wrap(["$colname	" for colname in nms])
  contents = map(rows) do row
    # Drop index.
    row = row[2:end]
    # Unpack the type and throw away mime info.
    elements = try
      first.(only(row))
    catch
      first.(first.(row))
    end
    elements = ["$e	" for e in elements]
    return _tr_wrap(elements)
  end
  content = join(contents, '\n')
  return """

      $headers
      $content

      """
end

abstract type Struct end

function symbol2type(s::Symbol)
  if s == :Tuple
    return Tuple
  elseif s == :Array
    return Array
  elseif s == :struct
    return Struct
  else
    @warn "Missing type: $s"
    return Missing
  end
end

"""
    _clean_tree(parent, element::Tuple{Any, Tuple{String, MIME}}, T)
Drop metadata.
For example, `(1, ("\"text\"", MIME type text/plain))` becomes "text".
"""
function _clean_tree(parent, element::Tuple{Any,Tuple{String,MIME}}, T)
  return first(last(element))
end

function _clean_tree(parent, element::Tuple{Any,Any}, T)
  embedded = first(last(element))
  if embedded isa String
    return embedded
  end
  struct_name = embedded[:prefix]
  elements = embedded[:elements]
  subelements = [_clean_tree(parent, e, Nothing) for e in elements]
  joined = join(subelements, ", ")
  return struct_name * '(' * joined * ')'
end

function _clean_tree(parent, elements::Tuple{Any,Tuple}, T)
  body = first(last(elements))
  T = symbol2type(body[:type])
  return _clean_tree(body, body[:elements], T)
end

function _clean_tree(parent, elements::AbstractVector, T::Type{Tuple})
  cleaned = [_clean_tree(parent, e, Nothing) for e in elements]
  joined = join(cleaned, ", ")
  return "($joined)"
end

function _clean_tree(parent, elements::AbstractVector, T::Type{Array})
  cleaned = [_clean_tree(parent, e, Nothing) for e in elements]
  joined = join(cleaned, ", ")
  return "[$joined]"
end

function _clean_tree(parent, elements::AbstractVector, T::Type{Struct})
  cleaned = [_clean_tree(parent, e, Nothing) for e in elements]
  joined = join(cleaned, ", ")
  return parent[:prefix] * '(' * joined * ')'
end

# Fallback. This shouldn't happen. Convert to string to avoid failure.
function _clean_tree(parent, elements, T)
  @warn "Couldn't convert $parent"
  return string(elements)::String
end

function _output2html(body::Dict{Symbol,Any}, ::MIME"application/vnd.pluto.tree+object", class)
  T = symbol2type(body[:type])
  cleaned = _clean_tree(body, body[:elements], T)
  return output_block(cleaned; class)
end

_output2html(body, ::MIME"text/plain", class) = output_block(body)
_output2html(body, ::MIME"text/html", class) = body
_output2html(body, T::MIME, class) = error("Unknown type: $T")

function _cell2html(cell::Cell, code_class, output_class, hide_md_code)
  code = _code2html(cell.code, code_class, hide_md_code)
  output = _output2html(cell.output.body, cell.output.mime, output_class)
  return """
      $code
      $output
      """
end

function run_notebook!(notebook, session)
  cells = [last(e) for e in notebook.cells_dict]
  update_run!(session, notebook, cells)
  return nothing
end

"""
    notebook2html(
        notebook::Notebook;
        session=ServerSession(),
        code_class="language-julia",
        output_class="code-output",
        hide_md_code=true
    )

Run the `notebook` and return the code and output as HTML.
"""
function notebook2html(
  notebook::Notebook;
  session = ServerSession(),
  code_class = "language-julia",
  output_class = "code-output",
  hide_md_code = true
)

  cells = [last(e) for e in notebook.cells_dict]
  update_run!(session, notebook, cells)
  order = notebook.cell_order
  outputs = map(order) do cell_uuid
    cell = notebook.cells_dict[cell_uuid]
    _cell2html(cell, code_class, output_class, hide_md_code)
  end
  return join(outputs, '\n')
end

function notebook2htmlSessionActions(input)
  session = ServerSession()
  notebook = SessionActions.open(session, input; run_async = false)

  println("Type of: ", typeof(notebook.cells))
  println("Cell order", notebook.cell_order)
end

"""
    notebook2html2(notebook::Notebook)

Run the Pluto notebook and return the code and output as HTML.
"""
function notebook2html2(
  notebook::Notebook;
  code_class = "language-julia",
  output_class = "code-output",
  hide_md_code = true
)

  cells = [last(e) for e in notebook.cells_dict]
  #update_run!(session, notebook, cells)
  order = notebook.cell_order
  outputs = map(order) do cell_uuid
    cell = notebook.cells_dict[cell_uuid]
    _cell2html(cell, code_class, output_class, hide_md_code)
  end
  return join(outputs, '\n')
end

function navbar(navbar::String; class = "nav-box")

end

function upperwrap()
  return """<!DOCTYPE html>

    <meta charset="UTF-8">
    <meta http-equiv="X-UA-Compatible" content="IE=edge">
    <meta name="viewport" content="width=device-width, initial-scale=1.0">

    <meta property="og:title" content="Genie Framework - Highly Productive Web Development with Julia"/>
    <meta property="og:description" content="Documentation for Genie Framework, the highly productive web development framework for Julia. It includes all you need to quickly build production ready web applications with Julia Lang."/>
    <meta property="og:image" itemprop="image" content="/assets/img/enviroment-website.svg">
    <link rel="icon" href="/assets/img/logo-genieframework.svg">
    <title>Documentation</title>

    <link href="/assets/external/bootstrap-5.1.3-dist/css/bootstrap.min.css"  rel="stylesheet"></link>

    <link href="/assets/external/fontawesome-free-6.0.0-beta3-web/css/fontawesome.min.css" rel="stylesheet"></link>
    <link href="/assets/external/fontawesome-free-6.0.0-beta3-web/css/solid.min.css" rel="stylesheet"></link>
    <link href="/assets/external/fontawesome-free-6.0.0-beta3-web/css/brands.min.css" rel="stylesheet"></link>

    <link rel="stylesheet" href="/assets/styles-docs.css">

    <link href="/assets/external/highlight/styles/xcode.min.css" rel="stylesheet"></link>

              GenieFramework

          <button class="navbar-toggler" type="button" data-bs-toggle="collapse" data-bs-target="#navbarNav" aria-controls="navbarNav" aria-expanded="false" aria-label="Toggle navigation">

          </button>

                Join our Community

Discord

Github

Gitter

            Tutorials
             »
            Welcome

Welcome

        <form class="docs-search" action="/docs/search.html"><input class="docs-search-query2 search" id="documenter-search-query2" name="q" type="text" placeholder="Search docs"/></form>

"""
end

#
function wrap_output(my_output; class = "container-fluid text")
  if my_output == ""
    return ""
  end
  return """
$my_output
"""
end

#

function lowerwrap()
  return """

        <form class="docs-search" action="/docs/search.html"><input class="docs-search-query" id="documenter-search-query" name="q" type="text" placeholder="Search docs"/></form>

                Genie

                  Tutorials

                        Welcome to Genie

                        Installing Genie

                        Getting Started

                        Creating a Web service

                        Developing MVC web applications

                        Handling URI/query params

                        Working with forms and POST payloads

                        Using JSON payloads

                        Uploading files

                        Adding your libraries into Genie

                        Loading and starting Genie apps

                        Managing Genie app's dependencies

                        Advanced routing

                        Auto-loading configuration code with initializers

                        The secrets file

                        Auto-loading user libraries

                        Using Genie with Docker

                        Working with Web Sockets

                        Force compiling route handlers

                        Deploying to Heroku with Buildpacks

                        Deploying to Heroku with Nginx

                  Guides

                        Working with Genie apps

                        Using Genie Plugins

                        Developing an API backend

                        Using Genie in an interactive environment

                  APIs

                          App

                          AppServer

                          Assets

                          Cache

                          Commands

                          Configuration

                          Cookies

                        Deploy Docker

                        Deploy Heroku

                          Encryption

                          Exceptions

                          FileTemplates

                          Flash

                          Generator

                          Genie

                          Headers

                          HttpUtils

                          Inflector

                          Input

                          Plugins

                        Renderer

                          HTML Renderer

                        JS Renderer

                          JSON Renderer

                          Requests

                          Responses

                          Router

                          Sessions

                          Toolbox

                          Util

                          WebChannels

GenieFramework

            ©2022 | All rights reserved

Genie

Stipple

SearchLight

Docs

Blog

App Gallery

Contact

Genie Cloud

Discover the Genie community through our social media channels 📣

"""
end

# Final page
function body(notebook::Notebook)
  return join(upperwrap() * wrap_output(notebook2html2(notebook::Notebook)) * lowerwrap())
end

export notebook2html2, run_notebook!, body

ext: jl


text:
using CLIMA.PlanetParameters: Omega
export Gravity, RayleighSponge, Subsidence, GeostrophicForcing, Coriolis

# kept for compatibility
# can be removed if no functions are using this
function atmos_source!(f::Function, m::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
  f(source, state, aux, t)
end
function atmos_source!(::Nothing, m::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
end
# sources are applied additively
function atmos_source!(stuple::Tuple, m::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
  map(s -> atmos_source!(s, m, source, state, aux, t), stuple)
end

abstract type Source
end

struct Gravity <: Source
end
function atmos_source!(::Gravity, m::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
  source.ρu -= state.ρ * aux.orientation.∇Φ
end

struct Subsidence <: Source
end
function atmos_source!(::Subsidence, m::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
  source.ρu -= state.ρ * m.radiation.D_subsidence
end

struct Coriolis <: Source
end
function atmos_source!(::Coriolis, m::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
  # note: this assumes a SphericalOrientation
  source.ρu -= SVector(0, 0, 2*Omega) × state.ρu
end

struct GeostrophicForcing{FT} <: Source
  f_coriolis::FT
  u_geostrophic::FT
  v_geostrophic::FT
end
function atmos_source!(s::GeostrophicForcing, m::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
  u = state.ρu / state.ρ
  u_geo = SVector(s.u_geostrophic, s.v_geostrophic, 0)
  source.ρu -= state.ρ * s.f_coriolis * (u - u_geo)
end

"""
  RayleighSponge{FT} <: Source
Rayleigh Damping (Linear Relaxation) for top wall momentum components
Assumes laterally periodic boundary conditions for LES flows. Momentum components
are relaxed to reference values (zero velocities) at the top boundary.
"""
struct RayleighSponge{FT} <: Source
  "Domain maximum height [m]"
  zmax::FT
  "Vertical extent at with sponge starts [m]"
  zsponge::FT
  "Sponge Strength 0 ⩽ c_sponge ⩽ 1"
  c_sponge::FT
end
function atmos_source!(s::RayleighSponge, m::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
  FT = eltype(state)
  z = aux.orientation.Φ / grav
  coeff = FT(0)
  if z >= s.zsponge
    coeff_top = s.c_sponge * (sinpi(FT(1/2)*(z - s.zsponge)/(s.zmax-s.zsponge)))^FT(4)
    coeff = min(coeff_top, 1.0)
  end
  source.ρu -= state.ρu * coeff
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using Test, Gurobi

@testset "GRB Parameters" begin
    env = Gurobi.Env()
    model = Gurobi.Model(env, "")

    @test_throws Gurobi.GurobiError Gurobi.get_str_param(model, "Foobar")
    @test_throws Gurobi.GurobiError Gurobi.set_str_param!(model, "Foobar", "")
    @test_throws Gurobi.GurobiError Gurobi.get_int_param(model, "Foobar")
    @test_throws Gurobi.GurobiError Gurobi.set_int_param!(model, "Foobar", 0)
    @test_throws Gurobi.GurobiError Gurobi.get_dbl_param(model, "Foobar")
    @test_throws Gurobi.GurobiError Gurobi.set_dbl_param!(model, "Foobar", 0.0)
end


text:
@testset "Canonisation and Compression" begin

D = 10
Dcut = 5

d = 2
sites = 5

T = Float32

ψ = CUDA.randn(CuMPS{T}, sites, D, d)
ϕ = CUDA.randn(CuMPS{T}, sites, D, d)
χ = CUDA.randn(CuMPS{T}, sites, D, d)
Φ = CUDA.randn(CuMPS{T}, sites, D, d)

@testset "Canonisation (left)" begin
    canonise!(ψ, :left)

    is_left_normalized = true
    for i ∈ 1:length(ψ)
        A = ψ[i]
        DD = size(A, 3)

        @cutensor Id[x, y] := conj(A[α, σ, x]) * A[α, σ, y] order = (α, σ)
        is_left_normalized *= norm(Id - cu(I(DD))) < 1e-5
    end

    @test is_left_normalized
    @test dot(ψ, ψ) ≈ 1
end

@testset "Canonisation (right)" begin
    canonise!(ϕ, :right)

    is_right_normalized = true
    for i ∈ 1:length(ϕ)
        B = ϕ[i]
        DD = size(B, 1)

        @cutensor Id[x, y] := B[x, σ, α] * conj(B[y, σ, α]) order = (α, σ)
        is_right_normalized *= norm(Id - cu(I(DD))) < 1e-5
    end

    @test is_right_normalized
    @test dot(ϕ, ϕ) ≈ 1
end

@testset "Cauchy-Schwarz inequality (after truncation)" begin
    @test abs(dot(ϕ, ψ)) <= norm(ϕ) * norm(ψ)
end

@testset "Canonisation (both)" begin
    canonise!(χ)
    @test dot(χ, χ) ≈ 1
end

@testset "Truncation (SVD, right)" begin
    truncate!(ψ, :right, Dcut)
    @test dot(ψ, ψ) ≈ 1
end

@testset "Truncation (SVD, left)" begin
    truncate!(ψ, :left, Dcut)
    @test dot(ψ, ψ) ≈ 1
end

@testset "Variational compression" begin
    Dcut = 5
    tol = 1E-4
    max_sweeps = 5

    canonise!(Φ, :right)
    @test dot(Φ, Φ) ≈ 1

    Ψ = compress(Φ, Dcut, tol, max_sweeps)

    @test dot(Ψ, Ψ) ≈ 1

    overlap = dot(Ψ, Φ)
    dist1 = 2 - 2 * real(overlap)
    dist2 = norm(Ψ)^2 + norm(Φ)^2 - 2 * real(overlap)

    @test abs(dist1 - dist2) < 1e-5
end

end


text:
###############################
# Weighted Geometric Chung Lu #
###############################
"""
wGCL(edges::Array{Int,2}, weights::Vector{Float64}, comm::Array{Int,2},
            embed::Array{Float64,2}, distances::Vector{Float64}, verbose::Bool = false)

Calculates Weighted Geometric Chung-Lu model and divergence score for graph and embedding.

**Arguments**
* `edges::Array{Int,2}` array with edges definition (two whitespace separated vertices ids)
* `eweights::Vector{Float64}` edges weights
* `comm::Array{Int,2}` assignment of vertices to communities
* `embed::Array{Float64,2}` array with vertices embeddings
* `distances::Vector{Float64}` distances between vertices
* `vweights::Vector{Float64}` landmarks total weights - used only with landmarks approximation
* `init_vweights::Vector{Float64}` vector with original (full) vertices weights - used only with landmarks approximation
* `v_to_l::Vector{Int}` mapping from vertices to landmarks (landmarks membership) - used only with landmarks approximation
* `init_edges::Array{Int,2}` array with original (full) graph edges - used only with landmarks approximation
* `init_eweights::Vector{Float64}` vector with original (full) edges weights - used only with landmarks approximation
* `init_embed::Matrix{Float64}` array with original embedding for full graph - used only with landmarks approximation
* `split::Bool` indicator for splitting JS divergence score (global score)
* `seed::Int` RNG seed for local measure score
* `auc_samples::Int` no. samples for local measure score
* `verbose::Bool` verbose switch, if true prints additional processing information
"""
function wGCL(edges::Array{Int,2}, eweights::Vector{Float64}, comm::Matrix{Int},
            embed::Matrix{Float64}, distances::Vector{Float64}, vweights::Vector{Float64},
            init_vweights::Vector{Float64}, v_to_l::Vector{Int}, init_edges::Array{Int,2},
            init_eweights::Vector{Float64}, init_embed::Matrix{Float64},
            split::Bool, seed::Int=-1, auc_samples::Int=10000, verbose::Bool=false)

    # Default parameters values
    epsilon = 0.25 # learning rate in Chung Lu model
    delta = 0.001 # desired precision of degree estimation in Chung Lu model
    AlphaMax = 10.0 # upper bound of alpha search
    AlphaStep = 0.25 # step in alpha search
    alpha_div_counter = alpha_auc_counter = 5 # early stopping threshold (iterations of alpha search without improvement)
    skip_div = skip_auc = false # skippig flags based on early stopping counters

    no_vertices = maximum(edges)
    no_edges = size(edges,1)
    verbose && println("auc_samples: $auc_samples")
    landmarks = !isempty(v_to_l)

    verbose && println("Graph has $no_vertices vertices and $no_edges edges")
    landmarks && verbose && println("Original graph has $(maximum(init_edges)) vertices and $(size(init_edges,1)) edges")

    # Read communities
    @assert size(comm,1) == no_vertices "No. communities not matching no. vertices"
    n_parts = maximum(comm)
    verbose && println("Graph has $n_parts communities")

    # Compute C-vector
    vect_len = Int(n_parts*(n_parts+1)/2)
    vect_C = zeros(Float64, vect_len)
    vect_B = zeros(Float64, vect_len)

    for i in 1:no_edges
        a = eweights[i]
        j,k = extrema([comm[edges[i,1]],comm[edges[i,2]]])
        vect_C[idx(n_parts, j, k)] += a
    end

    # Intra-community indicator
    vect_I = falses(vect_len)
    j = 1
    for i in 1:n_parts
        vect_I[j] = true
        j += (n_parts-i+1)
    end
    best_div = best_div_ext = best_div_int = best_auc_err = best_auc = typemax(Float64)
    best_alpha = best_alpha_auc = -1.0
    dim = size(embed,2)
    verbose && println("Embedding has $dim dimensions")
    # Loop over alpha's

    # Compute Euclidean distance vector D[] given the embedding
    p_len = no_vertices*(no_vertices + 1) ÷ 2
    D = zeros(Float64,p_len)
    @assert length(distances) == no_vertices "Distances vector length is not equal to no. vertices"
    for i in 1:no_vertices
        for j in i:no_vertices
            l = idx(no_vertices, i, j)
            if i == j
                D[l] = distances[i]
            else
                D[l] = dist(i,j,embed)
            end
        end
    end
    lo,hi = extrema(D)
    D = (D.-lo)./(hi-lo) # normalize to [0,1]

    adj_no_vertices = no_vertices
    adj_edges = edges
    adj_eweights = eweights
    if landmarks
        adj_no_vertices = length(init_vweights)
        adj_edges = init_edges
        adj_eweights = init_eweights
    end

    if landmarks
        adj_p_len = adj_no_vertices*(adj_no_vertices+1) ÷ 2
        full_graph_D = zeros(Float64, adj_p_len)
        for i in 1:adj_no_vertices
            for j in (i+1):adj_no_vertices
                l = idx(adj_no_vertices, i, j)
                full_graph_D[l] = dist(i,j,init_embed)
            end
        end
        lo,hi = extrema(full_graph_D)
        full_graph_D = (full_graph_D.-lo)./(hi-lo) # normalize to [0,1]
    end

    # Chung Lu model weights to be tuned
    T = ones(no_vertices)

    # Generate arrays of edges and non-edges
    NE = Tuple{Int64,Int64}[]
    for i in 1:adj_no_vertices
        for j in i:adj_no_vertices
            if i!=j
                push!(NE,(i,j))
            end
        end
    end

    ## Tuples of edges
    E = Tuple{Int64,Int64,Float64}[]
    for (i,e) in enumerate(eachrow(adj_edges))
        push!(E,(minimum(e),maximum(e),adj_eweights[i]))
    end

    ## Tuples of non-edges
    NE = collect(setdiff(Set(NE),Set([e[1:2] for e in E])));

    for alpha in AlphaStep:AlphaStep:(AlphaMax+delta)
        write(stderr, ".")
        # Apply kernel (g(dist))
        GD = zeros(Float64,p_len)
        for i in 1:no_vertices
            for j in i:no_vertices
                k = idx(no_vertices,i,j)
                GD[k] = (1-D[k])^alpha # transform w.r.t. alpha
            end
        end
        # Learn GCL model numerically
        diff = 1.0
        while diff > delta # stopping criterion
            S = zeros(no_vertices)
            for i in 1:no_vertices
                for j in i:no_vertices
                    tmp = T[i]*T[j]*GD[idx(no_vertices, i, j)]
                    S[i] += tmp
                    if i!=j S[j] += tmp end
                end
            end
            f = 0.0
            for i in 1:no_vertices
                move = epsilon*T[i]*(vweights[i]/S[i]-1.0)
                T[i] += move
                f = max(f,abs(vweights[i]-S[i])) # convergence w.r.t. degrees
            end
            diff = f
            # verbose && println("diff = $diff")
        end
        # Compute probas P[]
        P = zeros(Float64,p_len)
        for i in 1:no_vertices
            for j in i:no_vertices
                k = idx(no_vertices,i,j)
                P[k] = T[i]*T[j]*GD[k]
            end
        end

        if !skip_auc
            pos = zeros(auc_samples)
            neg = zeros(auc_samples)
            auc_weights = zeros(auc_samples)
            if landmarks
                ## Random positive cases
                seed != -1 && Random.seed!(seed)
                for (ind, edge) in enumerate(sample(E,auc_samples,replace=true))
                    i, j, w = edge
                    adj_Ti = T[v_to_l[i]]*init_vweights[i]/vweights[v_to_l[i]]
                    adj_Tj = T[v_to_l[j]]*init_vweights[j]/vweights[v_to_l[j]]
                    pos[ind] = adj_Ti*adj_Tj*((1-full_graph_D[idx(adj_no_vertices,i,j)])^alpha)
                    auc_weights[ind] = w
                end
                ## Random negative cases
                seed != -1 && Random.seed!(seed)
                for (ind, edge) in enumerate(sample(NE,auc_samples,replace=true))
                    i, j = edge
                    adj_T_i = T[v_to_l[i]]*init_vweights[i]/vweights[v_to_l[i]]
                    adj_T_j = T[v_to_l[j]]*init_vweights[j]/vweights[v_to_l[j]]
                    neg[ind] = adj_T_i*adj_T_j*((1-full_graph_D[idx(adj_no_vertices,i,j)])^alpha)
                end
            else
                ## Random positive cases
                seed != -1 && Random.seed!(seed)
                for (ind, edge) in enumerate(sample(E,auc_samples,replace=true))
                    i, j, w = edge
                    pos[ind] = P[idx(no_vertices,i,j)]
                    auc_weights[ind] = w
                end
                ## Random negative cases
                seed != -1 && Random.seed!(seed)
                neg = [P[idx(no_vertices,i,j)] for (i,j) in sample(NE,auc_samples,replace=true)]
            end
            # 1 - AUC
            auc = 1 - sum((pos .> neg).*auc_weights)/sum(auc_weights)

            if auc < best_auc
                best_auc = auc
                best_auc_err = 1.96 * sqrt(auc*(1-auc)/auc_samples) ## error from 95% CI
                best_alpha_auc = alpha
                alpha_auc_counter = 5
            else
                alpha_auc_counter -= 1
                skip_auc = alpha_auc_counter == 0
            end
        end

        if !skip_div
            # Compute B-vector given P[] and comm[]
            vect_B = zeros(vect_len)
            for i in 1:no_vertices
                for j in i:no_vertices
                    k,l = extrema([comm[i],comm[j]])
                    vect_B[idx(n_parts, k, l)] += P[idx(no_vertices,i,j)]
                end
            end
            if !split
                f = JS(vect_C, vect_B, Bool[], true)
            else
                x = JS(vect_C, vect_B, vect_I, true)
                y = JS(vect_C, vect_B, vect_I, false)
                f = (x+y)/2.0
            end
            if f < best_div
                best_div = f
                best_alpha = alpha
                best_div_ext = !split ? 0.0 : x
                best_div_int = !split ? 0.0 : y
                alpha_div_counter = 5
            else
                alpha_div_counter -= 1
                skip_div = alpha_div_counter == 0
            end
        end
        skip_div && skip_auc && break
    end
    write(stderr, "\n")
    return [best_alpha, best_div, best_div_ext, best_div_int, best_alpha_auc, best_auc, best_auc_err]
end

"""
wGCL_directed(edges::Array{Int,2}, weights::Vector{Float64}, comm::Array{Int,2},
            embed::Array{Float64,2}, distances::Vector{Float64}, verbose::Bool = false)

Calculates directed Weighted Geometric Chung-Lu model and divergence score for graph and embedding.

**Arguments**
* `edges::Array{Int,2}` array with edges definition (two whitespace separated vertices ids)
* `eweights::Vector{Float64}` edges weights
* `comm::Array{Int,2}` assignment of vertices to communities
* `embed::Array{Float64,2}` array with vertices embeddings
* `distances::Vector{Float64}` distances between vertices
* `vweights::Vector{Float64}` landmarks total weights - used only with landmarks approximation
* `init_vweights::Vector{Float64}` vector with original (full) vertices weights - used only with landmarks approximation
* `v_to_l::Vector{Int}` mapping from vertices to landmarks (landmarks membership) - used only with landmarks approximation
* `init_edges::Array{Int,2}` array with original (full) graph edges - used only with landmarks approximation
* `init_eweights::Vector{Float64}` vector with original (full) edges weights - used only with landmarks approximation
* `init_embed::Matrix{Float64}` array with original embedding for full graph - used only with landmarks approximation
* `split::Bool` indicator for splitting JS divergence score (global score)
* `seed::Int` RNG seed for local measure score
* `auc_samples::Int` no. samples for local measure score
* `verbose::Bool` verbose switch, if true prints additional processing information
"""
function wGCL_directed(edges::Array{Int,2}, eweights::Vector{Float64}, comm::Matrix{Int},
            embed::Matrix{Float64}, distances::Vector{Float64}, vweights::Vector{Float64},
            init_vweights::Vector{Float64}, v_to_l::Vector{Int}, init_edges::Array{Int,2},
            init_eweights::Vector{Float64}, init_embed::Matrix{Float64},
            split::Bool, seed::Int=-1, auc_samples::Int=10000, verbose::Bool=false)
    # Default values
    delta = 0.001 # desired precision of degree estimation in Chung Lu model
    AlphaMax = 10.0 # upper bound of alpha search
    AlphaStep = 0.25 # step in alpha search
    alpha_div_counter = alpha_auc_counter = 5 # early stopping threshold (iterations of alpha search without improvement)
    skip_div = skip_auc = false # skippig flags based on early stopping counters

    no_vertices = maximum(edges)
    no_edges = size(edges,1)
    verbose && println("auc_samples: $auc_samples")
    landmarks = !isempty(v_to_l)

    verbose && println("Graph has $no_vertices vertices and $no_edges edges")
    landmarks && verbose && println("Original graph has $(maximum(init_edges)) vertices and $(size(init_edges,1)) edges")

    # Read communities
    @assert size(comm,1) == no_vertices "No. communities not matching no. vertices"
    n_parts = maximum(comm)
    verbose && println("Graph has $n_parts communities")

    # Compute degrees
    degree_in = zeros(no_vertices)
    degree_out = zeros(no_vertices)
    star_check = zeros(Int64,no_vertices)
    for i in 1:no_edges
        a = eweights[i]
        v1 = edges[i,1]
        v2 = edges[i,2]
        degree_out[v1]+=a
        degree_in[v2]+=a
        star_check[v1]+=1
        star_check[v2]+=1
    end

    # Check for edge case - star graph
    is_star = false
    # Star based on either in or out edges
    if !isnothing(findfirst(==(no_vertices-1),star_check)) && sum(star_check) == 2*(no_vertices-1)
        is_star = true
    # Star based on both in and out edges
    elseif !isnothing(findfirst(==(2*(no_vertices-1)),star_check)) && sum(star_check .== 2) == no_vertices-1
        is_star = true
    end
    verbose && is_star && println("Graph is a star in respect to either in or out edges")

    if is_star
        return [-1.0,zeros(5)...]
    end

    # Compute C-vector
    vect_len = Int(n_parts*n_parts)
    vect_C = zeros(Float64, vect_len)
    vect_B = zeros(Float64, vect_len)

    for i in 1:no_edges
        j = comm[edges[i,1]]
        k = comm[edges[i,2]]
        vect_C[(j-1)*n_parts + k] += eweights[i]
    end

    # Indicator - internal to a community
    vect_I = falses(vect_len)
    for i in 1:(n_parts+1):vect_len
        vect_I[i] = true
    end
    best_div = best_div_ext = best_div_int = best_auc_err = best_auc = typemax(Float64)
    best_alpha = best_alpha_auc = -1.0
    dim = size(embed,2)
    verbose && println("Embedding has $dim dimensions")
    # Loop over alpha's

    # Compute Euclidean distance vector D[] given embed and alpha
    p_len = no_vertices*(no_vertices+1) ÷ 2
    D = zeros(Float64, p_len)

    #Read distances
    @assert size(distances,1) == no_vertices "Distances vector length is not equal to no. vertices"
    for i in 1:no_vertices
        for j in i:no_vertices
            l = idx(no_vertices, i, j)
            if i == j
                D[l] = distances[i]
            else
                D[l] = dist(i,j,embed)
            end
        end
    end
    lo,hi = extrema(D)
    D = (D.-lo)./(hi-lo) # normalize to [0,1]

    adj_no_vertices = no_vertices
    adj_edges = edges
    adj_eweights = eweights
    if landmarks
        adj_no_vertices = length(init_vweights)
        adj_edges = init_edges
        adj_eweights = init_eweights
    end

    if landmarks
        adj_p_len = adj_no_vertices*(adj_no_vertices+1) ÷ 2
        full_graph_D = zeros(Float64, adj_p_len)
        for i in 1:adj_no_vertices
            for j in (i+1):adj_no_vertices
                l = idx(adj_no_vertices, i, j)
                full_graph_D[l] = dist(i,j,init_embed)
            end
        end
        lo,hi = extrema(full_graph_D)
        full_graph_D = (full_graph_D.-lo)./(hi-lo) # normalize to [0,1]
    end
    # Loop here - exclude Alpha = 0
    Tin = ones(no_vertices)
    Tin[degree_in.==0] .= 0.0
    Tout = ones(no_vertices)
    Tout[degree_out.==0] .= 0.0

    # Generate arrays of edges and non-edges
    NE = Tuple{Int64,Int64}[]
    for i in 1:adj_no_vertices
        for j in 1:adj_no_vertices
            if i!=j
                push!(NE,(i,j))
            end
        end
    end

    ## tuples of edges
    E = Tuple{Int64,Int64,Float64}[]
    for (i,e) in enumerate(eachrow(adj_edges))
        push!(E,tuple(e...,adj_eweights[i]))
    end

    ## tuples of non-edges
    NE = collect(setdiff(Set(NE),Set([e[1:2] for e in E])))

    for alpha in AlphaStep:AlphaStep:(AlphaMax+delta)
        write(stderr, ".")
        # Apply kernel (g(dist))
        GD = zeros(Float64, p_len)
        for i in 1:no_vertices
            for j in i:no_vertices
                k = idx(no_vertices,i,j)
                GD[k] = (1-D[k])^alpha # transform w.r.t. alpha
            end
        end
        # Learn GCL model numerically
        diff = 1.0
        epsilon = 0.9
        while diff > delta # stopping criterion
            Sin = zeros(no_vertices)
            Sout = zeros(no_vertices)
            for i in 1:no_vertices
                for j in i:no_vertices
                    k = idx(no_vertices, i, j)
                    tmp1 = Tin[i]*Tout[j]*GD[k]
                    tmp2 = Tin[j]*Tout[i]*GD[k]
                    Sin[i] += tmp1
                    Sin[j] += tmp2
                    Sout[i] += tmp2
                    Sout[j] += tmp1
                end
            end

            f = 0.0
            for i in 1:no_vertices
                if degree_in[i]>0
                    Tin[i] += epsilon*Tin[i]*(degree_in[i]/Sin[i]-1.0)
                    f = max(f,abs(degree_in[i]-Sin[i])) # convergence w.r.t. degrees
                end
                if degree_out[i]>0
                    Tout[i] += epsilon*Tout[i]*(degree_out[i]/Sout[i]-1.0)
                    f = max(f,abs(degree_out[i]-Sout[i])) # convergence w.r.t. degrees
                end
            end
            if f > diff
                epsilon *= 0.99
            end
            diff = f
            # verbose && println("diff= $diff")
        end

        # Compute probas P[]
        P = zeros(Float64,2*p_len)
        for i in 1:no_vertices
            for j in 1:no_vertices
                a, b = extrema([i,j])
                P[no_vertices*(i-1)+j] = Tout[i]*Tin[j]*GD[idx(no_vertices,a,b)]
            end
        end

        if !skip_auc
            pos = zeros(auc_samples)
            neg = zeros(auc_samples)
            auc_weights = zeros(auc_samples)
            if landmarks
                ## random positive cases
                seed != -1 && Random.seed!(seed)
                for (ind, edge) in enumerate(sample(E,auc_samples,replace=true))
                    i, j, w = edge
                    a, b = extrema((i,j))
                    adj_Tout = Tout[v_to_l[i]]*init_vweights[i]/vweights[v_to_l[i]]
                    adj_Tin = Tin[v_to_l[j]]*init_vweights[j]/vweights[v_to_l[j]]
                    pos[ind] = adj_Tout*adj_Tin*((1-full_graph_D[idx(adj_no_vertices,a,b)])^alpha)
                    auc_weights[ind] = w
                end
                ## random negative cases
                seed != -1 && Random.seed!(seed)
                for (ind, edge) in enumerate(sample(NE,auc_samples,replace=true))
                    i, j = edge
                    a, b = extrema(edge)
                    adj_Tout = Tout[v_to_l[i]]*init_vweights[i]/vweights[v_to_l[i]]
                    adj_Tin = Tin[v_to_l[j]]*init_vweights[j]/vweights[v_to_l[j]]
                    neg[ind] = adj_Tout*adj_Tin*((1-full_graph_D[idx(adj_no_vertices,a,b)])^alpha)
                end
            else
                ## random positive cases
                seed != -1 && Random.seed!(seed)
                for (ind, edge) in enumerate(sample(E,auc_samples,replace=true))
                    i, j, w = edge
                    pos[ind] = P[no_vertices*(i-1)+j]
                    auc_weights[ind] = w
                end
                pos = [P[no_vertices*(i-1)+j] for (i,j) in sample(E,auc_samples,replace=true)]
                ## random negative cases
                seed != -1 && Random.seed!(seed)
                neg = [P[no_vertices*(i-1)+j] for (i,j) in sample(NE,auc_samples,replace=true)]
            end

            # 1 - AUC
            auc = 1 - sum((pos .> neg).*auc_weights)/sum(auc_weights)

            if auc < best_auc
                best_auc = auc
                best_auc_err = 1.96 * sqrt(auc*(1-auc)/auc_samples) ## error from 95% CI
                best_alpha_auc = alpha
                alpha_auc_counter = 5
            else
                alpha_auc_counter -= 1
                skip_auc = alpha_auc_counter == 0
            end
        end

        if !skip_div
            # Compute B-vector given P[] and comm[]
            vect_B = zeros(vect_len)
            for i in 1:no_vertices
                for j in 1:no_vertices
                    m = (comm[i]-1)*n_parts + comm[j]
                    vect_B[m] += P[no_vertices*(i-1)+j]
                end
            end
            if !split
                f = JS(vect_C, vect_B, Bool[], true)
            else
                x = JS(vect_C, vect_B, vect_I, true)
                y = JS(vect_C, vect_B, vect_I, false)
                f = (x+y)/2.0
            end
            if f < best_div
                best_div = f
                best_alpha = alpha
                best_div_ext = !split ? 0.0 : x
                best_div_int = !split ? 0.0 : y
                alpha_div_counter = 5
            else
                alpha_div_counter -= 1
                skip_div = alpha_div_counter == 0
            end
        end
        skip_div && skip_auc && break
    end
    write(stderr, "\n")
    return [best_alpha, best_div, best_div_ext, best_div_int, best_alpha_auc, best_auc, best_auc_err]
end


text:
using Serialization

using TodoGrabber

const cachefile = "./tests/serial_cache.bin"

function get_filenames(fpath::String)
    files = String[]
    for f in readdir(abspath(fpath), join=true)
        push!(files, f)
    end
    return files
end

function write_todos(tododict::Union{Dict{String,Vector{Todo}},Nothing},
                     targetfile::String)
    if tododict === nothing
        printstyled(stderr, "No todos found", color=:red)
        return nothing
    end
    open(targetfile, "w") do tfile
        for (fname, todo_vec) in tododict
            write(tfile, string("* ", fname, "\n"))
            for todo in todo_vec
                write(tfile, string("** [[TODO]] ", todo.ID, " ", todo.Suffix,
                                    "\n"))
            end
        end
    end
    return nothing
end

function save_dict(tododict::Dict{String,Vector{Todo}}; targetfile::String="")
    if isempty(targetfile)
        serialize(cachefile, tododict)
    else
        serialize(targetfile, tododict)
    end
    return nothing
end

function grab_dir(directory::String; target_file::String="")
    filenames = get_filenames(directory)
    todos = find_todos(filenames)
    save_dict(todos)
    if isempty(target_file)
        println(stdout, "[NOTE] No output file given. Writing on temp_stack.org")
        write_todos(todos, "./tests/temp_stack.org")
    else
        println(stdout, "Writing on $target_file")
        write_todos(todos, target_file)
    end
    return nothing
end

function find_todos(target_files::Vector{String})
    todo_dict = Dict{String,Vector{Todo}}()
    for filename in target_files
        captured = Vector{Todo}(undef, 0)
        open(filename) do file
            match_lines!(captured, filename, file, todo_dict)
        end # file loop
        if isempty(captured) # only push todos when they exist
            continue
        end
        todo_dict[filename] = captured
    end
    if isempty(todo_dict)
        return nothing
    end
    return todo_dict
end

function match_lines!(captured::Vector{Todo}, filename::String, file::IOStream,
                      todo_dict::Union{Dict{String,Vector{Todo}},Nothing})

    preffix = ""
    line::Int64 = 0
    sizet::Int64 = 0
    search_body::Bool = false
    temp_todo_vec = Vector{Todo}(undef, 256)
    println(stdout, "working on ", filename)
    while !eof(file)
        line += 1
        line_text = readline(file)
        m = match(r"^(.*)TODO: (.*)$", line_text)
        if m !== nothing
            # TODO: Rethink about stripping
            cleaned_capture = strip(m.captures[2], ['-','*','>',' ', '/'])
            if !haskey(todo_dict, filename) ||
                isempty(findall(x -> x == cleaned_capture, getindex(todo_dict,
                                                                    filename)))
                sizet += 1
                preffix = m.captures[1]
                temp_todo_vec[sizet] = Todo(preffix,
                                              String(cleaned_capture),
                                              filename,
                                              "",
                                              line,
                                              sizet)

                search_body = true
            end
        elseif !isempty(line_text) && !isempty(preffix) && search_body
            if preffix == line_text[1:sizeof(preffix)]
                temp_todo_vec[sizet].Body *= split(line_text, preffix)[2] * " "
            else
                preffix = ""
                search_body = false
            end
        end
    end # while loop
    copy!(captured, temp_todo_vec[1:sizet])
end # match_lines


text:
function molecule(::Molecule"phosphaethene")
  return """
C	0.0566650	1.0803320	0.0000000
P	0.0566650	-0.6006730	0.0000000
H	-0.8348350	1.6951630	0.0000000
H	1.0025300	1.6111320	0.0000000
H	-1.3576590	-0.7781910	0.0000000
  """
end


text:
using Flux: @treelike
using MacroTools: @forward

using ..Basic
using ..Basic: onehot, AbstractTransformer
using ..Stacks

struct Gpt <: AbstractTransformer
    ts::Stack
    drop::Dropout
end

@treelike Gpt

@forward Gpt.ts Base.getindex, Base.length

"""
    Gpt(size::Int, head::Int, ps::Int, layer::Int;
        act = gelu, pdrop = 0.1, attn_pdrop = 0.1)
    Gpt(size::Int, head::Int, hs::Int, ps::Int, layer::Int;
        act = gelu, pdrop = 0.1, attn_pdrop = 0.1)

the Generative Pretrained Transformer(GPT) model.

    (gpt::Gpt)(x::T, mask=nothing; all::Bool=false)

eval the gpt layer on input `x`. If length `mask` is given (in shape (1, seq_len, batch_size)), mask the attention with `mask`. Moreover, set `all` to `true` to get all
outputs of each transformer layer.
"""
function Gpt(size::Int, head::Int, ps::Int, layer::Int;
             act = gelu, pdrop = 0.1, attn_pdrop = 0.1)
    rem(size, head) != 0 && error("size not divisible by head")
    Gpt(size, head, div(size, head), ps, layer; act=act, pdrop=pdrop, attn_pdrop=attn_pdrop)
end

function Gpt(size::Int, head::Int, hs::Int, ps::Int, layer::Int;
             act = gelu, pdrop = 0.1, attn_pdrop = 0.1)
    Gpt(
        Stack(
            @nntopo_str("x':x => $layer"),
            [
                Transformer(size, head, hs, ps; future=false, act=act, pdrop=attn_pdrop)
                for i = 1:layer
            ]...
        ),
        Dropout(pdrop)
    )
end

function (gpt::Gpt)(x::T, mask=nothing; all::Bool=false) where T
    e = gpt.drop(x)
    t, ts = gpt.ts(e)
    t = mask === nothing ? t : t .* mask
    if all
        t, ts
    else
        t #size(t) == (size, seq_len, batch)
    end
end

"""
    lmloss(embed, onehot, encoding, mask)

compute the language modeling loss for Gpt, onehot is the onehot array of the origin
input sentence. encoding the output of Gpt, mask is the mask between input sentences.

"""
lmloss(embed::Embed{T}, o::OneHotArray, t::AbstractArray{T}, mask) where T = lmloss(embed, tofloat(T, o), t, mask)
function lmloss(embed::Embed{T}, et, t::AbstractArray{T, N}, mask) where {T,N}
    if N == 3
        t = t[:, 1:end-1, :]
        s = size(t)
        sim = logsoftmax(transpose(embed.embedding) * reshape(t, s[1], :)) #(vocab, seq_len*batch)
        sim = reshape(sim, :, s[2], s[3])
        return logcrossentropy(et[:, 2:end, :], sim, mask[:, 2:end, :])
    elseif N == 2
        t = t[:, 1:end-1]
        sim = logsoftmax(transpose(embed.embedding) * t)
        return logcrossentropy(et[:, 2:end], sim, mask[:, 2:end])
    end
end

function Base.show(io::IO, gpt::Gpt)
    hs = div(size(gpt.ts[1].mh.iqproj.W)[1], gpt.ts[1].mh.head)
    h, ps = size(gpt.ts[1].pw.dout.W)

    print(io, "Gpt(")
    print(io, "layers=$(length(gpt.ts)), ")
    print(io, "head=$(gpt.ts[1].mh.head), ")
    print(io, "head_size=$(hs), ")
    print(io, "pwffn_size=$(ps), ")
    print(io, "size=$(h))")
end


text:
# The now-defunct JuMPeR formulations of the problems are included here!

""" Robust facility location model, with JuMPeR. """
function robust_facility_model(c::Matrix, f::Vector, rho::Real, Gamm::Real)
    n, m = size(c)
    @assert length(f) == n
    model = RobustModel(solver = GurobiSolver())
    # model = Model(solver = GurobiSolver())

    # VARIABLES
    @variable(model, x[1:n], Bin)                 # Facility locations
    @variable(model, y[1:n, 1:m] >= 0)            # Flow variables
    @uncertain(model, -rho <= z[1:m] <= rho)      #  Demand uncertainties
    @constraint(model, norm(z, 1) <= Gamm)

    # CONSTRAINTS
    @constraint(model, [j=1:m], sum(y[i, j] for i = 1:n) >= d[j] + (P*z)[j]) # demand constraint
    @constraint(model, [i=1:n], sum(y[i, j] for j = 1:m) <= s[i] * x[i]) # capacity constraint
    # OBJECTIVE
    @objective(model, Min,
            sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * y[i, j] for i=1:n, j=1:m))
    return model, x, y
end

""" Adaptive facility location model, with JuMPeR and relaxation of binary variables. """
function adaptive_facility_model(c::Matrix, f::Vector, rho::Real, Gamm::Real)
    n, m = size(c)
    @assert length(f) == n
    model = RobustModel(solver = GurobiSolver())

    # VARIABLES
    @variable(model, 0 <= x[1:n] <= 1)                 # Facility locations
    @uncertain(model, -rho <= z[1:m] <= rho)      #  Demand uncertainties
    @constraint(model, norm(z, 1) <= Gamm)
    @adaptive(model, y[i=1:n, j=1:m] >= 0, policy = Affine, depends_on=z) # Flow variables

    # CONSTRAINTS
    @constraint(model, [j=1:m], sum(y[i, j] for i = 1:n) >= d[j] + (P*z)[j]) # demand constraint
    @constraint(model, [i=1:n], sum(y[i, j] for j = 1:m) <= s[i] * x[i]) # capacity constraint
    # OBJECTIVE
    @variable(model, F)
    @constraint(model, F >= sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * y[i, j] for i=1:n, j=1:m))
    @objective(model, Min, F)
    return model, x, y
end

# Adaptive solution with greedy binary heuristic (3.4 and 3.6)
function iterate_adaptive_model(c, f, rho, Gamm)
    ints = []
    xvals = zeros(size(c, 1))
    while !all([val in [0,1] for val in xvals]) || isempty(ints)
        model, x, y = adaptive_facility_model(c, f, rho, Gamm)
        for int in ints
            @constraint(model, x[int] == 1)
        end
        optimize!(model)
        xvals = value.(x)
        nonint = [!(var in [0,1]) for var in xvals]
        if sum(nonint) == 0
            return model, x, y, ints
        end
        max_val, max_idx = findmax(nonint.*xvals)
        push!(ints, max_idx)
    end
    return
end

# model, x, y, ints = iterate_adaptive_model(c, f, rho, Gamm)

# # Plot generation for adaptive solution (3.4 and 3.6)
# cost = objective_value(model)
# nonints = [val for val in 1:n if !(val in ints)]
# plot_model, x, y = facility_model(c, f)
# @constraint(plot_model, x[ints] .== 1)
# @constraint(plot_model, x[nonints] .== 0)
# optimize!(plot_model)
# plt = plot_solution(model, x, y, cost)
# @show plt


text:
# QP example
#
#    minimize 1/2*(2 x^2 + y^2 + xy) + x + y
#
#       s.t.  x, y >= 0
#             x + y = 1
#
#    solution: (0, 1), objv = 1.875
#
using Xpress, Compat.Test
using Compat.SparseArrays
using Compat.LinearAlgebra

@testset "QP 1" begin
    model = Xpress.Model("qp_02")

    add_cvars!(model, [1., 1.], 0., Inf)

    add_qpterms!(model, [1, 1, 2], [1, 2, 2], [2., 1., 1.])
    add_constr!(model, [1., 1.], '=', 1.)

    @test Xpress.getq_upper(model) == triu(sparse([2. 1.; 1. 1.]))
    @test Xpress.getq(model) == sparse([2. .5; .5 1.])
    @test Xpress.get_obj(model) == [1, 1]
    @test Xpress.get_constrmatrix(model) == sparse([1 1])
    @test Xpress.get_rhs(model) == [1]
    @test Xpress.get_lb(model) == [0, 0]
    @test Xpress.get_ub(model) == [1e20, 1e20]
    @test Xpress.get_sense(model) == [:(==)]

    optimize(model)

    ans = get_optiminfo(model)
    @test ans.status_lp == :optimal

    @test isapprox(get_solution(model), [0.0, 1]; atol = 1e-2)

    @test isapprox(get_objval(model), 1.5; atol = 1e-3)
end


text:
function check_open_orders_klines(config, open_orders, closed_orders, log, dt, x, high, low)

end

function check_open_orders(config, open_orders, closed_orders, log, current)

    if length(keys(open_orders)) > 0
        extra_order_count = parse(Int64, config["extra_order_count"])
        extra_order_exec  = log["extra_order_exec"]
        order_volume = config["order_volume"]
        extra_order_martingale = config["extra_order_martingale"]
        precision = config["precision"]
        extra_order_step = config["extra_order_step"]
        sell_precision = config["sell_precision"]
        stop_loss = parse(Float64, config["stop_loss"])
        fee = config["fee"]
        take_profit = config["take_profit"]

        if config["dry_mode"] == "true"

            if current >= parse(Float64, open_orders["SELL"]["price"]) ### TAKE PROFIT
#                 @info("Take profit was executed")
                open_orders["SELL"]["status"] = "FILLED"
                open_orders["SELL"]["updateTime"] = datetime2unix(now())*1000
                open_orders["SELL"]["executedQty"] = open_orders["SELL"]["origQty"]
                closed_orders[string(now())] = copy(open_orders["SELL"])
                open_orders["BUY"]["status"] == "CANCELLED"
                closed_orders[string(now() - Second(1))] = open_orders["BUY"]
                log["extra_order_exec"] = 0
                log["trade_ongoing"] = false

                open_orders = Dict()

                open("closed_orders.json", "w") do f
                    JSON.print(f, closed_orders, 4)
                end
                open("open_orders.json", "w") do f
                    JSON.print(f, open_orders, 4)
                end
                open("log.json", "w") do f
                    JSON.print(f, log, 4)
                end
                log = Pinguu.daily_log(closed_orders, config)
                #send_notification(config, "take profit", log);

            elseif extra_order_exec < extra_order_count
                if current <= parse(Float64, open_orders["BUY"]["price"]) ### EXTRA ORDER
#                 @info("Extra order was executed")
                    log["extra_order_exec"] += 1
                    closed_orders = JSON.parsefile("closed_orders.json");
                    open_orders = JSON.parsefile("open_orders.json");
                    executed_order = deepcopy(open_orders["BUY"])
                    executed_order["executedQty"] = executed_order["origQty"]
                    executed_order["status"] = "FILLED"
                    executed_order["updateTime"] = datetime2unix(now())*1000
                    failed_order = deepcopy(open_orders["SELL"])
                    failed_order["status"] = "CANCELLED"
                    closed_orders[string(unix2datetime(executed_order["updateTime"]/1000))] = executed_order
                    closed_orders[string(now() - Second(1))] = failed_order
                    open_orders["BUY"]["time"] = datetime2unix(now())*1000
                    open_orders["BUY"]["updateTime"] = datetime2unix(now())*1000
                    open_orders["BUY"]["orderId"] = rand(1:1:1e7)
                    old_quant = parse(Float64, open_orders["BUY"]["origQty"])
                    old_price = parse(Float64, open_orders["BUY"]["price"])
                    new_quant = floor(order_volume * (extra_order_martingale)^extra_order_exec / parse(Float64, open_orders["BUY"]["price"]), digits=precision)
                    open_orders["BUY"]["origQty"] = string(new_quant)
                    open_orders["BUY"]["executedQty"] = "0.0"
                    open_orders["BUY"]["status"] = "NEW"
                    if extra_order_exec != extra_order_count
                        open_orders["BUY"]["price"] = string(round(old_price * (1 - extra_order_step), digits=sell_precision))
                    else
                        open_orders["BUY"]["status"] = "STOPLOSSLIMIT"
                        open_orders["BUY"]["price"] = string(round(parse(Float64, open_orders["SELL"]["price"]) * (100 - stop_loss)/100, digits=sell_precision))
                    end

                    oriqQty_sell = parse(Float64, open_orders["SELL"]["origQty"])
                    newQty_sell  = floor(oriqQty_sell  + old_quant, digits=precision)
                    origPrice_sell = parse(Float64, open_orders["SELL"]["price"])
                    qty_old = 0
                    for i in 0:1:(extra_order_exec-1)
                        qty_old += extra_order_martingale^i
                    end
                    qty_new = extra_order_martingale^extra_order_exec
                    new_price = (qty_old * origPrice_sell + qty_new * old_price * (1 + fee + take_profit/100)) / (qty_old + qty_new)

                    open_orders["SELL"]["time"] = datetime2unix(now())*1000
                    open_orders["SELL"]["updateTime"] = datetime2unix(now())*1000
                    open_orders["SELL"]["orderId"] = rand(1:1:1e7)
                    open_orders["SELL"]["origQty"] = string(newQty_sell)
                    open_orders["SELL"]["price"] = string(round(new_price, digits=sell_precision))

                    open("closed_orders.json", "w") do f
                        JSON.print(f, closed_orders, 4)
                    end
                    open("open_orders.json", "w") do f
                        JSON.print(f, open_orders, 4)
                    end
                    open("log.json", "w") do f
                        JSON.print(f, log, 4)
                    end
                end
#                 send_notification(config, "extra order", log);
            elseif extra_order_exec == extra_order_count ### STOP LOSS
                if current <= parse(Float64, open_orders["SELL"]["price"]) * (100 - parse(Float64, stop_loss))/100
                    open_orders["SELL"]["status"] = "FILLED"
                    open_orders["SELL"]["updateTime"] = datetime2unix(now())*1000
                    open_orders["SELL"]["executedQty"] = open_orders["SELL"]["origQty"]
                    open_orders["SELL"]["price"] = string(round(current, digits=sell_precision))
                    closed_orders[string(now())] = copy(open_orders["SELL"])
                    log["extra_order_exec"] = 0
                    log["trade_ongoing"] = false

                    open_orders = Dict()

                    open("closed_orders.json", "w") do f
                        JSON.print(f, closed_orders, 4)
                    end
                    open("open_orders.json", "w") do f
                        JSON.print(f, open_orders, 4)
                    end
                    open("log.json", "w") do f
                        JSON.print(f, log, 4)
                    end
                    log = Pinguu.daily_log(closed_orders, config)
                    #send_notification(config, "stop loss", log);
                end
            end
        else
            order = client.get_order(
                symbol=config["symbol"],
                orderId=open_orders["SELL"]["orderId"])
            if order["status"] == "FILLED"
#                 @info("Take profit was executed")
                closed_orders[string(unix2datetime(order["time"]/1000))] = order
                result = client.cancel_order(
                    symbol=config["symbol"],
                    orderId=open_orders["BUY"]["orderId"]);
                open_orders = Dict()
            end
            #
            #
            #
            #
            #
            #
            #
            #
            #
            #
            #
            #
            #
            #
            #
            #
            #
        end
    end
    return open_orders, closed_orders, log
end


text:
# https://adventofcode.com/2020/day/24

const sides = Symbol[:N, :NE, :SE, :S, :SW, :NW]

function parse_line(line::AbstractString)
    coords = Symbol[]

    i = 1
    while i <= length(line)
        d = string(line[i])
        if d in ("s", "n")
            i += 1
            d *= line[i]
        end
        push!(coords, Symbol(uppercase(d)))
        i += 1
    end

    return coords
end

@assert all(parse_line("esenee") .== [:E, :SE, :NE, :E])

example = "
sesenwnenenewseeswwswswwnenewsewsw
neeenesenwnwwswnenewnwwsewnenwseswesw
seswneswswsenwwnwse
nwnwneseeswswnenewneswwnewseswneseene
swweswneswnenwsewnwneneseenw
eesenwseswswnenwswnwnwsewwnwsene
sewnenenenesenwsewnenwwwse
wenwwweseeeweswwwnwwe
wsweesenenewnwwnwsenewsenwwsesesenwne
neeswseenwwswnwswswnw
nenwswwsewswnenenewsenwsenwnesesenew
enewnwewneswsewnwswenweswnenwsenwsw
sweneswneswneneenwnewenewwneswswnese
swwesenesewenwneswnwwneseswwne
enesenwswwswneneswsenwnewswseenwsese
wnwnesenesenenwwnenwsewesewsesesew
nenewswnwewswnenesenwnesewesw
eneswnwswnwsenenwnwnwwseeswneewsenese
neswnwewnwnwseenwseesewsenwsweewe
wseweeenwnesenwwwswnew
"

function part1(input)

    for row in split(input, '\n', keepempty=false)
        coords = parse_line(row)
    end

    return
end

@assert part1(example) == 10

# test =
# println("Part 1: $(result1 = part1(test))")
# println("Part 2: $(result2 = part2(test))")

# @assert result1 ==
# @assert result2 ==

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
# Autogenerated wrapper script for LCIO_Julia_Wrapper_jll for i686-linux-gnu-cxx11-julia_version+1.6.0
export lciowrap

using libcxxwrap_julia_jll
using LCIO_jll
JLLWrappers.@generate_wrapper_header("LCIO_Julia_Wrapper")
JLLWrappers.@declare_library_product(lciowrap, "liblciowrap.so")
function __init__()
    JLLWrappers.@generate_init_header(libcxxwrap_julia_jll, LCIO_jll)
    JLLWrappers.@init_library_product(
        lciowrap,
        "lib/liblciowrap.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export qbf_operator

function qbf_operator(basis::Dictionary, degree::Int, oversamplingfactor)
    T = domaintype(basis)
    qbf = QuadQBF{T}(degree, oversamplingfactor*(degree+1))
    _qbf_operator(basis, qbf, oversamplingfactor)
end

function qbf_operator(basis::Dictionary, coef::AbstractVector, oversamplingfactor::Int)
    T = domaintype(basis)
    qbf = QuadQBF{T}(coef; oversamplingfactor=oversamplingfactor)
    _qbf_operator(basis, qbf, oversamplingfactor)
end

function _qbf_operator(basis::Dictionary, qbf::QuadQBF{T}, oversamplingfactor) where {T}
    grid = sampling_grid(basis, oversamplingfactor=oversamplingfactor)
    gb = GridBasis{coefficienttype(basis)}(grid)
    step = oversamplingfactor
    offset = -(length(qbf.x)>>1)
    HorizontalBandedOperator(gb, basis, qbf.w / sqrt(convert(T, length(basis))), step, offset)
end


text:
"""
    to_univariate_finite(scores::Scores)

Convert normalized scores to a vector of univariate finite distributions.

Parameters
----------
    scores::[`Scores`](@ref)

Returns
----------
    fit::UnivariateFiniteVector{OrderedFactor{2}}
The learned model of the given detector, which contains all the necessary information for later prediction and the
achieved outlier scores of the given input data `X`.
"""
function to_univariate_finite(scores::Scores)
    MLJ.UnivariateFinite([CLASS_NORMAL, CLASS_OUTLIER], scores; augment=true, pool=missing, ordered=true)
end

"""
    to_categorical(classes::Labels)

Convert a vector of classes (with possible missing values) to a categorical vector.

Parameters
----------
    classes::[`Labels`](@ref)

Returns
----------
    fit::CategoricalVector{Union{Missing, String},UInt32}
The learned model of the given detector, which contains all the necessary information for later prediction and the
achieved outlier scores of the given input data `X`.
"""
function to_categorical(classes::Labels)
    # explicit cast to Vector{Union{String, Missing}} in case only missing values are passed
    MLJ.categorical(Vector{Union{String, Missing}}(classes), ordered=true, levels=[CLASS_NORMAL, CLASS_OUTLIER])
end


text:
using CpelNano
using Test

@testset "Computations check" begin

    ## Marginal expectations

    # Transfer matrix objects
    @test CpelNano.get_u(0.0)==[1.0,1.0]
    @test CpelNano.get_W(0.0,0.0,0.0)==[1.0 1.0;1.0 1.0]

    # Partition function
    L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    u1 = CpelNano.get_u(αs[1]); uN = CpelNano.get_u(αs[end]);
    Ws = [CpelNano.get_W(αs[n],αs[n+1],βs[n]) for n=1:length(βs)];
    Z = CpelNano.get_Z(u1,uN,Ws)

    # Partition function
    @test Z==1024.0

    # Expected value of X
    @test CpelNano.get_E_X(u1,uN,Ws,Z)==zeros(10)

    # Expected value of XX
    @test CpelNano.get_E_XX(u1,uN,Ws,Z)==zeros(9)

    # Log-sum trick
    L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    logu1 = CpelNano.get_log_u(αs[1]);
    loguN = CpelNano.get_log_u(αs[end]);
    logWs = [CpelNano.get_log_W(αs[n],αs[n+1],βs[n]) for n=1:length(βs)];

    # Partition function
    logZ = CpelNano.get_log_Z(logu1,loguN,logWs)
    @test logZ==log(1024)

    # Expected value of X
    @test isapprox(CpelNano.get_E_X_log(logu1,loguN,logWs,logZ),zeros(10);atol=1e-5)

    # Expected value of XX
    @test isapprox(CpelNano.get_E_XX_log(logu1,loguN,logWs,logZ),zeros(9);atol=1e-5)

    ## Conditional expectations

    # Regular transfer matrix computation
    L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    obs = fill(CpelNano.MethCallCpgGrp(-50.0,-50.0),L);
    u1 = CpelNano.get_uc(αs[1],obs[1]); uN = CpelNano.get_uc(αs[end],obs[end]);
    Ws = [CpelNano.get_Wc(αs[n],αs[n+1],βs[n],obs[n],obs[n+1]) for n=1:length(βs)];

    # Partition function
    Zc = CpelNano.get_Zc(u1,uN,Ws);
    @test Zc==7.295566240503076e-215

    # Expected value of X
    @test CpelNano.get_Ec_X(u1,uN,Ws,Zc)==zeros(10)

    # Expected value of XX
    @test CpelNano.get_Ec_XX(u1,uN,Ws,Zc)==zeros(9)

    # Log-sum trick
    L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    obs = fill(CpelNano.MethCallCpgGrp(-100.0,-100.0),L);
    logu1 = CpelNano.get_log_uc(αs[1],obs[1]);
    loguN = CpelNano.get_log_uc(αs[end],obs[end]);
    logWs = [CpelNano.get_log_Wc(αs[n],αs[n+1],βs[n],obs[n],obs[n+1]) for n=1:length(βs)];

    # Partition function
    logZc = CpelNano.get_log_Zc(logu1,loguN,logWs)
    @test logZc==-993.0685281944009

    # Expected value of X
    @test isapprox(CpelNano.get_Ec_X_log(logu1,loguN,logWs,logZc),zeros(10);atol=1e-5)

    # Expected value of XX
    @test isapprox(CpelNano.get_Ec_XX_log(logu1,loguN,logWs,logZc),zeros(9);atol=1e-5)

end


text:
if VERSION < v"0.7.0"
  using Base.REPLCompletions: latex_symbols
else
  using REPL.REPLCompletions: latex_symbols
end

open(joinpath(dirname(@__FILE__), "..", "latex.ts"), "w") do f
  println(f, "export const latexSymbols = {")
  # sort by name length and name content
  symbols = sort(collect(latex_symbols))
  symbols = sort(symbols, by=x->length(x[1]))
  for (name, sym) in symbols
    println(f, "    '\\$name': '$sym',")
  end
  println(f, "};")
end

max_stars_repo_name: colinfang/unicode-latex

max_issues_repo_name: colinfang/unicode-latex

max_forks_repo_name: colinfang/unicode-latex


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

TANH_APPROX = nothing

"""
Load tanh approximator from bson if uninitiated.
"""
function load_tanh()
    if TANH_APPROX === nothing
        global TANH_APPROX = BSON.load(joinpath(@__DIR__, "tanh.bson"), @__MODULE__)[:network]
    end
    return TANH_APPROX
end

"""
Retrain neural network univariate tanh approximator. Ensures f(x) = f(-x) and |f(x)| <= 1.
"""
function retrain_tanh(; depth::Int64=1, width::Int64=10, num_iter::Int64=10000, num_samples::Int64=1000)
    # Create model
    layers = Dense[]
    push!(layers, Dense(1, width, relu))
    for _ in 1:(depth - 1)
        push!(layers, Dense(width, width, relu))
    end
    push!(layers, Dense(width, 1))
    model = Chain(layers...)
    ps = Flux.params(model)
    opt = AdaBelief(1e-3)

    # Train network
    p = Progress(num_iter)
    loss = Ref(0.0)
    for _ in 1:num_iter
        x = randn(1, num_samples)
        gs = gradient(ps) do
            loss[] = sum((tanh.(x) .- model(x)).^2)
        end
        Flux.update!(opt, ps, gs)
        ProgressMeter.next!(p; showvalues = [(:loss, loss[])])
    end

    # Antisymmetrize
    d1 = Dense(reshape([1.0f0, -1.0f0], 2, 1))
    d2 = Dense(reshape([0.5f0, -0.5f0], 1, 2))
    layers = Dense[]
    for d in m.layers
        W = [d.W zero(d.W); zero(d.W) d.W]
        b = [d.b; d.b]
        push!(layers, Dense(W, b, d.σ))
    end

    # Clamp
    u, v = -1.0f0, 1.0f0
    d3 = Dense(reshape([1.0f0, 1.0f0], 2, 1), [-u, -v], relu)
    d4 = Dense(reshape([1.0f0, -1.0f0], 1, 2), [u])

    global TANH_APPROX = Chain(d1, layers..., d2, d3, d4)
end

"""
Extended neural network containing metadata about transformations.
"""
Base.@kwdef mutable struct ExtendedNetwork
    nnet::Network
    step_up::Int64      # index of step-up layer
    step_down::Int64    # index of step-down layer
    scale::Int64 = 0.0  # distance scale (for reward shaping)
end

"""
Convert Flux activation functions to nnet functions. Limited to `identity` and `relu`.
"""
function convert_f(f::Function)
    fsym = Symbol(f)
    if fsym == :identity
        return NeuralVerification.Id()
    elseif fsym == :relu
        return NeuralVerification.ReLU()
    else
        error("Unable to convert $f.")
    end
end

"""
Pad front of matrix with zeros.
"""
function fpad(M::Matrix{T}, dims::Tuple{Int64,Int64}) where T
    i, j = size(M)
    ip, jp = dims
    Mp = hcat(zeros(T, i + ip, jp), vcat(zeros(T, ip, j), M))
    return Mp
end

"""
Pad front of vector with zeros.
"""
function fpad(v::Vector{T}, sz::Int64) where T
    return vcat(zeros(T, sz), v)
end

"""
Pad front of matrix with identity.
"""
function fpad_id(M::Matrix{T}, sz::Int64) where T
    Mp = fpad(M, (sz, sz))
    Mp[1:sz,1:sz] = Matrix{T}(I, sz, sz)
    return Mp
end

"""
Pad front of Flux.Dense layer and convert to nnet layer.
"""
function fpad(l::Dense, sz::Int64)
    W = fpad_id(l.weight, sz)
    b = fpad(l.bias, sz)
    σ = convert_f(l.σ)
    return Layer(W, b, σ)
end

"""
Construct block diagonal matrix from matrix arguments.
"""
function block_diag(matrices...)
    m = sum((M -> size(M, 1)).(matrices))
    n = sum((M -> size(M, 2)).(matrices))
    B = zeros(Float32, m, n)

    i, j = 1, 1
    for M in matrices
        mi, mj = size(M)
        B[i:i+mi-1, j:j+mj-1] = M
        i += mi
        j += mj
    end
    return B
end

"""
Construct block diagonal matrix repetition.
"""
function rep_diag(M::Matrix, n::Int64)
    A = [[i == j ? M : zero(M) for i in 1:n] for j in 1:n]
    return reduce(hcat, reduce(vcat, m) for m in A)
end

"""
Construct network that squashes input elementwise with tanh approximator, provided lower
and upper bounds.
"""
function squash_model(ls::Vector{Float64}, us::Vector{Float64})
    # Stack of tanh approximators
    n = length(ls)
    layers = Dense[]
    for d in load_tanh().layers
        W = rep_diag(d.W, n)
        b = repeat(d.b, n)
        push!(layers, Dense(W, b, d.σ))
    end

    # Rescale
    ls, us = Float32.(ls), Float32.(us)
    W = diagm((us - ls) / 2)
    b = (us + ls) / 2
    push!(layers, Dense(W, b))

    return Chain(layers...)
end

"""
Construct layer that computes mean of previous layer.
"""
function mean_layer(sz::Int64)
    W = Float32.(ones(1, sz) / sz)
    b = zeros(Float32, 1)
    σ = NeuralVerification.Id()
    return Layer(W, b, σ)
end

"""
Construct layer that computes elementwise absolute value of previous layer.
"""
function abs_layers(sz::Int64)
    M = Matrix{Float32}(I, sz, sz)
    W1 = vcat(M, -M)
    b1 = zeros(Float32, 2 * sz)
    σ1 = NeuralVerification.ReLU()

    W2 = hcat(M, M)
    b2 = zeros(Float32, sz)
    σ2 = NeuralVerification.Id()

    return (Layer(W1, b1, σ1), Layer(W2, b2, σ2))
end

"""
Coalesce artificially constructed network, incorporating layers without activation functions
into subsequent layers. Note that this is not guaranteed to reduce the number of parameters,
only the number of layers.
"""
function coalesce(nnet::Network)
    layers = []
    W, b, stored = nothing, nothing, false
    for l in nnet.layers
        W = stored ? l.weights * W : l.weights
        b = stored ? l.weights * b + l.bias : l.bias
        stored = l.activation isa NeuralVerification.Id && l != nnet.layers[end]
        stored || push!(layers, Layer(W, b, l.activation))
    end
    return Network(layers)
end

"""
Convert policy from actor-critic object to neural network. Concatenation and squashing are
represented by ReLU configurations.
"""
function policy_network(ac::GlobalResult; act_mins::Vector{Float64}, act_maxs::Vector{Float64})
    # Dimensions and parameters
    n_obs = size(ac.pi.net.layers[1].weight, 2)
    n_act = size(ac.pi.mu_layer.weight, 1)
    layers = []

    # Add duplication layer
    M = Matrix{Float32}(I, n_obs, n_obs)
    W = vcat(M, M)
    b = zeros(Float32, 2*n_obs)
    σ = NeuralVerification.Id()
    push!(layers, Layer(W, b, σ))
    step_up = length(layers)

    # Add policy net
    push!(layers, fpad.(ac.pi.net, n_obs)...)
    push!(layers, fpad(ac.pi.mu_layer, n_obs))

    # Add action squash
    squash = squash_model(act_mins, act_maxs)
    push!(layers, fpad.(squash.layers, n_obs)...)
    step_down = length(layers)

    return ExtendedNetwork(nnet=Network(layers), step_up=step_up, step_down=step_down)
end

"""
Convert actor-critic object to extended neural network representing ensemble of values.
"""
function values_network(ac::GlobalResult; act_mins::Vector{Float64}, act_maxs::Vector{Float64})
    # Compute policy network
    network = policy_network(ac; act_mins=act_mins, act_maxs=act_maxs)
    layers = network.nnet.layers

    # Replicate output
    sz = length(layers[end].bias)
    n_qs = length(ac.qs)
    W = repeat(Matrix{Float32}(I, sz, sz), n_qs)
    b = zeros(Float32, n_qs*sz)
    σ = NeuralVerification.Id()
    push!(layers, Layer(W, b, σ))

    # Concatenate critics
    for ls in zip((q -> q.q).(ac.qs)...)
        W = block_diag((l -> l.weight).(ls)...)
        b = vcat((l -> l.bias).(ls)...)
        σ = convert_f(ls[1].σ)
        push!(layers, Layer(W, b, σ))
    end
    return network
end

"""
Convert actor-critic object to extended neural network representing mean of ensemble.
"""
function mean_network(ac::GlobalResult; act_mins::Vector{Float64}, act_maxs::Vector{Float64}, s::Float64=0.0)
    network = values_network(ac; act_mins=act_mins, act_maxs=act_maxs)
    layers = network.nnet.layers

    sz = length(layers[end].bias)
    push!(layers, mean_layer(sz))

    # Compensation for reward shaping
    if s != 0.0
        for i in eachindex(layers)
            layer = layers[i]
            W = layer.weights
            if i > 1
                W = hcat(W, zeros(Float32, size(W, 1), 1))
            end
            W = vcat(W, zeros(Float32, 1, size(W, 2)))
            W[end, i == 1 ? 1 : end] = 1.0f0
            b = vcat(layer.bias, 0.0f0)
            layers[i] = Layer(W, b, layer.activation)
        end

        W = [1.0f0 -Float32(s)]
        b = [0.0f0]
        σ = NeuralVerification.Id()
        push!(layers, Layer(W, b, σ))
        network.scale = s
    end

    return network
end

"""
Convert actor-critic object to extended neural network representing spread of ensemble
(mean absolute deviation).
"""
function spread_network(ac::GlobalResult; act_mins::Vector{Float64}, act_maxs::Vector{Float64})
    network = values_network(ac; act_mins=act_mins, act_maxs=act_maxs)
    layers = network.nnet.layers

    # Distance to mean
    sz = length(layers[end].bias)
    M1 = Matrix{Float32}(I, sz, sz)
    M2 = Float32.(ones(sz, sz) / sz)
    W = M1 - M2
    b = zeros(Float32, sz)
    σ = NeuralVerification.Id()
    push!(layers, Layer(W, b, σ))

    push!(layers, abs_layers(sz)...)
    push!(layers, mean_layer(sz))
    return network
end

"""
    CrossSection

Cross-section type, defining linear mapping from reduced input space to network observation
space. Input can be specified as explicit matrix and bias or as functional mapping.

# Example
```
cs1 = CrossSection([1 -1; 0 1; 0 0], [0, 0, 5])
cs2 = CrossSection((x1, x2) -> (x1, -pi/2, -x2, pi/2, x1 + x2))
```
"""
mutable struct CrossSection
    W::Matrix{<:Real}
    b::Vector{<:Real}
end

function CrossSection(f::Function)
    unit = (n, i) -> [Float64(i == j) for j in 1:n]
    n = first(methods(f)).nargs - 1
    b = collect(f(zeros(n)...))
    W = reduce(hcat, collect(f(unit(n, i)...)) - b for i in 1:n)
    return CrossSection(W, b)
end

"""
    cross_section(network::ExtendedNetwork, cs::CrossSection, limits::NTuple{2, Vector{Float64}})

Produce cross-section of neural network to reduce input dimensionality.
"""
function cross_section(network::ExtendedNetwork, cs::CrossSection, limits::NTuple{2, Vector{Float64}})
    σ = NeuralVerification.Id()
    layer = [Layer(Float32.(cs.W), Float32.(cs.b), σ)]

    # Determine necessary step
    min_obs = compute_output(Network(layer), limits[1])
    max_obs = compute_output(Network(layer), limits[2])
    steps = -min.(min_obs, max_obs)

    # Duplicate and step network
    n_obs = length(cs.b)
    nnet = deepcopy(network.nnet)
    i1, i2 = network.step_up, network.step_down
    nnet.layers[i1].bias[1:n_obs] = Float32.(steps)
    nnet.layers[i2].bias[1:n_obs] = -Float32.(steps)

    # Add cross-section layer
    prepend!(nnet.layers, layer)

    return coalesce(nnet)
end


text:
function createtestvideo(filename; w = 300, h = 300, fps = 25, totalseconds = 10)
    totalframes = fps*totalseconds
    Δt = Millisecond(Second(1))/fps
    framenumber = Observable(1)
    msg = lift(framenumber) do fn
        t = Time(0) + (fn - 1)*Δt
        """
        Test video
        Frame #: $fn
        Time: $t"""
    end
    scene = Scene(resolution = (w, h))#, show_axis=false, backgroundcolor = :black)
    text!(scene, msg)#, color = :white)
    # scene
    record(scene, filename, 1:totalframes, framerate = fps) do i
        framenumber[] = i
    end
end


text:
@testset "indexing" begin
    i = Indices([1,2,3,4,5])

    @test issetequal(getindices(i, i)::Indices, i)
    @test view(i, i) === i

    i2 = Indices([1,2,3])

    @test issetequal(getindices(i, i2)::Indices, i2)
    @test view(i, i2) === i2

    d = zeros(Int64, i)::Dictionary

    @test isequal(getindices(d, i)::Dictionary, d)
    @test issetequal(pairs(getindices(d, i2)::Dictionary), [1=>0, 2=>0, 3=>0])
    @test issetequal(pairs(view(d, i2)::DictionaryView), [1=>0, 2=>0, 3=>0])

    setindices!(d, 1, i2)

    @test issetequal(pairs(d), [1=>1, 2=>1, 3=>1, 4=>0, 5=>0])

    d2 = Dictionary([3,4,2], [2,3,4])

    @test issetequal(pairs(getindices(i, d2)::Dictionary), [4=>3, 2=>4, 3=>2])
    @test issetequal(pairs(view(i, d2)::DictionaryView), [4=>3, 2=>4, 3=>2])

    @test issetequal(pairs(getindices(d, d2)::Dictionary), [2=>0, 3=>1, 4=>1])
    @test issetequal(pairs(view(d, d2)::DictionaryView), [2=>0, 3=>1, 4=>1])

    setindices!(d, 2, d2)

    @test issetequal(pairs(d), [1=>1, 2=>2, 3=>2, 4=>2, 5=>0])
end


text:
@system Water begin
    #FIXME check unit (w.r.t photosynthesis.ET_supply)
    #FIXME check name (probably not the same as photosynthesis.ET_supply)
    # ET_supply?
    water_supply => 0 ~ track(u"g/hr")
end


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 01-2018
#--------------------------------------------------------------------

function testoperator(Nx::Int, Ny::Int)::Array{Float64,4}
    """
    The main operator we're using
    sum(chebw(i,Nx)*chebw(j,Ny)*
         chebd(i,m,Nx)*delta(j,n)*phi[m,n]*
         delta(i,p)*chebd(j,q,Ny)*phi[p,q]
         for m in 1:Nx+1, n in 1:Ny+1, p in 1:Nx+1, q in 1:Ny+1,
         i in 1:Nx+1, j in 1:Ny+1)
    where we sum over i, j in the computation of the operator.

    And for computing the action of the operator on the vector field
    operator[m,n,p,q] phi[m,n] phi[p,q]
    and when you reshape it into a 2D array (m,n) and (p,q) are bunched
    together.
    """

    operator = zeros(Nx+1, Ny+1, Nx+1, Ny+1)
    for index in CartesianRange(size(operator))
        (m,n,p,q) = index.I
        operator[m,n,p,q] = chebw(p,Nx)*chebw(n,Ny)*chebd(p,m,Nx)*chebd(n,q,Ny)
    end
    return operator
end

function testboundaryOP(Nx::Int, Ny::Int)::Array{Float64,2}
    bnd = zeros(Nx+1, Ny+1)
    bnd[1,:] = bnd[:,1] = 1
    return diagm(vec(bnd))
end

function testRHS(fn::Function, Nx::Int, Ny::Int, M::Int, loc::Array{Int,1})
    patch = projectonPatchbyRestriction(fn, Nx, Ny, M, loc)
    for index in CartesianRange(size(patch))
        i = index.I[1]
        j = index.I[2]
        patch[index] = chebw(i,Nx)*chebw(j,Ny)*patch[index]
    end
    return patch
end

@test_broken 2*testoperator(2, 2) ≈ derivOP(2, 2)
@test_broken 2*testoperator(2, 4) ≈ derivOP(2, 4)
@test testboundaryOP(2,2) == shapeH2L(boundaryOP(2,2))
@test testboundaryOP(2,4) == shapeH2L(boundaryOP(2,4))
@test RHS((x,y)->x^8+y^9, 8, 4, 7, [2,1]) ≈ testRHS((x,y)->x^8+y^9, 8, 4, 7, [2,1])


text:
using Test
using StaticArrays
using PolygonArea
using PolygonArea: PolarHalfPlane, Reunion, Point, vertices, complement, rotate, translate, scale, square

@testset "Polygons" begin

    @testset "creation and contained points" begin
        unit_square = square((0.0, 0.0), 1.0)
        @test Point(0.5, 0.5) in unit_square
        @test Point(1.0, 1.0) in unit_square
        @test !(Point(1.2, 1.0) in unit_square)
        @test Point(0.0, 0.0) in vertices(unit_square)
        @test Point(1.0, 0.0) in vertices(unit_square)
        @test !(Point(0.5, 0.5) in vertices(unit_square))

        @test unit_square == rectangle(bottom_left=(0.0, 0.0), top_right=(1.0, 1.0))
        @test unit_square == rectangle(bottom_left=(0.0, 0.0), top_right=(1.0, 1.0))
        @test unit_square == rectangle((0.0, 0.0), (1.0, 1.0))

        @test !(PolygonArea.ConvexPolygon{Int}([]) == square((1, 1), 1))

        rec = PolygonArea.ConvexPolygon([(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0)])
        @test (0.5, 0.5) ∈ rec

        # Should return an error...
        inv_rec = PolygonArea.ConvexPolygon([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)])
        @test !((0.5, 0.5) ∈ inv_rec)
    end

    @testset "properties" begin
        @test isempty(PolygonArea.ConvexPolygon{Int}([]))

        @test PolygonArea.nb_sides(rectangle((0, 0), (2, 2))) == 4

        @test area(PolygonArea.square((0.0, 0.0), 1.0)) == 1.0
        @test area(rectangle((0, 0), (4, 3))) == 12
        @test isapprox(area(circle(Point(0.0, 0.0), 1.0, 100)), π, atol=1e3)
    end

    @testset "transformations" begin
        @test Point(0.0, 1.2) in rotate(rectangle(0.0, 0.0, 1.0, 1.0), π/4)
        @test Point(0.0, 1.2) in translate(rectangle(0.0, 0.0, 1.0, 1.0), SVector(-0.3, 0.5))
        @test Point(0.0, 1.2) in translate(rectangle(0.0, 0.0, 1.0, 1.0), [-0.3, 0.5])
        @test Point(0.0, 1.2) in translate(rectangle(0.0, 0.0, 1.0, 1.0), (-0.3, 0.5))
        @test Point(1.5, 1.5) in scale(rectangle(0.0, 0.0, 1.0, 1.0), 1.5)
        @test !(Point(1.5, 1.5) in scale(rectangle(0.0, 0.0, 1.0, 1.0), 1.5, center=(1.0, 1.0)))
        λ = rand(2)
        @test area(scale(rectangle(0.0, 0.0, 1.0, 1.0), λ)) ≈ prod(λ)*area(rectangle(0.0, 0.0, 1.0, 1.0))
    end

    @testset "union and intersection" begin
        # Intersection with HalfPlane
        unit_square = square((0.0, 0.0), 1.0)
        hp = HalfPlane(1.0, 1.0, -1.5)
        c1 = unit_square ∩ hp  # Cut the top-right corner
        @test Point(0.5, 0.5) in c1
        @test !(Point(0.9, 0.9) in c1)

        c2 = unit_square ∩ complement(hp)  # Keep only the top-right corner
        @test !(Point(0.5, 0.5) in c2)
        @test Point(0.9, 0.9) in c2

        @test area(c1) + area(c2) ≈ area(unit_square)

        # Intersection with union and intersection of HalfPlanes
        left = HalfPlane(1.0, 0.0, -0.5)
        top = HalfPlane(0.0, 1.0, -0.5)
        top_left = top ∩ left
        bottom_right = complement(top) ∩ complement(left)
        @test bottom_right == complement(top ∪ left)
        @test (unit_square ∩ left) ∩ top == unit_square ∩ (left ∩ top)

        @test area(unit_square ∩ left) == 0.5
        @test area((unit_square ∩ left) ∩ top) == 0.25
        @test area(unit_square ∩ (top ∪ left)) == 0.75

        @test area(unit_square ∩ (bottom_right ∪ top_left)) ≈ 0.5
        @test area(unit_square ∩ (bottom_right ∩ top_left)) ≈ 0.0

        @test (square((0.0, 0.0), 1.0) |> complement |> typeof) == Reunion{HalfPlane{Float64}}
        @test area(square((0.0, 0.0), 3.0) \ square((1.0, 1.0), 1.0)) == 8.0

        # Intersection with ConvexPolygon
        @test area(unit_square ∩ square((0.0, 0.5), 1.0)) == 0.5
        @test area(unit_square ∩ square((0.5, 0.5), 1.0)) == 0.25

        triangle = ConvexPolygon([(0.0, 0.0), (0.5, 1.0), (1.0, 0.0)])
        @test area(unit_square ∩ triangle) == area(triangle)

        # Intersection with union of ConvexPolygon
        r1 = rectangle(0.0, 0.0, 2.0, 2.0)
        r2 = rectangle(1.0, 1.0, 2.0, 2.0)
        r3 = rectangle(1.5, 1.5, 2.5, 2.5)
        @test area((unit_square ∪ r2) ∩ (r1 ∪ r3)) == 2.0

        # Intersection with circle
        c = circle(0.0, 0.0, 1.0, 100)
        r = rectangle(0.0, -1.0, 1.0, 1.0)
        @test area(c ∩ r) ≈ area(r ∩ c) ≈ area(c)/2

        # Disjoint union
        c = circle(0.0, 0.0, 1.0, 10)
        @test PolygonArea.disjoint(c) == c

        h = rectangle(0.0, 1/3, 1.0, 2/3)
        v = rectangle(1/3, 0.0, 2/3, 1.0)
        @test area(h ∪ v) ≈ 6/9
        @test area(PolygonArea.disjoint(h ∪ v)) ≈ 5/9

        # A floating point bug: intersection of a circle with a tangeant half plane
        c = circle((0.5, 0.5), 0.4, 60) |> PolygonArea.rotate(π/5, center=(0.5, 0.5))
        r = rectangle(0.45, 0.10, 0.46, 0.11)
        @test PolygonArea.area(r ∩ c) ≈ PolygonArea.area(c ∩ r)

        #
        horizontal_rectangle = PolygonArea.rectangle((0.4, 0.2), (0.6, 0.8))
        vertical_rectangle = PolygonArea.rectangle((0.2, 0.4), (0.8, 0.6))
        cross = horizontal_rectangle ∪ vertical_rectangle
        circ = PolygonArea.circle(0.5, 0.5, 0.35, 100)
        diff = circ \ cross
        @test area(diff) ≈ area(circ) - area(PolygonArea.disjoint(cross))
    end

    @testset "generic programming" begin
        using Unitful: m
        @test area(PolygonArea.ConvexPolygon([(1m, 0m), (1m, 1m), (0m, 1m), (0m, 0m)])) == 1m^2

        # Conversion and promotion
        T = Reunion{Intersection{HalfPlane{Float64}}}
        @test convert(T, HalfPlane(1f0, 1f0, 0f0)) isa T
        @test convert(T, rectangle(1f0, 1f0, 2f0, 2f0)) isa T

        @test (HalfPlane(1, 1, 1) ∩ HalfPlane(1f0, 2f0, 3f0)) isa Intersection{HalfPlane{Float32}}
        @test (rectangle(0, 0, 2, 2) ∩ rectangle(1f0, 1f0, 3f0, 3f0)) isa ConvexPolygon{Float32}
    end
end


text:
using AnalogDiscovery

function go()
  cDev = enumDevices()
  @show ("Available",cDev)
  myDev = Int32(0)
  if cDev > 0
    if enumDeviceIsOpened(myDev)
      error("Error opening device, already in use.")
    end
    hdwf = deviceOpen(myDev)
  else
    error("No devices available.")
  end

  digitalIOOutputEnableSet(hdwf, UInt32(0x0001))
  # Toggle DigitalIO pin 0 as fast as possible
  for i in 1:100000
    digitalIOOutputSet(hdwf, UInt32(0x01))
    digitalIOOutputSet(hdwf, UInt32(0x00))
  end

  deviceClose(hdwf)
end

go()


text:
"""
    glu(x, dim = 1)

The gated linear unit from the ["Language Modeling with Gated Convolutional Networks"](https://arxiv.org/abs/1612.08083) paper.

Calculates `a .* sigmoid(b)`, where `x` is split in half along given dimension `dim` to form `a` and `b`.
"""
function glu(x, dim = 1)
    maxdim = size(x, dim)
    @assert maxdim % 2 == 0 "Dimension must be even"
    half = maxdim ÷ 2
    a, b = selectdim(x, dim, 1:half), selectdim(x, dim, half+1:maxdim)
    a .* sigmoid.(b)
end


text:
### A Pluto.jl notebook ###
# v0.16.4

using Markdown
using InteractiveUtils

# ╔═╡ 4af33954-3bc9-11ec-0799-7570e35dcd61
begin
	import Pkg
	Pkg.activate(joinpath(@__DIR__))
	using SimLib, Plots, PlutoUI, Statistics, LinearAlgebra, LaTeXStrings, XXZNumerics, SpinSymmetry, LsqFit, SparseArrays
	gr()
	LinearAlgebra.BLAS.set_num_threads(8)
	#cd(joinpath(@__DIR__, ".."))
end

# ╔═╡ 46425991-0870-4b28-92c4-45307b1567c1
path = abspath(joinpath(@__DIR__, "../data/zero-field"))

# ╔═╡ 90d15710-4cb2-4021-978c-41e7fe33cbee
entdata = load_entropy(:box_pbc, 1, 12, 6, 5; prefix=path);

# ╔═╡ f36d56c2-6d96-4c84-b006-b35b2c5192aa
count(isnan, entdata.data), length(entdata.data)

# ╔═╡ 5c9c8eec-10d8-43bb-b680-0ff0704d34a8
plot(entdata.ρs, meandrop(entdata.data; dims=(1,2,3,4)))

# ╔═╡ be3b69f9-cdf3-4f7c-bd76-9fc3a9803385
posdata = load_positions(:box_pbc, 1, 12; prefix=path);

# ╔═╡ 5aefbe53-096a-41ef-8e32-83ee37be92ea
edd = EDDataDescriptor(posdata, 6, [0], :ensemble, symmetrized_basis(zbasis(12,5)))

# ╔═╡ 8c2b33f3-44d6-499c-90c7-377fa96eb30d
task = HalfChainEntropyZBlock(6)

# ╔═╡ 43e48e73-7f33-4a58-909c-631fc4ae9750
ED.initialize!(task, edd, ED._array_constructor);

# ╔═╡ f424dd93-d549-4687-bb29-262862fa14a4
task.data

# ╔═╡ 9de36247-b784-4bf9-8899-91e91808ea2f
ED.compute_task!(task, 1,1,1,[], rand(792,792))

# ╔═╡ Cell order:
# ╠═4af33954-3bc9-11ec-0799-7570e35dcd61
# ╠═46425991-0870-4b28-92c4-45307b1567c1
# ╠═90d15710-4cb2-4021-978c-41e7fe33cbee
# ╠═f36d56c2-6d96-4c84-b006-b35b2c5192aa
# ╠═5c9c8eec-10d8-43bb-b680-0ff0704d34a8
# ╠═be3b69f9-cdf3-4f7c-bd76-9fc3a9803385
# ╠═5aefbe53-096a-41ef-8e32-83ee37be92ea
# ╠═8c2b33f3-44d6-499c-90c7-377fa96eb30d
# ╠═43e48e73-7f33-4a58-909c-631fc4ae9750
# ╠═f424dd93-d549-4687-bb29-262862fa14a4
# ╠═9de36247-b784-4bf9-8899-91e91808ea2f