text:
### A Hamiltonian is composed of blocks (kpoints), which have a list of RealFourierOperator
# corresponding to each term
# This is the "high-level" interface, provided for convenience

struct HamiltonianBlock
    basis::PlaneWaveBasis
    kpoint::Kpoint

    # The operators are vectors of RealFourierOperator,
    # not typed because of type invariance issues.
    operators::Vector  # the original list of RealFourierOperator
                       # (as many as there are terms), kept for easier exploration
    optimized_operators::Vector  # the optimized list of RealFourierOperator, to be used for applying
    scratch  # Pre-allocated scratch arrays for fast application
end
function HamiltonianBlock(basis, kpt, operators, scratch)
    HamiltonianBlock(basis, kpt, operators, optimize_operators_(operators), scratch)
end
Base.eltype(block::HamiltonianBlock) = complex(eltype(block.basis))
Base.size(block::HamiltonianBlock) =
    (length(G_vectors(block.kpoint)), length(G_vectors(block.kpoint)))
Base.size(block::HamiltonianBlock, i::Integer) = i < 3 ? size(block)[i] : 1

struct Hamiltonian
    basis::PlaneWaveBasis
    blocks::Vector{HamiltonianBlock}
end

# Loop through bands, IFFT to get ψ in real space, loop through terms, FFT and accumulate into Hψ
# Do this band by band to conserve memory
# As an optimization we special-case nonlocal operators to apply them
# instead on the full block and benefit from BLAS3
@views @timing "Hamiltonian multiplication" function LinearAlgebra.mul!(Hψ::AbstractArray,
                                                                        H::HamiltonianBlock,
                                                                        ψ::AbstractArray)
    basis = H.basis
    kpt = H.kpoint
    nband = size(ψ, 2)

    # Allocate another scratch array:
    Hψ_fouriers = [similar(Hψ[:, 1]) for tid = 1:Threads.nthreads()]

    # take ψi, IFFT it to ψ_real, apply each term to Hψ_temp and Hψ_real, and add it to Hψi
    Threads.@threads for iband = 1:nband
        tid = Threads.threadid()
        ψ_real = H.scratch.ψ_reals[tid]
        Hψ_real = H.scratch.Hψ_reals[tid]
        Hψ_fourier = Hψ_fouriers[tid]

        Hψ_real .= 0
        Hψ_fourier .= 0
        G_to_r!(ψ_real, basis, kpt, ψ[:, iband])
        for op in H.optimized_operators
            if !(op isa NonlocalOperator)
                apply!((fourier=Hψ_fourier, real=Hψ_real),
                       op,
                       (fourier=ψ[:, iband], real=ψ_real))
            end
        end
        Hψ[:, iband] .= Hψ_fourier
        r_to_G!(Hψ_fourier, basis, kpt, Hψ_real)
        Hψ[:, iband] .+= Hψ_fourier
    end

    # Apply the nonlocal operators
    for op in H.optimized_operators
        if op isa NonlocalOperator
            apply!((fourier=Hψ, real=nothing), op, (fourier=ψ, real=nothing))
        end
    end

    Hψ
end
Base.:*(H::HamiltonianBlock, ψ) = mul!(similar(ψ), H, ψ)

function LinearAlgebra.mul!(Hψ, H::Hamiltonian, ψ)
    for ik = 1:length(H.basis.kpoints)
        mul!(Hψ[ik], H.blocks[ik], ψ[ik])
    end
    Hψ
end
# need `deepcopy` here to copy the elements of the array of arrays ψ (not just pointers)
Base.:*(H::Hamiltonian, ψ) = mul!(deepcopy(ψ), H, ψ)

# Get energies and Hamiltonian
# kwargs is additional info that might be useful for the energy terms to precompute
# (eg the density ρ)
@timing function energy_hamiltonian(basis::PlaneWaveBasis, ψ, occ; kwargs...)
    # it: index into terms, ik: index into kpoints
    ene_ops_arr = [ene_ops(term, ψ, occ; kwargs...) for term in basis.terms]
    energies    = [eh.E for eh in ene_ops_arr]
    operators   = [eh.ops for eh in ene_ops_arr]         # operators[it][ik]
    hks_per_k   = [[blocks[ik] for blocks in operators]  # hks_per_k[ik][it]
                   for ik = 1:length(basis.kpoints)]

    # Preallocated scratch arrays
    T = eltype(basis)
    scratch = (
        ψ_reals=[zeros(complex(T), basis.fft_size...) for tid = 1:Threads.nthreads()],
        Hψ_reals=[zeros(complex(T), basis.fft_size...) for tid = 1:Threads.nthreads()]
    )

    H = Hamiltonian(basis, [HamiltonianBlock(basis, kpt, hks, scratch)
                            for (hks, kpt) in zip(hks_per_k, basis.kpoints)])
    E = Energies(basis.model.term_types, energies)
    (E=E, H=H)
end
function Hamiltonian(basis::PlaneWaveBasis; ψ=nothing, occ=nothing, kwargs...)
    _, H = energy_hamiltonian(basis, ψ, occ; kwargs...)
    H
end

import Base: Matrix, Array
function Matrix(block::HamiltonianBlock)
    sum(Matrix(h) for h in block.optimized_operators)
end
Array(block::HamiltonianBlock) = Matrix(block)

"""
Get the total local potential of the given Hamiltonian, in real space.
"""
function total_local_potential(ham::Hamiltonian)
    @assert ham.basis.model.spin_polarization in (:none, :spinless)
    block = ham.blocks[1]  # all local potentials are the same
    rs = [o for o in block.optimized_operators if o isa RealSpaceMultiplication]
    @assert length(rs) == 1
    rs[1].potential
end


text:
function showtypetree(T, level = 0)
    println("\t"^level, last(split(string(T), ".")))
    for t in subtypes(T)
        showtypetree(t, level + 1)
    end
end

"""
Display the QuNet type tree
"""
function qunet_type_tree()
    showtypetree(QuNet.QObject)
end

ext: jl

max_stars_repo_name: drpeterrohde/QuNet

max_issues_repo_name: SamuelMarks/QuNet

max_forks_repo_name: SamuelMarks/QuNet


text:
using ArgParse
using JLD
using MAT
using JSON
using HDF5
using MFCC

include("util.jl")

function main(args)
    s = ArgParseSettings()
    s.description = "MFSC feature extraction from MAT file"

    @add_arg_table s begin
        ("--matfile"; required=true)
        ("--jsonfile"; required=true)
        ("--savefile"; required=true)
        ("--sr"; default=16000.0; arg_type=Float64)
        ("--preprocess"; action=:store_true)
        ("--period"; default=100; arg_type=Int64)
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); println(o); flush(STDOUT)

    # load mat file
    data = matread(o[:matfile])
    entries = JSON.parsefile(o[:jsonfile])
    # savefile = h5open(o[:savefile], "w")

    println("started.", now())
    features = Dict()
    @time for (i,entry) in enumerate(entries)
        longid = entry["longid"]
        sample = data[longid]
        feats = mfsc(sample)
        features[longid] = feats[1]
        # write(savefile, longid, feats[1])

        if i % o[:period] == 0
            println("$i samples processed so far...")
        end
    end
    # close(savefile)
    save(o[:savefile], "features", features)
    println("done.")
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ilkerkesen/phoneme-convnet

max_issues_repo_name: ilkerkesen/phoneme-convnet

max_forks_repo_name: ilkerkesen/phoneme-convnet


text:
function evaluate_matrix_for_w!(eval, template, w)
  for i in 1:length(template)
    if typeof(template[i]) <: Function
      eval[i] = template[i](w)
    end
  end
end

function get_prms_from_pairs(pairs)
  p = parameters()
  for pair in pairs
    setfield!(p, Symbol(pair[1]), float(pair[2]))
  end
  return p
end

function material_matrix_to_impedance(
            material_matrix::Array = [1 1 0 1; 0 1 0 1; 0 1 0 1; 0 1 0 1],
            prms_pairs = [];
            #
            f_list,
            #
            complex_type=ComplexF64,
            iterative_solver = false,
            verbose = false
            )

  if verbose
    @time params = get_prms_from_pairs(prms_pairs)
    @time header, sp_input, b_real = material_matrix_to_lin_sys(material_matrix, params)

    nz_el_count = length(sp_input[3])
    sp_input_vals_eval = Array{complex_type}(undef, nz_el_count)

    Z_list = []

    @time for i in 1:nz_el_count
      if !(typeof(sp_input[3][i]) <: Function)
        sp_input_vals_eval[i] = sp_input[3][i]
      end
    end

    b = convert.(complex_type, b_real)
    for f in f_list
      verbose && @show f
      @time evaluate_matrix_for_w!(sp_input_vals_eval, sp_input[3], 2*pi*f)

      A_eval = sparse(sp_input[1], sp_input[2], sp_input_vals_eval)

      #return A_eval, b
      @time if iterative_solver
        x = bicgstabl(A_eval, b)
      else
        x = A_eval \ b
      end
      push!(Z_list, 1/x[1])
    end
  else
    params = get_prms_from_pairs(prms_pairs)
    header, sp_input, b_real = material_matrix_to_lin_sys(material_matrix, params)

    nz_el_count = length(sp_input[3])
    sp_input_vals_eval = Array{complex_type}(undef, nz_el_count)

    Z_list = []

    for i in 1:nz_el_count
      if !(typeof(sp_input[3][i]) <: Function)
        sp_input_vals_eval[i] = sp_input[3][i]
      end
    end

    b = convert.(complex_type, b_real)
    for f in f_list
      verbose && @show f        s
      evaluate_matrix_for_w!(sp_input_vals_eval, sp_input[3], 2*pi*f)

      A_eval = sparse(sp_input[1], sp_input[2], sp_input_vals_eval)

      #return A_eval, b
      if iterative_solver
        x = bicgstabl(A_eval, b)
      else
        x = A_eval \ b
      end
      push!(Z_list, 1/x[1])
    end
  end

  return f_list, Z_list
end

max_stars_repo_name: Masicko/ImageToEIS

max_issues_repo_name: Masicko/ImageToEIS

max_forks_repo_name: Masicko/ImageToEIS


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto
#see LICENSE.txt
#this work was supported in part by DARPA Contract D15PC00135
#unum_solver.jl

#takes a function and finds the zeroes.

function scalc{ESS,FSS}(x::Unum{ESS,FSS})
  #string(calculate(x))
  string(convert(Float64, x))
end

function scalc{ESS,FSS}(x::Ubound{ESS,FSS})
  #string(calculate(x.lowbound), "->", calculate(x.highbound))
  string(convert(Float64, x.lowbound), "->", convert(Float64, x.highbound))
end

function spans_zero(x, y)
  if isa(x, Unum)
    #analyzes the input unum and the output and decides if it should report that it crosses zero.
    if is_exact(x)  #exact unums must result in exactly zero to qualify.
      isa(y, Ubound) && return false
      return is_zero(y)  #a unum must be exactly zero to qualify as spanning zero.
    end
  end
  #inexact unums and ubounds must span zero; and basically have to result in ubounds.
  isa(y, Unum) && return is_zero(y)  #a unum result must be exactly zero to span zero
  is_negative(y.lowbound) && is_positive(y.highbound)
end

function spy1(v)
  print("testing: ", describe(v))
  v
end

function spy2(v)
  println(v ? " crosses" : " does not cross")
  v
end

#runs an analysis on the result and decides if more exponent is necessary.
function more_exponent_necessary(f, a)
  #first check to see if any of the values in the array warrant increasing the exponent.
  for (idx = 1:length(a))
    #is the value outright mmr or sss?
    (is_mmr(a[idx]) || is_sss(a[idx])) && return true
    #next check the function evaluated at this value.
    fv = f(a[idx])
    if isa(fv, Ubound)
      (is_sss(fv.lowbound) || is_sss(fv.highbound) || is_mmr(fv.lowbound) || is_mmr(fv.highbound)) && return true
    else
      (is_sss(fv) || is_mmr(fv)) && return true
    end
  end
  return false
end

function bignum_width(x::Unum)
  calculate(width(x))
end

function default_keep{ESS,FSS}(x::Unum{ESS,FSS})
  (ESS < 4) && return true
  is_mmr(x) && return false
  is_sss(x) && return false
  return true
end

function solve(f, lim, sess = 0, sfss = 0; verbose = false, kp = default_keep)
  #walk the number line in the starting environment.
  #compose the passed function with the spans_zero assesment to create the boolean-valued
  #analyzer function.

  #rj is a function that rejects a solution unum.  rj(u) == true means it will be
  #rejected.  Usually this should be used to set bounds on just how far the solver will go.
  if verbose
    bf(v) = spy2(spans_zero(v, f(spy1(v))))
  else
    bf(v) = spans_zero(v, f(v))
  end

  bw(v) = bitwalk(bf, v, true, true)

  res = fullwalk(bf, sess, sfss)

  function report()
    #if verbose
      println("passed solutions:")
      map((d) -> println(describe(d), " evaluates as ", describe(f(d)), " width ", bignum_width(d) ), res)
    #end
  end

  while (true)
    #first, check to see if we have any results.
    (length(res) == 0) && return res
    #check to see if we meet the termination condition.
    if more_exponent_necessary(f, res)
      verbose && println("overflow error in environment {$sess, $sfss}")
      report()

      res = mapreduce((d)->promote_ess(bf, d), vcat, [], res)
      filter!(kp, res)
      sess += 1
      verbose && println("investigating environment {$sess, $sfss}")
      continue
    end

    #calculate the worst width within the set of solution ulps.
    worst_width = mapreduce(bignum_width, max, 0.0, res)
    verbose && println("worst width across solution set: ", worst_width)
    (worst_width < lim) && break
    verbose && println("does not satisfy criterion < ", lim)

    #check to see if we are already at the maximum size, in which case, promote the fss.
    if res[1].fsize == max_fsize(sfss)
      sfss += 1
      verbose && println("promoting fss across the result set, environment {$sess, $sfss}")
      res = map(promote_fss, res)
      filter!(kp, res)
    end

    report()
    #apply bitwalker across the entire set.
    res = mapreduce(bw, vcat, [], res)
    filter!(kp, res)
  end
  res
end

export solve

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
function loadGroups(data::DataFrame, sample::EzXML.Node)
    groups = filter(group-> group["name"] ≠ "All Samples", findall("""//SampleRefs/SampleRef[@sampleID='$(sample["sampleID"])']/../..""",sample))
    return length(groups) ≠ 0 ? DataFrame([ (group["name"]=>fill(true,size(data,1))) for group ∈ groups ]) : loadGroups(data,nothing)
end

function loadGroups(data::DataFrame, sample::Nothing)
    return DataFrame("Ungrouped"=>fill(true,size(data,1)))
end


text:
using PyCall
using Conda

Conda.add("matplotlib")

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
# example assuming * operation is always binary

# ENV["JULIA_DEBUG"] = Metatheory

using Metatheory
using Metatheory.Library
using TermInterface

EGraphs.make(::Val{:numberfold}, g::EGraph, n::ENodeLiteral) = n.value

# This should be auto-generated by a macro
function EGraphs.make(::Val{:numberfold}, g::EGraph, n::ENodeTerm)
  if exprhead(n) == :call && arity(n) == 2
    op = operation(n)
    args = arguments(n)
    l = g[args[1]]
    r = g[args[2]]
    ldata = getdata(l, :numberfold, nothing)
    rdata = getdata(r, :numberfold, nothing)

    # @show ldata rdata

    if ldata isa Number && rdata isa Number
      if op == :*
        return ldata * rdata
      elseif op == :+
        return ldata + rdata
      end
    end
  end

  return nothing
end

function EGraphs.join(an::Val{:numberfold}, from, to)
  if from isa Number
    if to isa Number
      @assert from == to
    else
      return from
    end
  end
  return to
end

function EGraphs.modify!(::Val{:numberfold}, g::EGraph, id::Int64)
  eclass = g.classes[id]
  d = getdata(eclass, :numberfold, nothing)
  if d isa Number
    newclass, _ = addexpr!(g, d)
    merge!(g, newclass.id, id)
  end
end

EGraphs.islazy(::Val{:numberfold}) = false

comm_monoid = @theory begin
  ~a * ~b --> ~b * ~a
  ~a * 1 --> ~a
  ~a * (~b * ~c) --> (~a * ~b) * ~c
end

G = EGraph(:(3 * 4))
analyze!(G, :numberfold)

# exit(0)

@testset "Basic Constant Folding Example - Commutative Monoid" begin
  @test (true == @areequalg G comm_monoid 3 * 4 12)

  @test (true == @areequalg G comm_monoid 3 * 4 12 4 * 3 6 * 2)
end

@testset "Basic Constant Folding Example 2 - Commutative Monoid" begin
  ex = :(a * 3 * b * 4)
  G = EGraph(ex)
  analyze!(G, :numberfold)
  addexpr!(G, :(12 * a))
  @test (true == @areequalg G comm_monoid (12 * a) * b ((6 * 2) * b) * a)
  @test (true == @areequalg G comm_monoid (3 * a) * (4 * b) (12 * a) * b ((6 * 2) * b) * a)
end

@testset "Basic Constant Folding Example - Adding analysis after saturation" begin
  G = EGraph(:(3 * 4))
  # addexpr!(G, 12)
  saturate!(G, comm_monoid)
  addexpr!(G, :(a * 2))
  analyze!(G, :numberfold)
  saturate!(G, comm_monoid)

  @test (true == areequal(G, comm_monoid, :(3 * 4), 12, :(4 * 3), :(6 * 2)))

  ex = :(a * 3 * b * 4)
  G = EGraph(ex)
  analyze!(G, :numberfold)
  params = SaturationParams(timeout = 15)
  @test areequal(G, comm_monoid, :((3 * a) * (4 * b)), :((12 * a) * b), :(((6 * 2) * b) * a); params = params)
end

@testset "Infinite Loops analysis" begin
  boson = @theory begin
    1 * ~x --> ~x
  end

  G = EGraph(:(1 * x))
  params = SaturationParams(timeout = 100)
  saturate!(G, boson, params)
  ex = extract!(G, astsize)

  boson = @theory begin
    (:c * :cdag) --> :cdag * :c + 1
    ~a * (~b + ~c) --> (~a * ~b) + (~a * ~c)
    (~b + ~c) * ~a --> (~b * ~a) + (~c * ~a)
    # 1 * x => x
    (~a * ~b) * ~c --> ~a * (~b * ~c)
    ~a * (~b * ~c) --> (~a * ~b) * ~c
  end

  g = EGraph(:(c * c * cdag * cdag))
  saturate!(g, boson)
  ex = extract!(g, astsize_inv)

end

@testset "Extraction" begin
  comm_monoid = @commutative_monoid (*) 1

  fold_mul = @theory begin
    ~a::Number * ~b::Number => ~a * ~b
  end

  t = comm_monoid ∪ fold_mul

  @testset "Extraction 1 - Commutative Monoid" begin
    G = EGraph(:(3 * 4))
    saturate!(G, t)
    @test (12 == extract!(G, astsize))

    ex = :(a * 3 * b * 4)
    G = EGraph(ex)
    params = SaturationParams(timeout = 15)
    saturate!(G, t, params)
    extr = extract!(G, astsize)
    @test extr == :((12 * a) * b) ||
          extr == :(12 * (a * b)) ||
          extr == :(a * (b * 12)) ||
          extr == :((a * b) * 12) ||
          extr == :((12a) * b) ||
          extr == :(a * (12b)) ||
          extr == :((b * (12a))) ||
          extr == :((b * 12) * a) ||
          extr == :((b * a) * 12) ||
          extr == :(b * (a * 12)) ||
          extr == :((12b) * a)
  end

  fold_add = @theory begin
    ~a::Number + ~b::Number => ~a + ~b
  end

  @testset "Extraction 2" begin
    comm_group = @commutative_group (+) 0 inv

    t = comm_monoid ∪ comm_group ∪ (@distrib (*) (+)) ∪ fold_mul ∪ fold_add

    # for i ∈ 1:20
    # sleep(0.3)
    ex = :((x * (a + b)) + (y * (a + b)))
    G = EGraph(ex)
    saturate!(G, t)
    # end

    extract!(G, astsize) == :((x + y) * (b + a))
  end

  @testset "Lazy Extraction 2" begin
    comm_group = @commutative_group (+) 0 inv

    t = comm_monoid ∪ comm_group ∪ (@distrib (*) (+)) ∪ fold_mul ∪ fold_add

    # for i ∈ 1:20
    # sleep(0.3)
    ex = :((x * (a + b)) + (y * (a + b)))
    G = EGraph(ex)
    saturate!(G, t)
    # end

    extract!(G, astsize) == :((x + y) * (b + a))
  end

  @testset "Extraction - Adding analysis after saturation" begin
    G = EGraph(:(3 * 4))
    addexpr!(G, 12)
    saturate!(G, t)
    addexpr!(G, :(a * 2))
    saturate!(G, t)

    saturate!(G, t)

    @test (12 == extract!(G, astsize))

    # for i ∈ 1:100
    ex = :(a * 3 * b * 4)
    G = EGraph(ex)
    analyze!(G, :numberfold)
    params = SaturationParams(timeout = 15)
    saturate!(G, comm_monoid, params)

    extr = extract!(G, astsize)

    @test extr == :((12 * a) * b) ||
          extr == :(12 * (a * b)) ||
          extr == :(a * (b * 12)) ||
          extr == :((a * b) * 12) ||
          extr == :((12a) * b) ||
          extr == :(a * (12b)) ||
          extr == :((b * (12a))) ||
          extr == :((b * 12) * a) ||
          extr == :((b * a) * 12) ||
          extr == :(b * (a * 12))
  end

  comm_monoid = @commutative_monoid (*) 1

  comm_group = @commutative_group (+) 0 inv

  powers = @theory begin
    ~a * ~a → (~a)^2
    ~a → (~a)^1
    (~a)^~n * (~a)^~m → (~a)^(~n + ~m)
  end
  logids = @theory begin
    log((~a)^~n) --> ~n * log(~a)
    log(~x * ~y) --> log(~x) + log(~y)
    log(1) --> 0
    log(:e) --> 1
    :e^(log(~x)) --> ~x
  end

  t = comm_monoid ∪ comm_group ∪ (@distrib (*) (+)) ∪ powers ∪ logids ∪ fold_mul ∪ fold_add

  @testset "Complex Extraction" begin
    G = EGraph(:(log(e) * log(e)))
    params = SaturationParams(timeout = 8)
    saturate!(G, t, params)
    @test extract!(G, astsize) == 1

    G = EGraph(:(log(e) * (log(e) * e^(log(3)))))
    params = SaturationParams(timeout = 7)
    saturate!(G, t, params)
    @test extract!(G, astsize) == 3

    @time begin
      G = EGraph(:(a^3 * a^2))
      saturate!(G, t)
      ex = extract!(G, astsize)
    end
    @test ex == :(a^5)

    @time begin
      G = EGraph(:(a^3 * a^2))
      saturate!(G, t)
      ex = extract!(G, astsize)
    end
    @test ex == :(a^5)

    function cust_astsize(n::ENodeTerm, g::EGraph)
      cost = 1 + arity(n)

      if operation(n) == :^
        cost += 2
      end

      for id in arguments(n)
        eclass = g[id]
        !hasdata(eclass, cust_astsize) && (cost += Inf; break)
        cost += last(getdata(eclass, cust_astsize))
      end
      return cost
    end

    cust_astsize(n::ENodeLiteral, g::EGraph) = 1

    @time begin
      G = EGraph(:((log(e) * log(e)) * (log(a^3 * a^2))))
      saturate!(G, t)
      @show getcost!(G, cust_astsize)
      ex = extract!(G, cust_astsize)
    end
    @test ex == :(5 * log(a)) || ex == :(log(a) * 5)
  end

  function costfun(n::ENodeTerm, g::EGraph)
    arity(n) != 2 && (return 1)
    left = arguments(n)[1]
    left_class = g[left]
    ENodeLiteral(:a) ∈ left_class.nodes ? 1 : 100
  end

  costfun(n::ENodeLiteral, g::EGraph) = 1

  moveright = @theory begin
    (:b * (:a * ~c)) --> (:a * (:b * ~c))
  end

  expr = :(a * (a * (b * (a * b))))
  res = rewrite(expr, moveright)

  g = EGraph(expr)
  saturate!(g, moveright)
  resg = extract!(g, costfun)

  @testset "Symbols in Right hand" begin
    @test resg == res == :(a * (a * (a * (b * b))))
  end

  co = @theory begin
    foo(~x ⋅ :bazoo ⋅ :woo) --> Σ(:n * ~x)
  end
  @testset "Consistency with classical backend" begin
    ex = :(foo(wa(rio) ⋅ bazoo ⋅ woo))
    g = EGraph(ex)
    saturate!(g, co)

    res = extract!(g, astsize)

    resclassic = rewrite(ex, co)

    @test res == resclassic
  end

  @testset "No arguments" begin
    ex = :(f())
    g = EGraph(ex)
    @test :(f()) == extract!(g, astsize)

    ex = :(f() + g())

    t = @theory begin
      f() + g() --> h()
    end

    gg = EGraph(ex)
    saturate!(gg, t)
    @show getcost!(gg, astsize)
    res = extract!(gg, astsize)

    @test res == :(h())
  end

  @testset "Symbol or function object operators in expressions in EGraphs" begin
    ex = :(($+)(x, y))
    t = [@rule a b a + b => 2]
    g = EGraph(ex)
    saturate!(g, t)
    @test extract!(g, astsize) == 2
  end
end


text:
# =====================================
# Flowpipe composition with a lazy map
# =====================================

