text:
#####
##### Nothing happens when your boundary condition is nothing
#####

fill_west_and_east_halo!(c,  ::Nothing, ::Nothing, args...; kwargs...) = NoneEvent()
fill_south_and_north_halo!(c,::Nothing, ::Nothing, args...; kwargs...) = NoneEvent()
fill_bottom_and_top_halo!(c, ::Nothing, ::Nothing, args...; kwargs...) = NoneEvent()

for dir in (:west, :east, :south, :north, :bottom, :top)
        fill_nothing! = Symbol( :fill_, dir, :_halo!)
    alt_fill_nothing! = Symbol(:_fill_, dir, :_halo!)
    @eval begin
        @inline     $fill_nothing!(c, ::Nothing, args...;  kwargs...)         = NoneEvent()
        @inline $alt_fill_nothing!(i, j, grid, c, ::Nothing, args...)         = nothing
        @inline $alt_fill_nothing!(i, j, grid, ::Nothing, ::Nothing, args...) = nothing
        @inline $alt_fill_nothing!(i, j, grid, ::Nothing, args...)            = nothing
    end
end


text:
function check_for_test_param(cur_item, cur_test_param)

  nested_files = get_nested_files(cur_item)
  clean_files_list(nested_files)

  cur_regex = Regex("$(cur_test_param)\\s*:\\s*true")

  for nested_file in nested_files

    opened_file = open(nested_file)

    for cur_line in eachline(opened_file)
      stripped_line = strip(cur_line)

      if ( stripped_line == "" )
        continue
      end

      if ( !startswith(stripped_line, "#") )
        break
      end

      if ismatch(cur_regex, stripped_line)
        return true
      end
    end

    close(opened_file)

  end

  return false

end

ext: jl


text:
using GRASS
using Documenter

DocMeta.setdocmeta!(GRASS, :DocTestSetup, :(using GRASS); recursive=true)

Introduction = "Introduction" => "index.md"

Getting_Started = "Getting Started" => ["installation.md", "examples/basic.md", "examples/advanced.md"]

License = "License" => "license.md"

Index = "Index" => "longlist.md"

pages = [Introduction, Getting_Started, License, Index]

makedocs(;
    modules=[GRASS],
    authors="Michael Palumbo",
    repo="https://github.com/palumbom/GRASS/blob/{commit}{path}#{line}",
    sitename="GRASS",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://palumbom.github.io/GRASS",
        assets=String[],
    ),
    pages=pages,
)

deploydocs(;
    repo="github.com/palumbom/GRASS",
    devbranch = "main",
)

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
include("building.jl")

# ================
# BLDF01 - BDS01
# ================

# general options
𝑂_BLDF01 = Options(:T => time_horizon, :mode => "check", :property => pBDS01)

# algorithm-specific options
𝑂_dense_BLDF01 = Options(:vars => [25], :partition => [1:24, [25], 26:48],
                         :δ => 0.003, :block_options_init => LazySets.LinearMap)
𝑂_discrete_BLDF01 =
    merge(𝑂_dense_BLDF01, Options(:discretization => "nobloating", :δ => 0.01))

# single run to verify that specification holds
sol_BLDF01_dense = solve(build_TV, 𝑂_BLDF01, op=BFFPSV18(𝑂_dense_BLDF01))
@assert sol_BLDF01_dense.satisfied
sol_BLDF01_discrete = solve(build_TV, 𝑂_BLDF01, op=BFFPSV18(𝑂_discrete_BLDF01))
@assert sol_BLDF01_discrete.satisfied

# benchmark
SUITE["Build"]["BLDF01-BDS01", "dense"] =
    @benchmarkable solve($build_TV, $𝑂_BLDF01, op=BFFPSV18($𝑂_dense_BLDF01))
SUITE["Build"]["BLDF01-BDS01", "discrete"] =
    @benchmarkable solve($build_TV, $𝑂_BLDF01, op=BFFPSV18($𝑂_discrete_BLDF01))

# ================
# BLDC01 - BDS01
# ================

# general options
𝑂_BLDC01 = Options(:T => time_horizon, :mode => "check", :property => pBLDC01)

# algorithm-specific options
𝑂_dense_BLDC01 = Options(:vars => [25], :partition => [1:24, [25], 26:48, [49]],
                         :δ => 0.005, :block_options_init => LazySets.LinearMap)
𝑂_discrete_BLDC01 =
    merge(𝑂_dense_BLDC01, Options(:discretization => "nobloating", :δ => 0.005))

# single run to verify that specification holds
sol_BLDC01_dense = solve(build_CONST, 𝑂_BLDC01, op=BFFPSV18(𝑂_dense_BLDC01))
@assert sol_BLDC01_dense.satisfied
sol_BLDC01_discrete = solve(build_CONST, 𝑂_BLDC01, op=BFFPSV18(𝑂_discrete_BLDC01))
@assert sol_BLDC01_discrete.satisfied

# register in benchmark suite
SUITE["Build"]["BLDC01-BDS01", "dense"] =
    @benchmarkable solve($build_CONST, $𝑂_BLDC01, op=BFFPSV18($𝑂_dense_BLDC01))
SUITE["Build"]["BLDC01-BDS01", "discrete"] =
    @benchmarkable solve($build_CONST, $𝑂_BLDC01, op=BFFPSV18($𝑂_discrete_BLDC01))

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
using Test
using Pkg.Artifacts
using NCDatasets

using RRTMGP
using RRTMGP.Device: array_type, array_device
using RRTMGP.Vmrs
using RRTMGP.LookUpTables
using RRTMGP.AtmosphericStates
using RRTMGP.Optics
using RRTMGP.Sources
using RRTMGP.BCs
using RRTMGP.Fluxes
using RRTMGP.AngularDiscretizations
using RRTMGP.RTE
using RRTMGP.RTESolver

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()
# overriding some parameters to match with RRTMGP FORTRAN code
CLIMAParameters.Planet.grav(::EarthParameterSet) = 9.80665
CLIMAParameters.Planet.molmass_dryair(::EarthParameterSet) = 0.028964
CLIMAParameters.Planet.molmass_water(::EarthParameterSet) = 0.018016

include("reference_files.jl")
include("read_rfmip_clear_sky.jl")
#---------------------------------------------------------------
function clear_sky(
    ::Type{OPC},
    ::Type{SRC},
    ::Type{VMR},
    ::Type{FT},
    ::Type{I},
    ::Type{DA},
) where {FT<:AbstractFloat,I<:Int,DA,OPC,SRC,VMR}
    opc = Symbol(OPC)
    lw_file = get_ref_filename(:lookup_tables, :clearsky, λ = :lw) # lw lookup tables for gas optics
    sw_file = get_ref_filename(:lookup_tables, :clearsky, λ = :sw) # sw lookup tables for gas optics
    input_file = get_ref_filename(:atmos_state, :clearsky)         # clear-sky atmos state
    # reference data files for comparison
    flux_up_file_lw = get_ref_filename(
        :comparison,
        :clearsky,
        λ = :lw,
        flux_up_dn = :flux_up,
        opc = opc,
    )
    flux_dn_file_lw = get_ref_filename(
        :comparison,
        :clearsky,
        λ = :lw,
        flux_up_dn = :flux_dn,
        opc = opc,
    )
    flux_up_file_sw = get_ref_filename(
        :comparison,
        :clearsky,
        λ = :sw,
        flux_up_dn = :flux_up,
        opc = opc,
    )
    flux_dn_file_sw = get_ref_filename(
        :comparison,
        :clearsky,
        λ = :sw,
        flux_up_dn = :flux_dn,
        opc = opc,
    )

    FTA1D = DA{FT,1}
    FTA2D = DA{FT,2}
    max_threads = Int(256)
    exp_no = 1
    n_gauss_angles = 1

    # reading longwave lookup data
    ds_lw = Dataset(lw_file, "r")
    lookup_lw, idx_gases = LookUpLW(ds_lw, I, FT, DA)
    close(ds_lw)

    # reading shortwave lookup data
    ds_sw = Dataset(sw_file, "r")
    lookup_sw, idx_gases = LookUpSW(ds_sw, I, FT, DA)
    close(ds_sw)

    # reading rfmip data to atmospheric state
    ds_lw_in = Dataset(input_file, "r")

    (as, sfc_emis, sfc_alb_direct, zenith, toa_flux, usecol) = setup_rfmip_as(
        ds_lw_in,
        idx_gases,
        exp_no,
        lookup_lw,
        FT,
        DA,
        VMR,
        max_threads,
    )
    close(ds_lw_in)

    ncol, nlay, ngpt_lw = as.ncol, as.nlay, lookup_lw.n_gpt
    nlev = nlay + 1
    op = OPC(FT, ncol, nlay, DA) # allocating optical properties object

    # setting up longwave problem
    ngpt_lw = lookup_lw.n_gpt
    src_lw = source_func_longwave(FT, ncol, nlay, opc, DA)         # allocating longwave source function object
    bcs_lw = LwBCs{FT,typeof(sfc_emis),Nothing}(sfc_emis, nothing) # setting up boundary conditions
    ang_disc = nothing#AngularDiscretization(opc, FT, n_gauss_angles, DA)  # initializing Angular discretization
    fluxb_lw = FluxLW(ncol, nlay, FT, DA)                          # flux storage for bandwise calculations
    flux_lw = FluxLW(ncol, nlay, FT, DA)                           # longwave fluxes
    # setting up shortwave problem
    ngpt_sw = lookup_sw.n_gpt
    src_sw = source_func_shortwave(FT, ncol, nlay, opc, DA)        # allocating shortwave source function object
    inc_flux_diffuse = nothing
    sfc_alb_diffuse = FTA2D(deepcopy(sfc_alb_direct))
    bcs_sw = SwBCs(
        zenith,
        toa_flux,
        sfc_alb_direct,
        inc_flux_diffuse,
        sfc_alb_diffuse,
    )

    fluxb_sw = FluxSW(ncol, nlay, FT, DA) # flux storage for bandwise calculations
    flux_sw = FluxSW(ncol, nlay, FT, DA)  # shortwave fluxes for band calculations
    #--------------------------------------------------
    # initializing RTE solver
    slv = Solver(
        as,
        op,
        src_lw,
        src_sw,
        bcs_lw,
        bcs_sw,
        fluxb_lw,
        fluxb_sw,
        flux_lw,
        flux_sw,
    )

    solve_lw!(slv, max_threads, lookup_lw)
    solve_sw!(slv, max_threads, lookup_sw)

    #for i = 1:10
    #    @time solve_lw!(slv, max_threads, lookup_lw)
    #    @time solve_sw!(slv, max_threads, lookup_sw)
    #end

    # comparing longwave fluxes with data from RRTMGP FORTRAN code
    flip_ind = nlev:-1:1

    ds_flux_up_lw = Dataset(flux_up_file_lw, "r")
    comp_flux_up_lw = ds_flux_up_lw["rlu"][:][flip_ind, :, exp_no]
    close(ds_flux_up_lw)

    ds_flux_dn_lw = Dataset(flux_dn_file_lw, "r")
    comp_flux_dn_lw = ds_flux_dn_lw["rld"][:][flip_ind, :, exp_no]
    close(ds_flux_dn_lw)

    flux_up_lw = Array(slv.flux_lw.flux_up)
    flux_dn_lw = Array(slv.flux_lw.flux_dn)

    max_err_flux_up_lw = maximum(abs.(flux_up_lw .- comp_flux_up_lw))
    max_err_flux_dn_lw = maximum(abs.(flux_dn_lw .- comp_flux_dn_lw))

    println("=======================================")
    println("Clear-sky longwave test - $opc")
    println("max_err_flux_up_lw = $max_err_flux_up_lw")
    println("max_err_flux_dn_lw = $max_err_flux_dn_lw")

    toler_lw = 1e-4
    #--------------------------------------------------------------
    # comparing shortwave fluxes with data from RRTMGP FORTRAN code
    ds_flux_up_sw = Dataset(flux_up_file_sw, "r")
    comp_flux_up_sw = ds_flux_up_sw["rsu"][:][flip_ind, :, exp_no]
    close(ds_flux_up_sw)

    ds_flux_dn_sw = Dataset(flux_dn_file_sw, "r")
    comp_flux_dn_sw = ds_flux_dn_sw["rsd"][:][flip_ind, :, exp_no]
    close(ds_flux_dn_sw)

    flux_up_sw = Array(slv.flux_sw.flux_up)
    flux_dn_sw = Array(slv.flux_sw.flux_dn)

    for i = 1:ncol
        if usecol[i] == 0
            flux_up_sw[:, i] .= FT(0)
            flux_dn_sw[:, i] .= FT(0)
        end
    end

    max_err_flux_up_sw = maximum(abs.(flux_up_sw .- comp_flux_up_sw))
    max_err_flux_dn_sw = maximum(abs.(flux_dn_sw .- comp_flux_dn_sw))

    println("Clear-sky shortwave test, opc  = $opc")
    println("max_err_flux_up_sw = $max_err_flux_up_sw")
    println("max_err_flux_dn_sw = $max_err_flux_dn_sw")
    println("=======================================")

    toler_sw = FT(0.001)

    @test max_err_flux_up_lw ≤ toler_lw
    @test max_err_flux_dn_lw ≤ toler_lw
    @test max_err_flux_up_sw ≤ toler_sw
    @test max_err_flux_dn_sw ≤ toler_sw
end

clear_sky(TwoStream, SourceLW2Str, VmrGM, Float64, Int, array_type())


text:
#####
##### constructor
#####

function rrule(::Type{T}, x::Number) where {T<:UniformScaling}
    UniformScaling_back(dx) = (NoTangent(), ProjectTo(x)(unthunk(dx).λ))
    return T(x), UniformScaling_back
