text:
# This example demonstrates the use of the medium-level R-ncdf-like interface
#
# First of all we create an array with top-of the atmosphere radiation data

using NetCDF
include("toa.jl")

# Define longitudes and latitudes, day and timesteps
lat=[-89:89]
lon=[0:359]
day=1
tim=[0:23]

# Create radiation array
rad = float64([g_pot(x2,x1,day,x3) for x1=lon, x2=lat, x3=tim])

# Define some attributes of the variable (optionlal)
varatts = @Compat.Dict("longname" => "Radiation at the top of the atmosphere",
           "units"    => "W/m^2")
lonatts = @Compat.Dict("longname" => "Longitude",
           "units"    => "degrees east")
latatts = @Compat.Dict("longname" => "Latitude",
           "units"    => "degrees north")
timatts = @Compat.Dict("longname" => "Time",
           "units"    => "hours since 01-01-2000 00:00:00")

#Here we start by defining the dimensions. This is done by creating NcDim objects:

latdim = NcDim("lat",lat,latatts)
londim = NcDim("lon",lon,lonatts)
timdim = NcDim("time",tim,timatts)

# Then we create an NcVar object, the data type is defined by the corresponding julia type:

radvar = NcVar("rad",[londim,latdim,timdim],varatts,Float32)

# Now we can finally create the netcdf-file and get a file handler in return:

isfile("radiation2.nc") ? rm("radiation2.nc") : nothing
nc = NetCDF.create("radiation2.nc",radvar)

# Writing data to the file is done using putvar

NetCDF.putvar(nc, "rad", rad )

# And we close the file

NetCDF.close( nc )

# We open a file for reading:

nc=NetCDF.open("radiation2.nc")

# we now have an NcFile object that can be browsed

println("Names of the NcFile fields: \n",names(nc))

println("Attributes of the variable rad: ",nc.vars["rad"].atts)

# we read the whole variable:

x = NetCDF.readvar(nc,"rad",[1,1,1],[-1,-1,-1])
println("Successfully read an array of size ",size(x))

# Additional
# Reading parts of a file, for example reading only time steps 5 and 6 can be done with:

x=NetCDF.readvar(nc,"rad",[1,1,5],[-1,-1,2])
println("Successfully read an array of size ",size(x))

#here the first array [1,1,5] gives the starting position for reading while the second array [1,1,2] gives the number of blocks to be read along each dimension.
# If the count is set to -1 the whole dimension is read.
NetCDF.close(nc)


text:
################################################################################
# NullSpace
################################################################################

mutable struct NullSpace{T,SVv,SVΔ,SVλ}
	probsize::ProblemSize
	mat::Matrix{T}
	vec::Vector{SVv}
	Δtraj::Vector{SVΔ}
	Δλ::Vector{SVλ}
end

function NullSpace(probsize)
	mat = zeros(0,0)
	vec = Vector([zeros(0)])
	Δtraj = Vector([zeros(0)])
	Δλ = Vector([zeros(0)])
	return NullSpace(probsize, mat, vec, Δtraj, Δλ)
end

function reset!(null::NullSpace{T,SVv,SVΔ,SVλ}) where {T,SVv,SVΔ,SVλ}
	null.mat = zeros(T,0,0)
	null.vec = Vector{SVv}()
	null.Δtraj = Vector{SVΔ}()
	null.Δλ = Vector{SVλ}()
	return nothing
end

function add_matrix!(null::NullSpace, mat::AbstractMatrix, hmask::Vector{Int})
	probsize = null.probsize
	N = probsize.N
	p = probsize.p
	S = probsize.S
	Sh = S + p*(p-1)*(N-1)
	m1, m2 = size(mat)
	@assert (S <= m1 <= Sh)

	null.mat = mat
	for l = 1:m2
		vec = zeros(Sh)
		vec[hmask] = mat[:,l]
		vec ./= mean(abs.(vec))
		Δtraj = vec[1:S]
		Δλ = vec[S+1:end]
		push!(null.vec, vec)
		push!(null.Δtraj, Δtraj)
		push!(null.Δλ, Δλ)
	end
	return nothing
end

################################################################################
# ActiveSetCore
################################################################################

mutable struct ActiveSetCore{Vr,Vrt,SMj,SVhi,SVvi,SAr,SAj,Sd}
	res::Vr               # residual vector
	res_tmp::Vrt          # holds a temporary copy of the residual vector
	jac::SMj              # residual sparse jacobian
	probsize::ProblemSize # size of the problem
	horiz_inds::SVhi      # indices for each variable in an horizontal block
	verti_inds::SVvi      # indices for each variable in an vertical block
	res_sub::SAr          # Residual views dictionary
	jac_sub::SAj          # Jacobian views dictionary
	dyn::Sd               # Dictionary of the Jacobian dynamics indices
	vmask::Vector{Int}    # Vertical mask of the active inequality constraints
	hmask::Vector{Int}    # Horizontal mask of the active inequality constraints
	null::NullSpace       # Structure holding nullspace of the active set Jacobian
end

function ActiveSetCore(probsize::ProblemSize)
	N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	Sv = probsize.S + Int(p*(p-1)*(N-1)/2)
	Sh = probsize.S + p*(p-1)*(N-1)
	res = zeros(Sv)
	res_tmp = deepcopy(res)
	jac = spzeros(Sv,Sh)
	verti_inds = complete_vertical_indices(probsize)
	horiz_inds = complete_horizontal_indices(probsize)
	res_sub = complete_residual_views(res, probsize, verti_inds)
	jac_sub = complete_jacobian_views(jac, probsize, verti_inds, horiz_inds)
	dyn = dynamics_indices(probsize)
	vmask = Vector{Int}(1:Sv)
	hmask = Vector{Int}(1:Sh)
	null = NullSpace(probsize)
	TYPE = typeof.((res, res_tmp, jac, horiz_inds, verti_inds, res_sub, jac_sub, dyn))
	return ActiveSetCore{TYPE...}(res, res_tmp, jac, probsize,
		horiz_inds, verti_inds, res_sub, jac_sub, dyn, vmask, hmask, null)
end

################################################################################
# Indices
################################################################################

function complete_vertical_indices(probsize::ProblemSize)
	N = probsize.N
	n = probsize.n
	p = probsize.p
	mi = probsize.mi
	verti_inds = Dict()
	off = 0
	for i = 1:p
		for k = 1:N-1
			stamp = stampify(:opt, i, :x, 1, k+1)
			verti_inds[stamp] = SVector{n,Int}(off .+ (1:n))
			off += n
			stamp = stampify(:opt, i, :u, i, k)
			verti_inds[stamp] = SVector{mi[i],Int}(off .+ (1:mi[i]))
			off += mi[i]
		end
	end
	for k = 1:N-1
		stamp = stampify(:dyn, 1, :x, 1, k)
		verti_inds[stamp] = SVector{n,Int}(off .+ (1:n))
		off += n
	end
	for k = 2:N
		for i = 1:p
			for j = i+1:p
				stamp = stampify(:v, :col, i, j, k)
				verti_inds[stamp] = SVector{1,Int}(off .+ (1:1))
				off += 1
			end
		end
	end
	return verti_inds
end

function complete_horizontal_indices(probsize::ProblemSize)
	N = probsize.N
	n = probsize.n
	p = probsize.p
	mi = probsize.mi

	horiz_inds = Dict()
	off = 0
	for k = 1:N-1
		stamp = stampify(:x, 1, k+1)
		horiz_inds[stamp] = SVector{n,Int}(off .+ (1:n))
		off += n
		for i = 1:p
			stamp = stampify(:u, i, k)
			horiz_inds[stamp] = SVector{mi[i],Int}(off .+ (1:mi[i]))
			off += mi[i]
		end
		for i = 1:p
			stamp = stampify(:λ, i, k)
			horiz_inds[stamp] = SVector{n,Int}(off .+ (1:n))
			off += n
		end
	end
	for k = 2:N
		for i = 1:p
			for j ∈ setdiff(1:p,i)
				stamp = stampify(:h, :col, i, j, k)
				horiz_inds[stamp] = SVector{1,Int}(off .+ (1:1))
				off += 1
			end
		end
	end
	return horiz_inds
end

function horizontal_idx(core::ActiveSetCore, stamp::CStamp)
	return core.horiz_inds[stamp]
end

function horizontal_idx(horiz_inds::Dict, stamp::CStamp)
	return horiz_inds[stamp]
end

function vertical_idx(core::ActiveSetCore, stamp::CStamp)
	return core.verti_inds[stamp]
end

function vertical_idx(verti_inds::Dict, stamp::CStamp)
	return verti_inds[stamp]
end

################################################################################
# SubArrays
################################################################################

function complete_residual_views(res::SV, probsize::ProblemSize, verti_inds::Dict) where {SV}
	N = probsize.N
	p = probsize.p

	# Residual Views
	res_sub = Dict{AbstractStamp,SubArray}()
	for prob ∈ (:dyn, :opt)
		for i0 = 1:p
			for n1 in (:x,:u)
				for i1 = 1:p
					for v1 = 1:N
						stamp = stampify(prob, i0, n1, i1, v1)
						if valid(stamp,N,p)
							res_sub[stamp] = view(res, vertical_idx(verti_inds, stamp))
						end
					end
				end
			end
		end
	end
	for k = 2:N
		for i = 1:p
			for j = i+1:p
				stamp = stampify(:v,:col, i, j, k)
				if valid(stamp,N,p)
					res_sub[stamp] = view(res, vertical_idx(verti_inds, stamp))
				end
			end
		end
	end
	return res_sub
end

function complete_jacobian_views(jac::SM, probsize::ProblemSize, verti_inds::Dict, horiz_inds::Dict) where {SM}
	N = probsize.N
	p = probsize.p

	# Jacobian Views
	jac_sub = Dict{Any,SubArray}()

	for prob ∈ (:opt,)
		for i0 = 1:p
			for n1 in (:x,)
				for i1 = 1:1
					for v1 = 1:N
						for j = 1:p
							vstamp = stampify(prob, i0, n1, i1, v1)
							hstamp = stampify(:h, :col, i0, j, v1)
							if valid(vstamp,N,p) && valid(hstamp,N,p)
								jac_sub[(vstamp, hstamp)] = view(jac,
									vertical_idx(verti_inds, vstamp),
									horizontal_idx(horiz_inds, hstamp))
							end
						end
					end
				end
			end
		end
	end
	for n2 in (:x,)
		for i2 = 1:1
			for v2 = 2:N
				for j = i2+1:p
					vstamp = stampify(:v, :col, i2, j, v2)
					hstamp = stampify(n2, i2, v2)
					if valid(vstamp,N,p) && valid(hstamp,N,p)
						jac_sub[(vstamp, hstamp)] = view(jac,
							vertical_idx(verti_inds, vstamp),
							horizontal_idx(horiz_inds, hstamp))
					end
				end
			end
		end
	end
	return jac_sub
end

ext: jl


text:
"""
```
imgb = binarize(Otsu(), img)
binarize!(Otsu(), img)
```

Under the assumption that the histogram is bimodal the threshold is
set so that the resultant inter-class variance is maximal.

# Arguments

The function arguments are described in more detail below.

##  `img`

An `AbstractGray` image which is to be binarized into background (`Gray(0)`)
and foreground (`Gray(1)`) pixels.

## Reference

Nobuyuki Otsu (1979). "A threshold selection method from gray-level histograms". *IEEE Trans. Sys., Man., Cyber.* 9 (1): 62–66. [doi:10.1109/TSMC.1979.4310076](http://dx.doi.org/doi:10.1109/TSMC.1979.4310076)
"""
function binarize(algorithm::Otsu,  img::AbstractArray{T,2}) where T <: AbstractGray
  binarize!(algorithm, copy(img))
end

"""
```
t = binarize!(Otsu(), img)
```
Same as [`binarize`](@ref binarize(::Otsu, ::AbstractArray{<:AbstractGray,2})) except that it modifies the image that was passed as an argument.
"""
function binarize!(algorithm::Otsu,  img::AbstractArray{T,2}) where T <: AbstractGray
  edges, counts = build_histogram(img,  256)
  t = find_threshold(HistogramThresholding.Otsu(), counts[1:end], edges)
  map!(img,img) do val
    val < t ? Gray(0) : Gray(1)
  end
end


text:
function autocorr(x::Vector, lag::Integer = 30)
    ## plot autocorrelation using Winston
    acVals = ones(lag)
    for ii=1:lag
        acVals[ii] = cor(x[1:(end-ii)], x[(1+ii):end])
    end
    Winston.plot(acVals, "-b")
    Winston.oplot(acVals, ".b")
    Winston.ylim(-0.1, 1)
    ## return acVals
end

function autocorr!(x::Vector, lag::Integer = 30)
    ## add autocorrelation plot to existing plot in blue
    acVals = ones(lag)
    for ii=1:lag
        acVals[ii] = cor(x[1:(end-ii)], x[(1+ii):end])
    end
    Winston.oplot(acVals, "-b")
    Winston.oplot(acVals, ".b")
end

function autocorr_red!(x::Vector, lag::Integer = 30)
    ## add autocorrelation plot to existing plot in red
    acVals = ones(lag)
    for ii=1:lag
        acVals[ii] = cor(x[1:(end-ii)], x[(1+ii):end])
    end
    Winston.oplot(acVals, "-r")
    Winston.oplot(acVals, ".r")
end


text:
# Tests for DirichletMultinomial

using Distributions
using Test, Random, SpecialFunctions

import SpecialFunctions: factorial
Random.seed!(123)

rng = MersenneTwister(123)

@testset "Testing DirichletMultinomial with $key" for (key, func) in
    Dict("rand(...)" => [rand, rand],
         "rand(rng, ...)" => [dist -> rand(rng, dist), (dist, n) -> rand(rng, dist, n)])

