text:
#module Contour

import .SimpleGraphics
#using NMR.Scaling
#using NMR.SimplePlot

export ContourPlot

const GridData = Array{Float64,2}

between(y0::Float64,y1::Float64) = y0/(y0-y1)

function TrInterp!(P::Array{Float64,1},ax,ay,az,bx,by,bz,cx,cy,cz)
  sa=sign(az)
  sb=sign(bz)
  sc=sign(cz)

  if sa==sb==sc return
  else   # contour intersects this triangle

    ic=between(az,bz)
    ia=between(bz,cz)
    ib=between(cz,az)

    # traverse the triangle in the direction A->B->C to
    # detect the end points of the contour.
    # This guarantees that the beginning -> end line
    # is aligned in the same sense of rotation.

    startx = 0.0
    starty = 0.0
    stopx = 0.0
    stopy = 0.0

    if sa==0
      if cz>bz
        (startx,starty) = (ax,ay)
      else
        (stopx,stopy) = (ax,ay)
      end
    end
    if 0<ic<1 # contour intersects A-B
      if az>bz
        (startx,starty)=((1-ic)*ax+ic*bx, (1-ic)*ay+ic*by)
      else
        (stopx,stopy)=((1-ic)*ax+ic*bx, (1-ic)*ay+ic*by)
      end
    end

    if sb==0
      if az>cz
        (startx,starty) = (bx,by)
      else
        (stopx,stopy) = (bx,by)
      end
    end
    if 0<ia<1 # contour intersects B-C
      if bz>cz
        (startx,starty)=((1-ia)*bx+ia*cx, (1-ia)*by+ia*cy)
      else
        (stopx,stopy)=((1-ia)*bx+ia*cx, (1-ia)*by+ia*cy)
      end
    end

    if sc==0
      if bz>az
        (startx,starty) = (cx,cy)
      else
        (stopx,stopy) = (cx,cy)
      end
    end
    if 0<ib<1 # contour intersects C-A
      if cz>az
        (startx,starty)=((1-ib)*cx+ib*ax, (1-ib)*cy+ib*ay)
      else
        (stopx,stopy)=((1-ib)*cx+ib*ax, (1-ib)*cy+ib*ay)
      end
    end

    # @printf("%10.5g%10.5g%10.5g%10.5g\n",startx,starty,stopx,stopy)
    push!(P,startx,starty,stopx,stopy)
  end
end

function triangles(z0::Array{Float64,2},l::Real=0)
   (n,m)=size(z0)
   z = z0-l
   P = Array{Float64,1}([]) # P is an array of line segments defined by four consecutive numbers: X1 Y1 X2 Y2

   for k = 1:n-1, l = 1:m-1
      ax=k
      ay=l
      az=z[ax,ay]
      bx=k
      by=l+1
      bz=z[bx,by]
      cx=k+1
      cy=l
      cz=z[cx,cy]

      TrInterp!(P,ax,ay,az,bx,by,bz,cx,cy,cz)

      ax=k+1
      ay=l+1
      az=z[ax,ay]
      bx=k+1
      by=l
      bz=z[bx,by]
      cx=k
      cy=l+1
      cz=z[cx,cy]

      TrInterp!(P,ax,ay,az,bx,by,bz,cx,cy,cz)

    end
    return P
end

const qtol2 = 1.e-20

function find_start(P::Array{Float64,1},used::Array{Bool,1},x::Float64,y::Float64)  # find segment with start point x,y in P
  nseg::Int64=Int(length(P)/4)
  for l=1:nseg
    if used[l] continue end
    qx::Float64=P[4*l-3]
    qy::Float64=P[4*l-2]
    if (x-qx)^2+(y-qy)^2<qtol2 return l end
  end
  return 0
end

function find_end(P::Array{Float64,1},used::Array{Bool,1},x::Float64,y::Float64)  # find segment with end point x,y in P
  nseg::Int64=Int(length(P)/4)
  for l=1:nseg
    if used[l] continue end
    qx::Float64=P[4*l-1]
    qy::Float64=P[4*l]
    if (x-qx)^2+(y-qy)^2<qtol2 return l end
  end
  return 0
end

function collect_branch(P::Array{Float64,1},attr=Dict{String,Any}(["fill"=>"blue" "stroke"=>"grey" "stroke-width"=>"0.01"])) # collect segments that share end points
  g=SimpleGraphics.Group()

  # We assume P to contain consecutive coordinates of beginning and
  # end points of segments. Each segment therefore uses 4 places in P
  # line segments are identified by their index number n/4

  n=length(P)
  nseg=Int(n/4)
  used=zeros(Bool,nseg)

  for k=1:nseg
    if used[k] continue end
    pts=Float64[] # points of contour line are accumulated in pts[]
    px=P[4*k-3]  # start point of k-th segment
    py=P[4*k-2]

    first=k

    # find the beginning of the contour. If it is circular, we go
    # through all segments until we re-encounter the present one (l==k)

    l=find_end(P,used,px,py)
    for j=1:nseg
      if l<=0  || l==k break end
      first=l
      (px,py) = (P[4*l-3],P[4*l-2])
      l=find_end(P,used,px,py)
    end

    px=P[4*first-1]
    py=P[4*first]

    # save the first segment in pts[]

    push!(pts,P[4*first-3],P[4*first-2],px,py)
    used[first]=true

    # find all subsequent segments and store them in pts[]. Stop when
    # no more segments found or when back to start (l==first)
    l=find_start(P,used,px,py)
    for j=1:nseg
      if l<=0 || l==first break end
      (px,py) = (P[4*l-1],P[4*l])
      push!(pts,px,py)
      used[l]=true
      l=find_start(P,used,px,py)
    end

    npts=length(pts)

    SimpleGraphics.push!(g,
      SimpleGraphics.Polygon(
        SimpleGraphics.VertexList(transpose(reshape(pts,2,Int(npts/2)))),
        attr))

  end
  return g
end

function ContourPlot(xrange,yrange,pts;
              levels=NiceScale(minimum(pts),maximum(pts),steps=20,padding=0),
              g=SimpleGraphics.Group(),
              frame=true,
              axes=false,
              show=true,
              yscale=NiceScale(minimum(yrange),maximum(yrange)),
              xscale=NiceScale(minimum(xrange),maximum(xrange)),
              Reverse=[false,false],
              style=Dict(["stroke"=>"black",
                        "fill"=>"none",
                        "stroke-width"=>"1"]),
              FrameW=1000,
              FrameH=1000)

    revrs=Reverse
    if frame g=Frame(xscale,yscale,FrameW=FrameW,FrameH=FrameH,Reverse=revrs) end
    if axes  g=Axes(xscale,yscale,FrameW=FrameW,FrameH=FrameH,Reverse=revrs) end

    (xreverse,yreverse)=Reverse

    if xreverse pts = pts[end:-1:1,:] end
    if yreverse pts = pts[:,end:-1:1] end

    (n,m)=size(pts)

    x0=scaled(minimum(xrange),xscale)
    x1=scaled(maximum(xrange),xscale)
    y0=scaled(minimum(yrange),yscale)
    y1=scaled(maximum(yrange),yscale)

    Q=[FrameW*(x1-x0)/(n-1) 0 ; 0 FrameH*(y1-y0)/(m-1)]

    for lvl in levels
        ut=triangles(pts,lvl)
        branch=collect_branch(ut,style)

        # brnch is on a scale of 1:size(pts). we need to scale and shift
        # so that it conforms into the axes

        SimpleGraphics.push!(g,SimpleGraphics.translate(SimpleGraphics.transform(branch,Q),(x0-1/n)*FrameW,(y0-1/m)*FrameH))

    end

    if show Show(g) end
    return(g)
end

#end # module Contour

max_stars_repo_name: marcel-utz/NMR

max_issues_repo_name: marcel-utz/NMR


text:
using Test
using HMatrices
using HMatrices: hilbert_linear_to_cartesian, hilbert_cartesian_to_linear

@testset "Hilbert points" begin
    # check that hilbert curves spans the lattice, and the inverse is correct
    for n in [2^p for p in 0:3]
        lattice = Set((i,j) for i in 0:n-1,j in 0:n-1)
        for d in 0:(n^2-1)
            x,y = hilbert_linear_to_cartesian(n,d)
            pop!(lattice,(x,y))
            @test hilbert_cartesian_to_linear(n,x,y) == d
        end
        @test isempty(lattice)
    end
end

ext: jl

max_stars_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
module PyLogging

using PyCall

