text:
__precompile__()

module KNITRO

    fn = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
    if isfile(fn)
        include(fn)
    else
        error("KNITRO not properly installed. Please run `] build KNITRO`")
    end

    using Libdl, SparseArrays
    import Base: show

    "A macro to make calling KNITRO's C API a little cleaner"
    macro ktr_ccall(func, args...)
        f = Base.Meta.quot(Symbol("KTR_$(func)"))
        args = [esc(a) for a in args]
        quote
            ccall(($f,libknitro), $(args...))
        end
    end
    "Load KNITRO version number via KTR API."
    function unsafe_get_release()
        len = 15
        out = zeros(Cchar,len)
        @ktr_ccall(get_release, Any, (Cint, Ptr{Cchar}), len, out)
        res = String(strip(String(convert(Vector{UInt8},out)), '\0'))
        return VersionNumber(split(res, " ")[2])
    end

    const KNITRO_VERSION = unsafe_get_release()

    # Wrapper of old API (soon deprecated)
    include("ktr_model.jl")
    include("ktr_callbacks.jl")
    include("ktr_functions.jl")
    include("ktr_defines.jl")
    include("ktr_params.jl")
    # wrapper with MathProgBase (only work with old API)
    include("MPBWrapper.jl")

    # Wrapper of new API (KNITRO's version > 11.0)
    # We load the new API only if KNITRO has correct version
    if KNITRO_VERSION >= v"11.0"
        include("kn_common.jl")
        include("kn_env.jl")

        include("kn_model.jl")
        include("kn_defines.jl")
        include("kn_params.jl")
        include("kn_variables.jl")
        include("kn_attributes.jl")
        include("kn_constraints.jl")
        include("kn_residuals.jl")
        include("kn_solve.jl")
        include("kn_callbacks.jl")

        # the MathOptInterface wrapper works only with the new API
        include("MOI_wrapper.jl")
    end
end


text:
using POMDPs
using POMDPToolbox
using Powseeker
using Plots

function plot_ro(policy::Policy)
    n = 100
    xs = linspace(-4000, 4000, 100)
    ys = linspace(-4000, 4000, 100)

    contour(xs, ys, (x,y)->simxy(policy, x, y), fill=true, title="$(typeof(policy))", zlim=(0,40))
end

function simxy(policy, x, y)
    mdp = PowseekerMDP()
    sim = RolloutSimulator(max_steps=mdp.duration, rng=MersenneTwister(1))
    s = SkierState(1, [x, y], deg2rad(90))
    return simulate(sim, mdp, policy, s)
end

mdp = PowseekerMDP()
plots = []
push!(plots, plot_ro(SkiOver(mdp)))
push!(plots, plot_ro(Downhill(mdp)))
push!(plots, plot_ro(Topout(mdp, 0.1)))
plot(plots...)
gui()

ext: jl


text:
using ArnoldiMethod, LinearAlgebra, NumericalIntegration, FFTW

h = LP.honeycomb() |>
    hamiltonian(hopping(1) + onsite(0., sublats = :A) - onsite(0., sublats = :B)) |>
    unitcell(region = RegionPresets.circle(20))

a = 0.1
eigs_dacp = DACPdiagonaliser(h, a,  maxdeg = 2, store_basis = false, d = 20)
eigs_dacp = DACPdiagonaliser(h, a,  maxdeg = 2, store_basis = true, d = 20)
eigs_dacp = DACPdiagonaliser(h, a,  maxdeg = 1, store_basis = true)
eigs_dacp = DACPdiagonaliser(h, a,  maxdeg = 1, store_basis = false)

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


text:
"""
    Value function approximation
"""
struct Value{N,NN}
    gradient::Vector{N}
    hessian::Vector{NN}
end

"""
    Action-value function approximation
"""
struct ActionValue{N,M,NN,MM,MN}
    gradient_state::Vector{N}
    gradient_action::Vector{M}
    hessian_state_state::Vector{NN}
    hessian_action_action::Vector{MM}
    hessian_action_state::Vector{MN}
end

"""
    Policy Data
"""
struct PolicyData{N,M,NN,MM,MN,NNN,MNN}
    # policy u = ū + K * (x - x̄) + k
    K::Vector{MN}
    k::Vector{M}

    K_candidate::Vector{MN}
    k_candidate::Vector{M}

    # value function approximation
    value::Value{N,NN}

    # action-value function approximation
    action_value::ActionValue{N,M,NN,MM,MN}

    # pre-allocated memory
	xx̂_tmp::Vector{NNN}
	ux̂_tmp::Vector{MNN}
	uu_tmp::Vector{MM}
	ux_tmp::Vector{MN}
end

function policy_data(dynamics::Vector{Dynamics{T}}) where T
    # policy
	K = [zeros(d.num_action, d.num_state) for d in dynamics]
    k = [zeros(d.num_action) for d in dynamics]

    K_candidate = [zeros(d.num_action, d.num_state) for d in dynamics]
    k_candidate = [zeros(d.num_action) for d in dynamics]

    # value function approximation
    P = [[zeros(d.num_state, d.num_state) for d in dynamics]...,
            zeros(dynamics[end].num_next_state, dynamics[end].num_next_state)]
    p =  [[zeros(d.num_state) for d in dynamics]...,
            zeros(dynamics[end].num_next_state)]

    value = Value(p, P)

    # action-value function approximation
    Qx = [zeros(d.num_state) for d in dynamics]
    Qu = [zeros(d.num_action) for d in dynamics]
    Qxx = [zeros(d.num_state, d.num_state) for d in dynamics]
    Quu = [zeros(d.num_action, d.num_action) for d in dynamics]
    Qux = [zeros(d.num_action, d.num_state) for d in dynamics]

    action_value = ActionValue(Qx, Qu, Qxx, Quu, Qux)

	xx̂_tmp = [zeros(d.num_state, d.num_next_state) for d in dynamics]
	ux̂_tmp = [zeros(d.num_action, d.num_next_state) for d in dynamics]
	uu_tmp = [zeros(d.num_action, d.num_action) for d in dynamics]
	ux_tmp = [zeros(d.num_action, d.num_state) for d in dynamics]

    PolicyData(K, k, K_candidate, k_candidate,
        value,
        action_value,
        xx̂_tmp, ux̂_tmp, uu_tmp, ux_tmp)
end


text:
# Replicate Huggett (1993)

# --------------------------------------------------------------------------- #
# Set parameters:
β = 0.9932
BoundSet = (-2.0, -4.0, -6.0, -8.0)

# Exogenous state variable grid
eH = 1.0
eL = 0.1
E = [eL; eH]

# Transition matrix: Pr(j|i) = Π(i,j) (columns sum to 1)
πHH = 0.925 # e_H | e_H
πHL = 0.5   # e_H | e_L
Π = [1 - πHL  πHL
    1 - πHH πHH]

# Guess initial price
q0 = 1.1

# --------------------------------------------------------------------------- #
# Solve for σ = 1.5
σ = 1.5

# Start loop over lower bound
for a_lb in BoundSet

    # Generate grid points for endogenous state variable
    a_max = 5
    a_size = 350
    A = range(a_lb, a_max, length = a_size)
    A = [A; ];

    (q, policy, Ψ) = step3(A, E, Π, q0, β, σ)

    println("- Solving the model for σ = $σ and lower bound of a = $a_lb.")
    println("  The optimal price is q = $q.")
    println("  ")

end

# --------------------------------------------------------------------------- #
# Solve for σ = 3
σ = 3

# Start loop over lower bound
for a_lb in BoundSet

    # Generate grid points for endogenous state variable
    a_max = 5
    a_size = 350
    A = range(a_lb, a_max, length = a_size)
    A = [A; ];

    (q, policy, Ψ) = step3(A, E, Π, q0, β, σ)

    println("- Solving the model for σ = $σ and lower bound of a = $a_lb.")
    println("  The optimal price is q = $q.")
    println("  ")

end

# --------------------------------------------------------------------------- #
# Graphs
q0 = 1.1
σ = 1.5
a_lb = -2.0
a_max = 5
a_size = 350
A = range(a_lb, a_max, length = a_size)
A = [A; ];

(q, policy, Ψ) = step3(A, E, Π, q0, β, σ)

# Find endogenous a upper bar (use aUB+20 in graphs to see the intersection)
aUB = findfirst(policy[:,2] - A[:] .<= 0.001)

# Unconstrained plot for Policy Function
p1 = plot(A, [policy A], linewidth = 2, title = "Policy function",
    label = ["eL" "eH" "45°"], legend =:topleft)

# Constrained  plot for Policy Function
p1B = plot(A[1:aUB+20], [policy[1:aUB+20,:] A[1:aUB+20]], linewidth = 2,
    title = "Policy function", label = ["eL" "eH" "45°"], legend =:topleft)

# Unconstrained plot for Stationary Distribution Ψ
p2 = plot(A, Ψ, linewidth = 2, title = "Stationary distribution",
    label = ["eL" "eH"], legend =:topleft)

# Constrained plot for Stationary Distribution Ψ
p2B = plot(A[1:aUB+20], Ψ[1:aUB+20,:], linewidth = 2,
    title = "Stationary distribution", label = ["eL" "eH"], legend =:topleft)

# Find rescaling value
scale = sum(Ψ[end, :])
F = Ψ ./ scale

# Unconstrained plot for Stationary Distribution Ψ
p3 = plot(A, F, linewidth = 2, title = "Stationary distribution",
    label = ["eL" "eH"], legend =:topleft)

# Constrained plot for Stationary Distribution Ψ
p3B = plot(A[1:aUB+20], F[1:aUB+20], linewidth = 2,
    title = "Stationary distribution", label = ["eL" "eH"], legend =:topleft)

include("Step3.jl")
include("Step2.jl")
include("Step1Grid.jl")


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module MultiplicativeInverses

import Base: div, divrem, rem, unsigned
using  Base: IndexLinear, IndexCartesian, tail
export multiplicativeinverse

unsigned(::Type{Int8}) = UInt8
unsigned(::Type{Int16}) = UInt16
unsigned(::Type{Int32}) = UInt32
unsigned(::Type{Int64}) = UInt64
unsigned(::Type{Int128}) = UInt128
unsigned(::Type{T}) where {T<:Unsigned} = T

abstract type  MultiplicativeInverse{T} end

# Computes integer division by a constant using multiply, add, and bitshift.

# The idea here is to compute floor(n/d) as floor(m*n/2^p) and then
# implement division by 2^p as a right bitshift.  The trick is finding
# m (the "magic number") and p. Roughly speaking, one can think of this as
#        floor(n/d) = floor((n/2^p) * (2^p/d))
# so that m is effectively 2^p/d.
#
# A few examples are illustrative:
# Division of Int32 by 3:
#   floor((2^32+2)/3 * n/2^32) = floor(n/3 + 2n/(3*2^32))
# The correction term, 2n/(3*2^32), is strictly less than 1/3 for any
# nonnegative n::Int32, so this divides any nonnegative Int32 by 3.
# (When n < 0, we add 1, and one can show that this computes
# ceil(n/d) = -floor(abs(n)/d).)
#
# Division of Int32 by 5 uses a magic number (2^33+3)/5 and then
# right-shifts by 33 rather than 32. Consequently, the size of the
# shift depends on the specific denominator.
#
# Division of Int32 by 7 would be problematic, because a viable magic
# number of (2^34+5)/7 is too big to represent as an Int32 (the
# unsigned representation needs 32 bits). We can exploit wrap-around
# and use (2^34+5)/7 - 2^32 (an Int32 < 0), and then correct the
# 64-bit product with an add (the `addmul` field below).
#
# Further details can be found in Hacker's Delight, Chapter 10.

struct SignedMultiplicativeInverse{T<:Signed} <: MultiplicativeInverse{T}
    divisor::T
    multiplier::T
    addmul::Int8
    shift::UInt8

    function SignedMultiplicativeInverse{T}(d::T) where T<:Signed
        d == 0 && throw(ArgumentError("cannot compute magic for d == $d"))
        signedmin = unsigned(typemin(T))
        UT = unsigned(T)

        # Algorithm from Hacker's Delight, section 10-4
        ad = unsigned(abs(d))
        t = signedmin + signbit(d)
        anc = t - one(UT) - rem(t, ad)   # absolute value of nc
        p = sizeof(d)*8 - 1
        q1, r1 = divrem(signedmin, anc)
        q2, r2 = divrem(signedmin, ad)
        while true
            p += 1                       # loop until we find a satisfactory p
            # update q1, r1 = divrem(2^p, abs(nc))
            q1 = q1<<1
            r1 = r1<<1
            if r1 >= anc                 # must be unsigned comparison
                q1 += one(UT)
                r1 -= anc
            end
            # update q2, r2 = divrem(2^p, abs(d))
            q2 = q2<<1
            r2 = r2<<1
            if r2 >= ad
                q2 += one(UT)
                r2 -= ad
            end
            delta = ad - r2
            (q1 < delta || (q1 == delta && r1 == 0)) || break
        end

        m = flipsign((q2 + one(UT)) % T, d)  # resulting magic number
        s = p - sizeof(d)*8                  # resulting shift
        new(d, m, d > 0 && m < 0 ? Int8(1) : d < 0 && m > 0 ? Int8(-1) : Int8(0), UInt8(s))
    end