# Test Constructors
d = DirichletMultinomial(10, ones(5))
d2 = DirichletMultinomial(10, ones(Int, 5))
@test typeof(d) == typeof(d2)
@test d == deepcopy(d)

α = func[1](5)
d = DirichletMultinomial(10, α)

# test parameters
@test length(d) == 5
@test ntrials(d) == 10
@test params(d) == (10, α)
@test d.α == α
@test d.α0 == sum(α)
@test partype(d) == eltype(α)

# test statistics
@test mean(d) ≈ α * (d.n / d.α0)
p = d.α / d.α0
@test var(d)  ≈ d.n * (d.n + d.α0) / (1 + d.α0) .* p .* (1.0 .- p)
x = func[2](d, 10_000)

# test statistics with mle fit
d = fit(DirichletMultinomial, x)
@test isapprox(mean(d), vec(mean(x, dims=2)), atol=.5)
@test isapprox(var(d) , vec(var(x, dims=2)) , atol=.5)
@test isapprox(cov(d) , cov(x, dims=2)      , atol=.5)

# test Evaluation
d = DirichletMultinomial(10, 5)
@test typeof(d) == DirichletMultinomial{Float64}
@test !insupport(d, func[1](5))
@test insupport(d, [2, 2, 2, 2, 2])
@test insupport(d, 2.0 * ones(5))
@test !insupport(d, 3.0 * ones(5))

for x in (2 * ones(5), [1, 2, 3, 4, 0], [3.0, 0.0, 3.0, 0.0, 4.0], [0, 0, 0, 0, 10])
    @test pdf(d, x) ≈
        factorial(d.n) * gamma(d.α0) / gamma(d.n + d.α0) * prod(gamma.(d.α + x) ./ factorial.(x) ./ gamma.(d.α))
    @test logpdf(d, x) ≈
        log(factorial(d.n)) + loggamma(d.α0) - loggamma(d.n + d.α0) + sum(loggamma, d.α + x) - sum(loggamma, d.α) - sum(log.(factorial.(x)))
end

# test Sampling
x = func[1](d)
@test isa(x, Vector{Int})
@test sum(x) == d.n
@test length(x) == length(d)
@test insupport(d, x)

x = func[2](d, 50)
@test all(x -> (x >= 0), x)
@test size(x, 1) == length(d)
@test size(x, 2) == 50
@test all(sum(x, dims=1) .== ntrials(d))
@test all(insupport(d, x))

# test MLE
x = func[2](d, 10_000)
ss = suffstats(DirichletMultinomial, x)
@test size(ss.s, 1) == length(d)
@test size(ss.s, 2) == ntrials(d)
mle = fit(DirichletMultinomial, x)
@test isapprox(mle.α, d.α, atol=.2)
mle = fit(DirichletMultinomial{Float64}, x)
@test isapprox(mle.α, d.α, atol=.2)

# test MLE with weights
for w in (.1 * ones(10_000), ones(10_000), 10 * ones(10_000))
    mle2 = fit(DirichletMultinomial, x, w)
    @test mle.α ≈ mle2.α
end

end


text:
import Base.convert
using OpenMDAO
import OpenMDAO: detect_apply_nonlinear, detect_guess_nonlinear, detect_apply_linear
using CCBlade: solve

struct CCBladeResidualComp <: OpenMDAO.AbstractImplicitComp
    num_nodes
    num_radial
    af
    B
    turbine
    rotors
    sections
    ops
end

function CCBladeResidualComp(; num_nodes, num_radial, af, B, turbine)
    # Check if the airfoil interpolation passed is a num_radial-length array.
    try
        num_af = length(af)
        if num_af == num_radial
            af = reshape(af, 1, num_radial)
        else
            throw(DomainError("af has length $num_af, but should have length $num_radial"))
        end
    catch e
        if isa(e, MethodError)
            # af is not an array of stuff, so assume it's just a single
            # function, and make it have shape (1, num_radial).
            af = fill(af, 1, num_radial)
        else
            # Some other error happened, so rethrow it.
            rethrow(e)
        end
    end

    Rhub = fill(0., num_nodes, 1)
    Rtip = fill(1., num_nodes, 1)
    pitch = fill(0., num_nodes, 1)
    precone = fill(0., num_nodes, 1)
    rotors = Rotor.(Rhub, Rtip, B, turbine, pitch, precone)

    r = fill(1., num_nodes, num_radial)
    chord = fill(1., num_nodes, num_radial)
    theta = fill(1., num_nodes, num_radial)
    sections = Section.(r, chord, theta, af)

    Vx = fill(1., num_nodes, num_radial)
    Vy = fill(1., num_nodes, num_radial)
    rho = fill(1., num_nodes, 1)
    mu = fill(1., num_nodes, 1)
    asound = fill(1., num_nodes, 1)
    ops = OperatingPoint.(Vx, Vy, rho, mu, asound)

    return CCBladeResidualComp(num_nodes, num_radial, af, B, turbine, rotors, sections, ops)
end

function OpenMDAO.setup(self::CCBladeResidualComp)
    num_nodes = self.num_nodes
    num_radial = self.num_radial

    input_data = [
        VarData("r", shape=(num_nodes, num_radial), val=1., units="m"),
        VarData("chord", shape=(num_nodes, num_radial), val=1., units="m"),
        VarData("theta", shape=(num_nodes, num_radial), val=1., units="rad"),
        VarData("Vx", shape=(num_nodes, num_radial), val=1., units="m/s"),
        VarData("Vy", shape=(num_nodes, num_radial), val=1., units="m/s"),
        VarData("rho", shape=(num_nodes, 1), val=1., units="kg/m**3"),
        VarData("mu", shape=(num_nodes, 1), val=1., units="N/m**2*s"),
        VarData("asound", shape=(num_nodes, 1), val=1., units="m/s"),
        VarData("Rhub", shape=(num_nodes, 1), val=1., units="m"),
        VarData("Rtip", shape=(num_nodes, 1), val=1., units="m"),
        VarData("pitch", shape=(num_nodes, 1), val=1., units="rad"),
        VarData("precone", shape=(num_nodes, 1), val=1., units="rad")]

    output_data = [
        VarData("phi", shape=(num_nodes, num_radial), val=1., units="rad"),
        VarData("Np", shape=(num_nodes, num_radial), val=1., units="N/m"),
        VarData("Tp", shape=(num_nodes, num_radial), val=1., units="N/m"),
        VarData("a", shape=(num_nodes, num_radial), val=1.),
        VarData("ap", shape=(num_nodes, num_radial), val=1.),
        VarData("u", shape=(num_nodes, num_radial), val=1., units="m/s"),
        VarData("v", shape=(num_nodes, num_radial), val=1., units="m/s"),
        VarData("W", shape=(num_nodes, num_radial), val=1., units="m/s"),
        VarData("cl", shape=(num_nodes, num_radial), val=1.),
        VarData("cd", shape=(num_nodes, num_radial), val=1.),
        VarData("cn", shape=(num_nodes, num_radial), val=1.),
        VarData("ct", shape=(num_nodes, num_radial), val=1.),
        VarData("F", shape=(num_nodes, num_radial), val=1.),
        VarData("G", shape=(num_nodes, num_radial), val=1.)]

    partials_data = Vector{PartialsData}()

    of_names = ["phi", "Np", "Tp", "a", "ap", "u", "v", "W", "cl", "cd", "cn", "ct", "F", "G"]

    ss_sizes = Dict(:i=>num_nodes, :j=>num_radial)
    rows, cols = get_rows_cols(ss_sizes=ss_sizes, of_ss=[:i, :j], wrt_ss=[:i])
    for name in of_names
        push!(partials_data, PartialsData(name, "Rhub", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "Rtip", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "pitch", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "precone", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "rho", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "mu", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "asound", rows=rows, cols=cols))
    end

    rows, cols = get_rows_cols(ss_sizes=ss_sizes, of_ss=[:i, :j], wrt_ss=[:i, :j])
    for name in of_names
        push!(partials_data, PartialsData(name, "r", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "chord", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "theta", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "Vx", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "Vy", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "phi", rows=rows, cols=cols))
    end

    for name in ["Np", "Tp", "a", "ap", "u", "v", "W", "cl", "cd", "cn", "ct", "F", "G"]
        push!(partials_data, PartialsData(name, name, rows=rows, cols=cols, val=1.0))
    end

    return input_data, output_data, partials_data
end

function OpenMDAO.linearize!(self::CCBladeResidualComp, inputs, outputs, partials)
    # Copy the data in `inputs` to the array of structs in `self`.
    set_inputs!(self, inputs)

    # Phi, the implicit variable.
    phis = outputs["phi"]

    num_nodes = self.num_nodes
    num_radial = self.num_radial

    # Get the derivatives of the residual. This will have size (num_nodes,
    # num_radial).
    residual_derivs = residual_partials.(phis, self.rotors, self.sections, self.ops)

    # Copy the derivatives of the residual to the partials dict. First do the
    # derivative of the phi residual wrt phi.
    deriv = transpose(reshape(partials["phi", "phi"], (num_radial, num_nodes)))
    deriv .= residual_derivs.phi

    # Copy the derivatives of the residual wrt each input.
    for str in keys(inputs)
        sym = Symbol(str)
        deriv = transpose(reshape(partials["phi", str], (num_radial, num_nodes)))
        @. deriv = getfield(residual_derivs, sym)
    end

    # Get the derivatives of the explicit outputs.
    output_derivs = output_partials.(phis, self.rotors, self.sections, self.ops)

    # Copy the derivatives of the explicit outputs into the partials dict.
    for of_str in keys(outputs)
        if of_str == "phi"
            continue
        else
            of_sym = Symbol(of_str)
            for wrt_str in keys(inputs)
                wrt_sym = Symbol(wrt_str)
                # reshape does not copy data: https://github.com/JuliaLang/julia/issues/112
                deriv = transpose(reshape(partials[of_str, wrt_str], (num_radial, num_nodes)))
                @. deriv = -getfield(getfield(output_derivs, of_sym), wrt_sym)
            end
            # Also need the derivative of each output with respect to phi, but
            # phi is an output, not an input, so we'll need to handle that
            # seperately.
            wrt_str = "phi"
            wrt_sym = Symbol(wrt_str)
            # reshape does not copy data: https://github.com/JuliaLang/julia/issues/112
            deriv = transpose(reshape(partials[of_str, wrt_str], (num_radial, num_nodes)))
            @. deriv = -getfield(getfield(output_derivs, of_sym), wrt_sym)
        end
    end

    return nothing
end

function OpenMDAO.solve_nonlinear!(self::CCBladeResidualComp, inputs, outputs)
    # Copy the data in `inputs` to the array of structs in `self`.
    set_inputs!(self, inputs)

    # When called this way, CCBlade.solve will return a 2D array of `Outputs`
    # objects.
    out = solve.(self.rotors, self.sections, self.ops)

    # Set the outputs.
    for out_str in keys(outputs)
        out_sym = Symbol(out_str)
        @. outputs[out_str] = getfield(out, out_sym)
    end

    return nothing
end

function set_inputs!(self::CCBladeResidualComp, inputs)
    # Rotor parameters.
    self.rotors.Rhub = inputs["Rhub"]
    self.rotors.Rtip = inputs["Rtip"]
    self.rotors.precone = inputs["precone"]

    # Blade section parameters.
    self.sections.r = inputs["r"]
    self.sections.chord = inputs["chord"]
    self.sections.theta = inputs["theta"]

    # Inflow parameters.
    self.ops.Vx = inputs["Vx"]
    self.ops.Vy = inputs["Vy"]
    self.ops.rho = inputs["rho"]
    self.ops.mu = inputs["mu"]
    self.ops.asound = inputs["asound"]

    return nothing
end

detect_apply_nonlinear(::Type{CCBladeResidualComp}) = false
detect_guess_nonlinear(::Type{CCBladeResidualComp}) = false
detect_apply_linear(::Type{CCBladeResidualComp}) = false


text:
"""
    Distributed schedule of when to increase evaluation
    resources and then compensate for the increase later.
"""
function evaluation_schedule(P, k, k_max, m, m_elite; p_UNUSED=0.2)
    N_max = k_max*m
    m_sched = floor(pdf(P, k)*N_max)
    if k == k_max # account for left-overs and over counting
        subtotal = sum(floor(pdf(P, i)*N_max) for i in 1:(k_max-1))
        m_sched = min(N_max - subtotal, N_max - m_sched)
    end
    mₑ = Int(m_sched)
    m_elite = min(m_elite, mₑ) # Clamp.
    return (mₑ, m_elite)
end

"""
    Schedule of when to increase evaluation resources and then
    compensate for the increase later.
"""
function evaluation_schedule_manual(k, m)
    if k == 1
        m+=10
    elseif k == 2
        m+=5
    elseif k == 3
        m-=1
    elseif 4 ≤ k ≤ 10
        m-=2
    end
    return m
end


text:
## MIT license (C) 2019 by Andrew Lyasoff

# Julia 1.1.0 code to illustrate the construction of univariate histograms
# and the Monte Carlo simulation technique for multivariate Gaussian samples with a given covariance matrix.
#See (3.67), p.91, in "Stochastic Methods in Asset Pricing."

begin
    using LinearAlgebra
    using SpecialFunctions
    using StatsBase
    using Random
    using Plots
    pyplot()
end

#We write our own histogram function. It takes as an input a single 1-dimensional array of data.
#The number of bins in the histogram is determined automatically by using the Diaconis-Friedman rule.
#The function returns two arrays: the mid-points of the bins and the (unnormalized) heights of the bars.

