text:
function irf(V::VAR, H::Int64, cholesky::Bool=true)
         K = size(V.Σ,1)
         cholesky == true ? Sigma = chol(V.Σ)' : Sigma = eye(K,K) # Cholesky or reduced form
         B0 = get_VAR1_rep(V)
         J = [eye(K,K) zeros(K,K*(V.p-1))]
         IRF = reshape((J*B0^0*J'*Sigma)',K^2,1)
         for i = 1:H
           IRF = [IRF reshape((J*B0^i*J'*Sigma)',K^2,1)]
         end
         return IRF
end

# Allow to specify a particular shock
function irf(VAR::VAR, H::Int64, shock::Int64, cholesky::Bool=true)
         K = size(V.Σ,1)
         cholesky == true ? Sigma = chol(V.Σ)' : Sigma = eye(K,K)
         abs(shock)>K && error("shock must be between 1 and $K")
         B0 = get_VAR1_rep(V)
         J = [eye(K,K) zeros(K,K*(V.p-1))]
         IRF = reshape((J*B0^0*J'*Sigma)',K^2,1) # before B0^0 why?
         for i = 1:H
             IRF = [IRF reshape((J*B0^i*J'*Sigma)',K^2,1)] #Cholesky here has also the intercept
         end
         IRFs = zeros(K,H+1)
         IRFs[:,:] = IRF[shock:K:(size(IRF,1)-K+shock),:]
         return IRFs
end


text:
"""
```
remove_presample!(Nt0, loglh, s_pred, P_pred, s_filt, P_filt)
```

Remove the first `Nt0` periods from all other input arguments and return.
"""
function remove_presample!(Nt0::Int, loglh::Vector{S},
                           s_pred::Matrix{S}, P_pred::Array{S, 3},
                           s_filt::Matrix{S}, P_filt::Array{S, 3};
                           outputs::Vector{Symbol} = [:loglh, :pred, :filt]) where {S<:Real}
    if Nt0 > 0
        if :loglh in outputs
            loglh  = loglh[(Nt0+1):end]
        end
        if :pred in outputs
            s_pred = s_pred[:,    (Nt0+1):end]
            P_pred = P_pred[:, :, (Nt0+1):end]
        end
        if :filt in outputs
            s_filt = s_filt[:,    (Nt0+1):end]
            P_filt = P_filt[:, :, (Nt0+1):end]
        end
    end
    return loglh, s_pred, P_pred, s_filt, P_filt
end

function remove_presample!(Nt0::Int, loglh::Vector{S}) where {S<:Real}
    if Nt0 > 0
        return loglh[(Nt0+1):end]
    end
    return loglh
end

"""
```
solve_discrete_lyapunov(A, B, max_it = 50)
```
Solves the discrete Lyapunov equation.

The problem is given by

```
    AXA' - X + B = 0
```

``X`` is computed by using a doubling algorithm. In particular, we iterate to
convergence on ``X_j`` with the following recursions for ``j = 1, 2, ...``
starting from ``X_0 = B, a_0 = A``:

```
    a_j = a_{j-1} a_{j-1}
    X_j = X_{j-1} + a_{j-1} X_{j-1} a_{j-1}'
```

This function is directly copied and pasted from the QuantEcon.jl package
and extended to handle Flux.Tracker through multiple dispatch.

##### Arguments

- `A::Matrix{Float64}` : An `n x n` matrix as described above.  We assume in order
  for  convergence that the eigenvalues of ``A`` have moduli bounded by unity
- `B::Matrix{Float64}` :  An `n x n` matrix as described above.  We assume in order
  for convergence that the eigenvalues of ``B`` have moduli bounded by unity
- `max_it::Int(50)` :  Maximum number of iterations

##### Returns

- `gamma1::Matrix{Float64}` Represents the value ``X``

"""
function solve_discrete_lyapunov(A::Matrix{S},
                                 B::Matrix{S},
                                 max_it::Int=50) where {S<:Real}
    # TODO: Implement Bartels-Stewardt
    n = size(A, 2)
    alpha0 = reshape([A;], n, n)
    gamma0 = reshape([B;], n, n)

    alpha1 = fill!(similar(alpha0), zero(eltype(alpha0)))
    gamma1 = fill!(similar(gamma0), zero(eltype(gamma0)))

    diff = 5
    n_its = 1

    while diff > 1e-15

        alpha1 = alpha0*alpha0
        gamma1 = gamma0 + alpha0*gamma0*alpha0'

        diff = maximum(abs, gamma1 - gamma0)
        alpha0 = alpha1
        gamma0 = gamma1

        n_its += 1

        if n_its > max_it
            error("Exceeded maximum iterations, check input matrices")
        end
    end

    return gamma1
end

function solve_discrete_lyapunov(A::TrackedArray{S},
                                 B::TrackedArray{S},
                                 max_it::Int=50) where {S<:Real}
    # TODO: Implement Bartels-Stewardt
    n = size(A, 2)
    alpha0 = reshape([A;], n, n)
    gamma0 = reshape([B;], n, n)

    alpha1 = fill!(similar(alpha0), zero(return_tracker_parameter_type(alpha0)))
    gamma1 = fill!(similar(gamma0), zero(return_tracker_parameter_type(gamma0)))

    diff = 5
    n_its = 1

    while diff > 1e-15

        alpha1 = alpha0*alpha0
        gamma1 = gamma0 + alpha0*gamma0*alpha0'

        diff = maximum(abs, gamma1 - gamma0)
        alpha0 = alpha1
        gamma0 = gamma1

        n_its += 1

        if n_its > max_it
            error("Exceeded maximum iterations, check input matrices")
        end
    end

    return gamma1
end

"""
```
return_tracker_parameter_type(x)
```
returns the underlying type for a Tracker variable.

```jldoctest
julia> return_tracker_type(param(1.0))
Float64
```
"""
return_tracker_parameter_type(x::Tracker.TrackedReal{S}) where S<:Real = S
return_tracker_parameter_type(x::TrackedArray{S}) where S<:Real = S

ext: jl


text:
using Reactive

function standard_push_test(non_input::Signal)
    m = map(x->2x, non_input)
    pval = number()
    push!(non_input, pval)
    step()

    @fact value(non_input) --> pval
    @fact value(m) --> 2pval
end

facts("Push to non-input nodes") do

    a = Signal(number(); name="a")
    b = map(x -> x*x, a; name="b")

    context("push to map") do
        m = map(x->2x, b)

        push!(b, 10.0)
        step()
        @fact value(b) --> 10.0
        @fact value(m) --> 2value(b)

        push!(a, 2.0)
        step()
        @fact value(b) --> 4.0
        @fact value(m) --> 2value(b)

        push!(b, 3.0)
        step()
        @fact value(b) --> 3.0
        @fact value(m) --> 2value(b)
    end

    context("push to merge") do
        ## Merge
        d = Signal(number(); name="d")
        e = merge(b, d, a; name="e")
        m = map(x->2x, e)
        # precedence to d
        @fact value(e) --> value(d)
        @fact value(m) --> 2value(e)

        standard_push_test(e)
    end

    context("push to foldp") do
        gc()
        x = Signal(number())
        f = foldp(+, 0, x)

        standard_push_test(f)

        # Check for correct continuation
        pval = value(f)
        sval = number()
        push!(x, sval)
        step()

        @fact value(f) --> pval + sval
    end

    context("push to filter") do
        g = Signal(0)
        pred = x -> x % 2 != 0
        h = filter(pred, 1, g)
        standard_push_test(h)
    end

    context("push to sampleon") do
        # sampleon
        g = Signal(0)
        nv = number()
        push!(g, nv)
        step()
        i = Signal(true)
        j = sampleon(i, g)
        standard_push_test(j)
    end

    context("push to droprepeats") do
        # droprepeats
        k = Signal(1)
        l = droprepeats(k)

        standard_push_test(l)
    end

    context("push to filterwhen") do
        # filterwhen
        b = Signal(false)
        n = Signal(1)
        dw = filterwhen(b, 0, n)
        standard_push_test(dw)
    end

    context("push to previous") do
        x = Signal(0)
        y = previous(x)
        standard_push_test(y)
    end

    context("push to delay") do
        x = Signal(0)
        y = delay(x)
        standard_push_test(y)
    end

    context("bind non-input") do
        s = Signal(1; name="sig 1")
        m = map(x->2x, s; name="m")
        s2 = Signal(3; name="sig 2")
        push!(m, 10)
        step()
        @fact value(m) --> 10

        bind!(m, s2) #two-way bind
        @fact value(m) --> 3
        @fact value(s2) --> 3

        push!(m, 6)
        step()
        @fact value(m) --> 6
        @fact value(s2) --> 6

        push!(s2, 10)
        step()
        @fact value(m) --> 10
        @fact value(s2) --> 10
    end

end


text:
using TimeZones.TZData: LATEST_FILE, read_latest
using TimeZones.TZData: tzdata_latest_version, tzdata_versions

@testset "tzdata_versions" begin
    versions = tzdata_versions()
    @test first(versions) == "93g"  # Earliest release
    @test "2016j" in versions
end

@testset "tzdata_latest_version" begin
    @test occursin(r"^(?:\d{2}){1,2}[a-z]?$", tzdata_latest_version())

    # Validate the contents of the LATEST_FILE which will be automatically created when
    # downloading the latest data.
    @test isfile(LATEST_FILE)
    version, retrieved = read_latest(LATEST_FILE)
    @test occursin(r"\A(?:\d{2}){1,2}[a-z]?\z", version)
    @test isa(retrieved, DateTime)
end


text:
@doc raw"""
    struct GeneralObs{Tag,D,T,FPT,K,Tlo,Tg,Td} <: Observation{D,T}
        lin_obs::Tlo
        g::Tg
        dist::Td
        obs::T
        t::Float64
        θ::Vector{K}
    end

General observation of the underlying process that is of the form:
```math
g(x)+ξ,
```
where $ξ$ is distributed according to `dist`, $g$ (corresponding to passed `g`)
is a function specified externally, `obs` is the observation made at time `t`,
`θ` is a container that may contain parameters and `lin_obs` is the
approximation to this general observation scheme via the linearization with
Gaussian noise.

    GeneralObs(
        t,
        obs,
        linearization;
        dist=_default_dist(obs),
        g=identity,
        fpt=NoFirstPassageTimes(),
        Tag=0,
        θ=[],
    )

Base constructor.
"""
struct GeneralObs{Tag,D,T,FPT,K,Tlo,Tg,Td} <: Observation{D,T}
    lin_obs::Tlo
    g::Tg
    dist::Td
    obs::T
    t::Float64
    θ::Vector{K}

    function GeneralObs(
            t, obs::T, lin_obs::Tlo, fpt::FPT, g::Tg, dist::Td, θ::Vector{K}, Tag,
        ) where {T,Tlo,FPT,Tg,Td,K}
        D = length(obs)
        new{Tag,D,T,FPT,K,Tlo,Tg,Td}(lin_obs, g, dist, obs, t, θ)
    end
end

function GeneralObs(
        t,
        obs,
        linearization;
        dist=_default_dist(obs),
        g=identity,
        fpt=NoFirstPassageTimes(),
        Tag=0,
        θ=[],
    )
    D = length(obs)
    GeneralObs(t, obs, linearization, fpt, g, dist, θ, Tag)
end

_default_dist(obs) = GaussianDistribtuions(_default_μ(obs), _default_Σ(obs))

fpt_info(::GeneralObs{Tag,D,T,FPT}) where {Tag,D,T,FPT} = FPT
get_tag(::GeneralObs{Tag}) where Tag = Tag

function Base.summary(io::IO, obs::GeneralObs)
    v, t = ν(obs), obs.t
    θ = parameters(obs)
    parametric_type = get_tag(obs) > 0
    println(io, repeat("⏤", 40 ))
    println(
        io,
        "|Observation `v = g(x)+ξ`, where `g` is an operator defined by a ",
        "function ", typeof(obs.g), ", and `ξ` is a random variable given by",
        typeof(obs.dist), "."
    )
    println(io, "|...")
    println(io, "|| ν: $v (observation),\n|| \t→ typeof(ν): ", typeof(v),",")
    println(io, "|| made at time $t.")
    println(io, "|...")
    println(io, "|This is NOT an exact observation.")
    println(io, "|...")
    println(
        io,
        "|It ",
        parametric_type ?
        "depends on additional parameters, which are set to: $θ." :
        "does not depend on any additional parameters."
    )
    println(io, "|...")
    _summary_fpt(io, fpt_info(obs); prepend="|")
    println(io, "|...")
    println(
        io,
        "|To inspect the linearized approximation to this observation ",
        "scheme please type in:"
    )
    println(io, "|\tsummary(<name-of-the-variable>.lin_obs)")
    println(io, "|and hit ENTER.")
    println(io, repeat("⋆ ", 10))
end

function Base.rand(rng::Random.AbstractRNG, o::GeneralObs, x)
    error("Not Implemented")
end

Base.rand(o::GeneralObs, x) = rand(Random.GLOBAL_RNG, o, x)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdcard

"""
    'feffery'_antdcard(;kwargs...)
    'feffery'_antdcard(children::Any;kwargs...)
    'feffery'_antdcard(children_maker::Function;kwargs...)

An AntdCard component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): The content of the tab - will only be displayed if this tab is selected
- `id` (String; optional)
- `bodyStyle` (Dict; optional)
- `bordered` (Bool; optional)
- `className` (String; optional)
- `coverImg` (optional): . coverImg has the following type: lists containing elements 'src', 'alt', 'style'.
Those elements have the following types:
  - `src` (String; optional)
  - `alt` (String; optional)
  - `style` (Dict; optional)
- `extraLink` (optional): . extraLink has the following type: lists containing elements 'content', 'href', 'target', 'className', 'style'.
Those elements have the following types:
  - `content` (String; optional)
  - `href` (String; optional)
  - `target` (String; optional)
  - `className` (String; optional)
  - `style` (Dict; optional)
- `headStyle` (Dict; optional)
- `hoverable` (Bool; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `size` (a value equal to: 'default', 'small'; optional)
- `style` (Dict; optional)
- `title` (String; optional)
"""
function 'feffery'_antdcard(; kwargs...)
        available_props = Symbol[:children, :id, :bodyStyle, :bordered, :className, :coverImg, :extraLink, :headStyle, :hoverable, :loading_state, :size, :style, :title]
        wild_props = Symbol[]
        return Component("'feffery'_antdcard", "AntdCard", "feffery_antd_components", available_props, wild_props; kwargs...)
end

'feffery'_antdcard(children::Any; kwargs...) = 'feffery'_antdcard(;kwargs..., children = children)
'feffery'_antdcard(children_maker::Function; kwargs...) = 'feffery'_antdcard(children_maker(); kwargs...)

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
macro sym_str(def)
    esc(:(Symbol($(Meta.parse("\"$(escape_string(def))\"")))))
end

join_args(args) = join(args, ", ")
splitargs(args) = filter(!isempty, strip.(split(args, ",")))
subindex(arr::AbstractArray, indices) = [arr[i] for i ∈ indices]

function splitjoin(str, removed_parts; delim=" ")
    split_str = split(str, delim)
    kept_parts = filter(x -> x ∉ removed_parts, 1:length(split_str))
    join(subindex(split_str, kept_parts), delim)
end

splitstrip(str; delim=" ", stripped_chars=[' ']) = strip.(split(str, delim), Ref(stripped_chars))

macro return_value_if_key_exists(collection, type)
    :($(esc(type)) ∈ keys($(esc(collection))) && return $(esc(collection))[$(esc(type))])
end

"""
Automatically create the pass! argument as an anonymous function from the corresponding pass! method with the other struct fields as kwargs.
"""
macro preprocess_pass(expr)
    # s_name = hasfield(expr.args[2])
    s_name = expr.args[2] isa Symbol ? expr.args[2] : first(expr.args[2].args)
    expr_nfields = length(expr.args[3].args)
    fields = getindex.(Ref(expr.args[3].args), 2:2:expr_nfields)
    (isempty(fields) || last(fields) ≠ :pass!) && error("The last field of the struct must be a pass! field")
    nopass_args = filter(x -> x ≠ :pass!, fields)
    quote
        struct $(esc(expr.args[2]))
            $(esc.(fields)...)
            $(esc(s_name))($(esc.(nopass_args)...)) = new($(esc.(nopass_args)...), preprocess_pass($(esc(s_name)); Core.NamedTuple{$(Tuple(nopass_args))}(Tuple([$(esc.(nopass_args)...)]))...))
        end
    end
end

reverse_dict(d::T) where {T <: AbstractDict} = T(v => k for (k, v) ∈ d)


text:
"""
Convolutional block
"""

# Re-define leakyrelu function
# https://github.com/FluxML/Flux.jl/issues/963
myleakyrelu(x::Real, a = oftype(x / one(x), 0.01)) = max(a * x, x / one(x))

conv_block(in, out) = [
        Conv((3, 3), in => out; init = glorot_normal, pad = (1, 1)),
        BatchNorm(out),
        x->myleakyrelu.(x, 0.2f0)
    ]

function build_layers(n_layers, in_chs, conv_chs, out_chs, σ)
    layers = conv_block(in_chs, conv_chs)
    for _ in 1:n_layers - 2
        push!(layers, conv_block(conv_chs, conv_chs)...)
    end
    tail_layer = Conv((3, 3), conv_chs => out_chs, σ;
        init = glorot_normal, pad = (1, 1))
    push!(layers, tail_layer)
    return Chain(layers...)
end

"""
    DiscriminatorPyramid
"""
mutable struct DiscriminatorPyramid{T <: Tuple}
    chains::T
    DiscriminatorPyramid(xs...) = new{typeof(xs)}(xs)
end

build_single_discriminator(n_layers, conv_chs) = build_layers(n_layers, 3, conv_chs, 1, identity)

function DiscriminatorPyramid(n_stage::Integer, n_layers::Integer)
    ds = build_single_discriminator.(n_layers, channel_pyramid(n_stage))
    return DiscriminatorPyramid(gpu.(ds)...)
end

function DiscriminatorPyramid(image_shapes::Vector{Tuple{Int64,Int64}}, n_layers::Integer)
    DiscriminatorPyramid(Base.length(image_shapes), n_layers)
end

function Base.show(io::IO, d::DiscriminatorPyramid)
    print(io, "DiscriminatorPyramid(")
    join(io, d.chains, ", \n")
    print(io, ")")
end

"""
    NoiseConnection
"""
mutable struct NoiseConnection
    layers
    pad::Int64
end

@Flux.functor NoiseConnection

function (nc::NoiseConnection)(prev::T, noise::T) where {T <: AbstractArray{Float32,4}}
    pad = nc.pad
    raw_output = nc.layers(noise + prev)::T + prev
    return raw_output[1 + pad:end - pad, 1 + pad:end - pad, :, :]
end

function Base.show(io::IO, nc::NoiseConnection)
    print(io, "NoiseConnection(", nc.layers, ", ", nc.pad, ")")
end

"""
    GeneratorPyramid
"""
mutable struct GeneratorPyramid{T <: Tuple}
    image_shapes::Vector{Tuple{Int64,Int64}}
    noise_shapes::Vector{Tuple{Int64,Int64}}
    pad::Int64
    chains::T
    GeneratorPyramid(image_shapes, noise_shapes, pad, xs...) = new{typeof(xs)}(image_shapes, noise_shapes, pad, xs)
end

build_single_gen_layers(n_layers, conv_chs) = build_layers(n_layers, 3, conv_chs, 3, tanh)
build_single_generator(n_layers, conv_chs, pad) = NoiseConnection(build_single_gen_layers(n_layers, conv_chs), pad)

function GeneratorPyramid(image_shapes::Vector{Tuple{Int64,Int64}}, n_layers::Integer, pad::Integer = 5)
    n_stage = Base.length(image_shapes)
    # receptive field = 11, floor(11/2) = 5
    noise_shapes = [2 * pad .+ s for s in image_shapes]
    ds = build_single_generator.(n_layers, channel_pyramid(n_stage), pad)
    return GeneratorPyramid(image_shapes, noise_shapes, pad, gpu.(ds)...)
end

function Base.show(io::IO, d::GeneratorPyramid)
    print(io, "GeneratorPyramid(")
    print(io, d.image_shapes, ", ")
    print(io, d.noise_shapes, ", ")
    println(io, d.pad, ", ")
    join(io, d.chains, ", \n")
    print(io, ")")
end

function (genp::GeneratorPyramid)(xs::AbstractVector{T}, st::Integer, resize::Bool) where {T <: AbstractArray{Float32,4}}
    if st == 0
        zeros_shape = resize ? first(genp.noise_shapes) : first(genp.image_shapes)
        return zeros_like(T, expand_dim(zeros_shape...))
    end
    prev = genp(xs, st - 1, true)
    out = genp.chains[st](prev, xs[st])
    return resize ? resize_and_padding(out, genp.image_shapes[st + 1], genp.noise_shapes[st + 1]) : out
end

"""
    HyperParams
"""
mutable struct HyperParams
    scale::Float64                  # progression scale, > 1
    min_size_x::Int64               # minimal image width
    min_size_y::Int64               # minimal image height
    img_size_x::Int64               # output image width
    img_size_y::Int64               # output image height
    n_layers::Int64                 # number of conv layers
    max_epoch::Int64                # training epochs
    reduce_lr_epoch::Int64          # reduce learining rate after training `redule_lr_epoch` epochs
    save_image_every_epoch::Int64   # save generated image every `save_image_every_epoch` epoch
    save_loss_every_epoch::Int64    # save loss every `save_loss_every_epoch` epoch
    loop_dscr::Int64                # training steps par descriminator training epoch
    loop_gen::Int64                 # training steps par generator training epoch
    lr_dscr::Float64                # discriminator learining rate
    lr_gen::Float64                 # generator learning rate
    alpha::Float32                  # rec loss coefficient
    amplifier_init::Float32         # noise amplifier
    HyperParams() = new(4/3, 25, 25, 128, 128, 5, 2000, 1600, 500, 100, 3, 3, 5e-4, 5e-4, 50f0, 1f0)
end

show_dict(hp::HyperParams) = OrderedDict(string(nm) => getfield(hp, nm) for nm in fieldnames(HyperParams))
image_shapes(hp::HyperParams) = size_pyramid(hp.scale, (hp.min_size_x, hp.min_size_y), (hp.img_size_x, hp.img_size_y))

function setup_models(hp::HyperParams)
    img_shapes = image_shapes(hp)
    dscrp = DiscriminatorPyramid(img_shapes, hp.n_layers) |> gpu
    genp = GeneratorPyramid(img_shapes, hp.n_layers) |> gpu
    return dscrp, genp
end

max_stars_repo_name: matsueushi/SinGAN

max_issues_repo_name: matsueushi/SinGAN

max_forks_repo_name: matsueushi/SinGAN


text:
using DataStructures: counter

"""
    tokenfrequencies(corpus)

Return a counter (dict) of the token frequencies for `corpus`.
"""
function tokenfrequencies(corpus::Array{Vector{T}}) where T
    h = counter(T)
    for tokens in corpus
        push!(h, tokens)
    end
    h
end


text:
module MyExample

function helper()
    println("Hello!")
end

function main_func()
    helper()
end
export main_func

end


text:
# # Using Atomistic-Compatible Packages

# There are three main steps to using an implementation of the Atomistic API.
# This simple toy example will walk you through each step using [NBodySimulator.jl](https://github.com/SciML/NBodySimulator.jl) for the dynamics.
# Some of the values from this example were taken from [this guide](https://ase.tufts.edu/chemistry/lin/images/FortranMD_TeachersGuide.pdf).

# ## Step 0: Load Dependencies

using Atomistic
using AtomsBase
using NBodySimulator
using Unitful
using UnitfulAtomic

# ## Step 1A: Configuring the System
# First, we must create an [AtomsBase.jl](https://github.com/JuliaMolSim/AtomsBase.jl)-style system.
# Here we will create a small Argon cluster with periodic boundary conditions.
# We could use the AtomsBase constructors directly, but in this case we will use NBodySimulator helper functions to initialize the system.

N = 864
element = :Ar
box_size = 3.47786u"nm"
reference_temp = 94.4u"K"

initial_bodies = generate_bodies_in_cell_nodes(N, element, box_size, reference_temp)
initial_system = FlexibleSystem(initial_bodies, CubicPeriodicBoundaryConditions(austrip(box_size)))

# ## Step 1B: Configuring the Simulator
# Second, we initialize our simulator. In this case we are using an AndersenThermostat which is provided by NBodySimulator.

Δt = 1e-2u"ps"
steps = 2000
thermostat = AndersenThermostat(austrip(reference_temp), 0.1 / austrip(Δt))
simulator = NBSimulator(Δt, steps, thermostat = thermostat)

# ## Step 1C: Configuring the Potential
# Lastly, we specify the interatomic potential that we will use for the simulation, Lennard-Jones in this case.

potential = LennardJonesParameters(1.657e-21u"J", 0.34u"nm", 0.765u"nm")

# ## Step 2: Running the Simulation

result = simulate(initial_system, simulator, potential)

# ## Step 3: Analyzing the Results
# We can now analyze the simulation results.
# The Atomistic API exposes a variety of quantities from each timestep of the simulation (more details [here](@ref MolecularDynamicsResult_Specification)).
# In this example we will look at the temperature, energy, and radial distribution function.

# ### Temperature
plot_temperature(result, simulator.steps ÷ 200)

# ### Energy
plot_energy(result, simulator.steps ÷ 200)

# ### Radial Distribution Function (RDF)
plot_rdf(result, potential.σ, 0.5)


text:
using Documenter, Dashboards

makedocs(
    sitename = "Dashboards",
    format = Documenter.HTML(),
    modules = [Dashboards],
    pages = [
        "index.md",
        "API Docs" => "api.md"
    ]
)

deploydocs(
    repo = "github.com/waralex/Dashboards.jl.git",
    target = "build"
)


text:
@testset "basic building blocks" begin
    mt = UltimateMortality([0.5, 0.5])

    # whole life insurance
    ins = Insurance(
        SingleLife(mortality = mt, issue_age = 0),
        Yields.Constant(0.05)
    )

    @test all(survival(SingleLife(mortality = mt, issue_age = 0)) .== [1.0, 0.5, 0.25])

    @test [timepoints(ins)...] == [1.0, 2.0]
    @test [survival(ins)...] == [1.0, 0.5]
    @test [discount(ins)...] == [1.0 / 1.05, 1 / 1.05^2]
    @test benefit(ins) == 1.0
    @test [probability(ins)...] == [0.5, 0.25]
    @test [cashflows(ins)...] == [0.5, 0.25]
    @test [cashflows(ins)...] == benefit(ins) .* probability(ins)
    @test present_value(ins) ≈ 0.5 / 1.05 + 0.5 * 0.5 / 1.05^2
    @test present_value(ins) ≈ LifeContingencies.A(ins.life, ins.int)

    # basic life contingency tests
    @test survival(LifeContingency(ins), 0, 1) ≈ 0.5
    @test survival(LifeContingency(ins), 0, 2) ≈ 0.25
    @test discount(LifeContingency(ins), 0, 1) ≈ 1 / 1.05
    @test discount(LifeContingency(ins), 0, 2) ≈ 1 / 1.05^2

    # term life insurance
    ins = Insurance(
        SingleLife(mortality = mt, issue_age = 0),
        Yields.Constant(0.05),
        1
    )

    @test [survival(ins)...] == [1.0]
    @test [discount(ins)...] == [1.0 / 1.05]
    @test benefit(ins) == 1.0
    @test [probability(ins)...] == [0.5]
    @test [cashflows(ins)...] == [0.5]
    @test [cashflows(ins)...] == benefit(ins) .* probability(ins)
    @test [timepoints(ins)...] == [1.0]
    @test present_value(ins) ≈ 0.5 / 1.05

    # annuity due
    ins = AnnuityDue(
        SingleLife(mortality = mt, issue_age = 0),
        Yields.Constant(0.05)
    )

    @test [survival(ins)...] == [1.0, 0.5, 0.25]
    @test [discount(ins)...] == [1.0, 1 / 1.05^1, 1 / 1.05^2]
    @test benefit(ins) == 1.0
    @test [timepoints(ins)...] == [0.0, 1.0, 2.0]
    @test [probability(ins)...] == [1.0, 0.5, 0.25]
    @test [cashflows(ins)...] == [1.0, 0.5, 0.25]
    @test [cashflows(ins)...] == benefit(ins) .* probability(ins)
    @test present_value(ins) ≈ 1 + 1 * 0.5 / 1.05 + 1 * 0.25 / 1.05^2

    ins = AnnuityDue(
        SingleLife(mortality = mt, issue_age = 0),
        Yields.Constant(0.05),
        2
    )
    @test present_value(ins) ≈ 1 + 1 * 0.5 / 1.05

    ins = AnnuityImmediate(
        SingleLife(mortality = mt, issue_age = 0),
        Yields.Constant(0.05),
        1
    )

    @test present_value(ins) ≈ 1 * 0.5 * 1 / 1.05

    ins = AnnuityDue(
        SingleLife(mortality = mt, issue_age = 0),
        Yields.Constant(0.05),
        1,
        start_time = 1
    )

    @test timepoints(ins) == [1.0]
    @test [probability(ins)...] == [0.5]
    @test present_value(ins) ≈ 1 * 0.5 * 1 / 1.05

end

@testset "one  year" begin
    ins = LifeContingency(
        SingleLife(
            mortality = UltimateMortality([0.5]),
            issue_age = 0
        ),
        Yields.Constant(0.05)
    )

    @test survival(ins, 1) ≈ 0.5
    @test survival(ins, 0) ≈ 1.0
    @test survival(ins.life, 0, 0.5) ≈ 1 - 0.5 * 0.5
    @test survival(ins, 0.5) ≈ 1 - 0.5 * 0.5

    @test omega(ins) ≈ 1
    @test present_value(AnnuityDue(ins)) ≈ 1 + 1 * 0.5 / 1.05
    @test present_value(AnnuityDue(ins, 1)) ≈ 1
    @test present_value(AnnuityDue(ins, 0)) == 0
    @test present_value(AnnuityImmediate(ins, 1)) ≈ 1 * 0.5 / 1.05
    @test present_value(AnnuityImmediate(ins, 0)) == 0

    @test present_value(Insurance(ins)) ≈ 0.5 / 1.05
    @test present_value(Insurance(ins, 1)) ≈ 0.5 / 1.05
    @test present_value(Insurance(ins, 0)) ≈ 0

    ins_jl = LifeContingency(
        JointLife((ins.life, ins.life), LastSurvivor(), Frasier()),
        Yields.Constant(0.05)
    )

    @test omega(ins_jl) ≈ 1
    @test present_value(AnnuityDue(ins_jl)) ≈ 1 + 1 * 0.75 / 1.05
    @test present_value(AnnuityDue(ins_jl, 1)) ≈ 1
    @test present_value(AnnuityImmediate(ins_jl, 1)) ≈ 1 * 0.75 / 1.05
    @test present_value(AnnuityDue(ins_jl, 2)) ≈ 1 + 1 * 0.75 / 1.05
    @test present_value(AnnuityDue(ins_jl, 2; certain = 2)) ≈ 1 + 1 / 1.05
    @test present_value(AnnuityDue(ins_jl, 0)) == 0

    @test survival(ins_jl, 1) ≈ 0.5 + 0.5 - 0.5 * 0.5
    @test all(survival(JointLife((ins.life, ins.life), LastSurvivor(), Frasier())) .== [1.0, 0.75])
    @test present_value(Insurance(ins_jl)) ≈ 0.25 / 1.05
    @test present_value(Insurance(ins_jl, 1)) ≈ 0.25 / 1.05
    @test present_value(Insurance(ins_jl, 0)) ≈ 0
end

@testset "two year no discount" begin
    ins = LifeContingency(
        SingleLife(
            mortality = UltimateMortality([0.0, 0.0]),
            issue_age = 0
        ),
        Yields.Constant(0.00)
    )

    @test omega(ins) ≈ 2
    @test present_value(AnnuityDue(ins)) ≈ 3
    @test present_value(AnnuityDue(ins; start_time = 1)) ≈ 2
    @test present_value(AnnuityDue(ins, 1)) ≈ 1
    @test present_value(AnnuityDue(ins, 2)) ≈ 2
    @test_throws BoundsError present_value(AnnuityDue(ins, 2; start_time = 2))
    @test present_value(AnnuityDue(ins, 3)) ≈ 3
    @test present_value(AnnuityDue(ins, 0)) == 0
    @test present_value(AnnuityImmediate(ins, 0)) == 0
    @test present_value(AnnuityImmediate(ins, 1)) ≈ 1
    @test present_value(AnnuityImmediate(ins, 1; start_time = 1)) ≈ 1
    @test_throws BoundsError present_value(AnnuityImmediate(ins; start_time = 2))
    @test present_value(AnnuityImmediate(ins, 2)) ≈ 2
    @test present_value(AnnuityImmediate(ins)) ≈ 2
    @test present_value(AnnuityImmediate(ins; certain = 0)) ≈ 2
    @test present_value(AnnuityImmediate(ins; certain = 2)) ≈ 2

    @test present_value(Insurance(ins)) ≈ 0
    @test present_value(Insurance(ins, 1)) ≈ 0
    @test present_value(Insurance(ins, 0)) ≈ 0

    ins_jl = LifeContingency(
        JointLife(
            lives = (ins.life, ins.life),
            contingency = LastSurvivor(),
            joint_assumption = Frasier()),
        Yields.Constant(0.00)
    )

    @test survival(ins_jl, 0) ≈ 1.0
    @test survival(ins_jl, 1) ≈ 1
    @test survival(ins_jl, 2) ≈ 1

end

@testset "two year with interest" begin
    ins = LifeContingency(
        SingleLife(
            mortality = UltimateMortality([0.5, 0.5]),
            issue_age = 0
        ),
        Yields.Constant(0.05)
    )

    @test omega(ins) ≈ 2
    @test present_value(AnnuityDue(ins)) ≈ 1 + 1 * 0.5 * 1 / 1.05 + 1 * 0.25 / 1.05^2
    @test present_value(AnnuityDue(ins, 1)) ≈ 1
    @test present_value(AnnuityDue(ins, 2)) ≈ 1 + 1 * 0.5 * 1 / 1.05
    @test present_value(AnnuityDue(ins, 3)) ≈ 1 + 1 * 0.5 * 1 / 1.05 + 1 * 0.25 / 1.05^2
    @test present_value(AnnuityDue(ins, 0)) ≈ 0
    @test present_value(AnnuityImmediate(ins, 0)) ≈ 0
    @test present_value(AnnuityImmediate(ins, 1)) ≈ 1 * 0.5 * 1 / 1.05

    @test present_value(Insurance(ins)) ≈ 0.5 / 1.05 + 0.5 * 0.5 / 1.05^2
    @test present_value(Insurance(ins, 1)) ≈ 0.5 / 1.05
    @test present_value(Insurance(ins, 0)) ≈ 0

    ins_jl = LifeContingency(
        JointLife(
            lives = (ins.life, ins.life),
            contingency = LastSurvivor(),
            joint_assumption = Frasier()),
        Yields.Constant(0.05)
    )

    @test survival(ins_jl, 0) ≈ 1.0
    @test survival(ins_jl, 1) ≈ 0.5 + 0.5 - 0.5 * 0.5
    @test survival(ins_jl, 2) ≈ (0.25 + 0.25 - 0.25 * 0.25)

end

@testset "two years" begin
    ins = LifeContingency(
        SingleLife(
            mortality = UltimateMortality([0.5, 0.5]),
            issue_age = 0
        ),
        Yields.Constant(0.05)
    )

    @test omega(ins) ≈ 2
    @test present_value(AnnuityDue(ins, 1)) ≈ 1
    @test present_value(AnnuityDue(ins, 2)) ≈ 1 + 1 * 0.5 * 1 / 1.05
    @test present_value(AnnuityDue(ins, 3)) ≈ 1 + 1 * 0.5 * 1 / 1.05 + 1 * 0.25 * 1 / 1.05^2
    @test present_value(AnnuityDue(ins, 0)) == 0
    @test present_value(AnnuityImmediate(ins, 0)) == 0
    @test present_value(AnnuityImmediate(ins, 1)) ≈ 1 * 0.5 * 1 / 1.05

    @test present_value(Insurance(ins)) ≈ 0.5 / 1.05 + 0.5 * 0.5 / 1.05^2
    @test present_value(Insurance(ins, 1)) ≈ 0.5 / 1.05
    @test present_value(Insurance(ins, 0)) ≈ 0

    ins_jl = LifeContingency(
        JointLife(
            lives = (ins.life, ins.life),
            contingency = LastSurvivor(),
            joint_assumption = Frasier()),
        Yields.Constant(0.05)
    )

    @test omega(ins_jl) ≈ 2
    @test present_value(AnnuityDue(ins_jl, 1)) ≈ 1
    @test present_value(AnnuityDue(ins_jl, 2)) ≈ 1 + 1 * 0.75 / 1.05
    @test present_value(AnnuityDue(ins_jl, 3)) ≈ 1 + 1 * 0.75 / 1.05 + 1 * survival(ins_jl, 2) / 1.05^2
    @test present_value(AnnuityDue(ins_jl, 0)) == 0
    @test present_value(AnnuityImmediate(ins_jl, 0)) == 0
    @test present_value(AnnuityImmediate(ins_jl, 1)) ≈ 1 * 0.75 / 1.05

end

# assumes embedded 'testMort' table in combination with MortalityTables
t = UltimateMortality(maleMort)

@testset "demo mortality" begin
    i = Yields.Constant(0.05)
    ins = LifeContingency(SingleLife(mortality = t, issue_age = 0), i)

    @test l(ins, 0) ≈ 1.0
    @test l(ins, 1) ≈ 0.993010000000000
    @test l(ins, 2) ≈ 0.992566124530000

    @test D(ins, 0) ≈ 1.0
    @test D(ins, 1) ≈ 0.945723809523809
    @test D(ins, 2) ≈ 0.900286734267574

    @test N(ins, 0) ≈ 20.113017073119200
    @test N(ins, 1) ≈ 19.113017073119200
    @test N(ins, 2) ≈ 18.167293263595400

    @test C(ins, 0) ≈ 0.006657142857142910
    @test C(ins, 1) ≈ 0.000402608136054441
    @test C(ins, 2) ≈ 0.000258082197156658

    @test M(ins, 0) ≈ 0.0422372822324182
    @test M(ins, 1) ≈ 0.0355801393752753
    @test M(ins, 2) ≈ 0.0351775312392208

    @test present_value(Insurance(ins)) ≈ 0.04223728223
    @test present_value(Insurance(ins, 0)) ≈ 0.0
    @test present_value(Insurance(ins, 1)) ≈ 0.0066571428571429
    @test present_value(AnnuityDue(ins, 0)) ≈ 0.0
    @test present_value(AnnuityDue(ins, 1)) ≈ 1.0
    @test present_value(AnnuityDue(ins, 2)) ≈ 1.0 + survival(ins, 1) / 1.05
    @test present_value(AnnuityImmediate(ins, 0)) ≈ 0.0
    @test present_value(AnnuityImmediate(ins, 1)) ≈ survival(ins, 1) / 1.05

    @test present_value(Insurance(ins, 30)) ≈ 0.0137761089686975

    @test N(ins, 26) ≈ 5.156762988852310
    @test D(ins, 26) ≈ 0.275358702015970
    @test present_value(AnnuityDue(ins, 26)) ≈ 14.9562540842669

end


text:
"""
K-D Tree index type for storing vectors. It is a wrapper
around a `KDTree` NN structure and performs a more efficient
search using a distance-based similarity between vectors.
"""
struct KDTreeIndex{A,D} <: AbstractIndex
    tree::KDTree{A,D}
end

KDTreeIndex(data, args...; kwargs...) = KDTreeIndex(KDTree(densify(data), args...; kwargs...))

# Nearest neighbor search method
function knn_search(index::KDTreeIndex{A,D},
                    point::AbstractVector,
                    k::Int,
                    keep::Vector{Int}=collect(1:length(index));
                    kwargs...
                   ) where {A<:AbstractArray, D<:Metric}
    # Uses Euclidean distance by default
    _k = min(k, length(keep))
    skip = idx->!in(idx, keep)
    idxs, scores = knn(index.tree, Vector(point), _k, true, skip)
    return idxs, scores
end

# Length method
length(index::KDTreeIndex) = length(index.tree.data)


text:
@testset "Power Residue" begin
    R = ResidueRing(ZZ,7)
    Rx,x = R["x"]
    b = x^4+x+1 #irreducible
    a = Base.rand(Rx, 10)
    a = div(a, gcd(a, b))
    d = fmpz(3)
    @test power_residue(b^2,b,d,fmpz(7)) == 0
    @test power_residue(a,b,d,fmpz(7)) == mod(a^(div(7^degree(b)-1, 3)), b)
end


text:
struct FINITO_LFinito_iterable{R<:Real,C<:RealOrComplex{R},Tx<:AbstractArray{C},Tf,Tg}
    F::Array{Tf}            # smooth term
    g::Tg                   # nonsmooth term
    x0::Tx                  # initial point
    N::Int                    # of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}}  # Lipschitz moduli of nabla f_i
    γ::Maybe{Union{Array{R},R}}  # stepsizes
    sweeping::Int8             # to only use one stepsize γ
    batch::Int                # batch size
    α::R                    # in (0, 1), e.g.: 0.99
end

mutable struct FINITO_LFinito_state{R<:Real,Tx}
    γ::Array{R}             # stepsize parameter
    hat_γ::R                # average γ
    av::Tx                  # the running average
    ind::Array{Array{Int}}  # running index set
    d::Int                  # number of batches
    # some extra placeholders
    z::Tx
    ∇f_temp::Tx             # placeholder for gradients
    z_full::Tx
    inds::Array{Int}        # needed for shuffled only!
end

function FINITO_LFinito_state(γ::Array{R}, hat_γ::R, av::Tx, ind, d) where {R,Tx}
    return FINITO_LFinito_state{R,Tx}(
        γ,
        hat_γ,
        av,
        ind,
        d,
        copy(av),
        copy(av),
        copy(av),
        collect(1:d),
    )
end

function Base.iterate(iter::FINITO_LFinito_iterable{R}) where {R}
    N = iter.N
    r = iter.batch # batch size
    # create index sets
    ind = Vector{Vector{Int}}(undef, 0)
    d = Int(floor(N / r))
    for i = 1:d
        push!(ind, collect(r*(i-1)+1:i*r))
    end
    r * d < N && push!(ind, collect(r*d+1:N))
    # updating the stepsize
    if iter.γ === nothing
        if iter.L === nothing
            @warn "--> smoothness parameter absent"
            return nothing
        else
            γ = zeros(R, N)
            for i = 1:N
                isa(iter.L, R) ? (γ = fill(iter.α * R(iter.N) / iter.L, (N,))) :
                (γ[i] = iter.α * R(N) / (iter.L[i]))
            end
        end
    else
        isa(iter.γ, R) ? (γ = fill(iter.γ, (N,))) : (γ = iter.γ) # provided γ
    end
    #initializing the vectors
    hat_γ = 1 / sum(1 ./ γ)
    av = copy(iter.x0)
    for i = 1:N
        ∇f, ~ = gradient(iter.F[i], iter.x0)
        ∇f .*= hat_γ / N
        av .-= ∇f
    end
    state = FINITO_LFinito_state(γ, hat_γ, av, ind, cld(N, r))

    return state, state
end

function Base.iterate(
    iter::FINITO_LFinito_iterable{R},
    state::FINITO_LFinito_state{R},
) where {R}
    # full update
    prox!(state.z_full, iter.g, state.av, state.hat_γ)
    state.av .= state.z_full
    for i = 1:iter.N
        gradient!(state.∇f_temp, iter.F[i], state.z_full) # update the gradient
        state.av .-= (state.hat_γ / iter.N) .* state.∇f_temp
    end
    iter.sweeping == 3 && (state.inds = randperm(state.d)) # shuffled

    for j in state.inds
        prox!(state.z, iter.g, state.av, state.hat_γ)
        for i in state.ind[j]
            gradient!(state.∇f_temp, iter.F[i], state.z_full) # update the gradient
            state.av .+= (state.hat_γ / iter.N) .* state.∇f_temp
            gradient!(state.∇f_temp, iter.F[i], state.z) # update the gradient
            state.av .-= (state.hat_γ / iter.N) .* state.∇f_temp
            state.av .+= (state.hat_γ / state.γ[i]) .* (state.z .- state.z_full)
        end
    end

    return state, state
end

solution(state::FINITO_LFinito_state) = state.z

max_issues_repo_name: kul-forbes/CIAOAlgorithms

max_forks_repo_name: kul-forbes/CIAOAlgorithms


text:
export randomhomotopy, randomsystem

"""
    randomhomotopy(::Type{AbstractPolynomialHomotopy{T}}, size::Int; kwargs...)

Create a total degree homotopy where the target system is a [`randomsystem(T, size, size; kwargs...)`](@ref).

## Example
```julia-repl
julia> H, solutions = randomhomotopy(StraightLineHomotopy{Complex128}, 2, mindegree=3, maxdegree=6);
julia> length(H)
3
julia> nvariables(H)
3
```
"""
function randomhomotopy(H::Type{<:AbstractPolynomialHomotopy{T}}, size::Int; kwargs...) where {T<:Complex}
    F = randomsystem(T, size, size; kwargs...)
    totaldegree(H, F)
end
function randomhomotopy(H::Type{<:AbstractPolynomialHomotopy}, size::Int; kwargs...)
    F = randomsystem(size, size; kwargs...)
    totaldegree(H, F)
end

"""
    randomsystem([T=Complex128,] nequations::Int, nvars::Int; mindegree=0, maxdegree=5, rng=Base.Random.GLOBAL_RNG, density=rand() * 0.8 + 0.1)

Creates a random polynomial system of `nequations` equations with `nvars` variables (named ``x_1``, ...``x_{nvars}``).
Each polynomial has a total degree uniformly drawn from ``[mindegree, maxdegree]``.
The coefficients are drawn independently from the given `rng`.
With `density` you can control how many coefficients are non-zero. A value of `1.0` creates
a dense polynomial (i.e. every coefficient is non-zero). A value of `0.5` creates a polynomial
where only half of all monomials are non zero.

    randomsystem([T=Complex128,] degrees::Vector{Int}, variables::Vector{Symbol}; rng=N(0,1))

Create a random polynomial system with the given `degrees` and `variables`.
"""
function randomsystem(::Type{T}, nequations::Int, nvars::Int; mindegree=2, maxdegree=5, rng=Base.Random.GLOBAL_RNG, density=(rand() * 0.8 + 0.1)) where {T<:Number}
    @assert (0 < density ≤ 1) "Expexted a density between 0.0 and 1.0"
    degrees = rand(mindegree:maxdegree, nequations)
    vars = map(i -> Symbol("x_$(i)"), 1:nvars)
    randomsystem(T, degrees, vars, rng, density)
end
function randomsystem(nequations::Int, nvars::Int; kwargs...)
    randomsystem(Complex128, nequations, nvars; kwargs...)
end
function randomsystem(::Type{T}, degrees::Vector{Int}, vars::Vector{Symbol}, rng, density) where {T<:Number}
    map(degrees) do degree
        exponents = create_exponents(degree, length(vars))
        nexponents = length(exponents)
        indices = rand(1:nexponents, ceil(Int, nexponents * density))
        exp_matrix = Matrix{Int}(length(vars), length(indices))
        for (j, index) in enumerate(indices)
            exp_matrix[:,j] .= exponents[index]
        end
        coeffs = randomcoeffs(rng, T, length(indices))
        FP.Polynomial(exp_matrix, coeffs, vars)
    end
end

randomcoeffs(rng, ::Type{Complex{BigFloat}}, n) = convert.(Complex{BigFloat}, rand(rng, Complex128, n))
randomcoeffs(rng, ::Type{BigFloat}, n) = convert.(BigFloat, rand(rng, Float64, n))
randomcoeffs(rng, ::Type{T}, n) where {T<:Number} = rand(rng, T, n)

function randomsystem(degrees::Vector{Int}, vars::Vector{Symbol}; kwargs...)
    randomsystem(Complex128, degrees, vars; kwargs...)
end

function exponents_helper!(results, partial_result, curr_sum::Int, target_sum::Int, remaining_elements::Int)::Vector{Vector{Int}}
    if remaining_elements == 0
        push!(results, partial_result)
        return results
    end
    if curr_sum == target_sum
        for _ = 1:remaining_elements
            push!(partial_result, 0)
        end
        push!(results, partial_result)
        return results
    end
    if remaining_elements == 1
        for x = 0:(target_sum - curr_sum) - 1
            push!(results, [partial_result; x])
        end
        push!(partial_result, target_sum - curr_sum)
        push!(results, partial_result)
        return results
    end

    for x=0:(target_sum-curr_sum)
        exponents_helper!(results, [partial_result; x], curr_sum + x, target_sum, remaining_elements - 1)
    end
    results
end

create_exponents(total_degree, nvars) = exponents_helper!(Vector{Vector{Int}}(), Int[], 0, total_degree, nvars)


text:
# Sample postprocessing script for customized figure.
#
# Reading 2D simulation data, zooming-in to a region of interest,
# adding streamlines and contour at specific levels on top of colored mesh.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, PyPlot

file = "bulk.0000501.vlsv"
nameρ = "rho"
nameV = "rho_v"

meta = load(file)

pcolormesh(meta, nameρ)
streamplot(meta, nameV, comp="xy", color="w", density=2.0)

f, ax = plt.gcf(), plt.gca()
cbar = ax.collections[end].colorbar
boxcoords = [0, 20, -15, 15]
ax.set_xlim([boxcoords[1],boxcoords[2]])
ax.set_ylim([boxcoords[3],boxcoords[4]])

# Contour line at a specific level
pArgs = Vlasiator.set_args(meta, nameρ, RE)
x, y = Vlasiator.get_axis(pArgs)
data = Vlasiator.prep2d(meta, nameρ)'
CS = plt.contour(x, y, data, levels = [1e7],
                 colors=("k",),linestyles=("-",),linewidths=(0.5,))


text:
"""
Metabolite structure.

# Fields
````
id :: String
formula :: String
charge :: Int
compartment :: String
notes :: Dict{String, Vector{String}}
annotation :: Dict{String, Union{Vector{String}, String}}
````
"""
mutable struct Metabolite
    id::String
    formula::Maybe{String}
    charge::Maybe{Int}
    compartment::Maybe{String}
    notes::Notes
    annotations::Annotations

    Metabolite(
        id = "";
        formula = nothing,
        charge = nothing,
        compartment = nothing,
        notes = Notes(),
        annotations = Annotations(),
    ) = new(id, formula, charge, compartment, notes, annotations)
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

run MOI tests
=#

using Test
using Printf
import Hypatia
import Hypatia.Solvers
include(joinpath(@__DIR__, "moicones.jl"))
include(joinpath(@__DIR__, "moi.jl"))

@testset "MathOptInterface wrapper tests" begin

@testset "MOI cone tests" begin
    println("starting MOI wrapper cone tests")
    real_types = [
        Float64,
        # Float32,
        BigFloat,
        ]
    for T in real_types
        @testset "$T" begin
            println(T, " ...")
            test_moi_cones(T)
        end
    end
end

@testset "MOI.Test tests" begin
    println("\nstarting MOI.Test tests")
    options = [
        (Float64, Solvers.CombinedStepper, false),
        # (Float64, Solvers.CombinedStepper, true),
        # (Float32, Solvers.CombinedStepper, true),
        (BigFloat, Solvers.PredOrCentStepper, true),
        ]
    for (T, stepper, use_dense_model) in options
        default_options = (
            # verbose = true,
            verbose = false,
            default_tol_relax = 4,
            stepper = stepper{T}(),
            )
        test_info = "$T, $use_dense_model"
        @testset "$test_info" begin
            println(test_info, " ...")
            test_time = @elapsed test_moi(T,
                use_dense_model = use_dense_model; default_options...)
            @printf("%8.2e seconds\n", test_time)
        end
    end
end

end
;


text:
__precompile__()
module REPLGameOfLife

using Compat
export gol, gameoflife

include("../Terminal.jl/Terminal.jl")
using .Terminal

include("braille_hd.jl")
include("presets.jl")

function next_gen!(board::Array{Int,2}, survive::Vector{Int}, birth::Vector{Int})
    newboard = copy(board)
    dx, dy = size(board)
    for i in 1:dx, j in 1:dy
        live_cells = live_neighbors(board, i, j)
        if board[i,j] == 0
            live_cells ∈ birth && (newboard[i,j] = 1)
        else
            live_cells ∈ survive ? (newboard[i,j] += 1) : (newboard[i,j] = 0)
        end
    end
    board[:,:] = newboard
end

function live_neighbors(board, i, j)
    dy, dx = size(board)
    wrap(idx,range,size) = (idx + range + size - 1) % size + 1
    sum(board[wrap.(i, -1:1, dy), wrap.(j, -1:1, dx)] .> 0) - (board[i,j] > 0)
end

function gol(board = rand([1,0,0,0,0,0,0,0,0,0,0], 80, 80); pause=0.1, survive=[2,3], birth=[3])
    rawmode() do
        abort=[false]
        clear_screen()
        update_board!(arr2braille(zeros(board)), arr2braille(board))
        sleep(pause)
        generation = 0
        @async while !abort[1]
            clear_screen()
            old = copy(board)
            next_gen!(board, survive, birth)
            generation += 1
            update_board!(arr2braille(old), arr2braille(board))
            sleep(pause)
        end
        while !abort[1]
            c = readKey()
            c in ["Ctrl-C"]     && (abort[1]=true)
        end
    end
end

function gol(preset::Preset, size, args...; kwargs...)
    x, y = Int.(floor.(size ./2))
    dy, dx = preset.size .- 1
    board = zeros(Int, (2y,2x))
    board[y:y+dy, x:x+dx] = preset.data
    gol(board, args...; kwargs...)
end

gameoflife(args...; kwargs...) = gol(args...; kwargs...)

end # module


text:
# asy_expansions.jl: contains all asymptotic expansions, needing results from the precomputations
# About
#   Author       - Peter Opsomer (Peter.Opsomer@cs.kuleuven.be)
#   History      - Created October 2014, last edit October 2015

# Evaluate the asymptotic expansion for the generalised Jacobi polynomial in the lens-shaped region with a certain normalisation. Derivative is only computed when dw is specified, default is monic normalisation without derivatives.
# Input
#   n,z          - Degree and point at which to evaluate
#   alpha,beta,h - Parts of the weight
#   psi(x)       - Function for analytic phase function psi(x)
#   nrT          - Number of terms. 1 gives leading order term
#   Dinf         - Limit of the Szego function
#   Uright,Uleft - Right and left U-matrices
#   [nor         - Normalisation: 'm' for monic(default), 'o' for orthonormal]
#   [dh,dpsi     - Derivatives of part of weight and phase function]
# Output
#   ipi          - Asymptotic expansion of the polynomial in the lens
#   [dipi        - Derivative of polynomial in the lens]
function asy_lens(n,z,alpha,beta,h,psi,nrT,Dinf,Uright,Uleft,nor = "m",dh = false,dpsi=false)

if ismatch(r"o",nor)
    if (nrT == 1) || (length(Uright) == 0) 	nmz = 1;     else
	# In Julia, 1:(nrT-1) already is a column vector instead of a row
	nmz = sqrt(real(1+2*1im*Dinf^2*sum(reshape(Uright[2,1,1:(nrT-1),1] + Uleft[2,1,1:(nrT-1),1], (nrT-1), 1)./(n+1).^(1:(nrT-1) ) ) ) );
    end
end
function w(x) return (0im+1-x).^alpha.*(0im+1+x).^beta.*h(0im+x); end

RI = (1+0im)*eye(2);
for k = 1:nrT-1
    for m = 1:ceil(Integer, k/2)
	RI = RI + (Uright[:,:,k,m]./(z-1).^m + Uleft[:,:,k,m]./(z+1).^m)/n.^k;
    end
end
if nor == "m"
    ipi = ([2.0.^(1/2-n)./(sqrt(0im+w(z) ).*(1.0+z+0.0im).^(1/4).*(1.0-z+0.0im).^(1/4) )   0;]*RI*[Dinf*cos((n+1/2)*acos(z) +psi(z) -pi/4); -1im/Dinf*cos((n-1/2)*acos(z) +psi(z) -pi/4)])[1,1];
elseif nor == "o"
    ipi = ([(2/pi)^(1/2)./(sqrt(0im+w(z) ).*(1.0+z+0.0im).^(1/4).*(1.0-z+0.0im).^(1/4) )   0;]*RI*[cos((n+1/2)*acos(z) +psi(z) -pi/4)  ;  -1im/Dinf^2*cos((n-1/2)*acos(z) +psi(z) -pi/4)]*nmz)[1,1];
else
    error("Wrong normalisation")
end
if dh == false return ipi end

function dw(x) return -alpha*(1-x).^(alpha-1).*(1+x).^beta.*h(x)  +(1-x).^alpha.*beta*(1+x).^(beta-1).*h(x) +(1-x).^alpha.*(1+x).^beta.*dh(x); end

dRI = (1+1im)*zeros(2,2);
for k = 1:nrT-1
    for m = 1:ceil(Integer, k/2)
	dRI = dRI + (Uright[:,:,k,m]./(z-1).^(m+1)*(-m) + Uleft[:,:,k,m]./(z+1).^(m+1)*(-m))/n.^k;
    end
end

if nor == "m"
    dipi = (([2.0.^(1/2-n).*(dw(z)*(-1/2)./(sqrt(0im+w(z) ).^3.*(1.0+z+0.0im).^(1/4).*(1.0+0.0im-z).^(1/4) ) -1/4*(-2*z)./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(5/4).*(1.0+0.0im-z).^(5/4) ) )    0;]*RI + [2.0.^(1/2-n)./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )   0;]*dRI)*[Dinf*cos((n+1/2)*acos(z) +psi(z) -pi/4); -1im/Dinf*cos((n-1/2)*acos(z) +psi(z) -pi/4)]+ [2.0.^(1/2-n)./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )     0;]*RI*[-Dinf*sin((n+1/2)*acos(z) +psi(z) -pi/4)*(dpsi(z)-(n+1/2)./sqrt(0im+1.0-z)./sqrt(0im+1.0+z) ); 1im/Dinf*sin((n-1/2)*acos(z) +psi(z) -pi/4)*(dpsi(z)-(n-1/2)./sqrt(0im+1.0-z)./sqrt(0im+1.0+z) )])[1,1];
elseif nor == "o"
    dipi = ((([sqrt(0im+2/pi).*(dw(z)*(-1/2)./(sqrt(0im+w(z) ).^3.*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) ) -1/4*(-2*z)./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(5/4).*(1.0+0.0im-z).^(5/4) ) )    0;]*RI+ [sqrt(0im+2/pi)./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )   0;]*dRI)*[cos((n+1/2)*acos(z) +psi(z) -pi/4); -1im/Dinf^2*cos((n-1/2)*acos(z) +psi(z) -pi/4)] + [sqrt(0im+2/pi)./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )     0;]*RI* [-sin((n+1/2)*acos(z) +psi(z) -pi/4)*(dpsi(z)-(n+1/2)./sqrt(0im+1.0-z)./sqrt(0im+1.0+z) ); 1im/Dinf^2*sin((n-1/2)*acos(z) +psi(z) -pi/4)*(dpsi(z)-(n-1/2)./sqrt(0im+1.0-z)./sqrt(0im+1.0+z) )])*nmz)[1,1];