end
SignedMultiplicativeInverse(x::Signed) = SignedMultiplicativeInverse{typeof(x)}(x)

struct UnsignedMultiplicativeInverse{T<:Unsigned} <: MultiplicativeInverse{T}
    divisor::T
    multiplier::T
    add::Bool
    shift::UInt8

    function UnsignedMultiplicativeInverse{T}(d::T) where T<:Unsigned
        d == 0 && throw(ArgumentError("cannot compute magic for d == $d"))
        u2 = convert(T, 2)
        add = false
        signedmin = one(d) << (sizeof(d)*8-1)
        signedmax = signedmin - one(T)
        allones = (zero(d) - 1) % T

        nc = allones - rem(convert(T, allones - d), d)
        p = 8*sizeof(d) - 1
        q1, r1 = divrem(signedmin, nc)
        q2, r2 = divrem(signedmax, d)
        while true
            p += 1
            if r1 >= convert(T, nc - r1)
                q1 = q1 + q1 + one(T)
                r1 = r1 + r1 - nc
            else
                q1 = q1 + q1
                r1 = r1 + r1
            end
            if convert(T, r2 + one(T)) >= convert(T, d - r2)
                add |= q2 >= signedmax
                q2 = q2 + q2 + one(T)
                r2 = r2 + r2 + one(T) - d
            else
                add |= q2 >= signedmin
                q2 = q2 + q2
                r2 = r2 + r2 + one(T)
            end
            delta = d - one(T) - r2
            (p < sizeof(d)*16 && (q1 < delta || (q1 == delta && r1 == 0))) || break
        end
        m = q2 + one(T)              # resulting magic number
        s = p - sizeof(d)*8 - add    # resulting shift
        new(d, m, add, s % UInt8)
    end
end
UnsignedMultiplicativeInverse(x::Unsigned) = UnsignedMultiplicativeInverse{typeof(x)}(x)

function div(a::T, b::SignedMultiplicativeInverse{T}) where T
    x = ((widen(a)*b.multiplier) >>> sizeof(a)*8) % T
    x += (a*b.addmul) % T
    ifelse(abs(b.divisor) == 1, a*b.divisor, (signbit(x) + (x >> b.shift)) % T)
end
function div(a::T, b::UnsignedMultiplicativeInverse{T}) where T
    x = ((widen(a)*b.multiplier) >>> sizeof(a)*8) % T
    x = ifelse(b.add, convert(T, convert(T, (convert(T, a - x) >>> 1)) + x), x)
    ifelse(b.divisor == 1, a, x >>> b.shift)
end

rem(a::T, b::MultiplicativeInverse{T}) where {T} =
    a - div(a, b)*b.divisor

function divrem(a::T, b::MultiplicativeInverse{T}) where T
    d = div(a, b)
    (d, a - d*b.divisor)
end

multiplicativeinverse(x::Signed) = SignedMultiplicativeInverse(x)
multiplicativeinverse(x::Unsigned) = UnsignedMultiplicativeInverse(x)

end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
function mimoUplink(Ns::Ti,Mc::Ti, M::Array{Ti,1}, K::Array{Ti,1},
     Ki::Array{Ti,1}, Tt::Array{Ti,1}, Td::Ti, rho::Array{Tf,1},
     gamma::Array{Tf,1}=[Inf],delta::Tf=0.0) where {Ti<:Integer,Tf<:AbstractFloat}
#mimoUplink(Ns,Mc,M,K,Tt,Td,rho)
#  Simulate uncoded transmission of QAM signals over a flat-fading Gaussian
#  uplink channel, for a base station with multiple receive
#  antennas. Constelations of size Mc are used with M receive antennas and K
#  txmitters, each with one transmit antenna. A power of "rho" per user is
#  utilized. Tt is the number of training samples to use, Td is the number
#  of data symbols, with Tt+Td being the training period. The resulting
#  symbol error rate is averaged over all users and Ns training periods.
#mimoUplink(Ns,Mc,M,K,Tt,Td,rho,gamma=Inf,delta=0)
#  If gamma is present, it is a CIR. If delta is present, it is
#  the pilot boost in db.
#
# Examples show SER as we deviate away from LTE PUSCH defaults
#  # Performance as a function of SNR
#  mimoUplink(100,4,4,[4],[0],[12],12,[-5.0:5:20;])
#  # Performance as a function of CIR
#  mimoUplink(400,4,4,[2],[1],[12],12,[15],[-15.0:5:25],0)
#  # Performance as a function of the number of antennas
#  mimoUplink(1000,4,[2:8],[2],[0],[4],8,[10])
#  # Performance as a function of Tt
#  mimoUplink(400,4,4,[2],[0],[2:2:16],12,[10])
#  # Performance as a function of the number of users
#  mimoUplink(400,4,4,[1:4],[0],[12],12,[10])
# Add an additional 'gamma' argument to the end of any of these functions to
# simulate with a single interferer of CIR gamma.

# By Christian Peel  (chris.peel@ieee.org)
# Last Modified: Fri 14 Nov 14, 11:57am by cpeel

println("   Mc     M     K    Ki    Tt    Td   rho  gamma")

#algs = getAlgs()
algs = []

# It would be great to print out the alg names here
#for ax = 1:min(length(algs),4)
#    @printf("%10s ",algs[ax][2])
#end
#println

out = []

xscale = :identity;
yscale = :log10;
if length(rho)>1
  for s = 1:length(rho)
    push!(out,simQAM(Ns,Mc,M[1],K[1],Ki[1],Tt[1],Td,rho[s],gamma[1],delta,algs));
  end
  xval = rho;
#  xlab = L"\rho"* " (SNR in dB)";
  xlab = "rho (SNR in dB)";
  tstr = @sprintf("%dQAM,Ns=%d,M=%d,K=%d,Ki=%d,Tt=%d,Td=%d,CIR=%4.1f",
                   Mc,Ns,M[1],K[1],Ki[1],Tt[1],Td,gamma[1]);
elseif length(gamma)>1
  for s = 1:length(gamma)
    push!(out,simQAM(Ns,Mc,M[1],K[1],Ki[1],Tt[1],Td,rho[1],gamma[s],delta,algs));
  end
  xval = gamma;
  xlab = "gamma (CIR in dB)";
  tstr = @sprintf("%dQAM,Ns=%d,M=%d,K=%d,Ki=%d,Tt=%d,Td=%d,SNR=%4.1f",
                 Mc,Ns,M[1],K[1],Ki[1],Tt[1],Td,rho[1]);
elseif length(K)>1
  for s = 1:length(K)
    push!(out,simQAM(Ns,Mc,M[1],K[s],Ki[1],Tt[1],Td,rho[1],gamma[1],delta,algs));
  end
  xval = K;
  xlab = "K (# users)";
  tstr = @sprintf("%dQAM,Ns=%d,M=%d,Ki=%d,Tt=%d,Td=%d,SNR=%4.1f,CIR=%4.1f",
                 Mc,Ns,M[1],Ki[1],Tt[1],Td,rho[1],gamma[1]);
elseif length(Tt)>1
  for s = 1:length(Tt)
    push!(out,simQAM(Ns,Mc,M[1],K[1],Ki[1],Tt[s],Td,rho[1],gamma[1],delta,algs));
  end
  xval = Tt;
  xlab = "Tt (training length)";
  tstr = @sprintf("%dQAM,Ns=%d,M=%d,K=%d,Ki=%d,Td=%d,SNR=%4.1f,CIR=%4.1f",
                 Mc,Ns,M[1],K[1],Ki[1],Td,rho[1],gamma[1]);
elseif length(M)>1
  for s = 1:length(M)
    push!(out,simQAM(Ns,Mc,M[s],K[1],Ki[1],Tt[1],Td,rho[1],gamma[1],delta,algs));
  end
  xval = M;
  xlab = "M (# Rx antennas)";
  tstr = @sprintf("%dQAM,Ns=%d,K=%d,Ki=%d,Tt=%d,Td=%d,SNR=%4.1f,CIR=%4.1f",
                 Mc,Ns,K[1],Ki[1],Tt[1],Td,rho[1],gamma[1]);
elseif length(K)==length(M) && length(K)>1
  for s = 1:length(K)
    push!(out,simQAM(NsMc,M[s],K[s],Ki[1],Tt[1],Td,rho[1],gamma[1],delta,algs));
  end
  xval = K;
  xlab = "M==K (# users, # antennas)";
  tstr = @sprintf("%dQAM,Ns=%d,M=%d,Ki=%d,Tt=%d,Td=%d,SNR=%4.1f,CIR=%4.1f",
                 Mc,Ns,M[1],Ki[1],Tt[1],Td,rho[1],gamma[1]);
elseif length(Mc)>1
  for s = 1:length(Mc)
    push!(out,simQAM(Ns,Mc[s],M[1],K[1],Ki[1],Tt[1],Td,rho[1],gamma[1],delta,algs));
  end
  xval = log2(Mc);
  xlab = "# constellation bits";
  tstr = "";
else
  out = simQAM(Ns,Mc,M[1],K[1],Ki[1],Tt[1],Td,rho[1],gamma[1],delta,algs);
end

pColor = ["r>-","bo--","kx-.","gd-","c^--","m*-.",
          "rs--","gp-.","bv-","kh--","c+-.","m.-"];
pIdx   = 1;

Ns = size(out,1);

pltT = plot(legend=:topleft);

ser = zeros(Ns);
for a=1:length(out[1][2])
    for k=1:Ns
        ser[k] = out[k][1][a];
#        mi(k) = out(k).mi(a);
    end
    plot!(pltT,xval,ser,xscale=xscale,yscale=yscale,label=out[1][2][a],linewidth=3);
#    hold(true);
##     figure(2)
##     plot(xval,mi,pColor[pIdx]);  hold(true);
    if pIdx == length(pColor)
         pIdx = 1;
    else
         pIdx = pIdx + 1;
    end
end

plot!(pltT,xlabel=xlab,ylabel="SER")

    # grid(which="both",axis="y")
# grid(which="major",axis="x")
# title(tstr)

## figure(2)
## xlabel(xlab);
## ylabel("Mut Info");
## legend(out(1).name,4);
## grid on;
## title(tstr)
display(pltT)

return
end

#######################################################################
function simQAM(Ns,Mc,M,K,Ki,Tt,Td,rdb,gdb,ddb,algs)

delta   = 10 .^(ddb/10); # pilot boost
rho     = 10 .^(rdb/10);
gamma   = 10 .^(gdb/10);
if K>M
    error("Can't have K>M!!")
end
T = Tt+Td;
Ka = K+Ki;
#
# Modulation class, modulated bits
#
Mp = convert(Int,round(sqrt(Mc)));
# Information bits, signals, Threshold.
Zreal = rand((0:Mp-1),K+Ki,Td*Ns);
Zimag = rand((0:Mp-1),K+Ki,Td*Ns);
Zi = Zreal*Mp+Zimag;
ZZ = Zi[1:K,:];
C,d,Cr = scodes(Mc,"QAM");
gam = sqrt(2/3*(Mc-1));

#
# Channel coefs, noise.
#
sigmaSq = 1;
rho     = rho;
# sigmaSq = 1/rho;
# rho     = 1;
HH = (randn(M,K+Ki,Ns) + im*randn(M,K+Ki,Ns))*sqrt(1/2);# Channel
HH[:,1:K,:] = HH[:,1:K,:]*sqrt(rho);                    # Signal power
HH[:,K+1:Ka,:] = HH[:,K+1:Ka,:]*sqrt(rho/gamma);        # Interference power
NN = (randn(M,T*Ns) + im*randn(M,T*Ns))*sqrt(sigmaSq/2);

#
# Training signal
#
IK = Matrix{Float64}(I,K,K)
IM = Matrix{Float64}(I,M,M)
St = repeat(IK,1,round(Int,ceil(Tt/K)));
St = St[:,1:Tt]*sqrt(delta)
SSi = sign.(rand(Ki,Tt,Ns) .- .5)/sqrt(Tt);

onesKT = ones(K,Td);
ri = [1:K;];
ii = [(K+1):(2*K);];
ria = [1:K;]';
iia = [(Ka+1):(Ka+K);]';

Nalgs = 5;
ZD = zeros(K,Td*Ns,Nalgs);
algNames = Array{Any}(missing,Nalgs)

