text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Behavior.Gherkin

"""
The location of a step definition (the Julia code of a test) is a filename and a line number.
"""
struct StepDefinitionLocation
    filename::String
    lineno::Int
end

"Thrown if there is no step definition matching a scenario step."
struct NoMatchingStepDefinition <: Exception end

"Thrown when more than one step definition matches a scenario step."
struct NonUniqueStepDefinition <: Exception
    locations::Vector{StepDefinitionLocation}
end

"A step definition matcher takes a scenario step and finds a matching step definition to execute."
abstract type StepDefinitionMatcher end

# A StepDefinitionMatcher should define a method
# findstepdefinition(::StepDefinitionMatcher, ::Gherkin.ScenarioStep)

#
# Matching utilities
#

"""
    symbolwithdefault(s::AbstractString) :: Symbol

If the string `s` is empty, then it returns the `:String` symbol. Otherwise
it returns `Symbol(s)`.

This function simply encapsulates the fact that an empty parameter `{}` in a
step defaults to the `String` type.
"""
function symbolwithdefault(s::AbstractString) :: Symbol
    if s == ""
        :String
    else
        Symbol(s)
    end
end

function makedescriptionregex(s::String) :: Tuple{Regex, AbstractVector{Symbol}}
    variablematches = eachmatch(r"{([^{}]*)}", s)
    # The variables here are the variable types if the step has parameters.
    # So if the step is
    #   Given some value {Int} and some text {String}
    # the `variables` will hold [:Int, :String]
    variables = [symbolwithdefault(m[1]) for m in variablematches]

    # Escaping the string here so that characters that have meaning in a regular expressios
    # is treated as their actual character, instead of as regular expression characters.
    # The escaped characters are the PCRE metacharacters.
    escaped_s = escape_string(s, "\$^()|.[]?*+")

    # Replace variables like {foo} with a regex expression for that
    # variable like (?<foo>.*)
    regex_s = replace(escaped_s, r"{([^{}]*)}" => s"(.*)")

    Regex("^$(regex_s)\$"), variables
end

"A step definition has a description, which is used to find it, a function to execute, and a location."
struct StepDefinition
    description::String
    descriptionregex::Regex
    definition::Function
    location::StepDefinitionLocation
    variabletypes::Vector{Symbol}

    function StepDefinition(description::String, definition::Function, location::StepDefinitionLocation)
        descriptionregex, variables = makedescriptionregex(description)
        new(description, descriptionregex, definition, location, variables)
    end
end

struct StepDefinitionMatch
    stepdefinition::StepDefinition
    variables::Vector{Any}

    StepDefinitionMatch(s::StepDefinition) = new(s, Any[])
    function StepDefinitionMatch(s::StepDefinition, variables::AbstractArray{<:Any})
        new(s, variables)
    end
end

"""
The context in which a step definition executes. This context is used to share data between
different step definitions. It is created newly for each scenario. Thus, two scenarios cannot share
data.
"""
mutable struct StepDefinitionContext
    storage::Dict{Symbol, Any}
    datatable::Gherkin.DataTable

    StepDefinitionContext() = new(Dict{Symbol, Any}(), Gherkin.DataTable())
end

"Find a variable value given a symbol name."
Base.getindex(context::StepDefinitionContext, sym::Symbol) = context.storage[sym]

"Set a variable value given a symbol name and a value."
Base.setindex!(context::StepDefinitionContext, value::Any, sym::Symbol) = context.storage[sym] = value

"Check for a mapping for a given key"
Base.haskey(context::StepDefinitionContext, sym::Symbol) = haskey(context.storage, sym)

#
# Global state
#

currentdefinitions = Vector{StepDefinition}()
currentfilename = ""

#
# Step definition macros
#

function step_definition_(definition::Expr, description::String)
    # The step definition function takes a context and executes a bit of code supplied by the
    # test writer. The bit of code is in $definition.
    definitionfunction = :((context, vars) -> $definition(context, vars...))
    descriptionregex = makedescriptionregex(description)
    quote
        # Push a given step definition to the global state so it can be found by the
        # `StepDefinitionMatcher`.
        push!(currentdefinitions,
            StepDefinition(
                $description,
                (context, vars) -> begin
                try
                    # Escape the step definition code so it gets the proper scope.
                    $(esc(definitionfunction))(context, vars)
                    # Any step definition that does not throw an exception is successful.
                    SuccessfulStepExecution()
                catch ex
                    # StepAssertFailures are turned into a failed result here, but all other exceptions
                    # are propagated.
                    if ex isa StepAssertFailure
                        StepFailed(ex.assertion, ex.evaluated)
                    else
                        rethrow()
                    end
                end
            end,
            StepDefinitionLocation(currentfilename, 0)))
    end
end

"Provide a more user friendly @given macro for a step definition."
macro given(description, definition)
    step_definition_(description, definition)
end

"Provide a more user friendly @when macro for a step definition."
macro when(description, definition)
    step_definition_(description, definition)
end

"Provide a more user friendly @then macro for a step definition."
macro then(description, definition)
    step_definition_(description, definition)
end

#
# Step Definition Matcher implementation
#

"""
Finds step definitions defined in a Julia file with the @given, @when, @then macros defined above.
Takes a source text as input and reads the code defined in it.
"""
struct FromMacroStepDefinitionMatcher <: StepDefinitionMatcher
    stepdefinitions::Vector{StepDefinition}
    filename::String

    function FromMacroStepDefinitionMatcher(source::AbstractString; filename::String="<no filename>")
        global currentdefinitions
        global currentfilename
        currentfilename = filename
        # Read all step definitions as Julia code.
        include_string(Main, source)
        # Save the step definitions found in the global variable `currentdefinitions` into a local
        # variable, so that we can clear the global state and read another file.
        mydefinitions = currentdefinitions
        this = new(mydefinitions, filename)
        currentdefinitions = Vector{StepDefinition}()
        this
    end
end

"""
    converttypes(typesymbol::Symbol, value) :: Any

Convert `value` to the type named by `typesymbol`.
This is necessary because `Number` types are converted from strings to its primitive
type using the `parse` method, while other types are converted from strings using the
`convert` method.

# Example
```julia-repl
julia> converttypes(:Int, "123")
123
```
"""
function converttypes(typesymbol::Symbol, value) :: Any
    t = eval(typesymbol)
    if t <: Number
        parse(t, value)
    else
        convert(t, value)
    end
end

function matchdefinition(stepdefinition::StepDefinition, description::String) :: Union{StepDefinitionMatch,Nothing}
    m = match(stepdefinition.descriptionregex, description)
    if m !== nothing
        variablestrings = String[String(x) for x in m.captures]
        # Combine the parameter values captured in m with their types, that
        # we have from the stepdefinition.
        varswithtypes = zip(stepdefinition.variabletypes, variablestrings)
        # Convert each parameter value to its expected type.
        # Example: if `varswithtypes = [(:Int, "17"), (:Bool, "true")]`
        # then `variables` will be `[17, true]`.
        variables = [converttypes(typesymbol, value)
                     for (typesymbol, value) in varswithtypes]
        StepDefinitionMatch(stepdefinition, variables)
    else
        nothing
    end
end

"""
    findstepdefinition(matcher::FromMacroStepDefinitionMatcher, step::Gherkin.ScenarioStep)

Find a step definition that has a description that matches the provided scenario step.
If no such step definition is found, throw a `NoMatchingStepDefinition`.
If more than one such step definition is found, throw a `NonUniqueStepDefinition`.
"""
function findstepdefinition(matcher::FromMacroStepDefinitionMatcher, step::Gherkin.ScenarioStep) :: StepDefinitionMatch
    allsteps = map(x -> matchdefinition(x, step.text), matcher.stepdefinitions)
    matches = filter(x -> x !== nothing, allsteps)
    if isempty(matches)
        throw(NoMatchingStepDefinition())
    end
    if length(matches) > 1
        locations = map(m -> StepDefinitionLocation(matcher.filename, 0),
                        matches)
        throw(NonUniqueStepDefinition(locations))
    end
    matches[1]
end

#
# Composite matcher
#
"""
Find step definitions from multiple other step definition matchers.
"""
mutable struct CompositeStepDefinitionMatcher <: StepDefinitionMatcher
    matchers::Vector{StepDefinitionMatcher}

    CompositeStepDefinitionMatcher(matchers...) = new([matchers...])
end

function findstepdefinition(composite::CompositeStepDefinitionMatcher, step::Gherkin.ScenarioStep)
    matching = StepDefinitionMatch[]
    nonuniquesfound = StepDefinitionLocation[]
    # Recursively call `findstepdefinition(...)` on all sub-matchers.
    # When they throw a `NonUniqueStepDefinition`, record the location so it can be shown to the
    # user where the step definitions are.
    # Ignore `NonUniqueStepDefinition` exceptions, as normally all but one of the matchers will
    # throw it.
    for matcher in composite.matchers
        try
            stepdefinitionmatch = findstepdefinition(matcher, step)
            push!(matching, stepdefinitionmatch)
        catch ex
            if ex isa NonUniqueStepDefinition
                append!(nonuniquesfound, ex.locations)
            end
        end
    end
    if length(matching) > 1 || !isempty(nonuniquesfound)
        locations = vcat(nonuniquesfound, [d.stepdefinition.location for d in matching])
        throw(NonUniqueStepDefinition(locations))
    end
    if isempty(matching)
        throw(NoMatchingStepDefinition())
    end
    matching[1]
end

function addmatcher!(composite::CompositeStepDefinitionMatcher, matcher::StepDefinitionMatcher)
    push!(composite.matchers, matcher)
end


text:
function winrate_table(rate, reps)
    d = Bernoulli(rate)
    ans = zeros(Int, reps)
    for i = 1:reps
        losses = 0
        wins = 0
        while losses < 3 && wins < 7
            if rand(d) == 0
                losses += 1
            else
                wins += 1
            end
        end
        ans[i] = wins
    end
    ans
end

ext: jl

max_stars_repo_name: asieev/ArenaSim

max_issues_repo_name: asieev/ArenaSim

max_forks_repo_name: asieev/ArenaSim


text:
include("Shor.jl")
using Test

"""Euler theorem states that the order is a devisor of Eulerφ (or the size of Z* group of `N`)"""
function check_Euler_theorem(N::Int)
    Z = NumberTheory.Z_star(N)
    Nz = length(Z)   # Eulerφ
    for x in Z
        @test powermod(x,Nz,N) == 1  # the order is a devisor of Eulerφ
    end
end

@testset "Euler" begin
    check_Euler_theorem(150)
end

@testset "shor_classical" begin
    Random.seed!(129)
    L = 35
    f = shor(L, Val(:classical))
    @test f == 5 || f == 7

    L = 25
    f = shor(L, Val(:classical))
    @test f == 5

    L = 7*11
    f = shor(L, Val(:classical))
    @test f == 7 || f == 11

    L = 14
    f = shor(L, Val(:classical))
    @test f == 2 || f == 7

    @test NumberTheory.factor_a_power_b(25) == (5, 2)
    @test NumberTheory.factor_a_power_b(15) == nothing
end

@testset "shor quantum" begin
    Random.seed!(129)
    L = 15
    f = shor(L, Val(:quantum))
    @test f == 5 || f == 3
end


text:
# so GRAPE methods use n_slices and so they can store that with themselves
# but dCRAB can be either continuous or pw defined and we can choose what to do
# we will start with something piecewise since we have the integrators written for that
using Parameters
abstract type Basis end

# think about what defines a basis and the fourier basis in particular
# its always going to have the form A1 * sin(omega1 * t + phi1) + A2 * cos(omega1 * t + phi2)
Base.@kwdef mutable struct Fourier2{F,C} <: Basis
    frequencies::F = [0.0]
    coefficients::C = [(0.0, 0.0)]
end

Base.@kwdef mutable struct Sigmoid{F} <: Basis
    random_stuff::F
end

# this is only true for the fourier basis, but the point is, we could define this for every basis individually
# that gives us a method to have a unified interface
function evaluate_basis!(basis::Fourier2, out, time_axis)
    # might also want to return the lambda sometimes?

    # lam = get_lambda(basis)
    lam = function (t)
        res = 0.0
        @unpack frequencies, coefficients = basis
        # lets imagine we store tuples
        n_freq = length(frequencies)
        for i = 1:n_freq
            coeffs_i = coefficients[i]
            res =
                res +
                coeffs_i[1] * sin(frequencies[i] * t) +
                coeffs_i[2] * cos(frequencies[i] * t)
        end
        return res
    end

    out .= lam.(time_axis)
    return
end

# now lets set up an example using the bases
my_fourier_components = Fourier2([1.0, 2.0], [(1.0, 2.0), (3.0, 4.0)])

my_fourier_components.frequencies[2] = rand()

append!(my_fourier_components.coefficients, [(4.0, 5.0)])
my_fourier_components.coefficients

# using Plots
# using BenchmarkTools

time_axis = 0:0.001:5.0
hold = similar(time_axis)
evaluate_basis!(my_fourier_components, hold, time_axis)

using Plots
plot(hold)
# @benchmark evaluate_basis($my_fourier_components, $time_axis)
# @code_warntype evaluate_basis(my_fourier_components, time_axis)

# @benchmark evaluate_basis($my_fourier_components, $hold, $time_axis)

Base.@kwdef mutable struct dCRAB{ITG,B,FD,IIP,OPTS,MI,SI}
    # characteristic of the integrator
    integrator_type::ITG # or if you want to do it continuously we have another method!
    basis::B = Fourier() # store bases
    freq_dist::FD = rand # distribution of frequencies
    isinplace::IIP = true
    optim_options::OPTS = Optim.Options()
    max_iters::MI = 100
    num_SI::SI = 5 # number of superiterations you want to do
