text:
# EXTRACT RELATIONS
# direct succession     >
# causality             ->
# parallel              ||
# choice                #

# FOOTPRINT MATRIX

# OUTLINE
# 1. Each activity corresponds to a transition in α(L)
# 2. Fix the set of start activities (first elements of each trace)
# 3. Fix the set of end activities (last elements of each trace)
# 4. Calculate pairs (A, B)
# 5. Delete non-maximal pairs (A, B)
# 6. Determine places p_(A, B) from pairs (A, B)
# 7. Determine the flow relation
# 8. Put everything together

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

# MAYBE NOT NEEDED
# mutable struct Footprint{T}
#     colnames::Array{T}
#     rownames::Array{T}
#     relation::Matrix{Symbol}
# end
# ------------------------

using ProcessMining
using Pipe

# 1) GET RAW TRACES AND ACTIVITY SET FROM THE LOG
function extract_event_traces(eventlog::EventLog)
    return [
        [event.name for event in trace.events] for trace in eventlog.traces
    ]
end

function extract_activity_set(extracted_event_traces::AbstractArray)
    activity_set = extracted_event_traces |>
        Iterators.flatten |>
        collect |>
        unique
    return activity_set
end
# ------------------------

# 2 + 3) SIMPLE STEPS: START AND END ACTIVITIES
function extract_start_activities(extracted_event_traces::AbstractArray)
    return unique([first(trace) for trace in extracted_event_traces])
end

function extract_end_activities(extracted_event_traces::AbstractArray)
    return unique([last(trace) for trace in extracted_event_traces])
end
# ------------------------

# 4a) EXTRACT RELATIONS FROM THE ACTIVITY LOG
function extract_direct_succession_relation(extracted_event_traces)
    direct_succession_relation = Set(Tuple{String, String}[])
    for trace in extracted_event_traces
        for i in 1:length(trace)
            if !(i == length(trace))
                push!(direct_succession_relation, (trace[i], trace[i + 1]))
            end
        end
    end
    return direct_succession_relation
end

function extract_causality_relation(direct_succession_relation)
    causality_relation = Set(Tuple{String, String}[])
    for tup in direct_succession_relation
        if !((tup[2], tup[1]) in direct_succession_relation)
            push!(causality_relation, tup)
        end
    end
    return causality_relation
end

function extract_parallel_relation(direct_succession_relation)
    parallel_relation = Set(Tuple{String, String}[])
    for tup in direct_succession_relation
        if !(tup in parallel_relation) && ((tup[2], tup[1]) in direct_succession_relation)
            push!(parallel_relation, tup)
            push!(parallel_relation, (tup[2], tup[1]))
        end
    end
    return parallel_relation
end

function extract_choice_relation(direct_succession_relation, activity_set)
    all_pairs = Set([(a, b) for a in activity_set for b in activity_set])
    choice_relation = setdiff(all_pairs, direct_succession_relation)
    return choice_relation
end
# ------------------------

# 4b) GET SET X_L
function iteratively_build_x(causality_relation, choice_relation)
    place_pairs = get_initial_place_pairs(causality_relation)
    while true
        size = length(place_pairs)
        build_place_pairs!(place_pairs, choice_relation, "forward")
        build_place_pairs!(place_pairs, choice_relation, "backward")
        new_size = length(place_pairs)
        if size == new_size
            break
        end
    end
    return place_pairs
end

function get_initial_place_pairs(causality_relation)
    place_pairs = Set([])
    for mapping in causality_relation
        push!(place_pairs, (Set([mapping[1]]), Set([mapping[2]])))
    end
    return place_pairs
end

function build_place_pairs!(place_pairs, choice_relation, mode::String)
    if mode == "forward"
        a = 1
        b = 2
    elseif mode == "backward"
        a = 2
        b = 1
    end
    place_pairs_addition = Set([])
    for p1 in place_pairs
        for p2 in place_pairs
            if p1 != p2 && (p1[a] == p2[a])
                combinations = get_combinations(p1, p2, b)
                flag = all_combinations_unrelated(combinations, choice_relation)
                if flag && (mode == "forward")
                    push!(place_pairs_addition, (p1[a], union(p1[b], p2[b])))
                elseif flag && (mode == "backward")
                    push!(place_pairs_addition, (union(p1[b], p2[b]), p1[a]))
                end
            end
        end
    end
    union!(place_pairs, place_pairs_addition)
    return place_pairs
end

function get_combinations(p1, p2, b)
    combinations = []
    for elem_p1 in collect(p1[b])
        for elem_p2 in collect(p2[b])
            push!(combinations, (elem_p1, elem_p2))
        end
    end
    return combinations
end

function all_combinations_unrelated(combinations, choice_relation)
    for comb in combinations
        if !(comb in choice_relation)
            return false
        end
    end
    return true
end
# ------------------------

# --------------------------- TESTING --------------------------- #
begin
    EVENT_LOG_FILE = "road_fines.xes"  # alternative: Performance.xes
    eventlog = read_xes(joinpath("data", EVENT_LOG_FILE))
    extracted_event_traces = extract_event_traces(eventlog)
    start_activities = extract_start_activities(extracted_event_traces)
    end_activities = extract_end_activities(extracted_event_traces)
    activity_set = extract_activity_set(extracted_event_traces)
    direct_succession_relation = extract_direct_succession_relation(extracted_event_traces)
    causality_relation = extract_causality_relation(direct_succession_relation)
    parallel_relation = extract_parallel_relation(direct_succession_relation)
    choice_relation = extract_choice_relation(direct_succession_relation, activity_set)
    X = iteratively_build_x(causality_relation, choice_relation)
end


text:
using Pkg
Pkg.activate("..")
using TSML
using TSML.TSMLTransformers
using DataFrames
using Dates
using CSV
using Plots

fname ="testdata.csv"
dat = CSV.read(fname)
rename!(dat,names(dat)[1]=>:Date,names(dat)[2]=>:Value)
dat[:Date] = DateTime.(dat[:Date],"d/m/y H:M")
orig = deepcopy(dat)
filter1 = DateValgator()
filter2 = DateValNNer(Dict(:nnsize=>1))

fit!(filter1,dat,[])
res1=transform!(filter1,dat)

fit!(filter2,res1,[])
res2=transform!(filter2,res1)

mypipeline = Pipeline(Dict(
		:transformers => [csvreader,filter1,filter2]
	)
)

fit!(mypipeline)
res = transform!(mypipeline)

Plots.plot(res[:Value][end-3000:end])

rfname = replace(fname,".csv"=>"-result.csv")
res |> CSV.write(rfname)

using TSML.TSMLTypes
import TSML.TSMLTypes.fit!
import TSML.TSMLTypes.transform!

mutable struct CSVDateValReader <: Transformer
	model
	args
	function CSVDateValReader(args=Dict())
		default_args = Dict(
			:filename => "",
			:dateformat => ""
		)
		new(nothing,mergedict(default_args,args))
	end
end

function fit!(csvrdr::CSVDateValReader,x::T=[],y::Vector=[]) where {T<:Union{DataFrame,Vector,Matrix}}
	fname = csvrdr.args[:filename]
	fmt = csvrdr.args[:dateformat]
	(fname != "" && fmt != "") || error("missing filename or date format")
	model = csvrdr.args
end

function transform!(csvrdr::CSVDateValReader,x::T=[]) where {T<:Union{DataFrame,Vector,Matrix}}
	fname = csvrdr.args[:filename]
	fmt = csvrdr.args[:dateformat]
	df = CSV.read(fname)
	ncol(df) == 2 || error("dataframe should have only two columns: Date,Value")
	rename!(df,names(df)[1]=>:Date,names(df)[2]=>:Value)
	df[:Date] = DateTime.(df[:Date],fmt)
	df
end

csvreader = CSVDateValReader(Dict(:filename=>"testdata.csv",:dateformat=>"d/m/y H:M"))

fit!(csvreader)
transform!(csvreader)

ext: jl

max_stars_repo_name: ppalmes/TSMLextra

max_issues_repo_name: ppalmes/TSMLextra

max_forks_repo_name: ppalmes/TSMLextra


text:
function solveRadau(rr::RadauIntegrator{T_object,NX,NC,NR,NSM}, x0::Vector{Float64}, t::Float64=0.0) where {T_object,NX,NC,NR,NSM}
    table = rr.table[rr.rule.s]
    calcJacobian!(rr, x0, t)
    return solveRadau_inner(rr, x0, table, t)
end

function solveRadau_inner(rr::RadauIntegrator{T_object,NX,NC,NR,NSM}, x0::Vector{Float64},
        table::RadauTable{NS}, t::Float64=0.0) where {T_object,NX,NC,NR,NSM,NS}

    simple_newton!(rr, x0, table, t)
    # update_x_err_norm!(rr, table, x0)
    calc_and_update_h!(rr, table, x0)
    update_rule!(rr)
    # h_new = calc_h_new(rr, table, x0, is_converge)
    # update_h!(rr, h_new)  # TODO: move this and the above line into one function
    # update_rule!(rr, is_converge)
    if is_exit_flag_success(rr)
        x_final = get_X_final(rr, table)
        h_taken = rr.step.hᵏ⁻¹
        return rr.step.hᵏ⁻¹, x_final, t + h_taken  # step actually taken, x_final
    else
        if rr.step.h < rr.step.h_min
            error("time step is too small, something is wrong")
        else
            table = get_table_from_current_s(rr)  # need to get new table because s may have changed
            return solveRadau_inner(rr, x0, table)
        end
    end
end

# function print_exit_flag(rr::RadauIntegrator)
#     exit_flag = rr.step.exit_flag
#     (exit_flag == -9999) && println("never run")
#     (exit_flag == 0) && println("step successfull")
#     (exit_flag == 1) && println("iteration limit exceeded")
#     (exit_flag == 2) && println("tolerance failure")
#     (exit_flag == 3) && println("divergence")
# end
#
# @inline set_exit_flag_success(rr::RadauIntegrator) = (rr.step.exit_flag = 0)
# @inline set_exit_flag_fail_iter(rr::RadauIntegrator) = (rr.step.exit_flag = 1)
# @inline set_exit_flag_fail_tol(rr::RadauIntegrator) = (rr.step.exit_flag = 2)
# @inline set_exit_flag_fail_diverge(rr::RadauIntegrator) = (rr.step.exit_flag = 3)

function simple_newton!(rr::RadauIntegrator{T_object,NX,NC,NR,NSM}, x0::Vector{Float64}, table::RadauTable{NS},
        t::Float64=0.0) where {T_object,NX,NC,NR,NSM,NS}

    # TODO: get rid of residual__ and replace with rr.rule.θ

    updateInvC!(rr, table)
    initialize_X_with_X0!(rr, table, x0)

    residual_prev = Inf
    k_iter_max = rr.rule.k_iter_max
    res_vec = SVector{3,Float64}(Inf,Inf,Inf)
    for k_iter = 1:k_iter_max
        rr.rule.k_iter = k_iter
        zeroFill!(rr.ct.Ew_stage, table)
        zeroFill!(rr.ct.delta_Z_stage, table)
        updateFX!(rr, table, x0, t)
        residual__ = calcEw!(rr, table, x0)
        if updateStageX!(rr, table)
            set_exit_flag_fail_diverge(rr)
            update_dense_unsuccessful!(rr)
            return nothing
        end

        if residual__ < rr.step.tol_newton
            set_exit_flag_success(rr)
            update_dense_successful!(rr, table)
            update_x_err_norm!(rr, table, x0)
            return nothing
            # return true
        end

        if k_iter != 1
            rr.rule.θᵏ⁻¹ = rr.rule.θ
            rr.rule.θ = sqrt(residual__)
            rr.rule.Ψ_k = sqrt(rr.rule.θᵏ⁻¹ * rr.rule.θ)
        else
            rr.rule.θ = sqrt(residual__)
            rr.rule.Ψ_k = rr.rule.θ
        end
        residual_prev = residual__  # update residual
        res_vec = SVector{3,Float64}(residual__, res_vec[1], res_vec[2])
        if res_vec[3] < res_vec[2] < res_vec[1]  # two consequitive issues
            set_exit_flag_fail_diverge(rr)
            update_dense_unsuccessful!(rr)
            return nothing
            # return false
        end
    end
    set_exit_flag_fail_iter(rr)
    update_dense_unsuccessful!(rr)
    return nothing
    # return false  # exceeded iteration limit
end


