text:
###
# Authors: Daniel Thuerck
# NEC Laboratories Europe GmbH, Copyright (c) 2021, All rights reserved.
#
# THIS HEADER MAY NOT BE EXTRACTED OR MODIFIED IN ANY WAY.
#
# SOFTWARE LICENSE AGREEMENT
# ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
# BY USING OR DOWNLOADING THE SOFTWARE, YOU ARE AGREEING TO THE TERMS
# OF THIS LICENSE AGREEMENT. IF YOU DO NOT AGREE WITH THESE TERMS, YOU MAY
# NOT USE OR DOWNLOAD THE SOFTWARE.
#
# This is a license agreement ("Agreement") between your academic institution or non-profit
# organization or self (called "Licensee" or "You" in this Agreement) and NEC Laboratories Europe
# GmbH (called "Licensor" in this Agreement). All rights not specifically granted to you in this
# Agreement are reserved for Licensor.
#
# RESERVATION OF OWNERSHIP AND GRANT OF LICENSE: Licensor retains exclusive
# ownership of any copy of the Software (as defined below) licensed under this Agreement and hereby
# grants to Licensee a personal, non-exclusive, non-transferable license to use the Software for
# noncommercial research purposes, without the right to sublicense, pursuant to the terms and
# conditions of this Agreement. NO EXPRESS OR IMPLIED LICENSES TO ANY OF LICENSOR’S
# PATENT RIGHTS ARE GRANTED BY THIS LICENSE. As used in this Agreement, the term
# "Software" means (i) the actual copy of all or any portion of code for program routines made
# accessible to Licensee by Licensor pursuant to this Agreement, inclusive of backups, updates, and/or
# merged copies permitted hereunder or subsequently supplied by Licensor, including all or any file
# structures, programming instructions, user interfaces and screen formats and sequences as well as any
# and all documentation and instructions related to it, and (ii) all or any derivatives and/or modifications
# created or made by You to any of the items specified in (i).
#
# CONFIDENTIALITY/PUBLICATIONS: Licensee acknowledges that the Software is proprietary to
# Licensor, and as such, Licensee agrees to receive all such materials and to use the Software only in
# accordance with the terms of this Agreement. Licensee agrees to use reasonable effort to protect the
# Software from unauthorized use, reproduction, distribution, or publication. All publication materials
# mentioning features or use of this software must explicitly include an acknowledgement the software
# was developed by NEC Laboratories Europe GmbH.
#
# COPYRIGHT: The Software is owned by Licensor.
#
# PERMITTED USES: The Software may be used for your own noncommercial internal research
# purposes. You understand and agree that Licensor is not obligated to implement any suggestionsand/or
# feedback you might provide regarding the Software, but to the extent Licensor does so, you are
# not entitled to any compensation related thereto.
#
# DERIVATIVES: You may create derivatives of or make modifications to the Software, however, You
# agree that all and any such derivatives and modifications will be owned by Licensor and become a part
# of the Software licensed to You under this Agreement. You may only use such derivatives and
# modifications for your own noncommercial internal research purposes, and you may not otherwise
# use, distribute or copy such derivatives and modifications in violation of this Agreement.
#
# BACKUPS: If Licensee is an organization, it may make that number of copies of the Software
# necessary for internal noncommercial use at a single site within its organization provided that all
# information appearing in or on the original labels, including the copyright and trademark notices are
# copied onto the labels of the copies.
#
# USES NOT PERMITTED: You may not distribute, copy or use the Software except as explicitly
# permitted herein. Licensee has not been granted any trademark license as part of this Agreement.
# Neither the name of NEC Laboratories Europe GmbH nor the names of its contributors may be used to
# endorse or promote products derived from this Software without specific prior written permission.
# You may not sell, rent, lease, sublicense, lend, time-share or transfer, in whole or in part, or provide
# third parties access to prior or present versions (or any parts thereof) of the Software.
#
# ASSIGNMENT: You may not assign this Agreement or your rights hereunder without the prior
# written consent of Licensor. Any attempted assignment without such consent shall be null and void.
#
# TERM: The term of the license granted by this Agreement is from Licensee's acceptance of this
# Agreement by downloading the Software or by using the Software until terminated as provided below.
# The Agreement automatically terminates without notice if you fail to comply with any provision of
# this Agreement. Licensee may terminate this Agreement by ceasing using the Software. Upon any
# termination of this Agreement, Licensee will delete any and all copies of the Software. You agree that
# all provisions which operate to protect the proprietary rights of Licensor shall remain in force should
# breach occur and that the obligation of confidentiality described in this Agreement is binding in
# perpetuity and, as such, survives the term of the Agreement.
#
# FEE: Provided Licensee abides completely by the terms and conditions of this Agreement, there is no
# fee due to Licensor for Licensee's use of the Software in accordance with this Agreement.
#
# DISCLAIMER OF WARRANTIES: THE SOFTWARE IS PROVIDED "AS-IS" WITHOUT
# WARRANTY OF ANY KIND INCLUDING ANY WARRANTIES OF PERFORMANCE OR
# MERCHANTABILITY OR FITNESS FOR A PARTICULAR USE OR PURPOSE OR OF NON-
# INFRINGEMENT. LICENSEE BEARS ALL RISK RELATING TO QUALITY AND
# PERFORMANCE OF THE SOFTWARE AND RELATED MATERIALS.
#
# SUPPORT AND MAINTENANCE: No Software support or training by the Licensor is provided as
# part of this Agreement.
#
# EXCLUSIVE REMEDY AND LIMITATION OF LIABILITY: To the maximum extent permitted
# under applicable law, Licensor shall not be liable for direct, indirect, special, incidental, or
# consequential damages or lost profits related to Licensee's use of and/or inability to use the Software,
# even if Licensor is advised of the possibility of such damage.
#
# EXPORT REGULATION: Licensee agrees to comply with any and all applicable export control laws,
# regulations, and/or other laws related to embargoes and sanction programs administered by law.
#
# SEVERABILITY: If any provision(s) of this Agreement shall be held to be invalid, illegal, or
# unenforceable by a court or other tribunal of competent jurisdiction, the validity, legality and
# enforceability of the remaining provisions shall not in any way be affected or impaired thereby.
#
# NO IMPLIED WAIVERS: No failure or delay by Licensor in enforcing any right or remedy under this
# Agreement shall be construed as a waiver of any future or other exercise of such right or remedy by
# Licensor.
#
# GOVERNING LAW: This Agreement shall be construed and enforced in accordance with the laws of
# Germany without reference to conflict of laws principles. You consent to the personal jurisdiction of
# the courts of this country and waive their rights to venue outside of Germany.
#
# ENTIRE AGREEMENT AND AMENDMENTS: This Agreement constitutes the sole and entire
# agreement between Licensee and Licensor as to the matter set forth herein and supersedes any
# previous agreements, understandings, and arrangements between the parties relating hereto.
#
# THIS HEADER MAY NOT BE EXTRACTED OR MODIFIED IN ANY WAY.
###

# list of symbols for primitive types
supported_continuous_type_symbols = [
    :Float16,
    :Float32,
    :Float64
]
supported_int_type_symbols = [
    :Int128,
    :Int16,
    :Int32,
    :Int64,
    :Int8,
    :UInt128,
    :UInt16,
    :UInt32,
    :UInt64,
    :UInt8,
    :Bool
]
unsupported_type_symbols = [
    :Integer,
    :Signed,
    :Unsigned,
    :Complex,
    :Real,
    :Number,
    :AbstractFloat,
    :BigFloat,
    :BigInt
]

# simpicistic type hierachy (just mentions, not, e.g., struct definitions)
abstract type type_t end
abstract type primitive_t <: type_t end
abstract type composite_t <: type_t end

mutable struct continuous_t <: primitive_t
    bit_width
    complex
end

mutable struct integer_t <: primitive_t
    bit_width
    unsigned
end

mutable struct void_t <: primitive_t
end

mutable struct ptr_t <: composite_t
    pointee_t :: type_t
end

mutable struct identifier_t
    name :: Symbol
    parts :: Vector{Symbol}
    fqname
end

mutable struct struct_ref_t <: composite_t
    name :: identifier_t
end

mutable struct array_t <: composite_t
    pointee_t :: type_t
    dims
end

mutable struct param
    name
    type :: type_t
    managed :: Bool

    replacement
end

mutable struct tuple_t <: composite_t
    members :: Array{param}
end

# structs are defined seperately (other then tuples, they are not anonymous)
mutable struct struct_def <: composite_t
    name :: identifier_t
    members :: Array{param}
end

# dumping back to expr
function expr_from_type(t::continuous_t)
    return Symbol("Float" * string(t.bit_width))
end

function expr_from_type(t::integer_t)
    if t.bit_width == 1
        return Symbol("Bool")
    end
    sign = ""
    if t.unsigned
        sign = "U"
    end
    return Symbol(sign * "Int" * string(t.bit_width))
end

function expr_from_type(t::void_t)
    return Symbol("Nothing")
end

function expr_from_type(t::ptr_t)
    return Expr(
        Symbol("curly"),
        Symbol("Ptr"),
        expr_from_type(t.pointee_t)
    )
end

function expr_from_type(t::struct_ref_t)
    return t.name.name
end

function expr_from_type(t::array_t)
    args_array = []
    if t.dims == 1
        push!(args_array, Symbol("Vector"))
    else
        push!(args_array, Symbol("Array"))
    end

    push!(args_array, expr_from_type(t.pointee_t))

    if t.dims > 1
        push!(args_array, t.dims)
    end

    return Expr(
        Symbol("curly"),
        args_array...
    )
end

function expr_from_type(t::tuple_t)
    member_args = []
    for m in t.members
        push!(member_args, expr_from_type(m))
    end

    return Expr(
        Symbol("curly"),
        Symbol("Tuple"),
        member_args...
    )
end

# mapping to C-Type
function cstring_from_type(t::continuous_t)
    if t.bit_width == 16
        return "half"
    elseif t.bit_width == 32
        return "float"
    end

    return "double"
end

function cstring_from_type(t::integer_t)
    s = ""

    if t.bit_width == 8
        s *= "char"
    elseif t.bit_width == 16
        s *= "short"
    elseif t.bit_width == 32
        s *= "int"
    else
        s *= "long long"
    end

    if t.unsigned
        s = "unsigned " * s
    end

    return s
end

function cstring_from_type(t::void_t)
    return "void"
end

function cstring_from_type(t::ptr_t)
    return cstring_from_type(t.pointee_t) * " *"
end

function cstring_from_type(t::struct_ref_t)
    return string(t.name.name)
end

function cstring_from_type(t::array_t)
    return cstring_from_type(t.pointee_t) * "[]"
end

function cstring_from_type(t::tuple_t)
    # tuples must be mapped to anonymous structs
    return hash_type(t)
end

# handles names with "." in it
function parse_id(expr :: Symbol)
    return identifier_t(expr, [expr], expr)
end

function name_parts_from_expr(expr)
    if isa(expr, Symbol)
        return [expr]
    elseif isa(expr, QuoteNode)
        return [expr.value]
    elseif expr.head == Symbol(".")
        names = []

        for part in expr.args
            append!(
                names,
                name_parts_from_expr(part)
            )
        end

        return names
    else
        error("Unknown name expr: ", expr)
    end
end

function parse_id(expr :: Expr)
    # get parts
    parts = name_parts_from_expr(expr)

    # FQN = expr of parts seperated by "."
    function expr_from_fqname(name)
        if length(name) == 1
            return name[1]
        else
            return Expr(
                Symbol("."),
                expr_from_fqname(name[1:(length(name) - 1)]),
                QuoteNode(name[length(name)])
            )
        end
    end

    fqname = expr_from_fqname(parts)

    # "use" name: replace "." by _ in a symbol
    use_name = string(parts[1])
    for p in parts[2:length(parts)]
        use_name *= "_" * string(p)
    end
    use_name = Symbol(use_name)

    return identifier_t(use_name, parts, fqname)
end

# descriptive strings for types
function hash_type(t::continuous_t)
    return ("continuous_" * string(t.bit_width) * (t.complex ? "_1" : "_0"))
end

function hash_type(t::integer_t)
    return ("integer_" * string(t.bit_width) * (t.unsigned ? "_1" : "_0"))
end

function hash_type(t::void_t)
    return "void"
end

function hash_type(t::ptr_t)
    return ("ptr{" * hash_type(t.pointee_t) * "}")
end

function hash_type(t::struct_ref_t)
    return t.name.name
end

function hash_type(t::array_t)
    return ("array{" * hash_type(t.pointee_t) * "_" * string(t.dims) * "}")
end

function hash_type(t::tuple_t)
    str = "tuple{"
    for (ix, elem) in enumerate(t.members)
        if ix > 1
            str = str * "_"
        end
        str = str * hash_type(elem)
    end

    return (str * "}")
end

function hash_type(t::param)
    return hash_type(t.type)
end

function hash_type(t::struct_def)
    hash = t.name[0]
    for n in t.name[2:length(t.name)]
        hash *= "." * string(n)
    end

    return hash
end

max_stars_repo_name: nec-research/jCall

max_issues_repo_name: nec-research/jCall

max_forks_repo_name: nec-research/jCall


text:
module TestPaging

  using GBIF
  using Test

  # No query
  set = occurrences()
  next!(set)
  @test length(set) == 40
  set.query["limit"] = 40
  next!(set)
  @test length(set) == 80

  # Queries
  qpars = Dict{String,Any}(
    "hasCoordinate" => true,
    "q" => "Lamellodiscus"
  )
  lam = occurrences(qpars)
  exp_count = lam.count
  complete!(lam)
  @test length(lam) == exp_count

end

ext: jl


text:
using Documenter, DM

makedocs(
    modules = [DM],
    clean = false,
    sitename="DM.jl",
    pages= [
        "Home" => "index.md",
    ],
)

deploydocs(
    repo = "github.com/neversakura/DM.jl.git",
    target = "build"
)


text:
using SinGAN

using Flux
using InteractiveUtils
using Test

include("utils.jl")
include("models.jl")
include("train.jl")

max_stars_repo_name: matsueushi/SinGAN

max_issues_repo_name: matsueushi/SinGAN

max_forks_repo_name: matsueushi/SinGAN


text:
export Scenario, BlackBox, BipartiteSignaling, LocalSignaling

export BipartiteNonSignaling

"""
A `Scenario` is an abstract type parent to all Bell scenarios.
Each child of this abstract type describes a distinct black-box system configuration.
"""
abstract type Scenario end

"""
    BlackBox(num_out :: Int64, num_in :: Int64) <: Scenario

A Bell scenario consisting of a single black-box device with `num_in` inputs and
`num_out` outputs.

![Black-Box Device](../assets/scenario_images/blackbox_device.png)

The black-box device computes the output ``y`` from the input ``x`` by performing
a stochastic map ``S(y|x)``.

### Errors
A `DomainError` is thrown if parameters `num_out` or `num_in` is less than 1.
"""
struct BlackBox <: Scenario
    num_out :: Int64
    num_in :: Int64
    BlackBox(
        num_out::Int64,
        num_in::Int64
    ) = ((num_out >= 1) && (num_in >= 1)) ? new(num_out, num_in) : throw(
        DomainError((num_out, num_in), "inputs must be ≥ 1")
    )
end

"""
    LocalSignaling(
        X :: Int64,
        Y :: Int64,
        d :: Int64,
    ) <: Scenario

A bipartite signaling scenario where information is passed from a transmitter
black-box to a receiver black-box using no more than `d` dits of communication.

![Local Signaling Scenario](../assets/scenario_images/local_signaling_scenario.png)

The transmitter device has `X` inputs and the receiver device has `Y` outputs and
shared randomness is held between the two devices.
When quantum communication is used instead of classical communication no Bell violations
occur.

### Errors
A `DomainError` is thrown if `X`, `Y`, or `d` is less than 1.
"""
struct LocalSignaling <: Scenario
    X :: Int64
    Y :: Int64
    d :: Int64
    LocalSignaling(X::Int64, Y::Int64, d::Int64) = (X ≥ 1 && Y ≥ 1 && d ≥ 1) ? new(X,Y,d) : throw(
        DomainError((X, Y, d), "parameters must be ≥ 1")
    )
end

"""
    BipartiteNonSignaling(
        A :: Int64,
        B :: Int64,
        X :: Int64,
        Y :: Int64
    ) <: Scenario

A bipartite non-signaling scenario where each device receives an input and produces
an output.
Let Alice be the device with `A` outputs and `X` inputs while Bob is the device
with `B` outputs and `Y` inputs.

![Bipartite Non-Signaling Scenario](../assets/scenario_images/bipartite_non_signaling_scenario.png)

Shared randomness is held between Alice and Bob.
When Alice and Bob share quantum entanglement, Bell violations are known to occur.

### Errors
A `DomainError` is thrown if `A`, `B`, `X`, or `Y` is less than 1.
"""
struct BipartiteNonSignaling <: Scenario
    A :: Int64
    B :: Int64
    X :: Int64
    Y :: Int64
    BipartiteNonSignaling(
        A::Int64, B::Int64, X::Int64, Y::Int64
    ) = (
            (A >= 1) && (B >= 1) && (X >= 1) && (Y >= 1)
        ) ? new(A, B, X, Y) : throw(
            DomainError((A, B, X, Y), "inputs must be ≥ 1")
        )
end

"""
    BipartiteSignaling(
        A :: Tuple{Int64, Int64},
        B :: Tuple{Int64, Int64};
        dits :: Int64 = 1,
        bidirectional :: Bool = false
    ) <: Scenario

A bipartite signaling scenario where each device can send a message to the other.

![Bipartite Signaling Scenario](../assets/scenario_images/bipartite_signaling_scenario.png)
"""
struct BipartiteSignaling <: Scenario
    A :: BlackBox
    B :: BlackBox
    dits :: Int64
    bidirectional :: Bool
    BipartiteSignaling(
        A::Tuple{Int64,Int64}, B::Tuple{Int64,Int64};
        dits::Int64=1, bidirectional::Bool=false
    ) = (dits >= 1) ? new(BlackBox(A...), BlackBox(B...), dits, bidirectional) : throw(
            DomainError(dits, "communication `dits` must be ≥ 1.")
        )
end


text:
#@variable(m, x )              # No bounds
#@variable(m, x >= lb )        # Lower bound only
#@variable(m, x <= ub )        # Upper bound only
#@variable(m, x == fixedval )  # Fixed to a value (lb == ub)

function calVarDict(varKey, label, nowLength)
  global varLength = varLength + nowLength
  varStartValue[varKey] = varStartLength
  varEndValue[varKey] = varLength
  varGroup[varKey] = label
  global varStartLength = varStartLength + nowLength
end

macro solveDEA_varibleInit(model)
  global varLength = 0
  global varStartLength = 1
  global varibleInitFlag = true
  global varStartValue = Dict()
  global varEndValue = Dict()
  global varGroup = Dict()
  global varSolution = Dict()
end

macro solveDEA_varible(m, xInequality, label)
  if xInequality.args[2] == :(>=) # Lower bound only
    tmpStr = string(xInequality.args[1])
    if contains(tmpStr, ":")
      start = search(tmpStr, '[')
      mid = search(tmpStr, ':')
      End = search(tmpStr, ']')
      scale = parse(Int32, tmpStr[mid+1:End-1])
      varStr = tmpStr[1:start-1]
      calVarDict(string(varStr), string(label), scale)
      quote
        @variable($m, $xInequality)
      end
    else
      calVarDict(string(tmpStr), string(label), 1)
      quote
        @variable($m, $xInequality)
      end
    end
  elseif xInequality.args[2] == :(<=) # Lower bound only or Lower and upper bounds
    tmpStr = string(xInequality.args[1])
    if contains(tmpStr, ":")
      start = search(tmpStr, '[')
      mid = search(tmpStr, ':')
      End = search(tmpStr, ']')
      scale = parse(Int32, tmpStr[mid+1:End-1])
      varStr = tmpStr[1:start-1]
      calVarDict(string(varStr), string(label), scale)
      quote
        @variable($m, $xInequality)
      end
    else
      calVarDict(string(tmpStr), string(label), 1)
      quote
        @variable($m, $xInequality)
      end
    end
  elseif xInequality.args[2] == :(==) # No bounds
    tmpStr = string(xInequality.args[1])
    if contains(tmpStr, ":")
      start = search(tmpStr, '[')
      mid = search(tmpStr, ':')
      End = search(tmpStr, ']')
      scale = parse(Int32, tmpStr[mid+1:End-1])
      varStr = tmpStr[1:start-1]
      calVarDict(string(varStr), string(label), scale)
      quote
        @variable($m, $xInequality)
      end
    else
      calVarLength(1)
      calVarDict(string(tmpStr), string(label), 1)
      quote
        @variable($m, $xInequality)
      end
    end
  else # No bounds
    tmpStr = string(xInequality.args[1])
    if contains(tmpStr, ":")
      start = search(tmpStr, '[')
      mid = search(tmpStr, ':')
      End = search(tmpStr, ']')
      scale = parse(Int32, tmpStr[mid+1:End-1])
      varStr = tmpStr[1:start-1]
      calVarDict(string(varStr), string(label), scale)
      quote
        @variable($m, $xInequality)
      end
    else
      calVarLength(1)
      calVarDict(string(tmpStr), string(label), 1)
      quote
        @variable($m, $xInequality)
      end
    end
  end
end


text:
module Test_dlyap

using LinearAlgebra
using MatrixEquations
using Test

@testset "Testing discrete Lyapunov equation solvers" begin

n = 100
Ty = Float64

@testset "Discrete Lyapunov equations" begin

