text:
using Catlab.CategoricalAlgebra.CSets
using PermutationGroups

# Color assigned to each elem of each compoennt
const CDict = Dict{Symbol, Vector{Int}}

# The maximum color of an empty color list is 0
max0(x::Vector{Int})::Int = isempty(x) ? 0 : maximum(x)

check_auto(x::CDict)::Bool = all(map(Base.isperm, values(x)))

"""
Construct permutation σ⁻¹ such that σσ⁻¹=id
"""
function invert_perms(x::CDict)::CDict
  return Dict([k=>Base.invperm(v) for (k, v) in collect(x)])
end

is_perms(x::CDict)::Bool = all([Base.isperm(v) for v in values(x)])

function compose_comp(xs::Vector{Int},ys::Vector{Int})::Vector{Int}
  return [ys[xs[i]] for i in eachindex(xs)]
end

"""Compose permutations"""
function compose_perms(x::CDict, y::CDict)::CDict
  Dict([k => compose_comp(v1,y[k]) for (k, v1) in collect(x)])
end

"""Convert sorted permutation into a single permutation"""
function to_perm(d::CDict, csum::Vector{Int})::Perm
  res = vcat([[vi + offset for vi in v]
              for ((_, v), offset) in zip(sort(collect(d)), csum)]...)
  return Perm(res)
end

"""Convert single permutation back into sorted permutation"""
function from_perm(p::Perm, syms::Vector{Symbol}, lens::Vector{Int},
                   csum::Vector{Int})::CDict
  res = CDict()
  for (k, l, off) in zip(syms, lens, csum)
    res[k] = [p[i+off]-off for i in 1:l]
  end
  return res
end

"""Enumerate the elements of a permutation group from its generators"""
function all_perms(perm_gens::Vector{CDict})::Set{CDict}
  syms = sort(collect(keys(perm_gens[1])))
  lens = map(length, [perm_gens[1][x] for x in syms])
  csum = vcat([0], cumsum(lens)[1:end-1])
  _,_,Cs = schreier_sims([to_perm(g, csum) for g in perm_gens])
  result = Set{CDict}()
  """
  Recursive function based on (sub)chain C and partial product r
  Algorithm from Fig II.1 of Alexander Hulpke's notes on Computational
  Group Theory (Spring 2010).
  """
  function enum(i::Int, r)::Nothing
    leaf = i == length(Cs)
    C = Cs[i]
    for d ∈ C.orb
      xr = C[d]*r
      leaf ? push!(result, from_perm(xr, syms, lens, csum)) : enum(i+1, xr)
    end
  end

  enum(1, Perm(1:sum(lens)))

  return result
end

"""Apply a coloring to a C-set to get an isomorphic cset"""
function apply_automorphism(c::StructACSet{S}, d::CDict)::StructACSet{S} where {S}
  check_auto(d) || error("received coloring that is not an automorphism: $d")
  new = deepcopy(c)
  tabs, arrs, srcs, tgts = ob(S), hom(S), dom(S), codom(S)
  for (arr, src, tgt) in zip(arrs,srcs,tgts)
    set_subpart!(new, d[src], arr, d[tgt][c[arr]])
  end
  return new
end

"""
Suppose we have a set of candidates for "best permutation"
  (results in lexicographic min ACSet so far)

We get a new candidate that is a refinement of a previous candidate

Handle ambiguity by computing best case and worst case scenarios.

Returns
 - nothing: the two are not comparable
 - true: a < b
 - false: b < a
"""
const VUNI = Vector{Union{Nothing, Int}}

function compare_perms(a::Vector{VUNI}, b::Vector{VUNI})::Union{Nothing, Bool}
  for (va, vb) ∈ zip(a,b)
    for (ia, ib) in zip(va, vb)
      if isnothing(ia) || isnothing(ib)
        return nothing
      elseif ia < ib
        return true
      elseif ib < ia
        return false
      end
    end
  end
end

"""
Compute as much of the canonical data that will be ordered given a partial
coloring of an ACSet.

Given the data of a morphism f: A→B, the i'th element of f is: σᵦ(f(σₐ(i))).
With an arbitrary partition, σₐ(i) a subset of elements in |A|. Likewise for |B|
If the partitioning is not discrete, there will be cases where σₐ is not defined
We put nothing in this case.

If the best-case-scenario is worse than SOME automorphism's worst-case-scenario,
then this branch can be pruned.
"""
function order_perms(c::StructACSet, d::CDict, symorder::Vector{Symbol}
                    )::Pair{Vector{VUNI}, Vector{VUNI}}
  best, worst = map(collect, zip([order_perm(c,d,s) for s in symorder]...))
  best == worst || error("best $best worst $worst")
  return best => worst
end

function order_perm(c::StructACSet{S}, coloring::CDict, sym::Symbol
                    )::Pair{VUNI, VUNI} where {S}
  ind = findfirst(==(sym), hom(S))
  d, cd = dom(S)[ind], codom(S)[ind]
  res = [let x=perm_best_worst(c,coloring, sym, i, d, cd);
         isnothing(x) ? nothing=>nothing : x end for i in parts(c, cd)]
  bst, wrst = isempty(res) ? VUNI()=>VUNI() : map(collect, zip(res...))
  return bst => wrst
end

function perm_best_worst(c::StructACSet, coloring::CDict, f::Symbol, i::Int,
                         dm::Symbol, cdm::Symbol)::Union{Nothing, Pair{Int,Int}}
  σᵦinv = findall(==(i), coloring[cdm])
  finv = Set(incident(c, σᵦinv, f))
  σₐinv = findall(∈(finv), coloring[dm])
  return isempty(σₐinv) ? nothing : minimum(σₐinv) => maximum(σₐinv)
end

"""
Get vector listing nontrivial colors (which component and which color index) as
well as how many elements have that color. E.g. for (V=[1,1,2], E=[1,2,2,2,3,3])
we would get `[2=>(:V,1), 3=>(:E,2), 2=>(:E, 3)]`
"""
function get_colors_by_size(coloring::CDict)::Vector{Pair{Int,Tuple{Symbol, Int}}}
  res = []
  for (k, v) in coloring
    for color in 1:max0(v)
      n_c = count(==(color), v)
      if n_c > 1
        # Store which table and which color
        push!(res, n_c => (k, color))
      end
    end
  end
  return res
end


text:
# module for representing robot model
# modeling as differential wheeled robot

using Plots, Random, Distributions
pyplot()

include(joinpath(split(@__FILE__, "src")[1], "src/slam/graph_based_slam/logger_agent.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/slam/graph_based_slam/psi_sensor.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/slam/fast_slam_2/fast_slam_2_agent.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/model/agent/agent.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/model/agent/puddle_ignore_agent.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/model/agent/dp_policy_agent.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/model/agent/q_agent.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/model/sensor/sensor.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/common/error_calculation/error_calculation.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/common/state_transition/state_transition.jl"))

mutable struct DifferentialWheeledRobot
  pose
  radius
  color
  agent
  delta_time
  noise_expon
  noise_norm
  dist_until_noise
  bias_rate_spd
  bias_rate_yr
  stuck_exporn
  escape_exporn
  time_until_stuck
  time_until_escape
  kidnap_exporn
  kidnap_x_uniform
  kidnap_y_uniform
  kidnap_theta_uniform
  time_until_kidnap
  sensor
  is_stuck
  traj_x
  traj_y

  # init
  function DifferentialWheeledRobot(pose::Array, radius::Float64,
                                    color::String, agent,
                                    delta_time::Float64; sensor=nothing,
                                    noise_per_meter::Int64=0,
                                    noise_std::Float64=0.0,
                                    bias_rate_stds::Array=[0.0, 0.0],
                                    exp_stuck_time::Float64=1e100,
                                    exp_escape_time::Float64=1e-100,
                                    exp_kidnap_time::Float64=1e100,
                                    kidnap_rx::Array=[-5.0, 5.0],
                                    kidnap_ry::Array=[-5.0, 5.0])
      self = new()
      self.pose = pose
      self.radius = radius
      self.color = color
      self.agent = agent
      self.delta_time = delta_time
      self.noise_expon = Exponential(1.0/(1e-100 + noise_per_meter))
      self.noise_norm = Normal(0.0, noise_std)
      self.dist_until_noise = rand(self.noise_expon)
      self.bias_rate_spd = rand(Normal(1.0, bias_rate_stds[1]))
      self.bias_rate_yr = rand(Normal(1.0, bias_rate_stds[2]))
      self.stuck_exporn = Exponential(exp_stuck_time)
      self.escape_exporn = Exponential(exp_escape_time)
      self.time_until_stuck = rand(self.stuck_exporn)
      self.time_until_escape = rand(self.escape_exporn)
      self.kidnap_exporn = Exponential(exp_kidnap_time)
      self.kidnap_x_uniform = Uniform(kidnap_rx[1], kidnap_rx[2])
      self.kidnap_y_uniform = Uniform(kidnap_ry[1], kidnap_ry[2])
      self.kidnap_theta_uniform = Uniform(0.0, 2 * pi)
      self.time_until_kidnap = rand(self.kidnap_exporn)
      self.sensor = sensor
      self.is_stuck = false
      self.traj_x = [pose[1]]
      self.traj_y = [pose[2]]
      return self
  end
end

function noise(self::DifferentialWheeledRobot, pose::Array,
               speed::Float64, yaw_rate::Float64, delta_time::Float64)
  self.dist_until_noise -= abs(speed) * delta_time + self.radius * abs(yaw_rate) * delta_time
  if self.dist_until_noise <= 0.0
    self.dist_until_noise += rand(self.noise_expon)
    pose[3] += rand(self.noise_norm)
  end
  return pose
end

function bias(self::DifferentialWheeledRobot, speed::Float64,
              yaw_rate::Float64)
  return speed * self.bias_rate_spd, yaw_rate * self.bias_rate_yr
end

function stuck(self::DifferentialWheeledRobot, speed::Float64,
             yaw_rate::Float64, delta_time::Float64)
  if self.is_stuck
    self.time_until_escape -= delta_time
    if self.time_until_escape <= 0.0
      self.time_until_escape += rand(self.escape_exporn)
      self.is_stuck = false
    end
  else
    self.time_until_stuck -= delta_time
    if self.time_until_stuck <= 0.0
      self.time_until_stuck += rand(self.stuck_exporn)
      self.is_stuck = true
    end
  end
  return speed * (!self.is_stuck), yaw_rate * (!self.is_stuck)
end

function kidnap(self::DifferentialWheeledRobot, pose::Array,
                delta_time::Float64)
  self.time_until_kidnap -= delta_time
  if self.time_until_kidnap <= 0.0
    self.time_until_kidnap += rand(self.kidnap_exporn)
    return [rand(self.kidnap_x_uniform), rand(self.kidnap_y_uniform),
            rand(self.kidnap_theta_uniform)]
  else
    return pose
  end
end

function circle(x, y, r)
  theta = LinRange(0, 2*pi, 500)
  x .+ r * sin.(theta), y .+ r * cos.(theta)
end

function draw!(self::DifferentialWheeledRobot)
  # current pose
  x, y, theta = self.pose[1], self.pose[2], self.pose[3]
  xn = x + self.radius * cos(theta)
  yn = y + self.radius * sin(theta)

  # draw robot and current pose
  plot!([x, xn], [y, yn], color=self.color,
        legend=false, aspect_ratio=true)
  plot!(circle(x, y, self.radius), linecolor=self.color,
        aspect_ratio=true)

  # draw trajectory
  plot!(self.traj_x, self.traj_y, color="black",
        legend=false, aspect_ratio=true)

  # draw observation
  observation = []
  if self.sensor !== nothing
    data(self.sensor, self.pose)
    draw!(self.sensor, self.pose)
    observation = self.sensor.last_data
  end

  # draw estimation
  spd, yr = draw_decision!(self.agent, observation)

  if self.agent.estimator !== nothing
    draw_lon_lat_error!(self.agent.estimator.estimated_pose, self.pose)
  end

  # next pose
  spd, yr = bias(self, spd, yr)
  spd, yr = stuck(self, spd, yr, self.delta_time)
  self.pose = state_transition(spd, yr, self.delta_time, self.pose)
  self.pose = noise(self, self.pose, spd, yr, self.delta_time)
  self.pose = kidnap(self, self.pose, self.delta_time)
  push!(self.traj_x, self.pose[1]), push!(self.traj_y, self.pose[2])
end

ext: jl

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
include("constants.jl")

function angulares(λ, D; arcs = true)
    if arcs == true
        θ = 2.5e5 * (λ / D) #Arcseconds
    else
        θ = 1.22 * (λ / D) #Radians
    end

    return θ
end

function rms(λ)
    return 0.03 * λ
end

export angulares, rms


text:
export gradDescent!, updateParams!, optimizeParams!

"""

    gradDescent!(pars::Vector{<:Real}, grads::Vector{<:Real}, η=0.001) ->
    pars::Vector{<:Real}

Default gradient descent method in used in Quiqbox.
"""
function gradDescent!(pars::Vector{<:Real}, grads::Vector{<:Real}, η=0.001)
    @assert length(pars) == length(grads) "The length of gradients and corresponding "*
                                          "parameters should be the same."
    pars .-= η*grads
end

"""

    updateParams!(pbs::Array{<:ParamBox, 1}, grads::Array{<:Real, 1};
                  method::F=gradDescent!) where {F<:Function} -> Array{<:ParamBox, 1}

Given a `Vector` of parameters::`ParamBox` and its gradients with respect to each
parameter, update the `ParamBox`s and return the updated values.
"""
function updateParams!(pbs::Vector{<:ParamBox}, grads::Vector{<:Real};
                       method::F=gradDescent!) where {F<:Function}
    parVals = [i[] for i in pbs]
    method(parVals, grads)
    for (m,n) in zip(pbs, parVals)
        m[] = n
    end
    parVals
end

"""

    defaultECmethod(HFtype::Symbol, nuc::Array{String, 1},
                    nucCoords::Array{<:AbstractArray, 1}) ->
    fEC::Function

    defaultECmethod(Ne::NTuple{2, Int}, nuc::Array{String, 1},
                    nucCoords::Array{<:AbstractArray, 1}) ->
    fEC::Function

Return the default `Function` in `optimizeParams!` that will be used to update Hartree-Fock
energy and coefficient matrix(s):

    fEC(Hcore, HeeI, S) -> E::Float64,
                           C::Union{Array{Float64, 1}, NTuple{2, Array{Float64, 2}}}

When `HFtype` is set to  `:RHF` or input argument is only a 2-element `Tuple`, the returned
function is for RHF methods.
"""
function defaultECmethod(HFtype::Symbol, nuc::Vector{String},
                         nucCoords::Vector{<:AbstractArray}, Ne::Int=getCharge(nuc))
    HFtype == :UHF && ( Ne = (Ne÷2, Ne-Ne÷2) )
    defaultECmethodCore(Ne, nuc, nucCoords)