for ix = 1:Ns
    zix = (ix-1)*Td+1:ix*Td
    N = NN[:,(ix-1)*T+1:ix*T]
    Z = Zi[:,zix]
    Si = SSi[:,:,ix]
    Sti = [St; Si]
    Sd = C[Z.+1]
    S = [Sti Sd]
    H = HH[:,:,ix]

    Y = H*S+N; # Channel
    Yt = Y[:,1:Tt];
    Yd = Y[:,Tt+1:Tt+Td];
    Ydc = [real(Yd); imag(Yd)];
    #
    # Genie channel
    #
    Hdes = H[:,1:K];
    Hc = [real(Hdes) -imag(Hdes); imag(Hdes) real(Hdes)];
    HcAll = [real(H) -imag(H); imag(H) real(H)];
    #
    # Training
    #
    Hhat = Yt*St'*inv(St*St');
    Hhatc = [real(Hhat) -imag(Hhat); imag(Hhat) real(Hhat)];
    #
    # additional Training
    #
    Hhat2= Hhat;
    #
    What = (Yt-Hhat2*St);
    Rhat = What*What';
    var2hat = sum(diag(Rhat))/(Tt-1)/M;
    powhat = sum(diag(Hhat2*Hhat2'))/K/M;
    Hhat2 = Yt*St'*inv(St*St'+ I*var2hat/powhat*10);

    ax=0;

    ax = ax+1;
    algNames[ax] = "Zero Forcing";
    Hhat = Yt*St'*inv(St*St');
    Wls = Hhat/(Hhat'*Hhat);
    Zd  = mimo_slice(Wls'*Yd,onesKT,C);
    ZD[:,zix,ax] = Zd;

    # ax = ax+1;
    # algNames[ax]= "MMSE size"
    # Bw, Tw = sizereduction([Hhat; sigmaSq*I]);
    # Ti = inv(Tw);
    # Wlr = Bw/(Bw'*Bw);
    # Ztmp = Wlr'*[Yd; zeros(K,Td)];
    # Stmp = roundc((Ztmp*gam-(1+im)*Ti*onesKT)/2)*2 +(1+im)*Ti*onesKT;
    # Zd = mimo_slice(Tw*Stmp/gam,onesKT,C);
    # ZD[:,zix,ax] = Zd;

    # ax = ax+1;
    # algNames[ax]= "MMSE brun"
    # Bw, Tw = brun([Hhat; sigmaSq*I]);
    # Ti = inv(Tw);
    # Wlr = Bw/(Bw'*Bw);
    # Ztmp = Wlr'*[Yd; zeros(K,Td)];
    # Stmp = roundc((Ztmp*gam-(1+im)*Ti*onesKT)/2)*2 +(1+im)*Ti*onesKT;
    # Zd = mimo_slice(Tw*Stmp/gam,onesKT,C);
    # ZD[:,zix,ax] = Zd;

    ax = ax+1;
    algNames[ax]= "MMSE LLL"
    Bw, Tw = lll([Hhat; sigmaSq*I]);
    Ti = inv(Tw);
    Wlr = Bw/(Bw'*Bw);
    Ztmp = Wlr'*[Yd; zeros(K,Td)];
    Stmp = roundc((Ztmp*gam-(1+im)*Ti*onesKT)/2)*2 +(1+im)*Ti*onesKT;
    Zd = mimo_slice(Tw*Stmp/gam,onesKT,C);
    ZD[:,zix,ax] = Zd;

    # ax = ax+1;
    # algNames[ax]= "MMSE Brun"
    # Bw, Tw = brun([Hhat; sigmaSq*complex(IK)])
    # Ti = inv(Tw);
    # Wlr = Bw/(Bw'*Bw);
    # Ztmp = Wlr'*[Yd; zeros(K,Td)];
    # Stmp = roundc((Ztmp*gam-(1+im)*Ti*onesKT)/2)*2 +(1+im)*Ti*onesKT;
    # Zd = mimo_slice(Tw*Stmp/gam,onesKT,C);
    # ZD[:,zix,ax] = Zd;

    ax = ax+1;
    algNames[ax]= "MMSE Seysen"
    Bw, Tw = seysen([Hhat; sigmaSq*I]);
    Ti = inv(Tw);
    Wlr = Bw/(Bw'*Bw);
    Ztmp = Wlr'*[Yd; zeros(K,Td)];
    Stmp = roundc((Ztmp*gam-(1+im)*Ti*onesKT)/2)*2 +(1+im)*Ti*onesKT;
    Zd = mimo_slice(Tw*Stmp/gam,onesKT,C);
    ZD[:,zix,ax] = Zd;

    ax = ax+1;
    algNames[ax] = "MMSE-VBLAST"
    (W,P,B) = vblast(Hhat,1/sigmaSq/K);
    Yp = W*Yd;
    Sp = zeros(Complex{Float64},K,Td);
    for tt=1:Td
        for kx=1:K
            Sp[kx,tt] = siso_demod(Yp[kx,tt] -transpose(B[kx,:])*Sp[:,tt],C);
        end
    end
    ZD[:,zix,ax] = mimo_slice(P*Sp,onesKT,C);

    # ax = ax+1;
    # algNames[ax] = "IRC";
    # Hhat = (Yt*St')/(St*St');
    # Wp = Yt-Hhat*St;     # To be used to estimate the noise variance
    # Rhat = Wp*Wp'/(max(Tt-K,1)*M);
    # Ri = inv(Rhat + tr(Rhat)/1e14*complex(IM));
    # W = Ri*Hhat/(Hhat'*Ri*Hhat);
    # Zd = mimo_slice(W'*Yd,W'*Hhat*onesKT,C);
    # ZD[:,zix,ax] = Zd;

    ## doesn't work for interference cases
    ax = ax+1;
    algNames[ax] = "Sphere Decoder";
    z_est = hardsphere(Ydc,Hhatc/gam,Mp);
    Zd = z_est[ri,:]*Mp + z_est[ii,:]
    Zd = Int.((Zd.+Mc.+1)./2 .-1)
    ZD[:,zix,ax] = Zd;

end
@printf("%5d %5d %5d %5d %5d %5d %5.1f %5.1f",
          Mc,  M,  K, Ki, Tt, Td,  rdb,  gdb)

ser = zeros(Nalgs);
#Calculate symbol error rate
for ax = 1:Nalgs
    ser[ax] = mean(map(!=,ZZ,ZD[:,:,ax]));
    if ser[ax]==0
        ser[ax] = 1e-6
    end
    if ax<6
        @printf("    %7.4f",ser[ax])
    end
end
@printf("\n")

return ser, algNames
end
#######################################################################

max_stars_repo_name: christianpeel/MUMIMO

max_issues_repo_name: christianpeel/MUMIMO

max_forks_repo_name: christianpeel/MUMIMO


text:
# Note: we use the 100s place of the enum to represent the number of dimensions
# And the 10s place to represent which dimension is the observation (N), 0 for no observation dimension
@enum ImageFormat begin
    PNG=0;
    L=100;
    CL=200; LC; HW; WH; NL=210; LN=220;
    HWC=300; WHC; CHW; CWH; NCL=310; NHW; NWH; NLC; CLN=330; LCN; HWN; WHN;
    NHWC=410; NWHC; NCHW; NCWH; HWCN=440; WHCN; CHWN; CWHN;
end

expected_ndims(x::ImageFormat) = Int(x) ÷ 100 # get the 100's digit
obs_dim(x::ImageFormat) = (Int(x) % 100) ÷ 10 # get the 10's digit
const strip_obs = Dict(
    NL=>L, LN=>L, NCL=>CL, CLN=>CL, NLC=>LC, LCN=>LC,
    NHW=>HW, HWN=>HW, NWH=>WH, WHN=>WH,
    NHWC=>HWC, NWHC=>WHC, NCHW=>CHW, NCWH=>CWH, HWCN=>HWC, WHCN=>WHC, CHWN=>CHW, CWHN=>CWH
)

function convert_to_CHW(imgArray::AbstractArray, format)
    converted =
        format == L   ? reshape(imgArray, (1, 1, size(imgArray, 1))) :
        format == CL  ? reshape(imgArray, (size(imgArray, 1), 1, size(imgArray, 2))) :
        format == LC  ? reshape(transpose(imgArray), (size(imgArray, 2), 1, size(imgArray, 1))) :
        format == HW  ? reshape(imgArray, (1, size(imgArray, 1), size(imgArray, 2))) :
        format == WH  ? reshape(transpose(imgArray), (1, size(imgArray, 2), size(imgArray, 1))) :
        format == HWC ? permutedims(imgArray, (3, 1, 2)) :
        format == WHC ? permutedims(imgArray, (3, 2, 1)) :
        format == CHW ? imgArray :
        format == CWH ? permutedims(imgArray, (1, 3, 2)) :
        #== else ==# throw("Invalid format")
    return converted
end

function ColorType_from_nchannels(channels)
    color =
        channels == 1 ? Gray :
        channels == 2 ? GrayA :
        channels == 3 ? RGB :
        channels == 4 ? RGBA :
        #== else ==# throw("Too many channels")
    return color
end


text:
"""

    EarlyStopper(c...; verbosity=0)

Instantiate an object for tracking whether one or more stopping
criterion `c` apply, given a sequence of losses.

For a list of possible criterion, do
`subtypes(EarlyStopping.StoppingCriterion)`.

### Sample usage

    stopper = EarlyStopper(Patience(1), NotANumber())
    done!(stopper, 0.123) # false
    done!(stopper, 0.234) # true

    julia> message(stopper)
    "Early stop triggered by Patience(1) stopping criterion. "

### Training losses

For criteria tracking both an "out-of-sample" loss and a "training"
loss (eg, stopping criterion of type `PQ`), specify `training=true` if
the update is for training, as in

    done!(stopper, 0.123, training=true)

In these cases, the out-of-sample update must always come after the
corresponding training update. Multiple training updates may preceed
the out-of-sample update.

"""
mutable struct EarlyStopper{S}
    criterion::S
    verbosity::Int
    state
    EarlyStopper(criterion::S; verbosity=0) where S =
        new{S}(criterion, verbosity)
end
EarlyStopper(criteria...; kwargs...) = EarlyStopper(sum(criteria); kwargs...)

for f in [:message, :done]
    eval(quote
         $f(stopper::EarlyStopper) = $f(stopper.criterion, stopper.state)
         end)
end

# defines 2 private functons `done_after_update!(stopper, loss)` and
# `done_after_training_update!(stopper, loss)`:
for f in [:update, :update_training]
    newf = Symbol(string("done_after_", f, "!"))
    eval(quote
         function $newf(stopper::EarlyStopper, loss)
             if isdefined(stopper, :state)
                 stopper.state = $f(stopper.criterion, loss, stopper.state)
             else
                 stopper.state = $f(stopper.criterion, loss)
             end
             return done(stopper)
         end
         end)
end

"""
    done!(stopper::EarlyStopper)
"""
function done!(stopper::EarlyStopper, loss; training=false)
    ret = if training
        done_after_update_training!(stopper, loss)
    else
        done_after_update!(stopper, loss)
    end
    if stopper.verbosity > 0
        suffix = training ? "training " : ""
        loss_str = suffix*"loss"
        @info "$loss_str: $loss\t state: $(stopper.state)"
    end
    return ret
end


text:
module _vscodeserver

@static if VERSION < v"0.7.0-DEV.357"
    function remlineinfo!(x)
        if isa(x, Expr)
            id = find(map(x -> (isa(x, Expr) && x.head == :line) || (isdefined(:LineNumberNode) && x isa LineNumberNode), x.args))
            deleteat!(x.args, id)
            for j in x.args
                remlineinfo!(j)
            end
        end
        x
    end
else
    function remlineinfo!(x)
        if isa(x, Expr)
            if x.head == :macrocall && x.args[2] != nothing
                id = find(map(x -> (isa(x, Expr) && x.head == :line) || (isdefined(:LineNumberNode) && x isa LineNumberNode), x.args))
                deleteat!(x.args, id)
                for j in x.args
                    remlineinfo!(j)
                end
                insert!(x.args, 2, nothing)
            else
                id = find(map(x -> (isa(x, Expr) && x.head == :line) || (isdefined(:LineNumberNode) && x isa LineNumberNode), x.args))
                deleteat!(x.args, id)
                for j in x.args
                    remlineinfo!(j)
                end
            end
        end
        x
    end
end

import Base: display, redisplay
global active_module = :Main

struct InlineDisplay <: Display end

pid = Base.ARGS[1]
Base.ENV["JULIA_EDITOR"] = Base.ARGS[2]

function change_module(newmodule::String, print_change = true)
    global active_module
    smods = Symbol.(split(newmodule, "."))

    val = Main
    for i = 1:length(smods)
        if isdefined(val, smods[i])
            val = getfield(val, smods[i])
        else
            println("Could not find module $newmodule")
            return
        end
    end
    expr = parse(newmodule)
    active_module = expr
    repl = Base.active_repl
    main_mode = repl.interface.modes[1]
    main_mode.prompt = string(newmodule,"> ")
    main_mode.on_done = Base.REPL.respond(repl,main_mode; pass_empty = false) do line
        if !isempty(line)
            ex = parse(line)
            if ex isa Expr && ex.head == :module
                ret = :( eval($expr, Expr(:(=), :ans, Expr(:toplevel, parse($line)))) )
            else
                ret = :( eval($expr, Expr(:(=), :ans, parse($line))) )
            end
        else
            ret = :(  )
        end
        out = connect(to_vscode)
        write(out, string("repl/variables,", getVariables(), "\n"))
        close(out)
        return ret
    end
    print(" \r ")
    print_change && println("Changed root module to $expr")
    print_with_color(:green, string(newmodule,"> "), bold = true)
end

function get_available_modules(m=Main, out = Module[])
    info("here")
    for n in names(m, true, true)
        if isdefined(m, n) && getfield(m, n) isa Module  && !(getfield(m, n) in out)
            M = getfield(m, n)
            push!(out, M)
            get_available_modules(M, out)
        end
    end
    out
end

function getVariables()
    M = current_module()
    variables = []
    msg = ""
    for n in names(M)
        !isdefined(M, n) && continue
        x = getfield(M, n)
        t = typeof(x)
        msg = string(msg, ";", n, "::", t)
    end
    return msg
end

function generate_pipe_name(name)
    if is_windows()
        "\\\\.\\pipe\\vscode-language-julia-$name-$pid"
    elseif is_unix()
        joinpath(tempdir(), "vscode-language-julia-$name-$pid")
    end
end

!(is_unix() || is_windows()) && error("Unknown operating system.")

global_lock_socket_name = generate_pipe_name("terminal")
from_vscode = generate_pipe_name("torepl")
to_vscode = generate_pipe_name("fromrepl")

if issocket(from_vscode)
    rm(from_vscode)
end

@async begin
    server = listen(from_vscode)
    while true
        sock = accept(server)
        cmd = readline(sock)
        !startswith(cmd, "repl/") && continue
        text = readuntil(sock, "repl/endMessage")[1:end-15]
        if cmd == "repl/getAvailableModules"
            oSTDERR = STDERR
            redirect_stderr()
            ms = get_available_modules(current_module())
            redirect_stderr(oSTDERR)
            names = unique(sort(string.(ms)))
            out = connect(to_vscode)
            write(out, string("repl/returnModules,", join(names, ","), "\n"))
            close(out)
        elseif cmd == "repl/changeModule"
            change_module(strip(text, '\n'))
        elseif cmd == "repl/include"
            cmod = eval(active_module)
            ex = Expr(:call, :include, strip(text, '\n'))
            cmod.eval(ex)
        elseif cmd == "repl/getVariables"
            out = connect(to_vscode)
            write(out, getVariables())
            close(out)
        end
    end
end

conn = connect(global_lock_socket_name)

function display(d::InlineDisplay, ::MIME{Symbol("image/png")}, x)
    payload = stringmime(MIME("image/png"), x)
    print(conn, "image/png", ":", endof(payload), ";")
    print(conn, payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("image/png")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("image/svg+xml")}, x)
    payload = stringmime(MIME("image/svg+xml"), x)
    print(conn, "image/svg+xml", ":", endof(payload), ";")
    print(conn, payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("image/svg+xml")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("text/html")}, x)
    payload = stringmime(MIME("text/html"), x)
    print(conn, "text/html", ":", endof(payload), ";")
    print(conn, payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("text/html")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("juliavscode/html")}, x)
    payload = stringmime(MIME("juliavscode/html"), x)
    print(conn, "juliavscode/html", ":", endof(payload), ";")
    print(conn, payload)
end

Base.Multimedia.istextmime(::MIME{Symbol("juliavscode/html")}) = true

displayable(d::InlineDisplay, ::MIME{Symbol("juliavscode/html")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v2+json")}, x)
    payload = stringmime(MIME("application/vnd.vegalite.v2+json"), x)
    print(conn, "application/vnd.vegalite.v2+json", ":", endof(payload), ";")
    print(conn, payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v2+json")}) = true

function display(d::InlineDisplay, x)
    if mimewritable("application/vnd.vegalite.v2+json", x)
        display(d,"application/vnd.vegalite.v2+json", x)
    elseif mimewritable("juliavscode/html", x)
        display(d,"juliavscode/html", x)
    # elseif mimewritable("text/html", x)
    #     display(d,"text/html", x)
    elseif mimewritable("image/svg+xml", x)
        display(d,"image/svg+xml", x)
    elseif mimewritable("image/png", x)
        display(d,"image/png", x)
    else
        throw(MethodError(display,(d,x)))
    end

end

atreplinit(i->Base.Multimedia.pushdisplay(InlineDisplay()))
@async while true
    if isdefined(Base, :active_repl)
        print("\r")
        change_module("Main", false)
        break
    end
end
end

max_stars_repo_name: kindlychung/julia-vscode

max_issues_repo_name: kindlychung/julia-vscode

max_forks_repo_name: kindlychung/julia-vscode


text:
@testset "incomplete objectives" begin
    import NLSolversBase: df, fdf, make_f, make_df, make_fdf
    function f(x)
        sum(x->x^2,x)
    end

    g!(G, x) = copyto!(G, 2 .* x)
    g(x) = 2 .* x

    h!(H, x) = copy!(H, Diagonal(fill(2, length(n))))
    h(x) = Diagonal(fill(2, length(n)))

    function fg!(G, x)
        copyto!(G, 2 .* x)
        sum(x->x^2,x)
    end
    function just_fg!(F, G, x)
        !(G == nothing) && copyto!(G, 2 .* x)
        !(F == nothing) && sum(x->x^2,x)
    end
    fg(x) = f(x), g(x)
    function just_fgh!(F, G, H, x)
        !(H == nothing) && copy!(H, Diagonal(fill(2, length(n))))
        !(G == nothing) && copyto!(G, 2 .* x)
        !(F == nothing) && return sum(x->x^2,x)
    end
    fgh(x) = f(x), g(x), h(x)
    fdf!_real = only_fg!(just_fg!)
    fdf_real = only_fg(fg)

    function just_hv!(Hv, x, v)
        copyto!(Hv, 2 .* v)
    end
    function just_fghv!(F, G, Hv, x, v)
        G  === nothing || copyto!(G, 2 .* x)
        Hv === nothing || copyto!(Hv, 2 .* v)
        F  === nothing || return sum(x->x^2, x)
    end

    df_fdf_real = only_g_and_fg(g, fg)
    Random.seed!(3259)
    x = rand(10)
    G_cache = similar(x)
    G_cache2 = similar(G_cache)

    @test df(fdf!_real) === nothing
    @test df(fdf_real) === nothing
    @test df(df_fdf_real) === g
    @test df(df_fdf_real)(x) == g(x)

    @test fdf(fdf!_real) === just_fg!
    @test fdf(fdf_real) === fg
    @test df(df_fdf_real) == g
    @test fdf(df_fdf_real) == fg
    @test df(df_fdf_real)(x) == g(x)
    @test fdf(df_fdf_real)(x) == fg(x)

    for FDF in (fdf_real, fdf!_real)
        @test make_f(FDF, x, x[1])(x) == f(x)
        make_df(FDF, x, x[1])(G_cache, x)
        g!(G_cache2, x)
        @test G_cache == G_cache2
        f1 = make_fdf(FDF, x, x[1])(G_cache, x.*2)
        f2 = fg!(G_cache2, x.*2)
        @test G_cache == G_cache2
        @test f1 == f2
    end

    nd_fg = NonDifferentiable(only_fg(fg), x)
    nd_fg! = NonDifferentiable(only_fg!(just_fg!), x)
    for ND in (nd_fg, nd_fg!)
        value!(ND, x)
        value(ND) == f(x)
    end
    od_fg = OnceDifferentiable(only_fg(fg), x)
    od_fg! = OnceDifferentiable(only_fg!(just_fg!), x)
    od_fgh! = TwiceDifferentiable(only_fgh!(just_fgh!), x)
    _F = zero(eltype(x))
    od_fgh! = TwiceDifferentiable(only_fgh!(just_fgh!), x, _F)
    od_fgh! = TwiceDifferentiable(only_fgh!(just_fgh!), x, _F, similar(x))
    od_fgh! = TwiceDifferentiable(only_fgh!(just_fgh!), x, _F, similar(x), NLSolversBase.alloc_DF(x, _F))
#    od_fgh = TwiceDifferentiable(only_fgh(fgh), x)
    for OD in (od_fg, od_fg!, od_fgh!)#, od_fgh)
        value!(OD, x)
        @test value(OD) == f(x)
        gradient!(OD, x)
        @test gradient(OD) == g(x)
        @test gradient(OD, x) == g(x)
        value_gradient!(OD, 2 .* x)
        @test value(OD) == f(2 .* x)
        @test gradient(OD) == g(2 .* x)
    end

    # Incomplete TwiceDifferentiableHv
    v = randn(10)
    od_fg_and_hv = TwiceDifferentiableHV(only_fg_and_hv!(just_fg!, just_hv!), x)
    od_fghv      = TwiceDifferentiableHV(only_fghv!(just_fghv!), x)
    for OD in (od_fg_and_hv, od_fghv)
        gradient!(OD, x)
        @test gradient(OD) == g(x)
        hv_product!(OD, x, v)
        @test OD.Hv == 2v
        OD.f(x) == f(x)
        _g = similar(x)
        OD.fdf(_g, x)
        @test _g == g(x)
        @test OD.fdf(_g, x) == f(x)
    end
end
@testset "incomplete objectives vectors" begin
    function tf(x)
        x.^2
    end
    function tf(F, x)
        copyto!(F, tf(x))
    end

    tj!(J, x) = copyto!(J, Matrix(Diagonal(x)))
    tj(x) = Matrix(Diagonal(x))

    function tfj!(F, J, x)
        copyto!(J, Matrix(Diagonal(x)))
        copyto!(F, tf(x))
    end
    function just_tfj!(F, J, x)
        !(J == nothing) && copyto!(J, Matrix(Diagonal(x)))
        !(F == nothing) && copyto!(F, tf(x))
    end
    tfj(x) = tf(x), tj(x)

    fdf!_real = only_fj!(just_tfj!)
    fdf_real = only_fj(tfj)

    df_fdf_real = only_j_and_fj(tj, tfj)
    Random.seed!(3259)
    x = rand(10)
    J_cache = similar(Matrix(Diagonal(x)))
    J_cache2 = similar(Matrix(Diagonal(x)))
    F_cache = similar(x)
    F_cache2 = similar(x)

    @test df(fdf!_real) === nothing
    @test df(fdf_real) === nothing
    @test df(df_fdf_real) === tj
    @test df(df_fdf_real)(x) == tj(x)

    @test fdf(fdf!_real) === just_tfj!
    @test fdf(fdf_real) === tfj
    @test df(df_fdf_real) == tj
    @test fdf(df_fdf_real) == tfj
    @test df(df_fdf_real)(x) == tj(x)
    @test fdf(df_fdf_real)(x) == tfj(x)

    for FDF in (fdf_real, fdf!_real)
        @test make_f(FDF, x, x)(F_cache, x) == tf(x)
        make_df(FDF, x, x)(J_cache, x)
        tj!(J_cache2, x)
        @test J_cache == J_cache2
        make_fdf(FDF, x, x)(F_cache, J_cache, x.*2)
        tfj!(F_cache2, J_cache2, x.*2)
        @test F_cache == F_cache2
        @test J_cache == J_cache2
    end

    nd_fj = NonDifferentiable(only_fj(tfj), x, x)
    nd_fj! = NonDifferentiable(only_fj!(just_tfj!), x, x)
    for ND in (nd_fj, nd_fj!)
        value!(ND, x)
        value(ND) == tf(x)
    end
    od_fj = OnceDifferentiable(only_fj(tfj), x, x)
    od_fj! = OnceDifferentiable(only_fj!(just_tfj!), x, x)
    for OD in (od_fj, od_fj!)
        value!(OD, x)
        @test value(OD) == tf(x)
        jacobian!(OD, x)
        @test jacobian(OD) == tj(x)
        @test jacobian(OD, x) == tj(x)
        value_jacobian!(OD, 2 .* x)
        @test value(OD) == tf(2 .* x)
        @test jacobian(OD) == tj(2 .* x)
    end

end


text:
#=
Foil Shape equation - equal contour
1-1-2021
Author: Brigham Ostergaard
=#

#=
  foilShape(x, m)

  foilShape is used to determine the thickness of the foil at a given position
    along its lateral axis. x is the position along the foil and m is the value
    of maximum thickness

  t = 10m * (0.2969√(x) - 0.1260x -0.3537x^2 + 0.2843x^3 - 0.1015x^4)
=#

function foilShape(x, m)
    #Determines foil thickness at length x from left side.

  a =  0.2969*sqrt(x);
  b = 0.1260*x;
  c = 0.3537*(x^2);
  d = 0.2843 * (x^3);
  e = 0.1015 * (x^4);

  t = (10*m) * (a - b - c + d - e)

  return t
end

function _init_()
  m = 0.1;                  #coefficient specific to wing shape.

  print("x   | y")          #title for table
  println()                 #formating

  for x = 0:0.1:1
    y = foilShape(x, m);
    print(x, " | ", y)
    #println()              #formating
  end

  #println()                #formating
end

max_stars_repo_name: brigoste/VortexRings

max_issues_repo_name: brigoste/VortexRings

max_forks_repo_name: brigoste/VortexRings


text:
## solver

immutable ProximalDescent{S<:DescentSolver} <: Solver
    dsolver::DescentSolver
end

ProximalDescent{S<:DescentSolver}(solver::S) = ProximalDescent{S}(solver)

typealias ProxGD ProximalDescent{GD}
typealias ProxAGD ProximalDescent{AGD}

ProxGD() = ProximalDescent(GD())::ProxGD
ProxAGD() = ProximalDescent(AGD())::ProxAGD

## higher level

function _solve{T<:AbstractFloat}(
    pb::Problem{T},
    reg::Regularizer,
    θ::Array{T},
    solver::ProximalDescent,
    options::Options,
    callback::Function)

    if has_bias(pb)
        rmodel = riskmodel(pred_with_bias(pb), loss(pb))
        f = RiskFun(rmodel, inputs(pb), outputs(pb))
        solve!(solver, f, reg, θ, options, callback)

    else
        rmodel = riskmodel(pred_without_bias(pb), loss(pb))
        f = RiskFun(rmodel, inputs(pb), outputs(pb))
        solve!(solver, f, reg, θ, options, callback)

    end::Solution{typeof(θ)}
end

## core skeleton

function solve!{T<:AbstractFloat}(solver::ProximalDescent,
    f::Functional{T},       # the objective function (smooth part)
    reg::Regularizer,       # the regularizer (non-smooth part)
    θ::Array{T},            # the solution (which would be updated inplace)
    options::Options,       # options to control the procedure
    callback::Function)     # callback function

    ## extract arguments and options

    dsolver = solver.dsolver
    maxiter = options.maxiter
    vbose = verbosity_level(options.verbosity)::Int
    β = options.beta

    ## prepare storage

    θ2 = similar(θ)    # tempoarily new parameter (linear search)
    g = similar(θ)     # gradient w.r.t. f
    g2 = similar(θ)    # another gradient w.r.t. f
    p = prep_dir(dsolver, g)

    ## main loop
    t = 0
    converged = false
    vf, _ = value_and_grad!(f, g, θ)
    vr = value(reg, θ)
    v = vf + vr

    if vbose >= VERBOSE_ITER
        print_iter_head(with_gnorm=false)
        print_iter(t, v)
    end

    states = init_states(dsolver, θ, g)
    α = one(T)
    α_cnt = 0

    while !converged && t < maxiter
        t += 1
        v_pre = v

        # accelerate solution
        if has_accelerate(dsolver)
            accelerate!(dsolver, t, states, θ)
            vf, _ = value_and_grad!(f, g, θ)
            vr = value(reg, θ)
            v = vf + vr
        end

        # descent direction: g
        descent_dir!(dsolver, t, states, p, θ, θ2, g, g2)

        # backtracking
        α_pre = α
        if α_cnt > 3  # if α stays there for a while, we may begin with a bigger step
            α_pre = max(α_pre / β, one(T))
            α_cnt = 0
        end

        # Note: proximal operator has been applied during backtrack
        α = prox_backtrack!(f, reg, θ2, θ, vf, vr, g, p, α_pre, options)
        if α == α_pre
            α_cnt += 1
        end
        θ, θ2 = θ2, θ  # swap current solution and previous solution

        # compute new gradient
        vf, _ = value_and_grad!(f, g2, θ)
        vr = value(reg, θ)
        v = vf + vr
        g, g2 = g2, g  # swap current gradient with previous gradient

        # test convergence
        converged = test_convergence(θ, θ2, v, v_pre, options)

        # print iteration
        if vbose >= VERBOSE_ITER
            print_iter(t, v, v_pre, α)
        end

        # invoke callback (when requested)
        callback(t, θ, v, g)
    end

    if vbose >= VERBOSE_FINAL
        print_final(t, v, converged)
    end

    return Solution(θ, v, t, converged)
end


text:
using LinearAlgebra

function ESPRIT(Y,K,m=0)
    # Estimation of signal parameters via rotational invariance techniques
    # (ESPRIT) algorithm for line spectral estimation from measurements.
    # Can handle multiple snapshots by constructing a low rank Hankel matrix
    # from the data.
    #
    # Inputs:
    # Y: Data matrix
    # K: number of sinusoids
    # m: Hankel parameter of the Hankel matrix of the data
    #
    # Parameters:
    # N = size(Y,1): number of measurements per snapshot
    # L = size(Y,2): number of snapshots
    #
    # Output:
    # f_hat: estimate of the frequencies of the K sinusoids in [rad]

    (N,L) = size(Y)

    #C = X*X'/L

    if m == 0
    m = ceil(Int,N/2)
    end

    C = zeros(m,m)
    for j = 1:L, i = m:N
        C = C + Y[i-m+1:i,j]*Y[i-m+1:i,j]'/(L*N)
    end
    N = m

    evdC = eigen(C)
    p = sortperm(evdC.values,rev=true)[1:K]

    S = evdC.vectors[:,p]

    phi = S[1:N-1,:]\S[2:N,:]

    f_hat = ((angle.(eigvals(phi)))/(2*pi) .+ 1).%1
    #f_hat = real(log.(eigvals(phi))/(1im*2*pi))
end


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

using AWS.Compat
using AWS.UUIDs
"""
    AssociateNode()

 Associates a new node with the server. For more information about how to disassociate a node, see DisassociateNode.  On a Chef server: This command is an alternative to knife bootstrap.  Example (Chef): aws opsworks-cm associate-node --server-name MyServer --node-name MyManagedNode --engine-attributes \"Name=CHEF_ORGANIZATION,Value=default\" \"Name=CHEF_NODE_PUBLIC_KEY,Value=public-key-pem\"   On a Puppet server, this command is an alternative to the puppet cert sign command that signs a Puppet node CSR.   Example (Puppet): aws opsworks-cm associate-node --server-name MyServer --node-name MyManagedNode --engine-attributes \"Name=PUPPET_NODE_CSR,Value=csr-pem\"   A node can can only be associated with servers that are in a HEALTHY state. Otherwise, an InvalidStateException is thrown. A ResourceNotFoundException is thrown when the server does not exist. A ValidationException is raised when parameters of the request are not valid. The AssociateNode API call can be integrated into Auto Scaling configurations, AWS Cloudformation templates, or the user data of a server's instance.

# Required Parameters
- `EngineAttributes`: Engine attributes used for associating the node.   Attributes accepted in a AssociateNode request for Chef     CHEF_ORGANIZATION: The Chef organization with which the node is associated. By default only one organization named default can exist.     CHEF_NODE_PUBLIC_KEY: A PEM-formatted public key. This key is required for the chef-client agent to access the Chef API.     Attributes accepted in a AssociateNode request for Puppet     PUPPET_NODE_CSR: A PEM-formatted certificate-signing request (CSR) that is created by the node.
- `NodeName`: The name of the node.
- `ServerName`: The name of the server with which to associate the node.

"""

associate_node(EngineAttributes, NodeName, ServerName; aws_config::AWSConfig=global_aws_config()) = opsworkscm("AssociateNode", Dict{String, Any}("EngineAttributes"=>EngineAttributes, "NodeName"=>NodeName, "ServerName"=>ServerName); aws_config=aws_config)
associate_node(EngineAttributes, NodeName, ServerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("AssociateNode", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("EngineAttributes"=>EngineAttributes, "NodeName"=>NodeName, "ServerName"=>ServerName), args)); aws_config=aws_config)

"""
    CreateBackup()

 Creates an application-level backup of a server. While the server is in the BACKING_UP state, the server cannot be changed, and no additional backup can be created.   Backups can be created for servers in RUNNING, HEALTHY, and UNHEALTHY states. By default, you can create a maximum of 50 manual backups.   This operation is asynchronous.   A LimitExceededException is thrown when the maximum number of manual backups is reached. An InvalidStateException is thrown when the server is not in any of the following states: RUNNING, HEALTHY, or UNHEALTHY. A ResourceNotFoundException is thrown when the server is not found. A ValidationException is thrown when parameters of the request are not valid.

# Required Parameters
- `ServerName`: The name of the server that you want to back up.

# Optional Parameters
- `Description`:  A user-defined description of the backup.
- `Tags`: A map that contains tag keys and tag values to attach to an AWS OpsWorks-CM server backup.   The key cannot be empty.   The key can be a maximum of 127 characters, and can contain only Unicode letters, numbers, or separators, or the following special characters: + - = . _ : /    The value can be a maximum 255 characters, and contain only Unicode letters, numbers, or separators, or the following special characters: + - = . _ : /    Leading and trailing white spaces are trimmed from both the key and value.   A maximum of 50 user-applied tags is allowed for tag-supported AWS OpsWorks-CM resources.
"""

create_backup(ServerName; aws_config::AWSConfig=global_aws_config()) = opsworkscm("CreateBackup", Dict{String, Any}("ServerName"=>ServerName); aws_config=aws_config)
create_backup(ServerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("CreateBackup", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ServerName"=>ServerName), args)); aws_config=aws_config)

"""
    CreateServer()

 Creates and immedately starts a new server. The server is ready to use when it is in the HEALTHY state. By default, you can create a maximum of 10 servers.   This operation is asynchronous.   A LimitExceededException is thrown when you have created the maximum number of servers (10). A ResourceAlreadyExistsException is thrown when a server with the same name already exists in the account. A ResourceNotFoundException is thrown when you specify a backup ID that is not valid or is for a backup that does not exist. A ValidationException is thrown when parameters of the request are not valid.   If you do not specify a security group by adding the SecurityGroupIds parameter, AWS OpsWorks creates a new security group.   Chef Automate: The default security group opens the Chef server to the world on TCP port 443. If a KeyName is present, AWS OpsWorks enables SSH access. SSH is also open to the world on TCP port 22.   Puppet Enterprise: The default security group opens TCP ports 22, 443, 4433, 8140, 8142, 8143, and 8170. If a KeyName is present, AWS OpsWorks enables SSH access. SSH is also open to the world on TCP port 22.  By default, your server is accessible from any IP address. We recommend that you update your security group rules to allow access from known IP addresses and address ranges only. To edit security group rules, open Security Groups in the navigation pane of the EC2 management console.  To specify your own domain for a server, and provide your own self-signed or CA-signed certificate and private key, specify values for CustomDomain, CustomCertificate, and CustomPrivateKey.

# Required Parameters
- `Engine`:  The configuration management engine to use. Valid values include ChefAutomate and Puppet.
- `InstanceProfileArn`:  The ARN of the instance profile that your Amazon EC2 instances use. Although the AWS OpsWorks console typically creates the instance profile for you, if you are using API commands instead, run the service-role-creation.yaml AWS CloudFormation template, located at https://s3.amazonaws.com/opsworks-cm-us-east-1-prod-default-assets/misc/opsworks-cm-roles.yaml. This template creates a CloudFormation stack that includes the instance profile you need.
- `InstanceType`:  The Amazon EC2 instance type to use. For example, m5.large.
- `ServerName`:  The name of the server. The server name must be unique within your AWS account, within each region. Server names must start with a letter; then letters, numbers, or hyphens (-) are allowed, up to a maximum of 40 characters.
- `ServiceRoleArn`:  The service role that the AWS OpsWorks CM service backend uses to work with your account. Although the AWS OpsWorks management console typically creates the service role for you, if you are using the AWS CLI or API commands, run the service-role-creation.yaml AWS CloudFormation template, located at https://s3.amazonaws.com/opsworks-cm-us-east-1-prod-default-assets/misc/opsworks-cm-roles.yaml. This template creates a CloudFormation stack that includes the service role and instance profile that you need.

# Optional Parameters
- `AssociatePublicIpAddress`:  Associate a public IP address with a server that you are launching. Valid values are true or false. The default value is true.
- `BackupId`:  If you specify this field, AWS OpsWorks CM creates the server by using the backup represented by BackupId.
- `BackupRetentionCount`:  The number of automated backups that you want to keep. Whenever a new backup is created, AWS OpsWorks CM deletes the oldest backups if this number is exceeded. The default value is 1.
- `CustomCertificate`: A PEM-formatted HTTPS certificate. The value can be be a single, self-signed certificate, or a certificate chain. If you specify a custom certificate, you must also specify values for CustomDomain and CustomPrivateKey. The following are requirements for the CustomCertificate value:   You can provide either a self-signed, custom certificate, or the full certificate chain.   The certificate must be a valid X509 certificate, or a certificate chain in PEM format.   The certificate must be valid at the time of upload. A certificate can't be used before its validity period begins (the certificate's NotBefore date), or after it expires (the certificate's NotAfter date).   The certificate’s common name or subject alternative names (SANs), if present, must match the value of CustomDomain.   The certificate must match the value of CustomPrivateKey.
- `CustomDomain`: An optional public endpoint of a server, such as https://aws.my-company.com. To access the server, create a CNAME DNS record in your preferred DNS service that points the custom domain to the endpoint that is generated when the server is created (the value of the CreateServer Endpoint attribute). You cannot access the server by using the generated Endpoint value if the server is using a custom domain. If you specify a custom domain, you must also specify values for CustomCertificate and CustomPrivateKey.
- `CustomPrivateKey`: A private key in PEM format for connecting to the server by using HTTPS. The private key must not be encrypted; it cannot be protected by a password or passphrase. If you specify a custom private key, you must also specify values for CustomDomain and CustomCertificate.
- `DisableAutomatedBackup`:  Enable or disable scheduled backups. Valid values are true or false. The default value is true.
- `EngineAttributes`: Optional engine attributes on a specified server.   Attributes accepted in a Chef createServer request:     CHEF_AUTOMATE_PIVOTAL_KEY: A base64-encoded RSA public key. The corresponding private key is required to access the Chef API. When no CHEF_AUTOMATE_PIVOTAL_KEY is set, a private key is generated and returned in the response.     CHEF_AUTOMATE_ADMIN_PASSWORD: The password for the administrative user in the Chef Automate web-based dashboard. The password length is a minimum of eight characters, and a maximum of 32. The password can contain letters, numbers, and special characters (!/@#%^&+=_). The password must contain at least one lower case letter, one upper case letter, one number, and one special character. When no CHEF_AUTOMATE_ADMIN_PASSWORD is set, one is generated and returned in the response.    Attributes accepted in a Puppet createServer request:     PUPPET_ADMIN_PASSWORD: To work with the Puppet Enterprise console, a password must use ASCII characters.    PUPPET_R10K_REMOTE: The r10k remote is the URL of your control repository (for example, ssh://git@your.git-repo.com:user/control-repo.git). Specifying an r10k remote opens TCP port 8170.    PUPPET_R10K_PRIVATE_KEY: If you are using a private Git repository, add PUPPET_R10K_PRIVATE_KEY to specify a PEM-encoded private SSH key.
- `EngineModel`:  The engine model of the server. Valid values in this release include Monolithic for Puppet and Single for Chef.
- `EngineVersion`:  The major release version of the engine that you want to use. For a Chef server, the valid value for EngineVersion is currently 2. For a Puppet server, the valid value is 2017.
- `KeyPair`:  The Amazon EC2 key pair to set for the instance. This parameter is optional; if desired, you may specify this parameter to connect to your instances by using SSH.
- `PreferredBackupWindow`:  The start time for a one-hour period during which AWS OpsWorks CM backs up application-level data on your server if automated backups are enabled. Valid values must be specified in one of the following formats:     HH:MM for daily backups    DDD:HH:MM for weekly backups    MM must be specified as 00. The specified time is in coordinated universal time (UTC). The default value is a random, daily start time.  Example: 08:00, which represents a daily start time of 08:00 UTC.  Example: Mon:08:00, which represents a start time of every Monday at 08:00 UTC. (8:00 a.m.)
- `PreferredMaintenanceWindow`:  The start time for a one-hour period each week during which AWS OpsWorks CM performs maintenance on the instance. Valid values must be specified in the following format: DDD:HH:MM. MM must be specified as 00. The specified time is in coordinated universal time (UTC). The default value is a random one-hour period on Tuesday, Wednesday, or Friday. See TimeWindowDefinition for more information.   Example: Mon:08:00, which represents a start time of every Monday at 08:00 UTC. (8:00 a.m.)
- `SecurityGroupIds`:  A list of security group IDs to attach to the Amazon EC2 instance. If you add this parameter, the specified security groups must be within the VPC that is specified by SubnetIds.   If you do not specify this parameter, AWS OpsWorks CM creates one new security group that uses TCP ports 22 and 443, open to 0.0.0.0/0 (everyone).
- `SubnetIds`:  The IDs of subnets in which to launch the server EC2 instance.   Amazon EC2-Classic customers: This field is required. All servers must run within a VPC. The VPC must have \"Auto Assign Public IP\" enabled.   EC2-VPC customers: This field is optional. If you do not specify subnet IDs, your EC2 instances are created in a default subnet that is selected by Amazon EC2. If you specify subnet IDs, the VPC must have \"Auto Assign Public IP\" enabled.  For more information about supported Amazon EC2 platforms, see Supported Platforms.
- `Tags`: A map that contains tag keys and tag values to attach to an AWS OpsWorks for Chef Automate or AWS OpsWorks for Puppet Enterprise server.   The key cannot be empty.   The key can be a maximum of 127 characters, and can contain only Unicode letters, numbers, or separators, or the following special characters: + - = . _ : / @    The value can be a maximum 255 characters, and contain only Unicode letters, numbers, or separators, or the following special characters: + - = . _ : / @    Leading and trailing white spaces are trimmed from both the key and value.   A maximum of 50 user-applied tags is allowed for any AWS OpsWorks-CM server.
"""

create_server(Engine, InstanceProfileArn, InstanceType, ServerName, ServiceRoleArn; aws_config::AWSConfig=global_aws_config()) = opsworkscm("CreateServer", Dict{String, Any}("Engine"=>Engine, "InstanceProfileArn"=>InstanceProfileArn, "InstanceType"=>InstanceType, "ServerName"=>ServerName, "ServiceRoleArn"=>ServiceRoleArn); aws_config=aws_config)
create_server(Engine, InstanceProfileArn, InstanceType, ServerName, ServiceRoleArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("CreateServer", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Engine"=>Engine, "InstanceProfileArn"=>InstanceProfileArn, "InstanceType"=>InstanceType, "ServerName"=>ServerName, "ServiceRoleArn"=>ServiceRoleArn), args)); aws_config=aws_config)

"""
    DeleteBackup()

 Deletes a backup. You can delete both manual and automated backups. This operation is asynchronous.   An InvalidStateException is thrown when a backup deletion is already in progress. A ResourceNotFoundException is thrown when the backup does not exist. A ValidationException is thrown when parameters of the request are not valid.

# Required Parameters
- `BackupId`: The ID of the backup to delete. Run the DescribeBackups command to get a list of backup IDs. Backup IDs are in the format ServerName-yyyyMMddHHmmssSSS.

"""

delete_backup(BackupId; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DeleteBackup", Dict{String, Any}("BackupId"=>BackupId); aws_config=aws_config)
delete_backup(BackupId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DeleteBackup", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("BackupId"=>BackupId), args)); aws_config=aws_config)

"""
    DeleteServer()

 Deletes the server and the underlying AWS CloudFormation stacks (including the server's EC2 instance). When you run this command, the server state is updated to DELETING. After the server is deleted, it is no longer returned by DescribeServer requests. If the AWS CloudFormation stack cannot be deleted, the server cannot be deleted.   This operation is asynchronous.   An InvalidStateException is thrown when a server deletion is already in progress. A ResourceNotFoundException is thrown when the server does not exist. A ValidationException is raised when parameters of the request are not valid.

# Required Parameters
- `ServerName`: The ID of the server to delete.

"""

delete_server(ServerName; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DeleteServer", Dict{String, Any}("ServerName"=>ServerName); aws_config=aws_config)
delete_server(ServerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DeleteServer", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ServerName"=>ServerName), args)); aws_config=aws_config)

"""
    DescribeAccountAttributes()

 Describes your OpsWorks-CM account attributes.   This operation is synchronous.

"""

describe_account_attributes(; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DescribeAccountAttributes"; aws_config=aws_config)
describe_account_attributes(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DescribeAccountAttributes", args; aws_config=aws_config)

"""
    DescribeBackups()

 Describes backups. The results are ordered by time, with newest backups first. If you do not specify a BackupId or ServerName, the command returns all backups.   This operation is synchronous.   A ResourceNotFoundException is thrown when the backup does not exist. A ValidationException is raised when parameters of the request are not valid.

# Optional Parameters
- `BackupId`: Describes a single backup.
- `MaxResults`: This is not currently implemented for DescribeBackups requests.
- `NextToken`: This is not currently implemented for DescribeBackups requests.
- `ServerName`: Returns backups for the server with the specified ServerName.
"""

describe_backups(; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DescribeBackups"; aws_config=aws_config)
describe_backups(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DescribeBackups", args; aws_config=aws_config)

"""
    DescribeEvents()

 Describes events for a specified server. Results are ordered by time, with newest events first.   This operation is synchronous.   A ResourceNotFoundException is thrown when the server does not exist. A ValidationException is raised when parameters of the request are not valid.

# Required Parameters
- `ServerName`: The name of the server for which you want to view events.

# Optional Parameters
- `MaxResults`: To receive a paginated response, use this parameter to specify the maximum number of results to be returned with a single call. If the number of available results exceeds this maximum, the response includes a NextToken value that you can assign to the NextToken request parameter to get the next set of results.
- `NextToken`: NextToken is a string that is returned in some command responses. It indicates that not all entries have been returned, and that you must run at least one more request to get remaining items. To get remaining results, call DescribeEvents again, and assign the token from the previous results as the value of the nextToken parameter. If there are no more results, the response object's nextToken parameter value is null. Setting a nextToken value that was not returned in your previous results causes an InvalidNextTokenException to occur.
"""

describe_events(ServerName; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DescribeEvents", Dict{String, Any}("ServerName"=>ServerName); aws_config=aws_config)
describe_events(ServerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DescribeEvents", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ServerName"=>ServerName), args)); aws_config=aws_config)

"""
    DescribeNodeAssociationStatus()

 Returns the current status of an existing association or disassociation request.   A ResourceNotFoundException is thrown when no recent association or disassociation request with the specified token is found, or when the server does not exist. A ValidationException is raised when parameters of the request are not valid.

# Required Parameters
- `NodeAssociationStatusToken`: The token returned in either the AssociateNodeResponse or the DisassociateNodeResponse.
- `ServerName`: The name of the server from which to disassociate the node.

"""

describe_node_association_status(NodeAssociationStatusToken, ServerName; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DescribeNodeAssociationStatus", Dict{String, Any}("NodeAssociationStatusToken"=>NodeAssociationStatusToken, "ServerName"=>ServerName); aws_config=aws_config)
describe_node_association_status(NodeAssociationStatusToken, ServerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DescribeNodeAssociationStatus", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("NodeAssociationStatusToken"=>NodeAssociationStatusToken, "ServerName"=>ServerName), args)); aws_config=aws_config)

"""
    DescribeServers()

 Lists all configuration management servers that are identified with your account. Only the stored results from Amazon DynamoDB are returned. AWS OpsWorks CM does not query other services.   This operation is synchronous.   A ResourceNotFoundException is thrown when the server does not exist. A ValidationException is raised when parameters of the request are not valid.

# Optional Parameters
- `MaxResults`: This is not currently implemented for DescribeServers requests.
- `NextToken`: This is not currently implemented for DescribeServers requests.
- `ServerName`: Describes the server with the specified ServerName.
"""

describe_servers(; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DescribeServers"; aws_config=aws_config)
describe_servers(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DescribeServers", args; aws_config=aws_config)

"""
    DisassociateNode()

 Disassociates a node from an AWS OpsWorks CM server, and removes the node from the server's managed nodes. After a node is disassociated, the node key pair is no longer valid for accessing the configuration manager's API. For more information about how to associate a node, see AssociateNode.  A node can can only be disassociated from a server that is in a HEALTHY state. Otherwise, an InvalidStateException is thrown. A ResourceNotFoundException is thrown when the server does not exist. A ValidationException is raised when parameters of the request are not valid.

# Required Parameters
- `NodeName`: The name of the client node.
- `ServerName`: The name of the server from which to disassociate the node.

# Optional Parameters
- `EngineAttributes`: Engine attributes that are used for disassociating the node. No attributes are required for Puppet.   Attributes required in a DisassociateNode request for Chef     CHEF_ORGANIZATION: The Chef organization with which the node was associated. By default only one organization named default can exist.
"""

disassociate_node(NodeName, ServerName; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DisassociateNode", Dict{String, Any}("NodeName"=>NodeName, "ServerName"=>ServerName); aws_config=aws_config)
disassociate_node(NodeName, ServerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("DisassociateNode", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("NodeName"=>NodeName, "ServerName"=>ServerName), args)); aws_config=aws_config)

"""
    ExportServerEngineAttribute()

 Exports a specified server engine attribute as a base64-encoded string. For example, you can export user data that you can use in EC2 to associate nodes with a server.   This operation is synchronous.   A ValidationException is raised when parameters of the request are not valid. A ResourceNotFoundException is thrown when the server does not exist. An InvalidStateException is thrown when the server is in any of the following states: CREATING, TERMINATED, FAILED or DELETING.

# Required Parameters
- `ExportAttributeName`: The name of the export attribute. Currently, the supported export attribute is Userdata. This exports a user data script that includes parameters and values provided in the InputAttributes list.
- `ServerName`: The name of the server from which you are exporting the attribute.

# Optional Parameters
- `InputAttributes`: The list of engine attributes. The list type is EngineAttribute. An EngineAttribute list item is a pair that includes an attribute name and its value. For the Userdata ExportAttributeName, the following are supported engine attribute names.    RunList In Chef, a list of roles or recipes that are run in the specified order. In Puppet, this parameter is ignored.    OrganizationName In Chef, an organization name. AWS OpsWorks for Chef Automate always creates the organization default. In Puppet, this parameter is ignored.    NodeEnvironment In Chef, a node environment (for example, development, staging, or one-box). In Puppet, this parameter is ignored.    NodeClientVersion In Chef, the version of the Chef engine (three numbers separated by dots, such as 13.8.5). If this attribute is empty, OpsWorks for Chef Automate uses the most current version. In Puppet, this parameter is ignored.
"""

export_server_engine_attribute(ExportAttributeName, ServerName; aws_config::AWSConfig=global_aws_config()) = opsworkscm("ExportServerEngineAttribute", Dict{String, Any}("ExportAttributeName"=>ExportAttributeName, "ServerName"=>ServerName); aws_config=aws_config)
export_server_engine_attribute(ExportAttributeName, ServerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("ExportServerEngineAttribute", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ExportAttributeName"=>ExportAttributeName, "ServerName"=>ServerName), args)); aws_config=aws_config)

"""
    ListTagsForResource()

Returns a list of tags that are applied to the specified AWS OpsWorks for Chef Automate or AWS OpsWorks for Puppet Enterprise servers or backups.

# Required Parameters
- `ResourceArn`: The Amazon Resource Number (ARN) of an AWS OpsWorks for Chef Automate or AWS OpsWorks for Puppet Enterprise server for which you want to show applied tags. For example, arn:aws:opsworks-cm:us-west-2:123456789012:server/test-owcm-server/EXAMPLE-66b0-4196-8274-d1a2bEXAMPLE.

# Optional Parameters
- `MaxResults`: To receive a paginated response, use this parameter to specify the maximum number of results to be returned with a single call. If the number of available results exceeds this maximum, the response includes a NextToken value that you can assign to the NextToken request parameter to get the next set of results.
- `NextToken`: NextToken is a string that is returned in some command responses. It indicates that not all entries have been returned, and that you must run at least one more request to get remaining items. To get remaining results, call ListTagsForResource again, and assign the token from the previous results as the value of the nextToken parameter. If there are no more results, the response object's nextToken parameter value is null. Setting a nextToken value that was not returned in your previous results causes an InvalidNextTokenException to occur.
"""

list_tags_for_resource(ResourceArn; aws_config::AWSConfig=global_aws_config()) = opsworkscm("ListTagsForResource", Dict{String, Any}("ResourceArn"=>ResourceArn); aws_config=aws_config)
list_tags_for_resource(ResourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("ListTagsForResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn), args)); aws_config=aws_config)

"""
    RestoreServer()

 Restores a backup to a server that is in a CONNECTION_LOST, HEALTHY, RUNNING, UNHEALTHY, or TERMINATED state. When you run RestoreServer, the server's EC2 instance is deleted, and a new EC2 instance is configured. RestoreServer maintains the existing server endpoint, so configuration management of the server's client devices (nodes) should continue to work.  Restoring from a backup is performed by creating a new EC2 instance. If restoration is successful, and the server is in a HEALTHY state, AWS OpsWorks CM switches traffic over to the new instance. After restoration is finished, the old EC2 instance is maintained in a Running or Stopped state, but is eventually terminated.  This operation is asynchronous.   An InvalidStateException is thrown when the server is not in a valid state. A ResourceNotFoundException is thrown when the server does not exist. A ValidationException is raised when parameters of the request are not valid.

# Required Parameters
- `BackupId`:  The ID of the backup that you want to use to restore a server.
- `ServerName`:  The name of the server that you want to restore.

# Optional Parameters
- `InstanceType`:  The type of instance to restore. Valid values must be specified in the following format: ^([cm][34]|t2).* For example, m5.large. Valid values are m5.large, r5.xlarge, and r5.2xlarge. If you do not specify this parameter, RestoreServer uses the instance type from the specified backup.
- `KeyPair`:  The name of the key pair to set on the new EC2 instance. This can be helpful if the administrator no longer has the SSH key.
"""

restore_server(BackupId, ServerName; aws_config::AWSConfig=global_aws_config()) = opsworkscm("RestoreServer", Dict{String, Any}("BackupId"=>BackupId, "ServerName"=>ServerName); aws_config=aws_config)
restore_server(BackupId, ServerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("RestoreServer", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("BackupId"=>BackupId, "ServerName"=>ServerName), args)); aws_config=aws_config)

"""
    StartMaintenance()

 Manually starts server maintenance. This command can be useful if an earlier maintenance attempt failed, and the underlying cause of maintenance failure has been resolved. The server is in an UNDER_MAINTENANCE state while maintenance is in progress.   Maintenance can only be started on servers in HEALTHY and UNHEALTHY states. Otherwise, an InvalidStateException is thrown. A ResourceNotFoundException is thrown when the server does not exist. A ValidationException is raised when parameters of the request are not valid.

# Required Parameters
- `ServerName`: The name of the server on which to run maintenance.

# Optional Parameters
- `EngineAttributes`: Engine attributes that are specific to the server on which you want to run maintenance.  Attributes accepted in a StartMaintenance request for Chef     CHEF_MAJOR_UPGRADE: If a Chef Automate server is eligible for upgrade to Chef Automate 2, add this engine attribute to a StartMaintenance request and set the value to true to upgrade the server to Chef Automate 2. For more information, see Upgrade an AWS OpsWorks for Chef Automate Server to Chef Automate 2.
"""

start_maintenance(ServerName; aws_config::AWSConfig=global_aws_config()) = opsworkscm("StartMaintenance", Dict{String, Any}("ServerName"=>ServerName); aws_config=aws_config)
start_maintenance(ServerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("StartMaintenance", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ServerName"=>ServerName), args)); aws_config=aws_config)

"""
    TagResource()

Applies tags to an AWS OpsWorks for Chef Automate or AWS OpsWorks for Puppet Enterprise server, or to server backups.

# Required Parameters
- `ResourceArn`: The Amazon Resource Number (ARN) of a resource to which you want to apply tags. For example, arn:aws:opsworks-cm:us-west-2:123456789012:server/test-owcm-server/EXAMPLE-66b0-4196-8274-d1a2bEXAMPLE.
- `Tags`: A map that contains tag keys and tag values to attach to AWS OpsWorks-CM servers or backups.   The key cannot be empty.   The key can be a maximum of 127 characters, and can contain only Unicode letters, numbers, or separators, or the following special characters: + - = . _ : /    The value can be a maximum 255 characters, and contain only Unicode letters, numbers, or separators, or the following special characters: + - = . _ : /    Leading and trailing white spaces are trimmed from both the key and value.   A maximum of 50 user-applied tags is allowed for any AWS OpsWorks-CM server or backup.

"""

tag_resource(ResourceArn, Tags; aws_config::AWSConfig=global_aws_config()) = opsworkscm("TagResource", Dict{String, Any}("ResourceArn"=>ResourceArn, "Tags"=>Tags); aws_config=aws_config)
tag_resource(ResourceArn, Tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("TagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn, "Tags"=>Tags), args)); aws_config=aws_config)

"""
    UntagResource()

Removes specified tags from an AWS OpsWorks-CM server or backup.

# Required Parameters
- `ResourceArn`: The Amazon Resource Number (ARN) of a resource from which you want to remove tags. For example, arn:aws:opsworks-cm:us-west-2:123456789012:server/test-owcm-server/EXAMPLE-66b0-4196-8274-d1a2bEXAMPLE.
- `TagKeys`: The keys of tags that you want to remove.

"""

untag_resource(ResourceArn, TagKeys; aws_config::AWSConfig=global_aws_config()) = opsworkscm("UntagResource", Dict{String, Any}("ResourceArn"=>ResourceArn, "TagKeys"=>TagKeys); aws_config=aws_config)
untag_resource(ResourceArn, TagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("UntagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceArn"=>ResourceArn, "TagKeys"=>TagKeys), args)); aws_config=aws_config)

"""
    UpdateServer()

 Updates settings for a server.   This operation is synchronous.

# Required Parameters
- `ServerName`: The name of the server to update.

# Optional Parameters
- `BackupRetentionCount`: Sets the number of automated backups that you want to keep.
- `DisableAutomatedBackup`: Setting DisableAutomatedBackup to true disables automated or scheduled backups. Automated backups are enabled by default.
- `PreferredBackupWindow`:
- `PreferredMaintenanceWindow`:
"""

update_server(ServerName; aws_config::AWSConfig=global_aws_config()) = opsworkscm("UpdateServer", Dict{String, Any}("ServerName"=>ServerName); aws_config=aws_config)
update_server(ServerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("UpdateServer", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ServerName"=>ServerName), args)); aws_config=aws_config)

"""
    UpdateServerEngineAttributes()

 Updates engine-specific attributes on a specified server. The server enters the MODIFYING state when this operation is in progress. Only one update can occur at a time. You can use this command to reset a Chef server's public key (CHEF_PIVOTAL_KEY) or a Puppet server's admin password (PUPPET_ADMIN_PASSWORD).   This operation is asynchronous.   This operation can only be called for servers in HEALTHY or UNHEALTHY states. Otherwise, an InvalidStateException is raised. A ResourceNotFoundException is thrown when the server does not exist. A ValidationException is raised when parameters of the request are not valid.

# Required Parameters
- `AttributeName`: The name of the engine attribute to update.
- `ServerName`: The name of the server to update.

# Optional Parameters
- `AttributeValue`: The value to set for the attribute.
"""

update_server_engine_attributes(AttributeName, ServerName; aws_config::AWSConfig=global_aws_config()) = opsworkscm("UpdateServerEngineAttributes", Dict{String, Any}("AttributeName"=>AttributeName, "ServerName"=>ServerName); aws_config=aws_config)
update_server_engine_attributes(AttributeName, ServerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = opsworkscm("UpdateServerEngineAttributes", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AttributeName"=>AttributeName, "ServerName"=>ServerName), args)); aws_config=aws_config)


text:
# Autogenerated wrapper script for GAP_lib_jll for any
JLLWrappers.@generate_wrapper_header("GAP_lib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using LowRankModels
using Base.Test

# Parameters
n = 200
m = 200
r = 5
eta = 0.01
delta = 1e-3

# Generate problem
srand(1)
Um = randn(r, n)
Vm = randn(r, m)
U = Um .+ sqrt(eta) * randn(r, n)
V = Vm .+ sqrt(eta) * randn(r, m)
Y = U' * V + sqrt(delta) * randn(n, m)

# Run algorithm
glrm = GLRM(Y, QuadLoss(), [RemQuadReg(50, Um[:, i]) for i = 1:n],
    [RemQuadReg(50, Vm[:, j]) for j = 1:m], r)
Uh, Vh, iter_info = fit!(glrm)

mseU, mseV, mseY = mean((U - Uh).^2), mean((V - Vh).^2), mean((Y - Uh' * Vh).^2)
@printf("MSE(U) = %.4g, MSE(V) = %.4g, MSE(Y) = %.4g\n", mseU, mseV, mseY)

# Perform some tests
@test mseU < 1e-3
@test mseV < 1e-3


text:
using Random
using LinearAlgebra

export spqr_ssp, normest, normestinv

#==
SPQR_SSP block power method or subspace iteration applied to A or A*N

 U,S,V,stats = spqr_ssp(A, N, k, opts)

 Uses the block power method or subspace iteration to estimate the largest k
 singular values, and left and right singular vectors an m-by-n matrix A or,
 if the n-by-p matrix N is included in the input, of A*N.

 Alternate input syntaxes (opts are set to defaults if not present):

       spqr_ssp(A)         k = 1, N is [ ]
       spqr_ssp(A,k)       k is a scalar
       spqr_ssp(A,N)       N is a matrix or struct
       spqr_ssp(A,N,k)     k is a scalar
       spqr_ssp(A,N,opts)  opts is a struct

 Input:
   A -- an m-by-n matrix
   N (optional) -- an n-by-p matrix representing the null space of A, stored
       either explicitly as a matrix or implicitly as a structure, and
       returned by spqr_basic, spqr_null, spqr_cod, or spqr_pinv.
       N is ignored if [ ].
   k (optional) -- the number of singular values to calculate. Also the block
       size in the block power method.  Also can be specified as opts.k.
       Default: 1.
   opts (optional) -- type 'help spqr_rank_opts' for details

 Output:
   U  -- m-by-k matrix containing estimates of the left singular vectors of A
       (or A*N if N is included in the input) corresponding to the singular
       values in S.
   S -- k-by-k diagonal matrix whose diagonal entries are the estimated
       singular values of A (or A*N). Also for i = 1:nsval, S[i,i] is a lower
       bound on singular value i of A (or A*N).
   V -- n-by-k matrix containing estimates of the right singular vectors of A
       corresponding to the singular values in S.
   stats -- type 'help spqr_rank_stats'

 Note that A*V == U*S (or A*N*V == U*S) and that U'*A (or U'*A*N)
 is approximately equal to S*V'.

 Output (for one or two output arguments): s,stats= spqr_ssp(...)
   s -- diag(S), with S as above

 Example:
    A = MatrixDepot.matrixdepot("kahan", 100))
    U,S,V = spqr_ssp(A)   # singular triple for largest singular value
    # or
    U,S,V,stats = spqr_ssp(A,4)  # same for largest 4 singular values
    # or
    s,stats = spqr_ssp(A,4)
    stats     # stats has information for largest 4 singular values
    # or
    N = spqr_null(A)   # N has a basis for the null space of A
    s = spqr_ssp(A,N)
    s     # s has estimate of largest singular value (or norm) of A*N

 See also spqr_basic, spqr_null, spqr_pinv, spqr_cod.

 Copyright 2012, Leslie Foster and Timothy A Davis.

 Outline of algorithm:
    Let B = A or B = A*N, if the optional input N is included
    where A is m-by-n and N is n-by-p. Note that the code is written so
    that B is not explicitly formed.  Let n = p if N is not input.
    Let k = block size = min(k,m,p)
    U = m-by-k random matrix with orthonormal columns
    repeat
       V1 = B' * U
       V, D1, X1= svd(V1) = the compact svd of V1
       U1 = B * V
       U, D2, X2 = svd(U1) = the compact svd of U1

       kth_est_error_bound = norm(B'*U[:,k] - V*(X2[:,k]*D2[k,k])) / sqrt(2)
       if kth_est_error_bound <= convergence_factor*s[k] && iters >= min_iters
            flag = 0
            exit loop
       elseif number of iteration > max number of iterations
            flag = 1
            exit loop
       else
            continue loop
       end
    end
    V = V*X2 # adjust V so that A*V = U*S
    calculate estimated error bounds for singular values 1:k-1
       as discussed in the code
    S = diag(s)
==#

"""
start block power method to estimate the largest singular values and the
    corresponding right singular vectors (in the n-by-k matrix V) and left
    singular vectors (in the m-by-k matrix U) of the m-by-n matrix A
"""
function spqr_ssp(A::AbstractMatrix{T}, N::Union{UniformScaling{T}, AbstractArray{T},NamedTuple}=one(T)I;
                  nargout=4, opts...) where T<:Number

##           k::Int = 1,   # the # of singular values to compute.
##           min_iters::Int = 4, # min # of iterations before checking convergence.
##           max_iters::Int = 10, # max # of iterations before stopping the iterations.
#       The default value = 10 appears, with the default value of opts.tol,
#       to provide sufficient accuracy to correctly determine the numerical
#       rank when spqr_ssi is called by spqr_basic, spqr_null, spqr_pinv or
#       spqr_cod in almost all cases, assuming that stats.flag is 0. For
#       values of opts.tol larger than the default value, a larger value of
#       opts.ssp_max_iters, for example 100, may be useful.
##            convergence_factor::Real = 0.1, # continue power method iterations until an
#       estimated bound on the relative error in the approximation
#       S[k] to singular value number k of A is <= convergence_factor.
#       The default value = 0.1 appears, with the default value of
#       opts.tol, to provide sufficient accuracy to correctly determine
#       the numerical rank in spqr_basic, spqr_null, spqr_pinv or spqr_cod
#       in almost all cases, assuming that stats.flag is 0.  For values
#       of opts.tol larger than the default value, a smaller value of
#       opts.ssp_convergence_factor, for example 0.01, may be useful.
##            get_details::Int = 1, # determine what statistics to return.
#       0: basic statistics
#       1: extensive statistics.
#       2: basic statistics and a few addional statistics.
#       See 'help spqr_rank_stats' for details.
##            repeatable::Bool = true # controls the random stream.
#       false: use the current random number stream.  The state of the stream will
#           be different when the method returns.  Repeated calls may generate
#           different results.  Faster for small matrices, however.
#       true: use a repeatable internal random number stream.  The current stream
#           will not be affected.  Repeated calls generate the same results.

    start_tic = time_ns()
    opts, get_details, repeatable, k, min_iters, max_iters, convergence_factor =
    get_opts(opts, :get_details, :repeatable, :k, :ssp_min_iters, :ssp_max_iters, :ssp_convergence_factor)
    stats = Statistics(real(T))
    stats.time_initialize = time_ns() - start_tic

    private_stream = repeatable ? MersenneTwister(1) : Random.GLOBAL_RNG
    #-------------------------------------------------------------------------------
    # initializations
    #-------------------------------------------------------------------------------

    m, n = size(A)
    n1 = N isa UniformScaling ? n : size(N, 1)
    n2 = N isa UniformScaling ? n : size(N, 2)
    k = max(k, 0)
    k = min(k, m, n2)     # number of singular values to compute
    max_iters = max(1, max_iters)

    stats.flag = 1
    stats.est_error_bounds = zeros(T, max(k, 1))
    #stats.sval_numbers_for_bounds = 1:k
    stats.opts_used = opts.data

    if k <= 0
        # quick return.  This is not an error condition.
        stats.flag = 0
        stats.est_svals = T[]
        stats.sval_numbers_for_bounds = 1:0
        U = zeros(m, 0)
        S = zeros(0)
        V = zeros(n, 0)
        return nargout == 1 ? S : nargout <= 2 ? (S, stats) : (U, S, V, stats)
    end

    U = randn(private_stream, m, k)
    U = qr(U).Q * Matrix{T}(I, m, k)  # iteration start - U has k orthonormal columns

    #-------------------------------------------------------------------------------
    # block power iterations
    #-------------------------------------------------------------------------------
    get_details >= 1 && ( t = time_ns() )

    D2 = T[]
    X2 = T[]
    kth_est_error_bound = T(0)

    ASU = similar(U, n, k)
    V1 = similar(U, n2, k)
    NV = similar(U, n1, k)
    U1 = similar(U)
    ASU1 = similar(U, n)
    V11 = similar(U, size(V1, 1))

    crit = eps(norm(A)) * min(m, n)

    iters = 0
    while iters < max_iters
        iters += 1
        mul!(ASU, A', U)
        mul!(V1, N', ASU)

        get_details == 1 && (time_svd = time_ns())
        V = svd(V1).U
        get_details == 1 && (stats.time_svd += time_ns() - time_svd)

        mul!(NV, N, V)
        mul!(U1, A, NV)

        get_details == 1 && (time_svd = time_ns())
        U, D2, X2 = svd(U1)
        get_details == 1 && (stats.time_svd += time_ns() - time_svd)

        # estimate error bound on kth singular value
        mul!(ASU1, A', U[:,k])
        mul!(V11, N', ASU1)
        V11 .-= V * ( X2[:,k] * D2[k] )
        kth_est_error_bound = norm(V11) / sqrt(2)

        # Note that
        #     [ 0     B ] [  U   ]  =    [       U * D2          ]
        #     [ B'    0 ] [ V*X2 ]       [     V * ( X2 * D2 )   ]
        # where B = A or B = A*N (if N is included in the input).  It follows that,
        # by Demmel's Applied Numerical Linear Algebra, Theorem 5.5, some
        # eigenvalue of
        #     C  =  [ 0    B ]
        #           [ B'   0 ]
        # will be within a distance of kth_est_error_bound of D2[k,k].  Typically
        # the eigenvalue of C is the kth singular value of B, although this is not
        # guaranteed.  kth_est_error_bound is our estimate of a bound on the error
        # in using D2[k,k] to approximate kth singular value of B.

        crit = max(convergence_factor * D2[k], crit)

        println("ssp: kth_est[$iters] = $kth_est_error_bound crit = $crit D2[$k] = $(D2[k])" )

        if kth_est_error_bound <= crit && iters >= min_iters
            # The test indicates that the estimated relative error in
            # the estimate, D2[k,k] for the kth singular value is smaller than
            # convergence_factor, which is the goal.
            stats.flag = 0
            break
        end
    end

    if get_details == 1
        stats.iters = iters
        stats.time_iters = time_ns() - t
    end

    S = D2[1:k]
    stats.est_error_bounds[k] = kth_est_error_bound
    stats.est_svals = S

    # adjust V so that A*V = U*S
    V = V*X2

    #-------------------------------------------------------------------------------
    # estimate error bounds for the 1st through (k-1)st singular values
    #-------------------------------------------------------------------------------

    if get_details >= 1
        t = time_ns()
    end

    if k > 1

        V1 = N' * (A' * U[:,1:k-1])
        U0 = V1 - V[:,1:k-1] * Diagonal(S[1:k-1])
        stats.est_error_bounds[1:k-1] .= vec(sqrt.(sum(abs2, U0, dims=1))) / sqrt(2)

        # Note (with V adjusted as above) that
        #    [ 0      B ]  [ U ]   -    [ U ] * S = [       0       ]
        #    [ B'     0 ]  [ V ]   -    [ V ]       [   B'*U - V*S  ]
        # where B = A or B = A*N (if N is included in the input).  It follows, by
        # Demmel's Applied Numerical Linear Algebra, Theorem 5.5, for i = 1:k, some
        # eigenvalue of
        #     C  =  [ 0    B ]
        #           [ B'   0 ]
        # will be within a distance of the norm of the ith column of [ B' * U - V *
        # S ] / sqrt(2) from S[i,i]. Typically this eigenvalue will be the ith
        # singular value number of B, although it is not guaranteed that the
        # singular value number is correct.  stats.est_error_bounds[i] is the norm
        # of the ith column of [ B' * U - V * S ] / sqrt(2).

    end

    #-------------------------------------------------------------------------------
    # return results
    #-------------------------------------------------------------------------------

    if get_details >= 1
        stats.time_est_error_bounds = time_ns() - t
        stats.time = time_ns() - start_tic
    end

    nargout  == 1 ? S : nargout <= 2 ? (S, stats) : (U, S, V, stats)
end

"""
    normest(A; [tol=convergence_factor], [max_iters=...])

Estimate 2-norm of matrix A, power-iterating `A'A` up to given tolerance.
"""
function normest(A::AbstractMatrix{T}; tol::AbstractFloat = 1e-6,
                 max_iters::Integer = 100) where T <:Number

    min_iters = 2
    max_iters = max(max_iters, min_iters)
    tol = max(tol, eps())
    ncols = size(A, 2)
    y = randn(T, ncols)
    k = argmax(sum(abs, A, dims=1))[2]
    y[k] += sign(y[k]) * 0.1
    lmul!(1/norm(y), y)
    i = 0
    n = 0.0
    n0 = n
    x = Vector{T}(undef, size(A,1))
    while ( abs(n - n0) > tol * n || i < min_iters ) && i < max_iters
        i += 1
        n0 = n
        mul!(x, A, y)
        normx = norm(x)
        if normx == 0
            x = rand(T, ncols)
        end
        lmul!(1/norm(x), x)
        mul!(y, A', x)
        n = norm(y)
    end
    n
end

"""
    normestinv(A; [tol=convergence_factor], [max_iters=...])

Estimate 2-norm of inverse of matrix A, power-iterating `A⁻¹'A⁻¹` up to given tolerance.
"""
function normestinv(A::AbstractMatrix{T}; tol::AbstractFloat = 1e-6,
                    max_iters::Integer = 100) where T <:Number

    if istriu(A)
        _normestinv(UpperTriangular(A), tol, max_iters)
    elseif istril(A)
        _normestinv(LowerTriangular(A), tol, max_iters)
    else
        _normestinv(A, tol, max_iters)
    end
end

function _normestinv(A::Union{Adjoint{<:Any,<:AbstractSparseMatrix},
                              Transpose{<:Any,<:AbstractSparseMatrix}},
                              tol::AbstractFloat, max_iters::Integer)

    _normestinv(A.parent, tol, max_iters)
end

function _normestinv(A::AbstractSparseMatrix{T}, tol::AbstractFloat = 1e-6,
                    max_iters::Integer = 100) where T <:Number

    m, n = size(A)
    tol2 = max(tol^2, eps())
    if m < n
        A = sparse(A')
    end
    qrf = qr(A, tol=tol^2) # TODO that is not clean

    R = qrf.R
    m1, n1 = size(R)
    if min(m1, n1) < min(m, n)
        return Inf
    elseif m1 > n1
        R = SparseMatrixCSC(n1, n1, R.colptr, R.rowval, R.nzval)
    elseif m1 < n1
        R = SparseMatrixCSC(m1, m1, R.colptr[1:m1+1], R.rowval, R.nzval)
    end
    _normestinv(UpperTriangular(R), tol, max_iters)
end

function _normestinv(A::AbstractMatrix{T}, tol::AbstractFloat = 1e-6,
                    max_iters::Integer = 100) where T <:Number

    m, n = size(A)
    if m < n
        A = A'
    end
    qrf = qr(A, Val(true))
    R = qrf.R
    m1, n1 = size(R)
    if min(m1, n1) < min(m, n)
        Inf
    else
        R = view(R, 1:m1, 1:m1)
    end
    _normestinv(UpperTriangular(R), tol, max_iters)
end

function normestinv(A::LinearAlgebra.AbstractTriangular{T}; tol::AbstractFloat = 1e-6,
                    max_iters::Integer = 100) where T <:Number
    _normestinv(A, tol, max_iters)
end

function _normestinv(A::LinearAlgebra.AbstractTriangular{T}, tol::AbstractFloat,
                    max_iters::Integer) where T <:Number

    min_iters = 2
    max_iters = max(max_iters, min_iters)
    tol = max(tol, eps())
    ncols = size(A, 2)
    y = randn(T, ncols)
    lmul!(1/norm(y), y)
    i = 0
    n = 0.0
    n0 = n
    x = Vector{T}(undef, size(A,1))
    while ( abs(n - n0) > tol * n || i < min_iters ) && i < max_iters
        i += 1
        n0 = n
        # ldiv!(x, A, y)
        x = A \ y
        normx = norm(x)
        if normx == 0
            x = rand(T, ncols)
        end
        lmul!(1/norm(x), x)
        #ldiv!(y, A', x)
        y = A' \ x
        n = norm(y)
    end
    n
end