end

#####
##### `+`
#####

function frule((_, Δx, ΔJ), ::typeof(+), x::AbstractMatrix, J::UniformScaling)
    return x + J, Δx + (zero(J) + ΔJ)  # This (0 + ΔJ) allows for ΔJ::Tangent{UniformScaling}
end

function frule((_, ΔJ, Δx), ::typeof(+), J::UniformScaling, x::AbstractMatrix)
    return J + x, (zero(J) + ΔJ) + Δx
end

function rrule(::typeof(+), x::AbstractMatrix, J::UniformScaling)
    project_x = ProjectTo(x)
    project_J = ProjectTo(J)
    function plus_back(dy)
        dx = unthunk(dy)
        return (NoTangent(), project_x(dx), project_J(I * tr(dx)))
    end
    return x + J, plus_back
end

function rrule(::typeof(+), J::UniformScaling, x::AbstractMatrix)
    y, back = rrule(+, x, J)
    function plus_back_2(dy)
        df, dx, dJ = back(dy)
        return (df, dJ, dx)
    end
    return y, plus_back_2
end

#####
##### `-`
#####

function frule((_, Δx, ΔJ), ::typeof(-), x::AbstractMatrix, J::UniformScaling)
    return x - J, Δx - (zero(J) + ΔJ)
end

function frule((_, ΔJ, Δx), ::typeof(-), J::UniformScaling, x::AbstractMatrix)
    return J - x, (zero(J) + ΔJ) - Δx
end

function rrule(::typeof(-), x::AbstractMatrix, J::UniformScaling)
    y, back = rrule(+, x, -J)
    project_J = ProjectTo(J)
    function minus_back_1(dy)
        df, dx, dJ = back(dy)
        return (df, dx, project_J(-dJ))  # re-project as -true isa Int
    end
    return y, minus_back_1
end

function rrule(::typeof(-), J::UniformScaling, x::AbstractMatrix)
    project_x = ProjectTo(x)
    project_J = ProjectTo(J)
    function minus_back_2(dy)
        dx = -unthunk(dy)
        return (NoTangent(), project_J(-tr(dx) * I), project_x(dx))
    end
    return J - x, minus_back_2
end

#####
##### `Matrix`
#####

function rrule(::Type{T}, I::UniformScaling{<:Bool}, (m, n)) where {T<:AbstractMatrix}
    Matrix_back_I(dy) = (NoTangent(), NoTangent(), NoTangent())
    return T(I, m, n), Matrix_back_I
end

function rrule(::Type{T}, J::UniformScaling, (m, n)) where {T<:AbstractMatrix}
    project_J = ProjectTo(J)
    function Matrix_back_I(dy)
        dJ = if m == n
            project_J(I * tr(unthunk(dy)))
        else
            project_J(I * sum(diag(unthunk(dy))))
        end
        return (NoTangent(), dJ, NoTangent())
    end
    return T(J, m, n), Matrix_back_I
end


text:
@testset "1579.remove-max-number-of-edges-to-keep-graph-fully-traversable.jl" begin
    @test max_num_edges_to_remove(
        4, [[3, 1, 2], [3, 2, 3], [1, 1, 3], [1, 2, 4], [1, 1, 2], [2, 3, 4]]
    ) == 2
    @test max_num_edges_to_remove(4, [[3, 2, 3], [1, 1, 2], [2, 3, 4]]) == -1
    @test max_num_edges_to_remove(4, [[3, 1, 2], [3, 2, 3], [1, 1, 4], [2, 1, 4]]) == 0
    @test max_num_edges_to_remove(4, [[3, 2, 3], [1, 1, 2], [1, 3, 4]]) == -1
end


text:
## test/frametools.jl -- whether some frame tools work as expected

### TODO: for these tests you need a toy frame to work on....
## I guess, it is fine to reload the data from extractions

@testset "Chain Dict" begin
    Trajectory(joinpath(dataloc, fname)) do traj
        ## stuff from Chemfiles
        fr = read(traj)
        top = Topology(fr)

        ## extracted from frame
        extracted = extractframe(fr)

        # now for the good part
        test_dict = getchaindict(extracted)

        ## check that chains are mutually exclusive (no index in diff chains)
        @test let rlist = extracted.res_list
            Q = true
            for i in eachindex(rlist)
                ## nonstandard (PDB) residues should not be in the chaindict
                res = @inbounds rlist[i]
                tst = !(res.standard_pdb)
                checker = @inbounds ifelse(tst, 0, 1)
                count = 0
                for (_, range_) in test_dict
                    if i in range_ count += 1 end
                end
                Q &= (count == checker)
                Q || break
            end
            Q
        end

        ## check that the chain ids of those in the extracted frame
        ## are the same for results in test_dict
        @test let rlist = extracted.res_list
            Q = true
            for i in eachindex(rlist)
                (; chainid, standard_pdb) = @inbounds rlist[i]
                ## remember to skip nonstandard residues
                standard_pdb || continue
                chid_ = Symbol(chainid)
                Q &= (i in test_dict[chid_])
                Q || break
            end
            Q
        end
    end
end

@testset "Backbone amide geometry, etc." begin
    Trajectory(joinpath(dataloc, fname)) do traj
        ## stuff from Chemfiles
        fr = read(traj)
        top = Topology(fr)

        ## extracted from frame; add protons to it
        extracted = extractframe(fr)
        has_nuc = fname == "3sn2.pdb"
        addprotons!(extracted; has_na=has_nuc)
        kAAN = ProtStructs.kAminoAcidNames

        innertest(res) = begin
            nuc = has_nuc && res.name ∉ kAAN
            ## should short circuit when possible; so it makes sense
            ## that the above holds when the residue is a nucleic acid
            !res.standard_pdb || res.name == :PRO || nuc
        end

        ## check the number of atoms; make sure length of positions is same
        @testset "Frame details" begin
            (; at_list, at_pos) = extracted
            ## lengths test
            @test length(at_list) == length(at_pos)
            ## if indexing works as expected
            @test let rlist = extracted.res_list
                Q = true
                ## drop one because proton is not added there
                for res in Iterators.drop(rlist, 1)
                    innertest(res) && continue
                    Q &= (at_list[getatom(res, :H)].name == :H)
                    Q || break
                end
                Q
            end
            ## number of protons added is same as standard protein res - 1
            @test let rlist = extracted.res_list
                ct = 0
                for res in Iterators.drop(rlist, 1)
                    innertest(res) && continue
                    for (name, _) in res.at_dict
                        if name == :H ct += 1 end
                    end
                end
                ct == count(Iterators.drop(rlist, 1)) do res
                    !innertest(res)
                end
            end
        end

        ## check that the bond length is 1A (as in DSSP HBond method)
        @test let (; at_pos, res_list) = extracted
            Q = true
            ## drop one because proton is not added there
            for res in Iterators.drop(res_list, 1)
                ## skip those that arent standard protein residues
                innertest(res) && continue
                Npos = at_pos[getatom(res, :N)]
                Hpos = at_pos[getatom(res, :H)]
                Q &= isapprox(ProtStructs.distance(Npos, Hpos), 1)
                Q || break
            end
            Q
        end

        ## check that the carbonyl bond is parallel to the amine bond
        ## (OCN) bond angle = (CNH) bond angle
        @test let (; at_pos, res_list) = extracted
            Q = true
            for i in Iterators.drop(eachindex(res_list), 1)
                res = @inbounds res_list[i]
                res_prev = @inbounds res_list[i-1]
                innertest(res) && continue
                res_prev.standard_pdb || continue
                Opos = at_pos[getatom(res_prev, :O)]
                Cpos = at_pos[getatom(res_prev, :C)]
                Npos = at_pos[getatom(res, :N)]
                Hpos = at_pos[getatom(res, :H)]

                OCN = anglespan(Opos, Cpos, Npos)
                CNH = anglespan(Cpos, Npos, Hpos)
                Q &= isapprox(OCN, CNH)
                Q || break
            end
            Q
        end

        ## check that the amide is planar (-π/π dihedral angle)
        @test let (; at_pos, res_list) = extracted
            Q = true
            for i in Iterators.drop(eachindex(res_list), 1)
                res = @inbounds res_list[i]
                res_prev = @inbounds res_list[i-1]
                innertest(res) && continue
                res_prev.standard_pdb || continue
                Opos = at_pos[getatom(res_prev, :O)]
                Cpos = at_pos[getatom(res_prev, :C)]
                Npos = at_pos[getatom(res, :N)]
                Hpos = at_pos[getatom(res, :H)]

                Q &= isapprox(abs(ProtStructs.dihedral(Opos, Cpos, Npos, Hpos)), π)
                Q || break
            end
            Q
        end
    end
end


text:
OFFSETS = (N = (0, 1), S = (0, -1), E = (1, 0), W = (-1, 0))

function solve(case, start, M)
    curr = start

    if sum(abs, curr) == 0
        println("Case #$(case): 0")
        return
    end

    for (i, x) in enumerate(M)
        curr .+= OFFSETS[Symbol(x)]
        if sum(abs, curr) <= i
            println("Case #$(case): $i")
            return
        end
    end

    if sum(abs, curr) < length(M) + 1
        println("Case #$(case): $(length(M))")
        return
    end

    println("Case #$(case): IMPOSSIBLE")
end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        _X, _Y, M = stdin |> readline |> split
        start = [parse(Int, _X), parse(Int, _Y)]
        solve(case, start, M)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
using BinaryBuilder, Pkg

name = "CUDA"
version = v"11.0.3"

dependencies = [BuildDependency(PackageSpec(name="CUDA_full_jll", version=version))]

script = raw"""
# First, find (true) CUDA toolkit directory in ~/.artifacts somewhere
CUDA_ARTIFACT_DIR=$(dirname $(dirname $(realpath $prefix/cuda/bin/ptxas${exeext})))
cd ${CUDA_ARTIFACT_DIR}

# Clear out our prefix
rm -rf ${prefix}

# license
install_license EULA.txt

# headers
mkdir -p ${prefix}/include
mv include/* ${prefix}/include
rm -rf ${prefix}/include/thrust

# binaries
mkdir -p ${bindir} ${libdir} ${prefix}/lib ${prefix}/share
if [[ ${target} == *-linux-gnu ]]; then
    # CUDA Runtime
    mv lib64/libcudart.so* lib64/libcudadevrt.a ${libdir}

    # CUDA FFT Library
    mv lib64/libcufft.so* lib64/libcufftw.so* ${libdir}

    # CUDA BLAS Library
    mv lib64/libcublas.so* lib64/libcublasLt.so* ${libdir}

    # NVIDIA "Drop-in" BLAS Library
    mv lib64/libnvblas.so* ${libdir}

    # CUDA Sparse Matrix Library
    mv lib64/libcusparse.so* ${libdir}

    # CUDA Linear Solver Library
    mv lib64/libcusolver.so* ${libdir}

    # CUDA Linear Solver Multi GPU Library
    mv lib64/libcusolverMg.so* ${libdir}

    # CUDA Random Number Generation Library
    mv lib64/libcurand.so* ${libdir}

    # NVIDIA Performance Primitives Library
    mv lib64/libnpp*.so* ${libdir}

    # NVIDIA Optimizing Compiler Library
    mv nvvm/lib64/libnvvm.so* ${libdir}

    # NVIDIA Common Device Math Functions Library
    mkdir ${prefix}/share/libdevice
    mv nvvm/libdevice/libdevice.10.bc ${prefix}/share/libdevice

    # CUDA Profiling Tools Interface (CUPTI) Library
    mv extras/CUPTI/lib64/libcupti.so* ${libdir}

    # NVIDIA Tools Extension Library
    mv lib64/libnvToolsExt.so* ${libdir}

    # Additional binaries
    mv bin/nvdisasm ${bindir}
    mv bin/cuda-memcheck ${bindir}
elif [[ ${target} == x86_64-w64-mingw32 ]]; then
    # CUDA Runtime
    mv bin/cudart64_*.dll ${bindir}
    mv lib/x64/cudadevrt.lib ${prefix}/lib

    # CUDA FFT Library
    mv bin/cufft64_*.dll bin/cufftw64_*.dll ${bindir}

    # CUDA BLAS Library
    mv bin/cublas64_*.dll bin/cublasLt64_*.dll ${bindir}

    # NVIDIA "Drop-in" BLAS Library
    mv bin/nvblas64_*.dll ${bindir}

    # CUDA Sparse Matrix Library
    mv bin/cusparse64_*.dll ${bindir}

    # CUDA Linear Solver Library
    mv bin/cusolver64_*.dll ${bindir}

    # CUDA Linear Solver Multi GPU Library
    mv bin/cusolverMg64_*.dll ${bindir}

    # CUDA Random Number Generation Library
    mv bin/curand64_*.dll ${bindir}

    # NVIDIA Performance Primitives Library
    mv bin/npp*64_*.dll ${bindir}

    # NVIDIA Optimizing Compiler Library
    mv nvvm/bin/nvvm64_*.dll ${bindir}

    # NVIDIA Common Device Math Functions Library
    mkdir ${prefix}/share/libdevice
    mv nvvm/libdevice/libdevice.10.bc ${prefix}/share/libdevice

    # CUDA Profiling Tools Interface (CUPTI) Library
    mv extras/CUPTI/lib64/cupti64_*.dll ${bindir}

    # NVIDIA Tools Extension Library
    mv bin/nvToolsExt64_1.dll ${bindir}

    # Additional binaries
    mv bin/nvdisasm.exe ${bindir}
    mv bin/cuda-memcheck.exe ${bindir}

    # Fix permissions
    chmod +x ${bindir}/*.{exe,dll}
fi
"""

products = [
    LibraryProduct(["libcudart", "cudart64_110"], :libcudart),
    FileProduct(["lib/libcudadevrt.a", "lib/cudadevrt.lib"], :libcudadevrt),
    LibraryProduct(["libcufft", "cufft64_10"], :libcufft),
    LibraryProduct(["libcufftw", "cufftw64_10"], :libcufftw),
    LibraryProduct(["libcublas", "cublas64_11"], :libcublas),
    LibraryProduct(["libcublasLt", "cublasLt64_11"], :libcublasLt),
    LibraryProduct(["libnvblas", "nvblas64_11"], :libnvblas),
    LibraryProduct(["libcusparse", "cusparse64_11"], :libcusparse),
    LibraryProduct(["libcusolver", "cusolver64_10"], :libcusolver),
    LibraryProduct(["libcusolverMg", "cusolverMg64_10"], :libcusolverMg),
    LibraryProduct(["libcurand", "curand64_10"], :libcurand),
    LibraryProduct(["libnppc", "nppc64_11"], :libnppc),
    LibraryProduct(["libnppial", "nppial64_11"], :libnppial),
    LibraryProduct(["libnppicc", "nppicc64_11"], :libnppicc),
    LibraryProduct(["libnppidei", "nppidei64_11"], :libnppidei),
    LibraryProduct(["libnppif", "nppif64_11"], :libnppif),
    LibraryProduct(["libnppig", "nppig64_11"], :libnppig),
    LibraryProduct(["libnppim", "nppim64_11"], :libnppim),
    LibraryProduct(["libnppist", "nppist64_11"], :libnppist),
    LibraryProduct(["libnppisu", "nppisu64_11"], :libnppisu),
    LibraryProduct(["libnppitc", "nppitc64_11"], :libnppitc),
    LibraryProduct(["libnpps", "npps64_11"], :libnpps),
    LibraryProduct(["libnvvm", "nvvm64_33_0"], :libnvvm),
    FileProduct("share/libdevice/libdevice.10.bc", :libdevice),
    LibraryProduct(["libcupti", "cupti64_2020.1.1"], :libcupti),
    LibraryProduct(["libnvToolsExt", "nvToolsExt64_1"], :libnvtoolsext),
    ExecutableProduct("nvdisasm", :nvdisasm),
]

build_tarballs(ARGS, name, version, [], script,
               [Platform("x86_64", "linux"), Platform("powerpc64le", "linux"), Platform("x86_64", "windows")],
               products, dependencies)

max_stars_repo_name: Affie/Yggdrasil

max_issues_repo_name: Affie/Yggdrasil

max_forks_repo_name: Affie/Yggdrasil


text:
# dnds.jl
# =======
#
# dNdS statistics computation.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/NaturalSelection.jl/blob/master/LICENSE.md

const CDN_POS_MASKS = (0xFFFFFFFFFFFFFFCF, 0xFFFFFFFFFFFFFFF3, 0xFFFFFFFFFFFFFFFC)
const SITE_PERMUTATIONS = [[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]]

function aligned_codons(x::BioSequence{T}, y::BioSequence{T}, start::Int = 1) where T <: NucAlphs
    xcdns = Vector{Kmer{eltype(T), 3}}()
    ycdns = Vector{Kmer{eltype(T), 3}}()
    pos = start
    while pos + 2 ≤ min(endof(x), endof(y))
        cdnx, okx = BioSequences.extract_kmer_impl(x, pos, 3)
        cdny, oky = BioSequences.extract_kmer_impl(y, pos, 3)
        if okx && oky
            push!(xcdns, convert(Kmer{eltype(T), 3}, cdnx))
            push!(ycdns, convert(Kmer{eltype(T), 3}, cdny))
        end
        pos += 3
    end
    return xcdns, ycdns
end

function aligned_codons(v::Vector{DNASequence}, start::Int = 1)
    nv = endof(v)
    cdns = Vector{Vector{Codon{DNA}}}(nv)
    oks = Vector{Vector{Bool}}(nv)
    @inbounds for i in eachindex(v)
        cdns[i] = Vector{Codon{DNA}}()
        oks[i] = Vector{Bool}()

        pos = start
        while pos + 2 ≤ minimum(endof.(v))
            cdn, ok = BioSequences.extract_kmer_impl(v[i], pos, 3)
            push!(cdns[i], convert(Codon{DNA}, cdn))
            push!(oks[i], ok)
            pos += 3
        end
    end
    finalbool = falses(oks[1])
    @inbounds for ok in oks
        for i in eachindex(ok)
            finalbool[i] |= ok[i]
        end
    end
    return Vector{Codon{DNA}}[cdn[finalbool] for cdn in cdns]
end

include("enumerator/enumerator.jl")
include("enumerator/allpaths.jl")
include("enumerator/shortestpath.jl")

include("lookups/lookups.jl")

include("NG86/NG86.jl")

include("mkt/mkt.jl")


text:
"""
    ArbNumerics

A performant package for numerics with settable extended precision.

primary exported types: `ArbFloat`, `ArbReal`, `ArbComplex`
"""
module ArbNumerics

export ArbNumber,
       ArbFloat, ArbReal, ArbComplex,
       interval, setinterval,
       lowerbound, upperbound,
       ball, setball,
       midpoint, radius,
       interval_abs,
       lowerbound_abs, upperbound_abs,
       increase_radius, decrease_radius,
       bits4digits, digits4bits,
       bits2digits, digits2bits,
       workingprecision, setworkingprecision,
       extrabits, setextrabits,

       stringall, showall,

       ispositive, isnegative, isnonpositive, isnonnegative,

       isspecial,     # zero, infinite or NaN
       isnonfinite,
       isinf, isposinf, isneginf,
       isnonzero,
       isnormal,      # nonzero finite
       isexact,       # zero radius
       isinexact,     # nonzero radius
       isinteger,     # integer valued
       isnoninteger,  # has nonzero fractional portion

       inf, posinf, neginf, nan,
       typemax, typemin, floatmax, floatmin,
       magnitude,     # complex magnitude, `angle` gives phase

       signs, signbits,
       significand_bits, relerror_bits, relaccuracy_bits,
       trim_bits, ulp, ulps, eps, epss,

       # other functions
       fmod, # inv of modf
       round!,
       add2!, sub2!, mul2!, div2!,
       addmul, submul, mulsub,
       square, cube, rsqrt, rcbrt, pow, root, loghypot,
       risingfactorial, doublefactorial,
       tanpi, cotpi,

       # special functions
       agm1, agm,
       erf, erfc, erfi, erfcx,
       logsinpi, barnesg, logbarnesg,
       ei, si, ci, shi, chi,
       airyai, airyaiprime, airybi, airybiprime,
       besselj, besselj0, besselj1, bessely, bessely0, bessely1, besseli, besselk,
       hypergeometric_0F1, hypergeometric_1F1, hypergeometric_2F1,
       regular_hypergeometric_0F1, regular_hypergeometric_1F1, regular_hypergeometric_2F1,
       F₀₁,  F₁₁, F₂₁, regularF₀₁,  regularF₁₁, regularF₂₁,
       elliptic_k, elliptic_e, elliptic_pi, elliptic_f,
       elliptic_k2, elliptic_e2, elliptic_pi2, elliptic_f2, # modulus^2
       elliptic_rf, elliptic_rg, elliptic_rj,
       weierstrass_p, weierstrass_invp, weierstrass_zeta, weierstrass_sigma,
       zeta, eta, xi,                  # Reimann
       lambertw, polylog,
       π, ℯ, γ, φ, catalan,

       # fft
       dft, inverse_dft,

       # special linear algebra
       determinant, inverse, fastmul, ⨱ # \timesbar

import Base: IEEEFloat,
       hash, promote_rule, promote_type, convert,
       copy, deepcopy, string, show,

       zero, one,

       Float16, Float32, Float64, float,
       UInt8, UInt16, UInt32, UInt64, UInt128,
       Int8, Int16, Int32, Int64, Int128,
       BigInt, BigFloat, Rational, Complex, real, imag, complex,
       floatmax, floatmin, typemax, typemin, maxintfloat,
       rationalize,

       iszero, isone, isinteger, isfinite, isinf, isnan, issubnormal,
       zero, one,
       precision, setprecision, eps, nextfloat, prevfloat, widen, isapprox,

       abs, abs2, signbit, sign, copysign, flipsign,

       (==), (!=), (<), (<=), (>=), (>), isequal, isless,
       (+), (-), (*), (/), (\), inv,
       div, fld, cld, rem, mod, divrem, fldmod,
       ceil, floor, trunc, round,
       muladd, fma,
       factorial, binomial,
       conj

import Base.Math:  max, min, minmax,
       exponent, significand,
       ldexp, frexp, modf,

       mod2pi, rem2pi,
       sqrt, cbrt, hypot,

       (^), exp, log, expm1, log1p, log2, log10,
       sin, cos, sincos, tan, csc, sec, cot, sinc, sinpi, cospi,
       asin, acos, atan, acsc, asec, acot,
       sinh, cosh, tanh, csch, sech, coth,
       asinh, acosh, atanh, acsch, asech, acoth,

       clamp, clamp!

import Random: rand, randn

import SpecialFunctions: gamma, lgamma, lfact, digamma, invdigamma, polygamma, trigamma,
     beta, lbeta, sinint, cosint,
     erf, erfc, erfcinv, erfcx, erfi, erfinv,
     airyai, airyaiprime, airybi, airybiprime,
     besselj, besselj0, besselj1, bessely, bessely0, bessely1, besseli, besselk,
     eta, zeta

using  GenericSVD

using  LinearAlgebra
import LinearAlgebra: tr, det, transpose, transpose!, norm, lu, ldlt,
                      cholesky, tril, triu, eigvals, svdvals, floatmin2,
                      mul!, rmul!, lmul!, eigvecs, svd, eigen, dot

export tr, det, transpose, transpose!, norm, mul!, lmul!, rmul!, lu, ldlt,
       cholesky, tril, tiru, eigvals, svdvals

import Base.MathConstants: π, ℯ, γ, φ, catalan

using Libdl
using Random
using Random: SamplerType, SamplerTrivial, CloseOpen01

using FLINT_jll
using Arb_jll

include("support/arblib.jl")
include("support/abstractions.jl")
include("support/matrices.jl")

include("libarb/roundingmodes.jl")
include("libarb/ArbMatrix.jl")
include("libarb/Flint2.jl")
include("libarb/Mag.jl")
include("libarb/ArbFloat.jl")
include("libarb/ArbReal.jl")
include("libarb/ArbComplex.jl")
include("support/complex.jl")

include("support/ArblibVector.jl")

include("support/random.jl")

const ArbNumber = Union{ArbFloat, ArbReal, ArbComplex}

include("libarb/ArbMatrix.jl")  # must preceed ArbRealMatrix
include("libarb/ArbRealMatrix.jl")  # must preceed ArbFloatMatrix
include("libarb/ArbFloatMatrix.jl") # must preceed ArbComplexMatrix
include("libarb/ArbComplexMatrix.jl")
include("libarb/ArbFFT.jl")
include("libarb/promote.jl")
include("libarb/string.jl")
include("libarb/parse.jl")
include("libarb/show.jl")

include("values/bitsdigits.jl")
include("values/intraconvert.jl")
include("values/interconvert.jl")
include("values/constructors.jl")
include("values/specialvals.jl")
include("values/compare.jl")
include("values/rounding.jl")
include("values/random.jl")

include("float/parametrics.jl")
include("float/prearith.jl")
include("float/arith.jl")
include("float/arith_inplace.jl")
include("float/morearith.jl")
include("float/muladd.jl")
include("float/dot.jl")
include("float/elementary.jl")
include("float/otherspecial.jl")
include("float/bessel.jl")
include("float/airy.jl")
include("float/elliptic.jl")
include("float/hypergeometric.jl")
include("float/export_fallbacks.jl")

include("float/float64.jl")

include("intervals/intervalparts.jl")
include("intervals/complex_intervalparts.jl")
include("intervals/eps_ulp.jl")

include("support/minprec.jl")
include("support/linearalgebra.jl")
include("support/printf.jl")
include("support/helptext.jl")

function __init__()
    ccall(@libflint(flint_set_num_threads), Cvoid, (Cint,), Sys.CPU_THREADS)
    return nothing
end

end # ArbNumerics


text:
import JuliaDB.ML
using JuliaDB
using Base.Test
using OnlineStats
using PooledArrays
using DataValues

