text:
module TransposeNoPacking

using KernelAbstractions.Extras: @unroll
using SIMD: Vec, vload, vstore
using StaticArrays: MArray, MVector
using MyGEMM.Types: TilingParams

export multiply_transpose_no_packing!

"""
  multiply_transpose_no_packing!(C, A, B, params)

Performs `C += A * B'` without packing optimization.
"""
@inline function multiply_transpose_no_packing!(C, A, B, params::TilingParams)
  # Get the matrix params
  m, n = size(C)
  k = size(B, 2)

  # Check the matrix params
  @assert size(C) == (m, n)
  @assert size(A) == (m, k)
  @assert size(B) == (n, k)

  # Get tiling parameters
  mr = params.mr
  nr = params.nr
  vl = params.vl

  # We assume that the matrices are even divisible by `nr` and `mr`,
  # that the blocking perfectly partitions the matrix,
  # and, vectors are perfectly divided by the chosen vector unit
  @assert mod(m, mr)  == 0
  @assert mod(n, nr)  == 0
  @assert mod(mr, vl) == 0

  # Execute C = C + A * B'
  loop5!(C, A, B, params)
end

@inline function loop5!(C, A, B, params::TilingParams)
  n = size(C, 2)
  nc = params.nc
  @inbounds for j1 = 1:nc:n
    j2 = min(n, j1 + nc - 1)
    C_j = @view C[:, j1:j2]
    B_j = @view B[j1:j2, :]
    loop4!(C_j, A, B_j, params)
  end
end

@inline function loop4!(C, A, B, params ::TilingParams)
  k = size(A, 2)
  kc = params.kc
  @inbounds for p1 = 1:kc:k
    p2 = min(k, p1 + kc - 1)
    A_p = @view A[:, p1:p2]
    B_p = @view B[:, p1:p2]
    loop3!(C, A_p, B_p, params
          )
  end
end

@inline function loop3!(C, A, B, params::TilingParams)
  m = size(C, 1)
  mc = params.mc
  @inbounds for i1 = 1:mc:m
    i2 = min(m, i1 + mc - 1)
    C_i = @view C[i1:i2, :]
    A_i = @view A[i1:i2, :]
    loop2!(C_i, A_i, B, params)
  end
end

@inline function loop2!(C, A, B, params::TilingParams)
  n = size(C, 2)
  nr = params.nr
  @inbounds for j1 = 1:nr:n
    j2 = min(n, j1 + nr - 1)
    C_j = @view C[:, j1:j2]
    B_j = @view B[j1:j2, :]
    loop1!(C_j, A, B_j, params
          )
  end
end

@inline function loop1!(C, A, B, params::TilingParams)
  m = size(C, 1)
  mr = params.mr
  @inbounds for i1 = 1:mr:m
    i2 = min(m, i1 + mr - 1)
    C_i = @view C[i1:i2, :]
    A_i = @view A[i1:i2, :]
    microknl!(C_i, A_i, B, params
      )
  end
end

@inline function microknl!(C, A, B, params::TilingParams)
  # Get relevant tile params
  m     = size(parent(C), 1)
  mr    = params.mr
  nr    = params.nr
  vl    = params.vl
  mr_vl = div(mr, vl)
  pend  = size(B, 2)

  # Set relevant types
  FT = eltype(C)
  VT = Vec{vl, FT}

  # Initialize vector storage
  C_tile  = MArray{Tuple{mr_vl, nr}, VT}(undef)
  A_panel = MVector{mr_vl, VT}(undef)

  # Load the columns of the microtile of C
  @inbounds @unroll for j = 1:nr
    @unroll for i = 1:mr_vl
      offset = ((i - 1) * vl + (j - 1) * m) * sizeof(FT)
      C_tile[i, j] = vload(VT, pointer(C) + offset)
    end
  end

  # Do the rank-one updates
  @inbounds @unroll for p = 1:pend
    # Load pieces of `A_panel`
    @unroll for i = 1:mr_vl
      offset = ((p - 1) * m + (i - 1) * vl) * sizeof(FT)
      A_panel[i] = vload(VT, pointer(A) + offset)
    end

    # Do outer product updates
    @unroll for j = 1:nr
      β = B[j, p]
      @unroll for i = 1:mr_vl
        α = A_panel[i]
        C_tile[i, j] = muladd(α, β, C_tile[i, j])
      end
    end
  end

  # Write back from vector storage
  @inbounds @unroll for j = 1:nr
    @unroll for i = 1:mr_vl
      offset = ((i - 1) * vl + (j - 1) * m) * sizeof(FT)
      vstore(C_tile[i, j], pointer(C) + offset)
    end
  end
end

end

max_stars_repo_name: bischtob/MyGEMM

max_issues_repo_name: bischtob/MyGEMM

max_forks_repo_name: bischtob/MyGEMM


text:
using Sundials
function exotherm(t, x, dx; n=1, a=1, b=1)
  p = x[2]^n * exp(x[1])
  dx[1] = p - a*x[1]
  dx[2] = -b*p
  return(dx)
end;

t = linspace(0.0,5.0,1001);
fexo(t,x,dx) = exotherm(t, x, dx, a=0.6, b=0.1);

x1 = Sundials.cvode(fexo, [0.9, 1.0], t);
using Winston
plot(t,x1[:,1]))

ext: jl

max_stars_repo_path: Module 3/Chapter06/exotherm1.jl

max_issues_repo_path: Module 3/Chapter06/exotherm1.jl

max_forks_repo_path: Module 3/Chapter06/exotherm1.jl


text:
#=
Handling of time in data as a physical quantity, and time-related data processing
=#
using Statistics, StatsBase
#########################################################################
# Datetime related
#########################################################################
using Dates
export monthday_indices, maxyearspan, daymonth, DAYS_IN_YEAR, time_in_days
export temporal_sampling

const DAYS_IN_YEAR = 365.26
millisecond2month(t) = Month(round(Int, t.value / 1000 / 60 / 60 / 24 / 30))
daymonth(t) = day(t), month(t)

maxyearspan(A::AbDimArray, tsamp = temporal_sampling(A)) =
maxyearspan(dims(A, Time).val, tsamp)

"""
    maxyearspan(A::ClimArray) = maxyearspan(dims(A, Time))
    maxyearspan(t::Vector{<:DateTime}) → i
Find the maximum index `i` of `t` so that `t[1:i]` covers exact(*) multiples of years.

(*) For monthly spaced data `i` is a multiple of `12` while for daily data we find
the largest possible multiple of `DAYS_IN_YEAR = 365.26`.
"""
function maxyearspan(times, tsamp = temporal_sampling(times))
    if tsamp == :monthly
        length(times) % 12 == 0 && return length(times)
        m = month(times[1])
        findlast(i -> month(times[i]) == m, 1:length(times)) - 1
    elseif tsamp == :yearly
        return length(times)
    elseif tsamp == :daily
        m, d = monthday(times[1])
        i = findlast(j -> monthday(t[j]) == (m, d), 1:length(t))
        return i-1
    else
        error("maxyearspan: not implemented yet for $tsamp data")
    end
end

"""
    monthday_indices(times, date = times[1])
Find the indices in `times` (which is a `Vector{Date}`) at which
the date in `times` gives the same day and month as `date`.
"""
function monthday_indices(times, date = times[1])
    d1, m1 = daymonth(date)
    a = findall(i -> daymonth(times[i]) == (d1, m1), 1:length(times))
end

"""
    monthamount(x::AbstractArray{<:AbstractDateTime})
Convert `x` to the amount of months (starting from 0).
"""
function monthamount(t)
    L = length(t)
    yd = [year(t[i]) - year(t[1]) for i in 1:L]
    md = [month(t[i]) - month(t[1]) for i in 1:L]
    [yd[i]*12 + md[i] for i in 1:L]
end

"""
    monthspan(t::TimeType)
Return a vector of daily spaced `Date`s that span the entire month that `t` belongs in.
"""
function monthspan(t::TimeType)
    m = mod1(month(t), 12)
    n = mod1(m+1, 12)
    y = year(t)
    u = m == 12 ? y+1 : y
    d = collect(Date(y, m, 1):Day(1):Date(u, n, 1))[1:end-1]
end

"""
    time_in_days(t::AbstractArray{<:TimeType}, T = Float32)
Convert a given date time array into measurement units of days:
a real-valued array which counts time in days, always increasing.
"""
function time_in_days(t::AbstractArray{<:TimeType}, T = Float32)
    ts = temporal_sampling(t)
    if ts == :monthly
        truetime = daysinmonth.(t)
        r = T.(cumsum(truetime))
    elseif ts == :yearly
        error("todo")
    elseif ts == :daily
        error("todo")
    end
    return r
end
time_in_days(t::AbstractArray{<:Real}) = t

"""
    temporal_sampling(x) → symbol
Return the temporal sampling type of `x`, which is either an array of `Date`s or
a dimensional array (with `Time` dimension).

Possible return values are:
- `:yearly`, where all dates have the same month+day, but different year.
- `:monthly`, where all dates have the same day, but different month.
- `:daily`, where the temporal difference between dates are exactly 1 day.
- `:hourly`, where the temporal difference between entries is exactly 1 hour.
- `:other`, which means that `x` doesn't fall to any of the above categories.

For vector input, only the first 3 entries of the temporal information are used
to deduce the sampling (while for ranges, checking the step is enough).
"""
temporal_sampling(A::AbDimArray) = temporal_sampling(dims(A, Time).val)
function temporal_sampling(t::AbstractVector{<:TimeType})
    #TODO: implement hourly!
    sampled_less_than_date(t) && error("Hourly sampling not yet implemented")
    d1 = daymonth(t[2]) .- daymonth(t[1])
    d2 = daymonth(t[3]) .- daymonth(t[2])
    samemonth = d1[2] == d1[2] == 0
    sameday = d1[1] == d2[1] == 0
    sameyear = year(t[1]) == year(t[2]) == year(t[3])
    if sameday && samemonth && !sameyear
        :yearly
    elseif sameday && !samemonth
        :monthly
    elseif !sameday && samemonth
        :daily
    # TODO: test that <0 works, add in test suite
    elseif !sameday && !samemonth && (d1[1] < 0 || d2[1] < 0)
        :daily
    else
        :other
    end
end
temporal_sampling(t::AbstractVector) = :other
temporal_sampling(t::StepRange{<:Any,Month}) = :monthly
temporal_sampling(t::StepRange{<:Any,Year}) = :yearly
temporal_sampling(t::StepRange{<:Any,Day}) = :daily
temporal_sampling(t::StepRange{<:Any,Hour}) = :hourly
temporal_sampling(t::StepRange{<:Any,<:Any}) = :other

"return true if hours or minutes are ≠ 0."
function sampled_less_than_date(t::AbstractVector{<:DateTime})
    r = 1:length(t)
    any(i -> Dates.hour(t[i]) ≠ 0, r) || any(i -> Dates.minute(t[i]) ≠ 0, r)