end
return (ipi, dipi)
end

# Evaluate the asymptotic expansion for the generalised Jacobi polynomial in the outer region with a certain normalisation. Derivative is only computed when dw is specified, default is monic normalisation without derivatives.
# Input
#   n,z          - Degree and point at which to evaluate
#   alpha,beta,h - Parts of the weight
#   psi(x)       - Function for analytic phase function psi(x)
#   nrT          - Number of terms. 1 gives leading order term
#   Dinf         - Limit of the Szego function
#   Uright,Uleft - Right and left U-matrices
#   [nor         - Normalisation: 'm' for monic(default), 'o' for orthonormal]
#   [dh,dpsi     - Derivatives of part of weight and phase function]
# Output
#   opi          - Asymptotic expansion of the outer polynomial
#   [dopi        - Derivative of outer polynomial]
function asy_outer(n,z,alpha,beta,h,psi,nrT,Dinf,Uright,Uleft,nor = "m",dh = false,dpsi=false)

if ismatch(r"o",nor)
    if (nrT == 1) || (length(Uright) == 0) nmz = 1; else
	nmz = sqrt(real(1+2*1im*Dinf^2*sum(reshape(Uright[2,1,1:(nrT-1),1] + Uleft[2,1,1:(nrT-1),1], (nrT-1), 1)./(n+1).^(1:(nrT-1) ) ) ) );
    end
