text:
module PBURF

import Polynomials
using Polynomials
import PyPlot;
#using PyPlot
import Colors
#using Colors

function homogeneousNormalization(
        twotuple::Tuple{Complex{Float64},Complex{Float64}},
        aproxprecision::Int64=15)
    tt1=complex(twotuple[1])
    tt2=complex(twotuple[2])
    if (abs(tt1)<1.0/10^aproxprecision && abs(tt2)<1.0/10^aproxprecision)
        hpoint=(complex(0.0),complex(0.0))
    else
        if abs(tt1)<= abs(tt2)
        hpoint=(tt1/tt2,complex(1.0))
        else
        hpoint=(complex(1.0),tt2/tt1)
        end
    end
    return hpoint
end

#Example:
homogeneousNormalization((1.0+0.0*im, 0.0+1.0*im))

function sphereBijection(twotuple::Tuple{Complex{Float64},Complex{Float64}},
                        aproxprecision::Int64=8)
    z=twotuple[1]
    t=twotuple[2]
    if (abs(z)<1.0/10^aproxprecision && abs(t)<1.0/10^aproxprecision)
    point=[0.0,0.0,0.0]
    else
    point=[real((conj(z)*t + conj(t)*z)/(conj(t)*t + conj(z)*z)),
            real((1im*(conj(z)*t - conj(t)*z))/(conj(t)*t + conj(z)*z)),
            real((-conj(t)*t + conj(z)*z)/(conj(t)*t + conj(z)*z))]
    end
    return point
end

sphereBijection((1.0+0.0*im, 0.0+1.0*im),9)

function chordalMetric(twotuple::Tuple{Complex{Float64},Complex{Float64}},
                        twotuple1::Tuple{Complex{Float64},Complex{Float64}},
                        aproxprecision::Int64=8)
    norma=vector->sqrt(vector[1]^2+vector[2]^2+vector[3]^2)
    return norma((sphereBijection(twotuple,aproxprecision)-
    sphereBijection(twotuple1,aproxprecision)))
end

chordalMetric((1.0+0.0*im, 0.0+1.0*im),(0.0+0.0*im, 1.0+0.0*im),9)

function bivariatepolyfunction(coefficientlist::Array{Complex{Float64},1},
     u::Complex{Float64}, t::Complex{Float64},d::Int64)
    ff=coefficientlist[1]*t^d
    for i in 2:length(coefficientlist)
        ff=ff+coefficientlist[i]*u^(i-1)*t^(d-i+1)
        +i
    end
    return ff
end

function bivariatepolyfunction(coefficientlist::Array{T,1},
    u::T, t::T,d::Int64) where {T<:Number}
    ff=coefficientlist[1]*t^d
    for i in 2:length(coefficientlist)
        ff=ff+coefficientlist[i]*u^(i-1)*t^(d-i+1)
        +i
    end
    return complex(ff)
end

function bivariatepolyfunction(coefficientlist::Array{T,1},
         u::Complex{Float64},
        t::Complex{Float64},d::Int64) where {T<:Number}
    coefficientlistcomplex=complex(coefficientlist)
    ff=coefficientlistcomplex[1]*t^d
    for i in 2:length(coefficientlist)
        ff=ff+coefficientlist[i]*u^(i-1)*t^(d-i+1)
    +i
    end
    return ff
end

function pairoffunctions(coefficientlistnum::Array{T,1},
    coefficientlistden::Array{T,1}) where {T<:Number}
    ln=length(coefficientlistnum)-1
    ld=length(coefficientlistden)-1
    d=max(ln,ld)
    fff(u::Complex{Float64}, t::Complex{Float64})=
    bivariatepolyfunction(coefficientlistnum,u,t,d)
    ggg(u::Complex{Float64}, t::Complex{Float64})=
    bivariatepolyfunction(coefficientlistden,u,t,d)
    return fff, ggg
end

coefficientlistnumcompila=complex([1.,0.,0.,2.])
coefficientlistdencompila=complex([0.,0.,3.,0.])

bivarej=bivariatepolyfunction(coefficientlistnumcompila,
    1.0+im* 2.3, 2.0+im*0.9,8)

hpaircompila=pairoffunctions(coefficientlistnumcompila,
    coefficientlistdencompila)

function rationalFunction(hpair::Tuple{Function,Function},
    twotuple::Tuple{Complex{Float64},Complex{Float64}},
    aproxprecision::Int64=8)
    c=twotuple[1]
    d=twotuple[2]
    F=hpair[1]
    G=hpair[2]
    cnew=F(c,d)
    dnew=G(c,d)
    hresult=homogeneousNormalization((cnew,dnew),aproxprecision)
    return  hresult
end

rationalFunction(hpaircompila,(1.0+0.0*im, 0.0+1.0*im),9)

