text:
export hoeffding_tree

import ..HoeffdingTreeModel

using mlpack._Internal.io

import mlpack_jll
const hoeffding_treeLibrary = mlpack_jll.libmlpack_julia_hoeffding_tree

# Call the C binding of the mlpack hoeffding_tree binding.
function hoeffding_tree_mlpackMain()
  success = ccall((:hoeffding_tree, hoeffding_treeLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module hoeffding_tree_internal
  import ..hoeffding_treeLibrary

import ...HoeffdingTreeModel

# Get the value of a model pointer parameter of type HoeffdingTreeModel.
function IOGetParamHoeffdingTreeModel(paramName::String)::HoeffdingTreeModel
  HoeffdingTreeModel(ccall((:IO_GetParamHoeffdingTreeModelPtr, hoeffding_treeLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type HoeffdingTreeModel.
function IOSetParamHoeffdingTreeModel(paramName::String, model::HoeffdingTreeModel)
  ccall((:IO_SetParamHoeffdingTreeModelPtr, hoeffding_treeLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeHoeffdingTreeModel(stream::IO, model::HoeffdingTreeModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeHoeffdingTreeModelPtr, hoeffding_treeLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeHoeffdingTreeModel(stream::IO)::HoeffdingTreeModel
  buffer = read(stream)
  HoeffdingTreeModel(ccall((:DeserializeHoeffdingTreeModelPtr, hoeffding_treeLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    hoeffding_tree(; [batch_mode, bins, confidence, info_gain, input_model, labels, max_samples, min_samples, numeric_split_strategy, observations_before_binning, passes, test, test_labels, training, verbose])

This program implements Hoeffding trees, a form of streaming decision tree
suited best for large (or streaming) datasets.  This program supports both
categorical and numeric data.  Given an input dataset, this program is able to
train the tree with numerous training options, and save the model to a file.
The program is also able to use a trained model or a model from file in order to
predict classes for a given test set.

The training file and associated labels are specified with the `training` and
`labels` parameters, respectively. Optionally, if `labels` is not specified, the
labels are assumed to be the last dimension of the training dataset.

The training may be performed in batch mode (like a typical decision tree
algorithm) by specifying the `batch_mode` option, but this may not be the best
option for large datasets.

When a model is trained, it may be saved via the `output_model` output
parameter.  A model may be loaded from file for further training or testing with
the `input_model` parameter.

Test data may be specified with the `test` parameter, and if performance
statistics are desired for that test set, labels may be specified with the
`test_labels` parameter.  Predictions for each test point may be saved with the
`predictions` output parameter, and class probabilities for each prediction may
be saved with the `probabilities` output parameter.

For example, to train a Hoeffding tree with confidence 0.99 with data `dataset`,
saving the trained tree to `tree`, the following command may be used:

```julia
julia> tree, _, _ = hoeffding_tree(confidence=0.99,
            training=dataset)
```

Then, this tree may be used to make predictions on the test set `test_set`,
saving the predictions into `predictions` and the class probabilities into
`class_probs` with the following command:

```julia
julia> _, predictions, class_probs =
            hoeffding_tree(input_model=tree, test=test_set)
```

# Arguments

 - `batch_mode::Bool`: If true, samples will be considered in batch
      instead of as a stream.  This generally results in better trees but at the
      cost of memory usage and runtime.  Default value `false`.

 - `bins::Int`: If the 'domingos' split strategy is used, this specifies
      the number of bins for each numeric split.  Default value `10`.

 - `confidence::Float64`: Confidence before splitting (between 0 and 1).
      Default value `0.95`.

 - `info_gain::Bool`: If set, information gain is used instead of Gini
      impurity for calculating Hoeffding bounds.  Default value `false`.

 - `input_model::HoeffdingTreeModel`: Input trained Hoeffding tree model.
 - `labels::Array{Int, 1}`: Labels for training dataset.
 - `max_samples::Int`: Maximum number of samples before splitting.
      Default value `5000`.

 - `min_samples::Int`: Minimum number of samples before splitting.
      Default value `100`.

 - `numeric_split_strategy::String`: The splitting strategy to use for
      numeric features: 'domingos' or 'binary'.  Default value `binary`.

 - `observations_before_binning::Int`: If the 'domingos' split strategy is
      used, this specifies the number of samples observed before binning is
      performed.  Default value `100`.

 - `passes::Int`: Number of passes to take over the dataset.  Default
      value `1`.

 - `test::Tuple{Array{Bool, 1}, Array{Float64, 2}}`: Testing dataset (may
      be categorical).
 - `test_labels::Array{Int, 1}`: Labels of test data.
 - `training::Tuple{Array{Bool, 1}, Array{Float64, 2}}`: Training dataset
      (may be categorical).
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output_model::HoeffdingTreeModel`: Output for trained Hoeffding tree
      model.
 - `predictions::Array{Int, 1}`: Matrix to output label predictions for
      test data into.
 - `probabilities::Array{Float64, 2}`: In addition to predicting labels,
      provide rediction probabilities in this matrix.

"""
function hoeffding_tree(;
                        batch_mode::Union{Bool, Missing} = missing,
                        bins::Union{Int, Missing} = missing,
                        confidence::Union{Float64, Missing} = missing,
                        info_gain::Union{Bool, Missing} = missing,
                        input_model::Union{HoeffdingTreeModel, Missing} = missing,
                        labels = missing,
                        max_samples::Union{Int, Missing} = missing,
                        min_samples::Union{Int, Missing} = missing,
                        numeric_split_strategy::Union{String, Missing} = missing,
                        observations_before_binning::Union{Int, Missing} = missing,
                        passes::Union{Int, Missing} = missing,
                        test::Union{Tuple{Array{Bool, 1}, Array{Float64, 2}}, Missing} = missing,
                        test_labels = missing,
                        training::Union{Tuple{Array{Bool, 1}, Array{Float64, 2}}, Missing} = missing,
                        verbose::Union{Bool, Missing} = missing,
                        points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, hoeffding_treeLibrary), Nothing, ());

  IORestoreSettings("Hoeffding trees")

  # Process each input argument before calling mlpackMain().
  if !ismissing(batch_mode)
    IOSetParam("batch_mode", convert(Bool, batch_mode))
  end
  if !ismissing(bins)
    IOSetParam("bins", convert(Int, bins))
  end
  if !ismissing(confidence)
    IOSetParam("confidence", convert(Float64, confidence))
  end
  if !ismissing(info_gain)
    IOSetParam("info_gain", convert(Bool, info_gain))
  end
  if !ismissing(input_model)
    hoeffding_tree_internal.IOSetParamHoeffdingTreeModel("input_model", convert(HoeffdingTreeModel, input_model))
  end
  if !ismissing(labels)
    IOSetParamURow("labels", labels)
  end
  if !ismissing(max_samples)
    IOSetParam("max_samples", convert(Int, max_samples))
  end
  if !ismissing(min_samples)
    IOSetParam("min_samples", convert(Int, min_samples))
  end
  if !ismissing(numeric_split_strategy)
    IOSetParam("numeric_split_strategy", convert(String, numeric_split_strategy))
  end
  if !ismissing(observations_before_binning)
    IOSetParam("observations_before_binning", convert(Int, observations_before_binning))
  end
  if !ismissing(passes)
    IOSetParam("passes", convert(Int, passes))
  end
  if !ismissing(test)
    IOSetParam("test", convert(Tuple{Array{Bool, 1}, Array{Float64, 2}}, test), points_are_rows)
  end
  if !ismissing(test_labels)
    IOSetParamURow("test_labels", test_labels)
  end
  if !ismissing(training)
    IOSetParam("training", convert(Tuple{Array{Bool, 1}, Array{Float64, 2}}, training), points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output_model")
  IOSetPassed("predictions")
  IOSetPassed("probabilities")
  # Call the program.
  hoeffding_tree_mlpackMain()

  return hoeffding_tree_internal.IOGetParamHoeffdingTreeModel("output_model"),
         IOGetParamURow("predictions"),
         IOGetParamMat("probabilities", points_are_rows)
end


text:
is_constant(t::Term) = is_constant(get(t))
is_constant(x) = typeof(x) ∉ [Expr, Symbolic, Variable]

is_ground(t::Term) = is_ground(get(t))
is_ground(ex::Expr) = all(is_ground, ex.args)
is_ground(::Symbolic) = false
is_ground(::Variable) = false
is_ground(x) = true

vars(t::Term) = vars(get(t))
vars(x::Variable) = [x]
vars(t::Expr) = [vars.(t.args)...;]
vars(x) = []

_image(x, S) = σ -> isvalid(Image(σ[x], S))

ext: jl


text:
# GFlops for vector triad
GFlops(N)=N*2.0/1.0e9

# Create arrays for test
function make_arrays(N)
    a = Array{Float64,1}(undef,N)
    b = Array{Float64,1}(undef,N)
    c = Array{Float64,1}(undef,N)
    d = Array{Float64,1}(undef,N)

    Threads.@threads for i=1:N
        a[i]=i
        b[i]=N-i
        c[i]=i
        d[i]=-i
    end
    return(a,b,c,d)
end

# Create shared arrays for test
function make_shared_arrays(N)
    (_a,_b,_c,_d)=make_arrays(N)
    a = SharedArray(_a)
    b = SharedArray(_b)
    c = SharedArray(_c)
    d = SharedArray(_d)
    return(a,b,c,d)
end

max_stars_repo_name: j-fu/julia-tests

max_issues_repo_name: j-fu/julia-tests

max_forks_repo_name: j-fu/julia-tests


text:
using Laplacians
include("compareSolvers.jl")
compareSolversOut(1000,10)


text:
lines = readlines("inputs/6-10/08.txt")

t = 0
for line in lines
  lin = split(line, ' ')[end-4:end]
  foo = filter(x -> length(x) ∈ [2,3,4,7], lin )
  t += length(foo)
end
@show t


text:
using LinearAlgebra
using FinEtools
using PGFPlotsX

objects = []

E = 200*phun("GPa");
nu = 0.3;
G = E/2/(1+nu)
a = 1.0*phun("m");
t = a/5
D = E*t^3/12/(1-nu^2)
P = 1.0
a = 1.0
r = collect(t/100000:t/10000:t)
w_b = @. P/16/pi/D * (2*r^2 * log(r/a) + (3+nu)/(1+nu)*(a^2 - r^2))
w_s = @. -P/4/pi/G/t * log(r/a)
max_w_b = maximum(w_b)

@pgf p = PGFPlotsX.Plot(
{
color = "black",
line_width  = 0.7,
style = "solid"
},
Coordinates([v for v in  zip(r./t, w_b./max_w_b)])
)
push!(objects, p)
push!(objects, LegendEntry("Bending only"))

@pgf p = PGFPlotsX.Plot(
{
color = "red",
line_width  = 0.7,
style = "solid"
},
Coordinates([v for v in  zip(r./t, (w_b+w_s)./max_w_b)])
)
push!(objects, p)
push!(objects, LegendEntry("Bending plus shear"))

@pgf ax = Axis(
    {
        xlabel = "Normalized Distance [ND]",
        ylabel = "Normalized Displacement [ND]",
        # xmin = range[1],
        # xmax = range[2],
        xmode = "linear",
        ymode = "linear",
        yminorgrids = "true",
        grid = "both",
        legend_style = {
            at = Coordinate(0.5, 1.05),
            anchor = "south",
            legend_columns = -1
        },
    },
    objects...
)

display(ax)
pgfsave("deflections.pdf", ax)


text:
if lowercase(get(ENV, "CI", "false")) == "true"
    include("install_dependencies.jl")
end

using StateSpaceReconstruction
using Statistics: std
using LinearAlgebra
using Test

if lowercase(get(ENV, "CI", "false")) == "true"
    include("install_dependencies.jl")
end

include("Embeddings/test_Embeddings.jl")
include("Embeddings/test_invariantize.jl")

include("RectangularPartitions/test_RectangularPartitions.jl")
include("Triangulations/test_Simplices.jl")
include("Triangulations/test_Delaunay.jl")
include("Triangulations/test_Triangulations.jl")


text:
# MvNormalCDF
# Copyright © 2019-2021 Vladimir Arnautov aka PharmCat <mail@pharmcat.net>, Andrew Gough
using MvNormalCDF
using Test, Distributions, StableRNGs

td = Array{Any}(undef,(14,9))
#  1-cov mtx 2-a 3-b 4-m 5-p 6-ptol 7-e 8-etol

# from MATLAB documenation 4 dim
td[1,1] = [4 3 2 1;3 5 -1 1;2 -1 4 2;1 1 2 5]  # Σ cov Matrix
td[1,2] = [-Inf; -Inf; -Inf; -Inf]             # a lower integration limit
td[1,3] = [1; 2; 3; 4 ]                        # b upper integration limit
td[1,4] = 5000
td[1,5] = 0.605653                            	# expected p value
td[1,6] = 0.001557374                          	# ± p tolerance
td[1,7] = 0.001394971                           # expected e (error) value
td[1,8] = 0.0009277058                        	# ± e tolerance

# 3 dim
td[2,1] =[1  3/5  1/3; 3/5 1  11/15; 1/3 11/15 1]
td[2,2] = [-Inf;-Inf;-Inf]
td[2,3] = [1;4;2]
td[2,4] = 3000
td[2,5] = 0.827985
td[2,6] = 3.529068e-5
td[2,7] = 2.8608229e-5
td[2,8] = 3.72604931e-5

# 3 dim
td[3,1] = [1 0.25 0.2; 0.25 1 0.333333333; 0.2 0.333333333 1]
td[3,2] = [-1;-4;-2]
td[3,3] = [1;4;2]
td[3,4] = 4000
td[3,5] = 0.65368
td[3,6] = 0.000002089699
td[3,7] = 0.000001799225
td[3,8] = 0.000001161755

# Genz book eq. 1.5 p. 4-5 & p. 63
# Genz gives wrong answer on p. 4-5
td[4,1] = [1/3 3/5 1/3; 3/5 1.0 11/15; 1/3 11/15 1.0]
td[4,2] = [-Inf; -Inf; -Inf]
td[4,3] = [1; 4; 2]
td[4,4] = 4000
td[4,5] = 0.943174
td[4,6] = 0.00006724871
td[4,7] = 0.0000509768
td[4,8] = 0.00003337435

# Genz p. 63 uses different r matrix
td[5,1] = [1 0 0; 3/5 1 0; 1/3 11/15 1]
td[5,2] = [-Inf; -Inf; -Inf]
td[5,3] = [1; 4; 2]
td[5,4] = 4000
td[5,5] = 0.827985
td[5,6] = 0.00001024322
td[5,7] = 0.000008933135
td[5,8] = 0.000005194712

# singular example
# problem reduces to univariate problem with
# p = cdf.(Normal(),1) = 0.84134476068543
td[6,1] = [1 1 1; 1 1 1; 1 1 1]
td[6,2] = [-Inf, -Inf, -Inf]
td[6,3] = [1, 1, 1]
td[6,4] = 3000
td[6,5] = 0.841345
td[6,6] = 4.440892E-16
td[6,7] = 0.0
td[6,8] = 2.220446E-16

# 5 dim example
td[7,1] = [1 1 1 1 1;
		 1 2 2 2 2;
		 1 2 3 3 3;
		 1 2 3 4 4;
		 1 2 3 4 5]
td[7,2] = [-1,-2,-3,-4,-5]
td[7,3] = [2,3,4,5,6]
td[7,4] = 6000
td[7,5] = 0.761243
td[7,6] = 0.0004715559
td[7,7] = 0.0003493157
td[7,8] = 0.0002319606

# Genz used wrong integration limits when computing
# above. see p. 63
td[8,1] = td[7,1]
td[8,2] = sort(td[7,2])
td[8,3] = 1 .- td[8,2]
td[8,4] = 6000
td[8,5] = 0.474128
td[8,6] = 0.000037335
td[8,7] = 0.00003169846
td[8,8] = 0.0000178166

# positive orthant probability of above
td[9,1] = td[7,1]
td[9,2] = [0,0,0,0,0]
td[9,3] = td[8,3]
td[9,4] = 6000
td[9,5] = 0.113537
td[9,6] = 0.0001740044
td[9,7] = 0.0001584365
td[9,8] = 0.0001053655

# test 7 Cov matrix, but now -Inf lower limit
td[10,1] = td[7,1]
td[10,2] = [-Inf,-Inf,-Inf,-Inf,-Inf]
td[10,3] = td[8,3]
td[10,4] = 6000
td[10,5] = 0.810315
td[10,6] = 0.00003136998
td[10,7] = 0.0000287006
td[10,8] = 0.00001567681

# eight dimensional test
td[11,1] = [1 1 1 1 1 1 1 1;
		 	1 2 2 2 2 2 2 2;
		 	1 2 3 3 3 3 3 3;
		 	1 2 3 4 4 4 4 4;
		 	1 2 3 4 5 5 5 5;
		 	1 2 3 4 5 6 6 6;
		 	1 2 3 4 5 6 7 7;
		 	1 2 3 4 5 6 7 8]
td[11,2] =  -1*[1,2,3,4,5,6,7,8]
td[11,3] = [2,3,4,5,6,7,8,9]
td[11,4] = 9000
td[11,5] = 0.759474
td[11,6] = 0.0006248131
td[11,7] = 0.0005161476
td[11,8] = 0.0003407714

# orthant probability of above
td[12,1] = td[11,1]
td[12,2] = [0,0,0,0,0,0,0,0]
td[12,3] = [Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf]
td[12,4] = 9000
td[12,5] = 0.19638
td[12,6] = 0.0002654318
td[12,7] = 0.0002384606
td[12,8] = 0.0001693536

# with dim with lower limit -Inf (exact result)
td[13,1] = td[11,1]
td[13,2] = -Inf*[1,1,1,1,1,1,1,1]
td[13,3] = td[12,3]
td[13,4] = 9000
td[13,5] = 1.0
td[13,6] = 2.220446E-16
td[13,7] = 0.0
td[13,8] = 2.220446E-16

# 25 dimensions
td[14,1] = [59.227 2.601 3.38 8.303 -0.334 11.029 10.908 0.739 4.703 7.075 8.049 1.403 9.838 5.46 11.949 2.272 7.234 15.215 -9.091 12.265 3.01 -3.199 10.608 8.464 -8.685;
 2.601 77.213 0.882 2.99 -2.536 -4.55 -3.874 -3.607 6.023 3.129 15.7 -7.271 -10.655 8.456 15.387 -5.764 6.617 -6.331 -2.244 -0.925 -7.516 13.836 9.243 -0.84 -3.781;
 3.38 0.882 79.72 -1.465 3.179 -1.799 9.842 9.165 -1.54 -8.03 0.778 9.053 -2.598 -8.844 15.857 13.613 -1.878 8.18 11.806 -0.242 4.711 2.258 9.554 1.184 -8.047;
 8.303 2.99 -1.465 65.041 3.803 2.13 -0.936 -5.996 2.719 -4.648 4.611 4.486 13.38 -0.376 0.179 14.654 -7.089 -1.194 9.357 5.12 4.943 -0.475 4.764 -8.56 2.337;
 -0.334 -2.536 3.179 3.803 64.128 4.85 -14.767 -10.044 12.437 5.065 8.191 0.391 0.068 9.754 -0.062 0.429 9.265 5.502 4.227 0.559 0.811 3.169 4.558 -1.878 -3.885;
 11.029 -4.55 -1.799 2.13 4.85 64.912 12.989 2.675 1.227 8.205 3.3 3.545 6.225 11.936 -2.956 6.188 -2.206 1.184 10.546 2.492 12.035 -15.789 4.296 1.086 7.225;
 10.908 -3.874 9.842 -0.936 -14.767 12.989 78.424 2.03 5.595 4.921 1.573 -7.355 -9.425 -4.024 4.912 -12.05 -2.034 -2.435 -10.355 0.985 8.23 2.806 6.254 -1.494 5.531;
 0.739 -3.607 9.165 -5.996 -10.044 2.675 2.03 93.184 8.84 3.843 -14.968 16.386 -0.223 4.398 -4.786 1.731 4.025 0.479 3.12 -15.591 12.614 -8.279 -3.582 6.597 -1.915;
 4.703 6.023 -1.54 2.719 12.437 1.227 5.595 8.84 73.981 4.019 -6.404 5.869 -4.305 5.936 2.016 5.987 10.051 -0.705 13.229 -1.715 7.102 12.89 10.967 5.262 15.954;
 7.075 3.129 -8.03 -4.648 5.065 8.205 4.921 3.843 4.019 68.048 7.346 7.412 9.956 6.743 2.547 0.177 0.844 2.147 -4.072 11.832 -3.55 -0.096 -1.96 -1.381 -3.249;
 8.049 15.7 0.778 4.611 8.191 3.3 1.573 -14.968 -6.404 7.346 66.708 -6.22 -8.952 0.647 3.039 -12.078 7.618 10.398 -5.055 1.577 -12.77 11.477 8.272 2.071 7.728;
 1.403 -7.271 9.053 4.486 0.391 3.545 -7.355 16.386 5.869 7.412 -6.22 55.025 0.273 -12.049 -2.71 11.531 12.604 1.294 2.791 -5.698 -2.231 15.025 12.229 -5.876 -3.374;
 9.838 -10.655 -2.598 13.38 0.068 6.225 -9.425 -0.223 -4.305 9.956 -8.952 0.273 82.184 3.466 -9.297 -4.347 12.586 4.372 13.705 -6.795 -5.818 -6.78 5.11 10.099 5.05;
 5.46 8.456 -8.844 -0.376 9.754 11.936 -4.024 4.398 5.936 6.743 0.647 -12.049 3.466 69.352 10.855 7.282 6.615 4.58 0.306 6.482 14.589 5.081 -9.141 -4.657 10.763;
 11.949 15.387 15.857 0.179 -0.062 -2.956 4.912 -4.786 2.016 2.547 3.039 -2.71 -9.297 10.855 87.988 2.504 4.226 10.461 9.703 -3.112 -13.348 0.944 -2.824 -4.498 10.551;
 2.272 -5.764 13.613 14.654 0.429 6.188 -12.05 1.731 5.987 0.177 -12.078 11.531 -4.347 7.282 2.504 63.611 3.024 2.35 -5.191 -6.101 -6.324 -0.483 9.899 5.768 2.382;
 7.234 6.617 -1.878 -7.089 9.265 -2.206 -2.034 4.025 10.051 0.844 7.618 12.604 12.586 6.615 4.226 3.024 78.3 1.54 -6.868 2.613 6.006 5.49 9.06 -4.229 -4.395;
 15.215 -6.331 8.18 -1.194 5.502 1.184 -2.435 0.479 -0.705 2.147 10.398 1.294 4.372 4.58 10.461 2.35 1.54 61.499 11.083 15.428 1.771 2.517 5.181 13.476 6.829;
 -9.091 -2.244 11.806 9.357 4.227 10.546 -10.355 3.12 13.229 -4.072 -5.055 2.791 13.705 0.306 9.703 -5.191 -6.868 11.083 76.538 4.287 6.564 -2.49 15.558 11.202 -16.964;
 12.265 -0.925 -0.242 5.12 0.559 2.492 0.985 -15.591 -1.715 11.832 1.577 -5.698 -6.795 6.482 -3.112 -6.101 2.613 15.428 4.287 59.057 -5.619 4.374 8.934 -6.203 14.687;
 3.01 -7.516 4.711 4.943 0.811 12.035 8.23 12.614 7.102 -3.55 -12.77 -2.231 -5.818 14.589 -13.348 -6.324 6.006 1.771 6.564 -5.619 69.833 -1.233 -7.922 -7.027 -2.315;
 -3.199 13.836 2.258 -0.475 3.169 -15.789 2.806 -8.279 12.89 -0.096 11.477 15.025 -6.78 5.081 0.944 -0.483 5.49 2.517 -2.49 4.374 -1.233 65.458 5.786 4.613 2.223;
 10.608 9.243 9.554 4.764 4.558 4.296 6.254 -3.582 10.967 -1.96 8.272 12.229 5.11 -9.141 -2.824 9.899 9.06 5.181 15.558 8.934 -7.922 5.786 69.574 3.527 -2.802;
 8.464 -0.84 1.184 -8.56 -1.878 1.086 -1.494 6.597 5.262 -1.381 2.071 -5.876 10.099 -4.657 -4.498 5.768 -4.229 13.476 11.202 -6.203 -7.027 4.613 3.527 65.98 3.495;
 -8.685 -3.781 -8.047 2.337 -3.885 7.225 5.531 -1.915 15.954 -3.249 7.728 -3.374 5.05 10.763 10.551 2.382 -4.395 6.829 -16.964 14.687 -2.315 2.223 -2.802 3.495 58.327]
 td[14,2] = vec(-Inf*fill(1,(25,1)))
 td[14,3] = [6.0; 9.0; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf]
 td[14,4] = 35000
 td[14,5] = 0.665342
 td[14,6] = 0.000001975944
 td[14,7] = 0.00000181746
 td[14,8] = 0.000001083068

 td[1,9] = 0.6054064156442294
 td[2,9] = 0.8279766144831429
 td[3,9] = 0.6536793504453967
 td[4,9] = 0.9431576518489296
 td[5,9] = 0.8279829100356028
 td[6,9] = 0.8413447460685427
 td[7,9] = 0.7612423324911164
 td[8,9] = 0.4741170481094123
 td[9,9] = 0.11351703757272677
 td[10,9] = 0.8103089978617527
 td[11,9] = 0.7593693258920418
 td[12,9] = 0.19640599118630003
 td[13,9] = 1.0
 td[14,9] = 0.6653426686040154

 @testset "MvNormalCDF test" begin
	for i in 1:14
	    r = td[i,1]
	    a = td[i,2]
	    b = td[i,3]
	    m = td[i,4]
	    pexpected = td[i,5]
	    eexpected = td[i,7]
	    ptol = td[i,6]
	    etol = td[i,8]

	    (p,e) = MvNormalCDF.mvnormcdf(r, a, b; m=m, rng = StableRNG(1234))
		v = p
		p = round(p, digits=6)

	    @test p ≈ pexpected atol=ptol
	    @test e ≈ eexpected atol=etol

		@test v ≈ td[i, 9] atol=1e-10
	end

	#just test not validation
	@test_nowarn MvNormalCDF.mvnormcdf(td[3,1], [0, 0, 0], [Inf, Inf, Inf]; m=m, rng = StableRNG(1234))

    # test warning on singular Σ
    r = td[6,1]
    a = td[6,2]
    b = td[6,3]
    m = td[6,4]

    #@test_logs (:warn,"covariance matrix Σ fails positive definite check") A.qsimvnv(r,a,b;m=m)

    # Σ dimension 1 throws error
    r = Array{Float64}(undef,(1,1))
    r[1,1] = 5
    @test_throws ErrorException MvNormalCDF.mvnormcdf(MvNormal(r),a,b)

    # a < b throws error
    r = td[3,1]
    a = [-Inf, 3, 0]
    b = [0, 0, 1]
    @test_throws ArgumentError MvNormalCDF.mvnormcdf(MvNormal(r),a,b)

    # non-square Σ throws error
    r = [1 2; 3 4; 5 6]
    a = [-Inf, -Inf, -Inf]
    b = [0, 0, 0]
    @test_throws DimensionMismatch MvNormalCDF.mvnormcdf(MvNormal(r),a,b)

	# mismatched dimensions throws errors
	r = [1 2; 3 4; 5 6]
    a = [-Inf, -Inf, -Inf]
	b = [0, 0]

	@test_throws DimensionMismatch MvNormalCDF.mvnormcdf(MvNormal(r),a,b)

 end

#=
using BenchmarkTools
mvn = MvNormal(td[14,1])
a = td[14,2]
b = td[14,3]
MvNormalCDF.mvnormcdf(mvn,a,b)
@benchmark  MvNormalCDF.mvnormcdf($mvn, $a, $b)

BenchmarkTools.Trial: 262 samples with 1 evaluation.
 Range (min … max):  18.809 ms …  29.450 ms  ┊ GC (min … max): 0.00% … 35.47%
 Time  (median):     19.086 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   19.140 ms ± 676.836 μs  ┊ GC (mean ± σ):  0.21% ±  2.19%

               ▁  ▁▂ ▁▃▃ ▂▂▂▄▁▂██▃▃▃▁▂▃ ▄  ▄   ▁
  ▃▁▃▁▆▁▁▃█▃▄▃▄█▆▇██▆███▆██████████████▇██▇█▇▇▇█▆█▄▆▆▄▄▃▃▄▁▃▃▃ ▄
  18.8 ms         Histogram: frequency by time         19.4 ms <

 Memory estimate: 892.28 KiB, allocs estimate: 2334.
 =#

#=
using BenchmarkTools
Σ = [4 3 2 1; 3 5 -1 1; 2 -1 4 2; 1 1 2 5]
μ = zeros(4)
a = [-Inf; -Inf; -Inf; -Inf]
b = [1; 2; 3; 4]
m = 5000
mvn = MvNormal(Σ)
@benchmark  MvNormalCDF.mvnormcdf($mvn, $a, $b)

BenchmarkTools.Trial: 7939 samples with 1 evaluation.
 Range (min … max):  602.700 μs …  18.335 ms  ┊ GC (min … max): 0.00% … 96.37%
 Time  (median):     621.000 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   627.835 μs ± 199.151 μs  ┊ GC (mean ± σ):  0.35% ±  1.08%

              █▇▁
  ▂▁▁▂▁▂▂▂▂▃▅▇███▆▄▄▄▄▅▄▄▄▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
  603 μs           Histogram: frequency by time          682 μs <

 Memory estimate: 38.27 KiB, allocs estimate: 98.

julia> @time MvNormalCDF.mvnormcdf(mvn, a, b)
  0.000663 seconds (98 allocations: 38.266 KiB)
(0.6055881147694031, 0.0009341283457742374)
@time MvNormalCDF.mvnormcdf(Σ, a, b)
  0.000761 seconds (98 allocations: 38.266 KiB)
(0.6064777984403792, 0.0012501174586413868)
julia> @time MvNormalCDF.mvnormcdf(μ, Σ, a, b)
  0.000668 seconds (98 allocations: 38.266 KiB)
(0.6057553050182675, 0.0019361170915866777)

=#

max_stars_repo_name: PharmCat/MvNormalCDF

max_issues_repo_name: PharmCat/MvNormalCDF

max_forks_repo_name: PharmCat/MvNormalCDF


text:
@testset "kron" begin
    m1 = GBMatrix(UInt64[1, 2, 3, 5], UInt64[1, 3, 1, 2], Int8[1, 2, 3, 5])
    n1 = GBMatrix(ones(UInt32, 4, 4))
    m2 = sparse([1, 2, 3, 5], [1, 3, 1, 2], Int8[1, 2, 3, 5])
    n2 = ones(Int32, 4, 4)
    o1 = kron(m1, n1)
    @test o1 == GBMatrix(kron(m2, n2)) #basic kron is equivalent
    mask = GBMatrix{Bool}(20, 12)
    mask[17:20, 5:8] = false #don't care value, using structural
    #mask out bottom chunk using structural complement
    o2 = kron(m1, n1; mask, desc=Descriptor(structural_mask=true, complement_mask=true))
    @test o2[20, 5] === nothing #We don't want values in masked out area
    @test o2[1:2:15, :] == o1[1:2:15, :] #The rest should match, test indexing too.
end


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

ToolInfo holds a <toolspecific> tag.

It wraps a vector of well formed elements parsed into [`AnyElement`](@ref)s.
for use by anything that understands toolname, version toolspecifics.
"""
@auto_hash_equals struct ToolInfo
    toolname::String
    version::String
    infos::Vector{AnyElement} #TODO specialize infos.
    xml::XMLNode
end

name(ti::ToolInfo) = ti.toolname
version(ti::ToolInfo) = ti.version
infos(ti::ToolInfo) = ti.infos

xmlnode(ti::ToolInfo) = ti.xml

###############################################################################

"""
$(TYPEDEF)
$(TYPEDFIELDS)

TokenGraphics is <toolspecific> content.
Combines the <tokengraphics> and <tokenposition> elements.
"""
struct TokenGraphics <: AbstractPnmlTool
    positions::Vector{Coordinate} #TODO: uses abstract type
end

# Empty TokenGraphics is allowed in spec.
TokenGraphics() = TokenGraphics(Coordinate[])


text:
import Base.Random: rand

"""
```julia
AbstractRNG{T} <: Base.Random.AbstractRNG
```

The abstract type of Random Number Generators. T indicates the original output type of a RNG.
"""
abstract AbstractRNG{T<:Number} <: Base.Random.AbstractRNG

typealias BitTypes Union{Bool, Signed, Unsigned}

# see https://github.com/sunoru/RNG.jl/issues/8
# TODO: find a better approach.
@inline function rand{T<:Union{UInt64, UInt128}}(rng::AbstractRNG{T}, ::Type{Float64}=Float64)
    u = rand(rng, UInt64)
    reinterpret(Float64, Base.exponent_one(Float64) | Base.significand_mask(Float64) & u) - 1.0
end
@inline function rand{T<:Union{UInt8, UInt16, UInt32}}(rng::AbstractRNG{T}, ::Type{Float64}=Float64)
    rand(rng, T) * exp2(-sizeof(T) << 3)
end

@inline function rand{T1<:BitTypes, T2<:BitTypes}(rng::AbstractRNG{T1}, ::Type{T2})
    s1 = sizeof(T1)
    s2 = sizeof(T2)
    t = rand(rng, T1) % T2
    s1 > s2 && return t
    for i in 2:(s2 ÷ s1)
        t |= (rand(rng, T1) % T2) << ((s1 << 3) * (i - 1))
    end
    t
end


text:
function get_rounded_zeros_and_ones!(v::Array{Float64}, tol_zero::Float64)

    for i in findall(abs.(v) .<= tol_zero)
        v[i] = 0
    end

    for i in findall(((1-tol_zero) .<= abs.(v) .<= (1+tol_zero)))
        if v[i] > 0
            v[i] = 1
        elseif v[i] < 0
            v[i] = -1
        end
    end

end

"""
    optimal_graph_edges(adjacency_matrix::Array{Float64})

Returns a vector of tuples of edges corresponding
to an input adjacency matrix of the graph.
"""
function optimal_graph_edges(adjacency_matrix::Array{Float64})

    edges = Vector{Tuple{Int64, Int64}}()

    num_nodes = size(adjacency_matrix)[1]

    if !LA.issymmetric(adjacency_matrix)
        Memento.error(_LOGGER, "Input adjacency matrix is asymmetric")
    end

    for i=1:num_nodes
        if !isapprox(adjacency_matrix[i,i], 0)
            Memento.error(_LOGGER, "Input adjacency matrix cannot have self loops")
        end
    end

    for i=1:(num_nodes-1)
        for j=(i+1):num_nodes
            if !isapprox(adjacency_matrix[i,j], 0)
                push!(edges, (i,j))
            end
        end
    end

    return edges
end

"""
    laplacian_matrix(adjacency_matrix::Array{Float64})

Returns the weighted Laplacian matrix for an input weighted adjacency matrix of the graph.
"""
function laplacian_matrix(adjacency_matrix::Array{Float64})

    if !LA.issymmetric(adjacency_matrix)
        Memento.error(_LOGGER, "Input adjacency matrix is asymmetric")
    end

    num_nodes = size(adjacency_matrix)[1]

    laplacian_matrix = zeros(Float64, num_nodes, num_nodes)

    for i = 1:num_nodes
        for j = i:num_nodes

            if i == j
                if !isapprox(adjacency_matrix[i,j], 0)
                    Memento.error(_LOGGER, "Input adjacency matrix cannot have self loops")
                end

                laplacian_matrix[i,j] = sum(adjacency_matrix[i,:])
            else
                if adjacency_matrix[i,j] <= -1E-6
                    Memento.error(_LOGGER, "Input adjacency matrix cannot have negative weights")
                end

                if !isapprox(adjacency_matrix[i,j], 0)
                    laplacian_matrix[i,j] = -adjacency_matrix[i,j]
                    laplacian_matrix[j,i] = -adjacency_matrix[i,j]
                end
            end

        end
    end

    return laplacian_matrix
end

"""
    fiedler_vector(adjacency_matrix::Array{Float64})

Returns the Fiedler vector or the eigenvector corresponding to the
second smallest eigenvalue of the Laplacian matrix for an input weighted adjacency matrix of the graph.
"""
function fiedler_vector(adjacency_matrix::Array{Float64})

    L_mat = LOpt.laplacian_matrix(adjacency_matrix)
    ac    = LOpt.algebraic_connectivity(adjacency_matrix)

    fiedler = LA.eigvecs(L_mat)[:,2]

    if !isapprox(fiedler' * L_mat * fiedler, ac, atol=1E-6)
        Memento.error(_LOGGER, "Evaluated eigenvector is not the Fiedler vector")
    end

    return fiedler

end

"""
    algebraic_connectivity(adjacency_matrix::Array{Float64})

Returns the algebraic connectivity or the
second smallest eigenvalue of the Laplacian matrix, for an input weighted adjacency matrix of the graph.
"""
function algebraic_connectivity(adjacency_matrix::Array{Float64})

    L_mat = LOpt.laplacian_matrix(adjacency_matrix)

    ac = sort(LA.eigvals(L_mat))[2]

    return ac
end


text:
using Distributed
@everywhere import Pkg
@everywhere Pkg.offline(true)  # Tell Julia package manager not to download an updated registry for
@everywhere Pkg.activate(".")
using NBInclude
println("# About to run Jupyter notebook")
flush(stdout)
@nbinclude("ex2.ipynb")
println("# Completed running notebook")
flush(stdout)


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.
module Lenslets
using LinearAlgebra
import Unitful
using Unitful:uconvert,ustrip
using Unitful.DefaultSymbols:mm,μm,nm #Unitful.DefaultSymbols exports a variable T which can cause major confusion with type declarations since T is a common parametric type symbol. Import only what is needed.
using Colors
using StaticArrays
import DataFrames
import ..Repeat
import SpecialFunctions
import Plots
import ...OpticSim
import ...OpticSim.Geometry

include("HexClusters.jl")
include("HexTilings.jl")
include("Analysis.jl")
include("DisplayGeneration.jl")
include("LensletAssignment.jl")

end # module
export Lenslets


text:
# TODO: These types aren't used anywhere for multiple dispatch, so why have them?
abstract type LyapunovMatrixEquationsMaps{T} <: LinearMaps.LinearMap{T} end
abstract type SylvesterMatrixEquationsMaps{T} <: LinearMaps.LinearMap{T} end
const MatrixEquationsMaps{T} = Union{LyapunovMatrixEquationsMaps{T},SylvesterMatrixEquationsMaps{T}}

abstract type ContinuousOrDiscrete end
struct Continuous <: ContinuousOrDiscrete end
struct Discrete <: ContinuousOrDiscrete end

"""
    M = trmatop(n, m)

Define the transposition operator `M: X -> X'` for all `n x m` matrices.
"""
struct trmatop{Int} <: LinearMaps.LinearMap{Int}
   size::Dims{2}
   function trmatop(dims::Dims{2})
      all(≥(0), dims) || throw(ArgumentError("dims must be non-negative"))
      return new{Int}(dims)
   end
   function trmatop(m::Int, n::Int)
      (m ≥ (0) & n ≥ (0))  || throw(ArgumentError("dimensions must be non-negative"))
      return new{Int}((m, n))
   end
end
trmatop(n::Int) = trmatop(n, n)
Base.size(A::trmatop) = (prod(A.size), prod(A.size))
LinearAlgebra.issymmetric(A::trmatop) = A.size[1] == A.size[2]
LinearAlgebra.ishermitian(A::trmatop) = A.size[1] == A.size[2]
function mul!(y::AbstractVector, A::trmatop, x::AbstractVector)
   X = reshape(x, A.size...)
   LinearMaps.check_dim_mul(y, A, x)
   copyto!(y, adjoint(X))
   return y
end

"""
    M = trmatop(A)

Define the transposition operator `M: X -> X'` of all matrices of the size of `A`.
"""
trmatop(A) = trmatop(size(A))

struct LyapunovMap{T,TA <: AbstractMatrix,CD <: ContinuousOrDiscrete} <: LyapunovMatrixEquationsMaps{T}
   A::TA
   her::Bool
   function LyapunovMap{T,TA,CD}(A::TA, her=false) where {T,TA<:AbstractMatrix{T},CD}
      LinearAlgebra.checksquare(A)
      return new{T,TA,CD}(A, her)
   end
end

LyapunovMap(A::TA, ::CD = Continuous(), her::Bool = false) where {T,TA<:AbstractMatrix{T},CD<:ContinuousOrDiscrete} =
   LyapunovMap{T,TA,CD}(A, her)

LinearAlgebra.adjoint(L::LyapunovMap{<:Any,<:Any,CD}) where {CD}   = LyapunovMap(L.A', CD(), L.her)
LinearAlgebra.transpose(L::LyapunovMap{<:Any,<:Any,CD}) where {CD} = LyapunovMap(L.A', CD(), L.her)

"""
    L = lyapop(A; disc = false, her = false)

Define, for an `n x n` matrix `A`, the continuous Lyapunov operator `L:X -> AX+XA'`
if `disc = false` or the discrete Lyapunov operator `L:X -> AXA'-X` if `disc = true`.
If `her = false` the Lyapunov operator `L:X -> Y` maps general square matrices `X`
into general square matrices `Y`, and the associated matrix `M = Matrix(L)` is
``n^2 \\times n^2``.
If `her = true` the Lyapunov operator `L:X -> Y` maps symmetric/Hermitian matrices `X`
into symmetric/Hermitian matrices `Y`, and the associated matrix `M = Matrix(L)` is
``n(n+1)/2 \\times n(n+1)/2``.
For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov
operators. Linear Algebra and its Applications 312:35–71, 2000.
"""
lyapop(A::AbstractMatrix; disc=false, her=false) =
   LyapunovMap(A, ifelse(disc, Discrete(), Continuous()), her)
lyapop(A::Schur; disc=false, her=false) =
   LyapunovMap(A.T, ifelse(disc, Discrete(), Continuous()), her)
lyapop(A::Number; disc=false, her=false) =
   LyapunovMap(fill(A, 1, 1), ifelse(disc, Discrete(), Continuous()), her)
function Base.size(L::LyapunovMap)
   n = size(L.A, 1)
   N = L.her ? Int(n * (n + 1) / 2) : n * n
   return (N, N)
end
function mul!(y::AbstractVector, L::LyapunovMap{T,TA,Discrete}, x::AbstractVector) where {T,TA}
   n = size(L.A, 1)
   T1 = promote_type(T, eltype(x))
   if L.her
      X = vec2triu(convert(AbstractVector{T1}, x), her=true)
      # (y .= triu2vec(utqu(X,L.A') - X))
      muldsym!(y, L.A, X)
   else
      X = reshape(convert(AbstractVector{T1}, x), n, n)
      # (y .= (L.A*X*L.A' - X)[:])
      mul!(y, -1, x)
      Y = reshape(y, n, n)
      mul!(Y, L.A, X * L.A', true, true)
   end
   return y
end
function mul!(y::AbstractVector, L::LyapunovMap{T,TA,Continuous}, x::AbstractVector) where {T,TA}
   n = size(L.A, 1)
   T1 = promote_type(T, eltype(x))
   if L.her
      X = vec2triu(convert(AbstractVector{T1}, x), her=true)
      mulcsym!(y, L.A, X)
   else
      X = reshape(convert(AbstractVector{T1}, x), n, n)
      # (y[:] = (L.A*X + X*L.A')[:])
      Y = reshape(y, n, n)
      mul!(Y, X, L.A')
      mul!(Y, L.A, X, true, true)
   end
   return y
end
function mulcsym!(y::AbstractVector, A::AbstractMatrix, X::AbstractMatrix)
   require_one_based_indexing(y, A, X)
   # A*X + X*A'
   n = size(A, 1)
   ZERO = zero(eltype(y))
   @inbounds begin
      k = 1
      for j = 1:n
         for i = 1:j
            temp = ZERO
            for l = 1:n
               temp += A[i,l] * X[l,j] + X[i,l] * conj(A[j,l])
            end
            y[k] = temp
            k += 1
         end
      end
   end
   return y
end
function muldsym!(y::AbstractVector, A::AbstractMatrix, X::AbstractMatrix)
   require_one_based_indexing(y, X)
   # A*X*A' - X
   n = size(A, 1)
   # t = triu(X)-diag(X)/2
   t = UpperTriangular(X) - Diagonal(X[diagind(X)] ./ 2)
   Y = similar(X, n, n)
   # Y = A*t*A'
   mul!(Y, A * t, A')
   # Y + Y' - X
   @inbounds begin
      k = 1
      for j = 1:n
         for i = 1:j
            y[k] = Y[i,j] + conj(Y[j,i]) - X[i,j]
            k += 1
         end
      end
   end
   return y
end

struct GeneralizedLyapunovMap{T,TA <: AbstractMatrix{T},TE <: AbstractMatrix{T},CD<:ContinuousOrDiscrete} <: LyapunovMatrixEquationsMaps{T}
   A::TA
   E::TE
   her::Bool
   function GeneralizedLyapunovMap{T,TA,TE,CD}(A::TA, E::TE, her=false) where {T,TA<:AbstractMatrix{T},TE<:AbstractMatrix{T},CD}
      n = LinearAlgebra.checksquare(A)
      n == LinearAlgebra.checksquare(E) ||
               throw(DimensionMismatch("E must be a square matrix of dimension $n"))
      return new{T,TA,TE,CD}(A, E, her)
   end
end
function GeneralizedLyapunovMap(A::AbstractMatrix, E::AbstractMatrix, ::CD = Continuous(), her::Bool = false) where {CD<:ContinuousOrDiscrete}
   n = LinearAlgebra.checksquare(A)
   n == LinearAlgebra.checksquare(E) ||
            throw(DimensionMismatch("E must be a square matrix of dimension $n"))
   T = promote_type(eltype(A), eltype(E))
   A = convert(AbstractMatrix{T}, A)
   E = convert(AbstractMatrix{T}, E)
   return GeneralizedLyapunovMap{T,typeof(A),typeof(E),CD}(A, E, her)
end

LinearAlgebra.adjoint(L::GeneralizedLyapunovMap{<:Any,<:Any,<:Any,CD}) where {CD} =
   GeneralizedLyapunovMap(L.A', L.E', CD(), L.her)
LinearAlgebra.transpose(L::GeneralizedLyapunovMap{<:Any,<:Any,<:Any,CD}) where {CD} =
   GeneralizedLyapunovMap(L.A', L.E', CD(), L.her)

"""
    L = lyapop(A, E; disc = false, her = false)

Define, for a pair `(A,E)` of `n x n` matrices, the continuous Lyapunov operator `L:X -> AXE'+EXA'`
if `disc = false` or the discrete Lyapunov operator `L:X -> AXA'-EXE'` if `disc = true`.
If `her = false` the Lyapunov operator `L:X -> Y` maps general square matrices `X`
into general square matrices `Y`, and the associated matrix `M = Matrix(L)` is
``n^2 \\times n^2``.
If `her = true` the Lyapunov operator `L:X -> Y` maps symmetric/Hermitian matrices `X`
into symmetric/Hermitian matrices `Y`, and the associated `M = Matrix(L)` is a
``n(n+1)/2 \\times n(n+1)/2``.
For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov
operators. Linear Algebra and its Applications 312:35–71, 2000.
"""
lyapop(A::AbstractMatrix, E::AbstractMatrix; disc=false, her=false) =
   GeneralizedLyapunovMap(A, E, ifelse(disc, Discrete(), Continuous()), her)
lyapop(F::GeneralizedSchur; disc=false, her=false) =
   GeneralizedLyapunovMap(F.S, F.T, ifelse(disc, Discrete(), Continuous()), her)
function lyapop(A::Number, E::Number; disc=false, her=false)
   A, E = promote(A, E)
   GeneralizedLyapunovMap(fill(A, 1, 1), fill(E, 1, 1), ifelse(disc, Discrete(), Continuous()), her)
end
function Base.size(L::GeneralizedLyapunovMap)
   n = size(L.A, 1)
   N = L.her ? Int(n * (n + 1) / 2) : n * n
   return (N, N)
end
function mul!(y::AbstractVector, L::GeneralizedLyapunovMap{T,<:Any,<:Any,Discrete}, x::AbstractVector) where {T}
   n = size(L.A, 1)
   T1 = promote_type(T, eltype(x))
   if L.her
      X = vec2triu(convert(AbstractVector{T1}, x), her=true)
      # (y .= triu2vec(utqu(X,L.A') - utqu(X,L.E')))
      muldsym!(y, L.A, L.E, X)
   else
      X = reshape(convert(AbstractVector{T1}, x), n, n)
      Y = reshape(y, n, n)
      temp = similar(Y, (n, n))
      # (y .= (L.A*X*L.A' - L.E*X*L.E')[:])
      # Y .= X
      mul!(temp, X, L.A')
      mul!(Y, L.A, temp)
      mul!(temp, X, L.E')
      mul!(Y, L.E, temp, -1, 1)
   end
   return y
end
function mul!(y::AbstractVector, L::GeneralizedLyapunovMap{T,<:Any,<:Any,Continuous}, x::AbstractVector) where {T}
   n = size(L.A, 1)
   T1 = promote_type(T, eltype(x))
   if L.her
      X = vec2triu(convert(AbstractVector{T1}, x), her=true)
      mulcsym!(y, L.A, L.E, X)
   else
      X = reshape(convert(AbstractVector{T1}, x), n, n)
      Y = reshape(y, n, n)
      temp = similar(Y, (n, n))
      # (y[:] = (L.A*X*L.E' + L.E*X*L.A')[:])
      mul!(temp, L.E, X)
      mul!(Y, temp, L.A')
      mul!(temp, X, L.E')
      mul!(Y, L.A, temp, 1, 1)
   end
   return y
end
function mulcsym!(y::AbstractVector, A::AbstractMatrix, E::AbstractMatrix, X::AbstractMatrix)
   require_one_based_indexing(y, A)
   # AXE' + EXA'
   n = size(A, 1)
   Y = similar(X, n, n)
   ZERO = zero(eltype(y))
   # Y = XE'
   mul!(Y, X, E')
   # AY + Y'A'
   @inbounds  begin
      k = 1
      for j = 1:n
         for i = 1:j
            temp = ZERO
            for l = 1:n
               temp += A[i,l] * Y[l,j] + conj(Y[l,i] * A[j,l])
            end
            y[k] = temp
            k += 1
         end
      end
   end
   return y
end
function muldsym!(y::AbstractVector, A::AbstractMatrix, E::AbstractMatrix, X::AbstractMatrix)
   require_one_based_indexing(y)
   # AXA' - EXE'
   n = size(A, 1)
   # t = triu(X)-diag(X)/2
   t = UpperTriangular(X) - Diagonal(X[diagind(X)] ./ 2)
   Y = similar(X, n, n)
   # Y = A*t*A' - E*t*E'
   mul!(Y, A * t, A')
   mul!(Y, E * t, E', -1, 1)
   # Y + Y'
   @inbounds  begin
      k = 1
      for j = 1:n
         for i = 1:j
            y[k] = Y[i,j] + conj(Y[j,i])
            k += 1
         end
      end
   end
   return y
end

struct InverseLyapunovMap{T,TA <: AbstractMatrix,adj,CD} <: LyapunovMatrixEquationsMaps{T}
   A::TA
   her::Bool
   sf::Bool
   function InverseLyapunovMap{T,TA,adj,CD}(A::TA, her::Bool, sf::Bool) where {T <: BlasFloat, TA <: AbstractMatrix{T}, adj, CD}
      LinearAlgebra.checksquare(A)
      return new{T,TA,adj,CD}(A, her, sf)
   end
end
function InverseLyapunovMap(A::AbstractMatrix{T}, ::CD = Continuous(), her::Bool = false) where {T <: BlasFloat,CD <: ContinuousOrDiscrete}
   schur_flag = isschur(A)
   return InverseLyapunovMap{T,typeof(A),false,CD}(A, her, schur_flag)
end
function InverseLyapunovMap(A::Adjoint{T,<:AbstractMatrix{T}}, ::CD = Continuous(), her::Bool = false) where {T <: BlasFloat,CD <: ContinuousOrDiscrete}
   Ap = A.parent
   schur_flag = isschur(Ap)
   return InverseLyapunovMap{T,typeof(Ap),true,CD}(Ap, her, schur_flag)
end

LinearAlgebra.adjoint(L::InverseLyapunovMap{<:Any,<:Any,adj,CD}) where {adj,CD} =
   InverseLyapunovMap(adj ? L.A : L.A', CD(), L.her)
LinearAlgebra.transpose(L::InverseLyapunovMap{<:Any,<:Any,adj,CD}) where {adj,CD} =
   InverseLyapunovMap(adj ? L.A : L.A', CD(), L.her)
LinearAlgebra.inv(L::LyapunovMap{T,TA,CD}) where {T,TA,CD} =
   InverseLyapunovMap(L.A, CD(), L.her)
LinearAlgebra.inv(L::InverseLyapunovMap{<:Any,<:Any,adj,CD}) where {adj,CD} =
   LyapunovMap(adj ? L.A' : L.A, CD(), L.her)

"""
    LINV = invlyapop(A; disc = false, her = false)

Define `LINV`, the inverse of the continuous Lyapunov operator `L:X -> AX+XA'` for `disc = false`
or the inverse of the discrete Lyapunov operator `L:X -> AXA'-X` for `disc = true`, where
`A` is an `n x n` matrix.
If `her = false` the inverse Lyapunov operator `LINV:Y -> X` maps general square matrices `Y`
into general square matrices `X`, and the associated matrix `M = Matrix(LINV)` is
``n^2 \\times n^2``.
If `her = true` the inverse Lyapunov operator `LINV:Y -> X` maps symmetric/Hermitian matrices `Y`
into symmetric/Hermitian matrices `X`, and the associated matrix `M = Matrix(LINV)` is
``n(n+1)/2 \\times n(n+1)/2``.
For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov
operators. Linear Algebra and its Applications 312:35–71, 2000.
"""
invlyapop(A::AbstractMatrix{<:BlasFloat}; disc=false, her=false) =
   InverseLyapunovMap(A, ifelse(disc, Discrete(), Continuous()), her)
function invlyapop(A::AbstractMatrix; disc=false, her=false)
   A = convert(AbstractMatrix{promote_type(eltype(A), Float64)}, A)
   InverseLyapunovMap(A, ifelse(disc, Discrete(), Continuous()), her)
end
invlyapop(A::Schur{<:BlasFloat}; disc=false, her=false) =
   InverseLyapunovMap(A.T, ifelse(disc, Discrete(), Continuous()), her)
function Base.size(L::InverseLyapunovMap)
   n = size(L.A, 1)
   N = L.her ? Int(n * (n + 1) / 2) : n * n
   return (N, N)
end
function mul!(y::AbstractVector, L::InverseLyapunovMap{T,<:Any,adj,Discrete}, x::AbstractVector) where {T <: BlasFloat,adj}
   n = size(L.A, 1)
   T1 = promote_type(T, eltype(x))
   try
      if L.sf && !(T <: Real && T1 <: Complex)
         if L.her
            Y = vec2triu(-convert(AbstractVector{T1}, x), her=true)
            lyapds!(L.A, Y, adj=adj)
            copyto!(y, triu2vec(Y))
         else
            Y = reshape(-convert(AbstractVector{T1}, x), n, n)
            adj ? sylvds!(-L.A, L.A, Y, adjA=true) : sylvds!(-L.A, L.A, Y, adjB=true)
            copyto!(y, Y)
         end
      else
         if L.her
            Y = vec2triu(-convert(AbstractVector{T1}, x), her=true)
            y .= triu2vec(lyapd(adj ? L.A' : L.A, Y))
         else
            Y = reshape(-convert(AbstractVector{T1}, x), n, n)
            copyto!(y, lyapd(adj ? L.A' : L.A, Y))
         end
      end
      return y
   catch err
      findfirst("SingularException", string(err)) === nothing &&
      findfirst("LAPACKException", string(err)) === nothing ? rethrow() : throw("ME:SingularException: Singular operator")
   end
end
function mul!(y::AbstractVector, L::InverseLyapunovMap{T,<:Any,adj,Continuous}, x::AbstractVector) where {T <: BlasFloat,adj}
   n = size(L.A, 1)
   T1 = promote_type(T, eltype(x))
   try
      if L.sf && !(T <: Real && T1 <: Complex)
         if L.her
            Y = vec2triu(-convert(AbstractVector{T1}, x), her=true)
            lyapcs!(L.A, Y; adj=adj)
            copyto!(y, triu2vec(Y))
         else
            Y = copy(reshape(convert(AbstractVector{T1}, x), n, n))
            adj ? (sylvcs!(L.A, L.A, Y, adjA=true)) : (sylvcs!(L.A, L.A, Y, adjB=true))
            copyto!(y, Y)
         end
      else
         if L.her
            Y = vec2triu(-convert(AbstractVector{T1}, x), her=true)
            y .= triu2vec(lyapc(adj ? L.A' : L.A, Y))
         else
            Y = reshape(-convert(AbstractVector{T1}, x), n, n)
            copyto!(y, lyapc(adj ? L.A' : L.A, Y))
         end
      end
      return y
   catch err
      findfirst("SingularException", string(err)) === nothing &&
      findfirst("LAPACKException", string(err)) === nothing ? rethrow() : throw("ME:SingularException: Singular operator")
   end
end

struct InverseGeneralizedLyapunovMap{T,TA <: AbstractMatrix,TE <: AbstractMatrix,adj,CD <: ContinuousOrDiscrete} <: LyapunovMatrixEquationsMaps{T}
  A::TA
  E::TE
  her::Bool
  sf::Bool
  function InverseGeneralizedLyapunovMap{T,TA,TE,adj,CD}(A::TA, E::TE, her::Bool, schur_flag::Bool) where {T<:BlasFloat,TA<:AbstractMatrix{T},TE<:AbstractMatrix{T},adj,CD}
      LinearAlgebra.checksquare(A)
      return new{T,TA,TE,adj,CD}(A, E, her, schur_flag)
   end
end

function InverseGeneralizedLyapunovMap(A::AbstractMatrix, E::AbstractMatrix, ::CD = Continuous(), her::Bool = false) where {CD <: ContinuousOrDiscrete}
   T = promote_type(eltype(A), eltype(E), Float64)
   A = convert(AbstractMatrix{T}, A)
   E = convert(AbstractMatrix{T}, E)
   schur_flag = isschur(A, E)
   return InverseGeneralizedLyapunovMap{T,typeof(A),typeof(E),false,CD}(A, E, her, schur_flag)
end
function InverseGeneralizedLyapunovMap(A::Adjoint{<:Any,<:AbstractMatrix}, E::Adjoint{<:Any,<:AbstractMatrix}, ::CD = Continuous(), her::Bool = false) where {CD <: ContinuousOrDiscrete}
   A = A.parent
   E = E.parent
   T = promote_type(eltype(A), eltype(E), Float64)
   A = convert(AbstractMatrix{T}, A)
   E = convert(AbstractMatrix{T}, E)
   schur_flag = isschur(A, E)
   return InverseGeneralizedLyapunovMap{T,typeof(A),typeof(E),true,CD}(A, E, her, schur_flag)
end

LinearAlgebra.adjoint(L::InverseGeneralizedLyapunovMap{<:Any,<:Any,<:Any,adj,CD}) where {adj,CD} =
   InverseGeneralizedLyapunovMap(adj ? L.A : L.A', adj ? L.E : L.E', CD(), L.her)
LinearAlgebra.transpose(L::InverseGeneralizedLyapunovMap{<:Any,<:Any,<:Any,adj,CD}) where {adj,CD} =
   InverseGeneralizedLyapunovMap(adj ? L.A : L.A', adj ? L.E : L.E', CD(), L.her)
LinearAlgebra.inv(L::GeneralizedLyapunovMap{<:Any,<:Any,<:Any,CD}) where {CD} =
   InverseGeneralizedLyapunovMap(L.A, L.E, CD(), L.her)
LinearAlgebra.inv(L::InverseGeneralizedLyapunovMap{<:Any,<:Any,<:Any,adj,CD}) where {adj,CD} =
   GeneralizedLyapunovMap(adj ? L.A' : L.A, adj ? L.E' : L.E, CD(), L.her)

"""
    LINV = invlyapop(A, E; disc = false, her = false)

Define `LINV`, the inverse of the continuous Lyapunov operator `L:X -> AXE'+EXA'` for `disc = false`
or the inverse of the discrete Lyapunov operator `L:X -> AXA'-EXE'` for `disc = true`, where
`(A,E)` is a pair of `n x n` matrices.
If `her = false` the inverse Lyapunov operator `LINV:Y -> X` maps general square matrices `Y`
into general square matrices `X`, and the associated matrix `M = Matrix(LINV)` is
``n^2 \\times n^2``.
If `her = true` the inverse Lyapunov operator `LINV:Y -> X` maps symmetric/Hermitian matrices `Y`
into symmetric/Hermitian matrices `X`, and the associated matrix `M = Matrix(LINV)` is
``n(n+1)/2 \\times n(n+1)/2``.
For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov
operators. Linear Algebra and its Applications 312:35–71, 2000.
"""
invlyapop(A::AbstractMatrix, E::AbstractMatrix; disc=false, her=false) =
   InverseGeneralizedLyapunovMap(A, E, ifelse(disc, Discrete(), Continuous()), her)
invlyapop(F::GeneralizedSchur; disc=false, her=false) =
   InverseGeneralizedLyapunovMap(F.S, F.T, ifelse(disc, Discrete(), Continuous()), her)
function Base.size(L::InverseGeneralizedLyapunovMap)
   n = size(L.A, 1)
   N = L.her ? Int(n * (n + 1) / 2) : n * n
   return (N, N)
end
function mul!(y::AbstractVector, L::InverseGeneralizedLyapunovMap{T,<:Any,<:Any,adj,Discrete}, x::AbstractVector) where {T,adj}
   n = size(L.A, 1)
   T1 = promote_type(T, eltype(x))
   try
      if L.sf && !(T <: Real && T1 <: Complex)
         if L.her
            Y = vec2triu(-convert(AbstractVector{T1}, x), her=true)
            lyapds!(L.A, L.E, Y; adj=adj)
            y .= triu2vec(Y)
         else
            copyto!(y, x)
            Y = reshape(y, (n, n))
            adj ? gsylvs!(L.A, L.A, -L.E, L.E, Y, adjAC=true) : gsylvs!(L.A, L.A, -L.E, L.E, Y, adjBD=true)
         end
      else
         if L.her
            Y = vec2triu(-convert(AbstractVector{T1}, x), her=true)
            y .= triu2vec(lyapd(adj ? L.A' : L.A, adj ? L.E' : L.E, Y))
         else
            Y = reshape(-convert(AbstractVector{T1}, x), n, n)
            # adj ? (y .= gsylv(-L.A',L.A,L.E',L.E,Y)[:]) : (y .= gsylv(-L.A,L.A',L.E,L.E',Y)[:])
            copyto!(y, lyapd(adj ? L.A' : L.A, adj ? L.E' : L.E, Y))
         end
      end
      return y
   catch err
      findfirst("SingularException", string(err)) === nothing &&
      findfirst("LAPACKException", string(err)) === nothing ? rethrow() : throw("ME:SingularException: Singular operator")
   end
end
function mul!(y::AbstractVector, L::InverseGeneralizedLyapunovMap{T,<:Any,<:Any,adj,Continuous}, x::AbstractVector) where {T,adj}
   n = size(L.A, 1)
   T1 = promote_type(T, eltype(x))
   try
      if L.sf && !(T <: Real && T1 <: Complex)
         if L.her
            Y = vec2triu(-convert(AbstractVector{T1}, x), her=true)
            lyapcs!(L.A, L.E, Y; adj=adj)
            y .= triu2vec(Y)
         else
            copyto!(y, x)
            Y = reshape(y, (n, n))
            adj ? (gsylvs!(L.A, L.E, L.E, L.A, Y, adjAC=true, DBSchur=true)) : (gsylvs!(L.A, L.E, L.E, L.A, Y, adjBD=true, DBSchur=true))
         end
      else
         if L.her
            Y = vec2triu(-convert(AbstractVector{T1}, x), her=true)
            y .= triu2vec(lyapc(adj ? L.A' : L.A, adj ? L.E' : L.E, Y))
         else
            Y = reshape(-convert(AbstractVector{T1}, x), n, n)
            # adj ? (y .= gsylv(L.A',L.E,L.E',L.A,Y)[:]) : (y .= gsylv(L.A,L.E',L.E,L.A',Y)[:])
            copyto!(y, lyapc(adj ? L.A' : L.A, adj ? L.E' : L.E, Y))
         end
      end
      return y
   catch err
      findfirst("SingularException", string(err)) === nothing &&
      findfirst("LAPACKException", string(err)) === nothing ? rethrow() : throw("ME:SingularException: Singular operator")
   end
end

struct SylvesterMap{T,TA <: AbstractMatrix,TB <: AbstractMatrix,CD <: ContinuousOrDiscrete} <: SylvesterMatrixEquationsMaps{T}
   A::TA
   B::TB
   function SylvesterMap{T,TA,TB,CD}(A::TA, B::TB) where {T,TA<:AbstractMatrix{T},TB<:AbstractMatrix{T},CD}
      LinearAlgebra.checksquare(A, B)
      return new{T,TA,TB,CD}(A, B)
   end
end
function SylvesterMap(A::AbstractMatrix, B::AbstractMatrix, ::CD = Continuous()) where {CD}
   T = promote_type(eltype(A), eltype(B))
   A = convert(AbstractMatrix{T}, A)
   B = convert(AbstractMatrix{T}, B)
   return SylvesterMap{T,typeof(A),typeof(B),CD}(A, B)
end

LinearAlgebra.adjoint(L::SylvesterMap{<:Any,<:Any,<:Any,CD}) where {CD} =
   SylvesterMap(L.A', L.B', CD())
LinearAlgebra.transpose(L::SylvesterMap{<:Any,<:Any,<:Any,CD}) where {CD} =
   SylvesterMap(L.A', L.B', CD())

"""
    M = sylvop(A, B; disc = false)

Define the continuous Sylvester operator `M: X -> AX+XB` if `disc = false`
or the discrete Sylvester operator `M: X -> AXB+X` if `disc = true`, where `A` and `B` are square matrices.
"""
sylvop(A::AbstractMatrix, B::AbstractMatrix; disc=false) =
   SylvesterMap(A, B, ifelse(disc, Discrete(), Continuous()))
sylvop(A::Schur, B::Schur; disc=false) = SylvesterMap(A.T, B.T, ifelse(disc, Discrete(), Continuous()))
function sylvop(A::Number, B::Number; disc=false)
   A, B = promote(A, B)
   SylvesterMap(fill(A, 1, 1), fill(B, 1, 1), ifelse(disc, Discrete(), Continuous()))
end
Base.size(L::SylvesterMap) = (N = size(L.A, 1) * size(L.B, 1); return (N, N))
function mul!(y::AbstractVector, L::SylvesterMap{T,<:Any,<:Any,Discrete}, x::AbstractVector) where T
   m = size(L.A, 1)
   n = size(L.B, 1)
   T1 = promote_type(T, eltype(x))
   X = reshape(convert(AbstractVector{T1}, x), (m, n))
   Y = reshape(y, (m, n))
   # Y = A * X * B + X
   copyto!(y, x)
   mul!(Y, L.A * X, L.B, true, true)
   return y
end
function mul!(y::AbstractVector, L::SylvesterMap{T,<:Any,<:Any,Continuous}, x::AbstractVector) where T
   m = size(L.A, 1)
   n = size(L.B, 1)
   T1 = promote_type(T, eltype(x))
   X = reshape(convert(AbstractVector{T1}, x), (m, n))
   Y = reshape(y, (m, n))
   # Y = A * X + X * B
   mul!(Y, L.A, X)
   mul!(Y, X, L.B, true, true)
   return y
end

struct InverseSylvesterMap{T,TA <: AbstractMatrix,TB <: AbstractMatrix,CD} <: SylvesterMatrixEquationsMaps{T}
   A::TA
   B::TB
   adjA::Bool
   adjB::Bool
   sf::Bool
   function InverseSylvesterMap{T,TA,TB,CD}(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where {T,TA,TB,CD}
      LinearAlgebra.checksquare(A, B)
      adjA = isa(A, Adjoint)
      adjB = isa(B, Adjoint)
      A = adjA ? A.parent : A
      B = adjB ? B.parent : B
      schur_flag = isschur(A) && isschur(B)
      return new{T,typeof(A),typeof(B),CD}(A, B, adjA, adjB, schur_flag)
   end
end
function InverseSylvesterMap(A::AbstractMatrix, B::AbstractMatrix, ::CD = Continuous) where {CD}
   T = promote_type(eltype(A), eltype(B))
   A = convert(AbstractMatrix{T}, A)
   B = convert(AbstractMatrix{T}, B)
   return InverseSylvesterMap{T,typeof(A),typeof(B),CD}(A, B)
end

LinearAlgebra.adjoint(L::InverseSylvesterMap{<:Any,<:Any,<:Any,CD}) where {CD} =
   InverseSylvesterMap(L.adjA ? L.A : L.A', L.adjB ? L.B : L.B', CD())
LinearAlgebra.transpose(L::InverseSylvesterMap{<:Any,<:Any,<:Any,CD}) where {CD} =
   InverseSylvesterMap(L.adjA ? L.A : L.A', L.adjB ? L.B : L.B', CD())
LinearAlgebra.inv(L::SylvesterMap{<:Any,<:Any,<:Any,CD}) where {CD} =
   InverseSylvesterMap(L.A, L.B, CD())
LinearAlgebra.inv(L::InverseSylvesterMap{<:Any,<:Any,<:Any,CD}) where {CD} =
   SylvesterMap(L.adjA ? L.A' : L.A, L.adjB ? L.B' : L.B, CD())

"""
    MINV = invsylvop(A, B; disc = false)

Define `MINV`, the inverse of the continuous Sylvester operator  `M: X -> AX+XB` if `disc = false`
or of the discrete Sylvester operator `M: X -> AXB+X` if `disc = true`, where `A` and `B` are square matrices.
"""
function invsylvop(A::AbstractMatrix, B::AbstractMatrix; disc=false)
   InverseSylvesterMap(A, B, ifelse(disc, Discrete(), Continuous()))
end
invsylvop(A::Schur, B::Schur; disc=false) =
   InverseSylvesterMap(A.T, B.T, ifelse(disc, Discrete(), Continuous()))
function invsylvop(A::Number, B::Number; disc=false)
   A, B = promote(A, B)
   InverseSylvesterMap(fill(A, 1, 1), fill(B, 1, 1), ifelse(disc, Discrete(), Continuous()))
end
Base.size(L::InverseSylvesterMap) = (N = size(L.A, 1) * size(L.B, 1); return (N, N))
function mul!(y::AbstractVector, L::InverseSylvesterMap{T,<:Any,<:Any,Discrete}, x::AbstractVector) where {T}
   m = size(L.A, 1)
   n = size(L.B, 1)
   T1 = promote_type(T, eltype(x))
   try
      if L.sf
         copyto!(y, x)
         Y = reshape(y, (m, n))
         sylvds!(L.A, L.B, Y, adjA=L.adjA, adjB=L.adjB)
      else
         Y = reshape(convert(AbstractVector{T1}, x), m, n)
         copyto!(y, sylvd(L.adjA ? L.A' : L.A, L.adjB ? L.B' : L.B, Y))
      end
      return y
   catch err
      findfirst("SingularException", string(err)) === nothing &&
      findfirst("LAPACKException", string(err)) === nothing ? rethrow() : throw("ME:SingularException: Singular operator")
   end
end
function mul!(y::AbstractVector, L::InverseSylvesterMap{T,<:Any,<:Any,Continuous}, x::AbstractVector) where T
   m = size(L.A, 1)
   n = size(L.B, 1)
   T1 = promote_type(T, eltype(x))
   try
      if L.sf
         copyto!(y, x)
         Y = reshape(y, (m, n))
         sylvcs!(L.A, L.B, Y, adjA=L.adjA, adjB=L.adjB)
      else
         Y = reshape(convert(AbstractVector{T1}, x), m, n)
         copyto!(y, sylvc(L.adjA ? L.A' : L.A, L.adjB ? L.B' : L.B, Y))
      end
      return y
   catch err
      findfirst("SingularException", string(err)) === nothing &&
      findfirst("LAPACKException", string(err)) === nothing ? rethrow() : throw("ME:SingularException: Singular operator")
   end
end

struct GeneralizedSylvesterMap{T,TA <: AbstractMatrix,TB <: AbstractMatrix,TC <: AbstractMatrix,TD <: AbstractMatrix} <: SylvesterMatrixEquationsMaps{T}
   A::TA
   B::TB
   C::TC
   D::TD
   function GeneralizedSylvesterMap(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix)
      m, n = LinearAlgebra.checksquare(A, B)
      [m, n] == LinearAlgebra.checksquare(C, D) ||
               throw(DimensionMismatch("A, B, C and D have incompatible dimensions"))
      T = promote_type(map(eltype, (A, B, C, D))...)
      A = convert(AbstractMatrix{T}, A)
      B = convert(AbstractMatrix{T}, B)
      C = convert(AbstractMatrix{T}, C)
      D = convert(AbstractMatrix{T}, D)
      return new{T,typeof(A),typeof(B),typeof(C),typeof(D)}(A, B, C, D)
   end
end

LinearAlgebra.adjoint(L::GeneralizedSylvesterMap{T}) where T =
   GeneralizedSylvesterMap(L.A', L.B', L.C', L.D')
LinearAlgebra.transpose(L::GeneralizedSylvesterMap{T}) where T =
   GeneralizedSylvesterMap(L.A', L.B', L.C', L.D')

"""
    M = sylvop(A, B, C, D)

Define the generalized Sylvester operator `M: X -> AXB+CXD`, where `(A,C)` and `(B,D)` a pairs of square matrices.
"""
sylvop(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix) =
   GeneralizedSylvesterMap(A, B, C, D)
Base.size(L::GeneralizedSylvesterMap) = (N = size(L.A, 1) * size(L.B, 1); return (N, N))
function mul!(y::AbstractVector, L::GeneralizedSylvesterMap{T}, x::AbstractVector) where T
   m = size(L.A, 1)
   n = size(L.B, 1)
   T1 = promote_type(T, eltype(x))
   X = reshape(convert(AbstractVector{T1}, x), (m, n))
   Y = reshape(y, (m, n))
   # Y = A * X * B + C * X * D
   temp = similar(Y, (m, n))
   mul!(temp, L.A, X)
   mul!(Y, temp, L.B, 1, 0)
   mul!(temp, L.C, X)
   mul!(Y, temp, L.D, 1, 1)
   return y
end

struct InverseGeneralizedSylvesterMap{T,TA <: AbstractMatrix,TB <: AbstractMatrix,TC <: AbstractMatrix,TD <: AbstractMatrix} <: SylvesterMatrixEquationsMaps{T}
   A::TA
   B::TB
   C::TC
   D::TD
   adjAC::Bool
   adjBD::Bool
   sf::Bool
   BDSchur::Bool
   DBSchur::Bool
   function InverseGeneralizedSylvesterMap(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix)
      m, n = LinearAlgebra.checksquare(A, B)
      [m, n] == LinearAlgebra.checksquare(C, D) ||
               throw(DimensionMismatch("A, B, C and D have incompatible dimensions"))
      adjAC = isa(A, Adjoint) && isa(C, Adjoint)
      adjBD = isa(B, Adjoint) && isa(D, Adjoint)
      ACSchur = (adjAC ? isschur(A.parent, C.parent) : isschur(A, C))
      BDSchur = (adjBD ? isschur(B.parent, D.parent) : isschur(B, D))
      DBSchur = (adjBD ? isschur(D.parent, B.parent) : isschur(D, B))
      schur_flag = ACSchur && (BDSchur || DBSchur)
      A = adjAC ? A.parent : A
      B = adjBD ? B.parent : B
      C = adjAC ? C.parent : C
      D = adjBD ? D.parent : D
      T = promote_type(map(eltype, (A, B, C, D))...)
      A = convert(AbstractMatrix{T}, A)
      B = convert(AbstractMatrix{T}, B)
      C = convert(AbstractMatrix{T}, C)
      D = convert(AbstractMatrix{T}, D)
    return new{T,typeof(A),typeof(B),typeof(C),typeof(D)}(A, B, C, D, adjAC, adjBD, schur_flag, BDSchur, DBSchur)
   end
end

LinearAlgebra.adjoint(L::InverseGeneralizedSylvesterMap) =
   InverseGeneralizedSylvesterMap(L.adjAC ? L.A : L.A', L.adjBD ? L.B : L.B',
                                    L.adjAC ? L.C : L.C', L.adjBD ? L.D : L.D')
LinearAlgebra.transpose(L::InverseGeneralizedSylvesterMap) =
   InverseGeneralizedSylvesterMap(L.adjAC ? L.A : L.A', L.adjBD ? L.B : L.B',
                                    L.adjAC ? L.C : L.C', L.adjBD ? L.D : L.D')
LinearAlgebra.inv(L::GeneralizedSylvesterMap) =
   InverseGeneralizedSylvesterMap(L.A, L.B, L.C, L.D)
LinearAlgebra.inv(L::InverseGeneralizedSylvesterMap) =
   GeneralizedSylvesterMap(L.adjAC ? L.A' : L.A, L.adjBD ? L.B' : L.B,
                           L.adjAC ? L.C' : L.C, L.adjBD ? L.D' : L.D)

"""
    MINV = invsylvop(A, B, C, D)

Define `MINV`, the inverse of the generalized Sylvester operator `M: X -> AXB+CXD`,
where (A,C) and (B,D) a pairs of square matrices.
"""
invsylvop(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix) =
   InverseGeneralizedSylvesterMap(A, B, C, D)
invsylvop(F::GeneralizedSchur, G::GeneralizedSchur) = invsylvop(F.S, G.S, F.T, G.T)
Base.size(L::InverseGeneralizedSylvesterMap) = (N = size(L.A, 1) * size(L.B, 1); return (N, N))
function mul!(y::AbstractVector, L::InverseGeneralizedSylvesterMap{T}, x::AbstractVector) where T
   m = size(L.A, 1)
   n = size(L.B, 1)
   T1 = promote_type(T, eltype(x))
   try
      if L.sf
         copyto!(y, x)
         Y = reshape(y, m, n)
         gsylvs!(L.A, L.B, L.C, L.D, Y, adjAC=L.adjAC, adjBD=L.adjBD, DBSchur=L.DBSchur && !L.BDSchur)
      else
         X = reshape(convert(AbstractVector{T1}, x), m, n)
         copyto!(y, gsylv(L.adjAC ? L.A' : L.A, L.adjBD ? L.B' : L.B,
                           L.adjAC ? L.C' : L.C, L.adjBD ? L.D' : L.D, X))
      end
   catch err
      findfirst("SingularException", string(err)) === nothing &&
      findfirst("LAPACKException", string(err)) === nothing ? rethrow() : throw("ME:SingularException: Singular operator")
   end
   return y
end

struct SylvesterSystemMap{T,TA <: AbstractMatrix,TB <: AbstractMatrix,TC <: AbstractMatrix,TD <: AbstractMatrix} <: SylvesterMatrixEquationsMaps{T}
   A::TA
   B::TB
   C::TC
   D::TD
   function SylvesterSystemMap(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix)
      m, n = LinearAlgebra.checksquare(A, B)
      [m, n] == LinearAlgebra.checksquare(C, D) ||
         throw(DimensionMismatch("A, B, C and D have incompatible dimensions"))
      T = promote_type(map(eltype, (A, B, C, D))...)
      A = convert(AbstractMatrix{T}, A)
      B = convert(AbstractMatrix{T}, B)
      C = convert(AbstractMatrix{T}, C)
      D = convert(AbstractMatrix{T}, D)
      return new{T,typeof(A),typeof(B),typeof(C),typeof(D)}(A, B, C, D)
   end
end

"""
    M = sylvsysop(A, B, C, D)

Define the operator `M: (X,Y) -> (AX+YB, CX+YD)`, where `(A,C)` and `(B,D)` are pairs of
square matrices.
"""
sylvsysop(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix) =
   SylvesterSystemMap(A, B, C, D)
Base.size(L::SylvesterSystemMap) = (N = size(L.A, 1) * size(L.B, 1) * 2; return (N, N))
function mul!(y::AbstractVector, L::SylvesterSystemMap{T}, x::AbstractVector) where T
   require_one_based_indexing(y, x)
   m = size(L.A, 1)
   n = size(L.B, 1)
   mn = m * n
   T1 = promote_type(T, eltype(x))
   X1 = reshape(convert(AbstractVector{T1}, view(x, 1:mn)), m, n)
   X2 = reshape(convert(AbstractVector{T1}, view(x, mn + 1:2 * mn)), m, n)
   U = reshape(view(y, 1:m * n), (m, n))
   V = reshape(view(y, mn + 1:2 * mn), (m, n))
   # U = A * X + Y * B
   # V = C * X + Y * D
   mul!(U, L.A, X1)
   mul!(U, X2, L.B, 1, 1)
   mul!(V, L.C, X1)
   mul!(V, X2, L.D, 1, 1)
   return y
end
for ttype in (LinearMaps.TransposeMap, LinearMaps.AdjointMap)
   @eval function mul!(y::AbstractVector, LT::$ttype{T,<:SylvesterSystemMap{T}}, x::AbstractVector) where T
      require_one_based_indexing(y, x)
      m = size(LT.lmap.A, 1)
      n = size(LT.lmap.B, 1)
      mn = m * n
      T1 = promote_type(T, eltype(x))
      X1 = reshape(convert(AbstractVector{T1}, view(x, 1:mn)), m, n)
      X2 = reshape(convert(AbstractVector{T1}, view(x, mn + 1:2 * mn)), m, n)
      U = reshape(view(y, 1:m * n), (m, n))
      V = reshape(view(y, mn + 1:2 * mn), (m, n))
      # U = A' * X + C' * Y
      # V = X * B' + Y * D'
      mul!(U, LT.lmap.A', X1)
      mul!(U, LT.lmap.C', X2, 1, 1)
      mul!(V, X1, LT.lmap.B')
      mul!(V, X2, LT.lmap.D', 1, 1)
      return y
   end
end

struct InverseSylvesterSystemMap{T,TA <: AbstractMatrix,TB <: AbstractMatrix,TC <: AbstractMatrix,TD <: AbstractMatrix} <: SylvesterMatrixEquationsMaps{T}
   A::TA
   B::TB
   C::TC
   D::TD
   sf::Bool
   function InverseSylvesterSystemMap(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix)
      m, n = LinearAlgebra.checksquare(A, B)
      [m, n] == LinearAlgebra.checksquare(C, D) ||
               throw(DimensionMismatch("A, B, C and D have incompatible dimensions"))
      T = promote_type(map(eltype, (A, B, C, D))...)
      A = convert(AbstractMatrix{T}, A)
      B = convert(AbstractMatrix{T}, B)
      C = convert(AbstractMatrix{T}, C)
      D = convert(AbstractMatrix{T}, D)
      schur_flag = isschur(A, C) && isschur(B, D)
      return new{T,typeof(A),typeof(B),typeof(C),typeof(D)}(A, B, C, D, schur_flag)
   end
end

LinearAlgebra.inv(L::SylvesterSystemMap) = InverseSylvesterSystemMap(L.A, L.B, L.C, L.D)
LinearAlgebra.inv(L::InverseSylvesterSystemMap) = SylvesterSystemMap(L.A, L.B, L.C, L.D)

"""
    MINV = invsylvsysop(A, B, C, D)

Define `MINV`, the inverse of the linear operator `M: (X,Y) -> (AX+YB, CX+YD )`,
where `(A,C)` and `(B,D)` a pairs of square matrices.
"""
invsylvsysop(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix) =
   InverseSylvesterSystemMap(A, B, C, D)
invsylvsysop(AC::GeneralizedSchur, BD::GeneralizedSchur) =
   InverseSylvesterSystemMap(AC.S, BD.S, AC.T, BD.T)
Base.size(L::InverseSylvesterSystemMap) = (N = size(L.A, 1) * size(L.B, 1) * 2; return (N, N))
function mul!(y::AbstractVector, L::InverseSylvesterSystemMap{T}, x::AbstractVector) where T
   require_one_based_indexing(y, x)
   m = size(L.A, 1)
   n = size(L.B, 1)
   T1 = promote_type(T, eltype(x))
   mn = m * n
   E = reshape(eltype(x) == T1 ? x[1:mn] : convert(Vector{T1}, view(x, 1:m * n)), m, n)
   F = reshape(eltype(x) == T1 ? x[mn + 1:2 * mn] : convert(Vector{T1}, view(x, mn + 1:2 * mn)), m, n)
   Y1 = view(y, 1:mn)
   Y2 = view(y, mn + 1:2mn)
   try
      if L.sf
         X1, X2 = sylvsyss!(L.A, L.B, E, L.C, L.D, F)
         copyto!(Y1, X1)
         copyto!(Y2, X2)
      else
         X1, X2  = sylvsys(L.A, L.B, E, L.C, L.D, F)
         copyto!(Y1, X1)
         copyto!(Y2, X2)
      end
   catch err
      findfirst("LAPACKException", string(err)) === nothing ? rethrow() : throw("ME:SingularException: Singular operator")
   end
   return y
end
for ttype in (LinearMaps.TransposeMap, LinearMaps.AdjointMap)
   @eval function mul!(y::AbstractVector, L::$ttype{T,<:InverseSylvesterSystemMap}, x::AbstractVector) where T
      require_one_based_indexing(y, x)
      m = size(L.lmap.A, 1)
      n = size(L.lmap.B, 1)
      mn = m * n
      T1 = promote_type(T, eltype(x))
      E = reshape(eltype(x) == T1 ? x[1:mn] : convert(Vector{T1}, view(x, 1:mn)), m, n)
      F = reshape(eltype(x) == T1 ? x[mn + 1:2 * mn] : convert(Vector{T1}, view(x, mn + 1:2 * mn)), m, n)
      Y1 = view(y, 1:mn)
      Y2 = view(y, mn + 1:2mn)
      try
         if L.lmap.sf
            X1, X2 = dsylvsyss!(L.lmap.A, L.lmap.B, E, L.lmap.C, L.lmap.D, F)
            copyto!(Y1, X1)
            copyto!(Y2, X2)
         else
            X1, X2 = dsylvsys(L.lmap.A', L.lmap.B', E, L.lmap.C', L.lmap.D', F)
            copyto!(Y1, X1)
            copyto!(Y2, X2)
         end
      catch err
         findfirst("LAPACKException", string(err)) === nothing ? rethrow() : throw("ME:SingularException: Singular operator")
      end
      return y
   end
end

max_stars_repo_name: andreasvarga/MatrixEquations

max_issues_repo_name: andreasvarga/MatrixEquations

max_forks_repo_name: andreasvarga/MatrixEquations


text:
module AWS

using Compat: Compat, @something
using Base64
using Dates
using Downloads: Downloads, Downloader, Curl
using HTTP
using MbedTLS
using Mocking
using OrderedCollections: LittleDict, OrderedDict
using Retry
using Sockets
using URIs
using UUIDs: UUIDs
using XMLDict

export @service
export _merge
export AbstractAWSConfig, AWSConfig, AWSExceptions, AWSServices, Request
export ec2_instance_metadata, ec2_instance_region
export generate_service_url, global_aws_config, set_user_agent
export sign!, sign_aws2!, sign_aws4!
export JSONService, RestJSONService, RestXMLService, QueryService, set_features

const DEFAULT_REGION = "us-east-1"

include(joinpath("utilities", "utilities.jl"))

include("AWSExceptions.jl")
include("AWSCredentials.jl")
include("AWSConfig.jl")
include("AWSMetadata.jl")

include(joinpath("utilities", "request.jl"))
include(joinpath("utilities", "response.jl"))
include(joinpath("utilities", "sign.jl"))
include(joinpath("utilities", "downloads_backend.jl"))

include("deprecated.jl")

using ..AWSExceptions
using ..AWSExceptions: AWSException

const user_agent = Ref("AWS.jl/1.0.0")
const aws_config = Ref{AbstractAWSConfig}()

"""
    FeatureSet

Allows end users to opt-in to new breaking behaviors prior before a major/breaking package
release. Each field of this struct contains a default which specifies uses the original
non-breaking behavior.

# Features
- `use_response_type::Bool=false`: When enabled, service calls will return an `AWS.Response`
  which provides streaming/raw/parsed access to the response. When disabled, the service
  call response typically be parsed but will vary depending on the following parameters:
  "return_headers", "return_stream", "return_raw", "response_dict_type".
"""
Base.@kwdef struct FeatureSet
    use_response_type::Bool = false
end

"""
    global_aws_config()

Retrieve the global AWS configuration.
If one is not set, create one with default configuration options.

# Keywords
- `kwargs...`: AWSConfig kwargs to be passed along if the global configuration is not already set

# Returns
- `AWSConfig`: The global AWS configuration
"""
function global_aws_config(; kwargs...)
    if !isassigned(aws_config) || !isempty(kwargs)
        aws_config[] = AWSConfig(; kwargs...)
    end

    return aws_config[]
end

"""
    global_aws_config(config::AbstractAWSConfig)

Set the global AWSConfig.

# Arguments
- `config::AWSConfig`: The AWSConfig to set in the global state

# Returns
- `AWSConfig`: Global AWSConfig
"""
function global_aws_config(config::AbstractAWSConfig)
    return aws_config[] = config
end

"""
    set_user_agent(new_user_agent::String)

Set the global user agent when making HTTP requests.

# Arguments
- `new_user_agent::String`: User agent to set when making HTTP requests

# Return
- `String`: The global user agent
"""
set_user_agent(new_user_agent::String) = return user_agent[] = new_user_agent

"""
    @service module_name feature=val...

Include a high-level service wrapper based off of the `module_name` parameter optionally
supplying a list of `features`.

When calling the macro you cannot match the predefined constant for the low-level API.
The low-level API constants are named in all lowercase, and spaces replaced with underscores.

Examples:
```julia
using AWS.AWSServices: secrets_manager
using AWS: @service

# This matches the constant and will error!
@service secrets_manager
> ERROR: cannot assign a value to variable AWSServices.secrets_manager from module Main

# This does NOT match the filename structure and will error!
@service secretsmanager
> ERROR: could not open file /.julia/dev/AWS.jl/src/services/secretsmanager.jl

# All of the examples below are valid!
@service Secrets_Manager
@service SECRETS_MANAGER
@service sECRETS_MANAGER

# Using a feature
@service Secrets_Manager use_response_type = true
```

# Arguments
- `module_name::Symbol`: Name of the module and service to include high-level API wrappers in your namespace
- `features=val...`: A list of features to enable/disable for this high-level API include.
  See `FeatureSet` for a list of available features.

# Return
- `Expression`: Module which embeds the high-level service API wrapper functions in your namespace
"""
macro service(module_name::Symbol, features...)
    service_name = joinpath(@__DIR__, "services", lowercase(string(module_name)) * ".jl")
    map(_assignment_to_kw!, features)

    module_block = quote
        using AWS: FeatureSet
        const SERVICE_FEATURE_SET = FeatureSet(; $(features...))
        include($service_name)
    end

    return Expr(:toplevel, Expr(:module, true, esc(module_name), esc(module_block)))
end

abstract type Service end

struct RestXMLService <: Service
    signing_name::String
    endpoint_prefix::String
    api_version::String
end

struct QueryService <: Service
    signing_name::String
    endpoint_prefix::String
    api_version::String
end

struct JSONService <: Service
    signing_name::String
    endpoint_prefix::String
    api_version::String

    json_version::String
    target::String
end

struct RestJSONService <: Service
    signing_name::String
    endpoint_prefix::String
    api_version::String

    service_specific_headers::LittleDict{String,String}
end

function RestJSONService(signing_name::String, endpoint_prefix::String, api_version::String)
    return RestJSONService(
        signing_name, endpoint_prefix, api_version, LittleDict{String,String}()
    )
end

struct ServiceWrapper{S<:Service}
    service::S
    feature_set::FeatureSet
end

function set_features(service::Service; features...)
    return ServiceWrapper(service, FeatureSet(; features...))
end

# Needs to be included after the definition of struct otherwise it cannot find them
include("AWSServices.jl")

function generate_service_url(aws::AbstractAWSConfig, service::String, resource::String)
    SERVICE_HOST = "amazonaws.com"
    reg = region(aws)

    regionless_services = ("iam", "route53")

    if service in regionless_services || (service == "sdb" && reg == DEFAULT_REGION)
        reg = ""
    end

    return string(
        "https://", service, ".", isempty(reg) ? "" : "$reg.", SERVICE_HOST, resource
    )
end

"""
    (service::RestXMLService)(
        request_method::String, request_uri::String, args::AbstractDict{String, <:Any}=Dict{String, String}();
        aws::AbstractAWSConfig=aws_config
    )

Perform a RestXML request to AWS.

# Arguments
- `request_method::String`: RESTful request type, e.g.: `GET`, `HEAD`, `PUT`, etc.
- `request_uri::String`: AWS URI for the endpoint
- `args::AbstractDict{String, <:Any}`: Additional arguments to be included in the request

# Keywords
- `aws::AbstractAWSConfig`: AWSConfig containing credentials and other information for fulfilling the request, default value is the global configuration
- `feature_set::FeatureSet`: Specifies opt-in functionality for this specific API call.

# Returns
- `Tuple` or `Dict`: If `return_headers` is passed in through `args` a Tuple containing the headers and response will be returned, otherwise just a `Dict`
"""
function (service::RestXMLService)(
    request_method::String,
    request_uri::String,
    args::AbstractDict{String,<:Any}=Dict{String,Any}();
    aws_config::AbstractAWSConfig=global_aws_config(),
    feature_set::FeatureSet=FeatureSet(),
)
    feature_set.use_response_type && _delete_legacy_response_kw_args!(args)

    return_headers = _pop!(args, "return_headers", nothing)

    request = Request(;
        _extract_common_kw_args(service, args)...,
        use_response_type=feature_set.use_response_type,
        request_method=request_method,
        content=_pop!(args, "body", ""),
    )

    if request.service == "s3"
        request_uri = _clean_s3_uri(request_uri)
    end
    request.resource = _generate_rest_resource(request_uri, args)
    query_str = HTTP.escapeuri(args)

    if !isempty(query_str)
        if occursin('?', request.resource)
            request.resource *= "&$query_str"
        else
            request.resource *= "?$query_str"
        end
    end

    request.url = generate_service_url(
        aws_config, service.endpoint_prefix, request.resource
    )

    return submit_request(aws_config, request; return_headers=return_headers)
end

"""
    (service::QueryService)(
        operation::String, args::AbstractDict{String, <:Any}=Dict{String, Any}();
        aws::AbstractAWSConfig=aws_config
    )

Perform a Query request to AWS.

# Arguments
- `operation::String`:
- `args::AbstractDict{String, <:Any}`: Additional arguments to be included in the request

# Keywords
- `aws::AbstractAWSConfig`: AWSConfig containing credentials and other information for fulfilling the request, default value is the global configuration
- `feature_set::FeatureSet`: Specifies opt-in functionality for this specific API call.

# Returns
- `Tuple` or `Dict`: If `return_headers` is passed in through `args` a Tuple containing the headers and response will be returned, otherwise just a `Dict`
"""
function (service::QueryService)(
    operation::String,
    args::AbstractDict{String,<:Any}=Dict{String,Any}();
    aws_config::AbstractAWSConfig=global_aws_config(),
    feature_set::FeatureSet=FeatureSet(),
)
    feature_set.use_response_type && _delete_legacy_response_kw_args!(args)

    POST_RESOURCE = "/"
    return_headers = _pop!(args, "return_headers", nothing)

    request = Request(;
        _extract_common_kw_args(service, args)...,
        use_response_type=feature_set.use_response_type,
        resource=POST_RESOURCE,
        request_method="POST",
        url=generate_service_url(aws_config, service.endpoint_prefix, POST_RESOURCE),
    )

    request.headers["Content-Type"] = "application/x-www-form-urlencoded; charset=utf-8"

    args["Action"] = operation
    args["Version"] = service.api_version
    request.content = HTTP.escapeuri(_flatten_query(service.signing_name, args))

    return submit_request(aws_config, request; return_headers=return_headers)
end

"""
    (service::JSONService)(
        operation::String, args::AbstractDict{String, <:Any}=Dict{String, Any}();
        aws::AbstractAWSConfig=aws_config
    )

Perform a JSON request to AWS.

# Arguments
- `operation::String`: Name of the operation to perform
- `args::AbstractDict{String, <:Any}`: Additional arguments to be included in the request

# Keywords
- `aws::AbstractAWSConfig`: AWSConfig containing credentials and other information for fulfilling the request, default value is the global configuration
- `feature_set::FeatureSet`: Specifies opt-in functionality for this specific API call.

# Returns
- `Tuple` or `Dict`: If `return_headers` is passed in through `args` a Tuple containing the headers and response will be returned, otherwise just a `Dict`
"""
function (service::JSONService)(
    operation::String,
    args::AbstractDict{String,<:Any}=Dict{String,Any}();
    aws_config::AbstractAWSConfig=global_aws_config(),
    feature_set::FeatureSet=FeatureSet(),
)
    feature_set.use_response_type && _delete_legacy_response_kw_args!(args)

    POST_RESOURCE = "/"
    return_headers = _pop!(args, "return_headers", nothing)

    request = Request(;
        _extract_common_kw_args(service, args)...,
        use_response_type=feature_set.use_response_type,
        resource=POST_RESOURCE,
        request_method="POST",
        content=json(args),
        url=generate_service_url(aws_config, service.endpoint_prefix, POST_RESOURCE),
    )

    request.headers["Content-Type"] = "application/x-amz-json-$(service.json_version)"
    request.headers["X-Amz-Target"] = "$(service.target).$(operation)"

    return submit_request(aws_config, request; return_headers=return_headers)
end

"""
    (service::RestJSONService)(
        request_method::String, request_uri::String, args::AbstractDict{String, <:Any}=Dict{String, String}();
        aws::AbstractAWSConfig=aws_config
    )

Perform a RestJSON request to AWS.

# Arguments
- `request_method::String`: RESTful request type, e.g.: `GET`, `HEAD`, `PUT`, etc.
- `request_uri::String`: AWS URI for the endpoint
- `args::AbstractDict{String, <:Any}`: Additional arguments to be included in the request

# Keywords
- `aws::AbstractAWSConfig`: AWSConfig containing credentials and other information for fulfilling the request, default value is the global configuration
- `feature_set::FeatureSet`: Specifies opt-in functionality for this specific API call.

# Returns
- `Tuple` or `Dict`: If `return_headers` is passed in through `args` a Tuple containing the headers and response will be returned, otherwise just a `Dict`
"""
function (service::RestJSONService)(
    request_method::String,
    request_uri::String,
    args::AbstractDict{String,<:Any}=Dict{String,String}();
    aws_config::AbstractAWSConfig=global_aws_config(),
    feature_set::FeatureSet=FeatureSet(),
)
    feature_set.use_response_type && _delete_legacy_response_kw_args!(args)

    return_headers = _pop!(args, "return_headers", nothing)

    request = Request(;
        _extract_common_kw_args(service, args)...,
        use_response_type=feature_set.use_response_type,
        request_method=request_method,
        resource=_generate_rest_resource(request_uri, args),
    )

    request.url = generate_service_url(
        aws_config, service.endpoint_prefix, request.resource
    )

    if !isempty(service.service_specific_headers)
        merge!(request.headers, service.service_specific_headers)
    end

    request.headers["Content-Type"] = "application/json"
    request.content = json(args)

    return submit_request(aws_config, request; return_headers=return_headers)
end

function (service::ServiceWrapper)(args...; feature_set=nothing, kwargs...)
    feature_set = something(feature_set, service.feature_set)
    return service.service(args...; feature_set=feature_set, kwargs...)
end

function __init__()
    DEFAULT_BACKEND[] = HTTPBackend()
    return nothing
end

end  # module AWS


text:
"""
    function_to_interpolated_function(f_xy; physwidth = 10.0m, physheight = 4.0m)

Input: f: (Q, Q)  → (Q, Q) Function taking a tuple quantities, outputs a tuple of quantities
    physwidth:  Q
    physheight: Q

Output: f: (Q, Q)  → (Q, Q) Function taking a tuple quantities, outputs a tuple of quantities

Convert the function f_xy into an interpolated function. The output function is linearly interpolated
between pixels, based on the currently defined sketch scale.

The output of the function can be otherwise; this is just the originally intended use.
"""
function function_to_interpolated_function(f_xy; physwidth = 10.0m, physheight = 4.0m, centered = true)
    # Position iterators
    xs, ys = x_y_iterators_at_pixels(;physwidth, physheight, centered)
    tuplemat = [f_xy(x, y) for x in xs, y in reverse(ys)]
    fxy_inter = interpolate((xs, reverse(ys)),
        tuplemat,
        Gridded(Linear()));
    # Extend the domain, using the same values as on the domain border. Don't know how the closest point
    # on domain boundary is found.
    extrapolate(fxy_inter, Flat())
end

"""
    matrix_to_function(matrix::Matrix, xmin, xmax, ymin, ymax, nx, ny)
    -> f: (Q, Q)  → (Q, Q)

Make a linear interpolation function from (complex) matrix.
This could be useful for e.g. a flow field, where the matrix contains complex
quantity values.

The imaginary component is mapped to the second parameter.
"""
function matrix_to_function(matrix::Matrix, xmin, xmax, ymin, ymax, nx, ny)
    xs = range(xmin, xmax, length = nx)
    ys = range(ymin, ymax, length = ny)
    # Adapt to Interpolations
    tuplemat = map( cmplx -> (real(cmplx), imag(cmplx)), transpose(matrix)[ : , end:-1:1])
    # Function that can interpolate between nodes.
    fxy_inter = interpolate((xs, reverse(ys)),
        tuplemat,
        Gridded(Linear()));
    # Extend the domain, using the same values as on the domain border. Don't know how the closest point
    # on domain boundary is found.
    extrapolate(fxy_inter, Flat())
end


text:
# https://oracle.github.io/odpi/doc/user_guide/data_types.html
# https://docs.oracle.com/cd/B28359_01/server.111/b28318/datatype.htm#CNCPT012
# https://docs.oracle.com/cd/B28359_01/server.111/b28320/limits001.htm

const SIZEOF_ORA_DATA = sizeof(Oracle.OraData)

@inline function oracle_type(v::AbstractOracleValue{O,N}) :: OraOracleTypeNum where {O,N}
    return O
end

@inline function native_type(v::AbstractOracleValue{O,N}) :: OraNativeTypeNum where {O,N}
    return N
end

@inline use_add_ref(::AbstractOracleValue) = false
@inline use_add_ref(v::ExternOracleValue) = v.use_add_ref

@inline parent(v::ExternOracleValue) = v.parent
@inline parent(v::JuliaOracleValue) = error("Not implemented.")

@inline function get_data_handle(val::ExternOracleValue, offset::Integer) :: Ptr{OraData}
    @assert offset > 0 "Invalid offset $offset."
    return val.data_handle + (offset-1)*SIZEOF_ORA_DATA
end

@inline function get_data_handle(val::ExternOracleValue) :: Ptr{OraData}
    # equivalent of get_data_handle(val, 0)
    return val.data_handle
end

@inline function get_data_handle(val::JuliaOracleValue, offset::Integer) :: Ptr{OraData}
    @assert offset > 0 "Invalid offset $offset."
    return pointer(val.buffer, offset)
end

@inline function get_data_handle(val::JuliaOracleValue) :: Ptr{OraData}
    # equivalent of get_data_handle(val, 0)
    return pointer(val.buffer)
end

@inline parse_oracle_value(val::AbstractOracleValue, offset::Integer) = parse_oracle_value_at(val, get_data_handle(val, offset))
@inline parse_oracle_value(val::AbstractOracleValue) = parse_oracle_value_at(val, get_data_handle(val))

@inline is_null(ptr::Ptr{OraData}) = Bool(dpiData_getIsNull(ptr))

@inline function parse_oracle_value_at(val::AbstractOracleValue, data_handle::Ptr{OraData})
    if is_null(data_handle)
        return missing
    else
        parse_non_null_oracle_value_at(val, data_handle)
    end
end

# O -> OraOracleTypeNum, N -> OraNativeTypeNum
@generated function parse_non_null_oracle_value_at(val::AbstractOracleValue{O, N}, data_handle::Ptr{OraData}) where {O, N}

    @assert isa(O, OraOracleTypeNum)
    @assert isa(N, OraNativeTypeNum)

    if O == ORA_ORACLE_TYPE_NATIVE_DOUBLE
        @assert N == ORA_NATIVE_TYPE_DOUBLE

        return quote
            dpiData_getDouble(data_handle)
        end
    end

    if O == ORA_ORACLE_TYPE_BOOLEAN
        @assert N == ORA_NATIVE_TYPE_BOOLEAN

        return quote
            Bool(dpiData_getBool(data_handle))
        end
    end

    if O == ORA_ORACLE_TYPE_NATIVE_FLOAT
        @assert N == ORA_NATIVE_TYPE_FLOAT

        return quote
            dpiData_getFloat(data_handle)
        end
    end

    if O == ORA_ORACLE_TYPE_NUMBER

        # DPI_NATIVE_TYPE_DOUBLE, DPI_NATIVE_TYPE_BYTES, DPI_NATIVE_TYPE_INT64, DPI_NATIVE_TYPE_UINT64

        if N == ORA_NATIVE_TYPE_DOUBLE
            return quote
                dpiData_getDouble(data_handle)
            end

        elseif N == ORA_NATIVE_TYPE_BYTES
            error("Numeric values with native type BYTES not supported.")

        elseif N == ORA_NATIVE_TYPE_INT64
            return quote
                dpiData_getInt64(data_handle)
            end

        elseif N == ORA_NATIVE_TYPE_UINT64
            return quote
                dpiData_getUint64(data_handle)
            end
# see issue #21
#        elseif N == ORA_NATIVE_TYPE_NUMBER
#            return quote
#                dpiData_getNumber(data_handle)
#            end
        else
            error("Native type $N not expected for value with numeric oracle type.")
        end
    end

    if N == ORA_NATIVE_TYPE_BYTES

        if O == ORA_ORACLE_TYPE_RAW || O == ORA_ORACLE_TYPE_LONG_RAW
            # binary data
            return quote
                ptr_bytes = dpiData_getBytes(data_handle) # get a Ptr{OraBytes}
                ora_bytes = unsafe_load(ptr_bytes) # get a OraBytes
                unsafe_array = unsafe_wrap(Vector{UInt8}, ora_bytes.ptr, ora_bytes.length)
                return copy(unsafe_array)
            end
        else
            # character data
            return quote
                ptr_bytes = dpiData_getBytes(data_handle) # get a Ptr{OraBytes}
                ora_bytes = unsafe_load(ptr_bytes) # get a OraBytes
                return unsafe_string(ora_bytes.ptr, ora_bytes.length)
            end
        end
    end

    if O == ORA_ORACLE_TYPE_DATE
        @assert N == ORA_NATIVE_TYPE_TIMESTAMP "Invalid combination: [ $O, $N ]."

        return quote
            ptr_native_timestamp = dpiData_getTimestamp(data_handle)
            local ts::OraTimestamp = unsafe_load(ptr_native_timestamp)
            return parse_datetime(ts)
        end
    end

    if O == ORA_ORACLE_TYPE_TIMESTAMP
        if N == ORA_NATIVE_TYPE_TIMESTAMP
            return quote
                ptr_native_timestamp = dpiData_getTimestamp(data_handle)
                local ts::OraTimestamp = unsafe_load(ptr_native_timestamp)
                return parse_timestamp(ts)
            end
        elseif N == ORA_NATIVE_TYPE_DOUBLE
            error("Not implemented.")
        else
            error("Invalid combination: [ $O, $N ].")
        end
    end

    if O == ORA_ORACLE_TYPE_TIMESTAMP_TZ
        @assert N == ORA_NATIVE_TYPE_TIMESTAMP "Invalid combination: [ $O, $N ]."

        return quote
            ptr_native_timestamp = dpiData_getTimestamp(data_handle)
            local ts::OraTimestamp = unsafe_load(ptr_native_timestamp)
            return parse_timestamp_tz(ts)
        end
    end

    if O == O == ORA_ORACLE_TYPE_TIMESTAMP_LTZ
        @assert N == ORA_NATIVE_TYPE_TIMESTAMP "Invalid combination: [ $O, $N ]."

        return quote
            ptr_native_timestamp = dpiData_getTimestamp(data_handle)
            local ts::OraTimestamp = unsafe_load(ptr_native_timestamp)
            return parse_timestamp_ltz(ts)
        end
    end

    if N == ORA_NATIVE_TYPE_LOB
        return quote
            ptr_native_lob = dpiData_getLOB(data_handle)
            return Lob(parent(val), ptr_native_lob, O; use_add_ref=use_add_ref(val))
        end
    end

    if O == ORA_ORACLE_TYPE_NATIVE_INT
        @assert N == ORA_NATIVE_TYPE_INT64

        return quote
            dpiData_getInt64(data_handle)
        end
    end

    if O == ORA_ORACLE_TYPE_NATIVE_UINT
        @assert N == ORA_NATIVE_TYPE_UINT64

        return quote
            dpiData_getUint64(data_handle)
        end
    end

    error("Couldn't parse value for oracle type $O, native type $N.")
end

@inline Base.getindex(val::AbstractOracleValue) = parse_oracle_value(val)
@inline Base.getindex(val::AbstractOracleValue, offset::Integer) = parse_oracle_value(val, offset)
@inline Base.setindex!(oracle_value::AbstractOracleValue, value, offset::Integer=1) = set_oracle_value!(oracle_value, value, offset)

@inline function set_oracle_value!(oracle_value::JuliaOracleValue{O,N,T}, val::T, offset::Integer=1) where {O,N,T}
    @assert offset > 0 "Invalid offset $offset."
    oracle_value.buffer[offset] = val
    nothing
end

@generated function set_oracle_value_at!(oracle_value::AbstractOracleValue{O,N}, val::T, at::Ref{OraData}) where {O,N,T}

    if val <: Missing
        return quote
            dpiData_setNull(at)
        end
    end

    if N == ORA_NATIVE_TYPE_BYTES
        if val <: String
            return quote
                dpiData_setBytes(at, val)
            end
        end

        if val <: Vector{UInt8}
            return quote
                dpiData_setBytes(at, pointer(val), UInt32(length(val)))
            end
        end
    end

    if N == ORA_NATIVE_TYPE_BOOLEAN
        return quote
            dpiData_setBool(at, val)
        end
    end

    if N == ORA_NATIVE_TYPE_DOUBLE
        return quote
            dpiData_setDouble(at, val)
        end
    end

    if N == ORA_NATIVE_TYPE_INT64
        return quote
            dpiData_setInt64(at, val)
        end
    end

    if N == ORA_NATIVE_TYPE_UINT64
        return quote
            dpiData_setUint64(at, val)
        end
    end

    if N == ORA_NATIVE_TYPE_TIMESTAMP

        # DPI_ORACLE_TYPE_DATE, DPI_ORACLE_TYPE_TIMESTAMP, DPI_ORACLE_TYPE_TIMESTAMP_LTZ, DPI_ORACLE_TYPE_TIMESTAMP_TZ

        if O == ORA_ORACLE_TYPE_DATE
            @assert val <: Date
        elseif O == ORA_ORACLE_TYPE_TIMESTAMP
            @assert (val <: Date) || (val <: DateTime) || (val <: Timestamp)
        elseif O == ORA_ORACLE_TYPE_TIMESTAMP_TZ
            @assert val <: TimestampTZ{false}
        elseif O == ORA_ORACLE_TYPE_TIMESTAMP_LTZ
            @assert val <: TimestampTZ{true}
        else
            error("Oracle type $O not supported.")
        end

        return quote
            ts = OraTimestamp(val)
            dpiData_setTimestamp(at, ts)
        end
    end

    if N == ORA_NATIVE_TYPE_LOB
        @assert val <: Lob "Value must be of type `Oracle.Lob`."
        return quote
            dpiData_setLOB(at, val.handle)
        end
    end

# see issue #21
#    if N == ORA_NATIVE_TYPE_NUMBER
#        @assert val <: OraNumber "Value must be of type `Oracle.OraNumber`."
#        return quote
#            dpiData_setNumber(at, val)
#        end
#    end

    error("Setting values to AbstractOracleValue{$O, $N} is not supported.")
end

encoding(ora_string::OraBytes) = unsafe_string(ora_string.encoding)

function encoding(ora_string_ptr::Ptr{OraBytes})
    ora_string = unsafe_load(ora_string_ptr)
    return encoding(ora_string)
end

@generated function encoding(val::AbstractOracleValue{O,N}, offset::Integer=1) where {O,N}
    @assert offset > 0 "Invalid offset $offset."
    @assert N == ORA_NATIVE_TYPE_BYTES "Native type must be Oracle.ORA_NATIVE_TYPE_BYTES. Found: $N."

    return quote
        ptr_bytes = dpiData_getBytes(get_data_handle(val, offset))
        return encoding(ptr_bytes)
    end
end

#
# Get/Set Implementation for JuliaOracleValue
#

@inline function parse_oracle_value(oracle_value::JuliaOracleValue{O,N,T}, offset::Integer=1) where {O,N,T}
    @assert offset > 0 "Invalid offset $offset."
    return oracle_value.buffer[offset]
end

@inline function set_oracle_value!(oracle_value::AbstractOracleValue{O,N}, val::T, offset::Integer) where {O,N,T}
    return set_oracle_value_at!(oracle_value, val, get_data_handle(oracle_value, offset))
end

struct OracleTypeTuple
    oracle_type::OraOracleTypeNum
    native_type::OraNativeTypeNum
end

# accept julia types as arguments
@inline infer_oracle_type_tuple(::Type{Bool}) = OracleTypeTuple(ORA_ORACLE_TYPE_BOOLEAN, ORA_NATIVE_TYPE_BOOLEAN)
@inline infer_oracle_type_tuple(::Type{Float64}) = OracleTypeTuple(ORA_ORACLE_TYPE_NATIVE_DOUBLE, ORA_NATIVE_TYPE_DOUBLE)
@inline infer_oracle_type_tuple(::Type{Int64}) = OracleTypeTuple(ORA_ORACLE_TYPE_NATIVE_INT, ORA_NATIVE_TYPE_INT64)
@inline infer_oracle_type_tuple(::Type{UInt64}) = OracleTypeTuple(ORA_ORACLE_TYPE_NATIVE_UINT, ORA_NATIVE_TYPE_UINT64)
@inline infer_oracle_type_tuple(::Type{Date}) = OracleTypeTuple(ORA_ORACLE_TYPE_DATE, ORA_NATIVE_TYPE_TIMESTAMP)
@inline infer_oracle_type_tuple(::Type{DateTime}) = infer_oracle_type_tuple(Timestamp)
@inline infer_oracle_type_tuple(::Type{Timestamp}) = OracleTypeTuple(ORA_ORACLE_TYPE_TIMESTAMP, ORA_NATIVE_TYPE_TIMESTAMP)

# see issue #21
#@inline infer_oracle_type_tuple(::Type{OraNumber}) = OracleTypeTuple(ORA_ORACLE_TYPE_NUMBER, ORA_NATIVE_TYPE_NUMBER)

# accept julia values as arguments
for type_sym in (:Bool, :Float64, :Int64, :UInt64, :Date, :DateTime, :Timestamp)
    @eval begin
        @inline infer_oracle_type_tuple(::$type_sym) = infer_oracle_type_tuple($type_sym)
    end
end

@inline infer_oracle_type_tuple(::Type{Lob{O,P}}) where {O,P} = OracleTypeTuple(O, ORA_NATIVE_TYPE_LOB)
@inline infer_oracle_type_tuple(::Lob{O,P}) where {O,P} = OracleTypeTuple(O, ORA_NATIVE_TYPE_LOB)

@inline infer_oracle_type_tuple(::Type{TimestampTZ{false}}) = OracleTypeTuple(ORA_ORACLE_TYPE_TIMESTAMP_TZ, ORA_NATIVE_TYPE_TIMESTAMP)
@inline infer_oracle_type_tuple(::TimestampTZ{false}) = OracleTypeTuple(ORA_ORACLE_TYPE_TIMESTAMP_TZ, ORA_NATIVE_TYPE_TIMESTAMP)
@inline infer_oracle_type_tuple(::Type{TimestampTZ{true}}) = OracleTypeTuple(ORA_ORACLE_TYPE_TIMESTAMP_LTZ, ORA_NATIVE_TYPE_TIMESTAMP)
@inline infer_oracle_type_tuple(::TimestampTZ{true}) = OracleTypeTuple(ORA_ORACLE_TYPE_TIMESTAMP_LTZ, ORA_NATIVE_TYPE_TIMESTAMP)

@inline function infer_oracle_type_tuple(s::String)
    # max VARCHAR2 size is 4000 bytes
    if sizeof(s) <= 4000
        return OracleTypeTuple(ORA_ORACLE_TYPE_NVARCHAR, ORA_NATIVE_TYPE_BYTES)
    else
        return OracleTypeTuple(ORA_ORACLE_TYPE_NCLOB, ORA_NATIVE_TYPE_LOB)
    end
end

@inline function infer_oracle_type_tuple(::Type{String})
    # without information about string length, will best guess as a NVARCHAR
    return OracleTypeTuple(ORA_ORACLE_TYPE_NVARCHAR, ORA_NATIVE_TYPE_BYTES)
end

@inline function infer_oracle_type_tuple(b::Vector{UInt8})
    if length(b) <= 2000
        return OracleTypeTuple(ORA_ORACLE_TYPE_RAW, ORA_NATIVE_TYPE_BYTES)
    else
        return OracleTypeTuple(ORA_ORACLE_TYPE_LONG_RAW, ORA_NATIVE_TYPE_BYTES)
    end
end

@inline function infer_oracle_type_tuple(::Type{Vector{UInt8}})
    # without information about string length, will best guess as a RAW
    return OracleTypeTuple(ORA_ORACLE_TYPE_RAW, ORA_NATIVE_TYPE_BYTES)
end

function JuliaOracleValue(scalar::T) where {T}
    ott = infer_oracle_type_tuple(scalar)
    val = JuliaOracleValue(ott.oracle_type, ott.native_type, T)
    val[] = scalar
    return val
end