text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using InteractiveDynamics, Random
using DynamicalBilliards
import GLMakie
using GLMakie: to_color, RGBf0, RGBAf0
using DynamicalSystems

fig = GLMakie.Figure(resolution = (figx*100*0.8, figy*100), fontsize = 32)
axhe = GLMakie.Axis3(fig[1,1])
axbi = GLMakie.Axis3(fig[1,2])
display(fig)

# %% Natural measure of henon map
ε = 0.001
x, y, p = begin
    ds = Systems.henon()
    X = trajectory(ds, 10^7; Ttr = 100)
    p, b = binhist(X, ε)
    x = [a[1] for a in b]
    y = [a[2] for a in b]
    x, y, p ./ maximum(p)
end

using GLMakie: Point3f0, Vec3f0, Rect3D

GLMakie.meshscatter!(axhe, vec(Point3f0.(x, y, 0.0)),
    markersize=Vec3f0.(2ε, 2ε, p), marker=Rect3D(Vec3f0(0), Vec3f0(1)),
    limits=Rect3D(Vec3f0(0), Vec3f0(1)),
    color = clamp.(p, 0, 0.25), alpha = 0.5
)
# GLMakie.ylims!(axhe, 1.5 .* extrema(y)...)
# GLMakie.xlims!(axhe, extrema(x)...)
GLMakie.zlims!(axhe, 0, 0.2)
# GLMakie.xticks!(fig; )
axhe.zlabel = "ρ"

axhe.azimuth = 7.175530633326986
axhe.elevation = 0.53
axhe.zticklabelsvisible = false
axhe.zlabeloffset = 20
axhe.xticklabelsize = 26
axhe.yticklabelsize = 26
axhe.title = "Hénon map"

# %% same but for mushroom billiard
w = 0.2
bd = billiard_mushroom(1, w, 1, 0)

pc = Particle(-0.01, 0.2, sqrt(3)/2)
pr = Particle(0.0, 1.2, 0.0)
pr2 = Particle(0.0, 1.9, 0.0)

colors = [to_color(COLORS[i]) for i in (1,2,4)]
cmaps = [[InteractiveDynamics.darken_color(c, 2), c, InteractiveDynamics.lighten_color(c, 2)] for c in colors]
particles = [pc, pr, pr2]

ε = 0.02

for i in 1:3
    x, y, p = begin
        bmap, = boundarymap(particles[i], bd, 10^7)
        X = Dataset(bmap)
        p, b = ChaosTools.binhist(X, ε)
        x = [a[1]/2 for a in b]
        y = [a[2] for a in b]

        x, y, p ./ maximum(p)
    end

    GLMakie.meshscatter!(axbi, vec(Point3f0.(x, y, 0.0)),
    markersize=Vec3f0.(ε, ε, p), marker=Rect3D(Vec3f0(0), Vec3f0(1)),
    limits=Rect3D(Vec3f0(0), Vec3f0(1)),
    color = clamp.(p, 0.25, 1), colormap = cmaps[i])
end
GLMakie.zlims!(axbi, 0, 1)

axbi.zlabel = "ρ"
axbi.zticklabelsvisible = false
axbi.xlabel = "ξ"
axbi.ylabel = "sin(θ)"
axbi.azimuth = 10.225530633326988
axbi.elevation = 0.53
axbi.zlabeloffset = 20
axbi.xticklabelsize = 26
axbi.yticklabelsize = 26
axbi.title = "mushroom billiard"

# %% Save
GLMakie.save(plotsdir("8", "natural_densities.png"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using Distributions
using FFTW
using Plots

# Create a simple signal with two frequencies
dt = 0.001
t = 0:dt:1
n = length(t)

f = sin.(50*2pi*t) + sin.(120*2pi*t) # Sum of 2 frequencies
f += 2.5*rand(Normal(), n) # Add some noise

# Compute the Fast Fourier Transform FFT
fhat = fft(f) # Compute the fast Fourier transform
PSD = abs2.(fhat)/n # Power spectrum (power per freq)
freq = 1/(dt*n)*(0:n) # Create x-axis of frequencies in Hz

# Use the PSD to filter out noise
indices = PSD .> 100 # Find all freqs with large power
PSDclean = PSD.*indices # Zero out all others
fhat = indices.*fhat # Zero out small Fourier coeffs. in Y
ffilt = real.(ifft(fhat)) # Inverse FFT for filtered time signal

# PLOTS
p1 = plot(t, f)
plot!(t, ffilt)

L = 1:Int(floor(n/2)) # Only plot the first half of freqs
p2 = plot(freq[L], PSD[L])
plot!(freq[L], PSDclean[L])

ext: jl


text:
abstract type SmoothingKernel end

struct IdentitySmoother <: SmoothingKernel
    S::Float64
    κ₁::Float64
    κ₂::Float64
    κ₃::Float64
end

struct TruncatedSmoother <: SmoothingKernel
    ξ::Integer
    S::Float64
    smoother::Function
    κ₁::Float64
    κ₂::Float64
    κ₃::Float64
end

struct BartlettSmoother <: SmoothingKernel
    ξ::Integer
    S::Float64
    smoother::Function
    κ₁::Float64
    κ₂::Float64
    κ₃::Float64
end

IdentitySmoother() = IdentitySmoother(1.0, 1.0, 1.0, 1.0)

function TruncatedSmoother(ξ::Integer)
    function smoother(G::Array{T, 2}) where T
        N, M = size(G)
        nG   = zeros(T, N, M)
        for m=1:M
            for t=1:N
                low = max((t-N), -ξ)
                high = min(t-1, ξ)
                for s = low:high
                    @inbounds nG[t, m] += G[t-s, m]
                end
            end
        end
        return(nG/(2.0*ξ+1.0))
    end
    TruncatedSmoother(ξ, (2.0*ξ+1.0)/2.0, smoother, 2.0, 2.0, 1.0)
end

function BartlettSmoother(ξ::Integer)
    function smoother(G::Array{T, 2}) where T
        N, M = size(G)
        nG   = zeros(T, N, M)
        St   = (2.0*ξ+1.0)/2.0
        for m=1:M
            for t=1:N
                low = max((t-N), -ξ)
                high = min(t-1, ξ)
                for s = low:high
                    κ = 1.0-s/St
                    @inbounds nG[t, m] += κ*G[t-s, m]
                end
            end
        end
        return(nG/(2*ξ+1))
    end
    BartlettSmoother(ξ, (2.0*ξ+1.0)/2.0, smoother, 1.0, 2.0/3.0, 0.5)
end

## Used to scale objective function
bw(k::SmoothingKernel) = k.S
κ₁(k::SmoothingKernel) = k.κ₁
κ₂(k::SmoothingKernel) = k.κ₂

smooth(g::Array{T, 2}, k::IdentitySmoother) where T = g
smooth(g::Array{T, 2}, k::SmoothingKernel) where T = k.smoother(g)


text:
using ReversibleJumpMCMC
using Test

@testset "ReversibleJumpMCMC.jl" begin
    # Write your tests here.
    @test 1==1
end


text:
using ITensors
using ITensorNetworkAD
using ITensorsVisualization

using ITensorNetworkAD.ITensorNetworks:
  inds_network, project_boundary, contract_approx, insert_projectors
using ITensorNetworkAD.ITensorNetworks: Models

model = Models.Model"ising"()
βc = Models.critical_point(model)
β = 1.001 * βc
@show β / βc

N = (3, 4)
ndims = length(N)

space = [QN(0) => 2]
#space = 2
tn_inds = inds_network(N...; linkdims=space, addtags="S=1/2")
A = Models.local_boltzmann_weight("ising", Val(ndims); β)

tn = itensor.((A,), tn_inds)

# Project periodic boundary indices
# onto the 1 state
state = 1
tn = project_boundary(tn, state)

_cutoff = 1e-15
_maxdim = 100

# Contract in every directions
boundary_mps = contract_approx(tn; maxdim=_maxdim, cutoff=_cutoff)

#
# Insert projectors horizontally (to measure e.g. properties
# in a row of the network)
#

row = 2
center = (row, :)
tn_projected = insert_projectors(tn, boundary_mps; center=center)
tn_split, Pl, Pr = tn_projected

Pl_flat = reduce(vcat, Pl)
Pr_flat = reduce(vcat, Pr)
tn_projected_flat = mapreduce(vec, vcat, (tn_split, Pl_flat, Pr_flat))

@show noncommoninds(tn_projected_flat...)
@visualize *(tn_projected_flat...) contract = false pause = true

@disable_warn_order begin
  @show contract(tn_projected_flat)[] / contract(vec(tn))[]
end

#
# Insert projectors vertically (to measure e.g. properties
# in a column of the network)
#

column = 2
center = (:, column)
tn_projected = insert_projectors(tn, boundary_mps; center=center)
tn_split, Pl, Pr = tn_projected

Pl_flat = reduce(vcat, Pl)
Pr_flat = reduce(vcat, Pr)
tn_projected_flat = mapreduce(vec, vcat, (tn_split, Pl_flat, Pr_flat))

@show noncommoninds(tn_projected_flat...)
@visualize *(tn_projected_flat...) contract = false

@disable_warn_order begin
  @show contract(tn_projected_flat)[] / contract(vec(tn))[]
end


text:
module IonBase

using Comonicon
using Comonicon.Tools
using PkgTemplates
using LocalRegistry
using RegistryTools
using GitHub
using StringDistances
using OrderedCollections
using Pkg
using Pkg.TOML

using RegistryTools: gitcmd
using PkgTemplates: @plugin, @with_kw_noshow
using Comonicon.Parse: default_name

Comonicon.set_brief_length!(120)

# GITHUB_TOKEN is used in github actions
# GITHUB_AUTH is suggested by GitHub.jl
const ENV_TOKEN_NAMES = ["GITHUB_TOKEN", "GITHUB_AUTH"]

module PATH

using ..IonBase

function default_ion_dir()
    if haskey(ENV, "DOT_ION_PATH")
        return ENV["DOT_ION_PATH"]
    else
        return joinpath(homedir(), ".ion")
    end
end

dot_ion(xs...) = joinpath(default_ion_dir(), xs...)

function templates(xs...)
    joinpath(dirname(dirname(pathof(IonBase))), "templates", xs...)
end

end

# Julia Pkg commands
include("project.jl")
include("registry.jl")

# extra Ion commands
include("create.jl")
include(joinpath("templates", "default.jl"))
include(joinpath("templates", "command.jl"))
include(joinpath("templates", "from_file.jl"))

include("install.jl")
include("clone.jl")
include("doc.jl")
include("release.jl")
include("search.jl")
include("utils.jl")

include("precompile.jl")
_precompile_()

end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using DocOpt

@inline function cndf2(in)
    out = 0.5 .+ 0.5 .* erf(0.707106781 .* in)
    return out
end

function blackscholes(sptprice, strike, rate, volatility, time)
    logterm = log10(sptprice ./ strike)
    powterm = .5 .* volatility .* volatility
    den = volatility .* sqrt(time)
    d1 = (((rate .+ powterm) .* time) .+ logterm) ./ den
    d2 = d1 .- den
    NofXd1 = cndf2(d1)
    NofXd2 = cndf2(d2)
    futureValue = strike .* exp(- rate .* time)
    c1 = futureValue .* NofXd2
    call = sptprice .* NofXd1 .- c1
    put  = call .- futureValue .+ sptprice

    # Please don't delete these.
    # We use them when we are working on some new feature and want easier blackscholes debugging.

    #    logterm = log10(sptprice ./ strike)
    #    powterm = .5 .* volatility .* volatility
    #    den = volatility .* sqrt(time)
    #    rp = rate .+ powterm
    #    rpt = rp .* time
    #    rptl = rpt .+ logterm
    #    d1 = rptl ./ den
    #    d2 = d1 .- den
    #    NofXd1 = cndf2(d1)
    #    NofXd2 = cndf2(d2)
    #    futureValue = strike .* exp(- rate .* time)
    #    c1 = futureValue .* NofXd2
    #    c2 = sptprice .* NofXd1
    #    call = c2 .- c1
    #    put  = call .- futureValue .+ sptprice;
    #    println("sptprice ",sptprice)
    #    println("strike ",strike)
    #    println("logterm ",logterm)
    #    println("volatility ",volatility)
    #    println("powterm ",powterm)
    #    println("time ",time)
    #    println("den ",den)
    #    println("rate ",rate)
    #    println("rp ",rp)
    #    println("rpt ",rpt)
    #    println("rptl ",rptl)
    #    println("d1 ",d1)
    #    println("d2 ",d2)
    #    println("NofXd1 ",NofXd1)
    #    println("NofXd2 ",NofXd2)
    #    println("futureValue ",futureValue)
    #    println("c1 ",c1)
    #    println("c2 ",c2)
    #    println("call ",call)
    #    println("put ",put)
    #    (powterm, den, rate, rp, rpt, rptl, d1, d2, NofXd1, NofXd2, futureValue, c1, c2, call, put)
    #    powterm
end

function run(iterations)
    sptprice   = Float64[ 42.0 for i = 1:iterations ]
    initStrike = Float64[ 40.0 + (i / iterations) for i = 1:iterations ]
    rate       = Float64[ 0.5 for i = 1:iterations ]
    volatility = Float64[ 0.2 for i = 1:iterations ]
    time       = Float64[ 0.5 for i = 1:iterations ]

    tic()
    put = blackscholes(sptprice, initStrike, rate, volatility, time)
    t = toq()
    println("checksum: ", sum(put))
    return t
end

function main()
    doc = """black-scholes.jl

Black-Scholes option pricing model.

Usage:
  black-scholes.jl -h | --help
  black-scholes.jl [--iterations=<iterations>]

Options:
  -h --help                  Show this screen.
  --iterations=<iterations>  Specify a number of iterations [default: 10000000].
"""
    arguments = docopt(doc)

    iterations = parse(Int, arguments["--iterations"])

    srand(0)

    println("iterations = ",iterations)

    tic()
    blackscholes(Float64[], Float64[], Float64[], Float64[], Float64[])
    println("SELFPRIMED ", toq())

    time = run(iterations)
    println("rate = ", iterations / time, " opts/sec")
    println("SELFTIMED ", time)

end

main()


text:
global breakFromLoop = false

"""
    stopVideoProducer()
Stop the running video producer and close the device and channel.
"""
function stopVideoProducer()
    global breakFromLoop = true
end

"""
    vidproducer(c::Channel, reader::T, devicename::String = "/dev/video0",  N::Int = 100) where T <: AbstractEncodings

"""
function videoproducer(c::Channel, decoder::T; devicename::String = "/dev/video0",  N::Int = 100, iomethod::IOMethods = Video4Linux.IO_METHOD_READ) where T <: AbstractEncodings
    global breakFromLoop = false
    ##
    set_io_method(iomethod)

    ## open device
    fid = open_device(devicename)

    ## init_device(fd, force_format);
    init_device(fid) != 0 && error("Erro in device init")

    ## start_capturing(fd);
    start_capturing(fid) != 0 && error("Erro in start_capturing")

    imy = (decoder.depth > 1) ? zeros(decoder.datatype, decoder.height, decoder.width, decoder.depth) :
                               zeros(decoder.datatype, decoder.height, decoder.width)

    NN = N > 0 ? N : Inf
    i = 0
    while i < NN
        i += 1
        mainloop( fid, 1 ) != 0 && error("Erro in mainloop")
        ## copy_buffer_bytes, copy the image buffer bytes to uint8 vector, the lenght will depend on the pixel format
        decoder(imy)
        put!(c, imy)
        if breakFromLoop
            break
        end
    end

    ## stop_capturing(fd);
    stop_capturing(fid) != 0 && error("Erro in stop_capturing")

    ## uninit_device();
    uninit_device(fid) != 0 && error("Erro in uninit_device")

    ## close device
    close_device(fid) != 0 && error("Erro in close_device")
end

##
#vidchan = Channel((c::Channel) -> videoproducer(c, yonly))


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

module SEIRmodel
using OrdinaryDiffEq

#Susceptible-exposed-infected-recovered model function
function seir_ode(dY,Y,p,t)
    #Infected per-Capita Rate
    β = p[1]
    #Incubation Rate
    σ = p[2]
    #Recover per-capita rate
    γ = p[3]
    #Death Rate
    μ = p[4]

    #Susceptible Individual
    S = Y[1]
    #Exposed Individual
    E = Y[2]
    #Infected Individual
    I = Y[3]
    #Recovered Individual
    #R = Y[4]

    dY[1] = μ-β*S*I-μ*S
    dY[2] = β*S*I-(σ+μ)*E
    dY[3] = σ*E - (γ+μ)*I
end

function main()
    #Pram (Infected Rate, Incubation Rate, Recover Rate, Death Rate)
    pram=[520/365,1/60,1/30,774835/(65640000*365)]
    #Initialize Param(Susceptible Individuals, Exposed Individuals, Infected Individuals)
    init=[0.8,0.1,0.1]
    tspan=(0.0,365.0)

    seir_prob = ODEProblem(seir_ode,init,tspan,pram)

    sol=solve(seir_prob, alg=Tsit5());
    return sol
end

# using Plots

# va = VectorOfArray(sol.u)
# y = convert(Array,va)
# R = ones(size(sol.t))' - sum(y,dims=1);

# plot(sol.t,[y',R'],xlabel="Time",ylabel="Proportion")

end


text:
export Do

function Do(m::Model, xs::Symbol...)
    theModule = getmodule(m)
    newargs = (arguments(m) ∪ xs)

    m_init = Model(theModule, newargs, NamedTuple(), NamedTuple(), nothing)
    m = foldl(setdiff(parameters(m), xs); init=m_init) do m0,v
        merge(m0, Model(theModule, findStatement(m, v)))
    end
end

Do(m::Model; kwargs...) = Do(m,keys(kwargs)...)(;kwargs...)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    Tests of alignments.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Alignments" begin
    # Left
    # ==========================================================================
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│ 1      │ false  │ 1.0    │ 1      │
│ 2      │ true   │ 2.0    │ 2      │
│ 3      │ false  │ 3.0    │ 3      │
│ 4      │ true   │ 4.0    │ 4      │
│ 5      │ false  │ 5.0    │ 5      │
│ 6      │ true   │ 6.0    │ 6      │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data; alignment = :l)
    @test result == expected

    # Center
    # ==========================================================================
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│   1    │ false  │  1.0   │   1    │
│   2    │  true  │  2.0   │   2    │
│   3    │ false  │  3.0   │   3    │
│   4    │  true  │  4.0   │   4    │
│   5    │ false  │  5.0   │   5    │
│   6    │  true  │  6.0   │   6    │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data; alignment = :c)
    @test result == expected

    # Per column configuration
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│ 1      │  false │  1.0   │      1 │
│ 2      │   true │  2.0   │      2 │
│ 3      │  false │  3.0   │      3 │
│ 4      │   true │  4.0   │      4 │
│ 5      │  false │  5.0   │      5 │
│ 6      │   true │  6.0   │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data; alignment = [:l,:r,:c,:r])
    @test result == expected

    # Cell override
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│ 1      │  false │  1.0   │ 1      │
│ 2      │   true │  2.0   │      2 │
│      3 │ false  │  3.0   │   3    │
│ 4      │   true │  4.0   │   4    │
│ 5      │  false │  5.0   │      5 │
│ 6      │   true │  6.0   │ 6      │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data;
                          alignment = [:l,:r,:c,:r],
                          cell_alignment = Dict( (3,1) => :r,
                                                 (3,2) => :l,
                                                 (1,4) => :l,
                                                 (3,4) => :c,
                                                 (4,4) => :c,
                                                 (6,4) => :l ))
    @test result == expected

    # Headers
    # ==========================================================================

    header = ["A" "B" "C" "D"
              "a" "b" "c" "d"]

    expected = """
┌───────────┬───────────┬───────────┬───────────┐
│ A         │     B     │         C │         D │
│ a         │     b     │         c │         d │
├───────────┼───────────┼───────────┼───────────┤
│ 1         │     false │    1.0    │ 1         │
│ 2         │      true │    2.0    │         2 │
│         3 │ false     │    3.0    │     3     │
│ 4         │      true │    4.0    │     4     │
│ 5         │     false │    5.0    │         5 │
│ 6         │      true │    6.0    │ 6         │
└───────────┴───────────┴───────────┴───────────┘
"""

    result = pretty_table(String, data, header;
                          alignment = [:l,:r,:c,:r],
                          cell_alignment = Dict( (3,1) => :r,
                                                 (3,2) => :l,
                                                 (1,4) => :l,
                                                 (3,4) => :c,
                                                 (4,4) => :c,
                                                 (6,4) => :l ),
                          columns_width = 9,
                          header_alignment = [:l,:c,:r,:r])
    @test result == expected

    expected = """
┌───────────┬───────────┬───────────┬───────────┐
│         A │     B     │         C │         D │
│ a         │ b         │         c │     d     │
├───────────┼───────────┼───────────┼───────────┤
│ 1         │     false │    1.0    │ 1         │
│ 2         │      true │    2.0    │         2 │
│         3 │ false     │    3.0    │     3     │
│ 4         │      true │    4.0    │     4     │
│ 5         │     false │    5.0    │         5 │
│ 6         │      true │    6.0    │ 6         │
└───────────┴───────────┴───────────┴───────────┘
"""

    result = pretty_table(String, data, header;
                          alignment = [:l,:r,:c,:r],
                          cell_alignment = Dict( (3,1) => :r,
                                                 (3,2) => :l,
                                                 (1,4) => :l,
                                                 (3,4) => :c,
                                                 (4,4) => :c,
                                                 (6,4) => :l ),
                          columns_width = 9,
                          header_alignment = [:l,:c,:r,:r],
                          header_cell_alignment = Dict( (1,1) => :r,
                                                        (2,2) => :l,
                                                        (2,4) => :c))
    @test result == expected
end


text:
module TestTransducers

using Test
using Transducers
using UnionArrays

using ..Utils: compile_enabled

function test()
    xs = UnionVector(Any[1, 2.0])
    @test foldl(+, Map(x -> 2x), xs) == 6

    ys = eduction(Map(x -> 2x), xs)
    if compile_enabled()
        @test (@inferred foldl(+, ys, init=0.0)) == 6
    end
end

end  # module


text:
# ------------ Dynamic addressing hint ------------ #

# Provides a warning to users when using dynamic addressing.

struct DynamicAddressingHint <: AddressingHint
    dynamic_addresses::Vector{Variable}
    ir_map::Dict{Variable, IRTools.Statement}
end

@inline Base.isempty(dah::DynamicAddressingHint) = isempty(dah.dynamic_addresses)

function Base.display(dah::DynamicAddressingHint)
    if isempty(dah)
        println("\u001b[32m ✓ (DynamicAddressHint): Compiler detected no dynamic addresses in your model code.\u001b[0m")
    else
        println("\u001b[33m ? (DynamicAddressHint): Compiler detected the following dynamic addresses in your model code.\u001b[0m")
        println(" ________________________\n")
        display(dah.ir_map)
        println(" ________________________\n")
        println("\u001b[33m    This will prevent the compiler from inferring a trace type for your model program, as well as statically checking the support.\u001b[0m\n")
        println("\u001b[32m ⇒ (Recommendation): transform control flow constructs into \u001b[3m\u001b[34;1mcombinators\u001b[0m\n\n      \u001b[34;1mhttps://www.gen.dev/dev/ref/combinators/\u001b[0m\n\n\u001b[32m    Then, provide static (e.g. \e[1mSymbol\u001b[0m\u001b[32m) addresses to combinator calls.\u001b[0m")
    end
end

function detect_dynamic_addresses(ir::IR)
    dynamic_addresses = Variable[]
    ir_map = Dict{Variable, IRTools.Statement}()
    for (v, st) in ir
        st.expr isa Expr || continue
        st.expr.head == :call || continue
        (st.expr.args[1] == trace || unwrap(st.expr.args[1]) == :trace) || continue
        st.expr.args[2] isa QuoteNode && continue
        addr = st.expr.args[2]
        if addr isa Variable
            st = ir[addr]
            vars = Set{Variable}([])
            if st.expr isa Expr
                MacroTools.postwalk(st.expr) do e
                    e isa Variable && push!(vars, addr)
                    e
                end
            end
            ir_map[addr] = st
            append!(dynamic_addresses, vars)
        end
    end
    dh = DynamicAddressingHint(dynamic_addresses, ir_map)
    dh
end

function detect_dynamic_addresses(fn, arg_types...)
    ir = lower_to_ir(fn, arg_types...)
    display(detect_dynamic_addresses(ir))
end

@inline dynamic_address_check(ir) = !isempty(detect_dynamic_addresses(ir))


text:
@testset "331.verify-preorder-serialization-of-a-binary-tree.jl" begin
    @test is_valid_serialization("9,3,4,#,#,1,#,#,2,#,6,#,#") == true
    @test is_valid_serialization("1,#") == false
    @test is_valid_serialization("9,#,#,1") == false
end


text:
import Tables

struct AlphaVantageResponse
    data::Vector{AbstractVector}
    names::Vector{AbstractString}

    AlphaVantageResponse(data::Vector{AbstractVector}, names::Vector{<:AbstractString}) = begin
        l1 = length(data[1])
        @assert all(t -> length(t) == l1, data)
        @assert length(data) == length(names)
        new(data, names)
    end
end

function narrow_types(v::AbstractVector)
    T = mapreduce(typeof, promote_type, v)
    convert(AbstractVector{T}, v)
end

AlphaVantageResponse(data::AbstractVector{<:AbstractVector{T}} where T, names::AbstractVector{<:AbstractString}) = begin
    AlphaVantageResponse(narrow_types(data), names)
end

AlphaVantageResponse(data::AbstractMatrix{T} where T, names::AbstractMatrix{<:AbstractString}) = begin
    v = AbstractVector[narrow_types(c) for c in eachcol(data)]
    n = vec(names)
    AlphaVantageResponse(v, n)
end

AlphaVantageResponse(raw::Tuple{AbstractMatrix{Any}, AbstractMatrix{<:AbstractString}}) = begin
    AlphaVantageResponse(raw[1], raw[2])
end

AlphaVantageResponse(d::Dict) = Tables.table(reshape(collect(values(d)), (1,:)), header=keys(d))

Tables.istable(::AlphaVantageResponse) = true
Tables.rowaccess(::AlphaVantageResponse) = false
Tables.columnaccess(::AlphaVantageResponse) = true
Tables.columns(t::AlphaVantageResponse) = t
Tables.getcolumn(t::AlphaVantageResponse, i::Int) = t.data[i]
Tables.getcolumn(t::AlphaVantageResponse, nm::Symbol) = begin
    ind = findfirst(==(nm), Symbol.(t.names))
    t.data[ind]
end
Tables.columnnames(t::AlphaVantageResponse) = Symbol.(t.names)


text:
function Base.serialize(S::AbstractSerializer, d::DArray)
    # Only send the ident for participating workers - we expect the DArray to exist in the
    # remote registry
    destpid = Base.worker_id_from_socket(S.io)
    Serializer.serialize_type(S, typeof(d))
    if (destpid in d.pids) || (destpid == d.identity[1])
        serialize(S, (true, d.identity))    # (identity_only, identity)
    else
        serialize(S, (false, d.identity))
        for n in [:dims, :pids, :indexes, :cuts]
            serialize(S, getfield(d, n))
        end
    end
end

function Base.deserialize{T<:DArray}(S::AbstractSerializer, t::Type{T})
    what = deserialize(S)
    identity_only = what[1]
    identity = what[2]

    if identity_only
        global registry
        if haskey(registry, (identity, :DARRAY))
            return registry[(identity, :DARRAY)]
        else
            # access to fields will throw an error, at least the deserialization process will not
            # result in worker death
            d = T()
            d.identity = identity
            return d
        end
    else
        # We are not a participating worker, deser fields and instantiate locally.
        dims = deserialize(S)
        pids = deserialize(S)
        indexes = deserialize(S)
        cuts = deserialize(S)
        return T(identity, dims, pids, indexes, cuts)
    end
end

# Serialize only those parts of the object as required by the destination worker.
type DestinationSerializer
    generate::Nullable{Function}     # Function to generate the part to be serialized
    pids::Nullable{Array}            # MUST have the same shape as the distribution

    deser_obj::Nullable{Any}         # Deserialized part

    DestinationSerializer(f,p,d) = new(f,p,d)
end

DestinationSerializer(f::Function, pids::Array) = DestinationSerializer(f, pids, Nullable{Any}())

# contructs a DestinationSerializer after verifying that the shape of pids.
function verified_destination_serializer(f::Function, pids::Array, verify_size)
    @assert size(pids) == verify_size
    return DestinationSerializer(f, pids)
end

DestinationSerializer(deser_obj::Any) = DestinationSerializer(Nullable{Function}(), Nullable{Array}(), deser_obj)

function Base.serialize(S::AbstractSerializer, s::DestinationSerializer)
    pid = Base.worker_id_from_socket(S.io)
    pididx = findfirst(get(s.pids), pid)
    Serializer.serialize_type(S, typeof(s))
    serialize(S, get(s.generate)(pididx))
end

function Base.deserialize{T<:DestinationSerializer}(S::AbstractSerializer, t::Type{T})
    lpart = deserialize(S)
    return DestinationSerializer(lpart)
end

function localpart(s::DestinationSerializer)
    if !isnull(s.deser_obj)
        return get(s.deser_obj)
    elseif  !isnull(s.generate) && (myid() in get(s.pids))
        # Handle the special case where myid() is part of s.pids.
        # In this case serialize/deserialize is not called as the remotecall is executed locally
        return get(s.generate)(findfirst(get(s.pids), myid()))
    else
        throw(ErrorException(string("Invalid state in DestinationSerializer.")))
    end
end


text:
# Lets use black-box optimization to do regression!
#
# To run this, just do:
#   julia regression_via_optimization.jl
# from a command line after having installed BlackBoxOptim.
#
# The basic idea is that a regression problem is a problem of deciding
# which set of coefficients (betas) minimize an objective function
# that compares the dependent variable y to a model based on the betas and the
# independent variables x.
#
# When implementing this, x will be a matrix d*n of floats where d is the dimension
# of the problem and n is the number of examples/cases. Thus y is an array/vector
# 1*n of floats and the coefficients is an array/vector (d+1)*1 of coefficients
# where the first is the intercept (beta0) while the rest are the coefficients.
#
# We get different types of regression depending on which objective function
# we select. Common to many of them is that they first calculate the deviances
# between the model and the actual values the model should predict:
function discrepancies(betas, x, y)
  y' .- ( betas[1] + sum(broadcast(*, betas[2:end], x), 1) )
end

# Given this setup we can now create an objective function for
# Ordinare Least Squares (OLS) regression. This is actually just the L2 norm:
function ols_regression_objective(beta, x, y)
  norm(discrepancies(beta, x, y), 2)
end

# And the L1 norm gives Least Absolute Deviations (LAD) (aka Robust regression) regression:
function lad_regression_objective(beta, x, y)
  norm(discrepancies(beta, x, y), 1)
end

# And we can do LASSO and Ridge regression by adding a penalty on large coefficients, but
# note that these also take a lambda constant:
function regularized_regression_objective(lambda, beta, x, y, p = 2, q = 1)
  norm(discrepancies(beta, x, y), p) + lambda * norm(beta[2:end], q)
end

# By selecting p = 2 and q = 1 we get LASSO regression:
function lasso_regression_objective(lambda, beta, x, y)
  regularized_regression_objective(lambda, beta, x, y, 2, 1)
end

# and by selecting p = 2 and q = 2 we get Ridge regression:
function ridge_regression_objective(lambda, beta, x, y)
  regularized_regression_objective(lambda, beta, x, y, 2, 2)
end

# Ok, we're ready to do some regression. We need some data, lets start with
# a way to generate data from models specified as julia functions. Lets start
# with a simple model of 3 vars: X1 + 2*X2 - X3
function m1(x)
  x[1,:] + 2*x[2,:] - x[3,:]
end

# and we generate random input for it:
x1 = rand(3, 100)
# and calc the outputs:
y1 = m1(x1)

# We can now search for the coefficients that minimize the OLS objective
# with a black-box optimization search like so (we allow coefficients to have
# a min value of -5.0 and a max value of 5.0, and the search is for 4
# coefficients, one intercept and three for each of the values of x):
using BlackBoxOptim
olsresult = bboptimize(b -> ols_regression_objective(b, x1, y1);
  SearchRange = (-10.0, 10.0), NumDimensions = 4, MaxSteps = 2e4)
@show (best_fitness(olsresult), best_candidate(olsresult))

# But the really nice thing is that we can easily consider other objectives such as the LAD:
ladresult = bboptimize(b -> lad_regression_objective(b, x1, y1);
  SearchRange = (-10.0, 10.0), NumDimensions = 4, MaxSteps = 2e4)
@show (best_fitness(ladresult), best_candidate(ladresult))

# For regularized regression we can optimize for different values of lambda so
# create a wrapper function that handles this:
function regularized_opt(lambda, func, x, y, dims, its = 2e4)
  bboptimize((b) -> func(lambda, b, x, y);
    SearchRange = (-10.0, 10.0), NumDimensions = dims, MaxSteps = its)
end

lassores1 = regularized_opt(1, lasso_regression_objective, x1, y1, 4)
lassores2 = regularized_opt(2, lasso_regression_objective, x1, y1, 4)
lassores3 = regularized_opt(3, lasso_regression_objective, x1, y1, 4)

ridgeres1 = regularized_opt(1, ridge_regression_objective, x1, y1, 4)
ridgeres2 = regularized_opt(2, ridge_regression_objective, x1, y1, 4)
ridgeres3 = regularized_opt(3, ridge_regression_objective, x1, y1, 4)

# Now lets create some support functions for printing models nicely.
linear_terms(num) = [@sprintf(" * X%d", i) for i in 1:num]
squared_terms(num) = [@sprintf(" * X%d^2", i) for i in 1:num]
linsq_terms(num) = vcat(linear_terms(num), squared_terms(num))

function sprint_predicted_model(fitresult, terms = nothing, skipIfLower = 1e-5)
  bestfit = best_candidate(fitresult)
  if terms == nothing
    terms = linear_terms(length(bestfit)-1)
  end
  if length(terms) < length(bestfit)
    terms = vcat([""], terms) # Put an empty term first which corresponds to intercept
  end

  signstr(value) = (value < 0.0) ? (@sprintf(" - %.3f", abs(value))) : (@sprintf(" + %.3f", value))

  elems = Any[]
  first_push = true
  for i in 1:length(bestfit)
    if abs(bestfit[i]) > skipIfLower
      str = join([signstr(bestfit[i]), terms[i]])
      if first_push
        # Strip away leading sign since this it the first push
        push!(elems, str[3:end])
        first_push = false
      else
        push!(elems, str)
      end
    end
  end

  join(elems)
end

# Let's try a model which involves a squared terms:
#   X1 + 4.13*X2*X2 - 3.14*X3
function m2(x)
  x[1,:] + (4.13 * x[2,:].^2) - (3.14 * x[3,:])
end

# and we generate random input for it:
x2 = rand(3, 100)
# and calc the outputs:
y2 = m2(x2)

# Before we regress we need to encode our beliefs about the general
# structure of the model. Let's say we believe there are squared terms but we
# do not know which ones. So we add one squared term per independent variable:
x2m = zeros(3+3, 100)
x2m[1:3,:] = x2
x2m[4,:] = x2[1,:].^2
x2m[5,:] = x2[2,:].^2
x2m[6,:] = x2[3,:].^2

