text:
function args_lesson(s)
  @add_arg_table s begin
    "--model-type", "-m"
      help = "model to be used - one of 'unetlike', 'multiscale3', 'fcrna'"
      arg_type = String
      default = "multiscale3"
    "--patch-mode", "-n"
      help = "number of patches used for training"
      default = nothing
    "--patch-size", "-s"
      help = "size of the patches used for training"
      arg_type = Int
      default = -1
    "--batch-size", "-b"
      help = "number of patches per mini-batch during training"
      arg_type = Int
      default = -1
    "--batch-normalization", "-B"
      help = "whether to activate batch normalization (0 or 1)"
      arg_type = Int
      default = -1
    "--directory", "-d"
      help = "Overwrite image/label folder suggested by the lesson file"
      arg_type = String
      default = ""
    "--optimizer", "-o"
      help = "Overwrite optimizer suggested by the lesson file"
      arg_type = String
      default = ""
    "--epochs", "-e"
      help = "number of training epochs"
      arg_type = Int
      default = -1
    "--learning-rate", "-l"
      help = "initial learning rate for the chosen optimizer [-o]"
      arg_type = Float64
      default = -1.
    "--proximity-kernel-size", "-k"
      help = "size of the peaks in the proximity map"
      arg_type = Int
      default = -1
    "--proximity-kernel-height", "-K"
      help = "height of the peaks in the proximity map"
      arg_type = Float64
      default = -1.
    "--imageop", "-I"
      help = "image operation to be performed on the patches"
      arg_type = String
      default = ""
    "--no-gpu"
      help = "prevent usage of gpu acceleration even if gpu support is detected"
      action = :store_true
    "lesson"
      help = "DCellC lesson (.dcct) file used for training"
      required = true
      arg_type = String
    "modelfile"
      help = "DCellC model (.dccm) file that the trained model is saved to"
      required = true
      arg_type = String
  end
end

function cmd_lesson(args)
  lesson = lessonload(args["lesson"])

  if args["directory"] != ""
    lesson.folder = args["directory"]
  end
  if args["optimizer"] != ""
    lesson.optimizer = args["optimizer"]
  end
  if args["epochs"] != -1
    lesson.epochs = args["epochs"]
  end
  if args["learning-rate"] > 0
    lesson.lr = args["learning-rate"]
  end
  if args["patch-mode"] != nothing
    lesson.patchmode = parse(Int, args["patch-mode"])
  end
  if args["patch-size"] != -1
    lesson.patchsize = args["patch-size"]
  end
  if args["batch-size"] != -1
    lesson.batchsize = args["batch-mode"]
  end
  if args["batch-normalization"] != -1
    lesson.batchnorm = Bool(args["batch-normalization"])
  end
  if args["proximity-kernel-size"] != -1
    lesson.kernelsize = args["proximity-kernel-size"]
  end
  if args["proximity-kernel-height"] > 0
    lesson.kernelheight = args["proximity-kernel-height"]
  end
  if args["imageop"] != ""
    lesson.imageop = eval(parse(args["imageop"]))
  end

  if args["no-gpu"] || gpu() < 0
    println("# CPU mode")
    at = Array{Float32}
  else
    println("# GPU mode")
    at = KnetArray{Float32}
  end

  train(lesson, modelpath = args["modelfile"], record = false)
end

max_stars_repo_name: tscode/DCellC

max_issues_repo_name: tscode/DCellC


text:
module Regional

using DrWatson
using Turing
using StatsBase
using MCMCChains: Chains, meanplot, gelmandiag
using PrettyTables
using NamedTupleTools
using BSON
# using Parameters
# using ArgCheck
# using RData
using CSV
using DataFrames
using OrderedCollections
using Dates
using StatsBase
using StatsFuns
using Base:@kwdef
using Base.Threads# using FillArrays
# using KernelFunctions
using Random
using StatsPlots
using Plots.PlotMeasures
using LaTeXStrings
using ArviZ
import PyPlot as plt
using Underscores
using WeakRefStrings
using InlineStrings

# using QuadGK
# using

using ..CovidSurvey#: NegativeBinomial2

include("models.jl")
include("utils.jl")
include("data.jl")
include("visualization.jl")

# const model = model_v1 # <= defines the "official" model for this sub-module

end

ext: jl

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
@testset "Deconstruct" begin
    let p = deconstruct(f)
        @test p[:category] == :function
        @test p[:name] == :f
        @test p[:args] == [:x, Expr(:kw, :y, 3)]
        @test p[:kwargs] == [:a, Expr(:kw, :b, 5)]
        @test p[:short]
        @test issubset([50], p[:body].args)
    end

    @test deconstruct(doc_f)[:docstring] == "I document function f"

    let p = deconstruct(f2)
        @test p[:category] == :function
        @test p[:name] == :f
        @test p[:args] == [:x]
        @test issubset([:(println(2)), :(return 5)], p[:body].args)
        @test !p[:short]
    end

    let p = deconstruct(s)
        @test p[:category] == :struct
        @test p[:name] == :HandleType
        @test p[:fields] == [:(handle::Ptr{Nothing})]
        @test p[:is_mutable] == true
        @test p[:constructors] == []
    end

    @test deconstruct(doc_s)[:docstring] == "I document structure s"

    let p = deconstruct(c)
        @test p[:category] == :const
        @test p[:value] == :val
    end

    let p = deconstruct(e)
        @test p[:category] == :enum
        @test p[:name] == :myenum
        @test p[:macro] == Symbol("@enum")
        @test p[:values] == [:a, :b, :c, :d]
    end

    let p = deconstruct(e2)
        @test p[:category] == :enum
        @test p[:name] == :myotherenum
        @test p[:macro] == Symbol("@cenum")
        @test p[:values] == [:(a = 1), :(b = 2), :(c = 3), :(d = 4)]
    end
end


text:
@testset "735.asteroid-collision.jl" begin
    @test asteroid_collision([5, 10, -5]) == [5, 10]
    @test asteroid_collision([8, -8]) == Int[]
    @test asteroid_collision([10, 2, -5]) == [10]
    @test asteroid_collision([-2, -1, 1, 2]) == [-2, -1, 1, 2]
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ dc111222-038e-479b-993e-5ed58755df27
    begin
       using Pkg
       Pkg.activate(".testenv")
	   Pkg.add("Revise")
	   using Revise
       Pkg.add(["PlutoUI","Triangulate","LaTeXStrings"])
	   Pkg.develop("PlutoVista")
     end

# ╔═╡ d6c0fb79-4129-444a-978a-bd2222b53df6
begin
    using PlutoUI
    using PlutoVista
    using Printf
	using LaTeXStrings
    import Triangulate
end

# ╔═╡ 93ca4fd0-8f61-4174-b459-55f5395c0f56
md"""
# Test Notebook for [PlutoVista](https://github.com/j-fu/PlutoVista.jl) and plotly.js
"""

# ╔═╡ 9072dcea-e634-493e-ba1a-890220737683
let X=collect(0:0.1:10); plot(X,sin.(X)) end

# ╔═╡ 763bdefc-200e-4afb-a6c9-71fb5ee16f58
let X=collect(0:0.1:10); plot(X,X.^2,xscale=:log,label="x", legend=:lt,legendfontsize=20) end

# ╔═╡ 4cd2182a-ea7c-49f6-a3f8-029b8727430f
let X=collect(0:0.1:100); plot(X,X.^2,yscale=:log,label="x",legend=:rb,tickfontsize=20) end

# ╔═╡ b9e1184b-62db-41c9-9340-b9db0aff5b78
 let X=collect(0:0.1:100);
	p=PlutoVistaPlot(xscale=:log,yscale=:log,ylabel="x",xlabel="ψ",title="title",legend=:lt,axisfontsize=20)
	plot!(p,X,X.^2;color=:red,label="2")
	plot!(p,X,X.^3,label="3")
	plot!(p,X,X.^4,label="4")
end

# ╔═╡ ad1da8b8-723f-4b18-ba50-fb5a6d5d1176
 let X=collect(0:0.1:1000);
	p=PlutoVistaPlot(dim=1,resolution=(600,300),title="title",titlefontsize=20,legend=:lt,xscale=:log,yscale=:log)
	plot!(p,X,X.^20,color=:red,label="2")
	plot!(p,X,X.^30,label="3")
plot!(p,X,X.^40,ylabel="y",xlabel="xxx",label="4")
end

# ╔═╡ 7c06fcf0-8c98-49f7-add8-435f57a9c9da
function maketriangulation(maxarea)

    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([-1.0 -1.0 ; 1.0 -1.0 ; 1.5  1.5 ; -1.0 1.0]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ; 3 4 ; 4 1 ]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4])
    area=@sprintf("%.15f",maxarea)
    (triout, vorout)=Triangulate.triangulate("pa$(area)DQ", triin)
    triout.pointlist, triout.trianglelist
end

# ╔═╡ db2823d9-aa6d-4be3-af5c-873c072cfd2b
md"""
Change grid resolution: $(@bind resolution Slider(20:200))
"""

# ╔═╡ 890710fe-dac0-4256-b1ba-79776f1ea7e5
(pts,tris)=maketriangulation(1/resolution^2)

# ╔═╡ b8a976e3-7fef-4527-ae6a-4da31c93a04f
func=0.5*[sin(10*pts[1,i])*cos(10*pts[2,i])*10*pts[1,i] for i=1:size(pts,2)]

# ╔═╡ 60dcfcf5-391e-418f-8e7c-3a0fe94f1e0d
p=let
	p=PlutoPlotlyPlot(resolution=(500,300))
	triplot!(p,pts,tris,func)
end

# ╔═╡ 401b36bd-fa8f-4a9c-9556-bbc82c3ddbca
md"""
Change time: $(@bind time Slider(0:0.1:10,show_value=true))
"""

# ╔═╡ e76f8a6a-ab91-454a-b200-cfc8b57eb331
triupdate!(p,pts,tris,0.5*[sin(10*pts[1,i]-time)*cos(10*pts[2,i]-time) for i=1:size(pts,2)])

# ╔═╡ bce0cfe7-4112-4bb8-aac6-43885f3746a9
md"""Number of gridpoints: $(size(pts,2)) """

# ╔═╡ 81046dcd-3cfb-4133-943f-61b9b3cdb183
let
	p=PlutoPlotlyPlot(resolution=(700,300))
	tricontour!(p,pts,tris,func,levels=11,colormap=:summer)
end

# ╔═╡ e900801e-2020-4aff-bfec-017ad6fcfdcf
md"""
shift: $(@bind shift Slider(-1:0.01:1,show_value=true))
"""

# ╔═╡ 5cb73674-c144-4fd9-8ff6-ac767548822e
f(x,y)=sin(3(x-shift))*cos(3(y-shift))*10

# ╔═╡ 34072263-1180-4bc0-a004-2f112ea4cbed
let
	X1=-2:0.1:2
	Y1=-1:0.1:1
	ff=[f(x,y) for x∈ X1, y∈ Y1]
	p=PlutoPlotlyPlot(resolution=(500,300))
	contour!(p,collect(X1),collect(Y1),ff,colormap=:hot)
end

# ╔═╡ 5c23809b-9d92-43da-a85d-6c8531c3b547
let X=collect(0:0.1:10); plot(X,sin.(X)) end

# ╔═╡ edfe1dd2-94a5-403f-b589-53bfca839057
"\$x^k\$"

# ╔═╡ 9709e98a-50f3-4ee5-912c-e370c8c93193
 function fpdens(x::AbstractFloat;sample_size=1000)
    xleft=x
    xright=x
    for i=1:sample_size
        xleft=prevfloat(xleft)
        xright=nextfloat(xright)
    end
    return prevfloat(2.0*sample_size/(xright-xleft))
end;

# ╔═╡ f8684eae-00c8-4e19-b638-5be32c045cfc
X=10.0.^collect(-40:0.1:40);

# ╔═╡ 7f5612f3-a0cb-4319-b4bf-4eefe5dc2a18
let
    p=plot(resolution=(600,300),
		title="Number of numbers per unit interval",
		xscale=:log,yscale=:log,xlabel="x",ylabel="n",legend=:rt)
    plot!(p,X,map(x->1.0e20/x,X), label="O(1/x)",linewidth=0.5,color=:black,linestyle=:dot)
    plot!(p,X,fpdens.(Float16.(X)),label="Float16")
    plot!(p,X,fpdens.(Float32.(X)),label="Float32")
    plot!(p,X,fpdens.(Float64.(X)),label="Float64")
    plot!(p,X,map(x-> x<typemax(Int8) ? 1  : 0,X),
		linestyle=:dash,label="Int8")
    plot!(p,X,map(x-> x<typemax(Int32) ? 1 : 0,X),
		linestyle=:dash,label="Int32")
    plot!(p,X,map(x-> x<typemax(Int64) ? 1 : 0,X),
		linestyle=:dash,label="Int64")
    plot!(p,X,map(x-> x<typemax(Int16) ? 1 : 0,X),
		linestyle=:dash,color=:black,label="Int16")
end

# ╔═╡ 88fd8e90-e054-4dd4-858d-7367b590b89d
html"""<hr>"""

