text:
#####
##### Momentum advection operators
#####

# Alternate names for advective fluxes
@inline _advective_momentum_flux_Uu(i, j, k, grid::APG, args...) = advective_momentum_flux_Uu(i, j, k, grid, args...)
@inline _advective_momentum_flux_Vu(i, j, k, grid::APG, args...) = advective_momentum_flux_Vu(i, j, k, grid, args...)
@inline _advective_momentum_flux_Wu(i, j, k, grid::APG, args...) = advective_momentum_flux_Wu(i, j, k, grid, args...)

@inline _advective_momentum_flux_Uv(i, j, k, grid::APG, args...) = advective_momentum_flux_Uv(i, j, k, grid, args...)
@inline _advective_momentum_flux_Vv(i, j, k, grid::APG, args...) = advective_momentum_flux_Vv(i, j, k, grid, args...)
@inline _advective_momentum_flux_Wv(i, j, k, grid::APG, args...) = advective_momentum_flux_Wv(i, j, k, grid, args...)

@inline _advective_momentum_flux_Uw(i, j, k, grid::APG, args...) = advective_momentum_flux_Uw(i, j, k, grid, args...)
@inline _advective_momentum_flux_Vw(i, j, k, grid::APG, args...) = advective_momentum_flux_Vw(i, j, k, grid, args...)
@inline _advective_momentum_flux_Ww(i, j, k, grid::APG, args...) = advective_momentum_flux_Ww(i, j, k, grid, args...)

@inline _advective_tracer_flux_x(i, j, k, grid::APG, args...) = _advective_tracer_flux_x(i, j, k, grid, args...)
@inline _advective_tracer_flux_y(i, j, k, grid::APG, args...) = _advective_tracer_flux_y(i, j, k, grid, args...)
@inline _advective_tracer_flux_z(i, j, k, grid::APG, args...) = _advective_tracer_flux_z(i, j, k, grid, args...)

"""
    div_Uu(i, j, k, grid, advection, U, u)

Calculate the advection of momentum in the x-direction using the conservative form, ∇·(Uu)

    1/Vᵘ * [δxᶠᵃᵃ(ℑxᶜᵃᵃ(Ax * u) * ℑxᶜᵃᵃ(u)) + δy_fca(ℑxᶠᵃᵃ(Ay * v) * ℑyᵃᶠᵃ(u)) + δz_fac(ℑxᶠᵃᵃ(Az * w) * ℑzᵃᵃᶠ(u))]

which will end up at the location `fcc`.
"""
@inline function div_Uu(i, j, k, grid, advection, U, u)
    return 1/Vᶠᶜᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, _advective_momentum_flux_Uu, advection, U[1], u) +
                                    δyᵃᶜᵃ(i, j, k, grid, _advective_momentum_flux_Vu, advection, U[2], u) +
                                    δzᵃᵃᶜ(i, j, k, grid, _advective_momentum_flux_Wu, advection, U[3], u))
end

"""
    div_Uv(i, j, k, grid, advection, U, v)

Calculate the advection of momentum in the y-direction using the conservative form, ∇·(Uv)

    1/Vʸ * [δx_cfa(ℑyᵃᶠᵃ(Ax * u) * ℑxᶠᵃᵃ(v)) + δyᵃᶠᵃ(ℑyᵃᶜᵃ(Ay * v) * ℑyᵃᶜᵃ(v)) + δz_afc(ℑxᶠᵃᵃ(Az * w) * ℑzᵃᵃᶠ(w))]

which will end up at the location `cfc`.
"""
@inline function div_Uv(i, j, k, grid, advection, U, v)
    return 1/Vᶜᶠᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, _advective_momentum_flux_Uv, advection, U[1], v) +
                                    δyᵃᶠᵃ(i, j, k, grid, _advective_momentum_flux_Vv, advection, U[2], v)    +
                                    δzᵃᵃᶜ(i, j, k, grid, _advective_momentum_flux_Wv, advection, U[3], v))
end

"""
    div_Uw(i, j, k, grid, advection, U, w)

Calculate the advection of momentum in the z-direction using the conservative form, ∇·(Uw)

    1/Vʷ * [δx_caf(ℑzᵃᵃᶠ(Ax * u) * ℑxᶠᵃᵃ(w)) + δy_acf(ℑzᵃᵃᶠ(Ay * v) * ℑyᵃᶠᵃ(w)) + δzᵃᵃᶠ(ℑzᵃᵃᶜ(Az * w) * ℑzᵃᵃᶜ(w))]

which will end up at the location `ccf`.
"""
@inline function div_Uw(i, j, k, grid, advection, U, w)
    return 1/Vᶜᶜᶠ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, _advective_momentum_flux_Uw, advection, U[1], w) +
                                    δyᵃᶜᵃ(i, j, k, grid, _advective_momentum_flux_Vw, advection, U[2], w) +
                                    δzᵃᵃᶠ(i, j, k, grid, _advective_momentum_flux_Ww, advection, U[3], w))
end

@inline div_Uu(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U, u) where FT = zero(FT)
@inline div_Uv(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U, v) where FT = zero(FT)
@inline div_Uw(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U, w) where FT = zero(FT)


text:
"""
Download seismic data from server

# Seisdata array structure

S[network_id][station_id][cnannel_id] gives a SeisData structure

May 29, 2019
Kurama Okubo
"""

using SeisIO, Noise, Printf, Dates, PlotlyJS, JLD2, FileIO, MPI, ProgressMeter

# Self-defined functions
include("./lib/utils.jl")
using .Utils

#------------------------------------------------------------------#
#For the time being, we need remove_response function from obspy
#This will be replaced by SeisIO modules in the near future.
#Please activate obspy enviroment before launching Julia.
include("./lib/remove_response_obspy.jl")
using .Remove_response_obspy
#------------------------------------------------------------------#

#---parameters---#
Network     = ["BP"]
#Station    = ["LCCB", "MMNB", "VCAB", "CCRB"]
Station     = ["LCCB", "MMNB"]
#Channels   = ["BP1", "BP2", "BP3"]
Channels    = ["BP1", "BP2", "BP3"]
src         = "NCEDC" #Data servise center

Starttime   = DateTime(2004,9,28,0,0,0)
Endtime     = DateTime(2004,9,28,24,0,0)
CC_time_unit = 3600 # minimum time unit for cross-correlation [s]

pre_filt    = (0.001, 0.002, 10.0, 20.0) #prefilter of remove_response: taper between f1 and f2, f3 and f4
downsample_fs = 20; #downsampling rate after filtering
MAXMPINUM   = 32 #Limit np for parallel downloading

IsRemoveStationXML = true #Removing all StationXML. (instrumental response is automatically removed before removing XML)
foname      = "./data/SeisData_BP" #output name
#----------------#

MPI.Init()

# establish the MPI communicator and obtain rank
comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

# print initial logo and info
if rank == 0 Utils.initlogo() end

#-------------------------------------------------------------------#
#NEVER CHANGE THIS THRESHOLD OTHERWISE IT OVERLOADS THE DATA SERVER
if size > MAXMPINUM throw(DomainError(size, "np must be smaller than $MAXMPINUM.")) end
#-------------------------------------------------------------------#

# calculate start time with each CC_time_unit
stlist = get_starttimelist(Starttime, Endtime, CC_time_unit)

S = Array{Any, 1}(undef, length(stlist));

mpiitrcount = 0

# progress bar
if rank == 0 prog = Progress(floor(Int, length(stlist)/size), 1.0) end

for stid = 1:length(stlist) #to be parallelized

    processID = stid - (size * mpiitrcount)

    if rank == processID-1

        S[stid] = Array{SeisData, 1}(undef, length(Station));

        for networkid = 1:length(Network)
        	for stationid = 1:length(Station)

                ns_id = (networkid-1)*(length(Station)) + stationid
                S[stid][ns_id] = SeisData(length(Channels)) #S[timeid][sta,netid][cha id]

                for channelsid = 1:length(Channels)

                    requeststr = @sprintf("%s.%s..%s", Network[networkid], Station[stationid], Channels[channelsid])

                    #---download data---#
                    Stemp = get_data("NCEDC", requeststr, s=stlist[stid],t=CC_time_unit, v=0, src=src, w=false, xf="$requeststr.$stid.xml")

                    #---remove response---#
                    Remove_response_obspy.remove_response_obspy!(Stemp, "$requeststr.$stid.xml", pre_filt=pre_filt, output="VEL")
                    if IsRemoveStationXML rm("$requeststr.$stid.xml") end

                    #---check for gaps---#
                    SeisIO.ungap!(Stemp)

                    #---remove earthquakes---#
                    # NOT IMPLEMENTED YET!

                    #---detrend---#
                    SeisIO.detrend!(Stemp)

                    #---bandpass filter---#
                    SeisIO.filtfilt!(Stemp,fl=0.01,fh=0.9*(0.5*Stemp.fs[1])) #0.9*Nyquist frequency

                    #---taper---#
                    SeisIO.taper!(Stemp,t_max=30.0,α=0.05)

                    #---sync starttime---#
                    SeisIO.sync!(Stemp,s=DateTime(stlist[stid]),t=DateTime(stlist[stid])+Dates.Second(CC_time_unit))

                    #---down sampling---#
                    # Note: filtering should be first then down sampling
                    Stemp = Noise.downsample(Stemp, float(downsample_fs))
                    SeisIO.note!(Stemp, "downsample!, downsample_fs=$downsample_fs")

                    #store Stemp to S
                    S[stid][ns_id][channelsid] = Stemp[1]

                end
            end
        end

        #println("stid:$stid done by rank $rank out of $size processors")
        global mpiitrcount += 1

        #progress bar
        if rank == 0 next!(prog) end

    end
end

#save struct
if rank == 0
    @save foname*".jld2" S
    println("Downloading data is successfully done.")
end

MPI.Finalize()

ext: jl

max_stars_repo_name: kura-okubo/SeisDownloadwJulia

max_issues_repo_name: kura-okubo/SeisDownloadwJulia

max_forks_repo_name: kura-okubo/SeisDownloadwJulia


text:
eta = 0.9

function f1(x)
    (1.0 / 2.0) * (x[1]^2 + eta * x[2]^2)
end

function g1(x, storage)
    storage[1] = x[1]
    storage[2] = eta * x[2]
end

function h1(x, storage)
    storage[1, 1] = 1.0
    storage[1, 2] = 0.0
    storage[2, 1] = 0.0
    storage[2, 2] = eta
end

results = optimize(f1, g1, h1, [127.0, 921.0])
@assert results.converged
@assert norm(results.minimum - [0.0, 0.0]) < 0.01

results = optimize(f1, g1, [127.0, 921.0])
@assert results.converged
@assert norm(results.minimum - [0.0, 0.0]) < 0.01

results = optimize(f1, [127.0, 921.0])
@assert results.converged
@assert norm(results.minimum - [0.0, 0.0]) < 0.01


text:
# small test with one tree quartet for the NLopt
# Claudia January 2015

treal = 0.1
obsCF=[1-2/3*exp(-treal),1/3*exp(-treal),1/3*exp(-treal)]
using NLopt

type ToyQuartet
    t::Float64
end

function objective(x::Float64,obsCF::Vector{Float64})
    length(obsCF) == 3 || error("obsCF must be size 3")
    val = obsCF[1]*log(1-2/3*exp(-x)) + obsCF[2]*log(1/3*exp(-x)) + obsCF[3]*log(1/3*exp(-x))
    return -val
end

function optimizar(obsCF::Vector{Float64}, q::ToyQuartet)
    t = q.t
    k = length(t)
    opt = NLopt.Opt(:LN_BOBYQA,k) # :LD_MMA if use gradient
    # criterion based on prof Bates code
    NLopt.ftol_rel!(opt,1e-12) # relative criterion
    NLopt.ftol_abs!(opt,1e-8) # absolute critetion
    NLopt.xtol_abs!(opt,1e-10) # criterion on parameter value changes
    NLopt.lower_bounds!(opt, zeros(k))
    NLopt.upper_bounds!(opt, Inf)
    count = 0
    function obj(x::Vector{Float64},g::Vector{Float64}) # added g::Vector{Float64} for gradient, ow error
        count += 1
        println("t is $(t) initially, x: $(x)")
        t = deepcopy(x)
        println("t is now $(t), x: $(x)")
        val = objective(x[1],obsCF)
        println("f_$count: $(round(val,5)), x: $(x)")
        return val
    end
    NLopt.min_objective!(opt,obj)
    fmin, xmin, ret = NLopt.optimize(opt,[t])
    return fmin,xmin
end

q = ToyQuartet(1.0)
fmin,xmin = optimizar(obsCF,q)


text:
# ---
# title: 93. Restore IP Addresses
# id: problem93
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Backtracking
# link: <https://leetcode.com/problems/restore-ip-addresses/description/>
# hidden: true
# ---
#
# Given a string `s` containing only digits, return all possible valid IP
# addresses that can be obtained from `s`. You can return them in **any** order.
#
# A **valid IP address** consists of exactly four integers, each integer is
# between `0` and `255`, separated by single dots and cannot have leading zeros.
# For example, "0.1.2.201" and "192.168.1.1" are **valid** IP addresses and
# "0.011.255.245", "192.168.1.312" and "192.168@1.1" are **invalid** IP
# addresses.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "25525511135"
#     Output: ["255.255.11.135","255.255.111.35"]
#
#
# **Example 2:**
#
#
#
#     Input: s = "0000"
#     Output: ["0.0.0.0"]
#
#
# **Example 3:**
#
#
#
#     Input: s = "1111"
#     Output: ["1.1.1.1"]
#
#
# **Example 4:**
#
#
#
#     Input: s = "010010"
#     Output: ["0.10.0.10","0.100.1.0"]
#
#
# **Example 5:**
#
#
#
#     Input: s = "101023"
#     Output: ["1.0.10.23","1.0.102.3","10.1.0.23","10.10.2.3","101.0.2.3"]
#
#
#
#
# **Constraints:**
#
#   * `0 <= s.length <= 3000`
#   * `s` consists of digits only.
#
#
## @lc code=start
using LeetCode

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


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -14.06559781669544 <= q <= 13.37751741699121)

                     add_NL_constraint(m, :(gelu(-0.7556338259975872 + -0.2850076752045205*gelu(0.7432962759699406 + 0.616365769922441*gelu(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + -0.9563519707270385*gelu(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.21036883565748976*gelu(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.5819748667402531*gelu(0.6367931081099356 + 0.08862555157308405*gelu(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.5030828808953012*gelu(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + -0.4316150506981855*gelu(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + 0.4871037750371179*gelu(0.12597161846431115 + 0.10770543162761159*gelu(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.46799585735795457*gelu(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.45029698767672954*gelu(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5])))) + gelu(0.9186311390475046 + -0.8559712370251158*gelu(0.7432962759699406 + 0.616365769922441*gelu(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + -0.9563519707270385*gelu(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.21036883565748976*gelu(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.2512416329037399*gelu(0.6367931081099356 + 0.08862555157308405*gelu(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.5030828808953012*gelu(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + -0.4316150506981855*gelu(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + 0.7592758779705528*gelu(0.12597161846431115 + 0.10770543162761159*gelu(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.46799585735795457*gelu(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.45029698767672954*gelu(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5])))) + gelu(-0.7151630404444922 + -0.15417874070690862*gelu(0.7432962759699406 + 0.616365769922441*gelu(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + -0.9563519707270385*gelu(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.21036883565748976*gelu(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.016895941959414795*gelu(0.6367931081099356 + 0.08862555157308405*gelu(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.5030828808953012*gelu(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + -0.4316150506981855*gelu(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.17047547651198514*gelu(0.12597161846431115 + 0.10770543162761159*gelu(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.46799585735795457*gelu(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.45029698767672954*gelu(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Suffix Array Construction Algorithms

# a wrapper type for a sequence that returns byte-convertible elements
immutable ByteSeq{S} <: AbstractVector{UInt8}
    data::S
end
Base.size(seq::ByteSeq) = (length(seq.data),)
@inline Base.getindex(seq::ByteSeq, i::Integer) = UInt8(seq.data[i])

# SuffixArrays.jl: https://github.com/quinnj/SuffixArrays.jl
function make_sa(T, seq, σ, mmap)
    n = length(seq)
    tmp_sa = mmap ? Mmap.mmap(Vector{Int}, n) : Vector{Int}(n)
    SuffixArrays.sais(ByteSeq(seq), tmp_sa, 0, n, nextpow2(σ), false)
    sa = mmap ? Mmap.mmap(Vector{T}, n) : Vector{T}(n)
    copy!(sa, tmp_sa)
    return sa
end

# pSAscan: https://www.cs.helsinki.fi/group/pads/pSAscan.html
function make_sa_pscan(T, seq, psascan, workdir, mmap)
    seqpath, io = mktemp(workdir)
    sapath = string(seqpath, ".sa5")
    try
        dump_seq(io, seq)
        run(`$psascan -o $sapath $seqpath`)
        return load_sa(T, sapath, mmap)
    catch
        rethrow()
    finally
        rm(seqpath)
        isfile(sapath) && rm(sapath)
    end
end

function dump_seq(io, seq)
    @inbounds for i in 1:length(seq)
        write(io, convert(UInt8, seq[i]))
    end
    close(io)
end

function load_sa(T, file, mmap)
    # load a 40-bit suffix array generated from psascan
    size = filesize(file)
    @assert size % 5 == 0 "file $file is not 40-bit integers"
    n = div(size, 5)
    sa = mmap ? Mmap.mmap(Vector{T}, n) : Vector{T}(n)
    open(file) do input
        load_sa!(input, sa)
    end
    return sa
end

function load_sa!{T}(input::IO, sa::Vector{T})
    # load a suffix array from the `input` into `sa`
    buf = Vector{UInt8}(5)
    i = 0
    while !eof(input)
        read!(input, buf)
        value = T(0)
        @inbounds for j in 1:5
            value |= convert(T, buf[j]) << 8 * (j - 1)
        end
        sa[i+=1] = value
    end
    @assert i == length(sa)
    return sa
end

# other utils

function index_type(n)
    n -= 1
    n ≤ typemax(UInt8)  ? UInt8  :
    n ≤ typemax(UInt16) ? UInt16 :
    n ≤ typemax(UInt32) ? UInt32 : UInt64
end

# suffix array sampling
function sample_sa{T}(sa::Vector{T}, r)
    n = length(sa)
    samples = Vector{T}(cld(n, r))
    sampled = falses(n)
    i′ = 0
    for i in 1:n
        @assert 0 ≤ sa[i] ≤ n - 1
        if sa[i] % r == 0
            samples[i′+=1] = sa[i]
            sampled[i] = true
        end
    end
    return samples, sampled
end

# Burrows-Wheeler Transform
function make_bwt(seq, sa)
    n = length(seq)
    @assert length(sa) == n
    ret = Vector{UInt8}(n)
    j = 1
    for i in 1:n
        # note that `sa` starts from zero
        p = sa[i]
        if p == 0
            ret[1] = seq[end]
        else
            ret[j+=1] = seq[p]
        end
    end
    return ret
end


text:
"""
Ledgers

This package provides support for general financial ledgers

See README.md for the full documentation

Copyright 2019-2020, Eric Forgy, Scott P. Jones and other contributors

Licensed under MIT License, see LICENSE.md
"""
module Ledgers

using UUIDs, StructArrays, AbstractTrees
import Instruments
using Instruments: Position
# Although importing Instruments is enough to extend the methods below, we still import
# them explicitly so we can use them like `Ledgers.position` instead of `Instruments.position`.
# Unfortunately, `Base` also exports `position`.
import Instruments: instrument, currency, symbol, position, amount

export Identifier, AccountId, AbstractAccount, AccountWrapper, AccountNumber
export LedgerAccount, Ledger, Entry, Account, AccountGroup, ledgeraccount, add_account, ledger
export id, balance, credit!, debit!, post!, add_entry!
# Although the methods below are exported by Instruments, we export them explicitly
# because we do not expect users to use `Instruments` directly.
export instrument, currency, symbol, position, amount

include("accounts.jl")

abstract type AbstractLedgerAccount{P} <: AbstractAccount{P} end

struct Entry{P,A<:AbstractLedgerAccount{P}}
    debit::A
    credit::A
    amount::P
end

struct LedgerId <: Identifier
    value::UUID
end

LedgerId() = LedgerId(uuid4())

struct Ledger{P <: Position,A<:AbstractLedgerAccount{P}}
    id::LedgerId
    accounts::Dict{AccountId,A}
    entries::Vector{Entry{P}}
end

function Ledger(accounts::AbstractVector{<:AbstractLedgerAccount{P}}, entries::AbstractVector{Entry{P}}; id=LedgerId()) where {P <: Position}
    ledger = Ledger(P)
    add_account!.(Ref(ledger), accounts)
    add_entry!.(Ref(ledger), entries)
    return ledger
end

Base.iterate(ledger::Ledger, state=1) = iterate(ledger.entries, state)

Ledger(::Type{P}) where {P <: Position} = Ledger{P,LedgerAccount{P}}(LedgerId(),Dict{AccountId,LedgerAccount{P}}(), Vector{Entry{P}}())

mutable struct LedgerAccount{P} <: AbstractLedgerAccount{P}
    id::AccountId
    balance::P
    ledger::Ledger{P}
end

function LedgerAccount(balance::P, ledger) where {P <: Position}
    acc = LedgerAccount{P}(AccountId(), balance, ledger)
    add_account!(ledger, acc)
    acc
end

LedgerAccount(::Type{P}, ledger) where {P <: Position} =
    LedgerAccount(P(0), ledger)

ledgeraccount(acc::LedgerAccount) = acc
ledgeraccount(acc::AbstractAccount) = ledgeraccount(acc.account)

ledger(acc::LedgerAccount) = acc.ledger
ledger(acc::AbstractAccount) = ledger(ledgeraccount(acc))

balance(acc::LedgerAccount) = acc.balance
balance(acc::AbstractAccount) = balance(ledgeraccount(acc))
id(acc::LedgerAccount) = acc.id

function add_account!(ledger::Ledger{P}, acc::AbstractAccount{P}) where {P <: Position}
    acc = ledgeraccount(acc)
    if id(acc) ∈ keys(ledger.accounts)
        warn("Account already in ledger.")
        return
    end
    ledger.accounts[acc.id] = acc
end

function add_accounts!(ledger::Ledger{P}, ag::AccountGroup) where P
    for acc in ag
        add_account!(ledger, acc)
    end
    ledger
end

function add_entry!(ledger::Ledger, entry::Entry)
    (id(entry.credit) ∈ keys(ledger.accounts)
        && id(entry.debit) ∈ keys(ledger.accounts)) || error("Unknown account in entry.")
    push!(ledger.entries, entry)
    post!(entry)
end

Instruments.symbol(::Type{A}) where {P, A<:AbstractAccount{P}} = symbol(P)

Instruments.currency(::Type{A}) where {P, A<:AbstractAccount{P}} = currency(P)

Instruments.instrument(::Type{A}) where {P, A<:AbstractAccount{P}} = instrument(P)

Instruments.position(::Type{A}) where {P, A<:AbstractAccount{P}} = P

Instruments.amount(acc::AbstractAccount) = amount(balance(acc))

debit!(acc::LedgerAccount, amt::Position) = (acc.balance -= amt)

credit!(acc::LedgerAccount, amt::Position) = (acc.balance += amt)

function post!(entry::Entry)
    debit!(ledgeraccount(entry.debit), entry.amount)
    credit!(ledgeraccount(entry.credit), entry.amount)
end

Base.getindex(ledger::Ledger, ix) = ledger.entries[ix]

Base.getindex(ledger::Ledger, id::AccountId) =
    ledger.accounts[id]

Base.getindex(grp::AccountGroup, id::AccountId) =
    grp.accounts[id]

# struct EntityId <: Identifier
#     value::UUID
# end

# EntityId() = EntityId(uuid4())

# struct Entity
#     id::EntityId
#     name::String
#     ledgers::Dict{Type{<:Position},Ledger{<:AbstractAccount}}
# end

# const chartofaccounts = Dict{String,AccountGroup{<:Cash}}()

# iscontra(a::AccountGroup) = !isequal(a,a.parent) && !isequal(a.parent,getledger(a)) && !isequal(a.parent.isdebit,a.isdebit)

# function loadchart(ledgername,ledgercode,csvfile)
#     data, headers = readdlm(csvfile,',',String,header=true)
#     nrow,ncol = size(data)

#     ledger = add(AccountGroup(ledgername,ledgercode))
#     for i = 1:nrow
#         row = data[i,:]
#         number = row[1]
#         name = row[2]
#         parent = chartofaccounts[row[3]]
#         isdebit = isequal(row[4],"Debit")
#         add(AccountGroup(parent,name,number,isdebit))
#     end
#     return ledger
# end

# function trim(a::AccountGroup,newparent::AccountGroup=AccountGroup(a.parent.name,a.parent.number))
#     newaccount = isequal(a,a.parent) ? newparent : AccountGroup(newparent,a.name,a.number,a.isdebit,a.balance)
#     for subaccount in a.accounts
#         balance(subaccount).a_print_account

function _print_account(io::IO, acc)
    iobuff = IOBuffer()
    isempty(acc.number.value) || print(iobuff,"[$(acc.number)] ")
    print(iobuff,acc.name,": ")
    acc.isdebit ? print(iobuff,balance(acc)) : print(iobuff,-balance(acc))
    print(io,String(take!(iobuff)))
end

Base.show(io::IO, id::Identifier) = print(io, id.value)

Base.show(io::IO, number::AccountNumber) = print(io, number.value)

Base.show(io::IO, acc::LedgerAccount) = print(io, "$(string(id(acc))): $(balance(acc))")

Base.show(io::IO, acc::AbstractAccount) = _print_account(io, acc)

Base.show(io::IO, acc::AccountGroup) = print_tree(io, acc)

Base.show(io::IO, entry::Entry) = print_tree(io, entry)

Base.show(io::IO, ledger::Ledger) = print_tree(io, ledger)

Base.show(io::IO, a::Vector{<:AbstractAccount}) = print_tree(io, a)

AbstractTrees.printnode(io::IO, acc::Ledger{P}) where {P <: Position} =
    print(io, "$(symbol(P)) Ledger: [$(acc.id)]")

AbstractTrees.children(acc::Ledger) = collect(values(acc.accounts))

AbstractTrees.printnode(io::IO, acc::AccountGroup) = _print_account(io, acc)

AbstractTrees.children(acc::AccountGroup) = vcat(collect(values(acc.subgroups)), collect(values(acc.accounts)))

AbstractTrees.printnode(io::IO, b::AbstractVector{<:AbstractAccount}) =
    isempty(b) ? print(io, "Accounts: None") : print(io, "Accounts:")

AbstractTrees.children(b::AbstractVector{<:AbstractAccount}) = b

AbstractTrees.children(entry::Entry) = [entry.debit, entry.credit]

AbstractTrees.printnode(io::IO, ::Entry) = print(io, "Entry:")

end # module Ledgers


text:
@testset "test_discrete" begin

C_111 = ss([-5], [2], [3], [0])
C_212 = ss([-5 -3; 2 -9], [1; 2], [1 0; 0 1], [0; 0])
C_221 = ss([-5 -3; 2 -9], [1 0; 0 2], [1 0], [0 0])
C_222_d = ss([-5 -3; 2 -9], [1 0; 0 2], [1 0; 0 1], [1 0; 0 1])

@test c2d_x0map(ss(4*[1 0; 0 1]), 0.5, :zoh) == (ss(4*[1 0; 0 1], 0.5), zeros(0, 2))
@test c2d_x0map(ss(4*[1 0; 0 1]), 0.5, :foh) == (ss(4*[1 0; 0 1], 0.5), zeros(0, 2))
@test_tupleapprox(c2d_x0map(C_111, 0.01, :zoh),
    ss([0.951229424500714], [0.019508230199714396], [3], [0], 0.01), [1 0])
@test_tupleapprox(c2d_x0map(C_111, 0.01, :foh),
    ss([0.951229424500714], [0.01902855227625244], [3], [0.029506188017136226], 0.01),
    [1 -0.009835396005712075])
@test_tupleapprox(c2d_x0map(C_212, 0.01, :zoh),
    ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
    [0.009466805409666932; 0.019219966830212765], [1 0; 0 1], [0; 0], 0.01),
    [1 0 0; 0 1 0])
@test_tupleapprox(c2d_x0map(C_212, 0.01, :foh),
    ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
    [0.008957940478201584; 0.018468989584974498], [1 0; 0 1], [0.004820885889482196;
    0.009738343195298675], 0.01), [1 0 -0.004820885889482196; 0 1 -0.009738343195298675])
@test_tupleapprox(c2d_x0map(C_221, 0.01, :zoh),
    ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
    [0.009753161420545834 -0.0002863560108789034; 9.54520036263011e-5 0.019124514826586465],
    [1.0 0.0], [0.0 0.0], 0.01), [1 0 0 0; 0 1 0 0])
@test_tupleapprox(c2d_x0map(C_221, 0.01, :foh),
    ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
    [0.009511049106772921 -0.0005531086285713394; 0.00018436954285711309 0.018284620042117387],
    [1 0], [0.004917457305816479 -9.657141633428213e-5], 0.01),
    [1 0 -0.004917457305816479 9.657141633428213e-5;
    0 1 -3.219047211142736e-5 -0.009706152723187249])
@test_tupleapprox(c2d_x0map(C_222_d, 0.01, :zoh),
    ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
    [0.009753161420545834 -0.0002863560108789034; 9.54520036263011e-5 0.019124514826586465],
    [1 0; 0 1], [1 0; 0 1], 0.01), [1 0 0 0; 0 1 0 0])
@test_tupleapprox(c2d_x0map(C_222_d, 0.01, :foh),
    ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
    [0.009511049106772921 -0.0005531086285713394; 0.00018436954285711309 0.018284620042117387],
    [1 0; 0 1], [1.0049174573058164 -9.657141633428213e-5; 3.219047211142736e-5 1.0097061527231872], 0.01),
    [1 0 -0.004917457305816479 9.657141633428213e-5; 0 1 -3.219047211142736e-5 -0.009706152723187249])

# Test c2d for transfer functions
G = tf([1, 1], [1, 3, 1])
Gd = c2d(G, 0.2)
@test typeof(Gd) <: TransferFunction{<:Discrete, <:ControlSystems.SisoRational}
@test Gd ≈ tf([0, 0.165883310712090, -0.135903621603238], [1.0, -1.518831946985175, 0.548811636094027], 0.2) rtol=1e-14

# Issue #391
@test c2d(tf(C_111), 0.01, :zoh) ≈ tf(c2d(C_111, 0.01, :zoh))
@test c2d(tf(C_111), 0.01, :foh) ≈ tf(c2d(C_111, 0.01, :foh))

# Issue #471 discussion, test c2d for zero row C
C_210 = ss(C_212.A, C_212.B, zeros(0, 2), zeros(0, 1))
@test c2d(C_210, 0.01).A ≈ c2d(C_212, 0.01).A

# c2d on a zpk model should arguably return a zpk model
Gdzpk = c2d(zpk(G), 0.2)
@test typeof(Gdzpk) <: TransferFunction{<:Discrete, <:ControlSystems.SisoZpk}
@test Gdzpk ≈ zpk([0.8192724211968178], [0.9264518519788194, 0.5923800950063556], 0.16588331071209, 0.2) rtol=1e-14

# ERRORS
@test_throws MethodError c2d(ss([1], [2], [3], [4], 0.01), 0.01)   # Already discrete
@test_throws MethodError c2d(ss([1], [2], [3], [4], -1), 0.01)     # Already discrete

# d2c
@test d2c(c2d(C_111, 0.01)) ≈ C_111
@test d2c(c2d(C_212, 0.01)) ≈ C_212
@test d2c(c2d(C_221, 0.01)) ≈ C_221
@test d2c(c2d(C_222_d, 0.01)) ≈ C_222_d
@test d2c(Gd) ≈ G

sys = ss([0 1; 0 0], [0;1], [1 0], 0)
sysd = c2d(sys, 1)
@test d2c(sysd) ≈ sys

# forward euler / tustin
@test c2d(C_111, 1, :fwdeuler).A == I + C_111.A
for method in (:fwdeuler, :tustin)
    @test d2c(c2d(C_111, 0.01, method), method) ≈ C_111 atol = sqrt(eps())
    @test d2c(c2d(C_212, 0.01, method), method) ≈ C_212 atol = sqrt(eps())
    @test d2c(c2d(C_221, 0.01, method), method) ≈ C_221 atol = sqrt(eps())
    @test d2c(c2d(C_222_d, 0.01, method), method) ≈ C_222_d atol = sqrt(eps())
    @test d2c(c2d(G, 0.01, method), method) ≈ G atol = sqrt(eps())
end

matlab_tustin = let
    A = [0.95094630230333 -0.02800401390866; 0.018669342605773 0.913607617091783]
    B = [0.009754731511517 -0.000280040139087; 9.3346713029e-5 0.019136076170918]
    C = [0.975473151151665 -0.01400200695433; 0.009334671302887 0.956803808545892]
    D = [1.004877365755758 -0.000140020069543; 4.6673356514e-5 1.009568038085459]
    ss(A,B,C,D,0.01)
end

sys_tustin = c2d(C_222_d, 0.01, :tustin)
@test sys_tustin ≈ matlab_tustin atol = 1e-12

matlab_prewarp = let
    A = [0.950906169509825 -0.028025790680969; 0.018683860453979 0.913538448601866]
    B = [0.009762667760265 -0.00028049168885; 9.3497229617e-5 0.019151346602063]
    C = [0.975453084754912 -0.014012895340485; 0.00934193022699 0.956769224300933]
    D = [1.004881333880133 -0.000140245844425; 4.6748614808e-5 1.009575673301031]
    sys = ss(A,B,C,D,0.01)
end

sys_prewarp = c2d(C_222_d, 0.01, :tustin, w_prewarp=10)
@test sys_prewarp ≈ matlab_prewarp atol = 1e-12

Cd = c2d(C_111, 0.001, :fwdeuler)
t = 0:0.001:2
y,_ = step(C_111, t)
yd,_ = step(Cd, t)
@test norm(yd-y) / norm(y) < 1e-3

# dab
import DSP: conv
ζ = 0.2
ω = 1
B = [1]
A   = [1, 2ζ*ω, ω^2]
P  = tf(B,A)
# Control design
ζ0 = 0.7
ω0 = 2
Am = [1, 2ζ0*ω0, ω0^2]
Ao = conv(2Am, [1/2, 1]) # Observer polynomial, add extra pole due to the integrator
AR = [1,0] # Force the controller to contain an integrator ( 1/(s+0) )

B⁺  = [1] # The process numerator polynomial can be facored as B = B⁺B⁻ where B⁻ contains the zeros we do not want to cancel (non-minimum phase and poorly damped zeros)
B⁻  = [1]
Bm  = conv(B⁺, B⁻) # In this case, keep the entire numerator polynomial of the process

R,S,T = rstc(B⁺,B⁻,A,Bm,Am,Ao,AR) # Calculate the 2-DOF controller polynomials

Gcl = tf(conv(B,T),zpconv(A,R,B,S)) # Form the closed loop polynomial from reference to output, the closed-loop characteristic polynomial is AR + BS, the function zpconv takes care of the polynomial multiplication and makes sure the coefficient vectores are of equal length

@test ControlSystems.isstable(Gcl)

p = poles(Gcl)
# Test that all desired poles are in the closed-loop system
@test norm(minimum(abs.((poles(tf(Bm,Am)) .- sort(p, by=imag)')), dims=2)) < 1e-6
# Test that the observer poles are in the closed-loop system
@test norm(minimum(abs.((poles(tf(1,Ao)) .- sort(p, by=imag)')), dims=2)) < 1e-6

end


text:
function readinput(filename)
    map(group -> filter(!isempty, split(group, " ")), readlines(filename))
end

lines = readinput("d02.in")

P1, P2 = 0, 0

for line in lines
    nums = parse.(Int, split(line[1], '-'))
    char = line[2][1]
    pass = line[3]
    x = count(x -> x == char, pass)
    if nums[1] <= x <= nums[2]
        global P1 += 1
    end
    if xor(pass[nums[1]] == char, pass[nums[2]] == char)
        global P2 += 1
    end
end

println("P1: ", P1)
println("P2: ", P2)

max_stars_repo_name: eduellery/adventofcode

max_issues_repo_name: eduellery/adventofcode

max_forks_repo_name: eduellery/adventofcode


text:
using Test
const desired_version = VersionNumber(ARGS[1])
include("../deps/deps.jl")
@testset "pyversion_build ≈ $desired_version" begin
    @test desired_version.major == pyversion_build.major
    @test desired_version.minor == pyversion_build.minor
    if desired_version.patch != 0
        @test desired_version.patch == pyversion_build.patch
    end
end


text:
using StatsBase
using LinearAlgebra

function get_freq(h, vals...)
    for (idx, val) ∈ enumerate(vals)
        if val < h.edges[idx][1] || val >= h.edges[idx][end]
            return eltype(h.weights)(0)
        end
    end
    proj = [searchsortedlast(h.edges[i], val) for (i, val) ∈ enumerate(vals)]
    h.weights[proj...]
end

function get_freq(h, vals::AbstractVector)
    for (idx, val) ∈ enumerate(vals)
        if val < h.edges[idx][1] || val >= h.edges[idx][end]
            return eltype(h.weights)(0)
        end
    end
    proj = [searchsortedlast(h.edges[i], val) for (i, val) ∈ enumerate(vals)]
    h.weights[proj...]
end

struct DPI{H<:Histogram, P}
    state_pdf::H
    policy::P
end

function DPI(states, policy; h_kwargs...)
    state_hist = fit(Histogram, Tuple(getindex.(states, i) for i in 1:length(states[1])); h_kwargs...)
    sh_norm = normalize(state_hist, mode=:probability)
    DPI(sh_norm, policy)
end

function (dpi::DPI)(state, action)
    s_p = get_freq(dpi.state_pdf, state)
    a_p = get(dpi.policy, state_t=state, action_t=action)
    s_p*a_p
end

struct HordeDPI{F}
    dpis::Vector{DPI}
    state_filter::F
end

function (hdpi::HordeDPI)(state, action)
    [dpi(hdpi.state_filter(state),action) for dpi in hdpi.dpis]
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 56807afc-164e-11ec-2cf4-354654176242
begin
	using Pkg
	Pkg.activate("../../")
	using Revise
end

# ╔═╡ 8a11c2b7-b6e5-4c26-9e06-ad3ab958db9c
using TestImages, ImageShow, Colors, ImgProcMic, Statistics, PlutoTest, PlutoUI

# ╔═╡ 3ccbf4fb-0707-4ce8-83d8-b1af555f3fe1
using PoissonRandom

# ╔═╡ 62c6e2f0-0007-40ec-81c5-61cda80f59e5
using ImageFiltering, IndexFunArrays

# ╔═╡ 31cc54e9-d47b-4df7-a785-dad0120e8cdb
md"## Activate our environment"

# ╔═╡ d8c83f2c-76ea-42b9-bd39-7fa5cafcd4d3
md"# Homework 01"

# ╔═╡ 334ebdc5-dc11-4e83-ad12-a1961420a97b
begin
	img = Float64.(testimage("mandril_gray"))
	gray_show(img)
end

# ╔═╡ 2122c12f-0832-4f19-9a77-3047d5311cba
md"## 1. Add Noise To Images

This homework is dedicated to introduce into the basics of Julia.
As examples, we want to apply and remove noise from images
"

# ╔═╡ a7d7a415-7464-43fb-a2bf-8e3f2d999c32
md"## Task 1.1  - Gaussian Noise
In the first step, we want to add Gaussian noise with a certain
standard deviation $\sigma$. Try to get yourself used to the Julia documentation
and check out whether such a function exists already and how you can apply it
to your image/array.

Codewords: **Julia, random, random normal distribution**
"

# ╔═╡ de104b0a-28a4-42de-ae8d-fa857a8f32e0
"""
	add_gauss_noise(img, σ=1)

This function adds normal distributed noise to `img`.
`σ` is an optional argument
"""
function add_gauss_noise(img, σ=one(img))
	# todo
	return similar(img)
end

# ╔═╡ f89bc1e4-d406-4550-988b-71496b64035a
md"
Now we want to use a for loop inside this function.
Try to find out how to write for loops to iterate through an array.
"

# ╔═╡ c0b5ff08-7342-4c25-9791-be3a4918c400
"""
	add_gauss_noise_fl!(img, σ=1)

This function adds normal distributed noise to `img`.
`σ` is an optional argument.
This function is memory efficient by using for loops.
`!` means that the input (`img`) is modified.
Therefore, don't return a new array but instead modify the existing one!
The bang (!) is a convention in Julia that a function modifies the input.
"""
function add_gauss_noise_fl!(img, σ=1)
	# todo
	return img
end

# ╔═╡ d56bc500-f5cf-40f3-8c0b-b4ae26a58367
md"### Tests 1.1
Don't modify but rather take those tests as input whether you are correct.
Green is excellent, red not :(
"

# ╔═╡ 1116455e-4698-489a-add2-4e1d8ecbb546
gray_show(add_gauss_noise(img, 0.15), set_one=true)

# ╔═╡ 8e343a3d-ab0f-45a1-9161-21fec07755c9
md"###### Type Stability of your functions.
That means that the output type is identical to the input type"

# ╔═╡ d865dddc-6d32-416b-ba12-8f6b3870e771
PlutoTest.@test Float32 == eltype(add_gauss_noise(ones(Float32, (512, 512)), 0.3f0))

# ╔═╡ 5c9a8c7b-e9c4-4177-86bb-b69113fce144
PlutoTest.@test Float32 == eltype(add_gauss_noise_fl!(ones(Float32, (512, 512)), 0.3))

# ╔═╡ cc7bbb55-b8dd-454b-b170-45137ab0c0b5
md"###### Check wether mean and standard deviation are correct"

# ╔═╡ c997c248-fba8-4ed8-ab75-179d248e6b83
PlutoTest.@test ≈(0.3, std(add_gauss_noise(ones((512, 512)), 0.3)), rtol=0.1)

# ╔═╡ 07bc00b7-81ad-4645-8fbe-8545ccdfe7e4
PlutoTest.@test ≈(0.3, std(add_gauss_noise_fl!(ones((512, 512)), 0.3)), rtol=0.1)

# ╔═╡ d317ed1a-1867-4293-8db4-51352212929c
PlutoTest.@test ≈(1, mean(add_gauss_noise(ones((512, 512)), 0.3)), rtol=0.1)

# ╔═╡ e69a8b77-833a-4498-8bc9-e4c970f6529d
md"## Task 1.2 - Poisson Noise
In microscopy and low light imaging situation, the dominant noise term is usually Poisson noise which we want to simulate here.

For adding Poisson noise we use: [PoissonRandom.jl](https://github.com/SciML/PoissonRandom.jl)

Read the documentation of this package how to generate Poisson Random numbers
"

# ╔═╡ 1e4b32cb-00ea-439f-900e-abadc850be95
"""
	add_poisson_noise!(img, scale_to=nothing)

This function adds poisson distributed noise to `img`.
Before adding noise, it scales the maximum value to `scale_to` and
divides by it afterwards.
With that we can set the number of events (like a photon count)

If `isnothing(scale_to) == true`, we don't modify/scale the array.

`!` means that the input is modified.
"""
function add_poisson_noise!(img, scale_to=nothing)
	# todo
	return img
end

# ╔═╡ 6cf8ff2a-4255-4c53-9376-5a0d2d372569
add_poisson_noise!(10 .* ones(Float32, (3, 3)))

# ╔═╡ e3fcd381-9652-4e1d-8585-a91f87b73ce5
md"### Test 1.2 - Poisson Noise

You probably encounter errors for Float32 types errors, try to put a `Float64(some_part)` at the right place.

Works like this:
* `typeof(1f0)` = $(typeof(1f0))
* `typeof(Float64(1f0))` = $(typeof(Float64(1f0)))

"

# ╔═╡ f3efcfcf-4516-43d2-9375-d8863b4b6f5b
gray_show(add_poisson_noise!(100 .* img), set_one=true)

# ╔═╡ 00c72e52-09fa-413b-acfd-83c760d42d7f
PlutoTest.@test Float32 == eltype(add_poisson_noise!(100 .* ones(Float32, (512, 512)), 0.3f0))

# ╔═╡ 5b68030f-557e-449e-8baf-cf7ae65e4425
mean(add_poisson_noise!(150 .* ones(Float64, (512, 512))))

# ╔═╡ 78a39a77-bc9d-427a-a113-4e64b0c9d311
PlutoTest.@test ≈(150, mean(add_poisson_noise!(150 .* ones(Float64, (512, 512)))), rtol=0.05)

# ╔═╡ 41118b3c-3e49-4fc7-8c80-7906b7cf91fe
PlutoTest.@test ≈(√(150), std(add_poisson_noise!(150 * ones(Float32, (512, 512)))), rtol=0.05)

# ╔═╡ f8be3a4b-d8f0-40cb-869e-81c17e58327a
md"###### Consider those tests as bonus"

# ╔═╡ 65a8981d-890b-4590-a913-6890ecf3817d
PlutoTest.@test ≈(150, 150 * mean(add_poisson_noise!(ones(Float32, (512, 512)), 150)), rtol=0.05)

# ╔═╡ 272c8afb-667f-49c7-af48-2feb1b9d06f7
PlutoTest.@test ≈(√(150), 150 * std(add_poisson_noise!(ones(Float32, (512, 512)), 150)), rtol=0.05)

# ╔═╡ 59482e49-d784-4303-9487-bf37f6a0462e
md"## 1.3 Hot Pixels

Another issues are hot pixels which show maximum value. This can be due to damaged pixels or some other noise (radioactivity, ...). Often this is called Salt (because white) noise

"

# ╔═╡ 6e050fc2-4db2-4e6e-abd4-2812b47c070f
"""
	add_hot_pixels(img, probability=0.1; max_value=one(eltype(img)))

Add randomly hot pixels. The probability for each pixel to be hot,
should be specified by `probability`.
`max_value` is a keyword argument which is the value the _hot_ pixel will have.
"""
function add_hot_pixels!(img, probability=0.1; max_value=one(eltype(img)))
	# todo
	return img
end

# ╔═╡ 22985061-9740-45e9-af56-e0787dadba7b
gray_show(add_hot_pixels!(copy(img)))

# ╔═╡ 42f92a8c-54bb-4033-9ef2-f5ffb78f4f3a
md"### 1.3 Test"

# ╔═╡ a9ad9d52-53c9-44ce-8edc-af8d7bfdc81f
PlutoTest.@test sum(map(x -> x ≈ 2, add_hot_pixels!(copy(img), 0.5, max_value=2))) ≈ length(img) .* 0.5 rtol=0.05

# ╔═╡ 373a3e39-fd8c-4cea-88bb-2b9c3734f14f
md"## 2. Remove Noise From Images"

# ╔═╡ 2c274829-ef33-4f34-9ba0-63d3b8cd7343
md"## 2.1 Remove Noise with Gaussian Blur Kernel

One option to remove noise is, is to blur the image with a Gaussian kernel.
We can convolve a small gaussian (odd sized) Kernel over the array.
For that we are using `ImageFiltering.imfilter` for the convolution.

To create a Gaussian shaped function, checkout `IndexFunArrays.normal`.
You probably want to normalize the sum of it again
"

# ╔═╡ 52ebebfc-940b-47e2-a059-40ec996f2f98
begin
	img_g = add_gauss_noise(img, 0.1)
	img_p = add_poisson_noise!(100 .* copy(img)) ./ 100
	img_h = add_hot_pixels!(copy(img))
end;

# ╔═╡ 2d6f25a5-b902-4e46-ac3f-cf452f525912
gray_show([img_g img_p img_h])

# ╔═╡ 75aa8212-8ba3-4fba-b8b2-af2dd9c82adf
# check out the help
imfilter

# ╔═╡ 6acf7ba6-445d-4ac0-a83f-6916318c783c
function gaussian_noise_remove(arr; kernel_size=(3,3), σ=1)
	# todo
	return similar(arr)
end

# ╔═╡ 4f683117-9a09-4837-b0af-5818ac2e62fd
md"σ = $(@bind σ Slider(0.01:0.1:10, show_value=true))"

# ╔═╡ 01203c2c-4e67-492c-9d9c-dc84a84afe15
md"kernel size = $(@bind ks Slider(3:2:21, show_value=true))"

# ╔═╡ daae9a7d-99f4-4d34-a7fd-886fa754d59b
begin
	img_p_gauss = gaussian_noise_remove(img_p, kernel_size=(ks,ks), σ=σ);
	img_g_gauss = gaussian_noise_remove(img_g, kernel_size=(ks,ks), σ=σ);
	img_h_gauss = gaussian_noise_remove(img_h, kernel_size=(ks,ks), σ=σ);
end;

# ╔═╡ 92d43328-29a2-45b2-ae71-e169dec02316
gray_show([img_p_gauss img_g_gauss img_h_gauss])

# ╔═╡ 512486eb-1919-446a-9f24-e9713546d237
md"### 2.1 Test"

# ╔═╡ 9b83d92d-f838-4329-b0b1-62bf7a18f123
arr_rand = add_gauss_noise(ones((500, 500)), 0.2);

# ╔═╡ bbbed1b7-a406-44a7-8c77-959a3eba0197
PlutoTest.@test std(gaussian_noise_remove(arr_rand, kernel_size=(8,8), σ=2)) < 0.05

# ╔═╡ 911c2e1e-f04a-4efa-9aff-7f4a424011ae
PlutoTest.@test sum(abs2, img_g .- img) > sum(abs2, img_g_gauss .- img)

# ╔═╡ 0f6ff33f-e50a-4317-a240-3b2f52de7a97
md"## 2.2 Noise Removal with Median Filter
The median filter slides with an quadratic shaped box over the arrays and
always takes the median value of this array.

"

# ╔═╡ 7bc0f845-d5d4-4192-be1d-b97750e95761
function median_noise_remove!(arr; kernel_size=(5,5))
	# todo
	return similar(arr)
end

# ╔═╡ ac6e7ea7-c277-4180-a809-a37531ab4f11
md"kernel\_size\_2 = $(@bind ks_2 Slider(3:2:9, show_value=true))"

# ╔═╡ 8208b936-ad5d-45a0-96d1-04fd131f1e29
begin
	img_p_median = median_noise_remove!(copy(img_p), kernel_size=(ks_2,ks_2));
	img_g_median = median_noise_remove!(copy(img_g), kernel_size=(ks_2,ks_2));
	img_h_median = median_noise_remove!(copy(img_h), kernel_size=(ks_2,ks_2));
end;

# ╔═╡ 06af0b94-f4e1-4114-b565-730b87567995
gray_show([img_p_median img_g_median img_h_median])

# ╔═╡ 707a3be6-ac83-421d-a8e1-31d75444aeef
md"### 2.2 Test"

# ╔═╡ 07f5e6da-7868-4e22-90aa-7a81c2b41ab3
PlutoTest.@test sum(abs2, img_p .- img) > sum(abs2, img_p_median .- img)

# ╔═╡ 3dbe2b38-8609-41e5-b3b9-d62da394b89d
PlutoTest.@test sum(abs2, img_g .- img) > sum(abs2, img_g_median .- img)

# ╔═╡ 867d5676-8dfc-4bb5-b555-73aa767f7e9e
md"## 3 Final Images"

# ╔═╡ 3605fda7-078e-4bff-9e33-91c4025e775a
gray_show([img_g img_p img_h])

# ╔═╡ c09b3470-95bd-4f49-b64b-26459c95ab19
gray_show([median_noise_remove!(copy(img_g)) median_noise_remove!(copy(img_p)) median_noise_remove!(copy(img_h))])

# ╔═╡ Cell order:
# ╟─31cc54e9-d47b-4df7-a785-dad0120e8cdb
# ╠═56807afc-164e-11ec-2cf4-354654176242
# ╠═8a11c2b7-b6e5-4c26-9e06-ad3ab958db9c
# ╠═d8c83f2c-76ea-42b9-bd39-7fa5cafcd4d3
# ╠═334ebdc5-dc11-4e83-ad12-a1961420a97b
# ╟─2122c12f-0832-4f19-9a77-3047d5311cba
# ╟─a7d7a415-7464-43fb-a2bf-8e3f2d999c32
# ╠═de104b0a-28a4-42de-ae8d-fa857a8f32e0
# ╟─f89bc1e4-d406-4550-988b-71496b64035a
# ╠═c0b5ff08-7342-4c25-9791-be3a4918c400
# ╟─d56bc500-f5cf-40f3-8c0b-b4ae26a58367
# ╠═1116455e-4698-489a-add2-4e1d8ecbb546
# ╟─8e343a3d-ab0f-45a1-9161-21fec07755c9
# ╠═d865dddc-6d32-416b-ba12-8f6b3870e771
# ╠═5c9a8c7b-e9c4-4177-86bb-b69113fce144
# ╟─cc7bbb55-b8dd-454b-b170-45137ab0c0b5
# ╠═c997c248-fba8-4ed8-ab75-179d248e6b83
# ╠═07bc00b7-81ad-4645-8fbe-8545ccdfe7e4
# ╠═d317ed1a-1867-4293-8db4-51352212929c
# ╟─e69a8b77-833a-4498-8bc9-e4c970f6529d
# ╠═3ccbf4fb-0707-4ce8-83d8-b1af555f3fe1
# ╠═1e4b32cb-00ea-439f-900e-abadc850be95
# ╠═6cf8ff2a-4255-4c53-9376-5a0d2d372569
# ╟─e3fcd381-9652-4e1d-8585-a91f87b73ce5
# ╠═f3efcfcf-4516-43d2-9375-d8863b4b6f5b
# ╠═00c72e52-09fa-413b-acfd-83c760d42d7f
# ╠═5b68030f-557e-449e-8baf-cf7ae65e4425
# ╠═78a39a77-bc9d-427a-a113-4e64b0c9d311
# ╠═41118b3c-3e49-4fc7-8c80-7906b7cf91fe
# ╠═f8be3a4b-d8f0-40cb-869e-81c17e58327a
# ╠═65a8981d-890b-4590-a913-6890ecf3817d
# ╠═272c8afb-667f-49c7-af48-2feb1b9d06f7
# ╟─59482e49-d784-4303-9487-bf37f6a0462e
# ╠═6e050fc2-4db2-4e6e-abd4-2812b47c070f
# ╠═22985061-9740-45e9-af56-e0787dadba7b
# ╟─42f92a8c-54bb-4033-9ef2-f5ffb78f4f3a
# ╠═a9ad9d52-53c9-44ce-8edc-af8d7bfdc81f
# ╟─373a3e39-fd8c-4cea-88bb-2b9c3734f14f
# ╟─2c274829-ef33-4f34-9ba0-63d3b8cd7343
# ╠═52ebebfc-940b-47e2-a059-40ec996f2f98
# ╠═2d6f25a5-b902-4e46-ac3f-cf452f525912
# ╠═62c6e2f0-0007-40ec-81c5-61cda80f59e5
# ╠═75aa8212-8ba3-4fba-b8b2-af2dd9c82adf
# ╠═6acf7ba6-445d-4ac0-a83f-6916318c783c
# ╟─4f683117-9a09-4837-b0af-5818ac2e62fd
# ╟─01203c2c-4e67-492c-9d9c-dc84a84afe15
# ╠═daae9a7d-99f4-4d34-a7fd-886fa754d59b
# ╠═92d43328-29a2-45b2-ae71-e169dec02316
# ╟─512486eb-1919-446a-9f24-e9713546d237
# ╠═9b83d92d-f838-4329-b0b1-62bf7a18f123
# ╠═bbbed1b7-a406-44a7-8c77-959a3eba0197
# ╠═911c2e1e-f04a-4efa-9aff-7f4a424011ae
# ╟─0f6ff33f-e50a-4317-a240-3b2f52de7a97
# ╠═7bc0f845-d5d4-4192-be1d-b97750e95761
# ╟─ac6e7ea7-c277-4180-a809-a37531ab4f11
# ╠═8208b936-ad5d-45a0-96d1-04fd131f1e29
# ╠═06af0b94-f4e1-4114-b565-730b87567995
# ╟─707a3be6-ac83-421d-a8e1-31d75444aeef
# ╠═07f5e6da-7868-4e22-90aa-7a81c2b41ab3
# ╠═3dbe2b38-8609-41e5-b3b9-d62da394b89d
# ╟─867d5676-8dfc-4bb5-b555-73aa767f7e9e
# ╠═3605fda7-078e-4bff-9e33-91c4025e775a
# ╠═c09b3470-95bd-4f49-b64b-26459c95ab19


text:
@testset "6.1 Hyperbolic sine" begin
    include("6.1.1 (c+d x)^m (a+b sinh)^n.jl")
    include("6.1.3 (e x)^m (a+b sinh(c+d x^n))^p.jl")
    include("6.1.4 (d+e x)^m sinh(a+b x+c x^2)^n.jl")
    include("6.1.5 Hyperbolic sine functions.jl")
    include("6.1.7 hyper^m (a+b sinh^n)^p.jl")
end

max_stars_repo_path: test/6 Hyperbolic functions/6.1 Hyperbolic sine/.jl


text:
using Proj
using BenchmarkTools
using StaticArrays

const ctrans = Proj.Transformation("EPSG:4326", "EPSG:28992", always_xy = true)

function timer(p)
    tot = 0.0
    for _ = 1:10_000
        tot += sum(ctrans(p))
    end
    tot
end

function timer(x, y)
    tot = 0.0
    for _ = 1:10_000
        tot += sum(ctrans(x, y))
    end
    tot
end

# test different input types
svector = SA[5.39, 52.16]
vector = [5.39, 52.16]
tupl = (5.39, 52.16)
x, y = 5.39, 52.16

# timed using Julia 1.8.0-beta3
@btime Proj.Coord($svector)  # 2.900 ns (0 allocations: 0 bytes)
@btime Proj.Coord($vector)   # 3.700 ns (0 allocations: 0 bytes)
@btime Proj.Coord($tupl)     # 2.600 ns (0 allocations: 0 bytes)
@btime Proj.Coord($x, $y)    # 2.800 ns (0 allocations: 0 bytes)

@btime timer($svector)       # 9.899 ms (0 allocations: 0 bytes)
@btime timer($vector)        # 9.899 ms (0 allocations: 0 bytes)
@btime timer($tupl)          # 9.908 ms (0 allocations: 0 bytes)
@btime timer($x, $y)         # 9.977 ms (0 allocations: 0 bytes)


text:
# In this tutorial, we present the use of Stopping to checkpointing.

# When using an optimizer for high-scale problems, the resolution process
# might be extremly long. In order to analyze the progress of the algorithm or save ongoing results,
# an idea is to introduce checkpointing, i.e. we save the output result in the file
# every `n`-steps. Using `Stopping` this operation is now very simple.

using ADNLPModels, FileIO, JLD2, LinearAlgebra, NLPModels, Printf, Random, Stopping

# In this tutorial, we will use the steepest descent method with a fixed stepsize.

function fixed_step_steepest_descent(stp :: NLPStopping; t = 1e-5)

  xk = stp.current_state.x

  OK = update_and_start!(stp, gx = grad(stp.pb, xk))

  @printf "%2s %7s\n" "k" "||∇f(x)||"
  @printf "%2d %7.1e\n" stp.meta.nb_of_stop norm(stp.current_state.current_score)
  while !OK
    xk -= t * stp.current_state.gx

    OK = update_and_stop!(stp, x = xk, gx = grad(stp.pb, xk))

    @printf "%2d %7.1e\n" stp.meta.nb_of_stop norm(stp.current_state.current_score)
  end
  return stp
end

# We now generate a regularized least squares problem using `ADNLPModels.jl`.

Random.seed!(1234)
m, n = 10_000, 10
A  = rand(m, n)
b  = A * ones(n)
f(x, A, b, λ) = norm(A * x - b)^2 + λ * norm(x)^2
pb = ADNLPModel(x -> f(x, A, b, 1e-2), zeros(n))

# The final step is now to initialize the `Stopping` and specify user-defined structures
# to store the parameter `n_save` set to 50 so that every 50 iterations the current
# stopping is saved using the package `JLD2.jl`.

save_check(stp, b) = begin
  if stp.meta.nb_of_stop % stp.stopping_user_struct[:n_save] == 0
    @save "checkpoint_stopping_$(stp.meta.nb_of_stop).jld2" stp
  end
end
n_save = 50
stp = NLPStopping(pb, user_struct = Dict(:n_save => n_save),
                      user_check_func! = save_check, max_iter = 99)
# Let's go
fixed_step_steepest_descent(stp)

# The algorithm has now generated two files `checkpoint_stopping_0.jld2` and `checkpoint_stopping_50.jld2` that
# can be analyzed.

stp0 = load("checkpoint_stopping_0.jld2")["stp"]
stp50 = load("checkpoint_stopping_50.jld2")["stp"]


text:
## Convert orignal Coconet training data into Julia format
using BSON
using PyCall
using Coconet

#===============================================================================
    SECTION 1: Python to Julia
===============================================================================#

path = "datasets"
name = "Tonnetz_Bach.bson"

let np = pyimport("numpy")
    # Load data as PyObject
    data = np.load(
        joinpath(path, "Jsb16thSeparated.npz"),
        allow_pickle = true,
        encoding = "latin1",
    )

    # Convert to Julia vectors
    test_data = [[timestep for timestep in song] for song in get(data, :test)]
    train_data = [[timestep for timestep in song] for song in get(data, :train)]
    valid_data = [[timestep for timestep in song] for song in get(data, :valid)]

    # Convert to arrays
    test_data = [hcat(song...) for song in test_data]
    train_data = [hcat(song...) for song in train_data]
    valid_data = [hcat(song...) for song in valid_data]

    # Transpose and re-convert to arrays
    test_data = [hcat(song') for song in test_data]
    train_data = [hcat(song') for song in train_data]
    valid_data = [hcat(song') for song in valid_data]

    # Save as dict
    bson(
        joinpath(path, name),
        test = test_data,
        train = train_data,
        valid = valid_data,
    )
end
#===============================================================================
    SECTION 2: Notes to tonnetz vectors
===============================================================================#

# Load in as dict
bach = BSON.load(joinpath(path, name))

notes = []
dataset = []
# for each component of the dataset
for test in [bach[v] for v in [:test, :train, :valid]]
    set = []

    # for each song
    for i = 1:length(test)
        song = zeros(Float32, 3, size(test[i])...)

        # for each instrument
        for instrument = 1:size(test[i])[2]
            # for each timestep
            for timestep = 1:size(test[i])[1]
                pitch = test[i][timestep, instrument]
                x, y, z = 0, 0, 0

                if !isnan(pitch)
                    z = 1
                    x, y = Coconet.pitch_to_tonnetz(Integer(pitch))
                end

                song[:, timestep, instrument] = [x,y,z]
                push!(notes, song[:, timestep, instrument])
            end # for each timestep
        end # for each instrument

        push!(set, song)
    end #for each song

    push!(dataset, set)
end # for each set

# Normalization
μ = sum(notes) / length(notes)
σ = sum([(note - μ).^2 for note in notes]) / length(notes)
for i in 1:length(dataset)
    for j in 1:length(dataset[i])
        for k in 1:size(dataset[i][j], 2)
            for l in 1:size(dataset[i][j], 3)
                dataset[i][j][:, k, l] = (dataset[i][j][:, k, l] - μ) ./ σ
            end
        end
    end
end

#===============================================================================
    SECTION 3: Save it
===============================================================================#
bson(
    joinpath(path, name),
    test = dataset[1],
    train = dataset[2],
    valid = dataset[3],
    μ = μ,
    σ = σ,
)

max_stars_repo_name: MacKenzieHnC/Coconet

max_issues_repo_name: MacKenzieHnC/Coconet

max_forks_repo_name: MacKenzieHnC/Coconet


text:
### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("FVheat2d");

useLog("FVheat2dlog")

# Configuration setup
domain(2)
solverType(FV)
timeStepper(RK4)

# Mesh
n = 20 # number of elements
mesh(QUADMESH, elsperdim=n)

# Variables and BCs
u = variable("u", location=CELL)
boundary(u, 1, FLUX, "0")

# Time interval and initial condition
T = 0.1;
timeInterval(T)
initial(u, "(sin(pi*x)*sin(2*pi*y))^4")

# The flux and source terms of the conservation equation
# F and S in the following equation:
# Dt(int(u dx)/A) = int(S dx) - int(F.n ds)
coefficient("D", 0.1) # Diffusion rate
flux(u, "-D*dot(grad(u),normal())")

printLatex(u)

#exportCode("fvheat2dcode") # uncomment to export generated code to a file
#importCode("fvheat2dcode") # uncomment to import code from a file

solve(u)

# output result to file
output_values(u, "fvheat2d", format="vtk");

@finalize()

##### Uncomment below to plot
# xy = Finch.fv_info.cellCenters

# using Plots
# pyplot();
# display(plot(xy[1,:], xy[2,:], u.values[:], st=:surface))

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
######### Stan program example  ###########

using Mamba, Stan

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  bernoullimodel = "
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
    y ~ bernoulli(theta);
  }
  "

  observeddata = [
    Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1]),
    Dict("N" => 10, "y" => [0, 1, 0, 0, 0, 0, 1, 0, 0, 1]),
    Dict("N" => 10, "y" => [0, 0, 0, 0, 0, 0, 1, 0, 1, 1]),
    Dict("N" => 10, "y" => [0, 0, 0, 1, 0, 0, 0, 1, 0, 1])
  ]

  initparms = Dict{String, Any}[
    Dict("theta" => 0.1),
    Dict("theta" => 0.4),
    Dict("theta" => 0.5),
    Dict("theta" => 0.9),
  ]

  monitor = ["theta", "lp__", "accept_stat__"]

  global stanmodel, rc, sim1
  stanmodel = Stanmodel(name="bernoulli", model=bernoullimodel, num_warmup=1);
  rc, sim1 = stan(stanmodel, observeddata, init=initparms, CmdStanDir=CMDSTAN_HOME)

  rc == 0 && sim1.value[1, end, :] |> display
end # cd


text:
export NDMComplex

struct NDMComplex{VT,MT} <: MatrixNeuralNetwork
    a::VT
    b::VT
    c::VT
    W::MT
    X::MT
end
@treelike NDMComplex

"""
    NDMComplex([T=STD_REAL_PREC], N, αₕ, αₐ, [initW, initb, inita])

Constructs a Neural Density Matrix with numerical precision `T` (Defaults to
Float32), `N` input neurons, N⋅αₕ hidden neurons and N⋅αₐ ancillary neurons.
This network ensure that the density matrix is always positive definite.

The number of input neurons `N` must match the size of the lattice.

This network is very similar to `NDM`, but instead of splitting the modulus and
phase of the density matrix in two disjoint networks, it uses a single network
with complex weights.

The initial parameters of the neurons are initialized with a rescaled normal
distribution of width 0.01 for the coupling matrix and 0.005 for the local
biases. The default initializers can be overriden by specifying

initW=(dims...)->rescaled_normal(T, 0.01, dims...),
initb=(dims...)->rescaled_normal(T, 0.005, dims...),
inita=(dims...)->rescaled_normal(T, 0.005, dims...))

Refs:
"""
NDMComplex(args...) = NDMComplex(Complex{STD_REAL_PREC}, args...)
NDMComplex(::Real, ::Int) = throw("NDMComplex needs complex type")
NDMComplex(T::Type{<:Complex}, in, αh, αa,
    initW=(dims...)->rescaled_normal(T, 0.01, dims...),
    initb=(dims...)->rescaled_normal(T, 0.005, dims...),
    inita=(dims...)->rescaled_normal(T, 0.005, dims...)) =
    NDMComplex(inita(in),
               initb(convert(Int,αh*in)),
               initb(convert(Int,αa*in)),
               initW(convert(Int,αa*in), in),
               initW(convert(Int,αh*in), in))

input_type(net::NDMComplex{VT,MT}) where {VT,MT} = real(eltype(VT))
weight_type(net::NDMComplex) = eltype(net.a)
out_type(net::NDMComplex{VT,MT}) where {VT,MT} = eltype(VT)
input_shape(net::NDMComplex) = (length(net.a), length(net.a))
random_input_state(net::NDMComplex{VT,MT}) where {VT,MT} =
    (input_type(net).([rand(0:1) for i=1:length(net.a)]), input_type(net).([rand(0:1) for i=1:length(net.a)]))
is_analytic(net::NDMComplex) = false

Base.show(io::IO, m::NDMComplex) = print(io,
    "NDMComplex($(eltype(m.a)), n=$(length(m.a)), αₕ=$(length(m.b)/length(m.a)), αₐ=$(length(m.c)/length(m.a)))")

Base.show(io::IO, ::MIME"text/plain", m::NDMComplex) = print(io,
    "NDMComplex($(eltype(m.a)), n=$(length(m.a)), α=$(length(m.b)/length(m.a)), αₐ=$(length(m.c)/length(m.a)))")

@inline (net::NDMComplex)(σ::State) = net(config(σ)...)
@inline (net::NDMComplex)(σ::Tuple) = net(σ...)
function (W::NDMComplex)(σ, σp)
    T=eltype(W.a)

    θ̃     = W.c +
              W.W * σ +
                conj(W.W) * σp

    θ_σ   = W.b +
              W.X*σ

    θ_σp   = W.b +
              W.X*σp

    logψ = dot(σ, W.a) + dot(W.a, σp) + #dot(σp, conj(W.a)) +
               sum(logℒ2.(θ̃))    +
                   sum(logℒ2.(θ_σ))  +
                        sum(logℒ2.(conj.(θ_σp))) + log(T(1/8))
    return logψ
end

# Cached version
mutable struct NDMComplexCache{VT,VTR} <: NNCache{NDMComplex}
    θ̃_r::VT
    θ̃::VT
    θ_σp::VT
    θ_σ::VT

    logℒθ̃::VT
    logℒθ_σp::VT
    logℒθ_σ::VT

    ∂logℒθ̃::VT
    ∂logℒθ_σp::VT
    ∂logℒθ_σ::VT

    σ::VT
    σp::VT
    ∑σ::VT
    Δσ::VT

    σ_row_cache::VTR
    i_σ_row_cache::Int

    valid::Bool # = false
end

cache(net::NDMComplex) =
    NDMComplexCache(similar(net.c),
              similar(net.c),
              similar(net.b),
              similar(net.b),

              similar(net.c),
              similar(net.b),
              similar(net.b),

              similar(net.c),
              similar(net.b),
              similar(net.b),

              similar(net.a),
              similar(net.a),
              similar(net.a),
              similar(net.a),

              similar(real.(net.a)),
              -1,

              false)

(net::NDMComplex)(c::NDMComplexCache, σ) = net(c, config(σ)...)
function (W::NDMComplex)(c::NDMComplexCache, σ_r, σp_r)
    T        = eltype(W.W)
    θ̃        = c.θ̃
    θ̃_r      = c.θ̃_r
    θ_σ      = c.θ_σ
    θ_σp     = c.θ_σp
    logℒθ̃    = c.logℒθ̃
    logℒθ_σ  = c.logℒθ_σ
    logℒθ_σp = c.logℒθ_σp

    # copy the states to complex valued states for the computations.
    σ  = c.σ;  copyto!(σ,  σ_r)
    σp = c.σp; copyto!(σp, σp_r)

    if !c.valid || c.σ_row_cache ≠ σ
        c.σ_row_cache .= σ
        c.valid = true

        #θ̃_r = W.c + W.W*σ
        mul!(θ̃_r, W.W, σ)
        θ̃_r .+= W.c

        #θ_σp   = W.b +  W.X*σ
        mul!(θ_σ, W.X, σ)
        θ_σ .+= W.b

        logℒθ_σ .= logℒ2.(θ_σ)
    end

    #θ̃ .+= conj(W.W) * σp
    mul!(θ̃, W.W, σp)
    conj!(θ̃)
    θ̃ .+= θ̃_r

    #θ_σ   = W.b + W.X*σp
    mul!(θ_σp, W.X, σp)
    θ_σp .+= W.b

    logℒθ̃    .= logℒ2.(θ̃)
    logℒθ_σp .= logℒ2.(conj.(θ_σp))

    dotσ  = dot(σ, W.a)
    dotσp = dot(W.a, σp)
    logψ  = log(T(1/8)) +  sum(logℒθ_σp) + sum(logℒθ_σ) + sum(logℒθ̃) +
              dotσ + dotσp

    return logψ
end

function logψ_and_∇logψ!(∇logψ, W::NDMComplex, c::NDMComplexCache, σ_r, σp_r)
    T        = eltype(W.W)
    θ̃        = c.θ̃
    θ̃_r      = c.θ̃_r
    θ_σ      = c.θ_σ
    θ_σp     = c.θ_σp
    logℒθ̃    = c.logℒθ̃
    logℒθ_σ  = c.logℒθ_σ
    logℒθ_σp = c.logℒθ_σp

    # copy the states to complex valued states for the computations.
    σ  = c.σ;  copyto!(σ,  σ_r)
    σp = c.σp; copyto!(σp, σp_r)

    if !c.valid || c.σ_row_cache ≠ σ
        c.σ_row_cache .= σ
        c.valid = true

        #θ̃_r = W.c + W.W*σ
        mul!(θ̃_r, W.W, σ)
        θ̃_r .+= W.c

        #θ_σp   = W.b +  W.X*σ
        mul!(θ_σ, W.X, σ)
        θ_σ .+= W.b

        logℒθ_σ .= logℒ2.(θ_σ)
    end

    #θ̃ .+= conj(W.W) * σp
    mul!(θ̃, W.W, σp)
    conj!(θ̃)
    θ̃ .+= θ̃_r

    #θ_σ   = W.b + W.X*σp
    mul!(θ_σp, W.X, σp)
    θ_σp .+= W.b

    logℒθ̃    .= logℒ2.(θ̃)
    logℒθ_σp .= logℒ2.(conj.(θ_σp))

    dotσ  = dot(σ, W.a)
    dotσp = dot(W.a, σp)
    logψ  = log(T(1/8)) +  sum(logℒθ_σp) + sum(logℒθ_σ) + sum(logℒθ̃) +
              dotσ + dotσp
    # --- End common terms with computation of ψ --- #
    ∑σ    = c.∑σ
    Δσ    = c.Δσ
    c.∑σ .= σ .+ σp
    c.Δσ .= σ .- σp

    # Compute additional terms for derivatives
    ∂logℒθ̃    = c.∂logℒθ̃;    ∂logℒθ̃    .= ∂logℒ2.(θ̃)
    ∂logℒθ_σ  = c.∂logℒθ_σ;  ∂logℒθ_σ  .= ∂logℒ2.(θ_σ)
    ∂logℒθ_σp = c.∂logℒθ_σp; ∂logℒθ_σp .= ∂logℒ2.(conj.(θ_σp))

    # Store the derivatives
    ∇logψ_r = real(∇logψ)
    ∇logψ_i = imag(∇logψ)

    ∇logψ_r.a .=     ∑σ
    ∇logψ_i.a .= im.*Δσ

    ∇logψ_r.b .=      ∂logℒθ_σ .+ ∂logℒθ_σp
    ∇logψ_i.b .= im.*(∂logℒθ_σ .- ∂logℒθ_σp)

    ∇logψ_r.X .=      ∂logℒθ_σ.*transpose(σ) .+ ∂logℒθ_σp.*transpose(σp)
    ∇logψ_i.X .= im.*(∂logℒθ_σ.*transpose(σ) .- ∂logℒθ_σp.*transpose(σp))

    ∇logψ_r.c .=      ∂logℒθ̃
    ∇logψ_i.c .=      ∂logℒθ̃

    ∇logψ_r.W .=      ∂logℒθ̃.*transpose(∑σ)
    ∇logψ_i.W .=  im.*∂logℒθ̃.*transpose(Δσ)

    return logψ
end


text:
using LinearAlgebra

"""
    KalmanObservation{T}(obs, obsCov)

The Kalman filter observation and covariance matrices

# Arguments
-`obs`: Observation `Matrix`
-`obsCov`: Observation `Matrix` covariance

See also [`kfsystem`](@ref), [`kfstate`](@ref)
"""
mutable struct KalmanObservation{T <: Real}

    # Measurement vector
    meas::Vector{T}
    # Observation matrix
    obs::Matrix{T}
    # Uncertainty on observation
    obsCov::Matrix{T}

end

KalmanObservation{T}(
        obs::Matrix{T},
        obsCov::Matrix{T}
) where {T<:Real} = KalmanObservation{T}(
    zeros(T, size(obs)[1]),
    obs,
    obsCov
);

KalmanObservation{T}(
        meas::Vector{T},
        obsCov::Matrix{T}
) where {T<:Real} = KalmanObservation{T}(
    meas,
    zeros(eltype(obsCov), size(obsCov)),
    obsCov
);

"""
    KalmanUpdate{T}(transition, control, noise)

The Kalman filter transition, control, and noise covariance matrices, collectively the update matrices

# Arguments
 -`transition`: State transition `Matrix`
 -`control`: Input/external control `Matrix `
 -`noise`: Estimated noise covariance `Matrix` of the state transition

See also [`kfobservation`](@ref), [`kfstate`](@ref)
"""
struct KalmanUpdate{T <: Real}
    # State transition matrix
    transition::Matrix{T}
    # Control matrix
    control::Matrix{T}
    # Input noise
    noise::Matrix{T}
end

"""
    KalmanState{T}(state, cov)

The Kalman filter state vector and associated covariance matrix

# Arguments
 -`state`: State `Vector`
 -`cov`: Covariance `Matrix` of the state

See also [`kfobservation`](@ref), [`kfsystem`](@ref)
"""
struct KalmanState{T <: Real}
    # State vector
    state::Vector{T}
    # Covariance matrix
    cov::Matrix{T}

end

KalmanState{T}(
        ks::KalmanState
) where {T<:Real} = begin
    KalmanState(ks.state, ks.cov)
end

"""
predict_state(state, system, input)

Predict the next state vector and covariance

# Arguments
 -`state`: [`kfstate`](@ref) instance
 -`update_matrices`: [`kfupdate`](@ref) instance
 -`input`: External input `Vector`

# Returns
[`kfstate`](@ref)

See also [`update_state`](@ref)
"""
function predict_state(
        state::KalmanState,
        update_matrices::KalmanUpdate,
        external_input::Vector
)::KalmanState
    # Update the state vector
    pred_state_vec::Vector = (update_matrices.transition * state.state
                              + update_matrices.control * external_input)
    pred_state_cov::Matrix = ((update_matrices.transition * state.cov
                               * transpose(update_matrices.transition))
                               + update_matrices.noise)

    pred_state = KalmanState{eltype(pred_state_vec)}(pred_state_vec, pred_state_cov)
    return pred_state
end

"""
    correct_state(predicted_state, current_filter, measurement)

Update the current state estimate with a measurement

# Arguments
 -`predicted_state`: [`kfstate`](@ref) instance after running [`extrapulate_state`](@ref)
 -`observation`: [`kfobservation`](@ref) instance dictating the observation of the system
 -`measurement`: The measurement `Vector`

# Returns
[`kfstate`](@ref)

See also [`extrapulate_state`](@ref)
"""
function correct_state(
        predicted_state::KalmanState,
        observation::KalmanObservation
)::KalmanState
    # Intermediate calculation (1)
    obs_transpose = transpose(observation.obs)

    # Compute the Kalman gain
    gain = (predicted_state.cov * obs_transpose
            * (observation.obs * predicted_state.cov * obs_transpose
               + observation.obsCov)^(-1))

    # Intermediate calculation (2)
    transformation = (I(size(predicted_state.cov)[1]) - gain * observation.obs)

    # Update the state estimate and covariance
    corr_state_vec = (predicted_state.state + gain
                      * (observation.meas - observation.obs * predicted_state.state))
    corr_cov = (transformation * predicted_state.cov * transpose(transformation)
                + gain * observation.obsCov * transpose(gain))
    corr_state = KalmanState{eltype(corr_state_vec)}(corr_state_vec, corr_cov)

    return corr_state
end;

max_stars_repo_name: hoganman/kalmanfilterjl

max_issues_repo_name: hoganman/kalmanfilterjl

max_forks_repo_name: hoganman/kalmanfilterjl


text:
#const CIS_LUT = SVector{64}(cis.((0:63) / 64 * 2π))
#function cis_fast(x)
#    @inbounds CIS_LUT[(floor(Int, x / 2π * 64) & 63) + 1]
#end

"""
$(SIGNATURES)

Get sine approximation very fast. It is valid between -π and π.
```julia-repl
julia> sin_vfast(π / 4)
```
"""
@inline function sin_vfast(x::T) where T <: Real
    B = T(4 / π)
    C = T(-4 / (π * π))
    B * x + C * x * abs(x)
end

"""
$(SIGNATURES)

Get cosine approximation very fast. It is valid between -π and π.
```julia-repl
julia> cos_vfast(π / 4)
```
"""
@inline function cos_vfast(x::T) where T <: Real
    x += T(π / 2)
    x -= (x > T(π)) * T(2π)
    sin_vfast(x)
end

"""
$(SIGNATURES)

Get sine approximation fast. It is valid between -π and π.
```julia-repl
julia> sin_fast(π / 4)
```
"""
@inline function sin_fast(x::T) where T <: Real
    P = T(0.225)
    y = sin_vfast(x)
    P * (y * abs(y) - y) + y
end

"""
$(SIGNATURES)

Get cosine approximation fast. It is valid between -π and π.
```julia-repl
julia> cos_fast(π / 4)
```
"""
@inline function cos_fast(x::T) where T <: Real
    x += T(π / 2)
    x -= (x > T(π)) * T(2π)
    sin_fast(x)
end

"""
$(SIGNATURES)

Get \exp(iz) approximation fast. It is valid between -π and π.
```julia-repl
julia> cos_fast(π / 4)
```
"""
function cis_fast(x)
    complex(cos_fast(x), sin_fast(x))
end

"""
$(SIGNATURES)

Get \exp(iz) approximation very fast. It is valid between -π and π.
```julia-repl
julia> cis_vfast(π / 4)
```
"""
function cis_vfast(x)
    complex(cos_vfast(x), sin_vfast(x))
end

"""
$(SIGNATURES)

Get \exp(iz) approximation fast. It is valid between -π and π.
```julia-repl
julia> get_carrier_fast_unsafe(π / 4)
```
"""
@inline function get_carrier_fast_unsafe(x)
    cis_fast(x)
end

"""
$(SIGNATURES)

Get \exp(iz) approximation very fast. It is valid between -π and π.
```julia-repl
julia> get_carrier_vfast_unsafe(π / 4)
```
"""
@inline function get_carrier_vfast_unsafe(x)
    cis_vfast(x)
end


text:
# a very simple MPI tag manager

"""
  Type that manages assigning MPI tags
"""
struct MPITagManager
  used_tags::Vector{Cint}
  start_tag::Int
end

"""
  Constructor

  **Inputs**

   * start_tag: tag to start counting from, default 5000
"""
function MPITagManager(start_tag=5000)
  used_tags = Array{Cint}(0)

  return MPITagManager(used_tags, start_tag)
end

"""
  Returns the next free tag
"""
function getNextTag(mgr::MPITagManager)

  if length(mgr.used_tags) == 0
    lasttag = mgr.start_tag
  else
    lasttag = mgr.used_tags[end]
  end

  newtag = lasttag + 1
  push!(mgr.used_tags, newtag)

  return newtag
end

#------------------------------------------------------------------------------
# default manager

global const TagManager = MPITagManager()


text:
# --------------------- #
# Post-estimation tools #
# --------------------- #

# Post-estimation methods           #
#-----------------------------------#

StatsBase.coef(e::MomentBasedEstimator{T}) where T = e.r.coef
StatsBase.coef(e::GenericMomentBasedEstimator) = e.coef

status(e::MomentBasedEstimator{T}) where T = e.r.status
objvalstatus(e::MomentBasedEstimator{T}) where T = e.r.objval

momentfunction(e::MomentBasedEstimator{T}) where T = momentfunction(e, Val{:smoothed})
momentfunction(e::MomentBasedEstimator{T}, ::Type{Val{:smoothed}}) where T = e.e.mf.s(coef(e))

momentfunction(e::MomentBasedEstimator{T}, ::Type{Val{:unsmoothed}}) where T = e.e.mf.g(coef(e))

momentfunction(e::MomentBasedEstimator{T}, theta) where T = e.e.mf.s(theta)
momentfunction(e::GenericMomentBasedEstimator, theta) where T = e.mf.s(theta)

function shrinkweight(p::Array{T}) where T
    mp = minimum(p)
    ϵ  = -min(mp, 0)
    (p .+ ϵ)./(1 .+ ϵ)    # n/(1 + e) + ne/(1+e) = n(1+e)/(1+e) = n
end

mdweights(e::MomentBasedEstimator{T}) where T<:MDEstimator= e.m.inner.x[1:nobs(e)]

function mdweights(e::MomentBasedEstimator{T}, ::Type{Val{:shrunk}}) where T<:MDEstimator
    shrinkweight(mdweights(e))::Array{Float64, 1}
end

function mdweights(e::MomentBasedEstimator{T}, ::Type{Val{:shrunk}}) where T<:GMMEstimator
    ones(first(size(e)))
end

# jacobian of moment function       #
#-----------------------------------#

function jacobian(e::MomentBasedEstimator{T}; weighted = true, shrinkweights = true) where T<:GenericMomentBasedEstimator
    t = weighted ? Val{:weighted} : Val{:unweighted}
    w = shrinkweights ? Val{:shrunk} : Val{:unshrunk}
    jacobian(e, t, w)
end

function jacobian(e::MomentBasedEstimator, t::Type{Val{:weighted}}, w)
    if isa(e.e.mf, MomentBasedEstimators.FADMomFun)
      p = mdweights(e, w)
      ws(theta) = momentfunction(e, theta)'*p
      theta = coef(e)
      ForwardDiff.jacobian(ws, theta)::Matrix{Float64}
    else
      e.e.mf.Dws(coef(e), mdweights(e, w))
    end
end

function jacobian(e::MomentBasedEstimator, t::Type{Val{:unweighted}}, w)
    if isa(e.e.mf, MomentBasedEstimators.FADMomFun)
      p = ones(first(size(e)))
      ws(theta) = At_mul_B(momentfunction(e, theta), p)
      theta = coef(e)
      ForwardDiff.jacobian(ws, theta)::Matrix{Float64}
    else
      e.e.mf.Dsn(coef(e))
    end
end

# smoothing of the moment function  #
#-----------------------------------#

κ₂(e::MomentBasedEstimator{T}) where T<:GenericMomentBasedEstimator = e.e.mf.kern.κ₂
κ₁(e::MomentBasedEstimator{T}) where T<:GenericMomentBasedEstimator = e.e.mf.kern.κ₁
κ₃(e::MomentBasedEstimator{T}) where T<:GenericMomentBasedEstimator = e.e.mf.kern.κ₃
bw(e::MomentBasedEstimator{T}) where T<:GenericMomentBasedEstimator = e.e.mf.kern.S

smoothing_kernel(e::MomentBasedEstimator{T}) where T<:GMMEstimator = e.e.mgr.k
smoothing_kernel(e::MomentBasedEstimator{T}) where T<:MDEstimator = e.e.mf.kern
iteration_manager(e::MomentBasedEstimator{T}) where T<:GMMEstimator = e.e.mgr

# covariance of the moment function #
#-----------------------------------#

function mfvcov(e::MomentBasedEstimator{T}, weighted::Bool = true, shrinkweights::Bool = true) where T<:MDEstimator
    t = weighted ? Val{:weighted} : Val{:unweighted}
    w = shrinkweights ? Val{:shrunk} : Val{:unshrunk}
    mfvcov(e, t, w)
end

function mfvcov(e::MomentBasedEstimator{T}, t::Type{Val{:weighted}}, w) where T<:MDEstimator
    mf = copy(momentfunction(e))
    p  = sqrt.(mdweights(e, w))
    broadcast!(*, mf, mf, p)
    S  = bw(e)
    k1 = κ₁(e)
    k2 = κ₂(e)
    Omega = vcov(mf, HC0())
    rmul!(Omega, k1^2*S/k2)
    return Omega
end

function mfvcov(e::MomentBasedEstimator{T}, t::Type{Val{:unweighted}}, w) where T<:MDEstimator
    mf = momentfunction(e)
    S  = bw(e)
    k1 = κ₁(e)
    k2 = κ₂(e)
    Omega = vcov(mf, HC0())
    Omega = scale!(k1^2*S/k2, Omega)
    return Omega
end

mfvcov(e::MomentBasedEstimator{T}) where T = vcov(momentfunction(e), smoothing_kernel(e))
mfvcov(e::MomentBasedEstimator{T}, k::RobustVariance) where T = vcov(momentfunction(e), k)

initial_weighting(e::MomentBasedEstimator{T}) where T<:GMMEstimator = e.e.W[end]

# covariance of the estimator       #
#-----------------------------------#

StatsBase.vcov(e::MomentBasedEstimator{T}) where T<:GMMEstimator = vcov(e, smoothing_kernel(e), iteration_manager(e))

function StatsBase.vcov(e::MomentBasedEstimator{T}, k::RobustVariance, mgr::TwoStepGMM) where T<:GMMEstimator
    n, p, m = size(e)
    G = jacobian(e)
    S = mfvcov(e, k)
    (n.^2/(n-p))*pinv(G'*pinv(S)*G)
end

function StatsBase.vcov(e::MomentBasedEstimator{T}, k::RobustVariance, mgr::OneStepGMM) where T<:GMMEstimator
    n, p, m = size(e)
    G = jacobian(e)
    S = mfvcov(e, k)
    W = initial_weighting(e)
    ## Use the general form of the variance covariance matrix
    ## that gives the correct covariance even when S \neq Var(\sqrt{N}g_N(\theta_0))
    A = pinv(G'*W*G)
    B = G'*W*S*W*G
    (n.^2/(n-p))*A*B*A
end

function StatsBase.vcov(e::MomentBasedEstimator{T}, k::RobustVariance) where T<:GMMEstimator
    vcov(e, k, iteration_manager(e))
end

function StatsBase.vcov(e::MomentBasedEstimator{T}, mgr::IterationManager) where T<:GMMEstimator
    vcov(e, smoothing_kernel(e), mgr)
end

function StatsBase.vcov(e::MomentBasedEstimator{T}; robust::Bool = false, weighted::Bool = true, shrinkweights = true) where T<:MDEstimator
    r = robust ? Val{:robust} : Val{:unrobust}
    t = weighted ? Val{:weighted} : Val{:unweighted}
    w = shrinkweights ? Val{:shrunk} : Val{:unshrunk}
    vcov(e, r, t, w)
end

function StatsBase.vcov(e::MomentBasedEstimator{T}, ::Type{Val{:robust}}, t, w) where T<:MDEstimator
    n, p, m = size(e)
    G = jacobian(e, t, w)
    S = mfvcov(e, t, w)
    V = G'pinv(S)*G/n
    H = inv(objhessian(e))
    V = H'*V*H
    sc = n/(n-p)
    return rmul!(V, sc)
end

function StatsBase.vcov(e::MomentBasedEstimator{T}, ::Type{Val{:unrobust}}, t, w) where T<:MDEstimator
    n, p, m = size(e)
    Ω = mfvcov(e, t, w)
    G = jacobian(e, t, w)
    V = pinv(G'pinv(Ω)*G/n)
    sc = n/(n-p)
    return rmul!(V, sc)
end

function StatsBase.vcov(e::MomentBasedEstimator{T}, k::RobustVariance) where T<:MDEstimator
    n, p, m = size(e)
    mf = momentfunction(e, Val{:unsmoothed})
    Ω = CovarianceMatrices.vcov(mf, k)
    G = jacobian(e, weighted = false, shrinkweights = false)
    V = pinv(G'pinv(Ω)*G/n)
    sc = n/(n-p)
    return rmul!(V, sc)
end

# hessian of the MD objective       #
#-----------------------------------#

function objhessian(m::MomentBasedEstimator{T}) where T<:MDEstimator
    @assert status(m) == :Optimal "the status of `::MDEstimator` is not :Optimal"
    n, _, p = size(m)
    ## In case KNITRO has ms_enable on disactivate it
    solver = m.s
    fixsolver!(solver)

    mdp = MinimumDivergenceProblem(Array{Float64}(undef, n,p), zeros(p), wlb = m.e.wlb,
                                   wub = m.e.wub, solver = solver, div = m.e.div)
    function obj(theta)
        mdp.e.mm.S[:] = m.e.mf.s(theta)
        solve!(mdp)
        mdp.m.inner.obj_val[1]
    end
    H = Calculus.hessian(obj, coef(m))
    k1 = κ₁(m)
    k2 = κ₂(m)
    S  = bw(m)
    k2*H/(S*k1^2)
end

# standard error of the estimator   #
#-----------------------------------#
#
# function StatsBase.stderror(e::MomentBasedEstimator{T}) where T<:MDEstimator
#     sqrt.(diag(vcov(e)))
# end
#
function StatsBase.stderror(e::MomentBasedEstimator{T}, mgr::IterationManager) where T<:GMMEstimator
    sqrt.(diag(vcov(e, smoothing_kernel(e), mgr)))
end
#
function StatsBase.stderror(e::MomentBasedEstimator{T}, k::RobustVariance) where T<:GMMEstimator
    sqrt.(diag(vcov(e, k)))
end

function StatsBase.stderror(e::MomentBasedEstimator{T}) where T<:GMMEstimator
    sqrt.(diag(vcov(e, e.e.mgr.k)))
end

#
function StatsBase.stderror(e::MomentBasedEstimator{T}; robust::Bool = false, weighted::Bool = true, shrinkweights::Bool = true) where T<:MDEstimator
    sqrt.(diag(vcov(e, robust = robust, weighted = weighted, shrinkweights = shrinkweights)))
end

# J test                            #
#-----------------------------------#

function J_test(e::MomentBasedEstimator)
    g = mean(momentfunction(e), dims = 1)
    S = pinv(mfvcov(e))
    j = (nobs(e)*(g*S*g'))[1]
    p = df(e) > 0 ? ccdf(Chisq(df(e)), j) : NaN
    # sometimes p is garbage, so we clamp it to be within reason
    return j, clamp(p, eps(), Inf)
end

df(e::MomentBasedEstimator{T}) where T = nmom(e) - npar(e)
z_stats(e::MomentBasedEstimator{T}) where T = coef(e) ./ stderror(e)
p_values(e::MomentBasedEstimator{T}) where T = 2*ccdf(Normal(), z_stats(e))

# coef table                        #
#-----------------------------------#

function StatsBase.coeftable(e::MomentBasedEstimator)
    cc = coef(e)
    se = stderror(e)
    zz = z_stats(e)
    CoefTable(hcat(cc, se, zz, 2.0*ccdf(Normal(), abs.(zz))),
              ["Estimate", "Std.Error", "z value", "Pr(>|z|)"],
              ["x$i" for i = 1:npar(e)],
              4)
end

# hack: fix solver                  #
#-----------------------------------#

fixsolver!(s::Ipopt.IpoptSolver) = s

# function fixsolver!(s::KnitroSolver)
#     opt = Array(Any, length(s.options))
#     for j in enumerate(s.options)
#         r = j[1]; o = j[2]
#         idx = (:outlev in o[1]) | (:ms_enable in o[1])
#         opt[r] = (o[1], idx ? o[2] : 0)
#     end
#     soptions[:] = opt
# end

# function Gadfly.plot{T <: MDEstimator}(m::MomentBasedEstimator{T}; xc = nothing, yc = nothing,  na...)
#     n, p, k = size(m)
#     @assert p <= 2 "no plotting for high dimensional MDEstimators"
#     solver = m.s
#     fixsolver!(solver)
#     mdp = MinimumDivergenceProblem(Array(Float64, (n, k)), zeros(k), wlb = m.e.wlb,
#                                    wub = m.e.wub, solver = solver, div = m.e.div)

#     cf = coef(m)
#     se = stderror(m)

#     xc = xc == nothing ? collect(cf[1]-5*se[1]:10*se[1]/20:cf[1]+5*se[1]) : xc
#     yc = yc == nothing && p > 1 ? collect(cf[2]-5*se[2]:10*se[2]/20:cf[2]+5*se[2]) : yc

#     if p == 1
#         objv = Array(Float64, length(xc))
#         for j in enumerate(xc)
#             r = j[1]
#             θ = j[2]
#             mdp.e.mm.S[:] = m.e.mf.s(θ)
#             solve!(mdp)
#             if mdp.m.inner.status == 0
#                 v = mdp.m.inner.obj_val[1]
#             else
#                 v = NaN
#             end
#             k1 = κ₁(m)
#             k2 = κ₂(m)
#             S  = bw(m)
#             objv[r] = k2*v/(S*k1^2)
#         end
#         plot(x = xc, y = objv, Geom.line,
#              Guide.xlabel("θ"), Guide.ylabel("P(θ)"))
#     else
#         objv = Array(Float64, (length(xc), length(yc)))
#         for j in enumerate(yc)
#             i  = j[1]
#             θ₂ = j[2]
#             for h in enumerate(xc)
#                 s  = h[1]
#                 θ₁ = h[2]
#                 mdp.e.mm.S[:] = m.e.mf.s([θ₁; θ₂])
#                 solve!(mdp)
#                 if mdp.m.inner.status == 0
#                     v = mdp.m.inner.obj_val[1]
#                 else
#                     v = NaN
#                 end
#                 k1 = MomentBasedEstimators.κ₁(m)
#                 k2 = MomentBasedEstimators.κ₂(m)
#                 S  = MomentBasedEstimators.bw(m)
#                 objv[s,i] = k2*v/(S*k1^2)
#             end
#         end
#         plot(z = objval, Geom.contour)
#         # fig = figure(figsize=(8,6))
#         # ax = fig[:gca](projection="3d")
#         # xgrid, ygrid = meshgrid(xc, yc)
#         # ax[:plot_surface](xgrid, ygrid, objv,
#         #                   rstride=2, cstride=2,
#         #                   cmap=ColorMap("magma"),
#         #                   alpha=0.7, linewidth=0.25)
#         # ax[:contour](xgrid, ygrid, objv, zdir = "z")
#     end
# end


text:
# Inspired by "Fizz Buzz in Tensorflow" blog by Joel Grus
# http://joelgrus.com/2016/05/23/fizz-buzz-in-tensorflow/
using Flux: Chain, Dense, params, logitcrossentropy, onehotbatch, ADAM, train!, softmax
using Test

# Data preparation
function fizzbuzz(x::Int)
    is_divisible_by_three = x % 3 == 0
    is_divisible_by_five = x % 5 == 0

    if is_divisible_by_three & is_divisible_by_five
        return "fizzbuzz"
    elseif is_divisible_by_three
        return "fizz"
    elseif is_divisible_by_five
        return "buzz"
    else
        return "else"
    end
end

const LABELS = ["fizz", "buzz", "fizzbuzz", "else"];

# Feature engineering
features(x) = float.([x % 3, x % 5, x % 15])
features(x::AbstractArray) = hcat(features.(x)...)

function getdata()

    @test fizzbuzz.([3, 5, 15, 98]) == LABELS

    raw_x = 1:100;
    raw_y = fizzbuzz.(raw_x);

    X = features(raw_x);
    y = onehotbatch(raw_y, LABELS);
    return X, y
end

function train()

    # Get Data
    X, y = getdata()

    # Model
    m = Chain(Dense(3, 10), Dense(10, 4))
    loss(x, y) = logitcrossentropy(m(x), y)

    # Helpers
    deepbuzz(x) = (a = argmax(m(features(x))); a == 4 ? x : LABELS[a])

    function monitor(e)
    	print("epoch $(lpad(e, 4)): loss = $(round(loss(X,y); digits=4)) | ")
        @show deepbuzz.([3, 5, 15, 98])
    end

    opt = ADAM()
    # Training
    for e in 0:500
        train!(loss, params(m), [(X, y)], opt)
        if e % 50 == 0
            monitor(e)
        end
    end
end

cd(@__DIR__)
train()

max_stars_repo_name: shawn-kern/model-zoo

max_issues_repo_name: shawn-kern/model-zoo

max_forks_repo_name: shawn-kern/model-zoo


text:
function test_aux()
  # test Givens reflector corner cases
  (c, s, ρ) = Krylov.sym_givens(0.0, 0.0)
  @test (c == 1.0) && (s == 0.0) && (ρ == 0.0)

  a = 3.14
  (c, s, ρ) = Krylov.sym_givens(a, 0.0)
  @test (c == 1.0) && (s == 0.0) && (ρ == a)
  (c, s, ρ) = Krylov.sym_givens(-a, 0.0)
  @test (c == -1.0) && (s == 0.0) && (ρ == a)

  b = 3.14
  (c, s, ρ) = Krylov.sym_givens(0.0, b)
  @test (c == 0.0) && (s == 1.0) && (ρ == b)
  (c, s, ρ) = Krylov.sym_givens(0.0, -b)
  @test (c == 0.0) && (s == -1.0) && (ρ == b)

  # test roots of a quadratic
  roots = Krylov.roots_quadratic(0.0, 0.0, 0.0)
  @test length(roots) == 1
  @test roots[1] == 0.0

  roots = Krylov.roots_quadratic(0.0, 0.0, 1.0)
  @test length(roots) == 0

  roots = Krylov.roots_quadratic(0.0, 3.14, -1.0)
  @test length(roots) == 1
  @test roots[1] == 1.0 / 3.14

  roots = Krylov.roots_quadratic(1.0, 0.0, 1.0)
  @test length(roots) == 0

  roots = Krylov.roots_quadratic(1.0, 0.0, 0.0)
  @test length(roots) == 2
  @test roots[1] == 0.0
  @test roots[2] == 0.0

  roots = Krylov.roots_quadratic(1.0, 3.0, 2.0)
  @test length(roots) == 2
  @test roots[1] ≈ -2.0
  @test roots[2] ≈ -1.0

  roots = Krylov.roots_quadratic(1.0e+8, 1.0, 1.0)
  @test length(roots) == 0

  # ill-conditioned quadratic
  roots = Krylov.roots_quadratic(-1.0e-8, 1.0e+5, 1.0, nitref=0)
  @test length(roots) == 2
  @test roots[1] == 1.0e+13
  @test roots[2] == 0.0

  # iterative refinement is crucial!
  roots = Krylov.roots_quadratic(-1.0e-8, 1.0e+5, 1.0, nitref=1)
  @test length(roots) == 2
  @test roots[1] == 1.0e+13
  @test roots[2] == -1.0e-05

  # not ill-conditioned quadratic
  roots = Krylov.roots_quadratic(-1.0e-7, 1.0, 1.0, nitref=0)
  @test length(roots) == 2
  @test isapprox(roots[1],  1.0e+7, rtol=1.0e-6)
  @test isapprox(roots[2], -1.0, rtol=1.0e-6)

  roots = Krylov.roots_quadratic(-1.0e-7, 1.0, 1.0, nitref=1)
  @test length(roots) == 2
  @test isapprox(roots[1], 1.0e+7, rtol=1.0e-6)
  @test isapprox(roots[2], -1.0, rtol=1.0e-6)

  # test trust-region boundary
  x = ones(5)
  d = ones(5); d[1:2:5] .= -1
  @test_throws ErrorException Krylov.to_boundary(x, d, -1.0)
  @test_throws ErrorException Krylov.to_boundary(x, d, 0.5)
  @test_throws ErrorException Krylov.to_boundary(x, zeros(5), 1.0)
  @test maximum(Krylov.to_boundary(x, d, 5.0)) ≈ 2.209975124224178
  @test minimum(Krylov.to_boundary(x, d, 5.0)) ≈ -1.8099751242241782
  @test maximum(Krylov.to_boundary(x, d, 5.0, flip=true)) ≈ 1.8099751242241782
  @test minimum(Krylov.to_boundary(x, d, 5.0, flip=true)) ≈ -2.209975124224178
end

test_aux()


text:
module TestTypeCheck
  using TypeCheck, FactCheck

  istype(t) = isa(t,TypeCheck.AType)

  facts("Check Return Types: Make Sure It Runs on Base") do
    for n in names(Base)
      if isdefined(Base,n)
        f = eval(Base,n)
        if isgeneric(f) && typeof(f) == Function
          context(string(n)) do
            @fact TypeCheck.checkreturntypes(f) => anything # => FunctionSignature([],Symbol)
            [@fact istype(TypeCheck.returntype(e)) => true for e in code_typed(f)]
            [@fact TypeCheck.returntype(e) => istype for e in code_typed(f)]
          end
        end
      else
        @fact n => x->isdefined(Base,x)
      end
    end
  end

  caught(x) = x[2] == true
  notcaught(x) = x[2] == false
  function checkreturn(f,check)
    @fact length(code_typed(f)) => 1
    @fact TypeCheck.checkreturntype(code_typed(f)[1]) => check
  end

  facts("Check Return Types: True Positives") do
    barr(x::Int) = isprime(x) ? x : false
    checkreturn(barr, caught)
  end

  facts("Check Return Types: False Negatives") do
    foo(x::Any) = isprime(x) ? x : false
    checkreturn(foo, notcaught)
  end

  facts("Check Loop Types: Make Sure It Runs on Base") do
    for n in names(Base)
      if isdefined(Base,n)
        f = eval(Base,n)
        if isgeneric(f) && typeof(f) == Function
          context(string(n)) do
            @fact TypeCheck.checklooptypes(f) => anything # => LoopResults
          end
        end
      else
        @fact n => x->isdefined(Base,x)
      end
    end
  end

  passed(x) = isempty(x.methods)
  failed(x) = !passed(x)
  function checkloops(f,check)
    @fact length(code_typed(f)) => 1
    @fact checklooptypes(f) => check
  end

  facts("Check Loop Types: True Positives") do
    function f1(x::Int)
      for n in 1:x
        x /= n
      end
      return x
    end
    checkloops(f1,failed)
  end

  facts("Check Loop Types: True Negatives") do
    function g1()
      x::Int = 5
      for i = 1:100
        x *= 2.5
      end
      return x
    end
    checkloops(g1,passed)
    function g2()
      x = 5
      x = 0.2
      for i = 1:10
        x *= 2
      end
      return x
    end
    checkloops(g2,passed)
  end

  facts("Check Method Calls: Make Sure It Runs on Base") do
    for n in names(Base)
      if isdefined(Base,n)
        f = eval(Base,n)
        if isgeneric(f) && typeof(f) == Function
          context(string(n)) do
            @fact TypeCheck.checkmethodcalls(f) => anything # => FunctionCalls
          end
        end
      else
        @fact n => x->isdefined(Base,x)
      end
    end
  end

  facts("Check Method Calls: True Positives") do
  end

  facts("Check Method Calls: False Negatives") do
  end

  facts("find_lhs_variables") do
    function lhs1(x::Int)
      x = 5
      return y
    end
    @fact TypeCheck.find_lhs_variables(code_typed(lhs1,(Int,))[1]) => Set{Symbol}([:x])

    function lhs2()
      y = 5
      true && (x = 42)
      foo(y) && (z = 55)
      return y
    end
    @fact TypeCheck.find_lhs_variables(code_typed(lhs2,())[1]) => Set{Symbol}([:x,:y,:z])
  end

  facts("find_rhs_variables") do
    function rhs1(x::Int)
      x = 5
      return y
    end
    @fact TypeCheck.find_rhs_variables(code_typed(rhs1,(Int,))[1]) => Set{Symbol}([:y])

    function rhs2()
      y = 5
      true && (x = 42)
      foo(y) && (z = 55)
      return y
    end
    @fact TypeCheck.find_rhs_variables(code_typed(rhs2,())[1]) => Set{Symbol}([:y])

    function rhs3()
      y = 5
      true && (x = y)
      foo(y) && (z += y)
      return w
    end
    @fact TypeCheck.find_rhs_variables(code_typed(rhs3,())[1]) => Set{Symbol}([:w,:y,:z])
  end

  facts("unused_locals") do
    function fl1(x::String)
      x *= "hello"
      y = x + "hi"
      return x
    end
    @fact TypeCheck.unused_locals(code_typed(fl1,(String,))[1]) => Set{Symbol}([:y])

    function fl2(x)
      answer = false
      if (x > 10)
        anwser = true
      end
      return answer
    end
    @fact TypeCheck.unused_locals(code_typed(fl2,(Int,))[1]) => Set{Symbol}([:anwser])
  end

  exitstatus()
end


text:
using IntervalArithmetic
if VERSION < v"0.7.0-DEV.2004"
    using Base.Test
else
    using Test
end

@test IntervalArithmetic.round_expr(:(a + b), RoundDown) == :($(Expr(:escape, :a)) + $(Expr(:escape, :b)) + $(RoundDown))

@test IntervalArithmetic.round_expr(:(sin(a)), RoundUp) ==
    :( sin( $(Expr(:escape, :a)), $(RoundUp) ) )


text:
# Structs
struct ChebyshevGrid{T <: AbstractFloat, P} <: ApproximationGrid
    a::T
    b::T
    n::Int
    p::P
end

# Constructors
function ChebyshevGrid(a::A, b::B, n) where {A, B}
    T = promote_type(float(A), float(B))
    p = FFTW.plan_r2r!(fill(zero(T), n), FFTW.REDFT10)
    return ChebyshevGrid{T, typeof(p)}(a, b, n, p)
end

# Nodes
function nodes(g::ChebyshevGrid{T}) where {T}
    n = g.n
    a = T(0.5) * (g.b - g.a)
    b = T(0.5) * (g.b + g.a)
    return [a * sinpi((2k + one(T) - n) / 2n) + b for k = 0:(n - 1)]
end

# Weights
function weights(g::ChebyshevGrid)
    n = g.n
    a = zeros(n)
    a[1:2:n] = (2 / n) ./ (1 .- ((1:2:n) .- 1).^2)
    return FFTW.r2r(a, FFTW.REDFT01)
end

# Jacobian
jacobian(g::ChebyshevGrid{T}) where {T} = T(0.5) * (g.b - g.a)

# Interpolation
function interpolate(g::ChebyshevGrid, v::AbstractVector, x′)
    @assert (g.a ≤ x′ ≤ g.b)
    x = 2 * (x′ - (g.b + g.a) / 2) / (g.b - g.a)

    n = length(v)

    c = reverse_plan(g.p, v, n)
    c[1] /= 2
    c   ./= n

    if n == 0
        return zero(eltype(c))
    elseif n == 1 # avoid issues with NaN x
        return first(c) * one(x)
    end

    x = 2x
    b₁ = b₂ = zero(eltype(c))
    @inbounds for k = n:-1:2
        b₂, b₁ = b₁, muladd(x, b₁, c[k] - b₂)
    end

    return muladd(x / 2, b₁, c[1] - b₂)
end

reverse_plan(p, v::AbstractVector{Float64}, n) = p * reverse(v)
function reverse_plan(p, v::AbstractVector{S}, n) where {S}
    r = reinterpret(eltype(S), reverse(v))
    m = length(r)
    d = div(m, n)
    w = permutedims(reshape(r, d, n))

    for i in 1:d
        p * view(w, :, i) # `p` is an in-place `FFTW` plan which mutates `w`
        r[i:d:m] = view(w, :, i)
    end

    return reinterpret(S, r)
end


text:
"""
Runs `ispc` on `code` and prints the resulting assembly to stdout.
"""
function ispc_native(code, options=``)
    ispc_cmd = `ispc -o - --emit-asm $options -`
    open(ispc_cmd, "w", STDOUT) do stdin
        write(stdin, code)
    end
end

"""
Runs `ispc` on `code` and prints the resulting LLVM IR to stdout.
"""
function ispc_llvm(code, options=``, llvm_dis="/usr/local/opt/llvm/bin/llvm-dis")
    ispc_cmd = pipeline(`ispc -o - --emit-llvm $options -`, `$llvm_dis -`)
    open(ispc_cmd, "w", STDOUT) do stdin
        write(stdin, code)
    end
end

# Detect how to link a shared library:
@osx_only begin
    libtool = "/usr/bin/libtool"
    println("Linker: $libtool")
    link(objfile, libfile) = run(`$libtool -dynamic -o "$libfile" "$objfile"`)
end

@linux_only begin
    gpp = strip(readall(`which g++`))
    gpp == "" && error("libtool or g++ is required")
    println("Linker: $gpp")
    link(objfile, libfile) = run(
        `$gpp -shared -Wl,-export-dynamic "$objfile" -o "$libfile"`)
end

@windows_only begin
    error("Not implemented: don't know how to create a shared lib on Windows")
end

_keep_files = false
_debug_generated = false

function keep_files(flag, debug=false)
    global _keep_files
    global _debug_generated
    _keep_files = flag
    _debug_generated = debug
end

"""
Runs `ispc` on `code`, creates a shared library and returns a library
handle `lib` for use with `Libdl.dlsym(lib, symbol)`.
"""
function load_ispc(code, options=``)

    function load(tmpdir)
        objfile = "$tmpdir/program.o"
        libfile = "$tmpdir/program.so"

        # Pipe the input program to ispc:
        ispc_cmd = `ispc -o "$objfile" --pic $options -`
        println(ispc_cmd)
        open(ispc_cmd, "w", STDERR) do stdin
            write(stdin, code)
        end

        if _debug_generated
            sourcefile = "$tmpdir/program.c"
            asmfile = "$tmpdir/program.s"
            bcfile = "$tmpdir/program.bc"
            open(sourcefile, "w") do f
                write(f, code)
            end
            run(`ispc -o $asmfile --emit-asm $options $sourcefile`)
            run(`ispc -o $bcfile --emit-llvm $options $sourcefile`)
        end

        # Create a shared library:
        link(objfile, libfile)
        return Libdl.dlopen(libfile)
    end

    if _keep_files
        load(mktempdir())
    else
        load(mktempdir())
#        mktempdir(load)
        # Note: the temp dir and library files will be deleted now.
        # That's OK because the library is already loaded.
    end
end


text:
using Printf

vandercorput(num::Integer, base::Integer) = sum(d * Float64(base) ^ -ex for (ex, d) in enumerate(digits(num, base)))

for base in 2:9
    @printf("%10s %i:", "Base", base)
    for num in 0:9 @printf("%7.3f", vandercorput(num, base)) end
    println(" [...]")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode