text:
# This code was originally taken and modifed from LineSearches.jl for the sake of
# improving performance. Attribution / credit belongs to the original authors.
# This was the attached liscense agreement:
#
# The original code has been separated out from the Optim.jl package,
# licensed under the MIT License:
# > Copyright (c) 2012: John Myles White and other contributors.
#
# The LineSearches.jl package is licensed under the MIT "Expat" License:
#
# > Copyright (c) 2016: Asbjørn Nilsen Riseth and other contributors.
# >
# > 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.

@with_kw struct BackTracking{TF, TI}
    c_1::TF = 1e-4
    ρ_hi::TF = 0.5
    ρ_lo::TF = 0.1
    iterations::TI = 1_000
    order::TI = 3
    maxstep::TF = Inf
end
BackTracking{TF}(args...; kwargs...) where TF = BackTracking{TF,Int}(args...; kwargs...)

function linesearch!(ls::BackTracking, df, x::AbstractArray{T}, s::AbstractArray{T},
                            α_0::Tα, x_new::AbstractArray{T}, ϕ_0, dϕ_0) where {T, Tα}

    alphamax = convert(real(T), Inf)
    # ϕ, dϕ = make_ϕ_dϕ(df, x_new, x, s)

    # α_0 = min(α_0, min(alphamax, ls.maxstep / norm(s, Inf)))
    α_0 = min(α_0, min(alphamax, ls.maxstep / SIMDArrays.maximum_abs(s)))
    # ls(ϕ, α_0, ϕ_0, dϕ_0)

    c_1, ρ_hi, ρ_lo, iterations, order = ls.c_1, ls.ρ_hi, ls.ρ_lo, ls.iterations, ls.order
    # @unpack c_1, ρ_hi, ρ_lo, iterations, order = ls

    iterfinitemax = -log2(eps(real(Tα)))

    # @assert order in (2,3)
    # Check the input is valid, and modify otherwise
    #backtrack_condition = 1.0 - 1.0/(2*ρ) # want guaranteed backtrack factor
    #if c_1 >= backtrack_condition
    #    warn("""The Armijo constant c_1 is too large; replacing it with
    #                   $(backtrack_condition)""")
    #   c_1 = backtrack_condition
    #end

    # Count the total number of iterations
    iteration = 0

    ϕx_0, ϕx_1 = ϕ_0, ϕ_0

    α_1, α_2 = α_0, α_0

    @fastmath @inbounds @simd for i ∈ 1:full_length(x_new)
        x_new[i] = x[i] + α_1 * s[i]
    end
    # Evaluate f(x+α*s)

    ϕx_1 = value!(df, x_new)
    # @show ϕx_1, α_1

    # Hard-coded backtrack until we find a finite function value
    iterfinite = 0
    while !isfinite(ϕx_1) && iterfinite < iterfinitemax
        iterfinite += 1
        α_1 = α_2
        α_2 = 0.5α_1

        @fastmath @inbounds @simd for i ∈ 1:full_length(x_new)
            x_new[i] = x[i] + α_2 * s[i]
        end
        # Evaluate f(x+α*s)

        ϕx_1 = value!(df, x_new)
        # @show ϕx_1, α_2
    end

    # Backtrack until we satisfy sufficient decrease condition
    while ϕx_1 > ϕ_0 + c_1 * α_2 * dϕ_0
        # Increment the number of steps we've had to perform
        iteration += 1

        # Ensure termination
        # if iteration > iterations
        #     throw(LineSearchException("Linesearch failed to converge, reached maximum iterations $(iterations).",
        #                               α_2))=
        # end
        # iteration > iterations && throw("Exceeded maximum number of iterations.")
        # Shrink proposed step-size:
        if order == 2 || iteration == 1
            # backtracking via quadratic interpolation:
            # This interpolates the available data
            #    f(0), f'(0), f(α)
            # with a quadractic which is then minimised; this comes with a
            # guaranteed backtracking factor 0.5 * (1-c_1)^{-1} which is < 1
            # provided that c_1 < 1/2; the backtrack_condition at the beginning
            # of the function guarantees at least a backtracking factor ρ.
            α_tmp = - (dϕ_0 * α_2^2) / ( 2 * (ϕx_1 - ϕ_0 - dϕ_0*α_2) )
        else
            div = one(Tα) / (α_1^2 * α_2^2 * (α_2 - α_1))
            a = (α_1^2*(ϕx_1 - ϕ_0 - dϕ_0*α_2) - α_2^2*(ϕx_0 - ϕ_0 - dϕ_0*α_1))*div
            b = (-α_1^3*(ϕx_1 - ϕ_0 - dϕ_0*α_2) + α_2^3*(ϕx_0 - ϕ_0 - dϕ_0*α_1))*div

            if isapprox(a, zero(a), atol=eps(real(Tα)))
                α_tmp = dϕ_0 / (2*b)
            else
                # discriminant
                d = max(b^2 - 3*a*dϕ_0, Tα(0))
                # quadratic equation root
                @fastmath α_tmp = (-b + sqrt(d)) / (3*a)
            end
        end

        α_1 = α_2

        # α_tmp = NaNMath.min(α_tmp, α_2*ρ_hi) # avoid too small reductions
        # α_2 = NaNMath.max(α_tmp, α_2*ρ_lo) # avoid too big reductions
        α_tmp = min(α_tmp, α_2*ρ_hi) # avoid too small reductions
        α_2 = max(α_tmp, α_2*ρ_lo) # avoid too big reductions

        @fastmath @inbounds @simd for i ∈ 1:full_length(x_new)
            x_new[i] = x[i] + α_2 * s[i]
        end

        # Evaluate f(x) at proposed position
        ϕx_0, ϕx_1 = ϕx_1, value!(df, x_new)
        # @show ϕx_1, α_2
    end

    return α_2, ϕx_1
end


text:
@testset "EDF Export" begin

    n_records = 100
    edf, edf_channel_indices = make_test_data(MersenneTwister(42), 256, 512, n_records)
    uuid = uuid4()

    onda_samples, nt = edf_to_onda_samples(edf)
    annotations = edf_to_onda_annotations(edf, uuid)

    signal_names = ["eeg", "eog", "ecg", "emg", "heart_rate", "tidal_volume",
                    "respiratory_effort", "snore", "positive_airway_pressure",
                    "pap_device_leak", "pap_device_cflow", "sao2", "ptaf"]
    samples_to_export = onda_samples[indexin(signal_names, getproperty.(getproperty.(onda_samples, :info), :kind))]
    exported_edf = onda_to_edf(samples_to_export, annotations)
    @test exported_edf.header.record_count == 200
    offset = 0
    for signal_name in signal_names
        samples = only(filter(s -> s.info.kind == signal_name, onda_samples))
        channel_names = samples.info.channels
        edf_indices = (1:length(channel_names)) .+ offset
        offset += length(channel_names)
        samples_data = Onda.decode(samples).data
        edf_samples = mapreduce(transpose ∘ EDF.decode, vcat, exported_edf.signals[edf_indices])
        @test isapprox(samples_data, edf_samples; rtol=0.02)
        for (i, channel_name) in zip(edf_indices, channel_names)
            s = exported_edf.signals[i]
            @test s.header.label == OndaEDF.export_edf_label(signal_name, channel_name)
            @test s.header.physical_dimension == OndaEDF.onda_to_edf_unit(samples.info.sample_unit)
        end
    end
    @testset "Record metadata" begin
        function change_sample_rate(samples; sample_rate)
            info = SamplesInfo(Tables.rowmerge(samples.info; sample_rate=sample_rate))
            new_data = similar(samples.data, 0, Onda.index_from_time(sample_rate, Onda.duration(samples)) - 1)
            return Samples(new_data, info, samples.encoded; validate=false)
        end

        eeg_samples = only(filter(row -> row.info.kind == "eeg", onda_samples))
        ecg_samples = only(filter(row -> row.info.kind == "ecg", onda_samples))

        massive_eeg = change_sample_rate(eeg_samples, sample_rate=5000.0)
        @test OndaEDF.edf_record_metadata([massive_eeg]) == (1000000, 1 / 5000)

        chunky_eeg = change_sample_rate(eeg_samples; sample_rate=9999.0)
        chunky_ecg = change_sample_rate(ecg_samples; sample_rate=425.0)
        @test_throws OndaEDF.RecordSizeException OndaEDF.edf_record_metadata([chunky_eeg, chunky_ecg])

        e_notation_eeg = change_sample_rate(eeg_samples; sample_rate=20_000_000.0)
        @test OndaEDF.edf_record_metadata([e_notation_eeg]) == (4.0e9, 1 / 20_000_000)

        too_big_and_thorny_eeg = change_sample_rate(eeg_samples; sample_rate=20_576_999.0)
        @test_throws OndaEDF.EDFPrecisionError OndaEDF.edf_record_metadata([too_big_and_thorny_eeg])

        floaty_eeg = change_sample_rate(eeg_samples; sample_rate=256.5)
        floaty_ecg = change_sample_rate(ecg_samples; sample_rate=340)
        @test OndaEDF.edf_record_metadata([floaty_eeg, floaty_ecg]) == (100, 2)

        resizable_eeg = change_sample_rate(eeg_samples; sample_rate=25.25)
        resizable_ecg = change_sample_rate(ecg_samples; sample_rate=10.4)
        @test OndaEDF.edf_record_metadata([resizable_eeg, resizable_ecg]) == (10, 20)

        @testset "Exception and Error handling" begin
            messages = ("RecordSizeException: sample rates [9999.0, 425.0] cannot be resolved to a data record size smaller than 61440 bytes",
                        "EDFPrecisionError: String representation of value 2.0576999e7 is longer than 8 ASCII characters")
            exceptions = (OndaEDF.RecordSizeException([chunky_eeg, chunky_ecg]), OndaEDF.EDFPrecisionError(20576999.0))
            for (message, exception) in zip(messages, exceptions)
                buffer = IOBuffer()
                showerror(buffer, exception)
                @test String(take!(buffer)) == message
            end
        end
    end
    @testset "annotation import/export via round trip" begin
        round_tripped = edf_to_onda_annotations(exported_edf, uuid)

        @test round_tripped isa Vector{<:Onda.Annotation}
        # annotations are sorted by start time on export
        ann_sorted = sort(annotations; by=row -> Onda.start(row.span))
        @test getproperty.(round_tripped, :span) == getproperty.(ann_sorted, :span)
        @test getproperty.(round_tripped, :value) == getproperty.(ann_sorted, :value)
        # same recording UUID passed as original:
        @test getproperty.(round_tripped, :recording) == getproperty.(ann_sorted, :recording)
        # new UUID for each annotation created during import
        @test all(getproperty.(round_tripped, :id) .!= getproperty.(ann_sorted, :id))
    end

    @testset "full service" begin
        # import annotations
        nt = store_edf_as_onda(exported_edf, mktempdir(), uuid)
        @test nt.annotations isa Vector{<:Onda.Annotation}
        # annotations are sorted by start time on export
        ann_sorted = sort(annotations; by=row -> Onda.start(row.span))
        @test getproperty.(nt.annotations, :span) == getproperty.(ann_sorted, :span)
        @test getproperty.(nt.annotations, :value) == getproperty.(ann_sorted, :value)
        # same recording UUID passed as original:
        @test getproperty.(nt.annotations, :recording) == getproperty.(ann_sorted, :recording)
        # new UUID for each annotation created during import
        @test all(getproperty.(nt.annotations, :id) .!= getproperty.(ann_sorted, :id))
        info_orig = first(onda_samples).info
        info_round_tripped = SamplesInfo(first(nt.signals))

        @test all(getproperty(info_orig, p) == getproperty(info_round_tripped, p) for p in propertynames(info_orig))

        # don't import annotations
        nt = store_edf_as_onda(exported_edf, mktempdir(), uuid; import_annotations=false)
        @test nt.annotations isa Vector{<:Onda.Annotation}
        @test length(nt.annotations) == 0

        # import empty annotations
        exported_edf2 = onda_to_edf(samples_to_export)
        @test_logs (:warn, r"No annotations found in") store_edf_as_onda(exported_edf2, mktempdir(), uuid; import_annotations=true)
    end

end

ext: jl


text:
module Raytracing

using Random

using GeometryBasics
using Rotations
using CoordinateTransformations
using StaticArrays
using LinearAlgebra

using FiniteDifferences
using ForwardDiff
using Zygote

using CUDA

using ColorTypes
using ColorVectorSpace
using Images
using ImageIO

include("ray.jl")
include("camera.jl")
include("sdf.jl")
include("utils.jl")
include("material.jl")
include("shape.jl")
include("scene.jl")
include("raytrace.jl")

end


text:
"""
`optimalPLD_createPath(machine, typeofpath, varargin...)`

Function to create a path to a particular folder; returns a string with the specified path.

`machine` = integer specifying computer where code is running. A value of 2 means to use the immediate subdirectory of the current working directory.

`typeofpath` = string specifying whether to create path to results, data, etc.

`varargin` = additional arguments that specify the directory (see code definition).

# Example

If the current directory is `/Users/example/Documents/optimalPLD/optimalPLD_exec_scripts/`, then

`optimalPLD_createPath(2, "data_dir", "testDataSet")`

will return

`/Users/example/Documents/optimalPLD/optimalPLD_data/testDataSet/`

"""
function optimalPLD_createPath(machine, typeofpath, varargin...)
	if (machine == 1)
		# VVM's macbook pro
		root = "/Users/vvmisic/Documents/optimalPLD/"
	elseif (machine == 2)
		# Catch-all -- go one directory up from the directory where .jl script
		# was executed.
		root = "../"
	else
		error("Unknown machine argument; use either 1 or 2")
	end

	if (typeofpath == "exp_dir")
		expname = varargin[1]
		path = "$root/optimalPLD_testing/$expname/";
	elseif (typeofpath == "data_dir")
		dataname = varargin[1]
		path = "$root/optimalPLD_data/$dataname/";
	else
		error("Unknown typeofpath! See optimalPLD_createPath code.")
	end

	return path