# ╔═╡ Cell order:
# ╠═93ca4fd0-8f61-4174-b459-55f5395c0f56
# ╠═d6c0fb79-4129-444a-978a-bd2222b53df6
# ╠═9072dcea-e634-493e-ba1a-890220737683
# ╠═763bdefc-200e-4afb-a6c9-71fb5ee16f58
# ╠═4cd2182a-ea7c-49f6-a3f8-029b8727430f
# ╠═b9e1184b-62db-41c9-9340-b9db0aff5b78
# ╠═ad1da8b8-723f-4b18-ba50-fb5a6d5d1176
# ╠═7c06fcf0-8c98-49f7-add8-435f57a9c9da
# ╠═890710fe-dac0-4256-b1ba-79776f1ea7e5
# ╠═b8a976e3-7fef-4527-ae6a-4da31c93a04f
# ╠═60dcfcf5-391e-418f-8e7c-3a0fe94f1e0d
# ╟─db2823d9-aa6d-4be3-af5c-873c072cfd2b
# ╟─401b36bd-fa8f-4a9c-9556-bbc82c3ddbca
# ╠═e76f8a6a-ab91-454a-b200-cfc8b57eb331
# ╟─bce0cfe7-4112-4bb8-aac6-43885f3746a9
# ╠═81046dcd-3cfb-4133-943f-61b9b3cdb183
# ╟─e900801e-2020-4aff-bfec-017ad6fcfdcf
# ╠═5cb73674-c144-4fd9-8ff6-ac767548822e
# ╠═34072263-1180-4bc0-a004-2f112ea4cbed
# ╠═5c23809b-9d92-43da-a85d-6c8531c3b547
# ╠═edfe1dd2-94a5-403f-b589-53bfca839057
# ╠═9709e98a-50f3-4ee5-912c-e370c8c93193
# ╠═f8684eae-00c8-4e19-b638-5be32c045cfc
# ╠═7f5612f3-a0cb-4319-b4bf-4eefe5dc2a18
# ╟─88fd8e90-e054-4dd4-858d-7367b590b89d
# ╠═dc111222-038e-479b-993e-5ed58755df27


text:
module Sigma

deps_dir = joinpath(dirname(@__FILE__),"..","deps")
prefix = joinpath(deps_dir,"usr")
src_dir = joinpath(prefix,"src")
bin_dir = joinpath(prefix,"bin")
lib_dir = joinpath(prefix,"lib")

using DReal
using Distributions
using AbstractDomains
using Lens
using DataStructures
# using Compat

import Base: ifelse, cond, isequal, isinf
import Base: sqrt, abs, promote_rule, convert, rand, getindex, string, size
import Base: show, print, showcompact
import Base: sum, dot, length, join, round
import Base: ndims, endof
import Base.isapprox
import Base: start, next, done
import Base: hash
import Base: ndims, isequal, union, push!, string, print, show, println
import Base.eltype
import Base.size
import Base.all
import Base: one, zero, norm, similar
import Base: real
import Base: var
import Base: in
import Base: collect
import Base: setindex!,
             hcat,
             vcat,
             .^,
             .*,
             .>,
             .>=,
             .<,
             .<=,
             .==,
             .!=,
             transpose,
             ctranspose,
             ==,
             ^,
             +,
             -,
             *,
             /,
             >,
             >=,
             <=,
             <,
             <,
             <,
             !=,
             |,
             &,
             !

import Base:  asin,
              sqrt,
              exp,
              log,
              cos,
              sin,
              tan,
              acos,
              asin,
              atan,
              cosh,
              sinh,
              tanh,
              acosh,
              asinh,
              atanh,
              abs,
              atan2,
              max,
              min,
              sign

# import Base: exponential
juliadir = joinpath(homedir(),".julia","v0.4")

# import Lens:benchmark

import AbstractDomains: dims, Interval, Boxes
import Distributions: quantile
import DReal: model, set_precision!

export
  # Random Variables
  RandVar,
  RandArray,
  dims,

  # Abstract Domains
  rangetype,

  LA,
  Lifted,
  liftedarray,
  LiftedArray,

  # Inference queries
  prob,
  sample_mean,

  setindex,
  model,

  # RandVar Types
  RandVar,
  AllRandVars,

  # Distributions
  normal,
  uniform,
  flip,
  exponential,
  logistic,
  geometric,
  discreteuniform,
  iid,

  dirichlet,
  mvnormal,
  mvuniform,
  mvexponential,
  mvlogistic,

  @noexpand,

  #utils
  sqr,
  ≊,
  pnormalize,

  #Plotting
  loadvis,
  plot_2d_boxes,
  plot_density,
  plot_cond_density,
  plot_volume_distribution,
  plot_performance,
  plot_sat_distribution,
  distinguished_colors,
  rand_color,
  plot_sample_cond_density,
  plot_sample_density,

  # RandVars
  DRealRandVar,

  # From abstract domains
  mid,

  isapprox,
  quantile

include("util.jl")
include("domains.jl")
include("omega.jl")
include("precision.jl")
include("randvar.jl")
include("interop.jl")
include("refinement/refinement.jl")
include("query.jl")
include("distributions.jl")
include("split.jl")
include("pmaplm.jl")
include("show.jl")

# Hack to avoid loading Gadfly each time
vispath = joinpath(juliadir, "Sigma","src","vis.jl")

loadvis() = include(vispath)

"""
A library for probabilistic programming.

API Overview
- randvar.jl contains random variable types
  `SymbolicRandVar' is the type used when constructed models, essentially an AST
- `randvars/elementary.jl` containts primitive random variables, uniform, normal
   etc.
- `randvars/symbolic.jl` has functions for constructing complex SymbolicRandVars
   pointwise, by implementing common functions.

"""
Sigma

end


text:
include("basics/basics.jl")
include("integrated_environments/integrated_environments.jl")


text:
# Basic Image(Shepp Logan Phantom) and

function test_kdata(N::Int64=32)
    @info "Testing simulating kdata with NFFT and exact evaluation"
    I = shepp_logan(N)
    tr = RadialTrajectory(N,N)
    println("Simulating kdata using NFFT")
    @time acqDataNFFT = simulation_fast(tr,I)
    @info "Simulating kdata rigorously"
    @time acqDataExplicit = simulation_explicit(tr,I)
    relError =  norm(acqDataExplicit.kdata[:]-acqDataNFFT.kdata[:]) / norm(acqDataExplicit.kdata[:])
    println("Relative error NFFT vs EXACT: ", relError)
    @test relError < 1e-2

end

function test_kdataMultipleSlices(N::Int64=32)
    @info "Testing simulating multiple 2d-slices with NFFT and exact evaluation"
    sh = ComplexF64.(shepp_logan(N))
    I = cat(sh,0.6*sh,0.3*sh,dims=3)
    tr = CartesianTrajectory(N,N)
    println("Simulating kdata using NFFT")
    @time acqDataNFFT = simulation(tr,I,opName="fast")
    @info "Simulating kdata rigorously"
    @time acqDataExplicit = simulation(tr,I,opName="explicit")
    relError =  norm(acqDataExplicit.kdata[:]-acqDataNFFT.kdata[:]) / norm(acqDataExplicit.kdata[:])
    println("Relative error NFFT vs EXACT: ", relError)
    @test relError < 1e-2

end

function test_kdata3d(N::Int64=16)
    @info "Testing simulating 3d-kdata with NFFT and exact evaluation"
    sh = ComplexF64.(shepp_logan(N))
    I = cat(sh,0.9*sh,0.8*sh,0.7*sh,0.6*sh,0.5*sh,0.4*sh,0.3*sh,dims=3)
    tr = CartesianTrajectory3D(N,N,numSlices=8)
    println("Simulating kdata using NFFT")
    @time acqDataNFFT = simulation(tr,I,opName="fast")
    @info "Simulating kdata rigorously"
    @time acqDataExplicit = simulation(tr,I,opName="explicit")
    relError =  norm(acqDataExplicit.kdata[:]-acqDataNFFT.kdata[:]) / norm(acqDataExplicit.kdata[:])
    println("Relative error NFFT vs EXACT: ", relError)
    @test relError < 1e-2

end

function test_kdataWithCorrection(N::Int64=32)
    # Testing generating kdata with fieldinhomogeneity
    I = shepp_logan(N)
    @info "Testing simulating kdata with correctionterm"
    fmap = quadraticFieldmap(N,N,125*2pi)[:,:,1]
    rmap = relaxationMap(I,5.0,50.0)
    cmap = rmap + 1im*fmap
    tr = SpiralTrajectory(N,2*N,TE=0.0,AQ=32e-3)

    println("Simulating kdata using NFFT to approx. correctionterm ...")
    @time acqDataNFFT = simulation_fast(tr,I,cmap;method="nfft")
    println("Simulating kdata using Least Squares to approx. correctionterm...")
    @time acqDataLSQR = simulation_fast(tr,I,cmap;method="leastsquare")
    println("Simulating kdata rigorously...")
    @time acqDataExplicit = simulation_explicit(tr,I,cmap)
    # Calculating and testing the relative error of kdata
    relErrorLeastSquare = norm(acqDataExplicit.kdata[:]-acqDataLSQR.kdata[:]) / norm(acqDataExplicit.kdata[:])
    println("Relative Error of leastsquare method: ", relErrorLeastSquare)
    @test relErrorLeastSquare < 1e-3

    relErrorNFFT = norm(acqDataExplicit.kdata[:] - acqDataNFFT.kdata[:]) / norm(acqDataExplicit.kdata[:])
    println("Relative Error of leastsquare method: ", relErrorNFFT)
    @test relErrorNFFT < 1e-3

end

function test_kdataMultiEcho(N=32)
    # image
    I = ComplexF64.(shepp_logan(N))
    @info "Testing simulating kdata with multiple echoes"
    rmap = 20.0*ones(N,N)

    # simulation
    params = Dict{Symbol, Any}()
    params[:simulation] = "fast"
    params[:trajName] = "Cartesian"
    params[:numProfiles] = floor(Int64, N)
    params[:numSamplingPerProfile] = N
    params[:TE] = 0.0

    acqData1 = simulation( real(I), params )
    kdata1 = acqData1.kdata

    params[:r2map] = rmap
    params[:T_echo] = [2.e-2, 4.e-2]
    params[:seqName] = "ME"
    params[:refocusingAngles] = Float64[pi,pi]

    acqData2 = simulation( real(I), params )
    kdata2 = reshape(acqData2.kdata,:,2)

    relErrorEcho1 = norm(exp(-20.0*2.e-2)*kdata1 - kdata2[:,1])/norm(exp(-20.0*2.e-2)*kdata1)
    @test relErrorEcho1 < 1e-3
    relErrorEcho2 = norm(exp(-20.0*4.e-2)*kdata1 - kdata2[:,2])/norm(exp(-20.0*4.e-2)*kdata1)
    @test relErrorEcho2 < 1e-3
end

function test_noise(N::Int64=32, snr::Float64=25.0)
    @info "Testing simulating kdata with NFFT and exact evaluation"
    I = shepp_logan(N)
    tr = RadialTrajectory(N,N)
    @time acqData = simulation_explicit(tr,I)
    acqDataNoisy = MRIReco.addNoise(acqData,snr)
    relError =  norm(acqData.kdata[:]-acqDataNoisy.kdata[:]) / norm(acqData.kdata[:])
    println("Relative error EXACT vs NOISY: ", relError)
    @test relError < 1e-1
end

function test_changeEncodingSize(N=32)
    @info "Testing reduction of encoding size"
    I = shepp_logan(2*N)
    tr = CartesianTrajectory(2*N,2*N)
    acqData = simulation_fast(tr,I)
    acqData.encodingSize = [2*N,2*N]
    acqData2 = MRIReco.changeEncodingSize2D(acqData, [N,N])
    kdata = reshape(acqData.kdata[1],2*N,2*N)[div(N,2)+1:end-div(N,2), div(N,2)+1:end-div(N,2)]
    kdata2 = reshape(acqData2.kdata[1],N,N)
    relError = norm(kdata[:]-4*kdata2[:]) / norm(kdata)
    @test relError < 1.e-3
end

function testAcqData(N=32)
    @test "testing AcquisitionData"
    # image
    x = shepp_logan(N)
    rmap = 20.0*ones(N,N)

    # coil sensitivites
    sensMaps = zeros(ComplexF64,N*N,2,1)
    sensMaps[1:floor(Int64, N*N/2),1,1] .= 1.0
    sensMaps[floor(Int64, N*N/2)+1:end,2,1] .= 1.0

    # simulation
    params = Dict{Symbol, Any}()
    params[:simulation] = "fast"
    params[:trajName] = "Cartesian"
    params[:numProfiles] = floor(Int64, N)
    params[:senseMaps] = reshape(sensMaps,N,N,1,2)
    params[:numSamplingPerProfile] = N

    #sequence parameters
    params[:r2map] = rmap
    params[:T_echo] = [2.e-2, 4.e-2]
    params[:seqName] = "ME"
    params[:refocusingAngles] = Float64[pi,pi]

    acqData = simulation(x,params)

    @test numContrasts(acqData) = 2
    @test numChannels(acqData) = 2
    @test numSlices(acqData) = 1
    @test numRepetitions(acqData) = 1

end

function testSimulation()
  @testset "simulations" begin
    test_kdata()
    test_kdataMultipleSlices()
    test_kdata3d()
    test_kdataWithCorrection(16)
    test_kdataMultiEcho()
    test_noise()
    test_changeEncodingSize(32)
  end
