text:
include("../src/autodiff_module.jl")
include("../../__lib__/math/common/numder/src/numder_module.jl")

module dtest
using PyPlot
PyPlot.pygui(true)
using ..autodiff
using ..numder

x = autodiff.DualNumber(1.0)
y = autodiff.DualNumber(3.0)

z = x + y

fcn1(x) = 1/3 * x^3 + 1/2 * x^2 + 5 * x + 10.0
fcn2_a(x) = x[1] + x[2]
fcn2_b(x) = x[1] * x[2]
fcn2_c(x) = x[1] / x[2]
fcn3(x) = [fcn2_a(x), fcn2_b(x), fcn2_c(x)]

x0 = 2.0
d1 = autodiff.derivative(fcn1, x0)

mode = autodiff.ForwardMode()

x0_a = [1.0, 3]
g2_a = autodiff.gradient(mode, fcn2_a, x0_a)
g2_b = autodiff.gradient(mode, fcn2_b, x0_a)
g2_c = autodiff.gradient(mode, fcn2_c, x0_a)

jac = autodiff.jacobian(mode, fcn3, x0_a)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
function initial_auth()
    @info "open the url in browser and copy the apikey back"
    auth_url = HTTP.URI(
        host="auth.tdameritrade.com",
        scheme="https",
        path="/auth",
        query = [
                 "response_type" => "code",
                 "redirect_uri" => AUTH_KEY.CALLBACK_URI,
                 "client_id"=>"$(AUTH_KEY.CONSUMER_KEY)@AMER.OAUTHAP"
                ]
   )
end

function access_token(code)
    url = "https://api.tdameritrade.com/v1/oauth2/token"
    params = HTTP.escape([
              "grant_type" => "authorization_code",
              "refresh_token" => "",
              "access_type" => "offline",
              "code" => code,
              "client_id" => AUTH_KEY.CONSUMER_KEY,
              "redirect_uri" => AUTH_KEY.CALLBACK_URI
             ])
    json_result = HTTP.post(
                            url,
                ["Content-Type"=>"application/x-www-form-urlencoded"],
                params
               ).body |> String |> JSON3.read
    json_result
end

function refresh(refresh_token)
    if now() - AUTH_KEY.LAST_REFRESH < Minute(29)
        @info "Refresh not needed"
        return AUTH_KEY.ACCESS_TOKEN
    end
    url = "https://api.tdameritrade.com/v1/oauth2/token"
    params = HTTP.escape([
              "grant_type" => "refresh_token",
              "refresh_token" => refresh_token,
              "access_type" => "",
              "code" => "",
              "client_id" => AUTH_KEY.CONSUMER_KEY,
              "redirect_uri" => ""
             ])
    json_result = HTTP.post(
                            url,
                ["Content-Type"=>"application/x-www-form-urlencoded"],
                params
               ).body |> String |> JSON3.read
    AUTH_KEY.LAST_REFRESH = now()
    json_result[:access_token]
end

"""
TD_auth()

run authentication procedure
1. already authenticated -> do nothing
2. access token expired && has refresh token -> refresh access token
3. refresh token expired -> get refresh token
4. no CONSUMER_KEY -> error
"""
function TD_auth()
    cache_path = joinpath(homedir(), ".JL_TD_TOKENS_CACHE")
    if isfile(cache_path)
        last_refresh, last_date = readdlm(cache_path, ',')
        if Date(last_date) > today() - Day(89) # refresh token expires every 90 days
            AUTH_KEY.REFRESH_TOKEN = last_refresh
        else
            @info "Cached token is too old, getting new one"
        end
    end

    AUTH_KEY.CONSUMER_KEY = try
        ENV["JL_TD_CONSUMER_KEY"]
    catch
        error("Key environmen variable JL_TD_CONSUMER_KEY not found")
    end

    if AUTH_KEY.REFRESH_TOKEN == ""
        @info "First time auth needs manual invervention,"
        println(initial_auth())
        @info "extract code=<copy this> from address bar to here"
        print("code:")
        AUTH_KEY.CODE = readline() |> HTTP.URIs.unescapeuri
        token = access_token(AUTH_KEY.CODE)
        AUTH_KEY.ACCESS_TOKEN, AUTH_KEY.REFRESH_TOKEN = token[:access_token], token[:refresh_token]

        AUTH_KEY.LAST_REFRESH = now()
        @info "creating cache at $cache_path for tokens"
        open(cache_path, "w") do io
            writedlm(io, [AUTH_KEY.REFRESH_TOKEN today()], ',')
        end
    else
        @info "REFRESH_TOKEN found, refreshing ACCESS_TOKEN"
        AUTH_KEY.ACCESS_TOKEN = refresh(AUTH_KEY.REFRESH_TOKEN)
    end

    @info "Authentication completed."
    nothing
end

ext: jl


text:
"""
    integerfeasibility(A,d)

Given a linear system `Ax=d`, return an integer vector `x` which satisfies the
system.

    integerfeasibility(A,d,true)

If the third argument is present and is `true`, then as well as returning a
solution `x`, also return a matrix `xNull` of vectors in the null space of
`A` which could be added to the `x` vector to find a solution which
satisfies a constraint such as `0 .≤ x .≤ u`; see the paper below.

This is not a robust tool, just a demo.

"Solving A System Of Diophantine Equations With Bounds On The Variables" by
Karen Aardal, Cor Hurkens, and Arjen Lenstra in Integer Programming and
Combinatorial Optimization, 6th International IPCO Conference, vol 1412, pp
229-242, 1998. See http://softlib.rice.edu/pub/CRPC-TRs/reports/CRPC-TR98782.pdf

# Examples
```jldoctest
julia> A=[10 1 -9; 1 8 8]; xtrue=[0; 2; 9]; d=A*xtrue;

julia> integerfeasibility(A,d)
3-element Vector{Int64}:
 0
 2
 9

julia> A=[10 1.1 -9.1; 1 8 8]; d=A*xtrue;

julia> integerfeasibility(A,d)
3-element Vector{Float64}:
 0.0
 2.0
 9.0

julia> n=20;m=30; A = rand(-10:10,n,m); xtrue = rand(0:10,m); d=A*xtrue;

julia> sum(abs.(xtrue - integerfeasibility(A,d) ))
0

```
"""
function integerfeasibility(A::AbstractArray{Td,2},
                            d::AbstractArray{Td,1}, flag=false) where {Td<:Number}
    m,n=size(A)
    # N1=10^3
    # N2=10^4
    N1=10^4
    N2=10^5
    B = [[I           zeros(Td,n,1)];
         zeros(Td,1,n) N1;
         N2*A          -N2*d]

    Bhat,_ = lll(B)
    #Bhat = l2(B)  # doesn't work?

    # the columns are not always sorted. A sort here helps sometimes.

    Bprime = Bhat[1:n+1,1:n-m+1]
    x_d = Bprime[1:n,end]

    if flag
        xNull=Bprime[1:n,1:n-m]
        return x_d,xNull
    else
        return x_d
    end
end

"""
    x = subsetsum(a,s)

For a vector of integers `a`, and an integer `s`, try to find a binary
vector `x` such that `x'*a=s`. We use the LLL algorithm to find the
solution. This is not a robust tool, just a demo.

This function tries first the technique in the `lagariasodlyzko` function,
and if it fails, a solution via `mdsubsetsum` is attempted.

It appears that this function can also solve some integer relations
problems. See the first example.

# Examples
```jldoctest
julia> a=[1.5;.5;0;.1;.2]; s=2.2; x,_=subsetsum(a,s,true); s-x'*a
A binary Lagarias-Odlyzko solution was found.
A solution was found via lagariasodlyzko
0.0

julia> a=[32771,65543,131101,262187,524387,1048759, # from Bremner p 117
          2097523,4195057,8390143,16780259,33560539,
          67121039,134242091,268484171,536968403];

julia> s=891221976; x,_=subsetsum(a,s,false); s-x'*a
0.0

julia> N=40;a=rand(1:2^BigInt(256),N);xtrue=rand(Bool,N); s=a'*xtrue;

julia> setprecision(BigFloat,300); x,_=subsetsum(a,s,false); s-x'*a
0.0

```
"""
function subsetsum(a::AbstractArray{Td,1},ss::Td,verbose=false) where {Td<:Number}

    xlo,binarylo = lagariasodlyzko(a,ss,0,verbose);
    if binarylo===true && ~any(ismissing.(xlo))
        verbose && print("A solution was found via lagariasodlyzko\n")
        return xlo,binarylo
    else
        x = mdsubsetsum(a,ss,.5,1)
        if ~ismissing(x)
            verbose && print("A solution was found via mdsubsetsum\n")
            return x,true
        elseif binarylo===false
            return xlo,binarylo
        end
    end
    return missing.*Bp[:,1]
end

"""
    x = lagariasodlyzko(a,s)

For a vector of integers `a`, and an integer `s`, try to find a binary
vector `x` such that `x'*a=s`. We use the LLL algorithm to find the
solution. This is not a robust tool, just a demo.

This follows the technique described by Lagarias and Odlyzko  in
"Solving Low-Density Subset Sum Problems"  in Journal of ACM, Jan 1985.
Code based on http://web.eecs.umich.edu/~cpeikert/lic15/lec05.pdf
We can likely get better results using techniques described and referenced in
https://www-almasty.lip6.fr/~joux/pages/papers/ToolBox.pdf

It's odd that permuting the `a` vector in the second example given below
causes the alg to often not find a binary solution. Apparently this is a
common oddity with lattice solvers.

# Examples
```jldoctest
julia> a=[1.5;.5;0;.1;.2]; s=2.2; x,_=lagariasodlyzko(a,s); s-x'*a
0.0

julia> a=[32771,65543,131101,262187,524387,1048759, # from Bremner p 117
          2097523,4195057,8390143,16780259,33560539,
          67121039,134242091,268484171,536968403];

julia> s=891221976; x,_=lagariasodlyzko(a,s); s-x'*a
0.0

julia> N=40;a=rand(1:2^BigInt(256),N);xtrue=rand(Bool,N); s=a'*xtrue;

julia> setprecision(BigFloat,300); x,_=lagariasodlyzko(a,s); s-x'*a
0.0

```
"""
function lagariasodlyzko(a::AbstractArray{Td,1},ss::Td,b=0,
                         verbose=false) where {Td<:Number}
    # page numbers below refer to lecture note above

    n = length(a)

    Ti = getIntType(Td)
    if ss<sum(a)/2
        flag=true
        s = sum(a)-ss
    else
        flag = false
        s = ss
    end
    if b≤0  # b is the "B" parameter defined at bottom of page 2
        nt = Ti(n)
        b = ceil(sqrt(nt*2^nt))
    end
    # BB is the 'B' matrix defined also at bottom of page 2
    BB = [Matrix{Ti}(I,n,n+1); -b*a' b*s]
    B,_ = lll(BB)
    #B = l2(BB)

    # Pick off the interesting parts of the B matrix
    Bp = B[1:n,1:n+1]
    # Check if any of the columns satisfy the requirement
    sums = abs.(a'*Bp)
    ixMatch = findall(sums.≈s)
    if length(ixMatch)==0
        @printf("A match that is close according to `isapprox` was not found.\n")
        err,ixx = findmin(sums .-s)
        if abs(err)<1e-4
            @printf("A solution w err %e will be used. Check that it's correct \n",
                    Float64(err))
            ixMatch = findall(sums.-s .≈ err) # lazy hack to find ixMatch
        else
            xr = zeros(Ti,n+1,1).*missing
            return xr, false
        end
    end

    binarySolution = false
    xr = zeros(Ti,n+1,1)
    for ix=1:length(ixMatch)
        # Pick off a column of B
        x = Bp[:,ixMatch[ix][2]]
        if !(any(x.<0)|any(x.>1))

            binarySolution= true
            xr = flag ? mod.(x .+ 1,2) : x
        end
        x = -x
        if !(any(x.<0)|any(x.>1))
            binarySolution= true
            xr = flag ? mod.(x .+ 1,2) : x
        end
    end

    if binarySolution
        verbose && print("A binary Lagarias-Odlyzko solution was found.\n")
    else
        if maximum(a)<2^(n^2/2) && verbose
            density = n/maximum(log2.(a))
            @printf("The density (%4.2f) of 'a' is not as low as required\n", density)
            @printf("(%4.2f) for Lagarias-Odlyzko to work. \n",2/n)
        end

        if length(ixMatch)>=1
            verbose && print("A non-binary Lagarias-Odlyzko solution was "*
                             "found; check that it's correct\n")
            xr=Bp[:,ixMatch[1][2]]
        else
            verbose && print("A Lagarias-Odlyzko solution was not found\n")
            xr= missing.*xr
        end
    end
    return xr, binarySolution
end

"""
    x = mdsubsetsum(a,sM,ratio=.5,Kpm=3)

For a vector of integers `a`, and an integer `sM`, try to find a binary
vector `x` such that `x'*a=s` using the technique from "Multidimensional
subset sum problem" [1][2]. A major goal of the technique is to solve
problems in which there are about 50% ones in `x`; other ratios of ones to zeros
can be specified in `ratio`.  The thesis also suggests searching `Kpm=3`
values around the nominal k. This technique is related to that in
[`subsetsum`](@ref) in that both use the LLL algorithm.  This is not a
robust tool, just a demo.

[1] https://scholarworks.rit.edu/theses/64/

[2] https://pdfs.semanticscholar.org/21a7/c2f9ff29507f1153aefcca04d1cd308e45c0.pdf

# Examples
```jldoctest
julia> using Random, LLLplus

julia> a=[1.5;.5;0;.1;.2]; s=2.2; x=mdsubsetsum(a,s); s-x'*a
0.0

julia> a=[32771,65543,131101,262187,524387,1048759, # from Bremner p 117
          2097523,4195057,8390143,16780259,33560539,
          67121039,134242091,268484171,536968403];

julia> sM=891221976; x=mdsubsetsum(a,sM); sM-x'*a
0

julia> setprecision(BigFloat,300); Random.seed!(0);

julia> N=40;a=rand(1:2^BigInt(256),N);xtrue=rand(Bool,N); s=a'*xtrue;

julia> x=mdsubsetsum(a,s); s-x'*a
0
```
"""
function mdsubsetsum(a::AbstractArray{Td,1},sM::Td,ratio=.5,Kpm=3,
                     verbose=false) where {Td<:Number}

    Ti = getIntType(Td)
    n = Ti(length(a))

    # r is heuristic value for r from end of Section 11 of thesis above
    r = Ti(ceil(2^round(.65*log2(maximum(a)))))
    c = 2^10 # a heuristic; look right below eq (10)

    p = a .÷ r # See start of Section 10
    s = a-p.*r
    k0 = Ti(round(sum(s)*ratio/r)) # Equation (9) with ratio=.5 and j=1.
    p0 = sM ÷ r
    m = sM - p0*r

    for k = k0-Kpm:k0+Kpm
        # Equation (14)
        Bt = [Matrix{Ti}(I,n,n)*2 c*s       c*p      zeros(Ti,n);
              ones(Ti,1,n)        c*(k*r+m) c*(p0-k) 1]
        BB = Bt'  # Notation in paper is row-based; my brain is column-based
        B,_ = lll(BB)
        #B = l2(BB)

        for nx = 1:n
            if abs(B[n+3,nx])==1 && B[n+2,nx]==0 && B[n+1,nx]==0 &&
                      all((B[1:n,nx] .==1) .| (B[1:n,nx] .==-1))
                xhat = abs.(B[1:n,nx] .- B[n+3,nx])/2
                return Ti.(xhat)
            end
        end
    end
    verbose && print("Solution not found")
    return missing
end

"""
    rationalapprox(x::AbstractArray{<:Real,1},M,Ti=BigInt,verbose=false)

For a vector of Reals `x`, and an integer `M`, find an integer q such that
`maximum(abs.(x*q-round.(x*q)))` is small; the vector `x` is approximated by
`round.(x*q)//q`.  The integer `q` is less than or equal to `M` and the
approximation satisfies `max(abs.(x*q-round.(x*q)))≤sqrt(5)*2^(n/4 -
5)*M^(-1/n)`; this equation comes from the paper below.  The LLL algorithm
reduction is used to find the solution. The approximation vector is returned.
This is also known as "simultaneous diophantine approximation"; see for
example the title of the Hanrot paper below.

This is not a robust tool, just a demo.

"LLL: A Tool for Effective Diophantine Approximation" by Guillaume Hanrot in
the book "The LLL Algorithm: Survey and Applications" edited by Phong
Q. Nguyen and Brigitte Vallée, Springer, Heidelberg, 2010.

See also Chapter 9 of M. R. Bremner, "Lattice Basis Reduction: An
Introduction to the LLL Algorithm and Its Applications" CRC Press, 2012.

# Examples
```jldoctest
julia> x = [0.3912641745333527; 0.5455179974014548; 0.1908698210882469];

julia> rationalapprox(x,1e4,Int64)
3-element Vector{Rational{Int64}}:
 43//110
  6//11
 21//110

```
"""
function rationalapprox(x::AbstractArray{<:Real,1},M,Ti=BigInt,verbose=false)
    # Follows the multivariate technique in Theorem 8 of the
    # Nguyen and Vallee paper above

    n = length(x)
    Q = (M/2^(1/4))^(1/n)
    C = Ti(ceil(Q^(n+1)))
    Cx = Ti.(round.(C*x))

    X = [one(Ti) zeros(Ti,1,n);
         [Cx     C*I] ]
    B,_ = lll(X)
    q = abs(B[1,1])
    maxErr = maximum(abs.(x-round.(x*q)/q))
    bound = sqrt(5)*2^(n/4 - 1)/Q /q
    if maxErr>bound
        warning("The max error in the approximation $(maxErr) is unexpectedly "*
                "greater than the bound $(bound).")
    end
    p = Ti.(round.(x*q));
    if verbose
        println("q        = $(q)\n"*
                "maxErr/q = $(maxErr)\n"*
                "bound/q  = $(bound)")
        display([x round.(x*q)/q abs.(x-round.(x*q)/q)])
    end
    return p//q
end

"""
    spigotBBPvec(Td::Type{Tr},s,b,n,K) where {Tr<:Number}

This is an auxiliary function for `spigotBBP`. It calculates the numbers
that `spigotBBP` uses to check for a BBP spigot infinite series.

If you want to try to get BBP coefficients without the aid of the
`spigotBBP` function, below is code you can play with.

# Example
```jldoctest
julia> vv=LLLplus.spigotBBPvec(Float64,1,16,8,20);

julia> v=[vv;-pi];

julia> b=1_000_000;

julia> A=[I; b* v'];

julia> B,T=lll(A);

julia> B[1:end-2,1]'*vv -pi <1e-14 # The error is small
true
```
"""
function spigotBBPvec(Td::Type{Tr},s,b,n,K) where {Tr<:Number}
    v = Vector{Td}(undef,n)
    k = Td(0):K
    for j =1:n
        v[j]= sum(1 ./(b.^k .*(k.*n .+j).^s))
    end
    return v
end

"""
    spigotBBP(α::Td,s,b,n,K,verbose=false) where {Td}

Check for a BBP-style [1] infinite series for the constant `α`.  These are
"spigot" formulas that can be used to generate (for example) the millionth
digit of the constant `α` without learning the previous
digits. Specifically, given the constant `α`, and parameters `b`, `n`, and
`s`, look for a vector of numbers `a_1` through `a_n` that satisfies the
following equation:

``\\alpha= \\sum_{k=0}^\\infty \\frac{1}{b^k} \\left( \\frac{a_1}{(nk+1)^s} + \\ldots + \\frac{a_n}{(nk+n)^s} \\right)``

Because it's hard to sum to infinity, the sum is stopped at K.
If a formula is found, it is printed to the screen in LaTeX and the
coefficents `a` are returned as a vector.  An online LaTeX viewer
such as https://www.latex4technics.com/ may be helpful.

This is not a robust tool, just a demo. For example, there may be a
problem with s≥2. See [2] for derivation of the technique used, and to
check whether a formula you find is new.

[1] David Bailey, Peter Borwein, and Simon Plouffe. "On the rapid
computation of various polylogarithmic constants." Mathematics of
Computation 66.218 (1997): 903-913.
https://www.ams.org/journals/mcom/1997-66-218/S0025-5718-97-00856-9/

[2] David Bailey, "A Compendium of BBP-Type Formulas for Mathematical
Constants". https://www.davidhbailey.com//dhbpapers/bbp-formulas.pdf

# Example
```jldoctest
julia> spigotBBP(BigFloat(pi),1,16,8,45,true);
A solution was found w error -4.728672e-60. In LaTeX form it is
\\alpha= \\sum_{k=0}^\\infty \\frac{1}{16^k} \\left(\\frac{4}{8k+1}-\\frac{2}{8k+4}-\\frac{1}{8k+5}-\\frac{1}{8k+6}\\right)
```

Other examples without output:
```julia
spigotBBP(Float64(pi),1,-4,4,22,true);
spigotBBP(log(2),1,2,2,30,true);
spigotBBP(9*log(3),1,9,2,30,true);
spigotBBP(atan(2)*8,1,16,8,30,true);
spigotBBP(8*sqrt(2)*log(1+sqrt(2)),1,16,8,25,true);
```

There is a formula for pi^2 which the following command should find, but it
does not find it. In fact the technique doesn't seem to work at all for
s>2; It's not obvious what the problem is
```julia
spigotBBP(BigFloat(pi)*pi,2,64,6,25,true);
```
"""
function spigotBBP(α::Td,s,b,n,K,verbose=false) where {Td}
    v = spigotBBPvec(Td,s,b,n,K)
    bN = 100_000_000
    av,_ = lagariasodlyzko(v,α,bN);
    if verbose && ~ismissing(av[1])
        @printf("A solution was found w error %e. In LaTeX form it is\n",av'*v-α)
        @printf("\\alpha= \\sum_{k=0}^\\infty \\frac{1}{%d^k} \\left(",b)
        for ix = 1:length(av)
            if abs(av[ix])>10*eps()
                if sign(av[ix])==1
                    ix>1 && @printf("+")
                else
                    @printf("-")
                end
                if s==1
                    @printf("\\frac{%d}{%dk+%d}",abs(Int64(av[ix])),n,ix)
                else
                    @printf("\\frac{%d}{(%dk+%d)^%d}",abs(Int64(av[ix])),n,ix,s)
                end
            end
        end
        @printf("\\right)\n")
    end
    if ismissing(av[1])
        verbose && @printf("A solution was not found.\n")
        return missing
    else
        return av
    end
end

max_stars_repo_name: christianpeel/LLLplus

max_issues_repo_name: christianpeel/LLLplus

max_forks_repo_name: christianpeel/LLLplus


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Segment(p1, p2)

A line segment with points `p1`, `p2`.
"""
struct Segment{Dim,T,V<:AbstractVector{Point{Dim,T}}} <: Polytope{1,Dim,T}
  vertices::V
end

measure(s::Segment) = norm(s.vertices[2] - s.vertices[1])

facets(s::Segment) = (v for v in s.vertices)


text:
function ODEStateAndDerivative(arg0::jdouble, arg1::Vector{jdouble}, arg2::Vector{jdouble})
    return ODEStateAndDerivative((jdouble, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2)
end

function ODEStateAndDerivative(arg0::jdouble, arg1::Vector{jdouble}, arg2::Vector{jdouble}, arg3::Vector{Vector{jdouble}}, arg4::Vector{Vector{jdouble}})
    return ODEStateAndDerivative((jdouble, Vector{jdouble}, Vector{jdouble}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}), arg0, arg1, arg2, arg3, arg4)
end

function get_complete_derivative(obj::ODEStateAndDerivative)
    return jcall(obj, "getCompleteDerivative", Vector{jdouble}, ())
end

function get_primary_derivative(obj::ODEStateAndDerivative)
    return jcall(obj, "getPrimaryDerivative", Vector{jdouble}, ())
end

function get_secondary_derivative(obj::ODEStateAndDerivative, arg0::jint)
    return jcall(obj, "getSecondaryDerivative", Vector{jdouble}, (jint,), arg0)
end


text:
using Test
using QuantumOpticsBase
using LinearAlgebra, Random

@testset "states" begin

Random.seed!(0)

D(x1::Number, x2::Number) = abs(x2-x1)
D(x1::StateVector, x2::StateVector) = norm(x2-x1)

b1 = GenericBasis(3)
b2 = GenericBasis(5)
b3 = GenericBasis(7)
b = b1 ⊗ b2

bra = Bra(b)
ket = Ket(b)

# Test creation
@test_throws DimensionMismatch Bra(b, [1, 2])
@test_throws DimensionMismatch Ket(b, [1, 2])
@test 0 ≈ norm(bra)
@test 0 ≈ norm(ket)
@test_throws QuantumOpticsBase.IncompatibleBases bra*Ket(b1)
@test bra == bra
@test length(bra) == length(bra.data) == 15
@test length(ket) == length(ket.data) == 15
@test QuantumOpticsBase.basis(ket) == b
@test QuantumOpticsBase.basis(bra) == b
@test bra != basisstate(b, 1)

# Test copy
psi1 = randstate(b1)
psi2 = copy(psi1)
@test psi1.data == psi2.data
@test !(psi1.data === psi2.data)
psi2.data[1] = complex(10.)
@test psi1.data[1] != psi2.data[1]

# Arithmetic operations
# =====================
bra_b1 = dagger(randstate(b1))
bra_b2 = dagger(randstate(b2))
bra_b3 = randstate(b3)'

ket_b1 = randstate(b1)
ket_b2 = randstate(b2)
ket_b3 = randstate(b3)

# Addition
@test_throws DimensionMismatch bra_b1 + bra_b2
@test_throws DimensionMismatch ket_b1 + ket_b2
@test 1e-14 > D(bra_b1 + Bra(b1), bra_b1)
@test 1e-14 > D(ket_b1 + Ket(b1), ket_b1)
@test 1e-14 > D(bra_b1 + dagger(ket_b1), dagger(ket_b1) + bra_b1)

# Subtraction
@test_throws DimensionMismatch bra_b1 - bra_b2
@test_throws DimensionMismatch ket_b1 - ket_b2
@test 1e-14 > D(bra_b1 - Bra(b1), bra_b1)
@test 1e-14 > D(ket_b1 - Ket(b1), ket_b1)
@test 1e-14 > D(bra_b1 - dagger(ket_b1), -dagger(ket_b1) + bra_b1)

# Multiplication
@test 1e-14 > D(-3*ket_b1, 3*(-ket_b1))
@test 1e-14 > D(0.3*(bra_b1 - dagger(ket_b1)), 0.3*bra_b1 - dagger(0.3*ket_b1))
@test 1e-14 > D(0.3*(bra_b1 - dagger(ket_b1)), bra_b1*0.3 - dagger(ket_b1*0.3))
@test 0 ≈ bra*ket
@test 1e-14 > D((bra_b1 ⊗ bra_b2)*(ket_b1 ⊗ ket_b2), (bra_b1*ket_b1)*(bra_b2*ket_b2))

# Tensor product
@test tensor(ket_b1) == ket_b1
@test 1e-14 > D((ket_b1 ⊗ ket_b2) ⊗ ket_b3, ket_b1 ⊗ (ket_b2 ⊗ ket_b3))
@test 1e-14 > D((bra_b1 ⊗ bra_b2) ⊗ bra_b3, bra_b1 ⊗ (bra_b2 ⊗ bra_b3))

ket_b1b2 = ket_b1 ⊗ ket_b2
shape = (ket_b1b2.basis.shape...,)
idx = LinearIndices(shape)[2, 3]
@test ket_b1b2.data[idx] == ket_b1.data[2]*ket_b2.data[3]
ket_b1b2b3 = ket_b1 ⊗ ket_b2 ⊗ ket_b3
@test ket_b1b2b3 == tensor(ket_b1, ket_b2, ket_b3)
shape = (ket_b1b2b3.basis.shape...,)
idx = LinearIndices(shape)[1, 4, 3]
@test ket_b1b2b3.data[idx] == ket_b1.data[1]*ket_b2.data[4]*ket_b3.data[3]

# Norm
bf = FockBasis(1)
bra = Bra(bf, [3.0im, -4])
ket = Ket(bf, [-4.0im, 3])
@test 5 ≈ norm(bra)
@test 5 ≈ norm(ket)

bra_normalized = normalize(bra)
ket_normalized = normalize(ket)
@test 5 ≈ norm(bra)
@test 5 ≈ norm(ket)
@test 1 ≈ norm(bra_normalized)
@test 1 ≈ norm(ket_normalized)

bra_copy = deepcopy(bra)
ket_copy = deepcopy(ket)
normalize!(bra_copy)
normalize!(ket_copy)
@test 5 ≈ norm(bra)
@test 5 ≈ norm(ket)
@test 1 ≈ norm(bra_copy)
@test 1 ≈ norm(ket_copy)

# Test basis state
b1 = GenericBasis(2)
b2 = GenericBasis(3)
b = b1 ⊗ b2
x1 = basisstate(b1, 2)
x2 = basisstate(b2, 1)

@test norm(x1) == 1
@test x1.data[2] == 1
@test basisstate(b, [2, 1]) == x1 ⊗ x2

# Test permutating systems
b1 = GenericBasis(2)
b2 = GenericBasis(5)
b3 = FockBasis(3)

psi1 = randstate(b1)
psi2 = randstate(b2)
psi3 = randstate(b3)

psi123 = psi1 ⊗ psi2 ⊗ psi3
psi132 = psi1 ⊗ psi3 ⊗ psi2
psi213 = psi2 ⊗ psi1 ⊗ psi3
psi231 = psi2 ⊗ psi3 ⊗ psi1
psi312 = psi3 ⊗ psi1 ⊗ psi2
psi321 = psi3 ⊗ psi2 ⊗ psi1

@test 1e-14 > D(psi132, permutesystems(psi123, [1, 3, 2]))
@test 1e-14 > D(psi213, permutesystems(psi123, [2, 1, 3]))
@test 1e-14 > D(psi231, permutesystems(psi123, [2, 3, 1]))
@test 1e-14 > D(psi312, permutesystems(psi123, [3, 1, 2]))
@test 1e-14 > D(psi321, permutesystems(psi123, [3, 2, 1]))

@test 1e-14 > D(dagger(psi132), permutesystems(dagger(psi123), [1, 3, 2]))
@test 1e-14 > D(dagger(psi213), permutesystems(dagger(psi123), [2, 1, 3]))
@test 1e-14 > D(dagger(psi231), permutesystems(dagger(psi123), [2, 3, 1]))
@test 1e-14 > D(dagger(psi312), permutesystems(dagger(psi123), [3, 1, 2]))
@test 1e-14 > D(dagger(psi321), permutesystems(dagger(psi123), [3, 2, 1]))

# Test Broadcasting
@test_throws QuantumOpticsBase.IncompatibleBases psi123 .= psi132
@test_throws QuantumOpticsBase.IncompatibleBases psi123 .+ psi132
bra123 = dagger(psi123)
bra132 = dagger(psi132)
@test_throws ArgumentError psi123 .+ bra123
@test_throws QuantumOpticsBase.IncompatibleBases bra123 .= bra132
@test_throws QuantumOpticsBase.IncompatibleBases bra123 .+ bra132
psi_ = copy(psi123)
psi_ .+= psi123
@test psi_ == 2*psi123
bra_ = copy(bra123)
bra_ .= 3*bra123
@test bra_ == 3*dagger(psi123)
@test_throws ErrorException cos.(psi_)
@test_throws ErrorException cos.(bra_)

end # testset


text:
# Julia wrapper for header: /usr/include/scip/pub_cutpool.h
# Automatically generated using Clang.jl wrap_c

function SCIPcutGetRow()
    ccall((:SCIPcutGetRow, libscip), Ptr{Cint}, ())
end

function SCIPcutGetAge(cut)
    ccall((:SCIPcutGetAge, libscip), Cint, (Ptr{SCIP_CUT},), cut)
end

function SCIPcutGetLPActivityQuot(cut)
    ccall((:SCIPcutGetLPActivityQuot, libscip), Cdouble, (Ptr{SCIP_CUT},), cut)
end

function SCIPcutpoolGetCuts(cutpool)
    ccall((:SCIPcutpoolGetCuts, libscip), Ptr{Ptr{SCIP_CUT}}, (Ptr{SCIP_CUTPOOL},), cutpool)
end

function SCIPcutpoolGetNCuts(cutpool)
    ccall((:SCIPcutpoolGetNCuts, libscip), Cint, (Ptr{SCIP_CUTPOOL},), cutpool)
end

function SCIPcutpoolGetMaxNCuts(cutpool)
    ccall((:SCIPcutpoolGetMaxNCuts, libscip), Cint, (Ptr{SCIP_CUTPOOL},), cutpool)
end

function SCIPcutpoolGetTime(cutpool)
    ccall((:SCIPcutpoolGetTime, libscip), Cdouble, (Ptr{SCIP_CUTPOOL},), cutpool)
end

function SCIPcutpoolGetNCalls(cutpool)
    ccall((:SCIPcutpoolGetNCalls, libscip), Clonglong, (Ptr{SCIP_CUTPOOL},), cutpool)
end

function SCIPcutpoolGetNCutsFound(cutpool)
    ccall((:SCIPcutpoolGetNCutsFound, libscip), Clonglong, (Ptr{SCIP_CUTPOOL},), cutpool)
end


text:
"""
    fcfs()

Implementation of first come first served scheduling algorithm

# Output
Tuple of vectors (process_id, burst_time, waiting_time, turnaround_time, avg_waiting_time, avg_turnaround_time)

# Example

```julia
n = 3 # number of processes
process_id = Any[1, 2, 3] # process ids
burst_times = Any[3, 4, 5] # burst times
fcfs(n, process_id, burst_times)
```
# Reference

https://en.wikipedia.org/wiki/Scheduling_(computing)#First_come,_first_served

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""
function fcfs(n, process_id, burst_time)
    waiting_time = []
    push!(waiting_time, 0)
    for i = 2:n
        # Calculates waiting time
        push!(waiting_time, burst_time[i-1] + waiting_time[i-1])
    end

    turnaround_time = []
    for i = 1:n
        # Calculates turnaround time
        push!(turnaround_time, burst_time[i] + waiting_time[i])
    end

    # Calculates Average waiting time
    avg_waiting_time = sum(waiting_time) / length(waiting_time)

    # Calculates Average turnaround time
    avg_turnaround_time = sum(turnaround_time) / length(turnaround_time)

    return process_id, burst_time, waiting_time, turnaround_time, avg_waiting_time, avg_turnaround_time
end

max_stars_repo_name: arubhardwaj/Julia

max_issues_repo_name: arubhardwaj/Julia

max_forks_repo_name: arubhardwaj/Julia


text:
## NODES

struct Node{T<:Union{Machine, Nothing}} <: AbstractNode

    operation   # eg, `predict` or a static operation, such as `exp`
    machine::T  # is `nothing` for static operations

    # nodes called to get args for `operation(model, ...) ` or
    # `operation(...)`:
    args::Tuple{Vararg{AbstractNode}}

    # sources of ancestor graph (training edges excluded)
    origins::Vector{Source}

    # all ancestors (training edges included) listed in
    # order consistent with extended graph, excluding self
    nodes::Vector{AbstractNode}

    function Node(operation,
                  machine::T,
                  args::AbstractNode...) where T<:Union{Machine, Nothing}

        # check the number of arguments:
        # if machine === nothing && isempty(args)
        #     error("`args` in `Node(::Function, args...)` must be non-empty. ")
        # end

        origins_ = unique(vcat([origins(arg) for arg in args]...))
        # length(origins_) == 1 ||
        #     @warn "A node referencing multiple origins when called " *
        #           "has been defined:\n$(origins_). "

        # initialize the list of upstream nodes:
        nodes_ = AbstractNode[]

        # merge the lists from arguments:
        nodes_ =
            vcat(AbstractNode[], (nodes(n) for n in args)...) |> unique

        # merge the lists from training arguments:
        if machine !== nothing
            nodes_ =
                vcat(nodes_, (nodes(n) for n in machine.args)...) |> unique
        end

        return new{T}(operation, machine, args, origins_, nodes_)
    end
end

"""

    nrows_at_source(N::node)

Return the number of rows of data wrapped at the source of `N`,
assumming this is unique.

Not to be confused with `J = nrows(N)`, which is a new node such that
`J() = nrows(N())`.

See also [`nrows`](@ref)

"""
function nrows_at_source(X::Node)
    ss = sources(X)
    length(ss) == 1 ||
        error("Node does not have a unique source. ")
    return nrows_at_source(first(ss))
end

"""
    origins(N)

Return a list of all origins of a node `N` accessed by a call `N()`.
These are the source nodes of ancestor graph of `N` if edges
corresponding to training arguments are excluded. A `Node` object
cannot be called on new data unless it has a unique origin.

Not to be confused with `sources(N)` which refers to the same graph
but without the training edge deletions.

See also: [`node`](@ref), [`source`](@ref).
"""
origins(X::Node) = X.origins

"""
    nodes(N)

Return all nodes upstream of a node `N`, including `N` itself, in an
order consistent with the extended directed acyclic graph of the
network. Here "extended" means edges corresponding to training
arguments are included.

*Warning.* Not the same as `N.nodes`, which may not include `N`
 itself.

"""
nodes(X::Node) = AbstractNode[X.nodes..., X]

color(N::Node{Nothing}) = :green
color(N::Node) = (N.machine.frozen ? :red : :green)

# constructor for static operations:
Node(operation, args::AbstractNode...) = Node(operation, nothing, args...)

_check(y::Node) = nothing
_check(y::Node{Nothing}) = length(y.origins) == 1 ? nothing :
    error("Node $y has multiple origins and cannot be called "*
              "on new data. ")

# make nodes callable:
(y::Node)(; rows=:) = _apply((y, y.machine); rows=rows)
(y::Node)(Xnew) = (_check(y); _apply((y, y.machine), Xnew))
(y::Node{Nothing})(; rows=:) = _apply((y, ); rows=rows)
(y::Node{Nothing})(Xnew)= (_check(y); _apply((y, ), Xnew))

function _apply(y_plus, input...; kwargs...)
    y = y_plus[1]
    mach = y_plus[2:end] # in static case this is ()
    raw_args = map(y.args) do arg
        arg(input...; kwargs...)
    end
    try
        (y.operation)(mach..., raw_args...)
    catch exception
        @error "Failed "*
        "to apply the operation `$(y.operation)` to the machine "*
        "$(y.machine), which receives it's data arguments from one or more "*
        "nodes in a learning network. Possibly, one of these nodes "*
        "is delivering data that is incompatible with the machine's model.\n"*
        diagnostics(y, input...; kwargs...)
        throw(exception)
    end
end

MLJScientificTypes.elscitype(N::Node) = Unknown
function MLJScientificTypes.elscitype(
    N::Node{<:Machine{<:Union{Deterministic,Unsupervised}}})
    if N.operation == MLJBase.predict
        return target_scitype(N.machine.model)
    elseif N.operation == MLJBase.transform
        return output_scitype(N.machine.model)
    elseif N.operation == MLJBase.inverse_transform
        return input_scitype(N.machine.model)
    end
    return Unknown
end

# TODO after
# https://github.com/alan-turing-institute/ScientificTypes.jl/issues/102 :
# Add Probabilistic case to above

MLJScientificTypes.scitype(N::Node) = CallableReturning{elscitype(N)}

## FITTING A NODE

# flush a (possibly remote) channel"
GenericChannel{T} = Union{Channel{T}, Distributed.RemoteChannel{<:Channel{T}}}
function flush!(c::GenericChannel{T}) where T
    ret = T[]
    while isready(c)
        push!(ret, take!(c))
    end
    return ret
end

"""
    fit!(N::Node;
         rows=nothing,
         verbosity=1,
         force=false,
         acceleration=CPU1())

Train all machines required to call the node `N`, in an appropriate
order.  These machines are those returned by
`machines(N)`.

"""
fit!(y::Node; acceleration=CPU1(), kwargs...) =
    fit!(y::Node, acceleration; kwargs...)

fit!(y::Node, ::AbstractResource; kwargs...) =
        error("Only `acceleration=CPU1()` currently supported")

function fit!(y::Node, ::CPU1; kwargs...)

    _machines = machines(y)

    # flush the fit_okay channels:
    @sync begin
        for mach in _machines
            @async flush!(mach.fit_okay)
        end
    end

    # fit the machines asynchronously;
    @sync begin
        for mach in _machines
            @async fit_only!(mach, true; kwargs...)
        end
    end

    return y
end
fit!(S::Source; args...) = S

# allow arguments of `Nodes` and `Machine`s to appear
# at REPL:
istoobig(d::Tuple{AbstractNode}) = length(d) > 10

# overload show method

_formula(stream::IO, X::AbstractNode, indent) =
    (print(stream, repeat(' ', indent));_formula(stream, X, 0, indent))
_formula(stream::IO, X::Source, depth, indent) = show(stream, X)
function _formula(stream, X::Node, depth, indent)
    operation_name = string(typeof(X.operation).name.mt.name)
    anti = max(length(operation_name) - INDENT)
    print(stream, operation_name, "(")
    n_args = length(X.args)
    if X.machine !== nothing
        print(stream, crind(indent + length(operation_name) - anti))
        printstyled(IOContext(stream, :color=>SHOW_COLOR),
#                        handle(X.machine),
                        X.machine,
                        bold=SHOW_COLOR)
        n_args == 0 || print(stream, ", ")
    end
    for k in 1:n_args
        print(stream, crind(indent + length(operation_name) - anti))
        _formula(stream, X.args[k],
                        depth + 1,
                        indent + length(operation_name) - anti )
        k == n_args || print(stream, ",")
    end
    print(stream, ")")
end

function Base.show(io::IO, ::MIME"text/plain", X::Node)
#    description = string(typeof(X).name.name)
    #    str = "$description $(handle(X))"
    println(io, "$X")
    println(io, "  args:")
    for i in eachindex(X.args)
        arg = X.args[i]
        println(io, "    $i:\t$arg")
    end
    print(io, "  formula:\n")
    _formula(io, X, 4)
    # print(io, " ")
    # printstyled(IOContext(io, :color=>SHOW_COLOR),
    #             handle(X),
    #             color=color(X))
end

## REPORTS AND FITRESULTS FOR NODES

# Both of these exposed but not intended for public use

# here `f` is `report` or `fitted_params`; returns a named tuple:
function item_given_machine(f, N)
    machs = machines(N) |> reverse
    items = map(machs) do m
        try
            f(m)
        catch exception
            if exception isa UndefRefError
                error("UndefRefError intercepted. Perhaps "*
                      "you forgot to `fit!` a machine or node?")
            else
                throw(exception)
            end
        end
    end
    key = f isa typeof(report) ?
        :report_given_machine :
        :fitted_params_given_machine
    dict = LittleDict(machs[j] => items[j] for j in eachindex(machs))
    return NamedTuple{(:machines, key)}((machs, dict))
end

report(N::Node) = item_given_machine(report, N)
report(::Source) = NamedTuple()

MLJModelInterface.fitted_params(N::Node) =
    item_given_machine(fitted_params, N)
MLJModelInterface.fitted_params(S::Source) = NamedTuple()

## SYNTACTIC SUGAR FOR LEARNING NETWORKS

"""
    N = node(f::Function, args...)

Defines a `Node` object `N` wrapping a static operation `f` and arguments
`args`. Each of the `n` elements of `args` must be a `Node` or `Source`
object. The node `N` has the following calling behaviour:

    N() = f(args[1](), args[2](), ..., args[n]())
    N(rows=r) = f(args[1](rows=r), args[2](rows=r), ..., args[n](rows=r))
    N(X) = f(args[1](X), args[2](X), ..., args[n](X))

    J = node(f, mach::Machine, args...)

Defines a dynamic `Node` object `J` wrapping a dynamic operation `f`
(`predict`, `predict_mean`, `transform`, etc), a nodal machine `mach` and
arguments `args`. Its calling behaviour, which depends on the outcome of
training `mach` (and, implicitly, on training outcomes affecting its
arguments) is this:

    J() = f(mach, args[1](), args[2](), ..., args[n]())
    J(rows=r) = f(mach, args[1](rows=r), args[2](rows=r), ..., args[n](rows=r))
    J(X) = f(mach, args[1](X), args[2](X), ..., args[n](X))

Generally `n=1` or `n=2` in this latter case.

    predict(mach, X::AbsractNode, y::AbstractNode)
    predict_mean(mach, X::AbstractNode, y::AbstractNode)
    predict_median(mach, X::AbstractNode, y::AbstractNode)
    predict_mode(mach, X::AbstractNode, y::AbstractNode)
    transform(mach, X::AbstractNode)
    inverse_transform(mach, X::AbstractNode)

Shortcuts for `J = node(predict, mach, X, y)`, etc.

Calling a node is a recursive operation which terminates in the call
to a source node (or nodes). Calling nodes on *new* data `X` fails unless the
number of such nodes is one.

See also: [`@node`](@ref), [`source`](@ref), [`origins`](@ref).

"""
const node = Node

"""
    @node f(...)

Construct a new node that applies the function `f` to some combination
of nodes, sources and other arguments.

*Important.* An argument not in global scope is assumed to be a node
 or source.

### Examples

```
X = source(π)
W = @node sin(X)
julia> W()
0

X = source(1:10)
Y = @node selectrows(X, 3:4)
julia> Y()
3:4

julia> Y([:one, :two, :three, :four])
2-element Array{Symbol,1}:
 :three
 :four

X1 = source(4)
X2 = source(5)
add(a, b, c) = a + b + c
N = @node add(X1, 1, X2)
julia> N()
10

```

See also [`node`](@ref)

"""
macro node(ex)
    ex.head == :call || error("@node syntax error")
    exs = ex.args
    f_ex = first(exs)
    arg_exs = exs[2:end]

    # build lambda expression lambda_left -> lambda_right
    stuff =
        first.(map(arg_exs) do ex
                   pair = (:nothing, false)
                   try
                       evaluated = __module__.eval(ex)
                       if evaluated isa AbstractNode
                           pair =  gensym("node"), true
                       else
                           pair = ex, false
                       end
                   catch e
                       if e isa UndefVarError
                           pair = gensym("node"), true
                       else
                           error()
                       end
                   end
               end |> zip)
    right = first.(stuff)
    mask = last.(stuff)
    left = right[mask]
    lambda_left = Expr(:tuple, left...)
    lambda_right = Expr(:call, f_ex, right...)
    lambda_ex = Expr(:->, lambda_left, lambda_right)

    # the node-only arguments:
    node_arg_exs = arg_exs[mask]

    esc(quote
        node($lambda_ex, $(node_arg_exs...))
    end)
end

"""
    glb(N1, N2, ...)

Construct a node `N` with the behaviour `N() = (N1(), N2(),
...)`. That is, `glb` is `tuple` overloaded for nodes.

Equivalent to `@tuple N1 N2 ...`

"""
glb(X::AbstractNode...) = node(tuple, X...)

"""
    @tuple N1 N2 ...

Construct a node `N` whose calling behaviour is `N() = (N1(), N2(), ...)`.

"""
macro tuple(ex...)
    esc(quote
        glb($(ex...))
        end)
end

"""
    nrows(X::AbstractNode)

Return a new node `N` such that `N() = nrows(X())` and `N(rows=rows) =
nrows(X(rows=rows))`. To obtain the number of rows of data at the
source of `X`, use `nrows_at_source(X)`.

"""
MLJModelInterface.nrows(X::AbstractNode) = node(nrows, X)

MMI.matrix(X::AbstractNode)      = node(matrix, X)
MMI.table(X::AbstractNode)       = node(table, X)
Base.vcat(args::AbstractNode...) = node(vcat, args...)
Base.hcat(args::AbstractNode...) = node(hcat, args...)
Statistics.mean(X::AbstractNode)   = node(v->mean.(v), X)
Statistics.median(X::AbstractNode) = node(v->median.(v), X)
StatsBase.mode(X::AbstractNode)    = node(v->mode.(v), X)
Base.log(X::AbstractNode) = node(v->log.(v), X)
Base.exp(X::AbstractNode) = node(v->exp.(v), X)

+(y1::AbstractNode, y2::AbstractNode) = node(+, y1, y2)
+(y1, y2::AbstractNode) = node(+, y1, y2)
+(y1::AbstractNode, y2) = node(+, y1, y2)
*(lambda::Real, y::AbstractNode) = node(y->lambda*y, y)

"""
    selectcols(X::AbstractNode, c)

Returns `Node` object `N` such that `N() = selectcols(X(), c)`.
"""
MMI.selectcols(X::AbstractNode, r) = node(XX->selectcols(XX, r), X)

"""
    selectrows(X::AbstractNode, r)

Returns a `Node` object `N` such that `N() = selectrows(X(), r)` (and
`N(rows=s) = selectrows(X(rows=s), r)`).

"""
MMI.selectrows(X::AbstractNode, r) = node(XX->selectrows(XX, r), X)

# for accessing and setting model hyperparameters at node:
getindex(n::Node{<:Machine{<:Model}}, s::Symbol) =
    getproperty(n.machine.model, s)
setindex!(n::Node{<:Machine{<:Model}}, v, s::Symbol) =
    setproperty!(n.machine.model, s, v)


text:
function MOI.get(
   model::Optimizer,
   ::MOI.ObjectiveSense
)
   return model.sense
end

function MOI.set(
   model::Optimizer,
   ::MOI.ObjectiveSense,
   sense::MOI.OptimizationSense
)
   model.sense = sense
   return
end

function MOI.set(
   model::Optimizer,
   ::MOI.ObjectiveFunction,
   func::Union{MOI.SingleVariable, MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction}
)
   check_inbounds(model, func)
   model.objective = func
   return
end

function MOI.get(
   model::Optimizer,
   ::MOI.ObjectiveValue
)
   return model.obj
end

function MOI.get(
   model::Optimizer,
   ::MOI.ObjectiveBound
)
   return model.obj_est
end


text:
using Random
using InteractBase, Observables, OrderedCollections, Colors, WebIO, CSSUtil
using Widgets
using Dates
import InteractBase: widgettype
import Widgets: components

using Test

include("test_observables.jl")
include("test_theme.jl")
include("test_deps.jl")


text:
#First delete any existing nc-files
fn1 = tempname2()
fn2 = tempname2()
fn3 = tempname2()
fn4 = tempname2()
fn5 = tempname2()

# Test Medium level Interface
# Test Dimension Creation
d1 = NcDim("Dim1",2;values=[5.0,10.0],atts=Dict("units"=>"deg C"));
d2 = NcDim("Dim2",collect(1:10));
d3 = NcDim("Dim3",20;atts=Dict("max"=>10));
d4 = NcDim("DimUnlim",0,unlimited=true)
d5 = NcDim("str_len",15)

# Test Variable creation
v1 = NcVar("v1",[d1,d2,d3],compress=5)  # With several dims in an Array, and compressed
v2 = NcVar("v2",[d1,d2,d3],atts=Dict("a1"=>"varatts"))  # with given attributes
v3 = NcVar("v3",d1)
vs = NcVar("vstr",d2,t=String)
vc = NcVar("vchar",[d5,d2],t=NetCDF.NC_CHAR)
vscal = NcVar("vscal",NcDim[])

tlist = [Float64, Float32, Int32, Int16, Int8]
vt = Array{NcVar}(undef,length(tlist))
for i= 1:length(tlist)
      vt[i]=NcVar("vt$i",d2,t = tlist[i])
end
vunlim = NcVar("vunlim",d4,t=Float64)

# Creating Files
nc1 = NetCDF.create(fn1,[v1,vs,vc],mode=NC_NETCDF4);
nc2 = NetCDF.create(fn2,NcVar[v2,v3,vscal],gatts=Dict("Some global attributes"=>2010));
nc3 = NetCDF.create(fn3,vt);
ncunlim = NetCDF.create(fn4,vunlim)

#Test Adding attributes
NetCDF.putatt(nc1,"v1",Dict("Additional String attribute"=>"att"))
NetCDF.putatt(nc1,"global",Dict("Additional global attribute"=>"gatt"))
NetCDF.putatt(nc1,"global",Dict("Additional Int8 attribute"=>Int8(20),
                                "Additional Int16 attribute"=>Int16(20),
                                "Additional Int32 attribute"=>Int32(20),
                                "Additional Float32 attribute"=>Float32(20),
                                "Additional Float64 attribute"=>Float64(20)))
NetCDF.putatt(nc1,"v1", Dict("Additional Int8 array attribute"=>Int8[i for i in 1:20],
                             "Additional Int16 array attribute"=>Int16[i for i in 1:20],
                             "Additional Int32 array attribute"=>Int32[i for i in 1:20],
                             "Additional Float32 array attribute"=>Float32[i for i in 1:20],
                             "Additional Float64 array attribute"=>Float64[i for i in 1:20]))

#Test writing data
#
#First generate the data
x1 = rand(2,10,20)
x2 = rand(2,10,20)
xt = [rand(tl,10) for tl in tlist]
xs = ["a","bb","ccc","dddd","eeeee","ffffff","ggggggg","hhhhhhhh","iiiiiiiii","jjjjjjjjjj"]
xscal = Array{Float64,0}(undef);xscal[1]=2.5
#
# And write it
#
NetCDF.putvar(nc1,"v1",x1)
NetCDF.putvar(nc1,"vstr",xs)
NetCDF.putvar(nc1,"vchar",nc_string2char(xs))
#Test sequential writing along one dimension
for i = 1:10
  NetCDF.putvar(nc2,"v2",x2[:,i,:],start=[1,i,1],count=[-1,1,-1])
end
#Test automatic type conversion
x4 = [1,2]
NetCDF.putvar(nc2,"v3",x4)

NetCDF.putvar(nc2,"vscal",xscal)

for i=1:length(tlist)
  NetCDF.putvar(nc3,"vt$i",xt[i])
end

NetCDF.putvar(ncunlim,"vunlim",collect(1:10))

NetCDF.close(nc1)
NetCDF.close(nc2)
NetCDF.close(nc3)

##Read the data back and compare
nc1 = NetCDF.open(fn1,mode=NC_NOWRITE);
nc2 = NetCDF.open(fn2,mode=NC_NOWRITE);
nc3 = NetCDF.open(fn3,mode=NC_NOWRITE);

@test x1 == NetCDF.readvar(nc1,"v1")
@test xs == NetCDF.readvar(nc1,"vstr")
@test x2 == NetCDF.readvar(nc2,"v2")
@test x4 == NetCDF.readvar(nc2,"v3")
@test xscal == NetCDF.readvar(nc2,"vscal")

@test xs == NetCDF.nc_char2string(NetCDF.readvar(nc1,"vchar"))

#Test -1 reading full dimension
NetCDF.readvar(nc1,"v1",start=[1,1,1],count=[-1,-1,-1])

#Test dimension variables type
dim1 = NcDim("dim1", 4, atts=Dict("units" => "m"), values=collect(Int32,1:4))
x5_1 = NcVar("x5_1", [dim1], atts=Dict("units" => "m"), t=Int32)

dim2 = NcDim("dim2", 4, atts=Dict("units" => "m"), values=collect(Float32, 1:4))
x5_2 = NcVar("x5_2", [dim2], atts=Dict("units" => "m"), t=Int32)

nc5 = NetCDF.create(fn5, Array{NcVar, 1}([x5_1, x5_2]))
NetCDF.close(nc5)

nc5 = NetCDF.open(fn5)
@test typeof(NetCDF.readvar(nc5, "dim1")) == Array{Int32, 1}
@test typeof(NetCDF.readvar(nc5, "dim2")) == Array{Float32, 1}


text:
using EffectiveWaves, Test
using LinearAlgebra

# this is for high volume and frequency
@testset "compare wienger-hopf and match method average scattering coefficient" begin

    medium = Medium(1.0,1.0+0.0im)

    ω = 1.0

    specie = Specie(ρ=0.0, r=0.5, c=0.0, volfrac = 0.30)
    radius_multiplier = 1.001

    tol = 1e-8
    hankel_order=0
    θ = pi/4

    k_effs = wavenumbers(ω, medium, [specie];
        radius_multiplier = radius_multiplier,
        tol=tol, hankel_order = hankel_order,
        num_wavenumbers = 50);

    wave_effs = [
        EffectiveWave(ω, k_eff, medium, [specie];
            hankel_order = hankel_order,
            radius_multiplier = radius_multiplier,
            tol = tol, extinction_rescale=false,
            method = :WienerHopf,
            θin=θ
        )
    for k_eff in k_effs];

    match_ws = MatchWave(ω, medium, specie;
        radius_multiplier = radius_multiplier,
        hankel_order = hankel_order,
        tol = tol, wave_effs = wave_effs[1:10],
        θin=θ,
        max_size = 800,
    );

    inds = Int.(round.(LinRange(1,length(wave_effs),5)))
    errors = map(inds) do i
        avg_WH = AverageWave(match_ws.average_wave.x[1:200], wave_effs[1:i]);
        norm(avg_WH.amplitudes - match_ws.average_wave.amplitudes[1:200,:,:])/norm(match_ws.average_wave.amplitudes[1:200,:,:])
    end

    # errors should be monotonically decreasing
    @test sort(errors; rev=true) == errors
    @test errors[end] < 0.05
end


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

struct WorkflowsApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_workflowsStepCompleted(_api::WorkflowsApi, workflow_step_execute_id::String; outputs=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/workflows.stepCompleted", ["slackAuth"])
    Swagger.set_param(_ctx.query, "workflow_step_execute_id", workflow_step_execute_id)  # type String
    Swagger.set_param(_ctx.query, "outputs", outputs)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Indicate that an app's step in a workflow completed execution.
Param: workflow_step_execute_id::String (required)
Param: outputs::String
Return: Dict{String, Any}
"""
function workflowsStepCompleted(_api::WorkflowsApi, workflow_step_execute_id::String; outputs=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_workflowsStepCompleted(_api, workflow_step_execute_id; outputs=outputs, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function workflowsStepCompleted(_api::WorkflowsApi, response_stream::Channel, workflow_step_execute_id::String; outputs=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_workflowsStepCompleted(_api, workflow_step_execute_id; outputs=outputs, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_workflowsStepFailed(_api::WorkflowsApi, workflow_step_execute_id::String, error::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/workflows.stepFailed", ["slackAuth"])
    Swagger.set_param(_ctx.query, "workflow_step_execute_id", workflow_step_execute_id)  # type String
    Swagger.set_param(_ctx.query, "error", error)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Indicate that an app's step in a workflow failed to execute.
Param: workflow_step_execute_id::String (required)
Param: error::String (required)
Return: Dict{String, Any}
"""
function workflowsStepFailed(_api::WorkflowsApi, workflow_step_execute_id::String, error::String; _mediaType=nothing)
    _ctx = _swaggerinternal_workflowsStepFailed(_api, workflow_step_execute_id, error; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function workflowsStepFailed(_api::WorkflowsApi, response_stream::Channel, workflow_step_execute_id::String, error::String; _mediaType=nothing)
    _ctx = _swaggerinternal_workflowsStepFailed(_api, workflow_step_execute_id, error; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_workflowsUpdateStep(_api::WorkflowsApi, workflow_step_edit_id::String; inputs=nothing, outputs=nothing, step_name=nothing, step_image_url=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/workflows.updateStep", ["slackAuth"])
    Swagger.set_param(_ctx.query, "workflow_step_edit_id", workflow_step_edit_id)  # type String
    Swagger.set_param(_ctx.query, "inputs", inputs)  # type String
    Swagger.set_param(_ctx.query, "outputs", outputs)  # type String
    Swagger.set_param(_ctx.query, "step_name", step_name)  # type String
    Swagger.set_param(_ctx.query, "step_image_url", step_image_url)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Update the configuration for a workflow extension step.
Param: workflow_step_edit_id::String (required)
Param: inputs::String
Param: outputs::String
Param: step_name::String
Param: step_image_url::String
Return: Dict{String, Any}
"""
function workflowsUpdateStep(_api::WorkflowsApi, workflow_step_edit_id::String; inputs=nothing, outputs=nothing, step_name=nothing, step_image_url=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_workflowsUpdateStep(_api, workflow_step_edit_id; inputs=inputs, outputs=outputs, step_name=step_name, step_image_url=step_image_url, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function workflowsUpdateStep(_api::WorkflowsApi, response_stream::Channel, workflow_step_edit_id::String; inputs=nothing, outputs=nothing, step_name=nothing, step_image_url=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_workflowsUpdateStep(_api, workflow_step_edit_id; inputs=inputs, outputs=outputs, step_name=step_name, step_image_url=step_image_url, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export workflowsStepCompleted, workflowsStepFailed, workflowsUpdateStep


text:
module TensorCore

using LinearAlgebra

export ⊙, hadamard, hadamard!
export ⊗, tensor, tensor!
export ⊡, boxdot, boxdot!

"""
    hadamard(a, b)
    a ⊙ b

For arrays `a` and `b`, perform elementwise multiplication.
`a` and `b` must have identical `axes`.

`⊙` can be passed as an operator to higher-order functions.

# Examples
```jldoctest; setup=:(using TensorCore)
julia> a = [2, 3]; b = [5, 7];

julia> a ⊙ b
2-element Array{$Int,1}:
 10
 21

julia> a ⊙ [5]
ERROR: DimensionMismatch("Axes of `A` and `B` must match, got (Base.OneTo(2),) and (Base.OneTo(1),)")
[...]
```

See also `hadamard!(y, a, b)`.
"""
function hadamard(A::AbstractArray, B::AbstractArray)
    @noinline throw_dmm(axA, axB) = throw(DimensionMismatch("Axes of `A` and `B` must match, got $axA and $axB"))

    axA, axB = axes(A), axes(B)
    axA == axB || throw_dmm(axA, axB)
    return map(*, A, B)
end
const ⊙ = hadamard

"""
    hadamard!(dest, A, B)

Similar to `hadamard(A, B)` (which can also be written `A ⊙ B`), but stores its results in
the pre-allocated array `dest`.
"""
function hadamard!(dest::AbstractArray, A::AbstractArray, B::AbstractArray)
    @noinline function throw_dmm(axA, axB, axdest)
        throw(DimensionMismatch("`axes(dest) = $axdest` must be equal to `axes(A) = $axA` and `axes(B) = $axB`"))
    end

    axA, axB, axdest = axes(A), axes(B), axes(dest)
    ((axdest == axA) & (axdest == axB)) || throw_dmm(axA, axB, axdest)
    @simd for I in eachindex(dest, A, B)
        @inbounds dest[I] = A[I] * B[I]
    end
    return dest
end

"""
    tensor(A, B)
    A ⊗ B

Compute the tensor product of `A` and `B`.
If `C = A ⊗ B`, then `C[i1, ..., im, j1, ..., jn] = A[i1, ... im] * B[j1, ..., jn]`.

For vectors `v` and `w`, the Kronecker product is related to the tensor product by
`kron(v,w) == vec(w ⊗ v)` or `w ⊗ v == reshape(kron(v,w), (length(w), length(v)))`.

# Examples
```jldoctest; setup=:(using TensorCore)
julia> a = [2, 3]; b = [5, 7, 11];

julia> a ⊗ b
2×3 Array{$Int,2}:
 10  14  22
 15  21  33
```
See also `tensor!(Y,A,B)`.
"""
tensor(A::AbstractArray, B::AbstractArray) = [a*b for a in A, b in B]
const ⊗ = tensor

const CovectorLike{T} = Union{Adjoint{T,<:AbstractVector},Transpose{T,<:AbstractVector}}
function tensor(u::AbstractArray, v::CovectorLike)
    # If `v` is thought of as a covector, you might want this to be two-dimensional,
    # but thought of as a matrix it should be three-dimensional.
    # The safest is to avoid supporting it at all. See discussion in #35150.
    error("`tensor` is not defined for co-vectors, perhaps you meant `*`?")
end
function tensor(u::CovectorLike, v::AbstractArray)
    error("`tensor` is not defined for co-vectors, perhaps you meant `*`?")
end
function tensor(u::CovectorLike, v::CovectorLike)
    error("`tensor` is not defined for co-vectors, perhaps you meant `*`?")
end

"""
    tensor!(dest, A, B)

Similar to `tensor(A, B)` (which can also be written `A ⊗ B`), but stores its results in
the pre-allocated array `dest`.
"""
function tensor!(dest::AbstractArray, A::AbstractArray, B::AbstractArray)
    @noinline function throw_dmm(axA, axB, axdest)
        throw(DimensionMismatch("`axes(dest) = $axdest` must concatenate `axes(A) = $axA` and `axes(B) = $axB`"))
    end

    axA, axB, axdest = axes(A), axes(B), axes(dest)
    axes(dest) == (axA..., axB...) || throw_dmm(axA, axB, axdest)
    if IndexStyle(dest) === IndexCartesian()
        for IB in CartesianIndices(axB)
            @inbounds b = B[IB]
            @simd for IA in CartesianIndices(axA)
                @inbounds dest[IA,IB] = A[IA]*b
            end
        end
    else
        i = firstindex(dest)
        @inbounds for b in B
            @simd for a in A
                dest[i] = a*b
                i += 1
            end
        end
    end
    return dest
end

export boxdot, ⊡, boxdot!

"""
    boxdot(A,B) = A ⊡ B    # \\boxdot

Generalised matrix multiplication: Contracts the last dimension of `A` with
the first dimension of `B`, for any `ndims(A)` & `ndims(B)`.
If both are vectors, then it returns a scalar `== sum(A .* B)`.

# Examples
```jldoctest; setup=:(using TensorCore)
julia> A = rand(3,4,5); B = rand(5,6,7);

julia> size(A ⊡ B)
(3, 4, 6, 7)

julia> typeof(rand(5) ⊡ rand(5))
Float64

julia> try B ⊡ A catch err println(err) end
DimensionMismatch("neighbouring axes of `A` and `B` must match, got Base.OneTo(7) and Base.OneTo(3)")
```
This is the same behaviour as Mathematica's function `Dot[A, B]`.
It is not identicaly to Python's `numpy.dot(A, B)`, which contracts with the second-last
dimension of `B` instead of the first, but both keep all the other dimensions.
Unlike Julia's `LinearAlgebra.dot`, it does not conjugate `A`, so these two agree only
for real-valued vectors.

When interacting with `Adjoint` vectors, this always obeys `(x ⊡ y)' == y' ⊡ x'`,
and hence may sometimes return another `Adjoint` vector. (And similarly for `Transpose`.)

```jldoctest; setup=:(using TensorCore)
julia> M = rand(5,5); v = rand(5);

julia> typeof(v ⊡ M')
Array{Float64,1}

julia> typeof(M ⊡ v')  # adjoint of the previous line
Adjoint{Float64,Array{Float64,1}}

julia> typeof(v' ⊡ M')  # same as *, and equal to adjoint(M ⊡ v)
Adjoint{Float64,Array{Float64,1}}

julia> typeof(v' ⊡ v)
Float64
```
See also `boxdot!(Y,A,B)`, which is to `⊡` as `mul!` is to `*`.
"""
function boxdot(A::AbstractArray, B::AbstractArray)
    Amat = _squash_left(A)
    Bmat = _squash_right(B)

    axA, axB = axes(Amat,2), axes(Bmat,1)
    axA == axB || _throw_dmm(axA, axB)

    return _boxdot_reshape(Amat * Bmat, A, B)
end

const ⊡ = boxdot

@noinline _throw_dmm(axA, axB) = throw(DimensionMismatch("neighbouring axes of `A` and `B` must match, got $axA and $axB"))

_squash_left(A::AbstractArray) = reshape(A, :,size(A,ndims(A)))
_squash_left(A::AbstractMatrix) = A

_squash_right(B::AbstractArray) = reshape(B, size(B,1),:)
_squash_right(B::AbstractVecOrMat) = B

function _boxdot_reshape(AB::AbstractArray, A::AbstractArray{T,N}, B::AbstractArray{S,M}) where {T,N,S,M}
    ax = ntuple(i -> i<N ? axes(A, i) : axes(B, i-N+2), Val(N+M-2))
    reshape(AB, ax) # some cases don't come here, so this doesn't really support OffsetArrays
end

# These can skip final reshape:
_boxdot_reshape(AB::AbstractVecOrMat, A::AbstractMatrix, B::AbstractVecOrMat) = AB

# These produce scalar output:
function boxdot(A::AbstractVector, B::AbstractVector)
    axA, axB = axes(A,1), axes(B,1)
    axA == axB || _throw_dmm(axA, axB)
    if eltype(A) <: Number
        return transpose(A)*B
    else
        return sum(a*b for (a,b) in zip(A,B))
    end
end

# Multiplication by a scalar:
boxdot(A::AbstractArray, b::Number) = A*b
boxdot(a::Number, B::AbstractArray) = a*B
boxdot(a::Number, b::Number) = a*b

using LinearAlgebra: AdjointAbsVec, TransposeAbsVec, AdjOrTransAbsVec

# Adjont and Transpose, vectors or almost (returning a scalar)
boxdot(A::AdjointAbsVec, B::AbstractVector) = A * B
boxdot(A::TransposeAbsVec, B::AbstractVector) = A * B

boxdot(A::AbstractVector, B::AdjointAbsVec) = A ⊡ vec(B)
boxdot(A::AbstractVector, B::TransposeAbsVec) = A ⊡ vec(B)

boxdot(A::AdjointAbsVec, B::AdjointAbsVec) = adjoint(adjoint(B) ⊡ adjoint(A))
boxdot(A::AdjointAbsVec, B::TransposeAbsVec) = vec(A) ⊡ vec(B)
boxdot(A::TransposeAbsVec, B::AdjointAbsVec) = vec(A) ⊡ vec(B)
boxdot(A::TransposeAbsVec, B::TransposeAbsVec) = transpose(transpose(B) ⊡ transpose(A))

# ... with a matrix (returning another such)
boxdot(A::AdjointAbsVec, B::AbstractMatrix) = A * B
boxdot(A::TransposeAbsVec, B::AbstractMatrix) = A * B

boxdot(A::AbstractMatrix, B::AdjointAbsVec) = (B' ⊡ A')'
boxdot(A::AbstractMatrix, B::TransposeAbsVec) = transpose(transpose(B) ⊡ transpose(A))

# ... and with higher-dim (returning a plain array)
boxdot(A::AdjointAbsVec, B::AbstractArray) = vec(A) ⊡ B
boxdot(A::TransposeAbsVec, B::AbstractArray) = vec(A) ⊡ B

boxdot(A::AbstractArray, B::AdjointAbsVec) = A ⊡ vec(B)
boxdot(A::AbstractArray, B::TransposeAbsVec) = A ⊡ vec(B)

"""
    boxdot!(Y, A, B, α=1, β=0)

In-place version of `boxdot`, i.e. `Y .= (A ⊡ B) .* β .+ Y .* α`.
Like 5-argument `mul!`, the use of `α, β` here requires Julia 1.3 or later.
"""
function boxdot! end

if VERSION < v"1.3" # Then 5-arg mul! isn't defined

    function boxdot!(Y::AbstractArray, A::AbstractArray, B::AbstractArray)
        szY = prod(size(A)[1:end-1]), prod(size(B)[2:end])
        mul!(reshape(Y, szY), _squash_left(A), _squash_right(B))
        Y
    end

    boxdot!(Y::AbstractArray, A::AbstractArray, B::AdjOrTransAbsVec) = boxdot!(Y, A, vec(B))

else

    function boxdot!(Y::AbstractArray, A::AbstractArray, B::AbstractArray, α::Number=true, β::Number=false)
        szY = prod(size(A)[1:end-1]), prod(size(B)[2:end])
        mul!(reshape(Y, szY), _squash_left(A), _squash_right(B), α, β)
        Y
    end

    # For boxdot!, only where mul! behaves differently:
    boxdot!(Y::AbstractArray, A::AbstractArray, B::AdjOrTransAbsVec,
        α::Number=true, β::Number=false) = boxdot!(Y, A, vec(B))

end

"""
    TensorCore._adjoint(A)

This extends `adjoint` to understand higher-dimensional arrays, always reversing the
order of dimensions. On Julia 1.5 and later, the symbol `'` can be overloaded locally
as `var"'"`, as shown below.

Then `(x ⊡ y)' == y' ⊡ x'` holds for `x` and `y` arrays of any dimension.

# Examples
```jldoctest; setup=:(using TensorCore)
julia> T3 = rand(3,4,5); v = rand(5);

julia> size(T3 ⊡ v')
(3, 4)

julia> let var"'" = TensorCore._adjoint
         v ⊡ T3' ≈ (T3 ⊡ v')'
       end
true
```
"""
_adjoint(x) = adjoint(x)
_adjoint(x::AbstractVecOrMat) = adjoint(x)
_adjoint(x::AbstractArray{T,N}) where {T<:Number,N} = conj(PermutedDimsArray(x, ntuple(i -> N-i+1, N)))
_adjoint(x::AbstractArray{T,N}) where {T,N} = adjoint.(PermutedDimsArray(x, ntuple(i -> N-i+1, N)))

_transpose(x) = transpose(x)
_transpose(x::AbstractVecOrMat) = transpose(x)
_transpose(x::AbstractArray{T,N}) where {T<:Number,N} = PermutedDimsArray(x, ntuple(i -> N-i+1, N))
_transpose(x::AbstractArray{T,N}) where {T,N} = transpose.(PermutedDimsArray(x, ntuple(i -> N-i+1, N)))

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*-------------------------------------------------------------------
#* Author: Rafael de Pelegrini Soares
#*         Maur좩cio Carvalho Maciel
#* $Id$
#*--------------------------------------------------------------------
type mixer3
	mixer3()=begin
		new(
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream 1",
				:PosX=>0,
				:PosY=>0.25,
				:Symbol=>"_{in1}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream 2",
				:PosX=>0,
				:PosY=>0.5059,
				:Symbol=>"_{in2}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream 3",
				:PosX=>0,
				:PosY=>0.75,
				:Symbol=>"_{in3}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>1,
				:PosY=>0.5059,
				:Symbol=>"_{out}"
			)),
			[
				:(Outlet.F = Inlet1.F + Inlet2.F + Inlet3.F),
				:(Outlet.F*Outlet.z = Inlet1.F*Inlet1.z + Inlet2.F*Inlet2.z + Inlet3.F*Inlet3.z),
				:(Outlet.F*Outlet.h = Inlet1.F*Inlet1.h + Inlet2.F*Inlet2.h + Inlet3.F*Inlet3.h),
				:(Outlet.P = min([Inlet1.P, Inlet2.P , Inlet3.P])),
			],
			[
				"Flow","Composition","Energy Balance","Pressure",
			],
			[:Inlet1,:Inlet2,:Inlet3,:Outlet,]
		)
	end
	Inlet1::stream
	Inlet2::stream
	Inlet3::stream
	Outlet::streamPH
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export mixer3
function setEquationFlow(in::mixer3)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
end
function atributes(in::mixer3,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/mixer"
	fields[:Brief]="Model of a mixer with 2 inputs"
	fields[:Info]="== Assumptions ==
* static
* adiabatic

== Specify ==
* the inlet streams"
	drive!(fields,_)
	return fields
end
mixer3(_::Dict{Symbol,Any})=begin
	newModel=mixer3()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using Cbc
using QuadraticToBinary

const CBC_OPTIMIZER = Cbc.Optimizer()
MOI.set(CBC_OPTIMIZER, MOI.Silent(), true)
const CBC_CACHE = MOIU.UniversalFallback(MOIU.Model{Float64}())
const CBC_CACHED = MOIU.CachingOptimizer(CBC_CACHE, CBC_OPTIMIZER)
const CBC_BRIDGED = MOI.Bridges.full_bridge_optimizer(CBC_CACHED, Float64)

push!(solvers, (opt = CBC_BRIDGED, mode = BilevelJuMP.SOS1Mode()))
push!(solvers_sos, (opt = CBC_BRIDGED, mode = BilevelJuMP.SOS1Mode()))
push!(solvers_indicator, (opt = CBC_BRIDGED, mode = BilevelJuMP.IndicatorMode()))

push!(solvers_bin_exp, (
    opt = QuadraticToBinary.Optimizer{Float64}(CBC_BRIDGED),
    mode = BilevelJuMP.StrongDualityInequalityMode{Float64}(1e-9)))
MOI.set(solvers_bin_exp[end].opt, QuadraticToBinary.GlobalVariablePrecision(), 1e-5)

push!(solvers_bin_exp, (
    opt = QuadraticToBinary.Optimizer{Float64}(CBC_BRIDGED),
    mode = BilevelJuMP.StrongDualityEqualityMode()))
MOI.set(solvers_bin_exp[end].opt, QuadraticToBinary.GlobalVariablePrecision(), 1e-5)

# push!(solvers_sos_quad, (
#     opt = QuadraticToBinary.Optimizer{Float64}(CBC_BRIDGED),
#     mode = BilevelJuMP.SOS1Mode()))
# MOI.set(solvers_bin_exp[end].opt, QuadraticToBinary.GlobalVariablePrecision(), 1e-5)


text:
roundable(num::Number) = isinteger(num) & (typemin(Int) <= num < typemax(Int))

ceil_neg_log10(x) = roundable(-log10(x)) ? ceil(Integer, -log10(x)) : floor(Integer, -log10(x))
round_neg_log10(x) = roundable(-log10(x)) ? round(Integer, -log10(x), RoundNearestTiesUp) : floor(Integer, -log10(x))
round_up_tick(x,m) = x == 0. ? 0.: (x > 0 ? ceil(x, ceil_neg_log10(m)) : -floor(-x, ceil_neg_log10(m)))
round_down_tick(x,m) = x == 0. ? 0.: (x > 0 ? floor(x, ceil_neg_log10(m)) : -ceil(-x, ceil_neg_log10(m)))
round_up_subtick(x,m) = x == 0. ? 0.: (x > 0 ? ceil(x, ceil_neg_log10(m)+1) : -floor(-x, ceil_neg_log10(m)+1))
round_down_subtick(x,m) = x == 0. ? 0.: (x > 0 ? floor(x, ceil_neg_log10(m)+1) : -ceil(-x, ceil_neg_log10(m)+1))
float_round_log10(x::F,m) where {F<:AbstractFloat} = x == 0. ? F(0) : (x > 0 ? round(x, ceil_neg_log10(m)+1)::F : -round(-x, ceil_neg_log10(m)+1)::F)
float_round_log10(x) = x > 0 ? float_round_log10(x,x) : float_round_log10(x,-x)

function plotting_range(xmin, xmax)
    diffX = xmax - xmin
    xmax = round_up_tick(xmax, diffX)
    xmin = round_down_tick(xmin, diffX)
    xmin, xmax
end

function plotting_range_narrow(xmin, xmax)
    diffX = xmax - xmin
    xmax = round_up_subtick(xmax, diffX)
    xmin = round_down_subtick(xmin, diffX)
    xmin, xmax
end

function extend_limits(vec, limits)
    mi = float(minimum(limits))
    ma = float(maximum(limits))
    if mi == 0. && ma == 0.
        mi = minimum(vec)
        ma = maximum(vec)
    end
    diff = ma - mi
    if diff == 0
        ma = mi + 1
        mi = mi - 1
    end
    (limits == [0.,0.]) ? plotting_range_narrow(mi, ma) : plotting_range_narrow(mi, ma)
end

const bordermap = Dict{Symbol,Dict{Symbol,String}}()
const border_solid  = Dict{Symbol,String}()
const border_bold   = Dict{Symbol,String}()
const border_none   = Dict{Symbol,String}()
const border_dashed = Dict{Symbol,String}()
const border_dotted = Dict{Symbol,String}()
const border_ascii  = Dict{Symbol,String}()
border_solid[:tl] = "┌"
border_solid[:tr] = "┐"
border_solid[:bl] = "└"
border_solid[:br] = "┘"
border_solid[:t] = "─"
border_solid[:l] = "│"
border_solid[:b] = "─"
border_solid[:r] = "│"
border_bold[:tl] = "┏"
border_bold[:tr] = "┓"
border_bold[:bl] = "┗"
border_bold[:br] = "┛"
border_bold[:t] = "━"
border_bold[:l] = "┃"
border_bold[:b] = "━"
border_bold[:r] = "┃"
border_none[:tl] = " "
border_none[:tr] = " "
border_none[:bl] = " "
border_none[:br] = " "
border_none[:t] = " "
border_none[:l] = " "
border_none[:b] = " "
border_none[:r] = " "
border_dashed[:tl] = "┌"
border_dashed[:tr] = "┐"
border_dashed[:bl] = "└"
border_dashed[:br] = "┘"
border_dashed[:t] = "╌"
border_dashed[:l] = "│"
border_dashed[:b] = "╌"
border_dashed[:r] = "│"
border_dotted[:tl] = "⡤"
border_dotted[:tr] = "⢤"
border_dotted[:bl] = "⠓"
border_dotted[:br] = "⠚"
border_dotted[:t] = "⠤"
border_dotted[:l] = "⡇"
border_dotted[:b] = "⠒"
border_dotted[:r] = "⢸"
border_ascii[:tl] = "+"
border_ascii[:tr] = "+"
border_ascii[:bl] = "+"
border_ascii[:br] = "+"
border_ascii[:t] = "-"
border_ascii[:l] = "|"
border_ascii[:b] = "-"
border_ascii[:r] = "|"
bordermap[:solid]  = border_solid
bordermap[:bold]   = border_bold
bordermap[:none]   = border_none
bordermap[:dashed] = border_dashed
bordermap[:dotted] = border_dotted
bordermap[:ascii]  = border_ascii

const color_cycle = [:blue, :red, :yellow, :magenta, :green, :cyan]
const color_encode = Dict{Symbol,UInt8}()
const color_decode = Dict{UInt8,Symbol}()
color_encode[:white]   = 0b000
color_encode[:blue]    = 0b001
color_encode[:red]     = 0b010
color_encode[:magenta] = 0b011
color_encode[:yellow]  = 0b100
color_encode[:green]   = 0b101
color_encode[:cyan]    = 0b110
for k in keys(color_encode)
    v = color_encode[k]
    color_decode[v] = k
end
color_decode[0b111] = :white

function print_color(color::UInt8, io::IO, args...)
    col = color_decode[color]
    str = string(args...)
    print_with_color(col, io, str)
end


text:
@testset "NumField/NfAbs" begin
  include("NfAbs/ConjugatesNS.jl")
end


text:
abstract type XTree end

struct XNode{T<:XTree, V} <: XTree
    tag::Symbol
    attributes::Dict{Symbol, V}
    children::Vector{T}
end

Base.show(io::IO, xnode::XNode) = print_tree(io, xnode, 3)

XNode(tag::Symbol) = XNode(tag, XTree[])
XNode(tag::Symbol, attributes::Dict) = XNode(tag, attributes, XTree[])
XNode(tag::Symbol, children::Vector) = XNode(tag, Dict{Symbol, Any}(), children)

AbstractTrees.children(xnode::XNode) = xnode.children
tag(xnode::XNode) = xnode.tag
attributes(xnode::XNode) = xnode.attributes

Base.iterate(xnode::XNode) = iterate(children(xnode))
Base.iterate(xnode::XNode, state) = iterate(children(xnode), state)
Base.IteratorSize(::Type{XNode{T}}) where T = Base.SizeUnknown()
Base.eltype(::Type{XNode{T}}) where T = T

Base.eltype(::Type{<:TreeIterator{XNode{T}}}) where T = T
Base.IteratorEltype(::Type{<:TreeIterator{XNode{T}}}) where T = Base.HasEltype()

withchildren(xnode::XNode, children) = XNode(tag(xnode), attributes(xnode), children)
withtag(xnode::XNode, tag) = XNode(tag, attributes(xnode), children(xnode))
withattributes(xnode::XNode, attributes) = XNode(tag(xnode), attributes, children(xnode))

struct XLeaf{T} <: XTree
    data::T
end

Base.getindex(xleaf::XLeaf) = xleaf.data

AbstractTrees.children(::XLeaf) = ()
AbstractTrees.printnode(io::IO, xleaf::XLeaf) = print(io, xleaf[])

# Catamorphism

"""
    catafold(f, xtree, state)

Fold a function `f : (node, state) -> (node', state')` post-order over `xtree` and
return a modified tree as well as the resulting state.

Use [`cata`](#) for a stateless catamorphism and [`fold`](#) if you don't want to
transform `xtree`.
"""
function catafold(f, xnode::XNode, state; T = XTree)
    cs = T[]
    for c in children(xnode)
        c, state = catafold(f, c, state; T = T)
        push!(cs, c)
    end
    return f(withchildren(xnode, cs), state)
end

catafold(f, xleaf::XLeaf, state; kwargs...) = f(xleaf, state)

"""
    cata(f, xtree)

Replace every node or leaf `v` in `xtree` with `f(v)`. Traverses in
post-order.
"""
cata(f, xnode::XNode) = f(withchildren(xnode, XTree[cata(f, c) for c in children(xnode)]))
cata(f, xleaf::XLeaf) = f(xleaf)

fold(f, xnode::XTree, init) = foldl(f, (l for l in PostOrderDFS(xnode)); init = init)
foldleaves(f, xnode::XTree, init) = foldl(f, (l[] for l in Leaves(xnode)); init = init)

function Base.:(==)(x1::XNode, x2::XNode)
    return ((tag(x1) == tag(x2)) &&
            (attributes(x1) == attributes(x2)) &&
            (length(children(x1)) == length(children(x2))) &&
            all(c1 == c2 for (c1, c2) in zip(children(x1), children(x2))))
end

Base.:(==)(::XNode, _) = false
Base.:(==)(_, ::XNode) = false

Base.:(==)(x1::XLeaf, x2::XLeaf) = x1[] == x2[]
Base.:(==)(::XLeaf, _) = false
Base.:(==)(::XNode, ::XLeaf) = false
Base.:(==)(::XLeaf, ::XNode) = false

# Showing

function AbstractTrees.printnode(io::IO, x::XNode)
    print(io, ":", x.tag)
    if !isempty(x.attributes)
        print(io, " [")
        for (i, (key, value)) in enumerate(x.attributes)
            if i != 1
                print(io, ", ")
            end
            print(io, key, " = ", value)
        end
        print(io, "]")
    end
end


text:
module Solvers

export
    PoissonSolver, PoissonBCs, solve_poisson_3d!,
    BatchedTridiagonalSolver, solve_batched_tridiagonal_system!

using Oceananigans.Grids

using Oceananigans: @hascuda
@hascuda using CUDAnative, CuArrays

abstract type PoissonBCs end

include("solver_utils.jl")
include("poisson_solver_cpu.jl")
include("poisson_solver_gpu.jl")
include("batched_tridiagonal_solver.jl")

end


text:
function optimize(f::Function,
                  initial_x::Array;
                  method = NelderMead(),
                  x_tol::Real = 1e-32,
                  f_tol::Real = 1e-32,
                  g_tol::Real = 1e-8,
                  iterations::Integer = 1_000,
                  store_trace::Bool = false,
                  show_trace::Bool = false,
                  extended_trace::Bool = false,
                  show_every::Integer = 1,
                  autodiff::Bool = false,
                  callback = nothing)
    options = OptimizationOptions(;
        x_tol = x_tol, f_tol = f_tol, g_tol = g_tol,
        iterations = iterations, store_trace = store_trace,
        show_trace = show_trace, extended_trace = extended_trace,
        callback = callback, show_every = show_every,
        autodiff = autodiff)
    optimize(f, initial_x, method, options)
end

function optimize(f::Function,
                  g!::Function,
                  initial_x::Array;
                  method = LBFGS(),
                  x_tol::Real = 1e-32,
                  f_tol::Real = 1e-32,
                  g_tol::Real = 1e-8,
                  iterations::Integer = 1_000,
                  store_trace::Bool = false,
                  show_trace::Bool = false,
                  extended_trace::Bool = false,
                  show_every::Integer = 1,
                  callback = nothing)
    options = OptimizationOptions(;
        x_tol = x_tol, f_tol = f_tol, g_tol = g_tol,
        iterations = iterations, store_trace = store_trace,
        show_trace = show_trace, extended_trace = extended_trace,
        callback = callback, show_every = show_every)
    optimize(f, g!, initial_x, method, options)
end

function optimize(f::Function,
                  g!::Function,
                  h!::Function,
                  initial_x::Array;
                  method = Newton(),
                  x_tol::Real = 1e-32,
                  f_tol::Real = 1e-32,
                  g_tol::Real = 1e-8,
                  iterations::Integer = 1_000,
                  store_trace::Bool = false,
                  show_trace::Bool = false,
                  extended_trace::Bool = false,
                  show_every::Integer = 1,
                  callback = nothing)
    options = OptimizationOptions(;
        x_tol = x_tol, f_tol = f_tol, g_tol = g_tol,
        iterations = iterations, store_trace = store_trace,
        show_trace = show_trace, extended_trace = extended_trace,
        callback = callback, show_every = show_every)
    optimize(f, g!, h!, initial_x, method, options)
end

function optimize(d::DifferentiableFunction,
                  initial_x::Array;
                  method = LBFGS(),
                  x_tol::Real = 1e-32,
                  f_tol::Real = 1e-32,
                  g_tol::Real = 1e-8,
                  iterations::Integer = 1_000,
                  store_trace::Bool = false,
                  show_trace::Bool = false,
                  extended_trace::Bool = false,
                  show_every::Integer = 1,
                  callback = nothing)
    options = OptimizationOptions(;
        x_tol = x_tol, f_tol = f_tol, g_tol = g_tol,
        iterations = iterations, store_trace = store_trace,
        show_trace = show_trace, extended_trace = extended_trace,
        callback = callback, show_every = show_every)
    optimize(d, initial_x, method, options)
end

function optimize(d::TwiceDifferentiableFunction,
                  initial_x::Array;
                  method = Newton(),
                  x_tol::Real = 1e-32,
                  f_tol::Real = 1e-32,
                  g_tol::Real = 1e-8,
                  iterations::Integer = 1_000,
                  store_trace::Bool = false,
                  show_trace::Bool = false,
                  extended_trace::Bool = false,
                  show_every::Integer = 1,
                  callback = nothing)
    options = OptimizationOptions(;
        x_tol = x_tol, f_tol = f_tol, g_tol = g_tol,
        iterations = iterations, store_trace = store_trace,
        show_trace = show_trace, extended_trace = extended_trace,
        callback = callback, show_every = show_every)
    optimize(d, initial_x, method, options)
end

function optimize(d,
                  initial_x::Array,
                  method::Optimizer,
                  options::OptimizationOptions = OptimizationOptions())
    optimize(d, initial_x, method, options)
end

function optimize(f::Function,
                  g!::Function,
                  initial_x::Array,
                  method::Optimizer,
                  options::OptimizationOptions = OptimizationOptions())
    d = DifferentiableFunction(f, g!)
    optimize(d, initial_x, method, options)
end

function optimize(f::Function,
                  g!::Function,
                  h!::Function,
                  initial_x::Array,
                  method::Optimizer,
                  options::OptimizationOptions = OptimizationOptions())
    d = TwiceDifferentiableFunction(f, g!, h!)
    optimize(d, initial_x, method, options)
end

function optimize{T}(f::Function,
                  initial_x::Array{T},
                  method::Optimizer,
                  options::OptimizationOptions)
    if !options.autodiff
        d = DifferentiableFunction(f)
    else
        g!(x, out) = ForwardDiff.gradient!(out, f, x)

        function fg!(x, out)
            gr_res = ForwardDiff.GradientResult(zero(T),out)
            ForwardDiff.gradient!(gr_res, f, x)
            ForwardDiff.value(gr_res)
        end
        d = DifferentiableFunction(f, g!, fg!)
    end
    optimize(d, initial_x, method, options)
end

function optimize{T}(f::Function,
                  initial_x::Array{T},
                  method::Union{Newton, NewtonTrustRegion},
                  options::OptimizationOptions)
    if !options.autodiff
        error("No gradient or Hessian was provided. Either provide a gradient and Hessian, set autodiff = true in the OptimizationOptions if applicable, or choose a solver that doesn't require a Hessian.")
    else
        g!(x, out) = ForwardDiff.gradient!(out, f, x)

        function fg!(x, out)
            gr_res = ForwardDiff.GradientResult(zero(T),out)
            ForwardDiff.gradient!(gr_res, f, x)
            ForwardDiff.value(gr_res)
        end

        h! = (x, out) -> ForwardDiff.hessian!(out, f, x)
        d = TwiceDifferentiableFunction(f, g!, fg!, h!)
    end
    optimize(d, initial_x, method, options)
end

function optimize(f::Function,
                  g!::Function,
                  initial_x::Array,
                  method::Union{Newton, NewtonTrustRegion},
                  options::OptimizationOptions)
    if !options.autodiff
        error("No Hessian was provided. Either provide a Hessian, set autodiff = true in the OptimizationOptions if applicable, or choose a solver that doesn't require a Hessian.")
    else
        function fg!(x, out)
            g!(x, out)
            f(x)
        end

        h! = (x, out) -> ForwardDiff.hessian!(out, f, x)
        d = TwiceDifferentiableFunction(f, g!, fg!, h!)
    end
    optimize(d, initial_x, method, options)
end

function optimize(d::DifferentiableFunction,
                  initial_x::Array,
                  method::Union{Newton, NewtonTrustRegion},
                  options::OptimizationOptions)
    if !options.autodiff
        error("No Hessian was provided. Either provide a Hessian, set autodiff = true in the OptimizationOptions if applicable, or choose a solver that doesn't require a Hessian.")
    else
        h! = (x, out) -> ForwardDiff.hessian!(out, d.f, x)
    end
    optimize(TwiceDifferentiableFunction(d.f, d.g!, d.fg!, h!), initial_x, method, options)
end

function optimize(d::DifferentiableFunction,
                  initial_x::Array,
                  method::Optimizer,
                  options::OptimizationOptions)
    optimize(d.f, initial_x, method, options)
end

function optimize(d::TwiceDifferentiableFunction,
                  initial_x::Array,
                  method::Optimizer,
                  options::OptimizationOptions)
    dn = DifferentiableFunction(d.f, d.g!, d.fg!)
    optimize(dn, initial_x, method, options)
end

function optimize{T <: AbstractFloat}(f::Function,
                                      lower::T,
                                      upper::T;
                                      method = Brent(),
                                      rel_tol::Real = sqrt(eps(T)),
                                      abs_tol::Real = eps(T),
                                      iterations::Integer = 1_000,
                                      store_trace::Bool = false,
                                      show_trace::Bool = false,
                                      callback = nothing,
                                      show_every = 1,
                                      extended_trace::Bool = false)
    show_every = show_every > 0 ? show_every: 1
    if extended_trace && callback == nothing
        show_trace = true
    end
    if show_trace
        @printf "Iter     Function value   Gradient norm \n"
    end
    optimize(f, Float64(lower), Float64(upper), method;
             rel_tol = rel_tol,
             abs_tol = abs_tol,
             iterations = iterations,
             store_trace = store_trace,
             show_trace = show_trace,
             show_every = show_every,
             callback = callback,
             extended_trace = extended_trace)
end

function optimize(f::Function,
                  lower::Real,
                  upper::Real;
                  kwargs...)
    optimize(f,
             Float64(lower),
             Float64(upper);
             kwargs...)
end


text:
using Test
using LinearAlgebra: eigen
using SparseArrays: SparseMatrixCSC
using QuantumLattices: PID, Point, Lattice, FID, Index, OID, Table, OIDToTuple, Metric, Hilbert, Spin, Fock
using QuantumLattices: Operator, OperatorSum, FockTerm, SpinTerm, Hopping, Onsite, Hubbard, Parameters
using QuantumLattices: ⊗, ⊕, dimension, contentnames, getcontent, parameternames, dtype, add!, matrix, idtype, kind, statistics, update!
using ExactDiagonalization

@testset "BinaryBasis" begin
    basis = BinaryBasis(5)
    @test basis==BinaryBasis((1, 3)) && isequal(basis, BinaryBasis{UInt}((1, 3)))
    @test basis<BinaryBasis(6) && isless(basis, BinaryBasis(6))
    @test one(basis) == one(typeof(basis)) == BinaryBasis(1)
    @test zero(basis) == zero(typeof(basis)) == BinaryBasis(0)
    @test string(basis) == "101"
    @test eltype(basis) == eltype(typeof(basis)) == Int
    @test collect(basis) == [1, 3]
    @test one(basis, 2) == BinaryBasis(7)
    @test isone(basis, 1)==true && isone(basis, 2)==false && isone(basis, 3)==true
    @test zero(basis, 1) == BinaryBasis(4)
    @test iszero(basis, 1)==false && iszero(basis, 2)==true && iszero(basis, 3)==false
    @test count(basis, 1, 1)==1 && count(basis, 1, 2)==1 && count(basis, 1, 3)==2 && count(basis, 2, 3)==1
    @test basis==BinaryBasis([1])⊗BinaryBasis([3])
end

@testset "BinaryBasisRange" begin
    bbr = BinaryBasisRange(UInt(0):UInt(4))
    @test issorted(bbr) == true
    @test dimension(bbr) == 5
    for i = 1:dimension(bbr)
        @test bbr[i] == BinaryBasis(UInt(i-1))
    end
end

@testset "BinaryBases" begin
    bs = BinaryBases(2)
    @test issorted(bs) == true
    @test dimension(bs) == length(bs) == 4
    @test bs==BinaryBases(1:2)
    @test isequal(bs, BinaryBases((2, 1)))
    for i = 1:dimension(bs)
        @test bs[i]==BinaryBasis(i-1)
        @test findfirst(bs[i], bs) == i
    end
    @test eltype(bs) == eltype(typeof(bs)) == BinaryBasis{UInt}
    @test collect(bs) == map(BinaryBasis, [0, 1, 2, 3])
    @test repr(bs) == "2^2"
    @test string(bs) == "2^2:\n  0\n  1\n  10\n  11\n"

    bs = BinaryBases(4, 2)
    @test collect(bs) == map(BinaryBasis, [3, 5, 6, 9, 10, 12])
    @test repr(bs) == "C(4, 2)"
    for i = 1:dimension(bs)
        @test findfirst(bs[i], bs) == i
    end

    bs = BinaryBases(1:2, 1) ⊗ BinaryBases(3:4, 1)
    @test collect(bs) == map(BinaryBasis, [5, 6, 9, 10])
    @test repr(bs) == "C(2, 1) ⊗ C(2, 1)"
end

@testset "TargetSpace" begin
    bs₁ = BinaryBases(1:4, 2)
    bs₂ = BinaryBases(1:4, 3)
    ts = TargetSpace(bs₁, bs₂)
    @test contentnames(typeof(ts)) == (:table,)
    @test getcontent(ts, :table) == ts.sectors
    @test add!(add!(TargetSpace(typeof(bs₁)[]), bs₁), bs₂) == add!(TargetSpace(typeof(bs₁)[]), ts) == ts
    @test bs₁⊕bs₂ == ts

    bs₃ = BinaryBases(5:8, 2)
    @test ts⊕bs₃ == TargetSpace(bs₁, bs₂, bs₃)
    @test bs₃⊕ts == TargetSpace(bs₃, bs₁, bs₂)
    @test ts⊕TargetSpace(bs₃) == TargetSpace(bs₁, bs₂, bs₃)
end

@testset "matrix" begin
    oids = [OID(Index(PID(i), FID{:f}(1, 1, 1)), [0.0, 0.0], [0.0, 0.0]) for i = 1:4]
    table = Table(oids, OIDToTuple(:site, :orbital, :spin))

    braket = (BinaryBases(1:4, 2), BinaryBases(1:4, 3))
    ops = [Operator(2.0, oid) for oid in oids]
    @test matrix(ops[1], braket, table) == SparseMatrixCSC(6, 4, [1, 2, 3, 4, 4], [3, 5, 6], [2.0, 2.0, 2.0])
    @test matrix(ops[2], braket, table) == SparseMatrixCSC(6, 4, [1, 2, 3, 3, 4], [2, 4, 6], [-2.0, -2.0, 2.0])
    @test matrix(ops[3], braket, table) == SparseMatrixCSC(6, 4, [1, 2, 2, 3, 4], [1, 4, 5], [2.0, -2.0, -2.0])
    @test matrix(ops[4], braket, table) == SparseMatrixCSC(6, 4, [1, 1, 2, 3, 4], [1, 2, 3], [2.0, 2.0, 2.0])

    braket = (BinaryBases(1:4, 2), BinaryBases(1:4, 2))
    @test matrix(ops[2]'*ops[1], braket, table) == SparseMatrixCSC(6, 6, [1, 1, 2, 2, 3, 3, 3], [3, 5], [4.0, 4.0])
    @test matrix(ops[3]'*ops[1], braket, table) == SparseMatrixCSC(6, 6, [1, 2, 2, 2, 3, 3, 3], [3, 6], [-4.0, 4.0])
end

@testset "EDMatrix" begin
    m = EDMatrix(BinaryBases(1:4, 2), SparseMatrixCSC(6, 6, [1, 1, 2, 2, 3, 3, 3], [3, 5], [4.0, 4.0]))
    @test parameternames(typeof(m)) == (:sector, :value)
    @test getcontent(m, :id) == (m.bra, m.ket)
    @test getcontent(m, :value) == m.matrix
    @test dtype(m) == dtype(typeof(m)) == Float64
    @test promote_type(typeof(m), ComplexF64) == EDMatrix{BinaryBases{BinaryBasis{UInt}, Vector{BinaryBasis{UInt}}}, SparseMatrixCSC{ComplexF64, Int}}
end

@testset "EDMatrixRepresentation && SectorFilter" begin
    oids = [OID(Index(PID(i), FID{:f}(1, 1, 1)), [0.0, 0.0], [0.0, 0.0]) for i = 1:4]
    table = Table(oids, OIDToTuple(:site, :orbital, :spin))
    op₁, op₂, op₃ = Operator(2.0, oids[2]', oids[1]), Operator(2.0, oids[3]', oids[2]), Operator(2.0, oids[4]', oids[3])
    ops = op₁ + op₂ + op₃
    target = BinaryBases(1:4, 1)⊕BinaryBases(1:4, 2)⊕BinaryBases(1:4, 3)
    M = EDMatrix{BinaryBases{BinaryBasis{UInt}, Vector{BinaryBasis{UInt}}}, SparseMatrixCSC{Float64, Int}}

    mr = EDMatrixRepresentation(target, table)
    @test valtype(typeof(mr), eltype(ops)) == valtype(typeof(mr), typeof(ops)) == OperatorSum{M, idtype(M)}

    ms = mr(ops)
    mr₁ = EDMatrixRepresentation(TargetSpace(target[1]), table)
    mr₂ = EDMatrixRepresentation(TargetSpace(target[2]), table)
    mr₃ = EDMatrixRepresentation(TargetSpace(target[3]), table)
    @test ms == mr₁(ops) + mr₂(ops) + mr₃(ops)
    @test mr₁(ops) == mr₁(op₁) + mr₁(op₂) + mr₁(op₃)
    @test mr₂(ops) == mr₂(op₁) + mr₂(op₂) + mr₂(op₃)
    @test mr₃(ops) == mr₃(op₁) + mr₃(op₂) + mr₃(op₃)

    sf = SectorFilter(target[1])
    @test sf == SectorFilter((target[1], target[1]))
    @test valtype(typeof(sf), typeof(ms)) == typeof(ms)
    @test sf(ms) == mr₁(ops)
end

@testset "ED" begin
    @test EDKind(FockTerm) == EDKind(:FED)
    @test EDKind(SpinTerm) == EDKind(:SED)
    @test EDKind(Tuple{Hopping, Onsite, Hubbard}) == EDKind(:FED)

    @test Metric(EDKind(:FED), Hilbert(PID(1)=>Fock{:f}(1, 2, 2))) == OIDToTuple(:spin, :site, :orbital)
    @test Metric(EDKind(:SED), Hilbert(PID(1)=>Spin{1//2}(1))) == OIDToTuple(:site, :orbital)

    lattice = Lattice(:L2P, [Point(PID(1), [0.0]), Point(PID(2), [1.0])])
    hilbert = Hilbert(pid=>Fock{:f}(1, 2, 2) for pid in lattice.pids)
    t = Hopping(:t, 1.0, 1)
    U = Hubbard(:U, 0.0, modulate=true)
    μ = Onsite(:μ, 0.0, modulate=true)

    bases = BinaryBases(1:2, 1)⊗BinaryBases(3:4, 1)
    ed = ED(lattice, hilbert, (t, U, μ), TargetSpace(bases))
    @test kind(ed) == kind(typeof(ed)) == EDKind(:FED)
    @test valtype(ed) == valtype(typeof(ed)) == Float64
    @test statistics(ed) == statistics(typeof(ed)) == :f
    @test Parameters(ed) == (t=1.0, U=0.0, μ=0.0)

    vector = [0.5; -0.5; -0.5; 0.5]
    eigensystem = eigen(matrix(ed); nev=1)
    @test isapprox(eigensystem.values, [-2.0]; atol=10^-10)
    @test isapprox(eigensystem.vectors, vector; atol=10^-10) || isapprox(eigensystem.vectors, -vector; atol=10^-10)

    vector = [-0.43516214649359913; 0.5573454101893041; 0.5573454101893037; -0.43516214649359913]
    update!(ed, U=1.0, μ=-0.5)
    eigensystem = eigen(matrix(ed); nev=1)
    @test isapprox(eigensystem.values, [-2.5615528128088303]; atol=10^-10)
    @test isapprox(eigensystem.vectors, vector; atol=10^-10) || isapprox(eigensystem.vectors, -vector; atol=10^-10)
end


text:
abstract type AbstractBeta{T1} <: AbstractDistributionLH{T1} end
abstract type AbstractBetaSwitches{T1} <: AbstractDistributionSwitches{T1} end

## -------------------  Generic

function draw(u :: AbstractBeta{T1}, nDims, rng :: AbstractRNG) where
    T1 <: AbstractFloat
    betaDistr = Distributions.Beta(u.alpha, u.beta);
    draws = min_value(u) .+ value_range(u) .* rand(rng, betaDistr, nDims...);
    return draws
end

function draw(u :: AbstractBeta{T1}, rng :: AbstractRNG)  where T1 <: AbstractFloat
    betaDistr = Distributions.Beta(u.alpha, u.beta);
    draws = min_value(u) + value_range(u) * rand(rng, betaDistr);
    return draws
end

function quantiles(u :: AbstractBeta{T1}, pctV)  where T1 <: AbstractFloat
    betaDistr = Distributions.Beta(u.alpha, u.beta);
    qV = min_value(u) .+ value_range(u) .* Distributions.quantile.(betaDistr, pctV);
    return qV
end

isbounded(u :: AbstractBeta{T1}) where T1 = true;
min_value(u :: AbstractBeta{T1}) where T1 = u.xMin;
max_value(u :: AbstractBeta{T1}) where T1 = u.xMin + u.xRange;
value_range(u :: AbstractBeta{T1}) where T1 = u.xRange;
alpha_value(u :: AbstractBeta{T1}) where T1 = u.alpha;
beta_value(u :: AbstractBeta{T1}) where T1 = u.beta;

Base.eltype(::AbstractBeta{T1}) where T1 <: AbstractFloat = T1;

Base.show(io :: IO, um :: AbstractBeta{T1}) where T1 =
    print(io, "Beta{$T1} over range $(min_value(um)) to $(max_value(um))");

# Validate draws
function check_draws(u :: AbstractBeta{T1}, drawM) where T1
    isValid = true;
    if any(x -> x < min_value(u), drawM)
        @warn "Below $(min_value(u)):  $(minimum(drawM))"
        isValid = false;
    end

    if any(x -> x > max_value(u), drawM)
        @warn "Above $(max_value(u)):  $(maximum(drawM))"
        isValid = false;
    end
    return isValid
end

init_distribution(objId :: ObjectId, switches :: AbstractBetaSwitches{T1}) where T1 =
    init_beta(objId, switches);

make_test_beta_switches() = BetaSwitches();
make_test_beta() = init_beta(ObjectId(:none),  make_test_beta_switches());

## ----------  Default Beta

"""
	$(SIGNATURES)

Switches for default Beta distribution object.
"""
Base.@kwdef mutable struct BetaSwitches{T1 <: AbstractFloat} <: AbstractBetaSwitches{T1}
    xMin :: T1 = 1.0
    xMinLb :: T1 = 0.1
    xMinUb :: T1 = 100.0
    xMinDescription :: String = "Lower bound"
    xMinLatex :: String = "x_{min}"
    xRange :: T1 = 3.0
    xRangeLb :: T1 = 0.1
    xRangeUb :: T1 = 100.0
    xRangeDescription :: String = "Range"
    xRangeLatex :: String = "dx"
    alpha :: T1 = 2.0
    alphaLb :: T1 = 1.0
    alphaUb :: T1 = 3.0
    alphaDescription :: String = "alpha"
    alphaLatex :: String = "\\alpha"
    beta :: T1 = 2.0
    betaLb :: T1 = 1.0
    betaUb :: T1 = 5.0
    betaDescription :: String = "beta"
    betaLatex :: String = "\\beta"
    calXMin :: Bool = true
    calXRange :: Bool = true
    calAlpha :: Bool = true
    calBeta :: Bool = true
end

"""
	$(SIGNATURES)

Beta marginal distribution. Characterized by lower bound, upper bound, and Beta parameters alpha and beta.
Cannot store the Beta distribution in the object because its parameters are not known.
"""
mutable struct Beta{T1} <: AbstractBeta{T1}
    objId :: ObjectId
    pvec :: ParamVector
	xMin :: T1
    xRange :: T1
    alpha :: T1
    beta :: T1
	# betaDistr :: Distributions.Beta{T1}
end

"""
	$(SIGNATURES)

Construct a Uniform distribution from its switches.
"""
function init_beta(objId :: ObjectId, switches :: BetaSwitches{T1}) where T1
    pMin = init_xmin(switches);
    pRange = init_xrange(switches);
    pAlpha = init_alpha_param(switches);
    pBeta = init_beta_param(switches);
    pvec = ParamVector(objId, [pMin, pRange, pAlpha, pBeta]);
    xMin = ModelParams.value(pMin);
    # betaDistr = Distributions.Beta(alpha, beta);
    return Beta(objId, pvec,
        ModelParams.value(pMin), ModelParams.value(pRange),
        ModelParams.value(pAlpha), ModelParams.value(pBeta));
end

init_xmin(s :: AbstractBetaSwitches{T1}) where T1 =
    Param(:xMin, s.xMinDescription, s.xMinLatex, s.xMin, s.xMin,
        s.xMinLb, s.xMinUb, s.calXMin);

init_xrange(s :: AbstractBetaSwitches{T1}) where T1 =
    Param(:xRange, s.xRangeDescription, s.xRangeLatex, s.xRange, s.xRange,
        s.xRangeLb, s.xRangeUb, s.calXRange);

init_alpha_param(s :: AbstractBetaSwitches{T1}) where T1 =
    Param(:alpha, s.alphaDescription, s.alphaLatex, s.alpha, s.alpha,
        s.alphaLb, s.alphaUb, s.calAlpha);

init_beta_param(s :: AbstractBetaSwitches{T1}) where T1 =
    Param(:beta, s.betaDescription, s.betaLatex, s.beta, s.beta,
        s.betaLb, s.betaUb, s.calBeta);

# --------------