text:
const DYNAMIC_DSL_TRACE = Symbol("@trace")

function make_dynamic_gen_function(name, args, body, return_type, annotations)
    escaped_args = map((arg) -> esc(arg.name), args)
    gf_args = [esc(state), escaped_args...]

    julia_fn_name = gensym(name)
    julia_fn_defn = Expr(:function,
        Expr(:call, esc(julia_fn_name), gf_args...),
        esc(body))
    arg_types = map((arg) -> esc(arg.typ), args)
    has_argument_grads = map(
        (arg) -> (DSL_ARG_GRAD_ANNOTATION in arg.annotations), args)
    accepts_output_grad = DSL_RET_GRAD_ANNOTATION in annotations

    quote
        # first define the underlying Julia function
        $julia_fn_defn

        # now wrap it in a DynamicDSLFunction value
        Core.@__doc__ $(esc(name)) = DynamicDSLFunction(
            Type[$(arg_types...)],
            $(esc(julia_fn_name)),
            $has_argument_grads,
            $(esc(return_type)),
            $accepts_output_grad)
    end
end

max_stars_repo_name: dendisuhubdy/Gen

max_issues_repo_name: kumsh/Gen

max_forks_repo_name: kumsh/Gen


text:
module EbitCloudSpatialExtend

using Revise
using QuadGK
using Interpolations
using Memoize

function phi(r, r_e, V_0)
    if r < r_e
        return V_0*(r/r_e)^2
    else
        return V_0*(2*log(r/r_e) + 1)
    end
end

function eta(qe_over_kT, r_stop, r_e, V_0, rtol=1e-6)
    f(r) = exp(-qe_over_kT*phi(r, r_e, V_0))*r
    quadgk(f, 0, r_stop, rtol=rtol, order=15)[1]
end

function ion_ion_overlap(qe_over_kT_i, qe_over_kT_j, r_stop, r_e, V_0, rtol=1e-6)
    retval = eta(qe_over_kT_i,  r_stop, r_e, V_0, rtol) / eta(qe_over_kT_j, r_stop, r_e, V_0, rtol)
    if isnan(retval)
        return 1.0
    else
        return min(1.0, retval)
    end
end

function electron_ion_overlap(qe_over_kT, r_stop, r_e, V_0, rtol=1e-6)
    retval = eta(qe_over_kT, r_e, r_e, V_0, rtol) / eta(qe_over_kT, r_stop, r_e, V_0, rtol)
    if isnan(retval)
        return 1.0
    else
        return min(1.0, retval)
    end
end

function heat_capacitance(qe_over_kT, r_stop, r_e, V_0; rtol=1e-6, nanval=2.5)
    norm = eta(qe_over_kT, r_stop, r_e, V_0)
    p2 = quadgk(r -> phi(r,r_e,V_0)^2*exp(-qe_over_kT*phi(r, r_e, V_0))*r, 0, r_stop, rtol=rtol, order=15)[1]
    p = quadgk(r -> phi(r,r_e,V_0)*exp(-qe_over_kT*phi(r, r_e, V_0))*r, 0, r_stop, rtol=rtol, order=15)[1]

    val = qe_over_kT^2 * (p2/norm - (p/norm)^2) + 1.5

    if isnan(val)
        return nanval
    else
        return val
    end
end

function effective_radius(qe_over_kT, r_stop, r_e, V_0; rtol=1e-6)
    f(r) = exp(-qe_over_kT*phi(r, r_e, V_0))*r
    g(r) = exp(-qe_over_kT*phi(r, r_e, V_0))*r*r

    return quadgk(g, 0, r_stop, rtol=rtol, order=21)[1] / quadgk(f, 0, r_stop, rtol=rtol, order=21)[1]
end

function simple_effective_radius(qe_over_kT, r_e, V_0)
    if ( V_0 * qe_over_kT > 1)
        return r_e * sqrt(1/(V_0 * qe_over_kT))
    else
        return r_e * exp(0.5*(1/(V_0*qe_over_kT) - 1))
    end
end

# comparison function
function create_interpolation(f, no_of_log_divisions=100,
                              q_over_kT_start=1e-12, q_over_kT_stop=1e12)

    _start = no_of_log_divisions*log10(q_over_kT_start)
    _end = no_of_log_divisions*log10(q_over_kT_stop)

    xs = vcat([0], [exp10(x/no_of_log_divisions) for x in _start:_end], [Inf])
    ys = [f(x) for x in xs]

    return LinearInterpolation(xs, ys)
end

function create_2d_interpolation(f, no_of_log_divisions=10,
                                 q_over_kT_start=1e-12, q_over_kT_stop=1e12)

    _start = no_of_log_divisions*log10(q_over_kT_start)
    _end = no_of_log_divisions*log10(q_over_kT_stop)

    xs = vcat([0], [exp10(x/no_of_log_divisions) for x in _start:_end], [Inf])
    zs = [f(x, y) for x in xs, y in xs]

    return LinearInterpolation((xs, xs), zs)
end

end

ext: jl


text:
##############################################################################
#
# ConvertionChecks.jl
#
# Part of CVortex.jl
# Checks that a type can be successfully converted to the types required
# for the underlying cvortex library, giving useful error messages.
#
# Copyright 2019 HJA Bird
#
# 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.
##############################################################################

"""
    Gives an assert message if something cannot be converted to something
    of type or eltype Float32.

Part of CVortex. Not for client use.
"""
function convertable_to_F32(arg :: Real, argname :: String)
    @assert(hasmethod(Float32, Tuple{typeof(arg)}),
        "Cannot convert "*argname*" of type "*
        string(typeof(arg))*" to a Float32.")
    return
end

function convertable_to_F32(arg :: Vector{<:Real}, argname :: String)
    @assert(hasmethod(Float32, Tuple{eltype(arg)}),
        "Cannot convert elements  of "*argname*" of type "*
        string(typeof(arg))*" to a Float32.")
    return
end

function convertable_to_F32(arg :: Matrix{<:Real}, argname :: String)
    @assert(hasmethod(Float32, Tuple{eltype(arg)}),
        "Cannot convert elements  of "*argname*" of type "*
		string(typeof(arg))*" to a Float32.")
    return
end

"""
    Gives an assert message if something cannot be converted to a
    CVortex.Vec3f or Vector{CVortex.Vec3f}.

Part of CVortex. Not for client use.
"""
function convertable_to_Vec3f_vect(arg :: Matrix{<:Real}, argname :: String)
    @assert(size(arg)[2]==3, "The argument "*argname*" should have size "*
        "N by 3, but actually has size "*string(size(arg))*".")
    convertable_to_F32(arg, argname)
    return
end

function convertable_to_Vec3f_vect(arg :: Vector{<:Real}, argname :: String)
    @assert(length(arg)==3, "The argument "*argname*" should have size "*
        "3, but actually has size "*string(length(arg))*".")
    convertable_to_F32(arg, argname)
    return
end

"""
    Gives an assert message if something cannot be converted to a
    CVortex.Vec2f or Vector{CVortex.Vec2f}.

Part of CVortex. Not for client use.
"""
function convertable_to_Vec2f_vect(arg :: Matrix{<:Real}, argname :: String)
    @assert(size(arg)[2]==2, "The argument "*argname*" should have size "*
        "N by 2, but actually has size "*string(size(arg))*".")
    convertable_to_F32(arg, argname)
    return
end

function convertable_to_Vec2f_vect(arg :: Vector{<:Real}, argname :: String)
    @assert(length(arg)==2, "The argument "*argname*" should have size "*
        "2, but actually has size "*string(length(arg))*".")
    convertable_to_F32(arg, argname)
    return
end

"""
    Gives an assert message if something cannot be converted to a
    CVortex.VortexFilament or Vector{CVortex.VortexFilament}.

Part of CVortex. Not for client use.
"""
function check_filament_definition(
    filament_start_coord :: Vector{<:Real},
    filament_end_coord :: Vector{<:Real},
    filament_strength :: Real)

    @assert(length(filament_start_coord)==3,
        "Filament start coordinate vector must have length 3.")
    @assert(length(filament_end_coord)==3,
        "Filament end coordinate vector must have length 3.")
    convertable_to_F32(filament_strength, "filament_strength")
    convertable_to_F32(filament_start_coord, "filament_start_coord")
    convertable_to_F32(filament_end_coord, "filament_end_coord")
    return
end

function check_filament_definition(
    filament_start_coord :: Matrix{<:Real},
    filament_end_coord :: Matrix{<:Real},
    filament_strength :: Vector{<:Real})

    @assert(size(filament_start_coord)[2]==3,
        "Filament start coordinate matrix must be N by 3. Actual size is "*
        string(size(filament_start_coord))*".")
    @assert(size(filament_end_coord)[2]==3,
        "Filament end coordinate matrix must be N by 3. Actual size is "*
        string(size(filament_end_coord))*".")
    @assert(size(filament_start_coord)[1]==size(filament_end_coord)[1],
        "Both filament start coordinate and filament end coordinate "*
        "matrices must be N by 3 (ie. of equal size), but the sizes "*
        "don't match. filament_start_coord matrix is "*
        string(size(filament_start_coord))*" and the end coordiante matrix is of"*
        "size "*string(size(filament_end_coord))*".")
    @assert(size(filament_start_coord)[1]==length(filament_strength),
        "The length of the filament strength vector does not match that "*
        "of the coordinate matrices. The filament coordinate matrices define"*
        " geometries for "* string(size(filament_start_coord)[1]) *
        " filaments whilst the strength vector has length " *
        string(length(filament_strength)) * ".")
    convertable_to_F32(filament_strength, "filament_strength")
    convertable_to_F32(filament_start_coord, "filament_start_coord")
    convertable_to_F32(filament_end_coord, "filament_end_coord")
    return
end

function check_filament_definition(
    filament_start_coord :: Matrix{<:Real},
    filament_end_coord :: Matrix{<:Real})

    @assert(size(filament_start_coord)[2]==3,
        "Filament start coordinate matrix must be N by 3. Actual size is "*
        string(size(filament_start_coord))*".")
    @assert(size(filament_end_coord)[2]==3,
        "Filament end coordinate matrix must be N by 3. Actual size is "*
        string(size(filament_end_coord))*".")
    @assert(size(filament_start_coord)[1]==size(filament_end_coord)[1],
        "Both filament start coordinate and filament end coordinate "*
        "matrices must be N by 3 (ie. of equal size), but the sizes "*
        "don't match. filament_start_coord matrix is "*
        string(size(filament_start_coord))*" and the end coordiante matrix is of"*
        "size "*string(size(filament_end_coord))*".")
    convertable_to_F32(filament_start_coord, "filament_start_coord")
    convertable_to_F32(filament_end_coord, "filament_end_coord")
    return
end

"""
    Gives an assert message if something cannot be converted to a
    CVortex.VortexParticle or Vector{CVortex.VortexParticle}.

Part of CVortex. Not for client use.
"""
function check_particle_definition_3D(
    particle_coords :: Matrix{<:Real},
    particle_vorts  :: Matrix{<:Real})

    convertable_to_Vec3f_vect(particle_coords, "particle_coords")
    convertable_to_Vec3f_vect(particle_vorts, "particle_vorts")
    @assert(size(particle_coords)==size(particle_vorts),
        "The number of particle defined by particle_coords does not "*
        "match that defined by particle_vorts. particle_coords defines "*
		string(size(particle_coords)[1])*" particle coordinates whilst "*
        "particle_vorts defines "*string(size(particle_vorts))*
        " particle vorticities.")
    return
end

function check_particle_definition_3D(
    particle_coords :: Matrix{<:Real},
    particle_vorts  :: Matrix{<:Real},
    particle_vols :: Vector{<:Real})

    convertable_to_Vec3f_vect(particle_coords, "particle_coords")
    convertable_to_Vec3f_vect(particle_vorts, "particle_vorts")
    convertable_to_F32(particle_vols, "particle_vols")
    @assert(size(particle_coords)==size(particle_vorts),
        "The number of particle defined by particle_coords does not "*
        "match that defined by particle_vorts. particle_coords defines "*
		string(size(particle_coords)[1])*" particle coordinates whilst "*
        "particle_vorts defines "*string(size(particle_vorts))*
        " particle vorticities.")
    @assert(length(particle_vols)==size(particle_coords)[1],
        "The number of particles defined by the particle volumes vector"*
        " does not match that of the coordinate and vorticity matrices. "*
        "The volume vector is of length "*string(length(particle_vols))*
        " and particle_coords suggests "*string(size(particle_coords)[1])*
        " particles.")
    return
end

function check_particle_definition_3D(
    particle_coords :: Vector{<:Real},
    particle_vorts  :: Vector{<:Real})

    convertable_to_Vec3f_vect(particle_coords, "particle_coords")
    convertable_to_Vec3f_vect(particle_vorts, "particle_vorts")
    return
end

function check_particle_definition_3D(
    particle_coords :: Vector{<:Real},
    particle_vorts  :: Vector{<:Real},
    particle_vol :: Real)

    convertable_to_Vec3f_vect(particle_coords, "particle_coords")
    convertable_to_Vec3f_vect(particle_vorts, "particle_vorts")
    convertable_to_F32(particle_vol, "particle_vol")
    return
end

"""
    Gives an assert message if something cannot be converted to a
    CVortex.VortexParticle2D or Vector{CVortex.VortexParticle2D}.

Part of CVortex. Not for client use.
"""
function check_particle_definition_2D(
    particle_coords :: Matrix{<:Real},
    particle_vorts  :: Vector{<:Real})

    convertable_to_Vec2f_vect(particle_coords, "particle_coords")
    convertable_to_F32(particle_vorts, "particle_vorts")
    @assert(size(particle_coords)[1]==length(particle_vorts),
        "The number of particle defined by particle_coords does not "*
        "match that defined by particle_vorts. particle_coords defines "*
		string(size(particle_coords)[1])*" particle coordinates whilst "*
        "particle_vorts defines "*string(length(particle_vorts))*
        " particle vorticities.")
    return
end

