text:
include("prelude.jl")

# Test simple matches
@testset "simple" begin
    interact(`printf 'hello\nworld\n'`, 1) do proc
        @require expect!(proc, "\n") == "hello"
        @require proc.match == "\n"
        @require expect!(proc, "\n") == "world"
        @require proc.match == "\n"
    end
end

# Test match/before
@testset "match" begin
    interact(`printf 'a\nbB'`, 1) do proc
        @require expect!(proc, [r"[AB]"]) == 1
        @require proc.match == "B"
        @require proc.before == "a\nb"
        @test_throws ExpectEOF expect!(proc, ["test"])
        @test_throws ExpectEOF expect!(proc, "test")
        @test_throws ExpectEOF expect!(proc, r"test")
        @require proc.match == nothing
    end
end

# expect! with coprocess
@testset "async" begin
    interact(`cat`, 1) do proc
        write(proc, "hello\nworld\n")
        @require expect!(proc, ["hello\n", "world\n"]) == 1
        @require proc.match == "hello\n"
        @require expect!(proc, ["hello\n", "world\n"]) == 2
        @require proc.match == "world\n"
        @test_throws ExpectTimeout expect!(proc, ["test"])
        @require proc.match == nothing
        @require process_running(proc)
    end
end


text:
using RDatasets

function col_normalize(X::Array{Float64})
    xmin = minimum(X, 1)
    xmax = maximum(X, 1)
    for i in 1:size(X, 2)
        X[:, i] = (X[:, i] .- xmin[i]) ./ (xmax[i] - xmin[i])
    end
    X
end

function get_iris()
    iris = dataset("datasets", "iris")

    iris_data = Array{Float64}(hcat(iris[:SepalLength], iris[:SepalWidth],
                                    iris[:PetalWidth], iris[:PetalLength]));
    X = col_normalize(iris_data)'

    sps = unique(iris[:Species])
    iris[:label] = indexin(iris[:Species], sps)
    Y = Array{Int64}(iris[:label])
    X, Y
end

function get_spirals()
    spirals = readtable("data/spiral.txt", header=false, separator=' ',
                        names=[:x, :y, :label])
    X = col_normalize(Array{Float64}(hcat(spirals[:x], spirals[:y])))'
    Y = Array{Int64}(spirals[:label])
    X, Y
end

function get_yeast()
    yeast = readtable("data/yeast.data", header=false, separator=' ',
                      names=[:seq, :mcg, :gvh, :alm, :mit, :erl, :pox, :vac,
                             :nuc, :local])
    yeast_data = Array{Float64}(hcat(
        map(i->yeast[i],[:mcg, :gvh, :alm, :mit, :erl, :pox, :vac, :nuc])...))
    X = col_normalize(yeast_data)'
    locs = unique(yeast[:local])
    yeast[:label] = indexin(yeast[:local], locs)
    Y = Array{Int64}(yeast[:label])
    X, Y
end

function get_data(problem::String)
    if problem == "iris"
        return get_iris()
    elseif problem == "spirals"
        return get_spirals()
    elseif problem == "yeast"
        return get_yeast()
    end
    error("Problem not found")
end

function train_test_split(X::Array{Float64}, Y::Array{Int64}, trainp::Float64=0.75)
    train = rand(length(Y)) .< trainp
    test = .~(train)
    X[:, train], Y[train], X[:, test], Y[test]
end

function onehot(Y::Array{Int64})
    oneh = zeros(Int64, length(Y), length(unique(Y)))
    for i in eachindex(Y)
        oneh[i, Y[i]] = 1
    end
    oneh
end

ext: jl


text:
module Archer_etal_2000
#=
This module serves to load the 6-box model matrix and grid.
The 6-box model is used because it makes it easier to build
operators for sinking particles in a regular grid,
and because the code to bin WOA data requires a regular grid.
These are the indices of the 6-box model:

┌───────┬────────────────────────┐
│ 1(H)  │         2(L)           │
├───────┼────────────────────────┤
│ 3(H)  │         4(D)           │
├───────┼────────────────────────┤
│       │                        │
│ 5(D)  │         6(D)           │
│       │                        │
└───────┴────────────────────────┘

The boxes are grouped as indicated:
- 1 + 3 represents the high-latitudes (surface) box
- 2 represents the low-latitudes (surface) box
- 4 + 5 + 6 represents the deep box
=#

using LinearAlgebra, SparseArrays
using ..GridTools
using Unitful, UnitfulAstro # for units
using Reexport
@reexport using OceanGrids            # To store the grid
using ..CirculationGeneration
CG = CirculationGeneration

build_wet3D() = trues(2, 1, 3)

function build_grid()
    wet3D = build_wet3D()
    # From Archer et al. [2000]
    vtot = 1.292e18u"m^3"
    area_ocean = 349e12u"m^2"
    depth_low = 100.0u"m"
    depth_high = 250.0u"m"
    fraction_area_high = 0.15 #
    # Infer other values
    depth_ocean = vtot / area_ocean
    # Build DY3d
    elon = [0,360] * u"°"
    elat = [-90, -90+fraction_area_high*180, 90] * u"°"
    edepth = [0u"m", depth_low, depth_high, depth_ocean]
    return OceanGrid(elon, elat, edepth)
end

function build_T(grid)
    # From Archer et al. [2000]
    v3D = array_of_volumes(grid)
    nb = length(v3D)

    # Mixing terms
    # between top boxes (1 and 2) (10 Sv = 10e6 cubic meters)
    T  = CG.flux_divergence_operator_from_advection(10e6u"m^3/s", [1, 2], v3D, nb)
    # between high-lat and deep (1 and 3)
    T += CG.flux_divergence_operator_from_advection(53e6u"m^3/s", [2, 5], v3D, nb)
    # between low-lat and deep (2 and 4)
    T += CG.flux_divergence_operator_from_advection( 1e6u"m^3/s", [2, 4], v3D, nb)
    # (trick) fast mixing between boxes of the same 3-box-model box
    T += CG.flux_divergence_operator_from_advection(1e10u"m^3/s", [1, 3], v3D, nb)
    T += CG.flux_divergence_operator_from_advection(1e10u"m^3/s", [4, 5], v3D, nb)
    T += CG.flux_divergence_operator_from_advection(1e10u"m^3/s", [5, 6], v3D, nb)
    T += CG.flux_divergence_operator_from_advection(1e10u"m^3/s", [4, 6], v3D, nb)

    # Overturning circulation
    T += CG.flux_divergence_operator_from_advection(19e6u"m^3/s", [1, 3, 5, 6, 4, 2], v3D, nb)

    return T
end

"""
    load

Returns wet3D, grid, and T (in that order).
"""
function load()
    print("Creating 6-box model")
    wet3D = build_wet3D()
    grid = build_grid()
    T = build_T(grid)
    println(" ✔")
    @info """
            You are about to use the 3-box model of Archer et al. [2000].
            If you use it for research, please cite:

            - Archer, D. E., Eshel, G., Winguth, A., Broecker, W., Pierrehumbert, R., Tobis, M., and Jacob, R. (2000), Atmospheric pCO2 sensitivity to the biological pump in the ocean, Global Biogeochem. Cycles, 14 (4), 1219–1230, doi:10.1029/1999GB001216.

            You can find the corresponding BibTeX entry in the CITATION.bib file
            at the root of the AIBECS.jl package repository.
            (Look for the "Archer_etal_2000" key.)

            Note that although this model represents the 3-box model of Archer et al. [2000],
            it effectively requires 6 boxes to represent in the AIBECS,
            so that particulate sinking flux divergence operators can be built.
            (See the comments at the start of the Archer_etal_2000.jl file for details.)
            """
    return wet3D, grid, T
end

end

export Archer_etal_2000


text:
using Molly
using DelimitedFiles

data_dir = normpath(dirname(pathof(Molly)), "..", "data")
ff_dir = joinpath(data_dir, "force_fields")
openmm_dir = joinpath(data_dir, "openmm_6mrr")

ff = OpenMMForceField(joinpath.(ff_dir, ["ff99SBildn.xml", "tip3p_standard.xml", "his.xml"])...)

atoms, atoms_data, specific_inter_lists, general_inters, neighbor_finder, coords, box_size = setupsystem(
    joinpath(data_dir, "6mrr_equil.pdb"), ff)

n_steps = 500
velocities = SVector{3}.(eachrow(readdlm(joinpath(openmm_dir, "velocities_300K.txt"))))u"nm * ps^-1"
simulator = VelocityVerlet(dt=0.0005u"ps", coupling=AndersenThermostat(300.0u"K", 1.0u"ps"))

s = System(
    atoms=atoms,
    atoms_data=atoms_data,
    general_inters=general_inters,
    specific_inter_lists=specific_inter_lists,
    coords=coords,
    velocities=velocities,
    box_size=box_size,
    neighbor_finder=neighbor_finder,
)

parallel = false
simulate!(s, simulator, 5; parallel=parallel)
@time simulate!(s, simulator, n_steps; parallel=parallel)

parallel = true
simulate!(s, simulator, 5; parallel=parallel)
@time simulate!(s, simulator, n_steps; parallel=parallel)


text:
using NLsolve
using ForwardDiff
using LinearAlgebra
using Printf
using Plots
include("../src/algorithms/linesearches.jl")

# ELEMENTAL FUNCTIONS