@testset "feature extraction" begin
    @testset "schema" begin
        @test ML.schema([1:10;]).series == fit!(Series(Variance()), [1:10;])
        x = repeat([1,2], inner=5)
        @test ML.schema(x, ML.Categorical).series == fit!(Series(CountMap(Int)), x) == ML.schema(PooledArray(x)).series
        m = ML.schema(DataValueArray(x))
        @test m isa ML.Maybe
        @test m.feature.series == ML.schema(x).series
        t = table(PooledArray([1,2,1,2]), [1,2,3,4], names=[:x, :y])
        sch = ML.schema(t)
        @test sch[:x].series == ML.schema(PooledArray([1,2,1,2])).series
        @test sch[:y].series == ML.schema([1,2,3,4]).series
    end

    @testset "featuremat" begin
        @test ML.featuremat([1,3,5]) == [-1.0, 0, 1.0]'
        Δ = 0.5 / std([1,2])
        @test ML.featuremat(DataValueArray([1,2,3], Bool[0,0,1])) ≈ [0 0 1; -Δ Δ 0]
        @test isempty(ML.featuremat(["x","y","x"]))
        @test ML.featuremat(PooledArray(["x","y","x"])) == [1  0  1; 0  1  0]

        t = table(PooledArray([1,1,2,2]), [1,2,3,4], DataValueArray([1,2,3,4]), ["x", "y", "z", "a"], names=[:a,:b,:c,:d])
        @test ML.featuremat(t) == collect(ML.featuremat(distribute(t, 2)))

        x = randn(100)
        @test ML.featuremat(x) ≈ ((x .- mean(x)) ./ std(x))'
    end

    @testset "Issue 149" begin
        x = rand(10), rand(1:4, 10), rand(1:4, 10)
        t  = table(x..., names = [:x1, :x2, :x3])
        td = table(x..., names = [:x1, :x2, :x3], chunks = 2)
        sch = ML.schema(t, hints=Dict(
            :x1 => ML.Continuous,
            :x2 => ML.Categorical,
            :x3 => ML.Categorical,
            )
        )
        schd = ML.schema(td, hints=Dict(
            :x1 => ML.Continuous,
            :x2 => ML.Categorical,
            :x3 => ML.Categorical,
            )
        )

        @test schd[:x1].series.stats[1].σ2 ≈ sch[:x1].series.stats[1].σ2
        @test schd[:x2].series == sch[:x2].series
        @test schd[:x3].series == sch[:x3].series
    end
end


text:
using Mimi

@defcomp grosseconomy begin
    YGROSS  = Variable(index=[time])    # Gross output
    K       = Variable(index=[time])    # Capital
    l       = Parameter(index=[time])   # Labor
    tfp     = Parameter(index=[time])   # Total factor productivity
    s       = Parameter(index=[time])   # Savings rate
    depk    = Parameter()               # Depreciation rate on capital - Note that it has no time index
    k0      = Parameter()               # Initial level of capital
    share   = Parameter()               # Capital share

    function run_timestep(p, v, d, t)
        # Define an equation for K
        if is_first(t)
            v.K[t]  = p.k0  # Note the use of v. and p. to distinguish between variables and parameters
        else
            v.K[t]  = (1 - p.depk)^5 * v.K[t-1] + v.YGROSS[t-1] * p.s[t-1] * 5
        end

        # Define an equation for YGROSS
        v.YGROSS[t] = p.tfp[t] * v.K[t]^p.share * p.l[t]^(1-p.share)
    end
end

@defcomp emissions begin
    E       = Variable(index=[time])    # Total greenhouse gas emissions
    sigma   = Parameter(index=[time])   # Emissions output ratio
    YGROSS  = Parameter(index=[time])   # Gross output - Note that YGROSS is now a parameter

    function run_timestep(p, v, d, t)
        # Define an eqation for E
        v.E[t] = p.YGROSS[t] * p.sigma[t]   # Note the p. in front of YGROSS
    end
end

function construct_model()
	m = Model()

	set_dimension!(m, :time, collect(2015:5:2110))

	# Order matters here. If the emissions component were defined first, the model would not run.
	add_comp!(m, grosseconomy)
	add_comp!(m, emissions)

	# Update parameters for the grosseconomy component
	update_param!(m, :grosseconomy, :l, [(1. + 0.015)^t *6404 for t in 1:20])
	update_param!(m, :grosseconomy, :tfp, [(1 + 0.065)^t * 3.57 for t in 1:20])
	update_param!(m, :grosseconomy, :s, ones(20).* 0.22)
	update_param!(m, :grosseconomy, :depk, 0.1)
	update_param!(m, :grosseconomy, :k0, 130.)
	update_param!(m, :grosseconomy, :share, 0.3)

	# Update and connect parameters for the emissions component
	update_param!(m, :emissions, :sigma, [(1. - 0.05)^t *0.58 for t in 1:20])
	connect_param!(m, :emissions, :YGROSS, :grosseconomy, :YGROSS)

	return m

end #end function

model = construct_model()

run(model)

# Show model results
model[:emissions, :E]

# Or, as a DataFrame
getdataframe(model, :emissions, :E)


text:
function updaterh!(system::CVSystem,n::Int64)
    # pulmonary, tricuspid valve states
    if system.heart.ra.P[n+1] > system.heart.rv.P[n+1]
        TV = "open";
    else
        TV = "closed";
    end

    if (system.heart.rv.P[n+1] > system.lungs.Pp[n+1] ||
        system.heart.rv.Q[n+1] > 0)
        PV = "open";
    else
        PV = "closed";
    end

    # right atrial V, Q (P updated elsewhere)
    system.heart.ra.V[n+2] = (system.heart.ra.V[n+1]+
        system.solverparams.h*(system.ivc.Q[n+1] + system.svc.Q[n+1]-
        system.heart.ra.Q[n+1]));
    if TV == "open"
        system.heart.ra.Q[n+2] = ((1 - system.heart.ra.R/system.heart.ra.L*
            system.solverparams.h)*system.heart.ra.Q[n+1]+
            system.solverparams.h/system.heart.ra.L*(system.heart.ra.P[n+1]-
            system.heart.rv.P[n+1]));
    else
        system.heart.ra.Q[n+2] = 0;
    end

    # right ventricular V, Q
    system.heart.rv.V[n+2] = (system.heart.rv.V[n+1]+
        system.solverparams.h*(system.heart.ra.Q[n+1]-
        system.heart.rv.Q[n+1]));
    if PV == "open"
        system.heart.rv.Q[n+2] = (system.heart.rv.Q[n+1]+
            system.solverparams.h/system.heart.rv.L*(system.heart.rv.P[n+1]-
            system.lungs.Pp[n+1]));
    else
        system.heart.rv.Q[n+2] = 0;
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
using Symbolics, Test
using ReferenceTests

@variables t a x(t) y(t)
expr = [a*x - x*y,-3y + x*y]
@test_reference "target_functions/1.stan" Symbolics.build_function(expr,[x,y],[a],t,target = Symbolics.StanTarget())

@test_reference "target_functions/1.c" Symbolics.build_function(expr,[x,y],[a],t,target = Symbolics.CTarget(),
                                     lhsname=:internal_var___du,
                                     rhsnames=[:internal_var___u,:internal_var___p,:t])

@test_reference "target_functions/1.m" Symbolics.build_function(expr,[x,y],[a],t,target = Symbolics.MATLABTarget())

@test Symbolics.build_function(expr,[x,y],[a],t,target = Symbolics.CTarget(),
                                     lhsname=:internal_var___du,
                                     rhsnames=[:internal_var___u,:internal_var___p,:t]) ==
    Symbolics.build_function(expr,[x,y],[a],t,target = Symbolics.CTarget(),
                                   lhsname=:internal_var___du,
                                   rhsnames=[:internal_var___u,:internal_var___p,:t])

@test Symbolics.build_function(expr,[x,y],[a],t,target = Symbolics.StanTarget()) ==
      Symbolics.build_function(expr,[x,y],[a],t,target = Symbolics.StanTarget())

@test Symbolics.build_function(expr,[x,y],[a],t,target = Symbolics.MATLABTarget()) ==
      Symbolics.build_function(expr,[x,y],[a],t,target = Symbolics.MATLABTarget())

@test Symbolics.build_function(expr,[x,y],[a],t,target = Symbolics.CTarget(),
                                     lhsname=:internal_var___du,
                                     rhsnames=[:internal_var___u,:internal_var___p,:t]) ==
      Symbolics.build_function(expr,vcat(x,y),[a],t,target = Symbolics.CTarget(),
                                     lhsname=:internal_var___du,
                                     rhsnames=[:internal_var___u,:internal_var___p,:t])

@test Symbolics.build_function(expr,[x,y],[a],t,target = Symbolics.StanTarget()) ==
      Symbolics.build_function(expr,vcat(x,y),[a],t,target = Symbolics.StanTarget())

@test Symbolics.build_function(expr,[x,y],[a],t,target = Symbolics.MATLABTarget()) ==
      Symbolics.build_function(expr,vcat(x,y),[a],t,target = Symbolics.MATLABTarget())

# Matrix CTarget test
let
    @variables x[1:4] y[1:4] z[1:4]
    expression = hcat(x,y,z) # Results in [x1 y1 z1; x2 y2 z2; ...]
    variables  = vcat(x,y,z) # Results in [x1, x2, x3, x4, y1, ...]
    cfunc      = build_function(expression, variables; target = Symbolics.CTarget(), expression = Val{true})
    @test_reference "target_functions/matrix.c" cfunc
end

# Scalar CTarget test
let
    @variables x y t
    expression = x + y + t
    cfunc = build_function(expression, [x], [y], t; target = Symbolics.CTarget(), expression = Val{true})
    @test_reference "target_functions/scalar1.c" cfunc
end

# Scalar CTarget test with scalar multiplication and powers
let
    @variables x y a t
    expression = x^2 + y^-1 + sin(a)^3.5 + 2t
    cfunc = build_function(expression, [x, y], [a], t; target = Symbolics.CTarget(), expression = Val{true})

    @test_reference "target_functions/scalar2.c" cfunc
end

# Matrix StanTarget test
let
    @variables x[1:4] y[1:4] z[1:4]
    expression = hcat(x,y,z) # Results in [x1 y1 z1; x2 y2 z2; ...]
    variables  = vcat(x,y,z) # Results in [x1, x2, x3, x4, y1, ...]
    sfunc      = build_function(expression, vcat(x,y), z, []; target = Symbolics.StanTarget(), expression = Val{true})

    @test_reference "target_functions/matrix.stan" sfunc
end

# Scalar StanTarget test
let
    @variables t x(t) y(t) z(t)
    expression = x + y + z
    sfunc = build_function(expression, vcat(x,y), [z], t; target = Symbolics.StanTarget(), expression = Val{true})

    @test_reference "target_functions/scalar1.stan" sfunc
end

# Matrix MATLABTarget test
let
    @variables x[1:4] y[1:4] z[1:4]
    expression = hcat(x,y,z) # Results in [x1 y1 z1; x2 y2 z2; ...]
    variables  = vcat(x,y,z) # Results in [x1, x2, x3, x4, y1, ...]
    mfunc      = build_function(expression, vcat(x,y,z); target = Symbolics.MATLABTarget(), expression = Val{true})

    @test_reference "target_functions/matrix.m" mfunc
end

# Scalar MATLABTarget test
let
    @variables x y z
    expression = x + y + z
    mfunc = build_function(expression, vcat(x,y,z); target = Symbolics.MATLABTarget(), expression = Val{true})

    @test_reference "target_functions/scalar1.m" mfunc
end


text:
module AppliGeneralLedger

include("./infrastructure/infrastructure.jl")

export create_journal_entry, process, read_from_file

end # module


text:
using Plots, JLD, Statistics, LaTeXStrings
number_list=[i for i in 0:9]
NList=[100,1000,10000,100000,1000000]
STDList=[]
final_res=[]
for N in NList
    res=rand(number_list, N)
    push!(final_res, res)
    frequency_list=[]
    for i in number_list
        push!(frequency_list, length(findall(x->x==i, res)))
    end
    push!(STDList,std(frequency_list)/N)
end

p1=histogram(final_res[1], bins=15, alpha=0.5, xlabel="The Random Number", ylabel="Number",c = :red, label="N=100")
p2=histogram(final_res[2], bins=15, alpha=0.5, xlabel="The Random Number", ylabel="Number",c = :green,label="N=1000")
p3=histogram(final_res[3], bins=15, alpha=0.5, xlabel="The Random Number", ylabel="Number",c = :yellow,label="N=10000")
p4=histogram(final_res[4], bins=15, alpha=0.5, xlabel="The Random Number", ylabel="Number",c = :blue,label="N=100000")
p5=histogram(final_res[5], bins=15, alpha=0.5, xlabel="The Random Number", ylabel="Number",c = :orange,label="N=1000000")
plot(p1,p2,p3,p4,p5, layout=5,dpi=600,guidefontsize=9, legendfont=6, legend = :outertop)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter6\\6.1\\Figs\\TheFrequencyCruve.png")


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

using Markdown
using InteractiveUtils

# ╔═╡ 822b4836-fe75-11ea-161f-35d589ab833f
using Pkg, DrWatson

# ╔═╡ 822b886e-fe75-11ea-3fd5-313bd9474e67
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanQuap
	using StatisticalRethinking
end

# ╔═╡ 416556d6-fe6e-11ea-338a-fb39b7c0f00f
md"## Clip-06-08s.jl"

# ╔═╡ 822c0d84-fe75-11ea-27b2-cb1e8fa46fb1
begin
	df = CSV.read(sr_datadir("milk.csv"), DataFrame; delim=';');
	scale!(df, [:kcal_per_g, :perc_fat, :perc_lactose])
end;

# ╔═╡ 823966f0-fe75-11ea-3c37-25437bd882cd
stan6_3 = "
data{
  int <lower=1> N;              // Sample size
  vector[N] K;
  vector[N] F;
}
parameters{
  real a;
  real bF;
  real<lower=0> sigma;
}
model{
  vector[N] mu;
  sigma ~ exponential( 1 );
  a ~ normal( 0 , 0.2 );
  bF ~ normal( 0 , 0.5 );
  mu = a + bF * F;
  K ~ normal( mu , sigma );
}
";

# ╔═╡ 82435d7e-fe75-11ea-3f52-7daa1879898d
begin
	data = (N = size(df, 1), F = df.perc_fat_s, K = df.kcal_per_g_s)
	init = (a=0.0, bF=1.0, sigma=1.0)
	q6_3s, m6_3s, o6_3s = stan_quap("m6.3s", stan6_3; data, init);
	if !isnothing(m6_3s)
		post6_3s_df = read_samples(m6_3s, :dataframe)
		PRECIS(post6_3s_df)
	end