text:
"""
   `interp_chunk_to_grid_linear!( flux_out, var_out, chunk_of_spectrum, wavelengths )`
Return spectra interpolated onto a grid of points using linear interpolation.
# Arguments:
- flux_out: (results stored into this array)
- var_out: (results stored into this array)
- chunk_of_spectrum
- wavelengths: AbstractRange or AbstractArray of locations where chunk is to be interpolated to
- boost_factor: divide wavelengths by boost_factor
Optional Arguments:
- Filter: Vector with pre-allocated workspace (if length>=1)
# Returns
- flux_out
"""
function interp_chunk_to_shifted_grid_linear!( flux_out::AbstractArray{T1,1}, var_out::AbstractArray{T2,1}, chunk::AC, grid::AR, boost_factor::Real ) where {
    T1<:Real, T2<:Real, AC<:AbstractChunkOfSpectrum, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    @assert size(flux_out) == size(var_out)
    @assert size(flux_out) == size(grid)
    lin_interp_flux = #LinearInterpolation.extrapolate(
        LinearInterpolation.make_interpolator_linear_flux(chunk) #, Flat() )
    lin_interp_var = #LinearInterpolation.extrapolate(
        LinearInterpolation.make_interpolator_linear_var(chunk)  #, Flat() )
    flux_out .= lin_interp_flux(grid./boost_factor)
    var_out .= lin_interp_var(grid./boost_factor)
    return flux_out
end

"""
   `interp_chunk_to_shifted_grid_linear( chunk_of_spectrum, wavelengths, boost_factor )
Return spectra interpolated onto a grid of points using linear interpolation.
# Arguments:
- chunk_of_spectrum
- wavelengths: AbstractRange or AbstractArray of locations where chunk is to be interpolated to
- boost_factor: divide wavelengths by boost_factor
Optional Arguments:
- Filter: Vector with pre-allocated workspace (if length>=1)
# Returns
- flux_out
"""
function interp_chunk_to_shifted_grid_linear( chunk::AC, grid::AR, boost_factor::Real ) where {  AC<:AbstractChunkOfSpectrum, T2<:Real, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    flux_out = Array{Float64,1}(undef,length(grid))
    var_out = Array{Float64,1}(undef,length(grid))
    interp_chunk_to_shifted_grid_linear!(flux_out, var_out, chunk, grid, boost_factor)
    return (flux=flux_out, var=var_out)
end

"""
   `interp_chunk_to_grid_linear!( flux_out, var_out, chunk_of_spectrum, wavelengths )`
Return spectra interpolated onto a grid of points using linear interpolation.
# Arguments:
- flux_out: (results stored into this array)
- var_out: (results stored into this array)
- chunk_of_spectrum
- wavelengths: AbstractRange or AbstractArray of locations where chunk is to be interpolated to
Optional Arguments:
- Filter: Vector with pre-allocated workspace (if length>=1)
# Returns
- flux_out
"""
function interp_chunk_to_grid_linear!( flux_out::AbstractArray{T1,1}, var_out::AbstractArray{T2,1}, chunk::AC, grid::AR ) where {
    T1<:Real, T2<:Real, AC<:AbstractChunkOfSpectrum, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    @assert size(flux_out) == size(var_out)
    @assert size(flux_out) == size(grid)
    lin_interp_flux = LinearInterpolation.make_interpolator_linear_flux(chunk)
    lin_interp_var = LinearInterpolation.make_interpolator_linear_var(chunk)
    flux_out .= lin_interp_flux(grid)
    var_out .= lin_interp_var(grid)
    return flux_out
end

"""
   `interp_chunk_to_grid_linear( chunk_of_spectrum, wavelengths )`
Return spectra interpolated onto a grid of points using linear interpolation.
# Arguments:
- chunk_of_spectrum
- wavelengths: AbstractRange or AbstractArray of locations where chunk is to be interpolated to
# Returns
- flux_out
"""
function interp_chunk_to_grid_linear( chunk::AC, grid::AR ) where {  AC<:AbstractChunkOfSpectrum, T2<:Real, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    flux_out = Array{Float64,1}(undef,length(grid))
    var_out = Array{Float64,1}(undef,length(grid))
    interp_chunk_to_grid_linear!(flux_out, var_out, chunk, grid)
    return (flux=flux_out, var=var_out)
end


text:
module mmrrigid
using FinEtools
using Compat.Test
import Arpack: eigs
import LinearAlgebra: norm
function test()

# println("""
#
# Example from Boundary element acoustics: Fundamentals and computer codes, TW Wu, page 123.
# Internal resonance problem. Reference frequencies: 90.7895, 181.579, 215.625, 233.959, 272.368, 281.895
# Quadrilateral mesh.
# """)

t0 = time()

rho=1.21*1e-9;# mass density
c =345.0*1000;# millimeters per second
bulk= c^2*rho;
Lx=1900.0;# length of the box, millimeters
Ly=800.0; # length of the box, millimeters
n=14;#
neigvs=18;
OmegaShift=10.0;

fens,fes = Q4block(Lx,Ly,n,n); # Mesh

geom = NodalField(fens.xyz)
P = NodalField(fill(zero(FFlt), size(fens.xyz,1),1))

numberdofs!(P)

femm = FEMMAcoust(IntegData(fes, GaussRule(2, 2)),
                     MatAcoustFluid(bulk,rho))

S = acousticstiffness(femm, geom, P);
C = acousticmass(femm, geom, P);

d,v,nev,nconv =eigs(C+OmegaShift*S, S; nev=neigvs, which=:SM)
d = d .- OmegaShift;
fs=real(sqrt.(complex(d)))/(2*pi)
# println("Eigenvalues: $fs [Hz]")
#
#
# println("Total time elapsed = ",time() - t0,"s")
# println("$fs[2]-90.7895")
@test abs(fs[2]-90.9801) < 1e-4

File =  "rigid_box.vtk"
scalarllist = Any[]
for n  = 1:15
  push!(scalarllist, ("Pressure_mode_$n", v[:,n]));
end
vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.Q4;
          scalars=scalarllist)
sleep(1.0)
try rm(File) catch end
# @async run(`"paraview.exe" $File`)

true
end
end
using .mmrrigid
mmrrigid.test()

module fahyL2example

using FinEtools
using Compat.Test
import Arpack: eigs
function test()
# println("""
# Example from Sound and Structural Vibration, Second Edition: Radiation, Transmission and Response [Paperback]
# Frank J. Fahy, Paolo Gardonio, page 483.
#
# 1D mesh.
# """)

t0  =  time()

rho = 1.21*1e-9;# mass density
c  = 343.0*1000;# millimeters per second
bulk =  c^2*rho;
L = 500.0;# length of the box, millimeters
A = 200.0; # cross-sectional area of the box
graphics =  true;# plot the solution as it is computed?
n = 40;#
neigvs = 8;
OmegaShift = 10.0;

fens,fes  =  L2block(L,n); # Mesh

geom = NodalField(fens.xyz)
P = NodalField(fill(zero(FFlt), size(fens.xyz,1),1))

numberdofs!(P)

femm = FEMMAcoust(IntegData(fes, GaussRule(1, 2)), MatAcoustFluid(bulk, rho))

S  =  acousticstiffness(femm, geom, P);
C  =  acousticmass(femm, geom, P);

d,v,nev,nconv  = eigs(C+OmegaShift*S, S; nev = neigvs, which = :SM)
d  =  d .- OmegaShift;
fs = real(sqrt.(complex(d)))/(2*pi)
# println("Eigenvalues: $fs [Hz]")

# println("Total time elapsed  =  ",time() - t0,"s")

# using Plots
# plotly()
# en = 2
# ix = sortperm(geom.values[:])
# plot(geom.values[:][ix], v[:,en][ix], color = "blue",
# title = "Fahy example, mode $en" , xlabel = "x", ylabel = "P")
# gui()

  @test (fs[2]-343.08)/343.08 < 0.001

  #println("Total time elapsed = ",time() - t0,"s")

  # using Winston
  # en=2
  # pl = FramedPlot(title="Fahy example, mode $en",xlabel="x",ylabel="P")
  # setattr(pl.frame, draw_grid=true)
  # ix=sortperm(geom.values[:])
  # add(pl, Curve(geom.values[:][ix],v[:,en][ix], color="blue"))

  # display(pl)

  true
end
end
using .fahyL2example
fahyL2example.test()

module mmfahyH8example
using FinEtools
using Compat.Test
import Arpack: eigs
function test()

# println("""
# Example from Sound and Structural Vibration, Second Edition: Radiation, Transmission and Response [Paperback]
# Frank J. Fahy, Paolo Gardonio, page 483.
#
# Hexahedral mesh.
# """)

t0 = time()

rho=1.21*1e-9;# mass density
c =343.0*1000;# millimeters per second
bulk= c^2*rho;
L=500.0;# length of the box, millimeters
A=200.0; # cross-sectional area of the box
n=40;#
neigvs=8;
OmegaShift=10.0;

fens,fes = H8block(L,sqrt(A),sqrt(A),n,1,1); # Mesh

geom = NodalField(fens.xyz)
P = NodalField(fill(zero(FFlt), size(fens.xyz,1),1))

numberdofs!(P)

femm = FEMMAcoust(IntegData(fes, GaussRule(3, 2)), MatAcoustFluid(bulk, rho))

S = acousticstiffness(femm, geom, P);
C = acousticmass(femm, geom, P);

d,v,nev,nconv = eigs(C+OmegaShift*S, S; nev=neigvs, which=:SM)
d = d .- OmegaShift;
fs = real(sqrt.(complex(d)))/(2*pi)
# println("Eigenvalues: $fs [Hz]")
#
#
# println("Total time elapsed = ",time() - t0,"s")

# File =  "fahy_H8.vtk"
# en = 5;
# vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H8;
#           scalars=[("Pressure_mode_$en", v[:,en])])
# println("Done")

  @test (fs[2]-343.08)/343.08 < 0.001

  true
end
end
using .mmfahyH8example
mmfahyH8example.test()

module mmfahyH27example
using FinEtools
using Compat.Test
import Arpack: eigs
function test()

# println("""
# Example from Sound and Structural Vibration, Second Edition: Radiation, Transmission and Response [Paperback]
# Frank J. Fahy, Paolo Gardonio, page 483.
#
# Hexahedral mesh.
# """)

t0 = time()

rho=1.21*1e-9;# mass density
c =343.0*1000;# millimeters per second
bulk= c^2*rho;
L=500.0;# length of the box, millimeters
A=200.0; # cross-sectional area of the box
n=40;#
neigvs=8;
OmegaShift=10.0;

fens,fes = H8block(L,sqrt(A),sqrt(A),n,1,1); # Mesh
fens,fes = H8toH27(fens,fes)

geom = NodalField(fens.xyz)
P = NodalField(fill(zero(FFlt), size(fens.xyz,1),1))

numberdofs!(P)

femm = FEMMAcoust(IntegData(fes, GaussRule(3, 3)), MatAcoustFluid(bulk, rho))

S = acousticstiffness(femm, geom, P);
C = acousticmass(femm, geom, P);

d,v,nev,nconv = eigs(C+OmegaShift*S, S; nev=neigvs, which=:SM)
d = d .- OmegaShift;
fs = real(sqrt.(complex(d)))/(2*pi)
# println("Eigenvalues: $fs [Hz]")
#
#
# println("Total time elapsed = ",time() - t0,"s")

# File =  "fahy_H8.vtk"
# en = 5;
# vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H8;
#           scalars=[("Pressure_mode_$en", v[:,en])])
# println("Done")

  @test (fs[2]-343.08)/343.08 < 0.001

  true
end
end
using .mmfahyH27example
mmfahyH27example.test()

module mstraight_duct_H8_1
using FinEtools
using Compat.Test
function test()
    t0  =  time()

    rho = 1.21*phun("kg/m^3");# mass density
    c  = 343.0*phun("m/s");# sound speed
    bulk =  c^2*rho;
    omega =  54.5901*phun("rev/s")
    vn0 =  -1.0*phun("m/s")
    Lx = 10.0*phun("m");# length of the box, millimeters
    Ly = 1.0*phun("m"); # length of the box, millimeters
    n = 20;#number of elements along the length

    # println("""
    #
    # Straight duct with anechoic termination.
    # Example from Boundary element acoustics: Fundamentals and computer codes, TW Wu, page 44.
    # Both real and imaginary components of the pressure should have amplitude of
    # rho*c = $(rho*c).
    #
    # Hexahedral mesh.
    # """)

    fens,fes  =  H8block(Lx,Ly,Ly,n,2,2); # Mesh
    bfes  =  meshboundary(fes)
    L0 = selectelem(fens,bfes, facing = true, direction = [-1.0 0.0 0.0])
    L10 = selectelem(fens,bfes, facing = true, direction = [+1.0 0.0 0.0])
    nLx = selectnode(fens,box = [0.0 Lx  0.0 0.0 0.0 0.0], inflate = Lx/1.0e5)

    geom  =  NodalField(fens.xyz)
    P = NodalField(fill(zero(Complex{Float64}), size(fens.xyz,1),1))

    numberdofs!(P)

    material = MatAcoustFluid(bulk,rho)
    femm  =  FEMMAcoust(IntegData(fes, GaussRule(3, 2)), material)

    S  =  acousticstiffness(femm, geom, P);
    C  =  acousticmass(femm, geom, P);

    E10femm  =  FEMMAcoustSurf(IntegData(subset(bfes,L10),GaussRule(2, 2)), material)
    D  =  acousticABC(E10femm, geom, P);

    E0femm  =  FEMMBase(IntegData(subset(bfes,L0), GaussRule(2,  2)))
    fi  =  ForceIntensity(-1.0im*omega*rho*vn0);
    F  =  distribloads(E0femm, geom, P, fi, 2);

    p = (-omega^2*S +omega*1.0im*D + C)\F
    scattersysvec!(P, p[:])

    # println("Pressure amplitude bounds: ")
    # println("  real $(minimum(real(P.values)))/$(maximum(real(P.values)))")
    # println("  imag $(minimum(imag(P.values)))/$(maximum(imag(P.values)))")
    #
    # println("Total time elapsed  =  ",time() - t0,"s")

    File  =   "straight_duct.vtk"
    scalars = real(P.values);
    vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.H8;
    scalars = [("Pressure", scalars)])
    try rm(File) catch end
    # @async run(`"paraview.exe" $File`)

    ref=rho*c
    @test abs(minimum(real(P.values))-(-ref))/ref < 0.02
    @test abs(minimum(imag(P.values))-(-ref))/ref < 0.02
    @test abs(maximum(real(P.values))-(+ref))/ref < 0.02
    @test abs(maximum(imag(P.values))-(+ref))/ref < 0.02

    # plotly()
    # ix = sortperm(geom.values[nLx,1])
    # plot(geom.values[nLx,1][ix], real(P.values)[nLx][ix], color = :blue, label = "real")
    # plot!(geom.values[nLx,1][ix], imag(P.values)[nLx][ix], color = :red, label  =  "imag")
    # plot!(title = "Straight duct with anechoic termination",
    # xlabel = "x", ylabel = "Pressure")
    # gui()

    true

end
end
using .mstraight_duct_H8_1
mstraight_duct_H8_1.test()

module mmsphere_dipole_1
using FinEtools
using Compat.Test
import LinearAlgebra: norm, lu, cross
function test()

# println("The interior sphere accelerates in the positive x-direction, generating
# positive pressure ahead of it, negative pressure behind
# ")
rho = 1.21*phun("kg/m^3");# mass density
c  = 343.0*phun("m/s");# sound speed
bulk =  c^2*rho;
a_amplitude=1.0*phun("mm/s^2");# amplitude of the  acceleration of the sphere
# omega = 2000*phun("rev/s");      # frequency of the incident wave
R = 5.0*phun("mm"); # radius of the interior sphere
Ro = 4*R # radius of the external sphere
P_amplitude = R*rho*a_amplitude; # pressure amplitude
nref=2;
nlayers=40;
tolerance = Ro/(nlayers)/100
frequency=2000; # Hz
omega=2*pi*frequency;
k = omega/c*[1.0 0.0 0.0];

# println("""
#         Spherical domain filled with acoustic medium, no scatterer.
# Hexahedral mesh.
#         """)

t0  =  time()

# Hexahedral mesh
fens,fes  =  H8sphere(R,nref);
fens1,fes1  =  mirrormesh(fens, fes, [-1.0, 0.0, 0.0], [0.0, 0.0, 0.0],
renumb =  r(c) = c[[1, 4, 3, 2, 5, 8, 7, 6]]);
fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
fes = cat(newfes1,fes2)

# Derive the finite element sets for the boundary
bfes  =  meshboundary(fes)
# Outer spherical boundary
function dout(xyz)
    return xyz/norm(xyz)
end

louter = selectelem(fens, bfes, facing = true, direction = dout)
outer_fes = subset(bfes, louter);

fens,fes  = H8extrudeQ4(fens, outer_fes,
nlayers, (xyz, layer)->(R+layer/nlayers*(Ro-R))*xyz/norm(xyz));
connected = findunconnnodes(fens, fes);
fens, new_numbering = compactnodes(fens, connected);
fess = renumberconn!(fes, new_numbering);

geom  =  NodalField(fens.xyz)
P = NodalField(fill(zero(FCplxFlt), size(fens.xyz,1),1))

bfes  =  meshboundary(fes)

# File  =   "Sphere.vtk"
# vtkexportmesh(File, bfes.conn, geom.values, FinEtools.MeshExportModule.Q4)
# @async run(`"paraview.exe" $File`)

# File  =   "Sphere.vtk"
# vtkexportmesh (File, fes.conn, fens.xyz, MeshExportModule.H8)
# @async run(`"C:/Program Files (x86)/ParaView 4.2.0/bin/paraview.exe" $File`)

linner = selectelem(fens, bfes, distance = R, from = [0.0 0.0 0.0],
  inflate = tolerance)
louter = selectelem(fens, bfes, facing = true, direction = dout)

# println("Pre-processing time elapsed  =  ",time() - t0,"s")

t1  =  time()

numberdofs!(P)

material = MatAcoustFluid(bulk,rho)
femm  =  FEMMAcoust(IntegData(fes, GaussRule(3, 2)), material)

S  =  acousticstiffness(femm, geom, P);
C  =  acousticmass(femm, geom, P);

abcfemm  =  FEMMAcoustSurf(IntegData(subset(bfes, louter), GaussRule(2, 2)), material)
D  =  acousticABC(abcfemm, geom, P);

# Inner sphere pressure loading
function dipole(dpdn, xyz, J, label)
    n = cross(J[:,1],J[:,2]);
    n = vec(n/norm(n));
    dpdn[1] = -rho*a_amplitude*n[1]
end

fi  =  ForceIntensity(FCplxFlt, 1, dipole);
#F  =  F - distribloads(pfemm, nothing, geom, P, fi, 2);
dipfemm  =  FEMMAcoustSurf(IntegData(subset(bfes, linner), GaussRule(2, 2)), material)
F  = distribloads(dipfemm, geom, P, fi, 2);

A = (1.0+0.0im)*(-omega^2*S + omega*1.0im*D + C)
K = lu(A) # We fake a complex matrix here
p = K\F  #

scattersysvec!(P, p[:])

# println(" Minimum/maximum pressure= $(minimum(real(p)))/$(maximum(real(p)))")
@test abs(-minimum(real(p))-maximum(real(p))) < 1.0e-9
@test abs(3.110989923540772e-6-maximum(real(p))) < 1.0e-9
# println("Computing time elapsed  =  ",time() - t1,"s")
# println("Total time elapsed  =  ",time() - t0,"s")

File  =   "sphere_dipole_1.vtk"
vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.H8;
 scalars = [( "realP", real(P.values))])
try rm(File) catch end
# @async run(`"paraview.exe" $File`)

true
end
end
using .mmsphere_dipole_1
mmsphere_dipole_1.test()

module mstraight_duct_T10_examplem
using FinEtools
using Compat.Test
function test()

t0  =  time()

rho = 1.21*phun("kg/m^3");# mass density
c  = 343.0*phun("m/s");# sound speed
bulk =  c^2*rho;
omega =  54.5901*phun("rev/s")
vn0 =  -1.0*phun("m/s")
Lx = 10.0*phun("m");# length of the box, millimeters
Ly = 1.0*phun("m"); # length of the box, millimeters
n = 20;#number of elements along the length

# println("""
#
# Straight duct with anechoic termination.
# Example from Boundary element acoustics: Fundamentals and computer codes, TW Wu, page 44.
# Both real and imaginary components of the pressure should have amplitude of
# rho*c = $(rho*c).
#
# Tetrahedral (quadratic) mesh.
# """)

fens,fes  =  T10block(Lx,Ly,Ly,n,2,2); # Mesh
bfes  =  meshboundary(fes)
L0 = selectelem(fens,bfes,facing = true, direction = [-1.0 0.0 0.0])
L10 = selectelem(fens,bfes,facing = true, direction = [+1.0 0.0 0.0])
nLx = selectnode(fens,box = [0.0 Lx  0.0 0.0 0.0 0.0], inflate = Lx/1.0e5)

geom  =  NodalField(fens.xyz)
P = NodalField(fill(zero(Complex{Float64}), size(fens.xyz,1),1))

numberdofs!(P)

material = MatAcoustFluid(bulk,rho)
femm  =  FEMMAcoust(IntegData(fes, TetRule(4)), material)

S  =  acousticstiffness(femm, geom, P);
C  =  acousticmass(femm, geom, P);

E10femm  =  FEMMAcoustSurf(IntegData(subset(bfes,L10), TriRule(3)), material)
D  =  acousticABC(E10femm, geom, P);

E0femm  =  FEMMBase(IntegData(subset(bfes,L0), TriRule(3)))
fi  =  ForceIntensity(-1.0im*omega*rho*vn0);
F  =  distribloads(E0femm, geom, P, fi, 2);

p = (-omega^2*S +omega*1.0im*D + C)\F
scattersysvec!(P, p[:])

# println("Pressure amplitude bounds: ")
# println("  real $(minimum(real(P.values)))/$(maximum(real(P.values)))")
# println("  imag $(minimum(imag(P.values)))/$(maximum(imag(P.values)))")
#
# println("Total time elapsed  =  ",time() - t0,"s")
@test abs(-414.0986190028914-minimum(imag(P.values))) < 1.e-6

File  =   "straight_duct.vtk"
scalars = real(P.values);
vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.T10;
scalars = [("Pressure", scalars)])
# @async run(`"paraview.exe" $File`)
try rm(File) catch end

