text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using SparseArrays
# using ADCMEKit
Random.seed!(233)

function sparse_factorization(A, s=99)
    ii, jj, vv = find(constant(A))
    d = size(A, 1)
    sparse_factorization_ = load_op_and_grad("./build/libfactorization","sparse_factorization")
    ii,jj,vv,d,s = convert_to_tensor([ii,jj,vv,d,s], [Int64,Int64,Float64,Int64,Int64])
    stop_gradient(sparse_factorization_(ii,jj,vv,d,s))
end

function sparse_solve(A,rhs,o)
    ii, jj, vv = find(constant(A))
    solve_ = load_op_and_grad("./build/libSolve","solve")
    rhs,ii, jj, vv,o = convert_to_tensor([rhs,ii, jj, vv,o], [Float64,Int64, Int64, Float64,Int64])
    solve_(rhs,ii, jj, vv,o)
end

function sparse_factorization_solve(A, rhs)
    o = sparse_factorization(A)
    sparse_solve(A, rhs, o)
end

# TODO: specify your input parameters
A = sprand(10,10,0.3)
rhs1 = rand(10)
rhs2 = rand(10)
u = sparse_factorization(A)

out1 = sparse_solve(A, rhs1, u)
out2 = sparse_solve(A, rhs2, u)
sess = Session(); init(sess)
u, o1, o2 = run(sess, [u, out1, out2])
@info u
@show o1 - A\rhs1
@show o2 - A\rhs2

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
ii, jj, vv = find(constant(A))
rhs = rand(10)

function scalar_function(m)
    A = SparseTensor(ii, jj, m, 10, 10)
    return sum(sparse_factorization_solve(A, rhs)^2)
end

# TODO: change `m_` and `v_` to appropriate values
k = length(vv)
m_ = constant(rand(k))
v_ = rand(k)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")

A = sprand(10, 10, 0.8)
ii, jj, vv = find(constant(A))
k = length(vv)
function while_loop_simulation(vv, rhs, ns = 10)
    A = SparseTensor(ii, jj, vv, 10, 10) + spdiag(10)*100.

    o = sparse_factorization(A)

    ta = TensorArray(ns)
    i = constant(2, dtype=Int32)
    ta = write(ta, 1, ones(10))
    function condition(i, ta)
        i<= ns
    end
    function body(i, ta)
        u = read(ta, i-1)
        res = sparse_solve(A, u + rhs, o)
        # res = u
        ta = write(ta, i, res)
        i+1, ta
    end
    _, out = while_loop(condition, body, [i, ta])
    sum(stack(out)^2)
end

vv_ = run(sess, vv)
pl = placeholder(vv_)

# test for vv
pl = placeholder(rand(k))
res = while_loop_simulation(pl, rhs , 100)
gradview(sess, pl, res, rand(k))

# test for right hand side
pl = placeholder(rand(10))
res = while_loop_simulation(vv_, pl , 100)
gradview(sess, pl, res, rand(10))


text:
"""
Describes a river that is diverted to a reservoir. It is akin to a natural
river, with the exception that the actual flow from this source can be decided,
as long as it respects some physical constraints:

  * at least some environmental flow remaining in the river.
  * at most the flow of the river is diverted.
  * at most the capacity of the pipe is diverted.

With respect to the abstract `River`, it defines two new variables to deal with
these new constraints: `environmental_flow` and `maximum_flow`.
"""
immutable DivertedRiver <: River
  name::AbstractString
  scenarios::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1} # 10^6 m^3 / period, indexed with scenarios.
  ext::Dict{AbstractString, Any}

  environmental_flow::Float64 # 10^6 m^3 / s
  maximum_flow::Float64 # 10^6 m^3 / s

  # Private constructor, without units nor defaut values.
  function DivertedRiver(name::AbstractString,
                scenarios::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1},
                ext::Dict{AbstractString, Any},
                environmental_flow::Float64, maximum_flow::Float64)
    o = new(name, scenarios, ext, environmental_flow, maximum_flow)
    _checkConsistency(o)
    return o
  end
end

# Public constructors, with units and keyword arguments.
DivertedRiver(;name::AbstractString=error("Missing reservoir name"),
               scenarios::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1}=error("Missing scenarios"),
               ext::Dict{AbstractString, Any}=Dict{AbstractString, Any}(),
               environmental_flow::Discharge=0m^3/s, maximum_flow::Discharge=(Inf)m^3/s) =
  DivertedRiver(name, scenarios, ext, _from_unitful(environmental_flow), _from_unitful(maximum_flow))

copy(r::DivertedRiver; name::AbstractString=r.name,
                       scenarios::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1}=r.scenarios,
                       ext::Dict{AbstractString, Any}=r.ext,
                       environmental_flow::Union{Discharge, Float64}=r.environmental_flow,
                       maximum_flow::Union{Discharge, Float64}=r.maximum_flow
    ) =
  # Always call the internal constructor!
  DivertedRiver(name, scenarios, ext, _from_unitful(environmental_flow), _from_unitful(maximum_flow))

Base.isempty(r::DivertedRiver) = false

getEnvironmentalFlow(r::DivertedRiver) = r.environmental_flow
getMaximumFlow(r::DivertedRiver) = r.maximum_flow
getEnvironmentalFlow(r::DivertedRiver, p::Period) = r.environmental_flow * period2seconds(p) # m^3/s -> m^3/period
getMaximumFlow(r::DivertedRiver, p::Period) = r.maximum_flow * period2seconds(p) # m^3/s -> m^3/period

function getMaximumAllowableFlow(r::DivertedRiver, scenario::Int, t::Int)
  period = getPeriod(r)

  actual_flow = getScenario(r, scenario, t)
  remaining_flow = getEnvironmentalFlow(r, period)
  maximum_diverted = getMaximumFlow(r, period)

  available_flow = actual_flow - remaining_flow
  available_flow = (available_flow > 0.) ? available_flow : 0.
  return (available_flow > maximum_diverted) ? maximum_diverted : available_flow
end

environmentalFlow(r::DivertedRiver) = getEnvironmentalFlow(r)
maxFlow(r::DivertedRiver) = getMaximumFlow(r)
environmentalFlow(r::DivertedRiver, p::Period) = getEnvironmentalFlow(r, p)
maxFlow(r::DivertedRiver, p::Period) = getMaximumFlow(r, p)
maxAllowableFlow(r::DivertedRiver, scenario::Int, t::Int) = getMaximumAllowableFlow(r, scenario, t)

ext: jl


text:
using LaserTypes, Test
using Unitful
using UnitfulAtomic
import PhysicalConstants.CODATA2018: c_0, m_e, μ_0, e, α

@testset "SI Unitful" begin
    c = c_0
    q = -e
    m = m_e
    λ = 0.8u"μm"
    w0 = 58.0u"μm"
    τ = 18.0u"fs"
    a₀ = 2.0
    p = GaussLaser(:SI_unitful; λ, w₀=w0, a₀=2.0, profile=GaussProfile(;τ,z₀=0.0u"μm"))
    s = setup_laser(GaussLaser, :SI_unitful; a₀=2.0, profile = GaussProfile, τ)
    @test p == s

    @test unit(s.constants.c) == u"m/s"
    @test unit(s.derived.k) == u"μm^-1"
    @test unit(s.derived.z_R) == u"μm"

    # Alternatives
    ω = 2π * c / λ
    a1 = setup_laser(GaussLaser, :SI_unitful; ω, w₀=w0, a₀=2.0, τ)
    @test a1 == s

    k = 2π / λ
    a2 = setup_laser(GaussLaser, :SI_unitful; ω, w₀=w0, a₀=2.0, profile = GaussProfile, τ)
    @test a2 == s

    z_R = w0^2 * k / 2
    a3 = setup_laser(GaussLaser, :SI_unitful; k, z_R, a₀=2.0, profile = GaussProfile(;τ, z₀=0.0u"μm"))
    @test a3 == s
end

@testset "Atomic Unitful" begin
    c = auconvert(c_0)
    q = auconvert(-e)
    m = auconvert(m_e)
    μ₀ = auconvert(μ_0)
    λ = auconvert(800u"nm")
    w0 = auconvert(58u"μm")
    τ = auconvert(18u"fs")
    a₀=2.0

    p = GaussLaser(:atomic_unitful; λ, w₀=w0, a₀=2.0, profile=GaussProfile(;τ,z₀=0.0u"a0_au"))
    s = setup_laser(GaussLaser, :atomic_unitful; a₀=2.0, profile = GaussProfile, τ)
    @test p == s

    @test unit(s.constants.c) == u"a0_au*Eh_au/ħ_au"
    @test unit(s.derived.k) == u"a0_au^-1"
    @test unit(s.derived.z_R) == u"a0_au"
    @test s.derived.ω ≈ 0.05695419u"Eh_au/ħ_au"
end

@testset "SI" begin
    c = ustrip(c_0)
    q = ustrip(-e)
    m = ustrip(m_e)
    μ₀ = ustrip(μ_0)
    λ = ustrip(u"m", 0.8u"μm")
    w0 = ustrip(u"m", 58.0u"μm")
    τ = ustrip(u"s", 18.0u"fs")
    p = GaussLaser(:SI; λ, w₀=w0, a₀=2.0, profile=GaussProfile(;τ, z₀=0.0))
    s = setup_laser(GaussLaser, :SI, a₀=2.0)
    @test p == s

    @test s.derived.ω ≈ 2.354564459e15
end

@testset "Atomic" begin
    c = ustrip(auconvert(c_0))
    q = ustrip(auconvert(-e))
    m = ustrip(auconvert(m_e))
    μ₀ = ustrip(auconvert(μ_0))
    λ = ustrip(auconvert(800u"nm"))
    w0 = ustrip(auconvert(58u"μm"))
    τ = ustrip(auconvert(18u"fs"))
    a₀=2.0

    p = GaussLaser(:atomic; λ, w₀=w0, a₀=2.0, profile=GaussProfile(τ=τ, z₀=0.0))
    s = setup_laser(GaussLaser, :atomic, a₀=2.0)
    @test p == s

    @test s.derived.ω ≈ 0.05695419
end


text:
export LTspiceSimulation
export circuitpath, logpath, executablepath
export parameternames, parametervalues
export measurementnames, measurementvalues
export stepnames, stepvalues
export run!

mutable struct Status
  ismeasurementsdirty :: Bool # true = need to run simulation
  timestamp :: DateTime; # timestamp from last simulation run
  duration :: Float64; # simulation time in seconds
  Status() = new(true,DateTime(1900),NaN)
end
mutable struct StepValues{Nstep}
  values ::NTuple{Nstep,Array{Float64,1}}
end
blankstepvalues(Nstep::Int) = StepValues{Nstep}(ntuple(d->[],Nstep))

"""
Access parameters and measurements of an LTspice simulation.  Runs simulation
as needed.

Access as a dictionary:
```julia
measurement_value = sim["measurement_name"]
parameter_value = sim["parameter_name"]
sim["parameter_name"] = new_parameter_value
```

Access as a function:
```julia
(m1,m2,m3) = sim(p1,p2,p3)  # simulation with three measurements and three parameters
```

Access as arrays or tuples:
```julia
pnames = parameternames(sim)
mnames = measurementnames(sim)
snames = stepnames(sim)
pvalues = parametervalues(sim)
mvalues = measurementvalues(sim)
svalues = stepvalues(sim)
```
"""
struct LTspiceSimulation{Nparam,Nmeas,Nmdim,Nstep}
  circuitpath :: AbstractString
  logpath :: AbstractString
  executablepath :: AbstractString
  circuitfilearray :: Array{AbstractString,1} # text of circuit file
  parameternames :: NTuple{Nparam,AbstractString}
  parametervalues :: ParameterValuesArray{Float64,1} # values in units A,V,W
  parametermultiplier :: NTuple{Nparam,Float64} # units
  parameterindex :: NTuple{Nparam,Int} # index into circuitfilearray
  parameterdict  :: Dict{AbstractString,Int} # index into name, value, multiplier, index arrays
  measurementnames :: NTuple{Nmeas,AbstractString}
  measurementdict :: Dict{AbstractString,Int} # index into measurmentnames
  measurementvalues :: MeasurementValuesArray{Float64,Nmdim}
  stepnames :: NTuple{Nstep,AbstractString}
  stepvalues :: StepValues{Nstep}
  status :: Status
  circuitfileencoding
  logfileencoding
end
NonSteppedSimulation{Nparam,Nmeas} = LTspiceSimulation{Nparam,Nmeas,1,0}

"""
    parametervalues(sim)

Returns an array of the parameters of `sim` in the order they appear in the
circuit file
"""
parametervalues

"""
    parameternames(sim)

Returns an tuple of the parameters names of `sim` in the order they appear in the
circuit file.
"""
parameternames

"""

    circuitpath(sim)

Returns path to the circuit file.

This is the path to the working circuit file.  If `tempdir=ture` was used
or if running under wine, this will not be the path given to the constructor.
"""
circuitpath

"""
    logpath(sim)

Returns path to the log file.
"""
logpath

"""
    executablepath(sim)

Returns path to the LTspice executable
"""
executablepath

"""
    measurementnames(sim)

Returns an tuple of the measurement names of `sim` in the order they appear in the
circuit file.
"""
measurementnames

"""
    stepnames(sim)

Returns an tuple of step names of `sim`.
"""
stepnames

"""
    measurementvalues(sim)

Retruns measurements of `sim` as an a array of Float64
values.

```julia
value = measurementvalues(sim)[inner_step,
                               middle_step,
                               outer_step,
                               measurement_name] # 3 nested steps
```
"""
measurementvalues

"""
    stepvalues(sim)

Returns the steps of `sim` as a tuple of three arrays of
the step values.
"""
stepvalues

circuitpath(x::LTspiceSimulation) = x.circuitpath
logpath(x::LTspiceSimulation) = x.logpath
executablepath(x::LTspiceSimulation) = x.executablepath
parameternames(x::LTspiceSimulation) = x.parameternames
parametervalues(x::LTspiceSimulation) = x.parametervalues
measurementnames(x::LTspiceSimulation) = x.measurementnames
function measurementvalues(x::LTspiceSimulation)
  run!(x)
  x.measurementvalues