"""
    MappedFlowpipe{FT<:AbstractFlowpipe, ST} <: AbstractFlowpipe

### Fields

- `F`    -- flowpipe
- `func` -- function representing the map
"""
struct MappedFlowpipe{FT<:AbstractFlowpipe, ST} <: AbstractFlowpipe
    F::FT
    func::ST
end

"""
    Projection(fp::AbstractFlowpipe, vars::NTuple{D, T}) where {D, T<:Integer}

Return the lazy projection of a flowpipe.

### Input

### Output

### Notes

The projection is lazy, and consists of mapping each set
`X` in the flowpipe to `MX`, where `M` is the projection matrix associated with
the given variables `vars`.
"""
function LazySets.Projection(fp::AbstractFlowpipe, vars::NTuple{D, T}) where {D, T<:Integer}
    # TODO: assert that vars belong to the variables of the flowpipe
    M = projection_matrix(collect(vars), dim(F), Float64)
    func = @map(x -> M*x)
    return MappedFlowpipe(fp, func)
end

function overapproximate(fp::Flowpipe, args...)
    return Flowpipe(map(R -> overapproximate(R, args...), fp), fp.ext)
end

function overapproximate(fp::VRT, args...) where {RT, VRT<:AbstractVector{RT}}
    return Flowpipe(map(R -> overapproximate(R, args...), fp))
end


text:
function _make_thres_gui()

    frame1_3=Frame("Threshold")
    vbox1_3_1=Box(:v)
    push!(frame1_3,vbox1_3_1)

    #sb=SpinButton(-300:300)
    #setproperty!(sb,:value,0)
    sb=Label("0")
    push!(vbox1_3_1,sb)

    button_thres_all = CheckButton()
    add_button_label(button_thres_all,"All Channels")
    Gtk.GAccessor.active(button_thres_all,false)
    push!(vbox1_3_1,button_thres_all)

    button_thres = CheckButton()
    add_button_label(button_thres,"Show")
    Gtk.GAccessor.active(button_thres,false)
    push!(vbox1_3_1,button_thres)

    vbox_slider=Box(:v)
    thres_slider = Scale(true, -300,300,1)
    adj_thres = Adjustment(thres_slider)
    Gtk.GAccessor.value(adj_thres,0)

    c_thres=Canvas(10,200)
    Gtk.GAccessor.vexpand(c_thres,false)

    Gtk.GAccessor.inverted(thres_slider,true)
    Gtk.GAccessor.draw_value(thres_slider,false)

    Gtk.GAccessor.vexpand(thres_slider,true)
    push!(vbox_slider,thres_slider)
    push!(vbox_slider,c_thres)

    thres_widgets=SpikeSorting.Thres_Widgets(sb,thres_slider,adj_thres,button_thres_all,button_thres)

    (thres_widgets,frame1_3,vbox_slider)
end

function add_thres_cb(sc_widgets)

    signal_connect(SpikeSorting.thres_cb,sc_widgets.thres_slider,"value-changed",Void,(),false,(sc_widgets,))
    signal_connect(SpikeSorting.thres_show_cb,sc_widgets.thres_widgets.show,"clicked",Void,(),false,(sc_widgets,))

    nothing
end

#=
Set Threshold for sorting equal to GUI threshold
=#
function thres_changed(han::Gui_Handles,s,fpga,backup)

    mythres=getproperty(han.sc.adj_thres,:value,Int)
    han.sc.thres=mythres

    update_thres(han,s,backup)

    send_thres_to_ic(han,fpga)

    han.sc.thres_changed=false

    nothing
end

#=
Functions used to update Sorting data structure with threshold from slider
=#
function update_thres(han::Gui_Handles,s::Array,backup)
    if (getproperty(han.sc.thres_widgets.all,:active,Bool))|(getproperty(han.sc.gain_widgets.all,:active,Bool))
        @inbounds for i=1:length(s)
            s[i].thres=-1*han.sc.thres/han.scale[i,1]+han.offset[i]
            f=open(string(backup,"thres/",i,".bin"),"w")
            write(f,s[i].thres)
            close(f)
        end
    else
        @inbounds s[han.sc.spike].thres=-1*han.sc.thres/han.sc.s+han.sc.o
        f=open(string(backup,"thres/",han.sc.spike,".bin"),"w")
        write(f,s[han.sc.spike].thres)
        close(f)
    end
end

function update_thres(han::Gui_Handles,s::DArray{T,1,Array{T,1}},backup) where T
    if (getproperty(han.sc.thres_widgets.all,:active,Bool))|(getproperty(han.sc.gain_widgets.all,:active,Bool))
        @sync begin
            for p in procs(s)
                @async remotecall_wait((ss)->set_multiple_thres(localpart(ss),han,localindexes(ss)),p,s)
            end
        end
        for i=1:size(han.scale,1)
            f=open(string(backup,"thres/",i,".bin"),"w")
            write(f,-1*han.sc.thres/han.scale[i,1]+h.offset[i])
            close(f)
        end
    else
        (nn,mycore)=get_thres_id(s,han.sc.spike)
        remotecall_wait(((x,h,num)->remote_set_thres(localpart(x)[num],h)),mycore,s,han,nn)
        f=open(string(backup,"thres/",han.sc.spike,".bin"),"w")
        write(f,-1*han.sc.thres/han.sc.s+han.sc.o)
        close(f)
    end
end

function remote_set_thres(x,h)

    x.thres = -1 * h.sc.thres/h.sc.s + h.sc.o

    nothing
end

function get_thres_id(s::DArray{T,1,Array{T,1}},ss::Int64) where T

    mycore=1
    mynum=ss

    for i=2:(length(procs(s))+1)
        if ss<s.cuts[1][i]
            mycore=i
            mynum=ss-s.cuts[1][i-1]+1
            break
        end
    end
    (mynum,mycore)
end

function set_multiple_thres(s::Array,han::Gui_Handles,inds)
    @inbounds for i=1:length(s)
        s[i].thres=-1*han.sc.thres/han.scale[inds[1][i],1]+han.offset[inds[1][i]]
    end
end

#=
Get threshold from Sorting data structure and set threshold in GUI handles equal to it.
=#

function get_thres(han::Gui_Handles,s::DArray{T,1,Array{T,1}}) where T<:Sorting
    (nn,mycore)=get_thres_id(s,han.sc.spike)

    mythres=remotecall_fetch(((x,h,num)->(localpart(x)[num].thres-h.sc.o)*h.sc.s*-1),mycore,s,han,nn)

    setproperty!(han.sc.adj_thres,:value,round(Int64,mythres)) #show threshold

    nothing
end

function get_thres(han::Gui_Handles,s::Array{T,1}) where T<:Sorting

    mythres=(s[han.sc.spike].thres-han.sc.o)*han.sc.s*-1
    setproperty!(han.sc.adj_thres,:value,round(Int64,mythres)) #show threshold

    nothing
end

#=
Threshold Callbacks
=#

function sb2_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han, rhd = user_data

    mygain=getproperty(han.sc.gain_widgets.all,:active,Bool)

    gainval=getproperty(han.sc.gain_widgets.gainbox,:value,Int)
    mythres=getproperty(han.sc.adj_thres,:value,Int)

    if mygain==true
        han.scale[:,1]=-1 .* gainval/1000
        han.scale[:,2]=-.2 * gainval/1000

        for i=1:size(han.scale,1)
            f=open(string(rhd.save.backup,"gain/",i,".bin"),"w")
            write(f,han.scale[i,1])
            close(f)
        end
    else
        han.scale[han.sc.spike,1]=-1 * gainval/1000
        han.scale[han.sc.spike,2]=-.2 * gainval/1000

        f=open(string(rhd.save.backup,"gain/",han.sc.spike,".bin"),"w")
        write(f,han.scale[han.sc.spike,1])
        close(f)
    end

    han.sc.s = han.scale[han.sc.spike,1]

    han.sc.thres_changed=true

    nothing
end


text:
# module for calculating particle motion
# considering shorter distance observation by occlusion

using Distributions, LinearAlgebra

include(joinpath(split(@__FILE__, "src")[1], "src/common/state_transition/state_transition.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/common/observation_function/observation_function.jl"))

mutable struct OcclusionFreeParticle
  pose
  weight

  # init
  function OcclusionFreeParticle(init_pose::Array, weight::Float64)
    self = new()
    self.pose = init_pose
    self.weight = weight
    return self
  end
end

function motion_update(self::OcclusionFreeParticle, speed,
                       yaw_rate, time_interval, noise_rate_pdf)
  ns = rand(noise_rate_pdf) # [nn, no, on, oo]
  noised_spd = speed + ns[1]*sqrt(abs(speed)/time_interval) + ns[2]*sqrt(abs(yaw_rate)/time_interval)
  noised_yr = yaw_rate + ns[3]*sqrt(abs(speed)/time_interval) + ns[4]*sqrt(abs(yaw_rate)/time_interval)
  self.pose = state_transition(noised_spd, noised_yr, time_interval, self.pose)
end

function observation_update(self::OcclusionFreeParticle, observation,
                            env_map, dist_dev_rate, dir_dev)
  for obs in observation
    obs_pose = obs[1] # [distance, direction]
    obs_id = obs[2]

    # calculate distance and direction from particle pos to object
    pos_on_map = env_map.objects[obs_id].pose
    obs_from_particle = observation_function(self.pose, pos_on_map) # [distance, direction]

    # calculate likelihood
    dist_dev = dist_dev_rate * obs_from_particle[1]
    cov = diagm(0 => [dist_dev^2, dir_dev^2])
    if obs_pose[1] < obs_from_particle[1]
      obs_pose[1] = obs_from_particle[1]
    end
    self.weight *= pdf(MvNormal(obs_from_particle, cov), obs_pose)
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text: Dict{String, Any}("local" => Dict{String, Any}("value" => "1987-07-05T17:45:00", "type" => "datetime-local"), "milli" => Dict{String, Any}("value" => "1977-12-21T10:32:00.555", "type" => "datetime-local"), "space" => Dict{String, Any}("value" => "1987-07-05T17:45:00", "type" => "datetime-local"))


text:
using CImGui
using CImGui.CSyntax
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using Printf
using ImageCore

function construct_B()
    B = repeat(1:500, 1, 401)
    low = colorant"navyblue";
    medium = colorant"limegreen";
    high = colorant"yellow";
    Bcolormap = colorsigned(low,medium,high) ∘ scalesigned(minimum(B),(minimum(B)+maximum(B)) / 2, maximum(B))
    Bcolormap.(B)
end

function example(B::AbstractArray)
    @static if Sys.isapple()
        # OpenGL 3.2 + GLSL 150
        glsl_version = 150
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
    else
        # OpenGL 3.0 + GLSL 130
        glsl_version = 130
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
        # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
    end

    # setup GLFW error callback
    error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
    GLFW.SetErrorCallback(error_callback)

    # create window
    window = GLFW.CreateWindow(1280, 720, "Demo")
    @assert window != C_NULL
    GLFW.MakeContextCurrent(window)
    GLFW.SwapInterval(1)  # enable vsync

    # setup Dear ImGui context
    ctx = CImGui.CreateContext()

    # setup Dear ImGui style
    CImGui.StyleColorsDark()

    # create texture for image drawing
    img₀ = B
    img₁= transpose(img₀)
    img_width, img_height = size(img₁)
    glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
    image_id = ImGui_ImplOpenGL3_CreateImageTexture(img_width, img_height, format = GL_RGBA)
    img′ = unsafe_wrap(Array{UInt8,3}, convert(Ptr{UInt8}, pointer(img₁)), (Cint(3), Cint(img_width), Cint(img_height)))
    ImGui_ImplOpenGL3_UpdateImageTexture(image_id, img′, img_width, img_height; format = GL_RGB)

    # setup Platform/Renderer bindings
    ImGui_ImplGlfw_InitForOpenGL(window, true)
    ImGui_ImplOpenGL3_Init(glsl_version)

    demo_open = true
    clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]
    while !GLFW.WindowShouldClose(window)
        demo_open # oh my global scope
        image_id, img_width, img_height

        GLFW.PollEvents()
        # start the Dear ImGui frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        CImGui.NewFrame()

        CImGui.Begin("Image Demo")
            CImGui.Image(Ptr{Cvoid}(image_id), (img_width, img_height))
        CImGui.End()

        # rendering
        CImGui.Render()
        GLFW.MakeContextCurrent(window)
        display_w, display_h = GLFW.GetFramebufferSize(window)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

        GLFW.MakeContextCurrent(window)
        GLFW.SwapBuffers(window)
    end

    # cleanup
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    CImGui.DestroyContext(ctx)
    GLFW.DestroyWindow(window)
end

# Run the demo
example(construct_B())


text:
"""Describes how to construct forecasts from raw timeseries data files."""
mutable struct TimeseriesFileMetadata
    simulation::String  # User description of simulation
    category::String  # String version of abstract type for the forecasted component.
                      # Calling module should determine the actual type.
    component_name::String  # Name of forecast component
    label::String  # Accessor function on component for source of timeseries
    scaling_factor::Union{String, Float64}  # Controls normalization of timeseries.
                                            # Use 1.0 for pre-normalized data.
                                            # Use 'Max' to divide the timeseries by the max
                                            #   value in the column.
                                            # Use any float for a custom scaling factor.
    data_file::String  # path to the timeseries data file
    percentiles::Vector{Float64}
    forecast_type::String
    component::Union{Nothing, InfrastructureSystemsType}  # Calling module must set.