function gfun!(F, x, u, p)

  # retrieve variables
  VM3 = x[1]
  VA3 = x[2]
  VA2 = x[3]
  P1  = x[4]

  VM1 = u[1]
  P2 = u[2]
  VM2 = u[3]

  VA1 = p[1]
  P3 = p[2]
  Q3 = p[3]

  # intermediate quantities
  VA23 = VA2 - VA3
  VA31 = VA3 - VA1
  VA32 = VA3 - VA2
  VA13 = VA1 - VA3

  F[1] = 4.0*VM2*VM2 + VM2*VM3*(-4*cos(VA23) + 10*sin(VA23)) - P2
  F[2] = (8.0*VM3*VM3 + VM3*VM1*(-4*cos(VA31) + 5*sin(VA31))
          + VM3*VM2*(-4*cos(VA32) + 10*sin(VA32)) + P3)
  F[3] = (15.0*VM3*VM3 + VM3*VM1*(-4*sin(VA31) - 5*cos(VA31))
          + VM3*VM2*(-4*sin(VA32) - 10*cos(VA32)) + Q3)
  F[4] = 4.0*VM1*VM1 + VM1*VM3*(-4*cos(VA13) + 5*sin(VA13)) - P1
end

function gfun(x, u, p, T=typeof(x))
  # Get Float64 of Vectors{Float64}, that is the first parameter
  F = zeros(T.parameters[1], 4)

  # retrieve variables
  VM3 = x[1]
  VA3 = x[2]
  VA2 = x[3]
  P1  = x[4]

  VM1 = u[1]
  P2 = u[2]
  VM2 = u[3]

  VA1 = p[1]
  P3 = p[2]
  Q3 = p[3]

  # intermediate quantities
  VA23 = VA2 - VA3
  VA31 = VA3 - VA1
  VA32 = VA3 - VA2
  VA13 = VA1 - VA3

  F[1] = 4.0*VM2*VM2 + VM2*VM3*(-4*cos(VA23) + 10*sin(VA23)) - P2
  F[2] = (8.0*VM3*VM3 + VM3*VM1*(-4*cos(VA31) + 5*sin(VA31))
          + VM3*VM2*(-4*cos(VA32) + 10*sin(VA32)) + P3)
  F[3] = (15.0*VM3*VM3 + VM3*VM1*(-4*sin(VA31) - 5*cos(VA31))
          + VM3*VM2*(-4*sin(VA32) - 10*cos(VA32)) + Q3)
  F[4] = 4.0*VM1*VM1 + VM1*VM3*(-4*cos(VA13) + 5*sin(VA13)) - P1
  return F
end

function cfun(x, u, p)

  VM3 = x[1]
  VA3 = x[2]
  P1  = x[4]

  VM1 = u[1]
  P2 = u[2]

  VA1 = p[1]

  # we fix generation weights inside the
  # function to simplify the script and
  # follow the paper closely.
  w1 = 1.0
  w2 = 1.0

  cost = P1 + P2

  return cost
end

# OPF COMPUTATION

function solve_pf(x, u, p, verbose=true)

  fun_pf!(F, x) = gfun!(F, x, u, p)
  x0 = copy(x)
  res = nlsolve(fun_pf!, x0)
  if verbose
    show(res)
    println("")
  end

  return res.zero

end

# initial parameters
x = zeros(4)
u = zeros(3)
p = zeros(3)

# this is an initial guess
x[1] = 1.0 #VM3
x[2] = 0.0 #VA3
x[3] = 0.0 #VA2
x[4] = 0.5 #P1

# this is given by the problem data, but might be "controlled" via OPF
u[1] = 1.0 #VM1
u[2] = 1.7 #P2
u[3] = 1.0 #VM2

# these parameters are fixed through the computation
p[1] = 0.0 #VA1, slack angle
p[2] = 2.0 #P3
p[3] = 1.0 #Q3

# print initial guesses
println(x)
println(u)

# copy to iteration vectors
xk = copy(x)
uk = copy(u)

iterations = 0

for i = 1:1000
  global xk
  global uk
  println("Iteration ", i)

  # solve power flow
  println("Solving power flow")
  xk = solve_pf(xk, uk, p, false)

  # jacobian
  gx_x(x) = gfun(x, uk, p, typeof(x))
  gx = x -> ForwardDiff.jacobian(gx_x, x)

  # gradient
  cfun_x(x) = cfun(x, uk, p)
  fx = x ->ForwardDiff.gradient(cfun_x, x)

  # lamba calculation
  println("Computing Lagrange multipliers")
  lambda = -inv(gx(xk)')*fx(xk)

  # compute g_u, g_u
  cfun_u(u) = cfun(xk, u, p)
  fu = u ->ForwardDiff.gradient(cfun_u, u)
  gx_u(u) = gfun(xk, u, p, typeof(u))
  gu = u -> ForwardDiff.jacobian(gx_u, u)

  # compute gradient of cost function
  Lu = u -> cfun_u(u) + gx_u(u)'*lambda
  grad_Lu = u -> (fu(u) + gu(u)'*lambda)
  grad_L = grad_Lu(uk)
  println("Norm of gradient ", norm(grad_L))

  # step
  println("Computing new control vector")
  c_par = 1.0
  # Optional linesearch
  c_par = ls(uk, grad_L, Lu, grad_Lu)
  @show c_par
  uk = uk - c_par*grad_L

  #@printf("VM3 %3.2f. VA3 %2.2f. VA2 %2.2f.\n", xk[1], xk[2], xk[3])
  #@printf("VM1 %3.2f. P2 %2.2f. VM2 %2.2f.\n", uk[1], uk[2], uk[3])
end


text:
using ArgParse

##################################################
# Parsing of the command-line arguments
##################################################
tabargs = ArgParseSettings()
@add_arg_table! tabargs begin
    "--parallel"
    help = "Parallel computation: yes/no"
    arg_type = String
    default = "yes"
    "--nbK"
    help = "Number of sampling points for the 3D Rosenbluth integrals"
    arg_type = Int64
    default = 100
    "--nbAvr"
    help = "Number of sampling points for the orbit-averaging integral"
    arg_type = Int64
    default = 100
    "--q"
    help = "Anisotropy parameter q for the Plummer model."
    arg_type = Float64
    default = 1.0
    ##################################################
    # Options to be removed in GitHub
    ##################################################
    "--gAni"
    help = "Anisotropy parameter gAni for the Isochrone model."
    arg_type = Float64
    default = 0.0
    "--logJrcutoff"
    help = "log10(Jr) lower cut-off"
    arg_type = Float64
    default = -Inf
    "--logLcutoff"
    help = "log10(L) lower cut-off"
    arg_type = Float64
    default = -Inf
end
parsed_args = parse_args(tabargs)

max_stars_repo_name: KerwannTEP/CAT

max_issues_repo_name: KerwannTEP/CAT

max_forks_repo_name: KerwannTEP/CAT


text:
#ДАНО: Робот - рядом с горизонтальной перегородкой (под ней), бесконечно продолжающейся в обе стороны,
#в которой имеется проход шириной в одну клетку.
#РЕЗУЛЬТАТ: Робот - в клетке под проходом
function peregorodka(r::Robot)
    side = 3
    counter = 1

    while true
        if !isborder(r, HorizonSide(0))
            break
        end
        for _ in 1:counter
            move!(r, HorizonSide(side))
        end
        side = (side+2)%4
        counter += 1
    end
end

max_stars_repo_name: listova-nastya/mirea-progs

max_issues_repo_name: listova-nastya/mirea-progs

max_forks_repo_name: listova-nastya/mirea-progs


text:
using AutoGrad: AutoGrad, @primitive1, value
import CUDA.CUDNN: cudnnConvolutionForwardAD
using CUDA.CUDNN:
    CUDNN_ACTIVATION_IDENTITY,
    CUDNN_ACTIVATION_RELU,
    CUDNN_NOT_PROPAGATE_NAN,
    cudnnActivationBackward,
    cudnnActivationDescriptor,
    cudnnConvolutionBackwardBias,
    cudnnConvolutionBackwardData,
    cudnnConvolutionBackwardFilter,
    cudnnConvolutionBwdDataAlgoPerf,
    cudnnConvolutionBwdFilterAlgoPerf,
    scalingParameter,
    with_workspace,
    handle

# Define gradients
@primitive1(
    (cudnnConvolutionForwardAD(w, x, bias, z; y, activation, convDesc, wDesc, xDesc, yDesc, zDesc, biasDesc, alpha, beta, dw, dx, dz, dbias, dready),_dy,_y),
    (dready[] || cudnnConvolutionBackward(_dy, _y, w, x, bias, z; y, activation, convDesc, wDesc, xDesc, yDesc, zDesc, biasDesc, alpha, beta, dw, dx, dz, dbias, dready); dw[]),
    (dready[] || cudnnConvolutionBackward(_dy, _y, w, x, bias, z; y, activation, convDesc, wDesc, xDesc, yDesc, zDesc, biasDesc, alpha, beta, dw, dx, dz, dbias, dready); dx[]),
    (dready[] || cudnnConvolutionBackward(_dy, _y, w, x, bias, z; y, activation, convDesc, wDesc, xDesc, yDesc, zDesc, biasDesc, alpha, beta, dw, dx, dz, dbias, dready); dbias[]),
    (dready[] || cudnnConvolutionBackward(_dy, _y, w, x, bias, z; y, activation, convDesc, wDesc, xDesc, yDesc, zDesc, biasDesc, alpha, beta, dw, dx, dz, dbias, dready); dz[]))

function cudnnConvolutionBackward(_dy, _y, w, x, bias, z; y, activation, convDesc, wDesc, xDesc, yDesc, zDesc, biasDesc, alpha, beta, dw, dx, dz, dbias, dready)
    # Make sure backward gets called only once
    dready[] && return
    dready[] = true
    # Read all relevant inputs
    (dy, y, w, x, bias, z) = value.((_dy, _y, w, x, bias, z))
    # Allocate gradient buffers if necessary
    if dw[] === nothing; dw[] = similar(w); end; @assert issimilar(w, dw[])
    if dx[] === nothing; dx[] = similar(x); end; @assert issimilar(x, dx[])
    if bias === nothing
        dbias[] = nothing
    elseif dbias[] === nothing
        dbias[] = similar(bias)
    end
    @assert issimilar(bias, dbias[])
    # Calculate pre-activation gradient if necessary, use dz[] for storage
    @assert isa(z,Nothing) || issimilar(z, y)
    if activation !== CUDNN_ACTIVATION_IDENTITY
        if dz[] === nothing; dz[] = similar(y); end; @assert issimilar(y, dz[])
        actback!(dz[],y,dy,y,activation)
        dy = dz[]
    end
    beta0, alpha1 = (a->scalingParameter(eltype(x),a)).((0,1))
    p = cudnnConvolutionBwdFilterAlgoPerf(xDesc, x, yDesc, dy, convDesc, wDesc, dw[]);
    with_workspace(p.memory) do workspace
        cudnnConvolutionBackwardFilter(handle(), alpha, xDesc, x, yDesc, dy, convDesc, p.algo, workspace, sizeof(workspace), beta0, wDesc, dw[])
    end
    p = cudnnConvolutionBwdDataAlgoPerf(wDesc, w, yDesc, dy, convDesc, xDesc, dx[]);
    with_workspace(p.memory) do workspace
        cudnnConvolutionBackwardData(handle(), alpha, wDesc, w, yDesc, dy, convDesc, p.algo, workspace, sizeof(workspace), beta0, xDesc, dx[])
    end
    if bias !== nothing
        cudnnConvolutionBackwardBias(handle(), alpha1, yDesc, dy, beta0, biasDesc, dbias[])
    end
    if beta[] == 0 || z === nothing
        dz[] = nothing
    elseif dz[] === nothing
        dz[] = beta[] .* dy
    else
        dz[] .= beta[] .* dy
    end
end

function actback!(dx,x,dy,y,activation)
    activationDesc = cudnnActivationDescriptor(activation, CUDNN_NOT_PROPAGATE_NAN, Cfloat(1))
    xDesc, yDesc = cudnnTensorDescriptor.((x,y))
    alpha, beta = (x->scalingParameter(eltype(y),x)).((1,0))
    cudnnActivationBackward(handle(), activationDesc, alpha, yDesc, y, yDesc, dy, xDesc, x, beta, xDesc, dx)
    return dx
end


text:
function meshpartition(dmd,p,t,f,t2t,tprd,elemtype,bcm,bndexpr,prdexpr,porder,nproc,metis)

# display("run facenumbering...");
# f, tprd = facenumbering(p,t,elemtype,bndexpr,prdexpr);

print("run elementpartition...\n");
if size(tprd,2)==size(t,1)
    dmd = elementpartition2(dmd,tprd,t2t,nproc,metis);
else
    dmd = elementpartition2(dmd,t,t2t,nproc,metis);
end

print("run facepartition...\n");
dmd = facepartition(dmd,p,t,tprd,f,bcm,elemtype,prdexpr,porder,nproc);

return dmd;

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
module PoissonFourier1D

using FFTW, LinearAlgebra

function generate_solver(ρ₀, Δx)
    PF = FFTW.plan_rfft(ρ₀)
    F = FFTW.rfft(ρ₀)
    ω = collect(2π*rfftfreq(length(ρ₀)))
    ω[1] = ω[2]
    ω⁻² = -1.0./(ω.^2)
    Pφ = FFTW.plan_irfft(F, length(ρ₀))
    φ = FFTW.irfft(F, length(ρ₀))
    function solve!(e::Vector{T}, ρ::Vector{T}) where {T}
        mul!(F, PF, ρ)
        F[1] = 0.0 + 0.0im
        for i = 2:length(F)
            F[i] = ω⁻²[i] * F[i]
        end
        mul!(φ, Pφ, F)
        e[1] = 0.5*(φ[2]-φ[end])/Δx
        for i = 2:length(φ)-1
            e[i] = 0.5*(φ[i+1] - φ[i-1])/Δx
        end
        e[end] = 0.5*(φ[1]-φ[end-1])/Δx
        return e
    end
    return solve!
end

end


text:
##############################################################################
#
# Basic DocumentTermMatrix type
#
##############################################################################

mutable struct DocumentTermMatrix
    dtm::SparseMatrixCSC{Int, Int}
    terms::Vector{String}
    column_indices::Dict{String, Int}
end

##############################################################################
#
# Construct a DocumentTermMatrix from a Corpus
#
##############################################################################

# create col index lookup dictionary from a (sorted?) vector of terms
function columnindices(terms::Vector{String})
    column_indices = Dict{String, Int}()
    for i in 1:length(terms)
        term = terms[i]
        column_indices[term] = i
    end
    column_indices
end

function DocumentTermMatrix(crps::Corpus, terms::Vector{String})
    column_indices = columnindices(terms)

    m = length(crps)
    n = length(terms)

    rows = Array{Int}(undef, 0)
    columns = Array{Int}(undef, 0)
    values = Array{Int}(undef, 0)
    for i in 1:m
        doc = crps.documents[i]
        ngs = ngrams(doc)
        for ngram in keys(ngs)
            j = get(column_indices, ngram, 0)
            v = ngs[ngram]
            if j != 0
                push!(rows, i)
                push!(columns, j)
                push!(values, v)
            end
        end
    end
    if length(rows) > 0
        dtm = sparse(rows, columns, values, m, n)
    else
        dtm = spzeros(Int, m, n)
    end
    DocumentTermMatrix(dtm, terms, column_indices)
end
DocumentTermMatrix(crps::Corpus) = DocumentTermMatrix(crps, lexicon(crps))

DocumentTermMatrix(crps::Corpus, lex::AbstractDict) = DocumentTermMatrix(crps, sort(collect(keys(lex))))

DocumentTermMatrix(dtm::SparseMatrixCSC{Int, Int},terms::Vector{String}) = DocumentTermMatrix(dtm, terms, columnindices(terms))

##############################################################################
#
# Access the DTM of a DocumentTermMatrix
#
##############################################################################

function dtm(d::DocumentTermMatrix, density::Symbol)
    if density == :sparse
        return d.dtm
    else
        return Matrix(d.dtm)
    end
end

function dtm(d::DocumentTermMatrix)
    return d.dtm
end

function dtm(crps::Corpus)
    dtm(DocumentTermMatrix(crps))
end

tdm(crps::DocumentTermMatrix, density::Symbol) = dtm(crps, density)' #'

tdm(crps::DocumentTermMatrix) = dtm(crps)' #'

tdm(crps::Corpus) = dtm(crps)' #'

##############################################################################
#
# Produce the signature of a DTM entry for a document
#
# TODO: Make this more efficieny by reusing column_indices
#
##############################################################################

function dtm_entries(d::AbstractDocument, lex::Dict{String, Int})
    ngs = ngrams(d)
    indices = Array{Int}(undef, 0)
    values = Array{Int}(undef, 0)
    terms = sort(collect(keys(lex)))
    column_indices = columnindices(terms)

    for ngram in keys(ngs)
        if haskey(column_indices, ngram)
            push!(indices, column_indices[ngram])
            push!(values, ngs[ngram])
        end
    end
    return (indices, values)
end

function dtv(d::AbstractDocument, lex::Dict{String, Int})
    p = length(keys(lex))
    row = zeros(Int, 1, p)
    indices, values = dtm_entries(d, lex)
    for i in 1:length(indices)
        row[1, indices[i]] = values[i]
    end
    return row
end

function dtv(d::AbstractDocument)
    error("Cannot construct a DTV without a pre-existing lexicon")
end

##############################################################################
#
# The hash trick: use a hash function instead of a lexicon to determine the
# columns of a DocumentTermMatrix-like encoding of the data
#
##############################################################################

function hash_dtv(d::AbstractDocument, h::TextHashFunction)
    p = cardinality(h)
    res = zeros(Int, 1, p)
    ngs = ngrams(d)
    for ng in keys(ngs)
        res[1, index_hash(ng, h)] += ngs[ng]
    end
    return res
end

hash_dtv(d::AbstractDocument) = hash_dtv(d, TextHashFunction())

function hash_dtm(crps::Corpus, h::TextHashFunction)
    n, p = length(crps), cardinality(h)
    res = zeros(Int, n, p)
    for i in 1:length(crps)
        doc = crps.documents[i]
        res[i, :] = hash_dtv(doc, h)
    end
    return res
end

hash_dtm(crps::Corpus) = hash_dtm(crps, hash_function(crps))

hash_tdm(crps::Corpus) = hash_dtm(crps)' #'

##############################################################################
#
# Produce entries for on-line analysis when DTM would not fit in memory
#
##############################################################################

mutable struct EachDTV
    crps::Corpus
end

start(edt::EachDTV) = 1

function next(edt::EachDTV, state::Int)
    return (dtv(edt.crps.documents[state], lexicon(edt.crps)), state + 1)
end

done(edt::EachDTV, state::Int) = state > length(edt.crps.documents)

mutable struct EachHashDTV
    crps::Corpus
end

start(edt::EachHashDTV) = 1

function next(edt::EachHashDTV, state::Int)
    (hash_dtv(edt.crps.documents[state]), state + 1)
end

done(edt::EachHashDTV, state::Int) = state > length(edt.crps.documents)

each_dtv(crps::Corpus) = EachDTV(crps)

each_hash_dtv(crps::Corpus) = EachHashDTV(crps)

##
## getindex() methods
##

Base.getindex(dtm::DocumentTermMatrix, k::AbstractString) = dtm.dtm[:, dtm.column_indices[k]]
Base.getindex(dtm::DocumentTermMatrix, i::Any) = dtm.dtm[i]
Base.getindex(dtm::DocumentTermMatrix, i::Any, j::Any) = dtm.dtm[i, j]


text:
# module Bukdu

import HttpCommon: Request, Response

module Server

function info()::String
    "Bukdu (commit $(commit_short())) with Julia $VERSION"
end

function commit_short()::String
    repo = LibGit2.GitRepo(Pkg.dir("Bukdu"))
    string(LibGit2.revparseid(repo, "HEAD"))[1:7]
end

end # module Bukdu.Server

function handler(req::Request, res::Response)
    if method_exists(before, (Request,Response))
        before(req, res)
    end
    verb = getfield(Bukdu, Symbol(lowercase(req.method)))
    local conn::Conn
    try
        conn = Routing.request(RouterRoute.routes, verb, req.resource, req.data) do route
            Base.function_name(route.verb) == Base.function_name(verb)
        end
    catch ex
        isa(ex, NoRouteError)
        conn = CONN_NOT_FOUND
    end
    for (key,value) in conn.resp_header
       res.headers[key] = value
    end
    res.headers["Server"] = Server.info()
    res.status = conn.status
    if isa(conn.resp_body, Vector{UInt8}) || isa(conn.resp_body, String)
       res.data = conn.resp_body
    else
       res.data = string(conn.resp_body)
    end
    if method_exists(after, (Request,Response))
        after(req, res)
    end
    res
end

module Farm

import HttpServer
servers = Vector{Tuple{HttpServer.Server,Task}}()

end # module Bukdu.Farm

import HttpServer

"""
    Bukdu.start(port::Int; host=getaddrinfo("localhost"))

Start Bukdu server with port.

```jula
julia> Bukdu.start(8080)
Listening on 127.0.0.1:8080...
```
"""
function start(port::Int; host=getaddrinfo("localhost"))
    start([port]; host=host)
end

"""
    Bukdu.start(ports::Vector{Int}; host=getaddrinfo("localhost"))

Start Bukdu server with multiple ports.

```jula
julia> Bukdu.start([8080, 8081])
Listening on 127.0.0.1:8080...
```
"""
function start(ports::Vector{Int}; host=getaddrinfo("localhost"))
    for port in ports
        server = HttpServer.Server(handler)
        server.http.events["listen"] = (port) -> Logger.info("Listening on $port..."; LF=!isdefined(:Juno))
        task = @async begin
            HttpServer.run(server, host=host, port=port)
        end
        if :queued == task.state
            push!(Farm.servers, (server,task))
        end
    end
end

"""
    Bukdu.stop()

Stop the Bukdu server.
"""
function stop()
    stopped = 0
    for (server,task) in Farm.servers
        try
            if Base.StatusActive == server.http.sock.status
                stopped += 1
            end
            close(server)
        catch e
        end
    end
    if stopped >= 1
        Logger.info("Stopped.")
        empty!(Farm.servers)
    end
    nothing
end


text:
mutable struct Nodes
    nodes::Array{Tuple{Integer, Integer}, 1}

    function Nodes(s::String)
        tiles = Tiles(s).data
        height, width = size(tiles)

        nodesRaw = Tuple{Char, Tuple{Integer, Integer}}[
            (tiles[y, x], (x * 8 - 4,  y * 8 - 4)) for y in 1:height, x in 1:width if tiles[y, x] != '0' && tiles[y, x] != ' '
        ]
        sort!(nodesRaw, by=(v -> v[1]))

        return new([node[2] for node in nodesRaw])
    end
end

max_stars_repo_name: CelestialCartographers/Maple

max_issues_repo_name: CelestialCartographers/Maple

max_forks_repo_name: CelestialCartographers/Maple


text:
# Copyright (c) 2019 Daniel Moura
# https://towardsdatascience.com/freeing-the-data-scientist-mind-from-the-curse-of-vectorization-11634c370107?source=friends_link&sk=e7b0086f3678261cb64b6dbb43744613
# linkedin.com/in/dmoura
# twitter.com/daniel_c_moura

# running this file:
# julia --optimize=3 --check-bounds=no distances_among_pairs.jl

using CPUTime
using LinearAlgebra

### functions under evaluation ###

function distances_among_pairs_loop(v)
    nin = length(v)
    nout = div(nin * nin - nin, 2) #length of the output vector
    dists = Vector{eltype(v)}(undef,nout) #preallocating output vector
    k = 1 #current position in the output vector
    @inbounds for i in 1:(nin-1)
        a = v[i]
        @inbounds for j in (i+1):nin
            b = v[j]
            dists[k] = abs(a-b)
            k += 1
        end
    end
    dists
end

distances_among_pairs_comp(v) =
    [abs(v[i]-v[j]) for i in eachindex(v), j in eachindex(v) if i<j]

function distances_among_pairs_outer(v)
    dists = v .- v'
    abs.(dists[tril!(trues(size(dists)), -1)])
end

### aux functions ###

function bench(label, f, v)
    GC.gc() #cleaning the garbage before each measurement
    print("$label:\n  ")
    @time @CPUtime f(v)
    println()
end

function readvec(run)
    open("vec$run.csv", "r") do f
        map(line -> parse(Float64, line), eachline(f))
    end
end

### main ###

function main()
    println("Warming up...")
    v = fill(1.0, 3)
    distances_among_pairs_loop(v)
    distances_among_pairs_comp(v)
    distances_among_pairs_outer(v)

    for r in 1:3
        println("\n----- RUN $r -----")
        v = readvec(r)
        bench("LOOP",           distances_among_pairs_loop,  v)
        bench("COMPREHENSION",  distances_among_pairs_comp,  v)
        bench("OUTER",          distances_among_pairs_outer, v)
    end
end

main()

max_stars_repo_name: dcmoura/blogposts

max_issues_repo_name: dcmoura/blogposts

max_forks_repo_name: dcmoura/blogposts


text:
"""
    fit(ts,Δ;model::Type{<:TimeSeriesModel},x₀,lowerΩcutoff,upperΩcutoff,x_lowerbounds,x_upperbounds,method,taper)
    fit(timeseries::TimeSeries;model::Type{<:TimeSeriesModel},x₀,lowerΩcutoff,upperΩcutoff,x_lowerbounds,x_upperbounds,method,taper)

Fit a time series model using the `IPNewton` method from Optim.jl.

# Arguments
- `ts`: `n` by `D` matrix containing the timeseries (or vector if `D=1`), where `n` is the number of observations and `D` is the number of series.
- `Δ`: The sampling rate, which should be a positive real number.
- `timeseries`: Can be provided in place of `ts` and `Δ`.
- `model`: The model which will be fitted. Should be a type (not a realisation of the model) e.g. JONSWAP{k} not JONSWAP{K}(x).
- `x₀`: The initial parameter guess.
- `lowerΩcutoff`: The lower cutoff for the frequency range to be used in fitting. Default is `0`.
- `upperΩcutoff`: The upper cutoff for the frequency range to be used in fitting. Default is `Inf`.
- `x_lowerbounds`: The lower bounds on the parameter space. If `nothing` is provided (the default) then these are set to default values based on the model.
- `x_upperbounds`: The upper bounds on the parameter space. If `nothing` is provided (the default) then these are set to default values based on the model.
- `method`: Either `:Whittle` or `:debiasedWhittle`.
- `taper`: The choice of tapering to be used. This should be `nothing` (in which case no taper is used) or `dpss_nw` where `nw` time-bandwith product (see DSP.dpss for more details).
- `options`: Options for optimisation of type `Optim.Options`.
"""
function fit(ts::VecOrMat{Float64},Δ::Real;model::Type{<:TimeSeriesModel},x₀,lowerΩcutoff::Real=0.0,upperΩcutoff::Real=Inf,
            x_lowerbounds=nothing,x_upperbounds=nothing,method=:debiasedWhittle,taper=nothing,options::Optim.Options=Optim.Options(iterations=1000))
    Δ > 0 || throw(ArgumentError("Δ should be a positive number."))
    model(x₀) # to check that the model can be constructed from the initial vector provided.
    size(ts,2) == ndims(model) || throw(ArgumentError("ts should by of size `(n,D)` where `D` is the dimension of the timeseries."))
    lowerΩcutoff < upperΩcutoff || throw(ArgumentError("lower frequency cutoff should be below upper frequency cutoff."))
    if x_lowerbounds === nothing
        x_lowerbounds = lowerbounds(model)
    else
        length(x_lowerbounds) == npars(model) || throw(ArgumentError("provided vector of lower bounds is not the same length as the number of parameters"))
    end
    if x_upperbounds === nothing
        x_upperbounds = upperbounds(model)
    else
        length(x_upperbounds) == npars(model) || throw(ArgumentError("provided vector of upper bounds is not the same length as the number of parameters"))
    end
    if method == :Whittle
        objective = WhittleLikelihood(model, ts, Δ; lowerΩcutoff = lowerΩcutoff, upperΩcutoff = upperΩcutoff, taper = maketaper(taper,size(ts,1)))
    elseif method == :debiasedWhittle
        objective = DebiasedWhittleLikelihood(model, ts, Δ; lowerΩcutoff = lowerΩcutoff, upperΩcutoff = upperΩcutoff, taper = maketaper(taper,size(ts,1)))
    else
        throw(ArgumentError("Method not recognised. Choose from :Whittle or :debiasedWhittle"))
    end
    constraints = Optim.TwiceDifferentiableConstraints(x_lowerbounds,x_upperbounds)
    obj = TwiceDifferentiable(Optim.only_fgh!(objective),x₀)
    res = optimize(obj, constraints, x₀, IPNewton(), options)
    Optim.converged(res) || @warn "Optimisation did not converge!"
    return res
end
function fit(ts::TimeSeries;model::Type{<:TimeSeriesModel},x₀,lowerΩcutoff::Real=0.0,upperΩcutoff::Real=Inf,
            x_lowerbounds=nothing,x_upperbounds=nothing,method=:debiasedWhittle,taper=nothing,options::Optim.Options=Optim.Options(iterations=1000))
    fit(ts.ts,ts.Δ,model=model,x₀=x₀,lowerΩcutoff=lowerΩcutoff,upperΩcutoff=upperΩcutoff,
    x_lowerbounds=x_lowerbounds,x_upperbounds=x_upperbounds,method=method,taper=taper,options=options)
end
maketaper(taper::Nothing,n) = nothing
maketaper(taper::Symbol,n) = maketaper(String(taper),n)
function maketaper(taper::String,n)
    tapername, nw = split(taper,'_')
    tapername == "dpss" || throw(ArgumentError("Taper should be nothing or dpss_nw where nw is the value of nw used for the taper. Other kinds of tapering are not yet supported."))
    return vec(dpss(n, parse(Float64,nw), 1))
end
lowerbounds(::Type{<:AdditiveTimeSeriesModel{M₁,M₂,D,T}}) where {M₁,M₂,D,T} = vcat(lowerbounds(M₁),lowerbounds(M₂))
upperbounds(::Type{<:AdditiveTimeSeriesModel{M₁,M₂,D,T}}) where {M₁,M₂,D,T} = vcat(upperbounds(M₁),upperbounds(M₂))


text:
################################################################################
###  Introduction to Scientific Programming and Machine Learning with Julia  ###
###                                                                          ###
### Run each script on a new clean Julia session                             ###
### GitHub: https://github.com/sylvaticus/IntroSPMLJuliaCourse               ###
### Licence (apply to all material of the course: scripts, videos, quizes,..)###
### Creative Commons By Attribution (CC BY 4.0), Antonello Lobianco          ###
################################################################################

# # 0106 Further Topics

# ## Some stuff to set-up the environment..

cd(@__DIR__)
using Pkg
Pkg.activate(".")
## If using a Julia version different than 1.7 please uncomment and run the following line (reproductibility guarantee will however be lost)
## Pkg.resolve()
## Pkg.instantiate() # run this if you didn't in Segment 01.01
using Random
Random.seed!(123)
using InteractiveUtils # loaded automatically when working... interactively

# ## Metaprogramming and macros

# "running" some code include the following passages (roughly):
# - parsing of the text defining the code and its translation in hierarchical expressions to the Abstract syntax Tree (AST) (syntax errors are caugth at this time)
# - on the first instance required ("just in time") compilation of the AST expressions into object code (using the LLVM compiler)
# - execution of the compiled object code

# "Macros" in many other language (e.g. C or C++) refer to the possibility to "pre-process" the textual representation of the code statements before it is parsed. In julia instead it refers to the possibility to alter the expression once has already being parsed in the AST, allowing a greater expressivity as we are no longer limited by the parsing syntax

# The AST is organised in a hierarchical tree of _expressions_ where each element (including the operators) is a _symbol_
# For variables, you can use symbols to refer to the actual identifiers instad to the variable's value

# Expressions themselves are objects representing unevaluated computer expressions

# ### Expressions and symbols

expr1 = Meta.parse("a = b + 2") # What the parser do when reading the source code. b doesn't need to actually been defined, it's just a namebinding without the reference to any object, not even `nothing`
typeof(expr1) # expressions are first class objects
expr2 = :(a = b + 1)
expr3 = quote a = b + 1 end
expr3
dump(expr1)   # The AST ! Note this is already a nested statement, an assignment of the result of an expression (the sum call between the symbol `:b` and 1) to the symbol `a`
expr4 = Expr(:(=),:a,Expr(:call,:+,:b,1)) # The AST using the "Expr" constructor

symbol1   = :(a)               # as for expressions
symbol2   = Meta.parse("a")    # as for expressions
symbol3   = Symbol("a")        # specific for symbols only
othsymbol = Symbol("aaa",10,"bbb")
typeof(symbol1)

# I can access any parts of my expression before evaluating it (indeed, that's what macro will do...)
myASymbol = expr1.args[2].args[1]
expr1.args[2].args[1] = :(*)
b = 2
## a # error, a not defined
eval(expr1)
a # here now is defined and it has an object associated... 4!

# !!! danger
#     The capability to evaluate expressions is very powerfull but due to obvious secutiry implications never evaluate expressions you aren't sure of their provenience. For example if you develop a Julia web app (e.g. using [Genie.jl](https://github.com/GenieFramework/Genie.jl)) never evaluate user provided expressions.

# Note that evaluation of expressions happens always at global scope, even if it done inside a function:
function foo()
    locVar = 1
    expr = :(locVar + 1)
    return eval(expr)
end
## a = foo() # error locVar not defined

# To refer to the _value_ of a variable rather than the identifier itself within an expression, interpolate the variable using the dollar sign:

expr = :($a + b) # here the identifier 'a' has been replaced with its numerical value, `4`
dump(expr)
eval(expr)
a = 10
eval(expr) # no changes
b = 100
eval(expr) # here it change, as it is at eval time that the identifier `b` is "replaced" with its value

# ### Macros

# One of the best usage of macros is they allow package developers to provide a very flexible API to their package, suited for the specific needs of the package, making life easier for the users. Compare for example [the API of JuMP](https://jump.dev/JuMP.jl/stable/manual/constraints/) with [those of Pyomo](https://pyomo.readthedocs.io/en/stable/pyomo_modeling_components/Constraints.html) to define model constraints !

# Some examples of macros:
# [MultiDimEquations.jl](https://github.com/sylvaticus/MultiDimEquations.jl)
# - from: `@meq par1[d1 in DIM1, d2 in DIM2, dfix3] =  par2[d1,d2]+par3[d1,d2]`
# - to:   `[par1[d1,d2,dfix3] =  par2[d1,d2]+par3[d1,d2] for d1 in DIM1, d2 in DIM2]`
# [Pipe.jl:](https://github.com/oxinabox/Pipe.jl)
# - from: `@pipe  10 |> foo(_,a) |> foo2(b,_,c) |> foo3(_)`
# - to:   `foo3(foo2(b,foo(10,a),c))`
# Brodcasting (Base):
# - from: `@. a + b * D^2`
# - to:   `a .+ b .* D.^2`

# Defining a macro...

# Like functions, but both the arguments and the returned output are expressions
macro customLoop(controlExpr,workExpr)
    return quote
      for i in $controlExpr
        $workExpr
      end
    end
end

# Invoking a macro....
a = 5
@customLoop 1:4 println(i) #note that "i" is in the macro
@customLoop 1:a println(i)
@customLoop 1:a if i > 3 println(i) end
@customLoop ["apple", "orange", "banana"]  println(i)
@customLoop ["apple", "orange", "banana"]  begin print("i: "); println(i)  end
@macroexpand @customLoop 1:4 println(i) # print what the macro does with the specific expressions provided

# String macros (aka "non-standard string literals")
# Invoked with the syntax `xxx" ...text..."` or `xxx""" ...multi-line text..."""` where `xxx` is the name of the macro and the macro must be defined as `macro  xxx_str`.
# Used to perform textual modification o nthe given text, for example this print the given text on a 8 characters

macro print8_str(mystr)                                 # input here is a string, not an expression
    limits = collect(1:8:length(mystr))
    for (i,j) in enumerate(limits)
      st = j
      en = i==length(limits) ? length(mystr) : j+7
      println(mystr[st:en])
    end
end

print8"123456789012345678"
print8"""This is a text that once printed in 8 columns with terminal will be several lines. Ok, no rammar rules relating to carriage returns are emploied here..."""

# While normally used to modify text, string macros are "true" macros:

macro customLoop_str(str)
    exprs = Meta.parse(str)
    controlExpr,workExpr = exprs.args[1],exprs.args[2]
    return quote
      for i in $controlExpr
        $workExpr
      end
    end
end

customLoop"""1:4; println(i)"""

# ## Interfacing with other languages

# There are 3 ways to interface Julia with programs or libraries wrote in other languages.
# At the lowest level, Julia allows to directly interface with C or Fortran libraries, and this means, aside using directly libraries written in C, to be able to interface with any programming language that offer also a C interface (R, Python...)
# Using this low level C Interface, users have created specific packages to interface many languages using a simple, Julian-way syntax. We will see these interfaces for R and Python.
# Finally, at the highest level, many common packages of other languages have been already "interfaced", so that the user can use the Julia Package without even knowing that this is an interface for an other package, for example `SymPy.jl` is a large interface to the Python package `SymPy`.

# ### Using C libraries

# Let's start by seing how to use a C library. For this example to work you will need to have the GCC compiler installed on your machine
# First let's write the header and source C files and write them to the disk:

cheader = """
extern int get5();
extern double mySum(float x, float y);
"""
csource = """
int get5(){
    return 5;
}

double mySum(float x, float y){
    return x+y;
}
"""

open(f->write(f,cheader),"myclib.h","w")  # We open a stream to file with the "w" parameter as for "writing", and we pass the stream to the anonymous function to actually write to the stream. If this funcitons is many lines of code, consider rewriting the `open` statement using a `do` block
open(f->write(f,csource),"myclib.c","w")

# Now let's run the command to compile the C code we saved as shared library using gcc, a C compiler.
# The following example assume that GCC is installed in the machine where this example is run and available as `gcc`.

compilationCommand1 = `gcc -o myclib.o -c myclib.c` # the actual compilation, note the backticks used to define a command
compilationCommand2 = `gcc -shared -o libmyclib.so myclib.o -lm -fPIC` # the linking into a shared library
run(compilationCommand1)
run(compilationCommand2)

# This should have created the C library `libmyclib.so` on disk. Let's gonna use it:
const myclib = joinpath(@__DIR__, "libmyclib.so")  # we need the full path
# ccall arguments:
# 1. A tuple with the funcion name to call and the library path. For both, if embedded in a variable, the variable must be set constant.
# 2. The Julia type that map to the C type returned by the function.
#    - `int` → `Int32` or `Int64` (or the easy-to remmeber `Cint` alias)
#    - `float` → `Float32` (or the `Cfloat` alias)
#    - `double` → `Float64` (or the `Cdouble` alias)
# 3. A tuple with the Julia types of the parameters passed to the C function
# 4. Any other argument are the values of the parameter passed
a = ccall((:get5,myclib), Int32, ())
b = ccall((:mySum,myclib), Float64, (Float32,Float32), 2.5, 1.5)

# More details on calling C or Fortran code can be obtained [in the official Julia documentation](https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/).

# ### Using Python in Julia

# The "default" way to use Python code in Julia is trough the [PyCall.jl](https://github.com/JuliaPy/PyCall.jl) package. It automatically take care of convert between Python types (including numpy arrays) and Julia types (types that can not be converted automatically are converted to the generic `PyObject` type).
ENV["PYTHON"] = "" # will force PyCall to download and use a "private to Julia" (conda based) version of Python. use "/path/to/python" if you want to reuse a version already installed on your system
## using Pkg
## Pkg.add("PyCall")
## Pkg.build("PyCall")
using PyCall

# #### Embed short python snippets in Julia
py"""
def sumMyArgs (i, j):
  return i+j
def getNthElement (vec,n):
  return vec[n]
"""
a = py"sumMyArgs"(3,4)             # 7 - here we call the Python object (a function) with Julia parameters
b = py"getNthElement"([1,2,3],1)   # 2 - attention to the diffferent convention for starting arrays!. Note the Julia Array ahas been converted automatically to a Python list
d = py"getNthElement([1,$a,3],1)"  # 7 - here we interpolate the Python call

# Alternativly, use `@pyinclude("pythonScript.py")`
pythonCode = """
def sumMyArgs (i, j, z):
  return i+j+z
"""
open(f->write(f,pythonCode),"pythonScript.py","w")
@pyinclude("pythonScript.py")
a = py"sumMyArgs"(3,4,5)

# !!! tip
#     Note thaat the 3 arguments definition of `sumMyArgs` has _replaced_ the 3-arguments one. This would now error `py"sumMyArgs"(3,4)`

# #### Use Python libraries

# Add a package to the local Python installation using Conda:
pyimport_conda("ezodf", "ezodf", "conda-forge") # pyimport_conda(module, package, channel)

const ez = pyimport("ezodf")  # Equiv. of Python `import ezodf as ez`
destDoc  = ez.newdoc(doctype="ods", filename="anOdsSheet.ods")
# Both `ez` and `destDoc` are `PyObjects` for which we can access attributes and call the methods using the usual `obj.method()` syntax as we would do in Python
sheet    = ez.Sheet("Sheet1", size=(10, 10))
destDoc.sheets.append(sheet)
## dcell1 = sheet[(2,3)] # This would error because the index is a tuple. Let's use directly the `get(obj,key)` function instead:
dcell1   = get(sheet,(2,3)) # Equiv. of Python `dcell1 = sheet[(2,3)]`. Attention again to Python indexing from zero: this is cell "D3", not "B3" !
dcell1.set_value("Hello")
get(sheet,"A9").set_value(10.5) # Equiv. of Python `sheet['A9'].set_value(10.5)`
destDoc.backup = false
destDoc.save()

# ### Using Julia in Python

# #### Installation of the Python package `PyJulia`

# [PyJulia](https://github.com/JuliaPy/pyjulia) can be installed using `pip`, taking note that its name using `pip` is `julia` not `PyJulia`:
# ```$ python3 -m pip install --user julia```

# We can now open a Python terminal and initialise PyJulia to work with our Julia version:
# ```python
# >>> import julia
# >>> julia.install() # Only once to set-up in julia the julia packages required by PyJulia
# ```

# If we have multiple Julia versions, we can specify the one to use in Python passing julia="/path/to/julia/binary/executable" (e.g. julia = "/home/myUser/lib/julia-1.1.0/bin/julia") to the install() function.

# #### Running Julia libraries and code in Python

# On each Python session we need to run the following code:
# ```python
# from julia import Julia
# Julia(compiled_modules=False)
# ```

# This is a workaround to the common situation when the Python interpreter is statically linked to libpython, but it will slow down the interactive experience, as it will disable Julia packages pre-compilation, and every time we will use a module for the first time, this will need to be compiled first. Other, more efficient but also more complicate, workarounds are given in the package documentation, under the [Troubleshooting section](https://pyjulia.readthedocs.io/en/stable/troubleshooting.html).

# We can now direcltly load a Julia module, including `Main`, the global namespace of Julia’s interpreter, with `from julia import ModuleToLoad` and access the module objects directly or using the `Module.evel()` interface.

# ##### Add a Julia package...
# ```python
# >>> from julia import Pkg
# >>> Pkg.add("BetaML")
# ```
# Of course we can add a package alternatively from within Julia

# ##### "Direct" calling of Julia functions...
# ```python
# >>> from julia import BetaML
# >>> import numpy as np
# >>> model = BetaML.buildForest([[1,10],[2,12],[12,1]],["a","a","b"])
# >>> predictions = BetaML.predict(model,np.array([[2,9],[13,0]]))
# >>> predictions
# [{'b': 0.36666666666666664, 'a': 0.6333333333333333}, {'b': 0.7333333333333333, 'a': 0.26666666666666666}]
# ```

# ##### Access using the `eval()` interface...
# If we are using the jl.eval() interface, the objects we use must be already known to julia. To pass objects from Python to Julia, we can import the julia Main module (the root module in julia) and assign the needed variables, e.g.

# ```python
# >>> X_python = [1,2,3,2,4]
# >>> from julia import Main
# >>> Main.X_julia = X_python
# >>> Main.eval('BetaML.gini(X_julia)')
# 0.7199999999999999
# >>> Main.eval("""
# ...   function makeProd(x,y)
# ...       return x*y
# ...   end
# ...   """
# ... )
# >>> Main.eval("makeProd(2,3)") # or Main.makeProd(2,3)
# ```
# For large scripts instead of using `eval()` we can equivalently use `Main.include("aJuliaScript.jl")`

# ### Using R in Julia

# To use R from within Julia we use the [RCall](https://github.com/JuliaInterop/RCall.jl) package.
ENV["R_HOME"] = "*" #  # will force RCall to download and use a "private to Julia" (conda based) version of R. use "/path/to/R/directory" (e.g. `/usr/lib/R`) if you want to reuse a version already installed on your system
## using Pkg
## Pkg.add("RCall")
## Pkg.build("RCall")
using RCall

R"""
sumMyArgs <- function(i,j) i+j
getNthElement <- function(vec,n) {
  return(vec[n])
}
"""
a = rcopy(R"sumMyArgs"(3,4))             # 7 - here we call the R object (a function) with Julia parameters
b = rcopy(R"getNthElement"([1,2,3],1))   # 1 - no differences in array indexing here
d = rcopy(R"as.integer(getNthElement(c(1,$a,3),2))")  # 7 - here we interpolate the R call
d = convert(Int64,R"getNthElement(c(1,$a,3),2)")

# While we don't have here the problem of different array indexing convention (both Julia and R start indexing arrays at 1), we have the "problem" that the output returned by using `R"..."` is not yet an exploitable Julia object but it remains as an `RObject` that we can convert with `rcopy()` or explicitly with `convert(T,obj)`. Also, R elements are all floats by default, so if we need an integer in Julia we need to explicitly convert it, either in R or in Julia.

# If the R code is on a script, we don't have here a sort of @Rinclude macro, so let's implement it ourselves by loading the file content as a file and evaluating it using the function `reval` provided by `RCall`:

macro Rinclude(fname)
    quote
        rCodeString = read($fname,String)
        reval(rCodeString)
        nothing
    end
end

RCode = """
sumMyArgs <- function(i, j, z) i+j+z
"""
open(f->write(f,RCode),"RScript.R","w")
@Rinclude("RScript.R")
a = rcopy(R"sumMyArgs"(3,4,5))  # 12
## a = rcopy(R"sumMyArgs"(3,4)) # error !  The 3-arguments version of `sumMyArgs` has _replaced_ the 2-arguments one

# ### Using Julia in R

# #### Installation of the R package `JuliaCall`

# [JuliaCall](https://github.com/Non-Contradiction/JuliaCall) can be installed from CRAN:

# ```{r}
# > install.packages("JuliaCall")
# > library(JuliaCall)
# install_julia()
# ```

# `install_julia()` will force the download of R and install a private copy of julia. If you prefer to use instead an existing version of julia and having R default to download a private version only if it can't find a version already installed, use `julia_setup(installJulia = TRUE)` instead of `install_julia()`, eventually passing the `JULIA_HOME = "/path/to/julia/binary/executable/directory"` (e.g. `JULIA_HOME = "/home/myUser/lib/julia-1.7.0/bin"`) parameter to the `julia_setup` call.

# `JuliaCall` depends for some things (like object conversion between Julia and R) from the Julia `RCall` package. If we don't already have it installed in Julia, it will try to install it automatically.

# #### Running Julia libraries and code in R

# On each R session we need to run the `julia_setup` function:
# ```R
# library(JuliaCall)
# julia_setup() # If we have already downloaded a private version of Julia for R it will be retrieved automatically
# ```

# We can now load a Julia module and access the module objects directly or using the `Module.evel()` interface.

# ##### Add a Julia package...

# ```{r}
# > julia_eval('using Pkg; Pkg.add("BetaML")')
# ```
# Of course we can add a package alternatively from within Julia

# Let's load some data from R and do some work with this data in Julia:

# ```{r}
# > library(datasets)
# > X <- as.matrix(sapply(iris[,1:4], as.numeric))
# > y <- sapply(iris[,5], as.integer)
# ```

# ##### Calling of Julia functions with `julia_call`...

# With `JuliaCall`, differently than `PyJulia`, we can't call direclty the julia functions but we need to employ the R function `julia_call("juliaFunction",args)`:

# ```{r}
# > julia_eval("using BetaML")
# > yencoded <- julia_call("integerEncoder",y)
# > ids      <- julia_call("shuffle",1:length(y))
# > Xs       <- X[ids,]
# > ys       <- yencoded[ids]
# > cOut     <- julia_call("kmeans",Xs,3L)    # kmeans expects K to be an integer
# > y_hat    <- sapply(cOut[1],as.integer)[,] # We need a vector, not a matrix
# > acc      <- julia_call("accuracy",y_hat,ys)
# > acc
# [1] 0.8933333
# ```
# ##### Access using the `eval()` interface...

# As alternative, we can embed Julia code directly in R using the `julia_eval()` function:

# ```{r}
# > kMeansR  <- julia_eval('
# +      function accFromKmeans(x,k,y_true)
# +        cOut = kmeans(x,Int(k))
# +        acc = accuracy(cOut[1],y_true)
# +        return acc
# +      end
# + ')
# ```

# We can then call the above function in R in one of the following three ways:
# 1. `kMeansR(Xs,3,ys)`
# 2. `julia_assign("Xs_julia", Xs); julia_assign("ys_julia", ys); julia_eval("accFromKmeans(Xs_julia,3,ys_julia)")`
# 3. `julia_call("accFromKmeans",Xs,3,ys)`.

# While other "convenience" functions are provided by the package, using  `julia_call` or `julia_assign` followed by `julia_eval` should suffix to accomplish most of the task we may need in Julia.

# ## Some performance tips

# ### Type stability

# "Type stable" functions guarantee to the compiler that given a certain method (i.e. with the arguments being of a given type) the object returned by the function is also of a certain fixed type. Type stability is fundamental to allow type inference continue across the function call stack.

function f1(x)    # Type unstable
    outVector = [1,2.0,"2"]
    if x < 0
        return outVector[1]
    elseif x == 0
        return outVector[2]
    else
        return outVector[3]
    end
end

function f2(x)   # Type stable
    outVector = [1,convert(Int64,2.0),parse(Int64,"2")]
    if x < 0
        return outVector[1]
    elseif x == 0
        return outVector[2]
    else
        return outVector[3]
    end
end

a = f1(0)
b = f1(1)
typeof(a)
typeof(b)

c = f2(0)
d = f2(1)
typeof(c)
typeof(d)

using BenchmarkTools
@btime f1(0) # 661 ns 6 allocations
@btime f2(0) #  55 ns 1 allocations

@code_warntype f1(0) # Body::Any
@code_warntype f2(0) # Body::Int64

# While in general is NOT important to annotate function parameters for performance, it is important to annotate struct fields with concrete types
abstract type Goo end
struct Foo <: Goo
    x::Number
end
struct Boo <: Goo
    x::Int64
end

function f1(o::Goo)
    return o.x +2
end

fobj = Foo(1)
bobj = Boo(1)
@btime f1($fobj) # 17.1 ns 0 allocations
@btime f1($bobj) #  2.8 ns 0 allocations

# Here the same function under some argument types is type stable, under other argument types is not
@code_warntype f1(fobj)
@code_warntype f1(bobj)

# #### Avoid (non-constant) global variables

g        = 2
const cg = 1   # we can't change the _type_ of the object binded to a constant variable
cg       = 2   # we can rebind to an other object of the same type, but we get a warning
## cg    = 2.5 # this would error !
f1(x,y) = x+y
f2(x)   = x + g
f3(x)   = x + cg

@btime f1(3,2)
@btime f2(3)    # 22 times slower !!!
@btime f3(3)    # as f1

# #### Loop arrays with the inner loop by rows

# Julia is column mayor (differently than Python) so arrays of bits types are contiguous in memory across the different rows of the same column

a = rand(1000,1000);
function f1(x)
    (R,C) = size(x)
    cum = 0.0
    for r in 1:R
        for c in 1:C
            cum += x[r,c]
        end
    end
    return cum
end
function f2(x)
    (R,C) = size(x)
    cum = 0.0
    for c in 1:C
        for r in 1:R
            cum += x[r,c]
        end
    end
    return cum
end
@btime f1($a) # 2.3 ms 0 allocations
@btime f2($a) # 1.3 ms 0 allocations

# #### Use low-level optimisation when possible

function f1(x)
    s = 0.0
    for i in 1:length(x)
        s += i * x[i]
    end
    return s
end

function f2(x)
    s = 0.0
    for i in 1:length(x)
        @inbounds s += i * x[i] # remove bound checks
    end
    return s
end

function f3(x)
    s = 0.0
    @simd for i in 1:length(x) # tell compiler it is allowed to run the loop in whatever order, allowing in-thread paralllelism of modern CPUs
        s += i * x[i]
    end
    return s
end

x = rand(10000);
@btime f1($x)
@btime f2($x)
@btime f3($x)

X = rand(100,20);
function f1(x)
    s = 0.0
    for i in 1:size(x,1)
        s += sum(x[i,:])
    end
    return s
end
function f2(x)
    s = 0.0
    @views for i in 1:size(x,1)
        s += sum(x[i,:])   # the slice operator copy the data.. the views macro force to have instead to have a view (reference)
    end
    return s
end

@btime f1($X)
@btime f2($X)

# !!! warning
#     Attention that while the `@views` macro "save time" by not copying the data, the resulting array has a pretty messy layout. If you need to use it for many subsequent operations it may be more efficient to "pay" the copy cost once and then have an array with a nicelly continuous block of memory..

function f1(x,y)
    if x+y > 100
        return x + y + 2
    else
        return x + 1
    end
end

@inline function f2(x,y)   # the function is "inlined", its whole definition copied at each calling place rather than being called
    if x+y > 100
        return x + y + 2
    else
        return x + 1
    end
end

function f3(y)
    s = 0.0
    for i in 2:y
       s += f1(i,i-1)
       s += f1(i,i)
    end
    return s
 end

function f4(y)
    s = 0.0
    for i in 2:y
       s += f2(i,i-1)
       s += f2(i,i)
    end
    return s
end

x = 1000
@btime f3($x)
@btime f4($x)

# But attention! Not always a good idea:
function f3(y)
    s = 0.0
    for i in 2:y
       s += sum(f1(a,a-1) for a in 2:i)
    end
    return s
 end

function f4(y)
    s = 0.0
    for i in 2:y
       s += sum(f2(a,a-1) for a in 2:i)
    end
    return s
end

x = 1000
@btime f3($x)
@btime f4($x)

# Note that the Julia compiles already inline small functions automatically when it thinks it will improve performances

# ## Profiling the code to discover bootlenecks

# We already see `@btime` and `@benchmark` from the package [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl)
# Remember to quote the global variables used as parameter of your function with the dollar sign to have accurate benchmarking of the function execution.
# Julia provide the macro `@time` but we should run on a second call to a given function (with a certain parameter types) or it will include compilation time in its output:

function fb(x)
    out = Union{Int64,Float64}[1,2.0,3]
    push!(out,4)
    if x > 10
        if ( x > 100)
            return [out[1],out[2]] |>  sum
        else
            return [out[2],out[3]] |>  sum
        end
    else
        return [out[1],out[3]] |>  sum
    end
end

@time fb(3)
@time fb(3)

# We can use `@profile function(x,y)` to use a sample-based profiling

using Profile # in the stdlib
function foo(n)
    a = rand(n,n)
    b = a + a
    c = b * b
    return c
end
@profile (for i = 1:100; foo(1000); end) # too fast otherwise
Profile.print() # on my pc: 243 rand, 174 the sum, 439 the matrix product
Profile.clear()

# ## Introspection and debugging

# To discover problems on the code more in general we can use several introspection functions that Julia provide us (some of which we have already saw):

## @less rand(3)  # Show the source code of the specific method invoked - use `q` to quit
## @edit rand(3)  # Like @loss but it opens the source code in an editor
methods(foo)
@which foo(2)          # which method am I using when I call foo with an integer?
typeof(a)
eltype(a)
fieldnames(Foo)
dump(fobj)
names(Main, all=false) # available (e.g. exported) identifiers of a given module
sizeof(2)              # bytes
typemin(Int64)
typemax(Int64)
bitstring(2)
# Various low-level interpretation of an expression
@code_native foo(3)
@code_llvm foo(3)
@code_typed foo(3)
@code_lowered foo(3)

# We can use a debugger, like e.g. the one integrated in Juno or VSCode.
# Graphical debuggers allow to put a _breakpoint_ on some specific line of code, run the code in debug mode (yes, it will be slower), let the program arrive to the breakpoint and inspect the state of the system at that point of the code, including local variables. In Julia we can also _change_ the program interactively !
# Other typycal functions are running a single line, running inside a function, running until the current function return, ecc..

# ## Runtime exceptions

# As many (all?) languages, Julia when "finds" an error issues an exception, that if it is not caugth at higher level in the call stack (i.e. recognised and handled) lead to an error and return to the prompt or termination of the script (and rarely with the Julia process crashing altogether).

# The idea is that we _try_ some potentially dangerous code and if some error is raised in this code we _catch_ it and handle it.

function customIndex(vect,idx;toReturn=0)
    try
        vect[idx]
    catch e
        if isa(e,BoundsError)
            return toReturn
        end
        rethrow(e)
    end
end

a = [1,2,3]
## a[4] # Error ("BoundsError" to be precise)
customIndex(a,4)

# Note that handling exceptions is computationally expensive, so do not use exceptions in place of conditional statements

# ## [Parallel computation](@id parallel_computation)

# Finally one note on parallel computation. We see only some basic usage of multithreading and multiprocesses in this course, but with Julia it is relativelly easy to parallelise the code either using multiple threads or multiple processes. What's the difference ?
# - **multithread**
#   - advantages: computationally "cheap" to create (the memory is shared)
#   - disadvantages: limited to the number of cores within a CPU, require attention in not overwriting the same memory or doing it at the intended order ("data race"), we can't add threads dynamically (within a script)
# - **multiprocesses**
#   - advantages: unlimited number, can be run in different CPUs of the same machine or differnet nodes of a cluster, even using SSH on different networks, we can add processes from within our code with `addprocs(nToAdd)`
#   - disadvantages: the memory being copied (each process will have its own memory) are computationally expensive (you need to have a gain higher than the cost on setting a new process) and require attention to select which memory a given process will need to "bring with it" for its functionality

# Note that if you are reading this document on the github pages, this script is compiled using GitHub actions where a single thread and process are available, so you will not see performance gains.

# ### Multithreading

# !!! warning
#     It is not possible to add threads dinamically, either we have to start Julia with the parameter `-t` (e.g. `-t 8` or `-t auto`) in the command line or use the VSCode Julia externsion setting `Julia: Num Threads`

function inner(x)
    s = 0.0
    for i in 1:x
        for j in 1:i
            if j%2 == 0
                s += j
            else
                s -= j
            end
        end
    end
    return s
end

function parentSingleThread(x,y)
    toTest = x .+ (1:y)
    out = zeros(length(toTest))
    for i in 1:length(toTest)
        out[i] = inner(toTest[i])
    end
    return out
end
function parentThreaded(x,y)
    toTest = x .+ (1:y)
    out = zeros(length(toTest))
    Threads.@threads for i in 1:length(toTest)
        out[i] = inner(toTest[i])
    end
    return out
end

x = 100
y = 20

str = parentSingleThread(x,y)
mtr = parentThreaded(x,y)

str == mtr # true
Threads.nthreads() # 4 in my case
Threads.threadid()
@btime parentSingleThread(100,20) # 140 μs on my machine
@btime parentThreaded(100,20)     #  47 μs

# ### Multiprocessing

# _**NOTE**_
# _The code on multiprocessing is commented out (not executed, so you don't see the output) as GitHub actions (that are used to run this code and render the web pages you are reading) have problems running multi-process functions:_

# ```julia
# using Distributed     # from the Standard Library
# addprocs(3)           # 2,3,4
# ```
# The first process is considered a sort of "master" process, the other one are the "workers"
# We can add processes on other machines by providing the SSH connection details directly in the `addprocs()` call (Julia must be installed on that machines as well)
# We can alternativly start Julia directly with _n_ worker processes using the armument `-p n` in the command line.

# ```julia
# println("Worker pids: ")
# for pid in workers()  # return a vector to the pids
#     println(pid)      # 2,3,4
# end
# rmprocs(workers()[1])    #  remove process pid 2
# println("Worker pids: ")
# for pid in workers()
#     println(pid) # 3,4 are left
# end
# @everywhere begin using Distributed end # this is needed only in GitHub action
# @everywhere println(myid()) # 4,3
# ```

# #### Run heavy tasks in parallel

# ```julia
# using Distributed, BenchmarkTools
# a = rand(1:35,100)
# @everywhere function fib(n)
#     if n == 0 return 0 end
#     if n == 1 return 1 end
#     return fib(n-1) + fib(n-2)
# end
# ```

# The macro `@everywhere` make available the given function (or functions with `@everywhere begin [shared function definitions] end` or `@everywhere include("sharedCode.jl")`) to all the current workers.

## result  = map(fib,a)

# The pmap function ("parallel" map) automatically pick up the free processes, assign them the job prom the "input" array and merge the results in the returned array. Note that the order is preserved:

# ```julia
# result2 = pmap(fib,a)
# result == result2
# @btime map(fib,$a)  # serialised:   median time: 514 ms    1 allocations
# @btime pmap(fib,$a) # parallelised: median time: 265 ms 4220 allocations # the memory of `a` need to be copied to all processes
# ```

# #### Divide and Conquer

# Rather than having a "heavy operation" and being interested in the individual results, here we have a "light" operation and we want to aggregate the results of the various computations using some aggreagation function.
# We can then use `@distributed (aggregationfunction) for [forConditions]` macro:

# ```julia
# using Distributed, BenchmarkTools
# function f(n)   # our single-process benchmark
#   s = 0.0
#   for i = 1:n
#     s += i/2
#   end
#     return s
# end
# function pf(n)
#   s = @distributed (+) for i = 1:n # aggregate using sum on variable s
#         i/2                        # the last element of the for cycle is used by the aggregator
#   end
#   return s
# end
# @btime  f(10000000) # median time: 11.1 ms   0 allocations
# @btime pf(10000000) # median time:  5.7 ms 145 allocations
# ```

# Note that also in this case the improvement is less than proportional with the number of processes we add

# Details on parallel comutation can be found [on the official documentation](https://docs.julialang.org/en/v1/manual/parallel-computing/), including information to run nativelly Julia on GPUs or TPUs.

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
test_LP() = CoreModel(
    zeros(4, 3),
    zeros(4),
    ones(3),
    ones(3),
    ones(3),
    ["r$x" for x = 1:3],
    ["m$x" for x = 1:4],
)

test_simpleLP() = CoreModel(
    [
        1.0 1.0
        -1.0 1.0
    ],
    [3.0, 1.0],
    [-0.25, 1.0],
    -ones(2),
    2.0 * ones(2),
    ["r$x" for x = 1:2],
    ["m$x" for x = 1:2],
)

test_simpleLP2() = CoreModel(
    zeros(2, 2),
    [0.0, 0.0],
    [-0.25, 1.0],
    -ones(2),
    2.0 * ones(2),
    ["r$x" for x = 1:2],
    ["m$x" for x = 1:2],
)

test_sparseLP() = CoreModel(
    sprand(4000, 3000, 0.5),
    sprand(4000, 0.5),
    sprand(3000, 0.5),
    sprand(3000, 0.5),
    sprand(3000, 0.5),
    ["r$x" for x = 1:3000],
    ["m$x" for x = 1:4000],
)

test_coupledLP() = CoreModelCoupled(
    CoreModel(
        sprand(4000, 3000, 0.5),
        sprand(4000, 0.5),
        sprand(3000, 0.5),
        sprand(3000, 0.5),
        sprand(3000, 0.5),
        ["r$x" for x = 1:3000],
        ["m$x" for x = 1:4000],
    ),
    sprand(2000, 3000, 0.5),
    sprand(2000, 0.5),
    sprand(2000, 0.5),
)


text:
using Documenter
using Pidfile

# The DOCSARGS environment variable can be used to pass additional arguments to make.jl.
# This is useful on CI, if you need to change the behavior of the build slightly but you
# can not change the .travis.yml or make.jl scripts any more (e.g. for a tag build).
if haskey(ENV, "DOCSARGS")
    for arg in split(ENV["DOCSARGS"])
        (arg in ARGS) || push!(ARGS, arg)
    end
end

makedocs(
    modules = [Pidfile],
    format = Documenter.HTML(
        prettyurls = true,
        canonical = "https://vtjnash.github.io/Pidfile.jl/",
    ),
    sitename = "Pidfile.jl",
    pages = [
        "Home" => "index.md",
    ],
    linkcheck = true,
    doctest = true,
)

#deploydocs(
#    repo = "github.com/vtjnash/Pidfile.jl.git",
#    target = "build",
#    push_preview = true,
#)


text:
# First verify that this Version of julia is compatible
# with this package.

if !isdefined(first(methods(sin)),:isstaged)
    error("This version of julia does not have staged functions.\n
           CXX.jl will not work!")
end

if !isdefined(Base,:llvmcall)
    error("This version of julia does not support llvmcall.\n
           CXX.jl will not work!")
end

println("Tuning for julia installation at: ",JULIA_HOME)

run(`make -f BuildBootstrap.Makefile JULIA_HOME=$JULIA_HOME`)


text:
using OrdinaryDiffEq, Test
function lorenz(du, u, p, t)
    du[1] = 10.0(u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
end
u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz, u0, tspan)
@test_nowarn sol = solve(prob, Tsit5(), reltol=1e-6)
sol = solve(prob, Tsit5(), rel_tol=1e-6)
@test_logs (:warn, DiffEqBase.KWARGWARN_MESSAGE) sol = solve(prob, Tsit5(), rel_tol=1e-6)
@test_throws DiffEqBase.CommonKwargError sol = solve(prob, Tsit5(), rel_tol=1e-6, kwargshandle=DiffEqBase.KeywordArgError)

prob = ODEProblem(lorenz, u0, tspan, test = 2.0)
@test_logs (:warn, DiffEqBase.KWARGWARN_MESSAGE) sol = solve(prob, Tsit5(), reltol=1e-6)


text:
function sparse_from_coo(::Type{<:SymSparseMatrixCSR{Bi}}, I,J,V,m,n) where Bi
  symsparsecsr(Val(Bi),I,J,V,m,n)
end

@inline function is_entry_stored(::Type{<:SymSparseMatrixCSR},i,j)
  (i<=j)
end

function finalize_coo!(T::Type{<:SymSparseMatrixCSR},I,J,V,m,n)
  m == n || throw(DimensionMismatch("matrix is not square: dimensions are ($m, $n)"))
  touched = zeros(Bool,m)
  for k in 1:length(I)
    Ik = I[k]
    Jk = J[k]
    if Ik == Jk
      touched[Ik] = true
    end
  end
  for k in 1:m
    if ! touched[k]
      push!(I,k)
      push!(J,k)
      push!(V,zero(eltype(V)))
    end
  end
end

function nz_index(A::SymSparseMatrixCSR,i,j)
  nz_index(A.uppertrian,i,j)
end

function add_entry!(combine::Function,A::SymSparseMatrixCSR,v::Number,i,j)
  if i<=j
    k = nz_index(A,i,j)
    nz = nonzeros(A)
    Aij = nz[k]
    nz[k] = combine(v,Aij)
  end
  A
end

@inline function push_coo!(::Type{<:SymSparseMatrixCSR},I,J,V,i,j,v)
  (i>j) && return nothing
  push!(I,i)
  push!(J,j)
  push!(V,v)
  nothing
end