end

# Direct specify number of α and β electrons for UHF.
defaultECmethod(Nes::NTuple{2, Int},
                nuc::Vector{String}, nucCoords::Vector{<:AbstractArray}) =
defaultECmethodCore(Nes, nuc, nucCoords)

function defaultECmethodCore(Ne, nuc::Vector{String}, nucCoords::Vector{<:AbstractArray})
    function (Hcore, HeeI, bs, S)
        X = getX(S)
        res = runHFcore(defaultSCFconfig, Ne, Hcore, HeeI, S, X,
                        guessC(:SAD, (length(Ne)==2), S, X,
                               Hcore, HeeI, bs, nuc, nucCoords);
                        printInfo=false)
        res.E0HF, res.C
    end
end

"""

    optimizeParams!(bs::Array{<:FloatingGTBasisFuncs, 1}, pbs::Array{<:ParamBox, 1},
                    nuc::Array{String, 1}, nucCoords::Array{<:AbstractArray, 1};
                    Etarget::Float64=NaN, threshold::Float64=1e-4, maxSteps::Int=2000,
                    printInfo::Bool=true, GDmethod::F1=gradDescent!,
                    ECmethod::F2=Quiqbox.defaultECmethod(:RHF, nuc, nucCoords)) where
                   {F1<:Function, F2<:Function} ->
    Es::Array{Float64, 1}, pars::Array{Float64, 2}, grads::Array{Float64, 2}

The main function to optimize the parameters of a given basis set.

=== Positional argument(s) ===

`bs::Array{<:FloatingGTBasisFuncs, 1}`: Basis set.

`pbs::Array{<:ParamBox, 1}`: The parameters to be optimized that are extracted from the
basis set.

`nuc::Array{String, 1}`: The nuclei of the molecule.

`nucCoords::Array{<:AbstractArray, 1}`: The nuclei coordinates.

`ECmethod::F2`: The `Function` used to update Hartree-Fock energy and coefficient matrix(s)
during the optimization iterations. Default setting is
`Quiqbox.defaultECmethod(:RHF, nuc, nucCoords)` for RHF and the number of electron is equal
to nuclei charge.

=== Keyword argument(s) ===

`Etarget::Float64`: The target Hartree-Hock energy intent to achieve.

`threshold::Float64`: The threshold for the convergence when evaluating difference between
the latest two energies.

`maxSteps::Int`: Maximum allowed iteration steps regardless of whether the optimization
iteration converges.

`printInfo::Bool`: Whether print out the information of each iteration step.

`GDmethod::F1`: Applied gradient descent `Function`.
"""
function optimizeParams!(bs::Vector{<:FloatingGTBasisFuncs}, pbs::Vector{<:ParamBox},
                         nuc::Vector{String}, nucCoords::Vector{<:AbstractArray},
                         ECmethod::F2=defaultECmethod(:RHF, nuc, nucCoords);
                         Etarget::Float64=NaN, threshold::Float64=1e-4, maxSteps::Int=2000,
                         printInfo::Bool=true, GDmethod::F1=gradDescent!) where
                        {F1<:Function, F2<:Function}
    tAll = @elapsed begin

        i = 0
        Es = Float64[]
        pars = zeros(0, length(pbs))
        grads = zeros(0, length(pbs))
        gap = max(5, maxSteps ÷ 200 * 5)

        if Etarget === NaN
            isConverged = (Es) -> isOscillateConverged(Es, threshold, leastCycles=3)
        else
            isConverged = Es -> (abs(Es[end] - Etarget) < threshold)
        end

        parsL = [i[] for i in pbs]

        while true
            S = overlaps(bs)
            Hcore = coreH(bs, nuc, nucCoords)
            HeeI = eeInteractions(bs)
            E, C = ECmethod(Hcore, HeeI, bs, S)

            t = @elapsed begin
                grad = gradHFenegy(bs, pbs, C, S, nuc, nucCoords)
            end

            push!(Es, E)
            pars = vcat(pars, parsL |> transpose)
            grads = vcat(grads, grad |> transpose)

            if i%gap == 0 && printInfo
                println(rpad("Step $i: ", 15), rpad("E = $(E)", 26))
                print(rpad("", 10), "params = ")
                println(IOContext(stdout, :limit => true), parsL)
                print(rpad("", 12), "grad = ")
                println(IOContext(stdout, :limit => true), grad)
                println("Step duration: ", t, " seconds.\n")
            end

            parsL = updateParams!(pbs, grad, method=GDmethod)

            !isConverged(Es) && i < maxSteps || break

            i += 1
        end

        popfirst!(Es)

    end

    if printInfo
        println("The iteration just ended at")
        println(rpad("Step $(i): ", 15), rpad("E = $(Es[end])", 26))
        print(rpad("", 10), "params = ")
        println(IOContext(stdout, :limit => true), pars[end, :])
        print(rpad("", 12), "grad = ")
        println(IOContext(stdout, :limit => true), grads[end, :])
        println("Optimization duration: ", tAll/60, " minutes.")
        !isConverged(Es) && println("The result has not converged.")
    end

    Es, pars, grads
end


text:
include("mul.jl")
include("lazymul.jl")
include("blasbroadcasting.jl")
include("blasmul.jl")
include("inv.jl")


text:
module Compress

# TODO: import TOML for Julia 1.6
import Pkg.TOML
const STDLIB_TOML = VERSION >= v"1.6.0-DEV.764"
import Pkg.Types: VersionSpec, VersionRange, VersionBound

"""
    compress_versions(pool::Vector{VersionNumber}, subset::Vector{VersionNumber})

Given `pool` as the pool of available versions (of some package) and `subset` as some
subset of the pool of available versions, this function computes a `VersionSpec` which
includes all versions in `subset` and none of the versions in its complement.
"""
function compress_versions(pool::Vector{VersionNumber}, subset::Vector{VersionNumber})
    # Explicitly drop prerelease/build numbers, as those can confuse this.
    # TODO: Rewrite all this to use VersionNumbers instead of VersionBounds
    drop_build_prerelease(v::VersionNumber) = VersionNumber(v.major, v.minor, v.patch)
    pool = drop_build_prerelease.(pool)
    subset = sort!(drop_build_prerelease.(subset))

    complement = sort!(setdiff(pool, subset))
    ranges = VersionRange[]
    @label again
    isempty(subset) && return VersionSpec(ranges)
    a = first(subset)
    for b in reverse(subset)
        a.major == b.major || continue
        for m = 1:3
            lo = VersionBound((a.major, a.minor, a.patch)[1:m]...)
            for n = 1:3
                hi = VersionBound((b.major, b.minor, b.patch)[1:n]...)
                r = VersionRange(lo, hi)
                if !any(v in r for v in complement)
                    filter!(!in(r), subset)
                    push!(ranges, r)
                    @goto again
                end
            end
        end
    end
end
function compress_versions(pool::Vector{VersionNumber}, subset)
    compress_versions(pool, filter(in(subset), pool))
end

function load_versions(path::String)
    versions_file = joinpath(dirname(path), "Versions.toml")
    versions_dict = TOML.parsefile(versions_file)
    sort!([VersionNumber(v) for v in keys(versions_dict)])
end

function load(path::String,
              versions::Vector{VersionNumber} = load_versions(path))
    compressed = TOML.parsefile(path)
    uncompressed = Dict{VersionNumber,Dict{Any,Any}}()
    for (vers, data) in compressed
        vs = VersionSpec(vers)
        for v in versions
            v in vs || continue
            uv = get!(uncompressed, v, Dict())
            for (key, value) in data
                if haskey(uv, key)
                    error("Overlapping ranges for $(key) in Compat. Detected for version $(v).")
                else
                    uv[key] = value
                end
            end
        end
    end
    return uncompressed
end

function compress(path::String, uncompressed::Dict,
                  versions::Vector{VersionNumber} = load_versions(path))
    inverted = Dict{Pair{String,Any},Vector{VersionNumber}}()
    for (ver, data) in uncompressed, (key, val) in data
        val isa Base.UUID && (val = string(val))
        push!(get!(inverted, key => val, VersionNumber[]), ver)
    end
    compressed = Dict{String,Dict{String,Any}}()
    for ((k, v), vers) in inverted
        for r in compress_versions(versions, sort!(vers)).ranges
            get!(compressed, string(r), Dict{String,Any}())[k] = v
        end
    end
    return compressed
end

function save(path::String, uncompressed::Dict,
              versions::Vector{VersionNumber} = load_versions(path))
    compressed = compress(path, uncompressed)
    open(path, write=true) do io
        if STDLIB_TOML
            TOML.print(string, io, compressed, sorted=true)
        else
            TOML.print(io, compressed, sorted=true)
        end
    end
end

#=
=#

end # module


text:
using NavAbilitySDK
using Test
using JSON
using LinearAlgebra
using Random

include("./testDistributions.jl")
include("./testFactors.jl")
include("./testVariables.jl")


text:
abstract type AbstractDevice end

include("LabStream.jl")
include("AbstractDevice.jl")

##### General interface for LabStream
## A stream of type T <: LabStream should define the methods
# serialize(::T, cmd)                   # Send data `cmd`
# init_devices!(::T, devs::AbstractDevice...)   # Initialize devices and connect them to stream
# send(::T)                            # Send set! commands for all devices to stream, using getwritecommand(stream::LabStream, dev::AbstractDevice, getsetvalue(dev::AbstractDevice))
# read(::T)                             # Send read commands for all devices to stream, sing get(dev::AbstractDevice)
# init_devices!                         # Initialize all connected devises, using initialize(::AbstractDevice)
# close                                # Close connection, call close(::AbstractDevice) on all connected devices

#Include the stream definitions
include("ComediStream.jl")
include("BeagleBoneStream.jl")

##### General interface for devices
## A device to type T<:AbstractDevice should define the methods
# set!(dev::T, val)    # Remember val for next send(::LabStream) command
# get(dev::T)          # Get val from last next read(::LabStream) command

# getwritecommand(stream::LabStream, dev::T, val)   #Get canonlical representation of send command
# getreadcommand(stream::LabStream, dev::T)         #Get canonlical representation of read command

## Default methods:
# initialize(::AbstractDevice) = nothing    #What to do when initializing
# close(::AbstractDevice) = nothing          #What to do when disconnecting
# getstream(dev::AbstractDevice) = dev.stream #Which stream is the Device connected to
# setstream!(dev::AbstractDevice, stream::LabStream) = dev.stream = stream #Set the stream the Device is connected to
##### END General interface for devices

#Include the device definitions
include("SysLED.jl")
include("10V.jl")


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

## basic task functions and TLS

Core.Task(@nospecialize(f), reserved_stack::Int=0) = Core._Task(f, reserved_stack, ThreadSynchronizer())

# Container for a captured exception and its backtrace. Can be serialized.
struct CapturedException <: Exception
    ex::Any
    processed_bt::Vector{Any}

    function CapturedException(ex, bt_raw::Vector)
        # bt_raw MUST be a vector that can be processed by StackTraces.stacktrace
        # Typically the result of a catch_backtrace()

        # Process bt_raw so that it can be safely serialized
        bt_lines = process_backtrace(bt_raw, 100) # Limiting this to 100 lines.
        CapturedException(ex, bt_lines)
    end

    CapturedException(ex, processed_bt::Vector{Any}) = new(ex, processed_bt)
end

function showerror(io::IO, ce::CapturedException)
    showerror(io, ce.ex, ce.processed_bt, backtrace=true)
end

"""
    capture_exception(ex, bt) -> Exception

Returns an exception, possibly incorporating information from a backtrace `bt`. Defaults to returning [`CapturedException(ex, bt)`](@ref).

Used in [`asyncmap`](@ref) and [`asyncmap!`](@ref) to capture exceptions thrown during
the user-supplied function call.
"""
capture_exception(ex, bt) = CapturedException(ex, bt)

"""
    CompositeException

Wrap a `Vector` of exceptions thrown by a [`Task`](@ref) (e.g. generated from a remote worker over a channel
or an asynchronously executing local I/O write or a remote worker under `pmap`) with information about the series of exceptions.
For example, if a group of workers are executing several tasks, and multiple workers fail, the resulting `CompositeException` will
contain a "bundle" of information from each worker indicating where and why the exception(s) occurred.
"""
struct CompositeException <: Exception
    exceptions::Vector{Any}
    CompositeException() = new(Any[])
    CompositeException(exceptions) = new(exceptions)
end
length(c::CompositeException) = length(c.exceptions)
push!(c::CompositeException, ex) = push!(c.exceptions, ex)
pushfirst!(c::CompositeException, ex) = pushfirst!(c.exceptions, ex)
isempty(c::CompositeException) = isempty(c.exceptions)
iterate(c::CompositeException, state...) = iterate(c.exceptions, state...)
eltype(::Type{CompositeException}) = Any

function showerror(io::IO, ex::CompositeException)
    if !isempty(ex)
        showerror(io, ex.exceptions[1])
        remaining = length(ex) - 1
        if remaining > 0
            print(io, "\n\n...and ", remaining, " more exception", remaining > 1 ? "s" : "", ".\n")
        end
    else
        print(io, "CompositeException()\n")
    end
end

"""
    TaskFailedException

This exception is thrown by a `wait(t)` call when task `t` fails.
`TaskFailedException` wraps the failed task `t`.
"""
struct TaskFailedException <: Exception
    task::Task
end

function showerror(io::IO, ex::TaskFailedException, bt = nothing; backtrace=true)
    print(io, "TaskFailedException")
    if bt !== nothing && backtrace
        show_backtrace(io, bt)
    end
    println(io)
    printstyled(io, "\n    nested task error: ", color=error_color())
    show_task_exception(io, ex.task)
end

function show_task_exception(io::IO, t::Task; indent = true)
    stack = current_exceptions(t)
    b = IOBuffer()
    if isempty(stack)
        # exception stack buffer not available; probably a serialized task
        showerror(IOContext(b, io), t.result)
    else
        show_exception_stack(IOContext(b, io), stack)
    end
    str = String(take!(b))
    if indent
        str = replace(str, "\n" => "\n    ")
    end
    print(io, str)
end

function show(io::IO, t::Task)
    print(io, "Task ($(t.state)) @0x$(string(convert(UInt, pointer_from_objref(t)), base = 16, pad = Sys.WORD_SIZE>>2))")
end

"""
    @task

Wrap an expression in a [`Task`](@ref) without executing it, and return the [`Task`](@ref). This only
creates a task, and does not run it.

# Examples
```jldoctest
julia> a1() = sum(i for i in 1:1000);

julia> b = @task a1();

julia> istaskstarted(b)
false

julia> schedule(b);

julia> yield();

julia> istaskdone(b)
true
```
"""
macro task(ex)
    :(Task(()->$(esc(ex))))
end

"""
    current_task()

Get the currently running [`Task`](@ref).
"""
current_task() = ccall(:jl_get_current_task, Ref{Task}, ())

# task states

const task_state_runnable = UInt8(0)
const task_state_done     = UInt8(1)
const task_state_failed   = UInt8(2)

const _state_index = findfirst(==(:_state), fieldnames(Task))
@eval function load_state_acquire(t)
    # TODO: Replace this by proper atomic operations when available
    @GC.preserve t llvmcall($("""
        %ptr = inttoptr i$(Sys.WORD_SIZE) %0 to i8*
        %rv = load atomic i8, i8* %ptr acquire, align 8
        ret i8 %rv
        """), UInt8, Tuple{Ptr{UInt8}},
        Ptr{UInt8}(pointer_from_objref(t) + fieldoffset(Task, _state_index)))
end

@inline function getproperty(t::Task, field::Symbol)
    if field === :state
        # TODO: this field name should be deprecated in 2.0
        st = load_state_acquire(t)
        if st === task_state_runnable
            return :runnable
        elseif st === task_state_done
            return :done
        elseif st === task_state_failed
            return :failed
        else
            @assert false
        end
    elseif field === :backtrace
        # TODO: this field name should be deprecated in 2.0
        return current_exceptions(t)[end][2]
    elseif field === :exception
        # TODO: this field name should be deprecated in 2.0
        return t._isexception ? t.result : nothing
    else
        return getfield(t, field)
    end
end

"""
    istaskdone(t::Task) -> Bool

Determine whether a task has exited.

# Examples
```jldoctest
julia> a2() = sum(i for i in 1:1000);

julia> b = Task(a2);

julia> istaskdone(b)
false

julia> schedule(b);

julia> yield();

julia> istaskdone(b)
true
```
"""
istaskdone(t::Task) = load_state_acquire(t) !== task_state_runnable

"""
    istaskstarted(t::Task) -> Bool

Determine whether a task has started executing.

# Examples
```jldoctest
julia> a3() = sum(i for i in 1:1000);

julia> b = Task(a3);

julia> istaskstarted(b)
false
```
"""
istaskstarted(t::Task) = ccall(:jl_is_task_started, Cint, (Any,), t) != 0

"""
    istaskfailed(t::Task) -> Bool

Determine whether a task has exited because an exception was thrown.

# Examples
```jldoctest
julia> a4() = error("task failed");

julia> b = Task(a4);

julia> istaskfailed(b)
false

julia> schedule(b);

julia> yield();

julia> istaskfailed(b)
true
```

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
istaskfailed(t::Task) = (load_state_acquire(t) === task_state_failed)

Threads.threadid(t::Task) = Int(ccall(:jl_get_task_tid, Int16, (Any,), t)+1)

task_result(t::Task) = t.result

task_local_storage() = get_task_tls(current_task())
function get_task_tls(t::Task)
    if t.storage === nothing
        t.storage = IdDict()
    end
    return (t.storage)::IdDict{Any,Any}
end

"""
    task_local_storage(key)

Look up the value of a key in the current task's task-local storage.
"""
task_local_storage(key) = task_local_storage()[key]

"""
    task_local_storage(key, value)

Assign a value to a key in the current task's task-local storage.
"""
task_local_storage(key, val) = (task_local_storage()[key] = val)

"""
    task_local_storage(body, key, value)

Call the function `body` with a modified task-local storage, in which `value` is assigned to
`key`; the previous value of `key`, or lack thereof, is restored afterwards. Useful
for emulating dynamic scoping.
"""
function task_local_storage(body::Function, key, val)
    tls = task_local_storage()
    hadkey = haskey(tls, key)
    old = get(tls, key, nothing)
    tls[key] = val
    try
        return body()
    finally
        hadkey ? (tls[key] = old) : delete!(tls, key)
    end
end

# just wait for a task to be done, no error propagation
function _wait(t::Task)
    if !istaskdone(t)
        lock(t.donenotify)
        try
            while !istaskdone(t)
                wait(t.donenotify)
            end
        finally
            unlock(t.donenotify)
        end
    end
    nothing
end

# have `waiter` wait for `t`
function _wait2(t::Task, waiter::Task)
    if !istaskdone(t)
        lock(t.donenotify)
        if !istaskdone(t)
            push!(t.donenotify.waitq, waiter)
            unlock(t.donenotify)
            # since _wait2 is similar to schedule, we should observe the sticky
            # bit, even if we aren't calling `schedule` due to this early-return
            if waiter.sticky && Threads.threadid(waiter) == 0
                # Issue #41324
                # t.sticky && tid == 0 is a task that needs to be co-scheduled with
                # the parent task. If the parent (current_task) is not sticky we must
                # set it to be sticky.
                # XXX: Ideally we would be able to unset this
                current_task().sticky = true
                tid = Threads.threadid()
                ccall(:jl_set_task_tid, Cint, (Any, Cint), waiter, tid-1)
            end
            return nothing
        else
            unlock(t.donenotify)
        end
    end
    schedule(waiter)
    nothing
end

function wait(t::Task)
    t === current_task() && error("deadlock detected: cannot wait on current task")
    _wait(t)
    if istaskfailed(t)
        throw(TaskFailedException(t))
    end
    nothing
end

fetch(@nospecialize x) = x

"""
    fetch(t::Task)

Wait for a Task to finish, then return its result value.
If the task fails with an exception, a `TaskFailedException` (which wraps the failed task)
is thrown.
"""
function fetch(t::Task)
    wait(t)
    return task_result(t)
end

## lexically-scoped waiting for multiple items

struct ScheduledAfterSyncException <: Exception
    values::Vector{Any}
end

function showerror(io::IO, ex::ScheduledAfterSyncException)
    print(io, "ScheduledAfterSyncException: ")
    if isempty(ex.values)
        print(io, "(no values)")
        return
    end
    show(io, ex.values[1])
    if length(ex.values) == 1
        print(io, " is")
    elseif length(ex.values) == 2
        print(io, " and one more ")
        print(io, nameof(typeof(ex.values[2])))
        print(io, " are")
    else
        print(io, " and ", length(ex.values) - 1, " more objects are")
    end
    print(io, " registered after the end of a `@sync` block")
end

function sync_end(c::Channel{Any})
    local c_ex
    while isready(c)
        r = take!(c)
        if isa(r, Task)
            _wait(r)
            if istaskfailed(r)
                if !@isdefined(c_ex)
                    c_ex = CompositeException()
                end
                push!(c_ex, TaskFailedException(r))
            end
        else
            try
                wait(r)
            catch e
                if !@isdefined(c_ex)
                    c_ex = CompositeException()
                end
                push!(c_ex, e)
            end
        end
    end
    close(c)

    # Capture all waitable objects scheduled after the end of `@sync` and
    # include them in the exception. This way, the user can check what was
    # scheduled by examining at the exception object.
    local racy
    for r in c
        if !@isdefined(racy)
            racy = []
        end
        push!(racy, r)
    end
    if @isdefined(racy)
        if !@isdefined(c_ex)
            c_ex = CompositeException()
        end
        # Since this is a clear programming error, show this exception first:
        pushfirst!(c_ex, ScheduledAfterSyncException(racy))
    end

    if @isdefined(c_ex)
        throw(c_ex)
    end
    nothing
end

const sync_varname = gensym(:sync)

"""
    @sync

Wait until all lexically-enclosed uses of `@async`, `@spawn`, `@spawnat` and `@distributed`
are complete. All exceptions thrown by enclosed async operations are collected and thrown as
a `CompositeException`.
"""
macro sync(block)
    var = esc(sync_varname)
    quote
        let $var = Channel(Inf)
            v = $(esc(block))
            sync_end($var)
            v
        end
    end
end

# schedule an expression to run asynchronously

"""
    @async

Wrap an expression in a [`Task`](@ref) and add it to the local machine's scheduler queue.

Values can be interpolated into `@async` via `\$`, which copies the value directly into the
constructed underlying closure. This allows you to insert the _value_ of a variable,
isolating the asynchronous code from changes to the variable's value in the current task.

!!! warning
    It is strongly encouraged to favor `Threads.@spawn` over `@async` always **even when no
    parallelism is required** especially in publicly distributed libraries.  This is
    because a use of `@async` disables the migration of the *parent* task across worker
    threads in the current implementation of Julia.  Thus, seemingly innocent use of
    `@async` in a library function can have a large impact on the performance of very
    different parts of user applications.

!!! compat "Julia 1.4"
    Interpolating values via `\$` is available as of Julia 1.4.
"""
macro async(expr)
    do_async_macro(expr)
end

# generate the code for @async, possibly wrapping the task in something before
# pushing it to the wait queue.
function do_async_macro(expr; wrap=identity)
    letargs = Base._lift_one_interp!(expr)

    thunk = esc(:(()->($expr)))
    var = esc(sync_varname)
    quote
        let $(letargs...)
            local task = Task($thunk)
            if $(Expr(:islocal, var))
                put!($var, $(wrap(:task)))
            end
            schedule(task)
            task
        end
    end
end

# task wrapper that doesn't create exceptions wrapped in TaskFailedException
struct UnwrapTaskFailedException
    task::Task
end

# common code for wait&fetch for UnwrapTaskFailedException
function unwrap_task_failed(f::Function, t::UnwrapTaskFailedException)
    try
        f(t.task)
    catch ex
        if ex isa TaskFailedException
            throw(ex.task.exception)
        else
            rethrow()
        end
    end
end

# the unwrapping for above task wrapper (gets triggered in sync_end())
wait(t::UnwrapTaskFailedException) = unwrap_task_failed(wait, t)

# same for fetching the tasks, for convenience
fetch(t::UnwrapTaskFailedException) = unwrap_task_failed(fetch, t)

# macro for running async code that doesn't throw wrapped exceptions
macro async_unwrap(expr)
    do_async_macro(expr, wrap=task->:(Base.UnwrapTaskFailedException($task)))
end

"""
    errormonitor(t::Task)

Print an error log to `stderr` if task `t` fails.
"""
function errormonitor(t::Task)
    t2 = Task() do
        if istaskfailed(t)
            local errs = stderr
            try # try to display the failure atomically
                errio = IOContext(PipeBuffer(), errs::IO)
                emphasize(errio, "Unhandled Task ")
                display_error(errio, scrub_repl_backtrace(current_exceptions(t)))
                write(errs, errio)
            catch
                try # try to display the secondary error atomically
                    errio = IOContext(PipeBuffer(), errs::IO)
                    print(errio, "\nSYSTEM: caught exception while trying to print a failed Task notice: ")
                    display_error(errio, scrub_repl_backtrace(current_exceptions()))
                    write(errs, errio)
                    flush(errs)
                    # and then the actual error, as best we can
                    Core.print(Core.stderr, "while handling: ")
                    Core.println(Core.stderr, current_exceptions(t)[end][1])
                catch e
                    # give up
                    Core.print(Core.stderr, "\nSYSTEM: caught exception of type ", typeof(e).name.name,
                            " while trying to print a failed Task notice; giving up\n")
                end
            end
        end
        nothing
    end
    t2.sticky = false
    _wait2(t, t2)
    return t
end

# Capture interpolated variables in $() and move them to let-block
function _lift_one_interp!(e)
    letargs = Any[]  # store the new gensymed arguments
    _lift_one_interp_helper(e, false, letargs) # Start out _not_ in a quote context (false)
    letargs
end
_lift_one_interp_helper(v, _, _) = v
function _lift_one_interp_helper(expr::Expr, in_quote_context, letargs)
    if expr.head === :$
        if in_quote_context  # This $ is simply interpolating out of the quote
            # Now, we're out of the quote, so any _further_ $ is ours.
            in_quote_context = false
        else
            newarg = gensym()
            push!(letargs, :($(esc(newarg)) = $(esc(expr.args[1]))))
            return newarg  # Don't recurse into the lifted $() exprs
        end
    elseif expr.head === :quote
        in_quote_context = true   # Don't try to lift $ directly out of quotes
    elseif expr.head === :macrocall
        return expr  # Don't recur into macro calls, since some other macros use $
    end
    for (i,e) in enumerate(expr.args)
        expr.args[i] = _lift_one_interp_helper(e, in_quote_context, letargs)
    end
    expr
end

# add a wait-able object to the sync pool
macro sync_add(expr)
    var = esc(sync_varname)
    quote
        local ref = $(esc(expr))
        put!($var, ref)
        ref
    end
end

# runtime system hook called when a task finishes
function task_done_hook(t::Task)
    # `finish_task` sets `sigatomic` before entering this function
    err = istaskfailed(t)
    result = task_result(t)
    handled = false

    donenotify = t.donenotify
    if isa(donenotify, ThreadSynchronizer)
        lock(donenotify)
        try
            if !isempty(donenotify.waitq)
                handled = true
                notify(donenotify)
            end
        finally
            unlock(donenotify)
        end
    end

    if err && !handled && Threads.threadid() == 1
        if isa(result, InterruptException) && isdefined(Base, :active_repl_backend) &&
            active_repl_backend.backend_task._state === task_state_runnable && isempty(Workqueue) &&
            active_repl_backend.in_eval
            throwto(active_repl_backend.backend_task, result) # this terminates the task
        end
    end
    # Clear sigatomic before waiting
    sigatomic_end()
    try
        wait() # this will not return
    catch e
        # If an InterruptException happens while blocked in the event loop, try handing
        # the exception to the REPL task since the current task is done.
        # issue #19467
        if Threads.threadid() == 1 &&
            isa(e, InterruptException) && isdefined(Base, :active_repl_backend) &&
            active_repl_backend.backend_task._state === task_state_runnable && isempty(Workqueue) &&
            active_repl_backend.in_eval
            throwto(active_repl_backend.backend_task, e)
        else
            rethrow()
        end
    end
end

## scheduler and work queue

struct InvasiveLinkedListSynchronized{T}
    queue::InvasiveLinkedList{T}
    lock::Threads.SpinLock
    InvasiveLinkedListSynchronized{T}() where {T} = new(InvasiveLinkedList{T}(), Threads.SpinLock())
end
isempty(W::InvasiveLinkedListSynchronized) = isempty(W.queue)
length(W::InvasiveLinkedListSynchronized) = length(W.queue)
function push!(W::InvasiveLinkedListSynchronized{T}, t::T) where T
    lock(W.lock)
    try
        push!(W.queue, t)
    finally
        unlock(W.lock)
    end
    return W
end
function pushfirst!(W::InvasiveLinkedListSynchronized{T}, t::T) where T
    lock(W.lock)
    try
        pushfirst!(W.queue, t)
    finally
        unlock(W.lock)
    end
    return W
end
function pop!(W::InvasiveLinkedListSynchronized)
    lock(W.lock)
    try
        return pop!(W.queue)
    finally
        unlock(W.lock)
    end
end
function popfirst!(W::InvasiveLinkedListSynchronized)
    lock(W.lock)
    try
        return popfirst!(W.queue)
    finally
        unlock(W.lock)
    end
end
function list_deletefirst!(W::InvasiveLinkedListSynchronized{T}, t::T) where T
    lock(W.lock)
    try
        list_deletefirst!(W.queue, t)
    finally
        unlock(W.lock)
    end
    return W
end

const StickyWorkqueue = InvasiveLinkedListSynchronized{Task}
global const Workqueues = [StickyWorkqueue()]
global const Workqueue = Workqueues[1] # default work queue is thread 1
function __preinit_threads__()
    nt = Threads.nthreads()
    if length(Workqueues) < nt
        resize!(Workqueues, nt)
        for i = 2:nt
            Workqueues[i] = StickyWorkqueue()
        end
    end
    Partr.multiq_init(nt)
    nothing
end

function enq_work(t::Task)
    (t._state === task_state_runnable && t.queue === nothing) || error("schedule: Task not runnable")
    tid = Threads.threadid(t)
    # Note there are three reasons a Task might be put into a sticky queue
    # even if t.sticky == false:
    # 1. The Task's stack is currently being used by the scheduler for a certain thread.
    # 2. There is only 1 thread.
    # 3. The multiq is full (can be fixed by making it growable).
    if t.sticky || Threads.nthreads() == 1
        if tid == 0
            # Issue #41324
            # t.sticky && tid == 0 is a task that needs to be co-scheduled with
            # the parent task. If the parent (current_task) is not sticky we must
            # set it to be sticky.
            # XXX: Ideally we would be able to unset this
            current_task().sticky = true
            tid = Threads.threadid()
            ccall(:jl_set_task_tid, Cint, (Any, Cint), t, tid-1)
        end
        push!(Workqueues[tid], t)
    else
        if !Partr.multiq_insert(t, t.priority)
            # if multiq is full, give to a random thread (TODO fix)
            if tid == 0
                tid = mod(time_ns() % Int, Threads.nthreads()) + 1
                ccall(:jl_set_task_tid, Cint, (Any, Cint), t, tid-1)
            end
            push!(Workqueues[tid], t)
        else
            tid = 0
        end
    end
    ccall(:jl_wakeup_thread, Cvoid, (Int16,), (tid - 1) % Int16)
    return t
end

schedule(t::Task) = enq_work(t)

"""
    schedule(t::Task, [val]; error=false)