end
sampled_less_than_date(t::AbstractVector{<:Date}) = false

"return the appropriate subtype of Dates.Period."
function tsamp2period(tsamp)
    tsamp == :monthly && return Month(1)
    tsamp == :yearly && return Year(1)
    tsamp == :daily && return Day(1)
    error("Don't know the period of $tsamp sampling!")
end

#########################################################################
# temporal statistics
#########################################################################
export timemean, timeagg

"""
    timemean(A::ClimArray [, w]) = timeagg(mean, A, w)
Temporal average of `A`.
"""
timemean(A::ClimArray, w = nothing) = timeagg(mean, A, w)

"""
    timeagg(f, A::ClimArray, W = nothing)
Perform a proper temporal aggregation of the function `f` (e.g. `mean, std`)
on `A` (assuming monthly spaced data) where:
* Only full year spans of `A` are included, see [`maxyearspan`](@ref)
  (because most processes are affected by yearly cycle,
  and averaging over an uneven number of cycles typically results in artifacts)
* Each month in `A` is weighted with its length in days (for monthly sampled data)

If you don't want these features, just do [`dropagg`](@ref)`(f, A, Time, W)`.

`W` are possible statistical weights that are used in conjuction to the temporal weighting,
to weight each time point differently.
If they are not a vector (a weight for each time point), then they have to be a dimensional
array of same dimensional layout as `A` (a weight for each data point).

    timeagg(f, t::Vector, x::Vector, w = nothing)
Same as above, but for arbitrary vector `x` accompanied by time vector `t`.
"""
function timeagg(f, A::AbDimArray, w = nothing)
    w isa AbDimArray && @assert dims(w) == dims(A)
    w isa Vector && @assert length(w) == size(A, Time)
    tsamp = temporal_sampling(A)
    r = if tsamp == :daily || tsamp == :other
        timeagg_daily(f, A, w)
    elseif tsamp == :monthly
        timeagg_monthly(f, A, w)
    elseif tsamp == :yearly
        timeagg_yearly(f, A, w)
    end
    n = A.name == "" ? "" : A.name*", temporally aggregated with $(string(f))"
    R = ClimArray(r, otherdims(A, Time()), n)
end

function timeagg_yearly(f, A, w)
     if isnothing(w)
        dropagg(f, A, Time)
    elseif w isa AbDimArray
        map(i -> f(view(A, i), weights(view(W, i))), otheridxs(A, Time()))
    elseif w isa Vector
        fw = weights(w)
        map(i -> f(view(A, i), fw), otheridxs(A, Time()))
    end
end

function timeagg_monthly(f, A::AbDimArray, w)
    t = dims(A, Time).val
    mys = maxyearspan(t, :monthly)
    tw = daysinmonth.(t)
    W = if isnothing(w)
        tw
    elseif w isa Vector
        tw .* w
    elseif w isa AbDimArray
        _w = dimwise(*, w, ClimArray(tw, (Time(t),)))
        @view _w[Time(1:mys)]
    end
    other = otherdims(A, Time)
    _A = @view A[Time(1:mys)]
    if w isa AbDimArray
        r = map(i -> f(view(_A, i), weights(view(W, i))), otheridxs(A, Time()))
    else
        fw = weights(view(W, 1:mys))
        r = map(i -> f(view(_A, i), fw), otheridxs(A, Time()))
    end
end

function timeagg_daily(f, A::AbDimArray, w)
    t = dims(A, Time).val
    mys = maxyearspan(t)
    _A = view(A, Time(1:mys))
    if w isa AbDimArray
        W = view(w, Time(1:mys))
        r = map(i -> f(view(_A, i), weights(view(W, i))), otheridxs(A, Time()))
    elseif w isa Vector
        fw = weights(view(w, 1:mys))
        r = map(i -> f(view(_A, i), fw), otheridxs(A, Time()))
    elseif isnothing(w)
        r = dropagg(f, _A, Time)
    end
end

function timeagg(f, a::AbDimArray{T, 1}, w = nothing) where {T} # version with only time dimension
    t = dims(a, Time).val
    return timeagg(f, t, Array(a), w)
end

function timeagg(f, T::AbstractVector{<:TimeType}, a::Vector, w = nothing) # version with just vectors
    tsamp = temporal_sampling(T)
    mys = maxyearspan(T, tsamp)
    t = view(T, 1:mys)
    if tsamp == :monthly
        dimw = float.(daysinmonth.(t))
        !isnothing(w) && (dimw .*= view(w, 1:mys))
        return f(view(a, 1:mys), weights(dimw))
    else
        if isnothing(w)
            return f(view(a, 1:mys))
        else
            return f(view(a, 1:mys), weights(view(w, 1:mys)))
        end
    end
end

#########################################################################
# Monthly/yearly/daily means
#########################################################################
export monthlyagg, yearlyagg, temporalrange

"""
    monthlyagg(A::ClimArray, f = mean) -> B
Create a new array where the temporal daily information has been aggregated to months
using the function `f`.
By convention, the dates of the new array always have day number of `15`.
"""
function monthlyagg(A::ClimArray, f = mean)
    t0 = dims(A, Time).val
    startdate = Date(year(t0[1]), month(t0[1]), 15)
    finaldate = Date(year(t0[end]), month(t0[end]), 16)
    t = startdate:Month(1):finaldate
    tranges = temporalrange(t0, Dates.month)
    return timegroup(A, f, t, tranges, "monthly")
end

"""
    yearlyagg(A::ClimArray, f = mean) -> B
Create a new array where the temporal information has been aggregated to years
using the function `f`.
By convention, the dates of the new array always have month and day number of `1`.
"""
function yearlyagg(A::ClimArray, f = mean)
    t0 = dims(A, Time).val
    startdate = Date(year(t0[1]), 1, 1)
    finaldate = Date(year(t0[end]), 2, 1)
    t = startdate:Year(1):finaldate
    tranges = temporalrange(t0, Dates.year)
    return timegroup(A, f, t, tranges, "yearly")
end

function timegroup(A, f, t, tranges, name)
    other = otherdims(A, Time)
    n = A.name == "" ? "" : A.name*", $(name) averaged"
    B = ClimArray(zeros(eltype(A), length.(other)..., length(t)), (other..., Time(t)), n)
    for i in 1:length(tranges)
        B[Time(i)] .= dropagg(f, view(A, Time(tranges[i])), Time)
    end
    return B
end

"""
    temporalrange(t::AbstractVector{<:TimeType}}, f = Dates.month) → r
Return a vector of ranges so that each range of indices are values of `t` that
belong in either the same month, year, or day, depending on `f`.
`f` can take the values: `Dates.year, Dates.month, Dates.day` (functions).

Used in e.g. [`monthlyagg`](@ref) or [`yearlyagg`](@ref).
"""
function temporalrange(t::AbstractArray{<:TimeType}, f = Dates.month)
    @assert issorted(t) "Sorted time required."
    L = length(t)
    r = Vector{UnitRange{Int}}()
    i, x = 1, f(t[1]) # previous entries
    for j in 2:L
        y = f(t[j])
        x == y && continue
        push!(r, i:(j-1))
        i, x = j, y
    end
    push!(r, i:L) # final range not included in for loop
    return r
end


text:
module ElecElecRepulsion

#LearnHatreeFock.jl Modules
using TypesBasis,TypesParticles
using UtilityFunctions

export buildelecelecrepulsion

@doc raw"""
function buildelecelecrepulsion()

description: Implement electron-electron repulsion within Hatree-Fock framework. The electron-electron repuslion is done within a mean-field approach. That is given the electron density what is the average repulsion energy an electron feels due to the charge density located at a different point. Recall that charge density is given by the overlap of two orbitals. Our operator is thus:

```
$ \int e^2 frac{\psi_{m}}{\psi_{n}} dr $
```

We need to have all possible combinations of orbitals for the electron-electron operator, N^4.

Gaussian orbitals centered at different locations. Break-down ad two set of orbitals using gaussian product theorem, so we get two gaussian at new centers called p and p-prime. The integral is then given by:

```
$g(l,m,n,o) = A_p\,A_{p^{`}} frac{2\pi^{\frac{5}{2}}}{p\,p^{`}\,\sqrt{p+p^{`}}} F_o\left(\alpha_{p,p^{`}}\,R_{p,p^{`}}^2 \right)$
```

The steps for solving the integral are:
1. Loop over all basis functions, i.e., l,m,n,o
2. Loop over all primitive Gaussian functions of l,m,n,o.
3. Apply Gaussian product theorem to p=lm and p`=no

input:
output:

"""
function buildelecelecrepulsion(numelec::Int,basisfunc::Array{GaussOrbitals},
                                basis::Basis)

    numbasisfunc = basis.nbasisfunc; #number of basis functions
    basissize = numelec*numbasisfunc;

    # Ne*Nbasis x Ne*Nbasis x Ne*Nbasis x Ne*Nbasis
    elecelecrepul = zeros(Float64,basissize,basissize,
                          basissize,basissize);

    flatbasisfunc = flattenbasisfunc(numelec,numbasisfunc,basisfunc);

    #NOTE: To try and comp. cost make min. call/ref.
    #Build electron-electron Coulomb intergral
    for l=1:basissize
        basisl = flatbasisfunc[l];
        primsl = length(basisl.coefs);
        for lp=1:primsl
            for m=1:basissize
                basism = flatbasisfunc[m];
                primsm = length(basism.coefs);
                for mp=1:primsm
                    #Build new Gaussian center from l,m Gaussian orbitals.
                    #TODO: Turn into function
                    gausslm = getnewcenter(basisl,lp,basism,mp);

                    #Now build new Gaussian center for n,o
                    for n=1:basissize
                        basisn = flatbasisfunc[n];
                        primsn = length(basisn.coefs);
                        for np=1:primsn
                            for o=1:basissize
                                basiso = flatbasisfunc[o];
                                primso = length(basiso.coefs);
                                for op=1:primso

                                    gaussno = getnewcenter(basisn,np,basiso,op);
                                    alpha = (gausslm.p*gaussno.p)/(gausslm.p+gaussno.p);
                                    dist = (gausslm.px-gaussno.px)^2 +
                                           (gausslm.py-gaussno.py)^2 +
                                           (gausslm.pz-gaussno.pz)^2;
                                    #NOTE: This calls a spherical implemented version
                                    boysfunc = boysfunction(alpha*dist,n=0);
                                    prefact = (2.00e0*pi^2.50e0)/(gausslm.p*gaussno.p*sqrt(gausslm.p+gaussno.p));
                                    elecelecrepul[l,m,n,o] += prefact*gausslm.ampli*gaussno.ampli*boysfunc;
                                end
                            end # o=1:basissize
                        end
                    end # n=1:basissize
                end
            end # m=1:basissize
        end
    end # l=1:basissize

    return elecelecrepul
end #buildelecelecrepulsion