end

# separate struct that defines the different solvers that are possible

prob = Problem(
    B = [Sx, Sy],
    A = Sz,
    Xi = ρinit,
    Xt = ρfin,
    T = 1.0,
    n_controls = 2,
    guess = rand(2, 10),
    sys_type = StateTransfer(),
)

# sol = solve(prob, GRAPE(n_slices = 10, isinplace = true))

# we need to get a functional f(x) which maps a sampled array to a figure of merit
# the array x will be the previous Fourier alongside the
function _get_functional(prob)
    @unpack B, A, Xi, Xt, T, n_controls, guess, sys_type = prob
    D = size(A, 1)
    u0 = typeof(A)(I(D))

    topt = function (x)
        U = pw_evolve(A, B, x, n_controls, T / 10, 10, u0)
        ev = U * Xi * U'
        return C1(Xt, ev)
    end
    return topt
end

func = _get_functional(prob)

algg = dCRAB(n_timeslices = 10)

@unpack n_timeslices, basis, freq_dist, isinplace, optim_options, max_iters, num_SI = algg
time_axis = collect(range(0.0, prob.T, length = n_timeslices))
# for n_si = 1:num_SI
n_si = 1
# draw a frequency
ω = freq_dist()
# now from the basis we need to know that we have two coefficients, we draw them at randomly
init_coeffs = Tuple(rand(2))
init_phases = Tuple(rand(2))

# now the point is that these are our parameters in the optimisation

# so we can update the basis function I guess, since we know that we're indexing into it
# at n_si every time
basis.frequencies[n_si] = ω
basis.coefficients[n_si] = init_coeffs
basis.phases[n_si] = init_phases

test_store = zeros(n_timeslices)

evaluate_basis!(basis, test_store, time_axis)

func(test_store)
#end

struct testnew3{T,K}
    name::T
    age::K
end

mytest = testnew3("Alastair", 11)
mytest

using Parameters

@unpack name, age = mytest

function getname(test::testnew)
    @show "hi world " * test.name
end

getname(testnew("Alastair"))


text:
function create_surface(cinst, window; glfw = true)
    if glfw
        surface_ptr = GLFW.CreateWindowSurface(cinst[], window, C_NULL)
    end

    surface_ptr
end

function has_presentation_support(physical_device, queue_index, surface)
    is_supported = Ref(VkBool32(true))
    @vkcheck vkGetPhysicalDeviceSurfaceSupportKHR(physical_device, queue_index, surface, is_supported)
    Bool(is_supported[])
end


text:
# MIT License
# Copyright (c) 2019: Antonio Saragga Seabra

# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

##################
# Settlement Cycle
##################

"""
    settle(transactiondate::Date, settlementperiod, :HolydayCalendar)

Computes the settlement day using:
- Transaction date
- Settlement cycle period measured in business days (e.g. T+3, or more simply just 3)
- Business days calendars (e.g. :USNYSE, :TARGET,...)
"""
function liqd(transdate::Dates.Date, settcycle::Int, calendar)
    return advancebdays(convert(HolidayCalendar, calendar), transdate, settcycle)
end

function liqd(transdate::Dates.Date, settcycle::String, calendar)
    return advancebdays(convert(HolidayCalendar, calendar), transdate, parse(Int,settcycle[3]))
end

###################################################################
# Additional Days Count Convention
# Note: in DayCounts.jl 'const ActualActualISMA = ActualActualICMA'
###################################################################

"""
    ActualActualISMA()

**Actual/Actual ISMA** day count convention (e.g. London Stock Exchange, NASDAQ).

The year fraction is computed as:
```math
\\frac{\\text{# days between previous coupon and settlement date}}{\\text{# days in coupon period} \\times \\text{# coupon periods in a year}}
```

# Reference
- Information for private investors, January 2015, published by the London Stock Exchange Group
  https://www.londonstockexchange.com/traders-and-brokers/security-types/retail-bonds/accrued-interest-corp-supra.pdf
- Guide to Calculation Methods for the FTSE Fixed Income Indexes, v1.8
  https://www.ftse.com/products/downloads/FTSE_Fixed_Income_Index_Guide_to_Calculation_new.pdf
"""
struct ActualActualISMA <: DayCount end
function frac(startdate::Dates.Date, enddate::Dates.Date, nextcoupondate::Dates.Date, ::ActualActualISMA)
    daysbetweencoupons = Dates.value(nextcoupondate-startdate)
    if startdate > enddate || enddate > nextcoupondate
        error("start date > end date, or enddate > nextcoupondate")
    end
    return Dates.value(enddate-startdate)/(daysbetweencoupons*round(365/daysbetweencoupons))
end

###########################
# Business Days Conventions
###########################

abstract type BusinessDayConvention end
"""
    AllBD

Cashflows can land on weekends or holidays.
"""
struct AllBD <: BusinessDayConvention end
function utild(transdate::Dates.Date, ::Type{AllBD})
    return transdate
end
function utild(transdate::Dates.Date, ::Type{AllBD}, calendar)
    return transdate
end

# (A.1) if "Following" is specified, that date will be the first following day that is a Business Day;
# With following, cashflows are adjusted to land on the next good business day.
"""
    Following

Timings of cashflows that fall on a non-business day are adjusted to be on the the first **following** day that is a business day
"""
struct Following <: BusinessDayConvention end
function utild(transdate::Dates.Date,::Type{Following}, calendar)
    return tobday(convert(HolidayCalendar, calendar), transdate; forward = true)
end

# (A.2) if "Modified Following" or "Modified" is specified, that date will be the first following day that is a Business Day unless that day falls in the next calendar month, in which case that date will be the first preceding day that is a Business Day
# This is probably the most common business day convention used.  Cash flows that fall on weekends or holidays are assumed to be distributed on the following business day. However if the following business day is in a different month, the previous business day is adopted instead.
"""
    ModifiedFollowing

Timings of cashflows that fall on a non-business day are adjusted to be on the first **following** day that is a business day unless that day falls in the next calendar month, in which case that date will be the first **preceding** day that is a business day.
"""
struct ModifiedFollowing <: BusinessDayConvention end
function utild(transdate::Dates.Date,::Type{ModifiedFollowing}, calendar)
    maybe = tobday(convert(HolidayCalendar, calendar), transdate; forward = true)
    if month(maybe) == month(transdate)
        return maybe
    else
        return tobday(convert(HolidayCalendar, calendar), transdate; forward = false)
    end
end

# (A.3) if "Preceding" is specified, that date will be the first preceding day that is a Business Day.
# Cash flows that fall on a non-business day are assumed to be distributed on the previous business day.
"""
    Preceding

Timings of cashflows that fall on a non-business day are adjusted to be on the first **preceding** day that is a business day.
"""
struct Preceding <: BusinessDayConvention end
function utild(transdate::Dates.Date,::Type{Preceding}, calendar)
    return tobday(convert(HolidayCalendar, calendar), transdate; forward = false)
end

# (5) Modified Preceding
# Cash flows that fall on a non-business day are assumed to be distributed on the previous business day. However if the previous business day is in a different month, the following business day is adopted instead.
struct ModifiedPreceding <: BusinessDayConvention end
function utild(transdate::Dates.Date,::Type{ModifiedPreceding}, calendar)
    maybe = tobday(convert(HolidayCalendar, calendar), transdate; forward = false)
    if month(maybe) == month(transdate)
        return maybe
    else
        return tobday(convert(HolidayCalendar, calendar), transdate; forward = true)
    end
end

# End of Month - unadjusted
# All cash flows are assumed to be distributed on the final day of the month, even if it is a non-business day.
"""
    EndOfMonth

All cash flows are assumed to be distributed on the **final day of the month**, even if it is a non-business day.
"""
struct EndOfMonth <: BusinessDayConvention end
function utild(transdate::Dates.Date, ::Type{EndOfMonth})
    return Dates.lastdayofmonth(transdate)
end

# End of Month - Following
# All cash flows are assumed to be distributed on the final day of the month. If the final day of the month is a non-business day, the following business day is adopted.

# End of Month - Preceding
# All cash flows are assumed to be distributed on the final day of the month. If the final day of the month is a non-business day, the previous business day is adopted.

"""
    bdarray(A, convention, calendar)
    bdarray(A, convention, calendar, subsetcolumns)
        A: nxk Array
        convention: BusinessDayConvention
        calender: HolydaysCalender
        subsetcolumns: number of the columns to be processed e.g. [1,4] - first and fouth columns processed, [1] - only first column processed. If omited, all columns will be processed.

Replace holidays and weekends in array A with business days according with a choosen business day convention and calander. The array A may have elements whose type is not Date (which will be skipped).

"""
function bdarray(A, convention, calendar)
    nlines = size(A,1)
    ncolumns = size(A,2)
    for j = 1:ncolumns
        for i = 1:nlines
            if typeof(A[i,j]) == Date
                A[i,j] = utild(A[i,j], convention, calendar)
            end
        end
    end
    return A
end

function bdarray(A, convention, calendar, subsetcolumns)
    nlines = size(A,1)
    for j in subsetcolumns
        for i = 1:nlines
            if typeof(A[i,j]) == Date
                A[i,j] = utild(A[i,j], convention, calendar)
            end
        end
    end
    return A
end


text:
using OrdinaryDiffEq, LinearSolve, Test

const N = 32
const xyd_brusselator = range(0,stop=1,length=N)
brusselator_f(x, y, t) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.
limit(a, N) = a == N+1 ? 1 : a == 0 ? N : a
function brusselator_2d_loop(du, u, p, t)
  A, B, alpha, dx = p
  alpha = alpha/dx^2
  @inbounds for I in CartesianIndices((N, N))
    i, j = Tuple(I)
    x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
    ip1, im1, jp1, jm1 = limit(i+1, N), limit(i-1, N), limit(j+1, N), limit(j-1, N)
    du[i,j,1] = alpha*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) +
                B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t)
    du[i,j,2] = alpha*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) +
                A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
    end
end
p = (3.4, 1., 10., step(xyd_brusselator))

function init_brusselator_2d(xyd)
  N = length(xyd)
  u = zeros(N, N, 2)
  for I in CartesianIndices((N, N))
    x = xyd[I[1]]
    y = xyd[I[2]]
    u[I,1] = 22*(y*(1-y))^(3/2)
    u[I,2] = 27*(x*(1-x))^(3/2)
  end
  u
end
u0 = init_brusselator_2d(xyd_brusselator)
prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
                                     u0,(0.,5.0),p)

integ1 = init(prob_ode_brusselator_2d,TRBDF2(),save_everystep=false)
integ2 = init(prob_ode_brusselator_2d,TRBDF2(linsolve=IterativeSolversJL_GMRES()),save_everystep=false)

nojac = @allocated init(prob_ode_brusselator_2d,TRBDF2(linsolve=IterativeSolversJL_GMRES()),save_everystep=false)
jac = @allocated init(prob_ode_brusselator_2d,TRBDF2(),save_everystep=false)
@test jac/nojac > 50
@test integ1.cache.nlsolver.cache.jac_config !== nothing
@test integ2.cache.nlsolver.cache.jac_config === nothing

## Test that no Jac Config is created

const k1=.35e0
const k2=.266e2
const k3=.123e5
const k4=.86e-3
const k5=.82e-3
const k6=.15e5
const k7=.13e-3
const k8=.24e5
const k9=.165e5
const k10=.9e4
const k11=.22e-1
const k12=.12e5
const k13=.188e1
const k14=.163e5
const k15=.48e7
const k16=.35e-3
const k17=.175e-1
const k18=.1e9
const k19=.444e12
const k20=.124e4
const k21=.21e1
const k22=.578e1
const k23=.474e-1
const k24=.178e4
const k25=.312e1

function pollu(dy,y,p,t)
 r1  = k1 *y[1]
 r2  = k2 *y[2]*y[4]
 r3  = k3 *y[5]*y[2]
 r4  = k4 *y[7]
 r5  = k5 *y[7]
 r6  = k6 *y[7]*y[6]
 r7  = k7 *y[9]
 r8  = k8 *y[9]*y[6]
 r9  = k9 *y[11]*y[2]
 r10 = k10*y[11]*y[1]
 r11 = k11*y[13]
 r12 = k12*y[10]*y[2]
 r13 = k13*y[14]
 r14 = k14*y[1]*y[6]
 r15 = k15*y[3]
 r16 = k16*y[4]
 r17 = k17*y[4]
 r18 = k18*y[16]
 r19 = k19*y[16]
 r20 = k20*y[17]*y[6]
 r21 = k21*y[19]
 r22 = k22*y[19]
 r23 = k23*y[1]*y[4]
 r24 = k24*y[19]*y[1]
 r25 = k25*y[20]

 dy[1]  = -r1-r10-r14-r23-r24+
          r2+r3+r9+r11+r12+r22+r25
 dy[2]  = -r2-r3-r9-r12+r1+r21
 dy[3]  = -r15+r1+r17+r19+r22
 dy[4]  = -r2-r16-r17-r23+r15
 dy[5]  = -r3+r4+r4+r6+r7+r13+r20
 dy[6]  = -r6-r8-r14-r20+r3+r18+r18
 dy[7]  = -r4-r5-r6+r13
 dy[8]  = r4+r5+r6+r7
 dy[9]  = -r7-r8
 dy[10] = -r12+r7+r9
 dy[11] = -r9-r10+r8+r11
 dy[12] = r9
 dy[13] = -r11+r10
 dy[14] = -r13+r12
 dy[15] = r14
 dy[16] = -r18-r19+r16
 dy[17] = -r20
 dy[18] = r20
 dy[19] = -r21-r22-r24+r23+r25
 dy[20] = -r25+r24