end

max_stars_repo_name: vvmisic/EFCPLO

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: vvmisic/EFCPLO

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: vvmisic/EFCPLO

max_forks_repo_licenses: MIT Unlicense


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module SeriesExpansion
using Nemo

export ModuleSeriesExpansion
export Coefficients, G000045, G000257, L000257
export G000032, L000032, G000073, L000073, G000108, L000108, G000957, L000957
export G001003, L001003, G001006, L001006, G001045, L001045, G002426, L002426
export G005043, L005043, G006318, G068875, L068875

"""
The generating functions of various combinatorial and number-theoretic functions.

* Coefficients, G000045, G000257, L000257, G000032, L000032, G000073, L000073, G000108, L000108, G000957, L000957, G001003, L001003, G001006, L001006, G001045, L001045, G002426, L002426, G005043, L005043, G006318, G068875, L068875
"""
const ModuleSeriesExpansion = ""

"""
Return the list of coefficients of the power series s.
"""
function Coefficients(s, len)
    R, x = PowerSeriesRing(ZZ, len + 2, "x")
    ser = s(x)
    [coeff(ser, k) for k in 0:len-1]
end

"""
The generating function of the Lucas numbers.
"""
G000032(x) = 1 + divexact(x*(1 + 2x), 1 - x - x^2)

"""
Return a list of Lucas numbers.
"""
L000032(n) = Coefficients(G000032, n)

"""
The generating function of the Fibonacci numbers.
"""
G000045(x) = divexact(x, 1 - x - x^2)

"""
The generating function of the Tribonacci numbers.
"""
G000073(x) = inv(1 - x - x^2 - x^3)

"""
Return a list of Tribonacci numbers.
"""
L000073(n) = Coefficients(G000073, n)

"""
The generating function of the Catalan numbers.
"""
G000108(x) = divexact(1 - sqrt(1 - 4x), 2x)

"""
Return a list of Catalan numbers.
"""
L000108(n) = Coefficients(G000108, n)

"""
The generating function of the number of rooted bicubic maps.
"""
G000257(x) = divexact(sqrt((1 - 8x)^3) + 8x^2 + 12x - 1, 32x^2)

"""
Return a list of the number of rooted bicubic maps.
"""
L000257(n) = Coefficients(G000257, n)

"""
The generating function of the Fine numbers (with a(0) = 1).
"""
G000957(x) = 1 + divexact(1 - sqrt(1 - 4x), 3 - sqrt(1 - 4x))

"""
Return a list of Fine numbers.
"""
L000957(n) = Coefficients(G000957, n)

"""
The generating function of the little Schröder numbers.
"""
G001003(x) = divexact(1 + x - sqrt(1 - 6x + x^2), 4x)

"""
Return a list of little Schröder numbers.
"""
L001003(n) = Coefficients(G001003, n)

"""
The generating function of the Motzkin numbers.
"""
G001006(x) = divexact(1 - x - sqrt(1 - 2x - 3x^2), 2x^2)

"""
Return a list of Motzkin numbers.
"""
L001006(n) = Coefficients(G001006, n)

"""
The generating function of the Jacobsthal numbers (with a(0) = 1).
"""
G001045(x) = divexact(2x^2 - 1, (x + 1)*(2x - 1))

"""
Return a list of Jacobsthal numbers.
"""
L001045(n) = Coefficients(G001045, n)

"""
The generating function of the central trinomial.
"""
G002426(x) = inv(sqrt(1 - 2x - 3x^2))

"""
Return a list of the central trinomials.
"""
L002426(n) = Coefficients(G002426, n)

"""
The generating function of the Riordan numbers with 1 prepended.
"""
G005043(x) = 1 + divexact(2x , 1 + x + sqrt(1 - 2x - 3x^2))

"""
Return a list of the Riordan numbers (1 prepended).
"""
L005043(n) = Coefficients(G005043, n)

"""
The generating function of the large Schröder numbers.
"""
G006318(x) = divexact(1 - x - sqrt(1 - 6x + x^2), 2x)

# already defined in module 'SelfConvolutive'.
#L006318(n) = Coefficients(G006318, n)

"""
The generating function of twice the Catalan numbers.
"""
G068875(x) = shift_right(1 - x - sqrt(1 - 4x), 1)

"""
Return a list of twice the Catalan numbers.
"""
L068875(n) = Coefficients(G068875, n)

#START-TEST-########################################################

using Test, SeqUtils

function test()
    @testset "SeriesCoefficients" begin

        a = L000032(14)
        b = [1, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, 199, 322, 521]
        @test all(a .== b)

        a = L000073(15)
        b = [1, 1, 2, 4, 7, 13, 24, 44, 81, 149, 274, 504, 927, 1705, 3136]
        @test all(a .== b)

        a = L000108(13)
        b = [1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786, 208012]
        @test all(a .== b)

        a = L000257(12)
        b = [1, 1, 3, 12, 56, 288, 1584, 9152, 54912, 339456, 2149888, 13891584]
        @test all(a .== b)

        a = L000957(13)
        b = [1, 1, 0, 1, 2, 6, 18, 57, 186, 622, 2120, 7338, 25724]
        @test all(a .== b)

        a = L001003(12)
        b = [1, 1, 3, 11, 45, 197, 903, 4279, 20793, 103049, 518859, 2646723]
        @test all(a .== b)

        a = L001006(14)
        b = [1, 1, 2, 4, 9, 21, 51, 127, 323, 835, 2188, 5798, 15511, 41835]
        @test all(a .== b)

        a = L001045(13)
        b = [1, 1, 1, 3, 5, 11, 21, 43, 85, 171,341, 683, 1365]
        @test all(a .== b)

        a = L002426(12)
        b = [1, 1, 3, 7, 19, 51, 141, 393, 1107, 3139, 8953, 25653]
        @test all(a .== b)

        a = L005043(12)
        b = [1, 1, 0, 1, 1, 3, 6, 15, 36, 91, 232, 603]
        @test all(a .== b)

        xL006318(n) = Coefficients(G006318, n)
        a = xL006318(12)
        b = [1, 2, 6, 22, 90, 394, 1806, 8558, 41586, 206098, 1037718, 5293446]
        @test all(a .== b)

        a = L068875(10)
        b = [1, 2, 4, 10, 28, 84, 264, 858, 2860, 9724]
        @test all(a .== b)
    end
end

function demo()
    L068875(12) |> println
end

"""
L068875(1000)
    0.025675 seconds (1.03 k allocations: 24.813 KiB)
"""
function perf()
    @time L068875(1000)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module


text:
include("C:\\Users\\daniel.coutinho\\github\\DSGE-HMC\\src\\priors.jl")

using Distributions
using Plots

x = randn(100)

y = 0.5*x + randn(100)

prior = Normal(0,1)

par = zeros(10000)

likelihood(y,x,par) = sum(logpdf.(Normal(0,1),y - x*par))

unif = Uniform(0,1)

for i in 2:10000
    kernel = Normal(par[i-1],0.5)
    novo_par = rand(kernel)
    kernel_novo = Normal(novo_par,0.5)

    num = likelihood(y,x,novo_par) + logpdf(prior,novo_par) #+ logpdf(kernel_novo,par[i-1])

    dem = likelihood(y,x,par[i-1]) + logpdf(prior,par[i-1]) #+ logpdf(kernel,novo_par)
    alpha = exp(min(0,num - dem))
    p = rand(unif)
    if alpha < p
        par[i] = par[i-1]
    else
        par[i] = novo_par
    end
end

histogram(par)

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
# Utilities for scan matching

export overlayScanMatcher

# this function uses the sum of squared differences between the two images.
# To use low-passed versions of the images, simpy set the kernel argument to
# Kernel.gaussian(10) (or an appropriate size)
getMismatch(a,b) = sqrt(sum((a-b).^2))

# sqrt(sum( imfilter( a.-b, Kernel.gaussian(5)).^2 ))
# sqrt(sum((imfilter(a, kernel).-imfilter(b, kernel)).^2))

# also do MMD version

# Next step is to define a function that applies a transform to the image. This
# transform consists of a translation and a rotation
function transformImage_SE2(img::AbstractMatrix,
                            trans::AbstractVector{<:Real},
                            rot::Real,
                            gridscale::Real=1  )
    #
    tf_ = Translation( gridscale.*trans... )
    img_r = imrotate(img, rot)
    tf_img = warp(img_r, tf_, degree=ImageTransformations.Constant())

    # replace NaN w/ 0
    mask = findall(x->isnan(x), tf_img)
    tf_img[mask] .= 0.0
    return tf_img
end

# now we can combine the two into an evaluation function
function evaluateTransform( a::AbstractMatrix,
                            b::AbstractMatrix,
                            trans::AbstractVector{<:Real},
                            rot::Real,
                            gridscale::Real=1 )
    #
    # transform image
    bp = transformImage_SE2(b, trans, rot, gridscale)

    # get matching padded views
    ap, bpp = paddedviews(0.0, a, bp)
    return getMismatch(ap,bpp)
end

#


text:
Pkg.add("FactCheck")
Pkg.add("GR")
Pkg.add("ProfileView")
Pkg.add("Formatting")
Pkg.checkout("FactCheck")

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
module ExactDiagonalization

using Reexport: @reexport

include("FED.jl")
include("FEDW.jl")
include("Perturbation.jl")
include("RPA.jl")

@reexport using .FED
@reexport using .FEDW
@reexport using .Perturbation
@reexport using .RPA
end # module


text:
function FindSol(Y,A,M,Mt)
    return A-Mt*(M*A-Y)
end

function InitSolution(MeasurementMatrix,MeasuredOutput,m)

    #Guess a random solution
    GuessedInput=zeros(m)
    #round., convert, and ones changed to v1.1 compatibility
    GuessedInput[round.(Int, rand(convert(Int, m/10))*(m-1)+ones(Int, convert(Int, m/10)))]=rand(convert(Int, m/10))

    #project the guessed solution into the space of exact answers
    GuessedInput=FindSol(MeasuredOutput,
        GuessedInput,
        MeasurementMatrix,
        pinv(MeasurementMatrix))

    return GuessedInput
end

function plusOp(x)
    max(x,0)
end

function soft(x,l)
    sign(x).*plusOp(abs(x)-l)
end


text:
module Jirachi

const VecI = AbstractVector; # Input  Vector
const VecO = AbstractVector; # Output Vector
const VecB = AbstractVector; # Buffer Vector
const MatI = AbstractMatrix; # Input  Matrix
const MatO = AbstractMatrix; # Output Matrix
const MatB = AbstractMatrix; # Buffer Matrix

const 𝚷 = 2.0 * π

fcall(fn::Function, x::VecI) = fn(x)

include("./BLAS/BLAS_Lv1.jl")
include("./BLAS/BLAS_Lv2.jl")
include("./DiffEq/RungeKutta.jl")
include("./DiffEq/BulirschStoer.jl")
include("./interpolation.jl")
include("./fft.jl")
include("./stats.jl")
include("./sorting.jl")
include("./QLearning.jl")
include("./optimizer/minimizer.jl")

end # module


text:
export facecurrents
export potential
export get_scatter_parameters

"""
	fcr, geo = facecurrents(coeffs, basis)

Compute the value of the function with the given collection of coeffient in the provided basis in all the centroids of the mesh underlying the basis. The mesh is returned together with the currents.
"""
function facecurrents(coeffs, basis)

	T = eltype(coeffs)
	RT = real(T)

	refs = refspace(basis)
	numrefs = numfunctions(refs)

	cells, tad = assemblydata(basis)

	mesh = geometry(basis)
	D = dimension(mesh)
	U = D+1

	# TODO: express relative to input types
	PT = SVector{U, T}
	fcr = zeros(PT, numcells(mesh))

	for (t,cell) in enumerate(cells)

		mps = neighborhood(cell, ones(RT,D)/(D+1))
        vals = refs(mps)

		# assemble in the right hand side vector
		for i in 1 : numrefs
			fx = vals[i][1]
            for (m,a) in tad[t,i]
				fcr[t] += coeffs[m] * a * fx
			end
		end
	end

	return fcr, geometry(basis)
end

function facecurrents(u, X1, Xs...)

    offset = 0

    n = numfunctions(X1)
    fcrs = [ facecurrents(u[offset + (1:n)], X1)[1] ]
    offset += n

    for i in 1:length(Xs)
        n = numfunctions(Xs[i])
        push!(fcrs, facecurrents(u[offset + (1:n)], Xs[i])[1])
        offset += n
    end

    fcrs
end

function facecurrents(u, X::DirectProductSpace)

	if length(X.factors) == 1
		return facecurrents(u, X.factors[1]), geometry(X.factors[1])[1]
	end

	fcrs = facecurrents(u, X.factors...)
	#fcr = append!(fcrs...)
	fcr = vcat(fcrs...)
	m = weld([geometry(x) for x in X.factors]...)

	fcr, m
end

function potential(op, points, coeffs, basis)
	T = SVector{3,Complex128}
	ff = zeros(T, length(points))
	store(v,m,n) = (ff[m] += v*coeffs[n])
	potential!(store, op, points, basis)
	return ff
end

function potential(op, points, coeffs, space::DirectProductSpace)
	T = SVector{3,Complex128}
	ff = zeros(T, length(points))

	@assert length(coeffs) == numfunctions(space)

	offset = 0
	for fct in space.factors
		store(v,m,n) = (ff[m] += v*coeffs[offset+n])
		potential!(store, op, points, fct)
		offset += numfunctions(fct)
	end

	ff
end

function potential!(store, op, points, basis)

	T = SVector{3,Complex128}
	z = zeros(T,length(points))

	els, ad = assemblydata(basis)
	rs = refspace(basis)

	zlocal = Array{T}(numfunctions(rs))
	qdata = quaddata(op,rs,els)

	#println("Computing nearfield.")
	print("dots out of 10: ")

	todo, done, pctg = length(points), 0, 0

	for (p,y) in enumerate(points)
		for (q,el) in enumerate(els)

			fill!(zlocal,zero(T))
			qr = quadrule(op,rs,p,y,q,el,qdata)
			farfieldlocal!(zlocal,op,rs,y,el,qr)

			# assemble from local contributions
			for (r,z) in enumerate(zlocal)
                for (n,b) in ad[q,r]
					store(z*b,p,n)
				end
			end
		end

		done += 1
		new_pctg = round(Int, done / todo * 100)
		if new_pctg > pctg + 9
				#println(todo," ",done," ",new_pctg)
				print(".")
				pctg = new_pctg
		end
	end

	println("")

end

function farfieldlocal!(zlocal,op,refspace,y,el,qr)

    for q in qr
        x = q.point
        F = q.value
        dx = q.weight

        krn = kernelvals(op, y, x)
        for r in 1 : length(zlocal)
            zlocal[r] += integrand(op,krn,y,F[r],x) * dx
        end

    end

end


text:
"""
`stepEuler!`: Perform an in place Euler step

### Fields
* `u` - Initial state
* `∇V` - Gradient of energy
* `Δt` - Time step
"""
function stepEuler!(u, ∇V::TGV, Δt) where {TGV}
    gradV = ∇V(u);
    @. u = u - Δt *gradV;

    u
end

"""
`stepRK4!`: Perform an in place RK4 step

### Fields
* `u` - Initial state
* `∇V` - Gradient of energy
* `Δt` - Time step
"""
function stepRK4!(u, ∇V::TGV, Δt) where {TGV}

    gradV1 = ∇V(u);
    gradV2 = ∇V(u - 0.5 * Δt * gradV1);
    gradV3 = ∇V(u - 0.5 * Δt * gradV2);
    gradV4 = ∇V(u - Δt * gradV3);

    @. u = u - Δt/6 * (gradV1 + 2 * gradV2 + 2 * gradV3 + gradV4);

    u

end

"""
`stepPCEuler!`: Perform an in place Euler step with preconditioning

### Fields
* `u` - Initial state
* `∇V` - Gradient of energy
* `P` - Preconditioner
* `Δt` - Time step
"""
function stepPCEuler!(u, ∇V::TGV, P::TP, Δt) where {TGV, TP}

    gradV = ∇V(u);
    Pu = P(u);
    s = Pu\gradV;
    @. u = u - Δt * s;

    u
end


text:
using LearnBase
"""
$(SIGNATURES)

Mostly useful for generating data for training offline.
"""
function batch_rollout(env::Env, init_xs, uss, T, B, dt)
    @assert size(init_xs)[end] == B
    xs = zeros(Float64, env.sys.n, T, B)
    @views xs[:, 1, :] .= init_xs
    for b in 1:B
        for k = 1:T-1
            @views TrajectoryOptimization.evaluate!(xs[:, k+1, b], env.sys, xs[:, k, b], uss[:, k, b], dt)
            @views xs[:, k+1, b] .= wraptopi(xs[:, k+1, b], env.thetamask)
        end
    end
    return xs, uss
end

"""
Rollout a trajectory using a policy in an environment for `T` timesteps.

$(SIGNATURES)
"""
function rollout(env::Env, policy, T, dt)
    xs = zeros(Float64, env.sys.n, T)
    us = zeros(Float64, env.sys.m, T-1)
    @views xs[:, 1] = env.x0
    for t in 1:T-1
        @views us[:, t] = policy(convert(Float64, t)*dt, xs[:, t], env.thetamask)
        @views TrajectoryOptimization.evaluate!(xs[:, t+1], env.sys, xs[:, t], us[:, t], dt)
        @views xs[:, t+1] .=  xs[:, t+1] # wraptopi(xs[:, t+1], env.thetamask)
    end
    return xs, us
end

"""
Converts a batch of trajectories of position, velocities and input forces into a Dataset

$(SIGNATURES)
"""
function trajbatch2dataset(xs, uss)
    A = slidingwindow(i->i+1, xs, 1, LearnBase.ObsDim.Constant(2))
    newxs = []
    newxsp = []
    newus = []
    for (i, (x, xp)) in enumerate(A)
        u = uss[:, i, :]
        push!(newus, u)
        push!(newxs, dropdims(x, dims=2))
        push!(newxsp, xp)
    end
    return Dataset(reduce(hcat, newxs), reduce(hcat, newus), reduce(hcat, newxsp))
end


text:
################################################################################
#                                 MvNormal                                     #
################################################################################

# gradients verified empirically using PyTorch version

function ∇logpdf_mu(dy, d::MvNormal, x)
    # device = iscuarray(x) ?
    iΣ = inv(d.Σ)
    iΣf = Matrix(iΣ)  # TODO: make device agnostic
    return dy .* iΣf * (x .- d.μ)
end