# With this we can fit models:
m2_olsres = bboptimize((b) -> ols_regression_objective(b, x2m, y2);
  SearchRange = (-10.0, 10.0), NumDimensions = 7, MaxSteps = 5e4)
m2_ladres = bboptimize((b) -> lad_regression_objective(b, x2m, y2);
  SearchRange = (-10.0, 10.0), NumDimensions = 7, MaxSteps = 5e4)
m2_lassores1 = regularized_opt(1, lasso_regression_objective, x2m, y2, 7, 5e4)
m2_lassores2 = regularized_opt(2, lasso_regression_objective, x2m, y2, 7, 5e4)
m2_lassores3 = regularized_opt(3, lasso_regression_objective, x2m, y2, 7, 5e4)
m2_ridgeres1 = regularized_opt(1, ridge_regression_objective, x2m, y2, 7, 5e4)
m2_ridgeres2 = regularized_opt(2, ridge_regression_objective, x2m, y2, 7, 5e4)
m2_ridgeres3 = regularized_opt(3, ridge_regression_objective, x2m, y2, 7, 5e4)

# And now lets print our models nicely so user can see the results...

println("Model1 = 1.000 * X1 + 2.000 * X2 - 1.000 * X3")
println("OLS best fit: ", sprint_predicted_model(olsresult))
println("LAD best fit: ", sprint_predicted_model(ladresult))
println("LASSO best fit, lambda = 1: ", sprint_predicted_model(lassores1))
println("LASSO best fit, lambda = 2: ", sprint_predicted_model(lassores2))
println("LASSO best fit, lambda = 3: ", sprint_predicted_model(lassores3))
println("Ridge best fit, lambda = 1: ", sprint_predicted_model(ridgeres1))
println("Ridge best fit, lambda = 2: ", sprint_predicted_model(ridgeres2))
println("Ridge best fit, lambda = 3: ", sprint_predicted_model(ridgeres3))

terms = linsq_terms(3)
println("")
println("Model2 = 1.000 * X1 - 3.140 * X3 + 4.130 * X2^2 ")
println("OLS best fit: ", sprint_predicted_model(m2_olsres, terms))
println("LAD best fit: ", sprint_predicted_model(m2_ladres, terms))
println("LASSO best fit, lambda = 1: ", sprint_predicted_model(m2_lassores1, terms))
println("LASSO best fit, lambda = 2: ", sprint_predicted_model(m2_lassores2, terms))
println("LASSO best fit, lambda = 3: ", sprint_predicted_model(m2_lassores3, terms))
println("Ridge best fit, lambda = 1: ", sprint_predicted_model(m2_ridgeres1, terms))
println("Ridge best fit, lambda = 2: ", sprint_predicted_model(m2_ridgeres2, terms))
println("Ridge best fit, lambda = 3: ", sprint_predicted_model(m2_ridgeres3, terms))

# Conclusion: With black-box optimization you can easily fit regression models
# from very different paradigms without having to implement very much code.


text:
function _majorVersion()
    ccall((:CSIPmajorVersion, libcsip), Cint, ())
end

function _minorVersion()
    ccall((:CSIPminorVersion, libcsip), Cint, ())
end

function _patchVersion()
    ccall((:CSIPpatchVersion, libcsip), Cint, ())
end

function _freeModel(model::SCIPMathProgModel)
    ccall((:CSIPfreeModel, libcsip), Cint, (Ptr{Void}, ), model.ptr_model)
end

function _addVar(model::SCIPMathProgModel, lowerbound::Cdouble,
                 upperbound::Cdouble, vartype::Cint, idx::Ptr{Cint})
    ccall((:CSIPaddVar, libcsip), Cint,
          (Ptr{Void}, Cdouble, Cdouble, Cint, Ptr{Cint}),
          model.ptr_model, lowerbound, upperbound, vartype, idx)
end

function _chgVarLB(model::SCIPMathProgModel, numindices::Cint,
                   indices::Vector{Cint}, lowerbounds::Vector{Cdouble})
    ccall((:CSIPchgVarLB, libcsip), Cint,
          (Ptr{Void}, Cint, Ptr{Cint}, Ptr{Cdouble}),
          model.ptr_model, numindices, indices, lowerbounds)
end

function _chgVarUB(model::SCIPMathProgModel, numindices::Cint,
                   indices::Vector{Cint}, upperbounds::Vector{Cdouble})
    ccall((:CSIPchgVarUB, libcsip), Cint,
          (Ptr{Void}, Cint, Ptr{Cint}, Ptr{Cdouble}),
          model.ptr_model, numindices, indices, upperbounds)
end

function _getVarType(model::SCIPMathProgModel, varindex::Cint)
    ccall((:CSIPgetVarType, libcsip), Cint,
          (Ptr{Void}, Cint),
          model.ptr_model, varindex)
end

function _chgVarType(model::SCIPMathProgModel, varindex::Cint, vartype::Cint)
    ccall((:CSIPchgVarType, libcsip), Cint,
          (Ptr{Void}, Cint, Cint),
          model.ptr_model, varindex, vartype)
end

function _addLinCons(model::SCIPMathProgModel, numindices::Cint,
                     indices::Vector{Cint}, coefs::Vector{Cdouble},
                     lhs::Cdouble, rhs::Cdouble, idx::Ptr{Cint})
    ccall((:CSIPaddLinCons, libcsip), Cint,
          (Ptr{Void}, Cint, Ptr{Cint}, Ptr{Cdouble}, Cdouble, Cdouble, Ptr{Cint}),
          model.ptr_model, numindices, indices, coefs, lhs, rhs, idx)
end

function _addQuadCons(model::SCIPMathProgModel, numlinindices::Cint,
                      linindices::Vector{Cint}, lincoefs::Vector{Cdouble},
                      numquadterms::Cint, quadrowindices::Vector{Cint},
                      quadcolindices::Vector{Cint}, quadcoefs::Vector{Cdouble},
                      lhs::Cdouble, rhs::Cdouble, idx::Ptr{Cint})
    ccall((:CSIPaddQuadCons, libcsip), Cint,
          (Ptr{Void}, Cint, Ptr{Cint}, Ptr{Cdouble}, Cint, Ptr{Cint}, Ptr{Cint},
           Ptr{Cdouble}, Cdouble, Cdouble, Ptr{Cint}),
          model.ptr_model, numlinindices, linindices, lincoefs, numquadterms,
          quadrowindices, quadcolindices, quadcoefs, lhs, rhs, idx)
end

function _addNonLinCons(model::SCIPMathProgModel, nops::Cint,
                     ops::Vector{Cint}, children::Vector{Cint},
                     beg::Vector{Cint}, values::Vector{Cdouble},
                     lhs::Cdouble, rhs::Cdouble, idx::Ptr{Cint})
    ccall((:CSIPaddNonLinCons, libcsip), Cint,
          (Ptr{Void}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Cdouble, Cdouble, Ptr{Cint}),
          model.ptr_model, nops, ops, children, beg, values, lhs, rhs, idx)
end

function _addSOS1(model::SCIPMathProgModel, numindices::Cint,
                  indices::Vector{Cint}, weights::Vector{Cdouble},
                  idx::Ptr{Cint})
    ccall((:CSIPaddSOS1, libcsip), Cint,
          (Ptr{Void}, Cint, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}),
          model.ptr_model, numindices, indices, weights, idx)
end

function _addSOS2(model::SCIPMathProgModel, numindices::Cint,
                  indices::Vector{Cint}, weights::Vector{Cdouble},
                  idx::Ptr{Cint})
    ccall((:CSIPaddSOS2, libcsip), Cint,
          (Ptr{Void}, Cint, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}),
          model.ptr_model, numindices, indices, weights, idx)
end

function _setObj(model::SCIPMathProgModel, numindices::Cint,
                 indices::Vector{Cint}, coefs::Vector{Cdouble})
    ccall((:CSIPsetObj, libcsip), Cint,
          (Ptr{Void}, Cint, Ptr{Cint}, Ptr{Cdouble}),
          model.ptr_model, numindices, indices, coefs)
end

function _setQuadObj(model::SCIPMathProgModel, numlinindices::Cint,
                     linindices::Vector{Cint}, lincoefs::Vector{Cdouble},
                     numquadterms::Cint, quadrowindices::Vector{Cint},
                     quadcolindices::Vector{Cint}, quadcoefs::Vector{Cdouble})
    ccall((:CSIPsetQuadObj, libcsip), Cint,
          (Ptr{Void}, Cint, Ptr{Cint}, Ptr{Cdouble}, Cint, Ptr{Cint}, Ptr{Cint},
           Ptr{Cdouble}),
          model.ptr_model, numlinindices, linindices, lincoefs, numquadterms,
          quadrowindices, quadcolindices, quadcoefs)
end

function _setNonlinearObj(model::SCIPMathProgModel, nops::Cint,
                          ops::Vector{Cint}, children::Vector{Cint},
                          beg::Vector{Cint}, values::Vector{Cdouble})
    ccall((:CSIPsetNonlinearObj, libcsip), Cint,
          (Ptr{Void}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
          model.ptr_model, nops, ops, children, beg, values)
end

function _setSenseMinimize(model::SCIPMathProgModel)
    ccall((:CSIPsetSenseMinimize, libcsip), Cint, (Ptr{Void}, ), model.ptr_model)
end

function _setSenseMaximize(model::SCIPMathProgModel)
    ccall((:CSIPsetSenseMaximize, libcsip), Cint, (Ptr{Void}, ), model.ptr_model)
end

function _solve(model::SCIPMathProgModel)
    ccall((:CSIPsolve, libcsip), Cint, (Ptr{Void}, ), model.ptr_model)
end

function _interrupt(model::SCIPMathProgModel)
    ccall((:CSIPinterrupt, libcsip), Cint, (Ptr{Void}, ), model.ptr_model)
end

function _getVarValues(model::SCIPMathProgModel, output::Vector{Cdouble})
    ccall((:CSIPgetVarValues, libcsip), Cint, (Ptr{Void}, Ptr{Cdouble}),
          model.ptr_model, output)
end

function _getObjValue(model::SCIPMathProgModel)
    ccall((:CSIPgetObjValue, libcsip), Cdouble, (Ptr{Void}, ), model.ptr_model)
end

function _getObjBound(model::SCIPMathProgModel)
    ccall((:CSIPgetObjBound, libcsip), Cdouble, (Ptr{Void}, ), model.ptr_model)
end

function _getStatus(model::SCIPMathProgModel)
    ccall((:CSIPgetStatus, libcsip), Cint, (Ptr{Void}, ), model.ptr_model)
end

function _getParamType(model::SCIPMathProgModel, name::Ptr{UInt8})
    ccall((:CSIPgetParamType, libcsip), Cint,
          (Ptr{Void}, Ptr{UInt8}),
          model.ptr_model, name)
end

function _setBoolParam(model::SCIPMathProgModel, name::Ptr{UInt8},
                              value::Cint)
    ccall((:CSIPsetBoolParam, libcsip), Cint,
          (Ptr{Void}, Ptr{UInt8}, Cint),
          model.ptr_model, name, value)
end

function _setIntParam(model::SCIPMathProgModel, name::Ptr{UInt8},
                              value::Cint)
    ccall((:CSIPsetIntParam, libcsip), Cint,
          (Ptr{Void}, Ptr{UInt8}, Cint),
          model.ptr_model, name, value)
end

function _setLongintParam(model::SCIPMathProgModel, name::Ptr{UInt8},
                              value::Clonglong)
    ccall((:CSIPsetLongintParam, libcsip), Cint,
          (Ptr{Void}, Ptr{UInt8}, Clonglong),
          model.ptr_model, name, value)
end

function _setRealParam(model::SCIPMathProgModel, name::Ptr{UInt8},
                              value::Cdouble)
    ccall((:CSIPsetRealParam, libcsip), Cint,
          (Ptr{Void}, Ptr{UInt8}, Cdouble),
          model.ptr_model, name, value)
end

function _setCharParam(model::SCIPMathProgModel, name::Ptr{UInt8},
                              value::Cchar)
    ccall((:CSIPsetCharParam, libcsip), Cint,
          (Ptr{Void}, Ptr{UInt8}, Cchar),
          model.ptr_model, name, value)
end

function _setStringParam(model::SCIPMathProgModel, name::Ptr{UInt8},
                         value::Ptr{UInt8})
    ccall((:CSIPsetStringParam, libcsip), Cint,
          (Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}),
          model.ptr_model, name, value)
