text:
using Dynesty
using Test
using Distributions

@testset "Dynesty.jl" begin
    # Write your tests here.
    # define a distribution
    d = MvNormal(ones(5))

    loglikelihood(x) = logpdf(d, x)

    prior_transform(p) = -5.0 .+ 10.0.*p

    smplr = NestedSampler(5; nlive=1000)
    dysmplr = DynamicNestedSampler(5)
    # sample using dynamic nested sampling with 500 initial live points
    res1 = sample(loglikelihood, prior_transform, smplr; dlogz=0.01, print_progress=false)
    res2 = sample(loglikelihood, prior_transform, smplr; dlogz=0.01, print_progress=false)
    res3 = sample(loglikelihood, prior_transform, smplr; dlogz=0.01, print_progress=false)

    mres = Dynesty.merge(res1, res2, res3)

    res = sample(loglikelihood, prior_transform, dysmplr; nlive_init=200, print_progress=false)

    runplot(res1)
    cornerplot(res1)
    traceplot(res1)
    boundplot(res1, (1,2),idx=10)
    cornerbound(res1, it=100)
    cornerpoints(res1)

    esamples = resample_equal(res3, Int(ceil(2*Dynesty.ess(res1))))

    @test isapprox(mean(getindex.(esamples, 1)), mean(d)[1], atol=1e-1)

end


text:
export sample!
using StatsBase: sample, Weights

"""
    sample!(model::ABM, n [, weight]; kwargs...)

Replace the agents of the `model` with a random sample of the current agents with
size `n`.

Optionally, provide a `weight`: Symbol (agent field) or function (input agent
out put number) to weight the sampling.
This means that the higher the `weight` of the agent, the higher the probability that
this agent will be chosen in the new sampling.

# Keywords
* `replace = true` : whether sampling is performed with replacement, i.e. all agents can
be chosen more than once.
* `rng = GLOBAL_RNG` : a random number generator to perform the sampling with.

Example usage in [Wright-Fisher model of evolution](@ref).
"""
function sample!(
    model::ABM,
    n::Int,
    weight = nothing;
    replace = true,
    rng::AbstractRNG = Random.GLOBAL_RNG,
)
    nagents(model) > 0 || return
    org_ids = collect(keys(model.agents))
    if weight != nothing
        weights = Weights([get_data(a, weight, identity) for a in values(model.agents)])
        newids = sample(rng, org_ids, weights, n, replace = replace)
    else
        newids = sample(rng, org_ids, n, replace = replace)
    end
    add_newids!(model, org_ids, newids)
end

"Used in sample!"
function add_newids!(model, org_ids, newids)
    n = nextid(model)
    for id in org_ids
        if !in(id, newids)
            kill_agent!(model.agents[id], model)
        else
            noccurances = count(x -> x == id, newids)
            for t in 2:noccurances
                newagent = deepcopy(model.agents[id])
                newagent.id = n
                add_agent_pos!(newagent, model)
                n += 1
            end
        end
    end
end

ext: jl


text:
# Solves the primal problem of a soft-margin SVM classifier using
# the ACCPM algorithm.

# The main function.
# This implementation is deterministic for illustration purpose.
# In practice we may want to use random starting point.
function svmaccpm(X::AbstractMatrix,
                  Y::AbstractVector,
                  C::Real,
                  ϵ::Real,
                  droppingstrategy::Symbol=:None;
                  α::Real=0.01, β::Real=0.05)
  w = zeros(size(X, 2))
  z = ones(Y)*2
  x0 = [w; z]
  A, b = initialpolytope(X, Y, C, w, z)
  wsize = size(w, 1)
  f(x) = sum(x[1: wsize].^2)/2 + C*sum(x[wsize+1: end])
  ∇f(x) = [x[1: size(w, 1)]; C*ones(Y)]
  x, optdistances, numconstraints = accpm(
    f, ∇f, A, b, x0, ϵ, droppingstrategy; α=α, β=β)
  x[1: wsize], optdistances, numconstraints
end

# Computes the initial polyhedron for the SVM problem with some given
# initial point.
# Notes that in the original problem formulation the polyhedron assocaited
# with the constraint is not bounded, but once given some possible value
# of the objective we can easily give an upper bound to |w| and z.
function initialpolytope(X::AbstractMatrix,
                         Y::AbstractVector,
                         C::Real,
                         w::AbstractVector,
                         z::AbstractVector)
  objvalue = sum(w.^2)/2 + C*sum(z)
  wupper = √(2*objvalue)
  zupper = objvalue/C
  A = [-Y.*X (-eye(size(Y, 1)));
       zeros(X) (-eye(size(Y, 1)));
       -eye(size(w, 1)) zeros(X)';
       eye(size(w, 1)+size(z, 1));]
  b = [-1.*ones(Y); zeros(Y);
       fill(wupper, 2*size(w, 1)); fill(zupper, size(z))]
  A, b
end

# The general ACCPM algorithm to minimize `f` using the first-order
# derivative information.
# We suppose a linear constraint Ax ≤ b exists for the original problem.
function accpm(f, ∇f,
               A::AbstractMatrix,
               b::AbstractVector,
               x0::AbstractVector,
               ϵ::Real,
               droppingstrategy::Symbol=:None;
               α::Real=0.01, β::Real=0.05)
  (all(A*x0.<b)
   || throw(DomainError("x0 must be strictly feasible.")))
  (droppingstrategy in [:None, :DropRedundant, :KeepConstant]
   || throw(ArgumentError(
        "droppingstrategy must be one of the three:"
        * ":None, :DropRedundant or :KeepConstant.")))
  x = copy(x0)
  ms = [size(A, 1)]
  u, l = f(x), -Inf
  optdistances::Vector{Float64} = []
  while true
    x, H = polytope_analcenter(A, b, x, 1e-6; α=α, β=β)
    ∇, fx = ∇f(x), f(x)
    u = min(u, fx)
    l = max(l, fx - ms[end]*√(∇'*(H\∇)))
    push!(optdistances, u - l)
    u - l < ϵ && break
    if droppingstrategy ≠ :None
      A, b = dropconstraint(A, b, x, H, droppingstrategy)
    end
    A, b = [A; ∇'], [b; ∇'*x]
    x = nextstartpoint(A, b, x)
    push!(ms, size(A, 1))
  end
  x, optdistances, ms
end

# In the ACCPM algorithm, at the end of each iteration, the vector `x`
# is found on the face of the polytope. However, to run the Newton method
# that gives the analytic center, a strictly feasible point must be given.
# This function supposes that we have Ax ≤ b but with eqaulity only
# on the last line and try to move `x` a little to get a new `x'` such
# that Ax' ≤ b (the last line to `A` is suppose to be not all zero).
function nextstartpoint(A::AbstractMatrix,
                        b::AbstractVector,
                        x::AbstractVector)
  a = A[end, :]
  upperbound = (b-A*x)[1: end-1]
  lowerbound = max.(0, -(A*a)[1: end-1])
  δ = minimum(upperbound./lowerbound)
  x - δ*a/2
end

# Drop the constraints according to different strategies.
#
# :DropRedundant:
#   Guarantees not to change the polytope in question.
#
# :KeepConstant:
#   Keeps a fixed number of constraints which are the most relevant.
#   Here I choose to keep 3n constraints where n is the variable dimension.
#   However, we shouldn't end up with something that doesn't satisfy the
#   initial constraints, so this implementation can have bugs in some
#   particular cases.
function dropconstraint(A::AbstractMatrix,
                        b::AbstractVector,
                        x::AbstractVector,
                        H::AbstractMatrix,
                        droppingstrategy::Symbol)
  m = size(A, 1)
  Hinv = inv(H)
  relevemeasure = (b-A*x)./[√((A*Hinv)[i, :]⋅A[i, :]) for i=1:m]
  if droppingstrategy == :DropRedundant
    tokeep = relevemeasure .< m
  end
  if droppingstrategy == :KeepConstant
    tokeep = sortperm(relevemeasure)[1: min(3*size(A, 2), end)]
  end
  A[tokeep, :], b[tokeep]
end


text:
@testset "half disk" begin
    ϕ = π/6  # start angle
    θ = ϕ-π/2  # outward normal to diameter
    c = 0.01 .* [cos(θ), sin(θ)]  # center of disk
    r = 2  # radius
    ∆ϕ = π  # half disk
    b = Sector(c, r, ϕ, ∆ϕ)  # note sector includes origin

    ∆vₛ = r * [cos(ϕ), sin(ϕ)]  # start of arc
    ∆vₑ = r * [cos(ϕ+∆ϕ), sin(ϕ+∆ϕ)]  # end of arc
    ∆vs = (∆vₛ, ∆vₑ)

    @test b == deepcopy(b)
    @test b ≈ Sector(c, r, ϕ+2π, ∆ϕ)
    @test b ≈ Sector(c, r, ϕ-2π, ∆ϕ)
    @test b ≈ Sector(c, r, ϕ+∆ϕ, -∆ϕ)
    @test b ≈ Sector(c, r, ϕ+∆ϕ-2π, -∆ϕ)
    @test hash(b) == hash(deepcopy(b))
    @test [0,0] ∈ b  # origin
    @test c ∈ b  # center
    @test (s = one⁻; s.*(∆vₛ+c) ∈ b && s.*(∆vₑ+c) ∈ b)  # slightly inside of end points of arc
    @test (s = one⁺; s.*(∆vₑ+c) ∉ b && s.*(∆vₑ+c) ∉ b)  # slightly outside of end points of arc

    @test all([(xn = surfpt_nearby(s.*c,b); xn ≈ (c, [cos(θ), sin(θ)])) for s = (one⁻⁻, 1, one⁺⁺)])  # handle point at center properly

    @test all([(p₀ = ∆vₛ/2 + c; p = s.*p₀; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1]) ≈ ϕ && n ≈ [cos(θ), sin(θ)]) for s = (one⁻⁻,1,one⁺⁺)])  # around start side
    @test all([(p₀ = ∆vₑ/2 + c; p = s.*p₀; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1]) ≈ ϕ+∆ϕ-2π && n ≈ [cos(θ), sin(θ)]) for s = (one⁻⁻,1,one⁺⁺)])  # around end side
    @test all([(α = ϕ+0.7∆ϕ; d̂ = [cos(α),sin(α)]; p = s.*r.*d̂ + c; surfpt_nearby(p,b) ≈ (r.*d̂.+c, d̂)) for s = (one⁻⁻,1,one⁺⁺)])  # around arc

    @test all([(p₀ = ∆v + c; p = one⁺⁺.*p₀; surfpt_nearby(p,b) ≈ (p₀, normalize(p-p₀))) for ∆v = ∆vs])  # outside end points of arc
    @test (p = ∆vₛ + c; surfpt_nearby(p,b) ≈ (∆vₛ+c, normalize([cos(ϕ),sin(ϕ)]+[cos(θ),sin(θ)])))   # at start point of arc
    @test (p = ∆vₑ + c; surfpt_nearby(p,b) ≈ (∆vₑ+c, normalize([cos(ϕ+∆ϕ),sin(ϕ+∆ϕ)]+[cos(θ),sin(θ)])))   # at end point of arc

    @test bounds(b) ≈ ([c[1]-r, c[2]+r*sin(ϕ+∆ϕ)], [c[1]+r*cos(ϕ), c[2]+r])
    @test checkbounds(b)

    @test (∆ = rand(2); translate(b,∆) ≈ Sector(c+∆, r, ϕ, ∆ϕ))
end  # @testset "half disk"

@testset "Sector, acute angle" begin
    ϕ = π/6  # start angle
    θ₁ = ϕ-π/2  # outward normal to start edge
    ∆ϕ = π/6  # half disk
    θ₂ = (ϕ+∆ϕ) + π/2  # outward normal to end edge

    ϕ₀ = ϕ + ∆ϕ/2
    c = -0.01 .* [cos(ϕ₀), sin(ϕ₀)]  # center
    r = 2  # radius
    b = Sector(c, r, ϕ, ∆ϕ)  # note sector includes origin

    ∆vₛ = r * [cos(ϕ), sin(ϕ)]  # start of arc
    ∆vₑ = r * [cos(ϕ+∆ϕ), sin(ϕ+∆ϕ)]  # end of arc
    ∆vs = (∆vₛ, ∆vₑ)

    @test b == deepcopy(b)
    @test b ≈ Sector(c, r, ϕ+2π, ∆ϕ)
    @test b ≈ Sector(c, r, ϕ-2π, ∆ϕ)
    @test b ≈ Sector(c, r, ϕ+∆ϕ, -∆ϕ)
    @test b ≈ Sector(c, r, ϕ+∆ϕ-2π, -∆ϕ)
    @test hash(b) == hash(deepcopy(b))
    @test [0,0] ∈ b  # origin
    @test c ∈ b  # center
    @test (s = one⁻; s.*(∆vₛ+c) ∈ b && s.*(∆vₑ+c) ∈ b)  # slightly outside of end points of arc
    @test (s = one⁺; s.*(∆vₑ+c) ∉ b && s.*(∆vₑ+c) ∉ b)  # slightly inside of end points of arc

    @test surfpt_nearby(c,b) ≈ (c, -[cos(ϕ₀), sin(ϕ₀)])  # handle point at center properly
    @test all([(p = c+0.1r.*[cos(θ),sin(θ)]; surfpt_nearby(p,b) ≈ (c, normalize(p-c))) for θ = (ϕ+π, ϕ+∆ϕ+π, ϕ₀+π, ϕ+π+0.1∆ϕ, ϕ+∆ϕ+π-0.1∆ϕ)])  # outside center

    @test all([(p₀ = ∆vₛ/2 + c; p = s.*p₀; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1]) ≈ ϕ && n ≈ [cos(θ₁), sin(θ₁)]) for s = (one⁻⁻,1,one⁺⁺)])  # around start side
    @test all([(p₀ = ∆vₑ/2 + c; p = s.*p₀; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1]) ≈ ϕ+∆ϕ && n ≈ [cos(θ₂), sin(θ₂)]) for s = (one⁻⁻,1,one⁺⁺)])  # around end side
    @test all([(α = ϕ+0.7∆ϕ; d̂ = [cos(α),sin(α)]; p = s.*r.*d̂ + c; surfpt_nearby(p,b) ≈ (r.*d̂.+c, d̂)) for s = (one⁻⁻,1,one⁺⁺)])  # around arc

    @test all([(p₀ = ∆v + c; p = one⁺⁺.*p₀; surfpt_nearby(p,b) ≈ (p₀, normalize(p-p₀))) for ∆v = ∆vs])  # outside end points of arc
    @test (p = ∆vₛ + c; surfpt_nearby(p,b) ≈ (∆vₛ+c, normalize([cos(ϕ),sin(ϕ)]+[cos(θ₁),sin(θ₁)])))   # at start point of arc
    @test (p = ∆vₑ + c; surfpt_nearby(p,b) ≈ (∆vₑ+c, normalize([cos(ϕ+∆ϕ),sin(ϕ+∆ϕ)]+[cos(θ₂),sin(θ₂)])))   # at end point of arc

    @test bounds(b) ≈ (c, [c[1]+r*cos(ϕ), c[2]+r*sin(ϕ+∆ϕ)])
    @test checkbounds(b)

    @test (∆ = rand(2); translate(b,∆) ≈ Sector(c+∆, r, ϕ, ∆ϕ))
