text:
data_array = [0.405503  0.707334  0.211959  0.2539  0.479168  0.161135  -0.0888367
0.415503  0.707334  0.211959  0.2539  0.479168  0.153381  -0.102154
0.405503  0.717334  0.211959  0.2539  0.479168  0.14144   -0.106894
0.405503  0.707334  0.221959  0.2539  0.479168  0.150405  -0.0882103
0.405503  0.707334  0.211959  0.2639  0.479168  0.135277  -0.108619
0.405503  0.707334  0.211959  0.2539  0.489168  0.162739  -0.0966177]

real_gd, imag_gd = get_val_gradient_from_data(data_array, δ=1e-2)

@test imag_gd[1, end-4] == (data_array[2, end] - data_array[1, end])*1e2
@test imag_gd[1, end-2] == (data_array[4, end] - data_array[1, end])*1e2
@test imag_gd[1, end] == (data_array[6, end] - data_array[1, end])*1e2
@test real_gd[1, end-4] == (data_array[2, end-1] - data_array[1, end-1])*1e2

# note that this value is wrong because the data was simulated with a δ=1e-2 pertubation
real_gd, imag_gd = get_val_gradient_from_data(data_array, δ=1e-4)

@test imag_gd[1, end-4] ≈ (data_array[2, end] - data_array[1, end])*1e4


text:
using BlackBoxOptim

# From Hans W. Borchers in
# https://groups.google.com/forum/#!topic/julia-opt/HltM-OGjueo
function trefethen(p)
  x = p[1]; y = p[2]
  return exp(sin(50 * x)) + sin(60 * exp(y)) + sin(70 * sin(x)) +
             sin(sin(80 * y)) - sin(10 * (x + y)) + (x^2 + y^2)/4
end

bboptimize(trefethen; SearchRange = (-1.0, 1.0), NumDimensions = 2,
    MaxTime = 2.0,
    Method = :generating_set_search)

# To compare several methods in several repetitions:
#BlackBoxOptim.repeated_bboptimize(5, p, 2, [
#  :generating_set_search,
#  :adaptive_de_rand_1_bin_radiuslimited,
#  :random_search,
#  :separable_nes],
#  5.0)

ext: jl


text:
@inline function Phi_s(Cell,s,z,Def,ϕ_tf,D,res0)
    """
    Solid Potential Transfer Function

    Phi_s(Cell,s,z,Def)

    """

 if Def == "Pos"
    Electrode = Cell.Pos #Electrode Length
 else
    Electrode = Cell.Neg #Electrode Length
 end

κ_eff = Cell.Const.κ*Electrode.ϵ_e^Electrode.κ_brug #Effective Electrolyte Conductivity
σ_eff = Electrode.σ*Electrode.ϵ_s^Electrode.σ_brug #Effective Electrode Conductivity
comb_cond_eff = κ_eff+σ_eff #combining into single variable

#Defining SOC
θ = Cell.Const.SOC * (Electrode.θ_100-Electrode.θ_0) + Electrode.θ_0

#Prepare for j0
cs0 = Electrode.cs_max * θ

#Current Flux Density
if Cell.Const.CellTyp == "Doyle_94"
   κ = Electrode.k_norm/Electrode.cs_max/Cell.Const.ce0^(1-Electrode.α)
   j0 = κ*(Cell.Const.ce0*(Electrode.cs_max-cs0))^(1-Electrode.α)*cs0^Electrode.α
else
   j0 = Electrode.k_norm*(Cell.Const.ce0*cs0*(Electrode.cs_max-cs0))^(1-Electrode.α)
end

#Resistance
Rtot = R*Cell.Const.T/(j0*F^2) + Electrode.RFilm
#Rtot = R*Cell.Const.T/(j0*Cell.Const.CC_A*F) + Electrode.RFilm

#∂Uocp_Def
∂Uocp_elc = Cell.Const.∂Uocp(Def,θ)/Electrode.cs_max

ν = @. Electrode.L*sqrt((Electrode.as/σ_eff+Electrode.as/κ_eff)/(Rtot+∂Uocp_elc*(Electrode.Rs/(F*Electrode.Ds))*(tanh(Electrode.β)/(tanh(Electrode.β)-Electrode.β)))) #Condensing Variable - eq. 4.13
ν_∞ = @. Electrode.L*sqrt((Electrode.as/κ_eff+Electrode.as/σ_eff)/(Rtot))

ϕ_tf .= @. (-Electrode.L*(κ_eff*(cosh(ν)-cosh((z-1)*ν)))-Electrode.L*(σ_eff*(1-cosh(z*ν)+z*ν*sinh(ν))))/(Cell.Const.CC_A*σ_eff*(comb_cond_eff)*ν*sinh(ν)) #Transfer Function - eq. 4.19
D .= @. (-Electrode.L*(κ_eff*(cosh(ν_∞)-cosh((z-1)*ν_∞)))-Electrode.L*(σ_eff*(1-cosh(z*ν_∞)+z*ν_∞*sinh(ν_∞))))/(Cell.Const.CC_A*σ_eff*(comb_cond_eff)*ν_∞*sinh(ν_∞)) # Contribution to D as G->∞
zero_tf = @. Electrode.L*(z-2)*z/(2*Cell.Const.CC_A*σ_eff)
ϕ_tf[:,findall(s.==0)] .= zero_tf[:,findall(s.==0)]
res0 .= zeros(length(z))

if Def == "Pos"
   ϕ_tf .= -ϕ_tf
   D .= -D
end

end

max_stars_repo_name: BradyPlanden/LiBRA

max_forks_repo_name: BradyPlanden/LiBRA


text:
#
# List of elements with properties
#
struct Element
    atomic_number::Int
    element::String
    pdb_name::String
    name::String
    mass::Float64
end
const elements = [
    Element(0, "X", "X", "NotFound", 0.00000)
    Element(1, "H", "H", "Hydrogen", 1.00797)
    Element(2, "He", "HEL", "Helium", 4.00260)
    Element(3, "Li", "LI", "Lithium", 6.941)
    Element(4, "Be", "BE", "Beryllium", 9.01218)
    Element(5, "B", "B", "Boron", 10.81)
    Element(6, "C", "C", "Carbon", 12.011)
    Element(7, "N", "N", "Nitrogen", 14.0067)
    Element(8, "O", "O", "Oxygen", 15.9994)
    Element(9, "F", "F", "Fluorine", 18.998403)
    Element(10, "Ne", "NEO", "Neon", 20.179)
    Element(11, "Na", "NA", "Sodium", 22.98977)
    Element(11, "Na", "SOD", "Sodium", 22.98977)
    Element(12, "Mg", "MG", "Magnesium", 24.305)
    Element(13, "Al", "AL", "Aluminum", 26.98154)
    Element(14, "Si", "SI", "Silicon", 28.0855)
    Element(15, "P", "P", "Phosphorus", 30.97376)
    Element(16, "S", "S", "Sulfur", 32.06)
    Element(17, "Cl", "CL", "Chlorine", 35.453)
    Element(17, "Cl", "CLA", "Chlorine", 35.453)
    Element(18, "Ar", "AR", "Argon", 39.948)
    Element(19, "K", "K", "Potassium", 39.0983)
    Element(19, "K", "POT", "Potassium", 39.0983)
    Element(20, "Ca", "CAL", "Calcium", 40.08)
    Element(21, "Sc", "SC", "Scandium", 44.9559)
    Element(22, "Ti", "TI", "Titanium", 47.90)
    Element(23, "V", "V", "Vanadium", 50.9415)
    Element(24, "Cr", "CR", "Chromium", 51.996)
    Element(25, "Mn", "MN", "Manganese", 54.9380)
    Element(26, "Fe", "FE", "Iron", 55.847)
    Element(27, "Co", "CO", "Cobalt", 58.9332)
    Element(28, "Ni", "NI", "Nickel", 58.70)
    Element(29, "Cu", "CU", "Copper", 63.546)
    Element(30, "Zn", "ZN", "Zinc", 65.38)
    Element(31, "Ga", "GA", "Gallium", 69.72)
    Element(32, "Ge", "GE", "Germanium", 72.59)
    Element(33, "As", "AS", "Arsenic", 74.9216)
    Element(34, "Se", "SE", "Selenium", 78.96)
    Element(35, "Br", "BR", "Bromine", 79.904)
    Element(36, "Kr", "KR", "Krypton", 83.80)
    Element(37, "Rb", "RB", "Rubidium", 85.4678)
    Element(38, "Sr", "SR", "Strontium", 87.62)
    Element(39, "Y", "Y", "Yttrium", 88.9059)
    Element(40, "Zr", "ZR", "Zirconium", 91.22)
    Element(41, "Nb", "NB", "Niobium", 92.9064)
    Element(42, "Mo", "MO", "Molybdenum", 95.94)
    Element(43, "Tc", "TC", "Technetium", 98.0)
    Element(44, "Ru", "RU", "Ruthenium", 101.07)
    Element(45, "Rh", "RH", "Rhodium", 102.9055)
    Element(46, "Pd", "PD", "Palladium", 106.4)
    Element(47, "Ag", "AG", "Silver", 107.868)
    Element(48, "Cd", "CAD", "Cadmium", 112.41)
    Element(49, "In", "IN", "Indium", 114.82)
    Element(50, "Sn", "SN", "Tin", 118.69)
    Element(51, "Sb", "SB", "Antimony", 121.75)
    Element(52, "Te", "TE", "Tellurium", 127.60)
    Element(53, "I", "I", "Iodine", 126.9045)
    Element(54, "Xe", "XE", "Xenon", 131.30)
    Element(55, "Cs", "CES", "Cesium", 132.9054)
    Element(56, "Ba", "BA", "Barium", 137.33)
    Element(57, "La", "LA", "Lanthanum", 138.9055)
    Element(58, "Ce", "CER", "Cerium", 140.12)
    Element(59, "Pr", "PR", "Praseodymium", 140.9077)
    Element(60, "Nd", "ND", "Neodymium", 144.24)
    Element(61, "Pm", "PM", "Promethium", 145)
    Element(62, "Sm", "SM", "Samarium", 150.4)
    Element(63, "Eu", "EU", "Europium", 151.96)
    Element(64, "Gd", "GD", "Gadolinium", 157.25)
    Element(65, "Tb", "TB", "Terbium", 158.9254)
    Element(66, "Dy", "DY", "Dysprosium", 162.50)
    Element(67, "Ho", "HLM", "Holmium", 164.9304)
    Element(68, "Er", "ER", "Erbium", 167.26)
    Element(69, "Tm", "TM", "Thulium", 168.9342)
    Element(70, "Yb", "YB", "Ytterbium", 173.04)
    Element(71, "Lu", "LU", "Lutetium", 174.967)
    Element(72, "Hf", "HAF", "Hafnium", 178.49)
    Element(73, "Ta", "TA", "Tantalum", 180.9479)
    Element(74, "W", "W", "Tungsten", 183.85)
    Element(75, "Re", "RE", "Rhenium", 186.207)
    Element(76, "Os", "OS", "Osmium", 190.2)
    Element(77, "Ir", "IR", "Iridium", 192.22)
    Element(78, "Pt", "PT", "Platinum", 195.09)
    Element(79, "Au", "AU", "Gold", 196.9665)
    Element(80, "Hg", "MER", "Mercury", 200.59)
    Element(81, "Tl", "TL", "Thallium", 204.37)
    Element(82, "Pb", "PB", "Lead", 207.2)
    Element(83, "Bi", "BI", "Bismuth", 208.9804)
    Element(84, "Po", "PO", "Polonium", 209)
    Element(85, "At", "AT", "Astatine", 210)
    Element(86, "Rn", "RN", "Radon", 222)
    Element(87, "Fr", "FR", "Francium", 223)
    Element(88, "Ra", "RA", "Radium", 226.0254)
    Element(89, "Ac", "AC", "Actinium", 227.0278)
    Element(90, "Th", "TH", "Thorium", 232.0381)
    Element(91, "Pa", "PA", "Protactinium", 231.0359)
    Element(92, "U", "U", "Uranium", 238.029)
]

