text:
using Distributed
addprocs(["compute-0-4","compute-0-5","compute-0-7","compute-0-9","compute-0-10"])
nprocs()
@everywhere println(Threads.nthreads())
@everywhere using FRGRealTime
@everywhere import FRGRealTime.ThresholdFunctions as TF
@everywhere using DifferentialEquations
@everywhere using HCubature
@everywhere using Dierckx
using LaTeXStrings
using CairoMakie
@everywhere using BenchmarkTools
@everywhere path="/home/tyy/Documents/FRGRealTime.jl/Applications/O(N) t-chanel Scheme2/"
@everywhere include(joinpath(path,"struct.jl"))
@everywhere include(joinpath(path,"flow.jl"))

include(joinpath(path,"solver.jl"))

@everywhere config_lpa = ODEConfig(dtmax = 100.0)
@everywhere config_spec = ODEConfig(dtmax = 0.1, progress = false, dense = false, save_on = false, save_start = false, save_end = true)

@everywhere lpakrang = RGscale(800.0, 1e-10)

@everywhere function tchanelLPASolve(
    Temper,
    krang::RGscale,
    massflowfun,
    lambdaflowfun;
    u0,
    config::ODEConfig
)
    function flow_LPA(du, u, p, k)
        #display(plot(x->-h2(x)[2],k,Λ))
        #println(k, " m2=", -u[2], " λ=", u[1])
        du[1] = massflowfun(k, u[1], p, u[2])
        du[2] = lambdaflowfun(k, u[1], p, u[2])
    end
    alg = Vern9()
    prob = ODEProblem(flow_LPA, u0, krang.kspan, Temper)
    solve(
        prob,
        alg,
        progress = config.progress,
        abstol = config.atol,
        reltol = config.rtol,
        adaptive = config.adaptive,
        dtmax = config.dtmax
    )
end

@time @everywhere sol1 = tchanelLPASolve(
    150.0, lpakrang,
    massflow2,
    lambdaflow2;
    u0 = [-116585.9, 6 * 8.0], config = config_lpa)

@everywhere solmassfun=Spline1D(reverse(sol1.t),reverse(hcat(sol1.u...)[1,:]))
@everywhere sollambdafun=Spline1D(reverse(sol1.t),reverse(hcat(sol1.u...)[2,:]))

testa=rand(1000)

@benchmark solmassfun.(testa)
myfun1=x->sol1(x)[1]
@benchmark myfun1.(testa)

plot(k -> sol1(k)[1], 0.001, 1.0)

@time tsolparallel1 = tchanelSolveParallel(
    145.0,
    1.0,
    1000.0,
    x -> sol1(x)[1],
    Imlambdaflow1,
    Relambdaflow1,
    4000,
    0.2,
    u0 = [-0.0001, 5 * 7.0],
    dtmax = 0.05,
)

@time tsolTwoPointparallel1 = tchanelSolveTwoPointParallel(
    150.0,
    1.0,
    800.0,
    x -> sol1(x)[1],
    x -> sol1(x)[2],
    Imlambdaflow2,
    Relambdaflow2,
    100,
    0.05,
    u0 = [-0.00001, 6 * 8.0],
    dtmax = 0.04,
)

@time tsolTwoPointparallel2 = tchanelSolveTwoPointParallel2(
    100.0,
    150.0,
    1.0,
    800.0,
    x -> sol1(x)[1],
    x -> sol1(x)[2],
    Imlambdaflow3,
    Relambdaflow3,
    1600,
    0.2,
    u0 = [-0.00001, 6 * 8.0],
)

@time tsolTwoPointparallel2 = tchanelSolveTwoPointParallel3(
    100.0,
    150.0, lpakrang,
    x -> sol1(x)[1],
    x -> sol1(x)[2],
    Imlambdaflow3,
    Relambdaflow3,
    800,
    0.2,
    u0 = [-0.00001, 6 * 8.0],
    config=config_spec,
)

config_spec.

config_spec

tsolTwoPointparallel2 = [[0.0, 0.0] for i = 1:400]

p0range = 1.0:1.0:400.0

@time Threads.@threads for i = 1:64
    tsolTwoPointparallel2[i] = tchanelSolveTwoPointParallel3(
        p0range[i],
        150.0, lpakrang,
        x -> sol1(x)[1],
        x -> sol1(x)[2],
        Imlambdaflow3,
        Relambdaflow3,
        800,
        0.2,
        u0 = [-0.00001, 6 * 8.0],)
end

config_spec.dtmax=50.0
config_spec.progress=false
config_spec.atol=1e-12
config_spec.rtol=1e-12
@time fourpointdatamk = tchanelSolveFourPointParallel(
    0.1,
    150.0,
    lpakrang,
    solmassfun,
    sollambdafun,
    Imlambdaflow3_ac,
    Relambdaflow3_ac,
    1600,
    0.005,
    u0 = [-0.000001, 6 * 8.0],
    config = config_spec,
)

writedlm(joinpath(path,"Data/k=100mk.dat"),hcat(fourpointdatamk...))

config_spec.maxiters=10^10

plot(fourpointdata[1],log.(10,abs.(-fourpointdata[2])),label =L"\mathrm{Im}\lambda",title="k=200MeV",xlabel=L"p_0")

testfun1=Spline1D(fourpointdata[1],log.(10,abs.(-fourpointdata[2])))
q0=Epi(5.0,sol1(5.0)[1])

plot!([q0],[testfun1(q0)],seriestype = :scatter,label = L"q_0")

plot!([3q0],[testfun1(3q0)],seriestype = :scatter,label = L"3q_0")

writedlm("/home/wjfu1/tyy/O(4)Data/vertex/full/k=200T=150.dat",hcat(fourpointdata...))

writedlm("/home/wjfu1/tyy/O(4)Data/vertex/full/k=200.cofig",config_spec)

config_spec

config_spec.maxiters=10^10

# yaxis = :log,
plot!(fourpointdata[1],-fourpointdata[2],label =L"\mathrm{Im}\lambda",title="k=200MeV",xlabel=L"p_0",yaxis = :log)

plot(p0->TF.ImFb2(p0 - 10.0, 10.0, sol1(10.0)[1], 150.0, 0.1),1.0,100.0)

testfun1=Spline1D(fourpointdata[1],fourpointdata[2])
q0=Epi(200.0,sol1(200.0)[1])

plot!([q0],[-testfun1(q0)],yaxis = :log,seriestype = :scatter,label = L"q_0")

plot!([3q0],[-testfun1(3q0)],yaxis = :log,seriestype = :scatter,label = L"3q_0")

tsolTwoPointparallel2