export @critical, @error, @warning, @info, @debug
export getLogger

const logging = pyimport("logging")

# 50, 40, ..., 0
const CRITICAL = convert(Int, logging["CRITICAL"])
const ERROR = convert(Int, logging["ERROR"])
const WARNING = convert(Int, logging["WARNING"])
const INFO = convert(Int, logging["INFO"])
const DEBUG = convert(Int, logging["DEBUG"])
const NOTSET = convert(Int, logging["NOTSET"])

immutable Logger
    logger::PyObject
    # The followings are cached so that we don't have to create many pyobjects.
    isEnabledFor::PyObject
    _log::PyObject

    function Logger(logger::PyObject)
        new(
            logger,
            logger["isEnabledFor"],
            logger["_log"]
        )
    end
end

function getLogger()
    logger = pycall(logging["getLogger"], PyObject)
    Logger(logger)
end

function getLogger(name::AbstractString)
    logger = pycall(logging["getLogger"], PyObject, name)
    Logger(logger)
end

function basicConfig(;kwargs...)
    pycall(logging["basicConfig"], PyObject; kwargs...)
end

function _log(logger::Logger, level::Int, msg::AbstractString)
    logger._log(level, msg, [])
end

function critical(logger::Logger, msg::AbstractString)
    if logger.isEnabledFor(CRITICAL)
        _log(logger, CRITICAL, msg)
    end
end

function error(logger::Logger, msg::AbstractString)
    if logger.isEnabledFor(ERROR)
        _log(logger, ERROR, msg)
    end
end

function warning(logger::Logger, msg::AbstractString)
    if logger.isEnabledFor(WARNING)
        _log(logger, WARNING, msg)
    end
end

function info(logger::Logger, msg::AbstractString)
    if logger.isEnabledFor(INFO)
        _log(logger, INFO, msg)
    end
end

function debug(logger::Logger, msg::AbstractString)
    if logger.isEnabledFor(DEBUG)
        _log(logger, DEBUG, msg)
    end
end

macro critical(logger, msg)
    quote
        x = $(esc(logger))
        if x.isEnabledFor(CRITICAL)
            _log(x, CRITICAL, $(esc(msg)))
        end
    end
end

macro error(logger, msg)
    quote
        x = $(esc(logger))
        if x.isEnabledFor(ERROR)
            _log(x, ERROR, $(esc(msg)))
        end
    end
end

macro warning(logger, msg)
    quote
        x = $(esc(logger))
        if x.isEnabledFor(WARNING)
            _log(x, WARNING, $(esc(msg)))
        end
    end
end

macro info(logger, msg)
    quote
        x = $(esc(logger))
        if x.isEnabledFor(INFO)
            _log(x, INFO, $(esc(msg)))
        end
    end
end

macro debug(logger, msg)
    quote
        x = $(esc(logger))
        if x.isEnabledFor(DEBUG)
            _log(x, DEBUG, $(esc(msg)))
        end
    end
end

end


text:
@vbtime "using TestRepo" begin
    using TestRepo
end

@vbtime "Heavy computation" begin
    TestRepo.heavy_computation()
end


text:
export WhithamGreenNaghdiGPU
using CUDAdrv, CuArrays, CuArrays.CUFFT,LinearAlgebra