function hstgram(data_sample::Array{Float64,1})
    data_sorted=sort(data_sample)
    first=data_sorted[1]
    last=data_sorted[end]
    nmb=length(data_sorted)
    IQR=percentile(data_sorted,75)-percentile(data_sorted,25)
    bin_size_loc = 2*IQR*(nmb^(-1.0/3))
    num_bins=Int(floor((last-first)/bin_size_loc))
    bin_size=(last-first)/(num_bins)
    bin_end_points=[first+(i-1)*bin_size for i=1:(num_bins+1)]
    ahist_val=[length(data_sorted[data_sorted .< u]) for u in bin_end_points]
    hist_val=[ahist_val[i+1]-ahist_val[i] for i=1:num_bins]
    mid_bins=[first-bin_size/2+i*bin_size for i=1:num_bins]
    return mid_bins, hist_val
end

## First, create data sampled from the standard univariate normal density.
val=(x->((2*π)^(-1/2)*exp(-x^2/2))).(-3.3:0.05:3.3);

Random.seed!(0xabcdef12); # if needed to generate the same samples

#simulate 10,000 standard normals and generate the histogram
begin
    nval=randn!(zeros(10000));
    U,V=hstgram(nval);
    VV=V/(sum(V)*(U[2]-U[1]));
end

#check the length of the bin
length(VV)

begin
    plot(U.+(U[2]-U[1])/2,VV,line=(:steppre,1),linewidth=0.05,label="histogram")
    xlabel!("samples")
    ylabel!("frequency")
    plot!(-3.3:0.05:3.3,val,label="normal density")
end

#generate another sample
begin
    nval=randn(10000);
    U,V=hstgram(nval);
    VV=V/(sum(V)*(U[2]-U[1]));
    plot(U.+(U[2]-U[1])/2,VV,line=(:steppre,1),linewidth=0.05,label="histogram")
    xlabel!("samples")
    ylabel!("frequency")
    plot!(-3.3:0.05:3.3,val,label="normal density")
end

# standard normals can be generated by sampling from the uniform distribution
begin
    uval=rand(10000);
    nval=(x->sqrt(2)*erfinv(2*x-1)).(uval);
    U,V=hstgram(nval);
    VV=V/(sum(V)*(U[2]-U[1]));
    plot(U.+(U[2]-U[1])/2,VV,line=(:steppre,1),linewidth=0.05,label="histogram")
    xlabel!("samples")
    ylabel!("frequency")
    plot!(-3.3:0.05:3.3,val,label="normal density")
end

#simulate 10,000 standard bi-variate normals
begin
    nval=randn!(zeros(10000));
    nnval=randn!(zeros(10000));
end

scatter(nval,nnval,ratio=1,markersize=1,label="")

# to simulate bi-variate normals with a given covariance matrix
#   first create a fictitios covariance matrix

begin
    A=rand(2,2);
    #Cov=A'A; # another alternative that always yields a positive definite matrix
    Cov=Symmetric(A) # may not produce a positive definite matrix
end

#check if Cov is positive definite; if not, repeat the last step
eigvals(Cov)

begin
    eigdCov=Diagonal(eigvals(Cov).^0.5)
    eigCov=eigvecs(Cov); # matrix of eigen vectors
    chlCov=cholesky(Cov); # Cholesky "square root" of Cov
    MM=eigdCov*eigCov; # Spectral "square root" of Cov
end

# check that the factorizations give what is expected
Cov-MM'*MM
Cov-chlCov.L*chlCov.U
Cov-(chlCov.U)'*chlCov.U
Cov-(chlCov.L)*(chlCov.L)'

#this should be an orthogonal matrix
eigCov'*eigCov

#Transform the randomly generated standard bi-variate sample through the "square root" of the covariance matrix.
#method 1
begin
    NN=hcat(nval,nnval)';
    data_2_dim=MM'NN;
    scatter(data_2_dim[1,:],data_2_dim[2,:],ratio=1,markersize=1,label="")
end

#method 2
begin
    data_2_dim=(chlCov.L)*NN;
    scatter(data_2_dim[1,:],data_2_dim[2,:],ratio=1,markersize=1,label="")
end

max_stars_repo_name: AndrewLyasoff/SMAP

max_issues_repo_name: lhyzh/SMAP

max_forks_repo_name: lhyzh/SMAP


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

