text:
function StokesPcCG_pardiso(U, P, KK, MM, GG, Rhs, ifree)

    rtol_Pat = 1e-8
    itmax_Pat = 200
    itmin_Pat = 10
    itnum = 0

    # allocate some stuff
    # fill!(U, 0.0)
    # fill!(P, 0.0)
    nU = length(Rhs)
    z = fill(0.0, nU)
    y = Vector{Float64}(undef, GG.m)
    Sq = Vector{Float64}(undef, GG.n)

    # forces resulting from pressure gradients
    # Rhs += spmv(GG, P)
    Rhs .+= GG*P

    # guess of the velocity field
    ps, A_pardiso = _MKLfactorize(KK, Rhs, ifree)
    _MKLsolve!(U, A_pardiso, ps, Rhs, ifree)

    # initial pressure residual vector
    GGtransp = GG'
    r = -GGtransp * U
    rm0space!(r)
    Prms = mynorm(r) # norm of r_i
    tol = Prms * rtol_Pat

    # get preconditioner
    # MM, pc = _preconditioner(:jacobi, MM)
    MM, pc = _preconditioner(:incomplete_chol, MM)
    # pc = AMGPreconditioner{RugeStuben}(MM)
    d = pc\r # precondition residual

    # Begin of Patera pressure iterations
    # d = _precondition(pc, MM, r) # precondition residual
    # d = _precondition(pc, r) # precondition residual
    q = deepcopy(d)  # define FIRST search direction q
    rlast = similar(r)

    for itPat = 1:itmax_Pat
        itnum +=1

        rm0space!(q)
        rd = mydot(r,d) # numerator for alpha

        #= Perform the S times q multiplication ======================
            S cannot be calculated explicitly, since Kinv cannot be formed
            Sq = S * q = (G' * Kinv * G) * q
            Hence, the muliplicatipon is done in 3 steps:
            (1) y   = G*q
            (2) K z = y
            (3) Sq  = G'*z
        =============================================================#
        # (1) y   = G*q
        # spmv!(y, GG, q)
        y .= GG*q
        # (2) Solve K z = y
        _MKLsolve!(z, A_pardiso, ps, y, ifree)
        # (3) Sq  = G'*z
        # spmv!(Sq, GGtransp, z)
        Sq .= GGtransp*z
        #=============================================================#
        qSq = mydot(q,Sq) # denominator to calculate alpha
        copyto!(rlast, r) # needed for Polak-Ribiere version 1
        α = rd/qSq # steps size in direction q

        # Update solution and residual
        _updatesolution!(P,U,r,q,z,Sq,α)
        # remove nullspace
        rm0space!(r)

        # Check convergence
        if mynorm(r) < tol && itPat>=itmin_Pat
            println("\n", itnum," CG iterations\n")
            break
        end

        # fill!(z, 0.0)

        d = pc\r # precondition residual
        # d  = _precondition(pc, r) # precondition residual
        # Make new search direction q S-orthogonal to all previous q's
        β  = mydot(r-rlast,d)/rd # Polak-Ribiere version 1
        xpy!(q, d, β)

    end

    _MKLrelease!(ps)

    return U,P

end

function StokesPcCG(U,P,KK,MM,GG,Rhs,ifree)

    rtol_Pat = 1e-8
    itmax_Pat = 200
    itmin_Pat = 10
    itnum = 0

    # allocate some stuff
    fill!(U, 0.0)
    fill!(P, 0.0)
    nU = length(Rhs)
    z = fill(0.0, nU)

    # forces resulting from pressure gradients
    # Rhs += spmv(GG, P)
    Rhs += GG*P

    # guess of the velocity field
    U,F = _CholeskyFactorizationSolve(U, KK, Rhs, ifree) # return factorization F to speed up next direct solvers

    # initial pressure residual vector
    GGtransp = sparse(GG')
    r = -GGtransp *U
    rm0space!(r)
    Prms = mynorm(r) # norm of r_i
    tol = Prms * rtol_Pat
    # get preconditioner
    MM, pc = _preconditioner(:jacobi, MM)
    # MM, pc = _preconditioner("lumped", MM)

    # Begin of Patera pressure iterations
    # d = _precondition(pc, MM, r) # precondition residual
    d = _precondition(pc, r) # precondition residual
    q = deepcopy(d)  # define FIRST search direction q
    rlast = similar(r)
    for itPat = 1:itmax_Pat
        itnum +=1

        rm0space!(q)
        rd = mydot(r,d) # numerator for alpha

        #= Perform the S times q multiplication ======================
            S cannot be calculated explicitly, since Kinv cannot be formed
            Sq = S * q = (Gᵀ * Kinv * G) * q
            Hence, the muliplicatipon is done in 3 steps:
            (1) y   = G*q
            (2) K z = y
            (3) Sq  = Gᵀ*z
        =============================================================#
        y = GG * q # (1) y   = G*q
        _CholeskyWithFactorization!(z, F, y, ifree) # (2) Solve K z = y
        Sq = GGtransp*z # (3) Sq  = G'*z
        #=============================================================#
        qSq = mydot(q, Sq) # denominator to calculate alpha
        copyto!(rlast, r) # needed for Polak-Ribiere version 1
        α = rd/qSq # steps size in direction q

        # Update solution and residual
        _updatesolution!(P,U,r,q,z,Sq,α)
        # remove nullspace
        rm0space!(r)

        # Check convergence
        if mynorm(r) < tol && itPat>=itmin_Pat
            println("\n", itnum," CG iterations\n")
            break
        end

        d  = _precondition(pc, r) # precondition residual
        # Make new search direction q S-orthogonal to all previous q's
        β  = mydot(r-rlast, d)/rd # Polak-Ribiere version 1
        xpy!(q, d, β)

    end

    return U,P

end

# Function takes out the mean of vector a
# (thereby removes the constant pressure mode, i.e. the nullspace)
function rm0space!(a)
    m = vmean(a)
    @tturbo for i in eachindex(a)
        a[i] -= m
    end
end

function vmean(a::AbstractArray{T}) where T
    m = zero(T)
    n = length(a)
    @tturbo for i in eachindex(a)
        m += a[i]
    end
    m/n
end

# PRECONDITIONERS ================================================================
@inline function _preconditioner(type, MM)

    MM = MM .+ tril(MM,-1)'

    if type == :diagonal
        #=
            Inverse of the diagonal of the mass matrix scaled by viscosity
        =#
        p = diagonalpreconditioner(MM)

    elseif type == :lumped
        #=
            Inverse of lumped mass matrix scaled by viscosity
        =#
        p  = lumpedpreconditioner(MM)

    elseif type == :jacobi
        #=
            Inverse of lumped mass matrix scaled by viscosity
        =#
        p  = lumpedpreconditioner(MM)

    elseif type == :incomplete_chol
        #=
            Inverse of lumped mass matrix scaled by viscosity
        =#
        p  = lltpreconditioner(MM)
    end
    return MM, p
end # END OF PRECONDITIONERS

struct DiagonalPreconditioner{T}
    p::T
end

function invdiagsp(A::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
    si, sj = A.m, A.n
    @assert si == sj
    I = collect(1:si)
    d = Vector{Tv}(undef, si)
    Threads.@threads for i in I
        @inbounds d[i] = 1/A[i,i]
    end
    sparse(I, I, d)
end

@inline function diagonalpreconditioner(MM)
    return (invdiagsp(MM))
end

@inline function lumpedpreconditioner(MM)
    return dropdims(sum(MM,dims=2), dims=2).^-1
end

@inline function lltpreconditioner(MM)
    return lldl(MM)
end

@inline function _updatesolution!(P,U,r,q,z,Sq,α)
    n1 = length(P)
    n2 = length(U)
    idx1 = 1:n1
    idx2 = (n1+1):n2
    @tturbo for i in idx1
        P[i] += α * q[i]
        U[i] += α * z[i]
        r[i] -= α * Sq[i]
    end
    @tturbo for i in idx2
        U[i] += α * z[i]
    end
end

# PRECONDITIONING: JACOBI ITERATIONS
@inline function _precondition(pc::Array{Float64,1}, MM::SparseMatrixCSC{Float64,Int64} ,r::Array{Float64,1})
    weight = [0.2500, 0.4000, 0.7500]
    d = r .* pc
    dummy = Vector{Float64}(undef, length(d))
    @inbounds for is in 1:3
        # d = d .+ weight[is] * pc .* (r .- MM*d)
        d = d .+ weight[is] * pc .* (r .- mul!(dummy, MM, d))
    end
    d
end

# PRECONDITIONING: lumped
@inline _precondition(pc::Vector,r::Vector) = r .* pc

function _CG!(T,KK,Rhs,ifree)
    # Load stiffness matrix and rhs
    A  = KK[ifree,ifree]
    b  = Vector{Float64}(undef,length(ifree))
    @turbo for  i ∈ 1:length(ifree)
        b[i] = Rhs[ifree[i]]
    end
    # Diagonal preconditioner
    p = DiagonalPreconditioner(A)
    # p = AMGPreconditioner{SmoothedAggregation}(A)
    # Conjugate gradients
    x = cg(A, b, Pl=p)

    # Rebuild solution array
    c   = 0
    @inbounds for i ∈ ifree
        c   +=1
        T[i] = x[c]
    end

end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
"""
Set Neumann boundary conditions
"""
function Neumann(t, xs, Ds, Fs, bc; bc_value = 0)
    bcs = []
    for F in Fs
        for i in eachindex(xs)
            temp = xs[i]
            xs[i] = bc
            push!(bcs,Ds[i](F(t, xs...)) ~ bc_value)
            xs[i] = temp
        end
    end

    return bcs
end

"""
Set Dirichlet boundary conditions

Dirichlet(t, xs, Es, lb, bc_value=2)

where Es must be the same length as xs.
"""
function Dirichlet(t, xs, Ds, Fs, bc; bc_value = 0)
    bcs = []
    for F in Fs
        for i in eachindex(xs)
            temp = xs[i]
            xs[i] = bc
            push!(bcs,F(t, xs...) ~ bc_value)
            xs[i] = temp
        end
    end

    return bcs
end

"""
Set reflective boundary conditions (for f)
"""
function Reflective(t, xs, vs, Fs, bc, a, g)
    bcs = []
    for F in Fs
        for i in eachindex(vs)
            temp = vs[i]
            vs[i] = bc
            lhs_value = F(t, xs..., vs...)
            vs[i] = -bc
            rhs_value = F(t, xs..., vs...)
            push!(bcs,lhs_value ~ a * rhs_value + g)
            vs[i] = temp
        end
    end

    return bcs
end

ext: jl


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# tests to check that code runs without crashing, but not checking that output is expected

# check that different sim configs can be opened, and sim can run
@testset "sim configs" begin
	@assert(isdir("data/cities/small/1/generated"))
	simConfigFolder = "data/cities/small/1/sim_configs"
	for configFilename in readdir(simConfigFolder)
		filename = joinpath(pwd(), simConfigFolder, configFilename)
		sim = initSim(filename);
		simulate!(sim)
		@test true
	end
end

# test the scripts in /example/scripts
if true
	@info("Skipped testing of example scripts") # default is to not run this test set
else
	@testset "example scripts" begin
		scriptsFolder = joinpath(dirname(pathof(JEMSS)), "..", "example", "scripts")

		# local search scripts, tested with a very small sim for speed
		cd("data/cities/small/6") do
			runGenConfig("gen_config.xml", overwriteOutputPath = true, doPrint = false)
			isdir("output") || mkdir("output")
			for script in ("deployment_local_search.jl", "deployment_search.jl", "deployment_test.jl",
					"priority_list_local_search.jl", "nested_comp_table_local_search.jl")
				@info(string("Running script: ", script))
				include(joinpath(scriptsFolder, script))
				@test true
				println()
			end
			rm("output", recursive = true)
		end

		# transient.jl
		cd("data/cities/small/2") do
			runGenConfig("gen_config.xml", overwriteOutputPath = true, doPrint = false)
			isdir("output") || mkdir("output")
			@info(string("Running script: transient.jl"))
			include(joinpath(scriptsFolder, "transient.jl"))
			@test true
			println()
			rm("output", recursive = true)
		end
	end
end

@testset "generate calls" begin
	cd("data/cities/small/3") do
		runGenConfig("gen_config_calls.xml", overwriteOutputPath = true, doPrint = false) # single calls file, limit calls by count
		runGenConfig("gen_config_calls_2.xml", overwriteOutputPath = true, doPrint = false) # single calls file, with maximum arrival time for last call
		runGenConfig("gen_config_calls_3.xml", overwriteOutputPath = true, doPrint = false) # single calls file, with minimum arrival time for last call
		genConfig = runGenConfig("gen_config_calls_multiple.xml", overwriteOutputPath = true, doPrint = false) # multiple calls files
		callSets = runGenConfigCalls(genConfig, doPrint = false, writeFile = false)
		@assert(isa(callSets, Vector{Vector{Call}}) && length(callSets) == genConfig.numCallsFiles)
		@test true
	end
end

@testset "set sim calls" begin
	sim = initSim("data/cities/small/1/sim_config.xml")
	genConfig = readGenConfig("data/cities/small/1/gen_config.xml")
	callSets = [makeCalls(genConfig) for i = 1:3]
	for calls in callSets
		for call in calls call.index = 0 end # test that setSimCalls! can handle this
		@assert(calls[1].nearestNodeIndex == nullIndex) # test that setSimCalls! can handle this
		setSimCalls!(sim, calls)
		@assert(sim.calls[1].index == 1)
		simulate!(sim)
	end
	@test true
end

@testset "sim replications" begin
	cd("data/cities/small/5") do
		runGenConfig("gen_config.xml", overwriteOutputPath = true, doPrint = false)
		sim = initSim("sim_config.xml")
		@assert(sim.numReps == 3)
		simulateReps!(sim.reps)
		@assert(all(rep -> rep.complete, sim.reps))
		writeStatsFiles(sim, sim.reps)
		resetReps!(sim.reps)
		@test true
	end
end

@testset "sim replications parallel" begin
	cd("data/cities/small/5") do
		runGenConfig("gen_config.xml", overwriteOutputPath = true, doPrint = false)
		sim = initSim("sim_config.xml")
		@assert(sim.numReps == 3)
		simulateReps!(sim.reps; parallel = true, numThreads = 2)
		@assert(all(rep -> rep.complete, sim.reps))
		writeStatsFiles(sim, sim.reps)
		resetReps!(sim.reps; parallel = true, numThreads = 2)
		@test true
	end
end

# test that mexclp can be solved, but not solution correctness
@testset "mexclp" begin
	@assert(isdir("data/cities/small/1/generated"))
	filename = joinpath(pwd(), "data/cities/small/1/mexclp/sim_config.xml")
	sim = initSim(filename);
	# solve mexclp with various kwarg values
	demandWeights = Dict([p => 1.0 for p in priorities])
	demandWeights[lowPriority] = 0
	stationsNumAmbs, deployment = solveMexclp!(sim; busyFraction =  0.4, demandWeights = demandWeights, stationCapacities = [2 for i = 1:sim.numStations])
	solveMexclp!(sim; busyFraction = -0.4) # test formulation where adding an ambulance at a station can reduce coverage
	# check that solution can be applied
	applyStationsNumAmbs!(sim, stationsNumAmbs)
	applyDeployment!(sim, deployment)
	@test true
end

@testset "cover bound" begin
	runGenConfig("data/cities/small/7/gen_config.xml", overwriteOutputPath = true, doPrint = false)
	sim = initSim("data/cities/small/7/sim_config.xml")
	genConfig = readGenConfig(sim.inputFiles["callGenConfig"].path)
	coverBoundSim = initCoverBoundSim(numReps = 30, numAmbs = sim.numAmbs, warmUpDuration = 0.1, minLastCallArrivalTime = 0.1 + 10,
		interarrivalTimeDistrRng = genConfig.interarrivalTimeDistrRng, ambBusyDurationSeed = 1)
	ambBusyDurationsToSample = [1:60;]/60/24 # 1 minute intervals, up to 1 hour
	queuedDurationsToSample = [0:maximum(sim.targetResponseDurations)*24*60;]/60/24 # 1 minute intervals
	coverBound = calcCoverBound!(sim; coverBoundSim = coverBoundSim, ambBusyDurationsToSample = ambBusyDurationsToSample, queuedDurationsToSample = queuedDurationsToSample)
	calcCoverBound!(coverBound; accountForQueuedDurations = false) # also calculate bound without accounting for call queueing durations (original cover bound)
	# @show coverBound.sim.bound # cover bound result
	simulateCoverBoundLowerBound!(coverBound) # also simulate a lower bound on the cover bound
	@test true
end


text:
using StructuresKit

#Verify Julia solution against Mathematica solution from
#Example 2 from Plaut, R.H., Moen, C.D.(2020). "Lateral-Torsional Deformations of Single-Span and Two-Span Thin-Walled Beams with Continuous Bracing".
#Journal of Constructional Steel Research.

#Z-section, no slope, single span, fixed-fixed
#kϕ=1500 N*mm/rad/mm, kx=0.1 N/mm^2, gravity load
#q = 1 kN/m
#Check again Figure 14 in the Plaut and Moen (2020) manuscript.
#  max ϕ

plaut_solution = 0.0164

#*********** Julia solution

#inputs

#Ix Iy Ixy J Cw
section_properties = [(3.230E6,449530,-865760, 397.09, 3.4104E9)]

#E  ν
material_properties = [(200,0.30)]

#ax ay
load_location = [(27.826,101.6)]

#kx kϕ
spring_stiffness = [(0.1/1000, 1500/1000)]

#ay_kx
spring_location = [(101.6)]

#member information
#L dL SectionProperties MaterialProperties LoadLocation SpringStiffness spring_location
member_definitions = [(7620,7620/48,1,1,1,1,1)]

#end boundary conditions
#type=1 u''=v''=ϕ''=0 (simply supported), type=2 u'=v'=ϕ'=0  (fixed), type=3 u''=v''=ϕ''=u'''=v'''=ϕ'''=0 (free end, e.g., a cantilever)
end_boundary_conditions = [2 2]

#supports
#location where u=v=ϕ=0
supports = [0.0 7620]

load = (0.0, 5.0/1000)  #kN/mm

z, u, v, ϕ, beam_properties = ThinWalledBeam.solve(member_definitions, section_properties, material_properties, spring_stiffness, spring_location, supports, load, load_location, end_boundary_conditions)

ϕmax = maximum(ϕ)

#calculate percent error as compared to Mathematica solution
error= abs.((ϕmax .- plaut_solution))./ plaut_solution

#accept 1% error from numerical solution
@test error <= 0.01


text:
@testset "Differential Common" begin
    @testset "Refine Differential" begin
        @test refine_differential(typeof(1.0 + 1im), Wirtinger(2,2)) == Wirtinger(2,2)
        @test refine_differential(typeof([1.0 + 1im]), Wirtinger(2,2)) == Wirtinger(2,2)

        @test refine_differential(typeof(1.2), Wirtinger(2,2)) == 4
        @test refine_differential(typeof([1.2]), Wirtinger(2,2)) == 4

        # For most differentials, in most domains, this does nothing
        for der in (DoesNotExist(), @thunk(23), @thunk(Wirtinger(2,2)), [1 2], One(), Zero(), 0.0)
            for 𝒟 in typeof.((1.0 + 1im, [1.0 + 1im], 1.2, [1.2]))
                @test refine_differential(𝒟, der) === der
            end
        end
    end
end


text:
module Remarkable
using HTTP
using JSON
using ZipFile
using Dates
using Crayons
using Parameters: @with_kw, type2dict
using UUIDs: uuid4
using AbstractTrees: AbstractTrees, print_tree
export  RemarkableClient,
        list_items,
        get_item,
        delete_item!,
        create_folder!,
        download_pdf,
        upload_pdf!
export print_tree

include("const.jl")
include("client.jl")
include("documents.jl")
include("api.jl")

end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible ideal GLM-MHD equations

equations = IdealGlmMhdEquations3D(1.4)

initial_condition = initial_condition_weak_blast_wave

surface_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
volume_flux  = (flux_hindenlang_gassner, flux_nonconservative_powell)
polydeg = 4
basis = LobattoLegendreBasis(polydeg)
indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=0.5,
                                         alpha_min=0.001,
                                         alpha_smooth=true,
                                         variable=density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)

solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux, volume_integral=volume_integral)

coordinates_min = (-2.0, -2.0, -2.0)
coordinates_max = ( 2.0,  2.0,  2.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=3,
                n_cells_max=10_000)

# create the semi discretization object
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

cfl = 1.4
stepsize_callback = StepsizeCallback(cfl=cfl)

glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        stepsize_callback,
                        glm_speed_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
using Test, Dates, Random, LinearAlgebra, Statistics
using StaticArrays
using Parameters
using OrdinaryDiffEq, DiffEqCallbacks
using JSON3, StructTypes, CSV
using Interpolations, NLsolve, FiniteDiff

using LongwaveModePropagator
using LongwaveModePropagator: QE, ME, distance, power
const LMP = LongwaveModePropagator

const TEST_RNG = MersenneTwister(1234)

const LWPC_PATH = "LWPC"

include("utils.jl")

#==
Scenarios
==#
const TEST_MODES = Dict{Any,Vector{EigenAngle}}()

const verticalB_scenario = (
    ea=EigenAngle(1.5 - 0.1im),
    bfield=BField(50e-6, π/2, 0),
    species=Species(QE, ME,
                    z->waitprofile(z, 75, 0.32; cutoff_low=40e3),
                    electroncollisionfrequency),
    ground=Ground(15, 0.001),
    tx=Transmitter(24e3),
    rx=GroundSampler(0:5e3:2000e3, Fields.Ez)
)

const multiplespecies_scenario = (
    ea=EigenAngle(1.5 - 0.1im),
    bfield=BField(50e-6, π/2, 0),
    species=(
        Species(QE, ME,
                z->waitprofile(z, 75, 0.32; cutoff_low=40e3),
                electroncollisionfrequency),
        Species(abs(QE), 58000*ME,  # used by LWPC
                z->waitprofile(z, 75, 0.32; cutoff_low=40e3),
                ioncollisionfrequency)
    ),
    ground=Ground(15, 0.001),
    tx=Transmitter(24e3),
    rx=GroundSampler(0:5e3:2000e3, Fields.Ez)
)

const isotropicB_resonant_scenario = (
    ea=EigenAngle(1.453098822238508 - 0.042008075239068944im),  # resonant
    bfield=BField(50e-6, 0, π/2),
    species=Species(QE, ME,
                    z->waitprofile(z, 75, 0.32; cutoff_low=40e3),
                    electroncollisionfrequency),
    ground=Ground(15, 0.001),
    tx=Transmitter(24e3),
    rx=GroundSampler(0:5e3:2000e3, Fields.Ez)
)

const resonant_scenario = (
    ea=EigenAngle(1.416127852502346 - 0.016482589477369265im),  # resonant
    bfield=BField(50e-6, deg2rad(68), deg2rad(111)),
    species=Species(QE, ME,
                    z->waitprofile(z, 75, 0.32; cutoff_low=40e3),
                    electroncollisionfrequency),
    ground=Ground(15, 0.001),
    tx=Transmitter(24e3),
    rx=GroundSampler(0:5e3:2000e3, Fields.Ez)
)

const resonant_elevatedrx_scenario = (
    ea=EigenAngle(1.416127852502346 - 0.016482589477369265im),  # resonant
    bfield=BField(50e-6, deg2rad(68), deg2rad(111)),
    species=Species(QE, ME,
                    z->waitprofile(z, 75, 0.32; cutoff_low=40e3),
                    electroncollisionfrequency),
    ground=Ground(15, 0.001),
    tx=Transmitter(24e3),
    rx=Sampler(0:5e3:2000e3, Fields.Ez, 8e3)
)

const resonant_horizontal_scenario = (
    ea=EigenAngle(1.416127852502346 - 0.016482589477369265im),  # resonant
    bfield=BField(50e-6, deg2rad(68), deg2rad(111)),
    species=Species(QE, ME,
                    z->waitprofile(z, 75, 0.32; cutoff_low=40e3),
                    electroncollisionfrequency),
    ground=Ground(15, 0.001),
    tx=Transmitter(24e3),
    rx=GroundSampler(0:5e3:2000e3, Fields.Ey)
)

const nonresonant_scenario = (
    ea=EigenAngle(1.5 - 0.1im),
    bfield=BField(50e-6, deg2rad(68), deg2rad(111)),
    species=Species(QE, ME,
                    z->waitprofile(z, 75, 0.32; cutoff_low=40e3),
                    electroncollisionfrequency),
    ground=Ground(5, 0.00005),
    tx=Transmitter(24e3),
    rx=GroundSampler(0:5e3:2000e3, Fields.Ez)
)

const interp_scenario = (
    ea=EigenAngle(1.5 - 0.1im),
    bfield=BField(50e-6, deg2rad(68), deg2rad(111)),
    species=Species(QE, ME,
        LinearInterpolation(0:500:110e3, waitprofile.(0:500:110e3, (75,), (0.32,); cutoff_low=40e3); extrapolation_bc=Line()),
        LinearInterpolation(0:500:110e3, electroncollisionfrequency.(0:500:110e3); extrapolation_bc=Line())),
    ground=Ground(5, 0.00005),
    tx=Transmitter(24e3),
    rx=GroundSampler(0:5e3:2000e3, Fields.Ez)
)

const homogeneousiono_scenario = (
    ea=EigenAngle(1.4161252139020892 - 0.016348911573820547im),  # resonant
    bfield=BField(50e-6, deg2rad(68), deg2rad(111)),
    species=Species(QE, ME,
                    z->2.65e6,
                    z->1e8),
    ground=Ground(15, 0.001),
    tx=Transmitter(24e3),
    rx=GroundSampler(0:5e3:2000e3, Fields.Ez)
)

const segmented_scenario = (
    ea=EigenAngle(1.5 - 0.1),
    bfield=[BField(50e-6, deg2rad(68), deg2rad(111)), BField(50e-6, deg2rad(68), deg2rad(111))],
    species=[Species(QE, ME,
                     z->waitprofile(z, 75, 0.32; cutoff_low=40e3),
                     electroncollisionfrequency),
             Species(QE, ME,
                     z->waitprofile(z, 80, 0.45; cutoff_low=40e3),
                     electroncollisionfrequency)],
    ground=[Ground(15, 0.001), Ground(15, 0.001)],
    tx=Transmitter(24e3),
    rx=GroundSampler(0:5e3:2000e3, Fields.Ez),
    distances=[0.0, 1000e3],
)

const θs = [complex(r,i) for r = range(deg2rad(60), deg2rad(89); length=30) for i =
            range(deg2rad(-10), deg2rad(0); length=11)]

maxabsdiff(a, b) = maximum(abs.(a - b))
meanabsdiff(a, b) = mean(abs.(a - b))

function findroots(scenario)
    @unpack bfield, species, ground, tx = scenario
    w = HomogeneousWaveguide(bfield, species, ground)
    me = PhysicalModeEquation(tx.frequency, w)
    origcoords = LMP.defaultmesh(tx.frequency)
    return findmodes(me, origcoords)
end

@testset "LongwaveModePropagator" begin
    include("EigenAngles.jl")
    include("Geophysics.jl")
    include("Waveguides.jl")

    include("magnetoionic.jl")
    include("TMatrix.jl")
    include("bookerquartic.jl")

    include("modefinder.jl")
    include("wavefields.jl")
    include("modeconversion.jl")
    include("modesum.jl")

    include("LongwaveModePropagator.jl")
    include("lwpc_comparisons.jl")

    include("IO.jl")
end


text:
# QuEST.jl/deps/build.jl
#
# Authors:
#  - Anabel Ovide, Uni Tartu
#  - Alejandro Villoria, Uni Tartu
#  - Dirk Oliver Theis, Ketita Labs & Uni Tartu
#
# MIT License
#
# (c) Ketita Labs, Uni Tartu, and the authors.
#
# 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.
#

#
# This script builds the package.
#

"""
True, if expert installation has been selected.  This variable changes the behavior of the
function `_auxBuild():
- if false, QuEST is built (cmake, make etc)
- if true, QuEST is not built.

In any case: QuEST is downloaded, and the C-interface is generated from the `.h`-files.
"""
const EXPERT_BUILD = haskey(ENV, "QUEST_JL_EXPERT_BUILD")

using Clang
using Clang.LibClang.Clang_jll

"""
Function `_bits_string()`

Transforms the precision *code* (1,2,4,...) to the number of bits as a string
("32,"64","2983",...).
"""
function _bits_string(;precision_code::Int) ::String
    @assert precision_code ∈ [1,2]
    precision_code==1  && return "32"
    precision_code==2  && return "64"
end

"""
Function `clang_from_h()` — Use Clang.jl to auto-generate C-interface frome `.h`-files

# Arguments:

"""
function _clang_from_h(makePrecision ::Int ; quest_root::String, isWindows ::Bool) ::Nothing
    precision = _bits_string(precision_code=makePrecision)

    questclang_include = joinpath(quest_root,"QuEST","include") |> normpath
    questclang_headers = [ joinpath(questclang_include,header)   for header in ["QuEST.h","QuEST_precision.h"] ] # ?readdir(questclang_include) if endswith(header, ".h")]
    wc = init(;
              headers           = questclang_headers,
              output_file       = joinpath(@__DIR__, "questclang_api_"*precision*".jl"),
              common_file       = joinpath(@__DIR__, "questclang_common_"*precision*".jl"),
              clang_includes    = vcat(questclang_include, CLANG_INCLUDE),
              clang_args        = ["-I", joinpath(questclang_include, ".."), "-DQuEST_PREC=$(makePrecision)"],
              header_wrapped    = (root, current)->root == current,
              header_library    = x->"questclang",
              clang_diagnostics = true,
              )
    run(wc)
end

"""
Execute commands to build QuEST
"""
function _auxBuild(makePrecision ::Int ; isWindows ::Bool) ::Nothing
    precision = _bits_string(precision_code=makePrecision)

    mkdir("build"*precision)
    cd("build"*precision)

    _clang_from_h(makePrecision ; quest_root="..", isWindows)

    if EXPERT_BUILD
        nothing
    else

        # `run()` with `wait=true` throws an error if anything goes wrong,
        # e.g., non-zero exit status
        if isWindows
            run(`cmake -DPRECISION=$makePrecision .. -G "MinGW Makefiles"`; wait=true)
        else
            run(`cmake -DPRECISION=$makePrecision ..` ; wait=true)
        end

        run(`make` ; wait=true)
    end

    cd("..")

    nothing
end

"""
Clone repository, call `_auxBuild()`
"""
function build(;isWindows::Bool) ::Nothing
    @info "Cloning QuEST..."
    run(`git clone https://github.com/TartuQC/QuEST.git` ; wait=true)
    cd("QuEST")
    @info "Build with 32-bit floats..."
    _auxBuild(1 ; isWindows)
    @info "Build with 64-bit floats..."
    _auxBuild(2 ; isWindows)
    @info "Build successful."
    cd("..")
end

#
# File "main" section
#

if EXPERT_BUILD
    @warn "Expert installation: Make sure the libraries `libQuEST_32` and `libQuEST_64` are loadable."
else
    @warn "Non-expert installation: Downloading & building QuEST with default settings"
end

if ispath("QuEST") && !isempty(readdir("QuEST"))
    rm("QuEST" ; force=true,recursive=true)
end

if Sys.isunix()
    build(isWindows=false)
elseif Sys.iswindows()
    build(isWindows=true)
else
    error("OS not supported.")
end

open("build_setup.jl", "w") do f
    write(f,
          "const EXPERT_BUILD = Bool($(EXPERT_BUILD))\n")
end

#EOF


text:
# -*- coding: utf-8 -*-
# DirectX

VERSION >= v"0.4.0-dev+6521" && __precompile__()
module DirectX

include("lib_fnc_defs.jl") # import DLL_Loader with import Relocator: _mf

import Relocator

export connect, close

const PSO_D3D, PSO_DEV, PSO_SPRITE, PSO_FONT, PSO_STRING, PSO_STRVBUF = 0:5
const res_default = (512, 512, "_string.png", "res")

type D3DMatrix # (fake to copy and read only access) in dx9adl.h
  aa::Float32; ba::Float32; ca::Float32; da::Float32
  ab::Float32; bb::Float32; cb::Float32; db::Float32
  ac::Float32; bc::Float32; cc::Float32; dc::Float32
  ad::Float32; bd::Float32; cd::Float32; dd::Float32
end

bitstype (8 * sizeof(D3DMatrix)) D3DMatrixBits # fake real byte size

D3DMatrix() = D3DMatrix(
  1., 0., 0., 0.,
  0., 1., 0., 0.,
  0., 0., 1., 0.,
  0., 0., 0., 1.)

function array_to(m::D3DMatrix, a::Array{Float32,2}) # 4x4 Array{Float32,2}
  m.aa = a[1, 1]; m.ba = a[2, 1]; m.ca = a[3, 1]; m.da = a[4, 1]
  m.ab = a[1, 2]; m.bb = a[2, 2]; m.cb = a[3, 2]; m.db = a[4, 2]
  m.ac = a[1, 3]; m.bc = a[2, 3]; m.cc = a[3, 3]; m.dc = a[4, 3]
  m.ad = a[1, 4]; m.bd = a[2, 4]; m.cd = a[3, 4]; m.dd = a[4, 4]
  m # load like a.'
end

type Q_D3DMatrix # in dx9adl.h
  tmp::Ptr{Void} # D3DMATRIX *
  world::Ptr{Void} # D3DMATRIX *
  view::Ptr{Void} # D3DMATRIX *
  proj::Ptr{Void} # D3DMATRIX *
end

type QQMatrix # in quaternion.h
  transform::Ptr{Void} # QUATERNIONIC_MATRIX *
  rotation::Ptr{Void} # QUATERNIONIC_MATRIX *
  scale::Ptr{Void} # QUATERNIONIC_MATRIX *
  translate::Ptr{Void} # QUATERNIONIC_MATRIX *
end

type VERTEX_GLYPH # in D3DxGlyph.h
  pQQM::Ptr{Void} # QQMATRIX * # in quaternion.h
  ppTexture::Ptr{Ptr{Void}} # LPDIRECT3DTEXTURE9 *
  ppVtxGlyph::Ptr{Ptr{Void}} # LPDIRECT3DVERTEXBUFFER9 *
  szGlyph::UInt32 # size_t
end

bitstype (8 * 4) FT_PosBits # fake real byte size # signed long
bitstype (8 * 4 * 2) FT_VectorBits # fake real byte size # {FT_Pos x, FT_Pos y}

type GLYPH_VTX # in D3DxFT2_types.h
  contour::Ptr{Void} # GLYPH_CONTOUR * # in D3DxFT2_types.h
  cmx::UInt32
  ctc::UInt32
  buf::Ptr{Void} # CUSTOMCUV * # in D3DxFT2_types.h
  bmx::UInt32
  btc::UInt32
  p::FT_VectorBits
  pz::FT_PosBits
  z::FT_PosBits
  col::UInt32
  bgc::UInt32
end

type GLYPH_TBL # in D3DxFT2_types.h
  pIS::Ptr{Void} # RENDERD3DITEMSSTATE * # in dx9adl.h
  pVG::Ptr{Void} # VERTEX_GLYPH * # in D3DxGlyph.h
  facename::Ptr{Cchar}
  utxt::Ptr{Cwchar_t}
  ratio::Float32
  angle::Float32
  reserved0::UInt32
  reserved1::UInt32
  mode::UInt32
  td::UInt32
  tw::UInt32
  th::UInt32
  rct::Ptr{Void} # CUSTOMRECT * # in D3DxFT2_types.h
  glyphBmp::Ptr{Void} # BOOL (*)(GLYPH_TBL *, FT_Bitmap *, FT_Int, FT_Int)
  matrix::Ptr{Void} # QUATERNIONIC_MATRIX * # in quaternion.h
  vtx::Ptr{Void} # GLYPH_VTX * # in D3DxFT2_types.h
  funcs::Ptr{Void} # FT_Outline_Funcs *
  glyphContours::Ptr{Void} # BOOL (*)(GLYPH_TBL *)
  glyphAlloc::Ptr{Void} # BOOL (*)(GLYPH_TBL *)
  glyphFree::Ptr{Void} # BOOL (*)(GLYPH_TBL *)
end

type D3DVector
  x::Float32; y::Float32; z::Float32
end

bitstype (8 * sizeof(D3DVector)) D3DVectorBits # fake real byte size

type D9F_Vecs # (construction values of matrices) in dx9adl.h
  eyePt::D3DVectorBits
  reserved0::UInt32
  lookatPt::D3DVectorBits
  reserved1::UInt32
  upVec::D3DVectorBits
  reserved2::UInt32
  fovY::Float32; aspect::Float32; zn::Float32; zf::Float32

#  function D9F_Vecs(
#    eyePt::D3DVector, reserved0::UInt32,
#    lookatPt::D3DVector, reserved1::UInt32,
#    upVec::D3DVector, reserved2::UInt32,
#    fovY::Float32, aspect::Float32, zn::Float32, zf::Float32
#  )
  function D9F_Vecs(eyePt, reserved0, lookatPt, reserved1, upVec, reserved2,
    fovY, aspect, zn, zf)
    return new(
      (@as_bits D3DVectorBits eyePt), reserved0,
      (@as_bits D3DVectorBits lookatPt), reserved1,
      (@as_bits D3DVectorBits upVec), reserved2,
      fovY, aspect, zn, zf)
  end
end

type D9Foundation # (fake to copy and read only access) in dx9adl.h
  pD3Dpp::Ptr{Ptr{Void}} # D3DPRESENT_PARAMETERS *
  pD3D::Ptr{Void} # LPDIRECT3D9
  pDev::Ptr{Void} # LPDIRECT3DDEVICE9
  pMenv::Ptr{Void} # Q_D3DMATRIX * # in dx9adl.h
  pVecs::Ptr{Void} # D9F_VECS * # in dx9adl.h
  imstring::Ptr{Cchar}
  imw::UInt32
  imh::UInt32

  function D9Foundation()
    return new()
  end
end

type RenderD3DItemsState # in dx9adl.h
  parent::Ptr{Void}
  d9fnd::Ptr{D9Foundation}
  sysChain::Ptr{Ptr{Void}}
  usrChain::Ptr{Ptr{Void}}
  smx::UInt32
  umx::UInt32
  width::UInt32
  height::UInt32
  bgc::UInt32
  fgc::UInt32
  mode::UInt32
  stat::UInt32
  hWnd::UInt32 # HWND
  fps::UInt32
  prevTime::UInt32
  nowTime::UInt32
end

m_tmp = D3DMatrix()
m_world = D3DMatrix()
m_view = D3DMatrix()
m_proj = D3DMatrix()
menv = Q_D3DMatrix(C_NULL, C_NULL, C_NULL, C_NULL) # re-set later
vecs = D9F_Vecs(
  D3DVector(5., 5., -5.), 0,
  D3DVector(0., 0., 0.), 0,
  D3DVector(0., 1., 0.), 0,
  3./4., 1., 1., 100.)
d9fnd = D9Foundation() # holder (must be out of struct Dx9adl ?)
istat = RenderD3DItemsState(C_NULL, C_NULL, C_NULL, C_NULL, # re-set later
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

type Dx9adl
  basepath::AbstractString # base path
  respath::AbstractString # resource path
  ims::AbstractString # to hold the pointer placing dynamic char[] (anti GC)
  d9fnd::D9Foundation # reference
  istat::RenderD3DItemsState # reference

  function Dx9adl(w::Int, h::Int, bp::AbstractString)
    res = Relocator.searchResDll(bp, res_default[4], true)
    ims = replace(res * "/" * res_default[3], "/", "\\") # only for Windows
    # OK pointer(ims) # AbstractString to Cchar
    # OK pointer(ims.data) # Array{UInt8,1} to Cchar
    # BAD @ptr ims
    # BAD @ptr ims.data
    # OK (@ptr ims.data) + 32 # OK but wrong way
    d9fnd.imstring = pointer(ims)
    d9fnd.imw = res_default[1]
    d9fnd.imh = res_default[2]
    menv.tmp = @ptr m_tmp
    menv.world = @ptr m_world
    menv.view = @ptr m_view
    menv.proj = @ptr m_proj
    d9fnd.pMenv = @ptr menv
    d9fnd.pVecs = @ptr vecs
    istat.d9fnd = @ptr d9fnd # or set C_NULL
    istat.width = w
    istat.height = h
    istat.bgc = 0xFFFF80FF
    istat.fgc = 0x80EE66CC
    istat.mode = 0x0CC00000
    # set mode 0 to skip debugalloc/debugfree
    d = new(bp, res, ims, d9fnd, istat) # set parent later
    d.istat.parent = @ptr d
    return d
  end
end

function connect(w::Int, h::Int, bp::AbstractString="")
  DLL_Loader.load(bp)
# debugalloc() # needless to call on Julia console ?
  return Dx9adl(w, h, bp)
end

function close(d9::Dx9adl)
# debugfree() # can not call on Julia console ?
  DLL_Loader.unload()
  return true
end

end


text:
function Vert(PP :: Array{VerticeType,1}, NewCut,CutNum)
    const EPS = 0;
    V_Plus = [];
    v_Plus = 0;
    V_Minus = [];
    v_Minus = 0;
    for i = 1:length(PP)
        if (sum(NewCut[1:end-1].*PP[i].Point) + NewCut[end] > EPS)
            push!(V_Plus, i)
            v_Plus = v_Plus + 1;
        else
            push!(V_Minus, i)
            v_Minus = v_Minus + 1;
        end
    end
    # println(V_Plus)
    # println(V_Minus)
# ---------------------------------------------------------------------------
    # Step I
    if v_Minus == 0
        error("The cut makes it infeasible!!!!")
        return []
    elseif v_Plus == 0
        println("The cut is redundant!!!!")
##############################################################
        return []
##########################################################3##
        return PP
    end

    if v_Plus <= v_Minus
        V1 = V_Plus;
        V2 = V_Minus;
    else
        V2 = V_Plus;
        V1 = V_Minus;
    end
# ---------------------------------------------------------------------------
    # Step II

    VS = [];
    counter = 1;
    for u in V1
        for v in PP[u].adj
            if v in V2
                W = FindIntersectingPoint(PP[u].Point,PP[v].Point,NewCut)
                NewPoint = VerticeType([],[],[],-1)
                NewPoint.Point = W;
                if v_Plus <= v_Minus
                    PP[v].adj = filter(i -> (i != u), PP[v].adj);
                    push!(PP[v].adj,length(PP)+counter);
                    counter = counter + 1;
                    push!(NewPoint.adj,v);
                else
                    PP[u].adj = filter(i -> (i != v), PP[u].adj);
                    push!(PP[u].adj,length(PP)+counter)
                    counter = counter + 1;
                    push!(NewPoint.adj,u)
                end
                temp = intersect(PP[u].J,PP[v].J)
                for ii in temp
                    push!(NewPoint.J,ii)
                end
                push!(NewPoint.J,CutNum)
                if isempty(VS)
                    VS = [NewPoint];
                else
                    push!(VS,NewPoint);
                end
            end
        end
    end
# ---------------------------------------------------------------------------
    IndTrans = Dict();
    for i = 1:v_Minus
        IndTrans[V_Minus[i]]=i;
        for j = 1:length(VS)
            VS[j].adj[VS[j].adj .== V_Minus[i]] =i;
        end
    end

    for i in V_Minus
        for j = 1:length(PP[i].adj)
            if PP[i].adj[j] > length(PP)
                PP[i].adj[j] = PP[i].adj[j]- v_Plus;
            else
                PP[i].adj[j] = IndTrans[PP[i].adj[j]];
            end
        end
    end
    PP = [PP[i] for i in V_Minus];

    # # Step III
    for u = 1:length(VS), v = 1:length(VS)
        if length(intersect(VS[u].J,VS[v].J)) == length(PP[1].Point)-1
            if !(length(V_Minus)+v in VS[u].adj)
                push!(VS[u].adj,length(V_Minus)+v)
            end
            if !(length(V_Minus)+u in VS[v].adj)
                push!(VS[v].adj,length(V_Minus)+u)
            end
        end
    end

    for i = 1:length(VS)
        VS[i].T = prod(VS[i].Point);
        push!(PP,VS[i])
    end
    return PP;

end


text:
using SuperFastGraphs
using Base.Test
using LightGraphs

gDWeakly = DiGraph(6)

add_edge!(gDWeakly, 1, 2)
add_edge!(gDWeakly, 1, 3)
add_edge!(gDWeakly, 1, 4)
add_edge!(gDWeakly, 5, 1)
add_edge!(gDWeakly, 2, 3)
add_edge!(gDWeakly, 3, 6)
add_edge!(gDWeakly, 4, 5)
println(is_directed(gDWeakly) && !is_strongly_connected(gDWeakly) && is_weakly_connected(gDWeakly))
println(preProcess(gDWeakly))

println("-----------------")
gDStrongly = DiGraph(3)
add_edge!(gDStrongly, 1, 2)
add_edge!(gDStrongly, 2, 3)
add_edge!(gDStrongly, 3, 1)
println(is_directed(gDStrongly) && is_strongly_connected(gDStrongly))
println(preProcess(gDStrongly))

println("---------------")

gNotD_NotC = Graph(6)
add_edge!(gNotD_NotC, 1, 2)
add_edge!(gNotD_NotC, 1, 3)
add_edge!(gNotD_NotC, 2, 4)
add_edge!(gNotD_NotC, 2, 5)
add_edge!(gNotD_NotC, 5, 4)
println("connesso: ",is_connected(gNotD_NotC))
println("componenti connesse: ", connected_components(gNotD_NotC))
println(!is_directed(gNotD_NotC) && !is_connected(gNotD_NotC))
println(preProcess(gNotD_NotC))

println("-------------------")
# g = loadgraph("testdata/karate.lg", "graph")
g = loadgraph("testdata/internet.lg", "graph")
println(g)
tic()
topKcc(g, 33)
toc()
# println(fastClosenessCentrality(g, 5)[1])
println("-----------------------")
tic()
closeness_centrality(g)
toc()
# println(sort!(closeness_centrality(g))[end-5])
#println(connected_components(g))

max_stars_repo_name: mameli/SuperFastGraphs

max_issues_repo_name: mameli/SuperFastGraphs

max_forks_repo_name: mameli/SuperFastGraphs


text:
import MPI

export MPIComm

"""
    MPIComm()
    MPIComm(comm::MPI.Comm)

An implementation of Comm using MPI
The no argument constructor uses MPI.COMM_WORLD
"""
struct MPIComm{GID <: Integer, PID <:Integer, LID <: Integer} <: Comm{GID, PID, LID}
    mpiComm::MPI.Comm
end

function MPIComm(GID::Type, PID::Type, LID::Type)
    MPIInit()
    comm = MPIComm{GID, PID, LID}(MPI.COMM_WORLD)

    comm
end

MPINeedsInitialization = true

"""
    MPIInit()

On the first call, initializes MPI and adds an exit hook to finalize MPI
Does nothing on subsequent calls
"""
function MPIInit()
    global MPINeedsInitialization
    if MPINeedsInitialization
        MPI.Init()
        atexit(() -> MPI.Finalize())

        MPINeedsInitialization = false
    end
end

function barrier(comm::MPIComm)
    MPI.Barrier(comm.mpiComm)
end

function broadcastAll(comm::MPIComm, myvals::AbstractArray{T}, root::Integer)::Array{T} where T
    vals = copy(myvals)
    result = MPI.Bcast!(vals, root-1, comm.mpiComm)
    result
end

function gatherAll(comm::MPIComm, myVals::AbstractArray{T})::Array{T} where T
    lengths = MPI.Allgather([convert(Cint, length(myVals))], comm.mpiComm)
    MPI.Allgatherv(myVals, lengths, comm.mpiComm)
end

function sumAll(comm::MPIComm, partialsums::AbstractArray{T})::Array{T} where T
    MPI.allreduce(partialsums, +, comm.mpiComm)
end

function maxAll(comm::MPIComm, partialmaxes::AbstractArray{T})::Array{T} where T
    MPI.allreduce(partialmaxes, max, comm.mpiComm)
end

function maxAll(comm::MPIComm, partialmaxes::AbstractArray{Bool})::Array{Bool}
    Array{Bool}(maxAll(comm, Array{UInt8}(partialmaxes)))
end

function minAll(comm::MPIComm, partialmins::AbstractArray{T})::Array{T} where T
    MPI.allreduce(partialmins, min, comm.mpiComm)
end

function minAll(comm::MPIComm, partialmins::AbstractArray{Bool})::Array{Bool}
    Array{Bool}(minAll(comm, Array{UInt8}(partialmins)))
end

function scanSum(comm::MPIComm, myvals::AbstractArray{T})::Array{T} where T
    MPI.Scan(myvals, length(myvals), MPI.SUM, comm.mpiComm)
end

function myPid(comm::MPIComm{GID, PID})::PID where {GID <: Integer, PID <: Integer}
    MPI.Comm_rank(comm.mpiComm) + 1
end

function numProc(comm::MPIComm{GID, PID})::PID where {GID <: Integer, PID <:Integer}
    MPI.Comm_size(comm.mpiComm)
end

function createDistributor(comm::MPIComm{GID, PID, LID})::MPIDistributor{GID, PID, LID}  where {GID <: Integer, PID <: Integer, LID <: Integer}
    MPIDistributor(comm)
end


text:
export LangC_MKL, LangC_OpenBLAS

# Data structures for the language.
"""
     LangC_MKL(gen_main::Bool)

Code generation in C using MKL implementation of BLAS.
"""
struct LangC_MKL
    gen_main::Any
    overwrite_input::Any
end

"""
     LangC_OpenBLAS(gen_main::Bool)

Code generation in C using OpenBLAS implementation of BLAS.
"""
struct LangC_OpenBLAS
    gen_main::Any
    overwrite_input::Any
end
LangC = Union{LangC_MKL,LangC_OpenBLAS}

# Generated constructors
# https://docs.julialang.org/en/v1/manual/metaprogramming/#Code-Generation
# to default to gen_main=false, overwrite_input=true
for L in [:LangC_MKL, :LangC_OpenBLAS]
    @eval $L() = $L(false, true)
    @eval $L(gen_main) = $L(gen_main, true)
end

# Language specific operations.
comment(::LangC, s) = "/* $s */"

slotname(::LangC, i) = "memslots[$(i-1)]"

# Variable declaration, initialization, and reference.
# In C, complex types are structures and are passed by reference.
function assign_coeff(::LangC_MKL, val::T, i) where {T<:Complex}
    assignment_string = [
        "coeff$i.real = " * string(real(val)) * ";",
        "coeff$i.imag = " * string(imag(val)) * ";",
    ]
    variable_string = "&coeff$i"
    return (variable_string, assignment_string)
end
function assign_coeff(::LangC_OpenBLAS, val::T, i) where {T<:Complex}
    assignment_string =
        ["coeff$i = " * string(real(val)) * " + " * string(imag(val)) * "*I;"]
    variable_string = "&coeff$i"
    return (variable_string, assignment_string)
end
function assign_coeff(::LangC, val::T, i) where {T<:Real}
    assignment_string = ["coeff$i = $val;"]
    variable_string = "coeff$i"
    return (variable_string, assignment_string)
end

# Constant declaration and reference.
function declare_constant(::LangC_MKL, val::T, id, type) where {T<:Complex}
    return "const $type $id = {.real = " *
           string(real(val)) *
           ", " *
           ".imag = " *
           string(imag(val)) *
           "};"
end
function declare_constant(::LangC_OpenBLAS, val::T, id, type) where {T<:Complex}
    return "const $type $id = " *
           string(real(val)) *
           " + " *
           string(imag(val)) *
           "*I;"
end
function declare_constant(::LangC, val::T, id, type) where {T<:Real}
    return "const $type $id = $val;"
end
reference_constant(::LangC, T, id) = (T <: Complex) ? "&$id" : "$id"

# Code generation.
# Return MKL/OpenBLAS-specific includes.
function get_blas_includes(::LangC_MKL)
    return "#include<mkl/mkl.h>"
end
function get_blas_includes(::LangC_OpenBLAS)
    return "#include<cblas.h>\n#include<lapacke.h>"
end

# Return type and BLAS/LAPACK prefix corresponding to type T.
function get_blas_type(::LangC, T::Type{Float32})
    return ("float", "s")
end
function get_blas_type(::LangC, T::Type{Float64})
    return ("double", "d")
end
function get_blas_type(::LangC_MKL, T::Type{Complex{Float32}})
    return ("MKL_Complex8", "c")
end
function get_blas_type(::LangC_MKL, T::Type{Complex{Float64}})
    return ("MKL_Complex16", "z")
end
function get_blas_type(::LangC_OpenBLAS, T::Type{Complex{Float32}})
    return ("openblas_complex_float", "c")
end
function get_blas_type(::LangC_OpenBLAS, T::Type{Complex{Float64}})
    return ("openblas_complex_double", "z")
end

function function_definition(lang::LangC, graph, T, funname, precomputed_nodes)
    (blas_type, blas_prefix) = get_blas_type(lang, T)
    code = init_code(lang)
    push_code!(code, get_blas_includes(lang), ind_lvl = 0)
    push_code!(code, "#include<assert.h>", ind_lvl = 0)
    push_code!(code, "#include<stdlib.h>", ind_lvl = 0)
    push_code!(code, "#include<string.h>", ind_lvl = 0)
    push_code!(code, "")
    push_comment!(code, "Code for matrix function evaluation.", ind_lvl = 0)
    input_variables =
        join(map(x -> "$blas_type *" * string(x), precomputed_nodes), ", ")
    push_code!(
        code,
        "void $blas_prefix$funname($input_variables, " *
        "const size_t n, $blas_type *output) {",
        ind_lvl = 0,
    )
    return code
end

function init_mem(lang::LangC, max_nof_nodes, precomputed_nodes)
    mem = CodeMem(max_nof_nodes, i -> slotname(lang, i))

    # Make sure the precomputed nodes have memslots
    for (i, n) in enumerate(precomputed_nodes)
        alloc_slot!(mem, i, n)
    end

    return mem
end

function function_init(lang::LangC, T, mem, graph, precomputed_nodes)
    (blas_type, blas_prefix) = get_blas_type(lang, T)
    code = init_code(lang)
    max_nodes = size(mem.slots, 1)

    # Initialization.
    graph_ops = values(graph.operations)
    # memalloc: matrices to allocate
    # memslots: numbers of pointers
    num_precomputed_nodes = size(precomputed_nodes, 1)
    num_memalloc =
        lang.overwrite_input ? max_nodes - num_precomputed_nodes : max_nodes
    num_preallocated_slots = max_nodes - num_memalloc
    push_code!(code, "size_t max_memalloc = $num_memalloc;")
    push_code!(code, "size_t max_memslots = $max_nodes;")
    push_code!(code, "")

    push_comment!(code, "Declarations and initializations.")
    # Allocate coefficients for linear combinations, if needed.
    if :lincomb in graph_ops
        push_code!(code, "$blas_type coeff1, coeff2;")
    end
    # Declare constants ZERO and ONE, if needed.
    if :mult in graph_ops
        push_code!(
            code,
            declare_constant(lang, convert(T, 0.0), "ZERO", blas_type),
        )
    end
    if :lincomb in graph_ops || :mult in graph_ops
        push_code!(
            code,
            declare_constant(lang, convert(T, 1.0), "ONE", blas_type),
        )
    end
    # Array ipiv for pivots of GEPP (needed only if graph has linear systems).
    if :ldiv in graph_ops
        push_code!(code, "lapack_int *ipiv = malloc(n*sizeof(*ipiv));")
    end
    push_code!(code, "size_t j;")
    push_code!(code, "")

    push_comment!(code, "Memory management.")
    push_code!(
        code,
        "$blas_type *master_mem = malloc(n*n*max_memalloc" *
        "*sizeof(*master_mem));",
    )
    push_code!(
        code,
        "$blas_type *memslots[max_memslots]; " *
        comment(lang, "As many slots as nodes"),
    )

    # Initialize the inputs
    for (i, n) in enumerate(precomputed_nodes)
        Ak_slot_name = get_slot_name(mem, n)
        if (lang.overwrite_input)
            # Just assign a pointer to the slot to allow overwrite
            push_code!(
                code,
                "$Ak_slot_name = $n; " * comment(lang, "Overwrite $n"),
            )
        else
            # Otherwise make a copy
            push_code!(
                code,
                "memcpy($Ak_slot_name,$n,n*n*sizeof(*master_mem));",
            )
        end
    end

    # Initialize pointers.
    push_comment!(code, "The other slots are pointers to allocated memory.")
    push_code!(code, "for (j=$(num_preallocated_slots-1); j<max_memalloc; j++)")
    push_code!(code, "memslots[j+1] = master_mem+j*n*n;", ind_lvl = 2)

    # Store identity explicitly only if graph has a linear combination of I.
    if has_identity_lincomb(graph)
        push_comment!(code, "Graph has linear combination of identities.")
        push_comment!(code, "The matrix I is explicitly allocated.")
        alloc_slot!(mem, num_precomputed_nodes + 1, :I)
        nodemem = get_slot_name(mem, :I)
        push_code!(code, "memset($nodemem, 0, n*n*sizeof(*master_mem));")
        push_code!(code, "for(j=0; j<n*n; j+=n+1)")
        push_code!(code, "*($nodemem+j) = ONE;", ind_lvl = 2)
    end

    return code
end

function function_end(lang::LangC, graph, mem)
    code = init_code(lang)
    retval_node = graph.outputs[end]
    retval = get_slot_name(mem, retval_node)
    push_code!(code, "")
    push_comment!(code, "Prepare output.")
    push_code!(code, "memcpy(output, $retval, n*n*sizeof(*output));")
    if :ldiv in values(graph.operations)
        push_code!(code, "free(ipiv);")
    end
    if size(mem.slots, 1) > 0
        push_code!(code, "free(master_mem);")
    end
    push_code!(code, "}", ind_lvl = 0)
    return code
end

function execute_operation!(lang::LangC, T, graph, node, dealloc_list, mem)
    (blas_type, blas_prefix) = get_blas_type(lang, T)

    op = graph.operations[node]
    parent1 = graph.parents[node][1]
    parent2 = graph.parents[node][2]

    # Keep deallocation list (used for smart memory management).
    dealloc_list = deepcopy(dealloc_list)
    setdiff!(dealloc_list, keys(mem.special_names))

    # Check if parents have a memory slot.
    if has_identity_lincomb(graph) # Linear combination of I.
        p1_is_identity = p2_is_identity = false
    else
        p1_is_identity = parent1 == :I
        p2_is_identity = parent2 == :I
        setdiff!(dealloc_list, [:I])
    end

    # Get memory slots of parents, if explicitly stored.
    if !p1_is_identity
        parent1mem = get_slot_name(mem, parent1)
    end
    if !p2_is_identity
        parent2mem = get_slot_name(mem, parent2)
    end

    rzero = reference_constant(lang, T, "ZERO")
    rone = reference_constant(lang, T, "ONE")

    code = init_code(lang)
    push_code!(code, "")
    push_comment!(code, "Computing $node with operation: $op")
    if op == :mult
        (nodemem_i, nodemem) = get_free_slot(mem)
        alloc_slot!(mem, nodemem_i, node)
        push_code!(
            code,
            "cblas_$blas_prefix" *
            "gemm(CblasColMajor, " *
            "CblasNoTrans, CblasNoTrans, n, n, n,\n" *
            "            $rone, $parent1mem, n, $parent2mem, n,\n" *
            "            $rzero, $nodemem, n);",
        )
    elseif op == :ldiv
        # Find a memory slot for the LU factors.
        # As ?detrs computes the LU decomposition of parent1 in-place, we need
        # to make a copy of $parent1mem, unless parent1 is on the deallocation
        # list.
        if (parent1 in dealloc_list)
            push_comment!(code, "Reusing memory of $parent1 for LU factors.")
            lhsmem = parent1mem
            lhsmem_i = get_slot_number(mem, parent1)
            # Remove parent 1 from the deallocation list, but only temporarily.
            # This is to make sure that the LU factors do not get overwritten by
            # the result of the sytem solve.
            setdiff!(dealloc_list, [parent1])
        else
            (lhsmem_i, lhsmem) = get_free_slot(mem)
            alloc_slot!(mem, lhsmem_i, :dummy)
            push_code!(
                code,
                "memcpy($lhsmem, $parent1mem, " * "n*n*sizeof(*master_mem));",
            )
        end

        # Compute LU decomposition of parent1.
        push_code!(
            code,
            "LAPACKE_$blas_prefix" *
            "getrf(LAPACK_COL_MAJOR, n, n, " *
            "$lhsmem, n, ipiv);",
        )

        # Solve the linear systems. We have two cases:
        if parent2 == :I
            # 1. If parent2 is the identity, we use the function ?getri, which
            # overwrites the LU factors with the matrix inverse.

            # Compute matrix inverse.
            push_code!(
                code,
                "LAPACKE_$blas_prefix" *
                "getri(LAPACK_COL_MAJOR, " *
                "n, $lhsmem, n, ipiv);",
            )

            # The LU factors should not be deallocated in this case, but we have
            # to point the slot lhsmem to the current node.
            alloc_slot!(mem, lhsmem_i, node)
        else
            # 2. If parent2 is not the identity, we use the function ?getrs,
            # which # overwrites B in AX = B. In this case, we need to copy
            # $parent2mem to # $nodemem first, unless parent2 is in the
            # deallocation list.

            # Allocate slot for result.
            if (parent2 in dealloc_list)
                push_comment!(code, "Reusing memory of $parent2 for solution.")
                nodemem = get_slot_name(mem, parent2)
                setdiff!(dealloc_list, [parent2])
                set_slot_number!(mem, get_slot_number(mem, parent2), node)
            else
                (nodemem_i, nodemem) = get_free_slot(mem)
                alloc_slot!(mem, nodemem_i, node)
                push_code!(
                    code,
                    "memcpy($nodemem, $parent2mem, " *
                    "n*n*sizeof(*master_mem));",
                )
            end

            # Solve linear system.
            push_code!(
                code,
                "LAPACKE_$blas_prefix" *
                "getrs(LAPACK_COL_MAJOR, " *
                "'N', n, n,\n" *
                "               $lhsmem, n, ipiv,\n" *
                "               $nodemem, n);",
            )

            # Deallocate LU factors.
            # TODO It might make sense to keep them if another system with the
            # same left-hand side is still in the graph.
            free!(mem, lhsmem_i)
        end

    elseif op == :lincomb
        coeff_vars = Vector{String}(undef, 2)

        # Coefficients of graph should have type T.
        (coeff1, coeff1_code) = assign_coeff(lang, graph.coeffs[node][1], 1)
        for s in coeff1_code
            push_code!(code, s)
        end
        (coeff2, coeff2_code) = assign_coeff(lang, graph.coeffs[node][2], 2)
        for s in coeff2_code
            push_code!(code, s)
        end

        # Reuse parent in deallocation list for output (saves memcopy).
        if ((parent1 in dealloc_list) || (parent2 in dealloc_list))
            recycle_parent = (parent1 in dealloc_list) ? parent1 : parent2
            push_comment!(code, "Smart lincomb recycle $recycle_parent")

            # Perform operation.
            nodemem = get_slot_name(mem, recycle_parent)
            if (recycle_parent == parent1)
                if p2_is_identity
                    push_code!(
                        code,
                        "cblas_$blas_prefix" *
                        "scal(n*n, $coeff1, " *
                        "$nodemem, 1);",
                    )
                    push_code!(
                        code,
                        "cblas_$blas_prefix" *
                        "axpby(n, " *
                        "$coeff2, &ONE, 0,\n" *
                        "             $rone, $nodemem, n+1);",
                    )
                else
                    push_code!(
                        code,
                        "cblas_$blas_prefix" *
                        "axpby(n*n, " *
                        "$coeff2, $parent2mem, 1,\n" *
                        "             $coeff1, $nodemem, 1);",
                    )
                end
            else
                if p1_is_identity
                    push_code!(
                        code,
                        "cblas_$blas_prefix" *
                        "scal(n*n, $coeff2, " *
                        "$nodemem, 1);",
                    )
                    push_code!(
                        code,
                        "cblas_$blas_prefix" *
                        "axpby(n, " *
                        "$coeff1, &ONE, 0,\n" *
                        "             $rone, $nodemem, n+1);",
                    )
                else
                    push_code!(
                        code,
                        "cblas_$blas_prefix" *
                        "axpby(n*n, " *
                        "$coeff1, $parent1mem, 1,\n" *
                        "             $coeff2, $nodemem, 1);",
                    )
                end
            end

            # Remove output from deallocation list.
            setdiff!(dealloc_list, [recycle_parent])

            # Update CodeMem.
            j = get_slot_number(mem, recycle_parent)
            set_slot_number!(mem, j, node)
        else
            # Allocate new slot for result.
            (nodemem_i, nodemem) = get_free_slot(mem)
            alloc_slot!(mem, nodemem_i, node)

            # Compute linear combination.
            if p1_is_identity
                push_code!(
                    code,
                    "memcpy($nodemem, $parent2mem, " *
                    "n*n*sizeof(*master_mem));",
                )
                push_code!(
                    code,
                    "cblas_$blas_prefix" *
                    "scal(n*n, $coeff2, " *
                    "$nodemem, 1);",
                )
                push_code!(
                    code,
                    "cblas_$blas_prefix" *
                    "axpby(n, " *
                    "$coeff1, &ONE, 0,\n" *
                    "             $rone, $nodemem, n+1);",
                )
            elseif p2_is_identity
                push_code!(
                    code,
                    "memcpy($nodemem, $parent1mem, " *
                    "n*n*sizeof(*master_mem));",
                )
                push_code!(
                    code,
                    "cblas_$blas_prefix" *
                    "scal(n*n, $coeff1, " *
                    "$nodemem, 1);",
                )
                push_code!(
                    code,
                    "cblas_$blas_prefix" *
                    "axpby(n*n, " *
                    "$coeff2, &ONE, 0,\n" *
                    "             $rone, $nodemem, n+1);",
                )
            else
                push_code!(
                    code,
                    "memcpy($nodemem, $parent2mem, " *
                    "n*n*sizeof(*master_mem));",
                )
                push_code!(
                    code,
                    "cblas_$blas_prefix" *
                    "axpby(n*n, " *
                    "$coeff1, $parent1mem, 1,\n" *
                    "             $coeff2, $nodemem, 1);",
                )
            end
        end
    else
        error("Unknown operation")
    end

    # Deallocate
    for n in dealloc_list
        i = get_slot_number(mem, n)
        push_comment!(code, "Deallocating $n in slot $i")
        free!(mem, i)
    end

    nodemem = 0
    return (code, nodemem)
end

function scalar_to_string(::LangC, z)
    return "$z"
end
function scalar_to_string(::LangC_OpenBLAS, z::T) where {T<:Complex}
    return "$(real(z)) + $(imag(z))*I"
end
function scalar_to_string(::LangC_MKL, z::T) where {T<:Complex}
    return "{$(real(z)), $(imag(z))}"
end
function print_indented_matrix(lang::LangC, code, A; ind_lvl = 1)
    (m, n) = size(A)
    for j = 1:n
        column_string = ""
        for i = 1:m
            column_string *=
                scalar_to_string(lang, A[i, j]) * (i != m ? ", " : "")
        end
        column_string *= (j != n ? "," : "")
        push_code!(code, column_string, ind_lvl = ind_lvl)
    end
end

function compilation_string(::LangC_OpenBLAS, fname)
    return "gcc -o main_compiled $fname -labials -llapacke"
end
function compilation_string(::LangC_MKL, fname)
    return "gcc -o main_compiled $fname -lmkl_rt"
end

function gen_main(lang::LangC, T, fname, funname; A = 10::Union{Integer,Matrix})
    code = init_code(lang)
    if (lang.gen_main)
        (blas_type, blas_prefix) = get_blas_type(lang, T)
        push_code!(code, "")
        push_code!(code, "")
        push_code!(code, "")
        push_code!(code, "typedef $blas_type blas_type;", ind_lvl = 0)
        push_code!(code, "")

        # Main function starts here.
        push_comment!(
            code,
            "Code snippet that calls $blas_prefix$funname().",
            ind_lvl = 0,
        )
        push_comment!(
            code,
            "With the GNU Compiler Collection, compile with:",
            ind_lvl = 0,
        )
        push_comment!(code, compilation_string(lang, fname), ind_lvl = 0)
        push_code!(code, "int main() {", ind_lvl = 0)
        push_code!(code, "size_t i;")

        # Generate matrix.
        if isa(A, Matrix)
            n = LinearAlgebra.checksquare(A)
            push_code!(code, "size_t n = $n;")
            push_code!(code, "blas_type A[$(n*n)] = {")
            print_indented_matrix(lang, code, A, ind_lvl = 2)
            push_code!(code, "};")
        else # A is an integer
            n = A
            push_code!(code, "size_t n = $n;")
            push_code!(code, "srand(0);")
            push_code!(code, "blas_type *A = malloc(n*n*sizeof(*A));")
            push_code!(code, "for(i=0; i<n*n; i++){")
            push_code!(code, "A[i] = rand() / (1.0*RAND_MAX);", ind_lvl = 2)
            push_code!(code, "}")
        end

        # Call polynomial evaluation function.
        push_code!(code, "blas_type *B = malloc(n*n*sizeof(*A));")
        push_code!(code, "$blas_prefix$funname(A,n,B);")
        push_code!(code, "return 0;")
        push_code!(code, "}", ind_lvl = 0)
    end

    return code
end


text:
using ControlSystems, RobustAndOptimalControl
using Plots

"""
This is a simple SISO example which was used for debugging the implementation,
as this exact example was use in the lecture notes of the "Principles of Optimal
Control" cours of the MIT OpenCourseWare [1], where the choice of weighting
functions and dynamics gave rise to an H-infinity optimal cotnroller with a
γ of approximately 1.36, where, in our case, we get a controller at 0.93

[1] https://ocw.mit.edu/courses/aeronautics-and-astronautics/16-323-principles-of-optimal-control-spring-2008/lecture-notes/lec15.pdf

The example can be set to visualize and save plots using the two variables
  ShowPlots - true/false (true if plots are to be generated, false for testing)
  filename  - Set to string if files are to be saved, otherwise set a empty list
"""
makeplots = true

# Define the process
G   = tf([200], [0.025,1.0025,10.1,1])

# Sensitivity weight function
M, wB, A = 1.5, 10, 1e-4
WS = tf([1/M, wB],[1, wB*A])

# Output sensitivity weight function
WU = ss(0.1)

# Complementary sensitivity weight function
WT = []

# Form augmented P dynamics in state-space
P = hinfpartition(G, WS, WU, WT)

# Check that the assumptions are satisfied
flag = hinfassumptions(P)

# Synthesize the H-infinity optimal controller
flag, C, γ = hinfsynthesize(P, γrel=1)

# Extract the transfer functions defining some signals of interest
Pcl, S, CS, T = hinfsignals(P, G, C)

## Plot the specifications
if makeplots
  plot(
    specificationplot([S, CS, T], [ss(WS), WU, WT], γ),
  ## Plot the closed loop gain from w to z
    specificationplot(Pcl, γ; s_labels=["\$\\sigma(P_{cl}(j\\omega))\$"], w_labels=["\$\\gamma\$"])
  )

end


text:
"""
    Ellipsoid(radius::Number, q::U, s::U, constant::Symbol) where U<:Real

Represents an ellipsoid with axis ratios `q` and `s`, formed by deforming a sphere
of a certain `radius` at constant `:volume` or major `:axis`.
"""
struct Ellipsoid{T<:Number,U<:Real}
    radius::T
    q::U
    s::U
    constant::Symbol
end

"""
    Sphere(radius::Number)

Returns a spherical `Ellipsoid` with the specified `radius`.
"""
Sphere(radius::Number) = Ellipsoid(radius, 1, 1, :axis)

"""
    rotation_matrix_edgeon(
        p::Particles;
        radius::Union{Number,Nothing}=nothing,
        mass_weighted::Bool=true,
        inertia_matrix_type=:unweighted,
        iterative::Bool=false,
        ellipsoidal_distance::Bool=false,
        constant_volume::Bool=true,
        algorithm::Union{Symbol,Nothing}=nothing,
        return_axes::Bool=false,
    )

Returns a tuple of the 3D rotation matrix for transforming the particles edge-on using the given shape
determination algorithm, and the axis ratios `q` and `s`.
If `return_axes` is `true`, returns the same tuple with an additional element: an array of
the ellipsoid axis vectors.

See [`inertia_matrix`](@ref) and [`inertia_matrix_iterative`](@ref) for the algorithms used for
computing the inertia matrix.

# Keyword Arguments
- `radius::Union{Number,Nothing}`: radius of the initial sphere for determining the inertia matrix (`nothing` to include all particles)
- `mass_weighted::Bool`: `true` if positions should be weighted by mass
- `inertia_matrix_type`: set to `:unweighted` or `:reduced`
- `iterative::Bool`: `true` to use [`inertia_matrix_iterative`](@ref), `false` for [`intertia_matrix`](@ref)
- `ellipsoidal_distance::Bool`: `true` only considered together with inertia matrix type `:reduced`
- `constant_volume::Bool`: `true` if the sphere should be deformed into an ellipsoid with constant volue, `false` for constant major axis
- `algorithm::Union{Symbol,Nothing}`: one of `:unw`, `:red`, `:unw_i_vol`, `:unw_i_ax`, `:red_i_vol`, `:red_i_ax`, `:red_ell_i_vol`, `:red_ell_i_ax`
- `return_axes::Bool`: `true` to return tuple `(rotmat, q, s)`, `false` to return only `rotmat`
"""
function rotation_matrix_edgeon(
    p::Particles;
    radius::Union{Number,Nothing}=nothing,
    mass_weighted::Bool=true,
    inertia_matrix_type=:unweighted,
    iterative::Bool=false,
    ellipsoidal_distance::Bool=false,
    constant_volume::Bool=true,
    algorithm::Union{Symbol,Nothing}=nothing,
    return_axes::Bool=false,
)
    # set algorithm options based on algorithm
    if !isnothing(algorithm)
        inertia_matrix_type, iterative, ellipsoidal_distance, constant_volume =
            get_algorithm_variables(algorithm)
    end

    mass = mass_weighted && haskey(p, "MASS") ? p.mass : nothing

    # select particles within radius if given
    if isnothing(radius)
        𝐈 = inertia_matrix(p.pos, mass; mass_weighted, inertia_matrix_type)
    else
        if iterative
            𝐈 = inertia_matrix_iterative(
                p.pos,
                mass,
                inertia_matrix_type;
                radius,
                mass_weighted,
                constant_volume,
                ellipsoidal_distance,
            )
        else
            mask_in_radius = dropdims(sum(abs2, p.pos; dims=1) .≤ radius^2; dims=1)
            𝐈 = inertia_matrix(
                p.pos[:, mask_in_radius],
                isnothing(mass) ? mass : mass[mask_in_radius];
                mass_weighted,
                inertia_matrix_type,
            )
        end
    end

    return rotation_matrix_axis_ratios(𝐈; return_axes)
end

"""
    rotation_matrix_2D(
        p::Particles;
        radius::Union{Number,Nothing}=nothing,
        mass_weighted::Bool=true,
        inertia_matrix_type=:unweighted,
        iterative::Bool=true,
        elliptical_distance::Bool=false,
        constant_area::Bool=true,
        algorithm::Union{Symbol,Nothing}=nothing,
        perspective::Symbol=:edgeon,
        matrix_2d::Bool=false,
    )

Returns a tuple of the 3D or 2D rotation matrix for rotating the particles in 2D, so that the fitted ellipse is
aligned with the axes, and the axis ratio `q`. If `matrix_2d` is `true`, will return the 2D rotation matrix,
otherwise the 3D rotation matrix.

See [`inertia_matrix`](@ref) and [`inertia_matrix_iterative`](@ref) for the algorithms used for
computing the inertia matrix (here only used in 2D, note that terminology is different than in 3D, i.e. elliptical instead of ellipsoidal and constant area instead of volume).

# Keyword Arguments
- `radius::Union{Number,Nothing}`: radius of the initial circle for determining the inertia matrix (`nothing` to include all particles)
- `mass_weighted::Bool`: `true` if positions should be weighted by mass
- `inertia_matrix_type`: set to `:unweighted` or `:reduced`
- `iterative::Bool`: `true` to use [`inertia_matrix_iterative`](@ref), `false` for [`intertia_matrix`](@ref)
- `elliptical_distance::Bool`: `true` only considered together with inertia matrix type `:reduced`
- `constant_area::Bool`: `true` if the circle should be deformed into an ellipse with constant area, `false` for constant major axis
- `algorithm::Union{Symbol,Nothing}`: one of `:unw`, `:red`, `:unw_i_vol`, `:unw_i_ax`, `:red_i_vol`, `:red_i_ax`, `:red_ell_i_vol`, `:red_ell_i_ax` (note the same use of `vol` like in the 3D case instead of `area`)
- `perspective::Symbol`: one of `:edgeon`, `:sideon`, `:faceon`
- `matrix_2d::Bool`: `true` to return 2D rotation matrix, `false` to return only 3D rotation matrix
"""
function rotation_matrix_2D(
    p::Particles;
    radius::Union{Number,Nothing}=nothing,
    mass_weighted::Bool=true,
    inertia_matrix_type=:unweighted,
    iterative::Bool=true,
    elliptical_distance::Bool=false,
    constant_area::Bool=true,
    algorithm::Union{Symbol,Nothing}=nothing,
    perspective::Symbol=:edgeon,
    matrix_2d::Bool=false,
)

    # set algorithm options based on algorithm
    if !isnothing(algorithm)
        inertia_matrix_type, iterative, elliptical_distance, constant_area =
            get_algorithm_variables(algorithm)
    end

    dims = get_dims(perspective)

    # select stars within radius if given
    mass = haskey(p, "MASS") ? p.mass : nothing
    if isnothing(radius)
        𝐈 = @views inertia_matrix(
            p.pos[dims, :],
            mass;
            mass_weighted,
            inertia_matrix_type,
        )
    else
        if iterative
            𝐈 = @views inertia_matrix_iterative(
                p.pos[dims, :],
                mass,
                inertia_matrix_type;
                radius,
                mass_weighted,
                constant_volume=constant_area,
                ellipsoidal_distance=elliptical_distance,
            )
        else
            mask_in_radius =
                @views dropdims(sum(abs2, p.pos[dims, :]; dims=1) .≤ radius^2; dims=1)
            𝐈 = inertia_matrix(
                p.pos[dims, mask_in_radius],
                isnothing(mass) ? mass : mass[mask_in_radius];
                mass_weighted,
                inertia_matrix_type,
            )
        end
    end

    Q⁻¹, q, s = rotation_matrix_axis_ratios(𝐈)

    # convert 2d rotation matrix to 3d
    if !matrix_2d
        Q⁻¹ = rotation_matrix_to_3d(Q⁻¹, dims)
    end
    return Q⁻¹, q
end

"""
    inertia_matrix(
        pos::AbstractMatrix{<:Number},
        mass::Union{AbstractVector{<:Number},Nothing};
        mass_weighted::Bool=true,
        inertia_matrix_type::Symbol=:unweighted,
    )

Create inertia matrix of different types.

# Keyword Arguments
- `mass_weighted::Bool`: `true` if positions should be weighted by mass
- `inertia_matrix_type::Symbol`: set to `:unweighted` or `:reduced`

## Inertia matrix types
- Unweighted:               Eq. 1 in [Joachimi et al. 2013](https://arxiv.org/abs/1203.6833),
- Mass-weighted unweighted: Eq. 2 in [Tenneti et al. 2015](https://arxiv.org/abs/1409.7297),
- Reduced:                  Eq. 3 in [Tenneti et al. 2015](https://arxiv.org/abs/1409.7297),
                            but without the mass weighting,
- Mass-weighted reduced:    Eq. 3 in [Tenneti et al. 2015](https://arxiv.org/abs/1409.7297).

Note that in the mass-weighted cases, it is not normalized to the total mass,
i.e. different than in Eq2/Eq3 of T2015.
"""
function inertia_matrix(
    pos::AbstractMatrix{<:Number},
    mass::Union{AbstractVector{<:Number},Nothing};
    mass_weighted::Bool=true,
    inertia_matrix_type::Symbol=:unweighted,
)
    # check dimensions
    @assert 2 ≤ size(pos, 1) ≤ 3
    if mass_weighted
        @assert !isnothing(mass)
        @assert length(mass) == size(pos, 2)
    end

    dim = size(pos, 1)

    # create empty matrix of correct float type (strips units)
    T = mass_weighted ? promote_type(typeof(ustrip(pos[1])), typeof(ustrip(mass[1]))) : typeof(ustrip(pos[1]))
    𝐈_up = zeros(T, dim, dim)

    fill_𝐈_up!(𝐈_up, pos, mass, mass_weighted, inertia_matrix_type)

    # produce symmetric matrix from upper triangle
    return Symmetric(𝐈_up, :U)
end

"""
    inertia_matrix_iterative(
        pos::AbstractMatrix{<:Number},
        mass::Union{AbstractVector{<:Number},Nothing},
        inertia_matrix_type::Symbol;
        radius::Union{Number,Nothing}=nothing,
        mass_weighted::Bool=true,
        ellipsoidal_distance::Bool=false,
        constant_volume::Bool=true,
        ϵ::Real=1e-4,
    )

Create inertia matrix of different types iteratively (generally more stable than [`inertia_matrix`](@ref)).

# Keyword Arguments
- `radius::Union{Number,Nothing}`: radius of the sphere to start the iteration with
- `mass_weighted::Bool`: `true` if positions should be weighted by mass
- `ellipsoidal_distance::Bool`: only in combination with `reduced` inertia matrix type; `true` to use distances corrected by ellipsoidal axes
- `constant_volume::Bool`: `true` if the sphere should be deformed into an ellipsoid with constant volume, `false` for constant major axis
- `ϵ::Real`: iterate until a relative change of the shape axis ratios less than `ϵ` is encountered

## Inertia matrix types
- Unweighted:               Eq. 1 in [Joachimi et al. 2013](https://arxiv.org/abs/1203.6833),
- Mass-weighted unweighted: Eq. 2 in [Tenneti et al. 2015](https://arxiv.org/abs/1409.7297),
- Reduced:                  Eq. 3 in [Tenneti et al. 2015](https://arxiv.org/abs/1409.7297),
                            but without the mass weighting,
- Mass-weighted reduced:    Eq. 3 in [Tenneti et al. 2015](https://arxiv.org/abs/1409.7297).
- Reduced with ellipsoidal distances:               TODO
- Mass-weighted reduced with ellipsoidal distances: TODO

Note that in the mass-weighted cases, it is not normalized to the total mass,
i.e. different than in Eq2/Eq3 of T2015.
"""
function inertia_matrix_iterative(
    pos::AbstractMatrix{<:Number},
    mass::Union{AbstractVector{<:Number},Nothing},
    inertia_matrix_type::Symbol;
    radius::Union{Number,Nothing}=nothing,
    mass_weighted::Bool=true,
    ellipsoidal_distance::Bool=false,
    constant_volume::Bool=true,
    ϵ::Real=1e-4,
)
    # check dimensions
    @assert 2 ≤ size(pos, 1) ≤ 3
    if mass_weighted
        @assert !isnothing(mass)
        @assert length(mass) == size(pos, 2)
    end

    dim = size(pos, 1)

    # create empty matrix of correct float type (strips units)
    T = mass_weighted ? promote_type(typeof(ustrip(pos[1])), typeof(ustrip(mass[1]))) : typeof(ustrip(pos[1]))
    𝐈_up = zeros(T, dim, dim)

    # initialize with spherical reduced inertia matrix
    q_old = s_old = one(T)
    if isnothing(radius)
        fill_𝐈_up!(𝐈_up, pos, mass, mass_weighted, inertia_matrix_type)
    else
        mask = dropdims(sum(abs2, pos; dims=1) .≤ radius^2; dims=1)
        fill_𝐈_up!(
            𝐈_up,
            pos[:, mask],
            isnothing(mass) ? mass : mass[mask],
            mass_weighted,
            inertia_matrix_type,
        )
    end
    𝐈 = Symmetric(𝐈_up, :U)
    Q⁻¹, q, s = rotation_matrix_axis_ratios(𝐈)

    # recompute inertia matrix by adapting ellipsoidal shape until convergence
    i = 0
    while (!isapprox(q_old, q; rtol=ϵ) || !(isnothing(s) || isapprox(s_old, s; rtol=ϵ))) && i < 50
        # ellipsoidal distance for masking
        r²_ell = r²_ellipsoid(rotate(pos, Q⁻¹), q, s)

        # ellipsoidal vs regular distance for inertia matrix
        if ellipsoidal_distance
            r² = r²_ell
        else
            r² = dropdims(sum(abs2, pos; dims=1); dims=1)
        end
        if isnothing(radius)
            fill_𝐈_up!(𝐈_up, pos, mass, mass_weighted, inertia_matrix_type, r²=r²)
        else
            factor = dim == 2 ? sqrt(q) : cbrt(q * s)
            r²_max = constant_volume ? (radius / factor)^2 : radius^2
            mask = reshape(r²_ell .≤ r²_max, :)
            fill_𝐈_up!(
                𝐈_up,
                pos[:, mask],
                isnothing(mass) ? mass : mass[mask],
                mass_weighted,
                inertia_matrix_type;
                r²=r²[mask],
            )
        end
        𝐈 = Symmetric(𝐈_up, :U)

        q_old, s_old = q, s
        Q⁻¹, q, s = rotation_matrix_axis_ratios(𝐈)

        i += 1
    end

    return 𝐈
end

"""
    r²_circle(pos::AbstractMatrix{<:Number})

Squared distances in 2D from origin.
"""
r²_circle(pos::AbstractMatrix{<:Number}) = @views @. pos[1, :]^2 + pos[2, :]^2

"""
    r²_sphere(pos::AbstractMatrix{<:Number})

Squared distances from origin in 2D (calls [`r²_circle`](@ref) in those cases) and 3D.
"""
function r²_sphere(pos::AbstractMatrix{<:Number})
    if size(pos, 1) == 2
        return r²_circle(pos)
    else
        return @views @. pos[1, :]^2 + pos[2, :]^2 + pos[3, :]^2
    end
end

"""
    r²_ellipsoid(pos::AbstractMatrix{<:Number}, q::Real, s::Nothing)

Squared ellipsoidal distances in 2D from origin with axis ratio `q`.

For 3D version, see Eq. 6 in Allgood et al. 2006, https://arxiv.org/abs/astro-ph/0508497
"""
function r²_ellipsoid(pos::AbstractMatrix{<:Number}, q::Real, s::Nothing)
    return @views @. pos[1, :]^2 + pos[2, :]^2 / q^2
end

"""
    r²_ellipsoid(pos::AbstractMatrix{<:Number}, q::Real, s::Real)

Squared ellipsoidal distances from origin with axis ratios `q` and `s`.

See Eq. 6 in Allgood et al. 2006, https://arxiv.org/abs/astro-ph/0508497
"""
function r²_ellipsoid(pos::AbstractMatrix{<:Number}, q::Real, s::Real)
    return @views @. pos[1, :]^2 + pos[2, :]^2 / q^2 + pos[3, :]^2 / s^2
end

"""
    ellipsoidal_mask(
        pos::AbstractMatrix{<:Number},
        ellipsoid::Ellipsoid;
        return_r²_ellipsoid::Bool=false,
    )

Create a `BitArray` for masking the particles within an `ellipsoid`.

The ellipsoid is centered around the origin of the coordinate system with its axes placed
on the x, y, and z axes, from smallest to largest.
If `return_r²_ellipsoid` is `true`, returns a tuple of the mask and the vector r² for all particles.
"""
function ellipsoidal_mask(
    pos::AbstractMatrix{<:Number},
    ellipsoid::Ellipsoid;
    return_r²_ellipsoid::Bool=false,
)
    r² = r²_ellipsoid(pos, ellipsoid.q, ellipsoid.s)
    r²_max = if ellipsoid.constant === :volume
        (ellipsoid.radius / cbrt(ellipsoid.q * ellipsoid.s))^2
    else
        ellipsoid.radius^2
    end

    mask = reshape(r² .≤ r²_max, :)

    return return_r²_ellipsoid ? (mask, r²) : mask
end

"""
    fill_𝐈_up!(
        𝐈_up::AbstractMatrix{<:Real},
        pos::AbstractMatrix{<:Number},
        mass::Union{AbstractVector{<:Number},Nothing},
        mass_weighted::Bool,
        inertia_matrix_type::Symbol;
        r²::Union{AbstractVector{<:Number},Nothing}=nothing,
    )

Fill upper triangle of the inertia matrix for [`inertia_matrix`](@ref) and [`inertia_matrix_iterative`](@ref)

# Arguments
- `mass_weighted::Bool`: `true` if positions should be weighted by mass
- `inertia_matrix_type::Symbol`: set to `:unweighted` or `:reduced`

# Keyword Arguments
- `r²::Union{AbstractVector{<:Number},Nothing}`: if precomputed, the absolute squared distances from the coordinate system's origin; only used when `:reduced` is used

Note that in the mass-weighted cases, it is not normalized to the total mass,
i.e. different than in Eq2/Eq3 of T2015.
"""
function fill_𝐈_up!(
    𝐈_up::AbstractMatrix{<:Real},
    pos::AbstractMatrix{<:Number},
    mass::Union{AbstractVector{<:Number},Nothing},
    mass_weighted::Bool,
    inertia_matrix_type::Symbol;
    r²::Union{AbstractVector{<:Number},Nothing}=nothing,
)
    # fill upper triangle
    for i in 1:size(𝐈_up, 1), j in i:size(𝐈_up, 1)
        if inertia_matrix_type === :unweighted
            if mass_weighted
                𝐈_up[i, j] = @views sum(mass .* pos[i, :] .* pos[j, :]) |> ustrip
            else
                𝐈_up[i, j] = @views sum(pos[i, :] .* pos[j, :]) |> ustrip
            end
        elseif inertia_matrix_type === :reduced
            r² = isnothing(r²) ? dropdims(sum(abs2, pos; dims=1); dims=1) : r²
            if mass_weighted
                𝐈_up[i, j] = @views sum(skipnan(mass .* pos[i, :] .* pos[j, :] ./ r²)) |> ustrip
            else
                𝐈_up[i, j] = @views sum(skipnan(pos[i, :] .* pos[j, :] ./ r²)) |> ustrip
            end
        end
    end

    return 𝐈_up
end

"""
    rotation_matrix_axis_ratios(𝐈::AbstractMatrix{<:Real})

Returns rotation matrix `Q⁻¹` and axis ratios `q` and `s` for an inertia matrix `𝐈`.
`s` is `nothing` if `𝐈` is a 2D matrix.
If `return_axes` is `true`, returns the same tuple with an additional element: an array of
the ellipsoid axis vectors.
"""
function rotation_matrix_axis_ratios(𝐈::AbstractMatrix{<:Real}; return_axes::Bool=false)
    # compute normalized and sorted (in descending order) eigenvalues and eigenvectors
    λ, Q = eigen(𝐈, sortby=λ -> -λ) # λ[k] -> Q[:,k]
    Q⁻¹ = inv(Q)

    # make rotation matrix always rotate to the same symmetry
    Q⁻¹ .*= replace(sign.(Q⁻¹[:, 1]), 0 => one(eltype(Q⁻¹)), -0.0 => one(eltype(Q⁻¹)))

    # determine axis ratios
    q = sqrt(λ[2] / λ[1])
    s = size(Q, 1) == 3 ? sqrt(λ[3] / λ[1]) : nothing

    if return_axes
        T = eltype(Q)

        # make vectors always point in same general direction
        Q .*= replace(sign.(Q[1, :]'), zero(T) => one(T), -zero(T) => one(T))

        return Q⁻¹, q, s, collect(eachcol(Q))
    end

    return Q⁻¹, q, s
end

"""
    triaxiality(q::Real, s::Real)

Triaxiality from the axis ratios `q` and `s`.

-   ``0 < T < 1/3``: oblate
- ``1/3 < T < 2/3``: triaxial
- ``2/3 < T < 1``:   prolate
"""
triaxiality(q::Real, s::Real) = (1 - q^2) / (1 - s^2)

@doc raw"""
    ellipticity(q::Real)

Ellipticity from axis ratio `q`: ``ϵ = 1 - q``
"""
ellipticity(q::Real) = 1 - q

@doc raw"""
    eccentricity(q::Real)

Eccentricity from axis ratio `q`: ``e = \sqrt{1 - q^2}``
"""
eccentricity(q::Real) = sqrt(1 - q^2)

"""
    rotation_matrix_to_3d(R::AbstractMatrix{T}, dims::AbstractVector{<:Integer}) where T<:Real

Converts a 2D rotation matrix `R` in the plane given by dimensions `dims` (e.g. `[1, 2]` or `[3, 1]`)
to a 3D rotation matrix.
"""
function rotation_matrix_to_3d(R::AbstractMatrix{T}, dims::AbstractVector{<:Integer}) where {T<:Real}
    @assert size(R) == (2, 2)
    @assert length(dims) == 2

    R3d = zeros(T, 3, 3)

    for i in 1:2, j in 1:2
        R3d[dims[i], dims[j]] = R[i, j]
    end

    R3d[(6 - sum(dims)), (6 - sum(dims))] = one(T)

    return R3d
end

function get_algorithm_variables(algorithm::Symbol)
    if algorithm === :unw
        :unweighted, false, false, false
    elseif algorithm === :red
        :reduced, false, false, false
    elseif algorithm === :unw_i_vol
        :unweighted, true, false, true
    elseif algorithm === :unw_i_ax
        :unweighted, true, false, false
    elseif algorithm === :red_i_vol
        :reduced, true, false, true
    elseif algorithm === :red_i_ax
        :reduced, true, false, false
    elseif algorithm === :red_ell_i_vol
        :reduced, true, true, true
    elseif algorithm === :red_ell_i_ax
        :reduced, true, true, false
    else
        error("""algorithm is none of the allowed values:
                  :unw, :red, :unw_i_vol, :unw_i_ax :red_i_vol, :red_i_ax,
                  :red_ell_i_vol, :red_ell_i_ax.""")
    end
end


text:
"""
   set up damping profile for PML boundary conditions, n is the model size include pml,
npml is the number of PML layers, omega is radian frequency.
"""
function get_gamma_profile(n::Ti, npml::Ti, h; apml=900.0, oneside_flag=false) where {Ti<:Int64}

    # constants
    pi_over_2 = 0.5 * pi   # pi/2

    # damping profile of PML
    gamma_integer = zeros(Float64, n)
    gamma_stagger = zeros(Float64, n)

    # PML width
    pml_width = npml * h

    # model width include PML layers
    model_width = (n - 1) * h

    # both side of gamma_integer and one side of gamma_stagger
    for i = 1 : npml

        # current coordinate of x
        x_integer = (i-1) * h
        x_stagger = (i-1+0.5) * h

        # distance of current grid to the edge of computation area
        dis_integer = pml_width - x_integer
        dis_stagger = pml_width - x_stagger

        # pml coefficients
        gamma_integer[i] = apml * (1.0 - cos(dis_integer / pml_width * pi_over_2))
        gamma_stagger[i] = apml * (1.0 - cos(dis_stagger / pml_width * pi_over_2))

        # symmetrical gamma_integer
        gamma_integer[n-i+1] = gamma_integer[i]
    end

    # one side PML obsorbing boundary
    if oneside_flag
       for i = 1 : npml
           gamma_integer[i] = 0.0
           gamma_stagger[i] = 0.0
       end
    end

    # the right side of gamma_stagger
    for i = 1 : npml + 1

        # the coordinate of the stagger grid (from right to left)
        x_stagger = model_width + 0.5 * h - (i-1)*h

        # distance of current grid to the edge of computation area
        dis_stagger = x_stagger - (model_width - pml_width)

        # pml coefficients
        gamma_stagger[n-i+1] = apml * (1.0 - cos(dis_stagger / pml_width * pi_over_2))
    end

    return gamma_integer, gamma_stagger
end

# """
#    struct for the coefficients of the split PML boundary
# """
# struct PMLCoefficients{Tv<:AbstractFloat, Ti<:Int64}
#     vz::SparseMatrixCSC{Tv, Ti}
#     vx::SparseMatrixCSC{Tv, Ti}
#     pz::SparseMatrixCSC{Tv, Ti}
#     px::SparseMatrixCSC{Tv, Ti}
# end
#
# """
#    constructor for PML obsorbing boundary coefficients
# """
# function PMLCoefficients(params::TdParams)
#
#     data_format = params.data_format
#
#     # coefficient depends on number of layers
#     R = data_format(0.01)
#     if params.npml == 10
#        R = data_format(0.001)
#     elseif params.npml >= 20
#        R = data_format(0.0001)
#     end
#
#     ll = data_format(2.)
#     vmin = minimum(params.vel)
#     alpha_z = data_format(-1.5) * vmin / (params.npml * params.dz) * log(R)
#     alpha_x = data_format(-1.5) * vmin / (params.npml * params.dx) * log(R)
#
#     # allocate memory for PML coefficient
#     vz = zeros(data_format, params.Nz, params.Nx)
#     vx = zeros(data_format, params.Nz, params.Nx)
#     pz = zeros(data_format, params.Nz, params.Nx)
#     px = zeros(data_format, params.Nz, params.Nx)
#
#     # expand model to accomdate PML layers
#     rho = model_padding(params.rho, params.npml, params.free_surface)
#
#     if params.free_surface
#        #  lower
#        for iz = params.nz+1 : params.Nz
#            vz[iz,:] = alpha_z * ((iz+0.5-params.nz) / params.npml)^ll * rho[iz,:]
#            pz[iz,:] = alpha_z * ((iz    -params.nz) / params.npml)^ll * ones(data_format, params.Nx)
#        end
#
#        #  left
#        for ix = 1 : params.npml
#            vx[:,ix] = alpha_x * ((params.npml+0.5-ix) / params.npml)^ll * rho[:,ix]
#            px[:,ix] = alpha_x * ((params.npml+1  -ix) / params.npml)^ll * ones(data_format, params.Nz)
#        end
#
#        #  right
#        for ix = params.nx + params.npml + 1 : params.Nx
#            vx[:,ix] = alpha_x * ((ix+0.5-params.nx-params.npml) / params.npml)^ll * rho[:,ix]
#            px[:,ix] = alpha_x * ((ix    -params.nx-params.npml) / params.npml)^ll * ones(data_format, params.Nz)
#        end
#     else
#        #  upper and left
#        for i = 1 : params.npml
#            vz[i,:] = alpha_z * ((params.npml+0.5-i) / params.npml)^ll * rho[i,:]
#            vx[:,i] = alpha_x * ((params.npml+0.5-i) / params.npml)^ll * rho[:,i]
#            pz[i,:] = alpha_z * ((params.npml+1  -i) / params.npml)^ll * ones(data_format, params.Nx)
#            px[:,i] = alpha_x * ((params.npml+1  -i) / params.npml)^ll * ones(data_format, params.Nz)
#        end
#
#        #  lower
#        for iz = params.nz + params.npml + 1 : params.Nz
#            vz[iz,:] = alpha_z * ((iz+0.5-params.nz-params.npml) / params.npml)^ll * rho[iz,:]
#            pz[iz,:] = alpha_z * ((iz    -params.nz-params.npml) / params.npml)^ll * ones(data_format, params.Nx)
#        end
#
#        #  right
#        for ix = params.nx + params.npml + 1 : params.Nx
#            vx[:,ix] = alpha_x * ((ix+0.5-params.nx-params.npml) / params.npml)^ll * rho[:,ix]
#            px[:,ix] = alpha_x * ((ix    -params.nx-params.npml) / params.npml)^ll * ones(data_format, params.Nz)
#        end
#     end
#
#     return PMLCoefficients(sparse(vz), sparse(vx), sparse(pz), sparse(px))
# end


text:
@testset "select" begin
    actual = Skans.remove_extra_href_info("href" => "www.bbc.com?bar")
    @test last(actual) == "www.bbc.com?[...]"

    actual = Skans.remove_extra_href_info("href" => "#header")
    @test last(actual) == "#[...]"

    html = """
        <!doctype html>

            <meta charset="UTF-8">
            <meta name="viewport" content="width=device-width, initial-scale=1">

            <link rel="stylesheet" href="/css/foo.css">
            <link rel="icon" href="/assets/favicon.png">
            <title> Home - Site.xyz</title>

            <meta property="og:title" content="Home" />

            <meta property="og:description" content="description" />

            <meta name="twitter:title" content="Home" />

NAME

About

        lorem ipsum

        """

    doc = parsehtml(html)
    # Smoke test.
    map!(identity, doc)

    expected = """

        <!DOCTYPE html>

            <meta/>
            <meta/>
            <link/>
            <link/>
            <title>
              Home - Site.xyz
            </title>
            <meta content="Home" property="og:title"/>
            <meta content="description" property="og:description"/>

            <meta content="Home" name="twitter:title"/>

                  NAME

                      About

              lorem ipsum

        """

    actual = Skans.clean(html)

    println("expected:")
    print(expected)
    println()
    println("actual:")
    print(actual)
    println()

    @test strip(expected) == strip(actual)
end