end

# ╔═╡ 82456bc6-fe75-11ea-3e9a-256e1fb9844f
!isnothing(m6_3s) && (part6_3s = Particles(post6_3s_df))

# ╔═╡ 82509064-fe75-11ea-16cb-e506f49efd72
if !isnothing(q6_3s)
	quap6_3s_df = sample(q6_3s)
	PRECIS(quap6_3s_df)
end

# ╔═╡ 8251194e-fe75-11ea-1403-d743608a5169
hpdi(part6_3s.bF.particles, alpha=0.11)

# ╔═╡ 825cb330-fe75-11ea-2df5-a92dca779683
md"## End of clip-06-08s.jl"

# ╔═╡ Cell order:
# ╟─416556d6-fe6e-11ea-338a-fb39b7c0f00f
# ╠═822b4836-fe75-11ea-161f-35d589ab833f
# ╠═822b886e-fe75-11ea-3fd5-313bd9474e67
# ╠═822c0d84-fe75-11ea-27b2-cb1e8fa46fb1
# ╠═823966f0-fe75-11ea-3c37-25437bd882cd
# ╠═82435d7e-fe75-11ea-3f52-7daa1879898d
# ╠═82456bc6-fe75-11ea-3e9a-256e1fb9844f
# ╠═82509064-fe75-11ea-16cb-e506f49efd72
# ╠═8251194e-fe75-11ea-1403-d743608a5169
# ╟─825cb330-fe75-11ea-2df5-a92dca779683


text:
# using Visualize: vert_linesegments, frag_linesegments, geom_linesegments, LineSegments, Vertices, LineAttributes
#
#
# function Drawable(w::AbstractGLWindow, primitive::LineSegments)
#     vbo = VertexArray(primitive[Vertices], face_type = Face{2, OffsetInteger{1, GLint}})
#     uniforms = UniformBuffer(LineAttributes(primitive))
#     args = (w[Scene], uniforms)
#     rasterizer(
#         vbo, args,
#         vert_linesegments, frag_linesegments;
#         geometryshader = geom_linesegments,
#         primitive_in = :lines
#     )
# end


text:
# This file plots the `findvolume.png` in docs
using Unitful, UnitfulAtomic, Plots
using EquationsOfState, EquationsOfState.Collections, EquationsOfState.Find

plotlyjs()

pressures = collect(0:20:200) .* u"GPa"
eos = BirchMurnaghan3rd(167 * u"angstrom^3", 2600 * u"kbar", 4.0)
volumes = map(
    p -> findvolume(PressureForm(), eos, p, (eps() * u"bohr^3", eos.v0 * 1.3)),
    pressures,
)
plot(ustrip.(volumes), ustrip.(pressures), label = "pressures")
scatter!(
    ustrip.(volumes),
    ustrip.(u"GPa", apply(PressureForm(), eos).(volumes)),
    label = "P(V)",
    dpi = 400,
)
xlabel!("volume (angstrom³)")
ylabel!("pressure (GPa)")
savefig("docs/src/assets/findvolume.png")


text:
@syms a b c

using Metatheory

@testset "Equality" begin
    @eqtest a == a
    @eqtest a != b
    @eqtest a*b == a*b
    @eqtest a*b != a
    @eqtest a != a*b
end

@testset "Literal Matcher" begin
    r = @rule 1 => 4
    @test r(1) === 4
    @test r(1.0) === 4
    @test r(2) === nothing
end

@testset "Slot matcher" begin
    @test @rule(~x => true)("?") === true
    @test @rule( ~x => ~x)(2) === 2
end

@testset "Term matcher" begin
    @test @rule(sin(~x) => ~x)(sin(a)) === a
    @eqtest @rule(sin(~x) => ~x)(sin(a^2)) == a^2
    @test @rule(sin(~x) => ~x)(sin(a)^2) === nothing
    @test @rule(sin(sin(~x)) => ~x)(sin(a^2)) === nothing
    @test @rule(sin(sin(~x)) => ~x)(sin(sin(a))) === a
    @test @rule(sin(~x)^2 => ~x)(sin(a)^2) === a
end

@testset "Equality matching" begin
    @test @rule((~x)^(~x) => ~x)(a^a) === a
    @test @rule((~x)^(~x) => ~x)(b^a) === nothing
    @test @rule((~x)^(~x) => ~x)(a+a) === nothing
    @eqtest @rule((~x)^(~x) => ~x)(sin(a)^sin(a)) == sin(a)
    @eqtest @rule((~x*~y + ~x*~z)  => ~x * (~y+~z))(a*b + a*c) == a*(b+c)

    @eqtest @rule(+(~~x) => ~~x)(a + b) == [a,b]
    @eqtest @rule(+(~~x) => ~~x)(term(+, a, b, c)) == [a,b,c]
    @eqtest @rule(+(~~x,~y, ~~x) => (~~x, ~y))(term(+,9,8,9,type=Any)) == ([9,],8)
    @eqtest @rule(+(~~x,~y, ~~x) => (~~x, ~y, ~~x))(term(+,9,8,9,9,8,type=Any)) == ([9,8], 9, [9,8])
    @eqtest @rule(+(~~x,~y,~~x) => (~~x, ~y, ~~x))(term(+,6,type=Any)) == ([], 6, [])
end

using SymbolicUtils: @capture

@testset "Capture form" begin

    ex = a^a

    #note that @test inserts a soft local scope (try-catch) that would gobble
    #the matches from assignment statements in @capture macro, so we call it
    #outside the test macro
    ret = @capture ex (~x)^(~x)
    @test ret
    @test @isdefined x
    @test x === a

    ex = b^a
    ret = @capture ex (~y)^(~y)
    @test !ret
    @test !(@isdefined y)

    ret = @capture (a + b) (+)(~~z)
    @test ret
    @test @isdefined z
    @test all(z .=== arguments(a + b))

    #a more typical way to use the @capture macro

    f(x) = if @capture x (~w)^(~w)
        w
    end

    @eqtest f(b^b) == b
    @test f(b+b) == nothing
end


text:
function compute_potential_vorticity!(q,H1h,H1hchol,dΩ,R,S,h,u,f,n)
  a(r,s) = ∫(s*h*r)dΩ
  c(s)   = ∫(perp(n,∇(s))⋅(u) + s*f)dΩ
  Gridap.FESpaces.assemble_matrix_and_vector!(a, c, H1h, get_free_dof_values(q), R, S)
  lu!(H1hchol, H1h)
  ldiv!(H1hchol, get_free_dof_values(q))
end

function compute_velocity!(u1,dΩ,dω,V,RTMMchol,u2,qAPVM,F,ϕ,n,dt1,dt2)
  b(v) = ∫(v⋅u2 - dt1*(qAPVM)*(v⋅⟂(F,n)))dΩ + ∫(dt2*DIV(v)*ϕ)dω
  Gridap.FESpaces.assemble_vector!(b, get_free_dof_values(u1), V)
  ldiv!(RTMMchol, get_free_dof_values(u1))
end

function compute_mass_flux!(F,dΩ,V,RTMMchol,u)
  b(v) = ∫(v⋅u)dΩ
  Gridap.FESpaces.assemble_vector!(b, get_free_dof_values(F), V)
  ldiv!(RTMMchol, get_free_dof_values(F))
end

function compute_depth!(h1,dΩ,dω,Q,L2MMchol,h2,F,dt)
  b(q)  = ∫(q*h2)dΩ - ∫(dt*q*DIV(F))dω
  Gridap.FESpaces.assemble_vector!(b, get_free_dof_values(h1), Q)
  ldiv!(L2MMchol, get_free_dof_values(h1))
end

function compute_bernoulli_potential!(ϕ,dΩ,Q,L2MMchol,uu,h,g)
  b(q)  = ∫(q*(0.5*uu + g*h))*dΩ
  Gridap.FESpaces.assemble_vector!(b, get_free_dof_values(ϕ), Q)
  ldiv!(L2MMchol, get_free_dof_values(ϕ))
end

function compute_diagnostic_vorticity!(w,dΩ,S,H1MMchol,u,n)
  b(s) = ∫(⟂(n,∇(s))⋅(u))dΩ
  Gridap.FESpaces.assemble_vector!(b, get_free_dof_values(w), S)
  ldiv!(H1MMchol, get_free_dof_values(w))
end

function shallow_water_explicit_time_step!(
     h₂, u₂, hₚ, uₚ, ϕ, F, q₁, q₂, H1h, H1hchol,      # in/out args
     model, dΩ, dω, V, Q, R, S, f, g, h₁, u₁, hₘ, uₘ, # in args
     RTMMchol, L2MMchol, dt, τ, leap_frog)            # more in args

  # energetically balanced explicit second order shallow water solver
  # reference: eqns (21-24) of
  # https://github.com/BOM-Monash-Collaborations/articles/blob/main/energetically_balanced_time_integration/EnergeticallyBalancedTimeIntegration_SW.tex

  n = get_normal_vector(model)
  # explicit step for provisional velocity, uₚ
  dt1 = dt
  if leap_frog
    dt1 = 2.0*dt
  end

  # 1.1: the mass flux
  compute_mass_flux!(F,dΩ,V,RTMMchol,u₁*h₁)
  # 1.2: the bernoulli function
  compute_bernoulli_potential!(ϕ,dΩ,Q,L2MMchol,u₁⋅u₁,h₁,g)
  # 1.3: the potential vorticity
  compute_potential_vorticity!(q₁,H1h,H1hchol,dΩ,R,S,h₁,u₁,f,n)
  # 1.4: solve for the provisional velocity
  compute_velocity!(uₚ,dΩ,dω,V,RTMMchol,uₘ,q₁-τ*u₁⋅∇(q₁),F,ϕ,n,dt1,dt1)
  # 1.5: solve for the provisional depth
  compute_depth!(hₚ,dΩ,dω,Q,L2MMchol,hₘ,F,dt1)

  # 2.1: the mass flux
  compute_mass_flux!(F,dΩ,V,RTMMchol,u₁*(2.0*h₁ + hₚ)/6.0+uₚ*(h₁ + 2.0*hₚ)/6.0)
  # 2.2: the bernoulli function
  compute_bernoulli_potential!(ϕ,dΩ,Q,L2MMchol,(u₁⋅u₁ + u₁⋅uₚ + uₚ⋅uₚ)/3.0,0.5*(h₁ + hₚ),g)
  # 2.3: the potential vorticity
  compute_potential_vorticity!(q₂,H1h,H1hchol,dΩ,R,S,hₚ,uₚ,f,n)
  # 2.4: solve for the final velocity
  compute_velocity!(u₂,dΩ,dω,V,RTMMchol,u₁,q₁-τ*u₁⋅∇(q₁)+q₂-τ*uₚ⋅∇(q₂),F,ϕ,n,0.5*dt,dt)
  # 2.5: solve for the final depth
  compute_depth!(h₂,dΩ,dω,Q,L2MMchol,h₁,F,dt)
end

function project_shallow_water_initial_conditions(dΩ, Q, V, S, L2MMchol, RTMMchol, H1MMchol, h₀, u₀, f₀)
  b₁(q)   = ∫(q*h₀)dΩ
  rhs1    = assemble_vector(b₁, Q)
  hn      = FEFunction(Q, copy(rhs1))

  b₂(v)   = ∫(v⋅u₀)dΩ
  rhs2    = assemble_vector(b₂, V)
  un      = FEFunction(V, copy(rhs2))

  b₃(s)   = ∫(s*f₀)*dΩ
  rhs3    = assemble_vector(b₃, S)
  f       = FEFunction(S, copy(rhs3))

  hnv,unv,fv=get_free_dof_values(hn,un,f)
  ldiv!(L2MMchol, hnv)
  ldiv!(RTMMchol, unv)
  ldiv!(H1MMchol, fv)

  hn, un, f, hnv, unv, fv
end