plot(hcat(tsolTwoPointparallel2...)'[:, 1])

imvec = hcat((tsolTwoPointparallel2)...)'[:, 1]
revec = hcat((tsolTwoPointparallel2)...)'[:, 2]

lambda_omega = range(0.0, stop = 400.0, length = 100)

plot(p0range, revec)

plot!(lambda_omega, abs.(revec), yaxis = :log)

plot(p0range, -imvec ./ (imvec .^ 2 + revec .^ 2), yaxis = :log)
tsolparallel1[3](40.0)

plot(k -> sol1(k)[1], 1.0, 800.0)

plot(k -> -tsol1(k)[1], 1, 1000.0)

plot(k -> tsol1(k)[2], 1.0, 800.0, yaxis = :log)


text:
using Documenter
using CuckooFilters

makedocs(
    sitename = "CuckooFilters",
    format = Documenter.HTML(),
    modules = [CuckooFilters]
)

deploydocs(
    repo = "github.com/kernelmethod/CuckooFilters.jl.git",
    devbranch = "main",
)

ext: jl


text:
using Mimi
using OptiMimi
using Base.Test

using ForwardDiff

include("testlp.jl")

# Differentiation-free optimization
# Create quadratic component
@defcomp quad1 begin
    regions = Index()

    # The x-value of the maximum of the quadratic
    maximum = Parameter(index=[regions])

    # The x-value to evaluate the quadratic
    input = Parameter(index=[regions])

    # The y-value of the quadratic at the x-value
    value = Variable(index=[regions])
end

function run_timestep(state::quad1, t)
    v = state.Variables
    p = state.Parameters

    v.value = -(p.input - p.maximum).^2
end

# Prepare model
model1 = Model()
setindex(model1, :time, 1)
setindex(model1, :regions, 2)
addcomponent(model1, quad1)
setparameter(model1, :quad1, :maximum, [2., 10.])
setparameter(model1, :quad1, :input, [0., 0.])

objective1(model::Model) = sum(model[:quad1, :value])

optprob = problem(model1, [:quad1], [:input], [0.], [100.0], objective1)
(maxf, maxx) = solution(optprob, () -> [0., 0.])

@test_approx_eq_eps maxf 0 1e-2
@test_approx_eq_eps maxx[1] 2 1e-2
@test_approx_eq_eps maxx[2] 10 1e-2

# Automatic differentiation

# Create quadratic component
@defcomp quad2 begin
    regions = Index()

    # The x-value of the maximum of the quadratic
    maximum = Parameter(index=[regions])

    # The x-value to evaluate the quadratic
    input = Parameter(index=[regions])

    # The y-value of the quadratic at the x-value
    value = Variable(index=[regions])
end

function run_timestep(state::quad2, t)
    v = state.Variables
    p = state.Parameters

    v.value = -(p.input - p.maximum).^2
end

# Prepare model
model2 = Model(Number)
setindex(model2, :time, 1)
setindex(model2, :regions, 2)
addcomponent(model2, quad2)
setparameter(model2, :quad2, :maximum, convert(Array{Number,1}, [2., 10.]))
setparameter(model2, :quad2, :input, convert(Array{Number,1}, [0., 0.]))

objective2(model::Model) = sum(model[:quad2, :value])

# Test the translation to a simple objective function
uo = unaryobjective(model2, [:quad2], [:input], objective2)
guo(xx) = ForwardDiff.gradient(uo, xx)
guos = guo([0., 0.])
@test_approx_eq_eps guos[1] 4 1e-2
@test_approx_eq_eps guos[2] 20 1e-2

# Optimize
optprob = problem(model2, [:quad2], [:input], [0.], [100.0], objective2)
(maxf, maxx) = solution(optprob, () -> [0., 0.])

@test_approx_eq_eps maxf 0 1e-2
@test_approx_eq_eps maxx[1] 2 1e-2
@test_approx_eq_eps maxx[2] 10 1e-2


text:
using QWNNM
using Test

nSig  = 10
Img_O = loadImg("lena.png")
Img_N = loadImg("lenaN.png")

@test psnr(addNoise(Img_O, 10), Img_O) ≈ psnr(Img_N, Img_O)


text:
module test_h3_api_coordinate_systems

using Test
using H3.API # Vec2d Vec3d CoordIJK FaceIJK ijToIjk ijkToHex2d ijkToIj ijkDistance ijkNormalize h3ToLocalIjk h3ToFaceIjk localIjkToH3 faceIjkToH3 hex2dToCoordIJK geoToVec3d geoToFaceIjk

@test ijkToHex2d(CoordIJK(0, 0, 0)) == Vec2d(0.0, 0.0)
@test ijkToHex2d(CoordIJK(0, 0, 1)) == Vec2d(-0.5, -0.8660254037844386)
@test ijkToHex2d(CoordIJK(0, 1, 0)) == Vec2d(-0.5, 0.8660254037844386)
@test ijkToHex2d(CoordIJK(0, 1, 1)) == Vec2d(-1.0, 0.0)
@test ijkToHex2d(CoordIJK(1, 0, 0)) == Vec2d(1.0, 0.0)
@test ijkToHex2d(CoordIJK(1, 0, 1)) == Vec2d(0.5, -0.8660254037844386)
@test ijkToHex2d(CoordIJK(1, 1, 0)) == Vec2d(0.5, 0.8660254037844386)

h = 0x85283473fffffff
ijk = CoordIJK(0, 37, 120)
@test ijkToHex2d(ijk) == Vec2d(-78.5, -71.88010851410841)
@test ijkToIj(ijk) == CoordIJ(-120, -83)
@test ijToIjk(CoordIJ(-120, -83)) == ijk
@test hex2dToCoordIJK(Vec2d(-78.5, -71.88010851410841)) == ijk
@test localIjkToH3(h, ijk) == 0x0851c29c7fffffff
@test h3ToLocalIjk(h, 0x0851c29c7fffffff) == ijk

faceijk = FaceIJK(7, ijk)
@test faceijk.face == 7
@test faceijk.coord == ijk
@test h3ToFaceIjk(h) == faceijk
@test faceIjkToH3(faceijk, 5) == h

geo = h3ToGeo(h)
@test geo == GeoCoord(0.6518070561696664, -2.128889370371519)
@test geoToH3(geo, 5) == h
@test geoToFaceIjk(geo, 5) == faceijk
@test geoToVec3d(geo) == Vec3d(-0.42100191417850585, -0.674362461213141, 0.6066239849895393)

ijk = CoordIJK(-1, 37, 120)
@test ijkDistance(ijk, CoordIJK(-1, 37, 120)) == 0
@test ijkDistance(ijk, CoordIJK(0, 37, 120))  == 1
@test ijkDistance(ijk, CoordIJK(1, 37, 120))  == 2
@test ijkNormalize(ijk) == CoordIJK(0, 38, 121)

end # module test_h3_api_coordinate_systems


text:
module CrossValidation

	include("./MatrixNetworks.jl")
	include("./NetworkBasedInference.jl")
	using Random
	using LinearAlgebra
	using .MatrixNetworks
	using .NetworkBasedInference

	export generateFolds, split_adjmat, alaimo2013_split_adjmat, cheng2012_split_adjmat
	export SDTNBI_LOOCV

	"""
		_10_fold_over_edges(A::Matrix, rng::Number = rand(1:typemax(Int64)))::Dict

	Generate 10 groups of edges for cross-validation.

	# Examples
	```julia-repl
	julia> A = [ 0 0 1 0 1 ; 0 0 0 0 1 ; 1 0 0 0 0 ; 0 0 0 0 0 ; 1 1 0 0 0 ]
	5×5 Array{Int64,2}:
	0  0  1  0  1
	0  0  0  0  1
	1  0  0  0  0
	0  0  0  0  0
	1  1  0  0  0
	julia> generateFolds_edges_10fold(A)
	Dict{Tuple{Int64,Int64},Int64} with 25 entries:
	  (3, 4) => 3
  	  (4, 3) => 3
  	  (2, 3) => 1
  	  (4, 4) => 7
  	  (2, 2) => 9
  	  (5, 4) => 6
  	  (2, 1) => 9
  	  (2, 5) => 2
  	  (3, 1) => 8
  	  (1, 4) => 2
  	  (3, 3) => 2
  	  (4, 5) => 0
  	  (1, 3) => 1
  	  (3, 2) => 3
  	  (5, 2) => 0
  	  (2, 4) => 5
  	  (1, 1) => 4
  	  (1, 2) => 5
  	  (5, 1) => 4
  	  (5, 5) => 5
  	  (5, 3) => 7
  	  (4, 1) => 4
  	  (1, 5) => 6
  	  (3, 5) => 1
  	  (4, 2) => 8
	```

	# Arguments
	- `M::Matrix`: rectangular M × N adjacency matrix.
	- `rng::Int`: seed for random number generator.

	"""
	function _10_fold_over_edges(A::Matrix, rng::Number = rand(1:typemax(Int64)))::Dict
		# Get all possible edges in graph
		M, N		= size(A)
		all_edges	= [(i,j) for i in 1:M for j in 1:N]

		# Assign fold randomly to all possible edges in graph
		seed			= Random.seed!(rng)
		shuffled_edges	= shuffle(seed, all_edges)
		edges_fold		= Dict((i,j) => mod(f,10) for (f,(i,j))
							   in enumerate(shuffled_edges))

		return edges_fold
	end

	"""
		generateFolds_edges_loo(A::Matrix, rng::Number = rand())::Dict

	Generate leave-one-out groups of edges for cross-validation.

	# Examples
	```julia-repl
	julia> A = [ 0 0 1 0 1 ; 0 0 0 0 1 ; 1 0 0 0 0 ; 0 0 0 0 0 ; 1 1 0 0 0 ]
	5×5 Array{Int64,2}:
	0  0  1  0  1
	0  0  0  0  1
	1  0  0  0  0
	0  0  0  0  0
	1  1  0  0  0
	julia> generateFolds_edges_loo(A, randint())
	Dict{Tuple{Int64,Int64},Int64} with 25 entries:
	  (3, 4) => 24
	  (4, 3) => 20
	  (2, 3) => 16
	  (2, 2) => 15
	  (4, 4) => 19
	  (5, 4) => 2
	  (2, 1) => 6
	  (2, 5) => 13
	  (3, 1) => 4
	  (1, 4) => 22
	  (3, 3) => 1
	  (4, 5) => 21
	  (1, 3) => 23
	  (5, 2) => 17
	  (3, 2) => 25
	  (2, 4) => 7
	  (1, 1) => 8
	  (1, 2) => 9
	  (5, 1) => 5
	  (5, 3) => 10
	  (4, 1) => 11
	  (5, 5) => 14
	  (1, 5) => 3
	  (4, 2) => 12
	  (3, 5) => 18
	```

	# Arguments
	- `M::Matrix`: rectangular M × N adjacency matrix.
	- `rng::Int`: seed for random number generator.

	"""
	function _leave_one_out_over_edges(A::Matrix, rng::Number = rand(1:typemax(Int64)))::Dict
		# Get all possible edges in graph
		M, N		= size(A)
		all_edges	= [(i,j) for i in 1:M for j in 1:N]

		# Assign fold randomly to all possible edges in graph
		seed			= Random.seed!(rng)
		shuffled_edges	= shuffle(seed, all_edges)
		edges_fold		= Dict((i,j) => f for (f,(i,j)) in enumerate(shuffled_edges))

		return edges_fold
	end

	"""
		generateFolds_nodes_loo(A::Matrix, rng::Number = rand())::Dict

	Generate leave-one-out groups of edges for cross-validation.

	# Examples
	```julia-repl
	julia> A = [ 0 0 1 0 1 ; 0 0 0 0 1 ; 1 0 0 0 0 ; 0 0 0 0 0 ; 1 1 0 0 0 ]
	5×5 Array{Int64,2}:
	0  0  1  0  1
	0  0  0  0  1
	1  0  0  0  0
	0  0  0  0  0
	1  1  0  0  0
	```

	# Arguments
	- `M::Matrix`: rectangular M × N adjacency matrix.
	- `rng::Int`: seed for random number generator.

	"""
	function _10_fold_over_nodes(A::Matrix, rng::Number = rand(1:typemax(Int64)))::Dict
		# Get all possible edges in graph
		M, N		= size(A)
		all_edges	= [(i,j) for i in 1:M for j in 1:N]

		# Assign fold randomly to all possible edges in graph
		seed			= Random.seed!(rng)
		shuffled_edges	= shuffle(seed, all_edges)
		edges_fold		= Dict((i,j) => f for (f,(i,j)) in enumerate(shuffled_edges))

		return edges_fold
	end

	function _leave_one_out_over_nodes(A::Matrix, rng::Number)::Dict
		# Get all possible edges in graph
		M, N		= size(A)
		all_edges	= [(i,j) for i in 1:M for j in 1:N]

		# Assign fold randomly to all possible edges in graph
		seed			= Random.seed!(rng)
		shuffled_edges	= shuffle(seed, all_edges)
		edges_fold		= Dict((i,j) => f for (f,(i,j)) in enumerate(shuffled_edges))

		return edges_fold
	end

	"""
		generateFolds(A::Matrix, method::String = "10-fold over edges")::Dict

	Wrapper function for generating fold groups
	"""
	function generateFolds(A::Matrix, method::String = "10-fold over edges", rng::Number = rand(1:typemax(Int64)))::Dict
		if method == "10-fold over edges"
			return _10_fold_over_edges(A, rng)
		elseif method == "Leave-one-out over edges"
			return _leave_one_out_over_edges(A, rng)
		elseif method == "10-fold over nodes"
			return _10_fold_over_nodes(A, rng)
		elseif method == "Leave-one-out over nodes"
			return _leave_one_out_over_nodes(A, rng)
		else
			error("$method not available. Available methods for fold generation:
				  - 10-fold over edges
				  - 10-fold over nodes
				  - Leave-one-out over edges
				  - Leave-one-out over nodes")
		end
	end

	"""

	"""
	function split_adjmat(A::Matrix, groups::Dict, fold::Int)
		# Create empty matrices
		M, N = size(A)
		training_A = deepcopy(A)
		test_A = zeros(Int64, M, N)

		# Get edges for given group
		test_edges = [e for (e,f) in groups if f == fold]

		# Populate matrices
		for (i,j) in test_edges
			training_A[i,j] = 0
			test_A[i,j] = A[i,j]
		end

		return training_A, test_edges
	end

	"""
		alaimo2013_split_adjmat(A::Matrix, groups::Dict, fold::Int)

	Generate test and training sets for cross-validation based on Alaimo et al. (2013)

	"""
	#=
	function alaimo2013_split_adjmat(A::Matrix, groups::Dict, fold::Int)
		# Create empty matrices
		M, N = size(A)
		training_A = deepcopy(A)
		test_A = zeros(Int64, M, N)

		# Get edges for given group
		test_edges = [e for (e,f) in groups if f == fold]

		# Populate matrices
		for (i,j) in test_edges
			j_adjacency = deepcopy(training_A[:,j])
			i_adjacency = deepcopy(training_A[i,:])
			j_adjacency[i] = 0
			i_adjacency[j] = 0

			if sum(j_adjacency) == 0 || sum(i_adjacency) == 0
				continue
				#filter!(edge->edge≠(i,j),test_edges)
			else
				training_A[i,j] = 0
				test_A[i,j] = A[i,j]
			end
		end

		return training_A, test_A, test_edges
	=#
	"""
	cheng2012_split_adjmat(A::Matrix, groups::Dict, fold::Int)

	Generate test and training sets for cross-validation based on Cheng et al. (2012)

	"""
	function cheng2012_split_adjmat(A::Matrix, groups::Dict, fold::Int)
		# Create empty matrices
		M, N		= size(A)
		training_A	= deepcopy(A)
		test_A		= zeros(Int64, M, N)

		# Get edges for given group
		test_edges_0 = [e for (e,f) in groups if f == fold]
		test_edges_1 = []

		# Populate matrices
		for (i,j) in test_edges_0
			training_A[i,j] = 0
			test_A[i,j] = A[i,j]
		end

		# Check if any vertex is isolated
		for i in 1:1:M
			if sum(training_A[i,:]) != 0
				edges_to_keep = [(i1, j1) for (i1, j1) in test_edges_0 if i1 == i]
				for e in edges_to_keep
					push!(test_edges_1, e)
				end
			end
		end

		for j in 1:1:N
			if sum(training_A[:,j]) != 0
				edges_to_keep = [(i1, j1) for (i1, j1) in test_edges_0 if j1 == j]
				for e in edges_to_keep
					push!(test_edges_1, e)
				end
			end
		end

		return training_A, test_A, test_edges_0, Set(test_edges_1)
	end

	function alaimo2013_split_adjmat(A::Matrix, groups::Dict, fold::Int)
		# Create empty matrices
		M, N		= size(A)
		training_A	= deepcopy(A)
		test_A		= zeros(Int64, M, N)

		# Get edges for given group
		test_edges = [e for (e,f) in groups if f == fold]

		# Populate matrices
		for (i,j) in test_edges
			training_A[i,j] = 0
			test_A[i,j] = A[i,j]
		end

		# Check if any vertex is isolated
		for i in 1:1:M
			if sum(training_A[i,:]) == 0
				i_neighbours = neighbours(A, i, "row")
				training_A[i, rand(i_neighbours)] = 1
			end
		end

		for j in 1:1:N
			if sum(training_A[:,j]) == 0
				j_neighbours = neighbours(A, j, "column")
				training_A[rand(j_neighbours), j] = 1
			end
		end

		return training_A, test_A, test_edges
	end

	function SDTNBI_LOOCV(DS::Matrix, DT::Matrix, C::Int, k::Int = 2)
		# Get number of nodes per type from adjacency matrices
		Nc, Nd, Ns, Nt = (1, size(DS, 1), size(DS, 2), size(DT, 2))

		# Create matrix A
		Mcc = zeros(Nc, Nc)
		Mcd = zeros(1, Nd-1)
		Mcs = DS[C,:]'
		Mct = zeros(1, Nt)

		Mdc = Mcd'
		Mdd = zeros(Nd-1, Nd-1)
		Mds = DS[1:end .!= C,:]
		Mdt = DT[1:end .!= C,:]

		Msc = Mcs'
		Msd = Mds'
		Mss = zeros(Ns, Ns)
		Mst = zeros(Ns, Nt)

		Mtc = Mct'
		Mtd = Mdt'
		Mts = Mst'
		Mtt = zeros(Nt, Nt)

		A = [ Mcc Mcd Mcs Mct ;
			  Mdc Mdd Mds Mdt ;
			  Msc Msd Mss Mst ;
			  Mtc Mtd Mts Mtt ]

		# Create matrix B
		B = deepcopy(A)
		B[1,:] = B[1,:] .= 0
		B[:,1] = B[:,1] .= 0

		# Calculate transfer matrix W
		W = SDTNBI(Symmetric(B))
		F = A * W^k
		R = F[1, Ns+Nd+1:end]

		return R
	end

end


text:
struct PositionAroundPoint{M,N,S,m} <: TaskMap{M,N,S}
    position_center::SVector{m,S}
end

PositionAroundPoint{M,N}(c::SVector{m,S}) where {M,N,S,m} = PositionAroundPoint{M,N,S,dim(M)}(c)

struct PositionAroundPointT{M,N,S,m} <: TaskMapT{M,N,S}
    base_map::PositionAroundPoint{M,N,S,m}
end
PositionAroundPointT{M,N}(c::SVector{m,S}) where {M,N,S,m} =
    PositionAroundPointT{M,N,S,m}(PositionAroundPoint{M,N}(c))

task_map_emb(::EmbRep, ::EmbRep, xm, task_map::PositionAroundPoint{ℝ{m},ℝ{m}}) where m =
    xm - task_map.position_center


text:
include("get_div_grad.jl")

function test_minres()
  minres_tol = 1.0e-6

  # 1. Symmetric and positive definite systems.
  #
  # Cubic spline matrix.
  n = 10;
  A = spdiagm(-1 => ones(n-1), 0 => 4*ones(n), 1 => ones(n-1))
  b = A * [1:n;]

  (x, stats) = minres(A, b)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("MINRES: Relative residual: %8.1e\n", resid)
  @test(resid <= minres_tol)
  @test(stats.solved)

  # radius = 0.75 * norm(x)
  # (x, stats) = minres(A, b, radius=radius, itmax=10)
  # show(stats)
  # @test(stats.solved)
  # @test(abs(radius - norm(x)) <= minres_tol * radius)

  # Symmetric indefinite variant.
  A = A - 3 * I
  b = A * [1:n;]
  (x, stats) = minres(A, b, itmax=10)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("MINRES: Relative residual: %8.1e\n", resid)
  @test(resid <= minres_tol)
  @test(stats.solved)

  # Code coverage.
  (x, stats) = minres(Matrix(A), b)
  show(stats)

  # Sparse Laplacian.
  A = get_div_grad(16, 16, 16)
  b = ones(size(A, 1))
  (x, stats) = minres(A, b)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("MINRES: Relative residual: %8.1e\n", resid)
  @test(resid <= minres_tol)
  @test(stats.solved)

  # radius = 0.75 * norm(x)
  # (x, stats) = minres(A, b, radius=radius, itmax=10)
  # show(stats)
  # @test(stats.solved)
  # @test(abs(radius - norm(x)) <= minres_tol * radius)

  # Symmetric indefinite variant, almost singular.
  A = A - 5 * I
  (x, stats) = minres(A, b)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("MINRES: Relative residual: %8.1e\n", resid)
  @test(resid <= 100 * minres_tol)
  @test(stats.solved)

  # Test b == 0
  (x, stats) = minres(A, zeros(size(A,1)))
  @test x == zeros(size(A,1))
  @test stats.status == "x = 0 is a zero-residual solution"

  # Test integer values
  A = spdiagm(-1 => ones(Int, n-1), 0 => 4*ones(Int, n), 1 => ones(Int, n-1))
  b = A * [1:n;]
  (x, stats) = minres(A, b)
  @test stats.solved
end

test_minres()


text:
import Base.isequal;
import Base.==;

function PQEcolaPoint.isequal(p1::PQContractPremiumLatticePoint, p2::PQContractPremiumLatticePoint)::Bool

    if ((p1.i == p2.i) && (p1.j == p2.j) && (p1.s == p2.s) && (p1.d == p2.d))
        return true
    end

    return false
end

function ==(p1::PQContractPremiumLatticePoint, p2::PQContractPremiumLatticePoint)::Bool

    if ((p1.i == p2.i) && (p1.j == p2.j) && (p1.s == p2.s) && (p1.d == p2.d))
        return true
    end

    return false
end

# == PUBLIC BELOW HERE ================================================================================================ #
function tabu(world::CRRJITContractPremiumLatticeModel, base::PQContractPremiumLatticePoint;
    fitness::Function, neighborhood::Function)::Tuple{Float64,PQContractPremiumLatticePoint}

    # initialize -
    max_number_of_iterations = 1000
    iteration_counter = 0
    length_of_simple_memory = 80
    should_stop_search = false
    current_point = base
    best_point = base
    best_fitness = fitness(world, best_point)

    # tabu list -
    tabu_list = Array{PQContractPremiumLatticePoint,1}()
    push!(tabu_list, best_point)

    # main loop -
    while (should_stop_search == false)

        # what is the current fitness?
        current_fitness = fitness(world, current_point)

        # generate a neighborhood of the current solution -> check: any good points there?
        list_of_candidate_points = neighborhood(current_point)
        for candidate_point ∈ list_of_candidate_points

            # is this candidate_point in the tabu_list?
            if (in(candidate_point, tabu_list) == false && fitness(world, candidate_point) > current_fitness)
                current_point = candidate_point
            end
        end

        # ok - is the current_point better than the best that we have?
        current_fitness = fitness(world, current_point)
        if (current_fitness > best_fitness)

            # update the best fitness -
            best_fitness = current_fitness

            # update -
            best_point = current_point

            # add this point to the tabu_list -
            push!(tabu_list, best_point)
        end

        # prune the tabu set - remove a random element -
        if (length(tabu_list) >= length_of_simple_memory)
            deleteat!(tabu_list,1) # remove a random element -
        end

        # update the iteration_counter -
        iteration_counter += 1
        if (iteration_counter > max_number_of_iterations)
            should_stop_search = true
        end
    end

    # return -
    return (best_fitness, best_point)
end
# == PUBLIC ABOVE HERE ================================================================================================ #


text:
using DFWannier
using Test
using LinearAlgebra
const DFW = DFWannier

assetfile(f) = joinpath("assets", f)

@time @testset "linalg" begin include("linalg.jl") end
@time @testset "hami_calcs" begin include("hami_calcs.jl") end
# @time @testset "berry" begin include("berry.jl") end
@time @testset "exchanges" begin include("exchanges.jl") end
@time @testset "wan_calcs" begin include("wan_calcs.jl") end


text:
"""
    σ()
σ, the scattering factor of a component.
See [here](https://archimed-platform.github.io/archimed-phys-user-doc/3-inputs/5-models/2-models_list/#translucent)
for more details
"""
Base.@kwdef struct σ{T} <: OpticalProperties
    PAR::T = 0.15
    NIR::T = 0.9
end

"""
Translucent model for light interception, see [here](https://archimed-platform.github.io/archimed-phys-user-doc//3-inputs/5-models/2-models_list/).
"""
Base.@kwdef struct Translucent{T} <: AbstractInterceptionModel
    transparency::T = 0.0
    optical_properties::σ = σ()
end


text:
# Invertible network based hyperbolic layers (Lensink et. al, 2019)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: February 2020

export NetworkHyperbolic, NetworkHyperbolic3D

"""
    H = NetworkHyperbolic(n_in, architecture; k=3, s=1, p=1, logdet=true, α=1f0)
    H = NetworkHyperbolic(n_in, architecture; k=3, s=1, p=1, logdet=true, α=1f0, ndims=2)
    H = NetworkHyperbolic3D(n_in, architecture; k=3, s=1, p=1, logdet=true, α=1f0)

 Create an invertible network based on hyperbolic layers. The network architecture is specified by a tuple
 of the form ((action_1, n_hidden_1), (action_2, n_hidden_2), ... ). Each inner tuple corresonds to an additional layer.
 The first inner tuple argument specifies whether the respective layer increases the number of channels (set to `1`),
 decreases it (set to `-1`) or leaves it constant (set to `0`).  The second argument specifies the number of hidden
 units for that layer.

 *Input*:

 - `n_in`: number of channels of input tensor.

 - `n_hidden`: number of hidden units in residual blocks

 - `architecture`: Tuple of tuples specifying the network architecture; ((action_1, n_hidden_1), (action_2, n_hidden_2))

 - `k`, `s`, `p`: Kernel size, stride and padding of convolutional kernels

 - `logdet`: Bool to indicate whether to return the logdet

 - `α`: Step size in hyperbolic network. Defaults to `1`

 - `ndims`: Number of dimension

 *Output*:

 - `H`: invertible hyperbolic network.

 *Usage:*

 - Forward mode: `Y_prev, Y_curr, logdet = H.forward(X_prev, X_curr)`

 - Inverse mode: `X_curr, X_new = H.inverse(Y_curr, Y_new)`

 - Backward mode: `ΔX_curr, ΔX_new, X_curr, X_new = H.backward(ΔY_curr, ΔY_new, Y_curr, Y_new)`

 *Trainable parameters:*

 - None in `H` itself

 - Trainable parameters in the hyperbolic layers `H.HL[j]`.

 See also: [`CouplingLayer!`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct NetworkHyperbolic <: InvertibleNetwork
    HL::AbstractArray{HyperbolicLayer, 1}
    logdet::Bool
end

@Flux.functor NetworkHyperbolic

# Constructor 2D
function NetworkHyperbolic(n_in::Int64, architecture::NTuple;
                           k=3, s=1, p=1, logdet=true, α=1f0, ndims=2)#, affine_layer=false)

    depth = length(architecture)
    HL = Array{HyperbolicLayer}(undef, depth)

    for j=1:depth

        # Hyperbolic layer at level j
        HL[j] = HyperbolicLayer(n_in, k, s, p; action=architecture[j][1], n_hidden=architecture[j][2], α=α, ndims=ndims)

        # adjust dimensions
        if architecture[j][1] == 1
            n_in = Int(n_in/2^ndims)
        elseif architecture[j][1] == -1
            n_in = Int(n_in*2^ndims)
        end
    end

    return NetworkHyperbolic(HL, logdet)
end

NetworkHyperbolic3D(args...; kw...) = NetworkHyperbolic(args...; kw..., ndims=3)

# Forward pass
function forward(X_prev, X_curr, H::NetworkHyperbolic)
    for j=1:length(H.HL)
        X_prev, X_curr = H.HL[j].forward(X_prev, X_curr)
    end
    return X_prev, X_curr, 1f0  # logdet is always 1
end

# Inverse pass
function inverse(Y_curr, Y_new, H::NetworkHyperbolic)
    for j=length(H.HL):-1:1
        Y_curr, Y_new = H.HL[j].inverse(Y_curr, Y_new)
    end
    return Y_curr, Y_new
end

# Backward pass
function backward(ΔY_curr, ΔY_new, Y_curr, Y_new, H::NetworkHyperbolic; set_grad::Bool=true)
    ~set_grad && (Δθ = Array{Parameter, 1}(undef, 0))
    for j=length(H.HL):-1:1
        if set_grad
            ΔY_curr, ΔY_new, Y_curr, Y_new = H.HL[j].backward(ΔY_curr, ΔY_new, Y_curr, Y_new)
        else
            ΔY_curr, ΔY_new, Δθ_HLj, Y_curr, Y_new = H.HL[j].backward(ΔY_curr, ΔY_new, Y_curr, Y_new; set_grad=set_grad)
            Δθ = cat(Δθ_HLj, Δθ; dims=1)
        end
    end
    set_grad ? (return ΔY_curr, ΔY_new, Y_curr, Y_new) : (return ΔY_curr, ΔY_new, Δθ, Y_curr, Y_new)
end

# Jacobian-related utils
function jacobian(ΔX_prev, ΔX_curr, Δθ::Array{Parameter, 1}, X_prev, X_curr, H::NetworkHyperbolic)
    npars_hl = Int64(length(Δθ)/length(H.HL))
    for j=1:length(H.HL)
        Δθj = Δθ[1+(j-1)*npars_hl:j*npars_hl]
        ΔX_prev, ΔX_curr, X_prev, X_curr = H.HL[j].jacobian(ΔX_prev, ΔX_curr, Δθj, X_prev, X_curr)
    end
    return ΔX_prev, ΔX_curr, X_prev, X_curr
end

adjointJacobian(ΔY_curr, ΔY_new, Y_curr, Y_new, H::NetworkHyperbolic) = backward(ΔY_curr, ΔY_new, Y_curr, Y_new, H; set_grad=false)

## Other utils

# Clear gradients
function clear_grad!(H::NetworkHyperbolic)
    depth = length(H.HL)
    for j=1:depth
        clear_grad!(H.HL[j])
    end
end

# Get parameters
function get_params(H::NetworkHyperbolic)
    depth = length(H.HL)
    p = get_params(H.HL[1])
    if depth > 1
        for j=2:depth
            p = cat(p, get_params(H.HL[j]); dims=1)
        end
    end
    return p
end


text:
Base.@kwdef struct CBSNode
    parent::Union{Nothing,CBSNode}
    constraint::Union{Nothing,Tuple{Int,Int,Int}}
    solution::Solution
    cost::Int
    conflict_heuristic::Int
end

function cbs_low_level(node, additional_constraint, mapf::MAPF)
    forbidden_vertices = Reservation()
    a, t, g = additional_constraint
    for v in mapf.conflict_groups[g]
        push!(forbidden_vertices, (t, v))
    end
    while true
        if isnothing(node.constraint)
            break
        else
            b, t, g = node.constraint
            if b == a
                for v in mapf.conflict_groups[g]
                    push!(forbidden_vertices, (t, v))
                end
            end
            node = node.parent
        end
    end
    graph, edge_weights = mapf.graph, mapf.edge_weights
    s, d, t0 = mapf.sources[a], mapf.destinations[a], mapf.starting_times[a]
    dist = mapf.distances_to_destinations[d]
    heuristic(v) = dist[v]
    path = temporal_astar(
        graph,
        s,
        d,
        t0;
        edge_weights=edge_weights,
        heuristic=heuristic,
        forbidden_vertices=forbidden_vertices,
    )
    return path
end

function conflict_based_search(mapf::MAPF; greedy=true)
    root_solution = independent_astar(mapf)
    root_cost = flowtime(root_solution, mapf)
    root_conflict_heuristic = nb_conflicting_pairs(root_solution, mapf)
    root = CBSNode(;
        parent=nothing,
        constraint=nothing,
        solution=root_solution,
        cost=root_cost,
        conflict_heuristic=root_conflict_heuristic,
    )

    if greedy
        open_queue = VectorPriorityQueue{CBSNode,Tuple{Int,Int}}()
        enqueue!(open_queue, root, (root.conflict_heuristic, root.cost))
    else
        open_queue = VectorPriorityQueue{CBSNode,Int}()
        enqueue!(open_queue, root, root.cost)
    end

    prog = ProgressUnknown("Nodes explored in the conflict tree:")
    while !isempty(open_queue)
        ProgressMeter.next!(prog)

        node = dequeue!(open_queue)
        has_empty_paths(node.solution) && continue
        conflict = find_conflict(node.solution, mapf)

        if isnothing(conflict)
            return node.solution
        else
            (; a, b, t, g) = conflict

            constraint_a = (a, t, g)
            solution_a = copy(node.solution)
            solution_a[a] = cbs_low_level(node, constraint_a, mapf)
            cost_a = flowtime(solution_a, mapf)
            conflict_heuristic_a = nb_conflicting_pairs(solution_a, mapf)
            node_a = CBSNode(;
                parent=node,
                constraint=(a, t, g),
                solution=solution_a,
                cost=cost_a,
                conflict_heuristic=conflict_heuristic_a,
            )

            constraint_b = (b, t, g)
            solution_b = copy(node.solution)
            solution_b[b] = cbs_low_level(node, constraint_b, mapf)
            cost_b = flowtime(solution_b, mapf)
            conflict_heuristic_b = nb_conflicting_pairs(solution_b, mapf)
            node_b = CBSNode(;
                parent=node,
                constraint=(b, t, g),
                solution=solution_b,
                cost=cost_b,
                conflict_heuristic=conflict_heuristic_b,
            )

            if greedy
                enqueue!(open_queue, node_a, (node_a.conflict_heuristic, node_a.cost))
                enqueue!(open_queue, node_b, (node_b.conflict_heuristic, node_b.cost))
            else
                enqueue!(open_queue, node_a, node_a.cost)
                enqueue!(open_queue, node_b, node_b.cost)
            end
        end
    end
    return nothing
end


text:
using Plots, Printf

"""
    shallow_water_1D(;nx=1024, do_visu=false)

1D shallow water equations solver for an instantaneous dam break.
The Lax-Friedrichs Method was applied to the continuity equation.
Geometry (length of 40 meters) and initial conditions (half of domain
have initial water level of 5 meters, other half is dry) match
BASEMENT version 2.8 Test Case H_1 "Dam break in a closed channel

# Arguments
    - `nx`: number of discretised cells.
    - `do_visu`: if true, each physical time step will be ploted.

# Return values
    - `H`: The solution arrray (Water surface height in m).
"""
@views function shallow_water_1D(;
    # Numerics
    nx = 1024,
    # Visualisation
    do_visu = false)
    # Physics
    Lx = 40.0
    g = 9.81
    u_max = 6 #from review of results with very small timesteps
    ttot = 20.0
    # Numerics
    nout = 100
    # Derived numerics
    dx = Lx / nx
    dt = dx / u_max / 10 #stability condition of Lax-Friedrichs Method, further divided by 10 for stability of initially discontinuous wave front
    #println(dt)
    nt = cld(ttot, dt)
    xc = LinRange(dx / 2, Lx - dx / 2, nx)
    # Array initialisation
    H = zeros(Float64, nx)
    H[1:round(Int64(nx / 2))] = 5.0 .+ zeros(Float64, round(Int64(nx / 2)))
    u = zeros(Float64, nx)
    dHdx = zeros(Float64, nx - 1)
    dudx = zeros(Float64, nx - 1)
    dudt = zeros(Float64, nx - 1)
    t_tic = 0.0
    niter = 0
    # Time loop
    for it = 0:nt
        if (it == 11)
            t_tic = Base.time()
            niter = 0
        end
        dHdx .= diff(H) ./ dx
        dudx .= diff(u) ./ dx
        dudt .= -1 / 2 .* (u[1:end-1] .+ u[2:end]) .* dudx .- (g .* dHdx) #momentum
        H[2:end-1] .= 1 / 2 .* (H[3:end] .+ H[1:end-2]) .- dt / 2 / dx .* ((H[3:end] .* u[3:end]) .- (H[1:end-2] .* u[1:end-2])) #continuity with Lax-Friedrichs Method
        H[1] = H[2]
        H[end] = H[end-1]
        u[2:end-1] .= u[2:end-1] .+ 1 / 2 .* (dudt[1:end-1] .+ dudt[2:end]) .* dt
        u[1] = -u[2]
        u[end] = -u[end-1]

        niter += 1
        if do_visu && (it % nout == 0)
            p1 = plot(xc, H, xlims = (xc[1], xc[end]), ylims = (0, 10),
                xlabel = "Lx (m)", ylabel = "water surface elevation (m)", label = "h",
                title = "time = $(round(it*dt, sigdigits=3)) s, stability: $(round(maximum(abs.(u))*dt/dx,sigdigits=3)) /1",
                linewidth = :1.0, framestyle = :box)
            plot!(xc, u, label = "u", linewidth = :1.0)
            display(p1)
        end
    end
    t_toc = Base.time() - t_tic
    A_eff = (2 * 2) / 1e9 * nx * sizeof(Float64)  # Effective main memory access per iteration [GB]
    t_it = t_toc / niter                      # Execution time per iteration [s]
    T_eff = A_eff / t_it                       # Effective memory throughput [GB/s]
    @printf("Time = %1.3f sec, T_eff = %1.2f GB/s (niter = %d)\n", t_toc, round(T_eff, sigdigits = 3), niter)
    return H
end

# shallow_water_1D(;do_visu=true)


text:
using Graphs
using GraphMakie, CairoMakie
using Chain
using LayeredLayouts

const ϕ = (1 + √5) / 2

macro map(col, steps...)
  quote
    map($col) do item
      @chain(item, $(steps...))
    end
  end
end

# The duration of a sequence of dots and dashes
height(morse) = @chain morse collect map(_) do c
  Dict('.' => 1, '-' => 2)[c]
end sum

# By treating each . as 0 and each - as 01 in base Fibonacci, e.g.
#
#   --.-. = 01010010 = F_8 + F_6 + F_3
#
# we can label all the sequences of the same duration with consecutive
# integers. This function uses base ϕ as a shortcut.
basefib(morse) = @chain morse collect map(_) do c
  Dict('.' => "0", '-' => "01")[c]
end join "0" * _ reverse enumerate map(_) do (i, c)
  c == '1' ? ϕ^(i - 1) : 0
end sum

morseplot(labelfunc = first, table = morsetable) = begin
  g = @chain table length SimpleDiGraph
  force_order = Pair{Int,Int}[]
  colors = Symbol[]

  for (i, l) in @chain table first.() enumerate
    for (j, r) in @chain table first.() enumerate
      if l * "." == r
        add_edge!(g, i => j)
        push!(colors, :steelblue)
      elseif l * "-" == r
        add_edge!(g, i => j)
        push!(colors, :brown)
      elseif endswith(l, ".") &&
             endswith(r, "-") &&
             l[1:end-1] == r[1:end-1]
        add_edge!(g, i => j)
        push!(colors, :darkkhaki)
      elseif endswith(l, "..") && endswith(r, "-") && l[1:end-2] == r[1:end-1]
        push!(force_order, i => j)
      end
    end
  end

  _, ax = graphplot(
    g,
    # curves = false,
    nlabels = labelfunc.(table),
    # size = (720, 720),
    # # method = :buchheim,
    # fontsize = 20,
    nlabels_textsize = 20,
    # nlabels_offset = Point(-0.1, -0.1),
    # nlabels_color = :blue,
    edge_color = colors,
    node_attr = (; marker = ' '),
    figure = (; resolution = (720, 720)),
    layout = g -> begin
      x, y, _ = solve_positions(
        Zarate(),
        g,
        force_order = force_order,
        force_layer = @chain table first.() height.() enumerate map(t -> Pair(t...), _)
      )
      Point.(y, -x)
    end
    # layout = _ -> Point.(
    #   (@map table first (basefib(_) + 1) / (ϕ^height(_) + 1)),
    #   (@map table first height -_)
    # ),
  )
  hidedecorations!(ax)
  hidespines!(ax)
  # ax.aspect = 0.1
  current_figure()
end

function make()
  morseplot(first, morsetablecore)
  save("ditsdahs.svg", Makie.current_scene(), resolution = (860, 600))
  morseplot(last, morsetablecore)
  save("alphabet.svg", Makie.current_scene(), resolution = (860, 600))
  morseplot(last, morsetable)
  save("extended.svg", Makie.current_scene(), resolution = (860, 600))
  nothing
end

morsetable = [
  "" => " ",
  "." => "E",
  "-" => "T",
  ".." => "I",
  "-." => "N",
  ".-" => "A",
  "--" => "M",
  "..." => "S",
  "-.." => "D",
  ".-." => "R",
  "--." => "G",
  "..-" => "U",
  "-.-" => "K",
  ".--" => "W",
  "---" => "O",
  "...." => "H",
  "-..." => "B",
  ".-.." => "L",
  "--.." => "Z",
  "..-." => "F",
  "-.-." => "C",
  ".--." => "P",
  "---." => "Ó",
  "...-" => "V",
  "-..-" => "X",
  ".-.-" => "Ä", # Also newline
  "--.-" => "Q",
  "..--" => "Ü",
  "-.--" => "Y",
  ".---" => "J",
  "----" => "Ĥ",
  "....." => "5",
  "-...." => "6",
  ".-..." => "&",
  "--..." => "7",
  # "..-.." => "***",
  "-.-.." => "Ć",
  # ".--.." => "***",
  "---.." => "8",
  # "...-." => "***",
  "-..-." => "/",
  ".-.-." => "+",
  # "--.-." => "***",
  # "..--." => "***",
  "-.--." => "(",
  # ".---." => "***",
  "----." => "9",
  "....-" => "4",
  "-...-" => "=",
  # ".-..-" => "***",
  # "--..-" => "***",
  # "..-.-" => "***",
  # "-.-.-" => "***",
  # ".--.-" => "***",
  # "---.-" => "***",
  "...--" => "3",
  # "-..--" => "***",
  # ".-.--" => "***",
  # "--.--" => "***",
  "..---" => "2",
  # "-.---" => "***",
  ".----" => "1",
  "-----" => "0",
]

morsetablecore = @chain morsetable filter(_) do x
  isascii(last(x)) && isletter(last(x)[1]) || last(x) == " "
end

;


text:
using Distributed
using JLD2

include("../BioClim/src/required.jl")

## Get & prepare data
@time @everywhere begin
    # Load data from CSV files
    df = CSV.read("/Users/kiristern/Documents/GitHub/SDM/data/pred_prey/Puma_concolor.csv", header=true)
    # Prepare data (select columns, arrange values)
    df = prepare_gbif_data(df)
    # Separate species
    taxa_occ = [df[df.species .== u,:] for u in unique(df.species)]
end

## Load predictions
@load "../SDM/data/predictions.jld2" predictions
pred = predictions[1]
prey1 = predictions[2]
prey2 = predictions[3]

# Combine layers
layers = [pred, prey1, prey2]

# Create function
function expand_layers(layers::Array{SDMLayer{Float64},1})
    # Get minimum coordinates
    min_lon = min(map(x -> x.left, layers)...)
    max_lon = max(map(x -> x.right, layers)...)
    min_lat = min(map(x -> x.bottom, layers)...)
    max_lat = max(map(x -> x.top, layers)...)

    # Get grid size (rounding should prevent problems with last decimal)
    grid_size_lons = 1/round(1/stride(layers[1],1))
    grid_size_lats = 1/round(1/stride(layers[1],2))

    # Get coordinate range of newlayer -> original layers must have same stride
    lons_newlayers = min_lon+grid_size_lons:2*grid_size_lons:max_lon-grid_size_lons
    lats_newlayers = min_lat+grid_size_lats:2*grid_size_lats:max_lat-grid_size_lats

    # Create expanded layers
    newlayers = []
    for i in 1:length(layers)
        # Get position of original layers in new layer
        ## find where the point = 0 (= min value; all others are the absolute value)
        m_lon = findmin(abs.(layers[i].left+grid_size_lons .- lons_newlayers))[2]
        M_lon = findmin(abs.(layers[i].right-grid_size_lons .- lons_newlayers))[2]
        m_lat = findmin(abs.(layers[i].bottom+grid_size_lats .- lats_newlayers))[2]
        M_lat = findmin(abs.(layers[i].top-grid_size_lats .- lats_newlayers))[2]
        # Create new grid
        newgrid = fill(NaN, length(lats_newlayers), length(lons_newlayers))
        # Fill in original values
        newgrid[(m_lat:M_lat), (m_lon:M_lon)] .= layers[i].grid
        # Convert to SDMLayer
        newlayer = SDMLayer(newgrid, min_lon, max_lon, min_lat, max_lat)
        # Export result
        push!(newlayers, newlayer)
    end
    return newlayers
end

@time newlayers = expand_layers(layers)

#find probability of species occuring together
probs = copy(newlayers[1].grid)
for j in 2:length(newlayers)
    for k in 1:length(newlayers[1].grid)
        if !isnan(probs[k]) && !isnan(newlayers[j][k])
            probs[k] = probs[k] .* newlayers[j][k]
        end
    end
end
probs

#check if probability of occurence changed (thus overlap in species occurence)
presence = probs .== newlayers[1].grid
#if no change, change to NaN
probs[presence] .= NaN

filter(!isnan, probs)

# Convert to SDMLayer
test = SDMLayer(probs, newlayers[1].left, newlayers[1].right, newlayers[1].bottom, newlayers[1].top)

plot_test = plotSDM(test)

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
# core.jl

using ImageGeoms
using FillArrays: Trues
using Unitful: m

using Test: @test, @testset, @test_throws, @inferred

@testset "construct" begin
    ig = @inferred ImageGeom()
    @test ig isa ImageGeom

    arg1 = ((2,), (3,), (0,))
    arg2 = ((2,3), (3m,4.), (0,1))
    arg3 = ((2,3,4), (3m,4.,1//2), (0,2//3,2.))

    ig1 = @inferred ImageGeom(arg1..., trues(arg1[1]...))
    ig2 = @inferred ImageGeom(arg2..., trues(arg2[1]...))
    ig3 = @inferred ImageGeom(arg3..., trues(arg3[1]...))
    @test ig3 isa ImageGeom

    ig1 = @inferred ImageGeom(arg1...)
    ig2 = @inferred ImageGeom(arg2...)
    ig3 = @inferred ImageGeom(arg3...)
    @test ig3 isa ImageGeom

    ig3a = @inferred ImageGeom(arg3..., Trues(arg3[1]))
    ig3b = @inferred ImageGeom(dims=arg3[1], deltas=arg3[2], offsets=arg3[3])
    ig3c = @inferred ImageGeom(dims=arg3[1], deltas=arg3[2], offsets=arg3[3],
        mask=Trues(arg3[1]))
    @test ig3a == ig3b == ig3c

    # components
    @test ig3.dims == arg3[1]
    @test ig3.deltas == arg3[2]
    @test all(ig3.offsets .≈ arg3[3])
    @test ig3.mask == Trues(ig3.dims)
end

@testset "methods" begin
    dims = (6,8)
    ig = @inferred ImageGeom(; dims)

    show(isinteractive() ? stdout : devnull, ig)
    show(isinteractive() ? stdout : devnull, MIME("text/plain"), ig)

    for f in (zeros, ones, falses, trues)
        @test f(ig.dims) == (@inferred f(ig))
    end
    @test ndims(ig) == 2
    @test size(ig) == dims
    @test size(ig,1) == dims[1]

    # _zero tests
    @test zero(ImageGeom((2,3), (3.0,4.0), (0,0))) === zero(Float64)
    @test zero(ImageGeom((2,3), (3.0,4.0f0), (0,0))) === zero(Float32)
    @test zero(ImageGeom((2,3), (3m,4.0), (0,0))) === zero(Int32) # alert
    @test zero(ImageGeom((2,3), (3.0m,4.0m), (0,0))) === zero(0.0m)

    ig_down = @inferred downsample(ig, (2,2))
    ig_over = @inferred oversample(ig_down, (2,2))
    @test ig_over == ig

    i3 = @inferred ImageGeom(MaskAllButEdge(); dims=(8,6,4))
    i3 = @inferred ImageGeom(MaskCircle(); dims=(8,6,4))
    i3_down = @inferred downsample(i3, 2)
    i3_over = @inferred oversample(i3_down, 2)
    @test i3_over.dims == i3.dims
    @test i3_over.deltas == i3.deltas
    @test i3_over.offsets == i3.offsets
    @test all(i3_over.mask[i3.mask])

    i3a = @inferred ImageGeom( dims=(2,3,4) )
    i3z = @inferred expand_nz(i3a, 1)
    @test size(i3z,3) == size(i3a,3) + 2

    i3z = @inferred expand_nz(i3, 4)
    @test size(i3z,3) == size(i3,3) + 8

    how = (args...; kw...) -> args[1] # trick to avoid Plots
    @test ImageGeoms.plot(ig, how) == axes(ig)[1]
    @test ImageGeoms.plot(i3, how) == axes(i3)[1]

    @test (@inferred axes(i3)) isa Tuple
    @test axis(i3, 2) isa StepRangeLen
    @test (@inferred grids(i3)) isa Tuple

    @test axisf(i3, 2) isa StepRangeLen
    @test (@inferred axesf(i3)) isa Tuple
    @test (@inferred gridf(i3)) isa Tuple
end

@testset "down" begin
    dims = (6,3)
    ig = @inferred ImageGeom(; dims)
    @test_throws String downsample(ig, (2,4))
end


text:
# Functions for parametric representations
@doc """
### para3D

Calculates the parametric locations for a 3D block Taken from GAETAN's code

**Inputs**

*  `X` : Coordinates, size(l,m,n,ndim)
*  `S` : The u, v, w parametric positions

**Outputs**

*  `u` : Averaged u parameters size(l)
*  `v` : Averaged v parameters size(m)
*  `w` : Averaged w parameters size(n)

REFERENCE: Gaetan parameterizations.f90

"""->

function para3D(X, S)

  # Zeros 3 low end faces (or edges if one plane is specified)
  S[1,:,:,1] = 0.0
  S[:,1,:,2] = 0.0
  S[:,:,1,3] = 0.0

  # Zeros out before working
  u = 0.0
  v = 0.0
  w = 0.0

  # Set up the low-edge lines since they are missed by the following loops over
  # most of the low end faces
  for i = 2:l
    S[i,1,1,1] = S[i-1,1,1,1] + del_i(i,1,1)
  end

  for j = 2:m
    S[1,j,1,2] = S[1,j-1,1,2] + del_j(1,j,1)
  end

  for k = 2:n
    S[1,1,k,3] = S[1,1,k-1,3] + del_k(1,1,k)
  end

  # Set up the rest of the low end face lines beacuse theu are missed by the
  # main loop over most of the volume
  for k = 2:n
    for j = 2:m
      S[1,j,k,2] = S[1,j-1,k,2] + del_j(1,j,k)
      S[1,j,k,3] = S[1,j,k-1,3] + del_k(1,j,k)
    end
  end

  for i = 2:l
    for k = 2:n
      S[i,1,k,1] = S[i-1,1,k,1] + del_i(i,1,k)
      S[i,1,k,3] = S[i,1,k-1,3] + del_k(i,1,k)
    end
  end

  for i = 2:l
    for j = 2:m
      S[i,j,1,1] = S[i-1,j,1,1] + del_i(i,j,1)
      S[i,j,1,2] = S[i,j-1,1,2] + del_j(i,j,1)
    end
  end

  # Traverse the block once for all lines except those within the low-end faces
  for i = 2:l
    for j = 2:m
      for k = 2:l
        S[i,j,k,1] = S[i-1,j,k,1] + del_i(i,j,k)
        S[i,j,k,2] = S[i,j-1,k,2] + del_j(i,j,k)
        S[i,j,k,3] = S[i,j,k-1,3] + del_k(i,j,k)
      end
    end
  end

  # Normalizing requires another pass through the volume.
  # Handle lines of zero length first by inserting uniform distributions. Then
  # the standard normalizations can be applied safely everywhere.

  for j = 1:m
    for k = 1:n
      if S[l,j,k,1] == 0.0 # zero length lines in `i` direction
        for i = 2:l
          S[i,j,k,1] = i - 1
        end
      end
    end
  end

  for i = 1:l
    for k = 1:n
      if S[i,m,k,2] == 0.0 # zero length lines in `j` direction
        for j = 2:m
          S[i,j,k,2] = j-1
        end
      end
    end
  end

  for i = 1:l
    for j = 1:m
      if S[i,j,l,3] == 0.0 # zero length lines in `k` direction
        for k = 2:l
          S[i,j,k,3] = k-1
        end
      end
    end
  end

  # Normalize
  for i = 1:l
    for j = 1:m
      for k = 1:n
        S[i,j,k,1] = S[i,j,k,1] / S[n,j,k,1]
        S[i,j,k,2] = S[i,j,k,2] / S[i,m,k,2]
        S[i,j,k,3] = S[i,j,k,3] / S[i,j,n,3]
      end
    end
  end

  # Finally, precise 1s for the three high end faces
  for j = 1:m
    for k = 1:n
      S[l,j,k,1] = 1.0
    end
  end

  for i = 1:l
    for k = 1:n
      S[i,m,k,2] = 1.0
    end
  end

  for i = 1:l
    for j = 1:m
      S[i,j,l,3] = 1.0
    end
  end

  # Get an average u, v, w
  # average u
  for j = 1:m
    for k = 1:n
      u += S[:,j,k,1]
    end
  end
  u = u/(m*n)

  # Average v
  for i = 1:l
    for k = 1:n
      v += S[i,:,k,2]
    end
  end
  v = v/(l*n)

  # Average w
  for i = 1:l
    for j = 1:m
      w += S[i,j,:,3]
    end
  end
  w = w/(l*m)

  return u, v, w
end

@doc """
### del_i, del_j, del_k

Functions to compute the distance between 2 points along a particular parametric
direction.

**Inputs**

*  `X` : Coordinates
*  `i`, `j`, `k` : Indices in the 3 directions

** Outputs**

*  `deli`, `delj`, `delk` : distances along the 3 directions
"""->

function del_i(X, i, j, k)

  deli = sqrt( (X[i,j,k,1] - X[i-1,j,k,1])^2 + (X[i,j,k,2] - X[i-1,j,k,2])^2 +
               (X[i,j,k,3] - X[i-1,j,k,3])^2 )

  return deli
end

function del_j(X, i, j, k)

  delj = sqrt( (X[i,j,k,1] - X[i,j-1,k,1])^2 + (X[i,j,k,2] - X[i,j-1,k,2])^2 +
               (X[i,j,k,3] - X[i,j-1,k,3])^2 )

  return delj
end

function del_k(X, i, j, k)

  delk = sqrt( (X[i,j,k,1] - X[i,j,k-1,1])^2 + (X[i,j,k,2] - X[i,j,k-1,2])^2 +
               (X[i,j,k,3] - X[i,j,k-1,3])^2 )

  return delk
end