@doc raw"""
function getnewcenter();

"""
function getnewcenter(basisn::GaussOrbitals,np::Int,
                      basism::GaussOrbitals,mp::Int)
    alpha_np = basisn.alphas[np]::Float64;
    alpha_mp = basism.alphas[mp]::Float64;
    coef_np = basisn.coefs[np]::Float64;
    coef_mp = basism.coefs[mp]::Float64;
    norm_np = basisn.norms[np]::Float64;
    norm_mp = basism.norms[mp]::Float64;
    p = alpha_np+alpha_mp;
    px = (basisn.ax*alpha_np + basism.ax*alpha_mp)/p;
    py = (basisn.ay*alpha_np + basism.ay*alpha_mp)/p;
    pz = (basisn.az*alpha_np + basism.az*alpha_mp)/p;
    spatialx_nm =gaussprod1D(basisn.ax,alpha_np,basism.ax,alpha_mp);
    spatialy_nm =gaussprod1D(basisn.ay,alpha_np,basism.ay,alpha_mp);
    spatialz_nm =gaussprod1D(basisn.az,alpha_np,basism.az,alpha_mp);
    spatial_nm = spatialx_nm*spatialy_nm*spatialz_nm;
    ampli_nm = norm_np * norm_mp * coef_np * coef_mp * spatial_nm;
    return NewGauss(p,px,py,pz,ampli_nm)

end #getnewcenter

"""
TODO: should I move this to BasisTypes.jl?

""" struct NewGauss
     p::Float64
     px::Float64
     py::Float64
     pz::Float64
     ampli::Float64
     function NewGauss(p,px,py,pz,ampli)
          new(p,px,py,pz,ampli);
      end
end #NewGauss

end #ElecElecRepulsion


text:
using Test
using LIBSVMFileIO
using SparseArrays
using Random
Random.seed!(0)

function readcorrectly(data, labels, d_read, l_read, precision=1e6)
	@test length(d_read) == length(data)
	@test length(l_read) == length(labels)
	data_correct = true
	labels_correct = true
	for i = 1:length(d_read)
		data_correct =
			all(isapprox.(d_read[i], round.(data[i].*precision)./precision)) &&
			data_correct
		labels_correct =
			all(isapprox.(l_read[i], round.(labels[i].*precision)./precision))&&
			labels_correct
	end
	@test data_correct
	@test labels_correct
end

@testset "Basics" begin
	println("Testing Basics")
	data = [20 .*randn(5) for _ = 1:10]
	labels = [rand(-10:10) for _ = 1:10]
	file = "tmp.libsvm"

	libsvmwrite(data, labels, file)
	for l in eachline(file)
		println(l)
	end

	nEx, nFeat = libsvmsize(file)
	@test nEx == 10
	@test nFeat == 5

	d_read, l_read = libsvmread(file)
	readcorrectly(data, labels, d_read, l_read)

	libsvmwrite(data, labels, file, precision=2)

	nEx, nFeat = libsvmsize(file)
	@test nEx == 10
	@test nFeat == 5

	d_read, l_read = libsvmread(file)
	readcorrectly(data, labels, d_read, l_read, 1e2)
end

@testset "Sparse Data" begin
	println("Testing Sparse Data")
	data = [20 .*sprandn(50, 0.4) for _ = 1:20]
	data[1][end] = 1.0
	labels = [rand(-10:10) for _ = 1:20]
	file = "tmp.libsvm"

	libsvmwrite(data, labels, file, precision=5)

	nEx, nFeat = libsvmsize(file)
	@test nEx == 20
	@test nFeat == 50

	d_read, l_read = libsvmread(file)
	readcorrectly(data, labels, d_read, l_read, 1e5)

	# Zeros out last data points
	nTrailing = 3
	for d in data
		d[(end-nTrailing):end] .= 0.0
	end

	libsvmwrite(data, labels, file)

	nEx, nFeat = libsvmsize(file)
	@test nEx == 20
	@test nFeat == 50-nTrailing-1

	data_ctrl = [d[1:(end-nTrailing-1)] for d in data]
	d_read, l_read = libsvmread(file)
	readcorrectly(data_ctrl, labels, d_read, l_read)
end

@testset "Parsing Types" begin
	println("Testing Parsing Types")
	data1 = [20 .*randn(5) for _ = 1:10]
	data2 = [rand(-10:10,5) for _ = 1:10]
	data = [data1;data2]
	labels = [rand(-10:10) for _ = 1:20]
	labels_f = randn(20)
	file = "tmp.libsvm"

	libsvmwrite(data, labels, file)
	d_read, l_read = libsvmread(file, valuetype=Float32)
	@test eltype(eltype(d_read)) == Float32
	@test d_read[1] isa SparseVector
	@test eltype(l_read) == Int
	readcorrectly(data, labels, d_read, l_read)

	d_read, l_read = libsvmread(file, labeltype=Float64)
	@test eltype(eltype(d_read)) == Float64
	@test d_read[1] isa SparseVector
	@test eltype(l_read) == Float64
	readcorrectly(data, labels, d_read, l_read)

	libsvmwrite(data, labels_f, file)
	for l in eachline(file)
		println(l)
	end
	d_read, l_read = libsvmread(file, labeltype=Float32, valuetype=Float32)
	@test eltype(eltype(d_read)) == Float32
	@test d_read[1] isa SparseVector
	@test eltype(l_read) == Float32
	readcorrectly(data, labels_f, d_read, l_read)

	d_read, l_read = libsvmread(file, labeltype=Float32, dense=true)
	@test eltype(eltype(d_read)) == Float64
	@test d_read[1] isa Vector
	@test eltype(l_read) == Float32
	readcorrectly(data, labels_f, d_read, l_read)
end

@testset "Selection" begin
	println("Testing Selection")
	data = [20 .*randn(10) for i = 1:30]
	labels = rand(-10:10, 30)
	file = "tmp.libsvm"

	libsvmwrite(data, labels, file)

	nEx, nFeat = libsvmsize(file)
	@test nEx == 30
	@test nFeat == 10

	selection = 1:nEx
	d_read, l_read = libsvmread(file, selection=selection, size=(nEx,nFeat))
	readcorrectly(data, labels, d_read, l_read)

	selection = 7:3:16
	d_read, l_read = libsvmread(file, selection=selection)
	readcorrectly(data[selection], labels[selection], d_read, l_read)

	selection = 13:13
	d_read, l_read = libsvmread(file, selection=selection, size=(nEx,nFeat))
	readcorrectly(data[selection], labels[selection], d_read, l_read)

	selection = 13:12
	d_read, l_read = libsvmread(file, selection=selection, size=(nEx,nFeat))
	readcorrectly(data[selection], labels[selection], d_read, l_read)

	selection = 2:2:21
	d_read, l_read = libsvmread(file, selection=selection, valuetype=Float32,
								labeltype=Float64, size=(nEx,nFeat))
	readcorrectly(data[selection], labels[selection], d_read, l_read)

end

@testset "Multi-label" begin
	println("Testing Multi-label")
	data = [20 .*randn(5) for i = 1:15]
	labels = [Tuple(randn(Float64,rand(1:6))) for _ = 1:15]
	file = "tmp.libsvm"

	libsvmwrite(data, labels, file, precision=2)

	for l in eachline(file)
		println(l)
	end

	nEx, nFeat = libsvmsize(file)
	@test nEx == 15
	@test nFeat == 5

	d_read, l_read = libsvmread(file, multilabel=true, labeltype=Float64)
	readcorrectly(data, labels, d_read, l_read, 1e2)
end


text:
using DifferentialEquations
f = (u,p,t) -> (p*u)
prob_ode_linear = ODEProblem(f,1/2,(0.0,1.0),1.01);

prob = prob_ode_linear
sol =solve(prob,Tsit5())
println(sol)

