text:
# # Introduction to GPU Programming

# The following tutorials assume that you have have setup [CUDA.jl](https://github.com/JuliaGPU/CUDA.jl). Detailed installation instructions can be found [here](https://juliagpu.gitlab.io/CUDA.jl/installation/overview/#InstallationOverview).

# You may check if your `CUDA.jl` installation is functional using

## using Pkg
## Pkg.add("CUDA")
using CUDA
CUDA.functional()

#=
If `CUDA.functional()` returns false then the package is in a non-functional state and you should follow the [documentation](https://juliagpu.gitlab.io/CUDA.jl/installation/overview/) to get it working.

Also explained in the [usage](https://juliagpu.gitlab.io/CUDA.jl/usage/overview/) section of the `CUDA.jl` [docs](https://juliagpu.gitlab.io/CUDA.jl/) is an overview of `CUDA.jl`'s functionality which can work at three distinct levels.

* Array Abstractions: With the help of the `CuArray` type we can use Base's array abstractions like broadcasting and mapreduce.
* Native Kernels: Write kernels which compiles to native GPU code directly from Julia.
* CUDA API wrappers: Call CUDA libraries directly from Julia for bleeding edge performance.

The purpose of these tutorials is to teach you effective GPU programming. The tutorials here complement other GPU programming resources such as the NVIDIA blogs, other online resources and formal textbooks. Using other resources in your study will complement these tutorials and is highly encouraged. Some suggested resources are listed at the end of this tutorial.

A GPU (graphical processing unit) is a device specially designed for graphics work. Graphical tasks are a good candidate for parallelization and  GPU's exploit it by having a large number of less powerful processors instead of a single very powerful processor. In 2007 NVIDIA released CUDA (Compute Unified Device Architecture), a parallel programming platform (hardware and software stack) which alongside graphics also focusses scientific computation. Modern GPU's are commonly called GPGPU (general purpose GPU) which shows their importance in scientific computation alongside graphics.

Programs which execute on the GPU are vastly different due to its different architecture. There are new paradigms and algorithms to learn. Understanding how a GPU works is crucial to maximizing the performance of your application.
=#

#=
## Parallelizing AXPY

[Basic Linear Algebra Subroutines(BLAS)](https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms) are subroutines for Linear Algebra operations. Linear algebra's importance in scientific computing makes BLAS essential to GPU computing. One of the most primitive BLAS operations is to scale a vector and add it to another vector. Given two vectors ($x$ and $y$) and a scalar ($\alpha$) we add $\alpha\cdot x$ to $y$. In BLAS libraries this manifests as the functions SAXPY, DAXPY and CAXPY. The difference between the three is that the data type of the vectors is `Float32`, `Float64` and `Complex{Float32}` respectively. However in this example we call our subroutine `axpy` and let Julia take care of the types.
=#

function axpy!(A, X, Y)
    for i in eachindex(Y)
        @inbounds Y[i] = A * X[i] + Y[i]
    end
end

N = 2^27
v1 = rand(Float32, N)
v2 = rand(Float32, N)
v2_copy = copy(v2) # maintain a copy of the original
α = rand()

axpy!(α, v1, v2)

#=
Alternatively, we can also use Julia's [broadcasting](https://docs.julialang.org/en/v1/manual/arrays/#Broadcasting) syntax which allows us to write it in simpler and equally performant version.
=#

v3 = copy(v2_copy)
v3 .+= α * v1

@show v2 == v3

#=
#### CPU multithreaded version

Consider parallelization on a CPU with `p` processors. We can divide our arrays into `p` subarrays of equal size and assign a processor to each subarray. This can theoretically make our parallel version `p` times faster. We say "theoretically" because there is an overhead of starting threads and synchronizing them. Our hope in parallel computing is that the cost will get amortized with the speedup of parallelization, but that may not be the case. Which is why measuring performance is extremely important. Nevertheless, the parallel version asymptotically scales linearly w.r.t `p` which is really good, so much so that these types of problems are classified as "embarassingly parallel". In other cases when processors need to communicate and synchronize frquently the benefit does not scale linearly with the number of processors.

We can use Julia's inbuilt multithreading functionality to use multiple CPU threads which is documented([here](https://docs.julialang.org/en/v1/manual/multi-threading/)). You need to ensure that Julia starts with the appropriate number of threads using the environment variable or startup option(`-t NUMTHREADS`), instructions for which are given in the docs.

A common theme in parallel computing is the concept of thread rank or id. Each thread has a unique id/rank which helps us identify them and map them to tasks easily.
=#

using Base.Threads

println("Number of CPU threads = ", nthreads())

## pseudocode for parallel saxpy
function parallel_axpy!(A, X, Y)
    len = cld(length(X), nthreads())

    ## Launch threads = nthreads()
    Threads.@threads for i in 1:nthreads()
        ## set id to thread rank/id
        tid = threadid()
        low = 1 + (tid - 1)*len
        high = min(length(X), len * tid) ## The last segment might have lesser elements than len

        ## Broadcast syntax, views used to avoid copying
        view(Y, low:high) .+= A.*view(X, low:high)
    end
    return
end

v4 = copy(v2_copy)
parallel_axpy!(α, v1, v4)

@show v2 == v4

#=
#### GPU version

Given below is the code for GPU
=#

function gpu_axpy!(A, X, Y)
    ## set tid to thread rank
    tid = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    tid > length(Y) && return
    @inbounds Y[tid] = A*X[tid] + Y[tid]
    return
end

## Transfer array to GPU memory
gpu_v1 = CuArray(v1)
gpu_v2 = CuArray(v2_copy)

numthreads = 256
numblocks = cld(N, numthreads)

@show numthreads
@show numblocks

## Launch the gpu_axpy! on the GPU
@cuda threads=numthreads blocks=numblocks gpu_axpy!(α, gpu_v1, gpu_v2)

## Copy back to RAM
v4 = Array(gpu_v2)

## Verify that the answers are the same
@show v2 == v4

#=
Compared to the CPU code there are a number of differences:

##### 1) Thread Indexing

The multithreaded CPU code used `threadid()` to get the current thread's rank whereas on the GPU a complicated expression `tid = (blockIdx().x - 1) * blockDim().x + threadIdx().x` computed rank. Furthermore, we are using two distinct terms, `blocks` and `threads` coupled with `idx`(index) and `Dim` (dimension).

##### 2) SIMT architecture

The multithreaded CPU code divided the array up into a handful of pieces equal to the number of processors. A modern consumer CPU has a handful of cores(4 - 8), hence each thread still works on a relatively large array whereas the GPU processes one element per thread. While both demonstrate parallelism their scales differ vastly.

[Flynn's taxonomy](https://hpc.llnl.gov/tutorials/introduction-parallel-computing/flynns-classical-taxonomy) is a popular way to classify parallel computer architectures.

|                      | single data | multiple data |
|:---------------------|:------------|:--------------|
| single instruction   | SISD        | SIMD          |
| multiple instruction | MISD        | MIMD          |

* SISD(single instruction single data) is the classical uniprocessor model. A single instruction stream executes, acting on a single data element at a time.
* SIMD(single instruction multiple data) incorporates a level of parallelism by having a single instruction stream acting on multiple data elements at a time. An example of this is vectorized CPU instructions which use large registers containing multiple data elements. Instructions that work with these large vector registers effectively work on multiple data elements in parallel with a single instruction by utilizing special hardware.
![SIMD AVX-2 4x64i addition](../assets/simd.jpg)
* MISD (multiple instruction single data) is currently only a theoretical model and no commercial machine has been built which uses it.
* MIMD (multiple instruction multiple data) is able to manage multiple instruction streams and acts on multiple data elements at the same time. The CPU multithreading model belongs to it. Each processor can work independently using a different instruction stream acting on different data as required.

To describe CUDA's parallel model NVIDIA coined the term SIMT (single instruction multiple threads) as an extention to SIMD classification. Just like a SIMD vector packs a fixed number of data elements in a wide register, a GPU packs a fixed number of threads in a single warp. Currently NVIDIA packs 32 threads in a single warp and AMD cards pack 64 threads. For more details refer to NVIDIA's docs [here](https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#hardware-implementation).

##### 3) Memory
CPU's memory (RAM) and GPU memory are distinct and is called *host* (CPU) and *device* (GPU) memory respectively. In Julia we need to explicitly transfer memory to and from GPU memory. The reason for this is that copying memory is an expensive operation with high latency. GPU's use PCIe lanes to transfer memory to and from RAM. Poor usage of memory transfers is detrimental to performace and can easily negate all benefits of using a GPU.

Device code referencing CPU memory will result in errors. Host code referencing device memory is produces a warning.
=#
arr = CUDA.rand(10);
arr[1]

#=
To disallow scalar operations altogether use the `CUDA.allowscalar()` function.
=#

CUDA.allowscalar(false)
try
    arr[1]
catch e
    println("Error Caught: ", e)
end

#=
To temporarily allow it in an experssion use the `@allowscalar` macro. However it is suggested that once your application executes correctly on the GPU, you should disallow scalar indexing and use GPU-friendly array operations instead. Accessing GPU memory in a scalar fashion is extremely detrimental to performance.
=#

CUDA.@allowscalar arr[1]

#=
A GPU also has different types of memory such as global memory, texture memory, constant memory which will be discussed later. In general what we call *global memory* is the GPUs DRAM which can be accessed by all threads and is what will be used most often. Memory transfers between the host and device involve the GPUs global memory.

You can check your GPU's memory using `CUDA.available_memory()` and `CUDA.total_memory()` which returns the number of bytes.
=#

@show CUDA.available_memory()
@show CUDA.total_memory();

#=
##### 4) Kernel

When we used the `@cuda` macro, it compiled the `gpu_saxpy!` function for execution on the GPU. A GPU has it's own [instruction set](https://simple.wikipedia.org/wiki/Instruction_set) just like a CPU. The compiled function is called the **kernel** and is sent to the GPU for execution. Once sent we can either wait for the GPU to complete execution or work on something different while it is executing. This can be done using the `blocking` option.

Although CUDA's native instruction set is proprietary there are other ways to inspect code at various stages of compilation. The [reflection page](https://juliagpu.gitlab.io/CUDA.jl/api/compiler/#Reflection) of the documentation should be consulted.

As an example consider `PTX` which resembles low level RISC-ISA like code. PTX is commonly used to inspect code and NVIDIA's [PTX docs](https://docs.nvidia.com/cuda/parallel-thread-execution/index.html#introduction) explains it well.
=#

@device_code_ptx @cuda threads=numthreads blocks=numblocks gpu_axpy!(α, gpu_v1, gpu_v2)

#=
#### Measuring Time
Since the primary inspiration for parallel programming is performance it is important to measure it effectively. Because the CPU and GPU can execute asynchronously there some nuance in their profiling. When we launch a CUDA kernel using `@cuda` after the kernel is launched, control is immediately returned back to the CPU. The CPU can continue executing other code until it's forced to synchronize with the GPU. Certain events like memory transfers and kernel launches can force synchronization. While measuring time and benchmarking we need to force synchronization otherwise we are measuring the time to launch kernels rather than the time it took to execute on the GPU.

Two simple ways to force synchronization are to use the `CUDA.@sync ex` where the CPU is blocked until `ex` finishes execution. The other is to use `CUDA.@time` which synchronizes before and after `ex`. Using `CUDA.@sync` is advisable when using a benchmarking package like [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl). Example `@benchmark CUDA.@sync ex`

Another way is to use [CUDA Events](https://juliagpu.gitlab.io/CUDA.jl/lib/driver/#Event-Management) which can be used in scenarios where a number of events and their statistics are to be collected.

Finally, having a look at NVIDIA's benchmarking tools like Nsight Systems and Nsight Compute can be very helpful in understanding an applications timeline and individual kernel performance. Both of these will be discussed in future tutorials.
=#

@time axpy!(α, v1, v2)
@time parallel_axpy!(α, v1, v2)
@time @cuda threads=numthreads blocks=numblocks gpu_axpy!(α, gpu_v1, gpu_v2)
sleep(0.1) ## Wait for the previous function to finish
@time CUDA.@sync @cuda threads=numthreads blocks=numblocks gpu_axpy!(α, gpu_v1, gpu_v2)
CUDA.@time @cuda threads=numthreads blocks=numblocks gpu_axpy!(α, gpu_v1, gpu_v2)

#=
Notice how the time with `@time @cuda` is much lesser than the `@time CUDA.@sync` and `CUDA.@time` counterparts.

## GPU Architecture

A GPU is made up of an array of *Streaming Multi-Processors*(SM) connected to *Global Memory*. Each streaming multiprocessor consists of warp schedulers, a register file and functional units like single/double precision ALU, Load-Store units,.etc to execute multiple warps concurrently. Effectively hundreds of threads can be executed concurrently on a single SM. Performance of a GPU scales with the number of SM's it has.

![GPU Architecture](../assets/GPU_diagram.png)

![SM Architecture](../assets/SM_diagram.png)

When a kernel is launched on a GPU we also specify a grid configuration using  the `blocks` and `threads` arguments. A grid is composed of "thread blocks": a logical collection of threads. The `blocks` argument defines the block configuration for the grid and the `threads` argument defines the thread configuration for the thread block.

The GPU schedules each thread block to any available SM with sufficient resources. Blocks can be processed in **any** order by the GPU. Multiple thread blocks may execute on a single SM if sufficient resources are available. As thread blocks complete execution other thread blocks take their place.

Each thread block contains a *cooperative thread array*(CTA) which is specified by the `threads` argument. Threads which belong to the same CTA can easily communicate and coordinate with each other because they belong to the same SM. They also have access to a shared memory which is much faster than global memory. The maximum size of a CTA is currently 1024 on NVIDIA hardware.

A small summary of some of the new terms we came across.
- *thread warp*: A set of threads with a fixed size(32). Instructions in a warp are executed together.
- *thread block*: A logical collection of threads which can communicate and coordinate easily.
- *grid*: A logical collection of thread blocks.

#### Resources for learning

Most learning material uses the C/C++ flavor of CUDA. However, there aren't any significant differences and most material can easily be translated from C/C++ to Julia. The main objective is to understand the programming model and how a GPU works rather than the syntax.

1. [CUDA C Programming Guide](https://docs.nvidia.com/cuda/archive/9.1/pdf/CUDA_C_Programming_Guide.pdf): A reference to the CUDA platform with details on hardware and CUDA C/C++ features.
2. [NVIDIA Developer Blog](https://developer.nvidia.com/blog/): Contains many educational blogposts such as [reduction](https://developer.nvidia.com/blog/faster-parallel-reductions-kepler/), [Matrix Multiply](https://developer.nvidia.com/blog/efficient-matrix-transpose-cuda-cc/),.etc
3. [Programming Massively Parallel Processors: A Hands-on Approach (David Kirk and Wen-mei Hwu)](https://www.elsevier.com/books/programming-massively-parallel-processors/kirk/978-0-12-811986-0): A formal textbook on parallel programming and GPU programming using CUDA.
4. Computer Architecture: A quantitative approach (6th Edition, Chapter 4): Exploring the hardware side of parallel technologies like vector architectures, vectorised SIMD instructions and GPU's.
=#


text:
#start boundary conditions
println("calculating right-hand-side...\n")

F = zeros(dofs,1)
#find indexes of edge that are to be bounded
nbc = count(x->x==-1,bc)+1
ix = findall(x->x==-1,bc)
ix = [ix[i][1] for i=1:nbc-1]
prepend!(ix,0)
ix = [collect(ix[i-1]+1:ix[i]-1) for i=2:nbc]
#ix now contains all indexes, seperated so you can now iterate ix[i], for each i you get the list of edges that are on boundarycondition group i
global bc,edge,elem,nodes,i,F
i=1
#iterate through groups
f(x) = [[0 0]', [0 1e6]']
while true
    if i>nbc-1
        break
    end
    global bc,edge,elem,nodes,rows,xedge,yedge,mapx,i,F
    edix = bc[ix[i]] #contains edge id
    ed = edge[edix,:]#contains edge nodes id
    rows = [i for k=1:length(ix[i]) for i=1:length(elem[:,1]) if issubset(ed[k,1],elem[i,:]) & issubset(ed[k,2],elem[i,:]) & issubset(ed[k,3],elem[i,:])] #contains superset id element
    xedge = nodes[elem[rows,:],:][:,:,1]; yedge = nodes[elem[rows,:],:][:,:,2] #contains x,y of edges
    #zero the non edge indexes of element
    superelem = elem[rows,:]
    mapx = [sum(dims=2,superelem[i,:] .== ed[i,:]') for i=1:length(superelem[:,1])]
    for k=1:length(mapx)
        global bc,edge,elem,nodes,rows,xedge,yedge,mapx,i,F
        g(x) = f(x)[i]
        edix = bc[ix[i]] #contains edge id
        ed = edge[edix,:]#contains edge nodes id
        a = h .* BC(g,xedge[k,:],yedge[k,:],mapx[k])
        id = [2 .*ed (2 .* ed .-1)]
        idlocal = [collect(zip(mapx[k].==1,mapx[k].==1))[i][j] for i=1:6 for j=1:2]
        F[id] .= F[id] .+ a[idlocal]'
    end
    i = i+1
end

ext: jl

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)
# The original type did not include beta_maxent vector
struct EPAlg{T<:Real}
    alpha::T
    beta_vec::AbstractVector{T}
    minvar::T
    maxvar::T
    epsconv::T
    damp::T
    maxiter::Int
    verbose::Bool
end


text:
using Distributed
using LightXML
import GZip
@everywhere import GZip
using Printf

mentalist_shared_data_path = string(ENV["HOME"], "/", ".mentalist")

function _get_first_line(node)
  return split(content(node), "\n")[1]
end

function _older_than_a_day(file)
  return (time() - stat(file).ctime)/86400 > 1
end

function _pubmlst_xml(download_dir=mentalist_shared_data_path)
  if !isdir(download_dir)
    mkdir(download_dir)
  end
  dbases_file = string(download_dir, "/", "dbases.xml")
  if !isfile(dbases_file) || _older_than_a_day(dbases_file)
    @info("Downloading the MLST database xml file...")
    download("https://pubmlst.org/data/dbases.xml", dbases_file)
  end
  # get
  return LightXML.parse_file(dbases_file)
end

function _cgmlst_http(download_dir=mentalist_shared_data_path)
  if !isdir(download_dir)
    mkdir(download_dir)
  end
  cgmlst_file = string(download_dir, "/", "cgmlst.html")
  if !isfile(cgmlst_file)
    @info("Downloading the cgmlist HTML to find schema...")
    download("www.cgmlst.org/ncs", cgmlst_file)
  end
  return cgmlst_file
end

function list_pubmlst_schema(prefix)
  xroot = root(_pubmlst_xml())
  scheme_list = Tuple{Int,String}[]
  for (id,species) in enumerate(xroot["species"])
    sp_name = _get_first_line(species)
    if prefix == nothing || startswith(sp_name, prefix)
      push!(scheme_list,(id,sp_name))
    end
  end
  println("#id\torganism")
  for (id, sp_name) in scheme_list
    @printf "%d\t%-30s\n" id sp_name
  end
  @info("$(length(scheme_list)) scheme(s) found.")
end

@everywhere function _download_to_folder(url, output_dir, overwrite=false, filename=nothing)
  println("Trying to download: $url")
  filepath = joinpath(output_dir, filename == nothing ? basename(url) : filename)
  if overwrite || (!isfile(filepath) || _older_than_a_day(filepath))
    mkpath(output_dir)
    download(url, filepath)
  end
  return filepath
end

function _find_publmst_species(xml_xroot, target_id)
  for (id, species) in enumerate(xml_xroot["species"])
    sp_name = _get_first_line(species)
    if target_id == "$id" || target_id == sp_name
      return(species)
    end
  end
end

function download_pubmlst_scheme(target_species, output_dir, overwrite=false)
  loci_files = String[]
  xroot = root(_pubmlst_xml())
  @info("Searching for the scheme ... ")
  species = _find_publmst_species(xroot, target_species)
  if species == nothing
    @warn("I did not find this scheme on pubmlst, please check the species spelling or the ID and try again.")
    exit(-1)
    return
  end
  @info("Downloading scheme for $(_get_first_line(species)) ... ")
  db = species["mlst"][1]["database"][1]
  @info("Downloading profile ...")
  profile_url = content(db["profiles"][1]["url"][1])
  profile_file = _download_to_folder(profile_url, output_dir, overwrite)
  for locus in db["loci"][1]["locus"]
    locus_name = _get_first_line(locus)
    @info("Downloading locus $locus_name ...")
    locus_url = content(locus["url"][1])
    filepath = _download_to_folder(locus_url, output_dir)
    push!(loci_files, filepath)
  end
  @info("Finished downloading.")
  return (loci_files, profile_file)
end

function available_cgmlst_schema()
  scheme_list = Tuple{String, String}[]
  open(_cgmlst_http()) do f
    for l in eachline(f)
      for row in split(l, "")
        m = match(r"<a href='https?://www.cgmlst.org/ncs/schema/(?<id>[0-9]+)/'>(?<species>.+)", row)
        if m != nothing
          id, species = m.captures
          push!(scheme_list,(id,species))
        end
      end
    end
  end
  return scheme_list
end

function list_cgmlst_schema(prefix)
  println("#id\torganism")
  count = 0
  for (id, species) in available_cgmlst_schema()
    if prefix == nothing || startswith(species,prefix)
      count += 1
      @printf "%s\t%-30s\n" id species
    end
  end
  @info("$count schema found.")
end

function _find_cgmlst_id(target_id)
  for (id,species) in available_cgmlst_schema()
    if id == target_id || species == target_id
      return id
    end
  end
  return nothing
end

@everywhere function _download_enterobase_locus(locus, output_dir, url_folder)
  fasta_locus = joinpath(output_dir, "$locus.fa")
  if !isfile(fasta_locus)
    gzip_locus = _download_to_folder("https://enterobase.warwick.ac.uk/schemes/$url_folder/$locus.fasta.gz", output_dir, false, "$locus.fa.gz")
    # gunzip to a FASTA and remove the gzip file;
    f_in = GZip.open(gzip_locus)
    f_out = open(fasta_locus, "w")
    while !eof(f_in)
      write(f_out, readline(f_in))
    end
    close(f_in)
    close(f_out)
    rm(gzip_locus)
  end
  return fasta_locus
end

function download_enterobase_scheme(scheme, s_type, output_dir, overwrite=false)
  verbose = Dict("S"=>"Salmonella", "Y"=>"Yersinia", "E"=>"Escherichia/Shigella", "C" => "clostridium")
  sp = Dict("E"=>"Escherichia", "S"=>"Salmonella", "Y"=>"Yersinia", "C" => "clostridium")

  if !haskey(sp, scheme)
    @info("Scheme has to be E, S, Y or C.")
    exit(-1)
  end
  # Version: empty for wg, v2 for Salmonella cg, v1 for other cg.
  version = s_type == "wg" ? "" : (scheme == "S" ? "v2" : "v1")
  scheme_name = "$(s_type)MLST$(version)"
  url_folder = "$(sp[scheme]).$(scheme_name)"
  filename = joinpath(dirname(@__FILE__), "../scripts/$(sp[scheme]).txt")

  loci = String[]
  open(filename) do f
    for l in eachline(f)
      values = split(strip(l),"\t")
      if startswith(values[1], s_type) # add locus if it is cg or wgmlst, according to input option
        push!(loci, split(values[2],' ')[1])
      end
    end
  end
  loci_files = pmap(locus->_download_enterobase_locus(locus, output_dir, url_folder), loci)
  return loci_files
end

function download_cgmlst_scheme(target_id, output_dir)
  id = _find_cgmlst_id(target_id)
  if id == nothing
    @warn("Id/species ($target_id) not found!")
    exit(-1)
  end
  @info("Downloading cgMLST scheme ...")
  scheme_zip_file = _download_to_folder("https://www.cgmlst.org/ncs/schema/$id/alleles", output_dir)
  locus_files = String[]
  locus = ""
  @info("Unzipping cgMLST scheme into individual FASTA files for each locus ...")
  scheme_dirname = dirname(scheme_zip_file)
  run(`unzip -oq $scheme_zip_file -d $scheme_dirname/tmp`)
  rm(scheme_zip_file)
  scheme_files = readdir(joinpath(scheme_dirname, "tmp"))

  for scheme_file in scheme_files
    # print to show progress:
    if length(locus_files) % 200 == 0
        print(".")
    end
    scheme_file_path = joinpath(scheme_dirname, scheme_file)
    push!(locus_files, scheme_file_path)
    # get locus ID from filename
    locus = split(scheme_file, ".")[1]
    fh = open(scheme_file_path, "w")
    for l in eachline(joinpath(scheme_dirname, "tmp", scheme_file))
      if length(strip(l)) == 0
        continue
      end
      if l[1] == '>'
        l = '>' * locus * "_" * l[2:end] * "\n"
        write(fh, l)
      else
        write(fh, l * "\n")
      end
    end
    close(fh)
  end
  println()
  total_loci = length(locus_files)
  @info("$total_loci loci found.")
  rm(joinpath(scheme_dirname, "tmp"), recursive=true)
  return locus_files
end

max_stars_repo_name: liberuser/MentaLiST

max_issues_repo_name: liberuser/MentaLiST

max_forks_repo_name: liberuser/MentaLiST


text:
module nucleus

const proton_mass = 938.27209 # [MeV/c^2]
const neutron_mass = 939.56563 # [MeV/c^2]
const c = 299792458 # [m/s]

struct Atom
    neutrons::Integer
    protons::Integer
    function Atom(n::Integer, p::Integer)
        n >= 0 || throw(ArgumentError("Number of neutrons $n must be non-negative"))
        p >= 0 || throw(ArgumentError("Number of protons $p must be non-negative"))
        new(n, p)
    end
end

const a_V = 15.75 # [MeV]
const a_S = 17.8 # [MeV]
const a_C = 0.711 # [MeV]
const a_A = 23.7 # [MeV]
const a_P = 11.18 # [MeV]
const k_P = -1 / 2

function δ(N::T, Z::T, A::T) where {T <: Integer}
    δ_0 = a_P * A^k_P

    if iseven(N) && iseven(Z)
        return +1 * δ_0
    end
    if isodd(A)
        return 0
    end
    if isodd(N) && isodd(Z)
        return -1 * δ_0
    end
end

δ(N::T, Z::T) where {T <: Integer} = δ(N, Z, N + Z)

function E_B(N::T, Z::T, A::T) where {T <: Integer}
    return a_V * A - a_S * A^(2 / 3) - a_C * (Z * (Z - 1)) / A^(1 / 3) - a_A * (A - 2 * Z)^2 / A - δ(N, Z)
end

E_B(N::T, Z::T) where {T <: Integer} = E_B(N, Z, N + Z)
E_B(a::Atom) = E_B(a.neutrons, a.protons)

function m(a::Atom)
    (N, Z) = [a.neutrons, a.protons]
    return N * neutron_mass + Z * proton_mass  - E_B(a) / c^2
end

mev_over_c2_to_u(x) = x * 0.00107354411

uranium_235 = Atom(143, 92)
println(mev_over_c2_to_u(m(uranium_235)))
println(E_B(uranium_235))

end # module

max_stars_repo_name: Pocket-titan/nucleus

max_issues_repo_name: Pocket-titan/nucleus

max_forks_repo_name: Pocket-titan/nucleus


text:
function StreamingStatistics()
    return StreamingStatistics(())
end

function StreamingStatistics(arg0::jboolean)
    return StreamingStatistics((jboolean,), arg0)
end

function accept(obj::StreamingStatistics, arg0::jdouble)
    return jcall(obj, "accept", void, (jdouble,), arg0)
end

function add_value(obj::StreamingStatistics, arg0::jdouble)
    return jcall(obj, "addValue", void, (jdouble,), arg0)
end

function aggregate(obj::StreamingStatistics, arg0::StreamingStatistics)
    return jcall(obj, "aggregate", void, (StreamingStatistics,), arg0)
end

function builder(::Type{StreamingStatistics})
    return jcall(StreamingStatistics, "builder", StreamingStatistics_StreamingStatisticsBuilder, ())
end

function clear(obj::StreamingStatistics)
    return jcall(obj, "clear", void, ())
end

function copy(obj::StreamingStatistics)
    return jcall(obj, "copy", StreamingStatistics, ())
end

function equals(obj::StreamingStatistics, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_geometric_mean(obj::StreamingStatistics)
    return jcall(obj, "getGeometricMean", jdouble, ())
end

function get_max(obj::StreamingStatistics)
    return jcall(obj, "getMax", jdouble, ())
end

function get_mean(obj::StreamingStatistics)
    return jcall(obj, "getMean", jdouble, ())
end

function get_median(obj::StreamingStatistics)
    return jcall(obj, "getMedian", jdouble, ())
end

function get_min(obj::StreamingStatistics)
    return jcall(obj, "getMin", jdouble, ())
end

function get_n(obj::StreamingStatistics)
    return jcall(obj, "getN", jlong, ())
end

function get_percentile(obj::StreamingStatistics, arg0::jdouble)
    return jcall(obj, "getPercentile", jdouble, (jdouble,), arg0)
end

function get_population_variance(obj::StreamingStatistics)
    return jcall(obj, "getPopulationVariance", jdouble, ())
end

function get_quadratic_mean(obj::StreamingStatistics)
    return jcall(obj, "getQuadraticMean", jdouble, ())
end

function get_second_moment(obj::StreamingStatistics)
    return jcall(obj, "getSecondMoment", jdouble, ())
end

function get_standard_deviation(obj::StreamingStatistics)
    return jcall(obj, "getStandardDeviation", jdouble, ())
end

function get_sum(obj::StreamingStatistics)
    return jcall(obj, "getSum", jdouble, ())
end

function get_sum_of_logs(obj::StreamingStatistics)
    return jcall(obj, "getSumOfLogs", jdouble, ())
end

function get_sum_of_squares(obj::StreamingStatistics)
    return jcall(obj, "getSumOfSquares", jdouble, ())
end

function get_summary(obj::StreamingStatistics)
    return jcall(obj, "getSummary", StatisticalSummary, ())
end

function get_variance(obj::StreamingStatistics)
    return jcall(obj, "getVariance", jdouble, ())
end

function hash_code(obj::StreamingStatistics)
    return jcall(obj, "hashCode", jint, ())
end

function to_string(obj::StreamingStatistics)
    return jcall(obj, "toString", JString, ())
end


text:
# ======================================================
# Type definitions
# ======================================================

abstract type AbstractValEdge{V <: Integer, E_VALS <: AbstractTuple} <: AbstractEdge{V} end

"""
    ValEdge{V, E_VALS} <: AbstractValEdge{V}

A data structure representing an undirected edge with multiple values.
"""
struct ValEdge{V<:Integer, E_VALS <: AbstractTuple} <: AbstractValEdge{V, E_VALS}
    src::V
    dst::V
    values::E_VALS

    function ValEdge{V, E_VALS}(src, dst, values) where {V, E_VALS}

        src, dst = minmax(src, dst) # TODO maybe use a branchless operator
        return new{V, E_VALS}(src, dst, values)
    end
end

"""
    ValDiEdge{V, E_VALS} <: AbstractValEdge{V}

A data structure representing a directed edge with multiple values.
"""
struct ValDiEdge{V<:Integer, E_VALS <: AbstractTuple} <: AbstractValEdge{V, E_VALS}
    src::V
    dst::V
    values::E_VALS

    function ValDiEdge{V, E_VALS}(src, dst, values) where {V, E_VALS}

        return new{V, E_VALS}(src, dst, values)
    end
end

# =========================================================
# Constructors
# =========================================================

"""
    ValEdge(s, d[, values])
    ValEdge{V}(s, d[, values])
    ValEdge{V, E_VALS}(s, d, values)

Create a `ValEdge` with source `s`, destination `d` and values `values`.

# Examples
```jldoctest
julia> e = ValEdge(1, 2, ('A',))
ValEdge 1 -- 2 with value A

julia> e = ValEdge(1, 2, ())
ValEdge 1 -- 2

julia> e = ValEdge(1,2, (first = 1, second = "2"))
ValEdge 1 -- 2 with values (first = 1, second = "2")
```
"""
function ValEdge(src::Integer, dst::Integer, values::E_VALS=()) where {E_VALS}

    src, dst = promote(src, dst)
    return ValEdge{typeof(src), E_VALS}(src, dst, values)
end

ValEdge{V}(s, d, values=()) where {V} = ValEdge{V, typeof(values)}(s, d, values)

"""
    ValDiEdge(s, d[, values])
    ValDiEdge{V}(s, d[, values])
    ValDiEdge{V, E_VALS}(s, d, values)

Create a `ValDiEdge` with source `s`, destination `d` and values `values`.

# Examples
```jldoctest
julia> e = ValDiEdge(1, 2, ('A',))
ValDiEdge 1 -> 2 with value A

julia> e = ValDiEdge(1, 2, ())
ValDiEdge 1 -> 2

julia> e = ValDiEdge(1,2, (first = 1, second = "2"))
ValDiEdge 1 -> 2 with values (first = 1, second = "2")
```
"""
function ValDiEdge(src::Integer, dst::Integer, values::E_VALS=()) where {E_VALS}

    src, dst = promote(src, dst)
    return ValDiEdge{typeof(src), E_VALS}(src, dst, values)
end

ValDiEdge{V}(s, d, values=()) where {V} = ValDiEdge{V, typeof(values)}(s, d, values)

# =========================================================
# Interface
# =========================================================

#  ------------------------------------------------------
#  src & dst
#  ------------------------------------------------------

Graphs.src(e::AbstractValEdge) = e.src
Graphs.dst(e::AbstractValEdge) = e.dst

#  ------------------------------------------------------
#  get_edgeval
#  ------------------------------------------------------

"""
    get_edgeval(e::AbstractValEdge, :)

Return the values attached to the edge `e`.

# Examples

```jldoctest
julia> e = ValEdge(1, 2, ("xyz", 123));

julia> get_edgeval(e, :)
("xyz", 123)

julia> e = ValDiEdge(1, 2, (weight=2.5,));

julia> get_edgeval(e, :)
(weight = 2.5,)
```
"""
get_edgeval(e::AbstractValEdge, ::Colon) = e.values

"""
    get_edgeval(e::AbstractValEdge, key)

Return the value attached to the edge `e` for the key `key`.

# Examples
```jldoctest
julia> e = ValEdge(1, 2, (a=11.0, ));

julia> get_edgeval(e, 1) # use integer key
11.0

julia> get_edgeval(e, :a) # use symbolic key
11.0

```
"""
get_edgeval(e::AbstractValEdge, key) = e.values[key]

#  ------------------------------------------------------
#  reverse
#  ------------------------------------------------------

Graphs.reverse(e::ValEdge) = e
Graphs.reverse(e::ValDiEdge) = ValDiEdge(dst(e), src(e), get_edgeval(e, :))

#  ------------------------------------------------------
#  is_directed
#  ------------------------------------------------------

Graphs.is_directed(::Type{<:ValEdge}) = false
Graphs.is_directed(::Type{<:ValDiEdge}) = true
Graphs.is_directed(e::AbstractValEdge) = is_directed(typeof(e))

#  ------------------------------------------------------
#  ==
#  ------------------------------------------------------

==(lhs::ValEdge, rhs::ValEdge) =
    lhs.src == rhs.src && lhs.dst == rhs.dst && lhs.values == rhs.values

==(lhs::ValDiEdge, rhs::ValDiEdge) =
    lhs.src == rhs.src && lhs.dst == rhs.dst && lhs.values == rhs.values

hash(e::ValEdge, h::UInt) = hash(e.src, hash(e.dst, hash(e.values, h)))
hash(e::ValDiEdge, h::UInt) = hash(e.src, hash(e.dst, hash(e.values, h)))

# ======================================================
# show
# ======================================================

function Base.show(io::IO, e::AbstractValEdge)
    isdir = is_directed(e)
    e_keys = keys(get_edgeval(e, :))
    has_symbol_keys = eltype(e_keys) === Symbol

    print(io, isdir ? "ValDiEdge" : "ValEdge")
    arrow = isdir ? "->" : "--"
    print(io, " $(src(e)) $arrow $(dst(e))")

    if length(e_keys) == 1
        print(io, " with value " * (has_symbol_keys ? "$(e_keys[1]) = $(get_edgeval(e, 1))" :
                                                      "$(get_edgeval(e, 1))"))
    elseif length(e_keys) > 1
        print(io, " with values $(get_edgeval(e, :))")
    end

    println(io)
end


text:
# Unit tests for Distances

function test_metricity(dist, x, y, z)
    @testset "Test metricity of $(typeof(dist))" begin
        dxy = evaluate(dist, x, y)
        dxz = evaluate(dist, x, z)
        dyz = evaluate(dist, y, z)
        if isa(dist, PreMetric)
            # Unfortunately small non-zero numbers (~10^-16) are appearing
            # in our tests due to accumulating floating point rounding errors.
            # We either need to allow small errors in our tests or change the
            # way we do accumulations...
            @test evaluate(dist, x, x) + one(eltype(x)) ≈ one(eltype(x))
            @test evaluate(dist, y, y) + one(eltype(y)) ≈ one(eltype(y))
            @test evaluate(dist, z, z) + one(eltype(z)) ≈ one(eltype(z))
            @test dxy ≥ zero(eltype(x))
            @test dxz ≥ zero(eltype(x))
            @test dyz ≥ zero(eltype(x))
        end
        if isa(dist, SemiMetric)
            @test dxy ≈ evaluate(dist, y, x)
            @test dxz ≈ evaluate(dist, z, x)
            @test dyz ≈ evaluate(dist, y, z)
        else # Not symmetric, so more PreMetric tests
            @test evaluate(dist, y, x) ≥ zero(eltype(x))
            @test evaluate(dist, z, x) ≥ zero(eltype(x))
            @test evaluate(dist, z, y) ≥ zero(eltype(x))
        end
        if isa(dist, Metric)
            # Again we have small rounding errors in accumulations
            @test dxz ≤ dxy + dyz || dxz ≈ dxy + dyz
            dyx = evaluate(dist, y, x)
            @test dyz ≤ dyx + dxz || dyz ≈ dyx + dxz
            dzy = evaluate(dist, z, y)
            @test dxy ≤ dxz + dzy || dxy ≈ dxz + dzy
        end
    end
end

@testset "PreMetric, SemiMetric, Metric on $T" for T in (Float64, F64)
    Random.seed!(123)
    n = 100
    x = rand(T, n)
    y = rand(T, n)
    z = rand(T, n)

    test_metricity(SqEuclidean(), x, y, z)
    test_metricity(Euclidean(), x, y, z)
    test_metricity(Cityblock(), x, y, z)
    test_metricity(TotalVariation(), x, y, z)
    test_metricity(Chebyshev(), x, y, z)
    test_metricity(Minkowski(2.5), x, y, z)

    test_metricity(CosineDist(), x, y, z)
    test_metricity(CorrDist(), x, y, z)

    test_metricity(ChiSqDist(), x, y, z)

    test_metricity(Jaccard(), x, y, z)
    test_metricity(SpanNormDist(), x, y, z)

    test_metricity(BhattacharyyaDist(), x, y, z)
    test_metricity(HellingerDist(), x, y, z)
    test_metricity(Bregman(x -> sqeuclidean(x, zero(x)), x -> 2*x), x, y, z);

    x₁ = rand(T, 2)
    x₂ = rand(T, 2)
    x₃ = rand(T, 2)

    test_metricity(Haversine(6371.0), x₁, x₂, x₃)

    k = rand(1:3, n)
    l = rand(1:3, n)
    m = rand(1:3, n)

    test_metricity(Hamming(), k, l, m)

    a = rand(Bool, n)
    b = rand(Bool, n)
    c = rand(Bool, n)

    test_metricity(RogersTanimoto(), a, b, c)
    test_metricity(BrayCurtis(), a, b, c)
    test_metricity(Jaccard(), a, b, c)

    w = rand(T, n)

    test_metricity(PeriodicEuclidean(w), x, y, z)
    test_metricity(WeightedSqEuclidean(w), x, y, z)
    test_metricity(WeightedEuclidean(w), x, y, z)
    test_metricity(WeightedCityblock(w), x, y, z)
    test_metricity(WeightedMinkowski(w, 2.5), x, y, z)
    test_metricity(WeightedHamming(w), a, b, c)

    Q = rand(T, n, n)
    Q = Q * Q'  # make sure Q is positive-definite

    test_metricity(SqMahalanobis(Q), x, y, z)
    test_metricity(Mahalanobis(Q), x, y, z)

    p = rand(T, n)
    q = rand(T, n)
    r = rand(T, n)
    p[p .< median(p)] .= 0
    p /= sum(p)
    q /= sum(q)
    r /= sum(r)

    test_metricity(KLDivergence(), p, q, r)
    test_metricity(RenyiDivergence(0.0), p, q, r)
    test_metricity(RenyiDivergence(1.0), p, q, r)
    test_metricity(RenyiDivergence(Inf), p, q, r)
    test_metricity(RenyiDivergence(0.5), p, q, r)
    test_metricity(RenyiDivergence(2), p, q, r)
    test_metricity(RenyiDivergence(10), p, q, r)
    test_metricity(JSDivergence(), p, q, r)
end

@testset "individual metrics" begin
    a = 1
    b = 2
    @test sqeuclidean(a, b) == 1.0

    @test euclidean(a, b) == 1.0
    @test cityblock(a, b) == 1.0
    @test totalvariation(a, b) == 0.5
    @test chebyshev(a, b) == 1.0
    @test minkowski(a, b, 2) == 1.0
    @test hamming(a, b) == 1
    @test peuclidean(a, b, 0.5) == 0
    @test peuclidean(a, b, 2) == 1.0

    bt = [true, false, true]
    bf = [false, true, true]
    @test rogerstanimoto(bt, bf) == 4.0 / 5.0
    @test braycurtis(bt, bf) == 0.5

    for T in (Float64, F64)

        for (_x, _y) in (([4.0, 5.0, 6.0, 7.0], [3.0, 9.0, 8.0, 1.0]),
                         ([4.0, 5.0, 6.0, 7.0], [3. 8.; 9. 1.0]))
            x, y = T.(_x), T.(_y)
            @test sqeuclidean(x, y) == 57.0
            @test euclidean(x, y) == sqrt(57.0)
            @test jaccard(x, y) == 13.0 / 28
            @test cityblock(x, y) == 13.0
            @test totalvariation(x, y) == 6.5
            @test chebyshev(x, y) == 6.0
            @test braycurtis(x, y) == 1.0 - (30.0 / 43.0)
            @test minkowski(x, y, 2) == sqrt(57.0)
            @test peuclidean(x, y, fill(10.0, 4)) == sqrt(37)
            @test peuclidean(x - vec(y), zero(y), fill(10.0, 4)) == peuclidean(x, y, fill(10.0, 4))
            @test peuclidean(x, y, [10.0, 10.0, 10.0, Inf]) == sqrt(57)
            @test_throws DimensionMismatch cosine_dist(1.0:2, 1.0:3)
            @test cosine_dist(x, y) ≈ (1.0 - 112. / sqrt(19530.0))
            x_int, y_int = Int64.(x), Int64.(y)
            @test cosine_dist(x_int, y_int) == (1.0 - 112.0 / sqrt(19530.0))
            @test corr_dist(x, y) ≈ cosine_dist(x .- mean(x), vec(y) .- mean(y))
            @test chisq_dist(x, y) == sum((x - vec(y)).^2 ./ (x + vec(y)))
            @test spannorm_dist(x, y) == maximum(x - vec(y)) - minimum(x - vec(y))

            @test gkl_divergence(x, y) ≈ sum(i -> x[i] * log(x[i] / y[i]) - x[i] + y[i], 1:length(x))

            @test meanad(x, y) ≈ mean(Float64[abs(x[i] - y[i]) for i in 1:length(x)])
            @test msd(x, y) ≈ mean(Float64[abs2(x[i] - y[i]) for i in 1:length(x)])
            @test rmsd(x, y) ≈ sqrt(msd(x, y))
            @test nrmsd(x, y) ≈ sqrt(msd(x, y)) / (maximum(x) - minimum(x))

            w = ones(4)
            @test sqeuclidean(x, y) ≈ wsqeuclidean(x, y, w)

            w = rand(Float64, size(x))
            @test wsqeuclidean(x, y, w) ≈ dot((x - vec(y)).^2, w)
            @test weuclidean(x, y, w) == sqrt(wsqeuclidean(x, y, w))
            @test wcityblock(x, y, w) ≈ dot(abs.(x - vec(y)), w)
            @test wminkowski(x, y, w, 2) ≈ weuclidean(x, y, w)
        end

        # Test ChiSq doesn't give NaN at zero
        @test chisq_dist([0.0], [0.0]) == 0.0

        # Test weighted Hamming distances with even weights
        a = T.([1.0, 2.0, 1.0, 3.0, 2.0, 1.0])
        b = T.([1.0, 3.0, 0.0, 2.0, 2.0, 0.0])
        w = rand(T, size(a))

        @test whamming(a, a, w) === T(0.0)
        @test whamming(a, b, w) === sum((a .!= b) .* w)

        # Minimal test of Jaccard - test return type stability.
        @inferred evaluate(Jaccard(), rand(T, 3), rand(T, 3))
        @inferred evaluate(Jaccard(), [1, 2, 3], [1, 2, 3])
        @inferred evaluate(Jaccard(), [true, false, true], [false, true, true])

        # Test Bray-Curtis. Should be 1 if no elements are shared, 0 if all are the same
        @test braycurtis([1,0,3],[0,1,0]) == 1.0
        @test braycurtis(rand(10), zeros(10)) == 1.0
        @test braycurtis([1,0],[1,0]) == 0.0

        # Test KL, Renyi and JS divergences
        r = rand(T, 12)
        p = copy(r)
        p[p .< median(p)] .= 0.0
        scale = sum(p) / sum(r)
        r /= sum(r)
        p /= sum(p)
        q = rand(T, 12)
        q /= sum(q)

        klv = 0.0
        for i = 1:length(p)
            if p[i] > 0
                klv += p[i] * log(p[i] / q[i])
            end
        end
        @test kl_divergence(p, q) ≈ klv
        @test typeof(kl_divergence(p, q)) == T

        @test renyi_divergence(p, r, 0) ≈ -log(scale)
        @test renyi_divergence(p, r, 1) ≈ -log(scale)
        @test renyi_divergence(p, r, 10) ≈ -log(scale)
        @test renyi_divergence(p, r, rand()) ≈ -log(scale)
        @test renyi_divergence(p, r, Inf) ≈ -log(scale)
        @test isinf(renyi_divergence([0.0, 0.5, 0.5], [0.0, 1.0, 0.0], Inf))
        @test renyi_divergence([0.0, 1.0, 0.0], [0.0, 0.5, 0.5], Inf) ≈ log(2.0)
        @test renyi_divergence(p, q, 1) ≈ kl_divergence(p, q)

        pm = (p + q) / 2
        jsv = kl_divergence(p, pm) / 2 + kl_divergence(q, pm) / 2
        @test js_divergence(p, q) ≈ jsv
    end
end # testset

@testset "NaN behavior" begin
    a = [NaN, 0]; b = [0, NaN]
    @test isnan(chebyshev(a, b)) == isnan(maximum(a - b))
    a = [NaN, 0]; b = [0, 1]
    @test isnan(chebyshev(a, b)) == isnan(maximum(a - b))
    @test isnan(renyi_divergence([0.5, 0.0, 0.5], [0.5, 0.5, NaN], 2))
end #testset

@testset "empty vector" begin
    for T in (Float64, F64)
        a = T[]
        b = T[]
        @test sqeuclidean(a, b) == 0.0
        @test isa(sqeuclidean(a, b), T)
        @test euclidean(a, b) == 0.0
        @test isa(euclidean(a, b), T)
        @test cityblock(a, b) == 0.0
        @test isa(cityblock(a, b), T)
        @test totalvariation(a, b) == 0.0
        @test isa(totalvariation(a, b), T)
        @test chebyshev(a, b) == 0.0
        @test isa(chebyshev(a, b), T)
        @test braycurtis(a, b) == 0.0
        @test isa(braycurtis(a, b), T)
        @test minkowski(a, b, 2) == 0.0
        @test isa(minkowski(a, b, 2), T)
        @test hamming(a, b) == 0.0
        @test isa(hamming(a, b), Int)
        @test renyi_divergence(a, b, 1.0) == 0.0
        @test isa(renyi_divergence(a, b, 2.0), T)
        @test braycurtis(a, b) == 0.0
        @test isa(braycurtis(a, b), T)

        w = T[]
        @test isa(whamming(a, b, w), T)
        @test peuclidean(a, b, w) == 0.0
        @test isa(peuclidean(a, b, w), T)
    end
end # testset

@testset "DimensionMismatch throwing" begin
    a = [1, 0]; b = [2]
    @test_throws DimensionMismatch sqeuclidean(a, b)
    a = [1, 0]; b = [2.0] ; w = [3.0]
    @test_throws DimensionMismatch wsqeuclidean(a, b, w)
    @test_throws DimensionMismatch peuclidean(a, b, w)
    a = [1, 0]; b = [2.0, 4.0] ; w = [3.0]
    @test_throws DimensionMismatch wsqeuclidean(a, b, w)
    @test_throws DimensionMismatch peuclidean(a, b, w)
    p = [0.5, 0.5]; q = [0.3, 0.3, 0.4]
    @test_throws DimensionMismatch bhattacharyya(p, q)
    @test_throws DimensionMismatch hellinger(q, p)
    Q = rand(length(p), length(p))
    Q = Q * Q'  # make sure Q is positive-definite
    @test_throws DimensionMismatch mahalanobis(p, q, Q)
    @test_throws DimensionMismatch mahalanobis(q, q, Q)
    mat23 = [0.3 0.2 0.0; 0.1 0.0 0.4]
    mat22 = [0.3 0.2; 0.1 0.4]
    @test_throws DimensionMismatch colwise!(mat23, Euclidean(), mat23, mat23)
    @test_throws DimensionMismatch colwise!(mat23, Euclidean(), mat23, q)
    @test_throws DimensionMismatch colwise!(mat23, Euclidean(), mat23, mat22)
    @test_throws DimensionMismatch colwise!(mat23, Bregman(x -> sqeuclidean(x, zero(x)), x -> 2*x), mat23, mat22)
    @test_throws DimensionMismatch evaluate(Bregman(x -> sqeuclidean(x, zero(x)), x -> 2*x), [1, 2, 3], [1, 2])
    @test_throws DimensionMismatch evaluate(Bregman(x -> sqeuclidean(x, zero(x)), x -> [1, 2]), [1, 2, 3], [1, 2, 3])
end # testset

@testset "mahalanobis" begin
    for T in (Float64, F64)
        x, y = T.([4.0, 5.0, 6.0, 7.0]), T.([3.0, 9.0, 8.0, 1.0])
        a = T.([1.0, 2.0, 1.0, 3.0, 2.0, 1.0])
        b = T.([1.0, 3.0, 0.0, 2.0, 2.0, 0.0])

        Q = rand(T, length(x), length(x))
        Q = Q * Q'  # make sure Q is positive-definite
        @test sqmahalanobis(x, y, Q) ≈ dot(x - y, Q * (x - y))
        @test eltype(sqmahalanobis(x, y, Q)) == T
        @test mahalanobis(x, y, Q) == sqrt(sqmahalanobis(x, y, Q))
        @test eltype(mahalanobis(x, y, Q)) == T
    end
end #testset

@testset "haversine" begin
    for T in (Float64, F64)
        @test haversine([-180.,0.], [180.,0.], 1.) ≈ 0 atol=1e-10
        @test haversine([0.,-90.],  [0.,90.],  1.) ≈ π atol=1e-10
        @test haversine((-180.,0.), (180.,0.), 1.) ≈ 0 atol=1e-10
        @test haversine((0.,-90.),  (0.,90.),  1.) ≈ π atol=1e-10
        @test haversine((1.,-15.625), (-179.,15.625), 6371.) ≈ 20015. atol=1e0
        @test_throws ArgumentError haversine([0.,-90., 0.25], [0.,90.], 1.)
    end
end

@testset "bhattacharyya / hellinger" begin
    for T in (Float64, F64)
        x, y = T.([4.0, 5.0, 6.0, 7.0]), T.([3.0, 9.0, 8.0, 1.0])
        a = T.([1.0, 2.0, 1.0, 3.0, 2.0, 1.0])
        b = T.([1.0, 3.0, 0.0, 2.0, 2.0, 0.0])
        p = rand(T, 12)
        p[p .< median(p)] .= 0.0
        q = rand(T, 12)

        # Bhattacharyya and Hellinger distances are defined for discrete
        # probability distributions so to calculate the expected values
        # we need to normalize vectors.
        px = x ./ sum(x)
        py = y ./ sum(y)
        expected_bc_x_y = sum(sqrt.(px .* py))
        @test Distances.bhattacharyya_coeff(x, y) ≈ expected_bc_x_y
        @test bhattacharyya(x, y) ≈ (-log(expected_bc_x_y))
        @test hellinger(x, y) ≈ sqrt(1 - expected_bc_x_y)

        pa = a ./ sum(a)
        pb = b ./ sum(b)
        expected_bc_a_b = sum(sqrt.(pa .* pb))
        @test Distances.bhattacharyya_coeff(a, b) ≈ expected_bc_a_b
        @test bhattacharyya(a, b) ≈ (-log(expected_bc_a_b))
        @test hellinger(a, b) ≈ sqrt(1 - expected_bc_a_b)

        pp = p ./ sum(p)
        pq = q ./ sum(q)
        expected_bc_p_q = sum(sqrt.(pp .* pq))
        @test Distances.bhattacharyya_coeff(p, q) ≈ expected_bc_p_q
        @test bhattacharyya(p, q) ≈ (-log(expected_bc_p_q))
        @test hellinger(p, q) ≈ sqrt(1 - expected_bc_p_q)

        # Ensure it is semimetric
        @test bhattacharyya(x, y) ≈ bhattacharyya(y, x)
    end
end #testset

function test_colwise(dist, x, y, T)
    @testset "Colwise test for $(typeof(dist))" begin
        n = size(x, 2)
        r1 = zeros(T, n)
        r2 = zeros(T, n)
        r3 = zeros(T, n)
        for j = 1:n
            r1[j] = evaluate(dist, x[:, j], y[:, j])
            r2[j] = evaluate(dist, x[:, 1], y[:, j])
            r3[j] = evaluate(dist, x[:, j], y[:, 1])
        end
        # ≈ and all( .≈ ) seem to behave slightly differently for F64
        @test all(colwise(dist, x, y) .≈ r1)
        @test all(colwise(dist, x[:, 1], y) .≈ r2)
        @test all(colwise(dist, x, y[:, 1]) .≈ r3)
    end
end

@testset "column-wise metrics on $T" for T in (Float64, F64)
    m = 5
    n = 8
    X = rand(T, m, n)
    Y = rand(T, m, n)
    A = rand(1:3, m, n)
    B = rand(1:3, m, n)

    P = rand(T, m, n)
    Q = rand(T, m, n)
    # Make sure not to remove all of the non-zeros from any column
    for i in 1:n
        P[P[:, i] .< median(P[:, i]) / 2, i] .= 0.0
    end

    test_colwise(SqEuclidean(), X, Y, T)
    test_colwise(Euclidean(), X, Y, T)
    test_colwise(Cityblock(), X, Y, T)
    test_colwise(TotalVariation(), X, Y, T)
    test_colwise(Chebyshev(), X, Y, T)
    test_colwise(Minkowski(2.5), X, Y, T)
    test_colwise(Hamming(), A, B, T)
    test_colwise(Bregman(x -> sqeuclidean(x, zero(x)), x -> 2*x), X, Y, T);

    test_colwise(CosineDist(), X, Y, T)
    test_colwise(CorrDist(), X, Y, T)

    test_colwise(ChiSqDist(), X, Y, T)
    test_colwise(KLDivergence(), P, Q, T)
    test_colwise(RenyiDivergence(0.0), P, Q, T)
    test_colwise(RenyiDivergence(1.0), P, Q, T)
    test_colwise(RenyiDivergence(Inf), P, Q, T)
    test_colwise(RenyiDivergence(0.5), P, Q, T)
    test_colwise(RenyiDivergence(2), P, Q, T)
    test_colwise(RenyiDivergence(10), P, Q, T)
    test_colwise(JSDivergence(), P, Q, T)
    test_colwise(SpanNormDist(), X, Y, T)

    test_colwise(BhattacharyyaDist(), X, Y, T)
    test_colwise(HellingerDist(), X, Y, T)
    test_colwise(BrayCurtis(), X, Y, T)

    w = rand(T, m)

    test_colwise(WeightedSqEuclidean(w), X, Y, T)
    test_colwise(WeightedEuclidean(w), X, Y, T)
    test_colwise(WeightedCityblock(w), X, Y, T)
    test_colwise(WeightedMinkowski(w, 2.5), X, Y, T)
    test_colwise(WeightedHamming(w), A, B, T)
    test_colwise(PeriodicEuclidean(w), X, Y, T)

    Q = rand(T, m, m)
    Q = Q * Q'  # make sure Q is positive-definite

    test_colwise(SqMahalanobis(Q), X, Y, T)
    test_colwise(Mahalanobis(Q), X, Y, T)
end

function test_pairwise(dist, x, y, T)
    @testset "Pairwise test for $(typeof(dist))" begin
        nx = size(x, 2)
        ny = size(y, 2)
        rxy = zeros(T, nx, ny)
        rxx = zeros(T, nx, nx)
        for j = 1:ny, i = 1:nx
            rxy[i, j] = evaluate(dist, x[:, i], y[:, j])
        end
        for j = 1:nx, i = 1:nx
            rxx[i, j] = evaluate(dist, x[:, i], x[:, j])
        end
        # As earlier, we have small rounding errors in accumulations
        @test pairwise(dist, x, y) ≈ rxy
        @test pairwise(dist, x) ≈ rxx
        @test pairwise(dist, x, y, dims=2) ≈ rxy
        @test pairwise(dist, x, dims=2) ≈ rxx
        @test pairwise(dist, permutedims(x), permutedims(y), dims=1) ≈ rxy
        @test pairwise(dist, permutedims(x), dims=1) ≈ rxx
    end
end

@testset "pairwise metrics on $T" for T in (Float64, F64)
    m = 5
    n = 8
    nx = 6
    ny = 8

    X = rand(T, m, nx)
    Y = rand(T, m, ny)
    A = rand(1:3, m, nx)
    B = rand(1:3, m, ny)

    P = rand(T, m, nx)
    Q = rand(T, m, ny)

    test_pairwise(SqEuclidean(), X, Y, T)
    test_pairwise(Euclidean(), X, Y, T)
    test_pairwise(Cityblock(), X, Y, T)
    test_pairwise(TotalVariation(), X, Y, T)
    test_pairwise(Chebyshev(), X, Y, T)
    test_pairwise(Minkowski(2.5), X, Y, T)
    test_pairwise(Hamming(), A, B, T)

    test_pairwise(CosineDist(), X, Y, T)
    test_pairwise(CorrDist(), X, Y, T)

    test_pairwise(ChiSqDist(), X, Y, T)
    test_pairwise(KLDivergence(), P, Q, T)
    test_pairwise(RenyiDivergence(0.0), P, Q, T)
    test_pairwise(RenyiDivergence(1.0), P, Q, T)
    test_pairwise(RenyiDivergence(Inf), P, Q, T)
    test_pairwise(RenyiDivergence(0.5), P, Q, T)
    test_pairwise(RenyiDivergence(2), P, Q, T)
    test_pairwise(JSDivergence(), P, Q, T)

    test_pairwise(BhattacharyyaDist(), X, Y, T)
    test_pairwise(HellingerDist(), X, Y, T)
    test_pairwise(BrayCurtis(), X, Y, T)
    test_pairwise(Bregman(x -> sqeuclidean(x, zero(x)), x -> 2*x), X, Y, T)

    w = rand(m)

    test_pairwise(WeightedSqEuclidean(w), X, Y, T)
    test_pairwise(WeightedEuclidean(w), X, Y, T)
    test_pairwise(WeightedCityblock(w), X, Y, T)
    test_pairwise(WeightedMinkowski(w, 2.5), X, Y, T)
    test_pairwise(WeightedHamming(w), A, B, T)
    test_pairwise(PeriodicEuclidean(w), X, Y, T)

    Q = rand(m, m)
    Q = Q * Q'  # make sure Q is positive-definite

    test_pairwise(SqMahalanobis(Q), X, Y, T)
    test_pairwise(Mahalanobis(Q), X, Y, T)
end

@testset "Euclidean precision" begin
    X = [0.1 0.2; 0.3 0.4; -0.1 -0.1]
    pd = pairwise(Euclidean(1e-12), X, X)
    @test pd[1, 1] == 0
    @test pd[2, 2] == 0
    pd = pairwise(Euclidean(1e-12), X)
    @test pd[1, 1] == 0
    @test pd[2, 2] == 0
    pd = pairwise(SqEuclidean(1e-12), X, X)
    @test pd[1, 1] == 0
    @test pd[2, 2] == 0
    pd = pairwise(SqEuclidean(1e-12), X)
    @test pd[1, 1] == 0
    @test pd[2, 2] == 0
end

@testset "Bregman Divergence" begin
    # Some basic tests.
    @test_throws ArgumentError bregman(x -> x, x -> 2*x, [1, 2, 3], [1, 2, 3])
    # Test if Bregman() correctly implements the gkl divergence between two random vectors.
    F(p) = LinearAlgebra.dot(p, log.(p));
    ∇(p) = map(x -> log(x) + 1, p)
    testDist = Bregman(F, ∇)
    p = rand(4)
    q = rand(4)
    p = p/sum(p);
    q = q/sum(q);
    @test evaluate(testDist, p, q) ≈ gkl_divergence(p, q)
    # Test if Bregman() correctly implements the squared euclidean dist. between them.
    @test bregman(x -> norm(x)^2, x -> 2*x, p, q) ≈ sqeuclidean(p, q)
    # Test if Bregman() correctly implements the IS distance.
    F(p) = -1 * sum(log.(p))
    ∇(p) = map(x -> -1 * x^(-1), p)
    function ISdist(p::AbstractVector, q::AbstractVector)
        return sum([p[i]/q[i] - log(p[i]/q[i]) - 1 for i in 1:length(p)])
    end
    @test bregman(F, ∇, p, q) ≈ ISdist(p, q)
end

#=
@testset "zero allocation colwise!" begin
    d = Euclidean()
    a = rand(2, 41)
    b = rand(2, 41)
    z = zeros(41)
    colwise!(z, d, a, b)
    # This fails when bounds checking is enforced
    bounds = Base.JLOptions().check_bounds
    if bounds == 0
        @test (@allocated colwise!(z, d, a, b)) == 0
    else
        @test_broken (@allocated colwise!(z, d, a, b)) == 0
    end
end
=#


text:
abstract Shape

immutable SequenceShape <: Shape
    children :: Array{Shape, 1}
end

function ==(a :: SequenceShape, b :: SequenceShape)
    return a.children == b.children
end

import Base.hash

const SequenceShape_seed = hash("SequenceShape")
function hash(a :: SequenceShape, h :: UInt)
    return hash(a.children, h + SequenceShape_seed)
end

function is_sequence(obj)
    if isa(obj, Tuple)
        return true
    end

    if !isa(obj, AbstractArray) || ndims(obj) != 1
        return false
    end

    if eltype(obj) <: Number
        return false
    end

    return true
end

immutable TensorShape <: Shape
    eltype :: Type
    size :: Array{Int64, 1}
end

function ==(a :: TensorShape, b :: TensorShape)
    return a.eltype == b.eltype && a.size == b.size
end

const TensorShape_seed = hash("TensorShape")
function hash(a :: TensorShape, h :: UInt)
    return hash(a.eltype, hash(a.size, h + TensorShape_seed))
end

function TensorShape(t :: Type, sizes :: Tuple)
    return TensorShape(t, collect(sizes))
end

function is_tensor(obj)
    return isa(obj, AbstractArray) && eltype(obj) <: Number
end

immutable TypeShape{T} <: Shape
end

function shape(obj)
    if is_sequence(obj)
        return SequenceShape(collect(map(shape, obj)))
    end
    if is_tensor(obj)
        return TensorShape(eltype(obj), collect(size(obj)))
    end
    return TypeShape{typeof(obj)}()
end

export shape


text:
# build the lookup
code = [0, 1, 2, 3, 4]

using Base.Iterators

code4 = product(code, code, code, code)

using DataFrames

function makeu(code4)
	c = [code4...]
	d = Dict{Int, Int}()
	upto = 1
	for i in 1:4
		if c[i] == 0
			## do nothing
		elseif haskey(d, c[i])
			c[i] = d[c[i]]
		else
			d[c[i]] = upto
			c[i] = upto
			upto += 1
		end
	end
	c
end

ucode4 = collect([makeu(c) for c in code4]) |> unique

res = [move!([ucode4...]) for ucode4 in ucode4]

a = mapreduce(x->transpose([x...]), vcat, ucode4)
b = mapreduce(x->transpose([x...]), vcat, res)

using CSV, DataFrames
df = DataFrame()
for (i, c) in enumerate(eachcol(hcat(a,b)))
	df[!, Symbol("ok"*string(i))] = c
end
df
CSV.write("d:/data/ok.csv", df)


text:
"""
    sv = ds_therm(v, M)

`sv` is an `M` by length(`v`) matrix of ±1s where sum(`sv`) = `v`
and the -1s are located in the bottom rows
"""
function ds_therm(v, M)
    sv = -ones(M, length(v))
    for i = 1:length(v)
        sv[1:((M+v[i])÷2), i] .= 1
    end
    return sv
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
export StokesVector, CoherencyMatrix, evpa, m̆, SingleStokes

"""
    $(TYPEDEF)
Static vector that holds the stokes parameters of a polarized
complex visibility

To convert between a `StokesVector` and `CoherencyMatrix` use the `convert`
function

```julia
convert(::CoherencyMatrix, StokesVector(1.0, 0.1, 0.1, 0.4))
```
"""
struct StokesVector{T} <: FieldVector{4,T}
    I::T
    Q::T
    U::T
    V::T
end

"""
    $(TYPEDEF)
Static matrix that holds construct the coherency matrix of a polarized
complex visibility

To convert between a `StokesVector` and `CoherencyMatrix` use the `convert`
function

```julia
convert(::StokesVector, CoherencyMatrix(1.0, 0.1, 0.1, 0.4))
```
"""
struct CoherencyMatrix{T} <: FieldMatrix{2,2,T}
    rr::T
    lr::T
    rl::T
    ll::T
end

@inline function Base.convert(::Type{CoherencyMatrix}, p::StokesVector)
    rr = p.I + p.V
    ll = p.I - p.V
    rl = p.Q + 1im*p.U
    lr = p.Q - 1im*p.U
    return CoherencyMatrix(rr, lr, rl, ll)
end

@inline function Base.convert(::Type{StokesVector}, p::CoherencyMatrix)
    i = (p.rr + p.ll)/2
    v = (p.rr - p.ll)/2
    q = (p.rl + p.lr)/2
    u = (p.rl - p.lr)/(2im)
    return StokesVector(i, q, u, v)
end

"""
    $(SIGNATURES)
Compute the fractional linear polarization of a stokes vector
or coherency matrix
"""
m̆(m::StokesVector) = (m.Q + 1im*m.U)/(m.I + eps())
m̆(m::CoherencyMatrix) = 2*m.rl/(m.rr+m.ll)

"""
    $(SIGNATURES)
Compute the evpa of a stokes vector or cohereny matrix.
"""
evpa(m::StokesVector) = 1/2*atan(m.U,m.Q)
evpa(m::StokesVector{<:Complex}) = 1/2*angle(m.U/m.Q)
evpa(m::CoherencyMatrix) = evpa(convert(StokesVector, m))

"""
    $(TYPEDEF)
Helper function that converts a model from something that compute polarized images
to just a single stokes parameter. This is useful if you just want to fit a single
stokes parameter.
"""
struct SingleStokes{M, S} <: ComradeBase.AbstractModel
    model::M
end
SingleStokes(m::M, param::Symbol) where {M} = SingleStokes{M, param}(m)

visanalytic(::Type{SingleStokes{M,S}}) where {M,S} = visanalytic(M)
imanalytic(::Type{SingleStokes{M,S}})  where {M,S} = imanalytic(M)
isprimitive(::Type{SingleStokes{M,S}}) where {M,S} = isprimitive(M)
@inline intensity_point(s::SingleStokes{M,S}, x,y) where {M,S} = getproperty(intensity_point(s.model, x,y), S)

# struct PolarizedMap{SI<:AbstractIntensityMap,
#                     SQ<:AbstractIntensityMap,
#                     SU<:AbstractIntensityMap,
#                     SV<:AbstractIntensityMap} <: AbstractPolarizedMap{SI,SQ,SU,SV}
#     I::SI
#     Q::SQ
#     U::SU
#     V::SV
#     function PolarizedMap(I::SI,Q::SQ,U::SU,V::SV) where {SI, SQ, SU, SV}
#         @assert size(I) == size(Q) == size(U) == size(V) "Image sizes must be equal in polarized map"
#         @assert fov(I) == fov(Q) == fov(U) == fov(V) "Image fov must be equal in polarized map"
#         new{SI,SQ,SU,SV}(I,Q,U,V)
#     end
# end

# Base.Base.@propagate_inbounds function Base.getindex(pimg::PolarizedMap, i...)
# return StokesVector(pimg.I[i...], pimg.Q[i...], pimg.U[i...], pimg.V[i...])
# end

# @inline stokes_parameter(pimg::PolarizedMap, p::Symbol) = getproperty(pimg, p)


text:
#Essa função de eliminação retorna U, depois de fazer a decomposição LU.
#Função faz pivotamento.

############
#Observação#
############

# Para usar essa função é necessario mandar uma matriz AA que seja do tipo float
# tentei fazer a conversão de tipo porem os resultados saiam como se fossem um número inteiro.

function eliminacao01(A)
	A
	epsilon = 0.00001;
	(m, n) = size(A);

	for i = 1:n-1
		pivot = A[i,i];

		if abs(pivot) < epsilon

			for k = i+1:m
				if abs(A[k,i] > epsilon)

					pivot = A[k,i];
					auxlinha = copy(A[k,:]);
					A[k,:] = copy(A[i,:]);
					A[i,:] = auxlinha;
					break

				end
			end
		end

		for j = i+1:m

			mult = A[j,i]/pivot;
			A[j,:] = A[j,:] - mult*A[i,:];

		end

	end

	return A;

end


text:
function _sites(t::ITensors.MPOTerm)
  return Tuple((only(o.site) for o in t.ops))
end

function opsum_sum(os::OpSum, nparts::Integer; in_partition=default_in_partition)
  oss = [OpSum() for _ in 1:nparts]
  for p in 1:nparts
    for n in 1:length(os)
      sites = _sites(os[n])
      if in_partition(sites, p, nparts)
        #oss[p] += os[n]
        push!(oss[p], os[n])
      end
    end
  end
  return oss
end

max_stars_repo_name: mtfishman/ITensorParallel

max_issues_repo_name: mtfishman/ITensorParallel

max_forks_repo_name: mtfishman/ITensorParallel


text:
using Base: AbstractFloat
import MLJModelInterface
import Tables
const MMI = MLJModelInterface

mutable struct RBFInterpolator{F <: AbstractFloat} <: MMI.Deterministic
	kernel_name :: Symbol
	kernel_args :: Vector{F}
	polynomial_degree :: Int
end

_precision( :: RBFInterpolator{F} ) where F = F

function MMI.clean!( m :: RBFInterpolator ) :: String
	warning = ""
	if m.kernel_name ∉ [:gaussian, :multiquadric, :inv_multiquadric, :cubic, :thin_plate_spline]
		warning *= "Parameter value $(m.kernel_name) for `kernel_name` not supported, using `:gaussian`.\n"
		m.kernel_name = :gaussian
	end

	φ_constructor = SymbolToRadialConstructor[ m.kernel_name ]
	try
		φ_constructor( m.kernel_args... )
	catch MethodError
		warning *= "Invalid kernel arguments. Using defaults for $(m.kernel_name).\n"
		m.kernel_args = collect(getfield.(φ_constructor(), fn for fn in fieldnames(φ_constructor)))
	end

	φ = φ_constructor( m.kernel_args... )
	if m.polynomial_degree < cpd_order(φ) - 1
		warning *= "Need a polynomial of degree at least $(cpd_order(φ) - 1).\n"
		m.polynomial_degree = cpd_order(φ) - 1
	end

	warning
end

function RBFInterpolator(; kernel_name :: Symbol = :gaussian, kernel_args :: Vector{<:Real} = [1,], polynomial_degree = 1)
	if !(eltype(kernel_args) <: AbstractFloat)
		kernel_args = Float16.(kernel_args)
	end
	model = RBFInterpolator(kernel_name, kernel_args,polynomial_degree)
	message = MMI.clean!( model )
	isempty(message) || @warn message
	return model
end

function MMI.fit( mod :: RBFInterpolator, verbosity, X :: AbstractVector{<:AbstractVector}, y :: AbstractVector )
	φ_constructor = SymbolToRadialConstructor[ mod.kernel_name ]
	φ = φ_constructor( mod.kernel_args... )

	inner_model = RBFInterpolationModel( X, y, φ , mod.polynomial_degree)
	return inner_model, nothing, nothing
end

function MMI.reformat( ::RBFInterpolator, X,  y :: AbstractVector )
	features = [ collect(Real,v)  for v in collect( Tables.rows(X) ) ]
	labels = Real.(y)
	return (features, labels)
end

function MMI.reformat( ::RBFInterpolator, X,  y )
	features = [ collect(Real,v)  for v in collect( Tables.rows(X) ) ]
	labels =  [ collect(Real,v)  for v in collect( Tables.rows(y) ) ]
	return (features, labels)
end

function MMI.reformat( ::RBFInterpolator, X )
	features = [ collect(Real,v)  for v in collect( Tables.rows(X) ) ]
	return (features,)
end

function MMI.selectrows( ::RBFInterpolator, I, X :: AbstractVector{<:AbstractVector} )
	return (X[I],)
end

function MMI.selectrows( ::RBFInterpolator, I, X :: AbstractVector{<:AbstractVector}, y :: Union{AbstractVector, Vector{<:AbstractVector}} )
	return (X[I],y[I])
end

function MMI.predict( :: RBFInterpolator, inner_model :: RBFInterpolationModel, Xnew :: AbstractVector{<:AbstractVector})
	return inner_model.(Xnew)
end

function MMI.fitted_params(model::RBFInterpolator, inner_model) ::NamedTuple
	return (
		rbf_weights = inner_model.model.rbf.weights,
		rbf_centers = [ k.c for k ∈ inner_model.model.rbf.kernels ],
		polynomial_weights = inner_model.model.psum.weights
	)
end

MMI.metadata_pkg(
	RBFInterpolator,
	name="RBFInterpolator",
        package_uuid="48790e7e-73b2-491a-afa5-62818081adcb",
        package_url="https://github.com/manuelbb-upb/RadialBasisFunctionModels.jl",
	package_license="MIT",
        is_pure_julia=true,
	is_wrapper=true,	# TODO is it a wrapper?
)

MMI.metadata_model(
    RBFInterpolator,
    input_scitype=MMI.Table(MMI.Continuous),
    target_scitype=Union{AbstractVector{MMI.Continuous}, MMI.Table(MMI.Continuous)},
    supports_weights=true,
    docstring="Radial Basis Function Interpolator",
    load_path="RadialBasisFunctionModels.RBFInterpolator"
)


text:
function examine_order(closed_order, id_str, client, config, bot_path)
    tmp_order = deepcopy(closed_order)
    order_check = Dict()
    order_check["orders"] = tmp_order

    tmp_pro    = 0
    order_check      = Dict()
    order_check["buys_exec"]  = 0
    order_check["total_vol"]  = 0
    order_check["use_vol"]    = 0
    order_check["commission"] = 0
    order_check["commission_coin"] = 0
    order_check["buy_vol"] = 0
    order_check["sell_vol"] = 0
    for (k2, v2) in tmp_order
        buy_sign = 1
        if v2["side"] == "BUY"
            order_check["buys_exec"] += 1
            buy_sign = -1
        end
        sleep(0.3)
        if !haskey(v2, "trades")
            trades = try client.get_my_trades(symbol=v2["symbol"], orderId=v2["orderId"])
                catch
                @info("Getting trades failed.")
                return
            end
            v2["trades"] = trades
            v2 = try add_fees_to_closed_order(v2, config, bot_path, client)
                catch
                @info("Adding fees failed.")
                return
            end
            closed_orders = JSON.parsefile(joinpath(bot_path, "closed_orders.json"));
            closed_orders[id_str][k2] = v2
            open(joinpath(bot_path, "closed_orders.json"), "w") do f
                JSON.print(f, closed_orders, 4)
            end
        else
            for trade in v2["trades"]
            @info("Trades already exist")
                if !haskey(trade, "com")
                    @info("No comissions stored")
                    v2 = try add_fees_to_closed_order(v2, config, bot_path, client)
                        catch
                        @info("Adding fees failed.")
                        return
                    end
                    closed_orders = JSON.parsefile(joinpath(bot_path, "closed_orders.json"));
                    closed_orders[id_str][k2] = v2
                    open(joinpath(bot_path, "closed_orders.json"), "w") do f
                        JSON.print(f, closed_orders, 4)
                    end
                end
            end
        end

        for t in v2["trades"]
            tmp_pro += buy_sign * parse(Float64, t["quoteQty"])
            tmp_pro -= t["com"]
            if buy_sign > 0
                order_check["sell_vol"] += parse(Float64, t["qty"])
            else
                order_check["buy_vol"] += parse(Float64, t["qty"])
            end
            order_check["total_vol"]  += -1 * buy_sign * parse(Float64, t["qty"])
            order_check["use_vol"]    += -1 * buy_sign * parse(Float64, t["qty"])
            if t["commissionAsset"] == config["farm_coin"]
                order_check["use_vol"] -= t["com_coin"]
            end
            order_check["commission"] += t["com"]
            order_check["commission_coin"] += t["com_coin"]
        end
        closed_orders = JSON.parsefile(joinpath(bot_path, "closed_orders.json"));
        if tmp_order != closed_orders[id_str]
            closed_orders[id_str] = tmp_order
            open(joinpath(bot_path, "closed_orders.json"), "w") do f
                JSON.print(f, closed_orders, 4)
            end
        end
    end
    order_check["current_profit"] = tmp_pro
    order_check["price_correction"]  = order_check["total_vol"] / order_check["use_vol"]
    return order_check
end
export examine_order


text:
using MPI, Test
include("../testhelpers.jl")

@testset "Ocean" begin
    runmpi(joinpath(@__DIR__, "HydrostaticBoussinesq/test_divergence_free.jl"))
    runmpi(joinpath(@__DIR__, "HydrostaticBoussinesq/test_ocean_gyre.jl"))
    # runmpi(joinpath(@__DIR__,"shallow_water/GyreDriver.jl"))
end


text:
export
    LateralDriverModel,
    ProportionalLaneTracker,
    LatLonSeparableDriver,
    Tim2DDriver,
    track_lane!

include("lateral_driving_models/lateral_driving_models.jl")
include("lane_change_models/main.jl")

include("lat_lon_separable_drivers.jl")
include("tim2d_drivers.jl")


text:
using SolidStateDetectors

include("Plotting.jl")
include("ReadGeant4Hits.jl")
include("EventSimulation.jl")

defaultDir = "plots/"

T = Float32

function BiasVariation(configFile::String, initVoltage::Real, finalVoltage::Real, length::Integer, CCDName::String="")::AbstractDict{Real, Simulation}
    biasRange = range(initVoltage, stop=finalVoltage, length=length)

    @info "Constructing base detector"
    sim_base = Simulation(configFile)
    sim_dict = Dict{Real, Simulation}()

    if CCDName != ""
        #include(CCDName)  This line needs to instead be put into the main script you are running
        sim_base.detector.semiconductors[1].charge_density_model = ccdm
    end

    @info "Calculating weighting potentials"
    for i in 1:size(sim_base.detector.contacts, 1)
        calculate_weighting_potential!(sim_base, i,
        max_refinements=2,
        convergence_limit=1e-3
        #init_grid_spacing = T.( (2e-4, 2e-4, 2e-4) )
        )
    end

    for i in 1:size(biasRange, 1)
        @info "Setting bias voltage"
        det = deepcopy(sim_base.detector)
        det.contacts[1].potential = -biasRange[i]

        sim = Simulation(det)
        calculate_electric_potential!(sim,
        depletion_handling=true,
        max_refinements=2,
        convergence_limit=1e-3
        #init_grid_spacing = T.( (2e-4, 2e-4, 2e-4) )
        )
        calculate_electric_field!(sim)

        sim.weighting_potentials = sim_base.weighting_potentials

        set_charge_drift_model!(sim, ADLChargeDriftModel())
        calculate_drift_fields!(sim)

        sim_dict[biasRange[i]] = sim
    end

    return sim_dict
end

function GetEligibleFiles(dataDir::String,
    minEnergy::Real=0., maxEnergy::Real=Inf,
    minAngle::Real=0., maxAngle::Real=Inf)::Vector{String}

    fileList = readdir(dataDir, join=true)

    filteredList = Vector{String}()

    regex = r"(\d+)+"

    for s in fileList
         m = collect(eachmatch(regex, s))
         E = parse(Int64, m[end-1].match)
         A = parse(Int64, m[end].match)
         if E >= minEnergy && E <= maxEnergy && A >= minAngle && A <= maxAngle
             push!(filteredList, s)
         end
    end

    return filteredList
end

function PlotSingleParticleDetectorResponse(sim::Simulation, fullPath::String,
    impactPosition::CartesianPoint=CartesianPoint{T}(0, 0, 0), prefix::String="",
    suffix::String="", dir::String=defaultDir, stepLimiter::Integer=1000, qe=x->1)::Vector{Event}

    fileName = split(fullPath, "/")[end]
    strippedName = split(fileName, ".")[1]
    suffix = "_" * strippedName * suffix

    gdf = GetHitInformation(fullPath)
    events = DriftGeant4Events(gdf, sim, impactPosition, stepLimiter=stepLimiter, qe=qe)
    PlotEvents(events, dir, prefix, suffix)
    return events
end

function PlotParticleDetectorResponse(sims::AbstractDict{Real, Simulation}, dataDir::String;
    impactPosition::CartesianPoint=CartesianPoint{T}(0, 0, 0), minEnergy::Real=0., maxEnergy::Real=Inf,
    minAngle::Real=0., maxAngle::Real=Inf,
    prefix::String="", dir::String=defaultDir, maxEvents::Integer=100, qe=x->1)::Nothing

    println(impactPosition)

    prefix = prefix * "particle_"
    biasRange = collect(keys(sims))
    fileList = GetEligibleFiles(dataDir, minEnergy, maxEnergy, minAngle, maxAngle)

    @showprogress for i in 1:size(biasRange, 1)
        sim = sims[biasRange[i]]
        suffix = @sprintf "_%.1fV" biasRange[i]
        @showprogress for s in fileList
            PlotSingleParticleDetectorResponse(sim, s, impactPosition, prefix, suffix, dir, maxEvents, qe)
        end
    end
end

function PlotDetectorPerformance(sims::AbstractDict{Real, Simulation}; prefix::String="", dir::String=defaultDir)::Nothing
    biasRange = collect(keys(sims))
    prefix = prefix * "detector_"
    PlotGeometry(sims[biasRange[1]], dir, prefix)
    Plot2DWeightingPotential(sims[biasRange[1]], dir, prefix)
    Plot1DWeightingPotential(sims[biasRange[1]], dir, prefix)
    PlotMaterialProperties(sims[biasRange[1]], dir, prefix)

    for i in 1:size(biasRange, 1)
        suffix = @sprintf "_%.1fV" biasRange[i]
        PlotElectricPotential(sims[biasRange[i]], dir, prefix, suffix)
        PlotElectricField(sims[biasRange[i]], dir, prefix, suffix)
        PlotPointType(sims[biasRange[i]], dir, prefix, suffix)
    end

    PlotCV(collect(values(sims)), dir, prefix)
end

max_stars_repo_name: leenderthayen/NDP

max_issues_repo_name: leenderthayen/NDP

max_forks_repo_name: leenderthayen/NDP


text:
module TestMultiReturningFunction5A

using Modia
@usingModiaPlot

mutable struct Mbs
    phi
    w
    derw
end

function setStates(phi, w)
    return Mbs(phi,w,0.0)
end

function setAccelerations(mbs::Mbs,derw)::Mbs
    mbs.derw = derw
    return mbs
end

function getForces(mbs::Mbs)
    # m*L^2*derw + m*g*l*sin(phi) ) = u
    #   m=L=g=1
    # derw + sin(phi) = 0
    u = mbs.derw + sin(mbs.phi)
    return u
end

Pendulum = Model(
    d   = 1.0,   # damping constant
    phi = Var(init=pi/2),
    w   = Var(init=0.0),
    equations = :[
        w    = der(phi)
        mbs1 = setStates(phi,w)
        mbs2 = setAccelerations(mbs1,der(w))
        u    = implicitDependency(getForces(mbs2),der(w))
        u    = -d*w
    ]
)

pendulum = @instantiateModel(Pendulum , unitless=true, log=false, logDetails=false, logCode=false, logStateSelection=true)

simulate!(pendulum, stopTime = 10.0, log=true, requiredFinalStates= [-0.012169911296941314, 0.0024087599260249094])

plot(pendulum, ("phi", "w"))

end


text:
# Check of handedness
@test CMBSim.rotationmatrix_normalized(CMBSim.qrotation_x(π/2)) * [0, 1, 0] ≈ [0, 0, 1]
@test CMBSim.rotationmatrix_normalized(CMBSim.qrotation_y(π/2)) * [1, 0, 0] ≈ [0, 0, -1]
@test CMBSim.rotationmatrix_normalized(CMBSim.qrotation_z(π/2)) * [1, 0, 0] ≈ [0, 1, 0]

# Simple composition of rotations
@test CMBSim.qrotation_x(π/3) * CMBSim.qrotation_x(π/3) ≈ CMBSim.qrotation_x(2π/3)
@test CMBSim.qrotation_y(π/3) * CMBSim.qrotation_y(π/3) ≈ CMBSim.qrotation_y(2π/3)
@test CMBSim.qrotation_z(π/3) * CMBSim.qrotation_z(π/3) ≈ CMBSim.qrotation_z(2π/3)

# More complex compositions
comp = CMBSim.qrotation_y(π/2) * CMBSim.qrotation_x(π/2)
@test CMBSim.rotationmatrix_normalized(comp) * [0, 1, 0] ≈ [1, 0, 0]

comp = CMBSim.qrotation_z(π/2) * CMBSim.qrotation_y(π/2)
@test CMBSim.rotationmatrix_normalized(comp) * [0, 0, 1] ≈ [0, 1, 0]


text:
export
    VE,
    ve

import Scruff.Operators.make_factors
import Scruff.SFuncs: Invertible, Serial
using Folds

"""
    VE(query_vars::Vector{Variable})

An instant algorithm that runs variable elimination.

# Arguments
- network
- `query_vars`: The variables to query, which are not eliminated
- depth: A depth of 1 means not to expand expanders, otherwise expands recursively to the given depth
- bounds: If true, return lower and upper bounds factors, otherwise just return a single factor
"""
struct VE <: InstantAlgorithm
    query_vars::Vector{<:Variable}
    depth::Int
    bounds::Bool
    default_range_size::Int
    VE(q; depth = 1, bounds = false, range_size = 1000) = new(q,depth,bounds,range_size)
end

function answer(::Marginal, alg::VE, runtime::Runtime, inst::VariableInstance{O}) where O
    inv = Invertible{Tuple{O},O}(x -> x[1], x -> (x,))
    jnt = joint(alg, runtime, [inst])
    if !alg.bounds
        s = Serial(Tuple{}, O, (jnt, inv))
        support(s, (), 10000, O[]) # need to precompute
        return s
    else
        (lower, upper) = jnt
        sl = Serial(Tuple{}, O, (lower, inv))
        su = Serial(Tuple{}, O, (upper, inv))
        support(sl, (), 10000, O[]) # need to precompute
        support(su, (), 10000, O[]) # need to precompute
        return (sl, su)
    end
end

function answer(query::ProbValue{O}, alg::VE, runtime::Runtime, inst::VariableInstance{O}) where O
    x = query.value
    if !alg.bounds
        m = answer(Marginal(), alg, runtime, inst)
        return cpdf(m, (), x)
    else
        (lm,um) = answer(Marginal(), alg, runtime, inst)
        l = cpdf(lm, (), x)
        u = cpdf(um, (), x)
        if isapprox(l,u)
            return(l)
        else
            error("Lower and upper bounds are not equal - use ProbabilityBounds")
        end
    end
end

function answer(query::ProbabilityBounds{O}, alg::VE, runtime::Runtime, inst::VariableInstance{O}) where O
    (lm,um) = answer(Marginal(), alg, runtime, inst)
    r = query.range
    n = length(r)
    ls = zeros(Float64, n)
    us = zeros(Float64, n)
    # Typically, when a value is not in the support of an sfunc, we interpret its probability as zero
    # However, for upper bounds, if it is in range, we want the upper bound to be 1
    usup = support(um, (), 1000, O[])
    for i in 1:n
        x = r[i]
        ls[i] = cpdf(lm, (), x)
        us[i] = x in usup ? cpdf(um, (), x) : 1.0
    end
    lsum = sum(ls)
    usum = sum(us)
    resultls = zeros(Float64, n)
    resultus = zeros(Float64, n)
    for i in 1:n
        # The probability of a value cannot be less than 1 - the upper bounds of the other values
        # or more than 1 - the lower bounds of the other values
        resultls[i] = max(ls[i], 1 - (usum - us[i]))
        resultus[i] = min(us[i], 1 - (lsum - ls[i]))
    end
    return (resultls, resultus)
end

function joint(alg::VE, runtime::Runtime, insts::Vector{<:Instance})
    vars::Vector{Variable} = [i.node for i in insts]
    if any(v -> !(v in alg.query_vars), vars)
        error("Cannot query eliminated variables")
    end
    (jnt, ids) = get_state(runtime, :joint_belief)
    if !alg.bounds
        # no bounds, just a single result
        return marginalize(runtime, normalize(jnt), ids, alg.query_vars, vars, alg.depth)
    else
        # jnt has lower and upper bounds to probabilities
        lower = marginalize(runtime, normalize(jnt[1]), ids, alg.query_vars, vars, alg.depth)
        upper = marginalize(runtime, normalize(jnt[2]), ids, alg.query_vars, vars, alg.depth)
        return (lower, upper)
    end
end

function marginalize(runtime, factor::Factor, keys::Dict{<:Node,Int}, query_vars::Vector{<:Variable}, vars_to_remain::Vector{<:Variable}, depth::Int)
    function get_var(k)
        for v in query_vars
            if keys[v] == k
                return v
            end
        end
    end
    index = [get_var(k) for k in factor.keys]
    ranges = [get_range(runtime, current_instance(runtime, v), depth) for v in index]
    selector = [v in vars_to_remain for v in index]
    combos = cartesian_product(ranges)
    result = Dict{Tuple, Float64}()
    for (combo, entry) in zip(combos, factor.entries)
        vals = []
        for i in 1:length(index)
            if selector[i]
                push!(vals, combo[i])
            end
        end
        tup = tuple(vals...)
        result[tup] = get(result, tup, 0.0) + entry
    end
    return Cat(result)
end

function infer(alg::VE, runtime::InstantRuntime,
    evidence::Dict{Symbol, Score} = Dict{Symbol, Score}(),
    interventions::Dict{Symbol, Dist} = Dict{Symbol, Dist}(),
    placeholder_beliefs = get_placeholder_beliefs(runtime))
    network = get_network(runtime)
    if !(isempty(interventions))
        error("VE cannot handle interventions")
    end
    order = topsort(get_initial_graph(network))
    ensure_all!(runtime)
    for pname in keys(placeholder_beliefs)
        ph = get_node(network, pname)
        inst = current_instance(runtime, ph)
        pi = placeholder_beliefs[pname]
        set_range!(runtime, inst, support(pi, NTuple{0,Vector}(),
            alg.default_range_size, output_type(ph)[]), alg.depth)
    end
    set_ranges!(runtime, evidence, alg.default_range_size, 1, order)
    for (n,e) in evidence
        node = get_node(network, n)
        inst = current_instance(runtime, node)
        post_evidence!(runtime, inst, e)
    end
    jnt = ve(runtime, order, alg.query_vars; depth = alg.depth, placeholder_beliefs = placeholder_beliefs, bounds = alg.bounds)
    set_state!(runtime, :joint_belief, jnt)
end

#############################################
#                                           #
# Helper functions for variable elimination #
#                                           #
#############################################

# Create a map from variable names to factor keys
function factorkeys(order) :: Dict{Node, Int}
    result = Dict{Node, Int}()
    for var in order
        result[var] = nextkey()
    end
    return result
end

# Create lower and upper bound factors for the given variable.
# Uses predetermined ranges for the variable and its parents.
#
# TODO: Change this to use an operator that returns a set of factors rather than probabilities.
# This will enable local factor decompositions for sfuncs.
function lower_and_upper_factors(runtime::Runtime, fn::Function,
                                 ids::Dict{Node, Int}, var::Variable{I,J,O}, depth) where {I,J,O}
    range = get_range(runtime, var, depth)
    pars = get_parents(get_network(runtime), var)
    parranges = []
    for v in pars
        push!(parranges, get_range(runtime, v, depth))
    end
    prtype = isempty(parranges) ? Vector{O} : typejoin([typeof(x) for x in parranges]...)
    parranges = convert(Vector{prtype}, parranges)

    inst = current_instance(runtime, var)
    sf = get_sfunc(inst)
    # @debug "lower_and_upper_factors" var=var.name sf=typeof(sf)
    parids = tuple([ids[p] for p in pars]...)
    # Since Expander needs a runtime to create factors, it has special purpose code
    # if is_fixed(var.model) && isa(make_sfunc(var, runtime), Expander)
    if isa(make_initial(var.model), Expander)
        (lowers, uppers) = expander_probs(runtime, fn, var, depth)
        keys = map(s -> ids[s], pars)
        dims = isempty(pars) ? Int[] : map(length, parranges)
        push!(keys, ids[var])
        push!(dims, length(range))
        kds = Tuple(dims)
        kts = Tuple(keys)
        return ([Factor(kds, kts, lowers)], [Factor(kds, kts, uppers)])
    else
        facts = make_factors(sf, range, Tuple(parranges), ids[var], parids)
        return facts
    end
end

function evidence_factor(runtime, var::Variable, ids::Dict)
    inst = current_instance(runtime, var)
    evidence = get_evidence(runtime, inst)
    range = get_range(runtime, inst)
    dims = Tuple(length(range))
    keys = Tuple(ids[var])
    entries = Array{Float64, 1}(undef, length(range))
    for i = 1:length(range)
        entries[i] = get_score(evidence, range[i])
    end
    return Factor(dims, keys, entries)
end

# Make all the initial factors for the given network
function produce_factors(runtime::Runtime, fn::Function,
                      order::Vector{Node}, ids::Dict, placeholder_beliefs, depth)
    lowers = []
    uppers = []
    for node in order
        if node isa Variable
            (lower, upper) =
                lower_and_upper_factors(runtime, fn, ids, node, depth)
            append!(lowers, lower)
            append!(uppers, upper)
            inst = current_instance(runtime, node)
            if has_evidence(runtime, inst)
                evfact = evidence_factor(runtime, node, ids)
                push!(lowers, evfact)
                push!(uppers, evfact)
            end
        else
            bel = placeholder_beliefs[node.name]
            range = get_range(runtime, node, depth)
            (phlower, phupper) = make_factors(bel, range, (), ids[node], ())
            append!(lowers, phlower)
            append!(uppers, phupper)
        end
    end
    return (lowers, uppers)
end

function make_graph(factors)
    result = Graph()
    for fact in factors
        ids = fact.keys
        for (i,id) in enumerate(ids)
            size = fact.dims[i]
            add_node!(result, id, size)
            for other in ids
                if other != id
                    add_undirected!(result, id, other)
                end
            end
        end
    end
    return result
end

# Eliminate the variable with the given id by multiplying all factors
# mentioning the variable and summing the variable out of the result
function eliminate(var_id, factors)
    relevant = filter(f -> var_id in f.keys, factors)
    remaining = filter(f -> !(var_id in f.keys), factors)
    if !isempty(relevant)
        prodfact = relevant[1]
        for i in 2:length(relevant)
            prodfact = product(prodfact, relevant[i])
        end
        sumfact = sum_over(prodfact, var_id)
        # If the node being eliminated is completely isolated from the rest of the network,
        # sumfact will be empty and shouldn't be added
        if !isnothing(sumfact.dims)
            push!(remaining, sumfact)
        end
    end
    return remaining
end

########################################################
#                                                      #
# Run the variable elimination algorithm               #
#                                                      #
# Works with both discrete and continuos variables,    #
# using previously determined ranges.                  #
#                                                      #
# The second argument is a topologically ordered       #
# list of variables to include in the computation.     #
# The code assumes that for any variable, its parents  #
# are present and precede it in the order.             #
# The third argument is a list of variables to query   #
# i.e. not to eliminate. This must be nonempty.        #
#                                                      #
# Returns lower and upper bound factors,               #
# as well as a key to variable names from factor keys. #
#                                                      #
########################################################

function ve(runtime::Runtime, order::Vector{<:Node},
            query_vars::Vector{<:Variable}; depth = 1, placeholder_beliefs = Dict{Symbol,Dist}(), bounds = false)
    @assert !isempty(query_vars)
    ids :: Dict{Node, Int} = factorkeys(order)
    # Making values for an expander takes a function to solve subnetworks.
    # Here we create a function that passes the bounds flag for ve.
    f(runtime, order, query_vars, depth) = ve(runtime, order, query_vars; depth = depth, placeholder_beliefs = placeholder_beliefs, bounds = bounds)
    (lowers, uppers) = produce_factors(runtime, f, order, ids, placeholder_beliefs, depth)
    ve_graph = make_graph(lowers) # Assumes that lower and upper factors have same structure
    elim_order = greedy_order(ve_graph, map(v -> ids[v], query_vars))
    for var_id in elim_order
        lowers = eliminate(var_id, lowers)
        if bounds
            uppers = eliminate(var_id, uppers)
        end
    end
    lowerprod = lowers[1]
    for i in 2:length(lowers)
        lowerprod = product(lowerprod, lowers[i])
    end
    if bounds
        # if (length(lowers[1].entries) > 45 && length(uppers[1].entries) > 45)
        #     @info("$(typeof(lowers[1])) && $(typeof(uppers[1]))",
        #         lowers=lowers[1].entries[40:1:45],
        #         uppers=uppers[1].entries[40:1:45])
        # end
        upperprod = uppers[1]
        for i in 2:length(uppers)
            upperprod = product(upperprod, uppers[i])
        end
        return ((lowerprod, upperprod), ids)
    else
        return (lowerprod, ids)
    end
end

function copy_graph(g :: Graph)
    ns = copy(g.nodes)
    es = Dict()
    ss = Dict()
    for n in ns
        es[n] = copy(g.edges[n])
        ss[n] = g.sizes[n]
    end
    return Graph(ns, es, ss)
end

##########################
#                        #
# Eliminating a variable #
#                        #
##########################

function unconnected_neighbors(g :: Graph, n :: Int)
    neighbors = g.edges[n]
    m = length(neighbors)
    result = []
    for i = 1:m
        n1 = neighbors[i]
        for j = i+1:m
            n2 = neighbors[j]
            if !(n1 in g.edges[n2]) || !(n2 in g.edges[n1])
                push!(result, (n1, n2))
            end
        end
    end
    return result
end

function eliminate(g :: Graph, n :: Int)
    ns = unconnected_neighbors(g, n)
    deleteat!(g.nodes, findfirst(m -> m == n, g.nodes))
    delete!(g.edges, n)
    delete!(g.sizes, n)
    for (m,ms) in g.edges
        if n in ms
            deleteat!(ms, findfirst(m -> m == n, ms))
        end
    end
    for (i, j) in ns
        add_undirected!(g, i, j)
    end
end

#############################################
#                                           #
# Greedily determining an elimination order #
# to minimize the number of edges added     #
#                                           #
#############################################

function cost(g :: Graph, n :: Int)
    return length(unconnected_neighbors(g, n))
end

function greedy_order(g :: Graph, to_leave :: Array{Int})
    candidates = filter(n -> !(n in to_leave), g.nodes)
    result = []
    h = copy_graph(g)
    while !isempty(candidates)
        costs = map(c -> cost(h, c), candidates)
        best = candidates[argmin(costs)]
        push!(result, best)
        eliminate(h, best)
        deleteat!(candidates, findfirst(n -> n == best, candidates))
    end
    return result
end

function greedy_order(g :: Graph)
    a :: Array{Int, 1} = []
    return greedy_order(g, a)
end


text:
module Coulter

    __precompile__()

    using KernelDensity
    using Distributions
    using Gadfly
    using StatsBase
    using DataStructures
    import Base.-, Base.deepcopy
    using Dates: Second, DateTime

    export CoulterCounterRun, loadZ2, -, volume, diameter,
        extract_peak, extract_peak!, extract_peak_interval

    include("utils.jl")

    """
    A simplified representation of a coulter counter run
    """
    mutable struct CoulterCounterRun
        filename::String
        sample::String
        timepoint::DateTime
        reltime::Union{Second, Nothing}
        binlims::Vector{Float64}
        binvols::Vector{Float64}
        binheights::Vector{Number}
        """Whether the `binheights` is the raw count or total volume"""
        yvariable::Symbol
        data::Vector{Float64}
        livepeak::Union{Float64, Nothing}
        allpeaks::Vector{Float64}
        params::Dict{String, Any}
    end

    include("analysis.jl")

    # copy constructor
    deepcopy(a::CoulterCounterRun) = CoulterCounterRun([deepcopy(getfield(a, field)) for field in fieldnames(CoulterCounterRun)]...)

    """
    loadZ2(filename::String, sample::String)

    Loads `filename` and assigns it a sample, returns a
    `CoulterCounterRun` object
    """
    function loadZ2(filepath::String, sample::String; yvariable=:count)
        open(filepath) do s
            # split windows newlines if present
            filebody = replace(read(s, String), "\r\n"=>"\n")
            # extract start time and date from body
            datetime = match(r"^StartTime= \d*\s*(?<time>\d*:\d*:\d*)\s*(?<date>\d*\s\w{3}\s\d{4})$"m, filebody)
            timepoint = DateTime("$(datetime[:date]) $(datetime[:time])", "dd uuu yyy HH:MM:SS")

            params = Dict{String, Any}()
            params["Current"] = parse(Float64, match(r"^Cur=([+-]?[0-9]*[.]?[0-9]+)$"m, filebody)[1])
            params["Pre-Amp Gain"] = parse(Float64, match(r"^PAGn=([+-]?[0-9]*[.]?[0-9]+)$"m, filebody)[1])
            params["Volume metered"] = parse(Float64, match(r"^Vol=([+-]?[0-9]*[.]?[0-9]+)$"m, filebody)[1])

            # extract data
            matcheddata = match(r"^\[#Bindiam\]\n(?<binvols>.*?)\n^\[Binunits\].*?\[#Binheight\]\n(?<binheight>.*?)\n^\[end\]"sm, filebody)
            binlims = [volume(parse(Float64, x)) for x in split(matcheddata[:binvols], "\n ")]
            # the Coulter software takes the mean the upper and lower limits to
            # compute the volume of everything that falls in the bin
            binvols = [(binlims[i] + binlims[i+1])/2 for i in 1:length(binlims)-1]
            binheights = Int[parse(Int, x) for x in split(matcheddata[:binheight], "\n ")]
            pop!(binheights) # remove the last extraneous value

            # unbin data, i.e. the inverse of the hist function
            data = repvec(binvols, binheights)

            if yvariable == :volume
                binheights = binvols .* binheights
            end

            CoulterCounterRun(basename(filepath), sample, timepoint, nothing,
                              binlims, binvols, binheights, yvariable, data,
                              nothing, Float64[], params)
        end
    end

    loadZ2(filepath::String) = loadZ2(filepath, "N/A")

    -(a::CoulterCounterRun, b::CoulterCounterRun) = a.timepoint - b.timepoint
    -(a::CoulterCounterRun, b::DateTime) = a.timepoint - b
end