end

function TimeseriesFileMetadata(simulation, category, component_name, label, scaling_factor,
                                data_file, percentiles, forecast_type)
    return TimeseriesFileMetadata(
        simulation, category, component_name, label, scaling_factor, data_file, percentiles,
        forecast_type, nothing,
    )
end

"""Reads forecast metadata and fixes relative paths to the data files."""
function read_time_series_metadata(file_path::AbstractString)
    if endswith(file_path, ".json")
        metadata = open(file_path) do io
            metadata = Vector{TimeseriesFileMetadata}()
            data = JSON.parse(io)
            for item in data
                category = _get_category(item["category"])
                scaling_factor = item["scaling_factor"]
                if !isa(scaling_factor, AbstractString)
                    scaling_factor = Float64(scaling_factor)
                end
                push!(metadata, TimeseriesFileMetadata(
                    item["simulation"],
                    item["category"],
                    item["component_name"],
                    item["label"],
                    scaling_factor,
                    item["data_file"],
                    # Use default values until CDM data is updated.
                    get(item, "percentiles", []),
                    get(item, "forecast_type", "DeterministicInternal"),
                ))
            end
            return metadata
        end
    elseif endswith(file_path, ".csv")
        csv = CSV.read(file_path)
        metadata = Vector{TimeseriesFileMetadata}()
        for row in eachrow(csv)
            category = _get_category(row.Category)
            push!(metadata, TimeseriesFileMetadata(row.Simulation,
                                                   row.Category,
                                                   row.Object,
                                                   row.Parameter,
                                                   row[Symbol("Scaling Factor")],
                                                   row[Symbol("Data File")],
                                                   # TODO: update CDM data for the next
                                                   # two fields.
                                                   [],
                                                   "DeterministicInternal",
                                                  )
            )
        end

    else
        error("file not supported")
    end

    directory = dirname(file_path)
    for ts_metadata in metadata
        ts_metadata.data_file = abspath(joinpath(directory, ts_metadata.data_file))
    end

    return metadata
end

function _get_category(category::String)
    # HACK alert because of wackiness in PowerSystems RTS data.
    if category == "LoadZone"
        category = "bus"
    end

    return lowercase(category)
end

"""
    set_component!(
                   metadata::TimeseriesFileMetadata,
                   components::Components,
                   mod::Module,
                  )

Set the component value in metadata by looking up the category in module.
This requires that category be a string version of a component's abstract type.
Modules can override for custom behavior.
"""
function set_component!(
                        metadata::TimeseriesFileMetadata,
                        components::Components,
                        mod::Module,
                       )

    # TODO: CDM data should change LoadZone to LoadZones.
    symbol = metadata.category == "LoadZone" ? :LoadZones : Symbol(metadata.category)
    category = getfield(mod, symbol)
    if isconcretetype(category)
        metadata.component = get_component(category, components, metadata.component_name)
        if isnothing(metadata.component)
            throw(DataFormatError(
                "no component category=$category name=$(metadata.component_name)"
            ))
        end
    else
        components = get_components_by_name(category, components, metadata.component_name)
        if length(components) == 0
            @warn "no component category=$category name=$(metadata.component_name)"
            metadata.component = nothing
        elseif length(components) == 1
            metadata.component = components[1]
        else
            throw(DataFormatError(
                "duplicate names type=$(category) name=$(metadata.component_name)"
            ))
        end
    end
end

struct ForecastInfo
    simulation::String
    component::InfrastructureSystemsType
    label::String  # Component field on which timeseries data is based.
    scaling_factor::Union{String, Float64}
    data::TimeSeries.TimeArray
    percentiles::Vector{Float64}
    file_path::String
    forecast_type::String

    function ForecastInfo(simulation, component, label, scaling_factor, data, percentiles,
                          file_path, forecast_type)
        new(simulation, component, label, scaling_factor, data, percentiles,
            abspath(file_path), forecast_type)
    end
end

function ForecastInfo(metadata::TimeseriesFileMetadata,
                      timeseries::TimeSeries.TimeArray)
    return ForecastInfo(metadata.simulation, metadata.component, metadata.label,
                        metadata.scaling_factor, timeseries, metadata.percentiles,
                        metadata.data_file, metadata.forecast_type)
end

function get_forecast_type(forecast_info::ForecastInfo)
    return getfield(InfrastructureSystems, Symbol(forecast_info.forecast_type))
end

struct ForecastCache
    forecasts::Vector{ForecastInfo}
    data_files::Dict{String, TimeSeries.TimeArray}
end

function ForecastCache()
    return ForecastCache(Vector{ForecastInfo}(),
                         Dict{String, TimeSeries.TimeArray}())
end

function handle_scaling_factor(timeseries::TimeSeries.TimeArray,
                               scaling_factor::Union{String, Float64})
    if scaling_factor isa String
        if lowercase(scaling_factor) == "max"
            max_value = maximum(TimeSeries.values(timeseries))
            timeseries = timeseries ./ max_value
            @debug "Normalize by max value" max_value
        else
            throw(DataFormatError("invalid scaling_factor=scaling_factor"))
        end
    elseif scaling_factor != 1.0
        timeseries = timeseries ./ scaling_factor
        @debug "Normalize by custom scaling factor" scaling_factor
    else
        @debug "forecast is already normalized"
    end

    return timeseries
end

function _add_forecast_info!(forecast_cache::ForecastCache, data_file::AbstractString,
                             component_name::Union{Nothing, String})
    if !haskey(forecast_cache.data_files, data_file)
        forecast_cache.data_files[data_file] = read_time_series(data_file, component_name)
        @debug "Added timeseries file" data_file
    end

    return forecast_cache.data_files[data_file]
end


text:
#=
Copyright 2017 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

# 1D blood flow within elastic vessels is ruled by a system of hyperbolic
# partial derivative equations
# $$
#   \begin{cases}
#     \dfrac{\partial A}{\partial t} + \dfrac{\partial Au}{\partial x} = 0, \\
#     \dfrac{\partial Au}{\partial t} + \dfrac{\partial Au^2}{\partial x}
#       \dfrac{A}{\rho}\dfrac{\partial P}{\partial x} = - Ru, \\
#     P(x, t) = P_e(x, t) + K_a \left[\left(\dfrac{A}{A_0} \right)^{1/2}
#     -1 \right],
#   \end{cases}
# $$
# By defining $\mathbf{U} = [A, u]^T$, the conservative form of the system
# can be written along with initial and boundary conditions as
# $$
#   \begin{cases}
#     \mathbf{U}_t + \mathbf{F}(\mathbf{U})_x = \mathbf{S}(\mathbf{U}),
#       \quad x \in \big[0, L \big], \\
#     \mathbf{U}(x, 0) = \mathbf{U}^{(0)}(x), \quad t>0, \\
#     \mathbf{U}(0, t) = \mathbf{U}_l(t), \mathbf{U}(L, t) = \mathbf{U}_r(t),
#   \end{cases}
# $$
# where $\mathbf{F}$ is the *flux* term, and $\mathbf{S}$ is the *source*
# term. This set of equations and initial and
# boundary conditions is also called Initial Boundary Value Problem (IBVP).
# The IBVP is solved by means of the Godunov's method [^1], a
# conservative first order finite volume method. The method is applied on
# a computational domain (red) subdivided in $M$ *cell*s $I_i$ (gray).

#

#

# Domain length is
# equal to vessel length $L$, and cell length is given by $\Delta x$ size.
# The total simulation time $T$ defines domain height, hence cell height is
# given by $\Delta t$. In each cell $I_i$ three points are defined: center point $x_i$, left
# extreme $x_{i-1/2}$, and right extreme $x_{i+1/2}$.

#

#

# Godunov's method assumes a piece-wise constant distribution of $\mathbf{U}$
# at each $t = t^n$.

#

#

# For a *small enough* $\Delta t$, the method updates the solution as
# $$
#   \mathbf{U}_i^{n+1} = \mathbf{U}_i^{n} + \frac{\Delta t}{\Delta x}
#     \left(\mathbf{F}_{i-1/2} - \mathbf{F}_{i+1/2} \right),
# $$
# where $\mathbf{F}_{i\pm 1/2} = \mathbf{F}\left(\mathbf{U}_{i\pm 1/2}(0)
# \right)$ is the *numerical* flux. The argument of the numerical flux,
# $\mathbf{U}_{i\pm 1/2}(0)$, is computed as the solution of a homogeneous (
# the source term will be added later on) quasi-linear local Riemann
# problem (RP) at the interface of two cells. For $\mathbf{U}_{i + 1/2}(0)$
# we have
# $$
#   \begin{cases}
#     \mathbf{U}_t + \mathbf{F}_\mathbf{U} \mathbf{U}_x = 0, \\
#     \mathbf{U}(x_i, 0) = \mathbf{U}^{(0)}(x_i) = \begin{cases}
#       \mathbf{U}_L = \mathbf{U}(x_{i-1/2}), \quad x<0, \\
#       \mathbf{U}_R = \mathbf{U}(x_{i+1/2}), \quad x>0. \\
#     \end{cases}
#   \end{cases}
# $$
# The solution of the RP is depicted below.

#

#

# The solution structure includes two families of waves associated to the
# two system Riemann invariants $\lambda_1$ and $\lambda_2$, respectively.
# There are three constant states $\mathbf{U}_L$, $\mathbf{U}^*$, and
# $\mathbf{U}_R$, of which $\mathbf{U}^*$ is unknown and refers
# to the solution in the star region (yellow). We need *jump*
# conditions to cross the waves and connect these states. Depending on the
# nature of each wave we have *shock* or *rarefaction* conditions.

# [`Godunov`](godunov.html#Godunov) function glues together all the functions
# needed to [solve the RP](godunov.html#riemannSolver)
# at each interface and to compute the [numerical flux](godunov.html#flux).

# *function* __`Godunov`__
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `i`              `::Int` cell index within the vessel numerical domain.
#
# `v`              `::Vessel` current vessel data structure.
#
# `dt`             `::Float` $\Delta t$, current time step.
#
# `b`              `::Blood` blood data structure.
# ----------------------------------------------------------------------------
#
# --------------------------------------------------------------------------
# Functioning:
# --------------------------------------------------------------------------
# In order to solve the local RP, local boundary conditions should be
# retrieved at both left ($i-1/2$) and right ($i+1/2$) interfaces.
# --------------------------------------------------------------------------
#
function Godunov(i :: Int64, v :: Vessel, dt :: Float64, b :: Blood)

  #

  #

  # First we compute the states $\mathbf{U}_L$ and $\mathbf{U}_R$ at the
  # left hand side interface ($i-1/2$). When we deal with the inlet node
  # ($i=1$), blue ghost cells are used to find left states. Right states are
  # always inside the numerical domain and they do not need ghost cells.
  #i-1/2 local left
  if i-1 == 0      #1st ghost cell
    UlA = v.U00A
    UlQ = v.U00Q

  elseif i-1 == -1 #2nd ghost cell
    UlA = v.U01A
    UlQ = v.U01Q

  else
    UlA = v.A[i-1]
    UlQ = v.Q[i-1]
  end

  #local right
  UrA = v.A[i]
  UrQ = v.Q[i]
  # Once the local states have been computed, the solution of the RP at the
  # left
  # interface is found by [`riemannSolver`](godunov.html#riemannSolver)
  # function.
  Al, Ql = riemannSolver(UlA, UlQ, UrA, UrQ, v)

  # The same process is applied to the local right interface.
  #i+1/2 local right
  if i+1 == v.M+1     #1st ghost cell
    UrA = v.UM1A
    UrQ = v.UM1Q

  elseif i+1 == v.M+2 #2nd ghost cell
    UrA = v.UM2A
    UrQ = v.UM2Q

  else
    UrA = v.A[i+1]
    UrQ = v.Q[i+1]
  end

  #local left
  UlA = v.A[i]
  UlQ = v.Q[i]

  Ar, Qr = riemannSolver(UlA, UlQ, UrA, UrQ, v)

  # When both interfaces are solved, left and right ($\mathbf{F}_{i\mp 1/2}$)
  # numerical fluxes are computed by [`flux`](godunov.html#flux) fucntion.
  Fl1, Fl2 = flux(Al, Ql, v.gamma)
  Fr1, Fr2 = flux(Ar, Qr, v.gamma)

  # Conservative variables `A` and `Q` are updated by Godunov's method
  # formula.
  v.A[i] += dt/v.dx * (Fl1-Fr1)
  v.Q[i] += dt/v.dx * (Fl2-Fr2)

  # The source term is solved by [`source`](godunov.html#source) function
  # and `Q` is directly updated. The source term is always zero for the
  # continuity equation, hence `A` is not updated.
  v.Q[i] += source(dt, v.A[i], v.u[i], b)

  # Eventually, primitive variables `P`, `u`, and `c` are updated accordingly
  # to their definitions by [`pressure`](converter.html#pressure) and
  # [`waveSpeed`](converter.html#waveSpeed) functions.
  v.P[i] = pressure(v.A[i], v.A0, v.beta, v.Pext)
  v.u[i] = v.Q[i]/v.A[i]
  v.c[i] = waveSpeed(v.A[i], v.gamma)

end

# The [RP solution](godunov.html#RP_solution) is herein coded. The Riemann
# solver computes the cross-sectional area in the star region
# $A^*$ starting from a guess solution (two rarefactions solution). The
# value of $A^*$ comes from the solution of the nonlinear equation
# $$
#   f(A^*) = f_L(A^*) + f_R(A^*) + \Delta u = 0,
# $$
# where $\Delta u = u_R - u_L$, and $f_L$ and $f_R$ are the local fluxes
# computed by [`fl`](godunov.html#fl) and [`fr`](godunov.html#fr),
# respectively.

# *function* __`riemannSolver`__ $\rightarrow$ `As::Float`, `Qs::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `Al`             `::Float` `A` value at the left side of the local
#                  cell interface.
#
# `Ql`             `::Float` `Q` value at the left side of the local
#                  cell interface.
#
# `Ar`             `::Float` `A` value at the right side of the local
#                  cell interface.
#
# `Qr`             `::Float` `Q` value at the right side of the local
#                  cell interface.
#
# `v`              `::Vessel` current vessel data structure.
# ----------------------------------------------------------------------------
#
function riemannSolver(Al :: Float64, Ql :: Float64,
                       Ar :: Float64, Qr :: Float64,
                        v :: Vessel)
  # --------------------------------------------------------------------------
  # Functioning:
  # --------------------------------------------------------------------------
  # Primitives are initialised by their definitions, and the two rarefaction
  # solution is applied to compute initial values for $c^*$ and $A^*$.
  # --------------------------------------------------------------------------
  #initialise solution
  ur = Qr/Ar
  ul = Ql/Al
  cl = waveSpeed(Al, v.gamma)
  cr = waveSpeed(Ar, v.gamma)
  du = ur - ul

#     #vacuum
#   if ul + 4*cl <= ur - 4*cr
# #     println("v")
# #     return (Al+Ar)*0.5, (Ql+Qr)*0.5
#     return v.A0, (Ql+Qr)*0.5
#   end

  #two rarefaction solution
  cs = 0.5 * (cl + cr) - (ur - ul)*0.125
  As = 4/9 * (cs^4)/(v.gamma*v.gamma)
  # The nonlinear equation for the RP is defined and solved by Newton's
  # method.
  fun(As) = fL(As, Al, cs, cl, v) + fR(As, Ar, cs, cr, v) + du

  As = newton(fun, As)

  # Eventually, $u^*$ is computed as
  # $$
  #   \frac{1}{2}\left(u_L + u_R\right) + \frac{1}{2}\left[f_R\left(A^*\right)
  #     -f_L\left(A^* \right) \right],
  # $$
  us = 0.5*((ur+ul) + (fR(As, Ar, cs, cr, v) - fL(As, Al, cs, cl, v)))
  # $Q^*$ is calculated from its definition and returned.
  Qs = As * us
  # --------------------------------------------------------------------------
  # Returns:
  # ---------- ---------------------------------------------------------------
  # `As`       `::Float` cross-sectional area in the star region, solution
  #            of the RP at the interface.
  #
  # `Qs`       `::Float` flow rate in the star region, solution
  #            of the RP at the interface.
  # --------------------------------------------------------------------------
  return As, Qs
end

# Local fluxes (left and right) are computed depending on the nature of the
# wave system located at the interface. In particular
# $$
#   f_L = \begin{cases}
#       4\left(c^* - c_L\right), & A^*\leq A_L \quad \text{(rarefaction)}, \\
#       \sqrt{\dfrac{\gamma\left(A^*-A_L\right)\left(A^{* 3/2}-
#         A_L^{3/2}\right)}{A_L A^*}}, & A^*>A_L \quad \text{(shock)},
#     \end{cases}
# $$
# $$
#   f_R = \begin{cases}
#       4\left(c^* - c_R\right), & A^*\leq A_R \quad \text{(rarefaction)}, \\
#       \sqrt{\dfrac{\gamma\left(A^*-A_R\right)\left(A^{* 3/2}-
#         A_R^{3/2}\right)}{A_R A^*}}, & A^*>A_R \quad \text{(shock)},
#     \end{cases}
# $$

# *function* __`fL`__ $\rightarrow$ `left_flux::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `A`           `::Float` current value of $A^*$ at the left interface.
#
# `Al`          `::Float` cross sectional area of the left cell.
#
# `c`           `::Float` current value of $c^*$ at the left interface.
#
# `cl`          `::Float` wave speed in the left cell.
#
# `v`           `::Vessel` current vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ------------- --------------------------------------------------------------
# `left_flux`   `::Float` local flux at the left interface.
# ----------------------------------------------------------------------------
#
function fL(A, Al :: Float64, c :: Float64, cl :: Float64, v :: Vessel)

  if A <= Al
    return 4*(c - cl)

  else
    return sqrt(
                (v.gamma*(A - Al)*(A^1.5 - Al^1.5)) /
                              (Al*A)
               )
  end

end

# *function* __`fR`__ $\rightarrow$ `right_flux::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `A`           `::Float` current value of $A^*$ at the right interface.
#
# `Ar`          `::Float` cross sectional area of the right cell.
#
# `c`           `::Float` current value of $c^*$ at the right interface.
#
# `cr`          `::Float` wave speed in the right cell.
#
# `v`           `::Vessel` current vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ------------- --------------------------------------------------------------
# `right_flux`   `::Float` local flux at the right interface.
# ----------------------------------------------------------------------------
#
function fR(A, Ar :: Float64, c :: Float64, cr :: Float64, v :: Vessel)

  if A <= Ar
    return 4*(c - cr)

  else
    return sqrt(
                (v.gamma*(A - Ar)*(A^1.5 - Ar^1.5)) /
                              (Ar*A)
               )
  end

end

# The flux in conservative form reads
# $$
#   \mathbf{F}(\mathbf{U}) = \left\{ \begin{array}{c}
#       Au \\
#       Au^2 + \gamma A^{3/2}
#     \end{array} \right\},
# $$
# and its numerical form is computed by [`flux`](godunov.html#flux) function.

# *function* __`flux`__ $\rightarrow$ `f1::Float`, `f2::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `A`           `::Float` cross sectional area along the cell.
#
# `Q`           `::Float` flow rate within the cell.
#
# `gamma`       `::Vessel` elastic constant for `c` calculation.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ---------- -----------------------------------------------------------------
# `f1`       `::Float` first component of the numerical flux.
#
# `f2`       `::Float` second component of the numerical flux.
# ----------------------------------------------------------------------------
#
function flux(A :: Float64, Q :: Float64, gamma :: Float64)

  f1 = Q
  f2 = Q*Q/A + gamma*A^1.5

  return f1, f2
end

# *function* __`source`__ $\rightarrow$ `S::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------ ----------------------------------------------------------------
# `dt`         `::Float` $\Delta t$, current time step.
#
# `A`          `::Float` cross sectional area along the cell.
#
# `u`          `::Vessel` longitudinal velocity along the cell.
#
# `b`          `::Blood` blood data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# The numerical source is updated by means of a first order implicit backward
# Euler method.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ---------- -----------------------------------------------------------------
# `S`        `::Float` source term.
# ----------------------------------------------------------------------------
#
function source(dt :: Float64, A :: Float64, u :: Float64, b :: Blood)

  return -b.Cf*u*dt

end

# *function* __`calculateDeltaT`__ $\rightarrow$ `DT::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `vessels`        `::Array{Vessel, 1}`
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# This function calculates the smallest $\Delta t$ allowed by all vessels
# geometries and properties.
# ----------------------------------------------------------------------------
#
function calculateDeltaT(vessels, dt :: Array{Float64, 1})

  # For each vessel the $\Delta t$ is computed as
  # $$
  #   \Delta t = \frac{\Delta x}{S_{max}} C_{cfl}
  # $$
  # where $S_{max}$ is the maximum between the forward and the backward
  # characteristics, and $C_{cfl}$ is the Courant-Friedrichs-Lewy condition
  # defined by the user.
  i = 1
  for v in vessels

    lambdap = zeros(Float64, v.M)
    for j in 1:v.M
      lambdap[j] = v.u[j] + v.c[j]
    end

    Smax = maximum(abs.(lambdap))

    dt[i] = v.dx*v.Ccfl/Smax

    i+=1
  end
  #
  # --------------------------------------------------------------------------
  # Returns:
  # ----------- --------------------------------------------------------------
  # `DT`        `::Float` $\Delta t$ computed as the smallest between
  #             the smallest $\Delta t$ associated with each vessel and
  #             the $\Delta t$ limiter $0.001$.
  # --------------------------------------------------------------------------
  return minimum(dt)

end

# *function* __`solveModel`__
#
# ----------------------------------------------------------------------------
# Parameters:
# --------------- ------------------------------------------------------------
# `grafo`         `::GenericGraph` graph data structure containing all the
#                 vessels.
#
# `edgess`        `::Array` collection of `edge` objects.
#
# `vessels`       `::Array{Vessel, 1}` collection of vessels in the system.
#
# `heart`         `::Heart` heart data structure.
#
# `blood`         `::Blood` blood data structure.
#
# `dt`            `::Float` $\Delta t$, current time step.
#
# `current_time`  `::Float` current simulation time.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# See the [detailed example](grafo.html).
#
# ----------------------------------------------------------------------------
#
function solveModel(vessels, heart :: Heart, edge_list,
                    blood :: Blood, dt :: Float64, current_time :: Float64)

  for j in 1:size(edge_list)[1]
    i = edge_list[j,1]
    s = edge_list[j,2]
    t = edge_list[j,3]
    v = vessels[i]
    lbl = v.label

    if size(find(edge_list[:,3] .== s))[1] == 0
      openBF.setInletBC(current_time, dt, v, heart)
    end

    openBF.MUSCL(v, dt, blood)

    if size(find(edge_list[:,2] .== t))[1] == 0
      openBF.setOutletBC(dt, v)
      # println("\t Outlet vessel - Compute outlet BC")

    elseif size(find(edge_list[:,2] .== t))[1] == 2
      d1_i = find(edge_list[:,2] .== t)[1]
      d2_i = find(edge_list[:,2] .== t)[2]

      openBF.joinVessels(blood, v, vessels[d1_i], vessels[d2_i])
      # println("\t\t Bifurcation at node $t between vessels $i, $d1_i, and $d2_i")

    elseif size(find(edge_list[:,3] .== t))[1] == 1
      d_i = find(edge_list[:,2] .== t)[1]

      openBF.joinVessels(blood, v, vessels[d_i])
      # println("\t\t Junction at node $t between vessels $i and $d_i")

    elseif size(find(edge_list[:,3] .== t))[1] == 2
      p1_i = find(edge_list[:,3] .== t)[1]
      p2_i = find(edge_list[:,3] .== t)[2]
      if maximum([p1_i, p2_i]) == i
        p2_i = minimum([p1_i, p2_i])
        d = find(edge_list[:,2] .== t)[1]
        openBF.solveAnastomosis(v, vessels[p2_i], vessels[d])
        # println("\t\t Anastomosis at node $s between vessels $i, $p2_i, and $d")
      end
    end
  end
end

# ### References
#
# [^1]: Toro, Eleuterio F. Riemann solvers and numerical methods for fluid
# dynamics: a practical introduction. Springer Science & Business Media, 2009.


text:
name = "nets"
for (_, _, files) in walkdir(joinpath(@__DIR__, name))
    for file in files
        endswith(file, ".jl") && include(joinpath(name, file))
    end
end


text: Dict{String,Any}("foo" => Any[Dict{String,Any}("bar" => Dict{String,Any}("value" => "\"{{baz}}\"","type" => "string"))])

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TestLU

using Test, LinearAlgebra, Random
using LinearAlgebra: ldiv!, BlasInt, BlasFloat, rdiv!

n = 10

# Split n into 2 parts for tests needing two matrices
n1 = div(n, 2)
n2 = 2*n1

Random.seed!(1234321)

areal = randn(n,n)/2
aimg  = randn(n,n)/2
breal = randn(n,2)/2
bimg  = randn(n,2)/2
creal = randn(n)/2
cimg  = randn(n)/2
dureal = randn(n-1)/2
duimg  = randn(n-1)/2
dlreal = randn(n-1)/2
dlimg  = randn(n-1)/2
dreal = randn(n)/2
dimg  = randn(n)/2

@testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
    a = eltya == Int ? rand(1:7, n, n) :
        convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
    d = if eltya == Int
        Tridiagonal(rand(1:7, n-1), rand(1:7, n), rand(1:7, n-1))
    elseif eltya <: Complex
        convert(Tridiagonal{eltya}, Tridiagonal(
            complex.(dlreal, dlimg), complex.(dreal, dimg), complex.(dureal, duimg)))
    else
        convert(Tridiagonal{eltya}, Tridiagonal(dlreal, dreal, dureal))
    end
    ε = εa = eps(abs(float(one(eltya))))

    if eltya <: BlasFloat
        @testset "LU factorization for Number" begin
            num = rand(eltya)
            @test (lu(num)...,) == (hcat(one(eltya)), hcat(num), [1])
            @test convert(Array, lu(num)) ≈ eltya[num]
        end
        @testset "Balancing in eigenvector calculations" begin
            A = convert(Matrix{eltya}, [ 3.0     -2.0      -0.9     2*eps(real(one(eltya)));
                                       -2.0      4.0       1.0    -eps(real(one(eltya)));
                                       -eps(real(one(eltya)))/4  eps(real(one(eltya)))/2  -1.0     0;
                                       -0.5     -0.5       0.1     1.0])
            F = eigen(A, permute=false, scale=false)
            @test F.vectors*Diagonal(F.values)/F.vectors ≈ A
            F = eigen(A)
            # @test norm(F.vectors*Diagonal(F.values)/F.vectors - A) > 0.01
        end
    end
    κ  = cond(a,1)
    @testset "(Automatic) Square LU decomposition" begin
        lua   = factorize(a)
        @test_throws ErrorException lua.Z
        l,u,p = lua.L, lua.U, lua.p
        ll,ul,pl = lu(a)
        @test ll * ul ≈ a[pl,:]
        @test l*u ≈ a[p,:]
        @test (l*u)[invperm(p),:] ≈ a
        @test a * inv(lua) ≈ Matrix(I, n, n)
        @test copy(lua) == lua
        if eltya <: BlasFloat
            # test conversion of LU factorization's numerical type
            bft = eltya <: Real ? LinearAlgebra.LU{BigFloat} : LinearAlgebra.LU{Complex{BigFloat}}
            bflua = convert(bft, lua)
            @test bflua.L*bflua.U ≈ big.(a)[p,:] rtol=ε
        end
        # compact printing
        lstring = sprint(show,l)
        ustring = sprint(show,u)
    end
    κd    = cond(Array(d),1)
    @testset "Tridiagonal LU" begin
        lud   = lu(d)
        @test LinearAlgebra.issuccess(lud)
        @test lu(lud) == lud
        @test_throws ErrorException lud.Z
        @test lud.L*lud.U ≈ lud.P*Array(d)
        @test lud.L*lud.U ≈ Array(d)[lud.p,:]
        @test AbstractArray(lud) ≈ d
    end
    @testset for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int)
        b  = eltyb == Int ? rand(1:5, n, 2) :
            convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal)
        c  = eltyb == Int ? rand(1:5, n) :
            convert(Vector{eltyb}, eltyb <: Complex ? complex.(creal, cimg) : creal)
        εb = eps(abs(float(one(eltyb))))
        ε  = max(εa,εb)
        @testset "(Automatic) Square LU decomposition" begin
            lua   = factorize(a)
            let Bs = copy(b), Cs = copy(c)
                for (bb, cc) in ((Bs, Cs), (view(Bs, 1:n, 1), view(Cs, 1:n)))
                    @test norm(a*(lua\bb) - bb, 1) < ε*κ*n*2 # Two because the right hand side has two columns
                    @test norm(a'*(lua'\bb) - bb, 1) < ε*κ*n*2 # Two because the right hand side has two columns
                    @test norm(a'*(lua'\a') - a', 1) < ε*κ*n^2
                    @test norm(a*(lua\cc) - cc, 1) < ε*κ*n # cc is a vector
                    @test norm(a'*(lua'\cc) - cc, 1) < ε*κ*n # cc is a vector
                    @test AbstractArray(lua) ≈ a
                    @test norm(transpose(a)*(transpose(lua)\bb) - bb,1) < ε*κ*n*2 # Two because the right hand side has two columns
                    @test norm(transpose(a)*(transpose(lua)\cc) - cc,1) < ε*κ*n
                end

                # Test whether Ax_ldiv_B!(y, LU, x) indeed overwrites y
                resultT = typeof(oneunit(eltyb) / oneunit(eltya))

                b_dest = similar(b, resultT)
                c_dest = similar(c, resultT)

                ldiv!(b_dest, lua, b)
                ldiv!(c_dest, lua, c)
                @test norm(b_dest - lua \ b, 1) < ε*κ*2n
                @test norm(c_dest - lua \ c, 1) < ε*κ*n

                ldiv!(b_dest, transpose(lua), b)
                ldiv!(c_dest, transpose(lua), c)
                @test norm(b_dest - transpose(lua) \ b, 1) < ε*κ*2n
                @test norm(c_dest - transpose(lua) \ c, 1) < ε*κ*n

                ldiv!(b_dest, adjoint(lua), b)
                ldiv!(c_dest, adjoint(lua), c)
                @test norm(b_dest - lua' \ b, 1) < ε*κ*2n
                @test norm(c_dest - lua' \ c, 1) < ε*κ*n

                if eltyb != Int && !(eltya <: Complex) || eltya <: Complex && eltyb <: Complex
                    p = Matrix(b')
                    q = Matrix(c')
                    p_dest = copy(p)
                    q_dest = copy(q)
                    rdiv!(p_dest, lua)
                    rdiv!(q_dest, lua)
                    @test norm(p_dest - p / lua, 1) < ε*κ*2n
                    @test norm(q_dest - q / lua, 1) < ε*κ*n
                end
            end
            if eltya <: BlasFloat && eltyb <: BlasFloat
                e = rand(eltyb,n,n)
                @test norm(e/lua - e/a,1) < ε*κ*n^2
            end
        end
        @testset "Tridiagonal LU" begin
            lud   = factorize(d)
            f = zeros(eltyb, n+1)
            @test_throws DimensionMismatch lud\f
            @test_throws DimensionMismatch transpose(lud)\f
            @test_throws DimensionMismatch lud'\f
            @test_throws DimensionMismatch LinearAlgebra.ldiv!(transpose(lud), f)
            let Bs = copy(b)
                for bb in (Bs, view(Bs, 1:n, 1))
                    @test norm(d*(lud\bb) - bb, 1) < ε*κd*n*2 # Two because the right hand side has two columns
                    if eltya <: Real
                        @test norm((transpose(lud)\bb) - Array(transpose(d))\bb, 1) < ε*κd*n*2 # Two because the right hand side has two columns
                        if eltya != Int && eltyb != Int
                            @test norm(LinearAlgebra.ldiv!(transpose(lud), copy(bb)) - Array(transpose(d))\bb, 1) < ε*κd*n*2
                        end
                    end
                    if eltya <: Complex
                        @test norm((lud'\bb) - Array(d')\bb, 1) < ε*κd*n*2 # Two because the right hand side has two columns
                    end
                end
            end
            if eltya <: BlasFloat && eltyb <: BlasFloat
                e = rand(eltyb,n,n)
                @test norm(e/lud - e/d,1) < ε*κ*n^2
                @test norm((transpose(lud)\e') - Array(transpose(d))\e',1) < ε*κd*n^2
                #test singular
                du = rand(eltya,n-1)
                dl = rand(eltya,n-1)
                dd = rand(eltya,n)
                dd[1] = zero(eltya)
                du[1] = zero(eltya)
                dl[1] = zero(eltya)
                zT = Tridiagonal(dl,dd,du)
                @test !LinearAlgebra.issuccess(lu(zT; check = false))
            end
        end
        @testset "Thin LU" begin
            lua   = @inferred lu(a[:,1:n1])
            @test lua.L*lua.U ≈ lua.P*a[:,1:n1]
        end
        @testset "Fat LU" begin
            lua   = lu(a[1:n1,:])
            @test lua.L*lua.U ≈ lua.P*a[1:n1,:]
        end
    end

    @testset "LU of Symmetric/Hermitian" begin
        for HS in (Hermitian(a'a), Symmetric(a'a))
            luhs = lu(HS)
            @test luhs.L*luhs.U ≈ luhs.P*Matrix(HS)
        end
    end
end

@testset "Singular matrices" for T in (Float64, ComplexF64)
    A = T[1 2; 0 0]
    @test_throws SingularException lu(A)
    @test_throws SingularException lu!(copy(A))
    @test_throws SingularException lu(A; check = true)
    @test_throws SingularException lu!(copy(A); check = true)
    @test !issuccess(lu(A; check = false))
    @test !issuccess(lu!(copy(A); check = false))
    @test_throws ZeroPivotException lu(A, Val(false))
    @test_throws ZeroPivotException lu!(copy(A), Val(false))
    @test_throws ZeroPivotException lu(A, Val(false); check = true)
    @test_throws ZeroPivotException lu!(copy(A), Val(false); check = true)
    @test !issuccess(lu(A, Val(false); check = false))
    @test !issuccess(lu!(copy(A), Val(false); check = false))
    F = lu(A; check = false)
    @test sprint((io, x) -> show(io, "text/plain", x), F) ==
        "Failed factorization of type $(typeof(F))"
end

@testset "conversion" begin
    Random.seed!(3)
    a = Tridiagonal(rand(9),rand(10),rand(9))
    fa = Array(a)
    falu = lu(fa)
    alu = lu(a)
    falu = convert(typeof(falu),alu)
    @test AbstractArray(alu) == fa
end

@testset "Rational Matrices" begin
    ## Integrate in general tests when more linear algebra is implemented in julia
    a = convert(Matrix{Rational{BigInt}}, rand(1:10//1,n,n))/n
    b = rand(1:10,n,2)
    @inferred lu(a)
    lua   = factorize(a)
    l,u,p = lua.L, lua.U, lua.p
    @test l*u ≈ a[p,:]
    @test l[invperm(p),:]*u ≈ a
    @test a*inv(lua) ≈ Matrix(I, n, n)
    let Bs = b
        for b in (Bs, view(Bs, 1:n, 1))
            @test a*(lua\b) ≈ b
        end
    end
    @test @inferred(det(a)) ≈ det(Array{Float64}(a))
end

@testset "Rational{BigInt} and BigFloat Hilbert Matrix" begin
    ## Hilbert Matrix (very ill conditioned)
    ## Testing Rational{BigInt} and BigFloat version
    nHilbert = 50
    H = Rational{BigInt}[1//(i+j-1) for i = 1:nHilbert,j = 1:nHilbert]
    Hinv = Rational{BigInt}[(-1)^(i+j)*(i+j-1)*binomial(nHilbert+i-1,nHilbert-j)*
        binomial(nHilbert+j-1,nHilbert-i)*binomial(i+j-2,i-1)^2
        for i = big(1):nHilbert,j=big(1):nHilbert]
    @test inv(H) == Hinv
    setprecision(2^10) do
        @test norm(Array{Float64}(inv(float(H)) - float(Hinv))) < 1e-100
    end
end

@testset "logdet" begin
    @test @inferred(logdet(ComplexF32[1.0f0 0.5f0; 0.5f0 -1.0f0])) === 0.22314355f0 + 3.1415927f0im
    @test_throws DomainError logdet([1 1; 1 -1])
end

@testset "REPL printing" begin
        bf = IOBuffer()
        show(bf, "text/plain", lu(Matrix(I, 4, 4)))
        seekstart(bf)
        @test String(take!(bf)) == """
LinearAlgebra.LU{Float64,Array{Float64,2}}
L factor:
4×4 Array{Float64,2}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0
U factor:
4×4 Array{Float64,2}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0"""
end

@testset "propertynames" begin
    names = sort!(collect(string.(Base.propertynames(lu(rand(3,3))))))
    @test names == ["L", "P", "U", "p"]
    allnames = sort!(collect(string.(Base.propertynames(lu(rand(3,3)), true))))
    @test allnames == ["L", "P", "U", "factors", "info", "ipiv", "p"]
end

include("trickyarithmetic.jl")

@testset "lu with type whose sum is another type" begin
    A = TrickyArithmetic.A[1 2; 3 4]
    ElT = TrickyArithmetic.D{TrickyArithmetic.C,TrickyArithmetic.C}
    B = lu(A, Val(false))
    @test B isa LinearAlgebra.LU{ElT,Matrix{ElT}}
end

@testset "Issue #30917. Determinant of integer matrix" begin
    @test det([1 1 0 0 1 0 0 0
               1 0 1 0 0 1 0 0
               1 0 0 1 0 0 1 0
               0 1 1 1 0 0 0 0
               0 1 0 0 0 0 1 1
               0 0 1 0 1 0 0 1
               0 0 0 1 1 1 0 0
               0 0 0 0 1 1 0 1]) ≈ 6
end

@testset "Issue #33177. No ldiv!(LU, Adjoint)" begin
    A = [1 0; 1 1]
    B = [1 2; 2 8]
    F = lu(B)
    @test (A  / F') * B == A
    @test (A' / F') * B == A'

    a = complex.(randn(2), randn(2))
    @test (a' / F') * B ≈ a'
    @test (transpose(a) / F') * B ≈ transpose(a)

    A = complex.(randn(2, 2), randn(2, 2))
    @test (A' / F') * B ≈ A'
    @test (transpose(A) / F') * B ≈ transpose(A)
end

@testset "0x0 matrix" begin
    A = ones(0, 0)
    F = lu(A)
    @test F.U == ones(0, 0)
    @test F.L == ones(0, 0)
    @test F.P == ones(0, 0)
    @test F.p == []
end

end # module TestLU

max_stars_repo_name: Dictino/julia

max_issues_repo_name: Dictino/julia

max_forks_repo_name: Dictino/julia


text:
function test_debug_envs()
    env = AutoEnvs.DeterministicSingleStepDebugEnv()
    x = reset(env)
    nx, r, t, _ = step(env, 1)
end

@time test_debug_envs()


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct HydroDispatch <: HydroGen
        name::String
        available::Bool
        bus::Bus
        active_power::Float64
        reactive_power::Float64
        rating::Float64
        prime_mover::PrimeMovers
        active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        reactive_power_limits::Union{Nothing, NamedTuple{(:min, :max), Tuple{Float64, Float64}}}
        ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
        time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
        base_power::Float64
        operation_cost::OperationalCost
        services::Vector{Service}
        dynamic_injector::Union{Nothing, DynamicInjection}
        ext::Dict{String, Any}
        time_series_container::InfrastructureSystems.TimeSeriesContainer
        internal::InfrastructureSystemsInternal
    end

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `active_power::Float64`
- `reactive_power::Float64`, validation range: `reactive_power_limits`, action if invalid: `warn`
- `rating::Float64`: Thermal limited MVA Power Output of the unit. <= Capacity, validation range: `(0, nothing)`, action if invalid: `error`
- `prime_mover::PrimeMovers`: Prime mover technology according to EIA 923
- `active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`, validation range: `(0, nothing)`, action if invalid: `warn`
- `reactive_power_limits::Union{Nothing, NamedTuple{(:min, :max), Tuple{Float64, Float64}}}`, action if invalid: `warn`
- `ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}`: ramp up and ramp down limits in MW (in component base per unit) per minute, validation range: `(0, nothing)`, action if invalid: `error`
- `time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}`: Minimum up and Minimum down time limits in hours, validation range: `(0, nothing)`, action if invalid: `error`
- `base_power::Float64`: Base power of the unit in MVA, validation range: `(0, nothing)`, action if invalid: `warn`
- `operation_cost::OperationalCost`: Operation Cost of Generation [`OperationalCost`](@ref)
- `services::Vector{Service}`: Services that this device contributes to
- `dynamic_injector::Union{Nothing, DynamicInjection}`: corresponding dynamic injection device
- `ext::Dict{String, Any}`
- `time_series_container::InfrastructureSystems.TimeSeriesContainer`: internal time_series storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct HydroDispatch <: HydroGen
    name::String
    available::Bool
    bus::Bus
    active_power::Float64
    reactive_power::Float64
    "Thermal limited MVA Power Output of the unit. <= Capacity"
    rating::Float64
    "Prime mover technology according to EIA 923"
    prime_mover::PrimeMovers
    active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    reactive_power_limits::Union{Nothing, NamedTuple{(:min, :max), Tuple{Float64, Float64}}}
    "ramp up and ramp down limits in MW (in component base per unit) per minute"
    ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    "Minimum up and Minimum down time limits in hours"
    time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    "Base power of the unit in MVA"
    base_power::Float64
    "Operation Cost of Generation [`OperationalCost`](@ref)"
    operation_cost::OperationalCost
    "Services that this device contributes to"
    services::Vector{Service}
    "corresponding dynamic injection device"
    dynamic_injector::Union{Nothing, DynamicInjection}
    ext::Dict{String, Any}
    "internal time_series storage"
    time_series_container::InfrastructureSystems.TimeSeriesContainer
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function HydroDispatch(name, available, bus, active_power, reactive_power, rating, prime_mover, active_power_limits, reactive_power_limits, ramp_limits, time_limits, base_power, operation_cost=TwoPartCost(0.0, 0.0), services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), )
    HydroDispatch(name, available, bus, active_power, reactive_power, rating, prime_mover, active_power_limits, reactive_power_limits, ramp_limits, time_limits, base_power, operation_cost, services, dynamic_injector, ext, time_series_container, InfrastructureSystemsInternal(), )
end

function HydroDispatch(; name, available, bus, active_power, reactive_power, rating, prime_mover, active_power_limits, reactive_power_limits, ramp_limits, time_limits, base_power, operation_cost=TwoPartCost(0.0, 0.0), services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), internal=InfrastructureSystemsInternal(), )
    HydroDispatch(name, available, bus, active_power, reactive_power, rating, prime_mover, active_power_limits, reactive_power_limits, ramp_limits, time_limits, base_power, operation_cost, services, dynamic_injector, ext, time_series_container, internal, )
end

# Constructor for demo purposes; non-functional.
function HydroDispatch(::Nothing)
    HydroDispatch(;
        name="init",
        available=false,
        bus=Bus(nothing),
        active_power=0.0,
        reactive_power=0.0,
        rating=0.0,
        prime_mover=PrimeMovers.HY,
        active_power_limits=(min=0.0, max=0.0),
        reactive_power_limits=nothing,
        ramp_limits=nothing,
        time_limits=nothing,
        base_power=0.0,
        operation_cost=TwoPartCost(nothing),
        services=Device[],
        dynamic_injector=nothing,
        ext=Dict{String, Any}(),
        time_series_container=InfrastructureSystems.TimeSeriesContainer(),
    )
end

"""Get [`HydroDispatch`](@ref) `name`."""
get_name(value::HydroDispatch) = value.name
"""Get [`HydroDispatch`](@ref) `available`."""
get_available(value::HydroDispatch) = value.available
"""Get [`HydroDispatch`](@ref) `bus`."""
get_bus(value::HydroDispatch) = value.bus
"""Get [`HydroDispatch`](@ref) `active_power`."""
get_active_power(value::HydroDispatch) = get_value(value, value.active_power)
"""Get [`HydroDispatch`](@ref) `reactive_power`."""
get_reactive_power(value::HydroDispatch) = get_value(value, value.reactive_power)
"""Get [`HydroDispatch`](@ref) `rating`."""
get_rating(value::HydroDispatch) = get_value(value, value.rating)
"""Get [`HydroDispatch`](@ref) `prime_mover`."""
get_prime_mover(value::HydroDispatch) = value.prime_mover
"""Get [`HydroDispatch`](@ref) `active_power_limits`."""
get_active_power_limits(value::HydroDispatch) = get_value(value, value.active_power_limits)
"""Get [`HydroDispatch`](@ref) `reactive_power_limits`."""
get_reactive_power_limits(value::HydroDispatch) = get_value(value, value.reactive_power_limits)
"""Get [`HydroDispatch`](@ref) `ramp_limits`."""
get_ramp_limits(value::HydroDispatch) = get_value(value, value.ramp_limits)
"""Get [`HydroDispatch`](@ref) `time_limits`."""
get_time_limits(value::HydroDispatch) = value.time_limits
"""Get [`HydroDispatch`](@ref) `base_power`."""
get_base_power(value::HydroDispatch) = value.base_power
"""Get [`HydroDispatch`](@ref) `operation_cost`."""
get_operation_cost(value::HydroDispatch) = value.operation_cost
"""Get [`HydroDispatch`](@ref) `services`."""
get_services(value::HydroDispatch) = value.services
"""Get [`HydroDispatch`](@ref) `dynamic_injector`."""
get_dynamic_injector(value::HydroDispatch) = value.dynamic_injector
"""Get [`HydroDispatch`](@ref) `ext`."""
get_ext(value::HydroDispatch) = value.ext
"""Get [`HydroDispatch`](@ref) `time_series_container`."""
get_time_series_container(value::HydroDispatch) = value.time_series_container
"""Get [`HydroDispatch`](@ref) `internal`."""
get_internal(value::HydroDispatch) = value.internal

"""Set [`HydroDispatch`](@ref) `name`."""
set_name!(value::HydroDispatch, val) = value.name = val
"""Set [`HydroDispatch`](@ref) `available`."""
set_available!(value::HydroDispatch, val) = value.available = val
"""Set [`HydroDispatch`](@ref) `bus`."""
set_bus!(value::HydroDispatch, val) = value.bus = val
"""Set [`HydroDispatch`](@ref) `active_power`."""
set_active_power!(value::HydroDispatch, val) = value.active_power = set_value(value, val)
"""Set [`HydroDispatch`](@ref) `reactive_power`."""
set_reactive_power!(value::HydroDispatch, val) = value.reactive_power = set_value(value, val)
"""Set [`HydroDispatch`](@ref) `rating`."""
set_rating!(value::HydroDispatch, val) = value.rating = set_value(value, val)
"""Set [`HydroDispatch`](@ref) `prime_mover`."""
set_prime_mover!(value::HydroDispatch, val) = value.prime_mover = val
"""Set [`HydroDispatch`](@ref) `active_power_limits`."""
set_active_power_limits!(value::HydroDispatch, val) = value.active_power_limits = set_value(value, val)
"""Set [`HydroDispatch`](@ref) `reactive_power_limits`."""
set_reactive_power_limits!(value::HydroDispatch, val) = value.reactive_power_limits = set_value(value, val)
"""Set [`HydroDispatch`](@ref) `ramp_limits`."""
set_ramp_limits!(value::HydroDispatch, val) = value.ramp_limits = set_value(value, val)
"""Set [`HydroDispatch`](@ref) `time_limits`."""
set_time_limits!(value::HydroDispatch, val) = value.time_limits = val
"""Set [`HydroDispatch`](@ref) `base_power`."""
set_base_power!(value::HydroDispatch, val) = value.base_power = val
"""Set [`HydroDispatch`](@ref) `operation_cost`."""
set_operation_cost!(value::HydroDispatch, val) = value.operation_cost = val
"""Set [`HydroDispatch`](@ref) `services`."""
set_services!(value::HydroDispatch, val) = value.services = val
"""Set [`HydroDispatch`](@ref) `ext`."""
set_ext!(value::HydroDispatch, val) = value.ext = val
"""Set [`HydroDispatch`](@ref) `time_series_container`."""
set_time_series_container!(value::HydroDispatch, val) = value.time_series_container = val


text:
using Judo
Judo.collate("Compose", template=Pkg.dir("Compose", "doc", "template"))


text:
module Morbit

using DocStringExtensions

using Printf: @sprintf

# steepest descent and directed search
using LinearAlgebra: norm, pinv
import JuMP;
import OSQP;

# LagrangeModels and PS
import NLopt;

using Parameters: @with_kw, @unpack, @pack!
using MathOptInterface;
const MOI = MathOptInterface;

import ThreadSafeDicts: ThreadSafeDict;
using Memoize;

import UUIDs;

import FiniteDiff#erences
const FD = FiniteDiff#erences

import ForwardDiff
const AD = ForwardDiff

import Logging: LogLevel, @logmsg
import Logging

export MixedMOP, RbfConfig, ExactConfig, TaylorConfig, TaylorApproximateConfig, TaylorCallbackConfig
export add_objective!, add_vector_objective!
export optimize
export SteepestDescentConfig, PascolettiSerafiniConfig

include("custom_logging.jl")

include("shorthands.jl");
include("Interfaces.jl");
include("diff_wrappers.jl");

# implementations (order should not matter)
include("VectorObjectiveFunction.jl");
include("MixedMOP.jl")
include("StaticMOP.jl")

include("ResultImplementation.jl")
include("DataBaseImplementation.jl")
include("IterDataImplementation.jl")
include("SurrogatesImplementation.jl");

include("ConfigImplementations.jl")

# utilities
include("adding_objectives.jl");
include("descent.jl")
#include("saving.jl")
include("utilities.jl")

# this method is used in the algorithm and delegates the work
function shrink_radius( ac :: AbstractConfig, Δ, steplength )
	return shrink_radius( Val( radius_update_method(ac) ), ac, Δ, steplength )
end

"Shrink radius according to `γ * Δ`."
function shrink_radius( ::Val{:standard}, ac, Δ, steplength )
	return Δ * _gamma_shrink( ac )
end

"Shrink radius according to `γ * ||s||`."
function shrink_radius( ::Val{:steplength}, ac, Δ, steplength )
	return steplength * _gamma_shrink( ac )
end

# this method is used in the algorithm and delegates the work
function shrink_radius_much( ac :: AbstractConfig, Δ, steplength )
	return shrink_radius_much( Val( radius_update_method(ac) ), ac, Δ, steplength )
end

"Shrink radius much according to `γ * Δ`."
function shrink_radius_much( ::Val{:standard}, ac, Δ, steplength )
	return Δ * _gamma_shrink_much( ac )
end

"Shrink radius according to `γ * ||s||`."
function shrink_radius_much( ::Val{:steplength}, ac, Δ, steplength )
	return steplength * _gamma_shrink_much( ac )
end

# this method is used in the algorithm and delegates the work
function grow_radius( ac :: AbstractConfig, Δ, steplength )
	return grow_radius( Val( radius_update_method(ac) ), ac, Δ, steplength )
end

"Grow radius according to `min( Δ_max, γ * Δ )`."
function grow_radius( ::Val{:standard}, ac, Δ, steplength )
	return min( get_delta_max(ac), _gamma_grow(ac) * Δ )
end

"Grow radius according to `min( Δ_max, (γ + ||s||/Δ) * Δ )`"
function grow_radius( ::Val{:steplength}, ac, Δ, steplength )
	return min( get_delta_max(ac), ( _gamma_grow(ac) .+ steplength ./ Δ ) .* Δ )
end

# we expect mop :: MixedMOP, but should work for static MOP if everything
# is set up properly
function initialize_data( mop :: AbstractMOP, x0 :: Vec, fx0 :: Vec = MIN_PRECISION[];
    algo_config :: Union{AbstractConfig, Nothing} = nothing,
    populated_db :: Union{AbstractDB, Nothing} = nothing )

	if num_objectives(mop) == 0
		error("`mop` has no objectives!")
	end

	@logmsg loglevel1 """\n
    |--------------------------------------------
    | Initialization
    |--------------------------------------------
	"""

	@logmsg loglevel1 "The evaluation counter of `mop` is reset."
	reset_evals!( mop )
	# initialize first iteration site
	@assert !isempty( x0 ) "Please provide a non-empty feasible starting point `x0`."

	# for backwards-compatibility with unconstrained problems:
	if num_vars(mop) == 0
		# will error if mop isa StaticMOP
		MOI.add_variables(mop, length(x0))
	end

	@assert length(x0) == num_vars( mop ) "Number of variables in `mop` does not match length of `x0`."

	# make problem static
	smop = StaticMOP(mop)

	# scale bound vars to unit hypercube
	x_scaled = scale( x0, smop )
	# ensure at least single-precision
	XTe = Base.promote_eltype( x_scaled, MIN_PRECISION )
	x = XTe.(x_scaled)
	XT = typeof(x)

	# initalize first objective vector
	if isempty( fx0 )
		# if no starting function value was provided, eval objectives
		fx_sorted = eval_all_objectives( smop, x )
	else
		fx_sorted = apply_internal_sorting( fx0, smop )
	end

	YTe = Base.promote_eltype( fx_sorted, MIN_PRECISION )
	fx = YTe.(fx_sorted)
	YT = typeof( fx )

	if isnothing( algo_config )
		ac = DefaultConfig()
	else
		ac = algo_config
	end

	# initialize iter data obj
	Δ_0 = let Δ = get_delta_0( ac );
		T = promote_type( typeof(Δ), XTe )
		T.(Δ)
	end
	DT = typeof(Δ_0)

	id = init_iter_data( IterData, x, fx, Δ_0 )

	# initialize database
	if !isnothing(populated_db)
		# has a database been provided? if yes, prepare (scale vars, sort values)
		data_base = populated_db
		transform!( data_base, mop )
	else
		result_type = Result{ XT, YT }
		data_base = init_db( ArrayDB, result_type, Nothing )
		set_transformed!(data_base, true)
	end

	# make sure, x & fx are in database
	x_id = ensure_contains_values!(data_base, x, fx)
	set_x_index!(id, x_id)

	# make the problem static
	# initialize surrogate models
	sc = init_surrogates(SurrogateContainer, smop, id, data_base, ac )

	# now that we have meta data available, we retrive the
	# right saveable type and make a new database, that can handle `CT`
	CT = get_saveable_type( sc )
	saveable_type = IterSaveable{DT, CT}
	new_data_base = copy_db(data_base; saveable_type )
	return (smop, id, new_data_base, sc, ac)
end

function iterate!( iter_data :: AbstractIterData, data_base :: AbstractDB, mop :: AbstractMOP,
	sc :: SurrogateContainer, algo_config :: AbstractConfig ) :: STOP_CODE

	# read iteration data
	x = get_x( iter_data )
	fx = get_fx( iter_data )
	Δ = get_delta( iter_data )

	# check (some) stopping conditions
	# (rest is done at end of this function, when trial point is known)
	if get_num_iterations(iter_data) >= max_iter(algo_config)
        @logmsg loglevel1 "Stopping. Maximum number of iterations reached."
        return MAX_ITER
    end

    if !_budget_okay(mop, algo_config)
        @logmsg loglevel1 "Stopping. Computational budget is exhausted."
        return BUDGET_EXHAUSTED
    end

    if Δ_abs_test( Δ, algo_config )
        return TOLERANCE
    end

	# read algorithm parameters from config
	ν_success = _nu_success( algo_config )
	ν_accept = _nu_accept( algo_config )
	μ = _mu( algo_config )
	β = _beta( algo_config )
	ε_c = _eps_crit( algo_config )
	γ_c = _gamma_crit( algo_config )

	SAVEABLE_TYPE = get_saveable_type( data_base )

    # set iteration counter
    if it_stat(iter_data) != MODELIMPROVING || count_nonlinear_iterations( algo_config )
        inc_num_iterations!(iter_data)
        set_num_model_improvements!(iter_data, 0);
    else
        inc_num_model_improvements!(iter_data)
    end

    @logmsg loglevel1 """\n
        |--------------------------------------------
        |Iteration $(get_num_iterations(iter_data)).$(get_num_model_improvements(iter_data))
        |--------------------------------------------
        |  Current trust region radius is $(Δ).
        |  Current number of function evals is $(num_evals(mop)).
        |  Iterate is $(_prettify(unscale(x, mop)))
        |  Values are $(_prettify(reverse_internal_sorting(fx,mop)))
        |--------------------------------------------
    """

    # update surrogate models
    if get_num_iterations(iter_data) > 1
        if it_stat(iter_data) == MODELIMPROVING
            improve_surrogates!( sc, mop, iter_data, data_base, algo_config; ensure_fully_linear = false );
        else
            update_surrogates!( sc, mop, iter_data, data_base, algo_config; ensure_fully_linear = false );
        end
    end

	# TODO: if problem is constrained by functions:
	# do "normal step here"

    # calculate descent step and criticality value
    ω, ω_data = get_criticality(mop, iter_data, data_base, sc, algo_config)
    @logmsg loglevel1 "Criticality is ω = $(ω)."

    # stop at the end of the this loop run?
    if ω_Δ_rel_test(ω, Δ, algo_config) || ω_abs_test( ω, algo_config )
        return CRITICAL
    end

    # Criticallity test
	_fully_linear_sc = fully_linear(sc)
	if ω <= ε_c && (!_fully_linear_sc|| all(Δ .> μ * ω))
		@logmsg loglevel1 "Entered Criticallity Test."
		DO_CRIT_LOOPS = true
        if !_fully_linear_sc
            @logmsg loglevel1 "Ensuring that all models are fully linear."
            update_surrogates!( sc, mop, iter_data, data_base, algo_config; ensure_fully_linear = true );

            ω, ω_data = get_criticality(mop, iter_data, data_base, sc, algo_config )
            if !fully_linear(sc)
				DO_CRIT_LOOPS = false
				@logmsg loglevel2 "Could not make all models fully linear. Trying one last descent step."
			end
		end
		if DO_CRIT_LOOPS
			num_critical_loops = 0
			while all(Δ .> μ * ω )
				@logmsg loglevel2 "Criticality loop $(num_critical_loops + 1)."
				if num_critical_loops >= max_critical_loops(algo_config)
					@logmsg loglevel1 "Maximum number ($(max_critical_loops(algo_config))) of critical loops reached. Exiting..."
					return CRITICAL
				end
				if !_budget_okay(mop, algo_config)
					@logmsg loglevel1 "Computational budget exhausted. Exiting…"
					return BUDGET_EXHAUSTED
				end
				# shrink radius
				Δ = γ_c .* Δ
				set_delta!(iter_data, Δ)

				# make models fully linear on smaller trust region
				update_surrogates!( sc, mop, iter_data, data_base, algo_config; ensure_fully_linear = true )

				# (re)calculate criticality
				ω, ω_data = get_criticality( mop, iter_data, data_base, sc, algo_config )
				num_critical_loops += 1

				if Δ_abs_test( Δ, algo_config ) ||
					ω_Δ_rel_test(ω, Δ, algo_config) || ω_abs_test( ω, algo_config )
					return CRITICAL
				end
				if !fully_linear(sc)
					@logmsg loglevel2 "Could not make all models fully linear. Trying one last descent step."
					break
				end
			end
			if num_critical_loops > 0
			end
		end
		@logmsg loglevel1 "Exiting after $(num_critical_loops) loops with ω = $(ω) and Δ = $(get_delta(iter_data))."
	end

	ω, x₊, mx₊, steplength = compute_descent_step(mop, iter_data, data_base, sc, algo_config, ω, ω_data)

    if ω_Δ_rel_test(ω, Δ, algo_config) || ω_abs_test( ω, algo_config )
        return CRITICAL
    end

	mx = eval_models(sc, x)
	fx₊ = eval_all_objectives(mop, x₊)

	if strict_acceptance_test( algo_config )
		_ρ = minimum( (fx .- fx₊) ./ (mx .- mx₊) )
	else
		_ρ = (maximum(fx) - maximum( fx₊ ))/(maximum(mx) - maximum(mx₊))
	end
	ρ = isnan(_ρ) ? -Inf : _ρ

	@logmsg loglevel2 """\n
	Attempting descent of length $steplength.
	| f(x)  | $(_prettify(fx))
	| f(x₊) | $(_prettify(fx₊))
	| m(x)  | $(_prettify(mx))
	| m(x₊) | $(_prettify(mx₊))
	The error betwenn f(x) and m(x) is $(sum(abs.(fx .- mx))).
	$(strict_acceptance_test(algo_config) ? "All" : "One") of the components must decrease.
	Thus, ρ is $ρ.
	"""

	# put new values in data base
	trial_index = new_result!(data_base, x₊, fx₊)

	Δ_old = Δ	# save for printing
	if ρ >= ν_success
		if Δ < β * ω
			Δ = grow_radius( algo_config, Δ, steplength )
		end
		it_stat!( iter_data, SUCCESSFULL )
	else
		if fully_linear(sc)
			if ρ < ν_accept
				Δ = shrink_radius_much( algo_config, Δ, steplength )
				it_stat!( iter_data, INACCEPTABLE )
			else
				Δ = shrink_radius(algo_config, Δ, steplength)
				it_stat!( iter_data, ACCEPTABLE )
			end
		else
			it_stat!( iter_data, MODELIMPROVING )
		end
	end

	# before updating `iter_data`, retrieve a saveable
	# (order is important here to have current index and next index)
	stamp_content = get_saveable( SAVEABLE_TYPE, iter_data;
		x_trial_index = trial_index,
		ρ, ω, sc,
		stepsize = steplength )
	stamp!( data_base, stamp_content )

	# update iter data
	set_delta!(iter_data, Δ)
	if it_stat(iter_data) in [SUCCESSFULL, ACCEPTABLE]
		set_x!(iter_data, x₊)
		set_fx!(iter_data, fx₊)
		set_x_index!( iter_data, trial_index )
	end

	@logmsg loglevel1 """\n
		The trial point was $( (it_stat(iter_data) == SUCCESSFULL) || (it_stat(iter_data) == ACCEPTABLE ) ? "" : "not ")accepted.
		The iteration is $(it_stat(iter_data)).
		Moreover, the radius was updated as below:
		old radius : $(Δ_old)
		new radius : $(get_delta(iter_data)) ($(round(get_delta(iter_data)/Δ_old * 100; digits=1)) %)
	"""

	if ( x_tol_rel_test( x, x₊, algo_config  ) || x_tol_abs_test( x, x₊, algo_config ) ||
		f_tol_rel_test( fx, fx₊, algo_config  ) || f_tol_abs_test( fx, fx₊, algo_config ) )
		return TOLERANCE
	end

	return CONTINUE
end

############################################
function optimize( mop :: AbstractMOP, x0 :: Vec, fx0 :: Vec = Float32[];
    algo_config :: Union{AbstractConfig, Nothing} = nothing,
    populated_db :: Union{AbstractDB, Nothing} = nothing )

	mop, iter_data, data_base, sc, ac = initialize_data(mop, x0, fx0; algo_config, populated_db)
    @logmsg loglevel1 _stop_info_str( ac, mop )

    @logmsg loglevel1 "Entering main optimization loop."
	ret_code = CONTINUE
    while ret_code == CONTINUE
        ret_code = iterate!(iter_data, data_base, mop, sc, ac)
    end# while

	ret_x, ret_fx = get_return_values(data_base, iter_data, mop )

    @logmsg loglevel1 _fin_info_str(data_base, iter_data, mop, ret_code)

    return ret_x, ret_fx, ret_code, data_base
end# function optimize

end