reltol = eps(float(100))
a = -2+im; b = 2; @time x = lyapd(a,b)
@test abs(a*x*a'-x+b) < reltol

reltol = eps(float(100f0))
a = 1f0-2f0im; b = 2f0; @time x = lyapd(a,b)
@test abs(a*x*a'-x+b)  < reltol

for Ty in (Float64, Float32)

ar = rand(Ty,n,n)
ac = rand(Ty,n,n)+im*rand(Ty,n,n)
c = rand(Ty,n,n)+im*rand(Ty,n,n)
qc = c'*c
Qr = real(qc)
Ty == Float64 ? reltol = eps(float(100)) : reltol = eps(100*n*one(Ty))

@time x = lyapd(ac,qc);
@test norm(ac*x*ac'-x+qc)/norm(x)/max(1.,norm(ac)^2) < reltol

α = 3+im; # α = 1; #SingularException
@time x = lyapd(α*I,qc);
@test norm(α*x*α'-x+qc)/norm(x) < reltol

try
   α = 1; #SingularException
   @time x = lyapd(α,qc);
catch
   @test true
end

α = 3. +im;
@time x = lyapd(α,qc);
@test norm(α*x*α'-x+qc)/norm(x) < reltol

@time x = lyapd(ac',qc);
@test norm(ac'*x*ac-x+qc)/norm(x)/max(1.,norm(ac)^2) < reltol

@time x = lyapd(ar,Qr)
@test norm(ar*x*ar'-x+Qr)/norm(x)/max(1.,norm(ar)^2) < reltol

@time x = lyapd(ar',Qr)
@test norm(ar'*x*ar-x+Qr)/norm(x)/max(1.,norm(ar)^2) < reltol

@time x = lyapd(ac,Qr);
@test norm(ac*x*ac'-x+Qr)/norm(x)/max(1.,norm(ac)^2)  < reltol

@time x = lyapd(ac',Qr);
@test norm(ac'*x*ac-x+Qr)/norm(x)/max(1.,norm(ac)^2)  < reltol

@time x = lyapd(ar,qc)
@test norm(ar*x*ar'-x+qc)/norm(x)/max(1.,norm(ar)^2) < reltol

@time x = lyapd(ar',qc)
@test norm(ar'*x*ar-x+qc)/norm(x)/max(1.,norm(ar)^2) < reltol
end
end

@testset "Discrete generalized Lyapunov equations" begin

reltol = eps(float(100))
ac = -2+im; ec = 4+im; b = 2; @time x = lyapd(ac,ec,b)
@test abs(ac*x*ac'-ec*x*ec'+b) < reltol

for Ty in (Float64, Float32)

ar = rand(Ty,n,n)
ac = rand(Ty,n,n)+im*rand(Ty,n,n)
er = rand(Ty,n,n)
ec = er+im*rand(Ty,n,n)

c = rand(Ty,n,n)+im*rand(Ty,n,n)
qc = c'*c
Qr = real(qc)
Ty == Float64 ? reltol = eps(float(100)) : reltol = eps(100*n*one(Ty))

@time x = lyapd(ac,ec,qc);
@test norm(ac*x*ac'-ec*x*ec'+qc)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

β = 3
@time x = lyapd(ac,β*I,qc);
@test norm(ac*x*ac'-β*x*β'+qc)/norm(x)/norm(ac) < reltol

β = (1+im);
@time x = lyapd(ac,β,qc);
@test norm(ac*x*ac'-β*x*β'+qc)/norm(x)/norm(ac) < reltol

α = 2+3im
@time x = lyapd(α*I,ec,qc);
@test norm(α*x*α'-ec*x*ec'+qc)/norm(x)/norm(ec) < reltol

α = 2+3im
@time x = lyapd(α,ec,qc);
@test norm(α*x*α'-ec*x*ec'+qc)/norm(x)/norm(ec) < reltol

try
   @time x = lyapd(0*ac,0*ec,qc);
   @test false
catch
   @test true
end

@test norm(α*x*α'-ec*x*ec'+qc)/norm(x)/norm(ec) < reltol

α = 2+3im; β = (1+im);
@time x = lyapd(α*I,β*I,qc);
@test norm(α*x*α'-β*x*β'+qc)/norm(x)/norm(ac)/norm(ec) < reltol

try
   α = 1; β = 1; #SingularException
   @time x = lyapd(α,β,qc);
   @test false
catch
   @test true
end

α = 2+3im; β = (1+im);
@time x = lyapd(α,β,qc);
@test norm(α*x*α'-β*x*β'+qc)/norm(x)/norm(ac)/norm(ec) < reltol

@time x = lyapd(ac',ec',qc);
@test norm(ac'*x*ac-ec'*x*ec+qc)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

β = (1+im);
@time x = lyapd(ac',β*I,qc);
@test norm(ac'*x*ac-β*x*β'+qc)/norm(x)/norm(ac) < reltol

α = 2+3im
@time x = lyapd(α*I,ec',qc);
@test norm(α*x*α'-ec'*x*ec+qc)/norm(x)/norm(ec) < reltol

@time x = lyapd(ac',ec,qc);
@test norm(ac'*x*ac-ec*x*ec'+qc)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

@time x = lyapd(ac,ec',qc);
@test norm(ac*x*ac'-ec'*x*ec+qc)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

@time x = lyapd(ar,er,Qr);
@test norm(ar*x*ar'-er*x*er'+Qr)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ar',er',Qr);
@test norm(ar'*x*ar-er'*x*er+Qr)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ar',er,Qr);
@test norm(ar'*x*ar-er*x*er'+Qr)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ar,er',Qr);
@test norm(ar*x*ar'-er'*x*er+Qr)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ac,ec,Qr);
@test norm(ac*x*ac'-ec*x*ec'+Qr)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

@time x = lyapd(ac',ec',Qr);
@test norm(ac'*x*ac-ec'*x*ec+Qr)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

@time x = lyapd(ac',ec,Qr);
@test norm(ac'*x*ac-ec*x*ec'+Qr)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

@time x = lyapd(ac,ec',Qr);
@test norm(ac*x*ac'-ec'*x*ec+Qr)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ar,er,qc);
@test norm(ar*x*ar'-er*x*er'+qc)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ar',er',qc);
@test norm(ar'*x*ar-er'*x*er+qc)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ar',er,qc);
@test norm(ar'*x*ar-er*x*er'+qc)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ar,er',qc);
@test norm(ar*x*ar'-er'*x*er+qc)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol
end
end

@testset "Discrete Lyapunov equations - Schur form" begin

for Ty in (Float64, Float32)

ar = rand(Ty,n,n)
ac = rand(Ty,n,n)+im*rand(Ty,n,n)
er = rand(Ty,n,n)
ec = er+im*rand(Ty,n,n)
as, es = schur(ar,er)
acs, ecs = schur(ac,ec)

c = rand(Ty,n,n)+im*rand(Ty,n,n)
qc = c'*c
Qr = real(qc)
Ty == Float64 ? reltol = eps(float(100)) : reltol = eps(100*n*one(Ty))

x = copy(Qr)
@time lyapds!(as,x);
@test norm(as*x*as'+Qr-x)/norm(x)/max(norm(as)^2,1.) < reltol

x = copy(qc)
@time lyapds!(as,x);
@test norm(as*x*as'+qc-x)/norm(x)/max(norm(as)^2,1.) < reltol

x = copy(Qr)
@time lyapds!(as,x,adj=true);
@test norm(as'*x*as+Qr-x)/norm(x)/max(norm(as)^2,1.) < reltol

x = copy(qc)
@time lyapds!(as,x,adj=true);
@test norm(as'*x*as+qc-x)/norm(x)/max(norm(as)^2,1.) < reltol

x = copy(qc)
@time lyapds!(acs,ecs,x);
@test norm(acs*x*acs'+qc-ecs*x*ecs')/norm(x)/max(norm(acs)^2,norm(ecs)^2) < reltol

x = copy(qc)
@time lyapds!(acs,I,x);
@test norm(acs*x*acs'+qc-x)/norm(x)/max(norm(acs)^2,norm(ecs)^2) < reltol

x = copy(qc)
@time lyapds!(acs,ecs,x,adj=true);
@test norm(acs'*x*acs+qc-ecs'*x*ecs)/norm(x)/max(norm(acs)^2,norm(ecs)^2) < reltol

x = copy(Qr)
@time lyapds!(as,es,x);
@test norm(as*x*as'+Qr-es*x*es')/norm(x)/max(norm(as)^2,norm(es)^2) < reltol

x = copy(Qr)
@time lyapds!(as,es,x,adj=true);
@test norm(as'*x*as+Qr-es'*x*es)/norm(x)/max(norm(as)^2,norm(es)^2) < reltol

x = copy(qc)
@time lyapds!(as,es,x);
@test norm(as*x*as'+qc-es*x*es')/norm(x)/max(norm(as)^2,norm(es)^2) < reltol

x = copy(qc)
@time lyapds!(as,es,x,adj=true);
@test norm(as'*x*as+qc-es'*x*es)/norm(x)/max(norm(as)^2,norm(es)^2) < reltol

end
end

end

end


text:
using SymmetricProducts
using Test

using LinearAlgebra: Diagonal, Hermitian, Symmetric

function test_symprod(A, ::Type{T}) where T
    A′A = SELF'*A
    @test A′A == A'*A
    @test A′A isa T
end

reals(dims...) = collect(reshape(range(1.0; length = prod(dims)), dims))

cplxs(dims...) =
    collect(reshape([Complex(Float64(i), i^2) for i in 1:prod(dims)], dims))

test_symprod(reals(3,3), Symmetric{Float64})
test_symprod(reals(3,3)', Symmetric{Float64})
test_symprod(reals(3)', Symmetric{Float64})
test_symprod(Diagonal(reals(3)), Diagonal{Float64})
test_symprod(Diagonal(reals(3))', Diagonal{Float64})
test_symprod(reals(3), Float64)

test_symprod(cplxs(3,3), Hermitian{Complex{Float64}})
test_symprod(cplxs(3,3)', Hermitian{Complex{Float64}})
test_symprod(cplxs(3)', Hermitian{Complex{Float64}})
test_symprod(Diagonal(cplxs(3)), Diagonal{Float64})
test_symprod(Diagonal(cplxs(3))', Diagonal{Float64})
test_symprod(cplxs(3), Float64)

max_stars_repo_name: tpapp/SymmetricProducts

max_issues_repo_name: tpapp/SymmetricProducts

max_forks_repo_name: tpapp/SymmetricProducts


text:
export ExpConstraint, conic_form!, vexity

### (Primal) exponential cone constraint ExpConstraint(x,y,z) => y exp(x/y) <= z & y>=0
struct ExpConstraint <: Constraint
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr, AbstractExpr, AbstractExpr} # (x, y, z)
    size::Tuple{Int, Int}
    dual::ValueOrNothing

    function ExpConstraint(x::AbstractExpr, y::AbstractExpr, z::AbstractExpr)
        @assert(x.size == y.size == z.size,
                "Exponential constraint requires x, y, and z to be of same size")
        # @assert(x.size == (1,1),
        #         "Exponential constraint requires x, y, and z to be scalar for now")
        sz = x.size
        id_hash = hash((x,y,z, :exp))
        return new(:exp, id_hash, (x, y, z), sz, nothing)
    end
end

ExpConstraint(x::AbstractExpr, y, z::AbstractExpr) = ExpConstraint(x, Constant(y), z)
ExpConstraint(x::AbstractExpr, y::AbstractExpr, z) = ExpConstraint(x, y, Constant(z))
ExpConstraint(x, y::AbstractExpr, z::AbstractExpr) = ExpConstraint(Constant(x), y, z)

function vexity(c::ExpConstraint)
    # TODO: check these...
    if vexity(c.x) == ConcaveVexity()
        error("Exponential constraint requires x to be convex")
    end
    if vexity(c.y)!=ConstVexity()
        error("Exponential constraint requires y to be constant")
    end
    if vexity(c.z) == ConvexVexity()
        error("Exponential constraint requires z to be concave")
    end
    return ConvexVexity()
end

function conic_form!(c::ExpConstraint, unique_conic_forms::UniqueConicForms=UniqueConicForms())
    if !has_conic_form(unique_conic_forms, c)
        conic_constrs = ConicConstr[]
        if c.size == (1, 1)
            objectives = Vector{ConicObj}(undef, 3)
            @inbounds for iobj = 1:3
                objectives[iobj] = conic_form!(c.children[iobj], unique_conic_forms)
            end
            push!(conic_constrs, ConicConstr(objectives, :ExpPrimal, [1, 1, 1]))
        else
            for i = 1:c.size[1]
                for j = 1:c.size[2]
                    objectives = Vector{ConicObj}(undef, 3)
                    @inbounds for iobj = 1:3
                        objectives[iobj] = conic_form!(c.children[iobj][i,j], unique_conic_forms)
                    end
                    push!(conic_constrs, ConicConstr(objectives, :ExpPrimal, [1, 1, 1]))
                end
            end
        end
        cache_conic_form!(unique_conic_forms, c, conic_constrs)
    end
    return get_conic_form(unique_conic_forms, c)
end


text:
using Gadfly

n = 10
xs, ys, labels = rand(n), rand(n), [randstring(6) for _ in 1:n]

l = plot(x=xs, y=ys, label=labels, Geom.point, Geom.label(position=:left))
r = plot(x=xs, y=ys, label=labels, Geom.point, Geom.label(position=:right))
a = plot(x=xs, y=ys, label=labels, Geom.point, Geom.label(position=:above))
b = plot(x=xs, y=ys, label=labels, Geom.point, Geom.label(position=:below))
c = plot(x=xs, y=ys, label=labels, Geom.point, Geom.label(position=:centered))

vstack(l, r, a, b, c)


text:
function *(n::Int, asset::FinancialAsset)
        asset.quantity = asset.quantity * n
    asset
end

function +(x::LongPut, y::ShortPut)
    if x.ticker != y.ticker
        error("Only same underlying spreads supported")
    elseif quantity != y.quantity
        error("Only same contract size supported")
    elseif x.expiry != y.expiry
        error("Only same expiry supported")
    elseif x.strike == y.strike
        error("Options cancel each other out")
    elseif x.strike > y.strike
        ShortPutVertical(x,y)
    else
        LongPutVertical(x,y)
    end
end

function +(x::Stock, y::Stock)
    if x.ticker != y.ticker
        error("Only same stock addition supported")
    elseif x.currency != y.currency
        error("Only same currency stock addition supported")
    else
        nothing
    end

    new_basis  = x.basis * x.quantity/(x.quantity + y.quantity) + y.basis * y.quantity/(x.quantity + y.quantity)
    new_quantity = x.quantity + y.quantity
    new_id     = isnull(x.id) ? y.id : x.id

    Stock(x.ticker, new_basis, new_quantity , x.currency, .01, 1., new_id)
end


text:
@testset "SlidingBuffer" begin

    @testset "Core Functionality" begin
        buf = SlidingBuffer{Int}(5)
        @testset "When empty" begin
            @test length(buf) == 0
            @test capacity(buf) == 5
            @test_throws BoundsError first(buf)
            @test isempty(buf) == true
            @test isfull(buf) == false
            @test eltype(buf) == Int
            @test eltype(typeof(buf)) == Int
        end

        @testset "With 1 element" begin
            push!(buf, 1)
            @test length(buf) == 1
            @test capacity(buf) == 5
            @test isfull(buf) == false
        end

        @testset "Appending many elements" begin
            # append!(buf, 2:8) for offset buffers you cannot in one operation append more elements than buffer length.
            append!(buf, 2:5)
            append!(buf, 6:8)
            @test length(buf) == capacity(buf)
            @test size(buf) == (length(buf),)
            @test isempty(buf) == false
            @test isfull(buf) == true
            @test convert(Array, buf) == Int[4,5,6,7,8]
        end

        @testset "getindex" begin
            @test buf[4] == 4
            @test buf[5] == 5
            @test buf[6] == 6
            @test buf[7] == 7
            @test buf[8] == 8
            @test_throws BoundsError buf[3]
            @test_throws BoundsError buf[3:6]
            @test_throws BoundsError buf[9]
            @test_throws BoundsError buf[6:9]
            @test buf[6:7] == Int[6,7]
            @test buf[[4,8]] == Int[4,8]
        end

        @testset "setindex" begin
            buf[6] = 999
            @test convert(Array, buf) == Int[4,5,999,7,8]
        end
    end

    @testset "other constructor" begin
        buf = SlidingBuffer(10)
        @test length(buf) == 0
        @test typeof(buf) <: SlidingBuffer{Any}
    end

    @testset "pushfirst" begin
        buf = SlidingBuffer{Int}(5)  # New, empty one for full test coverage
        for i in -5:5
            pushfirst!(buf, i)
        end
        arr = convert(Array, buf)
        @test arr == Int[5, 4, 3, 2, 1]
        for (idx, n) in enumerate(5:1)
            @test arr[idx] == n
        end
    end

    @testset "Issue 429" begin
        buf = SlidingBuffer{Int}(5)
        map(x -> pushfirst!(buf, x), 1:8)
        pop!(buf)
        pushfirst!(buf, 9)
        @test length(buf.buffer) == buf.capacity
        arr = convert(Array, buf)
        @test arr == Int[9, 8, 7, 6, 5]
    end

    @testset "Issue 379" begin
        buf = SlidingBuffer{Int}(5)
        pushfirst!(buf, 1)
        @test convert(Array, buf) == [1]
        pushfirst!(buf, 2)
        @test convert(Array, buf) == [2, 1]
    end

    @testset "empty!" begin
        buf = SlidingBuffer{Int}(5)
        push!(buf, 13)
        empty!(buf)
        @test length(buf) == 0
    end

    @testset "pop!" begin
        buf = SlidingBuffer{Int}(5)
        for i in 0:5    # one extra to force wraparound
            push!(buf, i)
        end
        for j in 5:-1:1
            @test pop!(buf) == j
            @test convert(Array, buf) == collect(1:j-1)
        end
        @test isempty(buf)
        @test_throws ArgumentError pop!(buf)
    end

    @testset "popfirst!" begin
        buf = SlidingBuffer{Int}(5)
        for i in 0:5    # one extra to force wraparound
            push!(buf, i)
        end
        for j in 1:5
            @test popfirst!(buf) == j
            @test convert(Array, buf) == collect(j+1:5)
        end
        @test isempty(buf)
        @test_throws ArgumentError popfirst!(buf)
    end

    @testset "fill!" begin
        @testset "fill an empty buffer" begin
            buf = SlidingBuffer{Int}(3)
            fill!(buf, 42)
            @test Array(buf) == [42, 42, 42]
        end
        @testset "fill a non empty buffer" begin
            buf = SlidingBuffer{Int}(3)
            push!(buf, 21)
            fill!(buf, 42)
            @test Array(buf) == [21, 42, 42]
        end
    end

    @testset "Growing indexing" begin
        @testset "push / pop with growing indexes" begin
            buf = SlidingBuffer{Int}(5)
            for i=1:8
                push!(buf, i)
            end

            @test firstindex(buf) == first(buf) == 4
            @test lastindex(buf) == last(buf) == 8

            @test Array(buf) == [4,5,6,7,8]

            for i=4:5
                popfirst!(buf)
            end

            @test firstindex(buf) == first(buf) == 6
            @test lastindex(buf) == last(buf) == 8

            @test length(buf) == 3
            @test Array(buf) == [6,7,8]
        end
        @testset "append with growing indexes" begin
            buf = SlidingBuffer{Int}(10)
            append!(buf, 1:10)

            append!(buf, 11:17)

            @test length(buf) == 10

            @test firstindex(buf) == first(buf) == 8
            @test lastindex(buf) == last(buf) == 17

            @test Array(buf) == 8:17

            for i = 8:12
                popfirst!(buf)
            end

            @test length(buf) == 5

            @test firstindex(buf) == first(buf) == 13
            @test lastindex(buf) == last(buf) == 17

            @test Array(buf) == 13:17
        end
    end

end


text:
# Autogenerated wrapper script for LibModbus_jll for x86_64-apple-darwin
export libmodbus

JLLWrappers.@generate_wrapper_header("LibModbus")
JLLWrappers.@declare_library_product(libmodbus, "@rpath/libmodbus.5.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libmodbus,
        "lib/libmodbus.5.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Pkg
Pkg.activate(".")
Pkg.add("Gadfly")
Pkg.add("Plots")
Pkg.update()

println("Making Figure 1")
include("figures/figureJ1.jl")
figureJ1()

println("Making Figure 2")
include("figures/figureJ2.jl")
figureJ2()


text:
module ParSpMatVec
using SparseArrays
using LinearAlgebra
using Libdl
const spmatveclib  = abspath(joinpath(splitdir(Base.source_path())[1],"..","deps","builds","ParSpMatVec"))

include("A_mul_B.jl")
include("Ac_mul_B.jl")

export isBuilt
function isBuilt()
	println(spmatveclib)
	return find_library([spmatveclib])!=""
end

end # module


text:
@testset "653.two-sum-iv-input-is-a-bst.jl" begin
    @test find_target_inorder_traversal(TreeNode{Int}([5, 3, 6, 2, 4, nothing, 7]), 9) ==
        true
    @test find_target_inorder_traversal(TreeNode{Int}([5, 3, 6, 2, 4, nothing, 7]), 28) ==
        false
    @test find_target_inorder_traversal(TreeNode{Int}([2, 1, 3]), 4) == true
    @test find_target_inorder_traversal(TreeNode{Int}([2, 1, 3]), 1) == false
    @test find_target_inorder_traversal(TreeNode{Int}([2, 1, 3]), 3) == true
end


text:
export SimulatedObservation!, ElasticPropagatorSolver, AcousticPropagatorSolver

function gather_ops(ipt,ii)
    gather_ops_ = load_op_and_grad("$(@__DIR__)/../deps/CustomOps/build/libADSeismic","gather_ops")
    ipt,ii = convert_to_tensor(Any[ipt,ii], [Float64,Int64])
    out = gather_ops_(ipt,ii)
    set_shape(out, (length(ii), ))
end

function Base.:getindex(o::PyObject, ii::Array{Int64,1})
    gather_ops(o, ii)
end

function scatter_add_op(ipt, ii, vv)
    scatter_add_ops_ = load_op_and_grad("$(@__DIR__)/../deps/CustomOps/build/libADSeismic","scatter_add_ops")
    ipt,ii,vv = convert_to_tensor(Any[ipt,ii,vv], [Float64,Int64,Float64])
    out = scatter_add_ops_(ipt,ii,vv)
    set_shape(out, (length(ipt), ))
    # ipt = scatter_add(ipt, ii, vv)
end

function scatter_nd_ops(ii,vv,m)
    scatter_nd_ops_ = load_op_and_grad("$(@__DIR__)/../deps/CustomOps/build/libADSeismic","scatter_nd_ops")
    ii,vv,m_ = convert_to_tensor(Any[ii,vv,m], [Int64,Float64,Int64])
    out = scatter_nd_ops_(ii,vv,m_)
    set_shape(out, (m,))
    # vector(ii, vv, m)
end

function ElasticPropagatorSolver(param::ElasticPropagatorParams, src::ElasticSource,
    ρ::Union{PyObject, Array{Float64, 2}}, λ::Union{PyObject, Array{Float64, 2}}, μ::Union{PyObject, Array{Float64, 2}})
    ρ = tf.reshape(convert_to_tensor(ρ),(-1,))
    λ = tf.reshape(convert_to_tensor(λ),(-1,))
    μ = tf.reshape(convert_to_tensor(μ),(-1,))
    NSTEP = param.NSTEP
    vx_arr = TensorArray(NSTEP+1); vx_arr = write(vx_arr, 1, constant(zeros((param.NX+2)*(param.NY+2))))
    vy_arr = TensorArray(NSTEP+1); vy_arr = write(vy_arr, 1, constant(zeros((param.NX+2)*(param.NY+2))))
    sigmaxx_arr = TensorArray(NSTEP+1); sigmaxx_arr = write(sigmaxx_arr, 1, constant(zeros((param.NX+2)*(param.NY+2))))
    sigmayy_arr = TensorArray(NSTEP+1); sigmayy_arr = write(sigmayy_arr, 1, constant(zeros((param.NX+2)*(param.NY+2))))
    sigmaxy_arr = TensorArray(NSTEP+1); sigmaxy_arr = write(sigmaxy_arr, 1, constant(zeros((param.NX+2)*(param.NY+2))))
    mem_arr = [TensorArray(NSTEP+1) for i = 1:8]
    for i = 1:8
        mem_arr[i] = write(mem_arr[i], 1, constant(zeros((param.NX+2)*(param.NY+2))))
    end
    ax,bx,kx,ay,by,ky = compute_PML_Params(param)

    function condition(i, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arr...)
        i<=NSTEP+1
    end

    function body(i, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arr...)
        mem = Array{PyObject}(undef, 8)
        σxx,σyy,σxy,vx,vy =
            read(sigmaxx_arr, i-1), read(sigmayy_arr, i-1), read(sigmaxy_arr, i-1),
            read(vx_arr, i-1), read(vy_arr, i-1)
        for k = 1:8
            mem[k] = read(mem_arr[k], i-1)
        end
        srci,srcj,srctype,srcv = ElasticSourceAtTimeT(src, i-1)
        σxx, σyy, σxy,vx, vy, mem = one_step(param, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
            ax, bx, ay, by, kx, ky, srci, srcj, srcv, srctype)

        sigmaxx_arr,sigmayy_arr,sigmaxy_arr,vx_arr,vy_arr =
            write(sigmaxx_arr, i, σxx), write(sigmayy_arr, i, σyy), write(sigmaxy_arr, i, σxy),
                    write(vx_arr, i, vx), write(vy_arr, i, vy)

        mem_arr_ = Array{PyObject}(undef, 8)
        for k = 1:8
            mem_arr_[k] = write(mem_arr[k], i, mem[k])
        end

        if param.IT_DISPLAY>0
            op = tf.cond(tf.equal(tf.math.floormod(i,param.IT_DISPLAY),1), ()->tf.print("time step = ", i-1, " out of $NSTEP"), ()->tf.no_op())
            i = bind(i, op)
        end

        i+1, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arr_...
    end
    i = constant(2, dtype=Int32)

    _, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arrs =
        while_loop(condition, body, [i, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arr...]; parallel_iterations = 1)

    vx = tf.reshape(stack(vx_arr), (param.NSTEP+1, param.NX+2, param.NY+2))
    vy = tf.reshape(stack(vy_arr), (param.NSTEP+1, param.NX+2, param.NY+2))
    σxx = tf.reshape(stack(sigmaxx_arr), (param.NSTEP+1, param.NX+2, param.NY+2))
    σyy = tf.reshape(stack(sigmayy_arr), (param.NSTEP+1, param.NX+2, param.NY+2))
    σxy = tf.reshape(stack(sigmaxy_arr), (param.NSTEP+1, param.NX+2, param.NY+2))
    ElasticPropagator(
        param, src, vx, vy, σxx, σyy, σxy, mem_arrs
    )
end

function fw1(param::ElasticPropagatorParams, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
        ax, bx, ay, by, kx, ky)
    NX, NY = param.NX, param.NY
    DELTAX, DELTAY, DELTAT = param.DELTAX, param.DELTAY, param.DELTAT
    i_1j = getid(param, 1:NX-1, 3:NY+1)
    ij = getid(param, 2:NX, 3:NY+1)
    i1j = getid(param, 3:NX+1, 3:NY+1)
    i2j = getid(param, 4:NX+2, 3:NY+1)

    ij_2 = getid(param, 2:NX, 1:NY-1)
    ij_1 = getid(param, 2:NX, 2:NY)
    ij1 = getid(param, 2:NX, 4:NY+2)

    λ_ = 0.5*(λ[i1j]+λ[ij])
    μ_ = 0.5*(μ[i1j]+μ[ij])
    λμ = λ_ + 2μ_
    value_dvx_dx = (27vx[i1j]-27vx[ij]-vx[i2j]+vx[i_1j])/(24*DELTAX)
    value_dvy_dy = (27*vy[ij]-27*vy[ij_1]-vy[ij1]+vy[ij_2]) / (24*DELTAY)

    # broadcast
    mem[1] =  makevector(mem[1], ij, adbroadcast(bx[2,1:end-1], mem[1][ij], 1) + adbroadcast(ax[2,1:end-1], value_dvx_dx, 1));
    mem[2] =  makevector(mem[2], ij, adbroadcast(by[1,2:end], mem[2][ij], 2) + adbroadcast(ay[1,2:end], value_dvy_dy, 2));

    value_dvx_dx = adbroadcast(value_dvx_dx, kx[2,1:end-1], 3) + mem[1][ij]
    value_dvy_dy = adbroadcast(value_dvy_dy, ky[1,2:end], 4) + mem[2][ij]

    σxx = scatter_add_op(σxx, ij, (λμ * value_dvx_dx + λ_ * value_dvy_dy) * DELTAT )
    σyy = scatter_add_op(σyy, ij, (λμ * value_dvy_dy + λ_ * value_dvx_dx) * DELTAT )

    return vx, vy, σxx, σyy,σxy, mem
end

function fw2(param::ElasticPropagatorParams, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
    ax, bx, ay, by, kx, ky)
    NX, NY = param.NX, param.NY
    DELTAX, DELTAY, DELTAT = param.DELTAX, param.DELTAY, param.DELTAT
    i_2j = getid(param, 1:NX-1, 2:NY)
    i_1j = getid(param, 2:NX, 2:NY)
    ij = getid(param, 3:NX+1, 2:NY)
    i1j = getid(param, 4:NX+2, 2:NY)

    ij_1 = getid(param, 3:NX+1, 1:NY-1)
    ij1 = getid(param, 3:NX+1, 3:NY+1)
    ij2 = getid(param, 3:NX+1, 4:NY+2)

    μ_ = 0.5*(μ[ij]+μ[ij1])
    value_dvy_dx = (27vy[ij]-27vy[i_1j]-vy[i1j]+vy[i_2j])/(24*DELTAX)
    value_dvx_dy = (27*vx[ij1]-27*vx[ij]-vx[ij2]+vx[ij_1]) / (24*DELTAY)

    # broadcast
    mem[3] = makevector(mem[3], ij,  adbroadcast(bx[1,2:end], mem[3][ij], 1) + adbroadcast(ax[1, 2:end], value_dvy_dx, 1))
    mem[4] = makevector(mem[4], ij,  adbroadcast(by[2,1:end-1], mem[4][ij], 2) + adbroadcast(ay[2, 1:end-1], value_dvx_dy, 2));

    value_dvy_dx = adbroadcast(value_dvy_dx, kx[1,2:end], 3) + mem[3][ij]
    value_dvx_dy = adbroadcast(value_dvx_dy, ky[1,1:end-1], 4) + mem[4][ij]

    σxy = scatter_add_op(σxy, ij, μ_*(value_dvy_dx + value_dvx_dy) * DELTAT )

    return vx, vy, σxx, σyy, σxy,mem
end

function fw3(param::ElasticPropagatorParams, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
    ax, bx, ay, by, kx, ky)
    NX, NY = param.NX, param.NY
    DELTAX, DELTAY, DELTAT = param.DELTAX, param.DELTAY, param.DELTAT
    i_2j = getid(param, 1:NX-1, 3:NY+1)
    i_1j = getid(param, 2:NX, 3:NY+1)
    ij = getid(param, 3:NX+1, 3:NY+1)
    i1j = getid(param, 4:NX+2, 3:NY+1)

    ij_2 = getid(param, 3:NX+1, 1:NY-1)
    ij_1 = getid(param,3:NX+1, 2:NY)
    ij1 = getid(param,3:NX+1, 4:NY+2)

    value_dsigmaxx_dx = (27*σxx[ij]-27*σxx[i_1j]-σxx[i1j]+σxx[i_2j]) / (24*DELTAX);
    value_dsigmaxy_dy = (27*σxy[ij]-27*σxy[ij_1]-σxy[ij1]+σxy[ij_2]) / (24*DELTAY);

    mem[5] = makevector(mem[5],ij, adbroadcast(bx[1,2:end], mem[5][ij], 1) + adbroadcast(ax[1,2:end], value_dsigmaxx_dx, 1))
    mem[6] = makevector(mem[6],ij, adbroadcast(by[1,2:end], mem[6][ij], 2) + adbroadcast(ay[1,2:end], value_dsigmaxy_dy, 2))

    value_dsigmaxx_dx = adbroadcast(value_dsigmaxx_dx, kx[1,2:end], 3) + mem[5][ij]
    value_dsigmaxy_dy = adbroadcast(value_dsigmaxy_dy, ky[1,2:end], 4) + mem[6][ij]

    vx = scatter_add_op(vx, ij, (value_dsigmaxx_dx + value_dsigmaxy_dy) * DELTAT / ρ[ij])
    return vx, vy, σxx, σyy, σxy, mem
end

function fw4(param::ElasticPropagatorParams, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
    ax, bx, ay, by, kx, ky)
    NX, NY = param.NX, param.NY
    DELTAX, DELTAY, DELTAT = param.DELTAX, param.DELTAY, param.DELTAT
    i_1j = getid(param, 1:NX-1, 2:NY)
    ij = getid(param, 2:NX, 2:NY)
    i1j = getid(param, 3:NX+1, 2:NY)
    i2j = getid(param, 4:NX+2, 2:NY)
    i1j1 = getid(param, 3:NX+1, 3:NY+1)

    ij_1 = getid(param, 2:NX, 1:NY-1)
    ij1 = getid(param, 2:NX, 3:NY+1)
    ij2 = getid(param, 2:NX, 4:NY+2)

    ρ_ = 0.25 * (ρ[ij] + ρ[i1j] + ρ[i1j1] + ρ[ij1]);

    value_dsigmaxy_dx = (27*σxy[i1j]-27*σxy[ij]-σxy[i2j]+σxy[i_1j]) / (24*DELTAX);
    value_dsigmayy_dy = (27*σyy[ij1]-27*σyy[ij]-σyy[ij2]+σyy[ij_1]) / (24*DELTAY);

    mem[7] = makevector(mem[7], ij, adbroadcast(bx[2,1:end-1], mem[7][ij], 1) + adbroadcast(ax[2,1:end-1] , value_dsigmaxy_dx, 1))
    mem[8] = makevector(mem[8], ij, adbroadcast(by[2,1:end-1], mem[8][ij], 2) + adbroadcast(ay[2,1:end-1] , value_dsigmayy_dy, 2))

    value_dsigmaxy_dx = adbroadcast(value_dsigmaxy_dx, kx[2,1:end-1], 3) + mem[7][ij]
    value_dsigmayy_dy = adbroadcast(value_dsigmayy_dy, ky[2,1:end-1], 4) + mem[8][ij]

    vy = scatter_add_op( vy, ij, (value_dsigmaxy_dx + value_dsigmayy_dy) * DELTAT / ρ_)
    return vx, vy, σxx, σyy, σxy,mem
end

function one_step(param::ElasticPropagatorParams, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
    ax, bx, ay, by, kx, ky, srci, srcj, srcv, srctype)
    vx, vy, σxx, σyy, σxy,mem = fw1(param, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
                    ax, bx, ay, by, kx, ky)
    vx, vy, σxx, σyy, σxy,mem = fw2(param, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
                    ax, bx, ay, by, kx, ky)
    vx, vy, σxx, σyy, σxy,mem = fw3(param, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
                    ax, bx, ay, by, kx, ky)
    vx, vy, σxx, σyy, σxy,mem = fw4(param, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
                    ax, bx, ay, by, kx, ky)
    add_source = load_op_and_grad("$(@__DIR__)/../deps/CustomOps/build/libADSeismic","add_source", multiple=true)
    σxx, σyy, σxy,vx, vy  = add_source(σxx, σyy, σxy,vx,vy,srci,srcj,srctype,constant(param.NX),constant(param.NY),srcv)
    return σxx, σyy, σxy,vx, vy, mem
end

function compute_PML_Params(param::ElasticPropagatorParams)
    ## define profile of absorption in PML region
    NPOINTS_PML = param.NPOINTS_PML
    NPOWER = param.NPOWER
    NX = param.NX
    NY = param.NY
    DELTAX = param.DELTAX
    DELTAY = param.DELTAY
    DELTAT = param.DELTAT
    K_MAX_PML = param.K_MAX_PML
    ALPHA_MAX_PML = param.ALPHA_MAX_PML
    USE_PML_XMIN = param.USE_PML_XMIN
    USE_PML_XMAX = param.USE_PML_XMAX
    USE_PML_YMIN = param.USE_PML_YMIN
    USE_PML_YMAX = param.USE_PML_YMAX
    Rcoef = param.Rcoef

    thickness_PML_x = NPOINTS_PML * DELTAX
    thickness_PML_y = NPOINTS_PML * DELTAY

    # reflection coefficient
    d0_x = - (NPOWER + 1) * param.vp_ref * log(Rcoef) / (2. * thickness_PML_x)
    d0_y = - (NPOWER + 1) * param.vp_ref * log(Rcoef) / (2. * thickness_PML_y)

    xoriginleft = thickness_PML_x
    xoriginright = (NX-0.5)*DELTAX - thickness_PML_x

    a_x = zeros(NX)
    a_y = zeros(NY)
    a_x_half = zeros(NX)
    a_y_half = zeros(NY)
    b_x = zeros(NX)
    b_y = zeros(NY)
    b_x_half = zeros(NX)
    b_y_half = zeros(NY)
    d_x = zeros(NX)
    d_y = zeros(NY)
    d_x_half = zeros(NX)
    d_y_half = zeros(NY)
    K_x = ones(NX)
    K_y = ones(NY)
    K_x_half = ones(NX)
    K_y_half = ones(NY)
    alpha_x = zeros(NX)
    alpha_y = zeros(NY)
    alpha_x_half = zeros(NX)
    alpha_y_half = zeros(NY)

    for i = 1:NX
        xval = DELTAX * (i-1)
        if (USE_PML_XMIN)

            abscissa_in_PML = xoriginleft - xval
            if (abscissa_in_PML >= 0.0)
                abscissa_normalized = abscissa_in_PML / thickness_PML_x
                d_x[i] = d0_x * abscissa_normalized^NPOWER
                K_x[i] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_x[i] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_x[i] = ALPHA_MAX_PML * abscissa_normalized
            end

            abscissa_in_PML = xoriginleft - (xval + DELTAX/2.)
            if (abscissa_in_PML >= 0.0)
                abscissa_normalized = abscissa_in_PML / thickness_PML_x
                d_x_half[i] = d0_x * abscissa_normalized^NPOWER
                K_x_half[i] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_x_half[i] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_x_half[i] = ALPHA_MAX_PML * abscissa_normalized
            end

        end

        if (USE_PML_XMAX)
            abscissa_in_PML = xval - xoriginright
            if (abscissa_in_PML >= 0)
                abscissa_normalized = abscissa_in_PML / thickness_PML_x
                d_x[i] = d0_x * abscissa_normalized^NPOWER
                K_x[i] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_x[i] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_x[i] = ALPHA_MAX_PML * abscissa_normalized
            end

            abscissa_in_PML = xval + DELTAX/2.0 - xoriginright
            if (abscissa_in_PML >= 0)
                abscissa_normalized = abscissa_in_PML / thickness_PML_x
                d_x_half[i] = d0_x * abscissa_normalized^NPOWER
                K_x_half[i] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_x_half[i] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_x_half[i] = ALPHA_MAX_PML * abscissa_normalized
            end
        end

        if (alpha_x[i] < 0)
            alpha_x[i] = 0
        end
        if (alpha_x_half[i] < 0)
            alpha_x_half[i] = 0
        end

        b_x[i] = exp(- (d_x[i] / K_x[i] + alpha_x[i]) * DELTAT)
        b_x_half[i] = exp(- (d_x_half[i] / K_x_half[i] + alpha_x_half[i]) * DELTAT)

        if (abs(d_x[i]) > 1e-6)
            a_x[i] = d_x[i] * (b_x[i] - 1.) / (K_x[i] * (d_x[i] + K_x[i] * alpha_x[i]))
        end
        if (abs(d_x_half[i]) > 1e-6)
            a_x_half[i] = d_x_half[i] * (b_x_half[i] - 1.) / (K_x_half[i] * (d_x_half[i] + K_x_half[i] * alpha_x_half[i]))
        end
    end

    yoriginbottom = thickness_PML_y
    yorigintop = (NY-0.5)*DELTAY - thickness_PML_y

    for j = 1:NY

        yval = DELTAY * (j-1)
        if param.USE_PML_YMIN
            abscissa_in_PML = yoriginbottom - yval
            if (abscissa_in_PML >= 0.0)
                abscissa_normalized = abscissa_in_PML / thickness_PML_y
                d_y[j] = d0_y * abscissa_normalized^NPOWER
                K_y[j] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_y[j] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_y[j] = ALPHA_MAX_PML * abscissa_normalized
            end
            abscissa_in_PML = yoriginbottom - (yval + DELTAY/2.)
            if abscissa_in_PML >= 0
                abscissa_normalized = abscissa_in_PML / thickness_PML_y
                d_y_half[j] = d0_y * abscissa_normalized^NPOWER
                K_y_half[j] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_y_half[j] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_y_half[j] = ALPHA_MAX_PML * abscissa_normalized
            end
        end

        if param.USE_PML_YMAX
            abscissa_in_PML = yval - yorigintop
            if abscissa_in_PML >= 0
                abscissa_normalized = abscissa_in_PML / thickness_PML_y
                d_y[j] = d0_y * abscissa_normalized^NPOWER
                K_y[j] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_y[j] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_y[j] = ALPHA_MAX_PML * abscissa_normalized
            end

            abscissa_in_PML = yval + DELTAY/2. - yorigintop
            if abscissa_in_PML >= 0
                abscissa_normalized = abscissa_in_PML / thickness_PML_y
                d_y_half[j] = d0_y * abscissa_normalized^NPOWER
                K_y_half[j] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_y_half[j] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_y_half[j] = ALPHA_MAX_PML * abscissa_normalized
            end
        end

        b_y[j] = exp(- (d_y[j] / K_y[j] + alpha_y[j]) * DELTAT)
        b_y_half[j] = exp(- (d_y_half[j] / K_y_half[j] + alpha_y_half[j]) * DELTAT)
        if abs(d_y[j]) > 1e-6
            a_y[j] = d_y[j] * (b_y[j] - 1.) / (K_y[j] * (d_y[j] + K_y[j] * alpha_y[j]))
        end
        if (abs(d_y_half[j]) > 1e-6)
            a_y_half[j] = d_y_half[j] * (b_y_half[j] - 1.) / (K_y_half[j] * (d_y_half[j] + K_y_half[j] * alpha_y_half[j]))
        end
    end
    # @show norm(d_x-d_y)

    ax = [a_x'; a_x_half']
    bx = [b_x'; b_x_half']
    kx = [K_x'; K_x_half']
    ay = [a_y'; a_y_half']
    by = [b_y'; b_y_half']
    ky = [K_y'; K_y_half']

    return ax,bx,kx,ay,by,ky

end

# # not in use any more
# function ElasticPropagator(src::ElasticSource, λ::Array{Float64},
#     ρ::Array{Float64}, μ::Array{Float64}, param::Union{ElasticPropagatorParams,Missing}=missing)
#     if ismissing(param)
#         param = ElasticPropagatorParams()
#     end
#     NSTEP = param.NSTEP
#     if size(src.srcv,1)<NSTEP
#         error("ADSeismic: length of sources=$(size(src.srci,1)) and NSTEP=$NSTEP does not match")
#     end
#     if size(src.srcv,1)!=NSTEP
#         @warn("ADSeismic: length of sources and NSTEP does not match $(size(src.srci,1))!=$NSTEP")
#     end
#     vx_arr = TensorArray(NSTEP+1); vx_arr = write(vx_arr, 1, zeros(param.NX+2, param.NY+2))
#     vy_arr = TensorArray(NSTEP+1); vy_arr = write(vy_arr, 1, zeros(param.NX+2, param.NY+2))
#     sigmaxx_arr = TensorArray(NSTEP+1); sigmaxx_arr = write(sigmaxx_arr, 1, zeros(param.NX+2, param.NY+2))
#     sigmayy_arr = TensorArray(NSTEP+1); sigmayy_arr = write(sigmayy_arr, 1, zeros(param.NX+2, param.NY+2))
#     sigmaxy_arr = TensorArray(NSTEP+1); sigmaxy_arr = write(sigmaxy_arr, 1, zeros(param.NX+2, param.NY+2))
#     mem_arr = TensorArray(NSTEP+1); mem_arr = write(mem_arr, 1, zeros(8, param.NX+2, param.NY+2))

#     ax,bx,kx,ay,by,ky = compute_PML_Params(param)

#     function condition(i, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arr)
#         i<=NSTEP+1
#     end

#     function body(i, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arr)
#         sigmaxx,sigmayy,sigmaxy,vx,vy,mem =
#             read(sigmaxx_arr, i-1), read(sigmayy_arr, i-1), read(sigmaxy_arr, i-1),
#             read(vx_arr, i-1), read(vy_arr, i-1), read(mem_arr, i-1)
#         srci,srcj,srctype,srcv = ElasticSourceAtTimeT(src, i-1)
#         # @show srci, srcj, srctype, srcv
#         # error()
#         sigmaxx, sigmayy, sigmaxy, vx,vy, mem = elastic_wave_op(λ,μ,ρ,sigmaxx,sigmayy,sigmaxy,vx,vy,
#                 mem,constant(ax),constant(bx),constant(kx),constant(zeros(length(ax))),constant(ay),constant(by),constant(ky),constant(zeros(length(bx))),constant(param.DELTAX),constant(param.DELTAY),constant(param.DELTAT),
#                 srci,srcj,srctype,srcv)
#         sigmaxx_arr,sigmayy_arr,sigmaxy_arr,vx_arr,vy_arr,mem_arr =
#             write(sigmaxx_arr, i, sigmaxx), write(sigmayy_arr, i, sigmayy), write(sigmaxy_arr, i, sigmaxy),
#                     write(vx_arr, i, vx), write(vy_arr, i, vy), write(mem_arr, i, mem)
#         op = tf.cond(tf.equal(tf.math.floormod(i,param.IT_DISPLAY),1), ()->tf.print("time step = ", i-1, " out of $NSTEP"), ()->tf.no_op())
#         i = bind(i, op)
#         i+1, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arr
#     end
#     i = constant(2, dtype=Int32)
#     _, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arr =
#         while_loop(condition, body, [i, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arr])
#     ElasticPropagator(
#             param,
#             src,
#             vx_arr,
#             vy_arr,
#             sigmaxx_arr,
#             sigmayy_arr,
#             sigmaxy_arr,
#             mem_arr
#     )
# end

################################################ ACOUSTIC ##################################################################

# todo:
# one_step(::AcousticPropagatorParams, ::PyCall.PyObject, ::PyCall.PyObject, ::PyCall.PyObject, ::PyCall.PyObject, ::PyCall.PyObject, ::PyCall.PyObject, ::PyCall.PyObject)
function acoustic_one_step_customop(param::AcousticPropagatorParams, w::PyObject, wold::PyObject, φ, ψ, σ::PyObject, τ::PyObject, c::PyObject)
    phi = φ
    psi = ψ
    sigma = σ
    tau = τ
    dt = param.DELTAT
    hx, hy = param.DELTAX, param.DELTAY
    nx, ny = param.NX, param.NY
    acoustic_one_step_ = load_op_and_grad("$(@__DIR__)/../deps/CustomOps/build/libADSeismic", "acoustic_one_step", multiple=true)
    w,wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny = convert_to_tensor(Any[w,wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny], [Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Int64,Int64])
    u, φ, ψ = acoustic_one_step_(w,wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny)
    set_shape(u, (param.NX+2)*(param.NY+2)), set_shape(φ, (param.NX+2)*(param.NY+2)), set_shape(ψ, (param.NX+2)*(param.NY+2))
end

function acoustic_one_step_customop_cpu(param::AcousticPropagatorParams, w::PyObject, wold::PyObject, φ, ψ, σ::PyObject, τ::PyObject, c::PyObject)
    phi = φ
    psi = ψ
    sigma = σ
    tau = τ
    dt = param.DELTAT
    hx, hy = param.DELTAX, param.DELTAY
    nx, ny = param.NX, param.NY
    acoustic_one_step_ = load_op_and_grad("$(@__DIR__)/../deps/CustomOps/build/libADSeismic", "acoustic_one_step_cpu", multiple=true)
    w,wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny = convert_to_tensor(Any[w,wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny], [Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Int64,Int64])
    u, φ, ψ = acoustic_one_step_(w,wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny)
    set_shape(u, (param.NX+2)*(param.NY+2)), set_shape(φ, (param.NX+2)*(param.NY+2)), set_shape(ψ, (param.NX+2)*(param.NY+2))
end

function acoustic_one_step_customop_ref(param::AcousticPropagatorParams, w::PyObject, wold::PyObject, φ, ψ, σ::PyObject, τ::PyObject, c::PyObject)
    Δt = param.DELTAT
    hx, hy = param.DELTAX, param.DELTAY
    IJ, IpJ, InJ, IJp, IJn, IpJp, IpJn, InJp, InJn =
        param.IJ, param.IpJ, param.InJ, param.IJp, param.IJn, param.IpJp, param.IpJn, param.InJp, param.InJn

    u = (2 - σ[IJ]*τ[IJ]*Δt^2 - 2*Δt^2/hx^2 * c[IJ] - 2*Δt^2/hy^2 * c[IJ]) * w[IJ] +
            c[IJ] * (Δt/hx)^2  *  (w[IpJ]+w[InJ]) +
            c[IJ] * (Δt/hy)^2  *  (w[IJp]+w[IJn]) +
            (Δt^2/(2hx))*(φ[IpJ]-φ[InJ]) +
            (Δt^2/(2hy))*(ψ[IJp]-ψ[IJn]) -
            (1 - (σ[IJ]+τ[IJ])*Δt/2) * wold[IJ]
    u = u / (1 + (σ[IJ]+τ[IJ])/2*Δt)
    u = scatter_nd_ops(IJ, u, (param.NX+2)*(param.NY+2))
    φ = (1. -Δt*σ[IJ]) * φ[IJ] + Δt * c[IJ] * (τ[IJ] -σ[IJ])/2hx *
        (w[IpJ]-w[InJ])
    ψ = (1. -Δt*τ[IJ]) * ψ[IJ] + Δt * c[IJ] * (σ[IJ] -τ[IJ])/2hy *
        (w[IJp]-w[IJn])
    φ = scatter_nd_ops(IJ, φ, (param.NX+2)*(param.NY+2))
    ψ = scatter_nd_ops(IJ, ψ, (param.NX+2)*(param.NY+2))
    u, φ, ψ
end

function one_step(param::AcousticPropagatorParams, w::PyObject, wold::PyObject, φ, ψ, σ::PyObject, τ::PyObject, c::PyObject)
    Δt = param.DELTAT
    hx, hy = param.DELTAX, param.DELTAY
    IJ, IpJ, InJ, IJp, IJn, IpJp, IpJn, InJp, InJn =
        param.IJ, param.IpJ, param.InJ, param.IJp, param.IJn, param.IpJp, param.IpJn, param.InJp, param.InJn

    u = (2 - σ[IJ]*τ[IJ]*Δt^2 - 2*Δt^2/hx^2 * c[IJ] - 2*Δt^2/hy^2 * c[IJ]) * w[IJ] +
            c[IJ] * (Δt/hx)^2  *  (w[IpJ]+w[InJ]) +
            c[IJ] * (Δt/hy)^2  *  (w[IJp]+w[IJn]) +
            (Δt^2/(2hx))*(φ[IpJ]-φ[InJ]) +
            (Δt^2/(2hy))*(ψ[IJp]-ψ[IJn]) -
            (1 - (σ[IJ]+τ[IJ])*Δt/2) * wold[IJ]
    u = u / (1 + (σ[IJ]+τ[IJ])/2*Δt)
    u = scatter_nd_ops(IJ, u, (param.NX+2)*(param.NY+2))
    φ = (1. -Δt*σ[IJ]) * φ[IJ] + Δt * c[IJ] * (τ[IJ] -σ[IJ])/2hx *
        (u[IpJ]-u[InJ])
    ψ = (1. -Δt*τ[IJ]) * ψ[IJ] + Δt * c[IJ] * (σ[IJ] -τ[IJ])/2hy *
        (u[IJp]-u[IJn])
    φ = scatter_nd_ops(IJ, φ, (param.NX+2)*(param.NY+2))
    ψ = scatter_nd_ops(IJ, ψ, (param.NX+2)*(param.NY+2))
    u, φ, ψ
end

"""
AcousticPropagatorSolver(param::AcousticPropagatorParams, src::AcousticSource, c::Union{PyObject, Array{Float64, 2}})

# Arguments
- 'param::AcousticPropagatorParams': model configuration parameters, such as NX, NT
- 'src::AcousticSource': source locations
- 'c::Union{PyObject, Array{Float64, 2}}': velocity model

# Return
- 'AcousticPropagator': forward calculation results, like displacement u
"""
function AcousticPropagatorSolver(param::AcousticPropagatorParams, src::AcousticSource, c::Union{PyObject, Array{Float64, 2}})

    c = tf.reshape(convert_to_tensor(c), (-1,))^2 ## convert to c^2
    compute_PML_Params!(param)

    σij = constant(param.Σx'[:])
    τij = constant(param.Σy'[:])

    tu = TensorArray(param.NSTEP+1; clear_after_read=false)
    tφ = TensorArray(param.NSTEP+1; clear_after_read=true)
    tψ = TensorArray(param.NSTEP+1; clear_after_read=true)

    function condition(i, ta, tφ, tψ)
        i<=param.NSTEP+1
    end

    function body(i, ta, tφ, tψ)

        if param.PropagatorKernel==0
            one_step_ = one_step
        elseif param.PropagatorKernel==1 && HAS_GPU
            @info "Use GPU custom one step..."
            one_step_ = acoustic_one_step_customop
        elseif param.PropagatorKernel==1 && !HAS_GPU
            @info "Use CPU custom one step..."
            one_step_ = acoustic_one_step_customop_cpu
        elseif param.PropagatorKernel==2
            @info "Use reference one step..."
            one_step_ = acoustic_one_step_customop_ref
        end
        u, φ, ψ = one_step_(param, read(ta, i-1), read(ta, i-2), read(tφ, i-1), read(tψ, i-1), σij, τij, c)
        srci, srcj, srcv = AcousticSourceAtTimeT(src, i-1)

        if param.IT_DISPLAY>0
            op = tf.cond(tf.equal(tf.math.floormod(i,param.IT_DISPLAY),1), ()->tf.print("time step = ", i-1, " out of $(param.NSTEP)"), ()->tf.no_op())
            i = bind(i, op)
        end

        src_index = (srci - 1) * (param.NY+2) + srcj
        u = scatter_add_op(u, src_index, srcv*param.DELTAT^2)

        ta_, tφ_, tψ_ = write(ta, i, u), write(tφ, i, φ), write(tψ, i, ψ)
        i+1, ta_, tφ_, tψ_
    end

    tu = write(tu, 1, constant(zeros((param.NX+2)*(param.NY+2))))
    tφ = write(tφ, 1, constant(zeros((param.NX+2)*(param.NY+2))))
    tψ = write(tψ, 1, constant(zeros((param.NX+2)*(param.NY+2))))
    tu = write(tu, 2, constant(zeros((param.NX+2)*(param.NY+2))))
    tφ = write(tφ, 2, constant(zeros((param.NX+2)*(param.NY+2))))
    tψ = write(tψ, 2, constant(zeros((param.NX+2)*(param.NY+2))))
    i = constant(3, dtype=Int32)
    _, tu, tφ, tψ = while_loop(condition, body, [i,tu,tφ,tψ]; swap_memory=true)
    tu = tf.reshape(stack(tu), (param.NSTEP+1, param.NX+2, param.NY+2))
    tφ = tf.reshape(stack(tφ), (param.NSTEP+1, param.NX+2, param.NY+2))
    tψ = tf.reshape(stack(tψ), (param.NSTEP+1, param.NX+2, param.NY+2))

    AcousticPropagator(param, src, tu, tφ, tψ)
end

function compute_PML_Params!(param::AcousticPropagatorParams)
    NX, NY = param.NX, param.NY
    # computing damping coefficient
    c, R = param.vp_ref, param.Rcoef
    Lx = param.NPOINTS_PML * param.DELTAX
    Ly = param.NPOINTS_PML * param.DELTAY
    param.damping_x = c/Lx*log(1/R)
    param.damping_y = c/Ly*log(1/R)
    # @show c, Lx, log(1/R), param.damping_x, param.damping_y

    X = (0:param.NX+1)*param.DELTAX
    Y = (0:param.NY+1)*param.DELTAY
    Σx = zeros(param.NX+2, param.NY+2)
    Σy = zeros(param.NX+2, param.NY+2)
    for i = 1:param.NX+2
        for j = 1:param.NY+2
            Σx[i,j], Σy[i,j] = pml_helper(X[i], Y[j], param)
        end
    end

    param.Σx = Σx
    param.Σy = Σy
    param.IJ = getid2(2:NX+1, 2:NY+1, NX, NY)
    param.IJn = getid2(2:NX+1, 1:NY, NX, NY)
    param.IJp = getid2(2:NX+1, 3:NY+2, NX, NY)
    param.IpJ = getid2(3:NX+2, 2:NY+1, NX, NY)
    param.IpJp = getid2(3:NX+2, 3:NY+2, NX, NY)
    param.IpJn = getid2(3:NX+2, 1:NY, NX, NY)
    param.InJ = getid2(1:NX, 2:NY+1, NX, NY)
    param.InJn = getid2(1:NX, 1:NY, NX, NY)
    param.InJp = getid2(1:NX, 3:NY+2, NX, NY)
    return param
end

# short functions #
function ElasticReceiver(rcvi::Union{Array{Int64, 1},PyObject}, rcvj::Union{Array{Int64, 1},PyObject},
    rcvtype::Union{Array{Int64, 1},PyObject})
    ElasticReceiver(rcvi, rcvj, rcvtype, missing)
end

function ElasticSourceAtTimeT(src::ElasticSource,i::PyObject)
    convert_to_tensor(src.srci), convert_to_tensor(src.srcj),
        convert_to_tensor(src.srctype), convert_to_tensor(src.srcv)[i]
end

function getid(param, a, b)
    idx = Int64[]
    for i = 1:length(a)
        for j = 1:length(b)
            push!(idx, (a[i]-1)*(param.NY+2)+b[j])
        end
    end
    idx
end

function adbroadcast(a, b, idx=1)
    if idx==2
        n = div(length(b), length(a))
        repeat(a, n)*b
    elseif idx==1
        n = div(length(b), length(a))
        repeat(a', n)[:]*b
    elseif idx==4
        n = div(length(a), length(b))
        a/repeat(b, n)
        return a
    elseif idx==3
        n = div(length(a), length(b))
        a/repeat(b', n)[:]
        return a
    end
end

function makevector(m::PyObject, ij::Array{Int64}, o::PyObject)
    scatter_nd_ops(ij, o, length(m))
end

function SimulatedObservation!(ep::ElasticPropagator, rcv::ElasticReceiver)
    N = (ep.param.NX+2)*(ep.param.NY+2)
    get_receive = load_op_and_grad("$(@__DIR__)/../deps/CustomOps/build/libADSeismic", "get_receive")
    u = get_receive(tf.reshape(ep.vx, (ep.param.NSTEP+1, N)),
                    tf.reshape(ep.vy, (ep.param.NSTEP+1, N)),
                    tf.reshape(ep.sigmaxx, (ep.param.NSTEP+1, N)),
                    tf.reshape(ep.sigmayy, (ep.param.NSTEP+1, N)),
                    tf.reshape(ep.sigmaxy, (ep.param.NSTEP+1, N)),
                    convert_to_tensor(rcv.rcvi), convert_to_tensor(rcv.rcvj), convert_to_tensor(rcv.rcvtype),
                        constant(ep.param.NX), constant(ep.param.NY))
    rcv.rcvv = u
end

function AcousticSourceAtTimeT(src::AcousticSource,i::PyObject)
    convert_to_tensor(src.srci, dtype=Int64), convert_to_tensor(src.srcj, dtype=Int64), convert_to_tensor(src.srcv)[i]
end

"""
    SimulatedObservation!(ap::AcousticPropagator, rcv::AcousticReceiver)

Extract and save simulated displacement u into rcv::AcousticReceiver.
"""
function SimulatedObservation!(ap::AcousticPropagator, rcv::AcousticReceiver)
    idx = (rcv.rcvi .- 1) *  (ap.param.NY+2) + rcv.rcvj
    u = reshape(ap.u, ap.param.NSTEP+1, (ap.param.NX+2)*(ap.param.NY+2))[:, idx]
    rcv.rcvv = u
end

"""
    pml_helper(x::Float64, y::Float64, param::AcousticPropagatorParams)

Computing the PML profile.
"""
function pml_helper(x::Float64, y::Float64, param::AcousticPropagatorParams)
    outx = 0.0; outy = 0.0
    ξx = param.damping_x
    Lx = param.NPOINTS_PML * param.DELTAX
    if x<Lx && param.USE_PML_XMIN
        d = abs(Lx-x)
        outx = ξx * (d/Lx - sin(2π*d/Lx)/(2π))
    elseif x>param.DELTAX*(param.NX+1)-Lx && param.USE_PML_XMAX
        d = abs(x-(param.DELTAX*(param.NX+1)-Lx))
        outx = ξx * (d/Lx - sin(2π*d/Lx)/(2π))
    end

    ξy = param.damping_y
    Ly = param.NPOINTS_PML * param.DELTAY
    if y<Ly && param.USE_PML_YMIN
        d = abs(Ly-y)
        outy = ξy * (d/Ly - sin(2π*d/Ly)/(2π))
    elseif y>param.DELTAY*(param.NY+1)-Ly && param.USE_PML_YMAX
        d = abs(y-(param.DELTAY*(param.NY+1)-Ly))
        outy = ξy * (d/Ly - sin(2π*d/Ly)/(2π))
    end

    return outx, outy
end

function getid2(a, b, NX, NY)
    idx = Int64[]
    for i = 1:length(a)
        for j = 1:length(b)
            push!(idx, (a[i]-1)*(NY+2)+b[j])
        end
    end
    return idx
end


text:
#   This unit allows to precompile the FourierAnalysis package
#   and its dependency to the PosDefManifold package locally.
#   v 0.O.1 - last update 10th of July  2019
#
#   MIT License
#   Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home
#
#   DIRECTIONS:
#   1) If you have installed the FourierAnalysis from github or Julia registry,
#      uninstall it.
#   2) Change the `juliaCodeDir` path here below to the path
#           where the top FourierAnalysis and PosDefMaifold folders
#           are located on your computer.
#   3) Under Linux, replace all '\\' with `/`
#   4) Put the cursor in this unit and hit SHIFT+CTRL+ENTER
#
#   Nota Bene: all you need is actually the 'push' lines and
#   the 'using' line. You can safely delete the rest once
#   you have identified the 'srcDir' to be used in the push command.

begin
    # change the 'juliaCodeDir' path to the folder where your projects are
    juliaCodeDir= homedir()*"\\Documents\\Code\\julia\\"

    # push!(LOAD_PATH, juliaCodeDir*"PosDefManifold"*"\\src\\")
    using Documenter, DocumenterTools, BenchmarkTools, Revise # PosDefManifold,
    # add simple modules to be used
    push!(LOAD_PATH, juliaCodeDir*"Modules")
    using IOtxt

    projectName = "FourierAnalysis"
    srcDir      = juliaCodeDir*projectName*"\\src\\"
    docsDir     = juliaCodeDir*projectName*"\\docs\\"
    push!(LOAD_PATH, srcDir)
    using FourierAnalysis

    # for compiling the documentation
    cd(docsDir)
    clipboard("""makedocs(sitename="$projectName", modules=[$projectName])""")
    @info("\nhit CTRL+V+ENTER on the REPL for building the documentation.");
end

max_stars_repo_path: src/START HERE.jl

max_issues_repo_path: src/START HERE.jl

max_forks_repo_path: src/START HERE.jl


text:
using Distributions
using LinearAlgebra

# Scene loader & Predictor parameters
scene_mode = "synthetic";                                                           # "synthetic" or "data"
prediction_mode = "gaussian";                                                       # "gaussian" or "trajectron" or "oracle"
prediction_device = "cpu";                                                          # "cpu" or "cuda"
prediction_steps = 12;                                                              # number of steps to look ahead in the future
ado_pos_init_dict = Dict("PEDESTRIAN/1" => [0.0, -5.0]);                            # initial ado positions [x, y] [m]
ado_vel_dict = Dict("PEDESTRIAN/1" => MvNormal([0.0, 1.0], Diagonal([0.01, 0.15])));# ado velocity distributions
dto = 0.4;                                                                          # observation update time interval [s]
prediction_rng_seed = 1;                                                            # random seed for prediction (and stochastic transition for "synthetic" scenes)
deterministic = false;                                                              # if true, a single, deterministic sample is drawn regardless of random seed. (num_samples = 1 is needed)
num_samples = 30;                                                                   # number of trajectory samples (per ado agent)
# Cost Parameters
include("params_cost.jl")
# Control Parameters
include("params_control.jl")
# Ego initial state
ego_pos_init_vec = [-5., 0.];                                                       # initial ego position [x, y] [m]
ego_pos_goal_vec = [5., 0.];                                                        # goal ego position [x, y] [m]
# Other parameters
pos_error_replan = 5.0;                                                             # position error for replanning target trajectory [m]
target_speed = 1.0;                                                                 # target speed [m/s]
sim_horizon = 10.0;                                                                 # simulation time horizon [s]


text:
@testset "mapmany" begin
    @test mapmany(i->1:i, 1:3) == [1, 1,2, 1,2,3]
    @test mapmany((i,j)->1:(i+j), 1:3, 1:3) == [1,2, 1,2,3,4, 1,2,3,4,5,6]
end

@testset "flatten" begin
    @test flatten([1:1, 1:2, 1:3]) == [1, 1,2, 1,2,3]
end

@testset "mapview" begin
    # Arrays
    a = [1,2,3]
    @test @inferred(mapview(-, a)) isa MappedArray{Int,1}
    b = mapview(-, a)
    @test b == [-1,-2,-3]
    a[1] = 11
    @test b == [-11,-2,-3]

    # Iterables
    iter = zip([1,2,3], [1,2,3])
    iter2 = mapview(x -> x[1] + x[2], iter)
    @test @inferred(first(iter2)) === 2
    @test collect(iter2) == [2, 4, 6]
end


text:
# QuantLib Time module
module Time

using Dates

import Dates.adjust

# Date helper funcs
within_next_week(d1::Date, d2::Date) = d2 >= d1 && d2 <= d1 + Dates.Day(7)
within_next_week(t1::Float64, t2::Float64) = t1 <= t2 && t2 <= t1 + (1.0/52.0)
within_previous_week(d1::Date, d2::Date) = d2 >= d1 - Dates.Day(7) && d2 <= d1

export within_next_week, within_previous_week

# Frequency.jl
export Frequency, NoFrequency, Once, Annual, Semiannual, EveryFourthMonth, Quarterly, Bimonthly, Monthly, EveryFourthWeek, Biweekly, Weekly, Daily, OtherFrequency,
value, period

# DayCount.jl
export DayCount, Actual360, Actual365, Thirty360, BondThirty360, EuroBondThirty360, ItalianThirty360, ActualActual, SimpleDayCount, ISMAActualActual, ISDAActualActual,
day_count, days_per_year, year_fraction

# BusinessCalendar.jl
export BusinessCalendar, WesternCalendar, OrthodoxCalendar, UnitedStatesCalendar, USSettlementCalendar, USNYSECalendar, USNERCCalendar,
USGovernmentBondCalendar, UnitedKingdomCalendar, UKSettlementCalendar, UKLSECalendar, UKLMECalendar, TargetCalendar, NullCalendar, JointCalendar,
BusinessDayConvention, Unadjusted, ModifiedFollowing, Following, is_business_day,
easter_date, is_holiday, advance, adjust

# tenor_period.jl
export TenorPeriod

# schedule.jl
export DateGenerationRule, DateGenerationForwards, DateGenerationBackwards, DateGenerationTwentieth, Schedule

# time_grid.jl
export TimeGrid, is_empty, closest_time, return_index

include("Frequency.jl")
include("DayCount.jl")
include("BusinessCalendar.jl")
include("tenor_period.jl")
include("schedule.jl")
include("time_grid.jl")

end


text:
module JCSA
import Base.getindex
import Base.setindex!
import Base.length
using Distributions
using ImageFiltering
using LightGraphs
using ImageSegmentation
include("ExponentialDistributions.jl")
include("combined_kmeans.jl")
import RGBDSegmentation: clusterize,
    CDNImage,
    rgb2lab!,
    colors,
    distances,
    normals

@enum ComponentType Color=1 Depth=2 Normal=3
struct Config
    kp::Float64
    thd::Float64
    thb::Float64
    thr::Float64
    max_iterations::Integer
end

"""
Construct a JCSA Config

# Arguments
- `kp::Float64`  Threshold  to decide a region as planar
- `thd::Float64` Threshold to decide the distance among two regions
- `thb::Float64` Threshold to decide the existence of boundary among two regions
- `thr::Float64` Threshold to determine the goodness of a plane fitting.
"""
function Config(; kp::Float64=5.0, thd::Float64=2.0, thb::Float64=3.0, thr::Float64=0.9, max_iterations::Integer=500)
    return Config(kp,thd,thb,thr,max_iterations)
end

"""
https://gist.github.com/manojpandey/f5ece715132c572c80421febebaf66ae
"""

const CDNComponentExponential = Tuple{MvGaussianExponential,
                                      MvGaussianExponential,
                                      MvWatsonExponential}
const  CDNComponent = Tuple{MvNormal,
                            MvNormal,
                            MvWatson}

const SufficientStatisticsGaussian  = Tuple{Array{<:AbstractFloat,2},
                                        Array{<:AbstractFloat, 3} }
const SufficientStatisticsNormal = Tuple{Matrix{<:AbstractFloat}}
struct CDNData
    color::SufficientStatisticsGaussian
    depth::SufficientStatisticsGaussian
    normal::SufficientStatisticsNormal
end
function length(d::CDNData)
    return size(d.color[1],2)
end
function getindex(d::CDNData,n::ComponentType)
    if n==Color
        return  d.color
    elseif n==Depth
        return d.depth
    elseif n==Normal
        return d.normal
    end
end
"""

"""
struct Model
    weights::Vector{Float64}
    components::Vector{CDNComponentExponential}
end
"""
```
function Model(k::Integer)

```
Construct a JCSA Model composed by a mixture of `k` components. Each component is a product of two multivariate gaussian and a multivariate watson distribution
"""
function Model(initial_assignments::Vector, xss::CDNData)
    (_, n) = size(xss.color[1])
    weights = []
    components = []
    for cluster in initial_assignments
        try
            push!(components, (
                MvGaussianExponential(
                    xss.color[1][:, cluster]),
                MvGaussianExponential(
                    xss.depth[1][:, findall(cluster)]),
                MvWatsonExponential(
                    vec(mean(xss.normal[1][:,cluster],dims=2)))
            ))
            push!(weights, length(cluster)/n)
        catch e
            @info "Model construction $e"
        end
    end
    return Model(weights, components)
end

include("RegionMerging.jl")

"""
```
function negative_log_likelihood(m::Model, data)::Float64
```
    """
function negative_log_likelihood(m::Model, data)::Float64
    k = length(m.components)
    n = size(data,2)
    pdfs = zeros(n,k,3)
    indices = [1:3, 4:6, 7:9]
    for j=1:k
        try
            distributions = to_source_parameters(m.components[j])
            for i=1:3
                Distributions.logpdf!(view(pdfs, :, j, i),
                                      distributions[i],
                                      data[indices[i], :])
            end
        catch e
        end
    end
    try
    logPDF = log.(m.weights)'.+ dropdims(sum(pdfs, dims=3), dims=3)
    y = maximum(logPDF, dims=2);
    x = y.-logPDF
    s = y + log.(sum(exp.(logPDF), dims=2));
    i = findall(.!isfinite.(y));
    if !isempty(i)
        s[i]= y[i]
    end
    return -sum(s)
    catch e
        return 500
    end
end
import Base.+
"""
```
function expectation(pij::Matrix{Float64}, DG::Array{<:AbstractFloat,3},  m::Model, x::CDNData)

```
"""
function expectation(pij::Matrix{Float64}, DG::Array{<:AbstractFloat,3},  m::Model, x::CDNData)
    k = length(m.weights)
    npatterns = length(x)
    div_sum = 0
    for j=1:k, r=1:3
        DG[j,r,:] = bregman_divergence(
            m.components[j][r],
            x[ComponentType(r)]...)
    end
    for j=1:k
        dg_k = dropdims(sum(DG[j,:,:], dims=1), dims=1)
        pij[j,:] = m.weights[j]*exp.(dg_k);
    end
    pij[isinf.(pij)] .= 1e7
    broadcast!(/,pij,pij,sum(pij,dims=1))
end
"""
```
function maximization(component::MvGaussianExponential , pij::Matrix{<:AbstractFloat}, data, j::Integer,sumpij::Vector )
```
Compute the maximization step of a multivariate gaussian
"""
function maximization(component::MvGaussianExponential , pij::Matrix{<:AbstractFloat}, data, j::Integer,sumpij::Vector )
    dim = size(data[2],1)
    n = size(data[2],3)
    data2 = reshape(data[2],dim*dim, n)
    n2 = sum(view(pij,j,:).*data2', dims=1)./sumpij[j]
    n1 = vec(sum(view(pij,j,:).*data[1]', dims=1)./sumpij[j])
    n2 = reshape(n2, (dim,dim))
    component.n =(n1,n2)
end
"""
```
function maximization(component::MvWatsonExponential , pij::Matrix{<:AbstractFloat}, data, j::Integer,sumpij::Vector )
```
Perform the maximization step of the directional data
"""
function maximization(component::MvWatsonExponential , pij::Matrix{<:AbstractFloat}, data, j::Integer,sumpij::Vector )
    component.n =  vec(sum(view(pij,j,:).*data[1]',dims=1)./sumpij[j])
end
"""
```
function maximization(m::Model,pij::Matrix{<:AbstractFloat}, x::CDNData)
```
"""

function maximization(m::Model,pij::Matrix{<:AbstractFloat}, x::CDNData)
    k = length(m.weights)
    nvectors = length(x)
    sumpij = vec(sum(pij, dims=2))
    #=Threads.@threads=# for j=1:k
        maximization(m.components[j][1], pij, x.color, j, sumpij)
        maximization(m.components[j][2], pij, x.depth, j, sumpij)
        maximization(m.components[j][3], pij, x.normal, j, sumpij)
        m.weights[j] = (sumpij[j]) /nvectors
    end
end

function update_parameters!(m::Model, pij)
    k = length(m.weights)
    components_to_remove = []
    invalid_component_lock = Threads.SpinLock()
   #= Threads.@threads=# for i=1:k
       try
           for j=1:3
               update_parameters!(m.components[i][j])
           end
       catch e
	   lock(invalid_component_lock)
           push!(components_to_remove, i)
           @warn("Component removed $e")
 	   unlock(invalid_component_lock)
       end
    end
    deleteat!(m.weights, components_to_remove)
    deleteat!(m.components, components_to_remove)
    for j in reverse(components_to_remove)
        pij = pij[1:end .!= j,: ]
    end
end

function get_extended_data(data::Array{T, 2}) where T<:Number
    color = sufficient_statistic(MvGaussianExponential, data[1:3, :])
    depth = sufficient_statistic(MvGaussianExponential, data[4:6, :])
    normal = sufficient_statistic(MvWatsonExponential, data[7:9, :])
    data = CDNData(color,depth,(normal, ))
    return data

end

function hard_clustering( pij::Matrix{<:AbstractFloat}, k::Integer, img_size)
    labels = [d[1] for d in findmax(pij, dims=1)[2]]
    labels = reshape(labels,img_size[1],img_size[2])
    return  round.(Integer,mapwindow(median!,labels, (3,3)))
end
"""
```
function clusterize(cfg::Config, img_a::CDNImage, k::Integer)
```
Entry point of the algorithm.

# Parameters
- `cfg::Config`: Config
- `img_a::CDNImage`: Color-Distance-Normals image of the scene.
- `k::Integer`: Number of clust

"""
function clusterize(cfg::Config, img::CDNImage, k::Integer)
    (model, pij, extended_data) = expectation_maximization(cfg, img, k)
    k = length(model.components)
    labels = hard_clustering(pij, k, (size(img,2), size(img,3)))
    rm_config = RegionMergingConfig()

    labels = process_labels(rm_config,
                            colors(img),
                            distances(img),
                            labels,
                            extended_data)
    return SegmentedImage(img, labels)
end
function expectation_maximization(cfg::Config,
                                  img::CDNImage,
                                  k::Integer) where T <:AbstractArray

    (nch,nr,nc) = size(img)
    npatterns = nr*nc
    p_data = reshape(img,
                     (nch, nr*nc))
    nk = k
    assignments = initial_clusters(p_data,k)
    xss = get_extended_data(p_data)
    model = Model([vec(assignments.==i) for i=1:k], xss)
    pij = zeros(Float64, k,npatterns)
    prev_nll = Inf
    current_nll = negative_log_likelihood(model, p_data)
    iteration = 1
    prev_pij = nothing
    DG = zeros(k,3,npatterns)
    current_nll = Inf
    best_nll = Inf
    while iteration < 15 || (iteration<cfg.max_iterations && current_nll<prev_nll)
        expectation(pij, DG, model, xss)
        maximization(model,pij,xss )
	update_parameters!(model, pij)
        prev_nll    = current_nll
        current_nll = negative_log_likelihood(model, p_data)
        if current_nll < best_nll
            best_nll = current_nll
        end
        @info """Iteration: $iteration
              | Current improvement NLL: $(prev_nll - current_nll)
              | Best: $best_nll
              """
        iteration = iteration + 1
    end
   return (model,pij, xss)
end
end


text:
# ---
# title: 342. Power of Four
# id: problem342
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/power-of-four/description/>
# hidden: true
# ---
#
# Given an integer `n`, return _`true` if it is a power of four. Otherwise,
# return `false`_.
#
# An integer `n` is a power of four, if there exists an integer `x` such that `n
# == 4x`.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 16
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: n = 5
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: n = 1
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `-231 <= n <= 231 - 1`
#
#
#
# **Follow up:** Could you solve it without loops/recursion?
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module TestMixins
using Test
using Mixins
using MLStyle: @match

@define_mixin f begin
  x::Int
  y::String
end

@define_mixin g begin
  z::Symbol
end

@use_mixins struct Cocktail
  @mixin(f)
  @mixin(g)
end

c = Cocktail(2,"gin",:tonic)
@test c.x == 2
@test c.y == "gin"
@test c.z == :tonic

macro modifyfields(structdef)
  modifieddef = @match structdef begin
    Expr(:struct, parent, name, body) => begin
      modifiedlines = map(body.args) do line
        @match line begin
          Expr(:(::), field, type) => :($(Symbol(string(field)*"prime")) :: $(type))
          _ => line
        end
      end
      Expr(:struct, parent, name, Expr(:block, modifiedlines...))
    end
    _ => error("expects a struct")
  end
  esc(modifieddef)
end

@use_mixins @modifyfields struct CocktailPrime
  @mixin(f)
  @mixin(g)
end

c = CocktailPrime(2,"gin",:tonic)
@test c.xprime == 2
@test c.yprime == "gin"
@test c.zprime == :tonic

end


text:
using JLD2, FileIO

function save(fname,args...;kwargs...)
     FileIO.save(fname,serialize.(args)...;kwargs...)
end
function load(fname,args...;kwargs...)
     serialize(FileIO.load(fname,args...;kwargs...))
end

macro save(filename, vars...)
    if isempty(vars)
        # Save all variables in the current module
        quote
            let
                m = $(__module__)
                f = jldopen($(esc(filename)), "w")
                wsession = JLD2.JLDWriteSession()
                try
                    for vname in names(m; all=true)
                        s = string(vname)
                        if !occursin(r"^_+[0-9]*$", s) # skip IJulia history vars
                            v = getfield(m, vname)
                            if !isa(v, Module)
                                try
                                    write(f, s, serialize(v), wsession)
                                catch e
                                    if isa(e, PointerException)
                                        @warn("skipping $vname because it contains a pointer")
                                    else
                                        rethrow(e)
                                    end
                                end
                            end
                        end
                    end
                finally
                    close(f)
                end
            end
        end
    else
        writeexprs = Vector{Expr}(undef, length(vars))
        for i = 1:length(vars)
            writeexprs[i] = :(write(f, $(string(vars[i])), serialize($(esc(vars[i]))), wsession))
        end

        quote
            jldopen($(esc(filename)), "w") do f
                wsession = JLD2.JLDWriteSession()
                $(Expr(:block, writeexprs...))
            end
        end
    end
end

macro load(filename, vars...)
    if isempty(vars)
        if isa(filename, Expr)
            throw(ArgumentError("filename argument must be a string literal unless variable names are specified"))
        end
        # Load all variables in the top level of the file
        readexprs = Expr[]
        vars = Symbol[]
        f = jldopen(filename)
        try
            for n in keys(f)
                if !JLD2.isgroup(f, JLD2.lookup_offset(f.root_group, n))
                    push!(vars, Symbol(n))
                end
            end
        finally
            close(f)
        end
    end
    return quote
        ($([esc(x) for x in vars]...),) = jldopen($(esc(filename))) do f
            ($([:(serialize(read(f, $(string(x))))) for x in vars]...),)
        end
        $(Symbol[v for v in vars]) # convert to Array
    end
end


text:
#Works for bar and grouped bar charts, add more cases if it makes sense
function coord_flip!(v::VegaVisualization)

    #Switch scales
    v.scales[findfirst([z.name == "x" for z in v.scales])].range = "height"
    v.scales[findfirst([z.name == "y" for z in v.scales])].range = "width"

    #Switch axes and title names
    v.axes[findfirst([z._type == "x" for z in v.axes])].scale = "y"
    v.axes[findfirst([z._type == "x" for z in v.axes])].title = "y"
    v.axes[findfirst([z._type == "y" for z in v.axes])].scale = "x"
    v.axes[findfirst([z._type == "y" for z in v.axes])].title = "x"

    #Change mark properties

    v.marks[1].properties.enter.y = VegaValueRef(scale = "x", field = "x")

    _offset = v.marks[1].properties.enter.width.offset
    v.marks[1].properties.enter.height = VegaValueRef(scale = "x", band = true, offset = _offset)
    v.marks[1].properties.enter.y2 = nothing

    if isempty(filter(x-> x.name == "stats", v.data))
        v.marks[1].properties.enter.x = VegaValueRef(scale = "y", field = "y")
        v.marks[1].properties.enter.x2 = VegaValueRef(scale = "y", value = 0)
    else  # stacked
        v.marks[1].properties.enter.x = VegaValueRef(scale = "y", field = "layout_start")
        v.marks[1].properties.enter.x2 = VegaValueRef(scale = "y", field = "layout_end")
    end

    #Grouped bar hack
    if v.name == "groupedbar"
        v.marks[1].marks[1].properties.enter.height = v.marks[1].marks[1].properties.enter.width
        v.marks[1].marks[1].properties.enter.width = nothing

        v.marks[1].marks[1].properties.enter.x2 = v.marks[1].marks[1].properties.enter.y2
        v.marks[1].marks[1].properties.enter.y2 = nothing

        v.marks[1].marks[1].properties.enter.y, v.marks[1].marks[1].properties.enter.x = v.marks[1].marks[1].properties.enter.x, v.marks[1].marks[1].properties.enter.y

        v.marks[1].properties.enter.width = nothing
        v.marks[1].properties.enter.x = nothing
        v.marks[1].properties.enter.x2 = nothing
        v.marks[1].properties.enter.y = VegaValueRef(scale = "x", field = "key")
        v.marks[1].scales[1].range = "height"
    end

    return v
end

function legend!(v::VegaVisualization; title::AbstractString="Group", show::Bool=true)

    show == true? v.legends = [VegaLegend(fill = "group", title = title)] : empty!(v.legends)
    return v

end

#Use ColorBrewer.jl scales
function colorscheme!(v::VegaVisualization; palette::Union{Tuple{AbstractString,Any}, AbstractString, Array} = "ordinal10", reversePalette::Bool = false)

    #See if group or color key exists
    i = findfirst([z.name == "group" for z in v.scales])
    i = i > 0? i: findfirst([z.name == "color" for z in v.scales]) #Was choropleth hack, standardized scale name to group

    #Copy scales property
    s = v.scales[i]

    #colorscheme, single color or array of arbitrary string colors
    if isa(palette, Tuple)
        reversePalette == true? s.range = reverse(colorpalettes[palette[1]]["$(palette[2])"]) : s.range = colorpalettes[palette[1]]["$(palette[2])"]
    elseif isa(palette, AbstractString)
        s.range = [palette]
    elseif isa(palette, Array)
        reversePalette == true? s.range = reverse(collect(palette)) : s.range = collect(palette)
    end

    return v

end

function stroke!(v::VegaVisualization; color::AbstractString = "Black", width::Real = 0.75, opacity::Real = 1, visible::Bool = true)

    if typeof(v.marks[1].marks) == Void
        if visible
            v.marks[1].properties.enter.stroke = VegaValueRef(value = color)
            v.marks[1].properties.enter.strokeWidth = VegaValueRef(value = width)
            v.marks[1].properties.enter.strokeOpacity = VegaValueRef(value = opacity)
        else
            v.marks[1].properties.enter.stroke = nothing
            v.marks[1].properties.enter.strokeWidth = nothing
            v.marks[1].properties.enter.strokeOpacity = nothing
        end
    #Group Marks
    else
        if visible
            v.marks[1].marks[1].properties.enter.stroke = VegaValueRef(value = color)
            v.marks[1].marks[1].properties.enter.strokeWidth = VegaValueRef(value = width)
            v.marks[1].marks[1].properties.enter.strokeOpacity = VegaValueRef(value = opacity)
        else
            v.marks[1].marks[1].properties.enter.stroke = nothing
            v.marks[1].marks[1].properties.enter.strokeWidth = nothing
            v.marks[1].marks[1].properties.enter.strokeOpacity = nothing
        end
    end

    return v
end

function hover!(v::VegaVisualization; opacity::Number = 1, color::AbstractString = "")

    #Hack due to boxplot having a different color setup than other visualizations
    v.name == "boxplot" ? field = "x" : field = "group"

    #Test if color is set in function, so that we know to reset the color on update
    color == ""? color_ = VegaValueRef(scale = "group", field = field) : color_ = VegaValueRef(value = color)

    #Grouped bar, line chart, popchart
    if v.marks[1].marks != nothing

        v.marks[1].marks[1].properties.hover = VegaMarkPropertySet(fillOpacity = VegaValueRef(value = opacity),
                                                                   strokeOpacity = VegaValueRef(value = opacity),
                                                                   fill = v.name == "lineplot"? nothing: color_,
                                                                   stroke = v.name == "lineplot"? color_ : nothing
                                                                   )
        v.marks[1].marks[1].properties.update = VegaMarkPropertySet(fillOpacity = VegaValueRef(value = 1),
                                                                    strokeOpacity = VegaValueRef(value = 1),
                                                                    fill = v.name == "lineplot"? nothing: VegaValueRef(scale = "group", field = field),
                                                                    stroke = v.name == "lineplot"? VegaValueRef(scale = "group", field = field) : nothing
                                                                    )

    #Area, Bar, Heatmap, Pie Chart, Scatterplot, Waterfall, WordCloud, Bubble Chart, Stemleaf, Choropleth
    else

        v.marks[1].properties.hover = VegaMarkPropertySet(fillOpacity = VegaValueRef(value = opacity),
                                                          strokeOpacity = VegaValueRef(value = opacity),
                                                          fill = color_
                                                        )
        v.marks[1].properties.update = VegaMarkPropertySet(fillOpacity = VegaValueRef(value = 1),
                                                           strokeOpacity = VegaValueRef(value = 1),
                                                           fill = v.name == "choropleth"? VegaValueRef(scale = "group", field = "table2.y"): VegaValueRef(scale = "group", field = field)
                                                        )

    end

    return v
end

function jitter!(v::VegaVisualization; pctXrange::Real = 0.05, pctYrange::Real = 0.05)

    #Calculate range, then use random pct of range
    _y = []
    _x = []

    for datum in v.data[1].values
        push!(_x, datum["x"])
        push!(_y, datum["y"])
    end

    xrange = maximum(_x) - minimum(_x)
    yrange = maximum(_y) - minimum(_y)

    #Iterate over VegaData to jitter dataset
    for datum in v.data[1].values
        if typeof(datum["x"]) <: Number
            datum["x"] = rand(0:0.001:pctXrange) * xrange + datum["x"]
        end
        if typeof(datum["y"]) <: Number
            datum["y"] = rand(0:0.001:pctYrange) * yrange + datum["y"]
        end
    end

    return v
end

max_stars_repo_name: jpthor/vega