text:
macro repeattimes(nusers::Int, ntimes::Int, ex::Expr)
    return quote
        @sync for usr=1:$nusers
            @async for itr=1:$ntimes
                $ex
            end
        end
    end
end

macro repeatduring(nusers::Int, seconds::Int, ex::Expr)
    return quote
        @sync for usr=1:$nusers
            start = time()
            @async while time() - start < $seconds
                $ex
            end
        end
    end
end

@doc """
Run expression by specieid number of users (threads) specified number of times
in parallel and collect results.
Example - run GET request by 100 users 5000 times each:

    using Requests
    timetable = @runtimes 100 5000 get("http://example.com")

""" ->
macro runtimes(nusers::Int, ntimes::Int, ex::Expr)
    return quote
        ctr = Controller(COLLECTOR_PORT)
        @sync for usr=1:$nusers
            r = Reporter(COLLECTOR_PORT, ctr.key)
            @async for itr=1:$ntimes
                success = true
                errorbuf = IOBuffer()
                now = time()
                t = @elapsed try
                    $ex
                catch e
                    success = false
                    bt = catch_backtrace()
                    showerror(errorbuf, e, bt)
                end
                err = split(bytestring(errorbuf), "\n")[1]
                report(r, "$now,$usr,$itr,$success,$err,$t")
            end
        end
        info("Ready, collecting results...")
        finalize(ctr)
    end
end

@doc """
Run expression by specified number of users during specified time in parallel
and collect results.
Example - run GET request by 100 users during 30 seconds:

    using Requests
    timetable = @runduring 100 30 get("http://example.com")

""" ->
macro runduring(nusers::Int, seconds::Int, ex::Expr)
    return quote
        ctr = Controller(COLLECTOR_PORT)
        @sync for usr=1:$nusers
            r = Reporter(COLLECTOR_PORT, ctr.key)
            start = time()
            @async begin
                itr = 0
                while time() - start < $seconds
                    itr += 1
                    success = true
                    errorbuf = IOBuffer()
                    now = time()
                    t = @elapsed try
                        $ex
                    catch e
                        success = false
                        bt = catch_backtrace()
                        showerror(errorbuf, e, bt)
                    end
                    err = split(bytestring(errorbuf), "\n")[1]
                    report(r, "$now,$usr,$itr,$success,$err,$t")
                end
            end
        end
        info("Ready, collecting results...")
        finalize(ctr)
    end
end


text:
## constants.jl -- necessary constants for setting up

"""
    AminoAcidNames

A `Dict` of three-letter amino acid residue codes and their corresponding
one-letter codes.
"""
const AminoAcidNames = Dict(
    :ALA => :A, :CYS => :C,
    :ASP => :D, :GLU => :E,
    :PHE => :F, :GLY => :G,
    :HIS => :H, :ILE => :I,
    :LYS => :K, :LEU => :L,
    :MET => :M, :ASN => :N,
    :PRO => :P, :GLN => :Q,
    :ARG => :R, :SER => :S,
    :THR => :T, :VAL => :V,
    :TRP => :W, :TYR => :Y
)
const kAminoAcidNames = keys(AminoAcidNames)

"""
    Protstructs.DSSPConstants

Constants used for calculations in DSSP source code. (will I get in trouble for
sharing this?). This `Dict` defines the following:
+ Minimum acceptable H--bond interaction distance: ``0.5`` (angstroms?)
+ Minimum distance between alpha carbons: ``9.0`` (angstroms?)
+ Maximum N--C distance in a peptide bond: ``2.5`` (angstroms?)
+ HBond energy range: ``[-9.9, -0.5]`` (kcal/mol)
+ Coupling constant: ``-332 * 0.42 * 0.2`` (appropriate units)

**Note**: I may get in trouble for sharing this. I thought it would be fine
since it's on the github hosted version.
"""
const DSSPConstants = Dict(
    :minCADistance => 9.0,                   ## between alpha carbons
    :minHBondEnergy => -9.9,                 ## constant for checking
    :maxHBondEnergy => -0.5,                 ## another for checking
    :coupleConstant => -332 * 0.42 * 0.2,    ## used in actual calculation
    :maxNCDistance => 2.5                    ## check length of peptide bond
)

ext: jl


text:
function get_evt_range!(filestr::String, idr::Regex, t::Array{Int64,1})
  f = h5open(filestr, "r")
  D = get_datasets(f)
  for i in D
    has(i, "event_id") || continue
    evid = read(i["event_id"])
    if occursin(idr, evid)
      t0 = read(i["starttime"])
      fs = read(i["sampling_rate"])
      nx = length(i)

      # convert fs to sampling interval in ns
      Δ = round(Int64, 1.0e9/fs)
      t1 = t0 + (nx-1)*Δ

      t[1] = min(t[1], t0-Δ)
      t[2] = max(t[2], t1+Δ)
    end
  end
  close(f)
  return nothing
end

@doc """
    EventCat = read_asdf_evt(filestr, event_id::Union{String, Regex}[, KWs])

Read data in seismic HDF5 format with ids matching `event_id` from files
matching pattern `filestr`. Returns an array of SeisEvent structures.

    EventCat = read_asdf_evt(filestr, [, KWs])

Read data in seismic HDF5 format from files matching pattern `filestr` into
SeisEvent structures. Matches any event ID in any matching file.

|KW         | Type    | Default | Meaning                                     |
|:---       |:---     |:---     |:---                                         |
| msr       | Bool    | true    | read full (MultiStageResp) instrument resp? |
| v         | Int64   | 0       | verbosity                                   |

See also: timespec, parsetimewin, read_data, read_hdf5
""" read_asdf_evt
function read_asdf_evt(filestr::String, event_id::Union{String, Regex};
  msr       ::Bool                  = true,                   # read multistage response?
  v         ::Int64                 = KW.v                    # verbosity
  )

  one_file = safe_isfile(filestr)
  S = SeisData()

  # Regex ID
  idr = isa(event_id, String) ? id_to_regex(event_id) : id

  # Time range to read
  t = zeros(Int64, 2)
  t[1] = typemax(Int64)
  t[2] = typemin(Int64)
  if one_file

    # Time range
    get_evt_range!(filestr, idr, t)
    t_start = string(u2d(t[1]*1.0e-9))
    t_end = string(u2d(t[2]*1.0e-9))
    (v > 0) && println("read range: ", t_start, " -- ", t_end)

    # Read data
    append!(S, read_asdf(filestr, "*", t_start, t_end, msr, v))
  else
    files = ls(filestr)

    # Time range
    for fname in files
      get_evt_range!(fname, idr, t)
    end
    t_start = string(u2d(t[1]*1.0e-9))
    t_end = string(u2d(t[2]*1.0e-9))
    (v > 0) && println("read range: ", t_start, " -- ", t_end)

    # Read data
    for fname in files
      (v > 0) && println("reading from ", fname)
      append!(S, read_asdf(fname, "*", t_start, t_end, msr, v))
    end
  end

  # List of event IDs for each channel
  ns = S.n
  sid = Array{String,1}(undef, ns)
  for i in 1:ns
    sid[i] = get(S.misc[i], "event_id", "")
  end

  (H,R) = asdf_rqml(filestr)

  # List of event IDs for each SeisHdr
  nh = length(H)
  hid = Array{String,1}(undef, nh)
  inds = Int64[]
  sizehint!(inds, nh)
  fill!(hid, "")
  for i in 1:nh
    evid = H[i].id
    hid[i] = evid
    if occursin(event_id, evid)
      push!(inds, i)
    end
  end

  # Form event catalog by matching IDs
  EC = Array{SeisEvent,1}(undef, length(inds))
  for (k,j) in enumerate(inds)
    chans = Int64[]
    for i in 1:S.n
      if sid[i] == hid[j]
        push!(chans, i)
      end
    end
    EC[k] = SeisEvent(hdr = H[j], source = R[j], data = S[chans])
  end
  return EC
end

read_asdf_evt(filestr::String;
  msr       ::Bool                  = true,                   # read multistage response?
  v         ::Int64                 = KW.v                    # verbosity
  ) = read_asdf_evt(filestr, "", msr=msr, v=v)


text:
using Raylib

function main()

    screenWidth = 800
    screenHeight = 450

    Raylib.InitWindow(screenWidth, screenHeight,
                         "raylib [core] example - basic window")

    Raylib.SetTargetFPS(60)

    while !Raylib.WindowShouldClose()
        Raylib.BeginDrawing()
        Raylib.ClearBackground(Raylib.RAYWHITE)
        Raylib.DrawText("Congrats! You created your first window!",
                           190, 200, 20, Raylib.LIGHTGRAY)
        Raylib.EndDrawing()
    end

    Raylib.CloseWindow()
end


text:
# adapted from https://gist.github.com/Keno/d0c2df947f67be543036238a0caeb1c6

let counter = 0
	global get_fast_sin_counter, reset_counter
	fast_sin(x::Float64) = (counter += 1; sin(x::Float64))
	get_fast_sin_counter() = counter
	reset_counter() = (counter = 0)
end

get_sin() = sin

function my_func(x)
	reset_counter()
	sin(x)
	@test get_fast_sin_counter() == 1
	get_sin()(x)
	@test get_fast_sin_counter() == 2
	get_sin()(Base.inferencebarrier(x))
	@test get_fast_sin_counter() == 3
end

with(FastSinPlugin) do
	my_func(1.0)
end


text:
using DocOpt

doc="""Training Gaussian Mixture models for Voice Conversion.

Usage:
    train_gmm.jl [options] <parallel_dir> <dst_jld>
    train_gmm.jl --version
    train_gmm.jl -h | --help

Options:
    -h --help           show this message
    --diff              use differencial
    --add_delta         add delta feature
    --n_components=MIX  number of mixtures [default: 16]
    --n_iter=ITER       number of iterations [default: 200]
    --n_init=N          number of initialization [default: 2]
    --max=MAX           maximum number of data we use [default: 100]
    --min_covar=MIN     minimum covariance [default: 1.0e-7]
    --refine            refine pre-trained GMM (<dst_jld>)
"""

using VoiceConversion
using HDF5, JLD
using PyCall

# TODO: remove this dependency
@pyimport sklearn.mixture as mixture

using Logging
@Logging.configure(level=DEBUG, output=STDOUT)

# pygmmmean2jl return transposed parameters because julia is column-major
pygmmmean2jl(means::AbstractMatrix) = means'

# shape (M, D, D) -> (D, D, M)
# where M is the # of mixtures and D is the dimention of feature vector
function pygmmcovar2jl{T}(py_covars::Array{T,3})
    n_components = size(py_covars, 1)
    order = size(py_covars, 2)
    covars = Array(eltype(py_covars), order, order, n_components)
    for m=1:n_components
        covars[:,:,m] = reshape(py_covars[m,:,:], order, order)
    end
    covars
end

jlgmmmean2py(means::AbstractMatrix) = means'

# shape (D, D, M) -> (M, D, D)
function jlgmmcovar2py{T}(jl_covars::Array{T,3})
    n_components = size(jl_covars, 3)
    order = size(jl_covars, 2)
    covars = Array(eltype(jl_covars), n_components, order, order)
    for m=1:n_components
        covars[m,:,:] = reshape(jl_covars[:,:,m], 1, order, order)
    end
    covars
end

function copy2pygmm!(jlgmm, pygmm)
    pygmm[:means_] = jlgmmmean2py(jlgmm["means"])
    pygmm[:covars_] = jlgmmcovar2py(jlgmm["covars"])
    pygmm[:weights_] = jlgmm["weights"]
    pygmm[:init_params] = ""
end