end
function w(x) return (1-x).^alpha.*(1+x).^beta.*h(x); end

RI = (1+0im)*eye(2);
for k = 1:nrT-1
    for m = 1:ceil(Integer, k/2)
	RI = RI + (Uright[:,:,k,m]./(z-1).^m + Uleft[:,:,k,m]./(z+1).^m)/n.^k;
    end
end

if nor == "m"
    opi = ([2.0.^(-1/2-n+0.0)./(sqrt(0im+w(z)).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )   0;]*RI*[Dinf*exp((-1.0+0im)^(angle(z) <= 0)*1im*((n+1/2)*acos(z) +psi(z) -pi/4) );-1im/Dinf*exp((-1.0+0im)^(angle(z) <= 0)*1im*((n-1/2)*acos(z) +psi(z) -pi/4))])[1,1];
elseif nor == "o"
    opi = ([(2*pi).^(-1/2+0.0)./(sqrt(0im+w(z)).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )    0;]*RI*[exp((-1.0+0im)^(angle(z) <= 0)*1im*((n+1/2)*acos(z) +psi(z) -pi/4) );   -1im/Dinf^2*exp((-1.0+0im)^(angle(z) <= 0)*1im*((n-1/2)*acos(z) +psi(z) -pi/4))]*nmz)[1,1];