prob = ODEProblem(f,1/2,(0//1,1//1),101//100);
sol = solve(prob,RK4(),dt=1//2^(6),adaptive=false)
println(sol)

prob = ODEProblem(f,BigInt(1)//BigInt(2),(0//1,1//1),101//100);
sol =solve(prob,RK4(),dt=1//2^(6),adaptive=false)
println(sol[end])

prob_ode_biglinear = ODEProblem(f,big(1.0)/big(2.0),(big(0.0),big(1.0)),big(1.01))
sol =solve(prob_ode_biglinear,Tsit5())
println(sol[end])

using DoubleFloats
prob_ode_doublelinear = ODEProblem(f,Double64(1)/Double64(2),(Double64(0),Double64(1)),Double64(1.01))
sol =solve(prob_ode_doublelinear,Tsit5())
println(sol[end])

using ArbNumerics
prob_ode_arbfloatlinear = ODEProblem(f,ArbFloat(1)/ArbFloat(2),(ArbFloat(0.0),ArbFloat(1.0)),ArbFloat(1.01))
sol =solve(prob_ode_arbfloatlinear,Tsit5())
println(sol[end])

using DecFP
prob_ode_decfplinear = ODEProblem(f,Dec128(1)/Dec128(2),(Dec128(0.0),Dec128(1.0)),Dec128(1.01))
sol =solve(prob_ode_decfplinear,Tsit5())
println(sol[end]); println(typeof(sol[end]))

using Decimals
prob_ode_decimallinear = ODEProblem(f,[decimal("1.0")]./[decimal("2.0")],(0//1,1//1),decimal(1.01))
sol =solve(prob_ode_decimallinear,RK4(),dt=1/2^(6)) #Fails
println(sol[end]); println(typeof(sol[end]))

using DiffEqTutorials
DiffEqTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
syms rho_p rho_m
syms E_p E_m
syms K_p K_m
syms Phi_p Phi_m
syms dz
p_p=E_p-rho_p*(K_p+Phi_p);
p_m=E_m-rho_m*(K_m+Phi_m);
D=2/(rho_p+rho_m)*(p_p-p_m)/dz;

dDdrho_p=simplify(diff(D,rho_p));
dDdrho_m=simplify(diff(D,rho_m));

aa=3;


text:
uniform_beta_mixture = HeterogeneousMixture([uniform, beta])

@testset "fixed mixture of different distributions" begin

    @test !is_discrete(uniform_beta_mixture)
    @test has_output_grad(uniform_beta_mixture)
    @test has_argument_grads(uniform_beta_mixture) == (true, true, true, true, true)

    w1 = 0.4
    w2 = 0.6
    a = 0.0
    b = 3.4
    alpha = 2.3
    beta = 1.0

    # random
    x = uniform_beta_mixture([w1, w2], a, b, alpha, beta)

    # logpdf
    x = 0.123
    actual = logpdf(uniform_beta_mixture, x, [w1, w2], a, b, alpha, beta)
    expected = log(w1 * exp(logpdf(uniform, x, a, b)) + w2 * exp(logpdf(Gen.Beta(), x, alpha, beta)))
    @test isapprox(actual, expected)

    # test logpdf_grad against finite differencing
    args = (x, [w1, w2], a, b, alpha, beta)
    (x_grad, weights_grad, a_grad, b_grad, alpha_grad, beta_grad) = logpdf_grad(
        uniform_beta_mixture, args...)
    f = (x, weights, a, b, alpha, beta) -> logpdf(uniform_beta_mixture, x, weights, a, b, alpha, beta)

    @test isapprox(x_grad, finite_diff(f, args, 1, dx))
    @test isapprox(weights_grad[1], finite_diff_vec(f, args, 2, 1, dx))
    @test isapprox(weights_grad[2], finite_diff_vec(f, args, 2, 2, dx))
    @test isapprox(a_grad, finite_diff(f, args, 3, dx))
    @test isapprox(b_grad, finite_diff(f, args, 4, dx))
    @test isapprox(alpha_grad, finite_diff(f, args, 5, dx))
    @test isapprox(beta_grad, finite_diff(f, args, 6, dx))
end

mixture_of_normals = HomogeneousMixture(normal, [0, 0])
mixture_of_binomials = HomogeneousMixture(binom, [0, 0])
mixture_of_mvnormals = HomogeneousMixture(mvnormal, [1, 2])

@testset "mixture of normals" begin

    @test !is_discrete(mixture_of_normals)
    @test has_output_grad(mixture_of_normals)
    @test has_argument_grads(mixture_of_normals) == (true, true, true)

    w1 = 0.4
    w2 = 0.6
    mu1 = 0.0
    mu2 = 1.0
    std1 = 2.3
    std2 = 1.0

    # random
    x = mixture_of_normals([w1, w2], [mu1, mu2], [std1, std2])

    # logpdf
    x = 1.123
    actual = logpdf(mixture_of_normals, x, [w1, w2], [mu1, mu2], [std1, std2])
    expected = log(w1 * exp(logpdf(normal, x, mu1, std1)) + w2 * exp(logpdf(normal, x, mu2, std2)))
    @test isapprox(actual, expected)

    # test logpdf_grad against finite differencing
    args = (x, [w1, w2], [mu1, mu2], [std1, std2])
    (x_grad, weights_grad, mus_grad, stds_grad) = logpdf_grad(
        mixture_of_normals, args...)
    @test size(mus_grad) == (2,)
    @test size(stds_grad) == (2,)
    f = (x, weights, mus, stds) -> logpdf(mixture_of_normals, x, weights, mus, stds)
    @test isapprox(x_grad, finite_diff(f, args, 1, dx))
    @test isapprox(weights_grad[1], finite_diff_vec(f, args, 2, 1, dx))
    @test isapprox(weights_grad[2], finite_diff_vec(f, args, 2, 2, dx))
    @test isapprox(mus_grad[1], finite_diff_vec(f, args, 3, 1, dx))
    @test isapprox(mus_grad[2], finite_diff_vec(f, args, 3, 2, dx))
    @test isapprox(stds_grad[1], finite_diff_vec(f, args, 4, 1, dx))
    @test isapprox(stds_grad[2], finite_diff_vec(f, args, 4, 2, dx))

    # test that logpdf can be differentiated by ReverseDiff.jl
    # and test against finite differencing
    tp = Gen.new_tape()
    x_tracked = Gen.track(x, tp)
    w1_tracked = Gen.track(w1, tp)
    w2_tracked = Gen.track(w2, tp)
    mu1_tracked = Gen.track(mu1, tp)
    mu2_tracked = Gen.track(mu2, tp)
    std1_tracked = Gen.track(std1, tp)
    std2_tracked = Gen.track(std2, tp)
    lpdf_tracked = logpdf(mixture_of_normals,
        x_tracked, [w1_tracked, w2_tracked], [mu1_tracked, mu2_tracked], [std1_tracked, std2_tracked])
    Gen.deriv!(lpdf_tracked, 1.0)
    Gen.reverse_pass!(tp)
    @test isapprox(Gen.deriv(x_tracked), finite_diff(f, args, 1, dx))
    @test isapprox(Gen.deriv(w1_tracked), finite_diff_vec(f, args, 2, 1, dx))
    @test isapprox(Gen.deriv(w2_tracked), finite_diff_vec(f, args, 2, 2, dx))
    @test isapprox(Gen.deriv(mu1_tracked), finite_diff_vec(f, args, 3, 1, dx))
    @test isapprox(Gen.deriv(mu2_tracked), finite_diff_vec(f, args, 3, 2, dx))
    @test isapprox(Gen.deriv(std1_tracked), finite_diff_vec(f, args, 4, 1, dx))
    @test isapprox(Gen.deriv(std2_tracked), finite_diff_vec(f, args, 4, 2, dx))
end

@testset "mixture of binomial" begin

    @test is_discrete(mixture_of_binomials)
    @test !has_output_grad(mixture_of_binomials)
    @test has_argument_grads(mixture_of_binomials) == (true, false, true)

    w1 = 0.4
    w2 = 0.6
    n1 = 10
    n2 = 20
    p1 = 0.2
    p2 = 0.3

    # random
    x = mixture_of_binomials([w1, w2], [n1, n2], [p1, p2])

    # logpdf
    x = 4
    actual = logpdf(mixture_of_binomials, x, [w1, w2], [n1, n2], [p1, p2])
    expected = log(w1 * exp(logpdf(binom, x, n1, p1)) + w2 * exp(logpdf(binom, x, n2, p2)))
    @test isapprox(actual, expected)

    # test logpdf_grad against finite differencing
    args = (x, [w1, w2], [n1, n2], [p1, p2])
    (x_grad, weights_grad, ns_grad, ps_grad) = logpdf_grad(
        mixture_of_binomials, args...)
    @test ns_grad == nothing
    f = (x, weights, ns, ps) -> logpdf(mixture_of_binomials, x, weights, ns, ps)
    @test isapprox(weights_grad[1], finite_diff_vec(f, args, 2, 1, dx))
    @test isapprox(weights_grad[2], finite_diff_vec(f, args, 2, 2, dx))
    @test isapprox(ps_grad[1], finite_diff_vec(f, args, 4, 1, dx))
    @test isapprox(ps_grad[2], finite_diff_vec(f, args, 4, 2, dx))
end

@testset "mixture of multivariate normals" begin

    @test !is_discrete(mixture_of_mvnormals)
    @test has_output_grad(mixture_of_mvnormals)
    @test has_argument_grads(mixture_of_mvnormals) == (true, true, true)

    w1 = 0.4
    w2 = 0.6
    mu1 = [0.0, 0.0]
    mu2 = [1.0, 1.0]
    cov1 = [2.3 0.0; 0.0 2.3]
    cov2 = [1.0 0.0; 0.0 1.0]

    # random
    mus = cat(mu1, mu2, dims=2)
    covs = cat(cov1, cov2, dims=3)
    x = mixture_of_mvnormals([w1, w2], mus, covs)

    # logpdf
    x = [1.123, -2.123]
    actual = logpdf(mixture_of_mvnormals, x, [w1, w2], mus, covs)
    expected = log(w1 * exp(logpdf(mvnormal, x, mu1, cov1)) + w2 * exp(logpdf(mvnormal, x, mu2, cov2)))
    @test isapprox(actual, expected)

    # test logpdf_grad against finite differencing
    args = (x, [w1, w2], mus, covs)
    (x_grad, weights_grad, mus_grad, covs_grad) = logpdf_grad(
        mixture_of_mvnormals, args...)
    @test size(mus_grad) == (2, 2)
    @test size(covs_grad) == (2, 2, 2)
    f = (x, weights, mus, covs) -> logpdf(mixture_of_mvnormals, x, weights, mus, covs)
    @test isapprox(x_grad[1], finite_diff_vec(f, args, 1, 1, dx))
    @test isapprox(x_grad[2], finite_diff_vec(f, args, 1, 2, dx))
    @test isapprox(weights_grad[1], finite_diff_vec(f, args, 2, 1, dx))
    @test isapprox(weights_grad[2], finite_diff_vec(f, args, 2, 2, dx))
    for idx in 1:length(mus_grad)
        @test isapprox(mus_grad[idx], finite_diff_arr(f, args, 3, idx, dx), rtol=1e-2)
    end

    # to test logpdf_grad with respect to covariance matrices using finite differencing
    # we can't directly introduce differences into the matrix itself, because then
    # it leaves the set of positive definite matrices
    # instead, we define the covariance matrix as a rotated version of some covariance matrix
    # and differentiate with respect to the scalar amount of rotation
    rotmat(theta) = [cos(theta) -sin(theta); sin(theta) cos(theta)]

    @gen function f_cov_model((grad)(theta::Float64))
        cov = [2.0 0.0; 0.0 1.0]
        rot = rotmat(theta)
        x ~ mixture_of_mvnormals([w1, w2], mus, cat(cov, rot * cov * rot', dims=3))
    end

    theta = pi/2
    trace1, = generate(f_cov_model, (theta-dx,), choicemap((:x, x)))
    trace2, = generate(f_cov_model, (theta+dx,), choicemap((:x, x)))
    expected = (get_score(trace2) - get_score(trace1)) / (2 * dx)
    trace, = generate(f_cov_model, (theta,), choicemap((:x, x)))
    ((actual,), _, _) = choice_gradients(trace, select())
    @test isapprox(actual, expected)
end

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
struct OrBind
  pat1
  pat2
end

or_(a, b) = OrBind(a, b)
or_(p...) = foldl(or_, p)

function match_inner(pat::OrBind, ex, env)
  env′ = trymatch(pat.pat1, ex)
  env′ === nothing ? match(pat.pat2, ex, env) : merge!(env, env′)
end

function isor(ex)
  if isexpr(ex, :call)
    arg1 = ex.args[1]
    return arg1 isa Symbol && arg1 in (:or_, :|)  # "isa Symbol" check improves inferrability (#166)
  end
  return false
end

function ornew(ex)
  isor(ex) || return ex
  or_(ex.args[2:end]...)
end

subor(s) = s
subor(s::Symbol) = s
subor(s::Expr) = isor(s) ? subor(ornew(s)) : Expr(s.head, map(subor, s.args)...)
subor(s::OrBind) = OrBind(subor(s.pat1), subor(s.pat2))


text:
module Peacock

using Base: Number
using LinearAlgebra: include
using PyPlot
using LinearAlgebra, FFTW
using CUDA

include("utils.jl")
export diff_yee2, eye, diag_R, mask

include("geometry.jl")
export Geometry, Mu, Eps, MaterialTensor, TE, TM, Polarisation

include("plane_wave_basis.jl")
export BrillouinZoneCoordinate, FDFDBasis

include("modes.jl")
export Mode, HilbertSpace, get_field, Mode_FDFD, get_field_FDFD, HilbertSpace_FDFD

include("solver.jl")
export Solver, solve

include("plotting.jl")
export plot

include("FDFD.jl")

include("band_diagrams.jl")
export plot_band_diagram

include("wilson_loops.jl")
export plot_wilson_loop_winding

# Submodule that defines my commonly used crystals
include("Zoo.jl")

include("Chern_number.jl")
export Chern_number

include("solver_EPWE.jl")
export solver_EPWE, abcd

end # module


max_stars_repo_name: Soyukke/julia-paralell-examples

max_issues_repo_name: Soyukke/julia-paralell-examples

max_forks_repo_name: Soyukke/julia-paralell-examples


text:
include("definitions.jl")

function saveSparseData(sparseDataFilename::AbstractString, documents::Array{SparseDocument,1})
#save a dataset in the sparse document format
#one document per line, consisting of pairs of token index and token count, space separated.
    f = open(sparseDataFilename, "w");
    for doc in documents
        for i = 1:length(doc.tokenIndex)
           write(f, "$(doc.tokenIndex[i]) $(doc.tokenCount[i])");
           if i < length(doc.tokenIndex)
               write(f, " ");
           end
        end
        write(f,"\n");
    end
    close(f);
end


text:
using DSP

include("constants.jl")
include("utils.jl")
include("watersheds.jl")

function calc_xmask(lon,lat)
    dd_res = mean([abs(lon[1]-lon[2]),abs(lat[1]-lat[2])])
    lons,lats = meshgrid(lon,lat)
    sides = dd_to_meters.(lons,lats;scale=dd_res)

    xdim,ydim = size(sides)

    xmask = zeros(Float64, xdim,ydim)

    for i in 1:xdim, j in 1:ydim
        xmask[i,j] = sqrt( +((sides[i,j].^2)...)) / 2
    end

    return xmask
end

function make_irf(xmask::Float64; diffusion::Float64=800.0,velocity::Float64=1.5, max_runoff_days::Int64=2)
    """
    Function for Calculate the impulse response function for a grid cells using equation 15 from Lohmann, et al. (1996) Tellus article
    Arguments:
        xmask -- longest path for any location within the grid cell to reach
                 the stream channel [m]
        diff -- diffusion value
        velo -- overland flow velocity value
    Returns:
        irf -- normalized impulse response function for a grid cell
    Keywords:
        diff -- the diffusivity parameter (default 800) [m^2/s]
        velo -- the flow velocity parameter (default 1.5) [m/s]
    """
    # get number of hours per impulse (typically 48 hours)
    nt = max_runoff_days * HOURSPERDAY
    step = collect(1:1:nt) #number of time steps

    # convert time steps to seconds
    time = @. step * SECSPERHOUR

    exponent = @. -1 * (velocity * time - xmask)^2 / (4.0 * diffusion * time)
    #impulse response function
    green = @. xmask / (2.0 * time * sqrt(pi * time * diffusion)) * exp(exponent)

    irf = green / sum(green) #normalize the output IRF

    return irf

end

function make_river_irf(irfs; max_runoff_days::Int64=2, max_basin_days::Int64=50)

    t_cell = Int(max_runoff_days * HOURSPERDAY)
    t_uh = Int(max_basin_days * HOURSPERDAY)

    n_cells = length(irfs)

    river_irf = zeros(Float64,t_uh,n_cells)

    river_irf[1:t_cell, 1] = irfs[1]

    for i in 2:n_cells
        irf_temp = conv(river_irf[:, i-1], irfs[i])

        river_irf[:,i] = irf_temp[1:t_uh] / sum(irf_temp[1:t_uh])
    end

    return river_irf

end

function make_grid_uh(river_irfs,uh_box; max_runoff_days::Int64=2, max_basin_days::Int64=50)

    t_cell = Int(max_runoff_days * HOURSPERDAY)
    t_uh = Int(max_basin_days * HOURSPERDAY)

    uh_shape = size(river_irfs)

    n_cells = uh_shape[end]

    # prealloc
    uh = zeros(Float64,uh_shape...)
    irf_temp = zeros(Float64, (t_uh + t_cell))

    # set init values
    irf_temp[1:length(uh_box)] = uh_box[:]
    uh[:,1] = irf_temp[1:t_uh]

    for i in 2:n_cells
        irf_temp = conv(uh_box, river_irfs[:,i-1])
        uh[:,i] = irf_temp[1:t_uh] / sum(irf_temp[1:t_uh])
    end

    return uh
end

function rout(runoff,unit_hydrograph,area)

    n_cells = size(unit_hydrograph)[end]
    days = size(runoff)[end]
    buffer = size(unit_hydrograph)[1]

    basinuh = zeros(Float64,n_cells,days) # ...basin wide UH
    griduh = zeros(Float64,days+buffer,buffer)

    for i in 1:n_cells

        factor = (area[i] / (SECSPERDAY * MMPERMETER))

        for u in 1:buffer
            # calculate the percent contributions from pixel to outlet at each time period
            griduh[u:days+u-1,u] = (runoff[i,:] * factor) * unit_hydrograph[u]
        end

        basinuh[i,:] = sum(griduh[1:days,:],dims=2)

    end

    finalq = dropdims(sum(basinuh,dims=1),dims=1)

    return finalq

end


text:
function writeJSON(json_string,directory,tileRef)
    open("$(directory)/$(tileRef).json","w") do f
        write(f,json_string)
    end
end
function readJSON(directory,tileRef)
    io = open("$(directory)/$(tileRef).json")
    content = readlines(io)
    close(io)
    data_dict = JSON.parse(content[1])
    data = DataFrame(data_dict)
    # As 'Missing's come back as 'Nothing's from JSON, this is corrected below
    prep_df(data,col) = data[findall(isnothing.(data[:,col])),col] .= missing
    prep_df.((data,),collect(1:ncol(data)))
    return data
end


text:
using Flux
using Flux: throttle, nfan, glorot_uniform, glorot_normal, kaiming_normal, kaiming_uniform, stack, unstack
using StatsBase: var, std
using Random
using Test

@testset "Throttle" begin
  @testset "default behaviour" begin
    a = []
    f = throttle(()->push!(a, time()), 1, leading=true, trailing=false)
    f()
    f()
    f()
    sleep(1.01)
    @test length(a) == 1
  end

  @testset "leading behaviour" begin
    a = []
    f = throttle(()->push!(a, time()), 1, leading=true, trailing=false)
    f()
    @test length(a) == 1
    f()
    @test length(a) == 1
    sleep(1.01)
    f()
    @test length(a) == 2
  end

  @testset "trailing behaviour" begin
    a = []
    f = throttle(()->push!(a, time()), 1, leading=false, trailing=true)
    f()
    @test length(a) == 0
    f()
    @test length(a) == 0
    sleep(1.01)
    @test length(a) == 1
  end

  @testset "arguments" begin
    a = []
    f = throttle((x)->push!(a, x), 1, leading=true, trailing=true)
    f(1)
    @test a == [1]
    f(2)
    @test a == [1]
    f(3)
    @test a == [1]
    sleep(1.01)
    @test a == [1, 3]
  end
end

@testset "Initialization" begin
  # Set random seed so that these tests don't fail randomly
  Random.seed!(0)

  @testset "Fan in/out" begin
    @test nfan() == (1, 1) #For a constant
    @test nfan(100) == (1, 100) #For vector
    @test nfan(100, 200) == (200, 100) == nfan((100, 200)) #For Dense layer
    @test nfan(2, 30, 40) == (2 * 30, 2 * 40) #For 1D Conv layer
    @test nfan(2, 3, 40, 50) == (2 * 3 * 40, 2 * 3 * 50) #For 2D Conv layer
    @test nfan(2, 3, 4, 50, 60) == (2 * 3 * 4 * 50, 2 * 3 * 4 * 60) #For 3D Conv layer
  end

  @testset "glorot" begin
    # glorot_uniform and glorot_normal should both yield a kernel with
    # variance ≈ 2/(fan_in + fan_out)
    for dims ∈ [(1000,), (100, 100), (100, 400), (2, 3, 32, 64), (2, 3, 4, 32, 64)]
      for init ∈ [glorot_uniform, glorot_normal]
        v = init(dims...)
        fan_in, fan_out = nfan(dims...)
        σ2 = 2 / (fan_in + fan_out)
        @test 0.9σ2 < var(v) < 1.1σ2
        @test eltype(v) == Float32
      end
    end
  end

  @testset "kaiming" begin
    # kaiming_uniform should yield a kernel in range [-sqrt(6/n_out), sqrt(6/n_out)]
    # and kaiming_normal should yield a kernel with stddev ~= sqrt(2/n_out)
    for (n_in, n_out) in [(100, 100), (100, 400)]
      v = kaiming_uniform(n_in, n_out)
      σ2 = sqrt(6/n_out)
      @test -1σ2  < minimum(v) < -0.9σ2
      @test 0.9σ2  < maximum(v) < 1σ2
      @test eltype(v) == Float32

      v = kaiming_normal(n_in, n_out)
      σ2 = sqrt(2/n_out)
      @test 0.9σ2 < std(v) < 1.1σ2
      @test eltype(v) == Float32
    end
  end
end

@testset "Params" begin
  m = Dense(10, 5)
  @test size.(params(m)) == [(5, 10), (5,)]
  m = RNN(10, 5)
  @test size.(params(m)) == [(5, 10), (5, 5), (5,), (5,)]

  # Layer duplicated in same chain, params just once pls.
  c = Chain(m, m)
  @test size.(params(c)) == [(5, 10), (5, 5), (5,), (5,)]

  # Self-referential array. Just want params, no stack overflow pls.
  r = Any[nothing,m]
  r[1] = r
  @test size.(params(r)) == [(5, 10), (5, 5), (5,), (5,)]
end

@testset "Basic Stacking" begin
  x = randn(3,3)
  stacked = stack([x, x], 2)
  @test size(stacked) == (3,2,3)
end

@testset "Precision" begin
  m = Chain(Dense(10, 5, relu), Dense(5, 2))
  x = rand(10)
  @test eltype(m[1].W) == Float32
  @test eltype(m(x)) == Float32
  @test eltype(f64(m)(x)) == Float64
  @test eltype(f64(m)[1].W) == Float64
  @test eltype(f32(f64(m))[1].W) == Float32
end

@testset "Stacking" begin
  stacked_array=[ 8 9 3 5; 9 6 6 9; 9 1 7 2; 7 4 10 6 ]
  unstacked_array=[[8, 9, 9, 7], [9, 6, 1, 4], [3, 6, 7, 10], [5, 9, 2, 6]]
  @test unstack(stacked_array, 2) == unstacked_array
  @test stack(unstacked_array, 2) == stacked_array
  @test stack(unstack(stacked_array, 1), 1) == stacked_array
end


text:
function check_skip_allele(a)
    if occursin(a, skip_chars)
        return true
    elseif occursin(a, valid_alleles)
        return false
    else
        throw(DomainError(a, string("$a not a valid allele")))
    end
end

function count_seq_diffs(seq1, seq2)
    if length(seq1) != length(seq2)
        throw(ErrorException(string(
            "contig lengths differ: ",
            "$(length(seq1)) != $(length(seq2))")))
    end

    total_bases = 0
    n_diffs = 0
    n_skipped = 0

    seq1 = uppercase(seq1)
    seq2 = uppercase(seq2)
    for (a1, a2) in zip(seq1, seq2)
        total_bases += 1
        if check_skip_allele(a1) || check_skip_allele(a2)
            n_skipped += 1
        elseif a1 != a2
            n_diffs += 1
        end
    end
    return [total_bases-n_skipped, n_diffs]
end

function fasta_edit_distances(fasta_name)
    names = []
    seqs = []
    open(fasta_name) do fasta
        for record in FASTA.Reader(fasta)
            name = FASTA.identifier(record)
            seq = String(sequence(LongCharSeq, record))
            push!(names, name)
            push!(seqs, seq)
        end
    end

    jobs = []
    for i = 1:(length(seqs)-1)
        for j = (i+1):length(seqs)
            push!(jobs, [i j])
        end
    end

    results = Array{Integer}(undef, length(jobs), 2)
    Threads.@threads for idx = 1:length(jobs)
        i = jobs[idx][1]
        j = jobs[idx][2]
        results[idx, :] = count_seq_diffs(seqs[i], seqs[j])
    end

    return DataFrame(
        Contig1 = [names[j[1]] for j in jobs],
        Contig2 = [names[j[2]] for j in jobs],
        SharedGenomeLen = results[:, 1],
        NumDiffs = results[:, 2]
    )
end

export fasta_edit_distances


text:
get_part_id(comm::MPI.Comm) = MPI.Comm_rank(comm)+1
num_parts(comm::MPI.Comm) = MPI.Comm_size(comm)

struct MPIBackend <: AbstractBackend end

const mpi = MPIBackend()

function get_part_ids(b::MPIBackend,nparts::Integer)
  comm = MPI.Comm_dup(MPI.COMM_WORLD)
  @notimplementedif num_parts(comm) != nparts
  MPIData(get_part_id(comm),comm,(nparts,))
end

function get_part_ids(b::MPIBackend,nparts::Tuple)
  comm = MPI.Comm_dup(MPI.COMM_WORLD)
  @notimplementedif num_parts(comm) != prod(nparts)
  MPIData(get_part_id(comm),comm,nparts)
end

function prun(driver::Function,b::MPIBackend,nparts)
  if !MPI.Initialized()
    MPI.Init()
  end
  #try
    part = get_part_ids(b,nparts)
    driver(part)
  #finally
  #  MPI.Finalize()
  #end
end

struct MPIData{T,N} <: AbstractPData{T,N}
  part::T
  comm::MPI.Comm
  size::NTuple{N,Int}
  function MPIData(
    part::T,
    comm::MPI.Comm,
    size::NTuple{N,<:Integer}) where {T,N}
    @assert num_parts(comm) == prod(size)
    new{T,N}(part,comm,size)
  end
end

Base.size(a::MPIData) = a.size
get_part_id(a::MPIData) = get_part_id(a.comm)
get_backend(a::MPIData) = mpi
i_am_main(a::MPIData) = get_part_id(a.comm) == MAIN

function get_part_ids(a::MPIData)
  MPIData(get_part_id(a.comm),a.comm,a.size)
end

function Base.iterate(a::MPIData)
  next = iterate(a.part)
  if next == nothing
    return next
  end
  item, state = next
  MPIData(item,a.comm,a.size), state
end

function Base.iterate(a::MPIData,state)
  next = iterate(a.part,state)
  if next == nothing
    return next
  end
  item, state = next
  MPIData(item,a.comm,a.size), state
end

function map_parts(task,args::MPIData...)
  @assert length(args) > 0
  @assert all(a->a.comm===first(args).comm,args)
  parts_in = map(a->a.part,args)
  part_out = task(parts_in...)
  MPIData(part_out,first(args).comm,first(args).size)
end

#TODO Tables are not displayed correctly
function Base.show(io::IO,k::MIME"text/plain",data::MPIData)
  MPI.Barrier(data.comm)
  str = """
  On part $(get_part_id(data)) of $(num_parts(data)):
  $(data.part)
  """
  print(io,str)
  MPI.Barrier(data.comm)
end

get_part(a::MPIData) = a.part

function get_part(a::MPIData,part::Integer)
  rcv = MPI.Bcast!(Ref(copy(a.part)),part-1,a.comm)
  rcv[]
end

function get_part(a::MPIData{<:AbstractVector},part::Integer)
  l = map_parts(length,a)
  l_at_part = get_part(l,part)
  if get_part_id(a) == part
    rcv = a.part
  else
    rcv = similar(a.part,eltype(a.part),l_at_part)
  end
  MPI.Bcast!(rcv,part-1,a.comm)
  rcv
end

function gather!(rcv::MPIData,snd::MPIData)
  @assert rcv.comm === snd.comm
  if get_part_id(snd) == MAIN
    @assert length(rcv.part) == num_parts(snd)
    rcv.part[MAIN] = snd.part
    MPI.Gather!(MPI.IN_PLACE,MPI.UBuffer(rcv.part,1),MAIN-1,snd.comm)
  else
    MPI.Gather!(Ref(snd.part),nothing,MAIN-1,snd.comm)
  end
  rcv
end

function gather!(rcv::MPIData{<:Table},snd::MPIData)
  @assert rcv.comm === snd.comm
  if get_part_id(snd) == MAIN
    @assert length(rcv.part) == num_parts(snd)
    kini = rcv.part.ptrs[MAIN]
    kend = rcv.part.ptrs[MAIN+1]-1
    rcv.part.data[kini:kend] = snd.part
    counts = ptrs_to_counts(rcv.part.ptrs)
    MPI.Gatherv!(MPI.IN_PLACE,MPI.VBuffer(rcv.part.data,counts),MAIN-1,snd.comm)
  else
    MPI.Gatherv!(snd.part,nothing,MAIN-1,snd.comm)
  end
  rcv
end

function gather_all!(rcv::MPIData,snd::MPIData)
  @assert rcv.comm === snd.comm
  @assert length(rcv.part) == num_parts(snd)
  MPI.Allgather!(Ref(snd.part),MPI.UBuffer(rcv.part,1),snd.comm)
  rcv
end

function gather_all!(rcv::MPIData{<:Table},snd::MPIData)
  @assert rcv.comm === snd.comm
  @assert length(rcv.part) == num_parts(snd)
  counts = ptrs_to_counts(rcv.part.ptrs)
  MPI.Allgatherv!(snd.part,MPI.VBuffer(rcv.part.data,counts),snd.comm)
  rcv
end

function scatter(snd::MPIData)
  if get_part_id(snd) == MAIN
    part = snd.part[MAIN]
    MPI.Scatter!(MPI.UBuffer(snd.part,1),MPI.IN_PLACE,MAIN-1,snd.comm)
  else
    rcv = Vector{eltype(snd.part)}(undef,1)
    MPI.Scatter!(nothing,rcv,MAIN-1,snd.comm)
    part = rcv[1]
  end
  MPIData(part,snd.comm,snd.size)
end

function scatter(snd::MPIData{<:Table})
  counts_main = map_parts(ptrs_to_counts∘get_ptrs,snd)
  counts_scat = scatter(counts_main)
  if get_part_id(snd) == MAIN
    buf = MPI.VBuffer(snd.part.data,counts_main.part)
    MPI.Scatterv!(buf,MPI.IN_PLACE,MAIN-1,snd.comm)
    rcv = snd.part[MAIN]
  else
    rcv = eltype(snd.part)(undef,counts_scat.part)
    MPI.Scatterv!(nothing,rcv,MAIN-1,snd.comm)
  end
  MPIData(rcv,snd.comm,snd.size)
end

function scatter(snd::MPIData{<:AbstractVector{<:AbstractVector}})
  snd_table = map_parts(Table,snd)
  scatter(snd_table)
end

function emit(snd::MPIData)
  MPIData(get_main_part(snd),snd.comm,snd.size)
end

function async_exchange!(
  data_rcv::MPIData,
  data_snd::MPIData,
  parts_rcv::MPIData,
  parts_snd::MPIData,
  t_in::MPIData)

  @check parts_rcv.comm === data_rcv.comm
  @check parts_rcv.comm === data_snd.comm
  @check parts_rcv.comm === t_in.comm
  comm = parts_rcv.comm
  s = parts_rcv.size

  t0 = t_in.part

  t1 = @async begin

    req_all = MPI.Request[]
    wait(schedule(t0))

    for (i,part_rcv) in enumerate(parts_rcv.part)
      rank_rcv = part_rcv-1
      buff_rcv = view(data_rcv.part,i:i)
      tag_rcv = part_rcv
      reqr = MPI.Irecv!(buff_rcv,rank_rcv,tag_rcv,comm)
      push!(req_all,reqr)
    end

    for (i,part_snd) in enumerate(parts_snd.part)
      rank_snd = part_snd-1
      buff_snd = view(data_snd.part,i:i)
      tag_snd = get_part_id(comm)
      reqs = MPI.Isend(buff_snd,rank_snd,tag_snd,comm)
      push!(req_all,reqs)
    end

    return req_all
  end

  t2 = @task begin
    req_all = fetch(t1)
    MPI.Waitall!(req_all)
  end

  t_out = MPIData(t2,comm,s)
  t_out
end

function async_exchange!(
  data_rcv::MPIData{<:Table},
  data_snd::MPIData{<:Table},
  parts_rcv::MPIData,
  parts_snd::MPIData,
  t_in::MPIData)

  @check parts_rcv.comm === data_rcv.comm
  @check parts_rcv.comm === data_snd.comm
  @check parts_rcv.comm === t_in.comm
  comm = parts_rcv.comm
  s = parts_rcv.size

  t0 = t_in.part

  t1 = @async begin

    req_all = MPI.Request[]
    wait(schedule(t0))

    for (i,part_rcv) in enumerate(parts_rcv.part)
      rank_rcv = part_rcv-1
      ptrs_rcv = data_rcv.part.ptrs
      buff_rcv = view(data_rcv.part.data,ptrs_rcv[i]:(ptrs_rcv[i+1]-1))
      tag_rcv = part_rcv
      reqr = MPI.Irecv!(buff_rcv,rank_rcv,tag_rcv,comm)
      push!(req_all,reqr)
    end

    for (i,part_snd) in enumerate(parts_snd.part)
      rank_snd = part_snd-1
      ptrs_snd = data_snd.part.ptrs
      buff_snd = view(data_snd.part.data,ptrs_snd[i]:(ptrs_snd[i+1]-1))
      tag_snd = get_part_id(comm)
      reqs = MPI.Isend(buff_snd,rank_snd,tag_snd,comm)
      push!(req_all,reqs)
    end

    return req_all
  end

  t2 = @task begin
    req_all = fetch(t1)
    MPI.Waitall!(req_all)
  end

  t_out = MPIData(t2,comm,s)
  t_out
end


text:
################################################################################
# MIT License
#
# Copyright (c) 2021 Hajime Nakagami<nakagami@gmail.com>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
################################################################################

mutable struct Arc4
    state::Vector{UInt8}
    x::UInt64
    y::UInt64

    function Arc4(key::Vector{UInt8})
        state::Vector{UInt8} = []
        for i in 0:255
            push!(state, i)
        end
        @assert length(state) == 256

        index1::UInt64 = 0
        index2::UInt64 = 0

        for i in 0:255
            index2 = (key[index1+1] + state[i+1] + index2) % 256
            (state[i+1], state[index2+1]) = (state[index2+1], state[i+1])
            index1 = (index1 + 1) % length(key)
        end

        new(state, 0, 0)
    end
end

function translate(arc4::Arc4, plain::Vector{UInt8})
    state = arc4.state
    enc::Vector{UInt8} = []
    for i in 1:length(plain)
        arc4.x = (arc4.x + 1) %256
        arc4.y = (arc4.y + state[arc4.x+1]) % 256
        state[arc4.x+1], state[arc4.y+1] = state[arc4.y+1], state[arc4.x+1]
        push!(enc, xor(plain[i], state[((state[arc4.x+1]+state[arc4.y+1]) % 256) + 1]))
    end
    enc
end


text:
mutable struct Shared
    df
    selectvalues::Vector{SpinBoxType}
    selectlist::Vector{Column}
    plotvalues::Vector{ComboBoxType}
    plotkwargs::TextBoxEntry
    smoother::SliderEntry
    splitting_var::ComboBoxType
    plt::Plots.Plot{Plots.GRBackend}
end
shared = Shared([], SpinBoxType[], Column[], ComboBoxType[],
    TextBoxEntry(""), SliderEntry(0.0), ComboBoxType("", [""]), plot())

"""
    build_window(; kwargs...)

Starts the GUI asking for a suitable csv file.
"""
function build_window(; kwargs...)
    @qmlapp joinpath(Pkg.dir("PlugAndPlot","src"), "QML", "choose_file.qml")
    path = TextBoxEntry("")
    @qmlset qmlcontext().choose = path
    exec()
    return build_window(path.value; kwargs...)
end

"""
    build_window(datafile::AbstractString; kwargs...)

Reads a csv file and starts `build_window` on the corresponding DataFrame
"""
build_window(datafile::AbstractString; kwargs...) = build_window(JuliaDB.loadtable(datafile); kwargs...)

"""
    build_window(dataset; nbox = 5)

Creates a GUI to analyze a data table interactively. Data can be selected either
on continuous columns, with SpinBoxes or on discrete columns with checkboxes,
provided there are less than `nbox` entries.
"""
build_window(dataset; kwargs...) = build_window(table(dataset); kwargs...)

function build_window(dataset::IndexedTables.AbstractIndexedTable; nbox = 5)
    shared.df = dataset
    shared.selectlist = [Column(string(name), string.(union(columns(shared.df, name))))
        for name in colnames(shared.df) if (1 < length(union(columns(shared.df, name))) < nbox)]

    shared.selectvalues = [SpinBoxType(string(name), Float64.([extrema(columns(shared.df, name))...]))
        for name in colnames(shared.df) if length(union(columns(shared.df, name))) >= nbox &&
        eltype(columns(shared.df, name)) <: Real]
    # run QML window
    qml_engine = init_qmlapplicationengine()
    @qmlset qmlcontext()._selectlist = ListModel(shared.selectlist)
    @qmlset qmlcontext()._selectvalues = ListModel(shared.selectvalues)
    shared.plotvalues = get_plotvalues(shared.df)
    @qmlset qmlcontext()._plotvalues = ListModel(shared.plotvalues)
    @qmlset qmlcontext().choose = shared.plotkwargs
    @qmlset qmlcontext().smoother = shared.smoother
    splitting_options = ["Don't split"; string.("splitby:",colnames(shared.df))]
    shared.splitting_var = ComboBoxType("spitting_var", splitting_options)
    @qmlset qmlcontext().splitting_var = shared.splitting_var
    qmlfunction("do_plot", PlugAndPlot.do_plot)
    qmlfunction("save_plot", PlugAndPlot.save_plot)
    qmlfunction("check_equality", PlugAndPlot.check_equality)
    qml_file = joinpath(Pkg.dir("PlugAndPlot","src"), "QML", "gui.qml")
    QML.load(qml_engine,qml_file)

    exec()
    return
end

function do_plot(d::JuliaDisplay, width, height, in_place = false)
    gr(grid = false, size=(Int64(round(width)),Int64(round(height))))
    get_plot!(shared, in_place)
    display(d, shared.plt)
end

save_plot(savename::AbstractString) =  savefig(shared.plt, savename)

function check_equality()
    (shared.plotvalues[1].chosen_value == shared.plotvalues[2].chosen_value) &&
    (shared.plotvalues[1].text_info == shared.plotvalues[2].text_info)
end


text:
using LinearAlgebra: svd, svdvals, cond, pinv, I, diagm
using Statistics
using Distributions
using Random
using JLD

include("../../../../Potentials.jl/src/Potentials.jl")

# Parameters
kb              = 8.617e-5 # eV/K
Tc              = 120.0    # K
μ               = 39.948   # u
ϵ               = Tc*kb    # 0.0103404 eV
σ               = 3.4      # Angstroms Å

# Use reduced units
ϵ               = 1.0
σ               = 1.0
μ               = 1.0
Tc_kb           = 1.0
Tc              = 1.0
Temp = 0.25
samples = load("examples/Argon7/lj_estimation/samples_lj_TEMP_25.jld")
map = samples["map"]
lj = Potentials.LennardJones(map.ϵ, map.σ)
c = Potentials.Configuration("examples/Argon7/DATA_training/DATA_start", atom_names = [:Ar], radii = [3.5], weights = [1.0])

Nt = 1000000
dT = 1000
dt = 0.01
path = "examples/Argon7/lj_estimation/TEMP_$(Int(100*Temp))/"
# mkdir(path)
# mkdir(path*"/MAP/")
for seed = 1:100
    if seed % 10 == 0
        println("seed $seed")
    end
    path = "examples/Argon7/lj_estimation/TEMP_$(Int(100*Temp))/MAP/"
    path = path * "seed_$seed/"
    mkdir(path)
    @time rdf, pe, msd = Potentials.run_md(c, lj, Nt, path; dt = dt, dim = 2, Temp = Temp, dT = dT, seed = 60)
end

# T = dT:dT:Nt
# num_configs = length(T)
# r = Vector{Potentials.Configuration}(undef, num_configs)
# for (i,t) in enumerate(T)
#     file = "DATA.$t"
#     c_temp = Potentials.Configuration("examples/Argon7/temp/"*file; atom_names = [:Ar],
#                     radii = [3.5], weights = [1.0])
#     r[i] = c_temp
# end

# for seed = 2:100
#     path = "examples/Argon7/DATA_md/LJ/seed_$seed/"
#     mkdir(path)
#     Nt = 50000
#     dT = 100
#     println("Seed $seed")
#     @time rdf, pe, msd = Potentials.run_md_lj(c, Nt, path; dim = 2, Temp = 0.75, dT = dT, seed = seed)
# end


text:
using ..Photon: getContext

"""
Datasets are repsonsible for loading a single sample. They need to have length and
index methods implemented. Hence, they can all be passed to a Dataloader which can load
multiple samples parallelly using threading.

Pleae note that in case your whole dataset would fit in memory, you can feed it
directly to the fit! function (so a dataset is no requirement)

```julia
X = [randn(10,10,16) for i in 1:100]
Y = [randn(1,16) for i in 1:100]
fit!(workout, zip(X,Y))
```

Otherwise the combination of dataset/dataloader is your best bet.
"""
abstract type Dataset end

"""
A dataset that contains random generated samples. Ideal for quick testing
of a model. The random values will be drawn from a normal distribution.

You need to provide the shape of X, Y and the number of samples that the
dataset should contain. Optionally you can specify a sleep value to simulate IO blocking.

# Usage

```julia
xshape = (28,28,1)
yshape = (10,)
ds = TestDataset(xshape, yshape, 60000)

ds = TestDataset((100,), (1,), 100, sleep=0.1)
```
"""
struct TestDataset <: Dataset
    x::Vector{Array}
    y::Vector{Array}
	sleep::Float64

	function TestDataset(shapeX::Tuple, shapeY::Tuple, samples::Int; sleep=0)
		dtype = getContext().dtype
		x = [randn(dtype, shapeX...) for _ in 1:samples]
		y = [randn(dtype, shapeY...) for _ in 1:samples]
		new(x,y,sleep)
	end
end

Base.length(ds::TestDataset) = length(ds.x)

function Base.getindex(ds::TestDataset, idx)
	ds.sleep > 0 && sleep(ds.sleep)
	return (ds.x[idx], ds.y[idx])
end

"""
Base dataset that contains two vectors.
"""
struct ArrayDataset{A,B} <: Dataset
    x::A
    y::B
end

Base.length(ds::ArrayDataset) = length(ds.x)
Base.getindex(ds::ArrayDataset, idx) = return (ds.x[idx], ds.y[idx])

"""
Dataset that loads data from a JLD(2) file. The format
that is stored is expected to be in the shape of: (X,Y)

Note: Tested with JLD2 library only.

Example:
========

```julia
using JLD2

jldopen("example.jld2", "w") do file
    file["image1"] = (randn(Float32,28,28,1), rand(0:9,1))
    file["image2"] = (randn(Float32,28,28,1), rand(0:9,1))
end

f = jldopen("example.jld2", "r")
ds = JLD2Dataset(f)
ds[1]
```
"""
struct JLDDataset <: Dataset
	f
    keys

	JLDDataset(f) = new(f,keys(f))
end

Base.length(ds::JLDDataset) = length(ds.keys)

function Base.getindex(ds::JLDDataset, idx)
	key = ds.keys[idx]
	ds.f[key]
end

"""
Dataset that loads data from a JuliaDB. Not yet implemented.
"""
struct JuliaDBDataset{A,B} <: Dataset
    filenames::A
    labels::B
end

Base.length(ds::JuliaDBDataset) = length(ds.filenames)

function Base.getindex(ds::JuliaDBDataset, idx)
	# TODO
end

"""
Dataset that loads an single image from a file and optionally resizes the image.
The labels are passed as is.

# Usage

```julia
ds = ImageDataset(filenames, labels, resize=(200,200))
```
"""
struct ImageDataset <: Dataset
    filenames::Vector{String}
    labels::Vector
	resize::Union{Nothing,Tuple}

	function ImageDataset(filenames, labels; resize=nothing)
		@assert length(filenames) == length(labels)
		try
			@eval import Images
			@eval import ImageMagick
		catch
			@warn "Package Images or ImageMagick not installed"
		end
		new(filenames, labels, resize)
	end

end

Base.length(ds::ImageDataset) = length(ds.filenames)

function Base.getindex(ds::ImageDataset, idx)
	filename = ds.filenames[idx]
	# img = Images.load(filename) has multi-thread issue
    img = ImageMagick.load(filename)
	img = Images.channelview(img)
	img = permutedims(img, [2,3,1])
	img = convert(Array{Float32}, img)
	if ds.resize !== nothing
		img = Images.imresize(img, ds.resize...)
	end
	return (img, ds.labels[idx])
end

"""
Dataset that retrieves is data from a dataframe. The provided column names
for X and Y can be either a single Symbol or a Vector of Symbols.

# Usage

```julia
df = DataFrame(randn(4,20))
ds = DFDataset(df, :x1, :x2)
ds = DFDataset(df, [:x1, :x3], [:x2,:x4])
```
"""
struct DFDataset <: Dataset
	df
    X::Vector{Symbol}
	Y::Vector{Symbol}

	DFDataset(df,X,Y) = new(df, makeArray(X), makeArray(Y))
end

Base.length(ds::DFDataset) = size(df,1)

function Base.getindex(ds::DFDataset, idx)
	df = ds.df
	(Vector(df[idx, ds.X]), Vector(df[idx, ds.Y]))
end


text:
export entropy,spectrum
"""
    Calculate the von-Neumann entanglement entropy of an MPS at bond_id
"""
function entropy(mps::AbstractMPS,bond_id::Int)
    S = spectrum(mps,bond_id)
    S = S./sqrt(sum(S.*S))  #normalize MPS
    S = S .+ 1E-100
    return sum(- S.*S .* log.(S.*S))
end

"""
    Calculate von-Neumann entanglement entropy of every bonds
"""
function entropy(mps::AbstractMPS,epsilon=1E-13)
    function spectrum2entropy(S::AbstractVector)
        t = S./sqrt(sum(S.*S))  #normalize MPS
        t = t .+ 1E-100
        return sum(- t.*t .* log.(t.*t))
    end
    # Calculate the von-Neumann entanglement entropy of an MPS
    #from left to right, qr
    v_entropy = zeros(Float64,mps.L-1)
    sweep!(mps,mps.L-1,LeftNormalization())
    #print (mps.bdim)
    #from right to left, svd
    for site = mps.L:-1:2
        S = single_tensor_spectrum!(mps,site)
        v_entropy[site-1] = spectrum2entropy(S)
    end
    return v_entropy
end

"""
Calculate Spectrum in MPS: from site-bond_id to site-bond_id+1
"""
function spectrum(mps::AbstractMPS,bond_id::Int)
    res = 0.0
    @assert bond_id>0 && bond_id <= mps.L
    normalization!(mps,bond_id,MixNormalization())
    U,S,V = svd(reshape(mps[bond_id],mps.bdim[bond_id-1]*mps.S,mps.bdim[bond_id]))
    return S
end

max_stars_repo_name: qiyang-ustc/JMPS

max_issues_repo_name: qiyang-ustc/JMPS

max_forks_repo_name: qiyang-ustc/JMPS


text:
mutable struct SpglibDataset
    spacegroup_number::Cint
    hall_number::Cint
    international_symbol::NTuple{11, UInt8}
    hall_symbol::NTuple{17, UInt8}
    choice::NTuple{6, UInt8}
    transformation_matrix::NTuple{9, Float64}
    origin_shift::NTuple{3, Float64}
    n_operations::Cint
    rotations::Ptr{NTuple{9, Cint}}
    translations::Ptr{NTuple{3, Float64}}
    n_atoms::Cint
    wyckoffs::Ptr{Cint}
    site_symmetry_symbols::Ptr{Tuple{7, UInt8}}
    equivalent_atoms::Ptr{Cint}
    mapping_to_primitive::Ptr{Cint}
    n_std_atoms::Cint
    std_lattice::NTuple{9, Float64}
    std_types::Ptr{Cint}
    std_positions::Ptr{NTuple{3, Float64}}
    std_rotation_matrix::NTuple{9, Float64}
    std_mapping_to_primitive::Ptr{Cint}
    pointgroup_symbol::NTuple{6, UInt8}
end

function spg_get_dataset(lattice::Array{Float64, 2},
                         positions::Array{Float64, 2},
                         types::Array{Int64, 1},
                         num_atom::Int64,
                         symprec::Float64)
    num_atom = Base.convert(Cint, num_atom)
    types = Base.convert(Array{Cint, 1}, types)

    return unsafe_load(ccall( (:spg_get_dataset, libsymspg), Ptr{SpglibDataset},
                            ( Ptr{Float64}, Ptr{Float64}, Ptr{Cint}, Cint, Float64 ),
                                lattice, positions, types, num_atom, symprec ))
end

function spgat_get_dataset(lattice::Array{Float64, 2},
                           positions::Array{Float64, 2},
                           types::Array{Int64, 1},
                           num_atom::Int64,
                           symprec::Float64,
                           angle_tolerance::Float64)
    num_atom = Base.convert(Cint, num_atom)
    types = Base.convert(Array{Cint, 1}, types)

    return unsafe_load(ccall( (:spgat_get_dataset, libsymspg), Ptr{SpglibDataset},
                                ( Ptr{Float64}, Ptr{Float64}, Ptr{Cint}, Cint, Float64, Float64 ),
                                    lattice, positions, types, num_atom, symprec, angle_tolerance ))
end

function spg_get_symmetry!(rotation::Array{Cint, 3},
                          translation::Array{Float64, 2},
                          max_size::Int64,
                          lattice::Array{Float64, 2},
                          position::Array{Float64, 2},
                          types::Array{Int64, 1},
                          num_atom::Int64,
                          symprec::Float64)
   #
   max_size = Base.convert(Cint, max_size)
   types = Base.convert(Array{Cint, 1}, types)
   num_atom = Base.convert(Cint, num_atom)

   return ccall((:spg_get_symmetry, libsymspg), Cint,
                (Ptr{Cint}, Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Float64}, Ptr{Cint}, Cint, Float64),
                 rotation, translation, max_size, lattice, position, types, num_atom, symprec)
end

function spgat_get_symmetry!(rotation::Array{Cint, 3},
                            translation::Array{Float64, 2},
                            max_size::Int64,
                            lattice::Array{Float64, 2},
                            position::Array{Float64, 2},
                            types::Array{Int64, 1},
                            num_atom::Int64,
                            symprec::Float64,
                            angle_tolerance::Float64)
   #
   max_size = Base.convert(Cint, max_size)
   types = Base.convert(Array{Cint, 1}, types)
   num_atom = Base.convert(Cint, num_atom)

   return ccall((:spgat_get_symmetry, libsymspg), Cint,
                (Ptr{Cint}, Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Float64}, Ptr{Cint}, Cint, Float64, Float64),
                 rotation, translation, max_size, lattice, position, types, num_atom, symprec, angle_tolerance)
end

# This method not stable at origin C-API
function spg_get_symmetry_with_collinear_spin!(rotation::Array{Cint, 3},
                                              translation::Array{Float64, 2},
                                              equivalent_atoms::Array{Cint, 1},
                                              max_size::Int64,
                                              lattice::Array{Float64, 2},
                                              position::Array{Float64, 2},
                                              types::Array{Int64, 1},
                                              spins::Array{Float64, 1},
                                              num_atom::Int64,
                                              symprec::Float64)
   #
   max_size = Base.convert(Cint, max_size)
   types = Base.convert(Array{Cint, 1}, types)
   num_atom = Base.convert(Cint, num_atom)

   return ccall((:spg_get_symmetry_with_collinear_spin, libsymspg), Cint,
                (Ptr{Cint}, Ptr{Float64}, Ptr{Cint}, Cint, Ptr{Float64}, Ptr{Float64}, Ptr{Cint}, Ptr{Float64}, Cint, Float64),
                 rotation, translation, equivalent_atoms, max_size, lattice, position, types, spins, num_atom, symprec)
end

function spgat_get_symmetry_with_collinear_spin!(rotation::Array{Cint, 3},
                                              translation::Array{Float64, 2},
                                              equivalent_atoms::Array{Cint, 1},
                                              max_size::Int64,
                                              lattice::Array{Float64, 2},
                                              position::Array{Float64, 2},
                                              types::Array{Int64, 1},
                                              spins::Array{Float64, 1},
                                              num_atom::Int64,
                                              symprec::Float64,
                                              angle_tolerance::Float64)
   #
   max_size = Base.convert(Cint, max_size)
   types = Base.convert(Array{Cint, 1}, types)
   num_atom = Base.convert(Cint, num_atom)

   return ccall((:spgat_get_symmetry_with_collinear_spin, libsymspg), Cint,
                (Ptr{Cint}, Ptr{Float64}, Ptr{Cint}, Cint, Ptr{Float64}, Ptr{Float64}, Ptr{Cint}, Ptr{Float64}, Cint, Float64, Float64),
                 rotation, translation, equivalent_atoms, max_size, lattice, position, types, spins, num_atom, symprec, angle_tolerance)
end

function spg_get_hall_number_from_symmetry(rotation::Array{Cint, 3},
                                  translation::Array{Float64, 2},
                                  num_operations::Int64,
                                  symprec::Float64)
    #
    num_operations = Base.convert(Cint, num_operations)

    return ccall((:spg_get_hall_number_from_symmetry, libsymspg), Cint,
                    (Ptr{Cint}, Ptr{Float64}, Cint, Float64),
                    rotation, translation, num_operations, symprec)
end

function spg_get_multiplicity(lattice::Array{Float64, 2},
                     position::Array{Float64, 2},
                     types::Array{Int64, 1},
                     num_atom::Int64,
                     symprec::Float64)
    #
    types = Base.convert(Array{Cint, 1}, types)
    num_atom = Base.convert(Cint, num_atom)

    return ccall((:spg_get_multiplicity, libsymspg), Cint,
                (Ptr{Float64}, Ptr{Float64}, Ptr{Cint}, Cint, Float64),
                 lattice, position, types, num_atom, symprec)
end

function spgat_get_multiplicity(lattice::Array{Float64, 2},
                     position::Array{Float64, 2},
                     types::Array{Int64, 1},
                     num_atom::Int64,
                     symprec::Float64,
                     angle_tolerance::Float64)
    #
    types = Base.convert(Array{Cint, 1}, types)
    num_atom = Base.convert(Cint, num_atom)

    return ccall((:spgat_get_multiplicity, libsymspg), Cint,
                (Ptr{Float64}, Ptr{Float64}, Ptr{Cint}, Cint, Float64, Float64),
                 lattice, position, types, num_atom, symprec, angle_tolerance)
end


text:
#!/usr/bin/env julia

"""
Usage

      ./CheckOptions.jl config_exp.toml

returns 0 if config_exp.toml is valid against harmonie.schema.json
   or a diagnostic error message if not
"""

import Harmonie, TOML

config = TOML.parsefile(ARGS[1])

Harmonie.isvalid(config) || error(Harmonie.diagnose(config))