function fixedPointsofaIrreduciblePair(
        coefficientlistnum::Union{Vector{Float64},Vector{Complex{Float64}}},
        coefficientlistden::Union{Vector{Float64},Vector{Complex{Float64}}},
        aproxprecision::Int64=8)
    coefficientlistnumcomplex=complex(coefficientlistnum)
    coefficientlistdencomplex=complex(coefficientlistden)
    if length(coefficientlistnum)!=length(coefficientlistden)
    println("The length of the first list
            is different of the length of the second")
    elseif  length(coefficientlistnum)==1 && length(coefficientlistden)==1
    finallistoffixedpoint=homogeneousNormalization((coefficientlistnumcomplex[1],
    coefficientlistdencomplex[1]),
    aproxprecision)
    return [finallistoffixedpoint]
    elseif abs(coefficientlistnum[1])==0.0 &&
        coefficientlistden[length(coefficientlistden)]==0.0 &&
        sum(i->abs(coefficientlistnum[i]-coefficientlistden[i-1]),
        2:length(coefficientlistden))==0.0
    println("All the points are fixed points")
    else
    ffgg=pairoffunctions(coefficientlistnumcomplex, coefficientlistdencomplex)
    Polynumerator=Poly(coefficientlistnumcomplex)
    Polydenominator=Poly(coefficientlistdencomplex)
    Polyx=Poly([0.0+ 0.0* im,1.0+0.0*im])
    lookingzeros= Polynumerator- Polydenominator * Polyx
    fix=roots(lookingzeros)
    le=length(fix)
        fixed=[(complex(fix[i]),1.0+0.0*im) for i in 1:le]
        if abs(ffgg[1](complex(0.0),complex(0.0)))==0.0 &&
            abs(ffgg[2](complex(0.0),complex(0.0)))==0.0
       newfixed=[(complex(0.0),complex(0.0))]
        else
       newfixed=[]
        end
        if abs(ffgg[1](complex(1.0),complex(0.0)))!=0.0 &&
            abs(ffgg[2](complex(1.0),complex(0.0)))==0.0
       morenewfixed=[(complex(1.0),complex(0.0))]
        else
       morenewfixed=[]
        end
        return newfixed, morenewfixed, fixed
    end
end

fixedPointsofaIrreduciblePairex=
    fixedPointsofaIrreduciblePair(coefficientlistnumcompila,
    coefficientlistdencompila,9)[3]

function newstep(hpair::Tuple{Function,Function},
                    iter::Int64,
                    iterprecision::Int64,
                    hpoint::Tuple{Complex{Float64},Complex{Float64}})
    point = hpoint
    tol=1.0/10^(iterprecision)
    number = 0
    imagepoint = rationalFunction(hpair,point,iterprecision)
    while
    (chordalMetric(point, imagepoint,iterprecision) > tol)&& (number < iter)
    point = imagepoint
    imagepoint = rationalFunction(hpair,point,iterprecision)
    number=number+1
    end
    return [imagepoint,number]
end

newstep(hpaircompila, 11, 3, (1.0+0.0*im, 0.0+1.0*im))

function rectangle(xinterval::Tuple{Float64,Float64}=(-1.5,1.5),
                    yinterval::Tuple{Float64,Float64}=(-1.5,1.5),
                    expprecision=10)
    tol=1.0/(2^expprecision)
    a=xinterval[1]
    b=xinterval[2]
    c=yinterval[1]
    d=yinterval[2]
    red=[(complex(r,i),complex(1.0)) for i=d:-tol:c, r=a:tol:b]
    return red
end

function newstep(hpair::Tuple{Function,Function}, iter::Int64,
        iterprecision::Int64,
        hrectangle::Array{Tuple{Complex{Float64},Complex{Float64}},2})
    size1=size(hrectangle)
    result=[newstep(hpair, iter, iterprecision, hrectangle[i,j])
        for i=1:size1[1], j=1:size1[2]]
    return result
end

#Example
rectanglecompila=rectangle((0.0,1.0),(0.0,1.0),1)
newstep(hpaircompila, 11, 3, rectanglecompila)

function positionuptotolerance(fixedPointList::Array{Tuple{Complex{Float64},
                            Complex{Float64}},1},
                            aproxprecision::Int64,
                            twotuple::Tuple{Complex{Float64},Complex{Float64}})
    pos=0
    it=1
    le = length(fixedPointList)
    tol=1/10^(aproxprecision)
    while (it < le+1)
        if (chordalMetric(twotuple, fixedPointList[it],aproxprecision) < tol)
            pos = it
        end
    it=it+1
    end
    return convert(Int64,pos)
end

positionuptotolerance(fixedPointsofaIrreduciblePairex,3,(-0.5 - 0.28867513459481275im,1.0 + 0.0im) )

function position_iteration_upto_tolerances(hpair::Tuple{Function,Function},
        fixedPointList::Array{Tuple{Complex{Float64},Complex{Float64}},1},
        iter::Int64, iterprecision::Int64,
        aproxprecision::Int64,
        twotuple::Tuple{Complex{Float64},Complex{Float64}})
    endpoint_iterations=newstep(hpair, iter, iterprecision, twotuple)
    endpoint=endpoint_iterations[1]
    iterations=convert(Int64,endpoint_iterations[2])
    pos=positionuptotolerance(fixedPointList, aproxprecision, endpoint)
    return (pos,   iterations)
end

position_iteration_upto_tolerances(hpaircompila,fixedPointsofaIrreduciblePairex, 11, 3, 2, (1.0+0.0*im, 0.0+1.0*im))

function position_iteration_upto_tolerances(hpair::Tuple{Function,Function},
        fixedPointList::Array{Tuple{Complex{Float64},Complex{Float64}},1},
        iter::Int64, iterprecision::Int64, aproxprecision::Int64,
        hrectangle::Array{Tuple{Complex{Float64},Complex{Float64}},2})
    size1=size(hrectangle)
    result=[position_iteration_upto_tolerances(hpair,fixedPointList, iter,
        iterprecision, aproxprecision,
        hrectangle[i,j]) for i=1:size1[1], j=1:size1[2]]
    return result
end

position_iteration_upto_tolerances(hpaircompila,fixedPointsofaIrreduciblePairex,11, 3, 2, rectanglecompila)

function fixedPointListex_matrixpositioninterations_RationalFunction(
        coefficientlistnum::Array{T,1},
        coefficientlistden::Array{T,1},
        preexpresolution::Int=8,
        preiteration_max::Int=25,
        preiterprecision::Int64=3,
        preaproxprecision::Int64=3,
        premodeldomain::AbstractString="localrectangle",
        prerectanglesidesdomain::Tuple{Float64,Float64,Float64,Float64}=
        (-1.5,1.5,-1.5,1.5)) where {T<:Number}
    xinterv=(prerectanglesidesdomain[1],prerectanglesidesdomain[2])
    yinterv=(prerectanglesidesdomain[3],prerectanglesidesdomain[4])
    rect=rectangle(xinterv,  yinterv, preexpresolution)
        hpair=pairoffunctions(coefficientlistnum, coefficientlistden)
        coefficientlistnumcomplex=complex(coefficientlistnum)
        coefficientlistdencomplex=complex(coefficientlistden)
    fixedPointListexcomplete=
        fixedPointsofaIrreduciblePair(coefficientlistnumcomplex,
        coefficientlistdencomplex,
        preaproxprecision)
    fixedPointListexcompletea=fixedPointListexcomplete[1]
    if length(fixedPointListexcomplete[2])==0
    fixedPointListexcompleteb=fixedPointListexcompletea
    else
    fixedPointListexcompleteb=union(fixedPointListexcompletea,
        [(complex(1.0),complex(0.0))])
    end
    fixedPointListexample=
    union(fixedPointListexcompleteb,fixedPointListexcomplete[3])
    positer=position_iteration_upto_tolerances(hpair,fixedPointListexample,
    preiteration_max, preiterprecision,
    preaproxprecision, rect)
    return   fixedPointListexample, positer, fixedPointListexcomplete
end

fmf=fixedPointListex_matrixpositioninterations_RationalFunction(
    coefficientlistnumcompila,
    coefficientlistdencompila,3)

function plot_matrixpositer_RationalFunction(
    fixedpointsmatrixpositer,iteration_max::Int=25,
    thecolorstrategy::AbstractString="positionplusiteration",
    themodel::AbstractString="localrectangle",
    therectanglesides::Tuple{Float64,Float64,Float64,Float64}=
    (-1.5,1.5,-1.5,1.5))
    numberoffixedPointListex=length(fixedpointsmatrixpositer[1])
    positer=fixedpointsmatrixpositer[2]
    ab=size(positer)
        if thecolorstrategy=="positionplusiteration"
        numberofcolors=numberoffixedPointListex+1
        integermatrix=
        [(positer[i,j][1]+(1.0-(positer[i,j][2])/iteration_max))+1.0  for i=1:ab[1], j=1:ab[2]];
        integermatrix[1,1]=1;integermatrix[ab[1],ab[2]]=numberofcolors-2
        elseif thecolorstrategy=="iteration"
        numberofcolors=iteration_max+1
        integermatrix=[positer[i,j][2]+0.25 for i=1:size(positer)[1],
            j=1:size(positer)[2]];integermatrix[1,1]=iteration_max
        else thecolorstrategy=="positionfixedpoints"
        numberofcolors=numberoffixedPointListex+1
        integermatrix=[positer[i,j][1]+0.25 for i=1:size(positer)[1],
            j=1:size(positer)[2]]; integermatrix[1,1]=0;integermatrix[ab[1],
            ab[2]]=numberoffixedPointListex
        end
    isinfinity=length(fixedpointsmatrixpositer[3][2])
    if  isinfinity==0
        seed1=[Colors.RGB(0.0,0.0,0.0),Colors.RGB(0.5,0.5,0.5)]
    else
        seed1=[Colors.RGB(0.0,0.0,0.0),Colors.RGB(0.5,0.5,0.5),Colors.RGB(1.0,1.0,0.0)]
    end
    seed=union(seed1,[Colors.RGB(1.0,0.0,0.0),Colors.RGB(0.0,1.0,0.0),Colors.RGB(0.0,0.0,1.0)])
    gseed=Colors.distinguishable_colors(numberofcolors, seed)
    fpcm=PyPlot.ColorMap(gseed)
    img=PyPlot.imshow(integermatrix, cmap=fpcm, extent=[therectanglesides[1],
        therectanglesides[2],
        therectanglesides[3],therectanglesides[4]])
    return  fixedpointsmatrixpositer[3],gseed, img
end

function aux_plottingBasinsUnivariateRationalFunctions(
    coefficientlistnum::Array{T,1},
    coefficientlistden::Array{T,1},
    expresolution::Int=8,
    iterationmax::Int=25,
    iterprecision::Int64=3,
    aproxprecision::Int64=3;
    colorstrategy::AbstractString="positionplusiteration",
    model::AbstractString="localrectangle",
    rectanglesides::Tuple{Float64,Float64,Float64,Float64}=
        (-1.5,1.5,-1.5,1.5)) where {T<:Number}
    maxdegree=max(length(coefficientlistnum),length(coefficientlistden))
    samedegreecoefficientlistnum=zeros(Complex{Float64}, maxdegree)
    for i in 1:length(coefficientlistnum)
        samedegreecoefficientlistnum[i]=coefficientlistnum[i]
    end
    samedegreecoefficientlistden=zeros(Complex{Float64}, maxdegree)
    for i in 1:length(coefficientlistden)
        samedegreecoefficientlistden[i]=coefficientlistden[i]
    end
    preexpresolution=expresolution
    preiterationmax=iterationmax
    preiterprecision=iterprecision
    preaproxprecision=aproxprecision
    premodeldomain=model
    prerectanglesidesdomain=rectanglesides
    fixedpointsmatrixpairpositerex=
    fixedPointListex_matrixpositioninterations_RationalFunction(
    samedegreecoefficientlistnum,
    samedegreecoefficientlistden, preexpresolution,preiterationmax,
    preiterprecision,preaproxprecision,
     premodeldomain,  prerectanglesidesdomain)
    fixedpointsmatrixpairpositer3seedimg=
    plot_matrixpositer_RationalFunction(fixedpointsmatrixpairpositerex,
                                        iterationmax,
                                        colorstrategy,
                                        model,
                                        rectanglesides)
    return  fixedpointsmatrixpairpositer3seedimg
end

aux_result=aux_plottingBasinsUnivariateRationalFunctions(coefficientlistnumcompila,
    coefficientlistdencompila)

function plottingBasinsUnivariateRationalFunctions(
    coefficientlistnum::Array{T,1},
    coefficientlistden::Array{T,1},
    expresolution::Int=8,
    iterationmax::Int=25,
    iterprecision::Int64=3,
    aproxprecision::Int64=3;
    colorstrategy::AbstractString="positionplusiteration",
    model::AbstractString="localrectangle",
    rectanglesides::Tuple{Float64,Float64,Float64,Float64}=
        (-1.5,1.5,-1.5,1.5)) where {T<:Number}
    maxdegree=max(length(coefficientlistnum),length(coefficientlistden))
    samedegreecoefficientlistnum=zeros(Complex{Float64}, maxdegree)
    for i in 1:length(coefficientlistnum)
        samedegreecoefficientlistnum[i]=coefficientlistnum[i]
    end
    samedegreecoefficientlistden=zeros(Complex{Float64}, maxdegree)
    for i in 1:length(coefficientlistden)
        samedegreecoefficientlistden[i]=coefficientlistden[i]
    end
    preexpresolution=expresolution
    preiterationmax=iterationmax
    preiterprecision=iterprecision
    preaproxprecision=aproxprecision
    premodeldomain=model
    prerectanglesidesdomain=rectanglesides
    fixedpointsmatrixpairpositerex=
    fixedPointListex_matrixpositioninterations_RationalFunction(
    samedegreecoefficientlistnum,
    samedegreecoefficientlistden, preexpresolution,preiterationmax,
    preiterprecision,preaproxprecision,
     premodeldomain,  prerectanglesidesdomain)
    fixedpointsmatrixpairpositer3seedimg=
    plot_matrixpositer_RationalFunction(fixedpointsmatrixpairpositerex,
                                        iterationmax,
                                        colorstrategy,
                                        model,
                                        rectanglesides)
    return  fixedpointsmatrixpairpositer3seedimg[3]
end

result=plottingBasinsUnivariateRationalFunctions(coefficientlistnumcompila,
    coefficientlistdencompila)
end # module

max_stars_repo_name: luisjavierhernandez/PBURF

max_issues_repo_name: luisjavierhernandez/PBURF

max_forks_repo_name: luisjavierhernandez/PBURF


text:
using Test, LinearAlgebra
using QBase

@testset "./src/constructors/measurements.jl" begin

@testset "mirror_symmetric_qubit_3povm" begin
    @testset "valid input θ = $θ" for θ in π/4:0.1:π/2
        Π = mirror_symmetric_qubit_3povm(θ)
        @test is_povm(Π)
        @test Π isa POVM{Float64}
    end

    @testset "invalid inputs" begin
        @test_throws DomainError mirror_symmetric_qubit_3povm(π/4 - 0.1)
        @test_throws DomainError mirror_symmetric_qubit_3povm(π/2 + 0.1)
    end
end

@testset "asymmetric_qubit_3povm" begin
    @testset "valid range of angles" begin
        @test asymmetric_qubit_3povm(π/2, -0.0001) isa POVM{Float64}
        @test asymmetric_qubit_3povm(0.00001, -π/2) isa POVM{Float64}

        for θ1 in π/24:π/24:π/2
            for θ2 in -π/2:π/24:(θ1-π/2 - 0.00001) # slight offset because θ2 cannot be 0
                Π1 = asymmetric_qubit_3povm(θ1, θ2)
                Π2 = asymmetric_qubit_3povm(θ2, θ1)

                @test Π1 isa POVM{Float64}
                @test Π2 isa POVM{Float64}
            end
        end

        for θ1 in -π/24:-π/24:-π/2
            for θ2 in π/2:-π/24:(π/2+θ1+0.00001)
                Π1 = asymmetric_qubit_3povm(θ1, θ2)
                Π2 = asymmetric_qubit_3povm(θ2, θ1)

                @test Π1 isa POVM{Float64}
                @test Π2 isa POVM{Float64}
            end
        end
    end

    @testset "invalid arguments" begin
        @test_throws DomainError asymmetric_qubit_3povm(π/2, 0)
        @test_throws DomainError asymmetric_qubit_3povm(π/2-0.0001, 0)
        @test_throws DomainError asymmetric_qubit_3povm(0, -π/2)
        @test_throws DomainError asymmetric_qubit_3povm(0, -π/2+0.001)
        @test_throws DomainError asymmetric_qubit_3povm(π/3, π/4)
    end
end

@testset "planar_symmetric_qubit_povm()" begin
    @testset "scanning over simple cases" begin
        for n in 2:50
            Π= planar_symmetric_qubit_povm(n)

            @test length(Π) == n
            @test Π isa POVM{Float64}
        end
    end
end

@testset "trine_qubit_povm()" begin
    trine_Π = trine_qubit_povm()
    @test trine_Π isa POVM{Float64}
    @test isapprox(trine_Π, mirror_symmetric_qubit_3povm(π/3), atol=1e-7)
end

@testset "sic_qubit_povm()" begin
    @test sic_qubit_povm() isa POVM{Complex{Float64}}
end

@testset "sqrt_povm()" begin
    priors = Probabilities([1/3;1/3;1/3])

    states = trine_qubit_states()
    Π = sqrt_povm(priors, states)

    @test Π isa POVM{Float64}
    @test isapprox(Π, trine_qubit_povm(), atol=1e-7)
end

@testset "_naimark_kraus_operators()" begin
    @testset "trine qubit povm" begin
        k = QBase._naimark_kraus_operators(trine_qubit_povm())
        trine_Π = trine_qubit_povm()

        @test k[1]'*k[1] ≈ trine_Π[1]
        @test k[2]'*k[2] ≈ trine_Π[2]
        @test k[3]'*k[3] ≈ trine_Π[3]
    end

    @testset "pentagram qubit povm" begin
        Π = planar_symmetric_qubit_povm(5)
        k = QBase._naimark_kraus_operators(Π)

        @test sum(k_i -> k_i'*k_i, k) ≈ [1 0;0 1]

        @test k[1]'*k[1] ≈ Π[1]
        @test k[2]'*k[2] ≈ Π[2]
        @test k[3]'*k[3] ≈ Π[3]
        @test k[4]'*k[4] ≈ Π[4]
        @test k[5]'*k[5] ≈ Π[5]
    end
end

@testset "naimark_dilation()" begin
    @testset "reconstructs trine probabilites" begin
        projector_dict = naimark_dilation(trine_qubit_povm())

        ρ_set = trine_qubit_states()

        @test 2/3 ≈ tr(Matrix{Complex{Float64}}(
            projector_dict["projectors"][1] * kron(ρ_set[1], projector_dict["ancilla"])
        ))
        @test 2/3 ≈ tr(Matrix{Complex{Float64}}(
            projector_dict["projectors"][2] * kron(ρ_set[2], projector_dict["ancilla"])
        ))
        @test 2/3 ≈ tr(Matrix{Complex{Float64}}(
            projector_dict["projectors"][3] * kron(ρ_set[3], projector_dict["ancilla"])
        ))

        @test 1/6 ≈ tr(Matrix{Complex{Float64}}(
            projector_dict["projectors"][1] * kron(ρ_set[2], projector_dict["ancilla"])
        ))
        @test 1/6 ≈ tr(Matrix{Complex{Float64}}(
            projector_dict["projectors"][2] * kron(ρ_set[1], projector_dict["ancilla"])
        ))
        @test 1/6 ≈ tr(Matrix{Complex{Float64}}(
            projector_dict["projectors"][2] * kron(ρ_set[3], projector_dict["ancilla"])
        ))
    end
end

end

ext: jl


text:
using ApproxTools

using Test
using LinearAlgebra

const BitsFloats = (Float32,Float64)
const Floats = (BitsFloats..., BigFloat)
const Reals = (Int, Floats...)
rnc(reals::Type) = (reals,complex.(reals))
rnc(reals::Tuple) = (reals...,complex.(reals)...)

include("tucker.jl")
include("grideval.jl")
include("LogNumbers.jl")
include("LinearCombinations.jl")
include("approximate.jl")
include("inner_products.jl")
include("points/ChebyshevPoints.jl")
include("points/Midpoints.jl")
include("points/TrapezoidalPoints.jl")
include("bases.jl")


text:
using Polynomials
using LsqFit

using Plots

function main()

    # nonlinear fit stuff
    p0 = [0.5, 0.5]                         # guess
    model(t, p) = p[1] * exp.(-p[2] * t)    # model trying to fit
    model1(t, p) = p[1] * sin.(p[2] * t) + p[3]

    xdata = range(0, stop=10, length=20)
    ydata = model(xdata, [1.0 2.0]) + 0.01*randn(length(xdata))

    # poly_fit(xdata, ydata)
    nlin_fit(model, xdata, ydata, p0)

end

function poly_fit(xdata, ydata)

    pfit1 = fit(xdata, ydata, 1)
    pfit2 = fit(xdata, ydata, 2)
    pfit3 = fit(xdata, ydata, 3)

    xlin = range(xdata[1], xdata[end], length=200)

    # plotting
    scatter(xdata, ydata, markersize=3, legend=:topright, label="data")
    # plot!(xlin, pfit1.(xlin), linestyle=:dash, label="f(x) = x")
    # plot!(xlin, pfit2.(xlin), label="f(x) = x²")
    # plot!(xlin, pfit3.(xlin), label="f(x) = x³", dpi=200)

    pfit1 = fit(xdata[1:5], ydata[1:5], 3)
    pfit2 = fit(xdata[5:end], ydata[5:end], 1)

    plot!(xlin[1:50], pfit1.(xlin[1:50]), label="f(x) = x³")
    plot!(xlin[50:end], pfit2.(xlin[50:end]), label="f(x) = x", dpi=200)

    xaxis!("x")
    yaxis!("y")
    title!("poly fit")

end

function nlin_fit(model, xdata, ydata, p0)

    nlinfit = curve_fit(model, xdata, ydata, p0)
    pfit = nlinfit.param
    print(pfit)
    xlin = range(xdata[1], xdata[end], length=200)

    scatter(xdata, ydata, markersize=3, legend=:topright, label="data")
    plot!(xlin, model(xlin, [p0[1], p0[2]]), label="initial model")
    plot!(xlin, model(xlin, [pfit[1], pfit[2]]), linestyle=:dash, label="fitted model", dpi=200)

    xaxis!("x")
    yaxis!("y")
    title!("nonlinear fit")

end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
# TODO: check effect of polar coordinates on ihomogeneous systems.
# TODO: loop over size of M to apply BC's deeper in the interior in order to do photonic crystal waveguides as open boundaries
# TODO: fix gradient

module DifferentialOperators

using ..Bravais,
..BoundaryConditions,
..CoordinateSystems,
LinearAlgebra,
SparseArrays

import ..BoundaryConditions: get_dim_side, get_dim, get_side, get_bc_type, reorder, apply_args

export laplacian

struct laplacian{DIM,POS} end
"""
    laplacian(N, coordinate_system, Δ, bcs[, bls, k, ka, kb, ind]; kwargs...)

    laplacian{DIM}(N, coordinate_system, Δ, bcs[, bls, k, ka, kb, ind]; kwargs...)
"""
function laplacian(N, coordinate_system::TCS, Δ, bcs, bls=((noBL{1,1}(),noBL{1,2}()),(noBL{2,1}(),noBL{2,2}())), args...; kwargs...) where TCS<:CoordinateSystem
    OD1 = define_op{2,1,TCS}(N,Δ,bcs,bls; kwargs...)
    OD2 = define_op{2,2,TCS}(N,Δ,bcs,bls; kwargs...)
    ∇₁² = diff_op(OD1,args...)
    ∇₂² = diff_op(OD2,args...)
    return sparse(∇₁², ∇₂²), (OD1,OD2)
end
function laplacian{DIM}(N, coordinate_system::TCS, Δ, bcs, bls=((noBL{1,1}(),noBL{1,2}()),(noBL{2,1}(),noBL{2,2}())), args...; kwargs...) where {DIM,TCS<:CoordinateSystem}
    N = DIM==1 ? (N[1],1) : (1,N[2])
    OD = define_op{2,DIM,TCS}(N,Δ,bcs,bls; kwargs...)
    ∇² = diff_op(OD,args...)
    return sparse(∇²), (OD,)
end

##############################################################################################
abstract type AbstractPolarity end
struct Backward<:AbstractPolarity end
struct Forward<:AbstractPolarity end
struct Central<:AbstractPolarity end

struct OperatorDefinition{ORD,DIM,TCS,TBC11,TBC12,TBC21,TBC22,TBL11,TBL12,TBL21,TBL22,TH1,TH2,TF,POL}
    N::NTuple{2,Int}
    Δ::NTuple{2,NTuple{2,Float64}}
    dx::NTuple{2,Float64}
    x::NTuple{2,LinRange{Float64}}
    bcs::Tuple{Tuple{TBC11,TBC12},Tuple{TBC21,TBC22}}
    bls::Tuple{Tuple{TBL11,TBL12},Tuple{TBL21,TBL22}}
    h::Tuple{Array{TH1,1},Array{TH2,1}}
    f::NTuple{2,Array{TF,1}}
    xd::NTuple{2,LinRange{Float64}}
    hd::Tuple{Array{TH1,1},Array{TH2,1}}
    fd::NTuple{2,Array{TF,1}}

    OperatorDefinition{ORD,DIM,CS}(N,Δ,bcs; kwargs...) where {ORD,DIM,CS} = OperatorDefinition{ORD,DIM,CS}(N,_oc_Δ(Δ),_oc_bcs(bcs),((noBL{1,1}(),noBL{1,2}()),(noBL{2,1}(),noBL{2,2}())); kwargs...)
    OperatorDefinition{ORD,DIM,CS}(N,Δ,bcs,bls; kwargs...) where {ORD,DIM,CS} = OperatorDefinition{ORD,DIM,CS}(N,_oc_Δ(Δ),_oc_bcs(bcs),_oc_bls(bls); kwargs...)
    function OperatorDefinition{ORD,DIM,CS}(N,Δ::Tuple{Tuple{Float64,Float64},Tuple{Float64,Float64}},bcs::Tuple{Tuple,Tuple},bls::Tuple{Tuple,Tuple}; kwargs...) where {ORD,DIM,CS,TD}
        bc1_bool = typeof(bcs[1][1])<:FloquetBC && typeof(bcs[1][2])<:FloquetBC
        bc2_bool = typeof(bcs[2][1])<:FloquetBC && typeof(bcs[2][2])<:FloquetBC
        if (bc1_bool || bc2_bool)
            lattice = haskey(kwargs,:lattice) ? kwargs[:lattice] : BravaisLattice()
        else
            lattice = BravaisLattice()
        end
        a, b, α, β = lattice.a, lattice.b, lattice.α, lattice.β

        dx1 = !isinf(a) ? a*sin(β-α)/sin(β)/N[1] : (Δ[1][2]-Δ[1][1])/N[1]
        dx2 = !isinf(b) ?  b*sin(β)/N[2] : (Δ[2][2]-Δ[2][1])/N[2]

        Δ = !isinf(a) ? ((Δ[1][1],Δ[1][1]+a*sin(β-α)/sin(β)),Δ[2]) : Δ
        Δ = !isinf(b) ? (Δ[1],(Δ[2][1],Δ[2][1]+b*sin(β))) : Δ

        xmin = float.((Δ[1][1],Δ[2][1]))
        dx = float.((dx1,dx2))

        x1 = LinRange(Δ[1][1]+dx[1]/2,Δ[1][2]-dx[1]/2,N[1])
        x2 = LinRange(Δ[2][1]+dx[2]/2,Δ[2][2]-dx[2]/2,N[2])
        x = (x1,x2)
        x1 = LinRange(Δ[1][1],Δ[1][2],N[1]+1)
        x2 = LinRange(Δ[2][1],Δ[2][2],N[2]+1)
        xd = (x1,x2)

        # bcs = apply_args(bcs, CS() ; N=N, dx=dx, xmin=xmin, lattice=lattice, kwargs...)
        # depth = haskey(kwargs,:depth) ? kwargs[:depth] : 0
        # bls = apply_args(bls; Δ=Δ, depth=depth)
        h1 = typeof(bls[1])<:Tuple{noBL,noBL} ? zeros(Int,size(x[1])) : 1im*(bls[1][1].(x[1]) + bls[1][2].(x[1]))
        h2 = typeof(bls[2])<:Tuple{noBL,noBL} ? zeros(Int,size(x[2])) : 1im*(bls[2][1].(x[2]) + bls[2][2].(x[2]))
        h = (h1,h2)
        h1 = typeof(bls[1])<:Tuple{noBL,noBL} ? zeros(Int,size(xd[1])) : 1im*(bls[1][1].(xd[1]) + bls[1][2].(xd[1]))
        h2 = typeof(bls[2])<:Tuple{noBL,noBL} ? zeros(Int,size(xd[2])) : 1im*(bls[2][1].(xd[2]) + bls[2][2].(xd[2]))
        hd = (h1,h2)
        f1 = typeof(bls[1])<:Tuple{noBL,noBL} ? zeros(Int,size(x[1])) : 1im*(-reverse(cumsum(bls[1][1].(x[1][end:-1:1]))) + cumsum(bls[1][2].(x[1])))*dx[1]
        f2 = zeros(Int,size(x[2]))
        f = (f1,f2)
        fd1 = typeof(bls[1])<:Tuple{noBL,noBL} ? zeros(Int,size(xd[1])) : 1im*(-reverse(cumsum(bls[1][1].(xd[1][end:-1:1]))) + cumsum(bls[1][2].(xd[1])))*dx[1]
        fd2 = zeros(Int,size(xd[2]))
        fd = (fd1,fd2)

        pol = haskey(kwargs,:polarity) ? typeof(kwargs[:polarity]) : (ORD==1 ? Central : Union{})

        tbc11 = typeof(bcs[1][1])
        tbc12 = typeof(bcs[1][2])
        tbc21 = typeof(bcs[2][1])
        tbc22 = typeof(bcs[2][2])

        tbl11 = typeof(bls[1][1])
        tbl12 = typeof(bls[1][2])
        tbl21 = typeof(bls[2][1])
        tbl22 = typeof(bls[2][2])

        return new{ORD,DIM,CS,tbc11,tbc12,tbc21,tbc22,tbl11,tbl12,tbl21,tbl22,eltype(h1),eltype(h2),eltype(f1),pol}(Tuple(N),Δ,dx,x,bcs,bls,h,f,xd,hd,fd)
    end

    Base.show(io::IO,OD::OperatorDefinition{ORD,DIM,TCS}) where {ORD,DIM,TCS} = print("OperatorDefinition{$ORD,$DIM,$TCS}")
end
_oc_Δ(Δ::Array{T,2}) where T = _oc_Δ(((Δ[1,1],Δ[2,1]),(Δ[1,2],Δ[2,2])))
_oc_Δ(Δ::Array{T,1}) where T = _oc_Δ(((Δ[1],Δ[2]),(Δ[3],Δ[4])))
_oc_Δ(Δ::Tuple) = _oc_Δ(((Δ[1],Δ[2]),(Δ[3],Δ[4])))
_oc_Δ(Δ::Tuple{Tuple,Tuple}) = ((float(Δ[1][1]),float(Δ[1][2])),(float(Δ[2][1]),float(Δ[2][2])))
_oc_Δ(Δ) = Δ

function _oc_bcs(bcs)
    bcs = reorder(bcs)
    return ((bcs[1],bcs[2]),(bcs[3],bcs[4]))
end
_oc_bcs(bcs::Tuple{Tuple,Tuple}) = bcs

_oc_bls(bls::Tuple) = _oc_bls(bls...)
function _oc_bls(bl...)
    bl = length(bl)==0 ? reorder(bl...,noBL{1,1}(),) : reorder(bl)
    bl = typeof(bl[1])<:AbstractBL{1,1} ? (bl...,) : (noBL{1,1}(),bl...)

    bl = length(bl)==1 ? reorder(bl...,noBL{1,2}(),) : reorder(bl)
    bl = typeof(bl[2])<:AbstractBL{1,2} ? (bl...,) : (noBL{1,2}(),bl...)

    bl = length(bl)==2 ? reorder(bl...,noBL{2,1}(),) : reorder(bl)
    bl = typeof(bl[3])<:AbstractBL{2,1} ? (bl...,) : (noBL{2,1}(),bl...)

    bl = length(bl)==3 ? reorder(bl...,noBL{2,2}(),) : reorder(bl)
    bl = typeof(bl[4])<:AbstractBL{2,2} ? (bl...,) : (noBL{2,2}(),bl...)
    return ((bl[1],bl[2]),(bl[3],bl[4]))
end
_oc_bls(bls::Tuple{Tuple,Tuple}) = bls
struct define_op{ORD,DIM,CS} end

"""
    define_op{ORD,DIM,CS}(args...; kwargs...)
"""
define_op{ORD,DIM,CS}(args...; kwargs...) where {ORD,DIM,CS} = OperatorDefinition{ORD,DIM,CS}(args...;kwargs...)

##################################################################

struct RowColVal{T}
    N::NTuple{2,Int}
    str::NTuple{2,Int}
    rows::Array{Int,1}
    cols::Array{Int,1}
    vals::Array{T,1}

    RowColVal(N,rows,cols,vals) = new{eltype(vals)}(N,Base.size_to_strides(1,N...),rows,cols,vals)
    RowColVal(N) = new{Union{}}(Tuple(N),Base.size_to_strides(1,N...))
    RowColVal(N,vec) = new{eltype(vec)}(N,Base.size_to_strides(1,N...),1:prod(N),1:prod(N),vec)
    function RowColVal(N,vec,dim)
        vec = dim==1 ? repeat(vec;outer=N[2]) : repeat(vec;inner=N[1])
        new{eltype(vec)}(N,Base.size_to_strides(1,N...),1:prod(N),1:prod(N),vec)
    end
    function (rcv::RowColVal{DIM})(i1, j1, i2, j2, vals) where DIM
        rows = 1 .+ rcv.str[1]*(i1.-1) + rcv.str[2]*(j1.-1)
        cols = 1 .+ rcv.str[1]*(i2.-1) + rcv.str[2]*(j2.-1)
        return new{eltype(vals)}(rcv.N,rcv.str,rows,cols,vals)
    end
    RowColVal(N,S::SparseMatrixCSC) = RowColVal(N,findnz(S)...)

    SparseArrays.sparse(rcv::RowColVal) = dropzeros!(sparse(rcv.rows,rcv.cols,rcv.vals,prod(rcv.N),prod(rcv.N),+))

    Base.promote_rule(::Type{RowColVal{T}}, ::Type{RowColVal{U}}) where {T,U} = RowColVal{promote_rule(T,U)}
    Base.convert(::Type{RowColVal{T}},rcv::RowColVal{U}) where {T,U} = RowColVal(rcv.N,rcv.rows,rcv.cols,convert(Array{T,1},rcv.vals))
    Base.convert(::Type{RowColVal{T}},rcv::RowColVal{T}) where T = rcv
end

##################################################################################

struct OperatorConstructor{ORD,DIM,CS,TBL1,TBL2,TBK,TB1,TB2,TOD}
    N::NTuple{2,Int}
    Δ::NTuple{2,NTuple{2,Float64}}
    dx::NTuple{2,Float64}
    bl1::TBL1
    bl2::TBL2
    bulk::RowColVal{TBK}
    bnd1::Array{RowColVal{TB1},1}
    bnd2::Array{RowColVal{TB2},1}
    definition::TOD

    function OperatorConstructor(OD::OperatorDefinition{ORD,DIM,CS,_1,_2,_3,_4,TBL1,TBL2,TBL3,TBL4}) where {ORD,DIM,CS,_1,_2,_3,_4,TBL1,TBL2,TBL3,TBL4}
        bulk, bnd1, bnd2 = RowColVal(OD.N), [RowColVal(OD.N)], [RowColVal(OD.N)]
        DIM==1 ? (tbl1,tbl2)=(TBL1,TBL2) : (tbl1,tbl2)=(TBL3,TBL4)
        return new{ORD,DIM,CS,tbl1,tbl2,Union{},Union{},Union{},typeof(OD)}(OD.N,OD.Δ,OD.dx,OD.bls[DIM][1],OD.bls[DIM][2],bulk,bnd1,bnd2,OD)
    end

    function (OC::OperatorConstructor{2,DIM,CS,_1,_2,Q,_3,_4,TZ})(args::Vararg{Number,M}) where {M,DIM,CS,_1,_2,_3,_4,Q,TZ}
        N,Δ,dx,def = OC.N,OC.Δ,OC.dx,OC.definition

        bulk_rows_on_site = 2:N[DIM]-1
        bulk_rows_nn = 1:N[DIM]-1
        bnd1_rows = [1]
        bnd2_rows = [N[DIM]]
        rows = vcat(bulk_rows_on_site,bulk_rows_nn,bulk_rows_nn.+1,bnd1_rows,bnd2_rows)
        cols = vcat(bulk_rows_on_site,bulk_rows_nn.+1,bulk_rows_nn,bnd1_rows,bnd2_rows)
        vals = get_vals(OC,args...)

        if CS<:Polar && DIM==2
            rf = def.x[1] .+ def.f[1]/args[1]
            rf⁻² = 1 ./rf.^2
            S = spdiagm(0=>rf⁻²)
        else
            S = DIM==1 ? sparse(1.0*I,N[2],N[2]) : sparse(1.0*I,N[1],N[1])
        end
        rows,cols,vals = DIM==1 ? findnz(kron(S,sparse(rows,cols,vals,N[1],N[1],+))) : findnz(kron(sparse(rows,cols,vals,N[2],N[2],+),S))
        rcv = RowColVal(N,rows,cols,vals)
        return new{2,DIM,CS,_1,_2,eltype(vals),_3,_4,TZ}(OC.N,OC.Δ,OC.dx,OC.bl1,OC.bl2,rcv,OC.bnd1,OC.bnd2,OC.definition)
    end

    # different dims
    (OC::OperatorConstructor{ORD,DIM1})(BC::AbstractBC{DIM2},args...) where {ORD,DIM1,DIM2} = OC

    # same dims, side=1
    function (OC::OperatorConstructor{_1,DIM,CS,_3,_4,_5,_6,_7,_8})(BC::AbstractBC{DIM,1},args...) where {_1,DIM,CS,_3,_4,_5,_6,_7,_8}
        if typeof(BC)<:FloquetBC
            @assert length(args)==4 "only $(length(args)) arguments passed. FloquetBC necessitates 4 args: k, ka, kb, ind"
            # ind_flag = args[4]==0 ? false : true
            # ind_flag ? nothing : args[4]=1
            # ind = args[4]
            args[4]==0 ? args=(args[1:3]...,1) : nothing
            ind = args[4]
            i1,j1,i2,j2 = BC.I1[ind],BC.J1[ind],BC.I2[ind],BC.J2[ind]
        else
            i1,j1,i2,j2 = BC.I1,BC.J1,BC.I2,BC.J2
        end

        vals = get_vals(OC,BC,args...)
        rcv = RowColVal(OC.N)
        rcv = rcv.(i1, j1, i2, j2, vals)
        N = OC.N
        if CS<:Polar && DIM==2
            rf = OC.definition.x[1] .+ OC.definition.f[1]/args[1]
            rf⁻² = 1 ./rf.^2
            S = kron(sparse(I,N[2],N[2]),spdiagm(0=>rf⁻²))
            RCV = Array{RowColVal{eltype(S)},1}(undef,length(rcv))
        else
            S = sparse(I,prod(N),prod(N))
            RCV = rcv
        end

        for i ∈ eachindex(rcv)
            RCV[i] = RowColVal(N,S*sparse(rcv[i]))
        end
        return new{_1,DIM,CS,_3,_4,_5,eltype(vals[1]),_7,_8}(OC.N,OC.Δ,OC.dx,OC.bl1,OC.bl2,OC.bulk,RCV,OC.bnd2,OC.definition)
    end
    # same dims, side=2
    function (OC::OperatorConstructor{_1,DIM,CS,_3,_4,_5,_6,_7,_8})(BC::AbstractBC{DIM,2},args...) where {_1,DIM,CS,_3,_4,_5,_6,_7,_8}
        if typeof(BC)<:FloquetBC
            @assert length(args)==4 "only $(length(args)) arguments passed. FloquetBC necessitates 4 args: k, ka, kb, ind"
            args[4]==0 ? args=(args[1:3]...,1) : nothing
            ind = args[4]
            i1,j1,i2,j2 = BC.I1[ind],BC.J1[ind],BC.I2[ind],BC.J2[ind]
        else
            i1,j1,i2,j2 = BC.I1,BC.J1,BC.I2,BC.J2
        end

        vals = get_vals(OC,BC,args...)
        rcv = RowColVal(OC.N)
        rcv = rcv.(i1, j1, i2, j2, vals)

        N = OC.N
        if CS<:Polar && DIM==2
            rf = OC.definition.x[1] .+ OC.definition.f[1]/args[1]
            rf⁻² = 1 ./rf.^2
            S = kron(sparse(I,N[2],N[2]),spdiagm(0=>rf⁻²))
            RCV = Array{RowColVal{eltype(S)},1}(undef,length(rcv))
        else
            S = sparse(I,prod(N),prod(N))
            RCV = rcv
        end
        for i ∈ eachindex(rcv)
            RCV[i] = RowColVal(N,S*sparse(rcv[i]))
        end
        return new{_1,DIM,CS,_3,_4,_5,_6,eltype(vals[1]),_8}(OC.N,OC.Δ,OC.dx,OC.bl1,OC.bl2,OC.bulk,OC.bnd1,RCV,OC.definition)
    end

    Base.show(io::IO,::OperatorConstructor{ORD,DIM,_1,_2,_3,_4,Union{},Union{},Union{}}) where {ORD,DIM,_1,_2,_3,_4} = print("OC{$ORD,$DIM}: 0/3")
    Base.show(io::IO,::OperatorConstructor{ORD,DIM,_1,_2,_3,_4,TBK,Union{},Union{}}) where {ORD,DIM,_1,_2,_3,_4,TBK} = print("OC{$ORD,$DIM}: 1/3")
    Base.show(io::IO,::OperatorConstructor{ORD,DIM,_1,_2,_3,_4,Union{},TB1,Union{}}) where {ORD,DIM,_1,_2,_3,_4,TB1} = print("OC{$ORD,$DIM}: 1/3")
    Base.show(io::IO,::OperatorConstructor{ORD,DIM,_1,_2,_3,_4,Union{},Union{},TB2}) where {ORD,DIM,_1,_2,_3,_4,TB2} = print("OC{$ORD,$DIM}: 1/3")
    Base.show(io::IO,::OperatorConstructor{ORD,DIM,_1,_2,_3,_4,TBK,TB1,Union{}}) where {ORD,DIM,_1,_2,_3,_4,TBK,TB1} = print("OC{$ORD,$DIM}: 2/3")
    Base.show(io::IO,::OperatorConstructor{ORD,DIM,_1,_2,_3,_4,TBK,Union{},TB2}) where {ORD,DIM,_1,_2,_3,_4,TBK,TB2} = print("OC{$ORD,$DIM}: 2/3")
    Base.show(io::IO,::OperatorConstructor{ORD,DIM,_1,_2,_3,_4,Union{},TB1,TB2}) where {ORD,DIM,_1,_2,_3,_4,TB1,TB2} = print("OC{$ORD,$DIM}: 2/3")
    Base.show(io::IO,::OperatorConstructor{ORD,DIM}) where {ORD,DIM} = print("OC{$ORD,$DIM}: 3/3")
end

function get_vals(OC::OperatorConstructor{2,DIM,CS,TBL1,TBL2},args::Vararg{Number,M}) where {M,DIM,CS,TD,TBL1,TBL2}
    N = OC.N

    h = 1 .+ OC.definition.hd[DIM]/args[1]
    vals = (1 ./h)/OC.dx[DIM]^2
    bulk_on_site = -(vals[2:end-2] + vals[3:end-1])
    bnd1_on_site = -[vals[1] + vals[2]]
    bnd2_on_site = -[vals[end-1] + vals[end]]
    bulk_nn = vals[2:N[DIM]]
    return vcat(bulk_on_site,bulk_nn,bulk_nn,bnd1_on_site,bnd2_on_site)
end
function get_vals(OC::OperatorConstructor{2,1,Polar,TBL1,TBL2},args::Vararg{Number,M}) where {M,TD,TBL1,TBL2}
    DIM=1

    N = OC.N

    r = OC.definition.xd[DIM]
    rf = r .+ OC.definition.fd[DIM]/args[1]
    h = 1 .+ OC.definition.hd[DIM]/args[1]

    vals = (rf./h)/OC.dx[DIM]^2
    bulk_on_site = -(vals[2:end-2] + vals[3:end-1])
    bnd1_on_site = -[vals[1] + vals[2]]
    bnd2_on_site = -[vals[end-1] + vals[end]]
    bulk_nn = vals[2:N[DIM]]
    return vcat(bulk_on_site,bulk_nn,bulk_nn,bnd1_on_site,bnd2_on_site)
end

function get_vals(OC::OperatorConstructor{2,DIM,CS,TBL1,TBL2},BC::T,args...) where {DIM,CS,T<:AbstractLocalBC{DIM,SIDE},TD,TBL1,TBL2} where SIDE
    @assert (TBL1<:noBL && TBL2<:noBL) || length(args)>0 "pml requires frequency argument. Try with additional argument"
    weights = BC.weights
    if CS<:Polar && DIM==1
        fd = SIDE==1 ? OC.definition.fd[DIM][1] : OC.definition.fd[DIM][end]
        rf = OC.Δ[DIM][SIDE] .+ fd/args[1]
    else
        rf = 1
    end
    ind = SIDE==1 ? 1 : OC.N[DIM]
    hd = 1 .+ OC.definition.hd[DIM][ind]/args[1]
    return [(rf./hd).*weights[1]/OC.dx[DIM]^2]
end

function get_vals(OC::OperatorConstructor{2,DIM,CS},BC::FloquetBC{DIM,SIDE},k::Number,ka::Number,kb::Number,ind::Int) where {DIM,SIDE,CS}
    @assert !(CS<:Polar && DIM==1) "combination DIM=1, CS=Polar, BC=Periodic not recognized"
    perp_ind = mod1(ind+1,2)
    I1, I2, J1, J2 = Int[], Int[], Int[], Int[]
    vals = ComplexF64[]
    ns = Int[]
    for i ∈ eachindex(BC.shifts[perp_ind])
        I1 = vcat(I1,BC.I1[perp_ind][i])
        I2 = vcat(I2,BC.I2[perp_ind][i])
        J1 = vcat(J1,BC.J1[perp_ind][i])
        J2 = vcat(J2,BC.J2[perp_ind][i])
        vals = vcat(vals,BC.weights[perp_ind][i])
        ns = vcat(ns,fill(BC.shifts[perp_ind][i],length(BC.weights[perp_ind][i])))
    end
    (K,κ) = ind==1 ? (ka,kb) : (kb,ka)
    (R,r) = ind==1 ? (BC.lattice.a,BC.lattice.b) : (BC.lattice.b,BC.lattice.a)
    if !isinf(r)
        ϕ=κ*r
    else
        ϕ=0.0
    end
    idx = 1
    weights = Array{Array{ComplexF64,1},1}(undef,length(BC.shifts[ind]))
    for i ∈ eachindex(BC.shifts[ind])
        rng = idx:(idx+length(BC.weights[ind][i])-1)
        vals[rng] .*= exp.(-1im*ns[rng]*ϕ)
        weights[i] = vals[rng]
        idx += length(BC.weights[ind][i])
    end
    h = SIDE==1 ? 1 .+ OC.definition.hd[DIM][1]/k : 1 .+ OC.definition.hd[DIM][end]/k
    for i ∈ eachindex(weights)
        weights[i] .*= (1/h)/OC.dx[DIM]^2
    end
    return weights
end

function get_vals(OC::OperatorConstructor{2,DIM,CS},BC::MatchedBC{DIM,SIDE},args...) where {DIM,CS,SIDE}
    vals = Array{Array{ComplexF64,1},1}(undef,length(BC.QNs))
    weights = BC.weights
    R = (CS<:Polar && DIM==1) ? OC.Δ[DIM][SIDE] : 1
    for i ∈ eachindex(vals)
        vals[i] = R.*weights[i]/OC.dx[DIM]^2
    end
    return vals
end

#######################################################################################

struct DifferentialOperator{TD,TS}
    N::NTuple{2,Int}
    d::Array{RowColVal{TD},1}
    f::Array{Function,1}
    s::RowColVal{TS}

    DifferentialOperator(d::RowColVal{TD},s::RowColVal{TS},f::TF) where {TD,TS,TF<:Function} = new{TD,TS}(d.N,[d],[f],s)
    DifferentialOperator(d::Array{RowColVal{TD},1},f::Array{Function,1},s::RowColVal{TS}) where {TD,TS} = new{TD,TS}(d[1].N,d,f,s)

    Base.show(io::IO,DO::DifferentialOperator) = print("Differential Operator")

    Base.vcat(D::DifferentialOperator) = D
    function Base.vcat(D1::DifferentialOperator{TD1,TS1},D2::DifferentialOperator{TD2,TS2}) where {TD1,TD2,TS1,TS2}
        @assert D1.N==D2.N "cannot combine Differential Operators with different N, here D1.N=$(D1.N), D2.N=$(D2.N)"
        d = vcat(D1.d,D2.d)
        f = vcat(D1.f,D2.f)
        s = D1.s
        new{promote_type(TD1,TD2),promote_type(TS1,TS2)}(D1.N,d,f,s)
    end
    Base.vcat(D::DifferentialOperator,Ds...) = vcat(D,vcat(Ds...))

    SparseArrays.sparse(DO::DifferentialOperator) = (sparse.(DO.d),DO.f,sparse(DO.s))
    function SparseArrays.sparse(D1::DifferentialOperator,D2::DifferentialOperator)
        d1,_,s1 = sparse(D1)
        d2,_,s2 = sparse(D2)
        T = promote_type(eltype(d1[1]),eltype(d2[1]))
        D = Array{SparseMatrixCSC{T,Int64},1}(undef,length(d1)+length(d2))
        for i ∈ eachindex(d1)
            D[i] = s2*d1[i]
        end
        for i ∈ eachindex(d2)
            D[length(d1) + i] = s1*d2[i]
        end
        return (D,vcat(D1.f,D2.f),s1*s2)
    end
end

"""
    diff_op(definition, args...)
"""
function diff_op(OD::OperatorDefinition, args...)
    if !(typeof(OD.bls)<:NTuple{2,NTuple{2,noBL}})
        @assert length(args)≥1 "must provide at least one argument for PMLs: k"
    elseif length(args)==0
        args = (1,)
    else
        nothing
    end
    OC = OperatorConstructor(OD)
    blk_op = bulk_op(OD,OC,args...)
    bnd_op = boundary_op(OD,OC,args...)
    return vcat(blk_op,bnd_op)
end

"""
    bulk_op(definition, args...)
"""
bulk_op(OD::OperatorDefinition, args...) = bulk_op(OD, OperatorConstructor(OD), args...)
function bulk_op(OD::OperatorDefinition{ORD,DIM,CS}, OC::OperatorConstructor, args...) where {ORD,DIM,CS}
    OC = OC(args...) # populate bulk nodes

    h = 1 .+ OC.definition.h[DIM]/args[1]
    if CS<:Polar && DIM==1
        rf = OC.definition.x[DIM] + OC.definition.f[DIM]/args[1]
        s = RowColVal(OD.N,rf.*h,DIM)
    else
        s = RowColVal(OD.N,h,DIM)
    end

    return DifferentialOperator(OC.bulk,s,(a...)->1)
end

"""
    boundary_op(definition, args...)
"""
boundary_op(OD::OperatorDefinition, args...) = boundary_op(OD,OperatorConstructor(OD),args...)
function boundary_op(OD::OperatorDefinition{ORD,DIM,CS},OC::OperatorConstructor, args...) where {ORD,DIM,CS}
    OC = OC(OD.bcs[DIM][1],args...)
    OC = OC(OD.bcs[DIM][2],args...)
    rcv = vcat(OC.bnd1,OC.bnd2)

    if typeof(OD.bcs[DIM][1])<:FloquetBC
        if args[4]==0
            fs1 = map(x->((k,_...)->OD.bcs[DIM][1](k,args[1],args[2],1)[x]),1:length(OC.bnd1))
            fs2 = map(x->((k,_...)->OD.bcs[DIM][2](k,args[1],args[2],1)[x]),1:length(OC.bnd2))
        elseif args[4]==1
            fs1 = map(x->((k,κ)->OD.bcs[DIM][1](k,κ,args[2],1)[x]),1:length(OC.bnd1))
            fs2 = map(x->((k,κ)->OD.bcs[DIM][2](k,κ,args[2],1)[x]),1:length(OC.bnd2))
        elseif args[4]==2
            fs1 = map(x->((k,κ)->OD.bcs[DIM][1](k,args[3],κ,2)[x]),1:length(OC.bnd1))
            fs2 = map(x->((k,κ)->OD.bcs[DIM][2](k,args[3],κ,2)[x]),1:length(OC.bnd2))
        end
    else
        fs1 = map(x->((k,_...)->OD.bcs[DIM][1](k)[x]),1:length(OC.bnd1))
        fs2 = map(x->((k,_...)->OD.bcs[DIM][2](k)[x]),1:length(OC.bnd2))
    end
    fs = vcat(fs1,fs2)

    h = 1 .+ OC.definition.h[DIM]/args[1]
    if CS<:Polar && DIM==1
        rf = OC.definition.x[DIM] + OC.definition.f[DIM]/args[1]
        s = RowColVal(OD.N,rf.*h,DIM)
    else
        s = RowColVal(OD.N,h,DIM)
    end
    return DifferentialOperator(rcv, fs, s)
end

end # module


text:
"""
    Schedule

Represents a scheduling queue in which events can be scheduled to occur at future time.
The schedule time is the time of the most recent scheduled event.
The schedule time starts from 0.0. The time -1.0 is the time before simulation.
The time after simulation is the time maxiumum Float64.
You can schedule your events at time greater then 0.0, such as 0.0 + eps().

The events are cloned before each simulation step, when an agent is scheduled,
there are no changing in its memory, until the  next simulation step.
This ensures that the agent memory remains unchanged during the simulation step,
and the agents can see the status of the other agents consistent with the current simulation time.

**Constructors**

    Schedule()

Construct `Schedule` with empty events list and initial time 0.0.

"""

mutable struct Schedule
    simstate::SimState
    events::PriorityQueue{Agent,Priority}
    endevents::Dict{Agent,Agent}
    steps::Int64
    time::Float64
    Schedule(simstate::SimState) =
        new(simstate,PriorityQueue{Agent,Priority}(), Dict{Agent,Agent}(), 0, 0.0)
end

"""
    scheduleOnce!(schedule::Schedule, agent::Agent, time::Float64=1.0; ordering::Int=0)

Schedules the `agent` agent to act at the `schedule` at time interval `time`
in the order for the given time  `ordering`.
"""
function scheduleOnce!(schedule::Schedule, agent::Agent, time::Float64=1.0; ordering::Int=0)
    push!(schedule.endevents,agent=>agent)
    enqueue!(schedule.events, agent, Priority(schedule.time+time,0))
end

function stop!(schedule::Schedule,agent::Agent)
    push!(schedule.endevents,agent=>agent)
end

"""
    scheduleRepeating!(schedule::Schedule, agent::Agent, time::Float64=1.0; ordering::Int=0)

Schedules the `agent` agent to repeatedly act at the `schedule` at time interval `time`
in the order for the given time  `ordering`.

"""
function scheduleRepeating!(schedule::Schedule, agent::Agent, time::Float64=1.0; ordering::Int=0)
    enqueue!(schedule.events, agent, Priority(time,ordering))
end

"""
"""
function simtime(schedule::Schedule)
    schedule.time
end

"""
"""
function step!(schedule::Schedule)
    schedule.steps+=1
    if length(schedule.events) == 0 return end
    schedule.time = peek(schedule.events).second.time
    ctime = schedule.time #current simulation time given by the event with small time
    cevents = Vector{Pair{Agent,Priority}}() #current events to be scheduled for this time
    while true
        if(isempty(schedule.events)) break end # no other events in the scheduling queue
        @inbounds event = peek(schedule.events)
        if(event.second.time > schedule.time) break end #no other events for currenttime
        ctime = event.second.time
        @inbounds dequeue!(schedule.events)
        push!(cevents,event)
    end #while loop

    #@sync @distributed
    for e in cevents
        #Double buffering on the event memory
    #    ce = deepcopy(e.first)
        ce = e.first
        ce.step(schedule.simstate,e.first)
        if(!haskey(schedule.endevents,e.first))
            enqueue!(schedule.events, ce, Priority(ctime+1.0, e.second.priority))
        end
    end
    #Reset the stopped events
    schedule.endevents = Dict{Agent,Agent}()
end


text:
using Test, LinearMaps, Aqua
import LinearMaps: FiveArg, ThreeArg

const matrixstyle = VERSION ≥ v"1.3.0-alpha.115" ? FiveArg() : ThreeArg()

const testallocs = VERSION ≥ v"1.4-"

@testset "code quality" begin
    Aqua.test_all(LinearMaps)
end

include("linearmaps.jl")

include("transpose.jl")

include("functionmap.jl")

include("scaledmap.jl")

include("composition.jl")

include("linearcombination.jl")

include("wrappedmap.jl")

include("uniformscalingmap.jl")

include("numbertypes.jl")

include("blockmap.jl")

include("kronecker.jl")

include("conversion.jl")

include("left.jl")

include("fillmap.jl")

if VERSION ≥ v"1.1"
    include("nontradaxes.jl")
end


text:
######################################################################
# DATA:          aminoAcidChainDict
######################################################################

#http://www.uniprot.org/manual/non_std;Selenocysteine (Sec) and pyrrolysine (Pyl)
aminoAcidDict=Dict{AbstractString,AbstractString}()
aminoAcidDict["ALA"]="A";
aminoAcidDict["ARG"]="R";
aminoAcidDict["ASN"]="N";
aminoAcidDict["ASP"]="D";
aminoAcidDict["CYS"]="C";
aminoAcidDict["GLU"]="E";
aminoAcidDict["GLN"]="Q";
aminoAcidDict["GLY"]="G";
aminoAcidDict["HIS"]="H";
aminoAcidDict["ILE"]="I";
aminoAcidDict["LEU"]="L";
aminoAcidDict["LYS"]="K";
aminoAcidDict["MET"]="M";
aminoAcidDict["PHE"]="F";
aminoAcidDict["PRO"]="P";
aminoAcidDict["SER"]="S";
aminoAcidDict["THR"]="T";
aminoAcidDict["TRP"]="W";
aminoAcidDict["TYR"]="Y";
aminoAcidDict["VAL"]="V";
aminoAcidDict["SEC"]="U";
aminoAcidDict["PYL"]="O";

max_stars_repo_name: christophfeinauer/PdbTool

max_issues_repo_name: christophfeinauer/PdbTool

max_forks_repo_name: christophfeinauer/PdbTool


text:
include(joinpath(@__DIR__, "..", "eap_chain.jl"));

const stepout = 500;
outpath(pargs) = "$(pargs["prefix"])_trajectory.csv";

[
  (chain::EAPChain, step::Int, startup::Bool, cleanup::Bool, pargs) -> begin;
    global outfile;
    if startup
      outfile = open(outpath(pargs), "a");
    end
  end,

  (chain::EAPChain, step::Int, startup::Bool, cleanup::Bool, pargs) -> begin;
    global outfile;
    if cleanup
      close(outfile);
    end
  end,

  (chain::EAPChain, step::Int, startup::Bool, cleanup::Bool, pargs) -> begin;
    global outfile;
    if step % stepout == 0
      writedlm(outfile,
               hcat(step, transpose(chain.r),
                    transpose(chain_μ(chain)), chain.U),
               ',');
    end
  end

]

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
############################################################
# judiWavefield ##############################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: June 2017

export judiWavefield, judiWavefieldException, judiDFTwavefield, muteWavefield, dump_wavefield, fft_wavefield

############################################################

mutable struct judiWavefield{vDT<:Number} <: joAbstractLinearOperator{vDT,vDT}
	name::String
	m::Integer
	n::Integer
	info::Info
    dt
	data
end

mutable struct judiWavefieldException <: Exception
	msg :: String
end

############################################################

## outer constructors

"""
judiWavefield
        name::String
        m::Integer
        n::Integer
        info::Info
        dt::Real
        data

Abstract vector for seismic wavefields.

Constructors
============

Construct wavefield vector from an info structure, a cell array of wavefields and the computational \\
time step dt:

    judiWavefield(info, nt, data)

"""
function judiWavefield(info,dt::Real,data::Union{Array, PyCall.PyObject, String};  vDT::DataType=Float32)

	# length of vector
	m = info.n * sum(info.nt)
	n = 1
	dataCell = [vDT.(data) for j=1:info.nsrc]

	return judiWavefield{vDT}("judiWavefield",m,n,info, Float32(dt), dataCell)
end

function judiWavefield(info, dt::Real, data::Union{Array{Any,1}, Array{Array{T, N}, 1}};vDT::DataType=Float32) where {T, N}
	# length of vector
	nsrc = length(data)
	nsrc != info.nsrc && throw("Different number of sources in info ($(info.nsrc)) and data array ($nsrc)")
	m = info.n * sum(info.nt)
	n = 1
	T != Float32 && (data = tof32.(data))
	return judiWavefield{vDT}("judiWavefield",m,n,info, Float32(dt), data)
end

####################################################################
## overloaded Base functions

# conj(jo)
conj(A::judiWavefield{vDT}) where vDT =
	judiWavefield{vDT}("conj("*A.name*")",A.m,A.n,A.info,A.dt,A.data)

# transpose(jo)
transpose(A::judiWavefield{vDT}) where vDT =
	judiWavefield{vDT}(""*A.name*".'",A.n,A.m,A.info,A.dt,A.data)

# adjoint(jo)
adjoint(A::judiWavefield{vDT}) where vDT = conj(transpose(A))

####################################################################

function vcat(ai::Vararg{judiWavefield{avDT}, N}) where {avDT, N}
	N == 1 && (return ai[1])
	N > 2 && (return vcat(ai[1], vcat(ai[2:end]...)))
	a, b = ai
	m = a.m + b.m
	n = 1
	nsrc = a.info.nsrc + b.info.nsrc
	data = Array{Array{avDT, 3}, 1}(undef, nsrc)
	nt = Array{Integer}(undef, nsrc)
	for j=1:a.info.nsrc
		data[j] = a.data[j]
		nt[j] = a.info.nt[j]
	end
	for j=a.info.nsrc+1:nsrc
		data[j] = b.data[j-a.info.nsrc]
		nt[j] = b.info.nt[j-a.info.nsrc]
	end
	info = Info(a.info.n, nsrc, nt)
	return judiWavefield{avDT}("judiWavefield",info.n * sum(info.nt), 1, info, a.dt ,data)
end

function push!(a::judiWavefield{T}, b::judiWavefield{T}) where T
	append!(a.data, b.data)
	a.m += b.m
	a.info.nsrc += b.info.nsrc
end

# DFT operator for wavefields, acts along time dimension
function fft_wavefield(x_in,mode)
	nsrc = x_in.info.nsrc
	nt = size(x_in.data[1], 1)
	if mode==1
		x = similar(x_in, Complex{Float32})
		for i=1:nsrc
			x.data[i] = fft(x_in.data[i], 1)/sqrt(nt)
		end
	elseif mode==-1
		x = similar(x_in, Float32)
		for i=1:nsrc
			x.data[i] = real(ifft(x_in.data[i], 1)) * sqrt(nt)
		end
	end
	return x
end

# Sampling mask to extract wavefields from full vector
subsample(u::judiWavefield,srcnum) = judiWavefield(u.info,u.data[srcnum];vDT=eltype(u))
similar(u::judiWavefield) = 0f0 * u
similar(u::judiWavefield, vDT::DataType) = judiWavefield{vDT}("judiWavefield",u.m, u.n, u.info, u.dt, Array{vDT}.(0f0.*u.data))

# norm
function norm(a::judiWavefield{avDT}, p::Real=2) where avDT
    if p == Inf
        return max([maximum(abs.(a.data[i])) for i=1:a.info.nsrc]...)
    end
    x = 0.f0
    for j=1:a.info.nsrc
        x += Float32(a.dt) * sum(abs.(vec(a.data[j])).^p)
    end
    return x^(1.f0/p)
end

# inner product
function dot(a::judiWavefield{avDT}, b::judiWavefield{bvDT}) where {avDT, bvDT}
	# Dot product for data containers
	size(a) == size(b) || throw(judiWavefieldException("dimension mismatch"))
	dotprod = 0f0
	for j=1:a.info.nsrc
		dotprod += Float32(a.dt) * dot(vec(a.data[j]),vec(b.data[j]))
	end
	return dotprod
end

# abs
function abs(a::judiWavefield{avDT}) where avDT
	b = deepcopy(a)
	for j=1:a.info.nsrc
		b.data[j] = abs.(a.data[j])
	end
	return b
end

function isapprox(x::judiWavefield, y::judiWavefield; rtol::Real=sqrt(eps()), atol::Real=0)
    x.info.nsrc == y.info.nsrc || throw(judiVectorException("Incompatible number of sources"))
    isapprox(x.data, y.data; rtol=rtol, atol=atol)
end

####################################################################################################

isfinite(x::judiWavefield) = all(all(isfinite.(x.data[i])) for i=1:length(x.data))

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
# https://plot.ly/javascript/getting-started

supportedArgs(::PlotlyBackend) = [
    :annotation,
    :background_color, :foreground_color, :color_palette,
    # :background_color_legend, :background_color_inside, :background_color_outside,
    # :foreground_color_legend, :foreground_color_grid, :foreground_color_axis,
    #     :foreground_color_text, :foreground_color_border,
    :group,
    :label,
    :linetype,
    :seriescolor, :seriesalpha,
    :linecolor, :linestyle, :linewidth, :linealpha,
    :markershape, :markercolor, :markersize, :markeralpha,
    :markerstrokewidth, :markerstrokecolor, :markerstrokealpha,
    :fillrange, :fillcolor, :fillalpha,
    :bins,
    :n, :nc, :nr, :layout,
    # :smooth,
    :title, :windowtitle, :show, :size,
    :x, :xlabel, :xlims, :xticks, :xscale, :xflip, :xrotation,
    :y, :ylabel, :ylims, :yticks, :yscale, :yflip, :yrotation,
    :z, :zlabel, :zlims, :zticks, :zscale, :zflip, :zrotation,
    :z,
    :tickfont, :guidefont, :legendfont,
    :grid, :legend, :colorbar,
    :marker_z, :levels,
    :xerror, :yerror,
    :ribbon, :quiver,
    :orientation,
    # :overwrite_figure,
    :polar,
    # :normalize, :weights, :contours, :aspect_ratio
  ]

supportedAxes(::PlotlyBackend) = [:auto, :left]
supportedTypes(::PlotlyBackend) = [:none, :line, :path, :scatter, :steppre, :steppost,
                                   :hist2d, :hist, :density, :bar, :contour, :surface, :path3d, :scatter3d,
                                   :pie, :heatmap] #,, :sticks, :hexbin, :hline, :vline]
supportedStyles(::PlotlyBackend) = [:auto, :solid, :dash, :dot, :dashdot]
supportedMarkers(::PlotlyBackend) = [:none, :auto, :ellipse, :rect, :diamond, :utriangle, :dtriangle, :cross, :xcross,
                                     :pentagon, :hexagon, :octagon, :vline, :hline] #vcat(_allMarkers, Shape)
supportedScales(::PlotlyBackend) = [:identity, :log10] #, :ln, :log2, :log10, :asinh, :sqrt]
subplotSupported(::PlotlyBackend) = true
stringsSupported(::PlotlyBackend) = true

# --------------------------------------------------------------------------------------

function _initialize_backend(::PlotlyBackend; kw...)
  @eval begin
    import JSON
    JSON._print(io::IO, state::JSON.State, dt::Union{Date,DateTime}) = print(io, '"', dt, '"')

    _js_path = Pkg.dir("Plots", "deps", "plotly-latest.min.js")
    _js_code = open(readall, _js_path, "r")

    # borrowed from https://github.com/plotly/plotly.py/blob/2594076e29584ede2d09f2aa40a8a195b3f3fc66/plotly/offline/offline.py#L64-L71 c/o @spencerlyon2
    _js_script = """

    """

    # if we're in IJulia call setupnotebook to load js and css
    if isijulia()
        display("text/html", _js_script)
    end

    # if isatom()
    #     import Atom
    #     Atom.@msg evaljs(_js_code)
    # end

  end
  # TODO: other initialization
end

# ---------------------------------------------------------------------------

function _create_plot(pkg::PlotlyBackend, d::KW)
  # TODO: create the window/canvas/context that is the plot within the backend (call it `o`)
  # TODO: initialize the plot... title, xlabel, bgcolor, etc
  Plot(nothing, pkg, 0, d, KW[])
end

function _add_series(::PlotlyBackend, plt::Plot, d::KW)
  # TODO: add one series to the underlying package
  push!(plt.seriesargs, d)
  plt
end

function _add_annotations{X,Y,V}(plt::Plot{PlotlyBackend}, anns::AVec{@compat(Tuple{X,Y,V})})
  # set or add to the annotation_list
  if haskey(plt.plotargs, :annotation_list)
    append!(plt.plotargs[:annotation_list], anns)
  else
    plt.plotargs[:annotation_list] = anns
  end
end

# ----------------------------------------------------------------

function _before_update_plot(plt::Plot{PlotlyBackend})
end

# TODO: override this to update plot items (title, xlabel, etc) after creation
function _update_plot(plt::Plot{PlotlyBackend}, d::KW)
end

function _update_plot_pos_size(plt::AbstractPlot{PlotlyBackend}, d::KW)
end

# ----------------------------------------------------------------

# accessors for x/y data

# function getxy(plt::Plot{PlotlyBackend}, i::Int)
#   d = plt.seriesargs[i]
#   d[:x], d[:y]
# end
#
# function setxy!{X,Y}(plt::Plot{PlotlyBackend}, xy::Tuple{X,Y}, i::Integer)
#   d = plt.seriesargs[i]
#   d[:x], d[:y] = xy
#   plt
# end

# ----------------------------------------------------------------

function _create_subplot(subplt::Subplot{PlotlyBackend}, isbefore::Bool)
  # TODO: build the underlying Subplot object.  this is where you might layout the panes within a GUI window, for example
  true
end

function _expand_limits(lims, plt::Plot{PlotlyBackend}, isx::Bool)
  # TODO: call expand limits for each plot data
end

function _remove_axis(plt::Plot{PlotlyBackend}, isx::Bool)
  # TODO: if plot is inner subplot, might need to remove ticks or axis labels
end

# ----------------------------------------------------------------

# TODO:
# _plotDefaults[:yrightlabel]       = ""
# _plotDefaults[:xlims]             = :auto
# _plotDefaults[:ylims]             = :auto
# _plotDefaults[:xticks]            = :auto
# _plotDefaults[:yticks]            = :auto
# _plotDefaults[:xscale]            = :identity
# _plotDefaults[:yscale]            = :identity
# _plotDefaults[:xflip]             = false
# _plotDefaults[:yflip]             = false

function plotlyfont(font::Font, color = font.color)
  KW(
      :family => font.family,
      :size   => round(Int, font.pointsize*1.4),
      :color  => webcolor(color),
    )
end

function get_annotation_dict(x, y, val::Union{AbstractString,Symbol})
  KW(
      :text => val,
      :xref => "x",
      :x => x,
      :yref => "y",
      :y => y,
      :showarrow => false,
    )
end

function get_annotation_dict(x, y, ptxt::PlotText)
  merge(get_annotation_dict(x, y, ptxt.str), KW(
      :font => plotlyfont(ptxt.font),
      :xanchor => ptxt.font.halign == :hcenter ? :center : ptxt.font.halign,
      :yanchor => ptxt.font.valign == :vcenter ? :middle : ptxt.font.valign,
      :rotation => ptxt.font.rotation,
    ))
end

function plotlyscale(scale::Symbol)
  if scale == :log10
    "log"
  else
    "-"
  end
end

use_axis_field(ticks) = !(ticks in (nothing, :none))

tickssym(letter) = symbol(letter * "ticks")
limssym(letter) = symbol(letter * "lims")
flipsym(letter) = symbol(letter * "flip")
scalesym(letter) = symbol(letter * "scale")
labelsym(letter) = symbol(letter * "label")
rotationsym(letter) = symbol(letter * "rotation")

function plotlyaxis(d::KW, letter)
  ax = KW(
      :title      => d[labelsym(letter)],
      :showgrid   => d[:grid],
      :zeroline   => false,
    )

  fgcolor = webcolor(d[:foreground_color])
  tsym = tickssym(letter)

  rot = d[rotationsym(letter)]
  if rot != 0
      ax[:tickangle] = rot
  end

  if use_axis_field(d[tsym])
    ax[:titlefont] = plotlyfont(d[:guidefont], fgcolor)
    ax[:type] = plotlyscale(d[scalesym(letter)])
    ax[:tickfont] = plotlyfont(d[:tickfont], fgcolor)
    ax[:tickcolor] = fgcolor
    ax[:linecolor] = fgcolor

    # xlims
    lims = d[limssym(letter)]
    if lims != :auto && limsType(lims) == :limits
      ax[:range] = lims
    end

    # xflip
    if d[flipsym(letter)]
      ax[:autorange] = "reversed"
    end

    # xticks
    ticks = d[tsym]
    if ticks != :auto
      ttype = ticksType(ticks)
      if ttype == :ticks
        ax[:tickmode] = "array"
        ax[:tickvals] = ticks
      elseif ttype == :ticks_and_labels
        ax[:tickmode] = "array"
        ax[:tickvals], ax[:ticktext] = ticks
      end
    end

    ax
  else
    ax[:showticklabels] = false
    ax[:showgrid] = false
  end

  ax
end

# function get_plot_json(plt::Plot{PlotlyBackend})
#   d = plt.plotargs
function plotly_layout(d::KW)
  d_out = KW()

  d_out[:width], d_out[:height] = d[:size]

  bgcolor = webcolor(d[:background_color])
  fgcolor = webcolor(d[:foreground_color])

  # set the fields for the plot
  d_out[:title] = d[:title]
  d_out[:titlefont] = plotlyfont(d[:guidefont], fgcolor)
  d_out[:margin] = KW(:l=>35, :b=>30, :r=>8, :t=>20)
  d_out[:plot_bgcolor] = bgcolor
  d_out[:paper_bgcolor] = bgcolor

  # TODO: x/y axis tick values/labels
  if is3d(d)
      d_out[:scene] = KW(
          :xaxis => plotlyaxis(d, "x"),
          :yaxis => plotlyaxis(d, "y"),
          :xzxis => plotlyaxis(d, "z"),
      )
  else
      d_out[:xaxis] = plotlyaxis(d, "x")
      d_out[:yaxis] = plotlyaxis(d, "y")
  end

  # legend
  d_out[:showlegend] = d[:legend] != :none
  if d[:legend] != :none
    d_out[:legend] = KW(
        :bgcolor  => bgcolor,
        :bordercolor => fgcolor,
        :font     => plotlyfont(d[:legendfont]),
      )
  end

  # annotations
  anns = get(d, :annotation_list, [])
  if !isempty(anns)
    d_out[:annotations] = [get_annotation_dict(ann...) for ann in anns]
  end

  if get(d, :polar, false)
      d_out[:direction] = "counterclockwise"
  end

  d_out
end

function get_plot_json(plt::Plot{PlotlyBackend})
  JSON.json(plotly_layout(plt.plotargs))
end

function plotly_colorscale(grad::ColorGradient, alpha = nothing)
  [[grad.values[i], webcolor(grad.colors[i], alpha)] for i in 1:length(grad.colors)]
end
plotly_colorscale(c, alpha = nothing) = plotly_colorscale(default_gradient(), alpha)

const _plotly_markers = KW(
    :rect       => "square",
    :xcross     => "x",
    :utriangle  => "triangle-up",
    :dtriangle  => "triangle-down",
    :star5      => "star-triangle-up",
    :vline      => "line-ns",
    :hline      => "line-ew",
  )

# get a dictionary representing the series params (d is the Plots-dict, d_out is the Plotly-dict)
function plotly_series(d::KW, plotargs::KW; plot_index = nothing)
    # dumpdict(d,"series",true)
  d_out = KW()

  x, y = collect(d[:x]), collect(d[:y])
  d_out[:name] = d[:label]

  lt = d[:linetype]
  isscatter = lt in (:scatter, :scatter3d)
  hasmarker = isscatter || d[:markershape] != :none
  hasline = !isscatter

  # set the "type"
  if lt in (:line, :path, :scatter, :steppre, :steppost)
    d_out[:type] = "scatter"
    d_out[:mode] = if hasmarker
      hasline ? "lines+markers" : "markers"
    else
      hasline ? "lines" : "none"
    end
    if d[:fillrange] == true || d[:fillrange] == 0
      d_out[:fill] = "tozeroy"
      d_out[:fillcolor] = webcolor(d[:fillcolor], d[:fillalpha])
    elseif !(d[:fillrange] in (false, nothing))
      warn("fillrange ignored... plotly only supports filling to zero. fillrange: $(d[:fillrange])")
    end
    d_out[:x], d_out[:y] = x, y

  elseif lt == :bar
    d_out[:type] = "bar"
    d_out[:x], d_out[:y] = x, y

  elseif lt == :hist2d
    d_out[:type] = "histogram2d"
    d_out[:x], d_out[:y] = x, y
    if isa(d[:bins], Tuple)
      xbins, ybins = d[:bins]
    else
      xbins = ybins = d[:bins]
    end
    d_out[:nbinsx] = xbins
    d_out[:nbinsy] = ybins

  elseif lt in (:hist, :density)
    d_out[:type] = "histogram"
    isvert = isvertical(d)
    d_out[isvert ? :x : :y] = y
    d_out[isvert ? :nbinsx : :nbinsy] = d[:bins]
    if lt == :density
      d_out[:histnorm] = "probability density"
    end

  elseif lt == :heatmap
    d_out[:type] = "heatmap"
    d_out[:x], d_out[:y] = x, y
    d_out[:z] = d[:z].surf
    d_out[:colorscale] = plotly_colorscale(d[:fillcolor], d[:fillalpha])

  elseif lt == :contour
    d_out[:type] = "contour"
    d_out[:x], d_out[:y] = x, y
    d_out[:z] = d[:z].surf
    # d_out[:showscale] = d[:colorbar] != :none
    d_out[:ncontours] = d[:levels]
    d_out[:contours] = KW(:coloring => d[:fillrange] != nothing ? "fill" : "lines")
    d_out[:colorscale] = plotly_colorscale(d[:linecolor], d[:linealpha])

  elseif lt in (:surface, :wireframe)
    d_out[:type] = "surface"
    d_out[:x], d_out[:y] = x, y
    d_out[:z] = d[:z].surf
    d_out[:colorscale] = plotly_colorscale(d[:fillcolor], d[:fillalpha])

  elseif lt == :pie
    d_out[:type] = "pie"
    d_out[:labels] = x
    d_out[:values] = y
    d_out[:hoverinfo] = "label+percent+name"

  elseif lt in (:path3d, :scatter3d)
    d_out[:type] = "scatter3d"
    d_out[:mode] = if hasmarker
      hasline ? "lines+markers" : "markers"
    else
      hasline ? "lines" : "none"
    end
    d_out[:x], d_out[:y] = x, y
    d_out[:z] = collect(d[:z])

  else
    warn("Plotly: linetype $lt isn't supported.")
    return KW()
  end

  # add "marker"
  if hasmarker
    d_out[:marker] = KW(
        :symbol => get(_plotly_markers, d[:markershape], string(d[:markershape])),
        :opacity => d[:markeralpha],
        :size => 2 * d[:markersize],
        :color => webcolor(d[:markercolor], d[:markeralpha]),
        :line => KW(
            :color => webcolor(d[:markerstrokecolor], d[:markerstrokealpha]),
            :width => d[:markerstrokewidth],
          ),
      )

    # gotta hack this (for now?) since plotly can't handle rgba values inside the gradient
    if d[:marker_z] != nothing
      # d_out[:marker][:color] = d[:marker_z]
      # d_out[:marker][:colorscale] = plotly_colorscale(d[:markercolor], d[:markeralpha])
      # d_out[:showscale] = true
      grad = ColorGradient(d[:markercolor], alpha=d[:markeralpha])
      zmin, zmax = extrema(d[:marker_z])
      d_out[:marker][:color] = [webcolor(getColorZ(grad, (zi - zmin) / (zmax - zmin))) for zi in d[:marker_z]]
    end

  end

  # add "line"
  if hasline
    d_out[:line] = KW(
        :color => webcolor(d[:linecolor], d[:linealpha]),
        :width => d[:linewidth],
        :shape => if lt == :steppre
          "vh"
        elseif lt == :steppost
          "hv"
        else
          "linear"
        end,
        :dash => string(d[:linestyle]),
        # :dash => "solid",
      )
  end

  # convert polar plots x/y to theta/radius
  if get(plotargs, :polar, false)
      d_out[:t] = rad2deg(pop!(d_out, :x))
      d_out[:r] = pop!(d_out, :y)
  end

  # # for subplots, we need to add the xaxis/yaxis fields
  # if plot_index != nothing
  #   d_out[:xaxis] = "x$(plot_index)"
  #   d_out[:yaxis] = "y$(plot_index)"
  # end

  d_out
end

# get a list of dictionaries, each representing the series params
function get_series_json(plt::Plot{PlotlyBackend})
  JSON.json(map(d -> plotly_series(d, plt.plotargs), plt.seriesargs))
end

function get_series_json(subplt::Subplot{PlotlyBackend})
  ds = KW[]
  for (i,plt) in enumerate(subplt.plts)
    for d in plt.seriesargs
      push!(ds, plotly_series(d, plt.plotargs, plot_index = i))
    end
  end
  JSON.json(ds)
end

# ----------------------------------------------------------------

function html_head(plt::AbstractPlot{PlotlyBackend})
  ""
end

function html_body(plt::Plot{PlotlyBackend}, style = nothing)
  if style == nothing
    w, h = plt.plotargs[:size]
    style = "width:$(w)px;height:$(h)px;"
  end
  uuid = Base.Random.uuid4()
  html = """

  """
  # @show html
  html
end

function js_body(plt::Plot{PlotlyBackend}, uuid)
    js = """
          PLOT = document.getElementById('$(uuid)');
          Plotly.plot(PLOT, $(get_series_json(plt)), $(get_plot_json(plt)));
    """
end

function html_body(subplt::Subplot{PlotlyBackend})
  w, h = subplt.plts[1].plotargs[:size]
  html = ["
"]
  nr = nrows(subplt.layout)
  ph = h / nr

  for r in 1:nr
    push!(html, "
")

    nc = ncols(subplt.layout, r)
    pw = w / nc

    for c in 1:nc
      plt = subplt[r,c]
      push!(html, html_body(plt, "float:left; width:$(pw)px; height:$(ph)px;"))
    end

    push!(html, "
")
  end
  push!(html, "
")

  join(html)
end

# ----------------------------------------------------------------

function Base.writemime(io::IO, ::MIME"image/png", plt::AbstractPlot{PlotlyBackend})
    writemime_png_from_html(io, plt)
end

function Base.writemime(io::IO, ::MIME"text/html", plt::AbstractPlot{PlotlyBackend})
    write(io, html_head(plt) * html_body(plt))
    # write(io, html_body(plt))
end

function Base.display(::PlotsDisplay, plt::AbstractPlot{PlotlyBackend})
  standalone_html_window(plt)
end

# function Base.display(::PlotsDisplay, plt::Subplot{PlotlyBackend})
#   # TODO: display/show the subplot
# end


text:
abstract type AbstractParser end

struct JSONParser <: AbstractParser end

struct TextParser <: AbstractParser end

struct App{P<:AbstractParser}
    router::HTTP.Router
    reqparser::P
    routelist::Dict{String, Vector{String}}
end

function App(; reqparser::AbstractParser = JSONParser())
    router = HTTP.Router()
    routelist = Dict(
        "GET" => Vector{String}(),
        "POST" => Vector{String}(),
        "PUT" => Vector{String}(),
        "PATCH" => Vector{String}(),
        "DELETE" => Vector{String}(),
        "OPTIONS" => Vector{String}(),
        "HEAD" => Vector{String}()
    )
    return App(router, reqparser, routelist)
end


text:
##################################
## Setting up agents properties ##
##################################
"""
`pick_random_node` function is used to get starting and ending node of agents.
Nodes are randomly chosen from set of provided areas.

**Input parameters**
* `OSMmap` : OpenStreetMapX mapData object with road network data
* `rects` : vector of Rect types interpreted as a set of rectangle areas
"""
function pick_random_node(OSMmap::OpenStreetMapX.MapData,
                          rects::Vector{Rect})
    nodes_in_rects = Vector{Int}()
    for rect in rects
        p1 = ENU(LLA(rect.p1[1], rect.p1[2]), OSMmap.bounds)
        p2 = ENU(LLA(rect.p2[1], rect.p2[2]), OSMmap.bounds)
        exE = extrema([p1.east, p2.east])
        exN = extrema([p1.north, p2.north])
        for key in keys(OSMmap.v)
            if (exE[1] <= OSMmap.nodes[key].east <= exE[2] &&
                exN[1] <= OSMmap.nodes[key].north <= exN[2])
                push!(nodes_in_rects, key)
            end
        end
    end
    unique_nodes = unique!(nodes_in_rects)
    return rand(unique_nodes)
end

"""
`get_nodes_set` function is used to get set of nodes within passed area.

**Input parameters**
* `OSMmap` : OpenStreetMapX mapData object with road network data
* `rects` : vector of Rect types interpreted as a set of rectangle areas
"""
function get_nodes_set(OSMmap::OpenStreetMapX.MapData,
                          rects::Vector{Rect})
    nodes_in_rects = Vector{Int}()
    for rect in rects
        p1 = ENU(LLA(rect.p1[1], rect.p1[2]), OSMmap.bounds)
        p2 = ENU(LLA(rect.p2[1], rect.p2[2]), OSMmap.bounds)
        exE = extrema([p1.east, p2.east])
        exN = extrema([p1.north, p2.north])
        for key in keys(OSMmap.v)
            if (exE[1] <= OSMmap.nodes[key].east <= exE[2] &&
                exN[1] <= OSMmap.nodes[key].north <= exN[2])
                push!(nodes_in_rects, key)
            end
        end
    end
    unique_nodes = unique!(nodes_in_rects)
    return unique_nodes
end

"""
`generate_agents` function creating vector of agents and returning travel time
for initial routes travelled with maximal speed

**Input parameters**
* `OSMmap` : OpenStreetMapX MapData object with road network data
* `N` : number of agents to be generated
* `StartArea` : vector of areas from which agents randomly pick starting point
* `EndArea` : vector of areas from which agents randomly pick ending point
* `α` : percentage of smart agents
* `k` : number of fastest routes used in k-path algorithm
* `T` : control variable for k-path algorithm probability distribution
* `AvgStartTime` : average agents start time in minutes
* `k_routes_dict` : dictionary with multiple shortest paths (values) between vertices (keys)
* `seed` : set seed for random function
"""
function generate_agents(OSMmap::OpenStreetMapX.MapData,
                        N::Int,
                        StartArea::Vector{Rect},
                        EndArea::Vector{Rect},
                        α::Float64,
                        k::Int = 3,
                        T::Float64 = 1.0,
                        AvgStartTime::Float64 = 0.0,
                        k_routes_dict::Dict{Tuple{Int,Int},Array{Vector{Int}}}=
                        Dict{Tuple{Int,Int},Array{Vector{Int}}}())
    #Initialize empty working variables
    AgentsArr = Vector{Agent}()
    #Indicate smart agents
    N_int= Int(ceil(N*α))
    smart_ind = [trues(N_int); falses(N-N_int)]
    #Obtaining nodes set for given starting and ending areas
    start_set = get_nodes_set(OSMmap, StartArea)
    end_set = get_nodes_set(OSMmap, EndArea)
    #Generating N agents
    for i in 1:N
        dist = Inf
        start_node = end_node = counter =  0
        init_route = Array{Int64,1}()
        while dist == Inf
            start_node = rand(start_set)
            end_node = rand(end_set)
            init_route, dist = OpenStreetMapX.fastest_route(OSMmap, start_node, end_node)[1:2]
            counter +=1
            counter == 100 && error("Route from starting to ending point can't be calculated.")
        end
        #First edge in vertices notation
        firstEdge = (OSMmap.v[init_route[1]], OSMmap.v[init_route[2]])
        sTime = rand()*AvgStartTime*60*2
        NewAgent = Agent(smart_ind[i], start_node, end_node,
                        init_route, sTime, 0.0, firstEdge,
                        ifelse(sTime == 0.0,true,false))
        push!(AgentsArr, NewAgent)
    end
    #Get averages speeds from base scenario - simulating short term memory
    AverageSpeeds = simulation_run(:base, OSMmap,AgentsArr; track_avg_speeds=true).AvgSpeeds
    #Recalculate agents paths with k-shortest path algorithm and avg speeds
    max_speeds = OpenStreetMapX.get_velocities(OSMmap)
    for a in AgentsArr
        k_shortest_path_rerouting!(OSMmap, k_routes_dict, a, AverageSpeeds,max_speeds, k, T, 0)
    end
    return AgentsArr
end


text:
module Day25

export day25

function powmod(a::T,n::T,m::T) where T<:Integer
    q = a
    product = T(1)
    while n > T(0)
        n,r = divrem(n,2)
        if r == T(1)
            product = mod(product * q, m)
        end
        q = mod(q^2, m)
    end
    return product
end

function invmod(a::T,m::T) where T<:Integer
    t = gcdx(a,m)
    if t[1] != T(1)
        return error("No inverse exists!")
    end
    return mod(t[2],m)
end

function dlog(g, y, p)
    q = Int(floor(sqrt(p)))
    table = Dict{Int,Int}([1 => 0])
    gq = powmod(g,q,p)
    current_power = gq
    for i in 1:div(p,q)
        table[current_power] = i
        current_power = mod(current_power*gq, p)
    end
    ginv = invmod(g,p)
    current_key = y
    for j in 0:q-1
        if haskey(table, current_key)
            return table[current_key] * q + j
        end
        current_key = mod(current_key*ginv, p)
    end
end

function day25()
    card = 10604480
    door = 4126658
    p = 20201227
    g = 7

    b = dlog(g, door, p)
    powmod(card,b,p)
end

end


text:
md2 = fill("metadata", 2)
md3 = fill("metadata", 3)
md4 = fill("metadata", 4)
a = BagNode(ArrayNode(rand(3, 4)), [1:4], md4)
b = BagNode(ArrayNode(rand(3, 4)), [1:2, 3:4], md4)
c = BagNode(ArrayNode(rand(3, 4)), [1:1, 2:2, 3:4], md4)
d = BagNode(ArrayNode(rand(3, 4)), [1:4, 0:-1], md4)
wa = WeightedBagNode(ArrayNode(rand(3, 4)), [1:4], rand(1:4, 4), md4)
wb = WeightedBagNode(ArrayNode(rand(3, 4)), [1:2, 3:4], rand(1:4, 4), md4)
wc = WeightedBagNode(ArrayNode(rand(3, 4)), [1:1, 2:2, 3:4], rand(1:4, 4), md4)
wd = WeightedBagNode(ArrayNode(rand(3, 4)), [1:4, 0:-1], rand(1:4, 4), md4)
e = ArrayNode(rand(2, 2), md2)

f = ProductNode((wb,b), md2)
g = ProductNode([c, wc], md3)
h = ProductNode((wc,c), md3)
i = ProductNode((
              b,
              ProductNode((
                        b,
                        BagNode(
                                BagNode(
                                        ArrayNode(rand(2, 4)),
                                        [1:1, 2:2, 3:3, 4:4]
                                       ),
                                [1:3, 4:4]
                               )
                       ))
             ), md2)
k = ProductNode((a = wb, b = b), md2)
l = ProductNode((a = wc, b = c), md3)
m = ProductNode((a = wc, c = c), md3)
n = ProductNode((a = c, c = wc), md3)

@testset "testing nobs" begin
    @test nobs(a) == nobs(wa) == 1
    @test nobs(b) == nobs(wb) == 2
    @test nobs(c) == nobs(wc) == 3
    @test nobs(d) == nobs(wd) == 2
    @test nobs(e) == 2
    @test nobs(f) == nobs(wb) == nobs(b) == nobs(k)
    @test nobs(g) == nobs(c) == nobs(wc) == nobs(l)
    @test nobs(h) == nobs(wc) == nobs(c)
    @test nobs(i) == nobs(b)
end

@testset "testing ArrayNode hcat and vcat" begin
    @test catobs(e, e).data == hcat(e.data, e.data) == reduce(catobs, [e,e]).data
    @test hcat(e, e).data == hcat(e.data, e.data)
    @test vcat(e, e).data == vcat(e.data, e.data)
    x = ArrayNode(randn(2,3), rand(2, 3))
    @test catobs(x, x[0:-1]) isa ArrayNode{Matrix{Float64}, Matrix{Float64}}
    @inferred catobs(x, x[0:-1])
    @test reduce(catobs, [x, x[0:-1]]) isa ArrayNode{Matrix{Float64}, Matrix{Float64}}
    @inferred reduce(catobs, [x, x[0:-1]])
    @test cat(e, e, dims = ndims(e)).data == hcat(e.data, e.data)
end

@testset "testing BagNode catobs" begin
    @test catobs(a, b, c).data.data == hcat(a.data.data, b.data.data, c.data.data)
    @test reduce(catobs, [a, b, c]).data.data == hcat(a.data.data, b.data.data, c.data.data)
    @test catobs(a, b, c).bags.bags == vcat(a.bags, b.bags, c.bags).bags
    @test reduce(catobs, [a, b, c]).bags.bags == vcat(a.bags, b.bags, c.bags).bags
    @test catobs(c, a).data.data == hcat(c.data.data, a.data.data)
    @test reduce(catobs, [c, a]).data.data == hcat(c.data.data, a.data.data)
    @test catobs(c, a).bags.bags == vcat(c.bags, a.bags).bags
    @test reduce(catobs, [c, a]).bags.bags == vcat(c.bags, a.bags).bags
    @test catobs(a, c).data.data == hcat(a.data.data, c.data.data)
    @test reduce(catobs, [a, c]).data.data == hcat(a.data.data, c.data.data)
    @test catobs(a, c).bags.bags == vcat(a.bags, c.bags).bags
    @test reduce(catobs, [a, c]).bags.bags == vcat(a.bags, c.bags).bags
    @test catobs(a, d).data.data == hcat(a.data.data, d.data.data)
    @test reduce(catobs, [a, d]).data.data == hcat(a.data.data, d.data.data)
    @test catobs(a, d).bags.bags == vcat(a.bags, d.bags).bags
    @test reduce(catobs, [a, d]).bags.bags == vcat(a.bags, d.bags).bags
    @test catobs(d, a).data.data == hcat(d.data.data, a.data.data)
    @test reduce(catobs, [d, a]).data.data == hcat(d.data.data, a.data.data)
    @test catobs(d, a).bags.bags == vcat(d.bags, a.bags).bags
    @test reduce(catobs, [d, a]).bags.bags == vcat(d.bags, a.bags).bags
    @test catobs(d).data.data == hcat(d.data.data)
    @test reduce(catobs, [d]).data.data == hcat(d.data.data)
    @test catobs(d).bags.bags == vcat(d.bags).bags
    @test reduce(catobs, [d]).bags.bags == vcat(d.bags).bags
    @test cat(a, b, dims = ndims(a)).data.data == hcat(a.data.data, b.data.data)
end

@testset "testing WeightedBagNode catobs" begin
    @test catobs(wa, wb, wc).data.data == hcat(wa.data.data, wb.data.data, wc.data.data)
    @test reduce(catobs, [wa, wb, wc]).data.data == hcat(wa.data.data, wb.data.data, wc.data.data)
    @test catobs(wa, wb, wc).bags.bags == vcat(wa.bags, wb.bags, wc.bags).bags
    @test reduce(catobs, [wa, wb, wc]).bags.bags == vcat(wa.bags, wb.bags, wc.bags).bags
    @test catobs(wc, wa).data.data == hcat(wc.data.data, wa.data.data)
    @test reduce(catobs, [wc, wa]).data.data == hcat(wc.data.data, wa.data.data)
    @test catobs(wc, wa).bags.bags == vcat(wc.bags, wa.bags).bags
    @test reduce(catobs, [wc, wa]).bags.bags == vcat(wc.bags, wa.bags).bags
    @test catobs(wa, wc).data.data == hcat(wa.data.data, wc.data.data)
    @test reduce(catobs, [wa, wc]).data.data == hcat(wa.data.data, wc.data.data)
    @test catobs(wa, wc).bags.bags == vcat(wa.bags, wc.bags).bags
    @test reduce(catobs, [wa, wc]).bags.bags == vcat(wa.bags, wc.bags).bags
    @test catobs(wa, wd).data.data == hcat(wa.data.data, wd.data.data)
    @test reduce(catobs, [wa, wd]).data.data == hcat(wa.data.data, wd.data.data)
    @test catobs(wa, wd).bags.bags == vcat(wa.bags, wd.bags).bags
    @test reduce(catobs, [wa, wd]).bags.bags == vcat(wa.bags, wd.bags).bags
    @test catobs(wd, wa).data.data == hcat(wd.data.data, wa.data.data)
    @test reduce(catobs, [wd, wa]).data.data == hcat(wd.data.data, wa.data.data)
    @test catobs(wd, wa).bags.bags == vcat(wd.bags, wa.bags).bags
    @test reduce(catobs, [wd, wa]).bags.bags == vcat(wd.bags, wa.bags).bags
    @test catobs(wd).data.data == hcat(wd.data.data)
    @test reduce(catobs, [wd]).data.data == hcat(wd.data.data)
    @test catobs(wd).bags.bags == vcat(wd.bags).bags
    @test reduce(catobs, [wd]).bags.bags == vcat(wd.bags).bags
end

@testset "testing hierarchical catobs on ProductNodes" begin
    @test catobs(f, h).data[1].data.data == reduce(catobs, [f, h]).data[1].data.data ==
        hcat(wb.data.data, wc.data.data)
    @test catobs(f, h).data[2].data.data == reduce(catobs, [f, h]).data[2].data.data ==
        hcat(b.data.data, c.data.data)
    @test catobs(f, h, f).data[1].data.data == reduce(catobs, [f, h, f]).data[1].data.data ==
        hcat(wb.data.data, wc.data.data, wb.data.data)
    @test nobs(catobs(f,h)) == nobs(f) + nobs(h)

    @test catobs(g, g).data[1].data.data == reduce(catobs, [g, g]).data[1].data.data ==
        hcat(c.data.data, c.data.data)
    @test catobs(g, g).data[2].data.data == reduce(catobs, [g, g]).data[2].data.data ==
        hcat(wc.data.data, wc.data.data)
    @test catobs(g, g, g).data[1].data.data == reduce(catobs, [g, g, g]).data[1].data.data ==
        hcat(c.data.data, c.data.data, c.data.data)
    @test nobs(catobs(g, g)) == 2nobs(g)

    @test catobs(k, l).data[1].data.data == hcat(wb.data.data, wc.data.data)
    @test catobs(k, l).data[2].data.data == hcat(b.data.data, c.data.data)
    @test nobs(catobs(k,l)) == nobs(k) + nobs(l)

    # correct length/keyset but different subtrees
    @test_throws MethodError catobs(f, i)
    @test_throws MethodError reduce(catobs, [f, i])
    @test_throws MethodError catobs(m, n)
    @test_throws MethodError reduce(catobs, [m, n])
    # different tuple length or keyset
    @test_throws ArgumentError catobs(g, i)
    @test_throws ArgumentError reduce(catobs, [g, i])
    @test_throws ArgumentError catobs(f, g)
    @test_throws ArgumentError reduce(catobs, [f, g])
    @test_throws ArgumentError catobs(k, m)
    @test_throws ArgumentError reduce(catobs, [k, m])
    @test_throws ArgumentError catobs(l, m)
    @test_throws ArgumentError reduce(catobs, [l, m])
end

@testset "testing catobs with missing values" begin
    @test catobs(a, b, c).data.data == catobs(a, b, missing, c).data.data
    @test catobs(a, b, c).bags.bags == catobs(a, b, missing, c).bags.bags
    @test catobs(wa, wb, wc).data.data == catobs(wa, wb, missing, wc).data.data
    @test catobs(wa, wb, wc).bags.bags == catobs(wa, wb, missing, wc).bags.bags
    @test catobs(wa, wb, wc).weights == catobs(wa, wb, missing, wc).weights
end

@testset "testing catobs stability" begin
    for n in [a, b, c, d, e, f, h, k, l, wa, wb, wc, wd]
        @inferred catobs(n, n)
        @inferred reduce(catobs, [n, n])
    end
end

@testset "testing BagNode indexing" begin
    @test a[1].data.data == a.data.data
    @test a[1].bags.bags == [1:4]
    @test b[1:2].data.data == b.data.data
    @test b[1:2].bags.bags == [1:2,3:4]
    @test b[2].data.data == b.data.data[:,3:4]
    @test b[2].bags.bags == [1:2]
    @test b[1].data.data == b.data.data[:,1:2]
    @test b[1].bags.bags == [1:2]
    @test c[1].data.data == c.data.data[:,1:1]
    @test c[1].bags.bags == [1:1]
    @test c[[1,2]].data.data == c.data.data[:,1:2]
    @test c[[1,2]].bags.bags == [1:1,2:2]
    @test c[[2,1]].data.data == c.data.data[:,[2,1]]
    @test c[[2,1]].bags.bags == [1:1,2:2]
    @test d[[2,1]].data.data == d.data.data
    @test d[[2,1]].bags.bags == [0:-1,1:4]
    @test d[1:2].data.data == d.data.data
    @test d[1:2].bags.bags == [1:4,0:-1]
    @test nobs(d[2].data) == 0
    @test d[2].bags.bags == [0:-1]
    @test isempty(a[2:1].bags.bags)
    @test nobs(a[2:1].data) == 0
end

@testset "testing WeightedBagNode indexing" begin
    @test wa[1].data.data == wa.data.data
    @test wa[1].bags.bags == [1:4]
    @test wb[1:2].data.data == wb.data.data
    @test wb[1:2].bags.bags == [1:2,3:4]
    @test wb[2].data.data == wb.data.data[:,3:4]
    @test wb[2].bags.bags == [1:2]
    @test wb[1].data.data == wb.data.data[:,1:2]
    @test wb[1].bags.bags == [1:2]
    @test wc[1].data.data == wc.data.data[:,1:1]
    @test wc[1].bags.bags == [1:1]
    @test wc[[1,2]].data.data == wc.data.data[:,1:2]
    @test wc[[1,2]].bags.bags == [1:1,2:2]
    @test wc[[2,1]].data.data == wc.data.data[:,[2,1]]
    @test wc[[2,1]].bags.bags == [1:1,2:2]
    @test wd[[2,1]].data.data == wd.data.data
    @test wd[[2,1]].bags.bags == [0:-1,1:4]
    @test wd[1:2].data.data == wd.data.data
    @test wd[1:2].bags.bags == [1:4,0:-1]
    @test nobs(wd[2].data) == 0
    @test wd[2].bags.bags == [0:-1]
end

@testset "testing ProductNode indexing" begin
    @test h[1].data[1].data.data ==  wc[1].data.data
    @test h[[3,2]].data[1].data.data ==  wc[[3,2]].data.data
    @test h[1:2].data[2].data.data ==  c[1:2].data.data
    @test h[end].data[2].data.data ==  c[end].data.data
    @test g[1].data[1].data.data ==  c[1].data.data
    @test g[[3,2]].data[1].data.data ==  c[[3,2]].data.data
    @test g[1:2].data[2].data.data ==  wc[1:2].data.data
    @test g[end].data[2].data.data ==  wc[end].data.data
    @test l[1].data[:a].data.data ==  wc[1].data.data
    @test l[[3,2]].data[:a].data.data ==  wc[[3,2]].data.data
    @test l[1:2].data[:b].data.data ==  c[1:2].data.data
    @test l[end].data[:b].data.data ==  c[end].data.data
end

@testset "testing nested ragged array" begin
    x = BagNode(ArrayNode(rand(3,10)),[1:2,3:3,0:-1,4:5,6:6,7:10])
    y = BagNode(x,[1:2,3:3,4:5,6:6])
    @test y[1].data.data.data == x.data.data[:,1:3]
    @test y[1].data.bags.bags == [1:2,3:3]
    @test y[1:2].data.data.data == x.data.data[:,1:3]
    @test y[1:2].data.bags.bags == [1:2,3:3,0:-1]
    @test y[2:3].data.data.data == x.data.data[:,4:6]
    @test y[2:3].data.bags.bags == [0:-1,1:2,3:3]
end

@testset "testing sparsify" begin
    @test sparsify(zeros(10, 10), 0.05) isa SparseMatrixCSC
    @test sparsify(randn(10, 10), 0.05) isa Matrix
    @test sparsify(randn(10), 0.05) isa Vector
end

@testset "testing sparsify and mapdata" begin
    x = ProductNode((ProductNode((ArrayNode(randn(5,5)), ArrayNode(zeros(5,5)))), ArrayNode(zeros(5,5))))
    xs = mapdata(i -> sparsify(i, 0.05), x)
    @test xs.data[2].data isa SparseMatrixCSC
    @test xs.data[1].data[2].data isa SparseMatrixCSC
    @test xs.data[1].data[1].data == x.data[1].data[1].data
end

@testset "testing missing mapdata" begin
    x = ProductNode((ProductNode((ArrayNode(randn(5,5)), ArrayNode(zeros(5,5)))), ArrayNode(zeros(5,5))), BagNode(missing, AlignedBags([0:-1]), nothing))
    xs = mapdata(i -> sparsify(i, 0.05), x)
    @test xs.data[2].data isa SparseMatrixCSC
    @test xs.data[1].data[2].data isa SparseMatrixCSC
    @test xs.data[1].data[1].data == x.data[1].data[1].data
end

@testset "testing equals and hash" begin
    a2 = deepcopy(a)
    i2 = deepcopy(i)
    k2 = ProductNode((a = wb, b = b), md2)
    metadata1 = fill("metadata", 4)
    metadata2 = "Oh, Hi Mark"
    r = rand(3,4)
    x1 = BagNode(ArrayNode(r),[1:4], metadata1)
    x2 = BagNode(ArrayNode(r),[1:4], metadata1)
    x3 = BagNode(ArrayNode(r),[1:4], metadata2)

    @test a != b
    @test a != i
    @test a == a2
    @test i == i2
    @test i != k
    @test k == k2
    @test x1 == x2
    @test x1 != x3
    @test hash(a) !== hash(b)
    @test hash(a) === hash(a2)
    @test hash(a) !== hash(i)
    @test hash(i) === hash(i2)
    @test hash(i) !== hash(k)
    @test hash(k) === hash(k2)
    @test hash(x1) === hash(x2)
    @test hash(x1) !== hash(x3)
end

@testset "testing equals with missings" begin
    a = ArrayNode([0.0 missing 0.0 0.0 1.0])
    b = ArrayNode([0.0 missing 0.0 0.0 2.0])
    c = ArrayNode([0.0 missing 0.0 missing 2.0])
    @test a != a
    @test isequal(a, a)
    @test a != b
    @test !isequal(a, b)
    @test b != b
    @test isequal(b, b)
    @test a != c
    @test !isequal(a, c)
    @test b != c
    @test !isequal(b, c)
    @test c != c
    @test isequal(c, c)

    d = ProductNode((; a=a, b=b))
    e = ProductNode((; a=a, b=b), [])
    @test isequal(d, d)
    @test !isequal(d, e)
    @test d != d
    @test d != e

    f = BagNode(a,[1:4])
    g = BagNode(b,[1:4])
    @test f != f
    @test f != g
    @test isequal(f, f)
    @test !isequal(f, g)

    wf = WeightedBagNode(a, [1:4], rand(1:4, 4))
    wg = WeightedBagNode(b, [1:4], rand(1:4, 4))
    @test wf != wf
    @test wf != wg
    @test isequal(wf, wf)
    @test !isequal(wf, wg)

    h = LazyNode(:Test, a.data)
    i = LazyNode(:Test, b.data)
    @test isequal(h == h, missing)
    @test h != i
    @test isequal(h, h)
    @test !isequal(h, i)
end

@testset "dropmeta" begin
    for n in [a, b, c, d, e, f, h, k, l, wa, wb, wc, wd]
        @test dropmeta(n) |> Mill.metadata |> isnothing
    end
end


text:
#=
ndgrid.jl
lazy and non-lazy versions of ndgrid
=#

using LazyArrays: Applied, ApplyArray

export ndgrid_array_r, ndgrid_lazy

"""
    grid = ndgrid_repeat(v::AbstractVector, dims::Dims{D}, d::Int)
Do the type of repeat needed for `ndgrid`, e.g., `repeat(x, 1, n)`,
for arbitrary dimensions, where `1 ≤ d ≤ D` and `dims[d] == length(v)`.
Output size is `dims`.
"""
function ndgrid_repeat(v::AbstractVector, dims::Dims{D}, d::Int) where D
    @boundscheck checkbounds([dims...], d)
    @boundscheck length(v) == dims[d] ||
        throw(DimensionMismatch("$d $(dims[d]) $(length(v))"))
    t1 = ntuple(i -> i == d ? length(v) : 1, Val(D)) # (1,…,1,n,1,…,1)
    t2 = ntuple(i -> i == d ? 1 : dims[i], Val(D)) # (?,…,?,1,?,…,?)
    repeat(reshape(v, t1), t2...)
end

"""
    ndgrid_array_r(args::AbstractVector...)
Returns tuple of `length(args)` arrays, each of size `tuple(length.(args)...)`.
This method is provided for convenience and testing,
but `ndgrid` is less efficient than broadcast so should be avoided.
The tuple returned here requires `prod(length.(args)) * length(args)` memory;
using `ndgrid_lazy` is an alternative that uses `O(length(args))` memory.
"""
function ndgrid_array_r(args::AbstractVector...)
    D = length(args)
    T = Tuple{ntuple(i -> Array{eltype(args[i]), D}, D)...} # return type
    fun = i -> ndgrid_repeat(args[i], length.(args), i)
    ntuple(fun, D)::T
end

# lazy array method for a single element of the ndgrid tuple:

const _Repeat = Applied{A, typeof(ndgrid_repeat), Tuple{V,Dim,Int}} where
    {A <: Any, V <: AbstractVector, Dim <: Dims{D} where D}
Base.ndims(r::_Repeat) = length(r.args[2])
Base.size(r::_Repeat) = r.args[2]
Base.eltype(r::_Repeat) = eltype(r.args[1])

const Repeat{T,N} = ApplyArray{T, N, typeof(ndgrid_repeat)}
# Base.IndexStyle(Repeat) = IndexCartesian() # default

@inline Base.@propagate_inbounds function Base.getindex(
    r::Repeat{T,N},
    i::Vararg{Int,N},
) where {T,N}
    @boundscheck checkbounds(r, i...)
    return @inbounds r.args[1][i[r.args[3]]]
end

"""
    ndgrid_lazy(args::AbstractVector...)
Returns tuple of lazy array `Repeat` items,
avoiding the memory issues of `ndgrid`.

# Examples
```jldoctest
julia> ndgrid_lazy(1:3, 1:2)
([1 1; 2 2; 3 3], [1 2; 1 2; 1 2])

julia> ndgrid_lazy(1:3, 1:2)[1]
ndgrid_repeat(3-element UnitRange{Int64}, Tuple{Int64, Int64}, Int64):
 1  1
 2  2
 3  3
```
"""
function ndgrid_lazy(args::AbstractVector...)
    fun = i -> ApplyArray(ndgrid_repeat, args[i], length.(args), i)
    return ntuple(fun, Val(length(args)))
end