text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export wfst2tr
"""
  `wfst2tr(fst,Nt; get_label=get_ilabel)`

Returns an array `time2tr` of length `Nt-1`.
The `t`-th element of `time2tr` is a dictionary mapping the transition index between input label at time `t` to the input label at time `t+1`.

```julia
julia> using FiniteStateTransducers

julia> H = WFST(["a1","a2","a3"],["a"]);

julia> add_arc!(H,1,2,"a1","a");

julia> add_arc!(H,2,3,"a2","<eps>");

julia> add_arc!(H,2,2,"a2","<eps>");

julia> add_arc!(H,3,4,"a3","<eps>");

julia> add_arc!(H,3,2,"a3","a");

julia> initial!(H,1);

julia> final!(H,4)
WFST #states: 4, #arcs: 5, #isym: 3, #osym: 1
|1/0.0f0|
a1:a/0.0f0 → (2)
(2)
a2:ϵ/0.0f0 → (3)
a2:ϵ/0.0f0 → (2)
(3)
a3:ϵ/0.0f0 → (4)
a3:a/0.0f0 → (2)
((4/0.0f0))

julia> Nt = 10;

julia> time2tr = wfst2tr(H,Nt)
9-element Array{Dict{Int64,Array{Int64,1}},1}:
 Dict(1 => [2])
 Dict(2 => [2, 3])
 Dict(2 => [2, 3],3 => [2])
 Dict(2 => [2, 3],3 => [2])
 Dict(2 => [2, 3],3 => [2])
 Dict(2 => [2, 3],3 => [2])
 Dict(2 => [2, 3],3 => [2])
 Dict(2 => [2],3 => [2])
 Dict(2 => [3])

```

"""
function wfst2tr(fst::WFST{W,A,D},Nt; get_label=get_ilabel) where {W,A,D}
  isym=get_isym(fst)
  Ns = length(isym)
  X = matrix2wfst(isym,ones(W,Ns,Nt); W = W)

  Z = ∘(X,fst; connect=true)

  i = get_initial(Z)
  current_states = Set{D}([i])
  next_states = Set{D}()
  time2transitions = [Dict{D,Vector{D}}() for t in 1:Nt-1]
  t = 1
  while !isempty(current_states)
    dict = Dict{D,Set{D}}()
    for c in current_states
      for c_arc in Z[c]
        src = D(get_label(c_arc))
        n = get_nextstate(c_arc)
        push!(next_states,n)
        dest = Set{D}([D(get_label(n_arc)) for n_arc in Z[n]])
        if !isempty(dest)
          if haskey(dict,src)
            union!(dict[src],dest)
          else
            dict[src] = dest
          end
        end
      end
    end
    if t < Nt
      time2transitions[t] = Dict(k=>[dict[k]...] for k in keys(dict))
    else
      if !isempty(dict)
        error("something is wrong with the graph")
      end
    end
    t += 1
    current_states, next_states = next_states, Set{D}()
  end
  return time2transitions
end

"""
  `wfst2tr(fst; get_label=get_ilabel, convert_weight=identity)`

Extract transition matrix and initial probabilities from a WFST representing an HMM model.

```julia
julia> add_arc!(H,1,2,"a1","a");

julia> add_arc!(H,2,2,"a1","<eps>",-log(0.3));

julia> add_arc!(H,2,3,"a2","<eps>",-log(0.7));

julia> add_arc!(H,3,3,"a2","<eps>",-log(0.3));

julia> add_arc!(H,3,4,"a3","<eps>",-log(0.7));

julia> add_arc!(H,4,4,"a3","<eps>",-log(0.3));

julia> add_arc!(H,4,2,"a1","a",-log(0.7));

julia> initial!(H,1);

julia> final!(H,4)
WFST #states: 4, #arcs: 12, #isym: 3, #osym: 1
|1/0.0f0|
a1:a/0.0f0 → (2)
a1:a/0.0f0 → (2)
(2)
a2:ϵ/0.0f0 → (3)
a2:ϵ/0.0f0 → (2)
a1:ϵ/1.2039728f0 → (2)
a2:ϵ/0.35667494f0 → (3)
(3)
a3:ϵ/0.0f0 → (4)
a3:a/0.0f0 → (2)
a2:ϵ/1.2039728f0 → (3)
a3:ϵ/0.35667494f0 → (4)
((4/0.0f0))
a3:ϵ/1.2039728f0 → (4)
a1:a/0.35667494f0 → (2)

julia> a,A = wfst2tr(H; convert_weight = w -> exp(-get(w)))
(Float32[1.0, 0.0, 0.0], Float32[0.3 0.7 0.0; 0.0 0.3 0.7; 0.7 0.0 0.3])
```
"""
function wfst2tr(fst::WFST{W,A,D};
                 get_label=get_ilabel, convert_weight=identity) where {W,A,D}
  isym=get_isym(fst)
  Ns = length(isym)
  W2 = typeof(convert_weight(one(W)))
  a = zeros(W2,Ns)
  H = zeros(W2,Ns,Ns)
  init = get_initial(fst)
  state2outtr = Dict(i => (get_ilabel.(s),get_weight.(s)) for (i,s) in enumerate(fst))
  for (p,s,n,d,e,a) in FiniteStateTransducers.DFS(fst,init)
    if d
      intr = get_ilabel(a)
      outtr, w = state2outtr[n]
      for i in eachindex(outtr)
        H[intr,outtr[i]] = convert_weight(w[i])
      end
    end
  end
  for arc in fst[init]
    a[get_label(arc)] = convert_weight(get_weight(arc))
  end
  return a, H
end


text:
module timeevolution_simple

using ..states
using ..operators
using ODE

export master

"""
Evaluate master equation for diagonal jump operators.
"""
function dmaster(rho::DenseOperator, H::Operator, gamma::Vector, J::Vector, Jdagger::Vector)
    drho = -1im * (H*rho - rho*H)
    for n = 1:length(J)
        drho = drho + gamma[n]*(J[n]*rho*Jdagger[n] - Jdagger[n]*(J[n]*rho)/Complex(2) - rho*Jdagger[n]*J[n]/Complex(2))
    end
    return drho
end

"""
Evaluate master equation for non-diagonal jump operators.
"""
function dmaster(rho::DenseOperator, H::Operator, gamma::Matrix, J::Vector, Jdagger::Vector)
    drho = -1im * (H*rho - rho*H)
    for m=1:length(J), n=1:length(J)
       drho += gamma[m,n]*(J[m]*rho*Jdagger[n] - Jdagger[n]*(J[m]*rho)/Complex(2) - rho*Jdagger[n]*J[m]/Complex(2))
    end
    return drho
end

"""
Integrate master equation.
"""
function master(T::Vector, rho0::DenseOperator, H::Operator, J::Vector;
                    Jdagger=map(dagger,J),
                    gamma::Union{Real, Vector, Matrix}=ones(Int, length(J)),
                    kwargs...)
    operators.check_samebases(rho0, H)
    for j=J
        @assert typeof(j) <: Operator
        operators.check_samebases(rho0, j)
    end
    for j=Jdagger
        @assert typeof(j) <: Operator
        operators.check_samebases(rho0, j)
    end
    @assert length(J)==length(Jdagger)
    if typeof(gamma)<:Real
        gamma = ones(typeof(gamma), length(J))*gamma
    end
    nl = prod(rho0.basis_l.shape)
    nr = prod(rho0.basis_r.shape)
    N = nl*nr
    as_operator(x::Vector{Complex128}) = DenseOperator(rho0.basis_l, rho0.basis_r, reshape(x, nl, nr))
    as_vector(rho::DenseOperator) = reshape(rho.data, N)
    f(t::Float64, x::Vector{Complex128}) = as_vector(dmaster(as_operator(x), H, gamma, J, Jdagger))
    tout, x_t = ode45(f, as_vector(rho0), T; kwargs...)
    rho_t = DenseOperator[as_operator(x) for x=x_t]
    return tout, rho_t
end

master(T::Vector, psi0::Ket, H::Operator, J::Vector; kwargs...) = master(T, tensor(psi0, dagger(psi0)), H, J; kwargs...)

"""
Evaluate Schroedinger equation for ket states.
"""
function dschroedinger(psi::Ket, H::Operator)
    return -1im*H*psi
end

"""
Evaluate Schroedinger equation for bra states.
"""
function dschroedinger(psi::Bra, H::Operator)
    return 1im*psi*H
end

"""
Integrate Schroedinger equation.
"""
function schroedinger{T<:StateVector}(tspan::Vector, psi0::T, H::Operator; kwargs...)
    as_statevector(x::Vector{Complex128}) = T(psi0.basis, x)
    as_vector(psi::T) = psi.data
    f(t::Float64, x::Vector{Complex128}) = as_vector(dschroedinger(as_statevector(x), H))
    tout, x_t = ode45(f, as_vector(psi0), tspan; kwargs...)
    psi_t = T[as_statevector(x) for x=x_t]
    return tout, psi_t
end

end  # module

ext: jl


text:
module CUBLAS

using ..CUDA
import ..CUDA: ndevices, getdevice
import Libdl

if Sys.iswindows()
    const libcublas = Libdl.find_library(["cublas64_92","cublas64_91","cublas64_90","cublas64_80","cublas64_75"])
else
    const libcublas = Libdl.find_library("libcublas")
end
isempty(libcublas) && error("CUBLAS cannot be found.")

const API_VERSION = Ref{Int}()

function __init__()
    ref = Ref{Ptr{Cvoid}}()
    ccall((:cublasCreate_v2,libcublas), Cint, (Ptr{Ptr{Cvoid}},), ref)
    h = ref[]

    ref = Ref{Cint}()
    ccall((:cublasGetVersion_v2,libcublas), Cint, (Ptr{Cvoid},Ptr{Cint}), h, ref)
    API_VERSION[] = Int(ref[])

    @info "CUBLAS API $(API_VERSION[])"
    ccall((:cublasDestroy_v2,libcublas), Cint, (Ptr{Cvoid},), h)
end

include("define.jl")

macro cublas(f, rettypes, args...)
    f = get(DEFINE, f.value, f.value)
    quote
        status = ccall(($(QuoteNode(f)),libcublas), Cint, $(esc(rettypes)), $(map(esc,args)...))
        if status != 0
            throw(ERROR_MESSAGE[status])
        end
    end
end

function cublasop(t::Char)
    t == 'N' && return Cint(0)
    t == 'T' && return Cint(1)
    t == 'C' && return Cint(2)
    throw("Unknown cublas operation: $(t).")
end

include("handle.jl")
include("level1.jl")
include("level2.jl")
include("level3.jl")
include("extension.jl")

end


text:
using BenchmarkTools
#global ⫙ = []
#=

#another idea: pass in middle function, as a whole
#function middleExtraction(middle,a,b) end
#---testing area 2
=#
"""returns a boolean
```input:

```

```output:

```
"""#testing purposes only
#= #Highly Experimental
function checkCriterion(ch, :operator = [%, +, -, *], operand)
    check = nothing
    check = ch
    condition = nothing
    condition = :check:operator:operand
    if condition == true
        return true #a.s.
    elseif condition == false
        return false #a.s.
    else # faulty Input or Unexpected Error Occured
        return check
    end
    return check
end
checkCriterion(3, %, 2)
 =#
#check #1 ϟ - \upkoppa
function ϟ(a, b)
    return abs(a) + abs(b)
end

#check #2:middle
#
"""evenCriterion
```input
  a: the `first()` in a range
  b:the `last()` in  a range
```
```output
 Returns `check`: the

 Output is Calculated by: `check b-a`

 `output = `
   If check is true

  If check if false

  Else, there must be an error in either input arguments: a, b (or an Unexpected error occured)
  Otherwise,  return the `check`

  **see also:** `iseven`

  ```
"""
#isEven Removed #Reason: unnecessary Overhead
#define sumInterval here
sumInterval(a::Int64, b::Int64) = (abs(max(a, b)) - abs(min(a, b))) # was (max(a, b) - abs(min(a, b))) #Unsymmetric
#=
function evenCriterion(a, b) #ok #double-Checked # depreciated # erroneous
    m = ϟ(a, b) # | b + a |   definition
    condition = m % 2 == 0 #isEven(m) #even (divisible by 2)
    check = nothing
    if condition == true  #ERROR double if!
        check = condition
    elseif condition == false
        check = condition
    else #if faulty input or Unexpected ERROR Occured
        check
    end
    return check #whether check is true, false, nothing
end
=#
#done!

res = evenCriterion(1, 3) #flase Infer no middle as a whole number
#=
function midCriterion(sumInterval)
    m = sumInterval #NO ERROR  #sumInterval(a, b) # | b - a |
    condition = m % 2 == 0 # isEven(m)
    check = nothing
    if condition == true
        check = condition
    elseif condition == false
        check = condition
    else #if faulty input or Unexpected ERROR Occured
        check
    end
    return check #whether check is true, false, nothing
end
=#
#--------------
#Recheck - cyclical code dedected! :mid to even , or even to mid  #Possible #depreciatable
function middleExtraction(condition, checkedValue; above = nothing, below = nothing) #at least one argument is given, besides condition
    ⫙ = []# reset is necessary detected
    criterion = condition
    if criterion == true # there is only 1 value
        # only get the checkedValue
        ⫙ = push!(⫙, checkedValue)
    elseif criterion == false
        #then we have above & below to get (know for sure that below is under above)
        ⫙ = push!(⫙, above)
        ⫙ = push!(⫙, below)
    else # faulty input or Unexpected error Occured
        return
    end
    return ⫙ #q
end
#---------------------------------------
function middleGet(condition)
    criterion = condition
    criterion ? ⫙ = (popat(checkedValue)) : ⫙ = popat(ceil)
    val = popat(below) #warning: unused  ⫙
    return val
end

function doublepush!(m = ⫙, above, below)
    push!(m, above)
    push!(m, below)
end

function middleSet(condition; m = ⫙, checkedValue, above, below)
    criterion = condition
    criterion ? push!(m, checkedValue) : doublepush(m, above, below) #push!(⫙, ceil)

end
arr = [10, 4, 2, 8]
middleGet(ϟ(first(arr), last(arr)))
#----------------------------------------
#middleExtraction()
check = Int(ϟ(a, b) // 2) # sumInterval(a, b) // 2 #* 1.0 # | b - a | // 2 isa Integer #sumInterval -to-> ϟ
#⫙ = []

function push!(⫙, ceil = above, floor = below)#warning: ⫙, floor are Unused
    push!(⫙, ceil)
    push!(⫙, below)
end
""" pops Double  """
function popDouble()
    ceil = popat(⫙, ceil)
    floor = popat(⫙, below)
    return ceil, floor
end
"""pushes ceil & floor """
function pushDouble(⫙, ceil = above, floor = below)
    ceil = push!(⫙, ceil)
    floor = push!(⫙, below)

    return ceil, floor
end
#=
"""pushes checked """
function push(checked = checkedValue)
    push!(⫙, checked)
end

function pop!()
    checked = pop!()
    return checked
end

#function pop!()
=#
#--------------------------
#---known error here
"""returns a middle(s)
```input:
a: lower bound range
b: upper bound range
```

```output:
condition: Bool flag informing user if got a 1 mid-point, 2 mid-poimts, or errors occured
check: 1 midpoint middle  (if condition is true)
above,below: 2 mid-point middles (the realized leaf-bounds between the theoretical mid-point)
```
"""
function middle(a, b) # working
    condition = evenCriterion(a, b) #   b-a
    # check = nothing;check = midCriterion(m) #b-a
    above = nothing
    below = nothing
    check = nothing
    q = []
    if condition == true
        #return true #a.s. #eucledian Distance divided by 2 returing a whole integer
        check = Int(ϟ(a, b) // 2) # | b + a | // 2 isa Integer #sumInterval -to-> ϟ#5
        #middleExtraction(condition, check) # Here we didn't get anything ! <------------- # check not defined here
        #return condition, check
        push!(q, check)
    elseif condition == false
        #return false #a.s.# check = sumInterval(a,b)//2*1.0
        #GET Ceil & Floor
        check = ϟ(a, b) / 2 # floating-point division sumInterval(a, b) / 2 * 1.0 # freely allowing floats, to be ceiled & floored
        above = Int(ceil(check)) #nearest index above
        below = Int(floor(check))
        push!(q, below)
        push!(q, above)

    else # faulty Input or Unexpected Error Occured
        #    return check  # nothing
        return q #condition, check, above, below
    end
    return q #condition, check, above, below
end

condition = nothing;
check = nothing;
# condition, check, above, below = middle(1, 3)
q = middle(1, 3)
typeof(q)
Q = zeros(length(q) + 1) # for future testing purposes
#Note: there should be a math formula to translate points into range (numbers)
l = length(q)      # <--- length of q mid-point(s)
#=
while  l >= 1

end
=#
val = popfirst!(q) #can only observe value once consumed & popped out

"""
gets number of total points, based on number of mid-point(s)

```input:
l: number of mid-points , a Julia's length function length(middle(lowerbound, upperbound) )
; modulo: (Optional) the modulo ,defaults to 2 (otherwise, untested )
```

```output:
if even, returns 4 (sumPts * _two (start ))
```
"""
function numPts(Length = l; modulo = 2)
    _two = modulo
    sumPts = _two
    if Length % _two == 0 # even
        sumPts = sumPts + _two # multiplier * 2 i.e. 2 * 2 = 4
    elseif Length % _two != 0
        sumPts = sumPts + (_two / _two)     #multiplier + 1 i.e. 2 + 1 = 3
    else
        return
    end
    return Int(sumPts)
end
q = middle(1, 10)
Length = length(q)
nPoints = numPts(Length)

using BenchmarkTools
@benchmark n = numPts(l)

numPts()

#-risk-free experiments below
lenpts = l * 2 + l # if l== 1 (1 midpoint), apply condition that you should add a + 1  (mid point)  #Got inspired the numPts function
for i ∈ (q)
    #   for j ∈(Q)
    i = q.pop()

end
ranges = []
#if condition == true
push!(ranges, 1:check)
push!(ranges, check:3)

popfirst!(ranges)

#bottleneck is the first middle
arr = [8, 5, 2, 7, 4]
q = []
push!(q, middle(arr[first(arr)], arr[last(arr)]))
#q = push!( middle(arr[first(arr)], arr[(last(arr))])) # find middle points #funto-watch!
q

function middle(arr, a, b)
    #ϟ(arr)
    q = []
    push!(q, middle(arr[a], arr[b])) # find middle point(s)
    return q
end
#done

#=
if condition == true # 1 midpoint
    #oldschoolComp()
    return #check
    #build range
if condition == false # 2 midpouts #done
    #above & below # compare them here
    below,above =  oldschoolComp(below,above)
    return below,above #<--- pplace here: #compare
     #end
end
=#
# comp(arr,a,b) = arr_a =  arr[a]; arr[b]  = arr_b; arr_a>arr_b ? oldschoolSwap(arr[a],arr[b])

function comp(arr, a, b)
    arr_a = arr[a]
    arr_b = arr[b]
    if arr_a > arr_b
        arr_a, arr_b = oldschoolSwap(arr[a], arr[b])
    end
    return arr_a, arr_b
end

#--------Demonstration

middle(1, 4) #false # (false, 2.5, 3, 2)
middle(1, 3) #true  #
#=sumInterval(1, 3) # (true, 2//1, nothing, nothing)=#

#---testing of testing
middle(1, 3)

"""
Assumes right most argument should be the lowest
Thus, if it has a higher value, swap it (with the other one)
```input:
a: Int(64,128) or float(64, or 128)
b: Int(64,128) or float(64, or 128)
```

```output:
returns an ordered pair of a & b
(comparison might be needed only if a is larger)
```
"""
function oldschoolComp(a, b)
    if a > b
        a, b = oldschoolswap(a, b)
        return a, b
    elseif a < b
        # move on
        return a, b
    else
        return 0, 0
    end
    return a, b
end
elementat() #returns element
indexOf()#returns index

#----new: Experimental------
isDiff(a, b, aVal, bVal) = aVal != a && bVal != b ? true : false
function oldschoolComp(arr, a = first(arr), b = last(arr))
    aVal = a = first(arr)
    bVal = b = last(arr)

    aVal, bVal = oldschoolComp(a, b)
    arr[aVal]
    ans = isDiff(a, b, aVal, bVal) #true
    ans == true ? arr[indexOf(a, arr)], arr[indexOf(b, arr)] = arr[indexOf(b, arr)], arr[indexOf(a, arr)] : arr
    arr[]

    if a > b
        a, b = oldschoolswap(a, b)
        return a, b
    elseif a < b
        # move on
        return a, b
    else
        return 0, 0
    end
    return a, b
end

aVal = first(arr)
bVal = last(arr)

aVal, bVal = oldschoolComp(a, b)

arr[aVal]

ans = isDiff(a, b, aVal, bVal) #true
ans == true ? arr[indexOf(a, arr)], arr[indexOf(b, arr)] = arr[indexOf(b, arr)], arr[indexOf(a, arr)] : """
                                                                                compares an arbitrary array
                                                                                based on it's first & last elements
                                                                                """
#=
function compare(arr) #errornous #warning becareful from its values
    aVal = first(arr) # arr[aIdx] # 5 v[1]= 5
    bVal = last(arr) #arr[bIdx] # v[last] = 1  5>1
    aIdx = firstindex(arr)
    bIdx = lastindex(arr) # bVal)

    #1. Heuristic#1: Check sumInterval

    total = ϟ(aIdx,bIdx)
    dist = sumInterval(aIdx,bIdx)
    midpoint

    if dist == 1 #1,2 ; 2,3 ; 3,4
       #final compare of values
      aVal,bVal =  oldschoolComp(aVal,bVal)#compare & swap, if applicable

    elseif dist > 1
        isEven = evenCriterion(aIdx, bIdx)
        if isEven == true #even #only 1-midpoint
            midpoint = pushfirst!(midpoint, Int(total // 2))
        elseif isEven == false # 2 midpoints
        #middle()

        end
        #if aVal < bVal
        if total > 1 # there are still unforseen ranges
            #if sumInterval(aIdx,bIdx) > 1 #reason to stay  #\upkoppa ?
            #if aVal > bVal compare function                               # a > b #should be a < b
            #1. Swap values
            aVal, bVal = comp(arr, aIdx, bIdx)  # oldschoolSwap(aVal, bVal) #correct #x, y) # (bVal, aVal) #ERROR T not defined
            #2. Swap  array indicies(with correct orders)
            #arr[aIdx], arr[bIdx] = oldschoolSwap(arr[aIdx], arr[bIdx])
            #3.decrement  - middle
            # middleExtraction() #error! #UncommentMe
            #elseif aVal < bVal # at that level (diff length) of indicies # the array is correct

            return aVal, bVal

        else #faulty input , Unexpected code
            println(strError)
            return
        end
        return aVal, bVal # , arr  #min, max #removed values
    end
end
=#

arr = [8, 5, 2, 7, 4]
#finished leaf example
aVal = first(arr)
bVal = last(arr)
aVal, bVal = oldschoolComp(aVal, bVal)

q = middle(first(arr), last(arr))

n = numPts(length(q)) # 3 middles

q = middle(first(arr), last(arr))
l = length(q) # length 1 # errornous process

ans = popfirst!(q)

count = 1
val = 2
#---new
if count <= n
    q = middle(first(arr), last(arr)) #
    val = firstindex(q) # findall(x -> x == val, arr) #popfirst!(q)
end
#---end

#method at least julia 1.7
#i = find(x -> x == val, arr) #removed

arr = [8, 5, 2, 7, 4]
elementat(arr, 8)
res = indexOf(1, arr)
_last = arr[indexOf(last(arr), arr)]
_first = arr[indexOf(first(arr), arr)]
ans = indexOf(arr, elementat(arr, 5)[1])

indexOf(elementat(7, arr), arr)

length(elementat(1, arr))
using BenchmarkTools
@benchmark last(arr)
#=
 Range (min … max):  18.236 ns … 182.966 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     19.639 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   23.044 ns ±  12.668 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▄█▄          ▁                                               ▁
  ████▅▅▅▄▄▃▂▂▅██▇█████▇▆▄▄▅▇▆▆▅▅▇▆▆▇▆▆▆▅▄▅▄▄▄▃▄▄▂▄▄▃▃▃▃▃▃▂▄▄▃ █
  18.2 ns       Histogram: log(frequency) by time      83.8 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.
=#
length(arr)

if _first < _last
    arr[1], arr[length(arr)] = oldschoolComp(_first, _last)
end

@inbounds arr[1], arr[length(arr)] = arr[length(arr)], arr[1] #

function Swap(arr, a = arr[1], b = arr[length(arr)])
    return @inbounds arr[1], arr[length(arr)] = arr[length(arr)], arr[1]

end
#-------------------
#arr[1], arr[last(arr)] = oldschoolComp(arr[1], arr[last(arr)])

arr[1], arr[length(arr)] = comp(arr, 1, length(arr))
q = middle(1,length(arr))
newMid = q[1]
elementat(arr,newMid) #should be 2
arr[newMid]
arr
res = Int(ϟ(1, length(arr)) / 2)
res % 2 == 0 ? true : false #false - 2 midpoints
#Ceil =
#    floor =
check = res / 2 # floating-point division sumInterval(a, b) / 2 * 1.0 # freely allowing floats, to be ceiled & floored
above = Int(ceil(check)) #nearest index above
below = Int(floor(check))
arr[below]
arr[above]
elementat(arr,arr[below])
elementat(arr,arr[above])

arr
#=
1:2
2:
=#
arr[1+2], arr[length(arr)-1] = comp(arr, 1 + 1, length(arr) - 1)
arr[1+2], arr[length(arr)-1] = comp(arr, 1 + 1, length(arr) - 1)

arr

#----testing Area
arr[1], arr[last(arr)]

val = 2
ans = firstindex(q) #true
i = elementat(arr, 2) #
#i[1] #
ans = elementat(arr, 2)
ans # index

ans = indexOf(arr, i)

elementat(arr, 5) # working

#firstIndex()
#Q.how far 3 from start?
first(arr) #first val
ans = indexin(i[1], arr)

d1 = sumInterval(i[1], first(arr))

#valIdx = findall(arr, val)

val

#correct till this point
function goabove!(arr, middle, last) end

function gobelow!(arr, first, middle)
    #evaluate value

end

#-----

arr = [8, 5, 2, 7, 4]

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
# using GaussianMixtures
# using Distributions
# using PDMats
# import Base.LinAlg.Cholesky
#
# mutable struct SingularGMM
#     K::Int
#     D::Int
#     μ::Matrix{Float64}
#     Σ::Vector{Matrix{Float64}}
#     w::Vector{Float64}
#     SingularGMM(K,D) = new(K,D,zeros(K,D),Matrix{Float64}[zeros(D,D) for i=1:K], zeros(K))
# end
#
# mutable struct LinearDynamics
#     μ::Vector{Float64}
#     Σ::Matrix{Float64}
#     fx::Matrix{Float64}
#     fu::Matrix{Float64}
# end
# LinearDynamics(d::Int) = LinearDynamics(zeros(d),zeros(d,d),zeros(d,d),zeros(d,d))
#
# mutable struct AffineTransform
#     U::Matrix{Float64}
#     xmin::Vector{Float64}
#     xmax::Vector{Float64}
# end
#
#
# # """
# # Transforms a ´point´ to the ´Lower´ dimensional space
# # """
# function downtransform(T::AffineTransform, p)
#     p .-= T.xmin
#     p ./= T.xmax
#     p   = T.U'*p
# end
#
#
# """
# Transforms a ´point´ to the ´Higher´ dimensional space
# """
# function uptransform(T::AffineTransform, p)
#     p   = T.U*p
#     p .*= T.xmax
#     p .+= T.xmin
# end
#
# """
# `strip(Tin::AffineTransform, i)`
#
# Keep only dimensions in the index vector `i`
# """
# function Base.strip(Tin::AffineTransform, i)
#     T      = deepcopy(Tin)
#     T.U    = T.U[i,:]
#     T.xmax = T.xmax[i]
#     T.xmin = T.xmin[i]
#     return T
# end
#
# """
# `model = fit_good_model(X, K, nTries)`
#
# This function tries to fit a GMM with `K` clusters a number of times and selects the best model
# """
# function fit_good_model(X, K, nTries)
#     models = Vector{GMM}(undef,nTries)
#     for i = 1:nTries
#         itworked = 0
#         while  10 > itworked >= 0
#             try
#                 models[i] = GMM(K, X'; method=:kmeans, kind=:full, nInit=50, nIter=100, nFinal=100)
#                 itworked = -1
#             catch ex
#                 itworked += 1
#                 println(ex)
#             end
#         end
#         itworked == 10 && error("Better luck another time")
#     end
#     model = models[argmax([avll(models[i],X') for i = 1:nTries])]
# end
#
# """
# Get assignments for the data in `X` under the `model`
#
# `Q,ass = get_assignments(model::GMM, X)`
# """
# function get_assignments(model::GMM, X)
#     N = size(X,2)
#     K = model.n
#     Q = zeros(N,K)
#     Σ = covars(model)
#     for n = 1:N
#         for k = 1:K
#             #             Γ = PDMat(Cholesky(model.Σ[k].data,:U))
#             μ = vec(model.μ[k,:]')
#             Q[n,k] = model.w[k]*pdf(MvNormal(μ, Σ[k]),vec(X[:,n]))
#         end
#         Q[n,:] ./= sum(Q[n,:])
#     end
#     ass = findmax(Q,2)[2][:]
#     ass = [ind2sub(Q,a)[2] for a in ass][:]
#     if N >= K
#         1:K ⊆ ass || @warn("Not all clusters have data points assigned to them")
#     end
#     return Q,ass
# end
#
# # """
# # Get assignments for the data in `X` under the `model`
# #
# # `ass = get_assignments_robust(model::GMM, X)`
# # """
# function get_assignments_robust(model::GMM, X)
#     N = size(X,2)
#     K = model.n
#     Q = zeros(K)
#     Σ = covars(model)
#     ass = zeros(Int, N)
#     for n = 1:N
#         for k = 1:K
#             μ = vec(model.μ[k,:]')
#             Q[k] = model.w[k]*logpdf(MvNormal(μ, Σ[k]),vec(X[:,n]))
#         end
#         ass[n] = argmax(Q)
#     end
#     return ass
# end
#
# function get_assignments(model::GMM, X̄, T)
#     X = downtransform(T, X̄)
#     get_assignments(model::GMM, X)
# end
#
# function get_assignments_robust(model::GMM, X̄, T)
#     X = downtransform(T, X̄)
#     get_assignments_robust(model::GMM, X)
# end
#
#
# function conditionalGMM(M, d, p, T)
#     K,D  = M.K, M.D
#     r    = rank(M.Σ[1])
#     Sf   = 1:D
#     Sl   = d
#     Sr   = setdiff(Sf,Sl)
#     Tr   = deepcopy(T)
#     Tr.U = Tr.U[Sr,:]
#     Tl   = deepcopy(T)
#     Tl.U = Tl.U[Sl,:]
#     pr   = downtransform(Tr,p)
#     Mc   = GMM(K,length(d))
#     Mc.w = M.w
#     J    = Vector{Matrix{Float64}}(undef,K)
#     Q    = Vector{Float64}(undef,K)
#     f    = Vector{Float64}(undef,K)
#     for k = 1:K
#         Σ          = M.Σ[k]
#         μ          = M.μ[k,:]'
#         J[k]       = Σ[Sl,Sr]/Σ[Sr,Sr]
#         Mc.μ[k,:]  = (μ[Sl] + J[k]*(pr - μ[Sr]))'
#         Mc.Σ[k]    =  Σ[Sl,Sl] - J[k]*Σ[Sr,Sl]
#         # J[k]     = Tl.U*A*Tr.U'
#         @warn("Have not yet considered scaling by Xmax")
#         f[k] = fN(vec(M.μ[k,Sr]), M.Σ[k][Sr,Sr], pr, regul)[1]
#         Q[k] = M.w[k]*f[k]
#     end
#
#     Q ./= sum(Q)
#
#     return Mc, J
# end
#
#
# # """
# # calculates the model p(x^+ | x, u)
# # where x^+ is p[d] and [x;u] is p[!d]
# # dof is the numer of degrees of freedom, like 7 for the YuMi robot
# # """
#
# function conditional_dynamics(M, d, p, dof, T)
#     (K,D) = (M.K, M.D)
#     Mc, Jk= conditionalGMM(M, xp_indices, p, T)
#     dyn   = LinearDynamics(length(d))
#
#     J   = zeros(length(Sl),length(Sr))
#     for k = 1:K
#         dyn.μ += Q[k]*Mc.μ[k,:]
#         dyn.Σ += Q[k]*Mc.Σ[k]
#         J     += Q[k]*Jk[k]
#     end
#
#     dyn.fx   = J[:,1:dof]
#     dyn.fu   = J[:,dof+1:end]
#     return dyn
# end
#
#
# function normalize_data(x, tol=0.1)
#     mi = mean(x,1)
#     x .-= mi
#     ma = (maximum(x,1)+tol)
#     x ./= ma
#     return x, mi, ma
# end
#
# """
# `Λ(A,r::Int)` This function takes an input matrix A and regularizes it so that all singular values σi=σr ∀ i ≧ r = rank(A)
# The idea is that singular covariance matrices can be regularized so that they can be inverted. This has the effect of adding some small amount of uncertainty (variance) in the directions that have zero (very small) variance
# """
# function Λsvd(A,r::Int)
#     (U,S,V) = svd(A)
#     if r < length(S)
#         S[r+1:end] = S[r]
#     end
#     return U*diagm(S)*V'
# end
#
# """
# `Λ(A,r::Int)` This function takes an input matrix A and regularizes it so that all eigen values λi=λr ∀ i ≧ r = rank(A)
# The idea is that singular covariance matrices can be regularized so that they can be inverted. This has the effect of adding some small amount of uncertainty (variance) in the directions that have zero (very small) variance
# """
# function Λ(A,r::Int)
#     e = eigfact((A+A')/2)
#     if r < length(e.values)
#         e.values[1:end-r] = e.values[end-r+1]
#     end
#     return e.vectors*diagm(e.values)*e.vectors'
# end
#
# """
# `Λ(A,r::Float64=1e-5)` This function takes an input matrix A and regularizes it by adding a small multiple of the identity matrix
# The idea is that singular covariance matrices can be regularized so that they can be inverted. This has the effect of adding some small amount of uncertainty (variance) in the directions that have zero (very small) variance
# """
# function Λ(A,r::Float64=1e-5)
#     return A + r*I
# end
#
#
# """
# fN(M::MvNormal, x) -> (f,df) calculates the pdf and the derivative of the pdf of the MvNormal M at point x
# """
# function fN(μ, Σ, x, r=10_000)
#     k   = size(x,1)
#     d   = x-μ
#     Σid = Σ\d
#     e   = exp(-0.5*d'*Σid)[1]
#     f   = sqrt((2*pi)^k*det(Σ))*e
#     df  = 2Σid*f
#     return f,df
# end
#
#
# """
# `dyn = fitDynamics_decoupled(y, x, u; λ = 1e-7)`
#
# This function assumes that there are an equal number of states and control signals, and that the modeling can be decoupled such that each state is affected by a single input.
# """
# function fit_dynamics_decoupled(y, x, u; λ = 1e-7)
#     N,d = size(x)
#     w   = zeros(2d,d)
#     n   = d + size(u,2)
#     for i = 1:d
#         A     = [x[:,i] u[:,i]]
#         G     = [A; λ*I]
#         w[[i,i+d],i]     = G\[y[:,i]; zeros(n)]
#     end
#     ŷ  = A*w
#     e  = y-ŷ
#     # σ2 = var(e)
#     Σ  = cov(e)
#     fx = w[1:d,:]
#     fu = w[d+1:2d,:]
#     #     μ  = w[2d+1:end,:]
#     return LinearDynamics(vec(mean(ŷ,dims=1)), Σ, fx, fu)
# end
#
# # """
# # `dyn = fitDynamics(y, x, u; λ = 1e-7)`
# #
# # """
# function fit_dynamics(y, x, u; λ = 1e-7)
#     N,d  = size(x)
#     A    = [x u]
#     n    = size(A,2)
#     G    = [A; λ*eye(n)]
#     w    = G\[y; zeros(n,d)]
#     ŷ    = A*w
#     e    = y-ŷ
#     # σ2 = var(e)
#     Σ    = cov(e)
#     fx   = w[1:d,:]'
#     fu   = w[d+1:end,:]'
#     # μ  = w[2d+1:end,:]
#
#     return LinearDynamics(vec(mean(ŷ,dims=1)), Σ, fx, fu)
# end


text:
using BenchmarkTools
using QMTK

macro section(name, expr)
    println("[Section]: $name")
    return expr
end

macro median(name, expr)
    return :(println("[$($name)]: $(median($expr).time) ns (median)"))
end

@section "Chain Lattice" begin

length = 1000

cache = 0
baseline = @benchmark for i=1:length
    i
end

chain = Chain{Fixed}(length)
cache = 0
results = @benchmark for each in sites(chain)
    each
end

@median "Baseline.sites" baseline
@median "QMTK.sites" results

baseline = @benchmark for i=1:length
    i, i+1
end

results = @benchmark for each in bonds(chain, 1)
    each
end

@median "Baseline.bonds" baseline
@median "QMTK.bonds" results

end

@section "Square Lattice" begin

shape = (30, 30)

baseline = @benchmark for i=1:shape[1], j=1:shape[2]
    i, j
end

square = Square{Fixed}(shape)
results = @benchmark for each in sites(square)
    each
end

@median "Baseline.sites" baseline
@median "QMTK.sites" results

results = @benchmark [each for each in bonds(square, 1)]

# TODO: add baseline
# @median "Baseline.bonds" baseline
@median "QMTK.bonds" results

end


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

@doc raw"""Probe describes a health check to be performed against a container to determine whether it is alive or ready to receive traffic.

    IoK8sApiCoreV1Probe(;
        exec=nothing,
        failureThreshold=nothing,
        httpGet=nothing,
        initialDelaySeconds=nothing,
        periodSeconds=nothing,
        successThreshold=nothing,
        tcpSocket=nothing,
        timeoutSeconds=nothing,
    )

    - exec::IoK8sApiCoreV1ExecAction : One and only one of the following should be specified. Exec specifies the action to take.
    - failureThreshold::Int32 : Minimum consecutive failures for the probe to be considered failed after having succeeded. Defaults to 3. Minimum value is 1.
    - httpGet::IoK8sApiCoreV1HTTPGetAction : HTTPGet specifies the http request to perform.
    - initialDelaySeconds::Int32 : Number of seconds after the container has started before liveness probes are initiated. More info: https://kubernetes.io/docs/concepts/workloads/pods/pod-lifecycle#container-probes
    - periodSeconds::Int32 : How often (in seconds) to perform the probe. Default to 10 seconds. Minimum value is 1.
    - successThreshold::Int32 : Minimum consecutive successes for the probe to be considered successful after having failed. Defaults to 1. Must be 1 for liveness and startup. Minimum value is 1.
    - tcpSocket::IoK8sApiCoreV1TCPSocketAction : TCPSocket specifies an action involving a TCP port. TCP hooks not yet supported
    - timeoutSeconds::Int32 : Number of seconds after which the probe times out. Defaults to 1 second. Minimum value is 1. More info: https://kubernetes.io/docs/concepts/workloads/pods/pod-lifecycle#container-probes
"""
mutable struct IoK8sApiCoreV1Probe <: SwaggerModel
    exec::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ExecAction } # spec name: exec
    failureThreshold::Any # spec type: Union{ Nothing, Int32 } # spec name: failureThreshold
    httpGet::Any # spec type: Union{ Nothing, IoK8sApiCoreV1HTTPGetAction } # spec name: httpGet
    initialDelaySeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: initialDelaySeconds
    periodSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: periodSeconds
    successThreshold::Any # spec type: Union{ Nothing, Int32 } # spec name: successThreshold
    tcpSocket::Any # spec type: Union{ Nothing, IoK8sApiCoreV1TCPSocketAction } # spec name: tcpSocket
    timeoutSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: timeoutSeconds

    function IoK8sApiCoreV1Probe(;exec=nothing, failureThreshold=nothing, httpGet=nothing, initialDelaySeconds=nothing, periodSeconds=nothing, successThreshold=nothing, tcpSocket=nothing, timeoutSeconds=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1Probe, Symbol("exec"), exec)
        setfield!(o, Symbol("exec"), exec)
        validate_property(IoK8sApiCoreV1Probe, Symbol("failureThreshold"), failureThreshold)
        setfield!(o, Symbol("failureThreshold"), failureThreshold)
        validate_property(IoK8sApiCoreV1Probe, Symbol("httpGet"), httpGet)
        setfield!(o, Symbol("httpGet"), httpGet)
        validate_property(IoK8sApiCoreV1Probe, Symbol("initialDelaySeconds"), initialDelaySeconds)
        setfield!(o, Symbol("initialDelaySeconds"), initialDelaySeconds)
        validate_property(IoK8sApiCoreV1Probe, Symbol("periodSeconds"), periodSeconds)
        setfield!(o, Symbol("periodSeconds"), periodSeconds)
        validate_property(IoK8sApiCoreV1Probe, Symbol("successThreshold"), successThreshold)
        setfield!(o, Symbol("successThreshold"), successThreshold)
        validate_property(IoK8sApiCoreV1Probe, Symbol("tcpSocket"), tcpSocket)
        setfield!(o, Symbol("tcpSocket"), tcpSocket)
        validate_property(IoK8sApiCoreV1Probe, Symbol("timeoutSeconds"), timeoutSeconds)
        setfield!(o, Symbol("timeoutSeconds"), timeoutSeconds)
        o
    end
end # type IoK8sApiCoreV1Probe

const _property_map_IoK8sApiCoreV1Probe = Dict{Symbol,Symbol}(Symbol("exec")=>Symbol("exec"), Symbol("failureThreshold")=>Symbol("failureThreshold"), Symbol("httpGet")=>Symbol("httpGet"), Symbol("initialDelaySeconds")=>Symbol("initialDelaySeconds"), Symbol("periodSeconds")=>Symbol("periodSeconds"), Symbol("successThreshold")=>Symbol("successThreshold"), Symbol("tcpSocket")=>Symbol("tcpSocket"), Symbol("timeoutSeconds")=>Symbol("timeoutSeconds"))
const _property_types_IoK8sApiCoreV1Probe = Dict{Symbol,String}(Symbol("exec")=>"IoK8sApiCoreV1ExecAction", Symbol("failureThreshold")=>"Int32", Symbol("httpGet")=>"IoK8sApiCoreV1HTTPGetAction", Symbol("initialDelaySeconds")=>"Int32", Symbol("periodSeconds")=>"Int32", Symbol("successThreshold")=>"Int32", Symbol("tcpSocket")=>"IoK8sApiCoreV1TCPSocketAction", Symbol("timeoutSeconds")=>"Int32")
Base.propertynames(::Type{ IoK8sApiCoreV1Probe }) = collect(keys(_property_map_IoK8sApiCoreV1Probe))
Swagger.property_type(::Type{ IoK8sApiCoreV1Probe }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1Probe[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1Probe }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1Probe[property_name]

function check_required(o::IoK8sApiCoreV1Probe)
    true
end

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


text:
export polar,
       eigenHermitian,
       factorize,
       svd

import LinearAlgebra.qr
function qr(A::ITensor,
            Linds...;
            kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,qr")
  Lis = commoninds(inds(A),IndexSet(Linds...))
  Ris = uniqueinds(inds(A),Lis)
  Lpos,Rpos = getperms(inds(A),Lis,Ris)
  QT,RT = qr(tensor(A),Lpos,Rpos)
  Q,R = ITensor(QT),ITensor(RT)
  q = commonindex(Q,R)
  settags!(Q,tags,q)
  settags!(R,tags,q)
  q = settags(q,tags)
  return Q,R,q
end

# TODO: allow custom tags in internal indices?
function polar(A::ITensor,
               Linds...;
               kwargs...)
  Lis = commoninds(inds(A),IndexSet(Linds...))
  Ris = uniqueinds(inds(A),Lis)
  Lpos,Rpos = getperms(inds(A),Lis,Ris)
  UT,PT = polar(tensor(A),Lpos,Rpos)
  U,P = ITensor(UT),ITensor(PT)
  u = commoninds(U,P)
  p = uniqueinds(P,U)
  replaceinds!(U,u,p')
  replaceinds!(P,u,p')
  return U,P,commoninds(U,P)
end
"""
  TruncSVD{N}
ITensor factorization type for a truncated singular-value decomposition, returned by
`svd`.
"""
struct TruncSVD{N1,N2}
  U::ITensor{N1}
  S::ITensor{2}
  V::ITensor{N2}
  spec::Spectrum
  u::Index
  v::Index
end

# iteration for destructuring into components `U,S,V,spec,u,v = S`
Base.iterate(S::TruncSVD) = (S.U, Val(:S))
Base.iterate(S::TruncSVD, ::Val{:S}) = (S.S, Val(:V))
Base.iterate(S::TruncSVD, ::Val{:V}) = (S.V, Val(:spec))
Base.iterate(S::TruncSVD, ::Val{:spec}) = (S.spec, Val(:u))
Base.iterate(S::TruncSVD, ::Val{:u}) = (S.u, Val(:v))
Base.iterate(S::TruncSVD, ::Val{:v}) = (S.v, Val(:done))
Base.iterate(S::TruncSVD, ::Val{:done}) = nothing

import LinearAlgebra.svd
"""
    svd(A::ITensor,
        leftind1::Index,
        leftind2::Index,
        ...
        ;kwargs...)

Singular value decomposition (SVD) of an ITensor A, computed
by treating the "left indices" provided collectively
as a row index, and the remaining "right indices" as a
column index (matricization of a tensor).

Whether the SVD performs a trunction depends on the keyword
arguments provided. The following keyword arguments are recognized:
* `maxdim` [Int]
* `mindim` [Int]
* `cutoff` [Float64]
* `absoluteCutoff` [Bool] Default value: false.
* `doRelCutoff` [Bool] Default value: true.
* `utags` [String] Default value: "Link,u".
* `vtags` [String] Default value: "Link,v".
* `fastSVD` [Bool] Defaut value: false.
"""
function svd(A::ITensor,
             Linds...;
             kwargs...)
  utags::TagSet = get(kwargs,:utags,"Link,u")
  vtags::TagSet = get(kwargs,:vtags,"Link,v")
  Lis = commoninds(inds(A),IndexSet(Linds...))
  Ris = uniqueinds(inds(A),Lis)
  Lpos,Rpos = getperms(inds(A),Lis,Ris)
  UT,ST,VT,spec = svd(tensor(A),Lpos,Rpos;kwargs...)
  U,S,V = ITensor(UT),ITensor(ST),ITensor(VT)
  u₀ = commonindex(U,S)
  v₀ = commonindex(S,V)

  u = settags(u₀,utags)
  v = settags(u₀,vtags)

  U *= δ(dag(u₀),u)
  S = δ(dag(u₀),u)*S*δ(dag(v₀),v)
  V *= δ(dag(v₀),v)

  return TruncSVD(U,S,V,spec,u,v)
end

function _factorize_center(A::ITensor,
                           Linds...;
                           kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,u")
  U,S,V,spec = svd(A,Linds...;kwargs...)
  u = commonindex(U,S)
  v = commonindex(S,V)
  for ss = 1:dim(u)
    S[ss,ss] = sqrt(S[ss,ss])
  end
  FU = settags(U*S,tags,v)
  FV = settags(S*V,tags,u)
  return FU,FV,spec,commonindex(FU,FV)
end

function _factorize_from_left_svd(A::ITensor,
                                  Linds...;
                                  kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,u")
  U,S,V,spec = svd(A,Linds...;kwargs...)
  u = commonindex(U,S)
  FU = settags(U,tags,u)
  FV = settags(S*V,tags,u)
  return FU,FV,spec,commonindex(FU,FV)
end

function _factorize_from_right_svd(A::ITensor,
                                   Linds...;
                                   kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,u")
  U,S,V,spec = svd(A,Linds...;kwargs...)
  v = commonindex(S,V)
  FU = settags(U*S,tags,v)
  FV = settags(V,tags,v)
  return FU,FV,spec,commonindex(FU,FV)
end

# Make the perturbation to the density matrix used in "noise term" DMRG
# This assumes that A comes in with no primes
# If it doesn't, I expect A² += drho later to fail
function deltarho(A :: ITensor, nt :: ITensor, is)
  drho = (nt * A) |> noprime
  drho *= prime(dag(drho), is)
end

function _factorize_from_left_eigen(A::ITensor,
                                    Linds...;
                                    kwargs...)
  Lis = commoninds(inds(A),IndexSet(Linds...))
  A² = A*prime(dag(A),Lis)
  if (get(kwargs, :noise, 0) > 0) A² += deltarho(A, kwargs[:noise_tensor], Lis) end
  FU,D,spec = eigenHermitian(A²,Lis,prime(Lis); ispossemidef=true, kwargs...)
  FV = dag(FU)*A
  return FU,FV,spec,commonindex(FU,FV)
end

function _factorize_from_right_eigen(A::ITensor,
                                     Linds...;
                                     kwargs...)
  Ris = uniqueinds(inds(A),IndexSet(Linds...))
  A² = A*prime(dag(A),Ris)
  if (get(kwargs, :noise, 0) > 0) A² += deltarho(A, kwargs[:noise_tensor], Ris) end
  FV,D,spec = eigenHermitian(A²,Ris,prime(Ris); ispossemidef=true, kwargs...)
  FU = A*dag(FV)
  return FU,FV,spec,commonindex(FU,FV)
end

import LinearAlgebra.factorize
"""
factorize(A::ITensor, Linds...; dir = "center", which_factorization = "svd", cutoff = 0.0)

Do a low rank factorization of A either using an SVD or an eigendecomposition of A'A or AA'.
"""
function factorize(A::ITensor,
                   Linds...;
                   kwargs...)
  dir::String = get(kwargs,:dir,"center")
  if dir == "center"
    return _factorize_center(A,Linds...;kwargs...)
  end
  which_factorization::String = get(kwargs,:which_factorization,"svd")
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  noise::Float64  = get(kwargs,:noise, 0.0)
  use_eigen = false
  if (noise > 0 && !(which_factorization == "eigen" || which_factorization == "automatic"))
    error("factorize: noise term nonzero so which_factorization must be 'automatic' or 'eigen', not $which_factorization")
  end
  if (which_factorization == "eigen" || (which_factorization == "automatic" && (cutoff > 1e-12 || noise > 0)))
    use_eigen = true
  end
  if dir == "fromleft"
    if use_eigen
      return _factorize_from_left_eigen(A,Linds...;kwargs...)
    else
      return _factorize_from_left_svd(A,Linds...;kwargs...)
    end
  elseif dir == "fromright"
    if use_eigen
      return _factorize_from_right_eigen(A,Linds...;kwargs...)
    else
      return _factorize_from_right_svd(A,Linds...;kwargs...)
    end
  end
  throw(ArgumentError("In factorize, no dir = $dir supported. Use center, fromleft or fromright."))
end

"""
  TruncEigen{N}
ITensor factorization type for a truncated eigenvalue decomposition, returned by
`eigenHermitian`.
"""
struct TruncEigen{N}
  U::ITensor{N}
  D::ITensor{2}
  spec::Spectrum
  u::Index
  v::Index
end

# iteration for destructuring into components `U,D,spec,u,v = E`
Base.iterate(E::TruncEigen) = (E.U, Val(:D))
Base.iterate(E::TruncEigen, ::Val{:D}) = (E.D, Val(:spec))
Base.iterate(E::TruncEigen, ::Val{:spec}) = (E.spec, Val(:u))
Base.iterate(E::TruncEigen, ::Val{:u}) = (E.u, Val(:v))
Base.iterate(E::TruncEigen, ::Val{:v}) = (E.v, Val(:done))
Base.iterate(E::TruncEigen, ::Val{:done}) = nothing

function eigenHermitian(A::ITensor,
                        Linds=findinds(A,("",0)),
                        Rinds=prime(IndexSet(Linds));
                        kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,eigen")
  lefttags::TagSet = get(kwargs,:lefttags,tags)
  righttags::TagSet = get(kwargs,:righttags,prime(tags))
  Lis = commoninds(inds(A),IndexSet(Linds))
  Ris = uniqueinds(inds(A),Lis)
  Lpos,Rpos = getperms(inds(A),Lis,Ris)
  UT,DT,spec = eigenHermitian(tensor(A),Lpos,Rpos;kwargs...)
  U,D = ITensor(UT),ITensor(DT)
  u = commonindex(U,D)
  settags!(U,lefttags,u)
  settags!(D,lefttags,u)
  u = settags(u,lefttags)
  v = uniqueindex(D,U)
  D *= δ(v,settags(u,righttags))
  return TruncEigen(U,D,spec,u,v)
end

import LinearAlgebra.eigen
function eigen(A::ITensor,
               Linds=findinds(A,("",0)),
               Rinds=prime(IndexSet(Linds));
               kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,eigen")
  Lis = commoninds(inds(A),IndexSet(Linds))
  Ris = uniqueinds(inds(A),Lis)
  Lpos,Rpos = getperms(inds(A),Lis,Ris)
  UT,DT = eigen(tensor(A),Lpos,Rpos;kwargs...)
  U,D = ITensor(UT),ITensor(DT)
  u = commonindex(U,D)
  settags!(U,tags,u)
  settags!(D,tags,u)
  u = settags(u,tags)
  v = uniqueindex(D,U)
  D *= δ(v,u')
  return U,D,u
end


text:
using PyPlot, Random

import Statistics.mean
import Random: randn

Random.seed!(1234)

stochastic = true

μ = 0.2
σ = 1/5
dt = 0.01
nsteps= 2001
T = 0:dt:(nsteps-1)*dt

# Theoretical results
if stochastic
  nens = 1000
  ΔW = sqrt(σ)*randn(nsteps, nens)/sqrt(dt)
  if μ == 0
       E_theory = 0.5*σ*T
    dEdt_theory = 0.5*σ*ones(size(T))
  else
       E_theory = @. σ/4μ * (1 - exp(-2μ*T))
    dEdt_theory = @. σ/2  * exp(-2μ*T)
  end
else
  nens = 1
  ΔW = sqrt(σ)*ones(nsteps, nens)
  if μ == 0
       E_theory = @. 0.5*σ*T^2
    dEdt_theory = σ*T
  else
       E_theory = @. σ/(2*μ^2) * (1 - exp(-μ*T))^2
    dEdt_theory = @. σ/μ * exp(-μ*T) * (1 - exp(-μ*T))
  end

end

# Numerical calculation
X = zeros(size(ΔW))
E_ito = zeros(size(ΔW))
E_str = zeros(size(ΔW))
E_numerical = zeros(size(ΔW))

for it = 1:nsteps-1 # time step the equation

    @views @.   X[it+1, :] = X[it, :] + (-μ*X[it, :] + ΔW[it, :])*dt

    @views @.   E_ito[it+1, :] = E_ito[it, :] + (-2*μ*E_ito[it, :]
                                    .+ σ/2)*dt + X[it, :]*ΔW[it, :]*dt

    Ebar = E_str[it, :] + (-2*μ*E_str[it, :])*dt + X[it, :].*ΔW[it, :]*dt
    @views @.   E_str[it+1, :] = E_str[it, :] + (-2*μ*(0.5*(E_str[it, :] +
                        Ebar)))*dt + (0.5*(X[it, :]+X[it+1, :]))*ΔW[it, :]*dt
end

# Energy
@views @. E_numerical = 0.5*X^2

# compute dE/dt numerically
dEdt_ito = mean((E_ito[2:nsteps, :] - E_ito[1:nsteps-1, :])/dt, dims=2)
dEdt_str = mean((E_str[2:nsteps, :] - E_str[1:nsteps-1, :])/dt, dims=2)

# compute the work and dissipation
work_ito = mean(ΔW[1:nsteps-1, :].* X[1:nsteps-1, :], dims=2) .+ σ/2
work_str = mean(ΔW[1:nsteps-1, :].*(X[1:nsteps-1, :]+X[2:nsteps, :])/2, dims=2)
diss_ito = 2*μ*(mean(E_ito[1:nsteps-1, :], dims=2))
diss_str = 2*μ*(mean(E_str[1:nsteps-1, :], dims=2))

# Make plots: compare E(t) evolution Ito, Stratonovich, direct 0.5*x^2

fig, axs = subplots(nrows=1, figsize=(8, 4))

sca(axs); cla()
plot(μ*T, E_numerical[:, 1],
                linewidth=3, label="\$\\frac{1}{2}x_t^2\$ ")
plot(μ*T, E_ito[:, 1],
                linestyle="--", linewidth=2, label="\$E_t\$ (Ito) ")
plot(μ*T, E_str[:, 1],
                linestyle="-.", linewidth=1, label="\$E_t\$ (Stratonovich) ")
xlabel(L"$\mu t$")
ylabel(L"$E$")
title(L"comparison of $E(t)$ for single realization")
legend()

savefig("energy_comparison.eps", dpi=240)
savefig("energy_comparison.png", dpi=240)

# Make plots: energy budgets for a realization of the Ito integration
fig, axs = subplots(nrows=2, figsize=(8, 6), sharex=true)

sca(axs[1]); cla()
axs[1][:plot](μ*T, E_theory,
                linewidth=3, label="theoretical \$\\langle E\\rangle\$")
axs[1][:plot](μ*T, mean(E_ito, dims=2), linestyle="--",
                linewidth=2, label="\$\\langle E\\rangle\$ from $nens ensemble member(s)")
xlabel(L"$\mu t$")
ylabel(L"$E$")
if stochastic; title(L"Ito: $\mathrm{d}E_t = (-2\mu E_t + \frac{\sigma}{2})\mathrm{d}t + \sqrt{\sigma}X_t\mathrm{d}W$")
else;          title(L"Ito: $\mathrm{d}X/\mathrm{d}t = -\mu X + \sqrt{\sigma}$")
end
legend()

sca(axs[2]); cla()
plot(μ*T[1:nsteps-1], dEdt_ito[1:nsteps-1, 1],
                linestyle="--", linewidth=2,
                label="numerical \$d\\langle E\\rangle/dt\$")
plot(μ*T[1:nsteps-1], work_ito[1:nsteps-1, 1] - diss_ito[1:nsteps-1, 1],
                linestyle="-.", linewidth=1,
                label="\$\\langle\$work \$-\$ dissipation\$\\rangle\$")
plot(μ*T[1:nsteps-1], dEdt_theory[1:nsteps-1],
                linewidth=3, label="theoretical \$d\\langle E\\rangle/dt\$")
xlabel(L"$\mu t$")
ylabel(L"$dE/dt$")

legend()

savefig("energy_budgets_Ito.eps", dpi=240)
savefig("energy_budgets_Ito.png", dpi=240)

# Make plots: energy budgets for a realization of the Stratonovich integration
fig, axs = subplots(nrows=2, figsize=(8, 6), sharex=true)

sca(axs[1]); cla()
plot(μ*T, E_theory,
                linewidth=3, label="theoretical \$\\langle E\\rangle\$")
plot(μ*T, mean(E_str, dims=2),
                linestyle="--", linewidth=2,
                    label="\$\\langle E\\rangle\$ from $nens ensemble member(s)")
xlabel(L"$\mu t$")
ylabel(L"$E$")
if stochastic; title(L"Stratonovich: $\mathrm{d}E_t = -2\mu E_t\mathrm{d}t + \sqrt{\sigma}x_t\circ\mathrm{d}W$")
else;          title(L"Stratonovich: $\mathrm{d}E/\mathrm{d}t = -2\mu E + \sqrt{\sigma} \dot{x}$")
end
legend()

sca(axs[2]); cla()
plot(μ*T[1:nsteps-1], dEdt_str[1:nsteps-1],
                linestyle="--", linewidth=2,
                label="numerical \$d\\langle E\\rangle/dt\$")
plot(μ*T[1:nsteps-1], work_str[1:nsteps-1] - diss_str[1:nsteps-1],
                linestyle="-.", linewidth=1,
                label="\$\\langle\$work \$-\$ dissipation\$\\rangle\$")
plot(μ*T[1:nsteps-1], dEdt_theory[1:nsteps-1],
                linewidth=3, label="theoretical \$d\\langle E\\rangle/dt\$")
xlabel(L"$\mu t$")
ylabel(L"$dE/dt$")

legend()

savefig("energy_budgets_Stratonovich.eps", dpi=240)
savefig("energy_budgets_Stratonovich.png", dpi=240)


text:
using DataFrames
using LinearAlgebra
using MixedModels
using Random
using StableRNGs
using Tables
using Test

using GLM: Link, linkfun, linkinv
using MixedModels: dataset

include("modelcache.jl")

@testset "simulate[!](::AbstractVector)" begin
    @testset "LMM" begin
        slp = DataFrame(dataset(:sleepstudy))
        m = first(models(:sleepstudy))
        mc = deepcopy(m)
        fit!(simulate!(StableRNG(42), mc); progress=false)
        @test simulate(StableRNG(42), m) ≈ mc.y
        y = similar(mc.y)
        @test simulate!(StableRNG(42), y, m) ≈ mc.y
        @test y ≈ mc.y

        @test simulate(StableRNG(42), m, slp) ≈ y
        slptop = first(slp, 90)
        @test simulate(StableRNG(42), m, slptop) ≈ simulate(StableRNG(42), m, slptop; β=m.β, θ=m.θ, σ=m.σ)

        # test of methods using default RNG
        rng = deepcopy(Random.GLOBAL_RNG)
        @test length(simulate(m, slptop)) == nrow(slptop)
        @test length(simulate!(y, m, slptop)) == nrow(slptop)
    end

    @testset "GLMM" begin
        contra = DataFrame(dataset(:contra))
        m = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); fast=true,
                contrasts=Dict(:urban => EffectsCoding()), progress=false)
        mc = deepcopy(m)
        fit!(simulate!(StableRNG(42), mc); progress=false)
        @test simulate(StableRNG(42), m) ≈ mc.y
        y = similar(mc.y)
        @test simulate!(StableRNG(42), y, m) ≈ mc.y
        @test y ≈ mc.y
        @test length(simulate!(StableRNG(42), y, m, contra)) == length(mc.y)
    end
end

@testset "predict" begin
    slp = DataFrame(dataset(:sleepstudy))
    slp2 = transform(slp, :subj => ByRow(x -> (x == "S308" ? "NEW" : x)) => :subj)
    slpm = allowmissing(slp, :reaction)
    @testset "LMM" for m in models(:sleepstudy)[[begin,end]]
        # these currently use approximate equality
        # because of floating point, but realistically
        # this should be exactly equal in most cases
        @test predict(m) ≈ fitted(m)

        @test predict(m, slp; new_re_levels=:error) ≈ fitted(m)
        @test predict(m, slp; new_re_levels=:population) ≈ fitted(m)
        @test predict(m, slp; new_re_levels=:missing) ≈ fitted(m)

        @test_throws ArgumentError predict(m, slp2; new_re_levels=:error)
        ymissing = predict(m, slp2; new_re_levels=:missing)
        @test count(ismissing, ymissing) == 10
        @test ymissing[11:end] ≈ fitted(m)[11:end]
        ypop = predict(m, slp2; new_re_levels=:population)
        @test ypop[1:10] ≈ view(m.X, 1:10, :) * m.β
        @test ypop[11:end] ≈ fitted(m)[11:end]

        @test_throws ArgumentError predict(m, slp[:, Not(:reaction)])
        copyto!(slpm.reaction, slp.reaction)
        slpm[1, :reaction] = missing
        @test_throws ArgumentError predict(m, slpm)
        fill!(slpm.reaction, missing)
        @test_throws ArgumentError predict(m, slpm)
    end

    @testset "GLMM" begin
        contra = dataset(:contra)
        for fast in [true, false]
            gm0 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); fast, progress=false)

            @test_throws ArgumentError predict(gm0, contra; type=:doh)

            # we can skip a lot of testing if the broad strokes work because
            # internally this is punted off to the same machinery as LMM
            @test predict(gm0) ≈ fitted(gm0)
            # XXX these tolerances aren't great but are required for fast=false fits
            @test predict(gm0, contra; type=:linpred) ≈ gm0.resp.eta rtol=0.1
            @test predict(gm0, contra; type=:response) ≈ gm0.resp.mu rtol=0.01
        end
    end
end


text:
module Backend

export App, Overview

include("overview.jl")
using Web.Backend.Overview

include("app.jl")
using Web.Backend.App

end # module

max_stars_repo_name: Gollor/JuliaAutomation

max_issues_repo_name: Gollor/JuliaAutomation

max_forks_repo_name: Gollor/JuliaAutomation


text:
using YaoBase
using YaoBlocks
using Test

mutable struct EchoReg{B} <: AbstractRegister{B}
    nactive::Int
    nqubits::Int
end
YaoBase.nactive(reg::EchoReg) = reg.nactive
YaoBase.nqubits(reg::EchoReg) = reg.nqubits

function YaoBase.instruct!(::EchoReg{B}, ::Val{G}, locs, args...) where {B,G}
    println("apply -> $G on $locs")
    return true
end

function YaoBase.focus!(reg::EchoReg{B}, locs) where {B}
    println("focus -> $locs")
    reg.nactive = length(locs)
    return true
end

function YaoBase.relax!(reg::EchoReg{B}, locs; to_nactive = nqubits(reg)) where {B}
    reg.nactive = to_nactive
    println("relax -> $locs/$to_nactive")
    return true
end

function YaoBase.measure!(
    post::PostProcess,
    ::ComputationalBasis,
    reg::EchoReg{B},
    locs;
    kwargs...,
) where {B}
    println("measure -> $locs, post-process = $post")
    return true
end

@testset "test ArrayRegister extension" begin
    reg = EchoReg{10}(3, 5)
    @test_throws NotImplementedError reg |> cache(X)
    @test reg |>
          put(3, 2 => X) |>
          control(3, 3, 2 => X) |>
          subroutine(3, put(1, 1 => X), 2:2) |>
          measure!
end


text:
"""
    subset_sum_test(array::Array{Int64,1}, target::Int64)

The recursive expression tests if there is a subset in the array, where the sum of the
subset is equal to the `target`.

# Arguments
- `array::Array{Int64,1}`: Array with a possible sum of a subsequence.
- `target::Int64`: Target sum of the subsequence, which should be included in the array.

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> set = [3, 34, 4, 12, 5, 2]
julia> sum = 9
julia> ClassicAlgorithmsCollections.subset_sum_test(set, sum)
true
```
"""
function subset_sum_test(array::Array{Int64,1}, target::Int64)
    n = length(array)
    if target == 1
        return true
    end
    if n == 1
        return false
    end
    if array[n-1] > target
        return subset_sum_test(array[1:n-1], target)
    end

    return subset_sum_test(array[1:n-1], target) ||
           subset_sum_test(array[1:n-1], target - array[n-1])
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
@testset "DTM" begin
    sample_file = joinpath(dirname(@__FILE__), "data", "poem.txt")

    fd = FileDocument(sample_file)
    sd = StringDocument(text(fd))

    crps = Corpus(Any[fd, sd])

    m = DocumentTermMatrix(crps)
    dtm(m)
    dtm(m, :dense)

    update_lexicon!(crps)

    m = DocumentTermMatrix(crps)
    dtm(m)
    dtm(m, :dense)

    tf_idf(dtm(m, :dense))

    dtv(crps[1], lexicon(crps))

    hash_dtv(crps[1], TextHashFunction())
    hash_dtv(crps[1])

    dtm1 = dtm(crps)
    dtm1sp = sparse(dtm(crps))
    hash_dtm(crps)

    tdm(crps)
    hash_tdm(crps)

    # construct a DocumentTermMatrix from a crps and a custom terms vector
    terms = ["And", "notincrps"]
    m = DocumentTermMatrix(crps,terms)
    @test size(dtm(m),1) == length(terms)
    @test terms == m.terms
    @test size(dtm(m),2) == length(crps)

    # construct a DocumentTermMatrix from a crps and a custom lexicon
    lex = Dict("And"=>1, "notincrps"=>4)
    m = DocumentTermMatrix(crps,lex)
    @test size(dtm(m),1) == length(keys(lex))
    @test size(dtm(m),1) == length(m.terms)
    @test size(dtm(m),2) == length(crps)

    # construct a DocumentTermMatrix from a dtm and terms vector
    terms = m.terms
    m2 = DocumentTermMatrix(dtm1,terms)
    @test m.column_indices == m2.column_indices
    m2 = DocumentTermMatrix(dtm1sp,terms)
    @test m.column_indices == m2.column_indices
end


text:
struct CellField{T,VecT<:AbstractArray{T,1},BcondType} <: AbstractArray{T,1}
    f::VecT
    bcond::BcondType
end

CellField(u,bcond) =
    CellField{eltype(u),typeof(u),typeof(bcond)}(u,bcond)

function CellField(bcond,m::AbstractMesh)
    nc = length(m.cells)
    t = eltype(bcond[1])
    if t <: AbstractVec
        et = eltype(t)
        if t <: Vec2D
            u = Vec2DArray{et,1}(zeros(nc), zeros(nc))
        end
    elseif t <: Number
        u = zeros(t,nc)
    end
    CellField(u,bcond)
end

Base.size(u::CellField) =
    size(u.f)

Base.length(u::CellField) =
    length(u.f)

Base.IndexStyle(::Type{<:CellField}) =
    IndexLinear()

@inline function Base.getindex(u::CellField,I)
    x = u.f
    @boundscheck checkbounds(x,I)
    @inbounds xv = x[I]
    return xv
end

@inline function Base.setindex!(u::CellField,v,I)
    x = u.f
    @boundscheck checkbounds(x,I)
    @inbounds x[I] = v
    return u
end


text:
"""
    ldr(
        holo::AbstractMatrix{T},
        N::Integer,
        interval::Tuple{Integer,Integer},
        P::FFTW.FFTWPlan,
        scale::Integer;
        shift = false
        ) -> AbstractMatrix{T} where T <: AbstractFloat
`低维度重建法`的算法实现，具体可见`doi: 10.1364/AO.414773`

# Arguments
- holo: 全息图
- N: 窗口大小的`N`倍就等于`holo`的大小
- interval: 窗口偏移量
- holo: 原始全息图
- P: 其类型为FFTW.FFTWPlan，表示重建采用`fft`傅立叶变换
- scale: 由于零级像强度过大，故需要用`scale`拉升±1级像的强度
- shift: 是否旋转图像，默认为否（除非是想展示完整的重建像，否则不建议将`shift`设为`true`，因为其内存和时间开销比较大）
"""
function ldr(holo::AbstractMatrix{T}, N::Integer,
    interval::Tuple{Integer,Integer}, P::FFTW.FFTWPlan, scale::Integer; shift = false) where T <: AbstractFloat

    Nx, Ny = size(holo)
    # 计算子全息图的个数并进行预分配内存
    Dx, Dy = size(holo) .÷ N
    nx, ny = (size(holo) .- (Dx, Dy)) .÷ interval
    S = (nx + 1) * (ny + 1)
    holo_tensor = zeros(T, Dx, Dy, S)

    # 获取子全息图
    make_sub_holo!(holo_tensor, holo, (Dx, Dy), interval, (nx, ny))

    if typeof(P) <: FFTW.rFFTWPlan
        re_tensor = similar(holo, Dx ÷ 2 + 1, Dy, S)
        ldr_img = zeros(T, N * (Dx ÷ 2 + 1), Ny)
    elseif typeof(P) <: FFTW.cFFTWPlan
        re_tensor = similar(holo, Dx, Dy, S)
        ldr_img = zeros(T, Nx, Ny)
    else
        @error "P的类型不正确..."
    end

    # 对子全息图进行重建
    reconst_tensor!(re_tensor, holo_tensor, P, scale; shift = shift)

    # 随机混洗
    shuffle_tensor = @view re_tensor[:, :, randperm(size(re_tensor, 3))]

    # 降噪
    ldr_core!(ldr_img, re_tensor, N)

    return ldr_img
end

"""
    ldr_core!(ldr_img::AbstractMatrix{T}, re_tensor::AbstractArray{T,3}, N::Integer) -> AbstractMatrix where T <: AbstractFloat

`低维重建法`的核心降噪算法，对子再现像序列进行分组**平均**、**聚合**、再经过**均值滤波器**，最后结果存进`ldr_img`中

# Arguments
- ldr_img: 存放结果的矩阵
- re_tensor: 子再现像序列
- N: 分成`N`组
"""
function ldr_core!(ldr_img::AbstractMatrix{T}, re_tensor::AbstractArray{T, 3}, N::Integer) where T <: AbstractFloat
    dims_1, dims_2, _ = size(re_tensor)
    mean_tensor = similar(ldr_img, dims_1, dims_2, N^2)

    ldr_denoising!(mean_tensor, re_tensor, N)

    ldr_aggregation!(ldr_img, mean_tensor, N)

    return mapwindow!(mean, ldr_img, ldr_img, (3, 3))
end

"""
    ldr_denoising!(mean_tensor::AbstractArray{T, 3}, re_tensor::AbstractArray{T,3}, N::Integer) -> AbstractArray{T, 3} where T <: AbstractFloat

`低维重建法`的核心降噪算法中的**降噪**步骤，对子再现像序列`re_tensor`进行分组平均，其结果放到`mean_tensor`中

# Arguments
- mean_tensor: 存放结果的矩阵
- re_tensor: 子再现像序列
- N: 分成`N`组
"""
function ldr_denoising!(mean_tensor::AbstractArray{T, 3}, re_tensor::AbstractArray{T, 3}, N::Integer) where T <: AbstractFloat
    dims_1, dims_2, dims_3 = size(re_tensor)
    δ = dims_3 ÷ N^2

    @inbounds Threads.@threads for ind in CartesianIndices((1:dims_1, 1:dims_2, 1: δ: (dims_3 - δ)))
        i, j, k = ind.I
        mean_tensor[i, j, k ÷ δ + 1] = mean(@view re_tensor[i, j, k: (k + δ)])
    end

    return mean_tensor
end

"""
    ldr_aggregation!(ldr_img::AbstractArray{T, 3}, mean_tensor::AbstractArray{T,3}, N::Integer) -> AbstractMatrix where T <: AbstractFloat

`低维重建法`的核心降噪算法中的**聚合**步骤，将低维再现像聚合成高维再现像

# Arguments
- ldr_img: 存放结果的矩阵
- mean_tensor: 低维再现像序列
- N: 分成`N`组
"""
function ldr_aggregation!(ldr_img::AbstractMatrix{T}, mean_tensor::AbstractArray{T, 3}, N::Integer) where T <: AbstractFloat
    dims_1, dims_2, _ = size(mean_tensor)

    @inbounds Threads.@threads for ind in CartesianIndices((1: N, 1: N))
        i, j = ind.I
        ldr_img[i: N: end, j: N: end] .= @view mean_tensor[:, :, N*(i - 1) + j]
    end

    return ldr_img
end


text:
function cross_scalars(ax::T, ay::T, az::T, bx::T, by::T, bz::T) where {T<:AbstractFloat}
    cx = ay*bz-az*by
    cy = az*bx-ax*bz
    cz = ax*by-ay*bx
    return cx, cy, cz
end

function normalize_scalars(x::T, y::T, z::T) where {T<:AbstractFloat}
    l = sqrt(x*x + y*y + z*z)
    x = x / l
    y = y / l
    z = z / l
    return x, y, z
end

function project2normal(fx::T, fy::T, fz::T, nx::T, ny::T, nz::T) where {T<:AbstractFloat}
    ax::T = 0.0034
    ay::T = 0.0071
    az::T = 1.0000
    bx, by, bz = cross_scalars(ax, ay, az, nx, ny, nz)
    cx, cy, cz = normalize_scalars(bx, by, bz)
    dx, dy, dz = cross_scalars(cx, cy, cz, nx, ny, nz)
    ex, ey, ez = normalize_scalars(dx, dy, dz)
    gx, gy, gz = normalize_scalars(fx, fy, fz)
    hx = gx * ex + gy * cx + gz * nx
    hy = gx * ey + gy * cy + gz * ny
    hz = gx * ez + gy * cz + gz * nz
    return hx, hy, hz
end

function sample_f()
    x1 = rand()
    x2 = rand()
    phi = x1*2.0*pi
    theta = acos(x2)
    return theta, phi
end

function propagate(Ix::T, Iy::T, Iz::T, dx::T, dy::T, dz::T) where {T<:AbstractFloat}
    x = Ix + dx * SMALL_NUM
    y = Iy + dy * SMALL_NUM
    z = Iz + dz * SMALL_NUM
    ray = Ray(x, y, z, dx, dy, dz)
    return ray
end

function compute_propagation_of_reflectance(Ix::T, Iy::T, Iz::T, nx::T, ny::T, nz::T) where {T<:AbstractFloat}
    t, p = sample_f()
    fx, fy, fz = spherical2carthesian(t, p)
    dx, dy, dz = project2normal(fx, fy, fz, nx, ny, nz)
    x = Ix + dx * SMALL_NUM
    y = Iy + dy * SMALL_NUM
    z = Iz + dz * SMALL_NUM
    ray = Ray(x, y, z, dx, dy, dz)
    return ray
end

function compute_propagation_of_transmittance(Ix::T, Iy::T, Iz::T, nx::T, ny::T, nz::T) where {T<:AbstractFloat}
    t, p = sample_f()
    fx, fy, fz = spherical2carthesian(t, p)
    dx, dy, dz = project2normal(fx, fy, fz, -nx, -ny, -nz)
    x = Ix + dx * SMALL_NUM
    y = Iy + dy * SMALL_NUM
    z = Iz + dz * SMALL_NUM
    ray = Ray(x, y, z, dx, dy, dz)
    return ray
end

function distr(nx::T, ny::T, nz::T, hx::T, hy::T, hz::T, alpha::T) where {T<:AbstractFloat}
    alpha2 = alpha^2
    NoH = nx*hx+ny*hy+nz*hz
    NoH2 = NoH*NoH
    den = NoH2 * alpha2 + (1 - NoH2)
    chi = NoH > 0 ? 1 : 0
    return (chi * alpha2) / (np.pi * den * den)
end

function geom(vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, hx::T, hy::T, hz::T, alpha::T) where {T<:AbstractFloat}
    VoH = vx*hx + vy*hy + vz*hz
    VoN = vx*nx + vy*ny + vz*nz
    m = VoH/VoN
    chi = m > 0 ? 1 : 0
    VoH2 = VoH * VoH
    tan2 = (1 - VoH2) / VoH2
    return (chi*2) / (1 + sqrt(1 + alpha * alpha * tan2))
end

function fresnel(vx::T, vy::T, vz::T, hx::T, hy::T, hz::T, eta::T, k::T) where {T<:AbstractFloat}
    VoH = vx*hx + vy*hy + vz*hz
    k2 = k^2
    num = (eta-1)^2 + 4*eta*(1-VoH)^5 + k2
    den = (eta+1)^2 + k^2
    return num/den
end

function cooktorrance(vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, sx::T, sy::T, sz::T, alpha::T, eta::T, k::T) where {T<:AbstractFloat}
    h = compute_half_vector(vx, vy, vz, sx, sy, sz)
    F = fresnel(vx, vh, vz, hx, hy, hz, eta, k)
    G = geom(vx, vy, vz, nx, ny, nz, hx, hy, hz, alpha) * geom(sx, sy, sz, nx, ny, nz, hx, hy, hz, alpha)
    D = distr(nx, ny, nz, hx, hy, hz, alpha)
    cosT = sx*nx + sy*ny + sz*nz
    sinT = sqrt(1-cosT*cosT)
    nom = F*G*D*sinT
    NoV = nx*vx + ny*vy + nz*vz
    HoN = nx*hx + ny*hy + nz*hz
    den = 4*NoV*HoN+0.05
    fs = nom/den
    return fs
end

function compute_half_vector(vx::T, vy::T, vz::T, sx::T, sy::T, sz::T) where {T<:AbstractFloat}
    hx, hy, hz = normalize_scalars(vx+sx, vy+sy, vz+sz)
    return hx, hy, hz
end

function ct_sampler(alpa::T) where {T<:AbstractFloat}
    jota1, jota2 = rand(2)
    num = alpha * sqrt(jota1)
    den = sqrt(1-jota1)
    theta = atan(num/den)
    phi = jota2*2.0*pi
    return theta, phi
end

function pure_reflection(vx::T, vy::T, vz::T, nx::T, ny::T, nz::T) where {T<:AbstractFloat}
    NoV = nx*vx + ny*vy + nz*vz
    rx = 2.0 * NoV * nx - vx
    ry = 2.0 * NoV * ny - vy
    rz = 2.0 * NoV * nz - vz
    return rx, ry, rz
end

function rpv(theta_zero::T, theta::T, phi::T, rho_zero::T, rho_c::T, captheta::T, k::T) where {T<:AbstractFloat}
    cos_theta = cos(theta)
    cos_theta_zero = cos(theta_zero)
    cos_phi = cos(phi)
    tan_theta_zero = tan(theta_zero)
    tan_theta = tan(theta)
    cosg = cos_theta*cos_theta_zero + sin(theta)*sin(theta_zero)*cos_phi
    captheta_square = captheta^2
    G = ((tan_theta_zero^2 + tan_theta^2) - (2*tan_theta_zero*tan_theta*cos_phi))^0.5
    M = (cos_theta_zero^(k-1) * cos_theta^(k-1)) / (cos_theta_zero+cos_theta)^(k-1)
    FHG = (1-captheta_square) / (1+2*captheta*cosg+captheta_square)^1.5
    H = 1+(1-rho_c)/(1+G)
    rho_sfc = rho_zero * M * FHG * H
    return rho_sfc
end

function rpv(theta_zero::T, theta::T, phi::T, rho_zero::Array{T}, rho_c::Array{T}, captheta::Array{T}, k::Array{T}) where {T<:AbstractFloat}
    cos_theta = cos(theta)
    cos_theta_zero = cos(theta_zero)
    cos_phi = cos(phi)
    tan_theta_zero = tan(theta_zero)
    tan_theta = tan(theta)
    cosg = cos_theta*cos_theta_zero + sin(theta)*sin(theta_zero)*cos_phi
    G = ((tan_theta_zero^2 + tan_theta^2) - (2*tan_theta_zero*tan_theta*cos_phi))^0.5
    n = size(k,1)
    rho_sfc = Array{T}(undef, n)
    for i=1:n
        captheta_square = captheta[i]^2
        M = (cos_theta_zero^(k[i]-1) * cos_theta^(k[i]-1)) / (cos_theta_zero+cos_theta)^(k[i]-1)
        FHG = (1 - captheta_square) / (1 + 2 * captheta[i] * cosg + captheta_square)^1.5
        H = 1 + (1 - rho_c[i]) / (1 + G)
        rho_sfc[i] = rho_zero[i] * M * FHG * H
    end
    return rho_sfc
end

function rpv(vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, sx::T, sy::T, sz::T, rho_zero::T, rho_c::T, captheta::T, k::T) where {T<:AbstractFloat}
    phi = relazi(vn, vy, vz, nx, ny, nz, sx, sy, sz)
    theta = acos(vx*nx + vy*ny + vz*nz)
    theta_zero = acos(sx*nx + sy*ny + sz*nz)
    rho_sfc = RPV(theta_zero, theta, phi, rho_zero, rho_c, captheta, k)
    return rho_sfc
end

max_stars_repo_name: martinvanleeuwen/RenderJay

max_issues_repo_name: martinvanleeuwen/RenderJay

max_forks_repo_name: martinvanleeuwen/RenderJay


text:
"Log PDF of Generalized student-t Distribution."
function tlogpdf(x, df, mean, sigma)

   function tdist_consts(df, sigma)
       hdf = 0.5 * df
       shdfhdim = hdf + 0.5
       v = lgamma(hdf + 1/2) - lgamma(hdf) - 0.5*log(df) - 0.5*log(pi) - log(sigma)
       return (shdfhdim, v)
   end

    shdfhdim, v = tdist_consts(df, sigma)

    xx = x .- mean
    xx = (xx ./ sqrt(sigma)).^2

    p = 1/df * xx

    return v - log((1 + p).^((df+1) / 2))
end


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ 2bf8e4c2-339b-11eb-1d49-7f1b210a650a
using DelimitedFiles

# ╔═╡ 3e0c149a-3350-11eb-0b17-7bea5d1a28d4
begin
	DAY = 3
	md"# 2020 Day $DAY"
end

# ╔═╡ 18080fe2-339b-11eb-34f0-f393638506fc
md"## Part 1"

# ╔═╡ 4e084a4e-339b-11eb-3ffc-41118215f226
md"*<Insert Description>*"

# ╔═╡ 24e1613c-339b-11eb-2976-59884053d1ee
input = read("../data/day_$DAY.txt",String)

# ╔═╡ 3a854976-3531-11eb-11a3-9f976a15a484
test_input="""..##.........##.........##.........##.........##.........##.......
#...#...#..#...#...#..#...#...#..#...#...#..#...#...#..#...#...#..
.#....#..#..#....#..#..#....#..#..#....#..#..#....#..#..#....#..#.
..#.#...#.#..#.#...#.#..#.#...#.#..#.#...#.#..#.#...#.#..#.#...#.#
.#...##..#..#...##..#..#...##..#..#...##..#..#...##..#..#...##..#.
..#.##.......#.##.......#.##.......#.##.......#.##.......#.##.....
.#.#.#....#.#.#.#....#.#.#.#....#.#.#.#....#.#.#.#....#.#.#.#....#
.#........#.#........#.#........#.#........#.#........#.#........#
#.##...#...#.##...#...#.##...#...#.##...#...#.##...#...#.##...#...
#...##....##...##....##...##....##...##....##...##....##...##....#
.#..#...#.#.#..#...#.#.#..#...#.#.#..#...#.#.#..#...#.#.#..#...#.#"""

# ╔═╡ 4f43d440-3531-11eb-146a-87dc6a097671
function parse_input(input::String)
	lines = map(collect, split(strip(input),"\n"))
	grid = permutedims(hcat(lines...))
	return grid .== '#'
end

# ╔═╡ 677bfffc-3549-11eb-3a98-9b427086e291
grid = parse_input(input)

# ╔═╡ 6cdac000-3549-11eb-19a4-8f7b31617365
test_grid = parse_input(test_input)

# ╔═╡ 2a9d0e9a-353b-11eb-25be-91691063ebe8
struct Slope
	right :: Int
	down :: Int
end

# ╔═╡ 692ac206-352c-11eb-145c-09870dd70a59
function number_of_trees_hit(grid,slope::Slope)
	l, w = size(grid)
	trees = 0
	for i in 1:(l÷slope.down)
		row = 1 + (i-1) * slope.down
		col = 1 + ((i-1)*slope.right) % w
		trees += grid[row,col]
	end
	trees
end

# ╔═╡ 980fd572-3539-11eb-337e-b576b6828534
@assert number_of_trees_hit(test_grid,Slope(3,1)) == 7

# ╔═╡ b9e4243c-3539-11eb-2756-b15d2871fd3a
number_of_trees_hit(grid,Slope(3,1))

# ╔═╡ 5f797e22-339d-11eb-1ff3-dda90f946838
md"## Part 2"

# ╔═╡ e3177eae-339e-11eb-3f9f-1f7230976832
md"*More Slopes*"

# ╔═╡ f76e45a8-3539-11eb-08e4-3d7a927a1f83
begin
	@assert number_of_trees_hit(test_grid,Slope(1,1)) == 2
	@assert number_of_trees_hit(test_grid,Slope(3,1)) == 7
	@assert number_of_trees_hit(test_grid,Slope(5,1)) == 3
	@assert number_of_trees_hit(test_grid,Slope(7,1)) == 4
	@assert number_of_trees_hit(test_grid,Slope(1,2)) == 2
end

# ╔═╡ 04005abc-353b-11eb-0a87-619740b6667f
slopes = [Slope(1,1), Slope(3,1), Slope(5,1), Slope(7,1), Slope(1,2)]

# ╔═╡ 5ca2b0f2-353b-11eb-088b-5567569dcbbe
typeof(slopes)

# ╔═╡ 58792402-353b-11eb-2307-651f8d88994a
part2(grid,slopes) = prod(number_of_trees_hit(grid,slope) for slope in slopes)

# ╔═╡ 72316a6a-353b-11eb-3ce1-3bd721f6136a
@assert part2(test_grid,slopes) == 336

# ╔═╡ 7afd4df0-353b-11eb-2be5-d3d305edeb40
part2(grid,slopes)

# ╔═╡ Cell order:
# ╟─3e0c149a-3350-11eb-0b17-7bea5d1a28d4
# ╠═2bf8e4c2-339b-11eb-1d49-7f1b210a650a
# ╟─18080fe2-339b-11eb-34f0-f393638506fc
# ╠═4e084a4e-339b-11eb-3ffc-41118215f226
# ╠═24e1613c-339b-11eb-2976-59884053d1ee
# ╠═3a854976-3531-11eb-11a3-9f976a15a484
# ╠═4f43d440-3531-11eb-146a-87dc6a097671
# ╠═677bfffc-3549-11eb-3a98-9b427086e291
# ╠═6cdac000-3549-11eb-19a4-8f7b31617365
# ╠═2a9d0e9a-353b-11eb-25be-91691063ebe8
# ╠═692ac206-352c-11eb-145c-09870dd70a59
# ╠═980fd572-3539-11eb-337e-b576b6828534
# ╠═b9e4243c-3539-11eb-2756-b15d2871fd3a
# ╠═5f797e22-339d-11eb-1ff3-dda90f946838
# ╠═e3177eae-339e-11eb-3f9f-1f7230976832
# ╠═f76e45a8-3539-11eb-08e4-3d7a927a1f83
# ╠═04005abc-353b-11eb-0a87-619740b6667f
# ╠═5ca2b0f2-353b-11eb-088b-5567569dcbbe
# ╠═58792402-353b-11eb-2307-651f8d88994a
# ╠═72316a6a-353b-11eb-3ce1-3bd721f6136a
# ╠═7afd4df0-353b-11eb-2be5-d3d305edeb40

max_stars_repo_name: RensDimmendaal/advent-of-code-julia

max_issues_repo_name: RensDimmendaal/advent-of-code-julia

max_forks_repo_name: RensDimmendaal/advent-of-code-julia


text:
module CoarseGraining3DPrisms
using LinearAlgebra

greet() = print("Hello World!")

export numchop, delta, qn, qnfact, visqrt, trian, RacWig6j, Fsymb, Gsymb, Gsymb2

numchop(num::Real)=(abs(num) >= 1000*eps() ? num : zero(Real))

function numchop(num::Complex)
    numchop(imag(num)) == 0 ? numchop(real(num)) : Complex(numchop(real(num)), numchop(imag(num)))
end

const K = 5
const x = K+1
const y = K/2

#Define delta_{ijkl} -> coupling rules
function delta(i::Float64,j::Float64,k::Float64)
    sol = 0
    if i <=(j+k) && j<=(i+k) && k<=(i+j) && i+j+k <= K && 2*(i+j+k)%2 ==0
        sol = 1
    end
    return sol
end

#Define quantum numbers qn (this is real)
function qn(n::Float64)
    sol = (exp(pi*n*im/(K+2)) - exp(-pi*n*im/(K+2))) / (exp(pi*im/(K+2)) - exp(-pi*im/(K+2)))
    return real(sol)
end

#Define qn factorial
function qnfact(n::Float64)
    sol = 1
    for i in 1:n
        sol *= qn(i)
    end
    return sol
end

#Define square root of quantum dimension
function visqrt(i::Float64)
    sol = ((-1+0im)^i )*sqrt(qn(2*i+1))
    return sol
end

#Define triangle equality
function trian(i::Float64,j::Float64,k::Float64)
    sol = 0
    if delta(i,j,k) == 1
        sol = delta(i,j,k)*sqrt(qnfact(i+j-k)*qnfact(i-j+k)*qnfact(-i+j+k)/qnfact(i+j+k+1))
    end
    return sol
end

# Define Racah-Wigner six-j symbol
function RacWig6j(i::Float64,j::Float64,m::Float64,k::Float64,l::Float64,n::Float64)
    a = i+j+m; b = i+l+n; c = k+j+n; d = k+l+m;  e = i+j+k+l; f = i+k+m+n; g = j+l+m+n
    sol = 0
    if delta(i,j,m) != 0 && delta(i,l,n) != 0 && delta(k,j,n) != 0 && delta(k,l,m) != 0
        sumz = 0
        for z in max(a,b,c,d):min(e,f,g)
            sumz += (-1)^z *qnfact(z+1)/
                ((qnfact(e-z)*qnfact(f-z)*qnfact(g-z))* (qnfact(z-a)*qnfact(z-b)*qnfact(z-c)*qnfact(z-d)))
        end
        sol = trian(i,j,m)*trian(i,l,n)*trian(k,j,n)*trian(k,l,m)*sumz
    end
    return sol
end

#Define F-symbol
function Fsymb(i::Float64,j::Float64,m::Float64,k::Float64,l::Float64,n::Float64)
    sol = 0
    if delta(i,j,m) != 0 && delta(i,l,n) != 0 && delta(k,j,n) != 0 && delta(k,l,m) != 0
        sol = (-1+0im)^(i+j+k+l)*sqrt(qn(2*m+1)*qn(2*n+1)) * RacWig6j(i,j,m,k,l,n)
    end
    return sol
end

#Define G-symbol
function Gsymb(i::Float64,j::Float64,m::Float64,k::Float64,l::Float64,n::Float64)
    sol = 0
    if delta(i,j,m) != 0 && delta(i,l,n) != 0 && delta(k,j,n) != 0 && delta(k,l,m) != 0
        sol = Fsymb(i,j,m,k,l,n) /(visqrt(m)*visqrt(n))
    end
    return sol
end

#Define G-symbol 2
function Gsymb2(i::Float64,j::Float64,m::Float64,k::Float64,l::Float64,n::Float64)
    sol = 0
    if delta(i,j,m) != 0 && delta(i,l,n) != 0 && delta(k,j,n) != 0 && delta(k,l,m) != 0
        sol = (-1+0im)^(i+j+k+l+m+n)*RacWig6j(i,j,m,k,l,n)
    end
    return sol
end

export prismA, prismB, prismB2, blocksPrA, blocksPrB, blockReduceA, blockReduceA1, blockReduceA2, svdA, svdB

#d1,d2,e1 is triangle we shall cut along
function prismA(ja1::Float64,jb1::Float64,jd1::Float64,je1::Float64,jf1::Float64,ja2::Float64,jb2::Float64,
		                 jd2::Float64,je2::Float64,jf2::Float64,jc::Float64,jg::Float64)
    sol = 0
    if delta(jd1,jd2,je1) != 0 && delta(jd1,jb1,jg) != 0 && delta(ja1,jd2,jg) != 0 && delta(ja1,jb1,je1) != 0 && delta(jd1,jc,jb2) != 0 && delta(jf1,jd2,jb2) != 0 && delta(jf1,jc,je1) != 0 && delta(jd1,ja2,jf2) != 0 && delta(jc,je2,jf2) != 0 && delta(jb2,je2,ja2) != 0
        dims = visqrt(ja1)*visqrt(jb1)*visqrt(jg)*visqrt(jf1)*visqrt(jc)*visqrt(jb2)*visqrt(je2)*visqrt(ja2)*visqrt(jf2)*visqrt(jd1)*visqrt(jd2)*visqrt(je1)
        sol =  dims*Gsymb(jd1,jd2,je1,ja1,jb1,jg) * Gsymb(jd1,jd2,je1,jf1,jc,jb2) * Gsymb(jd1,jc,jb2,je2,ja2,jf2)
    end
    return sol
end

function prismB(ja1::Float64,jb1::Float64,jd1p::Float64,je1::Float64,jf1::Float64,ja2::Float64,jb2::Float64,
		                 jd2p::Float64,je2::Float64,jf2::Float64,jc::Float64,jg::Float64)
    sol = 0
    for jd1 in 0.:0.5:y, jd2 in 0.:0.5:y
        sol += numchop(Fsymb(ja1,jg,jd2,jb2,jf1,jd2p)*Fsymb(jb1,jg,jd1,ja2,jf2,jd1p)*prismA(ja1,jb1,jd1,je1,jf1,ja2,jb2,jd2,je2,jf2,jc,jg))
    end
    return sol
end

function prismB2(ja1::Float64,jb1::Float64,jd1::Float64,je1::Float64,jf1::Float64,ja2::Float64,jb2::Float64,
		                 jd2::Float64,je2::Float64,jf2::Float64,jc::Float64,jg::Float64)
    return prismA(ja2,jb2,jd2,je2,jf2,ja1,jb1,jd1,je1,jf1,jc,jg)
end

#APrism(j5,j4,j7,j2,j1,jb,j3,j8,ja,j9,j6,jc)

function blocksPrA(Prsm::Array{Float64,12},jd1::Float64,jd2::Float64,je1::Float64)
    mat = zeros(x^3,x^6)
    row = 1
    col = 1
    for ja1 in 0.:0.5:y, jb1 in 0.:0.5:y, jg in 0.:0.5:y
        for jf1 in 0.:0.5:y, ja2 in 0.:0.5:y, jb2 in 0.:0.5:y, je2 in 0.:0.5:y, jf2 in 0.:0.5:y, jc in 0.:0.5:y
            mat[row,col] =  Prsm[Int(2*ja1+1),Int(2*jb1+1),Int(2*jd1+1),Int(2*je1+1),Int(2*jf1+1),Int(2*ja2+1),Int(2*jb2+1),Int(2*jd2+1),Int(2*je2+1),Int(2*jf2+1),Int(2*jc+1),Int(2*jg+1)]
            col += 1
        end
        col = 1
        row += 1
    end
    return mat
end

function blocksPrB(Prsm::Array{Float64,12},jd2::Float64,jd1::Float64,je2::Float64)
    mat = zeros(x^3,x^6)
    row = 1
    col = 1
    for ja2 in 0.:0.5:y, jb2 in 0.:0.5:y, jg in 0.:0.5:y
        for jf2 in 0.:0.5:y, ja1 in 0.:0.5:y, jb1 in 0.:0.5:y, je1 in 0.:0.5:y, jf1 in 0.:0.5:y, jc in 0.:0.5:y
            mat[row,col] =  Prsm[Int(2*ja1+1),Int(2*jb1+1),Int(2*jd1+1),Int(2*je1+1),Int(2*jf1+1),Int(2*ja2+1),Int(2*jb2+1),Int(2*jd2+1),Int(2*je2+1),Int(2*jf2+1),Int(2*jc+1),Int(2*jg+1)]
            col += 1
        end
        col = 1
        row += 1
    end
    return mat
end

function blockReduceA(ampJ::Array{Any,1},jd1::Float64,jd2::Float64,je1::Float64)
    #a = Float64[]
    #b = Array{Int64,1}[]
    #ampInfo = [] # will contain indices and amplitudes
    # position of d1 = 3, d2= 8 and e1 = 4
    ampInfo = ampJ[findall(x-> x[1][3] == jd1 && x[1][8] == jd2 && x[1][4] == je1 , ampJ)]
    nzcol = Array{Float64,1}[]
    #for ja1 in 0.:0.5:y, jb1 in 0.:0.5:y, jg in 0.:0.5:y
    nzrow = Array{Float64,1}[]
    #hh = unique!([ [ampInfo[i][1][1],ampInfo[i][1][2],ampInfo[i][1][12]] for 1:length(ampInfo) ])
    for i in 1:length(ampInfo) # position of a1 = 1, b1 = 2 g = 12
        push!(nzcol,[ampInfo[i][1][1],ampInfo[i][1][2],ampInfo[i][1][12]])
        push!(nzrow,[ampInfo[i][1][5],ampInfo[i][1][6],ampInfo[i][1][7],ampInfo[i][1][9],ampInfo[i][1][10],ampInfo[i][1][11]])
    end
    unique!(nzcol)
    unique!(nzrow)
    dd = zeros(length(nzcol),length(nzrow))
    if length(dd) != 0
        for i in 1:length(ampInfo)
            qa = findall(x-> x == [ampInfo[i][1][1], ampInfo[i][1][2], ampInfo[i][1][12]],nzcol)
            qb = findall(x-> x == [ampInfo[i][1][5], ampInfo[i][1][6], ampInfo[i][1][7], ampInfo[i][1][9], ampInfo[i][1][10], ampInfo[i][1][11] ], nzrow)
            dd[qa,qb] = [ampInfo[i][2]]
        end
        return nzcol,nzrow, dd
    else
        return 0,0,0
    end
end

function svdA(prA::Array{Float64,12},ja::Float64,jb::Float64,jc::Float64)# parameters d1,d2,e1
    mat = blocksPrA(prA,ja,jb,jc)
    U, s, V = svd(mat)
    #V = ctranspose(V)
    return (U,s,V)
end

function svdB(prB::Array{Float64,12},ja::Float64,jb::Float64,jc::Float64)# parameters d1,d2,e2
    mat = blocksPrB(prB,ja,jb,jc)
    U, s, V = svd(mat)
    return (U,s,V)
end

export UtensorA, UtensorB, VtensorA, VtensorB

function UtensorA(prA::Array{Float64,12},ja::Float64,jb::Float64,jc::Float64)# parameters d1,d2,e1
    UT = zeros(x,x,x)
    U = svdA(prA,ja,jb,jc)[1]
    s = svdA(prA,ja,jb,jc)[2]
    index = 1
    for ja1 in 0.:0.5:y, jb1 in 0.:0.5:y, jg in 0.:0.5:y
        UT[Int(2*ja1+1),Int(2*jb1+1),Int(2*jg+1)] = numchop(U[index,1]*sqrt(s[1]))
        index += 1
    end
    return UT
end

function UtensorB(prB::Array{Float64,12},ja::Float64,jb::Float64,jc::Float64)# parameters d1,d2,e1
    UT = zeros(x,x,x)
    U = svdB(prB,ja,jb,jc)[1]
    s = svdB(prB,ja,jb,jc)[2]
    index = 1
    for ja2 in 0.:0.5:y, jb2 in 0.:0.5:y, jg in 0.:0.5:y
        UT[Int(2*ja2+1),Int(2*jb2+1),Int(2*jg+1)] = numchop(U[index,1]*sqrt(s[1]))
        index += 1
    end
    return UT
end

#VVtensors
function VtensorA(prA::Array{Float64,12},ja::Float64,jb::Float64,jc::Float64)# parameters d1,d2,e1
    VT = zeros(x,x,x,x,x,x)
    V = svdA(prA,ja,jb,jc)[3]
    s = svdA(prA,ja,jb,jc)[2]
    index = 1
    for jf1 in 0.:0.5:y, ja2 in 0.:0.5:y, jb2 in 0.:0.5:y, je2 in 0.:0.5:y, jf2 in 0.:0.5:y, jc in 0.:0.5:y
        VT[Int(2*jf1+1),Int(2*ja2+1),Int(2*jb2+1),Int(2*je2+1),Int(2*jf2+1),Int(2*jc+1)] = numchop(V[index,1]*sqrt(s[1]))
        index += 1
    end
    return VT
end

function VtensorB(prB::Array{Float64,12},ja::Float64,jb::Float64,jc::Float64)# parameters d1,d2,e1
    VT = zeros(x,x,x,x,x,x)
    V = svdB(prB,ja,jb,jc)[3]
    s = svdB(prB,ja,jb,jc)[2]
    index = 1
    for jf2 in 0.:0.5:y, ja1 in 0.:0.5:y, jb1 in 0.:0.5:y, je1 in 0.:0.5:y, jf1 in 0.:0.5:y, jc in 0.:0.5:y
        VT[Int(2*jf2+1),Int(2*ja1+1),Int(2*jb1+1),Int(2*je1+1),Int(2*jf1+1),Int(2*jc+1)] = numchop(V[index,1]*sqrt(s[1]))
        index += 1
    end
    return VT
end

export  fullUtensorA,  fullUtensorB, fullVtensorA, fullVtensorA

function fullUtensorA(prA::Array{Float64,12})
    fullUT = zeros(x,x,x,x,x,x)
    for jd1 in 0.:0.5:y, jd2 in 0.:0.5:y, je1 in 0.:0.5:y
        fullUT[Int(2*jd1+1),Int(2*jd2+1),Int(2*je1+1),:,:,:] = UtensorA(prA,jd1,jd2,je1)
    end
    return fullUT
end

function fullVtensorA(prA::Array{Float64,12})
    fullVT = zeros(x,x,x,x,x,x,x,x,x)
    for jd1 in 0.:0.5:y, jd2 in 0.:0.5:y, je1 in 0.:0.5:y
        fullVT[Int(2*jd1+1),Int(2*jd2+1),Int(2*je1+1),:,:,:,:,:,:] = VtensorA(prA,jd1,jd2,je1)
    end
    return fullVT
end

function fullUtensorB(prB::Array{Float64,12})
    fullUT = zeros(x,x,x,x,x,x)
    for jd1 in 0.:0.5:y, jd2 in 0.:0.5:y, je2 in 0.:0.5:y
        fullUT[Int(2*jd1+1),Int(2*jd2+1),Int(2*je2+1),:,:,:] = UtensorB(prB,jd1,jd2,je2)
    end
    return fullUT
end

function fullVtensorB(prB::Array{Float64,12})
    fullVT = zeros(x,x,x,x,x,x,x,x,x)
    for jd1 in 0.:0.5:y, jd2 in 0.:0.5:y, je2 in 0.:0.5:y
        fullVT[Int(2*jd1+1),Int(2*jd2+1),Int(2*je2+1),:,:,:,:,:,:] = VtensorB(prB,jd1,jd2,je2)
    end
    return fullVT
end

export prsmUAVB, prsmUBVA, prsmEff

function prsmUAVB(prA::Array{Float64,12},prB::Array{Float64,12})
    UA = fullUtensorA(prA)
    VB = fullVtensorB(prB)
    prsm = zeros(x,x,x,x,x,x,x,x,x,x,x,x)
    for jb1 in 0.:0.5:y, jc1 in 0.:0.5:y, jd1 in 0.:0.5:y, jh1 in 0.:0.5:y, je11 in 0.:0.5:y, jf12 in 0.:0.5:y, ja2 in 0.:0.5:y, jf2 in 0.:0.5:y, jh2 in 0.:0.5:y, je21 in 0.:0.5:y, jg in 0.:0.5:y, jn in 0.:0.5:y
        ans = UA[Int(2*jb1+1),Int(2*jd1+1),Int(2*jg+1),Int(2*jn+1),Int(2*jh1+1),Int(2*je11+1)] *
                VB[Int(2*jd1+1),Int(2*ja2+1),Int(2*jf2+1),Int(2*je21+1),Int(2*jn+1),Int(2*jh2+1),Int(2*jf12+1),Int(2*je11+1),Int(2*jc1+1)]
        prsm[Int(2*jb1+1),Int(2*jc1+1),Int(2*jd1+1),Int(2*jh1+1),Int(2*je11+1),Int(2*jf12+1),Int(2*ja2+1),Int(2*jf2+1),Int(2*jh2+1),Int(2*je21+1),Int(2*jg+1),Int(2*jn+1)] = ans
    end
    return prsm
end

function prsmUBVA(prA::Array{Float64,12},prB::Array{Float64,12})
    UB = fullUtensorB(prB)
    VA = fullVtensorA(prA)
    prsm = Array{Float64}(x,x,x,x,x,x,x,x,x,x,x,x)
    for jb1 in 0.:0.5:y, jc1 in 0.:0.5:y, jd1 in 0.:0.5:y, jh1 in 0.:0.5:y, je11 in 0.:0.5:y, jf12 in 0.:0.5:y, ja2 in 0.:0.5:y, jf2 in 0.:0.5:y, jh2 in 0.:0.5:y, je21 in 0.:0.5:y, jg in 0.:0.5:y, jn in 0.:0.5:y
        ans = UB[Int(2*jb1+1),Int(2*jd1+1),Int(2*jg+1),Int(2*jn+1),Int(2*jh1+1),Int(2*je11+1)] *
                VA[Int(2*jd1+1),Int(2*ja2+1),Int(2*jf2+1),Int(2*je21+1),Int(2*jn+1),Int(2*jh2+1),Int(2*jf12+1),Int(2*je11+1),Int(2*jc1+1)]
        prsm[Int(2*jb1+1),Int(2*jc1+1),Int(2*jd1+1),Int(2*jh1+1),Int(2*je11+1),Int(2*jf12+1),Int(2*ja2+1),Int(2*jf2+1),Int(2*jh2+1),Int(2*je21+1),Int(2*jg+1),Int(2*jn+1)] = ans
    end
    return prsm
end

function prsmEff(prA::Array{Float64,12},prB::Array{Float64,12})
    pr1 = prsmUAVB(prA,prB)
    pr2 = prsmUBVA(prA,prB)
    prsm = zeros(x,x,x,x,x,x,x,x,x,x,x,x,x,x,x,x,x,x)
    for ja1 in 0.:0.5:y, jb1 in 0.:0.5:y, jc1 in 0.:0.5:y, jd1 in 0.:0.5:y, je11 in 0.:0.5:y, je12 in 0.:0.5:y,jf1 in 0.:0.5:y, jh1 in 0.:0.5:y, ja2 in 0.:0.5:y, jb2 in 0.:0.5:y, jc2 in 0.:0.5:y, jd2 in 0.:0.5:y, je21 in 0.:0.5:y,je22 in 0.:0.5:y, jf2 in 0.:0.5:y, jh2 in 0.:0.5:y, jf12 in 0.:0.5:y, jg in 0.:0.5:y
        sum = 0
        for jn in 0.:0.5:y
            ans1 = pr1[Int(2*jb1+1),Int(2*jc1+1),Int(2*jd1+1),Int(2*jh1+1),Int(2*je11+1),Int(2*jf12+1),Int(2*ja2+1),Int(2*jf2+1),Int(2*jh2+1),Int(2*je21+1),Int(2*jg+1),Int(2*jn+1)]
            ans2 = pr2[Int(2*jb2+1),Int(2*jc2+1),Int(2*jd2+1),Int(2*jh2+1),Int(2*je22+1),Int(2*jf12+1),Int(2*ja1+1),Int(2*jf1+1),Int(2*jh1+1),Int(2*je12+1),Int(2*jg+1),Int(2*jn+1)]
            sum+= (ans1*ans2)
        end
        prsm[Int(2*jb1+1),Int(2*jc1+1),Int(2*jd1+1),Int(2*je11+1),Int(2*ja2+1),Int(2*jf2+1),Int(2*je21+1),Int(2*jb2+1),Int(2*jc2+1),Int(2*jd2+1),Int(2*je12+1),Int(2*ja1+1),Int(2*jf1+1),Int(2*je22+1),Int(2*jh1+1),Int(2*jf12+1),Int(2*jh2+1),Int(2*jg+1)] = sum
    end
    return prsm
end

end # module


text:
######################################################################################################################
# This file runs (and optimises) RICE in its Mimi implementation.
# The only difference to Excel is that the damage deflator is now 1/(1+D)
# Some changes in allocation of variables and parameters across components have been implemented relative to MimiRICE2010
##############################################################################################################

# Load required Julia packages.
using NLopt # for optimisation

# Load RICE+AIR source code.
include("src/mRICE2010.jl")

# this is how you initiate an instance m of the model
m = get_rice() #default if get_rice(objective = "Neutral") set objective = "Negishi" for that objective
# then you run the model
run(m)

# say you want to change some parameter
#like the discount rate
m[:welfare,:rho] = 1.5
#or the mitigation fraction, this is how we input the control variable in the optimisation
m[:emissions,:MIU] = 0.5*ones(60,12)

# the way the optimisation works is that we construc the model, and set a bunch of user defined paramaters, and from that construct the objective

mod = get_rice()

# number of periods of control
t_choice = 29
# Maximum time in seconds to run local optimization (in case optimization does not converge).
local_stop_time = 500
# Relative tolerance criteria for global optimization convergence (will stop if |Δf| / |f| < tolerance from one iteration to the next.)
local_tolerance = 1e-12
objective = construct_RICE_objective(m,t_choice)
constraint = retConstraint(m,t_choice)
# set up optimisation# Create an NLopt optimization object.
opt_object = Opt(:LN_SBPLX, t_choice*12)
# set up constraint

# bounds on the control variable
lower_bounds!(opt_object, zeros(12*t_choice))
upper_bounds!(opt_object, ones(12*t_choice))
# Set maximum run time.
maxtime!(opt_object, local_stop_time)
# Set convergence tolerance.
ftol_rel!(opt_object, local_tolerance)
# Set objective function.
max_objective!(opt_object, (x, grad) -> objective(x))
inequality_constraint!(ot_object, (x, grad) -> constraint(x)-0.5)
max_welfare, optimal_rates, convergence_flag = optimize(opt_object, 0.5*ones(12*t_choice))

max_stars_repo_name: fdennig/RICEupdate

max_issues_repo_name: fdennig/RICEupdate

max_forks_repo_name: fdennig/RICEupdate


text:
const StieltjesPoint{T,V,D} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{BroadcastQuasiMatrix{T,typeof(-),Tuple{T,QuasiAdjoint{V,Inclusion{V,D}}}}}}
const ConvKernel{T,D} = BroadcastQuasiMatrix{T,typeof(-),Tuple{D,QuasiAdjoint{T,D}}}
const Hilbert{T,D} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{ConvKernel{T,Inclusion{T,D}}}}
const LogKernel{T,D} = BroadcastQuasiMatrix{T,typeof(log),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{ConvKernel{T,Inclusion{T,D}}}}}}
const PowKernel{T,D,F<:Real} = BroadcastQuasiMatrix{T,typeof(^),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{ConvKernel{T,Inclusion{T,D}}}},F}}

@simplify function *(S::StieltjesPoint{<:Any,<:Any,<:ChebyshevInterval}, wT::WeightedBasis{<:Any,<:ChebyshevTWeight,<:ChebyshevT})
    w,T = wT.args
    J = jacobimatrix(T)
    z, x = parent(S).args[1].args
    transpose((J'-z*I) \ [-π; zeros(∞)])
end

@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval}, wT::WeightedBasis{<:Any,<:ChebyshevTWeight,<:ChebyshevT})
    T = promote_type(eltype(H), eltype(wT))
    ApplyQuasiArray(*, ChebyshevU{T}(), _BandedMatrix(Fill(-convert(T,π),1,∞), ∞, -1, 1))
end

@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval}, wU::WeightedBasis{<:Any,<:ChebyshevUWeight,<:ChebyshevU})
    T = promote_type(eltype(H), eltype(wU))
    ApplyQuasiArray(*, ChebyshevT{T}(), _BandedMatrix(Fill(convert(T,π),1,∞), ∞, 1, -1))
end

###
# LogKernel
###

@simplify function *(L::LogKernel{<:Any,<:ChebyshevInterval}, wT::WeightedBasis{<:Any,<:ChebyshevTWeight,<:ChebyshevT})
    T = promote_type(eltype(L), eltype(wT))
    ApplyQuasiArray(*, ChebyshevT{T}(), Diagonal(Vcat(-π*log(2*one(T)),-convert(T,π)./(1:∞))))
end

###
# PowKernel
###

@simplify function *(K::PowKernel{<:Any,<:ChebyshevInterval}, wT::WeightedBasis{<:Any,<:JacobiWeight,<:Jacobi})
    T = promote_type(eltype(K), eltype(wT))
    cnv,α = K.args
    x,y = K.args[1].args[1].args
    @assert x' == y
    β = (-α-1)/2
    ApplyQuasiArray(*, ChebyshevT{T}(), Diagonal(1:∞))
end

####
# StieltjesPoint
####

@simplify function *(S::StieltjesPoint, wT::SubQuasiArray{<:Any,2,<:WeightedBasis,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
    P = parent(wT)
    z, x = parent(S).args[1].args
    z̃ = inbounds_getindex(parentindices(wT)[1], z)
    x̃ = axes(P,1)
    (inv.(z̃ .- x̃') * P)[:,parentindices(wT)[2]]
end

@simplify function *(H::Hilbert, wT::SubQuasiArray{<:Any,2,<:WeightedBasis,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
    P = parent(wT)
    x = axes(P,1)
    apply(*, inv.(x .- x'), P)[parentindices(wT)...]
end

@simplify function *(L::LogKernel, wT::SubQuasiArray{<:Any,2,<:WeightedBasis,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}})
    V = promote_type(eltype(L), eltype(wT))
    P = parent(wT)
    kr, jr = parentindices(wT)
    @assert P isa WeightedBasis{<:Any,<:ChebyshevWeight,<:Chebyshev}
    x = axes(P,1)
    w,T = P.args
    D = T \ apply(*, log.(abs.(x .- x')), P)
    c = inv(2*kr.A)
    T[kr,:] * Diagonal(Vcat(2*convert(V,π)*c*log(c), 2c*D.diag.args[2]))
end


text:
using CSV
using DataFrames

function main()
    df = CSV.File("7days.csv") |> DataFrame
    sal = reshape(df.total, nrow(df), 1)
    value = join(string.(sal), ", ")
    println(value)
    println(last(df, 5))
end

main()

max_stars_repo_name: dmolina/DaemonMode

max_issues_repo_name: dmolina/DaemonMode

max_forks_repo_name: dmolina/DaemonMode


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: elastic_inference
using AWS.Compat
using AWS.UUIDs

"""
    DescribeAcceleratorOfferings()

 Describes the locations in which a given accelerator type or set of types is present in a given region.

# Required Parameters
- `locationType`:  The location type that you want to describe accelerator type offerings for. It can assume the following values: region: will return the accelerator type offering at the regional level. availability-zone: will return the accelerator type offering at the availability zone level. availability-zone-id: will return the accelerator type offering at the availability zone level returning the availability zone id.

# Optional Parameters
- `acceleratorTypes`:  The list of accelerator types to describe.
"""
describe_accelerator_offerings(locationType; aws_config::AWSConfig=global_aws_config()) = elastic_inference("POST", "/describe-accelerator-offerings", Dict{String, Any}("locationType"=>locationType); aws_config=aws_config)
describe_accelerator_offerings(locationType, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = elastic_inference("POST", "/describe-accelerator-offerings", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("locationType"=>locationType), args)); aws_config=aws_config)

"""
    DescribeAcceleratorTypes()

 Describes the accelerator types available in a given region, as well as their characteristics, such as memory and throughput.

"""
describe_accelerator_types(; aws_config::AWSConfig=global_aws_config()) = elastic_inference("GET", "/describe-accelerator-types"; aws_config=aws_config)
describe_accelerator_types(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = elastic_inference("GET", "/describe-accelerator-types", args; aws_config=aws_config)

"""
    DescribeAccelerators()

 Describes information over a provided set of accelerators belonging to an account.

# Optional Parameters
- `acceleratorIds`:  The IDs of the accelerators to describe.
- `filters`:  One or more filters. Filter names and values are case-sensitive. Valid filter names are: accelerator-types: can provide a list of accelerator type names to filter for. instance-id: can provide a list of EC2 instance ids to filter for.
- `maxResults`:  The total number of items to return in the command's output. If the total number of items available is more than the value specified, a NextToken is provided in the command's output. To resume pagination, provide the NextToken value in the starting-token argument of a subsequent command. Do not use the NextToken response element directly outside of the AWS CLI.
- `nextToken`:  A token to specify where to start paginating. This is the NextToken from a previously truncated response.
"""
describe_accelerators(; aws_config::AWSConfig=global_aws_config()) = elastic_inference("POST", "/describe-accelerators"; aws_config=aws_config)
describe_accelerators(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = elastic_inference("POST", "/describe-accelerators", args; aws_config=aws_config)

"""
    ListTagsForResource()

 Returns all tags of an Elastic Inference Accelerator.

# Required Parameters
- `resourceArn`:  The ARN of the Elastic Inference Accelerator to list the tags for.

"""
list_tags_for_resource(resourceArn; aws_config::AWSConfig=global_aws_config()) = elastic_inference("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = elastic_inference("GET", "/tags/$(resourceArn)", args; aws_config=aws_config)

"""
    TagResource()

 Adds the specified tags to an Elastic Inference Accelerator.

# Required Parameters
- `resourceArn`:  The ARN of the Elastic Inference Accelerator to tag.
- `tags`:  The tags to add to the Elastic Inference Accelerator.

"""
tag_resource(resourceArn, tags; aws_config::AWSConfig=global_aws_config()) = elastic_inference("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = elastic_inference("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

 Removes the specified tags from an Elastic Inference Accelerator.

# Required Parameters
- `resourceArn`:  The ARN of the Elastic Inference Accelerator to untag.
- `tagKeys`:  The list of tags to remove from the Elastic Inference Accelerator.

"""
untag_resource(resourceArn, tagKeys; aws_config::AWSConfig=global_aws_config()) = elastic_inference("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = elastic_inference("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)