text:
SRC_DIR = joinpath(Base.source_dir(), "..", "src")
AUX_DIR = joinpath(SRC_DIR, "auxiliary")
OUT_DIR = joinpath(Base.source_dir(), "..", "output")
mkpath(OUT_DIR)

include(joinpath(SRC_DIR, "BridgeSDEInference.jl"))
using Main.BridgeSDEInference
using Distributions # to define priors
using Random        # to seed the random number generator
using DataFrames
using CSV
include(joinpath(AUX_DIR, "read_and_write_data.jl"))
include(joinpath(AUX_DIR, "transforms.jl"))

# decide if first passage time observations or partially observed diffusion
fptObsFlag = false

# pick dataset
filename = "path_part_obs_conj.csv"

# fetch the data
(df, x0, obs, obsTime, fpt,
      fptOrPartObs) = readData(Val(fptObsFlag), joinpath(OUT_DIR, filename))

param = :complexConjug
# Initial parameter guess.
θ₀ = [10.0, -8.0, 15.0, 0.0, 3.0]
# Target law
P˟ = FitzhughDiffusion(param, θ₀...)
# Auxiliary law
P̃ = [FitzhughDiffusionAux(param, θ₀..., t₀, u[1], T, v[1]) for (t₀,T,u,v)
     in zip(obsTime[1:end-1], obsTime[2:end], obs[1:end-1], obs[2:end])]
display(P̃[1])

L = @SMatrix [1. 0.]
Σdiagel = 10^(-10)
Σ = @SMatrix [Σdiagel]

Ls = [L for _ in P̃]
Σs = [Σ for _ in P̃]
τ(t₀,T) = (x) ->  t₀ + (x-t₀) * (2-(x-t₀)/(T-t₀))
numSteps=1*10^5
saveIter=3*10^2
tKernel = RandomWalk([3.0, 5.0, 5.0, 0.01, 0.5],
                     [false, false, false, false, true])
priors = Priors((MvNormal([0.0,0.0,0.0], diagm(0=>[1000.0, 1000.0, 1000.0])),
                 #ImproperPrior(),
                 ImproperPrior()))
𝔅 = NoBlocking()
blockingParams = ([], 0.1, NoChangePt())
changePt = NoChangePt()
#x0Pr = KnownStartingPt(x0)
x0Pr = GsnStartingPt(x0, x0, @SMatrix [20. 0; 0 20.])
warmUp = 100

Random.seed!(4)
start = time()
(chain, accRateImp, accRateUpdt,
    paths, time_) = mcmc(eltype(x0), fptOrPartObs, obs, obsTime, x0Pr, 0.0, P˟,
                         P̃, Ls, Σs, numSteps, tKernel, priors, τ;
                         fpt=fpt,
                         ρ=0.975,
                         dt=1/1000,
                         saveIter=saveIter,
                         verbIter=10^2,
                         updtCoord=(Val((true, true, true, false, false)),
                                    #Val((true, false, false, false, false)),
                                    Val((false, false, false, false, true)),
                                    ),
                         paramUpdt=true,
                         updtType=(ConjugateUpdt(),
                                   #MetropolisHastingsUpdt(),
                                   MetropolisHastingsUpdt(),
                                   ),
                         skipForSave=10^0,
                         blocking=𝔅,
                         blockingParams=blockingParams,
                         solver=Vern7(),
                         changePt=changePt,
                         warmUp=warmUp)
elapsed = time() - start
print("time elapsed: ", elapsed, "\n")

print("imputation acceptance rate: ", accRateImp,
      ", parameter update acceptance rate: ", accRateUpdt)

x0⁺, pathsToSave = transformMCMCOutput(x0, paths, saveIter; chain=chain,
                                       numGibbsSteps=2,
                                       parametrisation=param,
                                       warmUp=warmUp)

df2 = savePathsToFile(pathsToSave, time_, joinpath(OUT_DIR, "sampled_paths.csv"))
df3 = saveChainToFile(chain, joinpath(OUT_DIR, "chain.csv"))

include(joinpath(AUX_DIR, "plotting_fns.jl"))
set_default_plot_size(30cm, 20cm)
plotPaths(df2, obs=[Float64.(df.x1), [x0⁺[2]]],
          obsTime=[Float64.(df.time), [0.0]], obsCoords=[1,2])

plotChain(df3, coords=[1])
plotChain(df3, coords=[2])
plotChain(df3, coords=[3])
plotChain(df3, coords=[5])


text:
function powerlaw(x, k, G)

    # set negatives to zero -
    x = max.(0.0, x)

    # get system dimension -
    ℛ = length(k)
    ℳ = length(x)

    # initialize the rate array -
    rate_array = Array{Float64,1}(undef,ℛ)

    for i ∈ 1:ℛ

        # tmp -
        tmp_value = 1.0
        for j ∈ 1:ℳ
            tmp_value = tmp_value*(x[j])^(G[j,i])
        end

        # store the rate value -
        rate_array[i] = k[i]*tmp_value
    end

    # return -
    return rate_array
end

ext: jl

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
struct TestModel <: WhittleLikelihoodInference.TimeSeriesModel{2,Float64} end
struct TestModelUni <: WhittleLikelihoodInference.TimeSeriesModel{1,Float64} end
struct TestModel2 <: WhittleLikelihoodInference.TimeSeriesModel{2,Float64} end
struct TestModel2C <: WhittleLikelihoodInference.TimeSeriesModel{2,ComplexF64} end
struct TestModelUni2 <: WhittleLikelihoodInference.TimeSeriesModel{1,Float64} end
struct TestModelUni2C <: WhittleLikelihoodInference.TimeSeriesModel{1,ComplexF64} end
struct TestModel3 <: WhittleLikelihoodInference.TimeSeriesModel{2,Float64} end
struct TestModelUni3 <: WhittleLikelihoodInference.TimeSeriesModel{1,Float64} end
TestModel2(x) = TestModel2() # just to allow additive constructor to work
TestModelUni2(x) = TestModelUni2()
TestModel3(x) = TestModel3() # just to allow additive constructor to work
TestModelUni3(x) = TestModelUni3()
WhittleLikelihoodInference.npars(::Union{Type{TestModel2},Type{TestModel2C}}) = 4
WhittleLikelihoodInference.npars(::Union{Type{TestModelUni2},Type{TestModelUni2C}}) = 4
WhittleLikelihoodInference.npars(::Type{TestModel3}) = 4
WhittleLikelihoodInference.npars(::Type{TestModelUni3}) = 4
const TestDouble = TestModel2 + TestModel3
const TestDoubleUni = TestModelUni2 + TestModelUni3
WhittleLikelihoodInference.add_sdf!(out,::Union{TestModel2,TestModel2C},ω) = (out[2] += ω>0 ? 0.5+0.5im : 0.5-0.5im; nothing)
WhittleLikelihoodInference.grad_add_sdf!(out,::Union{TestModel2,TestModel2C},ω) = (out[2,:] .+= ω>0 ? 0.5+0.5im : 0.5-0.5im; nothing)
WhittleLikelihoodInference.hess_add_sdf!(out,::Union{TestModel2,TestModel2C},ω) = (out[2,:] .+= ω>0 ? 0.5+0.5im : 0.5-0.5im; nothing)
WhittleLikelihoodInference.sdf(::Union{TestModelUni2,TestModelUni2C},ω) = 0.5
WhittleLikelihoodInference.grad_add_sdf!(out,::Union{TestModelUni2,TestModelUni2C},ω) = (out .+= 0.5; nothing)
WhittleLikelihoodInference.hess_add_sdf!(out,::Union{TestModelUni2,TestModelUni2C},ω) = (out .+= 0.5; nothing)
WhittleLikelihoodInference.add_sdf!(out,::TestModel3,ω) = (out .+= 0.2; nothing)
WhittleLikelihoodInference.grad_add_sdf!(out,::TestModel3,ω) = (out .+= 0.2; nothing)
WhittleLikelihoodInference.hess_add_sdf!(out,::TestModel3,ω) = (out .+= 0.2; nothing)
WhittleLikelihoodInference.sdf(::TestModelUni3,ω) = 0.2
WhittleLikelihoodInference.grad_add_sdf!(out,::TestModelUni3,ω) = (out .+= 0.2; nothing)
WhittleLikelihoodInference.hess_add_sdf!(out,::TestModelUni3,ω) = (out .+= 0.2; nothing)
WhittleLikelihoodInference.nalias(::Union{TestModel2,TestModel2C}) = 2
WhittleLikelihoodInference.nalias(::Union{TestModelUni2,TestModelUni2C}) = 2
WhittleLikelihoodInference.nalias(::TestModel3) = 1
WhittleLikelihoodInference.nalias(::TestModelUni3) = 1
@testset "Spectral density" begin
    @testset "Function" begin
        @testset "Error handling" begin
            @test_throws ErrorException WhittleLikelihoodInference.add_sdf!(ones(2),TestModel(),1.0)
            @test_throws ErrorException WhittleLikelihoodInference.sdf(TestModelUni(),1.0)
        end
        @testset "sdf" begin
            @test sdf(TestModel2(),1.0) == SHermitianCompact(@SVector [0, complex(0.5,0.5), 0])
            @test sdf(TestModelUni2(),1.0) == 0.5
            @test sdf(TestDouble(zeros(8)),1.0) == SHermitianCompact(@SVector [0.2, complex(0.7,0.5), 0.2])
            @test sdf(TestDoubleUni(zeros(8)),1.0) == 0.7
        end
        @testset "asdf" begin
            @test asdf(TestModel2(),1.0,1.0) == SHermitianCompact(@SVector [0, complex(2.5,0.5), 0])
            @test asdf(TestModelUni2(),1.0,1.0) == 2.5 # two lots of aliasing
            @test asdf(TestDouble(zeros(8)),1.0,1.0) ≈ SHermitianCompact(@SVector [0.6, complex(3.1,0.5), 0.6])
            @test asdf(TestDoubleUni(zeros(8)),1.0,1.0) == 3.1 # two lots of aliasing of 0.5 plus one lot of 0.2 (but has floating point error)
        end
        @testset "Bounds check" begin
            store = allocate_memory_sdf_F(TestModel2, 1000, 1)
            storeuni = allocate_memory_sdf_F(TestModelUni2, 1000, 1)
            freqbad = FreqAcvEst(10,1)
            @test_throws ArgumentError WhittleLikelihoodInference.asdf!(store.funcmemory,TestModel2(),freqbad)
            @test_throws ArgumentError WhittleLikelihoodInference.asdf!(storeuni.funcmemory,TestModelUni2(),freqbad)
        end
        @testset "Complex vs real" begin
            for n in [100, 101]
                store = allocate_memory_sdf_F(TestModel2, n, 1.2)
                store_c = allocate_memory_sdf_F(TestModel2C, n, 1.2)
                storeuni = allocate_memory_sdf_F(TestModelUni2, n, 1.2)
                storeuni_c = allocate_memory_sdf_F(TestModelUni2C, n, 1.2)
                WhittleLikelihoodInference.asdf!(store,TestModel2())
                WhittleLikelihoodInference.asdf!(store_c,TestModel2C())
                WhittleLikelihoodInference.asdf!(storeuni,TestModelUni2())
                WhittleLikelihoodInference.asdf!(storeuni_c,TestModelUni2C())
                @test store.funcmemory.allocatedarray == store_c.funcmemory.allocatedarray
                @test storeuni.funcmemory.allocatedarray == storeuni_c.funcmemory.allocatedarray
            end
        end
        @testset "extract_asdf" begin
            @test_throws MethodError extract_asdf(allocate_memory_EI_F(TestModel2, 1000, 1).funcmemory)
            @test extract_asdf(allocate_memory_sdf_F(TestModel2, 1000, 1).funcmemory) isa Matrix{ComplexF64}
        end
    end
    @testset "Gradient" begin
        @testset "Error handling" begin
            @test_throws ErrorException WhittleLikelihoodInference.add_sdf!(ones(2),TestModel(),1.0)
            @test_throws ErrorException WhittleLikelihoodInference.sdf(TestModelUni(),1.0)
        end
        @testset "sdf" begin
            @test grad_sdf(TestModel2(),1.0) == repeat([0.0,complex(0.5,0.5),0.0],inner=(1,4))
            @test grad_sdf(TestModelUni2(),1.0) == fill(complex(0.5),4)
            @test grad_sdf(TestDouble(zeros(8)),1.0) == hcat(repeat([0.0,complex(0.5,0.5),0.0],inner=(1,4)),fill(complex(0.2),3,4))
            @test grad_sdf(TestDoubleUni(zeros(8)),1.0) == [fill(complex(0.5),4);fill(complex(0.2),4)]
        end
        @testset "Bounds check" begin
            store = allocate_memory_sdf_FG(TestModel2, 1000, 1)
            storeuni = allocate_memory_sdf_FG(TestModelUni2, 1000, 1)
            freqbad = FreqAcvEst(10,1)
            @test_throws ArgumentError WhittleLikelihoodInference.grad_asdf!(store.gradmemory,TestModel2(),freqbad)
            @test_throws ArgumentError WhittleLikelihoodInference.grad_asdf!(storeuni.gradmemory,TestModelUni2(),freqbad)
        end
        @testset "Complex vs real" begin
            for n in [100, 101]
                store = allocate_memory_sdf_FG(TestModel2, n, 1.2)
                store_c = allocate_memory_sdf_FG(TestModel2C, n, 1.2)
                storeuni = allocate_memory_sdf_FG(TestModelUni2, n, 1.2)
                storeuni_c = allocate_memory_sdf_FG(TestModelUni2C, n, 1.2)
                WhittleLikelihoodInference.grad_asdf!(store,TestModel2())
                WhittleLikelihoodInference.grad_asdf!(store_c,TestModel2C())
                WhittleLikelihoodInference.grad_asdf!(storeuni,TestModelUni2())
                WhittleLikelihoodInference.grad_asdf!(storeuni_c,TestModelUni2C())
                @test store.gradmemory.allocatedarray == store_c.gradmemory.allocatedarray
                @test storeuni.gradmemory.allocatedarray == storeuni_c.gradmemory.allocatedarray
            end
        end
    end
    @testset "Hessian" begin
        @testset "Error handling" begin
            @test_throws ErrorException WhittleLikelihoodInference.add_sdf!(ones(2),TestModel(),1.0)
            @test_throws ErrorException WhittleLikelihoodInference.sdf(TestModelUni(),1.0)
        end
        @testset "sdf" begin
            @test hess_sdf(TestModel2(),1.0) == repeat([0.0,complex(0.5,0.5),0.0],inner=(1,10))
            @test hess_sdf(TestModelUni2(),1.0) == fill(complex(0.5),10)
            @test hess_sdf(TestDouble(zeros(8)),1.0) == hcat(repeat([0.0,complex(0.5,0.5),0.0],inner=(1,10)),fill(complex(0.2),3,10))
            @test hess_sdf(TestDoubleUni(zeros(8)),1.0) == [fill(complex(0.5),10);fill(complex(0.2),10)]
        end
        @testset "Bounds check" begin
            store = allocate_memory_sdf_FGH(TestModel2, 1000, 1)
            storeuni = allocate_memory_sdf_FGH(TestModelUni2, 1000, 1)
            freqbad = FreqAcvEst(10,1)
            @test_throws ArgumentError WhittleLikelihoodInference.hess_asdf!(store.hessmemory,TestModel2(),freqbad)
            @test_throws ArgumentError WhittleLikelihoodInference.hess_asdf!(storeuni.hessmemory,TestModelUni2(),freqbad)
        end
        @testset "Complex vs real" begin
            for n in [100, 101]
                store = allocate_memory_sdf_FGH(TestModel2, n, 1.2)
                store_c = allocate_memory_sdf_FGH(TestModel2C, n, 1.2)
                storeuni = allocate_memory_sdf_FGH(TestModelUni2, n, 1.2)
                storeuni_c = allocate_memory_sdf_FGH(TestModelUni2C, n, 1.2)
                WhittleLikelihoodInference.hess_asdf!(store,TestModel2())
                WhittleLikelihoodInference.hess_asdf!(store_c,TestModel2C())
                WhittleLikelihoodInference.hess_asdf!(storeuni,TestModelUni2())
                WhittleLikelihoodInference.hess_asdf!(storeuni_c,TestModelUni2C())
                @test store.hessmemory.allocatedarray == store_c.hessmemory.allocatedarray
                @test storeuni.hessmemory.allocatedarray == storeuni_c.hessmemory.allocatedarray
            end
        end
    end