function check_particle_definition_2D(
    particle_coords :: Matrix{<:Real},
    particle_vorts  :: Vector{<:Real},
    particle_vols :: Vector{<:Real})

    convertable_to_Vec2f_vect(particle_coords, "particle_coords")
    convertable_to_F32(particle_vorts, "particle_vorts")
    convertable_to_F32(particle_vols, "particle_vols")
    @assert(size(particle_coords)[1]==length(particle_vorts),
        "The number of particle defined by particle_coords does not "*
        "match that defined by particle_vorts. particle_coords defines "*
		string(size(particle_coords)[1])* " particle coordinates whilst "*
        "particle_vorts defines "*string(length(particle_vorts))*
        " particle vorticities.")
    @assert(length(particle_vols)==size(particle_coords)[1],
        "The number of particles defined by the particle volumes vector"*
        " does not match that of the coordinate and vorticity matrices. "*
        "The volume vector is of length ", length(particle_vols),
        " and particle_coords suggests "*string(size(particle_coords)[1])*
        " particles.")
    return
end

function check_particle_definition_2D(
    particle_coords :: Vector{<:Real},
    particle_vort  :: Real)

    convertable_to_Vec2f_vect(particle_coords, "particle_coords")
    convertable_to_F32(particle_vort, "particle_vorts")
    return
end

function check_particle_definition_2D(
    particle_coords :: Vector{<:Real},
    particle_vort  :: Real,
    particle_vol :: Real)

    convertable_to_Vec2f_vect(particle_coords, "particle_coords")
    convertable_to_F32(particle_vort, "particle_vort")
    convertable_to_F32(particle_vol, "particle_vol")
    return
end


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__
cd(ProjDir) #do

fname = scm_path("..", "examples", "Shipley", "fig2.6", "fig2.6.dot")
#Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

# Read `=>` as `~` in regression models, or `<-` in causal models.
d = OrderedDict(
  :u => [:x, :v],
  :s1 => [:u],
  :w => [:v, :y],
  :s2 => [:w]
);

dag = DAG("fig2.6", d);

@testset "m_separation_02" begin

  @test to_ggm(dag) == "DAG(u ~ x + v, s1 ~ u, w ~ v + y, s2 ~ w)"
  @test m_separation(dag, :x, :v; debug=false) == true
  @test m_separation(dag, :x, [:v]; c=:u, debug=false) == false
  @test m_separation(dag, [:x], :v; c=:s1, debug=false) == false
  @test m_separation(dag, [:u], [:w]; debug=false) == false
  @test m_separation(dag, :u, [:w]; c=:v, debug=false) == true

  @test m_separation(dag, [:x], :y) == true
  @test m_separation(dag, :x, :y; c=[:u, :w], debug=false) == false
  @test m_separation(dag, :x, :y; c=[:s1, :s2], debug=false) == false
  @test m_separation(dag, :x, :y; c=[:u, :v, :w], debug=false) == true
  @test m_separation(dag, :x, [:v]; c=:u, debug=false) == false
  @test m_separation(dag, [:x], :v; c=:s1, debug=false) == false

end


text:
        # Exponential cone constrained problems

using Test, LinearAlgebra, SparseArrays
using COSMO

# This test is precision agnostic and will be run with T precision
if !@isdefined(UnitTestFloats)
    UnitTestFloats = [Float64] #if not run in full test setup, just do it for one float type
end

@testset "Exponential cone problems" begin

for T in UnitTestFloats

    @testset "Exponential cone problems (T = $(T))" begin

        @testset "Feasible" begin

            # solve the following exponential cone problem
            # max  x
            # s.t. y * exp(x / y) <= z
            #      y == 1, z == exp(5)

            P = spzeros(T, 3, 3)
            q = T[-1.; 0; 0]

            # y * exp( x / y) <= z
            A1 = spdiagm(0 => ones(T, 3))
            b1 = zeros(T, 3)

            # y == 1 and z == exp(5) = 148.4131591025766
            A2 = SparseMatrixCSC{T, Int}([0 1. 0; 0 0 1])
            b2 = [-one(T); -exp(T(5))]
            cs1 = COSMO.Constraint(A1, b1, COSMO.ExponentialCone)
            cs2 = COSMO.Constraint(A2, b2, COSMO.ZeroSet)

            model = COSMO.Model{T}()
            assemble!(model, P, q, [cs1; cs2], settings = COSMO.Settings{T}(eps_abs = T(1e-4), eps_rel = T(1e-4)))

            res = COSMO.optimize!(model)
            @test res.status == :Solved
            @test isapprox(res.obj_val, T(-5), atol=1e-2)
        end

        # Run infeasibility tests only if precision(T) is high enough
        if precision(T) >= precision(Float64)
            @testset "Primal infeasible 1" begin
                # min  x
                # s.t. y * exp( x / y) <= z
                #       y == 1
                #       z == -1

                P = spzeros(T, 3, 3)
                q = T[1.; 0; 0]

                # y * exp( x / y) <= z
                A1 = spdiagm(0 => ones(T, 3))
                b1 = zeros(T, 3)

                # y == 1
                A2 = T[0 -1. 0]
                b2 = T[-1.]

                # z == -1
                A3 = T[0 0 -1.]
                b3 = T[1]

                cs1 = COSMO.Constraint(A1, b1, COSMO.ExponentialCone)
                cs2 = COSMO.Constraint(A2, b2, COSMO.ZeroSet)
                cs3 = COSMO.Constraint(A3, b3, COSMO.ZeroSet)

                model = COSMO.Model{T}()
                assemble!(model, P, q, [cs1; cs2; cs3])

                res = COSMO.optimize!(model)
                @test res.status == :Primal_infeasible
            end

            @testset "Primal infeasible 2" begin
                # min  x
                # s.t. y * exp( x / y) <= z - 0.2
                #      y * exp( x / y) >= z + 0.2

                P = spzeros(T, 3, 3)
                q = T[1.; 0; 0]

                # y * exp( x / y) <= z - 0.2
                A1 = spdiagm(0 => ones(T, 3))
                b1 = T[0.; 0; -0.2]

                # y * exp( x / y) >= z + 0.3
                A2 = -spdiagm(0 => ones(T, 3))
                b2 = T[0.; 0; -0.3]

                cs1 = COSMO.Constraint(A1, b1, COSMO.ExponentialCone)
                cs2 = COSMO.Constraint(A2, b2, COSMO.ExponentialCone)

                model = COSMO.Model{T}()
                assemble!(model, P, q, [cs1; cs2])

                res = COSMO.optimize!(model)
                @test res.status == :Primal_infeasible
            end

            @testset "Dual infeasible" begin
                # max z
                # s.t. y * exp(x / y) <= z

                P = spzeros(T, 3, 3)
                q = T[0; 0; -1.0]

                # y * exp( x / y) <= z
                A1 = spdiagm(0 => ones(T, 3))
                b1 = zeros(T, 3)

                cs1 = COSMO.Constraint(A1, b1, COSMO.ExponentialCone)

                model = COSMO.Model{T}()
                assemble!(model, P, q, cs1)

                res = COSMO.optimize!(model)
                @test res.status == :Dual_infeasible
            end
        end
    end

    @testset "Dual Exponential cone problems (T = $(T))" begin
        if precision(T) >= precision(Float64)
        @testset "Feasible" begin

            # solve the following dual exponential cone problem
            # min  y
            # s.t. -x * exp(y / x) <= exp(1) * z
            #      x == -1, z == exp(5)

            P = spzeros(T, 3, 3)
            q = T[0; 1.; 0]

            # -x * exp( y / x) <= z
            A1 = spdiagm(0 => ones(T, 3))
            b1 = zeros(T, 3)
            cs1 = COSMO.Constraint(A1, b1, COSMO.DualExponentialCone)

            # y == 1 and z == exp(5)
            A2 = SparseMatrixCSC{T, Int}([1. 0 0; 0 0 1])
            b2 = [one(T); -exp(T(5))]
            cs2 = COSMO.Constraint(A2, b2, COSMO.ZeroSet)

            model = COSMO.Model{T}()
            assemble!(model, P, q, [cs1; cs2], settings = COSMO.Settings{T}())

            res = COSMO.optimize!(model)
            @test res.status == :Solved
            @test isapprox(res.obj_val, T(-6), atol=1e-3)
        end
    end
        # Run infeasibility tests only if precision(T) is high enough
        if precision(T) >= precision(Float64)
            @testset "Primal Infeasible" begin

                # solve the following dual exponential cone problem
                # min  u + v + w
                # s.t. -u * exp(v / u) <= exp(1) * w
                #      u == 1, v == 2

                P = spzeros(T, 3, 3)
                q = ones(T, 3)

                # -u * exp( v / u) <= exp(1) * w
                A1 = spdiagm(0 => ones(T, 3))
                b1 = zeros(T, 3)
                cs1 = COSMO.Constraint(A1, b1, COSMO.DualExponentialCone)

                #  u == 1 and b == 2
                A2 = SparseMatrixCSC{T, Int}([1. 0 0; 0 1 0])
                b2 = T[-1.; -2]
                cs2 = COSMO.Constraint(A2, b2, COSMO.ZeroSet)

                model = COSMO.Model{T}()
                assemble!(model, P, q, [cs1; cs2], settings = COSMO.Settings{T}())

                res = COSMO.optimize!(model)
                @test res.status == :Primal_infeasible

            end
        end
    end
end
end
nothing


text:
using Rimu
using Rimu.RMPI
using KrylovKit
using DistributedDictVectors

add = BoseFS((0,0,0,0,0,11,0,0,0,0,0))
H = HubbardMom1D(add; u=6.0)

if mpi_size() == 1 && Threads.nthreads() == 1
    dv = DVec(add => 1.0)
    println("crunch single")
    @time res = eigsolve(H, dv, 1, :SR; issymmetric=true)[1][1]
    println("result is $res")
elseif mpi_size() == 1
    tv = TVec(add => 1.0)
    println("crunch threads")
    @time res = eigsolve(H, tv, 1, :SR)[1][1]
    println("result is $res")
else
    tv = TVec(add => 1.0)
    if is_mpi_root()
        println("crunch mpi $(mpi_size())")
        @time res = eigsolve(H, tv, 1, :SR)[1][1]
        println("result is $res")
    else
        eigsolve(H, tv, 1, :SR)[1][1]
    end
end


text:
### Degree sequence functions

function find_startindex(degseq_rotated::Vector{<:Integer})
    min_value = Inf64
    min_index = 0
    current_sum = 0
    size = length(degseq_rotated)
    for i = 1:length(degseq_rotated)
        current_sum += degseq_rotated[i]-1
        if current_sum < min_value
            min_value = current_sum
            min_index = i
        end
    end
    return (min_index+1) % size
end

function fix_degseq(degseq_rotated::Vector{<:Integer})
    start_index = find_startindex(degseq_rotated)
    return circshift(degseq_rotated, -(start_index-1))
end

function degrees_count(size::Int, dist::DiscreteUnivariateDistribution, maxdegree::Int)
    buf = 16 # This is enough for almost all real applications
    degcount = Vector{Int}(undef, buf)

    while true
        cur = 0
        total_node = 0
        total_deg = 0
        p_cum = 0.0
        p_cur = 0.0
        p_para::Float64 = 0.0
        current_degcount = 0
        deg = 0

        while total_deg < size - 1 && total_node < size && deg <= maxdegree
            p_cur = pdf(dist, deg)
            p_para = p_cur/(1-p_cum)
            # This happens at the boundary
            if p_para > 1.0 || p_para < 0.0
                p_para = 1.0
            end
            current_degcount = rand(Binomial(size-total_node, p_para))
            total_node += current_degcount
            total_deg += deg*current_degcount
            if deg < length(degcount)
                degcount[deg+1] = current_degcount
            else
                push!(degcount, current_degcount)
            end
            p_cum += p_cur
            deg += 1
        end

        if total_deg == size - 1 && total_node == size
            return degcount[1:deg]
        end
    end
end

function degrees(size::Int, degcount::AbstractArray{T, 1}) where T <: Integer
    degseq = Array{Int, 1}(undef, size)
    pos = 1
    for deg in 0:length(degcount)-1
        n = degcount[deg+1]
        degseq[pos:pos+n-1] .= deg
        pos = pos + n
    end
    degseq
end

### Cayley Tree

"A conditional Galton-Watson tree with Poission(1) offspring distrbiution."
struct CayleyTree <: CondGWTree
    spec::TreeSpec
end
"""
    CayleyTree(size::Int)

Construct a `CayleyTree` of size `size`.
"""
CayleyTree(size::Int) = CayleyTree(TreeSpec(size, "Cayley Tree"))

function degrees_rotate(tree::CayleyTree, nbatch=2048)
    # See Svante Example 12.1 and Luc's simulation paper
    treesize = size(tree)
    degseq = zeros(Int32, treesize)
    batch = Vector{Int}(undef, nbatch)
    i = 0
    while i + nbatch < treesize
        rand!(batch, 1:treesize)
        for j = 1:nbatch
            @inbounds degseq[batch[j]] += 1
        end
        i += nbatch
    end
    for j = i+1:treesize-1
        box = rand(1:treesize)
        degseq[box] += 1
    end
    degseq
end

function degrees(tree::CondGWTree)
    degseq = degrees_rotate(tree)
    fix_degseq(degseq)
end

### Binary Tree

"A conditional Galton-Watson tree with `2 Bernoulli(1/2)` offspring distrbiution."
struct BinaryTree <: CondGWTree
    spec::TreeSpec
end

"""
    BinaryTree(size::Int)

Construct a `BinaryTree` of size `size`.
"""
function BinaryTree(size::Int)
    if iseven(size)
        size += 1
    end
    BinaryTree(TreeSpec(size, "Binary Tree"))
end

function degrees_rotate(tree::BinaryTree)
    treesize = size(tree)
    degseq = Array{Int32}(undef, treesize)
    lefanum = convert(Int, (treesize + 1)/2)

    # initialize the array
    degseq[1:lefanum] .= 0
    degseq[lefanum+1:end] .= 2

    shuffle!(degseq)
end

### General conditioned Galton-Watson tree

"A conditional Galton-Watson tree."
struct GeneralCondGWTree <: CondGWTree
    "Specification"
    spec::TreeSpec
    "Offspring distribution"
    dist::DiscreteUnivariateDistribution
end

distribution(tree::GeneralCondGWTree) = tree.dist

"""
    GeneralCondGWTree(size::Int, maxdegree::Int, dist::DiscreteUnivariateDistribution)

Construct a `GeneralCondGWTree` of size `size` and offspring distribution `dist`.
"""
function GeneralCondGWTree(size::Int, maxdegree::Int, dist::DiscreteUnivariateDistribution)
    GeneralCondGWTree(TreeSpec(size, "Conditional Galton-Watson Tree with offspring distribution $dist", maxdegree), dist)
end

"""
    GeneralCondGWTree(size::Int, maxdegree::Int, dist::DiscreteUnivariateDistribution)

Construct a `GeneralCondGWTree` of size `size`, maximum degree `maxdegree`, and offspring
distribution `dist`.
"""
function GeneralCondGWTree(size::Int, dist::DiscreteUnivariateDistribution)
    GeneralCondGWTree(size, size, dist)
end

"""
    CatalanTree(size::Int)

Construct a `GeneralCondGWTree` with `Geometric(1/2)` offspring distribution of size `size`.
"""
CatalanTree(size::Int) = GeneralCondGWTree(TreeSpec(size, "Catalan Tree"), Geometric(0.5))

"""
    DAryTree(size::Int, d::Int)

Construct a `GeneralCondGWTree` with `Binomial(d, 1/d)` offspring distribution of size `size`.
"""
DAryTree(size::Int, d::Int) = GeneralCondGWTree(TreeSpec(size, "$d-ary Tree", d), Binomial(d, 1/d))

"""
    MotzkinTree(size::Int)

Construct a `GeneralCondGWTree` with `DiscreteUniform(0, 2)` offspring distribution of size `size`.
"""
MotzkinTree(size::Int) = GeneralCondGWTree(TreeSpec(size, "Motzkin Tree", 2), DiscreteUniform(0, 2))

function degrees_rotate(tree::GeneralCondGWTree)
    degcount = degrees_count(size(tree), distribution(tree), maxdegree(tree))
    degseq = degrees(size(tree), degcount)
    shuffle!(degseq)
end

"""
    treegraph(tree)::FixedTreeGraph

Either convert a fixed `tree` to a `FixedTreeGraph`, or generate a random `FixedTreeGraph` according
to `tree`.
"""
function treegraph end

function treegraph(tree::FiniteTree)::FixedTreeGraph
    degseq = degrees(tree)
    walker = GraphWalker(size(tree))
    walk(degseq, walker, true)
    walker.tree_digraph
end

function treegraph(degseq::Vector{Int})::FixedTreeGraph
    walker = GraphWalker(length(degseq))
    walk(degseq, walker, true)
    walker.tree_digraph
end

### RandomRecursiveTree

"""
A random recursive tree.
"""
struct RandomRecursiveTree <: FiniteRandomTree
    spec::TreeSpec
end

"""
    RandomRecursiveTree(size::Int)

Construct a `RandomRecursiveTree` of size `size`.
"""
RandomRecursiveTree(size::Int) = RandomRecursiveTree(TreeSpec(size, "Random Recursive Tree"))

"""
    degrees(tree)

Return an iterator that traverses the degree sequence of `tree` in DFS order.
"""
function degrees end

function degrees(tree::RandomRecursiveTree)
    return treegraph(tree)
end

function treegraph(tree::RandomRecursiveTree)::FixedTreeGraph
    treesize = size(tree)
    graph = FixedTreeGraph(treesize)

    for current_node in 2:treesize
        parent_node = rand(1:current_node-1)
        addedge!(graph, parent_node, current_node)
    end

    graph
end


text:
include("UAVLanding.jl")
include("visual.jl")

POMCPOW_solver = POMCPOWSolver(max_depth = 1000, # the deeper, the more oscillations on heading angles.
                     tree_queries=1000,

                     criterion=MaxUCB(10.0),
                     enable_action_pw=true,
                     check_repeat_obs=false,
                     rng=MersenneTwister(2))

my_pomdp = UAVchasePOMDP()
policyPOMCPOW = POMCPOWPlanner(POMCPOW_solver, my_pomdp)
hr = HistoryRecorder(show_progress=true, max_steps=400)
up = SIRParticleFilter(my_pomdp, 400)
hist = simulate(hr, my_pomdp, policyPOMCPOW, up)

# view the simulation history:
# for (s, b, a, r, sp, o) in view(hist, 145:145)
#     @show s, a, r, o
# end
print(" Number of steps: ")
print(length(hist))

# plot a static plot showing the UAV and target trace:
# plt_static = plot3d(my_pomdp, hist)

# create a gif
gr()
frames = Frames(MIME("image/png"), fps=5)
@showprogress "Creating gif..." for i in 1:30
    push!(frames, plot3d(my_pomdp, view(hist, 1:i)))
end
write("out.gif", frames)


text:
using HighPrecisionNonsymmetricEigenproblem
using Test
using MultiFloats
using LinearAlgebra

using HighPrecisionNonsymmetricEigenproblem: pack!, do_qr!, unpack!, apply_tiny_Q_left!

@testset "HighPrecisionNonsymmetricEigenproblem.jl" begin
    T = Float64x2
    Q = rand(T, 32, 4)
    t = rand(T, 4)
    A = rand(T, 32, 9968)

    Q_packed, t_packed, A_packed = pack!(Q, t, A)
    do_qr!(Q_packed, t_packed, A_packed, T)
    A_fast_kernel = unpack!(similar(A), A_packed)
    A_slow_kernel = apply_tiny_Q_left!(Q, copy(A), t)

    @test norm(A_fast_kernel - A_slow_kernel) < 1e-20
end


text:
x = range(0, stop=2, length=25)
for α ∈ 0.01:0.05:0.8
    plot!(x, x.^α, lwidth=α/10, col=RGB(0.0,0.0,α), smooth=true)
end


text:
using PkgBenchmark
using Query
using QueryOperators
using DataTables

@benchgroup "Variable group columns" begin
    N = 100_000;
    A = rand(N);
    B = rand(1:100, N);
    C = rand(1:100, N);
    dt = DataTable([A, B, C], [:A, :B, :C]);

    @bench "one column" @from i in $dt begin
        @group i.A by i.B into g
        @select {m = mean(g)}
        @collect
    end

    @bench "two columns" @from i in $dt begin
        @group {i.A, i.B} by i.B into g
        @select {m = mean(g.A)}
        @collect
    end

    @bench "three columns" @from i in $dt begin
        @group {i.A, i.B, i.C} by i.B into g
        @select {m = mean(g.A)}
        @collect
    end
end

@benchgroup "Variable group by columns" begin
    N = 100_000;
    A = rand(N);
    B = rand(1:100, N);
    C = rand(1:100, N);
    dt = DataTable([A, B, C], [:A, :B, :C]);

    @bench "one column" @from i in $dt begin
        @group i.A by i.A into g
        @select {m = mean(g)}
        @collect DataTable
    end

    @bench "two columns" @from i in $dt begin
        @group i.A by {i.A, i.B} into g
        @select {m = mean(g)}
        @collect DataTable
    end

    @bench "three columns" @from i in $dt begin
        @group i.A by {i.A, i.B, i.C} into g
        @select {m = mean(g)}
        @collect DataTable
    end
end

@benchgroup "Variable datatypes" begin
    @benchgroup "Integer" begin
        N = 10_000;
        A = rand(1:10_000, N);
        B = rand(1:10_000, N);
        C = rand(1:10_000, N);
        dt = DataTable([A, B, C], [:A, :B, :C]);

        @bench "IntegerTest" @from i in $dt begin
            @group i.A by i.A into g
            @select g
            @collect
        end
    end

    @benchgroup "Float64" begin
        N = 10_000;
        A = rand(N);
        B = rand(N);
        C = rand(N);
        dt = DataTable([A, B, C], [:A, :B, :C]);

        @bench "Float64Test" @from i in $dt begin
            @group i.A by i.A into g
            @select g
            @collect
        end
    end

    @benchgroup "String" begin
        N = 10_000;
        A = [randstring(20) for i in 1:N];
        B = [randstring(20) for i in 1:N];
        C = [randstring(20) for i in 1:N];
        dt = DataTable([A, B, C], [:A, :B, :C]);

        @bench "StringTest" @from i in $dt begin
            @group i.A by i.A into g
            @select g
            @collect
        end
    end

    @benchgroup "Date" begin
        N = 10_000;
        A = [Dates.unix2datetime(time() - i) for i in rand(-1_000_000_000:1_000_000_000, N)];
        B = [Dates.unix2datetime(time() - i) for i in rand(-1_000_000_000:1_000_000_000, N)];
        C = [Dates.unix2datetime(time() - i) for i in rand(-1_000_000_000:1_000_000_000, N)];
        dt = DataTable([A, B, C], [:A, :B, :C]);

        @bench "DateTest" @from i in $dt begin
            @group i.A by i.A into g
            @select g
            @collect
        end
    end
end


text:
using TimeseriesSurrogates

export bbnue

# Keep this for when a common interface for optimized variable selection methods has been established
# """
#export BBNUE
#     BBNUE(est) <: TransferEntropyEstimator

# The bootstrap-based non-uniform embedding estimator (BB-NUE) for conditional transfer entropy (Montalto et al., 2014).
# Uses the estimator `est` to compute relevant marginal entropies. (e.g. `VisitationFrequency(RectangularBinning(3))`)

# [^Montalto2014]: Montalto, A.; Faes, L.; Marinazzo, D. MuTE: A MATLAB toolbox to compare established and novel estimators of the multivariate transfer entropy. PLoS ONE 2014, 9, e109462.
# """
# struct BBNUE{E} <: TransferEntropyEstimator
#     est::E
# end

"""
    bbnue(source, target, [cond], est;
        η = 1, include_instantaneous = true,
        method_delay = "ac_min", maxlag::Union{Int, Float64} = 0.05,
        surr::Surrogate = RandomShuffle(), nsurr = 100, α = 0.05)
        ) → te, js, τs, idxs_source, idxs_target, idxs_cond

Estimate transfer entropy (TE) from `source` to `target` (conditioned on `cond` if given)
for prediction lag `η`, using the bootstrap-based non-uniform embedding (BBNUE)
method from Montalta et al. (2004) [^Montalto2014].

## Implementation details

The BBNUE method uses a bootstrap-based criterion to identify the most relevant and
minimally redundant variables from the the past of `target`, present/past of `source`,
and (if given) the present/past of `cond`, that contribute most to `target`'s future.

This implementation uses a conditional entropy minimization criterion for selecting variables,
which is what Montalto et al. (2014)[^Montalto2014] uses for their bin-estimator. Here, any
estimator can be used, but variables will be selected using conditional entropy minimization
regardless of the choice of estimator.

Montalto et al.'s bin-estimator corresponds to using the `VisitationFrequency` estimator with bins
whose sides are equal-length, e.g. `VisitationFrequency(RectangularBinning(0.5))`.
In this implementation, any rectangular binning can be used.

## Input data

Multivariate `source`, `target` and `cond` (if given) can be given as univariate
`AbstractVector`s or as multivariate `Dataset`s or `Vector{AbstractVector}`.

For example, if you want to compute
the BBNUE-transfer entropy from a univariate source to a univariate target,
potentially conditioned on many different variables, you can do the following:

```julia
n = 1000
# Source and target variables
s, t = rand(n), rand(n)

# Variables that might potentially influence `t` along with `s`
c1, c2, c3 = rand(n), rand(n), rand(n)

est = NaiveKernel(0.3)
bbnue(s, t, Dataset([c1, c2, c3]), est)
```

## Variable selection and significance testing

In this implementation, the maximum lag for each embedding variable is determined using `estimate_delay`
from `DelayEmbeddings`. The keywords `method_delay` (the default is `"ac_min"`) controls the method
for estimating the delay, and `maxlag` is the maximum allowed delay (if `maxlag ∈ [0, 1]` is a fraction,
then the maximum lag is that fraction of the input time series length, and if `maxlag` is an integer,
then the maximum lag is `maxlag`).

If `instantaneous` is `true`, then instantaneous interactions are also considered, i.e. effects like
`source(t) → target(t)` are allowed. `η` is the forward prediction lag.

Significance testing is performed using a permutation test. At each iteration of the
variable selection routine, we first compute the transfer entropy using the new candidate
``c_k``. Then, the computation is repeated `nsurr` times, at each iteration replacing ``c_k``
with a surrogate of type `surr`. If transfer entropy using the original ``c_k`` exceeds the
the `1 - α`-quantile of that of the surrogate ensemble, then ``c_k`` is deemed significant
to the future of `target` and is included in the set of selected variables.

If no relevant variables pass the permutation test, then TE is not well-defined, and a value of `0.0`
is returned.

## Returns

A 6-tuple is returned, consisting of:
- `te`: The computed transfer entropy value.
- `js`: The indices of the selected variables. `js[i]` is the `i`-th entry in the array `[idxs_source..., idxs_target..., idxs_cond...,]`.
- `τs`: The embedding lags of the selected variables. `τs[i]` corresponds to `js[i]`.
- `idxs_source`: The indices of the source variables.
- `idxs_target`: The indices of the target variables.
- `idxs_cond`: The indices of the conditional variables (empty if `cond` is not given).

## Example

```julia
using CausalityTools, DynamicalSystems
sys = ExampleSystems.logistic2_unidir(c_xy = 0.8, r₁ = 3.78, r₂ = 3.92)
orbit = trajectory(sys, 10000, Ttr = 10000)
x, y = columns(orbit)

# Use a coarse-grained rectangular binning with subdivisions in each dimension,
# to keep computation costs low and to ensure the probability distributions
# over the bins don't approach the uniform distribution (need enough points
# to fill bins).
est = NaiveKernel(0.3)
te_xy = bbnue(x, y, est, surr = RandomShuffle(), nsurr = 100, include_instantaneous = true)
te_yx = bbnue(y, x, est, surr = RandomShuffle(), nsurr = 100, include_instantaneous = true)

te_xy, te_yx
```

[^Montalto2014]: Montalto, A.; Faes, L.; Marinazzo, D. MuTE: A MATLAB toolbox to compare established and novel estimators of the multivariate transfer entropy. PLoS ONE 2014, 9, e109462.
"""
function bbnue(source, target, cond, est; q = 1, base = 2,
        η = 1, include_instantaneous = true,
        method_delay = "ac_min", maxlag::Union{Int, Float64} = 0.05,
        α = 0.05, nsurr = 19, surr::TimeseriesSurrogates.Surrogate = RandomShuffle())

    Ω, Y⁺, τs, js, idxs_source, idxs_target, idxs_cond =
        candidate_embedding(
            process_input(source),
            process_input(target),
            process_input(cond),
            η = η,
            include_instantaneous = include_instantaneous,
            method_delay = method_delay,
            maxlag = maxlag)

    return optim_te(Ω, Y⁺, τs, js, idxs_source, idxs_target, idxs_cond, est,
        q = q, base = base, α = α, nsurr = nsurr, surr = surr)
end

function bbnue(source, target, est; q = 1, base = 2,
        η = 1, include_instantaneous = true,
        method_delay = "ac_min", maxlag::Union{Int, Float64} = 0.05,
        α = 0.05, nsurr = 19, surr::TimeseriesSurrogates.Surrogate = RandomShuffle())

    Ω, Y⁺, τs, js, idxs_source, idxs_target, idxs_cond =
        candidate_embedding(
            process_input(source),
            process_input(target),
            η = η,
            include_instantaneous = include_instantaneous,
            method_delay = method_delay,
            maxlag = maxlag)

    return optim_te(Ω, Y⁺, τs, js, idxs_source, idxs_target, idxs_cond, est,
        q = q, base = base, α = α, nsurr = nsurr, surr = surr)
end


text:
module ExampleModeling

include("model.jl")
export ExampleModel,
    Specs,
    Indices,
    Params,
    Variables,
    Objectives,
    model_to_dtype

include("io.jl")
export save_json,
    load_json

end # module

max_stars_repo_name: nnhjy/JuMP-project-template

max_issues_repo_name: nnhjy/JuMP-project-template

max_forks_repo_name: nnhjy/JuMP-project-template


text:
struct TimeLift{T} <: DynamicIterator
    iter::T
end

dyniterate(iter::TimeLift, ::Nothing) = dyniterate(iter, 0=>nothing)

function dyniterate(TL::TimeLift, (i, state)::Pair)
    x, state = @returnnothing dyniterate(TL.iter, state)
    i + 1 => x, (i + 1 => state)
end

function dyniterate(TL::TimeLift, (state, j)::NextKey)
    x, state = @returnnothing dyniterate(TL.iter, state)
    j => x, (j => state)
end
function dyniterate(TL::TimeLift, ((i, state), j)::NextKey{<:Pair})
    x, state = @returnnothing dyniterate(TL.iter, Steps(state, j-i))
    j => x, (j => state)
end
function dyniterate(TL::TimeLift, (state,)::NextKeys, key)
    u, state = @returnnothing dyniterate(TL, NextKey(state, key))
    u, NextKeys(state)
end
function dyniterate(TL::TimeLift, (nk,)::Start{<:NextKeys}, key)
    ϕ = dyniterate(TL, NextKey(nk[], key))
    u, state = @returnnothing ϕ
    u, NextKeys(state)
end

function dyniterate(TL::TimeLift, ((i,state), j)::NewKey{<:Pair})
    x, state = @returnnothing dyniterate(TL.iter, state)
    j => x, (j => state)
end

"""
    dyniterate(iter, state, (steps,)::Steps)

Advance the iterator `steps` times, and for negative
numbers, if implemented, rewind the iterator `-steps`
times.
"""
function dyniterate(iter, (state, n)::Steps)
    @assert n ≥ 1
    local x
    for k in 1:n
        x, state = @returnnothing dyniterate(iter, state)
    end
    x, state
end
function dyniterate(iter, (state,n)::Steps{Nothing})
    ϕ = @returnnothing dyniterate(iter)
    n == 1 && return ϕ
    dyniterate(iter, Steps(ϕ[2], n-1))
end
function dyniterate(iter::GEvolution, (x,n)::Steps)
    @assert n ≥ 0
    for k in 1:n
        x = @returnnothing evolve(iter, x)
    end
    dub(x)
end

Base.IteratorSize(tl::TimeLift) = Base.IteratorSize(tl.iter)


text:
# setting working directory
# cd("/Volumes/SSD Hans/Github/MacroFall2020/AdvMacro")
# mkpath("Assignment3")

# mkpath("graphs")
using Plots
# using LateXStrings # Pkg.add("LaTeXStrings") # https://github.com/stevengj/LaTeXStrings.jl
Pkg.add("Dierckx")
using Dierckx # Pkg.add("Dierckx") # https://github.com/kbarbary/Dierckx.jl
Pkg.add("Interpolations")
using Interpolations # Pkg.add("Interpolations") # https://github.com/JuliaMath/Interpolations.jl
Pkg.add("ForwardDiff")
using ForwardDiff # Pkg.add("ForwardDiff") # https://github.com/JuliaDiff/ForwardDiff.jl
Pkg.add("LaTeXStrings")
using LaTeXStrings
Pkg.add("Latexify")
using Latexify
Pkg.add("PyPlot")

cd("/Volumes/SSD Hans/Github/MacroFall2020/AdvMacro/Assignment3")

# define Functions

u1(x) = log(x)
u2(x) = (x)^(1/2)
u3(x,s) = x^(1-s)/(1-s)
u4(x) = u3(x,2)
u5(x) = u3(x,5)
u6(x) = u3(x,10)

# define derivatives
u1p(x) = 1/x
u2p(x) = (1/2)*x^(-1/2)
u3p(x,s) = x^(-s)
u4p(x) = x^(-2)
u5p(x) = x^(-5)
u6p(x) = x^(-10)

xax = range(0.05,2;length=1000)

u1.(xax)

# Polynomial Interpolation: Newton Basis
# Code from Giray Otken, First semester in numerical analysis with Julia

    function diff(x::Array,y::Array)
        m = length(x) #here m is the number of data points. #the degree of the polynomial is m-1 a=Array{Float64}(undef,m)
        a = [y[i] for i in 1:m]
        # for i in 1:m
        #     a[i]=y[i]
        # end
        #a = [(y[i]-y[i-1])/(x[i]-x[i-j+1]) for j in 2:m, i in reverse(collect(j:m))]
        for j in 2:m
            for i in reverse(collect(j:m))
                a[i]=(a[i]-a[i-1])/(x[i]-x[i-(j-1)])
            end
        end
        return a
    end

    function newton(x::Array,y::Array,z)
        m=length(x) #here m is the number of data points, not the degree # of the polynomial
        a=diff(x,y)
        sum=a[1]
        pr=1.0
        for j in 1:(m-1)
            pr=pr*(z-x[j])
            sum=sum+a[j+1]*pr
        end
        return sum
    end

    # Runge's Example with Newton Polynomial for different n

nom = ["Log","Square","CES σ=2", "CES σ=5", "CES σ=10"]
funs = [u1,u2,u4,u5,u6]
for i in 1:5
    name = nom[i]
    for n = (4,6,11,21)
        fn = funs[i].(xax)
        # Grid of nodes for interpolation
        xi = collect(range(0.05,2;length=n)) ; # Collect makes it an array instead of a collection
        yi = funs[i].(xi) # the corresponding y-coordinates
        # Interpolation
        interp=map(z->newton(xi,yi,z),xax) # Interpolating poly for the data
        # Plot
        gr()
        plot(title="Interpolation $name n=$n - Newton Polynomial")
        plot!(xax,fn,linewidth=3,label = "Function: $name",legend=(0.75,0.75),foreground_color_legend = nothing,background_color_legend = nothing)
        plot!(xax,interp,linewidth=3,label="Interpolation")
        plot!(xi,yi,linetype=:scatter,marker=(:diamond,9),markercolor=RGB(0.5,0.1,0.1),label = "Data")
        savefig("graphs/Newton $name _n_$n")

    end
end

function poly_inter_err(a,b,n,f::Function)
    a , b  = min(a,b), max(a,b)
    h = (b-a)/n
    x = [a+(i-1)*h for i in 1:(n+1)]
    Π = h^(n+1)*factorial(n)/4
    ξ = abs(maximum(f.(x)))
    err = Π*ξ/factorial(n+1)
    return err
end

err=[poly_inter_err.(0.05,2,i,funs) for i in [4,6,11]]
Err=[err[1] err[2] err[3]]

copy_to_clipboard(true)
mdtable(Err, side=nom, head=["n=4","n=6","n=11"],fmt = FancyNumberFormatter(4))|> print
mdtable(Err, side=nom, head=["n=4","n=6","n=11"],fmt = FancyNumberFormatter(4))|> display

## Cubic Splines - Natural

#The function CubicNatural takes the x and y-coordinates of the data as input, and
# computes the natural cubic spline interpolating the data, by solving the resulting matrix
# equation. The code is based on Algorithm 3.4 of Burden, Faires, Burden [4]. The output
# is the coefficients of the m 􀀀 1 cubic polynomials, ai; bi; ci; di; i = 1; :::;m 􀀀 1 where m is
# the number of data points. These coefficients are stored in the arrays a; b; c; d, which are
# declared global, so that we can access these arrays later to evaluate the spline for a given
# value w.

function CubicNatural(x::Array,y::Array)
    m=length(x) # m is the number of data points
    n=m-1
    a=Array{Float64}(undef,m)
    b=Array{Float64}(undef,n)
    c=Array{Float64}(undef,m)
    d=Array{Float64}(undef,n)
    a = y
    # for i in 1:m
    #     a[i]=y[i]
    # end
    h = [x[i+1]-x[i] for i in 1:n]
    # h=Array{Float64}(undef,n)
    # for i in 1:n
    #     h[i]=x[i+1]-x[i]
    # end
    u=Array{Float64}(undef,n)
    u[1]=0
    for i in 2:n
        u[i]=3*(a[i+1]-a[i])/h[i]-3*(a[i]-a[i-1])/h[i-1]
    end
    s=Array{Float64}(undef,m)
    z=Array{Float64}(undef,m)
    t=Array{Float64}(undef,n)
    s[1]=1
    z[1]=0
    t[1]=0
    for i in 2:n
        s[i]=2*(x[i+1]-x[i-1])-h[i-1]*t[i-1]
        t[i]=h[i]/s[i]
        z[i]=(u[i]-h[i-1]*z[i-1])/s[i]
    end
    s[m]=1
    z[m]=0
    c[m]=0
    for i in reverse(1:n)
        c[i]=z[i]-t[i]*c[i+1]
        b[i]=(a[i+1]-a[i])/h[i]-h[i]*(c[i+1]+2*c[i])/3
        d[i]=(c[i+1]-c[i])/(3*h[i])
    end
    return a, b, c, d
end

# CubicNaturalEval. The inputs are the value at which
# the spline is evaluated, w, and the x-coordinates of the data. The function first finds the
# interval [xi; xi+1]; i = 1; :::;m-1; w belongs to, and then evaluates the spline at w using the
# corresponding cubic polynomial.

function CubicNaturalEval(x::Array,y::Array,w)
    m=length(x)
    if w<x[1]||w>x[m]
        return print("error: spline evaluated outside its domain")
    end
    n=m-1
    p=1
    for i in 1:n
        if w<=x[i+1]
            break
        else
            p=p+1
        end
    end
    # p is the number of the subinterval w falls into, i.e., p=i means
    # w falls into the ith subinterval $(x_i,x_{i+1}), and therefore
    # the value of the spline at w is
    # a_i+b_i*(w-x_i)+c_i*(w-x_i)^2+d_i*(w-x_i)^3.
    a, b, c, d = CubicNatural(x,y)
    return a[p]+b[p]*(w-x[p])+c[p]*(w-x[p])^2+d[p]*(w-x[p])^3
end

nom = ["Log","Square","CES σ=2", "CES σ=5", "CES σ=10"]
funs = [u1,u2,u4,u5,u6]
for i in 1:5
    name = nom[i]
    for n = (4,6,11,21)
        fn = funs[i].(xax)
        # Grid of nodes for interpolation
        xi = collect(range(0.05,2;length=n)) ; # Collect makes it an array instead of a collection
        yi = funs[i].(xi) # the corresponding y-coordinates
        # Interpolation
        interp=map(z->CubicNaturalEval(xi,yi,z),xax) # Interpolating poly for the data
        # Plot
        gr()
        plot(title="Interpolation $name n=$n - Cubic Spline Natural")
        plot!(xax,fn,linewidth=3,label = "Function: $name",legend=(0.75,0.75),foreground_color_legend = nothing,background_color_legend = nothing)
        plot!(xax,interp,linewidth=3,label="Interpolation")
        plot!(xi,yi,linetype=:scatter,marker=(:diamond,9),markercolor=RGB(0.5,0.1,0.1),label = "Data")
        savefig("graphs/CSN_$name _n_$n")
    end
end

## Schumaker shape preserving Spline

##find_xi funtion that takes in x and y coordinates
##and slopes si at each xi and returns the ξ parameter
function find_xi(t::Array,z::Array,s::Array)
    #Based on Algorithm Tudd pag 233
    n = length(t)
    δ = [(z[i+1]-z[i])/(t[i+1]-t[i]) for i in 1:(n-1)]
    ξ = zeros(n-1)
    #Check if lemma 6.11.1 applies
    for i in 1:(n-1)
        if (s[i] + s[i+1])/2 == δ[i]
            ξ[i] = t[i+1]
        elseif (s[i]-δ[i])*(s[i+1]-δ[i]) ≥ 0
            ξ[i] = (t[i]+t[i+1])/2
        elseif abs(s[i+1]-δ[i]) < abs(s[i]-δ[i])
            ξub = t[i]+ (2*(t[i+1]-t[i])*(s[i+1]-δ[i]))/(s[i+1]-s[i])
            ξ[i] = (t[i]+ξub)/2
        else # abs(s[i+1]-δ[i]) ≥ abs(s[i]-δ[i])
            ξlb = t[i+1]+ (2*(t[i+1]-t[i])*(s[i]-δ[i]))/(s[i+1]-s[i])
            ξ[i] = (t[i+1]+ξlb)/2
        end
    end
    return ξ, δ
end

#Schumaker takes in x, y coordinates, slopes

function Schumaker(t::Array,z::Array,s::Array,w)
    m = length(t)
    if w < t[1] || w > t[m]
        return print("error: spline evaluated outside its domain")
    end
    n = m-1
    ξ, δ = find_xi(t,z,s)
    i = 1
    for j in 1:n
        if w ≤ t[j+1]
            break
        else
            i=i+1
        end
    end
    α = ξ[i]-t[i]
    β = t[i+1]-ξ[i]
    sb = (2*(z[i+1]-z[i])-(α*s[i]+β*s[i+1]))/(t[i+1]-t[i])
    A1, B1, C1 = z[i], s[i], (sb-s[i])/(2*α)
    A2, B2, C2 = A1 + α*B1 + C1*(α^2), sb, (s[i+1]-sb)/(2*β)
    if w ≤ ξ[i]
        w_t = A1 + B1*(w - t[i]) + C1*(w-t[i])^2
    else
        w_t = A2 + B2*(w - ξ[i]) + C2*(w-ξ[i])^2
    end
    return w_t, ξ, δ
end

# logfun = u1.(xax)
# # Grid of nodes for interpolation
# xi = collect(range(0.05,2;length=6)) ; # Collect makes it an array instead of a collection
# yi = u1.(xi) # the corresponding y-coordinates
# si = u1p.(xi) #the corresponding slopes
# xip = find_xi(xi,yi,si)[1]
#
# plot(xax,logfun)
# vline!(xi, color=:blue)
# vline!(xip, color=:red)
#
# # Interpolation
# interp=map(z->Schumaker(xi,yi,si,z)[1],xax) # Interpolating poly for the data
# # Plot
# gr()
# plot(xax,logfun,linewidth=3,legend=:bottomright,foreground_color_legend = nothing,background_color_legend = nothing)
# plot!(xax,interp,linewidth=3,label="Interpolation")
# plot!(xi,yi,linetype=:scatter,marker=(:diamond,9),markercolor=RGB(0.5,0.1,0.1),label = "Data")
#
# xi = [0.05,0.10]
# yi = u1.(xi)
# si = u1p.(xi)
#
# find_xi(xi,yi,si)
# Schumaker(xi,yi,si,0.075)

nom = ["Log","Square Root","CES σ=2", "CES σ=5", "CES σ=10"]
funs = [u1,u2,u4,u5,u6]
slopes = [u1p,u2p,u4p,u5p,u6p]
for i in 1:5
    name = nom[i]
    for n = (4,6,11,21)
        fn = funs[i].(xax)
        # Grid of nodes for interpolation
        global xi = collect(range(0.05,2;length=n)) ; # Collect makes it an array instead of a collection
        global yi = funs[i].(xi) # the corresponding y-coordinates
        global si = slopes[i].(xi) #the corresponding slopes
        # Interpolation
        interp=map(z->Schumaker(xi,yi,si,z)[1],xax) # Interpolating poly for the data
        # Plot
        gr()
        plot(title="Interpolation $name n=$n - Schumaker Shape Preserving Spline")
        plot!(xax,fn,linewidth=3,label = "Function: $name",legend=:bottomright,foreground_color_legend = nothing,background_color_legend = nothing)
        plot!(xax,interp,linewidth=3,label="Interpolation")
        plot!(xi,yi,linetype=:scatter,marker=(:diamond,9),markercolor=RGB(0.5,0.1,0.1),label = "Data")
        savefig("graphs/Schu_$name _n_$n")
    end
end

## graphs
##
fnom = ["Log","Square Root","CES "*latexinline("σ=2"), "CES "*latexinline("σ=5"), "CES "*latexinline("σ=10")]
funs = [u1,u2,u4,u5,u6]
slopes = [u1p,u2p,u4p,u5p,u6p]
method = [newton, CubicNaturalEval, Schumaker]
gridpts = [4,6,11,21]
nom = ["Log","SQR","CES_2", "CES_5", "CES_10"]
titlenom = ["Log","Square Root","CES σ=2", "CES σ=5", "CES σ=10"]
#titlenom = ["Log_n=6","Square_Root n=6","CES σ=2 n=6", "CES σ=5 n=6", "CES σ=10 n=6"]
for i in 1:length(funs)
    plotname = nom[i]
    title = titlenom[i]
    name = fnom[i]
    fn = funs[i].(xax)
    IE = Array{Float64}(undef, 3, 4)
    local interp = Array{Float64}(undef,1000,3)
    for n = gridpts
        for j in 1:length(method)
            met = method[j]
            # Grid of nodes for interpolation
            global xi = collect(range(0.05,2;length=n)) ; # Collect makes it an array instead of a collection
            global yi = funs[i].(xi) # the corresponding y-coordinates
            # Interpolation
            if j ≤ 2
                global interp[:,j] = map(z->method[j](xi,yi,z)[1],xax)
            else
                global si = slopes[i].(xi) #the corresponding slopes
                global interp[:,j]=map(z->Schumaker(xi,yi,si,z)[1],xax) # Interpolating poly for the data
            end
            #IE[j,n] = interp_err_SN(fn,interp)
            #IE[j,n] = sqrt(sum((fn .- interp).^2))
        end
        global fmet = ["$name","Newton Basis Polynomials","Natural Cubic Splines","Shape-preserving Schumaker Splines"]
        # Plot
        gr()
        plot(title="Interpolation $title Grid size n=$n")
        plot!(xax,fn,linewidth=3,label = "Function: $title",legend=:bottomright,foreground_color_legend = nothing,background_color_legend = nothing)
        plot!(xax,interp[:,1],linewidth=3,label=fmet[2], linestyle=:dash)
        plot!(xax,interp[:,2],linewidth=3,label=fmet[3],linestyle=:dot)
        plot!(xax,interp[:,3],linewidth=3,label=fmet[4], linestyle=:dashdot)
        plot!(xi,yi,linetype=:scatter,marker=(:circle,6),markercolor=:blue,label = "Data")
        savefig("graphs/$plotname $n.png")
    end
end

## Errors
##

fnom = ["Log","Square Root","CES "*latexinline("σ=2"), "CES "*latexinline("σ=5"), "CES "*latexinline("σ=10")]
funs = [u1,u2,u4,u5,u6]
slopes = [u1p,u2p,u4p,u5p,u6p]
method = [newton, CubicNaturalEval, Schumaker]
gridpts = [4,6,11,21]
nom = ["Log","SQR","CES_2", "CES_5", "CES_10"]
titlenom = ["Log","Square Root","CES σ=2", "CES σ=5", "CES σ=10"]
for i in 1:length(funs)
    plotname = nom[i]
    title = titlenom[i]
    name = fnom[i]
    fn = funs[i].(xax)
    IE = Array{Float64}(undef, 3, 4)
    for j in 1:length(method)
        for n in 1:length(gridpts)
            # Grid of nodes for interpolation
            global xi = collect(range(0.05,2;length=gridpts[n])) ; # Collect makes it an array instead of a collection
            global yi = funs[i].(xi) # the corresponding y-coordinates
            # Interpolation
            if j ≤ 2
                global interp = map(z->method[j](xi,yi,z)[1],xax)
            else
                global si = slopes[i].(xi) #the corresponding slopes
                global interp=map(z->Schumaker(xi,yi,si,z)[1],xax) # Interpolating poly for the data
            end
            IE[j,n] = sqrt(sum((fn .- interp).^2))
        end
    end
    global fmet = ["$name","Newton Basis Polynomials","Natural Cubic Splines","Shape-preserving Schumaker Splines"]
    mdtable(IE, side=fmet, head=["n=4","n=6","n=11","n=21"],fmt = FancyNumberFormatter(4))|> print
    # Plot
    gr()
    plot(title=title)
    plot!(gridpts,IE[1,:],linewidth=3,label = fmet[2], yscale=:log10, legend=:topright,foreground_color_legend = nothing,background_color_legend = nothing)
    plot!(gridpts,IE[2,:],linewidth=3,label = fmet[3], yscale=:log10)
    plot!(gridpts,IE[3,:],linewidth=3,label = fmet[4], yscale=:log10)
    savefig("graphs/IntError_$plotname")
end

## Grid spacing
##

# # No need for ScaledInterpolations if not using PolyRange
#grid = range(0,1,length=n)
# grid = a + (b-a)*grid^θ PolYnomial scaling
# a = 0.01, b = 0.5, θ = 2, interval = [0,1], number of points = 6
#    ci = 0.01 .+ (0.50-0.01).*range(0,1;length=6).^2

function polygrid(n::Integer,a::Float64,b::Float64,θ::Float64)
    grid = collect(range(0,1; length=n))
    xi = a .+ (b-a).*grid.^θ
    return xi
end

polygrid(6,0.05,2.0,2.0)

function interp_err_SN(F::Array,I::Array)
    error = maximum(abs.(I.-F))/maximum(abs.(F))
end

maximum(abs.(interp.-u1.(xax)))/maximum(abs.(u1.(xax)))

interp_err_SN(u1.(xax),interp)

fnom = ["Log","Square Root","CES "*latexinline("σ=2"), "CES "*latexinline("σ=5"), "CES "*latexinline("σ=10")]
funs = [u1,u2,u4,u5,u6]
slopes = [u1p,u2p,u4p,u5p,u6p]
method = [newton, CubicNaturalEval, Schumaker]
curvature = [1.0,2.0,3.0,4.0]
nom = ["Log","SQR","CES_2", "CES_5", "CES_10"]
titlenom = ["Log","Square Root","CES σ=2", "CES σ=5", "CES σ=10"]
#titlenom = ["Log_n=6","Square_Root n=6","CES σ=2 n=6", "CES σ=5 n=6", "CES σ=10 n=6"]
for i in 1:length(funs)
    plotname = nom[i]
    title = titlenom[i]
    name = fnom[i]
    fn = funs[i].(xax)
    IE = Array{Float64}(undef, 3, 4)
    for j in 1:length(method)
        met = method[j]
        for n in 1:length(curvature)
            curv = curvature[n]
            # Grid of nodes for interpolation
            global xi = polygrid(6,0.05,2.0,curvature[n]) ; # Collect makes it an array instead of a collection
            global yi = funs[i].(xi) # the corresponding y-coordinates
            # Interpolation
            if j ≤ 2
                global interp = map(z->method[j](xi,yi,z)[1],xax)
            else
                global si = slopes[i].(xi) #the corresponding slopes
                global interp=map(z->Schumaker(xi,yi,si,z)[1],xax) # Interpolating poly for the data
            end
            #IE[j,n] = interp_err_SN(fn,interp)
            IE[j,n] = sqrt(sum((fn .- interp).^2))
        end
    end
    global fmet = ["$name","Newton Basis Polynomials","Natural Cubic Splines","Shape-preserving Schumaker Splines"]
    mdtable(IE, side=fmet, head=latexinline(["θ=1","θ=2","θ=3","θ=4"]),fmt = FancyNumberFormatter(4))|> print
    # Plot
    gr()
    #pyplot()
    plot(title=title, xticks=curvature)
    plot!(curvature,IE[1,:],linewidth=3,label = fmet[2], yscale=:log10,  legend=:topleft,foreground_color_legend = nothing,background_color_legend = nothing)
    plot!(curvature,IE[2,:],linewidth=3,label = fmet[3],yscale=:log10,)
    plot!(curvature,IE[3,:],linewidth=3,label = fmet[4],yscale=:log10,)
    xlabel!("Curvature")
    savefig("graphs/IntError_Curv_$plotname")
end

fnom = ["Log","Square Root","CES "*latexinline("σ=2"), "CES "*latexinline("σ=5"), "CES "*latexinline("σ=10")]
funs = [u1,u2,u4,u5,u6]
slopes = [u1p,u2p,u4p,u5p,u6p]
method = [newton, CubicNaturalEval, Schumaker]
curvature = [1.0,1.5,2.0,3.0]
nom = ["Log","SQR","CES_2", "CES_5", "CES_10"]
titlenom = ["Log","Square Root","CES σ=2", "CES σ=5", "CES σ=10"]
#titlenom = ["Log_n=6","Square_Root n=6","CES σ=2 n=6", "CES σ=5 n=6", "CES σ=10 n=6"]
for i in 1:length(funs)
    plotname = nom[i]
    title = titlenom[i]
    name = fnom[i]
    fn = funs[i].(xax)
    IE = Array{Float64}(undef, 3, 4)
    local interp = Array{Float64}(undef,1000,3)
    for n in 1:length(curvature)
        curv = curvature[n]
        for j in 1:length(method)
            met = method[j]
            # Grid of nodes for interpolation
            global xi = polygrid(6,0.05,2.0,curvature[n]) ; # Collect makes it an array instead of a collection
            global yi = funs[i].(xi) # the corresponding y-coordinates
            # Interpolation
            if j ≤ 2
                global interp[:,j] = map(z->method[j](xi,yi,z)[1],xax)
            else
                global si = slopes[i].(xi) #the corresponding slopes
                global interp[:,j]=map(z->Schumaker(xi,yi,si,z)[1],xax) # Interpolating poly for the data
            end
            #IE[j,n] = interp_err_SN(fn,interp)
            #IE[j,n] = sqrt(sum((fn .- interp).^2))
        end
        global fmet = ["$name","Newton Basis Polynomials","Natural Cubic Splines","Shape-preserving Schumaker Splines"]
        # Plot
        gr()
        plot(title="Interpolation n=6 - $title, θ=$curv")
        plot!(log.(xax),fn,linewidth=3,label = "Function: $title",legend=:bottomright,foreground_color_legend = nothing,background_color_legend = nothing)
        plot!(log.(xax),interp[:,1],linewidth=3,label=fmet[2], linestyle=:dash)
        plot!(log.(xax),interp[:,2],linewidth=3,label=fmet[3],linestyle=:dot)
        plot!(log.(xax),interp[:,3],linewidth=3,label=fmet[4], linestyle=:dashdot)
        plot!(log.(xi),yi,linetype=:scatter,marker=(:circle,6),markercolor=:blue,label = "Data")
        savefig("graphs/$plotname $curv.png")
    end

    # mdtable(IE, side=fmet, head=latexinline(["θ=1","θ=2","θ=3","θ=4"]),fmt = FancyNumberFormatter(4))|> print
    # # Plot
    # gr()
    # #pyplot()
    # plot(title=title, xticks=curvature)
    # plot!(curvature,IE[1,:],linewidth=3,label = fmet[2], yscale=:log10,  legend=:topleft,foreground_color_legend = nothing,background_color_legend = nothing)
    # plot!(curvature,IE[2,:],linewidth=3,label = fmet[3],yscale=:log10,)
    # plot!(curvature,IE[3,:],linewidth=3,label = fmet[4],yscale=:log10,)
    # xlabel!("Curvature")
    # savefig("graphs/IntError_Curv_$plotname")
end

# function curv(n::Float64,a::Float64,b::Float64,F::function,T::function)
#     error = 1.0
#     θ = 1
#     xi = polygrid(n,a,b,θ)
#     yi = F.(xi)
#     fn = F.(xax)
#     interp = map(z->T(xi,yi,z),xax)
#     while error > 10E-5

max_stars_repo_name: py-r-hans/AdvMacro

max_issues_repo_name: py-r-hans/AdvMacro

max_forks_repo_name: py-r-hans/AdvMacro


text:
include(joinpath(@__DIR__, "../../src/FLS.jl"))

"""
CONSTANTS

Models are two constants from a set. A constant is an arrow from 1, the set with one element.
"""

s1 = Symbol("1")

const2schema = @acset LabeledGraph begin
    V = 2
    E = 2
    vlabel = [s1, :A]
    elabel = [:f, :g]
    src = [1,1]
    tgt = [2,2]
end

const2 = FLS(:const2, const2schema, [Cone(c2d, s1,[])], [])


text:
# Define variable
x = 1

# Define a function
"""
Add one to input (x)
"""
function addone(x)
  return x + 1
end

# Arrays
arr = [1, 2, 3]

# For comprehensions
arr_one_to_ten = 1:10
arr_one_to_ten = [i * 2 for i in 1:10]
one_to_ten_by_2 = 1:2:10

# For loop
for i in 1:10
  # Prints a newline after
  println(i)

  # Just prints i
  # print(i)
end

# Run command in terminal
run(`ls`)

"""
Add two to an integer
"""
function addtwo(x::Integer)
  return x + 2
end

# Get type of a variable
typeof('s')  # CHar
typeof("s")  # String

# Make your type
# Immutable is the default (without keyword mutable)
mutable struct Person
  name::String
  birth_year::Integer
  # age::Integer
end

person = Person("Arthur", 1989)

# person2 = Person(30, "Arthur")  # This throws error

# function increment_age(person::Person)
#   person.age += 1
# end

# increment_age(person)

function get_age(person, current_year::Integer)
  return current_year - person.birth_year
end

get_age(person, 2019)

# Generate random matrix
Z = randn(100, 5)

# element-wise mult
Z .* Z

Z'Z

# Greek letters
μ = 1

# Use a package
using Distributions

f(x) = x + 1

add_f(x, f) = f(x) + x

# Multiple dispatch (requires type annotation)
say_hi(x::String) = println("Hi, $(x)")
say_hi(x::Int) = println("Hi, $(x + 1)")

function readfile(fpath)
  x = open(fpath, "r") do f
    readlines(f)
  end
  return x
end

contents = readfile("names.txt")

using RCall
# Intro to RCall:
# http://luiarthur.github.io/usingrcall


text:
# https://www.hackerrank.com/challenges/30-data-types/problem

# Declare second integer, double, and String variables.
# Read and save an integer, double, and String to your variables.
inputs = readlines()
i = i + parse(Int, inputs[1])
d = d + parse(AbstractFloat, inputs[2])
s = s * " " * inputs[3]

# Print the sum of both integer variables on a new line.
println(i)
# Print the sum of the double variables on a new line.
println(d)
# Concatenate and print the String variables on a new line
# The 's' variable above should be printed first.
println(s)

max_stars_repo_path: Practice/Tutorials/30 Days of Code/Day 01 Data Types.jl

max_stars_repo_name: mtaufanr/hackerrank-code-julia

max_issues_repo_name: mtaufanr/hackerrank-code-julia

max_forks_repo_name: mtaufanr/hackerrank-code-julia


text:
export ISelectionStrategy, NthOfTypeStrategy
export IGroupAssemblyStrategy, NthGroupStrategy

abstract type ISelectionStrategy end

function selectFrom(this::ISelectionStrategy, phys::Vector{vk.VkPhysicalDevice})::vk.VkPhysicalDevice
    return vk.VK_NULL_HANDLE
end

mutable struct NthOfTypeStrategy <: ISelectionStrategy
    mType::vk.VkPhysicalDeviceType
    mN::UInt32

    NthOfTypeStrategy(type::vk.VkPhysicalDeviceType) = new(type, 0)
end
@class NthOfTypeStrategy

function selectFrom(this::NthOfTypeStrategy, phys::Vector{vk.VkPhysicalDevice})::vk.VkPhysicalDevice
    counter::UInt32 = 0
    for dev in phys
        props = VkExt.getProperties(dev)
        if (props.deviceType == this.mType)
            if (counter == this.mN)
                return dev;
            end
            counter += 1
        end
    end
    println("Failed to select physical device of ", this.mType, ".")
    return vk.VK_NULL_HANDLE
end

abstract type IGroupAssemblyStrategy end

function assembleFrom(this::IGroupAssemblyStrategy, groups::Vector{vk.VkPhysicalDeviceGroupProperties})::Vector{vk.VkPhysicalDevice}
    return []
end

mutable struct NthGroupStrategy <: IGroupAssemblyStrategy
    mN::UInt32

    # n start from 0
    NthGroupStrategy(n::Integer) = new(n)
end
@class NthGroupStrategy

function assembleFrom(this::NthGroupStrategy, groups::Vector{vk.VkPhysicalDeviceGroupProperties})::Vector{vk.VkPhysicalDevice}
    @assert(length(groups) > this.mN)

    group = groups[this.mN + 1]
    ret = []
    for i = 1 : group.physicalDeviceCount
        push!(ret, group.physicalDevices[i])
    end
    return ret
end


text:
using VITA2D
using Test
using Random
Random.seed!(42)

@testset "VITA2D.jl" begin
    @test half(4) == 2
    @test third(30) == 10

    L = 6
    P = 2
    θ = rand(3P)
    h = (ZZ = 1.0,
         X = 1.0,
         Z = 1.0)

    M = Model(L, θ, h)
    Obs = Thermo(P)
    Grad = Gradient(Obs, P)
    Sto = StochasticMatrix(Obs, P)
    C = Config((L, L, 2P+1))

    warmup!(C, M, 1)
    mc(C, M, Obs, Grad, 1, Sto)

    sgd_params = VITA2D.SGDParams()

    E = sgd(L=L,
              P=P,
              hX = 1.0,
              θ = θ,
              SGD=SGDParams(),
              path=false)

    E = sgd(L=L,
            P=P,
            hX = 3.0,
            θ = θ,
            SGD=sgd_params,
            print=true)
end


text:
using OrdinaryDiffEq, Plots, LaTeXStrings
using OpenQuantumTools

# Hamiltonian
H = DenseHamiltonian([(s)->1-s, (s)->s], -[σx, σz]/2, unit=:ħ)
# initial state
u0 = PauliVec[1][1]
# coupling
coupling = ConstantCouplings(["Z"], unit=:ħ)
# bath
bath = Ohmic(1e-4, 4, 16)
annealing = Annealing(H, u0; coupling=coupling, bath=bath)

tf = 60
U = solve_unitary(annealing, tf, alg=Tsit5(), abstol=1e-8, reltol=1e-8)
U = InplaceUnitary(U)

@time solr = solve_redfield(annealing, tf, U, alg=Tsit5())
# we set the integration error tolerance to 1e-5 for speed
@time solc = solve_cgme(annealing, tf, U, alg=Tsit5(), int_atol=1e-5, int_rtol=1e-5)
plot(solr, H, [0], 0:0.01:tf, linewidth=2, xlabel="t (ns)", ylabel="\$P_G(t)\$", label="Redfield")
plot!(solc, H, [0], 0:0.01:tf, linewidth=2, label="CGME")

using QuadGK

g(t) = quadgk((w)->sqrt(γ(w, bath))*exp(1.0im*w*t)/2/π, -Inf, Inf)[1]

t = range(-0.5,0.5,length=500)
g_value = g.(t)
c_value = [correlation(x, bath) for x in t];

plot(t, real.(g_value), label="Re[g(t)]", linewidth=2)
plot!(t, imag.(g_value), label="Im[g(t)]", linewidth=2)
plot!(t, real.(c_value), label="Re[C(t)]", linewidth=2)
plot!(t, imag.(c_value), label="Im[C(t)]", linewidth=2)
xlabel!("t (ns)")
ylabel!("correlation")

t = range(-4,4,length=2000)
g_value = g.(t)
gf = construct_interpolations(t, g_value, extrapolation = "flat")

t = range(-0.5,0.5,length=500)
g_value = gf.(t)
plot(t, real.(g_value), label="Re[g(t)]", linewidth=2)
plot!(t, imag.(g_value), label="Im[g(t)]", linewidth=2)
xlabel!("t (ns)")
ylabel!("fitted jump correlator")

ubath = ULEBath(gf)
annealing = Annealing(H, u0; coupling=coupling, bath=ubath)
@time solu = solve_ule(annealing, tf, U, alg=Tsit5(), int_atol=1e-5, int_rtol=1e-5)
plot(solr, H, [0], 0:0.01:tf, linewidth=2, xlabel="t (ns)", ylabel="\$P_G(t)\$", label="Redfield")
plot!(solc, H, [0], 0:0.01:tf, linewidth=2, label="CGME")
plot!(solu, H, [0], 0:0.01:tf, linewidth=2, label="ULE")


text:
using OffsetArrays

function reg_term(omega, B = size(omega, 2); alpha = 0.01)
  reg = 0.0
  for b = 1:B
    t = maximum(@views omega[:, b]) - omega[b, b]
    reg += t
  end
  return alpha * reg
end

function reg_term_turbo(omega, B = size(omega, 2); alpha = 0.01)
  reg = 0.0
  @turbo for b = 1:B
    m = -Inf
    for i in axes(omega, 1)
      m = max(m, omega[i, b])
    end
    t = m - omega[b, b]
    reg += t
  end
  return alpha * reg
end

function awmean_lv(x::AbstractArray{T1}, σ::AbstractArray{T2}) where {T1<:Number,T2<:Number}
  n = length(x)
  T3 = promote_type(T1, T2)
  T = sizeof(T3) ≤ 4 ? Float32 : Float64
  sum_of_values = sum_of_weights = χ2 = zero(T)
  @turbo for i = 1:n
    sum_of_values += x[i] / (σ[i] * σ[i])
    sum_of_weights += 1 / (σ[i] * σ[i])
  end
  wx = sum_of_values / sum_of_weights

  @turbo for i = 1:n
    χ2 += (x[i] - wx) * (x[i] - wx) / (σ[i] * σ[i])
  end
  mswd = χ2 / (n - 1)
  wσ = sqrt(one(T) / sum_of_weights)
  return wx, wσ, mswd
end
function awmean_simd(
  x::AbstractArray{T1},
  σ::AbstractArray{T2},
) where {T1<:Number,T2<:Number}
  n = length(x)
  T3 = promote_type(T1, T2)
  T = sizeof(T3) ≤ 4 ? Float32 : Float64
  sum_of_values = sum_of_weights = χ2 = zero(T)
  @fastmath @inbounds @simd for i = 1:n
    sum_of_values += x[i] / (σ[i] * σ[i])
    sum_of_weights += 1 / (σ[i] * σ[i])
  end
  wx = sum_of_values / sum_of_weights

  @fastmath @inbounds @simd for i = 1:n
    χ2 += (x[i] - wx) * (x[i] - wx) / (σ[i] * σ[i])
  end
  mswd = χ2 / (n - 1)
  wσ = sqrt(one(T) / sum_of_weights)
  return wx, wσ, mswd
end

function test_awmean(::Type{T}) where {T}
  for n ∈ 2:100
    if T <: Integer
      x = view(rand(T(-50):T(100), n + 32), 17:n+16)
      σ = view(rand(T(1):T(10), n + 32), 17:n+16)
    else
      x = view((randn(T, n + 32) .+= T(2)), 17:n+16)
      σ = view(rand(T, n + 32), 17:n+16)
    end
    wx, wσ, mswd = awmean_simd(x, σ)
    @test iszero(@allocated((wxlv, wσlv, mswdlv) = awmean_lv(x, σ)))
    wxlv, wσlv, mswdlv = awmean_lv(x, σ)
    isfinite(wx) && @test wx ≈ wxlv
    isfinite(wσ) && @test wσ ≈ wσlv
    isfinite(mswd) &&
      @test mswd ≈ mswdlv atol = eps(typeof(mswdlv)) rtol = sqrt(eps(typeof(mswdlv)))
  end
end

function logℒ_fast(α, β, t, c, x)
  eα = abs(α)
  n, k = size(x)

  (n == length(t) == length(c) && length(β) == k + 1) || throw(DimensionMismatch())
  s = zero(typeof(α))
  @inbounds for i = 1:n
    xb = sum(@inbounds(x[i, j] * β[j+1]) for j = 1:k) + β[1]
    ti = t[i]
    s += (1 - (c[i] == ti)) * (log(eα) + (eα - 1) * log(ti) + xb) - ti^eα * exp(xb)
  end
  return s
end
function logℒ_fast_turbo(α, β, t, c, x)
  eα = abs(α)
  n, k = size(x)

  (n == length(t) == length(c) && length(β) == k + 1) || throw(DimensionMismatch())
  s = zero(typeof(α))
  @turbo for i = 1:n
    xb0 = 0.0
    for j = 1:k
      xb0 += x[i, j] * β[j+1]
    end
    xb = xb0 + β[1]
    ti = t[i]
    s += (1 - (c[i] == ti)) * (log(eα) + (eα - 1) * log(ti) + xb) - ti^eα * exp(xb)
  end
  return s
end

function test_logℒ(n, k)
  t = rand(n)
  c = copy(t)
  b = rand(n) .> 0.5
  c[b] .= rand.()
  x = rand(n, k)
  α = 2.85
  β = rand(k + 1)
  @test logℒ_fast(α, β, t, c, x) ≈ logℒ_fast_turbo(α, β, t, c, x)
end

function not_an_outer_reduct!(r, N::Int, x = 2.0, y = nothing) # there was a bug where this was classified as one
  @turbo for i ∈ eachindex(r)
    acc = y === nothing ? x : r[i]
    for n ∈ 1:N
      acc += 0
    end
    r[i] = acc
  end
  r
end
function tk(v::Vector{T}, c_::Vector{T}, n::Int) where {T}
  l00 = l10 = l01 = l11 = zero(T)
  @turbo for j ∈ 1:n-1
    a00 = v[3n-3j+1]
    a01 = -v[3n-3j+2]
    a10 = v[3n-3j+2]
    a11 = v[3n-3j+3]
    b00 = c_[4j+1]
    b01 = c_[4j+2]
    b10 = c_[4j+3]
    b11 = c_[4j+4]
    l00 += a00 * b00 + a01 * b10
    l01 += a00 * b01 + a01 * b11
    l10 += a10 * b00 + a11 * b10
    l11 += a10 * b01 + a11 * b11
  end
  l00 + l01 + l10 + l11
end
function tk2(v::Vector{T}, c_::Vector{T}, n::Int) where {T}
  l00 = l10 = l01 = l11 = zero(T)
  @turbo for j ∈ 1:n-1
    i = n - j
    a00 = v[3i+1]
    a01 = -v[3i+2]
    a10 = v[3i+2]
    a11 = v[3i+3]
    b00 = c_[4j+1]
    b01 = c_[4j+2]
    b10 = c_[4j+3]
    b11 = c_[4j+4]
    l00 += a00 * b00 + a01 * b10
    l01 += a00 * b01 + a01 * b11
    l10 += a10 * b00 + a11 * b10
    l11 += a10 * b01 + a11 * b11
  end
  l00 + l01 + l10 + l11
end
function tk_base(v::Vector{T}, c_::Vector{T}, n::Int) where {T}
  l00 = l10 = l01 = l11 = zero(T)
  @inbounds @simd for j ∈ 1:n-1
    a00 = v[3n-3j+1]
    a01 = -v[3n-3j+2]
    a10 = v[3n-3j+2]
    a11 = v[3n-3j+3]
    i = n - j
    a00 = v[3i+1]
    a01 = -v[3i+2]
    a10 = v[3i+2]
    a11 = v[3i+3]
    b00 = c_[4j+1]
    b01 = c_[4j+2]
    b10 = c_[4j+3]
    b11 = c_[4j+4]
    l00 += a00 * b00 + a01 * b10
    l01 += a00 * b01 + a01 * b11
    l10 += a10 * b00 + a11 * b10
    l11 += a10 * b01 + a11 * b11
  end
  l00 + l01 + l10 + l11
end

function noturbosum(x, kern)
  ks = zero(eltype(x))
  @inbounds @fastmath for i = 2:size(x, 1)-1
    for j in axes(kern, 1), m in axes(kern, 2)
      ks += x[i+j, i+m] * kern[j, m]
    end
  end
  ks
end
function outersum(x, kern)
  ks = zero(eltype(x))
  @turbo for i = 2:size(x, 1)-1
    for j in axes(kern, 1), m in axes(kern, 2)
      ks += x[i+j, i+m] * kern[j, m]
    end
  end
  ks
end
function insum(x, kern)
  ks = zero(eltype(x))
  @turbo for i = 2:size(x, 1)-1
    s1 = zero(eltype(x))
    for j in axes(kern, 1), m in axes(kern, 2)
      s1 += x[i+j, i+m] * kern[j, m]
    end
    ks += s1
  end
  ks
end
function in2sum(x, kern)
  ks = zero(eltype(x))
  @turbo for i = 2:size(x, 1)-1
    for j in axes(kern, 1)
      s0 = zero(eltype(x))
      for m in axes(kern, 2)
        s0 += x[i+j, i+m] * kern[j, m]
      end
      ks += s0
    end
  end
  ks
end

@testset "Outer Reductions" begin
  for T ∈ [Float32, Float64, Int32, Int64]
    test_awmean(T)
    v = rand(T, 4 * 127)
    c = rand(T, 4 * 127)
    @test tk(v, c, 127) ≈ tk2(v, c, 127) ≈ tk_base(v, c, 127)
  end
  @test all(==(7.4), not_an_outer_reduct!(Vector{Float64}(undef, 5), 17, 7.4))
  for n ∈ 1:20, k ∈ 1:5
    test_logℒ(n, k)
  end
  omega = rand(87, 87)
  @test reg_term(omega) ≈ reg_term_turbo(omega)

  x = rand(100, 100)
  kern = OffsetArray(rand(3, 3), -2, -2)
  s = noturbosum(x, kern)
  @test s ≈ outersum(x, kern)
  @test s ≈ insum(x, kern)
  @test s ≈ in2sum(x, kern)
end


text:
using DashBootstrapComponents

offcanvas = html_div([
    dbc_button("Open Offcanvas", id = "open-offcanvas", n_clicks = 0),
    dbc_offcanvas(
        html_p(
            "This is the content of the Offcanvas. " *
            "Close it by clicking on the close button, or " *
            "the backdrop.",
        ),
        id = "offcanvas",
        title = "Title",
        is_open = false,
    ),
]);

callback!(
    app,
    Output("offcanvas", "is_open"),
    Input("open-offcanvas", "n_clicks"),
    State("offcanvas", "is_open"),
) do n1, is_open
    return n1 > 0 ? is_open == 0 : is_open
end;

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
## MNIST classification using Kernel machines
using MLDatasets
train_x, train_y = MNIST.traindata()
test_x,  test_y  = MNIST.testdata()
## There are 60,000 training images and 10,000 test images
# we are going to get only a subset of them:
# 5,000 training images and 500 test images
using Random
Random.seed!(0)
subset_train = randperm(length(train_y))[1:5000]
subset_test = randperm(length(test_y))[1:500]
train_x = train_x[:,:,subset_train]
test_x = test_x[:,:,subset_test]
train_y = train_y[subset_train]
test_y = test_y[subset_test]

## see how to use a Kernel machine, which involve solving a sym. pos. def linear system
target = 9
y = 2.0.*(train_y .== target) .- 1.0 # get a -1, 1 vector for target / not target

## Work through an example with Kernel machines
# Goal, train a 9 or not classifier based on a RBF Kernel
using LinearAlgebra
function build_kernel(train_x, sigma)
  N = size(train_x)[end]
  K = zeros(N,N)
  for i=1:N
    for j=1:N
      x = @view train_x[:,:,i]
      y = @view train_x[:,:,j]
      K[i,j] = exp(-norm(x-y)^2/(2*sigma))
    end
  end
  return K
end
sigma = 1.0
K = build_kernel(train_x, 1.0)
##
a = K\y
## Compute the prediction
function predict(x,a,sigma,train_x)
  N = size(train_x)[end]
  k = zeros(N)
  for i=1:N
    k[i] = exp(-norm(x-train_x[:,:,i])^2/(2*sigma))
  end
  return k'*a
end