end
stepnames(x::LTspiceSimulation) = x.stepnames
function stepvalues(x::LTspiceSimulation)
  run!(x) # step values can be a function of parameters
  x.stepvalues.values
end

LTspiceSimulation(circuitpath::AbstractString;
                  executablepath::AbstractString = defaultltspiceexecutable(),
                  tempdir::Bool = false,
                  librarysearchpaths = []) =
  LTspiceSimulation(circuitpath, executablepath, tempdir, librarysearchpaths)
function LTspiceSimulation(circuitpath::AbstractString,
                           executablepath::AbstractString,
                           istempdir::Bool,
                           librarysearhpaths)
  originalcircuitpath = circuitpath
  if istempdir
    circuitpath = preparetempdir(circuitpath, executablepath)
  end
  @static if Sys.islinux()
    circuitpath = linkintempdirectoryunderwine(circuitpath)
  end
  circuitparsed = parsecircuitfile(originalcircuitpath,
                                   circuitpath,
                                   executablepath,
                                   librarysearhpaths)
  Nparam = length(circuitparsed.parameternames)
  Nmeas = length(circuitparsed.measurementnames)
  Nstep = length(circuitparsed.stepnames)
  Nmdim = Nstep + 1
  parameterdict = Dict{AbstractString,Int}()
  for i in eachindex(circuitparsed.parameternames)
    parameterdict[circuitparsed.parameternames[i]] = i
  end
  measurementdict = Dict{AbstractString,Int}()
  for i in eachindex(circuitparsed.measurementnames)
    measurementdict[circuitparsed.measurementnames[i]] = i
  end
  return LTspiceSimulation{Nparam,Nmeas,Nmdim,Nstep}(
    circuitpath,
    logpath(circuitpath),
    executablepath,
    circuitparsed.circuitfilearray,
    (circuitparsed.parameternames...,),
    circuitparsed.parametervalues,
    (circuitparsed.parametermultiplier...,),
    (circuitparsed.parameterindex...,),
    parameterdict,
    (circuitparsed.measurementnames...,),
    measurementdict,
    MeasurementValuesArray{Float64,Nmdim}(fill(NaN,(ntuple(d->1,Nstep)...,Nmeas))), # measurementvalues
    (circuitparsed.stepnames...,),
    blankstepvalues(Nstep),
    Status(),
    circuitparsed.circuitfileencoding,
    PossibleEncodings([enc"UTF-16LE",enc"UTF-8",enc"windows-1252"],iscorrectencoding_logfile) # logfileencoding(executablepath) # LTspice changed?
  )
end

function preparetempdir(circuitpath::AbstractString, executablepath::AbstractString)
  td = mktempdir()
  atexit(()->ispath(td)&&rm(td,recursive=true)) # delete this on exit
  (d,f) = splitdir(circuitpath)
  workingcircuitpath = convert(AbstractString, joinpath(td,f))
  cp(circuitpath,workingcircuitpath)
  return workingcircuitpath
end

function Base.show(io::IO, x::LTspiceSimulation)
  println(io,"LTspiceSimulation:")
  showcircuitpath(io,x)
  showparameters(io,x)
  showmeasurements(io,x)
  showsteps(io,x)
  showtimeduration(io,x)
end

function showcircuitpath(io::IO, x::LTspiceSimulation)
  println(io,"circuit path = ",x.circuitpath)
end
showparameters(::IO, x::LTspiceSimulation{0,Nmeas,Nmdim,Nstep}) where {Nmeas,Nmdim,Nstep} = nothing
function showparameters(io::IO, x::LTspiceSimulation)
  println(io)
  println(io,"Parameters")
  for i in eachindex(x.parameternames)
    println(io,rpad(x.parameternames[i],25,' ')," = ",x.parametervalues[i])
  end
end
showmeasurments(::IO, x::LTspiceSimulation{0,0,1,0}) = nothing
showmeasurments(::IO, x::LTspiceSimulation{Nparam,0,1,0}) where Nparam= nothing
showmeasurments(::IO, x::LTspiceSimulation{Nparam,0,Nmdim,Nstep}) where {Nparam,Nmdim,Nstep} = nothing
function showmeasurements(io::IO, x::LTspiceSimulation{Nparam,Nmeas,1,0}) where {Nparam,Nmeas}
  if Nmeas!=0
    println(io)
    println(io,"Measurements")
    for i in eachindex(x.measurementnames)
      print(io,rpad(x.measurementnames[i],25,' '))
      if x.status.ismeasurementsdirty
        println(io)
      elseif isnan(x.measurementvalues[i])
        println(io," = measurement failed")
      else
        println(io," = ",x.measurementvalues[i])
      end
    end
  end
end
function showmeasurements(io::IO, x::LTspiceSimulation)
  println(io)
  println(io,"Measurements")
  totalsteps = 1
  s = size(x.measurementvalues)
  for i in 1:length(s)-1
    totalsteps *= s[i]
  end
  for i in eachindex(x.measurementnames)
    print(io,rpad(x.measurementnames[i],25,' '))
    if x.status.ismeasurementsdirty
      println(io)
    else
      println(io," ",totalsteps," values")
    end
  end
end
showsteps(::IO, x::LTspiceSimulation{Nparam,Nmeas,1,0}) where {Nparam,Nmeas}= nothing
function showsteps(io::IO, x::LTspiceSimulation)
  if length(x.stepnames)>0
    println(io)
    println(io,"Steps")
    for i in eachindex(x.stepnames)
      print(io,rpad(x.stepnames[i],25,' '))
      if x.status.ismeasurementsdirty
        println(io)
      else
        println(io," ",length(x.stepvalues.values[i])," steps")
      end
    end
  end
end
function showtimeduration(io::IO, x::LTspiceSimulation)
  if ~isnan(x.status.duration) # simulation was run at least once
    println(io)
    println(io,"Last Run")
    println(io,"time = ",x.status.timestamp)
    println(io,"duration = ",x.status.duration)
  end
end

Base.haskey(x::LTspiceSimulation, key::AbstractString) =
  haskey(x.parameterdict,key) || haskey(x.measurementdict,key)
Base.keys(x::LTspiceSimulation) =
  vcat(x.parameternames...,x.measurementnames...)
function Base.values(x::LTspiceSimulation)
  run!(x)
  allkeys = keys(x)
  allvalues = similar(allkeys, Float64)
  for i in eachindex(allkeys)
    allvalues[i] =  x[allkeys[i]]
  end
  return allvalues
end
@generated function Base.getindex(
        x::LTspiceSimulation{Nparam,Nmeas,Nmdim,Nstep},
        key::AbstractString) where {Nparam,Nmeas,Nmdim,Nstep}
  r = (
      :(x.measurementvalues[x.measurementdict[key]]),
      :(x.measurementvalues[:,x.measurementdict[key]]),
      :(x.measurementvalues[:,:,x.measurementdict[key]]),
      :(x.measurementvalues[:,:,:,x.measurementdict[key]])
      )
  return quote
    if haskey(x.parameterdict,key)
      return x.parametervalues[x.parameterdict[key]]
    elseif haskey(x.measurementdict,key)
      run!(x)
      return $(r[Nstep+1])
    else
      throw(KeyError(key))
    end
  end
end
function Base.get(x::LTspiceSimulation, key::AbstractString, default)
  if haskey(x,key)
    return(x[key])
  else
    return(default)
  end
end
function Base.setindex!(x::LTspiceSimulation, value::Float64, key::AbstractString)
  if haskey(x.parameterdict,key)
    x.parametervalues[x.parameterdict[key]] = value
  elseif haskey(x.measurementdict,key)
    error("measurements cannot be set.")
  else
    throw(KeyError(key))
  end
end
Base.eltype(x::LTspiceSimulation) = Float64
Base.length(x::LTspiceSimulation) = length(x.parametervalues) + length(x.measurementvalues)

(x::LTspiceSimulation)(args...) = throw(ArgumentError("number of arguments must match number of parameters"))
function (x::LTspiceSimulation{Nparam,Nmeas,Nmdim,Nstep})(args::Vararg{Any,Nparam}) where {Nparam,Nmeas,Nmdim,Nstep}
  for i in eachindex(args)
    x.parametervalues[i] = args[i]
  end
  measurementvalues(x)
end

"""
```julia
flush(sim)
```
Writes `sim`'s circuit file back to disk if any parameters have changed.  The
user does not usually need to call `flush`.  It will be called automatically
 when a measurement is requested and the log file needs to be updated.  It can be used
 to update a circuit file using julia for simulation with the LTspice GUI.
"""
function Base.flush(x::LTspiceSimulation, force=false)
	if x.parametervalues.ismodified || force
    updatecircuitfilearray!(x)
    writecircuitfilearray(x)
  	x.parametervalues.ismodified = false
    x.status.ismeasurementsdirty = true
  end
  return nothing
end
function updatecircuitfilearray!(x::LTspiceSimulation)
  for i in eachindex(x.parameternames)
    x.circuitfilearray[x.parameterindex[i]] =
      string(x.parametervalues[i]/x.parametermultiplier[i])
  end
end
function writecircuitfilearray(x::LTspiceSimulation)
  io = open(x.circuitpath, x.circuitfileencoding, "w")
  for text in x.circuitfilearray
    print(io,text)
  end
  close(io)
  return nothing
end

"""
```julia
run!(sim)
```
Writes circuit changes, calls LTspice to run `sim`, and reloads the log file.  The user
normally does not need to call this.
"""
function run!(x::LTspiceSimulation, force=false)
  flush(x,force)
  if x.status.ismeasurementsdirty || force
    if x.executablepath != ""  # so travis dosen't need to load LTspice
      @static if Sys.islinux()
        drive_c = "/home/$(ENV["USER"])/.wine/drive_c"
        winecircuitpath = joinpath("C:",relpath(x.circuitpath,drive_c))
        run(`wine $(x.executablepath) -b -Run $winecircuitpath`)
      else
        run(`$(x.executablepath) -b -Run $(x.circuitpath)`)
      end
    end
    parselog!(x)
    x.status.ismeasurementsdirty = false
  end
end