end


text:
module collective_modes

export Omega_k_chain, Gamma_k_chain, Omega_k_2D, Gamma_k_2D

# For a particularly elegant derivation of collective frequency shifts and decay rates in 1D
# atomic chains, see Asenjo-Garcia et al 10.1103/PhysRevX.7.031024.

import LinearAlgebra
import Polylogarithms

const la = LinearAlgebra
const pl = Polylogarithms
const k0 = 2pi
const Rot90 = [0 -1; 1 0]

"""
    polarization_renorm(polarization)

INTERNAL FUNCTION.
Normalizes polarization vector and separates in and out-of-lattice-plane components.

# Arguments
* `polarization`: 3D, complex vector of atomic polarization.
"""
function polarization_renorm(polarization)
    polarization = la.normalize(polarization)
    x_comp = abs(polarization[1])^2
    yz_comp = 1 - x_comp
    polarization = polarization[2:3]
    return polarization, x_comp, yz_comp
end

"""
    N_lightcone(rec, rec_mag)

INTERNAL FUNCTION.
Number of reciprocal lattice vectors that are within the lightcone.

# Arguments
* `rec`: Magnitude of reciprocal lattice vector overlap along axis.
* `rec_mag`: Magnitude of reciprocal lattice vector.
"""
function N_lightcone(rec, rec_mag)
    return floor(abs(k0 - rec)/rec_mag)
end

"""
    Omega_k_chain(k, a, polarization)

Collective frequency shift Omega_k of mode k for an infinite chain of atoms along x-axis.

WLOG, this calculation scales natural atomic frequency wavelength lambda0=1 and decay rate gamma0=1.

# Arguments
* `k`: x-axis quasimomentum of collective mode in first BZ such that |k|<= pi/a.
* `a`: Atomic lattice spacing.
* `polarization`: 3D, complex vector of atomic polarization.
"""
function Omega_k_chain(k::Real, a::Real, polarization::Array{<:Number, 1})
    polarization, x_comp, yz_comp = polarization_renorm(polarization)
    exp_sum = exp(im*(k0+k)*a)
    exp_diff = exp(im*(k0-k)*a)
    diag = pl.polylog(3, exp_sum) + pl.polylog(3, exp_diff) - im*k0*a*(pl.polylog(2, exp_sum) + pl.polylog(2, exp_diff))
    perp = (k0*a)^2*(log(1-exp_sum) + log(1-exp_diff))
    return 3/(4*(k0*a)^3) * (-2*x_comp*real(diag) + yz_comp*real(diag+perp))
end

"""
    Gamma_k_chain(k, a, polarization)

Collective decay rate Gamma_k of mode k for an infinite chain of atoms along x-axis.

WLOG, this calculation scales natural atomic frequency wavelength lambda0=1 and decay rate gamma0=1.

# Arguments
* `k`: x-axis quasimomentum of collective mode in first BZ such that |k|<= pi/a.
* `a`: Atomic lattice spacing.
* `polarization`: 3D, complex vector of atomic polarization.
"""
function Gamma_k_chain(k::Real, a::Real, polarization::Array{<:Number, 1})
    if abs(k) > k0
        return 0
    end
    polarization, x_comp, yz_comp = polarization_renorm(polarization)
    g = 2pi/a
    N = N_lightcone(k, g)
    Gamma_sum = 0
    for n = -N:N
        Gamma_sum += (k + n*g)^2
    end
    Gamma_sum = Gamma_sum/k0^2
    return 3pi/(4*k0*a) * (2*x_comp*(2*N+1-Gamma_sum) + yz_comp*(2*N+1+Gamma_sum))
end

# For a particularly elegant derivation of collective frequency shifts and decay rates in 2D
# atomic arrays, see the thesis of Dominik Wild, "Algorithms and Platforms for Quantum Science and Technology".

"""
    Bravais_consts(a_vec1, a_vec2)

INTERNAL FUNCTION.
Returns the BZ volume and reciprocal lattice vectors.

# Arguments
* `a_vec1, a_vec2`: 2D Bravais lattice vectors.
"""
function Bravais_consts(a_vec1::Array{<:Number, 1}, a_vec2::Array{<:Number, 1})
    V = abs(a_vec1[1]*a_vec2[2] - a_vec1[2]*a_vec2[1])
    b1 = 2pi * Rot90*a_vec2 / la.dot(a_vec1, Rot90*a_vec2)
    b2 = 2pi * Rot90*a_vec1 / la.dot(a_vec2, Rot90*a_vec1)
    return V, b1, b2
end

"""
    renorm_consts(V, b1, b2, Lambda, N1, N2)

INTERNAL FUNCTION.
Returns default cutoff frequency and number of reciprocal lattice terms to be summed if
keyword arguments not specified.

# Arguments
* `V`: Volume of BZ.
* `b1, b2`: Reciprocal lattice vectors from a_vec1 and a_vec2.
* `Lambda`: Cutoff frequency of renormalization.
* `N1`: Number of terms in a_vec1 reciprocal lattice direction.
* `N2`: Number of terms in a_vec2 reciprocal lattice direction.
"""
function renorm_consts(V, b1, b2, Lambda, N1, N2)
    if Lambda == nothing
        Lambda = 10/sqrt(V)
    end
    if N1 == nothing
        N1 = Int(ceil(5*Lambda/la.norm(b1)))
    end
    if N2 == nothing
        N2 = Int(ceil(5*Lambda/la.norm(b2)))
    end
    return Lambda, N1, N2
end

"""
    Omega_k_2D(k_vec, a_vec1, a_vec2, polarization)

Collective frequency shift Omega_k of in-plane mode k_vec for an 2D array of atoms in yz-plane.

WLOG, this calculation scales natural atomic frequency wavelength lambda0=1 and decay rate gamma0=1.

# Arguments
* `k_vec`: yz-axis quasimomentum of collective mode in first BZ such that |k|<= pi/a.
* `a_vec1, a_vec2`: 2D Bravais lattice vectors.
* `polarization`: 3D, complex vector of atomic polarization.
* `Lambda`: Cutoff frequency of renormalization.
* `N1`: Number of terms in a_vec1 reciprocal lattice direction.
* `N2`: Number of terms in a_vec2 reciprocal lattice direction.
"""
function Omega_k_2D(k_vec::Array{<:Number, 1}, a_vec1::Array{<:Number, 1}, a_vec2::Array{<:Number, 1}, polarization::Array{<:Number, 1}; Lambda=nothing, N1=nothing, N2=nothing)
    V, b1, b2 = Bravais_consts(a_vec1, a_vec2)
    Lambda, N1, N2 = renorm_consts(V, b1, b2, Lambda, N1, N2)
    polarization, x_comp, yz_comp = polarization_renorm(polarization)

    Omega_sum_x = 0
    Omega_sum_yz = 0
    for n1 = -N1:N1
        for n2 =-N2:N2
            g_eff = k_vec + n1*b1 + n2*b2
            discrim = sqrt(Complex(k0^2 - la.norm(g_eff)^2))
            if abs(discrim) > 0
                cut_discrim = im*exp(-(la.norm(g_eff)/Lambda)^2) / discrim
            else
                cut_discrim = 0
            end
            Omega_sum_x += real(la.norm(g_eff)^2 * cut_discrim)
            Omega_sum_yz += real((yz_comp*k0^2 - abs(la.dot(g_eff, polarization))^2) * cut_discrim)
        end
    end
    Omega0x = ( 2*k0^2*V*k0 / (32*sqrt(pi)) ) * ( 4*(Lambda/k0) + 2*(Lambda/k0)^3 ) * exp(-(k0/Lambda)^2)
    Omega0yz = ( 2*k0^2*V*k0 / (32*sqrt(pi)) ) * ( 2*(Lambda/k0) - (Lambda/k0)^3 ) * exp(-(k0/Lambda)^2)
    return -3pi/(2*k0^3*V) * (x_comp*(Omega_sum_x-Omega0x) + Omega_sum_yz - yz_comp*Omega0yz)
end

"""
    Gamma_k_2D(k_vec, a_vec1, a_vec2, polarization)

Collective decay rate Gamma_k of in-plane mode k_vec for an 2D array of atoms in yz-plane.

WLOG, this calculation scales natural atomic frequency wavelength lambda0=1 and decay rate gamma0=1.

# Arguments
* `k_vec`: yz-axis quasimomentum of collective mode in first BZ such that |k|<= pi/a.
* `a_vec1, a_vec2`: 2D Bravais lattice vectors.
* `polarization`: 3D, complex vector of atomic polarization.
"""
function Gamma_k_2D(k_vec::Array{<:Number, 1}, a_vec1::Array{<:Number, 1}, a_vec2::Array{<:Number, 1}, polarization::Array{<:Number, 1})
    if la.norm(k_vec) > k0
        return 0
    end
    V, b1, b2 = Bravais_consts(a_vec1, a_vec2)
    polarization, x_comp, yz_comp = polarization_renorm(polarization)

    N1 = N_lightcone(la.norm(la.dot(k_vec, b1)), la.norm(b1))
    N2 = N_lightcone(la.norm(la.dot(k_vec, b2)), la.norm(b2))
    Gamma_sum_x = 0
    Gamma_sum_yz = 0
    for n1 = -N1:N1
        for n2 =-N2:N2
            g_eff = k_vec + n1*b1 + n2*b2
            if la.norm(g_eff) < k0
                discrim = sqrt(k0^2 - la.norm(g_eff)^2)
                Gamma_sum_x += la.norm(g_eff)^2 / discrim
                Gamma_sum_yz += (yz_comp*k0^2 - abs(la.dot(g_eff, polarization))^2) / discrim
            end
        end
    end
    return 3pi/(k0^3*V) * (x_comp*Gamma_sum_x + Gamma_sum_yz)
end

end # module


text:
using BlockArrays, Test

@testset "foldl" begin
    x = mortar([rand(3), rand(2)])
    @test foldl(push!, x; init = []) == collect(x)

    x = PseudoBlockVector(rand(3), [1, 2])
    @test foldl(push!, x; init = []) == collect(x)
end

@testset "reduce" begin
    x = mortar([rand(Int, 3), rand(Int, 2)])
    @test reduce(+, x) == sum(collect(x))

    x = PseudoBlockVector(rand(Int, 3), [1, 2])
    @test reduce(+, x) == sum(collect(x))
end

@testset "sum (#141)" begin
    data = reshape(collect(1:20), 4, 5)
    A = BlockArray(data, [1,3], [2,3])
    @test sum(A) == sum(data)
    @test sum(A; dims=1) == sum(data; dims=1)
    @test sum(A; dims=2) == sum(data; dims=2)
    @test blockisequal(axes(A,2), axes(sum(A; dims=1),2))
    @test blockisequal(axes(A,1), axes(sum(A; dims=2),1))
end


text:
using Pkg.Artifacts
using polymake_jll

# adapted from https://github.com/JuliaLang/julia/pull/38797
function full_artifact_dir(m::Module)
   artifacts_toml = joinpath(dirname(dirname(Base.pathof(m))), "StdlibArtifacts.toml")

   # If this file exists, it's a stdlib JLL and we must download the artifact ourselves
   if isfile(artifacts_toml)
       # we need to remove the _jll for the artifact name
       meta = artifact_meta(string(m)[1:end-4], artifacts_toml)
       hash = Base.SHA1(meta["git-tree-sha1"])
       if !artifact_exists(hash)
           dl_info = first(meta["download"])
           download_artifact(hash, dl_info["url"], dl_info["sha256"])
       end
       return artifact_path(hash)
   else
      # Otherwise, we can just use the artifact directory given to us by the module
      return m.artifact_dir
   end
end

function prepare_deps_tree(targetdir::String)
   # create a directory tree for polymake with links to dependencies
   # looking similiar to the tree in the build environment
   # for compiling wrappers at run-time
   mkpath(joinpath(targetdir,"deps"))
   deps = [ polymake_jll.FLINT_jll,
            polymake_jll.GMP_jll,
            polymake_jll.MPFR_jll,
            polymake_jll.PPL_jll,
            polymake_jll.Perl_jll,
            polymake_jll.bliss_jll,
            polymake_jll.boost_jll,
            polymake_jll.cddlib_jll,
            polymake_jll.lrslib_jll,
            polymake_jll.normaliz_jll ]

   # dependencies
   for dep in deps
      target = joinpath(targetdir,"deps","$dep")
      rm(target, force=true)
      symlink(full_artifact_dir(dep), target)
   end

   # polymake prefix directories
   for dir in filter(d -> d!="bin", readdir(polymake_jll.artifact_dir))
      target = joinpath(targetdir,dir)
      rm(target, force=true)
      symlink(joinpath(polymake_jll.artifact_dir,dir), target)
   end

   # polymake perl scripts
   bindir(name) = joinpath(targetdir,"bin", name)
   mkpath(bindir(""))
   for file in ("polymake", "polymake-config")
      rm(bindir(file), force=true)
      symlink(joinpath(polymake_jll.artifact_dir, "bin", file), bindir(file))
   end

   # Point polymake to our custom tree
   return targetdir