end

testSimulation()


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1SecretKeySelector <: SwaggerModel
    key::Any # spec type: Union{ Nothing, String } # spec name: key
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    optional::Any # spec type: Union{ Nothing, Bool } # spec name: optional

    function IoK8sApiCoreV1SecretKeySelector(;key=nothing, name=nothing, optional=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1SecretKeySelector, Symbol("key"), key)
        setfield!(o, Symbol("key"), key)
        validate_property(IoK8sApiCoreV1SecretKeySelector, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1SecretKeySelector, Symbol("optional"), optional)
        setfield!(o, Symbol("optional"), optional)
        o
    end
end # type IoK8sApiCoreV1SecretKeySelector

const _property_map_IoK8sApiCoreV1SecretKeySelector = Dict{Symbol,Symbol}(Symbol("key")=>Symbol("key"), Symbol("name")=>Symbol("name"), Symbol("optional")=>Symbol("optional"))
const _property_types_IoK8sApiCoreV1SecretKeySelector = Dict{Symbol,String}(Symbol("key")=>"String", Symbol("name")=>"String", Symbol("optional")=>"Bool")
Base.propertynames(::Type{ IoK8sApiCoreV1SecretKeySelector }) = collect(keys(_property_map_IoK8sApiCoreV1SecretKeySelector))
Swagger.property_type(::Type{ IoK8sApiCoreV1SecretKeySelector }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1SecretKeySelector[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1SecretKeySelector }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1SecretKeySelector[property_name]

function check_required(o::IoK8sApiCoreV1SecretKeySelector)
    (getproperty(o, Symbol("key")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1SecretKeySelector }, name::Symbol, val)
end


text:
function _git_cmd(str::AbstractString;
                  adjust_PATH::Bool = true,
                  adjust_LIBPATH::Bool = true)
    cmd = git(; adjust_PATH, adjust_LIBPATH)
    return `$(cmd) $(split(str))`
end

macro git_cmd(ex)
    return _git_cmd(ex)
end


text:
function checked_mod1(x::T, y::T) where T<:Integer
    result = checked_mod(x, y)
    result = ifelse(result === zero(T), y, result)
    return result
end

function checked_fld1(x::T, y::T) where T<:Integer
    d = checked_div(x, y)
    return d + (!signbit(x ⊻ y) & (d * y !== x))
end

for (OP, CHK) in ((:(+), :checked_add), (:(-), :checked_sub),
                  (:(*), :checked_mul), (:div, :checked_div),
                  (:fld, :checked_fld), (:cld, :checked_cld),
                  (:rem, :checked_rem), (:mod, :checked_mod),
                  (:mod1, :checked_mod1), (:fld1, :checked_fld1))
    @eval begin
       @inline function $OP(x::T, y::T) where T<:SafeSigned
            ix = baseint(x)
            iy = baseint(y)
            result = $CHK(ix, iy)
            return safeint(result)
       end

       @inline function $OP(x::T, y::T) where T<:SafeUnsigned
            ix = baseint(x)
            iy = baseint(y)
            result = $CHK(ix, iy)
            return safeint(result)
       end

       @inline function $OP(x::T1, y::T2) where {T1<:SafeSigned, T2<:SafeSigned}
            xx, yy = promote(x, y)
            ix = baseint(xx)
            iy = baseint(yy)
            result = $CHK(ix, iy)
            return safeint(result)
        end

        @inline function $OP(x::T1, y::T2) where {T1<:SafeUnsigned, T2<:SafeUnsigned}
            xx, yy = promote(x, y)
            ix = baseint(xx)
            iy = baseint(yy)
            result = $CHK(ix, iy)
            return safeint(result)
        end

        @inline function $OP(x::T1, y::T2) where {T1<:SafeSigned, T2<:SafeUnsigned}
            xx, yy = promote(x, y)
            ix = baseint(xx)
            iy = baseint(yy)
            result = $CHK(ix, iy)
            return safeint(result)
        end

        @inline function $OP(x::T1, y::T2) where {T1<:SafeUnsigned, T2<:SafeSigned}
            xx, yy = promote(x, y)
            ix = baseint(xx)
            iy = baseint(yy)
            result = $CHK(ix, iy)
            return safeint(result)
        end

        @inline function $OP(x::T1, y::T2) where {T1<:SafeInteger, T2<:Integer}
            xx, yy = promote(x, y)
            return $OP(xx, yy)
        end

        @inline function $OP(x::T1, y::T2) where {T1<:Integer, T2<:SafeInteger}
            xx, yy = promote(x, y)
            return $OP(xx, yy)
        end

        @inline function $OP(x::T, y::Bool) where {T<:SafeInteger}
            xx, yy = promote(x, y)
            return $OP(xx, yy)
        end

        @inline function $OP(x::Bool, y::T) where {T<:SafeInteger}
            xx, yy = promote(x, y)
            return $OP(xx, yy)
        end
    end
end

function (/)(x::S, y::S) where S<:SafeInteger
    ix = baseint(x)
    iy = baseint(y)
    checked_div(ix, iy)
    result = ix / iy
    return result
end

function (\)(x::S, y::S) where S<:SafeInteger
    ix = baseint(y)
    iy = baseint(x)
    checked_div(iy, ix)
    result = ix / iy
    return result
end

function (/)(x::S1, y::S2) where {S1<:SafeInteger, S2<:SafeInteger}
   xx, yy = promote(x, y)
   return (/)(xx, yy)
end
function (/)(x::S, y::I) where {S<:SafeInteger, I<:Integer}
   xx, yy = promote(x, y)
   return (/)(xx, yy)
end
function (/)(x::I, y::S) where {I<:Integer, S<:SafeInteger}
   xx, yy = promote(x, y)
   return (/)(xx, yy)
end
function (/)(x::S, y::Bool) where {S<:SafeInteger}
   xx, yy = promote(x, y)
   return (/)(xx, yy)
end
function (/)(x::Bool, y::S) where {S<:SafeInteger}
   xx, yy = promote(x, y)
   return (/)(xx, yy)
end

function (\)(x::S1, y::S2) where {S1<:SafeInteger, S2<:SafeInteger}
   xx, yy = promote(x, y)
   return (\)(xx, yy)
end
function (\)(x::S, y::I) where {S<:SafeInteger, I<:Integer}
   xx, yy = promote(x, y)
   return (\)(xx, yy)
end
function (\)(x::I, y::S) where {I<:Integer, S<:SafeInteger}
   xx, yy = promote(x, y)
   return (\)(xx, yy)
end

function divrem(x::S, y::S) where S<:SafeInteger
    ix = baseint(x)
    iy = baseint(y)
    return safeint(div(ix, iy)), safeint(rem(ix, iy)) # div, rem already are checked
end

function divrem(x::S1, y::S2) where {S1<:SafeInteger, S2<:SafeInteger}
   xx, yy = promote(x, y)
   return divrem(xx, yy)
end
function divrem(x::S1, y::S2) where {S1<:SafeInteger, S2<:Integer}
   xx, yy = promote(x, y)
   return divrem(xx, yy)
end
function divrem(x::S1, y::S2) where {S2<:SafeInteger, S1<:Integer}
   xx, yy = promote(x, y)
   return divrem(xx, yy)
end
function divrem(x::S, y::Bool) where {S<:SafeInteger}
   xx, yy = promote(x, y)
   return divrem(xx, yy)
end
function divrem(x::Bool, y::S) where {S<:SafeInteger}
   xx, yy = promote(x, y)
   return divrem(xx, yy)
end

function fldmod(x::S, y::S) where S<:SafeInteger
    ix = baseint(x)
    iy = baseint(y)
    return safeint(fld(ix, iy)), safeint(mod(ix, iy)) # fld, mod already are checked
end

function fldmod(x::S1, y::S2) where {S1<:SafeInteger, S2<:SafeInteger}
   xx, yy = promote(x, y)
   return fldmod(xx, yy)
end
function fldmod(x::S1, y::S2) where {S1<:SafeInteger, S2<:Integer}
   xx, yy = promote(x, y)
   return fldmod(xx, yy)
end
function fldmod(x::S1, y::S2) where {S2<:SafeInteger, S1<:Integer}
   xx, yy = promote(x, y)
   return fldmod(xx, yy)
end
function fldmod(x::S, y::Bool) where {S<:SafeInteger}
   xx, yy = promote(x, y)
   return fldmod(xx, yy)
end
function fldmod(x::Bool, y::S) where {S<:SafeInteger}
   xx, yy = promote(x, y)
   return fldmod(xx, yy)
end

function fldmod1(x::S, y::S) where S<:SafeInteger
    ix = baseint(x)
    iy = baseint(y)
    return safeint(fld1(ix, iy)), safeint(mod1(ix, iy)) # fld1, mod1 already are checked
end

function fldmod1(x::S1, y::S2) where {S1<:SafeInteger, S2<:SafeInteger}
   xx, yy = promote(x, y)
   return fldmod1(xx, yy)
end
function fldmod1(x::S1, y::S2) where {S1<:SafeInteger, S2<:Integer}
   xx, yy = promote(x, y)
   return fldmod1(xx, yy)
end
function fldmod1(x::S1, y::S2) where {S2<:SafeInteger, S1<:Integer}
   xx, yy = promote(x, y)
   return fldmod1(xx, yy)
end
function fldmod1(x::S, y::Bool) where {S<:SafeInteger}
   xx, yy = promote(x, y)
   return fldmod1(xx, yy)
end
function fldmod1(x::Bool, y::S) where {S<:SafeInteger}
   xx, yy = promote(x, y)
   return fldmod1(xx, yy)
end

for F in (:gcd, :lcm)
  @eval begin
    function $F(x::S, y::S) where S<:SafeInteger
        ix = baseint(x)
        iy = baseint(y)
        return safeint($F(ix, iy))
    end

    function $F(x::S1, y::S2) where {S1<:SafeInteger, S2<:SafeInteger}
       xx, yy = promote(x, y)
       return $F(xx, yy)
    end

    function $F(x::S1, y::S2) where {S1<:SafeInteger, S2<:Integer}
       xx, yy = promote(x, y)
       return $F(xx, yy)
    end

    function $F(x::S1, y::S2) where {S2<:SafeInteger, S1<:Integer}
       xx, yy = promote(x, y)
       return $F(xx, yy)
    end

    function $F(x::S, y::Bool) where {S<:SafeInteger}
       xx, yy = promote(x, y)
       return $F(xx, yy)
    end

    function $F(x::Bool, y::S) where {S<:SafeInteger}
       xx, yy = promote(x, y)
       return $F(xx, yy)
    end
  end
end

function divgcd(x::S, y::S) where {S<:SafeInteger}
    g = gcd(x,y)
    return div(x,g), div(y,g)
end

divgcd(x::I, y::S) where {I<:Integer, S<:SafeInteger} = divgcd(promote(x,y)...)
divgcd(x::S, y::I) where {I<:Integer, S<:SafeInteger} = divgcd(promote(x,y)...)
divgcd(x::S1, y::S2) where {S1<:SafeInteger, S2<:SafeInteger} = divgcd(promote(x,y)...)
divgcd(x::S, y::Bool) where {S<:SafeInteger} = divgcd(promote(x,y)...)
divgcd(x::Bool, y::S) where {S<:SafeInteger} = divgcd(promote(x,y)...)


text:
using DynamicalSystemsBase, Test
using LinearAlgebra: norm
@testset "Projected integrator" begin

@testset "Lorenz" begin
    # Set β to 10, we have a stable fixed point
    ds = Systems.lorenz([0.0, 10.0, 1.0]; σ = 10.0, ρ = 28.0, β = 10)
    psys = projected_integrator(ds, 1:2, [0.0])
    @test get_state(psys) == get_state(ds)[1:2]

    reinit!(psys, [1.0, 1.0])
    for j in 1:5
      step!(psys, 1.0)
    end
    y = get_state(psys)
    @test length(y) == 2
    @test abs(sum(y - [16.43, 16.43])) < 0.01

    # Test complete state function
    complete = (y) -> [y[1], y[2], y[2] + 1]
    psys = projected_integrator(ds, 1:2, complete)
    reinit!(psys, [2.0, 1.0])
    y = get_state(psys)
    @test y[1] == 2
    @test length(y) == 2

    for j in 1:5
      step!(psys, 1.0)
    end
    u = get_state(psys)
    @test abs(sum(u - [16.43, 16.43])) < 0.01

    # Test projection function on a sphere of unit radius
    projection(u) = u/norm(u)
    complete = y -> 10y
    psys = projected_integrator(ds, projection, complete)
    @test norm(get_state(psys)) == 1
    reinit!(psys, ones(3))
    @test psys.integ.u[1] == 10

    for j in 1:5
      step!(psys, 1.0)
    end
    u = get_state(psys)
    @test abs(sum(u[1:2] - [0.461, 0.461])) < 0.01

    # Trajectory, continuous
    psys = projected_integrator(ds, projection, complete)
    tr = trajectory(psys, 10)
    @test dimension(tr) == 3
    @test norm(tr[end]) == 1

    # Trajectory, discrete
    ds = Systems.towel()
    psys = projected_integrator(ds, [1,2], [0.0])
    @test get_state(psys) == [0.085, -0.121]
    @test dimension(psys) == 2
    tr = trajectory(psys, 10)
    @test dimension(tr) == 2
end

end


text:
using BenchmarkTools
using JLD
using Test
using MIToS
using MIToS.MSA
using MIToS.Information
using MIToS.PDB
using MIToS.SIFTS
using MIToS.Pfam
using MIToS.Utils

# --------------------------------------------------------------------------- #
# MSA module benchmarks
const msa = BenchmarkGroup()

## Set up
### Files
const msafile_long_sth_gz = "../../data/PF00089.stockholm.gz"
const msafile_long_sth    = "../../data/PF00089.sth"
const msafile_long_fas_gz = "../../data/PF00089.fasta.gz"
const msafile_long_fas    = "../../data/PF00089.fasta"
const msafile_wide_sth_gz = "../../data/PF16957.stockholm.gz"
const msafile_wide_sth    = "../../data/PF16957.sth"
const msafile_wide_fas_gz = "../../data/PF16957.fasta.gz"
const msafile_wide_fas    = "../../data/PF16957.fasta"
### MSAs
const msa_long    = read(msafile_long_sth   , Stockholm)
const msa_wide    = read(msafile_wide_sth   , Stockholm)

#### Parse benchmarks

##### input
msa["input"] = BenchmarkGroup()

for (file,gzipped,shape,format) in [(msafile_long_sth_gz, "gzipped", "long", MIToS.MSA.Stockholm),
                                    (msafile_long_sth,  "ungzipped", "long", MIToS.MSA.Stockholm),
                                    (msafile_long_fas_gz, "gzipped", "long", MIToS.MSA.FASTA),
                                    (msafile_long_fas,  "ungzipped", "long", MIToS.MSA.FASTA),
                                    (msafile_wide_sth_gz, "gzipped", "wide", MIToS.MSA.Stockholm),
                                    (msafile_wide_sth,  "ungzipped", "wide", MIToS.MSA.Stockholm),
                                    (msafile_wide_fas_gz, "gzipped", "wide", MIToS.MSA.FASTA),
                                    (msafile_wide_fas,  "ungzipped", "wide", MIToS.MSA.FASTA)]

    # Default parser
    msa["input"][string(shape,"_",gzipped,"_",format)] = @benchmarkable read($file, $format)::MIToS.MSA.AnnotatedMultipleSequenceAlignment
    # With mapping
    msa["input"][string(shape,"_",gzipped,"_",format,"_mapping")] = @benchmarkable read($file, $format, generatemapping=true, useidcoordinates=true)::MIToS.MSA.AnnotatedMultipleSequenceAlignment
end

##### output
msa["output"] = BenchmarkGroup()

for (file,gzipped,shape,format) in [(msafile_long_sth_gz, "gzipped", "long", MIToS.MSA.Stockholm),
                                    (msafile_long_sth,  "ungzipped", "long", MIToS.MSA.Stockholm),
                                    (msafile_long_fas_gz, "gzipped", "long", MIToS.MSA.FASTA),
                                    (msafile_long_fas,  "ungzipped", "long", MIToS.MSA.FASTA),
                                    (msafile_wide_sth_gz, "gzipped", "wide", MIToS.MSA.Stockholm),
                                    (msafile_wide_sth,  "ungzipped", "wide", MIToS.MSA.Stockholm),
                                    (msafile_wide_fas_gz, "gzipped", "wide", MIToS.MSA.FASTA),
                                    (msafile_wide_fas,  "ungzipped", "wide", MIToS.MSA.FASTA)]
    outfile = string("./tmp/",split(file,"/")[end])
    aln = read(file, format)
    msa["output"][string(shape,"_",format,"_",gzipped)] = @benchmarkable write($outfile, $aln, $format)
end

##### Identity
msa["identity"] = BenchmarkGroup()

for (aln,label) in ((msa_long,"long"), (msa_wide,"wide"))
    for t in (Float16,Float32)
        msa["identity"][string("matrix_",t,"_",label)] = @benchmarkable percentidentity($aln, $t)
    end
    if label == "wide"
        msa["identity"][string("matrix_Float64_wide")] = @benchmarkable percentidentity($aln, Float64)
    end
    msa["identity"][string("mean_",label)] = @benchmarkable meanpercentidentity($aln)
end

##### Clustering
msa["hobohmI"] = BenchmarkGroup()

for (aln,label) in ((msa_long,"long"), (msa_wide,"wide"))
    for pid in 10:10:90
        msa["hobohmI"][string(pid,"_",label)] = @benchmarkable hobohmI($aln, $pid)
    end
end

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

# --------------------------------------------------------------------------- #
# Information module benchmarks
const information = BenchmarkGroup()

##### estimateincolumns
information["estimateincolumns"] = BenchmarkGroup()

for (aln,label) in ((msa_long,"long"), (msa_wide,"wide"))
    residues = getresidues(aln)
    information["estimateincolumns"][string("Entropy_Probability_",label)] = @benchmarkable estimateincolumns($residues, Int, ResidueProbability{Float64,1,false}, Entropy{Float64}())
    information["estimateincolumns"][string("Entropy_Count_",label)] = @benchmarkable estimateincolumns($residues, ResidueCount{Int,2,false}, Entropy{Float64}())
    information["estimateincolumns"][string("MI_Probability_",label)] = @benchmarkable estimateincolumns($residues, Int, ResidueProbability{Float64,2,false}, MutualInformation{Float64}())
    information["estimateincolumns"][string("MI_Count_",label)] = @benchmarkable estimateincolumns($residues, ResidueCount{Int,2,false}, MutualInformation{Float64}())
    information["estimateincolumns"][string("Entropy_Probability_Gapped_",label)] = @benchmarkable estimateincolumns($residues, Int, ResidueProbability{Float64,1,true}, Entropy{Float64}())
    information["estimateincolumns"][string("Entropy_Count_Gapped_",label)] = @benchmarkable estimateincolumns($residues, ResidueCount{Int,2,true}, Entropy{Float64}())
    information["estimateincolumns"][string("MI_Probability_Gapped_",label)] = @benchmarkable estimateincolumns($residues, Int, ResidueProbability{Float64,2,true}, MutualInformation{Float64}())
    information["estimateincolumns"][string("MI_Count_Gapped_",label)] = @benchmarkable estimateincolumns($residues, ResidueCount{Int,2,true}, MutualInformation{Float64}())
end

##### high level
information["highlevel"] = BenchmarkGroup()

for (aln,label) in ((msa_long,"long"), (msa_wide,"wide"))
    # default
    information["highlevel"][string("Buslje09_",label)] = @benchmarkable buslje09($aln)
    information["highlevel"][string("ZBLMIp_",label)] = @benchmarkable BLMI($aln)
end

##### low level
information["lowlevel"] = BenchmarkGroup()

const residues = getresidues(msa_long)
const column_i = residues[:,10]
const column_j = residues[:,20]
const column_k = residues[:,30]
const Pij = probabilities(column_i, column_j)
const Gij = ContingencyTable{Float64,  2, UngappedAlphabet}(UngappedAlphabet())
const nseq_msa_long = nsequences(msa_long)
const clusters_long = hobohmI(msa_long, 62)

information["lowlevel"]["count_col"] = @benchmarkable count($column_i)
information["lowlevel"]["count_col_col"] = @benchmarkable count($column_i, $column_j)
information["lowlevel"]["count_col_col_col"] = @benchmarkable count($column_i, $column_j, $column_k)
information["lowlevel"]["count_col_clusters"] = @benchmarkable count($column_i, weight=$clusters_long)
information["lowlevel"]["count_col_col_clusters"] = @benchmarkable count($column_i, $column_j, weight=$clusters_long)
information["lowlevel"]["count_col_col_col_clusters"] = @benchmarkable count($column_i, $column_j, $column_k, weight=$clusters_long)
information["lowlevel"]["probabilities_col"] = @benchmarkable probabilities($column_i)
information["lowlevel"]["probabilities_col_col"] = @benchmarkable probabilities($column_i, $column_j)
information["lowlevel"]["probabilities_col_col_col"] = @benchmarkable probabilities($column_i, $column_j, $column_k)

information["lowlevel"]["blosum_pseudofrequencies"] = @benchmarkable blosum_pseudofrequencies!($Gij, $Pij)
information["lowlevel"]["probabilities_blosum"] = @benchmarkable probabilities($column_i, $column_j, pseudofrequencies=BLOSUM_Pseudofrequencies(nsequences(msa), 8.512))

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

# --------------------------------------------------------------------------- #
# PDB module benchmarks
const mitos_pdb = BenchmarkGroup()

## Set up
### Files
const pdb_pdb_gz = "../../data/2XWB.pdb.gz"
const pdb_xml_gz = "../../data/2XWB.xml.gz"
const pdb_pdb    = "../../data/2XWB.pdb"
const pdb_xml    = "../../data/2XWB.xml"
### Residues
const pdb_residues = read(pdb_xml_gz, PDBML)

#### Parse benchmarks

##### input
mitos_pdb["input"] = BenchmarkGroup()

for (file,gzipped,label,format) in [(pdb_pdb_gz, "gzipped", "pdb", PDBFile),
                                    (pdb_xml_gz, "gzipped", "xml", PDBML),
                                    (pdb_pdb,  "ungzipped", "pdb", PDBFile),
                                    (pdb_xml,  "ungzipped", "xml", PDBML)]
    # Default parser
    mitos_pdb["input"][string(label,"_",format,"_",gzipped)] = @benchmarkable read($file, $format)
end

##### output
mitos_pdb["output"] = BenchmarkGroup()

for (file,gzipped,label,format) in [(pdb_pdb_gz, "gzipped", "pdb", PDBFile),
                                    (pdb_xml_gz, "gzipped", "xml", PDBML),
                                    (pdb_pdb,  "ungzipped", "pdb", PDBFile),
                                    (pdb_xml,  "ungzipped", "xml", PDBML)]
    outfile = string("./tmp/",split(file,"/")[end])
    mitos_pdb["output"][string(label,"_",format,"_",gzipped)] = @benchmarkable write($outfile, $pdb_residues, $format)
end

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

# --------------------------------------------------------------------------- #
# Pfam module (pipeline) benchmarks
const pfam = BenchmarkGroup()

# Set up
const aln = read(msafile_long_sth_gz, Stockholm, generatemapping=true, useidcoordinates=true)
const col2res = msacolumn2pdbresidue(aln, "CFAB_HUMAN/481-752", "2XWB", "F", "PF00089","../../data/2xwb.xml")
const pdb = read("../../data/2XWB.xml.gz", PDBML)
const resdict = @residuesdict pdb model "1" chain "F" group "ATOM" residue All
const cmap = msacontacts(aln, resdict, col2res)
const ZMIp, MIp = buslje09(aln)

pfam["read_pfam_gzipped"] = @benchmarkable read($msafile_long_sth_gz, Stockholm, generatemapping=true, useidcoordinates=true)
pfam["getseq2pdb"] = @benchmarkable getseq2pdb($aln)
pfam["msacolumn2pdbresidue_sifts"] = @benchmarkable msacolumn2pdbresidue($aln, "CFAB_HUMAN/481-752", "2XWB", "F", "PF00089","../../data/2xwb.xml")
pfam["msacolumn2pdbresidue_sifts_gzipped"] = @benchmarkable msacolumn2pdbresidue($aln, "CFAB_HUMAN/481-752", "2XWB", "F", "PF00089","../../data/2xwb.xml.gz")
pfam["read_PDBML_gzipped"] = @benchmarkable read("../../data/2XWB.xml.gz", PDBML)
pfam["residue_list_to_dict"] = @benchmarkable residuesdict($pdb,"1","F","ATOM",All)
pfam["msaresidues"] = @benchmarkable msaresidues($aln, $resdict, $col2res)
pfam["hasresidues"] = @benchmarkable hasresidues($aln, $col2res)
pfam["contact_map"] = @benchmarkable msacontacts($aln, $resdict, $col2res)
pfam["buslje09"] = @benchmarkable buslje09($aln)
pfam["AUC"] = @benchmarkable AUC($ZMIp, $cmap)

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

function SetUp!(;module_msa::Bool=true,module_information::Bool=true,module_pdb::Bool=true,module_pfam::Bool=true)
    if module_msa
        tune!(msa)
        JLD.save("msa.jld", "msa", params(msa))
    end
    if module_information
        tune!(information)
        JLD.save("information.jld", "information", params(information))
    end
    if module_pdb
        tune!(mitos_pdb)
        JLD.save("mitos_pdb.jld", "mitos_pdb", params(mitos_pdb))
    end
    if module_pfam
        tune!(pfam)
        JLD.save("pfam.jld", "pfam", params(pfam))
    end
end

function Run!(;module_msa::Bool=true,module_information::Bool=true,module_pdb::Bool=true,module_pfam::Bool=true)
    loadparams!(msa, JLD.load("msa.jld", "msa"), :evals, :samples);
    loadparams!(information, JLD.load("information.jld", "information"), :evals, :samples);
    loadparams!(mitos_pdb, JLD.load("mitos_pdb.jld", "mitos_pdb"), :evals, :samples);
    loadparams!(pfam, JLD.load("pfam.jld", "pfam"), :evals, :samples);
    bench = Dict{ASCIIString,BenchmarkGroup}()
    if module_msa
        bench["msa"] = run(msa)
    end
    if module_information
        bench["information"] = run(information)
    end
    if module_pdb
        bench["mitos_pdb"] = run(mitos_pdb)
    end
    if module_pfam
        bench["pfam"] = run(pfam)
    end
    bench
end

# @elapsed include("Benchmark.jl") # 486 s
# @elapsed SetUp!() # 5907 s
# @elapsed result = Run!() # 5642 s
# JLD.save("result_Benchmark.jld", "result", result);

max_stars_repo_name: diegozea/mitos-benchmarks

max_issues_repo_name: diegozea/mitos-benchmarks

max_forks_repo_name: diegozea/mitos-benchmarks


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Path

A path is a sequence of indices that can be
used to traverse a given domain/data object.
"""
abstract type Path end

"""
    traverse(object, path)

Traverse `object` with `path`.
"""
function traverse end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("paths/linear.jl")
include("paths/random.jl")
include("paths/source.jl")
include("paths/shifted.jl")


text:
# Utility operators for computing prox iteratively, e.g. using CG

import Base: *, size, eltype
import LinearAlgebra: mul!

abstract type LinOp end

infer_shape_of_y(Op, ::AbstractVector) = (size(Op, 1), )
infer_shape_of_y(Op, x::AbstractMatrix) = (size(Op, 1), size(x, 2))

function (*)(Op::LinOp, x)
    y = zeros(promote_type(eltype(Op), eltype(x)), infer_shape_of_y(Op, x))
    mul!(y, Op, x)
end

size(Op::LinOp, i::Integer) = i <= 2 ? size(Op)[i] : 1

# AAc (Gram matrix)

mutable struct AAc{M, T} <: LinOp
    A::M
    buf::Maybe{AbstractArray{T}}
    function AAc{M, T}(A::M) where {M, T}
        new(A, nothing)
    end
end

AAc(A::M) where {M} = AAc{M, eltype(A)}(A)

function mul!(y, Op::AAc, x)
    if Op.buf === nothing
        Op.buf = adjoint(Op.A) * x
    else
        mul!(Op.buf, adjoint(Op.A), x)
    end
    mul!(y, Op.A, Op.buf)
end

mul!(y, Op::Adjoint{AAc}, x) = mul!(y, adjoint(Op), x)

size(Op::AAc) = size(Op.A, 1), size(Op.A, 1)
eltype(Op::AAc) = eltype(Op.A)

# AcA (Covariance matrix)

mutable struct AcA{M, T} <: LinOp
    A::M
    buf::Maybe{AbstractArray{T}}
    function AcA{M, T}(A::M) where {M, T}
        new(A, nothing)
    end
end

AcA(A::M) where {M} = AcA{M, eltype(A)}(A)

function mul!(y, Op::AcA, x)
    if Op.buf === nothing
        Op.buf = Op.A * x
    else
        mul!(Op.buf, Op.A, x)
    end
    mul!(y, adjoint(Op.A), Op.buf)
end

mul!(y, Op::Adjoint{AcA}, x) = mul!(y, adjoint(Op), x)

size(Op::AcA) = size(Op.A, 2), size(Op.A, 2)
eltype(Op::AcA) = eltype(Op.A)

# Shifted symmetric linear operator

struct ScaleShift{M, T} <: LinOp
    alpha::T
    A::M
    rho::T
    function ScaleShift{M, T}(alpha::T, A::M, rho::T) where {M, T}
        if eltype(A) != T
            error("type of alpha, rho ($T) is different from that of A ($(eltype(A)))")
        end
        new(alpha, A, rho)
    end
end

ScaleShift(alpha::T, A::M, rho::T) where {M, T} = ScaleShift{M, T}(alpha, A, rho)

function mul!(y, Op::ScaleShift, x)
    mul!(y, Op.A, x)
    y .*= Op.alpha
    y .+= Op.rho .* x
end

function mul!(y, Op::Adjoint{ScaleShift}, x)
    mul!(y, adjoint(Op.A), x)
    y .*= Op.alpha
    y .+= Op.rho .* x
end

size(Op::ScaleShift) = size(Op.A, 2), size(Op.A, 2)
eltype(Op::ScaleShift) = eltype(Op.A)


text:
module TestGrad
    using Calculus2
    using Base.Test

    # multivariate test function
    function f(x::Vector)
        return sin(x[1]) + cos(x[2])
    end
    x = [0.0, 0.0]
    gr = similar(x)

    # Test in-place multivariate grad!
    grad!(gr, f, x)
    @test norm(gr - [cos(0.0), -sin(0.0)]) < 10e-8
    fill!(gr, 0.0)

    grad!(gr, f, x, method = :ad)
    @test norm(gr - [cos(0.0), -sin(0.0)]) < 10e-8
    fill!(gr, 0.0)

    ∇!(gr, f, x)
    @test norm(gr - [cos(0.0), -sin(0.0)]) < 10e-8
    fill!(gr, 0.0)

    # Test pure multivariate grad
    @test norm(grad(f, x) - [cos(0.0), sin(0.0)]) < 10e-8
    @test norm(∇(f, x) - [cos(0.0), sin(0.0)]) < 10e-8

    # Test in-place multivariate gradfunc!
    ∇f! = gradfunc!(f, x)
    ∇f!(gr, x)
    @test norm(gr - [cos(0.0), sin(0.0)]) < 10e-8

    # Test in-place multivariate gradfunc!
    ∇f! = gradfunc!(f, x, method = :ad)
    ∇f!(gr, x)
    @test norm(gr - [cos(0.0), sin(0.0)]) < 10e-8

    # Test pure multivariate gradfunc
    ∇f = gradfunc(f, x)
    @test norm(∇f(x) - [cos(0.0), sin(0.0)]) < 10e-8

    # Test pure univariate grad
    @test isapprox(grad(sin, 0.0), cos(0.0))
    @test isapprox(grad(sin, 0.0, method = :ad), cos(0.0))

    # Test pure univariate gradfunc
    @test isapprox(gradfunc(sin, 0.0)(0.0), cos(0.0))
    @test isapprox(gradfunc(sin, 0.0, method = :ad)(0.0), cos(0.0))

    # TODO: Test automatic differentation
    # TODO: Test symbolic differentiation
    # TODO: Test direction options
end


text:
# Advent of Code 2017 - day 2

f = open("input2_r")
input = readlines(f)
close(f)

function parse_nums(s)
  [parse(Int,n) for n in split(s)]
end

spreadsheet = [parse_nums(s) for s in input]

function diff(l)
  maximum(l)-minimum(l)
end

checksum = sum([diff(l) for l in spreadsheet])

println("part 1: ", checksum)

function divides(x,y)
  y !=x && rem(y,x)==0
end

function even_div(l)
  for x in l
    for y in l
      if divides(x,y) return div(y,x) end
    end
  end
end

div_sum = sum([even_div(l) for l in spreadsheet])

println("part 2: ", div_sum)

max_stars_repo_name: venanzio/AoC

max_issues_repo_name: venanzio/AoC

max_forks_repo_name: venanzio/AoC


text:
using ProgressLogging
using Setfield

"""
    GreeDS(alg=PCA(); threshold=0.0, progress=true)
    GreeDS(ncomps; threshold=0.0, progress=true)

Performs the greedy disk subtraction (GreeDS) algorithm.

This method is an iterative approach to standard ADI reduction which seeks to minimize over-subtraction by constraining the low-rank matrix approximation from `alg`.

For large data cubes the iteration can cause slowdowns, so a progress bar is provided using the [`ProgressLogging`](https://github.com/JunoLab/ProgressLogging.jl) API along with the `progress` keyword. It won't appear without a logging backend, such as [`TerminalLoggers`](https://github.com/c42f/TerminalLoggers.jl).

# Algorithms
Although the original paper explicitly uses PCA, we allow use of any linear ADI algorithm that is characterized by `ncomps`. By default, uses [`PCA`](@ref). If `ncomps` are provided, will default to [`PCA`](@ref) as well.

# References
1. [Pairet et al. 2018](https://ui.adsabs.harvard.edu/abs/2018arXiv181201333P) "Reference-less algorithm for circumstellar disks imaging"
2. [Pairet et al. 2020](https://ui.adsabs.harvard.edu/abs/2020arXiv200805170P) "MAYONNAISE: a morphological components analysis pipeline for circumstellar disks and exoplanets imaging in the near infrared"
"""
struct GreeDS{ALG<:LinearAlgorithm} <: LinearAlgorithm
    alg::ALG
    threshold::Float64
    progress::Bool
end

GreeDS(alg=PCA(); threshold=0.0, progress=true) = GreeDS(alg, threshold, progress)
GreeDS(ncomps::Int; kwargs...) = GreeDS(PCA(ncomps); kwargs...)

function decompose(alg::GreeDS, cube, angles, cube_ref=cube; kwargs...)
    target = expand(cube_ref)
    # get the number of components as a range from the underlying alg
    max_ncomps = isnothing(alg.alg.ncomps) ? size(target, 1) : alg.alg.ncomps
    # use the underlyhing algorithm with a lens for further processing
    _alg = alg.alg
    _alg = @set _alg.ncomps = 1
    reduced = _alg(target, angles)
    local A
    @progress "GreeDS" for n in 1:max_ncomps
        resid = target .- expand_rotate(reduced, angles, alg.threshold; kwargs...)
        # use lens to update number of components
        _alg = @set _alg.ncomps = n
        # use the `resid` cube as the reference frames for the next reduction
        A, w = decompose(_alg, target, angles, resid; kwargs...)
        reduced = collapse!(target .- expand(w * A), angles; kwargs...)
    end
    # project target cube onto reference basis
    w = flatten(cube) * A'
    return A, w
end

"""
    expand_rotate(frame, angles, threshold; kwargs...)

Takes a frame, expands it into a cube, rotates it clockwise by `angles`, and min-clips at `threshold`. Keyword arguments will be passed to `HCIToolbox.derotate!`.
"""
function expand_rotate(frame, angles, threshold; kwargs...)
    N = length(angles)
    _frame = @. ifelse(frame > threshold, frame, threshold)
    cube = similar(frame, N, size(frame)...)
    Threads.@threads for idx in axes(cube, 1)
        cube[idx, :, :] .= derotate(_frame, -angles[idx]; kwargs...)
    end
    return cube
end


text:
# using Base.Test
#include("../src/RunTests.jl")

using PhysicalConstant.CODATA2019, Measurements, Unitful
if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

import PhysicalConstant.CODATA2019: α, atm, c_0, e, ε_0, h, ħ, µ_0, m_e, m_H, C_0, ƛe, CmB, r_0

@testset "Base" begin
    @test ustrip(big(h)) == big"6.626070040e-34"
    @test setprecision(BigFloat, 768) do; precision(ustrip(big(c_0))) end == 768
    @test measurement(h) === measurement(h)
    @test iszero(measurement(α) - measurement(α))
    @test isone(measurement(BigFloat, atm) / measurement(BigFloat, atm))
    @test iszero(measurement(BigFloat, ħ) - (measurement(BigFloat, h) / 2big(pi)))
    @test isone(measurement(BigFloat, ħ) / (measurement(BigFloat, h) / 2big(pi)))
end

@testset "Utils" begin
    @test c_0.val === ustrip(float(c_0))
    @test C_0.val === ustrip(float(C_0))
    @test_throws ErrorException h.foo
end

@testset "Promotion" begin
    x = @inferred(inv(μ_0 * c_0 ^ 2))
    T = promote_type(typeof(ε_0), typeof(x))
    u = u"A^2 * kg^-1 * m^-3 * s^4"
    @test promote_type(typeof(α), BigInt) === BigFloat
    @test promote_type(typeof(α), typeof(1u"m/cm")) === Float64
    @test T === Quantity{Float64, dimension(x), typeof(u)}
    @test convert(T, ε_0) === (ε_0 / unit(ε_0)) * u
    @test convert(Float32, α) === float(Float32, α)
    @test uconvert(u"cm/s", c_0) === uconvert(u"cm/s", float(c_0))
end

@testset "Maths" begin
    @test α ≈ @inferred(e^2/(4 * pi * ε_0 * ħ * c_0))
    @test @inferred(ħ / (c_0 * m_e)) ≈ (ħ / float(c_0 * m_e))
#    @test ƛe ≈ @inferred(ħ / (c_0 * m_e))
    @test @inferred(α + 2) ≈ 2 + float(α)
    @test @inferred(5 + α) ≈ float(α) + 5
    @test @inferred(α + 2.718) ≈ 2.718 + float(α)
    @test @inferred(-3.1415 + α) ≈ float(α) - 3.1415
    @test ε_0 ≈ @inferred(1 / (μ_0 * c_0 ^ 2))
    @test @inferred(big(0) + α) == big(α)
    @test @inferred(α * 1.0) == float(α)
end

@testset "Show" begin
    @test repr(c_0) ==
        "Speed of light in vacuum (c_0)
Value                         = 2.99792458e8 m s^-1
Standard uncertainty          = (exact)
Relative standard uncertainty = (exact)
Reference                     = CODATA 2019"
    @test repr(α) ==
        "Fine-structure constant (α)
Value                         = 0.0072973525664
Standard uncertainty          = 1.7e-12
Relative standard uncertainty = 2.3e-10
Reference                     = CODATA 2019"
    @test repr(e) ==
        "Elementary charge (e)
Value                         = 1.6021766208e-19 C
Standard uncertainty          = 9.8e-28 C
Relative standard uncertainty = 6.1e-9
Reference                     = CODATA 2019"
    @test repr(ħ) ==
        "Planck constant over 2pi (ħ)
Value                         = 1.0545718001391127e-34 J s
Standard uncertainty          = 1.2891550390443523e-42 J s
Relative standard uncertainty = 1.2e-8
Reference                     = CODATA 2019"
    @test repr(r_0) ==
        "Bare Hydrogen radius (r_0)
Value                         = 5.291772103e-11 m
Standard uncertainty          = 1.2e-20 m
Relative standard uncertainty = 2.3e-10
Reference                     = CODATA 2019"
end


text:
function hofstQseq(n, typerst::Type=Int)
    nmax = maximum(n)
    r = Vector{typerst}(nmax)
    r[1] = 1
    if nmax ≥ 2 r[2] = 1 end
    for i in 3:nmax
        r[i] = r[i - r[i - 1]] + r[i - r[i - 2]]
    end
    return r[n]
end

println("First ten elements of sequence: ", join(hofstQseq(1:10), ", "))
println("1000-th element: ", hofstQseq(1000))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using ProjAssistant
@quickactivate "Chemostat_Rath2017"

@time begin

    using Plots
    using Statistics
    using ProgressMeter
    using Base.Threads

    import Chemostat
    import Chemostat.MetNets
    import Chemostat.MetLP
    import Chemostat.MetEP
    const Ch = Chemostat

    import Chemostat_Rath2017
    const ChR = Chemostat_Rath2017
    const Rd = ChR.RathData
    const H1 = ChR.Human1
    const HG = H1.HumanGEM
    const AG = H1.HartGEMs

end

## ---------------------------------------------------------------------
# # Load model
# let
#     fn = "/Users/Pereiro/University/Research/Metabolism/MaxEntEP2020/WIP/Chemostat_Rath2017/data/raw/Human1/Human1_Publication_Data_Scripts/ec_GEMs/models/HOP62/ecModel_batch.mat"
#     ec_model = MetNets.read_mat(fn)
#     MetNets.clampfields!(ec_model, [:lb, :ub]; abs_max = 999999, zeroth = 1e-8)
#     objidx = MetNets.rxnindex(ec_model, HG.HUMAN_BIOMASS_IDER)
#     ec_model = ChR.prepare_metnet(AG, ec_model; inf_medium = true);
#     ec_model = MetNets.force_dims(ec_model)
#     # fbaout = MetLP.fba!(ec_model, objidx)
#     sdat(AG, MetNets.compressed_model(ec_model), "ec_model", ".jls")
# end

## ---------------------------------------------------------------------
let

    met_readable = HG.load_met_readable_ids()
    met_map = HG.load_exch_met_map()

    @info("Loading model")
    # model = MetNets.uncompressed_model(ldat(AG, "ec_model", ".jls"))
    tissue = "GBM"
    modelid = "base"
    model = AG.load_model(;modelid, tissue, uncompress = true)
    objidx = MetNets.rxnindex(model, HG.HUMAN_BIOMASS_IDER)
    exglcidx = MetNets.rxnindex(model, HG.HUMAN_GLC_EX_IDER)

    @info("Finding open exchanges")
    exchis = filter(MetNets.exchanges(model)) do exchi
        MetNets.lb(model, exchi) < 0
    end
    @show length(exchis)

    # Open all
    @info("Fully opening all")
    for exchi in exchis
        MetNets.lb!(model, exchi, -1000)
    end
    fbaout = MetLP.fba!(model, objidx)
    biom0 = MetLP.objval(fbaout)
    @show biom0

    # find assentials
    @info("Finding assentials")
    cid = (:ESSENTIALS, exchis)
    # delcache(cid) # Reset
    essentialis = lcache(AG, cid) do
        _essentialis = Int[]
        for exchi in exchis
            lb_ = MetNets.lb(model, exchi)
            MetNets.lb!(model, exchi, 0.0)

            fbaout = MetLP.fba!(model, objidx)
            biom = MetLP.objval(fbaout)
            (biom == 0.0) && push!(_essentialis, exchi)

            exch = model.rxns[exchi]
            met = get(met_map, exch, exch)
            met = get(met_readable, met, met)

            MetNets.lb!(model, exchi, lb_)
        end
        return _essentialis
    end
    @show length(essentialis)

    @info("non_essentialis")
    non_essentialis = setdiff(exchis, essentialis)
    @show length(non_essentialis)

    @info("close non_essentials")
    for exchi in exchis
        MetNets.lb!(model, exchi, -1000.0)
    end
    for exchi in non_essentialis
        MetNets.lb!(model, exchi, 0.0)
    end

    @info("Test only essentials")
    fbaout = MetLP.fba!(model, objidx)
    biom = MetLP.objval(fbaout)
    @show biom

    @info("Test adding glucose")
    MetNets.lb!(model, exglcidx, -1000.0)
    fbaout = MetLP.fba!(model, objidx)
    biom = MetLP.objval(fbaout)
    @show biom

    # I will find the minimum medium with make glucose essential
    @threads for _ in collect(1:10000)
        thid = threadid()

        # local copy
        model_th = MetNets.MetNet(model;
            lb = copy(model.lb),
            ub = copy(model.ub),
            c = copy(model.c)
        )

        @info("Finding a greedy minimum medium", thid)
        minimum_mediumis = Set(exchis)
        for it in 1:(3*length(exchis))

            # Open minimum_mediumis
            for exchi in exchis
                MetNets.lb!(model_th, exchi, -0.0)
            end
            for exchi in minimum_mediumis
                MetNets.lb!(model_th, exchi, -1000.0)
            end
            MetNets.lb!(model_th, exglcidx, -1000.0)

            # close one
            rexchi = rand(setdiff(minimum_mediumis, essentialis))
            exch = model_th.rxns[rexchi]
            rmet = get(met_map, exch, exch)
            rmet = get(met_readable, rmet, rmet)
            rmet == "glucose[s]" && continue # protect glucose
            MetNets.lb!(model_th, rexchi, 0.0)

            fbaout = MetLP.fba!(model_th, objidx)
            biom = MetLP.objval(fbaout)

            # If rexch is not essential
            if abs(biom) > 0.9 * abs(biom0)

                @info("Met NOT essential", it, length(minimum_mediumis), rmet, biom0, biom, thid)

                # remove
                delete!(minimum_mediumis, rexchi)

                # Test glc essensiality
                MetNets.lb!(model_th, exglcidx, 0.0)
                fbaout = MetLP.fba!(model_th, objidx)
                biom = MetLP.objval(fbaout)
                # If glc is essential
                if abs(biom) < 0.1 * abs(biom0)

                    @info("GLC essential", it, length(minimum_mediumis), biom0, biom, thid)
                    @show minimum_mediumis

                    # save
                    medium_hash = hash((:MIN_MEDIUM, minimum_mediumis))
                    sdat(AG, minimum_mediumis,
                        "minimum_medium", (;medium_hash), ".jls";
                        verbose = true
                    )
                else
                    # @info("GLC NOT essential", it, length(minimum_mediumis), biom0, biom, thid)
                end
            else
                # If rexch is essential (open)
                # @info("Met essential", it, length(minimum_mediumis), rmet, biom0, biom, thid)
                MetNets.lb!(model_th, rexchi, -1000.0)
                continue
            end
            println()
        end

    end # threads

end

## ---------------------------------------------------------------------
## ---------------------------------------------------------------------
# ## ---------------------------------------------------------------------
# # histogram
# let
#     dir = datdir(AG)
#     lens = Int[]
#     medium_files = filter(readdir(dir; join = true)) do fn
#         contains(fn, "minimum_medium")
#     end
#     @show length(medium_files)
#     for fn in medium_files
#         push!(lens, length(ldat(fn)))
#     end
#     histogram(lens; bins = 20)
# end

# ## ---------------------------------------------------------------------
# let
#     dir = datdir(AG)
#     lens = Int[]
#     medium_files = filter(readdir(dir; join = true)) do fn
#         contains(fn, "minimum_medium")
#     end

#     met_readable = HG.load_met_readable_ids()
#     met_map = HG.load_exch_met_map()

#     @info("Loading model")
#     # model = MetNets.uncompressed_model(ldat(AG, "ec_model", ".jls"))
#     tissue = "GBM"
#     modelid = "base"
#     model = AG.load_model(;modelid, tissue, uncompress = true)
#     objidx = MetNets.rxnindex(model, HG.HUMAN_BIOMASS_IDER)

#     @info("Finding open exchanges")
#     exchis = filter(MetNets.exchanges(model)) do exchi
#         MetNets.lb(model, exchi) < 0
#     end
#     @show length(exchis)

#     # test
#     fbaout = MetLP.fba!(model, objidx)
#     biom0 = MetLP.objval(fbaout)
#     @show biom0

#     for medium_file in medium_files

#         println()
#         minimum_mediumis = ldat(medium_file)
#         @show minimum_mediumis
#         @show length(minimum_mediumis)

#         println()
#         @show biom0

#         # close all
#         for exchi in exchis
#             MetNets.lb!(model, exchi, 0.0)
#         end

#         # open minimum
#         for exchi in minimum_mediumis
#             MetNets.lb!(model, exchi, -1000.0)
#         end

#         # test
#         fbaout = MetLP.fba!(model, objidx)
#         biom = MetLP.objval(fbaout)
#         @show biom

#         # close glcose
#         MetNets.lb!(model, HG.HUMAN_GLC_EX_IDER, 0.0)

#         # test
#         fbaout = MetLP.fba!(model, objidx)
#         biom = MetLP.objval(fbaout)
#         @show biom

#     end

# end

# ## ---------------------------------------------------------------------
# ## ---------------------------------------------------------------------
# # let
# #     met_readable = HG.load_met_readable_ids()
# #     met_map = HG.load_exch_met_map()
# #     results = ldat(AG, "glc_limited_study", ".jls")

# #     li = 5
# #     biomis = Float64[]
# #     mets = String[]
# #     p = plot()
# #     for (exch, dat) in results
# #         met = get(met_map, exch, exch)
# #         met = get(met_readable, met, met)
# #         println("met: ", met)
# #         println("bioms: ", join(dat, ", "))
# #         println()
# #         # idx = findfirst(reverse(dat)) do biom
# #         #     biom < maximum(dat) * 0.9
# #         # end
# #         # idx = isnothing(idx) ? 0 : idx
# #         if any(iszero.(dat))
# #             # push!(biomis, idx)
# #             # push!(mets, met)
# #             plot!(p, dat; label = met)
# #         end
# #     end
# #     p
# #     # rang_ = 1:10
# #     # sis = sortperm(biomis; rev = true)
# #     # bar(mets[sis], biomis[sis]; xrotation = 45)

# #     # sfig(AG, ps,
# #     #     @fileid, "glc_limited_study", ".png";
# #     #     layout = (3,3)
# #     # )
# # end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -2.9922365365158297 <= q <= 1.3959867501971397)

                     add_NL_constraint(m, :(gelu(-0.7172866311917714 + 0.8530094465195428*gelu(0.49975691401319544 + 0.817850285997455*gelu(-0.9044530867886094 + -0.3504940249833268*gelu(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*gelu(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + -0.15011340339301915*gelu(-0.49866226232151334 + 0.31623027699891626*gelu(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*gelu(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4])))) + -0.2302663292873568*gelu(0.6056562038214484 + 0.017530899400252764*gelu(-0.9044530867886094 + -0.3504940249833268*gelu(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*gelu(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + 0.9141090581590978*gelu(-0.49866226232151334 + 0.31623027699891626*gelu(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*gelu(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))) + gelu(0.11802092343800519 + -0.23198020992542823*gelu(0.49975691401319544 + 0.817850285997455*gelu(-0.9044530867886094 + -0.3504940249833268*gelu(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*gelu(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + -0.15011340339301915*gelu(-0.49866226232151334 + 0.31623027699891626*gelu(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*gelu(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4])))) + 0.3694431087223533*gelu(0.6056562038214484 + 0.017530899400252764*gelu(-0.9044530867886094 + -0.3504940249833268*gelu(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*gelu(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + 0.9141090581590978*gelu(-0.49866226232151334 + 0.31623027699891626*gelu(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*gelu(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Use
#    @warnpcfail precompile(args...)
# if you want to be warned when a precompile directive fails
macro warnpcfail(ex::Expr)
    modl = __module__
    file = __source__.file === nothing ? "?" : String(__source__.file)
    line = __source__.line
    quote
        $(esc(ex)) || @warn """precompile directive
     $($(Expr(:quote, ex)))
 failed. Please report an issue in $($modl) (after checking for duplicates) or remove this directive.""" _file=$file _line=$line
    end
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(fetch_resource)),NamedTuple{(:servers,), Tuple{Vector{String}}},typeof(fetch_resource),String,String})
    Base.precompile(Tuple{typeof(get_registries),String,String})
    Base.precompile(Tuple{typeof(serve_json),HTTP.Streams.Stream{HTTP.Messages.Request, HTTP.ConnectionPool.Transaction{TCPSocket}},Dict{String, Any}})
    Base.precompile(Tuple{typeof(verify_registry_hash),String,String})
    Base.precompile(Tuple{typeof(wait_first),Task,Vararg{Task}})
    isdefined(PkgServer, Symbol("#22#26")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#22#26"))})
    isdefined(PkgServer, Symbol("#32#35")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#32#35"))})
    isdefined(PkgServer, Symbol("#40#42")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#40#42"))})
    isdefined(PkgServer, Symbol("#60#65")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#60#65"))})
    isdefined(PkgServer, Symbol("#61#66")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#61#66"))})
    isdefined(PkgServer, Symbol("#ServerConfig#57#58")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#ServerConfig#57#58")),InetAddr{IPv4},String,Dict{String, RegistryMeta},Vector{String},Vector{String},Int64,Type{ServerConfig}})
end

max_stars_repo_name: StefanKarpinski/PkgServer

max_issues_repo_name: StefanKarpinski/PkgServer

max_forks_repo_name: StefanKarpinski/PkgServer


text:
const DEFAULT_TRAIT_SUPERTYPE = Any

export trait

abstract type BinaryTrait{T} end

# This sub-module is used to keep standard prefix types
module Prefix
    import ..BinaryTraits
    using ..BinaryTraits: BinaryTrait

    # default types for both sides
    struct Positive{T} <: BinaryTrait{T} end
    struct Negative{T} <: BinaryTrait{T} end

    # Positive/Negative trait types are traits
    BinaryTraits.is_trait(::Type{Positive{T}}) where T = true
    BinaryTraits.is_trait(::Type{Negative{T}}) where T = true

    # Optional positive types that may be brought into user module namespace
    const Can = Positive
    const Has = Positive
    const Is = Positive

    # Optional negative types that may be brought into user module namespace
    const Cannot = Negative
    const IsNot = Negative
    const No = Negative
    const Not = Negative
end

"""
    trait(::Type{T}, x)

Returns the singleton positive/negative trait type for object `x`.
"""
function trait end

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

"""
    @trait <name> [as <category>] [prefix <positive>,<negative>] [with <trait1,trait2,...>]

Create a new trait type for `name` called `\$(name)Trait`:
* If the `as` clause is provided, then `category` (an abstract type) will be used as the
  super type of the trait type.
* If the `prefix` clause is provided, then it allows the user to choose different prefixes
  than the default ones (`Can` and `Cannot`) e.g. `prefix Is,Not` or `prefix Has,Not`.
* If the `with` clause is provided, then it defines a composite trait from existing traits.
Note that you must specify at least 2 traits to make a composite trait.
"""
macro trait(name::Symbol, args...)
    category, prefixes, underlying_traits = parse_trait_args(args)
    pos, neg = prefixes.args
    mod = __module__

    # Try to import the predefined prefix types to user module's namespace
    import_prefix_type(mod, pos)
    import_prefix_type(mod, neg)

    this_can_type = Expr(:curly, pos, name)
    this_cannot_type = Expr(:curly, neg, name)

    # Single traits - the default is "cannot".
    # Composite traits - the default is the AND-expression of all underlyings.
    default_expr = if underlying_traits !== nothing
        # Composite trait here:
        # Construct something like: trait(Fly,x) === Can{Fly}() && trait(Swim,x) === Can{Swim}()
        condition =
            Expr(:(&&),
                [let cap = t.args[1], trait = t.args[2]
                    :( BinaryTraits.trait($trait, x) === $(cap){$trait}())
                 end for t in underlying_traits.args]...)
        # Construct expression like: [condition] ? Can{FlySwim}() : Cannot{FlySwim}()
        Expr(:if, condition, Expr(:call, this_can_type), Expr(:call, this_cannot_type))
    else
        # Single trait here: default to "cannot" for every type
        Expr(:call, this_cannot_type)
    end

    # If it's composite trait, then I want to maintain a mapping from the
    # positive trait type to the underlying's positive trait types.
    # It is needed for interface checks later.
    composite_expr = if underlying_traits !== nothing
        underlying_types = underlying_traits.args
        quote
            BinaryTraits.push_composite_map!($mod, $this_can_type, Set([$(underlying_types...)]))
        end
    else
        nothing
    end

    expr = quote
        abstract type $name <: $category end
        BinaryTraits.trait(::Type{$name}, x::Type) = $default_expr
        BinaryTraits.is_trait(::Type{$name}) = true
        $composite_expr
    end

    display_expanded_code(expr)
    return esc(expr)
end

"""
Parse arguments for the @trait macro.
"""
function parse_trait_args(args)

    category = Symbol(DEFAULT_TRAIT_SUPERTYPE)
    prefixes = Expr(:tuple, :Positive, :Negative)
    traits = nothing

    usage = "Invalid @trait usage. See doc string for details."
    length(args) % 2 === 0 || throw(SyntaxError(usage))

    if length(args) > 0 && args[1] == :as
        category = args[2]
        category isa Symbol || throw(SyntaxError("Not a symbol: $category. $usage"))
        args = args[3:end]
    end

    if length(args) > 0 && args[1] == :prefix
        prefixes =  args[2]
        is_tuple_of_symbols(prefixes; n = 2) || throw(SyntaxError(usage))
        args = args[3:end]
    end

    if length(args) > 0 && args[1] == :with
        traits =  args[2]
        (traits isa Symbol || is_tuple_of_curly_expressions(traits; n = 2, op = >=)) ||
            throw(SyntaxError(usage))
        args = args[3:end]
    end

    return (category, prefixes, traits)
end

"""
    import_prefix_type(m::Module, prefix::Symbol)

Import predefined prexies `prefix` from BinaryTraits.Prefix module into
the specified module `m`. If the name is already bounded with a different value
in the module, an error is raised.
"""
function import_prefix_type(m::Module, prefix::Symbol)
    existing_names = names(m, all = true, imported = true)
    try
        prefix_type = Base.eval(m, prefix)
        parentmodule(prefix_type) == BinaryTraits.Prefix ||
            error("Unable to import prefix `$prefix` as it already exists in module $m")
    catch
        Base.eval(m, :(import BinaryTraits.Prefix: $prefix))
    end
end


text:
module Jasmin

using LinearAlgebra, Combinatorics

import Base.println

export AbstractStandard, BaseEnumeration, StandardSimplexe, solve!,  Simplexe,
     AbstractStatus, Optimal, Infeasible, Unknown, Unbounded,
     xstar, vstar

include("utils.jl")
include("Status.jl")
include("Standard/main.jl")

end # module


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 6538daa0-243a-11eb-2726-c9068913ab03
begin
	using PlutoUI, NeXLMatrixCorrection, Gadfly
	md"""
	# Plot the ϕ(ρz) Curve

	Plot the ϕ(ρz) curve as calculated using various common matrix correction models.
	"""
end

# ╔═╡ a3e07a00-243b-11eb-01e9-27ea5be6f71e
md"""
## Measurement Parameters

Material: $(@bind matstr TextField(default="Al"))

Beam energy (keV): $(@bind e0keV NumberField(0.0:0.1:100.0,default=29.0))

Characteristic X-ray: $(@bind cxrstr TextField(default="Mg K-L3"))

Take-Off Angle (°): $(@bind toa NumberField(1.0:0.1:89.0, default=40.0))
"""

# ╔═╡ 8f049f40-243a-11eb-2424-ef70b99f4423
plot([ XPP, CitZAF, Riveros1993, XPhi ], parse(Material,matstr), parse(CharXRay, cxrstr), 1.0e3*e0keV, deg2rad(toa))

# ╔═╡ c0ad8e42-2826-11eb-1b81-9f9a11b6c51d
begin
	cxr = parse(CharXRay,cxrstr)
	sh = inner(cxr)
md"""
  * E($cxr) = $(0.001*energy(cxr)) keV
  * Ec($sh) = $(0.001*energy(sh)) keV
  * U($sh @ $e0keV keV) = $(1.0e3*e0keV/energy(sh))
"""
end

# ╔═╡ Cell order:
# ╟─6538daa0-243a-11eb-2726-c9068913ab03
# ╟─a3e07a00-243b-11eb-01e9-27ea5be6f71e
# ╟─8f049f40-243a-11eb-2424-ef70b99f4423
# ╟─c0ad8e42-2826-11eb-1b81-9f9a11b6c51d

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
"""
    FileStack{X,K}

    FileStack{X,K}(filename, types, sizes, eachchunk, haschunks, write)

A wrapper object that holds file pointer and size/chunking
metadata for a multi-layered stack stored in a single file,
typically netcdf or hdf5.

`X` is a backend singleton like `GDALfile`, and `K` is a tuple
of `Symbol` keys.
"""
struct FileStack{X,K,F<:AbstractString,T<:NamedTuple{K},S<:NamedTuple{K},EC,HC}
    filename::F
    types::T
    sizes::S
    eachchunk::EC
    haschunks::HC
    write::Bool
end
function FileStack{X}(
    filename::F, types::T, sizes::S, eachchunk::EC, haschunks::HC, write::Bool
   ) where {X,K,F,T<:NamedTuple{K},S<:NamedTuple{K},EC,HC}
    FileStack{X,K,F,T,S,EC,HC}(filename, types, sizes, eachchunk, haschunks, write)
end

# FileStack has `X` parameter that is not recoverable from fields.
ConstructionBase.constructorof(::Type{<:FileStack{X}}) where {X} = FileStack{X}

filename(fs::FileStack) = fs.filename
Base.keys(fs::FileStack{<:Any,K}) where K = K
Base.values(fs::FileStack{<:Any}) = (fs[k] for k in keys(fs))
# Indexing FileStack returns a FileArray,
# referencing a specific key in the same file.
function Base.getindex(fs::FileStack{X}, key::Symbol) where X
    size = fs.sizes[key]
    eachchunk = fs.eachchunk[key]
    haschunks = fs.haschunks[key]
    T = fs.types[key]
    N = length(size)
    A = FileArray{X,T,N}(filename(fs), size, key, eachchunk, haschunks, fs.write)
end


text:
using Test
using DataStructures: OrderedDict
using Printf: @sprintf
using QuantumLattices.Mathematics.QuantumNumbers
using QuantumLattices.Interfaces: ⊕, ⊗, dimension, expand, permute, decompose, regularize
import Base.Iterators: Iterators

@testset "SQN" begin
    @test SQN|>fieldnames == (:Sz,)
    @test SQN|>periods == (Inf,)
end

@testset "PQN" begin
    @test PQN|>fieldnames == (:N,)
    @test PQN|>periods == (Inf,)
end

@testset "SPQN" begin
    @test SPQN|>fieldnames == (:N, :Sz)
    @test SPQN|>periods == (Inf, Inf)
end

@testset "Z2QN" begin
    @test Z2QN|>fieldnames == (:N,)
    @test Z2QN|>periods == (2,)
end

@testset "SQNS" begin
    @test SQNS(0.5) == AbelianNumbers('C', [SQN(-0.5), SQN(0.5)], [0, 1, 2], qnindptr)
end

@testset "PQNS" begin
    @test PQNS(1.0) == AbelianNumbers('C', [PQN(0.0), PQN(1.0)], [0, 1, 2], qnindptr)
end

@testset "SzPQNS" begin
    @test SzPQNS(0.5) == AbelianNumbers('C', [SPQN(0.0, 0.0), SPQN(1.0, 0.5)], [0, 1, 2], qnindptr)
end

@testset "SPQNS" begin
    @test SPQNS(0.5) == AbelianNumbers('C', [SPQN(0.0, 0.0), SPQN(1.0, -0.5), SPQN(1.0, 0.5), SPQN(2.0, 0.0)], [0, 1, 2, 3, 4], qnindptr)
end

@testset "Z2QNS" begin
    @test Z2QNS() == AbelianNumbers('C', [Z2QN(0.0), Z2QN(1.0)], [0, 1, 2], qnindptr)
end

@abeliannumber "CN" (:N,) (Inf,)
@abeliannumber "Z4" (:Z,) (4,)
@abeliannumber "CNZ4" (:N, :Z) (Inf, 4)

@testset "regularize" begin
    @test regularize(CNZ4, [1.5, 5.0]) == [1.5, 1.0]
    @test regularize(CNZ4, [4.0, -1.0]) == [4.0, 3.0]
    @test regularize(CNZ4, [1.5 4.0; 5.0 -1.0]) == [1.5 4.0; 1.0 3.0]
end

@testset "AbelianNumbers" begin
    qn1, qn2 = CNZ4(+1.0, +3.0), CNZ4(-1.0, +1.0)
    qns1 = AbelianNumbers('U', [qn1, qn2], [0, 2, 5], qnindptr)
    qns2 = AbelianNumbers('U', [qn1, qn2], [2, 3], qncounts)
    @test isequal(qns1, qns2)
    qns3 = AbelianNumbers(OrderedDict(qn1=>2, qn2=>3))
    qns4 = AbelianNumbers(OrderedDict(qn1=>1:2, qn2=>3:5))
    @test qns1 == qns2 == qns3 == qns4

    qns = AbelianNumbers('U', [qn1, qn2], [0, 3, 5], qnindptr)
    @test qns|>dimension == 5
    @test qns|>string == "QNS(2, 5)"
    @test qns|>length == 2
    @test qns|>eltype == CNZ4
    @test qns|>typeof|>eltype == CNZ4
    @test @sprintf("%s", qns) == "QNS(CNZ4(1.0, 3.0)=>1:3, CNZ4(-1.0, 1.0)=>4:5)"
    @test qns[1] == qn1
    @test qns[2] == qn2
    @test qns[1:2] == AbelianNumbers('G', [qn1, qn2], [0, 3, 5], qnindptr)
    @test qns[[2, 1]] == AbelianNumbers('G', [qn2, qn1], [0, 2, 5], qnindptr)
    @test qns|>collect == [qn1, qn2]
    @test qns|>Iterators.reverse|>collect == [qn2, qn1]
    @test qns|>keys|>collect == [qn1, qn2]
    @test values(qns, qnindptr)|>collect == [1:3, 4:5]
    @test values(qns, qncounts)|>collect == [3, 2]
    @test pairs(qns, qnindptr)|>collect == [qn1=>1:3, qn2=>4:5]
    @test pairs(qns, qncounts)|>collect == [qn1=>3, qn2=>2]
end

@testset "toordereddict" begin
    qn1, qn2 = CNZ4(1.0, 2.0), CNZ4(2.0, 3.0)
    qns = AbelianNumbers('U', [qn1, qn2], [2, 3], qncounts)
    @test toordereddict(qns, qnindptr) == OrderedDict(qn1=>1:2, qn2=>3:5)
    @test toordereddict(qns, qncounts) == OrderedDict(qn1=>2, qn2=>3)
end

@testset "arithmetic" begin
    qn = CNZ4(1.0, 3.0)
    @test qn|>dimension == 1
    @test qn|>typeof|>dimension == 1
    @test +qn == CNZ4(+1.0, +3.0)
    @test -qn == CNZ4(-1.0, +1.0)
    @test qn*4 == 4*qn == CNZ4(4.0, 0.0)
    @test qn^3 == CNZ4(3.0, 1.0)
    @test kron(CNZ4, CN(1.0), Z4(3.0)) == qn

    qn1, qn2 = CNZ4(1.0, 2.0), CNZ4(2.0, 3.0)
    @test qn1+qn2 == CNZ4(3.0, 1.0)
    @test qn1-qn2 == CNZ4(-1.0, 3.0)
    @test kron(qn1, qn2, signs=(+1, +1)) == CNZ4(+3.0, 1.0)
    @test kron(qn1, qn2, signs=(+1, -1)) == CNZ4(-1.0, 3.0)
    @test kron(qn1, qn2, signs=(-1, +1)) == CNZ4(+1.0, 1.0)
    @test kron(qn1, qn2, signs=(-1, -1)) == CNZ4(-3.0, 3.0)
    @test union(qn1, qn2, signs=(+1, +1)) == AbelianNumbers('G', [+qn1, +qn2], [0, 1, 2], qnindptr)
    @test union(qn1, qn2, signs=(+1, -1)) == AbelianNumbers('G', [+qn1, -qn2], [0, 1, 2], qnindptr)
    @test union(qn1, qn2, signs=(-1, +1)) == AbelianNumbers('G', [-qn1, +qn2], [0, 1, 2], qnindptr)
    @test union(qn1, qn2, signs=(-1, -1)) == AbelianNumbers('G', [-qn1, -qn2], [0, 1, 2], qnindptr)
    @test ⊗(qn1, qn2) == kron(qn1, qn2)
    @test ⊕(qn1, qn2) == union(qn1, qn2)

    qns = AbelianNumbers('U', [qn1, qn2], [0, 2, 4], qnindptr)
    @test +qns == AbelianNumbers('U', [+qn1, +qn2], [0, 2, 4], qnindptr)
    @test -qns == AbelianNumbers('U', [-qn1, -qn2], [0, 2, 4], qnindptr)
    @test qns*3 == 3*qns == AbelianNumbers('U', [qn1*3, qn2*3], [0, 2, 4], qnindptr)
    @test qns^2 == AbelianNumbers('G', [qn1+qn1, qn1+qn2, qn1+qn1, qn1+qn2, qn2+qn1, qn2+qn2, qn2+qn1, qn2+qn2], [2, 2, 2, 2, 2, 2, 2, 2], qncounts)

    @test qns+qn == qn+qns == AbelianNumbers('U', [qn1+qn, qn2+qn], [0, 2, 4], qnindptr)
    @test qns-qn == AbelianNumbers('U', [qn1-qn, qn2-qn], [0, 2, 4], qnindptr)
    @test qn-qns == AbelianNumbers('U', [qn-qn1, qn-qn2], [0, 2, 4], qnindptr)

    qns1 = AbelianNumbers('U', [+qn1, -qn2], [0, 2, 4], qnindptr)
    qns2 = AbelianNumbers('U', [-qn1, +qn2], [0, 3, 4], qnindptr)
    @test union(qns1, qns2, signs=(+1, +1)) == AbelianNumbers('G', [+qn1, -qn2, -qn1, +qn2], [2, 2, 3, 1], qncounts)
    @test union(qns1, qns2, signs=(+1, -1)) == AbelianNumbers('G', [+qn1, -qn2, +qn1, -qn2], [2, 2, 3, 1], qncounts)
    @test union(qns1, qns2, signs=(-1, +1)) == AbelianNumbers('G', [-qn1, +qn2, -qn1, +qn2], [2, 2, 3, 1], qncounts)
    @test union(qns1, qns2, signs=(-1, -1)) == AbelianNumbers('G', [-qn1, +qn2, +qn1, -qn2], [2, 2, 3, 1], qncounts)
    @test kron(qns1, qns2, signs=(+1, +1)) == AbelianNumbers('G', [+qn1-qn1, +qn1+qn2, +qn1-qn1, +qn1+qn2, -qn2-qn1, -qn2+qn2, -qn2-qn1, -qn2+qn2], [3, 1, 3, 1, 3, 1, 3, 1], qncounts)
    @test kron(qns1, qns2, signs=(+1, -1)) == AbelianNumbers('G', [+qn1+qn1, +qn1-qn2, +qn1+qn1, +qn1-qn2, -qn2+qn1, -qn2-qn2, -qn2+qn1, -qn2-qn2], [3, 1, 3, 1, 3, 1, 3, 1], qncounts)
    @test kron(qns1, qns2, signs=(-1, +1)) == AbelianNumbers('G', [-qn1-qn1, -qn1+qn2, -qn1-qn1, -qn1+qn2, +qn2-qn1, +qn2+qn2, +qn2-qn1, +qn2+qn2], [3, 1, 3, 1, 3, 1, 3, 1], qncounts)
    @test kron(qns1, qns2, signs=(-1, -1)) == AbelianNumbers('G', [-qn1+qn1, -qn1-qn2, -qn1+qn1, -qn1-qn2, +qn2+qn1, +qn2-qn2, +qn2+qn1, +qn2-qn2], [3, 1, 3, 1, 3, 1, 3, 1], qncounts)
    @test ⊕(qns1, qns2) == union(qns1, qns2)
    @test ⊗(qns1, qns2) == kron(qns1, qns2)
end

@testset "sort" begin
    qn1, qn2, qn3, qn4 = CNZ4(3.0, 2.0), CNZ4(4.0, 1.0), CNZ4(4.0, 2.0), CNZ4(3.0, 3.0)
    oldqns = AbelianNumbers('G', [qn1, qn2, qn3, qn4, qn2, qn3, qn1, qn4], [2, 1, 3, 4, 1, 1, 2, 3], qncounts)
    newqns, permutation = sort(oldqns)
    @test newqns == AbelianNumbers('C', [qn1, qn4, qn2, qn3], [4, 7, 2, 4], qncounts)
    @test permutation == [1, 2, 13, 14, 7, 8, 9, 10, 15, 16, 17, 3, 11, 4, 5, 6, 12]
end

@testset "findall" begin
    qn1, qn2, qn3 = CNZ4(2.0, 3.0), CNZ4(1.0, 2.0), CNZ4(0.0, 0.0)
    qns = AbelianNumbers('C', [qn2, qn1], [2, 3], qncounts)
    @test findall(qn1, qns, qncompression) == [2]
    @test findall(qn2, qns, qncompression) == [1]
    @test findall(qn3, qns, qncompression) == []
    @test findall(qn1, qns, qnexpansion) == [3, 4, 5]
    @test findall(qn2, qns, qnexpansion) == [1, 2]
    @test findall(qn3, qns, qnexpansion) == []

    qns = AbelianNumbers('G', [qn1, qn2, qn1], [2, 2, 1], qncounts)
    @test findall(qn1, qns, qncompression) == [1, 3]
    @test findall(qn2, qns, qncompression) == [2]
    @test findall(qn3, qns, qncompression) == []
    @test findall(qn1, qns, qnexpansion) == [1, 2, 5]
    @test findall(qn2, qns, qnexpansion) == [3, 4]
    @test findall(qn3, qns, qnexpansion) == []
end

@testset "filter" begin
    qn1, qn2 = CNZ4(1.0, 2.0), CNZ4(3.0, 0.0)
    qns = AbelianNumbers('G', [qn1, qn2, qn1, qn2], [1, 2, 3, 4], qncounts)
    @test filter((qn1, qn2), qns) == AbelianNumbers('G', [qn1, qn2, qn1, qn2], [1, 2, 3, 4], qncounts)
    @test filter(qn2, qns) == AbelianNumbers('G', [qn2, qn2], [2, 4], qncounts)
end

@testset "ukron" begin
    qn1, qn2 = CNZ4(1.0, 1.0), CNZ4(2.0, 3.0)
    qns1 = AbelianNumbers('U', [qn2, qn1], [4, 5], qncounts)
    qns2 = AbelianNumbers('U', [qn1, qn2], [2, 3], qncounts)
    qns, records = ukron(qns1, qns2, signs=(+1, -1))
    @test qns == AbelianNumbers('C', [qn1-qn2, zero(CNZ4), qn2-qn1], [15, 22, 8], qncounts)
    @test records == Dict(  (qn1-qn2) =>  Dict((qn1, qn2)=>1:15),
                            zero(CNZ4) => Dict((qn2, qn2)=>1:12, (qn1, qn1)=>13:22),
                            (qn2-qn1) =>  Dict((qn2, qn1)=>1:8)
                            )
end

@testset "expand" begin
    qn = CNZ4(3.0, 2.0)
    qns = AbelianNumbers(qn, 4)
    @test expand(qns, qnindices) == [1, 1, 1, 1]
    @test expand(qns, qncontents) == [qn, qn, qn, qn]
end

@testset "decompose" begin
    qn1, qn2 = CNZ4(0.0, 0.0), CNZ4(1.0, 1.0)
    qns = AbelianNumbers('U', [qn1, qn2], [0, 1, 2], qnindptr)
    target = CNZ4(2.0, 2.0)
    result = Set(((1, 1, 2, 2), (1, 2, 1, 2), (1, 2, 2, 1), (2, 1, 1, 2), (2, 1, 2, 1), (2, 2, 1, 1)))
    @test ⊆(Set(decompose((qns, qns, qns, qns), target, (1, 1, 1, 1), qnbruteforce, nmax=10)), result)
    @test ⊆(Set(decompose((qns, qns, qns, qns), target, (1, 1, 1, 1), qnmontecarlo, nmax=10)), result)
end

@testset "permute" begin
    qn1, qn2, qn3 = CNZ4(1.0, 2.0), CNZ4(3.0, 0.0), CNZ4(4.0, 1.0)
    qns = AbelianNumbers('G', [qn1, qn2, qn3], [2, 3, 4], qncounts)
    @test permute(qns, [3, 2, 1], qncompression) == AbelianNumbers('G', [qn3, qn2, qn1], [4, 3, 2], qncounts)
    @test permute(qns, [4, 6, 9, 8], qnexpansion) == AbelianNumbers('G', [qn2, qn3, qn3, qn3], [1, 1, 1, 1], qncounts)
end


text:
"""
cglorot_uniform([rng=GLOBAL_RNG], dims...)

A modification of the `glorot_uniform` function provided by `Flux` to accommodate Complex numbers.
This is necessary since the parameters of the global convolution operator in the Fourier Layer generally has complex weights.
"""
cglorot_uniform(rng::AbstractRNG, dims...) = (rand(rng, ComplexF32, dims...) .- 0.5f0) .* sqrt(24.0f0 / sum(nfan(dims...)))
cglorot_uniform(dims...) = cglorot_uniform(Random.GLOBAL_RNG, dims...)
cglorot_uniform(rng::AbstractRNG) = (dims...) -> cglorot_uniform(rng, dims...)

"""
cglorot_normal([rng=GLOBAL_RNG], dims...)

A modification of the `glorot_normal` function provided by `Flux` to accommodate Complex numbers.
This is necessary since the parameters of the global convolution operator in the Fourier Layer generally has complex weights.
"""
cglorot_normal(rng::AbstractRNG, dims...) = randn(rng, ComplexF32, dims...) .* sqrt(2.0f0 / sum(nfan(dims...)))
cglorot_normal(dims...) = cglorot_normal(Random.GLOBAL_RNG, dims...)
cglorot_normal(rng::AbstractRNG) = (dims...) -> cglorot_normal(rng, dims...)