else
    error("Wrong normalisation")
end

if dh == false return opi end

function dw(x) return -alpha*(1-x).^(alpha-1).*(1+x).^beta.*h(x)  +(1-x).^alpha.*beta*(1+x).^(beta-1).*h(x) +(1-x).^alpha.*(1+x).^beta.*dh(x); end

dRI = (1+1im)*zeros(2,2);
for k = 1:nrT-1
    for m = 1:ceil(Integer, k/2)
	dRI = dRI + (Uright[:,:,k,m]./(z-1).^(m+1)*(-m) + Uleft[:,:,k,m]./(z+1).^(m+1)*(-m))/n.^k;
    end
end
if nor=="m"
    dopi = (([2.0.^(-1/2-n+0.0).*(dw(z)*(-1/2)./(sqrt(0im+w(z) ).^3.*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) ) -1/4*(-2*z)./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(5/4).*(1.0+0.0im-z).^(5/4) ) )    0;]*RI  + [2.0.^(0.0-1/2-n)./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )    0;]*dRI)* [Dinf*exp((-1.0+0im)^(angle(z) <= 0)*1im*((n+1/2)*acos(z) +psi(z) -pi/4) );  -1im/Dinf*exp((-1.0+0im)^(angle(z) <= 0)*1im*((n-1/2)*acos(z) +psi(z) -pi/4))]    + [2.0.^(-1/2-n+0.0)./(sqrt(0im+w(z)).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )    0;]*RI*[Dinf*exp((-1.0+0im)^(angle(z) <= 0)*1im*((n+1/2)*acos(z) +psi(z) -pi/4) )*(-1.0+0im)^(angle(z) <= 0)*1im*(dpsi(0.0im+z)-(n+1/2)./sqrt(0im+1.0-z)./sqrt(0im+1.0+z) );-1im/Dinf*exp((-1.0+0im)^(angle(z) <= 0)*1im*((n-1/2)*acos(z) +psi(z) -pi/4) )*(-1.0+0im)^(angle(z) <= 0)*1im*(dpsi(0.0im+z)-(n-1/2)./sqrt(0im+1.0-z)./sqrt(0im+1.0+z) )])[1,1];
elseif nor=="o"
    dopi = ((([1/sqrt(0im+2*pi).*(dw(z)*(-1/2)./(sqrt(0im+w(z) ).^3.*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) ) -1/4*(-2*z)./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(5/4).*(1.0+0.0im-z).^(5/4) ) )    0;]*RI    + [1/sqrt(0im+2*pi)./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )    0;]*dRI)*   [exp((-1.0+0im)^(angle(z) <= 0)*1im*((n+1/2)*acos(z) +psi(z) -pi/4) );-1im/Dinf^2*exp((-1.0+0im)^(angle(z) <= 0)*1im*((n-1/2)*acos(z) +psi(z) -pi/4))] + [1/sqrt(0im+2*pi)./(sqrt(0im+w(z)).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )     0;]*RI* [exp((-1.0+0im)^(angle(z) <= 0)*1im*((n+1/2)*acos(z) +psi(z) -pi/4) )*(-1.0+0im)^(angle(z) <= 0)*1im*(dpsi(0.0im+z)-(n+1/2)./sqrt(0im+1.0-z)./sqrt(0im+1.0+z) );-1im/Dinf^2*exp((-1.0+0im)^(angle(z) <= 0)*1im*((n-1/2)*acos(z) +psi(z) -pi/4) )* (-1.0+0im)^(angle(z) <= 0)*1im*(dpsi(0.0im+z)-(n-1/2)./sqrt(0im+1.0-z)./sqrt(0im+1.0+z) )])*nmz)[1,1];
end
return (opi, dopi)
end

#Evaluate the asymptotic expansion for the generalised Jacobi polynomial in the right boundary region with a certain method. Derivative is only computed when assigned, default is monic normalisation without derivatives or Q-s. When using the Q-s, a series of the cosine transform of the asymptotic expansion is also used and theta=arccos(z).
# Input
#   n,z          - Degree and point at which to evaluate.
#   alpha,beta,h - Parts of the weight
#   psi(x)       - Function for analytic phase function psi(x)
#   nrT          - Number of terms. 1 gives leading order term
#   Dinf         - Limit of the Szego function
#   Uright,Uleft - Right and left U-matrices
#   [method        - 'm' for monic(default), 'o' for orthonormal, 'mQ' for monic using the Q-s, 'oQ' for orthonormal with Q-s]
#   [dh,dpsi     - Derivatives of part of weight and phase function]
#   [Qright      - Right Q-matrices for when ismatch(r"Q",method)]
#   [theta       - theta=arccos(z) when computing with Q-s to save digits]
#   [contpsi	 - Contour integral in psi, needed for series expansions]
# Output
#   bpiR         - Asymptotic expansion of the right boundary polynomial
#   [dbpiR       - Derivative of right boundary polynomial]
function asy_right(n,z,alpha,beta,h,psi,nrT,Dinf,Uright,Uleft,method = "m",dh=false,dpsi=false,Qright=false,theta=false,contpsi=false)

if ismatch(r"o",method)
    if (nrT == 1) || (length(Uright) == 0) nmz = 1; else
	nmz = sqrt(real(1+2*1im*Dinf^2*sum(reshape(Uright[2,1,1:(nrT-1),1] + Uleft[2,1,1:(nrT-1),1], (nrT-1), 1)./(n+1).^(1:(nrT-1) ) ) ) );
    end
end

if ismatch(r"Q",method)
    if (abs(z-1) >1e-2) || (abs(theta) > 1e-2)
	warn("Theta is not small so using Q-s will probably give an inaccurate result")
    end
    mo = size(Qright,4)-1;
    sbes = 0+0im;	sden = 0+0im;	sthosth = 0+0im;	ssinpl = 0+0im;	ssinmi = 0+0im;	sdbes = 0+0im;
    for nn = 0:mo # Replaced factorial(n) by gamma(1.0+n) to avoid overflow
	sbes = sbes + (-1)^nn*(n/2)^(alpha+2*nn)*theta^(2*nn)/gamma(1+nn+alpha)/gamma(1.0+nn);
	sden = sden + (-1)^nn*theta^(2*nn)/gamma(3+2*nn);
	sthosth = sthosth + (-1)^nn*theta^(2*nn)/gamma(2.0+2*nn);
	ssinpl = ssinpl + (-1)^nn*(psi(z)+alpha*pi/2+theta/2).^(2*nn)/gamma(2.0+2*nn);
	ssinmi = ssinmi + (-1)^nn*(psi(z)+alpha*pi/2-theta/2).^(2*nn)/gamma(2.0+2*nn);
	sdbes = sdbes + (-1)^nn/gamma(1.0+nn)*(n/2)^(2*nn+alpha+1)*(4/n^2*theta^(2*nn)/gamma(nn+alpha)-theta^(2*nn+2)/gamma(2+nn+alpha))/2;
    end
    plcoef = (alpha+beta+1)/2 +sthosth/(4im*pi)*contpsi(z,1);
    micoef = (alpha+beta-1)/2 +sthosth/(4im*pi)*contpsi(z,1);
    sbesowei = sbes*(sden)^(-alpha/2+0.0);
    splsindbes = plcoef*ssinpl*sdbes*sden^(-alpha/2+0.0);
    smisindbes = micoef*ssinmi*sdbes*sden^(-alpha/2+0.0);

    RR = (1+1im)*zeros(2,2);
    for k = (nrT-1):-1:1 # Counting down to avoid roundoff error
	tmp = zeros(2,2);
	for i = 1:size(Qright,4)
	    tmp = tmp + Qright[:,:,k,i]*(-theta^2*sden)^(i-1);
	end
	RR = (RR + tmp)/n;
    end
    RR = RR + eye(2);
else
    function w(x) return (1-x).^alpha.*(1+x).^beta.*h(x); end
    if abs(z-1) < eps()^(1/3)
	warn("z is close to 1 so use Q-s and series expansions in 'method'")
    end
    function brac(k) convert(FloatingPoint, (k == 0)) + convert(FloatingPoint, (k != 0))*prod(4*alpha^2 -(2*(1:k)-1).^2)/(2.0^(2*k)*(gamma(1.0+k) ) ); end

    function DeltaR(k,x) brac(k-1)/2.0^k./(log(x + sqrt(0im+x-1).*sqrt(0im+x+1)+0im).^k)./(2*sqrt(0im+x+1.0)*sqrt(0im+x-1.0))*[Dinf 0; 0 Dinf^(-1+0.0)]*		[sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )    1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) ); 		-1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )    sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) ) ]*		[exp( 1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2) )    0  ;   0    exp(-1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2))]*		[ ((-1.0)^k)/k*(alpha^2+k/2-1/4)    -1im*(k-1/2)   ;   (-1.0)^k*1im*(k-1/2)    (alpha^2+k/2-1/4)/k]*		[exp(-1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2))   0 ;   0   exp(1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2) )]*		[sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )     -1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )   ;   	1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )   sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) ) ]*[Dinf^(-1+0.0) 0; 0 Dinf]; end

    s = (1+1im)*zeros(2,2,nrT-1);
    R = (1+1im)*zeros(2,2,nrT-1);
    RR = eye(2) + (1+1im)*zeros(2,2);
    for k = 1:nrT-1
	for m = 1:ceil(Integer, k/2)
	    R[:,:,k] = R[:,:,k] + Uright[:,:,k,m]./(z-1).^m + Uleft[:,:,k,m]./(z+1).^m;
	end
	s[:,:,k] = DeltaR(k,z);
	if mod(k,2)==0 # k is even: add extra matrix
	    s[:,:,k] = s[:,:,k] -brac(k-1)*(4*alpha^2+2*k-1)/2.0^(k+1)/k/(log(z+sqrt(0im+z-1).*sqrt(0im+z+1) +0im))^k*eye(2);
	end
	RR = RR + (R[:,:,k]-s[:,:,k] )/n^k; # Avoids defining R_0^{right} = I
	for j = 1:(k-1)
	    RR = RR - R[:,:,k-j]*s[:,:,j]/n^k;
	end
    end
    b1 = cos(psi(z) +alpha*pi/2 + acos(z)/2)*besselj(alpha,n*acos(z) ) + sin(psi(z) +alpha*pi/2 + acos(z)/2)*(besselj(alpha-1,n*acos(z) ) -       besselj(alpha+1,n*acos(z) ) )/2;
    b2 = cos(psi(z) +alpha*pi/2 -acos(z)/2)*   besselj(alpha,n*acos(z) ) + sin(psi(z) +alpha*pi/2 -acos(z)/2)*   (besselj(alpha-1,n*acos(z) ) - besselj(alpha+1,n*acos(z) ) )/2;
end

if method == "m"
    bpiR = ([sqrt(0im+pi*n*acos(z))./(2.0^n.*sqrt(0im+w(z) ).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )     0;]*RR*[Dinf*b1; -1im/Dinf*b2])[1,1];
elseif method == "o"
    bpiR = ([sqrt(0im+n*acos(z))./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )    0;]*RR*[b1; (-1im/Dinf^2*b2)]*nmz)[1,1];
elseif method == "mQ"
    bpiR = ([sqrt(0im+pi*n/sthosth)./(2.0.^n.*sqrt(0im+(1.0+0im+z).^beta.*h(z) ) )      0;]*RR*[Dinf*(cos(0im+psi(0im+z) +alpha*pi/2 + theta/2)*sbesowei +splsindbes); -1im/Dinf*(cos(0im+psi(0im+z) +alpha*pi/2 -theta/2)*sbesowei + smisindbes)])[1,1];
elseif method == "oQ"
    bpiR = ([sqrt(0im+n/sthosth)./(sqrt(0im+(1.0+0im+z).^beta.*h(z) ) )     0;]*RR*[(cos(0im+psi(0im+z) +alpha*pi/2 + theta/2)*sbesowei +splsindbes); -1im/Dinf^2*(cos(0im+psi(0im+z) +alpha*pi/2 -theta/2)*sbesowei + smisindbes)]*nmz)[1,1];

else
    error("Wrong method")
end

if dh == false return bpiR end

dRR = (1.0+1.0im)*zeros(2,2);
if ismatch(r"Q",method)
    dsbes = 0+0im;	dsden = 0+0im;	dsthosth = 0+0im;	dssinpl = 0+0im;	dssinmi = 0+0im;
    dsdbes = (n/2)^(alpha+1.0)*(-2.0)/gamma(2.0+0im+alpha)/2;
    for nn = 1:mo
	dsbes = dsbes + (-1.0)^nn*(n/2)^(alpha+2.0*nn)*theta^(2*nn-2.0)*2.0*nn/gamma(1.0+nn+alpha)/gamma(1.0+nn);
	dsden = dsden + (-1.0)^nn*theta^(2*nn-2.0)*2.0*nn/gamma(3.0+2*nn);
	dsthosth = dsthosth + (-1.0)^nn*theta^(2*nn-2.0)*2.0*nn/gamma(2.0+2*nn);
	dssinpl = dssinpl + (-1.0)^nn*(psi(0.0im+z)+alpha*pi/2+theta/2).^(2*nn-2.0)*2.0*nn/gamma(2.0+2*nn)*plcoef^2;
	dssinmi = dssinmi + (-1.0)^nn*(psi(0.0im+z)+alpha*pi/2-theta/2).^(2*nn-2.0)*2.0*nn/gamma(2.0+2*nn)*micoef^2;
	dsdbes = dsdbes + (-1.0)^nn/gamma(1.0+nn)*(n/2.0)^(2*nn+alpha+1.0)*(4.0/n^2.0*theta^(2*nn-2.0)*2*nn/gamma(nn+alpha+0.0)-theta^(2.0*nn)*(2*nn+2.0)/gamma(2.0+nn+alpha))/2;
    end
    dsbesowei = dsbes*(sden)^(-alpha/2+0.0) -alpha/2*sbes*sden^(-alpha/2-1.0)*dsden;
    dscosbesoweipl = cos(psi(0.0im+z) +alpha*pi/2 +theta/2)*dsbesowei - sbesowei*ssinpl*plcoef^2;
    dscosbesoweimi = cos(psi(0.0im+z) +alpha*pi/2 -theta/2)*dsbesowei - sbesowei*ssinmi*micoef^2;

    common1 = -alpha/2*sdbes*(sden)^(-alpha/2-1.0)*dsden + (sden)^(-alpha/2+0.0)*dsdbes;
    common2 = (dsthosth*contpsi(z+0im,1.0)+sthosth^2*contpsi(z+0im,2.0) )/(4im*pi)*sdbes*(sden)^(-alpha/2+0.0);

    dssindbesoweipl = plcoef*(common1*ssinpl + sdbes*(sden)^(-alpha/2+0.0)*dssinpl) + common2*ssinpl;
    dssindbesoweimi = micoef*(common1*ssinmi + sdbes*(sden)^(-alpha/2+0.0)*dssinmi) + common2*ssinmi;

    for k = (nrT-1):-1:1
	tmp = (1.0+1.0im)*zeros(2,2);
	for i = 2:size(Qright,4) #Start from i=2 to avoid constant term
	    tmp = tmp + Qright[:,:,k,i]*(-theta^2*sden)^(i-2.0)*(i-1.0)*(-sthosth);
	end
	dRR = (dRR + tmp)/n;
    end
