text:
using DataFrames, ROOT, ROOTDataFrames
using Base.Test

dyndf = DataFrame(
    njet=Int64[3, 1],
    jet_pt=Vector{Float64}[Float64[1.0, 2.0, 3.0], Float64[4.0]],
    jet_eta=Vector{Float32}[Float32[3.0, 4.0, 5.0], Float32[6.0]],
    jet_id=Vector{Int64}[Int64[1,2,3], Float32[-4]],
    arr=Vector{Float64}[rand(5), rand(5)],
)

dtf = TreeDataFrame(dyndf, "dyn.root"; size_branches=Dict(:jet_pt=>:njet, :jet_eta=>:njet, :jet_id=>:njet))
@test dtf[:jet_pt] == dyndf[:jet_pt]
@test dtf[:jet_pt] == dyndf[:jet_pt]
@test dtf[:jet_id] == dyndf[:jet_id]
@test dtf[:arr] == dyndf[:arr]
print_with_color(:green, "small dynamic tree OK\n")

N = 10000
bigdf = similar(
    DataFrame(
        njet=Int64[],
        jet_pt=Vector{Float64}[],
        jet_eta=Vector{Float64}[]
    ),
    N
)

for iev=1:N
    nj = round(Int, rand()*10)
    bigdf[iev, :njet] = nj
    bigdf[iev, :jet_pt] = [rand()*200 for i=1:nj]
    bigdf[iev, :jet_eta] = [randn()*5 for i=1:nj]
end

dtf2 = TreeDataFrame(bigdf, "dyn2.root"; size_branches=Dict(:jet_pt=>:njet, :jet_eta=>:njet))
@test mean(map(sum, bigdf[:jet_pt])) == mean(map(sum, dtf2[:jet_pt]))
print_with_color(:green, "big dynamic tree OK\n")

df = TreeDataFrame([string(Pkg.dir(), "/ROOTDataFrames/dynamic.root")];
    treename="tthNtupleAnalyzer/events",
)

sumpt = 0.0
sumjet = 0
for i=1:nrow(df)
    load_row(df, i)
    n__jet = df.row.n__jet()
    jet__pt = df.row.jet__pt()[1:n__jet]
    sumjet += n__jet
    sumpt += sum(jet__pt)
end

@test_approx_eq_eps sumpt/sumjet 56.964209222875205 0.01
print_with_color(:green, "realistic dynamic tree OK\n")

r = with(df,
    row->Float64(sum(row.jet__pt()[1:row.n__jet()])),
    row->row.n__jet()>3, [:jet__pt, :n__jet], 1:length(df),
    Float64
)
@test_approx_eq_eps mean(r) 688.4442511628289 0.01
print_with_color(:green, "with() OK\n")

outdf = DataFrame()
v = Array(Vector{Float64}, 3)
v[1] = [1.0,2.0,3.0]
v[2] = [4.0,5.0]
v[3] = [6.0]

v2 = Array(Vector{Float64}, 3)
v2[1] = round(rand(2),2)
v2[2] = round(rand(4),2)
v2[3] = round(rand(6),2)

outdf[:b] = v2
outdf[:na] = [3,2,1]
outdf[:a] = v
writetree("dyntree.root", outdf)

df2 = TreeDataFrame(["dyntree.root"])
@test all(df2[:na] .== outdf[:na])

for i=1:length(df2)
    load_row(df2, i)
    @test all(df2[i, :a] .== outdf[i, :a])
end
print_with_color(:green, "dynamic tree with auto branch size OK\n")


text:
using NearestNeighbors
import NearestNeighbors: knn