end  # @testset "Sector, acute angle"


text: include("simplepolysos.jl")


text:
"""

# `check_dct_type(`

Check if dct == Dict[] and has length > 0.

### Method
```julia
check_dct_type((dct)
```

### Required arguments
```julia
* `dct::Dict{String, Any}`      : Observed data or parameter init data
```

"""
function check_dct_type(dct)
  if typeof(dct) <: Array && length(dct) > 0
    if typeof(dct) <: Array{Dict{S, T}} where {S <: AbstractString, T <: Any}
      return true
    end
  else
    if typeof(dct) <: Dict{S, T} where {S <: AbstractString, T <: Any}
      return true
    end
  end
  @error "Input $(dct) not of type Dict{String, Any}, but $(typeof(dct))"
  return false
end

"""

# update_R_file

Rewrite a dictionary of observed data or initial parameter values in R dump file
format to a file.

### Method
```julia
update_R_file{T<:Any}(file, dct)
```

### Required arguments
```julia
* `file::String`                : R file name
* `dct::Dict{String, T}`        : Dictionary to format in R
```

"""
function update_R_file(file::String, dct::Dict{String, T}) where T <: Any
  isfile(file) && rm(file)
	open(file, "w") do strmout
    str = ""
    for entry in dct
      str = "\""*entry[1]*"\" <- "
      val = entry[2]
      if length(val)==0 && length(size(val))==1          # Empty array
        str = str*"c()\n"
      elseif length(val)==1 && length(size(val))==0          # Scalar
        str = str*"$(val)\n"
          #elseif length(val)==1 && length(size(val))==1    # Single element vector
        #str = str*"$(val[1])\n"
      elseif length(val)>=1 && length(size(val))==1  # Vector
        str = str*"structure(c("
        write(strmout, str)
        str = ""
            writedlm(strmout, val', ',')
        str = str*"), .Dim=c($(length(val))))\n"
      elseif length(val)>1 && length(size(val))>1      # Array
        str = str*"structure(c("
        write(strmout, str)
        str = ""
            writedlm(strmout, val[:]', ',')
        dimstr = "c"*string(size(val))
        str = str*"), .Dim=$(dimstr))\n"
      end
      write(strmout, str)
    end
  end
end


text:
module TestJuliaPrograms

using Test

using Catlab, Catlab.Doctrines, Catlab.WiringDiagrams
using Catlab.Programs.JuliaPrograms

@present C(FreeCartesianCategory) begin
  W::Ob
  X::Ob
  Y::Ob
  Z::Ob
  f::Hom(X,Y)
  g::Hom(Y,Z)
  h::Hom(Z,W)
  l::Hom(otimes(Z,W),otimes(Y,X))
  m::Hom(otimes(X,Y),otimes(W,Z))
  n::Hom(otimes(W,Z),otimes(X,Y))
end

# Compilation
#############

X, Y, f, g, h = generators(C, [:X, :Y, :f, :g, :h])

@test compile(f) isa Function
@test compile(compose(f,g)) isa Function
@test compile(otimes(f,h)) isa Function
@test compile(compose(f, mcopy(Y), otimes(g,g))) isa Function

# Parsing
#########

f, g, m, n = generators(C, [:f, :g, :m, :n])

diagram = @parse_wiring_diagram(C, (x::X) -> f(x))
@test diagram == to_wiring_diagram(f)

diagram = @parse_wiring_diagram(C, (x::X) -> g(f(x)))
@test diagram == to_wiring_diagram(compose(f,g))

diagram = @parse_wiring_diagram C (x::X) begin
  y = f(x)
  z = g(y)
  return z
end
@test diagram == to_wiring_diagram(compose(f,g))

diagram = @parse_wiring_diagram C (v::X) begin
  v = f(v)
  v = g(v)
  v = h(v)
  v
end
@test is_permuted_equal(diagram, to_wiring_diagram(compose(f,g,h)), [2,3,1])

diagram = @parse_wiring_diagram C (x::X, y::Y) begin
  w, z = m(x, y)
  x2, y2 = n(w, z)
  return (x2, y2)
end
@test diagram == to_wiring_diagram(compose(m,n))

diagram = @parse_wiring_diagram(C, (x::X, y::Y) -> n(m(x,y)))
@test diagram == to_wiring_diagram(compose(m,n))

I = munit(FreeCartesianCategory.Ob)
c = add_generator!(C, Hom(:c, I, X))
d = add_generator!(C, Hom(:d, X, I))

diagram = @parse_wiring_diagram(C, () -> c())
@test diagram == to_wiring_diagram(c)

diagram = @parse_wiring_diagram(C, (x::X) -> c(d(x)))
@test diagram == to_wiring_diagram(compose(d,c))

diagram = @parse_wiring_diagram C (x::X) begin
  d(x)
  c()
end
@test diagram == to_wiring_diagram(compose(d,c))

# Roundtrip
###########

""" Convert morphism expression to Julia code and then back to wiring diagram.
"""
function roundtrip(f::HomExpr)::WiringDiagram
  expr = compile_expr(f, arg_types=first.(collect(dom(f))))
  parse_wiring_diagram(C, expr)
end

function test_roundtrip(f::HomExpr)
  @test roundtrip(f) == to_wiring_diagram(f)
end

f, g, h, l, m, n = generators(C, [:f, :g, :h, :l, :m, :n])

test_roundtrip(f)
test_roundtrip(compose(f,g))
test_roundtrip(otimes(f,h))
test_roundtrip(compose(m,n))
test_roundtrip(compose(l, braid(Y,X), m))
test_roundtrip(compose(mcopy(X), otimes(f,f)))

f = compose(f, mcopy(Y), otimes(g,g))
@test is_permuted_equal(roundtrip(f), to_wiring_diagram(f), [2,3,1])

end


text:
### A Pluto.jl notebook ###
# v0.18.1

using Markdown
using InteractiveUtils

# ╔═╡ 6cd99333-2b6f-470b-b140-07e4b2693e3e
using Unitful

# ╔═╡ 2ac4937e-7cab-48e3-9c76-8f860d2a2bd1
begin
	using CommonMark,
		LaTeXStrings,
		Lazy,
		PlutoUI,
		HypertextLiteral,
		MarkdownLiteral,
		Hyperscript,
		Random

	import MarkdownLiteral: @mdx

	macro css_str(stylesheet_url)
		the_css = open(download(stylesheet_url)) do io
			read(io, String)
		end
		return the_css
	end

	cm"Load support macros"

end

# ╔═╡ bde82185-3f3e-4633-849f-f1d65f063c11
cm"""
## Given

During steady-state operation, a gearbox receives 60 kW through
the input shaft and delivers power through the
output shaft.

{style="width: 65%; display:block; margin-inline:auto;"}
![gearbox](https://gitlab.com/tupeux/faire-thermodynamique-pluton/-/raw/main/notebooks/assets/ms6.x02.004.transp.png)

{ .fig_caption }
Funky little gearbox.

For the gearbox as the system, the rate of energy transfer by convection is

```math
\require{physics}
\begin{equation}
  \dot Q_{convection.out} = hA \qty( T_b - T_\infty )
\end{equation}
```

``\qquad`` where

```math
\begin{align}
	h &= \pu{0.171 kW/m2 \cdot K} \textsf{ is the heat transfer coefficient;} \\[3pt]
	A_{\it surface} &= \pu{1.0m2} \textsf{ is the outer surface area of the gearbox,} \\[3pt]
	T_{\it surface} &= \pu{27°C} \textsf{ is the temperature at the outer surface; and} \\[3pt]
	T_\infty &=  \pu{20°C} \textsf{ is the ambient temperature of the air.}
\end{align}
```
"""

# ╔═╡ ed521013-8bc1-485c-8574-fd2a7a9fc918
@mdx("""

## Find

For the gearbox, evaluate the heat transfer rate and the power delivered through the output shaft, each in kW.

""")

# ╔═╡ 58c14931-ff78-4a25-883e-4eb1dff0a15f
cm"""
## Get organized

These things you need to notice right away:
1. We were provided an expression for heat transfer that can be resolved completely from the information in **Given**{.Vollkorn .addblue}.
1. Steady state operation: _all time derivatives collapse to zero._{.killerorange}
2. All of the work is shaft work.
3. Multiple energy flows: check and verify you account for all of them.

Let's attack this problem with a _lazy_ approach. Lazy evaluation is a computer science concept that essentially means we will not calculate any expressions until we must. In other words, put away the calculator until one of these events calls you to action:
 -  you're terrified by imminent great harm; or
 -  you're intoxicated by the glee of great riches.

Apply your energy balance (1st Law)

 ```math
\begin{gather}
    \dot E_{in}
&-& \dot E_{out}
&=& {\color{#db5728}{\cancelto{\rlap{\textsf{ 0 : steady~state}}}{\class{hailmary}{\dv t}}}} E_{sys}

\\[6pt]

	\underbrace{\dot W_{\textit{input.shaft}}}_{\sf\large given}
&-& \underbrace{\dot W_{\textit{output.shaft}}}_{\sf\large target} -
    \underbrace{\dot Q_{\textit{convection.out}}}_{\sf\large ready.target}
&=& 0

\end{gather}
```

See? When you're organized, riches come quickly.

1. Remember, we can calculate ``\dot Q_{out}`` without the Energy Eqn, so we will:

   ```math
   \dot Q_{out} = hA \qty( T_b - T_\infty );
   ```

1. Finally, we close out the problem with a reality check,
   our energy equation, and little arithmetic:

   ```math
   \dot W_{\it output.shaft} = \dot W_{\it input.shaft} - \dot Q_{out}
   ```

   Pause for reality check. Our equation for the rate of energy (_power_ ≡ _rate of energy_ : your choice) leaving the generator suggests that the rate of energy brought in by the input shaft is equal to the sum of the rates of energy leaving the generator via heat transfer and the output shaft, respectively.

   Ask yourself often: _Am I making sense?_{.killerorange}

"""

# ╔═╡ 78f3f73c-b4fe-4492-9fc4-f74d083458da
cm"""
## Attack

It's all up to you (and Julia) now."""

# ╔═╡ 643b1d7d-2e4f-4a5d-9737-38111c07aa3a
# given

begin
	Ẇ_input_shaft_in = 60u"kW"
	h = 0.171u"kW/(m^2*K)"
	A_surface = 1.0u"m^2"
	T_surface = 27u"°C"
	T∞ = 20u"°C"
end;

# ╔═╡ becdc317-dd8f-47de-9804-70b695f3752e
Q̇_convection_out = h * A_surface * (T_surface - T∞)

# ╔═╡ ed727ea8-645f-4a9f-a3a9-dd428ad70563
Ẇ_output_shaft_out = Ẇ_input_shaft_in - Q̇_convection_out

# ╔═╡ 8348343b-ca7f-4663-9135-12c687689502
cm"## _finis_"

# ╔═╡ 3ded8340-bffe-421a-b6d1-777a4a907d88
L"""
\class{theanswer}{\begin{align}
\dot Q_{\it convection.out} &= \pu{%$Q̇_convection_out}   \\
\dot W_{\it output.shaft}   &= \pu{%$Ẇ_output_shaft_out}
\end{align}
}
"""

# ╔═╡ 4de68ffe-265f-4823-8cd9-df7b6ad20de6
cm"

"

# ╔═╡ c8ab6e9f-7c34-4ec4-a799-861741b18def
cm"

----

## Appendix"

# ╔═╡ 9270fc98-4ada-4e0a-ab45-255afe01c0e2
# TableOfContents()

# ╔═╡ 075e6adf-f818-47fa-aac1-60666df2f7d5
begin

	macro title_str(title::String, bookID::String="")
		cm_parser = CommonMark.Parser()
		enable!(cm_parser, AttributeRule)

		bookID == "" && return title |> cm_parser

		bib = Dict(
			"ms6" 	=> 	"""
						Moran, MJ & HN Shapiro (2008)
						Fundamentals of engineering thermodynamics 6e,
						978-0471787358
						https://gitlab.com/Good-Trouble/youCAN/doThermodynamics/-/raw/main/images/ms6cover.jpg
						""",
			"ms8" 	=> 	"""
						Moran, MJ Ed (2014)
						Fundamentals of engineering thermodynamics 8e,
						978-1118412930
						https://gitlab.com/Good-Trouble/youCAN/doThermodynamics/-/raw/main/images/ms8cover.jpg
						""",
			"ms9" 	=> 	"""
						Moran, MJ, _et al_ (2018)
						Fundamentals of engineering thermodynamics 9e,
						978-1119391388
						https://gitlab.com/Good-Trouble/youCAN/doThermodynamics/-/raw/main/images/ms9cover47x60.png
						""",
			"cb7" 	=> 	"""
						Çengel, YA & MA Boles (2011)
						Thermodynamics: An engineering approach 7e,
						978-0073529325
						https://gitlab.com/Good-Trouble/youCAN/doThermodynamics/-/raw/main/images/cb7cover.jpg
						""",
			"cb8" 	=> 	"""
						Çengel, YA & MA Boles (2015)
						Thermodynamics: An engineering approach 8e,
						978-0073398174
						https://gitlab.com/Good-Trouble/youCAN/doThermodynamics/-/raw/main/images/cb8cover.jpg
						""",
			"cb9" 	=> 	"""
						Çengel, YA, Boles, MA, & M Kanoglu (2019)
						Thermodynamics: An engineering approach 9e,
						978-1259822674
						https://gitlab.com/Good-Trouble/youCAN/doThermodynamics/-/raw/main/images/cb9cover48x60.jpg
						""",
		)

		haskey(bib, bookID) || return KeyError(bookID)

		the_bibliography, the_thumb =
			@as bb bib[bookID] begin
				split(bb, '\n', keepempty=false)
				strip.(bb)
				(
					join(bb[1:3], " _", "_ "),
					m("img", class="no-dark-invert", src=bb[4], title=join(bb[2:3]," "))

				)
			end

		out = @mdx("""

\n\n
			\n\n# $title\n\n{.source}\n$the_bibliography\n

			$the_thumb

			""")

		return out

	end

	  # ═══ load css ════════════════════════════════════════════*═

	@mdx("""

Load stylesheets

	{.appendix-item}
	Load Julia \"infrastructure\" pkgs
	""")

end

# ╔═╡ d5694646-e439-43ee-b2ef-29ef9741bb5a
title"Example 2.4"ms6

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CommonMark = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
Hyperscript = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
Lazy = "50d2b5c4-7a5e-59d5-8109-a42b560f39c0"
MarkdownLiteral = "736d6165-7244-6769-4267-6b50796e6954"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

[compat]
CommonMark = "~0.8.6"
Hyperscript = "~0.0.4"
HypertextLiteral = "~0.9.3"
LaTeXStrings = "~1.3.0"
Lazy = "~0.15.1"
MarkdownLiteral = "~0.1.1"
PlutoUI = "~0.7.37"
Unitful = "~1.11.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CommonMark]]
deps = ["Crayons", "JSON", "URIs"]
git-tree-sha1 = "4cd7063c9bdebdbd55ede1af70f3c2f48fab4215"
uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
version = "0.8.6"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.3.0"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"

[[deps.Lazy]]
deps = ["MacroTools"]
git-tree-sha1 = "1370f8202dac30758f3c345f9909b97f53d87d3f"
uuid = "50d2b5c4-7a5e-59d5-8109-a42b560f39c0"
version = "0.15.1"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MarkdownLiteral]]
deps = ["CommonMark", "HypertextLiteral"]
git-tree-sha1 = "0d3fa2dd374934b62ee16a4721fe68c418b92899"
uuid = "736d6165-7244-6769-4267-6b50796e6954"
version = "0.1.1"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "bf0a1121af131d9974241ba53f601211e9303a9e"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.37"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Unitful]]
deps = ["ConstructionBase", "Dates", "LinearAlgebra", "Random"]
git-tree-sha1 = "b649200e887a487468b71821e2644382699f1b0f"
uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
version = "1.11.0"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─d5694646-e439-43ee-b2ef-29ef9741bb5a
# ╠═bde82185-3f3e-4633-849f-f1d65f063c11
# ╟─ed521013-8bc1-485c-8574-fd2a7a9fc918
# ╟─58c14931-ff78-4a25-883e-4eb1dff0a15f
# ╟─78f3f73c-b4fe-4492-9fc4-f74d083458da
# ╠═6cd99333-2b6f-470b-b140-07e4b2693e3e
# ╠═643b1d7d-2e4f-4a5d-9737-38111c07aa3a
# ╠═becdc317-dd8f-47de-9804-70b695f3752e
# ╠═ed727ea8-645f-4a9f-a3a9-dd428ad70563
# ╟─8348343b-ca7f-4663-9135-12c687689502
# ╟─3ded8340-bffe-421a-b6d1-777a4a907d88
# ╟─4de68ffe-265f-4823-8cd9-df7b6ad20de6
# ╠═c8ab6e9f-7c34-4ec4-a799-861741b18def
# ╠═9270fc98-4ada-4e0a-ab45-255afe01c0e2
# ╟─2ac4937e-7cab-48e3-9c76-8f860d2a2bd1
# ╠═075e6adf-f818-47fa-aac1-60666df2f7d5
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: mdAshford/youCANdoThermodynamics

max_issues_repo_name: mdAshford/youCANdoThermodynamics

max_forks_repo_name: mdAshford/youCANdoThermodynamics


text:
function _computeVarElem(
    q, r,
    f::CDDirectDifferenceLoss,
    SxΔ, SyΔ,
    X, Y, row, col)

    nx, p = size(X)
    ny    = size(Y, 1)

    Sx = f.Σx
    Sy = f.Σy
    A = f.A

    fill!(q, 0.)
    fill!(r, 0.)

    if row == col
        # if diagonal element

        # compute t_aa
        @inbounds t_aa = A[row, row] - (Sy[row, row] - Sx[row, row])

        # compute qk
        for k=1:nx
            v = 0.
            for l=1:p
                @inbounds v += SyΔ[row, l] * X[k, l]
            end

            @inbounds q[k] = v * X[k, col] - Sy[row, row] + X[k, row]^2
        end

        # compute rk
        for k=1:ny

            v = 0.
            for l=1:p
                @inbounds v += SxΔ[row, l] * Y[k, l]
            end

            @inbounds r[k] = v * Y[k, row] - Y[k, row]^2 + Sx[row, row]
        end

        σ1 = sum(abs2, q) / (nx - 1) - nx / (nx - 1) * t_aa^2
        σ2 = sum(abs2, r) / (ny - 1) - ny / (ny - 1) * t_aa^2

        return σ1/nx + σ2/ny
    else
        # if off-diagonal element

        # compute t_ab
        @inbounds t_ab = A[row, col] + A[col, row] - (Sy[row, col] - Sx[row, col]) * 2.

        # compute qk
        for k=1:nx
            v1 = 0.
            v2 = 0.
            for l=1:p
                @inbounds v1 += SyΔ[row, l] * X[k, l]
                @inbounds v2 += SyΔ[col, l] * X[k, l]
            end

            @inbounds q[k] = (v1 * X[k, col] + v2 * X[k, row]) - (Sy[row, col] - X[k, row] * X[k, col]) * 2.
        end

        # compute rk
        for k=1:ny
            v1 = 0.
            v2 = 0.
            for l=1:p
                @inbounds v1 += SxΔ[row, l] * Y[k, l]
                @inbounds v2 += SxΔ[col, l] * Y[k, l]
            end

            @inbounds r[k] = (v2 * Y[k, row] + v1 * Y[k, col])  - (Y[k, row] * Y[k, col] - Sx[row, col]) * 2.
        end

        σ1 = sum(abs2, q) / (nx - 1) - nx / (nx - 1) * t_ab^2
        σ2 = sum(abs2, r) / (ny - 1) - ny / (ny - 1) * t_ab^2

        return (σ1/nx + σ2/ny) / 4.
    end
end

function _mul(S::Symmetric, θ::SymmetricSparseIterate)
    p = size(S, 1)

    out = zeros(p, p)
    for ci=1:p
        for ri=1:p
            for k=1:p
                @inbounds out[ri, ci] += S[ri, k] * θ[k, ci]
            end
        end
    end
    out
end

function _computeVar!(ω, f::CDDirectDifferenceLoss, X, Y, Δ)
    nx, p = size(X)
    ny    = size(Y, 1)

    q = zeros(nx)
    r = zeros(ny)

    SxΔ = _mul(f.Σx, Δ)
    SyΔ = _mul(f.Σy, Δ)

    for col=1:p
        for row=col:p
            linInd = sub2indLowerTriangular(p, row, col)
            ω[linInd] = sqrt(_computeVarElem(q, r, f, SxΔ, SyΔ, X, Y, row, col))
        end
    end

end

# our method
function diffEstimation(Sx::Symmetric, Sy::Symmetric, X, Y, λ, options=CDOptions())
    f = CDDirectDifferenceLoss(Sx, Sy)
    x = SymmetricSparseIterate(f.p)

    ##################
    #
    #  first stage
    #
    ##################
    # compute initial variance
    ω = Array{eltype(Sx)}(undef, length(x))
    _computeVar!(ω, f, X, Y, x)

    # compute initial estimate
    g = ProxL1(λ, ω)
    coordinateDescent!(x, f, g, options)

    ##################
    #
    #  second stage
    #
    ##################

    # recompute variance
    opt1 = CDOptions(;
        maxIter=options.maxIter,
        optTol=options.optTol,
        randomize=options.randomize,
        warmStart=true,
        numSteps=options.numSteps)

    _computeVar!(ω, f, X, Y, x)

    # recompute estimate
    g = ProxL1(λ, ω)
    coordinateDescent!(x, f, g, opt1)

    return x
end


text:
using Sobol, Test

# compare results with results from C++ code sobol.cc published on
# http://web.maths.unsw.edu.au/~fkuo/sobol/
# with new-joe-kuo-6.21201 file used as input.
#
# Command line used to generate output is
#  ./sobol $N 1024 new-joe-kuo-6.21201 > exp_results_$N
#
# For each of the dimensions below
# (except 21201 where only 64 samples are generated)

@testset "published results" begin
    dimensions = [5, 10, 20, 50, 100, 500, 21201]

    for dim in dimensions
        println("Testing dimension $(dim)")
        open(joinpath(dirname(@__FILE__), "results", "exp_results_$(dim)")) do exp_results_file
            s = SobolSeq(dim)
            x = zeros(dim)
            for line in eachline(exp_results_file)
                values = [parse(Float64, item) for item in split(line)]
                if length(values) > 0
                    @test x == values
                    next!(s, x)
                end
            end
        end
    end
end

using Base.Iterators: take
@testset "iterators" begin
    # issue #8
    @test [x[1] for x in collect(take(Sobol.SobolSeq(1),5))] == [0.5,0.75,0.25,0.375,0.875]
end

@testset "scaled" begin
    # ScaledSobolSeq constructors
    lb = [-1,0,0]
    ub = [1,3,2]
    N = length(lb)
    s = SobolSeq(lb,ub)
    @test s isa ScaledSobolSeq{3,Int}
    @test eltype(s) == Vector{Float64}
    @test eltype(SobolSeq(Float32.(lb),Float32.(ub))) == Vector{Float32}
    @test first(s) == [0,1.5,1]
    @test first(SobolSeq((x for x in lb), (x for x in ub))) == [0,1.5,1]
    @test SobolSeq(N,lb,ub) isa ScaledSobolSeq{3,Int}
    @test_throws BoundsError SobolSeq(2,lb,ub)
end

@testset "skip exactly n" begin
    # issue #22
    # want to get the nth element
    n = 10

    s1 = SobolSeq(1)
    x1 = [next!(s1)[1] for i in 1:n][n]

    s2 = SobolSeq(1)
    skip(s2, n-1, exact=true)
    x2 = next!(s2)[1]

    @test x1 == x2
end


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Plots
end;

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing6.11"

# ╔═╡ 1aa41898-3afa-11eb-2f49-570853ce6435
begin
	mu, sig = 2, 3
	eta = sqrt(3)*sig/pi
	n, N = 15, 10^7
	dNormal   = Normal(mu, sig)
	dLogistic = Logistic(mu, eta)
	xGrid = -8:0.1:12

	sNormal   = [var(rand(dNormal,n)) for _ in 1:N]
	sLogistic = [var(rand(dLogistic,n)) for _ in 1:N]
end;