# plotly()
# ix = sortperm(geom.values[nLx,1])
# plot(geom.values[nLx,1][ix], real(P.values)[nLx][ix], color = :blue, label = "real")
# plot!(geom.values[nLx,1][ix], imag(P.values)[nLx][ix], color = :red, label  =  "imag")
# plot!(title = "Straight duct with anechoic termination",
# xlabel = "x", ylabel = "Pressure")
# gui()

true
end
end
using .mstraight_duct_T10_examplem
mstraight_duct_T10_examplem.test()

module mmiintegrationmm
using FinEtools
using Compat.Test
import LinearAlgebra: norm
function test()
rho = 1000.0*phun("kg/m^3");# mass density of water
c  = 1.4491e+3*phun("m/s");# sound speed in water
bulk =  c^2*rho;
rhos = 2500.0*phun("kg/m^3");# mass density of the solid sphere
a_amplitude=1.0*phun("mm/s^2");# amplitude of the  acceleration of the sphere
R = 5.0*phun("mm"); # radius of the interior sphere
Ro = 3*R # radius of the external sphere
# Mesh parameters, geometrical tolerance
nperR = 6;
nlayers=10;
tolerance = Ro/(nlayers)/100
# Hexahedral mesh of the solid sphere
fens,fes  =  H8spheren(R, nperR);
r(c) = c[[1, 4, 3, 2, 5, 8, 7, 6]]
fens1,fes1  =  mirrormesh(fens, fes, [-1.0, 0.0, 0.0], [0.0, 0.0, 0.0],
          renumb =  r);
fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
fes = cat(newfes1,fes2)
fens1,fes1  =  mirrormesh(fens, fes, [0.0, -1.0, 0.0], [0.0, 0.0, 0.0],
          renumb =  r);
fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
fes = cat(newfes1,fes2)
fens1,fes1  =  mirrormesh(fens, fes, [0.0, 0.0, -1.0], [0.0, 0.0, 0.0],
          renumb =  r);
fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
fes = cat(newfes1,fes2)

# Find the inertial properties of the solid sphere
geom  =  NodalField(fens.xyz)
femm  =  FEMMBase(IntegData(fes, GaussRule(3, 2)))
V = integratefunction(femm, geom, (x) ->  1.0)
  # println("V=$(V/phun("mm^3"))")
Sx = integratefunction(femm, geom, (x) ->  x[1])
  # println("Sx=$(Sx/phun("mm^4"))")
Sy = integratefunction(femm, geom, (x) ->  x[2])
  # println("Sy=$(Sy/phun("mm^4"))")
Sz = integratefunction(femm, geom, (x) ->  x[3])
  # println("Sz=$(Sz/phun("mm^4"))")
CG = vec([Sx Sy Sz]/V)
  # println("CG=$(CG/phun("mm"))")
  I3 = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]