function ∇logpdf_sigma(dy, d::MvNormal, x)
    iΣ = inv(d.Σ)
    iΣf = Matrix(iΣ)  # TODO: make device agnostic
    diff = x .- d.μ
    return -0.5f0 .* (iΣf .- iΣ * diff * diff' * iΣf) .* dy
end

function ∇logpdf_x(dy::AbstractVector, d::MvNormal, x::AbstractMatrix)
    ret = similar(x)
    for j in 1:size(x, 2)
        ret[:, j] = Distributions.gradlogpdf(d, @view x[:, j]) .* dy[j]
    end
    return ret
end

function ∇logpdf_x(dy::Real, d::MvNormal, x::AbstractVector)
    return Distributions.gradlogpdf(d, x) .* dy
end

@ctor MvNormal(μ, Σ)
@diffrule logpdf(_d::MvNormal, x) _d.μ ∇logpdf_mu(dy, _d, x)
@diffrule logpdf(_d::MvNormal, x) _d.Σ ∇logpdf_sigma(dy, _d, x)
@diffrule logpdf(_d::MvNormal, x) x ∇logpdf_x(dy, _d, x)


text:
## GPU support
import GPUArrays: backend, gpu_call, linear_index, launch_heuristic, launch_configuration,
    AbstractGPUArrayStyle

# check the backend
import ArrayInterface: GPU, device, parent_type
devices(::Type{<:NamedTuple{<:Any,T}}) where {T} = devices(T)
devices(::Type{T}) where {T<:Tuple} = getsame(device, T.parameters...)
backends(::Type{<:NamedTuple{<:Any,T}}) where {T} = backends(T)
backends(::Type{T}) where {T<:Tuple} = getsame(backend, T.parameters...)

# Adapt fix
import Adapt: adapt_structure, adapt
@inline adapts(to, x, ys...) = (adapt(to, x), adapts(to, ys...)...)
@inline adapts(to) = ()

device(::Type{T}) where {T<:TupleDummy} = devices(T.parameters[4])
backend(::Type{T}) where {T<:TupleDummy} = backends(T.parameters[4])
adapt_structure(to, td::TupleDummy{T,N,L}) where {T,N,L} =
    TupleDummy{T,N,L}(adapts(to, td.arrays...), td.ax)

# general gpu_copyto!(modified from GPUArrays.jl's implement to support OffsetArrays.jl)

@inline function gpu_copyto!(dest::AbstractArray, bc::Broadcasted{Nothing})
    axes(dest) == axes(bc) || throwdm(axes(dest), axes(bc))
    isempty(dest) && return dest
    bc′ = preprocess(dest, bc)

    function broadcast_kernel(ctx, dest, bc, nelem)
        Inds = eachindex(bc)
        offset, Iᵉ = firstindex(Inds) - 1, lastindex(Inds)
        @inbounds for i = 1:nelem
            j = linear_index(ctx, i) + offset
            j > Iᵉ && return
            I = Inds[j]
            dest[I] = bc[I]
        end
    end

    elements = length(dest)
    elements_per_thread = typemax(Int)
    heuristic = launch_heuristic(backend(dest), broadcast_kernel, dest, bc′, 1; elements, elements_per_thread)
    config = launch_configuration(backend(dest), heuristic; elements, elements_per_thread)
    gpu_call(broadcast_kernel, dest, bc′, config.elements_per_thread;
             threads = config.threads, blocks = config.blocks)
    dest
end
## use style to dispatch
function dispatch_style(Style, need_more = false)
    @eval @inline mtb_materialize(bc::Broadcasted{<:$Style}, ::Integer) =
        copy(instantiate(bc))
    @eval @inline mtab_materialize(bc::Broadcasted{<:$Style}, ::Integer) =
        tab_copy(instantiate(bc))
    @eval @inline mtb_materialize!(::$Style, dest, bc::Broadcasted{Style}, ::Integer) where {Style} =
        gpu_copyto!(dest, instantiate(Broadcasted{Nothing}(bc.f, bc.args, axes(dest))))
    @eval @inline tab_copy(bc::Broadcasted{<:$Style}) =
        gpu_copyto!(toa_similar(bc), convert(Broadcasted{Nothing}, bc)) |> parent
end
dispatch_style(:AbstractGPUArrayStyle)
@inline tab_copy(bc::Broadcasted{<:AbstractGPUArrayStyle{0}}) = copy(bc)

## AbstractWrapper
function map_show_copy(WrapperType::Symbol)

    @eval @inline copyto!(dest::$WrapperType, bc::Broadcasted{Nothing}) = begin
        device(dest) isa GPU && return gpu_copyto!(dest, bc)
        invoke(copyto!, Tuple{AbstractArray,Broadcasted{Nothing}}, dest, bc)
    end

    @eval BroadcastStyle(::Type{Base.RefValue{AT}}) where {AT<:$WrapperType} =
        BroadcastStyle(AT) |> forcedim0
end

forcedim0(x) = x
forcedim0(::Style) where {Style<:AbstractArrayStyle} = Val(0) |> Style

@require OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" begin
    import .OffsetArrays: OffsetArray
    ## general
    ## device has been defined in ArrayInterface
    backend(::Type{T}) where {T<:OffsetArray} = parent_type(T) |> backend
    ## adapt_structure has been defined in OffsetArrays.jl
    map_show_copy(:OffsetArray)
    ## unique
    Base.collect(A::OffsetArray) = collect(parent(A))
    BroadcastStyle(::Type{OA}) where {OA<:OffsetArray} = parent_type(OA) |> BroadcastStyle
end

@require StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" begin
    import .StructArrays: StructArray, StructArrayStyle, components
    ## general
    device(::Type{T}) where {T<:StructArray} = T.parameters[3] |> devices
    backend(::Type{T}) where {T<:StructArray} = T.parameters[3] |> backends
    ## adapt_structure has been defined in StructArrays.jl
    map_show_copy(:StructArray)
    # unique
    dispatch_style(:(StructArrayStyle{<:AbstractGPUArrayStyle}))
    forcedim0(::StructArrayStyle{Style}) where {Style} = StructArrayStyle{typeof(forcedim0(Style()))}()

    function Base.similar(bc::Broadcasted{StructArrayStyle{S}}, ::Type{ElType}) where {S,ElType}
        bc′ = convert(Broadcasted{S}, bc)
        if isstructtype(ElType)
            return StructArrays.buildfromschema(T -> similar(bc′, T), ElType)
        else
            return similar(bc′, ElType)
        end
    end
end


text:
@testset "Utilities" begin
    @testset "Crossover/Crossunder" begin
        Random.seed!(SEED)
        x = cumsum(randn(N)) .+ X0
        y = x + randn(N)
        cxo = crossover(x, y)
        cxu = crossunder(x, y)
        @test any(cxo)
        @test any(cxu)
        @test !any(cxo .* cxu)  # ensure crossovers and crossunders never coincide
    end
end


text:
using DiffEqBase, StochasticDiffEq, DiffEqDevTools, Test, LinearAlgebra, Random

u0 = ones(2)
A = [-3/2 1/20
      1/20 -3/2]
B = [1/5 1/100
    1/100 1/5]

function f(du,u,p,t)
  mul!(du,A,u)
end
function σ(du,u,p,t)
  mul!(@view(du[:,1]),B,u)
  mul!(@view(du[:,2]),B,u)
end

function (::typeof(f))(::Type{Val{:analytic}},u0,p,t,W)
  tmp = (A-(B^2))*t + B*W[1] + B*W[2]
  exp(tmp)*u0
end

prob2 = SDEProblem(f,σ,u0,(0.0,1.0),noise_rate_prototype=rand(2,2))

sol2 = solve(prob2,EM(),dt=1/100)

dts = 1 ./ 2 .^ (14:-1:7)

println("First Test")
Random.seed!(100)
sim  = test_convergence(dts,prob2,EM(),numMonte=100)
@test abs(sim.𝒪est[:l2]-0.5) < 0.1

# using Plots; plot(sim)

u0 = rand(2)
A = [2.0 -1.0; -1.0 2.0]
B = [1/5 1/100
    1/100 1/5]

function f(du,u,p,t)
  mul!(du,A,u)
  du .+= 1.01u
end
function σ(du,u,p,t)
  mul!(@view(du[:,1]),B,u)
  mul!(@view(du[:,2]),B,u)
end

function (::typeof(f))(::Type{Val{:analytic}},u0,p,t,W)
 tmp = (A+1.01I-(B^2))*t + B*sum(W)
 exp(tmp)*u0
end

prob2 = SDEProblem(f,σ,u0,(0.0,1.0),noise_rate_prototype=rand(2,2))

sol2 = solve(prob2,EM(),dt=1/100)

dts = 1 ./ 2 .^ (14:-1:7)

println("Second Test")
Random.seed!(100)
sim  = test_convergence(dts,prob2,EM(),numMonte=50)
# Superconvergence
@test abs(sim.𝒪est[:l2]-1.0) < 0.1

# using Plots; plot(sim)


text:
"""
    get_spatial_profile_single_mode_and_fit(problem, mode_index)

return (data = ( x=DCM, y=ψ² ), fit=(x=x_fit, y=y_fit) )
"""
function get_spatial_profile_single_mode_and_fit(problem, mode_index)
    DCM, ψ² = get_spatial_profile_single_mode(problem, mode_index)
    x_fit, y_fit = get_decay_fit(problem, mode_index)
    return (DCM = DCM, ψ² = ψ², fit = (x = x_fit, y = y_fit))
end

function get_decay_fit(problem, mode_index)
    DCM, ψ² = get_spatial_profile_single_mode(problem, mode_index)

    x_use, y_use = select_points(DCM, ψ²)
    A_L₁, ξ_L₁, R1, y_fit_L₁ = get_coeffs_L1(x_use, y_use)

    if is_R1_localized(R1)
        return x_use, y_fit_L₁
    else
        A_L₂, ξ_L₂, R2, y_fit_L₂ = get_coeffs_L2(x_use, y_use)
        return x_use, y_fit_L₂
    end
end


text:
struct Satellite <: AbstractModel
    J::Diagonal{Float64,SVector{3,Float64}}
end

Satellite() = Satellite(Diagonal(@SVector ones(3)))

Base.size(::Satellite) = 7,3
Base.position(::Satellite, x::SVector) = @SVector zeros(3)
orientation(::Satellite, x::SVector) = @SVector [x[4], x[5], x[6], x[7]]

function dynamics(model::Satellite, x::SVector, u::SVector)
    ω = @SVector [x[1], x[2], x[3]]
    q = normalize(@SVector [x[4], x[5], x[6], x[7]])
    J = model.J

    ωdot = J\(u - ω × (J*ω))
    qdot = 0.5*Lmult(q)*Vmat()'ω
    return [ωdot; qdot]
end

function state_diff(model::Satellite, x::SVector, x0::SVector)::SVector{6}
    ω = @SVector [x[1], x[2], x[3]]
    q = @SVector [x[4], x[5], x[6], x[7]]
    ω0 = @SVector [x0[1], x0[2], x0[3]]
    q0 = @SVector [x0[4], x0[5], x0[6], x0[7]]

    δω = ω - ω0
    δq = Lmult(q0)'q
    ϕ = @SVector [δq[2]/δq[1], δq[3]/δq[1], δq[4]/δq[1]]
    return [δω; ϕ]
end

function state_diff_jacobian(model::Satellite, x::SVector)
    q = @SVector [x[4], x[5], x[6], x[7]]
    G = Lmult(q)*Vmat()'
    return @SMatrix [1 0 0 0 0 0;
                     0 1 0 0 0 0;
                     0 0 1 0 0 0;
                     0 0 0 G[1] G[5] G[ 9];
                     0 0 0 G[2] G[6] G[10];
                     0 0 0 G[3] G[7] G[11];
                     0 0 0 G[4] G[8] G[12];
                     ]
end

state_diff_size(::Satellite) = 6

function TrajectoryOptimization.state_diff_jacobian!(G, model::Satellite, Z::Traj)
    for k in eachindex(Z)
        G[k] = state_diff_jacobian(model, state(Z[k]))
    end
end

struct Satellite2{R,B} <: RigidBody{R}
    J::Diagonal{Float64,SVector{3,Float64}}
end

Satellite2() = Satellite2{UnitQuaternion{Float64,CayleyMap},true}(Diagonal(@SVector ones(3)))
Satellite2{R}(;use_rot=true) where R= Satellite2{R,use_rot}(Diagonal(@SVector ones(3)))

Base.size(::Satellite2{<:UnitQuaternion}) = 7,3
Base.size(::Satellite2) = 6,3
Base.position(::Satellite2, x::SVector, unitnorm=true) = @SVector zeros(3)
function orientation(::Satellite2{UnitQuaternion{T,D}}, x::SVector{N,T2},
        unitnorm=true) where {T,T2,D,N}
    q = UnitQuaternion{T2,D}(x[4], x[5], x[6], x[7])
    if unitnorm
        return normalize(q)
    end
    q
end
orientation(::Satellite2{R}, x::SVector{N,T2}, unitnorm=false) where {R,N,T2} =
    R(T2,x[4],x[5],x[6])

function build_state(model::Satellite2{R}, ω, q) where R
    q = SVector(R(q))
    return [ω; q]
end

function Base.rand(::Satellite2{R}) where R
    u = @SVector rand(3)
    ω = @SVector rand(3)
    q = rand(R)
    return [ω; SVector(q)], u
end

function dynamics(model::Satellite2{R}, x::SVector, u::SVector) where R
    ω = @SVector [x[1], x[2], x[3]]
    q = orientation(model, x)
    J = model.J

    ωdot = J\(u - ω × (J*ω))
    qdot = kinematics(q,ω)
    return [ωdot; qdot]
end

function state_diff(model::Satellite2, x::SVector, x0::SVector)::SVector{6}
    ω = @SVector [x[1], x[2], x[3]]
    ω0 = @SVector [x0[1], x0[2], x0[3]]
    q = orientation(model, x)
    q0 = orientation(model, x0)

    δω = ω - ω0
    ϕ = q ⊖ q0
    return [δω; ϕ]
end

function state_diff(model::Satellite2{R,false}, x::SVector, x0::SVector) where R
    return x - x0
end

function state_diff_jacobian(model::Satellite2{<:UnitQuaternion}, x::SVector)
    q = orientation(model, x)
    G = TrajectoryOptimization.∇differential(q)
    return @SMatrix [1 0 0 0 0 0;
                     0 1 0 0 0 0;
                     0 0 1 0 0 0;
                     0 0 0 G[1] G[5] G[ 9];
                     0 0 0 G[2] G[6] G[10];
                     0 0 0 G[3] G[7] G[11];
                     0 0 0 G[4] G[8] G[12];
                     ]
end

function state_diff_jacobian(model::Satellite2, x::SVector)
    q = orientation(model, x)
    G = TrajectoryOptimization.∇differential(q)
    return @SMatrix [1 0 0 0 0 0;
                     0 1 0 0 0 0;
                     0 0 1 0 0 0;
                     0 0 0 G[1] G[4] G[7];
                     0 0 0 G[2] G[5] G[8];
                     0 0 0 G[3] G[6] G[9];
                     ]
end

function state_diff_jacobian(model::Satellite2{R,false}, x::SVector) where R
    return I
end

import TrajectoryOptimization: ∇²differential
function ∇²differential(model::Satellite2, x::SVector, dx::SVector)
    q = orientation(model, x)
    dq = SVector(orientation(model, dx, false))  # gradient wrt the rotation
    G2 = ∇²differential(q, dq)
    return @SMatrix [0 0 0 0 0 0;
                     0 0 0 0 0 0;
                     0 0 0 0 0 0;
                     0 0 0 G2[1] G2[4] G2[7];
                     0 0 0 G2[2] G2[5] G2[8];
                     0 0 0 G2[3] G2[6] G2[9];
                     ]
end

function ∇²differential(model::Satellite2{R,false}, x::SVector, dx::SVector) where R
    return I*0
end

state_diff_size(::Satellite2) = 6
state_diff_size(::Satellite2{<:UnitQuaternion,false}) = 7

function TrajectoryOptimization.state_diff_jacobian!(G, model::Satellite2, Z::Traj)
    for k in eachindex(Z)
        G[k] = state_diff_jacobian(model, state(Z[k]))
    end
end


text:
using Geodesy
using Test

################################################
### Helpers for testing approximate equality ###
################################################

# Interesting that this isn't in Base...
Base.isapprox(a::T, b::T; kwargs...) where {T<:Tuple} = all(ntuple(i->isapprox(a[i],b[i]), length(a)); kwargs...)

@testset "Geodesy" begin
    include("points.jl")
    include("datums.jl")
    include("transverse_mercator.jl")
    include("polar_stereographic.jl")
    include("utm.jl")
    include("transformations.jl")
    include("conversion.jl")
    include("datum_transformations.jl")
end # @testset "Geodesy"


text:
if basename(pwd()) == "aoc"
    cd("2016/8")
end

function part1(filename::AbstractString, screensize::Tuple{Int,Int})
    screen = falses(screensize)
    for instr in eachline(filename)
        A, B = parse.(Int, match(r"(\d+)(?:x| by )(\d+)", instr))
        if startswith(instr, "rect")
            screen[1:B, 1:A] .= true
        elseif startswith(instr, "rotate row")
            screen[A+1, :] = circshift(screen[A+1, :], B)
        else
            screen[:, A+1] = circshift(screen[:, A+1], B)
        end
    end
    screen
end

function prettyprint(A)
    for row in eachrow(A)
        println((b ? '#' : '.' for b in row)...)
    end
    println()
end

# part1("example.txt", (3, 7))
result = part1("input.txt", (6, 50))
count(result)
prettyprint(result)

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
function resample(uvals::NTuple{N, Union{Real, AbstractUncertainValue}}) where N
    Tuple(convert(AbstractFloat, resample(uval)) for uval in uvals)
end

function resample(uvals::NTuple{N, Union{Real, AbstractUncertainValue}}, n::Int) where N
    [resample(uvals) for i = 1:n]
end

function resample_elwise(uvals::NTuple{N, Union{Real, AbstractUncertainValue}}) where N
    [resample(uval) for uval in uvals]
end

function resample_elwise(uvals::NTuple{N, Union{Real, AbstractUncertainValue}}, n::Int) where N
    [[resample(uval) for i = 1:n] for uval in uvals]
end


text:
using MLJModelInterface
using Random

const MMI = MLJModelInterface

const PLSRegressor_Desc  = "A Partial Least Squares Regressor. Contains PLS1, PLS2 (multi target) algorithms. Can be used mainly for regression."
const KPLSRegressor_Desc = "A Kernel Partial Least Squares Regressor. A Kernel PLS2 NIPALS algorithms. Can be used mainly for regression."

const MLJDICT = Dict(:Pls1 => PLS1Model,:Pls2 => PLS2Model,:Kpls => KPLSModel)

mutable struct PLS <: MMI.Deterministic
    n_factors::Int
end

mutable struct KPLS <: MMI.Deterministic
    n_factors::Integer
    kernel::String
    width::Real
end

function PLS(; n_factors=1)
    model   = PLS(n_factors)
    message = MLJModelInterface.clean!(model)
    isempty(message) || @warn message
    return model
end

function MMI.clean!(m::PLS)
    warning = ""
    if m.n_factors <= 0
        warning *= "Parameter `n_factors` expected to be positive, resetting to 1"
        m.n_factors = 1
    end
    return warning
end

function KPLS(; n_factors=1,kernel="rbf",width=1.0)
    model   = KPLS(n_factors,kernel,width)
    message = MLJModelInterface.clean!(model)
    isempty(message) || @warn message
    return model
end

function MMI.clean!(m::KPLS)
    warning = ""
    if m.n_factors <= 0
        warning *= "Parameter `n_factors` expected to be positive, resetting to 1"
        m.n_factors = 1
    end
    if m.kernel != "rbf"
        warning *= "Parameter `kernel` expected to be rbf, resetting to rbf"
        m.kernel = "rbf"
    end

    return warning
end

function MMI.fit(m::PLS, verbosity::Int, X,Y)

    X =  MMI.matrix(X)
    Y = (MMI.scitype(Y) == Table{AbstractArray{Continuous,1}} ? MMI.matrix(Y) : Y)
    Y = (MMI.scitype(Y) == Table{AbstractArray{Continuous,2}} ? MMI.matrix(Y) : Y)

    check_constant_cols(X)
    check_constant_cols(Y)
    check_params(m.n_factors, size(X,2),"linear")
    check_data(X, Y)

    model                    = PLSModel(X,Y,m.n_factors)

    (fitresult,cache,report) = fitting(model,X,Y)

    return (fitresult,cache,report)

end

function MMI.fit(m::KPLS, verbosity::Int, X,Y)

    X = MMI.matrix(X)
    Y = (MMI.scitype(Y) == Table{AbstractArray{Continuous,1}} ? MMI.matrix(Y) : Y)
    Y = (MMI.scitype(Y) == Table{AbstractArray{Continuous,2}} ? MMI.matrix(Y) : Y)

    check_constant_cols(X)
    check_constant_cols(Y)
    check_params(m.n_factors, size(X,2),m.kernel)
    check_data(X, Y)

    model = PLSModel(X,Y,
                 m.n_factors,
                 m.kernel,
                 m.width)

    (fitresult,cache,report) = fitting(model,X,Y)

    return (fitresult,cache,report)
end

function MMI.predict(m::Union{PLS,KPLS}, fitresult, X)

    X = MMI.matrix(X)
    check_data(X,fitresult.nfeatures)

    Y =  predicting(fitresult,X)
    Y = (length(size(Y)) > 1 ? MMI.table(Y) : Y)

    return Y
end

MMI.metadata_pkg.(
    (PLS, KPLS),
    name       = "Partial Least Squares Regressor",
    uuid       = "e010f91f-06b9-52b3-bed3-bb1da186bddc",
    url        = "https://github.com/lalvim/PLSRegressor.jl",
    julia      = true,
    license    = "MIT",
    is_wrapper = false)

MMI.metadata_model(PLS,
    input   = Table(Continuous),
    target  = Union{AbstractVector{<:Continuous},Table(Continuous)},
    weights = false,
    descr   = PLSRegressor_Desc)

MMI.metadata_model(KPLS,
    input   = Table(Continuous),
    target  = Union{AbstractVector{<:Continuous},Table(Continuous)},
    weights = false,
    descr   = KPLSRegressor_Desc)


text:
path = Base.source_dir()
for i in readdir(path * "/TestHelpers")
  if endswith(i, ".jl")
    include(joinpath(path, "TestHelpers", i))
  end
end


text:
 """
 Performs a convergence study of the discrete cubed sphere surface
 using a Gauss quadrature of degree degree for numerical integration.
 order is the order of the polynomial space used to approximate the
 surface of the cubed sphere.

 Returns a tuple with three entries.
  [1] Discrete cubed sphere surfaces for the associated values of n.
  [2] Relative errors among the approximate surface and the exact one.
  [3] Slope of the log(h)-log(err) relative error convergence curve.
 """
 function convergence_discrete_cubed_sphere_surface(hs,model_args_series,degree)
   approx_surfaces=Float64[]
   rel_errors=Float64[]
   exact_surface=4π
   current=1
   for args in model_args_series
      model=CubedSphereDiscreteModel(args...)
      approx_surface=compute_discrete_cubed_sphere_surface(model,degree)
      append!(approx_surfaces,approx_surface)
      append!(rel_errors,abs(approx_surface-exact_surface)/exact_surface)
      current=current*10
   end
   return approx_surfaces,rel_errors,slope(hs,rel_errors)
 end

 function compute_discrete_cubed_sphere_surface(model,degree)
   Ω = Triangulation(model)
   dΩ = Measure(Ω,degree)
   f(x)=1
   sum(∫(f)dΩ)
 end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
using IcanProj
using NetCDF
using DataFrames
using CSV
using Statistics
using PyPlot
using Dates

# Settings

path = "/data04/jmg/fsm_simulations/netcdf/fsmres_open"

cfg = 32

variable = "latmo"

resolution = "5km"

# Load and link custom results

file_fine = joinpath(path, "results_$(cfg)", "$(variable)_1km.nc")

file_coarse = joinpath(path, "results_$(cfg)", "$(variable)_$(resolution).nc")

df_links = link_results(file_fine, file_coarse)

var_coarse, var_ref, ngrids = unify_results(file_fine, file_coarse, df_links, variable)

(any(isnan.(var_ref)) | any(isnan.(var_coarse))) ? error("nans in results") : nothing

# Load and link swes

file_fine = joinpath(path, "results_$(cfg)", "swe_1km.nc")

file_coarse = joinpath(path, "results_$(cfg)", "swe_$(resolution).nc")

swe_coarse, swe_ref, ngrids = unify_results(file_fine, file_coarse, df_links, "swe")

(any(isnan.(swe_ref)) | any(isnan.(swe_coarse))) ? error("nans in results") : nothing

# Load time

ta = ncread(file_coarse, "time_array")

t = [DateTime(ta[1,i],ta[2,i],ta[3,i],ta[4,i]) for i in 1:size(ta, 2)]

# Compute scale error

scale_error = var_coarse - var_ref

scale_rmse = sqrt.(mean(scale_error.^2, dims = 1))

print("Number of time steps: $(size(scale_error, 1))")
print("Number of grid cells: $(size(scale_error, 2))")
print("Maximum rmse: $(maximum(scale_rmse))")

# Plots

if true

  fig, axes = plt[:subplots](nrows = 3, ncols = 1)

  axes[1][:plot](t, scale_error, color = "gray")

  axes[2][:plot](t, var_ref, color = "gray")

  axes[3][:plot](t, var_coarse, color = "gray")

  #axes[1][:plot](t, cumsum(scale_error, dims = 1), color = "gray")

  # mean_coarse = [mean(var_coarse[i, ikeep[i, :]]) for i in 1:size(ikeep, 1)]

  # mean_ref = [mean(var_ref[i, ikeep[i, :]]) for i in 1:size(ikeep, 1)]

end


text:
using PowerModels
using JuMP
using CPLEX
include("prob.jl")

function get_independent_load_shed(; file = "../data/nesta_case24_ieee_rts.m", solver = IpoptSolver(), model_constructor = SOCWRPowerModel, lines = final_lines)

    independent_load_shed = [i => 0.0 for i in lines]
    for line in lines
        data = PowerModels.parse_file(file)
        for branch in data["branch"]
            if Int(branch["index"]) == line
                branch["br_status"] = 0
            end
        end
        buses = [ Int(bus["index"]) => bus for bus in data["bus"] ]
        buses = filter((i, bus) -> bus["bus_type"] != 4, buses)
        branches = [ Int(branch["index"]) => branch for branch in data["branch"] ]
        branches = filter((i, branch) -> branch["br_status"] == 1 && branch["f_bus"] in keys(buses) && branch["t_bus"] in keys(buses), branches)
        branch_indexes = collect(keys(branches))
        pm = model_constructor(data; solver = solver)
        post_pfls(pm)
        status, solve_time = solve(pm)
        result = PowerModels.build_solution(pm, status, solve_time)
        independent_load_shed[line] = result["objective"]
    end
    return independent_load_shed;
end

max_stars_repo_name: kaarthiksundar/acopf-n-k

max_issues_repo_name: kaarthiksundar/acopf-n-k

max_forks_repo_name: kaarthiksundar/acopf-n-k


text:
using Base.Filesystem
using BinDeps: unpack_cmd, download_cmd

const meshcat_sha = "b34c2668c94675aa582c5dc2024fdff052db9511"
const meshcat_url = "https://github.com/rdeits/meshcat/archive/$meshcat_sha.zip"

const assets_dir = normpath(joinpath(@__DIR__, "..", "assets"))
const meshcat_dir = joinpath(assets_dir, "meshcat")
const stamp_file = joinpath(assets_dir, "meshcat.stamp")

function update_meshcat()
    if !isdir(assets_dir)
        mkpath(assets_dir)
    end

    if isdir(joinpath(meshcat_dir, ".git"))
        info("Meshcat assets in $meshcat_dir have been cloned with git, so they will not be automatically downloaded. To force a re-download, delete or rename the directory $meshcat_dir")
        return
    end
    if isdir(meshcat_dir) && isfile(stamp_file)
        stamped_sha = strip(open(readstring, stamp_file))
        if stamped_sha == meshcat_sha
            return
        else
            info("Updating meshcat assets in $meshcat_dir from SHA $stamped_sha to $meshcat_sha")
        end
    end

    if isdir(meshcat_dir)
        rm(meshcat_dir; recursive=true)
    end

    mktempdir() do download_dir
        download_path = joinpath(download_dir, "meshcat.zip")
        run(download_cmd(meshcat_url, download_path))
        run(unpack_cmd(download_path, download_dir, ".zip", nothing))
        mv(joinpath(download_dir, "meshcat-$meshcat_sha"), meshcat_dir; remove_destination=true)
    end
    open(stamp_file, "w") do file
        print(file, meshcat_sha)
    end
end

update_meshcat()


text:
module Conversion

using DataStructures
using PhaseSpaceIO
using Accessors
using ArgCheck

function compute_extra_ints(p::EGSParticle)::NTuple{1, Int32}
    ulatch = UInt32(p.latch)
    ilatch = Int32(ulatch)
    (ilatch,)
end

function compute_extra_floats(p::EGSParticle{Nothing})::NTuple{0, Float32}
    ()
end

function compute_extra_floats(p::EGSParticle{Float32})::NTuple{1, Float32}
    (p.zlast,)
end

function to_iaea(p::EGSParticle; z)
    IAEAParticle(
        typ = particle_type(p),
        x=p.x,
        y=p.y,
        z=z,
        u=p.u,
        v=p.v,
        w=p.w,
        E=p.E,
        new_history=p.new_history,
        weight=p.weight,
        extra_floats=compute_extra_floats(p),
        extra_ints=compute_extra_ints(p),
    )
end

function to_egs(p::IAEAParticle)
    typ = p.typ
    q = if typ == photon
        0
    elseif typ == electron
        -1
    elseif typ == positron
        1
    else
        msg = "Cannot convert particle of type = $typ to EGS"
        throw(ArgumentError(msg))
    end
    EGSParticle(
        latch=Latch(charge=q),
        new_history=p.new_history,
        x=p.x,
        y=p.y,
        u=p.u,
        v=p.v,
        w=p.w,
        E=p.E,
        weight=p.weight,
    )
end

function phsp_convert(src::AbstractString, dst::IAEAPath; z)
    phsp_convert([src], dst, z=z)
end

function phsp_convert(srcs::AbstractVector{<:AbstractString}, dst::IAEAPath;  z)
    iters = map(egs_iterator, srcs)
    ret = phsp_convert(iters, dst, z=z)
    foreach(close, iters)
    ret
end

function _create_iaea_header(iters; z)
    iter = first(iters)
    originalcount = sum(iters) do iter
        iter.header.originalcount
    end
    nt = (z=z,)
    _Nf(::Type{EGSParticle{Nothing}}) = 0
    _Nf(::Type{EGSParticle{Float32}}) = 1
    Nf = _Nf(eltype(iter))
    H = IAEAHeader{1,Nf}
    d = OrderedDict(:ORIG_HISTORIES => string(originalcount))
    header = H(nt, d)
    return header
end

@noinline function phsp_append(writer, iter, z)
    for p in iter
        write(writer, to_iaea(p, z=z))
    end
end

@noinline function phsp_convert(iters, dst::IAEAPath; z)
    if isempty(iters)
        return
    end
    h = _create_iaea_header(iters, z=z)
    iaea_writer(dst, h) do w
        for iter in iters
            phsp_append(w, iter, z)
        end
    end
end

end#module


text:
module NumDisplay

using PiGPIOC

# constructors
#, DisplayBCD
export DisplaySPI, DisplayDirect, DisplayBCD

# low level
export scan_rate, decode_mode, set_intensity, set_limit, shutdown_mode_on, shutdown_mode_off,
    test_mode_on, test_mode_off, write_digit,
    update

# high level
export write_number, write_symbols

include("decoding.jl")
include("extend.jl")
include("abstract-num-display.jl")
include("display-no-buffer.jl")
include("display-direct.jl")
include("display-bcd.jl")
include("display-spi.jl")

end # module


text:
### A Pluto.jl notebook ###
# v0.12.7

using Markdown
using InteractiveUtils

# ╔═╡ 5d84ea34-2202-11eb-0a5e-8153f59b5785
using Gumbo, HTTP,PyPlot,DataFrames,Dates,CSV

# ╔═╡ d02e25a2-2202-11eb-21b7-89a516182e0f
rawpage=HTTP.get("https://de.wikipedia.org/wiki/COVID-19-Pandemie_in_Deutschland/Statistik");

# ╔═╡ fb46d6f8-2202-11eb-1cda-8509c49cb0de
parsed_page=parsehtml(String(rawpage.body));

# ╔═╡ ce05b5c6-2205-11eb-3533-b3a8e0450f2e
function find_table(element,text)
	if isa(element,HTMLText)
     	return nothing
    end
	if length(element.children)>0 && element[1] isa Gumbo.HTMLElement{:caption} && occursin(text,String(element[1][1].text))
		return element
	else
		for i=1:length(element.children)
			 result=find_table(element.children[i],text)
             if result!=nothing
				return result
			end
		end
	end
end

# ╔═╡ aaf63af0-2206-11eb-0c77-df3ffa48b31b
t=find_table(parsed_page.root,"Infektionsfälle");

# ╔═╡ 05f3b3a0-2208-11eb-07f4-cbfd72d3aa4a
function header(t)
	hd=t[2][1]
	hdr=[hd[i][1][1].text for i=3:18]
	push!(hdr,"DE")
end

# ╔═╡ 226610d0-220a-11eb-1a51-e913fa1a9120
hd=header(t)

# ╔═╡ f85174da-220c-11eb-3bdc-0505db8fa860
function col(t,land)
	hd=header(t)
	for i=1:length(header(t))
		if land==hd[i]
			return i
		end
	end
	nothing
end

# ╔═╡ 32269a70-220b-11eb-0c7b-b389a74d0654
header(t)

# ╔═╡ f784c84c-220a-11eb-18c1-9f9ccbc76377
nrows(t)=length(t[2].children)-2

# ╔═╡ e2e5c45e-220a-11eb-1c25-fb79acd8d277
function row(t,i)
	tr=t[2][i+1]
	start=2
	if haskey(tr[1].attributes,"rowspan")
		start=3
	end
	[ replace(tr[i][1].text, r"\n|\t|\." => s"")  for i=start:length(tr.children)-2]
end

# ╔═╡ 9a7b1f6c-2213-11eb-2eed-0bfd1e360bc2
row(t,50)

# ╔═╡ 02467d42-2211-11eb-2b75-79953c5ca690
date_start=Date(2020,2,23)

# ╔═╡ bd80acd6-220f-11eb-223a-134dd4b501aa
num(s)= try parse(Int64,s) catch e 0 end

# ╔═╡ 406c112e-220f-11eb-22ff-e95bf1ed6128
begin
df=DataFrame(Datum=[date_start+Day(i) for i=1:nrows(t)])
for icol=1:length(header(t))
	df[Symbol(header(t)[icol])]=[num(row(t,i)[icol]) for i=1:nrows(t)]
end
end

# ╔═╡ 5455d774-220f-11eb-1038-8db3d3646f6f
df

# ╔═╡ 43270c24-2210-11eb-0789-e733ddad3fda
begin
	fig=figure(2)
	clf()
    xhd=header(t)
	for i=1:length(xhd)
	plot_date(df[:Datum],df[Symbol(xhd[i])],"-",label=xhd[i])
	end
	legend(loc="upper left")
	fig
end

# ╔═╡ edff59b8-2214-11eb-32a1-bd37af96c44b
CSV.write("rki_scraped.csv",df)

# ╔═╡ 90756dc8-2205-11eb-092b-e1f28bf76325
parsed_page.root[2][3][5][7][1][17];

# ╔═╡ Cell order:
# ╠═5d84ea34-2202-11eb-0a5e-8153f59b5785
# ╠═d02e25a2-2202-11eb-21b7-89a516182e0f
# ╠═fb46d6f8-2202-11eb-1cda-8509c49cb0de
# ╠═ce05b5c6-2205-11eb-3533-b3a8e0450f2e
# ╠═aaf63af0-2206-11eb-0c77-df3ffa48b31b
# ╠═226610d0-220a-11eb-1a51-e913fa1a9120
# ╠═f85174da-220c-11eb-3bdc-0505db8fa860
# ╠═05f3b3a0-2208-11eb-07f4-cbfd72d3aa4a
# ╠═32269a70-220b-11eb-0c7b-b389a74d0654
# ╠═f784c84c-220a-11eb-18c1-9f9ccbc76377
# ╠═e2e5c45e-220a-11eb-1c25-fb79acd8d277
# ╠═9a7b1f6c-2213-11eb-2eed-0bfd1e360bc2
# ╠═02467d42-2211-11eb-2b75-79953c5ca690
# ╠═bd80acd6-220f-11eb-223a-134dd4b501aa
# ╠═406c112e-220f-11eb-22ff-e95bf1ed6128
# ╠═5455d774-220f-11eb-1038-8db3d3646f6f
# ╠═43270c24-2210-11eb-0789-e733ddad3fda
# ╠═edff59b8-2214-11eb-32a1-bd37af96c44b
# ╠═90756dc8-2205-11eb-092b-e1f28bf76325

max_stars_repo_name: j-fu/coronaplot

max_issues_repo_name: j-fu/coronaplot

max_forks_repo_name: j-fu/coronaplot


text:
using Distributions, UncertainData, Combinatorics

# Test all combinations of different types of uncertain values
M = MixtureModel([Normal(3, 0.2), Normal(2, 1)])

r1 = UncertainValue(Normal, rand(), rand())
r2 = UncertainValue(rand(M, 10000))
r3 = UncertainValue(Normal, rand(Normal(4, 3.2), 10000))

uvals = [r1; r2; r3]

a = UncertainIndexDataset(uvals)
b = UncertainIndexDataset(uvals[end:-1:1])
N = length(a)

###########
# Addition
###########
@test a + b isa UncertainIndexDataset
@test 2 + b isa UncertainIndexDataset
@test [2 for i = 1:N] + b isa UncertainIndexDataset
@test 2.2 + b isa UncertainIndexDataset
@test [2.2 for i = 1:N] + b isa UncertainIndexDataset
@test a + [2 for i = 1:N] isa UncertainIndexDataset
@test a + [2.2 for i = 1:N] isa UncertainIndexDataset

#############
# Subtraction
#############
@test a - b isa UncertainIndexDataset
@test 2 - b isa UncertainIndexDataset
@test [2 for i = 1:N] - b isa UncertainIndexDataset
@test 2.2 - b isa UncertainIndexDataset
@test [2.2 for i = 1:N] - b isa UncertainIndexDataset
@test a - [2 for i = 1:N] isa UncertainIndexDataset
@test a - [2.2 for i = 1:N] isa UncertainIndexDataset

################
# Multiplication
################
@test a * b isa UncertainIndexDataset
@test 2 * b isa UncertainIndexDataset
@test [2 for i = 1:N] * b isa UncertainIndexDataset
@test 2.2 * b isa UncertainIndexDataset
@test [2.2 for i = 1:N] * b isa UncertainIndexDataset
@test a * [2 for i = 1:N] isa UncertainIndexDataset
@test a * [2.2 for i = 1:N] isa UncertainIndexDataset

################
# Division
################
@test a * b isa UncertainIndexDataset
@test 2 * b isa UncertainIndexDataset
@test [2 for i = 1:N] * b isa UncertainIndexDataset
@test 2.2 * b isa UncertainIndexDataset
@test [2.2 for i = 1:N] * b isa UncertainIndexDataset
@test a * [2 for i = 1:N] isa UncertainIndexDataset
@test a * [2.2 for i = 1:N] isa UncertainIndexDataset

################
# Exponentiation
################
@test a * b isa UncertainIndexDataset
@test 2 * b isa UncertainIndexDataset
@test [2 for i = 1:N] * b isa UncertainIndexDataset
@test 2.2 * b isa UncertainIndexDataset
@test [2.2 for i = 1:N] * b isa UncertainIndexDataset
@test a * [2 for i = 1:N] isa UncertainIndexDataset
@test a * [2.2 for i = 1:N] isa UncertainIndexDataset


text:
export Polynomial, exponents, coefficients, variables, nvariables, nterms,
    ishomogenized, ishomogenous, homogenize, dehomogenize, degree, evaluate, substitute,
    differentiate, ∇, weylnorm, weyldot, scale_coefficients!

"""
    Polynomial(p::MultivariatePolynomials.AbstractPolynomial [, variables [, homogenized=false]])

A structure for fast evaluation of multivariate polynomials.
The terms are sorted first by total degree, then lexicographically.
`Polynomial` has first class support for [homogenous polynomials](https://en.wikipedia.org/wiki/Homogeneous_polynomial).
This field indicates whether the first variable should be considered as the homogenization variable.

    Polynomial{T}(p::MultivariatePolynomials.AbstractPolynomial [, variables [, homogenized=false]])

You can force a coefficient type `T`. For optimal performance `T` should be same type as
the input to with which it will be evaluated.

    Polynomial(exponents::Matrix{Int}, coefficients::Vector{T}, variables, [, homogenized=false])

You can also create a polynomial directly. Note that in exponents each column represents the exponent of a term.

### Example
```
Poly([3 1; 1 1; 0 2 ], [-2.0, 3.0], [:x, :y, :z]) == 3.0x^2yz^2 - 2x^3y
```
"""
mutable struct Polynomial{T<:Number}
    exponents::Matrix{Int}
    coefficients::Vector{T}
    variables::Vector{Symbol}
    homogenized::Bool

    function Polynomial{T}(exponents::Matrix{Int}, coefficients::Vector{T}, variables::Vector{Symbol}, homogenized::Bool) where {T<:Number}
        sorted_cols = sort!([1:size(exponents,2);], lt=((i, j) -> lt_total_degree(exponents[:,i], exponents[:,j])), rev=true)
        exps = exponents[:, sorted_cols]
        coeffs = coefficients[sorted_cols]
        new(exps, coeffs, copy(variables), homogenized)
    end
end

function Polynomial(exponents::Matrix{Int}, coefficients::Vector{T}, variables::Vector{Symbol}, homogenized::Bool=false) where {T<:Number}
    Polynomial{T}(exponents, coefficients, variables, homogenized)
end

function Polynomial(exponents::Matrix{Int}, coefficients::Vector{T}, homogenized::Bool=false) where {T<:Number}
    Polynomial{T}(exponents, coefficients, [Symbol("x$(i)") for i=1:size(exponents, 1)], homogenized)
end
function Polynomial(p::MP.AbstractPolynomialLike, variables::Vector{<:MP.AbstractVariable}, homogenized::Bool=false)
    exps, coefficients = _coefficients_exponents(p, variables)
    Polynomial(exps, coefficients, Symbol.(variables), homogenized)
end
Polynomial(p::MP.AbstractPolynomialLike, homogenized::Bool=false) = Polynomial(p, MP.variables(p), homogenized)

function Polynomial{T}(p::MP.AbstractPolynomialLike, variables::Vector{<:MP.AbstractVariable}, homogenized::Bool=false) where {T<:Number}
    exps, coefficients = _coefficients_exponents(p, variables)
    Polynomial{T}(exps, convert(Vector{T}, coefficients), Symbol.(variables), homogenized)
end
Polynomial{T}(p::MP.AbstractPolynomialLike, homogenized::Bool=false) where {T<:Number} = Polynomial{T}(p, MP.variables(p), homogenized)

function _coefficients_exponents(poly::MP.AbstractPolynomialLike{T}, vars) where {T}
    terms = MP.terms(poly)
    nterms = length(terms)
    nvars = length(vars)
    exps = Matrix{Int}(undef, nvars, nterms)
    coefficients = Vector{T}(undef, nterms)
    for j = 1:nterms
        term = terms[j]
        coefficients[j] = MP.coefficient(term)
        for i = 1:nvars
            exps[i,j] = MP.degree(term, vars[i])
        end
    end
    exps, coefficients
end

"Sorts two vectory by total degree"
function lt_total_degree(a::Vector{T}, b::Vector{T}) where {T<:Real}
    sum_a = sum(a)
    sum_b = sum(b)
    if sum_a < sum_b
        return true
    elseif sum_a > sum_b
        return false
    else
        for i in eachindex(a)
            if a[i] < b[i]
                return true
            elseif a[i] > b[i]
                return false
            end
        end
    end
    false
end

"""
    exponents(p::Polynomial)

Returns the exponents matrix of `p`. Each column represents the exponents of a term of `p`.
"""
exponents(p::Polynomial) = p.exponents

"""
    coefficients(p::Polynomial)

Returns the coefficient vector of `p`.
"""
coefficients(p::Polynomial) = p.coefficients

"""
    variables(p::Polynomial)

Returns the variables of `p`.
"""
variables(p::Polynomial) = p.variables

"""
    ishomogenized(p::Polynomial)

Checks whether `p` was homogenized.
"""
ishomogenized(p::Polynomial) = p.homogenized

"""
    nterms(p::Polynomial)

Returns the number of terms of p
"""
nterms(p::Polynomial) = size(exponents(p), 2)

"""
    nvariables(p::Polynomial)

Returns the number of variables of p
"""
nvariables(p::Polynomial) = size(exponents(p), 1)

"""
    degree(p::Polynomial)

Returns the (total) degree of p.
"""
degree(p::Polynomial) = sum(exponents(p)[:,1])

==(p::Polynomial, q::Polynomial) = exponents(p) == exponents(q) && coefficients(p) == coefficients(q)
Base.isequal(p::Polynomial, q::Polynomial) = exponents(p) == exponents(q) && coefficients(p) == coefficients(q)

# ITERATOR
function Base.iterate(p::Polynomial, state...)
    n = nterms(p)
    istate = iterate(1:n, state...)
    istate === nothing && return nothing

    i, state = istate
    (coefficients(p)[i], exponents(p)[:,i]), state
end
Base.length(p::Polynomial) = nterms(p)
Base.eltype(p::Polynomial{T}) where {T} = T

Base.broadcastable(p::Polynomial) = Ref(p)

@inline pow(x::AbstractFloat, k::Integer) = Base.FastMath.pow_fast(x, k)
#@inline pow(x::Complex, k::Integer) = k == 1 ? x : x^k
# simplified from Base.power_by_squaring
@inline function pow(x::Number, p::Integer)
    if p == 1
        return copy(x)
    elseif p == 0
        return one(x)
    elseif p == 2
        return x*x
    end
    t = trailing_zeros(p) + 1
    p >>= t
    while (t -= 1) > 0
        x *= x
    end
    y = x
    while p > 0
        t = trailing_zeros(p) + 1
        p >>= t
        while (t -= 1) >= 0
            x *= x
        end
        y *= x
    end
    return y
end

"""
    evaluate(p::Polynomial{T}, x::AbstractVector{T})

Evaluates `p` at `x`, i.e. ``p(x)``.
`Polynomial` is also callable, i.e. you can also evaluate it via `p(x)`.
"""
function evaluate(p::Polynomial{S}, x::AbstractVector{T}) where {S<:Number, T<:Number}
    cfs = coefficients(p)
    exps = exponents(p)
    nvars, nterms = size(exps)
    res = zero(promote_type(S,T))
    for j = 1:nterms
        term = p.coefficients[j]
        for i = 1:nvars
            k = exps[i, j]
            term *= pow(x[i], k)
        end
        res += term
    end
    res
end

(p::Polynomial)(x) = evaluate(p, x)

evaluate(F::Vector{<:Polynomial}, x::AbstractVector) = map(f -> evaluate(f, x), F)
"""
    substitute(p::Polynomial, i, x)

In `p` substitute for the variable with index `i` the value `x`. You can use this for partial
evaluation of polynomial.

### Example
```julia-repl
julia> substitute(x^2+3y, 2, 5)
x^2+15
```
"""
function substitute(p::Polynomial{S}, varindex, x::T) where {S<:Number, T<:Number}
    cfs = coefficients(p)
    exps = exponents(p)
    nvars, nterms = size(exps)

    new_coefficients = Vector{promote_type(S,T)}()
    new_exps = Vector{Vector{Int}}()

    for j = 1:nterms
        coeff = cfs[j]
        exp = Vector{Int}(undef, nvars - 1)
        # first we calculate the new coefficient and remove the varindex-th row
        for i = 1:nvars
            if i == varindex
                coeff *= x^(exps[i, j])
            elseif i > varindex
                exp[i-1] = exps[i, j]
            else
                exp[i] = exps[i, j]
            end
        end
        # now we have to delete possible duplicates
        found_duplicate = false
        for k = 1:length(new_exps)
            if new_exps[k] == exp
                new_coefficients[k] += coeff
                found_duplicate = true
                break
            end
        end
        if !found_duplicate
            push!(new_coefficients, coeff)
            push!(new_exps, exp)
        end
    end

    # now we have to create a new matrix and return the poly
    Polynomial(hcat(new_exps...), new_coefficients, [variables(p)[1:varindex-1]; variables(p)[varindex+1:end]], p.homogenized)
end

"""
    differentiate(p::Polynomial, varindex::Int)

Differentiate `p` w.r.t the `varindex`th variable.

    differentiate(p::Polynomial)

Differentiate `p` w.r.t. all variables.
"""
function differentiate(p::Polynomial{T}, i_var) where T
    exps = copy(exponents(p))
    cfs = copy(coefficients(p))
    n_vars, n_terms = size(exps)

    zerocolumns = Int[]
    for j=1:n_terms
        k = exps[i_var, j]
        if k > 0
            exps[i_var, j] = max(0, k - 1)
            cfs[j] *= k
        else
            push!(zerocolumns , j)
        end
    end

    # now we have to get rid of all zeros
    nzeros = length(zerocolumns)
    if nzeros == 0
        return Polynomial(exps, cfs, variables(p), p.homogenized)
    end

    skipped_cols = 0
    new_exps = zeros(Int, n_vars, n_terms - nzeros)
    new_coefficients = zeros(T, n_terms - nzeros)

    for j = 1:n_terms
        # if we not yet have skipped all zero columns
        if skipped_cols < nzeros && j == zerocolumns[skipped_cols+1]
            skipped_cols += 1
            continue
        end

        new_exps[:, j - skipped_cols] = exps[:, j]
        new_coefficients[j - skipped_cols] = cfs[j]
    end

    Polynomial(new_exps, new_coefficients, variables(p), p.homogenized)
end
differentiate(poly::Polynomial) = map(i -> differentiate(poly, i), 1:nvariables(poly))

"""
    ∇(p::Polynomial)

Returns the gradient vector of `p`. This is the same as [`differentiate`](@ref).
"""
@inline ∇(poly::Polynomial) = differentiate(poly)

"""
    ishomogenous(p::Polynomial)

Checks whether `p` is a homogenous polynomial. Note that this is unaffected from the
value of `homogenized(p)`.
"""
function ishomogenous(p::Polynomial)
    monomials_degree = sum(exponents(p), dims=1)
    max_deg = monomials_degree[1]
    all(x -> x == max_deg, monomials_degree)
end

"""
    homogenize(p::Polynomial [, variable = :x0])

Makes `p` homogenous, if `ishomogenized(p)` is `true` this is just the identity.
The homogenization variable will always be considered as the first variable of the polynomial.
"""
function homogenize(p::Polynomial, variable::Symbol=:x0; respect_homogenous=true)
    if p.homogenized || (respect_homogenous && ishomogenous(p))
        p
    else
        monomials_degree = sum(exponents(p), dims=1)
        max_deg = monomials_degree[1]
        Polynomial([max_deg .- monomials_degree; exponents(p)], coefficients(p), [variable; variables(p)], true)
    end
end

"""
    dehomogenize(p::Polynomial)

Substitute `1` as for the first variable `p`, if `ishomogenized(p)` is `false` this is just
the identity.
"""
function dehomogenize(p::Polynomial)
    if !p.homogenized
        p
    else
        Polynomial(exponents(p)[2:end,:], coefficients(p), variables(p)[2:end], false)
    end
end

"Computes the multinomial coefficient (|k| \\over k)"
function multinomial(k::Vector{Int})
    s = 0
    result = 1
    @inbounds for i in k
        s += i
        result *= binomial(s, i)
    end
    result
end

"""
    weyldot(f::Polynomial, g::Polynomial)

Compute the [Bombieri-Weyl dot product](https://en.wikipedia.org/wiki/Bombieri_norm).
Note that this is only properly defined if `f` and `g` are homogenous.

    weyldot(f::Vector{Polynomial}, g::Vector{Polynomial})

Compute the dot product for vectors of polynomials.
"""
function weyldot(f::Polynomial,g::Polynomial)
    if (f === g)
        return sum(x -> abs2(x[1]) / multinomial(x[2]), f)
    end
    result = 0
    for (c_f, exp_f) in f
        normalizer = multinomial(exp_f)
        for (c_g, exp_g) in g
            if exp_f == exp_g
                result += (c_f * conj(c_g)) / normalizer
                break
            end
        end
    end
    result
end

function weyldot(F::Vector{Polynomial{T}}, G::Vector{Polynomial{S}}) where {T,S}
    res = zero(promote_type(S, T))
    for (f, g) in zip(F, G)
        res += weyldot(f, g)
    end
    res
end

"""
    weylnorm(f::Polynomial)

Compute the [Bombieri-Weyl norm](https://en.wikipedia.org/wiki/Bombieri_norm).
Note that this is only properly defined if `f` is homogenous.
"""
weylnorm(f::Polynomial) = √weyldot(f,f)
weylnorm(f::Vector{Polynomial{T}}) where T = √weyldot(f,f)

"""
    scale_coefficients!(f::Polynomial, λ)

Scale the coefficients of `f` with the factor `λ`.
"""
function scale_coefficients!(f::Polynomial, λ)
    f.coefficients .*= λ
end


text:
using Test
using JLD

include("core_julia1p3.jl")
###############################################################################
### This test script aims to verify correctness of our code,
### by rerunning the results reported in Table 4 of our paper.
### Note that we obtained the "optimal" value of these problems by running both
### the code of Berk+Bertsimas and our outer-approximation method, and verifying
### that we obtained the same optimal value from both methods.
### We have hard-coded some of the smaller test instances here, for simplicity.
###############################################################################

################################################################################
######## Test #1: Pitprops data set
################################################################################
@testset "Table 4" begin

pitprops=[[1,0.954,0.364,0.342,-0.129,0.313,0.496,0.424,0.592,0.545,0.084,-0.019,0.134];
       [0.954,1,0.297,0.284,-0.118,0.291,0.503,0.419,0.648,0.569,0.076,-0.036,0.144];
       [0.364,0.297,1,0.882,-0.148,0.153,-0.029,-0.054,0.125,-0.081,0.162,0.22,0.126];
       [0.342,0.284,0.882,1,0.22,0.381,0.174,-0.059,0.137,-0.014,0.097,0.169,0.015];
       [-0.129,-0.118,-0.148,0.22,1,0.364,0.296,0.004,-0.039,0.037,-0.091,-0.145,-0.208];
       [0.313,0.291,0.153,0.381,0.364,1,0.813,0.09,0.211,0.274,-0.036,0.024,-0.329];
       [0.496,0.503,-0.029,0.174,0.296,0.813,1,0.372,0.465,0.679,-0.113,-0.232,-0.424];
       [0.424,0.419,-0.054,-0.059,0.004,0.09,0.372,1,0.482,0.557,0.061,-0.357,-0.202];
       [0.592,0.648,0.125,0.137,-0.039,0.211,0.465,0.482,1,0.526,0.085,-0.127,-0.076];
       [0.545,0.569,-0.081,-0.014,0.037,0.274,0.679,0.557,0.526,1,-0.319,-0.368,-0.291];
       [0.084,0.076,0.162,0.097,-0.091,-0.036,-0.113,0.061,0.085,-0.319,1,0.029,0.007];
       [-0.019,-0.036,0.22,0.169,-0.145,0.024,-0.232,-0.357,-0.127,-0.368,0.029,1,0.184];
       [0.134,0.144,0.126,0.015,-0.208,-0.329,-0.424,-0.202,-0.076,-0.291,0.007,0.184,1];]
pitprops=reshape(pitprops, (13,13));

println("Running pitprops data set k=5 with Problem (18s-relax)")

k=5
lambda_true=3.406
@time UB, LB=getSDPUpperBound_gd(pitprops, k, true, false)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

println("Running pitprops data set k=5 with Problem (18s-relax-cuts)")

@time UB, LB=getSDPUpperBound_gd_highdim(pitprops, k, true, true)
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

k=10
println("Running pitprops data set k=10 with Problem (18s-relax)")
lambda_true=4.173
@time UB, LB=getSDPUpperBound_gd(pitprops, k, true, false)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

println("Running pitprops data set k=10 with Problem (18s-relax-cuts)")

@time UB, LB=getSDPUpperBound_gd_highdim(pitprops, k, true, true)
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

# ################################################################################
# #### Clear data
# ################################################################################
pitprops=nothing;
#
# ###############################################################################
# ####### Test #2: Norm wine data set.
# ##############################################################################
#
normwine=[[1.0       ,  0.0943969 ,  0.211545   , -0.310235  ,  0.270798  ,  0.289101 ,   0.236815 , -0.155929 ,  0.136698    , 0.546364   ,-0.0717472 ,  0.0723432 ,   0.64372   ];
[0.0943969 ,  1.0       ,  0.164045   ,  0.2885    , -0.0545751 , -0.335167 ,  -0.411007 ,  0.292977 , -0.220746    , 0.248985   ,-0.561296  , -0.36871   ,  -0.192011  ];
[0.211545  ,  0.164045  ,  1.0        ,  0.443367  ,  0.286587  ,  0.12898  ,   0.115077 ,  0.18623  ,  0.00965194  , 0.258887   ,-0.0746669 ,  0.00391123,   0.223626  ];
[-0.310235 ,   0.2885   ,   0.443367  ,   1.0      ,  -0.0833331,  -0.321113,   -0.35137 ,   0.361922,  -0.197327   ,  0.018732  , -0.273955 ,  -0.276769 ,   -0.440597 ];
[0.270798  , -0.0545751 ,  0.286587   , -0.0833331 ,  1.0       ,  0.214401 ,   0.195784 , -0.256294 ,  0.236441    , 0.19995    , 0.0553982 ,  0.0660039 ,   0.393351  ];
[0.289101  , -0.335167  ,  0.12898    , -0.321113  ,  0.214401  ,  1.0      ,   0.864564 , -0.449935 ,  0.612413    ,-0.0551364  , 0.433681  ,  0.699949  ,   0.498115  ];
[0.236815  , -0.411007  ,  0.115077   , -0.35137   ,  0.195784  ,  0.864564 ,   1.0      , -0.5379   ,  0.652692    ,-0.172379   , 0.543479  ,  0.787194  ,   0.494193  ];
[-0.155929 ,   0.292977 ,   0.18623   ,   0.361922 ,  -0.256294 ,  -0.449935,   -0.5379  ,   1.0     ,  -0.365845   ,  0.139057  , -0.26264  ,  -0.50327  ,   -0.311385 ];
[0.136698  , -0.220746  ,  0.00965194 , -0.197327  ,  0.236441  ,  0.612413 ,   0.652692 , -0.365845 ,  1.0         ,-0.0252499  , 0.295544  ,  0.519067  ,   0.330417  ];
[0.546364  ,  0.248985  ,  0.258887   ,  0.018732  ,  0.19995   , -0.0551364,  -0.172379 ,  0.139057 , -0.0252499   , 1.0        ,-0.521813  , -0.428815  ,   0.3161    ];
[-0.0717472,  -0.561296 ,  -0.0746669 ,  -0.273955 ,   0.0553982,   0.433681,    0.543479,  -0.26264 ,   0.295544   , -0.521813  ,  1.0      ,   0.565468 ,    0.236183 ];
[0.0723432 , -0.36871   ,  0.00391123 , -0.276769  ,  0.0660039 ,  0.699949 ,   0.787194 , -0.50327  ,  0.519067    ,-0.428815   , 0.565468  ,  1.0       ,   0.312761  ];
[0.64372   , -0.192011  ,  0.223626   , -0.440597  ,  0.393351  ,  0.498115 ,   0.494193 , -0.311385 ,  0.330417    , 0.3161     , 0.236183  ,  0.312761  ,   1.0	];]
normwine=reshape(normwine, (13,13))

println("Running wine data set k=5 with Problem (18s-relax)")

k=5
lambda_true=3.43978
@time UB, LB=getSDPUpperBound_gd(normwine, k, true, false)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

println("Running wine data set k=5 with Problem (18s-relax-cuts)")

@time UB, LB=getSDPUpperBound_gd_highdim(normwine, k, true, true)
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

k=10
println("Running wine data set k=10 with Problem (18s-relax)")
lambda_true=4.59429
@time UB, LB=getSDPUpperBound_gd(normwine, k, true, false)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

println("Running wine data set k=5 with Problem (18s-relax-cuts)")

@time UB, LB=getSDPUpperBound_gd_highdim(normwine, k, true, true)
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

# # ################################################################################
# # #### Clear data
# # ################################################################################
normwine=nothing;

################################################################################
######## Test #3: normMiniBooNE data set.
###############################################################################
normminiboone=load("data/miniBoone.jld",  "normMiniBooNE")

println("Running normMiniBooNE data set k=5 with Problem (18s-relax)")

k=5
lambda_true=5.0000
@time UB, LB=getSDPUpperBound_gd(normminiboone, k, true, false)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

println("Running normMiniBooNE data set k=5 with Problem (18s-relax-cuts)")

@time UB, LB=getSDPUpperBound_gd_highdim(normminiboone, k, true, true)
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

k=10
println("Running normMiniBooNE data set k=10 with Problem (18s)")
lambda_true=9.9999999
@time UB, LB=getSDPUpperBound_gd(normminiboone, k, true, false)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

println("Running normMiniBooNE data set k=10 with Problem (18s-relax-cuts)")

@time UB, LB=getSDPUpperBound_gd_highdim(normminiboone, k, true, true)
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

# # ################################################################################
# # #### Clear data
# # ################################################################################
normminiboone=nothing;

################################################################################
######## Test #4: normCommunities data set.
###############################################################################
normcommunities=load("data/communities.jld", "normCommunities")

println("Running normCommunites data set k=5 with Problem (18s-relax)")

k=5
lambda_true=4.86051
@time UB, LB=getSDPUpperBound_gd(normcommunities, k, true, false)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

println("Running normCommunites data set k=5 with Problem (18s-relax-cuts)")

@time UB, LB=getSDPUpperBound_gd_highdim(normcommunities, k, true, true)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

k=10
println("Running normCommunites data set k=10 with Problem (18s)")
lambda_true=8.8236
@time UB, LB=getSDPUpperBound_gd(normcommunities, k, true, false)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

println("Running normCommunites data set k=10 with Problem (18s-relax-cuts)")

@time UB, LB=getSDPUpperBound_gd_highdim(normcommunities, k, true, true)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

# # ################################################################################
# # #### Clear data
# # ################################################################################
normcommunities=nothing;

################################################################################
######## Test #4: normArrhythmia data set.
###############################################################################
normarrhythmia=load("data/arrhythmia.jld", "normArrhythmia");
@show size(normarrhythmia,1)
for i=1:size(normarrhythmia,1)
      for j=1:size(normarrhythmia,2)
            if isnan(normarrhythmia[i,j])
                  normarrhythmia[i,j]=0.0
            end
      end
end

println("Running normArrhythmia data set k=5 with Problem (18s-relax)")

k=5
lambda_true=4.2321
@time UB, LB=getSDPUpperBound_gd(normarrhythmia, k, true, false)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

println("Running normArrhythmia data set k=5 with Problem (18s-relax-cuts)")

@time UB, LB=getSDPUpperBound_gd_highdim(normarrhythmia, k, true, true)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

k=10
println("Running normArrhythmia data set k=10 with Problem (18s)")
lambda_true=7.53938
@time UB, LB=getSDPUpperBound_gd(normarrhythmia, k, true, false)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

println("Running normArrhythmia data set k=10 with Problem (18s-relax-cuts)")

@time UB, LB=getSDPUpperBound_gd_highdim(normarrhythmia, k, true, true)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

# # ################################################################################
# # #### Clear data
# # ################################################################################
normarrhythmia=nothing;

################################################################################
######## Test #4: normMicroMass data set.
###############################################################################
micromass=load("data/microMass.jld", "normMicroMass");

println("Running normMicroMass data set k=5 with Problem (18s-relax)")

k=5
lambda_true=4.936405
@time UB, LB=getSDPUpperBound_gd(micromass, k, true, false)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

println("Running normMicroMass data set k=5 with Problem (18s-relax-cuts)")

@time UB, LB=getSDPUpperBound_gd_highdim(micromass, k, true, true)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

k=10
println("Running normMicroMass data set k=10 with Problem (18s)")
lambda_true=9.045361166201573
@time UB, LB=getSDPUpperBound_gd(micromass, k, true, false)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

println("Running normMicroMass data set k=10 with Problem (18s-relax-cuts)")

@time UB, LB=getSDPUpperBound_gd_highdim(micromass, k, true, true)
# @show UB, LB
@show "R-gap is:" (UB-lambda_true)/lambda_true
@show "O-gap is:" (lambda_true-LB)/LB

# # ################################################################################
# # #### Clear data
# # ################################################################################
micromass=nothing;
end


text:
struct CompositeControl{A,B}
    a::A
    b::B
    function CompositeControl(a::A, b::B) where {A, B}
        # `_in` method extended post-facto
        a isa StoppingCriterion && _in(a, b) && return b
        b isa StoppingCriterion && _in(b, a) && return a
        return new{A,B}(a, b)
    end
end
CompositeControl() = Never()
CompositeControl(a) = a
CompositeControl(a, b, c...) = CompositeControl(CompositeControl(a,b), c...)

"""
    composite(controls...)

Construct an iteration control that applies the specified `controls`
in sequence.

"""
composite(controls...) = CompositeControl(controls...)

update!(c::CompositeControl, m, v, n) =
    (a = update!(c.a, m, v, n), b = update!(c.b, m, v, n))
update!(c::CompositeControl, m, v, n, state) =
    (a = update!(c.a, m, v, n, state.a),
     b = update!(c.b, m, v, n, state.b))

# # RECURSION TO FLATTEN A CONTROL OR ITS STATE

flat(state) = (state,)
flat(state::NamedTuple{(:a,:b)}) = tuple(flat(state.a)..., flat(state.b)...)
flat(d::CompositeControl) = tuple(flat(d.a)..., flat(d.b)...)

_in(::Never, ::Any) = true
_in(::Never, ::CompositeControl) = true
_in(c1::Any, c2::Any) = c1 == c2
_in(c::Any, d::CompositeControl) = c in flat(d)
_in(::CompositeControl, ::Any) = false

# # DISPLAY

function Base.show(io::IO, c::CompositeControl)
    list = join(string.(flat(c)), ", ")
    print(io, "CompositeControl($list)")
end

# # RECURSION TO DEFINE `done`

# fallback for atomic controls:
_done(control, state, old_done) = old_done || done(control, state)
%
# composite:
_done(c::CompositeControl, state, old_done) =
    _done(c.a, state.a, _done(c.b, state.b, old_done))

done(c::CompositeControl, state) = _done(c, state, false)

# # RECURSION TO DEFINE `takedown`

# fallback for atomic controls:
function _takedown(control, v, state, old_takedown)
    new_takedown = (control, takedown(control, v, state))
    return tuple(old_takedown..., new_takedown)
end

# composite:
_takedown(c::CompositeControl, v, state, old_takedown) =
    _takedown(c.b, v, state.b,
              _takedown(c.a, v, state.a, old_takedown))

takedown(c::CompositeControl, v, state) = _takedown(c, v, state, ())

# # TRAITS

for ex in [:needs_loss, :needs_training_losses]
    quote
        $ex(c::CompositeControl) = any($ex, flat(c))
    end |> eval
end


text:
using Base.Test
import ArchGDAL; const AG = ArchGDAL

@testset "Tests for field defn" begin
    AG.createfielddefn("fieldname",AG.OFTInteger) do fd
        println(fd)
        println("name: $(AG.getname(fd))")
        println("Setting name to \"newname\"")
        AG.setname!(fd, "newname")
        println(AG.getname(fd))
        println("type: $(AG.gettype(fd))")
        println("Setting name to $(AG.OFTDate)")
        AG.settype!(fd, AG.OFTDate)
        println(AG.gettype(fd))
        AG.settype!(fd, AG.OFTInteger)
        @test AG.getsubtype(fd) == AG.OFSTNone
        AG.setsubtype!(fd, AG.OFSTInt16)
        @test AG.getsubtype(fd) == AG.OFSTInt16
        AG.setsubtype!(fd, AG.OFSTBoolean)
        @test AG.getsubtype(fd) == AG.OFSTBoolean
        AG.setsubtype!(fd, AG.OFSTNone)
        @test AG.getjustify(fd) == AG.OJUndefined
        AG.setjustify!(fd, AG.OJLeft)
        @test AG.getjustify(fd) == AG.OJLeft
        println("width: $(AG.getwidth(fd))")
        AG.setwidth!(fd, 10)
        println("after setting width to 10: new width=$(AG.getwidth(fd))")
        println("precision: $(AG.getprecision(fd))")
        AG.setprecision!(fd, 20)
        println("after setting to 20: new precision=$(AG.getprecision(fd))")
        AG.setparams!(fd, "finalname", AG.OFTDate, nwidth=5, nprecision=2,
                      justify=AG.OJRight)
        println("type: $((AG.gettype(fd),AG.getname(fd),AG.getsubtype(fd),
                          AG.getjustify(fd),AG.getwidth(fd),AG.getprecision(fd)))")
        @test AG.isignored(fd) == false
        AG.setignored!(fd, true)
        @test AG.isignored(fd) == true
        AG.setignored!(fd, false)
        @test AG.isignored(fd) == false

        @test AG.isnullable(fd) == true
        AG.setnullable!(fd, false)
        @test AG.isnullable(fd) == false
        AG.setnullable!(fd, true)
        @test AG.isnullable(fd) == true

        @test AG.getdefault(fd) == ""
        AG.setdefault!(fd, "0001/01/01 00:00:00")
        @test AG.getdefault(fd) == "0001/01/01 00:00:00"
        @test AG.isdefaultdriverspecific(fd) == true
    end