# ╔═╡ 2eb0846a-3b2b-11eb-1154-ff7ca71f3c27
begin
	p1 = plot(xGrid, pdf.(dNormal,xGrid), c=:blue, label="Normal")
	p1 = plot!(xGrid, pdf.(dLogistic,xGrid), c=:red, label="Logistic",
		xlabel="x",ylabel="Density", xlims= (-8,12), ylims=(0,0.16))

	p2 = stephist(sNormal, bins=200, c=:blue, normed=true, label="Normal")
	p2 = stephist!(sLogistic, bins=200, c=:red, normed=true, label="Logistic",
		xlabel="Sample Variance", ylabel="Density", xlims=(0,30), ylims=(0,0.14))

	plot(p1, p2, size=(800, 400))
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing6.11"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═1aa41898-3afa-11eb-2f49-570853ce6435
# ╠═2eb0846a-3b2b-11eb-1154-ff7ca71f3c27
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NetworkInterfaceTapConfigurationPropertiesFormat2 <: SwaggerModel
    virtualNetworkTap::Any # spec type: Union{ Nothing, VirtualNetworkTap } # spec name: virtualNetworkTap
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function NetworkInterfaceTapConfigurationPropertiesFormat2(;virtualNetworkTap=nothing, provisioningState=nothing)
        o = new()
        validate_property(NetworkInterfaceTapConfigurationPropertiesFormat2, Symbol("virtualNetworkTap"), virtualNetworkTap)
        setfield!(o, Symbol("virtualNetworkTap"), virtualNetworkTap)
        validate_property(NetworkInterfaceTapConfigurationPropertiesFormat2, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type NetworkInterfaceTapConfigurationPropertiesFormat2

const _property_map_NetworkInterfaceTapConfigurationPropertiesFormat2 = Dict{Symbol,Symbol}(Symbol("virtualNetworkTap")=>Symbol("virtualNetworkTap"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_NetworkInterfaceTapConfigurationPropertiesFormat2 = Dict{Symbol,String}(Symbol("virtualNetworkTap")=>"VirtualNetworkTap", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ NetworkInterfaceTapConfigurationPropertiesFormat2 }) = collect(keys(_property_map_NetworkInterfaceTapConfigurationPropertiesFormat2))
Swagger.property_type(::Type{ NetworkInterfaceTapConfigurationPropertiesFormat2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkInterfaceTapConfigurationPropertiesFormat2[name]))}
Swagger.field_name(::Type{ NetworkInterfaceTapConfigurationPropertiesFormat2 }, property_name::Symbol) =  _property_map_NetworkInterfaceTapConfigurationPropertiesFormat2[property_name]

function check_required(o::NetworkInterfaceTapConfigurationPropertiesFormat2)
    true
end

function validate_property(::Type{ NetworkInterfaceTapConfigurationPropertiesFormat2 }, name::Symbol, val)
end


text:
module GeneralizedConjugateResidualSolver

export GeneralizedConjugateResidual

using ..LinearSolvers
const LS = LinearSolvers
using ..MPIStateArrays: device, realview

using LinearAlgebra
using LazyArrays
using StaticArrays
using KernelAbstractions

"""
    GeneralizedConjugateResidual(K, Q; rtol, atol)

This is an object for solving linear systems using an iterative Krylov method.
The constructor parameter `K` is the number of steps after which the algorithm
is restarted (if it has not converged), `Q` is a reference state used only
to allocate the solver internal state, and `tolerance` specifies the convergence
criterion based on the relative residual norm. The amount of memory
required by the solver state is roughly `(2K + 2) * size(Q)`.
This object is intended to be passed to the [`linearsolve!`](@ref) command.

This uses the restarted Generalized Conjugate Residual method of Eisenstat (1983).

### References

    @article{eisenstat1983variational,
      title={Variational iterative methods for nonsymmetric systems of linear equations},
      author={Eisenstat, Stanley C and Elman, Howard C and Schultz, Martin H},
      journal={SIAM Journal on Numerical Analysis},
      volume={20},
      number={2},
      pages={345--357},
      year={1983},
      publisher={SIAM}
    }
"""
mutable struct GeneralizedConjugateResidual{K, T, AT} <:
               LS.AbstractIterativeLinearSolver
    residual::AT
    L_residual::AT
    p::NTuple{K, AT}
    L_p::NTuple{K, AT}
    alpha::MArray{Tuple{K}, T, 1, K}
    normsq::MArray{Tuple{K}, T, 1, K}
    rtol::T
    atol::T

    function GeneralizedConjugateResidual(
        K,
        Q::AT;
        rtol = √eps(eltype(AT)),
        atol = eps(eltype(AT)),
    ) where {AT}
        T = eltype(Q)

        residual = similar(Q)
        L_residual = similar(Q)
        p = ntuple(i -> similar(Q), K)
        L_p = ntuple(i -> similar(Q), K)
        alpha = @MArray zeros(K)
        normsq = @MArray zeros(K)

        new{K, T, AT}(residual, L_residual, p, L_p, alpha, normsq, rtol, atol)
    end
end

const weighted = false

function LS.initialize!(
    linearoperator!,
    Q,
    Qrhs,
    solver::GeneralizedConjugateResidual,
    args...,
)
    residual = solver.residual
    p = solver.p
    L_p = solver.L_p

    @assert size(Q) == size(residual)
    rtol, atol = solver.rtol, solver.atol

    threshold = rtol * norm(Qrhs, weighted)
    linearoperator!(residual, Q, args...)
    residual .-= Qrhs

    converged = false
    residual_norm = norm(residual, weighted)
    if residual_norm < threshold
        converged = true
        return converged, threshold
    end

    p[1] .= residual
    linearoperator!(L_p[1], p[1], args...)

    threshold = max(atol, threshold)

    converged, threshold
end

function LS.doiteration!(
    linearoperator!,
    Q,
    Qrhs,
    solver::GeneralizedConjugateResidual{K},
    threshold,
    args...,
) where {K}

    residual = solver.residual
    p = solver.p
    L_residual = solver.L_residual
    L_p = solver.L_p
    normsq = solver.normsq
    alpha = solver.alpha

    residual_norm = typemax(eltype(Q))
    for k in 1:K
        normsq[k] = norm(L_p[k], weighted)^2
        beta = -dot(residual, L_p[k], weighted) / normsq[k]

        Q .+= beta * p[k]
        residual .+= beta * L_p[k]

        residual_norm = norm(residual, weighted)

        if residual_norm <= threshold
            return (true, k, residual_norm)
        end

        linearoperator!(L_residual, residual, args...)

        for l in 1:k
            alpha[l] = -dot(L_residual, L_p[l], weighted) / normsq[l]
        end

        if k < K
            rv_nextp = realview(p[k + 1])
            rv_L_nextp = realview(L_p[k + 1])
        else # restart
            rv_nextp = realview(p[1])
            rv_L_nextp = realview(L_p[1])
        end

        rv_residual = realview(residual)
        rv_p = realview.(p)
        rv_L_p = realview.(L_p)
        rv_L_residual = realview(L_residual)

        groupsize = 256
        T = eltype(alpha)

        event = Event(device(Q))
        event = LS.linearcombination!(device(Q), groupsize)(
            rv_nextp,
            (one(T), alpha[1:k]...),
            (rv_residual, rv_p[1:k]...),
            false;
            ndrange = length(rv_nextp),
            dependencies = (event,),
        )

        event = LS.linearcombination!(device(Q), groupsize)(
            rv_L_nextp,
            (one(T), alpha[1:k]...),
            (rv_L_residual, rv_L_p[1:k]...),
            false;
            ndrange = length(rv_nextp),
            dependencies = (event,),
        )
        wait(device(Q), event)
    end

    (false, K, residual_norm)
end

end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
cd(@__DIR__)
using Pkg
Pkg.activate(".")
##

include("sp_model.jl")
include("plot_utils.jl")
#include("results_handling.jl")
include("evaluation_utils.jl")
##

using DataFrames
using CSV

offset = 6531
timesteps = 1:168

pv = CSV.read("timeseries/basic_example.csv", DataFrame)[timesteps .+ offset, 3]
wind = CSV.read("timeseries/basic_example.csv", DataFrame)[timesteps .+ offset, 4]
demand = CSV.read("timeseries/basic_example.csv", DataFrame)[timesteps .+ offset, 2]

##

pars = copy(default_es_pars)

##

using Statistics
average_hourly_demand = mean(demand)

pars[:recovery_time] = 24
pars[:c_storage] = 100.
pars[:c_pv] = 300.
pars[:c_wind] = 800.
pars[:c_sto_op] = 0.00001

heatdemand = copy(demand)./300.
heatdemand0 = zeros(length(demand))
es = define_energy_system(pv, wind, demand, heatdemand; p = pars, strict_flex = true)

##

using Cbc

sp = instantiate(es, no_flex_pseudo_sampler(), optimizer = Cbc.Optimizer)

##

optimize!(sp)

##

println("Objective value: $(objective_value(sp))")

println("PV: $(round(value.(sp[1, :u_pv]); digits = 2)) -- Wind: $(round(value.(sp[1, :u_wind]); digits = 2)) -- Battery: $(round(value.(sp[1, :u_storage]); digits = 2))")
println("Heat storage: $(round(value.(sp[1, :u_heat_storage]); digits = 2))")
fixed_investment = Dict((
    :u_pv => value.(sp[1, :u_pv]),
    :u_wind => value.(sp[1, :u_wind]),
    :u_storage => value.(sp[1, :u_storage]),
    :u_heatpump => value.(sp[1, :u_heatpump]),
    :u_heat_storage => value.(sp[1, :u_heat_storage])
))

# for s in 1:minimum((n, 10))
#     println("Second stage objective value in scenario $s: $(objective_value(sp, s))")
# end

##

plot_results(sp, pv, wind, demand, hd = heatdemand, s = 1, stage_1 = [:gci, :gco], stage_2 = [:gci2, :gco2])
##

n = 100
F_max = average_hourly_demand * 0.1 # Have ancillary services equal to 10% of our typical demand
t_max = length(pv) - es.parameters[2].defaults[:recovery_time]

scens = simple_flex_sampler(n, F_max, t_max)
##
es_fixed = define_energy_system(pv, wind, demand, heatdemand; p = pars, strict_flex = true)
sp_f = instantiate(es_fixed, scens, optimizer = Cbc.Optimizer)
fix_investment!(sp_f, string.(keys(fixed_investment)), values(fixed_investment))

optimize!(sp_f)
##

optimal_decision(sp_f)
plot_results(sp_f, pv, wind, demand, hd = heatdemand, s = 8, stage_1 = [:gci, :gco], stage_2 = [:gci2, :gco2], inv_dict = fixed_investment)
##
cost_pos, pot_pos, cost_neg, pot_neg = test_decision(sp, 1:t_max)

plot_flexibility(1:t_max, cost_pos, pot_pos, cost_neg, pot_neg, objective_value(sp))

##
cost_pos_f, pot_pos_f, cost_neg_f, pot_neg_f = test_decision(sp_f, 1:t_max)

plot_flexibility(1:t_max, cost_pos_f, pot_pos_f, cost_neg_f, pot_neg_f, objective_value(sp))

#label_av = ["+ flexibility in scenario-aware system", "+ flexibility in scenario-unaware system"]
begin
plt_av = plot();
flexibility_availability(plt_av, pot_pos, label = "positive flexbility unaware");
flexibility_availability(plt_av, pot_pos_f, label = "positive fleibility aware");
flexibility_availability(plt_av, pot_neg, label = "negative flexibility unaware");
flexibility_availability(plt_av, pot_neg_f, label = "negative flexibility aware");
display(plt_av)
end
##


text:
# To create a vocabulary file
# for a given text file it counts all the words and writes result to a given save file
using ArgParse

function count_words(file::AbstractString; wfile=nothing)
    res = Dict{AbstractString, Int}()
    open(file) do f
        for line in eachline(f)
            words = split(line)
            for word in words
                if word in keys(res)
                    res[word] += 1
                else
                    res[word] = 1
                end
            end
        end
    end

    sorted_pairs = sort(collect(res), by=x->res[x[1]], rev=true)
    if wfile != nothing
        open(wfile, "w") do f
            for (k, v) in sorted_pairs
                write(f, "$v $k\n")
            end
        end
    else
        for (k, v) in sorted_pairs; println("$v $k");end;
    end

end

function main(args=ARGS)
    s = ArgParseSettings()
    s.exc_handler = ArgParse.debug_handler
    @add_arg_table s begin
        ("--textfile"; required=true; help="Raw text file containing all the information")
        ("--output"; help="Output count file contains the word frequency information")

    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    count_words(o[:textfile]; wfile=o[:output])

end
!isinteractive() && main(ARGS)

max_stars_repo_name: kirnap/KU-dependencyParser

max_issues_repo_name: kirnap/KU-dependencyParser

max_forks_repo_name: kirnap/KU-dependencyParser


text:
# Autogenerated wrapper script for CVRPSEP_jll for x86_64-apple-darwin
export libcvrpsep

JLLWrappers.@generate_wrapper_header("CVRPSEP")
JLLWrappers.@declare_library_product(libcvrpsep, "@rpath/libcvrpsep.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcvrpsep,
        "lib/libcvrpsep.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Sixel

export sixel_encode, sixel_decode

using ImageCore
using OffsetArrays
using IndirectArrays # sixel sequence is actually an indexed image format

import REPL: Terminals

include("interface.jl")
using .SixelInterface
const SEC = AbstractSixelEncoder
const SDC = AbstractSixelDecoder

include("terminaltools.jl")
using .TerminalTools

include("backend/libsixel/LibSixel.jl")
using .LibSixel

# Eventually we will rewrite everything in pure Julia :)
default_encoder(::AbstractArray) = LibSixel.LibSixelEncoder()
default_decoder() = LibSixel.LibSixelDecoder()

# The high-level API to deal with different Julia input types.
#
# We assume that encoder backends perfectly understands `io::IO` and `img::Array`. Any Julia-specific
# fancy array types (e.g., transpose, view) or lazy generators, if supported, should be handled here
# instead of in the backends.
include("encoder.jl")
include("decoder.jl")

# various frontend and input type supports
include("frontend/fileio.jl")

# Ref: https://invisible-island.net/xterm/ctlseqs/ctlseqs.html
"""
    is_sixel_supported(tty=stdout)::Bool

Check if given terminal `tty` supports sixel format.

!!! warning
    (Experiment) The return value is not fully tested on all terminals and all platforms.
"""
function is_sixel_supported(tty::Base.TTY=stdout)
    '4' in TerminalTools.query_terminal("\033[0c", tty)
end
is_sixel_supported(io::IO) = false

end # module


text:
#=
This file is auto-generated. Do not edit.
=#

mutable struct LoadZones <: Topology
    number::Int64
    name::String
    buses::Vector{Bus}
    maxactivepower::Float64
    maxreactivepower::Float64
    _forecasts::InfrastructureSystems.Forecasts
    internal::InfrastructureSystemsInternal
end

function LoadZones(number, name, buses, maxactivepower, maxreactivepower, _forecasts=InfrastructureSystems.Forecasts(), )
    LoadZones(number, name, buses, maxactivepower, maxreactivepower, _forecasts, InfrastructureSystemsInternal())
end

function LoadZones(; number, name, buses, maxactivepower, maxreactivepower, _forecasts=InfrastructureSystems.Forecasts(), )
    LoadZones(number, name, buses, maxactivepower, maxreactivepower, _forecasts, )
end

# Constructor for demo purposes; non-functional.

function LoadZones(::Nothing)
    LoadZones(;
        number=0,
        name="init",
        buses=[Bus(nothing)],
        maxactivepower=0.0,
        maxreactivepower=0.0,
        _forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get LoadZones number."""
get_number(value::LoadZones) = value.number
"""Get LoadZones name."""
get_name(value::LoadZones) = value.name
"""Get LoadZones buses."""
get_buses(value::LoadZones) = value.buses
"""Get LoadZones maxactivepower."""
get_maxactivepower(value::LoadZones) = value.maxactivepower
"""Get LoadZones maxreactivepower."""
get_maxreactivepower(value::LoadZones) = value.maxreactivepower
"""Get LoadZones _forecasts."""
get__forecasts(value::LoadZones) = value._forecasts
"""Get LoadZones internal."""
get_internal(value::LoadZones) = value.internal


text:
# Parameters

L = 5 # m
v = 5 # m/s

# Dynamics model functions

function next_state(x, θ, ϕ; dt = 0.05)
    # x is in meters and θ is in degrees
	# ϕ is steering angle input (deg)

	# Dynamics model
	ẋ = v * sind(θ)
	θ̇ = (v / L) * tand(ϕ)

	x′ = x + ẋ * dt
	θ′ = θ + rad2deg(θ̇) * dt

	return x′, θ′
end

function next_state(x, θ, y, ϕ; dt = 0.05)
    # for plotting purposes, we might also want the y value
    # x is in meters and θ is in degrees
	# ϕ is steering angle input (deg)

	# Dynamics model
	ẋ = v * sind(θ)
	ẏ = v * cosd(θ)
	θ̇ = (v / L) * tand(ϕ)

	x′ = x + ẋ * dt
	y′ = y + ẏ * dt
	θ′ = θ + rad2deg(θ̇) * dt

	return x′, θ′, y′
end

function next_state(x, θ, y, ϕ, cosval; dt = 0.05)
    # for plotting purposes, we might also want the y value
    # x is in meters and θ is in degrees
	# ϕ is steering angle input (deg)

	# Dynamics model
	ẋ = v * sind(θ)
	ẏ = v * cosval
	θ̇ = (v / L) * tand(ϕ)

	x′ = x + ẋ * dt
	y′ = y + ẏ * dt
	θ′ = θ + rad2deg(θ̇) * dt

	return x′, θ′, y′
end

function sim_steps(x, θ, ϕ; dt = 0.05, num_steps = 20)
    for i = 1:num_steps
		x, θ = next_state(x, θ, ϕ)
	end
    return x, θ
end

function sim_steps(x, θ, y, ϕ; dt = 0.05, num_steps = 20)
    for i = 1:num_steps
		x, θ, y = next_state(x, θ, y, ϕ)
	end
    return x, θ, y
end

function sim_steps(x, θ, y, ϕ, cosval; dt = 0.05, num_steps = 20)
    for i = 1:num_steps
		x, θ, y = next_state(x, θ, y, ϕ)
	end
    return x, θ, y
end

# NOTE: this function is only valid for θs between -90 and 90 degrees
function reachable_cell(lbs, ubs, ϕ_min, ϕ_max)
	next_lb_x, next_lb_θ = sim_steps(lbs[1], lbs[2], ϕ_min)
	next_ub_x, next_ub_θ = sim_steps(ubs[1], ubs[2], ϕ_max)
	return [next_lb_x, next_lb_θ], [next_ub_x, next_ub_θ]
end

function cosdextrema(xmin, xmax)
	b1, b2 = cosd(xmin), cosd(xmax)
	curr_min, curr_max = min(b1, b2), max(b1, b2)
	if xmin < 0 && 0 < xmax
		curr_max = 1
	end
	curr_min, curr_max
end

function reachable_cell_dtp(lbs, ubs, ϕ_min, ϕ_max; num_steps = 20)
	cosdmin, cosdmax = cosdextrema(ϕ_min, ϕ_max)

	next_lb_x, next_lb_θ, next_lb_y = sim_steps(lbs[1], lbs[2], lbs[3], ϕ_min, cosdmin, num_steps = num_steps)
	next_ub_x, next_ub_θ, next_ub_y = sim_steps(ubs[1], ubs[2], ubs[3], ϕ_max, cosdmax, num_steps = num_steps)
	return [next_lb_x, next_lb_θ, next_lb_y], [next_ub_x, next_ub_θ, next_ub_y]
end

max_stars_repo_name: sisl/VerifyGAN

max_issues_repo_name: sisl/VerifyGAN

max_forks_repo_name: sisl/VerifyGAN


text:
module Aswing
using OffsetArrays

import Base.eltype

function __init__()
    global UNITS = UNITS_s()
    global ASWING = ASWING_s()
    global EIGEN = EIGEN_s()
    global THCALC = THCALC_s()
    global CONLAW = CONLAW_s()
    initialize()

    global GlobalUNITS = UNITS
    global GlobalASWING = ASWING
    global GlobalEIGEN = EIGEN
    global GlobalTHCALC = THCALC
    global GlobalCONLAW = CONLAW
end

const depsfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("Aswing not properly installed. Please run Pkg.build(\"Aswing\") then restart Julia.")
end

function initialize()
    ccall((:init_, libaswing), Nothing, ())
    ccall((:setdef_, libaswing), Nothing, ())
    ccall((:iidef_, libaswing), Nothing, ())
    ccall((:nndefp_, libaswing), Nothing, ())
    ccall((:vginit_, libaswing), Nothing, ())
    ccall((:dscheck_, libaswing), Nothing, ())
    THCALC.CUSTOMINTERVALS[1] = false
    THCALC.INTERVALS[0:NBRKX+1] .= 0
end

# Global variable wrappers and fortran pointer definitions
include("globals/dimen.jl")
include("globals/indexb.jl")
include("globals/indexc.jl")
include("globals/indexp.jl")
include("globals/units.jl")
include("globals/aswing.jl")
include("globals/eigen.jl")
include("globals/thcalc.jl")
include("globals/conlaw.jl")

# defines ASWING geometry input
include("geom.jl")

# defines ASWING operating point definitions
include("pnt.jl")

# defines Aswing solution output
include("sol.jl")

# implements Aswing file i/o
include("io.jl")

# implements main Aswing calling methods
include("main.jl")

# deprecated functions
include("deprecated.jl")

end #modules


text:
using LightXML
using Base.Test

# document

docstr = """
<?xml version="1.0" encoding="UTF-8"?>
<bookstore>
  <book category="COOKING" tag="first">
    <title lang="en">Everyday Italian</title>
    <author>Giada De Laurentiis</author>
    <year>2005</year>
    <price>30.00</price>
  </book>
  <book category="CHILDREN">
    <title lang="en">Harry Potter</title>
    <author>J K. Rowling</author>
    <year>2005</year>
    <price>29.99</price>
  </book>
</bookstore>
"""

xdoc = parse_string(docstr)

@test version(xdoc) == "1.0"
@test encoding(xdoc) == "UTF-8"
@test standalone(xdoc) == -2

# root node

xroot = root(xdoc)

@test isa(xroot, XMLElement)
@test is_elementnode(xroot)
@test name(xroot) == "bookstore"
@test nodetype(xroot) == 1
@test !has_attributes(xroot)
@test has_children(xroot)

ras = collect(attributes(xroot))
@test isempty(ras)

# children of root (text nodes and books)

rcs = collect(child_nodes(xroot))
@test length(rcs) == 5  # text, book[1], text, book[1], text

@test is_textnode(rcs[1])
@test is_textnode(rcs[3])
@test is_textnode(rcs[5])

@test is_blanknode(rcs[1])
@test is_blanknode(rcs[3])
@test is_blanknode(rcs[5])

@test is_elementnode(rcs[2])
@test is_elementnode(rcs[4])

@test !is_blanknode(rcs[2])
@test !is_blanknode(rcs[4])

xb1 = XMLElement(rcs[2])

@test name(xb1) == "book"
@test nodetype(xb1) == 1
@test has_attributes(xb1)
@test has_children(xb1)
@test attribute(xb1, "category") == "COOKING"
@test attribute(xb1, "tag") == "first"

b1as = collect(attributes(xb1))
@test length(b1as) == 2

b1a1 = b1as[1]
@test isa(b1a1, XMLAttr)
@test name(b1a1) == "category"
@test value(b1a1) == "COOKING"

b1a2 = b1as[2]
@test isa(b1a2, XMLAttr)
@test name(b1a2) == "tag"
@test value(b1a2) == "first"

adct = attributes_dict(xb1)
@test length(adct) == 2
@test adct["category"] == "COOKING"
@test adct["tag"] == "first"

xb2 = XMLElement(rcs[4])

@test name(xb2) == "book"
@test nodetype(xb2) == 1
@test has_attributes(xb2)
@test has_children(xb2)
@test has_attribute(xb2, "category")
@test attribute(xb2, "category") == "CHILDREN"

@test !has_attribute(xb2, "wrongattr")
@test is(attribute(xb2, "wrongattr"), nothing)
@test_throws LightXML.XMLAttributeNotFound attribute(xb2, "wrongattr"; required=true)

rces = get_elements_by_tagname(xroot, "book")
@test length(rces) == 2
@test isa(rces, Vector{XMLElement})
@test attribute(rces[1], "category") == "COOKING"
@test attribute(rces[2], "category") == "CHILDREN"

# child elements of book[1]

ces = collect(child_elements(xb1))

@test length(ces) == 4
c1, c2, c3, c4 = ces[1], ces[2], ces[3], ces[4]

@test isa(c1, XMLElement)
@test name(c1) == "title"
@test has_attributes(c1)
@test attribute(c1, "lang") == "en"
@test content(c1) == "Everyday Italian"

@test has_children(c1)
c1cs = collect(child_nodes(c1))
@test length(c1cs) == 1
c1c = c1cs[1]
@test is_textnode(c1c)
@test content(c1c) == "Everyday Italian"

@test isa(c2, XMLElement)
@test name(c2) == "author"
@test !has_attributes(c2)
@test content(c2) == "Giada De Laurentiis"

@test isa(c3, XMLElement)
@test name(c3) == "year"
@test !has_attributes(c3)
@test content(c3) == "2005"

@test isa(c4, XMLElement)
@test name(c4) == "price"
@test !has_attributes(c4)
@test content(c4) == "30.00"

cy = find_element(xb1, "year")
@test isa(cy, XMLElement)
@test name(cy) == "year"
@test content(cy) == "2005"

cz = find_element(xb1, "abc")
@test is(cz, nothing)

free(xdoc)


text:
module Test03

import Automa
import Automa.RegExp: @re_str
using Test

@testset "Test03" begin
    re = Automa.RegExp
    header = re"[ -~]*"
    newline = re"\r?\n"
    sequence = re.rep(re.cat(re"[A-Za-z]*", newline))
    fasta = re.rep(re.cat('>', header, newline, sequence))

    machine = Automa.compile(fasta)

    for generator in (:table, :inline, :goto), checkbounds in (true, false), clean in (true, false)
        ctx = Automa.CodeGenContext(generator=generator, checkbounds=checkbounds, clean=clean)
        init_code = Automa.generate_init_code(ctx, machine)
        exec_code = Automa.generate_exec_code(ctx, machine)
        validate = @eval function (data)
            $(init_code)
            p_end = p_eof = lastindex(data)
            $(exec_code)
            return cs == 0
        end
        @test validate(b"") == true
        @test validate(b">\naa\n") == true
        @test validate(b">seq1\n") == true
        @test validate(b">seq1\na\n") == true
        @test validate(b">seq1\nac\ngt\n") == true
        @test validate(b">seq1\r\nacgt\r\n") == true
        @test validate(b">seq1\nac\n>seq2\ngt\n") == true
        @test validate(b"a") == false
        @test validate(b">") == false
        @test validate(b">seq1\na") == false
        @test validate(b">seq1\nac\ngt") == false
    end
end

end


text:
function Th=EnergyToPot(U,CG,Param)
nz=Param.Grid.nz;
RhoPos=Param.RhoPos;
uPos=Param.uPos;
vPos=Param.vPos;
wPos=Param.wPos;
ThPos=Param.ThPos;
W=[zeros(size(U,1),1)...
   ,U(:,:,wPos)];
WC=0.5*(W(:,1:nz)+W(:,2:nz+1));
KE=0.5*(U(:,:,uPos).*U(:,:,uPos)+U(:,:,vPos).*U(:,:,vPos)+WC.*WC);
Temp=(U(:,:,ThPos)-U(:,:,RhoPos).*(KE...
  +Param.Grav*repmat(Param.Grid.zP,1,size(U,1))'))/Param.Cvd;
Th=Temp.*(Param.p0./(Param.Rd.*Temp)).^(Param.Rd/Param.Cpd);
end


text:
using Distributions, Statistics, Random
rng = MersenneTwister(01012019);

"""
    kernel_scenarios(data::Union{Array{Int64,2},Array{Float64,2}}, kernelDistribution = Logistic; Markovian::Bool = true)

Returns an instance of a closure of the conditional density estimation method.

Args:
- data - a matrix of data in 2 dimension (NxT) where N is the number of trajectories and T is the number of stages.
- kernelDistribution - the distribution of the kernel functions. Default is Logistic distribution.

An optional keyword `Markovian` is passed to handle weights
for Markovian trajectories or non-Markovian ones.

To get a Markov trajectory from above,
`KernelScenarios(data,kernelDistribution;Markovian=true)()` gives the result for
Markovian samples, which is used to generate a scenario lattice.
"""
function kernel_scenarios(data::Union{Array{Int64,2},Array{Float64,2}}, kernelDistribution = Logistic; Markovian::Bool = true)
    function closure()
        N , T = size(data)                             # Dimensions of the data
        d = 1
        w = fill(1.0, N)                               # initialize the weights
        x = Array{Float64,2}(undef, T , 1)             # trajectory to be created
        for t = 1 : T
            w = w / sum(w)                             # normalized weights
            Nt = sum(w)^2 / sum(w.^2)                  # effective sample size
            σt = sqrt(sum(w .* ((data[:,t] .- sum(w .* data[:,t])).^2)))  # effective standard devistion
            ht = σt * Nt^(-1 / (d + 4)) + eps()
            # Composition method comes in here
            u = rand(rng, Uniform(0,1))
            _,jstar = findmin(cumsum(w) .< u * sum(w))
            x[t] = rand(rng, kernelDistribution(data[jstar,t], ht))
            if t<T
                if Markovian      # Markovian is used for scenario lattices
                    for j = 1:N
                        w[j] = pdf(kernelDistribution(data[j,t], ht), x[t])
                    end
                else             # Non-Markovian is used for scenario trees
                    for j = 1:N
                        w[j] *= pdf(kernelDistribution(data[j,t], ht), x[t])
                    end
                end
            end
        end
        return x        # The length of x is equal to the number stages = T
    end
end


text:
using AbstractTrees

mutable struct BinaryNode{T}
    data::T
    parent::BinaryNode{T}
    left::BinaryNode{T}
    right::BinaryNode{T}

    # Root constructor
    BinaryNode{T}(data) where T = new{T}(data)
    # Child node constructor
    BinaryNode{T}(data, parent::BinaryNode{T}) where T = new{T}(data, parent)
end
BinaryNode(data) = BinaryNode{typeof(data)}(data)

function leftchild(data, parent::BinaryNode)
    !isdefined(parent, :left) || error("left child is already assigned")
    node = typeof(parent)(data, parent)
    parent.left = node
end
function rightchild(data, parent::BinaryNode)
    !isdefined(parent, :right) || error("right child is already assigned")
    node = typeof(parent)(data, parent)
    parent.right = node
end

## Things we need to define
function AbstractTrees.children(node::BinaryNode)
    if isdefined(node, :left)
        if isdefined(node, :right)
            return (node.left, node.right)
        end
        return (node.left,)
    end
    isdefined(node, :right) && return (node.right,)
    return ()
end

## Things that make printing prettier
AbstractTrees.printnode(io::IO, node::BinaryNode) = print(io, node.data)

## Optional enhancements
# These next two definitions allow inference of the item type in iteration.
# (They are not sufficient to solve all internal inference issues, however.)
Base.eltype(::Type{<:TreeIterator{BinaryNode{T}}}) where T = BinaryNode{T}
Base.IteratorEltype(::Type{<:TreeIterator{BinaryNode{T}}}) where T = Base.HasEltype()

## Let's test it. First build a tree.
root = BinaryNode(0)
l = leftchild(1, root) # technically, this mutates root.
r = rightchild(2, root)
lr = rightchild(3, l)

print_tree(root)
collect(PostOrderDFS(root))
@static if isdefined(@__MODULE__, :Test)
    @testset "binarytree_easy.jl" begin
        @test [node.data for node in PostOrderDFS(root)] == [3, 1, 2, 0]
        @test [node.data for node in PreOrderDFS(root)] == [0, 1, 3, 2]
        @test [node.data for node in Leaves(root)] == [3, 2]
    end
end

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
mutable struct FitzHughNagumo <: SpecialDiffSys
    δ::Float64
    a::Float64
    b::Float64
    I::Float64
    noise::AbstractArray
end

name(::FitzHughNagumo) = "Fitz-Hugh Nagumo equation"
dimension(::FitzHughNagumo) = 2

FitzHughNagumo(; δ=0.1, a=1.0, b=0.8, I=0.8, noise=zeros(2)) = FitzHughNagumo(δ, a, b, I, noise)

function vec_field(fhn::FitzHughNagumo, x::AbstractArray)
    u, v = x
    δ = fhn.δ
    a = fhn.a
    b = fhn.b
    I = fhn.I
    du = δ*(v+a-b*u)
    dv = v-v^3/3.0-u+I
    return [du, dv]
end

function jacobian(fhn::FitzHughNagumo, x::AbstractArray)
    _, v = x
    δ = fhn.δ
    b = fhn.b
    J = [-δ*b δ; -1 1-v*v]
    return J
end


text:
using JLD
using PyCall
@pyimport corner
@pyimport numpy as py
@pyimport matplotlib.pyplot as plt

@load "tyc3559_mcmc_norv_fixedkernel_v03.jld"

#nsamp = 10000
nsamp = nwalkers*nsteps
npar = 8
#npar = 2
data = zeros(nsamp,npar)
ipar = [3,4,5,11,12,13,15,16]
#ipar = [3,4]
for i1=1:nwalkers
  for i2=1:nsteps
#  i1 = ceil(Int64,rand()*nwalkers)
#  i2 = ceil(Int64,rand()*nsteps)
    for j=1:npar
      data[(i1-1)*nsteps+i2,j]=par_mcmc[i1,i2,ipar[j]]
    end
  end
end

data[:,2] *= 100.
data[:,3] *= 24.
data[:,4] *= 1e6
data[:,5] *= 1e6
data[:,6] *= 1e6
data[:,7] *= 100.
data[:,8] *= 1e6

writedlm("tyc3559_mcmc_norv_fixedkernel_v03.csv",data,",")

#labname=[r"$b$",r"$D$",r"$T$",r"$a_{1c}$",r"$a_{1s}$",r"$a_{2c}$",r"$f_A/(f_A+f_B+f_{Bb})$",r"$f_{Bb}/(f_B+f_{Bb})$"]

#figure = corner.corner(data,labels=[r"\$b\$","$D$","$T$","$a_{1c}$","$a_{1s}$","$a_{2c}$","$f_A/(f_A+f_B+f_{Bb})$","$f_{Bb}/(f_B+f_{Bb})$"], quantiles=[0.1585, 0.5, 0.8415])
#figure = corner.corner(data, labels=PyAny([r"$b$",r"$D$"]), quantiles=[0.1585, 0.5, 0.8415])
#                       quantiles=[0.1585, 0.5, 0.8415],
#                       show_titles=True, title_kwargs={"fontsize": 12})

#figure = pycall(corner.corner,PyAny, data, labels=[r"$b$",r"$D$"], quantiles=[0.1585, 0.5, 0.8415])
#figure = pycall(corner.corner,PyAny, data, labels=['$b$','$D$'], quantiles=[0.1585, 0.5, 0.8415])
#plt.show(figure)


text:
import DataStructures
import StaticArrays
using BenchmarkTools

include("../src/soalight.jl")
include("../src/allocator.jl")
include("../src/queue.jl")

function testbase(n, p)
    q = DataStructures.Queue{Tuple{Int, typeof(p)}}()
    for j in 1:n[]
        DataStructures.enqueue!(q, (j, p))
    end
    for j in 1:n[]
        @assert j == DataStructures.dequeue!(q)[1]
    end
end

function testalloc(n, p, q)
    for j in 1:n[]
        push!(q, (j, p))
    end
    for j in 1:n[]
        @assert j == pop!(q)[1]
    end
    @assert isempty(q)
    emptyend!(q.alloc)
end

ans = nothing
for n in [100000, 10]
    println(n, " elements")
    for p in [1, 10, 100]
        println(" payload of ", p, " float(s)")

        Payload = StaticArrays.SVector{p ,Float64}

        GC.gc()

        print("  DataStructures.Queue           ")
        @btime ($ans = testbase(Ref($n), zeros($Payload)))

        GC.gc()

        print("  Queue with fixed allocator     ")
        alloc = Allocator{Tuple{Int, Payload}, Int}(n)
        q = Queue{Tuple{Int, Payload}, Int, Allocator{Tuple{Int, Payload}, Int}}(alloc)
        @btime ($ans = testalloc(Ref($n), zeros($Payload), $q))

        GC.gc()

        print("  Queue with resizable allocator ")
        alloc = Allocator{Tuple{Int, Payload}, Int}(nothing)
        q = Queue{Tuple{Int, Payload}, Int, Allocator{Tuple{Int, Payload}, Int}}(alloc)
        @btime ($ans = testalloc(Ref($n), zeros($Payload), $q))
    end
end


text:
using UnivariateFunctions
const tol = 10*eps()

# Make undefined function
und = Undefined_Function()
ismissing(und(5))
results = und.([5, 9]) # Test broadcasting
ismissing(results[1])
ismissing(results[2])

# pe functions
function test_pe(a::Float64,b::Float64,base::Float64,d::Int,x::Float64)
    return a * exp(b*(x-base)) * (x-base)^d
end
f1 = PE_Function(1.0, 2.0,3.0, 4)
f1_test_result = test_pe(1.0, 2.0,3.0, 4,5.0)
abs(evaluate(f1, 5.0) - f1_test_result) < tol
f2 = PE_Function(0.0, 2.0,3.0, 4)
f2_test_result = 0.0
abs(evaluate(f2, 5.0) - f2_test_result) < tol
f3 = PE_Function{Float64,Int64}(1.0, 8.0,7.0, 8)
f3_test_result = test_pe(1.0, 8.0,7.0, 8,5.0)
abs(f3(5.0) - f3_test_result) < tol
f4 = PE_Function(1.0)
f4_test_result = test_pe(1.0, 0.0,0.0, 0,5.0)
abs(f4(5.0) - f4_test_result) < tol

# Sum of functions
sum0 =  Sum_Of_Functions([])
typeof(sum0) == UnivariateFunctions.Sum_Of_Functions
sum1 =  Sum_Of_Functions([f1])
sum1 = Sum_Of_Functions(sum1)
typeof(sum1) == UnivariateFunctions.Sum_Of_Functions
sum3 =  Sum_Of_Functions([f2,f3])
typeof(sum3) == UnivariateFunctions.Sum_Of_Functions
sum4 =  Sum_Of_Functions([f1,f2,f3])
typeof(sum4) == UnivariateFunctions.Sum_Of_Functions
length(sum4.functions_) == 2
sum5 =  Sum_Of_Functions([sum3, sum4])
typeof(sum5) == UnivariateFunctions.Sum_Of_Functions
length(sum5.functions_) == 2
abs(sum1(5.0) - f1_test_result) < tol
abs(evaluate(sum5, 5.0) - 2*f3_test_result - f1_test_result) < 100* tol
abs(evaluate(sum0, 5.0)) < tol
sum6 =  Sum_Of_Functions([sum3, sum4, und])
typeof(sum6)  == Undefined_Function

fl = 8.9
integ = 7

function test_result(func, expected_type, eval_to, len = 1)
    val_test = abs(evaluate(func, 5.0) - eval_to) < 1e-09
    if (!val_test)
        print("Failed Val Test")
    end
    type_test = isa(func, expected_type)
    if (!type_test)
        print("Failed Type Test")
    end
    if typeof(func) == UnivariateFunctions.Sum_Of_Functions
        len_test = length(func.functions_) == len
        if (!len_test)
            print("Failed Length Test")
        end
        return all([val_test,type_test,len_test])
    else
        return all([val_test,type_test])
    end
end

# pe and Float
test_result(f1 + fl, UnivariateFunctions.Sum_Of_Functions, f1_test_result + fl,2 )
test_result(f1 - fl, UnivariateFunctions.Sum_Of_Functions, f1_test_result - fl,2 )
test_result(f1 * fl, UnivariateFunctions.PE_Function, (f1_test_result) * fl )
test_result(f1 / fl, UnivariateFunctions.PE_Function, (f1_test_result) / fl )
test_result(fl + f1, UnivariateFunctions.Sum_Of_Functions, f1_test_result + fl,2 )
test_result(fl - f1, UnivariateFunctions.Sum_Of_Functions, fl - f1_test_result,2 )
test_result(fl * f1, UnivariateFunctions.PE_Function, (f1_test_result) * fl )
# pe and Int
test_result(f1 + integ, UnivariateFunctions.Sum_Of_Functions, f1_test_result + integ,2 )
test_result(f1 - integ, UnivariateFunctions.Sum_Of_Functions, f1_test_result - integ,2 )
test_result(f1 * integ, UnivariateFunctions.PE_Function, (f1_test_result) * integ )
test_result(f1 / integ, UnivariateFunctions.PE_Function, (f1_test_result) / integ )
test_result(integ + f1, UnivariateFunctions.Sum_Of_Functions, f1_test_result + integ,2 )
test_result(integ - f1, UnivariateFunctions.Sum_Of_Functions, -f1_test_result + integ,2 )
test_result(integ * f1, UnivariateFunctions.PE_Function, (f1_test_result) * integ )

# Sum and Float
test_result(sum4 + fl, UnivariateFunctions.Sum_Of_Functions, f1_test_result + f3_test_result + fl, 3 )
test_result(sum4 - fl, UnivariateFunctions.Sum_Of_Functions, f1_test_result + f3_test_result - fl, 3 )
test_result(sum4 * fl, UnivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result) * fl, 2 )
test_result(sum4 / fl, UnivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result) / fl, 2 )
test_result(fl + sum4, UnivariateFunctions.Sum_Of_Functions, f1_test_result + f3_test_result + fl, 3 )
test_result(fl - sum4, UnivariateFunctions.Sum_Of_Functions, -f1_test_result -f3_test_result + fl, 3 )
test_result(fl * sum4, UnivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result) * fl, 2 )
# Sum and Int
test_result(sum4 + integ, UnivariateFunctions.Sum_Of_Functions, f1_test_result + f3_test_result + integ, 3 )
test_result(sum4 - integ, UnivariateFunctions.Sum_Of_Functions, f1_test_result + f3_test_result - integ, 3 )
test_result(sum4 * integ, UnivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result) * integ, 2 )
test_result(sum4 / integ, UnivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result) / integ, 2 )
test_result(integ + sum4, UnivariateFunctions.Sum_Of_Functions, f1_test_result + f3_test_result + integ, 3 )
test_result(integ - sum4, UnivariateFunctions.Sum_Of_Functions, -f1_test_result -f3_test_result + integ, 3 )
test_result(integ * sum4, UnivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result) * integ, 2 )

### Making sums with addition and subtraction.
evaluate(sum4, 5.0) == evaluate(f1 + f2 + f3, 5.0)
abs(evaluate(sum4, 5.0) - 2*evaluate(f3, 5.0) - evaluate(f1 + f2 - f3, 5.0)) < 0.001

### Changing of base
f1_unchanged = change_base_of_PE_Function(f1, 3.0)
isa(f1_unchanged, UnivariateFunctions.PE_Function)
abs(f1.base_ - f1_unchanged.base_) < tol

f1_changed = change_base_of_PE_Function(f1, 4.0)
typeof(f1_changed) == UnivariateFunctions.Sum_Of_Functions
abs(f1_changed.functions_[1].base_ - 4.0) < tol
abs(f1_changed.functions_[3].base_ - 4.0) < tol
abs(evaluate(f1_changed, 5.0) - f1_test_result) < 100*tol # Changing bases should not change this.

f1_changed_again = change_base_of_PE_Function(f1, -1.0)
typeof(f1_changed_again) == UnivariateFunctions.Sum_Of_Functions
abs(f1_changed_again.functions_[1].base_ + 1.0) < tol
abs(f1_changed_again.functions_[3].base_ + 1.0) < tol
abs(evaluate(f1_changed_again, 5.0) - f1_test_result) < 1000000*tol # Changing bases should not change this.

f3_changed = change_base_of_PE_Function(f3, 12.0)
typeof(f3_changed) == UnivariateFunctions.Sum_Of_Functions
abs(f3_changed.functions_[1].base_ - 12.0) < tol
abs(f3_changed.functions_[3].base_ - 12.0) < tol
abs(evaluate(f3_changed, 5.0) - f3_test_result) < 100*tol # Changing bases should not change this.

### multiplication of functions
test_result( f1 * f3 , UnivariateFunctions.Sum_Of_Functions, (f1_test_result * f3_test_result), 9 )
test_result( f1 * PE_Function(1.0,2.0,3.0,4) , UnivariateFunctions.PE_Function, (f1_test_result * f1_test_result), 1 )
test_result( f1 * sum4 , UnivariateFunctions.Sum_Of_Functions, f1_test_result * (f1_test_result + f3_test_result), 10 )
test_result( sum4 * f1 , UnivariateFunctions.Sum_Of_Functions, f1_test_result * (f1_test_result + f3_test_result), 10 )
test_result( sum4 * sum4 , UnivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result)^2, 11 )

### Powers pe
test_result( f1 ^ 0 , UnivariateFunctions.PE_Function, 1.0 )
test_result( f1 ^ 1 , UnivariateFunctions.PE_Function, (f1_test_result) )
test_result( f1 ^ 2 , UnivariateFunctions.PE_Function, (f1_test_result * f1_test_result) )
abs(evaluate(f1 ^ 4 ,5.0) - (f1_test_result * f1_test_result * f1_test_result * f1_test_result)) < 1e-03 # This fails normal tol.

### Powers sums
test_result( sum4 ^ 0 , UnivariateFunctions.PE_Function, 1.0 )
test_result( sum4 ^ 1 , UnivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result), 2 )
test_result( sum4 ^ 2 , UnivariateFunctions.Sum_Of_Functions, (f1_test_result + f3_test_result)^2, 11 )
# Higher powers give an underflow problem with base changes.

## Testing Calculus.
pe_const    = PE_Function(2.0,0.0,5.0,0)
pe_lin      = PE_Function(2.0,0.0,5.0,1)
pe_quad     = PE_Function(2.0,0.0,2.0,2)
pe_exp      = PE_Function(2.0,2.0,1.0,0)
pe_exp_quad = PE_Function(2.0,2.0,2.0,2)

# Linear gradient constant
abs(evaluate(derivative(pe_lin),5.0) - evaluate(derivative(pe_lin),1.0) ) < tol
typeof(derivative(pe_lin)) == UnivariateFunctions.Sum_Of_Functions
# This is also linear
abs(evaluate(derivative(derivative(pe_quad)),5.0) -evaluate(derivative(derivative(pe_quad)),9.0) ) < tol
# derivative into a sum of functions
typeof(derivative(pe_exp_quad)) == UnivariateFunctions.Sum_Of_Functions
abs(evaluate(derivative(pe_exp_quad),5.0) - ( 2.0*exp(2.0*(5.0-2.0))*(5.0-2.0)*(2 + 2.0*(5.0-2.0)) ) ) < tol
# Derivaitve of sum
abs(evaluate(derivative(sum4),5.0) - evaluate(derivative(f1),5.0) - evaluate(derivative(f3),5.0) ) < 1e-09

# integral of constant
indefinite_integral(pe_const).d_ == 1
# Integral of quadratic
abs(evaluate_integral(pe_quad,1.0,2.0) - ((2/3)*2^3 - 4*2^2 + 8*2) + ((2/3)*1^3 - 4*1^2 + 8*1)) < tol
# Integral of exponential
abs(evaluate_integral(pe_exp,0.2,3.8) - (exp(2.0*(3.8-1.0)) - exp(2.0*(0.2-1.0))  )) < tol
# Integral of combined.
combined_analytical_integral = PE_Function(1.0,2.0,2.0,2) - PE_Function(1.0,2.0,2.0,1) + PE_Function(0.5,2.0,2.0,0)
abs(evaluate_integral(pe_exp_quad,0.2,3.8) - ( evaluate(combined_analytical_integral, 3.8) - evaluate(combined_analytical_integral, 0.2))) < tol

# left and right integrals
l_int = left_integral(pe_exp_quad, 0.2)
(evaluate(l_int, 3.8) - evaluate_integral(pe_exp_quad,0.2,3.8)) < tol
r_int = right_integral(pe_exp_quad, 3.8)
(evaluate(r_int, 0.2) - evaluate_integral(pe_exp_quad,0.2,3.8)) < tol

## Testing linear rescaling
# For PE_Functions
test_func = PE_Function(0.5,1.2,0.9,2)
inp = 4.0
alpha = 0.435
beta = -1.52
rescaled_input = alpha*inp + beta
converted_test_func = convert_to_linearly_rescale_inputs(test_func, alpha, beta)
abs(evaluate(test_func, inp) - evaluate(converted_test_func, rescaled_input)) < 1e-10
# For undefined.
typeof(convert_to_linearly_rescale_inputs(Undefined_Function(), alpha, beta)) == UnivariateFunctions.Undefined_Function
# SumOfFunctions
sumFunc = test_func + pe_exp_quad + pe_exp
converted_sumFunc = convert_to_linearly_rescale_inputs(sumFunc, alpha, beta)
abs(evaluate(sumFunc, inp) - evaluate(converted_sumFunc, rescaled_input)) < 1e-10

### Test undefined
typeof(und + und)  == Undefined_Function
typeof(und - und)  == Undefined_Function
typeof(und / und)  == Undefined_Function
typeof(und * und)  == Undefined_Function

# PE_Function
typeof(f1 + und)  == Undefined_Function
typeof(f1 - und)  == Undefined_Function
typeof(f1 / und)  == Undefined_Function
typeof(f1 * und)  == Undefined_Function
typeof(und + f1)  == Undefined_Function
typeof(und - f1)  == Undefined_Function
typeof(und / f1)  == Undefined_Function
typeof(und * f1)  == Undefined_Function

# Sum_Of_Functions
typeof(sum0 + und)  == Undefined_Function
typeof(sum0 - und)  == Undefined_Function
typeof(sum0 / und)  == Undefined_Function
typeof(sum0 * und)  == Undefined_Function
typeof(und + sum0)  == Undefined_Function
typeof(und - sum0)  == Undefined_Function
typeof(und / sum0)  == Undefined_Function
typeof(und * sum0)  == Undefined_Function
typeof(und + und)  == Undefined_Function
typeof(und - und)  == Undefined_Function
typeof(und / und)  == Undefined_Function
typeof(und * und)  == Undefined_Function

# Numbers
typeof(und + 4)  == Undefined_Function
typeof(und - 4)  == Undefined_Function
typeof(und / 4)  == Undefined_Function
typeof(und * 4)  == Undefined_Function
typeof(und ^ 4)  == Undefined_Function
typeof(4 + und)  == Undefined_Function
typeof(4 - und)  == Undefined_Function
typeof(4 * und)  == Undefined_Function

typeof(und + 4.0)  == Undefined_Function
typeof(und - 4.0)  == Undefined_Function
typeof(und / 4.0)  == Undefined_Function
typeof(und * 4.0)  == Undefined_Function
typeof(4.0 + und)  == Undefined_Function
typeof(4.0 - und)  == Undefined_Function
typeof(4.0 * und)  == Undefined_Function


text:
export CUBLASError

struct CUBLASError <: Exception
    code::cublasStatus_t
    msg::AbstractString
end
Base.show(io::IO, err::CUBLASError) = print(io, "CUBLASError(code $(err.code), $(err.msg))")

function CUBLASError(code::cublasStatus_t)
    msg = status_message(code)
    return CUBLASError(code, msg)
end

function status_message(status)
    if status == CUBLAS_STATUS_SUCCESS
        return "the operation completed successfully"
    elseif status == CUBLAS_STATUS_NOT_INITIALIZED
        return "the library was not initialized"
    elseif status == CUBLAS_STATUS_ALLOC_FAILED
        return "the resource allocation failed"
    elseif status == CUBLAS_STATUS_INVALID_VALUE
        return "an invalid value was used as an argument"
    elseif status == CUBLAS_STATUS_ARCH_MISMATCH
        return "an absent device architectural feature is required"
    elseif status == CUBLAS_STATUS_MAPPING_ERROR
        return "an access to GPU memory space failed"
    elseif status == CUBLAS_STATUS_EXECUTION_FAILED
        return "the GPU program failed to execute"
    elseif status == CUBLAS_STATUS_INTERNAL_ERROR
        return "an internal operation failed"
    elseif status == CUBLAS_STATUS_NOT_SUPPORTED
        return "the requested feature is not supported"
    elseif status == CUBLAS_STATUS_LICENSE_ERROR
        return "error detected trying to check the license"
    else
        return "unknown status"
    end
end

macro check(blas_func)
    quote
        local err::cublasStatus_t
        err = $(esc(blas_func::Expr))
        if err != CUBLAS_STATUS_SUCCESS
            throw(CUBLASError(err))
        end
        err
    end
end


text:
abstract type AbstractBasis end

"""
    BasisCollocated()

Returns the singleton object corresponding to libCEED's `CEED_BASIS_COLLOCATED`.
"""
struct BasisCollocated <: AbstractBasis end
Base.getindex(::BasisCollocated) = C.CEED_BASIS_COLLOCATED[]

"""
    Basis

Wraps a `CeedBasis` object, representing a finite element basis. A `Basis` object can be
created using one of:

- [`create_tensor_h1_lagrange_basis`](@ref)
- [`create_tensor_h1_basis`](@ref)
- [`create_h1_basis`](@ref)
"""
mutable struct Basis <: AbstractBasis
    ref::RefValue{C.CeedBasis}
    function Basis(ref)
        obj = new(ref)
        finalizer(obj) do x
            # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x))
            destroy(x)
        end
        return obj
    end
end
destroy(b::Basis) = C.CeedBasisDestroy(b.ref) # COV_EXCL_LINE
Base.getindex(b::Basis) = b.ref[]
Base.show(io::IO, ::MIME"text/plain", b::Basis) = ceed_show(io, b, C.CeedBasisView)

@doc raw"""
    create_tensor_h1_lagrange_basis(ceed, dim, ncomp, p, q, qmode)

Create a tensor-product Lagrange basis.

# Arguments:
- `ceed`:  A [`Ceed`](@ref) object where the [`Basis`](@ref) will be created.
- `dim`:   Topological dimension of element.
- `ncomp`: Number of field components (1 for scalar fields).
- `p`:     Number of Gauss-Lobatto nodes in one dimension.  The polynomial degree of the
           resulting $Q_k$ element is $k=p-1$.
- `q`:     Number of quadrature points in one dimension.
- `qmode`: Distribution of the $q$ quadrature points (affects order of accuracy for the
           quadrature).
"""
function create_tensor_h1_lagrange_basis(c::Ceed, dim, ncomp, p, q, quad_mode::QuadMode)
    ref = Ref{C.CeedBasis}()
    C.CeedBasisCreateTensorH1Lagrange(c[], dim, ncomp, p, q, quad_mode, ref)
    Basis(ref)
end

@doc raw"""
    create_tensor_h1_basis(c::Ceed, dim, ncomp, p, q, interp1d, grad1d, qref1d, qweight1d)

Create a tensor-product basis for $H^1$ discretizations.

# Arguments:
- `ceed`:      A [`Ceed`](@ref) object where the [`Basis`](@ref) will be created.
- `dim`:       Topological dimension.
- `ncomp`:     Number of field components (1 for scalar fields).
- `p`:         Number of nodes in one dimension.
- `q`:         Number of quadrature points in one dimension
- `interp1d`:  Matrix of size `(q, p)` expressing the values of nodal basis functions at
               quadrature points.
- `grad1d`:    Matrix of size `(p, q)` expressing derivatives of nodal basis functions at
               quadrature points.
- `qref1d`:    Array of length `q` holding the locations of quadrature points on the 1D
               reference element $[-1, 1]$.
- `qweight1d`: Array of length `q` holding the quadrature weights on the reference element.
"""
function create_tensor_h1_basis(
    c::Ceed,
    dim,
    ncomp,
    p,
    q,
    interp1d,
    grad1d,
    qref1d,
    qweight1d,
)
    @assert size(interp1d) == (q, p)
    @assert size(grad1d) == (q, p)
    @assert length(qref1d) == q
    @assert length(qweight1d) == q

    # Convert from Julia matrices (column-major) to row-major format
    interp1d_rowmajor = collect(interp1d')
    grad1d_rowmajor = collect(grad1d')

    ref = Ref{C.CeedBasis}()
    C.CeedBasisCreateTensorH1(
        c[],
        dim,
        ncomp,
        p,
        q,
        interp1d_rowmajor,
        grad1d_rowmajor,
        qref1d,
        qweight1d,
        ref,
    )
    Basis(ref)
end

@doc raw"""
    create_h1_basis(c::Ceed, topo::Topology, ncomp, nnodes, nqpts, interp, grad, qref, qweight)

Create a non tensor-product basis for H^1 discretizations

# Arguments:
- `ceed`:    A [`Ceed`](@ref) object where the [`Basis`](@ref) will be created.
- `topo`:    [`Topology`](@ref) of element, e.g. hypercube, simplex, etc.
- `ncomp`:   Number of field components (1 for scalar fields).
- `nnodes`:  Total number of nodes.
- `nqpts`:   Total number of quadrature points.
- `interp`:  Matrix of size `(nqpts, nnodes)` expressing the values of nodal basis functions
             at quadrature points.
- `grad`:    Array of size `(dim, nqpts, nnodes)` expressing derivatives of nodal basis
             functions at quadrature points.
- `qref`:    Array of length `nqpts` holding the locations of quadrature points on the
             reference element $[-1, 1]$.
- `qweight`: Array of length `nqpts` holding the quadrature weights on the reference
             element.
"""
function create_h1_basis(
    c::Ceed,
    topo::Topology,
    ncomp,
    nnodes,
    nqpts,
    interp,
    grad,
    qref,
    qweight,
)
    dim = getdimension(topo)
    @assert size(interp) == (nqpts, nnodes)
    @assert size(grad) == (dim, nqpts, nnodes)
    @assert length(qref) == nqpts
    @assert length(qweight) == nqpts

    # Convert from Julia matrices and tensors (column-major) to row-major format
    interp_rowmajor = collect(interp')
    grad_rowmajor = permutedims(grad, [3, 2, 1])

    ref = Ref{C.CeedBasis}()
    C.CeedBasisCreateH1(
        c[],
        topo,
        ncomp,
        nnodes,
        nqpts,
        interp_rowmajor,
        grad_rowmajor,
        qref,
        qweight,
        ref,
    )
    Basis(ref)
end

"""
    apply!(b::Basis, nelem, tmode::TransposeMode, emode::EvalMode, u::AbstractCeedVector, v::AbstractCeedVector)

Apply basis evaluation from nodes to quadrature points or vice versa, storing the result in
the [`CeedVector`](@ref) `v`.

`nelem` specifies the number of elements to apply the basis evaluation to; the backend will
specify the ordering in CeedElemRestrictionCreateBlocked()

Set `tmode` to `CEED_NOTRANSPOSE` to evaluate from nodes to quadrature or to
`CEED_TRANSPOSE` to apply the transpose, mapping from quadrature points to nodes.

Set the [`EvalMode`](@ref) `emode` to:
- `CEED_EVAL_NONE` to use values directly,
- `CEED_EVAL_INTERP` to use interpolated values,
- `CEED_EVAL_GRAD` to use gradients,
- `CEED_EVAL_WEIGHT` to use quadrature weights.
"""
function apply!(
    b::Basis,
    nelem,
    tmode::TransposeMode,
    emode::EvalMode,
    u::AbstractCeedVector,
    v::AbstractCeedVector,
)
    C.CeedBasisApply(b[], nelem, tmode, emode, u[], v[])
end

"""
    apply(b::Basis, u::AbstractVector; nelem=1, tmode=NOTRANSPOSE, emode=EVAL_INTERP)

Performs the same function as the above-defined [`apply!`](@ref apply!(b::Basis, nelem,
tmode::TransposeMode, emode::EvalMode, u::AbstractCeedVector, v::AbstractCeedVector)), but
automatically convert from Julia arrays to [`CeedVector`](@ref) for convenience.

The result will be returned in a newly allocated array of the correct size.
"""
function apply(b::Basis, u::AbstractVector; nelem=1, tmode=NOTRANSPOSE, emode=EVAL_INTERP)
    ceed_ref = Ref{C.Ceed}()
    ccall((:CeedBasisGetCeed, C.libceed), Cint, (C.CeedBasis, Ptr{C.Ceed}), b[], ceed_ref)
    c = Ceed(ceed_ref)

    u_vec = CeedVector(c, u)

    len_v = (tmode == TRANSPOSE) ? getnumnodes(b) : getnumqpts(b)
    if emode == EVAL_GRAD
        len_v *= getdimension(b)
    end

    v_vec = CeedVector(c, len_v)

    apply!(b, nelem, tmode, emode, u_vec, v_vec)
    Vector(v_vec)
end

"""
    getdimension(b::Basis)

Return the spatial dimension of the given [`Basis`](@ref).
"""
function getdimension(b::Basis)
    dim = Ref{CeedInt}()
    C.CeedBasisGetDimension(b[], dim)
    dim[]
end

"""
    getdimension(t::Topology)

Return the spatial dimension of the given [`Topology`](@ref).
"""
function getdimension(t::Topology)
    return Int(t) >> 16
end

"""
    gettopology(b::Basis)

Return the [`Topology`](@ref) of the given [`Basis`](@ref).
"""
function gettopology(b::Basis)
    topo = Ref{Topology}()
    C.CeedBasisGetTopology(b[], topo)
    topo[]
end

"""
    getnumcomponents(b::Basis)

Return the number of components of the given [`Basis`](@ref).
"""
function getnumcomponents(b::Basis)
    ncomp = Ref{CeedInt}()
    C.CeedBasisGetNumComponents(b[], ncomp)
    ncomp[]
end

"""
    getnumnodes(b::Basis)

Return the number of nodes of the given [`Basis`](@ref).
"""
function getnumnodes(b::Basis)
    nnodes = Ref{CeedInt}()
    C.CeedBasisGetNumNodes(b[], nnodes)
    nnodes[]
end

"""
    getnumnodes1d(b::Basis)

    Return the number of 1D nodes of the given (tensor-product) [`Basis`](@ref).
"""
function getnumnodes1d(b::Basis)
    nnodes1d = Ref{CeedInt}()
    C.CeedBasisGetNumNodes1D(b[], nnodes1d)
    nnodes1d[]
end

"""
    getnumqpts(b::Basis)

Return the number of quadrature points of the given [`Basis`](@ref).
"""
function getnumqpts(b::Basis)
    nqpts = Ref{CeedInt}()
    C.CeedBasisGetNumQuadraturePoints(b[], nqpts)
    nqpts[]
end

"""
    getnumqpts1d(b::Basis)

Return the number of 1D quadrature points of the given (tensor-product) [`Basis`](@ref).
"""
function getnumqpts1d(b::Basis)
    nqpts1d = Ref{CeedInt}()
    C.CeedBasisGetNumQuadraturePoints1D(b[], nqpts1d)
    nqpts1d[]
end

"""
    getqref(b::Basis)

Get the reference coordinates of quadrature points (in `dim` dimensions) of the given
[`Basis`](@ref).
"""
function getqref(b::Basis)
    ref = Ref{Ptr{CeedScalar}}()
    C.CeedBasisGetQRef(b[], ref)
    copy(unsafe_wrap(Array, ref[], getnumqpts(b)))
end

"""
    getqref(b::Basis)

Get the quadrature weights of quadrature points (in `dim` dimensions) of the given
[`Basis`](@ref).
"""
function getqweights(b::Basis)
    ref = Ref{Ptr{CeedScalar}}()
    C.CeedBasisGetQWeights(b[], ref)
    copy(unsafe_wrap(Array, ref[], getnumqpts(b)))
end

"""
    getinterp(b::Basis)

Get the interpolation matrix of the given [`Basis`](@ref). Returns a matrix of size
`(getnumqpts(b), getnumnodes(b))`.
"""
function getinterp(b::Basis)
    ref = Ref{Ptr{CeedScalar}}()
    C.CeedBasisGetInterp(b[], ref)
    q = getnumqpts(b)
    p = getnumnodes(b)
    collect(unsafe_wrap(Array, ref[], (p, q))')
end

"""
    getinterp1d(b::Basis)

Get the 1D interpolation matrix of the given [`Basis`](@ref). `b` must be a tensor-product
basis, otherwise this function will fail. Returns a matrix of size `(getnumqpts1d(b),
getnumnodes1d(b))`.
"""
function getinterp1d(b::Basis)
    ref = Ref{Ptr{CeedScalar}}()
    C.CeedBasisGetInterp1D(b[], ref)
    q = getnumqpts1d(b)
    p = getnumnodes1d(b)
    collect(unsafe_wrap(Array, ref[], (p, q))')
end

"""
    getgad(b::Basis)

Get the gradient matrix of the given [`Basis`](@ref). Returns a tensor of size
`(getdimension(b), getnumqpts(b), getnumnodes(b))`.
"""
function getgrad(b::Basis)
    ref = Ref{Ptr{CeedScalar}}()
    C.CeedBasisGetGrad(b[], ref)
    dim = getdimension(b)
    q = getnumqpts(b)
    p = getnumnodes(b)
    permutedims(unsafe_wrap(Array, ref[], (p, q, dim)), [3, 2, 1])
end

"""
    getgrad1d(b::Basis)

Get the 1D derivative matrix of the given [`Basis`](@ref). Returns a matrix of size
`(getnumqpts(b), getnumnodes(b))`.
"""
function getgrad1d(b::Basis)
    ref = Ref{Ptr{CeedScalar}}()
    C.CeedBasisGetGrad1D(b[], ref)
    q = getnumqpts1d(b)
    p = getnumnodes1d(b)
    collect(unsafe_wrap(Array, ref[], (p, q))')
end

max_stars_repo_name: DiffeoInvariant/libCEED

max_issues_repo_name: DiffeoInvariant/libCEED

max_forks_repo_name: DiffeoInvariant/libCEED


text:
#
# This function computes the center of mass of a selection
#

function cm(selection,simulation,x,y,z)

  cm = [ 0., 0., 0. ]
  totmass = 0.
  for i in 1:length(selection)
    cm[1] = cm[1] + simulation.atoms[selection[i]].mass*x[selection[i]]
    cm[2] = cm[2] + simulation.atoms[selection[i]].mass*y[selection[i]]
    cm[3] = cm[3] + simulation.atoms[selection[i]].mass*z[selection[i]]
    totmass = totmass + simulation.atoms[selection[i]].mass
  end
  cm[1] = cm[1] / totmass
  cm[2] = cm[2] / totmass
  cm[3] = cm[3] / totmass

  return cm
end


text:
include("typs_and_structs.jl")

function calc_equilibrium_distribution(mom::moments1d_64, gravity=0)
    len = length(mom.height)
    v = 1.0
    vsquare = v^2
    equi_dist = distribution64_1d(zeros(len), zeros(len), zeros(len))

    for i in 1:len
        equi_dist.f0[i] = mom.height[i] - 1/(2*vsquare) * gravity * mom.height[i]^2 - 1/vsquare * mom.height[i] * mom.velocity[i]^2
        equi_dist.f1[i] = 1/(4*vsquare) * gravity * mom.height[i]^2 + 1/(2*v) * mom.height[i] * mom.velocity[i] + 1/(2*vsquare) * mom.height[i] * mom.velocity[i]^2
        equi_dist.f2[i] = 1/(4*vsquare) * gravity * mom.height[i]^2 - 1/(2*v) * mom.height[i] * mom.velocity[i] + 1/(2*vsquare) * mom.height[i] * mom.velocity[i]^2
    end
    return equi_dist
end

function calc_equilibrium_distribution(mom::moments1d_32, gravity=0)
    len = length(mom.height)
    v = 1.0
    vsquare = v^2
    equi_dist = distribution32_1d(zeros(len), zeros(len), zeros(len))

    for i in 1:len
        equi_dist.f0[i] = mom.height[i] - 1/(2*vsquare) * gravity * mom.height[i]^2 - 1/vsquare * mom.height[i] * mom.velocity[i]^2
        equi_dist.f1[i] = 1/(4*vsquare) * gravity * mom.height[i]^2 + 1/(2*v) * mom.height[i] * mom.velocity[i] + 1/(2*vsquare) * mom.height[i] * mom.velocity[i]^2
        equi_dist.f2[i] = 1/(4*vsquare) * gravity * mom.height[i]^2 - 1/(2*v) * mom.height[i] * mom.velocity[i] + 1/(2*vsquare) * mom.height[i] * mom.velocity[i]^2
    end
    return equi_dist
end