function shallow_water_explicit_time_stepper(model, order, degree,
                        h₀, u₀, f₀, g,
                        dt, τ, N;
                        write_diagnostics=true,
                        write_diagnostics_freq=1,
                        dump_diagnostics_on_screen=true,
                        write_solution=false,
                        write_solution_freq=N/10,
                        output_dir="nswe_eq_ncells_$(num_cells(model))_order_$(order)_explicit")

  # Forward integration of the shallow water equations
  Ω = Triangulation(model)
  dΩ = Measure(Ω, degree)
  dω = Measure(Ω, degree, ReferenceDomain())

  # Setup the trial and test spaces
  R, S, U, V, P, Q = setup_mixed_spaces(model, order)

  # assemble the mass matrices
  H1MM, _, L2MM, H1MMchol, RTMMchol, L2MMchol = setup_and_factorize_mass_matrices(dΩ, R, S, U, V, P, Q)

  # Project the initial conditions onto the trial spaces
  hn, un, f, hnv, unv, fv =  project_shallow_water_initial_conditions(dΩ, Q, V, S,
                               L2MMchol, RTMMchol, H1MMchol, h₀, u₀, f₀)

  # work arrays
  h_tmp = copy(hnv)
  w_tmp = copy(fv)
  # build the potential vorticity lhs operator once just to initialise
  bmm(a,b) = ∫(a*hn*b)dΩ
  H1h      = assemble_matrix(bmm, R, S)
  H1hchol  = lu(H1h)

  function run_simulation(pvd=nothing)
    diagnostics_file = joinpath(output_dir,"nswe_diagnostics.csv")

    hm1    = clone_fe_function(Q,hn)
    hm2    = clone_fe_function(Q,hn)
    hp     = clone_fe_function(Q,hn)
    ϕ      = clone_fe_function(Q,hn)

    um1    = clone_fe_function(V,un)
    um2    = clone_fe_function(V,un)
    up     = clone_fe_function(V,un)
    F      = clone_fe_function(V,un)

    wn     = clone_fe_function(S,f)
    q1     = clone_fe_function(S,f)
    q2     = clone_fe_function(S,f)

    # first step, no leap frog integration
    shallow_water_explicit_time_step!(hn, un, hp, up, ϕ, F, q1, q2, H1h, H1hchol,
                                      model, dΩ, dω, V, Q, R, S, f, g, hm1, um1, hm2, um2,
                                      RTMMchol, L2MMchol, dt, τ, false)

    if (write_diagnostics)
      initialize_csv(diagnostics_file,"time", "mass", "vorticity", "kinetic", "potential", "power")
    end

    if (write_diagnostics && write_diagnostics_freq==1)
      compute_diagnostic_vorticity!(wn, dΩ, S, H1MMchol, un, get_normal_vector(model))
      dump_diagnostics_shallow_water!(h_tmp, w_tmp,
                                      model, dΩ, dω, S, L2MM, H1MM,
                                      hn, un, wn, ϕ, F, g, 1, dt,
                                      diagnostics_file,
                                      dump_diagnostics_on_screen)
    end

    # subsequent steps, do leap frog integration
    # (now that we have the state at two previous time levels)
    for istep in 2:N
      h_aux = hm2
      hm2   = hm1
      hm1   = hn
      hn    = h_aux
      u_aux = um2
      um2   = um1
      um1   = un
      un    = u_aux

      shallow_water_explicit_time_step!(hn, un, hp, up, ϕ, F, q1, q2, H1h, H1hchol,
                                        model, dΩ, dω, V, Q, R, S, f, g, hm1, um1, hm2, um2,
                                        RTMMchol, L2MMchol, dt, τ, true)

      if (write_diagnostics && write_diagnostics_freq>0 && mod(istep, write_diagnostics_freq) == 0)
        compute_diagnostic_vorticity!(wn, dΩ, S, H1MMchol, un, get_normal_vector(model))
        dump_diagnostics_shallow_water!(h_tmp, w_tmp,
                                        model, dΩ, dω, S, L2MM, H1MM,
                                        hn, un, wn, ϕ, F, g, istep, dt,
                                        diagnostics_file,
                                        dump_diagnostics_on_screen)
      end
      if (write_solution && write_solution_freq>0 && mod(istep, write_solution_freq) == 0)
        compute_diagnostic_vorticity!(wn, dΩ, S, H1MMchol, un, get_normal_vector(model))
        pvd[dt*Float64(istep)] = new_vtk_step(Ω,joinpath(output_dir,"n=$(istep)"),["hn"=>hn,"un"=>un,"wn"=>wn])
      end
    end
    hn, un
  end
  if (write_diagnostics || write_solution)
    rm(output_dir,force=true,recursive=true)
    mkdir(output_dir)
  end
  if (write_solution)
    pvdfile=joinpath(output_dir,"nswe_eq_ncells_$(num_cells(model))_order_$(order)_explicit")
    paraview_collection(run_simulation,pvdfile)
  else
    run_simulation()
  end
end


text:
using JuMP
using IntervalArithmetic: Interval

# We vendor ConditionalJuMP (https://github.com/rdeits/ConditionalJuMP.jl/blob/e0c406077c0b07be76e02f72c3a7a7aa650df82f/src/ConditionalJuMP.jl)
# so that we can use JuMP >= 0.2.0

owner_model(x::JuMP.VariableRef)::Model = JuMP.owner_model(x)
function owner_model(x::JuMP.GenericAffExpr)::Union{Model,Nothing}
    if length(x.terms) == 0
        return nothing
    end
    return JuMP.owner_model(first(x.terms.keys))
end
function owner_model(
    xs::AbstractArray{T},
)::Model where {T<:Union{JuMP.VariableRef,JuMP.GenericAffExpr}}
    for x in xs
        m = owner_model(x)
        if m !== nothing
            return m
        end
    end
    return nothing
end

interval(x::Number) = Interval(x, x)
interval(x::JuMP.VariableRef) = Interval(lower_bound(x), upper_bound(x))
function interval(e::JuMP.GenericAffExpr)
    result = Interval(e.constant, e.constant)
    for (var, coef) in e.terms
        result += Interval(coef, coef) * Interval(lower_bound(var), upper_bound(var))
    end
    return result
end

lower_bound(x::Number) = x
upper_bound(x::Number) = x
lower_bound(x::JuMP.VariableRef) = JuMP.lower_bound(x)
upper_bound(x::JuMP.VariableRef) = JuMP.upper_bound(x)
lower_bound(e::JuMP.GenericAffExpr) = lower_bound(interval(e))
upper_bound(e::JuMP.GenericAffExpr) = upper_bound(interval(e))
lower_bound(i::Interval) = i.lo
upper_bound(i::Interval) = i.hi


text:
@testset "RFFTSpectrumArray" begin
    function cmparrays(a, b;
        rate_=rate(a), domain_=domain(a), nchannels_=nchannels(a),
        nframes_=nframes(a), names_=names(a), eltype_=eltype(a), typeof_=typeof(a), data_=data(a), ntimeframes_=ntimeframes(a))
        isnothing(rate_) || rate(b) == rate_ || error("rate $(rate(b)) != $(rate_)")
        isnothing(domain_) || domain(b) == domain_ || error("domain $(domain(b)) != $(domain_)")
        isnothing(nchannels_) || nchannels(b) == nchannels_ || error("nchannels $(nchannels(b)) != $(nchannels_)")
        isnothing(nframes_) || nframes(b) == nframes_ || error("nframes $(nframes(b)) != $(nframes_)")
        isnothing(ntimeframes_) || ntimeframes(b) == ntimeframes_ || error("ntimeframes $(ntimeframes(b)) != $(ntimeframes_)")
        isnothing(names_) || names(b) == names_ || error("names $(names(b)) != $(names_)")
        isnothing(eltype_) || eltype(b) == eltype_ || error("eltype $(eltype(b)) != $(eltype_)")
        isnothing(typeof_) || typeof(b) == typeof_ || error("typeof $(typeof(b)) != $(typeof_)")
        isnothing(data_) || data(b) == data_ || error("data not same!")
        true
    end

    sa = RFFTSpectrumArray(rand(Complex{Float64}, 4, 1), 10, 6)
    sc = RFFTSpectrumArray(rand(Complex{Float64}, 4, 2), 10Hz, 6)
    sd = RFFTSpectrumArray(rand(Complex{Float64}, 4, 2), 32Hz, 6)
    se = RFFTSpectrumArray(rand(Complex{Float32}, 4, 2), 32Hz, 6)
    sf = RFFTSpectrumArray([Complex(float(j), float(-j / i)) for i in 1:16, j in 1:8], 44100Hz, 30,
        [:front_left, :front_right, :rear_left, :rear_right,
            :front_center, :lfe, :side_left, :side_right])
    sg = copy(sf)

    @test sf == sg
    @test sf !== sg
    @test cmparrays(sf, sg)

    @test_throws DimensionMismatch RFFTSpectrumArray(rand(Complex{Float64}, 4, 2), 10Hz, 6, [:left, :right, :onemore])

    @testset "basic functions" begin
        @test domain(sa) == 0.0:1.6666666666666667:5.0
        @test domain_no0(sa) == 1.6666666666666667:1.6666666666666667:5.0
        @test data_no0(sa) == data(sa)[2:end, :]
        @test nchannels(sa) == 1
        @test nframes(sa) == 4
        @test ntimeframes(sc) == 6
        @test rate(sa) == 10.0
        @test rate(sd) == 32.0
        @test eltype(sa) == Complex{Float64}
        @test cmparrays(sa, sc, nchannels_=2, names_=[Symbol(1), Symbol(2)], data_=nothing)
        @test names(sf) == [:front_left, :front_right, :rear_left, :rear_right,
            :front_center, :lfe, :side_left, :side_right]

        sc2 = copy(sc)
        names!(sc2, :Right, 2)
        @test_throws ArgumentError names!(sc2, :Right, 1) # non unique name
        names!(sc2, :Left, 1)
        @test cmparrays(sc, sc2; names_=[:Left, :Right])
        names!(sc2, [:L, :R], 1:2)
        @test cmparrays(sc, sc2; names_=[:L, :R])
        names!(sc2, [:left, :right])
        @test_throws ArgumentError names!(sc2, [:same, :same]) # non unique name
        @test cmparrays(sc, sc2; names_=[:left, :right])
        names!(se, [:left, :right])
        @test cmparrays(sc2, se; rate_=32.0, domain_=nothing, data_=nothing, eltype_=Complex{Float32}, typeof_=RFFTSpectrumArray{Complex{Float32}})
        names!(sc2, [:RIGHT, :LEFT], [:right, :left])
        @test cmparrays(sc, sc2; names_=[:LEFT, :RIGHT])
    end

    @testset "indexing" begin
        # linear index
        @test sf[1] == data(sf)[1]
        @test typeof(sf[1]) == eltype(sf)

        @test_throws ArgumentError sf[1:10]
        @test_throws ArgumentError sf[1:10, :]
        @test_throws ArgumentError sf[0Hz..1.75Hz]

        # frame addressing, select all channels
        @test cmparrays(sf, sf[:])

        # channel addressing
        @test cmparrays(sf, sf[:, 2]; nchannels_=1, names_=[:front_right], data_=data(sf)[:, 2:2])
        @test cmparrays(sf[:, 2], sf[:, 2:2])
        @test cmparrays(sf[:, [3, 4, 5, 6, 7]], sf[:, 3:7])
        @test cmparrays(sf[:, 2], sf[:, [2]])
        @test cmparrays(sf[:, 2], sf[:, :front_right])
        @test cmparrays(sf[:, :front_right], sf[:, [:front_right]])
        @test_throws ArgumentError sf[:, :wrong_channel_name]
        @test cmparrays(sf[:, :front_right], sf[:, [:front_right]])
        @test cmparrays(sf[:, :rear_left..:side_left], sf[:, 3:7])
        @test cmparrays(sf[:, [name for name in names(sf) if occursin("left", String(name))]], sf[:, [1, 3, 7]])

        # getting integer indices
        @test [SampleArrays.toindex(sa, t) for t in [0Hz, 0.8Hz, 1Hz, 10Hz]] == [1, 1, 2, 7]
        @test [SampleArrays.toindex(sa, t) for t in [0, 0.8, 1, 10]] == [1, 1, 2, 7]

        @test_throws ArgumentError SampleArrays.toframeidx(sa, 0s..1s)
        @test SampleArrays.toframeidx(sa, 0Hz..3.4Hz) == 1:3
        @test SampleArrays.toframeidx(sa, 0Hz..4.3Hz) == 1:4

        # zero channels
        sc2 = sc[:, []]
        @test nchannels(sc2) == 0
        sc3 = sc[:, [false, false]]
        @test nchannels(sc3) == 0
        @test cmparrays(sc2, sc3)

        # setindex!
        sc[1] = -2
        @test data(sc)[1] == -2

        @test_throws ArgumentError sc[:, :] = -2
        sc[:] .= -2
        @test all(data(sc)[:, :] .== -2)

        sc[:, 1:2] .= -5
        @test all(data(sc)[:, 1:2] .== -5)

        sc[:, 1] .= -13
        @test all(data(sc[:, 1]) .== -13)

        se[:, :left] .= -14
        @test all(data(se[:, 1]) .== -14)

        se[:, [:left, :right]] .= -15
        @test all(data(se[:, 1:2]) .== -15)

        sc[:, [1]] .= -20
        @test all(data(sc)[:, [1]] .== -20)

        sc[:, [false, true]] .= -50
        @test all(data(sc)[:, [false, true]] .== -50)
    end

    @testset "similar" begin
        @test cmparrays(sc, similar(sc); data_=nothing)
        @test cmparrays(sc, similar(sc, Float32); data_=nothing, eltype_=Float32, typeof_=RFFTSpectrumArray{Float32})
        @test cmparrays(sc, similar(sc, eltype(sc), (4, 2)); domain_=nothing, data_=nothing)
        @test cmparrays(sc, similar(sc, Complex{Float32}, (4, 2)); domain_=nothing,
            data_=nothing, eltype_=Complex{Float32}, typeof_=RFFTSpectrumArray{Complex{Float32}})
        names!(sc, [:left, :right])
        @test cmparrays(sc, similar(sc, Complex{Float64}, (4, 1)); domain_=nothing, data_=nothing, nchannels_=1, names_=[:left])
        @test cmparrays(sc, similar(sc, Complex{Float64}, (4, 3)); domain_=nothing, data_=nothing,
            nchannels_=3, names_=SampleArrays._default_channel_names(3))
    end

    @testset "view" begin
        sc2 = copy(sc)
        vsc2 = view(sc2, :, 1)
        typeof(vsc2) <: SubArray

        vsc2[:] .= -5
        @test all(sc2[:, 1:1] .== -5)
    end

    @testset "broadcasting" begin
        @test sa + sa == 2 * sa
        @test data(sa .+ 5) == data(sa) .+ 5

        sa32 = Complex{Float32}.(sa)
        ssum = sa + sa32
        @test ssum ≈ 2sa

        @test cmparrays(ssum, sa32 .+ sc; nchannels_=2, names_=nothing, data_=nothing)
    end

    @testset "concatenation" begin
        @test data(hcat(sa, sc)) == hcat(data(sa), data(sc))
        se2 = hcat(se, se)
        @test cmparrays(se, se2[:, 1:2])
        @test cmparrays(se, se2[:, 3:4]; names_=[:left_2, :right_2])
        se2 = hcat(se, Complex{Float32}.(se))
        @test cmparrays(se, se2[:, 1:2])
        @test cmparrays(se, se2[:, 3:4]; names_=[:left_2, :right_2], data_=nothing)
        @test_throws ArgumentError hcat(sa, se) # different rate

        @test_throws ErrorException vcat(sc, sc)
        @test_throws ErrorException cat(sa, sa; dims=1)
    end

    @testset "phase" begin
        sdz = zerophase(sd)
        @test sdz !== sd
        @test all(imag.(sdz) .== 0)
        sdz2 = copy(sd)
        zerophase!(sdz2)
        @test sdz == sdz2
        @test sdz !== sdz2
    end

    @testset "FFT" begin
        smparray = SampleArray(rand(6, 2), 10Hz)
        sa = rfft(smparray)
        @show sa |> typeof
        @test irfft(sa) ≈ smparray
        @test cmparrays(smparray, irfft(sa); data_=nothing, ntimeframes_=nothing)

        smparray = SampleArray(rand(7, 2), 10Hz)
        sa = rfft(smparray)
        @test irfft(sa) ≈ smparray
        @test cmparrays(smparray, irfft(sa); data_=nothing, ntimeframes_=nothing)

        @test_throws DimensionMismatch irfft(RFFTSpectrumArray(rand(Complex{Float64}, 4, 2), 10Hz, 1000)) # wrong number of original frames

        isf = irfft(sf, nframes(sf) * 2 - 2)
        irfft(sf, nframes(sf) * 2 - 1)

        sf2 = MagPhase.(sf)
        isf_mp = irfft(sf2, nframes(sf2) * 2 - 2)
        @test isf ≈ isf_mp
        @test cmparrays(isf, isf_mp, ntimeframes_=nothing, data_=nothing)
    end