end

function fjac(J,y,p,t)
      J .= 0.0
      J[1,1]   = -k1-k10*y[11]-k14*y[6]-k23*y[4]-k24*y[19]
      J[1,11]  = -k10*y[1]+k9*y[2]
      J[1,6]   = -k14*y[1]
      J[1,4]   = -k23*y[1]+k2*y[2]
      J[1,19]  = -k24*y[1]+k22
      J[1,2]   = k2*y[4]+k9*y[11]+k3*y[5]+k12*y[10]
      J[1,13]  = k11
      J[1,20]  = k25
      J[1,5]   = k3*y[2]
      J[1,10]  = k12*y[2]

      J[2,4]   = -k2*y[2]
      J[2,5]   = -k3*y[2]
      J[2,11]  = -k9*y[2]
      J[2,10]  = -k12*y[2]
      J[2,19]  = k21
      J[2,1]   = k1
      J[2,2]   = -k2*y[4]-k3*y[5]-k9*y[11]-k12*y[10]

      J[3,1]   = k1
      J[3,4]   = k17
      J[3,16]  = k19
      J[3,19]  = k22
      J[3,3]   = -k15

      J[4,4]   = -k2*y[2]-k16-k17-k23*y[1]
      J[4,2]   = -k2*y[4]
      J[4,1]   = -k23*y[4]
      J[4,3]   = k15

      J[5,5]   = -k3*y[2]
      J[5,2]   = -k3*y[5]
      J[5,7]   = 2k4+k6*y[6]
      J[5,6]   = k6*y[7]+k20*y[17]
      J[5,9]   = k7
      J[5,14]  = k13
      J[5,17]  = k20*y[6]

      J[6,6]   = -k6*y[7]-k8*y[9]-k14*y[1]-k20*y[17]
      J[6,7]   = -k6*y[6]
      J[6,9]   = -k8*y[6]
      J[6,1]   = -k14*y[6]
      J[6,17]  = -k20*y[6]
      J[6,2]   = k3*y[5]
      J[6,5]   = k3*y[2]
      J[6,16]  = 2k18

      J[7,7]   = -k4-k5-k6*y[6]
      J[7,6]   = -k6*y[7]
      J[7,14]  = k13

      J[8,7]   = k4+k5+k6*y[6]
      J[8,6]   = k6*y[7]
      J[8,9]   = k7

      J[9,9]   = -k7-k8*y[6]
      J[9,6]   = -k8*y[9]

      J[10,10] = -k12*y[2]
      J[10,2]  = -k12*y[10]+k9*y[11]
      J[10,9]  = k7
      J[10,11] = k9*y[2]

      J[11,11] = -k9*y[2]-k10*y[1]
      J[11,2]  = -k9*y[11]
      J[11,1]  = -k10*y[11]
      J[11,9]  = k8*y[6]
      J[11,6]  = k8*y[9]
      J[11,13] = k11

      J[12,11] = k9*y[2]
      J[12,2]  = k9*y[11]

      J[13,13] = -k11
      J[13,11] = k10*y[1]
      J[13,1]  = k10*y[11]

      J[14,14] = -k13
      J[14,10] = k12*y[2]
      J[14,2]  = k12*y[10]

      J[15,1]  = k14*y[6]
      J[15,6]  = k14*y[1]

      J[16,16] = -k18-k19
      J[16,4]  = k16

      J[17,17] = -k20*y[6]
      J[17,6]  = -k20*y[17]

      J[18,17] = k20*y[6]
      J[18,6]  = k20*y[17]

      J[19,19] = -k21-k22-k24*y[1]
      J[19,1]  = -k24*y[19]+k23*y[4]
      J[19,4]  = k23*y[1]
      J[19,20] = k25

      J[20,20] = -k25
      J[20,1]  = k24*y[19]
      J[20,19] = k24*y[1]

      return
end

u0 = zeros(20)
u0[2]  = 0.2
u0[4]  = 0.04
u0[7]  = 0.1
u0[8]  = 0.3
u0[9]  = 0.01
u0[17] = 0.007
prob = ODEProblem(ODEFunction(pollu, jac=fjac),u0,(0.0,60.0))

integ = init(prob,Rosenbrock23(),abstol=1e-6,reltol=1e-6)
@test integ.cache.jac_config === nothing
integ = init(prob,Rosenbrock23(linsolve=SimpleLUFactorization()),abstol=1e-6,reltol=1e-6)
@test integ.cache.jac_config === nothing
integ = init(prob,Rosenbrock23(linsolve=GenericLUFactorization()),abstol=1e-6,reltol=1e-6)
@test integ.cache.jac_config === nothing
integ = init(prob,Rosenbrock23(linsolve=RFLUFactorization(),chunk_size = Val{3}()),abstol=1e-6,reltol=1e-6)
@test integ.cache.jac_config === nothing


text:
#!/usr/bin/env julia

using Luxor, Colors

function arrow_test(fname)
  pagewidth, pageheight = 2000, 2000
  Drawing(pagewidth, pageheight, fname)
  origin() # move 0/0 to center
  background("ivory")
  setopacity(0.8)
  setline(2)
  pagetiles = Tiler(pagewidth, pageheight, 40, 8, margin=50)
  for (pos, n) in pagetiles
    randomhue()
    gsave()
    arrow(pos, Point(pos.x + 200, pos.y), linewidth=rand(1:4), arrowheadlength=rand(1:10), arrowheadangle=rand(pi/20:pi/12:pi/3))
    # same length as a non-arrow line?
    line(Point(pos.x, pos.y+10), Point(pos.x + 200, pos.y+10), :stroke)
    grestore()
  end
  finish()
  println("arrow-line-test saved in $(fname)")
end

srand(42)
fname = "/tmp/arrow-line-test.pdf"
arrow_test(fname)


text:
#=
This file is part of the replication code for: Hasenzagl, T., Pellegrino, F., Reichlin, L., & Ricco, G. (2020). A Model of the Fed's View on Inflation.
Please cite the paper if you are using any part of the code for academic work (including, but not limited to, conference and peer-reviewed papers).
=#

#=
Name: tc_main.jl
Description: Execution manager
=#

# ----------------------------------------------------------------------------------------------------------------------
# Execution: run_type == 1
# - Single iteration: it executes the code using the most updated datapoints
# ----------------------------------------------------------------------------------------------------------------------

if run_type == 1

    # ----- Standardize the data -----

    # Mixed-frequency
    if nM > 0 && nQ > 0
         σʸ_monthly   = hcat([std(collect(skipmissing(diff(data[:, i]))), dims=1) for i=1:nM]...);
         σʸ_quarterly = hcat([std(collect(skipmissing(diff(data[3:3:end, i]))), dims=1) for i=nM+1:nM+nQ]...);
         σʸ           = [σʸ_monthly σʸ_quarterly];

    # Mono frequency
    else
         σʸ = hcat([std(collect(skipmissing(diff(data[:, i]))), dims=1) for i=1:nM+nQ]...);
    end

    data = data ./ σʸ; # Standardize the data
    data = [data; missing.*ones(h, nM+nQ)]; # add "missings" for forecast

    # ----- Run the Metropolis-Within-Gibbs -----

    distr_α, distr_fcst, chain_θ_unb, chain_θ_bound, mwg_const, acc_rate, par, par_ind, par_size, distr_par =
         tc_mwg(data, h, nDraws, burnin, mwg_const, σʸ);

    # Save res in jld format
    save(string("./res_", res_name, ".jld"), Dict("distr_α" => distr_α, "distr_fcst" => distr_fcst, "chain_θ_unb" => chain_θ_unb,
                "chain_θ_bound" => chain_θ_bound, "mwg_const" => mwg_const, "acc_rate" => acc_rate, "par" => par,
                "nDraws" => nDraws, "burnin" => burnin, "data" => data, "date" => date, "nM" => nM, "nQ" => nQ,
                "MNEMONIC" => MNEMONIC, "par_ind" => par_ind, "par_size" => par_size, "distr_par" => distr_par, "σʸ" => σʸ));

# ----------------------------------------------------------------------------------------------------------------------
# Execution: run_type == 2
# - Conditional forecasts
# ----------------------------------------------------------------------------------------------------------------------

elseif run_type == 2

    # ----- Load in-sample output -----

    # Load jld output
    res_iis = load(string("./res_", res_name_iis, ".jld"));

    # Minimum output to compute the conditional forecasts
    data      = res_iis["data"];
    date      = res_iis["date"];
    σʸ        = res_iis["σʸ"];
    distr_par = res_iis["distr_par"];
    MNEMONIC  = res_iis["MNEMONIC"];

    # Remove h closing missing values
    data  = data[1:end-(size(data)[1] - size(date)[1]), :];

    # ----- Conditional forecasts -----

    for i=1:length(cond)

        # Size
        observations, n = size(data);

        # Conditioning path
        keys_ith   = collect(keys(cond[i]));
        values_ith = collect(values(cond[i]));
        pos_ith    = vcat([findall(MNEMONIC .== j) for j in keys_ith]...);
        h_ith      = length.(values_ith);

        # Generate conditioning data
        data_ith = copy(data);
        data_ith = [data; missing .* ones(maximum([h; h_ith]), n)];             # add "missings" for conditional forecast
        for j=1:length(pos_ith)
            data_ith[observations+1:observations+h_ith[j], pos_ith[j]] = values_ith[j] ./ σʸ[pos_ith[j]];
        end

        # Pre-allocate memory
        k               = size(distr_par[1].T)[1];
        m, n            = size(data_ith);
        distr_fcst_cond = zeros(m, n, nDraws[2]-burnin[2]);
        distr_α_cond    = zeros(k, m, nDraws[2]-burnin[2]);

        # Run the conditional forecast
        for draw=1:nDraws[2]-burnin[2]
            if draw > 1 && mod(draw, 100) == 0
                print("Conditional forecast $i (out of $(size(cond)[1])) > $draw-th iteration (out of $(nDraws[2]-burnin[2])) \n");
            end

            # Draw
            par_draw                 = distr_par[draw];
            par_draw.y               = permutedims(data_ith);
            α_draw, _                = kalman_diffuse!(par_draw, 0, 1, 1);
            distr_α_cond[:, :, draw] = α_draw;

            # Forecast
            distr_fcst_cond[:, :, draw] = (par_draw.Z * α_draw)' .* σʸ;
        end

        print("\n");

        # Save res in jld format
        save(string("./res_", res_name, "_cond$(i).jld"), Dict("data_ith" => data_ith, "distr_fcst_cond" => distr_fcst_cond,
                                                                "distr_α_cond" => distr_α_cond, "conditional_path" => cond[i]));
    end

# ----------------------------------------------------------------------------------------------------------------------
# Execution: run_type == 3
# -  Out-of-sample: out-of-sample exercise, forecasting period starts after end_presample_vec
# ----------------------------------------------------------------------------------------------------------------------