end

@testset "Tests for Geom Field Defn" begin
    AG.creategeomfielddefn("geomname", AG.wkbPolygon) do gfd
        @test AG.getname(gfd) == "geomname"
        AG.setname!(gfd, "my name!")
        @test AG.getname(gfd) == "my name!"

        @test AG.gettype(gfd) == AG.wkbPolygon
        AG.settype!(gfd, AG.wkbPolyhedralSurface)
        @test AG.gettype(gfd) == AG.wkbPolyhedralSurface

        println(AG.getspatialref(gfd))
        # AG.setspatialref!(gfd, AG.unsafe_fromEPSG(4326))
        # println(AG.getspatialref(gfd))

        @test AG.isignored(gfd) == false
        AG.setignored!(gfd, true)
        @test AG.isignored(gfd) == true
        AG.setignored!(gfd, false)
        @test AG.isignored(gfd) == false

        @test AG.isnullable(gfd) == true
        AG.setnullable!(gfd, false)
        @test AG.isnullable(gfd) == false
        AG.setnullable!(gfd, true)
        @test AG.isnullable(gfd) == true
    end
end

@testset "Tests for Feature Defn" begin
    AG.createfeaturedefn("new_feature") do fd
        println(AG.nreference(fd))
        AG.reference(fd)
        println(AG.nreference(fd))
        AG.reference(fd)
        println(AG.nreference(fd))
        AG.reference(fd)
        println(AG.nreference(fd))
        AG.release(fd)
        println(AG.nreference(fd))
        AG.dereference(fd)
        println(AG.nreference(fd))
        AG.dereference(fd)
        AG.createfielddefn("fieldname",AG.OFTInteger) do fielddef
            @test AG.nfield(fd) == 0
            AG.addfielddefn!(fd, fielddef)
            @test AG.nfield(fd) == 1
            AG.addfielddefn!(fd, fielddef)
            @test AG.nfield(fd) == 2
            AG.addfielddefn!(fd, fielddef)
            @test AG.nfield(fd) == 3
            AG.createfielddefn("newfield",AG.OFTInteger) do fielddef2
                AG.addfielddefn!(fd, fielddef2)
                @test AG.nfield(fd) == 4
                for i in 0:3
                    println("$i : $(AG.getname(AG.getfielddefn(fd,i)))")
                end
            end
            AG.deletefielddefn!(fd, 0)
            @test AG.nfield(fd) == 3
            for i in 0:2
                println("$i : $(AG.getname(AG.getfielddefn(fd,i)))")
            end
            println("After reordering (in reverse):")
            AG.reorderfielddefns!(fd, Cint[2,1,0])
            @test AG.nfield(fd) == 3
            for i in 0:2
                println("$i : $(AG.getname(AG.getfielddefn(fd,i)))")
            end
        end
        @test AG.ngeomfield(fd) == 1
        @test AG.getgeomtype(fd) == AG.wkbUnknown
        AG.setgeomtype!(fd, AG.wkbPolygon)
        @test AG.getgeomtype(fd) == AG.wkbPolygon
        @test AG.ngeomfield(fd) == 1

        @test AG.isgeomignored(fd) == false
        AG.setgeomignored!(fd, true)
        @test AG.isgeomignored(fd) == true
        AG.setgeomignored!(fd, false)
        @test AG.isgeomignored(fd) == false

        @test AG.isstyleignored(fd) == false
        AG.setstyleignored!(fd, true)
        @test AG.isstyleignored(fd) == true
        AG.setstyleignored!(fd, false)
        @test AG.isstyleignored(fd) == false

        @test AG.getgeomfieldindex(fd) == 0
        gfd0 = AG.getgeomfielddefn(fd, 0)
        @test AG.ngeomfield(fd) == 1
        AG.addgeomfielddefn!(fd, gfd0)
        @test AG.ngeomfield(fd) == 2
        gfd1 = AG.getgeomfielddefn(fd, 1)
        AG.setname!(gfd0, "name0")
        AG.setname!(gfd1, "name1")
        @test AG.getgeomfieldindex(fd, "") == -1
        @test AG.getgeomfieldindex(fd, "name0") == 0
        @test AG.getgeomfieldindex(fd, "name1") == 1
        AG.deletegeomfielddefn!(fd, 0)
        @test AG.ngeomfield(fd) == 1
        @test AG.getgeomfieldindex(fd, "") == -1
        @test AG.getgeomfieldindex(fd, "name0") == -1
        @test AG.getgeomfieldindex(fd, "name1") == 0
        @test AG.nreference(fd) == 0
        AG.createfeature(fd) do f
            @test AG.nreference(fd) == 2 # artificially inflated
            @test AG.issame(AG.getfeaturedefn(f), fd) == true
        end
        @test AG.nreference(fd) == 0
    end
end