else
    function dw(x) return -alpha*(1-x).^(alpha-1).*(1+x).^beta.*h(x)  +(1-x).^alpha.*beta*(1+x).^(beta-1).*h(x)+(1-x).^alpha.*(1+x).^beta.*dh(x); end
    function dSqrtPhi(x) 1/2./sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) ).*(1+x./sqrt(0im+x-1.0)./sqrt(0im+x+1.0) ); end

    function dDeltaR(k,x) DeltaR(k,x).*(log(x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0)+0im).^k).*2*sqrt(0im+x+1.0)*sqrt(0im+x-1.0).*( (-k)*log(x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0)+0im).^(-k-1.0)./(x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0)).*(1 +x./sqrt(0im+x-1.0)./sqrt(0im+x+1.0) )./(2*sqrt(0im+x+1.0)*sqrt(0im+x-1.0))	+1./(log(x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0)+0im).^k).*(-x/2).*(x+1.0+0im)^(-3/2+0.0).*(x-1.0+0im)^(-3/2+0.0) )    +     # End of derivative of term in front
		brac(k-1)/2.0^k./(log(x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0)+0im).^k)./(2*sqrt(0im+x+1.0)*sqrt(0im+x-1.0))*[Dinf 0; 0 Dinf^(-1+0.0)]*( [dSqrtPhi(x)    -1im*dSqrtPhi(x)./(x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )   ;    1im*dSqrtPhi(x)./(x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )    dSqrtPhi(x)]*[exp( 1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2) )   0  ;   0   exp(-1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2))]  +  [sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )     1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )   ;   -1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )     sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) ) ]*[exp( 1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2) )*1im*(-1.0+0im)^(angle(x-1) <= 0)*dpsi(x)    0   ;   0   exp(-1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2))*(-1im)*(-1.0+0im)^(angle(x-1) <= 0)*dpsi(x)])*	[ ((-1)^k)/k*(alpha^2+k/2-1/4)    -1im*(k-1/2)    ;   (-1.0)^k*1im*(k-1/2)     (alpha^2+k/2-1/4)/k]*	[exp(-1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2))   0  ;   0   exp(1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2) )]*[sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )     -1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )   ;   1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )   sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) ) ]*[Dinf^(-1+0.0) 0; 0 Dinf]   +  # End of derivative of Lambda, now derivative of Lambda^(-1)
		brac(k-1)/2.0^k./(log(x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0)+0im).^k)./(2*sqrt(0im+x+1.0)*sqrt(0im+x-1.0))*[Dinf 0; 0 Dinf^(-1+0.0)]*	[sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )     1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )  ;   -1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )    sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) ) ]*	[exp( 1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2) )   0  ;  0   exp(-1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2))]*[ ((-1.0)^k)/k*(alpha^2+k/2-1/4)    -1im*(k-1/2)   ;   (-1.0)^k*1im*(k-1/2)   (alpha^2+k/2-1/4)/k]*	(		[-1im*(-1.0+0im)^(angle(x-1) <= 0)*dpsi(x)*exp(-1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2))   0  ;   0   1im*(-1.0+0im)^(angle(x-1) <= 0)*dpsi(x)*exp(1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2) )]*		[sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )    -1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )   ;    	1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )    sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) ) ] + 	[exp(-1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2)) 0  ;   0   exp(1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) +alpha*pi/2) )]*	[dSqrtPhi(x)    1im*dSqrtPhi(x)./(x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )   ;    -1im*dSqrtPhi(x)./(x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )    dSqrtPhi(x)])*[Dinf^(-1+0.0) 0; 0 Dinf]; end

    ds = (1.0+1.0im)*zeros(2,2,nrT-1);
    dR = (1.0+1.0im)*zeros(2,2,nrT-1);

    for k = 1:nrT-1
	for m = 1:ceil(Integer, k/2)
	    dR[:,:,k] = dR[:,:,k] + Uright[:,:,k,m]./(z-1).^(m+1)*(-m) + Uleft[:,:,k,m]./(z+1).^(m+1)*(-m);
	end

	ds[:,:,k] = dDeltaR(k,z);
	if mod(k,2)==0 # k is even: add extra matrix
	    ds[:,:,k] = ds[:,:,k] -brac(k-1)*(4*alpha^2+2*k-1)/2.0^(k+1)/k*(-k)*log(z + sqrt(0im+z-1.0).*sqrt(0im+z+1.0)+0im).^(-k-1.0)./(z + sqrt(0im+z-1.0).*sqrt(0im+z+1.0)).*(1.0 +z./sqrt(0im+z-1.0)./sqrt(0im+z+1.0) )*eye(2);
	end
	dRR = dRR + (dR[:,:,k]-ds[:,:,k] )/n^k; # Avoids defining R_0^{right} = I
	for j = 1:(k-1)
	    dRR = dRR - (dR[:,:,k-j]*s[:,:,j]+R[:,:,k-j]*ds[:,:,j])/n^k;
	end
    end
    db1 = -sin(psi(z) +alpha*pi/2 + acos(z)/2)*(dpsi(z) -1/2/sqrt(0im+1.0-z)/sqrt(0im+1.0+z))*besselj(alpha,n*acos(z) ) +cos(psi(z) +alpha*pi/2 + acos(z)/2)*(besselj(alpha-1,n*acos(z) )    -besselj(alpha+1,n*acos(z) ) )/2*n*(-1.0)/sqrt(0im+1.0-z)/sqrt(0im+1.0+z) + # Above derivative of cos-term, below sine term
        cos(psi(z) +alpha*pi/2 + acos(z)/2)*(dpsi(z) -1/2/sqrt(0im+1.0-z)/sqrt(0im+1.0+z))*   (besselj(alpha-1,n*acos(z) ) - besselj(alpha+1,n*acos(z) ) )/2  +sin(psi(z) +alpha*pi/2 + acos(z)/2)*(besselj(alpha-2,n*acos(z) ) -2*besselj(alpha,n*acos(z) ) + besselj(alpha+2,n*acos(z) ) )/4*   n*(-1)/sqrt(0im+1.0-z)/sqrt(0im+1.0+z);
    db2 = -sin(psi(z) +alpha*pi/2 -acos(z)/2)* (dpsi(z) +1/2/sqrt(0im+1.0-z)/sqrt(0im+1.0+z) )*besselj(alpha,n*acos(z) )+cos(psi(z) +alpha*pi/2 -acos(z)/2)* (besselj(alpha-1,n*acos(z) ) - besselj(alpha+1,n*acos(z) ) )/2.0*(-n)./sqrt(0im+1.0-z)./sqrt(0im+1.0+z)  + cos(psi(z) +alpha*pi/2 -acos(z)/2)*(dpsi(z) +1/2/sqrt(0im+1.0-z)/sqrt(0im+1.0+z) )* (besselj(alpha-1,n*acos(z) ) - besselj(alpha+1,n*acos(z) ) )/2  + sin(psi(z) +alpha*pi/2 -acos(z)/2)*(besselj(alpha-2,n*acos(z) ) -2*besselj(alpha,n*acos(z) ) + besselj(alpha+2,n*acos(z) ) )/4*n*(-1)/sqrt(0im+1.0-z)/sqrt(0im+1.0+z) ;
end

if method == "m"
    dbpiR = (([sqrt(0im+pi*n)/2.0^n*(1/2/sqrt(0im+acos(z) )*(-1.0)./sqrt(0im+1.0-z)./sqrt(0im+1.0+z)./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) ) + sqrt(0im+acos(z) ).*(dw(z)*(-1/2)./(sqrt(0im+w(z) ).^3.*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) ) -1/4*(-2*z)./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(5/4).*(1.0+0.0im-z).^(5/4) ) ) )     0;]*RR + [sqrt(0im+pi*n*acos(z))./(2.0^n.*sqrt(0im+w(z) ).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )     0;]*dRR        )*[Dinf*b1  ;   -1im/Dinf*b2]  + # Above derivative of row vector and matrix, below of column vector
        [sqrt(0im+pi*n*acos(z))./(2.0^n.*sqrt(0im+w(z) ).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )    0;]*RR*[Dinf*db1;-1im/Dinf*db2])[1,1];
elseif method == "o"
    dbpiR = (( ([sqrt(0im+n)*(1/2/sqrt(0im+acos(z) )*(-1.0)./sqrt(0im+1.0-z)./sqrt(0im+1.0+z)./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) ) + sqrt(0im+acos(z) ).*(dw(z)*(-1/2)./(sqrt(0im+w(z) ).^3.*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) ) -1/4*(-2*z)./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(5/4).*(1.0+0.0im-z).^(5/4) ) ) )     0;]*RR + [sqrt(0im+n*acos(z))./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )     0;]*dRR)*[b1; -1im/Dinf^2*b2]  + #Above derivative of row vector and matrix, below of column vector
        [sqrt(0im+n*acos(z))./(sqrt(0im+w(z) ).*(1.0+0.0im+z).^(1/4).*(1.0+0.0im-z).^(1/4) )    0;]*RR*[db1; -1im/Dinf^2*db2] )*nmz )[1,1];
elseif method == "mQ"
    dbpiR = (([sqrt(0im+pi*n)/2.0^n*(-1/2*sthosth^(-3/2+0.0)*dsthosth./(sqrt((0.0im+1+z).^beta.*h(0.0im+z) ) ) -  (dh(0.0im+z)*(-1/2)./(sqrt(h(0.0im+z).^3.*(0.0im+1+z).^beta) ) -beta/2./(sqrt(h(0.0im+z).*(0.0im+1+z).^beta).*(0.0im+1+z) ) ) )     0;]*RR + [sqrt(0im+pi*n/sthosth)./(2.0.^n.*sqrt(h(0im+z).*(0im+1+z).^beta) )    0;]*dRR)*[Dinf*(cos(psi(0im+z) +alpha*pi/2 + theta/2)*sbesowei + splsindbes ); -1im/Dinf*(cos(psi(0im+z) +alpha*pi/2 -theta/2)* sbesowei + smisindbes)]  +
#Above derivative of row vector and matrix, below of column vector
        [sqrt(0im+pi*n/sthosth)./(2.0.^n.*sqrt(h(0im+z).*(0im+1+z).^beta) )      0;]*RR*[Dinf*(dscosbesoweipl +dssindbesoweipl);-1im/Dinf*(dscosbesoweimi+dssindbesoweimi)])[1,1];
	dbpiR = -dbpiR/sthosth;
elseif method == "oQ"
    dbpiR = (([sqrt(0im+n)*(-1/2*sthosth^(-3/2+0.0)*dsthosth./(sqrt((0.0im+1+z).^beta.*h(0.0im+z) ) ) -  (dh(0.0im+z)*(-1/2)./(sqrt(h(0.0im+z).^3.*(0.0im+1+z).^beta) ) -beta/2./(sqrt(h(0.0im+z).*(0.0im+1+z).^beta).*(0.0im+1+z) ) ) )     0;]*RR + [sqrt(0im+n/sthosth)./(sqrt(h(0im+z).*(0im+1+z).^beta) )    0;]*dRR)*[(cos(psi(0im+z) +alpha*pi/2 + theta/2)*sbesowei + splsindbes ); -1im/Dinf^2*(cos(psi(0im+z) +alpha*pi/2 -theta/2)* sbesowei + smisindbes)]  +
#Above derivative of row vector and matrix, below of column vector
        [sqrt(0im+n/sthosth)./(sqrt(h(0im+z).*(0im+1+z).^beta) )    0;]*RR*[(dscosbesoweipl +dssindbesoweipl);-1im/Dinf^2*(dscosbesoweimi+dssindbesoweimi)])[1,1];
	dbpiR = -dbpiR/sthosth*nmz;
end
return (bpiR, dbpiR)
end

# Evaluate the asymptotic expansion for the generalised Jacobi polynomial in the left boundary region with a certain method. Derivative is only computed when dw is specified, default is monic normalisation without Q-s or derivatives. When using the Q-s, a series of the cosine transform of the asymptotic expansion is also used and theta=arccos(-z).
# Input
#   n,z            - Degree and point at which to evaluate
#   alpha,beta,h   - Parts of the weight
#   psi(x)         - Function for analytic phase function psi(x)
#   nrT            - Number of terms. 1 gives leading order term
#   Dinf           - Limit of the Szego function
#   Uright,Uleft   - Right and left U-matrices
#   [method        - 'm' for monic(default), 'o' for orthonormal, 'mQ' for monic using the Q-s, 'oQ' for orthonormal with Q-s]
#   [dh,dpsi       - Derivatives of part of weight and phase function]
#   [Qleft         - Left Q-matrices for when regexp(method,'Q') == 1]
#   [theta	   - theta=arccos(-z) when computing with Q-s to save digits]
#   [contpsi	   - Contour integral in psi, needed for series expansions]
# Output
#   bpiL           - Asymptotic expansion of the left boundary polynomial
#   [dbpiL         - Derivative of left boundary polynomial]
function asy_left(n,z, alpha,beta,h,psi,nrT,Dinf,Uright,Uleft, method="m",dh=false, dpsi=false,Qleft=false,theta=false,contpsi=false)

if ismatch(r"o",method)
    if (nrT == 1) || (length(Uright) == 0) nmz = 1; else
	nmz = sqrt(real(1+2*1im*Dinf^2*sum(reshape(Uright[2,1,1:(nrT-1),1] + Uleft[2,1,1:(nrT-1),1], (nrT-1), 1)./(n+1).^(1:(nrT-1) ) ) ) );
    end
end

if ismatch(r"Q",method)
    if (abs(z+1) >1e-2) || (abs(theta) > 1e-2)
	warn("Theta is not small so using Q-s will probably give an inaccurate result")
    end
    mo = size(Qleft,4)-1;
    sbes = 0im+0;	sden = 0im+0;	sthosth = 0im+0;	ssinpl = 0im+0;	ssinmi = 0im+0;	sdbes = 0im+0;
    for nn = 0:mo
	sbes = sbes + (-1)^nn*(n/2)^(beta+2*nn)*theta^(2*nn)/gamma(1+nn+beta)/gamma(1.0+nn);
	sden = sden + (-1)^nn*theta^(2*nn)/gamma(3+2*nn);
	sthosth = sthosth + (-1)^nn*theta^(2*nn)/gamma(2.0+2*nn);

	ssinpl = ssinpl + (-1)^(nn+1)*(psi(0im+z) -beta*pi/2 -theta/2).^(2*nn)/gamma(2.0+2*nn);
	ssinmi = ssinmi + (-1)^nn*(psi(0im+z) -beta*pi/2 +theta/2).^(2*nn)/gamma(2.0+2*nn);
	sdbes = sdbes + (-1)^nn/gamma(1.0+nn)*(n/2)^(2*nn+beta+1)*(4/n^2*theta^(2*nn)/gamma(nn+beta)-theta^(2*nn+2)/gamma(2+nn+beta))/2;
    end
    plcoef = (-alpha-beta-1)/2 +sthosth/(4im*pi)*contpsi(0im+z,1);
    micoef = (-alpha-beta+1)/2 +sthosth/(4im*pi)*contpsi(0im+z,1);
    sbesowei = sbes*(sden)^(-beta/2+0.0);
    splsindbes = plcoef*ssinpl*sdbes*sden^(-beta/2+0.0);
    smisindbes = micoef*ssinmi*sdbes*sden^(-beta/2+0.0);
    RL = (1+0im)*zeros(2,2);
    for k = (nrT-1):-1:1 # Counting down to avoid roundoff error
	tmp = zeros(2,2);
	for i = 1:size(Qleft,4)
	    tmp = tmp + Qleft[:,:,k,i]*(-theta^2*sden)^(i-1);
	end
	RL = (RL + tmp)/n;
    end
    RL = RL + eye(2);
else
    function w(x) return (1-x).^alpha.*(1+x).^beta.*h(x); end
    if abs(z+1) < eps()^(1/3)
	warn("z is close to -1 so use Q-s and series expansions in 'method'")
    end
    function brac(k) convert(FloatingPoint, (k == 0)) + convert(FloatingPoint, (k != 0))*prod(4*beta^2 -(2*(1:k)-1).^2)/(2.0^(2*k)*gamma(1.0+k) ); end # = (beta,k)

    function DeltaL(k,x) brac(k-1)/2.0^k./(log(-x -sqrt(0im+x-1.0).*sqrt(0im+x+1.0)+0im).^k)./(2*sqrt(0im+x+1.0)*sqrt(0im+x-1.0))*[Dinf 0; 0 (Dinf+0.0)^(-1+0.0)]*[sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )      1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) ); -1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )      sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) ) ]*	[exp( 1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2) ) 0; 0 exp(-1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2))]*[ ((-1.0)^k)/k*(beta^2+k/2-1/4)     1im*(k-1/2)    ;    (-1.0+0im)^(k+1)*1im*(k-1/2)       (beta^2+k/2-1/4)/k]*	[exp(-1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2)) 0; 0 exp(1im*(-1.0+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2) )]*[sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )     -1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) ); 1im/sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) )     sqrt(0im+x + sqrt(0im+x-1.0).*sqrt(0im+x+1.0) ) ]*[(Dinf+0.0)^(-1+0.0) 0; 0 Dinf]; end

    s = (1.0+1.0im)*zeros(2,2,nrT-1);
    R = (1.0+1.0im)*zeros(2,2,nrT-1);
    RL = (1+0im)*eye(2);
    for k = 1:nrT-1
	s[:,:,k] = DeltaL(k,z);
	for m = 1:ceil(Integer, k/2)
	    R[:,:,k] = R[:,:,k] + Uright[:,:,k,m]./(z-1).^m + Uleft[:,:,k,m]./(z+1).^m;
	end
	if mod(k,2)==0 # k is even: add extra matrix
	    s[:,:,k] = s[:,:,k] -brac(k-1)*(4*beta^2+2*k-1)/2.0^(k+1)/k/(log(-z-sqrt(0im+z-1.0).*sqrt(0im+z+1.0) +0im) )^k*eye(2);
	end
	RL = RL + (R[:,:,k]-s[:,:,k] )/n^k;
	for j = 1:(k-1)
	    RL = RL - R[:,:,k-j]*s[:,:,j]/n^k;
	end
    end
    b1 = sin(psi(z) -beta*pi/2 + acos(z)/2)*besselj(beta,n*acos(-z) ) +   cos(psi(z) -beta*pi/2 + acos(z)/2)*(besselj(beta-1,n*acos(-z) ) - besselj(beta+1,n*acos(-z) ) )/2;
    b2 = sin(psi(z) -beta*pi/2 -acos(z)/2)*  besselj(beta,n*acos(-z) ) + cos(psi(z) -beta*pi/2 -acos(z)/2)*    (besselj(beta-1,n*acos(-z) ) - besselj(beta+1,n*acos(-z) ) )/2;