end

max_stars_repo_name: cometscome/Yggdrasil

max_issues_repo_name: cometscome/Yggdrasil

max_forks_repo_name: cometscome/Yggdrasil


text:
"""
    _get_lp_relaxation_vertices(univariate_function_data::UnivariateFunctionData)::Vector{Vertex2d}

Return vertices of the LP relaxation of the given function.
"""
function _get_lp_relaxation_vertices(
    univariate_function_data::UnivariateFunctionData,
)::Vector{Vertex2d}
    sec_vs, tan_vs = _collect_vertices(univariate_function_data)
    vertices = Vertex2d[]
    push!(vertices, sec_vs[1])
    append!(vertices, tan_vs)
    push!(vertices, sec_vs[end])
    return vertices
end

"""
    _build_univariate_lp_relaxation!(m, x, y, univariate_function_data, pre_base_name)

Build LP relaxation for ``y=f(x)`` given the univariate function data.
"""
function _build_univariate_lp_relaxation!(
    m::JuMP.Model,
    x::JuMP.VariableRef,
    y::JuMP.VariableRef,
    univariate_function_data::UnivariateFunctionData,
    pre_base_name::AbstractString,
)::FormulationInfo
    vertices = _get_lp_relaxation_vertices(univariate_function_data)
    num_vars = length(vertices)
    formulation_info = FormulationInfo()

    # add variables
    lambda =
        formulation_info.variables[:lambda] = @variable(
            m,
            [1:num_vars],
            lower_bound = 0.0,
            upper_bound = 1.0,
            base_name = pre_base_name * "_lambda"
        )
    formulation_info.variables[:lambda] = lambda

    # add constraints
    formulation_info.constraints[:sum_lambda] = @constraint(m, sum(lambda) == 1)
    formulation_info.constraints[:x] =
        @constraint(m, x == sum(lambda[i] * vertices[i][1] for i in 1:num_vars))
    formulation_info.constraints[:y] =
        @constraint(m, y == sum(lambda[i] * vertices[i][2] for i in 1:num_vars))

    return formulation_info
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 16ba1b52-6ca5-4cdb-8657-5c0ec912bb56
using SparkSQL, DataFrames

# ╔═╡ 1a82a4df-fd0c-4c5f-b240-9672b20d017d
using GLM

# ╔═╡ 010e59a4-0933-11ec-005b-3314aac96bb8
# SparkSQL.jl machine Learning tutorial.

# ╔═╡ f44c2335-bed9-4bf2-ba4e-a30a7e26fe3e
# Example shows how to build a logistic regression machine learning model with  Julia and SparkSQL.jl. A sample advertising use case is used to provide a real world example. The example model predicts which ad campaign and media format combination produces a "click".

# ╔═╡ b2356d86-4c41-4e13-b1e6-c9fc64f9993c
SparkSQL.initJVM()

# ╔═╡ f68d5f13-89c6-4755-b8ec-c123b87c2c35
sprk = SparkSQL.SparkSession("spark://localhost:7077", "Julia & Spark Machine Learning Example")

# ╔═╡ 8f224f26-af8b-45db-b007-b542189ae20f
data = sql(sprk, "SELECT CAST(_c0 AS INTEGER) AS click, _c1 AS media, _c2 as ad_shown FROM CSV.`sampleDataML.csv`");

# ╔═╡ 196603e4-270a-45b8-afb9-4649599e78f2
# Create a temporary database view in Spark with Julia to allow queries on loaded data.

# ╔═╡ 1392d02d-0258-4183-b756-969a14dfd033
createOrReplaceTempView(data,"advertisement_data");

# ╔═╡ 35950c0c-f87a-425f-ae96-f7555dd8eaba
# Create Training data set. Get a random sample of the data for training.

# ╔═╡ 9d228706-000d-41c3-8162-9e8894eed559
trainingData = sql(sprk, "SELECT * FROM advertisement_data TABLESAMPLE(40 PERCENT)")

# ╔═╡ b22ed63c-b646-4e43-88bf-d1e5c754333a
trainingDF = toJuliaDF(trainingData)

# ╔═╡ b5ed1dd2-d1af-4d22-a6da-8b5586ccac9d
# Create test data set to validate the model. Get a random sample of the data for testing the model.

# ╔═╡ 41acbbe9-2c68-493d-b7e9-9811817a4f5a
testData = sql(sprk,"SELECT * FROM advertisement_data TABLESAMPLE(25 PERCENT)")

# ╔═╡ 4ac78ed0-ab59-4a05-b2d3-2e5143e781a9
testDF = toJuliaDF(testData)

# ╔═╡ a5c2fe1d-aa59-4048-91ad-ed00262e0f57
# The GLM (Generalized Linear Model) Julia package.

# ╔═╡ 75eafc32-fc5b-44aa-9d06-cf85c8b1d338
# Create the Logistic regression machine learning model. Ad was clicked based on the media and ad shown.

# ╔═╡ a8dfa116-2677-4df1-acc6-738b27a864bf
f = @formula(click ~ media + ad_shown)

# ╔═╡ 0254392c-61b0-4198-860f-f7450f1e04e1
model = glm(f, trainingDF, Binomial(), LogitLink());

# ╔═╡ 6b5efee9-ceb5-488c-ae8f-5adf107bf7f5
coeftable(model);

# ╔═╡ 1ec2a9c1-0033-4320-84eb-74d665bcec51
# Run the model against the test data.

# ╔═╡ 4a4fe208-f98a-4952-b0b8-f32e8a63be61
prediction = round.(predict(model, testDF), digits=3);

# ╔═╡ 15265b71-0392-411d-a9b9-77b4e72a638e
# Compare model's performance against test data.

# ╔═╡ 31a75e65-c4ab-43e0-8566-954b7274add2
insertcols!(testDF, 1, :model_predicted => [if p >= 0.5 1 else 0 end for p in prediction])

# ╔═╡ cbdc3eef-cd65-4101-bf7e-b649f632939a
# The model shows a "1" if predicting the ad was clicked and a "0" if not.

# ╔═╡ 69ab1393-ae18-4761-a188-14f6c8464152
# Model validation is done using a confusion matrix, and ROC (Receiver Operating Characteristic) curve.

# ╔═╡ Cell order:
# ╠═010e59a4-0933-11ec-005b-3314aac96bb8
# ╠═f44c2335-bed9-4bf2-ba4e-a30a7e26fe3e
# ╠═16ba1b52-6ca5-4cdb-8657-5c0ec912bb56
# ╠═b2356d86-4c41-4e13-b1e6-c9fc64f9993c
# ╠═f68d5f13-89c6-4755-b8ec-c123b87c2c35
# ╠═8f224f26-af8b-45db-b007-b542189ae20f
# ╠═196603e4-270a-45b8-afb9-4649599e78f2
# ╠═1392d02d-0258-4183-b756-969a14dfd033
# ╠═35950c0c-f87a-425f-ae96-f7555dd8eaba
# ╠═9d228706-000d-41c3-8162-9e8894eed559
# ╠═b22ed63c-b646-4e43-88bf-d1e5c754333a
# ╠═b5ed1dd2-d1af-4d22-a6da-8b5586ccac9d
# ╠═41acbbe9-2c68-493d-b7e9-9811817a4f5a
# ╠═4ac78ed0-ab59-4a05-b2d3-2e5143e781a9
# ╠═a5c2fe1d-aa59-4048-91ad-ed00262e0f57
# ╠═1a82a4df-fd0c-4c5f-b240-9672b20d017d
# ╠═75eafc32-fc5b-44aa-9d06-cf85c8b1d338
# ╠═a8dfa116-2677-4df1-acc6-738b27a864bf
# ╠═0254392c-61b0-4198-860f-f7450f1e04e1
# ╠═6b5efee9-ceb5-488c-ae8f-5adf107bf7f5
# ╠═1ec2a9c1-0033-4320-84eb-74d665bcec51
# ╠═4a4fe208-f98a-4952-b0b8-f32e8a63be61
# ╠═15265b71-0392-411d-a9b9-77b4e72a638e
# ╠═31a75e65-c4ab-43e0-8566-954b7274add2
# ╠═cbdc3eef-cd65-4101-bf7e-b649f632939a
# ╠═69ab1393-ae18-4761-a188-14f6c8464152

max_stars_repo_name: propelledanalytics/Tutorials

max_issues_repo_name: propelledanalytics/Tutorials

max_forks_repo_name: propelledanalytics/Tutorials


text:
@testset "LinearRegression" begin
    mean = LossFunction.MSE([1 2 5 10 11], [1 3 5 7 9])
    @test mean == 2.8
    x = [15.43 23.01 5.0 12.56 8.67 7.31 9.66 13.64 14.92 18.47 15.48 22.13 10.11 26.95 5.68 21.76]
    t = [170.91 160.68 129.0 159.7 155.46 140.56 153.65 159.43 164.7 169.65 160.71 173.29 159.31 171.52 138.96 165.87]
    w =LinearRegression.SGD.fit(x, t)
    p = LinearRegression.SGD.predict(x, w)
    @test 1 < w[1] < 2
    @test 100 < w[2]
    x =[15.43 23.01 5.0 12.56 8.67 7.31 9.66 13.64 14.92 18.47 15.48 22.13 10.11 26.95 5.68 21.76; 70.43 58.15 37.22 56.51 57.32 40.84 57.79 56.94 63.03 65.69 62.33 64.95 57.73 66.89 46.68 61.08]
    t = [170.91 160.68 129.0 159.7 155.46 140.56 153.65 159.43 164.7 169.65 160.71 173.29 159.31 171.52 138.96 165.87]
    w = LinearRegression.MR.fit(x, t)
    p = LinearRegression.MR.predict(x, w)
    print(p)
    w = LinearRegression.RR.fit(x, t)
    p = LinearRegression.RR.predict(x, w)
    w = LinearRegression.LR.fit(x, t)
    p = LinearRegression.LR.predict(x, w)
    println(p)
    w, m, s = LinearRegression.BFM.fit(x[1, :], t, 4)
    println("w:", w, "m:", m, "s:", s)
    p = LinearRegression.BFM.predict(x[1, :], w, m, s)
    println(p)
    mse = LossFunction.MSE(p, t)
    println(mse)
    w, m, s = LinearRegression.BFM.fit(x[1, :], t, 4, alpha = 0.1, f = "polynomial")
    println("p(polynomial):", p)
    p = LinearRegression.BFM.predict(x[1, :], w, m, s, f = "polynomial")
    println("p:", p)
    mse = LossFunction.MSE(p, t)
end


text:
##  Jacobi Operator

recA{T}(::Type{T},S::Ultraspherical,k) = (2*(k+order(S)))/(k+one(T))   # one(T) ensures we get correct type
recB{T}(::Type{T},::Ultraspherical,::) = zero(T)
recC{T}(::Type{T},S::Ultraspherical,k) = (k-one(T)+2order(S))/(k+one(T))   # one(T) ensures we get correct type

# x p_k
recα{T}(::Type{T},::Ultraspherical,::) = zero(T)
recβ{T}(::Type{T},S::Ultraspherical,k) = k/(2*(k-one(T)+order(S)))   # one(T) ensures we get correct type
recγ{T}(::Type{T},S::Ultraspherical,k) = (k-2+2order(S))/(2*(k-one(T)+order(S)))   # one(T) ensures we get correct type

## Multiplication
# these are special cases

Base.stride{U<:Chebyshev,V<:Ultraspherical}(M::ConcreteMultiplication{U,V}) =
    stride(M.f)
Base.stride{U<:Ultraspherical,V<:Chebyshev}(M::ConcreteMultiplication{U,V}) =
    stride(M.f)
Base.stride{U<:Ultraspherical,V<:Ultraspherical}(M::ConcreteMultiplication{U,V}) =
    stride(M.f)

function Multiplication{C<:Chebyshev}(f::Fun{C},sp::Ultraspherical{Int})
    if order(sp) == 1
        cfs = f.coefficients
        MultiplicationWrapper(f,
            SpaceOperator(SymToeplitzOperator(cfs/2) +
                                HankelOperator(@compat(view(cfs,3:length(cfs)))/(-2)),
                          sp,sp))

    else
        ConcreteMultiplication(f,sp)
    end
end

## Derivative

#Derivative(k::Integer,d::IntervalDomain)=Derivative(k-1:k,d)
#Derivative(d::IntervalDomain)=Derivative(1,d)

Derivative{LT,DD<:Interval}(sp::Ultraspherical{LT,DD},m::Integer) =
    ConcreteDerivative(sp,m)
function Integral{LT,DD<:Interval}(sp::Ultraspherical{LT,DD},m::Integer)
    λ = order(sp)
    if m ≤ λ
        ConcreteIntegral(sp,m)
    else # Convert up
        nsp = Ultraspherical(λ+1,domain(sp))
        IntegralWrapper(Integral(nsp,m)*Conversion(sp,nsp),m)
    end
end

rangespace{LT,DD<:Interval}(D::ConcreteDerivative{Ultraspherical{LT,DD}}) =
    Ultraspherical(order(domainspace(D))+D.order,domain(D))
bandinds{LT,DD<:Interval}(D::ConcreteDerivative{Ultraspherical{LT,DD}}) = 0,D.order
bandinds{LT,DD<:Interval}(D::ConcreteIntegral{Ultraspherical{LT,DD}}) = -D.order,0
Base.stride{LT,DD<:Interval}(D::ConcreteDerivative{Ultraspherical{LT,DD}}) = D.order

function getindex{TT,DD<:Interval,K,T}(D::ConcreteDerivative{Ultraspherical{TT,DD},K,T},
                                      k::Integer,j::Integer)
    m=D.order
    d=domain(D)
    λ=order(domainspace(D))

    if j==k+m
        (pochhammer(one(T)*λ,m)*(4./(d.b-d.a)).^m)::T
    else
        zero(T)
    end