"""
    knn(xgrid,ygrid::MeshArray,x,y::Array{T,1},k::Int)

Find k nearest neighbors to each point in x,y on xgrid,ygrid

```
lon=collect(0.1:0.5:2.1); lat=collect(0.1:0.5:2.1);
(f,i,j,c)=knn(Γ.XC,Γ.YC,lon,lat)
```
"""
function knn(xgrid::MeshArray,ygrid::MeshArray,
        xvec::Array{T,1},yvec::Array{T,1},k=1::Int) where {T}

        #ancillary variables
        γ=xgrid.grid
        a_f=MeshArray(γ,Int); [a_f[ii][:,:].=ii for ii=1:γ.nFaces]
        a_i=MeshArray(γ,Int); [a_i[ii]=collect(1:γ.fSize[ii][1])*ones(Int,1,γ.fSize[ii][2]) for ii=1:γ.nFaces]
        a_j=MeshArray(γ,Int); [a_j[ii]=ones(Int,γ.fSize[ii][1],1)*collect(1:γ.fSize[ii][2])' for ii=1:γ.nFaces]

        #convert to flat Array format
        a_x=write(xgrid)
        a_y=write(ygrid)
        a_f=write(a_f)
        a_i=write(a_i)
        a_j=write(a_j)

        #vector of grid points in Cartesian, 3D, coordinates
        kk=findall(isfinite.(a_x))
        x=sin.(pi/2 .-a_y[kk]*pi/180).*cos.(a_x[kk]*pi/180)
        y=sin.(pi/2 .-a_y[kk]*pi/180).*sin.(a_x[kk]*pi/180)
        z=cos.(pi/2 .-a_y[kk]*pi/180);

        #vector of target points in Cartesian, 3D, coordinates
        xx=sin.(pi/2 .-yvec*pi/180).*cos.(xvec*pi/180);
        yy=sin.(pi/2 .-yvec*pi/180).*sin.(xvec*pi/180);
        zz=cos.(pi/2 .-yvec*pi/180);

        #define tree
        kdtree = KDTree([x y z]')

        #find nearest neighbors
        idxs, _ = knn(kdtree, [xx yy zz]', k, true)
        idxs=[idxs[i][j] for i=1:length(idxs),j=1:k]

        return a_f[kk[idxs]],a_i[kk[idxs]],a_j[kk[idxs]],kk[idxs]
end

knn(xgrid::MeshArray,ygrid::MeshArray,lon::Number,lat::Number) = knn(xgrid::MeshArray,ygrid::MeshArray,[lon],[lat])

"""
    Interpolate(z_in::MeshArray,f,i,j,w)

```
using MeshArrays
lon=[i for i=-179.5:1.0:179.5, j=-89.5:1.0:89.5]
lat=[j for i=-179.5:1.0:179.5, j=-89.5:1.0:89.5]

Γ=GridLoad(GridSpec("LatLonCap",MeshArrays.GRID_LLC90); option="full")
(f,i,j,w,j_f,j_x,j_y)=InterpolationFactors(Γ,vec(lon),vec(lat))
DD=Interpolate(Γ.Depth,f,i,j,w)

using Plots
contourf(vec(lon[:,1]),vec(lat[1,:]),DD,clims=(0.,6000.))
```
"""
function Interpolate(z_in::MeshArray,f,i,j,w)
    z_out=NaN*similar(f[:,1])
    for jj=1:size(f,1)
        if !isnan(sum(w[jj,:]))
            x=[z_in[f[jj,ii]][i[jj,ii],j[jj,ii]] for ii=1:4]
            kk=findall(isfinite.(x))
            ~isempty(kk) ? z_out[jj]=sum(w[jj,kk].*x[kk])/sum(w[jj,kk]) : nothing
        end
    end
    return z_out
end

"""
    InterpolationFactors(Γ,lon::Array{T,1},lat::Array{T,1})

Compute interpolation coefficients etc from grid `Γ` to `lon,lat`

```jldoctest; output = false
using MeshArrays
γ=GridSpec("CubeSphere",MeshArrays.GRID_CS32)
Γ=GridLoad(γ; option="full")
lon=collect(45.:0.1:46.); lat=collect(60.:0.1:61.)
(f,i,j,w,j_f,j_x,j_y)=InterpolationFactors(Γ,lon,lat)
YC=Interpolate(Γ.YC,f,i,j,w)
extrema(i)==(9,10)

# output

true
```
"""
function InterpolationFactors(Γ,lon::Array{T,1},lat::Array{T,1}) where {T}
#to match gcmfaces test case (`interp=gcmfaces_interp_coeffs(0.1,0.1);`)
#set: iiTile=17; XC0=6.5000; YC0=-0.1994

        #main loop
        i_f=fill(0,length(lon),4)
        i_i=fill(0,length(lon),4)
        i_j=fill(0,length(lon),4)
        i_w=fill(NaN,length(lon),4)
        j_f=fill(0,length(lon),1)
        j_x=fill(0.0,length(lon),1)
        j_y=fill(0.0,length(lon),1)

        #ancillary variables
        (f,i,j,c)=knn(Γ.XC,Γ.YC,lon,lat)

        #1. tiles and τ
        fs=Γ.XC.fSize
        s=fill(0,2*length(fs))
        [s[collect(1:2) .+ (i-1)*2]=collect(fs[i]) for i in 1:length(fs)]
        ni=gcd(s); nj=gcd(s); γ=Γ.XC.grid
        τ=Tiles(γ,ni,nj); tiles=MeshArray(γ,Int);
        [tiles[τ[ii].face][τ[ii].i,τ[ii].j].=ii for ii in 1:length(τ)]

        #2. t_XC, t_XC, t_f, t_i, t_j
        t=vec(write(tiles)[c])
        t_list=unique(t)
        t_XC=Tiles(τ,exchange(Γ.XC))
        t_YC=Tiles(τ,exchange(Γ.YC))

        t_f=MeshArray(γ,Int); [t_f[ii][:,:].=ii for ii=1:γ.nFaces]
        t_i=MeshArray(γ,Int); [t_i[ii]=collect(1:γ.fSize[ii][1])*ones(Int,1,γ.fSize[ii][2]) for ii=1:γ.nFaces]
        t_j=MeshArray(γ,Int); [t_j[ii]=ones(Int,γ.fSize[ii][1],1)*collect(1:γ.fSize[ii][2])' for ii=1:γ.nFaces]

        t_f=Tiles(τ,exchange(t_f))
        t_i=Tiles(τ,exchange(t_i))
        t_j=Tiles(τ,exchange(t_j))

        x_q=fill(0.0,1,4)
        y_q=fill(0.0,1,4)
        tmpx=fill(0.0,1,4)
        tmpy=fill(0.0,1,4)
        w=fill(0.0,1,4)

        #main loop
        for ii=1:length(t_list)
                tt=t_list[ii]

                ff=τ[tt].face
                ii0=minimum(τ[tt].i)+Int(ni/2)
                jj0=minimum(τ[tt].j)+Int(nj/2)
                XC0=Γ.XG.f[ff][ii0,jj0]
                YC0=Γ.YG.f[ff][ii0,jj0]
                #
                (x_grid,y_grid)=StereographicProjection(XC0,YC0,t_XC[tt],t_YC[tt])
                (x_quad,y_quad,i_quad,j_quad)=QuadArrays(x_grid,y_grid)
                #
                x=minimum(τ[tt].i)-0.5 .+collect(-1:ni)*ones(Int,1,nj+2)
                y=minimum(τ[tt].j)-0.5 .+ones(Int,ni+2,1)*collect(-1:nj)'
                #
                angsum=fill(0.0,size(x_quad,1))

                for pp in findall(t.==tt)
                        (x_trgt,y_trgt)=StereographicProjection(XC0,YC0,lon[pp],lat[pp])
                        PolygonAngle(x_quad,y_quad,x_trgt,y_trgt,angsum)
                        #angsum=PolygonAngle_deprecated(x_quad,y_quad,x_trgt,y_trgt)
                        if sum(angsum.>180.)>0
                                kk=findall(angsum.>180.)[end]
                                i_f[pp,:].=[t_f[tt][i_quad[kk,i]+1,j_quad[kk,i]+1] for i=1:4]
                                i_i[pp,:].=[t_i[tt][i_quad[kk,i]+1,j_quad[kk,i]+1] for i=1:4]
                                i_j[pp,:].=[t_j[tt][i_quad[kk,i]+1,j_quad[kk,i]+1] for i=1:4]
                                x_q[:].=x_quad[kk,:]
                                y_q[:].=y_quad[kk,:]
                                QuadCoeffs(x_q,y_q,x_trgt,y_trgt,w)
                                i_w[pp,:].=w[:]
                                #
                                [tmpx[i]=x[i_quad[kk,i]+1,j_quad[kk,i]+1] for i=1:4]
                                [tmpy[i]=y[i_quad[kk,i]+1,j_quad[kk,i]+1] for i=1:4]
                                #
                                j_f[pp]=ff
                                j_x[pp]=sum(tmpx[:].*i_w[pp,:])
                                j_y[pp]=sum(tmpy[:].*i_w[pp,:])
                        end
                end
        end

        return i_f,i_i,i_j,i_w,j_f,j_x,j_y
end

InterpolationFactors(Γ,lon::Number,lat::Number) = InterpolationFactors(Γ,[lon],[lat])

"""
    StereographicProjection(XC0,YC0,XC,YC)

Apply stereographic projection that puts `XC0,YC0` at `0.0,0.0`
to target point(s) `XC,YC`

```
lon=collect(45.:0.1:46.); lat=collect(60.:0.1:61.)
x,y=StereographicProjection(45.,60.,lon,lat)
```
"""
function StereographicProjection(XC0::Number,YC0::Number,XC,YC)
        #compute spherical coordinates:
        phi=XC; theta=90 .-YC;
        phi0=XC0; theta0=90-YC0;

        #compute cartesian coordinates:
        X=sind.(theta).*cosd.(phi);
        Y=sind.(theta).*sind.(phi);
        Z=cosd.(theta);

        x=X; y=Y; z=Z;

        #bring chosen point to the north pole:
        xx=x; yy=y; zz=z;
        x=cosd(phi0).*xx+sind(phi0).*yy;
        y=-sind(phi0).*xx+cosd(phi0).*yy;
        z=zz;

        xx=x; yy=y; zz=z;
        x=cosd(theta0)*xx-sind(theta0)*zz;
        y=yy;
        z=sind(theta0)*xx+cosd(theta0)*zz;

        #stereographic projection from the south pole:
        xx=x./(1 .+z);
        yy=y./(1 .+z);

        #nrm=sqrt(xx.^2+yy.^2);
        #msk=1+0*nrm; msk(nrm>tan(pi/4/2))=NaN;%mask points outside of pi/4 cone

        return xx,yy
end

"""
    PolygonAngle_deprecated(px,py,x=[],y=[])

Compute sum of interior angles for polygons or points-to-polygons (when
`px,py,x,y` is provided as input). `px,py` are `MxN` matrices where each line
specifies one polygon. (optional) `x,y` are position vectors.
"""
function PolygonAngle_deprecated(px,py,x=[],y=[])

        M=size(px,1); N=size(px,2); P=1;
        doPointsInPolygon=false
        if length(x)>0;
                doPointsInPolygon=true
                x=reshape(x,1,length(x))
                y=reshape(y,1,length(y))
                P=length(x)
        end;

        angsum=zeros(M,P)
        for ii=0:N-1
                ppx=circshift(px,[0 -ii])
                ppy=circshift(py,[0 -ii])
                if ~doPointsInPolygon
                        #compute sum of corner angles
                        v1x=ppx[:,2]-ppx[:,1]
                        v1y=ppy[:,2]-ppy[:,1]
                        v2x=ppx[:,4]-ppx[:,1]
                        v2y=ppy[:,4]-ppy[:,1]
                else;
                        #compute sum of sector angles
                        v1x=ppx[:,1]*ones(1,P)-ones(M,1)*x
                        v1y=ppy[:,1]*ones(1,P)-ones(M,1)*y
                        v2x=ppx[:,2]*ones(1,P)-ones(M,1)*x
                        v2y=ppy[:,2]*ones(1,P)-ones(M,1)*y
                end
                tmp=( v1x.*v2x+v1y.*v2y )./sqrt.( v1x.*v1x+v1y.*v1y )./sqrt.( v2x.*v2x+v2y.*v2y )
                g_acos=acos.( min.(max.(tmp,-1.0),1.0) )
                g_sin= ( v1x.*v2y-v1y.*v2x )./sqrt.( v1x.*v1x+v1y.*v1y )./sqrt.( v2x.*v2x+v2y.*v2y )
                angsum=angsum+rad2deg.(g_acos.*sign.(g_sin));
        end;

        return angsum
end

"""
    PolygonAngle(px::Array,py::Array,angsum::Array)
"""
function PolygonAngle(px::Array,py::Array,angsum::Array)
        M=size(px,1)
        N=size(px,2)
        angsum .= 0.0
        for ii=0:N-1
                i1=mod1(ii+1,N)
                i2=mod1(ii+2,N)
                i4=mod1(ii+4,N)
                for jj=1:M
                #compute sum of sector angles
                v1x=px[jj,i2]-px[jj,i1]
                v1y=py[jj,i2]-py[jj,i1]
                v2x=px[jj,i4]-px[jj,i1]
                v2y=py[jj,i4]-py[jj,i1]

                tmp=( v1x.*v2x+v1y.*v2y )./sqrt.( v1x.*v1x+v1y.*v1y )./sqrt.( v2x.*v2x+v2y.*v2y )
                g_acos=acos.( min.(max.(tmp,-1.0),1.0) )
                g_sin= ( v1x.*v2y-v1y.*v2x )./sqrt.( v1x.*v1x+v1y.*v1y )./sqrt.( v2x.*v2x+v2y.*v2y )
                angsum[jj] += rad2deg(g_acos*sign(g_sin))
                end
        end
end

"""
    function PolygonAngle(px::Array,py::Array,x::Array,y::Array,angsum)

Compute sum of interior angles (`angsum`) for polygons or points-to-polygons (when
`px,py,x,y,angsum` is provided as input). `px,py` are `MxN` matrices where each line
specifies one polygon. (optional) `x,y` are position vectors.

```jldoctest; output = false
using MeshArrays
px=[0. 0. 1. 1.]; py=[0. 1. 1. 0.];
x=collect(-1.0:0.25:2.0); y=x;
angsum=fill(0.0,1,length(x))
MeshArrays.PolygonAngle(px,py,x,y,angsum)

isa(angsum,Array)

# output

true
```

"""
function PolygonAngle(px::Array,py::Array,x::Array,y::Array,angsum)
        for ii in 1:length(x)
                PolygonAngle(px,py,x[ii],y[ii],view(angsum,:,ii))
        end
end

"""
    PolygonAngle(px::Array,py::Array,x::Number,y::Number,angsum)
"""
function PolygonAngle(px::Array,py::Array,x::Number,y::Number,angsum)
        M=size(px,1)
        N=size(px,2)
        angsum .= 0.0
        for ii=0:N-1
                i1=mod1(ii+1,N)
                i2=mod1(ii+2,N)
                for jj=1:M
                #compute sum of sector angles
                v1x=px[jj,i1] -x
                v1y=py[jj,i1] -y
                v2x=px[jj,i2] -x
                v2y=py[jj,i2] -y

                tmp=( v1x.*v2x+v1y.*v2y )./sqrt.( v1x.*v1x+v1y.*v1y )./sqrt.( v2x.*v2x+v2y.*v2y )
                g_acos=acos.( min.(max.(tmp,-1.0),1.0) )
                g_sin= ( v1x.*v2y-v1y.*v2x )./sqrt.( v1x.*v1x+v1y.*v1y )./sqrt.( v2x.*v2x+v2y.*v2y )
                angsum[jj] += rad2deg(g_acos*sign(g_sin))
                end
        end
end

"""
    QuadArrays(x_grid,y_grid)

Transform x_grid,y_grid (size ni+2,nj+2) into x_quad,y_quad,i_quad,j_quad
quadrilaterals (size ni+1*nj+1,4) where i_quad,j_quad are point indices
"""
function QuadArrays(x_grid::Array{T,2},y_grid::Array{T,2}) where {T}
        ni,nj=size(x_grid) .-2

        x_quad=Array{Float64,2}(undef,(ni+1)*(nj+1),4)
        y_quad=Array{Float64,2}(undef,(ni+1)*(nj+1),4)
        i_quad=Array{Int64,2}(undef,(ni+1)*(nj+1),4)
        j_quad=Array{Int64,2}(undef,(ni+1)*(nj+1),4)

        didj=[[0 0];[1 0];[1 1];[0 1]]
        for pp=1:4
                di=didj[pp,1]
                dj=didj[pp,2]

                #note the shift in indices due to exchange above
                tmp=x_grid[1+di:ni+1+di,1+dj:nj+1+dj]
                x_quad[:,pp]=vec(tmp)
                tmp=y_grid[1+di:ni+1+di,1+dj:nj+1+dj]
                y_quad[:,pp]=vec(tmp)

                tmp=collect(0+di:ni+di)*ones(1,nj+1)
                i_quad[:,pp]=vec(tmp)
                tmp=ones(ni+1,1)*transpose(collect(0+dj:nj+dj));
                j_quad[:,pp]=vec(tmp)
        end

        return x_quad,y_quad,i_quad,j_quad
end

"""
    QuadCoeffs(px,py,ox=[],oy=[],ow=[])

Compute bilinear interpolation coefficients for `ox,oy` within `px,py`
by remapping these quadrilaterals to the `unit square`.
- `px,py` are `Mx4` matrices where each line specifies one quadrilateral.
- `ox,oy` are `MxP` position matrices
- `ow` (output) are the `MxPx4` bilinear interpolation weights

```
QuadCoeffs([-1., 8., 13., -4.]',[-1., 3., 11., 8.]',0.,6.)
QuadCoeffs([0., 2., 3., 1.]',[0., 0., 1., 1.]',0.1,0.1)
```
"""
function QuadCoeffs(px,py,ox=[],oy=[],ow=[])
        #test case from https://www.particleincell.com/2012/quad-interpolation/
        #  QuadCoeffs([-1., 8., 13., -4.]',[-1., 3., 11., 8.]',0.,6.)
        #However the case of a perfect parallelogram needs special treatment
        #  QuadCoeffs([0., 2., 3., 1.]',[0., 0., 1., 1.]',0.1,0.1)
        #Deals with this situtation by falling back to ParaCoeffs
        #  ParaCoeffs([0., 2., 3., 1.]',[0., 0., 1., 1.]',0.1,0.1)

        #1. solve linear problem (`a,b` vectors from `px,py`)
        #  A=[1 0 0 0;1 1 0 0;1 1 1 1;1 0 1 0]; AI = inv(A);
        #  AI=[1 0 0 0;-1 1 0 0;-1 0 0 1; 1 -1 1 -1];
        #  a = AI*px';
        #  b = AI*py';
        #This defines the mapping from logical `l,m` to physical `x,y` as
        #  x=a(1)+a(2)*l+a(3)*m+a(2)*l*m;
        #  y=b(1)+b(2)*l+b(3)*m+b(2)*l*m;

        a=[px[1] -px[1]+px[2] -px[1]+px[4] px[1]-px[2]+px[3]-px[4]]
        a[findall(abs.(a).<1e-8)].=0.0

        b=[py[1] -py[1]+py[2] -py[1]+py[4] py[1]-py[2]+py[3]-py[4]]
        b[findall(abs.(b).<1e-8)].=0.0

        #2. select between the two solutions (to 2nd order
        #non-linear problem below) using polygon interior angles
        angsum=fill(0.0,1)
        PolygonAngle(px,py,angsum)
        sgn=fill(NaN,1)
        isapprox(angsum[1],360.0,rtol=0.01) ? sgn[1]=1.0 : nothing
        isapprox(angsum[1],-360.0,rtol=0.01) ? sgn[1]=-1.0 : nothing
        #ii=findall(isnan.(angsum))
        #length(ii)>0 ? println("warning: edge point was found") : nothing

        #3. solve non-linear problem for `pl,pm` from `px,py` & `a,b`
        # This defines the mapping from physical `x,y` to logical `l,m`

        a=reshape(a,(size(a,1),1,size(a,2)));
        b=reshape(b,(size(b,1),1,size(b,2)));

        # quadratic equation coeffs, `aa*mm^2+bb*m+cc=0`
        if ~isempty(ox);
                x=ox; y=oy;
        else;
                x=px; y=py;
                a=repeat(a,1,size(x,2),1);
                b=repeat(b,1,size(x,2),1);
                sgn=repeat(sgn,1,size(x,2));
        end;
        #
        det=fill(0.0,size(x))
        pm=fill(0.0,size(x))
        pl=fill(0.0,size(x))
        for ii=1:size(x,1), jj=1:size(x,2)
                aa = a[ii,jj,4]*b[ii,jj,3]-a[ii,jj,3]*b[ii,jj,4]
                bb = a[ii,jj,4]*b[ii,jj,1]-a[ii,jj,1]*b[ii,jj,4]+a[ii,jj,2]*b[ii,jj,3]-a[ii,jj,3]*b[ii,jj,2]+x[ii,jj]*b[ii,jj,4]-y[ii,jj].*a[ii,jj,4]
                cc = a[ii,jj,2]*b[ii,jj,1]-a[ii,jj,1]*b[ii,jj,2]+x[ii,jj]*b[ii,jj,2]-y[ii,jj].*a[ii,jj,2]
                #compute `pm = (-b+sqrt(b^2-4ac))/(2a)`
                det[ii,jj] = sqrt(bb*bb - 4.0*aa*cc)
                pm[ii,jj]  = (-bb+sgn[ii,jj]*det[ii,jj])/(2.0*aa)
                #compute `pl` by substitution in equation system
                pl[ii,jj]  = (x[ii,jj]-a[ii,jj,1]-a[ii,jj,3]*pm[ii,jj])/(a[ii,jj,2]+a[ii,jj,4]*pm[ii,jj])
        end

        if ~isempty(ox);
                tmp1=(1 .-pl).*(1 .-pm)
                tmp2=pl.*(1 .-pm)
                tmp3=pl.*pm
                tmp4=(1 .-pl).*pm
                ow[:]=cat(tmp1,tmp2,tmp3,tmp4; dims=3)
                #ow[:].=cat(tmp1,tmp2,tmp3,tmp4; dims=3)[:]
                #treat pathological cases if needed
                tmp=ParaCoeffs(px,py,[ox],[oy])
                for kk=1:length(ow)
                        !isfinite(ow[kk]) ? ow[kk]=tmp[kk] : nothing
                end
        end

end

"""
    ParaCoeffs(px,py,ox=[],oy=[])

Compute bilinear interpolation coefficients for `ox,oy` within `px,py`
by remapping these parallelograms to the `unit square`.
- `px,py` are `Mx4` matrices where each line specifies one quadrilateral.
- `ox,oy` are `MxP` position matrices
- `pw` (output) are the `MxPx4` bilinear interpolation weights

```
x=1.0; y=1.0 #Try send the corners to unit square corners?
println(vec(ParaCoeffs([0., 2., 3., 1.]',[0., 0., 1., 1.]',x,y)))
println(vec(QuadCoeffs([0., 2.01, 3., 1.]',[0., 0., 1., 1.]',x,y)))
```
"""
function ParaCoeffs(px,py,ox=[],oy=[])

        tmp1=px[:,1];
        tmp2=-px[:,1]+px[:,2];
        tmp3=-px[:,2]+px[:,3];
        a=[tmp1 tmp2 tmp3];

        tmp1=py[:,1];
        tmp2=-py[:,1]+py[:,2];
        tmp3=-py[:,2]+py[:,3];
        b=[tmp1 tmp2 tmp3];

#        (m,l)=inv([a[1,2] a[1,3];b[1,2] b[1,3]])*[ox[1]-a[1,1]; oy[1]-b[1,1]]
        m=( b[:,3].*(ox-a[:,1])-a[:,3].*(oy-b[:,1]) ) ./(a[:,2].*b[:,3]-a[:,3].*b[:,2])
        l=( -b[:,2].*(ox-a[:,1])+a[:,2].*(oy-b[:,1]) ) ./(a[:,2].*b[:,3]-a[:,3].*b[:,2])

        ow=[];
        tmp1=(1 .-l).*(1 .-m)
        tmp4=l.*(1 .-m)
        tmp3=l.*m
        tmp2=(1 .-l).*m
        ow=cat(tmp1,tmp2,tmp3,tmp4; dims=3)

        return ow
end

ext: jl


text:
include("../examples/rijke.jl")
using OrdinaryDiffEq
using DiffEqSensitivity, Zygote
using JLD
using Distributed
#using PyPlot
using SharedArrays
function ad_solve(u0, eps, n, v, ind)
	d = N
	s = [7.0, 0.2]
	un = u0 + eps*v
	un = Rijke_ODE(un, s, n)
	return un[ind]
end
function ad_adj_solve(u0, n, v)
	d = N
	s = [7.0, 0.2]
	un = Rijke_ODE(u0, s, n)
	return dot(un, v)
end

function compute_perturbations()
d = N
u = rand(d)
s = [7.0, 0.2]
nSteps = 2500
ad_norm = zeros(nSteps)
ad_adj_norm = zeros(nSteps)
ad = SharedArray{Float64}(d)
ad_adj = zeros(d)
eps_ad = 0.
v_norm, w_norm, fd_norm = zeros(nSteps),
				zeros(nSteps), zeros(nSteps)
v, w, fd = rand(d), rand(d),
				rand(d)
u_trj = zeros(d, nSteps)
eps = 1.e-4
u_p = u + eps*fd
for n = 1:nSteps
	ans = @distributed for j = 1:d
			ad[j] = Zygote.gradient(eps -> ad_solve(u, eps, 1, v, j), eps)[1]
	end
	wait(ans)
	ad_norm[n] = norm(ad)

	du = dRijke(u, s, 1.e-6)
	v = du*v
	v_norm[n] = norm(v)

	u_trj[:,n] = u
	u = Rijke_ODE(u, s, 1)
	u_p = Rijke_ODE(u_p, s, 1)

	fd = (u_p - u)./eps
	fd_norm[n] = norm(fd)

end
for n = nSteps:-1:1
	u = u_trj[:,n]
	ad_adj .= Zygote.gradient(u -> ad_adj_solve(u, 1, w), u)[1]
	ad_adj_norm[n] = norm(ad_adj)
	du = dRijke(u, s, 1.e-6)
	du = du'
	w = du*w
	w_norm[n] = norm(w)
end
save("../data/rijke_perturbations/v_w_fd_ada_norms.jld",
	 "v_norm", v_norm,
	 "w_norm", w_norm,
	 "fd_norm", fd_norm,
	 "ad_norm", ad_norm,
	 "ad_adj_norm", ad_adj_norm)
end
function plot_perturbations()
	X = load("../data/rijke_perturbations/v_w_fd_ada_norms.jld")
	v_norm = X["v_norm"]
	w_norm = X["w_norm"]
	fd_norm = X["fd_norm"]
	#X = load("../data/rijke_perturbations/ad_norms.jld")
	ad_norm = X["ad_norm"]
	ad_adj_norm = X["ad_adj_norm"]
	@show maximum(ad_adj_norm), maximum(w_norm)
	nSteps = div(size(v_norm)[1],1)
	v_norm = v_norm[1:nSteps]
	fd_norm = fd_norm[1:nSteps]
	w_norm = w_norm[1:nSteps]
	ad_norm = ad_norm[1:nSteps]
	fig, ax = subplots(1,1)
	ax.semilogy(dt*(1:nSteps), v_norm, "v", ms=15.0,
				label="tangent")
	ax.semilogy(dt*(1:nSteps), w_norm, ".", ms=20.0,
				label="adjoint")
	ax.semilogy(dt*(1:nSteps), fd_norm, "1", ms=4.0,
				label="FD")
	ax.semilogy(dt*(1:nSteps), ad_norm, "^", ms=4.0,
				label="forward AD")
	ax.semilogy(dt*(1:nSteps), ad_adj_norm, "P", ms=4.0,
				label="reverse AD")

	ax.xaxis.set_tick_params(labelsize=28)
	ax.yaxis.set_tick_params(labelsize=28)
	ax.set_xlabel("time", fontsize=28)
	ax.set_ylabel("Perturbation norms", fontsize=28)
	ax.grid(true)
	lgnd = fig.legend(loc=(0.3,0.6),
					  fontsize=28)
	lgnd.legendHandles[1]._legmarker.set_markersize(20)
	lgnd.legendHandles[2]._legmarker.set_markersize(20)
	lgnd.legendHandles[3]._legmarker.set_markersize(20)
	lgnd.legendHandles[4]._legmarker.set_markersize(20)
	lgnd.legendHandles[5]._legmarker.set_markersize(20)
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
module  ProximalFunctionTest

using Test
using ProximalBase
using Random
using LinearAlgebra
using Distributions
import JuMP, Ipopt

Random.seed!(1)

@testset "prox zero" begin
  f = ProxZero()

  p = 10
  x = randn(p)
  z = zeros(p)

  @test prox!(f, z, x) == x
  @test prox(f, x) == z
end

@testset "proximal_l1" begin
  @testset "shrink to zero" begin
    x = randn(10)
    lambda = maximum(abs, x) + 0.1
    g = ProxL1(lambda)

    hat_x = randn(10)
    prox!(g, hat_x, x)
    @test maximum(abs, hat_x) == 0.

    @test value(g, x) ≈ lambda * sum(abs, x)
  end

  @testset "shrink towards zero" begin
    x = [1.0, 1.4, -3.2]
    lambda = 0.1
    g = ProxL1(lambda)

    hat_x = randn(3)
    prox!(g, hat_x, x)
    @test  hat_x ≈ [0.9, 1.3, -3.1]

    @test value(g, x) ≈ lambda * sum(abs, x)
    @test value(g, hat_x) ≈ lambda * sum(abs, hat_x)
  end

  @testset "shrink with vector" begin
    x = [1.0, 1.4, -3.2]
    lambda = [0.1, 1.5, 1]
    g = ProxL1(1., lambda)

    hat_x = randn(3)
    prox!(g, hat_x, x)
    @test hat_x ≈ [0.9, 0., -2.2]
    @test prox(g, x) ≈ [0.9, 0., -2.2]

    @test value(g, x) ≈ dot(lambda, abs.(x))
    @test value(g, hat_x) ≈ dot(lambda, abs.(hat_x))
  end

end

@testset "proximal_l1_fused" begin

    @testset "random" begin

    m = JuMP.Model(JuMP.with_optimizer(Ipopt.Optimizer, print_level=0))

    JuMP.@variable(m, z1)
    JuMP.@variable(m, z2)
    JuMP.@variable(m, t[1:3] >= 0)
    JuMP.@constraint(m, z1 <= t[1] )
    JuMP.@constraint(m, -t[1] <= z1 )
    JuMP.@constraint(m, z2 <= t[2] )
    JuMP.@constraint(m, -t[2] <= z2 )
    JuMP.@constraint(m, z1 - z2 <= t[3])
    JuMP.@constraint(m, -t[3] <= z1 - z2 )

    for i=1:1000
      x1 = randn()
      x2 = randn()
      λ1 = rand(Uniform(0.,1.))
      λ2 = rand(Uniform(0.,1.))

      JuMP.@objective(m, Min, ((x1-z1)^2+(x2-z2)^2)/2. + λ1 * (t[1]+t[2]) + λ2 * t[3])
      JuMP.optimize!(m)

      zp1, zp2 = ProximalBase.proxL1Fused(x1, x2, λ1, λ2)

      @test abs(JuMP.value.(z1) - zp1) + abs(JuMP.value.(z2) - zp2)  ≈ 0. atol=2e-4
    end

    m = JuMP.Model(JuMP.with_optimizer(Ipopt.Optimizer, print_level=0))
    JuMP.@variable(m, z1)
    JuMP.@variable(m, z2)
    JuMP.@variable(m, t >= 0)
    JuMP.@constraint(m, z1 - z2 <= t)
    JuMP.@constraint(m, -t <= z1 - z2 )

    for i=1:1000
      x1 = randn()
      x2 = randn()
      λ = rand(Uniform(0.,1.))

      JuMP.@objective(m, Min, ((x1-z1)^2+(x2-z2)^2)/2. + λ * t)
      JuMP.optimize!(m)

      zp1, zp2 = ProximalBase.proxL1Fused(x1, x2, 0., λ)

      @test abs(JuMP.value.(z1) - zp1) + abs(JuMP.value.(z2) - zp2)  ≈ 0. atol=1e-4
    end
  end # random testset

  @testset "nonrandom" begin
      x1 = randn()
      x2 = randn()
      λ1 = rand(Uniform(0.,1.))

      zp1, zp2 = ProximalBase.proxL1Fused(x1, x2, λ1, 0.)
      @test ProximalBase.shrink(x1, λ1) ≈ zp1
      @test ProximalBase.shrink(x2, λ1) ≈ zp2

      λ1 = 1000.
      zp1, zp2 = ProximalBase.proxL1Fused(x1, x2, λ1, 0.)
      @test ProximalBase.shrink(x1, λ1) ≈ 0.
      @test ProximalBase.shrink(x2, λ1) ≈ 0.
  end
end # proximal_l1

@testset "proximal_l2" begin

  @testset "shrink to zero" begin
    p = 10
    x = randn(p)
    hat_x = randn(p)

    lambda = norm(x) + 0.1
    g = ProxL2(lambda)
    # check the norm
    @test value(g, x) ≈ lambda * norm(x)
    ln = 0.
    for j=1:p
      ln += x[j]^2
    end
    @test value(g, x) ≈ lambda * sqrt(ln)

    # check shrinkage
    prox!(g, hat_x, x, 1.)
    @test hat_x ≈ zeros(Float64, p)
    @test norm(hat_x) ≈ 0.
  end

  @testset "shrink not to zero" begin
    p = 2
    x = [1., 2.]
    hat_x = randn(p)

    lambda = 1.
    # check the norm
    g = ProxL2(lambda)

    @test value(g, x) ≈ lambda * norm(x)
    ln = 0.
    for j=1:p
      ln += x[j]^2
    end
    @test value(g, x) ≈ lambda * sqrt(ln)

    # check shrinkage
    prox!(g, hat_x, x, 1.)

    @test hat_x ≈ (1. - 1. / sqrt(5.))*x
    @test norm(hat_x) ≈ sqrt(5.) - 1.
  end

end

@testset "proximal_l2sq" begin

  @testset "prox operator" begin
    x = [1., 2.]
    lambda = 1.

    # norm value
    g = ProxL2Sq(0.)
    @test value(g, x) == 0.

    g = ProxL2Sq(1.)
    @test value(g, x) == sum(abs2, x)

    g = ProxL2Sq(2.)
    @test value(g, x) ≈ 2. * norm(x)^2

    # prox
    hat_x = similar(x)

    g = ProxL2Sq(0.)
    @test prox!(g, hat_x, x) ≈ x

    g = ProxL2Sq(1.)
    @test prox!(g, hat_x, x) ≈ x / 3.

    g = ProxL2Sq(2.)
    @test prox!(g, hat_x, x) ≈ x / 5.

  end

end

@testset "proximal_nuclear" begin

  y = randn(30, 10)
  y = y' * y

  U, S, V = svd(y)
  for k=1:10
    S[k] = max(S[k] - 0.1, 0.)
  end
  tx = U * (Diagonal(S)*V')

  x = randn(10, 10)
  g = ProxNuclear(0.1)
  prox!(g, x, y)

  @test x ≈ tx
end

@testset "proximal_l1l2" begin

  @testset "shrink to zero" begin
    p = 10
    x = AtomIterate(p, 2)
    hat_x = AtomIterate(p, 2)

    x .= randn(p)
    hat_x .= randn(p)

    lambda = norm(x) + 0.1

    # check the norm
    g = ProxL2(lambda, [1., 1.])
    @test value(g, x) ≈ lambda * norm(x[1:5]) + lambda * norm(x[6:10])

    # check shrinkage
    prox!(g, hat_x, x, 1.)

    @test hat_x.storage ≈ zeros(Float64, p)
    @test norm(hat_x) ≈ 0.
  end

  @testset "shrink not to zero" begin
    p = 10
    x = AtomIterate(p, 2)
    hat_x = AtomIterate(p, 2)

    x .= randn(p)
    hat_x .= randn(p)

    lambda = min(norm(x[1:5]), norm(x[6:10])) - 0.2

    # check the norm
    g = ProxL2(lambda, [1., 1.])
    @test value(g, x) ≈ lambda * norm(x[1:5]) + lambda * norm(x[6:10])

    # check shrinkage
    prox!(g, hat_x, x, 1.)

    @test hat_x.storage[1:5] ≈  (1. - lambda/norm(x[1:5])) * x[1:5]
    @test hat_x.storage[6:10] ≈ (1. - lambda/norm(x[6:10])) * x[6:10]
  end

end

@testset "proximal_logdet" begin
  p = 10

  out = zeros(p,p)

  # identity
  Σ = Symmetric(Matrix{Float64}(I, p, p))
  V = Matrix{Float64}(I, p, p)
  g = ProxGaussLikelihood(Σ)
  @test prox(g, V) ≈ Matrix{Float64}(I, p, p)
  @test prox!(g, out, V) ≈ Matrix{Float64}(I, p, p)
  @test value(g, V) ≈ p

  γ = 0.5
  ρ = 2.
  sol = Matrix{Float64}(I, p, p) * (-(1. - ρ)+sqrt((1. - ρ)^2 + 4. * ρ)) / (2. * ρ)
  @test prox!(g, out, V, γ) ≈ sol

end

end


text:
using Literate
using Weave

# https://github.com/baggepinnen/LiterateWeave.jl/
function literateweave(source, args...; kwargs...)
    tmpname = tempname()
    Literate.markdown(source, tmpname, documenter=false)
    if source[end-1:end] == "jl"
      sourcename = source[1:end-2] * "md"
    else
      @error "Need .jl file!"
    end
    sourcename = match(r"(\w+.md)", sourcename)[1]
    sourcename = joinpath(tmpname,sourcename)
    jmdsource = replace(sourcename,".md"=>".jmd")
    run(`cp $(sourcename) $(jmdsource)`)
    weave(jmdsource, args...; kwargs...)
end

tmp_dir = @__DIR__
#Literate.markdown(joinpath(tmp_dir,"linear_estimation.jl"), tmp_dir,
#                   documenter=false)

literateweave(joinpath(tmp_dir,"linear_estimation.jl"), out_path=tmp_dir)
literateweave(joinpath(tmp_dir,"data_analysis.jl"), out_path=tmp_dir)


text:
using Dates, HTTP
"""
    latest_beagle()
---
Return the latest beagle URL
"""
function latest_beagle()
    function later(a, b)
        mon = Dict("Jan" => 1, "Feb" => 2, "Mar" => 3,
                   "Apr" => 4, "May" => 5, "Jun" => 6,
                   "Jul" => 7, "Aug" => 8, "Sep" => 9,
                   "Oct" => 10, "Nov" => 11, "Dec" => 12, "mmm" => 0)
        da, db = a[8:9], b[8:9]
        ma, mb = mon[a[10:12]], mon[b[10:12]]
        ya, yb = a[13:14], b[13:14]
        if ya == yb
            if ma == mb
                return da > db ? a : b
            else
                return ma > mb ? a : b
            end
        else
            return ya > yb ? a : b
        end
    end

    BeagleURL = "https://faculty.washington.edu/browning/beagle/"
    # Note Beagle is not version controlled by git.
    # It seems that names for new versions are following the pattern like
    # beagle.ddMmmyy.xxx.jar
    # so I just grab all the file names in $BeagleURL
    # There also must be better way to get the URL.
    beagle = "beagle.00mmm00.xxx.jar"
    web = HTTP.request("GET", BeagleURL)
    txt = split(String(web.body), "\n")
    for line in txt
        r = match(r"beagle\.\d{2}\w{3}\d{2}.{5}jar", line)
        r == nothing || (beagle = later(beagle, r.match))
    end
    message("\tUpdating Beagle to $beagle.")
    return joinpath(BeagleURL, beagle)
end

"""
    latest_plink()
---
Return URL of the latest plink version 1 for Linux x86_64.
"""
function update_plink()
    message("\tUpdating plink to the latest")
    plinkURL = "http://s3.amazonaws.com/plink1-assets"
    latest = "plink_linux_x86_64_latest.zip"
    odir = pwd()
    cd(bin_dir)
    download(joinpath(plinkURL, latest), latest)
    run(`unzip -u $latest`)
    cd(odir)
    # using LightXML
    # Below uses "LightXML", may serve for future usage.
    # fxml = download("plinkURL")  # A tmp file store the XML result
    # xdoc = parse_file(fxml)
    # xroot = root(xdoc)
    # cnt = get_elements_by_tagname(xroot, "Contents")
    # for c in cnt
    #     t = find_element(c, "Key")
    #     f = content(t)
    # end
end

"""
    make_bins()
---
Compile Xijiang's C++ codes into binaries.
"""
function make_bins()
    bins = ["mrg2bgl", "raw2gt"]
    for b in bins
        if !isfile(joinpath(bin_dir, b))
            print(lpad("g++ -O3 -Wall -std=c++17 $b.cpp -o $b", 60))
            run(`g++ -O3 -Wall -std=c++17 -o $bin_dir/$b $cpp_dir/$b.cpp`)
            done("Compiled")
        end
    end
end

"""
    Update()
---
Every time one starts this package, the package will automatically check if `plink` and `beagle.jar` exist.
If not, the package will download the latest version of these two files.
The package will update, or download them on 17th every month anyway.
"""
function Update(force = false)
    title("Initialization")
    # Welcome and copyright message
    msg = begin
        w = 60
        a = "NSG Julia Package"
        b = repeat('=', length(a)+2)
        c = "Developed by Xijiang Yu @NMBU"
        d = "Copyright © 2020"
        join([lpad(a, w), lpad(b, w+1), lpad(c, w), lpad(d, w)], '\n')
    end
    message(msg)

    isdir(bin_dir) || mkdir(bin_dir)
    if Dates.day(now()) == 17 || force
        rm(beagle, force=true)
        rm(plink,  force=true)
    end

    subtitle("Updating the binaries")

    item("Beagle related")
    beagle2vcfURL = "https://faculty.washington.edu/browning/beagle_utilities/"
    beagle2vcf = "beagle2vcf.jar"

    isfile(beagle) || download(latest_beagle(), beagle)
    if !isfile(joinpath(bin_dir, beagle2vcf))
        message("\tDownloading beagle2vcf.jar")
        download(joinpath(beagle2vcfURL, beagle2vcf), joinpath(bin_dir, beagle2vcf))
        done("Newly downloaded")
    else
        done("Exists")
    end

    item("plink")
    if isfile(plink)
        done("Exists")
    else
        update_plink()
        done("Newly downloaded")
    end

    item("C++ binaries")
    make_bins()
    done()
end


text:
using PastaQ
using ITensors
using Test
using LinearAlgebra
using Random

@testset "quantum state fidelity: normalized input states" begin
  N = 4
  circuit1 = randomcircuit(N; depth = 3)
  circuit2 = randomcircuit(N; depth = 3)

  sites = siteinds("Qubit", N)
  # MPS wavefunction
  ψ1 = runcircuit(sites, circuit1)
  ψ2 = runcircuit(sites, circuit2)
  # MPO density matrix
  ρ1 = runcircuit(sites, circuit1; noise=("DEP", (p=0.01,)))
  ρ2 = runcircuit(sites, circuit2; noise=("DEP", (p=0.01,)))
  # LPDO density matrix
  ϱ1 = normalize!(randomstate(sites; χ = 2, ξ = 3))
  ϱ2 = normalize!(randomstate(sites; χ = 2, ξ = 3))

  ψ1vec = PastaQ.array(ψ1)
  ρ1mat = PastaQ.array(ρ1)
  ϱ1mat = PastaQ.array(ϱ1)

  ψ2vec = PastaQ.array(ψ2)
  ρ2mat = PastaQ.array(ρ2)
  ϱ2mat = PastaQ.array(ϱ2)

  @test fidelity(ψ1, ψ2) ≈ abs2(ψ1vec' * ψ2vec)
  @test fidelity(ψ1, ρ2) ≈ ψ1vec' * ρ2mat * ψ1vec
  @test fidelity(ψ1, ϱ2) ≈ real((ψ1vec' * ϱ2mat * ψ1vec))

  @test fidelity(ρ1, ρ2) ≈ real(tr(sqrt(sqrt(ρ1mat) * ρ2mat * sqrt(ρ1mat))))^2 atol = 1e-7
  @test fidelity(ρ1, ϱ2) ≈ real(tr(sqrt(sqrt(ρ1mat) * ϱ2mat * sqrt(ρ1mat))))^2 atol = 1e-7
  @test fidelity(ϱ1, ϱ2) ≈ real(tr(sqrt(sqrt(ϱ1mat) * ϱ2mat * sqrt(ϱ1mat))))^2 atol = 1e-7

  # ITensors

  ψ1prod = prod(ψ1)
  ψ2prod = prod(ψ2)
  ρ1prod = prod(ρ1)
  ρ2prod = prod(ρ2)
  ϱ1prod = prod(ϱ1)
  ϱ2prod = prod(ϱ2)

  @test fidelity(ψ1prod, ψ2prod) ≈ fidelity(ψ1, ψ2)
  @test fidelity(ψ1prod, ψ2) ≈ fidelity(ψ1, ψ2)
  @test fidelity(ψ1, ψ2prod) ≈ fidelity(ψ1, ψ2)

  @test fidelity(ψ1prod, ρ2prod) ≈ fidelity(ψ1, ρ2)
  @test fidelity(ψ1prod, ρ2) ≈ fidelity(ψ1, ρ2)
  @test fidelity(ψ1, ρ2prod) ≈ fidelity(ψ1, ρ2)

  @test fidelity(ψ1prod, ϱ2prod) ≈ fidelity(ψ1, ϱ2)
  @test fidelity(ψ1prod, ϱ2) ≈ fidelity(ψ1, ϱ2)
  @test fidelity(ψ1, ϱ2prod) ≈ fidelity(ψ1, ϱ2)

  @test fidelity(ρ1prod, ρ2prod) ≈ fidelity(ρ1, ρ2)
  @test fidelity(ρ1, ρ2prod) ≈ fidelity(ρ1, ρ2)
  @test fidelity(ρ1prod, ρ2) ≈ fidelity(ρ1, ρ2)

  @test fidelity(ρ1prod, ϱ2prod) ≈ fidelity(ρ1, ϱ2)
  @test fidelity(ρ1, ϱ2prod) ≈ fidelity(ρ1, ϱ2)
  @test fidelity(ρ1prod, ϱ2) ≈ fidelity(ρ1, ϱ2)

  @test fidelity(ϱ1prod, ϱ2prod) ≈ fidelity(ϱ1, ϱ2)
  @test fidelity(ϱ1prod, ϱ2) ≈ fidelity(ϱ1, ϱ2)
  @test fidelity(ϱ1, ϱ2prod) ≈ fidelity(ϱ1, ϱ2)

end

@testset "quantum state fidelity: unnormalized input states" begin
  N = 4
  sites = siteinds("Qubit",N)
  # MPS wavefunction
  ψ1 = randomstate(sites; χ = 4)
  ψ2 = randomstate(sites; χ = 5)
  # LPDO density matrix
  ϱ1 = randomstate(sites; χ = 5, ξ = 2)
  ϱ2 = randomstate(sites; χ = 5, ξ = 3)

  # MPO density matrix
  ρ1 = MPO(ϱ1)
  ρ2 = MPO(ϱ2)

  ψ1vec = PastaQ.array(normalize!(copy(ψ1)))
  ρ1mat = PastaQ.array(normalize!(copy(ρ1)))
  ϱ1mat = PastaQ.array(normalize!(copy(ϱ1)))

  ψ2vec = PastaQ.array(normalize!(copy(ψ2)))
  ρ2mat = PastaQ.array(normalize!(copy(ρ2)))
  ϱ2mat = PastaQ.array(normalize!(copy(ϱ2)))

  @test fidelity(ψ1, ψ2) ≈ abs2(ψ1vec' * ψ2vec)
  @test fidelity(ψ1, ρ2) ≈ ψ1vec' * ρ2mat * ψ1vec
  @test fidelity(ψ1, ϱ2) ≈ real((ψ1vec' * ϱ2mat * ψ1vec))

  @test fidelity(ρ1, ρ2) ≈ real(tr(sqrt(sqrt(ρ1mat) * ρ2mat * sqrt(ρ1mat))))^2 atol = 1e-7
  @test fidelity(ρ1, ϱ2) ≈ real(tr(sqrt(sqrt(ρ1mat) * ϱ2mat * sqrt(ρ1mat))))^2 atol = 1e-7
  @test fidelity(ϱ1, ϱ2) ≈ real(tr(sqrt(sqrt(ϱ1mat) * ϱ2mat * sqrt(ϱ1mat))))^2 atol = 1e-7

  # ITensors

  ψ1prod = prod(ψ1)
  ψ2prod = prod(ψ2)
  ρ1prod = prod(ρ1)
  ρ2prod = prod(ρ2)
  ϱ1prod = prod(ϱ1)
  ϱ2prod = prod(ϱ2)

  @test fidelity(ψ1prod, ψ2prod) ≈ fidelity(ψ1, ψ2)
  @test fidelity(ψ1prod, ψ2) ≈ fidelity(ψ1, ψ2)
  @test fidelity(ψ1, ψ2prod) ≈ fidelity(ψ1, ψ2)

  @test fidelity(ψ1prod, ρ2prod) ≈ fidelity(ψ1, ρ2)
  @test fidelity(ψ1prod, ρ2) ≈ fidelity(ψ1, ρ2)
  @test fidelity(ψ1, ρ2prod) ≈ fidelity(ψ1, ρ2)

  @test fidelity(ψ1prod, ϱ2prod) ≈ fidelity(ψ1, ϱ2)
  @test fidelity(ψ1prod, ϱ2) ≈ fidelity(ψ1, ϱ2)
  @test fidelity(ψ1, ϱ2prod) ≈ fidelity(ψ1, ϱ2)

  @test fidelity(ρ1prod, ρ2prod) ≈ fidelity(ρ1, ρ2)
  @test fidelity(ρ1, ρ2prod) ≈ fidelity(ρ1, ρ2)
  @test fidelity(ρ1prod, ρ2) ≈ fidelity(ρ1, ρ2)

  @test fidelity(ρ1prod, ϱ2prod) ≈ fidelity(ρ1, ϱ2)
  @test fidelity(ρ1, ϱ2prod) ≈ fidelity(ρ1, ϱ2)
  @test fidelity(ρ1prod, ϱ2) ≈ fidelity(ρ1, ϱ2)

  @test fidelity(ϱ1prod, ϱ2prod) ≈ fidelity(ϱ1, ϱ2)
  @test fidelity(ϱ1prod, ϱ2) ≈ fidelity(ϱ1, ϱ2)
  @test fidelity(ϱ1, ϱ2prod) ≈ fidelity(ϱ1, ϱ2)

end

@testset "quantum process fidelity" begin
  N = 3
  sites = siteinds("Qubit",N)

  circuit1 = randomcircuit(N; depth =  3)
  circuit2 = randomcircuit(N; depth =  3)
  # MPO unitary
  U1 = runcircuit(sites, circuit1; process=true)
  U2 = randomprocess(sites)

  # MPO Choi matrix
  ρ1 = PastaQ.choimatrix(sites, circuit1; noise=("DEP", (p=0.01,)))
  ρ2 = PastaQ.choimatrix(sites, circuit2; noise=("DEP", (p=0.01,)))
  # LPDO Choi matrix
  ϱ1 = normalize!(randomprocess(sites; mixed=true))
  ϱ2 = normalize!(randomprocess(sites; mixed=true))

  @disable_warn_order begin
    ϕ1 = PastaQ.unitary_mpo_to_choi_mps(U1)
    normalize!(ϕ1)
    ϕ1vec = PastaQ.array(ϕ1)
    ρ1mat = PastaQ.array(ρ1)
    ρ1mat = ρ1mat / tr(ρ1mat)
    ϱ1mat = PastaQ.array(ϱ1)
    ϱ1mat = ϱ1mat / tr(ϱ1mat)

    ϕ2 = PastaQ.unitary_mpo_to_choi_mps(U2)
    normalize!(ϕ2)
    ϕ2vec = PastaQ.array(ϕ2)
    ρ2mat = PastaQ.array(ρ2)
    ρ2mat = ρ2mat / tr(ρ2mat)
    ϱ2mat = PastaQ.array(ϱ2)
    ϱ2mat = ϱ2mat / tr(ϱ2mat)

    @test fidelity(U1, U2; process=true) ≈ abs2(ϕ1vec' * ϕ2vec)
    @test fidelity(U1, ρ2; process=true) ≈ ϕ1vec' * ρ2mat * ϕ1vec
    @test fidelity(U1, ϱ2; process=true) ≈ (ϕ1vec' * ϱ2mat * ϕ1vec)

    @test fidelity(ρ1, ρ2; process=true) ≈
          real(tr(sqrt(sqrt(ρ1mat) * ρ2mat * sqrt(ρ1mat))))^2 atol = 1e-7
    @test fidelity(ρ1, ϱ2; process=true) ≈
          real(tr(sqrt(sqrt(ρ1mat) * ϱ2mat * sqrt(ρ1mat))))^2 atol = 1e-7
    @test fidelity(ϱ1, ϱ2; process=true, cutoff = 1e-14) ≈
          real(tr(sqrt(sqrt(ϱ1mat) * ϱ2mat * sqrt(ϱ1mat))))^2 atol = 1e-7

    # ITensors
    U1prod = prod(U1)
    U2prod = prod(U2)
    ρ1prod = prod(ρ1)
    ρ2prod = prod(ρ2)
    ϱ1prod = prod(ϱ1)
    ϱ2prod = prod(ϱ2)

    @test fidelity(U1prod, U2prod; process=true) ≈ fidelity(U1, U2; process=true)
    @test fidelity(U1, U2prod; process=true) ≈ fidelity(U1, U2; process=true)
    @test fidelity(U1prod, U2; process=true) ≈ fidelity(U1, U2; process=true)

    @test fidelity(U1prod, ρ2prod; process=true) ≈ fidelity(U1, ρ2; process=true)
    @test fidelity(U1, ρ2prod; process=true) ≈ fidelity(U1, ρ2; process=true)
    @test fidelity(U1prod, ρ2; process=true) ≈ fidelity(U1, ρ2; process=true)

    @test fidelity(U1prod, ϱ2prod; process=true) ≈ fidelity(U1, ϱ2; process=true)
    @test fidelity(U1prod, ϱ2; process=true) ≈ fidelity(U1, ϱ2; process=true)
    @test fidelity(U1, ϱ2prod; process=true) ≈ fidelity(U1, ϱ2; process=true)
  end
end

@testset "quantum process fidelity: unnormalized states" begin
  N = 3
  sites = siteinds("Qubit",N)

  circuit1 = randomcircuit(N; depth =  3)
  circuit2 = randomcircuit(N; depth =  3)
  # MPO unitary
  U1 = randomprocess(sites; χ = 3)
  U2 = randomprocess(sites; χ = 4)

  # LPDO Choi matrix
  ϱ1 = randomprocess(sites; ξ = 3, χ = 2)
  ϱ2 = randomprocess(sites; ξ = 3, χ = 3)

  @disable_warn_order begin
    ϕ1 = PastaQ.unitary_mpo_to_choi_mps(U1)
    normalize!(ϕ1)
    ϕ1vec = PastaQ.array(ϕ1)
    ϱ1mat = PastaQ.array(ϱ1)
    ϱ1mat = ϱ1mat / tr(ϱ1mat)

    ϕ2 = PastaQ.unitary_mpo_to_choi_mps(U2)
    normalize!(ϕ2)
    ϕ2vec = PastaQ.array(ϕ2)
    ϱ2mat = PastaQ.array(ϱ2)
    ϱ2mat = ϱ2mat / tr(ϱ2mat)

    @test fidelity(U1, U2; process=true) ≈ abs2(ϕ1vec' * ϕ2vec)
    @test fidelity(U1, ϱ2; process=true) ≈ (ϕ1vec' * ϱ2mat * ϕ1vec)
    @test fidelity(ϱ1, ϱ2; process=true) ≈
          real(tr(sqrt(sqrt(ϱ1mat) * ϱ2mat * sqrt(ϱ1mat))))^2 atol = 1e-7

    ## ITensors
    U1prod = prod(U1)
    U2prod = prod(U2)
    ϱ1prod = prod(ϱ1)
    ϱ2prod = prod(ϱ2)

    @test fidelity(U1prod, U2prod; process=true) ≈ fidelity(U1, U2; process=true)
    @test fidelity(U1, U2prod; process=true) ≈ fidelity(U1, U2; process=true)
    @test fidelity(U1prod, U2; process=true) ≈ fidelity(U1, U2; process=true)

    @test fidelity(U1prod, ϱ2prod; process=true) ≈ fidelity(U1, ϱ2; process=true)
    @test fidelity(U1, ϱ2prod; process=true) ≈ fidelity(U1, ϱ2; process=true)
  end
end

@testset "frobenius distance" begin
  N = 4
  Random.seed!(1111)
  ψ1 = randomstate(N; χ=2)
  Random.seed!(2222)
  ψ2 = randomstate(ψ1; χ=2)

  ρ_mpo = MPO(ψ1)
  σ_mpo = MPO(ψ2)

  ρ_mat = PastaQ.array(ρ_mpo)
  σ_mat = PastaQ.array(σ_mpo)
  Kρ = tr(ρ_mat)
  Kσ = tr(σ_mat)

  δ = ρ_mat / Kρ - σ_mat / Kσ

  T = sqrt(tr(conj(transpose(δ)) * δ))

  F = frobenius_distance(ρ_mpo, σ_mpo)
  @test T ≈ F
  @test F ≈ frobenius_distance(ψ1, σ_mpo)
  @test F ≈ frobenius_distance(ρ_mpo, ψ2)
  @test F ≈ frobenius_distance(ψ1, ψ2)

  Random.seed!(1111)
  ρ = randomstate(ψ1; mixed=true, χ=2, ξ=2)

  ρ_mpo = MPO(ρ)
  σ_mpo = MPO(ψ2)

  ρ_mat = PastaQ.array(ρ_mpo)
  σ_mat = PastaQ.array(σ_mpo)
  Kρ = tr(ρ_mat)
  Kσ = tr(σ_mat)

  δ = ρ_mat / Kρ - σ_mat / Kσ

  T = sqrt(tr(conj(transpose(δ)) * δ))

  F = frobenius_distance(ρ, σ_mpo)
  @test T ≈ F
  @test F ≈ frobenius_distance(ρ_mpo, ψ2)

  Random.seed!(1111)
  σ = randomstate(ψ1; mixed=true, χ=2, ξ=2)

  ρ_mpo = MPO(ψ1)
  σ_mpo = MPO(σ)

  ρ_mat = PastaQ.array(ρ_mpo)
  σ_mat = PastaQ.array(σ_mpo)
  Kρ = tr(ρ_mat)
  Kσ = tr(σ_mat)

  δ = ρ_mat / Kρ - σ_mat / Kσ

  T = sqrt(tr(conj(transpose(δ)) * δ))

  F = frobenius_distance(ρ_mpo, σ)
  @test T ≈ F
  @test F ≈ frobenius_distance(ψ1, σ_mpo)

  Random.seed!(1111)
  ρ = randomstate(N; mixed=true, χ=2, ξ=2)
  Random.seed!(1111)
  σ = randomstate(ρ; mixed=true, χ=2, ξ=2)

  ρ_mpo = MPO(ρ)
  σ_mpo = MPO(σ)

  ρ_mat = PastaQ.array(ρ_mpo)
  σ_mat = PastaQ.array(σ_mpo)
  Kρ = tr(ρ_mat)
  Kσ = tr(σ_mat)

  δ = ρ_mat / Kρ - σ_mat / Kσ

  T = sqrt(tr(conj(transpose(δ)) * δ))

  F = frobenius_distance(ρ, σ)
  @test T ≈ F
end

@testset "fidelity bound" begin
  N = 4
  Random.seed!(1111)
  ψ1 = randomstate(N; χ=2)
  Random.seed!(2222)
  ψ2 = randomstate(ψ1; χ=2)

  ρ_mpo = MPO(ψ1)
  σ_mpo = MPO(ψ2)

  ρ_mat = PastaQ.array(ρ_mpo)
  σ_mat = PastaQ.array(σ_mpo)
  Kρ = tr(ρ_mat)
  Kσ = tr(σ_mat)

  f = tr(conj(transpose(ρ_mat / Kρ)) * (σ_mat / Kσ))
  F̃ = fidelity_bound(ρ_mpo, σ_mpo)
  @test f ≈ F̃
  @test F̃ ≈ fidelity(ψ1, σ_mpo)
  @test F̃ ≈ fidelity(ρ_mpo, ψ2)
  @test F̃ ≈ fidelity(ψ1, ψ2)

  Random.seed!(1111)
  ρ = randomstate(ψ1; mixed=true, χ=2, ξ=2)

  ρ_mpo = MPO(ρ)
  σ_mpo = MPO(ψ2)

  ρ_mat = PastaQ.array(ρ_mpo)
  σ_mat = PastaQ.array(σ_mpo)
  Kρ = tr(ρ_mat)
  Kσ = tr(σ_mat)

  f = tr(conj(transpose(ρ_mat / Kρ)) * (σ_mat / Kσ))
  F̃ = fidelity_bound(ρ, σ_mpo)
  @test f ≈ F̃
  @test F̃ ≈ fidelity(ρ_mpo, ψ2)

  Random.seed!(1111)
  σ = randomstate(ψ1; mixed=true, χ=2, ξ=2)

  ρ_mpo = MPO(ψ1)
  σ_mpo = MPO(σ)

  ρ_mat = PastaQ.array(ρ_mpo)
  σ_mat = PastaQ.array(σ_mpo)
  Kρ = tr(ρ_mat)
  Kσ = tr(σ_mat)
  f = tr(conj(transpose(ρ_mat / Kρ)) * (σ_mat / Kσ))
  F̃ = fidelity_bound(ρ_mpo, σ)
  @test f ≈ F̃
  @test F̃ ≈ fidelity(ψ1, σ_mpo)

  Random.seed!(1111)
  ρ = randomstate(N; mixed=true, χ=2, ξ=2)
  Random.seed!(1111)
  σ = randomstate(ρ; mixed=true, χ=2, ξ=2)

  ρ_mpo = MPO(ρ)
  σ_mpo = MPO(σ)

  ρ_mat = PastaQ.array(ρ_mpo)
  σ_mat = PastaQ.array(σ_mpo)
  Kρ = tr(ρ_mat)
  Kσ = tr(σ_mat)

  f = tr(conj(transpose(ρ_mat / Kρ)) * (σ_mat / Kσ))
  F̃ = fidelity_bound(ρ, σ)
  @test f ≈ F̃
end


text:
# Create example data
x = randn(100_000)
y = randn(100_000)
# Draw the hexbin plot
hexbin(x, y)


text:
"A container for the bedGraph header."
mutable struct BedgraphHeader{T} #TODO: determine what and how this will be.
    data::T
end
BedgraphHeader{T}() where T = BedgraphHeader{T}(T())
BedgraphHeader() = BedgraphHeader{Vector{String}}()

"Push data into [`BedgraphHeader`](@ref) container."
function Base.push!(sink::BedgraphHeader, data)
    push!(sink.data, data) #Note: converts data to sink.data's eltype.
end

"Convert [`BedgraphHeader`](@ref) to type `String`."
function Base.convert(::Type{String}, header::BedgraphHeader{<:AbstractVector{<:AbstractString}})

    str = ""
    for line in header.data
        str = string(str, line, '\n')
    end

    return str
end

"Generate a basic bedGraph header given a vector of [`Record`](@ref)s."
function generate_basic_header(records::AbstractVector{<:Record}; bump_forward=true) #Note: we assume that records are sorted by chrom and left position.

    chrom = records[1].chrom

    pos_start = records[1].first
    pos_end = records[end].last

    if bump_forward
        pos_start = pos_start + 1
        pos_end = pos_end + 1
    end

    return BedgraphHeader(["browser position $chrom:$pos_start-$pos_end", "track type=bedGraph"])
end

generate_basic_header(chrom::AbstractString, pos_start::Int, pos_end::Int; bump_forward=true) = generate_basic_header([Record(chrom, pos_start, pos_end, 0)], bump_forward=bump_forward)

"Seek and then read bedGraph header into sink."
function read_header(io::IO, sink)
    position(io) == 0 || seekstart(io)

    line = readline(io)

    while !eof(io) && !is_like_record(line) # TODO: seek more rebust check.
        push!(sink, line)
        line = readline(io)
    end

    return sink
end

"Read bedGraph header into sink."
function Base.read(io::IO, sink::Type{<:BedgraphHeader})
    return read_header(io, sink())
end

"Write bedGraph header to `IO` as `String`."
function Base.write(io::IO, header::BedgraphHeader)
    return Base.write(io, convert(String, header))
end


text:
#
# @description Implements the MD5 hash algorithm as defined in RFC 1321.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2015 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
export Md5Hasher
export MD5_DIGEST_SIZE, MD5_BLOCK_SIZE
export digest_size, reset!, update!, digest

# @description Defines the MD5 algorithm data struture.
mutable struct Md5Hasher <: HashFunction
  scratch::Vector{UInt32}
  partial_block::Vector{UInt8}
  block_size::UInt32
  data_length::UInt64

  Md5Hasher() = reset!(new())
end

# @description The MD5 resulting message digest size in bytes.
const MD5_DIGEST_SIZE = 16

# @description The MD5 operation block size in bytes.
const MD5_BLOCK_SIZE = 64

# @description Gets the size of the resulting message digest in bytes.
#
# @param {Md5Hasher} self The hash algorithm.
#
# @return {Integer} The size of the resulting message digest in bytes.
function digest_size(self::Md5Hasher)
  return MD5_DIGEST_SIZE
end

# @description Resets the state of the MD5 data structure.
#
# @param {Md5Hasher} self The hash algorithm.
#
# @return The reference to the resetted hash algorithm.
function reset!(self::Md5Hasher)
  self.scratch       = UInt32[ 0x67452301, 0xEFCDAB89, 0x98BADCFE, 0x10325476 ]
  self.partial_block = zeros(UInt8, MD5_BLOCK_SIZE)
  self.block_size    = 0
  self.data_length   = 0

  return self
end

# @description Continue hashing of a message by consuming the next chunk of data.
#
# @param  {Md5Hasher}     self The hash algorithm.
# @param  {Vector{UInt8}} data The next message chunk.
#
# @return The reference to the updated hash algorithm.
function update!(self::Md5Hasher, data::Vector{UInt8})
  index            = 0
  data_length      = length(data)
  self.data_length = self.data_length + data_length

  if self.block_size > 0
    data = vcat(self.partial_block[1 : self.block_size], data)
    self.block_size = 0
    data_length     = length(data)
  end

  while data_length >= (index + MD5_BLOCK_SIZE)
    process_block(self, data[(index + 1):(index + MD5_BLOCK_SIZE)])
    index += MD5_BLOCK_SIZE
  end

  if data_length > index
    setindex!(self.partial_block, data[(index + 1) : end], 1 : (data_length - index))
    self.block_size = data_length - index
  end

  return self
end

# @description Gets the digest of the data passed to the update!() method so far.
#
# @param {Md5Hasher} self The hash algorithm.
#
# @return {Vector{UInt8}} The message digest.
function digest(self::Md5Hasher)
  algorithm   = deepcopy(self)
  r           = algorithm.data_length % MD5_BLOCK_SIZE
  extra       = (r > 56) ? MD5_BLOCK_SIZE : 0
  last_block  = zeros(UInt8, extra + MD5_BLOCK_SIZE - r)
  len::UInt64 = algorithm.data_length << 3

  last_block[1]       = 0x80
  last_block[end - 7] = convert(UInt8, len & 0xFF)

  for i in 1:7
    last_block[end - 7 + i] = convert(UInt8, (len >> (8 * i)) & 0xFF)
  end

  update!(algorithm, last_block)
  @assert (algorithm.block_size == 0) "Block should be empty."

  result = zeros(UInt8, MD5_DIGEST_SIZE)

  for i in 1:length(algorithm.scratch)
    value = bitstring(algorithm.scratch[i])
    result[4 * (i - 1) + 1] = parse(UInt8, value[(end -  7) : (end -  0)], base = 2)
    result[4 * (i - 1) + 2] = parse(UInt8, value[(end - 15) : (end -  8)], base = 2)
    result[4 * (i - 1) + 3] = parse(UInt8, value[(end - 23) : (end - 16)], base = 2)
    result[4 * (i - 1) + 4] = parse(UInt8, value[(end - 31) : (end - 24)], base = 2)
  end

  return result
end

# PRIVATE IMPLEMENTATION #######################################################

# @description MD5 block processing shift constants.
const MD5_SHIFTS = UInt32[
  7, 12, 17, 22, 7, 12, 17, 22,  7, 12, 17, 22, 7, 12, 17, 22,
  5,  9, 14, 20, 5,  9, 14, 20,  5,  9, 14, 20, 5,  9, 14, 20,
  4, 11, 16, 23, 4, 11, 16, 23,  4, 11, 16, 23, 4, 11, 16, 23,
  6, 10, 15, 21, 6, 10, 15, 21,  6, 10, 15, 21, 6, 10, 15, 21
]

# @description MD5 block processing bytes constants.
const MD5_CONSTANTS_TABLE = UInt32[
  0xd76aa478, 0xe8c7b756, 0x242070db, 0xc1bdceee, 0xf57c0faf, 0x4787c62a, 0xa8304613, 0xfd469501,
  0x698098d8, 0x8b44f7af, 0xffff5bb1, 0x895cd7be, 0x6b901122, 0xfd987193, 0xa679438e, 0x49b40821,
  0xf61e2562, 0xc040b340, 0x265e5a51, 0xe9b6c7aa, 0xd62f105d, 0x02441453, 0xd8a1e681, 0xe7d3fbc8,
  0x21e1cde6, 0xc33707d6, 0xf4d50d87, 0x455a14ed, 0xa9e3e905, 0xfcefa3f8, 0x676f02d9, 0x8d2a4c8a,
  0xfffa3942, 0x8771f681, 0x6d9d6122, 0xfde5380c, 0xa4beea44, 0x4bdecfa9, 0xf6bb4b60, 0xbebfbc70,
  0x289b7ec6, 0xeaa127fa, 0xd4ef3085, 0x04881d05, 0xd9d4d039, 0xe6db99e5, 0x1fa27cf8, 0xc4ac5665,
  0xf4292244, 0x432aff97, 0xab9423a7, 0xfc93a039, 0x655b59c3, 0x8f0ccc92, 0xffeff47d, 0x85845dd1,
  0x6fa87e4f, 0xfe2ce6e0, 0xa3014314, 0x4e0811a1, 0xf7537e82, 0xbd3af235, 0x2ad7d2bb, 0xeb86d391
]

# @description Process a block.
#
# @param {Md5Hasher}     self  The hash algorithm.
# @param {Vector{UInt8}} block The block data buffer.
function process_block(self::Md5Hasher, block::Vector{UInt8})
  a0::UInt32 = a::UInt32 = self.scratch[1]
  b0::UInt32 = b::UInt32 = self.scratch[2]
  c0::UInt32 = c::UInt32 = self.scratch[3]
  d0::UInt32 = d::UInt32 = self.scratch[4]

  buffer = [ bytes_to_int(block[i], block[i + 1], block[i + 2], block[i + 3]) for i in 1 : 4 : 64 ]

  for j in 0:63
    f::UInt32 = 0
    position  = j
    round     = j >> 4

    if round == 0
      f = (b & c) | (~b & d)
    elseif round == 1
      f = (b & d) | (c & ~d)
      position = (position * 5 + 1) & 0x0F
    elseif round == 2
      f = b ⊻ c ⊻ d
      position = (position * 3 + 5) & 0x0F
    else
      f = c ⊻ (b | ~d)
      position = (position * 7) & 0x0F
    end

    sa  = MD5_SHIFTS[j + 1]
    a  += f + buffer[position + 1] + MD5_CONSTANTS_TABLE[j + 1]

    a, d, c, b = (d, c, b, (a << sa | a >> (32 - sa) + b))
  end

  self.scratch[1] = a0 + a
  self.scratch[2] = b0 + b
  self.scratch[3] = c0 + c
  self.scratch[4] = d0 + d
end

# @description Packs the bytes into a integer value.
#
# @param {UInt8} a, b, c, d The bytes.
#
# @return {UInt32} The integer.
function bytes_to_int(a, b, c, d)
  return convert(UInt32, a) | convert(UInt32, b) << 8 | convert(UInt32, c) << 16 | convert(UInt32, d) << 24
end


text:
# deprecations
@deprecate iparameters(args...) getiparams(args...)
@deprecate setiparameters!(args...) setiparams!(args...)
@deprecate niparameters(args...) niparams(args...)
@deprecate parameter_type(args...) parameters_eltype(args...)
@deprecate PreserveStyle(args...) PropertyTrait(args...)
@deprecate Sum(args...) Add(args...)
@deprecate Add(blocks::AbstractVector{<:AbstractBlock{N}}) where {N} Add{N}(blocks)
@deprecate mathgate(f; nbits) mathgate(nbits, f)


text:
mutable struct MCMC{T <: EpidemicModel}
  event_observations::EventObservations{T}
  event_extents::EventExtents{T}
  population::Population
  starting_states::Vector{DiseaseState}
  risk_functions::RiskFunctions{T}
  risk_priors::RiskPriors{T}
  transmission_network_prior::Union{Nothing, TNPrior}
  markov_chains::Vector{MarkovChain{T}}

  function MCMC(obs::EventObservations{T},
                ee::EventExtents{T},
                pop::Population,
                states::Vector{DiseaseState},
                rf::RiskFunctions{T},
                rp::RiskPriors{T};
                tnprior::Union{Nothing, TNPrior}=nothing) where T <: EpidemicModel
    return new{T}(obs, ee, pop, states, rf, rp, tnprior, MarkovChain{T}[])
  end
end

function MCMC(obs::EventObservations{T},
              ee::EventExtents{T},
              pop::Population,
              rf::RiskFunctions{T},
              rp::RiskPriors{T};
              tnprior::Union{Nothing, TNPrior}=nothing) where T <: EpidemicModel
  return MCMC(obs, ee, pop, fill(State_S, pop.individuals), rf, rp, tnprior=tnprior)
end

function Base.show(io::IO, x::MCMC{T}) where T <: EpidemicModel
  return print(io, "$T model MCMC with $(length(x.markov_chains)) chains")
end

function TransmissionNetworkDistribution(x::MCMC; burnin::Int64=0, thin::Int64=1)
  return TNDistribution([TNDistribution(y, burnin=burnin, thin=thin) for y in x.markov_chains])
end


text:
using Test
using Random
using Distributions: Bernoulli, LogNormal, Poisson
using Statistics:std

using BayesianExperiments
include("utils.jl")

include("util.jl")
include("data.jl")
include("experiment.jl")
include("conjugate.jl")
include("bayesfactor.jl")
include("simulation.jl")


text:
export generate_psf
export generate_downsample, my_interpolate
export center_extract, center_set!, get_indices_around_center, center_pos
export conv_psf, conv_otf, conv_otf_r, plan_conv_r

"""
    conv_psf(obj, psf [, dims])

Convolve `obj` with `psf` over `dims` dimensions.
Based on FFT convolution.
"""
function conv_psf(obj, psf, dims=[1, 2])
    return real.(ifft(fft(obj, dims) .* fft(psf, dims), dims))
end

"""
    conv_otf(obj, otf [, dims])

Performs a FFT-based convolution of an `obj`
with an `otf`. `otf = fft(psf)`. The 0 frequency of the `otf` must be located
at position [1, 1, 1].
The `obj` can be of arbitrary dimension but `ndims(psf) ≥ ndims(otf)`.
The convolution happens over the `dims` array. Any further dimensions are broadcasted.
Per default `dims = [1, 2]`.
"""
function conv_otf(obj, otf, dims=[1, 2])
    return real.(ifft(fft(obj, dims) .* otf, dims))
end

"""
    conv_otf_r(obj, otf [, dims])

Performs a FFT-based convolution of an `obj`
with an `otf`.
Same arguments as `conv_otf` but with `obj` being real and `otf=rfft(psf)`.
"""
function conv_otf_r(obj, otf, dims=[1, 2])
    return real.(irfft(rfft(obj, dims) .* otf, size(obj)[dims[1]], dims))
end

"""
    plan_conv_r(psf [, dims])

Pre-plan an optimized convolution for array shaped like `psf` (based on pre-plan FFT)
along the given dimenions `dims`.
`dims = [1, 2]` per default.
The 0 frequency of `psf` must be located at [1, 1, 1].
We return first the `otf` (obtained by `rfft(psf))`.
The second return is the convolution function `conv`.
`conv` itself has two arguments. `conv(obj, otf)` where `obj` is the object and `otf` the otf.

This function achieves faster convolution than `conv_psf(obj, psf)`.
"""
function plan_conv_r(psf, rec0, dims=[1, 2])
    # do the preplanning step
    P = plan_rfft(rec0, dims)
    rec0_fft = P * rec0
    P_inv = plan_irfft(rec0_fft, size(psf)[dims[1]], dims)

    # obtain the otf by real based fft
    otf = rfft(psf, dims)
    # construct the efficient conv function
    conv(obj, otf) = real.(P_inv * ((P * obj) .* otf))

    return otf, conv
end

"""
    generate_downsample(num_dim, downsample_dims, factor)

Generate a function (based on Tullio.jl) which can be used to downsample arrays.
`num_dim` (Integer) are the dimensions of the array.
`downsample_dims` is a list of which dimensions should be downsampled.
`factor` is a downsampling factor. It needs to be an integer number.

# Examples
```jldoctest
julia> ds = generate_downsample(2, [1, 2], 2)
[...]
julia> ds([1 2; 3 4; 5 6; 7 8])
2×1 Array{Float64,2}:
 2.5
 6.5

julia> ds = generate_downsample(2, [1], 2)
[...]
julia> ds([1 2; 3 5; 5 6; 7 8])
2×2 Array{Float64,2}:
 2.0  3.5
 6.0  7.0
```
"""
function generate_downsample(num_dim, downsample_dims, factor)
    @assert num_dim ≥ length(downsample_dims)
    # create unit cell with Cartesian Index
    # dims_units containts every where a 1 where the downsampling should happen
    dims_units = zeros(Int, num_dim)
    # here we set which dimensions should be downsamples
    dims_units[downsample_dims] .= 1
    # the unit cell expressed in CartesianIndex
    one = CartesianIndex(dims_units...)

    # create a list of symbols
    # these list represents the symbols to access the arrays
    ind = :i
    inds_out = map(1:num_dim) do di
        i = Symbol(ind, di)
    end

    # output list for the add commands
    add = []
    # via CartesianIndex we can loop over all rectangular neighbours
    # we loop only over the neighbours in the downsample_dims
    for n = one:one * factor
        # for each index calculate the offset to the neighbour
        inds = map(1:num_dim) do di
                i = Symbol(ind, di)
                if n[di] == 0
                    di = i
                else
                    expr = :($factor * $i)
                    diff = -factor + n[di]
                    di = :($expr + $diff)
                end
            end
        # push this single neighbour to add list
        push!(add, :(arr[$(inds...)]))
    end
    # combine the different parts and divide for averaging
    expr = [:(@tullio res[$(inds_out...)] := (+($(add...))) / $factor ^ $(length(downsample_dims)))]
    #= return expr =#
    # evaluate to function
    @eval f = arr -> ($(expr...))
    return f
end

"""
    my_interpolate(arr, size_n, [interp_type])

Interpolates `arr` to the sizes provided in `size_n`.
Therefore it holds `ndims(arr) == length(size_n)`.
`interp_type` specifies the interpolation type.
See Interpolations.jl for all options
"""
function my_interpolate(arr, size_n, interp_type=BSpline(Linear()))

    # we construct a arr which includes the interpolation
    # type for each dimension
    interp = []
    for s in size_n
        # if the outpute size is of the s-th dimension=1,
        # do NoInterp
        if s == 1
            push!(interp, NoInterp())
        else
            push!(interp, interp_type)
        end
    end
    # prepare the interpolation
    arr_n = interpolate(arr, Tuple(interp))

    # interpolate introduces fractional indices
    # via LinRange we access these fractional indices
    inds = []
    for d = 1:ndims(arr)
        push!(inds, LinRange(1, size(arr)[d], size_n[d]))
    end

    # return the new array sampled at the positions of inds
    # this accessing actually interpolates the data
    return arr_n(inds...)
end

"""
    get_indices_around_center(i_in, i_out)
A function which provides two output indices `i1` and `i2`
where `i2 - i1 = i_out`
The indices are chosen in a way that the set `i1:i2`
cuts the interval `1:i_in` in a way that the center frequency
stays at the center position.
Works for both odd and even indices
"""
function get_indices_around_center(i_in, i_out)
    if (mod(i_in, 2) == 0 && mod(i_out, 2) == 0
     || mod(i_in, 2) == 1 && mod(i_out, 2) == 1)
        x = (i_in - i_out) ÷ 2
        return 1 + x, i_in - x
    elseif mod(i_in, 2) == 1 && mod(i_out, 2) == 0
        x = (i_in - 1 - i_out) ÷ 2
        return 1 + x, i_in - x - 1
    elseif mod(i_in, 2) == 0 && mod(i_out, 2) == 1
        x = (i_in - (i_out - 1)) ÷ 2
        return 1 + x, i_in - (x - 1)
    end
end

"""
    center_extract(arr, new_size_array)
Extracts a center of an array.
`new_size_array` must be list of sizes indicating the output
size of each dimension. Centered means that a center frequency
stays at the center position. Works for even and uneven.
If `length(new_size_array) < length(ndims(arr))` the remaining dimensions
are untouched and copied.
# Examples
```jldoctest
julia> DeconvOptim.center_extract([1 2; 3 4], [1])
1×2 Array{Int64,2}:
 3  4
julia> DeconvOptim.center_extract([1 2; 3 4], [1, 1])
1×1 Array{Int64,2}:
 4
julia> DeconvOptim.center_extract([1 2 3; 3 4 5; 6 7 8], [2 2])
2×2 Array{Int64,2}:
 1  2
 3  4
```
"""
function center_extract(arr::AbstractArray, new_size_array)
    new_size_array = collect(new_size_array)

    # we construct two lists
    # the reason is, that we don't change higher dimensions which are not
    # specified in new_size_array
    out_indices1 = [get_indices_around_center(size(arr)[x], new_size_array[x])
                    for x = 1:length(new_size_array)]

    out_indices1 = [x[1]:x[2] for x = out_indices1]

    # out_indices2 contains just ranges covering the full size of each dimension
    out_indices2 = [1:size(arr)[i] for i = (1 + length(new_size_array)):ndims(arr)]
    return arr[out_indices1..., out_indices2...]
end

"""
    center_set!(arr_large, arr_small)
Puts the `arr_small` central into `arr_large`.
The convention, where the center is, is the same as the definition
as for FFT based centered.
Function works both for even and uneven arrays.
# Examples
```jldoctest
julia> DeconvOptim.center_set!([1, 1, 1, 1, 1, 1], [5, 5, 5])
6-element Array{Int64,1}:
 1
 1
 5
 5
 5
 1
```
"""
function center_set!(arr_large, arr_small)
    out_is = []
    for i = 1:ndims(arr_large)
        a, b = get_indices_around_center(size(arr_large)[i], size(arr_small)[i])
        push!(out_is, a:b)
    end

    #rest = ones(Int, ndims(arr_large) - 3)
    arr_large[out_is...] = arr_small

    return arr_large
end

"""
    center_pos(x)
Calculate the position of the center frequency.
Size of the array is `x`
# Examples
```jldoctest
julia> DeconvOptim.center_pos(3)
2
julia> DeconvOptim.center_pos(4)
3
```
"""
function center_pos(x::Integer)
    # integer division
    return div(x, 2) + 1
end

"""
    generate_psf(psf_size, radius)

Generation of an approximate 2D PSF.
`psf_size` is the output size of the PSF. The PSF will be centered
around the point [1, 1],
`radius` indicates the pupil diameter in pixel from which the PSF is generated.

# Examples
```julia-repl
julia> generate_psf([5, 5], 2)
5×5 Array{Float64,2}:
 0.36       0.104721    0.0152786    0.0152786    0.104721
 0.104721   0.0304627   0.00444444   0.00444444   0.0304627
 0.0152786  0.00444444  0.000648436  0.000648436  0.00444444
 0.0152786  0.00444444  0.000648436  0.000648436  0.00444444
 0.104721   0.0304627   0.00444444   0.00444444   0.0304627
```
"""
function generate_psf(psf_size, radius)
    mask = rr_2D(psf_size) .<= radius
    mask_ft = fft(mask)
    psf = abs2.(mask_ft)
    return psf ./ sum(psf)
end

function rr_3D(s)
    rarr = zeros((s...))
    for k = 1:s[3]
        for j = 1:s[2]
            for i = 1:s[1]
                rarr[i, j, k] = sqrt( (i-center_pos(s[1]))^2 + (j-center_pos(s[2]))^2 + (k-center_pos(s[3]))^2)
            end
        end
    end
    return rarr
end

"""
    rr_2D(s)

Generate a image with values being the distance to the center pixel.
`s` specifies the output size of the 2D array.

# Examples
```julia-repl
julia> DeconvOptim.rr_2D((6, 6))
6×6 Array{Float64,2}:
 4.24264  3.60555  3.16228  3.0  3.16228  3.60555
 3.60555  2.82843  2.23607  2.0  2.23607  2.82843
 3.16228  2.23607  1.41421  1.0  1.41421  2.23607
 3.0      2.0      1.0      0.0  1.0      2.0
 3.16228  2.23607  1.41421  1.0  1.41421  2.23607
 3.60555  2.82843  2.23607  2.0  2.23607  2.82843
```
"""
function rr_2D(s)
    rarr = zeros((s...))
    for j = 1:s[2]
        for i = 1:s[1]
               rarr[i, j] = sqrt( (i - center_pos(s[1]))^2 + (j - center_pos(s[2]))^2)
        end
    end
    return rarr
end


text:
# .................... cubic b-spline functions ...............................

# n / 0 == 0
#
function zero_divide(a, b)
    x = a / b
    return ifelse(isfinite(x), x, 0)
end

# uniform-knot p:degree n:number_of_points
#
function knot(p, n)
    base = [ifelse(k <= p, 0.0, 1.0) for k in (0:p + n)]
    if n - p == 1
        return base
    else
        for k in 0:n-p
            base[p + k + 1] = k / (n - p)
        end
    end
    return base
end

# b-spline basis p:degree i:number u:parameter k:knot
#
function bbasis(p, i, u, k)
    if p == 0
        return ifelse(k[i] <= u <= k[i+1], 1, 0)
    else
        m0 = zero_divide(u - k[i], k[i+p] - k[i])
        m1 = zero_divide(k[i+p+1] - u, k[i+p+1] - k[i+1])
        return m0 * bbasis(p-1, i, u, k) + m1 * bbasis(p-1, i+1, u, k)
    end
end

# b-spline basis first-derivative
#
function bbasis_fd(p, i, u, k)
    return zero_divide(p, k[i+p] - k[i]) * bbasis(p-1, i, u, k) -
        zero_divide(p, k[i+p+1] - k[i+1]) * bbasis(p-1, i+1, u, k)
end

# b-spline curve, points: as columns, t: 1d vector
#
function bcurve(points, u, bas=bbasis)
    knt = Ref(knot(3, size(points, 2)))
    tmp = zeros(size(points, 1), length(u), size(points, 2))
    for n in axes(points, 2)
        tmp[:,:,n] = points[:,n] .* bas.(3, n, u, knt)'
    end
    return dropdims(sum(tmp, dims=3), dims=3)
end

# b-spline surface, points: size(cp_u, cp_v, 3), uv: size(nu, nv, 2)
#
function bsurface(points, uv, bu=bbasis, bv=bbasis)
    ku, kv = knot.(3, size(points))
    uv_tmp = zeros(size(uv, 1), size(uv, 2), size(points, 3))

    for I in CartesianIndices(uv[:,:,1])
        for J in CartesianIndices(points[:,:,1])
            uv_tmp[I,:] += points[J,:] .*
                bu(3, J[1], uv[I, 1], ku) .*
                bv(3, J[2], uv[I, 2], kv)
        end
    end
    return uv_tmp
end

function mgrid(nu, nv, a=0:1, b=0:1)
    uv = zeros(nu, nv, 2)
    uv[:,:,1] .= first(a) : ((last(a)-first(a)) / (nu-1)) : last(a)
    uv[:,:,2] .= (first(b) : ((last(b)-first(b)) / (nv-1)) : last(b))'
    return uv
end

function rnd_grid(x, y)
    uv = mgrid(x, y, 0:x, 0:y)
    az = rand(Float32, x, y)
    return cat(uv, az, dims = 3)
end

max_stars_repo_name: zaika-generative/ShapeFactory

max_issues_repo_name: zaika-generative/ShapeFactory

max_forks_repo_name: zaika-generative/ShapeFactory


text:
# Analyze convergence of Monte Carlo tree search for various horizon lenghts

using SubmodularMaximization
using PyPlot
using Statistics
using Base.Threads
using Base.Iterators
using JLD2
using Printf
using HistogramFilters

close("all")

experiment_name = "evaluate_solver_objectives"
data_folder = "./data"

reprocess = false

# Values taken from the entropy_by_solver script
num_robots = 16
trials = 1:20
solver_ind = 4 # dist-4

# Drop the first n and keep every few after
sample_inds = 20:100

library_solver_ind = 4 # dist-4

# Range limits
communication_range = 20
robot_target_range_limit = 12
range_limit_partitions = 4
range_limit_solver(x) = solve_communication_range_limit(x,
                                                        num_partitions=
                                                        range_limit_partitions,
                                                        communication_range=
                                                        communication_range)

solver_rounds = [2, 4, 8]
solvers = [(solve_random, Inf),
           (solve_myopic, Inf),
           map(num->(prob->solve_n_partitions(num, prob), Inf), solver_rounds)...,
           (range_limit_solver, robot_target_range_limit),
           (solve_sequential, Inf)]
solver_strings = ["random",
                  "myopic",
                  map(x->string("dist. ",x," rnds."), solver_rounds)...,
                  string("rng. lim. ", range_limit_partitions, " rnds."),
                  "sequential"]
solver_inds = 1:length(solvers)

#
# Load prior trials
#

println("Loading rollouts")
library_file = "./data/entropy_by_solver.jld2"

let
  @load library_file results
  global library_results = results
end

# Create a map from number of robots to a vector of configs sample pairs
samples = vcat(map(trials) do trial
                 rollout = library_results[library_solver_ind, num_robots, trial]

                 [(configs=rollout.configs, sample=sample)
                  for sample in rollout.data[sample_inds]]
               end...)

#
# Evaluate results on solvers
#

all_tests = product(solver_inds, 1:length(samples))

function trial_fun(x)
  solver_ind, sample_ind = x
  solver, robot_target_range_limit = solvers[solver_ind]

  test = samples[sample_ind]

  sample = test.sample
  configs = MultiRobotTargetTrackingConfigs(test.configs,
                                            robot_target_range_limit=
                                            robot_target_range_limit
                                           )

  problem = MultiRobotTargetTrackingProblem(sample.robot_states,
                                            sample.target_filters,
                                            configs)
  solution = solver(problem)

  solution.value
end
function print_summary(x)
  solver_ind, sample_ind = x
  _, robot_target_range_limit = solvers[solver_ind]

  println("Solver: ", solver_strings[solver_ind],
          " Robot-Target Range Limit: ", robot_target_range_limit,
          " Sample: ", sample_ind)
end

@time results = run_experiments(all_tests,
                                trial_fun=trial_fun,
                                print_summary=print_summary,
                                experiment_name=experiment_name,
                                data_folder=data_folder,
                                reprocess=reprocess
                               )

println("Analyzing results")

#
# Results by objective values
#

results_matrix = map(x->results[x], all_tests)
boxplot(results_matrix', notch=false, vert=false)

title("Solver Objective Performance")
xlabel("Mutual Information Objective (bits)")
yticks(solver_inds, solver_strings)

save_fig("fig", "solver_objectives")

#
# Results normalized by maximum objective values on each problem
#

# Max values over different solvers for each problem instance
figure()

maximum_values = maximum(results_matrix, dims=1)
normalized_results = results_matrix ./ maximum_values

boxplot(normalized_results', notch=false, vert=false)

title("Normalized Objective Performance")
xlabel("Objective (fraction of max)")
yticks(solver_inds, solver_strings)

save_fig("fig", "normalized_solver_objectives")

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
# start the netcat (nc) tool at the prompt to make a connection with the Julia server on port 8080:
# nc localhost 8080
server = listen(8080)
while true
  conn = accept(server)
  @async begin
    try
      while true
        line = readline(conn)
        println(line)
        write(conn,line)
      end
    catch ex
      print("connection ended with error $ex")
    end
  end # end coroutine block
end

max_stars_repo_path: Module 1/Chapter08/echoserver.jl

max_issues_repo_path: Module 1/Chapter08/echoserver.jl

max_forks_repo_path: Module 1/Chapter08/echoserver.jl


text:
"""
    Spherical <: Kernel

TODO
"""
struct Spherical <: Kernel
    # TODO
end


text:
# This file was generated, do not modify it. # hide
using MLJ
import DataFrames
import Statistics
using PrettyPrinting
using StableRNGs

MLJ.color_off() # hide

rng = StableRNG(512)
Xraw = rand(rng, 300, 3)
y = exp.(Xraw[:,1] - Xraw[:,2] - 2Xraw[:,3] + 0.1*rand(rng, 300))
X = DataFrames.DataFrame(Xraw)

train, test = partition(eachindex(y), 0.7);


text:
#=
constants meaning

τ = the time grid
Et = electric field
Ip = ionization energy for the gas
ns, cn2, wt = constants from adk paper
W = ionization probability
Na = nutral atoms density
Ne = electrons density

p = pressure
Tk = temperture in kelven
c = speed of light
KB = boltezman constant
J = quantum number
h = plancks' constant
e = electron charge
hbar = reduced planck constant
N0 = particles density at atomsperic pressure
m_e = electron mass
ϵ0 = vaccume premitivity
nJ = number of quantum states to be used in Raman QSM
pre = Raman normalization factor obtaied from getpre() function
Bk =The rotational constants
τrot_k = The relaxation factors
=#

const e = 1.60217662e-19
const hbar = 1.0545718e-34
const N0 =  2.5000013629442974e25 # number of gas atoms(molecules) in 1 m^3 air at standard condition
const m_e = 9.10938356e-31
const ϵ0 = 8.854e-12
const Tk = 293.0
const c = 3e8
const KB = 1.38064852e-23
const h = 6.62607004e-34

const p_0 = 1 #standered pressure
const Tk_0 = 272.15 #standered temperture
const u = 2.40482555769577 #first order bessel solution http://wwwal.kuicr.kyoto-u.ac.jp/www/accelerator/a4/besselroot.htmlx
const eh = 4.35974434e-18      # Hartree energy
const a_b = 0.52917721092e-10  # Bohr radius
const a_E = 5.14220652e11      # atomic unit of electric field
const a_t = 2.418884326502e-17 # atomic unit of time

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev


text:
include("test_utils.jl")
include("test_curves.jl")
include("test_roadways.jl")
include("test_roadway_generation.jl")


text:
"""
            Objects for SPH mapping.

    Author: Ludwig Böss
    Contact: lboess@usm.lmu.de
    Created: 2020-07-03

"""

# using Distributed

struct mappingParameters{T}

    x_lim::Vector{T}
    y_lim::Vector{T}
    z_lim::Vector{T}
    center::Vector{T}
    halfsize::Vector{T}
    len2pix::T
    pixelSideLength::T
    Npixels::Vector{Int64}
    boxsize::T
    periodic::Bool

end

"""
    mappingParameters( T::DataType=Float64;
                       x_lim::Vector{Float64}   = [-1.0, -1.0],
                       y_lim::Vector{Float64}   = [-1.0, -1.0],
                       z_lim::Vector{Float64}   = [-1.0, -1.0],
                       center::Vector{Float64}  = [-1.0, -1.0, -1.0],
                       x_size::Float64          =  -1.0,
                       y_size::Float64          =  -1.0,
                       z_size::Float64          =  -1.0,
                       pixelSideLength::Float64 =  -1.0,
                       Npixels::Int64           =   0)

Parameter object for sph to grid mapping. Define either `*_lim`, or `center` and `*_size`.
Resolution is defined by `pixelSideLength` or `Npixels`.
"""
function mappingParameters( T::DataType=Float64;
                            x_lim::Vector{<:Real}   = [-1.0, -1.0],
                            y_lim::Vector{<:Real}   = [-1.0, -1.0],
                            z_lim::Vector{<:Real}   = [-1.0, -1.0],
                            center::Vector{<:Real}  = [-1.0, -1.0, -1.0],
                            x_size::Real            =  -1.0,
                            y_size::Real            =  -1.0,
                            z_size::Real            =  -1.0,
                            pixelSideLength::Real   =  -1.0,
                            Npixels::Integer        =   0,
                            boxsize::Real           =  -1.0)

    # calculate limits if center position and sizes are given
    if ( x_lim == [-1.0, -1.0] && ( y_lim == [-1.0, -1.0] && z_lim == [-1.0, -1.0] ))

        if (center != [-1.0, -1.0, -1.0]) && ( x_size != -1.0 && (y_size != -1.0 && z_size != -1.0) )

            x_lim = [ center[1] - 0.5x_size, center[1] + 0.5x_size ]
            y_lim = [ center[2] - 0.5y_size, center[2] + 0.5y_size ]
            z_lim = [ center[3] - 0.5z_size, center[3] + 0.5z_size ]

        else
            error("Giving a center position requires extent in x, y and z direction.")
        end
    end

    # calculate side lengths from limits
    if x_size == -1.0
        x_size = x_lim[2] - x_lim[1]
    end
    if y_size == -1.0
        y_size = y_lim[2] - y_lim[1]
    end
    if z_size == -1.0
        z_size = z_lim[2] - z_lim[1]
    end

    if center == [-1.0, -1.0, -1.0]
        center[1] = x_lim[1] + 0.5*x_size
        center[2] = y_lim[1] + 0.5*y_size
        center[3] = z_lim[1] + 0.5*z_size
    end

    # find the maximum extent of the map
    max_size = max(x_size, y_size)

    if (pixelSideLength == -1.0) & (Npixels != 0)
        pixelSideLength = max_size/Npixels
    elseif (pixelSideLength != -1.0) & (Npixels == 0)
        Npixels = floor(Int64, max_size/pixelSideLength)
        # recalculate pixelSideLenght to account for rounding
        pixelSideLength = max_size/Npixels
    else
        error("Please specify pixelSideLenght or number of pixels!")
    end

    # area of a pixel in code units
    pixelArea = pixelSideLength^2

    # number of pixels in each dimension
    Npix = [ Npixels, Npixels, Npixels ]

    periodic = false

    if boxsize != -1.0
        periodic = true
    end

    halfsize = 0.5 .* [ x_size, y_size, z_size ]

    len2pix = 1.0/pixelSideLength

    mappingParameters{T}(x_lim, y_lim, z_lim,
                         center,
                         halfsize,
                         len2pix,
                         pixelSideLength,
                         Npix,
                         boxsize, periodic)

end


text:
const OUTPUT_VERSIONS = ["1.04", "1.05"]

"""
# Scanner

Examples
```julia
NMAP.Scanner(
    NMAP.ports("80"),
    NMAP.targets("127.0.0.1"),
    NMAP.timingtemplate(NMAP.paranoid)
)

NMAP.Scanner(
    NMAP.ports("80"),
    NMAP.targets("127.0.0.1"),
    NMAP.timingtemplate(NMAP.paranoid),
    NMAP.os_detection(
        NMAP.osscan_guess()
    )
)

# Mixed option definition
NMAP.Scanner(
    NMAP.fastmode(),
    "-sV",
    NMAP.targets("127.0.0.1"),
)

# Options as strings
NMAP.Scanner("-p1-80", "-sV", "127.0.01")
```
"""
mutable struct Scanner
    cmd::Cmd
    args::Vector{String}

    binpath::String
    stderr
    stdout

    debug::Bool # When true,  will not attempt to generate Scan after completion
    xml::Bool   # When false, will not attempt to generate Scan after completion

    # create kwarg constructor to prevent some constructor and dispatch issues that appeared when using Base.@kwdef
    function Scanner(;
            cmd     = Cmd(``),
            args    = Vector{String}(),
            binpath = "nmap",
            stderr  = Base.stderr,
            stdout  = Base.stdout,
            debug   = false,
            xml     = true
        )
        return new(cmd, args, binpath, stderr, stdout, debug, xml)
    end
end

Base.push!(scanner::Scanner,    arg::String) = Base.push!(scanner.args, arg)
Base.push!(scanner::Scanner,    arg)         = Base.push!(scanner.args, string(arg))
Base.append!(scanner::Scanner,  args::Vector{String}) = append!(scanner.args, args)
Base.push!(scanner::Scanner,    args::Vector{String}) = append!(scanner, args)

function Base.Cmd(scanner::Scanner) ::Cmd
    return Cmd([scanner.binpath, scanner.args...])
end

"""
    function Scanner(options....; binpath::String = "nmap", kwargs...) ::Scanner

Build a Scanner from provided `options`

Input `options` can be strings, functions, Option objects or a mix of both

Examples

```julia
NMAP.Scanner(
    NMAP.ports("80"),
    NMAP.targets("127.0.0.1"),
    NMAP.timingtemplate(NMAP.paranoid))

NMAP.Scanner("-p1-80", "-sV", "127.0.01")

NMAP.Scanner(
    NMAP.fastmode(),
    "-sV",
    NMAP.targets("127.0.0.1"))
```
"""
function Scanner(options...;
    args::Vector{String} = Vector{String}(),
    binpath::String      = "nmap",
    stdout               = Base.stdout,
    stderr               = Base.stderr,
    debug::Bool          = false,
    kwargs...)

    this = Scanner(; args=args, binpath=binpath, stdout=stdout, stderr=stderr, debug=debug)

    for option in options
        if option isa Function || option isa Option
            option(this)
        elseif option isa String
            push!(this.args, option)
        elseif option isa Vector
            append!(this.args, option)
        end
    end
    this.cmd = Cmd(this)

    return this
end

"""
    Option

Options are used to group scanner options.

An Option will apply or remove options from scanner arguments. For some options, it can/will also mutate other scanner properties.

## Examples

### Defining an Option
```julia
opt = Option(function(scanner::Scanner) push!(scanner.args, "-someOpt") end)
```
### Using Option
```julia
opt = Option("-someOpt")

opts = Option("-someOpt", "someValue", anotherValue)

append!(scanner, [opt, opts])
```
"""
mutable struct Option
    fn::Function
end
function Option(option::String)
    return Option(
        function(scan::Scanner)
            push!(scan.args, option)
        end)
end
function Option(options...)
    return Option(
        function(scan::Scanner)
            [push!(scan, option) for option in options]
        end
    )
end

"""
    Apply `option` to `scanner`

Invokes `option.fn(scanner)`
"""
function(option::Option)(scanner::Scanner)
    return option.fn(scanner)
end

"""
    Apply `option` to `scanner`

Invokes `option(scanner)`
"""
function (scanner::Scanner)(option::Option)
    return option(scanner)
end

Base.push!(scanner::Scanner,    option::Option)          = option(scanner)
Base.append!(scanner::Scanner,  options::Vector{Option}) = [x(scanner) for x in options]
Base.push!(scanner::Scanner,    options::Vector{Option}) = append!(scanner, options)

"""
    forcexml(out::String="-") :: Option

Enable XML output in stdout
"""
function _forcexml(out::String = "-") ::Option
    return Option(
        function(scanner::Scanner)
            # Enable XML Output
            push!(scanner.args, "-oX")
            # Write output to output stream
            push!(scanner.args, out)
        end)
end

# Utility Type, as JSON is a Module, not a Type
abstract type Json end
const Sinks = Union{Scan, Dict, Json, String}

Marshalling.unmarshall(::Type{String}, xml::XMLDict.XMLDictElement; kwargs...) = String(XMLDict.dict_xml(xml))
Marshalling.unmarshall(::Type{Dict},   xml::XMLDict.XMLDictElement; kwargs...) = XMLDict.xml_dict(xml)
Marshalling.unmarshall(::Type{Json},   xml::XMLDict.XMLDictElement; kwargs...) = JSON.json(Marshalling.unmarshall(Dict, xml))
Marshalling.unmarshall(::Type{Scan},   xml::XMLDict.XMLDictElement; kwargs...) = Scan(xml)

Marshalling.unmarshall(::Type{String}, xml::String; kwargs...) = xml
Marshalling.unmarshall(::Type{Dict},   xml::String; kwargs...) = Marshalling.unmarshall(Dict, XMLDict.parse_xml(xml))
Marshalling.unmarshall(::Type{Json},   xml::String; kwargs...) = JSON.json(Marshalling.unmarshall(Dict, xml))
Marshalling.unmarshall(::Type{Scan},   xml::String; kwargs...) = Scan(xml)

"""
    run!(scan::Scanner)

run! the `scanner` synchronously and return scan results

Invoking `run!` will mutate the `scanner`

Default behaviour is to return a `Scan type`. This can be configured via `sink` kwarg

WIP: To redirect formatted output to both an output stream and Base.stdout, set `interactive` to false
"""
function run!(
        scanner::Scanner;
        autoremove::Bool             = true,
        interactive::Bool            = true,
        file::Union{Nothing, String} = nothing,
        replace::Bool                = false,
        sink::Type{S}                = Scan) where S<:Sinks

    if scanner.debug || !scanner.xml
        @warn "Scanning in debug mode, scan outputs will not be automatically parsed to a Scan"
    end

    if scanner.xml
        if !interactive
            # deactive interactive output and write results to standard output stream ("-") or given `file`
            path = isnothing(file) ? "-" : file
        else
            # set output format to XML and rediret xml to tmp or given file
            if isnothing(file)
                path, io = mktemp(tempdir(), cleanup=true)
                close(io)
            else
                path = file
            end
        end
        scanner(_forcexml(path))
    end

    # recreate Cmd in case it has been altered since initial creation
    scanner.cmd = Cmd(scanner)
    @info scanner.cmd

    runctx = (stdout = scanner.stdout, stderr = scanner.stderr)
    if !interactive
        runctx = (stdout = IOBuffer(), stderr = scanner.stderr)
    end
    raw = Base.run(pipeline(scanner.cmd; runctx...), wait=true)

    output = interactive ? read(path, String) : String(take!(runctx.stdout))
    if !scanner.debug && scanner.xml
        # build scan output as type `sink`, default is Scan
        output = Marshalling.unmarshall(sink, output)
    end
    if (scanner.xml && path != "-") && autoremove
        @async rm(path)
    end
    return output
end


text:
using DifferentialEquations
using PolynomialRoots
using Printf
# using MATLAB
# using ProgressMeter

include("leap_frog_definitions.jl")
max_hit_q1 = 6
max_hit_p2 = 10000
barrier = 5
H=.25

condition_max_hits(u, t, integrator) = u[3] >max_hit_q1 || u[6] > max_hit_p2 || maximum([abs(u[1]),abs(u[2]),abs(u[4]),abs(u[5])]) > barrier
affect_stop!(integrator) = terminate!(integrator)
function condition_hits_PSS_q1(u, t, integrator) # Event when event_f(u,t) == 0
    u[1]
end

function condition_hits_PSS_p2f(u, t, integrator) # Event when event_f(u,t) == 0
    u[5]+100
end

function affect_update_iterator_q1!(integrator)
    q2=integrator.u[2]
    p2=integrator.u[5]
    p1=integrator.u[4]
    P1=P1_poly(q2,p2, H)
    if !isempty(p1) && abs(p1-P1)<1e-5
        integrator.u[3] = integrator.u[3] + 1
    end
end
function affect_update_iterator_p2!(integrator)
    integrator.u[6] = integrator.u[6] + 1
end

callback_max_hits = DiscreteCallback(condition_max_hits, affect_stop!)
callback_hits_PSS_q1_forward = ContinuousCallback(condition_hits_PSS_q1, nothing,  affect_update_iterator_q1!,rootfind=true)
callback_hits_PSS_q1_backward = ContinuousCallback(condition_hits_PSS_q1,  affect_update_iterator_q1!,nothing,rootfind=true)

callback_hits_PSS_p2_forward = ContinuousCallback(condition_hits_PSS_p2f, affect_update_iterator_p2!)
callback_hits_PSS_p2_backward = ContinuousCallback(condition_hits_PSS_p2f, affect_update_iterator_p2!)

# callback_hits_PSS_q1_forward= ContinuousCallback(condition_hits_PSS_q1, affect_update_iterator_q1!)
# callback_hits_PSS_q1_backward= ContinuousCallback(condition_hits_PSS_q1,affect_update_iterator_q1!)

# callback_hits_PSS_p2_forward = ContinuousCallback(condition_hits_PSS_p2,  affect_update_iterator_p2!)
# callback_hits_PSS_p2_backward = ContinuousCallback(condition_hits_PSS_p2, affect_update_iterator_p2!)

cb_forward = CallbackSet(callback_hits_PSS_q1_forward, callback_hits_PSS_p2_forward, callback_max_hits)
cb_backward = CallbackSet(callback_hits_PSS_q1_backward, callback_hits_PSS_p2_backward, callback_max_hits)

function Q1_trajectories(Q2,P2, t_end, H,flag)
    P1 = P1_poly(Q2, P2, H)
    if ~isempty(P1)
        # println(Q1)
        q0, p0 = [zeros(3) for i in 1:2]
        q0[1] = 0
        q0[2] = Q2
        q0[3] = 0
        p0[1] = P1
        p0[2] = P2
        p0[3] = 0
       #constructor for ODE
        prob_f = HamiltonianProblem{true}(Hamiltonian_Dimer, q0, p0, (0., t_end));
        #solve ode , save_everystep=false is important to prevent sol to include all points, not just event points
        # sol_f = solve(prob_f,  RK4(),maxiters=1e20,callback = cb_forward, save_start = false, save_end = true, save_everystep = flag)
        sol_f = solve(prob_f, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback = cb_forward, save_start = false, save_end = true, save_everystep = flag)

        # sol_f=solve(prob, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb_forward,save_start=false,save_end=true,save_everystep=false)

        prob_b = HamiltonianProblem{true}(Hamiltonian_Dimer_Backwards, q0, p0, (0., t_end));
        #solve ode , save_everystep=false is important to prevent sol to include all points, not just event points
        # sol_b = solve(prob_b,  Vern9(), maxiters = 1e20, reltol = 1e-13, abstol = 1e-15,  callback = cb_backward, save_start = false, save_end = true, save_everystep = flag)
        # sol_b = solve(prob_b,  RK4(), maxiters = 1e20, callback = cb_backward, save_start = false, save_end = true, save_everystep = flag)
        sol_b=solve(prob_b, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb_backward,save_start=false,save_end=true,save_everystep=flag)
        point=[sol_f[1,2],sol_f[2,2],sol_f[4,2],sol_f[5,2]]
        return sol_f.t, sol_f[1,:],sol_f[3,end],sol_b.t, sol_b[1,:],sol_b[3,end]

    else
        return -2
    end

end


text:
export find_path, dist

"""
`find_path(G::ImplicitGraph,s,t)` finds a shortest path from `s` to `t`.
"""
function find_path(G::ImplicitGraph{T}, s::T, t::T) where {T}

    if !has(G, s) || !has(G, t)
        error("Source and/or target vertex is not in this graph")
    end
    if s == t
        return [s]
    end

    # set up a queue for vertex exploration
    Q = Queue{T}()
    enqueue!(Q, s)

    # set up trace-back dictionary
    tracer = Dict{T,T}()
    tracer[s] = s

    while length(Q) > 0
        v = dequeue!(Q)
        Nv = G[v]
        for w in Nv
            if haskey(tracer, w)
                continue
            end
            tracer[w] = v
            enqueue!(Q, w)

            if w == t  # success!
                path = Array{T}(undef, 1)
                path[1] = t
                while path[1] != s
                    v = tracer[path[1]]
                    pushfirst!(path, v)
                end
                return path

            end
        end
    end
    return T[]   # return empty array if no path found
end

dist(G::ImplicitGraph{T}, s::T, t::T) where {T} = length(find_path(G, s, t)) - 1


text:
"""
Instead of intrucing fresh variables, DPNeurifyZono relaxes some neurons with a
zonotope.
"""
@with_kw struct DPNeurifyZono <: Solver
    max_zono::Integer = 10
    zono_frac::Float64 = 0.5
    get_zono_idxs = get_fresh_var_idxs_largest_range
end

function forward_linear(solver::DPNeurifyZono, L::LayerNegPosIdx, input::SymbolicIntervalBounds)
    Low, Up = interval_map(L.W_neg, L.W_pos, input.sym.Low, input.sym.Up)
    Low[:, end] .+= L.bias
    Up[:, end] .+= L.bias
    sym = SymbolicInterval(Low, Up, domain(input))
    return SymbolicIntervalBounds(sym, input.lbs, input.ubs)
end

function forward_act(solver::DPNeurifyZono, L::LayerNegPosIdx{ReLU}, input::SymbolicIntervalBounds)
    n = n_nodes(L)
    current_n_vars = size(input.sym.Low, 2) - 1  # TODO: not really n_vars, but n_sym, add field in SymbolicInterval!

    low_low = lower_bounds(input.sym.Low, domain(input), input.lbs[L.index], input.ubs[L.index])
    low_up = upper_bounds(input.sym.Low, domain(input), low_low, input.ubs[L.index])
    up_up = upper_bounds(input.sym.Up, domain(input), low_up, input.ubs[L.index])
    up_low = lower_bounds(input.sym.Up, domain(input), low_low, up_up)

    zono_idxs = solver.get_zono_idxs(solver.max_zono, current_n_vars, low_low, up_up, solver.zono_frac)
    n_zonos = length(zono_idxs)

    λ_l = relaxed_relu_gradient_lower.(low_low, low_up)
    λ_u = relaxed_relu_gradient.(up_low, up_up)
    β_u = -up_low

    n, m = size(input.sym.Low)

    Low = λ_l .* input.sym.Low
    Up  = λ_u .* input.sym.Up
    Up[:, end] .+= λ_u .* max.(β_u, 0)

    G = zeros(n, n_zonos)  # generator matrix for new error terms
    for (i, z) in enumerate(zono_idxs)
        λ = relaxed_relu_gradient(low_low[z], up_up[z])
        β = -low_low[z]
        Low[z,:] = λ * input.sym.Low[z,:]
        Up[z,:] = λ * input.sym.Low[z,:]
        G[z, i] = λ * max(β, 0)  # we can remove the max, as we only introduce zonos for non-fixed ReLUs
    end

    # add all the error terms ϵᵢ ∈ [0, 1] to the input domain
    dom = Hyperrectangle(low=[low(domain(input)); zeros(n_zonos)], high=[high(domain(input)); ones(n_zonos)])
    # add the generator matrix for the error terms to the symbolic bounds
    sym = SymbolicInterval([Low[:, 1:m-1] G Low[:, end]], [Up[:, 1:m-1] G Up[:, end]], dom)

    output = SymbolicIntervalBounds(sym, input.lbs, input.ubs)
    output.lbs[L.index] .= low_low
    output.ubs[L.index] .= up_up
    return output
end

function forward_act(solver::DPNeurifyZono, L::LayerNegPosIdx{Id}, input::SymbolicIntervalBounds)
    low_low = lower_bounds(input.sym.Low, domain(input), input.lbs[L.index], input.ubs[L.index])
    up_up = upper_bounds(input.sym.Up, domain(input), low_low, input.ubs[L.index])
    input.lbs[L.index] .= low_low
    input.ubs[L.index] .= up_up

    return input
end


text:
using QPSReader, QuadraticModels, QuadraticModelsGurobi
using Test

@testset "QuadraticModelsGurobi.jl" begin
    qps1 = readqps("QAFIRO.SIF") #lower bounds
    stats1 = gurobi(QuadraticModel(qps1))
    @test isapprox(stats1.objective, -1.59078179, atol=1e-2)

    qps2 = readqps("HS21.SIF") # low/upp bounds
    stats2 = gurobi(QuadraticModel(qps2))
    @test isapprox(stats2.objective, -9.99599999e1, atol=1e-2)

    qps3 = readqps("HS52.SIF") # free bounds
    stats3 = gurobi(QuadraticModel(qps3))
    @test isapprox(stats3.objective, 5.32664756, atol=1e-2)
end


text:
# Lebesgue measure

export Lebesgue

struct Lebesgue{T} <: AbstractMeasure
    support::T
end

Lebesgue() = Lebesgue(ℝ)

gentype(::Lebesgue) = Float64

testvalue(d::Lebesgue) = testvalue(d.support)

proxy(d::Lebesgue) = restrict(in(d.support), LebesgueMeasure())

Base.:∘(::typeof(basemeasure), ::Type{Lebesgue}) = LebesgueMeasure()

Base.show(io::IO, d::Lebesgue) = print("Lebesgue(",d.support,")")

insupport(μ::Lebesgue, x) = x ∈ μ.support

insupport(::Lebesgue{RealNumbers}) = Returns(true)


text:
using Primes

struct primePower
    prime::Integer
    power::Int
end
function maxpow(n::Int, p::Int)
    pow = 0
    while true
        rem = n%p
        if rem != 0
            return pow
        end
        pow+=1
        n /=p
    end
end

oddq(n) =  maxpow(n,2)==0

function evenpart(n::Integer)
    ep = 1
    while true
        rem = n%2
        if rem != 0
            return ep
        end
        ep *= 2
        n /= 2
    end
end

function jacobisymbol(a::Integer, b::Integer)
    mp = maxpow(a, 2)
    if mp%2 != 0
        ared =div(a, 2^mp)
        jtemp = jacobisymbol(ared, b)
        b8 = b%8
        if b8 == 1 || b8 == 7
            return jtemp
        else
            return -jtemp
        end
    end
    if a == b-1
        b4 = b%4
        return b4==1?1:3
    end
    if a==1
        return 1
    end
    if a>b
        return jacobisymbol(a%b, b)
    end
    a4 = a%4
    b4 = b%4
    jtemp = jacobisymbol(b, a)
    if a4 == 3 && b4 ==3
         return -jtemp
    else
        return jtemp
    end
end

max_stars_repo_name: igorrivin/juliant

max_issues_repo_name: igorrivin/juliant

max_forks_repo_name: igorrivin/juliant