end

if method == "m"
    bpiL = ([sqrt(0im+pi*n*acos(-z))./((-2.0)^n.*sqrt(0im+w(z)).*(1.0+z+0.0im).^(1/4).*(1.0-z+0.0im).^(1/4) )      0;]*RL*[Dinf*b1; -1im/Dinf*b2])[1,1]; # *(-1.0+0.0im)^(abs(angle(z+1))>pi/2);
elseif method == "o"
    bpiL = ([sqrt(0im+n*acos(-z))./((-1.0)^n.*sqrt(0im+w(z)).*(1.0+z+0.0im).^(1/4).*(1.0-z+0.0im).^(1/4) )    0;]*RL*[b1; -1im/Dinf^2*b2]*nmz)[1,1]; # *(-1.0+0.0im)^(abs(angle(z+1))>pi/2);
elseif method == "mQ"
    bpiL = ([sqrt(0im+pi*n/sthosth)./((-2.0).^n*sqrt((0im+1-z).^alpha.*h(0im+z) ) )     0;]*RL*[Dinf*(sin(psi(0im+z) -beta*pi/2 + (pi-theta)/2)*sbesowei +splsindbes); -1im/Dinf*(sin(psi(z) -beta*pi/2 -(pi-theta)/2)*sbesowei + smisindbes)])[1,1]; # *(-1.0+0.0im)^(abs(angle(z+1))>pi/2);
elseif method == "oQ"
    bpiL = ([sqrt(0im+n/sthosth)./((-1).^n*sqrt((0im+1-z).^alpha.*h(0im+z) ) )     0;]*RL*[(sin(psi(0im+z) -beta*pi/2 + (pi-theta)/2)*sbesowei +splsindbes); -1im/Dinf^2*(sin(psi(z) -beta*pi/2 -(pi-theta)/2)*sbesowei + smisindbes)]*nmz)[1,1]; # *(-1.0+0.0im)^(abs(angle(z+1))>pi/2);
else
    error("Wrong method")
end

if dh == false return bpiL end

dRL = (1.0+1im)*zeros(2,2);
if ismatch(r"Q",method)
    dsbes = 0im+0;	dsden = 0im+0;	dsthosth = 0im+0;	dssinpl = 0im+0;	dssinmi = 0im+0;
    dsdbes = (n/2)^(beta+1)*(-2)/gamma(2.0+beta)/2;
    for nn = 1:mo
	dsbes = dsbes + (-1)^nn*(n/2)^(beta+2*nn)*theta^(2*nn-2)*2*nn/gamma(1.0+nn+beta)/gamma(1.0+nn);
	dsden = dsden + (-1.0)^nn*theta^(2.0*nn-2.0)*2*nn/gamma(3.0+2*nn);
	dsthosth = dsthosth + (-1)^nn*theta^(2.0*nn-2.0)*2.0*nn/gamma(2.0+2*nn);
	dssinpl = dssinpl + (-1)^(nn+1.0)*(psi(0im+z)-beta*pi/2 -theta/2).^(2*nn-2)*2*nn/gamma(2.0+2*nn)*plcoef^2;
	dssinmi = dssinmi + (-1.0)^nn*(psi(0im+z)-beta*pi/2 +theta/2).^(2*nn-2)*2*nn/gamma(2.0+2*nn)*micoef^2;
	dsdbes = dsdbes + (-1.0)^nn/gamma(1.0+nn)*(n/2.0)^(2*nn+beta+1.0)*(4.0/n^2.0*theta^(2*nn-2.0)*2*nn/gamma(nn+beta+0.0)-theta^(2.0*nn)*(2*nn+2.0)/gamma(2.0+nn+beta))/2;
    end
    dsbesowei = dsbes*(sden)^(-beta/2+0.0) -beta/2*sbes*sden^(-beta/2-1.0)*dsden;
    dscosbesoweipl = sin(psi(0im+z) -beta*pi/2 +(pi-theta)/2)*dsbesowei + sbesowei*ssinpl*plcoef^2;
    dscosbesoweimi = sin(psi(0im+z) -beta*pi/2 -(pi-theta)/2)*dsbesowei + sbesowei*ssinmi*micoef^2;

    common1 = -beta/2*sdbes*(sden)^(-beta/2-1.0)*dsden + (sden)^(-beta/2+0.0)*dsdbes;
    common2 = (dsthosth*contpsi(0im+z,1.0)+sthosth^2*contpsi(0im+z,2.0) )/(4im*pi)*sdbes*(sden)^(-beta/2+0.0);

    dssindbesoweipl = plcoef*(common1*ssinpl + sdbes*(sden)^(-beta/2+0.0)*dssinpl) + common2*ssinpl;
    dssindbesoweimi = micoef*(common1*ssinmi + sdbes*(sden)^(-beta/2+0.0)*dssinmi) + common2*ssinmi;
    for k = (nrT-1):-1:1
	tmp = zeros(2,2);
	for i = 2:size(Qleft,4)
	    tmp = tmp + Qleft[:,:,k,i]*(theta^2*sden)^(i-2.0)*(i-1.0)*sthosth;
	end
	dRL = (dRL + tmp)/n;
    end
else
    function dw(x) return -alpha*(1-x).^(alpha-1).*(1+x).^beta.*h(x)  +(1-x).^alpha.*beta*(1+x).^(beta-1).*h(x)+(1-x).^alpha.*(1+x).^beta.*dh(x);end
    function dSqrtPhi(x) 1/2./sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) ).*(1+x./sqrt(0im+x-1)./sqrt(0im+x+1) ); end

    function dDeltaL(k,x) DeltaL(k,x).*(log(-x -sqrt(0im+x-1).*sqrt(0im+x+1)+0im).^k).*2*sqrt(0im+x+1)*sqrt(0im+x-1).*(  (-k)*log(-x -sqrt(0im+x-1).*sqrt(0im+x+1)+0im).^(-k-1.0)./(-x -sqrt(0im+x-1).*sqrt(0im+x+1)).*(-1 -x./sqrt(0im+x-1)./sqrt(0im+x+1) )./(2*sqrt(0im+x+1)*sqrt(0im+x-1))+1./(log(-x -sqrt(0im+x-1).*sqrt(0im+x+1)+0im).^k).*(-x/2).*(x+1+0im)^(-3/2+0.0).*(x-1+0im)^(-3/2+0.0) )   +   # End of derivative of term in front
		brac(k-1)/2.0^k./(log(-x -sqrt(0im+x-1).*sqrt(0im+x+1)+0im).^k)./(2*sqrt(0im+x+1)*sqrt(0im+x-1))*[Dinf 0; 0 Dinf^(-1+0.0)]*( [dSqrtPhi(x)    -1im*dSqrtPhi(x)./(x + sqrt(0im+x-1).*sqrt(0im+x+1) );  1im*dSqrtPhi(x)./(x + sqrt(0im+x-1).*sqrt(0im+x+1) )     dSqrtPhi(x)]*	[exp( 1im*(-1+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2) ) 0; 0 exp(-1im*(-1+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2))]	+[sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) )     1im/sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) ); 	-1im/sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) )     sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) ) ]*	[exp( 1im*(-1+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2) )*1im*(-1+0im)^(angle(x-1) <= 0)*dpsi(x) 0; 0 exp(-1im*(-1+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2))*(-1im)*(-1+0im)^(angle(x-1) <= 0)*dpsi(x)] )*[ ((-1)^k)/k*(beta^2+k/2-1/4)       1im*(k-1/2) ; (-1+0im)^(k+1)*1im*(k-1/2)        (beta^2+k/2-1/4)/k]* [exp(-1im*(-1+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2)) 0; 0 exp(1im*(-1+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2) )]*[sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) )       -1im/sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) ); 		1im/sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) )      sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) ) ]*[Dinf^(-1.0) 0; 0 Dinf]	+   # End of derivative of Lambda, now derivative of Lambda^(-1)
		brac(k-1)/2.0^k./(log(-x -sqrt(0im+x-1).*sqrt(0im+x+1)+0im).^k)./(2*sqrt(0im+x+1)*sqrt(0im+x-1))*[Dinf 0; 0 Dinf^(-1.0)]*	[sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) )       1im/sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) ); 	-1im/sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) )      sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) ) ]*[exp( 1im*(-1+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2) ) 0; 0 exp(-1im*(-1+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2))]*	[ ((-1)^k)/k*(beta^2+k/2-1/4)       1im*(k-1/2) ; (-1)^(k+1)*1im*(k-1/2)        (beta^2+k/2-1/4)/k]*	(	[-1im*(-1+0im)^(angle(x-1) <= 0)*dpsi(x)*exp(-1im*(-1+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2)) 0; 	0 1im*(-1+0im)^(angle(x-1) <= 0)*dpsi(x)*exp(1im*(-1+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2) )]*	[sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) )       -1im/sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) ); 	1im/sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) )      sqrt(0im+x + sqrt(0im+x-1).*sqrt(0im+x+1) ) ] + 	[exp(-1im*(-1+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2)) 0; 0 exp(1im*(-1+0im)^(angle(x-1) <= 0)*(psi(x) -beta*pi/2) )]*	[dSqrtPhi(x)      1im*dSqrtPhi(x)./(x + sqrt(0im+x-1).*sqrt(0im+x+1) )    ;     -1im*dSqrtPhi(x)./(x + sqrt(0im+x-1).*sqrt(0im+x+1) )     dSqrtPhi(x)])*[Dinf^(-1.0) 0; 0 Dinf]; end

    ds = (1.0+1im)*zeros(2,2,nrT-1);
    dR = (1.0+1im)*zeros(2,2,nrT-1);
    for k = 1:nrT-1
	for m = 1:ceil(Integer, k/2)
	    dR[:,:,k] = dR[:,:,k] + Uright[:,:,k,m]./(z-1).^(m+1)*(-m) + Uleft[:,:,k,m]./(z+1).^(m+1)*(-m);
	end
	ds[:,:,k] = dDeltaL(k,z);
	if mod(k,2)==0 # k is even: add extra matrix
	    ds[:,:,k] = ds[:,:,k] -brac(k-1)*(4*beta^2+2*k-1)/2.0^(k+1)/k		*(-k)*log(-z -sqrt(0im+z-1).*sqrt(0im+z+1)+0im).^(-k-1.0)./(-z -sqrt(0im+z-1).*sqrt(0im+z+1))	.*(-1 -z./sqrt(0im+z-1)./sqrt(0im+z+1) )*eye(2);
	end
	dRL = dRL + (dR[:,:,k]-ds[:,:,k] )/n^k; # Avoids defining R_0^{right} = I
	for j = 1:(k-1)
	    dRL = dRL - (dR[:,:,k-j]*s[:,:,j]+R[:,:,k-j]*ds[:,:,j])/n^k;
	end
    end
    db1 = cos(psi(z) -beta*pi/2 + acos(z)/2)*(dpsi(z)-1/2/sqrt(0im+1-z)/sqrt(0im+1+z))*    besselj(beta,n*acos(-z) ) +    sin(psi(z) -beta*pi/2 + acos(z)/2)*(besselj(beta-1,n*acos(-z) )    -besselj(beta+1,n*acos(-z) ))/2*n*(+1)/sqrt(0im+1.0-z)/sqrt(0im+1.0+z)   +  # Above derivative of sine-term, below cos term
        -sin(psi(z) -beta*pi/2 + acos(z)/2)*(dpsi(z)-1/2/sqrt(0im+1.0-z)/sqrt(0im+1+z) )*   (besselj(beta-1,n*acos(-z) ) - besselj(beta+1,n*acos(-z) ) )/2  +cos(psi(z) -beta*pi/2 + acos(z)/2)*(besselj(beta-2,n*acos(-z) )  -2*besselj(beta,n*acos(-z) ) +besselj(beta+2,n*acos(-z) ) )/4*   n*(+1)/sqrt(0im+1.0-z)/sqrt(0im+1.0+z);
    db2 = cos(psi(z) -beta*pi/2 -acos(z)/2)*(dpsi(z)-1/2*(-1)/sqrt(0im+1.0-z)/sqrt(0im+1.0+z) )*    besselj(beta,n*acos(-z) )  + sin(psi(z) -beta*pi/2 -acos(z)/2)*(besselj(beta-1,n*acos(-z) )   -besselj(beta+1,n*acos(-z) ) )/2*n*(+1)/sqrt(0im+1.0-z)/sqrt(0im+1.0+z)   -sin(psi(z) -beta*pi/2 -acos(z)/2)*(dpsi(z)+1/2/sqrt(0im+1.0-z)/sqrt(0im+1.0+z))*   (besselj(beta-1,n*acos(-z) ) - besselj(beta+1,n*acos(-z) ) )/2 + cos(psi(z) -beta*pi/2 -acos(z)/2)*   (besselj(beta-2,n*acos(-z) ) -2*besselj(beta,n*acos(-z) )   +besselj(beta+2,n*acos(-z) ) )/4*n/sqrt(0im+1.0-z)/sqrt(0im+1.0+z);
end

if method == "m"
	dbpiL = (([sqrt(0im+pi*n)*(1/2./sqrt(0im+acos(-z) )./(sqrt(0im+1-z).*sqrt(0im+1+z).*(-2.0)^n.*sqrt(0im+w(z)).*	(1.0+z+0.0im).^(1/4).*(1.0-z+0.0im).^(1/4) ) 	+sqrt(0im+acos(-z))./(-2.0)^n*(dw(z)*(-1/2)./(sqrt(0im+w(z) ).^3.*(1.0+z+0.0im).^(1/4).*(1.0-z+0.0im).^(1/4) ) 	-1/4*(-2*z)./(sqrt(0im+w(z) ).*(1.0+z+0.0im).^(5/4).*(1.0-z+0.0im).^(5/4) ) ) )        0;]*RL +       [sqrt(0im+pi*n*acos(-z))./((-2.0)^n.*sqrt(0im+w(z)).*(1.0+z+0.0im).^(1/4).*(1.0-z+0.0im).^(1/4) )       0;]*dRL)*      [Dinf*b1; -1im/Dinf*b2]  +  # Above derivative of row vector and matrix, below of column vector
        [sqrt(0im+pi*n*acos(-z))./((-2.0)^n.*sqrt(0im+w(z)).*(1.0+z+0.0im).^(1/4).*(1.0-z+0.0im).^(1/4) )      0;]*RL*[Dinf*db1 ; -1im/Dinf*db2])[1,1]; # *(-1.0+0.0im)^(abs(angle(z+1))>pi/2);