#
# Retrive index of element in elements list from name. Returns 1 (element "X" of list) if not found
#
function element_index(name::String)
    # Try to find if there is any exact match
    i = findfirst(
        el -> (el.element == name || el.name == name || el.pdb_name == name),
        elements,
    )
    (i != nothing) && return i
    #
    # Try to find by PDB name, note that NT2 and 2NT2 must match N, for example
    #
    # First, check if the first char is a number
    i0 = try
        parse(Int, name[1:1])
        2
    catch
        1
    end
    # Now check if the start of this name matches some PDB code
    for (iel, el) in pairs(elements)
        lpdb = length(el.pdb_name)
        # If the name is shorter than the PDB name, no match
        (length(name) - i0 + 1) < lpdb && continue
        if name[i0:i0+lpdb-1] == el.pdb_name
            return iel
        end
    end
    return 1
end

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
using Distributed
@everywhere using MCPhylo

## Data
mice = Dict{Symbol, Any}(
  :t =>
    [12  1  21 25 11  26  27  30 13 12  21 20  23  25  23  29 35 NaN 31 36
     32 27  23 12 18 NaN NaN  38 29 30 NaN 32 NaN NaN NaN NaN 25  30 37 27
     22 26 NaN 28 19  15  12  35 35 10  22 18 NaN  12 NaN NaN 31  24 37 29
     27 18  22 13 18  29  28 NaN 16 22  26 19 NaN NaN  17  28 26  12 17 26],
  :tcensor =>
    [0 0  0 0 0  0  0  0 0 0  0 0  0  0  0  0 0 40 0 0
     0 0  0 0 0 40 40  0 0 0 40 0 40 40 40 40 0  0 0 0
     0 0 10 0 0  0  0  0 0 0  0 0 24  0 40 40 0  0 0 0
     0 0  0 0 0  0  0 20 0 0  0 0 29 10  0  0 0  0 0 0]
)
mice[:M] = size(mice[:t], 1)
mice[:N] = size(mice[:t], 2)

## Model Specification
model = Model(

  t = Stochastic(2,
    (r, beta, tcensor, M, N) ->
      UnivariateDistribution[(
        lambda = exp(-beta[i] / r);
        0 < lambda < Inf ?
          Truncated(Weibull(r, lambda), tcensor[i, j], Inf) :
          Uniform(0, Inf)) for i in 1:M, j in 1:N
      ],
    false
  ),

  r = Stochastic(
    () -> Exponential(1000)
  ),

  beta = Stochastic(1,
    () -> Normal(0, 10),
    false
  ),

  median = Logical(1,
    (beta, r) -> exp.(-beta / r) * log(2)^(1 / r)
  ),

  veh_control = Logical(
    beta -> beta[2] - beta[1]
  ),

  test_sub = Logical(
    beta -> beta[3] - beta[1]
  ),

  pos_control = Logical(
    beta -> beta[4] - beta[1]
  )

)

## Initial Values
inits = [
  Dict(:t => mice[:t], :beta => fill(-1, mice[:M]), :r => 1.0),
  Dict(:t => mice[:t], :beta => fill(-2, mice[:M]), :r => 1.0)
]

## Sampling Scheme
scheme = [MISS(:t),
          Slice(:beta, 1.0, Univariate),
          Slice(:r, 0.25)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, mice, inits, 20000, burnin=2500, thin=2, chains=2)
describe(sim)


text:
"""
    @traitfn <function definition>

Define a function that contains traits arguments.  The macro expands the function
definition using the Holy Traits pattern such that two functions are defined:

1. A dispatch function that uses [`trait`](@ref) function to determine the type of trait args
2. An implementation function that contains trait type arguments in the front of the signature

# Example

```
@traitfn second(v::Is{Indexable}) = v[2]
```

is expanded to:
```
second(v::T) where T = second(BinaryTraits.trait(Indexable, T), v)
second(::Positive{Indexable}, v::T) where T = v[2]
```
"""
macro traitfn(ex::Expr)
    mod = __module__
    def = splitdef(ex)

    trait_args = find_trait_args(mod, def[:args])
    length(trait_args) > 0 || error("Unable to find any trait arguments")

    my_arg_names = arg_names(def[:args])
    holy_def = move_trait_args(mod, def, trait_args)

    dfunc = make_dispatch_function(holy_def, trait_args, my_arg_names)
    ifunc = make_implementation_function(holy_def, trait_args)
    expr = quote
        $dfunc
        $ifunc
    end
    display_expanded_code(expr)
    return esc(expr)
end

"""
    make_dispatch_function

Make the dispatch function.  It needs to make `trait` method call and pass keyword arguments
to the implementation function.
"""
function make_dispatch_function(def::AbstractDict, trait_args::AbstractVector, my_arg_names::Vector{Symbol})

    def = copy(def)  # make a copy to avoid side effects

    # Make sure that all arguments has a name
    for (i, arg) in enumerate(def[:args])
        if is_typed_arg(arg) && length(arg.args) == 1
            pushfirst!(arg.args, my_arg_names[i])
        end
    end

    # Start building the body expression
    body = Expr(:call, def[:name])

    # Add kwargs if needed.
    # Translate orignal form (k1=v1, k2=v2) into (k1=k1, k2=k2) since we're passing
    # to the implementation function.
    if haskey(def, :kwargs)
        kwargs = Any[Expr(:kw, k.args[1], k.args[1]) for k in def[:kwargs]]
        push!(body.args, Expr(:parameters, kwargs...))
    end

    # Make BinaryTraits.trait() method calls
    trait_exprs = [make_trait_expr(nt.type, nt.param) for nt in trait_args]
    append!(body.args, trait_exprs)

    # Append the regular function arguments
    append!(body.args, my_arg_names)

    def[:body] = body
    return combinedef(def)
end

"""
    make_implementation_function

Make the implementation function, which contains trait type args in the front.
"""
function make_implementation_function(def::AbstractDict, trait_args::AbstractVector)
    def = copy(def)  # make a copy to avoid side effects
    arg_exprs = [:(::$(nt.type)) for nt in trait_args]
    prepend!(def[:args], arg_exprs)
    return combinedef(def)
end

"""
    move_trait_args(mod::Module, def::AbstractDict, trait_args::AbstractVector)

Move the trait args to the where-clause.
"""
function move_trait_args(mod::Module, def::AbstractDict, trait_args::AbstractVector)
    # Make a copy so that we don't mess with the original
    def = deepcopy(def)

    # replace trait arg types with parameter names e.g. x::Can{Fly} ==> x::T1
    for nt in trait_args
        arg = def[:args][nt.pos]                # locate arg that needs to be updated
        argpos = length(arg.args) > 1 ? 2 : 1   # x::T or ::T
        arg.args[argpos] = nt.param             # replace RHS with the param
    end

    # add where parameters e.g. where {T1,T2,...}
    if !haskey(def, :whereparams)
        def[:whereparams] = []
    end
    for nt in trait_args
        push!(def[:whereparams], nt.param)
    end

    return def
end

"""
    arg_names(args::AbstractVector)

Return the function argument names.  If name is missing, a default name is
generated.
"""
function arg_names(args::AbstractVector)
    return map(enumerate(args)) do (i, arg)
        if arg isa Symbol
            arg
        elseif is_typed_arg(arg)
            if length(arg.args) > 1
                arg.args[1]
            else
                Symbol("#ARG", i)
            end
        else
            error("Impossible path, perhaps an upstream language bug? arg=$arg")
        end
    end
end

"""
    find_trait_args(mod::Module, args::Vector)

Find trait arguments from the provided arguments.
Return an array of named tuples with the `name` and `type`,
and position `pos` of the trait arguments.
"""
function find_trait_args(mod::Module, args::Vector)
    trait_args = NamedTuple{(:name, :type, :pos, :param),Tuple{Symbol,DataType,Int,Symbol}}[]
    for (i, arg) in enumerate(args)
        if is_typed_arg(arg)                     # typed argument e.g. x::T, ::T
            if length(arg.args) > 1
                name = arg.args[1]               # extract x from x::T
                type = arg.args[2]               # extract T from x::T
            else
                name = Symbol("#ARG", i)         # derive argument name as ARG<pos>
                type = arg.args[1]               # extract T from ::T
            end
            try
                T = Base.eval(mod, type)         # eval at the user module to find the data type
                if T <: BinaryTrait
                    push!(trait_args, (name = name, type = T, pos = i, param = Symbol("#T", i)))
                end
            catch
                # It could be a user-provided type parameter (e.g. x::S) and S is
                # in the where-clause later.  No need to do anything here since we
                # know it is not a trait arg.
            end
        end
    end
    return trait_args
end

"""
    make_trait_expr

Make an expression that calls `BinaryTraits.trait` function with the specified
trait type `T`.  This expression will be used in a parametric method so a
parameter name is needed.

# Example

```
julia> make_trait_expr(Can{Fly}, :T1)
:(BinaryTraits.trait(Fly, var"#T1"))
```
"""
function make_trait_expr(T::DataType, param::Symbol)
    trait = T.parameters[1]   # extract Fly from Can{Fly}
    return :( BinaryTraits.trait($trait, $param) )
end

"""
    is_typed_arg(arg::Union{Symbol,Expr})

Returns `true` for typed argument expression.
"""
is_typed_arg(arg::Union{Symbol,Expr}) = arg isa Expr && arg.head == :(::)


text:
import NPZ
import MLDataPattern

"""
    load_numpy(path)

Load a numpy array and transpose it so that it is
in column-major form. This lets us make use of the
same data files for our models in Python.

Say you created a numpy array created by, and saved
it as `my_array.npy`:

```python
np.arange(10).reshape((2, 5))

>>> array([[0, 1, 2, 3, 4],
           [5, 6, 7, 8, 9]])
```

loading that array with `load_numpy("my_array.npy")`,
we would get

```julia
5×2 LinearAlgebra.Transpose{Int64,Array{Int64,2}}:
 0  5
 1  6
 2  7
 3  8
 4  9
```
"""
function load_numpy(path)
    arry = NPZ.npzread(path)
    transpose(arry)
end

"""
    split_data(x, y, ratio=0.85)

Split a dataset of x, y pairs into two sets with a given ratio.
This will return a pair of tuples with the `ratio` percent of
of the data belonging to the first tuple:

```
(x_split1, y_split1), (x_split2, y_split2)
```

# Examples
```julia-repl
julia> x = transpose(reshape(collect(1:10), (2,5)))
julia> y = transpose(reshape(collect(0:1), (2,1)))
julia> split_data(x, y, 0.5)
(([1; 3; … ; 7; 9], [0]), ([2; 4; … ; 8; 10], [1]))
```
"""
function split_data(x, y, ratio=0.85)
    MLDataPattern.splitobs((x, y); at = ratio)
end


text:
module MassInstallAction

include("types.jl")

include("public.jl")

include("default_workflows.jl")
include("repo_url.jl")
include("utils.jl")

end # module


text:
using Omega
using UnicodePlots

weight = β(2.0,2.0)
beta_samples = rand(weight, 10000)

UnicodePlots.histogram(beta_samples)

nflips = 4
coinflips_ = [bernoulli(weight, Bool) for i = 1:nflips]


text:
# Routines for Laguerre polynomials

# Compute the expansion of the orthonormal polynomial without e^(qm*x^m/2) nor a constant factor based on some heuristics.
function polyAsyRHgen(np, y, alpha, T::Int64, qm, m::Int64, UQ)
    if (qm == 1) && (m == 1)
        z = y/4/np
        mnxi = 2*np*( sqrt(z)*sqrt(1-z) - acos(sqrt(z) ) ) # = -n*xin/i
    else
        A = zeros(m+1)
        for k =0:m
            A[k+1] = prod((2*(1:k)-1)/2/(1:k))
        end
        z = y/(np*2/m/qm/A[m+1] )^(1/m)
        # Also correct but much slower: Hn = 4*m/(2*m-1)*double(hypergeom([1, 1-m], 3/2-m, z))/m
        Hn = 2/A[m+1]*sum(z.^(0:m-1).*A[m-(0:m-1)])/m
        mnxi = np*(sqrt(z+0im).*sqrt(1-z+0im).*Hn/2 -2*acos(sqrt(z+0im)))
    end
    # We could avoid these tests by splitting the loop k=1:mn into three parts with heuristics for the bounding indices.
    if y < sqrt(np)
        # The fixed delta in the Riemann-Hilbert Problem would mean this bound has to be proportional to n, but x(1:k) are O(1/n) so choose the bound in between them to make more use of the (cheap) expansion in the bulk.
        return asyBesselgen(np, z, alpha, T, qm, m, UQ, mnxi + pi*np, true)
    elseif y > 3.7*np
        # Use the expansion in terms of the (expensive) Airy function, although the corresponding weights will start underflowing for n >= 186 for standard associated Laguerre polynomials.
        return asyAirygen(np, z, alpha, T, qm, m, UQ, (mnxi*3im/2)^(2/3), true)
    end
    asyBulkgen(np, z, alpha, T, qm, m, UQ, mnxi)