where `toml` is a Comby configuration.
"""
struct Pattern
    id::Int
    title::String
    tags::Vector{Any}
    description::String
    toml::String
end

unary_operators = ['!', '~', '+', '-']

val_rx = raw"[\w_\[\]:\.\\'*+-^′]+"
extend_val_rx_left = raw"(?<![+-\/*] )"
extend_val_rx_right = raw"(?! [+-\/*])"
fn_rx = raw"[\w_\.]+"

patterns = [
    Pattern(1, "Use rand() instead of rand(1)[1]", ["julia"],
        """
        According to [DNF](https://stackoverflow.com/questions/65403410/), many people use
        ```
        julia> using Random

        julia> Random.seed!(1234); rand(1)[1]
        0.5908446386657102
        ```
        instead of (the more performant)
        ```
        julia> Random.seed!(1234); rand()
        0.5908446386657102
        ```
        """,
        """
        [p1a]
        match="rand(1)[1]"
        rewrite="rand()"

        [p1b]
        match="rand(:[a],:[space~[ ]*]1)[1]"
        rewrite="rand(:[a])"
        """
    ),
    Pattern(2, "Use ismissing", ["julia", "performance-decrease"],
        """
        This pattern is disabled for now by default, due to discussion at
        <https://discourse.julialang.org/t/ismissing-x-versus-x-missing/52171>.
        """,
        raw"""
        [p2]
        match=':[[var]] :[first~[=|!]]== missing'

        rule='where rewrite :[first] { "=" -> "" }'

        rewrite=':[first]ismissing(:[var])'
        """
    ),
    Pattern(3, "Avoid x -> f(x)", ["julia"],
        """
        From the [Julia Style Guide](https://docs.julialang.org/en/v1/manual/style-guide/).
        For unary operators, such as `x -> !x`, see [pattern 8](#8).

        **Example**
        ```
        julia> f(x) = 2x
        f (generic function with 1 method)

        julia> map(x -> f(x), [1, 2])
        2-element Array{Int64,1}:
         2
         4

        julia> map(f, [1, 2])
        2-element Array{Int64,1}:
         2
         4
        ```
        """,
        raw"""
        [p3]
        match=':[[x]] -> :[f~[\w_]*](:[x]):[end~(,|\n)]'

        rewrite=':[f]:[end]'
        """
    ),
    Pattern(4, "Omit comparison with boolean constant", ["generic"],
        "From [staticcheck](https://staticcheck.io/docs/checks).",
        """
        [p4]
        match=':[[var]] :[first~(=|!)]= :[bool~(true|false)]'

        rule='where
            rewrite :[first] { "=" -> "" },
            rewrite :[bool] { "false" -> "!" },
            rewrite :[bool] { "true" -> "" }'

        rewrite=':[first]:[bool]:[var]'
        """
    ),
    Pattern(5, "Omit a == a and a != a", ["julia"],
        raw"""
        SA4000 in [staticcheck](https://staticcheck.io/docs/checks).

        Be careful when applying this pattern to situations where `a` can be of type `Missing`.
        For example, `a == a` returns `missing` and not `true`.
        On the other hand, `if missing` returns an error, so maybe this isn't a valid argument.

        **Examples**
        ```
        julia> a == 1
        true

        julia> a == a
        true

        julia> a != a
        false
        ```
        """,
        # Keep this pattern specific to avoid massive slowdown.
        # So, don't do ":[x] :[bool] :[x]" without "if".
        """
        [p5a]
        match='if :[x] :[bool~(=|!)]= :[x]'
        rule='where
            rewrite :[bool] { "=" -> "true" },
            rewrite :[bool] { "!" -> "false" }'
        rewrite='if :[bool]'

        [p5b]
        match=''':[x~$val_rx]:[bool~(=|!)]=:[x~$val_rx]'''
        rule='where
            rewrite :[bool] { "=" -> "true" },
            rewrite :[bool] { "!" -> "false" }'
        rewrite=':[bool]'
        """
    ),
    Pattern(6, "Avoid comparing findfirst to nothing", ["julia"],
        """
        For example, instead of `findfirst('a', \"ab\") === nothing`, use `occursin('a', \"ab\")` or `contains` (Julia ≥1.5).
        """,
        """
        [p6]
        match='findfirst(:[a], :[b]) :[bool~(=|!)]== nothing'

        rule='where
            rewrite :[bool] { "!" -> "" },
            rewrite :[bool] { "=" -> "!" }'

        rewrite=':[bool]occursin(:[a], :[b])'
        """
    ),
    Pattern(7, "Avoid findall(x -> x == false, Y)", ["julia"],
        """
        Instead, use `findall(.!Y)`.
        """,
        """
        [p7]
        match='findall(:[x] -> :[x] == false, :[Y])'

        rewrite='findall(.!:[Y])'
        """
    ),
    Pattern(8, "Use unary operator without anonymous function", ["julia"],
        """
        This is an extension of [pattern 3](#3) for the unary operators:
        `$(join(unary_operators, "`, `"))`.

        **Example**
        ```
        julia> map(x -> !x, [true, missing])
        2-element Array{Union{Missing, Bool},1}:
         false
              missing

        julia> map(!, [true, missing]
        2-element Array{Union{Missing, Bool},1}:
         false
              missing
        ```
        """,
        """
        [p8]
        match=':[[x]] -> :[unary~[$(join(unary_operators))]]:[x~[\\w\\_\\[\\]\\.]+]'

        rewrite=':[unary]'
        """
    ),
    Pattern(9, "Abbreviate keyword argument values", ["julia", 1.5],
        """
        From [Julia 1.5 onwards](https://julialang.org/blog/2020/08/julia-1.5-highlights)
        it is possible to abbreviate
        ```
        printstyled("text"; color = color)
        ```
        to
        ```
        printstyled("text"; color)
        ```
        """,
        """
        [p9]
        match="(:[a]; :[c] = :[c])"
        rewrite="(:[a]; :[c])"
        """
    ),
    Pattern(10, "Replace P && P with P", ["julia", "disabled"],
        """
        *Disabled because it has many false-positives and no true-positives.*

        The statement `P && P` is logically equivalent to `P`.
        This is called a micro-clone \\citep{pTonder2016}.

        **Example**
        ```
        julia> P = 3 < 1
        false

        julia> P && P
        false

        julia> P
        false
        ```
        """,
        """
        [p10a]
        match="if :[a] == :[b] && :[a] == :[b]"
        rewrite="if :[a] == :[b]"
        [p10b]
        match=":[a~[^ ]+] && :[a~[^ ]+]"
        rewrite=":[a]"
        """
    ),
    Pattern(11, "Replace dirname(dirname(pathof(M))) with pkgdir(M)", ["julia", 1.4],
        """
        From Julia 1.4 onwards it is possible to use `pkgdir(M)`.

        **Example**
        ```
        julia> using DataFrames

        julia> pathof(DataFrames)
        "/home/rik/.julia/packages/DataFrames/oQ5c7/src/DataFrames.jl"

        julia> dirname(dirname(pathof(DataFrames)))
        "/home/rik/.julia/packages/DataFrames/oQ5c7"

        julia> pkgdir(DataFrames)
        "/home/rik/.julia/packages/DataFrames/oQ5c7"
        ```
        """,
        """
        [p11]
        match="dirname(dirname(pathof(:[M])))"
        rewrite="pkgdir(:[M])"
        """
    )
]

function patterns_have_valid_indexes()::Bool
    collect(1:length(patterns)) == [p.id for p in patterns]
end

max_stars_repo_name: rikhuijzer/display-lint

max_issues_repo_name: rikhuijzer/display-lint

max_forks_repo_name: rikhuijzer/ShowLint


text:
using RvSpectraKitLearn
using MultivariateStats

# Likely need to update this to point to directory containing the spectra you wish to analyze
include(joinpath(Pkg.dir("RvSpectraKitLearn"),"path_to_spectra.jl"))
cd(path_to_spectra)

(datafiles, phases) = get_filenames(); # num_max_files=4);
(lambda, obs_noisefree) = read_filelist(datafiles);

sampling = 3.0
planet_amplitude = 10.0
rvs_true = -planet_amplitude*cos.(2pi*phases)

snr = 100.0
obs = RvSpectraKitLearn.make_noisy_spectra(obs_noisefree,snr,sampling=sampling)

pca_out = fit_pca_default(obs);
pca_scores = transform(pca_out,obs)   # compute scores using Julia's PCA routine for comparison's sake
pca_eford_out = fit_pca_eford(obs);   # compute first few scores using Eric's itterative PCA routine

doppler_comp_simple = calc_doppler_component_simple(lambda,obs);
genpca_simple_out = fit_gen_pca_rv(obs,doppler_comp_simple)
rvs_out_simple = est_rvs_from_pc(obs,genpca_simple_out[1],genpca_simple_out[2][:,1])

doppler_comp_gp = calc_doppler_component_gp(lambda,obs);
genpca_gp_out = fit_gen_pca_rv(obs,doppler_comp_gp)
rvs_out_gp = est_rvs_from_pc(obs,genpca_gp_out[1],genpca_gp_out[2][:,1])
rms_simple = sqrt(sumabs2(rvs_out_simple.-rvs_true)/length(rvs_true))
rms_gp = sqrt(sumabs2(rvs_out_gp.-rvs_true)/length(rvs_true))
println("snr = ", snr, " rms = ", rms_simple, "  ", rms_gp)

using PyPlot
plot(phases,rvs_true,"b-")
 plot(phases,rvs_out_simple,"r.")
 plot(phases,rvs_out_gp,"g.")
 plot(phases,rvs_out_simple.-rvs_true,"r.")
 plot(phases,rvs_out_gp.-rvs_true,"g.")
 xlabel("Phase")
 ylabel("Velocity (m/s)")

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
#EM500N3LO

## LO con  10^4 GeV^-2
const C0_1S0  = -0.147167
const C0_3S1  = -0.11897839697675265
const C_CSB = 0.00049950
const C_CIB = 0.00069075
## NLO con 10^4 GeV^-4
const C2_3S1  = 0.76
const C2_3P0  = 1.487
const C2_1P1  = 0.656
const C2_3P1  = -0.630
const C2_1S0  = 2.380
const C2_3SD1 = 0.826
const C2_3P2  = -0.538
## N3LO con 10^4 GeV^-6
const hD_1S0 = -2.545
const D_1S0 = -16.0
const D_1P1 = 5.25
const D_3P0 = 0.245
const D_3P1 = 2.35
const D_3P2 = 2.295
const hD_3S1 = 7.0
const D_3S1 = 6.55
const hD_3SD1 = 2.25
const D_3SD1 = 6.61
const D_3D1 = -2.80
const D_1D2 = -1.77
const D_3D2 = -1.46
const D_3PF2 = -0.465
const D_3D3 = 5.66

## For pion exchange contribution
## NNLO GeV^-1
const c1_NNLO =  -0.81
const c2_NNLO =  2.80
const c3_NNLO = -3.20
const c4_NNLO =  5.40
## N3LO GeV^-2
const d12 = 3.06
const d3 = -3.27
const d5 = 0.45
const d145 = -5.65
### NLOvs
const c_vs_1 =  0.19162
const c_vs_2 = -0.28374
const c_vs_3 =  0.02680
const c_vs_4 = -0.34499
const c_vs_5 = -0.16235

## for eff3nf
const cD =  1.621
const cE = -0.143


text:
# This file is part of Kpax3. License is MIT.

abstract type PriorColPartition end

mutable struct AminoAcidPriorCol <: PriorColPartition
  logγ::Vector{Float64}
  logω::Vector{Vector{Float64}}

  A::Matrix{Float64}
  B::Matrix{Float64}
end

function AminoAcidPriorCol(data::Matrix{UInt8},
                           γ::Vector{Float64},
                           r::Float64)
  if length(γ) != 3
    throw(KInputError("Argument 'γ' does not have length 3."))
  elseif γ[1] < 0
    throw(KDomainError("Argument 'γ[1]' is negative."))
  elseif γ[2] < 0
    throw(KDomainError("Argument 'γ[2]' is negative."))
  elseif γ[3] < 0
    throw(KDomainError("Argument 'γ[3]' is negative."))
  end

  if r <= 0
    throw(KDomainError("Argument 'r' is negative."))
  end

  (m, n) = size(data)

  # probabilities must sum to one
  tot = γ[1] + γ[2] + γ[3]
  γ[1] /= tot
  γ[2] /= tot
  γ[3] /= tot

  # log(0) is -Inf
  # since we are going to multiply this matrix with a positive scalar (1 / k)
  # no NaN can be produced even if some γ[i] are zero
  # attributes 3 and 4 are both from property 3 -> use γ[3] twice
  logγ = [log(γ[1]); log(γ[2]); log(γ[3])]

  logω = Vector{Float64}[[log(k - 1) - log(k); -log(k)] for k in 1:n]

  n1s = zeros(Float64, m)
  for a in 1:n, b in 1:m
    n1s[b] += float(data[b, a])
  end

  A = zeros(Float64, 4, m)
  B = zeros(Float64, 4, m)

  if n > r
    for b in 1:m
      # uninformative attributes
      # If n > r, these two parameters sum to (r+1), i.e. the theoretical
      # sample size for the characteristic hyperparameters. This is done so
      # because we don't want them to be overwhelmed by the data.
      # The mean is the same to the one obtained with a Jeffreys prior
      A[1, b] = (r + 1) * (n1s[b] + 0.5) / (n + 1)
      B[1, b] = (r + 1) - A[1, b]

      # informative but not characteristic for any cluster
      A[2, b] = 1
      B[2, b] = 1

      # informative and characteristic... but for another cluster
      A[3, b] = 1
      B[3, b] = r

      # informative and characteristic for this cluster
      A[4, b] = r
      B[4, b] = 1
    end
  else
    for b in 1:m
      # uninformative attributes
      A[1, b] = n1s[b] + 0.5
      B[1, b] = n - n1s[b] + 0.5

      # informative but not characteristic for any cluster
      A[2, b] = 1
      B[2, b] = 1

      # informative and characteristic... but for another cluster
      A[3, b] = 1
      B[3, b] = r

      # informative and characteristic for this cluster
      A[4, b] = r
      B[4, b] = 1
    end
  end

  AminoAcidPriorCol(logγ, logω, A, B)
end


text:
struct T1
    x::NTuple{1000, Int}
    y::Int
end

mutable struct T2
    x::NTuple{1000, Int}
    y::Int
end

function worker1()
    p = T1(ntuple(x->1, 1000), 0)
    for i in 1:10^6
        p = T1(p.x, p.y+1)
    end
    p
end

function worker2()
    p = T2(ntuple(x->1, 1000), 0)
    for i in 1:10^6
        p.y += 1
    end
    p
end

max_stars_repo_path: Chapter06/07. The efficiency of mutable versus immutable types/work.jl


text:
module Harmonie_namelists

import JSONSchema, YAML
using OrderedCollections, TOML
import Base.parse

# import Base.isvalid

# export schema, isvalid, diagnose
export NAMELIST_DIR
export dict2namelist, read_namelist, merge_namelists, replace_env!

const moduledir=@__DIR__
# const schemafile="$moduledir/namelists.schema.json"
#const schema = JSONSchema.Schema(read(schemafile,String),parentFileDirectory="$moduledir/")
# const schema = JSONSchema.Schema(read(schemafile,String))
# const NAMELIST_DIR = "$moduledir/namelists/ifs"
const NAMELIST_DIR = "$moduledir/namelists/"
const TYPE_INFO = YAML.load(open("$NAMELIST_DIR/type_info.yaml"))
const DEFAULTS  = YAML.load(open("$NAMELIST_DIR/defaults.yaml"))
const MYENV =  merge((env, def) -> env, ENV,DEFAULTS)

include("dict2namelist.jl")
include("replace_env.jl")
# include("jsonschema.jl")

# Why is this not in Base?
Base.parse(::Type{String},str::String) = str
Base.parse(::Type{Int64}, int::Int64)  = int
Base.parse(::Type{Bool}, bool::Bool)  = bool

read_namelist(name) = YAML.load(open("$NAMELIST_DIR/yaml/ifs/$name.yaml"))
# read_namelist(name) = TOML.parsefile("$NAMELIST_DIR/toml/ifs/$name.toml")

"""
    merge_namelists(dicts)

Returns a dictionary based on an arrays of dicts
"""
merge_namelists(dicts) = merge!(merge,OrderedDict{String,Any}(),dicts...)

end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

abstract type MCMCConvergenceTestResult end

function check_convergence!(
    ct::MCMCConvergenceTest,
    chains::AbstractVector{<:MCMCIterator},
    samples::AbstractVector{<:DensitySampleVector}
)
    result = check_convergence(ct, samples)
    for chain in chains
        chain.info = MCMCIteratorInfo(chain.info, converged = result.converged)
    end
    result
end

"""
    gr_Rsqr(stats::AbstractVector{<:MCMCBasicStats})
    gr_Rsqr(samples::AbstractVector{<:DensitySampleVector})

*BAT-internal, not part of stable public API.*

Gelman-Rubin ``\$R^2\$`` for all DOF.
"""
function gr_Rsqr end

function gr_Rsqr(stats::AbstractVector{<:MCMCBasicStats})
    m = totalndof(first(stats))
    W = mean([cs.param_stats.cov[i,i] for cs in stats, i in 1:m], dims=1)[:]
    B = var([cs.param_stats.mean[i] for cs in stats, i in 1:m], dims=1)[:]
    (W .+ B) ./ W
end

function gr_Rsqr(samples::AbstractVector{<:DensitySampleVector})
    gr_Rsqr(MCMCBasicStats.(samples))
end

"""
    struct GelmanRubinConvergence <: MCMCConvergenceTest

Gelman-Rubin maximum R^2 convergence test.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct GelmanRubinConvergence <: MCMCConvergenceTest
    threshold::Float64 = 1.1
end

export GelmanRubinConvergence

struct GRConvergenceResult <: MCMCConvergenceTestResult
    converged::Bool
    max_Rsqr::Float64
end

function check_convergence(ct::GelmanRubinConvergence, samples::AbstractVector{<:DensitySampleVector})
    max_Rsqr = maximum(gr_Rsqr(samples))
    converged = max_Rsqr <= ct.threshold
    @debug begin
        success_str = converged ? "have" : "have *not*"
        "Chains $success_str converged, max(R^2) = $(max_Rsqr), threshold = $(ct.threshold)"
    end
    GRConvergenceResult(converged, max_Rsqr)
end

@doc doc"""
    bg_R_2sqr(stats::AbstractVector{<:MCMCBasicStats}; corrected::Bool = false)
    bg_R_2sqr(samples::AbstractVector{<:DensitySampleVector}; corrected::Bool = false)

*BAT-internal, not part of stable public API.*

Brooks-Gelman R_2^2 for all DOF.
If normality is assumed, 'corrected' should be set to true to account for the sampling variability.
"""
function bg_R_2sqr(stats::AbstractVector{<:MCMCBasicStats}; corrected::Bool = false)
    p = totalndof(first(stats))
    m = length(stats)
    n = mean(Float64.(nsamples.(stats)))

    σ_W = var([cs.param_stats.cov[i,i] for cs in stats, i in 1:p], dims = 1)[:]
    B  = var([cs.param_stats.mean[i] for cs in stats, i in 1:p], dims = 1)[:]
    W = mean([cs.param_stats.cov[i,i] for cs in stats, i in 1:p], dims = 1)[:]

    σ_sq = m * (n - 1) / (m*n - 1) * W + n * (m - 1) / (m*n - 1) * B

    R_unc = σ_sq ./ W

    if corrected == false
        return R_unc
    end

    σ_ij = [cs.param_stats.cov[i,i] for cs in stats, i in 1:p]
    x_ij = [cs.param_stats.mean[i] for cs in stats, i in 1:p]

    cov_σx = [cov(σ_ij[:,j], x_ij[:,j]) for j in 1:p]
    cov_σx_sq = [cov(σ_ij[:,j], x_ij[:,j].^2) for j in 1:p]

    N = (n-1)/n
    M = (m-1)/m
    V = N*σ_sq + M*B

    σ_V = N^2/m*σ_W + 2*M/(m-1)*B.^2 + 2*M*N/m*(cov_σx_sq - 2*B.*cov_σx)
    d = 2 * V.^2 ./ σ_V

    R_unc.*(d.+3)./(d.+1)
end

function bg_R_2sqr(samples::AbstractVector{<:DensitySampleVector}; corrected::Bool = false)
    bg_R_2sqr(MCMCBasicStats.(samples), corrected = corrected)
end

"""
    struct BrooksGelmanConvergence <: MCMCConvergenceTest

Brooks-Gelman maximum R^2 convergence test.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct BrooksGelmanConvergence <: MCMCConvergenceTest
    threshold::Float64 = 1.1
    corrected::Bool = false
end

export BrooksGelmanConvergence

struct BGConvergenceResult <: MCMCConvergenceTestResult
    converged::Bool
    max_Rsqr::Float64
end

function check_convergence(ct::BrooksGelmanConvergence, samples::AbstractVector{<:DensitySampleVector})
    max_Rsqr = maximum(bg_R_2sqr(samples, corrected = ct.corrected))
    converged = max_Rsqr <= ct.threshold
    @debug begin
        success_str = converged ? "have" : "have *not*"
        "Chains $success_str converged, max(R^2) = $(max_Rsqr), threshold = $(ct.threshold)"
    end
    GRConvergenceResult(converged, max_Rsqr)
end

function check_convergence(ct::BrooksGelmanConvergence, samples::DensitySampleVector)

    # create a vector of chains
    chains_ind = unique([i.chainid for i in samples.info])
    vector_chains = DensitySampleVector[]
    for i in chains_ind
        mask_chain = [j.chainid == i for j in samples.info]
        push!(vector_chains, samples[mask_chain])
    end

    check_convergence(ct, vector_chains)

end


text:
# Write the output of psim_vrf to files
function writedlm_vrf(results::Dict{ASCIIString, Any}, fpath::ASCIIString, dlm)
  for i in keys(results)
    if isa(results[i], AbstractArray)
      writedlm(fpath*"_"*i*".txt", results[i], dlm)
    end
  end
end

# Read file and rewrite it into a new one with specified format
function writedlm_reformat(infile::ASCIIString, outfile::ASCIIString, outfmt::Array{ASCIIString, 1};
  indlm::ASCIIString=" ", outdlm::ASCIIString=" ", outendline::ASCIIString="\n")
  instream, outstream = open(infile), open(outfile, "w")

  while !eof(instream)
    inline, outline = readline(instream), string()
    if ismatch(r"[^\s\t\n\r]+", inline) && !ismatch(r"^#", inline)
      tokens = split(rstrip(inline), indlm)
      for i in 1:length(tokens)
        outline *= (@eval @sprintf($(outfmt[i]), parse(utf8($(tokens[i])))))*(i==length(tokens) ? outendline : outdlm)
      end
      write(outstream, outline)
    end
  end

  close(instream)
  close(outstream)
end


text:
using Mocking

add_attributes_to_findings_1_patch = @patch function add_attributes_to_findings(a...; b...)
    return Dict{String, Any}("failedItems" => Dict{String, Any}())
end

create_assessment_target_1_patch = @patch function create_assessment_target(a...; b...)
    return Dict{String, Any}("assessmentTargetArn" => "arn:aws:inspector:us-west-2:123456789012:target/0-nvgVhaxX")
end

create_assessment_template_1_patch = @patch function create_assessment_template(a...; b...)
    return Dict{String, Any}("assessmentTemplateArn" => "arn:aws:inspector:us-west-2:123456789012:target/0-nvgVhaxX/template/0-it5r2S4T")
end

create_resource_group_1_patch = @patch function create_resource_group(a...; b...)
    return Dict{String, Any}("resourceGroupArn" => "arn:aws:inspector:us-west-2:123456789012:resourcegroup/0-AB6DMKnv")
end

describe_assessment_runs_1_patch = @patch function describe_assessment_runs(a...; b...)
    return Dict{String, Any}("failedItems" => Dict{String, Any}(), "assessmentRuns" => Any[Dict{String, Any}("findingCounts" => Dict{String, Any}("Low" => 0, "Medium" => 2, "Informational" => 0, "High" => 14, "Undefined" => 0), "assessmentTemplateArn" => "arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq/template/0-4r1V2mAw", "createdAt" => "1458680170.035", "state" => "COMPLETED", "name" => "Run 1 for ExampleAssessmentTemplate", "dataCollected" => true, "arn" => "arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq/template/0-4r1V2mAw/run/0-MKkpXXPE", "completedAt" => "1458680301.4", "durationInSeconds" => 3600, "stateChanges" => Any[Dict{String, Any}("stateChangedAt" => "1458680170.035", "state" => "CREATED"), Dict{String, Any}("stateChangedAt" => "1458680170.065", "state" => "START_DATA_COLLECTION_PENDING"), Dict{String, Any}("stateChangedAt" => "1458680170.096", "state" => "START_DATA_COLLECTION_IN_PROGRESS"), Dict{String, Any}("stateChangedAt" => "1458680170.161", "state" => "COLLECTING_DATA"), Dict{String, Any}("stateChangedAt" => "1458680239.883", "state" => "STOP_DATA_COLLECTION_PENDING"), Dict{String, Any}("stateChangedAt" => "1458680299.847", "state" => "DATA_COLLECTED"), Dict{String, Any}("stateChangedAt" => "1458680300.099", "state" => "EVALUATING_RULES"), Dict{String, Any}("stateChangedAt" => "1458680301.4", "state" => "COMPLETED")], "notifications" => Any[], "rulesPackageArns" => Any["arn:aws:inspector:us-west-2:758058086616:rulespackage/0-X1KXtawP"], "startedAt" => "1458680170.161", "userAttributesForFindings" => Any[], "stateChangedAt" => "1458680301.4")])
end

describe_assessment_targets_1_patch = @patch function describe_assessment_targets(a...; b...)
    return Dict{String, Any}("failedItems" => Dict{String, Any}(), "assessmentTargets" => Any[Dict{String, Any}("name" => "ExampleAssessmentTarget", "updatedAt" => "1458074191.459", "createdAt" => "1458074191.459", "arn" => "arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq", "resourceGroupArn" => "arn:aws:inspector:us-west-2:123456789012:resourcegroup/0-PyGXopAI")])
end

describe_assessment_templates_1_patch = @patch function describe_assessment_templates(a...; b...)
    return Dict{String, Any}("assessmentTemplates" => Any[Dict{String, Any}("rulesPackageArns" => Any["arn:aws:inspector:us-west-2:758058086616:rulespackage/0-X1KXtawP"], "name" => "ExampleAssessmentTemplate", "userAttributesForFindings" => Any[], "durationInSeconds" => 3600, "assessmentRunCount" => 0, "createdAt" => "1458074191.844", "arn" => "arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq/template/0-4r1V2mAw", "assessmentTargetArn" => "arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq")], "failedItems" => Dict{String, Any}())
end

describe_cross_account_access_role_1_patch = @patch function describe_cross_account_access_role(a...; b...)
    return Dict{String, Any}("roleArn" => "arn:aws:iam::123456789012:role/inspector", "valid" => true, "registeredAt" => "1458069182.826")
end

describe_findings_1_patch = @patch function describe_findings(a...; b...)
    return Dict{String, Any}("failedItems" => Dict{String, Any}(), "findings" => Any[Dict{String, Any}("indicatorOfCompromise" => false, "recommendation" => "No remediation needed.", "schemaVersion" => 1, "createdAt" => "1458680301.37", "confidence" => 10, "service" => "Inspector", "description" => "Amazon Inspector did not find any potential security issues during this assessment.", "assetAttributes" => Dict{String, Any}("ipv4Addresses" => Any[], "schemaVersion" => 1), "arn" => "arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq/template/0-4r1V2mAw/run/0-MKkpXXPE/finding/0-HwPnsDm4", "serviceAttributes" => Dict{String, Any}("assessmentRunArn" => "arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq/template/0-4r1V2mAw/run/0-MKkpXXPE", "rulesPackageArn" => "arn:aws:inspector:us-west-2:758058086616:rulespackage/0-X1KXtawP", "schemaVersion" => 1), "severity" => "Informational", "updatedAt" => "1458680301.37", "userAttributes" => Any[], "numericSeverity" => 0, "assetType" => "ec2-instance", "attributes" => Any[], "title" => "No potential security issues found")])
end

describe_resource_groups_1_patch = @patch function describe_resource_groups(a...; b...)
    return Dict{String, Any}("failedItems" => Dict{String, Any}(), "resourceGroups" => Any[Dict{String, Any}("tags" => Any[Dict{String, Any}("key" => "Name", "value" => "example")], "createdAt" => "1458074191.098", "arn" => "arn:aws:inspector:us-west-2:123456789012:resourcegroup/0-PyGXopAI")])
end

describe_rules_packages_1_patch = @patch function describe_rules_packages(a...; b...)
    return Dict{String, Any}("rulesPackages" => Any[Dict{String, Any}("name" => "Security Best Practices", "provider" => "Amazon Web Services, Inc.", "version" => "1.1", "arn" => "arn:aws:inspector:us-west-2:758058086616:rulespackage/0-JJOtZiqQ", "description" => "The rules in this package help determine whether your systems are configured securely.")], "failedItems" => Dict{String, Any}())
end

get_telemetry_metadata_1_patch = @patch function get_telemetry_metadata(a...; b...)
    return Dict{String, Any}("telemetryMetadata" => Any[Dict{String, Any}("dataSize" => 345, "count" => 2, "messageType" => "InspectorDuplicateProcess"), Dict{String, Any}("dataSize" => 255, "count" => 3, "messageType" => "InspectorTimeEventMsg"), Dict{String, Any}("dataSize" => 1082, "count" => 4, "messageType" => "InspectorNetworkInterface"), Dict{String, Any}("dataSize" => 349, "count" => 2, "messageType" => "InspectorDnsEntry"), Dict{String, Any}("dataSize" => 2514, "count" => 11, "messageType" => "InspectorDirectoryInfoMsg"), Dict{String, Any}("dataSize" => 179, "count" => 1, "messageType" => "InspectorTcpV6ListeningPort"), Dict{String, Any}("dataSize" => 10949, "count" => 101, "messageType" => "InspectorTerminal"), Dict{String, Any}("dataSize" => 5916, "count" => 26, "messageType" => "InspectorUser"), Dict{String, Any}("dataSize" => 32148, "count" => 282, "messageType" => "InspectorDynamicallyLoadedCodeModule"), Dict{String, Any}("dataSize" => 10172, "count" => 18, "messageType" => "InspectorCreateProcess"), Dict{String, Any}("dataSize" => 8001, "count" => 3, "messageType" => "InspectorProcessPerformance"), Dict{String, Any}("dataSize" => 360, "count" => 1, "messageType" => "InspectorOperatingSystem"), Dict{String, Any}("dataSize" => 546, "count" => 6, "messageType" => "InspectorStopProcess"), Dict{String, Any}("dataSize" => 1553, "count" => 1, "messageType" => "InspectorInstanceMetaData"), Dict{String, Any}("dataSize" => 434, "count" => 2, "messageType" => "InspectorTcpV4Connection"), Dict{String, Any}("dataSize" => 2960322, "count" => 474, "messageType" => "InspectorPackageInfo"), Dict{String, Any}("dataSize" => 2235, "count" => 3, "messageType" => "InspectorSystemPerformance"), Dict{String, Any}("dataSize" => 46048, "count" => 105, "messageType" => "InspectorCodeModule"), Dict{String, Any}("dataSize" => 182, "count" => 1, "messageType" => "InspectorUdpV6ListeningPort"), Dict{String, Any}("dataSize" => 371, "count" => 2, "messageType" => "InspectorUdpV4ListeningPort"), Dict{String, Any}("dataSize" => 8362, "count" => 18, "messageType" => "InspectorKernelModule"), Dict{String, Any}("dataSize" => 48788, "count" => 29, "messageType" => "InspectorConfigurationInfo"), Dict{String, Any}("dataSize" => 79, "count" => 1, "messageType" => "InspectorMonitoringStart"), Dict{String, Any}("dataSize" => 0, "count" => 5, "messageType" => "InspectorSplitMsgBegin"), Dict{String, Any}("dataSize" => 4593, "count" => 51, "messageType" => "InspectorGroup"), Dict{String, Any}("dataSize" => 184, "count" => 1, "messageType" => "InspectorTcpV4ListeningPort"), Dict{String, Any}("dataSize" => 3146579, "count" => 1159, "messageType" => "Total"), Dict{String, Any}("dataSize" => 0, "count" => 5, "messageType" => "InspectorSplitMsgEnd"), Dict{String, Any}("dataSize" => 612, "count" => 1, "messageType" => "InspectorLoadImageInProcess")])
end

list_assessment_run_agents_1_patch = @patch function list_assessment_run_agents(a...; b...)
    return Dict{String, Any}("assessmentRunAgents" => Any[Dict{String, Any}("assessmentRunArn" => "arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq/template/0-4r1V2mAw/run/0-MKkpXXPE", "agentHealth" => "HEALTHY", "telemetryMetadata" => Any[Dict{String, Any}("dataSize" => 345, "count" => 2, "messageType" => "InspectorDuplicateProcess"), Dict{String, Any}("dataSize" => 255, "count" => 3, "messageType" => "InspectorTimeEventMsg"), Dict{String, Any}("dataSize" => 1082, "count" => 4, "messageType" => "InspectorNetworkInterface"), Dict{String, Any}("dataSize" => 349, "count" => 2, "messageType" => "InspectorDnsEntry"), Dict{String, Any}("dataSize" => 2514, "count" => 11, "messageType" => "InspectorDirectoryInfoMsg"), Dict{String, Any}("dataSize" => 179, "count" => 1, "messageType" => "InspectorTcpV6ListeningPort"), Dict{String, Any}("dataSize" => 10949, "count" => 101, "messageType" => "InspectorTerminal"), Dict{String, Any}("dataSize" => 5916, "count" => 26, "messageType" => "InspectorUser"), Dict{String, Any}("dataSize" => 32148, "count" => 282, "messageType" => "InspectorDynamicallyLoadedCodeModule"), Dict{String, Any}("dataSize" => 10172, "count" => 18, "messageType" => "InspectorCreateProcess"), Dict{String, Any}("dataSize" => 8001, "count" => 3, "messageType" => "InspectorProcessPerformance"), Dict{String, Any}("dataSize" => 360, "count" => 1, "messageType" => "InspectorOperatingSystem"), Dict{String, Any}("dataSize" => 546, "count" => 6, "messageType" => "InspectorStopProcess"), Dict{String, Any}("dataSize" => 1553, "count" => 1, "messageType" => "InspectorInstanceMetaData"), Dict{String, Any}("dataSize" => 434, "count" => 2, "messageType" => "InspectorTcpV4Connection"), Dict{String, Any}("dataSize" => 2960322, "count" => 474, "messageType" => "InspectorPackageInfo"), Dict{String, Any}("dataSize" => 2235, "count" => 3, "messageType" => "InspectorSystemPerformance"), Dict{String, Any}("dataSize" => 46048, "count" => 105, "messageType" => "InspectorCodeModule"), Dict{String, Any}("dataSize" => 182, "count" => 1, "messageType" => "InspectorUdpV6ListeningPort"), Dict{String, Any}("dataSize" => 371, "count" => 2, "messageType" => "InspectorUdpV4ListeningPort"), Dict{String, Any}("dataSize" => 8362, "count" => 18, "messageType" => "InspectorKernelModule"), Dict{String, Any}("dataSize" => 48788, "count" => 29, "messageType" => "InspectorConfigurationInfo"), Dict{String, Any}("dataSize" => 79, "count" => 1, "messageType" => "InspectorMonitoringStart"), Dict{String, Any}("dataSize" => 0, "count" => 5, "messageType" => "InspectorSplitMsgBegin"), Dict{String, Any}("dataSize" => 4593, "count" => 51, "messageType" => "InspectorGroup"), Dict{String, Any}("dataSize" => 184, "count" => 1, "messageType" => "InspectorTcpV4ListeningPort"), Dict{String, Any}("dataSize" => 3146579, "count" => 1159, "messageType" => "Total"), Dict{String, Any}("dataSize" => 0, "count" => 5, "messageType" => "InspectorSplitMsgEnd"), Dict{String, Any}("dataSize" => 612, "count" => 1, "messageType" => "InspectorLoadImageInProcess")], "agentHealthCode" => "RUNNING", "agentId" => "i-49113b93")], "nextToken" => "1")
end

list_assessment_runs_1_patch = @patch function list_assessment_runs(a...; b...)
    return Dict{String, Any}("assessmentRunArns" => Any["arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq/template/0-4r1V2mAw/run/0-MKkpXXPE", "arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq/template/0-4r1V2mAw/run/0-v5D6fI3v"], "nextToken" => "1")
end

list_assessment_targets_1_patch = @patch function list_assessment_targets(a...; b...)
    return Dict{String, Any}("assessmentTargetArns" => Any["arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq"], "nextToken" => "1")
end

list_assessment_templates_1_patch = @patch function list_assessment_templates(a...; b...)
    return Dict{String, Any}("assessmentTemplateArns" => Any["arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq/template/0-4r1V2mAw", "arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq/template/0-Uza6ihLh"], "nextToken" => "1")
end

list_event_subscriptions_1_patch = @patch function list_event_subscriptions(a...; b...)
    return Dict{String, Any}("subscriptions" => Any[Dict{String, Any}("resourceArn" => "arn:aws:inspector:us-west-2:123456789012:target/0-nvgVhaxX/template/0-7sbz2Kz0", "eventSubscriptions" => Any[Dict{String, Any}("event" => "ASSESSMENT_RUN_COMPLETED", "subscribedAt" => "1459455440.867")], "topicArn" => "arn:aws:sns:us-west-2:123456789012:exampletopic")], "nextToken" => "1")
end

list_findings_1_patch = @patch function list_findings(a...; b...)
    return Dict{String, Any}("findingArns" => Any["arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq/template/0-4r1V2mAw/run/0-MKkpXXPE/finding/0-HwPnsDm4", "arn:aws:inspector:us-west-2:123456789012:target/0-0kFIPusq/template/0-4r1V2mAw/run/0-v5D6fI3v/finding/0-tyvmqBLy"], "nextToken" => "1")
end

list_rules_packages_1_patch = @patch function list_rules_packages(a...; b...)
    return Dict{String, Any}("rulesPackageArns" => Any["arn:aws:inspector:us-west-2:758058086616:rulespackage/0-9hgA516p", "arn:aws:inspector:us-west-2:758058086616:rulespackage/0-H5hpSawc", "arn:aws:inspector:us-west-2:758058086616:rulespackage/0-JJOtZiqQ", "arn:aws:inspector:us-west-2:758058086616:rulespackage/0-vg5GGHSD"], "nextToken" => "1")
end

list_tags_for_resource_1_patch = @patch function list_tags_for_resource(a...; b...)
    return Dict{String, Any}("tags" => Any[Dict{String, Any}("key" => "Name", "value" => "Example")])
end

preview_agents_1_patch = @patch function preview_agents(a...; b...)
    return Dict{String, Any}("agentPreviews" => Any[Dict{String, Any}("agentId" => "i-49113b93")], "nextToken" => "1")
end

remove_attributes_from_findings_1_patch = @patch function remove_attributes_from_findings(a...; b...)
    return Dict{String, Any}("failedItems" => Dict{String, Any}())
end

start_assessment_run_1_patch = @patch function start_assessment_run(a...; b...)
    return Dict{String, Any}("assessmentRunArn" => "arn:aws:inspector:us-west-2:123456789012:target/0-nvgVhaxX/template/0-it5r2S4T/run/0-jOoroxyY")
end


text:
using Revise # reduce recompilation time
using Plots
using LinearAlgebra
using BenchmarkTools
using UnPack

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D
using Basis2DTri
using UniformTriMesh

using SetupDG

push!(LOAD_PATH, "./examples/EntropyStableEuler")
using EntropyStableEuler

"Approximation parameters"
N = 2 # The order of approximation
K1D = 3
CFL = 2 # CFL goes up to 2.5ish
T = .01 #1.0 # endtime

"Mesh related variables"
Kx = convert(Int,4/3*K1D)
Ky = K1D
VX, VY, EToV = uniform_tri_mesh(Kx, Ky)
@. VX = 15*(1+VX)/2
@. VY = 5*VY

# initialize ref element and mesh
# specify lobatto nodes to automatically get DG-SEM mass lumping
rd = init_reference_tri(N)
md = init_mesh((VX,VY),EToV,rd)

# Make domain periodic
@unpack Nfaces,Vf = rd
@unpack xf,yf,K,mapM,mapP,mapB = md
LX,LY = (x->maximum(x)-minimum(x)).((VX,VY)) # find lengths of domain
mapPB = build_periodic_boundary_maps(xf,yf,LX,LY,Nfaces*K,mapM,mapP,mapB)
mapP[mapB] = mapPB
@pack! md = mapP

# construct hybridized SBP operators
@unpack M,Dr,Ds,Vq,Pq,Vf,wf,nrJ,nsJ = rd
Qr = Pq'*M*Dr*Pq
Qs = Pq'*M*Ds*Pq
Ef = Vf*Pq
Br = diagm(wf.*nrJ)
Bs = diagm(wf.*nsJ)
Qrh = .5*[Qr-Qr' Ef'*Br;
         -Br*Ef  Br]
Qsh = .5*[Qs-Qs' Ef'*Bs;
         -Bs*Ef  Bs]

Vh = [Vq;Vf]
Ph = M\transpose(Vh)
VhP = Vh*Pq

# make sparse skew symmetric versions of the operators"
# precompute union of sparse ids for Qr, Qs
Qrhskew = .5*(Qrh-transpose(Qrh))
Qshskew = .5*(Qsh-transpose(Qsh))

# define initial conditions at nodes
@unpack x,y = md
rho,u,v,p = vortex(x,y,0)
Q = primitive_to_conservative(rho,u,v,p)

# interpolate geofacs to both vol/surf nodes
@unpack rxJ, sxJ, ryJ, syJ = md
rxJ, sxJ, ryJ, syJ = (x->Vh*x).((rxJ, sxJ, ryJ, syJ)) # interp to hybridized points
@pack! md = rxJ, sxJ, ryJ, syJ

# pack SBP operators into tuple
@unpack LIFT = rd
ops = (Qrhskew,Qshskew,VhP,Ph,LIFT,Vq)

#=
"Time integration"
rk4a,rk4b,rk4c = rk45_coeffs()
CN = (N+1)*(N+2)/2  # estimated trace constant for CFL
h  = 2/K1D
dt = CFL * h / CN
Nsteps = convert(Int,ceil(T/dt))
dt = T/Nsteps

"dense version - speed up by prealloc + transpose for col major "
function dense_hadamard_sum(Qhe,ops,vgeo,flux_fun)

    (Qr,Qs) = ops
    (rxJ,sxJ,ryJ,syJ) = vgeo
    # transpose for column-major evals
    QxTr = transpose(rxJ*Qr + sxJ*Qs)
    QyTr = transpose(ryJ*Qr + syJ*Qs)

    # precompute logs for logmean
    (rho,u,v,beta) = Qhe
    Qlog = (log.(rho), log.(beta))
    Qlogi = zeros(length(Qlog))
    Qlogj = zeros(length(Qlog))

    n = size(Qr,1)
    nfields = length(Qhe)

    QF = ntuple(x->zeros(n),nfields)
    QFi = zeros(nfields)
    for i = 1:n
        Qi = getindex.(Qhe,i)
        Qlogi = getindex.(Qlog,i)

        fill!(QFi,0)
        for j = 1:n
            Qj = getindex.(Qhe,j)
            Qlogj = getindex.(Qlog,j)

            Fx,Fy = flux_fun(Qi,Qj,Qlogi,Qlogj)
            @. QFi += QxTr[j,i]*Fx + QyTr[j,i]*Fy
        end

        for fld in eachindex(Qhe)
            QF[fld][i] = QFi[fld]
        end
    end

    return QF
end

"Qh = (rho,u,v,beta), while Uh = conservative vars"
function rhs(Q,md::MeshData,ops,flux_fun,compute_rhstest=false)

    @unpack rxJ,sxJ,ryJ,syJ,J,wJq = md
    @unpack nxJ,nyJ,sJ,mapP,mapB,K = md
    Qrh,Qsh,VhP,Ph,Lf,Vq = ops
    Nh,Nq = size(VhP)

    # entropy var projection
    VU = v_ufun((x->Vq*x).(Q)...)
    VU = (x->VhP*x).(VU)
    Uh = u_vfun(VU...)

    # convert to rho,u,v,beta vars
    (rho,rhou,rhov,E) = Uh
    beta = betafun(rho,rhou,rhov,E)
    Qh = (rho, rhou./rho, rhov./rho, beta) # redefine Q = (rho,u,v,β)

    # compute face values
    QM = (x->x[Nq+1:Nh,:]).(Qh)
    QP = (x->x[mapP]).(QM)

    # simple lax friedrichs dissipation
    Uf =  (x->x[Nq+1:Nh,:]).(Uh)
    (rhoM,rhouM,rhovM,EM) = Uf
    rhoUM_n = @. (rhouM*nxJ + rhovM*nyJ)/sJ
    lam = abs.(wavespeed(rhoM,rhoUM_n,EM))
    LFc = .25*max.(lam,lam[mapP]).*sJ

    fSx,fSy = flux_fun(QM,QP)
    normal_flux(fx,fy,u) = fx.*nxJ + fy.*nyJ - LFc.*(u[mapP]-u)
    flux = normal_flux.(fSx,fSy,Uf)
    rhsQ = (x->Lf*x).(flux)

    # compute volume contributions using flux differencing
    for e = 1:K
        Qhe = tuple(getindex.(Qh,:,e)...) # force tuples for fast splatting
        vgeo_local = getindex.((rxJ,sxJ,ryJ,syJ),1,e) # assumes affine elements for now

        Qops = (Qrh,Qsh)
        QFe = dense_hadamard_sum(Qhe,Qops,vgeo_local,flux_fun)

        mxm_accum!(X,x) = X[:,e] += 2*Ph*x
        mxm_accum!.(rhsQ,QFe)
    end

    rhsQ = (x -> -x./J).(rhsQ)

    rhstest = 0
    if compute_rhstest
        for fld in eachindex(rhsQ)
            VUq = VU[fld][1:Nq,:]
            rhstest += sum(wJq.*VUq.*(Vq*rhsQ[fld]))
        end
    end

    return rhsQ,rhstest # scale by Jacobian
end

Q = collect(Q) # make Q,resQ arrays of arrays for mutability
resQ = [zeros(size(x)) for i in eachindex(Q)]

# # testing
# rhsQ,rhstest = rhs(Q,md,ops,euler_fluxes,true)
# println("Testing: rhstest = $rhstest")

for i = 1:Nsteps

    rhstest = 0
    for INTRK = 1:5
        compute_rhstest = true # INTRK==5
        rhsQ,rhstest = rhs(Q,md,ops,euler_fluxes,compute_rhstest)

        @. resQ = rk4a[INTRK]*resQ + dt*rhsQ
        @. Q += rk4b[INTRK]*resQ
    end

    if i%10==0 || i==Nsteps
        println("Time step: $i out of $Nsteps with rhstest = $rhstest")
    end
end

# "project solution back to GLL nodes"
# Q = (x->Pq*x).(Q)

# use a higher degree quadrature for error evaluation
@unpack VDM = rd
@unpack J = md
rq2,sq2,wq2 = quad_nodes_2D(N+2)
Vq2 = vandermonde_2D(N,rq2,sq2)/VDM
wJq2 = diagm(wq2)*(Vq2*J)
xq2,yq2 = (x->Vq2*x).((x,y))

Qq = (x->Vq2*x).(Q)
Qex = primitive_to_conservative(vortex(xq2,yq2,T)...)
L2err = 0.0
for fld in eachindex(Q)
    global L2err
    L2err += sum(@. wJq2*(Qq[fld]-Qex[fld])^2)
end
L2err = sqrt(L2err)
println("L2err at final time T = $T is $L2err\n")

#plotting nodes
@unpack Vp = rd
gr(aspect_ratio=:equal,legend=false,
   markerstrokewidth=0,markersize=2)

vv = Vp*Q[1]
scatter(Vp*x,Vp*y,vv,zcolor=vv,camera=(0,90))
=#

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
include("../util/numeric.jl")

function solve_modal_eigen(structure,loadcase,restrainedDOFs::Vector{Int};path=pwd())
    @info "-------------------------- 求解自振模态特征值 --------------------------" 工况=loadcase.id 前一步基础工况=loadcase.plc
    K=structure.K
    P=loadcase.P
    u=zero(P)
    F₀=zero(P)
    if loadcase.plc!=""
        u=read_vector(path*"/.analysis",loadcase.plc*"_u.v")
        F₀=read_vector(path*"/.analysis",loadcase.plc*"_F.v")
        Kσ=calc_Kσ(structure,u)
        K=structure.K+Kσ #初始位移及初始刚度来自前一步结果
    end
    K̄=introduce_BC(K,restrainedDOFs)
    M̄=introduce_BC(structure.M,restrainedDOFs)

    nev=min(loadcase.nev,size(M̄,1))
    tol=loadcase.tolev
    maxiter=loadcase.maxiterev

    if size(K̄,1)>512
        ω²,ϕ=eigs(K̄,M̄,nev=nev,which=:SM,tol=tol,maxiter=maxiter)
    else
        eig=eigen(Array(K̄),Array(M̄))
        idx=collect(1:length(eig.values))
        sort!(idx,by=x->eig.values[x])
        ω²=zero(eig.values)
        ϕ=zero(eig.vectors)
        for i in 1:length(idx)
            ω²[i]=eig.values[idx[i]]
            ϕ[:,i]=eig.vectors[:,idx[i]]
        end
        ω²,ϕ=ω²[1:nev],ϕ[:,1:nev]
    end

    @info "------------------------------ 求解完成 ------------------------------"
    return ω²,ϕ
end

function solve_modal_Ritz(structure,loadcase,restrainedDOFs::Vector{Int};path=pwd())
    @info "-------------------------- 求解自振模态Ritz值 -------------------------" 工况=loadcase.id 前一步基础工况=loadcase.plc
    K=structure.K
    P=loadcase.P
    u=zero(P)
    F₀=zero(P)
    if loadcase.plc!=""
        u=read_vector(path*"/.analysis",loadcase.plc*"_u.v")
        F₀=read_vector(path*"/.analysis",loadcase.plc*"_F.v")
        Kσ=calc_Kσ(structure,u)
        K=structure.K+Kσ #初始位移及初始刚度来自前一步结果
    end
    K̄=introduce_BC(K,restrainedDOFs)
    F̄=introduce_BC(P,restrainedDOFs)
    M̄=introduce_BC(structure.M,restrainedDOFs)

    if F̄==zero(F̄)
        F̄.=1
    end

    nev=min(loadcase.nev,length(F̄))
    tol=loadcase.tolev
    maxiter=loadcase.maxiterev

    X=zeros(length(F̄),nev)
    if USE_PARDISO
        ps=Pardiso.MKLPardisoSolver()
        Pardiso.set_matrixtype!(ps,1) #实对称矩阵
        x̂=Pardiso.solve(ps,K̄,F̄)
        x̂=Pardiso.solve(ps,K̄,F̄)
    else
        x̂=Symmetric(K̄) \ F̄
    end
    β=sqrt(x̂'*M̄*x̂)
    X[:,1]=x̂/β
    for i in 2:nev
        if USE_PARDISO
            x̃=Pardiso.solve(ps,K̄,Array(M̄*X[:,i-1]))
        else
            x̃=Symmetric(K̄) \ Array(M̄*X[:,i-1])
        end
        x̂=copy(x̃)
        for j in 1:i-1
            α=x̃'*M̄*X[:,j]
            x̂-=α*X[:,j]
        end
        β=sqrt(x̂'*M̄*x̂)
        X[:,i]=x̂/β
    end
    ω²,ϕ=eigs(X'*K̄*X,X'*M̄*X,nev=nev,which=:SM,tol=tol,maxiter=maxiter) #X'*M̄*X=I?
    @info "------------------------------ 求解完成 ------------------------------"
    return ω²,ϕ
end

function solve_newmark_beta(structure,loadcase,restrainedDOFs::Vector{Int};path=pwd())
    T=loadcase.t
    Δt=(T[end]-T[1])/(length(T)-1)

    β=loadcase.β
    γ=loadcase.γ
    @info "----------------------- 求解Newmark-β法逐步积分 -----------------------" Δt=Δt β=β γ=γ
    K̄=introduce_BC(structure.K,restrainedDOFs)
    M̄=introduce_BC(structure.M,restrainedDOFs)
    C̄=introduce_BC(structure.C,restrainedDOFs)
    P̄=introduce_BC(loadcase.P,restrainedDOFs)
    Q̄=loadcase.f'.*P̄
    u₀=zeros(size(K̄,1))
    v₀=zeros(size(K̄,1))

    if !(β>=0.25*(0.5+γ)^2 && γ>=0.5)
        @warn "算法可能不稳定"
    end
    c₀=1/(β*Δt^2)
    c₁=γ/(β*Δt)
    c₂=1/(β*Δt)
    c₃=1/(2*β)-1
    c₄=γ/β-1
    c₅=Δt/2*(γ/β-2)
    c₆=Δt*(1-γ)
    c₇=γ*Δt

    K̂=K̄+c₀*M̄+c₁*C̄
    LDLᵀ=ldlt(Symmetric(K̂))

    ū=Array{Float64,2}(undef,length(u₀),length(T))
    v̄=Array{Float64,2}(undef,length(v₀),length(T))
    ā=Array{Float64,2}(undef,length(v₀),length(T))
    ū[:,1]=u₀
    v̄[:,1]=v₀
    if USE_PARDISO
        ps=Pardiso.MKLPardisoSolver()
        Pardiso.set_matrixtype!(ps,1)
        a₀=Pardiso.solve(ps,M̄,(Q̄[:,1]-C̄*v₀-K̄*u₀)[:,1])
    else
        a₀=Symmetric(M̄) \ (Q̄[:,1]-C̄*v₀-K̄*u₀)[:,1]
    end

    ā[:,1]=a₀
    for t in 1:length(T)-1
        Q̂=Q̄[:,t+1]+M̄*(c₀*ū[:,t]+c₂*v̄[:,t]+c₃*ā[:,t])+C̄*(c₁*ū[:,t]+c₄*v̄[:,t]+c₅*ā[:,t])
        Q̂=reshape(Q̂,length(Q̂))
        if USE_PARDISO
            ū[:,t+1]=Pardiso.solve(ps,K̂,Q̂) #May LDLT here
        else
            ū[:,t+1]=Symmetric(K̂) \ Q̂
        end
        ā[:,t+1]=c₀*(ū[:,t+1]-ū[:,t])-c₂*v̄[:,t]-c₃*ā[:,t]
        v̄[:,t+1]=v̄[:,t]+c₆*ā[:,t]+c₇*ā[:,t+1]
    end
    @info "------------------------------ 求解完成 ------------------------------"

    u=zeros(size(structure.K,1),length(T))
    v=zeros(size(structure.K,1),length(T))
    a=zeros(size(structure.K,1),length(T))
    for t in 1:length(T)
        u[:,t]=resolve_BC(ū[:,t],restrainedDOFs)
        v[:,t]=resolve_BC(v̄[:,t],restrainedDOFs)
        a[:,t]=resolve_BC(ā[:,t],restrainedDOFs)
    end
    return u,v,a
end

function solve_wilson_theta(structure,loadcase,restrainedDOFs::Vector{Int};path=pwd())
    @info "----------------------- 求解Wilson-θ法逐步积分 -----------------------"
    T=loadcase.t
    Δt=(T[end]-T[1])/(length(T)-1)

    β=loadcase.β
    γ=loadcase.γ
    θ=loadcase.θ

    if !(θ<1.37)
        @warn "算法可能不稳定"
    end

    K̄=introduce_BC(structure.K,restrainedDOFs)
    M̄=introduce_BC(structure.M,restrainedDOFs)
    C̄=introduce_BC(structure.C,restrainedDOFs)
    P̄=introduce_BC(loadcase.P,restrainedDOFs)
    Q̄=loadcase.f'.*P̄
    u₀=zeros(size(K̄,1))
    v₀=zeros(size(K̄,1))

    Δt̄=θ*Δt
    b₁=1/(β*Δt^2)
    b₂=-1/(β*Δt̄)
    b₃=(1/2-β)/β
    b₄=γ*Δt̄*b₁
    b₅=1+γ*Δt̄*b₂
    b₆=Δt̄*(1+γ*b₃-γ)

    K̂=K̄+b₁*M̄+b₄*C̄
    LDLᵀ=ldlt(Symmetric(K̂))

    ū=Array{Float64,2}(undef,length(u₀),length(T))
    v̄=Array{Float64,2}(undef,length(v₀),length(T))
    ā=Array{Float64,2}(undef,length(v₀),length(T))
    ū[:,1]=u₀
    v̄[:,1]=v₀

    if USE_PARDISO
        ps=Pardiso.MKLPardisoSolver()
        Pardiso.set_matrixtype!(ps,1)
        a₀=Pardiso.solve(ps,M̄,(Q̄[:,1]-C̄*v₀-K̄*u₀)[:,1])
    else
        a₀=Symmetric(M̄) \ (Q̄[:,1]-C̄*v₀-K̄*u₀)[:,1]
    end

    ā[:,1]=a₀
    for t in 1:length(T)-1
        Q̂=Q̄[:,t+1]+M̄*(b₁*ū[:,t]+b₂*v̄[:,t]+b₃*ā[:,t])+C̄*(b₄*ū[:,t]+b₅*v̄[:,t]+b₆*ā[:,t])
        Q̂=reshape(Q̂,length(Q̂))
        if USE_PARDISO
            ū[:,t+1]=Pardiso.solve(ps,K̂,Q̂) #May LDLT here
        else
            ū[:,t+1]=Symmetric(K̂) \ Q̂
        end
        v̄[:,t+1]=b₄*(ū[:,t+1]-ū[:,t])+b₅*v̄[:,t]+b₆*ā[:,t]
        ā[:,t+1]=b₁*(ū[:,t+1]-ū[:,t])+b₂*v̄[:,t]+b₃*ā[:,t]

        ā[:,t+1]=ā[:,t]+1/θ*(ā[:,t+1]-ā[:,t])
        v̄[:,t+1]=v̄[:,t]+((1-γ)*ā[:,t]+γ*ā[:,t+1])*Δt
        ū[:,t+1]=ū[:,t]+v̄[:,t]*Δt
    end
    @info "------------------------------ 求解完成 ------------------------------"

    u=zeros(size(structure.K,1),length(T))
    v=zeros(size(structure.K,1),length(T))
    a=zeros(size(structure.K,1),length(T))
    for t in 1:length(T)
        u[:,t]=resolve_BC(ū[:,t],restrainedDOFs)
        v[:,t]=resolve_BC(v̄[:,t],restrainedDOFs)
        a[:,t]=resolve_BC(ā[:,t],restrainedDOFs)
    end
    return u,v,a
end

#not finished

function solve_HHT_alpha(structure,loadcase,restrainedDOFs::Vector{Int};path=pwd())
    T=loadcase.t
    Δt=(T[end]-T[1])/(length(T)-1)
    α=loadcase.α
    β=(2-α^2)/4
    γ=1.5-α
    @info "----------------------- 求解Hilber-Hughes-Taylor法逐步积分 -----------------------"
    K̄=introduce_BC(structure.K,restrainedDOFs)
    M̄=introduce_BC(structure.M,restrainedDOFs)
    C̄=introduce_BC(structure.C,restrainedDOFs)
    P̄=introduce_BC(loadcase.P,restrainedDOFs)
    Q̄=loadcase.f'.*P̄
    u₀=zeros(size(K̄,1))
    v₀=zeros(size(K̄,1))

    if !(0.5<α<1.0)
        @warn "算法可能不精确"
    end
    c₀=1/(β*Δt^2)
    c₁=α*γ/(β*Δt)
    c₂=1/(β*Δt)
    c₃=1/(2*β)-1
    c₄=γ/β-1
    c₅=Δt/2*(γ/β-2)
    c₆=Δt*(1-γ)
    c₇=γ*Δt

    K̂=α*K̄+c₀*M̄+c₁*C̄

    LDLᵀ=ldlt(Symmetric(K̂))

    ū=Array{Float64,2}(undef,length(u₀),length(T))
    v̄=Array{Float64,2}(undef,length(v₀),length(T))
    ā=Array{Float64,2}(undef,length(v₀),length(T))
    ū[:,1]=u₀
    v̄[:,1]=v₀
    if USE_PARDISO
        ps=Pardiso.MKLPardisoSolver()
        Pardiso.set_matrixtype!(ps,1)
        a₀=Pardiso.solve(ps,M̄,(Q̄[:,1]-C̄*v₀-K̄*u₀)[:,1])
    else
        a₀=Symmetric(M̄) \ (Q̄[:,1]-C̄*v₀-K̄*u₀)[:,1]
    end

    ā[:,1]=a₀
    for t in 1:length(T)-1
        Q̂=Q̄[:,t+1]+M̄*(c₀*ū[:,t]+c₂*v̄[:,t]+c₃*ā[:,t])+C̄*(c₁*ū[:,t]+c₄*v̄[:,t]+c₅*ā[:,t])

        Q̂=reshape(Q̂,length(Q̂))
        if USE_PARDISO
            ū[:,t+1]=Pardiso.solve(ps,K̂,Q̂) #May LDLT here
        else
            ū[:,t+1]=Symmetric(K̂) \ Q̂
        end
        ā[:,t+1]=c₀*(ū[:,t+1]-ū[:,t])-c₂*v̄[:,t]-c₃*ā[:,t]
        v̄[:,t+1]=v̄[:,t]+c₆*ā[:,t]+c₇*ā[:,t+1]
    end
    @info "------------------------------ 求解完成 ------------------------------"

    u=zeros(size(structure.K,1),length(T))
    v=zeros(size(structure.K,1),length(T))
    a=zeros(size(structure.K,1),length(T))
    for t in 1:length(T)
        u[:,t]=resolve_BC(ū[:,t],restrainedDOFs)
        v[:,t]=resolve_BC(v̄[:,t],restrainedDOFs)
        a[:,t]=resolve_BC(ā[:,t],restrainedDOFs)
    end
    return u,v,a
end

function duhamel(ω,ζ,t,p)
    dt=(t[end]-t[1])/(length(t)-1)
    N=length(t)
    ω̄=ω*sqrt(1-ζ^2)
    q=exp.(-ζ*ω*t).*sin.(ω*t)
    d=(1/ω)*conv(p,q)[1:N]*dt
end

#差分代微分
function diffasdiff(x₀,y₀)
    d=diff(y₀,dims=2)./diff(x₀)'
    dydx=zero(y₀)
    dydx[:,1:end-1]=d
    dydx[:,2:end]+=d
    dydx[:,1]+=d[:,1]
    dydx[:,end]+=d[:,end]
    dydx./2
end

function solve_modal_decomposition(structure,loadcase,restrainedDOFs::Vector{Int};path=pwd())
    modal_case=loadcase.modal_case
    ω²=read_vector(path*"/.analysis",modal_case*"_o.v")
    ϕ=read_matrix(path*"/.analysis",modal_case*"_p.m")

    @info "------------------------ 求解振型分解法频域积分 ------------------------" 模态工况=modal_case 模态数量=length(ω²)

    K̄=introduce_BC(structure.K,restrainedDOFs)
    M̄=introduce_BC(structure.M,restrainedDOFs)
    C̄=introduce_BC(structure.C,restrainedDOFs)
    P̄=introduce_BC(loadcase.P,restrainedDOFs)
    T=loadcase.t #should be a nx1 vector
    P̄=loadcase.f'.*P̄ #should be a nfDOF*nt matrix
    ω̂=sqrt.(ω²)
    K̂=diag(ϕ'*K̄*ϕ)
    M̂=diag(ϕ'*M̄*ϕ)
    Ĉ=diag(ϕ'*C̄*ϕ)
    P̂=[(ϕ'*P̄)[i,:] for i in 1:size(ϕ,2)]
    X=zeros(length(ω²),length(loadcase.f))
    for i in 1:length(ω²)
        ζ=Ĉ[i]/(2*sqrt(M̂[i]*K̂[i]))
        X[i,:]=duhamel(ω̂[i],ζ,T,P̂[i])
    end
    ū=ϕ*X
    v̄=diffasdiff(T,ū)
    ā=diffasdiff(T,v̄)
    @info "------------------------------ 求解完成 ------------------------------"
    u=zeros(size(structure.K,1),length(T))
    v=zeros(size(structure.K,1),length(T))
    a=zeros(size(structure.K,1),length(T))
    for t in 1:length(T)
        u[:,t]=resolve_BC(ū[:,t],restrainedDOFs)
        v[:,t]=resolve_BC(v̄[:,t],restrainedDOFs)
        a[:,t]=resolve_BC(ā[:,t],restrainedDOFs)
    end
    return u,v,a
end

function solve_central_diff(structure,loadcase,restrainedDOFs::Vector{Int};path=pwd(),Δtcr=0)
    T=loadcase.t
    Δt=(T[end]-T[1])/(length(T)-1)
    @info "----------------------- 求解中心差分法逐步积分 -----------------------" Δt=Δt
    if Δtcr==0
        @warn "未考虑临界Δt，算法可能不稳定"
    end
    c₀=1/Δt^2
    c₁=1/2Δt
    c₂=2c₀
    c₃=1/c₂

    K̄=introduce_BC(structure.K,restrainedDOFs)
    M̄=introduce_BC(structure.M,restrainedDOFs)
    C̄=introduce_BC(structure.C,restrainedDOFs)
    P̄=introduce_BC(loadcase.P,restrainedDOFs)
    Q̄=loadcase.f'.*P̄

    u₀=zeros(size(K̄,1))
    v₀=zeros(size(K̄,1))
    a₀=zeros(size(K̄,1))

    ū=Array{Float64,2}(undef,length(u₀),length(T)+1)
    v̄=Array{Float64,2}(undef,length(v₀),length(T)+1)
    ā=Array{Float64,2}(undef,length(v₀),length(T)+1)

    ū[:,2]=u₀
    v̄[:,2]=v₀
    ā[:,2]=a₀

    ū[:,1]=u₀-Δt*v₀+c₃*a₀
    M̂=c₀*M̄+c₁*C̄

    LDLᵀ=ldlt(Symmetric(M̂))
    if USE_PARDISO
        ps=Pardiso.MKLPardisoSolver()
        Pardiso.set_matrixtype!(ps,1)
    end
    for t in 2:length(T)
        Q̂=Q̄[:,t]-(K̄-c₂*M̄)*ū[:,t]-(c₀*M̄-c₁*C̄)*ū[:,t-1]
        if USE_PARDISO
            ū[:,t+1]=Pardiso.solve(ps,M̂,Q̂)
        else
            ū[:,t+1]=Symmetric(M̂) \ Q̂
        end
        v̄[:,t]=c₁*(-ū[:,t-1]+ū[:,t+1])
        ā[:,t]=c₀*(ū[:,t-1]-2*ū[:,t]+ū[:,t+1])
    end
    @info "------------------------------ 求解完成 ------------------------------"

    u=zeros(size(structure.K,1),length(T))
    v=zeros(size(structure.K,1),length(T))
    a=zeros(size(structure.K,1),length(T))
    for t in 1:length(T)
        u[:,t]=resolve_BC(ū[:,t+1],restrainedDOFs)
        v[:,t]=resolve_BC(v̄[:,t+1],restrainedDOFs)
        a[:,t]=resolve_BC(ā[:,t+1],restrainedDOFs)
    end
    return u,v,a
end

function solve_response_spectrum(structure,loadcase,restrainedDOFs::Vector{Int};path=pwd())
    modal_case=loadcase.modal_case
    ω²=read_vector(path*"/.analysis",modal_case*"_o.v")
    ϕ=read_matrix(path*"/.analysis",modal_case*"_p.m")

    K̄=introduce_BC(structure.K,restrainedDOFs)
    M̄=introduce_BC(structure.M,restrainedDOFs)
    C̄=introduce_BC(structure.C,restrainedDOFs)
    P̄=introduce_BC(loadcase.P,restrainedDOFs)
    T=loadcase.t #should be a nx1 vector
    P̄=loadcase.f'.*P̄ #should be a nfDOF*nt matrix
    ω̂=sqrt.(ω²)
    K̂=diag(ϕ'*K̄*ϕ)
    M̂=diag(ϕ'*M̄*ϕ)
    Ĉ=diag(ϕ'*C̄*ϕ)
    P̂=[(ϕ'*P̄)[i,:] for i in 1:size(ϕ,2)]
    X=zeros(length(ω²),length(loadcase.f))
    for i in 1:length(ω²)
        ζ=Ĉ[i]/(2*sqrt(M̂[i]*K̂[i]))
        X[i,:]=duhamel(ω̂[i],ζ,T,P̂[i])
    end
    ū=ϕ*X
    v̄=ū
    ā=ū
    @info "------------------------------ 求解完成 ------------------------------"
    u=zeros(size(structure.K,1),length(T))
    v=zeros(size(structure.K,1),length(T))
    a=zeros(size(structure.K,1),length(T))
    for t in 1:length(T)
        u[:,t]=resolve_BC(ū[:,t],restrainedDOFs)
        v[:,t]=resolve_BC(v̄[:,t],restrainedDOFs)
        a[:,t]=resolve_BC(ā[:,t],restrainedDOFs)
    end
    return u,v,a
end

#
# """
# spec: a {'T':period,'a':acceleration} dictionary of spectrum\n
# mdd: a list of modal damping ratio\n
# comb: combination method, 'CQC' or 'SRSS'
# """
# function solve_response_spectrum(assembly:assembly,T::Vector,α::Vector,n=60,comb='CQC'):
#     K=assembly.K̄
#     M=assembly.M̄
#     DOF=size(K̄)[1]
#     w,f,T,mode=eigen_mode(assembly,DOF)
#     mode=mode[:,1:n]#use n modes only.
#     M̂=ϕ'*M*ϕ#generalized mass
#     K̂=ϕ'*K*ϕ#generalized stiffness
#     L̂=M*ϕ
#     px=[]
#     Vx=[]
#     Xm=[]
#     γ=[]
#     mx=diag(M)
#     for i in range(len(mode)):
#         #mass participate factor
#         push!(px,-(ϕ[:,i])'*mx)
#         push!(Vx,px[end]^2)
#         push!(Xm,Vx[end]/3/m)
#         #modal participate factor
#         push!(γ,L̂[i]/M̂[i,i])
#     end
#     S=zeros(DOF,n)
#
#     for i in 1:n:
#         ξ=assembly.C̄[i]
#         y=linear_interp(T[i],T,α)
#         y/=ω²[i]
#         S[:,i]=γ[i]*y*mode[:,i]
#     end
#
#     if comb=='CQC':
#         cqc=0
#         rho=Matrix(,n,n)
#         for i in range(mode.shape[1]):
#             for j in range(mode.shape[1]):
#                 if i!=j:
#                     r=T[i]/T[j]
#                     rho[i,j]=8*xi^2*(1+r)*r^1.5/((1-r^2)^2+4*xi^2*r*(1+r)^2)
#                 end
#                 cqc+=rho[i,j]*S[:,i]*S[:,j]
#             end
#         end
#         cqc=sqrt(cqc)
#         print(cqc)
#     elseif comb=='SRSS':
#         SRSS(S)
#     end
# end
#
#


text:
using Test, StructTypes, Tables, Strapping

tbl = (a=[10], b=[3.14])
tbl2 = (id=[10, 10, 10], floats=[3.14, 3.15, 3.16])

struct AB
    a::Int
    b::Float64
end
StructTypes.StructType(::Type{AB}) = StructTypes.Struct()

ab = AB(10, 3.14)
@test Strapping.construct(AB, tbl) == ab
@test Strapping.construct(Vector{AB}, tbl) == [ab]

mutable struct ABM
    a::Int
    b::Float64
    ABM() = new()
    ABM(a::Int, b::Float64) = new(a, b)
end
Base.:(==)(a::ABM, b::ABM) = a.a == b.a && a.b == b.b
StructTypes.StructType(::Type{ABM}) = StructTypes.Mutable()

abm = ABM(10, 3.14)
@test Strapping.construct(ABM, tbl) == abm
@test Strapping.construct(Vector{ABM}, tbl) == [abm]

@test Strapping.construct(NamedTuple, tbl) == (a=10, b=3.14)
@test Strapping.construct(NamedTuple{(:a, :b), Tuple{Int, Float64}}, tbl) == (a=10, b=3.14)
@test Strapping.construct(NamedTuple{(:a, :b)}, tbl) == (a=10, b=3.14)
@test Strapping.construct(Dict, tbl) == Dict("a" => 10, "b" => 3.14)
@test Strapping.construct(Dict{String, Float64}, tbl) == Dict("a" => 10.0, "b" => 3.14)
@test Strapping.construct(Dict{Symbol, Any}, tbl) == Dict(:a => 10, :b => 3.14)
@test Strapping.construct(Array, tbl) == [10, 3.14]
@test Strapping.construct(Vector{Float64}, tbl) == [10.0] # because Float64 is scalar, only 1st field of result is used, other fields are ignored
@test Strapping.construct(Set, tbl) == Set([10, 3.14])
@test Strapping.construct(Tuple, tbl) == (10, 3.14)
@test Strapping.construct(Tuple{Int, Float64}, tbl) == (10, 3.14)
@test Strapping.construct(Vector{Any}, tbl) == [Dict("a" => 10, "b" => 3.14)]

struct AB2
    id::Int
    floats::Vector{Float64}
end

Base.:(==)(a::AB2, b::AB2) = a.id == b.id && a.floats == b.floats
StructTypes.StructType(::Type{AB2}) = StructTypes.Struct()
StructTypes.idproperty(::Type{AB2}) = :id

ab2 = AB2(10, [3.14, 3.15, 3.16])
@test Strapping.construct(AB2, tbl2) == ab2
@test Strapping.construct(Vector{AB2}, tbl2) == [ab2]

struct AB3
    a::Int
    ab::AB
end

StructTypes.StructType(::Type{AB3}) = StructTypes.Struct()

tbl3 = (a=[10], ab_a=[10], ab_b=[3.14])

ab3 = AB3(10, AB(10, 3.14))
@test Strapping.construct(AB3, tbl3) == ab3
@test Strapping.construct(Vector{AB3}, tbl3) == [ab3]

mutable struct AB4
    a::Int
    b::Float64
    abs::Vector{AB}
    AB4() = new()
    AB4(a::Int, b::Float64, abs::Vector{AB}) = new(a, b, abs)
end

Base.:(==)(a::AB4, b::AB4) = a.a == b.a && a.b == b.b && a.abs == b.abs
StructTypes.StructType(::Type{AB4}) = StructTypes.Mutable()
StructTypes.idproperty(::Type{AB4}) = :a

tbl4 = (a=[10, 10, 10], b=[3.14, 3.14, 3.14], abs_a=[10, 10, 10], abs_b=[3.14, 3.15, 3.16])

ab4 = AB4(10, 3.14, AB[AB(10, 3.14), AB(10, 3.15), AB(10, 3.16)])
@test Strapping.construct(AB4, tbl4) == ab4
@test Strapping.construct(Vector{AB4}, tbl4) == [ab4]

# multiple levels of nesting
struct AB5
    a::Int
    ab::AB3
end

StructTypes.StructType(::Type{AB5}) = StructTypes.Struct()

tbl5 = (a=[10], ab_a=[10], ab_ab_a=[10], ab_ab_b=[3.14])

ab5 = AB5(10, AB3(10, AB(10, 3.14)))
@test Strapping.construct(AB5, tbl5) == ab5
@test Strapping.construct(Vector{AB5}, tbl5) == [ab5]

# deconstruction
@test columntable(Strapping.deconstruct(ab)) == tbl
@test columntable(Strapping.deconstruct(ab2)) == tbl2
@test columntable(Strapping.deconstruct(ab3)) == tbl3
@test columntable(Strapping.deconstruct(ab4)) == tbl4
@test columntable(Strapping.deconstruct(ab5)) == tbl5

struct AB6
    a::Int
    b::Vector{Float64}
    c::AB
    d::Vector{AB}
end

Base.:(==)(a::AB6, b::AB6) = a.a == b.a && a.b == b.b && a.c == b.c && a.d == b.d
StructTypes.StructType(::Type{AB6}) = StructTypes.Struct()
StructTypes.idproperty(::Type{AB6}) = :a

tbl6 = (a=[1, 1, 1], b=[3.14, 3.15, 3.16], c_a=[2, 2, 2], c_b=[0.01, 0.01, 0.01], d_a=[10, 11, 12], d_b=[1.1, 2.2, 3.3])

ab6 = AB6(1, [3.14, 3.15, 3.16], AB(2, 0.01), [AB(10, 1.1), AB(11, 2.2), AB(12, 3.3)])
@test Strapping.construct(AB6, tbl6) == ab6
@test Strapping.construct(Vector{AB6}, tbl6) == [ab6]
@test columntable(Strapping.deconstruct(ab6)) == tbl6

# https://github.com/JuliaData/Strapping.jl/issues/3
struct TestStruct
    a::Float64
    b::Float64
    id::Int
end

StructTypes.StructType(::Type{TestStruct}) = StructTypes.Struct()
StructTypes.idproperty(::Type{TestStruct}) = :id

data = [ TestStruct(rand(2)..., n) for n = 1:5]

tbl = Strapping.deconstruct(data)
@test length(Tables.columntable(tbl)[1]) == 5

struct AA1
    a::Int64
    b::Dict{Symbol, Any}
end

StructTypes.StructType(::Type{AA1}) = StructTypes.Struct()

data = [ AA1(1, Dict{Symbol, Any}(:aa => 2, :bb => 3)) ]
tbl = Strapping.deconstruct(data)
tbl2 = Tables.columntable(tbl)
@test tbl2.a == [1]
@test tbl2.b_aa == [2]
@test tbl2.b_bb == [3]


text:
__precompile__(true)
module ParticleSwarmOptimizer

include("TestFunctions.jl")
include("helper.jl")
include("neighbourhood.jl")
include("pso.jl")

export TestFunctions

end # module