text:
function OffPolicyGAIL(;π,
                        S,
                        𝒟_demo,
                        𝒟_ndas::Array{ExperienceBuffer} = ExperienceBuffer[],
                        normalize_demo::Bool=true,
                        D::ContinuousNetwork,
                        solver=SAC,
                        d_opt::NamedTuple=(epochs=5,),
                        log::NamedTuple=(;),
                        kwargs...)

    # Define the training parameters for the desciminator
    d_opt = TrainingParams(;name="discriminator_", loss=()->nothing, d_opt...)

    # Setup NDA parameters
    N_nda = length(𝒟_ndas)
    λ_nda = Float32(-1 / N_nda)
    N_datasets = 2 + N_nda

    # Normalize and/or change device of expert and NDA data
    dev = device(π)
    A = action_space(π)
    normalize_demo && (𝒟_demo = normalize!(deepcopy(𝒟_demo), S, A))
    𝒟_demo = 𝒟_demo |> dev
    for i in 1:N_nda
        𝒟_ndas[i] = normalize_demo && normalize!(deepcopy(𝒟_ndas[i]), S, A)
        𝒟_ndas[i] = 𝒟_ndas[i] |> dev
    end

    # Build the solver
    𝒮 = solver(;π=π,
                S=S,
                log=(dir="log/offpolicygail", period=500, log...),
                kwargs...)

    # Setup the training of the discriminator
    B = d_opt.batch_size

    # These are minibatch buffers
    𝒟_batch = buffer_like(𝒮.buffer, capacity=B, device=dev)
    𝒟_demo_batch = deepcopy(𝒟_batch)
    𝒟_demo_π_batch = deepcopy(𝒟_batch)
    𝒟_ndas_batch = [deepcopy(𝒟_batch) for 𝒟_nda in 𝒟_ndas]
    𝒟_ndas_π_batch = [deepcopy(𝒟_batch) for 𝒟_nda in 𝒟_ndas]

    function GAIL_callback(𝒟; info=Dict(), kwargs...)
        for i=1:d_opt.epochs

            # Sample minibatchs
            rand!(𝒟_demo_batch, 𝒟_demo)
            # rand!(𝒟_demo_π_batch, 𝒟_demo)
            # 𝒟_demo_π_batch.data[:a] = action(π, 𝒟_demo_π_batch[:s])
            rand!(𝒟_batch, 𝒮.buffer)
            for i in 1:N_nda
                rand!(𝒟_ndas_batch[i], 𝒟_ndas[i])
                # rand!(𝒟_ndas_π_batch[i], 𝒟_ndas[i])
                # 𝒟_ndas_π_batch[i].data[:a] = action(π, 𝒟_ndas_π_batch[i][:s])
            end
            # Concatenate minibatches into on buffer
            𝒟_full = hcat(𝒟_demo_batch, 𝒟_batch, 𝒟_ndas_batch...)

            # concat inputs
            x = cat(flatten(𝒟_full[:s]), 𝒟_full[:a], dims=1)

            # Add some noise (Optional)
            # x .+= Float32.(rand(Normal(0, 0.2f0), size(x))) |> dev

            # Create labels
            y_demo = Flux.onehotbatch(ones(Int, B), 1:N_datasets)
            # y_demo_π = Flux.onehotbatch(2*ones(Int, B), 1:N_datasets)
            y_π = Flux.onehotbatch(2*ones(Int, B), 1:N_datasets)
            y_ndas = [Flux.onehotbatch((i+2)*ones(Int, B), 1:N_datasets) for i=1:N_nda]
            # y_ndas_π = [Flux.onehotbatch(2*ones(Int, B), 1:N_datasets) for i=1:N_nda]

            y = cat(y_demo, y_π, y_ndas..., dims=2) |> dev

            # println(gradient_penalty(D, x_demo, x_π))
            # + 10f0 * gradient_penalty(D, x_demo, x_π)
            train!(D, (;kwargs...) -> Flux.Losses.logitcrossentropy(D(x), y), d_opt, info=info)
        end

        # ## replace the bufffer
        # rand!(𝒟_demo_batch, 𝒟_demo)
        # # rand!(𝒟_demo_π_batch, 𝒟_demo)
        # # 𝒟_demo_π_batch.data[:a] = action(π, 𝒟_demo_π_batch[:s])
        # rand!(𝒟_batch, 𝒮.buffer)
        # for i in 1:N_nda
        #     rand!(𝒟_ndas_batch[i], 𝒟_ndas[i])
        #     # rand!(𝒟_ndas_π_batch[i], 𝒟_ndas[i])
        #     # 𝒟_ndas_π_batch[i].data[:a] = action(π, 𝒟_ndas_π_batch[i][:s])
        # end
        #
        # 𝒟_full = hcat(𝒟_demo_batch, 𝒟_batch, 𝒟_ndas_batch...)
        #
        # for k in keys(𝒟)
        #     𝒟.data[k] = 𝒟_full.data[k]
        # end
        # 𝒟.elements = 𝒟_full.elements
        # 𝒟.next_ind = 𝒟_full.next_ind

        ## Compute the rewards
        D_out = Flux.softmax(value(D, flatten(𝒟[:s]), 𝒟[:a]))
        w = [1f0, 0f0, λ_nda*ones(Float32, N_nda)...] |> dev

        𝒟[:r] .= sum((Base.log.(D_out .+ 1f-5) .- Base.log.(1f0 .- D_out .+ 1f-5)) .* w, dims=1)
    end

    𝒮.post_batch_callback = GAIL_callback
    𝒮
end


text:
using TightlyBound
using Test
using Suppressor

TESTDIR=TightlyBound.TESTDIR

@testset "test atomicproj" begin

    @suppress begin

        units_old = TightlyBound.set_units()
        TightlyBound.set_units(both="atomic")

        filname = "$TESTDIR/data_forces/znse.in_vnscf_vol_2/projham.xml.gz"
        tbc_ref = read_tb_crys(filname)

        fil = "$TESTDIR/data_forces/znse.in_vnscf_vol_2/"
        dft = TightlyBound.QE.loadXML(fil*"/qe.save/")

        tbc, tbck, proj_warn = TightlyBound.AtomicProj.projwfx_workf(dft, directory=fil, nprocs=1, writefile=missing, writefilek=missing, skip_og=true, skip_proj=true, freeze=true, localized_factor = 0.15, only_kspace = false, screening = 1.0 )

        energy_tot, tbc_, conv_flag = scf_energy(tbc)

        @test (energy_tot - dft.atomize_energy) < 1e-4

        energy_tot_ref, tbc_, conv_flag = scf_energy(tbc_ref)
        @test (energy_tot_ref - dft.atomize_energy) < 1e-4

        TightlyBound.set_units(energy = units_old[1], length=units_old[2])

    end

end

ext: jl


text:
using HTTP
using Sockets
using Logging

global_logger(NullLogger()) # Shut off logging

const HOST = "127.0.0.1"
const PORT = 8080
const STATIC_DIR = "./static/"

function file_server()
    web_server = HTTP.Servers.Server(handler, ratelimit=1000000//1)
    @async HTTP.Servers.serve(web_server, HOST, PORT, false)
end

function handler(req::HTTP.Request)
    (req.method != "GET") && return HTTP.Response(404, "Method $(req.method) not supported!")
    try
        return HTTP.Response(get_file(req))
    catch e
        return HTTP.Response(404, "Error: $e")
    end
end

function get_file(req::HTTP.Request) :: Vector{UInt8}
    t = (req.target == "/") ? "/index.html" : req.target
    read(STATIC_DIR * t)
end

# An alternate implementation
function __start()
    n::Int64 = 0
    errors::Int64 = 0
    @async begin
        HTTP.listen(pipeline_limit = 1000000) do request::HTTP.Request
            n += 1
            try
                @info "$n requests with $errors failures."
                return HTTP.Response("Hello There!")
            catch e
                errors += 1
                @info "$n requests with $errors failures."
                return HTTP.Response(404, "Error: $e")
            end
        end
    end
end


text:
# # Classification of penguin species
#
# ## Packages

using CairoMakie
using CalibrationErrors
using DataFrames
using Distances
using Distributions
using PalmerPenguins
using Query

using Random

using CairoMakie.AbstractPlotting.ColorSchemes: Dark2_8

# ## Data
#
# In this example we study the calibration of different models that classify three penguin
# species based on measurements of their bill and flipper lengths.
#
# We use the [Palmer penguins dataset](https://allisonhorst.github.io/palmerpenguins/) to
# to train and validate the models.

penguins = dropmissing(DataFrame(PalmerPenguins.load()))

f = Figure()
ax = Axis(f[1, 1]; xlabel="bill length [mm]", ylabel="flipper length [mm]")
for (i, (key, df)) in enumerate(pairs(groupby(penguins, :species)))
    scatter!(
        df.bill_length_mm, df.flipper_length_mm; color=(Dark2_8[i], 0.8), label=key.species
    )
end
Legend(f[1, 2], ax, "species")
#!jl save("./figures/penguins.svg", f);

#!jl # ![](./figures/penguins.svg)

# We split the data randomly into a training and validation dataset. The training dataset
# contains around 60% of the samples.

Random.seed!(1234)
idxs = shuffle(axes(penguins, 1))
k = floor(Int, 0.6 * length(idxs))
train_idxs = @view idxs[1:k]
val_idxs = @view idxs[(k + 1):end]

train_penguins = penguins[train_idxs, :]
val_penguins = penguins[val_idxs, :];

#-

f = Figure()
ax = Axis(f[1, 1]; xlabel="bill length [mm]", ylabel="flipper length [mm]")
for (i, df) in enumerate((train_penguins, val_penguins))
    for (j, (key, subdf)) in enumerate(pairs(groupby(df, :species)))
        scatter!(
            subdf.bill_length_mm,
            subdf.flipper_length_mm;
            color=(Dark2_8[j], 0.8),
            marker=i == 1 ? :circle : :diamond,
        )
    end
end
group_marker = [
    MarkerElement(; marker=m, color=:black, strokecolor=:transparent, markersize=20) for
    m in (:circle, :diamond)
]
group_color = [PolyElement(; color=Dark2_8[i], strokecolor=:transparent) for i in 1:3]
Legend(
    f[1, 2],
    [group_marker, group_color],
    [["training", "validation"], string.(levels(penguins.species))],
    ["dataset", "species"],
)
#!jl save("./figures/penguins_datasets.svg", f);

#!jl # ![](./figures/penguins_datasets.svg)

# ## Fitting normal distributions
#
# For each species, we fit independent normal distributions to the observations of the bill
# and flipper length in the training data, using maximum likelihood estimation.

penguins_fit = @from i in train_penguins begin
    @group i by i.species into g
    @select {
        species = key(g),
        proportion = length(g) / nrow(train_penguins),
        bill = fit(Normal, g.bill_length_mm),
        flipper = fit(Normal, g.flipper_length_mm),
    }
    @collect DataFrame
end

# We plot the estimated normal distributions.

function xrange(dists, alpha=0.0001)
    xmin = minimum(Base.Fix2(quantile, alpha), dists)
    xmax = maximum(Base.Fix2(quantile, 1 - alpha), dists)
    return range(xmin, xmax; length=1_000)
end

function plot_normal_fit(dists, species, xlabel)
    f = Figure()
    Axis(f[1, 1]; xlabel=xlabel, ylabel="density")
    xs = xrange(dists)
    plots = map(enumerate(dists)) do (i, dist)
        ys = pdf.(dist, xs)
        l = lines!(xs, ys; color=Dark2_8[i])
        b = band!(xs, 0, ys; color=(Dark2_8[i], 0.2))
        return [l, b]
    end
    Legend(f[1, 2], plots, species, "species")
    return f
end

plot_normal_fit(penguins_fit.bill, penguins_fit.species, "bill length [mm]")
#!jl save("./figures/normal_fit_bill.svg", current_figure());

#!jl # ![](./figures/normal_fit_bill.svg)

plot_normal_fit(penguins_fit.flipper, penguins_fit.species, "flipper length [mm]")
#!jl save("./figures/normal_fit_flipper.svg", current_figure());

#!jl # ![](./figures/normal_fit_flipper.svg)

# ## Naive Bayes classifier
#
# Let us assume that the bill and flipper length are conditionally independent given the
# penguin species. Then Bayes' theorem implies that
# ```math
# \begin{aligned}
# \mathbb{P}(\mathrm{species} \,|\, \mathrm{bill}, \mathrm{flipper})
# &= \frac{\mathbb{P}(\mathrm{species}) \mathbb{P}(\mathrm{bill}, \mathrm{flipper} \,|\, \mathrm{species})}{\mathbb{P}(\mathrm{bill}, \mathrm{flipper})} \\
# &= \frac{\mathbb{P}(\mathrm{species}) \mathbb{P}(\mathrm{bill} \,|\, \mathrm{species}) \mathbb{P}(\mathrm{flipper} \,|\, \mathrm{species})}{\mathbb{P}(\mathrm{bill}, \mathrm{flipper})}.
# \end{aligned}
# ```
# This predictive model is known as
# [naive Bayes classifier](https://en.wikipedia.org/wiki/Naive_Bayes_classifier).
#
# In the section above, we estimated $\mathbb{P}(\mathrm{species})$,
# $\mathbb{P}(\mathrm{bill} \,|\, \mathrm{species})$, and
# $\mathbb{P}(\mathrm{flipper} \,|\, \mathrm{species})$ for each penguin species from
# the training data. For the conditional distributions we used a Gaussian approximation.

function predict_naive_bayes_classifier(fit, data)
    ## Compute unnormalized probabilities
    z =
        log.(permutedims(fit.proportion)) .+
        logpdf.(permutedims(fit.bill), data.bill_length_mm) .+
        logpdf.(permutedims(fit.flipper), data.flipper_length_mm)

    ## Normalize probabilities
    u = maximum(z; dims=2)
    z .= exp.(z .- u)
    sum!(u, z)
    z ./= u

    return DataFrame(z, fit.species)
end

train_predict = predict_naive_bayes_classifier(penguins_fit, train_penguins)
val_predict = predict_naive_bayes_classifier(penguins_fit, val_penguins);

# ## Evaluation
#
# We evaluate the probabilistic predictions of the naive Bayes classifier that we just
# trained. It is easier to work with a numerical encoding of the true penguin species and a
# corresponding vector of predictions.

train_species = convert(Vector{Int}, indexin(train_penguins.species, names(train_predict)))
train_probs = collect(Vector{Float64}, eachrow(train_predict))

val_species = convert(Vector{Int}, indexin(val_penguins.species, names(val_predict)))
val_probs = collect(Vector{Float64}, eachrow(val_predict));

# ### Log-likelihood
#
# We compute the average log-likelihood of the training and validation data. It is
# equivalent to the negative cross-entropy.

function mean_loglikelihood(species, probs)
    return mean(log(p[s]) for (s, p) in zip(species, probs))
end

mean_loglikelihood(train_species, train_probs)

#-

mean_loglikelihood(val_species, val_probs)

# ### Brier score
#
# The average log-likelihood is also equivalent to the
# [logarithmic score](https://sites.stat.washington.edu/raftery/Research/PDF/Gneiting2007jasa.pdf).
# The Brier score is another strictly proper scoring rule that can be used for evaluating
# probabilistic predictions.

function brier_score(species, probs)
    return mean(
        sum(abs2(pi - (i == s)) for (i, pi) in enumerate(p)) for
        (s, p) in zip(species, probs)
    )
end

brier_score(train_species, train_probs)

#-

brier_score(val_species, val_probs)

# ### Expected calibration error
#
# As all proper scoring rules, the logarithmic and the Brier score can be [decomposed in
# three terms that quantify the sharpness and calibration of the predictive model and the
# irreducible uncertainty of the targets that is inherent to the prediction
# problem](https://doi.org/10.1002/qj.456). The calibration term in this decomposition is
# the expected calibration error (ECE)
# ```math
# \mathbb{E} d\big(P_X, \mathrm{law}(Y \,|\, P_X)\big)
# ```
# with respect to the score divergence $d$.
#
# Scoring rules, however, include also the sharpness and the uncertainty term. Thus models
# can trade off calibration for sharpness and therefore scoring rules are not suitable for
# specifically evaluating calibration of predictive models.
#
# The score divergence to the logarithmic and the Brier score are the Kullback-Leibler (KL)
# divergence
# ```math
# d\big(P_X, \mathrm{law}(Y \,|\, P_X)\big) = \sum_{y} \mathbb{P}(Y = y \,|\, P_X)
# \log\big(\mathbb{P}(Y = y \,|\, P_X) / P_X(\{y\})\big)
# ```
# and the squared Euclidean distance
# ```math
# d\big(P_X, \mathrm{law}(Y \,|\, P_X)\big) = \sum_{y} \big(P_X - \mathrm{law}(Y \,|\, P_X)\big)^2(\{y\}),
# ```
# respectively. The KL divergence is defined only if $\mathrm{law}(Y \,|\, P_X)$ is
# absolutely continuous with respect to $P_X$, i.e., if $P_X(\{y\}) = 0$ implies
# $\mathbb{P}(Y = y \,|\, P_X) = 0$.

# We estimate the ECE by binning the probability simplex of predictions $P_X$ and computing
# the weighted average of the distances between the mean prediction and the distribution of
# targets in each bin.
#
# One approach is to use bins of uniform size.

ece = ECE(UniformBinning(10), (μ, y) -> kl_divergence(y, μ))
ece(train_probs, train_species)

#-

ece(val_probs, val_species)

# For the squared Euclidean distance we obtain:

ece = ECE(UniformBinning(10), SqEuclidean())
ece(train_probs, train_species)

#-

ece(val_probs, val_species)

# Alternatively, one can use a data-dependent binning scheme that tries to split the
# predictions in a way that minimizes the variance in each bin.
#
# With the KL divergence we get:

ece = ECE(MedianVarianceBinning(5), (μ, y) -> kl_divergence(y, μ))
ece(train_probs, train_species)

#-

ece(val_probs, val_species)

# For the squared Euclidean distance we obtain:

ece = ECE(MedianVarianceBinning(5), SqEuclidean())
ece(train_probs, train_species)

#-

ece(val_probs, val_species)

# We see that the estimates (of the same theoretical quantity!) are highly dependent on the
# chosen binning scheme.

# ### Kernel calibration error
#
# As an alternative to the ECE, we estimate the kernel calibration error (KCE). We keep it
# simple here, and use the tensor product kernel
# ```math
# k\big((\mu, y), (\mu', y')\big) = \delta_{y,y'} \exp{\bigg(-\frac{{\|\mu - \mu'\|}_2^2}{2\nu^2} \bigg)}
# ```
# with length scale $\nu > 0$ for predictions $\mu,\mu'$ and corresponding targets $y, y'$.
# For simplicity, we estimate length scale $\nu$ with the median heuristic.

distances = pairwise(SqEuclidean(), train_probs)
λ = sqrt(median(distances[i] for i in CartesianIndices(distances) if i[1] < i[2]))
kernel = (GaussianKernel() ∘ ScaleTransform(inv(λ))) ⊗ WhiteKernel();

# We obtain the following biased estimates of the squared KCE (SKCE):

skce = BiasedSKCE(kernel)
skce(train_probs, train_species)

#-

skce(val_probs, val_species)

# Similar to the biased estimates of the ECE, the biased estimates of the SKCE are always
# non-negative. The unbiased estimates can be negative as well, in particular if the model
# is (close to being) calibrated:

skce = UnbiasedSKCE(kernel)
skce(train_probs, train_species)

#-

skce(val_probs, val_species)

# When the datasets are large, the quadratic sample complexity of the standard biased and
# unbiased estimators of the SKCE can become prohibitive. In these cases, one can resort to
# an estimator that averages estimates of non-overlapping blocks of samples. This estimator
# allows to trade off computational cost for increased variance.
#
# Here we consider the extreme case of blocks with two samples, which yields an estimator
# with linear sample complexity:

skce = BlockUnbiasedSKCE(kernel, 2)
skce(train_probs, train_species)

#-

skce(val_probs, val_species)


text:
# This script is written for calculation the Hover Result of Helicopter Rotor

@everywhere function ui2hover(ψ,rb,T_hover=T,A_hover=A)
    # this function is for solving uniform induced velocity and total velocity of
    # blade elements in rotation coordination in hover condition
    vind_hover = [0.,0.,-sqrt(T_hover/2/ρ/A_hover)] # 动量理论求解悬停诱导速度

    # ---Initialize Induced Velocity Array & Total Velocity Array---
    vindhover_r = Array{Vector}(NR,Nbe)
    vallhover_r = Array{Vector}(NR,Nbe)

    # ---Calculate Induced Velocity and Total Velocity in rotation coordination---
    for k in 1:NR
        ψk = ψ+(k-1)*2*π/NR
        for i in 1:Nbe
            vindhover_r[k,i] = systoro(vind_hover,ψk)
            vOmega_r = [0.,-Ω*rb[k,i],0.]
            vallhover_r[k,i] = vindhover_r[k,i]+vOmega_r
        end
    end

    return vindhover_r,vallhover_r,vind_hover[3]
end

@everywhere function betaang(θcp,θlon,θlat,λ,γ)
    # calculate current flap angle
    beta0 = γ/2*θcp+γ/2*λ
    betalon = θlon
    betalat = -θlat
    return beta0,betalon,betalat
end

@everywhere function binaryse(rot, θcp, dθ, T_hover=T, simθ=10.)
    if abs(rot-T_hover)<simθ
        return true, θcp, dθ
    elseif rot<T_hover
        θcp = θcp-dθ
        return false, θcp, dθ/2.
    else
        θcp = θcp+dθ
        return false, θcp, dθ/2.
    end
end

@everywhere function uitest2hover(θ7,thelat,thelon,rb)
    # this function is for running the whole hovering calculation of rotor

    ψ = 0.0
    θcp = zeros(Float64,NR)
    θlat = zeros(Float64,NR)
    θlon = zeros(Float64,NR)
    θ0 = Array{Float64}(NR,Nbe)
    for k in 1:NR
      θcp[k] = θ7
      θlat[k] = thelat
      θlon[k] = thelon
      for i in 1:Nbe
        θ0[k,i] = θcp[k]+θtw*(rb[k,i]/R-0.7)
      end
    end
    betalon = thelon
    betalat = -thelat
    index = 1
    trimind = 1
    rot = 0.
    MQ  = 0.
    dθ  = 0.2*π/180
    trimrot = zeros(Float64,2)

    # ---Initialize Induced Velocity Field---
    uitmp = ui2hover(ψ,rb)

    while true  # trim iteration
        vindhover_r = uitmp[1]
        vallhover_r = uitmp[2]
        λ           = uitmp[3]

        # ---calculate flap coefficient
        betatmp = betaang(θcp[1],thelon,thelat,λ,(ρ*8.2*(0.062)*R^4)/Iβ)
        β0      = betatmp[1]
        βlon    = betatmp[2]
        βlat    = betatmp[3]

        # ---calculate current flap angle
        betaval = β0-βlon*cos(ψ)-βlat*sin(ψ)

        # ---calculate the a.o.a of every blade element---
        aoatmp = aoaget(vallhover_r,ψ,θ0,θlat,θlon,betaval)
        α_aero = aoatmp[1]
        θ      = aoatmp[2]

        # ---calculate the lift&drag coefficient of each blade element
        clcdtmp = clcdget(α_aero,vallhover_r)
        cl = clcdtmp[1]
        cd = clcdtmp[2]

        # calculate the whole aerodynamic force and moment of the rotor in hover
        rftmp = rotorforce(ψ,vallhover_r,α_aero,θ,cl,cd)
        rot = rot+rftmp[3]
        MQ  = MQ+rftmp[4]
        ψ += dψ

        # judge if trim condition is satisfied
        if index%npsi==0
            rot = rot/npsi
            if (trimind-1)%3!=0
                trimrot[(trimind-1)%3] = rot
            else
                trimrot0 = rot
            end
            MQ  = MQ/npsi
            print("---The Trust now is : $(rot) N\n")
            if abs(rot-T)<10||index>=1000*npsi
                print("---The collective Pitch is : $(θcp/π*180)\n")
                break
            else
                print("\n$(trimind)\n")
                if trimind%3 == 1
                    for k in 1:NR
                        θcp[k] = θcp[k]+dθ
                        for i in 1:Nbe
                            θ0[k,i] = θcp[k]+θtw*(rb[k,i]/R-0.7)
                        end
                    end
                elseif trimind%3 == 2
                    for k in 1:NR
                        θcp[k] = θcp[k]-dθ
                        for i in 1:Nbe
                            θ0[k,i] = θcp[k]+θtw*(rb[k,i]/R-0.7)
                        end
                    end
                else
                    for k in 1:NR
                        θcp[k] = θcp[k]+(T-trimrot0)*2*dθ/(trimrot[1]-trimrot[2])
                        for i in 1:Nbe
                            θ0[k,i] = θcp[k]+θtw*(rb[k,i]/R-0.7)
                        end
                    end
                end
            end

            ψ = 0.  # Initialize the azimuth every circle
            trimind += 1
            rot = 0.
            MQ  = 0.
            # # ---run the binary search---
            # bsetmp = binaryse(rot, θcp[1], dθ)
            # if bsetmp[1]||index==100*npsi
            #     print("---The collective Pitch is : $(θcp/π*180)\n")
            #     break
            # else
            #     dθ = bsetmp[3]
            #     print("\n$(bsetmp)\n\n")
            #     for k in 1:NR
            #         θcp[k] = bsetmp[2]
            #         for i in 1:Nbe
            #             θ0[k,i] = θcp[k]+θtw*(rb[k,i]/R-0.7)
            #         end
            #     end
            #     rot = 0.
            #     MQ  = 0.
            # end
        end

        uitmp = ui2hover(ψ,rb)
        index += 1
    end

    if index>=1000*npsi
        return false,abs(MQ*Ω*10)
    else
        return true,abs(MQ*Ω)
    end
end

max_stars_repo_name: dingzhw/RotorTrim

max_issues_repo_name: helicopterliu/RotorTrim

max_forks_repo_name: helicopterliu/RotorTrim


text:
# Based on an example from
#	https://jump.dev/JuMP.jl/previews/PR2476/examples/Semidefinite%20programs/max_cut_sdp/
# which in turn uses methods from
#	Goemans, M. X., & Williamson, D. P. (1995). Improved approximation
#	algorithms for maximum cut and satisfiability problems using semidefinite
#	programming. Journal of the ACM (JACM), 42(6), 1115-1145.

using SparseArrays
using LinearAlgebra

# Model the max-cut SDP relaxation --------------------------------------------

cut_value(C, cut) =	sum(C .* (cut * cut'))
laplacian(W) = spdiagm(vec(sum(W, dims=1))) - W

function extract_cut(X)
	rnd = normalize(rand(size(X, 1)))

	factorization = cholesky(Hermitian(X), Val(true); check=false)
	V = (factorization.P * factorization.L)'
	vs = map(normalize, eachcol(V))

	map(x -> dot(rnd, x) |> sign, vs)
end

function max_cut_sdp_model(weights)
	n = size(weights, 1)

	C = laplacian(weights) / 4
	As = [sparse([i], [i], [1], n, n) for i in 1:n]
	b = ones(n) / n

	C, As, b
end

function solve_maxcut(W)
	n = size(W, 1)

	C, As, b = max_cut_sdp_model(W)
	X = sketchy_cgal(normalize(-C), As, b; R=n, info_io=devnull)

	cut = extract_cut(X)
	value = cut_value(C, cut)

	cut, value
end

# Examples --------------------------------------------------------------------

function maxcut_example_1()
	W = [
		0 5;
		5 0;
	]

	cut, value = solve_maxcut(W)

	@test isapprox(value, 5)
	@test cut[1] != cut[2]
end

function maxcut_example_2()
	W = [
		0 5 7 6;
		5 0 0 1;
		7 0 0 1;
		6 1 1 0;
	]

	cut, value = solve_maxcut(W)

	@test isapprox(value, 18)
	@test cut[1] != cut[2]
	@test cut[2] == cut[3] == cut[4]
end

function maxcut_example_3()
	W = [
		0 1 5 0;
		1 0 0 9;
		5 0 0 2;
		0 9 2 0;
	]

	cut, value = solve_maxcut(W)

	@test isapprox(value, 17)
	@test cut[1] == cut[4]
	@test cut[2] == cut[3]
	@test cut[1] != cut[2]
end


text:
immutable Inplace{N} end

immutable Transpose{T}
    obj::T
end
immutable CTranspose{T}
    obj::T
end

immutable Operator{S} end

cat_tuple(t1::Tuple, t2::Tuple) = tuple(t1...,t2...)


text:
import LightGraphs as LG
import Random
import Base.Threads
include(joinpath(pwd(), "src", "Heuristics.jl"))

mutable struct GraphData
    g::LG.AbstractGraph
    e_width::Int
    ub::Int
    best_order::Vector{Int}
    nums::Vector{Int}
    finishTime::Float64 # epoch
    numNodes::Int # num nodes visited
    timeout::Bool
    seed::Int
    randomContext::Random.MersenneTwister
end

function branch_bound(g::LG.AbstractGraph, duration::Float64, seed=time()::Int)
    G = deepcopy(g)
    # best_order, ub = min_width(G)
    best_order = Int[]
    ub = 10000
    x = Int[]
    println(best_order, ub)
    nums = genNumberList(LG.nv(g))
    finishTime = time() + duration
    gd = GraphData(G, 0, ub, best_order, nums, finishTime, 1, false, seed, Random.MersenneTwister(seed))
    dfs(gd, x)
    gd
end

@inline function process(gd::GraphData, x::Vector{Int}, v::Int, nVertices::Int)

    # Connecting edges
    gm = joinVerts!(gd.g, v)
    temp_ewidth = gd.e_width
    gd.e_width = max(gd.e_width, LG.degree(gd.g, v))
    LG.rem_vertex!(gd.g, v)

    # Handling nums
    temp = gd.nums[v]
    gd.nums[v] = gd.nums[nVertices]

    if gd.e_width < gd.ub
        # gd_new = GraphData(gd.g, gd.e_width, gd.ub, gd.best_order, nums_old, gd.finishTime, 1, gd.timeout)
        dfs(gd, [x; temp])
    end

    gd.e_width = temp_ewidth
    unmodify_graph!(gd.g, gm, v)
    gd.nums[v] = temp
end

function dfs(gd::GraphData, x::Vector{Int})
    gd.numNodes += 1
    nVertices = LG.nv(gd.g)

    if time() >= gd.finishTime
        # println("timeout!")
        gd.timeout = true
        return
    end

    if nVertices < gd.e_width
        if gd.e_width < gd.ub
            gd.ub = gd.e_width
            gd.best_order = [x; gd.nums[1:nVertices]]
        end
        return
    end

    # If parallel
    if length(x) < 1
        gds = [deepcopy(gd) for _ = 1:Threads.nthreads()]
        for v in gds
            v.randomContext = Random.MersenneTwister(v.seed + Threads.threadid())
        end

        Threads.@threads for v in Random.shuffle(gds[1].randomContext, LG.vertices(gd.g))
            if gd.timeout break end
            process(gds[Threads.threadid()], x, v, nVertices)
        end

        # Pick out best order here
        best_gd = argmin([_gd.ub for _gd in gds])
        println("best ub: ", gds[best_gd].ub)
        println("Did timeout: ", gds[best_gd].timeout)
        gd.numNodes = sum([_gd.numNodes for _gd in gds])
        gd.best_order = gds[best_gd].best_order
        gd.ub = gds[best_gd].ub
        gd.timeout = gds[best_gd].timeout

    else
        for v in Random.shuffle(gd.randomContext, LG.vertices(gd.g))
            if gd.timeout break end
            process(gd, x, v, nVertices)
        end
    end
end

# revert changes made to the graph
@inline function unmodify_graph!(g::LG.AbstractGraph, gm::GraphModification, v::Int)
    LG.add_vertex!(g)

    nb = copy(LG.neighbors(g, v))

    for i in nb
        LG.add_edge!(g, i, LG.nv(g))
        LG.rem_edge!(g, i, v)
    end

    for i in gm.nb
        LG.add_edge!(g, i, v)
    end

    for i in gm.addedEdges
        LG.rem_edge!(g, i[1], i[2])
    end
end

max_stars_repo_name: OleeL/treewidth-heuristics

max_issues_repo_name: OleeL/treewidth-heuristics

max_forks_repo_name: OleeL/treewidth-heuristics


text:
using Redpitaya
using PyPlot
using FFTW

rp = RedPitaya("192.168.178.66")

freq = 25000.0
dec = optimalDecimation(rp,freq)
#freq=25000
numPeriods = 5
freqR = roundFreq(rp,dec,freq)
numSampPerPeriod = numSamplesPerPeriod(rp,dec,freqR)
numSamp = numSampPerPeriod*numPeriods
if numSamp>16384
#	error("To many sampling points for RedPitaya buffer. Reduce Decimation or reduce number of periods")
end

println("Frequency = $freqR Hz")
println("Number Sampling Points per Period: $numSampPerPeriod")

# start sending
send(rp,"GEN:RST")
#send(rp,"SOUR1:BURS:STAT ON")
#change this to SQUARE to get a SQUARE excitation
sendAnalogBurstSignal(rp,1,"SINE",freqR,0.4,offset=0.2,ncycBurst=numPeriods+1,repetitions=2,repetitionCycleTime=1/freq*10*1e6)
#send(rp,"SOUR1:VOLT:OFFS 0")
# receive data
trigger = "AWG_PE" # "CH1_PE"
u1 = receiveAnalogSignalWithTrigger(rp, 1, 0, numSamp, dec=dec, delay=0.2, typ="OLD",
                      trigger=trigger, triggerLevel=-0.0, binary=true, raw=false,
                      triggerDelay=numSampPerPeriod)

figure(3)
subplot(2,1,1)
plot(u1)
subplot(2,1,2)
semilogy(abs.(rfft(u1)),"o-b",lw=2)

#Lets check if the data is periodic
figure(4)
plot(1:length(u1),u1,"x-r",lw=2)
plot(length(u1)+1:2*length(u1),u1,"x-b",lw=2)
send(rp,"GEN:RST")

trapez=trapeziod(0.9)
sendAnalogBurstSignal(rp,1,"ARBITRARY",freqR,0.4,offset=0.2,ncycBurst=5+1,repetitions=2,repetitionCycleTime=1/freq*10*1e6,AWG_data=trapez)
#send(rp,"SOUR1:VOLT:OFFS 0")
# receive data
trigger = "AWG_PE" # "CH1_PE"
u1 = receiveAnalogSignalWithTrigger(rp, 1, 0, numSamp*4, dec=dec, delay=0.0, typ="OLD",
                      trigger=trigger, triggerLevel=-0.1, binary=true, raw=false,
                      triggerDelay=0*numSampPerPeriod)

figure(5)
subplot(2,1,1)
plot(u1)


text:
import Distributions.MvNormalStats, Distributions.lpgamma, Distributions.suffstats, Distributions.mean

function suffstats{Tf <: FloatingPoint}(D::Type{MvNormal}, x::DenseArray{Tf, 2},
        w::DenseArray{Float64}, z::DenseArray{Tf, 2})
    d = size(x, 1)
    n = size(x, 2)

    tw = sum(w)
    s = zeros(d)
    for i=1:n
        @devec s[:] += x[:, i] .* w[i]
    end
    m = s * inv(tw)
    for i in 1:n
        z_i = view(z, :, i)
        x_i = view(x, :, i)
        @devec z_i[:] = (x_i .- m) .* sqrt(w[i])
    end
    s2 = A_mul_Bt(z, z)

    MvNormalStats(s, m, s2, tw)
end

function posterior_cool(prior::NormalWishart, ss::MvNormalStats)
    if ss.tw < eps() return prior end

    mu0 = prior.mu
    kappa0 = prior.kappa
    TC0 = prior.Tchol
    nu0 = prior.nu

    kappa = kappa0 + ss.tw
    nu = nu0 + ss.tw
    mu = (kappa0.*mu0 + ss.s) ./ kappa

    Lam0 = inv(TC0)
    z = prior.zeromean ? ss.m : ss.m - mu0
    Lam = inv(Lam0 + ss.s2 + kappa0*ss.tw/kappa*(z*z'))

    return NormalWishart(mu, kappa, cholfact(Lam), nu)
end

function expected_logdet(nw::NormalWishart)
    logd = 0

    for i=1:nw.dim
        logd += digamma(0.5 * (nw.nu + 1 - i))
    end

    logd += nw.dim * log(2)
    logd += logdet(nw.Tchol)

    return logd
end

function expected_truncation_level(nw::NormalWishart)
    return nw.Tchol[:U]' * nw.Tchol[:U] * nw.nu
end

function mean(nw::NormalWishart)
    T = expected_truncation_level(nw)
    return MvNormalCanon(T * nw.mu, T)
end

function lognorm(nw::NormalWishart)
    return (nw.dim / 2) * (log(2 * pi) - log(nw.kappa)) + (nw.nu / 2) * logdet(nw.Tchol) + (nw.dim * nw.nu / 2) * log(2.) + lpgamma(nw.dim, nw.nu / 2)
end

function entropy(nw::NormalWishart)
    en = 0.

    logd = expected_logdet(nw)

    en -= logd
    en -= nw.dim * (log(nw.kappa) - log(2*pi))
    en += nw.dim
    en /= 2

    en -= 0.5 * (nw.nu - nw.dim - 1) * logd
    en += nw.nu * nw.dim / 2
    en += 0.5 * nw.nu * logdet(nw.Tchol) + 0.5 * nw.nu * nw.dim * log(2.) + lpgamma(nw.dim, nw.nu/2)

    return en
end

function marginal_loglikelihood(prior::NormalWishart, posterior::NormalWishart, n::Float64)
    d = prior.dim
    return lognorm(posterior) - lognorm(prior) - (n*d/2) * log(2 * pi)
end

function gaussian_mixture{Tf <: FloatingPoint}(prior::NormalWishart, T::Int64, alpha::Float64, x::DenseArray{Tf, 2})
    dim, N = size(x)
    theta = Array(NormalWishart, T)

    buffer = similar(x)

    for k=1:T
        theta[k] = prior
    end

    function cluster_update(k::Int64, z::DenseArray{Float64})
        nw = posterior_cool(prior, suffstats(MvNormal, x, z, buffer))
        theta[k] = nw
    end

    function cluster_entropy(k::Int64)
        return entropy(theta[k])
    end

    function cluster_loglikelihood(k::Int64, z::DenseArray{Float64})
        post = posterior_cool(prior, suffstats(MvNormal, x, z, buffer))
        n = sum(z)

        ll = -lognorm(prior)

        nw = theta[k]
        ll -= nw.dim * n / 2 * log(2*pi)

        ll += (post.nu - post.dim)/2 * expected_logdet(nw)
        W = nw.Tchol[:U]' * nw.Tchol[:U]
        ll -= nw.nu * trace(inv(post.Tchol) * W) / 2

        dm = nw.Tchol[:U] * (nw.mu - post.mu)
        ll -= post.kappa * (nw.dim / nw.kappa + nw.nu * dot(dm, dm)) / 2

        return ll
    end

    function object_loglikelihood(k::Int64, i::Int64)
        nw = theta[k]
        ll = -nw.dim/2 * log(2*pi) + expected_logdet(nw)/2
        dx = nw.Tchol[:U] * (x[:, i] - nw.mu)
        ll -= (nw.dim / nw.kappa + nw.nu * dot(dx, dx)) / 2

        return ll
    end

    mm = TSBPMM(N, T, alpha, cluster_update,
            cluster_loglikelihood,
            object_loglikelihood,
            cluster_entropy; random_init=true)

    function predictive_loglikelihood{Tf <: FloatingPoint}(xt::DenseArray{Tf, 2})
        nt = size(xt, 2)
        logp = zeros(nt)

        π = zeros(T)
        pi!(π, mm)

        pred = map(predictive, theta)

        for i=1:nt
            for k=1:T
                logp[i] += π[k] * pdf(pred[k], x[:, i])
            end
        end
        log!(logp)

        return logp
    end

    return mm, theta, predictive_loglikelihood
end

export gaussian_mixture


text:
## random number generation

using Random
import RandomNumbers

# global state

# we provide 128 bytes of random state, shared across the thread block using shared memory.
# the same amount of seeding data is available in constant memory, set during initial load.

@eval @inline function global_random_seed()
    ptr = Base.llvmcall(
        $("""@global_random_seed = weak addrspace($(AS.Constant)) externally_initialized global [32 x i32] zeroinitializer, align 32
             define i8 addrspace($(AS.Constant))* @entry() #0 {
                 %ptr = getelementptr inbounds [32 x i32], [32 x i32] addrspace($(AS.Constant))* @global_random_seed, i64 0, i64 0
                 %untyped_ptr = bitcast i32 addrspace($(AS.Constant))* %ptr to i8 addrspace($(AS.Constant))*
                 ret i8 addrspace($(AS.Constant))* %untyped_ptr
             }
             attributes #0 = { alwaysinline }
          """, "entry"), LLVMPtr{UInt32, AS.Constant}, Tuple{})
    CuDeviceArray((32,), ptr)
end

@eval @inline function global_random_state()
    ptr = Base.llvmcall(
        $("""@global_random_state = weak addrspace($(AS.Shared)) global [32 x i32] zeroinitializer, align 32
             define i8 addrspace($(AS.Shared))* @entry() #0 {
                 %ptr = getelementptr inbounds [32 x i32], [32 x i32] addrspace($(AS.Shared))* @global_random_state, i64 0, i64 0
                 %untyped_ptr = bitcast i32 addrspace($(AS.Shared))* %ptr to i8 addrspace($(AS.Shared))*
                 ret i8 addrspace($(AS.Shared))* %untyped_ptr
             }
             attributes #0 = { alwaysinline }
          """, "entry"), LLVMPtr{UInt32, AS.Shared}, Tuple{})
    CuDeviceArray((32,), ptr)
end

function initialize_random_seeds!(mod)
    seed = CuGlobal{NTuple{32,UInt32}}(mod, "global_random_seed")
    seed[async=true] = Tuple(rand(UInt32, 32))
end

@device_override Random.make_seed() = clock(UInt32)

@generated function make_full_seed(seed::Integer)
    quote
        s_0 = seed % UInt32
        Base.@nexprs 32 i->s_i = xorshift(s_{i-1})
        Base.@ncall 32 tuple s
    end
end

# helpers

function xorshift(x::UInt32)::UInt32
    x = xor(x, x << 13)
    x = xor(x, x >> 17)
    x = xor(x, x << 5)
    return x
end

# generators

"""
    SharedTauswortheGenerator()

A maximally equidistributed combined Tausworthe generator.

The generator uses 32 bytes of random state stored in shared memory. This memory is
zero-initialized; When the first random number is generated, each block will derive an
initial state from the seed provided during module compilation, and the block identifier.
Each block will then indepentendly, but deterministically use that state to generate random
numbers. Within a block, the 32-bytes of random state are identical, so the warp identifier
will be mixed in to ensure uniqueness across warps. Finally, the first warp of each block
updates the shared state.

!!! warning

    Although the numbers obtained from this generator "look OK", they do not pass the
    SmallCrush test suite, so the generator should be deemed experimental.
"""
struct SharedTauswortheGenerator <: RandomNumbers.AbstractRNG{UInt32}
end

@inline function Base.getproperty(rng::SharedTauswortheGenerator, field::Symbol)
    if field === :seed
        global_random_seed()
    elseif field === :state
        global_random_state()
    end
end

@device_override Random.default_rng() = SharedTauswortheGenerator()

function Random.seed!(rng::SharedTauswortheGenerator, seed::Integer)
    Random.seed!(rng, make_full_seed(seed))
    return
end

"""
    Random.seed!(rng::SharedTauswortheGenerator, seeds)

Seed the on-device Tausworthe generator with a 32-element tuple or array of UInt32 seeds.
This function only needs to be called by a single warp (i.e. 32 threads) from each block,
but it does not hurt to call it from all threads (as long as threads are synchronized after
that). It should always be called by at least 32 threads to ensure the random state is fully
initialized, even if you will be using the generator from fewer threads!
"""
@inline Base.@propagate_inbounds function Random.seed!(rng::SharedTauswortheGenerator, seed)
    state = initial_state(seed)
    @inbounds rng.state[laneid()] = state
    return
end

@inline Base.@propagate_inbounds function initial_state(seeds)
    z = seeds[laneid()]

    # mix-in the block id to ensure unique values across blocks
    # XXX: is this OK? shouldn't we use a generator that allows skipping ahead?
    #      https://stackoverflow.com/questions/11692785/efficient-xorshift-skip-ahead
    blockId = blockIdx().x + (blockIdx().y - 1) * gridDim().x +
                             (blockIdx().z - 1) * gridDim().x * gridDim().y
    z = xorshift(z ⊻ blockId%UInt32)

    return z
end

# NOTE: @propagate_inbounds for when we know the passed seed contains at least 32 elements
# XXX: without @inline @propagate_inbounds, we get local memory accesses (e.g. with `rand!`)

#const TausShift1 = CuConstantMemory{UInt32}((6, 2, 13, 3))
#const TausShift2 = CuConstantMemory{UInt32}((13, 27, 21, 12))
#const TausShift3 = CuConstantMemory{UInt32}((18, 2, 7, 13))
#const TausOffset = CuConstantMemory{UInt32}((4294967294, 4294967288, 4294967280, 4294967168))
# XXX: constant memory isn't supported yet, so let's resort to llvmcall
for (name, vals) in ("TausShift1" => (6, 2, 13, 3),
                     "TausShift2" => (13, 27, 21, 12),
                     "TausShift3" => (18, 2, 7, 13),
                     "TausOffset" => (4294967294, 4294967288, 4294967280, 4294967168))
    @eval @inline function $(Symbol(name))()
        ptr = Base.llvmcall(
            $("""@$(name) = weak addrspace($(AS.Constant)) global [4 x i32] [i32 $(vals[1]), i32 $(vals[2]), i32 $(vals[3]), i32 $(vals[4])], align 4
                define i8 addrspace($(AS.Constant))* @entry() #0 {
                    %ptr = getelementptr inbounds [4 x i32], [4 x i32] addrspace($(AS.Constant))* @$(name), i64 0, i64 0
                    %untyped_ptr = bitcast i32 addrspace($(AS.Constant))* %ptr to i8 addrspace($(AS.Constant))*
                    ret i8 addrspace($(AS.Constant))* %untyped_ptr
                }
                attributes #0 = { alwaysinline }
            """, "entry"), LLVMPtr{UInt32, AS.Constant}, Tuple{})
        CuDeviceArray((4,), ptr)
    end
end

function TausStep(z::Unsigned, S1::Integer, S2::Integer, S3::Integer, M::Unsigned)
    b = (((z << S1) ⊻ z) >> S2)
    return (((z & M) << S3) ⊻ b)
end

"""
    Random.rand(rng::SharedTauswortheGenerator, UInt32)

Generate a byte of random data using the on-device Tausworthe generator.

!!! warning

    This functions performs synchronization, so should only be called uniformly (i.e. by all
    threads in a block). Do not call it from a branch that diverges within a block, or your
    kernel may deadlock.
"""
function Random.rand(rng::SharedTauswortheGenerator, ::Type{UInt32})
    @inline pow2_mod1(x, y) = (x-1)&(y-1) + 1

    threadId = UInt32(threadIdx().x + (threadIdx().y - 1) * blockDim().x +
                                      (threadIdx().z - 1) * blockDim().x * blockDim().y)
    warpId = (threadId-UInt32(1)) >> 5 + UInt32(1)  # fld1
    i = pow2_mod1(threadId, 32)
    j = pow2_mod1(threadId, 4)

    @inbounds begin
        # get state
        z = rng.state[i]
        if z == 0
            z = initial_state(rng.seed)
        end

        # mix-in the warp id to ensure unique values across blocks.
        # we have max 1024 threads per block, so can safely shift by 16 bits.
        # XXX: see comment in `initial_state`
        z = xorshift(z ⊻ (warpId << 16))

        sync_threads()

        # advance & update state
        S1, S2, S3, M = TausShift1()[j], TausShift2()[j], TausShift3()[j], TausOffset()[j]
        state = TausStep(z, S1, S2, S3, M)
        if warpId == 1
            rng.state[i] = state
        end

        sync_threads()

        # generate based on 4 bytes of state
        mask = active_mask()
        i1 = pow2_mod1(threadId+1, 32)
        i2 = pow2_mod1(threadId+2, 32)
        i3 = pow2_mod1(threadId+3, 32)
        if active_mask() == typemax(UInt32)
            # we have a full warp, so we can safely shuffle.
            # get the warp-local states from neighbouring threads.
            s1 = shfl_sync(FULL_MASK, state, i1)
            s2 = shfl_sync(FULL_MASK, state, i2)
            s3 = shfl_sync(FULL_MASK, state, i3)
        else
            # can't shuffle, so fall back to fetching block-global state.
            # this results is less entropy, as it reuses data from the first warp.
            s1 = rng.state[i1]
            s2 = rng.state[i2]
            s3 = rng.state[i3]
        end
        state ⊻ s1 ⊻ s2 ⊻ s3
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess

@testset "Tet10 + convection" begin
    # For some reason Tet10 fails, maybe because of convection.
    mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
    mesh = aster_read_mesh(mesh_file, "TETRA_TET10_1")
    prob = Problem(Heat, "tet", 1)
    face = Problem(Heat, "face 4", 1)
    fixed = Problem(Dirichlet, "fixed face 3", 1, "temperature")
    prob.elements = create_elements(mesh, "TET")
    update!(prob, "temperature thermal conductivity", 50.0)
    face.elements = create_elements(mesh, "FACE4")
    update!(face, "temperature external temperature", 20.0)
    update!(face, "temperature heat transfer coefficient", 60.0)
    fixed.elements = create_elements(mesh, "FACE2")
    info("# of elements in fixed set: $(length(fixed))")
    update!(fixed, "temperature 1", 0.0)
    solver = LinearSolver(prob, face, fixed)
    solver()
    T = prob.assembly.u
    info("Solution: $T")
    T_expected = [ # using code aster
        1.45606533688540E+01
        5.01315339269860E-17
        3.02236827927507E-17
       -2.01049663215778E-16
        1.05228712963739E+01
        0.00000000000000E+00
        9.44202309239159E+00
        1.05228712963739E+01
        4.44089209850063E-16
        0.00000000000000E+00]
    @test isapprox(T, T_expected; rtol=1.0e-6)
end

@testset "2d heat problem (one element)" begin

    X = Dict{Int, Vector{Float64}}(
        1 => [0.0,0.0],
        2 => [1.0,0.0],
        3 => [1.0,1.0],
        4 => [0.0,1.0])

    # define volume element
    el1 = Element(Quad4, [1, 2, 3, 4])

    update!(el1, "geometry", X)
    update!(el1, "temperature thermal conductivity", 6.0)
    update!(el1, "temperature load", 12.0)

    # define boundary element for flux
    el2 = Element(Seg2, [1, 2])
    update!(el2, "geometry", X)
    # linear ramp from 0 -> 6 in time 0 -> 1
    update!(el2, "temperature flux", 0.0 => 0.0)
    update!(el2, "temperature flux", 1.0 => 6.0)

    # define heat problem and push elements to problem
    problem = Problem(Heat, "one element heat problem", 1)
    problem.properties.formulation = "2D"
    push!(problem, el1, el2)

    # Set constant source f=12 with k=6. Accurate solution is
    # T=1 on free boundary, u(x,y) = -1/6*(1/2*f*x^2 - f*x)
    # when boundary flux not active (at t=0)
    assemble!(problem, 0.0)
    A = full(problem.assembly.K)
    b = full(problem.assembly.f)
    A_expected = [
         4.0 -1.0 -2.0 -1.0
        -1.0  4.0 -1.0 -2.0
        -2.0 -1.0  4.0 -1.0
        -1.0 -2.0 -1.0  4.0]
    free_dofs = [1, 2]
    @test isapprox(A, A_expected)
    @test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [1.0, 1.0])

    # Set constant flux g=6 on boundary. Accurate solution is
    # u(x,y) = x which equals T=1 on boundary.
    # at time t=1.0 all loads should be on.
    empty!(problem)
    assemble!(problem, 1.0)
    A = full(problem.assembly.K)
    b = full(problem.assembly.f)
    @test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [2.0, 2.0])
end

function T_acc(x)
    # accurate solution
    a = 0.01
    L = 0.20
    k = 50.0
    Tᵤ = 20.0
    h = 10.0
    P = 4*a
    A = a^2
    α = h
    β = sqrt((h*P)/(k*A))
    T̂ = 100.0
    C = [1.0 1.0; (α+k*β)*exp(β*L) (α-k*β)*exp(-β*L)] \ [T̂-Tᵤ, 0.0]
    return dot(C, [exp(β*x), exp(-β*x)]) + Tᵤ
end

#=
@testset "test 1d heat problem" begin
    X = Dict{Int, Vector{Float64}}(
        1 => [0.0, 0.0, 0.0],
        2 => [0.1, 0.0, 0.0],
        3 => [0.2, 0.0, 0.0])
    e1 = Element(Seg2, [1, 2])
    e2 = Element(Seg2, [2, 3])
    e3 = Element(Poi1, [3])

    p1 = Problem(Heat, "1d heat problem", 1)
    p1.properties.formulation = "1D"
    push!(p1, e1, e2, e3)
    update!(p1, "geometry", X)
    a = 0.010
    update!(p1, "cross-section area", a^2)
    update!(p1, "cross-section perimeter", 4*a)
    update!(p1, "temperature thermal conductivity", 50.0) # k [W/(m∘C)]
    update!(p1, "temperature heat transfer coefficient", 10.0) # h [W/(m²∘C)]
    update!(p1, "temperature external temperature", 20.0)

    p2 = Problem(Dirichlet, "left boundary", 1, "temperature")
    e3 = Element(Poi1, [1])
    update!(e3, "geometry", X)
    update!(e3, "temperature 1", 100.0)
    push!(p2, e3)

    solver = LinearSolver(p1, p2)
    solver()
    T_min = minimum(p1.assembly.u)
    @test isapprox(T_max, T_acc(0.2); rtol=4.5e-2)
end
=#

@testset "compare simple 3d heat problem to code aster solution" begin
    fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med"
    mesh = aster_read_mesh(fn, "Hex8")
    element_sets = join(keys(mesh.element_sets), ", ")
    info("element sets: $element_sets")

    p1 = Problem(Heat, "rod", 1)
    rod = create_elements(mesh, "ROD")
    face2 = create_elements(mesh, "FACE2")
    face3 = create_elements(mesh, "FACE3")
    face4 = create_elements(mesh, "FACE4")
    face5 = create_elements(mesh, "FACE5")
    face6 = create_elements(mesh, "FACE6")
    update!(rod, "temperature thermal conductivity", 50.0)
    update!(face2, "temperature external temperature", 20.0)
    update!(face2, "temperature heat transfer coefficient", 60.0)
    update!(face3, "temperature external temperature", 30.0)
    update!(face3, "temperature heat transfer coefficient", 50.0)
    update!(face4, "temperature external temperature", 40.0)
    update!(face4, "temperature heat transfer coefficient", 40.0)
    update!(face5, "temperature external temperature", 50.0)
    update!(face5, "temperature heat transfer coefficient", 30.0)
    update!(face6, "temperature external temperature", 60.0)
    update!(face6, "temperature heat transfer coefficient", 20.0)
    push!(p1, rod, face2, face3, face4, face5, face6)

    p2 = Problem(Dirichlet, "left support T=100", 1, "temperature")
    push!(p2, create_elements(mesh, "FACE1"))
    update!(p2, "temperature 1", 100.0)

    solver = LinearSolver(p1, p2)
    solver()

    # fields extracted from Code Aster .resu file
    TEMP = Dict{Int64, Float64}(
        1 => 1.00000000000000E+02,
        2 => 1.00000000000000E+02,
        3 => 1.00000000000000E+02,
        4 => 1.00000000000000E+02,
        5 => 3.01613322896279E+01,
        6 => 3.01263406641066E+01,
        7 => 3.02559777927923E+01,
        8 => 3.02209215997131E+01)
    FLUX_ELGA = Dict{Int64, Vector{Float64}}(
        1 => [1.74565160615448E+04, -9.99903237329079E+01, -3.69874201221677E+01],
        2 => [1.74565160615448E+04, -3.73168968436642E+02, -1.38038931136833E+02],
        3 => [1.74428571293096E+04, -9.99903237329079E+01, -3.70268090662933E+01],
        4 => [1.74428571293096E+04, -3.73168968436642E+02, -1.38185932677561E+02],
        5 => [1.74615686370955E+04, -9.99509347888079E+01, -3.69874201221677E+01],
        6 => [1.74615686370955E+04, -3.73021966895897E+02, -1.38038931136833E+02],
        7 => [1.74479150854902E+04, -9.99509347888065E+01, -3.70268090662933E+01],
        8 => [1.74479150854901E+04, -3.73021966895874E+02, -1.38185932677561E+02])
    FLUX_NOEU = Dict{Int64, Vector{Float64}}(
        1 => [1.74596669275930E+04,  7.55555618070503E-11,  3.68594044175552E-12],
        2 => [1.74684148339734E+04,  1.10418341137120E-11,  3.48876483258209E-12],
        3 => [1.74360055518019E+04,  7.91828824731056E-11,  1.95399252334028E-13],
        4 => [1.74447696000717E+04, -3.49587025993969E-12,  3.55271367880050E-13],
        5 => [1.74596669275931E+04, -4.73227515822099E+02, -1.74958127606525E+02],
        6 => [1.74684148339733E+04, -4.72904678032251E+02, -1.74958127606524E+02],
        7 => [1.74360055518019E+04, -4.73227515822118E+02, -1.75280965396335E+02],
        8 => [1.74447696000717E+04, -4.72904678032179E+02, -1.75280965396335E+02])

    postprocessor = Postprocessor(p1)
    flux = full(postprocessor())
    fluxd = Dict{Int64, Vector{Float64}}()
    for j=1:8
        fluxd[j] = vec(flux[j,:])
    end

    T = p1("temperature")

    for j in sort(collect(keys(T)))
        T1 = T[j][1]
        T2 = TEMP[j]
        rtol = norm(T1-T2)/max(T1,T2)*100.0
        @printf "node %i temp, JF: %e, CA: %e, rtol: %10.6f %%\n" j T1 T2 rtol
        @test rtol < 1.0e-9
    end

    for j=1:8
        q1 = get_integration_points(first(rod))[j]("heat flux", 0.0)
        q2 = FLUX_ELGA[j]
        rtol = norm(q1-q2)/max(norm(q1),norm(q2))*100.0
        @printf "ip   %i flux, JF: (% e,% e,% e), CA: (% e,% e,% e), rtol: %10.6f %%\n" j q1... q2... rtol
        # @test rtol < 0.05
        # testing in integration points makes no sense because they are in different order in CA
    end

    for j in sort(collect(keys(fluxd)))
        q1 = fluxd[j]
        q2 = FLUX_NOEU[j]
        rtol = norm(q1-q2)/max(norm(q1),norm(q2))*100.0
        @printf "node %i flux, JF: (% e,% e,% e), CA: (% e,% e,% e), rtol: %10.6f %%\n" j q1... q2... rtol
        @test rtol < 1.0e-9
    end

end

@testset "compare simple 3d heat problem to analytical solution" begin
    function calc_3d_heat_model(mesh_name)
        fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med"
        mesh = aster_read_mesh(fn, mesh_name)
        p1 = Problem(Heat, "rod", 1)
        p2 = Problem(Dirichlet, "left support T=100", 1, "temperature")
        p1.elements = create_elements(mesh, "ROD", "FACE2")
        p2.elements = create_elements(mesh, "FACE1")
        update!(p1, "temperature thermal conductivity", 100.0)
        update!(p1, "temperature external temperature", 0.0)
        update!(p1, "temperature heat transfer coefficient", 1000.0)
        update!(p2, "temperature 1", 100.0)
        solver = LinearSolver(p1, p2)
        solver()
        T_min = minimum(p1.assembly.u)
        return T_min
    end
    for model in ["Tet4", "Tet10", "Hex8", "Hex20", "Hex27"]
        Tmin = calc_3d_heat_model(model)
        Tacc = 100/3
        rtol = norm(Tmin-Tacc)/max(Tmin,Tacc)*100.0
        @printf "%-10s : Tmin = % g, Tacc = % g, rtol = %g %%\n" model Tmin Tacc rtol
        @test isapprox(Tmin, 100/3)
    end
end

@testset "compare simple 3d heat problem to code aster solution" begin

    function calc_3d_heat_model(mesh_name)
        fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med"
        mesh = aster_read_mesh(fn, mesh_name)
        element_sets = join(keys(mesh.element_sets), ", ")
        info("element sets: $element_sets")
        # x -> FACE1 ... FACE2
        # y -> FACE3 ... FACE4
        # z -> FACE5 ... FACE6
        # rod has longer dimension in x direction, first face comes
        # first in corresponding axis direction
        p1 = Problem(Heat, "rod", 1)
        rod = create_elements(mesh, "ROD")
        face2 = create_elements(mesh, "FACE2")
        face3 = create_elements(mesh, "FACE3")
        face4 = create_elements(mesh, "FACE4")
        face5 = create_elements(mesh, "FACE5")
        face6 = create_elements(mesh, "FACE6")
        update!(rod, "temperature thermal conductivity", 50.0)
        update!(face2, "temperature external temperature", 20.0)
        update!(face2, "temperature heat transfer coefficient", 60.0)
        update!(face3, "temperature external temperature", 30.0)
        update!(face3, "temperature heat transfer coefficient", 50.0)
        update!(face4, "temperature external temperature", 40.0)
        update!(face4, "temperature heat transfer coefficient", 40.0)
        update!(face5, "temperature external temperature", 50.0)
        update!(face5, "temperature heat transfer coefficient", 30.0)
        update!(face6, "temperature external temperature", 60.0)
        update!(face6, "temperature heat transfer coefficient", 20.0)
        push!(p1, rod, face2, face3, face4, face5, face6)
        p2 = Problem(Dirichlet, "left support T=100", 1, "temperature")
        p2.elements = create_elements(mesh, "FACE1")
        update!(p2, "temperature 1", 100.0)
        solver = LinearSolver(p1, p2)
        solver()
        return p1.assembly.u
    end

    CA_sol = Dict(
        "Tet4" => 3.01872246268290E+01,
        "Hex8" => 3.01263406641066E+01,
        "Tet10" => 4.38924023356612E+01,
        "Hex20" => 4.57539800177123E+01,
        "Hex27" => 4.57760386068096E+01)

    models = ["Tet4", "Hex8", "Hex20", "Hex27", "Tet10"]

    for model in models
        T = calc_3d_heat_model(model)
        T_min = minimum(T)
        T_ca = CA_sol[model]
        rtol = norm(T_min-T_ca)/max(T_min,T_ca)*100.0
        @printf "%-10s : T_min = % g, T_ca = % g, rtol = %g %%\n" model T_min T_ca rtol
        if rtol > 1.0e-9
            info("Solution vector")
            dump(T)
        end
        @test rtol < 1.0e-9
    end

end