let
    args = docopt(doc, version=v"0.0.2")
    diff = args["--diff"]
    n_components = parse(Int, args["--n_components"])
    n_iter = parse(Int, args["--n_iter"])
    n_init = parse(Int, args["--n_init"])
    min_covar = parse(Float64, args["--min_covar"])
    refine = args["--refine"]

    dataset = ParallelDataset(args["<parallel_dir>"],
                              joint=true,
                              diff=diff,
                              add_delta=args["--add_delta"],
                              nmax=parse(Int, args["--max"]))

    savepath = args["<dst_jld>"]

    gmm = mixture.GMM(n_components=n_components,
                      covariance_type="full",
                      n_iter=n_iter,
                      n_init=n_init,
                      min_covar=min_covar
                      )

    savepath = args["<dst_jld>"]
    if refine
        if !isfile(savepath)
            error("$(savepath) not found. Cannot refine model.")
        end
        println("refine pre-trained GMM")
        pretrained_jlgmm = load(savepath)
        copy2pygmm!(pretrained_jlgmm, gmm)
    end

    @show gmm
    # pass transposed matrix because python is row-major language
    elapsed = @elapsed gmm[:fit](dataset.X')
    @info("Elapsed time in training is $(elapsed) sec.")

    save(savepath,
         "description", "Parameters of Gaussian Mixture Model",
         "diff", diff,
         "n_components", length(gmm[:weights_]),
         "weights", gmm[:weights_],
         "means", pygmmmean2jl(gmm[:means_]),
         "covars", pygmmcovar2jl(gmm[:covars_])
         )

    @info("Dumped to $(savepath)")

    println("Finished")
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Area <: AggregationTopology
        name::String
        maxactivepower::Float64
        maxreactivepower::Float64
        load_response::Float64
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

A collection of buses for control purposes.

# Arguments
- `name::String`
- `maxactivepower::Float64`
- `maxreactivepower::Float64`
- `load_response::Float64`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct Area <: AggregationTopology
    name::String
    maxactivepower::Float64
    maxreactivepower::Float64
    load_response::Float64
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function Area(name, maxactivepower=0.0, maxreactivepower=0.0, load_response=0.0, forecasts=InfrastructureSystems.Forecasts(), )
    Area(name, maxactivepower, maxreactivepower, load_response, forecasts, InfrastructureSystemsInternal(), )
end

function Area(; name, maxactivepower=0.0, maxreactivepower=0.0, load_response=0.0, forecasts=InfrastructureSystems.Forecasts(), )
    Area(name, maxactivepower, maxreactivepower, load_response, forecasts, )
end

# Constructor for demo purposes; non-functional.
function Area(::Nothing)
    Area(;
        name="init",
        maxactivepower=0.0,
        maxreactivepower=0.0,
        load_response=0.0,
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::Area) = value.name
"""Get Area maxactivepower."""
get_maxactivepower(value::Area) = value.maxactivepower
"""Get Area maxreactivepower."""
get_maxreactivepower(value::Area) = value.maxreactivepower
"""Get Area load_response."""
get_load_response(value::Area) = value.load_response

InfrastructureSystems.get_forecasts(value::Area) = value.forecasts
"""Get Area internal."""
get_internal(value::Area) = value.internal

InfrastructureSystems.set_name!(value::Area, val::String) = value.name = val
"""Set Area maxactivepower."""
set_maxactivepower!(value::Area, val::Float64) = value.maxactivepower = val
"""Set Area maxreactivepower."""
set_maxreactivepower!(value::Area, val::Float64) = value.maxreactivepower = val
"""Set Area load_response."""
set_load_response!(value::Area, val::Float64) = value.load_response = val

InfrastructureSystems.set_forecasts!(value::Area, val::InfrastructureSystems.Forecasts) = value.forecasts = val
"""Set Area internal."""
set_internal!(value::Area, val::InfrastructureSystemsInternal) = value.internal = val


text:
__precompile__(true)

module BridgeDiffEq

using Reexport
@reexport using DiffEqBase

using StaticArrays
import Bridge

import DiffEqBase: solve

const warnkeywords =
    (:save_idxs, :d_discontinuities, :unstable_check, :save_everystep,
     :save_end, :initialize_save, :adaptive, :abstol, :reltol, :dtmax,
     :dtmin, :force_dtmin, :internalnorm, :gamma, :beta1, :beta2,
     :qmax, :qmin, :qsteady_min, :qsteady_max, :qoldinit, :failfactor,
     :maxiters, :isoutofdomain, :unstable_check,
     :calck, :progress, :timeseries_steps, :tstops, :saveat, :dense)

function __init__()
    const global warnlist = Set(warnkeywords)
end

include("algorithms.jl")
include("solve.jl")

export BridgeAlgorithm, BridgeEuler, BridgeHeun, BridgeSRK, BridgeR3, BridgeBS3

end # module


text:
using DataFrames
using Rimu
using LinearAlgebra
using SafeTestsets
using StaticArrays
using Statistics
using Suppressor
using Test
using Rimu.StatsTools, Rimu.RimuIO

# assuming VERSION ≥ v"1.6"
# the following is needed because random numbers of collections are computed
# differently after version 1.6, and thus the results of many tests change
# for Golden Master Testing (@https://en.wikipedia.org/wiki/Characterization_test)
@assert VERSION ≥ v"1.6"

@safetestset "Interfaces" begin
    include("Interfaces.jl")
end

@safetestset "StatsTools" begin
    include("StatsTools.jl")
end

@safetestset "BitStringAddresses" begin
    include("BitStringAddresses.jl")
end

using Rimu.ConsistentRNG
@testset "ConsistentRNG" begin
    seedCRNG!(127) # uses `RandomNumbers.Xorshifts.Xoshiro256StarStar()`
    @test cRand(UInt128) == 0x50f0f296b239b257a8c2ac2f11d6d2cb

    @test rand(ConsistentRNG.CRNGs[][1],UInt128) == 0xba97314c00e092e448993a2bef41d28d
    # Only looks at first element of the `NTuple`. This should be reproducible
    # regardless of `numthreads()`.
    @test rand(trng(),UInt16) == 0x4e65
    @test rand(newChildRNG(),UInt16) == 0x03aa
    @test ConsistentRNG.check_crng_independence(0) == Threads.nthreads()
end

@safetestset "StochasticStyles" begin
    include("StochasticStyles.jl")
end

@safetestset "DictVectors" begin
    include("DictVectors.jl")
end

@testset "Hamiltonians" begin
    include("Hamiltonians.jl")
end

@safetestset "lomc!" begin
    include("lomc.jl")
end

@testset "MemoryStrategy" begin
    # Define the initial Fock state with n particles and m modes
    n = m = 9
    add = near_uniform(BoseFS{n,m})
    H = HubbardReal1D(add; u = 6.0, t = 1.0)
    dv = DVec(add => 1; style=IsStochasticWithThreshold(1.0))
    s_strat = DoubleLogUpdate(targetwalkers=100)

    @testset "NoMemory" begin
        seedCRNG!(12345)
        df = lomc!(
            H, copy(dv);
            laststep=100, s_strat, m_strat=NoMemory(), maxlength=2*dimension(H)
        ).df
        @test sum(df[:,:norm]) ≈ 2698 atol=1
    end

    @testset "DeltaMemory" begin
        seedCRNG!(12345)
        df = lomc!(
            H, copy(dv);
            laststep=100, s_strat, m_strat=DeltaMemory(1), maxlength=2*dimension(H)
        ).df
        @test sum(df[:,:norm]) ≈ 2698 atol=1

        seedCRNG!(12345)
        df = lomc!(
            H, copy(dv);
            laststep=100, s_strat, m_strat=DeltaMemory(10), maxlength=2*dimension(H)
        ).df
        @test sum(df[:,:norm]) ≈ 2005 atol=1
    end

    @testset "DeltaMemory2" begin
        seedCRNG!(12345)
        df = lomc!(
            H, copy(dv);
            laststep=100, s_strat, m_strat=Rimu.DeltaMemory2(1), maxlength=2*dimension(H)
        ).df
        @test sum(df[:,:norm]) ≈ 2698 atol=1

        seedCRNG!(12345)
        df = lomc!(
            H, copy(dv);
            laststep=100, s_strat, m_strat=Rimu.DeltaMemory2(10), maxlength=2*dimension(H)
        ).df
        @test sum(df[:,:norm]) ≈ 1848 atol=1
    end

    @testset "ShiftMemory" begin
        seedCRNG!(12345)
        df = lomc!(
            H, copy(dv);
            laststep=100, s_strat, m_strat=ShiftMemory(1), maxlength=2*dimension(H)
        ).df
        @test sum(df[:,:norm]) ≈ 2698 atol=1

        seedCRNG!(12345)
        df = lomc!(
            H, copy(dv);
            laststep=100, s_strat, m_strat=ShiftMemory(10), maxlength=2*dimension(H)
        ).df
        @test sum(df[:,:norm]) ≈ 2811 atol=1
    end
end

@testset "IsDeterministic with Vector" begin
    ham = HubbardReal1D(BoseFS((1, 1, 1, 1)))
    dim = dimension(ham)
    sm, basis = Rimu.Hamiltonians.build_sparse_matrix_from_LO(ham, starting_address(ham))
    @test dim == length(basis)
    # run lomc! in deterministic mode with Matrix and Vector
    a = lomc!(sm, ones(dim); threading=true).df # no actual threading is done, though
    b = lomc!(sm, ones(dim); threading=false).df
    @test a.shift ≈ b.shift
    # run lomc! in deterministic mode with Hamiltonian and DVec
    v = DVec(k=>1.0 for k in basis) # corresponds to `ones(dim)`
    c = lomc!(ham, v).df
    @test a.shift ≈ c.shift
end

@testset "helpers" begin
    @testset "walkernumber" begin
        v = [1,2,3]
        @test walkernumber(v) == norm(v,1)
        dvc = DVec(:a => 2-5im)
        @test StochasticStyle(dvc) isa StochasticStyles.IsStochastic2Pop
        @test walkernumber(dvc) == 2.0 + 5.0im
        Rimu.purge_negative_walkers!(dvc)
        @test walkernumber(dvc) == 2.0 + 0.0im
        dvi= DVec(:a=>Complex{Int32}(2-5im))
        @test StochasticStyle(dvi) isa StochasticStyles.IsStochastic2Pop
        dvr = DVec(i => cRandn() for i in 1:100; capacity = 100)
        @test walkernumber(dvr) ≈ norm(dvr,1)
    end
    @testset "MultiScalar" begin
        a = Rimu.MultiScalar(1, 1.0, SVector(1))
        b = Rimu.MultiScalar(SVector(2, 3.0, SVector(4)))
        c = Rimu.MultiScalar((3, 4.0, SVector(5)))
        @test a + b == c
        @test_throws MethodError a + Rimu.MultiScalar(1, 1, 1)

        @test Rimu.combine_stats(a) == a
        @test Rimu.combine_stats([a, b]) == c
    end
end

using Rimu.Blocking
@testset "Blocking" begin
    n=10
    a = rand(n)
    m = mean(a)
    @test m == sum(a)/n
    myvar(a,m) = sum((a .- m).^2)/n
    @test var(a) == sum((a .- m).^2)/(n-1)
    @test var(a, corrected=false) == sum((a .- m).^2)/n == myvar(a,m)
    @test var(a, corrected=false) == var(a, corrected=false, mean = m)
    # @benchmark myvar($a, $m)
    # @benchmark var($a, corrected=false, mean = $m)
    # evaluating the above shows that the library function is faster and avoids
    # memory allocations completely

    # test se
    a = collect(1:10)
    @test Rimu.Blocking.se(a) ≈ 0.9574271077563381
    @test Rimu.Blocking.se(a;corrected=false) ≈ 0.9082951062292475
    # test autocovariance
    @test autocovariance(a,1) ≈ 6.416666666666667
    @test autocovariance(a,1;corrected=false) ≈ 5.775
    # test covariance
    b = collect(2:11)
    @test covariance(a,b) ≈ 9.166666666666666
    @test covariance(a,b;corrected=false) ≈ 8.25
    c = collect(2:20) # should be truncated
    @test covariance(a,b) == covariance(a,c)
    @test covariance(a,b;corrected=false) == covariance(a,c;corrected=false)

    # Define the initial Fock state with n particles and m modes
    n = m = 9
    aIni = near_uniform(BoseFS{n,m})
    ham = HubbardReal1D(aIni; u = 6.0, t = 1.0)
    pa = RunTillLastStep(laststep = 1000)

    # standard fciqmc
    s = DoubleLogUpdate(targetwalkers = 100)
    svec = DVec(Dict(aIni => 2))
    StochasticStyle(svec)
    vs = copy(svec)
    post_step = ProjectedEnergy(ham, svec)
    τ_strat = ConstantTimeStep()

    seedCRNG!(12345) # uses RandomNumbers.Xorshifts.Xoroshiro128Plus()
    # @time rdfs = fciqmc!(vs, pa, ham, s, r_strat, τ_strat, similar(vs))
    @time rdfs = lomc!(
        ham, vs; params = pa, s_strat = s, post_step, τ_strat, wm = similar(vs),
    ).df
    r = autoblock(rdfs, start=101)
    @test r.s̄ ≈ -5.36 atol=0.1
    @test r.σs ≈ 0.27 atol=0.1
    @test r.ē ≈ -7.46 atol=0.1
    @test r.σe ≈ 0.58 atol=0.1
    @test r.k == 6

    g = growthWitness(rdfs, b=50)
    # @test sum(g) ≈ -5725.3936298329545
    @test length(g) == nrow(rdfs)
    g = growthWitness(rdfs, b=50, pad = :false)
    @test length(g) == nrow(rdfs) - 50
    @test_throws AssertionError growthWitness(rdfs.norm, rdfs.shift[1:end-1],rdfs.dτ[1])
end

@testset "RimuIO" begin
    @testset "save_df, load_df" begin
        file = joinpath(@__DIR__, "tmp.arrow")
        rm(file; force=true)

        df = DataFrame(a=[1, 2, 3], b=Complex{Float64}[1, 2, 3+im], d=rand(Complex{Int}, 3))
        RimuIO.save_df(file, df)
        df2 = RimuIO.load_df(file)
        @test df == df2

        rm(file)
    end
    @testset "save_dvec, load_dvec" begin
        # BSON is currently broken on 1.8
        if VERSION ≤ v"1.7"
            file1 = joinpath(@__DIR__, "tmp1.bson")
            file2 = joinpath(@__DIR__, "tmp2.bson")
            rm(file1; force=true)
            rm(file2; force=true)

            add = BoseFS2C((1,1,0,1), (1,1,0,0))
            dv = InitiatorDVec(add => 1.0, style=IsDynamicSemistochastic(abs_threshold=3.5))
            H = BoseHubbardMom1D2C(add)

            _, state = lomc!(H, dv; replica=NoStats(2))
            RimuIO.save_dvec(file1, state.replicas[1].v)
            RimuIO.save_dvec(file2, state.replicas[2].v)

            dv1 = RimuIO.load_dvec(file1)
            dv2 = RimuIO.load_dvec(file2)

            @test dv1 == state.replicas[1].v
            @test typeof(dv2) == typeof(state.replicas[1].v)
            @test StochasticStyle(dv1) == StochasticStyle(state.replicas[1].v)
            @test storage(dv2) == storage(state.replicas[2].v)

            rm(file1; force=true)
            rm(file2; force=true)
        end
    end
end

using Rimu.EmbarrassinglyDistributed
@testset "EmbarrassinglyDistributed" begin
    add = BoseFS((1,1,0,1))
    v = DVec(add => 2; capacity = 200)
    ham = HubbardReal1D(add, u=4.0)
    @test setup_workers(4) == 4 # add workers and load code (Rimu and its modules)
    seedCRNGs_workers!(127)     # seed rgns on workers deterministically
    nt = d_lomc!(ham, v; eqsteps = 1_000, laststep = 21_000) # perform parallel lomc!
    @test [size(df)[1] for df in nt.dfs] == [6000, 6000, 6000, 6000]
    ntc = combine_dfs(nt) # combine results into one DataFrame
    @test size(ntc.df)[1] == 20997
    energies = autoblock(ntc) # perform `autoblock()` discarding `eqsteps` time steps
    # in a single line:
    # energies = d_lomc!(ham, v; eqsteps = 1_000, laststep = 21_000) |> combine_dfs |> autoblock
    @test ismissing(energies.ē) && ismissing(energies.σe)
    # golden master test on results because qmc evolution is deterministic
    @test energies.s̄ ≈ -4.1 atol=0.1
    @test energies.σs ≈ 0.006 atol=1e-3
end

@testset "BoseFS2C" begin
    bfs2c = BoseFS2C(BoseFS((1,2,0,4)),BoseFS((4,0,3,1)))
    @test typeof(bfs2c) <: BoseFS2C{7,8,4}
    @test num_occupied_modes(bfs2c.bsa) == 3
    @test num_occupied_modes(bfs2c.bsb) == 3
    @test onr(bfs2c.bsa) == [1,2,0,4]
    @test onr(bfs2c.bsb) == [4,0,3,1]
    @test Hamiltonians.bose_hubbard_2c_interaction(bfs2c) == 8 # n_a*n_b over all sites
end

@testset "TwoComponentBosonicHamiltonian" begin
    aIni2cReal = BoseFS2C(BoseFS((1,1,1,1)),BoseFS((1,1,1,1))) # real space two-component
    Ĥ2cReal = BoseHubbardReal1D2C(aIni2cReal; ua = 6.0, ub = 6.0, ta = 1.0, tb = 1.0, v= 6.0)
    hamA = HubbardReal1D(BoseFS((1,1,1,1)); u=6.0, t=1.0)
    hamB = HubbardReal1D(BoseFS((1,1,1,1)); u=6.0)
    @test hamA == Ĥ2cReal.ha
    @test hamB == Ĥ2cReal.hb
    @test num_offdiagonals(Ĥ2cReal,aIni2cReal) == 16
    @test num_offdiagonals(Ĥ2cReal,aIni2cReal) == num_offdiagonals(Ĥ2cReal.ha,aIni2cReal.bsa)+num_offdiagonals(Ĥ2cReal.hb,aIni2cReal.bsb)
    @test dimension(Ĥ2cReal) == 1225
    @test dimension(Float64, Ĥ2cReal) == 1225.0

    hp2c = offdiagonals(Ĥ2cReal,aIni2cReal)
    @test length(hp2c) == 16
    @test hp2c[1][1] == BoseFS2C(BoseFS((0,2,1,1)), BoseFS((1,1,1,1)))
    @test hp2c[1][2] ≈ -1.4142135623730951
    @test diagonal_element(Ĥ2cReal,aIni2cReal) ≈ 24.0 # from the V term

    aIni2cMom = BoseFS2C(BoseFS((0,4,0,0)),BoseFS((0,4,0,0))) # momentum space two-component
    Ĥ2cMom = BoseHubbardMom1D2C(aIni2cMom; ua = 6.0, ub = 6.0, ta = 1.0, tb = 1.0, v= 6.0)
    @test num_offdiagonals(Ĥ2cMom,aIni2cMom) == 9
    @test dimension(Ĥ2cMom) == 1225
    @test dimension(Float64, Ĥ2cMom) == 1225.0

    hp2cMom = offdiagonals(Ĥ2cMom,aIni2cMom)
    @test length(hp2cMom) == 9
    @test hp2cMom[1][1] == BoseFS2C(BoseFS((1,2,1,0)), BoseFS((0,4,0,0)))
    @test hp2cMom[1][2] ≈ 2.598076211353316

    smat2cReal, adds2cReal = Hamiltonians.build_sparse_matrix_from_LO(Ĥ2cReal,aIni2cReal)
    eig2cReal = eigen(Matrix(smat2cReal))
    smat2cMom, adds2cMom = Hamiltonians.build_sparse_matrix_from_LO(Ĥ2cMom,aIni2cMom)
    eig2cMom = eigen(Matrix(smat2cMom))
    @test eig2cReal.values[1] ≈ eig2cMom.values[1]
end

@safetestset "KrylovKit" begin
    include("KrylovKit.jl")
end
@safetestset "RMPI" begin
    include("RMPI.jl")
end

@testset "example script" begin
    include("../scripts/BHM-example.jl")
    dfr = load_df("fciqmcdata.arrow")
    qmcdata = last(dfr,steps_measure)
    (qmcShift,qmcShiftErr) = mean_and_se(qmcdata.shift)
    @test qmcShift ≈ -4.11255936332424 rtol=0.01

    # clean up
    rm("fciqmcdata.arrow", force=true)
end

@safetestset "doctests" begin
    include("doctests.jl")
end

# Note: This test is only for local testing, as MPI is tested separately on CI
@testset "MPI" begin
    # read name of mpi executable from environment variable if defined
    # necessary for allow-run-as root workaround for Pipelines
    mpiexec = haskey(ENV, "JULIA_MPIEXEC") ? ENV["JULIA_MPIEXEC"] : "mpirun"
    is_local = !haskey(ENV, "CI")

    juliaexec = Base.julia_cmd()

    if is_local
        mpi_test_filename = isfile("mpi_runtests.jl") ?  "mpi_runtests.jl" : "test/mpi_runtests.jl"
        if isfile(mpi_test_filename)
            rr = run(`$mpiexec -np 2 $juliaexec -t 1 $mpi_test_filename`)
            @test rr.exitcode == 0
        else
            @warn "Could not find mpi_runtests.jl. Not testing MPI."
        end
    else
        @info "not testing MPI on CI"
    end
end


text:
# constructs a semantic theory about a commutative monoid
# A monoid whose operation is commutative is called a
# commutative monoid (or, less commonly, an abelian monoid).

include("../docstrings.jl")

module Library

@info "NOTE: The current implementation of the Metatheory library currently works
correctly only with the EGraphs backend."

include("algebra.jl")

using ..Metatheory
# include("../util.jl")
# include("../rule.jl")

# theory generation macros
export commutativity
export associativity
export identity_left
export identity_right
export distrib_left
export distrib_right
export distrib
export monoid
export commutative_monoid
export commutative_group
export abelian_group
export @commutative_monoid
export @commutative_group
export @abelian_group

end


text:
"""
    SkewedExponentialPower(μ, σ, p, α)

The *Skewed exponential power distribution*, with location `μ`, scale `σ`, shape `p`, and skewness `α`
has the probability density function [1]
```math
f(x; \\mu, \\sigma, p, \\alpha) =
\\begin{cases}
\\frac{1}{\\sigma 2p^{1/p}\\Gamma(1+1/p)} \\exp \\left\\{ - \\frac{1}{2p}\\Big| \\frac{x-\\mu}{\\alpha \\sigma} \\Big|^p \\right\\}, & \\text{if } x \\leq \\mu \\\\
\\frac{1}{\\sigma 2p^{1/p}\\Gamma(1+1/p)} \\exp \\left\\{ - \\frac{1}{2p}\\Big| \\frac{x-\\mu}{(1-\\alpha) \\sigma} \\Big|^p \\right\\}, & \\text{if } x > \\mu
\\end{cases}.
```
The Skewed exponential power distribution (SEPD) incorporates the Laplace (``p=1, \\alpha=0.5``),
normal (``p=2, \\alpha=0.5``), uniform (``p\\rightarrow \\infty, \\alpha=0.5``), asymmetric Laplace (``p=1``), skew normal (``p=2``),
and exponential power distribution (``\\alpha = 0.5``) as special cases.

[1] Zhy, D. and V. Zinde-Walsh (2009). Properties and estimation of asymmetric exponential power distribution. _Journal of econometrics_, 148(1):86-96, 2009.

```julia
SkewedExponentialPower()            # SEPD with shape 2, scale 1, location 0, and skewness 0.5 (the standard normal distribution)
SkewedExponentialPower(μ, σ, p, α)  # SEPD with location μ, scale σ, shape p, and skewness α
SkewedExponentialPower(μ, σ, p)     # SEPD with location μ, scale σ, shape p, and skewness 0.5 (the exponential power distribution)
SkewedExponentialPower(μ, σ)        # SEPD with location μ, scale σ, shape 2, and skewness 0.5 (the normal distribution)
SkewedExponentialPower(μ)           # SEPD with location μ, scale 1, shape 2, and skewness 0.5 (the normal distribution)

params(d)       # Get the parameters, i.e. (μ, σ, p, α)
shape(d)        # Get the shape parameter, i.e. p
location(d)     # Get the location parameter, i.e. μ
scale(d)        # Get the scale parameter, i.e. σ
```
"""
struct SkewedExponentialPower{T <: Real} <: ContinuousUnivariateDistribution
    μ::T
    σ::T
    p::T
    α::T
    SkewedExponentialPower{T}(μ::T, σ::T, p::T, α::T) where {T} = new{T}(μ, σ, p, α)
end

function SkewedExponentialPower(µ::T, σ::T, p::T, α::T; check_args::Bool=true) where {T <: Real}
    if check_args
        @check_args(SkewedExponentialPower, σ > zero(σ))
        @check_args(SkewedExponentialPower, p > zero(p))
        @check_args(SkewedExponentialPower, zero(α) < α < one(α))
    end
    return SkewedExponentialPower{T}(µ, σ, p, α)
end

function SkewedExponentialPower(μ::Real, σ::Real, p::Real=2, α::Real=1//2; check_args::Bool=true)
    return SkewedExponentialPower(promote(μ, σ, p, α)...; check_args=check_args)
end
SkewedExponentialPower(μ::Real=0) = SkewedExponentialPower(μ, 1, 2, 1//2; check_args=false)

@distr_support SkewedExponentialPower -Inf Inf

### Conversions
convert(::Type{SkewedExponentialPower{T}}, μ::S, σ::S, p::S, α::S) where {T <: Real, S <: Real} = SkewedExponentialPower(T(μ), T(σ), T(p), T(α))
convert(::Type{SkewedExponentialPower{T}}, d::SkewedExponentialPower{S}) where {T <: Real, S <: Real} = SkewedExponentialPower(T(d.μ), T(d.σ), T(d.p), T(d.α), check_args=false)
convert(::Type{SkewedExponentialPower{T}}, d::SkewedExponentialPower{T}) where {T<:Real} = d

### Parameters
@inline partype(d::SkewedExponentialPower{T}) where {T<:Real} = T

params(d::SkewedExponentialPower) = (d.μ, d.σ, d.p, d.α)
location(d::SkewedExponentialPower) = d.μ
shape(d::SkewedExponentialPower) = d.p
scale(d::SkewedExponentialPower) = d.σ

### Statistics

#Calculates the kth central moment of the SEPD
function m_k(d::SkewedExponentialPower, k::Integer)
    _, σ, p, α = params(d)
    inv_p = inv(p)
    return  k * (logtwo + inv_p * log(p) + log(σ)) + loggamma((1 + k) * inv_p) -
        loggamma(inv_p) + log(abs((-1)^k * α^(1 + k) + (1 - α)^(1 + k)))
end

# needed for odd moments on log-scale
sgn(d::SkewedExponentialPower) = d.α > 1//2 ? -1 : 1

mean(d::SkewedExponentialPower) = d.α == 1//2 ? float(d.μ) : sgn(d)*exp(m_k(d, 1)) + d.μ
mode(d::SkewedExponentialPower) =  mean(d)
var(d::SkewedExponentialPower) = exp(m_k(d, 2)) - exp(2*m_k(d, 1))
skewness(d::SkewedExponentialPower) = d.α == 1//2 ? float(zero(partype(d))) : sgn(d)*exp(m_k(d, 3)) / (std(d))^3
kurtosis(d::SkewedExponentialPower) = exp(m_k(d, 4))/var(d)^2 - 3

function logpdf(d::SkewedExponentialPower, x::Real)
    μ, σ, p, α = params(d)
    a = x < μ ? α : 1 - α
    inv_p = inv(p)
    return -(logtwo + log(σ) + inv_p * log(p) + loggamma(1 + inv_p) + inv_p * (abs(μ - x) / (2 * σ * a))^p)
end

function cdf(d::SkewedExponentialPower, x::Real)
    μ, σ, p, α = params(d)
    inv_p = inv(p)
    if x <= μ
        α * ccdf(Gamma(inv_p), inv_p * (abs((x-μ)/σ) / (2*α))^p)
    else
        α + (1-α) * cdf(Gamma(inv_p), inv_p * (abs((x-μ)/σ) / (2*(1-α)))^p)
    end
end

function quantile(d::SkewedExponentialPower, p::Real)
    μ, σ, _, α = params(d)
    inv_p = inv(d.p)
    if p <= α
        μ - 2*α*σ * (d.p * quantile(Gamma(inv_p), 1-p/α))^(inv_p)
    else
        μ + 2*(1-α)*σ * (d.p * quantile(Gamma(inv_p), 1-(1-p)/(1-α)))^(inv_p)
    end
end

function rand(rng::AbstractRNG, d::SkewedExponentialPower)
    μ, σ, p, α = params(d)
    inv_p = inv(d.p)
    if rand(rng) < d.α
        μ - σ * 2*p^(inv_p) * α * rand(Gamma(inv_p, 1))^(inv_p)
    else
        μ + σ * 2*p^(inv_p) * (1-α) * rand(Gamma(inv_p, 1))^(inv_p)
    end
end


text:
function StokesPcCG_gpu(U, P, KK, MM, GG, Rhs, ifree)

    rtol_Pat = 1e-8
    itmax_Pat = 200
    itmin_Pat = 10
    itnum = 0

    nU = length(Rhs)

    # allocate some stuff
    Ud = CUDA.zeros(nU)
    z = CUDA.zeros(nU)
    Pd = CUDA.zeros(length(P))

    # move FEM matrices to device
    KKd = CuSparseMatrixCSC{Float32}(KK[ifree,ifree])
    GGd = CuSparseMatrixCSC{Float32}(GG)
    GGtransp = CuSparseMatrixCSC{Float32}(sparse(GG'))
    negGGtransp = CuSparseMatrixCSC{Float32}(sparse(-GG'))
    Rhsd = CuArray{Float32}(Rhs)

    # guess of the velocity field
    sol, = Krylov.cg(KKd, Rhsd[ifree])
    Ud[ifree] .= sol

    # initial pressure residual vector
    r = negGGtransp * Ud
    rm0space!(r)
    Prms = norm(r) # norm of r_i
    tol = Prms * rtol_Pat
    # get preconditioner
    MM, pc = _preconditioner(:jacobi, MM)
    MMd = CuSparseMatrixCSC{Float32}(MM)
    pcd = CuArray{Float32}(pc)

    # Begin of Patera pressure iterations
    # d = _precondition(pc, MM, r) # precondition residual
    d = _precondition(pcd, r) # precondition residual
    q = deepcopy(d)  # define FIRST search direction q
    rlast = similar(r)
    for itPat = 1:itmax_Pat
        itnum +=1

        rm0space!(q)
        rd = dot(r,d) # numerator for alpha

        #= Perform the S times q multiplication ======================
            S cannot be calculated explicitly, since Kinv cannot be formed
            Sq = S * q = (Gᵀ * Kinv * G) * q
            Hence, the muliplicatipon is done in 3 steps:
            (1) y   = G*q
            (2) K z = y
            (3) Sq  = Gᵀ*z
        =============================================================#
        y = GGd * q # (1) y = G*q
        sol, = Krylov.cg(KKd, y[ifree])
        z[ifree] .= sol
        # _CholeskyWithFactorization!(z, F, y, ifree) # (2) Solve K z = y
        Sq = GGtransp*z # (3) Sq  = G'*z
        #=============================================================#
        qSq = dot(q, Sq) # denominator to calculate alpha
        copyto!(rlast, r) # needed for Polak-Ribiere version 1
        α = rd/qSq # steps size in direction q

        # Update solution and residual
        Pd .+= α .* q
        Ud .+= α .* z
        r  .-= α .* Sq

        # remove nullspace
        rm0space!(r)

        # Check convergence
        # if norm(r) < tol && itPat>=itmin_Pat
        #     println("\n", itnum," CG iterations\n")
        #     break
        # end

        d  .= _precondition(pcd, r) # precondition residual
        # Make new search direction q S-orthogonal to all previous q's
        β  = dot(r-rlast, d)/rd # Polak-Ribiere version 1
        q .= xpy(q, d, β)

    end

    return Array(Ud), Array(Pd)

end

# Function takes out the mean of vector a
# (thereby removes the constant pressure mode, i.e. the nullspace)
function rm0space!(a::CuArray)
    m = sum(a)/length(a)
    a .=  a.-m
end

# PRECONDITIONING: JACOBI ITERATIONS
@inline function _precondition(pcd::CuArray, MMd::CuSparseMatrixCSC, r::CuArray)
    weight = (0.2500, 0.4000, 0.7500)
    # d = r .* pcd
    d = r .* pcd .+ weight[1] .* pcd .* (r .- MMd*d)
    d .+= weight[2] .* pcd .* (r .- MMd*d)
    d .+= weight[3] .* pcd .* (r .- MMd*d)
    return d
end

@inline _precondition(pc::CuArray, r::CuArray) = r .* pc

@inline function _updatesolution_gpu!(P,U,r,q,z,Sq,α)
    n1 = length(P)
    n2 = length(U)
    idx1 = 1:n1
    idx2 = (n1+1):n2
    @tturbo for i in idx1
        P[i] += α * q[i]
        U[i] += α * z[i]
        r[i] -= α * Sq[i]
    end
    @tturbo for i in idx2
        U[i] += α * z[i]
    end
end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


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

# ╔═╡ 4d9f28c0-6f82-11eb-1df0-f97522f733b1
begin
	#import Pkg;
	#Pkg.add("PlutoUI")
	#Pkg.add("Random")
	using PlutoUI
	using Random
end

# ╔═╡ 1d4cdc00-706b-11eb-3b0d-f7e15cd4ea6b
md"""
## Many bind objects in a markdown interpolation: **multiparameter function does not work**
"""

# ╔═╡ 17f84e80-6f83-11eb-0a6c-133363e073d6
md"Choose seed 👉 $(@bind seed Slider(1:1000, default = 1, show_value = true)) for random numbers."

# ╔═╡ dd4da650-6f81-11eb-2565-f5fbd7f13c8c
md"Choose number of 👉 $(@bind row_slider Slider(1:10,default=2, show_value = true)) rows of random matrix."

# ╔═╡ 4bbfe170-6f82-11eb-10eb-89f060f3148f
md"Choose number of 👉 $(@bind columns_slider Slider(1:10,default=1, show_value = true)) columns of random matrix."

# ╔═╡ 1c6417a0-6f84-11eb-3bdd-c1dbff2b5c63
md"""
Choose the multipliaction factor 👉 $(@bind mult_fact Select(["1"=>"One","10"=>"Ten","100"=>"Hundred","1000"=>"Thousand"]))
"""

# ╔═╡ af0fec20-6f82-11eb-385a-0d297bfc5ba4
A = round.(Int,tryparse(Int, mult_fact)*rand(MersenneTwister(seed), row_slider, columns_slider))

# ╔═╡ 3e6b0f20-706a-11eb-2eb6-e5a0b72f57dc
md"""
## Many bind objects with a **workaround of local variables in different cells** (Using a seed for the random variable)
"""

# ╔═╡ 7b626010-72d4-11eb-1f7f-0d63820d8b02
begin
	local bond_1_1 = @bind seed_2 Slider(1:1000, default = 1, show_value = true)
	md"""Choose seed 👉 $(bond_1_1) for random numbers.
"""
end

# ╔═╡ 9cbb5cd0-72d4-11eb-0e42-e79af8807c2d
begin
local bond_1 = @bind row_slider_2 Slider(1:10,default=2, show_value = true)
	md"""
Choose number of 👉 $(bond_1) rows of random matrix.
	"""
end

# ╔═╡ 934dd280-7069-11eb-29cd-2156919e9c30
begin
local bond_1 = @bind columns_slider_2 Slider(1:10,default=1, show_value = true)
md"""
Choose number of 👉 $(bond_1) columns of random matrix.
"""
end

# ╔═╡ d4c62240-72d4-11eb-0a92-15bd1bcbb249
begin
local bond_1 = @bind mult_fact_2 Select(["1"=>"One","10"=>"Ten","100"=>"Hundred","1000"=>"Thousand"])
md"""
Choose the multipliaction factor 👉 $(bond_1)
"""
end

# ╔═╡ 31eafc60-706a-11eb-0e11-6bac82af32fa
B = round.(Int,tryparse(Int, mult_fact_2)*rand(MersenneTwister(seed_2), row_slider_2, columns_slider_2))

# ╔═╡ 05c8a820-706b-11eb-3329-138b145aa667
md"""
## **Using rand():**
"""

# ╔═╡ 595bba00-706a-11eb-2ff7-e36aeec39df7
begin
local bond_1 = @bind row_slider_3 Slider(1:10,default=2, show_value = true)
local bond_2 = @bind columns_slider_3 Slider(1:10,default=1, show_value = true)
local bond_3 = @bind mult_fact_3 Select(["1"=>"One","10"=>"Ten","100"=>"Hundred","1000"=>"Thousand"])
md"""
Choose number of 👉 $(bond_1) rows of random matrix.

Choose number of 👉 $(bond_2) columns of random matrix.

Choose the multipliaction factor 👉 $(bond_3)
"""
end

# ╔═╡ c7a10e30-7069-11eb-1b2d-e5e5f3846ac3
C = round.(Int,tryparse(Int, mult_fact_3)*rand(row_slider_3, columns_slider_3))

# ╔═╡ Cell order:
# ╟─4d9f28c0-6f82-11eb-1df0-f97522f733b1
# ╟─1d4cdc00-706b-11eb-3b0d-f7e15cd4ea6b
# ╟─17f84e80-6f83-11eb-0a6c-133363e073d6
# ╟─dd4da650-6f81-11eb-2565-f5fbd7f13c8c
# ╟─4bbfe170-6f82-11eb-10eb-89f060f3148f
# ╟─1c6417a0-6f84-11eb-3bdd-c1dbff2b5c63
# ╟─af0fec20-6f82-11eb-385a-0d297bfc5ba4
# ╟─3e6b0f20-706a-11eb-2eb6-e5a0b72f57dc
# ╟─7b626010-72d4-11eb-1f7f-0d63820d8b02
# ╟─9cbb5cd0-72d4-11eb-0e42-e79af8807c2d
# ╟─934dd280-7069-11eb-29cd-2156919e9c30
# ╟─d4c62240-72d4-11eb-0a92-15bd1bcbb249
# ╟─31eafc60-706a-11eb-0e11-6bac82af32fa
# ╟─05c8a820-706b-11eb-3329-138b145aa667
# ╟─595bba00-706a-11eb-2ff7-e36aeec39df7
# ╟─c7a10e30-7069-11eb-1b2d-e5e5f3846ac3

max_stars_repo_name: fonsp/bayes-imoox

max_issues_repo_name: fonsp/bayes-imoox

max_forks_repo_name: fonsp/bayes-imoox


text:
using CartesianGrids
using Test
##using TestSetExtensions

#@test isempty(detect_ambiguities(ViscousFlow))
include("fields.jl")
include("points.jl")
include("generatedfields.jl")

#@testset ExtendedTestSet "All tests" begin
#    @includetests ARGS
#end

#if isempty(ARGS)
#    include("../docs/make.jl")
#end


text:
#=
This file contains code for separating a timeseries into its seasonal and residual
fluctuations. It interfaces SignalSeparation.jl and also uses some analytic knowledge
of the form of the incoming solar radiation
=#
using SignalDecomposition

export decompose, Fourier, MeanAlbedo, ProductInversion, fielddecompose, Sinusoidal
export nrmse, rmse, TimeAnomaly

# my own decomposition method that dispatches on dimensional array for input method

function fielddecompose(A::AbDimArray, method::Decomposition)
    t = dims(A, Time)
    fielddecompose(t, A, method)
end
function fielddecompose(A::AbDimArray, method::Fourier)
    t = 1:size(A, Time)
    fielddecompose(t, A, method)
end

function fielddecompose(t::AbstractVector, A::AbDimArray, method::Decomposition)
    seasonal = DimensionalArray(copy(Array(A)), dims(A), A.name*"x")
    residual = DimensionalArray(copy(Array(A)), dims(A), A.name*"r")
    for i in spatialidxs(A)
        y = Array(A[i...])
        sea, res = SignalDecomposition.decompose(t, y, method)
        seasonal[i...] .= sea
        residual[i...] .= res
    end
    return seasonal, residual
end

SignalDecomposition.decompose(t::Vector{<:TimeType}, s, m::Fourier) =
decompose(1:length(t), s, m)

SignalDecomposition.decompose(s::AbDimArray, m::Sinusoidal) =
decompose(dims(s, Time).val, Array(s), m)

function SignalDecomposition.decompose(t::AbstractVector{<:TimeType}, s, m::Sinusoidal)
    truetime = realtime_days(t)
    decompose(truetime, s, m)
end

#########################################################################
# Mean Albedo
#########################################################################
"""
    MeanAlbedo(S::Vector) <: Decomposition
    MeanAlbedo(t::Vector, φ::Real)
Decomposition that splits a timeseries `R` into its mean albedo component
``\\bar{\\alpha}S`` with ``S`` the Solar
insolation (either given directly as `S` or calculated analytically with given time
`t` and latitude `φ`) and ``\\bar{\\alpha} = `` `timeaverage(R)/timeaverage(S)`,
and the residuals of that.

It is supposed to approximately purely internal dynamical fluctuations.
"""
struct MeanAlbedo{T<:Real} <: Decomposition
    solar::Vector{T}
end
MeanAlbedo(t::Vector, φ::Real) = MeanAlbedo(insolation.(t, φ))
MeanAlbedo(x) = MeanAlbedo(Array(x))

# Extention for convenient usage in decompose below
ClimateBase.timeagg(f, t::Vector, a::AbDimArray{<:Any, 1}) = timeagg(f, t, a.data)

function SignalDecomposition.decompose(t, x, method::MeanAlbedo)
    @assert length(x) == length(method.solar)
    if eltype(t) <: TimeType
        meanα = timeagg(mean, t, x)/timeagg(mean, t, method.solar) # average albedo
    else
        meanα = mean(x)/mean(method.solar)
    end
    seasonal = meanα .* method.solar
    return seasonal, x .- seasonal
end

function meanalbedo_decompose(A::AbDimArray, S)
    solar = DimensionalArray(copy(Array(A)), dims(A), "MeanAlbedo αS")
    residual = DimensionalArray(copy(Array(A)), dims(A), "MeanAlbedo r")
    for i in spatialidxs(A)
        y = Array(A[i...])
        s = Array(S[i...])
        method = MeanAlbedo(s)
        sea, res = SignalDecomposition.decompose(y, method)
        solar[i...] .= sea
        residual[i...] .= res
    end
    return solar, residual
end

max_stars_repo_name: Datseris/EarthAlbedoSymmetry

max_issues_repo_name: Datseris/EarthAlbedoSymmetry

max_forks_repo_name: Datseris/EarthAlbedoSymmetry


text:
# select a list of hard case, perturb the demand, and then test
# those cases are selected as the two solvers' pg are very different from each other >= 1%

# compare the solvers for all cases
using Distributed;

addprocs(20);
@everywhere using PowerModels, Gurobi, Ipopt, Random, Distributions;
@everywhere using JLD,HDF5;
@everywhere include("solver_run.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();

list_file = readdir("./data");
result_Gurobi = Dict();
result_Ipopt = Dict();

itema = "case5.m";
network_data = parse_file("./data/"*itema);
load_keys = [i for i in keys(network_data["load"])];
σ_load = Dict()
load_distrn = Dict()
for i in 1:length(load_keys)
    if network_data["load"][load_keys[i]]["pd"] != 0
        σ_load[load_keys[i]] = abs(network_data["load"][load_keys[i]]["pd"]*0.03)
        load_distrn[load_keys[i]] = Normal(0,σ_load[load_keys[i]])
    else
        σ_load[load_keys[i]] = 0;
        load_distrn[load_keys[i]] = Uniform(0,0);
    end
end
pmap(ω -> run_diffd_para(network_data, load_distrn, load_keys), 1:20);

list_file_hard = ["case1888rte.m"
                 "case6515rte.m"
                 "case89pegase.m"
                 "case1354pegase.m"
                 "pglib_opf_case2746wp_k__api.m"
                 "pglib_opf_case24_ieee_rts__api.m"
                 "case6470rte.m"
                 "case6468rte.m"
                 "case1951rte.m"
                 "pglib_opf_case2746wop_k__api.m"
                 "pglib_opf_case3012wp_k__api.m"
                 "case13659pegase.m"
                 "case2848rte.m"
                 "case30pwl.m"
                 "pglib_opf_case2383wp_k__api.m"
                 "pglib_opf_case3375wp_k__api.m"
                 "case2869pegase.m"
                 "case6495rte.m"
                 "pglib_opf_case30000_goc__api.m"
                 "case9241pegase.m"
                 "pglib_opf_case60_c.m"
                 "case2868rte.m"
                  ];
perturbResult_dict = Dict();

for itema in list_file_hard
    if occursin(".m",itema)
        network_data = parse_file("./data/"*itema);
        load_keys = [i for i in keys(network_data["load"])]
        σ_load = Dict()
        load_distrn = Dict()
        for i in 1:length(load_keys)
            if network_data["load"][load_keys[i]]["pd"] != 0
                σ_load[load_keys[i]] = abs(network_data["load"][load_keys[i]]["pd"]*0.03)
                load_distrn[load_keys[i]] = Normal(0,σ_load[load_keys[i]])
            else
                σ_load[load_keys[i]] = 0;
                load_distrn[load_keys[i]] = Uniform(0,0.0000001);
            end
        end
        perturbResult_dict[itema] = pmap(ω -> run_diffd_para(network_data, load_distrn, load_keys), 1:1000);
        save("stochastic.jld","results",perturbResult_dict);
    end
end


text:
function solve()
    n, m, x = [parse(Int, x) for x in split(readline())]
    c = zeros(Int, n)
    a = zeros(Int, n, m)
    for i=1:n
        l = [parse(Int, x) for x in split(readline())]
        c[i] = l[1]
        a[i, :] = l[2:end]
    end

    min_p = sum(c) + 1
    for v=0:2^n-1
        p = 0
        r = zeros(Int, m)
        for i=0:n-1
            if v & (1 << i) != 0
                p += c[i+1]
                r += a[i+1, :]
            end
        end
        if minimum(r) >= x
            min_p = min(min_p, p)
        end
    end
    min_p > sum(c) ? -1 : min_p
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# auxiliary function to perform timings for action
function time_kernels!(
    to :: TimerOutput,
    r  :: Reduced_lattice,
    m  :: Mesh,
    a  :: Action,
    ap :: Action
    )  :: Nothing

    # get timer parameters
    Λ   = rand()
    v   = rand()
    dv  = rand()
    wc  = rand()
    vc  = rand()
    vcp = rand()

    # init buffers for evaluation of rhs
    buff = zeros(Float64, length(a.Γ), length(r.sites))
    temp = zeros(Float64, length(r.sites), length(a.Γ), 2)

    # fill self energy with random values
    a.Σ  .= rand(31)
    ap.Σ .= rand(31)

    # fill channels with random values
    for i in eachindex(a.Γ)
        a.Γ[i].ch_s.q1 .= rand(length(r.sites), 31)
        a.Γ[i].ch_t.q1 .= rand(length(r.sites), 31)
        a.Γ[i].ch_u.q1 .= rand(length(r.sites), 31)

        a.Γ[i].ch_s.q2_1 .= rand(length(r.sites), 31, 31)
        a.Γ[i].ch_t.q2_1 .= rand(length(r.sites), 31, 31)
        a.Γ[i].ch_u.q2_1 .= rand(length(r.sites), 31, 31)

        a.Γ[i].ch_s.q2_2 .= rand(length(r.sites), 31, 31)
        a.Γ[i].ch_t.q2_2 .= rand(length(r.sites), 31, 31)
        a.Γ[i].ch_u.q2_2 .= rand(length(r.sites), 31, 31)

        a.Γ[i].ch_s.q3 .= rand(length(r.sites), 31, 31, 31)
        a.Γ[i].ch_t.q3 .= rand(length(r.sites), 31, 31, 31)
        a.Γ[i].ch_u.q3 .= rand(length(r.sites), 31, 31, 31)

        ap.Γ[i].ch_s.q1 .= rand(length(r.sites), 31)
        ap.Γ[i].ch_t.q1 .= rand(length(r.sites), 31)
        ap.Γ[i].ch_u.q1 .= rand(length(r.sites), 31)

        ap.Γ[i].ch_s.q2_1 .= rand(length(r.sites), 31, 31)
        ap.Γ[i].ch_t.q2_1 .= rand(length(r.sites), 31, 31)
        ap.Γ[i].ch_u.q2_1 .= rand(length(r.sites), 31, 31)

        ap.Γ[i].ch_s.q2_2 .= rand(length(r.sites), 31, 31)
        ap.Γ[i].ch_t.q2_2 .= rand(length(r.sites), 31, 31)
        ap.Γ[i].ch_u.q2_2 .= rand(length(r.sites), 31, 31)

        ap.Γ[i].ch_s.q3 .= rand(length(r.sites), 31, 31, 31)
        ap.Γ[i].ch_t.q3 .= rand(length(r.sites), 31, 31, 31)
        ap.Γ[i].ch_u.q3 .= rand(length(r.sites), 31, 31, 31)
    end

    for rep in 1 : 10
        # time parquet equations
        @timeit to "=> parquet" begin
            # time SDE
            @timeit to "-> SDE" compute_Σ_kernel(Λ, v, wc, r, m, a, ap)

            # time BSE in all channels
            @timeit to "-> BSE s channel" compute_s_BSE!(Λ, buff, v, dv, wc, vc, vcp, r, m, a, temp)
            @timeit to "-> BSE t channel" compute_t_BSE!(Λ, buff, v, dv, wc, vc, vcp, r, m, a, temp)
            @timeit to "-> BSE u channel" compute_u_BSE!(Λ, buff, v, dv, wc, vc, vcp, r, m, a, temp)
        end

        # time flow equations
        @timeit to "=> FRG" begin
            # time self energy
            @timeit to "-> Σ one loop" compute_dΣ_kernel(Λ, v, wc, r, m, a)
            @timeit to "-> Σ corr 1"   compute_dΣ_kernel_corr1(Λ, v, wc, r, m, a, ap)
            @timeit to "-> Σ corr 2"   compute_dΣ_kernel_corr2(Λ, v, wc, r, m, a, ap)

            # time Katanin part in all channels
            @timeit to "-> Katanin s channel" compute_s_kat!(Λ, buff, v, dv, wc, vc, vcp, r, m, a, ap, temp)
            @timeit to "-> Katanin t channel" compute_t_kat!(Λ, buff, v, dv, wc, vc, vcp, r, m, a, ap, temp)
            @timeit to "-> Katanin u channel" compute_u_kat!(Λ, buff, v, dv, wc, vc, vcp, r, m, a, ap, temp)

            # time left part in all channels
            @timeit to "-> left s channel" compute_s_left!(Λ, buff, v, dv, wc, vc, vcp, r, m, a, ap, temp)
            @timeit to "-> left t channel" compute_t_left!(Λ, buff, v, dv, wc, vc, vcp, r, m, a, ap, temp)
            @timeit to "-> left u channel" compute_u_left!(Λ, buff, v, dv, wc, vc, vcp, r, m, a, ap, temp)

            # time central part in all channels
            @timeit to "-> central s channel" compute_s_central!(Λ, buff, v, dv, wc, vc, vcp, r, m, a, ap, temp)
            @timeit to "-> central t channel" compute_t_central!(Λ, buff, v, dv, wc, vc, vcp, r, m, a, ap, temp)
            @timeit to "-> central u channel" compute_u_central!(Λ, buff, v, dv, wc, vc, vcp, r, m, a, ap, temp)
        end
    end

    return nothing
end

"""
    get_flow_timers() :: Nothing

Time current implementation of flow equations by running the different integration kernels.
"""
function get_flow_timers() :: Nothing

    println("NOTE: TIMERS BETWEEN DIFFERENT SYMMETRIES SHALL NOT BE COMPARED (RESULTS ARE MODEL DEPENDENT).")
    println()

    # init test dummys
    list = get_mesh(rand(), 1.0, 30, 0.4)
    m    = Mesh(31, 31, 31, 31, list, list, list, list, list, list)

    # init timer
    to = TimerOutput()

    # time evals of integration kernels for action_su2
    @timeit to "=> action_su2" begin
        # generate action dummy for square lattice Heisenberg model
        l  = get_lattice("square", 5, verbose = false)
        r  = get_reduced_lattice("heisenberg", [[1.0]], l, verbose = false)
        a  = get_action_empty("su2", r, m); init_action!(l, r, a)
        ap = get_action_empty("su2", r, m)

        # time integration kernels
        time_kernels!(to, r, m, a, ap)
    end

    # time evals of integration kernels for action_u1_dm
    @timeit to "=> action_u1_dm" begin
        # generate action dummy for triangular lattice dm-c3 model
        l  = get_lattice("triangular", 5, verbose = false)
        r  = get_reduced_lattice("triangular-dm-c3", [[1.0, 1.0, 1.0]], l, verbose = false)
        a  = get_action_empty("u1-dm", r, m); init_action!(l, r, a)
        ap = get_action_empty("u1-dm", r, m)

        # time integration kernels
        time_kernels!(to, r, m, a, ap)
    end

    show(to)

    return nothing
end


text:
export FacElem, FacElemMon, simplify, factor_coprime

export transform

################################################################################
#
#  Multiplicative representation
#
################################################################################

# abstract nonsense

const FacElemMonDict = ObjectIdDict()

function (x::FacElemMon{S})() where S
  z = FacElem{elem_type(S), S}()
  z.fac = Dict{elem_type(S), fmpz}()
  z.parent = x
  return z
end
doc"""
    FacElem{B}(base::Array{B, 1}, exp::Array{fmpz, 1}) -> FacElem{B}
> Returns the element $\prod b_i^{e_i}$, un-expanded.
"""
function FacElem(base::Array{B, 1}, exp::Array{fmpz, 1}) where B

  length(base) == 0 && error("Array must not be empty")

  length(base) != length(exp) && error("not the same length")

  z = FacElem{B, typeof(parent(base[1]))}()

  for i in 1:length(base)
    if exp[i] == 0
      continue
    end
    if haskey(z.fac, base[i])
      z.fac[base[i]] += exp[i]
    else
      z.fac[base[i]] = exp[i]
    end
  end

  z.parent = FacElemMon(parent(base[1]))
  return z
end

doc"""
    FacElem{B}(d::Dict{B, fmpz}) -> FacElem{B}
    FacElem{B}(d::Dict{B, Integer}) -> FacElem{B}
> Returns the element $\prod b^{d[p]}$, un-expanded.
"""
function FacElem(d::Dict{B, fmpz}) where B

  length(d) == 0 && error("Dictionary must not be empty")

  z = FacElem{B, typeof(parent(first(keys(d))))}()
  z.fac = d

  z.parent = FacElemMon(parent(first(keys(d))))
  return z
end

function FacElem(d::Dict{B, T}) where {B, T <: Integer}

  length(d) == 0 && error("Dictionary must not be empty")

  z = FacElem{B, typeof(parent(first(keys(d))))}()
  z.fac = Dict{B, fmpz}((k,fmpz(v)) for (k,v) = d)

  z.parent = FacElemMon(parent(first(keys(d))))
  return z
end

parent(x::FacElem) = x.parent

base_ring(x::FacElemMon) = x.base_ring

base_ring(x::FacElem) = base_ring(parent(x))

base(x::FacElem) = keys(x.fac)

function Base.deepcopy_internal(x::FacElem{B, S}, dict::ObjectIdDict) where {B, S}
  z = FacElem{B, S}()
  z.fac = Base.deepcopy_internal(x.fac, dict)
  if isdefined(x, :parent)
    z.parent = x.parent
  end
  z.hash = x.hash
  return z
end

function Base.copy(x::FacElem{B, S}) where {B, S}
  z = FacElem{B, S}()
  z.fac = Base.copy(x.fac)
  if isdefined(x, :parent)
    z.parent = x.parent
  end
  z.hash = x.hash
  return z
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, x::FacElemMon)
  print(io, "Factored elements over $(x.base_ring)")
end

function show(io::IO, x::FacElem)
  print(io, "Factored element with data\n$(x.fac)")
end

################################################################################
#
#  Inverse
#
################################################################################
function inv!(x::FacElem)
  for a in base(x)
    x.fac[a] = -x.fac[a]
  end
end

function inv(x::FacElem)
  y = copy(x)
  inv!(y)
  return y
end

################################################################################
#
#  Exponentiation
#
################################################################################

function pow!(z::FacElem, x::FacElem, y::T) where T <: Union{fmpz, Integer}
  z.fac = copy(x.fac)
  for a in base(x)
    # this should be inplace ... not sure anymore: using copy, inplace is bad
    z.fac[a] = y*x.fac[a]
  end
end

# ^(x::FacElem, y::Union{fmpz, Integer}) is ambiguous
for T in [:Integer, fmpz]
  @eval begin
    function ^(x::FacElem, y::($T))
      z = parent(x)()
      if y == 0
        return z
      end
      if y == 1
        return copy(x)
      else
        z.fac = copy(x.fac)
        for a in base(x)
          # this should be inplac
          z.fac[a] = y*x.fac[a]
        end
        return z
      end
    end
  end
end

################################################################################
#
#  Multiplication
#
################################################################################

function mul!(z::FacElem{B, S}, x::FacElem{B, S}, y::FacElem{B, S}) where {B, S}
  z.fac = copy(x.fac)
  for a in base(y)
    if haskey(x.fac, a)
      z.fac[a] = z.fac[a] + y.fac[a]
    else
      z.fac[a] = y.fac[a]
    end
  end
end

function *(x::FacElem{B, S}, y::FacElem{B, S}) where {B, S}
  if length(x.fac) == 0
    return copy(y)
  end

  if length(y.fac) == 0
    return copy(x)
  end

  z = copy(x)
  for a in base(y)
    if haskey(x.fac, a)
      z.fac[a] = z.fac[a] + y.fac[a]
    else
      z.fac[a] = y.fac[a]
    end
  end
  return z
end

function *(x::FacElem{B}, y::B) where B
  z = copy(x)
  if haskey(x.fac, y)
    z.fac[y] = z.fac[y] + 1
  else
    z.fac[y] = 1
  end
  return z
end

function *(y::B, x::FacElem{B}) where B
  z = copy(x)
  if haskey(x.fac, y)
    z.fac[y] = z.fac[y] + 1
  else
    z.fac[y] = 1
  end
  return z
end

function div(x::FacElem{B}, y::FacElem{B}) where B
  z = copy(x)
  for a in base(y)
    if haskey(x.fac, a)
      z.fac[a] = z.fac[a] - y.fac[a]
    else
      z.fac[a] = -y.fac[a]
    end
  end
  return z
end

################################################################################
#
#  Apply transformation
#
################################################################################

# return (x1,...,xr)*y
function _transform(x::Array{FacElem{T, S}, 1}, y::fmpz_mat) where {T, S}
  length(x) != rows(y) &&
              error("Length of array must be number of rows of matrix")

  z = Array{FacElem{T, S}}(cols(y))

  t = parent(x[1])()

  for i in 1:cols(y)
    z[i] = x[1]^y[1,i]
    for j in 2:rows(y)
      if y[j, i] == 0
        continue
      end
      pow!(t, x[j], y[j, i])
      mul!(z[i], z[i], t)
    end
  end
  return z
end

function transform(x::Array{FacElem{S, T}, 1}, y::fmpz_mat) where {S, T}
  return _transform(x, y)
end

function transform_left!(x::Array{FacElem{S, T}, 1}, y::TrafoSwap{fmpz}) where {S, T}
  x[y.i], x[y.j] = x[y.j], x[y.i]
  nothing
end

function transform_left!(x::Array{FacElem{S, T}, 1}, y::TrafoAddScaled{fmpz}) where {S, T}
  x[y.j] = x[y.j] * x[y.i]^y.s
  nothing
end

function transform_left!(x::Array{FacElem{S, T}, 1}, y::TrafoPartialDense{R}) where {S, T, R}
  z = view(deepcopy(x), y.rows)
  xx = view(x, y.rows)
  for i in 1:rows(y.U)
    xx[i] = z[1]^Int(y.U[i, 1])
    for j in 2:cols(y.U)
      xx[i] *= z[j]^Int(y.U[i, j])
    end
  end
end

function transform_left!(x::Array{T, 1}, t::TrafoDeleteZero{S}) where {S, T}
  # move ith position to the back
  for j in t.i:length(x)-1
    r = x[j]
    x[j] = x[j + 1]
    x[j + 1] = r
  end
end

function transform_left!(x::Array{FacElem{S, T}, 1}, y::TrafoParaAddScaled{R}) where {S, T, R}
  # [ Ri; Rj] <- [a b; c d]* [ Ri; Rj ]
  ri = deepcopy(x[y.i])
  rj = deepcopy(x[y.j])
  x[y.i] = ri^y.a * rj^y.b
  x[y.j] = ri^y.c * rj^y.d
end

################################################################################
#
#  Evaluate
#
################################################################################

function evaluate(x::FacElem{NfOrdIdl, NfOrdIdlSet}; coprime::Bool = false)
  O = order(base_ring(x))
  if !coprime
    x = simplify(x) # the other method won't work due to one()
  end
  if length(x.fac)==0
    return frac_ideal(O, O(1))
  end
  # still doesn't work

  return prod([(p//1)^Int(k) for (p,k) = x.fac])
end

function _ev(d::Dict{T, fmpz}, oe::T) where T
  z = copy(oe)
  if length(d)==0
    return z
  elseif length(d)==1
    x = first(d)
    return x[1]^x[2]
  end
  b = similar(d)
  for (k,v) in d
    if v>-10 && v<10
      z *= k^Int(v)
    else
      r = isodd(v) ? 1 :0
      vv = div(v-r, 2)
      if vv!=0
        b[k] = vv
      end
      if r!=0
        z*= k
      end
    end
  end
  return _ev(b, oe)^2*z
end

function one(A::NfOrdFracIdlSet)
  return ideal(order(A), 1)//1
end

function copy(A::NfOrdFracIdl)
  return deepcopy(A)
end

function ^(A::NfOrdFracIdl, d::fmpz)
  return A^Int(d)
end

doc"""
***
  evaluate{T}(x::FacElem{T}) -> T

> Expands or evaluates the factored element, i.e. actually computes the
> value.
> Does "square-and-multiply" on the exponent vectors.
"""
function evaluate(x::FacElem{T}) where T
  return _ev(x.fac, one(base_ring(x)))
end

doc"""
***
    evaluate(x::FacElem{fmpq}) -> fmpq
    evaluate(x::FacElem{fmpz}) -> fmpz

> Expands or evaluates the factored element, i.e. actually computes the
> the element.
> Works by first obtaining a simplified version of the power product
> into coprime base elements.
"""
function evaluate(x::FacElem{fmpq})
  return evaluate_naive(simplify(x))
end

function evaluate(x::FacElem{fmpz})
  return evaluate_naive(simplify(x))
end
doc"""
***
  simplify(x::FacElem{fmpq}) -> FacElem{fmpq}
  simplify(x::FacElem{fmpz}) -> FacElem{fmpz}

> Simplfies the factored element, i.e. arranges for the base to be coprime.
"""
function simplify(x::FacElem{fmpq})
  y = deepcopy(x)
  simplify!(y)
  return y
end

function simplify(x::FacElem{fmpz})
  y = deepcopy(x)
  simplify!(y)
  return y
end

function simplify!(x::FacElem{fmpq})
  if length(x.fac) <= 1
    return
  end
  cp = coprime_base(vcat([denominator(y) for y = base(x)], [numerator(y) for y=base(x)]))
  ev = Dict{fmpq, fmpz}()
  for p = cp
    if p == 1 || p == -1
      continue
    end
    v = fmpz(0)
    for b = base(x)
      v += valuation(b, abs(p))*x.fac[b]
    end
    if v != 0
      ev[fmpq(abs(p))] = v
    end
  end
  s = prod(map(b -> b < 0 && isodd(x.fac[b]) ? -1 : 1, base(x)))
  if s == -1
    ev[fmpq(-1)] = 1
  else
    if length(ev) == 0
      ev[fmpq(1)] = 0
    end
  end
  x.fac = ev
  nothing
end

function simplify!(x::FacElem{fmpz})
  if length(x.fac) <= 1
    return
  end
  cp = coprime_base(collect(base(x)))
  ev = Dict{fmpz, fmpz}()
  for p = cp
    if p == 1 || p == -1
      continue
    end
    v = fmpz(0)
    for b = base(x)
      v += valuation(b, abs(p))*x.fac[b]
    end
    if v < 0
      throw(DomainError())
    end
    if v != 0
      ev[abs(p)] = v
    end
  end
  s = prod(map(b -> b < 0 && isodd(x.fac[b]) ? -1 : 1, base(x)))
  if s == -1
    ev[-1] = 1
  else
    if length(ev) == 0
      ev[fmpz(1)] = 0
    end
  end
  x.fac = ev
  nothing
end

doc"""
***
    isone(x::FacElem{fmpq}) -> Bool
    isone(x::FacElem{fmpz}) -> Bool
> Tests if $x$ represents $1$ without an evaluation.
"""
function isone(x::FacElem{fmpq})
  y = simplify(x)
  return all(iszero, values(y.fac))
end

function isone(x::FacElem{fmpz})
  y = simplify(x)
  return all(iszero, values(y.fac))
end

#TODO: expand the coprime stuff to automatically also get the exponents
function simplify(x::FacElem{NfOrdIdl, NfOrdIdlSet})
  z = deepcopy(x)
  simplify!(z)
  return z
end

function simplify!(x::FacElem{NfOrdIdl, NfOrdIdlSet})
  if length(x.fac) <= 1
    p = first(keys(x.fac))
    x.fac =  Dict(p => x.fac[p])
    return
  end
  cp = coprime_base(collect(base(x)))
  ev = Dict{NfOrdIdl, fmpz}()
  for p = cp
    if isone(p)
      continue
    end
    assure_2_normal(p)
    v = fmpz(0)
    for b = base(x)
      v += valuation(b, p)*x.fac[b]
    end
    if v != 0
      ev[p] = v
    end
  end
  x.fac = ev
end

function simplify(x::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  z = deepcopy(x)
  simplify!(z)
  return z
end

function simplify!(x::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  de = factor_coprime(x)
  if length(de)==0
    de = Dict(ideal(order(base_ring(parent(x))), 1) => fmpz(1))
  end
  x.fac = Dict((i//1, k) for (i,k) = de)
end

function factor_coprime(x::FacElem{NfOrdIdl, NfOrdIdlSet})
  x = deepcopy(x)
  simplify!(x)
  return Dict(p=>Int(v) for (p,v) = x.fac)
end

function factor_coprime(x::FacElem{fmpz})
  x = deepcopy(x)
  simplify!(x)
  d = Dict(abs(p) => Int(v) for (p,v) = x.fac)
  if haskey(d, fmpz(-1))
    delete!(d, fmpz(-1))
    return Fac(fmpz(-1), d)
  else
    return Fac(fmpz(1), d)
  end
end

doc"""
***
  evaluate_naive{T}(x::FacElem{T}) -> T

> Expands or evaluates the factored element, i.e. actually computes the
> value. Uses the obvious naive algorithm. Faster for input in finite rings.
"""
function evaluate_naive(x::FacElem{T}) where T
  z = one(base_ring(x))
  d = x.fac
  if length(d)==0
    return z
  end
  for (k,v) in d
    z *= k^v
  end
  return z
end

function ^(a::fmpz, k::fmpz)
  if a == 0
    if k == 0
      return fmpz(1)
    end
    return fmpz(0)
  end

  if a == 1
    return fmpz(1)
  end
  if a == -1
    if isodd(k)
      return fmpz(-1)
    else
      return fmpz(1)
    end
  end
  return a^Int(k)
end

function ^(a::fmpq, k::fmpz)
  if a == 0
    if k == 0
      return fmpq(1)
    end
    return fmpq(0)
  end

  if a == 1
    return fmpq(1)
  end
  if a == -1
    if isodd(k)
      return fmpq(-1)
    else
      return fmpq(1)
    end
  end
  return a^Int(k)
end

#################################################################################
doc"""
    max_exp(a::FacElem)
> Finds the largest exponent in the factored element $a$
"""
function max_exp(a::FacElem)
  return maximum(values(a.fac))
end

doc"""
    min_exp(a::FacElem)
> Finds the smallest exponent in the factored element $a$
"""
function min_exp(a::FacElem)
  return minimum(values(a.fac))
end

doc"""
    maxabs_exp(a::FacElem)
> Finds the largest exponent by absolute value the factored element $a$
"""
function maxabs_exp(a::FacElem)
  return maximum(abs, values(a.fac))
end

function Base.hash(a::FacElem, u::UInt)
  if a.hash == UInt(0)
    h = hash(u, UInt(3127))
    for (k,v) = a.fac
      h = hash(k, hash(v, h))
    end
    a.hash = h
  end
  return a.hash
end

#used (hopefully) only inside the class group
function FacElem(A::Array{nf_elem_or_fac_elem, 1}, v::Array{fmpz, 1})
  if typeof(A[1]) == nf_elem
    B = FacElem(A[1])
  else
    B = A[1]
  end
  B = B^v[1]
  for i=2:length(A)
    if typeof(A[i]) == nf_elem
      C = FacElem(A[i])^v[i]
    else
      C = A[i]^v[i]
    end
    B *= C
  end
  return B
end

#################################################################################
##
##  Auxillary deep copy functions
##
#################################################################################
#
#function _deepcopy{K, V}(x::Dict{K, V})
#  z = typeof(x)()
#  for a in keys(x)
#    z[a] = deepcopy(x[a])
#  end
#  return z
#end

################################################################################
#
#  Parent object overloading
#
################################################################################

function (F::FacElemMon{T})(a::T) where T
  z = F()
  z.fac[a] = fmpz(1)
  return z
end


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    vecdot(vec1::Matrix{F},vec2::Matrix{F})::Matrix{F} where {F<:AbstractFloat}

or

    vecdot(vec1::Vector{F},vec2::Vector{F})::F where {F<:AbstractFloat}

Determines scalar dot product of two vectors.

# Arguments:
- `vec1`, `vec2`: An array of column vectors represented by a matrix of size
  `(dim,nv)`, where `dim` is the dimension of the vector and `nv` the number
  of vectors.

# Output:
- `dot`: Row vector of scalars, each element corresponding to the dot product
  of the respective components in `vec1` and `vec2`.

# Examples:
```julia
julia> dot = vecdot(vec2,vec2)
```
Determine the dot product of

    [2.3;3.4;5.6] and [1.2;4.5;1.2]
    [5.1;0.0;2.3] and [2.5;3.2;4.0]

```julia
julia> dot = vecdot([2.3 5.1; 3.4 0.0; 5.6 2.3],[1.2 2.5; 4.5 3.2; 1.2 4.0])
```
"""
function vecdot(vec1::Matrix{F},vec2::Matrix{F}
               )::Matrix{F} where {F<:AbstractFloat}
return sum(vec1 .* vec2,dims=1);
end # vecdot

# TODO: look if this is even necessary for this toolbox (using LinearAlgebra)
# overloaded function for vector input to get a float instead of a vector
function vecdot(vec1::Vector{F},vec2::Vector{F})::F where {F<:AbstractFloat}
return vec1' * vec2;
end # vecdot


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

using Markdown
using InteractiveUtils

# ╔═╡ ac31c1f4-0709-11eb-09ab-41e27b54a47b
using Pkg, DrWatson

# ╔═╡ a42e0abc-0709-11eb-27eb-412153f97b8f
begin
	@quickactivate "SR2TuringPluto"
	using Turing
	using StatisticalRethinking
    Turing.setprogress!(false)
end

# ╔═╡ b2616f98-0709-11eb-32a2-514954fe4fa5
md"## Clip-04-42-47t.jl"

# ╔═╡ a131667e-0709-11eb-28be-7906b12ad346
begin
	df = CSV.read(sr_datadir("Howell1.csv"), DataFrame)
	df = df[df.age .>= 18, :]
	x̄ = mean(df.weight)
	x = range(minimum(df.weight), maximum(df.weight), length = 100)     # either
end;

# ╔═╡ caba36be-070a-11eb-3536-3d17bead1bcd
md"### snippet 4.42"

# ╔═╡ d2b9c2ee-070a-11eb-09f8-3390f676f682
@model function m4_3(weights, heights)
    a ~ Normal(178, 20)
    b ~ LogNormal(0, 1)
    σ ~ Uniform(0, 50)
    μ = a .+ b .* (weights .- x̄)
    heights ~ MvNormal(μ, σ)
end

# ╔═╡ d2ba4cfa-070a-11eb-394f-57a4000ca479
m4_3t = m4_3(df.weight, df.height);

# ╔═╡ d2bae2e4-070a-11eb-30fc-638d5366be34
q4_3t = quap(m4_3t, NelderMead())

# ╔═╡ d2cb37f4-070a-11eb-2aee-af17621d1840
md"### snippet 4.43"

# ╔═╡ d2dc3644-070a-11eb-0cd6-19687545a003
@model function m4_3l(weights, heights)
    a ~ Normal(178, 20)
    log_b ~ Normal(0, 1)
    σ ~ Uniform(0, 50)
    μ = a .+ exp(log_b) .* (weights .- mean(weights))
    heights .~ Normal.(μ, σ)
end

# ╔═╡ d2e087b2-070a-11eb-3089-7d6511ec927f
m4_3tl = m4_3l(df.weight, df.height)

# ╔═╡ d2e9d2c4-070a-11eb-2130-6397de3611e4
q4_3tl = quap(m4_3tl, NelderMead())

# ╔═╡ d2f3e980-070a-11eb-0ac4-a9f0c01b3453
md"### snippets 4.44, 4.45"

# ╔═╡ d30746de-070a-11eb-2f08-d76d77bf0853
round.(q4_3t.vcov, digits = 3)

# ╔═╡ d308808c-070a-11eb-080c-27e7c48d4fbe
md"### snippet 4.46"

# ╔═╡ d3222082-070a-11eb-2804-b9341057fa9a
begin
	scatter(df.weight, df.height, leg=false)
	post = DataFrame(rand(q4_3t.distr, 10_000)', q4_3t.params)
	a_map = mean(post.a)
	b_map = mean(post.b)
	plot!(x, a_map .+ b_map .* (x .- x̄))
end

# ╔═╡ d327a962-070a-11eb-10af-0d4f260b4aef
md"### snippets 4.47"

# ╔═╡ d33379f4-070a-11eb-14b7-39e2ba66bd11
Text(precis(post; io=String))

# ╔═╡ c45fa368-0709-11eb-3013-83becc440a56
md"## End of clip-04-42-47t.jl"

# ╔═╡ Cell order:
# ╠═b2616f98-0709-11eb-32a2-514954fe4fa5
# ╠═ac31c1f4-0709-11eb-09ab-41e27b54a47b
# ╠═a42e0abc-0709-11eb-27eb-412153f97b8f
# ╠═a131667e-0709-11eb-28be-7906b12ad346
# ╟─caba36be-070a-11eb-3536-3d17bead1bcd
# ╠═d2b9c2ee-070a-11eb-09f8-3390f676f682
# ╠═d2ba4cfa-070a-11eb-394f-57a4000ca479
# ╠═d2bae2e4-070a-11eb-30fc-638d5366be34
# ╟─d2cb37f4-070a-11eb-2aee-af17621d1840
# ╠═d2dc3644-070a-11eb-0cd6-19687545a003
# ╠═d2e087b2-070a-11eb-3089-7d6511ec927f
# ╠═d2e9d2c4-070a-11eb-2130-6397de3611e4
# ╟─d2f3e980-070a-11eb-0ac4-a9f0c01b3453
# ╠═d30746de-070a-11eb-2f08-d76d77bf0853
# ╟─d308808c-070a-11eb-080c-27e7c48d4fbe
# ╠═d3222082-070a-11eb-2804-b9341057fa9a
# ╟─d327a962-070a-11eb-10af-0d4f260b4aef
# ╠═d33379f4-070a-11eb-14b7-39e2ba66bd11
# ╟─c45fa368-0709-11eb-3013-83becc440a56


text:
@enum(ChannelAssignment,
      Independent,
      LeftSide,
      RightSide,
      MidSide)

@enum(FrameNumberType,
      FrameNumber,
      SampleNumber)

immutable FrameHeader
    blocksize::Cuint
    sample_rate::Cuint
    channels::Cuint
    channel_assignment::ChannelAssignment
    bits_per_sample::Cuint
    typ::FrameNumberType
    num::UInt32
    crc::UInt8
end

const FrameHeaderSync = 0x00003ffe

const MAX_CHANNELS = 0x00000008

immutable FrameFooter
    crc::UInt16
end

@enum(EntropyCodingMethodType,
      PartitionedRice,
      PartitionedRice2)

immutable PartitionedRiceContents
    parameters::Ptr{Cuint}
    raw_bits::Ptr{Cuint}
    capacity_by_order::Cuint
end

immutable PartitionedRiceT
    order::Cuint
    contents::Ptr{PartitionedRiceContents}
end

immutable EntropyCodingMethod
    typ::EntropyCodingMethodType
    data::PartitionedRiceT
end

@enum(SubframeType,
      SubframeConstant,
      SubframeVerbatim,
      SubframeFixed,
      SubframeLPC)

abstract Subframe

immutable Subframe_Constant <: Subframe
    value::Int32
end

immutable Subframe_Verbatim <: Subframe
    data::Ptr{Int32}
end

const MAX_FIXED_ORDER = 0x00000004

immutable SubFrame_Fixed <: Subframe
    method::EntropyCodingMethod
    order::Cuint
    warmup::NTuple{Int(MAX_FIXED_ORDER),Int32}
    residual::Ptr{Int32}
end

const MAX_LPC_ORDER = 0x00000020

immutable SubFrame_LPC <: Subframe
    method::EntropyCodingMethod
    order::Cuint
    qlp_coeff_precision::Cuint
    quantization_level::Cint
    qlp_coeff::NTuple{Int(MAX_LPC_ORDER),Int32}
    warmup::NTuple{Int(MAX_LPC_ORDER),Int32}
    residual::Ptr{Int32}
end

immutable SubFrame{SF <: Subframe}
    typ::SubframeType
    data::SF
    wasted_bits::Cuint
end

immutable Frame
    header::FrameHeader
    subframes::NTuple{Int(MAX_CHANNELS),SubFrame}
    footer::FrameFooter
end


text:
# include("Move.jl")
# include("Wavefunction.jl")
# include("SwapWavefunction.jl")

abstract Walk

is_nontrivial_transition(walk::Walk) = walk.transition_is_nontrivial

function get_probability_ratio(walk::Walk)
    if !walk.transition_is_nontrivial
        warn("Asking for the probability ratio ( = 1) of a trivial transition.")
        return 1.
    else
        return walk.probability_ratio
    end
end

# Lots of code copying below since Julia does't really have abstract base classes.  There has to be a better way..
# I could have one class with finish_move! and cancel_move! functions passed to the constructor and held as members,
# but that seems less safe and less explicit.

# standard walk for ordinary wavefunctions / measurements
type StandardWalk <: Walk
    wf::Wavefunction
    transition_function::Function
    probability_ratio::Float64
    transition_is_nontrivial::Bool
    # transition_in_progress::Bool
    StandardWalk(wf::Wavefunction, transition_function::Function) = new(wf, transition_function)
end

function propose_and_assess_random_transition!(walk::StandardWalk)
    move = walk.transition_function(walk.wf)

    if !is_nontrivial(move)
        walk.transition_is_nontrivial = false
    else
        probrat = assess_move!(walk.wf, move)
        walk.transition_is_nontrivial = true
        walk.probability_ratio = probrat
    end
end

accept_transition!(walk::StandardWalk) = finish_move!(walk.wf)
reject_transition!(walk::StandardWalk) = cancel_move!(walk.wf)

check_for_numerical_error!(walk::StandardWalk) = check_for_numerical_error!(walk.wf)

# walk for Renyi mod measurement
type RenyiModWalk <: Walk
    wf::Union{SwapWavefunction, SwapWavefunctionNS}
    transition_function::Function
    probability_ratio::Float64
    transition_is_nontrivial::Bool
    # transition_in_progress::Bool
    RenyiModWalk(wf::Union{SwapWavefunction, SwapWavefunctionNS}, transition_function::Function) = new(wf, transition_function)
end

function propose_and_assess_random_transition!(walk::RenyiModWalk)
    move = walk.transition_function(walk.wf)
    # move here is either a 2-tuple of moves for the different copies (SwapWavefunction) or a Tuple{ParticlesMove{LatticeParticleMove},Int64} (SwapWavefunctionNS),
    # where the Int here is the copy_to_move
    @assert length(move) == 2

    if isa(move[2], Int) # for SwapWavefunctionNS
        trivial_condition = !is_nontrivial(move[1])
    else # for SwapWavefunction
        trivial_condition = all(map(!, map(is_nontrivial, move)))
    end

    if trivial_condition
        walk.transition_is_nontrivial = false
    else
        probrat = assess_move_for_mod!(walk.wf, move...)
        walk.transition_is_nontrivial = true
        walk.probability_ratio = probrat
    end
end

accept_transition!(walk::RenyiModWalk) = finish_move_for_mod!(walk.wf)
reject_transition!(walk::RenyiModWalk) = cancel_move_for_mod!(walk.wf)

check_for_numerical_error!(walk::RenyiModWalk) = check_for_numerical_error_for_mod!(walk.wf)

# walk for Renyi sign measurement
type RenyiSignWalk <: Walk
    wf::Union{SwapWavefunction, SwapWavefunctionNS}
    transition_function::Function
    probability_ratio::Float64
    transition_is_nontrivial::Bool
    # transition_in_progress::Bool
    RenyiSignWalk(wf::Union{SwapWavefunction, SwapWavefunctionNS}, transition_function::Function) = new(wf, transition_function)
end

function propose_and_assess_random_transition!(walk::RenyiSignWalk)
    move = walk.transition_function(walk.wf)
    @assert length(move) == 2 # a 2-tuple of moves for the different copies

    if isa(move[2], Int) # for SwapWavefunctionNS
        trivial_condition = !is_nontrivial(move[1])
    else # for SwapWavefunction
        trivial_condition = all(map(!, map(is_nontrivial, move)))
    end

    if trivial_condition
        walk.transition_is_nontrivial = false
    else
        probrat = assess_move_for_sign!(walk.wf, move...)
        walk.transition_is_nontrivial = true
        walk.probability_ratio = probrat
    end
end

accept_transition!(walk::RenyiSignWalk) = finish_move_for_sign!(walk.wf)
reject_transition!(walk::RenyiSignWalk) = cancel_move_for_sign!(walk.wf)

check_for_numerical_error!(walk::RenyiSignWalk) = check_for_numerical_error_for_sign!(walk.wf)

# type RenyiWalk <: Walk
#     wf::SwapWavefunction
#     transition_function::Function
#     probability_ratio::Float64
#     transition_is_nontrivial::Bool
#     # transition_in_progress::Bool
#     RenyiWalk(wf::SwapWavefunction, transition_function::Function) = new(wf, transition_function)
# end

# type RenyiModWalk
#     walk::RenyiWalk
#     RenyiModWalk(wf::SwapWavefunction, transition_function::Function) = new(RenyiWalk(wf, transition_function))
# end

# How to then make, e.g., RenyiModWalk.wf invoke RenyiModWalk.walk.wf?


text:
mutable struct wasi_instance_t end

mutable struct WasiConfig
    wasi_config_ptr::Ptr{wasi_config_t}
    program_name::String
    owned::Bool

    WasiConfig(wasi_config_ptr, name) =
        finalizer(new(wasi_config_ptr, name, false)) do config
            config.owned || wasi_config_delete(config)
        end
end
WasiConfig(program_name::String) = WasiConfig(wasi_config_new(), program_name)

Base.unsafe_convert(::Type{Ptr{wasi_config_t}}, wasi_config::WasiConfig) =
    wasi_config.wasi_config_ptr

function set_argv!(config::WasiConfig, argv)
    @assert !config.owned "This WASI configuration is already owned by a store, it can't be modified"
    arr_ptr = Base.unsafe_convert.(Cstring, argv)

    GC.@preserve arr_ptr argv wasi_config_set_argv(config, length(argv), pointer(arr_ptr))
end


text:
export velocity

function velocity(v::Array{Float64,2}, mask, velocity::Array{Float64,1})
    v[:, mask] .*= 0.0
    v[:, mask] .+= vec(velocity)
    return v
end

function velocity(v::Array{Float64,2}, mask, velocity::Array{Float64,2})
    v[:, mask] .= velocity
    return v
end

function velocity(obj::T, f::Function, vel) where T <: Union{Shape, PostOpObj}
    # x, v, y, vol, type
    function func(out)
        v = out[:v]
        mask = vec(f(out))
        return repack!(out, [:v], [v])
    end
    if isa(obj, Shape)
        return PostOpObj(obj, func)
    elseif isa(obj, PostOpObj)
        push!(obj.operations, func)
        return obj
    else
        error("Not allowed.")
    end
end

max_stars_repo_name: sajidmannan/PDMesh

max_issues_repo_name: sajidmannan/PDMesh

max_forks_repo_name: sajidmannan/PDMesh


text:
using MLStyle

mlstyle(node) = @match node begin
    :({ $(kvs...) }) =>
        let f = @λ begin
                :($k : $v) -> Expr(:call, =>, string(k), mlstyle(v))
                :($k) -> Expr(:call, =>, string(k), mlstyle(k))
            end

            Expr(:call, Dict, (f(kv) for kv in kvs if !(kv isa LineNumberNode))...)
        end
    :[$(elts...)] => Expr(:vect, map(mlstyle, elts)...)
    a :: Symbol || a :: T where T <: Number || a :: T where T <: AbstractString => a
end

macro mlstyle(expr)
    mlstyle(expr) |> esc
end


text:
const c = Center()
const f = Face()

const XBoundedGrid = AbstractGrid{<:Any, <:Bounded}
const YBoundedGrid = AbstractGrid{<:Any, <:Any, <:Bounded}
const ZBoundedGrid = AbstractGrid{<:Any, <:Any, <:Any, <:Bounded}

const XYBoundedGrid = AbstractGrid{<:Any, <:Bounded, <:Bounded}
const XZBoundedGrid = AbstractGrid{<:Any, <:Bounded, <:Any, <:Bounded}
const YZBoundedGrid = AbstractGrid{<:Any, <:Any, <:Bounded, <:Bounded}

const FullBoundedGrid =  AbstractGrid{<:Any, <:Bounded, <:Bounded, <:Bounded}

# Fallback for general grids
@inline solid_node(i, j, k, grid) = false
@inline solid_node(i, j, k, grid::XBoundedGrid) = ifelse((i < 1) | (i > grid.Nx), true, false)
@inline solid_node(i, j, k, grid::YBoundedGrid) = ifelse((j < 1) | (j > grid.Ny), true, false)
@inline solid_node(i, j, k, grid::ZBoundedGrid) = ifelse((k < 1) | (k > grid.Nz), true, false)

@inline solid_node(i, j, k, grid::XYBoundedGrid) = ifelse((i < 1) | (i > grid.Nx) | (j < 1) | (j > grid.Ny), true, false)
@inline solid_node(i, j, k, grid::XZBoundedGrid) = ifelse((i < 1) | (i > grid.Nx) | (k < 1) | (k > grid.Nz), true, false)
@inline solid_node(i, j, k, grid::YZBoundedGrid) = ifelse((j < 1) | (j > grid.Ny) | (k < 1) | (k > grid.Nz), true, false)

@inline solid_node(i, j, k, grid::FullBoundedGrid) = ifelse((i < 1) | (i > grid.Nx) |
                                                            (j < 1) | (j > grid.Ny) |
                                                            (k < 1) | (k > grid.Nz), true, false)

@inline solid_node(LX, LY, LZ, i, j, k, grid)      = solid_node(i, j, k, grid)
@inline solid_interface(LX, LY, LZ, i, j, k, grid) = solid_node(i, j, k, grid)

@inline solid_node(::Face, LY, LZ, i, j, k, grid) = solid_node(i, j, k, grid) & solid_node(i-1, j, k, grid)
@inline solid_node(LX, ::Face, LZ, i, j, k, grid) = solid_node(i, j, k, grid) & solid_node(i, j-1, k, grid)
@inline solid_node(LX, LY, ::Face, i, j, k, grid) = solid_node(i, j, k, grid) & solid_node(i, j, k-1, grid)

@inline solid_node(::Face, ::Face, LZ, i, j, k, grid) = solid_node(c, f, c, i, j, k, grid) & solid_node(c, f, c, i-1, j, k, grid)
@inline solid_node(::Face, LY, ::Face, i, j, k, grid) = solid_node(c, c, f, i, j, k, grid) & solid_node(c, c, f, i-1, j, k, grid)
@inline solid_node(LX, ::Face, ::Face, i, j, k, grid) = solid_node(c, f, c, i, j, k, grid) & solid_node(c, f, c, i, j, k-1, grid)

@inline solid_node(::Face, ::Face, ::Face, i, j, k, grid) = solid_node(c, f, f, i, j, k, grid) & solid_node(c, f, f, i-1, j, k, grid)

@inline solid_interface(::Face, LY, LZ, i, j, k, grid) = solid_node(i, j, k, grid) | solid_node(i-1, j, k, grid)
@inline solid_interface(LX, ::Face, LZ, i, j, k, grid) = solid_node(i, j, k, grid) | solid_node(i, j-1, k, grid)
@inline solid_interface(LX, LY, ::Face, i, j, k, grid) = solid_node(i, j, k, grid) | solid_node(i, j, k-1, grid)

@inline solid_interface(::Face, ::Face, LZ, i, j, k, grid) = solid_interface(c, f, c, i, j, k, grid) | solid_interface(c, f, c, i-1, j, k, grid)
@inline solid_interface(::Face, LY, ::Face, i, j, k, grid) = solid_interface(c, c, f, i, j, k, grid) | solid_interface(c, c, f, i-1, j, k, grid)
@inline solid_interface(LX, ::Face, ::Face, i, j, k, grid) = solid_interface(c, f, c, i, j, k, grid) | solid_interface(c, f, c, i, j, k-1, grid)

@inline solid_interface(::Face, ::Face, ::Face, i, j, k, grid) = solid_interface(c, f, f, i, j, k, grid) | solid_interface(c, f, f, i-1, j, k, grid)


text:
# Implement external sort for 4-byte integers
#
# For more information https://en.wikipedia.org/wiki/External_sorting

"""
External sort.
"""
function external_sort(max_memory,
                       input_filename="sample.bin",
                       output_filename="sorted-sample.bin")
    # Phase 1
    #
    # 1. Read max_memory numbers of the data in main memory
    # 2. Sort the data in main memory by some conventional method.
    # 3. Write the sorted data to disk.
    # 4. Repeat 1-3 until end of file.
    input_file = open(input_filename, "r")
    bin_names = Array(ASCIIString,0)
    number_of_bins = 0
    while ~ eof(input_file)
        temporary_vector = Array(Int32, 0)
        # The loop is needed because otherwise it can raise
        # "EOFError: read end of file"
        for i in 1:max_memory
            if eof(input_file)
                break
            else
                push!(temporary_vector, read(input_file, Int32))
            end
        end
        sort!(temporary_vector)

        number_of_bins += 1
        bin_name = string("bin-", number_of_bins, ".bin")
        push!(bin_names, bin_name)
        output_file = open(bin_name, "w")
        write(output_file, temporary_vector)
        close(output_file)
    end
    close(input_file)

    # if number_of_bins is 1 we already have the solution
    if number_of_bins == 1
        mv(bin_names[1], output_filename, remove_destination=true)
        return
    end

    external_sort_phase2(max_memory,
                        [string("bin-", i, ".bin") for i in 1:number_of_bins],
                        output_filename)
end

"""
Second phase of external sort

We reuse the function at the parallel version.
"""
function external_sort_phase2(max_memory,
                              input_filename_list,
                              output_filename="sorted-sample.bin")
    # Phase 2
    #
    # 1. Read the first max_memory / (number_of_bins + 1) of each sorted chunk into
    #    input buffers in main memory and allocate the remaining memory for an
    #    output buffer.
    # 2. Perform a k-way merge and store the result in the output buffer.
    #    Whenever the output buffer fills, write it to the final sorted file and empty it.
    #    Whenever any of the 9 input buffers empties, fill it with the next
    #    sorted chunk until no more data from the chunk is available.
    number_of_bins = length(input_filename_list)

    size_of_sorted_chunk = floor(Integer, max_memory / (number_of_bins + 1))
    if size_of_sorted_chunk == 0
        error("Number of bins is too large. Need more memory.")
    end

    output_file = open(output_filename, "w")

    # Open all bins
    bin_files = Array(IOStream, number_of_bins)
    for i in 1:number_of_bins
        bin_files[i] = open(input_filename_list[i], "r")
    end

    # Read first chunk
    input_buffer = Array(Array{Int32}, number_of_bins)
    for i in 1:number_of_bins
        input_buffer[i] = Array(Int32, 0)
        get_new_chunk(input_buffer[i], bin_files[i], size_of_sorted_chunk)
    end

    while ~ input_buffer_is_empty(input_buffer)
        selected_bin = get_first_non_empty_input_buffer(input_buffer)
        for i in selected_bin + 1:number_of_bins
            if length(input_buffer[i]) > 0 && input_buffer[i][1] < input_buffer[selected_bin][1]
                selected_bin = i
            end
        end
        write(output_file, shift!(input_buffer[selected_bin]))

        if length(input_buffer[selected_bin]) == 0
            get_new_chunk(input_buffer[selected_bin], bin_files[selected_bin], size_of_sorted_chunk)
        end
    end
    close(output_file)

    return
end

"""
Write chunck of data to file.
"""
function write_chunck_to_file(chunck,
                              output_filename,
                              sort=false,
                              append=true)
    if sort
        sort!(chunck)
    end

    if append
        mode = "a"
    else
        mode = "w"
    end
    output_file = open(output_filename, mode)
    write(output_file, chunck)
    close(output_file)
end

"""
Get new chunck of sorted data.
"""
function get_new_chunk(input_buffer_array,
                       input_file_stream,
                       size_of_sorted_chunk)
    for j in 1:size_of_sorted_chunk
        if eof(input_file_stream)
            close(input_file_stream)
        else
            push!(input_buffer_array, read(input_file_stream, Int32))
        end
    end
end

"""
Check if input buffer is empty.
"""
function input_buffer_is_empty(input_buffer)
    for chunck in input_buffer
        if length(chunck) != 0
            return false
        end
    end
    return true
end

"""
Get first non empty buffer.
"""
function get_first_non_empty_input_buffer(input_buffer)
    for i in 1:length(input_buffer)
        if length(input_buffer[i]) != 0
            return i
        end
    end
    return 0
end

max_stars_repo_name: rgaiacs/external-sort

max_issues_repo_name: rgaiacs/external-sort

max_forks_repo_name: rgaiacs/external-sort


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "nPth"
version = v"1.6"

# Collection of sources required to build libgcrypt
sources = [
    ArchiveSource("https://gnupg.org/ftp/gcrypt/npth/npth-1.6.tar.bz2",
                  "1393abd9adcf0762d34798dc34fdcf4d0d22a8410721e76f1e3afcd1daa4e2d1"),
]

# Bash recipe for building across all platforms

# Tried -no-undefined but still couldn't build for windows
script = raw"""
cd $WORKSPACE/srcdir/npth-*/
if [[ "${target}" == powerpc64le-* ]]; then
    autoreconf -vi
fi
update_configure_scripts

./configure --prefix=${prefix} --host=${target} --build=${MACHTYPE} --disable-static
make -j${nproc}
make install

install_license ${WORKSPACE}/srcdir/npth-*/COPYING.LIB
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.  We are manually disabling
# many platforms that do not seem to work.
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["libnpth", "libnpth6"], :libnpth),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using LatticeQCD

betas = ["1.40","1.65"]
println("# beta plaquette")
for beta in betas
        dir = "samples/HMC_L08080808_su2nf8_beta$(beta)_Staggered_mass0.015"
        p = get_plaquette_average(dir)
        println("$beta $p #$dir")
end
println("# ")
println("# beta Polyakov loop")
for beta in betas
        dir = "samples/HMC_L08080808_su2nf8_beta$(beta)_Staggered_mass0.015"
        p = get_polyakov_average(dir)
        println("$beta $p #$dir")
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "MPIR"
version = v"3.0.0"

# Collection of sources required to build MPFRBuilder
sources = [
    ArchiveSource("http://mpir.org/mpir-$(version).tar.bz2",
                  "52f63459cf3f9478859de29e00357f004050ead70b45913f2c2269d9708675bb"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/mpir-*

# We need `yasm`
apk add yasm

./configure --enable-cxx --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-static --enable-shared
make -j
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [p for p in supported_platforms() if arch(p) == "x86_64" && !Sys.isfreebsd(p)]

# The products that we will ensure are always built
products = [
    LibraryProduct("libmpir", :libmpir)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)