end

function asyBulkgen(np, z, alpha, T::Int64, qm, m::Int64, UQ, mnxi)
    if T == 1
        return real( 2/(z+0im)^(1/4 + alpha/2)/(1 - z+0im)^(1/4)*cos(acos(2*z - 1+0im)*(1/2 + alpha/2) - mnxi - pi/4) )
    end
    R = [1 0]
    for k = 1:T-1
        for i = 1:ceil(Int64, 3*T/2)
            R = R + reshape((UQ[1,:,k,i,1]/(z-1)^i + UQ[1,:,k,i,2]/z^i)/np^k,(1,2))
        end
    end
    p = real( 2/(z+0im)^(1/4 + alpha/2)*(cos(acos(2*z-1+0im)*(1/2+alpha/2) - mnxi-pi/4)*R[1]-cos(acos(2*z-1+0im)*(-1/2+alpha/2)-mnxi-pi/4)*R[2]*1im*4^alpha)/(1 - z+0im)^(1/4) )
end

function asyBesselgen(np, z, alpha, T::Int64, qm, m::Int64, UQ, npb, useQ::Bool)
    if T == 1
        return real( sqrt(2*pi)*(-1)^np*sqrt(npb)/(z+0im)^(1/4+alpha/2)/(1 - z+0im)^(1/4)*(sin( (alpha + 1)/2*acos(2*z - 1+0im) - pi*alpha/2)*besselj(alpha,npb) + cos( (alpha + 1)/2*acos(2*z - 1+0im) - pi*alpha/2)*(besselj(alpha-1,npb) - alpha/(npb)*besselj(alpha, npb) ) ) )
    end
    R = (1+0im)*[1 0]
    # Use the series expansion of R because it is faster and we use asyBessel only very close to zero to have less calls to besselj.
    if useQ
        for k = 1:T-1
            for i = 1:min(size(UQ,4),9-k)
                R = R + reshape(UQ[1, :, k, i, 4]*z^(i-1)/np^k,(1,2))
            end
        end
    else
        d = z-1+0im
        phi = 2*z-1+2*sqrt(z)*sqrt(d)
        Rko = (1+0im)zeros(T-1,2)
        sL = (1+0im)zeros(2,2,T-1)
		for m = 1:T-1
            for i = 1:ceil(Int64,3*m/2)
                Rko[m,:] += UQ[1, :, m, i, 1]/d^i+UQ[1, :, m, i, 2]/z^i
            end
            sL[:,:,m] = brac(m-1,alpha)/2^(1+2*m)/(npb/2im/np)^m*( [2^(-alpha)  0 ; 0   2^(alpha)]*[sqrt(phi)    1im/sqrt(phi)  ;  -1im/sqrt(phi)   sqrt(phi)]/2/z^(1/2)/d^(1/2)*[(-phi)^(alpha/2)  0 ; 0   (-phi)^(-alpha/2) ]*[((-1)^m)/m*(alpha^2+m/2-1/4)     (m-1/2)*1im  ;  -((-1)^m)*(m-1/2)*1im    (alpha^2+m/2-1/4)/m]*[(-phi)^(-alpha/2)  0 ; 0   (-phi)^(alpha/2) ]*[sqrt(phi)     -1im/sqrt(phi) ;  1im/sqrt(phi)    sqrt(phi)]*[2^(alpha)   0 ; 0 2^(-alpha)] -mod(m+1,2)*(4*alpha^2+2*m-1)/m*eye(2,2) )
        end
        for k = 1:T-1
            R += reshape((Rko[k,:] - sL[1,:,k])/np^k,(1,2))
            for m = 1:k-1
                R -= reshape(reshape(Rko[k-m,:],(1,2))*sL[:,:,m]/np^k,(1,2))
            end
        end
    end
    p = real( sqrt(2*pi)*(-1)^np*sqrt(npb)/(z+0im)^(1/4+alpha/2)/(1 - z+0im)^(1/4)*( (sin( (alpha + 1)/2*acos(2*z - 1+0im) - pi*alpha/2)*R[1] -sin( (alpha - 1)/2*acos(2*z - 1+0im) - pi*alpha/2)*R[2]*1im*4^alpha)*besselj(alpha, npb) + (cos( (alpha + 1)/2*acos(2*z - 1+0im)- pi*alpha/2)*R[1] - cos( (alpha - 1)/2*acos(2*z - 1+0im) - pi*alpha/2)*R[2]*1im*4^alpha)*(besselj(alpha-1, npb) - alpha/npb*besselj(alpha, npb) ) ) )
end

function asyAirygen(np, z, alpha, T::Int64, qm, m::Int64, UQ, fn, useQ::Bool, xin=NaN+NaN*1im)
    d = z - 1.0 +0im
    if T == 1
        return real( 4*sqrt(pi)/(z+0im)^(1/4+alpha/2)/d^(1/4)*(cos( (alpha + 1)/2*acos(2*z - 1+0im) )*fn^(1/4)*airyai(fn) -1im*sin( (alpha + 1)/2*acos(2*z - 1+0im) )*ifelse(angle(z-1) <= 0, -one(z), one(z) )*fn^(-1/4)*airyaiprime(fn) ) )
    end
    R = (1+0im)*[1 0]
    if useQ
        for k = 1:T-1
            for i = 1:min(size(UQ,4),9-k)
                R = R + reshape(UQ[1, :, k, i, 3]*d^(i-1)/np^k,(1,2))
            end
        end
    else
        phi = 2*z-1+2*sqrt(z)*sqrt(d)
        Rko = (1+0im)zeros(T-1,2)
        sR = (1+0im)zeros(2,2,T-1)
	for m = 1:T-1
            for i = 1:ceil(Int64,3*m/2)
                Rko[m,:] += UQ[1, :, m, i, 1]/d^i+UQ[1, :, m, i, 2]/z^i
            end
            sR[:,:,m] = nuk(m)/xin^m*( [2^(-alpha)  0 ; 0   2^(alpha)]*[sqrt(phi)    1im/sqrt(phi)  ;  -1im/sqrt(phi)   sqrt(phi)]/8/z^(1/2)/d^(1/2)*[phi^(alpha/2)  0 ; 0   phi^(-alpha/2) ]*[(-1.0)^m  -m*6im ; 6im*m*(-1)^m    1.0]*[phi^(-alpha/2)  0 ; 0   phi^(alpha/2) ]*[sqrt(phi)     -1im/sqrt(phi) ;  1im/sqrt(phi)    sqrt(phi)]*[2^(alpha)   0 ; 0 2^(-alpha)] - mod(m+1,2)*eye(2,2) )
        end
        for k = 1:T-1
            R += reshape((Rko[k,:] -sR[1,:,k])/np^k,(1,2))
            for m = 1:k-1
                R -= reshape(reshape(Rko[k-m,:],(1,2))*sR[:,:,m]/np^k,(1,2))
            end
        end
    end
    p = real( 4*sqrt(pi)/(z+0im)^(1/4+alpha/2)/d^(1/4)*( (R[1]*cos( (alpha + 1)/2*acos(2*z - 1+0im) ) -cos( (alpha - 1)/2*acos(2*z - 1+0im) )*R[2]*1im*4^alpha)*fn^(1/4)*airyai(fn) + 1im*(-sin( (alpha + 1)/2*acos(2*z - 1+0im) )*R[1] +sin( (alpha - 1)/2*acos(2*z - 1+0im) )*R[2]*1im*4^alpha)*ifelse(angle(z-1) <= 0, -one(z), one(z) )*fn^(-1/4)*airyaiprime(fn) ) )
end

# Additional short functions
# pochhammer does not seem to exist yet in Julia
poch(x,n) = prod(x .+ (0:(n-1)))

function binom(x,n) # binomial only works for integer x
    b = one(x)
    for i = 1:n
        b *= (x-(n-i))/i
    end
    b
end

nuk(n) = -gamma(3*n-1/2)*2^n/27^n/2/n/sqrt(pi)/gamma(n*2)
brac(n,alpha) = prod(4*alpha^2 .- (2*(1:n) .- 1).^2 )/(2^(2*n)*gamma(1.0+n))