"""
    WhithamGreenNaghdiGPU(param;kwargs)

    GPU version of [`WhithamGreenNaghdiGPU`](@ref)

"""
struct WhithamGreenNaghdiGPU <: AbstractModel

    label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	mapfrofull	:: Function
	param	:: NamedTuple
	kwargs  :: NamedTuple

    function WhithamGreenNaghdiGPU(
        param::NamedTuple;
        SGN = false,
        ktol = 0,
        dealias = 0,
		verbose = true
    )

        @show CUDAdrv.name(CuDevice(0))

		if SGN == true
			label = string("Serre-Green-Naghdi")
		else
			label = string("Whitham-Green-Naghdi")
		end
		if verbose @info string("model ",label) end

		kwargs = (SGN=SGN,dealias=dealias,ktol=ktol,verbose=verbose)
		μ 	= param.μ
		ϵ 	= param.ϵ
		mesh = Mesh(param)
		param = ( ϵ = ϵ, μ = μ, xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N )

        k = mesh.k
        x = mesh.x
        x₀ = mesh.x[1]

        ∂ₓ = 1im * mesh.k
        F₁ = tanh.(sqrt(μ) * abs.(k)) ./ (sqrt(μ) * abs.(k))
        F₁[1] = 1                 # Differentiation

        if SGN == true
            F₀ = sqrt(μ) * ∂ₓ
        else
            F₀ = 1im * sqrt.(3 * (1 ./ F₁ .- 1)) .* sign.(k)
        end

        K = mesh.kmax * (1 - dealias / (2 + dealias))
        Π⅔ = abs.(mesh.k) .<= K # Dealiasing low-pass filter
		if dealias == 0
			if verbose @info "no dealiasing" end
			Π⅔ 	= ones(size(mesh.k))
		elseif verbose
			@info string("dealiasing : spectral scheme for power ", dealias + 1," nonlinearity ")
		end
        if verbose @info "LU decomposition" end

        FFT = exp.(-1im * k * (x .- x₀)')
        IFFT = exp.(1im * k * (x .- x₀)') / length(x)
        IFFTF₀ = IFFT * Diagonal(F₀)
        M₀ = IFFT * Diagonal(F₀) * FFT
        Id = Diagonal(ones(size(x)))

        function f!( U::Array{ComplexF64,2} )

            d_fftη = CuArray(U[:, 1])
            d_fftu = CuArray(U[:, 2])
            d_fftv = copy(d_fftu)
            d_h = similar(d_fftu)
            d_u = similar(d_fftu)
            d_hdu = similar(d_fftu)
            d_mo = CuArray(M₀)
            d_pi = CuArray(Π⅔)
            d_id = CuArray(Id)
            d_fo = CuArray(F₀)
            d_L = CuArray(FFT)
            d_dx = CuArray(∂ₓ)
            d_fft = CuArray(FFT)
            d_ifft0 = CuArray(IFFTF₀)

            fw = CUFFT.plan_fft!(d_fftη)
            bw = CUFFT.plan_ifft!(d_fftv)

            d_fftv .= d_fftu
            d_h .= d_fftη
            bw * d_h
            d_h .*= ϵ
            d_h .+= 1

            d_L .= Diagonal(d_h .* d_h .* d_h) * d_ifft0
            d_L .= d_mo * d_L
            d_L .= Diagonal(1 ./ d_h) * d_L
            d_L .= d_fft * d_L
            d_L .= d_id - 1 / 3 * d_L

            d_L, ipiv = CuArrays.CUSOLVER.getrf!(d_L)
            CuArrays.CUSOLVER.getrs!('N', d_L, ipiv, d_fftu)

            d_u .= d_fftu
            d_hdu .= d_u

            bw * d_u

            d_hdu .*= d_fo
            d_hdu .*= d_pi
            bw * d_hdu
            d_hdu .*= d_h
            d_hdu .*= d_hdu
            d_hdu .*= 1 / 2

            bw * d_fftv
            d_fftv .*= d_u
            d_fftv .-= 1 / 2 .* d_u .* d_u
            d_fftv .-= d_hdu

            fw * d_fftv
            d_fftv .*= ϵ
            d_fftv .+= d_fftη
            d_fftv .*= d_pi
            d_fftv .*= d_dx

            bw * d_fftη
            d_fftη .*= d_u
            fw * d_fftη
            d_fftη .*= ϵ
            d_fftη .+= d_fftu
            d_fftη .*= d_pi
            d_fftη .*= d_dx

            d_fftη .*= -1
            d_fftv .*= -1

            U[:, 1] .= Array(d_fftη)
            U[:, 2] .= Array(d_fftv)

        	U[abs.(U) .< ktol] .= 0.0

        end

        function mapto(data::InitialData)

            U = [Π⅔ .* fft(data.η(x)) Π⅔ .* fft(data.v(x))]
            U[abs.(U).<ktol] .= 0
            return U

        end

        function mapfro(U)
            real(ifft(U[:, 1])), real(ifft(U[:, 2]))
        end

        function mapfrofull(U)
            fftη = U[:, 1]
            h = 1 .+ ϵ * ifft(fftη)
            L = Id - 1 / 3 * FFT * Diagonal(1 ./ h) * M₀ * Diagonal(h .^ 3) * IFFTF₀
            real(ifft(U[:, 1])), real(ifft(U[:, 2])), real(ifft(L \ U[:, 2]))
        end

        new( label, f!, mapto, mapfro, mapfrofull, param, kwargs )
    end
end


text:
# The number, 197, is called a circular prime because all rotations of the
# digits: 197, 971, and 719, are themselves prime.
#
# There are thirteen such primes below 100: 2, 3, 5, 7, 11, 13, 17, 31, 37, 71,
# 73, 79, and 97.
#
# How many circular primes are there below one million?

using ProjectEulerSolutions

# Helper function to generate an array of circular shifts given an input
# integer.
function circularnumbers(n::Integer)::Array{Integer}
    digs = digits(n)
    numdigits = length(digs)
    multipliers = 10 .^ (0:numdigits-1)

    circular_nums = Integer[]
    for i in 1:numdigits-1
        push!(circular_nums, sum(multipliers .* circshift(digs, i)))
    end
    return circular_nums
end

# Use Sieve of Eratosthenes to get primes, store in set, then cycle through
# them creating circular shifts and checking whether all shifts are prime.  If
# so, perform some set magic to move those primes into the verified set and
# remove from the potential circular primes set.
function p035solution(n::Integer=100)::Integer
    primes = Set{Integer}(sieve_eratosthenes(n))
    verified = Set{Integer}()

    for prime in primes
        if ndigits(prime) == 1
            push!(verified, prime)
        else
            if mapreduce(x-> x % 2, *, digits(prime)) == 0
                # If any digit is even, ignore number
                continue
            end

            circnums = circularnumbers(prime)
            if mapreduce(prime -> in(prime, primes), &, circnums)
                push!(verified, prime)
                union!(verified, circnums)
                delete!(primes, prime)
                setdiff!(primes, circnums)
            end
        end
    end

    return length(verified)
end

p035 = Problems.Problem(p035solution)

Problems.benchmark(p035, 1_000_000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct Entity <: AbstractEntity
    iid::String
    type::EntityType
end

function Entity(t::Proto.Thing)
    iid = bytes2hex(t.iid)
    isempty(iid) && throw(TypeDBClientException(CONCEPT_MISSING_IID))
    return Entity(iid, instantiate(t._type))
end


text:
"""
	`SaveAtEvent(positions::Tuple)`

This event implements the detection of when the parameter values, used during continuation, equals one of the values in `positions`. This state is then saved in the branch.

For example, you can use it like `continuation(args...; event = SaveAtEvent((1., 2., -3.)))`
"""
SaveAtEvent(positions::Tuple) = ContinuousEvent(length(positions), (it, state) -> map(x -> x - getp(state), positions), ntuple(x->"save-$x",length(positions)))
####################################################################################################
# detection of Fold bifurcation, should be based on BorderedPred
"""
	`FoldDetectEvent`

This event implements the detection of Fold points based on the p-component of the tangent vector to the continuation curve. It is designed to work with the predictor `BorderedPred()` that you pass to `continuation` with the keyword argument `tangentAlgo`.
"""
FoldDetectEvent = ContinuousEvent(1, (it, state) -> state.tau.p, ("fold",))
####################################################################################################
# detection of codim 1 bifurcation
struct BifEvent{Tcb} <: AbstractDiscreteEvent
	nb::Int
	condition::Tcb
end

computeEigenElements(::BifEvent) = true
@inline length(eve::BifEvent) = eve.nb
@inline hasCustomLabels(::BifEvent) = true

function detectBifurcationEVE(iter, state)
	# Note that the computation of eigen-elements should have occured before events are called
	# state should be thus up to date at this stage
	@assert state.n_unstable[1] >=0 "Issue with `detectBifurcationEVE`. Please open an issue on https://github.com/rveltz/BifurcationKit.jl/issues."
	# put the max because n_unstable is initialized at -1 at the beginning of the continuation
	return max(0, state.n_unstable[1])
end

"""
	`BifDetectEvent`

This event implements the detection of bifurcations points along a continuation curve. The detection is based on monitoring  the number of unstable eigenvalues. More details are given at [Detection of bifurcation points](@ref).
"""
BifDetectEvent = BifEvent(1, detectBifurcationEVE)

function getEventType(event::BifEvent, iter::AbstractContinuationIterable, state, verbosity, status::Symbol, interval::Tuple{T, T}, ind = :; typeE = :user) where T
	return getBifurcationType(iter.contParams, state, iter.normC, iter.printSolution, verbosity, status, interval)
end


text:
import RCall  # slow

"""
Auto-Correlation function by RCall.
"""
function acf(X::AbstractVector{T}, lagmax::Int) where {T<:Real}
    res = RCall.rcopy(RCall.rcall(:acf, X, lagmax, plot=false, na_action=:na_pass))
    return res[:acf][2:end]
end

"""
Partial Auto-Correlation function by RCall.
"""
function pacf(X::AbstractVector{T}, lagmax::Int) where {T<:Real}
    res = RCall.rcopy(RCall.rcall(:pacf, X, lagmax, plot=false, na_action=:na_pass))
    return res[:acf][:]
end

"""
Auto-Correlation function of increment process.
"""
function acf_incr(X::AbstractVector{T}, dlags::Union{Int, AbstractVector{Int}}, lagmax::Int; method::Symbol=:acf) where {T<:Real}
    # for single value of dlag: convert to a list
    if typeof(dlags) <: Integer
        dlags = [dlags]
    end

    A = if method==:acf
        [acf((X[l+1:end]-X[1:end-l])[100:end-100], lagmax) for l in dlags]
    elseif method==:pacf
        [pacf((X[l+1:end]-X[1:end-l])[100:end-100], lagmax) for l in dlags]
    else
        error("Invalid method: $(method)")
    end

    # A = []  # output of acf
    # for l in dlags
    #     dX = X[l+1:end]-X[1:end-l]
    #     if method==:acf
    #         push!(A, acf(dX, lagmax))
    #         # push!(A, [cor(dX[t+1:end], dX[1:end-t]) for t in tidx])
    #     elseif method==:pacf
    #         push!(A, pacf(dX, lagmax))
    #     else
    #         error("Invalid method: $(method)")
    #     end
    # end
    return hcat(A...)
end

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
@test adjacency_matrix(g3)[3,2] == 1
@test adjacency_matrix(g3)[2,4] == 0
@test laplacian_matrix(g3)[3,2] == -1
@test laplacian_matrix(g3)[1,3] == 0
@test laplacian_spectrum(g3)[5] == 3.6180339887498945
@test adjacency_spectrum(g3)[1] == -1.732050807568878
@test laplacian_spectrum(g5)[3] == laplacian_spectrum(g5,:both)[3] == 3.0
@test laplacian_spectrum(g5,:in)[3] == 1.0
@test laplacian_spectrum(g5,:out)[3] == 1.0

# check adjacency matrices with self loops
g = copy(g3)
add_edge!(g,1,1)
@test adjacency_matrix(g)[1,1] == 2

g10 = CompleteGraph(10)
B, em = non_backtracking_matrix(g10)
@test length(em) == 2*ne(g10)
@test size(B) == (2*ne(g10),2*ne(g10))
for i=1:10
    @test sum(B[:,i]) == 8
    @test sum(B[i,:]) == 8
end
@test !issymmetric(B)

v = ones(Float64, ne(g10))
z = zeros(Float64, nv(g10))
n10 = Nonbacktracking(g10)
@test size(n10) == (2*ne(g10), 2*ne(g10))
@test eltype(n10) == Float64
@test !issymmetric(n10)

LightGraphs.contract!(z, n10, v)

zprime = contract(n10, v)
@test z == zprime
@test z == 9*ones(Float64, nv(g10))

@test_approx_eq_eps(adjacency_spectrum(g5)[3],0.311, 0.001)

@test adjacency_matrix(g3) ==
    adjacency_matrix(g3, :out) ==
    adjacency_matrix(g3, :in) ==
    adjacency_matrix(g3, :both)

@test_throws ErrorException adjacency_matrix(g3, :purple)

#that call signature works
inmat   = adjacency_matrix(g5, :in, Int)
outmat  = adjacency_matrix(g5, :out, Int)
bothmat = adjacency_matrix(g5, :both, Int)

#relations that should be true
@test inmat' == outmat
@test all((bothmat - outmat) .>= 0)
@test all((bothmat - inmat)  .>= 0)

#check properties of the undirected laplacian carry over.
for dir in [:in, :out, :both]
    amat = adjacency_matrix(g5, dir, Float64)
    lmat = laplacian_matrix(g5, dir, Float64)
    @test isa(amat, SparseMatrixCSC{Float64, Int64})
    @test isa(lmat, SparseMatrixCSC{Float64, Int64})
    evals = eigvals(full(lmat))
    @test all(evals .>= -1e-15) # positive semidefinite
    @test_approx_eq_eps minimum(evals) 0 1e-13
end

# testing incidence_matrix, first directed graph
@test size(incidence_matrix(g4)) == (5,4)
@test incidence_matrix(g4)[1,1] == -1
@test incidence_matrix(g4)[2,1] == 1
@test incidence_matrix(g4)[3,1] == 0
# now undirected graph
@test size(incidence_matrix(g3)) == (5,4)
@test incidence_matrix(g3)[1,1] == 1
@test incidence_matrix(g3)[2,1] == 1
@test incidence_matrix(g3)[3,1] == 0

# TESTS FOR Nonbacktracking operator.

n = 10; k = 5
pg = CompleteGraph(n)
# ϕ1 = nonbacktrack_embedding(pg, k)'

nbt = Nonbacktracking(pg)
B, emap = non_backtracking_matrix(pg)
Bs = sparse(nbt)
@test sparse(B) == Bs
@test_approx_eq_eps(eigs(nbt, nev=1)[1], eigs(B, nev=1)[1], 1e-5)

# check that matvec works
x = ones(Float64, nbt.m)
y = nbt * x
z = B * x
@test norm(y-z) < 1e-8

#check that matmat works and full(nbt) == B
@test norm(nbt*eye(nbt.m) - B) < 1e-8

#check that matmat works and full(nbt) == B
@test norm(nbt*eye(nbt.m) - B) < 1e-8

#check that we can use the implicit matvec in nonbacktrack_embedding
@test size(y) == size(x)

# END tests for Nonbacktracking

# spectral distance checks
for n=3:10
  polygon = random_regular_graph(n, 2)
  @test isapprox(spectral_distance(polygon, polygon), 0, atol=1e-8)
  @test isapprox(spectral_distance(polygon, polygon, 1), 0, atol=1e-8)
end


text:
#= MIT License

Copyright (c) 2020, 2021, 2022 Uwe Fechner

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE. =#

#= Model of a kite-power system in implicit form: residual = f(y, yd)

This model implements a 3D mass-spring system with reel-out. It uses six tether segments (the number can be
configured in the file data/settings.yaml). The kite is modelled using 4 point masses and 3 aerodynamic
surfaces. The spring constant and the damping decrease with the segment length. The aerodynamic kite forces
are acting on three of the four kite point masses.

Four point kite model, included from KiteModels.jl.

Scientific background: http://arxiv.org/abs/1406.6218 =#

# Array of connections of bridlepoints.
# First point, second point, unstressed length.
const SPRINGS_INPUT = [0.    1.  150.
                       1.    2.   -1. # s1, p7, p8
                       4.    2.   -1. # s2, p10, p8
                       4.    5.   -1. # s3, p10, p11
                       3.    4.   -1. # s4, p9, p10
                       5.    1.   -1. # s5, p11, p7
                       4.    1.   -1. # s6, p10, p7
                       3.    5.   -1. # s7, p9, p11
                       5.    2.   -1. # s8, p11, p8
                       2.    3.   -1.] # s9, p8, p9

# struct, defining the phyical parameters of one spring
@with_kw struct Spring{I, S}
    p1::I = 1         # number of the first point
    p2::I = 2         # number of the second point
    length::S = 1.0   # current unstressed spring length
    c_spring::S = 1.0 # spring constant [N/m]
    damping::S  = 0.1 # damping coefficent [Ns/m]
end

const SP = Spring{Int16, Float64}
const KITE_PARTICLES = 4
const KITE_SPRINGS = 9
const KITE_ANGLE = 3.83 # angle between the kite and the last tether segment due to the mass of the control pod
const PRE_STRESS  = 0.9998   # Multiplier for the initial spring lengths.
const KS = deg2rad(16.565 * 1.064 * 0.875 * 1.033 * 0.9757 * 1.083)  # max steering
const DRAG_CORR = 0.93       # correction of the drag for the 4-point model
function zero(::Type{SP})
    SP(0,0,0,0,0)
end

"""
    mutable struct KPS4{S, T, P, Q, SP} <: AbstractKiteModel

State of the kite power system, using a 4 point kite model. Parameters:
- S: Scalar type, e.g. SimFloat
  In the documentation mentioned as Any, but when used in this module it is always SimFloat and not Any.
- T: Vector type, e.g. MVector{3, SimFloat}
- P: number of points of the system, segments+1
- Q: number of springs in the system, P-1
- SP: struct type, describing a spring
Normally a user of this package will not have to access any of the members of this type directly,
use the input and output functions instead.

$(TYPEDFIELDS)
"""
@with_kw mutable struct KPS4{S, T, P, Q, SP} <: AbstractKiteModel
    "Reference to the settings struct"
    set::Settings = se()
    "Reference to the KCU model (Kite Control Unit as implemented in the package KitePodModels"
    kcu::KCU = KCU()
    "Reference to the atmospheric model as implemented in the package AtmosphericModels"
    am::AtmosphericModel = AtmosphericModel()
    "Reference to winch model as implemented in the package WinchModels"
    wm::AbstractWinchModel = AsyncMachine()
    "Iterations, number of calls to the function residual!"
    iter:: Int64 = 0
    "Function for calculation the lift coefficent, using a spline based on the provided value pairs."
    calc_cl = Spline1D(se().alpha_cl, se().cl_list)
    "Function for calculation the drag coefficent, using a spline based on the provided value pairs."
    calc_cd = Spline1D(se().alpha_cd, se().cd_list)
    "wind vector at the height of the kite"
    v_wind::T =           zeros(S, 3)
    "wind vector at reference height"
    v_wind_gnd::T =       zeros(S, 3)
    "wind vector used for the calculation of the tether drag"
    v_wind_tether::T =    zeros(S, 3)
    "apparent wind vector at the kite"
    v_apparent::T =       zeros(S, 3)
    "drag force of kite and bridle; output of calc_aero_forces!"
    drag_force::T =       zeros(S, 3)
    "lift force of the kite; output of calc_aero_forces!"
    lift_force::T =       zeros(S, 3)
    "spring force of the current tether segment, output of calc_particle_forces!"
    spring_force::T =     zeros(S, 3)
    "last winch force"
    last_force::T =       zeros(S, 3)
    "a copy of the residual one (pos,vel) for debugging and unit tests"
    res1::SVector{P, KVec3} = zeros(SVector{P, KVec3})
    "a copy of the residual two (vel,acc) for debugging and unit tests"
    res2::SVector{P, KVec3} = zeros(SVector{P, KVec3})
    "a copy of the actual positions as output for the user"
    pos::SVector{P, KVec3} = zeros(SVector{P, KVec3})
    "velocity vector of the kite"
    vel_kite::T =          zeros(S, 3)
    "unstressed segment length [m]"
    segment_length::S =           0.0
    "lift coefficient of the kite, depending on the angle of attack"
    param_cl::S =         0.2
    "drag coefficient of the kite, depending on the angle of attack"
    param_cd::S =         1.0
    "azimuth angle in radian; inital value is zero"
    psi::S =              zero(S)
    "depower angle [deg]"
    alpha_depower::S =     0.0
    "relative start time of the current time interval"
    t_0::S =               0.0
    "reel out speed of the winch"
    v_reel_out::S =        0.0
    "reel out speed at the last time step"
    last_v_reel_out::S =   0.0
    "unstretched tether length"
    l_tether::S =          0.0
    "air density at the height of the kite"
    rho::S =               0.0
    "actual relative depower setting,  must be between    0 .. 1.0"
    depower::S =           0.0
    "actual relative steering setting, must be between -1.0 .. 1.0"
    steering::S =          0.0
    "multiplier for the stiffniss of tether and bridle"
    stiffness_factor::S =  1.0
    "initial masses of the point masses"
    initial_masses::MVector{P, S} = ones(P)
    "current masses, depending on the total tether length"
    masses::MVector{P, S}         = zeros(P)
    "vector of the springs, defined as struct"
    springs::MVector{Q, SP}       = zeros(SP, Q)
    "vector of the forces, acting on the particles"
    forces::SVector{P, KVec3} = zeros(SVector{P, KVec3})
    "synchronous speed of the motor/ generator"
    sync_speed::S =        0.0
    "x vector of kite reference frame"
    x::T =                 zeros(S, 3)
    "y vector of kite reference frame"
    y::T =                 zeros(S, 3)
    "z vector of kite reference frame"
    z::T =                 zeros(S, 3)
end

"""
    clear!(s::KPS4)

Initialize the kite power model.
"""
function clear!(s::KPS4)
    s.t_0 = 0.0                              # relative start time of the current time interval
    s.v_reel_out = 0.0
    s.last_v_reel_out = 0.0
    s.v_wind_gnd    .= [s.set.v_wind, 0, 0]    # wind vector at reference height
    s.v_wind_tether .= [s.set.v_wind, 0, 0]
    s.v_apparent    .= [s.set.v_wind, 0, 0]
    height = sin(deg2rad(s.set.elevation)) * (s.set.l_tether)
    s.v_wind .= s.v_wind_gnd * calc_wind_factor(s.am, height)

    s.l_tether = s.set.l_tether
    s.segment_length = s.l_tether / s.set.segments
    init_masses!(s)
    init_springs!(s)
    for i in 1:se().segments + KiteModels.KITE_PARTICLES + 1
        s.forces[i] .= zeros(3)
    end
    s.drag_force .= [0.0, 0, 0]
    s.lift_force .= [0.0, 0, 0]
    s.rho = s.set.rho_0
    s.calc_cl = Spline1D(s.set.alpha_cl, s.set.cl_list)
    s.calc_cd = Spline1D(s.set.alpha_cd, s.set.cd_list)
    s.kcu.depower = s.set.depower/100.0
    s.kcu.set_depower = s.kcu.depower
    KiteModels.set_depower_steering!(s, get_depower(s.kcu), get_steering(s.kcu))
end

function KPS4(kcu::KCU)
    s = KPS4{SimFloat, KVec3, kcu.set.segments+KITE_PARTICLES+1, kcu.set.segments+KITE_SPRINGS, SP}()
    s.set = kcu.set
    s.kcu = kcu
    s.calc_cl = Spline1D(s.set.alpha_cl, s.set.cl_list)
    s.calc_cd = Spline1D(s.set.alpha_cd, s.set.cd_list)
    clear!(s)
    return s
end

"""
    calc_particle_forces!(s::KPS4, pos1, pos2, vel1, vel2, spring, segments, d_tether, rho, i)

Calculate the drag force of the tether segment, defined by the parameters pos1, pos2, vel1 and vel2
and distribute it equally on the two particles, that are attached to the segment.
The result is stored in the array s.forces.
"""
@inline function calc_particle_forces!(s, pos1, pos2, vel1, vel2, spring, segments, d_tether, rho, i)
    l_0 = spring.length # Unstressed length
    k = spring.c_spring * s.stiffness_factor  # Spring constant
    c = spring.damping                        # Damping coefficient
    segment = pos1 - pos2
    rel_vel = vel1 - vel2
    av_vel = 0.5 * (vel1 + vel2)
    norm1 = norm(segment)
    unit_vector = segment / norm1

    k1 = 0.25 * k # compression stiffness kite segments
    k2 = 0.1 * k  # compression stiffness tether segments
    c1 = 6.0 * c  # damping kite segments
    spring_vel   = unit_vector ⋅ rel_vel
    if (norm1 - l_0) > 0.0
        if i > segments  # kite springs
             s.spring_force .= (k *  (norm1 - l_0) + (c1 * spring_vel)) * unit_vector
        else
             s.spring_force .= (k *  (norm1 - l_0) + (c * spring_vel)) * unit_vector
        end
    elseif i > segments # kite spring
        s.spring_force .= (k1 *  (norm1 - l_0) + (c * spring_vel)) * unit_vector
    else
        s.spring_force .= (k2 *  (norm1 - l_0) + (c * spring_vel)) * unit_vector
    end

    s.v_apparent .= s.v_wind_tether - av_vel
    if s.set.version == 1
        area = norm1 * d_tether
    else
        area = norm1 * s.set.d_line * 0.001
    end
    v_app_perp = s.v_apparent - s.v_apparent ⋅ unit_vector * unit_vector
    half_drag_force = (-0.25 * rho * s.set.cd_tether * norm(v_app_perp) * area) * v_app_perp

    @inbounds s.forces[spring.p1] .+= half_drag_force + s.spring_force
    @inbounds s.forces[spring.p2] .+= half_drag_force - s.spring_force
    if i == 1 s.last_force .= s.forces[spring.p1] end
    nothing
end

"""
    calc_aero_forces!(s::KPS4, pos, vel, rho, alpha_depower, rel_steering)

Calculates the aerodynamic forces acting on the kite particles.

Parameters:
- pos:              vector of the particle positions
- vel:              vector of the particle velocities
- rho:              air density [kg/m^3]
- rel_depower:      value between  0.0 and  1.0
- alpha_depower:    depower angle [degrees]
- rel_steering:     value between -1.0 and +1.0

Updates the vector s.forces of the first parameter.
"""
function calc_aero_forces!(s::KPS4, pos, vel, rho, alpha_depower, rel_steering)
    rel_side_area = s.set.rel_side_area/100.0    # defined in percent
    K = 1 - rel_side_area                        # correction factor for the drag
    # pos_B, pos_C, pos_D: position of the kite particles B, C, and D
    # v_B,   v_C,   v_D:   velocity of the kite particles B, C, and D
    pos_B, pos_C, pos_D = pos[s.set.segments+3], pos[s.set.segments+4], pos[s.set.segments+5]
    v_B,   v_C,   v_D   = vel[s.set.segments+3], vel[s.set.segments+4], vel[s.set.segments+5]
    va_2,  va_3,  va_4  = s.v_wind - v_B, s.v_wind - v_C, s.v_wind - v_D

    pos_centre = 0.5 * (pos_C + pos_D)
    delta = pos_B - pos_centre
    z = -normalize(delta)
    y = normalize(pos_C - pos_D)
    x = y × z
    s.x .= x; s.y .= y; s.z .= z # save the kite reference frame in the state

    va_xz2 = va_2 - (va_2 ⋅ y) * y
    va_xy3 = va_3 - (va_3 ⋅ z) * z
    va_xy4 = va_4 - (va_4 ⋅ z) * z

    alpha_2 = rad2deg(π - acos2(normalize(va_xz2) ⋅ x) - alpha_depower)     + s.set.alpha_zero
    alpha_3 = rad2deg(π - acos2(normalize(va_xy3) ⋅ x) - rel_steering * KS) + s.set.alpha_ztip
    alpha_4 = rad2deg(π - acos2(normalize(va_xy4) ⋅ x) + rel_steering * KS) + s.set.alpha_ztip

    CL2, CD2 = calc_cl(alpha_2), DRAG_CORR * calc_cd(alpha_2)
    CL3, CD3 = calc_cl(alpha_3), DRAG_CORR * calc_cd(alpha_3)
    CL4, CD4 = calc_cl(alpha_4), DRAG_CORR * calc_cd(alpha_4)
    L2 = (-0.5 * rho * (norm(va_xz2))^2 * s.set.area * CL2) * normalize(va_2 × y)
    L3 = (-0.5 * rho * (norm(va_xy3))^2 * s.set.area * rel_side_area * CL3) * normalize(va_3 × z)
    L4 = (-0.5 * rho * (norm(va_xy4))^2 * s.set.area * rel_side_area * CL4) * normalize(z × va_4)
    D2 = (-0.5 * K * rho * norm(va_2) * s.set.area * CD2) * va_2
    D3 = (-0.5 * K * rho * norm(va_3) * s.set.area * rel_side_area * CD3) * va_3
    D4 = (-0.5 * K * rho * norm(va_4) * s.set.area * rel_side_area * CD4) * va_4
    s.lift_force .= L2
    s.drag_force .= D2 + D3 + D4

    s.forces[s.set.segments + 3] .+= (L2 + D2)
    s.forces[s.set.segments + 4] .+= (L3 + D3)
    s.forces[s.set.segments + 5] .+= (L4 + D4)
end

"""
    inner_loop!(s::KPS4, pos, vel, v_wind_gnd, segments, d_tether)

Calculate the forces, acting on all particles.

Output:
- s.forces
- `s.v_wind_tether`
"""
@inline function inner_loop!(s::KPS4, pos, vel, v_wind_gnd, segments, d_tether)
    for i in eachindex(s.springs)
        p1 = s.springs[i].p1  # First point nr.
        p2 = s.springs[i].p2  # Second point nr.
        height = 0.5 * (pos[p1][3] + pos[p2][3])
        rho = calc_rho(s.am, height)
        @assert height > 0

        s.v_wind_tether .= calc_wind_factor(s.am, height) * v_wind_gnd
        calc_particle_forces!(s, pos[p1], pos[p2], vel[p1], vel[p2], s.springs[i], segments, d_tether, rho, i)
    end
    nothing
end

"""
    loop!(s::KPS4, pos, vel, posd, veld)

Calculate the vectors s.res1 and calculate s.res2 using loops
that iterate over all tether segments.
"""
function loop!(s::KPS4, pos, vel, posd, veld)
    L_0      = s.l_tether / s.set.segments
    if s.set.version == 1
        # for compatibility with the python code and paper
        mass_per_meter = 0.011
    else
        mass_per_meter = s.set.rho_tether * π * (s.set.d_tether/2000.0)^2
    end
    s.res1[1] .= pos[1]
    s.res2[1] .= vel[1]
    particles = s.set.segments + KITE_PARTICLES + 1
    for i in 2:particles
        s.res1[i] .= vel[i] - posd[i]
    end
    # Compute the masses and forces
    m_tether_particle = mass_per_meter * s.segment_length
    s.masses[s.set.segments+1] = s.set.kcu_mass + 0.5 * m_tether_particle
    # TODO: check if the next two lines are correct
    damping  = s.set.damping / L_0
    c_spring = s.set.c_spring/L_0
    for i in 1:s.set.segments
        @inbounds s.masses[i] = m_tether_particle
        @inbounds s.springs[i] = SP(s.springs[i].p1, s.springs[i].p2, s.segment_length, c_spring, damping)
    end
    inner_loop!(s, pos, vel, s.v_wind_gnd, s.set.segments, s.set.d_tether/1000.0)
    for i in 2:particles
        s.res2[i] .= veld[i] - (SVector(0, 0, -G_EARTH) - s.forces[i] / s.masses[i])
    end
end

"""
    residual!(res, yd, y::MVector{S, SimFloat}, s::KPS4, time) where S

    N-point tether model, four points for the kite on top:
    Inputs:
    State vector y   = pos1,  pos2, ... , posn,  vel1,  vel2, . .., veln,  length, v_reel_out
    Derivative   yd  = posd1, posd2, ..., posdn, veld1, veld2, ..., veldn, lengthd, v_reel_outd
    Output:
    Residual     res = res1, res2 = vel1-posd1,  ..., veld1-acc1, ...,

    Additional parameters:
    s: Struct with work variables, type KPS4
    S: The dimension of the state vector
The number of the point masses of the model N = S/6, the state of each point
is represented by two 3 element vectors.
"""
function residual!(res, yd, y::MVector{S, SimFloat}, s::KPS4, time) where S
    T = S-2 # T: three times the number of particles excluding the origin
    segments = div(T,6) - KITE_PARTICLES

    # Reset the force vector to zero.
    for i in 1:segments+KITE_PARTICLES+1
        s.forces[i] .= SVector(0.0, 0, 0)
    end
    # extract the data for the winch simulation
    length,  v_reel_out  = y[end-1],  y[end]
    lengthd, v_reel_outd = yd[end-1], yd[end]
    # extract the data of the particles
    y_  = @view y[1:end-2]
    yd_ = @view yd[1:end-2]
    # unpack the vectors y and yd
    part  = reshape(SVector{T}(y_),  Size(3, div(T,6), 2))
    partd = reshape(SVector{T}(yd_), Size(3, div(T,6), 2))
    pos1, vel1 = part[:,:,1], part[:,:,2]
    pos = SVector{div(T,6)+1}(if i==1 SVector(0.0,0,0) else SVector(pos1[:,i-1]) end for i in 1:div(T,6)+1)
    vel = SVector{div(T,6)+1}(if i==1 SVector(0.0,0,0) else SVector(vel1[:,i-1]) end for i in 1:div(T,6)+1)
    posd1, veld1 = partd[:,:,1], partd[:,:,2]
    posd = SVector{div(T,6)+1}(if i==1 SVector(0.0,0,0) else SVector(posd1[:,i-1]) end for i in 1:div(T,6)+1)
    veld = SVector{div(T,6)+1}(if i==1 SVector(0.0,0,0) else SVector(veld1[:,i-1]) end for i in 1:div(T,6)+1)
    @assert ! isnan(pos[2][3])

    # core calculations
    s.l_tether = length
    s.segment_length = length / segments
    calc_aero_forces!(s, pos, vel, s.rho, s.alpha_depower, s.steering)
    loop!(s, pos, vel, posd, veld)

    # winch calculations
    res[end-1] = lengthd - v_reel_out
    res[end] = v_reel_outd - calc_acceleration(s.wm, s.sync_speed, v_reel_out, norm(s.forces[1]), true)

    # copy and flatten result
    for i in 2:div(T,6)+1
        for j in 1:3
            @inbounds res[3*(i-2)+j]              = s.res1[i][j]
            @inbounds res[3*(div(T,6))+3*(i-2)+j] = s.res2[i][j]
        end
    end
    # copy the position vector for easy debugging
    for i in 1:div(T,6)+1
        @inbounds s.pos[i] .= pos[i]
    end
    s.vel_kite .= vel[end-2]
    s.v_reel_out = v_reel_out

    @assert ! isnan(norm(res))
    s.iter += 1

    nothing
end

# =================== getter functions ====================================================

"""
    calc_height(s::KPS4)

Determine the height of the topmost kite particle above ground.
"""
function calc_height(s::KPS4)
    pos_kite(s)[3]
end

"""
    pos_kite(s::KPS4)

Return the position of the kite (top particle).
"""
function pos_kite(s::KPS4)
    s.pos[end-2]
end

"""
    kite_ref_frame(s::KPS4)

Returns a tuple of the x, y, and z vectors of the kite reference frame.
"""
function kite_ref_frame(s::KPS4)
    s.x, s.y, s.z
end

"""
    winch_force(s::KPS4)

Return the absolute value of the force at the winch as calculated during the last timestep.
"""
function winch_force(s::KPS4) norm(s.last_force) end

# ==================== end of getter functions ================================================

function spring_forces(s::KPS4)
    forces = zeros(SimFloat, s.set.segments+KITE_SPRINGS)
    for i in 1:s.set.segments
        forces[i] =  s.springs[i].c_spring * (norm(s.pos[i+1] - s.pos[i]) - s.segment_length) * s.stiffness_factor
        if forces[i] > 4000.0
            println("Tether raptures for segment $i !")
        end
    end
    for i in 1:KITE_SPRINGS
        p1 = s.springs[i+s.set.segments].p1  # First point nr.
        p2 = s.springs[i+s.set.segments].p2  # Second point nr.
        pos1, pos2 = s.pos[p1], s.pos[p2]
        spring = s.springs[i+s.set.segments]
        l_0 = spring.length # Unstressed length
        k = spring.c_spring * s.stiffness_factor       # Spring constant
        segment = pos1 - pos2
        norm1 = norm(segment)
        k1 = 0.25 * k # compression stiffness kite segments
        if (norm1 - l_0) > 0.0
            spring_force = k *  (norm1 - l_0)
        else
            spring_force = k1 *  (norm1 - l_0)
        end
        forces[i+s.set.segments] = spring_force
        if norm(s.spring_force) > 4000.0
            println("Bridle brakes for spring $i !")
        end
    end
    forces
end

"""
    find_steady_state!(s::KPS4; prn=false, delta = 0.0, stiffness_factor=0.035)

Find an initial equilibrium, based on the inital parameters
`l_tether`, elevation and `v_reel_out`.
"""
function find_steady_state!(s::KPS4; prn=false, delta = 0.0, stiffness_factor=0.035)
    s.stiffness_factor = stiffness_factor
    res = zeros(MVector{6*(s.set.segments+KITE_PARTICLES)+2, SimFloat})
    iter = 0

    # helper function for the steady state finder
    function test_initial_condition!(F, x::Vector)
        x1 = copy(x)
        y0, yd0 = init(s, x1)
        residual!(res, yd0, y0, s, 0.0)
        for i in 1:s.set.segments+KITE_PARTICLES-1
            if i != s.set.segments+KITE_PARTICLES-1
                j = i
            else
                j = i + 1
            end
            # copy the x-component of the residual res2 (acceleration)
            F[i]                               = res[1 + 3*(j-1) + 3*(s.set.segments+KITE_PARTICLES)]
            # copy the z-component of the residual res2
            F[i+s.set.segments+KITE_PARTICLES] = res[3 + 3*(j-1) + 3*(s.set.segments+KITE_PARTICLES)]
        end
        # copy the acceleration of point KCU in x direction
        i = s.set.segments+1
        F[end-1]                               = res[1 + 3*(i-1) + 3*(s.set.segments+KITE_PARTICLES)]
        # copy the acceleration of point C in y direction
        i = s.set.segments+3
        F[end]                                 = res[2 + 3*(i-1) + 3*(s.set.segments+KITE_PARTICLES)]
        iter += 1
        return nothing
    end
    if prn println("\nStarted function test_nlsolve...") end
    X00 = zeros(SimFloat, 2*(s.set.segments+KITE_PARTICLES-1)+2)
    results = nlsolve(test_initial_condition!, X00, autoscale=true, xtol=2e-7, ftol=2e-7, iterations=s.set.max_iter)
    if prn println("\nresult: $results") end
    init(s, results.zero)
end


text:
using ComplexMixtures
using PDBTools
using Plots
using LaTeXStrings
using EasyFit

# Plot defaults
function fig()

plot_font = "Computer Modern"
default(
    fontfamily=plot_font,
    linewidth=2.5,
    framestyle=:box,
    label=nothing,
    grid=false,
    palette=:tab10
)
scalefontsizes()
#scalefontsizes(1.3)

# Load system PDB file
system = readPDB("../simulation/equilibrated.pdb")

# Select the atoms corresponding to glycerol and water
dmf = select(system,"resname DMF")
acr = select(system,"resname FACR or resname ACR or resname LACR")

# Set ComplexMixtures selections
solute = Selection(acr,nmols=1)
solvent = Selection(dmf,natomspermol=12)

# Load previously computed data
mddf_dmf_acr = load("./mddf_dmf_acr.json")

l = @layout [ grid(2,1); grid(2,1); c{0.3h} ]
plot(layout=l)

# Plot the MDDF of DMF relative to PolyACR and its corresponding KB integral
x = mddf_dmf_acr.d
y = movavg(mddf_dmf_acr.mddf,n=10).x
plot!(x,y,ylabel="MDDF",subplot=1)
plot!(xlim=(0,20),subplot=1,xticks=:none)

y = movavg(mddf_dmf_acr.kb,n=10).x
plot!(x,y,xlabel=L"\textrm{Distance / \AA}",ylabel=L"\textrm{KB~/~cm^3~mol^{-1}}",subplot=2)
plot!(xlim=(0,20),subplot=2)

# Plot DMF group contributions to the MDDF. We use a vector of tuples here,
# the first element of the tuple is a list of atom names, and the second element
# is the label of the group, to be used in the plot.
groups = [
    (["C","O"],"CO"), # carbonyl
    (["HA"],"HA"),
    (["N"],"N"),
    (["CC","CT","HC1","HC2","HC3","HT1","HT2","HT3"],"Methyl groups")
]
x = mddf_dmf_acr.d
y = movavg(mddf_dmf_acr.mddf,n=10).x
plot!(x,y,label="Total",subplot=3,xticks=:none)
for group in groups
    # Retrieve the contributions of the atoms of this group
    group_contrib = contrib(solvent,mddf_dmf_acr.solvent_atom,group[1])
    y = movavg(group_contrib,n=10).x
    plot!(x,y,label=group[2],subplot=3)
end
plot!(
    xlim=(1,5),
    ylabel="MDDF",
    subplot=3
)

# Plot ACR group contributions to the MDDF
groups = [
    (["CF","HF1","HF2","HF3",
      "CL","HL1","HL2","HL3"],L"\textrm{CH_3}"), # terminal methyles
    (["OE1","CD"],"CO"), # carbonyl
    (["NE2","HE22","HE21"],L"\textrm{NH_2}"), # amine
    (["C","H2","H1","CA","HA"],L"\textrm{CHCH_2}"), # backbone
]
x = mddf_dmf_acr.d
y = movavg(mddf_dmf_acr.mddf,n=10).x
plot!(x,y,label="Total",subplot=4)
for group in groups
    group_contrib = contrib(solute,mddf_dmf_acr.solute_atom,group[1])
    y = movavg(group_contrib,n=10).x
    plot!(x,y,label=group[2],subplot=4)
end
plot!(
    xlim=(1,5),
    xlabel=L"\textrm{Distance / \AA}",
    ylabel="MDDF", subplot=4
)

# 2D plot of group contributions
groups = [
    (["CF","HF1","HF2","HF3"],L"\textrm{CH_3}"),
    (["OE1","CD"],"CO"),
    (["NE2","HE22","HE21"],L"\textrm{NH_2}"),
    (["C","H2","H1","CA","HA"],L"\textrm{CHCH_2}"),
    (["CL","HL1","HL2","HL3"],L"\textrm{CH_3}")
]

# Here we split the polymer in residues, to extract the contribution of
# each chemical group of each polymer mer independently
mers = collect(eachresidue(acr))
group_contribs = zeros(length(mddf_dmf_acr.d),0)
labels = String[]
for (imer, mer) in pairs(mers)
    # The mer.atoms array is just the original array of atoms,
    # and the atoms of this residue are retrieved by choosing the
    # corresponding range (the eachresidue function does not copy)
    mer_atoms = mer.atoms[mer.range]
    for igroup in 1:length(groups)
        group = groups[igroup]
        (imer != 1 && igroup == 1) && continue # only first residue has a first CH3
        (imer != 5 && igroup == 5) && continue # only last residue has a terminal CH3
        # And from the mer_atoms atoms, filter the ones corresponding to this group
        atoms = filter(at -> at.name in group[1], mer_atoms)
        # Retrive the contribution of these atoms to the MDDF
        contribs = movavg(contrib(solute,mddf_dmf_acr.solute_atom,atoms),n=10).x
        # Concatenate the results to build the 2D matrix
        group_contribs = hcat(group_contribs,contribs)
        # Push label to label list
        push!(labels,group[2])
    end
end

# Find the indices of the limits of the map we want
idmin = findfirst( d -> d > 1.5, mddf_dmf_acr.d)
idmax = findfirst( d -> d > 3.2, mddf_dmf_acr.d)

# Plot contour map
contourf!(
    1:length(labels),
    mddf_dmf_acr.d[idmin:idmax],
    group_contribs[idmin:idmax,:],
    color=cgrad(:tempo),linewidth=1,linecolor=:black,
    colorbar=:none,levels=10,
    xlabel="PolyACR group",ylabel=L"\textrm{Distance/\AA}",xrotation=60,
    xticks=(1:length(labels),labels),
    subplot=5,
)

for (y,lab) in [(8.2,"A)"),
                (7.1,"B)"),
                (5.4,"C)"),
                (4.4,"D)"),
                (2.4,"E)")]
    annotate!(
        -2.3, y,
        text(lab,plot_font,12),
        subplot=5,
    )
end

annotate!(1.1, 2.9, text("DMF\ngroups", plot_font, 10, :left), subplot=3)
annotate!(1.1, 2.9, text("PolyACR\ngroups", plot_font, 10, :left), subplot=4)

plot!(
    size=(550,800),
    leftmargin=10Plots.Measures.mm # adjust margin
)

savefig("./polyacr.png")
savefig("./figure5.pdf")

end; fig()


text:
using Revise
Revise.revise()

using Term
import Term: AbstractRenderable
import Term: install_stacktrace

install_stacktrace()

# ------------------------------- inspect types ------------------------------ #
inspect(Panel)
inspect(AbstractRenderable;)
inspect(String;)

# ----------------------------- inspect variables ---------------------------- #
inspect("test")
inspect(1)

# ------------------------------ inspect methods ----------------------------- #
# inspect(print)
inspect(inspect)


text:
"""Julia program to implement Quick sort
Quick Sort is a Divide and Conquer strategy based sorting algoithm, that partitions the array
around a selected pivot element. In the below implemented version of Quick Sort we always took the
last element as the pivot element.

"""
function  partition(arr, low, high)
    # Select the last element as pivot
    pivot = arr[high]
    i = low - 1
    for j in low:high
        # Check if there is any element less than the pivot element
        if(arr[j] < pivot)
            i = i + 1
            # Swap such an element with the to the start position
            arr[i], arr[j] = arr[j], arr[i]
        end
    end
    # Now place the pivot to its original position
    arr[i+1], arr[high] = arr[high], arr[i+1]
    return (i+1)
end

function quick_sort(arr, low, high)
    if (low < high)
        pi = partition(arr, low, high)
        # Recursively call quick_sort function
        quick_sort(arr, low, pi - 1)
        quick_sort(arr, pi + 1, high)
    end
end

print("How many numbers do you want to sort? ")
n = readline()
n = parse(Int, n)
if (n <= 0)
    println("No numbers to sort!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
quick_sort(arr, 1, n)
print("The numbers in sorted order is ")
print(arr)

"""
Time Complexity - O(n^2), where 'n' is the size of the array
Space Complexity - O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE I
How many numbers do you want to sort? 5
Enter the numbers: 4 2 5 1 3
The numbers in sorted order is [1, 2, 3, 4, 5]

"""


text:
@testset "Impute" begin
    ds = KeyedDataset(
        flatten([
            :train => [
                :temp => KeyedArray([1.0 1.1; missing 2.2; 3.0 3.3]; time=1:3, id=[:a, :b]),
                :load => KeyedArray([7.0 7.7; 8.0 missing; 9.0 9.9]; time=1:3, loc=[:x, :y]),
            ],
            :predict => [
                :temp => KeyedArray([1.0 missing; 2.0 2.2; 3.0 3.3]; time=1:3, id=[:a, :b]),
                :load => KeyedArray([7.0 7.7; 8.1 missing; 9.0 9.9]; time=1:3, loc=[:x, :y]),
            ]
        ])...
    );

    @testset "validate" begin
        @test_throws ThresholdError Impute.threshold(ds; limit=0.1)
        r = Impute.threshold(ds; limit=1.0)
        @test isequal(r, ds)
    end

    @testset "declaremissings" begin
        nonmissing = KeyedDataset(
           flatten([
               :train => [
                   :temp => KeyedArray(allowmissing([1.0 1.1; -9999.0 2.2; 3.0 3.3]); time=1:3, id=[:a, :b]),
                   :load => KeyedArray(allowmissing([7.0 7.7; 8.0 NaN; 9.0 9.9]); time=1:3, loc=[:x, :y]),
                ],
                :predict => [
                   :temp => KeyedArray(allowmissing([1.0 -9999.0; 2.0 2.2; 3.0 3.3]); time=1:3, id=[:a, :b]),
                   :load => KeyedArray(allowmissing([7.0 7.7; 8.1 NaN; 9.0 9.9]); time=1:3, loc=[:x, :y]),
                ]
            ])...
        );

        r = Impute.declaremissings(nonmissing; values=(NaN, -9999.0))
        @test isequal(r, ds)
    end

    @testset "impute" begin
        expected = KeyedDataset(
            flatten([
                :train => [
                    :temp => KeyedArray(allowmissing([1.0 1.1; 2.2 2.2; 3.0 3.3]); time=1:3, id=[:a, :b]),
                    :load => KeyedArray(allowmissing([7.0 7.7; 8.0 8.0; 9.0 9.9]); time=1:3, loc=[:x, :y]),
                ],
                :predict => [
                    :temp => KeyedArray(allowmissing([1.0 1.0; 2.0 2.2; 3.0 3.3]); time=1:3, id=[:a, :b]),
                    :load => KeyedArray(allowmissing([7.0 7.7; 8.1 8.1; 9.0 9.9]); time=1:3, loc=[:x, :y]),
                ]
            ])...
        );

        r = Impute.substitute(ds; dims=:time)
        @test isequal(r, expected)

        expected = KeyedDataset(
            flatten([
                :train => [
                    :temp => KeyedArray([1.0 1.1; missing 2.2; 3.0 3.3]; time=1:3, id=[:a, :b]),
                    :load => KeyedArray(allowmissing([7.0 7.7; 8.0 8.8; 9.0 9.9]); time=1:3, loc=[:x, :y]),
                ],
                :predict => [
                    :temp => KeyedArray([1.0 missing; 2.0 2.2; 3.0 3.3]; time=1:3, id=[:a, :b]),
                    :load => KeyedArray(allowmissing([7.0 7.7; 8.1 8.8; 9.0 9.9]); time=1:3, loc=[:x, :y]),
                ]
            ])...
        );

        r = Impute.substitute(ds; dims=:loc)
        @test isequal(r, expected)
    end

    @testset "filter" begin
        # All time axis are constraint, so they must match and removing a row from 1
        # component will remove it from all
        expected = KeyedDataset(
            flatten([
                :train => [
                    :temp => KeyedArray([3.0 3.3]; time=3:3, id=[:a, :b]),
                    :load => KeyedArray([9.0 9.9]; time=3:3, loc=[:x, :y]),
                ],
                :predict => [
                    :temp => KeyedArray([3.0 3.3]; time=3:3, id=[:a, :b]),
                    :load => KeyedArray([9.0 9.9]; time=3:3, loc=[:x, :y]),
                ]
            ])...
        );

        r = Impute.filter(ds; dims=:time)
        @test isequal(r, expected)

        # All time axis are constraint but only :load is checked. :temp is updated to match
        expected = KeyedDataset(
            flatten([
                :train => [
                    :temp => KeyedArray([1.0 1.1; 3.0 3.3]; time=[1, 3], id=[:a, :b]),
                    :load => KeyedArray([7.0 7.7; 9.0 9.9]; time=[1, 3], loc=[:x, :y]),
                ],
                :predict => [
                    :temp => KeyedArray([1.0 missing; 3.0 3.3]; time=[1, 3], id=[:a, :b]),
                    :load => KeyedArray([7.0 7.7; 9.0 9.9]; time=[1, 3], loc=[:x, :y]),
                ]
            ])...
        );
        r = Impute.filter(ds; dims=Pattern(:__, :load, :time))
        @test isequal(r, expected)

        r = Impute.filter(ds; dims=(:__, :load, :time))
        @test isequal(r, expected)

        # Pattern must end in a dimname
        @test_throws ArgumentError Impute.filter(ds; dims=Pattern(:__, :load, :_))

        # Only :load has a shared :loc axis, so we see that that :y location is dropped from
        # both train and predict.
        expected = KeyedDataset(
            flatten([
                :train => [
                    :temp => KeyedArray([1.0 1.1; missing 2.2; 3.0 3.3]; time=1:3, id=[:a, :b]),
                    :load => KeyedArray([7.0, 8.0, 9.0][:, :]; time=1:3, loc=[:x]),
                ],
                :predict => [
                    :temp => KeyedArray([1.0 missing; 2.0 2.2; 3.0 3.3]; time=1:3, id=[:a, :b]),
                    :load => KeyedArray([7.0, 8.1, 9.0][:, :]; time=1:3, loc=[:x]),
                ]
            ])...
        );

        r = Impute.filter(ds; dims=:loc)
        @test isequal(r, expected)

        constrained_ds = KeyedDataset(
            flatten([
                :train => [
                    :temp => KeyedArray([1.0 1.1; missing 2.2; 3.0 3.3]; time=1:3, id=[:a, :b]),
                    :load => KeyedArray([7.0 7.7; 8.0 missing; 9.0 9.9]; time=1:3, loc=[:x, :y]),
                ],
                :predict => [
                    :temp => KeyedArray([1.0 missing; 2.0 2.2; 3.0 3.3]; time=1:3, id=[:a, :b]),
                    :load => KeyedArray([7.0 7.7; 8.1 missing; 9.0 9.9]; time=1:3, loc=[:x, :y]),
                ]
            ])...,
            constraints = Pattern[(:__, :time), (:__, :loc), (:train, :__, :id)]
        );

        # :id for :train is included in the constraints and :predict is unconstrained
        # Both drop only the :ids missing in their own table
        expected = KeyedDataset(
            flatten([
                :train => [
                    :temp => KeyedArray([1.1; 2.2; 3.3][:, :]; time=1:3, id=[:b]),
                    :load => KeyedArray([7.0 7.7; 8.0 missing; 9.0 9.9]; time=1:3, loc=[:x, :y]),
                ],
                :predict => [
                    :temp => KeyedArray([1.0; 2.0; 3.0][:, :]; time=1:3, id=[:a]),
                    :load => KeyedArray([7.0 7.7; 8.1 missing; 9.0 9.9]; time=1:3, loc=[:x, :y]),
                ]
            ])...,
            constraints = Pattern[(:__, :time), (:__, :loc),  (:train, :__, :id)]
        );

        r = Impute.filter(constrained_ds; dims=:id)
        @test isequal(r, expected)

    end
end