elseif method == "o"
	dbpiL = ((([sqrt(0im+n)*(1/2./sqrt(0im+acos(-z) )./(sqrt(0im+1.0-z).*sqrt(0im+1.0+z).*(-1.0)^n.*sqrt(0im+w(z)).*	(1.0+z+0.0im).^(1/4).*(1.0-z+0.0im).^(1/4) ) +sqrt(0im+acos(-z))./(-1.0)^n*(dw(z)*(-1/2)./(sqrt(0im+w(z) ).^3.*(1.0+z+0.0im).^(1/4).*(1.0+0.0im-z).^(1/4) ) -1/4*(-2*z)./(sqrt(0im+w(z) ).*(1.0+z+0.0im).^(5/4).*(1.0-z+0.0im).^(5/4) ) ) )      0;]*RL +     [sqrt(0im+n*acos(-z))./((-1.0)^n.*sqrt(0im+w(z)).*(1.0+z+0.0im).^(1/4).*(1.0-z+0.0im).^(1/4) )       0;]*dRL)*    [b1; -1im/Dinf^2*b2]      +   # Above derivative of row vector and matrix, below of column vector
        [sqrt(0im+n*acos(-z))./((-1.0)^n.*sqrt(0im+w(z)).*(1.0+z+0.0im).^(1/4).*(1.0-z+0.0im).^(1/4) )      0;]*RL*[db1;  -1im/Dinf^2*db2])*nmz)[1,1]; # *(-1.0+0.0im)^(abs(angle(z+1))>pi/2);
elseif method == "mQ"
    dbpiL = (([sqrt(0im+pi*n).*(0im-2.0)^(-n+0.0)*(-1/2*(0.0im+sthosth)^(-3/2+0.0)*dsthosth./(sqrt(0im+(0im+1-z).^alpha.*h(0im+z) ) ) + (dh(0im+z)*(-1/2)./(sqrt(h(0im+z).^3.*(0im+1-z).^alpha) ) -alpha/2./(sqrt(h(0im+z).*(0im+1-z).^alpha).*(0im+1-z).*(-1) ) ) )     0;]*RL + [sqrt(0im+n*pi/sthosth).*(0im-2.0)^(-n+0.0)./(sqrt(h(0im+z) .*(0im+1-z).^alpha ) )     0;]*dRL )*[Dinf*(sin(psi(0im+z) -beta*pi/2 + (pi-theta)/2)*sbesowei +splsindbes ); -1im/Dinf*(sin(psi(0im+z) -beta*pi/2 -(pi-theta)/2)*sbesowei + smisindbes)] + [sqrt(0im+n*pi/sthosth).*(0im-2.0)^(-n+0.0)./(sqrt(h(0im+z).*(0im+1-z).^alpha) )     0;]*RL*[Dinf*(dscosbesoweipl +dssindbesoweipl);-1im/Dinf*(dscosbesoweimi+dssindbesoweimi)])[1,1];
    dbpiL = +dbpiL/sthosth; # *(-1.0+0.0im)^(abs(angle(z+1))>pi/2);
elseif method == "oQ"
    dbpiL = (([sqrt(0im+n).*(0.0im-1)^(-n+0.0)*(-1/2*(0.0im+sthosth)^(-3/2+0.0)*dsthosth./(sqrt(0im+(0im+1-z).^alpha.*h(0im+z) ) ) + (dh(0im+z)*(-1/2)./(sqrt(h(0im+z).^3.*(0im+1-z).^alpha) ) -alpha/2./(sqrt(h(0im+z).*(0im+1-z).^alpha).*(0im+1-z).*(-1) ) ) )     0;]*RL + [sqrt(0im+n/sthosth).*(-1)^(-n+0.0)./(sqrt(h(0im+z) .*(0im+1-z).^alpha ) )     0;]*dRL )*[(sin(psi(0im+z) -beta*pi/2 + (pi-theta)/2)*sbesowei +splsindbes ); -1im/Dinf^2*(sin(psi(0im+z) -beta*pi/2 -(pi-theta)/2)*sbesowei + smisindbes)] + [sqrt(0im+n/sthosth).*(-1)^(-n+0.0)./(sqrt(h(0im+z).*(0im+1-z).^alpha) )     0;]*RL*[(dscosbesoweipl +dssindbesoweipl);-1im/Dinf^2*(dscosbesoweimi+dssindbesoweimi)])[1,1];
    dbpiL = +dbpiL/sthosth*nmz; # *(-1.0+0.0im)^(abs(angle(z+1))>pi/2);
end

return (bpiL, dbpiL)
end

# Evaluate the asymptotic expansion for the gamma_n of the generalised Jacobi polynomial.
# Input
#   n            - Degree at which to evaluate
#   nrT          - Number of terms. 1 gives leading order term
#   Dinf         - Limit of the Szego function
#   Uright,Uleft - Right and left U-matrices
# Output
#   gamma        - Asymptotic expansion of the gamma_n
function gamman(n,nrT,Dinf,Uright,Uleft)
    if (nrT == 1) || (length(Uright) == 0)
	return 2.0.^n/Dinf/sqrt(pi);
    else
	return 2.0^(n)/Dinf/sqrt(pi)*sqrt(real(1+2*1im*Dinf^2*sum(reshape(Uright[2,1,1:(nrT-1),1] + Uleft[2,1,1:(nrT-1),1], (nrT-1), 1)./(n+1).^(1:(nrT-1) ) ) ));
    end
end

# Evaluate the asymptotic expansion for the recursion coefficient alpha_n of the generalised Jacobi polynomial.
# Input
#   n            - Degree at which to evaluate
#   nrT          - Number of terms. 1 gives leading order term
#   Uright,Uleft - Right and left U-matrices
# Output
#   alphan       - Asymptotic expansion of the alpha_n
function alphan(n,nrT,Uright,Uleft)
    if (nrT == 1) || (length(Uright) == 0)
	return 0;
    else
	return -sum(reshape(Uright[1,1,1:(nrT-1),1] + Uleft[1,1,1:(nrT-1),1], (nrT-1), 1)./(n+1).^(1:(nrT-1))) - sum(reshape(Uright[2,2,1:(nrT-1),1] + Uleft[2,2,1:(nrT-1),1],(nrT-1), 1)./n.^(1:(nrT-1)));
    end
end

# Evaluate the asymptotic expansion for the recursion coefficient beta_n of the generalised Jacobi polynomial.
# Input
#   n            - Degree at which to evaluate
#   nrT          - Number of terms. 1 gives leading order term
#   Dinf         - Limit of the Szego function
#   Uright,Uleft - Right and left U-matrices
# Output
#   betan       - Asymptotic expansion of the beta_n
function betan(n,nrT,Dinf,Uright,Uleft)
    if (nrT == 1) || (length(Uright) == 0)
	return 1/4;
    else
	return (1/(2im*Dinf^2) + sum(reshape(Uright[2,1,1:(nrT-1),1] + Uleft[2,1,1:(nrT-1),1],(nrT-1), 1)./n.^(1:(nrT-1))) )*(-Dinf^2/(2im) + sum(reshape(Uright[1,2,1:(nrT-1),1] + Uleft[1,2,1:(nrT-1),1],(nrT-1), 1)./n.^(1:(nrT-1))) );
    end
end


text:
"""
    EventSet

A set of PAPI native or preset events
"""
mutable struct EventSet
    val::Cint
    events::Vector{Event}

    function EventSet()
        evtset = new(PAPI_NULL, Event[])
        @papichk ccall((:PAPI_create_eventset, libpapi), Cint, (Ptr{Cint},), evtset)
        REFCOUNT[] += 1

        finalizer(evtset) do x
            @papichk ccall((:PAPI_cleanup_eventset, libpapi), Cint, (Cint,), evtset)
            @papichk ccall((:PAPI_destroy_eventset, libpapi), Cint, (Ptr{Cint},), evtset)
            deref_shutdown()
        end
        evtset
    end
end

# allows us to pass EventSet objects directly into ccall signatures
function Base.cconvert(::Type{Cint}, evtset::EventSet)
    evtset.val
end

# allows us to pass EventSet objects directly into Ptr{} ccall signatures
function Base.unsafe_convert(::Type{Ptr{Cint}}, evtset::EventSet)
    convert(Ptr{Cint}, pointer_from_objref(evtset))
end

"""
    try_add_event(evtset::EventSet, evt::Event)

Tries to add the event `evt` to the set. Returns success or not.
"""
function try_add_event(evtset::EventSet, evt::Event)
    success = ccall((:PAPI_add_event, libpapi), Cint, (Cint, Cuint), evtset, evt) == PAPI_OK
    if success
        push!(evtset.events, evt)
    end
    success
end

Base.eltype(::Type{EventSet}) = Event
Base.length(evtset::EventSet) = length(evtset.events)
Base.iterate(evtset::EventSet, i=1) = iterate(evtset.events, i)

"""
    push!(evtset::EventSet, evt::Event)

Adds one PAPI event to the set. The event can be either a native or preset event.
"""
function Base.push!(evtset::EventSet, evt::Event)
    @papichk ccall((:PAPI_add_event, libpapi), Cint, (Cint, Cuint), evtset, evt)
    push!(evtset.events, evt)
    evtset
end

"""
    append!(evtset::EventSet, evt::Event)

Adds multiple PAPI events to the set. The events can be either a native or preset events.
"""
function Base.append!(evtset::EventSet, evts::Vector{Event})
    errcode = ccall((:PAPI_add_events, libpapi), Cint, (Cint, Ptr{Cuint}, Cint), evtset, evts, length(evts))
    if errcode < 0
        throw(PAPIError(errcode))
    end

    if errcode != PAPI_OK
        @papichk ccall((:PAPI_remove_events, libpapi), Cint, (Cint, Ptr{Cuint}, Cint), evtset, evts, errcode)
        throw(PAPIError(PAPI_ECNFLCT))
    end

    append!(evtset.events, evts)
    evtset
end
Base.append!(evtset::EventSet, evts::AbstractVector{Event}) = append!(evtset, collect(evts))

function remove!(collection::AbstractVector, item)
    index = findfirst(isequal(item), collection)
    if index === nothing
        error("$item is not in $collection")
    else
        deleteat!(collection, index)
    end
end

"""
    delete!(evtset::EventSet, evt::Event)

Remove a PAPI events from the set. The event should have previously been added using `push!` or `append!`.
"""
function Base.delete!(evtset::EventSet, evt::Event)
    @papichk ccall((:PAPI_remove_event, libpapi), Cint, (Cint, Cuint), evtset, evt)
    remove!(evtset.events, evt)
    evtset
end

"""
    empty!(evtset::EventSet)

Remove all PAPI events from the set.
"""
function Base.empty!(evtset::EventSet)
    @papichk ccall((:PAPI_remove_events, libpapi), Cint, (Cint, Ptr{Cuint}, Cint), evtset, evtset.events, length(evtset))
    empty!(evtset.events)
    evtset
end

"""
    EventSet(evts::AbstractVector{Events})

Create an eventset containing the specified events. The events can be either native or preset PAPI events.
"""
function EventSet(evts::AbstractVector{Event})
    evtset = EventSet()
    append!(evtset, evts)
    evtset
end

"""
    EventSet(evts::Events...)

Create an eventset containing the specified events. The events can be either native or preset PAPI events.
"""
EventSet(evts::Event...) = EventSet(collect(Event, evts))

function find_add_group!(evt::Event, groups::Vector{EventSet})
    for g in groups
        if try_add_event(g, evt)
            return g
        end
    end

    push!(groups, EventSet(evt))
    last(groups)
end

"""
    measurement_groups(events::Vector{Event})

Split the list of events into subsets that can be counted at the same time.
A greedy algorithm is used, so the subsets might not be optimal. The subsets computed depends on the order of the events.
"""
function measurement_groups(events::Vector{Event})
    groups = EventSet[]

    for evt in events
        find_add_group!(evt, groups)
    end

    groups
end

"""
    collect_events(evtsets::Vector{EventSet})

Collect the event in the groups in order
"""
function collect_events(evtsets::Vector{EventSet})
    if length(evtsets) == 1
        evtsets[1].events
    else
        num_events = sum(length, evtsets)
        events = Vector{Event}(undef, num_events)

        last = 0
        for evtset in evtsets
            first = last + 1
            last = last + length(evtset)
            events[first:last] = evtset.events
        end

        events
    end
end


text:
function solvesplits!(iters::Vector{Int64},Q::Vector{Float64},A::Vector{Float64},
    children::Vector{Int64},Qold::Vector{Float64},Aold::Vector{Float64},W::Vector{Float64},
    beta::Vector{Float64},A0::Vector{Float64},c0::Vector{Float64},rho::Float64,
    f::Function,J::Function,maxiter::Int16,maxval::Float64,epsJ::Float64,epsN::Float64,hemoflag="no",
    injured = false,h = 0.,Ph = 0.,Vloss = [0.])
    if hemoflag == "no"
        CVModule.newton!(iters,Q,A,children,Qold,Aold,W,beta,A0,c0,rho,f,J,
            maxiter,maxval,epsJ,epsN);
    elseif hemoflag == "yes"
        if injured == false
            CVModule.newton!(iters,Q,A,children,Qold,Aold,W,beta,A0,c0,rho,f,J,
                maxiter,maxval,epsJ,epsN);
        elseif injured == true
            yout = zeros(3*(length(children)+1)+1);
            CVModule.model1dhemo!(yout,children,Ph,beta,rho,c0,A0,h,W[1],W[2:end])
            A[1] = yout[1];
            W[1] = yout[2];
            Q[1] = yout[3];
            Vloss[1] = yout[4];
            for i = 1:length(children)
                A[i+1] = yout[5+3*(i-1)];
                W[i+1] = yout[6+3*(i-1)];
                Q[i+1] = yout[7+3*(i-1)];
            end
        end
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
module Julz

  using Compat
  using EnglishText
  using Mustache
  using PkgDev
  using Revise

  using Reexport

  export @reexport

  @reexport using DataStructures

  cd("$(dirname(@__FILE__))/..") do
    include("../config/bootload.jl")
  end

  function main()
    load_input("input.jl", true)

    open("output.jl", "w") do cur_file
      write(cur_file, "output = []")
    end

    println("done.")
  end

end


text:
# Autogenerated wrapper script for AzStorage_jll for powerpc64le-linux-gnu
export libAzStorage

using CompilerSupportLibraries_jll
using LibCURL_jll
JLLWrappers.@generate_wrapper_header("AzStorage")
JLLWrappers.@declare_library_product(libAzStorage, "libAzStorage.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, LibCURL_jll)
    JLLWrappers.@init_library_product(
        libAzStorage,
        "lib/libAzStorage.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using HiQGA.transD_GP, PyPlot
## model fixed parts, i.e., air, but only to be compatible with AEM ...
zfixed   = [-1e5]
ρfixed   = [1e12]
nmax = 200
# z grid spec starts, the first z and first ρ will be unused in MT
zstart = 0.0
extendfrac, dz = 1.22, 15
zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=50, showplot=true, atol=1e-3)
z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)
## COPROD from Occam paper
T = [28.5 38.5 52.0 70.5 95.5 129.0 174.6 236.2 319.6 432.5 585.1 791.7 1071.1 1449.2 1960.7][:]
d_log10_ρ = [2.315  2.254  2.229  2.188  2.18  2.162  2.151  2.208  2.194  2.299  2.338  2.42  2.405  2.308  2.397][:]
d_phase_deg = -[57.19  58.19  61.39  59.09  59.89  51.19  46.89  42.79  36.89  32.0  44.0  32.0  37.59  45.29  50.09][:]
σ_log10_ρ = [0.0721  0.0425  0.0244  0.021  0.0164  0.0173  0.0287  0.0328  0.0193  0.027  0.0591  0.0506  0.0825  0.1233  0.0927][:]
σ_phase_deg = [22.95  22.95  4.96  4.46  5.96  22.95  22.95  2.46  1.65  22.95  6.37  2.46  22.95  4.15  22.95][:]
F = transD_GP.MT1DInversion.MT1DZ_nodepthprior(d_log10_ρ = d_log10_ρ, d_phase_deg = d_phase_deg,
                                    σ_log10_ρ = σ_log10_ρ, σ_phase_deg = σ_phase_deg,
                                    freqs = 1 ./T, zboundaries = zboundaries, irxlayer = 1)
transD_GP.MT1DInversion.plotdata(F)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA