text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script creating monthly data files
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, CSV, XLSX, LinearAlgebra, Statistics

# other scripts
dir_home = "x:/VIVES/1-Personal/Florian/git/2022_01_comext/src/"
include(dir_home * "functions.jl")

# location of data input/output (io)
dir_io = "C:/Users/u0148308/data/comext/"
dir_dropbox = "C:/Users/u0148308/Dropbox/BREXIT/"

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# timespan
years = string.(2001:2021)
months = lpad.(1:12, 2, '0')

# subsets string
g(x) = x[end-1:end]

periods = repeat(years, inner=12) .* repeat(months, outer=length(years))
periods = parse.(Int64, periods) # convert to Int64 (actually needed? think we can do all with String)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# create monthly files
# Notes:
#   - import data for 3 months, compute prices and shares
#       + could use 2 months only but then we can only compute MOM. In this setting we obtain also UNIT_PRICE_CHANGE between Jan and Mar,
#         this might be useful later if we realize MOM is too restrictive
#   - subset to middle month data export (i.e. Feb if Jan, Feb, Mar)
#       + use Jan & Feb to calculate MOM change
#       + use Feb & Mar to calculate PRICE_CHANGE (actually when subset and exported we lose this again => use 2 months only and drop PRICE_CHANGE)
#   - lose first and last month (i.e. cannot compute changes since no previous/following data)

for (i, time) in enumerate(periods[1:end-3])

    # specify data to import
    year = string.(time)[1:4]
    month = g.(string.(periods[[i, i+1, i+2]]))
    sub_month = month[2]

    # initialize dataframe
    df = DataFrame(DECLARANT_ISO=String[], PARTNER_ISO=String[], TRADE_TYPE=String[], PRODUCT_NC=String[], FLOW=String[], PERIOD=Int64[],
                    VALUE_IN_EUROS=Union{Missing, Float64}[], QUANTITY_IN_KG=Union{Missing, Float64}[],
                    PRODUCT_SITC=String[], PRODUCT_CPA2008=String[], PRODUCT_CPA2002=String[], PRODUCT_CPA2_1=String[],
                    PRODUCT_BEC=String[], PRODUCT_BEC5=String[], PRODUCT_SECTION=String[])

    # add 3 months to dataframe and do initial cleaning
    append!(df, initial_cleaning(year, month[1]))
    append!(df, initial_cleaning(year, sub_month))
    append!(df, initial_cleaning(year, month[3]))

    #sort!(df)
    sort!(df, [:DECLARANT_ISO, :PARTNER_ISO, :PRODUCT_NC, :FLOW, :TRADE_TYPE, :PERIOD]) # do sorting just in case (unexplainable problem below)

    # -------------------------------------------------------------------------------------------------------------------------------------------------------------
    # Unit prices
    # -------------------------------------------------------------------------------------------------------------------------------------------------------------

    # compute UNIT_PRICE
    transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)

    # compute UNIT_PRICE_CHANGE
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "PRODUCT_NC", "FLOW"]
    gdf = groupby(df, cols_grouping)
    df = transform(gdf, :UNIT_PRICE => pct_change => :UNIT_PRICE_CHANGE)

    # compute MOM_PRICE_CHANGE
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "PRODUCT_NC", "FLOW"]
    gdf = groupby(df, cols_grouping)
    df = transform(gdf, [:PERIOD, :UNIT_PRICE] => mom_change => :MOM)

    # export data
    df_export = subset(df, :PERIOD => ByRow( x-> string(x) == year*sub_month))
    CSV.write(dir_io * "clean/prices/" * "df_prices_" * year * sub_month * ".csv", df_export)

    println(" ✓ Price data for $sub_month/$year was successfully created.")

    # -------------------------------------------------------------------------------------------------------------------------------------------------------------
    # Share of Belgium
    # -------------------------------------------------------------------------------------------------------------------------------------------------------------

    # compute EU_WORLD (EU imports/exports to the entire world)
    cols_grouping = ["PRODUCT_NC", "FLOW", "PERIOD"]
    prod_class_other = ["PRODUCT_SITC", "PRODUCT_CPA2002", "PRODUCT_CPA2008", "PRODUCT_CPA2_1", "PRODUCT_BEC", "PRODUCT_BEC5"]
    cols_grouping = [cols_grouping; prod_class_other]

    gdf = groupby(df, cols_grouping)
    df_sum_EU_WORLD = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)
    df_sum_EU_WORLD.TRADE_TYPE .= "total"

    # compute EU_EXTRA and EU_INTRA (EU imports/exports to extra/intra EU countries)
    cols_grouping = ["TRADE_TYPE", "PRODUCT_NC", "FLOW", "PERIOD"]
    cols_grouping = [cols_grouping; prod_class_other]
    gdf = groupby(df, cols_grouping)
    df_sum_EU_EXTRA_INTRA = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

    # compute BE_WORLD (BE imports/exports to the entire world)
    df_BE = subset(df, :DECLARANT_ISO => ByRow(x -> x == "BE"))
    cols_grouping = ["PRODUCT_NC", "FLOW", "PERIOD"]
    cols_grouping = [cols_grouping; prod_class_other]

    gdf = groupby(df_BE, cols_grouping)
    df_sum_BE_WORLD = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)
    df_sum_BE_WORLD.TRADE_TYPE .= "total"

    # compute BE_EXTRA and BE_INTRA (BE imports/exports to extra/intra EU countries)
    cols_grouping = ["TRADE_TYPE", "PRODUCT_NC", "FLOW", "PERIOD"]
    cols_grouping = [cols_grouping; prod_class_other]
    gdf = groupby(df_BE, cols_grouping)
    df_sum_BE_EXTRA_INTRA = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

    # merge tables
    df_sum_EU = vcat(df_sum_EU_EXTRA_INTRA, df_sum_EU_WORLD)
    df_sum_BE = vcat(df_sum_BE_EXTRA_INTRA, df_sum_BE_WORLD)

    # join, rename and compute shares in percentages (?)
    cols_join = ["TRADE_TYPE", "PRODUCT_NC", "FLOW", "PERIOD"]
    cols_join = [cols_join; prod_class_other]
    df_join = leftjoin(df_sum_EU, df_sum_BE, on=cols_join, makeunique=true)

    rename!(df_join, ["VALUE_IN_EUROS", "QUANTITY_IN_KG", "VALUE_IN_EUROS_1", "QUANTITY_IN_KG_1"] .=> ["VALUE_EU", "QUANITY_EU", "VALUE_BE", "QUANITY_BE"])

    df_join.VALUE_SHARE = df_join.VALUE_BE ./ df_join.VALUE_EU .* 100
    df_join.QUANTITY_SHARE = df_join.QUANITY_BE ./ df_join.QUANITY_EU .* 100

    # possibly need to sort before grouping, seems that grouping sorts incorrectly sometimes?
    sort!(df_join, [:PRODUCT_NC, :FLOW, :TRADE_TYPE, :PERIOD]) # dont understand why I need the sorting here? otherwise some MOM are not taken?

    # compute SHARE_CHANGE
    cols_grouping = ["TRADE_TYPE", "PRODUCT_NC", "FLOW"]
    gdf = groupby(df_join, cols_grouping)
    df_share_BE = transform(gdf, [:VALUE_SHARE, :QUANTITY_SHARE] .=> pct_change .=> [:VALUE_SHARE_CHANGE, :QUANTITY_SHARE_CHANGE])

    # compute MOM_CHANGE
    cols_grouping = ["TRADE_TYPE", "PRODUCT_NC", "FLOW"]
    gdf = groupby(df_share_BE, cols_grouping)
    df_share_BE = transform(gdf, [:PERIOD, :VALUE_SHARE] => mom_change => :MOM_VALUE,
                                 [:PERIOD, :QUANTITY_SHARE] => mom_change => :MOM_QUANTITY)

    # export data
    df_export = subset(df_share_BE, :PERIOD => ByRow( x-> string(x) == year*sub_month))
    CSV.write(dir_io * "clean/shares/" * "df_share_BE_" * year * sub_month * ".csv", df_export)

    println(" ✓ Belgian shares data for $sub_month/$year were successfully created. \n")

end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------


text:
# ---
# title: 1004. Max Consecutive Ones III
# id: problem1004
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Two Pointers, Sliding Window
# link: <https://leetcode.com/problems/max-consecutive-ones-iii/description/>
# hidden: true
# ---
#
# Given an array `A` of 0s and 1s, we may change up to `K` values from 0 to 1.
#
# Return the length of the longest (contiguous) subarray that contains only 1s.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [1,1,1,0,0,0,1,1,1,1,0], K = 2
#     Output: 6
#     Explanation:
#     [1,1,1,0,0, _ **1** ,1,1,1,1, **1**_ ]
#     Bolded numbers were flipped from 0 to 1.  The longest subarray is underlined.
#
# **Example 2:**
#
#
#
#     Input: A = [0,0,1,1,0,0,1,1,1,0,1,1,0,0,0,1,1,1,1], K = 3
#     Output: 10
#     Explanation:
#     [0,0, _1,1, **1** , **1** ,1,1,1, **1** ,1,1_,0,0,0,1,1,1,1]
#     Bolded numbers were flipped from 0 to 1.  The longest subarray is underlined.
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 20000`
#   2. `0 <= K <= A.length`
#   3. `A[i]` is `0` or `1`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
module S3

using AWSCore
include("../Utils.jl")

aws = aws_config()
endpoint = Utils.Endpoint("s3", Utils.rest_xml, "s3.amazonaws.com", "2006-03-01")

# https://github.com/aws/aws-sdk-js/blob/v2.591.0/apis/s3-2006-03-01.min.json#L2421
function list_buckets(args...)
    method = "GET"
    request_uri = "/"

    AWSCore.service_rest_xml(
        aws;
        service      = endpoint.service,
        version      = endpoint.api_version,
        verb         = method,
        resource     = request_uri,
        args         = args
    )
end

# https://github.com/aws/aws-sdk-js/blob/v2.591.0/apis/s3-2006-03-01.min.json#L2648
function list_objects(bucket, args...)
    method = "GET"
    request_uri = "/$bucket"

    AWSCore.service_rest_xml(
        aws;
        service      = endpoint.service,
        version      = endpoint.api_version,
        verb         = method,
        resource     = request_uri,
        args         = args
    )
end

end  # End Module

max_stars_repo_name: mattBrzezinski/AWSCore-Prototype

max_issues_repo_name: mattBrzezinski/AWSCore-Prototype

max_forks_repo_name: mattBrzezinski/AWSCore-Prototype


text:
spectral_resolutions = (31,42,85,170,341)

@testset "FFT of geopotential" begin

    # Test for variable resolution
    for trunc in spectral_resolutions

        P = Parameters(NF=Float64;trunc=trunc)
        G = GeoSpectral(P)
        B = Boundaries(P,G)

        geopot_surf_spectral = B.geopot_surf
        geopot_surf_grid = gridded(geopot_surf_spectral,G)

        @test all(geopot_surf_grid .≈ fourier_inverse(fourier(geopot_surf_grid,G),G))
    end
end

@testset "Legendre transform of geopotential" begin

    # Test for variable resolution
    for trunc in spectral_resolutions

        P = Parameters(NF=Float64;trunc=trunc)
        G = GeoSpectral(P)
        B = Boundaries(P,G)

        geopot_surf_spectral = B.geopot_surf
        geopot_surf_grid = gridded(geopot_surf_spectral,G)

        # using the already spectrally truncated geopotential
        # geopot_surf_grid2 = gridded(spectral(geopot_surf_grid,G),G)

        @test geopot_surf_grid ≈ gridded(spectral(geopot_surf_grid,G),G)

    end
end

@testset "Spectral transform of spectral noise" begin

    # Test for variable resolution
    for trunc in spectral_resolutions

        P = Parameters(NF=Float64;trunc=trunc)
        G = GeoSpectral(P)

        mx = G.spectral.mx
        nx = G.spectral.nx

        A = rand(mx,nx) + im*rand(mx,nx)
        At = spectral(gridded(A,G),G)   # the first transform includes truncation

        # the next should be exact to machine precision
        At2 = spectral(gridded(At,G),G)
        At3 = spectral(gridded(At2,G),G)
        @test all(At .≈ At2)
        @test all(At2 .≈ At3)
    end
end


text:
@testset "Recoverable problem tests" begin
    @test recoverableProblem([1, 2], [5, 3], getKnapsackConstraints([1, 1], 1), 1.0, KnapsackProblemDescriptor(2)) == ([1.0, 0.0], [0.0, 1.0], 4.0)
end


text:
create(c::Client, ::Type{Emoji}, g::AbstractGuild; kwargs...) = create_guild_emoji(c, g.id; kwargs...)

retrieve(c::Client, ::Type{Emoji}, g::AbstractGuild, e::Emoji) = get_guild_emoji(c, g.id, e.id)
retrieve(c::Client, ::Type{Emoji}, g::AbstractGuild) = list_guild_emojis(c, g.id)

update(c::Client, e::Emoji, g::AbstractGuild; kwargs...) = modify_guild_emoji(c, g.id, e.id; kwargs...)

delete(c::Client, e::Emoji, g::AbstractGuild) = delete_guild_emoji(c, g.id, e.id)


text:
# Benchmark suite for BioStructures
# Run with something like:
#   using BioStructures
#   using PkgBenchmark
#   results = benchmarkpkg(".")
#   export_markdown(out_file, results)

using BioStructures
using BenchmarkTools
import BioCore # Imported to avoid clash with BioGenerics distance
using Random

# Use files in BioFmtSpecimens
fmtdir = BioCore.Testing.get_bio_fmt_specimens()
testfilepath(path::AbstractString...) = joinpath(fmtdir, path...)

# All writing is done to one temporary file
temp_filename, io = mktemp()
close(io)

pdbids = ["1AKE", "1EN2", "1SSU"]
formats = Dict("PDB"=> PDB, "mmCIF"=> MMCIF, "MMTF"=> MMTF)
writefunctions = Dict("PDB"=> writepdb, "mmCIF"=> writemmcif, "MMTF"=> writemmtf)

const SUITE = BenchmarkGroup(
        [],
        "read"   => BenchmarkGroup([], [f=> BenchmarkGroup() for f in keys(formats)]...),
        "write"  => BenchmarkGroup([], [f=> BenchmarkGroup() for f in keys(formats)]...),
        "dict"   => BenchmarkGroup(),
        "model"  => BenchmarkGroup(),
        "collect"=> BenchmarkGroup(),
        "spatial"=> BenchmarkGroup(),
)

struc = Dict{String, ProteinStructure}()
for pdbid in pdbids
    struc[pdbid] = read(testfilepath("PDB", "$pdbid.pdb"), PDB)
end

for pdbid in pdbids
    for f in keys(formats)
        SUITE["read"][f][pdbid] = @benchmarkable read(
                $(testfilepath(f, "$pdbid.$(pdbextension[formats[f]])")), $(formats[f]))
        SUITE["write"][f][pdbid] = @benchmarkable $(writefunctions[f])(
                $temp_filename, $(struc[pdbid])) teardown=(rm(temp_filename, force=true))
    end
end

SUITE["dict"]["mmCIF"] = @benchmarkable MMCIFDict($(testfilepath("mmCIF", "1AKE.cif" )))
SUITE["dict"]["MMTF" ] = @benchmarkable MMTFDict( $(testfilepath("MMTF" , "1AKE.mmtf")))

SUITE["model"]["atomaccess"  ] = @benchmarkable $(struc["1EN2"])["A"][20]["CA"]
SUITE["model"]["sortatoms"   ] = @benchmarkable sort(ats) setup=(ats = shuffle(collectatoms(   struc["1EN2"])))
SUITE["model"]["sortresidues"] = @benchmarkable sort(res) setup=(res = shuffle(collectresidues(struc["1EN2"])))
SUITE["model"]["iterate"     ] = @benchmarkable for mod in $(struc["1AKE"]) for ch in mod for res in ch for at in res end end end end
SUITE["model"]["sequence"    ] = @benchmarkable LongAminoAcidSeq($(collectresidues(struc["1EN2"])))
SUITE["model"]["pdbline"     ] = @benchmarkable pdbline($(struc["1EN2"]["A"][20]["CA"]))

SUITE["collect"]["atoms"       ] = @benchmarkable collectatoms(   $(struc["1EN2"]))
SUITE["collect"]["atomssel"    ] = @benchmarkable collectatoms(   $(struc["1EN2"]), calphaselector)
SUITE["collect"]["atomsdis"    ] = @benchmarkable collectatoms(   $(struc["1EN2"]), expand_disordered=true)
SUITE["collect"]["residues"    ] = @benchmarkable collectresidues($(struc["1EN2"]))
SUITE["collect"]["residuessel" ] = @benchmarkable collectresidues($(struc["1EN2"]), standardselector)
SUITE["collect"]["residuesdis" ] = @benchmarkable collectresidues($(struc["1EN2"]), expand_disordered=true)
SUITE["collect"]["chains"      ] = @benchmarkable collectchains(  $(struc["1AKE"]))
SUITE["collect"]["models"      ] = @benchmarkable collectmodels(  $(struc["1SSU"]))

SUITE["spatial"]["coordarray"    ] = @benchmarkable coordarray($(collectatoms(struc["1AKE"])))
SUITE["spatial"]["transformation"] = @benchmarkable Transformation($(collectresidues(struc["1SSU"][5], standardselector)), $(collectresidues(struc["1SSU"][10], standardselector)))
SUITE["spatial"]["rmsd"          ] = @benchmarkable rmsd($(coordarray(struc["1SSU"][5], heavyatomselector)), $(coordarray(struc["1SSU"][10], heavyatomselector)))
SUITE["spatial"]["distance"      ] = @benchmarkable distance($(struc["1AKE"]["A"][50]), $(struc["1AKE"]["A"][60]))
SUITE["spatial"]["bondangle"     ] = @benchmarkable bondangle($(struc["1AKE"]["A"][20]["N"]), $(struc["1AKE"]["A"][20]["CA"]), $(struc["1AKE"]["A"][20]["C"]))
SUITE["spatial"]["dihedralangle" ] = @benchmarkable dihedralangle($(struc["1AKE"]["A"][20]["N"]), $(struc["1AKE"]["A"][20]["CA"]), $(struc["1AKE"]["A"][20]["C"]), $(struc["1AKE"]["A"][21]["N"]))
SUITE["spatial"]["ramachandran"  ] = @benchmarkable ramachandranangles($(collectresidues(struc["1AKE"]["A"], standardselector)))
SUITE["spatial"]["contactmap"    ] = @benchmarkable ContactMap($(collectatoms(struc["1AKE"]["A"], cbetaselector)), 8.0)
SUITE["spatial"]["distancemap"   ] = @benchmarkable DistanceMap($(collectresidues(struc["1AKE"]["A"], standardselector)))
SUITE["spatial"]["contactgraph"  ] = @benchmarkable MetaGraph($(collectatoms(struc["1AKE"]["A"], cbetaselector)), 8.0)


text:
using FactCheck
facts("testing ols.jl") do
    # OLS
    x = [ones(10,1) (1:10)]
    srand(1)
    y = rand(10,1)
    b = ols(y,x,silent=true)[1]
    @fact b[1,1] --> roughly(0.07655377367377458)
    # samin
    include("../examples/samin_example.jl");
    @fact xopt[2] -->roughly(3.25)
    # npreg
    include("../examples/npreg_example.jl");
    @fact yhat[1,1] -->roughly(-0.8108308248103517)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
using GitHub

name = "julia"

if !any(a->startswith(a, "--branch"), ARGS)
    version = v"1.0.3"
    commit_hash = "099e826241fca365a120df9bac9a9fede6e7bae4"
else
    idx = findfirst(a->startswith(a, "--branch"), ARGS)
    branch_name = ARGS[idx+1]
    deleteat!(ARGS, idx:idx+1)
    sanitized_branch = replace(branch_name, "/" => ".")
    version = VersionNumber("1.2.0-$(sanitized_branch)")
    commit_hash = GitHub.reference("JuliaLang/julia", "heads/$(branch_name)").object["sha"]
end

# Collection of sources required to build julia
sources = [
    "https://github.com/JuliaLang/julia.git" =>
    commit_hash,
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
mount -t devpts -o newinstance jrunpts /dev/pts
mount -o bind /dev/pts/ptmx /dev/ptmx
cd julia
cat > Make.user <<EOF
LLVM_ASSERTIONS=1
LIBSSH2_ENABLE_TESTS=0
JULIA_CPU_TARGET=generic;sandybridge,-xsaveopt,clone_all;haswell,-rdrnd,base(1)
USE_BINARYBUILDER_OPENBLAS=1
EOF
make -j${nproc}
make install
cp LICENSE.md ${prefix}
contrib/fixup-libgfortran.sh ${prefix}/lib/julia
contrib/fixup-libstdc++.sh ${prefix}/lib ${prefix}/lib/julia
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "sys", :sys),
    ExecutableProduct(prefix, "julia", :julia)
]

# Dependencies that must be installed before this package can be built
dependencies = []

# Build the tarballs, and possibly a `build.jl` as well.
product_hashes = build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)
tarball, hash = product_hashes["x86_64-linux-gnu"]

version_stanza = """
["$version"]
file = "$tarball"
sha = "$hash"
"""
if isfile(joinpath(@__DIR__, "..", "deps", "downloads", tarball))
  @warn "$tarball already exists in deps/downloads folder. Not copying."
  println("You may manually copy the file from products/ and add the following stanza to Versions.toml:")
  println(version_stanza)
else
  mkpath(joinpath(@__DIR__, "..", "deps", "downloads"))
  cp(joinpath(@__DIR__, "products", tarball), joinpath(@__DIR__, "..", "deps", "downloads", tarball))
  open(joinpath(@__DIR__, "..", "deps", "Versions.toml"); append=true) do f
    println(f, version_stanza)
  end
end


text:
#####
##### `AbstractSynthController`
#####

abstract type AbstractSynthController end

#####
##### `PitchController`
#####

struct PitchController{
    S<:AbstractFloat,
    T<:NTuple{<:Any,Pair{<:AbstractOscillator,<:Any}},
} <: AbstractSynthController
    base_frequency::Ref{S}
    to_set::T

    function PitchController(
        base_frequency::Real,
        to_set::Pair{<:AbstractOscillator,<:Any}...,
    )
        base_frequency = float(base_frequency)
        _setfrequency!(to_set, base_frequency)
        S = typeof(base_frequency)
        T = typeof(to_set)
        return new{S,T}(Ref(base_frequency), to_set)
    end
end

function PitchController(
    base_frequency::Real,
    to_set::Pair{<:AbstractOscillator,<:Real}...,
)
    to_set = map(((osc, c),) -> (osc => Base.Fix1(*, c)), to_set)
    return PitchController(base_frequency, to_set...)
end

function _setfrequency!(
    to_set::NTuple{<:Any,Pair{<:AbstractOscillator,<:Any}},
    frequency::Real
)
    for (osc, f) in to_set
        setfrequency!(osc, f(frequency))
    end
end

getfrequency(p::PitchController) = p.base_frequency[]
function setfrequency!(p::PitchController, frequency::Real)
    _setfrequency!(p.to_set, frequency)
    p.base_frequency[] = float(frequency)
end

@enum KeyState::Bool PRESSED=true RELEASED=false

struct Envelope{F,S<:Real,T<:AbstractFloat}
end


text:
"""
   out = adaptive_pmmh(theta0, prior, state, n; kwargs...)

Generic particle marginal Metropolis-Hastings (PMMH) with adaptive
Metropolis proposal on the parameters.

# Arguments
- `theta0`: Initial parameter vector
- `prior`: Function returning prior log density values for parameters
- `state`: `SMCState` data structure.
- `n`: Number of iterations

# Optional keyword arguments
- `b`: Burn-in length; default `0.1n`
- `thin`: Thinning interval; default `1` (no thinning)
- `save_paths`: Whether to save (samples of) latent trajectories; default `false`
- `show_progress`: Show progress every `show_progress` seconds (default: `false`).

The hidden Markov model and SMC state is defined in `state`; see `?SMCState`.

The output `out` is a `NamedTuple` containing the simulation output: `out.acc`
is the mean acceptance rate, `out.Theta` contains the simulated parameter
values (each column is a parameter vector). If requested (by `save_paths=true`),
`out.X[i][k]` contains the simulated state corresponding to `out.Theta[:,i]`
at time `k`.
"""
function adaptive_pmmh(theta0::ParamT, prior::Function, state, n::Int;
    b::Int=Int(ceil(0.1n)), thin::Int=1, save_paths::Bool=false,
    show_progress::Real=false) where {FT<:AbstractFloat, ParamT<:AbstractVector{FT}}

    # Run the particle filter for HMM w/ param theta0:
    _set_model_param!(state, theta0)
    _run_smc!(state)
    # The initial log-target density:
    p = _log_likelihood(state) + prior(theta0)

    # Initialise random walk Metropolis state
    r = RWMState(theta0)
    # Initialise adaptive Metropolis state
    s = AdaptiveMetropolis(theta0)
    # Allocate space for simulated parameters
    nsim = Int(floor((n-b)/thin))
    Theta = zeros(FT, length(theta0), nsim)
    D = zeros(nsim)
    acc = zero(FT)
    if save_paths
        _pick_particle!(state)
        X = _init_path_storage(state, nsim)
    else
        X = Missing
    end
    if show_progress > 0
        progress = Progress(n, Float64(show_progress))
    else
        progress = NullProgress()
    end
    for k = 1:n
        # PMMH: Draw proposal for parameter, run SMC:
        draw!(r, s); pr_ = prior(r.y)
        if pr_ > -Inf
            _set_model_param!(state, r.y)
            _run_smc!(state)
            p_ = _log_likelihood(state) + pr_
            #alpha = min(one(FT), exp(p_ - p))
            alpha = _accept_prob(p, p_, FT)
            if rand() <= alpha
                accept!(r)
                p = p_
                acc += 1
                if save_paths
                    _pick_particle!(state)
                    _save_reference!(state)
                end
            end
        else
            alpha = 0.0
        end
        adapt!(s, r, alpha, k)

        # Save output:
        @inbounds if k>b && rem(k-b, thin)==0
            i = Int((k-b)/thin)
            Theta[:,i] = r.x
            D[i] = p
            if save_paths
                _copy_reference!(X[i], state)
            end
        end
        next!(progress)
    end
    (Theta=Theta, D=D, acc=acc/n, X=X, am=s, theta0=theta0)
end


text:
function PowellOptimizer(arg0::jdouble, arg1::jdouble)
    return PowellOptimizer((jdouble, jdouble), arg0, arg1)
end

function PowellOptimizer(arg0::jdouble, arg1::jdouble, arg2::ConvergenceChecker)
    return PowellOptimizer((jdouble, jdouble, ConvergenceChecker), arg0, arg1, arg2)
end

function PowellOptimizer(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return PowellOptimizer((jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function PowellOptimizer(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::ConvergenceChecker)
    return PowellOptimizer((jdouble, jdouble, jdouble, jdouble, ConvergenceChecker), arg0, arg1, arg2, arg3, arg4)
end


text:
using DelimitedFiles
include("GB14.jl")
import GLM: fit, predict, LinearModel
using Distributions,SpecialFunctions
using LinearAlgebra, MultivariateStats
using StatsBase, Random, Distributed
using DelimitedFiles, DataStructures
using CSV
using DataFrames

rng = Random.seed!(2019)

gbdt = GB14.GBDT(;
  sampling_rate = 0.75,                # Sampling rate
  learning_rate = 0.1,                # Learning rate
  num_iterations = 100,               # NuGBer of iterations
  tree_options = Dict(
      :nsubfeatures => 0,
      :min_samples_leaf => 1,
      :max_depth => 5))
gbl = GB14.GBLearner(gbdt, :regression)

######## σ² = 1 #########
Pval = []
for i in 1:1000
    if (10 <= i <= 19)
        push!(Pval, i)
    elseif (20 <= i <= 99)
        if (i % 10 == 0)
            push!(Pval, i)
        end
    else
        if (i % 100 == 0)
            push!(Pval, i)
        end

    end
end

Pval
temp_storedim = size(collect(Pval),1)

rmsefitted2 = []
HyperS = Array{Union{Missing, Nothing, Int, Float64}, 2}(nothing, temp_storedim, 1000)

STORE = []

Pval
for pval in Pval
    # reproducibility
    rng = Random.seed!(2019)

    ### simulate the training and testing data
    X, Y = GB14.test_data(500; p = pval, sigma = 1.0)
    X_test, Y_test = GB14.test_data(500; p = pval, sigma = 1.0)
    μ = GB14.fried(X)
    μ_test = GB14.fried(X_test)

    hypers = GB14.Hypers(X, Y, 1, 1, 0.95)
    hypers.σ = sqrt(1)
    hypers.τ = 3 / 10 * var(Y) / gbdt.num_iterations
    hypers.α = 0.95
    hypers.s = ones(size(X, 2)) / size(X, 2)

    # initialize opt
    opts = GB14.Opts()
    opts.num_burn = 50

    # train the model and get the final model as well as rmse2
    rmse2, Store_s = GB14.TBD_predict_and_rmse!(
        X,
        Y,
        X_test,
        Y_test,
        μ_test,
        gbdt,
        hypers,
        opts,
        rng)

    println("The rmse2 b.w. μ_test and pred for p = $(pval)  is $(rmse2).")

    push!(rmsefitted2, rmse2)
    push!(STORE, Store_s)
    # store the hyper.s for the ensembel tree
    j = findall(x -> x == pval, Pval)
    HyperS[j, 1:size(hypers.s, 1)] = hypers.s
    println("After $(opts.num_burn) times of burning, the final sum of the first 5 elements of the prob vec is:\n$(sum(hypers.s[1:5])) ")
end

rmsefitted2
HyperS

HyperS_copy = copy(HyperS)
for i in 1 : size(HyperS, 1)
    for j in 1 : size(HyperS, 2)
        if HyperS_copy[i, j] == nothing
            HyperS_copy[i, j] = 0.0
        end
    end
end
HyperS_copy

rmsefitted2 = DataFrame(rmsefitted2 = rmsefitted2)
CSV.write("_rmse2_sigma1.csv", rmsefitted2)
HyperS = DataFrame(HyperS)
HyperS_copy = DataFrame(HyperS_copy)
CSV.write("_HyperS_copy_sigma1.csv", HyperS_copy)
##### not runned yet

######## σ² = 10 #########
Pval = []
for i in 1:1000
    if (10 <= i <= 19)
        push!(Pval, i)
    elseif (20 <= i <= 99)
        if (i % 10 == 0)
            push!(Pval, i)
        end
    else
        if (i % 100 == 0)
            push!(Pval, i)
        end

    end
end

Pval
temp_storedim = size(collect(Pval),1)
rmsefitted2_sigma10 = []
HyperS_sigma10 = Array{Union{Missing, Nothing, Int, Float64}, 2}(nothing, temp_storedim, 1000)
STORE = []
Pval

##### TODO:Pval = [10, 11,12,13,14,15,16,17,18,19,29,39,49,59,69,79,89,99,209,309,409,509,609,709,809,909,1000]
s
for pval in Pval[1]
    # reproducibility
    rng = Random.seed!(2019)
    ### simulate the training and testing data
    X, Y = GB14.test_data(500; p = pval, sigma = sqrt(10.0))
    X_test, Y_test = GB14.test_data(500; p = pval, sigma = sqrt(10.0))
    μ = GB14.fried(X)
    μ_test = GB14.fried(X_test)

    hypers = GB14.Hypers(X, Y, 1, 1, 0.95)
    hypers.σ = sqrt(1)  ### TODO: I am using
    hypers.τ = 3 / 10 * var(Y) / gbdt.num_iterations
    hypers.α = 0.95
    hypers.s = ones(size(X, 2)) / size(X, 2)

    # initialize opt
    opts = GB14.Opts()
    opts.num_burn = 50

    # train the model and get the final model as well as rmse2
    rmse2, Store_s = GB14.TBD_predict_and_rmse!(
        X,
        Y,
        X_test,
        Y_test,
        μ_test,
        gbdt,
        hypers,
        opts,
        rng)

    println("The rmse2 b.w. μ_test and pred for p = $(pval)  is $(rmse2).")

    push!(rmsefitted2, rmse2)
    push!(STORE, Store_s)
    # store the hyper.s for the ensembel tree
    j = findall(x -> x == pval, Pval)
    HyperS[j, 1:size(hypers.s, 1)] = hypers.s
    println("After $(opts.num_burn) times of burning, the sum of the first 5 elements in the prob vec is: $(sum(hypers.s[1:5])).")
end

rmsefitted2_sigma10
HyperS_sigma10

rmsefitted_sigma10 = DataFrame(rmsefitted_sigma10=rmsefitted_sigma10)
rmsefitted2_sigma10 = DataFrame(rmsefitted2_sigma10 = rmsefitted2_sigma10)
CSV.write("_rmse2_sigma10.csv", rmsefitted2_sigma10)
HyperS_sigma10 = DataFrame(HyperS_sigma10)
HyperS_copy_sigma10 = copy(HyperS_sigma10)

for i in 1 : size(HyperS_sigma10, 1)

    for j in 1 : size(HyperS_sigma10, 2)

        if HyperS_copy_sigma10[i, j] == nothing

            HyperS_copy_sigma10[i, j] = 0.0

        end

    end

end

HyperS_copy_sigma10
CSV.write("_HyperS_sigma10.csv", HyperS_sigma10)

max_stars_repo_name: YinpuLi/SparsityBoosting

max_issues_repo_name: YinpuLi/SparsityBoosting

max_forks_repo_name: YinpuLi/SparsityBoosting


text:
function scan_entry!(result, idx, stmt)
    # NewVarNodes count as defs for the purpose
    # of liveness analysis (i.e. they kill use chains)
    if isa(stmt, NewvarNode)
        result[slot_id(stmt.slot)].any_newvar[] = true
        push!(result[slot_id(stmt.slot)].defs, idx)
        return
    elseif isexpr(stmt, :(=))
        if isa(stmt.args[1], SlotNumber)
            push!(result[slot_id(stmt.args[1])].defs, idx)
        end
        stmt = stmt.args[2]
    end
    if isa(stmt, Union{SlotNumber, TypedSlot})
        push!(result[slot_id(stmt)].uses, idx)
        return
    end
    for op in userefs(stmt)
        val = op[]
        if isa(val, Union{SlotNumber, TypedSlot})
            push!(result[slot_id(val)].uses, idx)
        end
    end
end

struct SlotInfo
    defs::Vector{Int}
    uses::Vector{Int}
    any_newvar::Ref{Bool}
end
SlotInfo() = SlotInfo(Int[], Int[], Ref{Bool}(false))

function lift_defuse(cfg, defuse)
    map(defuse) do slot
        SlotInfo(
            map(x->block_for_inst(cfg, x), slot.defs),
            map(x->block_for_inst(cfg, x), slot.uses),
            Ref{Bool}(slot.any_newvar[])
        )
    end
end

@inline slot_id(s) = isa(s, SlotNumber) ? (s::SlotNumber).id : (s::TypedSlot).id
function scan_slot_def_use(nargs, ci)
    nslots = length(ci.slotnames)
    result = [SlotInfo() for i = 1:nslots]
    # Set defs for arguments
    for var in result[1:(1+nargs)]
        push!(var.defs, 0)
    end
    for (idx, stmt) in enumerate(ci.code)
        scan_entry!(result, idx, stmt)
    end
    result
end

function renumber_ssa(stmt::SSAValue, ssanums, new_ssa=false, used_ssa = nothing)
    id = stmt.id + (new_ssa ? 0 : 1)
    if id > length(ssanums)
        return stmt
    end
    val = ssanums[id]
    if isa(val, SSAValue) && used_ssa !== nothing
        used_ssa[val.id] += 1
    end
    return val
end

function renumber_ssa!(stmt, ssanums, new_ssa=false, used_ssa = nothing)
    isa(stmt, SSAValue) && return renumber_ssa(stmt, ssanums, new_ssa, used_ssa)
    return ssamap(val->renumber_ssa(val, ssanums, new_ssa, used_ssa), stmt)
end

function make_ssa!(ci, idx, slot, typ)
    (idx == 0) && return Argument(slot)
    stmt = ci.code[idx]
    @assert isexpr(stmt, :(=))
    push!(ci.ssavaluetypes, typ)
    ssa = length(ci.ssavaluetypes)-1
    stmt.args[1] = SSAValue(ssa)
    ssa
end

struct UndefToken
end
const undef_token = UndefToken()

function new_to_regular(stmt)
    if isa(stmt, NewSSAValue)
        return SSAValue(stmt.id)
    end
    urs = userefs(stmt)
    urs === () && return stmt
    for op in urs
        val = op[]
        if isa(val, NewSSAValue)
            op[] = SSAValue(val.id)
        end
    end
    urs[]
end

function fixup_slot!(ir, ci, idx, slot, stmt::Union{SlotNumber, TypedSlot}, ssa)
    # We don't really have the information here to get rid of these.
    # We'll do so later
    if ssa === undef_token
        insert_node!(ir, idx, Any, Expr(:throw_undef_if_not, ci.slotnames[slot], false))
        return undef_token
    end
    if !isa(ssa, Argument) && !(ssa === nothing) && ((ci.slotflags[slot] & NI.SLOT_USEDUNDEF) != 0)
        insert_node!(ir, idx, Any, Expr(:undefcheck, ci.slotnames[slot], ssa))
    end
    if isa(stmt, SlotNumber)
        return ssa
    elseif isa(stmt, TypedSlot)
        return NewSSAValue(insert_node!(ir, idx, stmt.typ, PiNode(ssa, stmt.typ)).id)
    end
end

function fixemup!(cond, rename, ir, ci, idx, stmt)
    if isa(stmt, Union{SlotNumber, TypedSlot}) && cond(stmt)
        return fixup_slot!(ir, ci, idx, slot_id(stmt), stmt, rename(stmt))
    end
    if isexpr(stmt, :(=))
        stmt.args[2] = fixemup!(cond, rename, ir, ci, idx, stmt.args[2])
        return stmt
    end
    if isa(stmt, PhiNode)
        for i = 1:length(stmt.edges)
            isassigned(stmt.values, i) || continue
            val = stmt.values[i]
            isa(val, Union{SlotNumber, TypedSlot}) || continue
            cond(val) || continue
            bb_idx = block_for_inst(ir.cfg, stmt.edges[i])
            from_bb_terminator = last(ir.cfg.blocks[bb_idx].stmts)
            stmt.values[i] = fixup_slot!(ir, ci, from_bb_terminator, slot_id(val), val, rename(val))
        end
        return stmt
    end
    urs = userefs(stmt)
    urs === () && return stmt
    for op in urs
        val = op[]
        if isa(val, Union{SlotNumber, TypedSlot}) && cond(val)
            x = fixup_slot!(ir, ci, idx, slot_id(val), val, rename(val))
            # We inserted an undef error node. Delete subsequent statement
            # to avoid confusing the optimizer
            if x === undef_token
                return nothing
            end
            op[] = x
        end
    end
    urs[]
end

function fixup_uses!(ir, ci, uses, slot, ssa)
    for use in uses
        ci.code[use] = fixemup!(stmt->slot_id(stmt)==slot,stmt->ssa,ir,ci,use,ci.code[use])
    end
end

function rename_uses!(ir, ci, idx, stmt, renames)
    fixemup!(stmt->true,stmt->renames[slot_id(stmt)],ir,ci,idx,stmt)
end

function strip_trailing_junk(code)
    # Remove `nothing`s at the end, we don't handle them well
    # (we expect the last instruction to be a terminator)
    for i = length(code):-1:1
        if code[i] !== nothing && !isexpr(code[i], :meta) &&
            !isa(code[i], LineNumberNode) && !isexpr(code[i], :line)
            resize!(code, i)
            break
        end
    end
    # If the last instruction is not a terminator, add one. This can
    # happen for implicit return on dead branches.
    if !isexpr(code[end], :return) && !isexpr(code[end], :gotoifnot) && !isa(code, GotoNode) &&
            !isa(code[end], ReturnNode) && !isa(code[end], GotoIfNot)
        push!(code, ReturnNode{Any}())
    end
    return code
end

struct DelayedTyp
    phi::NewSSAValue
end

function typ_for_val(val, ir, ci)
    isa(val, Expr) && return val.typ
    isa(val, GlobalRef) && return typeof(getfield(val.mod, val.name))
    isa(val, SSAValue) && return ci.ssavaluetypes[val.id+1]
    isa(val, Argument) && return ci.slottypes[val.n]
    isa(val, NewSSAValue) && return DelayedTyp(val)
    return typeof(val)
end

function construct_ssa!(ci, mod, cfg, domtree, defuse, nargs)
    left = Int[]
    defuse_blocks = lift_defuse(cfg, defuse)
    phi_slots = [Vector{Int}() for _ = 1:length(cfg.blocks)]
    phi_nodes = [Vector{Pair{Int,PhiNode}}() for _ = 1:length(cfg.blocks)]
    phi_ssas = SSAValue[]
    code = Any[nothing for _ = 1:length(ci.code)]
    ir = IRCode(code, cfg, ci.slottypes[1:(nargs+1)], mod)
    for (idx, slot) in enumerate(defuse)
        # No uses => no need for phi nodes
        isempty(slot.uses) && continue
        if length(slot.defs) == 1 && slot.any_newvar[]
            if slot.defs[] == 0
                typ = ci.slottypes[idx]
                ssaval = Argument(idx)
                fixup_uses!(ir, ci, slot.uses, idx, ssaval)
            elseif isa(ci.code[slot.defs[]], NewvarNode)
                typ = NI.MaybeUndef(Union{})
                ssaval = nothing
                for use in slot.uses[]
                    insert_node!(ir, use, Union{}, Expr(:throw_undef_if_not, ci.slotnames[idx], false))
                end
                fixup_uses!(ir, ci, slot.uses, idx, nothing)
            else
                val = ci.code[slot.defs[]].args[2]
                typ = typ_for_val(val, ir, ci)
                ssaval = SSAValue(make_ssa!(ci, slot.defs[], idx, typ))
                fixup_uses!(ir, ci, slot.uses, idx, ssaval)
            end
            continue
        end
        # TODO: Perform liveness here to eliminate dead phi nodes
        phiblocks = idf(cfg, defuse_blocks, domtree, idx)
        for block in phiblocks
            push!(phi_slots[block], idx)
            node = PhiNode()
            ssa = insert_node!(ir, first_insert_for_bb(ci.code, cfg, block), Union{}, node)
            push!(phi_nodes[block], ssa.id=>node)
        end
        push!(left, idx)
    end
    # Perform SSA renaming
    initial_incoming_vals = map(1:length(ci.slotnames)) do x
        if 0 in defuse[x].defs
            return Argument(x)
        elseif !defuse[x].any_newvar[]
            return undef_token
        else
            return SSAValue(-1)
        end
    end
    worklist = Any[(1, 0, initial_incoming_vals)]
    visited = Set{Int}()
    type_refine_phi = Set{Int}()
    while !isempty(worklist)
        (item, pred, incoming_vals) = pop!(worklist)
        # Insert phi nodes if necessary
        for (idx, slot) in enumerate(phi_slots[item])
            ssaval, node = phi_nodes[item][idx]
            incoming_val = incoming_vals[slot]
            if incoming_val == SSAValue(-1)
                # Optimistically omit this path.
                # Liveness analysis would probably have prevented us from inserting this phi node
                continue
            end
            push!(node.edges, pred)
            if incoming_val == undef_token
                resize!(node.values, length(node.values)+1)
            else
                push!(node.values, incoming_val)
            end
            # TODO: Remove the next line, it shouldn't be necessary
            push!(type_refine_phi, ssaval)
            if isa(incoming_val, NewSSAValue)
                push!(type_refine_phi, ssaval)
            end
            typ = incoming_val == undef_token ? NI.MaybeUndef(Union{}) : typ_for_val(incoming_val, ir, ci)
            new_node_id = ssaval - length(ir.stmts)
            old_insert, old_typ, _ = ir.new_nodes[new_node_id]
            if isa(typ, DelayedTyp)
                push!(type_refine_phi, ssaval)
            end
            ir.new_nodes[new_node_id] = (old_insert,
                isa(typ, DelayedTyp) ? Union{} : NI.tmerge(old_typ, typ), node)
            incoming_vals[slot] = NewSSAValue(ssaval)
        end
        (item in visited) && continue
        push!(visited, item)
        for idx in cfg.blocks[item].stmts
            stmt = ci.code[idx]
            if isa(stmt, NewvarNode)
                incoming_vals[slot_id(stmt.slot)] = undef_token
                ci.code[idx] = nothing
            else
                ci.code[idx] = rename_uses!(ir, ci, idx, stmt, incoming_vals)
                # Record a store
                if isexpr(stmt, :(=)) && isa(stmt.args[1], SlotNumber)
                    id = slot_id(stmt.args[1])
                    val = stmt.args[2]
                    typ = typ_for_val(val, ir, ci)
                    incoming_vals[id] = SSAValue(make_ssa!(ci, idx, id, typ))
                end
            end
        end
        for succ in cfg.blocks[item].succs
            push!(worklist, (succ, item, copy(incoming_vals)))
        end
    end
    # Delete any instruction in unreachable blocks
    for bb in setdiff(Set{Int}(1:length(cfg.blocks)), visited)
        for idx in cfg.blocks[bb].stmts
            ci.code[idx] = nothing
        end
    end
    # Convert into IRCode form
    ssavalmap = fill(SSAValue(-1), length(ci.ssavaluetypes)+1)
    types = Any[Any for _ = 1:(length(code)+length(ir.new_nodes))]
    # Detect statement positions for assignments and construct array
    for (idx, stmt) in enumerate(ci.code)
        if isexpr(stmt, :(=)) && isa(stmt.args[1], SSAValue)
            ssavalmap[stmt.args[1].id + 1] = SSAValue(idx)
            types[idx] = ci.ssavaluetypes[stmt.args[1].id + 1]
            stmt = stmt.args[2]
            if isa(stmt, PhiNode)
                edges = map(stmt.edges) do edge
                    block_for_inst(cfg, edge)
                end
                code[idx] = PhiNode(convert(Vector{Any}, edges), stmt.values)
            else
                code[idx] = stmt
            end
        # Convert GotoNode/GotoIfNot/PhiNode to BB addressing
        elseif isa(stmt, GotoNode)
            code[idx] = GotoNode(block_for_inst(cfg, stmt.label))
        elseif isa(stmt, GotoIfNot)
            code[idx] = GotoIfNot{Any}(stmt.cond, block_for_inst(cfg, stmt.dest))
        else
            code[idx] = stmt
        end
    end
    # This is a bit awkward, because it basically duplicates what type
    # inference does. Ideally, we'd just use this representation earlier
    # to make sure phi nodes have accurate types
    changed = true
    while changed
        changed = false
        for phi in type_refine_phi
            new_idx = phi - length(ir.stmts)
            old_insert, old_typ, node = ir.new_nodes[new_idx]
            new_typ = Union{}
            for i = 1:length(node.values)
                if !isassigned(node.values, i)
                    if !isa(new_typ, NI.MaybeUndef)
                        new_typ = NI.MaybeUndef(new_typ)
                    end
                    continue
                end
                typ = typ_for_val(node.values[i], ir, ci)
                if isa(typ, DelayedTyp)
                    typ = ir.new_nodes[typ.phi.id - length(ir.stmts)][2]
                end
                new_typ = NI.tmerge(new_typ, typ)
            end
            if !(NI.:⊑)(old_typ, new_typ) || !(NI.:⊑)(new_typ, old_typ)
                #ccall(:jl_, Cvoid, (Any,), (phi, old_typ, new_typ))
                ir.new_nodes[new_idx] = (old_insert, new_typ, node)
                changed = true
            end
        end
    end
    types = map(types) do typ
        isa(typ, DelayedTyp) ? ir.new_nodes[typ.phi.id - length(ir.stmts)][2] : typ
    end
    new_nodes = map(ir.new_nodes) do (pos, typ, node)
        typ = isa(typ, DelayedTyp) ? ir.new_nodes[typ.phi.id - length(ir.stmts)][2] : typ
        (pos, typ, node)
    end
    # Renumber SSA values
    code = map(stmt->new_to_regular(renumber_ssa!(stmt, ssavalmap)), code)
    new_nodes = map(((pt,typ,stmt),)->(pt, typ, new_to_regular(renumber_ssa!(stmt, ssavalmap))), new_nodes)
    IRCode(code, types, ci.slottypes[1:(nargs+1)], cfg, new_nodes, mod)
end


text:
using aoc2019.computer

input = joinpath(@__DIR__, "input")
p = computer.load_program(input)

stdin, stdout = computer.io()
task = @async computer.run(copy(p), stdin, stdout)

while !istaskdone(task)
    sleep(.05)
    while isready(stdout)
        c = take!(stdout)
        print(convert(Char, c))
    end
    print("\$ ")
    cmd = readline()
    for c in cmd put!(stdin, c) end
    put!(stdin, '\n')
end


text:
export Zero

"""
$(TYPEDEF)

Always outputs exactly zero.
"""
struct Zero <: Layer end

Base.hash(a::Zero, h::UInt) = hash(:Zero, h)

function Base.show(io::IO, p::Zero)
    print(io, "Zero()")
end

(p::Zero)(x::Array{<:JuMPReal}) = 0


text:
# mnist_7_or_9_sgfs
#
# This script compares the Stein discrepancy of sample sequences drawn
# from SGFS as done in Experiment 5.1 of (http://arxiv.org/abs/1206.6380).
# We got the Matlab code from the authors and have used their
# code to re-generate the same samples as used in the paper.
# By plotting the Stein discrepancy results against
# a measure of computation required per sample (e.g., the number of
# datapoint likelihood evaluations needed), we can quantify
# which of the samplers is indeed best for representing the target.
#
# The target distribution is a Bayesian logistic regression posterior
# under independent Gaussian priors based on a the mnist dataset
# with either 7 or 9 as the labels.

# parallelize now!
addprocs(CPU_CORES - 1)

using DataFrames

using SteinDistributions:
    SteinLogisticRegressionGaussianPrior,
    gradlogdensity
using SteinDiscrepancy:
    ksd,
    SteinGaussianKernel,
    SteinGaussianPowerKernel,
    SteinMaternRadialKernel,
    SteinInverseMultiquadricKernel

BASEPATH = "src/experiments/mnist"

include("experiment_utils.jl")

# Set random number generator seed
@parseintcli seed "s" "seed" 7
# set the sampler
@parsestringcli sampler "r" "sampler" "SGFS-f"
# set the kernel
@parsestringcli kernelname "k" "kernel" "inversemultiquadric"
kernel = nothing
# Setup the data
data_set = "mnist"
Xtrain = convert(
    Array,
    readtable(joinpath(BASEPATH, "mnist79Xtrain50proj.csv"), header=false)
)
ytrain = convert(
    Array,
    readtable(joinpath(BASEPATH, "mnist79ytrain50proj.csv"), header=false)
)
# get the training sizes
ntrain = size(Xtrain, 1)
# give X an intercept term
X = hcat(Xtrain, ones(ntrain))
# make y be -1 or 1
y = vec(2.0 * ytrain - 1.0)
# Parameter dimension
d = size(X, 2)
# name of prior + posterior
distname = "logisticgaussianprior"
# sigma2 term for prior
sigma2 = 1.0
# setup distribution
target = SteinLogisticRegressionGaussianPrior(X, y, zeros(d), sigma2 .* ones(d))
# define gradlogp
function gradlogp(x::Array{Float64,1})
    gradlogdensity(target, x)
end
# number of burnin samples
b = 50000

if sampler in ["SGFS-d", "SGFS-f"]
    sample_filename = joinpath(
        BASEPATH,
        @sprintf("mnist79_%s_n=100000.csv", sampler)
    )
    betas = convert(Array, readtable(sample_filename, header=false))
    betas = betas[(b+1):end,:]
    n = size(betas,1)

    # checkpoints for the kernel
    checkpoints = vcat(100:100:min(1000,n),
                       2000:1000:min(n,10000),
                       20000:10000:min(n,100000))
    if checkpoints[end] != n
        push!(checkpoints, n)
    end
    # get the kernel
    if kernelname == "gaussian"
        kernel = SteinGaussianKernel()
    elseif kernelname == "gaussianpower"
        kernel = SteinGaussianPowerKernel(2.0)
    elseif kernelname == "maternradial"
        kernel = SteinMaternRadialKernel()
    elseif kernelname == "inversemultiquadric"
        kernel = SteinInverseMultiquadricKernel(0.5)
    end
    # Compute Stein discrepancy
    res = @setseed seed ksd(
        points=betas,
        gradlogdensity=gradlogp,
        kernel=kernel,
        checkpoints=checkpoints)

    betahat = vec(mean(betas, 1))
    sigmahat = cov(betas)

    instance_data = Dict{Any, Any}(
        "distname" => distname,
        "n" => checkpoints,
        "d" => d,
        "betahat" => betahat,
        "sigmahat" => sigmahat,
        "kernel" => kernelname,
        "betas" => betas,
        "discrepancy" => sqrt(res.discrepancy2),
        "solvetime" => res.solvetime,
        "sampler" => sampler,
        "ncores" => nprocs(),
        "burnin" => b,
    )

    save_json(
        instance_data;
        dir="mnist_7_or_9_sgfs",
        distname=distname,
        dataset=data_set,
        kernel=kernelname,
        sampler=sampler,
        n=n,
        d=d,
        ncores=nprocs()
    )
end


text:
module CUDAnative

using CUDAapi
using CUDAdrv

using LLVM
using LLVM.Interop

using GPUCompiler

using Adapt

## deferred initialization

# CUDA packages require complex initialization (discover CUDA, download artifacts, etc)
# that can't happen at module load time, so defer that to run time upon actual use.

const configured = Ref{Union{Nothing,Bool}}(nothing)

"""
    functional(show_reason=false)

Check if the package has been configured successfully and is ready to use.

This call is intended for packages that support conditionally using an available GPU. If you
fail to check whether CUDA is functional, actual use of functionality might warn and error.
"""
function functional(show_reason::Bool=false)
    if configured[] === nothing
        _functional(show_reason)
    end
    configured[]::Bool
end

const configure_lock = ReentrantLock()
@noinline function _functional(show_reason::Bool=false)
    lock(configure_lock) do
        if configured[] === nothing
            if __configure__(show_reason)
                configured[] = true
                try
                    __runtime_init__()
                catch
                    configured[] = false
                    rethrow()
                end
            else
                configured[] = false
            end
        end
    end
end

# macro to guard code that only can run after the package has successfully initialized
macro after_init(ex)
    quote
        @assert functional(true) "CUDAnative.jl did not successfully initialize, and is not usable."
        $(esc(ex))
    end
end

## source code includes

include("init.jl")
include("compatibility.jl")
include("bindeps.jl")

include("cupti/CUPTI.jl")
include("nvtx/NVTX.jl")

# needs to be loaded _before_ the compiler infrastructure, because of generated functions
include("device/pointer.jl")
include("device/array.jl")
include("device/cuda.jl")
include("device/llvm.jl")
include("device/runtime.jl")

include("compiler.jl")
include("execution.jl")
include("exceptions.jl")
include("reflection.jl")
include("array.jl")

include("deprecated.jl")

export CUPTI, NVTX

## initialization

# device compatibility
const __target_support = Ref{Vector{VersionNumber}}()
const __ptx_support = Ref{Vector{VersionNumber}}()

target_support() = @after_init(__target_support[])
ptx_support() = @after_init(__ptx_support[])

function __init__()
    if Base.libllvm_version != LLVM.version()
        error("LLVM $(LLVM.version()) incompatible with Julia's LLVM $(Base.libllvm_version)")
    end

    # enable generation of FMA instructions to mimic behavior of nvcc
    LLVM.clopts("-nvptx-fma-level=1")

    resize!(thread_contexts, Threads.nthreads())
    fill!(thread_contexts, nothing)

    resize!(thread_tasks, Threads.nthreads())
    fill!(thread_tasks, nothing)

    CUDAdrv.initializer(prepare_cuda_call)
end

function __configure__(show_reason::Bool)
    # if any dependent GPU package failed, expect it to have logged an error and bail out
    if !CUDAdrv.functional(show_reason)
        show_reason && @warn "CUDAnative.jl did not initialize because CUDAdrv.jl failed to"
        return false
    end

    return __configure_dependencies__(show_reason)
end

function __runtime_init__()
    if release() < v"9"
        @warn "CUDAnative.jl only supports CUDA 9.0 or higher (your toolkit provides CUDA $(release()))"
    elseif release() > CUDAdrv.release()
        @warn """You are using CUDA toolkit $(release()) with a driver that only supports up to $(CUDAdrv.release()).
                 It is recommended to upgrade your driver, or switch to automatic installation of CUDA."""
    end

    # device compatibility

    llvm_support = llvm_compat()
    cuda_support = cuda_compat()

    __target_support[] = sort(collect(llvm_support.cap ∩ cuda_support.cap))
    isempty(__target_support[]) && error("Your toolchain does not support any device capability")

    __ptx_support[] = sort(collect(llvm_support.ptx ∩ cuda_support.ptx))
    isempty(__ptx_support[]) && error("Your toolchain does not support any PTX ISA")

    @debug("Toolchain with LLVM $(LLVM.version()), CUDA driver $(CUDAdrv.version()) and toolkit $(CUDAnative.version()) supports devices $(verlist(__target_support[])); PTX $(verlist(__ptx_support[]))")
end

end


text:
module PlutoUtils

using Comonicon

# run(host, port::Integer; launchbrowser::Bool=false, session=ServerSession())

"""
Pluto CLI - Lightweight reactive notebooks for Julia

# Arguments

- `host`: default is 127.0.0.1 (localhost)

# Options

- `-p,--port <int>`: port you want to specify, default is 1234

# Flags

- `-l,--launchbrowser`: add this flag to launch browser.
"""
@main function pluto(;host="127.0.0.1", port::Int=1234, launchbrowser::Bool=false)
    # workaround the package environment problem
    # might be resolved by Pluto#142
    julia = joinpath(Sys.BINDIR::String, Base.julia_exename())
    script = """
    try
        import Pluto
    catch e
        if e isa ArgumentError
            print("Pluto not found, install Pluto? [Y/n]")
            x = read(stdin, Char)
            if x in ['Y', 'y', '\n']
                using Pkg
                Pkg.add("Pluto")
                import Pluto
            else
                exit(0)
            end
        end
    end

    Pluto.run("$host", $port; launchbrowser=$launchbrowser)
    """

    PLUTO_ENV = filter(x->x.first!="JULIA_PROJECT", ENV)
    run(setenv(`$julia -e "$script"`, PLUTO_ENV))
end

end


text:
using Particula, Unitful, Test

@testset "Particula.jl" begin
    # Write your tests here.
    @test Particula.R |> ustrip >= 0

end

include("test_constants.jl")
include("test_inputting.jl")
include("test_dviscosity.jl")
include("test_meanfreepath.jl")


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

using Markdown
using InteractiveUtils

# ╔═╡ acd7da78-1d69-11eb-2517-e1bafdc14f4a
using Pkg, DrWatson

# ╔═╡ acd82a8c-1d69-11eb-10c9-0902b833d79e
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Plots, LaTeXStrings, Measures
end

# ╔═╡ 2e6c02cc-1d69-11eb-283a-1bcf9e0044bf
md"## Listing 1.10"

# ╔═╡ acd8c1cc-1d69-11eb-3164-5b9e66fc3664
begin
	f(x,y) = x^2 + y^2
	f0(x) = f(x,0)
	f2(x) = f(x,2)
end

# ╔═╡ ace3c478-1d69-11eb-0a4d-f191cf032e78
begin
	x, y = -5:0.1:5 , -5:0.1:5
	plot(x, [f0.(x), f2.(x)],
		c=[:blue :red], xlims=(-5, 5), legend=:top,
		ylims=(-5, 25), ylabel=L"f(x,⋅)", label=["f(x, 0)" "f(x, 2)"])
	p1 = annotate!(-1, -0.2, text("Minimum = (0,0)", :left, :top, 10))
end

# ╔═╡ ace462fa-1d69-11eb-2b31-c1c9efe2376f
begin
	z = [ f(x,y) for y in y, x in x ]
	p2 = surface(x, y, z, c=cgrad([:blue, :red]), legend=:none, ylabel="y", zlabel="f(x,y)")
end

# ╔═╡ acf332f0-1d69-11eb-09df-43d188f0b58d
begin
	M = z[1:10,1:10]
	p3 = heatmap(M, c=cgrad([:blue, :red]), yflip=true, ylabel="y",
		xticks=([1:10;], x), yticks=([1:10;], y))
end

# ╔═╡ ad012874-1d69-11eb-24a0-e936c7c266d5
md"## End of listing 1.10"

# ╔═╡ Cell order:
# ╟─2e6c02cc-1d69-11eb-283a-1bcf9e0044bf
# ╠═acd7da78-1d69-11eb-2517-e1bafdc14f4a
# ╠═acd82a8c-1d69-11eb-10c9-0902b833d79e
# ╠═acd8c1cc-1d69-11eb-3164-5b9e66fc3664
# ╠═ace3c478-1d69-11eb-0a4d-f191cf032e78
# ╠═ace462fa-1d69-11eb-2b31-c1c9efe2376f
# ╠═acf332f0-1d69-11eb-09df-43d188f0b58d
# ╟─ad012874-1d69-11eb-24a0-e936c7c266d5


text:
using PackedASCIITypes
using Test

const TYPEINFO = [(PackedASCII1, 1),
                  (PackedASCII2, 2),
                  (PackedASCII4, 4),
                  (PackedASCII8, 8),
                  (PackedASCII17, 17)]

@testset "non-standard literals and creation" begin
    @test pasc1"a"::PackedASCII1 ≡ pasc"a"::PackedASCII1 ≡ PackedASCII1("a") ≡
        PackedASCII("a")
    @test pasc2"ab"::PackedASCII2 ≡ pasc"ab"::PackedASCII2 ≡
        PackedASCII2("ab") ≡ PackedASCII("ab")
    @test pasc4"abcd"::PackedASCII4 ≡ pasc"abcd"::PackedASCII4 ≡
        PackedASCII4("abcd") ≡ PackedASCII("abcd")
    @test pasc8"abcdefgh"::PackedASCII8 ≡ pasc"abcdefgh"::PackedASCII8 ≡
        PackedASCII8("abcdefgh") ≡ PackedASCII("abcdefgh")
    @test pasc17"fourtytwo"::PackedASCII17 ≡ pasc"fourtytwo"::PackedASCII17 ≡
        PackedASCII17("fourtytwo") ≡ PackedASCII("fourtytwo")
end

@testset "small tests" begin
    for (T, l) in TYPEINFO
        @test PackedASCIITypes.maxlength(T) == l
        str = reduce(*, rand('a':'z', l))
        p = T(str)
        @test p isa T
        @test String(p) == str
        @test_throws InexactError T(str * 'x')
        a = T("a")
        b = T("b")
        @test a ≡ a
        @test a == a
        @test b ≠ a
        if l ≥ 2
            aa = T("aa")
            @test a < aa < b
        else
            @test a < b
        end
        @test String(typemin(T)) == ""
        @test repr(a) == "pasc$(l)\"a\""
    end
end

@testset "mixing types, equality, hashing" begin
    a = pasc2"a"
    aa = pasc8"aa"
    b = pasc4"b"
    @test aa == pasc4"aa"
    @test a < b
    @test a < aa
    for T in first.(TYPEINFO)
        a2 = T(a)
        @test a == a2
        @test a2 < aa < T(b)
        @test hash(a) == hash(a2)
        h = rand(UInt)
        @test hash(a, h) == hash(a2, h)
    end
end

@testset "constructors" begin
    @test_throws InexactError PackedASCII("1"^18)
    @test_throws InexactError PackedASCII("pillangó")
end

@testset "I/O" begin
    io = IOBuffer()
    a = pasc1"a"
    for T in first.(TYPEINFO)
        write(io, T(a))
    end
    b = take!(io)
    @test length(b) == 2^5-1
    i = 0
    for T in first.(TYPEINFO)
        j = sizeof(T)
        a2 = reinterpret(T, b[i.+(1:j)])[1]
        @test T(a) ≡ a2::T
        i += j
    end
end


text:
"""
    posterior(Cidx::Int, Pidxs::Vector{Int}, mcr::ModelClusteringResult, X::AbstractMatrix) -> Float64

Calculate posterior of the clusrer `C` given the PLM cluster `P`, p(C|P).
Clusters are given by their respective indices in the model class presented as model-base clustering 'mcr'.
"""
function posterior(Cidx::Int, Pidxs::Vector{Int},
                   mcr::ModelClusteringResult,
                   X::AbstractMatrix,
                   marginal::MixtureModel)
    Ppts = findall(i-> i ∈ Pidxs, assignments(mcr))
    # Ppts = mapreduce(i->clusters(mcr)[i].idx, vcat, Pidxs)
    Cdist = components(marginal)[Cidx]
    pc = probs(marginal)[Cidx]
    return sum(i->
        let x = X[:,i]
            pcx = pdf(Cdist, x)
            px = pdf(marginal, x)
            pcx*pc/px
        end, Ppts)/length(Ppts)
end

function posterior(Cidx::Int, marginal::MixtureModel, X::AbstractMatrix)
    pn = size(X,2)
    Cdist = components(marginal)[Cidx]
    pc = probs(marginal)[Cidx]
    pc*sum(pdf(Cdist, X) ./ pdf(marginal, X))/pn
end

posterior(Cidx::Int, Pidxs::Vector{Int}, cnt::Vector{Int}) =
    Cidx ∉ Pidxs ? 0.0 : cnt[Cidx]/sum(i->cnt[i], Pidxs)

function posterior(Cidx::Int, priorC::AbstractVector{T}, logCs::AbstractMatrix{T},
                   mgnl::AbstractVector{T}, idxs::Vector{Bool}) where {T <: AbstractFloat}
    n = sum(idxs)
    logCP = view(logCs, idxs, :)
    priorC[Cidx]*sum(exp.(view(logCP,:,Cidx))./view(mgnl,idxs))/n
end

function ibdiff(P′::Vector{Vector{Int}},
                mcr::ModelClusteringResult,
                X::AbstractMatrix; kwargs...)
    logps = hcat((logpdf(p, X) for p in models(mcr))...)
    return ibdiff(P′, logps, assignments(mcr); kwargs...)
end

function ibdiff(P′::Vector{Vector{Int}}, logCs::AbstractMatrix{T},
                assign::Vector{Int}; hard=true,
                marginal::Union{AbstractVector{T}, Nothing}=nothing,
                kwargs...)  where {T <: AbstractFloat}
    # precompute distributions
    n, k = size(logCs)
    cs = zeros(Int, k)
    for i in assign
        cs[i] += 1
    end
    pC = map(c->c/n, cs)
    pP′ = sum(sum(pC[P]) for P in P′)
    if marginal === nothing && !hard
        marginal = sum(exp.(logCs' .+ log.(pC)), dims=1)'
    end
    # calculate IB difference
    P′idxs = vcat(P′...)
    ibd = 0.0
    for C in 1:k
        pCP′ = if hard
            posterior(C, P′idxs, cs)
        else
            P′pts = map(i-> i ∈ P′idxs, assign)
            posterior(C, pC, logCs, marginal, P′pts)
        end
        Dₖₗ2 = pCP′ == 0.0 ? 0.0 : pCP′*log(pCP′/pP′)
        for P in P′
            pP  = sum(pC[P])
            pCP = if hard
                posterior(C, P, cs)
            else
                Ppts = map(i-> i ∈ P, assign)
                posterior(C, pC, logCs, marginal, Ppts)
            end
            Dₖₗ1 = pCP == 0.0 ? 0.0 : pCP*log(pCP/pP)
            ibd += pC[C]*(Dₖₗ1 - Dₖₗ2)
        end
        # @debug "IB Difference" cluster=C IB=ibd
    end
    return ibd
end

max_stars_repo_name: wildart/PLMC

max_issues_repo_name: wildart/PLMC

max_forks_repo_name: wildart/PLMC


text:
#Test against Convex.jl

using Convex

pkgdir(pkg::String) = abspath(joinpath(dirname(Base.find_package(pkg)), ".."))

single_solver_test_file = joinpath(pkgdir("Convex"),"test/runtests_single_solver.jl")
set_default_solver(DR(eps=1e-8, verbose=0, debug=0))

include(single_solver_test_file)

# TODO TODO
# #TODO Throws error, integrate with Base.Test
# FactCheck.exitstatus()


text:
export Dirichlet

import Distributions: Dirichlet
import SpecialFunctions: digamma

vague(::Type{ <: Dirichlet }, dims::Int) = Dirichlet(ones(dims))

prod_analytical_rule(::Type{ <: Dirichlet }, ::Type{ <: Dirichlet }) = ProdAnalyticalRuleAvailable()

function prod(::ProdAnalytical, left::Dirichlet, right::Dirichlet)
    mvec = probvec(left) .+ probvec(right)
    mvec = mvec .- one(eltype(mvec))
    return Dirichlet(mvec)
end

probvec(dist::Dirichlet) = params(dist)[1]

# probvec is not normalised
mean(::typeof(log), dist::Dirichlet) = digamma.(probvec(dist)) .- digamma(sum(probvec(dist)))

# Variate forms promotion

promote_variate_type(::Type{ Multivariate },  ::Type{ <: Dirichlet })          = Dirichlet
promote_variate_type(::Type{ Matrixvariate }, ::Type{ <: Dirichlet })          = MatrixDirichlet

promote_variate_type(::Type{ Multivariate },  ::Type{ <: MatrixDirichlet }) = Dirichlet
promote_variate_type(::Type{ Matrixvariate }, ::Type{ <: MatrixDirichlet }) = MatrixDirichlet


text:
using Statistics
using NCDatasets

using ArgParse
using JSON

include("nanop.jl")
include("PCA.jl")

using .PCA

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--data-file"
            help = "Ocean SSTA data file."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

Dataset(parsed["domain-file"], "r") do ds
    global lon = replace(ds["xc"][:, 1], missing=>NaN)
    global lat = replace(ds["yc"][1, :], missing=>NaN)
    global mask = replace(ds["mask"][:], missing=>NaN)

    global Nx = length(lon)
    global Ny = length(lat)

    global valid_grids = []

    for i = 1:Nx, j = 1:Ny
        if mask[i, j] == 0.0 && ( 115.0 <= lon[i] <= 255.0 ) && ( 20.0 <= lat[j] <= 65.0 )
            push!(valid_grids, (i, j))
        end
    end
end

Dataset(parsed["data-file"], "r") do ds

    global SSTA  = replace(ds["T_ML_MA"][:, :, 1, :], missing=>NaN)
    global months  = ds.dim["time"]

    global EOF_input = zeros(Float64, length(valid_grids), months)
    for i = 1:length(valid_grids)
        EOF_input[i, :] = SSTA[valid_grids[i][1], valid_grids[i][2], :]
    end
end

modes = 2

println("Solving for PCA...")
eigen_vectors = PCA.findPCAs(EOF_input, num=modes)
println("done.")

PCAs = zeros(Float64, Nx, Ny, modes)
PCAs .= NaN

for i = 1:length(valid_grids)
    PCAs[valid_grids[i][1], valid_grids[i][2], :] = eigen_vectors[i, :]
end

PDO = zeros(Float64, months)

PDO_mode = view(PCAs, :, :, 1)
for i = 1:length(PDO)
    PDO[i] = nansum(PDO_mode .* SSTA[:, :, i])
end

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "Nx", Nx)
    defDim(ds, "Ny", Ny)
    defDim(ds, "time", Inf)
    defDim(ds, "modes", modes)

    for (varname, vardata, vardim, attrib) in [
        ("PCAs",  PCAs, ("Nx", "Ny", "modes",), Dict()),
        ("PDO",   PDO, ("time",), Dict()),
    ]

        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
# %%

for i in 1:10
    println(i)
end

# %%

a = split("A,B,C,D",",")
@show typeof(a)
@show a

@show [1,2,3]
@show typeof([1,2,3])
@show [1 2 3]
@show typeof([1 2 3])

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# ---
# title: 303. Range Sum Query - Immutable
# id: problem303
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/range-sum-query-immutable/description/>
# hidden: true
# ---
#
# Given an integer array `nums`, find the sum of the elements between indices
# `i` and `j` `(i ≤ j)`, inclusive.
#
# Implement the `NumArray` class:
#
#   * `NumArray(int[] nums)` Initializes the object with the integer array `nums`.
#   * `int sumRange(int i, int j)` Return the sum of the elements of the `nums` array in the range `[i, j]` inclusive (i.e., `sum(nums[i], nums[i + 1], ... , nums[j])`)
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["NumArray", "sumRange", "sumRange", "sumRange"]
#     [[[-2, 0, 3, -5, 2, -1]], [0, 2], [2, 5], [0, 5]]
#     **Output**
#     [null, 1, -1, -3]
#
#     **Explanation**
#     NumArray numArray = new NumArray([-2, 0, 3, -5, 2, -1]);
#     numArray.sumRange(0, 2); // return 1 ((-2) + 0 + 3)
#     numArray.sumRange(2, 5); // return -1 (3 + (-5) + 2 + (-1))
#     numArray.sumRange(0, 5); // return -3 ((-2) + 0 + 3 + (-5) + 2 + (-1))
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 104`
#   * `-105 <= nums[i] <= 105`
#   * `0 <= i <= j < nums.length`
#   * At most `104` calls will be made to `sumRange`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
info("Testing ", Bound)
for T in (FloatingPointTypes..., IntegerTypes...)
    for isopen in (true, false)
        B = Bound(one(T), isopen)
        @test B.value == one(T)
        @test B.isopen == isopen
        @test B == Bound(one(T), isopen ? :open : :closed)

        if T in FloatingPointTypes
            @test eltype(convert(Bound{Float32}, B)) == Float32
            @test eltype(convert(Bound{Float64}, B)) == Float64
        else
            @test eltype(convert(Bound{Int32}, B)) == Int32
            @test eltype(convert(Bound{Int64}, B)) == Int64
        end

        @test_throws ErrorException Bound(one(T), :test)
    end
end

info("Testing ", Interval)
for T in (FloatingPointTypes..., IntegerTypes...)
    for lisopen in (true, false), uisopen in (true, false)
        Bl = Bound(zero(T), lisopen)
        Br = Bound(convert(T,2),  uisopen)

        I = Interval(Bl, Br)

        @test eltype(I) == T

        @test get(I.left) == Bl
        @test get(I.right) == Br

        @test MOD.checkbounds(I, -one(T)) == false
        @test MOD.checkbounds(I, convert(T,3)) == false

        @test MOD.checkbounds(I, one(T))

        if lisopen
            @test MOD.checkbounds(I, zero(T)) == false
        else
            @test MOD.checkbounds(I, zero(T)) == true
        end

        if uisopen
            @test MOD.checkbounds(I, convert(T,2)) == false
        else
            @test MOD.checkbounds(I, convert(T,2)) == true
        end

        for lnull in (true, false), unull in (true, false)
            I = Interval(lnull ? Nullable{Bound{T}}() : Nullable(Bl),
                         unull ? Nullable{Bound{T}}() : Nullable(Br))
            @test lnull ? isnull(I.left) : get(I.left) == Bl
            @test unull ? isnull(I.right) : get(I.right) == Br

            if T in FloatingPointTypes
                I2 = convert(Interval{Float32}, I)
                @test lnull ? isnull(I2.left) : get(I2.left) == convert(Bound{Float32}, Bl)
                @test unull ? isnull(I2.right) : get(I2.right) == convert(Bound{Float32}, Br)
                I2 = convert(Interval{Float64}, I)
                @test lnull ? isnull(I2.left) : get(I2.left) == convert(Bound{Float64}, Bl)
                @test unull ? isnull(I2.right) : get(I2.right) == convert(Bound{Float64}, Br)
            end

            # Test that output does not create error
            show(DevNull, I)
        end

        Br = Bound(zero(T), uisopen)

        if lisopen || uisopen
            @test_throws ErrorException Interval(Bl, Br)
        else
            I = Interval(Bl, Br)
            @test get(I.left) == Bl
            @test get(I.right) == Br
        end
    end
    for isopen in (true, false)
        B = Bound(convert(T,2), isopen)

        I = leftbounded(B)
        @test isnull(I.right)
        @test get(I.left) == B
        @test MOD.checkbounds(I, one(T)) == false
        @test MOD.checkbounds(I, convert(T,3)) == true
        @test MOD.checkbounds(I, convert(T,2)) == (isopen ? false : true)

        I = leftbounded(convert(T,2), isopen ? :open : :closed)
        @test isnull(I.right)
        @test get(I.left) == B

        @test_throws ErrorException leftbounded(one(T), :test)

        I = rightbounded(B)
        @test isnull(I.left)
        @test get(I.right) == B
        @test MOD.checkbounds(I, one(T)) == true
        @test MOD.checkbounds(I, convert(T,3)) == false
        @test MOD.checkbounds(I, convert(T,2)) == (isopen ? false : true)

        I = rightbounded(convert(T,2), isopen ? :open : :closed)
        @test isnull(I.left)
        @test get(I.right) == B

        @test_throws ErrorException rightbounded(zero(T), :test)

        I = unbounded(T)
        @test MOD.checkbounds(I, -one(T)) == true
        @test MOD.checkbounds(I, zero(T)) == true
        @test MOD.checkbounds(I,  one(T)) == true
    end
end

info("Testing ", Variable)
for T in (FloatingPointTypes..., IntegerTypes...)

    v = Variable(one(T))
    @test v.value   == one(T)
    @test v.isfixed == false

    @test Variable(one(T), true).value   == one(T)
    @test Variable(one(T), true).isfixed == false

    @test fixed(one(T)).value   == one(T)
    @test fixed(one(T)).isfixed == true

    if T in FloatingPointTypes
        @test eltype(convert(Variable{Float32}, v)) == Float32
        @test eltype(convert(Variable{Float64}, v)) == Float64
    else
        @test eltype(convert(Variable{Int32}, v)) == Int32
        @test eltype(convert(Variable{Int64}, v)) == Int64
    end
end

info("Testing ", HyperParameter)
for T in (FloatingPointTypes..., IntegerTypes...)
    I = rightbounded(one(T), :open)

    for fixedvar in (true, false)
        P = HyperParameter(zero(T), I, fixedvar)
        @test P.value == zero(T)
        @test P.isfixed == fixedvar
    end

    @test_throws ErrorException HyperParameter(one(T), I, true)

    show(DevNull, HyperParameter(zero(T), I, true))

    P = HyperParameter(convert(T,2), unbounded(T), false)

    @test_approx_eq 3*P 3*2
    @test_approx_eq P*3 3*2

    @test_approx_eq 3/P 3/2
    @test_approx_eq P/3 2/3

    if !(T <: Unsigned)
        @test_approx_eq -P  -2
        @test_approx_eq 3-P 3-2
        @test_approx_eq P-3 2-3
    end

    @test_approx_eq 3^P 3^2
    if !(T <: Integer)
        @test_approx_eq P^convert(T,0.5) 2^convert(T,0.5)
    end
    @test_approx_eq P^3 2^3

    @test_approx_eq besselk(P, 1) besselk(2,1)
    @test_approx_eq exp(P) exp(2)
    @test_approx_eq gamma(P) gamma(2)
    @test_approx_eq tanh(P) tanh(2)
end


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

# run download lib svm before running this code
include("../src/read-libsvm/read-libsvm.jl")

using Test
@testset "load_libsvm_file" begin
  # Checks that E2006.train is loaded and there are no errors in the data.
  file_name = joinpath(@__DIR__, "../data/E2006.train")
  result = load_libsvm_file(file_name)
  @test length(result.target_values) == 16087
  @test size(result.coefficient_matrix, 1) == 16087
  @test size(result.coefficient_matrix, 2) == 150360
  @test result.target_values[1] == -3.58943963121738
  @test result.target_values[2] == -3.52663816315417
  @test result.coefficient_matrix[1,376] == 0.000187741243486602
  @test result.coefficient_matrix[1,435] == 0.00012679437327207
  @test result.coefficient_matrix[1,435] == 0.00012679437327207
  @test result.coefficient_matrix[2,667] == 0.000139063750019391
  preprocessed_data = preprocess_regression_problem(result)
  @test size(preprocessed_data.coefficient_matrix, 2) == 150349
end


text:
include("read_line_list.jl")
include("extract_orders.jl")
include("ccf_total.jl")
include("rvs_from_ccf_total.jl")
include("ccf_orders.jl")


text:
module TablesTests

using InplaceArrays.Arrays

data = Float64[2,3,1,3,6,7,3,2,5,6,3,4]
ptrs = [1,4,4,7,13]
a = Table(data,ptrs)
b = [ data[ptrs[i]:ptrs[i+1]-1] for i in 1:length(ptrs)-1]
test_array(a,b)

data = Float64[]
ptrs = [1,]
a = Table(data,ptrs)
b = [ data[ptrs[i]:ptrs[i+1]-1] for i in 1:length(ptrs)-1]
test_array(a,b)

end # module


text:
using Test
using HTTP, HTTP.IOExtras
using BufferedStreams

# For more information see: https://github.com/JuliaWeb/HTTP.jl/pull/288
@testset "Chunking" begin
    sz = 90
    port = 8095
    hex(n) = string(n, base=16)
    encoded_data = "$(hex(sz + 9))\r\n" * "data: 1$(repeat("x", sz))\n\n" * "\r\n" *
                   "$(hex(sz + 9))\r\n" * "data: 2$(repeat("x", sz))\n\n" * "\r\n" *
                   "$(hex(sz + 9))\r\n" * "data: 3$(repeat("x", sz))\n\n" * "\r\n"
    decoded_data = "data: 1$(repeat("x", sz))\n\n" *
                   "data: 2$(repeat("x", sz))\n\n" *
                   "data: 3$(repeat("x", sz))\n\n"
    split1 = 106
    split2 = 300

    t = @async HTTP.listen("127.0.0.1", port) do http
        startwrite(http)
        tcp = http.stream.c.io

        write(tcp, encoded_data[1:split1])
        flush(tcp)
        sleep(1)

        write(tcp, encoded_data[split1+1:split2])
        flush(tcp)
        sleep(1)

        write(tcp, encoded_data[split2+1:end])
        flush(tcp)
    end

    sleep(1)
    @assert !istaskdone(t)

    r = HTTP.get("http://127.0.0.1:$port")

    @test String(r.body) == decoded_data

    for wrap in (identity, BufferedInputStream)
        r = ""

        # Ignore byte-by-byte read warning
        CL = Base.CoreLogging
        CL.with_logger(CL.SimpleLogger(stderr, CL.Error)) do
            HTTP.open("GET", "http://127.0.0.1:$port") do io
                io = wrap(io)
                x = split(decoded_data, "\n")

                for i in 1:6
                    l = readline(io)
                    @test l == x[i]
                    r *= l * "\n"
                end
            end
        end

        @test r == decoded_data
    end
end


text:
#module ListItem

#export ListIt, Node, EmptyNode

abstract type Node end

mutable struct EmptyNode <: Node
  next::Node
  previous::Node
  EmptyNode() = (x=new(); x.next=x; x.previous=x; x)
end
# aa=EmptyNode()

mutable struct ListIt <: Node
    next           :: Node
    previous       :: Node
    bucketSizeRow  :: Int
    MAXBUCKETS     :: Int
    bucketTotal    :: Array{Float64,1}
    bucketVariance :: Array{Float64,1}

    function ListIt()
        x = new()
        x.next = x
        x.previous = x
        x.bucketSizeRow = 1
        x.MAXBUCKETS = 1000
        x.bucketTotal = zeros(x.MAXBUCKETS+1)
        x.bucketVariance = zeros(x.MAXBUCKETS+1)
        x
    end
end

function insertBucket!(lst_it::Node, value::Float64, variance::Float64)
    k = lst_it.bucketSizeRow
    lst_it.bucketSizeRow += 1
    setTotal!(lst_it, value, k)
    setVariance!(lst_it, variance, k)
end

function removeBucket(lst_it::ListIt)
    # Removes the first Bucket
    compressBucketsRow(lst_it, 1);
end

function compressBucketsRow!(lst_it::ListIt, numItemsDeleted::Int)
    for k = numItemsDeleted : lst_it.MAXBUCKETS
        lst_it.bucketTotal[k-numItemsDeleted] = lst_it.bucketTotal[k]
        lst_it.bucketVariance[k-numItemsDeleted]=lst_it.bucketVariance[k]
    end
    for k = 1 : numItemsDeleted
        clearBucket!(lst_it, lst_it.MAXBUCKETS - k + 1)
    end
    lst_it.bucketSizeRow -= numItemsDeleted
end

function measureByteSize(lst_it::ListIt)
    # int size = (int) SizeOf.sizeOf(this); return size;
end

function total(lst_it::ListIt, k::Int)
    lst_it.bucketTotal[k]
end

function variance(lst_it::ListIt, k::Int)
    lst_it.bucketVariance[k]
end

function setTotal!(lst_it::ListIt, val::Float64, k::Int)
    lst_it.bucketTotal[k] = val
end

function setVariance!(lst_it::ListIt, val::Float64, k::Int)
    lst_it.bucketVariance[k] = val
end

function clearBucket!(lst_it::ListIt, k::Int)
    setTotal!(lst_it, 0, k)
    setVariance!(lst_it, 0, k)
end

function clear!(lst_it::ListIt)
    lst_it.bucketSizeRow = 0
    for k = 1 : lst_it.MAXBUCKETS
        clearBucket!(lst_it, k)
    end
end
#end # module
###############################################################################

#using ListItem

""" Pointers to the head and tail of the linked list  """
mutable struct List
    count :: Int
    head  :: Node
    tail  :: Node
    function List()
        x = new()
        x.count = 0
        x.head = ListIt() # EmptyNode()
        x.tail = ListIt() # EmptyNode()
        x
    end
end

function clearLst!(lst::List)
  # public void clear()
  # post: clears the list so that it contains no elements.
  lst.head = EmptyNode()
  lst.tail = EmptyNode()
  lst.count = 0
end

function addToHead!(lst::List)
  a = ListIt()
  lst.head.previous = a
  a.next = lst.head
  lst.head = a
  if typeof(lst.tail) <: EmptyNode
    println("Yep: typeof(lst.tail) <: EmptyNode")
    lst.tail = lst.head
    println("Nope: typeof(lst.tail) <: EmptyNode")
  end
  lst.count += 1
end

function removeFromHead!(lst::List)
  # pre: list is not empty
  # post: removes and returns first object from the list
  lst.head = lst.head.next
  if !(typeof(lst.head) <: EmptyNode)
    lst.head.previous = lst.head.previous #EmptyNode()
  else
    lst.tail = EmptyNode()
  end
  lst.count -= 1
  lst.head
end

function addToTail!(lst::List)
  # line 94:
  # public void addToTail() {
  #  pre: anObject is non-null
  # post: the object is added at the end of the list
  b = ListIt()
  lst.tail.next = b
  b.previous = lst.tail
  lst.tail = b
  if typeof(lst.head) <: EmptyNode
    lst.head = lst.tail
  end
  lst.count += 1
end

function removeFromTail!(lst::List)
  #   public void removeFromTail() {
  # //			pre: list is not empty
  # //			post: the last object in the list is removed and returned
  # //			ListItem temp;
  # //			temp = this.tail;
  lst.tail = lst.tail.previous
  #if typeof(lst.tail) <: EmptyNode
  #  lst.head = EmptyNode()
  #else
  lst.tail.next = lst.tail.next
  #end
  lst.count -= 1
end
function myCounter(lst_it::Node, counter=1)

  # println("... myCounter: counter = $counter")
  if typeof(lst_it) <: EmptyNode || lst_it.next == lst_it
    return counter
  else
    if lst_it == lst_it.next.previous
      myCounter(lst_it.next, counter+1)
    else
      return counter
    end
  end
end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
# area plots

# Vega
using Vega, KernelDensity, Distributions

x = rand(Beta(3.0, 2.0), 10)
k = kde(x)

areaplot(x = k.x, y = k.density)

# stacked area

x = [0,1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6,7,8,9]
y = [28, 43, 81, 19, 52, 24, 87, 17, 68, 49, 55, 91, 53, 87, 48, 49, 66, 27, 16, 15]
g = [0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1]

a = areaplot(x = x, y = y, group = g, stacked = true)

# Plotly
using PlotlyJS

x = ["Winter", "Spring", "Summer", "Fall"]

plot([
    scatter(
        x=x, y=[40, 60, 40, 10],
        stackgroup="one", mode="lines", hoverinfo="x+y",
        line=attr(width=0.5, color="rgb(131, 90, 241)")
    ),
    scatter(
        x=x, y=[20, 10, 10, 60],
        stackgroup="one", mode="lines", hoverinfo="x+y",
        line=attr(width=0.5, color="rgb(111, 231, 219)")
    ),
    scatter(
        x=x, y=[40, 30, 50, 30],
        stackgroup="one", mode="lines", hoverinfo="x+y",
        line=attr(width=0.5, color="rgb(184, 247, 2121)")
    ),
], Layout(yaxis_range=(0, 100)))

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
using Bebop
using Base.Test

if length(ARGS) < 1
    error("missing argument")
end

const duration = 5

openFile(ARGS[1]) do file
    sampRate = file.sfinfo.samplerate
    channels = file.sfinfo.channels
    frames = readFrames(file, sampRate * duration)
    letsjam(sampRate, channels, Int16) do
        play(frames)
    end
end


text:
using Test

function add!(a,b,c)
    local tid = 1
    while (tid <= min(length(a), length(b), length(c)))
        c[tid] = a[tid] + b[tid]
        tid += 1
    end
end

function main()

    # Make three vectors
    a = Vector{Any}(fill(0, 10))
    b = Vector{Any}(fill(0, 10))
    c = Vector{Any}(fill(0, 10))

    # Fill a and b with values
    for i in 1:10
        a[i] = i
        b[i] = i * 2
    end

    # Fill c with values
    add!(a,b,c)

    # Do a sanity check
    for i in 1:length(a)
        @test a[i] + b[i] ≈ c[i]
    end
end

main()


text:
function test_hello(name::AbstractString)
    num_greetings = length(GREETINGS)
    rand_greeting = floor(Int, (num_greetings-1)*rand()) + 1
    string(GREETINGS[rand_greeting], " ", name)
end

function test_exception()
    ex = thriftbuild(InvalidOperation, Dict(:oper => "test_exception"))
    throw(ex)
end

function test_oneway()
    println("server received oneway method call")
end

function ping()
    println("server received ping method call")
end

function test_enum(enum_val::Int32)
    if enum_val == TestEnum.ONE
        return TestEnum.TEN
    elseif enum_val == TestEnum.TWO
        return TestEnum.TWENTY
    end
    ex = thriftbuild(InvalidOperation, Dict(:oper => "test_enum"))
    throw(ex)
end

function _test_types(types)
    set_field!(types, :bool_val,     !get_field(types, :bool_val))
    set_field!(types, :byte_val,     UInt8(get_field(types, :byte_val) + 1))
    set_field!(types, :i16_val,      Int16(get_field(types, :i16_val) + 1))
    set_field!(types, :i32_val,      Int32(get_field(types, :i32_val) + 1))
    set_field!(types, :i64_val,      Int64(get_field(types, :i64_val) + 1))
    set_field!(types, :double_val,   -(get_field(types, :double_val)))
    set_field!(types, :string_val,   uppercase(get_field(types, :string_val)))

    d = Dict{Int32,Int16}()
    for (k,v) in get_field(types, :map_val)
        d[k] = Int16(2*v)
    end
    set_field!(types, :map_val, d)

    l = Int16[]
    for v in get_field(types, :list_val)
        push!(l, Int16(v+10))
    end
    set_field!(types, :list_val, l)

    s = Set{UInt8}()
    for v in get_field(types, :set_val)
        push!(s, UInt8(v+10))
    end
    set_field!(types, :set_val, s)
    types
end

test_types(types::AllTypes) = _test_types(types)
test_types_default(types::AllTypesDefault) = _test_types(types)


text:
"""
`row` and `col` must be in [1,size(mat)], i.e indexing starts at 1 (Julia).

# Implementation
For some unkwnown reason, calling `MatSetValue` fails.
"""
function Base.setindex!(mat::PetscMat, value::Number, row::Integer, col::Integer)
    MatSetValues(mat, PetscInt[row - 1], PetscInt[col - 1], PetscScalar[value], INSERT_VALUES)
end

# This is stupid but I don't know how to do better yet
Base.setindex!(mat::PetscMat, values, row::Integer, cols) = MatSetValues(mat, [row-1], collect(cols) .- 1, values, INSERT_VALUES)
Base.setindex!(mat::PetscMat, values, rows, col::Integer) = MatSetValues(mat, collect(rows) .- 1, [col-1], values, INSERT_VALUES)

Base.ndims(::Type{PetscMat}) = 2

"""
    create_matrix(nrows, ncols, nrows_loc = PETSC_DECIDE, ncols_loc = PETSC_DECIDE; auto_setup = false)

Create a `PetscMat` matrix of global size `(nrows, ncols)`.

Use `auto_setup = true` to immediatly call `set_from_options!` and `set_up!`.
"""
function create_matrix(nrows, ncols, nrows_loc = PETSC_DECIDE, ncols_loc = PETSC_DECIDE; auto_setup = false, comm::MPI.Comm = MPI.COMM_WORLD)
    mat = MatCreate()
    MatSetSizes(mat::PetscMat, nrows_loc, ncols_loc, nrows, ncols)

    if (auto_setup)
        set_from_options!(mat)
        set_up!(mat)
    end
    return mat
end

"""
Wrapper to `MatCreateComposite` using the "alternative construction" from the PETSc documentation.
"""
function create_composite_add(matrices)
    N, M = MatGetSize(matrices[1])
    n, m = MatGetLocalSize(matrices[1])
    mat = create_matrix(N, M, n, m; auto_setup = false, comm = matrices[1].comm)
    MatSetType(mat, "composite")
    for m in matrices
        MatCompositeAddMat(mat, m)
    end
    assemble!(mat)
    return mat
end

set_global_size!(mat::PetscMat, nrows, ncols) = MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, nrows, ncols)
set_local_size!(mat::PetscMat, nrows, ncols) = MatSetSizes(mat, nrows, ncols, PETSC_DECIDE, PETSC_DECIDE)

set_from_options!(mat::PetscMat) = MatSetFromOptions(mat)

set_up!(mat::PetscMat) = MatSetUp(mat)

"""
    get_range(mat::PetscMat)

Wrapper to `MatGetOwnershipRange`

However, the result `(rstart, rend)` is such that `mat[rstart:rend]` are the rows handled by the local processor.
This is different from the default `PETSc.MatGetOwnershipRange` result where the indexing starts at zero and where
`rend-1` is last row handled by the local processor.
"""
function get_range(mat::PetscMat)
    rstart, rend = MatGetOwnershipRange(mat)
    return (rstart + 1, rend)
end

"""
    get_urange(mat::PetscMat)

Provide a `UnitRange` from the method `get_range`.
"""
function get_urange(mat::PetscMat)
    rstart, rend = MatGetOwnershipRange(mat)
    return rstart+1:rend
end

"""
    Wrapper to `MatAssemblyBegin` and `MatAssemblyEnd` successively.
"""
function assemble!(mat::PetscMat, type::MatAssemblyType = MAT_FINAL_ASSEMBLY)
    MatAssemblyBegin(mat, type)
    MatAssemblyEnd(mat, type)
end

"""
    set_values!(mat::PetscMat, I, J, V, mode = ADD_VALUES)

Set values of `mat` in `SparseArrays` fashion : using COO format:
`mat[I[k], J[k]] = V[k]`.
"""
function set_values!(mat::PetscMat, I::Vector{PetscInt}, J::Vector{PetscInt}, V::Vector{PetscScalar}, mode = ADD_VALUES)
    for (i,j,v) in zip(I, J, V)
        MatSetValue(mat, i - 1, j - 1, v, mode)
    end
end

set_values!(mat, I, J, V, mode = ADD_VALUES) = set_values!(mat, PetscInt.(I), PetscInt.(J), PetscScalar.(V), mode)

function preallocate_MPIAIJ(mat::PetscMat; dnz::PetscInt = 0, dnnz::Vector{PetscInt} = C_NULL, onz::PetscInt = 0, onnz::Vector{PetscInt} = C_NULL)
    MatMPIAIJSetPreallocation(mat, dnz, dnnz, onz, onnz)
end

"""
    mat2file(mat::PetscMat, filename::String, format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, type::String = "ascii")

Write a PetscMat to a file.
"""
function mat2file(mat::PetscMat, filename::String, format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, type::String = "ascii")
    viewer = PetscViewer(mat.comm, filename, format, type)
    MatView(mat, viewer)
    destroy!(viewer)
end

Base.show(::IO, mat::PetscMat) = MatView(mat)

destroy!(mat::PetscMat) = MatDestroy(mat)

max_stars_repo_name: bmxam/JuliaPetsc

max_issues_repo_name: bmxam/JuliaPetsc

max_forks_repo_name: bmxam/JuliaPetsc


text:
using COSMO
using Plots
using PortfolioOpt
using PortfolioOpt.TestUtils:
    sequential_backtest_market, get_test_data, mean_variance,
    percentchange, timestamp, rename!
using Statistics

include("./examples/example_utils.jl")

############ Read Prices #############
Prices = get_test_data()
numD, numA = size(Prices) # A: Assets    D: Days

############# Calculating returns #####################
returns_series = percentchange(Prices)
returns = values(returns_series)
# risk free asset
rf = fill(0.0, numD) # fill(3.2e-4, numD) # readcsv(".\\rf.csv")

############ Efficient frontier ###################

DEFAULT_SOLVER = optimizer_with_attributes(
    COSMO.Optimizer, "verbose" => false, "max_iter" => 9000000, "eps_abs" => 1e-6, "eps_rel" => 1e-6
)

# General Parameters
t = 120
k_back = 60
Σ, r̄ = mean_variance(returns[(t - k_back):(t - 1), :])

# Parameter Strategies
δ = 0.028 * 6.05
Δ = std(returns[(end - k_back):end, :]; dims=1)'[:,1] / 5

formulation_markowitz = MeanVariance(;
    predicted_mean = r̄,
    predicted_covariance = Σ,
)
formulation_soyster = RobustBertsimas(;
    predicted_mean = r̄,
    predicted_covariance = Σ,
    uncertainty_delta = Δ,
    bertsimas_budget = 6.0,
)
formulation_bertsimas = RobustBertsimas(;
    predicted_mean = r̄,
    predicted_covariance = Σ,
    uncertainty_delta = Δ,
    bertsimas_budget = 2.8,
)
formulation_bental = RobustBenTal(;
    predicted_mean = r̄,
    predicted_covariance = Σ,
    uncertainty_delta = δ
)

# Loop range and log (no robustness)
range_R = 0:0.00005:0.0015; # Used for `po_min_variance_limit_return`
len = size(range_R, 1)
R_noR = zeros(len);
E_noR = zeros(len);
σ_noR = zeros(len);
R_soy = zeros(len);
x_noR = zeros(len, size(r̄, 1));

# Loop
iter = 0
for R_0 in range_R
    global iter += 1
    model = po_min_variance_limit_return(formulation_markowitz, R_0; rf = rf[t])
    x_noR[iter, :], R_noR[iter] = compute_solution(model, DEFAULT_SOLVER), value.(model[:R])
    E_noR[iter] = sum(x_noR[iter, :]'r̄)
    σ_noR[iter] = sqrt(sum(x_noR[iter, :]'Σ * x_noR[iter, :]))
end

# Loop range and log
range_R = 0:0.00005:0.0006; # Used for `po_min_variance_limit_return`
len = size(range_R, 1)
R_ber = zeros(len);
E_ber = zeros(len);
σ_ber = zeros(len);
R_soy = zeros(len);
E_soy = zeros(len);
σ_soy = zeros(len);
R_ben_tal = zeros(len);
E_ben_tal = zeros(len);
σ_ben_tal = zeros(len);
x_soy = zeros(len, size(r̄, 1));
x_ber = zeros(len, size(r̄, 1));
x_ben_tal = zeros(len, size(r̄, 1));

# Loop robust
iter = 0
for R_0 in range_R
    global iter += 1
    model = po_min_variance_limit_return(formulation_markowitz, R_0; rf = rf[t])
    x_noR[iter, :], R_noR[iter] = compute_solution(model, DEFAULT_SOLVER), value.(model[:R])
    E_noR[iter] = sum(x_noR[iter, :]'r̄)
    σ_noR[iter] = sqrt(sum(x_noR[iter, :]'Σ * x_noR[iter, :]))

    model = po_min_variance_limit_return(formulation_soyster, R_0; rf = rf[t])
    x_soy[iter, :], R_soy[iter] = compute_solution(model, DEFAULT_SOLVER), value.(model[:R])
    model = po_min_variance_limit_return(formulation_bertsimas, R_0; rf = rf[t])
    x_ber[iter, :], R_ber[iter] = compute_solution(model, DEFAULT_SOLVER), value.(model[:R])
    model = po_min_variance_limit_return(formulation_bental, R_0; rf = rf[t])
    x_ben_tal[iter, :], R_ben_tal[iter] = compute_solution(model, DEFAULT_SOLVER), value.(model[:R])

    E_ber[iter] = sum(x_ber[iter, :]'r̄)
    E_soy[iter] = sum(x_soy[iter, :]'r̄)
    E_ben_tal[iter] = sum(x_ben_tal[iter, :]'r̄)

    σ_ber[iter] = sqrt(sum(x_ber[iter, :]'Σ * x_ber[iter, :]))
    σ_soy[iter] = sqrt(sum(x_soy[iter, :]'Σ * x_soy[iter, :]))
    σ_ben_tal[iter] = sqrt(sum(x_ben_tal[iter, :]'Σ * x_ben_tal[iter, :]))
end

# plot optimistic frontier
plt = plot(
    σ_noR,
    E_noR;
    title="Efficient Frontier Mean-Variance",
    xlabel="Standard Deviation",
    ylabel="Expected Return",
    label="Markowitz",
    legend=:outertopright,
    ylims = (0.0, 0.0015),
    xlims = (0.0, 0.01)
    # markershape=:circle
);
plot!(plt, σ_ber, E_ber; label="Bertsimas");
plot!(plt, σ_soy, E_soy; label="Soyster");
# plot!(plt, σ_ben_tal, E_ben_tal; label="BenTal");
plt

# plot robust frontier
plt = plot(
    σ_noR,
    R_noR;
    title="Problem-Specific Efficient Frontier Mean-Variance",
    xlabel="Standard Deviation",
    ylabel="Robust Return (R)",
    label="Markowitz",
    legend=:outertopright,
    ylims = (0.0, 0.0006),
    xlims = (0.0, 0.015)
);
plot!(plt, σ_ber, R_ber; label="Bertsimas");
plot!(plt, σ_soy, R_soy; label="Soyster");
# plot!(plt, σ_ben_tal, R_ben_tal; label="BenTal");
plt

# plot frontier regret bertsimas
min_return_noR_bertsimas = [worst_case_return(x_noR[i, :], formulation_bertsimas, DEFAULT_SOLVER) for i = 1:length(R_noR)]
max_return_noR_bertsimas = [best_case_return(x_noR[i, :], formulation_bertsimas, DEFAULT_SOLVER) for i = 1:length(R_noR)]
max_return_bertsimas = [best_case_return(x_ber[i, :], formulation_bertsimas, DEFAULT_SOLVER) for i = 1:length(R_ber)]

plt = plot(
    σ_noR,
    R_noR;
    title="Uncertain Efficient Frontier Mean-Variance (Bertsimas)",
    xlabel="Standard Deviation",
    ylabel="Return",
    label="Markowitz",
    legend=:outertopright,
    ribbon=(R_noR-min_return_noR_bertsimas,max_return_noR_bertsimas - R_noR),
    xlims = (0.0, 0.005),
    ylims = (0.0, 0.003),
);
plot!(plt, σ_ber, (R_ber+max_return_bertsimas)/2;
    label="Bertsimas",
    ribbon=((R_ber+max_return_bertsimas)/2-R_ber,max_return_bertsimas - (R_ber+max_return_bertsimas)/2)
);
plt

# plot frontier regret soyster
min_return_noR_soyster = [worst_case_return(x_noR[i, :], formulation_soyster, DEFAULT_SOLVER) for i = 1:length(R_noR)]
max_return_noR_soyster = [best_case_return(x_noR[i, :], formulation_soyster, DEFAULT_SOLVER) for i = 1:length(R_noR)]
max_return_bertsimas_soyester = [best_case_return(x_ber[i, :], formulation_soyster, DEFAULT_SOLVER) for i = 1:length(R_ber)]
min_return_bertsimas_soyester = [worst_case_return(x_ber[i, :], formulation_soyster, DEFAULT_SOLVER) for i = 1:length(R_ber)]
max_return_soyester = [best_case_return(x_soy[i, :], formulation_soyster, DEFAULT_SOLVER) for i = 1:length(R_soy)]

plt = plot(
    σ_noR,
    R_noR;
    title="Uncertain Efficient Frontier Mean-Variance (Soyster)",
    xlabel="Standard Deviation",
    ylabel="Return",
    label="Markowitz",
    legend=:outertopright,
    ribbon=(R_noR-min_return_noR_soyster,max_return_noR_soyster - R_noR),
    xlims = (0.0, 0.005),
    ylims = (-0.001, 0.004),
);
plot!(plt, σ_ber, (min_return_bertsimas_soyester+max_return_bertsimas_soyester)/2;
    label="Bertsimas",
    ribbon=((min_return_bertsimas_soyester+max_return_bertsimas_soyester)/2-min_return_bertsimas_soyester,max_return_bertsimas_soyester - (min_return_bertsimas_soyester+max_return_bertsimas_soyester)/2)
);
plot!(plt, σ_soy, (R_soy+max_return_soyester)/2;
    label="Soyster",
    ribbon=((R_soy+max_return_soyester)/2-R_soy,max_return_soyester - (R_soy+max_return_soyester)/2)
);
plt

# plot frontier regret bental
min_return_noR_bental = [worst_case_return(x_noR[i, :], formulation_bental, DEFAULT_SOLVER) for i = 1:length(R_noR)]
max_return_noR_bental = [best_case_return(x_noR[i, :], formulation_bental, DEFAULT_SOLVER) for i = 1:length(R_noR)]
max_return_bental = [best_case_return(x_ben_tal[i, :], formulation_bental, DEFAULT_SOLVER) for i = 1:length(R_ben_tal)]

plt = plot(
    σ_noR,
    R_noR;
    title="Uncertain Efficient Frontier Mean-Variance (Ben-Tal)",
    xlabel="Standard Deviation",
    ylabel="Return",
    label="Markowitz",
    legend=:outertopright,
    ribbon=(R_noR-min_return_noR_bental,max_return_noR_bental - R_noR),
);
plot!(plt, σ_ben_tal, (R_ben_tal+max_return_bental)/2;
    label="Ben-Tal",
    ribbon=((R_ben_tal+max_return_bental)/2-R_ben_tal,max_return_bental - (R_ben_tal+max_return_bental)/2)
);
plt

max_stars_repo_name: andrewrosemberg/PortfolioOpt

max_issues_repo_name: andrewrosemberg/PortfolioOpt

max_forks_repo_name: andrewrosemberg/PortfolioOpt


text:
# src/runtests.jl
# Test Suite for {PROJECT_NAME}

using {PROJECT_NAME}

max_stars_repo_name: JuliaTagBot/julia-template

max_issues_repo_name: JuliaTagBot/julia-template

max_forks_repo_name: JuliaTagBot/julia-template