anim=  @animate for i=1:50
  yhat = predict(test_x[:,:,i],a,sigma,train_x)
  ytrue = test_y[i]
  color = :Greys_r
  if (ytrue == target) != (yhat .>= 0)
    color=(:Reds_r)
  end
  heatmap(Float64.(test_x[:,:,i]'),yflip=true,framestyle=:none,
    colorbar=false,color=color)
  plot!(size=(300,300),dpi=300)
  title!("Predicted $(target) = $(yhat .>= 0)")
end
gif(anim, "mnist-kernel.gif", fps=1)
## Give the number of errors
ntest = size(test_x)[end]
yhat = zeros(Bool, ntest)
for i=1:ntest
  yhat[i] = predict(test_x[:,:,i],a,sigma,train_x) .>= 0
end
## Show the confusion matrix
tp = sum( (yhat .== true) .& (test_y .== target) )
tn = sum( (yhat .== false) .& (test_y .!= target) )
fp = sum( (yhat .== true) .& (test_y .!= target) )
fn = sum( (yhat .== false) .& (test_y .== target) )
C = [tp tn ; fp fn]
accuracy = (tp+tn)/(sum(C))
## Look at accuracy across values of sigma
function accuracy_of_sigma(target, sigma)
  K = build_kernel(train_x, sigma)
  y = 2.0.*(train_y .== target) .- 1.0 # get a -1, 1 vector for target / not target
  a = K\y

  ntest = size(test_x)[end]
  yhat = zeros(Bool, ntest)
  for i=1:ntest
    yhat[i] = predict(test_x[:,:,i],a,sigma,train_x) .>= 0
  end

  tp = sum( (yhat .== true) .& (test_y .== target) )
  tn = sum( (yhat .== false) .& (test_y .!= target) )
  fp = sum( (yhat .== true) .& (test_y .!= target) )
  fn = sum( (yhat .== false) .& (test_y .== target) )
  C = [tp tn ; fp fn]
  accuracy = (tp+tn)/(sum(C))
  return accuracy, C
end
@time accuracy_of_sigma(9, 1.0)
##
@time accuracy_of_sigma(9, 1.2)
##
@time accuracy_of_sigma(9, 2.0)


text:
# Type Operators

#------------------------------------------------------------------------------
# isa operator
#= REPL
julia> 1 isa Int
true

julia> 1 isa Float64
false

julia> 1 isa Real
true
=#

#------------------------------------------------------------------------------
# What's Int type's heritage?
#= REPL
julia> supertype(Int)
Signed

julia> supertype(Signed)
Integer

julia> supertype(Integer)
Real
=#

#------------------------------------------------------------------------------
# check if Int is a subtype of Real
#= REPL
julia> Int <: Real
true
=#

#------------------------------------------------------------------------------
# documentation
#= REPL
help?> isa
search: isa isascii isapprox isabspath isassigned isabstracttype disable_sigint isnan ispath isvalid ismarked istaskdone

  isa(x, type) -> Bool

  Determine whether x is of the given type. Can also be used as an infix operator, e.g. x isa type.

help?> <:
search: <:

  <:(T1, T2)

  Subtype operator: returns true if and only if all values of type T1 are also of type T2.
=#


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

using Markdown
using InteractiveUtils

# ╔═╡ 8a4bdc28-f9b7-11ea-068e-eb70e542a539
using Pkg, DrWatson

# ╔═╡ 8a4c20a2-f9b7-11ea-2d67-a1e2531dcbd6
begin
	@quickactivate "StatisticalRethinkingTuring"
	using Turing
	using StatisticalRethinking
	using StatsPlots
end

# ╔═╡ 2a9eafae-1246-11eb-0073-f7df8df64915
using XGBoost

# ╔═╡ 4eb6e25c-f9b7-11ea-3e43-05d4d6f124f5
md"## m0.2t.jl"

# ╔═╡ 9de23082-04c7-11eb-3b7c-393a11a68cf4
md"##### This notebook takes a look at the [rstar()](https://arxiv.org/pdf/2003.07900.pdf) diagnostic. Currently I can't get MLJ & MLJModels to work (serialization issue?)."

# ╔═╡ 8a4c985e-f9b7-11ea-22cb-717b854951d6
md"##### Define a simple Normal model with unknown mean and variance."

# ╔═╡ 8a53fb6a-f9b7-11ea-1a81-5b84e8ec7e71
@model gdemo(x, y) = begin
  s ~ InverseGamma(2, 3)
  m ~ Normal(0, sqrt(s))
  x ~ Normal(m, sqrt(s))
  y ~ Normal(m, sqrt(s))
end

# ╔═╡ 8a5473b8-f9b7-11ea-3b00-75cd16c17551
md"#####  Run sampler, collect results."

# ╔═╡ 8a5e1b36-f9b7-11ea-36ef-ebe92bc31e40
chns = mapreduce(c -> sample(gdemo(1.5, 2), NUTS(0.65), 2000), chainscat, 1:4)

# ╔═╡ 8a5ea05e-f9b7-11ea-2cbe-e15b867548b3
md"##### Plot the chains."

# ╔═╡ 8a6d9e26-f9b7-11ea-1d0f-f95263880beb
plot(chns; seriestype=:traceplot)

# ╔═╡ 45425d72-f9b8-11ea-1faa-41e93626b8b5
plot(chns; seriestype=:density)

# ╔═╡ 5001e4fc-1246-11eb-38f0-a5435c180d1a
"""
    rstar(chains::Chains; subset = 0.8, niter = 1_000, eta = 0.5, XGBoostParams)
    rstar(chains::Chains, iterations::Int; subset = 0.8, niter = 1_000, eta = 0.5, XGBoostParams)
    rstar(x::AbstractMatrix, y::AbstractVector, nchains::Int, iterations::Int; subset = 0.8, niter = 1_000, eta = 0.5, XGBoostParams)

Compute the R* statistic for convergence diagnostic of MCMC. This implementation is an adaption of Algorithm 1 & 2, described in [Lambert & Vehtari]. Note that the correctness of the statistic depends on the convergence of the classifier used internally in the statistic. You can track if the training of the classifier converged by inspection of the printed RMSE values from the XGBoost backend. To adjust the number of iterations used to train the classifier set `niter` accordingly.

# Usage

```julia-repl
using XGBoost
# You need to load XGBoost before using MCMCChains.Rstar

...

chn = ...

# Compute R⋆ using defaults settings for the  gradient boosting classifier used to compute the statistic.
# This is the recomended use.
R = rstar(chn)

# Compute 100 samples of the R⋆ statistic using sampling from according to the prediction probabilities.
# This approach can be slow and results in a less accurate estimation of the R* statistic.
# See discussion in Section 3.1.3 in the paper.
Rs = rstar(chn, 100)

# estimate Rstar
R = mean(Rs)

# visualize distribution
histogram(Rs)
```

## References:
[Lambert & Vehtari] Ben Lambert and Aki Vehtari. "R∗: A robust MCMC convergence diagnostic with uncertainty using gradient-boostined machines." Arxiv 2020.
"""
function rstar(x::AbstractMatrix, y::AbstractVector, nchains::Int, iterations::Int; subset = 0.8, niter = 1_000, eta = 0.5, xgboostparams...)

    N = length(y)

    # randomly sub-select training and testing set
    Ntrain = round(Int, N*subset)
    Ntest = N - Ntrain

    ids = Random.randperm(N)
    train_ids = view(ids, 1:Ntrain)
    test_ids = view(ids, (Ntrain+1):N)

    @assert length(test_ids) == Ntest

    # use predicted probabilities?
    mode = iterations > 1 ? "multi:softprob" : "multi:softmax"

    @info "Training classifier"
    # train classifier using XGBoost
    classif = XGBoost.xgboost(x[train_ids,:], niter; label = y[train_ids],
                      objective = mode, num_class = nchains,
                      xgboostparams...)

    @info "Computing R* statistics"
    Rstats = ones(iterations) * Inf
    for i in 1:iterations
        # predict labels for "test" data
        p = XGBoost.predict(classif, x[test_ids,:])

        pred = if length(p) == Ntest*nchains
            probs = reshape(p, Ntest, nchains)
            map(s -> rand(Categorical(s / sum(s))), eachrow(probs))
        else
            p
        end

        # compute statistic
        a = mean(pred .== y[test_ids])
        Rstats[i] = nchains*a
    end

    return Rstats
end

# ╔═╡ 5002183c-1246-11eb-2867-adbcc9adab3c
function rstar(chn::Chains, iterations; kwargs...)
    nchains = size(chn, 3)
    @assert nchains > 1

    # collect data
    x = mapreduce(c -> Array(chn[:,:,c]), vcat, chains(chn))
    y = mapreduce(c -> ones(Int, length(chn))*c, vcat, chains(chn)) .- 1

    return rstar(x, y, nchains, iterations, kwargs...)
end

# ╔═╡ 5002de04-1246-11eb-0cfb-e931786d8c7f
rstar(chn::Chains; kwargs...) = first(rstar(chn, 1; kwargs...))

# ╔═╡ 4d6193a0-1246-11eb-2b25-1d555f611e5c
rstar(chns, 1)

# ╔═╡ 727fa3e8-1246-11eb-3057-6fe955204f01
rs = rstar(chns, 100)

# ╔═╡ 7280a8e2-1246-11eb-3cd5-8b6f9981eada
begin
	density(rs, lab="rs density")
	vline!([mean(rs)], lab="mean(rs)")
end

# ╔═╡ 8a6e37c6-f9b7-11ea-2abd-5db63b061a7d
md"## End m0.2t.jl"

# ╔═╡ Cell order:
# ╟─4eb6e25c-f9b7-11ea-3e43-05d4d6f124f5
# ╟─9de23082-04c7-11eb-3b7c-393a11a68cf4
# ╠═8a4bdc28-f9b7-11ea-068e-eb70e542a539
# ╠═8a4c20a2-f9b7-11ea-2d67-a1e2531dcbd6
# ╟─8a4c985e-f9b7-11ea-22cb-717b854951d6
# ╠═8a53fb6a-f9b7-11ea-1a81-5b84e8ec7e71
# ╟─8a5473b8-f9b7-11ea-3b00-75cd16c17551
# ╠═8a5e1b36-f9b7-11ea-36ef-ebe92bc31e40
# ╟─8a5ea05e-f9b7-11ea-2cbe-e15b867548b3
# ╠═8a6d9e26-f9b7-11ea-1d0f-f95263880beb
# ╠═45425d72-f9b8-11ea-1faa-41e93626b8b5
# ╠═2a9eafae-1246-11eb-0073-f7df8df64915
# ╠═5001e4fc-1246-11eb-38f0-a5435c180d1a
# ╠═5002183c-1246-11eb-2867-adbcc9adab3c
# ╠═5002de04-1246-11eb-0cfb-e931786d8c7f
# ╠═4d6193a0-1246-11eb-2b25-1d555f611e5c
# ╠═727fa3e8-1246-11eb-3057-6fe955204f01
# ╠═7280a8e2-1246-11eb-3cd5-8b6f9981eada
# ╟─8a6e37c6-f9b7-11ea-2abd-5db63b061a7d


text:
abstract PlotItem

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

type Series <: PlotItem
  axis::Symbol
  label::AbstractString
  idx::Int
  color
  plt
  x::Vector{Float64}
  y::Vector{Float64}

  function Series(axis::Symbol, label::AbstractString, idx::Int, color, plt)
    @assert axis in (:left, :right)
    new(axis, label, idx, color, plt, zeros(0), zeros(0))
  end
end

Base.string(line::Series) = "Series{axis=$(line.axis) label=$(line.label) idx=$(line.idx) npoints=$(length(line.x))}"
Base.print(io::IO, line::Series) = print(io, string(line))
Base.show(io::IO, line::Series) = print(io, string(line))

isleft(line::Series) = line.axis == :left
isright(line::Series) = line.axis == :right

function Base.push!(line::Series, x::Float64, y::Float64)
  push!(line.x, x)
  push!(line.y, y)
  nothing
end

function Base.push!(line::Series, x::Number, y::Number)
  push!(line, convert(Float64, x), convert(Float64, y))
end

function setdata(line::Series, x::Vector{Float64}, y::Vector{Float64})
  @assert length(x) == length(y)
  line.x = x
  line.y = y
end

function Base.empty!(line::Series)
  setdata(line, zeros(0), zeros(0))
end

updateWidgetData(line::Series, pyobj::PyObject) = pyobj[:setData](line.x, line.y)

getreglinecolor(line::Series) = line.color

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

type HeatMap <: PlotItem
  axis::Symbol
  label::AbstractString
  idx::Int
  x::Vector{Float64}
  y::Vector{Float64}
  recalcOnUpdate::Bool
  nx::Int
  ny::Int
  plt

  function HeatMap(axis::Symbol, label::AbstractString, idx::Int, nx::Int, ny::Int, plt)
    @assert axis in (:left, :right)
    new(axis, label, idx, zeros(0), zeros(0), true, nx, ny, plt)
  end
end

Base.string(heatmap::HeatMap) = "HeatMap{axis=$(heatmap.axis) label=$(heatmap.label) idx=$(heatmap.idx) npoints=$(length(heatmap.x))}"
Base.print(io::IO, heatmap::HeatMap) = print(io, string(heatmap))
Base.show(io::IO, heatmap::HeatMap) = print(io, string(heatmap))

isleft(heatmap::HeatMap) = heatmap.axis == :left
isright(heatmap::HeatMap) = heatmap.axis == :right

function Base.push!(heatmap::HeatMap, x::Float64, y::Float64)
  push!(heatmap.x, x)
  push!(heatmap.y, y)
  nothing
end

function Base.push!(heatmap::HeatMap, x::Number, y::Number)
  push!(heatmap, convert(Float64, x), convert(Float64, y))
end

function setdata(heatmap::HeatMap, x::Vector{Float64}, y::Vector{Float64})
  @assert length(x) == length(y)
  heatmap.x = x
  heatmap.y = y
end

function updateWidgetData(heatmap::HeatMap, pyobj::PyObject)
  if heatmap.recalcOnUpdate
    heatMapData = PLOT.HeatMapData(heatmap.x, heatmap.y, heatmap.nx, heatmap.ny)
    pyobj[:setData](heatMapData)
    recalcOnUpdate = false
  end
end

getreglinecolor(hm::HeatMap) = :black


text:
import StatsFuns.logsumexp
import Statistics.var

@inline function logmeanexp(array::Vector{Float64})
  return logsumexp(array) - log(length(array))
end

function testEstimator(estimator, trials)
  vals::Vector{Float64} = Vector{Float64}(undef, trials)
  totalTime = @elapsed for trial in 1:trials
    vals[trial] = estimator()
  end
  logMean::Float64 = logmeanexp(vals)
  relativeVariance::Float64 = var(exp.(vals .- logMean))
  averageTime::Float64 = totalTime / trials * Threads.nthreads()
  return averageTime, logMean, relativeVariance
end


text:
# SPDX-License-Identifier: X11
# 2020-12-07
# Day 1, Part 1

function main()
  input = parse.(Int, readlines())
  nums = falses(2020)

  for x ∈ input
    nums[x] = true
  end

  for x ∈ input
    if 2020 - x > 0 && nums[2020 - x]
      println(x * (2020 - x))
      return
    end
  end
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
function clip(l, rect::Rectangle)
    A = VoronoiDelaunay.geta(l)
    B = VoronoiDelaunay.getb(l)

    clip(A, B, rect)
end

function clip(A::VoronoiDelaunay.AbstractPoint2D, B::VoronoiDelaunay.AbstractPoint2D, rect::Rectangle)
    if isinside(A, rect) && isinside(B, rect)
        VoronoiDelaunay.Line(A, B)
        # GeometryBasics.Line(A, B)
	end

	t0 = 0.0
	t1 = 1.0
	D = B - A
	# left, right, bottom, top
	# Parametrization of line segment from A to B
	p = [-getx(D) ; getx(D) ; -gety(D) ; gety(D)]
	q = [
		getx(A) - left(rect) ;
		right(rect) - getx(A) ;
		gety(A) - lower(rect) ;
		upper(rect) - gety(A)
	]

	for k in 1:4
		if p[k] == 0.0
			# Line parallel with k'th edge
			if q[k] < 0.0
				return nothing
			end
		elseif p[k] < 0.0
			# Outside to inside
			u = q[k] / p[k]
			if u > t1
				return nothing
			end
			t0 = max(u, t0)
		else # p[k] > 0.0
			# Inside to outside
			u = q[k] / p[k]
			if u < t0
				return nothing
			end
			t1 = min(u, t1)
		end
	end

	return VoronoiDelaunay.Line(A + t0*D, A + t1*D)
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL; GL = ViewerGL

# input of primitive shapes
W,(WW,EW,FW) = Lar.cuboid([1,1], true)
square = (W,EW)
# hierarchical assembly
model = Lar.Struct([
            Lar.Struct([
                square,
                Lar.t(.175,.175), Lar.s(.25,.25),
                square
            ]),
            Lar.t(.825,.825), Lar.s(-.25,-.25),
            square
        ])
# visualization of generated wire-frame model
V,EV = Lar.struct2lar(model)
VV = [[k] for k=1:size(V,2)]
GL.VIEW(GL.numbering(.4)((V,[VV, EV])));

triangles = Lar.triangulate2d(V, EV)
GL.VIEW([ GL.GLGrid(V,triangles) ]);

# generate edges
ev = map(sort,cat([[[u,v], [v,w], [w,u]] for (u,v,w) in triangles]))
# remove duplicated edges from triangulation
nodupev = collect(Set(ev))
GL.VIEW(GL.numbering(.25)((V,[VV, nodupev, triangles]),GL.COLORS[1],0.2));


text:
printstyled("  poly\n", color=:light_green)
import SeisIO:polyfit, polyval
N = 1000
C = 2.0
P = 3
at = 1.0e-5

t = collect(1.0:1.0:Float64(N))
x = C.*t

pf = polyfit(t, x, 1)
@test isapprox(pf[1], C, atol=at)
@test isapprox(pf[2], 0.0, atol=at)

pf = polyfit(t, x, 3)
@test isapprox(pf[1], 0.0, atol=at)
@test isapprox(pf[2], 0.0, atol=at)
@test isapprox(pf[3], C, atol=at)
@test isapprox(pf[4], 0.0, atol=at)

pv = polyval(pf, t)
@test isapprox(pv, x)

y = t.^P
pf = polyfit(t, y, P)
@test isapprox(pf[1], 1.0, atol=at)
for i = 2:P+1
  @test isapprox(pf[i], 0.0, atol=at)
end

# this is the current response of CC.VALT..BHE as of 2014,079,00:00:00
resp = PZResp(a0 = 9.092142f11, f0 = 1.0f0,
  p = Complex{Float32}[-0.1486+0.1486im, -0.1486-0.1486im, -391.96+850.69im, -391.96-850.69im, -471.24+0.0im, -2199.1+0.0im],
  z = Complex{Float32}[0.0+0.0im, 0.0+0.0im])

# test: recover a0 from poles and zeroes in response file using above poly, polyval, polyfit functions
T = typeof(resp.a0)
Z = poly(resp.z)
P = poly(resp.p)
s = Complex{T}(2*pi*im*resp.f0)
a0 = one(T)/T(abs(polyval(Z, s)/polyval(P, s)))
@test isapprox(a0, resp.a0)

# test: adapted from Octave
r = Float64.(collect(0:10:50))
p = poly(r)
p ./= maximum(abs.(p))
x = Float64.(collect(0:5:50))
y = polyval(p,x) + 0.25*sin.(100.0*x)
y2 = similar(y)
for i = 1:length(y)
  y2[i] = polyval(p, x[i]) + 0.25*sin(100.0*x[i])
end
pf = polyfit(x, y, length(r))
y3 = polyval(pf, x)
y_expect = [0.00000, -1.34741,  0.20672,  0.16168,  0.23251, -0.45550,  0.05480,  0.47582, -0.17088, -0.99408, -0.24699]
y3_expect = [-0.021348, -1.234641, -0.027165,  0.412554,  0.035480, -0.273558, -0.034148,  0.358771,  0.031322, -1.105561, -0.225033]

@test isapprox(y, y_expect, atol=1.0e-5)
@test isapprox(y2, y_expect, atol=1.0e-5)
@test isapprox(y3, y3_expect, atol=1.0e-5)