text:
int_rules_4_1_4_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.4.1*(a+b*sin)^m*(A+B*sin+C*sin^2) =#
    @apply_utils Antiderivative((~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(B') * sin(~(e') + ~(f') * ~x) + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (1 / ~b) * Antiderivative((~b * sin(~e + ~f * ~x)) ^ (~m + 1) * (~B + ~C * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~b, ~e, ~f, ~B, ~C, ~m], ~x)
    @apply_utils Antiderivative((~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~A + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (~A * cos(~e + ~f * ~x) * (~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * (~m + 1)) <-- FreeQ([~b, ~e, ~f, ~A, ~C, ~m], ~x) && EqQ(~A * (~m + 2) + ~C * (~m + 1), 0)
    @apply_utils Antiderivative((~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~A + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (~A * cos(~e + ~f * ~x) * (~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * (~m + 1)) + ((~A * (~m + 2) + ~C * (~m + 1)) / ((~b) ^ 2 * (~m + 1))) * Antiderivative((~b * sin(~e + ~f * ~x)) ^ (~m + 2), ~x) <-- FreeQ([~b, ~e, ~f, ~A, ~C], ~x) && LtQ(~m, -1)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(m') * (~A + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (-1 / ~f) * Subst(Antiderivative((1 - (~x) ^ 2) ^ ((~m - 1) / 2) * ((~A + ~C) - ~C * (~x) ^ 2), ~x), ~x, cos(~e + ~f * ~x)) <-- FreeQ([~e, ~f, ~A, ~C], ~x) && IGtQ((~m + 1) / 2, 0)
    @apply_utils Antiderivative((~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~A + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(~C) * cos(~e + ~f * ~x) * (~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * (~m + 2)) + ((~A * (~m + 2) + ~C * (~m + 1)) / (~m + 2)) * Antiderivative((~b * sin(~e + ~f * ~x)) ^ ~m, ~x) <-- FreeQ([~b, ~e, ~f, ~A, ~C, ~m], ~x) && Not(LtQ(~m, -1))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x) + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (1 / (~b) ^ 2) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * Simp((~b * ~B - ~a * ~C) + ~b * ~C * sin(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C, ~m], ~x) && EqQ((~A * (~b) ^ 2 - ~a * ~b * ~B) + (~a) ^ 2 * ~C, 0)
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(A') + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (~C / (~b) ^ 2) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * Simp(-(~a) + ~b * sin(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C, ~m], ~x) && EqQ(~A * (~b) ^ 2 + (~a) ^ 2 * ~C, 0)
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x) + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (~A - ~C) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (1 + sin(~e + ~f * ~x)), ~x) + ~C * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (1 + sin(~e + ~f * ~x)) ^ 2, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C, ~m], ~x) && (EqQ((~A - ~B) + ~C, 0) && Not(IntegerQ(2 * ~m)))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(A') + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (~A - ~C) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (1 + sin(~e + ~f * ~x)), ~x) + ~C * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (1 + sin(~e + ~f * ~x)) ^ 2, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C, ~m], ~x) && (EqQ(~A + ~C, 0) && Not(IntegerQ(2 * ~m)))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x) + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (((~A * ~b - ~a * ~B) + ~b * ~C) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m) / (~a * ~f * (2 * ~m + 1)) + (1 / ((~a) ^ 2 * (2 * ~m + 1))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * Simp(~a * ~A * (~m + 1) + ~m * (~b * ~B - ~a * ~C) + ~b * ~C * (2 * ~m + 1) * sin(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C], ~x) && (LtQ(~m, -1) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (~b * (~A + ~C) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m) / (~a * ~f * (2 * ~m + 1)) + (1 / ((~a) ^ 2 * (2 * ~m + 1))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * Simp((~a * ~A * (~m + 1) - ~a * ~C * ~m) + ~b * ~C * (2 * ~m + 1) * sin(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C], ~x) && (LtQ(~m, -1) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x) + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(((~A * (~b) ^ 2 - ~a * ~b * ~B) + (~a) ^ 2 * ~C)) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * (~m + 1) * ((~a) ^ 2 - (~b) ^ 2)) + (1 / (~b * (~m + 1) * ((~a) ^ 2 - (~b) ^ 2))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * Simp(~b * ((~a * ~A - ~b * ~B) + ~a * ~C) * (~m + 1) - ((~A * (~b) ^ 2 - ~a * ~b * ~B) + (~a) ^ 2 * ~C + ~b * ((~A * ~b - ~a * ~B) + ~b * ~C) * (~m + 1)) * sin(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C], ~x) && (LtQ(~m, -1) && NeQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (-((~A * (~b) ^ 2 + (~a) ^ 2 * ~C)) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * (~m + 1) * ((~a) ^ 2 - (~b) ^ 2)) + (1 / (~b * (~m + 1) * ((~a) ^ 2 - (~b) ^ 2))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * Simp(~a * ~b * (~A + ~C) * (~m + 1) - (~A * (~b) ^ 2 + (~a) ^ 2 * ~C + (~b) ^ 2 * (~A + ~C) * (~m + 1)) * sin(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C], ~x) && (LtQ(~m, -1) && NeQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x) + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(~C) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * (~m + 2)) + (1 / (~b * (~m + 2))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * Simp(~A * ~b * (~m + 2) + ~b * ~C * (~m + 1) + (~b * ~B * (~m + 2) - ~a * ~C) * sin(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B, ~C, ~m], ~x) && Not(LtQ(~m, -1))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(A') + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(~C) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * (~m + 2)) + (1 / (~b * (~m + 2))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * Simp((~A * ~b * (~m + 2) + ~b * ~C * (~m + 1)) - ~a * ~C * sin(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~C, ~m], ~x) && Not(LtQ(~m, -1))
    @apply_utils Antiderivative((~(b') * sin(~(e') + ~(f') * ~x) ^ ~p) ^ ~m * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x) + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => ((~b * sin(~e + ~f * ~x) ^ ~p) ^ ~m / (~b * sin(~e + ~f * ~x)) ^ (~m * ~p)) * Antiderivative((~b * sin(~e + ~f * ~x)) ^ (~m * ~p) * (~A + ~B * sin(~e + ~f * ~x) + ~C * sin(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~b, ~e, ~f, ~A, ~B, ~C, ~m, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(b') * cos(~(e') + ~(f') * ~x) ^ ~p) ^ ~m * (~(A') + ~(B') * cos(~(e') + ~(f') * ~x) + ~(C') * cos(~(e') + ~(f') * ~x) ^ 2), ~x) => ((~b * cos(~e + ~f * ~x) ^ ~p) ^ ~m / (~b * cos(~e + ~f * ~x)) ^ (~m * ~p)) * Antiderivative((~b * cos(~e + ~f * ~x)) ^ (~m * ~p) * (~A + ~B * cos(~e + ~f * ~x) + ~C * cos(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~b, ~e, ~f, ~A, ~B, ~C, ~m, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(b') * sin(~(e') + ~(f') * ~x) ^ ~p) ^ ~m * (~(A') + ~(C') * sin(~(e') + ~(f') * ~x) ^ 2), ~x) => ((~b * sin(~e + ~f * ~x) ^ ~p) ^ ~m / (~b * sin(~e + ~f * ~x)) ^ (~m * ~p)) * Antiderivative((~b * sin(~e + ~f * ~x)) ^ (~m * ~p) * (~A + ~C * sin(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~b, ~e, ~f, ~A, ~C, ~m, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(b') * cos(~(e') + ~(f') * ~x) ^ ~p) ^ ~m * (~(A') + ~(C') * cos(~(e') + ~(f') * ~x) ^ 2), ~x) => ((~b * cos(~e + ~f * ~x) ^ ~p) ^ ~m / (~b * cos(~e + ~f * ~x)) ^ (~m * ~p)) * Antiderivative((~b * cos(~e + ~f * ~x)) ^ (~m * ~p) * (~A + ~C * cos(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~b, ~e, ~f, ~A, ~C, ~m, ~p], ~x) && Not(IntegerQ(~m))
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.4.1 (a+b sin)^m (A+B sin+C sin^2).jl


text:
function get_test_spectrum(seed)
    seed!(seed)
    x = collect(0:0.1:200)
    y = @. exp(-(x-15)^2) + 2 * exp(-(x-30)^2)
    @. y += 1.0 + x*1.5e-2 - (x-50)^3*3e-7 # add polynomial baseline
    n = length(x)
    δx = x[2] - x[1]
    return Curve(x, y .+ (get_cov(δx, n, 0.1, 0.5) |> MvNormal |> rand))
end

@testset "regression tests" begin
    @testset "integration of single peak" begin
        c = crop(get_test_spectrum(1), 10, 40)
        uc = add_noise(c, MvGaussianNoiseModel(0.1, 0.1, 0.5))
        ub = UncertainBound(15.0, scale_shift_beta(2.0, 2.0, 3.0, 4.0), uc)
        area = mc_integrate(uc, ub, subtract_baseline=true)
        @test mean(area) ≈ 1.808961690 atol = 1e-8
        @test std(area) ≈ 0.2488841198 atol = 1e-8
    end
    @testset "integration of two peaks" begin
        c = crop(get_test_spectrum(1), 10, 40)
        uc = add_noise(c, MvGaussianNoiseModel(0.1, 0.1, 0.5))
        ubs = UncertainBound([15., 30.], scale_shift_beta(2.0, 2.0, 3.0, 4.0), uc)
        area1, area2 = mc_integrate(uc, ubs, subtract_baseline=true)
        @test mean(area1) ≈ 1.808961690 atol = 1e-8
        @test std(area1) ≈ 0.2488841198 atol = 1e-8
        @test mean(area2) ≈ 3.442573414 atol = 1e-8
        @test std(area2) ≈ 0.3122904192 atol = 1e-8
    end
end

@testset "pass custom integration function" begin
    f(args...; kw...) = 0.0 # this should make integration yield area = 0 ± 0
    c = crop(get_test_spectrum(1), 10, 40)
    uc = add_noise(c, MvGaussianNoiseModel(0.1, 0.1, 0.5))
    ub = UncertainBound(15.0, scale_shift_beta(2.0, 2.0, 3.0, 4.0), uc)
    area = mc_integrate(uc, ub; intfun=f)
    @test mean(area) == 0.0
    @test std(area) == 0.0
end

@testset "integrate using Simpsons rule from external package" begin

    using NumericalIntegration: integrate, SimpsonEven

    function simps(x, y, a, b; subtract_baseline=true)
        idx = [(xᵢ >= a) & (xᵢ <= b) for xᵢ in x]
        return integrate(x[idx], y[idx], SimpsonEven())
    end
    c = crop(get_test_spectrum(1), 10, 40)
    uc = add_noise(c, MvGaussianNoiseModel(0.1, 0.1, 0.5))
    ub = UncertainBound(15.0, scale_shift_beta(2.0, 2.0, 3.0, 4.0), uc)
    area_s = mc_integrate(uc, ub; intfun=simps)
    @test mean(area_s) ≈ 5.6868358 atol = 1e-7
    @test std(area_s) ≈ 0.32897085 atol = 1e-7
end

@testset "error if local_baseline and subtract_baseline true" begin
    c = NoisySignalIntegration.testdata_1()
    uc = add_noise(c, MvGaussianNoiseModel(0.1, 0.1, 0.5))
    ub = UncertainBound(15.0, scale_shift_beta(2.0, 2.0, 3.0, 4.0), uc)
    @test_throws ErrorException area = mc_integrate(uc, ub; local_baseline=true, subtract_baseline=true)
end

@testset "test _local_baseline" begin
    xs = 1.0:10.0 |> collect
    b = UncertainBound(Uniform(1.5, 2.5), Uniform(4.5, 6.5))
    @testset "test _local_baseline with linear dataset" begin
        ys = xs
        xₗ, xᵣ, yₗ, yᵣ = nsi._local_baseline(xs, ys, 1.75, 5.25, b)
        yyₗ = nsi.lininterp(1.75, xs, ys)
        yyᵣ = nsi.lininterp(5.25, xs, ys)
        @test yyₗ == yₗ
        @test yyᵣ == yᵣ
    end
    @testset "test _local_baseline with parabolic dataset" begin
        ys = xs.^2
        xₗ, xᵣ, yₗ, yᵣ = nsi._local_baseline(xs, ys, 2.25, 5.25, b)
        @test yₗ ≈ 6.1329 atol=1e-4  # results derived by doing computation by hand in Numpy
        @test yᵣ ≈ 28.8471 atol=1e-4
    end
    @testset "test _local_baseline with normal and uniform distribution on noise" begin
        noise = [
            1 1.51
            2 1.41
            3 1.2
            4 0.88
            5 0.5
            6 0.09
            7 -0.32
            8 -0.67
            9 -0.94
            10 -1.09
            11 -1.12
            12 -1.04
            13 -0.86
            14 -0.63
            15 -0.39
            16 -0.18
            17 -0.04
            18 0.02
            19 -0.01
            20 -0.12
            21 -0.25
            22 -0.36
            23 -0.43
            24 -0.43
            25 -0.4
            26 -0.36
            27 -0.35
            28 -0.39
            29 -0.48
            30 -0.6
            31 -0.73
            32 -0.8
            33 -0.8
            34 -0.71
            35 -0.53
            36 -0.3
            37 -0.06
            38 0.15
            39 0.27
            40 0.29
        ]
        bnd = UncertainBound(Normal(4.0, 2.0), Uniform(30.0, 35.0))
        xₗ, xᵣ, yₗ, yᵣ = nsi._local_baseline(noise[:,1], noise[:,2], 4.0, 32.5, bnd)
        @test yₗ ≈ 0.77916842211701 # manually derived in LibreCalc
        @test yᵣ ≈ -0.7210156
    end
end


text:
function rms(A)
   s = 0.0
   for a in A
      s += a*a
   end
   return sqrt(s / length(A))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function data33_14()
    row = [1, 9, 11, 18, 2, 9, 10, 11, 19, 20, 30, 31, 9, 10, 21, 33, 9, 10, 11, 20, 21, 27, 9, 11, 12, 18, 19, 26, 3, 11, 12, 14, 15, 16, 24, 26,
           4, 13, 27, 28, 29, 5, 13, 28, 17, 22, 29, 14, 23, 12, 14, 23, 6, 12, 15, 16, 24, 7, 12, 15, 16, 17, 25, 8, 16, 17, 22, 25, 32]
    col = [1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 9,
           9, 9, 10, 10, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14]
    cof = [4.0, -16.90045631, -4.483500717, 4.483500717, 14.0, 33.37432577, -5.051270395, -5.751092708, 5.751092708, 5.671506352, 10.0, 30.0,
           -5.051270395, 10.89819783, 5.84692744, 1.0, -5.671506352, -5.84692744, -23.74732843, -5.671506352, -5.84692744, 4.88951266, -5.751092708,
           38.23936719, -4.257445335, -4.483500717, -5.751092708, 4.257445335, 1.0, -4.257445335, 20.87061322, -5.027652086, -3.909151323, -7.676364474,
           3.909151323, -4.257445335, 1.0, -5.676979847, -4.88951266, 5.676979847, 9.09008272, 1.0, 5.676979847, -5.676979847, -3.69849841, 3.69849841,
           -9.09008272, -5.206435154, 5.206435154, -5.027652086, 10.23408724, -5.206435154, 1.0, -3.909151323, 8.912153124, -5.003001801, -3.909151323, 1.0,
           -7.676364474, -5.003001801, 15.55276436, -2.873398081, 2.873398081, 1.0, -2.873398081, 6.57189649, -3.69849841, -2.873398081, 6.0]
    H = sparse(row, col, cof)

    z = [0.0, -0.087476928, -0.259210842, -0.242730471, -0.242719726, -0.278676774, -0.281695364, -0.299976206, 0.183367375, -0.94065623, -0.077346635,
     -0.112748813, -0.000215539, -0.0356358, -0.061994423, -0.134358588, -0.148770477, 0.712008594, 0.408969228, 0.552256282, -0.242603781, 0.096901144,
     -0.032970493, 0.076149458, 0.053230571, 0.427620284, 0.284565166, 4.68e-5, 0.283247793, 4.0, 8.0, 2.0, 0.0]

    v = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

    return H, z, v
end

function discrete6_4()
       probability1 = [1, 2]
       table1 = [0.0 0.8; 0.1 0.7; 0.3 0.2]

       probability2 = [3]
       table2 = [0.2, 0.1, 0.6]

       probability3 = [2, 3, 4, 5]
       table3 = zeros(2, 3, 2, 3)
       table3[1, 1, 1, 1] =	0.2
       table3[2, 1, 1, 1] =	0.1
       table3[1, 2, 1, 1] = 0.5
       table3[2, 2, 1, 1] = 0.6
       table3[1, 3, 1, 1] = 0.2
       table3[2, 3, 1, 1] = 0.4
       table3[1, 1, 2, 1] = 0.4
       table3[2, 1, 2, 1] = 0.2
       table3[1, 2, 2, 1] = 0.9
       table3[2, 2, 2, 1] = 0.4
       table3[1, 3, 2, 1] = 0.9
       table3[2, 3, 2, 1] = 0.4
       table3[1, 1, 1, 2] = 0.8
       table3[2, 1, 1, 2] = 0.2
       table3[1, 2, 1, 2] = 0.7
       table3[2, 2, 1, 2] = 0.1
       table3[1, 3, 1, 2] = 0.6
       table3[2, 3, 1, 2] = 0.5
       table3[1, 1, 2, 2] = 0.1
       table3[2, 1, 2, 2] = 0.1
       table3[1, 2, 2, 2] = 0.5
       table3[2, 2, 2, 2] = 0.5
       table3[1, 3, 2, 2] = 0.7
       table3[2, 3, 2, 2] = 0.6
       table3[1, 1, 1, 3] = 0.1
       table3[2, 1, 1, 3] = 0.4
       table3[1, 2, 1, 3] = 0.1
       table3[2, 2, 1, 3] = 0.3
       table3[1, 3, 1, 3] = 0.1
       table3[2, 3, 1, 3] = 0.2
       table3[1, 1, 2, 3] = 0.2
       table3[2, 1, 2, 3] = 0.1
       table3[1, 2, 2, 3] = 0.2
       table3[2, 2, 2, 3] = 0.8
       table3[1, 3, 2, 3] = 0.9
       table3[2, 3, 2, 3] = 0.9

       probability4 = [4]
       table4 = [0.6, 0.2]

       probability5 = [6, 5]
       table5 = [0.5 0.3 0.2; 0.1 0.7 0.1; 0.8 0.1 0.8; 0.3 0.4 0.6]

       probability6 = [5]
       table6 = [0.9, 0.5, 0.1]

       probability = [probability1, probability2, probability3, probability4, probability5, probability6]
       table = [table1, table2, table3, table4, table5, table6]

       return probability, table
end


text:
@deprecate error(obs::Observable) std_error(obs::Observable)
export error

@deprecate binning_error(ts) std_error(ts, method=:full)
export binning_error

@deprecate jackknife_error(g, obs...) jackknife(g, obs...)[2]

@deprecate rename(obs::Observable, name) rename!(obs::Observable, name)
export rename

@deprecate add!(obs::Observable, measurement; kw...) push!(obs::Observable, measurement; kw...)


text:
import Profile

import SpecialFunctions # for triggering invalidations
import Static # for triggering invalidations

x = rand(1000)

function foo()
    s = 0.0
    for i in x
        s += i - rand()
    end
    return s
end

for i in 1:100
    foo()
end
Profile.clear_malloc_data()
for i in 1:100
    foo()
end


text:
Base.@kwdef struct ModelParameters{A<:AbstractArray{<:AbstractFloat}, B<:AbstractArray{<:Bool}}
    Δt::A
    Tconn::A
    Pi::A
    Bwi::A
    Boi::A
    cw::A
    co::A
    cf::A
    Swi::A
    Vpi::A
    Tconst::A
    Gw::A
    Jp::A
    M::A
    Jinj::A
    λ::A
    Pmin::A
    Pmax::A
    Qliq_h::A
    Qinj_h::A
    Pcalc::A
    Qliq::A
    Qinj::A
    jac_next::A
    Pbhp::A
    Pinj::A
    Tupd::B
    Vupd::B
    cupd::B
end

Base.@kwdef struct ProblemCache{T<:AbstractFloat}
    Vw::Vector{T}
    Vo::Vector{T}
    Vwprev::Vector{T}
    Voprev::Vector{T}
    Vwi::Vector{T}
    Voi::Vector{T}
    cwf::Vector{T}
    cof::Vector{T}
    Qsum::Vector{T}
    CTC::SparseMatrixCSC{T, Int}
    Cbuf::SparseMatrixCSC{T, Int}
    diagJ::Vector{T}
end

Base.@kwdef struct NonlinearProblem{T<:AbstractFloat}
    params::ModelParameters{Matrix{T}, BitMatrix}
    pviews::Vector{ModelParameters{ColumnSlice{T}, ColumnSliceBool}}
    C::SparseMatrixCSC{T, Int}
    cache::ProblemCache{T}
end

function build_incidence_matrix(::Type{T}, df) where {T}

    # Tanks numbering
    dict = @with df begin
        @_ :Tank::Vector{String} |>
            unique(__) |>
            Dict(__ .=> 1:length(__))
    end

    # Numbers of neighbouring tanks
    df_view = getparams(df, Val(:Tconn))
    N = @with df_view begin
        @_ [:Tank :Neighb]::Matrix{String} |>
            unique(__, dims=1) |>
            getindex.(Ref(dict), __)
    end

    # Tank adjacency matrix
    C = spzeros(T, size(N, 1), length(dict))
    C[CartesianIndex.(axes(C, 1), N)] .= [-1 1]

    return C
end

function build_parameter_matrix(::Type{T}, df, name::Symbol, N, M) where {T}
    A = Array{T}(undef, N, M)
    df_view = getparams(df, Val(name))
    @with df_view @byrow begin
        cols = UnitRange(:Jstart::Int, :Jstop::Int)
        A[:Istart::Int, cols] .= convert(T, :Init_value)::T
    end
    return A
end

function build_update_matrix(::Type{T}, df, names::NTuple{N, Symbol}, M) where {T, N}
    A = falses(1, M)
    for name ∈ names
        df_view = getparams(df, Val(name))
        A[df_view.Jstart] .= true
    end
    return A
end

function build_rate_matrix(::Type{T}, data, N, M) where {T}
    @_ data |>
        convert(Vector{T}, __)::Vector{T} |>
        reshape(__, M, N) |>
        permutedims(__)
end

function params_and_views(::Type{T}, data) where {T}
    params = T(; data...)
    pviews = map(zip(eachcol.(values(data))...)) do x
        T(; (keys(data) .=> x)...)
    end
    return params, pviews
end

function NonlinearProblem{T}(df_rates::AbstractDataFrame, df_params::AbstractDataFrame) where {T}

    # Adjacency matrix building
    C = build_incidence_matrix(T, df_params)
    # Number of connections and tanks
    Nc, Nt = size(C)
    # Number of time steps
    dates = unique(df_rates.Date::Vector{Date})
    Nd = length(dates)

    # Form matrices of parameters and history of production/injection
    kwargs = (
        # Input variables of model
        Δt = build_rate_matrix(T, daysinmonth.(dates), 1, Nd),
        Tconn = build_parameter_matrix(T, df_params, :Tconn, Nc, Nd),
        Pi = build_parameter_matrix(T, df_params, :Pi, Nt, Nd),
        Bwi = build_parameter_matrix(T, df_params, :Bwi, Nt, Nd),
        Boi = build_parameter_matrix(T, df_params, :Boi, Nt, Nd),
        cw = build_parameter_matrix(T, df_params, :cw, Nt, Nd),
        co = build_parameter_matrix(T, df_params, :co, Nt, Nd),
        cf = build_parameter_matrix(T, df_params, :cf, Nt, Nd),
        Swi = build_parameter_matrix(T, df_params, :Swi, Nt, Nd),
        Vpi = build_parameter_matrix(T, df_params, :Vpi, Nt, Nd),
        Tconst = build_parameter_matrix(T, df_params, :Tconst, Nt, Nd),
        λ = build_parameter_matrix(T, df_params, :λ, Nt, Nd),
        Pmin = build_parameter_matrix(T, df_params, :Pmin, Nt, Nd),
        Pmax = build_parameter_matrix(T, df_params, :Pmax, Nt, Nd),
        Qliq_h = build_rate_matrix(T, df_rates.Qliq, Nt, Nd),
        Qinj_h = build_rate_matrix(T, df_rates.Qinj, Nt, Nd),
        M = build_rate_matrix(T, df_rates.Total_mobility, Nt, Nd),

        # Auxiliary buffer update flags
        Tupd = build_update_matrix(T, df_params, (:Tconn,), Nd),
        Vupd = build_update_matrix(T, df_params, (:Vpi, :Bwi, :Boi, :Swi,), Nd),
        cupd = build_update_matrix(T, df_params, (:cw, :co, :cf,), Nd),

        # Output variables of model
        Pcalc = build_parameter_matrix(T, df_params, :Pi, Nt, Nd),
        Qliq = build_rate_matrix(T, df_rates.Qliq, Nt, Nd),
        Qinj = build_rate_matrix(T, df_rates.Qinj, Nt, Nd),
        Gw = build_parameter_matrix(T, df_params, :Gw, Nt, Nd),
        Jp = build_parameter_matrix(T, df_params, :Jp, Nt, Nd),
        Jinj = build_parameter_matrix(T, df_params, :Jinj, Nt, Nd),
        Pbhp = build_parameter_matrix(T, df_params, :Pi, Nt, Nd),
        Pinj = build_parameter_matrix(T, df_params, :Pi, Nt, Nd),
        jac_next = Array{T}(undef, Nt, Nd),
    )

    params, pviews = params_and_views(ModelParameters, kwargs)
    cache = ProblemCache{T}(Nt, Nc)
    NonlinearProblem{T}(; params, C, pviews, cache)
end

function ProblemCache{T}(Nt, Nc) where {T}
    kwargs = (
        Vw = Array{T}(undef, Nt),
        Vo = Array{T}(undef, Nt),
        Vwprev = Array{T}(undef, Nt),
        Voprev = Array{T}(undef, Nt),
        Vwi = Array{T}(undef, Nt),
        Voi = Array{T}(undef, Nt),
        cwf = Array{T}(undef, Nt),
        cof = Array{T}(undef, Nt),
        Qsum = Array{T}(undef, Nt),
        CTC = spzeros(T, Nt, Nt),
        Cbuf = spzeros(T, Nc, Nt),
        diagJ = Array{T}(undef, Nt),
    )
    ProblemCache{T}(; kwargs...)
end

function val_and_jac!(r, J, P, Δt, prob::NonlinearProblem, n)

    @unpack Vwprev, Voprev, Vw, Vo, CTC = prob.cache
    @unpack Vwi, Voi, cwf, cof, Qsum, diagJ = prob.cache
    @unpack Pi, Tconst = @inbounds prob.pviews[n]

    # Residual and Jacobian initialization
    mul!(r, CTC, P)
    copyto!(J, CTC)

    @turbo for i = 1:length(P)
        # Pore volumes of water and oil in surface conditions
        Vw[i] = Vwi[i] * exp(cwf[i] * (P[i] - Pi[i]))
        Vo[i] = Voi[i] * exp(cof[i] * (P[i] - Pi[i]))

        # Update the residual vector until the condition is met ∥r∥ ≈ 0
        r[i] += (Vw[i] - Vwprev[i] + Vo[i] - Voprev[i]) / Δt
        r[i] += Tconst[i] * (P[i] - Pi[i]) + Qsum[i]

        # In general, the Jacobian is static except for the diagonal elements
        diagJ[i] = Tconst[i] + (Vw[i] * cwf[i] + Vo[i] * cof[i]) / Δt
    end

    # Update the main diagonal of the Jacobian
    # TODO: We took it out from previous loop, because macro '@turbo' does not work with 'Sparse Arrays'
    @inbounds @simd for i = 1:length(diagJ)
        J[i, i] += diagJ[i]
    end

    return r, J
end


text:
using BinaryBuilder, Pkg

name = "BVLS"
version = v"1990.03.19"
sources = [
    ArchiveSource("https://cran.r-project.org/src/contrib/bvls_1.4.tar.gz", "a22e1951d280c7d281ffb6d4440938dad83b59a810739e8a1a7ad632e7714aa1")
]

script = raw"""
cd ${WORKSPACE}/srcdir/bvls
mkdir -p "${libdir}"
gfortran -shared -fPIC -o $libdir/libbvls.${dlext} src/bvls.f
install_license inst/COPYRIGHTS
sed -n 10,19p src/bvls.f > LICENSE
install_license LICENSE
"""

platforms = supported_platforms()
platforms = expand_gfortran_versions(platforms)

products = [
    LibraryProduct("libbvls", :libbvls),
]

dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll",
        uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using FileTools
using Base.Test

# array.json
@test JSONFile("array.json") == f"array.json"
@test JSONFile == typeof(f"array.json")
@test 20 == length(f"array.json")
@test 413 == filesize(f"array.json")
@test f"array.json"[1] == json"""[[1,"hello","world"]]"""[1]
@test [18,19,20] == map(first, f"array.json"[end-2:end])
@test 210 == sum(map(first, f"array.json"))

# dict.json
@test JSONFile("dict.json") == f"dict.json"
@test JSONFile == typeof(f"dict.json")
@test 2 == length(f"dict.json")
@test 25 == filesize(f"dict.json")
@test parse(f"dict.json") == json"""{"id": 1, "value": true}"""
@test 1 == f"dict.json"["id"]
@test true == f"dict.json"["value"]
@test 1 == json"""{"id": 1, "value": true}"""["id"]

@test ["id", "value"] == keys(f"dict.json")
@test [1, true] == values(f"dict.json")

# int-array.json
save(f"int-array.json", [1,2,3])
@test [1,2,3] == parse(f"int-array.json")

open(f"int-array.json") do a
  push!(a, 5)
  push!(a, 6)
end
@test [1,2,3,5,6] == parse(f"int-array.json")

push!(f"int-array.json", 7)
@test [1,2,3,5,6,7] == parse(f"int-array.json")

# int-dict.json
save(f"int-dict.json", Dict(2=>3))
@test Dict("2"=>3) == parse(f"int-dict.json")
open(f"int-dict.json") do d
  push!(d, 5=>6)
end
@test Dict("2"=>3, "5"=>6) == parse(f"int-dict.json")

# 0.5
# for a ∈ f"array.json"
for a in f"array.json"
  (idx,hello,world) = a
  #println(idx, hello, world)
  break
end

for (k,v) in f"dict.json"
  #println(k, v)
  break
end


text:
function ll{T<:Real}(Θᵏ::Vector{T}, D, c, J)
    update_Pᵈ_at_xᵈ!(c.Pᵈ, D.xᵈ, Θᵏ, c.Δz̃ᵈ, c.Δẽᵈ, J)

    ll = zero(T)
    lengthxd = length(D.xᵈ)
    @inbounds for j = 1:J
        for i = 1:lengthxd
            ll -=  D.nxjᵈ[i,j]*log(c.Pᵈ[j][i])
        end
    end
    ll/D.nobs
end

∇ll(Θᵏ, grad, D, c, J) =
                  ∇ll(Θᵏ,            grad, D.x,  J, D.nobs, D.xᵈ, c, D)
function ∇ll{T<:Real}(Θᵏ::Vector{T}, grad,   xᵈ, J, n, Dxd, c, D)
    update_Pᵈ_at_xᵈ!(c.Pᵈ, Dxd, Θᵏ, c.Δz̃ᵈ, c.Δẽᵈ, J)

    z̄ = sum(c.Pᵈ[j].*c.z̃ᵈ[j] for j = 1:J)
    grad .= zero(T)
    @views for j in 1:J
        grad[:] -= sum(D.nxjᵈ[:, j].*(c.z̃ᵈ[j]-z̄),1)'
    end
    grad .= grad/n
end

ll_∇ll(Θᵏ, grad, D, c, J) =
                  ll_∇ll(Θᵏ,           grad, D.x,  J, D.nobs, c, D.nxjᵈ, D.xᵈ)
function ll_∇ll{T<:Real}(x::Vector{T}, grad,   xᵈ, J,   n,    c,   nxjᵈ, Dxd)

    K = length(x)
    lengthxd = length(Dxd)
    update_Pᵈ_at_xᵈ!(c.Pᵈ, Dxd, x, c.Δz̃ᵈ, c.Δẽᵈ, J)
    ll = zero(T)
    @inbounds for j = 1:J
        for i = 1:lengthxd
            ll -=  nxjᵈ[i,j]*log(c.Pᵈ[j][i])
        end
    end

    z̄ = sum(c.Pᵈ[j].*c.z̃ᵈ[j] for j = 1:J)
    grad .= zero(T)
    @views for j in 1:J
        grad[:] -= sum(nxjᵈ[:, j].*(c.z̃ᵈ[j]-z̄),1)'
    end
    grad .= grad/n
    ll/n
end

function ∇²ll{T<:Real}(Θᵏ::Vector{T}, hessian, D, c, J)
    update_Pᵈ_at_xᵈ!(c.Pᵈ, D.xᵈ, Θᵏ, c.Δz̃ᵈ, c.Δẽᵈ, J)
    K = size(c.z̃ᵈ[1],2)

    z̄ = sum(c.Pᵈ[j].*c.z̃ᵈ[j] for j = 1:J)
    nxx = sum(D.nxjᵈ,2)
    hessian .= zero(T)
    for j in 1:J
        hessian .+= (nxx.*c.Pᵈ[j].*(c.z̃ᵈ[j]-z̄))'*(c.z̃ᵈ[j]-z̄)
    end

    hessian .= hessian/D.nobs
end

update_Pᵈ_at_xᵈ!{T}(c, D, Θᵏ::Vector{T}, J) =
          update_Pᵈ_at_xᵈ!(c.Pᵈ, D.xᵈ, Θᵏ,          c.Δz̃ᵈ, c.Δẽᵈ, J)
function update_Pᵈ_at_xᵈ!{T}(Pᵈ,   xᵈ, Θᵏ::Vector{T}, Δz̃ᵈ,   Δẽᵈ, J)
    K = size(Δz̃ᵈ[1],2) # nvar
    nxᵈ = length(xᵈ)
    @inbounds for j = 1:J-1
        Pᵈ[j] .= Δẽᵈ[j]
        zj = Δz̃ᵈ[j]
        for k = 1:K
            @simd for i = 1:nxᵈ
                Pᵈ[j][i] += zj[i, k]*Θᵏ[k] # hoist
            end
        end
        Pᵈ[j] .= exp.(Pᵈ[j])
    end
    denom = 1+sum(Pᵈ[1:J-1])
    Pᵈ[end] .= one(T)
    for j = 1:J-1
        Pᵈ[j] .= Pᵈ[j]./denom
        Pᵈ[end] .-= Pᵈ[j]
    end
end


text:
using MolarWeight
using Test

@testset "MolarWeight.jl" begin
    # Write your tests here.
end


text:
export XOR, SWAP, NEG, CONJ
export ROT, IROT

function NEG(a!::Number)
    -a!
end
@selfdual NEG

function CONJ(x::Number)
    x'
end
@selfdual CONJ

function XOR(a!::Number, b::Number)
    a!⊻b, b
end
@selfdual XOR

function SWAP(a!::Number, b!::Number)
    b!, a!
end
@selfdual SWAP

function ROT(i::Number, j::Number, θ::Number)
    a, b = rot(value(i), value(j), value(θ))
    @assign value(i) a
    @assign value(j) b
    i, j, θ
end
function IROT(i::Number, j::Number, θ::Number)
    i, j, _ = ROT(i, j, -θ)
    i, j, θ
end
@dual ROT IROT

for F1 in [:NEG, :CONJ]
    @eval @i function $F1(a!)
        $F1(value(a!))
    end
    @eval NiLangCore.nouts(::typeof($F1)) = 1
    @eval NiLangCore.nargs(::typeof($F1)) = 1
end

for F2 in [:XOR, :SWAP]
    @eval @i function $F2(a, b)
        $F2(value(a), value(b))
    end
    @eval NiLangCore.nouts(::typeof($F2)) = $(F2 == :SWAP ? 2 : 1)
    @eval NiLangCore.nargs(::typeof($F2)) = 2
end

for F3 in [:ROT, :ITOR]
    @eval @i function $F3(a, b, c)
        $F3(value(a), value(b), value(c))
    end
    @eval NiLangCore.nouts(::typeof($F3)) = 2
    @eval NiLangCore.nargs(::typeof($F3)) = 3
end

for (TP, OP) in [(:PlusEq, :+), (:MinusEq, :-), (:XorEq, :⊻)]
    @eval @i function (inf::$TP)(out!, args...; kwargs...)
        inf(value(out!), value.(args)...; kwargs...)
    end
    @eval NiLangCore.nouts(::$TP) = 1

    for SOP in [:*, :/, :^]
        @eval NiLangCore.nargs(::$TP{typeof($SOP)}) = 3
    end
    for SOP in [:sin, :cos, :log, :exp]
        @eval NiLangCore.nargs(::$TP{typeof($SOP)}) = 2
    end
end


text:
module DependencyWalker

using ObjectFile
import Libdl

export Library

struct Library{OH<:Union{ObjectHandle,Missing}}
    path::String
    handle::OH
    level::Int
    deps::Vector{Library}
end

function Library(path::String, level::Int = 0)
    if !isfile(path)
        # TODO: should check also if it can be opened?
        return Library(path, missing, level, Library[])
    end
    io = open(path, "r")
    oh = try
        readmeta(io)
    catch
        close(io)
        missing
    end
    lib = Library(path, oh, level, dependency_tree(oh, level))
    isopen(io) && close(io)
    return lib
end

Library(path::String, nil::Missing, level::Int) =
    Library{Missing}(path, nil, level, [])

Library(path::AbstractString, handle, level) =
    Library(String(path), handle, level)

function dependency_tree(oh::ObjectHandle, level::Int; dlext::String = Libdl.dlext)
    # Get the list of needed libraries
    deps_names = keys(find_libraries(oh))
    # Get list of already dlopen'ed libraries
    open_dls = Libdl.dllist()
    # Initialise list of dependencies
    deps = Library[]
    for dep in deps_names
        split_dep = split(dep, '.')
        # Get rid of the soversion.  TODO: this is only for GNU/Linux and
        # FreeBSD, no idea what we have to do for the other operating systems.
        idx = findlast(isequal(dlext), split_dep)
        if !isnothing(idx)
            dep = join(split_dep[1:idx], '.')
        end
        # Get the first dlopen'ed library matching the needed library, if any
        idx = findfirst(d -> occursin(dep, basename(d)), open_dls)
        if isnothing(idx)
            # Push a missing library to the list
            push!(deps, Library(dep, missing, level + 1))
        else
            # Push the found library to the list
            push!(deps, Library(open_dls[idx], level + 1))
        end
    end
    return deps
end

dependency_tree(::Missing, level) = Library[]

function Base.show(io::IO, lib::Library{<:ObjectHandle})
    print(io, repeat("  ", lib.level), "* ", lib.path)
    for dep in lib.deps
        println(io)
        show(io, dep)
    end
end

Base.show(io::IO, lib::Library{Missing}) =
    printstyled(io, repeat("  ", lib.level), "* ", lib.path, " (NOT FOUND)", color = :yellow)

end # module


text:
module ApproxManifoldProducts

using Reexport
@reexport using KernelDensityEstimate
@reexport using TransformUtils  # likely to be deprecated
using Random

import ManifoldsBase
import ManifoldsBase: AbstractManifold
const MB = ManifoldsBase
using Manifolds

using DocStringExtensions

using NLsolve
import Optim
using CoordinateTransformations
using Requires
# using SLEEFPirates
using LinearAlgebra
using TensorCast
using StaticArrays
using Logging
using Statistics

import Base: *, isapprox, convert
import LinearAlgebra: rotate!
import Statistics: mean
import KernelDensityEstimate: getPoints, getBW

const AMP = ApproxManifoldProducts
const KDE = KernelDensityEstimate
const TUs = TransformUtils
const CTs = CoordinateTransformations

# TODO temporary for initial version of on-manifold products
KDE.setForceEvalDirect!(true)

export
  # new local features
  AMP,
  MKD,
  AbstractManifold,
  ManifoldKernelDensity,
  get2DLambda,
  get2DMu,
  get2DMuMin,
  resid2DLinear,
  solveresid2DLinear!,
  solveresid2DLinear,
  *,
  isapprox,

  # APi and util functions
  buildHybridManifoldCallbacks,
  getKDEManifoldBandwidths,
  manifoldProduct,
  manikde!,
  calcCovarianceBasic,
  isPartial,
  mean,
  calcProductGaussians

# internal features not exported
include("_BiMaps.jl")

# AMP types and some legacy support
include("entities/ManifoldKernelDensity.jl")
include("entities/ManifoldDefinitions.jl")
include("Legacy.jl")
include("services/ManifoldPartials.jl")
include("Interface.jl")

# regular features
include("CommonUtils.jl")
include("services/ManifoldKernelDensity.jl")
include("services/Euclidean.jl")
include("services/CircularUtils.jl")
include("services/Circular.jl")
include("KernelHilbertEmbeddings.jl")

include("TrackingLabels.jl")

# include("Serialization.jl") # moved downstream to IIF to use InferenceVariable serialized types instead
include("API.jl")

include("Deprecated.jl")

function __init__()
  @require Gadfly="c91e804a-d5a3-530f-b6f0-dfbca275c004" begin
    @require Colors="5ae59095-9a9b-59fe-a467-6f913c188581" include("plotting/CircularPlotting.jl")
  end
end

end


text:
export plot

if Pkg.installed("MATLAB") == Nothing()

	function missingPackage()

		warn("Install Julia package MATLAB to enable TensorMesh3D mesh plotting.")
		return

	end

	plot(mesh::TensorMesh3D)                      = missingPackage()
	plot(mesh::TensorMesh3D, u::Array{Int64,1})   = missingPackage()
	plot(mesh::TensorMesh3D, u::Array{Int64,3})   = missingPackage()
	plot(mesh::TensorMesh3D, u::Array{Float64,1}) = missingPackage()
	plot(mesh::TensorMesh3D, u::Array{Float64,3}) = missingPackage()

else

	using MATLAB

	function plot(mesh::TensorMesh3D; f = 1)

		v = diag(getVolume(mesh))
		u = log10(v ./ minimum(v)) .+ 1.0

		plot(mesh, u; f = f)

		return

	end

	function plot(mesh::TensorMesh3D, u::Array{Int64,1}; f = 1)

		u = reshape(u, (mesh.n[1], mesh.n[2], mesh.n[3]))

		plot(mesh, u; f = f)

	end

	function plot(mesh::TensorMesh3D, u::Array{Float64,1}; f = 1)

		u = reshape(u, (mesh.n[1], mesh.n[2], mesh.n[3]))

		plot(mesh, u; f = f)

	end

	function plot(mesh::TensorMesh3D, u::Array{Int64,3}; f = 1)

		if (mesh.n[1],mesh.n[2],mesh.n[3]) != size(u)
			error("Incompatible array sizes")
		end

		x,y,z = getNodalAxes(mesh)

		# set path to Matlab function plotTensorMesh3D.m which resides in the same directory like this Julia function
		(dname,fname) = splitdir(@__FILE__())
		mxcall(:addpath, 0, dname)

		# plot
		mxcall(:plotTensorMesh3D, 0, f, x, y, z, u)

		return

	end

	function plot(mesh::TensorMesh3D, u::Array{Float64,3}; f = 1)

		if (mesh.n[1],mesh.n[2],mesh.n[3]) != size(u)
			error("Incompatible array sizes")
		end

		x,y,z = getNodalAxes(mesh)

		# set path to Matlab function plotTensorMesh3D.m which resides in the same directory like this Julia function
		(dname,fname) = splitdir(@__FILE__())
		mxcall(:addpath, 0, dname)

		# plot
		mxcall(:plotTensorMesh3D, 0, f, x, y, z, u)

		return

	end

end


text:
module GibbsTypePriors

export Pkn_NGG, Pkn_2PD, Pkn_PY, Pkn_Dirichlet, Pkn_Dirichlet_mult, Pkn_PYM, Pkn_NGGM, E_2PD, E_PY, E_stable, E_Dirichlet

include("common_functions.jl")
include("Cnk.jl")
include("Cnk_variable_prec.jl")
include("Vnk.jl")
include("Vnk_variable_prec.jl")
include("Pkn.jl")
include("Pkn_variable_prec.jl")
include("Expect_Kn.jl")
end # module

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
__precompile__()

module OREnvironment

using ORInterface

include("Status.jl")
include("Solution.jl")
include("Constraints.jl")
include("Problem.jl")
include("Movements.jl")

end


text:
export DoubleLearner

using Random

"""
    DoubleLearner(;L1, L2, rng=Random.GLOBAL_RNG)

This is a meta-learner, it will randomly select one learner and update another learner.
The estimation of an observation is the sum of result from two learners.
"""
Base.@kwdef struct DoubleLearner{T1<:Any,T2<:Any,R<:AbstractRNG} <: Any
    L1::T1
    L2::T2
    rng::R = Random.GLOBAL_RNG
end

(learner::DoubleLearner)(env) = learner.L1(env) .+ learner.L2(env)

function RLBase.update!(
    L::DoubleLearner{<:TDLearner},
    t::Any,
    ::AbstractEnv,
    ::PostEpisodeStage,
)
    if rand(L.rng, Bool)
        L, Lₜ = L.L1, L.L2
    else
        L, Lₜ = L.L2, L.L1
    end

    S = t[:state]
    A = t[:action]
    R = t[:reward]
    n, γ, Q, Qₜ = L.n, L.γ, L.approximator, Lₜ.approximator
    G = 0.0
    for i in 1:min(n + 1, length(R))
        G = R[end-i+1] + γ * G
        s, a = S[end-i], A[end-i]
        update!(Q, (s, a) => Q(s, a) - G)
    end
end

function RLBase.update!(L::DoubleLearner{<:TDLearner}, t::Any, ::AbstractEnv, ::PreActStage)
    if rand(L.rng, Bool)
        L, Lₜ = L.L1, L.L2
    else
        L, Lₜ = L.L2, L.L1
    end

    S = t[:state]
    A = t[:action]
    R = t[:reward]
    n, γ, Q, Qₜ = L.n, L.γ, L.approximator, Lₜ.approximator

    if length(R) >= n + 1
        s, a, s′ = S[end-n-1], A[end-n-1], S[end]
        G =
            discount_rewards_reduced(@view(R[end-n:end]), γ) +
            γ^(n + 1) * Qₜ(s′, argmax(Q(s′)))
        update!(Q, (s, a) => Q(s, a) - G)
    end
end

function RLBase.update!(
    t::Any,
    # not very elegant
    ::Union{
        QBasedPolicy{<:DoubleLearner{<:TDLearner}},
        NamedPolicy{<:QBasedPolicy{<:DoubleLearner{<:TDLearner}}},
    },
    ::AbstractEnv,
    ::PreEpisodeStage,
)
    empty!(t)
end


text:
# # GARTEUR SM-AG19 Testbed: Construction of the geometry

# Source code: [`garteur_geometry_tut.jl`](garteur_geometry_tut.jl)

# ## Description

# This virtual test application is based on the test article
# used by the GARTEUR Structures & Materials Action Group 19
# which organized a Round Robin exercise where 12 European laboratories
# tested a single structure between 1995 and 1997. The benchmark structure
# was a laboratory structure built to simulate the dynamic behaviour
# of an aeroplane. The structure was initially built for a benchmark
# study on experimental modal analysis conducted by the
# Structures and Materials Action Group (SM-AG19) of the Group
# for Aeronautical Research and Technology in EURope (GARTEUR).
# The test-bed was designed and manufactured by ONERA, France.

# ![](IMAC97photo.png)

# ### References

# - [GARTEUR] Ground Vibration Test Techniques, compiled by A Gravelle, GARTEUR
#   Structures & Materials Action Group 19 Technical report TP-115, 1999.
# - [BW] Etienne Balmes, Jan R. Wright, GARTEUR GROUP ON GROUND VIBRATION
#   TESTING | RESULTS FROM THE TEST OF A SINGLE STRUCTURE BY 12 LABORATORIES IN
#   EUROPE, Proceedings of DETC'97, 1997 ASME Design Engineering Technical
#   Conferences, September 14-17, 1997, Sacramento, California.

# ## Goals

# - Show how to construct model from multiple connected beams.
# - Demonstrate the use of massless connectors.
# - Visualize the structure interactively.
#

# ## Geometry of the testbed airplane.

# The aluminum testbed was a rather simple structure which was reasonably
# dynamically representative of a simple airplane structure [GARTEUR](@ref
# References). It was composed of several beams simulating a fuselage with wings
# and a tail. Wing tip drums allowed to adjust bending and torsion frequencies
# similarly to airplane ones, with some very close modal frequencies.

# ![](garteur-geom.png)

# The structural-element finite element code relies on the basic functionality
# implemented in this general FEM package.

using FinEtools

# This is the characteristic length. The dimensions of the aircraft frame are
# expressed in terms of multiples of this characteristic unit.
L = 0.1*phun("m");

# ## Cross-section

# Cross-sectional properties are incorporated in the cross-section property.
# There are several rectangular cross-sections in the model: the fuselage, the
# wing, the tail. There are also three massless connectors: connections between
# the fuselage and the wing, between the wing structure and the viscoelastic
# damping layer, and between the fuselage and the tail.

using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionRectangle
# Body of the frame (fuselage).
cs_body = CrossSectionRectangle(s -> 1.5*L, s -> L/2, s -> [1.0, 0.0, 1.0]; label = 1)
# Wing beam.
cs_wing = CrossSectionRectangle(s -> L/10, s -> L, s -> [0.0, 0.0, 1.0]; label = 2)
# Wing drums.
cs_drum = CrossSectionRectangle(s -> L/10, s -> L, s -> [0.0, 0.0, 1.0]; label = 3)
# Vertical part of the tail.
cs_tailv = CrossSectionRectangle(s -> L, s -> L/10, s -> [1.0, 0.0, 1.0]; label = 4)
# Horizontal part of the tail.
cs_tailh = CrossSectionRectangle(s -> L/10, s -> L, s -> [0.0, 0.0, 1.0]; label = 5)
# Constraining plate on top of the viscoelastic tape.
cs_constrp = CrossSectionRectangle(s -> L*(1.1/100), s -> L*(76.2/100), s -> [0.0, 0.0, 1.0]; label = 6)
# Massless connectors of the structural parts of the wing: the main beam and the
# constraining plate.
cs_connw2c = CrossSectionRectangle(s -> L/2, s -> L/2, s -> [1.0, 0.0, 1.0]; label = 7)
# Massless connectors of the wing and the drums.
cs_connw2d = CrossSectionRectangle(s -> L/2, s -> L/2, s -> [1.0, 0.0, 1.0]; label = 8)
# Massless connector between the body and the wings.
cs_connb2w = CrossSectionRectangle(s -> L/5, s -> L, s -> [0.0, 1.0, 0.0]; label = 9)
# Massless connector between the body and the tail.
cs_conntb2t = CrossSectionRectangle(s -> L, s -> L/3, s -> [1.0, 0.0, 1.0]; label = 10)
# Massless connector between the structure and the sensors and point masses.
cs_connta2p = CrossSectionRectangle(s -> L/5, s -> L/5, s -> [1.0, 0.0, 1.0]; label = 11)

# ## Mesh

# We shall use the `frame_member` utility function to generate the mesh of the individual
# parts. This will result in a number of separate meshes for the members. These
# separate meshes will then be glued together (merged) based on the tolerance on
# the location of the nodes.
using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member
tolerance = L/10000;
# Number of intervals from 0.25*L to 8.5*L (the extent of the constraining plate).
nc = 8

# We start with an empty array, and for each part that we mesh, we append the
# mesh to the array. All the meshes in the array will then be considered for
# merging based on geometrical tolerances.
meshes = Tuple{FENodeSet, AbstractFESet}[]

# Define the constituent parts of the body of the aircraft.
push!(meshes, frame_member([-9*L 0 0; -8.5*L 0 0], 1, cs_body; label = cs_body.label))
push!(meshes, frame_member([-8.5*L 0 0; -8.0*L 0 0], 1, cs_body; label = cs_body.label))
push!(meshes, frame_member([-8.0*L 0 0; -2.0*L 0 0], 2, cs_body; label = cs_body.label))
push!(meshes, frame_member([-2.0*L 0 0; 0 0 0], 1, cs_body; label = cs_body.label))
push!(meshes, frame_member([0 0 0; 6*L 0 0], 2, cs_body; label = cs_body.label))

# Define the aluminum parts of the wings.
push!(meshes, frame_member([0 0 0.805*L;  0 0.25*L 0.805*L], 1, cs_wing; label = cs_wing.label))
push!(meshes, frame_member([0 0 0.805*L;  0 -0.25*L 0.805*L], 1, cs_wing; label = cs_wing.label))
push!(meshes, frame_member([0 0.25*L 0.805*L;  0 8.5*L 0.805*L], nc, cs_wing; label = cs_wing.label))
push!(meshes, frame_member([0 -0.25*L 0.805*L;  0 -8.5*L 0.805*L], nc, cs_wing; label = cs_wing.label))
push!(meshes, frame_member([0 8.5*L 0.805*L;  0 9.5*L 0.805*L], 1, cs_wing; label = cs_wing.label))
push!(meshes, frame_member([0 -8.5*L 0.805*L;  0 -9.5*L 0.805*L], 1, cs_wing; label = cs_wing.label))
push!(meshes, frame_member([0 9.5*L 0.805*L;  0 10.0*L 0.805*L], 1, cs_wing; label = cs_wing.label))
push!(meshes, frame_member([0 -9.5*L 0.805*L;  0 -10.0*L 0.805*L], 1, cs_wing; label = cs_wing.label))

# Define the drums at the ends of the wings.
push!(meshes, frame_member([0 +9.5*L +0.91*L; +2*L +9.5*L +0.91*L], 1, cs_drum; label = cs_drum.label))
push!(meshes, frame_member([0 +9.5*L +0.91*L; -2*L +9.5*L +0.91*L], 1, cs_drum; label = cs_drum.label))
push!(meshes, frame_member([0 -9.5*L +0.91*L; +2*L -9.5*L +0.91*L], 1, cs_drum; label = cs_drum.label))
push!(meshes, frame_member([0 -9.5*L +0.91*L; -2*L -9.5*L +0.91*L], 1, cs_drum; label = cs_drum.label))

# Define the horizontal and vertical parts of the tail.
push!(meshes, frame_member([-8*L 0 .75*L; -8*L 0 3.35*L], 2, cs_tailv; label = cs_tailv.label))
push!(meshes, frame_member([-8*L 0 3.35*L; -8*L 0 3.75*L], 2, cs_tailv; label = cs_tailv.label))
push!(meshes, frame_member([-8*L 0 3.8*L; -8*L 2*L 3.8*L], 2, cs_tailh; label = cs_tailh.label))
push!(meshes, frame_member([-8*L 0 3.8*L; -8*L -2*L 3.8*L], 2, cs_tailh; label = cs_tailh.label))

# Define the parts of the aluminum constraining plate for the viscoelastic layer.
push!(meshes, frame_member([-.119*L 0 0.8665*L;  -.119*L 0.25*L 0.8665*L], 1, cs_constrp; label = cs_constrp.label))
push!(meshes, frame_member([-.119*L 0 0.8665*L;  -.119*L -0.25*L 0.8665*L], 1, cs_constrp; label = cs_constrp.label))
push!(meshes, frame_member([-.119*L 0.25*L 0.8665*L;  -.119*L 8.5*L 0.8665*L], nc, cs_constrp; label = cs_constrp.label))
push!(meshes, frame_member([-.119*L -0.25*L 0.8665*L;  -.119*L -8.5*L 0.8665*L], nc, cs_constrp; label = cs_constrp.label))

# Define the massless connectors between:
# Wing - Wingdrum
push!(meshes, frame_member([0 +9.5*L +0.805*L;0 +9.5*L +0.91*L], 1, cs_connw2d; label = cs_connw2d.label))
push!(meshes, frame_member([0 -9.5*L +0.805*L;0 -9.5*L +0.91*L], 1, cs_connw2d; label = cs_connw2d.label))
# Body-Wing
push!(meshes, frame_member([0 0 0; 0 0.25*L .805*L], 1, cs_connb2w; label = cs_connb2w.label))
push!(meshes, frame_member([0 0 0; 0 -0.25*L .805*L], 1, cs_connb2w; label = cs_connb2w.label))
# Body-Tail
push!(meshes, frame_member([-8*L 0 0; -8*L 0 .75*L], 1, cs_conntb2t; label = cs_conntb2t.label))
# Tail-Taildrum
push!(meshes, frame_member([-8*L 0 3.75*L; -8*L 0 3.8*L], 1, cs_connw2d; label = cs_connw2d.label))

# Wing-Constraining plate for the viscoelastic layer
# Middle connector created individually
push!(meshes, frame_member([0 0 .805*L; -.119*L 0 0.8665*L], 1, cs_connw2c; label = cs_connw2c.label))
# Connectors alongside both wings
for i in 1:nc+1
    push!(meshes, frame_member([0 (0.25+(i-1)*8.25/nc)*L .805*L;  -.119*L (0.25+(i-1)*8.25/nc)*L 0.8665*L], 1, cs_connw2c; label = cs_connw2c.label))
    push!(meshes, frame_member([0 -(0.25+(i-1)*8.25/nc)*L .805*L;  -.119*L -(0.25+(i-1)*8.25/nc)*L 0.8665*L], 1, cs_connw2c; label = cs_connw2c.label))
end

# Massless Sensor Connectors
# Tail Sensors
push!(meshes, frame_member([-8*L 2*L 3.8*L; -(153/20)*L (37/20)*L 3.85*L], 1, cs_connta2p; label = cs_connta2p.label))# 303
push!(meshes, frame_member([-8*L -2*L 3.8*L; -(153/20)*L -(37/20)*L 3.85*L], 1, cs_connta2p; label = cs_connta2p.label))# 301

# Wingdrum Sensors
push!(meshes, frame_member([0 9.5*L .91*L ; 0 9.8*L .96*L], 1, cs_connta2p; label = cs_connta2p.label))# 101
push!(meshes, frame_member([-2*L 9.5*L .91*L ; -1.8*L 9.8*L .96*L], 1, cs_connta2p; label = cs_connta2p.label))# 112
push!(meshes, frame_member([2*L 9.5*L .91*L ; 1.8*L 9.8*L .96*L], 1, cs_connta2p; label = cs_connta2p.label))# 111

push!(meshes, frame_member([0 -9.5*L .91*L ; 0 -9.8*L .96*L], 1, cs_connta2p; label = cs_connta2p.label))# 1
push!(meshes, frame_member([-2*L -9.5*L .91*L ; -1.8*L -9.8*L .96*L], 1, cs_connta2p; label = cs_connta2p.label))# 12
push!(meshes, frame_member([2*L -9.5*L .91*L ; 1.8*L -9.8*L .96*L], 1, cs_connta2p; label = cs_connta2p.label))# 11

# Wingdrum complementary masses
push!(meshes, frame_member([2*L 9.5*L .91*L ; 1.8*L 9.2*L .96*L], 1, cs_connta2p; label = cs_connta2p.label))# added mass
push!(meshes, frame_member([2*L -9.5*L .91*L ; 1.8*L -9.2*L .96*L], 1, cs_connta2p; label = cs_connta2p.label))# added mass

# ## Merge all the individual members into one coherent structure

# Merge all the meshes of individual parts. This will glue together nodes which
# are in the "same" location. The parts of the mesh can be distinguished based
# on the label.
fens, fesa = mergenmeshes(meshes, tolerance)

# The result is the set of the finite element nodes, and an array of the finite element sets.
# The sets of the finite elements can be distinguished based on the label.

# The number of nodes in the mesh:
@show count(fens)

# The number of the finite element sets:
@show length(fesa)

# The labels of the finite element sets
@show [s.label[1] for s in fesa]'

# End of the tutorial
nothing


text:
@testset "tarray" begin
    @testset "tzeros and tfill" begin
        a = tzeros(Float64, 10)
        @test eltype(a) === Float64
        @test all(iszero, a)
        @test size(a) == (10,)

        a = tzeros(10)
        @test eltype(a) === Float64
        @test all(iszero, a)
        @test size(a) == (10,)

        a = tfill(5, 10)
        @test eltype(a) === Int
        @test all(==(5), a)
        @test size(a) == (10,)

        a = tfill(19, (5, 5, 5, 5))
        @test eltype(a) === Int
        @test all(==(19), a)
        @test size(a) == (5, 5, 5, 5)
    end

    @testset "push! and pop!" begin
        ta1 = TArray{Int}(4)
        push!(ta1, 1)
        push!(ta1, 2)
        @test pop!(ta1) == 2

        # another constructor
        ta1_2 = TArray{Int, 1}(4)
        push!(ta1_2, 1)
        push!(ta1_2, 2)
        @test pop!(ta1_2) == 2
    end

    @testset "constructors and conversion" begin
        ta2 = TArray{Int}(4, 4)
        @test ta2 isa TArray{Int,2}
        @test size(ta2) == (4, 4)

        ta2 = TArray{Int}(undef, 4, 4)
        @test ta2 isa TArray{Int,2}
        @test size(ta2) == (4, 4)

        ta2 = TArray{Int,2}(4, 4)
        @test ta2 isa TArray{Int,2}
        @test size(ta2) == (4, 4)

        ta2 = TArray{Int,2}(undef, 4, 4)
        @test ta2 isa TArray{Int,2}
        @test size(ta2) == (4, 4)

        @test_throws MethodError TArray{Int,2}(4)
        @test_throws MethodError TArray{Int,2}(undef, 4)

        ta3 = TArray{Int, 4}(4, 3, 2, 1)
        ta4 = Libtask._get_local_storage(ta3)
        @test ta3[3] == ta4[3]

        ta5 = TArray{Int}(4)
        @test ta5 isa TArray{Int,1}
        @test size(ta5) == (4,)

        ta5 = TArray{Int}(undef, 4)
        @test ta5 isa TArray{Int,1}
        @test size(ta5) == (4,)

        ta5 = TArray{Int,1}(4)
        @test ta5 isa TArray{Int,1}
        @test size(ta5) == (4,)

        ta5 = TArray{Int,1}(undef, 4)
        @test ta5 isa TArray{Int,1}
        @test size(ta5) == (4,)

        @test_throws MethodError TArray{Int,1}(4, 4)
        @test_throws MethodError TArray{Int,1}(undef, 4, 4)

        for i in 1:4
            ta5[i] = i
        end
        @test Array(ta5) == [1, 2, 3, 4]
        @test convert(Array, ta5) == [1, 2, 3, 4]
        @test convert(Array{Int, 1}, ta5) == [1, 2, 3, 4]
        @test ta5 == convert(TArray, [1, 2, 3, 4])
        @test ta5 == convert(TArray{Int, 1}, [1, 2, 3, 4])
        @test_throws MethodError convert(TArray{Int, 2}, [1, 2, 3, 4])
        @test_throws MethodError convert(Array{Int, 2}, ta5)

        @test Array(tzeros(4)) == zeros(4)

        ta6 = TArray{Float64}(4)
        for i in 1:4
            ta6[i] = i / 10
        end
        @test ta6[1] == 0.1
        @test Array(ta6) == [0.1, 0.2, 0.3, 0.4]

        # TODO: add test for use this multi-dim array
        ta7 = TArray{Int, 2}((2, 2))
    end

    @testset "stdlib functions" begin
        ta = TArray{Int}(4, 4)

        @test view(ta, 3:5) isa TArray{Int, 1}
        @test view(ta, 3:5) == ta[3:5]

        @test -ta isa TArray{Int, 2}
        @test -(-ta) == ta

        @test transpose(ta)[2, 3] == ta[3, 2]

        @test repeat(ta, 2) == vcat(ta, ta)

        @test repeat(ta, 1, 2) == hcat(ta, ta)

        @test ta .+ ta == Libtask._get_local_storage(ta) .+ Libtask._get_local_storage(ta)
    end

    @testset "task copy" begin
        function f()
            t = TArray(Int, 1)
            t[1] = 0
            while true
                produce(t[1])
                t[1]
                t[1] = 1 + t[1]
            end
        end

        ctask = CTask(f)

        consume(ctask)
        consume(ctask)
        a = copy(ctask)
        consume(a)
        consume(a)

        @test consume(ctask) == 2
        @task consume(a) == 4

        DATA = Dict{Task, Array}()
        function g()
            ta = tzeros(UInt64, 4)
            for i in 1:4
                ta[i] = hash(Libtask._current_task())
                DATA[Libtask._current_task()] = convert(Array, ta)
                produce(ta[i])
            end
        end

        ctask = CTask(g)
        @test consume(ctask) == hash(ctask.task) # index = 1
        @test consume(ctask) == hash(ctask.task) # index = 2

        a = copy(ctask)
        @test consume(a) == hash(a.task) # index = 3
        @test consume(a) == hash(a.task) # index = 4

        @test consume(ctask) == hash(ctask.task) # index = 3

        @test DATA[ctask.task] == [hash(ctask.task), hash(ctask.task), hash(ctask.task), 0]
        @test DATA[a.task] == [hash(ctask.task), hash(ctask.task), hash(a.task),
                               hash(a.task)]
    end

    @testset "Issue: PR-86 (DynamicPPL.jl/pull/261)" begin
        function f()
            t = TArray(Int, 1)
            t[1] = 0
            while true
                produce(t[1])
                t[1]
                t[1] = 1 + t[1]
            end
        end

        ctask = CTask(f)

        ex = try
            for _ in 1:10000
                consume(ctask)
                consume(ctask)
                a = copy(ctask)
                consume(a)
                consume(a)
            end
        catch ex
            ex
        end
        @test ex === nothing
    end

end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

see "A Direct Formulation for Sparse PCA Using Semidefinite Programming" by
Alexandre d’Aspremont, Laurent El Ghaoui, Michael I. Jordan, Gert R. G. Lanckriet
=#

import Distributions
using SparseArrays

struct SparsePCANative{T <: Real} <: ExampleInstanceNative{T}
    p::Int
    k::Int
    use_epinorminfdual::Bool # use dual of epinorminf cone, else nonnegative cones
    noise_ratio::T
end

function build(inst::SparsePCANative{T}) where {T <: Real}
    (p, k, noise_ratio) = (inst.p, inst.k, inst.noise_ratio)
    @assert 0 < k <= p

    # sample components that will carry the signal
    signal_idxs = Distributions.sample(1:p, k, replace = false)
    if noise_ratio <= 0
        # noiseless model
        x = zeros(T, p)
        x[signal_idxs] = rand(T, k)
        sigma = x * x'
        sigma ./= tr(sigma)
    else
        # simulate some observations with noise
        x = randn(T, p, 100)
        sigma = x * x'
        y = rand(Distributions.Normal(zero(T), noise_ratio), k)
        sigma[signal_idxs, signal_idxs] .+= y * y'
        sigma ./= 100
    end

    dimx = Cones.svec_length(p)
    # x is the svec (lower triangle, row-wise) of the matrix solution
    c = Cones.smat_to_svec!(zeros(T, dimx), -sigma, sqrt(T(2)))
    b = T[1]
    A = zeros(T, 1, dimx)
    for i in 1:p
        A[sum(1:i)] = 1
    end
    hpsd = zeros(T, dimx)
    cones = Cones.Cone{T}[Cones.PosSemidefTri{T, T}(dimx)]

    if inst.use_epinorminfdual
        # l1 cone
        # double off-diagonals, which are already scaled by rt2
        Gl1vec = fill(-one(T), dimx)
        Cones.scale_svec!(Gl1vec, sqrt(T(2)))
        G = [
            sparse(-one(T) * I, dimx, dimx);
            spzeros(T, 1, dimx);
            Diagonal(Gl1vec);
            ]
        h = vcat(hpsd, T(k), zeros(T, dimx))
        push!(cones, Cones.EpiNormInf{T, T}(1 + dimx, use_dual = true))
    else
        l1 = Cones.scale_svec!(ones(T, dimx), sqrt(T(2)))
        G = [
            -I    spzeros(T, dimx, 2 * dimx);
            spzeros(T, 2 * dimx, dimx)    -I;
            spzeros(T, 1, dimx)    repeat(l1', 1, 2);
            ]
        A = [
            sparse(A)    spzeros(T, 1, 2 * dimx);
            -I    -I    I;
            ]
        c = vcat(c, zeros(T, 2 * dimx))
        b = vcat(b, zeros(T, dimx))
        h = vcat(hpsd, zeros(T, 2 * dimx), k)
        push!(cones, Cones.Nonnegative{T}(2 * dimx + 1))
    end

    model = Models.Model{T}(c, A, b, G, h, cones)
    return model
end

function test_extra(
    inst::SparsePCANative{T},
    solve_stats::NamedTuple,
    ::NamedTuple,
    ) where T
    @test solve_stats.status == Solvers.Optimal
    if solve_stats.status == Solvers.Optimal && iszero(inst.noise_ratio)
        # check objective value is correct
        tol = eps(T)^0.25
        @test solve_stats.primal_obj ≈ -1 atol=tol rtol=tol
    end
end


text:
export OrthoCube
export find_generators
export generate_coordinates
export volume
export isequiv

"""
    OrthoCube(shape)

Represent a hypercubic (Bravais) lattice.

# Fields
* `shape_matrix`: a matrix whose columns are the lattice vectors of the
* `inverse_shape_matrix`: `(shape_matrix)`
* `wrap`: periodic wrapping function which takes an integer array, and maps it onto a
  site in the Bravais lattice. Analogous to `fldmod`.
"""
struct OrthoCube
    shape_matrix::Matrix{Int}
    inverse_shape_matrix::Matrix{Rational{Int}}
    wrap::Function

    function OrthoCube(shape_matrix::AbstractMatrix{<:Integer})
        dim, dim2 = size(shape_matrix)
        if dim != dim2
            msg = "scale_matrix is not square: dimensions are ($dim, $dim2)"
            throw(DimensionMismatch(msg))
        end
        det = ExactLinearAlgebra.determinant(shape_matrix)
        det == 0 && throw(ArgumentError("scale matrix null"))

        inverse_shape_matrix = ExactLinearAlgebra.inverse(shape_matrix)
        function wrap(r::AbstractArray{<:Integer}, mode::RoundingMode=RoundDown)
            rnd = (x) -> round(Int, x, mode)
            R = rnd.(inverse_shape_matrix * r)
            r2 = r - shape_matrix * R
            return R, r2
        end
        return new(shape_matrix, inverse_shape_matrix, wrap)
    end
end

"""
    dimension(ortho::OrthoCube)

Return the spatial dimension of the orthocube.
"""
dimension(ortho::OrthoCube) = size(ortho.shape_matrix, 1)

"""
    volume(ortho::OrthoCube)

Return the signed volume of the orthocube, defined by the determinant of the shape.
"""
volume(ortho::OrthoCube) = ExactLinearAlgebra.determinant(ortho.shape_matrix)

function Base.:(==)(lhs::OrthoCube, rhs::OrthoCube)
    return (lhs.shape_matrix == rhs.shape_matrix)
end

"""
    isequiv(lhs::OrthoCube, rhs::OrthoCube)

Check whether the two orthocubes are equivalent.
"""
function isequiv(lhs::OrthoCube, rhs::OrthoCube)
    R, r = lhs.wrap(rhs.shape_matrix)
    if abs(ExactLinearAlgebra.determinant(R)) != 1 || !iszero(r)
        return false
    end
    return true
end

# function isequiv_old(lhs::OrthoCube, rhs::OrthoCube)
#     det_lhs = TightBindingLattice.ExactLinearAlgebra.determinant(lhs.shape_matrix)
#     det_rhs = TightBindingLattice.ExactLinearAlgebra.determinant(rhs.shape_matrix)
#     det_lhs != det_rhs && return false
#     inv_lhs = TightBindingLattice.ExactLinearAlgebra.inverse(lhs.shape_matrix)
#     inv_rhs = TightBindingLattice.ExactLinearAlgebra.inverse(rhs.shape_matrix)
#     return all(isinteger.(inv_lhs * rhs.shape_matrix)) && all(isinteger.(inv_rhs * lhs.shape_matrix))
# end

"""
    find_generators(ortho::OrthoCube)

Find generators of an OrthoCube. For one-dimensional lattice, which is isomorphic to Zₙ, the
generator is +1. For two-dimensional lattices, this function invokes
[`find_generators_2d`](@ref). Higher dimensions are currently not supported.
"""
function find_generators(ortho::OrthoCube)
    if dimension(ortho) == 1
        return ones(Int, 1, 1)
    elseif dimension(ortho) == 2
        return find_generators_2d(ortho)
    else
        error("Currently, only 1D and 2D lattices are supported")
    end
end

"""
    find_generators_2d(ortho::OrthoCube)

Find translation generators of an OrthoCube.
"""
function find_generators_2d(ortho::OrthoCube)
    dimension(ortho) != 2 && throw(ArgumentError("shape matrix should be 2x2"))

    orthovolume = abs(volume(ortho))
    allowed_pairs = [[1,0]]
    append!(
        allowed_pairs,
        [x, y] for y in 1:orthovolume, x in 0:orthovolume if gcd(x,y) == 1
    )

    function make_loop(t::Vector{Int})
        axis = Vector{Int}[zero(t)]
        r = ortho.wrap(t)[2]
        while !iszero(r)
            push!(axis, r)
            r = ortho.wrap(r .+ t)[2]
        end
        @assert allunique(axis)
        return axis
    end

    for r1p in allowed_pairs
        r1p[1] == 0 && continue
        r1 = [r1p[1], -r1p[2]]
        axis1 = make_loop(r1)
        mod(orthovolume, length(axis1)) != 0 && continue
        for r2 in allowed_pairs
            ExactLinearAlgebra.determinant(hcat(r1, r2)) != 1 && continue
            axis2 = make_loop(r2)
            orthovolume != length(axis1) * length(axis2) && continue
            coordinates = vec([ortho.wrap(sum(x))[2]
                                   for x in Iterators.product(axis1, axis2)])
            allunique(coordinates) && return hcat(r1, r2)
        end
    end
    error("Failed to find decomposition") # COV_EXCL_LINE
end

"""
    generate_coordinates(ortho, generators)

Generate a list of coordinates of the orthocube

# Arguments
* `ortho::OrthoCube`
* `generator_translations::AbstractMatrix{<:Integer}`
"""
function generate_coordinates(
    ortho::OrthoCube,
    generator_translations::AbstractMatrix{<:Integer},
)
    let dim = dimension(ortho)
        if size(generator_translations) != (dim, dim)
            throw(DimensionMismatch("OrthoCube and generator_translations have different dimensions"))
        elseif abs(ExactLinearAlgebra.determinant(generator_translations)) != 1
            throw(ArgumentError("generator_translations is not unimodular"))
        end
    end
    function make_loop(t::AbstractVector{<:Integer})
        axis = Vector{Int}[zero(t)]
        r = ortho.wrap(t)[2]
        while !iszero(r)
            push!(axis, r)
            r = ortho.wrap(r .+ t)[2]
        end
        return axis
    end
    axes = [make_loop(t) for t in eachcol(generator_translations)]
    coordinates = vec([ortho.wrap(sum(x))[2] for x in Iterators.product(axes...)])

    if length(coordinates) != abs(volume(ortho))
        msg = "translations $generator_translations generates $(length(coordinates)) " *
            "coordinates, while the volume is $(abs(volume(ortho)))"
        throw(ArgumentError(msg))
    end
    @assert allunique(coordinates)
    return coordinates
end


text:
import SwitchOnSafety: isinfeasible, isfeasible, isdecided, usestep
export sdpbound_γ, sdpbound_gamma

function sdpbisection_alg(s::StateDepDiscreteSwitchedLinearSystem, optimizer, lo, hi; verbose=1, tol=1e-5, step=0.5)
    @assert hi > lo
    while hi - lo ≥ tol
        γ_mid = midpoint(hi, lo, step)
        status = sdplyap_feasible_prog(γ_mid, s, optimizer, verbose=verbose-1)
        if verbose ≥ 1
            println("$γ_mid")
        end
        if !isdecided(status, true)
            if usestep(lo, hi)
                step *= 2
                continue
            end
            γ_midlo = min(midpoint(lo, γ_mid, step), shift(γ_mid, -tol/2))
            statuslo = sdplyap_feasible_prog(γ_midlo, s, optimizer, verbose=verbose-1)
            if isdecided(statuslo, true)
                γ_mid = γ_midlo
                status = statuslo
            else
                γ_midhi = max(midpoint(γ_mid, hi, step), shift(γ_mid, tol/2))
                statushi = sdplyap_feasible_prog(γ_midhi, s, optimizer, verbose=verbose-1)
                if isdecided(statushi, true)
                    γ_mid = γ_midhi
                    status = statushi
                end
            end
        end
        if isinfeasible(status, true)
            lo = γ_mid
        end
        if isfeasible(status, true)
            hi = γ_mid
        end
    end
    return lo, hi
end

"""
    sdpbound_γ(s::StateDepDiscreteSwitchedLinearSystem)

Compute minimum value of ``\\gamma`` for which the SDP is feasible.

For keywords, see [`sosbound_γ`](@ref)
"""
function sdpbound_γ(s::StateDepDiscreteSwitchedLinearSystem; optimizer=optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true), tol=1e-5, verbose=0, initstep=1.1)
    lo_init, hi_init = 0, initstep; hi_stat = (MOI.INFEASIBLE, MOI.INFEASIBLE_POINT, MOI.INFEASIBILITY_CERTIFICATE, "")
    while isinfeasible(hi_stat, true)
        hi_stat = sdplyap_feasible_prog(hi_init, s, optimizer, verbose=verbose-1)
        if isinfeasible(hi_stat, true)
            lo_init = copy(hi_init)
            hi_init += initstep
        end
        initstep *= 2
    end
    lo, hi = sdpbisection_alg(s, optimizer, lo_init, hi_init, verbose=verbose)
    return hi
end

"""
    sdpbound_gamma(s::StateDepDiscreteSwitchedLinearSystem; optimizer=nothing, tol=1e-5, verbose=0, initstep=1.1)

Alias for [`sdpbound_γ`](@ref).
"""
sdpbound_gamma(s::StateDepDiscreteSwitchedLinearSystem; optimizer=optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true), tol=1e-5, verbose=0, initstep=1.1) = sdpbound_γ(s::StateDepDiscreteSwitchedLinearSystem; optimizer=optimizer, tol=tol, verbose=verbose, initstep=initstep)


text:
#############################################################################
#Copyright (c) 2016 Charles Le Losq
#
#The MIT License (MIT)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the #Software without restriction, including without limitation the rights to use, copy, #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, #and to permit persons to whom the Software is furnished to do so, subject to the #following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
# functions.jl contains several mathematic functions
#
#
#############################################################################

"""
Bootstrap data generation function

    bootstrap(x::Array{Float64}, y::Array{Float64};boottype = "np")

Bootstrap of data. We generate new datapoints with the basis of existing data. Two options: non-parametric and parametric bootstrapping.

type = "np": Non-parametric bootstrap is made by re-sampling with replacement the data.

type = "p": Parametric bootstrapping is made by re-sampling the data from randowmly picking new values in their probability density distribution.
For now, only normal PDF are supported.

RETURN: b_x_f, b_y_f, the arrays with bootstrapped x and y values.

TODO: proposing other error PDF for parametric bootstrapping.

"""
function bootsample(x::Array{Float64}, y::Array{Float64}; boottype::ASCIIString = "np", ese::Array{Float64} = [0.0])
    testx = size(x,1)
    testy = size(y,1)
    teste = size(ese,1)
    if boottype == "np" && testx == testy
        vect = collect(1:size(x)[1]) # for real bootstrapping
        idx = sample(vect,size(vect)[1],replace=true) #resampling data with replacement...
        b_x_f = x[idx,:] # we pick the right x
        b_y_f = y[idx,:] # we pick the right y
    elseif boottype == "p" && testx == testy && testx == teste
	b_x_f = x
	b_y_f = ones(size(y))
	b_y_f = randn!(b_y_f[:,:]) .* ese[:,:] + y[:,:]
    else
	error("Something is wrong. Check size of x, y and ese as well as the boottype (either p or np). Providing an ese array is mandatory for parametric bootstrapping")
    end

    return b_x_f, b_y_f
end

"""
bootperf: reading the bootstrap parameter array for peak fitting

    bootperf(params_boot::Array{Float64}; plotting::ASCIIString = "True", parameter::Int64 = 0, feature::Int64 = 0, histogram_step::Int64 = 100, savefigures::ASCIIString = "False", save_bootrecord::ASCIIString = "Boot_record.pdf", save_histogram::ASCIIString = "Boot_histogram.pdf")

params_boot[i,j] or [i,j,k]: array with i the bootstrrap experiment, j the parameter and k the feature being calculated (e.g., a peak during peak fitting)

plotting: switch to plotting mode ("True") or not ("False"). If true, parameter and feature must be provided, otherwise an error message is returned.

histogram_step is an integer value to control the histogram X axis division.

savefigures: "True" or "False", explicit, save in the current working directory.

save_bootrecord: Name for the graphic showing the bootstrap mean and std evolutions, with extension.

save_histogram: Name for the graphic showing the histogram for the parameter and feature of interest, with extension.

RETURN: std_record, mean_record, the arrays recording how the standard deviation and mean of the parameters as a function of the bootstrap advance.

"""
function bootperf(params_boot::Array{Float64}; plotting::ASCIIString = "True", parameter::Int64 = 1, feature::Int64 = 1, histogram_step::Int64 = 100, savefigures::ASCIIString = "False", save_bootrecord::ASCIIString = "Boot_record.pdf", save_histogram::ASCIIString = "Boot_histogram.pdf")

    # test
    if plotting != "True" && plotting != "False"
        error("Error: plotting should be set to True or False")
    end

    nb_boot::Int64 = size(params_boot)[1]
    nb_param::Int64 = size(params_boot)[2]

   if ndims(params_boot) == 2 #### TWO DIMENSIONS PARAMS_BOOT

       mean_record::Array{Float64} = zeros(size(params_boot))
       std_record::Array{Float64} = zeros(size(params_boot))
       for k = 1:nb_boot
           for j = 1:size(params_boot)[2]
       	        if k == 1
       	            mean_record[k,j] = params_boot[k,j]
       		    std_record[k,j] = 0
       	        else
       	            mean_record[k,j] = sum(mean(params_boot[1:k,j],1))
       		    std_record[k,j] = sum(std(params_boot[1:k,j],1))
       	        end
       	    end
       	end

       if plotting == "True" # Graphics
           # tests:
           if savefigures != "True" && savefigures != "False"
               error("savefigure should be set to True or False.")
           end
           if isa(save_bootrecord ,ASCIIString) != true
               error("Error: save_bootrecord is not a valid savename. Check it is a string of characters.")
           end
           if isa(save_histogram ,ASCIIString) != true
               error("Error: save_histogram is not a valid savename. Check it is a string of characters.")
           end

           if parameter == 0 || histogram_step == 0
               error("parameter or histogram_step appear to have been set to 0...")
           elseif parameter > size(params_boot)[2]
               error("Value entered for parameter is too high.")
           end

           fig = figure()
           #title("Mean and Std of parameter $(parameter), feature $(feature)",fontsize = 10, fontweight = "bold", fontname="Arial")
           subplot(121)
           plot(1:nb_boot, mean_record[:,parameter],label="Mean")
           xlabel("Number of iterations during bootstrap",fontsize = 10, fontweight = "bold", fontname="Arial")
           ylabel("Mean value",fontsize = 10, fontweight = "bold", fontname="Arial")

           subplot(122)
           plot(1:nb_boot, std_record[:,parameter],label="Std")
           xlabel("Number of iterations during bootstrap",fontsize = 10, fontweight = "bold", fontname="Arial")
           ylabel("Standard deviation",fontsize = 10, fontweight = "bold", fontname="Arial")

           if savefigures == "True"
               savefig(save_bootrecord)
           end

           fig2= figure()
           h = plt[:hist](params_boot[:,parameter],100)
           #title("Histogram of parameter $(parameter), feature $(feature)",fontsize = 10, fontweight = "bold", fontname="Arial")
           xlabel("Parameter value",fontsize = 10, fontweight = "bold", fontname="Arial")
           ylabel("Number of bootstrap values",fontsize = 10, fontweight = "bold", fontname="Arial")

           if savefigures == "True"
               savefig(save_histogram)
       	   end

       end # end if plotting

       return mean_record, std_record

   elseif ndims(params_boot) == 3 #### THREE DIMENSIONS PARAMS_BOOT
       mean_record = zeros(size(params_boot))
       std_record = zeros(size(params_boot))
       for k = 1:nb_boot
           for j = 1:size(params_boot)[2]
   	    for i = 1:size(params_boot)[3]
   	        if k == 1
   	            mean_record[k,j,i] = params_boot[k,j,i]
   		    std_record[k,j,i] = 0
   	        else
   	            mean_record[k,j,i] = sum(mean(params_boot[1:k,j,i],1))
   		    std_record[k,j,i] = sum(std(params_boot[1:k,j,i],1))
   	        end
   	    end
   	end
       end

       if plotting == "True" # Graphics
           # tests:
           if savefigures != "True" && savefigures != "False"
               error("savefigure should be set to True or False.")
           end
           if isa(save_bootrecord ,ASCIIString) != true
               error("Error: save_bootrecord is not a valid savename. Check it is a string of characters.")
           end
           if isa(save_bootrecord ,ASCIIString) != true
               error("Error: save_histogram is not a valid savename. Check it is a string of characters.")
           end

           if parameter == 0 || feature == 0 || histogram_step == 0
               error("parameter, feature or histogram_step appear to have been set to 0...")
           elseif parameter > size(params_boot)[2]
               error("Value entered for parameter is too high.")
           elseif feature > size(params_boot)[3]
   	    error("Value entered for feature is too high.")
           end

           fig = figure()
           #title("Mean and Std of parameter $(parameter), feature $(feature)",fontsize = 10, fontweight = "bold", fontname="Arial")
           subplot(121)
           plot(1:nb_boot, mean_record[:,parameter,feature],label="Mean")
           xlabel("Number of iterations during bootstrap",fontsize = 10, fontweight = "bold", fontname="Arial")
           ylabel("Mean value",fontsize = 10, fontweight = "bold", fontname="Arial")

           subplot(122)
           plot(1:nb_boot, std_record[:,parameter,feature],label="Std")
           xlabel("Number of iterations during bootstrap",fontsize = 10, fontweight = "bold", fontname="Arial")
           ylabel("Standard deviation",fontsize = 10, fontweight = "bold", fontname="Arial")

           if savefigures == "True"
               savefig(save_bootrecord)
           end

           fig2= figure()
           h = plt[:hist](params_boot[:,parameter,feature],100)
           #title("Histogram of parameter $(parameter), feature $(feature)",fontsize = 10, fontweight = "bold", fontname="Arial")
           xlabel("Parameter value",fontsize = 10, fontweight = "bold", fontname="Arial")
           ylabel("Number of bootstrap values",fontsize = 10, fontweight = "bold", fontname="Arial")

           if savefigures == "True"
               savefig(save_histogram)
   	end

       end # end if plotting

       return mean_record, std_record

    elseif ndims(params_boot) > 3
	error("Not implemented, params_boot should have 3 dimensions, maximum.")
    end

end # end function