# Compute the W or V-matrices to construct the asymptotic expansion of R.
# Input
#   alpha, qm, m - Factors in the weight function w(x) = x^alpha*exp(-qm*x^m)
#   maxOrder     - The maximum order of the error
#   r            - 1 when computing Wright, -1 when computing Wleft
#   isW          - Whether we compute W(iso V)-matrices
# Output
#   WV           - Coefficient matrices for (z + 1/2 \pm 1/2)^m of Delta_k(z) or s_k(z)
function getV(alpha,qm,m::Int64,maxOrder::Int64,r)
    mo = ceil(Int64, 3*maxOrder/2) + 4
    ns = 0:mo
    f = NaN*zeros(mo+1) # Coefficients in the expansion of \bar{phi}_n(z) or \xi_n(z)
    g = NaN*zeros(maxOrder-1,mo+1)

    A = zeros(m+1)
    for k =0:m
        A[k+1] = prod((2*(1:k) .- 1.0)/2/(1:k))
    end
    if (r == 1) # Right disk: near z=1
        f = NaN*zeros(mo+2)
        ns = [ns; ns[mo+1]+1] # Extend by one because f(1) = 0 while not for left
        u = zeros(ns[mo+2]+1,ns[mo+2]+2)
        v = zeros(ns[mo+2]+1,ns[mo+2]+2)
        u[1,1] = 1.0
        v[1,1] = 1.0
        for n = [ns; ns[mo+2]+1]
            u[2,n+1] = binom(1/2,n+1)
            v[2,n+1] = binom(1/2,n+2)
        end
        for kt = 2:ns[mo+2]
            for n = ns
                u[kt+1,n+1] = sum(u[kt,(0:n) .+ 1].*u[2,n .- (0:n) .+ 1])
                v[kt+1,n+1] = sum(v[kt,(0:n) .+ 1].*v[2,n .- (0:n) .+ 1])
            end
        end
        q = zeros(ns[mo+2]+1)
        rr = zeros(ns[mo+2]+1) # Coeffs in the expansion of sqrt(2-2*sqrt(1-w))
        for kt = ns
            for l = 0:kt
                q[kt+1] = q[kt+1] + poch(1/2,kt-l)*u[kt-l+1,l+1]/(-2)^(kt-l)/gamma(1.0+kt-l)/(1+2*(kt-l) )
                rr[kt+1] = rr[kt+1] + binom(1/2,kt-l)*v[kt-l+1,l+1]*2^(kt-l)
            end
        end
        if (m == 1)
            for n = ns
                f[n+1,1] = -2*binom(1/2,n)
                for l = 0:n
                    f[n+1,1] = f[n+1,1] + 2*q[l+1]*rr[n-l+1]
                end
            end
        else
            for j = ns
                f[j+1,1] = 0
                for i=0:min(j,m-1)
                    f[j+1,1] = f[j+1,1] + binom(1/2,j-i)*(-1)^(m-i-1)*gamma(-1/2-i)/gamma(1/2-m)/gamma(m-i)
                end
                f[j+1,1] = -f[j+1,1]/m/A[m+1]
                for l = 0:j
                    f[j+1,1] = f[j+1,1] + 2*q[l+1]*rr[j-l+1]
                end
            end
        end
        if(abs(f[1]) > 10*eps(Float64) )
            error("xi_n should be O( (z-1)^(3/2) ): Expected f[1] to be zero")
        end
        ns = ns[1:mo+1] # Reset ns to its value before computing f's
        g[1,1,1] = -1/f[2,1]
        for n = 1:mo
            g[1,n+1,1] = -sum(reshape(g[1,1:n,1],(n,1)).*reshape(f[(n+2):-1:3,1],(n,1) ) )/f[2,1]
        end
    else # Left disk: near z=0
        if (m == 1)
            for n = ns
                f[n+1,1] = -(binom(1/2,n)*(-1)^n + poch(1/2,n)./(1+2*n)./gamma(1.0+n))
            end
        else
            for n = ns
                f[n+1,1] = 0.0
                for k = 0:min(m-1,n)
                    f[n+1,1] = f[n+1,1] + binom(1/2,n-k)*(-1)^(n-k)*A[m-k]
                end
                f[n+1,1] = -f[n+1,1]/2/m/A[m+1]-poch(1/2,n)./(1+2*n)./gamma(1.0+n)
            end
        end
        g[1,1,1] = 1/f[1,1]
        for n = 1:mo
            g[1,n+1,1] = -sum(reshape(g[1,1:n,1],(n,1)).*reshape(f[(n+1):-1:2,1],(n,1)) )/f[1,1]
        end
    end
    rho = (1+1im)*zeros(2*mo+3,mo+1)
    for n = ns
        rho[2,n+1] = poch(1/2,n)/gamma(1.0+n)/(1+2*n)*(-r)^n
    end
    rho[1,1] = 1
    for i = 2:(maxOrder-1)
        for n = ns
            g[i,n+1] = sum(g[i-1,1:(n+1) ].*g[1,(n+1):-1:1] )
        end
    end
    for i = 2:(mo*2+2)
        for n = ns
            rho[i+1,n+1] = sum(rho[i,1:(n+1) ].*rho[2,(n+1):-1:1] )
        end
    end
    OmOdd = (1+1im)*zeros(mo+1); OmEven = (1+1im)*zeros(mo+1)
    XiOdd = (1+1im)*zeros(mo+1); XiEven = (1+1im)*zeros(mo+1)
    ThOdd = (1+1im)*zeros(mo+1); ThEven = (1+1im)*zeros(mo+1)
    OmO = (1+1im)*zeros(mo+1); OmE = (1+1im)*zeros(mo+1)
    XiO = (1+1im)*zeros(mo+1); XiE = (1+1im)*zeros(mo+1)
    ThO = (1+1im)*zeros(mo+1); ThE = (1+1im)*zeros(mo+1)
    for n = ns
        js = 0:n
        for j = js
            OmOdd[n+1] = OmOdd[n+1] + (-1)^j/gamma(1.0+2.0*j)*(-2*alpha/sqrt(-r+0.0im))^(2*j)*rho[2*j+1,n-j+1]
            XiOdd[n+1] = XiOdd[n+1] + (-1)^j/gamma(1.0+2.0*j)*(-2*(alpha+1)/sqrt(-r+0.0im))^(2*j)*rho[2*j+1,n-j+1]
            ThOdd[n+1] = ThOdd[n+1] + (-1)^j/gamma(1.0+2.0*j)*(-2*(alpha-1)/sqrt(-r+0.0im))^(2*j)*rho[2*j+1,n-j+1]
            OmEven[n+1] = OmEven[n+1] + (-1)^j/gamma(1.0+2*j+1.0)*(-2*alpha/sqrt(-r+0.0im))^(2*j+1)*rho[2*j+2,n-j+1]
            XiEven[n+1] = XiEven[n+1] + (-1)^j/gamma(1.0+2*j+1.0)*(-2*(alpha+1)/sqrt(-r+0.0im))^(2*j+1)*rho[2*j+2,n-j+1]
            ThEven[n+1] = ThEven[n+1] + (-1)^j/gamma(1.0+2*j+1.0)*(-2*(alpha-1)/sqrt(-r+0.0im))^(2*j+1)*rho[2*j+2,n-j+1]
        end
        for j = js
            OmO[n+1] = OmO[n+1] + binom(-1/2,j)*(r)^j*OmOdd[n-j+1]
            XiO[n+1] = XiO[n+1] + binom(-1/2,j)*(r)^j*XiOdd[n-j+1]
            ThO[n+1] = ThO[n+1] + binom(-1/2,j)*(r)^j*ThOdd[n-j+1]
            OmE[n+1] = OmE[n+1] + binom(-1/2,j)*(r)^j*OmEven[n-j+1]
            XiE[n+1] = XiE[n+1] + binom(-1/2,j)*(r)^j*XiEven[n-j+1]
            ThE[n+1] = ThE[n+1] + binom(-1/2,j)*(r)^j*ThEven[n-j+1]
        end
    end
    Ts = zeros(ComplexF64, 2, 2, mo+1) # = G_{k,n}^{odd/even} depending on k, overwritten on each new k
    WV = zeros(ComplexF64, 2, 2, maxOrder-1, mo+1)
    for k = 1:(maxOrder-1)
        Ts[:,:,:] .= 0
        if r == 1
            if mod(k,2) == 1
                for n = 0:mo
                    Ts[:,:,n+1] = nuk(k)*[-2*(2*binom(-1/2,n-1)*(n>0)+binom(-1/2,n))     2im*4^(-alpha)*binom(-1/2,n)    ;    2im*4^(alpha)*binom(-1/2,n)    (2*(2*binom(-1/2,n-1)*(n>0) +binom(-1/2,n)))] -6*k*nuk(k)*[-2*OmO[n+1]   4^(-alpha)*2im*XiO[n+1]  ;   4^(alpha)*2im*ThO[n+1]    2*OmO[n+1]]
                    WV[:,:,k,n+1] = sum(repeat(reshape(g[k,1:(n+1)], (1,1,n+1) ), outer=(2,2,1)) .* Ts[:,:,(n+1):-1:1], dims=3)/8
                end
            else
                for n = 0:mo
                     Ts[:,:,n+1] = nuk(k)*4*(n==0)*I +6*k*nuk(k)*[-2im*OmE[n+1]    -2*4^(-alpha)*XiE[n+1]  ;   -2*4^alpha*ThE[n+1]   2im*OmE[n+1]]
                     WV[:,:,k,n+1] = sum(repeat(reshape(g[k,1:(n+1) ], (1,1,n+1) ), outer=[2,2,1]).*Ts[:,:,(n+1):-1:1], dims=3)/8
                end
            end
        else
            if mod(k,2) == 1
                for n = 0:mo
                    Ts[:,:,n+1] = -(alpha^2+k/2-1/4)/k*[-(-1)^n*(2*binom(-1/2,n-1)*(n>0)+binom(-1/2,n))*2    -1im*4^(-alpha)*2*(-1)^n*binom(-1/2,n)  ;  -1im*4^(alpha)*2*(-1)^n*binom(-1/2,n)     ( (-1)^n*(2*binom(-1/2,n-1)*(n>0) +binom(-1/2,n))*2)] - (k-1/2)*[2*OmO[n+1]   4^(-alpha)*2im*XiO[n+1]  ;   4^(alpha)*2im*ThO[n+1]   -2*OmO[n+1]] # binom(-1/2,-1) should be zero
                    WV[:,:,k,n+1] = -(-1)^(ceil(Int64, k/2)+1)*(1im*sqrt(2))^k*(-2+0im)^(-k/2)/4^(k+1)*brac(k-1,alpha)*sum(repeat(reshape(g[k,1:(n+1) ], (1,1,n+1) ), outer=[2,2,1]).*Ts[:,:,(n+1):-1:1], dims=3)
                end
            else
                for n = 0:mo
                    Ts[:,:,n+1] = (alpha^2+k/2-1/4)/k*4*(n==0)*I  -2*(k-1/2)*[ OmE[n+1]   4^(-alpha)*1im*XiE[n+1]  ;   4^alpha*1im*ThE[n+1]   -OmE[n+1] ]
                    WV[:,:,k,n+1] = -(-1)^(ceil(Int64, k/2)+1)*(1im*sqrt(2))^k*(-2)^(-k/2)/4^(k+1)*brac(k-1,alpha)*sum(repeat(reshape(g[k,1:(n+1) ], (1,1,n+1) ), outer=[2,2,1]).*Ts[:,:,(n+1):-1:1],dims=3)
                end
            end
        end
    end
    WV
end

# Get the U-matrices to construct the asymptotic expansion of R using the procedure with the convolutions with a specified method.
# Input
#   alpha, qm, m - Parts of the weight function
#   maxOrder     - The maximal order of the error
# Output
#   UQ           - Coefficient matrices of R_k(z) for (z-1)^(-m) [Uright], or z^(-m) [Uleft] of R_k^{right}(z) for (z-1)^n [Qright] and of R_k^{left}(z) for z^n [Qleft]
function getUQ(alpha, qm, m::Int64, maxOrder::Int64)
    Vr = getV(alpha, qm, m, maxOrder, 1)
    Vl = getV(alpha, qm, m, maxOrder, -1)
    UQ = (1+1im)*zeros(2,2,maxOrder-1,ceil(Int64,3*maxOrder/2)+2, 4)
    for kt = 0:(maxOrder-2)
        # Uright(:,:,(maxOrder-1)+1,:) will not be used later on because first term in expansions is without U's
        for mt = 0:(ceil(Int64,3*(kt+1)/2)-1)
            UQ[:,:,kt+1,mt+1,1] = Vr[:,:,kt+1,ceil(Int64,3*(kt+1)/2)-mt]
            for j = 0:(kt-1)
                for l = 0:(ceil(Int64,3*(j+1)/2)-mt-1)
                    UQ[:,:,kt+1,mt+1,1] = UQ[:,:,kt+1,mt+1,1] + UQ[:,:,kt-j,l+1,3]*Vr[:,:,j+1,ceil(Int64,3*(j+1)/2)-l-mt]
                end
            end
        end
        for mt = 0:(ceil(Int64,(kt+1)/2)-1)
            UQ[:,:,kt+1,mt+1,2] = Vl[:,:,kt+1,ceil(Int64,(kt+1)/2)-mt]
            for j= 0:(kt-1)
                for l = 0:(ceil(Int64,(j+1)/2)-mt-1)
                    UQ[:,:,kt+1,mt+1,2] = UQ[:,:,kt+1,mt+1,2] + UQ[:,:,kt-j,l+1,4]*Vl[:,:,j+1,ceil(Int64,(j+1)/2)-l-mt]
                end
            end
        end
        for n = 0:(ceil(Int64,3*(maxOrder-kt+1)/2)-1)
            UQ[:,:,kt+1,n+1,3] = -Vr[:,:,kt+1,ceil(Int64,3*(kt+1)/2)+1+n]
            UQ[:,:,kt+1,n+1,4] = -Vl[:,:,kt+1,ceil(Int64,(kt+1)/2)+1+n]
            for i = 0:(ceil(Int64,(kt+1)/2)-1)
                UQ[:,:,kt+1,n+1,3] = UQ[:,:,kt+1,n+1,3] + binom(-i-1,n)*UQ[:,:,kt+1,i+1,2]
            end
            for i = 0:(ceil(Int64,3*(kt+1)/2)-1)
                UQ[:,:,kt+1,n+1,4] = UQ[:,:,kt+1,n+1,4] + binom(-i-1,n)*(-1.0)^(-i-1-n)*UQ[:,:,kt+1,i+1,1]
            end
            for j = 0:(kt-1)
                for l = 0:(ceil(Int64, (j+1)/2)+n)
                    UQ[:,:,kt+1,n+1,4] = UQ[:,:,kt+1,n+1,4] -UQ[:,:,kt-j,l+1,4]*Vl[:,:,j+1,n-l+1+ceil(Int64,(j+1)/2) ]
                end
                for l = 0:(ceil(Int64, 3*(j+1)/2)+n)
                    UQ[:,:,kt+1,n+1,3] = UQ[:,:,kt+1,n+1,3] -UQ[:,:,kt-j,l+1,3]*Vr[:,:,j+1,n-l+1+ceil(Int64, 3*(j+1)/2) ]
                end
            end
        end
    end
    UQ
end


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

mutable struct KeyVaultSecretReference <: SwaggerModel
    secretUrl::Any # spec type: Union{ Nothing, String } # spec name: secretUrl
    sourceVault::Any # spec type: Union{ Nothing, SubResource } # spec name: sourceVault

    function KeyVaultSecretReference(;secretUrl=nothing, sourceVault=nothing)
        o = new()
        validate_property(KeyVaultSecretReference, Symbol("secretUrl"), secretUrl)
        setfield!(o, Symbol("secretUrl"), secretUrl)
        validate_property(KeyVaultSecretReference, Symbol("sourceVault"), sourceVault)
        setfield!(o, Symbol("sourceVault"), sourceVault)
        o
    end
end # type KeyVaultSecretReference

const _property_map_KeyVaultSecretReference = Dict{Symbol,Symbol}(Symbol("secretUrl")=>Symbol("secretUrl"), Symbol("sourceVault")=>Symbol("sourceVault"))
const _property_types_KeyVaultSecretReference = Dict{Symbol,String}(Symbol("secretUrl")=>"String", Symbol("sourceVault")=>"SubResource")
Base.propertynames(::Type{ KeyVaultSecretReference }) = collect(keys(_property_map_KeyVaultSecretReference))
Swagger.property_type(::Type{ KeyVaultSecretReference }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_KeyVaultSecretReference[name]))}
Swagger.field_name(::Type{ KeyVaultSecretReference }, property_name::Symbol) =  _property_map_KeyVaultSecretReference[property_name]

function check_required(o::KeyVaultSecretReference)
    (getproperty(o, Symbol("secretUrl")) === nothing) && (return false)
    (getproperty(o, Symbol("sourceVault")) === nothing) && (return false)
    true
end

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


text:
module RepeatedSolvers

abstract Solver
include("direct_solver.jl")
include("jacobi_solver.jl")

end


text:
module Inspect
include("Utils.jl")
include("Constants.jl")
include("SurfaceReactions.jl")
using LightXML
using Printf

"""
This function is to test a surface reaction mechanism.
Ideally the sum of fluxes must vanish at steasy state.
The function also gives the surface coverages for a
given gasphase concentration and T,p conditions.
The program generates only a screen outout
#   Usage:
-   run(input_file, lib_dir)
-   input_file: name of the input file (String)
-   lib_dir: path to the folder where the data files are present
"""
function  run(input_file::AbstractString, lib_dir::AbstractString)

    local gasphase = Array{String,1}()

    xmldoc = parse_file(input_file)
    xmlroot = root(xmldoc)
    #Read the gasphase
    #tag_gasphase = get_elements_by_tagname(xmlroot,"gasphase")
    gasphase = get_collection_from_xml(xmlroot,"gasphase")
    thermo_file = lib_dir*"therm.dat"
    thermo_obj = SurfaceReactions.IdealGas.create_thermo(gasphase, thermo_file)

    #Get the molefractions
    mole_fracs = get_molefraction_from_xml(xmlroot,thermo_obj.molwt,gasphase)

    #Read the temperature xml
    local T = get_value_from_xml(xmlroot,"T")
    #Get the pressure from xml
    local p = get_value_from_xml(xmlroot,"p")

    #Get the mechanism file from xml
    local mech_file = get_text_from_xml(xmlroot,"mech")
    mech_file = lib_dir*"/"*mech_file

    md = SurfaceReactions.compile_mech(mech_file,thermo_obj,gasphase)
    covg = md.sm.si.ini_covg
    n_species = length(gasphase)+length(md.sm.species)
    rate = zeros(n_species)
    #Create the state object
    state = SurfaceReactions.State(mole_fracs,covg,T,p,rate)

    #calculate the molar production rates
    t, state = SurfaceReactions.calculate_ss_molar_production_rates!(state,thermo_obj,md,10.0)
    println("\nInitial coverage: ")
    println("-----------------------")
    for i in eachindex(covg)
        @printf("%12s \t %.4e \n", md.sm.species[i], md.sm.si.ini_covg[i])
    end
    println("\nT(K): \t",T)
    println("p(Pa): \t", p)
    println("\nFinal coverage and flux after $t(s): \n")
    println("     Species \t coverage \t flux(mol/m2-s)")
    println("-------------------------------------------------")
    ng = length(state.mole_frac)

    for i in eachindex(covg)
        @printf("%12s \t %.4e \t %+.4e \n",md.sm.species[i], state.covg[i], state.s_rate[ng+i])
    end

    println("\nSum of fluxes (mol/m2-s): ", sum(state.s_rate[ng+1:lastindex(state.s_rate)]))
    println("Sum of coverages: ", sum(state.covg))

    return Symbol("Success")

end

end


text:
@testset "power" begin

    p = Power(3)
    @test p isa Transform
    @test cardinality(p) == OneToOne()

    @test FeatureTransforms.apply([1, 2, 3], p) == [1, 8, 27]

end


text:
using SFML

"""
Interface to image file containing tiles.
"""
type TileSet
    texture::Texture
    size::Tuple{Int, Int}  # number of tiles in each direction
    tilesize::Vector2i # in pixel
    function TileSet(texture::Texture, siz::Tuple{Int, Int}, tilesize::Vector2i)
        texture_size = texture |> get_size |> Vector2i
        calc_size = Vector2i(siz...) * tilesize
        @assert map(<=, calc_size, texture_size) |> all
        new(texture, siz, tilesize)
    end
end

tileset(t::TileSet) = t
tilesize(ts::TileSet) = ts.tilesize
tilesize(t) = t |> tileset |> tilesize
Base.size(ts::TileSet) = ts.size
Base.size(ts::TileSet, i) = size(ts)[i]
tilewidth(ts) = tilesize(ts)[1]
tileheight(ts) = tilesize(ts)[2]

function TileSet(texture::Texture, siz, tilesiz)
    TileSet(texture, tuple(siz...), Vector2i(tilesiz...))
end

function TileSet(filename::String, size_, tilesize)
    TileSet(Texture(filename), size_, tilesize)
end

immutable Tile
    corners::NTuple{4, Vector2f} # pixel coords of corners of tile in TileSet
end

corners(t::Tile) = t.corners

_colors4() = (SFML.white,SFML.white,SFML.white,SFML.white)

function Base.convert(::Type{VertexArray}, t::Tile)
    vertices = map(Vertex, corners(t), _colors4(), corners(t))
    va = VertexArray()
    set_primitive_type(va, SFML.quads)
    for v in vertices
        append(va, v)
    end
    va
end

type TileMap
    tilemat::Matrix{Tile}
    globalbounds::FloatRect
    ts::TileSet
end

function TileMap(tmat::Matrix{Tile}, ts::TileSet)
    width, height = Vector2i(size(tmat)...) * tilesize(ts)
    r = FloatRect(0,0,width, height)
    return TileMap(tmat, r, ts)
end

tilematrix(tm::TileMap) = tm.tilemat
Base.size(tm::TileMap, dims...) = size(tm |> tilematrix, dims...)

function subrect(tm::TileMap,i,j)
    r = tm.globalbounds
    n,m = size(tm)
    width = r.width/n
    height = r.height/m
    left = r.left + (i-1) * width
    top = r.top + (j-1) * height
    FloatRect(left, top, width, height)
end

function Base.convert(::Type{VertexArray}, tm::TileMap)
    va = VertexArray()
    set_primitive_type(va, SFML.quads)
    for i in 1:size(tm,1), j in 1:size(tm,2)
        t = tilematrix(tm)[i,j]
        r = subrect(tm, i, j)
        pushtile!(va, r, t)
    end
    va
end

function pushtile!(va::VertexArray, r::FloatRect, t::Tile)
    # colors = (SFML.white, SFML.white, SFML.white, SFML.white)
    verts = map(Vertex, corners(r), _colors4(), corners(t))
    for v in verts
        append(va, v)
    end
    va
end

function Base.getindex(ts::TileSet, i::Int, j::Int)
    c = setcoords2tilerect(ts, Vector2i(i,j)) |> corners |> Tile
end

function corners(r::Union{FloatRect, IntRect})
    x_lower = r.left
    x_upper = x_lower + r.width
    y_lower = r.top
    y_upper = r.top + r.height
    return (
        Vector2(x_lower, y_upper),
        Vector2(x_lower, y_lower),
        Vector2(x_upper, y_lower),
        Vector2(x_upper, y_upper)
    )
end

function setcoords2tilerect(ts::TileSet, setcoords::Vector2i)
    left, top = (ts |> tilesize) * (setcoords - Vector2i(1,1))
    width, height = (ts |> tilesize)
    return FloatRect(left, top, width, height)
end

ts = TileSet("tiles.png", (7, 5), (17,17))
water = ts[1,1]
wall = ts[2,1]
lava = ts[3,1]
tube = ts[4,5]

m = [water water water;
    lava   tube  tube ;
    lava   lava  lava ;
    wall   wall  wall ;
]

vertex_array = TileMap(m, ts) |> VertexArray
states = RenderStates(ts.texture)

window=RenderWindow("quick window", 800, 600)
event = Event()
while isopen(window)
    while pollevent(window, event)
        (get_type(event) == EventType.CLOSED) && close(window)
    end
    clear(window, SFML.black)
    draw(window, vertex_array, states)
    display(window)
end


text:
function normcol_equal(matin)
# solve the proximal problem
# matout = argmin||matout-matin||_F^2, s.t. matout(:,i)=1
    matin ./ repmat(sqrt(sum(matin .^ 2, 1) + eps()), size(matin, 1), 1)
end


text:
module Data

using Statistics
using Random

import Base.size

using ..Utils:
    flatmap, copy_with_eltype, issomething

export AbstractData, WXData, PlainXData, XData, XYData, XBatches, XYBatches, Dataset, UnlabeledDataset,
LabeledDataset, XDataset, XYDataset, BatchedXDataset, BatchedXYDataset, weights,
num_examples, total_example_weight, max_batch_size, num_features, num_labels, num_batches,
feature_matrix, unweighted_data, feature_data, labels, weights, aggr_weight_type, feature_type, label_type,
train, valid, test,
shuffle, batch, threshold, fully_factorized_likelihood,
ll_per_example, bits_per_pixel, vslice, unbatch

#####################
# Types
#####################

"""
Any form of data where features are of type X
"""
abstract type AbstractData{X} end

"""
Abstract unlabeled data (X-values)
"""
abstract type XData{X,M<:AbstractMatrix{<:X}} <: AbstractData{X} end

"""
Unlabeled data (X-values) with a direct matrix representation
"""
struct PlainXData{X,M<:AbstractMatrix{<:X}} <: XData{X,M} # repeating the type bound here allows for automatic inference of X in constructor
    x::M
end

"""
Unlabeled data (X-values) with weights for each example
"""
struct WXData{X,W,M} <: XData{X,M}
    xd::PlainXData{X,M}
    w::Vector{W}
    WXData(x::PlainXData{X,M}, w::AbstractVector{<:W}) where {X,W,M} =
        (num_examples(x) == length(w)) ? new{X,W,M}(x,w) : error("x and w have different numbers of examples ($(num_examples(x)) vs. $(length(w)))")
end

"Labeled data (X-values with Y-labels)"
struct XYData{X,Y,XD<:XData{<:X}, V<:AbstractVector{<:Y}} <: AbstractData{X}
    xd::XD
    y::V
    XYData(xd::XD,y::V) where {XD,V} =
        (num_examples(xd) == size(y)[1]) ? new{feature_type(xd),eltype(V),XD,V}(xd,y) : error("x and y have different numbers of examples")
end

# batches of data data

"Batches of unlabeled `XData`"
const XBatches{X,XD<:XData{X}} = AbstractVector{<:XD}

"Batches of labeled `XYData`"
const XYBatches{X,Y, XYD<:XYData{X,Y}} = AbstractVector{<:XYD}

# datasets with splits

"A dataset with train/validation/test splits"
abstract type Dataset{X} <: AbstractData{X} end;

"An unlabeled dataset with train/validation/test splits"
abstract type UnlabeledDataset{X} <: Dataset{X} end;

"A labeled dataset with train/validation/test splits"
abstract type LabeledDataset{X,Y} <: Dataset{X} end;

"An single unlabeled dataset with train/validation/test splits"
struct XDataset{X} <: UnlabeledDataset{X}
    train::XData{<:X}
    valid::Union{Nothing,<:XData{<:X}}
    test::XData{<:X}
end

"An single labeled dataset with train/validation/test splits"
struct XYDataset{X,Y} <: LabeledDataset{X,Y}
    train::XYData{<:X,<:Y}
    valid::Union{Nothing,<:XYData{X,Y}}
    test::XYData{<:X,<:Y}
end

"Batches of unlabeled data with train/validation/test splits"
struct BatchedXDataset{X} <: UnlabeledDataset{X}
    train::XBatches{<:X}
    valid::Union{Nothing,<:XBatches{<:X}}
    test::XBatches{<:X}
end

"Batches of labeled data with train/validation/test splits"
struct BatchedXYDataset{X,Y} <: LabeledDataset{X,Y}
    train::XYBatches{<:X,<:Y}
    valid::Union{Nothing,<:XYBatches{<:X,<:Y}}
    test::XYBatches{<:X,<:Y}
end

#####################
# Constructors and conversions
#####################

WXData(x::M, w::AbstractVector{<:W}) where {M,W} = WXData(PlainXData(x), w)
WXData(x::M) where M = WXData(x,ones(num_examples(x)))

XData(x::M) where M = PlainXData(x)

XDataset(xyd::XYDataset{X}) where X = XDataset{X}(feature_data(xyd.train),
                                        map_something(v -> feature_data(v), xyd.valid),
                                        feature_data(xyd.test))

Base.convert(::Type{XBatches}, d::XData) = [d]
Base.convert(::Type{XBatches{X}}, d::XData{X}) where X = [d]
Base.convert(::Type{XYBatches{X,Y}}, d::XYBatches{X,Y}) where {X,Y} = [d]

Base.convert(::Type{BatchedXDataset}, d::XDataset) = BatchedXDataset(d.train, d.valid, d.test)
Base.convert(::Type{BatchedXDataset{X}}, d::XDataset{X}) where X = BatchedXDataset(d.train, d.valid, d.test)
Base.convert(::Type{BatchedXYDataset{X,Y}}, d::XYDataset{X,Y}) where {X,Y} = BatchedXYDataset(d.train, d.valid, d.test)

#####################
# Methods
#####################

"Replace the feature matrix in data by a given matrix"
replace_features(xd::PlainXData, x::M) where M = PlainXData(x)
replace_features(xd::WXData, x::M) where M = WXData(x,weights(xd))
replace_features(xyd::XYData, x::M) where M = XYData(replace_features(feature_data(xyd)),labels(xyd))

"Number of examples in data"
num_examples(::Nothing) = 0
num_examples(xd::Union{XData,XYData}) = size(feature_matrix(xd))[1]
num_examples(bs::Union{XBatches,XYBatches}) = sum(b -> num_examples(b), bs)
num_examples(ds::Dataset) = num_examples(train(ds)) + num_examples(valid(ds)) + num_examples(test(ds))

"Count the number of example in a weighted dataset"
total_example_weight(::Nothing)::Float64 = 0.0
total_example_weight(d::PlainXData)::Float64 = num_examples(d)
total_example_weight(d::WXData)::Float64 = sum(weights(d))
total_example_weight(d::XYData)::Float64 = total_example_weight(feature_data(d))
total_example_weight(bs::Union{XBatches,XYBatches})::Float64 = sum(b -> total_example_weight(b), bs)
total_example_weight(ds::Dataset)::Float64 =
    total_example_weight(train(ds)) + total_example_weight(valid(ds)) + total_example_weight(test(ds))

"Number of batches in the dataset"
num_batches(bs::Union{XBatches,XYBatches}) = length(bs)
num_batches(ds::Dataset) = num_batches(train(ds)) + num_batches(valid(ds)) + num_batches(test(ds))

"Size of the largest batch"
max_batch_size(batch::Union{XData,XYData}) = num_examples(batch)
max_batch_size(batches::Union{XBatches,XYBatches}) = maximum(b -> max_batch_size(b), batches)

"Number of features in the data"
num_features(wxd::WXData) = size(feature_matrix(wxd))[2]
num_features(xd::Union{XData,XYData}) = size(feature_matrix(xd))[2]
num_features(bs::Union{XBatches,XYBatches}) = num_features(bs[1]) # assume constructor enforces consistency between batches
num_features(ds::Dataset) = num_features(train(ds)) # assume constructor enforces consistency between folds

"Number of labels in the data"
num_labels(xyd::XYData) = maximum(xyd.y)+1 # assume labels start at 0
num_labels(bs::XYBatches) = num_labels(bs[1])
num_labels(ds::LabeledDataset) = num_labels(train(ds))

"Get the unlabeled feature portion of the data"
feature_data(xyd::XYData)::XData = xyd.xd

"Get the data as a feature matrix"
feature_matrix(xd::XData) = xd.x
feature_matrix(wxd::WXData) = wxd.xd.x
feature_matrix(xyd::XYData) = feature_matrix(feature_data(xyd))

"Get the data as plain unweighted data"
unweighted_data(xd::PlainXData) = xd
unweighted_data(wxd::WXData) = wxd.xd

"Get the weight vector of the data"
weights(xd::WXData) = xd.w
weights(xd::PlainXData) = ones(num_examples(xd))

"Get the label vector of the data"
labels(xyd::XYData) = xyd.y

"What type should instance weights aggregate to?"
aggr_weight_type(x) = aggr_weight_type(typeof(x))
aggr_weight_type(::Type{<:XData}) = Int
aggr_weight_type(::Type{<:WXData{X,<:Integer}}) where X = Int # also covers Bool!
aggr_weight_type(::Type{<:WXData{X,<:AbstractFloat}}) where X = Float64
aggr_weight_type(::Type{B}) where {B<:Union{XBatches,XYBatches}} = aggr_weight_type(eltype(B))

"Get the type of the feature matrix elements"
feature_type(::Type{<:AbstractData{X}}) where X = X
feature_type(x::AbstractData) = feature_type(typeof(x))
feature_type(::AbstractVector{XD}) where {XD<:AbstractData} = feature_type(XD)

"Get the type of the label vector elements"
label_type(::XYData{X,Y}) where {X,Y} = Y
label_type(::Type{XYData{X,Y}}) where {X,Y} = Y
label_type(::AbstractVector{XYD}) where {XYD<:XYData} = label_type(XYD)

"Get the training fold of the data"
train(ds::Dataset) = ds.train

"Get the validation fold of the data"
valid(ds::Dataset) = ds.valid

"Get the testing fold of the data"
test(ds::Dataset) = ds.test

function assert_matching_folds(tr,v,te)
    @assert (num_features(tr) == num_features(te))
    @assert (isnothing(v) || num_features(tr)== num_features(v))
end

# basic data manipulations

function shuffle(xd::PlainXData, perm=nothing)
    issomething(perm) || (perm = randperm(num_examples(xd)))
    PlainXData(feature_matrix(xd)[perm,:])
end

function shuffle(wxd::WXData, perm=nothing)
    issomething(perm) || (perm = randperm(num_examples(wxd)))
    WXData(shuffle(wxd.xd, perm), weights(wxd)[perm])
end

function shuffle(xyd::XYData)
    perm = randperm(num_examples(xyd))
    XYData(shuffle(feature_data(xyd),perm),labels(xyd)[perm])
end

shuffle(::Nothing) = nothing
shuffle(d::D) where {D <: Union{XDataset,XYDataset}} =
    D(shuffle(train(d)), shuffle(valid(d)), shuffle(test(d)))
#no point in shuffling when already turned into mini batches, don't implement it

slice(xd::PlainXData, first, last) = PlainXData(feature_matrix(xd)[first:last, :])
slice(wxd::WXData, first, last) = WXData(feature_matrix(wxd)[first:last, :], weights(wxd)[first:last])
slice(xyd::XYData, first, last) = XYData(slice(feature_data(xyd),first,last), labels(xyd)[first:last])

"Get slices of data (rows and columns ids) as a view (no memory allocation)"
vslice(x::PlainXData, rows, cols) = PlainXData(@view feature_matrix(x)[rows, cols])
function vslice(xd::WXData, rows, cols)
    x = feature_matrix(xd)
    w = weights(xd)
    xv = @view x[rows, cols]
    wv = @view w[rows]
    return WXData(xv, wv)
end

function batch(xd::Union{XData,XYData}, batch_size::Integer)
    start_indices = 1:batch_size:num_examples(xd)
    # note that the last batch may be less than batch_size
    [slice(xd, start, min(num_examples(xd),start+batch_size-1)) for start in start_indices]
end
batch(::Nothing, ::Any) = nothing
batch(d::XDataset{X}, batch_size::Integer) where {X} =
    BatchedXDataset{X}(batch(train(d), batch_size),batch(valid(d), batch_size),batch(test(d), batch_size))
batch(d::XYDataset{X,Y}, batch_size::Integer) where {X,Y} =
    BatchedXYDataset{X,Y}(batch(train(d), batch_size), batch(valid(d), batch_size), batch(test(d), batch_size))

"Undo `batch` operation, combine batches into one"
unslice(xds::Vector{PlainXData{X, M}}) where {X, M} =
    PlainXData(vcat([feature_matrix(xd) for xd in xds]...))
unslice(wxds::AbstractVector{WXData{X, W, M}}) where {X, W, M} =
    WXData(vcat([feature_matrix(wxd) for wxd in wxds]...), vcat([weights(wxd) for wxd in wxds]...))
# unslice(xyds::AbstractVector{XYData}) = XYData(unslice(feature_data(xyds)), vcat([labels(xyds) for xyd in xyds]...)) # untest

unbatch(d::PlainXData{X, M}) where {X, M} = d
unbatch(d::XBatches{X}) where {X} = unslice(d)
unbatch(d::BatchedXDataset) where {X} =
    unbatch(train(d)), unbatch(valid(d)), unbatch(test(d)) # TODO: calling `XDataset` here triger signal (4) illegal hardware instruction, fix later

reweigh(wxdata::WXData, w::AbstractVector)::WXData =
    WXData(unweighted_data(wcdata), w)
reweigh(wxbatches::XBatches, ws::AbstractVector{<:AbstractVector})::XBatches =
    map((b,w) -> reweigh(b,w), wxbatches, ws)

threshold(xd) = threshold(xd, 0.05) # default threshold offset (used for MNIST)

function threshold(xd::XData, offset)
    x = feature_matrix(xd)
    means = Statistics.mean(x, dims=1)
    std = Statistics.std(x, dims=1)
    threshold = means .+ offset .* std
    th_x = BitArray(x .> threshold)
    (replace_features(xd,th_x), threshold)
end

threshold(xyd::XYData, offset) = replace_features(xyd,threshold(feature_matrix(xyd),offset)[1])

function threshold(ds::XDataset, offset)
    (tr,thr) = threshold(train(ds), offset)
    valid_m = map_something(v -> BitArray(feature_matrix(v) .> thr), valid(ds))
    test_m = BitArray(feature_matrix(test(ds)) .> thr)
    v = map_something(v -> replace_features(valid(ds),v), valid_m)
    te = replace_features(test(ds),test_m)
    XDataset{Bool}(tr,v,te)
end

function threshold(xyds::XYDataset{X,Y}, offset) where {X,Y}
    xds = threshold(XDataset(xyds), offset)
    v = map_something(v -> XYData(valid(xds),labels(v)), valid(xyds))
    XYDataset{Bool,Y}(XYData(xds.train,labels(train(xyds))), v, XYData(xds.test,labels(test(xyds))))
end

# basic inference

function fully_factorized_likelihood(batches::XBatches; pseudocount)
    counts = mapreduce(+, batches) do batch
        sum(feature_matrix(batch), dims=1)
    end
    smoothed_counts = counts .+ pseudocount/2.0
    log_estimates = log.(smoothed_counts ./ (num_examples(batches) + pseudocount))
    ll = sum(log_estimates .* counts)
    ll += sum(log.(1 .- exp.(log_estimates)) .* (num_examples(batches) .- counts))
    ll / num_examples(batches)
end

ll_per_example(ll::Float64, data) = ll / total_example_weight(data)
bits_per_pixel(ll::Float64, data) = -(ll_per_example(ll, data)  / num_features(data)) / log(2)

end # module


text:
# execute tests inside the cluster as a JuliaParBatch job
using JuliaRunClient
using Base.Test

const SUCCESS_INDICATOR = joinpath(ENV["JRUN_DIR"], "success.txt")

isfile(SUCCESS_INDICATOR) && rm(SUCCESS_INDICATOR)
ctx = Context()

# let's see if we can connect to the cluster
@test @result getSystemStatus(ctx)

# get status of self job
myself = self()
mystatus = @result getJobStatus(ctx, myself)
@test mystatus["created"]
@test !mystatus["succeeded"]

myscale = @result getJobScale(ctx, myself)
@test myscale[1] == myscale[2] == 0

initParallel()
sleep(10)
# scale up to 2 workers
@test @result setJobScale(ctx, myself, 2)
waitForWorkers(2)
@test length(workers()) == 2

@test @result setJobScale(ctx, myself, 0)
while true
    w = workers()
    (length(w) == 1) && (w[1] == 1) && break
    sleep(5)
end
touch(SUCCESS_INDICATOR)
chmod(SUCCESS_INDICATOR, 0o777)


text:
using Plots, ApproxFun

x=Fun()
u0=0.x

N=u->[u(-1.),u(1.),0.001u''+6*(1-x^2)*u'+u^2-1.]
u=newton(N,u0)
plot(u)

max_stars_repo_path: examples/Newton iteration for Nonlinear BVPs.jl


text:
module TestBackground
using Test
using BackgroundSubtraction
using LinearAlgebra
function synthetic_data(n = 128, m = 16)
    f(x) = sin(2π*x) + 1.1
    x = range(0, stop = 1, length = n)
    fx = f.(x)
    A = repeat(fx, 1, m)

    # add noise
    σ = 1e-2
    @. A += σ * randn()
    # add positive outliers (i.e. peaks)
    npeaks = n*m ÷ 4
    outind = rand(eachindex(A), npeaks)
    @. A[outind] += exp(randn())
    return A, x, fx
end

function synthetic_images(nx = 128, ny = 128, m = 16)
    f(x, y) = sin(2π*(x+y)) + 1.1
    x = range(0, stop = 1, length = nx)
    y = range(0, stop = 1, length = ny)
    fxy = f.(x, y')
    A = repeat(fxy, 1, 1, m)
    # add noise
    σ = 1e-2
    @. A += σ * randn()
    # add positive outliers (i.e. peaks)
    npeaks = nx*ny*m ÷ 4
    outind = rand(eachindex(A), npeaks)
    @. A[outind] += .1exp(randn())
    return A, x, y, fxy
end

@testset "background" begin
    n, m = 128, 16
    A, x, fx = synthetic_data(n, m)
    k = 2
    l = 1/2
    background = mcbl(A, k, x, l)
    @test size(background) == size(A)
    @test all(≥(0), background)
    tol = 5e-2
    @test maximum(abs, background[:, 1]-fx) < tol

    doplot = false
    if doplot
        using Plots
        plotly()
        plot(x, A[:, i], label = "data")
        plot!(x, background[:, i], label = "inferred background")
        plot!(x, fx, label = "background")
        gui()
    end

    # testing with only vector input: less powerful, and needs stronger l regularization
    i = 1
    background = mcbl(A[:,i], x, l, maxiter = 16)

    @test background isa Vector
    @test length(background) == length(A[:, i])
    tol = 5e-2
    @test maximum(abs, background[:, 1]-fx) < tol

    # testing smooth 2D projection
    lc = 1. # composition length scale
    c = randn(2, m) # 2 dimensional composition dimension
    a = A
    background = mcbl(a, x, l, c, lc)
    tol = 5e-2
    @test maximum(abs, background[:, 1]-fx) < tol

    doplot = false
    if doplot
        using Plots
        plotly()
        plot(x, A[:, i], label = "data")
        plot!(x, background[:, i], label = "inferred background")
        plot!(x, fx, label = "background")
        gui()
    end

    # testing 3D mcbl
    nx, ny, m = 128, 128, 16
    k = 1
    lx, ly = .2, .2
    A, x, y, fxy = synthetic_images(nx, ny, m)
    background = mcbl(A, k, x, lx, y, ly)
    tol = 5e-2
    @test maximum(abs, background[:, :, 1]-fxy) < tol

    # testing 3D smooth projection
    background = mcbl(A, x, lx, y, ly, c, lc)
    tol = 5e-2
    @test maximum(abs, background[:, :, 1]-fxy) < tol

    doplot = false
    if doplot
        using Plots
        plotly()
        surface(x, y, A[:, :, i], label = "data")
        gui()
        surface(x, y, background[:, :, i], label = "inferred background")
        gui()
        surface(x, y, fxy[:, :, i], label = "background")
        gui()
    end

    # testing 3D smooth projection
    # update default parameters

end

end # TestBackground


text:
abstract type AbstractIpState end

"""This structure saves the state of the material."""
struct IpStateSingle{T<:Real} <:AbstractIpState
    data::T
end

"""This structure saves the state of the material."""
struct IpStateArray{T<:Real, N} <:AbstractIpState
    data::AbstractArray{T, N}
end

"""
This function has to be implemented by the user for different types that they wish to implement.
This function gets the state of the material, If they exist in the Dictionary
for the given material/integration point in the given element,
it updates the data with the available data in stateDict.
If they don't exist, it just fills the state varibles with zeros.

    data = getState!(data, stateDict, elementNo, integrationPt)
"""
function getIpState!(data::T, stateDict::Dict{Tuple{Int64, Int64}, IpStateSingle},
    elementNo::Int64= 1, integrationPt::Int64=1) where {T <:Real}

    if (elementNo, integrationPt) ∈ keys(stateDict)
        data = stateDict[elementNo, integrationPt].data
    return data
    else
        zero = zeros(typeof(data), 1)
        data = zero[1]
        return data
    end
    return data
end

function getIpState!(data::AbstractArray{T,N}, stateDict::Dict{Tuple{Int64, Int64}, IpStateArray},
    elementNo::Int64= 1, integrationPt::Int64=1) where {T <:Real, N}

    if (elementNo, integrationPt) ∈ keys(stateDict)
        data .= stateDict[elementNo, integrationPt].data
    return data
    else
        zero = zeros(T, 1)
        fill!(data, zero[1])
        return data
    end
    return data
end

"""
This function updates the StateDict according to the passed data a specific element number and
an integration point within the given element.

    updateStateDict!(data, stateDict, elementNo, integrationPt)
"""
function updateIpStateDict!(data::T1, stateDict::Dict{Tuple{Int64, Int64}, T2},
    elementNo::Int64= 1, integrationPt::Int64=1) where {T1, T2}
    stateDict[elementNo, integrationPt] = T2(data)
    return nothing
end

"""
This function creates a Dictionary of a give Type to store the state of the material.

    stateDict = createIpStateDict(type)
"""

function createIpStateDict(T::Type)
    return Dict{Tuple{Int64, Int64}, T}()
end

function updateStateDict!(stateDict::Dict{Tuple{Int64, Int64}, T},
    stateDictBuffer::Dict{Tuple{Int64, Int64}, T}) where T
    merge!(stateDict, stateDictBuffer)
end

function updateIpStateDict4rmBuffer!(stateDict::Dict{Tuple{Int64, Int64}, T},
    stateDictBuffer::Dict{Tuple{Int64, Int64}, T}) where T
    updateStateDict!(stateDict, stateDictBuffer)
end


text:
using TensorBoardLogger #import the TensorBoardLogger package
using Logging #import Logging package
using TestImages
using Flux: Data
using Random
logger = TBLogger("imagelogs", tb_append) #create tensorboard logger

################log images example: mri################
mri = testimage("mri")

#using logger interface
with_logger(logger) do
    @info "image/mri/loggerinterface" mri
end
#using explicit function interface
log_image(logger, "image/mri/explicitinterface", mri, step = 0)

################log images example: MNIST data################
images = shuffle(Data.MNIST.images())[1:5]

#using logger interface
with_logger(logger) do
    @info "image/mnist/loggerinterface" images = TBImages(images, HW)
end
#using explicit function interface
log_images(logger, "image/mnist/explicitinterface", images, step = 0)

################log images example: random arrays################
noise = rand(16, 16, 3, 4) #Format is HWCN

#using logger interface
with_logger(logger) do
    @info "image/noise/loggerinterface" random = TBImage(noise, HWCN)
end
#using explicit function interface
log_image(logger, "image/noise/explicitinterface", noise, HWCN, step = 0)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# helper type for concise definitions
const PointOrGeometry{Dim,T} = Union{Point{Dim,T},Geometry{Dim,T}}

"""
    Collection(items)

A collection of `items` seen as a single entity.
"""
struct Collection{Dim,T,I<:PointOrGeometry{Dim,T}} <: Domain{Dim,T}
  items::Vector{I}
end

# -----------------
# DOMAIN INTERFACE
# -----------------

Base.getindex(c::Collection, ind::Int) = c.items[ind]

nelements(c::Collection) = length(c.items)

# ------------------------
# SPECIAL CASE: POINT SET
# ------------------------

const PointSet{Dim,T} = Collection{Dim,T,Point{Dim,T}}

"""
    PointSet(points)

A set of `points` (a.k.a. point cloud) seen as a single entity.

## Examples

All point sets below are the same and contain two points in R³:

```julia
julia> PointSet([Point(1,2,3), Point(4,5,6)])
julia> PointSet(Point(1,2,3), Point(4,5,6))
julia> PointSet([(1,2,3), (4,5,6)])
julia> PointSet((1,2,3), (4,5,6))
julia> PointSet([[1,2,3], [4,5,6]])
julia> PointSet([1,2,3], [4,5,6])
julia> PointSet([1 4; 2 5; 3 6])
```
"""
PointSet(points::AbstractVector{P}) where {P<:Point} = PointSet{embeddim(P),coordtype(P)}(points)
PointSet(points::Vararg{P}) where {P<:Point} = PointSet(collect(points))
PointSet(coords::AbstractVector{TP}) where {TP<:Tuple} = PointSet(Point.(coords))
PointSet(coords::Vararg{TP}) where {TP<:Tuple} = PointSet(collect(coords))
PointSet(coords::AbstractVector{V}) where {V<:AbstractVector} = PointSet(Point.(coords))
PointSet(coords::Vararg{V}) where {V<:AbstractVector} = PointSet(collect(coords))
PointSet(coords::AbstractMatrix) = PointSet(Point.(eachcol(coords)))

centroid(pset::PointSet, ind::Int) = pset[ind]

function Base.show(io::IO, pset::PointSet{Dim,T}) where {Dim,T}
  nelm = nelements(pset)
  print(io, "$nelm PointSet{$Dim,$T}")
end

# ---------------------------
# SPECIAL CASE: GEOMETRY SET
# ---------------------------

const GeometrySet{Dim,T,G<:Geometry{Dim,T}} = Collection{Dim,T,G}

"""
    GeometrySet(geometries)

A set of `geometries` seen as a single entity.
"""
GeometrySet(geometries::AbstractVector{G}) where {G<:Geometry} =
  GeometrySet{embeddim(G),coordtype(G),G}(geometries)

function Base.show(io::IO, gset::GeometrySet{Dim,T}) where {Dim,T}
  nelm = nelements(gset)
  print(io, "$nelm GeometrySet{$Dim,$T}")
end


text:
# ---
# title: 662. Maximum Width of Binary Tree
# id: problem662
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/maximum-width-of-binary-tree/description/>
# hidden: true
# ---
#
# Given a binary tree, write a function to get the maximum width of the given
# tree. The maximum width of a tree is the maximum width among all levels.
#
# The width of one level is defined as the length between the end-nodes (the
# leftmost and right most non-null nodes in the level, where the `null` nodes
# between the end-nodes are also counted into the length calculation.
#
# It is **guaranteed** that the answer will in the range of 32-bit signed
# integer.
#
# **Example 1:**
#
#
#
#     Input:
#
#                1
#              /   \
#             3     2
#            / \     \
#           5   3     9
#
#     Output: 4
#     Explanation: The maximum width existing in the third level with the length 4 (5,3,null,9).
#
#
# **Example 2:**
#
#
#
#     Input:
#
#               1
#              /
#             3
#            / \
#           5   3
#
#     Output: 2
#     Explanation: The maximum width existing in the third level with the length 2 (5,3).
#
#
# **Example 3:**
#
#
#
#     Input:
#
#               1
#              / \
#             3   2
#            /
#           5
#
#     Output: 2
#     Explanation: The maximum width existing in the second level with the length 2 (3,2).
#
#
# **Example 4:**
#
#
#
#     Input:
#
#               1
#              / \
#             3   2
#            /     \
#           5       9
#          /         \
#         6           7
#     Output: 8
#     Explanation: The maximum width existing in the fourth level with the length 8 (6,null,null,null,null,null,null,7).
#
#
#
#
# **Constraints:**
#
#   * The given binary tree will have between `1` and `3000` nodes.
#
#
## @lc code=start
using LeetCode

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


text:
export HeatMap ,saveplot , LinePlot, Axes

import Gaston: plot, heatmap , Axes , save

function saveplot end
saveoptions =  "font 'Consolas,10' size 1280,900 lw 1 background 'white'"
function (sim::Simulation)(::typeof(saveplot);term = "png",output = "./figure.png",saveopts = saveoptions)
    save(term = term , output = output,saveopts = saveopts)
end

function HeatMap end

heatmapax = Axes(auto = "fix", palette = :balance)
function (sim::Simulation)(::typeof(HeatMap), q::Symbol; apply::Function = real, axes = heatmapax, args...)
  x = range(sim.grid,x̂)
  y = range(sim.grid,ŷ)
  z  = getproperty(sim,q)
  heatmap(x,y,apply.(z)',axes; args...)
end

function (sim::Simulation)(::typeof(HeatMap), q::Symbol ,dir::Direction; apply::Function = real , axes = heatmapax ,  args...)
  x = range(sim.grid,x̂)
  y = range(sim.grid,ŷ)
  z  = sim(q,dir)
  heatmap(x,y,apply.(z)',axes; args...)
end

function (sim::Simulation)(::typeof(HeatMap), z::AbstractArray ; apply::Function = real,axes = heatmapax, args...)
  x = range(sim.grid,x̂)
  y = range(sim.grid,ŷ)
  heatmap(x,y,apply.(z)',axes; args...)
end

function LinePlot end

function (sim::Simulation)(::typeof(LinePlot), q::Symbol; apply::Function = real, axes = Axes(auto = "fix"), args...)
  x = range(sim.grid,x̂)
  z  = getproperty(sim,q)
  plot(x,apply.(z),axes,w = :l ; args...)
end

function (sim::Simulation)(::typeof(LinePlot), z::AbstractArray ; apply::Function = real,axes = Axes(auto = "fix"), args...)
  x = range(sim.grid,x̂)
  plot(x,apply.(z),axes,w = :l ; args...)
end

function (sim::Simulation)(::typeof(LinePlot), q::Symbol ,dir::Direction; apply::Function = real , axes = Axes(auto = "fix") ,  args...)
  x = range(sim.grid,x̂)
  z  = sim(q,dir)
  plot(x,apply.(z),axes,w = :l ; args...)
end


text:
## ----------------------------  Linear tuition function
# In gpa, parental pct, quality
# Optional additive term by year in college. Purpose is to penalize late graduation.
# Can be calibrated by multiplying entire time profile by a number.

# When constructing these switches, set `gpaGradient` etc to data values.
Base.@kwdef mutable struct TuitionLinearSwitches <: AbstractTuitionSwitches
    nColleges :: CollInt
    qualBaseV :: Vector{Double}
    calQualBase :: Bool = true
    gpaGradient :: Double = 0.0
    calGpaGradient :: Bool = true
    parentalGradient :: Double = 0.0
    calParentalGradient :: Bool = true
    # Add-on pattern by year in college
    byYearV :: Vector{Double} = fill(0.0, 6)
    byYearFactor :: Double = 0.0
    calByYear :: Bool = false
    # For experiments: specify which colleges are free
    # Cannot be used with calibration (then qualBaseV no longer matters for free colleges)
    freeIdxV :: Vector{CollInt} = Vector{CollInt}()
end

"""
	$(SIGNATURES)

Linear Tuition function. A base value for each quality plus linear terms in HS GPA percentile and parental income percentile. No time variation.
Coefficients may be calibrated or fixed. Stored in model dollars.
"""
mutable struct TuitionFunctionLinear{T1 <: AbstractFloat} <: AbstractTuitionFunction
    objId :: ObjectId
    switches :: TuitionLinearSwitches
    pvec :: ParamVector

    # Base tuition by quality
    qualityBaseV :: Vector{T1}
    # Gradient w.r.to HS GPA percentile
    gpaGradient :: T1
    # Gradient w.r.to parental percentile
    parentalGradient :: T1
    byYearFactor :: T1
end

function StructLH.describe(switches :: TuitionLinearSwitches)
    gpaStr = calibrated_string(switches.calGpaGradient;
        fixedValue = round.(switches.gpaGradient; digits = 2));
    ypStr = calibrated_string(switches.calParentalGradient;
        fixedValue = round.(switches.parentalGradient; digits = 2));
    qualityStr = calibrated_string(switches.calQualBase;
        fixedValue = round.(switches.qualBaseV; digits = 2));
    return [
        "Tuition "  "linear in quality, gpa, parental pct.";
        "GPA gradient"  gpaStr;
        "Parental gradient "  ypStr;
        "Quality dummies"  qualityStr
    ]
end

## -----   Convenience

cal_gpa_gradient(s :: TuitionLinearSwitches, doCal :: Bool) =
    s.calGpaGradient = doCal;
cal_parental_gradient(s :: TuitionLinearSwitches, doCal :: Bool) =
    s.calParentalGradient = doCal;
cal_qual_base(s :: TuitionLinearSwitches, doCal :: Bool) =
    s.calQualBase = doCal;
cal_year_add_on!(s :: TuitionLinearSwitches, doCal :: Bool) =
    s.calByYear = doCal;

quality_base(tf :: TuitionFunctionLinear, qual) = tf.qualityBaseV[qual];
gpa_gradient(tf :: TuitionFunctionLinear) = tf.gpaGradient;
parental_gradient(tf :: TuitionFunctionLinear) = tf.parentalGradient;
year_add_on(tf :: TuitionFunctionLinear, t :: Integer) =
    tf.byYearFactor * tf.switches.byYearV[t];

isfree(tf :: TuitionFunctionLinear, qual :: Integer) =
    isfree(tf.switches, qual);

isfree(switches :: TuitionLinearSwitches, qual :: Integer) =
    qual ∈ switches.freeIdxV;

function make_free!(switches :: TuitionLinearSwitches, iCollege :: Integer)
    !isfree(switches, iCollege)  &&  push!(switches.freeIdxV, iCollege);
end

## -----  Initialize

test_tuition_linear_switches(nc :: Integer; freeIdxV = Vector{CollInt}()) =
    TuitionLinearSwitches(
        nColleges = nc,
        gpaGradient = 1.0,
        parentalGradient = 2.0,
        qualBaseV = collect(LinRange(3.0, 5.0, nc)),
        freeIdxV = freeIdxV);

# Initialize tuition function, taking coefficients from data regression.
function init_tuition_fct(objId :: ObjectId,  switches :: TuitionLinearSwitches)
    # rt = load_tuition_regr(ds, modelUnits = true);
    pQualBase = init_qual_base(switches);
    pGpaGrad = init_gpa_gradient(switches);
    pParentalGrad = init_parental_gradient(switches);
    pByYear = init_by_year(switches);
    pvec = ParamVector(objId, [pQualBase, pGpaGrad, pParentalGrad, pByYear]);

    tf = TuitionFunctionLinear(objId, switches, pvec,
        ModelParams.value(pQualBase), ModelParams.value(pGpaGrad),
        ModelParams.value(pParentalGrad), ModelParams.value(pByYear));
    @assert validate_tf(tf)
    return tf
end

function init_qual_base(switches :: TuitionLinearSwitches)
    # nameV = [regressor_name(:quality, ic)  for ic = 2 : n_colleges(switches)];
    # dataV, _ = get_coeff_se_multiple(rt, vcat(intercept_name(), nameV));

    dataV = switches.qualBaseV;
    p = Param(:qualityBaseV,  ldescription(:tuitionBaseV), lsymbol(:tuitionBaseV),
        dataV, dataV, dataV .- 5.0, dataV .+ 5.0, switches.calQualBase);
    return p
end

# There is no good way of using the regression, which is for GPA quartiles.
function init_gpa_gradient(switches :: TuitionLinearSwitches)
    # Gradient = last gpa coefficient - 0 (first coefficient; omitted)
    # nGpa = find_n_regressors(rt, :afqt);
    # grad = get_coefficient(rt, regressor_name(:gpa, nGpa));
    grad = switches.gpaGradient;
    # Units are: model dollars as gpa percentile goes from 0 to 100
    lb = grad - dollars_data_to_model(2000.0, :perYear);
    ub = grad + dollars_data_to_model(2000.0, :perYear);
    p = Param(:gpaGradient,  ldescription(:tuitionGpaGrad), lsymbol(:tuitionGpaGrad),
        grad, grad, lb, ub, switches.calGpaGradient);
    return p
end

function init_parental_gradient(switches :: TuitionLinearSwitches)
    # nParental = find_n_regressors(rt, :parental);
    # grad = get_coefficient(rt, regressor_name(:parental, nParental));
    grad = switches.parentalGradient;
    # Units are: model dollars as gpa percentile goes from 0 to 100
    lb = grad - dollars_data_to_model(2000.0, :perYear);
    ub = grad + dollars_data_to_model(2000.0, :perYear);
    p = Param(:parentalGradient,
         ldescription(:tuitionYpGrad), lsymbol(:tuitionYpGrad),
        grad, grad, lb, ub, switches.calParentalGradient);
    return p
end

function init_by_year(switches :: TuitionLinearSwitches)
    byFactor = switches.byYearFactor;
    ub = dollars_data_to_model(5000.0, :perYear);
    p = Param(:byYearFactor, "Tuition scale for add on by year", "byYear",
        byFactor, byFactor, 0.0, ub, switches.calByYear);
    return p
end

##  ----  Compute tuition

"""
	$(SIGNATURES)

Tuition for given characteristics. Inputs can be scalar or array as long as dimensions conform. `qual` is expected to be `Integer`.
"""
function tuition(tf :: TuitionFunctionLinear, qual, gpa, parental, t :: Integer;
    modelUnits :: Bool = true)
    if dbgHigh
        @assert all_at_most(gpa, 1.0)
        @assert all_at_least(gpa, 0.0)
        @assert all_at_most(parental, 1.0)
        @assert all_at_least(parental, 0.0)
    end

    tuitionV = (
        quality_base(tf, qual) .+
        gpa_gradient(tf) .* gpa .+
        parental_gradient(tf) .* parental .+
        year_add_on(tf, t))  .* free_factor(tf, qual);
    # tuitionV = zero_out_free_colleges(tf, tuitionV, qual);
    if !modelUnits
        tuitionV = dollars_model_to_data(tuitionV, :perYear);
    end
    return tuitionV
end

# -----------


text:
#= extract_object_from_radar.jl
extract the object_list from raw-data of radar =#

using DataFrames, XLSX;
include("visualize_radar_points.jl")
using .visualize

prj_path = dirname(dirname(@__FILE__));
xls_file_path = prj_path * raw"\data\radar_process\radar_points.xlsx";

#=
for id_sample in 1:38
    sh_name = "$id_sample.pcd"
    df=DataFrame(XLSX.readtable(
        xls_file_path,
        sh_name));
    println(df)
end =#

data, col_names = XLSX.readtable(xls_file_path, "1.pcd", header=true);

df = DataFrame(data, col_names);

x = convert(Array{Float64,1}, df.x);
y = convert(Array{Float64,1}, df.y);
vx_comp = convert(Array{Float64,1}, df.vx_comp);
vy_comp = convert(Array{Float64,1}, df.vy_comp);
v_comp = @. abs(vx_comp + im * vy_comp)

# set the v_comp-threshold(10km/s)
v_comp_threshold = 10 / 3.6;
id_v_comp = findall(x -> x >= v_comp_threshold, v_comp);

println(names(df))