#     @testset "resamplefreq" begin
#     end
end


text:
using Documenter, ImageContainers

makedocs(;
    modules=[ImageContainers],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/Lirimy/ImageContainers.jl/blob/{commit}{path}#L{line}",
    sitename="ImageContainers.jl",
    authors="Lirimy <31124605+Lirimy@users.noreply.github.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/Lirimy/ImageContainers.jl",
)


text:
module LowLevelFloatFunctions

export bitwidth, signbit,
       sign, precision, exponent, biased_exponent, significand,
       sign_bits, exponent_bits, significand_bits,
       exponent_max, exponent_min, exponent_field_max,
       exponent_bias, unbias_exponent, bias_exponent,
       sign_field, exponent_field, significand_field,
       unbiased_exponent_field, biased_exponent_field,
       sign_and_exponent_fields, exponent_and_significand_fields,
       mask_sign, mask_exponent, mask_significand,
       mask_exponent_and_significand, mask_sign_and_significand,
       floating

import Base: signed, unsigned
import Base.Math: precision, significand_bits, exponent_bits, exponent_bias,
                  sign_mask, exponent_mask, significand_mask

const SysFloat = Union{Float64, Float32, Float16}

biased_exponent(x::T) where T<:SysFloat = Int(biased_exponent_field(x))
unbias_exponent(x::T) where T<:SysFloat = Int(x - exponent_bias(T))

# extend coverage to Unsigneds for field processing functions

for F in (:precision, :significand_bits, :exponent_bits)
    for (T,U) in ((:Float64, :UInt64), (:Float32, :UInt32), (:Float16, :UInt16))
        @eval begin
            @inline $F(::Type{$U}) = $F($T)
        end
    end
end

@inline sign_bits(::Type{T}) where T<:SysFloat = 1
@inline sign_bits(::Type{T}) where T<:Union{UInt64, UInt32, UInt16} = 1

@inline bitwidth(::Type{T}) where T = sizeof(T) * 8

@inline exponent_max(::Type{Float16})  =     15
@inline exponent_max(::Type{Float32})  =    127
@inline exponent_max(::Type{Float64})  =   1023

@inline exponent_min(::Type{T}) where T<:SysFloat = 1 - exponent_max(T)

@inline exponent_bias(::Type{UInt16})  =     15%UInt16
@inline exponent_bias(::Type{UInt32})  =    127%UInt32
@inline exponent_bias(::Type{UInt64})  =   1023%UInt64

@inline exponent_field_max(::Type{T}) where T<:SysFloat = exponent_max(T) + one(convert(Unsigned, T))

# extend coverage to Unsigneds for field processing functions

for F in (:exponent_max, :exponent_min, :exponent_field_min)
    for (T,U) in ((:Float64, :UInt64), (:Float32, :UInt32), (:Float16, :UInt16))
        @eval begin
            @inline $F(::Type{$U}) = $F($T)
        end
    end
end

for (F,G) in ((:mask_sign, :sign_mask), (:mask_exponent, :exponent_mask), (:mask_significand, :significand_mask))
    for (T,U) in ((:Float64, :UInt64), (:Float32, :UInt32), (:Float16, :UInt16))
        @eval begin
            @inline $F(x::$T) = $G($T) & reinterpret($U,x)
            @inline $F(x::$U) = $G($T) & x
        end
    end
end

include("convert.jl")
include("def_fields.jl")
include("getset_fields.jl")

end # module LowLevelFloatFunctions


text:
"""
    register!(
        mi::MutualInformationContainer,
        full_image::AbstractArray{T,2},
        fixed::AbstractArray{T,2},
        moving_bbox::AbstractVector{Int},
        max_shift_x::Int,
        max_shift_y::Int,
        buffer::AbstractArray{T,2};
        set_buffer! = (buffer, current_frame, moving_bbox) -> set_buffer!(buffer, current_frame, moving_bbox, max_shift_x, max_shift_y),
        get_buffer_crop = (buffer, moving_bbox, shift_x, shift_y) -> get_buffer_crop(buffer, moving_bbox, shift_x, shift_y, max_shift_x, max_shift_y),
        prefilter_frame_crop! = x -> nothing,
        start_shift_x = 3,
        start_shift_y = 3,
        expand_border = 1,
        expand_increment = 1,
    ) where {T<:Integer}

Calculates the shift that best aligns the `moving_bbox` to the `fixed` image within the `full_image`. At a high level,
this attempts to best match the view of `moving_bbox` inside `full_image` to the `fixed` image. This only considers
shifts along the x and y axes; no rotation is considered. This function incrementally expands the maximum shift horizon
to evaluate as few shifts as possible. At a minimum, all shifts within `± start_shift_x` and `± start_shift_y` are
considered. If the optimal shift falls within `expand_border` of the horizon, the horizon will be expanded by
`expand_increment` and all new shifts will be evaluated. This process repeats until the optimal shift does not fall
within `expand_border` of the horizon, or until the horizon has reached `max_shift_x` and `max_shift_y`.

Adding some padding to `moving_bbox` is a good idea to improve registration stability.
E.g. `my_bbox .+ [-10, -10, 10, 10]`. You will need to determine the best padding value for your data.

The parameter `buffer` is required for temporary storage. Generally, you can define `buffer` as
`Array{T}(undef, (size(fixed) .+ (MAX_SHIFT * 2))...)`. The buffer must have a size which is at least the size of the
`fixed` image expanded by the maximum and minimum x and y shifts on each respective axis. The parameters `set_buffer!`
and `get_buffer_crop` are used to write to and read from this buffer, respectively. There is typically no need to change
these from their defaults.

The parameter `prefilter_frame_crop!` can be specified if you want to apply image filtering before computing the
mutual information between the two images. This function must mutate the image it is given with whatever filtering
operation you implement. Also, the `fixed` image must have the filtering pre-applied. See this package's tests for an
example.
"""
function register!(
    mi::MutualInformationContainer,
    full_image::AbstractArray{T,2},
    fixed::AbstractArray{T,2},
    moving_bbox::AbstractVector{Int},
    max_shift_x::Int,
    max_shift_y::Int,
    buffer::AbstractArray{T,2};
    set_buffer! = (buffer, current_frame, moving_bbox) ->
        set_buffer!(buffer, current_frame, moving_bbox, max_shift_x, max_shift_y),
    get_buffer_crop = (buffer, moving_bbox, shift_x, shift_y) ->
        get_buffer_crop(buffer, moving_bbox, shift_x, shift_y, max_shift_x, max_shift_y),
    prefilter_frame_crop! = x -> nothing,
    start_shift_x = 3,
    start_shift_y = 3,
    expand_border = 1,
    expand_increment = 1,
) where {T<:Integer}
    allowjumps = [-start_shift_x, -start_shift_y, start_shift_x, start_shift_y]
    shift = (0, 0)
    best_mi = 0
    prev_mis::Union{Missing,AbstractArray{Float32,2}} = missing

    while allowjumps[1] >= -max_shift_x &&
              allowjumps[3] <= max_shift_x &&
              allowjumps[2] >= -max_shift_y &&
              allowjumps[4] <= max_shift_y
        shift, best_mi, prev_mis = register!(
            mi,
            full_image,
            fixed,
            moving_bbox,
            allowjumps[1]:allowjumps[3],
            allowjumps[2]:allowjumps[4],
            buffer,
            prev_mis;
            set_buffer!,
            get_buffer_crop,
            prefilter_frame_crop!,
        )

        tooclose = 0

        # Expand the maximum shift by expand_increment on any side that is within the expand_border
        if shift[1] <= allowjumps[1] + expand_border
            allowjumps[1] -= expand_increment
            tooclose += 1
        end
        if shift[1] >= allowjumps[3] - expand_border
            allowjumps[3] += expand_increment
            tooclose += 1
        end
        if shift[2] <= allowjumps[2] + expand_border
            allowjumps[2] -= expand_increment
            tooclose += 1
        end
        if shift[2] >= allowjumps[4] - expand_border
            allowjumps[4] += expand_increment
            tooclose += 1
        end

        # If the shift was not within the expand_border, the shift is a global maximum, so we are done
        if tooclose == 0
            break
        end
    end

    return shift, best_mi, prev_mis
end

"""
    register!(
        mi::MutualInformationContainer,
        full_image::AbstractArray{T,2},
        fixed::AbstractArray{T,2},
        moving_bbox::AbstractVector{Int},
        range_x::AbstractVector{Int},
        range_y::AbstractVector{Int},
        buffer::AbstractArray{T,2},
        prev_mis::Union{Missing,AbstractArray{Float32,2}};
        set_buffer! = (buffer, current_frame, moving_bbox) -> set_buffer!(buffer, current_frame, moving_bbox, maximum(range_x), maximum(range_y)),
        get_buffer_crop = (buffer, moving_bbox, shift_x, shift_y) -> get_buffer_crop(buffer, moving_bbox, shift_x, shift_y, maximum(range_x), maximum(range_y)),
        prefilter_frame_crop! = x -> nothing,
    ) where {T<:Integer}

Calculates the shift that best aligns the `moving_bbox` to the `fixed` image within the `full_image`. At a high level,
this attempts to best match the view of `moving_bbox` inside `full_image` to the `fixed` image. This only considers
shifts along the x and y axes; no rotation is considered. All combinations of shifts within `range_x` and `range_y` are
considered.

Adding some padding to `moving_bbox` is a good idea to improve registration stability.
E.g. `my_bbox .+ [-10, -10, 10, 10]`. You will need to determine the best padding value for your data.

The parameter `buffer` is required for temporary storage between calls to this function because it is assumed that
you will call this function in a loop to register many similar images. Generally, you can define `buffer` as
`Array{T}(undef, (size(fixed) .+ (MAX_SHIFT * 2))...)`. The buffer must have a size which is at least the size of the
`fixed` image expanded by the maximum and minimum x and y shifts on each respective axis. The parameters `set_buffer!`
and `get_buffer_crop` are used to write to and read from this buffer, respectively. There is typically no need to change
these from their defaults.

The parameter `prev_mis` is used to memoize the mutual information calculations if this function is being called
"incrementally" with an expanding shift horizon. If you are calling this function directly (and not in a loop to
gradually expand a maximum shift horizon), you should set this to `missing`, which would cause all combinations of
shifts within `range_x` and `range_y` to be considered.

The parameter `prefilter_frame_crop!` can be specified if you want to apply image filtering before computing the
mutual information between the two images. This function must mutate the image it is given with whatever filtering
operation you implement. Also, the `fixed` image must have the filtering pre-applied. See this package's tests for an
example.
"""
function register!(
    mi::MutualInformationContainer,
    full_image::AbstractArray{T,2},
    fixed::AbstractArray{T,2},
    moving_bbox::AbstractVector{Int},
    range_x::AbstractVector{Int},
    range_y::AbstractVector{Int},
    buffer::AbstractArray{T,2},
    prev_mis::Union{Missing,AbstractArray{Float32,2}};
    set_buffer! = (buffer, current_frame, moving_bbox) ->
        set_buffer!(buffer, current_frame, moving_bbox, maximum(range_x), maximum(range_y)),
    get_buffer_crop = (buffer, moving_bbox, shift_x, shift_y) -> get_buffer_crop(
        buffer,
        moving_bbox,
        shift_x,
        shift_y,
        maximum(range_x),
        maximum(range_y),
    ),
    prefilter_frame_crop! = x -> nothing,
) where {T<:Integer}
    mis = mutual_information!(
        mi,
        fixed,
        buffer,
        full_image,
        moving_bbox,
        range_x,
        range_y,
        prev_mis;
        set_buffer!,
        get_buffer_crop,
        prefilter_frame_crop!,
    )
    best_mi, idx = findmax(mis)
    return idx.I, best_mi, mis
end

function set_buffer!(buffer, current_frame, moving_bbox, max_shift_x, max_shift_y)
    x_inds = (moving_bbox[1]-max_shift_x):(moving_bbox[3]+max_shift_x)
    y_inds = (moving_bbox[2]-max_shift_y):(moving_bbox[4]+max_shift_y)
    buffer .= view(current_frame, x_inds, y_inds)
    nothing