end

function _getNumVars(model::SCIPMathProgModel)
    ccall((:CSIPgetNumVars, libcsip), Cint, (Ptr{Void}, ), model.ptr_model)
end

function _getNumConss(model::SCIPMathProgModel)
    ccall((:CSIPgetNumConss, libcsip), Cint, (Ptr{Void}, ), model.ptr_model)
end

function _setInitialSolution(model::SCIPMathProgModel, values::Vector{Cdouble})
    ccall((:CSIPsetInitialSolution, libcsip), Cint, (Ptr{Void}, Ptr{Cdouble}),
          model.ptr_model, values)
end

function _lazyGetContext(lazydata::Ptr{Void})
    ccall((:CSIPlazyGetContext, libcsip), Cint, (Ptr{Void}, ), lazydata)
end

function _lazyGetVarValues(lazydata::Ptr{Void}, output::Vector{Cdouble})
    ccall((:CSIPlazyGetVarValues, libcsip), Cint, (Ptr{Void}, Ptr{Cdouble}),
          lazydata, output)
end

function _lazyAddLinCons(lazydata::Ptr{Void}, numindices::Cint,
                       indices::Vector{Cint}, coefs::Vector{Cdouble},
                       lhs::Cdouble, rhs::Cdouble, islocal::Cint)
    ccall((:CSIPlazyAddLinCons, libcsip), Cint,
          (Ptr{Void}, Cint, Ptr{Cint}, Ptr{Cdouble}, Cdouble, Cdouble, Cint),
          lazydata, numindices, indices, coefs, lhs, rhs, islocal)
end

function _addLazyCallback(model::SCIPMathProgModel, lazycb::Ptr{Void}, userdata)
    ccall((:CSIPaddLazyCallback, libcsip), Cint,
          (Ptr{Void}, Ptr{Void}, Any),
          model.ptr_model, lazycb, userdata)
end

function _heurGetVarValues(heurdata::Ptr{Void}, output::Vector{Cdouble})
    ccall((:CSIPheurGetVarValues, libcsip), Cint, (Ptr{Void}, Ptr{Cdouble}),
          heurdata, output)
end

function _heurAddSolution(heurdata::Ptr{Void}, values::Vector{Cdouble})
    ccall((:CSIPheurAddSolution, libcsip), Cint, (Ptr{Void}, Ptr{Cdouble}),
          heurdata, values)
end

function _addHeuristicCallback(model::SCIPMathProgModel, heur::Ptr{Void}, userdata)
    ccall((:CSIPaddHeuristicCallback, libcsip), Cint, (Ptr{Void}, Ptr{Void}, Any),
          model.ptr_model, heur, userdata)
end

function _getInternalSCIP(model::SCIPMathProgModel)
    ccall((:CSIPgetInternalSCIP, libcsip), Ptr{Void}, (Ptr{Void}, ),
          model.ptr_model)
end


text:
struct Linear
end
struct Left
end
struct Right
end

"""
    interpolate(method, X::Trajectory, s)

Interpolate trajectory `X` using `method in [Linear(), Left(), ...]`
at time `s`.
"""
function interpolate
end

function interpolate(::Linear, X::Trajectory, s)
    i = searchsorted(keys(X), s)
    t, x = Pair(X)
    n = length(X)
    f, l = first(i), last(i)
    if f <= 1
        return x[1]
    elseif l >= n
        return x[end]
    elseif f <= l
        return x[f]
    else
        λ = (t[f] - s)/(t[f]-t[l])
        return λ*x[l] + (1-λ)*x[f]
    end
end

function interpolate(::Left, X::Trajectory, s)
    i = searchsorted(keys(X), s)
    t, x = Pair(X)
    n = length(X)
    f, l = first(i), last(i)
    return x[l]
end

function interpolate(::Right, X::Trajectory, s)
    i = searchsorted(keys(X), s)
    t, x = Pair(X)
    n = length(X)
    f, l = first(i), last(i)
    return x[f]
end


text:
## test_show.jl : test show commands

@testset "display commands         " begin

    str = string(Matern(1,1))
    @test str isa String
    str = string(Exponential(1))
    @test str isa String
    str = string(Whittle(1))
    @test str isa String
    str = string(Linear(1))
    @test str isa String
    str = string(Spherical(1))
    @test str isa String
    str = string(SquaredExponential(1))
    @test str isa String
    str = string(AnisotropicExponential([1000 0; 0 1000]))
    @test str isa String
    str = string(CovarianceFunction(2,Matern(1.0,2)))
    @test str isa String
    str = string(Cholesky())
    @test str isa String
    str = string(Spectral())
    @test str isa String
    str = string(KarhunenLoeve(10))
    @test str isa String
    str = string(CirculantEmbedding())
    @test str isa String
    str = string(GaussianRandomField(CovarianceFunction(1,Exponential(1)),Cholesky(),0:0.1:1))
    @test str isa String
    str = string(SeparableCovarianceFunction(Exponential(1),Exponential(1)))
    @test str isa String
    str = string(GaussianRandomField(CovarianceFunction(2,Exponential(0.1)),Cholesky(),0:0.1:1,0:0.1:1))
    @test str isa String
    str = string(GaussianRandomField(CovarianceFunction(2,Exponential(0.1)),KarhunenLoeve(10),rand(25,2)))
    @test str isa String
    str = string(GaussianRandomField(CovarianceFunction(2,Exponential(0.1)),KarhunenLoeve(10),rand(3,2),collect([1,2,3]')))
    @test str isa String

end


text:
#Script for including all the other functions

  include("getalpha.jl")
  include("numerovfunc.jl")
  include("GetAtomParams.jl")
  include("Rb87numbers.jl")
  include("filenamemake.jl")
  include("SIunits.jl")
  include("krondelt.jl")
  include("triangle_coeff.jl")
  include("Wigner6j.jl")
  include("Wigner6jcheck.jl")
  include("fung.jl")
  include("tfunction3j.jl")
  include("Wigner3j.jl")
  include("ClebschGord.jl")
  include("radiel.jl")
  include("matrix_elements.jl")
  include("Radiative_lifetimes.jl")
  include("Rb87DecayPaths.jl")
  include("Rb87blackbody.jl")
  include("firstexcitedrabiRb.jl")
  include("matrixel2p.jl")
  include("angcalc.jl")
  include("BlockadeShift.jl")
  include("rectpulse.jl")
  include("combvechack.jl")
  include("envalfunc.jl")

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


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

"""
    CreateCluster()

Creates a DAX cluster. All nodes in the cluster run the same DAX caching software.

# Required Parameters
- `ClusterName`: The cluster identifier. This parameter is stored as a lowercase string.  Constraints:    A name must contain from 1 to 20 alphanumeric characters or hyphens.   The first character must be a letter.   A name cannot end with a hyphen or contain two consecutive hyphens.
- `IamRoleArn`: A valid Amazon Resource Name (ARN) that identifies an IAM role. At runtime, DAX will assume this role and use the role's permissions to access DynamoDB on your behalf.
- `NodeType`: The compute and memory capacity of the nodes in the cluster.
- `ReplicationFactor`: The number of nodes in the DAX cluster. A replication factor of 1 will create a single-node cluster, without any read replicas. For additional fault tolerance, you can create a multiple node cluster with one or more read replicas. To do this, set ReplicationFactor to a number between 3 (one primary and two read replicas) and 10 (one primary and nine read replicas). If the AvailabilityZones parameter is provided, its length must equal the ReplicationFactor.  AWS recommends that you have at least two read replicas per cluster.

# Optional Parameters
- `AvailabilityZones`: The Availability Zones (AZs) in which the cluster nodes will reside after the cluster has been created or updated. If provided, the length of this list must equal the ReplicationFactor parameter. If you omit this parameter, DAX will spread the nodes across Availability Zones for the highest availability.
- `Description`: A description of the cluster.
- `NotificationTopicArn`: The Amazon Resource Name (ARN) of the Amazon SNS topic to which notifications will be sent.  The Amazon SNS topic owner must be same as the DAX cluster owner.
- `ParameterGroupName`: The parameter group to be associated with the DAX cluster.
- `PreferredMaintenanceWindow`: Specifies the weekly time range during which maintenance on the DAX cluster is performed. It is specified as a range in the format ddd:hh24:mi-ddd:hh24:mi (24H Clock UTC). The minimum maintenance window is a 60 minute period. Valid values for ddd are:    sun     mon     tue     wed     thu     fri     sat    Example: sun:05:00-sun:09:00   If you don't specify a preferred maintenance window when you create or modify a cache cluster, DAX assigns a 60-minute maintenance window on a randomly selected day of the week.
- `SSESpecification`: Represents the settings used to enable server-side encryption on the cluster.
- `SecurityGroupIds`: A list of security group IDs to be assigned to each node in the DAX cluster. (Each of the security group ID is system-generated.) If this parameter is not specified, DAX assigns the default VPC security group to each node.
- `SubnetGroupName`: The name of the subnet group to be used for the replication group.  DAX clusters can only run in an Amazon VPC environment. All of the subnets that you specify in a subnet group must exist in the same VPC.
- `Tags`: A set of tags to associate with the DAX cluster.
"""
create_cluster(ClusterName, IamRoleArn, NodeType, ReplicationFactor; aws_config::AWSConfig=global_aws_config()) = dax("CreateCluster", Dict{String, Any}("ClusterName"=>ClusterName, "IamRoleArn"=>IamRoleArn, "NodeType"=>NodeType, "ReplicationFactor"=>ReplicationFactor); aws_config=aws_config)
create_cluster(ClusterName, IamRoleArn, NodeType, ReplicationFactor, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("CreateCluster", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClusterName"=>ClusterName, "IamRoleArn"=>IamRoleArn, "NodeType"=>NodeType, "ReplicationFactor"=>ReplicationFactor), args)); aws_config=aws_config)

"""
    CreateParameterGroup()

Creates a new parameter group. A parameter group is a collection of parameters that you apply to all of the nodes in a DAX cluster.

# Required Parameters
- `ParameterGroupName`: The name of the parameter group to apply to all of the clusters in this replication group.

# Optional Parameters
- `Description`: A description of the parameter group.
"""
create_parameter_group(ParameterGroupName; aws_config::AWSConfig=global_aws_config()) = dax("CreateParameterGroup", Dict{String, Any}("ParameterGroupName"=>ParameterGroupName); aws_config=aws_config)
create_parameter_group(ParameterGroupName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("CreateParameterGroup", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ParameterGroupName"=>ParameterGroupName), args)); aws_config=aws_config)

"""
    CreateSubnetGroup()

Creates a new subnet group.

# Required Parameters
- `SubnetGroupName`: A name for the subnet group. This value is stored as a lowercase string.
- `SubnetIds`: A list of VPC subnet IDs for the subnet group.

# Optional Parameters
- `Description`: A description for the subnet group
"""
create_subnet_group(SubnetGroupName, SubnetIds; aws_config::AWSConfig=global_aws_config()) = dax("CreateSubnetGroup", Dict{String, Any}("SubnetGroupName"=>SubnetGroupName, "SubnetIds"=>SubnetIds); aws_config=aws_config)
create_subnet_group(SubnetGroupName, SubnetIds, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("CreateSubnetGroup", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("SubnetGroupName"=>SubnetGroupName, "SubnetIds"=>SubnetIds), args)); aws_config=aws_config)

"""
    DecreaseReplicationFactor()

Removes one or more nodes from a DAX cluster.  You cannot use DecreaseReplicationFactor to remove the last node in a DAX cluster. If you need to do this, use DeleteCluster instead.

# Required Parameters
- `ClusterName`: The name of the DAX cluster from which you want to remove nodes.
- `NewReplicationFactor`: The new number of nodes for the DAX cluster.

# Optional Parameters
- `AvailabilityZones`: The Availability Zone(s) from which to remove nodes.
- `NodeIdsToRemove`: The unique identifiers of the nodes to be removed from the cluster.
"""
decrease_replication_factor(ClusterName, NewReplicationFactor; aws_config::AWSConfig=global_aws_config()) = dax("DecreaseReplicationFactor", Dict{String, Any}("ClusterName"=>ClusterName, "NewReplicationFactor"=>NewReplicationFactor); aws_config=aws_config)
decrease_replication_factor(ClusterName, NewReplicationFactor, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("DecreaseReplicationFactor", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClusterName"=>ClusterName, "NewReplicationFactor"=>NewReplicationFactor), args)); aws_config=aws_config)

"""
    DeleteCluster()

Deletes a previously provisioned DAX cluster. DeleteCluster deletes all associated nodes, node endpoints and the DAX cluster itself. When you receive a successful response from this action, DAX immediately begins deleting the cluster; you cannot cancel or revert this action.

# Required Parameters
- `ClusterName`: The name of the cluster to be deleted.

"""
delete_cluster(ClusterName; aws_config::AWSConfig=global_aws_config()) = dax("DeleteCluster", Dict{String, Any}("ClusterName"=>ClusterName); aws_config=aws_config)
delete_cluster(ClusterName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("DeleteCluster", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClusterName"=>ClusterName), args)); aws_config=aws_config)

"""
    DeleteParameterGroup()

Deletes the specified parameter group. You cannot delete a parameter group if it is associated with any DAX clusters.

# Required Parameters
- `ParameterGroupName`: The name of the parameter group to delete.

"""
delete_parameter_group(ParameterGroupName; aws_config::AWSConfig=global_aws_config()) = dax("DeleteParameterGroup", Dict{String, Any}("ParameterGroupName"=>ParameterGroupName); aws_config=aws_config)
delete_parameter_group(ParameterGroupName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("DeleteParameterGroup", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ParameterGroupName"=>ParameterGroupName), args)); aws_config=aws_config)

"""
    DeleteSubnetGroup()

Deletes a subnet group.  You cannot delete a subnet group if it is associated with any DAX clusters.

# Required Parameters
- `SubnetGroupName`: The name of the subnet group to delete.

"""
delete_subnet_group(SubnetGroupName; aws_config::AWSConfig=global_aws_config()) = dax("DeleteSubnetGroup", Dict{String, Any}("SubnetGroupName"=>SubnetGroupName); aws_config=aws_config)
delete_subnet_group(SubnetGroupName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("DeleteSubnetGroup", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("SubnetGroupName"=>SubnetGroupName), args)); aws_config=aws_config)

"""
    DescribeClusters()

Returns information about all provisioned DAX clusters if no cluster identifier is specified, or about a specific DAX cluster if a cluster identifier is supplied. If the cluster is in the CREATING state, only cluster level information will be displayed until all of the nodes are successfully provisioned. If the cluster is in the DELETING state, only cluster level information will be displayed. If nodes are currently being added to the DAX cluster, node endpoint information and creation time for the additional nodes will not be displayed until they are completely provisioned. When the DAX cluster state is available, the cluster is ready for use. If nodes are currently being removed from the DAX cluster, no endpoint information for the removed nodes is displayed.

# Optional Parameters
- `ClusterNames`: The names of the DAX clusters being described.
- `MaxResults`: The maximum number of results to include in the response. If more results exist than the specified MaxResults value, a token is included in the response so that the remaining results can be retrieved. The value for MaxResults must be between 20 and 100.
- `NextToken`: An optional token returned from a prior request. Use this token for pagination of results from this action. If this parameter is specified, the response includes only results beyond the token, up to the value specified by MaxResults.
"""
describe_clusters(; aws_config::AWSConfig=global_aws_config()) = dax("DescribeClusters"; aws_config=aws_config)
describe_clusters(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("DescribeClusters", args; aws_config=aws_config)

"""
    DescribeDefaultParameters()

Returns the default system parameter information for the DAX caching software.

# Optional Parameters
- `MaxResults`: The maximum number of results to include in the response. If more results exist than the specified MaxResults value, a token is included in the response so that the remaining results can be retrieved. The value for MaxResults must be between 20 and 100.
- `NextToken`: An optional token returned from a prior request. Use this token for pagination of results from this action. If this parameter is specified, the response includes only results beyond the token, up to the value specified by MaxResults.
"""
describe_default_parameters(; aws_config::AWSConfig=global_aws_config()) = dax("DescribeDefaultParameters"; aws_config=aws_config)
describe_default_parameters(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("DescribeDefaultParameters", args; aws_config=aws_config)

"""
    DescribeEvents()

Returns events related to DAX clusters and parameter groups. You can obtain events specific to a particular DAX cluster or parameter group by providing the name as a parameter. By default, only the events occurring within the last 24 hours are returned; however, you can retrieve up to 14 days' worth of events if necessary.

# Optional Parameters
- `Duration`: The number of minutes' worth of events to retrieve.
- `EndTime`: The end of the time interval for which to retrieve events, specified in ISO 8601 format.
- `MaxResults`: The maximum number of results to include in the response. If more results exist than the specified MaxResults value, a token is included in the response so that the remaining results can be retrieved. The value for MaxResults must be between 20 and 100.
- `NextToken`: An optional token returned from a prior request. Use this token for pagination of results from this action. If this parameter is specified, the response includes only results beyond the token, up to the value specified by MaxResults.
- `SourceName`: The identifier of the event source for which events will be returned. If not specified, then all sources are included in the response.
- `SourceType`: The event source to retrieve events for. If no value is specified, all events are returned.
- `StartTime`: The beginning of the time interval to retrieve events for, specified in ISO 8601 format.
"""
describe_events(; aws_config::AWSConfig=global_aws_config()) = dax("DescribeEvents"; aws_config=aws_config)
describe_events(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("DescribeEvents", args; aws_config=aws_config)

"""
    DescribeParameterGroups()

Returns a list of parameter group descriptions. If a parameter group name is specified, the list will contain only the descriptions for that group.

# Optional Parameters
- `MaxResults`: The maximum number of results to include in the response. If more results exist than the specified MaxResults value, a token is included in the response so that the remaining results can be retrieved. The value for MaxResults must be between 20 and 100.
- `NextToken`: An optional token returned from a prior request. Use this token for pagination of results from this action. If this parameter is specified, the response includes only results beyond the token, up to the value specified by MaxResults.
- `ParameterGroupNames`: The names of the parameter groups.
"""
describe_parameter_groups(; aws_config::AWSConfig=global_aws_config()) = dax("DescribeParameterGroups"; aws_config=aws_config)
describe_parameter_groups(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("DescribeParameterGroups", args; aws_config=aws_config)

"""
    DescribeParameters()

Returns the detailed parameter list for a particular parameter group.

# Required Parameters
- `ParameterGroupName`: The name of the parameter group.

# Optional Parameters
- `MaxResults`: The maximum number of results to include in the response. If more results exist than the specified MaxResults value, a token is included in the response so that the remaining results can be retrieved. The value for MaxResults must be between 20 and 100.
- `NextToken`: An optional token returned from a prior request. Use this token for pagination of results from this action. If this parameter is specified, the response includes only results beyond the token, up to the value specified by MaxResults.
- `Source`: How the parameter is defined. For example, system denotes a system-defined parameter.
"""
describe_parameters(ParameterGroupName; aws_config::AWSConfig=global_aws_config()) = dax("DescribeParameters", Dict{String, Any}("ParameterGroupName"=>ParameterGroupName); aws_config=aws_config)
describe_parameters(ParameterGroupName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("DescribeParameters", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ParameterGroupName"=>ParameterGroupName), args)); aws_config=aws_config)

"""
    DescribeSubnetGroups()

Returns a list of subnet group descriptions. If a subnet group name is specified, the list will contain only the description of that group.

# Optional Parameters
- `MaxResults`: The maximum number of results to include in the response. If more results exist than the specified MaxResults value, a token is included in the response so that the remaining results can be retrieved. The value for MaxResults must be between 20 and 100.
- `NextToken`: An optional token returned from a prior request. Use this token for pagination of results from this action. If this parameter is specified, the response includes only results beyond the token, up to the value specified by MaxResults.
- `SubnetGroupNames`: The name of the subnet group.
"""
describe_subnet_groups(; aws_config::AWSConfig=global_aws_config()) = dax("DescribeSubnetGroups"; aws_config=aws_config)
describe_subnet_groups(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("DescribeSubnetGroups", args; aws_config=aws_config)

"""
    IncreaseReplicationFactor()

Adds one or more nodes to a DAX cluster.

# Required Parameters
- `ClusterName`: The name of the DAX cluster that will receive additional nodes.
- `NewReplicationFactor`: The new number of nodes for the DAX cluster.

# Optional Parameters
- `AvailabilityZones`: The Availability Zones (AZs) in which the cluster nodes will be created. All nodes belonging to the cluster are placed in these Availability Zones. Use this parameter if you want to distribute the nodes across multiple AZs.
"""
increase_replication_factor(ClusterName, NewReplicationFactor; aws_config::AWSConfig=global_aws_config()) = dax("IncreaseReplicationFactor", Dict{String, Any}("ClusterName"=>ClusterName, "NewReplicationFactor"=>NewReplicationFactor); aws_config=aws_config)
increase_replication_factor(ClusterName, NewReplicationFactor, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("IncreaseReplicationFactor", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClusterName"=>ClusterName, "NewReplicationFactor"=>NewReplicationFactor), args)); aws_config=aws_config)

"""
    ListTags()

List all of the tags for a DAX cluster. You can call ListTags up to 10 times per second, per account.

# Required Parameters
- `ResourceName`: The name of the DAX resource to which the tags belong.

# Optional Parameters
- `NextToken`: An optional token returned from a prior request. Use this token for pagination of results from this action. If this parameter is specified, the response includes only results beyond the token.
"""
list_tags(ResourceName; aws_config::AWSConfig=global_aws_config()) = dax("ListTags", Dict{String, Any}("ResourceName"=>ResourceName); aws_config=aws_config)
list_tags(ResourceName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("ListTags", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceName"=>ResourceName), args)); aws_config=aws_config)

"""
    RebootNode()

Reboots a single node of a DAX cluster. The reboot action takes place as soon as possible. During the reboot, the node status is set to REBOOTING.   RebootNode restarts the DAX engine process and does not remove the contents of the cache.

# Required Parameters
- `ClusterName`: The name of the DAX cluster containing the node to be rebooted.
- `NodeId`: The system-assigned ID of the node to be rebooted.

"""
reboot_node(ClusterName, NodeId; aws_config::AWSConfig=global_aws_config()) = dax("RebootNode", Dict{String, Any}("ClusterName"=>ClusterName, "NodeId"=>NodeId); aws_config=aws_config)
reboot_node(ClusterName, NodeId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("RebootNode", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClusterName"=>ClusterName, "NodeId"=>NodeId), args)); aws_config=aws_config)

"""
    TagResource()

Associates a set of tags with a DAX resource. You can call TagResource up to 5 times per second, per account.

# Required Parameters
- `ResourceName`: The name of the DAX resource to which tags should be added.
- `Tags`: The tags to be assigned to the DAX resource.

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

"""
    UntagResource()

Removes the association of tags from a DAX resource. You can call UntagResource up to 5 times per second, per account.

# Required Parameters
- `ResourceName`: The name of the DAX resource from which the tags should be removed.
- `TagKeys`: A list of tag keys. If the DAX cluster has any tags with these keys, then the tags are removed from the cluster.

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

"""
    UpdateCluster()

Modifies the settings for a DAX cluster. You can use this action to change one or more cluster configuration parameters by specifying the parameters and the new values.

# Required Parameters
- `ClusterName`: The name of the DAX cluster to be modified.

# Optional Parameters
- `Description`: A description of the changes being made to the cluster.
- `NotificationTopicArn`: The Amazon Resource Name (ARN) that identifies the topic.
- `NotificationTopicStatus`: The current state of the topic.
- `ParameterGroupName`: The name of a parameter group for this cluster.
- `PreferredMaintenanceWindow`: A range of time when maintenance of DAX cluster software will be performed. For example: sun:01:00-sun:09:00. Cluster maintenance normally takes less than 30 minutes, and is performed automatically within the maintenance window.
- `SecurityGroupIds`: A list of user-specified security group IDs to be assigned to each node in the DAX cluster. If this parameter is not specified, DAX assigns the default VPC security group to each node.
"""
update_cluster(ClusterName; aws_config::AWSConfig=global_aws_config()) = dax("UpdateCluster", Dict{String, Any}("ClusterName"=>ClusterName); aws_config=aws_config)
update_cluster(ClusterName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("UpdateCluster", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClusterName"=>ClusterName), args)); aws_config=aws_config)

"""
    UpdateParameterGroup()

Modifies the parameters of a parameter group. You can modify up to 20 parameters in a single request by submitting a list parameter name and value pairs.

# Required Parameters
- `ParameterGroupName`: The name of the parameter group.
- `ParameterNameValues`: An array of name-value pairs for the parameters in the group. Each element in the array represents a single parameter.

"""
update_parameter_group(ParameterGroupName, ParameterNameValues; aws_config::AWSConfig=global_aws_config()) = dax("UpdateParameterGroup", Dict{String, Any}("ParameterGroupName"=>ParameterGroupName, "ParameterNameValues"=>ParameterNameValues); aws_config=aws_config)
update_parameter_group(ParameterGroupName, ParameterNameValues, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("UpdateParameterGroup", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ParameterGroupName"=>ParameterGroupName, "ParameterNameValues"=>ParameterNameValues), args)); aws_config=aws_config)

"""
    UpdateSubnetGroup()

Modifies an existing subnet group.

# Required Parameters
- `SubnetGroupName`: The name of the subnet group.

# Optional Parameters
- `Description`: A description of the subnet group.
- `SubnetIds`: A list of subnet IDs in the subnet group.
"""
update_subnet_group(SubnetGroupName; aws_config::AWSConfig=global_aws_config()) = dax("UpdateSubnetGroup", Dict{String, Any}("SubnetGroupName"=>SubnetGroupName); aws_config=aws_config)
update_subnet_group(SubnetGroupName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dax("UpdateSubnetGroup", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("SubnetGroupName"=>SubnetGroupName), args)); aws_config=aws_config)


text:
module Example
export hello, testing

"""
    hello(who::String)

Return "Hello, `who`".
"""
hello(who::String) = "Hello, $who"

testing() = "boink"

end


text:
using Test

# Category
##########

A, B = Ob(FreeCategory, :A), Ob(FreeCategory, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)

# Expression types
@test isa(A, FreeCategory.Ob) && isa(f, FreeCategory.Hom)
@test isa(A, ObExpr) && isa(f, HomExpr)
@test isa(A, CategoryExpr) && isa(f, CategoryExpr)

# Domains and codomains
@test dom(f) == A
@test codom(f) == B
@test dom(compose(f,g)) == A
@test codom(compose(f,g)) == A
@test_throws SyntaxDomainError compose(f,f)

# Associativity
@test compose(compose(f,g),f) == compose(f,compose(g,f))

# Extra syntax
@test compose(f,g,f) == compose(compose(f,g),f)
@test compose([f,g,f]) == compose(compose(f,g),f)
@test f⋅g == compose(f,g)
@test f⋅g⋅f⋅g == compose(f,g,f,g)

# String format
@test string(A) == "A"
@test string(f) == "f"
@test string(compose(f,g)) == "compose(f,g)"
@test string(compose(f,g,f)) == "compose(f,g,f)"
@test string(Ob(FreeCategory, nothing)) != ""

# S-expressions
@test sexpr(A) == ":A"
@test sexpr(f) == ":f"
@test sexpr(compose(f,g)) == "(compose :f :g)"
@test sexpr(compose(f,g,f)) == "(compose :f :g :f)"

# Infix notation (Unicode)
@test unicode(A) == "A"
@test unicode(f) == "f"
@test unicode(id(A)) == "id{A}"
@test unicode(compose(f,g)) == "f⋅g"

# Infix notation (LaTeX)
@test latex(A) == "A"
@test latex(f) == "f"
@test latex(id(A)) == "\\mathrm{id}_{A}"
@test latex(compose(f,g)) == "f \\cdot g"

@test latex(Ob(FreeCategory, "x")) == "x"
@test latex(Ob(FreeCategory, "sin")) == "\\mathrm{sin}"
@test latex(Ob(FreeCategory, "\\alpha")) == "\\alpha"

# 2-category
############

A, B, C, D = Ob(FreeCategory2, :A, :B, :C, :D)
f, g, F, G = [ Hom(sym, A, B) for sym in [:f,:g,:F,:G] ]
h, k, H, K = [ Hom(sym, B, C) for sym in [:h,:k,:H,:K] ]

# Domains and codomains
α, β = Hom2(:α, f, g), Hom2(:β, g, h)
@test dom(α) == f
@test codom(α) == g
@test dom(dom(α)) == A
@test codom(dom(α)) == B
@test dom(compose(α,β)) == f
@test codom(compose(α,β)) == h
@test_throws SyntaxDomainError compose2(α,β)

α, β = Hom2(:α, f, g), Hom2(:β, h, k)
@test dom(compose2(α,β)) == compose(f,h)
@test codom(compose2(α,β)) == compose(g,k)

# Infix notation (Unicode)
α, β = Hom2(:α, f, g), Hom2(:β, g, h)
@test unicode(compose(f,h)) == "f⋅h"
@test unicode(compose(α,β)) == "α⋅β"

α, β = Hom2(:α, f, g), Hom2(:β, h, k)
@test unicode(compose2(α,β)) == "α*β"

# Incompatible theories
@test_throws MethodError Hom(:f, Ob(FreeCategory, :A), Ob(FreeCategory2, :B))
f = Hom(:f, Ob(FreeCategory, :A), Ob(FreeCategory, :B))
g = Hom(:g, Ob(FreeCategory2, :B), Ob(FreeCategory2, :C))
@test_throws MethodError compose(f,g)


text:
using DoubleExponentialFormulas
using LinearAlgebra: norm
using Test

quadde32 = QuadDE(Float32)
quadde64 = QuadDE(Float64)
quaddeBF = QuadDE(BigFloat)

# Test Base.show works without error
let
    io = IOBuffer()
    @test (Base.show(io, MIME("text/plain"), quadde32); true)
    @test (Base.show(io, MIME("text/plain"), quadde64); true)
    @test (Base.show(io, MIME("text/plain"), quaddeBF); true)
end

# Abel transform
let
    f(r::AbstractFloat) = exp(-(r^2 - 1))
    expect = sqrt(BigFloat(π))

    y = 1.0
    g = r -> 2*f(r)*r/sqrt(r^2 - y^2)

    I, E = quadde32(g, 1, Inf)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(g, 1, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(g, 1, Inf, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end


text:
# ---
# title: 374. Guess Number Higher or Lower
# id: problem374
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Binary Search
# link: <https://leetcode.com/problems/guess-number-higher-or-lower/description/>
# hidden: true
# ---
#
# We are playing the Guess Game. The game is as follows:
#
# I pick a number from `1` to `n`. You have to guess which number I picked.
#
# Every time you guess wrong, I will tell you whether the number I picked is
# higher or lower than your guess.
#
# You call a pre-defined API `int guess(int num)`, which returns 3 possible
# results:
#
#   * `-1`: The number I picked is lower than your guess (i.e. `pick < num`).
#   * `1`: The number I picked is higher than your guess (i.e. `pick > num`).
#   * `0`: The number I picked is equal to your guess (i.e. `pick == num`).
#
# Return _the number that I picked_.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 10, pick = 6
#     Output: 6
#
#
# **Example 2:**
#
#
#
#     Input: n = 1, pick = 1
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: n = 2, pick = 1
#     Output: 1
#
#
# **Example 4:**
#
#
#
#     Input: n = 2, pick = 2
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 231 - 1`
#   * `1 <= pick <= n`
#
#
## @lc code=start
using LeetCode

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


text:
Base.include(Main, "opfdata.jl")
using JuMP
using Ipopt
using Printf
using SparseArrays
using DelimitedFiles

mutable struct SCOPFData
  opfdata::OPFData
  lines_off::Array
  #Float64::gener_ramp #generator ramp limit for contingency (percentage)
end

function scopf_solve(opfmodel, scopf_data)
  #
  # Initial point - needed especially for pegase cases
  #
  Pg0,Qg0,Vm0,Va0 = acopf_initialPt_IPOPT(scopf_data.opfdata)
  setvalue(getindex(opfmodel, :Pg), Pg0)
  setvalue(getindex(opfmodel, :Qg), Qg0)
  extra_jump=getindex(opfmodel, :extra)
  Vm_jump=getindex(opfmodel, :Vm)
  Va_jump=getindex(opfmodel, :Va)

  setvalue(Vm_jump[:,0], Vm0)
  setvalue(Va_jump[:,0], Va0)
  setvalue(extra_jump[:,0], 0.025*Pg0)
  ncont=length(scopf_data.lines_off); nbus=length(scopf_data.opfdata.buses)

  #println(opfmodel)

  for c in 1:ncont
    opfm1=opf_loaddata(ARGS[1], scopf_data.lines_off[c])
    Pg0,Qg0,Vm0,Va0 = acopf_initialPt_IPOPT(opfm1)
    setvalue(extra_jump[:,c], 0.025*Pg0)
    setvalue(Vm_jump[:,c], Vm0)
    setvalue(Va_jump[:,c], Va0)
  end

  status = solve(opfmodel)

  if status != :Optimal
    println("Could not solve the model to optimality.")
  end

  #scopf_outputAll(opfmodel, scopf_data)

  return opfmodel,status
end

function scopf_model(scopf_data)
  sd=scopf_data
  ncont=length(sd.lines_off); nbus=length(sd.opfdata.buses); ngen=length(sd.opfdata.generators)

  generators=sd.opfdata.generators; buses=sd.opfdata.buses; baseMVA = sd.opfdata.baseMVA

  opfmodel = Model(solver=IpoptSolver(print_level=5))

  println("Considering ", ncont, " contingengies")

  @variable(opfmodel, generators[i].Pmin <= Pg[i=1:ngen] <= generators[i].Pmax)
  @variable(opfmodel, 0<=extra[i=1:ngen,0:ncont]<=0.05*generators[i].Pmax)
  #@variable(opfmodel,Pg[i=1:ngen])
  @variable(opfmodel, generators[i].Qmin <= Qg[i=1:ngen] <= generators[i].Qmax)

  @NLobjective(opfmodel, Min, sum( generators[i].coeff[generators[i].n-2]*(baseMVA*(Pg[i]+extra[i,c]))^2
			             +generators[i].coeff[generators[i].n-1]*(baseMVA*(Pg[i]+extra[i,c]))
				     +generators[i].coeff[generators[i].n  ] for i=1:ngen, c=0:ncont) / (1+ncont))

  @variable(opfmodel, sd.opfdata.buses[i].Vmin <= Vm[i=1:nbus, c=0:ncont] <= sd.opfdata.buses[i].Vmax)
  @variable(opfmodel, Va[1:nbus, c=0:ncont])

  #fix the voltage angle at the reference bus
  for c in 0:ncont
    setlowerbound(Va[sd.opfdata.bus_ref,c], buses[sd.opfdata.bus_ref].Va)
    setupperbound(Va[sd.opfdata.bus_ref,c], buses[sd.opfdata.bus_ref].Va)
  end

  @constraint(opfmodel, ex[i=1:ngen,co=0:ncont], generators[i].Pmin <= Pg[i] + extra[i,co] <= generators[i].Pmax)
  for co=0:ncont
    #branch admitances
    if co==0
      YffR,YffI,YttR,YttI,YftR,YftI,YtfR,YtfI,YshR,YshI = computeAdmitances(sd.opfdata.lines, sd.opfdata.buses, sd.opfdata.baseMVA)
      lines=sd.opfdata.lines
      busIdx=sd.opfdata.BusIdx;FromLines=sd.opfdata.FromLines; ToLines=sd.opfdata.ToLines; BusGeners=sd.opfdata.BusGenerators
    else
      opfm1=opf_loaddata(ARGS[1], sd.lines_off[co])
      YffR,YffI,YttR,YttI,YftR,YftI,YtfR,YtfI,YshR,YshI = computeAdmitances(opfm1.lines, opfm1.buses, opfm1.baseMVA)
      lines=opfm1.lines
      busIdx = opfm1.BusIdx; FromLines = opfm1.FromLines; ToLines = opfm1.ToLines; BusGeners = opfm1.BusGenerators
    end
    nline=length(lines)

    # power flow balance
    for b in 1:nbus
      #real part
      @NLconstraint( opfmodel,
        ( sum( YffR[l] for l in FromLines[b]) + sum( YttR[l] for l in ToLines[b]) + YshR[b] ) * Vm[b,co]^2
        + sum(  Vm[b,co]*Vm[busIdx[lines[l].to],  co]*( YftR[l]*cos(Va[b,co]-Va[busIdx[lines[l].to],  co])
              + YftI[l]*sin(Va[b,co]-Va[busIdx[lines[l].to],co]  )) for l in FromLines[b] )
        + sum(  Vm[b,co]*Vm[busIdx[lines[l].from],co]*( YtfR[l]*cos(Va[b,co]-Va[busIdx[lines[l].from],co])
              + YtfI[l]*sin(Va[b,co]-Va[busIdx[lines[l].from],co])) for l in ToLines[b]   )
        -(sum(baseMVA*(Pg[g]+extra[g,co]) for g in BusGeners[b]) - buses[b].Pd) / baseMVA      # Sbus part
        ==0)

      #imaginary part
      @NLconstraint( opfmodel,
        ( sum(-YffI[l] for l in FromLines[b]) + sum(-YttI[l] for l in ToLines[b]) - YshI[b] ) * Vm[b,co]^2
        + sum(  Vm[b,co]*Vm[busIdx[lines[l].to],co]  *(-YftI[l]*cos(Va[b,co]-Va[busIdx[lines[l].to],  co])
              + YftR[l]*sin(Va[b,co]-Va[busIdx[lines[l].to],co]  )) for l in FromLines[b] )
        + sum( Vm[b,co]*Vm[busIdx[lines[l].from],co] *(-YtfI[l]*cos(Va[b,co]-Va[busIdx[lines[l].from],co])
              + YtfR[l]*sin(Va[b,co]-Va[busIdx[lines[l].from],co])) for l in ToLines[b]   )
        -(sum(baseMVA*Qg[g] for g in BusGeners[b]) - buses[b].Qd) / baseMVA      #Sbus part
        ==0)
    end # of for: power flow balance loop

    # branch/lines flow limits

    nlinelim=0
    for l in 1:nline
      if lines[l].rateA!=0 && lines[l].rateA<1.0e10
        nlinelim += 1
        flowmax=(lines[l].rateA/baseMVA)^2

        #branch apparent power limits (from bus)
        Yff_abs2=YffR[l]^2+YffI[l]^2; Yft_abs2=YftR[l]^2+YftI[l]^2
        Yre=YffR[l]*YftR[l]+YffI[l]*YftI[l]; Yim=-YffR[l]*YftI[l]+YffI[l]*YftR[l]
        @NLconstraint( opfmodel,
	  Vm[busIdx[lines[l].from],co]^2 *
	  ( Yff_abs2*Vm[busIdx[lines[l].from],co]^2 + Yft_abs2*Vm[busIdx[lines[l].to],co]^2
	  + 2*Vm[busIdx[lines[l].from],co]*Vm[busIdx[lines[l].to],co]
               *(Yre*cos(Va[busIdx[lines[l].from],co]-Va[busIdx[lines[l].to],co])-Yim*sin(Va[busIdx[lines[l].from],co]-Va[busIdx[lines[l].to],co]))
	  )
          - flowmax <=0)

        #branch apparent power limits (to bus)
        Ytf_abs2=YtfR[l]^2+YtfI[l]^2; Ytt_abs2=YttR[l]^2+YttI[l]^2
        Yre=YtfR[l]*YttR[l]+YtfI[l]*YttI[l]; Yim=-YtfR[l]*YttI[l]+YtfI[l]*YttR[l]
        @NLconstraint(
          opfmodel,
    	  Vm[busIdx[lines[l].to],co]^2 *
          ( Ytf_abs2*Vm[busIdx[lines[l].from],co]^2 + Ytt_abs2*Vm[busIdx[lines[l].to],co]^2
          + 2*Vm[busIdx[lines[l].from],co]*Vm[busIdx[lines[l].to],co]
               *(Yre*cos(Va[busIdx[lines[l].from],co]-Va[busIdx[lines[l].to],co])-Yim*sin(Va[busIdx[lines[l].from],co]-Va[busIdx[lines[l].to],co]))
          )
          - flowmax <=0)
      end
    end # of for: branch power limits

    @printf("Contingency %d -> Buses: %d  Lines: %d  Generators: %d\n", co, nbus, nline, ngen)
    println("     lines with limits:  ", nlinelim)

  end # of for: contingencies

  # minimize active power
#  @NLobjective(opfmodel,
#		  Min,
#		  sum( generators[i].coeff[generators[i].n] +
#		       sum(generators[i].coeff[generators[i].n-k]*(baseMVA*Pg[i])^k for k=1:generators[i].n-1)
#		       for i=1:ngen)
#		 )

  return opfmodel
end

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

  #values = zeros(2*nbus+2*ngen)
  ## values[1:2*nbus+2*ngen] = readdlm("/sandbox/petra/work/installs/matpower5.1/vars2.txt")
  #values[1:2*nbus+2*ngen] = readdlm("/sandbox/petra/work/installs/matpower5.1/vars3_9241.txt")
  #d = JuMP.NLPEvaluator(opfmodel)
  #MathProgBase.initialize(d, [:Jac])

  #g = zeros(2*nbus+2*nlinelim)
  #MathProgBase.eval_g(d, g, values)
  #println("f=", MathProgBase.eval_f(d,values))

  #gmat=zeros(2*nbus+2*nlinelim)
  #gmat[1:end] = readdlm("/sandbox/petra/work/installs/matpower5.1/cons3_9241.txt")
  #println("diff: ", norm(gmat-g))

  #println(opfmodel)

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

function scopf_outputAll(opfmodel, scopf_data)
  #shortcuts for compactness
  sd=scopf_data; opf_data=sd.opfdata
  lines = opf_data.lines; buses = opf_data.buses; generators = opf_data.generators; baseMVA = opf_data.baseMVA
  busIdx = opf_data.BusIdx; FromLines = opf_data.FromLines; ToLines = opf_data.ToLines; BusGeners = opf_data.BusGenerators;

  nbus  = length(buses); nline = length(lines); ngen  = length(generators)

  # OUTPUTING
  println("Objective value: ", getobjectivevalue(opfmodel), "USD/hr")
  VM=getvalue(getindex(opfmodel,:Vm)); VA=getvalue(getindex(opfmodel,:Va));
  PG=getvalue(getindex(opfmodel,:Pg)); QG=getvalue(getindex(opfmodel,:Qg));

  VM=VM[:,0]; VA=VA[:,0]; #base case

  EX=getvalue(getindex(opfmodel,:extra));
  EX=EX[:,0];

  # printing the first stage variables
  println("============================= BUSES ==================================")
  println("  Generator  |  extra ")   # |    P (MW)     Q (MVAr)")  #|         (load)   ")
  println("----------------------------------------------------------------------")
  for i in 1:ngen
      @printf("  %10d | %6.2f \n",generators[i].bus, EX[i])
  end
  println("\n")

  println("============================= BUSES ==================================")
  println("  BUS    Vm     Va   |   Pg (MW)    Qg(MVAr) ")   # |    P (MW)     Q (MVAr)")  #|         (load)   ")

  println("                     |     (generation)      ")
  println("----------------------------------------------------------------------")
  for i in 1:nbus
    @printf("%4d | %6.2f  %6.2f | %s  | \n",
	    buses[i].bus_i, VM[i], VA[i]*180/pi,
	    length(BusGeners[i])==0 ? "   --          --  " : @sprintf("%7.2f     %7.2f", baseMVA*PG[BusGeners[i][1]], baseMVA*QG[BusGeners[i][1]]))
  end
  println("\n")

  within=20 # percentage close to the limits

  nflowlim=0
  for l in 1:nline
    if lines[l].rateA!=0 && lines[l].rateA<1.0e10
      nflowlim += 1
    end
  end
  return

  if nflowlim>0
    println("Number of lines with flow limits: ", nflowlim)

    optvec=zeros(2*nbus+2*ngen)
    optvec[1:ngen]=PG
    optvec[ngen+1:2*ngen]=QG
    optvec[2*ngen+1:2*ngen+nbus]=VM
    optvec[2*ngen+nbus+1:2*ngen+2*nbus]=VA

    d = JuMP.NLPEvaluator(opfmodel)
    MathProgBase.initialize(d, [:Jac])

    consRhs = zeros(2*nbus+2*nflowlim)
    MathProgBase.eval_g(d, consRhs, optvec)

    #println(consRhs)

    @printf("================ Lines within %d %s of flow capacity ===================\n", within, "%")
    println("Line   From Bus    To Bus    At capacity")

    nlim=1
    for l in 1:nline
      if lines[l].rateA!=0 && lines[l].rateA<1.0e10
        flowmax=(lines[l].rateA/baseMVA)^2
        idx = 2*nbus+nlim

        if( (consRhs[idx]+flowmax)  >= (1-within/100)^2*flowmax )
          @printf("%3d      %3d      %3d        %5.3f%s\n", l, lines[l].from, lines[l].to, 100*sqrt((consRhs[idx]+flowmax)/flowmax), "%" )
          #@printf("%7.4f   %7.4f    %7.4f \n", consRhs[idx], consRhs[idx]+flowmax,  flowmax)
        end
        nlim += 1
      end
    end
  end

  return
end

# Compute initial point for IPOPT based on the values provided in the case data
function acopf_initialPt_IPOPT(opfdata)
  Pg=zeros(length(opfdata.generators)); Qg=zeros(length(opfdata.generators)); i=1
  for g in opfdata.generators
    # set the power levels in in between the bounds as suggested by matpower
    # (case data also contains initial values in .Pg and .Qg - not used with IPOPT)
    Pg[i]=0.5*(g.Pmax+g.Pmin)
    Qg[i]=0.5*(g.Qmax+g.Qmin)
    i=i+1
  end
  @assert i-1==length(opfdata.generators)

  Vm=zeros(length(opfdata.buses)); i=1;
  for b in opfdata.buses
    # set the ini val for voltage magnitude in between the bounds
    # (case data contains initials values in Vm - not used with IPOPT)
    Vm[i]=0.5*(b.Vmax+b.Vmin);
    i=i+1
  end
  @assert i-1==length(opfdata.buses)

  # set all angles to the angle of the reference bus
  Va = opfdata.buses[opfdata.bus_ref].Va * ones(length(opfdata.buses))

  return Pg,Qg,Vm,Va
end


text:
function leaky_fixed_rnn(activation, alpha, W, W_in, x, y)
    return (1-alpha).*x + alpha*activation.((W*x)+(W_in*y))
end

function double_leaky_fixed_rnn(alpha, f_activation, s_activation, f_lambda, s_lambda, W, W_in, x, y)
    return (1-alpha).*x + f_lambda*f_activation.((W*x)+(W_in*y)) + s_lambda*s_activation.((W*x)+(W_in*y))
end