end

## Integral

linesum{LT,DD<:Interval}(f::Fun{Ultraspherical{LT,DD}}) =
    sum(setcanonicaldomain(f))*arclength(d)/2

rangespace{LT,DD<:Interval}(D::ConcreteIntegral{Ultraspherical{LT,DD}}) =
    order(domainspace(D)) == 1 ? Chebyshev() : Ultraspherical(order(domainspace(D))-D.order,domain(D))

function getindex{LT,DD<:Interval,T}(D::ConcreteIntegral{Ultraspherical{LT,DD},T},k::Integer,j::Integer)
    m=D.order
    d=domain(D)
    λ=order(domainspace(D))
    @assert m<=λ

    if λ == 1 && k==j+1
        C = .5(d.b-d.a)
        C./(k-one(T))
    elseif λ > 1 && k==j+m
        pochhammer(one(T)*λ,-m)*(.25(d.b-d.a))^m
    else
        zero(T)
    end
end

## Conversion Operator

function Conversion(A::Chebyshev,B::Ultraspherical)
    if order(B) ≤ 1
        ConcreteConversion(A,B)
    else
        d=domain(A)
        US=Ultraspherical(order(B)-1,d)
        ConversionWrapper(TimesOperator(Conversion(US,B),
                                        Conversion(Chebyshev(d),US)))
    end
end

function Conversion(A::Ultraspherical,B::Chebyshev)
    if order(A) == 1//2
        ConcreteConversion(A,B)
    else
        error("Not implemented")
    end
end

maxspace_rule(A::Ultraspherical,B::Chebyshev) = A

function Conversion(A::Ultraspherical,B::Ultraspherical)
    a=order(A); b=order(B)
    if b==a
        ConversionWrapper(eye(A))
    elseif a<b≤a+1  || b<a≤b+1
        ConcreteConversion(A,B)
    else
        d=domain(A)
        ConversionWrapper(TimesOperator(Conversion(Ultraspherical(b-1,d),B),
                                        Conversion(A,Ultraspherical(b-1,d))))
    end
end

maxspace_rule(A::Ultraspherical,B::Ultraspherical) = order(A) > order(B) ? A : B

function getindex{DD,C<:Chebyshev,T}(M::ConcreteConversion{C,Ultraspherical{Int,DD},T},
                                     k::Integer,j::Integer)
   # order must be 1
    if k==j==1
        one(T)
    elseif k==j
        one(T)/2
    elseif j==k+2
        -one(T)/2
    else
        zero(T)
    end
end

function getindex{DD,T}(M::ConcreteConversion{Ultraspherical{Int,DD},Ultraspherical{Int,DD},T},
                            k::Integer,j::Integer)
    #  we can assume that λ==m+1
    λ=order(rangespace(M))
    c=λ-one(T)  # this supports big types
    if k==j
        c/(k - 2 + λ)
    elseif j==k+2
        -c/(k + λ)
    else
        zero(T)
    end
end

function getindex{LT,DD,T}(M::ConcreteConversion{Ultraspherical{LT,DD},Ultraspherical{LT,DD},T},
                            k::Integer,j::Integer)
    λ=order(rangespace(M))
    if order(domainspace(M))+1==λ
        c=λ-one(T)  # this supports big types
        if k==j
            c/(k - 2 + λ)
        elseif j==k+2
            -c/(k + λ)
        else
            zero(T)
        end
    else
        error("Not implemented")
    end
end

bandinds{DD}(C::ConcreteConversion{Chebyshev{DD},Ultraspherical{Int,DD}}) = 0,2  # order == 1
bandinds{DD}(C::ConcreteConversion{Ultraspherical{Int,DD},Ultraspherical{Int,DD}}) = 0,2

bandinds{LT,DD}(C::ConcreteConversion{Chebyshev{DD},Ultraspherical{LT,DD}}) =
    0,order(rangespace(C))==1?2:∞
bandinds{LT,DD}(C::ConcreteConversion{Ultraspherical{LT,DD},Chebyshev{DD}}) =
    0,order(domainspace(C))==1?2:∞

bandinds{LT1,LT2,DD}(C::ConcreteConversion{Ultraspherical{LT1,DD},Ultraspherical{LT2,DD}}) =
    0,order(domainspace(C))+1==order(rangespace(C))?2:∞

Base.stride{DD}(C::ConcreteConversion{Chebyshev{DD},Ultraspherical{Int,DD}}) = 2
Base.stride{LT1,LT2,DD}(C::ConcreteConversion{Ultraspherical{LT1,DD},Ultraspherical{LT2,DD}}) = 2

## coefficients

# return the space that has banded Conversion to the other
conversion_rule(a::Chebyshev,b::Ultraspherical{Int}) =
    if domainscompatible(a,b)
        a
    else
        NoSpace()
    end

conversion_rule{LT}(a::Ultraspherical{LT},b::Ultraspherical{LT}) =
    if domainscompatible(a,b) && isapproxinteger(order(a)-order(b))
        order(a) < order(b)?a:b
    else
        NoSpace()
    end

function coefficients(g::Vector,sp::Ultraspherical{Int},C::Chebyshev)
    if order(sp) == 1
        ultraiconversion(g)
    else
        # do one at a time
        coefficients(g,sp,Ultraspherical(1,domain(sp)),C)
    end
end
function coefficients(g::Vector,C::Chebyshev,sp::Ultraspherical)
    if order(sp) == 1
        ultraconversion(g)
    else
        # do one at a time
        coefficients(g,C,Ultraspherical(1,domain(sp)),sp)
    end
end

# TODO: include in getindex to speed up
Integral{DD<:Interval}(sp::Chebyshev{DD},m::Integer) =
    IntegralWrapper(TimesOperator([Integral(Ultraspherical(m,domain(sp)),m),
                                   Conversion(sp,Ultraspherical(m,domain(sp)))]),m)

## Non-banded conversions

function getindex{DD,LT,C<:Chebyshev,T}(M::ConcreteConversion{C,Ultraspherical{LT,DD},T},
                                     k::Integer,j::Integer)
    λ = order(rangespace(M))
    if λ == 1
        if k==j==1
            one(T)
        elseif k==j
            one(T)/2
        elseif j==k+2
            -one(T)/2
        else
            zero(T)
        end
    elseif λ == 0.5
        # Cheb-to-Leg
        if j==k==1
            one(T)
        elseif j==k
            T(sqrt(π)/(2FastTransforms.Λ(k-1)))
        elseif k < j && iseven(k-j)
            T(-(j-1)*(k-0.5)*(FastTransforms.Λ((j-k-2)/2)/(j-k))*
                            (FastTransforms.Λ((j+k-3)/2)/(j+k-1)))
        else
            zero(T)
        end
    else
        error("Not implemented")
    end
end

function getindex{DD,LT,C<:Chebyshev,T}(M::ConcreteConversion{Ultraspherical{LT,DD},C,T},
                                     k::Integer,j::Integer)
    λ = order(domainspace(M))
    if λ == 1
        # order must be 1
        if k==j==1
            one(T)
        elseif k==j
            one(T)/2
        elseif j==k+2
            -one(T)/2
        else
            zero(T)
        end
    elseif λ == 0.5
        if k==1 && isodd(j)
            T(FastTransforms.Λ((j-1)/2)^2/π)
        elseif k ≤ j && iseven(k-j)
            T(FastTransforms.Λ((j-k)/2)*FastTransforms.Λ((k+j-2)/2)*2/π)
        else
            zero(T)
        end
    else
        error("Not implemented")
    end
end

function getindex{DD,LT,LT2,T}(M::ConcreteConversion{Ultraspherical{LT,DD},
                                                     Ultraspherical{LT2,DD},T},
                                     k::Integer,j::Integer)
    λ1 = order(domainspace(M))
    λ2 = order(rangespace(M))
    if abs(λ1-λ2) < 1
        if j ≥ k && iseven(k-j)
            gamma(λ2)*(k-1+λ2)/(gamma(λ1)*gamma(λ1-λ2))*
                (gamma((j-k)/2+λ1-λ2)/gamma((j-k)/2+1))*
                (gamma((k+j-2)/2+λ1)/gamma((k+j-2)/2+λ2+1))
        else
            zero(T)
        end
    else
        error("Not implemented")
    end
end


text:
"""
    ZeroKernel()

Zero kernel.

# Definition

For inputs ``x, x'``, the zero kernel is defined as
```math
k(x, x') = 0.
```
The output type depends on ``x`` and ``x'``.

See also: [`ConstantKernel`](@ref)
"""
struct ZeroKernel <: SimpleKernel end

kappa(κ::ZeroKernel, d::T) where {T<:Real} = zero(T)

metric(::ZeroKernel) = Delta()

Base.show(io::IO, ::ZeroKernel) = print(io, "Zero Kernel")

"""
    WhiteKernel()

White noise kernel.

# Definition

For inputs ``x, x'``, the white noise kernel is defined as
```math
k(x, x') = \\delta(x, x').
```
"""
struct WhiteKernel <: SimpleKernel end

"""
    EyeKernel()

Alias of [`WhiteKernel`](@ref).
"""
const EyeKernel = WhiteKernel

kappa(κ::WhiteKernel, δₓₓ::Real) = δₓₓ

metric(::WhiteKernel) = Delta()

Base.show(io::IO, ::WhiteKernel) = print(io, "White Kernel")

"""
    ConstantKernel(; c::Real=1.0)

Kernel of constant value `c`.

# Definition

For inputs ``x, x'``, the kernel of constant value ``c \\geq 0`` is defined as
```math
k(x, x') = c.
```

See also: [`ZeroKernel`](@ref)
"""
struct ConstantKernel{Tc<:Real} <: SimpleKernel
    c::Vector{Tc}

    function ConstantKernel(; c::Real=1.0)
        @check_args(ConstantKernel, c, c >= zero(c), "c ≥ 0")
        return new{typeof(c)}([c])
    end
end

@functor ConstantKernel

kappa(κ::ConstantKernel, x::Real) = first(κ.c) * one(x)

metric(::ConstantKernel) = Delta()

Base.show(io::IO, κ::ConstantKernel) = print(io, "Constant Kernel (c = ", first(κ.c), ")")


text:
using Printf
import PyPlot

const plt = PyPlot

include("gen_rmesh_exp.jl")
include("lagrange_interp.jl")
include("radial_interp.jl")

function myfunc(x)
    return exp(-x)
end

function main()
    r_min = 0.0
    r_max = 10.0
    a = 1e9
    N = 50

    rmesh = gen_rmesh_exp(r_min, r_max, a, N)
    fr = myfunc.(rmesh)

    #NptsPlot = 5*N
    #rmesh_dense = gen_rmesh_exp(r_min, r_max, a, NptsPlot)
    #fr_interp = zeros(NptsPlot)
    #for i in 1:NptsPlot
    #    fr_interp[i] = radial_interp(rmesh, fr, )
    #end

    Vmid = zeros(N-1)
    rmid = zeros(N-1)
    for i in 1:N-1
        rmid[i] = 0.5*(rmesh[i] + rmesh[i+1])
        Vmid[i] = radial_interp( rmesh, fr, rmid[i], i+1)
        Vtrue = myfunc(rmid[i])
        err = abs(Vmid[i] - Vtrue)
        @printf("%18.10f %18.10f %18.10e\n", Vmid[i], Vtrue, err)
    end

    plt.clf()
    plt.plot(rmesh, fr, marker="o", label="data")
    plt.plot(rmid, Vmid, marker="o", label="midval")
    plt.xlim(0.0, 2.0)
    plt.legend()
    plt.savefig("IMG_test_radial_interp.png")
end

main()


text:
function parse_all(doc, server)
    ps = CSTParser.ParseState(doc._content)
    if endswith(doc._uri, ".jmd")
        doc.code.cst, ps = parse_jmd(ps, doc._content)
    else
        doc.code.cst, ps = CSTParser.parse(ps, true)
    end
    # update_includes(doc, server)
    empty!(doc.diagnostics)
    ls_diags = []

    for err in ps.errors
        if err.description == "Expected end."
            rng2 = max(0, first(err.loc)-1):last(err.loc)
            stack, offsets = StaticLint.get_stack(doc.code.cst, first(rng2))
            for i = length(stack):-1:1
                if stack[i] isa CSTParser.EXPR{T} where T <: Union{CSTParser.Begin,CSTParser.Quote,CSTParser.ModuleH,CSTParser.Function,CSTParser.Macro,CSTParser.For,CSTParser.While,CSTParser.If} && last(stack[i].args) isa CSTParser.EXPR{CSTParser.ErrorToken} && stack[i].args[end].args[1] isa CSTParser.KEYWORD
                    rng1 = offsets[i] .+ (1:stack[i].args[1].span)
                    push!(ls_diags, Diagnostic(Range(doc, rng1), 1, "Parsing error", "Julia language server", "Closing end is missing.", nothing))
                end
            end
            push!(ls_diags, Diagnostic(Range(doc, rng2), 1, "Parsing error", "Julia language server", err.description, nothing))
        else
            rng = max(0, first(err.loc)-1):last(err.loc)
            push!(ls_diags, Diagnostic(Range(doc, rng), 1, "Parsing error", "Julia language server", err.description, nothing))
        end
    end

    if server.runlinter
        if doc._runlinter
            StaticLint.pass(doc.code)
            state = StaticLint.build_bindings(server, find_root(doc, server).code);
            empty!(doc.code.rref)
            empty!(doc.code.uref)
            StaticLint.resolve_refs(doc.code.state.refs, state, doc.code.rref, doc.code.uref);

            for r in doc.code.uref
                r isa StaticLint.Reference{CSTParser.BinarySyntaxOpCall} && continue
                push!(ls_diags ,Diagnostic(r, doc))
            end
            for err in doc.code.state.linterrors
                push!(ls_diags ,Diagnostic(err, doc))
            end
        end
    end
    send(JSONRPC.Request{Val{Symbol("textDocument/publishDiagnostics")},PublishDiagnosticsParams}(nothing, PublishDiagnosticsParams(doc._uri, ls_diags)), server)
end

Diagnostic(b::StaticLint.Binding, doc) = Diagnostic(Range(doc, b.loc.offset .+ (0:b.val.span)), 2, "Unused variable", "Julia language server", "Variable declared but not used: $(string(Expr(b.val)))", nothing)
Diagnostic(r::StaticLint.Reference, doc) = Diagnostic(Range(doc, r.loc.offset .+ (0:r.val.span)), 2, "Missing variable", "Julia language server", "Use of possibly undeclared variable: $(string(Expr(r.val)))", nothing)
Diagnostic(err, doc) = Diagnostic(Range(doc, err.loc.offset .+ (0:err.val.span)), 2, "Lint error", "Julia language server", get(StaticLint.LintMessages, err.code, ""), nothing)

StaticLint.getfile(server::LanguageServerInstance, path) = server.documents[URI2(filepath2uri(path))].code
StaticLint.setfile(server::LanguageServerInstance, path, x) = server.documents[URI2(filepath2uri(path))] = x
StaticLint.is_loaded(server::LanguageServerInstance, path) = haskey(server.documents, URI2(filepath2uri(path)))
function StaticLint.load_file(server::LanguageServerInstance, path::String, index, nb, parent)
    code = read(path, String)
    uri = filepath2uri(path)
    doc = Document(uri, code, true, server, index, nb, parent)
    StaticLint.setfile(server, path, doc)
    return doc.code
end

function convert_diagnostic(h::LSDiagnostic{T}, doc::Document) where {T}
    rng = Range(doc, h.loc)
    code = 1
    Diagnostic(rng, code, string(T), string(typeof(h).name), isempty(h.message) ? string(T) : h.message, nothing)
end

function publish_diagnostics(doc::Document, server)
    ls_diags = map(doc.diagnostics) do diag
        convert_diagnostic(diag, doc)
    end
    publishDiagnosticsParams = PublishDiagnosticsParams(doc._uri, ls_diags)
    response =  JSONRPC.Request{Val{Symbol("textDocument/publishDiagnostics")},PublishDiagnosticsParams}(nothing, publishDiagnosticsParams)
    send(response, server)
end

function clear_diagnostics(uri::URI2, server)
    doc = server.documents[uri]
    empty!(doc.diagnostics)
    response =  JSONRPC.Request{Val{Symbol("textDocument/publishDiagnostics")},PublishDiagnosticsParams}(nothing, PublishDiagnosticsParams(doc._uri, Diagnostic[]))
    send(response, server)
end

function clear_diagnostics(server)
    for (uri, doc) in server.documents
        clear_diagnostics(uri, server)
    end
end


text:
module InducingPoints

import Base: rand, show

using AbstractGPs: AbstractGPs
using Clustering: kmeans!
using DataStructures
using DeterminantalPointProcesses: DPP
using Distances
using KernelFunctions
using KernelFunctions: ColVecs, RowVecs, vec_of_vecs
using LinearAlgebra#: Symmetric, Eigen, eigen, eigvals, I, logdet, diag, norm
using Random: rand, bitrand, AbstractRNG, MersenneTwister, GLOBAL_RNG
using StatsBase: Weights, sample, wsample

export AbstractInducingPointsSelectionAlg

## Generic functions

export inducingpoints
export initZ
export updateZ!, updateZ

## Offline algorithms
export KmeansAlg
export RandomSubset
export StdDPP, kDPP
export Greedy

## Online algorithms
export OIPS
export SeqDPP
export StreamKmeans
export UniGrid
export Webscale

## Custom output object
export UniformGrid

const jitt = 1e-5

abstract type AbstractInducingPointsSelectionAlg end

const AIPSA = AbstractInducingPointsSelectionAlg

abstract type OfflineInducingPointsSelectionAlg <: AIPSA end

const OffIPSA = OfflineInducingPointsSelectionAlg

abstract type OnlineInducingPointsSelectionAlg <: AIPSA end

const OnIPSA = OnlineInducingPointsSelectionAlg

## Wrapper for matrices
"""
     inducingpoints([rng::AbstractRNG], alg::AIPSA, X::AbstractVector; [kwargs...])
     inducingpoints([rng::AbstractRNG], alg::AIPSA, X::AbstractMatrix; obsdim=1, [kwargs...])

Select inducing points according to the algorithm `alg`. For some algorithms, additional keyword arguments are required.
"""
inducingpoints

function inducingpoints(
    rng::AbstractRNG, alg::AIPSA, X::AbstractMatrix; obsdim=1, kwargs...
)
    return inducingpoints(rng, alg, vec_of_vecs(X; obsdim=obsdim); kwargs...)
end

function inducingpoints(alg::AIPSA, X::AbstractMatrix; obsdim=1, kwargs...)
    return inducingpoints(GLOBAL_RNG, alg, X; obsdim=obsdim, kwargs...)
end

## Wrapper for the RNG
function inducingpoints(alg::AIPSA, X::AbstractVector; kwargs...)
    return inducingpoints(GLOBAL_RNG, alg, X; kwargs...)
end

## Online IP selection functions
@doc raw"""
    updateZ!([rng::AbstractRNG], Z::AbstractVector, alg::OnIPSA, X::AbstractVector; [kwargs...])

Update inducing points `Z` with data `X` and algorithm `alg`. Requires additional keyword arguments
for some algorithms. Also see `InducingPoints`.
"""
updateZ!

function updateZ!(Z::AbstractVector, alg::OnIPSA, X::AbstractVector; kwargs...)
    return updateZ!(GLOBAL_RNG, Z, alg, X; kwargs...)
end
# Default behavior is to simply add points
function updateZ!(
    rng::AbstractRNG, Z::AbstractVector, alg::OnIPSA, X::AbstractVector; kwargs...
)
    return add_point!(rng, Z, alg, X; kwargs...)
end

@doc raw"""
    updateZ([rng::AbstractRNG], Z::AbstractVector, alg::OnIPSA, X::AbstractVector; kwargs...)

Return new vector of inducing points `Z` with data `X` and algorithm `alg` without changing the original one
"""
updateZ

function updateZ(Z::AbstractVector, alg::OnIPSA, X::AbstractVector; kwargs...)
    return updateZ(GLOBAL_RNG, Z, alg, X; kwargs...)
end
# Default behavior is to simply add points
function updateZ(
    rng::AbstractRNG, Z::AbstractVector, alg::OnIPSA, X::AbstractVector; kwargs...
)
    return add_point(rng, Z, alg, X; kwargs...)
end

## Offline algorithms
include("offline/kmeans.jl")
include("offline/randomsubset.jl")
include("offline/stddpp.jl")
include("offline/kdpp.jl")
include("offline/greedyip.jl")

## Online algorithms
include("online/seqdpp.jl")
include("online/streamkmeans.jl")
include("online/webscale.jl")
include("online/oips.jl")
include("online/unigrid.jl")

## Utilities
include("utils.jl")

@deprecate initZ inducingpoints

end


text:
using ITensorParallel
using Test

@testset "ITensorParallel.jl" begin
    # Write your tests here.
end

max_stars_repo_name: mtfishman/ITensorParallel

max_issues_repo_name: mtfishman/ITensorParallel

max_forks_repo_name: mtfishman/ITensorParallel


text:
function NormalDistribution()
    return NormalDistribution(())
end

function NormalDistribution(arg0::jdouble, arg1::jdouble)
    return NormalDistribution((jdouble, jdouble), arg0, arg1)
end

function cumulative_probability(obj::NormalDistribution, arg0::jdouble)
    return jcall(obj, "cumulativeProbability", jdouble, (jdouble,), arg0)
end

function density(obj::NormalDistribution, arg0::jdouble)
    return jcall(obj, "density", jdouble, (jdouble,), arg0)
end

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

function get_numerical_mean(obj::NormalDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::NormalDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

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

function get_support_lower_bound(obj::NormalDistribution)
    return jcall(obj, "getSupportLowerBound", jdouble, ())
end

function get_support_upper_bound(obj::NormalDistribution)
    return jcall(obj, "getSupportUpperBound", jdouble, ())
end

function inverse_cumulative_probability(obj::NormalDistribution, arg0::jdouble)
    return jcall(obj, "inverseCumulativeProbability", jdouble, (jdouble,), arg0)
end

function is_support_connected(obj::NormalDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function log_density(obj::NormalDistribution, arg0::jdouble)
    return jcall(obj, "logDensity", jdouble, (jdouble,), arg0)
end

function probability(obj::NormalDistribution, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "probability", jdouble, (jdouble, jdouble), arg0, arg1)
end


text:
###########################################################################################################
###########################################################################################################
############################ OVERLAP LIKE G VİA CUBA NUMERICAL INTEGRATION ################################
###########################################################################################################
###########################################################################################################
export CuhreAuxiliaryG, VegasAuxiliaryG, SuaveAuxiliaryG
###########################################################################################################
function CuhreAuxiliaryG(n1, q, n2, n3, p1, p2, p3, acc)
    coef = ((p1^(n1))//Gamma(n1+1))
    function integrand(x, f)

        f[1] = ( (((1/(1 - x[1]))*(2*x[2] - 1))^q)*
                 (((1/(1 - x[1])) + (2*x[2] - 1))^n2)*
                 (((1/(1 - x[1])) - (2*x[2] - 1))^n3)*
                 exp(-p2*(1/(1 - x[1])) - p3*(2*x[2] - 1))*
                 (2/((1 - x[1])^2))
        )

    end
    res, err = cuhre(integrand, 2, 1, atol=acc, rtol=acc)
    result = NO(coef) * res

    return result
end

function VegasAuxiliaryG(n1, q, n2, n3, p1, p2, p3, acc)
    coef = ((p1^(n1))//Gamma(n1+1))
    function integrand(x, f)

        f[1] = ( (((1/(1 - x[1]))*(2*x[2] - 1))^q)*
                 (((1/(1 - x[1])) + (2*x[2] - 1))^n2)*
                 (((1/(1 - x[1])) - (2*x[2] - 1))^n3)*
                 exp(-p2*(1/(1 - x[1])) - p3*(2*x[2] - 1))*
                 (2/((1 - x[1])^2))
        )

    end
    res, err = vegas(integrand, 2, 1, atol=acc, rtol=acc)
    result = NO(coef) * res

    return result
end

function SuaveAuxiliaryG(n1, q, n2, n3, p1, p2, p3, acc)
    coef = ((p1^(n1))//Gamma(n1+1))
    function integrand(x, f)

        f[1] = ( (((1/(1 - x[1]))*(2*x[2] - 1))^q)*
                 (((1/(1 - x[1])) + (2*x[2] - 1))^n2)*
                 (((1/(1 - x[1])) - (2*x[2] - 1))^n3)*
                 exp(-p2*(1/(1 - x[1])) - p3*(2*x[2] - 1))*
                 (2/((1 - x[1])^2))
        )

    end
    res, err = suave(integrand, 2, 1, atol=acc, rtol=acc)
    result = NO(coef) * res

    return result
end


text:
# read a PBF format file
# format documentation: https://wiki.openstreetmap.org/wiki/PBF_Format, https://github.com/mapbox/osmpbf-tutorial

import .OSMPBF_pb
using ProtoBuf
using Mmap
using CodecZlib

# read a blob header (they occur many times in the file). off points to the offset into the PBF of the
# header length (which is four raw bytes before the protobuf header)
# all of these read_* functions return their data, and the offset in the file for the next read
function read_blob_header(pbf::Vector{UInt8}, off::Int64)::Tuple{OSMPBF_pb.BlobHeader,Int64}
    # read the length of the header
    hdr_len::Int32 = convert(Int32, ntoh(reinterpret(Int32, @view pbf[off:off + 3])[1]))
    off += 4

    buf = IOBuffer(@view pbf[off:off + hdr_len - 1])
    out = readproto(buf, OSMPBF_pb.BlobHeader())
    off += hdr_len
    return out, off
end

# This reads the header block which occurs once per file
function read_file_header!(pbf::Vector{UInt8}, off::Int64, blhdr::OSMPBF_pb.BlobHeader, out::OSMPBF_pb.HeaderBlock)::Int64
    # read the blob
    off = read_blob!(pbf, off, blhdr, out)
    return off + b
end

function read_blob(pbf::Vector{UInt8}, off::Int64, hdr::OSMPBF_pb.BlobHeader, T::Type)::Tuple{T, Int64}
    buf = IOBuffer(@view pbf[off : off + hdr.datasize - 1])
    blob = readproto(buf, OSMPBF_pb.Blob())

    if hasproperty(blob, :raw)
        buf = IOBuffer(blob.raw)
        out = readproto(buf, T())
    elseif hasproperty(blob, :zlib_data)
        buf = ZlibDecompressorStream(IOBuffer(blob.zlib_data))
        out = readproto(buf, T())
    else
        error("Blob not in raw or zlib format")
    end

    return out, off + hdr.datasize
end

# Parsing raw nodes (i.e. not packed into dense nodes)
# This function is untested because PBF files with regular rather than dense nodes are few and far between
function parse_nodes(nodelist::Vector{OSMPBF_pb.Node}, block::OSMPBF_pb.PrimitiveBlock, strtab::Vector{String}, cb::Function)
    @info "got $(length(nodelist)) nodes"
    for raw_node in nodelist
        nodeid = raw_node.id::Int64
        lat = 1e-9 * (block.lat_offset + (block.granularity * raw_node.lat))::Float64
        lon = 1e-9 * (block.lon_offset + (block.granularity * raw_node.lon))::Float64
        tags = map(Iterators.zip(raw_node.keys, raw_node.vals)) do t
            return strtab[t[1]] => strtab[t[2]]
        end |> Dict
        node = Node(
            nodeid,
            lat,
            lon,
            tags
        )
        cb(node)
    end

    return nothing
end

function parse_dense_nodes(dense::OSMPBF_pb.DenseNodes, block::OSMPBF_pb.PrimitiveBlock, strtab::Vector{String}, cb::Function)
    # first, parse and de-delta-code
    ids = Vector{Int64}()
    lats = Vector{Float64}()
    lons = Vector{Float64}()
    n = length(dense.id)
    sizehint!(ids, n)
    sizehint!(lats, n)
    sizehint!(lons, n)

    # read and de-delta-code ids
    current_id = 0
    for id in dense.id
        current_id += id
        push!(ids, current_id)
    end

    current_lat = 0
    for lat in dense.lat
        current_lat += lat
        push!(lats, 1e-9 * (block.lat_offset + (block.granularity * current_lat)))
    end

    current_lon = 0
    for lon in dense.lon
        current_lon += lon
        push!(lons, 1e-9 * (block.lon_offset + (block.granularity * current_lon)))
    end

    if length(dense.keys_vals) > 0
        # parse tags and yield simultaneously to save memory
        off = 1
        current_node = 1
        current_tags = Dict{String, String}()
        while off <= length(dense.keys_vals)
            if dense.keys_vals[off] == 0
                # yield this node
                Node(ids[current_node], lats[current_node], lons[current_node], current_tags) |> cb

                # prepare for next node
                current_tags = Dict{String, String}()  # can't just call empty! as caller may have a reference to tags from current node
                current_node += 1
                off += 1
                continue
            end

            key = strtab[dense.keys_vals[off]]
            val = strtab[dense.keys_vals[off + 1]]
            current_tags[key] = val
            off += 2
        end

        # yield final node if we have not already
        # the spec doesn't say whether the last element of keys_vals needs to be a zero. osmium does
        # zero-terminate these arrays, but I don't want to count on this behavior
        if current_node == n
            Node(ids[current_node], lats[current_node], lons[current_node], current_tags) |> cb
            current_node += 1
        end

        # make sure we have yielded everything
        @assert current_node == n + 1
    else
        # special case of no tags
        for node in 1:n
            Node(ids[current_node], lats[current_node], lons[current_node], Dict{String, String}()) |> cb
        end
    end

    return nothing
end

function parse_ways(ways::Vector{OSMPBF_pb.Way}, strtab::Vector{String}, cb::Function)
    for raw_way in ways
        wayid = raw_way.id::Int64

        nodes = Vector{Int64}()
        sizehint!(nodes, length(raw_way.refs))

        current_nodeid = 0
        for nodeid in raw_way.refs
            current_nodeid += nodeid
            push!(nodes, current_nodeid)
        end

        tags = map(Iterators.zip(raw_way.keys, raw_way.vals)) do t
            return strtab[t[1]] => strtab[t[2]]
        end |> Dict

        Way(wayid, nodes, tags) |> cb
    end

    return nothing
end

function parse_relations(relations::Vector{OSMPBF_pb.Relation}, strtab::Vector{String}, cb::Function)
    for raw_rel in relations
        relid = raw_rel.id::Int64
        tags = map(Iterators.zip(raw_rel.keys, raw_rel.vals)) do t
            return strtab[t[1]] => strtab[t[2]]
        end |> Dict
        members = Vector{RelationMember}()
        sizehint!(members, length(raw_rel.memids))

        memid = 0
        for i in 1:length(raw_rel.memids)
            memid += raw_rel.memids[i]  # de-delta-code
            role = strtab[raw_rel.roles_sid[i]]
            raw_type = raw_rel.types[i]
            type = raw_type == OSMPBF_pb.Relation_MemberType.NODE ? node :
                (raw_type == OSMPBF_pb.Relation_MemberType.WAY ? way : relation)
            push!(members, RelationMember(memid, type, role))
        end

        Relation(relid, members, tags) |> cb
    end
end

# hoping that the compiler is smart enough that the Any return types won't cause type instability, since the return values are not used...
function scan_pbf(pbffile; nodes::Union{Function, Missing}=missing, ways::Union{Function, Missing}=missing, relations::Union{Function, Missing}=missing)
    # read the file
    open(pbffile, "r") do stream
        pbf::Vector{UInt8} = Mmap.mmap(stream)

        off::Int64 = 1

        blhdr, off = read_blob_header(pbf, off)

        @assert blhdr._type == "OSMHeader"  "Malformed OSM PBF file, expected OSMData block but found $(blhdr._type)"

         file_header, off = read_blob(pbf, off, blhdr, OSMPBF_pb.HeaderBlock)

        @info "Reading file written by $(file_header.writingprogram)"

        # track these to avoid spammy log messages, only log once
        has_relations = false
        has_changesets = false

        # <= b/c julia starts array indices at 1, so pbf[length(pbf)] is the last byte of the file
        while (off <= length(pbf))
            # TODO should be okay to reuse blob headers right
            blhdr, off = read_blob_header(pbf, off)
            @assert blhdr._type == "OSMData"  "Malformed OSM PBF file, expected OSMData block but found $(blhdr._type)"

            block, off = read_blob(pbf, off, blhdr, OSMPBF_pb.PrimitiveBlock)

            # the primitiveblock contains several primitivegroups which may contain nodes, ways, etc.
            # but they reference string tables
            # drop index 1 (in Julia land)/index 0 (everywhere) which is always empty string and used as a delimiter
            # this also means that if a node etc say string x, strtab[x] is the correct string - no off-by-one errors
            # due to Julia starting arrays at 1.
            strtab::Vector{String} = map(String, block.stringtable.s[2:length(block.stringtable.s)])

            for grp in block.primitivegroup
                if (length(grp.nodes) > 0)
                    if !ismissing(nodes)
                        parse_nodes(grp.nodes, block, strtab, nodes)
                    end
                elseif hasproperty(grp, :dense) && length(grp.dense.id) > 0
                    if !ismissing(nodes)
                        parse_dense_nodes(grp.dense, block, strtab, nodes)
                    end
                elseif (length(grp.ways) > 0)
                    if !ismissing(ways)
                        parse_ways(grp.ways, strtab, ways)
                    end
                elseif (length(grp.relations) > 0)
                    if !ismissing(relations)
                        parse_relations(grp.relations, strtab, relations)
                    end
                elseif (length(grp.changesets) > 0)
                    if !has_changesets
                        @warn "changesets are present in $pbffile, but are not yet supported"
                        has_changesets = true  # don't warn again
                    end
                else
                    error("Malformed PBF file, PrimitiveGroup has no ways, nodes, dense nodes, relations, or changesets")
                end
            end
        end
    end
end


text:
module Spark

export SparkConfig,
       SparkContext,
       RDD,
       JuliaRDD,
       JavaRDD,
       text_file,
       parallelize,
       map,
       map_partitions,
       map_partitions_with_index,
       reduce,
       collect,
       count,
       close,
       typehint!

include("core.jl")

end


text:
"""
    approximate(f, S) -> p::LinearCombination

Approximate `f` using the approximation scheme `S`.

Possible approximation schemes are:
 - Basis with associated canonical interpolation points (e.g. `S = Chebyshev(5)`).
 - Pair of interpolation point and basis (e.g. `S = ([-1,0,1], Monomials(3))`).
 - More to come.

# Examples
```
coeffs(approximate(x->x^2, Monomials(3))) ≈ [0,0,1]
coeffs(approximate(sin, ([-π,0,π],Monomials(3))) ≈ [0,0,0]
coeffs(approximate((x1,x2)->x1*x2, Monomials(2))) ≈ [0 0; 0 1]
```
"""
function approximate end

# Broadcast single argument to number of dimensions
approximate(f, S) = approximate(f, ntuple(i->S,Val(fndims(f))))
approximate(f, S::Tuple{AbstractVector,AbstractBasis}) = invoke(approximate, Tuple{Any,Any}, f,S)

# Generic algorithm for linear approximation schemes
approximate(f, S::NTuple{N,Any}) where {N} =
    LinearCombination(
        tucker(
            grideval(f,approxpoints.(S)),
            map(Si->(f->approxtransform(Si,f)), S)
        ),
        basis.(S)
    )

approxpoints((x,B)::Tuple{AbstractVector,AbstractBasis}) = x
approxtransform((x,B)::Tuple{AbstractVector{<:Number},AbstractBasis}, f) = Matrix(B,x)\f

basis(B::AbstractBasis) = B
basis((x,B)::Tuple{AbstractVector,AbstractBasis}) = B


text:
using StaticModules, Test

using StaticModules

@staticmodule Foo begin
    x = 1
    f(y) = x^2 + 2y
end

@test :x ∈ propertynames(Foo)
@test :f ∈ propertynames(Foo)

@test Foo.x == 1

@test @with Foo begin
    f(1) == 3x
end

@test @with (;f = x -> x + 1, x = 2) begin
    f(1) == x
end

struct  Bar
    a
    b
end

@test @with Bar(1, 2) begin
    (a^2, b^2) == (1, 4)
end

@const_staticmodule X begin
    @staticmodule X begin
        X = 1
    end
end

@test @with X begin
    @with X begin
        X^2 + 1 == 2
    end
end

macro foo()
    quote
        x = 1
        y = 2
    end |> esc
end

@staticmodule Foo begin
    @foo
end

@test Foo.x == 1
@test Foo.y == 2

macro bar()
    esc(:(y))
end

@test @with Foo begin
     @bar() == 2
end

@test @with ((a=1, b=2), (b=3, c=4),(c=5, d=6)) begin
    (a,b,c,d) == (1,2,4,6)
end


text:
using LinearAlgebra
import Plots.plot as plot
# using Plots
using Random

include("PlottingWrappers.jl")
include("PointsSubstitution.jl")

"""
	function expand_matrix(input_matrix, expansion_size, last_components;do_plot=false)

Takes 'input_matrix' (an ordering matrix used for creating cliques) and and adds
2×'expansion_size' number of rows. 'last_components' are the values in original
matrix that are added last to the clique.

Results may be be plotted by setting 'do_plot=true'.
"""
function expand_matrix(input_matrix, expansion_size, last_components;do_plot=false)
	new_comp = last_components
	matrix_size = size(input_matrix,1)
	for mat_sizes = matrix_size:2:(matrix_size+2expansion_size)
		input_matrix, new_comp = add_step_to_matrix(input_matrix, new_comp)
	end
	if do_plot
    # TODO separate plotting from processing
		expand_plt_ref = plot_square_heatmap(input_matrix, 1,size(input_matrix,1);
											plt_title = "Original, size:$(matrix_size)",
											color_palete=:lightrainbow)
		display(expand_plt_ref)
	end
	return input_matrix, expand_plt_ref
end

# Shuffle matrix entries
"""
	function shuffle_matrix(input_matrix, shuffles; do_plot=false)

Takes symmetric 'input_matrix' and randomly swaps rows 'shuffles' many times.

Results may be plotted by setting 'do_plot=true'.
"""

function shuffle_matrix(input_matrix, shuffles; do_plot=false)
	matrix_size = size(input_matrix,1)
	rows = randcycle(matrix_size)
	shuffled_ord_mat = copy(input_matrix)

	for k = 1:shuffles
		# global shuffled_ord_mat, rows
		srcs, trgts = rand(rows,2)
		swap_rows!(shuffled_ord_mat, srcs, trgts)
	end

	if do_plot
    # TODO separate plotting from processing
		shuff_plt_ref = plot_square_heatmap(shuffled_ord_mat, 1,size(shuffled_ord_mat,1);
											plt_title = "Shuffled, size:$(matrix_size)",
											color_palete=:lightrainbow)
		display(shuff_plt_ref)
	end
	return input_matrix, shuff_plt_ref
end

"""
	function organize_shuff_matrix(input_matrix; do_plots=false)

Reorganizes 'input_matrix' so that values highest values in a row are positioned
next to the diagonal.

Results may be plotted by setting 'do_plot=true'.
"""
function organize_shuff_matrix(input_matrix; do_plots=false)
	unscrambled_matrix = copy(input_matrix)
	matrix_size = size(input_matrix,1)
	for k = matrix_size:-2:2
		max_row_val = findmax(unscrambled_matrix[k,:])[2]
		# put to the previous last position
		swap_rows!(unscrambled_matrix, max_row_val, k-1)
		# skip 1 row and work on next one
	end
	if do_plots
    # TODO separate plotting from processing
		reorganized_plt_ref = plot_square_heatmap(unscrambled_matrix, 1,size(unscrambled_matrix,1);
									plt_title = "unscrambled_matrix, size:$(matrix_size)",
									color_palete=:lightrainbow)
		display(reorganized_plt_ref)
		return unscrambled_matrix, reorganized_plt_ref
	else
		return unscrambled_matrix
	end
end

"""
		function order_max_vals_near_diagonal(input_matrix; do_plots=false, direction=:descending)

Orders values in 'input_matrix' so that values next to diagonal are descending
(by default).

TODO- not working-  Optionally, ascending order can be used by setting 'direction' to
':ascending'.

Results may be plotted by setting 'do_plot=true'.
"""
function order_max_vals_near_diagonal(input_matrix; do_plots=false, direction=:descending)
	# Find max values next to the diagonal
	matrix_size = size(input_matrix,1)

	if direction == :descending
		# ordering_function = findmax
		new_ord_value = -1
		iteration_values = matrix_size:-2:2
		# iteration_values = 2:2:matrix_size

	elseif direction == :ascending
		# ordering_function = findmin
		new_ord_value = findmax(input_matrix)[1]*2
		iteration_values = 2:2:matrix_size
	else
		# @error "Unknow ordering was given"
		throw("Unknow ordering was given")
	end

	reordered_matrix = copy(input_matrix)
	row_indices = 1:2:matrix_size
	col_indices = 2:2:matrix_size
	coord_set = [CartesianIndex(row_indices[k], col_indices[k]) for k=1:matrix_size÷2]
	diag_max_values = reordered_matrix[coord_set]

	for k = iteration_values
		max_val, max_ind = findmax(diag_max_values)
		(direction == :descending) ? (position = floor(k÷2)) : (position = floor(k÷2))
		diag_max_values[max_ind] = diag_max_values[position]
		diag_max_values[position] = new_ord_value
		max_ind *= 2

		swap_rows!(reordered_matrix, k, max_ind)
		swap_rows!(reordered_matrix, k-1, max_ind-1)
	end

	if do_plots
    # TODO separate plotting from processing
		reorganized_plt_ref = plot_square_heatmap(reordered_matrix, 1,size(reordered_matrix,1);
									plt_title = "reordered_matrix, size:$(matrix_size)",
									color_palete=:lightrainbow)
		display(reorganized_plt_ref)
		return reordered_matrix, reorganized_plt_ref
	else
		return reordered_matrix
	end
end

"""
	function fine_tune_matrix(input_matrix; do_plots=false)

Check if velues next to the maximal values are organized in descending order.
"""
function fine_tune_matrix(input_matrix; do_plots=false)#, direction=:descending)
	# Find max values next to the diagonal
	matrix_size = size(input_matrix,1)
	fine_tune_matrix = copy(input_matrix)

	# if direction == :descending
	# 	# ordering_function = findmax
	# 	new_ord_value = -1
	# 	iteration_values = matrix_size:-2:2
	# 	# iteration_values = 2:2:matrix_size
	#
	# elseif direction == :ascending
	# 	# ordering_function = findmin
	# 	new_ord_value = findmax(input_matrix)[1]*2
	# 	iteration_values = 2:2:matrix_size
	# else
	# 	# @error "Unknow ordering was given"
	# 	throw("Unknow ordering was given")
	# end

	for k = 2:2:matrix_size-1
		if fine_tune_matrix[k-1,k+1] > fine_tune_matrix[k,k+1]
			swap_rows!(fine_tune_matrix, k, k-1)
		end
	end

	if do_plots
    # TODO separate plotting from processing
		fine_tuned_plt_ref = plot_square_heatmap(fine_tune_matrix, 1,size(reordered_matrix,1);
										plt_title = "fine_tuned, size:$(matrix_size)",
										color_palete=:lightrainbow)
		display(fine_tuned_plt_ref)
	end
	return fine_tune_matrix, fine_tuned_plt_ref
end

# TODO separate plotting from processing
function order_max_vals_by_row_avg(input_matrix; do_plots=false)
	# Find max values next to the diagonal
	matrix_size = size(input_matrix,1)

	# Row average
	row_avg = reshape(mean(input_matrix, dims=1),(matrix_size, 1))

	# Using funciton below, because sortperm is not working on Array{Float64,2}
	sorted_rows_indexes = [findall(x->x==sort(row_avg, dims=1)[k], row_avg)[1][1] for k=1:matrix_size]
	matrix_indices = collect(range(1,matrix_size))
	# Sort indices by values (highest to lowest)
	# Create a list of indices, which corresponding valeus are ordered
	sorted_indices = sort!([1:matrix_size;],
						by=i->(sorted_rows_indexes[i],matrix_indices[i]))

	sorted_matrix = copy(input_matrix)
	for k = 1:matrix_size÷2 #iteration_values
		max_ind = sorted_indices[k]

		sorted_indices[k] = k
		sorted_indices[max_ind] = max_ind

		swap_rows!(sorted_matrix, k, max_ind)
		# swap_rows!(sorted_matrix, k-1, max_ind-1)
	end
    # TODO separate plotting from processing
	reorganized_plt_ref = plot_square_heatmap(sorted_matrix, 1,size(reordered_matrix,1);
							plt_title = "reordered_matrix, size:$(matrix_size)",
							color_palete=:lightrainbow)

    # TODO separate plotting from processing
		input_mat_plt_ref = plot_square_heatmap(input_matrix, 1,size(reordered_matrix,1);
									plt_title = "input_matrix, size:$(matrix_size)",
									color_palete=:lightrainbow)

		common_plot1 = plot(input_mat_plt_ref, reorganized_plt_ref, layout=(1,2),
																size=(800,400))

	# reordered_matrix = copy(input_matrix)
	# row_indices = 1:2:matrix_size
	# col_indices = 2:2:matrix_size
	# coord_set = [CartesianIndex(row_indices[k], col_indices[k]) for k=1:matrix_size÷2]
	#
	#
	# for k = iteration_values
	# 	max_val, max_ind = findmax(diag_max_values)
	# 	position = floor(k÷2)
	# 	diag_max_values[max_ind] = diag_max_values[position]
	# 	diag_max_values[position] = new_ord_value
	# 	max_ind *= 2
	#
	# 	swap_rows!(reordered_matrix, k, max_ind)
	# 	swap_rows!(reordered_matrix, k-1, max_ind-1)
	# end

	if do_plots
    # TODO separate plotting from processing
		reorganized_plt_ref = plot_square_heatmap(sorted_matrix, 1,size(reordered_matrix,1);
									plt_title = "reordered_matrix, size:$(matrix_size)",
									color_palete=:lightrainbow)
		display(reorganized_plt_ref)
	end
	return reordered_matrix, reorganized_plt_ref
end

function order_max_vals_near_diagonal2(input_matrix; do_final_plot=false, do_all_plots = false, direction=:descending)
	# Find max values next to the diagonal
	matrix_size = size(input_matrix,1)
	reordered_matrix = copy(input_matrix)
	reorganized_plt_ref = []

	# for every row in matrix
	for k = 1:2:matrix_size-1
		# global reordered_matrix
    # TODO separate plotting from processing
		reorganized_plt_ref_pt0 = plot_square_heatmap(reordered_matrix, 1,size(reordered_matrix,1);
								plt_title = "reordered_matrix, size:$(matrix_size)",
								color_palete=:lightrainbow)
		max_val, max_ind = findmax(reordered_matrix[k:end, k:end])
		# Take the smaller coordinate
		# (max_ind[1] < max_ind[2]) ? (target_row = max_ind[1]) : (target_row = max_ind[2])
		target_row = max_ind[1]+k-1

		reordered_matrix = swap_rows(reordered_matrix, k, target_row)
    # TODO separate plotting from processing
		reorganized_plt_ref_pt1 = plot_square_heatmap(reordered_matrix, 1,size(reordered_matrix,1);
									plt_title = "reordered_matrix, size:$(matrix_size)",
									color_palete=:lightrainbow)
	    val, second_target = findmax(reordered_matrix[k,k:end])
		second_target = second_target+k-1
		reordered_matrix = swap_rows(reordered_matrix, k+1, second_target)
		# end
		#
		#
		if do_all_plots
    # TODO separate plotting from processing
			reorganized_plt_ref_pt2 = plot_square_heatmap(reordered_matrix, 1,size(reordered_matrix,1);
										plt_title = "reordered_matrix, size:$(matrix_size)",
										color_palete=:lightrainbow)
			reorganized_plt_ref = plot(reorganized_plt_ref_pt0, reorganized_plt_ref_pt1, reorganized_plt_ref_pt2, layout=(1,3), size=(1400,400))
			display(reorganized_plt_ref)
		end
	end
	if do_final_plot
    # TODO separate plotting from processing
		reorganized_plt_ref = plot_square_heatmap(reordered_matrix, 1,size(reordered_matrix,1);
									plt_title = "reordered_matrix, size:$(matrix_size)",
									color_palete=:lightrainbow)
		# display(reorganized_plt_ref)
	else
		reorganized_plt_ref=[]
	end
	return reordered_matrix, reorganized_plt_ref
end

"""
    function get_key_for_value(d::Dict, target_value)

Returns key of the dictionary which corresponds to the given target value.
"""
function get_key_for_value(d::Dict, target_value)

    for (key, value) in d
        if value == target_value
            return key
        end
    end
end

function order_max_vals_near_diagonal3(input_matrix, ordering; direction=:descending)
    # Find max values next to the diagonal
    matrix_size = size(input_matrix,1)
    reordered_matrix = deepcopy(input_matrix)
    new_ordering = Dict()

    # for every row in matrix
    for m = 1:2:matrix_size-1
        # global reordered_matrix

        max_val, max_ind = findmax(reordered_matrix[m:end, m:end])
        # Take the smaller coordinate
        first_target = max_ind[1]+m-1
        reordered_matrix = swap_rows(reordered_matrix, m, first_target)

        # check for duplicates
        val, second_target = findmax(reordered_matrix[m,m:end])
        second_target = second_target+m-1
        if first_target == second_target
            @debug "are same"
            second_target -= 1
        end

        reordered_matrix = swap_rows(reordered_matrix, m+1, second_target)

        # find key which initially had region = first_target
        region1 = get_key_for_value(ordering, first_target)
        region2 = get_key_for_value(ordering, second_target)

        if region1 in keys(new_ordering)
            @warn "repeated"
        end
        if region2 in keys(new_ordering)
            @warn "repeated2"
        end

        new_ordering[region1] = m
        new_ordering[region2] = m +1
        println("Replaced $(region1) fom $(ordering[region1]) to $(m)")
        println("Replaced $(region2) fom $(ordering[region2]) to $(m+1)")
    end

    return reordered_matrix, new_ordering
end

##
# """
#    matrix_poling!(input_matrix; method = "avg_pooling")
#
# Takes a matrix and changes it's values to the same value, according to 'method'.
# Possible methods are:
# - 'max_pooling'- finds maximal value and replaces all values with the maximal
# 	value.
# - 'avg_pooling'- changes values to the average value
# - 'gauss_pooling'- uses gausian kernel as weights to the values in the matrix
# """
# function matrix_poling!(input_matrix::Array; method::String = "max_pooling")
# 	if method == "max_pooling"
# 		max_val = findmax(input_matrix)[1]
# 		input_matrix .= max_val
# 	end
# 	return input_matrix
# end

# matrix_poling(mat[1:3,1:3]; method = "gauss_pooling")
function matrix_poling(input_matrix::Array; method = "avg_pooling", kernel_size=3,gauss_sigma=1)
	out_matrix = copy(input_matrix)
	if method == "max_pooling"
		max_val = findmax(out_matrix)[1]
		out_matrix .= max_val
	elseif method == "avg_pooling"
		avg_val = mean(out_matrix)
		out_matrix .= floor(Int,avg_val)
	elseif method == "gauss_pooling"
		@debug "Gauss pooling"

		# gauss_kernel = ones(kernel_size,kernel_size)
		# gauss_kernel[kernel_size÷2+1,kernel_size÷2+1] = 2
		filtering_kernel = Kernel.gaussian(gauss_sigma)
		# gauss_kernel2 = imfilter(gauss_kernel, filtering_kernel)
		# gauss_kernel3 = gauss_kernel2.-findmin(gauss_kernel2)[1]/findmax(gauss_kernel2)[1]
		# gauss_kernel3 = gauss_kernel3./sum(gauss_kernel3)
		out_matrix = imfilter(out_matrix, filtering_kernel)
		# out_matrix += out_matrix.*gauss_kernel3
		out_matrix .= floor.(Int,out_matrix)
	end
	return out_matrix
end

function subsample_matrix(square_matrix::Array; subsamp_size::Int=2, method="max_pooling")
	if !issymmetric(square_matrix)
		error("Input matrix is not square")
		return
	end

	if subsamp_size == 2
		reorganize_matrix(square_matrix; subsamp_size=subsamp_size, method=method)
		return reorganize_matrix(square_matrix; subsamp_size=subsamp_size, method=method)
	elseif subsamp_size%2 == 0
		new_matrix = reorganize_matrix(square_matrix; subsamp_size=subsamp_size, method=method)
		return reorganize_matrix(new_matrix; subsamp_size=2, method=method)
	else
		new_matrix = reorganize_matrix(square_matrix; subsamp_size=subsamp_size, method=method)
		return new_matrix
	end
end

#= Should the result be overlapping or not? Options:
	- filter it as a whole image, return diagonal
	- filter subimages- the problem is that htere will be edge effect at every border of cells
	- filter subimages and assign midde value to whole patch
	- filter whole upper diagonal matrix

Plot gaussian kernel
Should we care about
=#
function reorganize_matrix(square_matrix::Array; subsamp_size::Int=2, method="max_pooling", overlap::Int=0,gauss_sigma=1)
	if method == "gauss_pooling"
		(subsamp_size >= 3) || error("Can not do gaussian pooling for area smaller than 3x3")
	end
	@debug method
	# Subsample upper half
	square_matrix2 = Float64.(copy(square_matrix))
	total_rows, total_cols = size(square_matrix)
	size_mismatch_flag = false

	# if total_rows%2 != 0
	# 	total_rows -= 1
	# end
	# if total_cols%2 != 0
	# 	total_cols -= 1
	# end

	if method == "gauss_pooling"
		square_matrix2 = zeros(Int,size(square_matrix))
		square_matrix2[1:end-1,2:end] = UpperTriangular(square_matrix[1:end-1,2:end])

		# flip matrix
		do_matrix_flip = true
		if do_matrix_flip
			square_matrix3 = zeros(Float64,size(square_matrix))
			for row in 0:total_rows-1
				for col in 0:total_cols-1
					square_matrix3[row+1,col+1] = square_matrix[end-row,end-col]
				end
			end
			square_matrix3[1:end-1,:] = square_matrix3[2:end,:]
			square_matrix3[:,2:end] = square_matrix3[:,1:end-1]
		else
			square_matrix3 = copy(square_matrix2)
			square_matrix3[1:end-1,:] = square_matrix2[2:end,:]
			square_matrix3[:,1:end-1] = square_matrix2[:,2:end]
		end

		for row in 1:total_rows
			for col in 1:row
				square_matrix2[row,col] = square_matrix3[row,col]
			end
		end

		filtering_kernel = Kernel.gaussian(gauss_sigma)
		square_matrix2 = imfilter(square_matrix2, filtering_kernel)
		square_matrix2 .= Int.(floor.(Int,square_matrix2))
	elseif method == "row_pooling"
		function gauss_func(σ,len;μ=0)
			maxv = len÷2
			 minv= -len÷2
			if len%2 == 0
				maxv-=1
			end
			x = collect(minv:1:maxv)
			return exp.(-(((x.-μ)./σ)./2).^2)./(σ*sqrt(2π))
		end
		# Take 'subsamp_size'in total in horizontal and in vertical line from
		# current matrix element
		# subsamp_size = 5
		val_range = subsamp_size÷2
		r = (subsamp_size÷2)*2
		# total_rows = 266
		# total_cols = 266
		# row = 3

		for row = 1:1:(total_rows-1)
			for col = (row+1):1:total_cols
				if row < r  &&  col <= r
					row_range = row - 1
					col_range = val_range + (val_range-row_range÷2)
				else
					row_range = val_range
				end
				if row > total_rows-r  &&  col >= total_cols-r
					col_range = total_cols - row -1
					row_range = val_range + (val_range-col_range)
				else
					col_range = val_range
				end

				r_beg = row - row_range
				r_end = row + row_range
				c_beg = col - col_range
				c_end = col + col_range

				# if r_beg < 1 && r_end > total_rows
				if r_beg < 1
					r_end += abs(r_beg)+1
					r_beg = 1
				end
				if r_end > col
					r_beg -= abs(r_end-col)
					if r_beg <1
						r_beg=1
					end

					r_end = col-1
				end
				# end # if both

				# if c_beg < row+1 && c_end > total_cols
				if c_beg < row+1
					c_end += abs(c_beg-(row+1))
					c_beg = row+1
				end
				if c_end > total_cols
					c_beg -= abs(total_rows-c_end)
					c_end = total_cols
				end
				vrange = r_beg:r_end
				try
					square_matrix2[row,col] += sum(
												square_matrix[vrange,col]
												.* gauss_func(gauss_sigma,length(vrange))
												)
												vrange = c_beg:c_end
					square_matrix2[row,col] += sum(
													square_matrix[row,c_beg:c_end] .*
													 gauss_func(gauss_sigma,length(vrange))
													)
				catch e
					@error "Failed to compute row pooling"
					@error "row" row
					@error "col" col
					square_matrix2[row,col] = 0
					break
					# error(e)
				end

			end # for col
		end # for rows

	else
		step = subsamp_size-overlap
		for row = 1:step:(total_rows-2)
			for col = (row+subsamp_size):step:total_cols
				r_beg = row
				r_end = row+subsamp_size-1
				c_beg = col
				c_end = col+subsamp_size-1

				if r_end > total_rows || c_end > total_cols
					size_mismatch_flag = true
					continue
				end
				square_matrix2[r_beg:r_end,c_beg:c_end] =
							matrix_poling(square_matrix[r_beg:r_end,c_beg:c_end]; method=method,kernel_size=subsamp_size, gauss_sigma=gauss_sigma)
			end # for col
			size_mismatch_flag && continue
		end # for rows
	end # if method

	# Copy over lower half
	for row in 2:total_rows
		for col in 1:row-1
			square_matrix2[row,col] = square_matrix2[col,row]
		end
	end

	# keep same values on diagonal
	for row in 1:total_rows
		square_matrix2[row,row] = square_matrix[row,row]
	end

	return square_matrix2
end

function pool_matrix(square_matrix::Array; method="max_pooling")
	out_matrix = copy(square_matrix)
	pool_matrix!(square_matrix; method=method)
	return out_matrix
end

"""
	add_random_patch(input_matrix; patch_size=1, total_patches=1, locations)

Takes a matrix and replaces some values with random values. Returns a new matrix
with replaced values and indicies where replacement took place.

Values can be
replaced by setting 'patch_size' to values bigger than 1. If the input matrix
is symmetric, then output matrix will be symmetric as well (values from above
diagnoal will be copied over values from below diagonal).
"""
function add_random_patch(input_matrix::Matrix; patch_size=1, total_patches=1, locations=CartesianIndex(0))
	total_rows, total_cols = size(input_matrix)
	max_row = total_rows-patch_size+1
	max_col = total_cols-patch_size+1

	output_matrix = copy(input_matrix)
	max_val = findmax(output_matrix)[1]
	min_val = findmin(output_matrix)[1]
	matrix_type = typeof(output_matrix[1])

	if patch_size>total_rows || patch_size>total_cols
		error(DimensionMismatch,": Patch size is bigger than the matrix!")
	end

	# ===
	issymmetric(input_matrix) ? (symmetrize_matrix = true) : (symmetrize_matrix = false)

	if locations == CartesianIndex(0)
		@debug "Locations were not specified- random locations will be used"
		if symmetrize_matrix
			possible_indices = findall(x->true,UpperTriangular(output_matrix))
			possible_indices = possible_indices[findall(x->x[1]<=x[2], possible_indices)]
			possible_indices = possible_indices[findall(x->x[1]<=max_row, possible_indices)]
			possible_indices = possible_indices[findall(x->x[2]<=max_col, possible_indices)]
		else
			possible_indices = possible_indices = findall(x->true,output_matrix)
		end

		tartget_indices = possible_indices[randcycle(length(possible_indices))]

	else
		wrong_indices = findall(x->x[1]>max_row || x[2]>max_col, locations)
		if isempty(wrong_indices)
			tartget_indices = locations
			total_patches = size(locations)[1]
		else
			error(DimensionMismatch,": Given indices are bigger than the matrix dimensions!")
		end
	end

	changed_indices = CartesianIndex[]
	for replacement=1:total_patches
		row = tartget_indices[replacement][1]
		col = tartget_indices[replacement][2]
		r_range = row:row+patch_size-1
		c_range = col:col+patch_size-1

		for ind in CartesianIndices((r_range,c_range))
			push!(changed_indices,ind)
		end

		new_rand_matrix = floor.(matrix_type, rand(patch_size,patch_size) .* (max_val-min_val+1) .+ min_val)

		output_matrix[r_range,c_range] .= new_rand_matrix
	end

	if symmetrize_matrix
		# Inverse second column
		changed_indices2 = [changed_indices changed_indices]
		for ind = 1:size(changed_indices)[1]
			c_ind = changed_indices2[ind,2]
			changed_indices2[ind,2] = CartesianIndex(c_ind[2],c_ind[1])
		end

		# Copy over lower half
		for row in 2:total_rows
			for col in 1:row-1
				output_matrix[row,col] = output_matrix[col,row]
			end
		end

		@debug "Returned symmetric matrix" output_matrix
		return output_matrix, changed_indices2
	else
		return output_matrix, changed_indices
	end
end

function scramble_matrix(in_matrix::Array; k::Int=2, max_iterations=-1)
	out_matrix = copy(in_matrix)
	total_rows, total_cols = size(in_matrix)
	counter = 0
	if max_iterations < 1
		max_iterations = (total_cols*(total_cols-1))/2
	end

	for row = 1:k:total_rows-k
		# @info "row:" row
		for col=total_cols:-k:row+1
			# @info "col:" col
			if row == col-k+1
				# @info "shoulbreak"
				continue
			end
			indices = collect(CartesianIndices((row:row+k-1,col-k+1:col)))
			permut_indices = shuffle(indices)
			out_matrix[indices] .= in_matrix[permut_indices]
			counter +=1
			if counter >= max_iterations
				break
			end
		end
		if counter >= max_iterations
			break
		end
	end
	return out_matrix
end


text:
module Periodize

using Cubature: hcubature
using Cuba
using JSON
using Mea.Green

II = [1.0 + 0.0im 0.0 0.0 0.0
      0.0 1.0 0.0 0.0
      0.0 0.0 1.0 0.0
      0.0 0.0 0.0 1.0]

function buildmodelvec(finsE::String, finparams::String)
    params = JSON.parsefile(finparams)
    t = 1.0
    tp = params["tp"][1]
    mu = params["mu"][1]
    (wvec, sEvec_c) = Green.readgreen_c(finsE)
    return ModelVector(t, tp, mu, wvec, sEvec_c)
end

type ModelVector
   t_::Float64 ; tp_::Float64 ; mu_::Float64
   wvec_::Array{Float64, 1} ; sEvec_c_::Array{Complex{Float64}, 3}
   cumulants_::Array{Complex{Float64}, 3}

   function ModelVector(t::Float64, tp::Float64, mu::Float64,
                        wvec::Array{Float64, 1}, sEvec_c::Array{Complex{Float64}, 3})

        cumulants = build_cumulants(wvec, mu, sEvec_c)
        return new(t, tp, mu, wvec, sEvec_c, cumulants)
    end

end

function build_cumulants(wvec::Array{Float64, 1}, mu::Float64, sEvec_c::Array{Complex{Float64}, 3})

      cumulants = zeros(Complex{Float64}, size(sEvec_c))

      for (ii, ww) in enumerate(wvec)
          tmp = zeros(Complex{Float64}, (4, 4))
          tmp[1, 1] = tmp[2, 2] = tmp[3, 3] = tmp[4, 4] = (ww + mu)
          tmp -= sEvec_c[ii, :, :]
          cumulants[ii, :, :] = inv(tmp)
      end

      return cumulants
end

type Model
    t_::Float64  ; tp_::Float64 ; mu_::Float64
    w_::Float64 ; sE_::Array{Complex{Float64}, 2}
    cumulant_::Array{Complex{Float64}, 2}

    function Model(modelvec::ModelVector, ii::Integer)
        (t, tp, mu, w, sE_c, cumulant) = (modelvec.t_, modelvec.tp_, modelvec.mu_, modelvec.wvec_[ii], modelvec.sEvec_c_[ii, :, :], modelvec.cumulants_[ii, :, :])
        return new(t, tp, mu, w, sE_c, cumulant)
    end

end

function t_value(model::Model, kx::Float64, ky::Float64) # This is t_ij(k_tilde)
    t = model.t_; tp = model.tp_
    t_val = zeros(Complex{Float64}, (4, 4))
    ex = exp(-2.0im*kx) ; emx = conj(ex)
    ey = exp(-2.0im*ky) ; emy = conj(ey)
    tloc = [0.0 -t -tp -t
            -t 0.0 -t 0.0
            -tp -t 0.0 -t
            -t 0.0 -t 0.0]

    t_val += tloc

    t_val[1, 1] += 0.0;               t_val[1, 2] += -t*ex;               t_val[1, 3] += -tp*ex*ey;    t_val[1, 4] += -t*ey
    t_val[2, 1] += -t*emx;            t_val[2, 2] += 0.0;                 t_val[2, 3] += -t*ey;        t_val[2, 4] += -tp*(emx + ey)
    t_val[3, 1] += -tp*emx*emy;       t_val[3, 2] += -t*emy;              t_val[3, 3] += 0.0;          t_val[3, 4] += -t*emx
    t_val[4, 1] += -t*emy;            t_val[4, 2] += -tp*(ex + emy);      t_val[4, 3] += -t*ex;        t_val[4, 4] += 0.0

    return (t_val)
end

function exp_k(kx::Float64, ky::Float64)

    expk = zeros(Complex{Float64}, 4) # Here exp_k[i] is in fact e**(-j*dot(k, r_i)) where r_i is a site of the cluster
    expk[1] = 1.0
    expk[2] = exp(1.0im*kx)
    expk[3] = exp(1.0im*(kx+ky))
    expk[4] = exp(1.0im*ky)
    return expk
end

function eps_0(model::Model, kx::Float64, ky::Float64)
    return (-2.0*model.t_*(cos(kx) + cos(ky)) - 2.0*model.tp_*cos(kx + ky) )
end

function hopping_test(model::Model, kx::Float64, ky::Float64)
    t = model.t_ ; tp = model.tp_
    k = [kx, ky]
    N_c = 4
    r_sites = [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
    K_sites = pi*deepcopy(r_sites)
    t_array = zeros(Complex{Float64}, (N_c, N_c))

    for i in 1:N_c
        for j in 1:N_c
            for K in K_sites
                t_array[i, j] += 1.0/N_c * exp(1.0im*dot(K + k, r_sites[i] - r_sites[j])) * eps_0(model, (K + k)...)
            end
        end
    end

    return t_array

end

function build_gf_ktilde(model::Model, kx::Float64, ky::Float64)
    return inv((model.w_ + model.mu_) * II - t_value(model, kx, ky) - model.sE_)
end

function build_gf_ktilde_inverse(model::Model, kx::Float64, ky::Float64)
    return( (model.w_ + model.mu_) * II - t_value(model, kx, ky) - model.sE_)
end

function periodize(arg::Array{Complex{Float64}, 2}, kx::Float64, ky::Float64)
    N_c = 4.0
    expk = exp_k(kx, ky)
    return (1.0/N_c * dot(expk, (arg * expk)) )
end

function make_akwgreen(model::Model)
  function akw(kk::Array{Float64, 1}) # periodize the imaginary part (Ak_w)
      #N_c = 4.0
      gf_ktilde = build_gf_ktilde(model, kk[1], kk[2])
      #expk = exp_k(kx, ky)
      #return imag(-2.0/N_c * dot(expk, (gf_ktilde * expk)) )
      return imag(-2.0*periodize(gf_ktilde, kk[1], kk[2]))
  end
  return akw
end

function make_akwcum(model::Model)
  function akwcum(kk::Array{Float64, 1}) # periodize the imaginary part (Ak_w)
      cump = periodize(model.cumulant_, kk[1], kk[2])
      return imag(-2.0*inv(inv(cump) - eps_0(model, kk[1], kk[2])))
  end
  return akwcum
end

function make_akwtrace(model::Model)
    function akwtrace(kk::Array{Float64, 1}) # periodize the imaginary part (Ak_w)
        N_c = 4.0
        return 1.0/N_c*trace(-2.0*imag(build_gf_ktilde(model, kk[1], kk[2])))
    end
    return akwtrace
end

function make_akw2green(model::Model)
    akw = make_akwgreen(model)
    function akw2(kk::Array{Float64, 1})
        return (akw(kk)^2.0)
    end
    return akw2
end

function make_akw2cum(model::Model)
    akwcum = make_akwcum(model)
    function akw2cum(kk::Array{Float64, 1})
        return (akwcum(kk)^2.0)
    end
    return akw2cum
end

function make_akw2trace(model::Model)
    function akw2trace(kk::Array{Float64, 1}) # periodize the imaginary part (Ak_w)
        N_c = 4.0
        gfktilde = build_gf_ktilde(model, kk[1], kk[2])
        return 1.0/N_c*trace(-2.0*imag(gfktilde)*-2.0*imag(gfktilde) )
    end
    return akw2trace
end

function calcintegral(modelvector::ModelVector, fct; fout_name::String="out.dat", maxevals::Int64=100000, kwargs...)

    len_sEvec_c = size(modelvector.sEvec_c_)[1]
    result = zeros(Float64, len_sEvec_c)

    #println("in calcintegral, kwargs = ", kwargs...)
    #fct = getfield(Mea.Periodize, Symbol(fctname))
    for n in 1:len_sEvec_c
        model = Model(modelvector, n)
        result[n] = (2.0*pi)^(-2.0)*hcubature(fct(model; kwargs...), (-pi, -pi), (pi, pi), reltol=1.49e-8, abstol=1.49e-8, maxevals=maxevals)[1]
    end

    result_out = hcat(modelvector.wvec_, result)
    writedlm(fout_name, result_out, " ")
    return result_out
end

function calcintegral_cuba(modelvector::ModelVector, fct; fout_name::String="out.dat", maxevals::Int64=100000, kwargs...)

    len_sEvec_c = size(modelvector.sEvec_c_)[1]
    result = zeros(Float64, len_sEvec_c)

    #fct = getfield(Mea.Periodize, Symbol(fctname))
    for n in 1:len_sEvec_c
        model = Model(modelvector, n)
        result[n] = (2.0*pi)^(-2.0)*divonne(fct(model; kwargs...), 2, 1, reltol=1.49e-8, abstol=1.49e-8, maxevals=maxevals).integral[1]
    end

    result_out = hcat(modelvector.wvec_, result)
    writedlm(fout_name, result_out, " ")
    return result_out
end

 function calcdos(modelvector::ModelVector; fout_name::String="dos.dat", maxevals::Int64=100000)
     dos_out = calcintegral(modelvector, make_akwgreen, fout_name=fout_name, maxevals=maxevals)
     return dos_out
 end

 function calcdos2(modelvector::ModelVector; fout_name::String="dos2.dat", maxevals::Int64=100000)
     dos2_out = calcintegral(modelvector, make_akw2green, fout_name=fout_name, maxevals=maxevals)
     return dos2_out
 end

end