elseif run_type == 3

    # ----- Initialise -----

    data_full = copy(data);

    # end_presample_vec
    date_vec      = vcat([[Dates.day(date[i]) Dates.month(date[i]) Dates.year(date[i])] for i=1:length(date)]...);
    end_presample = findall(sum(date_vec .== end_presample_vec', dims=2)[:] .== 3)[1];
    end_oos       = size(date_vec)[1];

    if end_presample == end_oos
        error("end_presample_vec is not set correctly");
    end

    oos_length = end_oos-end_presample+1;

    # ----- Run the out-of-sample -----

    for t=end_presample:end_oos

        # data and data_cond
        data = data_full[1:t, :];

        # Mixed-frequency
        if nM > 0 && nQ > 0
            σʸ_monthly   = hcat([std(collect(skipmissing(diff(data[:, i]))), dims=1) for i=1:nM]...);
            σʸ_quarterly = hcat([std(collect(skipmissing(diff(data[3:3:end, i]))), dims=1) for i=nM+1:nM+nQ]...);
            σʸ           = [σʸ_monthly σʸ_quarterly];

        # Mono frequency
        else
            σʸ = hcat([std(collect(skipmissing(diff(data[:, i]))), dims=1) for i=1:nM+nQ]...);
        end

        data = data ./ σʸ;
        data = [data; missing .* ones(h, nM+nQ)];

        # Run the Metropolis-Within-Gibbs
        global mwg_const;
        distr_α, distr_fcst, chain_θ_unb, chain_θ_bound, mwg_const, acc_rate, par, par_ind, par_size, distr_par =
             tc_mwg(data, h, nDraws, burnin, mwg_const, σʸ);

        # Re-attribute standard deviation of the delta to the forecasts
        for draw=1:nDraws[2]-burnin[2]
            distr_fcst[:, :, draw] = distr_fcst[:, :, draw] .* σʸ;
        end

        # Save res for time t in jld format
        save(string("./res_", res_name, "_chunk", t-end_presample+1, ".jld"), Dict("distr_α" => distr_α, "distr_fcst" => distr_fcst,
                    "chain_θ_unb" => chain_θ_unb, "chain_θ_bound" => chain_θ_bound, "mwg_const" => mwg_const,
                    "acc_rate" => acc_rate, "par" => par, "par_ind" => par_ind, "par_size" => par_size,
                    "distr_par" => distr_par, "data" => data, "σʸ" => σʸ));
    end

    # Save res in jld format
    save(string("./res_", res_name, "_chunk0.jld"), Dict("end_presample" => end_presample, "end_oos" => end_oos,
                "oos_length" => oos_length, "nDraws" => nDraws, "burnin" => burnin, "date" => date,
                "nM" => nM, "nQ" => nQ, "MNEMONIC" => MNEMONIC, "data_full" => data_full));
end


text:
# using TrussMorph::SectionProperties

function compute_round_section_properties(radius::Float64)
    # @assert(radius > 0)
    A = pi * radius^2
    Jx = 0.5 * pi * radius^4
    Iy = pi * radius^4 / 4
    return SectionProperties(A, Jx, Iy, Iy) # assuming Iz = Iy
end


text:
@doc raw"""
    hsbondenv(p::PEPOTensor)

Returns the Hilbert-Schmidt bond evironment of a tensor
"""
function hsbondenv(vp::VectorPEPOTensor{S}, bond::S) where {S<:IndexSpace}
    hs, wu, eu = hilbertschmidt(vp)
    guess = InfiniteMPS([codomain(hs)[2]], [bond])
    ts, envs, _ = leading_boundary(
        guess, InfiniteMPO(hs), Vumps(; maxiter=100, verbose=true)
    )
    FL = permute(leftenv(envs, 1, 1, convert(MPSMultiline, ts)), Tuple(1), (2, 3))
    FR = rightenv(envs, 1, 1, convert(MPSMultiline, ts))
    AL = ts.AL[1]
    AR = ts.AR[1]
    C = ts.CR[1]
    # original γ[w2 ; e1]
    # without flipping after split:
    # 3. e1u |-<-  -<-| w2u 1.
    # 2. e1d |->-  ->-| w2d 4.
    γ = _construct_hsbondenv(hs, wu, eu, ts, envs)
    return γ
end
function hsbondenv(vhs::TensorMap{ComplexSpace,6,6}, D::Integer)
    sp = hcat(codomain(vhs)..., domain(vhs)...)
    A = TensorMap(randn, ComplexF64, ℂ^D * sp[6], sp[12] * ℂ^D)
    λ, AL, C, AR, FL, FR = vumps(A, vhs; maxiter=400,verbose=true)
    return _construct_hsbondenv(vhs, AL, C, AR, FL, FR)
end

function _construct_hsbondenv(hs, wu, eu, ts, envs)
    wui = inv(wu)
    eui = inv(eu)
    γ = TensorMap(
        undef,
        ComplexF64,
        codomain(wui)[1] * domain(eui)[2]',
        domain(eui)[1] * codomain(wui)[2]',
    )
    FL = permute(leftenv(envs, 1, 1, convert(MPSMultiline, ts)), Tuple(1), (2, 3))
    FR = rightenv(envs, 1, 1, convert(MPSMultiline, ts))
    AL = ts.AL[1]
    AR = ts.AR[1]
    C = ts.CR[1]
    @time begin
        # use @tensoropt macro, otherwise uses GiB of memory
        @tensoropt γ[w2u e1d; e1u w2d] =
            FL[8; w1 1] *
            AL[1 n1; 2] *
            C[2; 3] *
            AR[3 n2; 4] *
            FR[4 e2; 5] *
            (AR')[5; 6 s2] *
            (C')[6; 7] *
            (AL')[7; 8 s1] *
            hs[w1 s1; n1 e1] *
            eui[e1; e1u e1d] *
            hs[w2 s2; n2 e2] *
            wui[w2u w2d; w2]
    end
    return γ
end
function _construct_hsbondenv(
    M::TensorMap{S,6,6},
    AL::TensorMap{S,2,2},
    C::TensorMap{S,1,1},
    AR::TensorMap{S,2,2},
    FL::TensorMap{S,3,3},
    FR::TensorMap{S,3,3},
) where {S<:IndexSpace}
    @info "Constructing bond environment..."
    @time begin
    @tensoropt γ[-1 -2 -3 -4; -5 -6 -7 -8] :=
        M[l1 l2 l3 -2 -1 l6; l7 l8 l9 -6 -5 l12] *
        M[-4 -3 r3 r4 r5 r6; -8 -7 r9 r10 r11 r12] *
        FL[sw l8 l7; l2 l1 nw] *
        FR[r11 r10 ne; r5 r4 se] *
        AL[nw l12; l6 au] *
        C[au; bu] *
        AR[bu r12; r6 ne] *
        (AL')[l9 ad; sw l3] *
        (C')[bd; ad] *
        (AR')[r9 se; bd r3]
    end
    return γ
end

# PERF: this uses loads of memory. Probably need to come up with a better way of reshaping
# that avoids tensor contractions.
function hilbertschmidt(vp1::TensorMap{S,3,2}, vp2::TensorMap{S,3,2}) where {S<:IndexSpace}
    # @tensor vphs[wu wd su sd; nu nd eu ed] := vp1[wu su ph; nu eu] * (vp2')[wd sd; nd ed ph] old
    @info "Constructing HS tensor..."
    @time begin
        @tensoropt vphs[wu wd su sd; nu nd eu ed] :=
            vp1[wu su ph; nu eu] * (vp2')[nd ed; wd sd ph]
        @info "Size of vphs: $(Base.summarysize(vphs)/1024/1024) MiB"
    end
    cod = codomain(vphs)
    dom = domain(vphs)
    # codomain
    wu = unitary(fuse(cod[1] ⊗ cod[2]), cod[1] ⊗ cod[2])
    su = unitary(fuse(cod[3] ⊗ cod[4]), cod[3] ⊗ cod[4])
    # domain
    nu = unitary(dom[1] ⊗ dom[2], fuse(dom[1] ⊗ dom[2]))
    eu = unitary(dom[3] ⊗ dom[4], fuse(dom[3] ⊗ dom[4]))
    @info "fusing indices..."
    @time begin
        @tensoropt rv[w s; n e] :=
            wu[w; wu wd] *
            su[s; su sd] *
            vphs[wu wd su sd; nu nd eu ed] *
            nu[nu nd; n] *
            eu[eu ed; e]
    end

    # @time @tensoropt rv1[w su sd; nu nd eu ed]:=vphs[wu wd su sd; nu nd eu ed]*wu[w; wu wd]
    # @time @tensoropt rv2[w s; nu nd eu ed]:=rv1[w su sd; nu nd eu ed]*su[s; su sd]
    # @time @tensoropt rv3[w s; n eu ed]:=rv2[w s; nu nd eu ed]*nu[nu nd ; n]
    # @time @tensoropt rv[w s; n e]:=rv3[w s; n eu ed]*eu[eu ed; e]
    return rv, wu, eu
end

function hilbertschmidt(vp1::TensorMap{S,4,3}, vp2::TensorMap{S,4,3}) where {S<:IndexSpace}
    @tensoropt vhs[1 2 3 4 5 6; 7 8 9 10 11 12] :=
        vp1[1 2 3 s; 10 11 12] * (vp2')[4 5 6; 7 8 9 s]
    return vhs
end

function hilbertschmidt(
    vp1::VectorPEPOTensor{S}, vp2::VectorPEPOTensor{S}
) where {S<:EuclideanSpace}
    return hilbertschmidt(vp1.tensor, vp2.tensor)
end

hilbertschmidt(vp::TensorMap{<:IndexSpace,3,2}) = hilbertschmidt(vp, vp)
hilbertschmidt(vp::VectorPEPOTensor) = hilbertschmidt(vp.tensor, vp.tensor)

max_stars_repo_name: jack-dunham/PEPOKit

max_issues_repo_name: jack-dunham/PEPOKit

max_forks_repo_name: jack-dunham/PEPOKit


text:
using Colors
using ImageView

function show_points(im::Image, keypoints::Array{Float64, 2})
    imgc, imsl = ImageView.view(im)
    for i = 1:size(keypoints, 2)
        x, y = keypoints[1:2, i]
        ImageView.annotate!(imgc, imsl, ImageView.AnnotationPoint(x, y, shape='+', size=20, color=Colors.RGB(1,0,0)))
    end
end

# Places images side by side and draws corresponces between them
# matches; array of [ax, ay, bx, by] points
function show_matches(a::Image, b::Image, matches)
    im, sizeA = sidebyside(a, b)
    offsetB = sizeA[1]  # X dim of A is offset for b coordinates
    imgc, imsl = ImageView.view(im)
    for i = 1:size(matches, 2)
        xa, ya, xb, yb = matches[1:4, i]
        xb += offsetB
        ImageView.annotate!(imgc, imsl, ImageView.AnnotationPoint(xa, ya, shape='+', size=5, color=Colors.RGB(1,0,0)))
        ImageView.annotate!(imgc, imsl, ImageView.AnnotationPoint(xb, yb, shape='+', size=5, color=Colors.RGB(1,0,0)))
        ImageView.annotate!(imgc, imsl, ImageView.AnnotationLine(xa, ya, xb, yb, linewidth=1, color=Colors.RGB(1,0,0)))
    end
end

# Helper function to go from images to correspondences
function show_matches(a::Image, b::Image)
    ca, da = vision.sift(a)
    cb, db = vision.sift(b)
    matches, _ = find_matches(da, db, 0.5)
    coords = indices_to_coordinates(ca[1:2, :], cb[1:2, :], matches)
    show_matches(a, b, coords)
end

# Concatenate images so we can use ImageView annotations to draw lines
function sidebyside(a::Image, b::Image)
    prop = a.properties
    maxx = max(size(a, 1), size(b, 1))
    maxy = max(size(a, 2), size(b, 2))
    a = padarray(a, [0,0], [maxx - size(a, 1), maxy - size(a, 2)], "replicate")
    b = padarray(b, [0,0], [maxx - size(b, 1), maxy - size(b, 2)], "replicate")
    println(size(a))
    println(size(b))
    c = Images.Image(vcat(Images.data(a), Images.data(b)), prop)
    return (c, size(a))
end


text:
using Test
using LinearAlgebra
using QBase

@testset "./src/math/matrices.jl" begin

@testset "partial_trace()" begin
    @testset "invalid inputs" begin
        @test_throws DomainError partial_trace([1 0;0 1], [2,2], 1)
        @test_throws DomainError partial_trace([1 0 0 0 0;0 0 0 0 0;0 0 0 0 0;0 0 0 0 0], [2,2], 1)
    end

    @testset "diagonal density operator" begin
        ρ = (1/6)*Matrix(1I,6,6)

        @test partial_trace(ρ, [2,3],2) == [0.5 0; 0 0.5]
        @test partial_trace(ρ, [2,3],1) == (1/3)*I
        @test partial_trace(ρ, [3,2],2) == (1/3)*I
    end

    @testset "maximally entangled bell state" begin
        ρ = 0.5*[1 0 0 1; 0 0 0 0; 0 0 0 0; 1 0 0 1]

        @test partial_trace(ρ,[2,2],2) == [0.5 0; 0 0.5]
        @test partial_trace(ρ,[2,2],1) == [0.5 0; 0 0.5]
    end
end

@testset "computational_basis_vectors()" begin
    @test computational_basis_vectors(2) == [[1,0],[0,1]]
    @test computational_basis_vectors(5) == [
        [1,0,0,0,0],[0,1,0,0,0],[0,0,1,0,0],[0,0,0,1,0],[0,0,0,0,1]
    ]
end

@testset "commutes()" begin
    @test_throws DomainError commutes([1 0;0 1;1 1],[0 1;1 0])

    @test commutes([1 0;0 0],[0 0;0 1])
    @test !commutes([0 1;1 0],[1 0;0 -1])

    @test commutes.([[0 1;1 0],[1 0;0 -1]],[[0 im;-im 0],[0 im;-im 0]]) == [false,false]
end

@testset "is_hermitian()" begin
    @testset "invalid cases" begin
        @test !is_hermitian([1 0;1 1])
        @test !is_hermitian([1 1;-1 1])
    end
    @testset "valid cases" begin
        @test is_hermitian([1 0;0 1])
        @test is_hermitian([1 0;0 -1])
        @test is_hermitian([1 im; -im 1])
    end
    @testset "using atol" begin
        @test !is_hermitian([1 1;1+1e-6 1])
        @test is_hermitian([1 1;1+1e-6 1], atol=1e-5)

        @test !is_hermitian([1 1e-6 + 1im;-1im 1])
        @test is_hermitian([1 1e-6 + 1im;-1im 1], atol=1e-5)
    end
end

@testset "is_positive_semidefinite()" begin
    @testset "invalid cases" begin
        @test !is_positive_semidefinite([1 0 ; 0 -1])
        @test !is_positive_semidefinite([0 1;1 0])
        @test !is_positive_semidefinite([1 1;-1 1])
        @test !is_positive_semidefinite([1 0 0;0 1 0])
    end

    @testset "valid cases" begin
        @test is_positive_semidefinite([1 0; 0 1])
    end

    @testset "using atol" begin
        @test !is_positive_semidefinite([1 0;0 -1e-6])
        @test is_positive_semidefinite([1 0;0 -1e-6], atol=1e-6)
    end
end

@testset "is_orthonormal_basis()" begin
    @testset "valid cases" begin
        @test is_orthonormal_basis([[1,0],[0,1]])
        @test is_orthonormal_basis(Ket.([[1,0],[0,1]]))
        @test is_orthonormal_basis(computational_basis_vectors(5))
        @test is_orthonormal_basis([[1,1im]/sqrt(2),[1,-1im]/sqrt(2)])
    end

    @testset "invalid cases" begin
        @test !is_orthonormal_basis([[1,1]/2,[1,-1]/2])
        @test !is_orthonormal_basis([[0.5,0.5],[0.5,0.5]])
    end

    @testset "atol" begin
        @test !is_orthonormal_basis([[1,0.1],[0,1]])
        @test is_orthonormal_basis([[1,0.1],[0,1]],atol=0.2)
    end
end

@testset "n_product_id()" begin
    m0 = [1 0;0 0]
    m1 = [0 0;1 0]
    m2 = [0 0 0;0 0 0;0 1 0;0 0 0]
    m3 = [0 0 0 0;0 0 0 0;0 0 0 1]

    col_id = n_product_id([1,1],[2,2])
    row_id = n_product_id([1,2],[2,2])
    @test col_id == 1
    @test row_id == 2
    @test kron(m0,m1)[row_id,col_id] == 1

    col_id = n_product_id([1,1,2,4],[2,2,3,4])
    row_id = n_product_id([1,2,3,3],[2,2,4,3])
    @test col_id == 8
    @test row_id == 21
    @test kron(m0,m1,m2,m3)[row_id,col_id] == 1

    @testset "errors" begin
        @test_throws DomainError n_product_id([1,1,1],[2,2])
        @test_throws DomainError n_product_id([1,3],[2,2])
    end
end

@testset "is_complete()" begin
    @testset "valid cases" begin
        @test is_complete([[1 0;0  0],[0 0;0 1]])
        @test is_complete([[0.25 -0.25im;0.25im 0.25],[0.25 0.25im;-0.25im 0.25],[0.5 0;0 0.5]])
        @test is_complete(computational_basis_states(3))
        @test is_complete([[1 0;0.5 0],[0 0;0 1]], atol=0.51)
    end

    @testset "invalid cases" begin
        @test !is_complete([[1 0;0.5 0],[0 0;0 1]])
    end
end

end


text:
export getBoundaryFuncVal

"""
gdir, gneu = function getBoundaryFuncVal(gd::Function, gn::Function, M::RegularMesh,bc=(:dir,:neu,:neu,:neu))

evaluates functions gd and gn on boundary faces. Function input arguments must be arrays that store points rowwise, e.g.,

gd = X -> X[:,1].*X[:,2].*exp(X[:,3])
"""
function getBoundaryFuncVal(gd::Function, gn::Function, M::RegularMesh,bc=(:dir,:neu,:neu,:neu))
	# get boundary conditions
	gdx  = getBoundaryCondition(gd, M)
	gnx  = getBoundaryCondition(gn, M)
	gdir = zeros(0)
	gneu = zeros(0)
	for k=1:length(bc)
		if bc[k]==:dir
			gdir = [gdir; vec(gdx[k])]
		elseif bc[k]==:neu
			inout = (k<=M.dim) ? -1 : 1
			comp = mod.((1:2*M.dim) .- 1,M.dim) .+ 1
			gneu = [gneu; inout*M.h[comp[k]]*vec(gnx[k])]
		end
	end
	return gdir, gneu
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
module TestRobustLQ

using QuantEcon
using Base.Test
using FactCheck

rough_kwargs = Dict(:atol => 1e-4, :rtol => 1e-4)

# set up
a_0     = 100
a_1     = 0.5
ρ     = 0.9
sigma_d = 0.05
β    = 0.95
c       = 2
γ   = 50.0
θ   = 0.002
ac      = (a_0 - c) / 2.0

R = [0   ac    0
     ac  -a_1  0.5
     0.  0.5   0]

R = -R
Q = γ / 2

A = [1. 0. 0.
     0. 1. 0.
     0. 0. ρ]
B = [0.0 1.0 0.0]'
C = [0.0 0.0 sigma_d]'

rblq = RBLQ(Q, R, A, B, C, β, θ)
lq = LQ(Q, R, A, B, C, β)

Fr, Kr, Pr = robust_rule(rblq)

facts("Testing robustlq") do
    # test stuff
    context("test robust vs simple") do
        Fs, Ks, Ps = robust_rule_simple(rblq, Pr; tol=1e-12)

        @fact Fr --> roughly(Fs; rough_kwargs...)
        @fact Kr --> roughly(Ks; rough_kwargs...)
        @fact Pr --> roughly(Ps; rough_kwargs...)
    end

    context("test f2k and k2f") do
        K_f2k, P_f2k = F_to_K(rblq, Fr)
        F_k2f, P_k2f = K_to_F(rblq, Kr)

        @fact K_f2k --> roughly(Kr; rough_kwargs...)
        @fact F_k2f --> roughly(Fr; rough_kwargs...)
        @fact P_f2k --> roughly(P_k2f; rough_kwargs...)
    end

    context("test evaluate F") do
        Kf, Pf, df, Of, of =  evaluate_F(rblq, Fr)

        @fact Pf --> roughly(Pr; rough_kwargs...)
        @fact Kf --> roughly(Kr; rough_kwargs...)
    end

    context("test no run-time error in robust_rule_simple") do
        # this will just print out a warning
        robust_rule_simple(rblq, ones(Pr); tol=eps(), max_iter=1)
    end
end  # facts
end  # module


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    bspderiv(d::I,c::Matrix{F},k::Vector{F}
            ) where {I<:Integer,F<:AbstractFloat}

Evaluate the control points and knot sequence of the derivative of a univariate
B-Spline.

# Arguments:
- `d`: Degree of the B-Spline.
- `c`: Control points, matrix of size `(dim,nc)`.
- `k`: Knot sequence, vector of size `nk`.

# Output:
- `dc`: Control points of the derivative.
- `dk`: Knot sequence of the derivative.

# Examples:
```julia
julia> dc,dk = bspderiv(d,c,k)
```

# Reference:
- derived from Algorithm A3.3 from 'The NURBS BOOK' pg98 ([1997Piegl](@cite))
"""
function bspderiv(d::I,c::Matrix{F},k::Vector{F}
                   ) where {I<:Integer,F<:AbstractFloat}
# TODO: the original C function passed c, k, dc, and dk as reference
# TODO: this function has optimization potential

mc,nc = size(c);
nk = length(k);

# control points of the derivative
dc = zeros(F,(mc,nc-1));
for i in 0:nc-2
  tmp = d / (k[i+d+2] - k[i+2]);
  for j in 0:mc-1
    # TODO yeah 1-1 is zero have a look if that is indented to have weights set
    #      to zero
    dc[j+1,i+1] = tmp * (c[j+1,i+2] - c[j+1,i+1]);
  end # for j
end # for i

# knots of the derivative
dk = zeros(F,nk-2);
for i in 1:nk-2
  dk[i] = k[i+1];
end # for i

return dc,dk;
end # bspderiv


text:
module GoogleMaps

    import JSON, HTTP, ImageMagick
    include("utils.jl")
    include("geocode.jl")
end


text:
# export ProductOperation

# export apply_operation
# export canonize
# export iscanonical
# export domaintype
# export isidentity

# struct ProductOperation{S<:Real, F<:Tuple}
#     dimension::Int
#     factors::F
#     function ProductOperation(factors::AbstractSpaceSymmetryOperation...)
#         F = typeof(factors)
#         dim = 0
#         S = Bool
#         for f in factors
#             S = promote_type(domaintype(f), S)
#             df = dimension(f)
#             if df != 0
#                 if dim == 0
#                     dim = df
#                 else
#                     dim != df && throw(ArgumentError("All factors must have the same dimension"))
#                 end
#             else
#                 # this should not happen since `identity * something` is always `something`
#             end
#         end
#         new{S, F}(dim, factors)
#     end

#     function ProductOperation{S}(factors::AbstractSpaceSymmetryOperation...) where {S<:Real}
#         F = typeof(factors)
#         dim = 0
#         for f in factors
#             df = dimension(f)
#             if df != 0
#                 if dim == 0
#                     dim = df
#                 else
#                     dim != df && throw(ArgumentError("All factors must have the same dimension"))
#                 end
#             else
#                 # this should not happen since `identity * something` is always `something`
#             end
#         end
#         new{S, F}(dim, factors)
#     end
# end

# ## properties
# dimension(arg::ProductOperation) = arg.dimension
# domaintype(arg::ProductOperation{S,F}) where {S,F} = S
# combinable(arg::ProductOperation) = false
# isidentity(arg::ProductOperation) = isa(canonize(arg), IdentityOperation) # <- TODO: Think about whether this is a good idea

# ## operators
# import Base.==
# """
#     ==(lhs::ProductOperation, rhs::ProductOperation)

# Test for LITERAL equality.
# """
# function ==(lhs::ProductOperation, rhs::ProductOperation)
#     dimension(lhs) != dimension(rhs) && return false
#     length(lhs.factors) != length(rhs.factors) && return false
#     return all(l == r for (l, r) in zip(lhs.factors, rhs.factors))
# end

# import Base.*
# function (*)(lhs::AbstractSpaceSymmetryOperation, rhs::AbstractSpaceSymmetryOperation)
#     return ProductOperation(lhs, rhs)
# end

# function (*)(lhs::ProductOperation, rhs::AbstractSpaceSymmetryOperation)
#     return ProductOperation(lhs.factors..., rhs)
# end

# function (*)(lhs::AbstractSpaceSymmetryOperation, rhs::ProductOperation)
#     return ProductOperation(lhs, rhs.factors...)
# end

# function (*)(lhs::ProductOperation, rhs::ProductOperation)
#     return ProductOperation(lhs.factors..., rhs.factors...)
# end

# import Base.^
# function (^)(lhs::ProductOperation{S,F}, rhs::Integer)::ProductOperation where {S,F<:Tuple}
#     if rhs == 0
#         return ProductOperation()
#     elseif rhs > 0
#         return ProductOperation(vcat([collect(lhs.factors) for i in 1:rhs]...)...)
#     else
#         lhs_inv = inv(lhs)
#         return ProductOperation(vcat([collect(lhs_inv.factors) for i in 1:(-rhs)]...)...)
#     end
# end

# import Base.inv
# function inv(arg::ProductOperation)
#     return ProductOperation(reverse(inv.(arg.factors))...)
# end

# ## apply
# function apply_operation(symop::ProductOperation{S,F}, coord::AbstractVector{S}) where {S<:Real, F}
#     return foldr(apply_operation, symop.factors; init=coord)
# end

# function (symop::ProductOperation{S,F})(coord::AbstractVector{S}) where {S, F}
#     return foldr(apply_operation, symop.factors; init=coord)
# end

# ## canonical
# iscanonical(pop::ProductOperation) = false
# function iscanonical(
#             pop::ProductOperation{<:Real,
#                                   <:Tuple{<:PointOperation, <:TranslationOperation}})
#     return iscanonical(pop.factors[1]) && iscanonical(pop.factors[2])
# end

# function canonize(arg::ProductOperation)
#     _reorder(lhs::IdentityOperation, rhs::AbstractSpaceSymmetryOperation) = (rhs, lhs)
#     _reorder(lhs::AbstractSpaceSymmetryOperation, rhs::IdentityOperation) = (lhs, rhs)
#     _reorder(lhs::TranslationOperation, rhs::TranslationOperation) = (lhs, rhs)
#     _reorder(lhs::PointOperation, rhs::PointOperation) = (lhs, rhs)
#     _reorder(lhs::PointOperation, rhs::TranslationOperation) = (lhs, rhs)
#     function _reorder(lhs::TranslationOperation, rhs::PointOperation)
#         rhs_inv = inv(rhs)
#         (rhs, TranslationOperation(rhs_inv.matrix * lhs.displacement))
#     end

#     factors = collect(arg.factors)
#     # first, bubble sort
#     # order: Point, Translation, Identity (implemented in _reorder)
#     n = length(factors)
#     for m in n:-1:2, i in 1:(m-1)
#         factors[i], factors[i+1] = _reorder(factors[i], factors[i+1])
#     end

#     new_factors = AbstractSpaceSymmetryOperation[]
#     op = IdentityOperation()
#     for f in factors
#         if combinable(op, f)
#             op *= f
#         else
#             op = canonize(op)
#             op != IdentityOperation() && push!(new_factors, op)
#             op = f
#         end
#     end
#     op = canonize(op)
#     op != IdentityOperation() && push!(new_factors, op)

#     if length(new_factors) == 0
#         return IdentityOperation()
#     elseif length(new_factors) == 1
#         return first(new_factors)
#     else
#         return ProductOperation(new_factors...)
#     end
# end


text:
"""
    AnovaFactor <: AnovaEffect

A set of values for an Anova effect which has a mean square.

`name` - the name of this factor

`ss` - sum of squares

`df` - degrees of freedom

`ms` - mean square
"""
struct AnovaFactor <: AnovaEffect
    name::String
    ss::Float64
    df::Float64
    ms::Float64
end

AnovaFactor(name, ss, df) = AnovaFactor(name, ss, df, ss / df)
AnovaFactor(name, av::AnovaEffect) = AnovaFactor(name, av.ss, av.df)
AnovaFactor(av::AnovaEffect) = AnovaFactor(av.name, av.ss, av.df)

import Base.isapprox
isapprox(x::AnovaFactor, y::AnovaFactor) =
    x.name == y.name &&
    x.ss ≈ y.ss &&
    x.df == y.df &&
    (isnan(x.ms) && isnan(y.ms) || x.ms ≈ y.ms)

import Base.zero
zero(::Type{AnovaFactor}) = AnovaFactor("", 0, 1, 0)

import Base.-, Base.+
-(x::AnovaFactor, y::AnovaEffect) = AnovaFactor(x.name, x.ss - y.ss, x.df - y.df)
+(x::AnovaFactor, y::AnovaEffect) = AnovaFactor(x.name, x.ss + y.ss, x.df + y.df)


text:
module InternalForces

using DiffEqOperators: CenteredDifference

export moment, shear, torsion, bimoment, calculateDerivativeOperators

#############################################################
# Fornberg algorithm

# This implements the Fornberg (1988) algorithm (https://doi.org/10.1090/S0025-5718-1988-0935077-0)
# to obtain Finite Difference weights over arbitrary points to arbitrary order

function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real
    #=
        order: The derivative order for which we need the coefficients
        x0   : The point in the array 'x' for which we need the coefficients
        x    : A dummy array with relative coordinates, eg. central differences
               need coordinates centred at 0 while those at boundaries need
               coordinates starting from 0 to the end point
        The approximation order of the stencil is automatically determined from
        the number of requested stencil points.
    =#
    N = length(x)
    @assert order < N "Not enough points for the requested order."
    M = order
    c1 = one(T)
    c4 = x[1] - x0
    C = zeros(T, N, M+1)
    C[1,1] = 1
    @inbounds for i in 1 : N-1
        i1 = i + 1
        mn = min(i, M)
        c2 = one(T)
        c5 = c4
        c4 = x[i1] - x0
        for j in 0 : i-1
            j1 = j + 1
            c3 = x[i1] - x[j1]
            c2 *= c3
            if j == i-1
                for s in mn : -1 : 1
                    s1 = s + 1
                    C[i1,s1] = c1*(s*C[i,s] - c5*C[i,s1]) / c2
                end
                C[i1,1] = -c1*c5*C[i,1] / c2
           end
            for s in mn : -1 : 1
                s1 = s + 1
                C[j1,s1] = (c4*C[j1,s1] - s*C[j1,s]) / c3
            end
            C[j1,1] = c4 * C[j1,1] / c3
        end
        c1 = c2
    end
    #=
        This is to fix the problem of numerical instability which occurs when the sum of the stencil_coefficients is not
        exactly 0.
        https://scicomp.stackexchange.com/questions/11249/numerical-derivative-and-finite-difference-coefficients-any-update-of-the-fornb
        Stack Overflow answer on this issue.
        http://epubs.siam.org/doi/pdf/10.1137/S0036144596322507 - Modified Fornberg Algorithm
    =#
    _C = C[:,end]
    _C[div(N,2)+1] -= sum(_C)
    return _C
end

function operatorjumps(z, dm, A, order)

    diffdm = diff(dm)
    jumps = findall(x-> x > 0.0, diffdm)

    for i in eachindex(jumps)

        #left pointing single sided stencil, left of jump
        x = z[jumps[i]-4-1:jumps[i]-1]
        x0 = z[jumps[i]-1]
        coeffs=calculate_weights(order, x0, x)
        A[jumps[i]-1, 1:end] .= 0.0
        A[jumps[i]-1, jumps[i]-4-1:jumps[i]-1] = coeffs

        #left pointing single sided stencil, at jump
        x = z[jumps[i]-4:jumps[i]]
        x0 = z[jumps[i]]
        coeffs=calculate_weights(order, x0, x)
        A[jumps[i], 1:end] .= 0.0
        A[jumps[i], jumps[i]-4:jumps[i]] = coeffs

        # right pointing single sided stencil, right of jump
        x = z[jumps[i]+1:jumps[i]+1 + 4]
        x0 = z[jumps[i]+1]
        coeffs=calculate_weights(order, x0, x)
        A[jumps[i]+1, 1:end] .= 0.0
        A[jumps[i]+1, jumps[i]+1:jumps[i]+4+1] = coeffs

        # right pointing single sided stencil, two right of jump
        x = z[jumps[i]+2:jumps[i]+2 + 4]
        x0 = z[jumps[i]+2]
        coeffs=calculate_weights(order, x0, x)
        A[jumps[i]+2, 1:end] .= 0.0
        A[jumps[i]+2, jumps[i]+2:jumps[i]+4+2] = coeffs

    end

    return A

end

function calculateDerivativeOperators(z, dm)

    NumberOfNodes=length(z)
    dz=diff(z)
    dz=[dz[1]; dz; dz[end]]
    dL=dz

    #first derivative

    NthDerivative = 1
    DerivativeOrder = 2
    Az = CenteredDifference(NthDerivative, DerivativeOrder, dL, NumberOfNodes)
    Az=Array(Az)
    Az=Az[:,2:end-1]  #remove ghost nodes
    Az[1,:]=zeros(NumberOfNodes)
    Az[2,:]=zeros(NumberOfNodes)
    Az[(end-1),:]=zeros(NumberOfNodes)
    Az[end,:]=zeros(NumberOfNodes)

    #update boundary condition with single sided stencils
    order=1

    x=z[1:5]
    x0=x[1]
    coeffs=calculate_weights(order, x0, x)
    Az[1,1:5]=coeffs

    x=z[1:5]
    x0=x[2]
    coeffs=calculate_weights(order, x0, x)
    Az[2,1:5]=coeffs

    x=z[end-4:end]
    x0=x[end-1]
    coeffs=calculate_weights(order, x0, x)
    Az[end-1,end-4:end]=coeffs

    x=z[end-4:end]
    x0=x[end]
    coeffs=calculate_weights(order, x0, x)
    Az[end,end-4:end]=coeffs

    #place singled sided stencils at discontinuities
    Az = operatorjumps(z, dm, Az, order)

    #second derivative

    NthDerivative = 2
    DerivativeOrder = 2
    Azz = CenteredDifference(NthDerivative, DerivativeOrder, dL, NumberOfNodes)
    Azz=Array(Azz)
    Azz=Azz[:,2:end-1]
    Azz[1,:]=zeros(NumberOfNodes)
    Azz[2,:]=zeros(NumberOfNodes)
    Azz[(end-1),:]=zeros(NumberOfNodes)
    Azz[end,:]=zeros(NumberOfNodes)

    order=2   #update stencil to consider boundary conditions without ghost nodes

    x=z[1:5]
    x0=x[1]
    coeffs=calculate_weights(order, x0, x)
    Azz[1,1:5]=coeffs

    x=z[1:5]
    x0=x[2]
    coeffs=calculate_weights(order, x0, x)
    Azz[2,1:5]=coeffs

    x=z[end-4:end]
    x0=x[end-1]
    coeffs=calculate_weights(order, x0, x)
    Azz[end-1,end-4:end]=coeffs

    x=z[end-4:end]
    x0=x[end]
    coeffs=calculate_weights(order, x0, x)
    Azz[end,end-4:end]=coeffs

    #place singled sided stencils at discontinuities
    Azz = operatorjumps(z, dm, Azz, order)

    #third derivative

    NthDerivative = 3
    DerivativeOrder = 2
    Azzz = CenteredDifference(NthDerivative, DerivativeOrder, dL, NumberOfNodes)
    Azzz=Array(Azzz)
    Azzz=Azzz[:,2:end-1]
    Azzz[1,:]=zeros(NumberOfNodes)
    Azzz[2,:]=zeros(NumberOfNodes)
    Azzz[(end-1),:]=zeros(NumberOfNodes)
    Azzz[end,:]=zeros(NumberOfNodes)

    order=3

    x=z[1:5]
    x0=x[1]
    coeffs=calculate_weights(order, x0, x)
    Azzz[1,1:5]=coeffs

    x=z[1:5]
    x0=x[2]
    coeffs=calculate_weights(order, x0, x)
    Azzz[2,1:5]=coeffs

    x=z[end-4:end]
    x0=x[end-1]
    coeffs=calculate_weights(order, x0, x)
    Azzz[end-1,end-4:end]=coeffs

    x=z[end-4:end]
    x0=x[end]
    coeffs=calculate_weights(order, x0, x)
    Azzz[end,end-4:end]=coeffs

    #place singled sided stencils at discontinuities
    Azzz = operatorjumps(z, dm, Azzz, order)

    return Az, Azz, Azzz

end

function moment(z, dm, Δ, E, I)

    Az, Azz, Azzz = calculateDerivativeOperators(z, dm)

    M = E .* I .* Azz * Δ

    return M

end

function shear(z, dm, Δ, E, I)

    Az, Azz, Azzz = calculateDerivativeOperators(z, dm)

    V = E .* I .* Azzz * Δ

    return V

end

function torsion(z, dm, ϕ, E, G, J, Cw)

    Az, Azz, Azzz = calculateDerivativeOperators(z, dm)

    #T=GJ*dϕ/dz - ECw*dϕ3/dz3
    T=G.*J.*Az*ϕ .- E.*Cw.*Azzz*ϕ

    return T

end

function bimoment(z, dm, ϕ, E, Cw)

    Az, Azz, Azzz = calculateDerivativeOperators(z, dm)

    #B=ECwϕ''
    B=E.*Cw.*Azz*ϕ

    return B

end

end #module


text:
# # Heat equation tutorial

# In this tutorial, we'll be solving the [heat
# equation](https://en.wikipedia.org/wiki/Heat_equation):

# ``
# \frac{∂ ρcT}{∂ t} + ∇ ⋅ (-α ∇ρcT) = 0
# ``

# where
#  - `t` is time
#  - `α` is the thermal diffusivity
#  - `T` is the temperature
#  - `ρ` is the density
#  - `c` is the heat capacity
#  - `ρcT` is the thermal energy

# To put this in the form of ClimateMachine's [`BalanceLaw`](@ref
# ClimateMachine.DGMethods.BalanceLaw), we'll re-write the equation as:

# ``
# \frac{∂ ρcT}{∂ t} + ∇ ⋅ (F(α, ρcT, t)) = 0
# ``

# where
#  - ``F(α, ρcT, t) = -α ∇ρcT`` is the second-order flux

# with boundary conditions
#  - Fixed temperature ``T_{surface}`` at ``z_{min}`` (non-zero Dirichlet)
#  - No thermal flux at ``z_{min}`` (zero Neumann)

# Solving these equations is broken down into the following steps:
# 1) Preliminary configuration
# 2) PDEs
# 3) Space discretization
# 4) Time discretization
# 5) Solver hooks / callbacks
# 6) Solve
# 7) Post-processing

# # Preliminary configuration

# ## Loading code

# First, we'll load our pre-requisites:
#  - load external packages:
using MPI
using OrderedCollections
using Plots
using StaticArrays

#  - load CLIMAParameters and set up to use it:

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

#  - load necessary ClimateMachine modules:
using ClimateMachine
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods: BalanceLaw, LocalGeometry
using ClimateMachine.MPIStateArrays
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.SingleStackUtils

#  - import necessary ClimateMachine modules: (`import`ing enables us to
#  provide implementations of these structs/methods)
import ClimateMachine.DGMethods:
    vars_state_auxiliary,
    vars_state_conservative,
    vars_state_gradient,
    vars_state_gradient_flux,
    source!,
    flux_second_order!,
    flux_first_order!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    update_auxiliary_state!,
    nodal_update_auxiliary_state!,
    init_state_auxiliary!,
    init_state_conservative!,
    boundary_state!

# ## Initialization

# Define the float type (`Float64` or `Float32`)
FT = Float64;
# Initialize ClimateMachine for CPU.
ClimateMachine.init(; disable_gpu = true);

const clima_dir = dirname(dirname(pathof(ClimateMachine)));

# Load some helper functions for plotting
include(joinpath(clima_dir, "docs", "plothelpers.jl"));

# # Define the set of Partial Differential Equations (PDEs)

# ## Define the model

# Model parameters can be stored in the particular [`BalanceLaw`](@ref
# ClimateMachine.DGMethods.BalanceLaw), in this case, a `HeatModel`:

Base.@kwdef struct HeatModel{FT} <: BalanceLaw
    "Parameters"
    param_set::AbstractParameterSet = param_set
    "Heat capacity"
    ρc::FT = 1
    "Thermal diffusivity"
    α::FT = 0.01
    "Initial conditions for temperature"
    initialT::FT = 295.15
    "Bottom boundary value for temperature (Dirichlet boundary conditions)"
    T_bottom::FT = 300.0
    "Top flux (α∇ρcT) at top boundary (Neumann boundary conditions)"
    flux_top::FT = 0.0
end

# Create an instance of the `HeatModel`:
m = HeatModel{FT}();

# This model dictates the flow control, using [Dynamic Multiple
# Dispatch](https://en.wikipedia.org/wiki/Multiple_dispatch), for which
# kernels are executed.

# ## Define the variables

# All of the methods defined in this section were `import`ed in # [Loading
# code](@ref) to let us provide implementations for our `HeatModel` as they
# will be used by the solver.

# Specify auxiliary variables for `HeatModel`
vars_state_auxiliary(::HeatModel, FT) = @vars(z::FT, T::FT);

# Specify state variables, the variables solved for in the PDEs, for
# `HeatModel`
vars_state_conservative(::HeatModel, FT) = @vars(ρcT::FT);

# Specify state variables whose gradients are needed for `HeatModel`
vars_state_gradient(::HeatModel, FT) = @vars(ρcT::FT);

# Specify gradient variables for `HeatModel`
vars_state_gradient_flux(::HeatModel, FT) = @vars(α∇ρcT::SVector{3, FT});

# ## Define the compute kernels

# Specify the initial values in `aux::Vars`, which are available in
# `init_state_conservative!`. Note that
# - this method is only called at `t=0`
# - `aux.z` and `aux.T` are available here because we've specified `z` and `T`
# in `vars_state_auxiliary`
function init_state_auxiliary!(m::HeatModel, aux::Vars, geom::LocalGeometry)
    aux.z = geom.coord[3]
    aux.T = m.initialT
end;

# Specify the initial values in `state::Vars`. Note that
# - this method is only called at `t=0`
# - `state.ρcT` is available here because we've specified `ρcT` in
# `vars_state_conservative`
function init_state_conservative!(
    m::HeatModel,
    state::Vars,
    aux::Vars,
    coords,
    t::Real,
)
    state.ρcT = m.ρc * aux.T
end;

# The remaining methods, defined in this section, are called at every
# time-step in the solver by the [`BalanceLaw`](@ref
# ClimateMachine.DGMethods.BalanceLaw) framework.

# Overload `update_auxiliary_state!` to call `heat_eq_nodal_update_aux!`, or
# any other auxiliary methods
function update_auxiliary_state!(
    dg::DGModel,
    m::HeatModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    nodal_update_auxiliary_state!(heat_eq_nodal_update_aux!, dg, m, Q, t, elems)
    return true # TODO: remove return true
end;

# Compute/update all auxiliary variables at each node. Note that
# - `aux.T` is available here because we've specified `T` in
# `vars_state_auxiliary`
function heat_eq_nodal_update_aux!(
    m::HeatModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    aux.T = state.ρcT / m.ρc
end;

# Since we have second-order fluxes, we must tell `ClimateMachine` to compute
# the gradient of `ρcT`. Here, we specify how `ρcT` is computed. Note that
#  - `transform.ρcT` is available here because we've specified `ρcT` in
#  `vars_state_gradient`
function compute_gradient_argument!(
    m::HeatModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.ρcT = state.ρcT
end;

# Specify where in `diffusive::Vars` to store the computed gradient from
# `compute_gradient_argument!`. Note that:
#  - `diffusive.α∇ρcT` is available here because we've specified `α∇ρcT` in
#  `vars_state_gradient_flux`
#  - `∇transform.ρcT` is available here because we've specified `ρcT`  in
#  `vars_state_gradient`
function compute_gradient_flux!(
    m::HeatModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.α∇ρcT = -m.α * ∇transform.ρcT
end;

# We have no sources, nor non-diffusive fluxes.
function source!(m::HeatModel, _...) end;
function flux_first_order!(
    m::HeatModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) end;

# Compute diffusive flux (``F(α, ρcT, t) = -α ∇ρcT`` in the original PDE).
# Note that:
# - `diffusive.α∇ρcT` is available here because we've specified `α∇ρcT` in
# `vars_state_gradient_flux`
function flux_second_order!(
    m::HeatModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    flux.ρcT += diffusive.α∇ρcT
end;

# ### Boundary conditions

# Second-order terms in our equations, ``∇⋅(F)`` where ``F = -α∇ρcT``, are
# internally reformulated to first-order unknowns.
# Boundary conditions must be specified for all unknowns, both first-order and
# second-order unknowns which have been reformulated.

# The boundary conditions for `ρcT` (first order unknown)
function boundary_state!(
    nf,
    m::HeatModel,
    state⁺::Vars,
    aux⁺::Vars,
    n⁻,
    state⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
    _...,
)
    if bctype == 1 # bottom
        state⁺.ρcT = m.ρc * m.T_bottom
    elseif bctype == 2 # top
        nothing
    end
end;

# The boundary conditions for `ρcT` are specified here for second-order
# unknowns
function boundary_state!(
    nf,
    m::HeatModel,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n⁻,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
    _...,
)
    if bctype == 1 # bottom
        state⁺.ρcT = m.ρc * m.T_bottom
    elseif bctype == 2 # top
        diff⁺.α∇ρcT = n⁻ * m.flux_top
    end
end;

# # Spatial discretization

# Prescribe polynomial order of basis functions in finite elements
N_poly = 5;

# Specify the number of vertical elements
nelem_vert = 10;

# Specify the domain height
zmax = FT(1);

# Establish a `ClimateMachine` single stack configuration
driver_config = ClimateMachine.SingleStackConfiguration(
    "HeatEquation",
    N_poly,
    nelem_vert,
    zmax,
    param_set,
    m,
    numerical_flux_first_order = CentralNumericalFluxFirstOrder(),
);

# # Time discretization

# Specify simulation time (SI units)
t0 = FT(0)
timeend = FT(40)

# We'll define the time-step based on the [Fourier
# number](https://en.wikipedia.org/wiki/Fourier_number)
Δ = min_node_distance(driver_config.grid)

given_Fourier = FT(0.08);
Fourier_bound = given_Fourier * Δ^2 / m.α;
dt = Fourier_bound

# # Configure a `ClimateMachine` solver.

# This initializes the state vector and allocates memory for the solution in
# space (`dg` has the model `m`, which describes the PDEs as well as the
# function used for initialization). This additionally initializes the ODE
# solver, by default an explicit Low-Storage
# [Runge-Kutta](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods)
# method.

solver_config =
    ClimateMachine.SolverConfiguration(t0, timeend, driver_config, ode_dt = dt);
grid = solver_config.dg.grid;
Q = solver_config.Q;
aux = solver_config.dg.state_auxiliary;

# ## Inspect the initial conditions

# Let's export a plot of the initial state
output_dir = @__DIR__;

mkpath(output_dir);

z_scale = 100 # convert from meters to cm
z_key = "z"
z_label = "z [cm]"
z = get_z(grid, z_scale)
state_vars = SingleStackUtils.get_vars_from_nodal_stack(
    grid,
    Q,
    vars_state_conservative(m, FT),
)
aux_vars = SingleStackUtils.get_vars_from_nodal_stack(
    grid,
    aux,
    vars_state_auxiliary(m, FT),
)
all_vars = OrderedDict(state_vars..., aux_vars...);
export_plot_snapshot(
    z,
    all_vars,
    ("ρcT",),
    joinpath(output_dir, "initial_condition.png"),
    z_label,
);
# ![](initial_condition.png)

# It matches what we have in `init_state_conservative!(m::HeatModel, ...)`, so
# let's continue.

# # Solver hooks / callbacks

# Define the number of outputs from `t0` to `timeend`
const n_outputs = 5;

# This equates to exports every ceil(Int, timeend/n_outputs) time-step:
const every_x_simulation_time = ceil(Int, timeend / n_outputs);

# Create a nested dictionary to store the solution:
all_data = Dict([k => Dict() for k in 0:n_outputs]...)
all_data[0] = all_vars # store initial condition at ``t=0``

# The `ClimateMachine`'s time-steppers provide hooks, or callbacks, which
# allow users to inject code to be executed at specified intervals. In this
# callback, the state and aux variables are collected, combined into a single
# `OrderedDict` and written to a NetCDF file (for each output step `step`).
step = [1];
callback = GenericCallbacks.EveryXSimulationTime(
    every_x_simulation_time,
    solver_config.solver,
) do (init = false)
    state_vars = SingleStackUtils.get_vars_from_nodal_stack(
        grid,
        Q,
        vars_state_conservative(m, FT),
    )
    aux_vars = SingleStackUtils.get_vars_from_nodal_stack(
        grid,
        aux,
        vars_state_auxiliary(m, FT);
        exclude = [z_key],
    )
    all_vars = OrderedDict(state_vars..., aux_vars...)
    all_data[step[1]] = all_vars

    step[1] += 1
    nothing
end;

# # Solve

# This is the main `ClimateMachine` solver invocation. While users do not have
# access to the time-stepping loop, code may be injected via `user_callbacks`,
# which is a `Tuple` of callbacks in [`GenericCallbacks`](@ref ClimateMachine.GenericCallbacks).
ClimateMachine.invoke!(solver_config; user_callbacks = (callback,));

# # Post-processing

# Our solution is stored in the nested dictionary `all_data` whose keys are
# the output interval. The next level keys are the variable names, and the
# values are the values along the grid:

# To get `T` at ``t=0``, we can use `T_at_t_0 = all_data[0]["T"][:]`
@show keys(all_data[0])

# Let's plot the solution:

export_plot(
    z,
    all_data,
    ("ρcT",),
    joinpath(output_dir, "solution_vs_time.png"),
    z_label,
);
# ![](solution_vs_time.png)

# The results look as we would expect: a fixed temperature at the bottom is
# resulting in heat flux that propagates up the domain. To run this file, and
# inspect the solution in `all_data`, include this tutorial in the Julia REPL
# with:

# ```julia
# include(joinpath("tutorials", "Land", "Heat", "heat_equation.jl"))
# ```

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
@testset "Sampling Tests" begin

    model = load_model(StandardModel, model_paths["e_coli_core.json"])

    warmup, lbs, ubs = warmup_from_variability(model, Tulip.Optimizer; workers = W)

    samples = affine_hit_and_run(
        warmup,
        lbs,
        ubs;
        sample_iters = 10 * (1:3),
        workers = W,
        chains = length(W),
    )

    @test size(samples, 1) == size(warmup, 1)
    @test size(samples, 2) == size(warmup, 2) * 3 * length(W)

    @test all(samples .>= lbs)
    @test all(samples .<= ubs)
    @test all(stoichiometry(model) * samples .< TEST_TOLERANCE)
end


text:
struct CTKAlg <: SciMLBase.AbstractSteadyStateAlgorithm end

"""
    solve(prob::SciMLBase.AbstractSteadyStateProblem,
          alg::CTKAlg;
          nrm=norm,
          τstop=1e12*365*24*60*60,
          preprint="",
          maxItNewton=50)

Solves `prob` using the modified C.T.Kelley Shamanskii algorithm.
"""
function SciMLBase.solve(prob::SciMLBase.AbstractSteadyStateProblem,
                          alg::CTKAlg;
                          nrm=norm,
                          τstop=ustrip(u"s", 1e6u"Myr"),
                          preprint="",
                          maxItNewton=50)
    # Define the functions according to SciMLBase.SteadyStateProblem type
    p = prob.p
    x0 = copy(prob.u0)
    dx = copy(x0)
    function F(x)
        if SciMLBase.isinplace(prob)
            prob.f(dx, x, p)
            dx
        else
            prob.f(x, p)
        end
    end
    ∇ₓF(x) = prob.f.jac(x, p)
    # Compute `u_steady` and `resid` as per SciMLBase using my algorithm
    x_steady = NewtonChordShamanskii(F, ∇ₓF, nrm, x0, τstop; preprint=preprint, maxItNewton=maxItNewton)
    resid = F(x_steady)
    # Return the common SciMLBase solution type
    SciMLBase.build_solution(prob, alg, x_steady, resid; retcode=:Success)
end

export solve, SteadyStateProblem, CTKAlg


text:
"""
	SimplexGraph{T<:AbstractMatrix{Bool}, S<:Union{Number,AbstractVector{<:Number}}}

Struct describing a set of simplices with vertices numbered from 1:N.
Optionally, each simplex can be assigned a weight.

Usage
--------------

	SimplexGraph(G, w=true)

Construct a SimplexGraph. `G` is a `MxN` matrix of booleans, where `M` is the number of vertices and `N` is the number of simplices.
Each column in `G` represents one simplex. True means that a vertex is part of the simplex and false that it is not.
Optionally, a vector `w` of length N can be used to specify a weight (total mass) for each simplex. The weight defaults to 1 (true).

See also `pma`, `groupsimplices`, `timeseriessimplices`, `neighborsimplices`.
"""
struct SimplexGraph{T<:AbstractMatrix{Bool}, S<:Union{Number,AbstractVector{<:Number}}}
	G::T
	w::S

	function SimplexGraph(G::T, w::S=true) where {T<:AbstractMatrix{Bool}, S<:Union{Number,AbstractVector{<:Number}}}
		@assert ndims(w)==0 || (size(w,1)==size(G,2))
		new{T,S}(G,w)
	end
end

"""
	groupsimplices(groupby::AbstractVector)

Create SimplexGraph connecting elements with identical values in the groupby vector.

# Examples
```
julia> sg = groupsimplices(["A","A","B","C","B"]);

julia> sg.G
5×3 BitArray{2}:
 1  0  0
 1  0  0
 0  1  0
 0  0  1
 0  1  0

julia> sg.w
3-element Array{Int64,1}:
 2
 2
 1
```
"""
function groupsimplices(groupby::AbstractVector)
	u = unique(groupby)
	ind = indexin(groupby,u)
	I = Vector{Int}(undef,length(groupby))
	J = Vector{Int}(undef,length(groupby))
	w = zeros(Int,length(u))
	for (i,j) in enumerate(ind)
		I[i] = i
		J[i] = j
		w[j] += 1
	end
	SimplexGraph(sparse(I,J,true),w)
end

"""
	timeseriessimplices(time::AbstractVector; groupby::AbstractVector)

Create SimplexGraph connecting elements adjacent in time.
In case of ties, **all** elements at a unique timepoint will be connected to the **all** elements at the previous, current and next timepoints.
If `groupby` is specified, the elements are first divided by group, and then connected by time.

# Examples
```
julia> sg = timeseriessimplices([0.5, 1.0, 4.0]);

julia> sg.G
3×3 BitArray{2}:
 1  1  0
 1  1  1
 0  1  1

julia> sg = timeseriessimplices([0.5, 1.0, 1.0, 4.0]);

julia> sg.G
4×3 BitArray{2}:
 1  1  0
 1  1  1
 1  1  1
 0  1  1

julia> sg.w
3-element Array{Int64,1}:
 1
 2
 1

julia> sg = timeseriessimplices([2, 4, 6, 2, 4, 8]; groupby=["A","A","A","B","B","B"]);

julia> sg.G
6×6 BitArray{2}:
 1  1  0  0  0  0
 1  1  1  0  0  0
 0  1  1  0  0  0
 0  0  0  1  1  0
 0  0  0  1  1  1
 0  0  0  0  1  1
```
"""
function timeseriessimplices(time::AbstractVector; groupby::AbstractVector=falses(length(time)))
	N = length(groupby)

	I,J = Int[],Int[]
	w = Int[]

	# TODO: Rewrite. Can be simplified and optimized.
	col = 1
	for g in unique(groupby)
		ind = findall( groupby.==g )
		time2 = time[ind]
		uniqueTimes = sort(unique(time2))

		currInd = Int[]
		nextInd = ind[time2.==uniqueTimes[1]]

		for t in Iterators.drop(uniqueTimes,1)
			prevInd = currInd
			currInd = nextInd
			nextInd = ind[time2.==t]

			append!(I, prevInd)
			append!(I, currInd)
			append!(I, nextInd)
			append!(J, repeat([col],outer=length(prevInd)+length(currInd)+length(nextInd)))
			push!(w, length(currInd))
			col += 1
		end
		if length(uniqueTimes)==2 # Two timepoints means the same simplex is used for both.
			w[end] += length(nextInd)
		else
			append!(I, currInd)
			append!(I, nextInd)
			append!(J, repeat([col],outer=length(currInd)+length(nextInd)))
			push!(w, length(nextInd))
			col += 1
		end
	end
	SimplexGraph(sparse(I,J,true),w)
end

"""
	neighborsimplices(A::AbstractMatrix; k, r, dim, symmetric, normalizedist, groupby)

Create simplex graph connecting nearest neighbor samples.

# Inputs
* `A`: Data matrix (variables × samples).
* `k`: Number of nearest neighbors to connect. Default: `0`.
* `r`: Connected all neighbors with disctance `≤r`. Default: `0.0`.
* `dim`: Reduce the dimension to `dim` before computing distances. Useful to reduce noise. Default: Disabled.
* `symmetric`: Make the simplex graph symmetric. Default: `false`.
* `normalizedist`: Normalize distances to the scale [0.0,1.0] such that the maximal distance from a point to the origin is 0.5. Affects the `r` parameter. Default: `true`.
* `groupby`: Only connected samples within the specified groups. Default: Disabled.

# Examples
```
julia> sg = neighborsimplices([0 0 2 2; 0 1 1 0]; k=1); sg.G
4×4 BitArray{2}:
 1  1  0  0
 1  1  0  0
 0  0  1  1
 0  0  1  1

julia> sg = neighborsimplices([0 0 2 2; 0 1 1 0]; r=0.45); sg.G
4×4 BitArray{2}:
 1  1  0  1
 1  1  1  0
 0  1  1  1
 1  0  1  1
```
"""
function neighborsimplices(A::AbstractMatrix; dim::Integer=typemax(Int), kwargs...)
	if dim<minimum(size(A))
		F = svdbyeigen(A; nsv=dim)
		A = Diagonal(F.S)*F.Vt
	end
	_neighborsimplices(A; kwargs...)
end

function _neighborsimplices(A::AbstractMatrix; k::Integer=0, r::Real=0.0, symmetric=false, normalizedist=true, groupby=nothing)
	k==0 && r==0.0 && return SimplexGraph(sparse(Diagonal(trues(size(A,2))))) # trivial case, avoid some computations.
	r>0.0 && normalizedist && (r *= 2*sqrt(maximum(sum(x->x^2, A, dims=1))))

	N = size(A,2)
	rowInd,colInd = Int[],Int[]

	if groupby===nothing
		_neighborsimplices!(rowInd, colInd, A, identity, k, r)
	else
		for g in unique(groupby)
			ind = findall( groupby.==g )
			_neighborsimplices!(rowInd, colInd, A[:,ind], x->ind[x], k, r)
		end
	end

	SimplexGraph(symmetric ? sparse([rowInd;colInd],[colInd;rowInd],true) : sparse(rowInd,colInd,true))
end

function _neighborsimplices!(rowInd, colInd, A, indfun, k, r)
	tree = BallTree(A)

	if k>0
		indices,_ = knn(tree, A, min(k+1,size(A,2)))
		for (j,I) in enumerate(indices)
			append!(rowInd,indfun.(I))
			append!(colInd,repeat([indfun(j)],outer=length(I)))
		end
	end

	if r>0.0
		indices = inrange(tree, A, r)
		for (j,I) in enumerate(indices)
			append!(rowInd,indfun.(I))
			append!(colInd,repeat([indfun(j)],outer=length(I)))
		end
	end
end


text:
"""
mcppv(ts::AbstractVector; shuffle::Int = 500)
    return p-value of a change point of mean in ts
    shuffle is number of randomization used to generate distribution
    uses global random number generator
"""
function mcppv(ts::AbstractVector; shuffle::Int = 500)
    length(ts) <= 1 && return 1.                    # large pvalue is no change

    dev = ts .- mean(ts)                            # centered
    cusum = Vector{Float64}(undef, length(ts))
    cumsum!(cusum, dev)                             # cumsum
    cmax_obs = maximum(abs, cusum)                  # max abs cumsum observed

    n = 0
    for _ in 1:shuffle
        shuffle!(dev)
        cumsum!(cusum, dev)
        cmax = maximum(abs, cusum)
        if cmax >= cmax_obs
            n += 1
        end
    end

    round((n + 1) / (shuffle + 1), digits = 5)      # p-value
end

"""
mcptime(ts::AbstractVector; minlen::Int = 1)
    return index of start of new segment via minimum ssq
    0 is returned if minlen is not least twice length of ts
    minlen is minimum length of segments to consider
"""
function mcptime(ts::AbstractVector; minlen::Int = 3)
    minlen < 1 && return 0
    len = length(ts)
    len < 2*minlen && return 0

    chgpoint = 0
    minssq = Inf
    @inbounds for brk = minlen+1:len+1-minlen
        @views ssqtot = ssq(ts[1:brk-1]) + ssq(ts[brk:len])
        if ssqtot < minssq
            minssq = ssqtot
            chgpoint = brk
        end
    end

    chgpoint
end

"""
mcptimechk(ts::AbstractVector; minlen::Int = 1)
    return index of start of new segment via minimum ssq
    0 is returned if minlen is not least twice length of ts
    minlen is minimum length of segments to consider
"""
function mcptimechk(ts::AbstractVector; minlen::Int = 3)
    minlen < 1 && return 0
    len = length(ts)
    len < 2*minlen && return 0

    chgpoint = 0
    minssq = Inf
    @inbounds for brk = minlen+1:len+1-minlen
        @views ssq_l = ssq(ts[1:brk-1])
        @views ssq_r = ssq(ts[brk:len])
        ssqtot = ssq_l + ssq_r
        println("Index $brk   $ssqtot $ssq_l $ssq_r")
        if ssqtot < minssq
            minssq = ssqtot
            chgpoint = brk
        end
    end

    chgpoint
end

"""
ssq(v) returns sum of squares
"""
@inline ssq(v::AbstractVector) = var(v) * (length(v) - 1)

"""
mcpoint(ts::AbstractVector; pcut=0.05, shuffle=500, minlen=3)
    return index of change point in mean if it exists, 0 otherwise
    pcut is p-value cut-off, shuffle is number of randomizations
    minlen is minimum length to consider for a segment
"""
function mcpoint(ts::AbstractVector; pcut=0.05, shuffle::Int=500, minlen::Int=3)
    minlen < 1 && return 0
    len = length(ts)
    len < 2*minlen && return 0

    chgpoint = 0
    if mcppv(ts; shuffle = shuffle) <= pcut
        brk = mcptime(ts; minlen = minlen)
        if minlen < brk <= (len - minlen + 1)
            chgpoint = brk
        end
    end

    chgpoint
end

"""
mcplast(ts::AbstractVector; pcut=0.05, shuffle=500, minlen=3)
    return index of rightmost change point if any, 0 otherwise
    pcut is p-value cut-off, shuffle is number of randomizations
    minlen is the minimum length of a segment
"""
function mcplast(ts::AbstractVector; pcut=0.05, shuffle::Int=500, minlen::Int=3)
    minlen < 1 && return 0
    len = length(ts)
    len < 2*minlen && return 0

    head = 1
    brk = mcpoint(ts[head:len]; pcut=pcut, minlen=minlen, shuffle=shuffle)
    while brk > 0
        head += brk - 1         # keep checking right hand segment
        brk = mcpoint(ts[head:len]; pcut=pcut, minlen=minlen, shuffle=shuffle)
    end

    head == 1 ? 0 : head
end

##---   all change points if they exist

"""
mcpall(ts::AbstractVector; pcut = 0.05, shuffle = 500, minlen = 3)
    return vector of index of all change points via recursive partitioning
    pcut is p-value cut-off, shuffle is number of randomizations
    minlen is the minimum length of a segment
"""
function mcpall(ts::AbstractVector; pcut=0.05, shuffle::Int=500, minlen::Int=3)
    chgpts = Vector{Int}(undef, 0)
    minlen < 1 && return chgpts
    len = length(ts)
    len < 2*minlen && return chgpts

    _mcpall!(chgpts, 1, ts; pcut=pcut, shuffle=shuffle, minlen=minlen)
    sort!(chgpts)
end

function _mcpall!(chgpts::AbstractVector, head::Int, ts::AbstractVector;
                pcut=0.05, shuffle::Int=500, minlen::Int=3)
    len = length(ts)
    chg = mcpoint(ts; pcut=pcut, shuffle=shuffle, minlen=minlen)
    chg == 0 && return                      # no change point, stop
    append!(chgpts, head + chg - 1)         # remember change point
    _mcpall!(chgpts, head, ts[1:chg-1]; pcut=pcut, shuffle=shuffle, minlen=minlen)
    _mcpall!(chgpts, head+chg-1, ts[chg:len]; pcut=pcut, shuffle=shuffle, minlen=minlen)

    nothing
end

##---   mcplot

"""
mcplot(ts::AbstractVector; chgpts=Int[], palette=:seaborn_bright)
    returns plot of ts with chgpts as change points (default is none)
    palette is the default color scheme to use

    e.g., mcplot(ts) is a plain plot of the time series
    mcplot(ts, chgpts=mcpall(ts)) is plot of ts with mcpall(ts) change points
"""
function mcplot(ts::AbstractVector; chgpts=Int[], palette=:seaborn_bright)
    len = length(ts)
    plt = plot(size=(600, 400), bar_width=0.9, xlims=(0.5, len+0.5),
        titlefont=(8), xtickfontsize=6, palette=palette)

    vline!(-0.5 .+ chgpts, lc=:blue, ls=:dash, lw=1, label=false)

    n = 0
    for (t1,t2) in zip(reverse([1; chgpts]), reverse([chgpts .- 1; length(ts)]))
        n += 1
        m = mean(ts[t1:t2])
        plot!(t1:t2, ts[t1:t2], st=:bar, lc=nothing, label=false, fc=n)
        plot!([t1-0.5, t2+0.5], [m, m], lc=:blue, ls=:dot, lw=2, label=false)
    end

    plt
end


text:
# The liblzma Interfaces
# ======================

# Return code
const LZMA_OK                = Cint(0)
const LZMA_STREAM_END        = Cint(1)
const LZMA_NO_CHECK          = Cint(2)
const LZMA_UNSUPPORTED_CHECK = Cint(3)
const LZMA_GET_CHECK         = Cint(4)
const LZMA_MEM_ERROR         = Cint(5)
const LZMA_MEMLIMIT_ERROR    = Cint(6)
const LZMA_FORMAT_ERROR      = Cint(7)
const LZMA_OPTIONS_ERROR     = Cint(8)
const LZMA_DATA_ERROR        = Cint(9)
const LZMA_BUF_ERROR         = Cint(10)
const LZMA_PROG_ERROR        = Cint(11)

# Action code
const LZMA_RUN          = Cint(0)
const LZMA_SYNC_FLUSH   = Cint(1)
const LZMA_FULL_FLUSH   = Cint(2)
const LZMA_FULL_BARRIER = Cint(4)
const LZMA_FINISH       = Cint(3)

# Flag
const LZMA_TELL_NO_CHECK          = UInt32(0x01)
const LZMA_TELL_UNSUPPORTED_CHECK = UInt32(0x02)
const LZMA_TELL_ANY_CHECK         = UInt32(0x04)
const LZMA_IGNORE_CHECK           = UInt32(0x10)
const LZMA_CONCATENATED           = UInt32(0x08)

# Check
const LZMA_CHECK_NONE     = Cint(0)
const LZMA_CHECK_CRC32    = Cint(1)
const LZMA_CHECK_CRC64    = Cint(4)
const LZMA_CHECK_SHA256   = Cint(10)

mutable struct LZMAStream
    next_in::Ptr{UInt8}
    avail_in::Csize_t
    total_in::UInt64

    next_out::Ptr{UInt8}
    avail_out::Csize_t
    total_out::UInt64

    allocator::Ptr{Cvoid}

    internal::Ptr{Cvoid}

    reserved_ptr::NTuple{4,Ptr{Cvoid}}
    reserved_uint::NTuple{2,UInt64}
    reserved_size::NTuple{2,Csize_t}
    reserved_enum::NTuple{2,Cint}
end

function LZMAStream()
    return LZMAStream(
        C_NULL, 0, 0,
        C_NULL, 0, 0,
        C_NULL,
        C_NULL,
        (C_NULL, C_NULL, C_NULL, C_NULL),
        (0, 0), (0, 0),
        (0, 0))
end

function lzmaerror(stream::LZMAStream, code::Cint)
    error(lzma_error_string(code))
end

function lzma_error_string(code::Cint)
    return "lzma error: code = $(code)"
end

function easy_encoder(stream::LZMAStream, preset::Integer, check::Integer)
    return ccall(
        (:lzma_easy_encoder, liblzma),
        Cint,
        (Ref{LZMAStream}, UInt32, Cint),
        stream, preset, check)
end

function stream_decoder(stream::LZMAStream, memlimit::Integer, flags::Integer)
    return ccall(
        (:lzma_stream_decoder, liblzma),
        Cint,
        (Ref{LZMAStream}, UInt64, UInt32),
        stream, memlimit, flags)
end

function code(stream::LZMAStream, action::Integer)
    return ccall(
        (:lzma_code, liblzma),
        Cint,
        (Ref{LZMAStream}, Cint),
        stream, action)
end

function free(stream::LZMAStream)
    ccall((:lzma_end, liblzma), Cvoid, (Ref{LZMAStream},), stream)
end


text:
function bogosort!(arr::AbstractVector)
    while !issorted(arr)
        shuffle!(arr)
    end
    return arr
end

v = rand(-10:10, 10)
println("# unordered: $v\n -> ordered: ", bogosort!(v))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode