text:
# Disable this because of broken dependencies
# using FastBroadcast: @..

####
#### Transformations to our representations
####

## Caution, we choose that this differs from OpenFermion by a factor of 1/2 in the two-body integral.
## In OpenFermion, they include the factor when constructing InteractionOperator.
## See spinorb_from_spatial in molecular_data.py.
"""
    spin_orbital_from_spatial(one_body_integrals, two_body_integrals; block_spin=false, index_order=:physicist)

Create coefficients for spin-orbitals from coefficients for spatial orbitals only.

If `block_spin` is `true` the spins are in two sectors. Otherwise, the spin variable varies
more rapidly than the spatial variable, i.e. the spins are interleaved; orbitals
corresponding to the same spatial orbital, but with differing spins are adjacent.  Note that
Qiskit uses the former ordering and OpenFermion uses the latter.

`index_order` must be one of `:physicist`, `:chemist`, or `:intermediate`, and must correspond to the
actual symmetry of `two_body_integrals` in order that the spin-orbital tensor be constructed correctly.
This symmetry may be determined via `find_index_order`.
"""
function spin_orbital_from_spatial(one_body_integrals, two_body_integrals; block_spin=false, index_order=:physicist)
    spatial_dim = size(one_body_integrals)[1]
    spin_orb_dim = 2 * spatial_dim
    if block_spin
        r1 = 1:spatial_dim  # first block of orbitals
        r2 = (spatial_dim+1):spin_orb_dim  # second block of orbitals
    else  # interleaved spin order
        r1 = 1:2:spin_orb_dim # odd indices
        r2 = 2:2:spin_orb_dim # even indices
    end
    return _spin_orbital_from_spatial(one_body_integrals, two_body_integrals, r1, r2, index_order)
end

## Create spin-orbital coefficients from spatial-orbital coeffcients. r1 and r2 are indices
## encoding whether spin variable varies faster or slower than spatial-orbital variable.
function _spin_orbital_from_spatial(one_body_integrals, two_body_integrals, r1, r2, index_order=:physicist)
    new_dim = 2 * size(one_body_integrals)[1]

    one_body_coefficients = zeros(eltype(one_body_integrals), (new_dim, new_dim))
    @. one_body_coefficients[r1, r1] = one_body_integrals
    @. one_body_coefficients[r2, r2] = one_body_integrals

    two_body_coefficients = zeros(eltype(two_body_integrals), (fill(new_dim, 4)...,))
    if index_order == :physicist
        index_set = ((r1, r2, r2, r1), (r2, r1, r1, r2), (r1, r1, r1, r1), (r2, r2, r2, r2))
    elseif index_order == :intermediate
        index_set = ((r1, r2, r1, r2), (r2, r1, r2, r1), (r1, r1, r1, r1), (r2, r2, r2, r2))
    elseif index_order == :chemist
        index_set = ((r1, r1, r2, r2), (r2, r2, r1, r1), (r1, r1, r1, r1), (r2, r2, r2, r2))
    else
        throw(ArgumentError("Unknown index order: $index_order"))
    end
    @inbounds for inds in index_set
        @. two_body_coefficients[inds...] = two_body_integrals
    end
    @. two_body_coefficients = two_body_coefficients / 2

    # Should we do zchop ?
    # return ZChop.zchop!.((one_body_coefficients, two_body_coefficients))
    return (one_body_coefficients, two_body_coefficients)
end

"""
    struct InteractionOperator

Structure to represent a molecular Hamiltonian.
"""
struct InteractionOperator
    nuclear_repulsion::Float64
    one_body_tensor::Array{Float64, 2}
    two_body_tensor::Array{Float64, 4}
end

Base.:(==)(iop1::InteractionOperator, iop2::InteractionOperator) =
    iop1.nuclear_repulsion == iop2.nuclear_repulsion &&
    iop1.one_body_tensor == iop2.one_body_tensor &&
    iop1.two_body_tensor == iop2.two_body_tensor

Base.isapprox(iop1::InteractionOperator, iop2::InteractionOperator) =
    isapprox(iop1.nuclear_repulsion, iop2.nuclear_repulsion) &&
    isapprox(iop1.one_body_tensor, iop2.one_body_tensor) &&
    isapprox(iop1.two_body_tensor, iop2.two_body_tensor)

## NOTE !!!! Conversion of InteractionOperator to FermionOperator is in OF conversions.py
## !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
"""
    InteractionOperator(mol_data::MolecularData; spin_order=:block, index_order=:default)

Create an InteractionOperator from `mol_data`.

`mol_data` contains integrals for spatial orbitals. `InteractionOrbitals` creates coefficients for
spin orbitals, doubling the number of orbitals.

If `spin_order` is `:block` the two-body tensor is constructed with spins in two
sectors. If `spin_order` is `:alternating`, the spin variable varies more rapidly than the spatial variable,
i.e. adjacent orbitals have differing spins. Note that Qiskit uses the former ordering and
OpenFermion uses the latter. If `spin_order` is `:none`, then the operator includes only the input
spatial oribtals; no spin orbitals are computed.

If `index_order` is `:default` then it is assumed that the two-body tensor in `mol_data` is
in chemists' index order and that the desired order for the spin-orbital tensor is
physicists' order.  Otherwise, `index_order` must be one of `:chemist`, `:physicist`, or
`:intermediate` and an attempt will be made to construct the spin-orbital tensor in the
corresponding index order.

The default values of `block_spin` and `index_order` agree with the `InteractionOperator` from OpenFermion.
"""
function InteractionOperator(mol_data::MolecularData; spin_order=:block, index_order=:default)
    interaction_operator(mol_data.nuclear_repulsion, mol_data.one_body_integrals, mol_data.two_body_integrals;
                         spin_order=spin_order, index_order=index_order)
end

function interaction_operator(nuclear_repulsion, one_body_integrals, two_body_integrals;
                              spin_order=:block, index_order=:default)
    tb = two_body_integrals
    if index_order == :default
        tens_tmp = chem_to_phys(tb)
        index_order = :physicist
    else
        tens_tmp = to_index_order(tb, index_order)
    end
    if spin_order == :none
        return InteractionOperator(nuclear_repulsion, one_body_integrals, tens_tmp)
    end
    if spin_order == :block
        block_spin = true
    else
        block_spin = false
    end
    tens1, tens2 = spin_orbital_from_spatial(one_body_integrals, tens_tmp, block_spin=block_spin, index_order=index_order)
    return InteractionOperator(nuclear_repulsion, tens1, tens2)
end

"""
    zchop!(iop::InteractionOperator, zeps=ZChop.ZEPS)

Remove terms in `iop` that are close to zero.
"""
function ZChop.zchop!(iop::InteractionOperator, zeps::Real=ZChop.ZEPS)
    ZChop.zchop!(iop.one_body_tensor, zeps)
    ZChop.zchop!(iop.two_body_tensor, zeps)
    return iop
end

function ZChop.nchop(iop::InteractionOperator)
    return InteractionOperator(ZChop.nchop(iop.nuclear_repulsion),
                               ZChop.nchop(iop.one_body_tensor),
                               ZChop.nchop(iop.two_body_tensor))
end

"""
    zchop(iop::InteractionOperator, zeps=ZChop.ZEPS)

See `zchop!(iop::InteractionOperator)`
"""
function ZChop.zchop(iop::InteractionOperator, zeps::Real=ZChop.ZEPS)
    return InteractionOperator(ZChop.zchop(iop.nuclear_repulsion, zeps),
                               ZChop.zchop(iop.one_body_tensor, zeps),
                               ZChop.zchop(iop.two_body_tensor, zeps))
end

"""
    find_index_order(iop::InteractionOperator)

Return the index-ordering convention of the two-body tensor in `iop`.
"""
find_index_order(iop::InteractionOperator) = find_index_order(iop.two_body_tensor)


text:
push!(LOAD_PATH,"../src/")
using Documenter, Forecast

makedocs(
    sitename="Forecast.jl",
    modules = [Forecast],
    pages = [
        "index.md",
        "man/quickstart.md",
        "man/methods.md",
        "man/datasets.md",
        "man/docstrings.md",
    ]
)

deploydocs(
    repo = "github.com/viraltux/Forecast.jl.git",
    devbranch = "main",
)

ext: jl


text:
module Pkg442

using Dep442A
using Dep442B

export check442

check442() = true

end # module


text:
@hydromodel Stochastic EmptyReservoirs = begin
    # First stage
    # ========================================================
    @stage 1 begin
        @parameters begin
            horizon
            indices
            data
        end
        @unpack plants = indices
        @unpack hydrodata = data
        # Variables
        # ========================================================
        @variable(model, M₀[p = plants], lowerbound = 0, upperbound = hydrodata[p].M̄)
    end
    # Second stage
    # =======================================================
    @stage 2 begin
        @parameters begin
            horizon
            indices
            data
        end
        @unpack hours, plants, segments = indices
        @unpack hydrodata = data
        @uncertain ρ from EmptyReservoirsScenario
        # Variables
        # =======================================================
        # First stage
        @decision M₀
        # -------------------------------------------------------
        @variable(model, Q[p = plants, s = segments, t = hours], lowerbound = 0, upperbound = hydrodata[p].Q̄[s])
        @variable(model, S[p = plants, t = hours] >= 0)
        @variable(model, M[p = plants, t = hours], lowerbound = 0, upperbound = hydrodata[p].M̄)
        @variable(model, H[t = hours] >= 0)
        @variable(model, Qf[p = plants, t = hours] >= 0)
        @variable(model, Sf[p = plants, t = hours] >= 0)
        # Objectives
        # ========================================================
        # Net profit
        @expression(model, net_profit,
                    sum(ρ[t]*H[t]
                        for t = hours))
        @objective(model, Max, net_profit)
        # Constraints
        # ========================================================
        # Hydrological balance
        @constraint(model, hydro_constraints[p = plants, t = hours],
                    # Previous reservoir content
                    M[p,t] == (t > 1 ? M[p,t-1] : M₀[p])
                    # Inflow
                    + sum(Qf[i,t]
                          for i = intersect(hydrodata.Qu[p],plants))
                    + sum(Sf[i,t]
                          for i = intersect(hydrodata.Su[p],plants))
                    # Outflow
                    - sum(Q[p,s,t]
                          for s = segments)
                    - S[p,t]
                    )
        # Production
        @constraint(model, production[t = hours],
                    H[t] == sum(hydrodata[p].μ[s]*Q[p,s,t]
                                for p = plants, s = segments)
                    )
        # Water flow: Discharge + Spillage
        @constraintref Qflow[1:length(plants),1:nhours(horizon)]
        @constraintref Sflow[1:length(plants),1:nhours(horizon)]
        for (pidx,p) = enumerate(plants)
            for t = hours
                if t - hydrodata[p].Rqh > 1
                    Qflow[pidx,t] = @constraint(model,
                                                Qf[p,t] == (hydrodata[p].Rqm/60)*sum(Q[p,s,t-(hydrodata[p].Rqh+1)]
                                                                                for s = segments)
                                                + (1-hydrodata[p].Rqm/60)*sum(Q[p,s,t-hydrodata[p].Rqh]
                                                                          for s = segments)
                                                )
                elseif t - hydrodata[p].Rqh > 0
                    Qflow[pidx,t] = @constraint(model,
                                                Qf[p,t] == (1-hydrodata[p].Rqm/60)*sum(Q[p,s,t-hydrodata[p].Rqh]
                                                                                for s = segments)
                                                )
                else
                    Qflow[pidx,t] = @constraint(model,
                                                Qf[p,t] == 0
                                                )
                end
                if t - hydrodata[p].Rsh > 1
                    Sflow[pidx,t] = @constraint(model,
                                                Sf[p,t] == (hydrodata[p].Rsm/60)*S[p,t-(hydrodata[p].Rsh+1)]
                                                + (1-hydrodata[p].Rsm/60)*S[p,t-hydrodata[p].Rsh]
                                                )
                elseif t - hydrodata[p].Rsh > 0
                    Sflow[pidx,t] = @constraint(model,
                                                Sf[p,t] == (1-hydrodata[p].Rsm/60)*S[p,t-hydrodata[p].Rsh]
                                                )
                else
                    Sflow[pidx,t] = @constraint(model,
                                                Sf[p,t] == 0
                                                )
                end
            end
        end
    end
end

function EmptyReservoirsModel(data::AbstractModelData, args...; kw...)
    return EmptyReservoirsModel(maximum_horizon(data), data, args...; kw...)
end


text:
#############################
# Directed acyclic graph IR #
#############################

include("dag.jl")

##############################
# generative functions types #
##############################

struct StaticIRGenerativeFunctionOptions
    track_diffs::Bool
    cache_julia_nodes::Bool
end

# trace code generation
include("trace.jl")

"""
    StaticIRGenerativeFunction{T,U} <: GenerativeFunction{T,U}

Abstact type for a static IR generative function with return type T and trace type U.

Contains an intermediate representation based on a directed acyclic graph.
Most generative function interface methods are generated from the intermediate representation.
"""
abstract type StaticIRGenerativeFunction{T,U} <: GenerativeFunction{T,U} end

function get_ir end
function get_gen_fn_type end
function get_options end

function generate_generative_function(ir::StaticIR, name::Symbol; track_diffs=false, cache_julia_nodes=true)
    options = StaticIRGenerativeFunctionOptions(track_diffs, cache_julia_nodes)
    generate_generative_function(ir, name, options)
end

function generate_generative_function(ir::StaticIR, name::Symbol, options::StaticIRGenerativeFunctionOptions)

    (trace_defns, trace_struct_name) = generate_trace_type_and_methods(ir, name, options)

    gen_fn_type_name = gensym("StaticGenFunction_$name")
    return_type = ir.return_node.typ
    trace_type = trace_struct_name
    has_argument_grads = tuple(map((node) -> node.compute_grad, ir.arg_nodes)...)
    accepts_output_grad = ir.accepts_output_grad

    gen_fn_defn = quote
        struct $gen_fn_type_name <: $(QuoteNode(StaticIRGenerativeFunction)){$return_type,$trace_type}
            params_grad::Dict{Symbol,Any}
            params::Dict{Symbol,Any}
        end
        (gen_fn::$gen_fn_type_name)(args...) = $(GlobalRef(Gen, :propose))(gen_fn, args)[3]
        $(GlobalRef(Gen, :get_ir))(::Type{$gen_fn_type_name}) = $(QuoteNode(ir))
        $(GlobalRef(Gen, :get_trace_type))(::Type{$gen_fn_type_name}) = $trace_struct_name
        $(GlobalRef(Gen, :has_argument_grads))(::$gen_fn_type_name) = $(QuoteNode(has_argument_grads))
        $(GlobalRef(Gen, :accepts_output_grad))(::$gen_fn_type_name) = $(QuoteNode(accepts_output_grad))
        $(GlobalRef(Gen, :get_gen_fn))(trace::$trace_struct_name) = $(Expr(:(.), :trace, QuoteNode(static_ir_gen_fn_ref)))
        $(GlobalRef(Gen, :get_gen_fn_type))(::Type{$trace_struct_name}) = $gen_fn_type_name
        $(GlobalRef(Gen, :get_options))(::Type{$gen_fn_type_name}) = $(QuoteNode(options))
    end
    Expr(:block, trace_defns, gen_fn_defn, Expr(:call, gen_fn_type_name, :(Dict{Symbol,Any}()), :(Dict{Symbol,Any}())))
end

include("render_ir.jl")

###########################
# generative function API #
###########################

# variable names used in generated code
const trace = gensym("trace")
const weight = gensym("weight")
const subtrace = gensym("subtrace")

include("simulate.jl")
include("generate.jl")
include("project.jl")
include("update.jl")
include("backprop.jl")


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

mutable struct IoK8sApiAdmissionregistrationV1MutatingWebhook <: SwaggerModel
    admissionReviewVersions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: admissionReviewVersions
    clientConfig::Any # spec type: Union{ Nothing, IoK8sApiAdmissionregistrationV1WebhookClientConfig } # spec name: clientConfig
    failurePolicy::Any # spec type: Union{ Nothing, String } # spec name: failurePolicy
    matchPolicy::Any # spec type: Union{ Nothing, String } # spec name: matchPolicy
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    namespaceSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: namespaceSelector
    objectSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: objectSelector
    reinvocationPolicy::Any # spec type: Union{ Nothing, String } # spec name: reinvocationPolicy
    rules::Any # spec type: Union{ Nothing, Vector{IoK8sApiAdmissionregistrationV1RuleWithOperations} } # spec name: rules
    sideEffects::Any # spec type: Union{ Nothing, String } # spec name: sideEffects
    timeoutSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: timeoutSeconds

    function IoK8sApiAdmissionregistrationV1MutatingWebhook(;admissionReviewVersions=nothing, clientConfig=nothing, failurePolicy=nothing, matchPolicy=nothing, name=nothing, namespaceSelector=nothing, objectSelector=nothing, reinvocationPolicy=nothing, rules=nothing, sideEffects=nothing, timeoutSeconds=nothing)
        o = new()
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("admissionReviewVersions"), admissionReviewVersions)
        setfield!(o, Symbol("admissionReviewVersions"), admissionReviewVersions)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("clientConfig"), clientConfig)
        setfield!(o, Symbol("clientConfig"), clientConfig)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("failurePolicy"), failurePolicy)
        setfield!(o, Symbol("failurePolicy"), failurePolicy)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("matchPolicy"), matchPolicy)
        setfield!(o, Symbol("matchPolicy"), matchPolicy)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("namespaceSelector"), namespaceSelector)
        setfield!(o, Symbol("namespaceSelector"), namespaceSelector)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("objectSelector"), objectSelector)
        setfield!(o, Symbol("objectSelector"), objectSelector)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("reinvocationPolicy"), reinvocationPolicy)
        setfield!(o, Symbol("reinvocationPolicy"), reinvocationPolicy)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("rules"), rules)
        setfield!(o, Symbol("rules"), rules)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("sideEffects"), sideEffects)
        setfield!(o, Symbol("sideEffects"), sideEffects)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("timeoutSeconds"), timeoutSeconds)
        setfield!(o, Symbol("timeoutSeconds"), timeoutSeconds)
        o
    end
end # type IoK8sApiAdmissionregistrationV1MutatingWebhook

const _property_map_IoK8sApiAdmissionregistrationV1MutatingWebhook = Dict{Symbol,Symbol}(Symbol("admissionReviewVersions")=>Symbol("admissionReviewVersions"), Symbol("clientConfig")=>Symbol("clientConfig"), Symbol("failurePolicy")=>Symbol("failurePolicy"), Symbol("matchPolicy")=>Symbol("matchPolicy"), Symbol("name")=>Symbol("name"), Symbol("namespaceSelector")=>Symbol("namespaceSelector"), Symbol("objectSelector")=>Symbol("objectSelector"), Symbol("reinvocationPolicy")=>Symbol("reinvocationPolicy"), Symbol("rules")=>Symbol("rules"), Symbol("sideEffects")=>Symbol("sideEffects"), Symbol("timeoutSeconds")=>Symbol("timeoutSeconds"))
const _property_types_IoK8sApiAdmissionregistrationV1MutatingWebhook = Dict{Symbol,String}(Symbol("admissionReviewVersions")=>"Vector{String}", Symbol("clientConfig")=>"IoK8sApiAdmissionregistrationV1WebhookClientConfig", Symbol("failurePolicy")=>"String", Symbol("matchPolicy")=>"String", Symbol("name")=>"String", Symbol("namespaceSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("objectSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("reinvocationPolicy")=>"String", Symbol("rules")=>"Vector{IoK8sApiAdmissionregistrationV1RuleWithOperations}", Symbol("sideEffects")=>"String", Symbol("timeoutSeconds")=>"Int32")
Base.propertynames(::Type{ IoK8sApiAdmissionregistrationV1MutatingWebhook }) = collect(keys(_property_map_IoK8sApiAdmissionregistrationV1MutatingWebhook))
Swagger.property_type(::Type{ IoK8sApiAdmissionregistrationV1MutatingWebhook }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAdmissionregistrationV1MutatingWebhook[name]))}
Swagger.field_name(::Type{ IoK8sApiAdmissionregistrationV1MutatingWebhook }, property_name::Symbol) =  _property_map_IoK8sApiAdmissionregistrationV1MutatingWebhook[property_name]

function check_required(o::IoK8sApiAdmissionregistrationV1MutatingWebhook)
    (getproperty(o, Symbol("admissionReviewVersions")) === nothing) && (return false)
    (getproperty(o, Symbol("clientConfig")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("sideEffects")) === nothing) && (return false)
    true
end

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


text:
using DataStructures

lines = readlines("inputs/11-15/15.txt")
# lines = readlines("inputs/11-15/15ex.txt")

M = zeros(Int, length(lines),length(lines[1]))
for (i,l) in enumerate(lines)
  for (j,c) in enumerate(split(l,""))
    # @show c
    M[i,j] = parse(Int,String(c))
  end
end

# is (x,y) a valid coord in matrix A?
inrange((x,y),A) = x in 1:size(A,1) && y in 1:size(A,2)

q = [(0,1,1)]
costs = Dict()
while true
  global q
  cost,x,y = q[1]
  if x==size(M)[1] && y==size(M)[2]
    @show cost
    break
  end
  popfirst!(q)
  for (dx,dy) ∈ [(-1,0), (1,0), (0,1), (0,-1)]
    if inrange((x+dx,y+dy),M)
      ncost = cost + M[x+dx,y+dy]
      if (x+dx,y+dy) ∈ keys(costs) && costs[(x+dx,y+dy)] <= ncost
        continue
      end
      costs[(x+dx,y+dy)] = ncost
      push!(q,(ncost,x+dx,y+dy))
    end
  end
  q = sort(q)
end


text:
# Tests
# #+HTML:  Source
# #+HTML:

# [[file:~/.julia/dev/MatrixProductStates/README.org::*Tests][Tests:1]]
using Test, MatrixProductStates, SparseArrays, Arpack

@testset "TFIM   " begin
    g = 1.0; L = 7

    function H_TFIM(g, L)
        id = [1  0;
              0  1]
        σˣ = [0  1;
              1  0]
        σᶻ = [1  0;
              0 -1]
        W_tnsr = zeros(Complex{Float64}, 3, 3, 2, 2)
        W_tnsr[1, 1, :, :] = id
        W_tnsr[2, 1, :, :] = -σᶻ
        W_tnsr[3, 1, :, :] = -g*σˣ
        W_tnsr[3, 2, :, :] = σᶻ
        W_tnsr[3, 3, :, :] = id

        return MPO(W_tnsr, L)
    end
    H = H_TFIM(g, L)
    ψ = randn(MPS{L, Complex{Float64}}, 100, 2)

    ψ̃ = compress(ψ, left, Dcut=80)[1] # Note: no actual information is lost in this
    # compression because of the small size of the chain

    @test              ψ̃'ψ̃ ≈ 1
    @test          ψ'ψ/ψ'ψ ≈ ψ̃'ψ̃
    @test ((ψ'*(H*ψ))/ψ'ψ) ≈ (ψ̃' * (H * ψ̃))/ψ̃'ψ̃
    @test ((ψ'*(H*ψ))/ψ'ψ) ≈ (ψ̃' * (H * ψ))/ψ̃'ψ

    ϕ, E₀ = ground_state(ψ, H, quiet=true)
    @test ϕ' * H * H * ϕ ≈ (ϕ'*H*ϕ)^2
end

@testset "Hubbard" begin

    id = [1 0
          0 1]
    c  = [0 0
          1 0] #Anti commuting matrix
    c_up = c  ⊗ id
    c_dn = id ⊗ c
    id²  = id ⊗ id
    n_up = c_up' * c_up
    n_dn = c_dn' * c_dn

    P_up = (id² - 2c_up'*c_up) # Spin up parity operator
    P_dn = (id² - 2c_dn'*c_dn) # Spin down parity operator

    function H_hub(U, μ, L)
        W_tnsr = zeros(Complex{Float64}, 6, 6, 4, 4)
        W_tnsr[1, 1, :, :] = id²
        W_tnsr[2, 1, :, :] = c_up'
        W_tnsr[3, 1, :, :] = c_dn'
        W_tnsr[4, 1, :, :] = c_up
        W_tnsr[5, 1, :, :] = c_dn
        W_tnsr[6, 1, :, :] = U*(n_up * n_dn) - μ*(n_up + n_dn)
        W_tnsr[6, 2, :, :] =  c_up  * P_up  # Must multiply by the parity operator to get
        W_tnsr[6, 3, :, :] =  c_dn  * P_dn  # correct off-site commutation relations!
        W_tnsr[6, 4, :, :] = -c_up' * P_up
        W_tnsr[6, 5, :, :] = -c_dn' * P_dn
        W_tnsr[6, 6, :, :] = id²
        MPO(W_tnsr, L)
    end

    function solve_hub(U, μ, L; retfull=true, quiet=true)
        H = H_hub(U, μ, L)
        ψ = randn(MPS{L, Complex{Float64}}, 100, 4)
        (ϕ, E₀), t, bytes = @timed ground_state(ψ, H, ϵ=1e-5, quiet=quiet)

        (ϕ=ϕ, E₀=E₀, H=H, t=t, Gbytes=bytes/1e9)
    end

    function Hub_ED(U, μ, L,)
        Û = U*(n_up * n_dn) - μ*(n_up + n_dn)
        c_dg_up(i) = foldl(⊗, sparse.([i==j ? c_up' : id² for j in 1:L]))
        cup(i)     = foldl(⊗, sparse.([i==j ? c_up  : id² for j in 1:L]))
        c_dg_dn(i) = foldl(⊗, sparse.([i==j ? c_dn' : id² for j in 1:L]))
        cdn(i)     = foldl(⊗, sparse.([i==j ? c_dn  : id² for j in 1:L]))
        Ûf(i)      = foldl(⊗, sparse.([i==j ? Û     : id² for j in 1:L]))
        function c_dg_c(i)
            out = c_dg_up(i)*cup(i+1) + c_dg_dn(i)*cdn(i+1)
            out + out'
        end
        H = -sum(c_dg_c, 1:(L-1)) + sum(Ûf, 1:L)

        λ, ϕ = eigs(H, nev=1, which=:SR)
        (ϕ'H*ϕ)[]
    end

    U = 3.0; μ = -1.0; L = 4
    H = H_hub(U, μ, L)

    ϕ, E₀ = solve_hub(U, μ, L, retfull=true, quiet=true)
    @test ϕ' * H * H * ϕ ≈ (ϕ'*H*ϕ)^2  # Make sure energy is eigenvalue
    @test ϕ' * H * ϕ ≈ E₀              # make sure eigenvalue matches one produced by alogrithm
    @test ϕ' * H * ϕ ≈ Hub_ED(U, μ, L) # check against exact diagonalization
end
# Tests:1 ends here


text:
using Distributed
addprocs(24; exeflags = "--project")

println("Loading packages")
@everywhere using Pkg
@everywhere Pkg.instantiate()
@everywhere using gcSolver

println("Starting fitting")
gcSolver.runFit()


text:
function hargreaves(
    forcings::DataFrame;
    tminCol::Symbol=:tmin,
    tmaxCol::Symbol=:tmax,
    dateCol::Symbol=:datetime
)

    dts = forcings[!,dateCol]
    tmin = forcings[!,tminCol]
    tmax = forcings[!,tmaxCol]
    len = length(tmax)
    Gsc = 367
    lhov = 2.257

    doy = map(dayofyear, dts)

    tavg = map(mean,zip(tmin,tmax))

    eto = zeros(len)

    for (i,t) in enumerate(doy)
        b = 2 * pi * (Int16(t)/365)
        Rav = 1.00011 + 0.034221*cos(b) + 0.00128*sin(b) + 0.000719*cos(2*b) + 0.000077*sin(2*b)
        Ho = ((Gsc * Rav) * 86400)/1e6

        eto[i] = (0.0023 * Ho * (tmax[i]-tmin[i])^0.5 * (tavg[i]+17.8))
    end

    return eto
end

function hamon(forcings::DataFrame, par::Float64, lat::Float64; tavgCol::Symbol=:tavg, dateCol::Symbol=:dataetime)

    dts = forcings[!,dateCol]
    tavg = forcings[!,tavgCol]
    len = length(tavg)

    doy = map(dayofyear, dts)

    θ = @. 0.2163108 + 2 * atan(0.9671396 * tan(0.0086 * (doy - 186)))
	pi_v = @. asin(0.39795 * cos(θ))
    daylighthr = @. 24 - 24/pi * acos((sin(0.8333 * pi/180) + sin(lat * pi/180) * sin(pi_v))/(cos(lat * pi/180) * cos(pi_v)))

    esat = @. 0.611 * exp(17.27 * tavg/(237.3 + tavg))

    eto = @. par * 29.8 * daylighthr * (esat/(tavg + 273.2))

  return eto
end


text:
using Compose
import Cairo, Fontconfig

imgs = [SVG("transformations.svg", 7cm, 7cm),
        PDF("transformations.pdf", 7cm, 7cm)]

img = compose(context(),
    (context(0,0,0.5,0.5), xgon(0.2,0.2,0.1,3,0.3),
    (context(mirror=Mirror(-π/4, 0.4,0.4)), xgon(0.2,0.2,0.1,3,0.3)) ),
    (context(0.5,0,0.5,0.5), xgon(0.5,0.5,0.2,3,0.3), fill("silver"),
    (context(shear=Shear(0.8, 0, 0.5,0.5)), xgon(0.5,0.5,0.2,3,0.3), fill(nothing), stroke("black")) )
)

draw.(imgs, [img])


text:
int_rules_4_7_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.7.2*Tangent*normalization*rules =#
    @apply_utils Antiderivative(~u * (~(c') * cot(~(a') + ~(b') * ~x)) ^ ~(m') * (~(d') * tan(~(a') + ~(b') * ~x)) ^ ~(n'), ~x) => (~c * cot(~a + ~b * ~x)) ^ ~m * (~d * tan(~a + ~b * ~x)) ^ ~m * Antiderivative(ActivateTrig(~u) * (~d * tan(~a + ~b * ~x)) ^ (~n - ~m), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && KnownTangentIntegrandQ(~u, ~x)
    @apply_utils Antiderivative(~u * (~(c') * tan(~(a') + ~(b') * ~x)) ^ ~(m') * (~(d') * cos(~(a') + ~(b') * ~x)) ^ ~(n'), ~x) => (((~c * tan(~a + ~b * ~x)) ^ ~m * (~d * cos(~a + ~b * ~x)) ^ ~m) / (~d * sin(~a + ~b * ~x)) ^ ~m) * Antiderivative((ActivateTrig(~u) * (~d * sin(~a + ~b * ~x)) ^ ~m) / (~d * cos(~a + ~b * ~x)) ^ (~m - ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && KnownCotangentIntegrandQ(~u, ~x)
    @apply_utils Antiderivative(~u * (~(c') * cot(~(a') + ~(b') * ~x)) ^ ~(m'), ~x) => (~c * cot(~a + ~b * ~x)) ^ ~m * (~c * tan(~a + ~b * ~x)) ^ ~m * Antiderivative(ActivateTrig(~u) / (~c * tan(~a + ~b * ~x)) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~m], ~x) && (Not(IntegerQ(~m)) && KnownTangentIntegrandQ(~u, ~x))
    @apply_utils Antiderivative(~u * (~(c') * tan(~(a') + ~(b') * ~x)) ^ ~(m'), ~x) => (~c * cot(~a + ~b * ~x)) ^ ~m * (~c * tan(~a + ~b * ~x)) ^ ~m * Antiderivative(ActivateTrig(~u) / (~c * cot(~a + ~b * ~x)) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~m], ~x) && (Not(IntegerQ(~m)) && KnownCotangentIntegrandQ(~u, ~x))
    @apply_utils Antiderivative(~u * (~(c') * tan(~(a') + ~(b') * ~x)) ^ ~(n') * (~A + ~(B') * cot(~(a') + ~(b') * ~x)), ~x) => ~c * Antiderivative(ActivateTrig(~u) * (~c * tan(~a + ~b * ~x)) ^ (~n - 1) * (~B + ~A * tan(~a + ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~A, ~B, ~n], ~x) && KnownTangentIntegrandQ(~u, ~x)
    @apply_utils Antiderivative(~u * (~(c') * cot(~(a') + ~(b') * ~x)) ^ ~(n') * (~A + ~(B') * tan(~(a') + ~(b') * ~x)), ~x) => ~c * Antiderivative(ActivateTrig(~u) * (~c * cot(~a + ~b * ~x)) ^ (~n - 1) * (~B + ~A * cot(~a + ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~A, ~B, ~n], ~x) && KnownCotangentIntegrandQ(~u, ~x)
    @apply_utils Antiderivative(~u * (~A + ~(B') * cot(~(a') + ~(b') * ~x)), ~x) => Antiderivative((ActivateTrig(~u) * (~B + ~A * tan(~a + ~b * ~x))) / tan(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~A, ~B], ~x) && KnownTangentIntegrandQ(~u, ~x)
    @apply_utils Antiderivative(~u * (~A + ~(B') * tan(~(a') + ~(b') * ~x)), ~x) => Antiderivative((ActivateTrig(~u) * (~B + ~A * cot(~a + ~b * ~x))) / cot(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~A, ~B], ~x) && KnownCotangentIntegrandQ(~u, ~x)
    @apply_utils Antiderivative(~(u') * (~(c') * tan(~(a') + ~(b') * ~x)) ^ ~(n') * (~(A') + ~(B') * cot(~(a') + ~(b') * ~x) + ~(C') * cot(~(a') + ~(b') * ~x) ^ 2), ~x) => (~c) ^ 2 * Antiderivative(ActivateTrig(~u) * (~c * tan(~a + ~b * ~x)) ^ (~n - 2) * (~C + ~B * tan(~a + ~b * ~x) + ~A * tan(~a + ~b * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~A, ~B, ~C, ~n], ~x) && KnownTangentIntegrandQ(~u, ~x)
    @apply_utils Antiderivative(~(u') * (~(c') * cot(~(a') + ~(b') * ~x)) ^ ~(n') * (~(A') + ~(B') * tan(~(a') + ~(b') * ~x) + ~(C') * tan(~(a') + ~(b') * ~x) ^ 2), ~x) => (~c) ^ 2 * Antiderivative(ActivateTrig(~u) * (~c * cot(~a + ~b * ~x)) ^ (~n - 2) * (~C + ~B * cot(~a + ~b * ~x) + ~A * cot(~a + ~b * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~A, ~B, ~C, ~n], ~x) && KnownCotangentIntegrandQ(~u, ~x)
    @apply_utils Antiderivative(~(u') * (~(c') * tan(~(a') + ~(b') * ~x)) ^ ~(n') * (~A + ~(C') * cot(~(a') + ~(b') * ~x) ^ 2), ~x) => (~c) ^ 2 * Antiderivative(ActivateTrig(~u) * (~c * tan(~a + ~b * ~x)) ^ (~n - 2) * (~C + ~A * tan(~a + ~b * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~A, ~C, ~n], ~x) && KnownTangentIntegrandQ(~u, ~x)
    @apply_utils Antiderivative(~(u') * (~(c') * cot(~(a') + ~(b') * ~x)) ^ ~(n') * (~A + ~(C') * tan(~(a') + ~(b') * ~x) ^ 2), ~x) => (~c) ^ 2 * Antiderivative(ActivateTrig(~u) * (~c * cot(~a + ~b * ~x)) ^ (~n - 2) * (~C + ~A * cot(~a + ~b * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~A, ~C, ~n], ~x) && KnownCotangentIntegrandQ(~u, ~x)
    @apply_utils Antiderivative(~u * (~(A') + ~(B') * cot(~(a') + ~(b') * ~x) + ~(C') * cot(~(a') + ~(b') * ~x) ^ 2), ~x) => Antiderivative((ActivateTrig(~u) * (~C + ~B * tan(~a + ~b * ~x) + ~A * tan(~a + ~b * ~x) ^ 2)) / tan(~a + ~b * ~x) ^ 2, ~x) <-- FreeQ([~a, ~b, ~A, ~B, ~C], ~x) && KnownTangentIntegrandQ(~u, ~x)
    @apply_utils Antiderivative(~u * (~(A') + ~(B') * tan(~(a') + ~(b') * ~x) + ~(C') * tan(~(a') + ~(b') * ~x) ^ 2), ~x) => Antiderivative((ActivateTrig(~u) * (~C + ~B * cot(~a + ~b * ~x) + ~A * cot(~a + ~b * ~x) ^ 2)) / cot(~a + ~b * ~x) ^ 2, ~x) <-- FreeQ([~a, ~b, ~A, ~B, ~C], ~x) && KnownCotangentIntegrandQ(~u, ~x)
    @apply_utils Antiderivative(~u * (~A + ~(C') * cot(~(a') + ~(b') * ~x) ^ 2), ~x) => Antiderivative((ActivateTrig(~u) * (~C + ~A * tan(~a + ~b * ~x) ^ 2)) / tan(~a + ~b * ~x) ^ 2, ~x) <-- FreeQ([~a, ~b, ~A, ~C], ~x) && KnownTangentIntegrandQ(~u, ~x)
    @apply_utils Antiderivative(~u * (~A + ~(C') * tan(~(a') + ~(b') * ~x) ^ 2), ~x) => Antiderivative((ActivateTrig(~u) * (~C + ~A * cot(~a + ~b * ~x) ^ 2)) / cot(~a + ~b * ~x) ^ 2, ~x) <-- FreeQ([~a, ~b, ~A, ~C], ~x) && KnownCotangentIntegrandQ(~u, ~x)
    @apply_utils Antiderivative(~u * (~(A') + ~(B') * tan(~(a') + ~(b') * ~x) + ~(C') * cot(~(a') + ~(b') * ~x)), ~x) => Antiderivative((ActivateTrig(~u) * (~C + ~A * tan(~a + ~b * ~x) + ~B * tan(~a + ~b * ~x) ^ 2)) / tan(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~A, ~B, ~C], ~x)
    @apply_utils Antiderivative(~u * (~(A') * tan(~(a') + ~(b') * ~x) ^ ~(n') + ~(B') * tan(~(a') + ~(b') * ~x) ^ ~n1 + ~(C') * tan(~(a') + ~(b') * ~x) ^ ~n2), ~x) => Antiderivative(ActivateTrig(~u) * tan(~a + ~b * ~x) ^ ~n * (~A + ~B * tan(~a + ~b * ~x) + ~C * tan(~a + ~b * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~A, ~B, ~C, ~n], ~x) && (EqQ(~n1, ~n + 1) && EqQ(~n2, ~n + 2))
    @apply_utils Antiderivative(~u * (~(A') * cot(~(a') + ~(b') * ~x) ^ ~(n') + ~(B') * cot(~(a') + ~(b') * ~x) ^ ~n1 + ~(C') * cot(~(a') + ~(b') * ~x) ^ ~n2), ~x) => Antiderivative(ActivateTrig(~u) * cot(~a + ~b * ~x) ^ ~n * (~A + ~B * cot(~a + ~b * ~x) + ~C * cot(~a + ~b * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~A, ~B, ~C, ~n], ~x) && (EqQ(~n1, ~n + 1) && EqQ(~n2, ~n + 2))
end

max_stars_repo_path: src/rules/4 Trig functions/4.7 Miscellaneous/4.7.2 Tangent normalization rules.jl


text:
export blockdiag, contract, cp_als, diagt, dirsum, dropdims, hosvd, innerprod, krontm, matten, mkrontm, mkrontv, mrank, mttkrp, neye, nrank, nvecs
export sthosvd, tenmat, tkron, ttm, ttt, ttv

"""
   blockdiag(X,Y)

Create block diagonal tensor where tensors X and Y are block elements. If X and Y are matrices, equal to blkdiag for sparse matrices.
"""
function blockdiag(X1::Array{<:Number,N},X2::Array{<:Number,N}) where N
  sz=tuple([size(X1)...]+[size(X2)...]...)
  Xd=zeros(sz)
  R1=CartesianIndices(size(X1))
  I1=last(R1)
  R2=CartesianIndices(size(X2))
  for In in R1
    Xd[In]=X1[In]
  end
  for In in R2
    Xd[I1+In]=X2[In]
  end
  Xd
end
"""
   dirsum(X,Y)

Direct sum od tensors. Equal to blockdiag.
"""
dirsum(X1::AbstractArray{<:Number},X2::AbstractArray{<:Number}) = blockdiag(X1,X2)

"""
    contract(X,Y)
    contract(X,indX,Y,indY[,perm])
    contract(X::TensorCell)

Contracted product of tensors. Contract indX modes of array X to indY modes of array Y and permute the result by vector perm.
Default: indX=[ndims(X)], indY=[1].
"""
function contract(X1::AbstractArray{<:Number},ind1::Vector{Int},X2::AbstractArray{<:Number},ind2::AbstractVector{<:Integer},perm=[])
	if length(X1) == 0
		return X2
	elseif length(X2)==0
		return X1
	end
    sz1=[size(X1)...]
    sz2=[size(X2)...]
    @assert(sz1[ind1]==sz2[ind2],"Dimension mismatch.")
    sz=[sz1[setdiff(1:ndims(X1),ind1)];sz2[setdiff(1:ndims(X2),ind2)]]
    #Xres=reshape(transpose(tenmat(X1,row=ind1))*tenmat(X2,row=ind2),tuple(sz...))
	Xres=reshape(transpose(tenmat(X1,row=ind1))*tenmat(X2,row=ind2),tuple(sz...))
    if perm!=[]
      Xres=permutedims(Xres,perm)
    end
    Xres
end
contract(X1::AbstractArray{<:Number},X2::AbstractArray{<:Number},ind::AbstractVector{<:Integer})=contract(X1,ind,X2,ind)
contract(X1::AbstractArray{<:Number},ind1::Int,X2::AbstractArray{<:Number},ind2::Int,perm=[]) =contract(X1,[ind1],X2,[ind2],perm)
contract(X1::AbstractArray{<:Number},X2::AbstractArray{<:Number}) = contract(X1,ndims(X1),X2,1)
function contract(X1::AbstractArray{<:Number},X2::AbstractArray{<:Number},Xk...)
	Xres=contract(X1,X2)
	for k=1:length(Xk)
		Xres=contract(Xres,Xk[k])
	end
	Xres
end
function contract(X::Array{AbstractArray{<:Number,D} where D,1},squeeze=true)
    N=length(X)
    if N==1
        return X[1]
    end
    @assert(any([size(X[n])[end]==size(X[n+1])[1] for n=1:N-1]),"Dimensions mismatch.")
    Xcontr=contract(X[1],X[2])
    for n=3:N
        Xcontr=contract(Xcontr,X[n])
    end
    squeeze == true ? dropdims(Xcontr) : Xcontr
end
function contract(X::Array{AbstractArray{<:Number,D},1},squeeze=true) where D
	N=length(X)
	if N==1
		return X[1]
	end
	@assert(any([size(X[n])[end]==size(X[n+1])[1] for n=1:N-1]),"Dimensions mismatch.")
	Xcontr=contract(X[1],X[2])
	for n=3:N
		Xcontr=contract(Xcontr,X[n])
	end
	squeeze == true ? dropdims(Xcontr) : Xcontr
end

"""
    cp_als(X,R;init,tol,maxit,dimorder)

Compute a CP decomposition with R components of a tensor X .

## Arguments:
- `init` ∈ {MatrixCell,"rand","nvecs","eigs"}. Initial guess for factor matrices. If init="nvecs" (same as "eigs") initialize matrices with function nvecs.
- `tol`: Tolerance. Defualt: 1e-4.
- `maxit`: Maximal number of iterations. Default: 1000.
- `dimorder': Order of dimensions. Default: 1:ndims(A).
"""
function cp_als(X::AbstractArray{<:Number},R::Integer;init="rand",tol=1e-4,maxit=1000,dimorder=[])
    N=ndims(X)
    nr=norm(X)
    K=ktensor
    if length(dimorder) == 0
        dimorder=collect(1:N)
    end
    fmat=MatrixCell(undef,N)
    if isa(init,Vector) || isa(init,MatrixCell)
        @assert(length(init)==N,"Wrong number of initial matrices.")
        for n in dimorder[2:end]
            @assert(size(init[n])==(size(X,n),R),"$(n)-th initial matrix is of wrong size.")
            fmat[n]=init[n]
        end
    elseif init=="rand"
        [fmat[n]=rand(size(X,n),R) for n in dimorder[2:end]]
    elseif init=="eigs" || init=="nvecs"
        [fmat[n]=nvecs(X,n,R) for n in dimorder[2:end]]
    else
        error("Initialization method wrong.")
    end
    G = zeros(R,R,N); #initalize gramians
    for n in dimorder[2:end]
      if !isempty(fmat[n])
        G[:,:,n]=fmat[n]'*fmat[n]
      end
    end
    fit=0
    for k=1:maxit
        fitold=fit
        lambda=[]
        for n in dimorder
            fmat[n]=mttkrp(X,fmat,n)
            W=reshape(prod(G[:,:,setdiff(collect(1:N),n)],dims=3),Val(2))
            fmat[n]=fmat[n]/W
            if k == 1
                lambda = sqrt.(sum(fmat[n].^2,dims=1))[:] #2-norm
            else
                lambda = maximum(maximum(abs.(fmat[n]),dims=1),dims=1)[:] #max-norm
            end
            fmat[n] = fmat[n]./lambda'
            G[:,:,n] = fmat[n]'*fmat[n]
        end
        K=ktensor(lambda,fmat)
        if nr==0
            fit=norm(K)^2-2*innerprod(X,K)
        else
            nr_res=sqrt.(abs.(nr^2+norm(K)^2 .-2*innerprod(X,K)))
            fir=1 .-nr_res/nr
        end
        fitchange=abs.(fitold-fit)
        if k>1 && fitchange<tol
            break
        end
    end
    arrange!(K)
    fixsigns!(K)
    K
end

"""
    diagt(v[,dims])

Create a diagonal tensor for a given vector of diagonal elements. Generalization of diagm.
"""
function diagt(v::AbstractVector{<:Number})
  N=length(v)
  sz=tuple(repeat([N],N,1)[:]...)
  D=zeros(sz)
  R=CartesianIndices(sz)
  In=first(R)
  for i=1:N
    D[In]=v[i]
    In=In+1*one(In)
  end
  D
end

function diagt(v::AbstractVector{<:Number},dims::AbstractVector{:Integer})
  Dt=diagt(v)
  sz=dims-[size(Dt)...]
  blockdiag(Dt,zeros(sz...))
end

"""
    hosvd(X; <keyword arguments>)

Higher-order singular value decomposition.

## Arguments:
- `X`: Tensor (multidimensional array) or ttensor.
- `method` ∈ {"svd","lanczos","randsvd"} Method for SVD. Default: "svd".
- `reqrank::Vector`: Requested mutlilinear rank. Optional.
- `eps_abs::Number/Vector`: Drop singular values (of mode-n matricization) below eps_abs. Default: 1e-8.
- `eps_rel::Number/Vector`: Drop singular values (of mode-n matricization) below eps_rel*sigma_1. Optional.
- `p::Integer`: Oversampling parameter. Defaul p=10.
"""
function hosvd(X::AbstractArray{<:Number,N};method="svd",reqrank=[],eps_abs=[],eps_rel=[],p=10) where N
  fmat=MatrixCell(undef,N)

  reqrank=check_vector_input(reqrank,N,0)
  eps_abs=check_vector_input(eps_abs,N,1e-8)
  eps_rel=check_vector_input(eps_rel,N,0)

  for n=1:N
    Xn=tenmat(X,n)
    if method == "lanczos"
      fmat[n],S=lanczos(Xn,tol=eps_abs[n],reqrank=reqrank[n],p=p)
    elseif method == "randsvd"
      fmat[n],S=randsvd(Xn,tol=eps_abs[n],reqrank=reqrank[n],p=p)
    else
      fmat[n],S,V=svd(Xn)
    end
    if reqrank[n]!=0 && size(fmat[n],2)>reqrank[n]
      fmat[n]=fmat[n][:,1:reqrank[n]]
    else
      eps_rel[n] != 0 ? tol=eps_rel[n]*S[1] : tol=eps_abs[n]
      K=findall(x-> x>tol ? true : false,S)
      fmat[n]=fmat[n][:,K]
    end
  end
  ttensor(ttm(X,fmat,'t'),fmat)
end

"""
   innerprod(X,Y)
   innerprod(X::ttensor,Y::ttensor)
   innerprod(X::ktensor,Y::ktensor)

Inner product of two tensors.
"""
function innerprod(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N}) where N
	@assert(size(X1) == size(X2),"Dimension mismatch")
	sum(X1.*conj(X2))
end

"""
   krontm(X,Y,M[,modes,t='n'])

Kronecker product of two tensors times matrix (n-mode product): (X ⊗ Y) x₁ M₁ x₂ M₂ x₃ ⋯ xₙ Mₙ.

## Arguments:
- `X::Array`
- `Y::Array`
- `M::Matrix/MatrixCell`
- `modes::Integer/Vector` : Modes for multiplication. Default: 1:length(M).
- `t='t'`: Transpose matrices from M.
"""
function krontm(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::MatrixCell,modes::AbstractVector{<:Integer},t='n') where N
  if t=='t'
    [M[n]=M[n]' for n=1:length(M)]
	end
	@assert(length(modes)<=length(M)<=N,"Dimension mismatch.")
  sz=[size(X1)...].*[size(X2)...]
  R=copy(sz)
  if length(modes) != length(M)
    M=M[modes] #discard matrices not needed for multiplication
  end
  for n=1:length(modes)
      R[modes[n]]=size(M[n],1) #vector of rₖ
  end
  #Order of multiplication - if tkron(X₁,X₂) is i₁ × i₂ × ... × iₙ and Mₖ are rₖ × iₖ, sort by largest possible dimension reduction iₖ-rₖ
  p=sortperm(sz[modes]-R[modes],rev=true)
  M=M[p]
  modes=modes[p]
  @assert(sz[modes[1]] == size(M[1],2),"Dimensions mismatch")
  Xn=copy(mkrontv(X1,X2,copy(M[1]'),modes[1],'t')')
  sz[modes[1]]=size(M[1],1)
  X=matten(Xn,modes[1],sz)
  for n=2:length(M)
	   @assert(sz[modes[n]] == size(M[n],2),"Dimensions mismatch")
     Xn=tenmat(X,modes[n])
	   sz[modes[n]]=size(M[n],1)
	   X=matten(M[n]*Xn,modes[n],sz)
	end
#  end
  X
end
krontm(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::AbstractMatrix{<:Number},n::Integer,t='n') where N=krontm(X1,X2,MatrixCell([M]),[n],t)
krontm(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::MatrixCell,t::Char) where N = krontm(X1,X2,M,1:length(M),t)
krontm(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::MatrixCell) where N = krontm(X1,X2,M,1:length(M))
krontm(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::MatrixCell,modes::AbstractRange{<:Integer},t::Char) where N = krontm(X1,X2,M,collect(modes),t)
krontm(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::MatrixCell,modes::AbstractRange{<:Integer}) where N = krontm(X1,X2,M,collect(modes))
function krontm(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::MatrixCell,n::Integer,t='n') where N
 	if n>0
 		krontm(X1,X2,M[n],n,t)
 	else
 		modes=setdiff(1:N,-n)
 		krontm(X1,X2,M,modes,t)
 	end
end
#If array of matrices isn't defined as MatrixCell, but as M=[M1,M2,...,Mn]:
krontm(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::Array{Matrix{T},1},modes::AbstractVector{<:Integer},t='n') where {T<:Number,N}=krontm(X1,X2,MatrixCell(M),modes,t)
krontm(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::Array{Matrix{T},1},t::Char) where {T<:Number,N}=krontm(X1,X2,MatrixCell(M),t)
krontm(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::Array{Matrix{T},1}) where {T<:Number,N}=krontm(X1,X2,MatrixCell(M))
krontm(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::Array{Matrix{T},1},modes::AbstractRange{<:Integer},t::Char) where {T<:Number,N}=krontm(X1,X2,MatrixCell(M),modes,t)
krontm(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::Array{Matrix{T},1},modes::AbstractRange{<:Integer}) where {T<:Number,N}=krontm(X1,X2,MatrixCell(M),modes)
krontm(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::Array{Matrix{T},1},n::Integer,t='n') where {T<:Number,N}=krontm(X1,X2,MatrixCell(M),n,t)

"""
    matten(A,n,dims)
    matten(A,R,C,dims)

Fold matrix A into a tensor of dimension dims by mode n or by row and column vectors R and C.
"""
function matten(A::AbstractMatrix{<:Number},n::Integer,dims::AbstractVector{<:Integer})
	@assert(dims[n]==size(A,1),"Dimensions mismatch")
	m = setdiff(1:length(dims), n)
	@assert prod(dims[m])==size(A,2)
	X = reshape(A,[dims[n];dims[m]]...)
	permutedims(X,invperm([n;m]))
end

function matten(A::AbstractMatrix{<:Number},row::Vector{<:Integer},col::AbstractVector{<:Integer},dims::Vector{<:Integer})
	@assert(prod(dims[row])==size(A,1) && prod(dims[col])==size(A,2),"Dimensions mismatch")
	X = reshape(A,[dims[row];dims[col]]...)
	permutedims(X,invperm([row;col]))
end

"""
    mkrontv(X,Y,v,n,t='n')

Matricized Kronecker product of tensors X and Y times vector v (n-mode multiplication): (X ⊗ Y)ₙv.
If t='t', transpose matricized Kronecker product: (X ⊗ Y)ᵀₙv.
If v is a matrix, multiply column by column.
"""
function mkrontv(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},v::AbstractVector{<:Number},n::Integer,t='n') where N
  I1=size(X1)
  I2=size(X2)
  kronsize=tuple(([I1...].*[I2...])...);
  ind=setdiff(1:N,n) #all indices but n
  X1n=tenmat(X1,n);
  X2n=tenmat(X2,n);
  perfect_shuffle=Int[ [2*k-1 for k=1:N-1]; [2*k for k=1:N-1] ]
  if t=='n'
    @assert(length(v) == prod(kronsize[ind]),"Vector is of inapropriate size.")
    tenshape=vec([[I2[ind]...] [I1[ind]...]]');
    vperm=permutedims(reshape(v,tenshape...),perfect_shuffle);
    vec(X2n*reshape(vperm,size(X2n,2),size(X1n,2))*X1n')
  elseif t=='t'
    @assert(length(v) == kronsize[n],"Vector is of inapropriate size.")
    if I1[n]*prod(I2)+prod(I1)*prod(I2[ind]) > I2[n]*prod(I1)+prod(I2)*prod(I1[ind])
      W=(X2n'*reshape(v,size(X2n,1),size(X1n,1)))*X1n
    else
      W=X2n'*(reshape(v,size(X2n,1),size(X1n,1))*X1n)
    end
    tenshape=[[I2[ind]...];[I1[ind]...]];
    vec(permutedims(reshape(W,tenshape...),invperm(perfect_shuffle)))
    end
end

#Matricized Kronecker product times matrix - column by column.
function mkrontv(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::AbstractMatrix{<:Number},n::Integer,t='n') where N
  if sort(collect(size(vec(M))))[1]==1
    return mkrontv(A,B,vec(M));
  end
  I1=size(X1);
  I2=size(X2);
  kronsize=([I1...].*[I2...]);
  ind=setdiff(1:N,n) #all indices but n
  if t=='n'
    Mprod=zeros(kronsize[n],size(M,2))
  else
    Mprod=zeros(prod(kronsize[ind]),size(M,2))
  end
  for i=1:size(M,2)
    Mprod[:,i]=mkrontv(X1,X2,M[:,i],n,t);
  end
  Mprod
end
function mkrontm(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N},M::AbstractMatrix{<:Number},n::Integer,t='n') where N
  @warn "Function mkrontm is depricated. Use mkrontv."
  mkrontv(X1,X2,M,n,t)
end

"""
    mrank(X[,tol])
    mrank(X::ttensor[,tol])

Multilinear rank of a tensor with optionally given tolerance.
"""
function mrank(X::AbstractArray{<:Number,N}) where N
   ntuple(n->nrank(X,n),N)
end
function mrank(X::AbstractArray{<:Number,N},tol::Number) where N
   ntuple(n->nrank(X,n,tol),N)
end

"""
    mttkrp(X,M,n)
    mttkrp(X::ttensor,M,n)
    mttkrp(X::ktensor,M,n)

Mode-n matricized tensor X times Khatri-Rao product of matrices from M (except nth) in reverse order.
"""
function mttkrp(X::AbstractArray{<:Number,N},M::MatrixCell,n::Integer) where N
  @assert(N-1<=length(M)<=N,"Wrong number of matrices")
  if length(M)==N-1  #if nth matrix not defined
    push!(M,M[end])
    [M[m]=M[m-1] for m=N-1:-1:n+1]
  end
  modes=setdiff(1:N,n)
  sz=[size(X)...]
  K=size(M[modes[1]],2)
  @assert(!any(map(Bool,[size(M[m],2)-K for m in modes])),"Matrices must have the same number of columns.")
  @assert(!any(map(Bool,[size(M[m],1)-sz[m] for m in modes])),"Matrices are of wrong size.")
  Xn=tenmat(X,n)
  Xn*khatrirao(reverse(M[modes]))
end
mttkrp(X::AbstractArray{<:Number},M::Array{Matrix{T},1},n::Integer) where {T<:Number}=mttkrp(X,MatrixCell(M),n)

"""
    neye(dims)

Identity tensor of a given dimension. Generalization of eye.
"""
function neye(dims::Vector{<:Integer})
  dims=tuple(dims...)
  A=zeros(dims)
  R=CartesianIndices(dims)
  Ifirst=first(R)
  Iend=last(R)
  while Ifirst<=Iend
      A[Ifirst]=1
      Ifirst=Ifirst+1*one(Ifirst)
  end
  A
end
function neye(d1,d2...)
  if isempty(d2)
    return neye([d1...])
  end
  dims=[d1]
  for d in d2
    push!(dims,d)
  end
  neye(dims)
end
function neye(d::Integer;order=0)
  @assert(order>0,"Wrong input.")
  neye(repeat([d],order,1)[:])
end

"""
    nrank(X,n[,tol])
    nrank(X::ttensor,n[,tol])

Rank of the n-mode matricization of a tensor X (n-rank).
"""
function nrank(X::AbstractArray{<:Number},n::Integer)
  rank(tenmat(X,n))
end
function nrank(X::AbstractArray{<:Number},n::Integer,tol::Number)
  rank(tenmat(X,n),tol)
end

"""
    nvecs(X,n,r=0;flipsign=false,svds=false)
    nvecs(X::ttensor,n,r=0;flipsign=false)
    nvecs(X::ktensor,n,r=0;flipsign=false)

Computes the r leading singular vectors of mode-n matricization of a tensor X.
Works with XₙXₙᵀ.

## Arguments:
- `flipsign=true`: Make the largest magnitude element be positive.
- `svds=true`: Use svds on Xₙ rather than eigs on XₙXₙᵀ.
"""
function nvecs(X::AbstractArray{<:Number},n::Integer,r=0;flipsign=false,svds=false)
  if r==0
    r=size(X,n)
  end
  Xn=tenmat(X,n)
  if svds
    #U=svds(Xn,nsv=r)[1][:U]
    #if size(U,2)<r
      U=svdfact(Xn)[:U][1:r]
    #end
  else
    G=Symmetric(Xn*Xn') #Gramian
    #U=eigs(G,nev=r,which=:LM)[2] #has bugs!
    #if size(U,2)<r
       U=eigen(G).vectors[:,end:-1:end-r+1]
    #end
  end
  if flipsign
      maxind = findmax(abs.(U),1)[2]
      for i = 1:r
          ind=ind2sub(size(U),maxind[i])
          if U[ind...] < 0
             U[:,ind[2]] = U[:,ind[2]] * -1
          end
      end
  end
  U
end

#Squeeze all singleton dimensions. **Documentation in Base.jl.
function dropdims(A::AbstractArray{<:Number})
  sz=size(A)
  sdims=findall(sz.==1) #singleton dimensions
  dropdims(A,dims=tuple(sdims...))
end

"""
    sthosvd(X;reqrank::Vector,order::Vector,tol::Number)

Sequentially truncated HOSVD of a tensor X of predifined rank and processing order p.
"""
function sthosvd(X::AbstractArray{<:Number,N};reqrank=[],order=[],tol=1e-8) where N
    if order==[]
        order=collect(1:N)
    end
    @assert(N==length(order),"Dimensions mismatch")
    if reqrank!=[]
        @assert(N==length(reqrank),"Dimensions mismatch")
    end
	sz=[size(X)...]
	fmat=MatrixCell(undef,N)
	for n in order
		Xn=tenmat(X,n)
		U,S,V=svd(Xn)
        if reqrank!=[]
            fmat[n]=U[:,1:reqrank[n]]
		    Xn=Diagonal(S[1:reqrank[n]])*V'[1:reqrank[n],:]
            sz[n]=reqrank[n]
        else
            K=findall(x-> x>tol ? true : false,S)
            fmat[n]=U[:,K]
		    Xn=Diagonal(S[K])*V[:,K]'
            sz[n]=length(K)
        end
		X=matten(Xn,n,sz)
	end
	ttensor(X,fmat)
end

"""
    tenmat(X,n)
    tenmat(X,row=[],col=[])
    tenmat(X::ttensor,n)
    tenmat(X::ktensor,n)

Mode-n matricization of a tensor or matricization by row and column vectors R and C.
"""
function tenmat(X::AbstractArray{<:Number,N},n::Integer) where N
	@assert(n<=ndims(X),"Mode exceedes number of dimensions")
	sz=size(X)
	m=setdiff(1:N,n)
    if [n;m]!=collect(1:N)
        X=permutedims(X,[n;m])
    end
	reshape(X,sz[n],prod(sz[m]))
end

function tenmat(X::AbstractArray{<:Number,N};row=[],col=[]) where N
    @assert(row!=[] || col!=[],"Al least one of row and col needs to be specified.")
    if row!=[] && col!=[]
        @assert(sort([row;col])==collect(1:N),"Incorrect mode partitioning.")
    elseif row==[]
        @assert(!(false in [c in collect(1:N) for c in col]),"Incorrect modes.")
        if isa(col,Integer)
            col=[col]
        end
        row=collect(1:N)
        deleteat!(row,sort(col))
    else
        @assert(!(false in [r in collect(1:N) for r in row]),"Incorrect modes.")
        if isa(row,Integer)
            row=[row]
        end
        col=collect(1:N)
        deleteat!(col,sort(row))
    end
	sz=size(X)
    J=prod(sz[row]);K=prod(sz[col])
    if [row;col]!=collect(1:N)
        X=permutedims(X,[row;col])
    end
	reshape(X,J,K)
end

"""
    tkron(X,Y)

Kronecker product of two tensors X and Y. Direct generalization of Kronecker product of matrices.
"""
function tkron(X1::AbstractArray{<:Number,N},X2::AbstractArray{<:Number,N}) where N
  if N<3
    return kron(X1,X2)
  end
  s1=size(X1)
  s2=size(X2)
  Xk=zeros(s1.*s2)
  R1=CartesianIndices(s1)
  R2=CartesianIndices(s2)
  Il=last(R2)
  i=0
  for I1 in R1
    for I2 in R2
      Xk[I2+ewprod((I1-1*one(I1)),Il)]=X1[I1]*X2[I2]
    end
    i+=1
  end
  Xk
end

"""
    ttm(X,M[,modes,t='n'])
    ttm(X::ttensor,M[,modes,t='n'])

Tensor times matrix (n-mode product):  X x₁ M₁ x₂ M₂ x₃ ⋯ xₙ Mₙ
Default modes: 1:length(M).
If t='t', transpose matrices from M.
"""
function ttm(X::AbstractArray{<:Number,N},A::MatrixCell,modes::AbstractVector{<:Integer},t='n') where N
  M=deepcopy(A)
  if t=='t'
    [M[n]=transpose(M[n]) for n=1:length(M)]
  end
  @assert(length(modes)<=length(M),"Too few matrices.")
  @assert(length(M)<=N,"Too many matrices.")
  sz=[size(X)...]
  if length(modes) < length(M)
    M=M[modes] #discard matrices not needed for multiplication
  end
  R=copy(sz)
  for n=1:length(modes)
    R[modes[n]]=size(M[n],1) #vector of rₖ
  end
  #Order of multiplication - if X is i₁ × i₂ × ... × iₙ and Mₖ is rₖ × iₖ, sort by largest possible dimension reduction iₖ-rₖ
  p=sortperm(sz[modes]-R[modes],rev=true)
  M=M[p]
  modes=modes[p]
  for n=1:length(M)
	 @assert(sz[modes[n]] == size(M[n],2),"Dimensions mismatch")
     Xn=tenmat(X,modes[n])
	 sz[modes[n]]=size(M[n],1)
	 X=matten(M[n]*Xn,modes[n],sz)
  end
  X
end
ttm(X::AbstractArray{<:Number},M::AbstractMatrix{<:Number},n::Integer,t='n')=ttm(X,MatrixCell([M]),[n],t)
ttm(X::AbstractArray{<:Number},M::MatrixCell,t::Char)=ttm(X,M,1:length(M),t)
ttm(X::AbstractArray{<:Number},M::MatrixCell)=ttm(X,M,1:length(M))
ttm(X::AbstractArray{<:Number},M::MatrixCell,R::AbstractRange{<:Integer},t::Char)=ttm(X,M,collect(R),t)
ttm(X::AbstractArray{<:Number},M::MatrixCell,R::AbstractRange{<:Integer})=ttm(X,M,collect(R))
function ttm(X::AbstractArray{<:Number,N},M::MatrixCell,n::Integer,t='n') where N
	if n>0
		ttm(X,M[n],n,t)
	else
		modes=setdiff(1:N,-n)
		ttm(X,M,modes,t)
	end
end
#If array of matrices isn't defined as MatrixCell, but as M=[M1,M2,...,Mn]:
ttm(X::AbstractArray{<:Number},M::Array{Matrix{T},1},modes::Array{<:Integer,1},t='n') where{T<:Number}=ttm(X,MatrixCell(M),modes,t)
ttm(X::AbstractArray{<:Number},M::Array{Matrix{T},1},t::Char) where{T<:Number}=ttm(X,MatrixCell(M),t)
ttm(X::AbstractArray{<:Number},M::Array{Matrix{T},1}) where{T<:Number}=ttm(X,MatrixCell(M))
ttm(X::AbstractArray{<:Number},M::Array{Matrix{T},1},R::AbstractRange{<:Integer},t::Char) where{T<:Number}=ttm(X,MatrixCell(M),R,t)
ttm(X::AbstractArray{<:Number},M::Array{Matrix{T},1},R::AbstractRange{<:Integer}) where{T<:Number}=ttm(X,MatrixCell(M),R)
ttm(X::AbstractArray{<:Number},M::Array{Matrix{T},1},n::Integer,t='n') where{T<:Number}=ttm(X,MatrixCell(M),n,t)

"""
    ttt(X,Y)

Outer product of two tensors.
"""
function ttt(X1::AbstractArray{<:Number},X2::AbstractArray{<:Number})
  sz=tuple([[size(X1)...];[size(X2)...]]...);
  Xprod=vec(X1)*vec(X2)';
  reshape(Xprod,sz)
end

"""
    ttv(X,V[,modes])

Tensor times vectors (n-mode product):  X x₁ V₁ x₂ V₂ x₃ ⋯ xₙ Vₙ.
Default modes: 1:length(M).
"""
function ttv(X::AbstractArray{<:Number,N},V::VectorCell,modes::AbstractVector{<:Integer}) where N
  remmodes=setdiff(1:N,modes)'
  if N > 1
    X=permutedims(X,[remmodes modes'])
  end
  sz=size(X)
  if length(modes) < length(V)
    V=V[modes]
  end
  M=N
  for n=length(modes):-1:1
    X=reshape(X,prod(sz[1:M-1]),sz[M])
    X=X*V[n]
    M-=1
  end
  if M>0
    X=reshape(X,sz[1:M])
  end
  X
end
#ttv(X::AbstractArray{<:Number,N},v::AbstractVector{<:Number},n::Integer) where N =ttv(X,[v],[n])
ttv(X::AbstractArray{<:Number,N},v::AbstractVector{<:Number},n::Integer) where N =ttv(X,VectorCell([v]),[n])
ttv(X::AbstractArray{<:Number,N},V::VectorCell) where N = ttv(X,V,collect(1:length(V)))
function ttv(X::AbstractArray{<:Number,N},V::VectorCell,n::Integer) where N
	if n>0
		ttv(X,V[n],n)
	else
		modes=setdiff(1:N,-n)
		ttv(X,V,modes)
	end
end
#If array of vectors isn't defined as VectorCell, but as V=[v1,v2,...,vn]:
ttv(X::AbstractArray{<:Number,N},V::Array{Vector{T},1},modes::AbstractVector{<:Integer}) where {T<:Number,N} = ttv(X,VectorCell(V),modes)
ttv(X::AbstractArray{<:Number,N},V::Array{Vector{T},1}) where {T<:Number,N} = ttv(X,VectorCell(V))
ttv(X::AbstractArray{<:Number,N},V::Array{Vector{T},1},n::Integer) where {T<:Number,N} = ttv(X,VectorCell(V),n)


text:
# Julia wrapper for header: /opt/ffmpeg/include/libswresample/swresample.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    swr_get_class,
    swr_alloc,
    swr_init,
    swr_is_initialized,
    swr_alloc_set_opts,
    swr_free,
    swr_close,
    swr_convert,
    swr_next_pts,
    swr_set_compensation,
    swr_set_channel_mapping,
    swr_set_matrix,
    swr_drop_output,
    swr_inject_silence,
    swr_get_delay,
    #swresample_version,
    swresample_configuration,
    swresample_license

function swr_get_class()
    ccall((:swr_get_class,libswresample),Ptr{AVClass},())
end

function swr_alloc()
    ccall((:swr_alloc,libswresample),Ptr{SwrContext},())
end

function swr_init(s)
    ccall((:swr_init,libswresample),Cint,(Ptr{SwrContext},),s)
end

function swr_is_initialized(s)
    ccall((:swr_is_initialized,libswresample),Cint,(Ptr{SwrContext},),s)
end

function swr_alloc_set_opts(s,out_ch_layout::Int64,out_sample_fmt::AVSampleFormat,out_sample_rate::Integer,in_ch_layout::Int64,in_sample_fmt::AVSampleFormat,in_sample_rate::Integer,log_offset::Integer,log_ctx)
    ccall((:swr_alloc_set_opts,libswresample),Ptr{SwrContext},(Ptr{SwrContext},Int64,AVSampleFormat,Cint,Int64,AVSampleFormat,Cint,Cint,Ptr{Void}),s,out_ch_layout,out_sample_fmt,out_sample_rate,in_ch_layout,in_sample_fmt,in_sample_rate,log_offset,log_ctx)
end

function swr_free(s)
    ccall((:swr_free,libswresample),Void,(Ptr{Ptr{SwrContext}},),s)
end

function swr_close(s)
    ccall((:swr_close,libswresample),Void,(Ptr{SwrContext},),s)
end

function swr_convert(s,out,out_count::Integer,_in,in_count::Integer)
    ccall((:swr_convert,libswresample),Cint,(Ptr{SwrContext},Ptr{Ptr{UInt8}},Cint,Ptr{Ptr{UInt8}},Cint),s,out,out_count,_in,in_count)
end

function swr_next_pts(s,pts::Int64)
    ccall((:swr_next_pts,libswresample),Int64,(Ptr{SwrContext},Int64),s,pts)
end

function swr_set_compensation(s,sample_delta::Integer,compensation_distance::Integer)
    ccall((:swr_set_compensation,libswresample),Cint,(Ptr{SwrContext},Cint,Cint),s,sample_delta,compensation_distance)
end

function swr_set_channel_mapping(s,channel_map)
    ccall((:swr_set_channel_mapping,libswresample),Cint,(Ptr{SwrContext},Ptr{Cint}),s,channel_map)
end

function swr_set_matrix(s,matrix,stride::Integer)
    ccall((:swr_set_matrix,libswresample),Cint,(Ptr{SwrContext},Ptr{Cdouble},Cint),s,matrix,stride)
end

function swr_drop_output(s,count::Integer)
    ccall((:swr_drop_output,libswresample),Cint,(Ptr{SwrContext},Cint),s,count)
end

function swr_inject_silence(s,count::Integer)
    ccall((:swr_inject_silence,libswresample),Cint,(Ptr{SwrContext},Cint),s,count)
end

function swr_get_delay(s,base::Int64)
    ccall((:swr_get_delay,libswresample),Int64,(Ptr{SwrContext},Int64),s,base)
end

function swresample_version()
    ccall((:swresample_version,libswresample),UInt32,())
end

function swresample_configuration()
    ccall((:swresample_configuration,libswresample),Ptr{UInt8},())
end

function swresample_license()
    ccall((:swresample_license,libswresample),Ptr{UInt8},())
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
using Test
using StaticArrays
using HMatrices
using LinearAlgebra

include(joinpath(HMatrices.PROJECT_ROOT,"test","testutils.jl"))

@testset "Assemble" begin
    m,n  = 1_000,1_000
    X  = Y   = rand(SVector{3,Float64},m)
    splitter  = CardinalitySplitter(nmax=20)
    Xclt = Yclt = ClusterTree(X,splitter)
    adm       = StrongAdmissibilityStd(eta=3)
    rtol      = 1e-5
    comp      = PartialACA(rtol=rtol)
    # Laplace
    for threads in (true,false)
        K         = laplace_matrix(X,Y)
        H         = assemble_hmat(K,Xclt,Yclt;adm,comp,threads)
        @test norm(Matrix(K) - Matrix(H;global_index=true)) < rtol*norm(Matrix(K))
        H         = assemble_hmat(K;threads,distributed=false)
        @test norm(Matrix(K) - Matrix(H;global_index=true)) < rtol*norm(Matrix(K))
        adjH      = adjoint(H)
        H_full    = Matrix(H;global_index=false)
        adjH_full = adjoint(H_full)
        @testset "getindex" begin
            HMatrices.enable_getindex()
            @test H_full[8,64] ≈ H[8,64]
            HMatrices.disable_getindex()
            @test_throws ErrorException H_full[8,64] ≈ H[8,64]
            HMatrices.enable_getindex()
            @test H_full[:,666] ≈ H[:,666]
            @test adjH_full[:,666] ≈ adjH[:,666]
        end
        # Elastostatic
        K         = elastostatic_matrix(X,Y,1.1,1.2)
        H         = assemble_hmat(K,Xclt,Yclt;adm,comp,threads)
        @test norm(Matrix(K) - Matrix(H;global_index=true)) < rtol*norm(Matrix(K))
    end
end

max_stars_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
@testset "Example 1" begin
    Mod.@polyvar x y
    p = 2x + 3.0x*y^2 + y
    @test differentiate(p, x) == 2 + 3y^2
    @test differentiate.(p, (x, y)) == (2 + 3y^2, 6x*y + 1)
    @test p((x, y)=>(y, x)) == 2y + 3y*x^2 + x
    @test subs(p, y=>x^2) == 2x + 3x^5 + x^2
    @test p(x=>1, y=>2) == 2 + 3*4 + 2
end


text:
module Builder

using ..Configs: Configs
using ..ZSHCompletions
using ..Comonicon: get_version, Tools
using PackageCompiler
using Logging
using TOML: TOML
using Libdl: Libdl
using Pkg: Pkg
using UUIDs: UUID
using Scratch: get_scratch!

const COMONICON_URL = "https://github.com/comonicon/Comonicon.jl"

include("cli.jl")
include("install.jl")
include("sysimg.jl")
include("app.jl")
include("tarball.jl")
include("rcfile.jl")

end


text:
#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the methods in `jltest/utils.jl`.

Notes
-----
* For the unit tests in this files, failures and errors are expected.
"""

# --- Imports

# Standard library
using Logging
using Test
using Test: DefaultTestSet

# External packages
using Coverage: Coverage
using Suppressor

# Local modules
using TestTools.jltest

# --- Tests

@testset EnhancedTestSet "jltest.run_tests(): basic tests" begin
    # --- Preparations

    # Construct path to test directory
    test_dir = joinpath(@__DIR__, "data-basic-tests")
    test_dir_relpath = relpath(test_dir)

    # Precompute commonly used values
    some_tests_file = joinpath(test_dir, "some_tests.jl")
    expected_output_some_tests = "$(joinpath(test_dir_relpath, "some_tests")): .."

    some_tests_no_testset_file = joinpath(test_dir, "some_tests_no_testset.jl")
    expected_output_some_tests_no_testset = "$(joinpath(test_dir_relpath, "some_tests_no_testset")): .."

    failing_tests_file = joinpath(test_dir, "failing_tests.jl")
    expected_output_failing_tests = strip(
        """
        $(joinpath(test_dir_relpath, "failing_tests")): .
        =====================================================
        failing tests: Test Failed at $(failing_tests_file):27
          Expression: 2 == 1
           Evaluated: 2 == 1

        Stacktrace:
         [1]
        """
    )

    failing_tests_no_testset_file = joinpath(test_dir, "failing_tests_no_testset.jl")
    expected_output_failing_tests_no_testset = strip(
        """
        $(joinpath(test_dir_relpath, "failing_tests_no_testset")): .
        =====================================================
        test set: Test Failed at $(failing_tests_no_testset_file):26
          Expression: 2 == 1
           Evaluated: 2 == 1

        Stacktrace:
          [1]
        """
    )

    more_tests_file = joinpath(test_dir, "subdir", "more_tests.jl")
    expected_output_more_tests = "$(joinpath(test_dir_relpath, "subdir", "more_tests")): .."

    # --- Tests for run_tests(tests::Vector{<:AbstractString})

    # `tests` contains tests named with ".jl" extension
    tests = [some_tests_file]
    output = strip(@capture_out begin
        run_tests(tests)
    end)
    @test output == expected_output_some_tests

    # `tests` contains tests named without ".jl" extension
    some_tests_no_testset_file_without_extension = some_tests_no_testset_file[1:(end - 3)]
    tests = [some_tests_no_testset_file_without_extension]
    output = strip(@capture_out begin
        run_tests(tests)
    end)
    @test output == expected_output_some_tests_no_testset

    # `tests` contains only a directory
    tests = [test_dir]
    output = strip(@capture_out begin
        run_tests(tests)
    end)

    expected_output_lines = [
        expected_output_some_tests,
        expected_output_some_tests_no_testset,
        expected_output_failing_tests,
        expected_output_failing_tests_no_testset,
        expected_output_more_tests,
    ]
    for line in expected_output_lines
        @test occursin(line, output)
    end

    # `tests` contains both directories and files
    tests = [test_dir, some_tests_file]
    output = strip(@capture_out begin
        run_tests(tests)
    end)

    expected_output_lines = [
        expected_output_some_tests,
        expected_output_some_tests_no_testset,
        expected_output_failing_tests,
        expected_output_failing_tests_no_testset,
        expected_output_more_tests,
    ]
    for line in expected_output_lines
        @test occursin(line, output)
    end

    output_lines = split(output, '\n')
    @test count(i -> (i == expected_output_some_tests), output_lines) == 2

    # `tests` is empty list
    tests = Vector{String}()
    @test_throws ArgumentError run_tests(tests)

    # --- Tests for run_tests(tests::AbstractString)

    # `tests` is a string
    tests = some_tests_file
    output = strip(@capture_out begin
        run_tests(tests)
    end)

    @test output == expected_output_some_tests

    # `tests` is empty string
    tests = ""
    @test_throws ArgumentError run_tests(tests)

    # --- Keyword arguments tests

    # desc
    tests = [failing_tests_no_testset_file]
    desc = "test description"
    output = strip(@capture_out begin
        run_tests(tests; desc=desc)
    end)
    output_line_three = split(output, '\n')[3]
    @test startswith(output_line_three, desc)

    # test_set_type
    tests = [failing_tests_file]
    output = strip(@capture_out begin
        run_tests(tests; test_set_type=DefaultTestSet)
    end)
    expected_prefix =
        "$(joinpath(test_dir_relpath, "failing_tests")): failing tests" *
        ": Test Failed at $(joinpath(test_dir, "failing_tests.jl")):27"
    @test startswith(output, expected_prefix)

    # test_set_type = nothing
    tests = [failing_tests_file]
    output = strip(@capture_out begin
        run_tests(tests; test_set_type=nothing)
    end)
    @test startswith(output, expected_output_failing_tests)

    # recursive = false
    tests = [test_dir]
    output = strip(@capture_out begin
        run_tests(tests; recursive=false)
    end)

    expected_output_lines = [
        expected_output_some_tests,
        expected_output_some_tests_no_testset,
        expected_output_failing_tests,
        expected_output_failing_tests_no_testset,
    ]
    for line in expected_output_lines
        @test occursin(line, output)
    end
    @test !occursin(expected_output_more_tests, output)

    # exclude_runtests = false
    tests = [test_dir]
    output = strip(@capture_out begin
        run_tests(tests; exclude_runtests=false)
    end)

    expected_output_lines = [
        expected_output_some_tests,
        expected_output_some_tests_no_testset,
        expected_output_failing_tests,
        expected_output_failing_tests_no_testset,
        "$(joinpath(test_dir_relpath, "failing_tests")):",
    ]
    for line in expected_output_lines
        @test occursin(line, output)
    end
end

@testset EnhancedTestSet "jltest.run_tests(): log message tests" begin
    # --- Preparations

    # Construct path to test directory
    test_dir = joinpath(@__DIR__, "data-log-message-tests")
    test_dir_relpath = relpath(test_dir)

    # Set up Julia environment
    push!(LOAD_PATH, test_dir)

    # Temporarily enable logging at all levels
    env_julia_debug_save = get(ENV, "JULIA_DEBUG", nothing)
    ENV["JULIA_DEBUG"] = "all"

    # Precompute commonly used values
    log_message_tests_file = joinpath(test_dir, "log_message_tests.jl")
    expected_output_log_message_tests = "$(joinpath(test_dir_relpath, "log_message_tests")):"

    location_prefix =
        "TestTools.jltest.##$(joinpath(test_dir_relpath, "log_message_tests"))#[0-9]+ " *
        "$(Base.contractuser(log_message_tests_file))"
    if Sys.iswindows()
        location_prefix = replace(location_prefix, "\\" => "\\\\")
    end

    expected_log_messages_log_message_tests = [
        "[ Warning: Single line @warn message test",
        Regex(strip("""
                    ┌ Warning: Multi-line @warn message test.
                    │ Second line.
                    │ Third line.
                    └ @ $(location_prefix):[0-9]+
                    """)),
        "[ Info: Single line @info message test",
        strip("""
              ┌ Info: Multi-line @info message test.
              │ Second line.
              └ Third line.
              """),
        "[ Debug: Single line @debug message test",
        Regex(strip("""
                    ┌ Debug: Multi-line @debug message test.
                    │ Second line.
                    │ Third line.
                    └ @ $(location_prefix):[0-9]+
                    """)),
    ]

    missing_dependencies_tests_file = joinpath(test_dir, "missing_dependencies_tests.jl")
    expected_output_missing_dependencies_tests = "$(joinpath(test_dir_relpath, "missing_dependencies_tests")):"
    expected_log_messages_missing_dependencies_tests = "[ Info: Non-missing dependency log message"

    # --- Tests

    local output

    # Case: "Package TestTools does not have ... in its dependencies" warning suppressed
    tests = [missing_dependencies_tests_file]
    log_msg = strip(@capture_err begin
        output = strip(@capture_out begin
            run_tests(tests)
        end)
    end)
    @test output == expected_output_missing_dependencies_tests
    @test occursin(expected_log_messages_missing_dependencies_tests, log_msg)

    # Case: only non-missing dependency log messages
    tests = [log_message_tests_file]
    log_msg = strip(@capture_err begin
        output = strip(@capture_out begin
            run_tests(tests)
        end)
    end)

    @test output == expected_output_log_message_tests

    for message in expected_log_messages_log_message_tests
        @test occursin(message, log_msg)
    end

    # --- Clean up

    # Remove Manifest.toml
    rm(joinpath(test_dir, "Manifest.toml"); force=true)

    # Restore LOAD_PATH
    filter!(x -> x != test_dir, LOAD_PATH)

    # Restore state of ENV
    if isnothing(env_julia_debug_save)
        delete!(ENV, "JULIA_DEBUG")
    else
        ENV["JULIA_DEBUG"] = env_julia_debug_save
    end
end

@testset EnhancedTestSet "jltest.run_tests(): current directory checks" begin
    # --- Preparations

    # Get current directory
    cwd = pwd()

    # Construct path to test directory
    test_dir = joinpath(@__DIR__, "data-directory-change-tests")

    # Change to test directory
    cd(test_dir)

    # Precompute commonly used values
    change_dir_file = joinpath(test_dir, "change_dir.jl")
    expected_output_change_dir = "$(relpath(joinpath(test_dir, "change_dir"), test_dir)): "

    # Delete old coverage data files
    @suppress begin
        Coverage.clean_folder(test_dir)
    end

    # --- Tests

    # Case: change directory before running test file
    check_dir_file = joinpath(test_dir, "check_dir.jl")
    tests = [change_dir_file, check_dir_file]
    output = strip(@capture_out begin
        run_tests(tests)
    end)

    expected_output = join(
        [
            expected_output_change_dir,
            "$(relpath(joinpath(test_dir, "check_dir"), test_dir)): .",
        ],
        '\n',
    )
    @test output == expected_output

    # Case: change directory before running tests in a directory
    tests = [change_dir_file, joinpath(test_dir, "subdir")]
    output = strip(@capture_out begin
        run_tests(tests)
    end)

    expected_output = join(
        [
            expected_output_change_dir,
            "$(relpath(joinpath(test_dir, "subdir", "check_dir"), test_dir)): .",
        ],
        '\n',
    )
    @test output == expected_output

    # --- Clean up

    # Restore current directory
    cd(cwd)
end

@testset EnhancedTestSet "jltest.run_tests(): JLTEST_FAIL_FAST tests" begin
    # --- Preparations

    # Construct path to test directory
    test_dir = joinpath(@__DIR__, "data-basic-tests")
    test_dir_relpath = relpath(test_dir)

    # Save original value of JLTEST_FAIL_FAST environment variable
    env_jltest_fail_fast_original = get(ENV, "JLTEST_FAIL_FAST", nothing)

    # Precompute commonly used values
    some_tests_file = joinpath(test_dir, "some_tests.jl")
    expected_output_some_tests = "$(joinpath(test_dir_relpath, "some_tests")): .."

    failing_tests_file = joinpath(test_dir, "failing_tests.jl")
    expected_output_failing_tests = strip(
        """
        $(joinpath(test_dir_relpath, "failing_tests")): .
        =====================================================
        failing tests: Test Failed at $(failing_tests_file):27
          Expression: 2 == 1
           Evaluated: 2 == 1

        Stacktrace:
         [1]
        """
    )

    expected_output_failing_tests_fail_fast = strip(
        """
        $(joinpath(test_dir_relpath, "failing_tests")): .
        =====================================================
        Test Failed at $(failing_tests_file):27
          Expression: 2 == 1
           Evaluated: 2 == 1

        =====================================================
        Error During Test at
        """
    )

    more_tests_file = joinpath(test_dir, "subdir", "more_tests.jl")
    expected_output_more_tests = "$(joinpath(test_dir_relpath, "subdir", "more_tests")): .."

    # --- Tests

    # Case: ENV["JLTEST_FAIL_FAST"] = false
    ENV["JLTEST_FAIL_FAST"] = false
    tests = [failing_tests_file, some_tests_file]
    local error = nothing
    output = strip(@capture_out begin
        try
            jltest.run_tests(tests)
        catch error
        end
    end)

    if isnothing(env_jltest_fail_fast_original)  # Restore ENV["JLTEST_FAIL_FAST"]
        delete!(ENV, "JLTEST_FAIL_FAST")
    else
        ENV["JLTEST_FAIL_FAST"] = env_jltest_fail_fast_original
    end

    @test startswith(output, expected_output_failing_tests)
    @test occursin(expected_output_some_tests, output)
    @test isnothing(error)

    # Case: ENV["JLTEST_FAIL_FAST"] is undefined
    # Expect: same behavior as ENV["JLTEST_FAIL_FAST"] = false
    delete!(ENV, "JLTEST_FAIL_FAST")
    tests = [failing_tests_file, some_tests_file]
    local error = nothing
    output = strip(@capture_out begin
        try
            jltest.run_tests(tests)
        catch error
        end
    end)

    if isnothing(env_jltest_fail_fast_original)  # Restore ENV["JLTEST_FAIL_FAST"]
        delete!(ENV, "JLTEST_FAIL_FAST")
    else
        ENV["JLTEST_FAIL_FAST"] = env_jltest_fail_fast_original
    end

    @test startswith(output, expected_output_failing_tests)
    @test occursin(expected_output_some_tests, output)
    @test isnothing(error)

    # Case: ENV["JLTEST_FAIL_FAST"] = true
    ENV["JLTEST_FAIL_FAST"] = "true"
    tests = [failing_tests_file, some_tests_file]
    local error = nothing
    output = strip(@capture_out begin
        try
            jltest.run_tests(tests)
        catch error
        end
    end)

    if isnothing(env_jltest_fail_fast_original)  # Restore ENV["JLTEST_FAIL_FAST"]
        delete!(ENV, "JLTEST_FAIL_FAST")
    else
        ENV["JLTEST_FAIL_FAST"] = env_jltest_fail_fast_original
    end

    @test startswith(output, expected_output_failing_tests_fail_fast)
    @test !occursin(expected_output_some_tests, output)
    @test error isa Test.FallbackTestSetException
    @test error.msg == "There was an error during testing"
end

@testset EnhancedTestSet "jltest.find_tests()" begin

    # --- directory without subdirectories

    test_dir = joinpath(@__DIR__, "data-find-tests", "subdir")
    tests = Set(find_tests(test_dir))
    expected_tests = Set([joinpath(test_dir, "more_tests.jl")])
    @test tests == expected_tests

    # --- directory with subdirectories

    test_dir = joinpath(@__DIR__, "data-find-tests")
    tests = Set(find_tests(test_dir))
    expected_tests = Set([
        joinpath(test_dir, file) for
        file in ["some_tests.jl", joinpath("subdir", "more_tests.jl")]
    ])
    @test tests == expected_tests

    # --- Keyword arguments tests

    # recursive = false
    test_dir = joinpath(@__DIR__, "data-find-tests")
    tests = Set(find_tests(test_dir; recursive=false))
    expected_tests = Set([joinpath(test_dir, "some_tests.jl")])
    @test tests == expected_tests

    # exclude_runtests = false
    test_dir = joinpath(@__DIR__, "data-find-tests")
    tests = Set(find_tests(test_dir; exclude_runtests=false))
    expected_tests = Set([
        joinpath(test_dir, file) for file in [
            "runtests.jl",
            "some_tests.jl",
            joinpath("subdir", "more_tests.jl"),
            joinpath("subdir", "runtests.jl"),
        ]
    ])
    @test tests == expected_tests
end

@testset EnhancedTestSet "jltest: Pkg.test() tests" begin
    # --- Preparations

    # Get current directory
    cwd = pwd()

    test_pkg_dir = joinpath(@__DIR__, "data-missing-package-dependency", "TestPackage")
    cmd_options = `--startup-file=no --project=@. -O0`
    cmd = `julia $(cmd_options) -e 'import Pkg; Pkg.test(coverage=true)'`

    # --- Tests

    test_error = strip(@capture_out begin
        try
            @suppress_err begin
                Base.run(Cmd(cmd; dir=test_pkg_dir); wait=true)
            end
        catch process_error
            @test process_error isa ProcessFailedException
        end
    end)

    src_error_file = abspath(
        joinpath(dirname(dirname(@__DIR__)), "src", "jltest", "utils.jl")
    )
    test_error_file = joinpath(test_pkg_dir, "test", "missing_dependency_tests.jl")
    expected_test_error = strip(
        """
        missing_dependency_tests:
        =====================================================
        test set: Error During Test at $(src_error_file):285
          Got exception outside of a @test
          The test environment is missing InteractiveUtils from its dependencies.
          Error occurred at $(test_error_file):22
          Stacktrace:
            [1]
        """
    )
    @test startswith(test_error, expected_test_error)

    # --- Clean up

    # Remove Manifest.toml
    rm(joinpath(test_pkg_dir, "Manifest.toml"); force=true)
end

# --- Emit message about expected failures and errors

println()
@info "For $(basename(@__FILE__)), 13 failures and 0 errors are expected."


text:
using TcpInstruments
using Test

m = initialize(KeysightDMM34465A)

V = get_voltage(m)
@info get_voltage V
@test V isa Float64

A = get_current(m)
@info get_current A
@test A isa Float64

Ω = get_resistance(m; wire=2)
@info "get_resistance(; wire=2)" Ω
@test Ω isa Float64

Ω = get_resistance(m; wire=4)
@info "get_resistance(; wire=4)" Ω
@test Ω isa Float64

channel = get_channel(m)
@info get_channel channel
@test channel in ["FRON", "REAR"]

@info "°" get_tc_temperature(m) get_temp_unit(m)

set_temp_unit_kelvin(m)

@info get_tc_temperature(m) get_temp_unit(m)

# set_tc_type(m; type="K")

set_temp_unit_farenheit(m)

@info get_tc_temperature(m) get_temp_unit(m)

set_temp_unit_celsius(m)

@info get_tc_temperature(m) get_temp_unit(m)

terminate(m)


text:
#module Coordinatesm
#export Coordinates, PossibleCoordinates, Coordinates2d, Coordinates2dIndex, polar_to_2d, polar_to_2d_index, calc_possible_coordinates, get_random_coordinates

using Random

struct PossibleCoordinates
    ϕs::Array{Float64}
    θs::Array{Float64}
end

struct Coordinates
    ϕ::Float64
    θ::Float64
end
struct Coordinates2d
    x::Float64
    y::Float64
end
struct Coordinates2dIndex
    x::Int
    y::Int
end
function polar_to_2d(p::Coordinates)
    x=sin(p.θ)*cos(p.ϕ)
    y=sin(p.θ)*sin(p.ϕ)
    return Coordinates2d(x,y)
end
function polar_to_2d_index(p::Coordinates,width::Int,height::Int)
    c=polar_to_2d(p)
    if p.θ < pi/2
        xindex=1+floor(Int,((1.0+c.x)/2.0)*(width/2.0)) #left
    else
        xindex=floor(Int,width/2.0)+1+floor(Int,((1.0+c.x)/2.0)*(width/2.0)) #right
    end
    yindex=1+floor(Int,((1.0+c.y)/2.0)*height)
    return Coordinates2dIndex(xindex,yindex)
end

function calc_possible_coordinates(number_cells::Int)
    nθ=floor(Int,sqrt(number_cells/2))
    nϕ=2*nθ
    ϕs=range(0,2π,length=nϕ)
    θs=range(0,π,length=nθ)
    PossibleCoordinates(ϕs,θs)
end

function get_random_coordinates(possible_coordinates::PossibleCoordinates)
    return Coordinates(rand(possible_coordinates.ϕs),rand(possible_coordinates.θs))
end

function align_to_possible_coordinates(c::Coordinates, possible_coordinates::PossibleCoordinates)
    aligned_ϕ=-1.0
    if c.ϕ < possible_coordinates.ϕs[1]
        aligned_ϕ = possible_coordinates.ϕs[1]
    elseif c.ϕ >= possible_coordinates.ϕs[end]
        aligned_ϕ = possible_coordinates.ϕs[1]
    else
        for ϕ_index in 1:(length(possible_coordinates.ϕs)-1)
            if c.ϕ >= possible_coordinates.ϕs[ϕ_index] && c.ϕ < possible_coordinates.ϕs[ϕ_index+1]
                aligned_ϕ = possible_coordinates.ϕs[ϕ_index]
            end
        end
    end
    aligned_θ=-1.0
    if c.θ < possible_coordinates.θs[1]
        aligned_θ = possible_coordinates.θs[1]
    elseif c.θ >= possible_coordinates.θs[end]
        aligned_θ = possible_coordinates.θs[1]
    else
        for θ_index in 1:(length(possible_coordinates.θs)-1)
            if c.θ >= possible_coordinates.θs[θ_index] && c.θ < possible_coordinates.θs[θ_index+1]
                aligned_θ = possible_coordinates.θs[θ_index]
            end
        end
    end
    return Coordinates(aligned_ϕ,aligned_θ)
end

#end


text:
"""
AstroImage(fits::FITS, ext::Int=1)

Given an open FITS file from the FITSIO library,
load the HDU number `ext` as an AstroImage.
"""
AstroImage(fits::FITS, ext::Int=1, args...; kwargs...) = AstroImage(fits[ext], args...; kwargs...)

"""
AstroImage(hdu::HDU)

Given an open FITS HDU, load it as an AstroImage.
"""
AstroImage(hdu::HDU, args...; kwargs...) = AstroImage(read(hdu), args..., read_header(hdu); kwargs...)

"""
img = AstroImage(filename::AbstractString, ext::Integer=1)

Load an image HDU `ext` from the  FITS file at `filename` as an AstroImage.
"""
function AstroImage(filename::AbstractString, ext::Integer, args...; kwargs...)
    return FITS(filename, "r") do fits
        return AstroImage(fits[ext], args...; kwargs...)
    end
end
"""
img1, img2 = AstroImage(filename::AbstractString, exts)

Load multiple image HDUs `exts` from an FITS file at `filename` as an AstroImage.
`exts` must be a tuple, range, :, or array of Integers.
All listed HDUs in `exts` must be image HDUs or an error will occur.

Example:
```julia
img1, img2 = AstroImage("abc.fits", (1,3)) # loads the first and third HDU as images.
imgs = AstroImage("abc.fits", 1:3) # loads the first three HDUs as images.
imgs = AstroImage("abc.fits", :) # loads all HDUs as images.
```
"""
function AstroImage(filename::AbstractString, exts::Union{NTuple{N,<:Integer},AbstractArray{<:Integer}}, args...; kwargs...) where {N}
    return FITS(filename, "r") do fits
        return map(exts) do ext
            return AstroImage(fits[ext], args...; kwargs...)
        end
    end
end
function AstroImage(filename::AbstractString; kwargs...) where {N}
    return FITS(filename, "r") do fits
        ext = indexer(fits)
        return AstroImage(fits[ext]; kwargs...)
    end
end
function AstroImage(filename::AbstractString, ::Colon, args...; kwargs...) where {N}
    return FITS(filename, "r") do fits
        return map(fits) do hdu
            return AstroImage(hdu, args...; kwargs...)
        end
    end
end

"""
load(fitsfile::String)

Read and return the data from the first ImageHDU in a FITS file
as an AstroImage. If no ImageHDUs are present, an error is returned.

load(fitsfile::String, ext::Int)

Read and return the data from the HDU `ext`. If it is an ImageHDU,
as AstroImage is returned. If it is a TableHDU, a plain Julia
column table is returned.

load(fitsfile::String, :)

Read and return the data from each HDU in an FITS file. ImageHDUs are
returned as AstroImage, and TableHDUs are returned as column tables.

load(fitsfile::String, exts::Union{NTuple, AbstractArray})

Read and return the data from the HDUs given by `exts`. ImageHDUs are
returned as AstroImage, and TableHDUs are returned as column tables.

!!! Currently any header on TableHDUs are not supported and are ignored.
"""
function fileio_load(f::File{format"FITS"}, ext::Union{Int,Nothing}=nothing, args...; kwargs...) where {N}
    return FITS(f.filename, "r") do fits
        if isnothing(ext)
            ext = indexer(fits)
        end
        _loadhdu(fits[ext], args...; kwargs...)
    end
end
function fileio_load(f::File{format"FITS"}, exts::Union{NTuple{N,<:Integer},AbstractArray{<:Integer}}, args...; kwargs...) where {N}
    return FITS(f.filename, "r") do fits
        map(exts) do ext
            _loadhdu(fits[ext], args...; kwargs...)
        end
    end
end
function fileio_load(f::File{format"FITS"}, ::Colon, args...; kwargs...) where {N}
    return FITS(f.filename, "r") do fits
        exts_resolved = 1:length(fits)
        map(exts_resolved) do ext
            _loadhdu(fits[ext], args...; kwargs...)
        end
    end
end

_loadhdu(hdu::FITSIO.TableHDU) = Tables.columntable(hdu)
function _loadhdu(hdu::FITSIO.ImageHDU, args...; kwargs...)
    if size(hdu) != ()
        return AstroImage(hdu, args...; kwargs...)
    else
        # Sometimes files have an empty data HDU that shows up as an image HDU but has headers.
        # Fallback to creating an empty AstroImage with those headers.
        emptydata = fill(missing)
        return AstroImage(emptydata, (), (), read_header(hdu), [emptywcs(emptydata)], Ref(false), ())
    end
end
function indexer(fits::FITS)
    ext = 0
    for (i, hdu) in enumerate(fits)
        if hdu isa ImageHDU && length(size(hdu)) >= 1 # check if Image is atleast 1D
            ext = i
            break
        end
    end
    if ext > 1
        @info "Image was loaded from HDU $ext"
    elseif ext == 0
        error("There are no ImageHDU extensions in '$(fits.filename)'")
    end
    return ext
end
indexer(fits::NTuple{N,FITS}) where {N} = ntuple(i -> indexer(fits[i]), N)

# Fallback for saving arbitrary arrays
function fileio_save(f::File{format"FITS"}, args...)
    return writefits(f.filename, args...)
end
"""
    writefits("abc.fits", img1, img2, table1,...)

Write arguments to a FITS file.

See also `Fileio.save`
"""
function writefits(fname, args...)
    FITS(fname, "w") do fits
        for arg in args
            writearg(fits, arg)
        end
    end
end
writearg(fits, img::AstroImage) = write(fits, parent(img), header=header(img))
# Fallback for writing plain arrays
writearg(fits, arr::AbstractArray) = write(fits, arr)
# For table compatible data.
# This allows users to round trip: dat = load("abc.fits", :); write("abc", dat)
# when it contains FITS tables.
function writearg(fits, table)
    if !Tables.istable(table)
        error("Cannot save argument to FITS file. Value is not an AbstractArray or table.")
    end
    # FITSIO has fairly restrictive input types for writing tables (assertions for documentation only)
    colname_strings = string.(collect(Tables.columnnames(table)))::Vector{String}
    columns = collect(Tables.columns(table))::Vector
    write(
        fits,
        colname_strings,
        columns;
        hdutype=TableHDU
        # TODO: In future, we want to be able to access and round-trip coments
        # on table HDUs
        # header=nothing
    )
end


text:
"""
Quite crude way of checking that a number is a probability

The two steps are

1. The number should be of the `Float64` type
2. The number should belong to [0.0,1.0]
"""
function _value_is_a_probability(p::T) where {T<:AbstractFloat}
	(0.0 <= p <= 1.0) || throw(ArgumentError("The value p=$(p) is not a probability"))
	return nothing
end

"""
Expected value of a single Bernoulli event

Simply f(p): p
"""
function _single_bernoulli_expectation(p::Float64)
    _value_is_a_probability(p)
	return p
end

"""
Variance of a single Bernoulli event

f(p): p(1-p)
"""
function _single_bernoulli_variance(p::Float64)
    _value_is_a_probability(p)
	return p*(1.0-p)
end

"""
Variance of a series of additive Bernoulli events

f(p): ∑(p(1-p))
"""
function _additive_bernoulli_variance(p::Array{Float64})
	return sum(_single_bernoulli_variance.(p))
end

"""
Variance of a series of multiplicative Bernoulli events
"""
function _multiplicative_bernoulli_variance(p::Array{Float64})
	v = _single_bernoulli_variance.(p)
	e = _single_bernoulli_expectation.(p)
	return prod(v.+e.^2.0)-prod(e.^2.0)
end


text:
function onpatch!(x::AbstractVector, v̄::PVector{T,N}) where {T,N}
    ranges = ProjectiveVectors.dimension_indices(v̄)
    for range in ranges
        λ = zero(eltype(x))
        @inbounds for i in range
            λ += v̄[i] * x[i]
        end
        λ⁻¹ = @fastmath inv(λ)
        for i in range
            x[i] *= λ⁻¹
        end
    end
    x
end

function evaluate_patch!(u, v̄::PVector{S,N}, x::PVector{T,N}) where {S,T,N}
    ranges = ProjectiveVectors.dimension_indices(v̄)
    n = length(u) - N
    for (k, range) in enumerate(ranges)
        out = -one(eltype(x))
        for i in range
            out += v̄[i] * x[i]
        end
        u[n+k] = out
    end
    nothing
end

function jacobian_patch!(U, v̄::PVector{S,N}, x::PVector) where {S,T,N}
    ranges = ProjectiveVectors.dimension_indices(v̄)
    n = size(U, 1) - N
    for j = 1:size(U, 2), i = (n+1):size(U, 1)
        U[i, j] = zero(eltype(U))
    end
    for (k, range) in enumerate(ranges)
        for j in range
            U[n+k, j] = v̄[j]
        end
    end
    nothing
end


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

mutable struct IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource <: SwaggerModel
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    partition::Any # spec type: Union{ Nothing, Int32 } # spec name: partition
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    volumeID::Any # spec type: Union{ Nothing, String } # spec name: volumeID

    function IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource(;fsType=nothing, partition=nothing, readOnly=nothing, volumeID=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource, Symbol("partition"), partition)
        setfield!(o, Symbol("partition"), partition)
        validate_property(IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource, Symbol("volumeID"), volumeID)
        setfield!(o, Symbol("volumeID"), volumeID)
        o
    end
end # type IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource

const _property_map_IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource = Dict{Symbol,Symbol}(Symbol("fsType")=>Symbol("fsType"), Symbol("partition")=>Symbol("partition"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("volumeID")=>Symbol("volumeID"))
const _property_types_IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource = Dict{Symbol,String}(Symbol("fsType")=>"String", Symbol("partition")=>"Int32", Symbol("readOnly")=>"Bool", Symbol("volumeID")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource)
    (getproperty(o, Symbol("volumeID")) === nothing) && (return false)
    true
end

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