end

function get_buffer_crop(buffer, moving_bbox, shift_x, shift_y, max_shift_x, max_shift_y)
    # Extract the bbox + (shift_x, shift_y). We are extracting it from an array of
    # bbox ± max_shift so we need to transform the coordinates into the right frame by
    # subtracting the origin of the bbox ± max_shift frame.
    x_inds =
        ((moving_bbox[1]:moving_bbox[3]) .+ shift_x) .- (moving_bbox[1] - max_shift_x) .+ 1
    y_inds =
        ((moving_bbox[2]:moving_bbox[4]) .+ shift_y) .- (moving_bbox[2] - max_shift_y) .+ 1
    view(buffer, x_inds, y_inds)
end


text:
using AdaGram
vm, dict = load_model("1bil.test.bin")

println(expected_pi(vm, dict.word2id["new"]))
println(nearest_neighbors(vm, dict, "new", 1, 10))
println(nearest_neighbors(vm, dict, "new", 2, 10))
println(disambiguate(vm, dict, "new", split("york")))
println(disambiguate(vm, dict, "new", split("create")))


text:
using Dates

"""
Write a win file at the indicated prefix.
Parameters to Wannier90 can be added as kwargs: e.g. `num_iter=500`.
"""
function write_w90_win(fileprefix::String, basis::PlaneWaveBasis;
                       bands_plot=false, wannier_plot=false, kwargs...)
    @assert :num_bands in keys(kwargs)  # Required parameter
    @assert :num_wann  in keys(kwargs)  # Required parameter
    isnothing(basis.kgrid) && error("The basis must be constructed from a MP grid.")

    open(fileprefix *".win", "w") do fp
        println(fp, "! Generated by DFTK.jl at $(now())\n")

        # Drop parameters
        for (key, value) in pairs(kwargs)
            @printf fp "%-20s =   %-30s\n" key value
        end
        if wannier_plot
            println(fp, "wvfn_formatted = True")
            println(fp, "wannier_plot   = True")
        end

        # Delimiter for system
        println(fp, "\n", "!"^20 * " System \n\n")

        # Lattice vectors are in rows in Wannier90
        println(fp, "begin unit_cell_cart\nbohr")
        for vec in eachcol(basis.model.lattice)
            @printf fp "%10.6f %10.6f %10.6f \n" vec...
        end
        println(fp, "end unit_cell_cart \n")

        println(fp, "begin atoms_frac")
        for (element, positions) in basis.model.atoms
            for position in positions
                @printf fp "%-2s %10.6f %10.6f %10.6f \n" element.symbol position...
            end
        end
        println(fp, "end atoms_frac\n")

        # Delimiter for k-points block
        println(fp, "!"^20 * " k_points\n")

        if bands_plot
            kpathdata = high_symmetry_kpath(basis.model)
            length(kpathdata.kpath) > 1 || @warn(
                "Only first kpath branch considered in write_w90_win")
            path = kpathdata.kpath[1]

            println(fp, "begin kpoint_path")
            for i in 1:length(path)-1
                A, B = path[i:i+1]  # write segment A -> B
                @printf(fp, "%s %10.6f %10.6f %10.6f  ", A, round.(kpathdata.klabels[A], digits=5)...)
                @printf(fp, "%s %10.6f %10.6f %10.6f\n", B, round.(kpathdata.klabels[B], digits=5)...)
            end
            println(fp, "end kpoint_path")
            println(fp, "bands_plot = true\n")
        end

        println(fp, "mp_grid : $(basis.kgrid[1]) $(basis.kgrid[2]) $(basis.kgrid[3])\n")
        println(fp, "begin kpoints")
        for kpt in basis.kpoints
            @printf fp  "%10.6f %10.6f %10.6f\n" kpt.coordinate...
        end
        println(fp, "end kpoints")
    end
    nothing
end

"""
Read the .nnkp file provided by the preprocessing routine of Wannier90
(i.e. "wannier90.x -pp prefix")
Returns:
1) the array 'nnkpts' of k points, their respective nearest neighbors and associated
   shifing vectors (non zero if the neighbor is located in another cell).
2) the number 'nntot' of neighbors per k point.

TODO: add the possibility to exclude bands
"""
function read_w90_nnkp(fileprefix::String)
    fn = fileprefix * ".nnkp"
    !isfile(fn) && error("Expected file $fileprefix.nnkp not found.")

    # Extract nnkp block
    lines = readlines(fn)
    ibegin_nnkp = findfirst(lines .== "begin nnkpts")
    iend_nnkp   = findfirst(lines .== "end nnkpts")
    @assert !isnothing(ibegin_nnkp) && !isnothing(iend_nnkp)

    nntot  = parse(Int, lines[ibegin_nnkp + 1])
    nnkpts = map(ibegin_nnkp+2:iend_nnkp-1) do iline
        splitted = parse.(Int, split(lines[iline], ' ', keepempty=false))
        @assert length(splitted) == 5

        # Meaning of the entries:
        # 1st: Index of k-point
        # 2nd: Index of periodic image of k+b k-point
        # 3rd: Shift vector to get k-point of ikpb to the actual k+b point required
        (ik=splitted[1], ikpb=splitted[2], G_shift=splitted[3:5])
    end
    (; nntot, nnkpts)
end

"""
Write the eigenvalues in a format readable by Wannier90.
"""
function write_w90_eig(fileprefix::String, eigenvalues; n_bands)
    open("$fileprefix.eig", "w") do fp
        for (k, εk) in enumerate(eigenvalues)
            for (n, εnk) in enumerate(εk[1:n_bands])
                @printf fp "%3i  %3i   %25.18f \n" n k auconvert(Unitful.eV, εnk).val
            end
        end
    end
    nothing
end

@timing function write_w90_unk(fileprefix::String, basis, ψ; n_bands, spin=1)
    fft_size = basis.fft_size

    for ik in krange_spin(basis, spin)
        open(dirname(fileprefix) * (@sprintf "/UNK%05i.%i" ik spin), "w") do fp
            println(fp, "$(fft_size[1]) $(fft_size[2]) $(fft_size[3]) $ik $n_bands")
            for iband in 1:n_bands
                ψnk_real = G_to_r(basis, basis.kpoints[ik], @view ψ[ik][:, iband])
                for iz in 1:fft_size[3], iy in 1:fft_size[2], ix in 1:fft_size[1]
                    println(fp, real(ψnk_real[ix, iy, iz]), " ", imag(ψnk_real[ix, iy, iz]))
                end
            end
        end
    end
    nothing
end

@doc raw"""
Computes the matrix ``[M^{k,b}]_{m,n} = \langle u_{m,k} | u_{n,k+b} \rangle``
for given k, kpb = k+b.

`G_shift` is the "shifting" vector, correction due to the periodicity conditions
imposed on k -> ψ_k.
It is non zero if kpb is taken in another unit cell of the reciprocal lattice.
We use here that : ``u_{n(k + G_shift)}(r) = e^{-i*\langle G_shift,r \rangle} u_{nk}``
"""
@views function overlap_Mmn_k_kpb(basis::PlaneWaveBasis, ψ, ik, ikpb, G_shift, n_bands)
    # Search for common Fourier modes and their resp. indices in Bloch states k and kpb
    # TODO Check if this can be improved using the G vector mapping in the kpoints
    k   = basis.kpoints[ik]
    kpb = basis.kpoints[ikpb]
    equivalent_G_vectors = [(iGk, DFTK.index_G_vectors(basis, kpb, Gk + G_shift))
                            for (iGk, Gk) in enumerate(G_vectors(basis, k))]
    iGk   = [eqG[1] for eqG in equivalent_G_vectors if !isnothing(eqG[2])]
    iGkpb = [eqG[2] for eqG in equivalent_G_vectors if !isnothing(eqG[2])]

    # Compute overlaps
    # TODO This should be improved ...
    Mkb = zeros(ComplexF64, (n_bands, n_bands))
    for n in 1:n_bands
        for m in 1:n_bands
            # Select the coefficient in right order
            Mkb[m, n] = dot(ψ[ik][iGk, m], ψ[ikpb][iGkpb, n])
        end
    end
    iszero(Mkb) && return Matrix(I, n_bands, n_bands)
    Mkb
end

@timing function write_w90_mmn(fileprefix::String, basis::PlaneWaveBasis, ψ, nnkp; n_bands)
    open(fileprefix * ".mmn", "w") do fp
        println(fp, "Generated by DFTK at $(now())")
        println(fp, "$n_bands  $(length(ψ))  $(nnkp.nntot)")
        for (ik, ikpb, G_shift) in nnkp.nnkpts
            @printf fp "%i  %i  %i  %i  %i \n" ik ikpb G_shift...
            for ovlp in overlap_Mmn_k_kpb(basis, ψ, ik, ikpb, G_shift, n_bands)
                @printf fp "%22.18f %22.18f \n" real(ovlp) imag(ovlp)
            end
        end
    end
    nothing
end

@doc raw"""
Compute the matrix ``[A_k]_{m,n} = \langle ψ_m^k | g^{\text{per}}_n \rangle``

``g^{per}_n`` are periodized gaussians whose respective centers are given as an
 (num_bands,1) array [ [center 1], ... ].

Centers are to be given in lattice coordinates and G_vectors in reduced coordinates.
The dot product is computed in the Fourier space.

Given an orbital ``g_n``, the periodized orbital is defined by :
 ``g^{per}_n =  \sum\limits_{R \in {\rm lattice}} g_n( \cdot - R)``.
The  Fourier coefficient of ``g^{per}_n`` at any G
is given by the value of the Fourier transform of ``g_n`` in G.
"""
function compute_Ak_gaussian_guess(basis::PlaneWaveBasis, ψk, kpt, centers, n_bands)
    n_wannier = length(centers)
    # TODO This function should be improved in performance

    # associate a center with the fourier transform of the corresponding gaussian
    fourier_gn(center, qs) = [exp( 2π*(-im*dot(q, center) - dot(q, q) / 4) ) for q in qs]
    qs = vec(map(G -> G .+ kpt.coordinate, G_vectors(basis)))  # all q = k+G in reduced coordinates
    Ak = zeros(eltype(ψk), (n_bands, n_wannier))

    # Compute Ak
    for n in 1:n_wannier
        # Functions are l^2 normalized in Fourier, in DFTK conventions.
        norm_gn_per = norm(fourier_gn(centers[n], qs), 2)
        # Fourier coeffs of gn_per for k+G in common with ψk
        coeffs_gn_per = fourier_gn(centers[n], qs[kpt.mapping]) ./ norm_gn_per
        # Compute overlap
        for m in 1:n_bands
            # TODO Check the ordering of m and n here!
            Ak[m, n] = dot(ψk[:, m], coeffs_gn_per)
        end
    end
    Ak
end

@timing function write_w90_amn(fileprefix::String, basis::PlaneWaveBasis, ψ; n_bands, centers)
    open(fileprefix * ".amn", "w") do fp
        println(fp, "Generated by DFTK at $(now())")
        println(fp, "$n_bands   $(length(basis.kpoints))  $(length(centers))")

        for (ik, (ψk, kpt)) in enumerate(zip(ψ, basis.kpoints))
            Ak = compute_Ak_gaussian_guess(basis, ψk, kpt, centers, n_bands)
            for n in 1:size(Ak, 2)
                for m in 1:size(Ak, 1)
                    @printf(fp, "%3i %3i %3i  %22.18f %22.18f \n",
                            m, n, ik, real(Ak[m, n]), imag(Ak[m, n]))
                end
            end
        end
    end
    nothing
end

"""
Default random Gaussian guess for maximally-localised wannier functions
generated in reduced coordinates.
"""
default_wannier_centres(n_wannier) = [rand(1, 3) for _ in 1:n_wannier]

@timing function run_wannier90(scfres;
                               n_bands=size(scfres.ψ[1], 2) - scfres.n_ep_extra,
                               n_wannier=n_bands,
                               centers=default_wannier_centres(n_wannier),
                               fileprefix=joinpath("wannier90", "wannier"),
                               wannier_plot=false, kwargs...)
    # TODO None of the routines consider spin at the moment
    @assert scfres.basis.model.spin_polarization in (:none, :spinless)
    @assert length(centers) == n_wannier

    # Undo symmetry operations to get full k-point list
    scfres_unfold = unfold_bz(scfres)
    basis = scfres_unfold.basis
    ψ = scfres_unfold.ψ

    # Make wannier directory ...
    dir, prefix = dirname(fileprefix), basename(fileprefix)
    mkpath(dir)

    # Write input file and launch Wannier90 preprocessing
    write_w90_win(fileprefix, basis;
                  num_wann=length(centers), num_bands=n_bands, wannier_plot, kwargs...)
    wannier90_jll.wannier90(exe -> run(Cmd(`$exe -pp $prefix`; dir)))
    nnkp = read_w90_nnkp(fileprefix)

    # Files for main wannierization run
    write_w90_eig(fileprefix, scfres_unfold.eigenvalues; n_bands)
    write_w90_amn(fileprefix, basis, ψ; n_bands, centers)
    write_w90_mmn(fileprefix, basis, ψ, nnkp; n_bands)

    # Writing the unk files is expensive (requires FFTs), so only do if needed.
    if wannier_plot
        write_w90_unk(fileprefix, basis, ψ; n_bands, spin=1)
    end

    # Run Wannierisation procedure
    @timing "Wannierization" wannier90_jll.wannier90(exe -> run(Cmd(`$exe $prefix`; dir)))
    fileprefix
end