Add a [`Task`](@ref) to the scheduler's queue. This causes the task to run constantly when the system
is otherwise idle, unless the task performs a blocking operation such as [`wait`](@ref).

If a second argument `val` is provided, it will be passed to the task (via the return value of
[`yieldto`](@ref)) when it runs again. If `error` is `true`, the value is raised as an exception in
the woken task.

!!! warning
    It is incorrect to use `schedule` on an arbitrary `Task` that has already been started.
    See [the API reference](@ref low-level-schedule-wait) for more information.

# Examples
```jldoctest
julia> a5() = sum(i for i in 1:1000);

julia> b = Task(a5);

julia> istaskstarted(b)
false

julia> schedule(b);

julia> yield();

julia> istaskstarted(b)
true

julia> istaskdone(b)
true
```
"""
function schedule(t::Task, @nospecialize(arg); error=false)
    # schedule a task to be (re)started with the given value or exception
    t._state === task_state_runnable || Base.error("schedule: Task not runnable")
    if error
        t.queue === nothing || Base.list_deletefirst!(t.queue, t)
        setfield!(t, :result, arg)
        setfield!(t, :_isexception, true)
    else
        t.queue === nothing || Base.error("schedule: Task not runnable")
        setfield!(t, :result, arg)
    end
    enq_work(t)
    return t
end

"""
    yield()

Switch to the scheduler to allow another scheduled task to run. A task that calls this
function is still runnable, and will be restarted immediately if there are no other runnable
tasks.
"""
function yield()
    ct = current_task()
    enq_work(ct)
    try
        wait()
    catch
        ct.queue === nothing || list_deletefirst!(ct.queue, ct)
        rethrow()
    end
end

@inline set_next_task(t::Task) = ccall(:jl_set_next_task, Cvoid, (Any,), t)

"""
    yield(t::Task, arg = nothing)

A fast, unfair-scheduling version of `schedule(t, arg); yield()` which
immediately yields to `t` before calling the scheduler.
"""
function yield(t::Task, @nospecialize(x=nothing))
    (t._state === task_state_runnable && t.queue === nothing) || error("yield: Task not runnable")
    t.result = x
    enq_work(current_task())
    set_next_task(t)
    return try_yieldto(ensure_rescheduled)
end

"""
    yieldto(t::Task, arg = nothing)

Switch to the given task. The first time a task is switched to, the task's function is
called with no arguments. On subsequent switches, `arg` is returned from the task's last
call to `yieldto`. This is a low-level call that only switches tasks, not considering states
or scheduling in any way. Its use is discouraged.
"""
function yieldto(t::Task, @nospecialize(x=nothing))
    # TODO: these are legacy behaviors; these should perhaps be a scheduler
    # state error instead.
    if t._state === task_state_done
        return x
    elseif t._state === task_state_failed
        throw(t.result)
    end
    t.result = x
    set_next_task(t)
    return try_yieldto(identity)
end

function try_yieldto(undo)
    try
        ccall(:jl_switch, Cvoid, ())
    catch
        undo(ccall(:jl_get_next_task, Ref{Task}, ()))
        rethrow()
    end
    ct = current_task()
    if ct._isexception
        exc = ct.result
        ct.result = nothing
        ct._isexception = false
        throw(exc)
    end
    result = ct.result
    ct.result = nothing
    return result
end

# yield to a task, throwing an exception in it
function throwto(t::Task, @nospecialize exc)
    t.result = exc
    t._isexception = true
    set_next_task(t)
    return try_yieldto(identity)
end

function ensure_rescheduled(othertask::Task)
    ct = current_task()
    W = Workqueues[Threads.threadid()]
    if ct !== othertask && othertask._state === task_state_runnable
        # we failed to yield to othertask
        # return it to the head of a queue to be retried later
        tid = Threads.threadid(othertask)
        Wother = tid == 0 ? W : Workqueues[tid]
        pushfirst!(Wother, othertask)
    end
    # if the current task was queued,
    # also need to return it to the runnable state
    # before throwing an error
    list_deletefirst!(W, ct)
    nothing
end

function trypoptask(W::StickyWorkqueue)
    while !isempty(W)
        t = popfirst!(W)
        if t._state !== task_state_runnable
            # assume this somehow got queued twice,
            # probably broken now, but try discarding this switch and keep going
            # can't throw here, because it's probably not the fault of the caller to wait
            # and don't want to use print() here, because that may try to incur a task switch
            ccall(:jl_safe_printf, Cvoid, (Ptr{UInt8}, Int32...),
                "\nWARNING: Workqueue inconsistency detected: popfirst!(Workqueue).state != :runnable\n")
            continue
        end
        return t
    end
    return Partr.multiq_deletemin()
end

function checktaskempty()
    return Partr.multiq_check_empty()
end

@noinline function poptask(W::StickyWorkqueue)
    task = trypoptask(W)
    if !(task isa Task)
        task = ccall(:jl_task_get_next, Ref{Task}, (Any, Any, Any), trypoptask, W, checktaskempty)
    end
    set_next_task(task)
    nothing
end

function wait()
    GC.safepoint()
    W = Workqueues[Threads.threadid()]
    poptask(W)
    result = try_yieldto(ensure_rescheduled)
    process_events()
    # return when we come out of the queue
    return result
end

if Sys.iswindows()
    pause() = ccall(:Sleep, stdcall, Cvoid, (UInt32,), 0xffffffff)
else
    pause() = ccall(:pause, Cvoid, ())
end


text:
abstract type AbstractCatabolism end

"""
    catabolism!(o, u)

Where `o` is an `Organ` and `u` its state variables

Catabolism for E, C and N, or C, N and E reserves.
"""
catabolism!(o, u) = catabolism!(catabolism_pars(o), o, u)

kC(p::AbstractCatabolism) = p.kC
kN(p::AbstractCatabolism) = p.kN
kE(p::AbstractCatabolism) = p.k

abstract type AbstractCatabolismCN <: AbstractCatabolism end

catabolism!(p::AbstractCatabolismCN, o, u) = begin
    v, J = vars(o), flux(o)
    turnover = (kC(p), kN(p)) .* tempcorrection(v) .* scaling(v)
    rel_reserve = (u[:C], u[:N]) ./ u[:V]
    corrected_j_E_mai = j_E_mai(o) * tempcorrection(v)

    # Calculate the growth rate
    r, alive = calc_rate(su_pars(o), rel_reserve, turnover, corrected_j_E_mai,
                         y_E_C(o), y_E_N(o), y_V_E(o), κsoma(o), tstep(o))
    set_rate!(o, r)

    # Catabolise reserve depending on reserve state and turnover minus the growth rate
    C_ctb, N_ctb = (u[:C], u[:N]) .* (turnover .- r)
    C_rej, N_rej, E_ctb = stoich_merge(su_pars(o), C_ctb, N_ctb, y_E_C(o), y_E_N(o))
    # Set catabolism var for use in other components
    set_E_ctb!(o, E_ctb)
    # Set rejected catabolised reserves to be translocated, or lost
    J[:C, :rej] = -C_rej
    J[:N, :rej] = -N_rej
    alive
end

"""
    CatabolismCN(kC, kN)

2-pareter catabolism where reserve turnover rate can be different for C and N reserves.

$(FIELDDOCTABLE)
"""
@columns struct CatabolismCN{KC,KN} <: AbstractCatabolismCN
    kC::KC | 0.2 | mol*mol^-1*d^-1 | (0.0,1.0) | _   | "C-reserve turnover rate"
    kN::KN | 0.2 | mol*mol^-1*d^-1 | (0.0,1.0) | _   | "N-reserve turnover rate"
end

"""
    CatabolismCNshared(k)

1-pararameter catabolism where reserve turnover rate is the same for both reserves.

$(FIELDDOCTABLE)
"""
@columns struct CatabolismCNshared{K} <: AbstractCatabolismCN
    k::K  | 0.2 | mol*mol^-1*d^-1 | (0.0,1.0) | _   | "Reserve turnover rate"
end

kC(p::CatabolismCNshared) = p.k
kN(p::CatabolismCNshared) = p.k

abstract type AbstractCatabolismCNE <: AbstractCatabolism end

catabolism!(p::AbstractCatabolismCNE, o, u) = begin
    v, J, J1 = vars(o), flux(o), flux1(o)
    turnover = (kC(p), kN(p), kE(p)) .* tempcorrection(v) .* scaling(v)
    reserve = (u[:C], u[:N], u[:E])
    rel_reserve = reserve ./ u[:V]
    corrected_j_E_mai = j_E_mai(o) * tempcorrection(o)

    r, isalive = calc_rate(su_pars(o), rel_reserve, turnover, corrected_j_E_mai(o), y_E_C(o), y_E_N(o), y_V_E(o), κsoma(o), tstep(o))
    set_rate!(o, r)

    J1_C_ctb, J1_N_ctb, J_EE_ctb = reserve .* (turnover .- r)
    C_rej, N_rej, J1_CN_ctb = stoich_merge(su_pars(o), J_C_ctb, J1_N_ctb, y_E_C(o), y_E_N(o))
    E_ctb = J1_EE_ctb + J1_CN_ctb # Total catabolic flux

    # Set rejected flux
    J[:N,:rej] = -C_rej
    J[:C,:rej] = -N_rej
    # Set catabolism vars for use in other components
    set_θE!(o, J_EE_ctb / J_E_ctb) # Proportion of general reserve flux in total catabolic flux
    set_E_ctb!(o, J_E_ctb)
    isalive
end

"""
    CatabolismCNshared(k, kC, kN)

3-pareter catabolism where reserve turnover rate can be different for C, N and E reserves.

$(FIELDDOCTABLE)
"""
@columns struct CatabolismCNE{KE,KC,KN} <: AbstractCatabolismCNE
    k::KE  | 0.2 | mol*mol^-1*d^-1 | (0.0,1.0) | _   | "Reserve turnover rate"
    kC::KC | 0.2 | mol*mol^-1*d^-1 | (0.0,1.0) | _   | "C-reserve turnover rate"
    kN::KN | 0.2 | mol*mol^-1*d^-1 | (0.0,1.0) | _   | "N-reserve turnover rate"
end


text:
module BlueNoiseGeneration

using FFTW
using Plots
using Random
using StaticArrays

struct Parameters
    starting_density::Float64
    σ::Float64
    neighborhood_size::Int
end

Parameters() = Parameters(1 / 10, 2.0, 128)

function voidandcluster_energy_kernel(
    parameters::Parameters, dims::NTuple{D,Int}
)::Array{Float64,D} where {D}
    result = zeros(Float64, dims)
    for other_coordinate in CartesianIndices(dims)
        axis_distances = abs.(Tuple(other_coordinate) .- 1)
        # "Toroidal distance", i.e. pick which is closer between going to the
        # point the normal way or going to the point by wrapping around the
        # other direction:
        #  x-----x    distance: 6
        # -x     x-   distance: 3
        # 123456789
        axis_distances = min.(axis_distances, dims .- axis_distances)
        distance_squared = sum(axis_distances.^2)
        # Gaussian filter.
        result[other_coordinate] = exp(distance_squared / (-2 * parameters.σ^2))
    end
    result
end

const REDUCE_THRESHOLD = 1e-5

# struct ReducedKernel{D,T}
#     taps::SVector{T,Tuple{CartesianIndex{D},Float64}}
# end

struct ReducedKernel{D,T}
    taps::Vector{Tuple{CartesianIndex{D},Float64}}
end

function reduce_kernel(full::Array{Float64,D}) where D
    taps::Vector{Tuple{CartesianIndex{D},Float64}} = []
    for index in CartesianIndices(full)
        value_here = full[index]
        if value_here >= REDUCE_THRESHOLD
            push!(taps, (CartesianIndex(Tuple(index) .- 1), value_here))
        end
    end
    reduction = prod(size(full)) ÷ length(taps)
    println("Reduced kernel size x$reduction (now $(length(taps)) taps)")
    # num_taps = length(taps)
    # staps = SVector{num_taps,Tuple{CartesianIndex{D},Float64}}(taps)
    # ReducedKernel(staps)
    ReducedKernel{D,0}(taps)
end

function apply_kernel!(
    target::Array{Float64,D},
    kernel::ReducedKernel{D,T},
    offset::CartesianIndex{D},
    ::Val{A}
) where {D,T,A}
    target_size = size(target)
    for (tap_offset::CartesianIndex{D}, tap_value::Float64) in kernel.taps
        target_pos = mod1.(Tuple(tap_offset) .+ Tuple(offset), target_size)
        target_pos = CartesianIndex(target_pos)
        if A
            target[target_pos] += tap_value
        else
            target[target_pos] -= tap_value
        end
    end
end

add_kernel!(target::Array{Float64,D}, kernel::ReducedKernel{D}, offset::CartesianIndex{D}) where D =
    apply_kernel!(target, kernel, offset, Val(true))

sub_kernel!(target::Array{Float64,D}, kernel::ReducedKernel{D}, offset::CartesianIndex{D}) where D =
    apply_kernel!(target, kernel, offset, Val(false))

function generate(parameters::Parameters, size::NTuple{D,Int}) where {D}
    binary_pattern = Array{Bool,D}(undef, size...)
    energy = zeros(Float64, size)
    kernel = reduce_kernel(voidandcluster_energy_kernel(parameters, size))
    indices = CartesianIndices(binary_pattern)
    ns = parameters.neighborhood_size ÷ 2^D
    println("Neighborhood size is $ns^$D")
    neighborhood_origins = @view indices[ntuple(x -> 1:ns:size[x], ndims(binary_pattern))...]
    neighborhood_element_offsets = CartesianIndices(ntuple(x -> ns, ndims(binary_pattern)))

    function find_tightest_cluster()::Int
        raw_index = 1
        highest_value = energy[raw_index]
        for candidate in 2:length(energy)
            if binary_pattern[candidate] == false
                continue
            end
            value_here = energy[candidate]
            if value_here > highest_value
                highest_value = value_here
                raw_index = candidate
            end
        end
        raw_index
    end

    function find_largest_void()::Int
        raw_index = 1
        lowest_value = energy[raw_index]
        for candidate in 2:length(energy)
            if binary_pattern[candidate] == true
                continue
            end
            value_here = energy[candidate]
            if value_here < lowest_value
                lowest_value = value_here
                raw_index = candidate
            end
        end
        raw_index
    end

    function find_largest_void_in_neighborhood(neighborhood_origin::CartesianIndex{D})::CartesianIndex{D}
        pos = neighborhood_origin
        lowest_value = energy[pos]
        for offset in neighborhood_element_offsets
            candidate = CartesianIndex(Tuple(neighborhood_origin + offset) .- 1)
            if binary_pattern[candidate] == true
                continue
            end
            value_here = energy[candidate]
            if value_here < lowest_value
                lowest_value = value_here
                pos = candidate
            end
        end
        pos
    end

    # Set all pixels to be zero
    for index in indices
        binary_pattern[index] = false
    end
    # Set up to (total_pixels * ONES_DENSITY) pixels to be ones.
    for _ in 1:floor(Int64, prod(size) * parameters.starting_density^length(size))
        index = rand(indices)
        if binary_pattern[index] == true
            # Don't add energy to the same pixel twice.
            continue
        end
        binary_pattern[index] = true
        add_kernel!(energy, kernel, indices[index])
    end

    # Distribute them more regularly.
    iteration_limit = prod(size)
    while iteration_limit > 0
        tightest_cluster = find_tightest_cluster()
        binary_pattern[tightest_cluster] = false
        sub_kernel!(energy, kernel, indices[tightest_cluster])
        largest_void = find_largest_void()
        if largest_void === tightest_cluster
            # Restore the '1' we just removed.
            binary_pattern[tightest_cluster] = true
            add_kernel!(energy, kernel, indices[tightest_cluster])
            break
        end
        # @assert binary_pattern[largest_void] == false
        binary_pattern[largest_void] = true
        add_kernel!(energy, kernel, indices[largest_void])
        iteration_limit -= 1
    end
    if iteration_limit == 0
        throw(ErrorException("Iteration limit exceeded."))
    end

    initial_binary_pattern = copy(binary_pattern)
    ones_in_ibp = sum(Int64, initial_binary_pattern)
    result = zeros(Int64, size...)

    # Phase I from paper, write a value to the result for each 1 in the IBP.
    rank = ones_in_ibp - 1
    while rank >= 0
        tightest_cluster = find_tightest_cluster()
        # @assert binary_pattern[tightest_cluster] == true
        binary_pattern[tightest_cluster] = false
        sub_kernel!(energy, kernel, indices[tightest_cluster])
        result[tightest_cluster] = rank
        rank -= 1
    end

    # Phase II and III, add pixels to largest voids. The paper splits this into
    # two steps that are actually mathematically equivalent :P oh well.
    binary_pattern .= initial_binary_pattern
    # Assert all the energies are basically zero.
    # @assert prod(abs.(energy .- zeros(size)) .< 1f-4)
    rank = ones_in_ibp
    ## Set up energy to be correct.
    @inbounds for raw_index in eachindex(binary_pattern)
        if binary_pattern[raw_index] != true
            continue
        end
        index = indices[raw_index]
        add_kernel!(energy, kernel, index)
    end
    neighborhood_order = shuffle(neighborhood_origins)
    neighborhood_index = 1
    while rank < Int64(prod(size))
        neighborhood = neighborhood_order[neighborhood_index]
        largest_void = find_largest_void_in_neighborhood(neighborhood)
        # @assert binary_pattern[largest_void] == false
        binary_pattern[largest_void] = true
        add_kernel!(energy, kernel, indices[largest_void])
        result[largest_void] = rank
        rank += 1
        neighborhood_index += 1
        if neighborhood_index > length(neighborhood_order)
            neighborhood_index = 1
        end
    end

    result
end

function serialize(data::Array{Int64,D})::String where D
    buffer = Array{UInt8}(undef, 3 + length(data) * 16)
    buffer[1] = '('
    buffer[2] = 'b'
    buffer[3] = '"'
    pos = 4
    for element in data
        for digit in 1:16
            buffer[pos + 16 - digit] = b"0123456789ABCDEF"[(element >> (digit * 4 - 4)) & 0xF + 1]
        end
        pos += 16
    end
    String(buffer) * "\", $(size(data)))"
end

function make_noise(dim::Val{D}, size::Val{S}, channels::Int) where {D,S}
    dims = ntuple(x -> S, dim)
    data = Array{Int64,D + 1}(undef, (dims..., channels))
    total_time = 0.0
    for channel in 1:channels
        index = (ntuple(x -> :, dim)..., channel)
        start = time()
        data[index...] .= generate(Parameters(), dims)
        total_time += time() - start
    end
    println("Generated in $(ceil(Int64, total_time)) seconds")
    serialized = serialize(data)
    open(
        f -> write(f, serialized),
        "BlueNoise/src/generated_$(D)d$(channels)v_res$(S).jl",
        "w"
    )
end

function frequency_plot(data::Array{Float64,D}) where {D}
    heatmap(abs.(FFTW.dct(data .- 0.5)), aspect_ratio=1, clim=(0.0, 1.0))
end

export Parameters, generate, frequency_plot, make_noise

# begin
#     params = Parameters()
#     d1 = generate(params, (64, 64))
#     p1 = frequency_plot(d1)
#     params = Parameters(1//30, 1.5)
#     d2 = generate(params, (64, 64))
#     p2 = frequency_plot(d2)
#     plot(p1, p2, heatmap(d1), heatmap(d2))
# end

end # module

max_stars_repo_name: joshua-maros/julia-blue-noise

max_issues_repo_name: joshua-maros/julia-blue-noise

max_forks_repo_name: joshua-maros/julia-blue-noise


text:
module TestStaticArrays
using Test
using Setfield
using StaticArrays
using StaticNumbers

@testset "StaticArrays" begin
    obj = StaticArrays.@SMatrix [1 2; 3 4]
    @testset for l in [
            (@lens _[2,1]),
        ]
        @test get(obj, l) == 3
        @test set(obj, l, 5) == StaticArrays.@SMatrix [1 2; 5 4]
        @test setindex(obj, 5, 2, 1) == StaticArrays.@SMatrix [1 2; 5 4]
    end

    v = @SVector [1,2,3]
    @test (@set v[1] = 10) === @SVector [10,2,3]
    @test_broken (@set v[1] = π) === @SVector [π,2,3]

    @testset "Multi-dynamic indexing" begin
        two = 2
        plusone(x) = x + 1
        l1 = @lens _.a[2, 1].b
        l2 = @lens _.a[plusone(end) - two, end÷2].b
        m_orig = @SMatrix [
            (a=1, b=10) (a=2, b=20)
            (a=3, b=30) (a=4, b=40)
            (a=5, b=50) (a=6, b=60)
        ]
        m_mod = @SMatrix [
            (a=1, b=10) (a=2, b=20)
            (a=3, b=3000) (a=4, b=40)
            (a=5, b=50) (a=6, b=60)
        ]
        obj = (a=m_orig, b=4)
        @test get(obj, l1) === get(obj, l2) === 30
        @test set(obj, l1, 3000) === set(obj, l2, 3000) === (a=m_mod, b=4)
    end
end
end


text:
module NetworkAlign

using MatrixNetworks

"""
Module ``NetworkAlign``: Documentation on the module

You can check the readme file here: \n
"https://github.com/nassarhuda/NetworkAlign.jl/blob/master/README.md"
"""
NetworkAlign

include("normout.jl")
include("othermaxplus.jl")
include("othersum.jl")
include("round_messages.jl")
include("column_maxmatchsum.jl")
include("load_netalign_problem.jl")
include("netalignmr.jl")
include("isorank.jl")
include("netalignbp.jl")
include("count_overlap.jl")
include("make_squares.jl")
include("netalign_setup.jl")

export load_netalign_problem, netalignmr, isorank, netalignbp,
        count_overlap, make_squares, netalign_setup, build_Si

end # end module

max_stars_repo_name: nassarhuda/NetworkAlign

max_issues_repo_name: nassarhuda/NetworkAlign

max_forks_repo_name: nassarhuda/NetworkAlign


text:
module macro.main
using base: @kwdf
using pest
using pest.parse

macro compute_euler()
    return :(nil)
end

macro compute_algol()
    return :(nil)
end

macro compute_mcarlo()
    return :(nil)
end

@compute_euler
@compute_algol
@compute_mcarlo

# macro with one method

export macro.main
export main.compute()

max_stars_repo_name: timo-cmd/DANCE

max_issues_repo_name: timo-cmd/Twister

max_forks_repo_name: timo-cmd/DANCE


text:
using SimpleDrawing, Plots, Clines
#resize_gr_window()

C = Circle(0, 0, 15)
D = Circle(23, 0, 8)

newdraw()
draw(C, linewidth = 2);
draw(D, linewidth = 2);

x = 15 * cos(pi / 6)
y = 15 * sin(pi / 6)

draw_vector(x, y, linecolor = :black)
annotate!(x / 2, y / 2 + 1, "a")

x = 8 * cos(pi / 6)
y = 8 * sin(pi / 6)
draw_vector(x, y, 23, 0, color = :black)

annotate!(x + 19, y / 2 + 0.7, "b")

#draw_point(23+offset,0,marker=4,color=:red,linecolor=:white)

offset = 6
draw(Circle(23 + offset + 0.2, 0, 0.2), fill = true, linecolor = :red, color = :red)

draw_vector(offset, 0, 23, 0, linecolor = :red)
annotate!(26, -0.5, "offset", annotationfontsize = 10)

draw_point(0, 0, color = :black, marker = 4)
draw_point(23, 0, color = :black, marker = 4)

finish()


text:
@testset "Energy minimization" begin
    coords = [
        SVector(1.0, 1.0, 1.0)u"nm",
        SVector(1.6, 1.0, 1.0)u"nm",
        SVector(1.4, 1.6, 1.0)u"nm",
    ]
    sys = System(
        atoms=[Atom(σ=(0.4 / (2 ^ (1 / 6)))u"nm", ϵ=1.0u"kJ * mol^-1") for i in 1:3],
        pairwise_inters=(LennardJones(),),
        coords=coords,
        box_size=SVector(5.0, 5.0, 5.0)u"nm",
    )
    sim = SteepestDescentMinimizer(tol=1.0u"kJ * mol^-1 * nm^-1")

    simulate!(sys, sim)
    dists = distances(sys.coords, sys.box_size)
    dists_flat = dists[triu(trues(3, 3), 1)]
    @test all(x -> isapprox(x, 0.4u"nm"; atol=1e-3u"nm"), dists_flat)
    @test isapprox(potential_energy(sys), -3.0u"kJ * mol^-1";
                    atol=1e-4u"kJ * mol^-1")

    # No units
    coords = [
        SVector(1.0, 1.0, 1.0),
        SVector(1.6, 1.0, 1.0),
        SVector(1.4, 1.6, 1.0),
    ]
    sys = System(
        atoms=[Atom(σ=0.4 / (2 ^ (1 / 6)), ϵ=1.0) for i in 1:3],
        pairwise_inters=(LennardJones(force_units=NoUnits, energy_units=NoUnits),),
        coords=coords,
        box_size=SVector(5.0, 5.0, 5.0),
        force_units=NoUnits,
        energy_units=NoUnits,
    )
    sim = SteepestDescentMinimizer(step_size=0.01, tol=1.0)

    simulate!(sys, sim)
    dists = distances(sys.coords, sys.box_size) * u"nm"
    dists_flat = dists[triu(trues(3, 3), 1)]
    @test all(x -> isapprox(x, 0.4u"nm"; atol=1e-3u"nm"), dists_flat)
    @test isapprox(potential_energy(sys) * u"kJ * mol^-1", -3.0u"kJ * mol^-1";
                    atol=1e-4u"kJ * mol^-1")

    if run_gpu_tests
        coords = cu([
            SVector(1.0, 1.0, 1.0)u"nm",
            SVector(1.6, 1.0, 1.0)u"nm",
            SVector(1.4, 1.6, 1.0)u"nm",
        ])
        sys = System(
            atoms=cu([Atom(σ=(0.4 / (2 ^ (1 / 6)))u"nm", ϵ=1.0u"kJ * mol^-1") for i in 1:3]),
            pairwise_inters=(LennardJones(),),
            coords=coords,
            box_size=SVector(5.0, 5.0, 5.0)u"nm",
        )
        sim = SteepestDescentMinimizer(tol=1.0u"kJ * mol^-1 * nm^-1")

        simulate!(sys, sim)
        dists = distances(sys.coords, sys.box_size)
        dists_flat = dists[triu(trues(3, 3), 1)]
        @test all(x -> isapprox(x, 0.4u"nm"; atol=1e-3u"nm"), dists_flat)
        neighbors = find_neighbors(sys)
        @test isapprox(potential_energy(sys, neighbors), -3.0u"kJ * mol^-1";
                        atol=1e-4u"kJ * mol^-1")
    end
end


text:
#= functions for the power series exterior map from unit disk =#

struct PowerSeries{N,T}
  ccoeff :: Vector{T}
  dcoeff :: Vector{T}
end

function PowerSeries(ccoeff::Vector{T},dcoeff::Vector{T}) where T
  ncoeff = length(ccoeff)-2
  PowerSeries{ncoeff,T}(ccoeff,dcoeff)
end

function Base.show(io::IO,ps::PowerSeries{N,T}) where {N,T}
  println(io, "multipole coefficients:")
  println(io, "  c₁ = $(ps.ccoeff[1]), ")
  println(io, "  c₀ = $(ps.ccoeff[2]), ")
  print(io,   "  c₋ᵢ = ")
  for i = 1:N
    print(io,"$(ps.ccoeff[2+i]), ")
  end
  println(io, "i = 1:$(N)")
end

function (ps::PowerSeries)(ζ::T) where T<:Number
  ζⁿ = ζ
  z = zero(ζ)
  for c in ps.ccoeff
    z += c*ζⁿ
    ζⁿ /= ζ
  end
  return z
end

(ps::PowerSeries)(ζ::Vector{T}) where T<:Number = ps.(ζ)

struct PowerSeriesDerivatives
  ps :: PowerSeries
end

function (dps::PowerSeriesDerivatives)(ζ::T) where T<:Number
  C = dps.ps.ccoeff
  dz = C[1]
  ζⁿ = 1/ζ^2
  ddz = ComplexF64(0)
  for n in 1:length(C)-2
    dz -= n*C[n+2]*ζⁿ
    ζⁿ /= ζ
    ddz += n*(n+1)*C[n+2]*ζⁿ
  end
  return dz, ddz
end

function (dps::PowerSeriesDerivatives)(ζs::Vector{T}) where T<:Number
  dz = zero(ζs)
  ddz = zero(ζs)
  for (i,ζ) in enumerate(ζs)
    dz[i], ddz[i] = dps(ζ)
  end
  return dz, ddz
end

function evalinv_exterior(z::Vector{ComplexF64},ps::PowerSeries,
                                dps::PowerSeriesDerivatives)
#=
Evaluates the inverse of the exterior power series mapping, using a combination
of integration and Newton iteration.
=#

   ζ = zeros(ComplexF64,size(z))
   lenz = length(z)
   ζ0 = []
   maxiter = 10
   tol = 1.0e-8

   # Find z values close to vertices and set ζ to the corresponding
   # prevertices
   done = zeros(Bool,size(z))

   # Now, for remaining z values, first try to integrate
   #  dζ/dt = (z - z(ζ₀))/z'(ζ) from t = 0 to t = 1,
   # with the initial condition ζ(0) = ζ₀.
   if isempty(ζ0)
     # choose a point on the unit circle
     ζ0 = exp.(im*zeros(lenz))
     ζ0[isapprox.(angle.(z),π;atol=eps())] .= exp(im*π)
     dz0,ddz0 = dps(ζ0)
     # check for starting points on edges of the body, and rotate them
     # a bit if so
     onedge = isapprox.(abs.(dps(ζ0)[1]),0.0;atol=eps())
     ζ0[onedge] .*= exp(im*π/20)
     z0 = ps(ζ0)
   else
     z0 = ps(ζ0)
     if length(ζ0)==1 && lenz > 1
       ζ0 = repeat(transpose(ζ0), lenz)
       z0 = repeat(transpose(z0), lenz)
     end
     z0 = z0[.~done]
     ζ0 = ζ0[.~done]
   end
   odetol = max(tol,1e-3)
   scale = z[.~done] - z0

   ζ0 = [real(ζ0);imag(ζ0)]

   f(ζ,p,t) = invfunc(ζ,scale,dps)
   tspan = (0.0,1.0)
   prob = ODEProblem(f,ζ0,tspan)
   sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
   lenu = length(ζ0)
   ζ[.~done] = sol.u[end][1:lenz]+im*sol.u[end][lenz+1:lenu];
   out = abs.(ζ) .> 1
   ζ[out] = sign.(ζ[out])

   # Now use Newton iterations to improve the solution
   ζn = ζ
   k = 0
   while ~all(done) && k < maxiter
     F = z[.~done] - ps(ζn[.~done])
     M = length(F)
     dF, ddz = dps(ζn[.~done])
     ζn[.~done] = ζn[.~done] + F./dF

     done[.~done] = abs.(F).< tol
     k += 1
   end
   F = z[.~done] - ps(ζn[.~done])
   if any(abs.(F).> tol)
     error("Check solution")
   end
   ζ = ζn

end

function invfunc(u,scale,dps::PowerSeriesDerivatives)
    lenu = length(u)
    lenzp = Int(lenu/2)
    ζ = u[1:lenzp]+im*u[lenzp+1:lenu]

    dz, ddz = dps(ζ)
    f = scale./dz
    zdot = [real(f);imag(f)]
end


text:
module EOS
using PolynomialRoots
using Unitful

∑(x)=sum(x)

const R = 8.31446261815324u"MPa*cm^3/(K*mol)"

@doc raw"""
    Van_Der_Waals(T,P,Tc,Pc)

**Ecuación de la forma:**

```math
P=\frac{RT}{ν-b}-\frac{a}{ν^2}
```
Donde:
```math
\begin{aligned}
a&=\frac{27(RT_c)^2}{64P_c} & b&=\frac{RT_c}{8P_c}\\
\end{aligned}
```
```math
\frac{U-U^{ig}}{RT}=\frac{a}{RTν}
```
```math
\frac{G-G^{ig}}{RT}=z-1-\frac{a}{RTν}-\log_{10} (z-B)
```
**Los campos de entrada son:**
- `T :: Float` Es la temperatura del sistema
- `P :: Float` Es la presión del sistema
- `Tc :: Float` Es la temperatura critica del sistema
- `Pc :: Float` Es la presión critica del sistema
- `𝜔 :: Float` Es el factor acéntrico

**Salida:**

Regresa un conjunto de propiedades termodinámicas.

- `ℍ :: Float` Entalpía residual
- `𝕌 :: Float` Energía interna residual
- `𝕊 :: Float` Entropía residual
- `𝔾 :: Float` Energía libre de Gibbs residual
- `𝔸 :: Float` Energía libre de helmholtz residual
- `φ :: Float` Coeficiente de fugacidad
- `f :: Float` Fugacidad
- `z :: Float` Factor de compresibilidad
- `ν :: Float` Volumen molar

"""
function Van_Der_Waals(T,P,Tc,Pc)
    T = uconvert(u"K",T)
    P = uconvert(u"MPa",P)
    Tc = uconvert(u"K",Tc)
    Pc = uconvert(u"MPa",Pc)

    a = (27/64)*R^2 *Tc^2 /Pc
    b = R*Tc/(8*Pc)
    A = a*P/(R*T)^2
    B = b*P/(R*T)

    𝒰 = -(1+B)
    𝒮 = A
    𝒯 = -A*B

    z = roots([𝒯, 𝒮, 𝒰, 1])
    index=findall(x->abs(imag(x))>1e-5,z)
    deleteat!(z,index)
    z = real(z)
    index=findall(x->x<0,z)
    deleteat!(z,index)
    if length(z)==1
        z=z[1]
    end

    v = z*R*T/P

    U_rt = -a./(R*T*v)
    H_rt = U_rt+z.-1
    G_rt = z.-1-log.(z.-B)-a./(R*T*v)
    S_r = H_rt-G_rt
    A_rt = U_rt-S_r
    φ = exp.(G_rt)
    f = φ*P

    𝕌 = U_rt*R*T
    ℍ = H_rt*R*T
    𝕊 = S_r*R
    𝔾 = G_rt*R*T
    𝔸 = A_rt*R*T

    𝕌 = uconvert.(u"J/mol",𝕌)
    ℍ = uconvert.(u"J/mol",ℍ)
    𝕊 = uconvert.(u"J/(K*mol)",𝕊)
    𝔾 = uconvert.(u"J/mol",𝔾)
    𝔸 = uconvert.(u"J/mol",𝔸)

    ℍ = ustrip(ℍ)
    𝕌 = ustrip(𝕌)
    𝕊 = ustrip(𝕊)
    𝔾 = ustrip(𝔾)
    𝔸 = ustrip(𝔸)
    f = ustrip(f)
    v = ustrip(v)

    ℍ = round.(ℍ; digits=6)
    𝕌 = round.(𝕌; digits=6)
    𝕊 = round.(𝕊; digits=6)
    𝔾 = round.(𝔾; digits=6)
    𝔸 = round.(𝔸; digits=6)
    φ = round.(φ; digits=6)
    f = round.(f; digits=6)
    v = round.(v; digits=6)
    z = round.(z; digits=6)

    return ℍ,𝕌,𝕊,𝔾,𝔸,φ,f,z,v
end

@doc raw"""
    Redlich_Kwong(T,P,Tc,Pc)

**Ecuación de la forma:**

```math
P=\frac{RT}{ν-b}-\frac{a}{(T^{1.5} ν(ν+b))}
```
Donde:
```math
\begin{aligned}
a&=\frac{0.42748 R^2 T_c^{2.5}}{P_c} & b&=\frac{0.08664 RT_c}{P_c}\\
\end{aligned}
```
```math
\frac{U-U^{ig}}{RT}=\frac{3a\log_{10} (1+z/B)}{2bR T^{1.5}}
```
```math
\frac{G-G^{ig}}{RT}=z-1+\frac{a\log_{10} (1+z/B)}{bR T^{1.5}}-\log_{10} (z-B)
```

**Los campos de entrada son:**
- `T :: Float` Es la temperatura del sistema
- `P :: Float` Es la presión del sistema
- `Tc :: Float` Es la temperatura critica del sistema
- `Pc :: Float` Es la presión critica del sistema
- `𝜔 :: Float` Es el factor acéntrico

**Salida:**

Regresa un conjunto de propiedades termodinámicas.

- `ℍ :: Float` Entalpía residual
- `𝕌 :: Float` Energía interna residual
- `𝕊 :: Float` Entropía residual
- `𝔾 :: Float` Energía libre de Gibbs residual
- `𝔸 :: Float` Energía libre de helmholtz residual
- `φ :: Float` Coeficiente de fugacidad
- `f :: Float` Fugacidad
- `z :: Float` Factor de compresibilidad
- `v :: Float` Volumen molar
"""
function Redlich_Kwong(T,P,Tc,Pc)
    T = uconvert(u"K",T)
    P = uconvert(u"MPa",P)
    Tc = uconvert(u"K",Tc)
    Pc = uconvert(u"MPa",Pc)

    Tᵣ = T/Tc
    a = (0.42748*R^2 *Tc^2.5)/Pc
    b = 0.08664*R*Tc/Pc

    A = a*P/(R^2 *T^2.5)
    B = b*P/(R*T)
    𝒰 = -1
    𝒮 = A-B-B^2
    𝒯 = -A*B

    z = roots([𝒯, 𝒮, 𝒰, 1])
    index=findall(x->abs(imag(x))>1e-5,z)
    deleteat!(z,index)
    z = real(z)
    index=findall(x->x<0,z)
    deleteat!(z,index)
    if length(z)==1
        z=z[1]
    end

    v = z*R*T/P
    U_rt = 3*a*log.(z./(B.+z))/(2*b*R*T^1.5)
    H_rt = U_rt+z.-1
    G_rt = z.-1-log.(z.-B)+(a*log.(z./(z.+B)))./(b*R*T^1.5)
    S_r = H_rt-G_rt
    A_rt = U_rt-S_r
    φ = exp.(G_rt)
    f = φ*P
    𝕌 = U_rt*R*T
    ℍ = H_rt*R*T
    𝕊 = S_r*R
    𝔾 = G_rt*R*T
    𝔸 = A_rt*R*T

    𝕌 = uconvert.(u"J/mol",𝕌)
    ℍ = uconvert.(u"J/mol",ℍ)
    𝕊 = uconvert.(u"J/(K*mol)",𝕊)
    𝔾 = uconvert.(u"J/mol",𝔾)
    𝔸 = uconvert.(u"J/mol",𝔸)

    ℍ = ustrip(ℍ)
    𝕌 = ustrip(𝕌)
    𝕊 = ustrip(𝕊)
    𝔾 = ustrip(𝔾)
    𝔸 = ustrip(𝔸)
    f = ustrip(f)
    v = ustrip(v)

    ℍ = round.(ℍ; digits=6)
    𝕌 = round.(𝕌; digits=6)
    𝕊 = round.(𝕊; digits=6)
    𝔾 = round.(𝔾; digits=6)
    𝔸 = round.(𝔸; digits=6)
    φ = round.(φ; digits=6)
    f = round.(f; digits=6)
    v = round.(v; digits=6)
    z = round.(z; digits=6)
    #return ℍ,𝕌,𝕊,z,v,φ,f
    return ℍ,𝕌,𝕊,𝔾,𝔸,φ,f,z,v

end # function

@doc raw"""
    Soave_Redlich_Kwong(T,P,Tc,Pc,𝜔)

**Ecuación de la forma:**

```math
P=\frac{RT}{V-b}-\frac{aα}{V(V+b)}
```
Donde:
```math
\begin{aligned}
a&=\frac{0.42748 R^2 T_c^{2}}{P_c} & b&=\frac{0.08664 RT_c}{P_c} & α&=(1 + κ(1 -T_r^{0.5}))^2\\
\end{aligned}
```
```math
κ=0.48+1.574w-0.176w^2
```
```math
\frac{U-U^{ig}}{RT}=-\frac{(aα+a\sqrt{T_rα}κ)\log_{10} (1+B/z)}{bRT}
```
```math
\frac{S-S^{ig}}{R}=\log_{10}(z-B)-\frac{a\sqrt{T_rα}κ\log_{10}(1+B/z)}{bRT}
```

**Los campos de entrada son:**
- `T :: Float` Es la temperatura del sistema
- `P :: Float` Es la presión del sistema
- `Tc :: Float` Es la temperatura critica del sistema
- `Pc :: Float` Es la presión critica del sistema
- `𝜔 :: Float` Es el factor acéntrico

**Salida:**

Regresa un conjunto de propiedades termodinámicas.

- `ℍ :: Float` Entalpía residual
- `𝕌 :: Float` Energía interna residual
- `𝕊 :: Float` Entropía residual
- `𝔾 :: Float` Energía libre de Gibbs residual
- `𝔸 :: Float` Energía libre de helmholtz residual
- `φ :: Float` Coeficiente de fugacidad
- `f :: Float` Fugacidad
- `z :: Float` Factor de compresibilidad
- `v :: Float` Volumen molar
"""
function Soave_Redlich_Kwong(T,P,Tc,Pc,𝜔)
    T = uconvert(u"K",T)
    P = uconvert(u"MPa",P)
    Tc = uconvert(u"K",Tc)
    Pc = uconvert(u"MPa",Pc)

    Tᵣ = T/Tc
    κ = 0.48 + 1.574*𝜔 -0.176*𝜔^2
    α = (1 + κ*(1 -Tᵣ^.5))^2
    ac = 0.42748*((R^2)*(Tc^2)/Pc)
    a = ac*α
    b = 0.08664*R*Tc/Pc

    A = a*P/(R*T)^2
    B = b*P/(R*T)
    𝒰 = -1
    𝒮 = A-B-B^2
    𝒯 = -A*B

    z = roots([𝒯, 𝒮, 𝒰, 1])
    index=findall(x->abs(imag(x))>1e-5,z)
    deleteat!(z,index)
    z = real(z)
    index=findall(x->x<0,z)
    deleteat!(z,index)
    if length(z)==1
        z=z[1]
    end

    v = z*R*T/P
    U_rt = -((a.+ac*√Tᵣ*√α*κ)*log.(1 .+B./z))/(b*R*T)
    H_rt = U_rt+z.-1
    S_r = log.(z.-B)-(ac*√Tᵣ*√α*κ*log.(1 .+B./z))/(b*R*T)
    G_rt = H_rt-S_r
    A_rt = U_rt-S_r
    φ = exp.(G_rt)
    f = φ*P
    𝕌 = U_rt*R*T
    ℍ = H_rt*R*T
    𝕊 = S_r*R
    𝔾 = G_rt*R*T
    𝔸 = A_rt*R*T

    𝕌 = uconvert.(u"J/mol",𝕌)
    ℍ = uconvert.(u"J/mol",ℍ)
    𝕊 = uconvert.(u"J/(K*mol)",𝕊)
    𝔾 = uconvert.(u"J/mol",𝔾)
    𝔸 = uconvert.(u"J/mol",𝔸)

    ℍ = ustrip(ℍ)
    𝕌 = ustrip(𝕌)
    𝕊 = ustrip(𝕊)
    𝔾 = ustrip(𝔾)
    𝔸 = ustrip(𝔸)
    f = ustrip(f)
    v = ustrip(v)

    ℍ = round.(ℍ; digits=6)
    𝕌 = round.(𝕌; digits=6)
    𝕊 = round.(𝕊; digits=6)
    𝔾 = round.(𝔾; digits=6)
    𝔸 = round.(𝔸; digits=6)
    φ = round.(φ; digits=6)
    f = round.(f; digits=6)
    v = round.(v; digits=6)
    z = round.(z; digits=6)
    #return ℍ,𝕌,𝕊,z,v,φ,f
    return ℍ,𝕌,𝕊,𝔾,𝔸,φ,f,z,v

end # function

@doc raw"""
    Peng_Robinson(T,P,Tc,Pc,𝜔)

**Ecuación de la forma:**

```math
P=\frac{RT}{V-b}-\frac{aα}{V^2 +2bV+b^2}
```
Donde:
```math
\begin{aligned}
a&=\frac{0.45724 R^2 T_c^{2}}{P_c} & b&=\frac{0.0778 RT_c}{P_c} & α&=(1 + κ(1 -T_r^{0.5}))^2\\
\end{aligned}
```
```math
κ=0.37464 + 1.54226*w -0.26992*w^2
```
```math
\frac{U-U^{ig}}{RT}=-(aα+a\sqrt{T_rα}κ)*(\coth^{-1}{\sqrt{2}}+\tanh^{-1}\frac{(-1+B/z)/\sqrt{2}}{bRT\sqrt{2}})
```
```math
\frac{G-G^{ig}}{R}=z-1-\log(z-B)-aα*(\coth^{-1}{\sqrt{2}}+\tanh^{-1}\frac{(B-z)/z\sqrt{2}}{bRT\sqrt{2}})
```

**Los campos de entrada son:**
- `T :: Float` Es la temperatura del sistema
- `P :: Float` Es la presión del sistema
- `Tc :: Float` Es la temperatura critica del sistema
- `Pc :: Float` Es la presión critica del sistema
- `𝜔 :: Float` Es el factor acéntrico

**Salida:**

Regresa un conjunto de propiedades termodinámicas.

- `ℍ :: Float` Entalpía residual
- `𝕌 :: Float` Energía interna residual
- `𝕊 :: Float` Entropía residual
- `𝔾 :: Float` Energía libre de Gibbs residual
- `𝔸 :: Float` Energía libre de helmholtz residual
- `φ :: Float` Coeficiente de fugacidad
- `f :: Float` Fugacidad
- `z :: Float` Factor de compresibilidad
- `v :: Float` Volumen molar
"""
function Peng_Robinson(T,P,Tc,Pc,𝜔)
    T = uconvert(u"K",T)
    P = uconvert(u"MPa",P)
    Tc = uconvert(u"K",Tc)
    Pc = uconvert(u"MPa",Pc)

    Tᵣ = T/Tc
    κ = 0.37464 + 1.54226*𝜔 -0.26992*𝜔^2
    α = (1 + κ*(1 -Tᵣ^.5))^2
    ac = 0.45724*((R^2)*(Tc^2)/Pc)
    a = ac*α
    b = 0.0778*R*Tc/Pc

    𝒰 = b*P/(R*T) -1
    𝒮 = a*P/(R*T)^2 -3*(b*P/(R*T))^2 -2*b*P/(R*T)
    𝒯 = (b*P/(R*T))^3 +(b*P/(R*T))^2 -a*P*b*P/(R*T)^3

    z = roots([𝒯, 𝒮, 𝒰, 1])
    index=findall(x->abs(imag(x))>1e-5,z)
    deleteat!(z,index)
    z = real(z)
    index=findall(x->x<0,z)
    deleteat!(z,index)
    if length(z)==1
        z=z[1]
    end

    v = z*R*T/P
    A1 = a*P/(R*T)^2
    B = b*P/(R*T)

    U_rt = -(a.+ac*√Tᵣ*√α*κ)*(acoth(√2).+atanh.((-1 .+B./z)./√2))/(√2*b*R*T)
    H_rt = U_rt+z.-1
    G_rt = z.-1-log.(z.-B)-a*(acoth(√2).+atanh.((B.-z)./(z*√2)))/(√2*b*R*T)
    S_r = H_rt-G_rt
    A_rt = U_rt-S_r
    φ = exp.(G_rt)
    f = φ*P
    𝕌 = U_rt*R*T
    ℍ = H_rt*R*T
    𝕊 = S_r*R
    𝔾 = G_rt*R*T
    𝔸 = A_rt*R*T

    𝕌 = uconvert.(u"J/mol",𝕌)
    ℍ = uconvert.(u"J/mol",ℍ)
    𝕊 = uconvert.(u"J/(K*mol)",𝕊)
    𝔾 = uconvert.(u"J/mol",𝔾)
    𝔸 = uconvert.(u"J/mol",𝔸)

    ℍ = ustrip(ℍ)
    𝕌 = ustrip(𝕌)
    𝕊 = ustrip(𝕊)
    𝔾 = ustrip(𝔾)
    𝔸 = ustrip(𝔸)
    f = ustrip(f)
    v = ustrip(v)

    ℍ = round.(ℍ; digits=6)
    𝕌 = round.(𝕌; digits=6)
    𝕊 = round.(𝕊; digits=6)
    𝔾 = round.(𝔾; digits=6)
    𝔸 = round.(𝔸; digits=6)
    φ = round.(φ; digits=6)
    f = round.(f; digits=6)
    v = round.(v; digits=6)
    z = round.(z; digits=6)

    return ℍ,𝕌,𝕊,𝔾,𝔸,φ,f,z,v
end

function Van_Der_Waals_M(T,P,Tc,Pc,y)
    n = length(y)
    T = uconvert(u"K",T)
    P = uconvert(u"MPa",P)
    Tc = uconvert.(u"K",Tc)
    Pc = uconvert.(u"MPa",Pc)

    a = (27/64)*R^2 *Tc.^2 ./Pc
    b = R*Tc./(8*Pc)
    aa(i,j) =√(a[i]*a[j])
    aₘ = ∑(∑(y[i]*y[j]*aa(i,j) for j=1:n) for i=1:n)
    bₘ = ∑(y[i]*b[i] for i=1:n)

    A = aₘ*P/(R*T)^2
    B = bₘ*P/(R*T)

    𝒰 = -(1+B)
    𝒮 = A
    𝒯 = -A*B

    z = roots([𝒯, 𝒮, 𝒰, 1])
    index=findall(x->abs(imag(x))>1e-5,z)
    deleteat!(z,index)
    z = real(z)
    index=findall(x->x<0,z)
    deleteat!(z,index)
    if length(z)==1
        z=z[1]
    end

    v = z*R*T/P

    U_rt = -aₘ./(R*T*v)
    H_rt = U_rt+z.-1
    G_rt = z.-1-log.(1 .-bₘ./v)-log.(z)-aₘ./(R*T*v)
    S_r = H_rt-G_rt
    A_rt = U_rt-S_r
    φ = exp.(G_rt)
    f = φ*P

    𝕌 = U_rt*R*T
    ℍ = H_rt*R*T
    𝕊 = S_r*R
    𝔾 = G_rt*R*T
    𝔸 = A_rt*R*T

    𝕌 = uconvert.(u"J/mol",𝕌)
    ℍ = uconvert.(u"J/mol",ℍ)
    𝕊 = uconvert.(u"J/(K*mol)",𝕊)
    𝔾 = uconvert.(u"J/mol",𝔾)
    𝔸 = uconvert.(u"J/mol",𝔸)

    ℍ = ustrip(ℍ)
    𝕌 = ustrip(𝕌)
    𝕊 = ustrip(𝕊)
    𝔾 = ustrip(𝔾)
    𝔸 = ustrip(𝔸)
    f = ustrip(f)
    v = ustrip(v)

    ℍ = round.(ℍ; digits=6)
    𝕌 = round.(𝕌; digits=6)
    𝕊 = round.(𝕊; digits=6)
    𝔾 = round.(𝔾; digits=6)
    𝔸 = round.(𝔸; digits=6)
    φ = round.(φ; digits=6)
    f = round.(f; digits=6)
    v = round.(v; digits=6)
    z = round.(z; digits=6)

    return ℍ,𝕌,𝕊,𝔾,𝔸,φ,f,z,v
end

function Redlich_Kwong_M(T,P,Tc,Pc,y)
    n = length(y)
    T = uconvert(u"K",T)
    P = uconvert(u"MPa",P)
    Tc = uconvert.(u"K",Tc)
    Pc = uconvert.(u"MPa",Pc)

    Tᵣ = T./Tc
    a = (0.42748*R^2 *Tc.^2.5)./Pc
    b = 0.08664*R*Tc./Pc

    aa(i,j) =√(a[i]*a[j])
    #aa(i,j) =(1-k)*√(a[i]*a[j])
    aₘ = ∑(∑(y[i]*y[j]*aa(i,j) for j=1:n) for i=1:n)
    bₘ = ∑(y[i]*b[i] for i=1:n)

    A = aₘ*P/(R^2 *T^2.5)
    B = bₘ*P/(R*T)
    𝒰 = -1
    𝒮 = A-B-B^2
    𝒯 = -A*B

    z = roots([𝒯, 𝒮, 𝒰, 1])
    index=findall(x->abs(imag(x))>1e-5,z)
    deleteat!(z,index)
    z = real(z)
    index=findall(x->x<0,z)
    deleteat!(z,index)
    if length(z)==1
        z=z[1]
    end

    v = z*R*T/P

    U_rt = 3*aₘ*log.(z./(B.+z))/(2*bₘ*R*T^1.5)
    H_rt = U_rt+z.-1
    G_rt = z.-1-log.(z.-B)+(aₘ*log.(z./(z.+B)))./(bₘ*R*T^1.5)
    S_r = H_rt-G_rt
    A_rt = U_rt-S_r
    φ = exp.(G_rt)
    f = φ*P
    𝕌 = U_rt*R*T
    ℍ = H_rt*R*T
    𝕊 = S_r*R
    𝔾 = G_rt*R*T
    𝔸 = A_rt*R*T

    𝕌 = uconvert.(u"J/mol",𝕌)
    ℍ = uconvert.(u"J/mol",ℍ)
    𝕊 = uconvert.(u"J/(K*mol)",𝕊)
    𝔾 = uconvert.(u"J/mol",𝔾)
    𝔸 = uconvert.(u"J/mol",𝔸)

    ℍ = ustrip(ℍ)
    𝕌 = ustrip(𝕌)
    𝕊 = ustrip(𝕊)
    𝔾 = ustrip(𝔾)
    𝔸 = ustrip(𝔸)
    f = ustrip(f)
    v = ustrip(v)

    ℍ = round.(ℍ; digits=6)
    𝕌 = round.(𝕌; digits=6)
    𝕊 = round.(𝕊; digits=6)
    𝔾 = round.(𝔾; digits=6)
    𝔸 = round.(𝔸; digits=6)
    φ = round.(φ; digits=6)
    f = round.(f; digits=6)
    v = round.(v; digits=6)
    z = round.(z; digits=6)
    #return ℍ,𝕌,𝕊,z,v,φ,f
    return ℍ,𝕌,𝕊,𝔾,𝔸,φ,f,z,v

end # function

function Soave_Redlich_Kwong_M(T,P,Tc,Pc,y,𝜔)
    n = length(y)
    T = uconvert(u"K",T)
    P = uconvert(u"MPa",P)
    Tc = uconvert.(u"K",Tc)
    Pc = uconvert.(u"MPa",Pc)

    Tᵣ = T./Tc
    κ = 0.48508 .+ 1.55171*𝜔 -0.15613*𝜔.^2
    α = (1 .+ κ.*(1 .-Tᵣ.^.5)).^2
    ac = 0.42747*((R^2)*(Tc.^2)./Pc)

    a = ac.*α
    b = 0.08664*R*Tc./Pc
    A = a*P/(R*T)^2
    B = b*P/(R*T)

    aa(i,j) =√(a[i]*a[j])
    AA(i,j) =√(A[i]*A[j])
    #aa(i,j) =(1-k)*√(a[i]*a[j])
    aₘ = ∑(∑(y[i]*y[j]*aa(i,j) for j=1:n) for i=1:n)
    bₘ = ∑(y[i]*b[i] for i=1:n)
    Aₘ = aₘ*P/(R*T)^2
    Bₘ = bₘ*P/(R*T)

    𝒰 = -1
    𝒮 = Aₘ-Bₘ-Bₘ^2
    𝒯 = -Aₘ*Bₘ

    z = roots([𝒯, 𝒮, 𝒰, 1])
    index=findall(x->abs(imag(x))>1e-5,z)
    deleteat!(z,index)
    z = real(z)
    index=findall(x->x<0,z)
    deleteat!(z,index)
    if length(z)==1
        z=z[1]
    end

    D = (-1/sqrt(T))*∑(∑(y[i]*y[j]*aa(i,j)*(κ[j]/sqrt(Tc[j]*α[j])+κ[i]/sqrt(Tc[i]*α[i]))/2 for j=1:n) for i=1:n)

    v = z*R*T/P

    U_rt = log(ustrip(R*T))*(aₘ-T*D)/(bₘ*R*T) .-log.(ustrip(R*T.+Bₘ*R*T./z))*(aₘ-T*D)/(bₘ*R*T)
    H_rt = z.-1 + U_rt
    A_rt = -log.(z) +(-R*T*log.(1 .-Bₘ./z)-aₘ*log.(1 .+Bₘ./z)/bₘ)/(R*T)
    S_r = -A_rt +U_rt
    G_rt = H_rt-S_r

    φ = exp.(G_rt)
    f = φ*P
    𝕌 = U_rt*R*T
    ℍ = H_rt*R*T
    𝕊 = S_r*R
    𝔾 = G_rt*R*T
    𝔸 = A_rt*R*T

    𝕌 = uconvert.(u"J/mol",𝕌)
    ℍ = uconvert.(u"J/mol",ℍ)
    𝕊 = uconvert.(u"J/(K*mol)",𝕊)
    𝔾 = uconvert.(u"J/mol",𝔾)
    𝔸 = uconvert.(u"J/mol",𝔸)

    ℍ = ustrip(ℍ)
    𝕌 = ustrip(𝕌)
    𝕊 = ustrip(𝕊)
    𝔾 = ustrip(𝔾)
    𝔸 = ustrip(𝔸)
    f = ustrip(f)
    v = ustrip(v)

    ℍ = round.(ℍ; digits=6)
    𝕌 = round.(𝕌; digits=6)
    𝕊 = round.(𝕊; digits=6)
    𝔾 = round.(𝔾; digits=6)
    𝔸 = round.(𝔸; digits=6)
    φ = round.(φ; digits=6)
    f = round.(f; digits=6)
    v = round.(v; digits=6)
    z = round.(z; digits=6)
    #return ℍ,𝕌,𝕊,z,v,φ,f
    return ℍ,𝕌,𝕊,𝔾,𝔸,φ,f,z,v

end # function

"""
    Peng_Robinson_M(T,P,Tc,Pc,y,𝜔)

**Ecuación de la forma:**

```math
P=\\frac{RT}{V-b}-\\frac{aα}{V^2 +2bV+b^2}
```

**Los campos de entrada son:**
- `T :: Float` Es la temperatura del sistema
- `P :: Float` Es la presión del sistema
- `Tc :: Float` Es la temperatura critica del sistema
- `Pc :: Float` Es la presión critica del sistema
- `𝜔 :: Float` Es el factor acéntrico

**Salida:**

Regresa un conjunto de propiedades termodinámicas.

- `ℍ :: Float` Entalpía residual
- `𝕌 :: Float` Energía interna residual
- `𝕊 :: Float` Entropía residual
- `𝔾 :: Float` Energía libre de Gibbs residual
- `𝔸 :: Float` Energía libre de helmholtz residual
- `φ :: Float` Coeficiente de fugacidad
- `f :: Float` Fugacidad
- `z :: Float` Factor de compresibilidad
- `v :: Float` Volumen molar
"""
function Peng_Robinson_M(T,P,Tc,Pc,y,𝜔)
    n = length(y)
    T = uconvert(u"K",T)
    P = uconvert(u"MPa",P)
    Tc = uconvert.(u"K",Tc)
    Pc = uconvert.(u"MPa",Pc)

    Tᵣ = T./Tc
    κ = 0.37464 .+ 1.54226*𝜔 -0.26993*𝜔.^2
    α = (1 .+ κ.*(1 .-Tᵣ.^.5)).^2
    ac = 0.457235529*((R^2)*(Tc.^2)./Pc)

    a = ac.*α
    b = 0.077796074*R*Tc./Pc
    A = a*P/(R*T)^2
    B = b*P/(R*T)

    aa(i,j) =√(a[i]*a[j])
    AA(i,j) =√(A[i]*A[j])
    #aa(i,j) =(1-k)*√(a[i]*a[j])
    aₘ = ∑(∑(y[i]*y[j]*aa(i,j) for j=1:n) for i=1:n)
    bₘ = ∑(y[i]*b[i] for i=1:n)
    Aₘ = aₘ*P/(R*T)^2
    Bₘ = bₘ*P/(R*T)

    𝒰 = +Bₘ-1
    𝒮 = Aₘ-2*Bₘ-3*Bₘ^2
    𝒯 = -(Aₘ*Bₘ-Bₘ^2-Bₘ^3)

    z = roots([𝒯, 𝒮, 𝒰, 1])
    index=findall(x->abs(imag(x))>1e-5,z)
    deleteat!(z,index)
    z = real(z)
    index=findall(x->x<0,z)
    deleteat!(z,index)
    if length(z)==1
        z=z[1]
    end

    v = z*R*T/P

    D = (-1/√T)*∑(∑(y[i]*y[j]*aa(i,j)*(κ[j]/√(Tc[j]*α[j])+κ[i]/√(Tc[i]*α[i]))/2 for j=1:n) for i=1:n)
    H_rt = -1 .+z+(acoth(√2).+atanh.((-1 .+Bₘ./z)/(√2)))*(-aₘ+T*D)/(√2*bₘ*R*T)
    U_rt = H_rt-(z.-1)
    A_rt = -log.(z.-Bₘ)-aₘ*(acoth(√2).+atanh.((-1 .+Bₘ./z)/(√2)))/(√2*bₘ*R*T)
    S_r = U_rt-A_rt
    G_rt = H_rt-S_r

    φ = exp.(G_rt)
    f = φ*P
    𝕌 = U_rt*R*T
    ℍ = H_rt*R*T
    𝕊 = S_r*R
    𝔾 = G_rt*R*T
    𝔸 = A_rt*R*T

    𝕌 = uconvert.(u"J/mol",𝕌)
    ℍ = uconvert.(u"J/mol",ℍ)
    𝕊 = uconvert.(u"J/(K*mol)",𝕊)
    𝔾 = uconvert.(u"J/mol",𝔾)
    𝔸 = uconvert.(u"J/mol",𝔸)

    ℍ = ustrip(ℍ)
    𝕌 = ustrip(𝕌)
    𝕊 = ustrip(𝕊)
    𝔾 = ustrip(𝔾)
    𝔸 = ustrip(𝔸)
    f = ustrip(f)
    v = ustrip(v)

    ℍ = round.(ℍ; digits=6)
    𝕌 = round.(𝕌; digits=6)
    𝕊 = round.(𝕊; digits=6)
    𝔾 = round.(𝔾; digits=6)
    𝔸 = round.(𝔸; digits=6)
    φ = round.(φ; digits=6)
    f = round.(f; digits=6)
    v = round.(v; digits=6)
    z = round.(z; digits=6)
    return ℍ,𝕌,𝕊,𝔾,𝔸,φ,f,z,v
end

end

max_stars_repo_name: EmilioAlvizo/TermoQuimica

max_issues_repo_name: EmilioAlvizo/TermoQuimica

max_forks_repo_name: EmilioAlvizo/TermoQuimica


text:
function _not_EAGO_error!(m::JuMP.Model)
    if JuMP.solver_name(m) !== "EAGO: Easy Advanced Global Optimization"
        error("Solve attached to model must be EAGO.Optimizer")
    end
end

#=
Reference for auxillary variables
=#
struct AuxiliaryVariableRef <: JuMP.AbstractVariableRef
    idx::Int
    model::JuMP.Model
end

Base.@kwdef mutable struct _AuxVarData
    aux::Dict{JuMP.VariableRef, AuxiliaryVariableRef} = Dict{JuMP.VariableRef, AuxiliaryVariableRef}()
    #mimo_expr::Vector{MIMOExpr} = MIMOExpr[]
    last_hook::Union{Nothing,Function} = nothing
end
function is_auxilliary_variable(m::_AuxVarData, i::Int)
    false
end
is_auxilliary_variable(::Nothing, ::Int) = false

#=
function aux_variable_optimizehook(model::JuMP.Model)
    initialize_auxillary_variables!(model)
    model.optimize_hook = model.ext[:aux_var].last_hook
    optimize!(model)
    model.optimize_hook = aux_variable_optimizehook
    return
end
function _initialize_auxillary_variable_data(model::Model)
    model.ext[:aux_var] = _AuxVarData()
    model.ext[:aux_var].last_hook = model.optimize_hook
    model.optimize_hook = aux_variable_optimizehook
    return model
end
function enable_auxillary_variables(model::Model)
    haskey(model.ext, :aux_var) && error("Model has auxillary variables parameter enabled")
    return _initialize_auxillary_variable_data(model)
end
EAGOModel(m::JuMP.Model) = enable_auxillary_variables(m)
EAGOModel() = EAGOModel(Model(EAGO.Optimizer))

_getmodel(v::AuxiliaryVariableRef) = v.model
function _getauxdata(v::AuxiliaryVariableRef)::_AuxVarData
    return _getauxdata(_getmodel(v))::_AuxVarData
end
function _getauxdata(model::Model)::_AuxVarData
    auxvar = get(model.ext, :aux_var, nothing)
    if auxvar !== nothing
        return auxvar
    end
    return enable_auxillary_variables(model)
end

JuMP.index(p::AuxiliaryVariableRef) = p.idx
Base.iszero(::AuxiliaryVariableRef) = false
Base.copy(p::AuxiliaryVariableRef) = ParameterRef(p.idx, p.model)

struct AuxiliaryVariableNotOwned <: Exception
    aux::AuxiliaryVariableRef
end
JuMP.owner_model(p::AuxiliaryVariableRef) = p.model
function JuMP.check_belongs_to_model(p::AuxiliaryVariableRef, model::AbstractModel)
    if owner_model(p) !== model
        throw(AuxiliaryVariableNotOwned(p))
    end
end
function JuMP.is_valid(model::Model, aux::AuxiliaryVariableRef)
    return model === owner_model(aux)
end

function Base.hash(p::AuxiliaryVariableRef, h::UInt)
    return hash(objectid(owner_model(p)), hash(p.ind, h))
end
function Base.isequal(p1::AuxiliaryVariableRef, p2::AuxiliaryVariableRef)
    return owner_model(p1) === owner_model(p2) && p1.ind == p2.ind
end

JuMP.name(p::AuxiliaryVariableRef) = get(_getauxdata(p).names, p, "")
JuMP.set_name(p::AuxiliaryVariableRef, s::String) = _getauxdata(p).names[p] = s

struct AuxVar end

_aux_msg(msg) = "Invalid initialization of auxillary variable. " * msg * " not supported."
function JuMP.build_variable(_error::Function, info::JuMP.VariableInfo, ::AuxVar)
    info.has_lb    && _error(_aux_msg("Lower bound"))
    info.has_ub    && _error(_aux_msg("Upper bound"))
    info.binary    && _error(_aux_msg("Binary"))
    info.integer   && _error(_aux_msg("Integer"))
    info.has_start && _error(_aux_msg("Initial value"))
    info.has_fix   && _error(_aux_msg("Fixed value"))
    return AuxVar()
end
function JuMP.add_variable(m::JuMP.Model, v::AuxVar, name::String="")
    vref = _add_auxillary_variable(m)
    if !isempty(name)
        JuMP.set_name(vref, name)
    end
    return vref
end

macro auxiliary_variable(m, args...)
    esc(quote
            @show ($args)
            @show ($args...)
            @show $(args...)
            #@variable($m, ($args...), AuxVar())
        end
    )
end

macro mimo_expression(m, name, f!, y, x) end

macro mimo_expression(m, f!, y, x)
    esc(quote
            name = Symbol($f!)
            @mimo_expression($m, name, $f!, $y, $x)
        end
    )
end
=#

#=
struct MIMOExpr
    m::Model
    f!::Function
    y::Vector{Union{AuxiliaryVariableRef, JuMP.VariableRef}}
    x::Vector{Union{AuxiliaryVariableRef, JuMP.VariableRef}}
end

function initialize_auxillary_variables!(m::JuMP.Model)
    set_optimizer_attribute(m, "_auxillary_variable_info", m.ext[:aux_var])
    return
end

function aux_variable_optimize!(m::JuMP.Model)
    initialize_auxillary_variables!(m)
    m.optimize_hook = m.ext[:aux_var].last_hook
    optimize!(m)
    m.optimize_hook = aux_variable_optimize!
    return
end

function JuMP.add_variable(m::JuMP.Model, v::AuxVar, name::String="")
    x = JuMP.add_variable(m, JuMP.ScalarVariable(v.info), name)
    m.ext[:aux_var].var_to_aux[x] = v
    return x
end

JuMP.build_variable(_err::Function, info::JuMP.VariableInfo, ::Type{AuxVar}; kwargs...) = AuxVar(info)
function JuMP.add_variable(m::JuMP.Model, v::AuxVar, name::String)
    x = JuMP.add_variable(m, JuMP.ScalarVariable(v.info), name)
    m.ext[:aux_var].var_to_aux[x] = v
    return x
end

"
m = EAGOModel()
@variable(m, -2 <= x[i=1:2] <= 2)
"
function add_auxillary_variable(m::JuMP.Model, x, l, u)
end

"""
"""
function add_mimo_expression(m::JuMP.Model, name::String, f!::Function, y, x)
    MIMOExpr
    m = mimo_expr
end
=#

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
const single_qubit_ops = (
    (:X, :sX), (:Y, :sY), (:Z, :sZ), (:H, :sHadamard),
    (:S, :sPhase), (:Sdag, :sInvPhase),
    (:I2, :sId1)
)
for (sym, op) in single_qubit_ops
    @eval function Yao.instruct!(r::StabilizerReg, ::Val{$(QuoteNode(sym))}, locs)
        for l in locs
            QuantumClifford.apply!(r.st, QuantumClifford.$(op)(l))
        end
        return r
    end
end

function Yao.instruct!(r::StabilizerReg, ::Val{:SWAP}, locs)
    length(locs) == 2 || error("SWAP gate should be applied on 2 qubits")
    QuantumClifford.apply!(r.st, QuantumClifford.sSWAP(locs[1], locs[2]))
    return r
end

function Yao.instruct!(r::StabilizerReg, ::Val{:X}, locs, ctrl_locs, ctrl_configs)
    length(ctrl_locs) == 1 || error("Gates with multiple control qubits are not supported")
    for l in locs
        ctrl_configs !== (1,) && QuantumClifford.apply!(r.st, QuantumClifford.sX(ctrl_locs))
        QuantumClifford.apply!(r.st, QuantumClifford.sCNOT(ctrl_locs[1], l))
        ctrl_configs !== (1,) && QuantumClifford.apply!(r.st, QuantumClifford.sX(ctrl_locs))
    end
    return r
end
function Yao.instruct!(r::StabilizerReg, ::Val{:Z}, locs, ctrl_locs, ctrl_configs)
    Yao.instruct!(r, Val(:H), locs)
    Yao.instruct!(r, Val(:X), locs, ctrl_locs, ctrl_configs)
    Yao.instruct!(r, Val(:H), locs)
    return r
end