function Iinteg(x)
  (norm(x-CG)^2*I3-(x-CG)*(x-CG)')
end
Ixx = integratefunction(femm, geom, (x) ->  Iinteg(x)[1, 1])
Ixy = integratefunction(femm, geom, (x) ->  Iinteg(x)[1, 2])
Ixz = integratefunction(femm, geom, (x) ->  Iinteg(x)[1, 3])
Iyy = integratefunction(femm, geom, (x) ->  Iinteg(x)[2, 2])
Iyz = integratefunction(femm, geom, (x) ->  Iinteg(x)[2, 3])
Izz = integratefunction(femm, geom, (x) ->  Iinteg(x)[3, 3])
I = [Ixx Ixy Ixz; Ixy Iyy Iyz; Ixz Iyz Izz]
  # println("I=$(I/phun("mm^4"))")
@test abs(Ixx/phun("mm^4")-4.9809) < 0.001
mass =  V*rhos;
Inertia = I*rhos;
end
end
using .mmiintegrationmm
mmiintegrationmm.test()

module mmtransientsphere
using FinEtools
using Compat.Test
import LinearAlgebra: norm, cross
function test()
  # println("The interior sphere accelerates in the alternately in the positive
  # and negative x-direction, generating positive pressure ahead of it, negative
  # pressure behind. Time-dependent simulation.
  # ")
  rho = 1.21*phun("kg/m^3");# mass density
  c  = 343.0*phun("m/s");# sound speed
  bulk =  c^2*rho;
  a_amplitude=1.0*phun("mm/s^2");# amplitude of the  acceleration of the sphere
  # omega = 2000*phun("rev/s");      # frequency of the incident wave
  R = 50.0*phun("mm"); # radius of the interior sphere
  Ro = 8*R # radius of the external sphere
  P_amplitude = R*rho*a_amplitude; # pressure amplitude
  nref=2;
  nlayers=40;
  tolerance = Ro/(nlayers)/100
  frequency=1200.; # Hz
  omega=2*pi*frequency;
  dt=1.0/frequency/20;
  tfinal=90*dt;
  nsteps=round(tfinal/dt)+1;

  # Hexahedral mesh
  fens,fes  =  H8sphere(R,nref);
  fens1,fes1  =  mirrormesh(fens, fes, [-1.0, 0.0, 0.0], [0.0, 0.0, 0.0],
  renumb =  r(c) = c[[1, 4, 3, 2, 5, 8, 7, 6]]);
  fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
  fes = cat(newfes1,fes2)

  # Derive the finite element sets for the boundary
  bfes  =  meshboundary(fes)
  # Outer spherical boundary
  function dout(xyz)
      return xyz/norm(xyz)
  end

  louter = selectelem(fens, bfes, facing = true, direction = dout)
  outer_fes = subset(bfes, louter);

  fens,fes  = H8extrudeQ4(fens, outer_fes,
  nlayers, (xyz, layer)->(R+layer/nlayers*(Ro-R))*xyz/norm(xyz));
  connected = findunconnnodes(fens, fes);
  fens, new_numbering = compactnodes(fens, connected);
  fess = renumberconn!(fes, new_numbering);

  geom  =  NodalField(fens.xyz)
  P  =  NodalField(fill(zero(FCplxFlt), size(fens.xyz,1),1))

  bfes  =  meshboundary(fes)

  linner = selectelem(fens, bfes, distance = R, from = [0.0 0.0 0.0],
    inflate = tolerance)
  louter = selectelem(fens, bfes, facing = true, direction = dout)

  numberdofs!(P)

  material = MatAcoustFluid(bulk,rho)
  femm  =  FEMMAcoust(IntegData(fes, GaussRule(3, 2)), material)

  S  =  acousticstiffness(femm, geom, P);
  C  =  acousticmass(femm, geom, P);

  abcfemm  =  FEMMAcoustSurf(IntegData(subset(bfes, louter), GaussRule(2, 2)), material)
  D  =  acousticABC(abcfemm, geom, P);

  # Inner sphere pressure loading
  function dipole(dpdn, xyz, J, label, t)
      n = cross(J[:,1],J[:,2]);
      n = vec(n/norm(n));
      dpdn[1] = -rho*a_amplitude*sin(omega*t)*n[1]
  end

  dipfemm  =  FEMMAcoustSurf(IntegData(subset(bfes, linner), GaussRule(2, 2)), material)

  # Solve
  P0 = deepcopy(P)
  P0.values[:] .= 0.0; # initially all pressure is zero
  vP0 = gathersysvec(P0);
  vP1 = fill!(similar(vP0), zero(eltype(vP0)))
  vQ0 = fill!(similar(vP0), zero(eltype(vP0)))
  vQ1 = fill!(similar(vP0), zero(eltype(vP0)))
  t = 0.0;
  P1 = deepcopy(P0);

  fi  =  ForceIntensity(FCplxFlt, 1, (dpdn, xyz, J, label)->dipole(dpdn, xyz, J, label, t));
  La0 = distribloads(dipfemm, geom, P1, fi, 2);

  A = (2.0/dt)*S + D + (dt/2.)*C;

  step =0;
  while t <= tfinal
    step = step  + 1;
    # println("Time $t ($(step)/$(round(tfinal/dt)+1))")
    t = t+dt;
    fi  =  ForceIntensity(FCplxFlt, 1, (dpdn, xyz, J, label)->dipole(dpdn, xyz, J, label, t));
    La1 = distribloads(dipfemm, geom, P1, fi, 2);
    vQ1 = A\((2/dt)*(S*vQ0)-D*vQ0-C*(2*vP0+(dt/2)*vQ0)+La0+La1);
    vP1 = vP0 + (dt/2)*(vQ0+vQ1);
    vP0 = vP1;
    vQ0 = vQ1;
    P1 = scattersysvec!(P1, vec(vP1));
    P0 = P1;
    La0 = La1;
  end

  @test abs(vP1[1] - (3.0254474104972927e-6 + 0.0im)) < 1.e-10

  # File  =   "sphere_dipole_1.vtk"
  # vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H8;
  #  scalars = [( "realP", real(P1.values))])
  # @async run(`"paraview.exe" $File`)

  true

end
end
using .mmtransientsphere
mmtransientsphere.test()

module mmhhemispheremm
using FinEtools
using Compat.Test
using SparseArrays
import LinearAlgebra: norm, dot, lu, diff, cross
function test()

  # println("Rigid movable hemisphere in  water. Time-dependent simulation.
  # ")
  rho = 1000.0*phun("kg/m^3");# mass density of water
  c  = 1.4491e+3*phun("m/s");# sound speed in water
  bulk =  c^2*rho;
  rhos = 2500.0*phun("kg/m^3");# mass density of the solid sphere
  a_amplitude=1.0*phun("mm/s^2");# amplitude of the  acceleration of the sphere
  R = 5.0*phun("mm"); # radius of the interior sphere
  Ro = 3*R # radius of the external sphere
  P_amplitude = 1000*phun("Pa"); # pressure amplitude
  frequency=200.; # Hz
  omega=2*pi*frequency;
  wavenumber=omega/c;
  angl = 0.0
  front_normal=vec([cos(pi/180*angl),sin(pi/180*angl),0]);
  wavevector=wavenumber*front_normal;
  tshift = 0.0
  uincA=P_amplitude/rho/c/omega; # amplitude of the displacement
  # This is the analytical solution of Hickling and Wang
  uamplHW=P_amplitude/(rho*c^2)/(omega/c)*3/(2*rhos/rho+1);
  # Mesh parameters, geometrical tolerance
  nperR = 8;
  nlayers=nperR;
  tolerance = Ro/(nlayers)/100
  # Time step
  dt = 1.0/frequency/20;
  tfinal = 9/frequency;
  nsteps = Int(round(tfinal/dt)+1);

  t0  =  time()

  # Hexahedral mesh of the solid sphere
  fens,fes  =  H8spheren(R, nperR);
  r(c) = c[[1, 4, 3, 2, 5, 8, 7, 6]]
  fens1,fes1  =  mirrormesh(fens, fes, [-1.0, 0.0, 0.0], [0.0, 0.0, 0.0],
            renumb = r);
  fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
  fes = cat(newfes1,fes2)
  fens1,fes1  =  mirrormesh(fens, fes, [0.0, -1.0, 0.0], [0.0, 0.0, 0.0],
            renumb = r);
  fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
  fes = cat(newfes1,fes2)
  fens1,fes1  =  mirrormesh(fens, fes, [0.0, 0.0, -1.0], [0.0, 0.0, 0.0],
            renumb = r);
  fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
  fes = cat(newfes1,fes2)

  # Retain only half for the hemisphere
  Solidl = selectelem(fens, fes, box=[-Inf,Inf,-Inf,Inf,0,Inf], inflate=(Ro-R)/nlayers/100);
  # This is the half sphere which is filled with water
  fens1w, fes1w = deepcopy(fens), subset(fes, setdiff(collect(1:count(fes)), Solidl))

  # Find the inertial properties of the solid sphere
  geom  =  NodalField(fens.xyz)
  femm  =  FEMMBase(IntegData(subset(fes, Solidl), GaussRule(3, 2)))
  V = integratefunction(femm, geom, (x) ->  1.0)
    # println("V=$(V/phun("mm^3"))")
  Sx = integratefunction(femm, geom, (x) ->  x[1])
    # println("Sx=$(Sx/phun("mm^4"))")
  Sy = integratefunction(femm, geom, (x) ->  x[2])
    # println("Sy=$(Sy/phun("mm^4"))")
  Sz = integratefunction(femm, geom, (x) ->  x[3])
    # println("Sz=$(Sz/phun("mm^4"))")
  CG = vec([Sx Sy Sz]/V)
    # println("CG=$(CG/phun("mm"))")
    I3 = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]
  function Iinteg(x)
    (norm(x-CG)^2*I3-(x-CG)*(x-CG)')
  end
  Ixx = integratefunction(femm, geom, (x) ->  Iinteg(x)[1, 1])
  Ixy = integratefunction(femm, geom, (x) ->  Iinteg(x)[1, 2])
  Ixz = integratefunction(femm, geom, (x) ->  Iinteg(x)[1, 3])
  Iyy = integratefunction(femm, geom, (x) ->  Iinteg(x)[2, 2])
  Iyz = integratefunction(femm, geom, (x) ->  Iinteg(x)[2, 3])
  Izz = integratefunction(femm, geom, (x) ->  Iinteg(x)[3, 3])
  It = [Ixx Ixy Ixz; Ixy Iyy Iyz; Ixz Iyz Izz]
    # println("It=$(It/phun("mm^4"))")
  mass =  V*rhos;
  Inertia = It*rhos;

  # The boundary surface of the solid sphere will now be extruded to form  the
  # layers of fluid around it
  # Outer spherical boundary
  function dout(xyz)
      return xyz/norm(xyz)
  end
  bfes  =  meshboundary(fes)
  louter = selectelem(fens, bfes, facing = true, direction = dout)
  outer_fes = subset(bfes, louter);
  fens, fes  = H8extrudeQ4(fens, outer_fes,
            nlayers, (xyz, layer)->(R+layer/nlayers*(Ro-R))*xyz/norm(xyz));
  fens, newfes1, fes2 =  mergemeshes(fens, fes, fens1w, fes1w, tolerance)
  fes = cat(newfes1, fes2)
  connected = findunconnnodes(fens, fes);
  fens, new_numbering = compactnodes(fens, connected);
  fes = renumberconn!(fes, new_numbering);

  # File  =   "Sphere1.vtk"
  # vtkexportmesh(File, fes.conn, fens.xyz, FinEtools.MeshExportModule.H8)
  # @async run(`"paraview.exe" $File`)

  geom  =  NodalField(fens.xyz)
  P  =  NodalField(fill(zero(FFlt), size(fens.xyz,1),1))

  bfes  =  meshboundary(fes)

  linner = selectelem(fens, bfes, distance = R, from = [0.0 0.0 0.0],
    inflate = tolerance)
  louter = selectelem(fens, bfes, facing = true, direction = dout)

  # println("Pre-processing time elapsed  =  ",time() - t0,"s")

  t1  =  time()

  numberdofs!(P)

  material = MatAcoustFluid(bulk,rho)
  femm  =  FEMMAcoust(IntegData(fes, GaussRule(3, 2)), material)

  S  =  acousticstiffness(femm, geom, P);
  C  =  acousticmass(femm, geom, P);

  abcfemm  =  FEMMAcoustSurf(IntegData(subset(bfes, louter), GaussRule(2, 2)), material)
  D  =  acousticABC(abcfemm, geom, P);

  # ABC surface pressure loading
  function abcp(dpdn, xyz, J, label, t)
      n = cross(J[:,1],J[:,2]);
      n = vec(n/norm(n));
      arg = (-dot(vec(xyz),vec(wavevector))+omega*(t+tshift));
      dpdn[1] = P_amplitude*cos(arg)*(-dot(vec(n),vec(wavevector)));
  end

  targetfemm  =  FEMMAcoustSurf(IntegData(subset(bfes, linner), GaussRule(2, 2)), material)

  ForceF = GeneralField([zero(FFlt) for i=1:3, j=1:1])
  numberdofs!(ForceF)
  TorqueF = GeneralField([zero(FFlt) for i=1:3, j=1:1])
  numberdofs!(TorqueF)

  GF = pressure2resultantforce(targetfemm, geom, P, ForceF)
  GT = pressure2resultanttorque(targetfemm, geom, P, TorqueF, CG)

  H = transpose(GF)*(rho/mass)*GF + transpose(GT)*(sparse(rho*inv(Inertia)))*GT;

  Ctild = C + H;

  # Solve
  P0 = deepcopy(P)
  P0.values[:] .= 0.0; # initially all pressure is zero
  vP0 = gathersysvec(P0);
  vP1 = fill!(similar(vP0), zero(eltype(vP0)))
  vQ0 = fill!(similar(vP0), zero(eltype(vP0)))
  vQ1 = fill!(similar(vP0), zero(eltype(vP0)))
  t = 0.0;
  P1 = deepcopy(P0);
  tTa1 = [zero(FFlt) for i=1:3]
  tAa1 = [zero(FFlt) for i=1:3]
  tTa0 = [zero(FFlt) for i=1:3]
  tAa0 = [zero(FFlt) for i=1:3]
  tTv1 = [zero(FFlt) for i=1:3]
  tAv1 = [zero(FFlt) for i=1:3]
  tTv0 = [zero(FFlt) for i=1:3]
  tAv0 = [zero(FFlt) for i=1:3]
  tTu1 = [zero(FFlt) for i=1:3]
  tAu1 = [zero(FFlt) for i=1:3]
  tTu0 = [zero(FFlt) for i=1:3]
  tAu0 = [zero(FFlt) for i=1:3]
  tTa_store = fill(zero(FFlt), 3, nsteps+1)
  tAa_store = fill!(similar(tTa_store), zero(eltype(tTa_store)))
  tTv_store = fill!(similar(tTa_store), zero(eltype(tTa_store)))
  tAv_store = fill!(similar(tTa_store), zero(eltype(tTa_store)))
  tTu_store = fill!(similar(tTa_store), zero(eltype(tTa_store)))
  tAu_store = fill!(similar(tTa_store), zero(eltype(tTa_store)))
  t_store = fill(zero(FFlt), 1, nsteps+1)

  pinc = deepcopy(P)
  pincdd = deepcopy(P)
  pincv = fill(zero(FFlt), P.nfreedofs)
  pincddv = deepcopy(pincv)
  p1v = deepcopy(pincv)
  p1 = deepcopy(P)

  function recalculate_incident!(t, pinc, pincdd)
    for j = 1:count(fens)
      arg=(-dot(vec(fens.xyz[j,:]),vec(wavevector))+omega*(t+tshift));# wave along the wavevector
      pinc.values[j] = P_amplitude*sin(arg)
      pincdd.values[j] = -omega^2*P_amplitude*sin(arg)
    end
    return pinc, pincdd
  end

  pinc, pincdd = recalculate_incident!(t, pinc, pincdd)
  L0 = -S*gathersysvec!(pincdd, pincddv) - Ctild*gathersysvec!(pinc, pincv);
  fabcp(dpdn, xyz, J, label) = abcp(dpdn, xyz, J, label, t)
  fi  =  ForceIntensity(FFlt, 1, fabcp);
  La0 = distribloads(abcfemm, geom, P1, fi, 2);

  A = lu((2.0/dt)*S + D + (dt/2.)*Ctild);

  step = 1;
  while step <= nsteps
    # println("Time $t ($(step)/$(Int(round(tfinal/dt))+1))")
    # Store for output
    tTa_store[:, step] = tTa1
    tAa_store[:, step] = tAa1
    tTv_store[:, step] = tTv1
    tAv_store[:, step] = tAv1
    tTu_store[:, step] = tTu1
    tAu_store[:, step] = tAu1
    t_store[1, step] = t
    # New loads and update
    step = step  + 1;
    t = t+dt;
    pinc, pincdd = recalculate_incident!(t, pinc, pincdd)
    L1 = -S*gathersysvec!(pincdd, pincddv) - Ctild*gathersysvec!(pinc, pincv);
    fi  =  ForceIntensity(FFlt, 1, fabcp);
    La1 = distribloads(abcfemm, geom, P1, fi, 2);
    vQ1 = A\((2/dt)*(S*vQ0)-D*vQ0-Ctild*(2*vP0+(dt/2)*vQ0)+L0+L1+La0+La1);
    vP1 = vP0 + (dt/2)*(vQ0+vQ1);
    P1 = scattersysvec!(P1, vec(vP1));
    p1.values[:] = P1.values[:] + pinc.values[:]
    p1v = gathersysvec!(p1, p1v);
    Fresultant = GF*p1v;
    Tresultant = GT*p1v;
    # println("Fresultant = $Fresultant")
    tTa1 = Fresultant/mass; tAa1 = Inertia\Tresultant;
    # Update  the velocity and position of the rigid body
    tTv1 = tTv0 + dt/2*(tTa0+tTa1);
    tTu1 = tTu0 + dt/2*(tTv0+tTv1);
    tAv1 = tAv1 + dt/2*(tAa0+tAa1);
    tAu1 = tAu0 + dt/2*(tAv0+tAv1);
    # show(tTa1)
    # Reset for the next time step
    copyto!(vP0, vP1);
    copyto!(vQ0, vQ1);
    P0 = deepcopy(P1);
    copyto!(L0, L1);
    copyto!(La0, La1);
    copyto!(tTa0, tTa1); copyto!(tAa0, tAa1)
    copyto!(tTv0, tTv1); copyto!(tAv0, tAv1)
    copyto!(tTu0, tTu1); copyto!(tAu0, tAu1)
  end
  tTa_store[:, step] = tTa1
  tAa_store[:, step] = tAa1
  tTv_store[:, step] = tTv1
  tAv_store[:, step] = tAv1
  tTu_store[:, step] = tTu1
  tAu_store[:, step] = tAu1
  t_store[1, step] = t

  function remove_bias(it, iy)
    local t = deepcopy(reshape(it,length(it),1))
    local y = deepcopy(reshape(iy,length(iy),1))
    t = t/mean(t);
    t = reshape(t,length(t),1);
    y = reshape(y,length(y),1);
    local B = [t [1.0 for m in t]];
    local p = (B'*B)\(B'*y);
    return y.-p[1]*t.-p[2];
  end

  ux = remove_bias(t_store, tTu_store[1,:])

  function find_peaks(t,y)
      p=Float64[];
      dy=diff(y, dims=1)
      for j3=1:length(dy)-1
          if (dy[j3]*dy[j3+1]<0)
              push!(p, y[j3+1]);
          end
      end
      return p
  end
  function amplitude(t,y)
      p= find_peaks(t,y) ;
      A=mean(abs.(p)) ;
      return A
  end
  # println("Normalized amplitude  = $(amplitude(t_store,ux)/uamplHW)")
@test abs(amplitude(t_store,ux)/uamplHW-0.9361240886218495) < 1.0e-5
  # using Plots
  # plotly()
  # plot(transpose(t_store), transpose(tTu_store)/phun("mm"),
  #  xlabel = "t", ylabel = "Displacement", label = "Displacements")
  # plot(transpose(t_store), transpose(tAu_store),
  #  xlabel = "t", ylabel = "Displacement", label = "Rotations")
  # # plot(transpose(t_store), (ux)/uamplHW,
  # #  xlabel = "t", ylabel = "Normalized displacement", label = "Rectified")
  # gui()
  #
  # #
  # File  =   "hemisphere_final_td1.vtk"
  # vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H8;
  #  scalars = [( "realP", real(P1.values))])
  # @async run(`"paraview.exe" $File`)

end
end
using .mmhhemispheremm
mmhhemispheremm.test()

module mmbbaffledmm
using FinEtools
using Compat.Test
import LinearAlgebra: norm
function test()
rho = 1.21*phun("kg/m^3");# mass density
c  = 343.0*phun("m/s");# sound speed
bulk =  c^2*rho;
omega =  7500*phun("rev/s");      # frequency of the piston
a_piston =  -1.0*phun("mm/s")     # amplitude of the piston acceleration
R = 50.0*phun("mm");# radius of the piston
Ro = 150.0*phun("mm"); # radius of the enclosure
nref = 3;#number of refinements of the sphere around the piston
nlayers = 25;                     # number of layers of elements surrounding the piston
tolerance = R/(2^nref)/100

# println("""
#
# Baffled piston in a half-sphere domain with ABC.
#
# Hexahedral mesh. Algorithm version.
# """)

t0  =  time()

# Hexahedral mesh
fens,fes  =  H8sphere(R,nref);
bfes  =  meshboundary(fes)
# File  =   "baffledabc_boundary.vtk"
# vtkexportmesh(File, bfes.conn, fens.xyz, FinEtools.MeshExportModule.Q4)
#  @async run(`"paraview.exe" $File`)

l = selectelem(fens,bfes,facing = true,direction = [1.0 1.0  1.0], dotmin= 0.001)
ex(xyz, layer) = (R+layer/nlayers*(Ro-R))*xyz/norm(xyz)
fens1,fes1  =  H8extrudeQ4(fens, subset(bfes,l), nlayers, ex);
fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
fes = cat(newfes1,fes2)

# Piston surface mesh
bfes  =  meshboundary(fes)
l1 = selectelem(fens, bfes, facing = true, direction = [-1.0 0.0 0.0])
l2 = selectelem(fens, bfes, distance = R, from = [0.0 0.0 0.0], inflate = tolerance)
piston_fes = subset(bfes,intersect(l1,l2));

# Outer spherical boundary
louter = selectelem(fens, bfes, facing = true, direction = [1.0 1.0  1.0], dotmin= 0.001)
outer_fes = subset(bfes,louter);

# println("Pre-processing time elapsed  =  ",time() - t0,"s")

t1  =  time()

material = MatAcoustFluid(bulk, rho)
# Region of the fluid
region1 =  FDataDict("femm"=>FEMMAcoust(IntegData(fes, GaussRule(3, 2)), material))

# Surface for the ABC
abc1  =  FDataDict("femm"=>FEMMAcoustSurf(IntegData(outer_fes, GaussRule(2, 2)),
          material))

# Surface of the piston
flux1  =  FDataDict("femm"=>FEMMAcoustSurf(IntegData(piston_fes, GaussRule(2, 2)),
          material),  "normal_flux"=> -rho*a_piston+0.0im);

# Make model data
modeldata =  FDataDict("fens"=>  fens,
                 "omega"=>omega,
                 "regions"=>[region1],
                 "flux_bcs"=>[flux1], "ABCs"=>[abc1])

# Call the solver
modeldata = FinEtools.AlgoAcoustModule.steadystate(modeldata)

# println("Computing time elapsed  =  ",time() - t1,"s")
# println("Total time elapsed  =  ",time() - t0,"s")

geom = modeldata["geom"]
P = modeldata["P"]

@test abs(P.values[1]-(4.355914465396856e-6 - 1.2061599990585114e-6im)) < 1.e-10
# File  =   "baffledabc.vtk"
# vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H8;
# scalars = [("absP", abs.(P.values))])
# @async run(`"paraview.exe" $File`)

# using Winston
# pl  =  FramedPlot(title = "Matrix",xlabel = "x",ylabel = "Re P, Im P")
# setattr(pl.frame, draw_grid = true)
# add(pl, Curve([1:length(C[:])],vec(C[:]), color = "blue"))

# # pl = plot(geom.values[nLx,1][ix],scalars[nLx][ix])
# # xlabel("x")
# # ylabel("Pressure")
# display(pl)

true
end
end
using .mmbbaffledmm
mmbbaffledmm.test()

module mmtransientmm1mm
using FinEtools
using Compat.Test
import LinearAlgebra: norm
function test()
    rho = 1.21*phun("kg/m^3");# mass density
    c  = 343.0*phun("m/s");# sound speed
    bulk =  c^2*rho;
    frequency=2000; # Hz
    omega=2*pi*frequency;
    P_piston = 1.0e-3*phun("MPa"); # amplitude of the piston pressure
    R = 50.0*phun("mm");# radius of the piston
    Ro = 400.0*phun("mm"); # radius of the enclosure
    nref = 2;#number of refinements of the sphere around the piston
    nlayers = 48;                     # number of layers of elements surrounding the piston
    tolerance = R/(2^nref)/10
    dt=1.0/frequency/20;
    tfinal=50*dt;
    nsteps=round(tfinal/dt)+1;

    # println("""
    #
    # Baffled piston in a half-sphere domain with hard boundary condition.
    #
    # Hexahedral mesh. Algorithm version.
    # """)

    t0  =  time()

    # Hexahedral mesh
    fens,fes  =  H8sphere(R,nref);
    bfes  =  meshboundary(fes)
    # File  =   "baffledabc_boundary.vtk"
    # vtkexportmesh(File, bfes.conn, fens.xyz, FinEtools.MeshExportModule.Q4)
    #  @async run(`"paraview.exe" $File`)

    l = selectelem(fens,bfes,facing = true,direction = [1.0 1.0  1.0], dotmin= 0.001)
    ex(xyz, layer) = (R+layer/nlayers*(Ro-R))*xyz/norm(xyz)
    fens1,fes1  =  H8extrudeQ4(fens, subset(bfes,l), nlayers, ex);
    fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
    fes = cat(newfes1,fes2)

    # Piston surface mesh
    bfes  =  meshboundary(fes)
    l1 = selectelem(fens, bfes, facing = true, direction = [-1.0 0.0 0.0])
    l2 = selectelem(fens, bfes, distance = R, from = [0.0 0.0 0.0], inflate = tolerance)
    piston_fes = subset(bfes,intersect(l1,l2));
    # File  =   "piston_fes.vtk"
    # vtkexportmesh(File, piston_fes.conn, fens.xyz, FinEtools.MeshExportModule.Q4)
    #  @async run(`"paraview.exe" $File`)

    # Outer spherical boundary
    louter = selectelem(fens, bfes, facing = true, direction = [1.0 1.0  1.0], dotmin= 0.001)
    outer_fes = subset(bfes,louter);
    # File  =   "outer_fes.vtk"
    # vtkexportmesh(File, outer_fes.conn, fens.xyz, FinEtools.MeshExportModule.Q4)
    #  @async run(`"paraview.exe" $File`)

    # println("Pre-processing time elapsed  =  ",time() - t0,"s")

    t1  =  time()

    material = MatAcoustFluid(bulk, rho)
    # Region of the fluid
    femm = FEMMAcoust(IntegData(fes, GaussRule(3, 2)), material)

    geom = NodalField(fens.xyz);
    P = NodalField(fill(zero(FFlt), nnodes(geom),1));
    piston_fenids = connectednodes(piston_fes);
    setebc!(P, piston_fenids, true, 1, 0.0);
    applyebc!(P);
    numberdofs!(P);

    S  =  acousticstiffness(femm, geom, P);
    C  =  acousticmass(femm, geom, P);
    D= (2.0/dt)*S +(dt/2.0)*C;

    P0 = deepcopy(P)
    Pdd0 = deepcopy(P)
    P1 = deepcopy(P)
    Pdd1 = deepcopy(P)
    TMPF = deepcopy(P)
    P = nothing # we don't need this field anymore
    vP0 = fill(zero(FFlt), P0.nfreedofs)
    vP1 = deepcopy(vP0)
    vPd0 = deepcopy(vP0)
    vPd1 = deepcopy(vP0)

    t = 0.0
    P0.fixed_values[piston_fenids,1] .= P_piston*sin(omega*t)
    Pdd0.fixed_values[piston_fenids,1] .= P_piston*(-omega^2)*sin(omega*t)
    vP0 = gathersysvec!(P0, vP0)

    nh = selectnode(fens, nearestto = [R+Ro/2, 0.0, 0.0] )
    Pnh = [P1.values[nh, 1][1]]
    ts = [t]
    step =0;
    while t <=tfinal
        step = step  +1;
        t=t+dt;
        P1.fixed_values[piston_fenids,1] .= P_piston*sin(omega*t)
        Pdd1.fixed_values[piston_fenids,1] .= P_piston*(-omega^2)*sin(omega*t)
        TMPF.fixed_values = P0.fixed_values + P1.fixed_values
        F = nzebcloadsacousticmass(femm, geom, TMPF);
        TMPF.fixed_values = Pdd0.fixed_values + Pdd1.fixed_values
        F = F + nzebcloadsacousticstiffness(femm, geom, TMPF);
        # println("$(norm(F))")
        vPd1 = D\((2/dt)*(S*vPd0) - C*(2*vP0+(dt/2)*vPd0) + F);
        vP1 = vP0 + (dt/2)*(vPd0+vPd1);
        scattersysvec!(P1, vP1); # store current pressure
        # Swap variables for the next step
        copyto!(vP0, vP1)
        copyto!(vPd0, vPd1)
        copyto!(P0, P1)
        copyto!(Pdd0, Pdd1)
        # Graphics output
        # File  =   "baffled_piston-$(step).vtk"
        # vtkexportmesh(File, fes.conn, fens.xyz, FinEtools.MeshExportModule.H8;
        #     scalars = [(  "P", P1.values)])
        #  @async run(`"paraview.exe" $File`)
        push!(Pnh, P1.values[nh, 1][1])
        push!(ts, t)
        # println("step = $( step )")
    end

    # using Plots
    # plotly()
    # plot(vec(ts), vec(Pnh))
    # gui()
    # println("$(Pnh[end])")
    # println("$(length(Pnh))")
    # println("$(Pnh[end])")
@test abs(Pnh[end]-221.11319820621947) < 1.0e-3
end
end
using .mmtransientmm1mm
mmtransientmm1mm.test()

module mmAnnularmm
using FinEtools
using FinEtools.AlgoAcoustModule
using Compat.Test
function test()

    # println("""
    # Annular region, pressure BC + rigid wall.
    # This version uses the FinEtools algorithm module.
    # Version: 08/21/2017
    # """)

    t0 = time()

    rho = 1001*phun("kg/m^3");# mass density
    c  = 1500.0*phun("m/s");# sound speed
    bulk =  c^2*rho;
    omega =  2000*phun("rev/s");      # frequency of the piston
    rin =  1.0*phun("m");#internal radius

    rex =  2.0*phun("m"); #external radius
    nr = 20; nc = 120;
    Angle = 2*pi;
    thickness =  1.0*phun("m/s");
    tolerance = min(rin/nr,  rin/nc/2/pi)/10000;

    fens, fes = Q4annulus(rin, rex, nr, nc, Angle)
    fens, fes = mergenodes(fens,  fes,  tolerance);
    edge_fes = meshboundary(fes);

    # The pressure boundary condition
    l1 = selectelem(fens, edge_fes, box=[-1.1*rex -0.9*rex -0.5*rex 0.5*rex]);
    ebc1 = FDataDict("node_list"=>connectednodes(subset(edge_fes, l1)),
     "pressure"=>x -> cos(2*pi*x[2]/rin)+1im*sin(2*pi*x[2]/rin)) # entering the domain

    material = MatAcoustFluid(bulk, rho)
    femm = FEMMAcoust(IntegData(fes,  GaussRule(2, 2)),  material)
    region1 = FDataDict("femm"=>femm)

    # Make model data
    modeldata = FDataDict("fens"=>fens,
    "omega"=>omega,
    "regions"=>[region1], "essential_bcs"=>[ebc1]);

    # Call the solver
    modeldata = FinEtools.AlgoAcoustModule.steadystate(modeldata)
    geom=modeldata["geom"]
    P=modeldata["P"]
    # println("Minimum/maximum pressure, real= $(minimum(real(P.values)))/$(maximum(real(P.values))))")
    # println("Minimum/maximum pressure, imag= $(minimum(imag(P.values)))/$(maximum(imag(P.values))))")
    @test abs(minimum(real(P.values)) - -26.02026924534437) < 1.0e-5
    @test abs(maximum(real(P.values)) - 32.15639028990587) < 1.0e-5
    @test abs(minimum(imag(P.values)) - -3.083023162317023) < 1.0e-5
    @test abs(maximum(imag(P.values)) - 3.083023162317045) < 1.0e-5
    # println("Total time elapsed = ",time() - t0,"s")

    # # Postprocessing
    # File = "acou_annulusmod.vtk"
    # vtkexportmesh(File, fes.conn, geom.values,
    # FinEtools.MeshExportModule.Q4; scalars=[("Pre", real(P.values)), ("Pim", imag(P.values))])
    # @async run(`"paraview.exe" $File`)

end
end
using .mmAnnularmm
mmAnnularmm.test()


text:
#---------------------------------------------------------
# # Ideal mean age
#---------------------------------------------------------

#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`ideal_mean_age.ipynb`](@__NBVIEWER_ROOT_URL__/examples/generated/ideal_mean_age.ipynb)

#---------------------------------------------------------
# ## The model
#---------------------------------------------------------

# We will simulate the ideal mean age of water.
# That is, the average amount of time since a water parcel had last contact with the surface.

# ### Tracer equation
#
# The ideal mean age is transported with water, is equal to $0$ at the surface, and increases by one second every second everywhere.
# In other words, the 3D field of the age, $a$, is governed by the tracer equation
#
# $$\frac{\partial a}{\partial t} + \nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \cdot \nabla \right] a = 1,$$
#
# where $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \cdot \nabla \right]$ is a differential operator that represents the transport by the ocean circulation.
# ($\boldsymbol{u}$ is the 3D vector field for the advection and $\mathbf{K}$ is the diffusivity matrix.)
# In the equation above, we also assume that there is the boundary condition that $a=0$ at the surface.

#---------------------------------
# ### Discretized tracer equation
#---------------------------------

# In AIBECS, the linear differential operator defined by $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \cdot \nabla \right]$ is approximated by a constant matrix $\mathbf{T}$ when discretizing the continuous 3D ocean onto the model grid.
# This matrix can be small (e.g., for models with a few boxes), or large, like for the OCIM (more on the OCIM later).
# Similarly, the continuous 3D field of the age, $a$, is discretized into a column-vector, $\boldsymbol{a}$.
# (We represent scalars in italic, vectors in bold italic, and matrices in upstraight bold.)

# In the discrete case, we replace the boundary condition (that $a = 0$ at the surface) by imposing $\boldsymbol{a} = 0$ in the surface layer of the model grid.
# In practice, this is done by restoring $\boldsymbol{a}$ to $0$ with a very short timescale.
# The tracer equation thgus takes the form of
#
# $$\frac{\partial\boldsymbol{a}}{\partial t} = -\mathbf{T} \, \boldsymbol{a} + 1 - \boldsymbol{\Lambda} \, \boldsymbol{a},$$
#
# where $\boldsymbol{\Lambda}$ is a diagonal matrix with entries equal to $1 / \tau$ in the surface layer of the model grid and 0 otherwise.
# The timescale $\tau$ is chosen to be very small, ensuring that $\boldsymbol{a}$ is very close to $0$ at the surface.
# The first term represents the transport by the ocean circulation, the second term the source of 1 second per second everywhere, and the last term the fast relaxation.

#---------------------------------
# ### Steady-state
#---------------------------------

# The steady-state is the equilibrium that would be reached if we wait long enough for $a$ to not change anymore.
# Mathematically, the steady-state is also the state for which
#
# $$\frac{\partial a}{\partial t} = 0.$$

# Computationally, in the discrete case, this means that we just need to solve
#
# $$0 = -\mathbf{T} \, \boldsymbol{a} + 1 - \boldsymbol{\Lambda} \, \boldsymbol{a}$$
#
# to find $\boldsymbol{a}$.
# More specifically, we need to solve
#
# $$(\mathbf{T} + \boldsymbol{\Lambda}) \, \boldsymbol{a} = 1.$$

# Now that we have the equations laid down, let us chose the circulation transport matrix, $\mathbf{T}$.

#---------------------------------------------------------
# ## Using AIBECS
#---------------------------------------------------------

#md # !!! note
#md #     If this is the first time you are trying AIBECS, make sure you go through the prerequisites!
#nb # > **Note**
#nb # > If this is the first time you are trying AIBECS, make sure you go through the prerequisites!

# AIBECS can interpret tracer equations as long as you arrange them under the generic form:
#
# $$\frac{\partial \boldsymbol{x}}{\partial t} = \boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p}),$$
#
# where $\boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p})$ is the rate of change of the state and $\boldsymbol{p}$ is the vector of model parameters.
# We only track the age here, so that the entire state of the system is determined by the age itself.
# In other words, here, $\boldsymbol{x} = \boldsymbol{a}$.

# We will use AIBECS to find the steady-state of the system.
# For AIBECS, this translates into finding the solution of $\boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p}) = 0$.

# We start by telling Julia that we want to use the AIBECS package via

using AIBECS

#md # !!! note
#md #     If it's the first time you are running this line, the package will need precompiling.
#md #     This may take a minute or two.
#md #     (Just be patient... Or read on while you wait!)
#
#nb # > **Note**
#nb # > If it's the first time you are running this line, the package will need precompiling.
#nb # > This may take a minute or two.
#nb # > (Just be patient... Or read on while you wait!)
#
#md # !!! note
#md #     You should see a `Warning` for the `Flatten` package — just disregard it...
#md #     If you get an error though, please send me a copy of the output/error message, and I will try to troubleshoot it.
#nb # > **Note**
#nb # > You should see a `Warning` for the `Flatten` package — just disregard it...
#nb # > If you get an error though, please send me a copy of the output/error message, and I will try to troubleshoot it.

#---------------------------------
# ### The circulation
#---------------------------------

# We will use the circulation output from the Ocean Circulation Inverse Model (OCIM1).
# Basically, the OCIM provides researchers and oceanographers with a big sparse matrix that represents the global ocean circulation (advection and diffusion), which allows them to efficiently simulate the transport of passive tracers, like the age.
# (For more details, see Tim DeVries's [website](https://tdevries.eri.ucsb.edu/models-and-data-products/) and references therein.)
# With AIBECS, the OCIM0.1 and OCIM1 circulations can be loaded really easily, by simply typing

wet3D, grd, T_OCIM = AIBECS.OCIM1.load()
typeof(T_OCIM), size(T_OCIM)

#
#md # !!! note
#md #     Julia may ask you to download the OCIM matrix for you, in which case you should say yes (i.e., type `y`).
#md #     Once downloaded, AIBECS will remember where it downloaded the file and it will only load it from your laptop.
#nb # > **Note**
#nb # > Julia may ask you to download the OCIM matrix for you, in which case you should say yes (i.e., type `y`).
#nb # > Once downloaded, AIBECS will remember where it downloaded the file and it will only load it from your laptop.
#
# Additionally to downloading the OCIM file, the `load()` command loads 3 variables in the Julia workspace:
# - `wet3D` — a 3D array of the model grid, filled with `1`'s at "wet" grid boxes and `0`'s and "land" grid boxes.
# - `grd` — a `OceanGrid` object containing information about the 3D grid of the OCIM circulation, like the latitude, longitude, and depth of each grid boxes.
# - `T_OCIM` — the transport matrix representing advection and diffusion.
#
# The second line in command above tells you the type and the size of `T_OCIM`.
# It is a sparse matrix (CSC just means that it is stored in Compressed Sparse Column format) and is quite big!
#md # !!! note
#md #     A sparse matrix is just a matrix with very few non-zero entries.
#md #     Computationally, sparse matrices are stored differently than full matrices to save memory (no need to save all those zeros), and are much faster to use too!
#nb # > **Note**
#nb # > A sparse matrix is just a matrix with very few non-zero entries.
#nb # > Computationally, sparse matrices are stored differently than full matrices to save memory (no need to save all those zeros), and are much faster to use too!

# Anyway, this looks good, so let's move on with setting up the model!
#
# We have already loaded the transport matrix, `T_OCIM`, for the ocean circulation, but we must tell AIBECS that it applies to the age.
# To do that, we define a function of the parameters (although there are no parameters involved in this case, this is just the way AIBECS works for the moment).

T_age(p) = T_OCIM

# (Functions in Julia can be created in one line, just as above.)
# That's it for the circulation.
# Now, let's define the local sources and sinks.

#---------------------------------
# ### The local sources and sinks
#---------------------------------

# We will denote the age, $\boldsymbol{a}$, by the variable `age` in Julia.
# (It's good practice to use explicit names!)
# We need to translate the local sources and sinks in our discretized state function $\boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p})$ into Julia code.

# #### The source
#--------------------------

# Remember the age increases by $1$ second every second and everywhere.
# So its source function is equal to, well, `1`! (seconds per seconds means it is unitless).
# Let's create the local source function:

source_age(age, p) = 1

# #### The sink
#--------------------------

# Recall that the age must also be $0$ at the surface.
# And that we implement this boundary condition by restoring the age very fast to $0$ in the surface layer.
# This will act as the sink for the age.
# But first, we need to figure out where "the surface layer" is.
# To do that, AIBECS can generate a number of useful constants for you.
# (You can see the list of functions by typing `varinfo(AIBECS)` at the REPL.)
# Here we will use the vector of grid box depths, `z`, which AIBECS can generate for us via

z = vector_of_depths(wet3D, grd)

# So what is the top layer?
# Let's investigate what's the minimum depth:

minimum(z)

# The surface layer in the OCIM grid has its center at about $18\,$m depth.
# We can create a mask of the surface layer via `z .< 20`.
# (This will return a vector of `0`s and `1`s, depending on whether the depth, `z`, is less than `20`.)
#
#md # !!! note
#md #     In Julia (like in MATLAB), placing a dot, `.`, in front of operators is a convenient way to do element-wise operations.
#nb # > **Note**
#nb # > In Julia (like in MATLAB), placing a dot, `.`, in front of operators is a convenient way to do element-wise operations.
#
# Then, we implement the local sink by restoring the age to `0` with a timescale `τ`, via

function sink_age(age, p)
    τ = p.τ
    return age .* (z .< 20u"m") / τ
end

#md # !!! note
#md #     Julia allows you to use unicode for your functions and variables, like for `τ`.
#nb # > **Note**
#nb # > Julia allows you to use unicode for your functions and variables, like for `τ`.
#
# Here, we have defined a Julia function using the `function` keyword because the sink is a bit more complicated, so that we needed two lines to define it.
# The first line unpacks the model parameters, which is just the restoring timescale, `τ`, in this case.
# We will chose the value for `τ` later.

# #### Net sources and sinks
#--------------------------

# The sources minus the sinks are simply defined by

sms_age(age, p) = source_age(age, p) .- sink_age(age, p)

# #### Model parameters
#--------------------------

# We must define the parameters... And AIBECS comes with an API for that!

t = empty_parameter_table()    # initialize table of parameters
add_parameter!(t, :τ, 1u"s")   # add the parameter we want (τ = 1s)
initialize_Parameters_type(t, "IdealAgeParameters")  # Generate the parameter type
t

# Note, in particular, that we gave our parameter `τ` a unit.
# Yes, Julia comes with some nice functionality to deal with units directly!
# The lines above created a table that contains all the info for generating the parameters vector, $\boldsymbol{p}$.
# To generate the parameters in AIBECS we do:

p₀ = IdealAgeParameters()

# where we have used the constructor `IdealAgeParameters`, whose name we defined in the previous cell.
# Here we did not really need to create `p₀` as a parameters vector, since it has only one element, `τ`, in it.
# However, we are here to learn, and this structure and functionality comes in very handy when one deals with many parameters.
# (And as you can imagine, having all the parameters in a nice table ready for being used in a publication comes quite handy!)

# #### State function and Jacobian
#--------------------------

# Similarly to `p₀`, let's create a state `x₀` to start with.
# The vector `x₀` will be our initial guess for the state.
# Let's assume that the age is `1` (seconds) everywhere (as an initial guess):

nb = number_of_wet_boxes(wet3D)  # number of wet boxes
x₀ = ones(nb)

# The first line above defines the number of wet grid boxes, `nb`.
# Here, this is also the length of the state vector `x`, because there is only one tracer, `age`.
# In the second line, the `ones` function creates a vector of `1`s of the size you give it (the number of wet grid boxes, `nb`, here, which we defined earlier).

# Finally, the last step for the set up is to define $\boldsymbol{F}$.
# Using AIBECS, this is done via

T_matrices = (T_age,)           # bundles all the transport matrices in a tuple
sources_minus_sinks = (sms_age,) # bundles all the source-sink functions in a tuple
F, ∇ₓF = state_function_and_Jacobian(T_matrices, sources_minus_sinks, nb) # generates the state function (and its Jacobian!)
F(x₀,p₀)

# That's it!
# We have just created a model of the mean age.
#
#
# Lines 2 and 3 are just telling AIBECS
# - what transport matrices it should use for the transport of these tracers, and
# - and what local sources and sinks should be appplied to these tracers
#
#md # !!! note
#md #     The `(x,)` syntax returns a tuple of one element — the comma is necessary because without it, `(x)` would be just like `x` with brackets around it.
#md #     This interface of AIBECS was developed for case with multiple tracers in mind, and might look a bit odd for a single tracer.
#md #     But in the future, this might be cleaned up to be easier to work with single tracers.
#nb # > **Note**
#nb # > The `(x,)` syntax returns a tuple of one element — the comma is necessary because without it, `(x)` would be just like `x` with brackets around it.
#nb # > This interface of AIBECS was developed for case with multiple tracers in mind, and might look a bit odd for a single tracer.
#nb # > But in the future, this might be cleaned up to be easier to work with single tracers.

# The fourth line creates two functions:
# - `F` — the numerical version of the **state function**, $\boldsymbol{F}$, of our model of the mean age, and
# - `∇ₓF` — the **Jacobian matrix** of the state function, i.e., $\nabla_{\boldsymbol{x}}\boldsymbol{F}$.
# Yes, AIBECS just automatically created an exact derivative of your input, using autodifferentiation via dual numbers.
# (I'd be very excited to detail how this is implemented here, but it is an entirely different discussion.)

# The last line just checks that our generated `F` works with our initial guess `x₀` and parameter vector `p₀`.

# #### Solving for the steady-state
#--------------------------

# The Jacobian, `∇ₓF` is essential to solving the steady-state equation $\boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p}) = 0$ fast.
# Specifically, solving $\boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p}) = 0$ is done via Newton's method.
# By starting from an initial guess, that you will have to provide, it will iterate over this recursion relation
#
# $$\boldsymbol{x}_{k+1} = \boldsymbol{x}_{k} - \nabla_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{x}_{k}, \boldsymbol{p})^{-1} \boldsymbol{F}(\boldsymbol{x}_{k}, \boldsymbol{p})$$
#
# until $\boldsymbol{F}(\boldsymbol{x}_{k}, \boldsymbol{p})$ is sufficiently small.
# Now I should note that here, our age model is linear in $x$ (or `age` in our code), so that the solution will be found in a single iteration, or a sinle "matrix inversion", as could be seen from our steady-state equation for $\boldsymbol{a}$.

#md # !!! note
#md #     AIBECS comes with a built-in algorithm and an API to solve for the steady-state, so you don't have to worry about all these details!
#nb # > **Note**
#nb # > AIBECS comes with a built-in algorithm and an API to solve for the steady-state, so you don't have to worry about all these details!

# ##### Define the Steady-state problem in AIBECS

# First, we create an instance of the steady-state problem, via

prob = SteadyStateProblem(F, ∇ₓF, x₀, p₀)

# where we have simply provided the state function, $\boldsymbol{F}$, the Jacobian, $\nabla_{\boldsymbol{x}}\boldsymbol{F}$, the initial guess and the parameters.
# The `SteadyStateProblem` function is a standard "DiffEqBase" constructor that I have overloaded in my package so that you can easily generate the model here.

# ##### Solve for the steady-state with AIBECS

# Finally, we can find the solution in litterally one line, via the `solve` function:

age = solve(prob, CTKAlg())

# Here, I have provided the `solve` function with two things:
# - the problem, `prob`, which we just defined, and
# - the quasi-Newton algorithm that I wrote in Julia, denoted by `CTKAlg()` after C.T. Kelley, who originally wrote it in MATLAB.
#
# The last line should take about 10 seconds to 1 minute, depending on your laptop.
# That's it!
# We solved for the steady state!
# Everyone here deserves a nice tap on the shoulder — Good job!
# Now let's see what this age looks like on a map

#---------------------------------------------------------
# ## Figures
#---------------------------------------------------------

# We will plot a horizontal slice of the age at about 1000m depth using Cartopy.

# First, we must rearrange `age` into the 3D model grid.
# For that we will need the vector of the indices of wet points in the 3D grid, which we will denote by `iwet`, and which AIBECS generates via

iwet = indices_of_wet_boxes(wet3D)

# We then rearrange the column vector `age` into a 3D array via

age_3D = fill(NaN, size(wet3D)) # creates a 3D array of NaNs of the same size as `wet3D`
age_3D[iwet] = age              # Fills the wet grid boxes with the age values
size(age_3D)                    # Just to check the size of age_3D

# The last line just shows you the size of `age_3D`, which is a 3D-array as expected.

# Now let us find the index of the depth that is closest to 1000m.
# To do that we must use the depth information contained in `grd`.
# Let us first create a small vector of the depths of the grid:

depth = grd.depth

# We could count the index of the entry we want, but here we will use the `findfirst` function to find the first depth index that is greater than 1000m.
#nb # (Feel free to change the value of `iz` if you want to see a slice at another depth.)

iz = findfirst(depth .> 1000u"m")
iz, depth[iz]

# We get `iz = 13`, which is a layer that lies at 1104m, close to 1000m like we wanted.
#
# Finally, we need the latitude and longitudes of the grid, contained in `grd`.
# As for `depth`, we can use the OCIM's `grd` output:

lat, lon = ustrip.(grd.lat), ustrip.(grd.lon)

# So these are the latitudes and longitudes of the map we are about to plot.

# A last thing we can do is convert the age from seconds, `u"s"`, to years, `u"yr"`, because the age is large.
# This can be done via the Unitful package (loaded automatically by AIBECS).

age_3D_1000m_yr = age_3D[:,:,iz] * ustrip(1.0u"s" |> u"yr")

# Finally! Let's have a look at this ideal mean age!
# To make figures, here, we use Cartopy.
#nb # (You should have installed Cartopy if you went through the prerequisites correctly).
# To use it we simply type

ENV["MPLBACKEND"]="qt5agg"
using PyPlot, PyCall

#md # !!! note
#md #     The first line is needed for Mac users.
#md #     It's a bug that should eventually be resolved, but for now this seems to make it work.
#nb # > **Note**
#nb # > The first line is needed for Mac users.
#nb # > It's a bug that should eventually be resolved, but for now this seems to make it work.

# We import cartopy, define a new plot, add some coastlines because they are pretty, and add our slice of age at 1000m depth to it via

clf()
ccrs = pyimport("cartopy.crs")
ax = subplot(projection=ccrs.EqualEarth(central_longitude=-155.0))
ax.coastlines()
lon_cyc = [lon; 360+lon[1]] # making it cyclic for Cartopy
age_cyc = hcat(age_3D_1000m_yr, age_3D_1000m_yr[:,1])
p = contourf(lon_cyc, lat, age_cyc, levels=0:100:1200, transform=ccrs.PlateCarree(), zorder=-1)
colorbar(p, orientation="horizontal")
gcf() # gets the current figure to display

# That's it!
# Good job!
#
# At 1000m, the age ranges from a few years below deep water formation regions (Wedell Sea, North Atlantic), and reaches a dozen of centuries in the North Pacific!
# This is pretty good for so little work!


text:
using StateSpaceReconstruction
using Test
import DelayEmbeddings.Dataset

# Some example data
pts = rand(4, 25)
E = cembed(pts)
D = Dataset(transpose(pts))

# Test constructors
@test DelaunayTriangulation(pts) isa DelaunayTriangulation
@test DelaunayTriangulation(transpose(pts)) isa DelaunayTriangulation
@test DelaunayTriangulation(E) isa DelaunayTriangulation
@test DelaunayTriangulation(D) isa DelaunayTriangulation

# Test indexing
DT = DelaunayTriangulation(pts)
@test DT[1] isa Vector{Int32}
@test length(DT[1]) == 5
@test DT[1:2] isa Vector{Vector{Int32}}
@test length(DT[1:2]) == 2
@test DT[:, 1] isa Array{Int32, 2}
@test size(DT[:, 1]) == (5, 1)

@test DT[:, 1:4] isa Array{Int32, 2}
@test size(DT[:, 1:4]) == (5, 4)
@test DT[:, [1, 2]] isa Array{Int32, 2}
@test DT[:, :] isa Array{Int32, 2}


text:
# Author: Lukasz Bratos

# Funkcja przyjmuje zakres [a, b].
# Jesli w danym zakresie liczby sa rownomiernie rozmieszczone,
# to zwraca roznice pomiedzy dwoma kolejnymi liczbami,
# w przeciwnym przypadku zwraca 0.
function calcualateSpread(a :: Float64, b :: Float64)
    exponentA = SubString(bitstring(a), 2:12)
    exponentB = SubString(bitstring(prevfloat(b)), 2:12)

    if exponentA != exponentB
        return 0.0
    end

    exponent = parse(Int, exponentA, base = 2)

    return (2.0 ^ (exponent - 1023)) * (2.0 ^ (-52))
end

println("Spread for [0.5, 1]: ", calcualateSpread(0.5, 1.0))
println("Spread for [1, 2]: ", calcualateSpread(1.0, 2.0))
println("Spread for [2, 4]: ", calcualateSpread(2.0, 4.0))


text:
 using CUDAnative
device!(0)
using CuArrays
using CuArrays.CUBLAS: libcublas_handle, cublasSetStream_v2, CuDefaultStream
using CUDAnative.CUDAdrv: CuStream
using GPUArrays
using BlockArrays: _BlockArray, PseudoBlockArray, BlockArray, BlockMatrix, BlockVector,
                  nblocks, Block, cumulsizes, AbstractBlockVector
using BlockBandedMatrices: BandedBlockBandedMatrix, _BandedBlockBandedMatrix,
                           blockbandwidths, subblockbandwidths, blockbandwidth,
                           BandedBlockBandedSizes
using LinearAlgebra: BLAS
using BandedMatrices: _BandedMatrix
using SharedArrays
using ArrayLayouts
import Distributed

import Adapt: adapt
import LinearAlgebra

############### Loot and plunder
# BlockArrays
adapt(T::Type{<:AbstractArray}, b::BlockArray) =
    _BlockArray(T.(b.blocks), b.block_sizes)
adapt(T::Type{<:AbstractArray}, b::PseudoBlockArray) =
    PseudoBlockArray(T(b.blocks), b.block_sizes)
adapt(T::Type{<:PseudoBlockArray}, b::BlockArray) = T(b.blocks, b.block_sizes)
adapt(T::Type{<:BlockArray}, b::PseudoBlockArray) = T(b.blocks, b.block_sizes)
# CuArrays and BlockArrays
if @isdefined CuArray
  adapt(T::Type{<:CuArray}, b::PseudoBlockArray) = adapt(T, BlockArray(b))
end
###############

  adapt(T::Type, b::BandedBlockBandedMatrix) =
    _BandedBlockBandedMatrix(adapt(T, b.data), b.block_sizes)

function LinearAlgebra.mul!(c::BlockVector{T},
                            A::BandedBlockBandedMatrix{T, <: BlockMatrix},
                            x::BlockVector{T}) where T
    @assert nblocks(A, 1) == nblocks(c, 1)
    @assert cumulsizes(A, 1) == cumulsizes(c, 1)
    @assert nblocks(A, 2) == nblocks(x, 1)
    @assert cumulsizes(A, 2) == cumulsizes(x, 1)

    for block in c.blocks
      fill!(block, zero(eltype(block)))
    end
    l, u = blockbandwidths(A)
    λ, μ = subblockbandwidths(A)
    N,M = nblocks(A)

    @inbounds for i = 1:N, j = max(1,i-l):min(M,i+u)
        BLAS.gbmv!('N', size(view(A, Block(i, j)), 1), λ, μ, one(T),
                   A.data.blocks[i - j + u + 1, j],
                   x.blocks[j], one(T), c.blocks[i])
    end

    # mostly for symmetry with streamed version
    # synchronizing could be left to user
    synchronize(c)
    c
end

function streamed_mul!(c::BlockVector{T},
                       A::BandedBlockBandedMatrix{T, <: BlockMatrix},
                       x::BlockVector{T};
                       nstreams :: Union{Integer, Nothing} = nothing) where T
    @assert nblocks(A, 1) == nblocks(c, 1)
    @assert cumulsizes(A, 1) == cumulsizes(c, 1)
    @assert nblocks(A, 2) == nblocks(x, 1)
    @assert cumulsizes(A, 2) == cumulsizes(x, 1)

    if nstreams isa Nothing
      nₛ = nblocks(A, 1)
    else
      nₛ = min(nstreams, nblocks(A, 1))
    end
    streams = [CuStream() for u in 1:nₛ]
    for (i, block) in enumerate(c.blocks)
      cublasSetStream_v2(libcublas_handle[], streams[(i - 1) % nₛ + 1])
      fill!(block, zero(eltype(block)))
    end
    l, u = blockbandwidths(A)
    λ, μ = subblockbandwidths(A)
    N,M = nblocks(A)

    @inbounds for i = 1:N, j = max(1,i-l):min(M,i+u)
        cublasSetStream_v2(libcublas_handle[], streams[(i - 1) % nₛ + 1])
        BLAS.gbmv!('N', size(view(A, Block(i, j)), 1), λ, μ, one(T),
                   A.data.blocks[i - j + u + 1, j],
                   x.blocks[j], one(T), c.blocks[i])
    end

    synchronize(c)
    cublasSetStream_v2(libcublas_handle[], CuDefaultStream())
    c
end

function banded_mul!(c::BlockVector{T},
                     A::BandedBlockBandedMatrix{T, <: BlockMatrix},
                     x::AbstractBlockVector{T}) where T
    @assert nblocks(A, 1) == nblocks(c, 1)
    @assert cumulsizes(A, 1) == cumulsizes(c, 1)
    @assert nblocks(A, 2) == nblocks(x, 1)
    @assert cumulsizes(A, 2) == cumulsizes(x, 1)

    for block in c.blocks
      fill!(block, zero(eltype(block)))
    end
    l, u = blockbandwidths(A)
    λ, μ = subblockbandwidths(A)
    N, M = nblocks(A)

    @inbounds for i = 1:N, j = max(1,i-l):min(M,i+u)
      B = _BandedMatrix(A.data.blocks[i - j + u + 1, j],
                        size(view(A, Block(i, j)), 1),
                        λ, μ)
        muladd!(one(T), B, x.blocks[j], one(T), c[Block(i)])
    end

    c
end

using Test

function testme()
  @testset "block-banded on NVIDIA gpus" begin

    @testset "BlockArray Adapters" begin
      bmat = BlockArray{Float64}(undef, [1, 1], [2, 2])
      @test adapt(JLArray, bmat) isa BlockArray{T, 2, JLArray{T, 2}} where T
      @test eltype(adapt(JLArray, bmat)) === Float64
      if @isdefined CuArray
        @test cu(bmat) isa BlockArray{T, 2, CuArray{T, 2}} where T
        @test eltype(cu(bmat)) === Float32
      end
    end

    @testset "PseudoBlockArray Adapters" begin
      bmat = PseudoBlockArray{Float64}(undef, [1, 1], [2, 2])
      @test eltype(adapt(JLArray, bmat)) === Float64
      @test adapt(JLArray, bmat) isa PseudoBlockArray
      if @isdefined CuArray
        @test !(adapt(CuArray, bmat) isa PseudoBlockArray)
        @test adapt(CuArray, bmat) isa BlockArray{T, 2, CuArray{T, 2}} where T
        @test cu(bmat) isa BlockArray{T, 2, CuArray{T, 2}} where T
        @test eltype(cu(bmat)) === Float32
      end
    end

    @testset "PseudoBlockArray Adapters" begin
      bmat = BandedBlockBandedMatrix{Float64}(undef, ([1, 1], [2, 2]), (1, 2), (1, 1))
      @test adapt(JLArray, bmat) isa BandedBlockBandedMatrix
      @test adapt(JLArray, bmat).data isa PseudoBlockArray{T, 2, JLArray{T, 2}} where T
      @test eltype(adapt(JLArray, bmat)) === Float64
      if @isdefined CuArray
        @test adapt(CuArray, bmat).data isa BlockArray{T, 2, CuArray{T, 2}} where T
        @test cu(bmat) isa BandedBlockBandedMatrix
        @test cu(bmat).data isa BlockArray{T, 2, CuArray{T, 2}} where T
        @test eltype(cu(bmat)) === Float32
      end
    end

    @testset "Multiplication" begin
       N, M = rand(1:20, 2)
       l, u, λ, μ = rand(0:2, 4)
       n, m = rand(max(l, u, λ, μ):20, N), rand(max(l, u, λ, μ):20, M)
       A = BandedBlockBandedMatrix{Float64}(undef, (n, m), (l, u), (λ, μ))
       A.data .= rand.()
       Ablock = adapt(BlockArray, A)
       cblock = BlockArray(Array{Float64, 1}(undef, size(A, 1)), n)
       cblock .= rand.()
       x = PseudoBlockArray(Array{Float64, 1}(undef, size(A, 2)), m)
       x .= rand.()
       xblock = adapt(BlockArray, x)

       @test LinearAlgebra.mul!(cblock, Ablock, xblock) ≈ A * x
       cblock .= 0
       @test banded_mul!(cblock, Ablock, xblock) ≈ A * x
       cgpu = cu(cblock)
       @test streamed_mul!(cgpu, cu(A), cu(x)) ≈ A * x
    end

  end
end

using BenchmarkTools
using Statistics

function benchmarks()
  suite = BenchmarkGroup()
  # suite["viabm"] = BenchmarkGroup()
  suite["pseudo"] = BenchmarkGroup()
  suite["block"] = BenchmarkGroup()
  if @isdefined CuArrays
    suite["gpu"] = BenchmarkGroup()
    suite["streams"] = BenchmarkGroup()
  end
  for N in [10, 100, 500, 1000], n = [N ÷ 5, N, 5N, 10N]
    l, u, λ, μ = rand(0:2, 4)
    M, m = N, n

    A = BandedBlockBandedMatrix{Float64}(
             undef, (repeat([n], N), repeat([m], M)), (l, u), (λ, μ))
    A.data .= rand.()
    c = PseudoBlockArray(Array{Float64, 1}(undef, size(A, 1)), repeat([n], N))
    c .= rand.()
    x = PseudoBlockArray(Array{Float64, 1}(undef, size(A, 2)), repeat([m], M))
    x .= rand.()

    suite["pseudo"]["N=$N n=$n"] = @benchmarkable begin
      $c .= Mul($A, $x)
    end
    suite["block"]["N=$N n=$n"] = @benchmarkable begin
      LinearAlgebra.mul!($(adapt(BlockArray, c)), $(adapt(BlockArray, A)),
                         $(adapt(BlockArray, x)))
    end
    if @isdefined CuArrays
      gpus = Dict(:c => adapt(CuArray, c),
                  :x => adapt(CuArray, x),
                  :A => adapt(CuArray, A))

      suite["gpu"]["N=$N n=$n"] = @benchmarkable begin
        LinearAlgebra.mul!($(gpus[:c]), $(gpus[:A]), $(gpus[:x]))
      end
      suite["streams"]["N=$N n=$n"] = @benchmarkable begin
        gpuc = streamed_mul!($(gpus[:c]), $(gpus[:A]), $(gpus[:x]))
      end
    end
  end
  suite
end

block_ratio(result, name; method=median) =
    ratio(method(result["block"][name]), method(result["pseudo"][name]))
viabm_ratio(result, name; method=median) =
    ratio(method(result["viabm"][name]), method(result["block"][name]))


text:
using Base.Test

include("apitest.jl")
include("testexamples.jl")
include("testshow.jl")
include("test_mathprogbase.jl")


text:
using LinearAlgebra
using Random
using Statistics
using DataFrames
using FreqTables

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 1
#:::::::::::::::::::::::::::::::::::::::::::::::::::
using Optim
f(x) = -x[1]^4-10x[1]^3-2x[1]^2-3x[1]-2
minusf(x) = x[1]^4+10x[1]^3+2x[1]^2+3x[1]+2
startval = rand(1)   # random starting value
result = optimize(minusf, startval, BFGS())

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 2
#:::::::::::::::::::::::::::::::::::::::::::::::::::
using DataFrames
using CSV
using HTTP
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS1-julia-intro/nlsw88.csv"
df = CSV.read(HTTP.get(url).body,DataFrame)
X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.married.==1

function ols(beta, X, y)
    ssr = (y.-X*beta)'*(y.-X*beta)
    return ssr
end

beta_hat_ols = optimize(b -> ols(b, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
println(beta_hat_ols.minimizer)

using GLM
bols = inv(X'*X)*X'*y
df.white = df.race.==1
bols_lm = lm(@formula(married ~ age + white + collgrad), df)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 3
#:::::::::::::::::::::::::::::::::::::::::::::::::::
function logit(alpha, X, d)
    loglike = sum([d[i]*log(exp(sum(alpha.*X[i,:]))/(1+exp(sum(alpha.*X[i,:]))))+(1-d[i])*log(1/(1+exp(sum(alpha.*X[i,:])))) for i=1:size(X,1)])
    return loglike
end

startval = rand(size(X,2))
obj = optimize(alpha->-logit(alpha, X, y), startval, LBFGS())
println(obj.minimizer)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 4
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# see Lecture 3 slides for example
res = glm(X, y, Binomial(), LogitLink())
println(res)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 5
#:::::::::::::::::::::::::::::::::::::::::::::::::::
freqtable(df, :occupation) # note small number of obs in some occupations
df = dropmissing(df, :occupation)
df[df.occupation.==8 ,:occupation] .= 7
df[df.occupation.==9 ,:occupation] .= 7
df[df.occupation.==10,:occupation] .= 7
df[df.occupation.==11,:occupation] .= 7
df[df.occupation.==12,:occupation] .= 7
df[df.occupation.==13,:occupation] .= 7
freqtable(df, :occupation) # problem solved

X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.occupation

function mlogit(alpha, X, d)
    loglike = 0
    for i = 1:size(X,1)
        term1 = sum([convert(Float64,d[i]==j)*sum(alpha[:,j].*X[i,:]) for j =1:6])
        term2 = log(1 + sum([exp(sum(alpha[:,j].*X[i,:])) for j =1:6]))
        loglike += term1 - term2
    end
    return loglike
end

#startval = zeros(size(X,2),6)
#startval = rand(size(X,2),6)
startval = 2*rand(size(X,2),6).-1
obj = optimize(alpha->-mlogit(alpha, X, y), startval, LBFGS(), Optim.Options(g_tol=1e-5, iterations=100000))
println(obj.minimizer)


text:
__precompile__()

module ResettableStacks

  const FULL_RESET_COUNT = 10000

  using StaticArrays

  import Base: isempty, length, push!, pop!, iterate

  include("core.jl")

  export ResettableStack, reset!
end # module


text:
const CEG = SumOfSquares.Certificate.ChordalExtensionGraph

"""
    struct MonomialSparsity{C<:CEG.AbstractCompletion} <: Sparsity
        completion::C
        k::Int
        use_all_monomials::Bool
    end

Monomial or term sparsity as developed in [WML20a, WML20b].
The `completion` field should be `ClusterCompletion()` [default] for the block-closure or cluster completion [WML20a],
and `ChordalCompletion()` for chordal completion [WML20b].
The integer `k` [default=0] corresponds to `Σ_k` defined in [(3.2), WML20a]
and `k = 0` corresponds to `Σ_*` defined in [(3.3), WML20a].
If `use_all_monomials` is `false` then some monomials of the basis
might be dropped from the basis if not needed.

[WML20a] Wang, Jie, Victor Magron, and Jean-Bernard Lasserre.
*TSSOS: A Moment-SOS hierarchy that exploits term sparsity*.
arXiv preprint arXiv:1912.08899 (2020).

[WML20b] Wang, Jie, Victor Magron, and Jean-Bernard Lasserre.
*Chordal-TSSOS: a moment-SOS hierarchy that exploits term sparsity with chordal extension*.
arXiv preprint arXiv:2003.03210 (2020).
"""
struct MonomialSparsity{C<:CEG.AbstractCompletion} <: Sparsity
    completion::C
    k::Int
    use_all_monomials::Bool
    function MonomialSparsity(
        completion::CEG.AbstractCompletion=CEG.ClusterCompletion(),
        k::Int=0,
        use_all_monomials::Bool=false,
    )
        return new{typeof(completion)}(completion, k, use_all_monomials)
    end
end

# Note: For some implementation of MultivariatePolynomials such as
# https://github.com/blegat/CondensedMatterSOS.jl,
# the product of monomials may be a term so wrap any multiplication of monomials by
# `MP.monomial`.

const MP = SumOfSquares.MP
function monomial_sparsity_graph(monos, P, use_all_monomials::Bool)
    g = CEG.LabelledGraph{eltype(monos)}()
    if use_all_monomials
        squares = nothing # No need to record this
    else
        squares = Set{eltype(monos)}()
    end
    for mono in monos
        CEG.add_node!(g, mono)
    end
    for a in monos
        for b in monos
            if MP.monomial(a * b) in P
                if a != b
                    CEG.add_edge!(g, a, b)
                elseif squares !== nothing
                    push!(squares, a)
                end
            end
        end
    end
    return g, squares
end
function _add_monos(add, monos, H, squares)
    for a in monos
        neighbors = H.graph.neighbors[H.n2int[a]]
        # If `!isempty!(neighbors)`, the monomial part of a clique with another
        # monomial so the diagonal entry should be nonzero.
        # If `a in squares`, the monomial square is needed so it is added.
        # e.g. if the polynomial is x^2 but the basis is `[x, 1]`,
        # `squares` is `Set([x])` so `x^2` will be added but not `1^2`.
        if squares === nothing || !isempty(neighbors) || a in squares
            add(a^2)
        end
        for bi in neighbors
            b = H.int2n[bi]
            add(MP.monomial(a * b))
        end
    end
end
function completion_with_squares(g, squares, completion)
    H, cliques = CEG.completion(g, completion)
    if squares !== nothing
        cliques = filter(monos -> length(monos) > 1 || (monos[1] in squares), collect(cliques))
    end
    return H, cliques
end
function monomial_sparsity_iteration(P, completion, use_all_monomials::Bool, monos)
    P_next = Set{eltype(P)}()
    g, squares = monomial_sparsity_graph(monos, P, use_all_monomials)
    H, cliques = completion_with_squares(g, squares, completion)
    _add_monos(mono -> push!(P_next, mono), monos, H, squares)
    return P_next, cliques
end
struct GeneratorP{PT, GT}
    P::PT
    generator_monos::GT
end
function Base.in(mono, g::GeneratorP)
    return any(g.generator_monos) do g_mono
        MP.monomial(mono * g_mono) in g.P
    end
end
function monomial_sparsity_iteration(P, completion, use_all_monomials::Bool, monos, multiplier_generator_monos)
    P_next, cliques = monomial_sparsity_iteration(P, completion, use_all_monomials, monos)
    multiplier_cliques = map(multiplier_generator_monos) do m
        multiplier_monos, generator_monos = m
        g, squares = monomial_sparsity_graph(multiplier_monos, GeneratorP(P, generator_monos), use_all_monomials)
        H, _cliques = completion_with_squares(g, squares, completion)
        _add_monos(multiplier_monos, H, squares) do mono
            for b in generator_monos
                push!(P_next, MP.monomial(mono * b))
            end
        end
        return _cliques
    end
    return P_next, (cliques, multiplier_cliques)
end
_monovec(cliques::AbstractVector{<:MP.AbstractMonomial}) = MP.monovec(cliques)
_monovec(cliques) = _monovec.(cliques)
function sparsity(monos::AbstractVector{<:MP.AbstractMonomial}, sp::MonomialSparsity,
                  gram_monos::AbstractVector = SumOfSquares.Certificate.monomials_half_newton_polytope(monos, tuple()),
                  args...)
    P = Set(monos)
    if sp.use_all_monomials
        for mono in gram_monos
            push!(P, mono^2)
        end
    end
    cliques = nothing
    iter = 0
    while iter < sp.k || iszero(sp.k)
        P_prev = P
        P, cliques = monomial_sparsity_iteration(P_prev, sp.completion, sp.use_all_monomials, gram_monos, args...)
        if iszero(iter)
            # If gram_monos + gram_monos !⊆ monos, then it's possible that P_prev !⊆ P
            P == P_prev && break
        else
            length(P) >= length(P_prev) || error("Set of monomials should be increasing in monomial sparsity iterations.")
            length(P) == length(P_prev) && break
        end
        iter += 1
    end
    return _monovec(cliques)
end
# This also checks that it is indeed a monomial basis
_monos(basis::MB.MonomialBasis) = basis.monomials
function _gram_monos(vars, certificate::MaxDegree{CT, MB.MonomialBasis}) where CT
    return _monos(maxdegree_gram_basis(MB.MonomialBasis, vars, certificate.maxdegree))
end
function sparsity(poly::MP.AbstractPolynomial, domain::BasicSemialgebraicSet, sp::MonomialSparsity, certificate::AbstractPreorderCertificate)
    gram_monos = _gram_monos(
        reduce((v, q) -> unique!(sort!([v..., variables(q)...], rev=true)),
                  domain.p, init = variables(poly)),
        get(certificate, IdealCertificate())
    )
    processed = get(certificate, PreprocessedDomain(), domain, poly)
    multiplier_generator_monos = [
        (_monos(get(certificate, MultiplierBasis(), index, processed)),
         monomials(get(certificate, Generator(), index, processed)))
        for index in get(certificate, PreorderIndices(), processed)
    ]
    cliques, multiplier_cliques = sparsity(monomials(poly), sp, gram_monos, multiplier_generator_monos)
    return MB.MonomialBasis.(cliques), [MB.MonomialBasis.(clique) for clique in multiplier_cliques]
end