text:
mutable struct Descriptor <: AbstractDescriptor
    p::libgb.GrB_Descriptor
    function Descriptor(p::libgb.GrB_Descriptor)
        d = new(p)
        function f(descriptor)
            libgb.GrB_Descriptor_free(Ref(descriptor.p))
        end
        return finalizer(f, d)
    end
end

function Descriptor()
    return Descriptor(libgb.GrB_Descriptor_new())
end

Base.unsafe_convert(::Type{libgb.GrB_Descriptor}, d::Descriptor) = d.p

function Base.getproperty(d::Descriptor, s::Symbol)
    if s == :p
        return getfield(d, s)
    elseif s == :output
        f = libgb.GrB_OUTP
    elseif s == :mask
        f = libgb.GrB_MASK
    elseif s == :input1
        f = libgb.GrB_INP0
    elseif s == :input2
        f = libgb.GrB_INP1
    #elseif s == :nthreads
        #f = libgb.GxB_DESCRIPTOR_NTHREADS
    #elseif s == :chunk
        #f = libgb.GxB_DESCRIPTOR_CHUNK
    #elseif s == :axb_method
        #f = libgb.GxB_AxB_METHOD
    #elseif s == :sort
        #f = libgb.GxB_SORT
    else
        error("type Descriptor has no field $s")
    end
    return libgb.GxB_Descriptor_get(d, f)
end

function Base.setproperty!(d::Descriptor, s::Symbol, x)
    if s == :p
        setfield!(d, s, v)
    elseif s == :output
        f = libgb.GrB_OUTP
    elseif s == :mask
        f = libgb.GrB_MASK
    elseif s == :input1
        f = libgb.GrB_INP0
    elseif s == :input2
        f = libgb.GrB_INP1
    #elseif s == :nthreads
        #f = libgb.GxB_DESCRIPTOR_NTHREADS
    #elseif s == :chunk
        #f = libgb.GxB_DESCRIPTOR_CHUNK
    #elseif s == :axb_method
        #f = libgb.GxB_AxB_METHOD
    #elseif s == :sort
        #f = libgb.GxB_SORT
    else
        error("type Descriptor has no field $s")
    end
    libgb.GrB_Descriptor_set(d, f, x)
end

function Base.:+(d1::Descriptor, d2::Descriptor)
    d = Descriptor()
    for f ∈ propertynames(d)
        if f == :input1
            if getproperty(d1, f) == Descriptors.TRANSPOSE && getproperty(d2, f) == Descriptors.TRANSPOSE
                setproperty!(d, f, Descriptors.DEFAULT)
            elseif getproperty(d1, f) == Descriptors.TRANSPOSE || getproperty(d2, f) == Descriptors.TRANSPOSE
                setproperty!(d, f, Descriptors.TRANSPOSE)
            end
        elseif f == :input2
            if getproperty(d1, f) == Descriptors.TRANSPOSE && getproperty(d2, f) == Descriptors.TRANSPOSE
                setproperty!(d, f, Descriptors.DEFAULT)
            elseif getproperty(d1, f) == Descriptors.TRANSPOSE || getproperty(d2, f) == Descriptors.TRANSPOSE
                setproperty!(d, f, Descriptors.TRANSPOSE)
            end
        else
            if getproperty(d1, f) != Descriptors.DEFAULT
                setproperty!(d, f, getproperty(d1, f))
            end
            if getproperty(d2, f) != Descriptors.DEFAULT
                setproperty!(d, f, getproperty(d2, f))
            end
        end
    end
    return d
end

#This is probably not ideal. Perhaps kwargs = nothing by default is better
Base.:+(d1::Descriptor, ::Nothing) = d1
Base.:+(::Nothing, d2::Descriptor) = d2

function Base.propertynames(d::Descriptor)
    return (
    :output,
    :mask,
    :input1,
    :input2,
    #:nthreads,
    #:chunk,
    #:axb_method,
    #:sort,
    )
end
baremodule Descriptors
import ..SSGrB: load_global, Descriptor
import ..libgb:
GB_Descriptor_opaque,
GrB_Descriptor,
GxB_DEFAULT,
GrB_REPLACE,
GrB_COMP,
GrB_STRUCTURE,
GrB_TRAN,
GxB_GPU_ALWAYS,
GxB_GPU_NEVER,
GxB_AxB_GUSTAVSON,
GxB_AxB_DOT,
GxB_AxB_HASH,
GxB_AxB_SAXPY
import ..Types
import Base.C_NULL

const DEFAULT = GxB_DEFAULT
const REPLACE = GrB_REPLACE
const COMPLEMENT = GrB_COMP
const STRUCTURE = GrB_STRUCTURE
const TRANSPOSE = GrB_TRAN
const GUSTAVSON = GxB_AxB_GUSTAVSON
const DOT = GxB_AxB_DOT
const HASH = GxB_AxB_HASH
const SAXPY = GxB_AxB_SAXPY
end

function loaddescriptors()
    builtins = ["GrB_DESC_T1",
    "GrB_DESC_T0",
    "GrB_DESC_T0T1",
    "GrB_DESC_C",
    "GrB_DESC_CT1",
    "GrB_DESC_CT0",
    "GrB_DESC_CT0T1",
    "GrB_DESC_S",
    "GrB_DESC_ST1",
    "GrB_DESC_ST0",
    "GrB_DESC_ST0T1",
    "GrB_DESC_SC",
    "GrB_DESC_SCT1",
    "GrB_DESC_SCT0",
    "GrB_DESC_SCT0T1",
    "GrB_DESC_R",
    "GrB_DESC_RT1",
    "GrB_DESC_RT0",
    "GrB_DESC_RT0T1",
    "GrB_DESC_RC",
    "GrB_DESC_RCT1",
    "GrB_DESC_RCT0",
    "GrB_DESC_RCT0T1",
    "GrB_DESC_RS",
    "GrB_DESC_RST1",
    "GrB_DESC_RST0",
    "GrB_DESC_RST0T1",
    "GrB_DESC_RSC",
    "GrB_DESC_RSCT1",
    "GrB_DESC_RSCT0",
    "GrB_DESC_RSCT0T1"]
    for name ∈ builtins
        simple = Symbol(string(name[10:end]))
        constquote = quote
            const $simple = Descriptor(load_global($name, GB_Descriptor_opaque))
        end
        @eval(Descriptors, $constquote)
    end

end
isloaded(d::AbstractDescriptor) = d.p !== C_NULL

function printtest(::IO, M::Descriptor)
    str = mktemp() do _, f
        cf = Libc.FILE(f)
        libgb.GxB_Descriptor_fprint(M, "", libgb.GxB_SHORT, cf)
        flush(f)
        seekstart(f)
        x = read(f, String)
        close(cf)
        x
    end
    print(str)
end

function Base.show(io::IO, d::Descriptor)
    printtest(io, d)
end


text:
module Roxygen
    export roxygenize

    include("constants.jl")
    include("types.jl")
    include("chunker.jl")
    include("lexer.jl")
    include("semantics.jl")
    include("markdown.jl")
    include("roxygenize.jl")
end

ext: jl


text:
## ExoplanetsSysSim/examples/run_abc.jl
## (c) 2015 Eric B. Ford

include("abc_setup.jl")

using SysSimABC

abc_plan = setup_abc(3)
@time output = SysSimABC.run_abc(abc_plan)

#@time output2 = SysSimABC.run_abc(abc_plan,output)
#@time output3 = SysSimABC.run_abc(abc_plan,output2)


text:
# Gillespie simulations wrapper
struct Gillespie end

"""
#Remarks
We could think about a generalization of this to pass also custom events, threhold, etc. But then this exactly becomes the event_handler function ... So I would keep this as a specialization.

We should however move events to the function argument! Because this is relevant for the update function

#See also direct construction of event handlers
[`ListEventRateSimple`](@ref), [`ListEventRateActiveMask`](@ref)
"""
function initialize(alg::Gillespie, rates::Vector{Float64}; type_event_handler::String = "ListEventRateSimple")
  events = collect(1:length(rates))
  noevent = 0
  threshold = 0.0
  if type_event_handler == "ListEventRateSimple"
   event_handler = ListEventRateSimple{Int}(events, rates, threshold, noevent)
  elseif type_event_handler == "ListEventRateActiveMask"
    event_handler = ListEventRateActiveMask{Int}(events, rates, threshold, noevent)
  else
    throw(UndefVarError(:type_event_handler))
  end

  return event_handler
end

function next(alg::Gillespie, rng::AbstractRNG, event_handler::AbstractEventHandlerRate)::Tuple{Float64,Int}
  next(KineticMonteCarlo(), rng, event_handler)
end

function advance!(alg::Gillespie, rng::AbstractRNG, event_handler::AbstractEventHandlerRate, update!::Function, total_time::T)::T where {T <: AbstractFloat}
  advance!(KineticMonteCarlo(), rng, event_handler, update!, total_time)
end

function advance!(alg::Gillespie, event_handler::AbstractEventHandlerRate, update!::Function, total_time::T)::T where {T <: AbstractFloat}
  advance!(alg, Random.GLOBAL_RNG, event_handler, update!, total_time)
end


text:
module EvoTrees

export init_evotree, grow_evotree!, grow_tree, fit_evotree, predict,
    fit_evotree_gpu, predict_gpu,
    EvoTreeRegressor, EvoTreeCount, EvoTreeClassifier, EvoTreeGaussian,
    EvoTreeRModels, importance

using Statistics
using Base.Threads: @threads
using StatsBase: sample, quantile
using Random
using StaticArrays
using Distributions
using CategoricalArrays
using CUDA
import MLJModelInterface
import MLJModelInterface: fit, update
import MLJModelInterface: predict

include("models.jl")
include("structs.jl")
include("loss.jl")
include("eval.jl")
include("predict.jl")
include("find_split.jl")
include("fit.jl")
include("importance.jl")
include("MLJ.jl")

include("structs_gpu.jl")
include("loss_gpu.jl")
include("eval_gpu.jl")
include("predict_gpu.jl")
include("find_split_gpu.jl")
include("fit_gpu.jl")

end # module


text:
baremodule ShowCode

export @sc_intel,
    @sc_ircode,
    @sc_llvm,
    @sc_native,
    sc_intel,
    sc_ircode,
    sc_llvm,
    sc_native

#=
macro sc_sexpr end
macro sc_lowered end
macro sc_typed end
=#
macro sc_ircode end
macro sc_llvm end
macro sc_native end
macro sc_intel end

function sc_ircode end
function sc_llvm end
function sc_native end
function sc_intel end

function from_llvm end
function from_native end

module Implementations

using ..ShowCode: ShowCode
import ..ShowCode:
    @sc_intel,
    @sc_ircode,
    @sc_llvm,
    @sc_native,
    sc_intel,
    sc_ircode,
    sc_llvm,
    sc_native

import AbstractTrees
import LLVM_jll
using Accessors: @set
using InteractiveUtils:
    InteractiveUtils,
    code_llvm,
    code_lowered,
    code_native,
    code_typed,
    gen_call_with_extracted_types_and_kwargs,
    print_llvm,
    print_native
using Requires: @require
using UnPack: @unpack

const DetachNode = try
    Core.DetachNode
catch
    Union{}
end

const ReattachNode = try
    Core.ReattachNode
catch
    Union{}
end

const SyncNode = try
    Core.SyncNode
catch
    Union{}
end

include("compilerutils.jl")
include("core.jl")
include("code_ircode.jl")
include("ircode.jl")
include("llvm.jl")
include("native.jl")

function __init__()
    @require Cthulhu="f68482b8-f384-11e8-15f7-abe071a5a75f" include("cthulhu.jl")
end

end

const CONFIG = Implementations.ShowCodeConfig()

Implementations.finalize_module()

end


text:
# Bogumił Kamiński, 2019-03-25

using DataFrames

lines = readlines("sample_dataframes.log")

author = r"^Author: ([^<]*) <"
insc = r"^.+changed, ([0-9]+) insertion"
delc = r"^.+changed.*, ([0-9]+) deletion"
authordata = DataFrame(author=String[], action=String[], count=Int[])
curauthor = ""

for line in lines
    m = match(author, line)
    if m !== nothing
        global curauthor = m[1]
    end
    m = match(insc, line)
    m === nothing || push!(authordata,
                           [curauthor, "insertion", parse(Int, m[1])])
    m = match(delc, line)
    m === nothing || push!(authordata,
                           [curauthor, "deletion", parse(Int, m[1])])
end

authorstats = by(authordata, [:author, :action], count = :count=>sum)
df = unstack(authorstats, :action, :count)
df.all = coalesce.(df.deletion, 0) .+ coalesce.(df.insertion, 0)
display(sort!(df, :all, rev=true))

max_stars_repo_path: 5. Logparser/logparser_1.jl

max_issues_repo_path: 5. Logparser/logparser_1.jl

max_forks_repo_path: 5. Logparser/logparser_1.jl


text:
module Opus
using Ogg, Opus_jll
using FileIO
import Base: convert, show, write

export OpusDecoder, OpusEncoder, OpusArray, load, save

include("defines.jl")
include("decoder.jl")
include("encoder.jl")
include("opusarray.jl")

end # module


text:
#import Requests
import HTTPClient

dict2t(d) = map(p -> (p.first, p.second), d)

headers = Dict("User-Agent" => "EVE.jl/$EVE_JL_VERSION (https://github.com/Daagr/EVE.jl)")
function eget(url; query = Dict())
    r = HTTPClient.get(url; headers = dict2t(headers), query_params = dict2t(query))

    if r.http_code ≠ 200
        error("Request failed $r.http_code")
    end
    bytestring(r.body)
end


text:
module Properties
using Base.FastMath: add_fast, mul_fast, min_fast, max_fast

export return_type, initial, eltype_bound, instance
export Guard, Assume

"""
    initial(op, T)

Return the default initial value for the reduction operator `op` on the type
`T`. Return `nothing` if you such a value does not exist.

See also: [`zero`](@ref), [`oneunit`](@ref).

# Examples
```jldoctest
julia> initial(+, Int)
0
julia> initial(+, Real)
false
```
"""

@inline initial(f, T) = nothing

@inline initial(::typeof(+), T::Type{<:Number}) = Some(zero(T))
@inline initial(::typeof(add_fast), T::Type{<:Number}) = Some(zero(T))
@inline initial(::typeof(*), T::Type{<:Number}) = Some(oneunit(T))
@inline initial(::typeof(mul_fast), T::Type{<:Number}) = Some(oneunit(T))
@inline initial(::typeof(max), T::Type{<:Number}) = Some(typemin(T))
@inline initial(::typeof(max_fast), T::Type{<:Number}) = Some(typemin(T))
@inline initial(::typeof(min), T::Type{<:Number}) = Some(typemax(T))
@inline initial(::typeof(min_fast), T::Type{<:Number}) = Some(typemax(T))

"""
    instance(T::Type)

Return an instance of the type if there is only one instance.
Return `nothing` otherwise.

# Examples
```jldoctest
julia> instance(IndexCartesian)
Some(IndexCartesian())
julia> instance(AbstractArray)
nothing
julia> instance(Nothing)
Some(nothing)
```
"""
@inline function instance(T::Type)
    if isbitstype(T) && sizeof(T)==0 && :instance in propertynames(T)
        return Some(T.instance)
    else
        return nothing
    end
end

"""
    return_type(f, args...)

Returns a type `T` containing the return type of the function `f` called with
arguments of type `args...`. Should run in time linear with the description of
the types of its arguments.

See also: [`Core.Compiler.return_type`](@ref).

# Examples
```jldoctest
julia> return_type(+, Int, Int)
Int64
julia> return_type(+, Real, Real)
Any
```
"""
@inline function return_type(f, args...)
    return Core.Compiler.return_type(f, args)
end

Rational, Complex

const AtomNumber = Union{BigFloat,
                           Float16,
                           Float32,
                           Float64,
                           Bool,
                           BigInt,
                           Int128,
                           Int16,
                           Int32,
                           Int64,
                           Int8,
                           UInt128,
                           UInt16,
                           UInt32,
                           UInt64,
                           UInt8}

const MoleculeNumber = Union{AtomNumber, Complex{<:AtomNumber}, Rational{<:AtomNumber}}

@inline return_type(::typeof(+), a::Type{<:MoleculeNumber}, b::Type{<:MoleculeNumber}) = promote_type(a, b)
@inline return_type(::typeof(add_fast), a, b) = return_type(+, a, b)
@inline return_type(::typeof(*), a::Type{<:MoleculeNumber}, b::Type{<:MoleculeNumber}) = promote_type(a, b)
@inline return_type(::typeof(mul_fast), a, b) = return_type(*, a, b)
@inline return_type(::typeof(max), a, b) = Union{a, b}
@inline return_type(::typeof(max_fast), a, b) = Union{a, b}
@inline return_type(::typeof(min), a, b) = Union{a, b}
@inline return_type(::typeof(min_fast), a, b) = Union{a, b}

"""
    eltype_bound(A)

Return the strictest bound on the element type of `A` that is possible to obtain
in constant time. This may be stricter than `eltype(A)`, which corresponds to
the declared element type of `A`

See also: [`eltype`](@ref).
"""
@inline eltype_bound(arg) = eltype(arg)

#declare opposite

struct Guard{Op}
    op::Op
end

(op::Guard)(x::Nothing, y) = y
(op::Guard)(x, y) = op.op(x, y)
@inline return_type(op::Guard, T, S) = return_type(op.op, T, S)
@inline return_type(op::Guard, ::Type{Union{Nothing, T}}, S) where {T} = return_type(op.op, T, S)
@inline return_type(op::Guard, ::Type{Nothing}, S) = S

@inline initial(::Guard, ::Any) = Some(nothing)

struct Assume{T, Op}
    op::Op
end

@inline ((op::Assume{T})(x, y)::T) where {T} = op.op(x, y)
@inline function return_type(op::Assume{T}, S, R)::T where {T}
    Q = Properties.return_type(op.op, S, R)
    if Q <: T
        return Q
    else
        return T
    end
end

@inline initial(op::Assume, T) = initial(op.op, T)

end


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

using ValueShapes
using Test

using LinearAlgebra

@testset "functions" begin
    input_shape = NamedTupleShape(
        a = ArrayShape{Real}(3,2),
        b = ArrayShape{Real}(2)
    )

    output_shape = NamedTupleShape(
        x = ScalarShape{Real}(),
        y = ArrayShape{Real}(3)
    )

    @testset "UnshapedFunction" begin
        f(x) = (x = norm(x.a)^2 + norm(x.b)^2, y = vec(sum(x.a, dims = 2)))
        ux = [2, 7, 4, 6, 8, 5, 3, 1]
        x = input_shape(ux)[]
        y = f(x)
        uy = unshaped(y, output_shape)
        @test @inferred(unshaped(f, input_shape)(ux)) == y
        @test @inferred(unshaped(f, input_shape, nothing)(ux)) == y
        @test @inferred(unshaped(f, input_shape, output_shape)(ux)) == uy
        @test @inferred(unshaped(f, nothing, output_shape)(x)) == uy
        @test @inferred(unshaped(f, nothing)(x)) == y
        @test @inferred(unshaped(f, nothing, nothing)(x)) == y
    end

    @testset "FuncWithVarShape" begin
        x_flat = rand(@inferred totalndof(input_shape))
        x = @inferred stripscalar(input_shape(x_flat))

        g(x) = sum(x.a * x.b)

        @test @inferred(input_shape >> g) isa ValueShapes.FuncWithVarShape
        @test @inferred(varshape((input_shape >> g))) == input_shape
        @test @inferred(vardof((input_shape >> g))) == totalndof(input_shape)

        fws = input_shape >> g
        @test g(x) == @inferred(fws(x))
        @test g(x) == @inferred(fws(x_flat))
    end
end


text:
"""
    convert_samplepath(X)

Useful for storage of a samplepath of states
Ordering is as follows:
1) time
2) landmark nr
3) for each landmark: q1, q2 p1, p2

With m time-steps, n landmarks, this entails a vector of length m * n * 2 * d
"""
function convert_samplepath(X)
    VA = VectorOfArray(map(x->deepvec(x),X.yy))
    vec(convert(Array,VA))
end

"""
    convert_samplepath(Xvec::Vector)

Useful for storage of a samplepath of states
Ordering is as follows:
0) shape
1) time
2) landmark nr
3) for each landmark: q1, q2 p1, p2

With m time-steps, n landmarks, this entails a vector of length m * n * 2 * d
"""
function convert_samplepath(Xvec::Vector)
    nshapes = length(Xvec)
    out = [convert_samplepath(Xvec[1])]
    for k in 2:nshapes
        push!(out,convert_samplepath(Xvec[k]))
    end
    vcat(out...)
end

"""
    extract_initial_and_endstate(iteratenr::Int64, Xpath)

Extract from a sample path its states at time zero and time T (endtime)

iteratnr:: integer that indicates iteration.nr in mcmc-algorithm
Xpath:: Bridge.Samplepath

Returns a matrix with first column iteratenr repeated, second column state at time zero, third column state a time T

The second and third column as ordered as follows (in case of two dimensional landmarks):
landmark 1 pos1
landmark 1 pos2
landmark 1 mom1
landmark 1 mom2

landmark 2 pos1
landmark 2 pos2
landmark 2 mom1
landmark 2 mom2

...
"""
function extract_initial_and_endstate(iteratenr::Int64, Xpath)
    X0 = deepvec(Xpath.yy[1])  # state at time zero
    XT = deepvec(Xpath.yy[end]) # state at time T
    ℓ = length(XT)
    out = [fill(iteratenr, ℓ) X0 XT]
    out
end

function write_output(xobs0, xobsT, parsave, Xsave, initendstates_save, elapsed, accinfo,tt, n,nshapes, subsamples,ITER, updatescheme, Σobs, pars, ρ, δ,P, outdir)
    println("Elapsed time: ",round(elapsed/60;digits=2), " minutes")
    ave_acc = [mean(x) for x in eachcol(accinfo[!,1:end-1])]
    println("Average acceptance of updatesteps: ", round.(ave_acc;digits=2))
    write_mcmc_iterates(Xsave, tt, n, nshapes, subsamples, outdir)
    write_info(ITER, n, tt, updatescheme, Σobs, pars, ρ, δ , ave_acc, elapsed, outdir)
    write_observations(xobs0, xobsT, n, nshapes, outdir)
    write_acc(accinfo,outdir)
    write_params(parsave, 0:ITER, outdir)
    write_noisefields(P, outdir)
    write_initendstates(initendstates_save, n, outdir)
end

"""
    write_mcmc_iterates(Xsave, tt_, n, nshapes, subsamples, outdir)

Write bridge iterates to file

## Arguments:
- `Xsave`: contains values of bridges
- `tt_`: grid on which bridges are simulated
- `n`: nr of bridges
- `nshapes`: nr of shapes
- `subsamples`: indices of iterates that are kept and saved in Xsave
- `outdir`: output directory
"""
function write_mcmc_iterates(Xsave, tt_, n, nshapes, subsamples, outdir)
    outdir[end] == "/" && error("provide pathname without trailing '/'")
    iterates = reshape(vcat(Xsave...),2*d*length(tt_)*n*nshapes, length(subsamples)) # each column contains samplepath of an iteration
    # total number of rows is nshapes * length(tt_) * n * (2d)
    # Ordering in each column is as follows:
    # 0) shape
    # 1) time
    # 2) landmark nr
    # 3) for each landmark: q1, q2 p1, p2
    #if !(d==2) error("pqtype in write_mcmc_iterates only implemented in case d=2") end
    shapes = repeat(1:nshapes, inner=length(tt_)*2d*n)
        #times = repeat(tt_,inner=2d*n*nshapes)  # original buggy version when nshapes=1
    times = repeat(tt_, inner=2*d*n, outer=nshapes)
    landmarkid = repeat(1:n, inner=2d, outer=length(tt_)*nshapes)
    if d==1
        pqtype = repeat(["pos1", "mom1"], outer=length(tt_)*n*nshapes)
    elseif d==2
        pqtype = repeat(["pos1", "pos2", "mom1", "mom2"], outer=length(tt_)*n*nshapes)
    elseif d==3
        pqtype = repeat(["pos1", "pos2", "pos3", "mom1", "mom2", "mom3"], outer=length(tt_)*n*nshapes)
    end
    out = hcat(times,pqtype,landmarkid,shapes,iterates)
    headline = "time " * "pqtype " * "landmarkid " * "shapes " * prod(map(x -> "iter"*string(x)*" ",subsamples))
    headline = chop(headline,tail=1) * "\n"

    fn = joinpath(outdir, "iterates.csv")
    f = open(fn,"w")
    write(f, headline)
    writedlm(f,out)
    close(f)
end

"""
    Write info to txt file
"""
function write_info(ITER, n, tt_, updatescheme, Σobs, tp, ρ, δ, perc_acc_pcn, elapsed, outdir)
    outdir[end] == "/" && error("provide pathname without trailing '/'")

    f = open(joinpath(outdir,"info.txt"),"w")
    write(f, "Number of iterations: ",string(ITER),"\n")
    write(f, "Number of landmarks: ",string(n),"\n")
    write(f, "Length time grid: ", string(length(tt_)),"\n")
    write(f, "Endpoint: ",string(tt_[end]),"\n")
    write(f, "updatescheme: ", string(updatescheme),"\n")
    write(f, "Average acceptance rates update steps: ",string(perc_acc_pcn),"\n\n")
    write(f, "Noise Sigma: ",string(Σobs),"\n")
    write(f, "tuningpars_mcmc: ", string(tp),"\n")
    write(f, "Final value  of rho (Crank-Nicholsen parameter: ",string(ρ),"\n")
    write(f, "Final value of MALA parameter (delta): ",string(δ),"\n")
    write(f, "skip in evaluation of loglikelihood: ",string(sk),"\n")
    write(f, "Elapsed time: ", string(elapsed),"\n")
    close(f)
end

"""
    write_acc(accinfo,outdir)

Write dataframe with acceptance numbers of updatesteps to a csv file
"""
function write_acc(accinfo,outdir)
    outdir[end] == "/" && error("provide pathname without trailing '/'")
    CSV.write(joinpath(outdir, "accdf.csv"), accinfo; delim=";")
end

"""
    write_observations(xobs0, xobsT, n, nshapes, outdir)

Write observations to csv file
"""
function write_observations(xobs0, xobsT, n, nshapes, outdir)
    outdir[end] == "/" && error("provide pathname without trailing '/'")

    valueT = vcat(map(x->deepvec(x), xobsT)...) # merge all observations at time T in one vector
    if d==1
        posT = repeat(["pos1"], n*nshapes)
    elseif d==2
        posT = repeat(["pos1","pos2"], n*nshapes)
    elseif d==3
        posT = repeat(["pos1","pos2","pos3"], n*nshapes)
    end
    shT = repeat(1:nshapes, inner=d*n)
    obsTdf = DataFrame(pos=posT,shape=shT, value=valueT,landmark=repeat(1:n,inner=d,outer=nshapes))

    if d==1
        obs0df = DataFrame(pos1=extractcomp(xobs0,1),landmark=1:n)
    elseif d==2
        obs0df = DataFrame(pos1=extractcomp(xobs0,1), pos2=extractcomp(xobs0,2),landmark=1:n)
    elseif d==3
        obs0df = DataFrame(pos1=extractcomp(xobs0,1), pos2=extractcomp(xobs0,2), pos3=extractcomp(xobs0,3),landmark=1:n)
    end

    CSV.write(joinpath(outdir, "obs0.csv"), obs0df; delim=";")
    CSV.write(joinpath(outdir, "obsT.csv"), obsTdf; delim=";")
end

"""
    write_params(parsave,subsamples,outdir)

Write parameter iterates to file
"""
function write_params(parsave,subsamples,outdir)
    outdir[end] == "/" && error("provide pathname without trailing '/'")

    parsdf = DataFrame(a=extractcomp(parsave,1),#c=extractcomp(parsave,2),
            gamma=extractcomp(parsave,2), iterate=subsamples)
    CSV.write(joinpath(outdir, "parameters.csv"), parsdf; delim=";")
end

"""
    write_noisefields(P,outdir)

Write noisefields to file (empty in case of MS-model)
"""
function write_noisefields(P,outdir)
    outdir[end] == "/" && error("provide pathname without trailing '/'")

    if isa(P,Landmarks)
        nfsloc = [P.nfs[j].δ for j in eachindex(P.nfs)]
        if d==2
            nfsdf = DataFrame(locx =  extractcomp(nfsloc,1),
                          locy =  extractcomp(nfsloc,2),
                          nfstd=fill(P.nfstd,length(P.nfs)))
        elseif d==1
            nfsdf = DataFrame(locx =  extractcomp(nfsloc,1),
                          nfstd=fill(P.nfstd,length(P.nfs)))
        end
    elseif isa(P,MarslandShardlow)
        nfsdf =DataFrame(locx=Int64[], locy=Int64[], nfstd=Int64[])
    end
    CSV.write(joinpath(outdir, "noisefields.csv"), nfsdf; delim=";")
end

"""
    write_initendstates(initendstates_save)

initendstates_save is an Array of matrices, each matrix corresponding to one iteration,
with ordering as specified in `extract_initial_and_endstate`
"""
function write_initendstates(initendstates_save, n, outdir)
    a = initendstates_save

    A = DataFrame(vcat(a...))
    rename!(A, ["iterate", "time0", "time1"])
    N = div(nrow(A),2d)
    if d==1
        lab = repeat(["pos1", "mom1"], N)
    elseif d==2
        lab = repeat(["pos1","pos2", "mom1", "mom2"], N)
    elseif d==3
        lab = repeat(["pos1","pos2","pos3", "mom1", "mom2", "mom3"], N)
    end
    A[!,:type] = lab
    A[!,:iterate] = Int.(A[!,:iterate])
    landmarkid = repeat(1:n, inner=2d, outer=length(a))
    A[!,:landmarkid] = landmarkid
    CSV.write(joinpath(outdir, "initendstates.csv"), A; delim=";")
end


text:
struct FibonacciIterable{I}
    s0::I
    s1::I
end

@testset "IterationTools" begin

    using Printf
    using ProximalAlgorithms: IterationTools
    using Random

    Random.seed!(0)

    Base.iterate(iter::FibonacciIterable, state=(iter.s0, iter.s1)) = state[1], (state[2], sum(state))

    @testset "Looping" begin
        iter = rand(Float64, 10)
        last = @inferred IterationTools.loop(iter)
        @test last == iter[end]
    end

    @testset "Halting" begin
        iter = [
            0,
            1,
            1,
            2,
            3,
            5,
            8,
            13,
            21,
            34,
            55,
            89,
            144,
            233,
            377,
            610,
            987,
            1597,
            2584,
            4181,
            6765,
        ]
        truncated = IterationTools.halt(iter, x -> x >= 1000)
        @test Base.IteratorSize(truncated) == Base.IteratorSize(iter)
        @test eltype(truncated) == eltype(iter)
        last = @inferred IterationTools.loop(truncated)
        @test last == 1597
    end

    @testset "Side effects" begin
        disp(s) = @info "$s"

        iter = IterationTools.tee(FibonacciIterable(BigInt(0), BigInt(1)), disp)

        fib = [0, 1, 1, 2, 3, 5, 8, 13, 21, 34]

        _, state = @test_logs (:info, "$(fib[1])") iterate(iter)
        for k = 2:10
            _, state = @test_logs (:info, "$(fib[k])") iterate(iter, state)
        end
    end

    @testset "Sampling" begin
        iter = randn(Float64, 147)
        sample = IterationTools.sample(iter, 10)
        @test eltype(sample) == Float64
        @test length(sample) == 15
        k = 0
        for x in sample
            idx = min(147, (k + 1) * 10)
            @test x == iter[idx]
            k = k + 1
        end
    end

    @testset "Timing" begin
        iter = randn(Float64, 10)
        timed = IterationTools.stopwatch(iter)
        @test length(timed) == length(iter)
        @test axes(timed) == axes(iter)
        @test eltype(timed) == (UInt64, eltype(iter))
        k = 0
        for (t, x) in timed
            @test x == iter[k+1]
            @test t >= k * 1e8 * 0.9 # 1e8 ns = 0.1 s, factor 0.9 is for safety
            sleep(0.1)
            k = k + 1
        end
    end

end


text:
using Pkg
Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
using CSV, MLJ, DataFrames, Random

function generate(; option = "drop", std = "false", valid = "true", test = "true")

    data = CSV.read(joinpath(@__DIR__, "data", "trainingdata.csv"), DataFrame);
    test_all = CSV.read(joinpath(@__DIR__, "data", "testdata.csv"), DataFrame)

    Random.seed!(2809)

    if option == "drop"
        new_data = dropmissing(data);
        coerce!(new_data, :precipitation_nextday => Multiclass);
        idxs = randperm(size(new_data, 1));
        dim = size(new_data)[1];
        idx = Vector{Int}([1, floor(dim/3), floor(2*dim/3), dim]);
        if std == "true"
            std_mach = machine(Standardizer(), select(new_data, Not([:ZER_sunshine_1, :ALT_sunshine_4])))
            fit!(std_mach)
            new_data = MLJ.transform(std_mach, select(new_data, Not([:ZER_sunshine_1, :ALT_sunshine_4])))
            if valid == "true"
                if test == "true"
                    new_data = (train = new_data[idxs[idx[1]:idx[2]], :], valid = new_data[idxs[idx[2]:idx[3]], :], test = new_data[idxs[idx[3]:idx[4]], :]);
                    new_data
                else
                    new_data = DataFrame()
                    new_data
                end
            else
                if test == "true"
                    new_data = (train = new_data[idxs[idx[1]:idx[3]], :], test = new_data[idxs[idx[3]:idx[4]], :]);
                    new_data
                else
                    new_test = MLJ.transform(std_mach, select(test_all, Not([:ZER_sunshine_1, :ALT_sunshine_4])));
                    new_data, new_test
                end
            end
        else
            if valid == "true"
                if test == "true"
                    new_data = (train = new_data[idxs[idx[1]:idx[2]], :], valid = new_data[idxs[idx[2]:idx[3]], :], test = new_data[idxs[idx[3]:idx[4]], :]);
                    new_data
                else
                    new_data = DataFrame()
                    new_data
                end
            else
                if test == "true"
                    new_data = (train = new_data[idxs[idx[1]:idx[3]], :], test = new_data[idxs[idx[3]:idx[4]], :]);
                    new_data
                else
                    new_data, test_all
                end
            end
        end
    elseif option == "med"
        new_data = MLJ.transform(fit!(machine(FillImputer(), select(data, Not(:precipitation_nextday)))), select(data, Not(:precipitation_nextday)));
        insertcols!(new_data,:precipitation_nextday=>data.precipitation_nextday[:]);
        coerce!(new_data, :precipitation_nextday => Multiclass);
        idxs = randperm(size(new_data, 1));
        dim = size(new_data)[1];
        idx = Vector{Int}([1, floor(dim/3), floor(2*dim/3), dim]);
        if std == "true"
            std_mach = machine(Standardizer(), select(new_data, Not([:ZER_sunshine_1, :ALT_sunshine_4])))
            fit!(std_mach)
            new_data = MLJ.transform(std_mach, select(new_data, Not([:ZER_sunshine_1, :ALT_sunshine_4])))
            if valid == "true"
                if test == "true"
                    #If valid is true then a test set is also created
                    new_data = (train = new_data[idxs[idx[1]:idx[2]], :], valid = new_data[idxs[idx[2]:idx[3]], :], test = new_data[idxs[idx[3]:idx[4]], :]);
                    new_data
                else
                    empty = DataFrame()
                    empty
                end
            else valid == "false"
                if test == "true"
                    idx_test = Vector{Int}([1, floor(2*dim/3), dim]);
                    new_data = (train = new_data[idxs[idx[1]:idx[3]], :], test = new_data[idxs[idx[3]:idx[4]], :]);
                    new_data
                else
                    new_test = MLJ.transform(std_mach, select(test_all, Not([:ZER_sunshine_1, :ALT_sunshine_4])))
                    new_data, new_test
                end
            end
        else
            if valid == "true"
                if test == "true"
                    #If valid is true then a test set is also created
                    new_data = (train = new_data[idxs[idx[1]:idx[2]], :], valid = new_data[idxs[idx[2]:idx[3]], :], test = new_data[idxs[idx[3]:idx[4]], :]);
                    new_data
                else
                    empty = DataFrame()
                    empty
                end
            else
                if test == "true"
                    idx_test = Vector{Int}([1, floor(2*dim/3), dim]);
                    new_data = (train = new_data[idxs[idx[1]:idx[3]], :], test = new_data[idxs[idx[3]:idx[4]], :]);
                    new_data
                else
                    new_data, test_all
                end
            end
        end
    end
end

#Commonly used losses for binary classification
function losses(machine, input, response)
    (loglikelihood = -sum(log_loss(predict(machine, input), response)),
     misclassification_rate = mean(predict_mode(machine, input) .!= response),
     accuracy = accuracy(predict_mode(machine, input), response),
     auc = MLJ.auc(predict(machine, input), response))
end

max_stars_repo_name: clara-rossignol/Will-it-rain-tomorrow

max_issues_repo_name: clara-rossignol/Will-it-rain-tomorrow

max_forks_repo_name: clara-rossignol/Will-it-rain-tomorrow


text:
using LinearAlgebra

using Test
using PrettyTables
using ProgressMeter
using UnicodePlots

using RollingFunctions

function test_derivatives(R::AbstractQuasiMatrix,
                          f::Function, g::Function, h::Function)
    r = axes(R,1)

    D = Derivative(r)

    ∇ = R' * D * R
    ∇² = R' * D' * D * R

    S = R'R
    nonortho = bandwidths(S) != (0,0)
    S⁻¹ = nonortho ? qr(S) : I

    fv = R \ f.(r)

    gv = similar(fv)
    mul!(gv, ∇, fv)
    nonortho && ldiv!(S⁻¹, gv)
    δgv = R \ g.(r) - gv

    hv = similar(gv)
    mul!(hv, ∇, gv)
    nonortho && ldiv!(S⁻¹, hv)
    δhv = R \ h.(r) - hv

    hv′ = similar(fv)
    mul!(hv′, ∇², fv)
    nonortho && ldiv!(S⁻¹, hv′)
    δhv′ = R \ h.(r) - hv′

    R,fv,gv,hv,hv′,δgv,δhv,δhv′
end

function estimate_convergence_rate(ρ, ϵ)
    i = findall(!iszero, ϵ)
    x = log10.(ρ[i])
    y = log10.(abs.(ϵ[i]))

    f = (x,y) -> ([x ones(length(x))] \ y)[1]

    # Find all slopes for consecutive groups of three errors
    slopes = rolling(f, x, y, 3)
    e = findfirst(s -> s < 0, slopes)
    if e !== nothing && e > 5
        slopes = slopes[1:max(1,e-1)]
    end

    maximum(slopes)
end

function error_slopes(hs, ϵ, names, verbosity=0)
    loghs = log10.(hs)
    logϵ = log10.(abs.(ϵ))
    n = size(logϵ,2)

    slopes = zeros(n)
    for j = 1:n
        slopes[j] = estimate_convergence_rate(hs, ϵ[:,j])
    end

    if verbosity > 0
        pretty_table([hs ϵ], vcat("h", ["$(name) [$(slope)]"
                                        for (name,slope) in zip(names,slopes)]))
        plt = lineplot(loghs, logϵ[:,1], name=names[1], xlabel="log10(h)", ylabel="log10(error)")
        for j = 2:n
            lineplot!(plt, loghs, logϵ[:,j], name=names[j])
        end
        # for (q,o) in enumerate(["linear", "quadratic", "cubic", "quartic", "quintic"])
        #     lineplot!(plt, loghs, q*loghs, name=o)
        # end
        display(plt)
        println()
        println()
    end

    slopes
end

function derivative_test_functions(d)
    a,b = √d*[-1,1]

    # The functions vanish at the boundaries; Dirichlet0 boundary
    # conditions. They /cannot/ be evaluated at the boundaries, since
    # that results in NaNs.
    f = x -> exp(-1/(d-x^2))
    g = x -> -2*exp(-1/(d-x^2))*x/((d-x^2)^2)
    h = x -> -2*exp(-1/(d-x^2))*(d^2 + 2*(d-1)*x^2-3x^4)/((d-x^2)^4)

    f,g,h,a,b
end

function compute_derivative_errors(fun::Function, Ns,
                                   f::Function, g::Function, h::Function,
                                   verbosity=0)
    p = Progress(length(Ns))
    errors = map(Ns) do N
        R = fun(N)
        R,fv,gv,hv,hv′,δgv,δhv,δhv′ = test_derivatives(R, f, g, h)
        S = R'R
        ProgressMeter.next!(p)
        [√(dot(δgv, S, δgv)) √(dot(δhv, S, δhv)) √(dot(δhv′, S, δhv′))]
    end |> e -> vcat(e...)

    ϵg = errors[:,1]
    ϵh = errors[:,2]
    ϵh′ = errors[:,3]

    hs = 1.0 ./ Ns
    pg,ph,ph′ = error_slopes(hs, errors, ["δg", "δh", "δh′"], verbosity)

    hs,ϵg,ϵh,ϵh′,pg,ph,ph′
end

function norm_rot!(v)
    normalize!(v)
    vc = v.args[2]
    vc[:] *= sign(vc[1])
    v
end

function diagonalize_hamiltonian(H, R::B, nev, σ; method=:arnoldi_shift_invert) where {B<:AbstractQuasiMatrix}
    A,target = if method == :arnoldi_shift_invert
        ShiftAndInvert(H, R, σ), LR()
    else
        H, SR()
    end
    schurQR,history = partialschur(A, nev=nev, which=target)

    ϕ = [norm_rot!(R*schurQR.Q[:,j]) for j = 1:nev]

    θ = schurQR.eigenvalues
    λ = if method == :arnoldi_shift_invert
        σ .+ inv.(θ)
    else
        θ
    end

    λ,ϕ
end

function get_kinetic_operator(R::B, Z=0, ℓ=0) where {B<:AbstractQuasiMatrix}
    D̃ = Derivative(Base.axes(R,1))
    D = !iszero(Z) ? CoulombDerivative(D̃, Z, ℓ) : D̃
    apply(*, R', D', D, R) / -2
end

function test_particle_in_a_box(R, L, nev; kwargs...) where {B<:AbstractQuasiMatrix}
    Tm = get_kinetic_operator(R)
    λ,ϕ = diagonalize_hamiltonian(Tm, R, nev, 0.0; kwargs...)

    r̃ = CompactBases.locs(R)

    n = 1:nev
    δλ = abs.(λ - n.^2*π^2/(2L^2))

    # Could/should also test eigenvectors

    λ,ϕ,r̃,R,δλ
end

function test_singular_scheme(R, Z, ℓ, nev; kwargs...) where {B<:AbstractQuasiMatrix}
    r = axes(R,1)

    n = 1:nev
    λₐ = -inv.(2(n .+ ℓ).^2)

    coulomb(r) = -inv(r)
    centrifugal(r) = ℓ*(ℓ+1)/2r^2

    Tm = get_kinetic_operator(R, Z, ℓ) + R'*QuasiDiagonal(centrifugal.(r))*R
    V = R'*QuasiDiagonal(coulomb.(r))*R
    H = Tm + V
    λ,ϕ = diagonalize_hamiltonian(H, R, nev, 1.1λₐ[1]; kwargs...)

    r̃ = CompactBases.locs(R)

    δλ = abs.(λ - λₐ)

    # Could/should also test eigenvectors

    λ,ϕ,r̃,R,δλ
end

function compute_diagonalization_errors(fun::Function, problem::Function, Ns, args...; verbosity=0, kwargs...) where B
    p = Progress(length(Ns))
    errors = map(Ns) do N
        t = time()
        λ,ϕ,r̃,R,δλ = problem(fun(N), args...; kwargs...)
        elapsed = time()-t
        ProgressMeter.next!(p)
        vcat(δλ,elapsed)'
    end |> e -> vcat(e...)

    errors,elapsed = errors[:,1:end-1],errors[:,end]

    hs = 1.0 ./ Ns
    slopes = error_slopes(hs, errors, ["$j" for j in 1:size(errors,2)], verbosity)

    errors,slopes,elapsed
end


text:
#Code written by Skipper af Rosenborg and Esben Bager
##### MEAN-VARIANCE OPTIMIZATION
#This is not supposed to produce great results, since mean and covariance
#is supposed to be very exact before good weights are found

#Load packages and support functions
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
include("SupportFunction.jl")
include("DataLoad.jl")
include("CEOSupportFunctions.jl")
println("Leeeeroooy Jenkins")

#Set path for dataloading
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/Data/IndexDataDiff/"
path = "/zhome/9f/d/88706/SpecialeCode/Thesis/Data/IndexDataDiff/"

#Set parameters
trainingSize = 240
possibilities = 5
industries = ["NoDur", "Durbl", "Manuf", "Enrgy", "HiTec", "Telcm", "Shops", "Hlth", "Utils", "Other"]
industriesTotal = length(industries)

#Define the input and transformations
modelMatrix = zeros(industriesTotal, possibilities)
noDurModel = [1 0 0 1 1]
testModel = [1 0 0 1 1]
modelMatrix[1, :] = noDurModel
for i=2:industriesTotal
    modelMatrix[i, :] = noDurModel
end

#Create and fill X and Y arrays
XArrays = Array{Array{Float64, 2}}(industriesTotal)
YArrays = Array{Array{Float64, 2}}(industriesTotal)
XArrays, YArrays = generateXandYs(industries, modelMatrix)

rfRates = loadRiskFreeRate("NoDur", path)
rfRates = rfRates[:,1]

#Set parameter
nGammas = 5
standY = YArrays[1]
nRows = size(standY)[1]
amountOfModels = nGammas^4

#Initialization of parameters
w1N = repeat([1/11], outer = 11) #1/N weights
gamma = 2 #risk aversion
validationPeriod = 0 # Number of periods to do parameter training
PMatrix = zeros(nRows-trainingSize, amountOfModels)

bestModelAmount = 5
bestModelConfigs = zeros(bestModelAmount, 4)
bestModelIndexes = zeros(bestModelAmount)

return1NMatrix = zeros(nRows-trainingSize)
returnCEOMatrix = zeros(nRows-trainingSize, bestModelAmount)
returnPerfectMatrix = zeros(nRows-trainingSize)

weightsPerfect = zeros(nRows-trainingSize, 11)
weightsCEO     = zeros(nRows-trainingSize, 11, bestModelAmount)
expectedReturnMatrix = zeros(nRows-trainingSize, 11)
forecastErrors = zeros(nRows-trainingSize, 11)

#Load perfect results
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/RFR/VIXTimeTA2.4/"
path = "/zhome/9f/d/88706/SpecialeCode/Results/CEORFR/"
weightsPerfect = Array{Float64}(CSV.read(path*"weightsPerfect.csv",header=false, datarow=1, nullable=false))
returnPerfectMatrix = Array{Float64}(CSV.read(path*"returnPerfectMatrix.csv", header=false, datarow=1, nullable=false))

#Set models to investigate
bestModelConfigs[1,:] = [1 0.5  0.25  1]
bestModelConfigs[2,:] = [1 0.5  0.25  2.5]
bestModelConfigs[3,:] = [1 0.5  0.25  10]
bestModelConfigs[4,:] = [1 0.5  0.25  50]
bestModelConfigs[5,:] = [1 0.5  0.25  100]

for i = 1:bestModelAmount
    bestModelIndexes[i] = i
end

inputArg1 = 0 #83 breaks the trainingSize
inputArg1 = parse(Int64,ARGS[1])
if 10+validationPeriod+(10*inputArg1) <= nRows-trainingSize-2
    println("Starting CEO Validation loop")
    for t=validationPeriod+1+(10*inputArg1):10+validationPeriod+(10*inputArg1)#(nRows-trainingSize-2)
        println("Time $t/50")
        trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY = createDataSplits(XArrays, YArrays, t, trainingSize)

        for m = 1:bestModelAmount
            expectedReturns = zeros(11)

            rfRatesVec = rfRates[t:(t+trainingSize-1)]

            #Build and solve CEO
            betaArray, U = @time(runCEORFR(trainingXArrays, trainingYArrays, bestModelConfigs[m, :], gamma, rfRatesVec))

            # Forecast next period
            expectedReturns[1:10] = generateExpectedReturns(betaArray, trainingXArrays, trainingYArrays, validationXRows)
            expectedReturns[11] = rfRates[t+trainingSize]
            if m == 1
                expectedReturnMatrix[t, 1:10] = (exp.(expectedReturns[1:10])-1)*100
                expectedReturnMatrix[t, 11] = (exp.(expectedReturns[11])-1)
            end

            valY = zeros(11)
            for i = 1:10
                valY[i] = validationY[i][1]
            end
            valY[11] = rfRates[t+trainingSize]

            #Set Sigma for MV
            rfRatesVec = rfRates[t:(t+trainingSize-1)]
            trainX = hcat(trainingXArrays[1][:,1:10], rfRatesVec)
            Sigma =  cov(trainX)
            F = lufact(Sigma)
            U = F[:U]  #Cholesky factorization of Sigma

            return1N, returnCEO, wStar, forecastRow = performMVOptimizationRISK(expectedReturns, U, gamma, valY, valY)
            weightsCEO[t, 1:11, Int64(m)] = wStar
            return1NMatrix[t]      = return1N
            returnCEOMatrix[t, Int64(bestModelIndexes[m])] = returnCEO
            returnPerfect = returnPerfectMatrix[t]
            if m == 1
                forecastErrors[t, 1:11]  = forecastRow-expectedReturnMatrix[t,:]
            end
            println("1N returns is $return1N, returnPerfect is $returnPerfect and returnCEO is $returnCEO")
            PMatrix[t, Int64(bestModelIndexes[m])] = calculatePvalue(return1N, returnPerfect, returnCEO)
        end
    end
elseif validationPeriod+(10*inputArg1) <= nRows-trainingSize-2
    for t=validationPeriod+1+(10*inputArg1):(nRows-trainingSize-2)
        println("Time $t/50")
        trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY = createDataSplits(XArrays, YArrays, t, trainingSize)

        for m = 1:bestModelAmount
            expectedReturns = zeros(11)
            rfRatesVec = rfRates[t:(t+trainingSize-1)]
            betaArray, U = @time(runCEORFR(trainingXArrays, trainingYArrays, bestModelConfigs[m, :], gamma, rfRatesVec))

            expectedReturns[1:10] = generateExpectedReturns(betaArray, trainingXArrays, trainingYArrays, validationXRows)
            expectedReturns[11] = rfRates[t+trainingSize]
            if m == 1
                expectedReturnMatrix[t, 1:10] = (exp.(expectedReturns[1:10])-1)*100
                expectedReturnMatrix[t, 11] = (exp.(expectedReturns[11])-1)
            end

            valY = zeros(11)
            for i = 1:10
                valY[i] = validationY[i][1]
            end
            valY[11] = rfRates[t+trainingSize]

            rfRatesVec = rfRates[t:(t+trainingSize-1)]
            trainX = hcat(trainingXArrays[1][:,1:10], rfRatesVec)
            Sigma =  cov(trainX)
            F = lufact(Sigma)
            U = F[:U]  #Cholesky factorization of Sigma

            return1N, returnCEO, wStar, forecastRow = performMVOptimizationRISK(expectedReturns, U, gamma, valY, valY)
            weightsCEO[t, 1:11, Int64(m)] = wStar
            return1NMatrix[t]      = return1N
            returnCEOMatrix[t, Int64(bestModelIndexes[m])]  = returnCEO
            returnPerfect = returnPerfectMatrix[t]
            if m == 1
                forecastErrors[t, 1:11]  = forecastRow-expectedReturnMatrix[t,:]
            end
            println("1N returns is $return1N, returnPerfect is $returnPerfect and returnCEO is $returnCEO")
            PMatrix[t, Int64(bestModelIndexes[m])] = calculatePvalue(return1N, returnPerfect, returnCEO)
        end
    end
end

#Output data
combinedPortfolios = hcat(returnPerfectMatrix[1:nRows-trainingSize-2, 1], return1NMatrix[1:nRows-trainingSize-2, 1],
    returnCEOMatrix[1:nRows-trainingSize-2, Array{Int64}(bestModelIndexes)], PMatrix[1:nRows-trainingSize-2, Array{Int64}(bestModelIndexes)])
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/"
path = "/zhome/9f/d/88706/SpecialeCode/Results/CEORFR/"
writedlm(path*string(inputArg1)*"_returnPvalueOutcome1to200.csv", combinedPortfolios, ",")
for i = 1:bestModelAmount
    writedlm(path*"wightsCEO_Model"*string(i)*"_"*string(inputArg1)*".csv",weightsCEO[:,:,i], ",")
end

combinedPortfolios2 = hcat(expectedReturnMatrix, forecastErrors)
writedlm(path*string(inputArg1)*"_forecastsAndErrors1to200.csv", combinedPortfolios2, ",")
println("Finished Everything")

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_path: ML/Code to hand in/CEORFR.jl

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_path: ML/Code to hand in/CEORFR.jl

max_forks_repo_name: skipperafrosenborg/Thesis


text:
# # Mesh grid for mode finding - Part 1
#
# The [Global complex Roots and Poles Finding](https://doi.org/10.1109/TAP.2018.2869213)
# (GRPF) algorithm searches for roots and poles of complex-valued functions by sampling
# the function at the nodes of a regular mesh and analyzing the function's complex phase.
# The mesh search region can be arbitrarily shaped and guarantees that none of
# the roots/poles can be missed as long as the mesh step is sufficiently small.
# Specifically, the phase change of the function value cannot exceed three quadrants
# across mesh edges. Effectively this means there cannot be equal numbers of roots and
# poles within the boundary of each identified candidate region.
#
# Because the GRPF algorithm uses an adaptive meshing process, the initial mesh
# node spacing can actually be greater than the distance between the roots/poles
# of the function, but there is no established method for _a priori_ estimation
# of the required initial sampling of an arbitrary function.
# The GRPF algorithm is implemented in
# [RootsAndPoles.jl](https://github.com/fgasdia/RootsAndPoles.jl) and used by
# LongwaveModePropagator to identify eigenangles of [`HomogeneousWaveguide`](@ref)'s.
# Therefore, we must experimentally determine a sufficient initial mesh grid spacing to
# identify roots of the mode equation.
#
# ## Mode equation
#
# The criteria for mode propagation in the Earth-ionosphere waveguide is described by the
# fundamental equation of mode theory
#
# ```math
# \det(\overline{\bm{R}}(\theta)\bm{R}(\theta) - \bm{I}) = 0
# ```
#
# where ``\overline{\bm{R}}(\theta)`` is the reflection coefficient matrix for the
# ground and ``\bm{R}(\theta)`` is the reflection coefficient of the ionosphere.
# Both are functions of the complex angle ``\theta``. The discrete values of
# ``\theta`` for which the fundamental equation is satisfied are known as
# _eigenangles_.
#
# First, import the necessary packages.

using Plots
using Plots.Measures
using RootsAndPoles

using LongwaveModePropagator
using LongwaveModePropagator: QE, ME, solvemodalequation, trianglemesh

# ## Exploring the mode equation phase
#
# The U.S. Navy's Long Wavelength Propagation Capability (LWPC) searches the region of the
# complex plane from 30° to 90° in the real axis and 0° to -10° in the
# imaginary axis. The lowest order modes are nearest 90° - i0°, excluding 90°.
# Modes with large negative imaginary components are highly attenuated and have
# relatively little affect on the total field in the waveguide.
#
# Let's begin by exploring the modal equation phase on a fine grid across the bottom
# right complex quadrant.
#
# We'll define four different [`PhysicalModeEquation`](@ref)'s that we'll use
# throughout this example.

lowfrequency = Frequency(10e3)
midfrequency = Frequency(20e3)
highfrequency = Frequency(100e3)

day = Species(QE, ME, z->waitprofile(z, 75, 0.35), electroncollisionfrequency)
night = Species(QE, ME, z->waitprofile(z, 85, 0.9), electroncollisionfrequency)

daywaveguide = HomogeneousWaveguide(BField(50e-6, π/2, 0), day, GROUND[5])
nightwaveguide = HomogeneousWaveguide(BField(50e-6, π/2, 0), night, GROUND[5])

day_low_me = PhysicalModeEquation(lowfrequency, daywaveguide)
day_mid_me = PhysicalModeEquation(midfrequency, daywaveguide)
day_high_me = PhysicalModeEquation(highfrequency, daywaveguide)

night_mid_me = PhysicalModeEquation(midfrequency, nightwaveguide)

day_low_title = "10 kHz\nh′: 75, β: 0.35"
day_mid_title = "20 kHz\nh′: 75, β: 0.35"
day_high_title = "100 kHz\nh′: 75, β: 0.35"
night_mid_title = "20 kHz\nh′: 85, β: 0.9"
nothing  #hide

# Here we define a dense rectangular mesh.

Δr = 0.2
x = 0:Δr:90
y = -40:Δr:0
mesh = x .+ im*y';

# Now we simply iterate over each node of the mesh, evaluating the modal equation with
# [`LongwaveModePropagator.solvemodalequation`](@ref), explicitly imported from
# LongwaveModePropagator above.
# We use Julia's `Threads.@threads` multithreading capability to speed up the computation.

function modeequationphase(me, mesh)
    phase = Vector{Float64}(undef, length(mesh))
    Threads.@threads for i in eachindex(mesh)
        f = solvemodalequation(deg2rad(mesh[i]), me)
        phase[i] = rad2deg(angle(f))
    end
    return phase
end

phase = modeequationphase(day_mid_me, mesh);

# If an instability occurs in the integration of the reflection coefficient,
# it likely happened at the angle 90° + i0°, which is not a valid mode.
# This isn't a problem for plotting though!
#
# Plotting the results, we can see that there are clearly identifiable locations where
# white, black, blue, and orange, each representing a different quadrant of the
# complex plane, all meet.
# Each of these locations are either a root or pole in the daytime ionosphere.

heatmap(x, y, reshape(phase, length(x), length(y))';
        color=:twilight, clims=(-180, 180),
        xlims=(0, 90), ylims=(-40, 0),
        xlabel="real(θ)", ylabel="imag(θ)",
        title=day_mid_title,
        right_margin=2mm)
#md savefig("meshgrid_20kday.png"); nothing # hide
#md # ![](meshgrid_20kday.png)

# We can zoom in to the upper right corner of the plot to see the lowest order modes:

heatmap(x, y, reshape(phase, length(x), length(y))';
        color=:twilight, clims=(-180, 180),
        xlims=(30, 90), ylims=(-10, 0),
        xlabel="real(θ)", ylabel="imag(θ)",
        title=day_mid_title,
        right_margin=2mm)
#md savefig("meshgrid_20kdayzoom.png"); nothing # hide
#md # ![](meshgrid_20kdayzoom.png)

# If we switch to a nighttime ionosphere with a high Wait β parameter, we see that the
# roots move closer to the axes.
# A perfectly reflecting conductor has eigenangles along the real and complex axes.

phase = modeequationphase(night_mid_me, mesh);

heatmap(x, y, reshape(phase, length(x), length(y))';
        color=:twilight, clims=(-180, 180),
        xlims=(0, 90), ylims=(-40, 0),
        xlabel="real(θ)", ylabel="imag(θ)",
        title=night_mid_title,
        right_margin=2mm)
#md savefig("meshgrid_20knight.png"); nothing # hide
#md # ![](meshgrid_20knight.png)

# At lower frequencies, the roots/poles move further apart.

phase = modeequationphase(day_low_me, mesh);

heatmap(x, y, reshape(phase, length(x), length(y))';
        color=:twilight, clims=(-180, 180),
        xlims=(0, 90), ylims=(-40, 0),
        xlabel="real(θ)", ylabel="imag(θ)",
        title=day_low_title,
        right_margin=2mm)
#md savefig("meshgrid_10kday.png"); nothing # hide
#md # ![](meshgrid_10kday.png)

# ## Global complex roots and poles finding
#
# The global complex roots and poles finding (GRPF) algorithm is most efficient when
# the initial mesh grid consists of equilateral triangles.
# Such a grid can be produced with the `rectangulardomain` function from
# RootsAndPoles, but as seen in the evaluation of the modal equation above, no
# roots or poles appear in the lower right diagonal half of the domain.
# Even if they did, they would correspond to highly attenuated modes.
# Therefore, to save compute time, we can exclude the lower right triangle of the domain
# from the initial mesh.
#
# The function [`LongwaveModePropagator.trianglemesh`](@ref) is built into
# LongwaveModePropagator for this purpose.
# The inputs are specified by the complex bottom left corner `zbl`, the top right corner
# `ztr`, and the mesh spacing `Δr` in _radians_.

zbl = deg2rad(complex(30.0, -10.0))
ztr = deg2rad(complex(89.9, 0.0))
Δr = deg2rad(0.5)

mesh = trianglemesh(zbl, ztr, Δr);

# We convert back to degrees just for plotting.
# Here's a zoomed in portion of the upper right of the domain.

meshdeg = rad2deg.(mesh)

img = plot(real(meshdeg), imag(meshdeg); seriestype=:scatter,
           xlims=(80, 90), ylims=(-10, 0),
           xlabel="real(θ)", ylabel="imag(θ)",
           legend=false, size=(450,375));
plot!(img, [80, 90], [0, 0]; color="red");
plot!(img, [0, 90], [-90, 0]; color="red")
#md savefig(img, "meshgrid_trianglemesh.png"); nothing # hide
#md # ![](meshgrid_trianglemesh.png)

#
# Now let's apply `grpf` to the modal equation on the triangle mesh.
# `grpf` adaptively refines the mesh to obtain a more accurate estimate of the position
# of the roots and poles.
#
# [`LMPParams`](@ref) is a struct use for passing parameters across LongwaveModePropagator.
# Default values are automatically inserted when instantiating the struct.
# We'll extract the `grpfparams` field which contains a `GRPFParams` struct for passing
# to `grpf`.
#
# We will also pass RootsAndPoles' `PlotData()` argument to obtain additional
# information on the function phase for plotting.

params = LMPParams().grpfparams
roots, poles, quads, phasediffs, tess, g2f = grpf(θ->solvemodalequation(θ, day_mid_me),
                                                  mesh, PlotData(), params);

# The `getplotdata` function, provided by RootsAndPoles, is a convenience function
# for plotting a color-coded tesselation from `grpf` based on the phase of the function.

z, edgecolors = getplotdata(tess, quads, phasediffs, g2f)

rootsdeg = rad2deg.(roots)
polesdeg = rad2deg.(poles)
zdeg = rad2deg.(z)

twilightquads = [
    colorant"#9E3D36",
    colorant"#C99478",
    colorant"#7599C2",
    colorant"#5C389E",
    colorant"#404040",
    RGB(0.0, 0.0, 0.0)
]

img = plot(real(zdeg), imag(zdeg); group=edgecolors, palette=twilightquads, linewidth=1.5,
           xlims=(30, 90), ylims=(-10, 0),
           xlabel="real(θ)", ylabel="imag(θ)", legend=false,
           title=day_mid_title);
plot!(img, real(rootsdeg), imag(rootsdeg); color="red",
      seriestype=:scatter, markersize=5);
plot!(img, real(polesdeg), imag(polesdeg); color="red",
      seriestype=:scatter, markershape=:utriangle, markersize=5)
#md savefig(img, "meshgrid_20kdaymesh.png"); nothing # hide
#md # ![](meshgrid_20kdaymesh.png)

# In the plot above, roots are marked with red circles and poles are marked
# with red triangles.
# The automatic refinement of the mesh is clearly visible.
#
# Here are similar plots for the other three scenarios.

## Daytime ionosphere, low frequency
roots, poles, quads, phasediffs, tess, g2f = grpf(θ->solvemodalequation(θ, day_low_me),
                                                  mesh, PlotData(), params);
z, edgecolors = getplotdata(tess, quads, phasediffs, g2f)

rootsdeg = rad2deg.(roots)
polesdeg = rad2deg.(poles)
zdeg = rad2deg.(z)

img = plot(real(zdeg), imag(zdeg); group=edgecolors, palette=twilightquads, linewidth=1.5,
           xlims=(30, 90), ylims=(-10, 0),
           xlabel="real(θ)", ylabel="imag(θ)", legend=false,
           title=day_low_title);
plot!(img, real(rootsdeg), imag(rootsdeg); color="red",
      seriestype=:scatter, markersize=5);
plot!(img, real(polesdeg), imag(polesdeg); color="red",
      seriestype=:scatter, markershape=:utriangle, markersize=5)
#md savefig(img, "meshgrid_10kdaymesh.png"); nothing # hide
#md # ![](meshgrid_10kdaymesh.png)

# At 100 kHz, `grpf` requires more mesh refinements and takes considerably
# more time to run.

## Daytime ionosphere, high frequency
roots, poles, quads, phasediffs, tess, g2f = grpf(θ->solvemodalequation(θ, day_high_me),
                                                  mesh, PlotData(), params);
z, edgecolors = getplotdata(tess, quads, phasediffs, g2f)

rootsdeg = rad2deg.(roots)
polesdeg = rad2deg.(poles)
zdeg = rad2deg.(z)

img = plot(real(zdeg), imag(zdeg); group=edgecolors, palette=twilightquads, linewidth=1.5,
           xlims=(30, 90), ylims=(-10, 0),
           xlabel="real(θ)", ylabel="imag(θ)", legend=false,
           title=day_high_title);
plot!(img, real(rootsdeg), imag(rootsdeg); color="red",
      seriestype=:scatter, markersize=5);
plot!(img, real(polesdeg), imag(polesdeg); color="red",
      seriestype=:scatter, markershape=:utriangle, markersize=5)
#md savefig(img, "meshgrid_100kdaymesh.png"); nothing # hide
#md # ![](meshgrid_100kdaymesh.png)

#
## Nighttime ionosphere, 20 kHz
roots, poles, quads, phasediffs, tess, g2f = grpf(θ->solvemodalequation(θ, night_mid_me),
                                                  mesh, PlotData(), params);
z, edgecolors = getplotdata(tess, quads, phasediffs, g2f)

rootsdeg = rad2deg.(roots)
polesdeg = rad2deg.(poles)
zdeg = rad2deg.(z)

img = plot(real(zdeg), imag(zdeg); group=edgecolors, palette=twilightquads, linewidth=1.5,
           xlims=(30, 90), ylims=(-10, 0),
           xlabel="real(θ)", ylabel="imag(θ)", legend=false,
           title=night_mid_title);
plot!(img, real(rootsdeg), imag(rootsdeg); color="red",
      seriestype=:scatter, markersize=5);
plot!(img, real(polesdeg), imag(polesdeg); color="red",
      seriestype=:scatter, markershape=:utriangle, markersize=5)
#md savefig(img, "meshgrid_20knightmesh.png"); nothing # hide
#md # ![](meshgrid_20knightmesh.png)

# The example continues in [Mesh grid for mode finding - Part 2](@ref).


text:
module SetIntersectionProjection

#packages that used to be in base in Julia 0.6
using Distributed
using LinearAlgebra
using SparseArrays
using Printf
using Statistics

#other packages
using Parameters
using Interpolations
using DistributedArrays
using JOLI
using JOLI.FFTW, JOLI.Wavelets
using SortingAlgorithms
using TimerOutputs
# using Flux
# using NNlib
# using CUDA

export log_type_PARSDMM, set_properties, PARSDMM_options, set_definitions

#main scripts
include("PARSDMM.jl")
include("PARSDMM_multi_level.jl")
include("PARSDMM_precompute_distribute.jl")

#algorithms and solver scripts
include("argmin_x.jl");
include("adapt_rho_gamma.jl");
include("update_y_l.jl");
include("stop_PARSDMM.jl");
include("PARSDMM_initialize.jl")
include("Q_update!.jl")
include("rhs_compose.jl")

#linear algebra functions
include("cg.jl");
include("mat2CDS.jl")
include("CDS_MVp.jl")
include("CDS_MVp_MT.jl")
include("CDS_MVp_MT_subfunc.jl")
include("CDS_scaled_add!.jl")
#include("CDS2stencil.jl")

#scripts for parallelism
include("update_y_l_parallel.jl")
include("adapt_rho_gamma_parallel.jl")
include("compute_relative_feasibility.jl")

#multi-threaded scripts
include("a_is_b_min_c_MT!.jl")

#multi level related scripts
include("setup_multi_level_PARSDMM.jl")
include("constraint2coarse.jl")
include("interpolate_y_l.jl")

#scripts for setting up constraints, projectors, linear operators
include("default_PARSDMM_options.jl")
include("convert_options!.jl")
include("get_discrete_Grad.jl");
include("get_TD_operator.jl");
include("get_projector.jl");
include("get_bound_constraints.jl");
include("setup_constraints.jl");

#projectors
include("projectors/project_bounds!.jl")
include("projectors/project_cardinality!.jl");
include("projectors/project_rank!.jl");
include("projectors/project_nuclear!.jl");
include("projectors/project_l1_Duchi!.jl");
include("projectors/project_l2!.jl");
include("projectors/project_annulus!.jl");
include("projectors/project_subspace!.jl");
include("projectors/project_histogram_relaxed.jl");

#other proximal maps
include("prox_l2s!.jl")
include("prox_l1!.jl")

#scripts that are required to run examples
include("constraint_learning_by_observation.jl")

#Minkowski sum related scripts
include("PARSDMM_precompute_distribute_Minkowski.jl")

#define types
mutable struct log_type_PARSDMM
      set_feasibility   :: Array{Real,2}
      r_dual            :: Array{Real,2}
      r_pri             :: Array{Real,2}
      r_dual_total      :: Vector{Real}
      r_pri_total       :: Vector{Real}
      obj               :: Vector{Real}
      evol_x            :: Vector{Real}
      rho               :: Array{Real,2}
      gamma             :: Array{Real,2}
      cg_it             :: Vector{Integer}
      cg_relres         :: Vector{Real}
      timing            :: TimerOutput
end

@with_kw mutable struct PARSDMM_options
  x_min_solver          :: String       = "CG_normal" #what algorithm to use for the x-minimization (CG applied to normal equations)
  maxit                 :: Integer      = 200         #max number of PARSDMM iterations
  evol_rel_tol          :: Real         = 1e-3        #stop PARSDMM if ||x^k - X^{k-1}||_2 / || x^k || < options.evol_rel_tol AND options.feas_tol is reached
  feas_tol              :: Real         = 5e-2        #stop PARSDMM if the transform-domain relative feasibility error is < options.feas_tol AND options.evol_rel_tol is reached
  obj_tol               :: Real         = 1e-3        #optional stopping criterion for change in distance from point that we want to project
  rho_ini               :: Vector{Real} = [10.0]      #initial values for the augmented-Lagrangian penalty parameters. One value in array or one value per constraint set in array
  rho_update_frequency  :: Integer      = 2           #update augmented-Lagrangian penalty parameters and relaxation parameters every X number of PARSDMM iterations
  gamma_ini             :: Real         = 1.0         #initial value for all relaxation parameters (scalar)
  adjust_rho            :: Bool         = true        #adapt augmented-Lagrangian penalty parameters or not
  adjust_gamma          :: Bool         = true        #adapt relaxation parameters in PARSDMM
  adjust_feasibility_rho:: Bool         = true        #adapt augmented-Lagrangian penalty parameters based on constraint set feasibility errors (can be used in combination with options.adjust_rho)
  Blas_active           :: Bool         = true        #use direct BLAS calls, otherwise the code will use Julia loop-fusion where possible
  feasibility_only      :: Bool         = false       #drop distance term and solve a feasibility problem
  FL                    :: DataType     = Float32     #type of Float: Float32 or Float64
  parallel              :: Bool         = false       #compute proximal mappings, multiplier updates, rho and gamma updates in parallel
  zero_ini_guess        :: Bool         = true        #zero initial guess for primal, auxiliary, and multipliers
  Minkowski             :: Bool         = false       #the intersection of sets includes a Minkowski set
end

#save properties of the constraint set and its linear operator (not all are currently used in the code)
#each entry in the following vectors contain the information corresponding to one set
mutable struct set_properties
           ncvx       ::Vector{Bool}                                          #is the set non-convex? (true/false)
           AtA_diag   ::Vector{Bool}                                          #is A^T A a diagonal matrix?
           dense      ::Vector{Bool}                                          #is A a dense matrix?
           TD_n       ::Vector{Tuple}                                         #the grid dimensions in the transform-domain.
           tag        ::Vector{Tuple{String,String,String,String}}            #(constraint type, linear operator description, application mode,application direction if applicable)
           banded     ::Vector{Bool}                                          #is A a banded matrix?
           AtA_offsets::Union{Vector{Vector{Int32}},Vector{Vector{Int64}}}    #only required if A is banded. A vector of indices of the non-zero diagonals, where the main diagonal is index 0
end

mutable struct set_definitions
  set_type          ::String #"bounds","l1","l2","annulus","subspace","nuclear","rank","histogram","cardinality"
  TD_OP             ::String
  min               ::Union{Vector,Real}
  max               ::Union{Vector,Real}
  app_mode          ::Tuple{String,String} #("tensor/"matrix"/"slice"/"fiber" , "x","y","z") , x,y,z apply only to slice and fibers of a tensor or matrix
  custom_TD_OP      ::Tuple{Union{Array{Any},SparseMatrixCSC{Float32,Int64},SparseMatrixCSC{Float64,Int64},Array{Float64,2},Array{Float32,2}, joAbstractLinearOperator},Bool} #custom matrix for subspace, and boolean to indicate orthogonality
end

end # module

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
using CellularAutomata
using Test

# Function to test Conway's blinker oscillator and toad oscillator
function blinker_test()
    init = Array{Int}(zeros(10,10))

    # 'Toad Oscillator'
    init[6, 4] = 1
    init[6, 5] = 1
    init[6, 6] = 1

    init[7, 5] = 1
    init[7, 6] = 1
    init[7, 7] = 1

    ca = CA2d([3], [2,3], init, 9);

    @test ca.cells[:, :, 3] == ca.cells[:, :, 1]
    @test ca.cells[:, :, 5] == ca.cells[:, :, 3]
    @test ca.cells[:, :, 4] == ca.cells[:, :, 2]
    @test ca.cells[:, :, 8] == ca.cells[:, :, 6]

    # 'Blinker Oscillator'
    init[6, 4] = 0
    init[6, 5] = 0
    init[6, 6] = 0

    ca = CA2d([3], [2,3], init, 9)

    @test ca.cells[:, :, 3] == ca.cells[:, :, 1]
    @test ca.cells[:, :, 5] == ca.cells[:, :, 3]
    @test ca.cells[:, :, 4] == ca.cells[:, :, 2]
    @test ca.cells[:, :, 8] == ca.cells[:, :, 6]
end

# Function to test Conway's turbine automota which repeats every ninth step
function turbine_test()
    init = zeros(Int, 20, 20)
    init[4, 4:9] .= 1
    init[5, 4:9] .= 1

    init[4:9, 11] .= 1
    init[4:9, 12] .= 1

    init[11, 7:12] .= 1
    init[12, 7:12] .= 1

    init[7:12, 4] .= 1
    init[7:12, 5] .= 1

    ca = CA2d([3], [2,3], init, 18)

    for i in 9:18
        @test ca.cells[:, :, i] == ca.cells[:, :, i-8]
    end
end

turbine_test()

blinker_test()


text:
# This file contains the definition for the Lines struct.

"""
calc_lineparams(lines_df::DataFrame)
	Calculates line parameters with regard to voltage level, number of circuits,
	and Thermal Reliability Margin (TRM).
	Returns reactance (p.u.), resistance (p.u.), and pmax (MW)
"""
function calc_lineparams(lines_df::DataFrame, Sbase::Number, TRM::Number)

	zbase(voltage::Number) = (voltage*1E3)^2 / (Sbase*1E6)

	# L = Symbol.(lines_df[!, 1])
	L = lines_df[!, 1]
	zbase_vlevel = map(zbase, lines_df[!, :voltage])	# Calculate reactance in p.u. given voltage level

	reactance = Dict(zip(L, lines_df[!, :reactance] ./ lines_df[!, :circuits] ./ zbase_vlevel))
	resistance = Dict(zip(L,  lines_df[!, :resistance] ./ lines_df[!, :circuits] ./ zbase_vlevel))

	pmax = Dict(zip(L, lines_df[!, :pmax] .* lines_df[!, :circuits] .* (1-TRM) ))
	return reactance, resistance, pmax
end

mutable struct Lines
    id
    from
    to
	voltage
	reactance
	resistance
	pmax
	circuits

    function Lines(lines_df::DataFrame, Sbase::Number, TRM::Number)

		# L = Symbol.(lines_df[!, 1])
		L = lines_df[!, 1]
		reactance, resistance, pmax = calc_lineparams(lines_df, Sbase, TRM)                    # U = unit

		function lines_df_to_dict_with_id(df::DataFrame)
		    Dict(col[1] => Dict(zip(df[!, 1], col[2])) for col in eachcol(df, true))
		end

		lines_dict = lines_df_to_dict_with_id(lines_df)

        return new(L,
                   lines_dict[:from],
                   lines_dict[:to],
				   lines_dict[:voltage],
				   reactance,
				   resistance,
				   pmax,
				   lines_dict[:circuits]
				   )
    end
end

"""
    create_incidence(lines::Lines, nodes::Nodes)
Create the incidence matrix from 'lines' and 'nodes'.
"""
function create_incidence(lines::Lines, nodes::Nodes)
	incidence = zeros(Int, length(lines.id), length(nodes.id))
	for (i,l) in enumerate(lines.id)
	    incidence[i, findfirst(nodes.id .== lines.from[l])] = 1
	    incidence[i, findfirst(nodes.id .== lines.to[l])] = -1
	end
	return incidence
end

"""
    calc_h_b(lines::Lines, nodes::Nodes)
Calculate the H and B matrices from 'lines' and 'nodes'.
"""
function calc_h_b(lines::Lines, nodes::Nodes)

	incidence = create_incidence(lines, nodes)
	bvector = [lines.reactance[l] /
		(lines.reactance[l].^2 .+ lines.resistance[l].^2) for l in lines.id]
	h_matrix = bvector .* incidence
	b_matrix = h_matrix' *incidence

	h = Dict((l,n) => h_matrix[i,j]
		for (i,l) in enumerate(lines.id), (j,n) in enumerate(nodes.id))

	b = Dict((n,nn) => b_matrix[i,j]
		for (i,n) in enumerate(nodes.id), (j,nn) in enumerate(nodes.id))

	return h,b
end

max_issues_repo_name: bobbyxiong/redispatch-ptg

max_forks_repo_name: bobbyxiong/redispatch-ptg


text:
function crop(img::AbstractImage, x::Range, y::Range)
    @assert ndims(img) == 2
    xidx = dimindex(img, "x")
    ndata = if xidx == 1
        w, h = size(img)
        Base.getindex(img, x, y)
    else
        h, w = size(img)
        Base.getindex(img, y, x)
    end
    copyproperties(img, ndata)
end


text:
using StaticArrays

struct Polynomial{N,T}
    data::SVector{N,T}
end
Polynomial(x::T...) where {T} = Polynomial{length(x),T}(SVector{length(x),T}(x))
Polynomial(a::T) where {T<:Number} = Polynomial{1,T}(SVector{1,T}(a))

Base.copy(p::Polynomial) = p
Base.one(p::Polynomial{N,T}) where{N,T} = Polynomial(@SVector[one(T)])

Base.length(p::Polynomial) = length(p.data)
Base.eltype(p::Polynomial) = eltype(p.data)
Base.getindex(p::Polynomial,i::Int) = p.data[i+1]
Base.setindex!(p::Polynomial,v,i::Int) = (p.data[i+1] = v)

degree(p::Polynomial) = length(p.data)-1

function derive(p)
    x = eltype(p)[d*p[d] for d in 1:degree(p)]
    Polynomial(x...)
end

function integrate(p::Polynomial{N,T},x0,y0) where {N,T}

    coeffs = similar(p.data, length(p.data)+1)
    fill!(coeffs, 0)
    for i in 2:length(coeffs)
        coeffs[i] = p.data[i-1]/(i-1)
    end
    temp = Polynomial(SVector(coeffs...))
    coeffs[1] = y0 - temp(x0)
    # TODO: Remove high order terms with zero coefficient
    return Polynomial(SVector(coeffs...))
end

function evaluate(p::Polynomial, t)
    d = p.data
    T = promote_type(typeof(t), eltype(d))
    r = zero(T)
    for i in 1:length(d)
        r += d[i] * t^(i-1)
    end
    return r
end

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

function integrate(p)
    x = zeros(eltype(p), degree(p)+2)
    for i in 1:degree(p)+1
        x[i+1] = p[i-1]/i
    end
    Polynomial(x...)
end

import Base: *,+,-,promote_rule

function *(p::Polynomial,q::Polynomial)
    T = promote_type(eltype(p), eltype(q))
    D = degree(p)+degree(q)+1
    x = zeros(T, D)
    for d in 0 : D
        for k in 0 : D
            k <= degree(p) || continue
            0 <= d-k <= degree(q) || continue
            x[d+1] += p[k] * q[d-k]
        end
    end
    Polynomial(x...)
end

# function *(p::Polynomial{1}, q::Polynomial)
#     Polynomial(p.data[1] * q.data)
# end
#
# function *(p::Polynomial, q::Polynomial{1})
#     Polynomial(q.data[1] * p.data)
# end

function *(a::Number, q::Polynomial)
    T = promote_type(typeof(a), eltype(q))
    Polynomial{length(q),T}(a*q.data)
end

function +(p::Polynomial, q::Polynomial)
    T = promote_type(eltype(p), eltype(q))
    D = max(length(p), length(q))
    x = zeros(T,D)
    for i in 1 : length(p); x[i] += p.data[i]; end
    for i in 1 : length(q); x[i] += q.data[i]; end
    Polynomial(x...)
end

+(p::Polynomial, a::Number) = p + Polynomial(a)
-(p::Polynomial, a::Number) = p - Polynomial(a)
+(a::Number, p::Polynomial) = p + a
-(a::Number, p::Polynomial) = Polynomial(a) - p

-(p::Polynomial, q::Polynomial) = p + (-one(eltype(q))*q)

function substitute(p::Polynomial, q::Polynomial)
    T = promote_type(eltype(p), eltype(q))
    r = Polynomial(p[degree(p)])
    for d = degree(p)-1 : -1 : 0
        r = q*r + p[d]
    end
    return r
end

mutable struct PieceWisePolynomial
end

function Bernstein(n,u)
    basis = zeros(n+1,1)
    for i = 0:n
        basis[i+1] = binomial(n,i) * u^i * (1.0-u)^(n-i)
    end
    basis
end


text: multiply(a, b) = a * b

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# julia graphs

cd("G:\\My Drive\\18. Github Repo\\DataScience");

# load library
using LightGraphs, Plots, VegaLite, LinearAlgebra, SparseArrays, DataFrames, VegaDatasets, MatrixNetworks

# load data
airports = dataset("airports");
flightsairport = dataset("flights-airport");

# convert to airport
airportsdf = DataFrame(airports);
flightsairportdf = DataFrame(flightsairport);

# show head
first(airportsdf, 10);
first(flightsairportdf, 10);

# only include airport thats in the :origin and :destination
uairports = unique([flightsairportdf[!, :origin]; flightsairportdf[!, :destination]]);
subidx = map(x -> findall(airportsdf[!, :iata] .== uairports[x])[1], 1:length(uairports));
airportsdf_sub = airportsdf[subidx, :];
first(airportsdf_sub, 10)

# build the adjacency matrix
ei_ids = findfirst.(isequal.(flightsairportdf[!,:origin]), [uairports]);
ej_ids = findfirst.(isequal.(flightsairportdf[!,:destination]), [uairports]);
edgeweights = flightsairportdf[!,:count]

A = sparse(ei_ids,ej_ids,1,length(uairports),length(uairports))
A = max.(A,A')

# grapg
G = SimpleGraph(10) #SimpleGraph(nnodes,nedges)
add_edge!(G,7,5) #modifies graph in place.
add_edge!(G,3,5)
add_edge!(G,5,2)

cc = scomponents(A)

degrees = sum(A,dims=2)[:]
p1 = plot(sort(degrees,rev=true),ylabel="log degree",legend=false,yaxis=:log)
p2 = plot(sort(degrees,rev=true),ylabel="degree",legend=false)
plot(p1,p2,size=(600,300))

maxidx = argmax(degrees)
airportsdf_sub[maxidx, :] # ATL

# plot
us10m = dataset("us-10m")
@vlplot(width=500, height=300) +
@vlplot(
    mark={
        :geoshape,
        fill=:lightgray,
        stroke=:white
    },
    data={
        values=us10m,
        format={
            type=:topojson,
            feature=:states
        }
    },
    projection={type=:albersUsa},
) +
@vlplot(
    :circle,
    data=airportsdf_sub,
    projection={type=:albersUsa},
    longitude="longitude:q",
    latitude="latitude:q",
    size={value=10},
    color={value=:steelblue}
)+
@vlplot(
    :rule,
    data=flightsairport,
    transform=[
        {filter={field=:origin,equal=:ATL}},
        {
            lookup=:origin,
            from={
                data=airportsdf_sub,
                key=:iata,
                fields=["latitude", "longitude"]
            },
            as=["origin_latitude", "origin_longitude"]
        },
        {
            lookup=:destination,
            from={
                data=airportsdf_sub,
                key=:iata,
                fields=["latitude", "longitude"]
            },
            as=["dest_latitude", "dest_longitude"]
        }
    ],
    projection={type=:albersUsa},
    longitude="origin_longitude:q",
    latitude="origin_latitude:q",
    longitude2="dest_longitude:q",
    latitude2="dest_latitude:q"
)

# shortest path
ATL_paths = dijkstra(A, maxidx);
ATL_paths[1][maxidx]
maximum(ATL_paths[1])
@show stop1 = argmax(ATL_paths[1])
@show uairports[stop1]

airportsdf_sub[airportsdf_sub[!, :iata] .== "GST", :];

@show stop2 = ATL_paths[2][stop1]
@show uairports[stop2]
airportsdf_sub[airportsdf_sub[!, :iata] .== "JNU", :]

@show stop3 = ATL_paths[2][stop2]
@show uairports[stop3]
airportsdf_sub[airportsdf_sub[!, :iata] .== "SEA", :]

@show stop4 = ATL_paths[2][stop3]
@show uairports[stop4]
airportsdf_sub[airportsdf_sub[!, :iata] .== "ATL", :]

us10m = dataset("us-10m")
airports = dataset("airports")

@vlplot(width=800, height=500) +
@vlplot(
    mark={
        :geoshape,
        fill="#eee",
        stroke=:white
    },
    data={
        values=us10m,
        format={
            type=:topojson,
            feature=:states
        }
    },
    projection={type=:albersUsa},
) +
@vlplot(
    :circle,
    data=airportsdf_sub,
    projection={type=:albersUsa},
    longitude="longitude:q",
    latitude="latitude:q",
    size={value=5},
    color={value=:gray}
) +
@vlplot(
    :line,
    data={
        values=[
            {airport=:ATL,order=1},
            {airport=:SEA,order=2},
            {airport=:JNU,order=3},
            {airport=:GST,order=4}
        ]
    },
    transform=[{
        lookup=:airport,
        from={
            data=airports,
            key=:iata,
            fields=["latitude","longitude"]
        }
    }],
    projection={type=:albersUsa},
    longitude="longitude:q",
    latitude="latitude:q",
    order={field=:order,type=:ordinal}
)

# Minimum Spanning Tree (MST)
@time ti,tj,tv,nverts = mst_prim(A);

df_edges = DataFrame(:ei=>uairports[ti],:ej=>uairports[tj])

@vlplot(width=800, height=500) +
@vlplot(
    mark={
        :geoshape,
        fill="#eee",
        stroke=:white
    },
    data={
        values=us10m,
        format={
            type=:topojson,
            feature=:states
        }
    },
    projection={type=:albersUsa},
) +
@vlplot(
    :circle,
    data=airportsdf_sub,
    projection={type=:albersUsa},
    longitude="longitude:q",
    latitude="latitude:q",
    size={value=20},
    color={value=:gray}
) +
@vlplot(
    :rule,
    data=df_edges, #data=flightsairport,
    transform=[
        {
            lookup=:ei,
            from={
                data=airportsdf_sub,
                key=:iata,
                fields=["latitude", "longitude"]
            },
            as=["originx", "originy"]
        },
        {
            lookup=:ej,
            from={
                data=airportsdf_sub,
                key=:iata,
                fields=["latitude", "longitude"]
            },
            as=["destx", "desty"]
        }
    ],
    projection={type=:albersUsa},
    longitude="originy:q",
    latitude="originx:q",
    longitude2="desty:q",
    latitude2="destx:q"
)

# PageRank
v = MatrixNetworks.pagerank(A,0.85);
insertcols!(airportsdf_sub,7,:pagerank_value => v)

@vlplot(width=500, height=300) +
@vlplot(
    mark={
        :geoshape,
        fill="#eee",
        stroke=:white
    },
    data={
        values=us10m,
        format={
            type=:topojson,
            feature=:states
        }
    },
    projection={type=:albersUsa},
) +
@vlplot(
    :circle,
    data=airportsdf_sub,
    projection={type=:albersUsa},
    longitude="longitude:q",
    latitude="latitude:q",
    size="pagerank_value:q",
    color={value=:steelblue}
)

max_stars_repo_path: 08. Graphs.jl

max_issues_repo_path: 08. Graphs.jl

max_forks_repo_path: 08. Graphs.jl


text:
"""
    ShowHelpMarker(desc)
Helper to display a little (?) mark which shows a tooltip when hovered.
"""
function ShowHelpMarker(desc)
    CImGui.TextDisabled("(?)")
    if CImGui.IsItemHovered()
        CImGui.BeginTooltip()
        CImGui.PushTextWrapPos(CImGui.GetFontSize() * 35.0)
        CImGui.TextUnformatted(desc)
        CImGui.PopTextWrapPos()
        CImGui.EndTooltip()
    end
end

"""
    ShowUserGuide()
Helper to display basic user controls.
"""
function ShowUserGuide()
    CImGui.BulletText("Double-click on title bar to collapse window.")
    CImGui.BulletText("Click and drag on lower right corner to resize window\n(double-click to auto fit window to its contents).")
    CImGui.BulletText("Click and drag on any empty space to move window.")
    CImGui.BulletText("TAB/SHIFT+TAB to cycle through keyboard editable fields.")
    CImGui.BulletText("CTRL+Click on a slider or drag box to input value as text.")
    unsafe_load(CImGui.GetIO().FontAllowUserScaling) && CImGui.BulletText("CTRL+Mouse Wheel to zoom window contents.")
    CImGui.BulletText("Mouse Wheel to scroll.")
    CImGui.BulletText("While editing text:\n")
    CImGui.Indent()
    CImGui.BulletText("Hold SHIFT or use mouse to select text.")
    CImGui.BulletText("CTRL+Left/Right to word jump.")
    CImGui.BulletText("CTRL+A or double-click to select all.")
    CImGui.BulletText("CTRL+X,CTRL+C,CTRL+V to use clipboard.")
    CImGui.BulletText("CTRL+Z,CTRL+Y to undo/redo.")
    CImGui.BulletText("ESCAPE to revert.")
    CImGui.BulletText("You can apply arithmetic operators +,*,/ on numerical values.\nUse +- to subtract.")
    CImGui.Unindent()
end

"""
    ShowAboutWindow(p_open::Ref{Bool})
Access from ImGui Demo -> Help -> About
"""
function ShowAboutWindow(p_open::Ref{Bool})
    CImGui.Begin("About Dear ImGui", p_open, CImGui.ImGuiWindowFlags_AlwaysAutoResize) || (CImGui.End(); return)

    CImGui.Text("Dear ImGui $(CImGui.IMGUI_VERSION)")
    CImGui.Separator()
    CImGui.Text("By Omar Cornut and all dear imgui contributors.")
    CImGui.Text("Dear ImGui is licensed under the MIT License, see LICENSE for more information.")

@cstatic show_config_info=false begin
    @c CImGui.Checkbox("Config/Build Information", &show_config_info)
    if show_config_info
        io = CImGui.GetIO()
        style = CImGui.GetStyle()

        copy_to_clipboard = CImGui.Button("Copy to clipboard")
        CImGui.BeginChildFrame(CImGui.GetID("cfginfos"), ImVec2(0, CImGui.GetTextLineHeightWithSpacing() * 18), CImGui.ImGuiWindowFlags_NoMove)
        copy_to_clipboard && CImGui.LogToClipboard()

        CImGui.Text("Dear ImGui $(CImGui.IMGUI_VERSION)")
        CImGui.Separator()
        CImGui.Text("sizeof(size_t): $(sizeof(Csize_t)), sizeof(ImDrawIdx): $(sizeof(CImGui.ImDrawIdx)), sizeof(ImDrawVert): $(sizeof(CImGui.ImDrawVert))")
        CImGui.Separator()
        CImGui.Text(@sprintf("io.BackendPlatformName: %s", unsafe_load(io.BackendPlatformName) ? unsafe_load(io.BackendPlatformName) : "NULL"))
        CImGui.Text(@sprintf("io.BackendRendererName: %s", unsafe_load(io.BackendRendererName) ? unsafe_load(io.BackendRendererName) : "NULL"))
        CImGui.Text(@sprintf("io.ConfigFlags: 0x%08X", unsafe_load(io.ConfigFlags)))
        unsafe_load(io.ConfigFlags) & CImGui.ImGuiConfigFlags_NavEnableKeyboard    != 0 && CImGui.Text(" NavEnableKeyboard")
        unsafe_load(io.ConfigFlags) & CImGui.ImGuiConfigFlags_NavEnableGamepad     != 0 && CImGui.Text(" NavEnableGamepad")
        unsafe_load(io.ConfigFlags) & CImGui.ImGuiConfigFlags_NavEnableSetMousePos != 0 && CImGui.Text(" NavEnableSetMousePos")
        unsafe_load(io.ConfigFlags) & CImGui.ImGuiConfigFlags_NavNoCaptureKeyboard != 0 && CImGui.Text(" NavNoCaptureKeyboard")
        unsafe_load(io.ConfigFlags) & CImGui.ImGuiConfigFlags_NoMouse              != 0 && CImGui.Text(" NoMouse")
        unsafe_load(io.ConfigFlags) & CImGui.ImGuiConfigFlags_NoMouseCursorChange  != 0 && CImGui.Text(" NoMouseCursorChange")
        unsafe_load(io.MouseDrawCursor)                   && CImGui.Text("io.MouseDrawCursor")
        unsafe_load(io.ConfigMacOSXBehaviors)             && CImGui.Text("io.ConfigMacOSXBehaviors")
        unsafe_load(io.ConfigInputTextCursorBlink)        && CImGui.Text("io.ConfigInputTextCursorBlink")
        unsafe_load(io.ConfigWindowsResizeFromEdges)      && CImGui.Text("io.ConfigWindowsResizeFromEdges")
        unsafe_load(io.ConfigWindowsMoveFromTitleBarOnly) && CImGui.Text("io.ConfigWindowsMoveFromTitleBarOnly")
        CImGui.Text(@sprintf("io.BackendFlags: 0x%08X", unsafe_load(io.BackendFlags)))
        unsafe_load(io.BackendFlags) & ImGuiBackendFlags_HasGamepad      != 0  && CImGui.Text(" HasGamepad")
        unsafe_load(io.BackendFlags) & ImGuiBackendFlags_HasMouseCursors != 0  && CImGui.Text(" HasMouseCursors")
        unsafe_load(io.BackendFlags) & ImGuiBackendFlags_HasSetMousePos  != 0  && CImGui.Text(" HasSetMousePos")
        CImGui.Separator()
        # CImGui.Text("io.Fonts: %d fonts, Flags: 0x%08X, TexSize: %d,%d", io.Fonts->Fonts.Size, io.Fonts->Flags, io.Fonts->TexWidth, io.Fonts->TexHeight);
        CImGui.Text(@sprintf("io.DisplaySize: %.2f,%.2f", unsafe_load(io.DisplaySize).x, unsafe_load(io.DisplaySize).y))
        CImGui.Text(@sprintf("io.DisplayFramebufferScale: %.2f,%.2f", unsafe_load(io.DisplayFramebufferScale).x, unsafe_load(io.DisplayFramebufferScale).y))
        CImGui.Separator()
        CImGui.Text(@sprintf("style.WindowPadding: %.2f,%.2f", style.WindowPadding.x, style.WindowPadding.y))
        CImGui.Text(@sprintf("style.WindowBorderSize: %.2f", style.WindowBorderSize))
        CImGui.Text(@sprintf("style.FramePadding: %.2f,%.2f", style.FramePadding.x, style.FramePadding.y))
        CImGui.Text(@sprintf("style.FrameRounding: %.2f", style.FrameRounding))
        CImGui.Text(@sprintf("style.FrameBorderSize: %.2f", style.FrameBorderSize))
        CImGui.Text(@sprintf("style.ItemSpacing: %.2f,%.2f", style.ItemSpacing.x, style.ItemSpacing.y))
        CImGui.Text(@sprintf("style.ItemInnerSpacing: %.2f,%.2f", style.ItemInnerSpacing.x, style.ItemInnerSpacing.y))

        copy_to_clipboard && CImGui.LogFinish()
        CImGui.EndChildFrame()
    end
end
    CImGui.End()
end


text:
# Bogumił Kamiński, 2019-03-25

function uppersum(M::AbstractMatrix{T}) where {T<:Number}
    n = size(M, 1)
    @assert size(M, 1) == size(M, 1)
    s = zero(T)
    for i in 2:n
        @simd for j in 1:i
            @inbounds s += M[j, i]
        end
    end
    s
end

@time m = randn(10_000, 10_000);
@time m = randn(10_000, 10_000);

m

@time sum(m)
@time sum(m)

@time uppersum(M)
@time uppersum(M)

max_stars_repo_path: 3. Uppersum/uppersum_2.jl

max_issues_repo_path: 3. Uppersum/uppersum_2.jl

max_forks_repo_path: 3. Uppersum/uppersum_2.jl


text:
"""
````
simple!(u,A,b;
                 abstol::Real = zero(real(eltype(b))),
                 reltol::Real = sqrt(eps(real(eltype(b)))),
                 log=false,
                 maxiter=100,
                 P=nothing
                 ) -> solution, [history]
````

Simple iteration scheme ``u_{i+1}= u_i - P^{-1} (A u_i -b)`` with similar API as the methods in IterativeSolvers.jl.

"""
function simple!(u,A,b;
                 abstol::Real = zero(real(eltype(b))),
                 reltol::Real = sqrt(eps(real(eltype(b)))),
                 log=false,
                 maxiter=100,
                 Pl=nothing
                 )
    res=A*u-b # initial residual
    upd=similar(res)
    r0=norm(res) # residual norm
    history=[r0]  # intialize history recording
    for i=1:maxiter
        ldiv!(upd,Pl,res)
	u.-=upd # solve preconditioning system and update solution
	mul!(res,A,u)    # calculate residual
        res.-=b
	r=norm(res)       # residual norm
	push!(history,r)  # record in history
	if (r/r0)<reltol || r<abstol    # check for tolerance
            break
	end
    end
    if log
	return u,Dict( :resnorm =>history )
    else
        return u
    end
end

simple(A,b;kwargs...)=simple!(zeros(eltype(b),length(b)),A,b;kwargs...)


text:
using PGFPlotsX
using CSV
using FinEtools

objects = []

f = CSV.File("$(@__DIR__)/spherical_cap_expl_22-KE-1.0.res")
xdata, ydata = f["t"][1:2:end],  f["v"][1:2:end];
@pgf p = PGFPlotsX.Plot(
{
color = "black",
line_width  = 1.0,
},
Coordinates([v for v in  zip(xdata, ydata)])
)
push!(objects, p)
push!(objects, LegendEntry("1.0"))

f = CSV.File("$(@__DIR__)/spherical_cap_expl_22-KE-0.1.res")
xdata, ydata = f["t"][1:71:end],  f["v"][1:71:end];
@pgf p = PGFPlotsX.Plot(
{
only_marks,
        color = "black",
        line_width  = 1.0,
        mark = "square",
},
Coordinates([v for v in  zip(xdata, ydata)])
)
push!(objects, p)
push!(objects, LegendEntry("0.1"))

f = CSV.File("$(@__DIR__)/spherical_cap_expl_22-KE-10.0.res")
xdata, ydata = f["t"][17:71:end],  f["v"][17:71:end];
@pgf p = PGFPlotsX.Plot(
{
only_marks,
        color = "black",
        line_width  = 1.0,
        mark = "diamond",
},
Coordinates([v for v in  zip(xdata, ydata)])
)
push!(objects, p)
push!(objects, LegendEntry("10.0"))

@pgf ax = Axis(
{
xlabel = "Time [ms]",
ylabel = "Kinetic energy [lbf*in]",
# ymin = -1.1e-5,
# ymax = 1.1e-5,
xmin = 0.0,
xmax = 1.0,
xmode = "linear",
ymode = "linear",
yminorgrids = "true",
grid = "both",
legend_style = {
at = Coordinate(0.5, 1.05),
anchor = "south",
legend_columns = -1
},
},
objects...
)

display(ax)
pgfsave("spherical_cap_expl_22-KE.pdf", ax)


text:
using JLD2

function oscillator()
    A = [0.0 1; -4π 0]
    X0 = BallInf([0.0, 10.0], 0.1)
    X = Universe(2)
    P = @ivp(x' = A*x, x ∈ X, x(0) ∈ X0)
    return P
end

function freedom()
    m1 = 1.0
    m2 = 1.0
    k1 = 1e4
    k2 = 1.0

    M = [m1 0; 0 m2]
    K = [(k1+k2) -k2; -k2 k2]

    sys = SecondOrderLinearContinuousSystem(M, zeros(size(M)), K)

    n = 2
    A = [zeros(n, n) Matrix(1.0I, n, n)
         -inv(M)*K    zeros(n, n)       ]

    X0 = Hyperrectangle([1.0, 10.0, 0.0, 0.0], [0.1, 0.5, 0.5, 0.5])
    P = @ivp(x' = Ax, x(0) ∈ X0)
    P = normalize(P)
    return P
end

function iss()
    ISS_path = joinpath(@__DIR__, "ISS.jld2")
    @load ISS_path A B

    A = Matrix(A)  # some methods cannot deal with sparse matrices
    U = Hyperrectangle(low=[0.0, 0.8, 0.9], high=[0.1, 1., 1.])
    U = ConstantInput(linear_map(B, U))
    X = Universe(270)
    X0 = BallInf(zeros(size(A, 1)), 0.0001)  # -0.0001 <= xi <= 0.0001 for all i
    P = @ivp(x' = A*x + u, x(0) ∈ X0, u ∈ U, x ∈ X)
    return P
end


text:
X = [        0.659027622374092  0.231933368553031
0.659027622374092  0.109039009072877
0.231933368553031  0.659027622374092
0.231933368553031  0.109039009072877
0.109039009072877  0.659027622374092
0.109039009072877  0.231933368553031]
w = [ 0.16666666666666666667
0.16666666666666666667
0.16666666666666666667
0.16666666666666666667
0.16666666666666666667
0.16666666666666666667  ]

function integrate(f)
    s = 0.0
    for i = 1:length(w)
        s += f(X[i,1], X[i,2]) * w[i]
    end
    0.5 * s
end

function calc_basis(ξ, η, x1=1.0, y1=0.0, x2=0.0, y2=1.0, x3=0.0, y3=0.0)
    λ = [ξ;η;1-ξ-η]
    # λ = [1-ξ-η;ξ;η]
    K = 0.5 * abs(det([
        x2-x1 x3-x1
        y2-y1 y3-y1
    ]))
    ∇λ = [
        y2-y3 x3-x2
        y3-y1 x1-x3
        y1-y2 x2-x1
    ]/2K
    ∇λ = [∇λ[:,2] -∇λ[:,1]]
    return λ, ∇λ
end

function φ1(ξ, η)
    λ, ∇λ = calc_basis(ξ, η)
    return λ[1]*∇λ[2,:]
end

function φ2(ξ, η)
    λ, ∇λ = calc_basis(ξ, η)
    return -λ[2]*∇λ[1,:]
end

function φ3(ξ, η)
    λ, ∇λ = calc_basis(ξ, η)
    return λ[2]*∇λ[3,:]
end

function φ4(ξ, η)
    λ, ∇λ = calc_basis(ξ, η)
    return -λ[3]*∇λ[2,:]
end

function φ5(ξ, η)
    λ, ∇λ = calc_basis(ξ, η)
    return λ[3]*∇λ[1,:]
end

function φ6(ξ, η)
    λ, ∇λ = calc_basis(ξ, η)
    return -λ[1]*∇λ[3,:]
end

Basis = [φ1, φ2, φ3, φ4, φ5, φ6]

S = zeros(6,6)
for i = 1:6
    for j = 1:6
        p = (x,y)->Basis[i](x,y)' * Basis[j](x,y)
        S[i,j] = integrate(p)
    end
end
S[ abs.(S) .< 1e-5 ] .= 0.0
@info sum(S), sum(S.^2), sum(S.^3)