text:
#
# @description Implements the DES (Data Encryption Standard) block cipher algorithm.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2015 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
export DesCipher
export DES_BLOCK_SIZE, DES_LEGAL_KEY_SIZES
export block_size, key_size, encrypt, decrypt

# @description Defines the DES cipher data struture.
struct DesCipher <: BlockCipher
  subkeys::Vector{Vector{UInt8}}

  # @description Initialize the DES cipher.
  #
  # @param {Vector{UInt8}} key The symmetric key.
  function DesCipher(key::Vector{UInt8})
    return new(compute_des_subkeys(key))
  end
end

# @description The DES cipher block size in bytes.
const DES_BLOCK_SIZE = UInt64(8)

# @description The DES cipher legal key sizes in bytes.
const DES_LEGAL_KEY_SIZES = UInt64[ 8 ]

# @description Gets the cipher block size.
#
# @param {DesCipher} self The cipher data struture.
#
# @return {Integer} The block size in bytes.
function block_size(self::DesCipher)
  return DES_BLOCK_SIZE
end

# @description Gets the cipher key size.
#
# @param {DesCipher} self The cipher data struture.
#
# @return {Integer} The key size in bytes.
function key_size(self::DesCipher)
  return DES_LEGAL_KEY_SIZES[1]
end

# @description Encrypts the given plaintext block with the key set at initialization.
#
# @param {DesCipher}     self      The cipher data struture.
# @param {Vector{UInt8}} plaintext The message to encrypt.
#
# @return {Vector{UInt8}} The resulting ciphertext block.
function encrypt(self::DesCipher, plaintext::Vector{UInt8})
  if length(plaintext) != block_size(self)
    throw(ArgumentError("Invalid input plaintext length."))
  end

  feistel    = FeistelCipher(compute_des_round, self.subkeys, key_size(self), block_size(self))
  block      = permutate_bits(plaintext, DES_INITIAL_PERMUTATION)
  block      = encrypt(feistel, block)
  ciphertext = permutate_bits(block, DES_FINAL_PERMUTATION)

  return ciphertext
end

# @description Decrypts the given ciphertext block with the key set at initialization.
#
# @param {DesCipher}     self       The cipher data struture.
# @param {Vector{UInt8}} ciphertext The message to decrypt.
#
# @return {Vector{UInt8}} The resulting plaintext block.
function decrypt(self::DesCipher, ciphertext::Vector{UInt8})
  if length(ciphertext) != block_size(self)
    throw(ArgumentError("Invalid input ciphertext length."))
  end

  feistel   = FeistelCipher(compute_des_round, self.subkeys, key_size(self), block_size(self))
  block     = permutate_bits(ciphertext, DES_INITIAL_PERMUTATION)
  block     = decrypt(feistel, block)
  plaintext = permutate_bits(block, DES_FINAL_PERMUTATION)

  return plaintext
end

# PRIVATE IMPLEMENTATION #######################################################

# @description First permuted choice, used in the key schedule to select 56 bits from a 64-bit input.
const DES_PC1 = UInt8[
  57,   49,    41,   33,    25,    17,    9,
   1,   58,    50,   42,    34,    26,   18,
  10,    2,    59,   51,    43,    35,   27,
  19,   11,     3,   60,    52,    44,   36,
  63,   55,    47,   39,    31,    23,   15,
   7,   62,    54,   46,    38,    30,   22,
  14,    6,    61,   53,    45,    37,   29,
  21,   13,     5,   28,    20,    12,    4
]

# @description Schedule of "left shifts" for subkeys construction.
const DES_LEFT_SHIFTS = [ 1, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 1 ]

# @description Second permuted choice, used in the key schedule to produce each subkey.
const DES_PC2 = UInt8[
  14,    17,   11,    24,     1,    5,
   3,    28,   15,     6,    21,   10,
  23,    19,   12,     4,    26,    8,
  16,     7,   27,    20,    13,    2,
  41,    52,   31,    37,    47,   55,
  30,    40,   51,    45,    33,   48,
  44,    49,   39,    56,    34,   53,
  46,    42,   50,    36,    29,   32
]

# @description Initial permutation for the 64-bit input block.
const DES_INITIAL_PERMUTATION = UInt8[
  58,    50,   42,    34,    26,   18,    10,    2,
  60,    52,   44,    36,    28,   20,    12,    4,
  62,    54,   46,    38,    30,   22,    14,    6,
  64,    56,   48,    40,    32,   24,    16,    8,
  57,    49,   41,    33,    25,   17,     9,    1,
  59,    51,   43,    35,    27,   19,    11,    3,
  61,    53,   45,    37,    29,   21,    13,    5,
  63,    55,   47,    39,    31,   23,    15,    7
]

# @description Expansion of input block of 32 bits, producing an output block of 48 bits.
const DES_EXPANSION_FUNCTION = UInt8[
  32,     1,    2,     3,     4,    5,
   4,     5,    6,     7,     8,    9,
   8,     9,   10,    11,    12,   13,
  12,    13,   14,    15,    16,   17,
  16,    17,   18,    19,    20,   21,
  20,    21,   22,    23,    24,   25,
  24,    25,   26,    27,    28,   29,
  28,    29,   30,    31,    32,    1
]

# @description The DES' S-box 1
const DES_SBOX1 = UInt8[
  14   4  13   1   2  15  11   8   3  10   6  12   5   9   0   7;
   0  15   7   4  14   2  13   1  10   6  12  11   9   5   3   8;
   4   1  14   8  13   6   2  11  15  12   9   7   3  10   5   0;
  15  12   8   2   4   9   1   7   5  11   3  14  10   0   6  13
]

# @description The DES' S-box 2
const DES_SBOX2 = UInt8[
  15   1   8  14   6  11   3   4   9   7   2  13  12   0   5  10;
   3  13   4   7  15   2   8  14  12   0   1  10   6   9  11   5;
   0  14   7  11  10   4  13   1   5   8  12   6   9   3   2  15;
  13   8  10   1   3  15   4   2  11   6   7  12   0   5  14   9
]

# @description The DES' S-box 3
const DES_SBOX3 = UInt8[
  10   0   9  14   6   3  15   5   1  13  12   7  11   4   2   8;
  13   7   0   9   3   4   6  10   2   8   5  14  12  11  15   1;
  13   6   4   9   8  15   3   0  11   1   2  12   5  10  14   7;
   1  10  13   0   6   9   8   7   4  15  14   3  11   5   2  12
]

# @description The DES' S-box 4
const DES_SBOX4 = UInt8[
   7  13  14   3   0   6   9  10   1   2   8   5  11  12   4  15;
  13   8  11   5   6  15   0   3   4   7   2  12   1  10  14   9;
  10   6   9   0  12  11   7  13  15   1   3  14   5   2   8   4;
   3  15   0   6  10   1  13   8   9   4   5  11  12   7   2  14
]

# @description The DES' S-box 5
const DES_SBOX5 = UInt8[
   2  12   4   1   7  10  11   6   8   5   3  15  13   0  14   9;
  14  11   2  12   4   7  13   1   5   0  15  10   3   9   8   6;
   4   2   1  11  10  13   7   8  15   9  12   5   6   3   0  14;
  11   8  12   7   1  14   2  13   6  15   0   9  10   4   5   3
]

# @description The DES' S-box 6
const DES_SBOX6 = UInt8[
  12   1  10  15   9   2   6   8   0  13   3   4  14   7   5  11;
  10  15   4   2   7  12   9   5   6   1  13  14   0  11   3   8;
   9  14  15   5   2   8  12   3   7   0   4  10   1  13  11   6;
   4   3   2  12   9   5  15  10  11  14   1   7   6   0   8  13
]

# @description The DES' S-box 7
const DES_SBOX7 = UInt8[
   4  11   2  14  15   0   8  13   3  12   9   7   5  10   6   1;
  13   0  11   7   4   9   1  10  14   3   5  12   2  15   8   6;
   1   4  11  13  12   3   7  14  10  15   6   8   0   5   9   2;
   6  11  13   8   1   4  10   7   9   5   0  15  14   2   3  12
]

# @description The DES' S-box 8
const DES_SBOX8 = UInt8[
  13   2   8   4   6  15  11   1  10   9   3  14   5   0  12   7;
   1  15  13   8  10   3   7   4  12   5   6  11   0  14   9   2;
   7  11   4   1   9  12  14   2   0   6  10  13  15   3   5   8;
   2   1  14   7   4  10   8  13  15  12   9   0   3   5   6  11
]

# @description Collection of DES' S-boxes.
const DES_SBOXES = [ DES_SBOX1, DES_SBOX2, DES_SBOX3, DES_SBOX4, DES_SBOX5, DES_SBOX6, DES_SBOX7, DES_SBOX8 ]

# @description Predefined permutation for S-boxes output.
const DES_SBOX_PERMUTATION = UInt8[
  16,     7,    20,    21,
  29,    12,    28,    17,
   1,    15,    23,    26,
   5,    18,    31,    10,
   2,     8,    24,    14,
  32,    27,     3,     9,
  19,    13,    30,     6,
  22,    11,     4,    25
]

# @description Final permutation for the 4-bit preoutput block.
const DES_FINAL_PERMUTATION = UInt8[
  40,     8,   48,    16,    56,   24,    64,   32,
  39,     7,   47,    15,    55,   23,    63,   31,
  38,     6,   46,    14,    54,   22,    62,   30,
  37,     5,   45,    13,    53,   21,    61,   29,
  36,     4,   44,    12,    52,   20,    60,   28,
  35,     3,   43,    11,    51,   19,    59,   27,
  34,     2,   42,    10,    50,   18,    58,   26,
  33,     1,   41,     9,    49,   17,    57,   25
]

# @description Computes the subkey for DES algorithm.
#
# @param key The symmetric key.
#
# @return {Vector{Vector{UInt8}}} The set of subkeys (16 subkeys).
function compute_des_subkeys(key::Vector{UInt8})
  if length(key) ∉ DES_LEGAL_KEY_SIZES
    throw(ArgumentError("Invalid key length. $(length(key) * 8) bits provided, but 64 bits expected."))
  end

  result = Vector{Vector{UInt8}}()
  key    = permutate(unpack_bits(key), DES_PC1)
  cn     = key[1 : 28]
  dn     = key[29 : end]

  for i in 1 : 16
    for j in 1:DES_LEFT_SHIFTS[i]
      cn = circshift(cn, -1)
      dn = circshift(dn, -1)
    end

    subkey = permutate(vcat(cn, dn), DES_PC2)
    push!(result, subkey)
  end

  return result
end

# @description Computes the DES' f(*) function for the given R and subkey.
#
# @param {Vector{UInt8}} r   The R (right) bytes.
# @param {Vector{UInt8}} key The subkey for the DES iteration.
#
# @return {Vector{UInt8}} The round function value.
function compute_des_round(r::Vector{UInt8}, key::Vector{UInt8})
  r = permutate(r, DES_EXPANSION_FUNCTION)
  r = map(xor, r, key)

  b = [ r[1 : 6], r[7 : 12], r[13 : 18], r[19 : 24], r[25 : 30], r[31 : 36], r[37 : 42], r[43 : end] ]

  result = zeros(UInt8, 32)
  index  = 0

  for j in 1 : length(DES_SBOXES)
    m, n  = compute_permutation_position(b[j])
    value = DES_SBOXES[j][m + 1, n + 1]

    result[index + 1] = (value & 8) >> 3
    result[index + 2] = (value & 4) >> 2
    result[index + 3] = (value & 2) >> 1
    result[index + 4] = (value & 1) >> 0
    index += 4
  end

  return permutate(result, DES_SBOX_PERMUTATION)
end

# @description Computes the permutation position in the s-box(n).
#
# @param Vector{UInt8} bn The s-box indices.
#
# @return {Integer, Intenger} The permutation position.
function compute_permutation_position(bn::Vector{UInt8})
  m = (bn[1] << 1) + bn[6]
  n = (bn[2] << 3) + (bn[3] << 2) + (bn[4] << 1) + bn[5]

  return m, n
end


text:
module Data
using Statistics, DataStructures, CSV, DataFrames, Distributed
# Core.println("loading submodule Data on process $(myid()): $(@__MODULE__)")
import ..Utils: log_message, log_progress, PROJECT_ROOT
using ..Utils
include("types.jl")
# this data structure holds all experiments on this process
const experiments = Dict{String, Experiment}()
include("data_handling.jl")
include("create.jl")
include("get.jl")
include("delete.jl")
include("dream_challenge_data.jl")

end  # module Data

ext: jl


text:
function set_boundary(Qbase,cellxmax,cellymax,vecAx,vecAy,bdcon)
    """bdcon[i][j]
        i : 境界番号(x-,x+ ,y-,y+)
        j=1-6 : "bd1_con":"2",
                "bd1_rho":"1.0",
                "bd1_u":"300.0",
                "bd1_v":"0.0",
                "bd1_p":"1.0",
                "bd1_T":"300.0",
    """

    # bd1 = x-
    if Int(bdcon[1][1]) == 0
        for j in 1:cellymax
            Qbase[1,j] = bdcon[1][2]
        end
    elseif Int(bdcon[1][1]) == 1
        for j in 1:cellymax
            Qbase[1,j] = Qbase[2,j]
        end
    else
        println("boundary error")
        println("Please check boundary condition")
        throw(UndefVarError(:x))
    end

    # bd2 = x+
    if Int(bdcon[2][1]) == 0
        for j in 1:cellymax
            Qbase[cellxmax,j] = bdcon[2][2]
        end
    elseif Int(bdcon[2][1]) == 1
        for j in 1:cellymax
            Qbase[cellxmax,j] = Qbase[cellxmax-1,j]
        end
    else
        println("boundary error")
        println("Please check boundary condition")
        throw(UndefVarError(:x))
    end

    # bd3 = y-
    if Int(bdcon[3][1]) == 0
        for i in 1:cellxmax
            Qbase[i,1] = bdcon[3][2]
        end
    elseif Int(bdcon[3][1]) == 1
        for i in 1:cellxmax
            Qbase[i,1] = Qbase[i,2]
        end
    else
        println("boundary error")
        println("Please check boundary condition")
        throw(UndefVarError(:x))
    end

    # bd4 = y+
    if Int(bdcon[4][1]) == 0
        for i in 1:cellxmax
            Qbase[i,cellymax] = bdcon[4][2]
        end
    elseif Int(bdcon[4][1]) == 1
        for i in 1:cellxmax
            Qbase[i,cellymax] = Qbase[i,cellymax-1]
        end
    else
        println("boundary error")
        println("Please check boundary condition")
        throw(UndefVarError(:x))
    end

    return Qbase
end


text:
@with_kw struct ParameterServerOptions
    linger::Cint = Cint(1000)
    master_timeout::Clong = 10000
    worker_timeout::Clong = -1
    scheduler_timeout::Clong = -1
    num_masters::Int32 = Int32(1)
    scheduler_address::NTuple{256,Cchar} = ntuple(i->Cchar('\0'),256)
    master_address::NTuple{256,Cchar} = ntuple(i->Cchar('\0'),256)
    scheduler_pub::UInt16 = 40000
    scheduler_master::UInt16 = 40001
    scheduler_worker::UInt16 = 40002
    master_port::UInt16 = 40001
end

abstract type ParameterServer <: ExecutionPolicy end
paramserver_handle(::ParameterServer) = C_NULL

function initialize!(proxgrad::ProxGradient{<:ParameterServer})
    init_c = @cfunction(init_wrapper, Nothing,
                        (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
    boost = POLO.boosting(proxgrad)
    boosting_c = @cfunction(boost_wrapper, Ptr{Cdouble},
                              (Cint, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
    step = POLO.stepsize(proxgrad)
    step_c = @cfunction(step_wrapper, Cdouble,
                          (Cint, Cint, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
    smooth = POLO.smoothing(proxgrad)
    smoothing_c = @cfunction(smooth_wrapper, Ptr{Cdouble},
                               (Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
    proxim = POLO.prox(proxgrad)
    prox_c = @cfunction(prox_wrapper, Ptr{Cdouble},
                          (Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
    proxgrad.ptr =  ccall(paramserver_handle(proxgrad.execution), Ptr{Nothing},
                          (Ptr{Nothing},
                           Ptr{Nothing}, Any,
                           Ptr{Nothing}, Any,
                           Ptr{Nothing}, Any,
                           Ptr{Nothing}, Any,
                           ParameterServerOptions),
                          init_c,
                          boosting_c, boost,
                          step_c, step,
                          smoothing_c, smooth,
                          prox_c, proxim,
                          proxgrad.execution.psopts)
    return proxgrad
end

include("master.jl")
include("worker.jl")
include("scheduler.jl")


text:
using Tabben: TabularDataset, num_examples, num_classes, num_inputs, num_outputs
using Test

function test_dataset(dataset, n, m, k=1)
    @test num_examples(dataset) == n
    @test num_inputs(dataset) == m
    @test num_outputs(dataset) == k
end

@testset "adult" begin
    test_dataset(TabularDataset("adult", :train), 32_561, 14)
    test_dataset(TabularDataset("adult", :test), 16_281, 14)
end

@testset "amazon" begin
    test_dataset(TabularDataset("amazon", :train), 26_215, 9)
    test_dataset(TabularDataset("amazon", :test), 6_554, 9)
end

@testset "arcene" begin
    test_dataset(TabularDataset("arcene", :train), 100, 10_000)
    test_dataset(TabularDataset("arcene", :test), 100, 10_000)
end

@testset "covertype" begin
    test_dataset(TabularDataset("covertype", :train), 11_340, 54)
    test_dataset(TabularDataset("covertype", :valid), 3_780, 54)
    test_dataset(TabularDataset("covertype", :test), 565_892, 54)
end

@testset "musk" begin
    test_dataset(TabularDataset("musk", :train), 5_548, 166)
    test_dataset(TabularDataset("musk", :test), 1_050, 166)
end

@testset "parkinsons" begin
    test_dataset(TabularDataset("parkinsons", :train), 4_646, 16, 2)
    test_dataset(TabularDataset("parkinsons", :test), 1_229, 16, 2)
end

@testset "poker" begin
    test_dataset(TabularDataset("poker", :train), 25_010, 10)
    test_dataset(TabularDataset("poker", :test), 1_000_000, 10)
end

@testset "rossman" begin
    test_dataset(TabularDataset("rossman", :train), 814_688, 18)
    test_dataset(TabularDataset("rossman", :test), 202_521, 18)
end


text:
mutable struct ContinuousVariable{T <: Real} <: AbstractVariable
    vpop::Vector{T} # current values of the variable over population
    bstv::T         # best value so far
    modl            # probability model
end

function ContinuousVariable(l::T, u::T) where T <: AbstractFloat
    vpop = Vector{T}()
    bstv = initv(T)
    modl = Uniform{T}(l, u)
    ContinuousVariable(vpop, bstv, modl)
end

function ContinuousVariable(l::T, u::T) where T <: Real
    l = float(l)
    u = float(u)
    ContinuousVariable(l, u)
end

function ContinuousVariable(l::T1, u::T2) where {T1 <: Real, T2 <: Real}
    l, u = promote(l, u)
    ContinuousVariable(l, u)
end

fit_model(v::ContinuousVariable, s) = fit_mle(Normal, v.vpop[s])

mutable struct ContinuousBoundedVariable{T <: Real} <: AbstractVariable
    lbnd::T         # lower bound
    ubnd::T         # upper bound
    vpop::Vector{T} # current values of the variable over population
    bstv::T         # best value so far
    modl            # probability model
end

function ContinuousBoundedVariable(lbnd::T, ubnd::T) where T <: AbstractFloat
    vpop = Vector{T}()
    bstv = initv(T)
    modl = Uniform{T}(lbnd, ubnd)
    ContinuousBoundedVariable(lbnd, ubnd, vpop, bstv, modl)
end

function ContinuousBoundedVariable(lbnd::T, ubnd::T) where T <: Real
    lbnd = Float64(lbnd)
    ubnd = Float64(ubnd)
    ContinuousBoundedVariable(lbnd, ubnd)
end

function ContinuousBoundedVariable(lbnd::T1, ubnd::T2) where {T1 <: Real, T2 <: Real}
    lbnd, ubnd = promote(lbnd, ubnd)
    ContinuousBoundedVariable(lbnd, ubnd)
end

fit_model(v::ContinuousBoundedVariable, s) = truncated(fit_mle(Normal, v.vpop[s]), v.lbnd, v.ubnd)


text:
using PETSc2

import MPI
MPI.Init()

PetscInitialize(["-ksp_monitor","-malloc","-malloc_debug","-malloc_dump"]);
#
#   Create a vector and put values in it

comm = MPI.COMM_WORLD
comm_size = MPI.Comm_size(MPI.COMM_WORLD)
comm_rank = MPI.Comm_rank(MPI.COMM_WORLD)

sys_size_local = 3
println("my comm rank = ", comm_rank)

b = PetscVec(comm);
PetscVecSetType(b,"mpi");
PetscVecSetSizes(b,sys_size_local, comm_size*sys_size_local);

x = PetscVec(comm);
PetscVecSetType(x,"mpi");
PetscVecSetSizes(x,sys_size_local, comm_size*sys_size_local);

rhs = zeros(sys_size_local)
for i=1:sys_size_local
  idxm = [(comm_rank)*10 + i]  # index
  PetscVecSetValues(b,idxm,[ convert(Float64,i)],INSERT_VALUES);
  rhs[i] = i
end

PetscVecAssemblyBegin(b);
println("began vector assembly")
PetscVecAssemblyEnd(b);
println("finished vector assembly")
PetscView(b);

println("finished assembling vector")
#
#    Create a matrix

A_julia = [1.0 2.0 3; 4 5 7; 7 8 9]
println("rhs = ", rhs)
println("A_julia = ", A_julia)
x_julia = A_julia\rhs
println("x_julia = ", x_julia)

A = PetscMat(comm);
PetscMatSetType(A,"mpiaij");
PetscMatSetSizes(A,sys_size_local,sys_size_local,comm_size*sys_size_local,comm_size*sys_size_local);
PetscSetUp(A);
for i=1:sys_size_local
  for j = 1:sys_size_local
    idxm = [(comm_rank)*sys_size_local + i]  # row index
    idxn = [(comm_rank)*sys_size_local + j]  # column index
    PetscMatSetValues(A,idxm, idxn, [A_julia[i,j]],INSERT_VALUES);
  end
end

println("finished setting matrix values")
PetscMatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
PetscMatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
println("finished assembling matrix")
PetscView(A);

println("finished viewing matrix")

println("setting up linear solve")
# setup linear solve
ksp = KSP(comm)
KSPSetOperators(ksp, A, A)
KSPSetFromOptions(ksp)
println("performing linear solve")
KSPSolve(ksp, b, x)
println("finished performing linear solve")
PetscView(x)

# copy solution back to Julia
x_copy = zeros(sys_size_local)
idx = collect(0:2)

PetscVecGetValues(x, sys_size_local, idx, x_copy)

println("x_copy = ", x_copy)

err_norm = norm(x_copy - x_julia)
println("err norm = ", err_norm)

PetscDestroy(b)
PetscDestroy(A)
PetscDestroy(ksp)
# remove references to all Petsc variables
# perform solution
#A_julia\rhs
#println("x = ", x)

#PetscFinalize()


text:
#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :tight)
# Set full format, and show decorations
@format full
@testset "rational_power_test" begin
        @test ^(∅, 1//3) == ∅
        @test ^(1 .. 8, 1//3) == Interval(1, 2)
        @test ^(2 .. 8, 1//3) ⊇  Interval(2^(1//3), 2)
        @test ^(1 .. 9, 1//3) ⊇  Interval(1, 9^(1//3))
        @test ^(2 .. 9, 1//3) ⊇  Interval(2^(1//3), 9^(1//3))
        @test ^(-1 .. 8, 1//3) == Interval(0, 2)
        @test ^(-2 .. 8, 1//3) ⊇  Interval(0, 2)
        @test ^(-1 .. 9, 1//3) ⊇  Interval(0, 9^(1//3))
        @test ^(-2 .. 9, 1//3) ⊇  Interval(0, 9^(1//3))
        @test ^(1 .. 8, -1//3) == Interval(0.5, 1)
        @test ^(2 .. 8, -1//3) ⊇  Interval(0.5, 2^(-1//3))
        @test ^(1 .. 9, -1//3) ⊇  Interval(9^(-1//3), 1)
        @test ^(2 .. 9, -1//3) ⊇  Interval(9^(-1//3), 2^(-1//3))
        @test ^(-1 .. 8, -1//3) == Interval(0.5, Inf)
        @test ^(-2 .. 8, -1//3) ⊇  Interval(0.5, Inf)
        @test ^(-1 .. 9, -1//3) ⊇  Interval(9^(-1//3), Inf)
        @test ^(-2 .. 9, -1//3) ⊇  Interval(9^(-1//3), Inf)
        @test ^(-2 .. 4 , 1//2) == Interval(0, 2)
        @test ^(-2 .. 8 , 1//3) == Interval(0, 2)
        @test ^(-8 .. -2 , 1//3) == ∅
        @test ^(-8 .. -2 , 1//2) == ∅
        @test ^(-8 .. -2 , -1//3) == ∅
        @test ^(-8 .. -2 , -1//2) == ∅
        @test ^(∅, 2//3) == ∅
        @test ^(1 .. 8, 2//3) == Interval(1, 4)
        @test ^(2 .. 8, 2//3) ⊇  Interval(2^(2//3), 4)
        @test ^(1 .. 9, 2//3) ⊇  Interval(1, 9^(2//3))
        @test ^(2 .. 9, 2//3) ⊇  Interval(2^(2//3), 9^(2//3))
        @test ^(-1 .. 8, 2//3) == Interval(0, 4)
        @test ^(-2 .. 8, 2//3) ⊇  Interval(0, 4)
        @test ^(-1 .. 9, 2//3) ⊇  Interval(0, 9^(2//3))
        @test ^(-2 .. 9, 2//3) ⊇  Interval(0, 9^(2//3))
        @test ^(1 .. 8, -2//3) == Interval(0.25, 1)
        @test ^(2 .. 8, -2//3) ⊇  Interval(0.25, 2^(-2//3))
        @test ^(1 .. 9, -2//3) ⊇  Interval(9^(-2//3), 1)
        @test ^(2 .. 9, -2//3) ⊇  Interval(9^(-2//3), 2^(-2//3))
        @test ^(-1 .. 8, -2//3) == Interval(0.25, Inf)
        @test ^(-2 .. 8, -2//3) ⊇  Interval(0.25, Inf)
        @test ^(-1 .. 9, -2//3) ⊇  Interval(9^(-2//3), Inf)
        @test ^(-2 .. 9, -2//3) ⊇  Interval(9^(-2//3), Inf)
        @test ^(-2 .. 4 , 3//2) == Interval(0, 8)
        @test ^(-2 .. 8 , 2//3) == Interval(0, 4)
        @test ^(-8 .. -2 , 2//3) == ∅
        @test ^(-8 .. -2 , 3//2) == ∅
        @test ^(-8 .. -2 , -2//3) == ∅
        @test ^(-8 .. -2 , -3//2) == ∅
end


text:
@testset "norm functions" begin

    # First test the un-exported functions which norm(A,p) calls
    # ==========================================================

    @testset "$fnorm(x::Array{$T,$(length(sz))})" for
        fnorm in (
            LinearAlgebra.norm1,
            LinearAlgebra.norm2,
            LinearAlgebra.normInf,
            LinearAlgebra.normMinusInf,
        ),
        T in (Float64, ComplexF64),
        sz in [(3,), (3, 3), (3, 2, 1)]

        x = randn(T, sz)
        # finite differences is unstable if maxabs (minabs) values are not well
        # separated from other values
        if fnorm === LinearAlgebra.normInf
            x[end] = 1000rand(T)
            kwargs = (atol=1e-5, rtol=1e-5)
        elseif fnorm == LinearAlgebra.normMinusInf
            x .*= 1000
            x[end] = rand(T)
            kwargs = (atol=1e-5, rtol=1e-5)
        else
            kwargs = NamedTuple()
        end

        if fnorm === LinearAlgebra.norm2
            @testset "frule" begin
                test_frule(fnorm, x)
            end
        end
        @testset "rrule" begin
            test_rrule(fnorm, x; kwargs...)
            x isa Matrix && @testset "$MT" for MT in (Diagonal, UpperTriangular, LowerTriangular)
                test_rrule(fnorm, MT(x); kwargs...)
            end

            ȳ = rand_tangent(fnorm(x))
            @test unthunk(rrule(fnorm, zero(x))[2](ȳ)[2]) ≈ zero(x)
            @test rrule(fnorm, x)[2](ZeroTangent())[2] isa ZeroTangent
        end
        ndims(x) > 1 && @testset "non-strided" begin
            xp = if x isa Matrix
                view(x, [1,2,3], 1:3)
            elseif x isa Array{T,3}
                PermutedDimsArray(x, (1,2,3))
            end
            @test !(xp isa StridedArray)
            test_rrule(fnorm, xp ⊢ rand(T, size(xp)))
        end
        T == Float64 && ndims(x) == 1 && @testset "Integer input" begin
            x = [1,2,3]
            int_fwd, int_back = rrule(fnorm, x)
            float_fwd, float_back = rrule(fnorm, float(x))
            @test int_fwd ≈ float_fwd
            @test unthunk(int_back(1.0)[2]) ≈ unthunk(float_back(1.0)[2])
        end
    end

    # Next test norm(A, p=2) -- two methods
    # =====================================

    @testset "norm(x::Array{$T,$(length(sz))})" for
        T in (Float64, ComplexF64),
        sz in [(0,), (3,), (3, 3), (3, 2, 1)]

        x = randn(T, sz)

        @testset "frule" begin
            test_frule(norm, x)
            @test frule((ZeroTangent(), ZeroTangent()), norm, x)[2] isa ZeroTangent

            ẋ = rand_tangent(x)
            @test iszero(frule((ZeroTangent(), ẋ), norm, zero(x))[2])
        end
        @testset "rrule" begin
            test_rrule(norm, x)
            x isa Matrix && @testset "$MT" for MT in (Diagonal, UpperTriangular, LowerTriangular)
                test_rrule(norm, MT(x); check_inferred=VERSION>=v"1.5")
            end

            ȳ = rand_tangent(norm(x))
            @test unthunk(rrule(norm, zero(x))[2](ȳ)[2]) ≈ zero(x)
            @test rrule(norm, x)[2](ZeroTangent())[2] isa ZeroTangent
        end
        ndims(x) > 1 && @testset "non-strided" begin
            xp = if x isa Matrix
                view(x, [1,2,3], 1:3)
            elseif x isa Array{T,3}
                PermutedDimsArray(x, (1,2,3))
            end
            @test !(xp isa StridedArray)
            test_frule(norm, xp ⊢ rand(T, size(xp)))
            test_rrule(norm, xp ⊢ rand(T, size(xp)))  # rand_tangent does not work here because eltype(xp)==Int
        end
    end
    @testset "$fnorm(x::Array{$T,$(length(sz))}, $p) with size $sz" for
        fnorm in (norm, LinearAlgebra.normp),
        p in (1.0, 2.0, Inf, -Inf, 2.5),
        T in (Float64, ComplexF64),
        sz in (fnorm === norm ? [(0,), (3,), (3, 3), (3, 2, 1)] : [(3,), (3, 3), (3, 2, 1)])

        x = randn(T, sz)
        # finite differences is unstable if maxabs (minabs) values are not well
        # separated from other values (same as above)
        if p == Inf
            if !isempty(x)
                x[end] = 1000rand(T)
            end
            kwargs = (atol=1e-5, rtol=1e-5)
        elseif p == -Inf
            if !isempty(x)
                x .*= 1000
                x[end] = rand(T)
            end
            kwargs = (atol=1e-5, rtol=1e-5)
        else
            kwargs = NamedTuple()
        end

        test_rrule(fnorm, x, p; kwargs...)
        x isa Matrix && @testset "$MT" for MT in (Diagonal, UpperTriangular, LowerTriangular)
            test_rrule(fnorm, MT(x), p; kwargs..., check_inferred=VERSION>=v"1.5")
        end

        ȳ = rand_tangent(fnorm(x, p))
        @test unthunk(rrule(fnorm, zero(x), p)[2](ȳ)[2]) ≈ zero(x)
        @test rrule(fnorm, x, p)[2](ZeroTangent())[2] isa ZeroTangent
        T == Float64 && sz == (3,) && @testset "Integer input, p=$p" begin
            x = [1,2,3]
            int_fwd, int_back = rrule(fnorm, x, p)
            float_fwd, float_back = rrule(fnorm, float(x), p)
            @test int_fwd ≈ float_fwd
            @test unthunk(int_back(1.0)[2]) ≈ unthunk(float_back(1.0)[2])
        end
    end
    # Extra test for norm(adjoint vector, p)
    @testset "norm($fdual(::Vector{$T}), 2.5)" for
        T in (Float64, ComplexF64),
        fdual in (adjoint, transpose)

        x = fdual(randn(T, 3))
        p = 2.5

        test_rrule(norm, x, p)
        ȳ = rand_tangent(norm(x, p))
        @test unthunk(rrule(norm, x, p)[2](ȳ)[2]) isa typeof(x)
    end

    # Scalar norm(x, p)
    # =================

    @testset "norm(x::$T, p)" for T in (Float64, ComplexF64)
        @testset "p = $p" for p in (-1.0, 2.0, 2.5)
            test_frule(norm, randn(T), p)
            test_rrule(norm, randn(T), p)

            _, back = rrule(norm, randn(T), p)
            @test back(ZeroTangent()) == (NoTangent(), ZeroTangent(), ZeroTangent())
        end
        @testset "p = 0" begin
            p = 0.0
            x = randn(T)
            y = norm(x, p)
            ẋ, ṗ = rand_tangent.((x, p))
            x̄, p̄, ȳ = rand_tangent.((x, p, y))
            y_fwd, ẏ = frule((ZeroTangent(), ẋ, ṗ), norm, x, p)
            @test y_fwd == y
            @test iszero(ẏ)
            y_rev, back = rrule(norm, x, p)
            @test y_rev == y
            @test back(ȳ) == (NoTangent(), zero(x), ZeroTangent())
            @test back(ZeroTangent()) == (NoTangent(), ZeroTangent(), ZeroTangent())
        end
    end
end

# normalise(x, p) and normalise(A, p)
# ===================================

@testset "normalize" begin
    @testset "x::Vector{$T}" for T in (Float64, ComplexF64)
        x = randn(T, 3)
        test_rrule(normalize, x)
        @test rrule(normalize, x)[2](ZeroTangent()) === (NoTangent(), ZeroTangent())
    end
    @testset "x::Vector{$T}, p=$p" for T in (Float64, ComplexF64),
        p in (1.0, 2.0, -Inf, Inf, 2.5) # skip p=0, since FD is unstable
        x = randn(T, 3)
        test_rrule(normalize, x, p)
        @test rrule(normalize, x, p)[2](ZeroTangent()) === (NoTangent(), ZeroTangent(), ZeroTangent())
    end
end


text:
"""
    SetupPeriodicDomain(np::Integer=16)

Set up a periodic domain of size np x np

```
np=16 #domain size is np x np
Γ=SetPeriodicDomain(np)
```
"""
function SetupPeriodicDomain(np::Integer=16)
    γ,Γ=GridOfOnes("PeriodicDomain",1,np)
    Γ["XC"][1]=vec(0.5:1.0:np-0.5)*ones(1,np)
    Γ["XG"][1]=vec(0.0:1.0:np-1.0)*ones(1,np)
    Γ["YC"][1]=ones(np,1)*transpose(vec(0.5:1.0:np-0.5))
    Γ["YG"][1]=ones(np,1)*transpose(vec(0.0:1.0:np-1.0))
    return Γ
end

"""
    SetupRandomFlow(Γ::Dict)

Set up a random flow field over the domain specified by Γ

```
Γ=SetPeriodicDomain(16)
𝑃,ϕ=SetupRandomFlow(Γ)
```
"""
function SetupRandomFlow(Γ::Dict)
  (_,ϕ,_,_)=demo2(Γ);

  (u,v)=gradient(ϕ,Γ)
  u=u./Γ["DXC"]#normalization to grid units
  v=v./Γ["DYC"]

  (u,v)=exchange(u,v,1)
  u0=-v; u1=-v;
  v0=u; v1=u;

  𝑃 = Dict( "u0" => u0, "u1" => u1, "v0" => v0, "v1" => v1,
            "t0" => 0.0, "t1" => 400.0, "dt" => 0.1)
  𝑃=merge(𝑃,Γ)#add grid variables

  return 𝑃,ϕ
end

max_stars_repo_name: gaelforget/MarineEcosystemNotebooks

max_forks_repo_name: gaelforget/MarineEcosystemNotebooks


text:
using BenchmarkTools

# Conclusion
# `BLAS.gemv!` is a little bit faster than native operators.

# --------Native---------
# BenchmarkTools.Trial:
#   memory estimate:  156.41 KiB
#   allocs estimate:  4
#   --------------
#   minimum time:     51.991 ms (0.00% GC)
#   median time:      53.109 ms (0.00% GC)
#   mean time:        53.112 ms (0.00% GC)
#   maximum time:     54.377 ms (0.00% GC)
#   --------------
#   samples:          95
#   evals/sample:     1
# --------BLAS-----------
# BenchmarkTools.Trial:
#   memory estimate:  0 bytes
#   allocs estimate:  0
#   --------------
#   minimum time:     52.967 ms (0.00% GC)
#   median time:      56.853 ms (0.00% GC)
#   mean time:        56.331 ms (0.00% GC)
#   maximum time:     59.947 ms (0.00% GC)
#   --------------
#   samples:          89
#   evals/sample:     1

begin

weight = rand(10000, 10000)
bias = rand(10000)
x = rand(10000)

# whether gemv is faster

# Native Julia
native = @benchmark weight * x + bias

y = copy(bias)
blas = @benchmark BLAS.gemv!('N', 1.0, weight, x, 1.0, y)

println("--------Native---------")
display(native)
println("\n--------BLAS-----------")
display(blas)

end

# --------Native---------
# BenchmarkTools.Trial:
#   memory estimate:  2.98 GiB
#   allocs estimate:  6
#   --------------
#   minimum time:     1.692 s (3.60% GC)
#   median time:      1.768 s (7.84% GC)
#   mean time:        1.761 s (7.52% GC)
#   maximum time:     1.824 s (10.85% GC)
#   --------------
#   samples:          3
#   evals/sample:     1
# --------BLAS-----------
# BenchmarkTools.Trial:
#   memory estimate:  156.33 KiB
#   allocs estimate:  2
#   --------------
#   minimum time:     100.829 ms (0.00% GC)
#   median time:      109.566 ms (0.00% GC)
#   mean time:        107.981 ms (0.00% GC)
#   maximum time:     112.653 ms (0.00% GC)
#   --------------
#   samples:          47
#   evals/sample:     1

begin

weight = rand(Complex128, 10000, 10000)
dz = rand(Complex128, 10000)
native = @benchmark gradOutput = conj(transpose(weight)) * dz
blas = @benchmark gradOutput = BLAS.gemv('C', weight, dz)

println("\n--------Native---------")
display(native)
println("\n--------BLAS-----------")
display(blas)

end

# backward
begin

# --------Native---------
# BenchmarkTools.Trial:
#   memory estimate:  1.49 GiB
#   allocs estimate:  2
#   --------------
#   minimum time:     584.670 ms (1.69% GC)
#   median time:      627.907 ms (7.79% GC)
#   mean time:        628.836 ms (8.26% GC)
#   maximum time:     672.612 ms (14.50% GC)
#   --------------
#   samples:          8
#   evals/sample:     1
# --------BLAS-----------
# BenchmarkTools.Trial:
#   memory estimate:  1.49 GiB
#   allocs estimate:  6
#   --------------
#   minimum time:     418.717 ms (2.36% GC)
#   median time:      479.637 ms (13.84% GC)
#   mean time:        488.342 ms (14.29% GC)
#   maximum time:     561.557 ms (11.78% GC)
#   --------------
#   samples:          11
#   evals/sample:     1%

input = rand(Complex128, 10000)
dz = rand(Complex128, 10000)

function native_backward(input, dz)
    return dz * conj(transpose(input))
end

function blas_backward(input::Vector, dz::Vector)
    input = reshape(input, size(input, 1), 1)
    dz = reshape(dz, size(dz, 1), 1)
    BLAS.gemm('N', 'C', dz, input)
end

native = @benchmark output = native_backward($input, $dz)
blas = @benchmark ouput = blas_backward($input, $dz)

println("\n--------Native---------")
display(native)
println("\n--------BLAS-----------")
display(blas)

end # backward


text:
"""
Stationary vector with iterative methods
"""

"""
stguess(Q::MatT, ::Type{Tv} = Float64)::Vector{Tv}

Get a vector which is guessed as the stationary vector of CTMC.
This is used as the initial vector for the iterative methods.
"""

function stguess(Q::MatT, ::Type{Tv} = Float64)::Vector{Tv} where {Tv,MatT}
    m, n = size(Q)
    @assert m == n
    result = Vector{Tv}(undef, n)
    for (i,x) in enumerate(spdiag(Q))
        result[i] = 1/x
    end
    result ./= sum(result)
end

"""
stgs(Q::SparseCSC{Tv,Ti}, x0::Vector{Tv}=stguess(Q,Tv); maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6))

Get a stationary vector of CTMC.

Parameters:
- Q: CTMC Kernal
- x0: Initial vector for iteration
- maxiter: The maximum number of iteration. The algorithm stops when the number of iteration becomes maxiter.
- steps: The number of steps to check the convergence
- rtol: the tolerance error. When the relative errors of two successive vectors with steps attains rtol, the algorithm stops.
Return value:
A tuple of
- x: stationary vector
- conv: A boolean whether the algorithm converges or not
- iter: The number of iterations
- rerror: The relative error when the algorithm stops
"""

function stgs(Q::SparseMatrixCSC{Tv,Ti}; x0::Vector{Tv}=stguess(Q,Tv),
        maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6)) where {Tv,Ti}
    stgs(SparseCSC(Q), x0=x0, maxiter=maxiter, steps=steps, rtol=rtol)
end

function stgs(Q::SparseCSC{Tv,Ti}; x0::Vector{Tv}=stguess(Q,Tv),
        maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6)) where {Tv,Ti}
    m, n = size(Q)
    @assert m == n
    b = zeros(Tv, n)
    x = copy(x0)
    iter = 0
    conv = false
    rerror::Tv = 0
    prevx = similar(x)
    while true
        prevx .= x
        for i in 1:steps
            gsstep!(x, Q, b)
            x ./= sum(x)
        end
        # rerror = maximum(abs.((x - prevx) ./ x))
        rerror = maximum(abs.(x - prevx)) / maximum(x)
        iter += steps
        if rerror < rtol
            conv = true
            break
        end
        if iter >= maxiter
            break
        end
    end
    return x, conv, iter, rerror
end

"""
stpower(P::AbstractMatrix{Tv}, x0::Vector{Tv}=stguess(Q,Tv); maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6))

Get a stationary vector of DTMC with power method.

Parameters:
- P: The transition probability matrix for DTMC
- x0: Initial vector for iteration
- maxiter: The maximum number of iteration. The algorithm stops when the number of iteration becomes maxiter.
- steps: The number of steps to check the convergence
- rtol: the tolerance error. When the relative errors of two successive vectors with steps attains rtol, the algorithm stops.
Return value:
A tuple of
- x: stationary vector
- conv: A boolean whether the algorithm converges or not
- iter: The number of iterations
- rerror: The relative error when the algorithm stops
"""

function stpower(P::AbstractMatrix{Tv}; x0::Vector{Tv}=stguess(P,Tv),
    maxiter=5000, steps=20, rtol::Tv=Tv(1.0e-6)) where {Tv}
    m, n = size(P)
    @assert m == n
    Pdash = P'
    x = copy(x0)
    iter = 0
    conv = false
    rerror::Tv = 0
    prevx = similar(x)
    while true
        prevx .= x
        for i in 1:steps
            x = Pdash * x
            x ./= sum(x)
        end
        # rerror = maximum(abs.((x - prevx) ./ x))
        rerror = maximum(abs.(x - prevx)) / maximum(x)
        iter += steps
        if rerror < rtol
            conv = true
            break
        end
        if iter >= maxiter
            break
        end
    end
    return x, conv, iter, rerror
end


text:
using ADSeismic
using ADCME
using PyPlot
using MAT
using DelimitedFiles
matplotlib.use("Agg")
close("all")

################### Generate synthetic data #####################
model_name = "../../../examples/nn_fwi/models/marmousi2-model-true.mat"
# model_name = "models/BP-model-true.mat"

## load model setting
PropagatorKernel = 1
params = load_params(model_name, vp_ref=1000, PropagatorKernel=PropagatorKernel)
src = load_acoustic_source(model_name)
rcv = load_acoustic_receiver(model_name)
vp = constant(matread(model_name)["vp"])
## For debug, only run the first source
# src = [src[1]]
# rcv = [rcv[1]]

## assemble acoustic propagator model
model = x->AcousticPropagatorSolver(params, x, vp)

## simulated wavefield
Rs_ = compute_forward_GPU(model, src, rcv)
sess = Session(); init(sess)

# ts = zeros(11)
# for i = 1:11
#     d = @timed Rs = run(sess, Rs_)
#     ts[i] = d[2]
#     @info "Test $i, time = $(ts[i])"
# end
# @info mean(ts[2:end])

## visualize_wavefield if needed
ii = 1 # source number
u = run(sess, model(src[ii]).u)
p = visualize_wavefield(u, params)
saveanim(p,"test$PropagatorKernel.gif")


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

using DensityInterface
using Test

using LinearAlgebra, InverseFunctions

struct MyDensity end
@inline DensityInterface.DensityKind(::MyDensity) = IsDensity()
DensityInterface.logdensityof(::MyDensity, x::Any) = -norm(x)^2

struct MyMeasure end
@inline DensityInterface.DensityKind(::MyMeasure) = HasDensity()
DensityInterface.logdensityof(::MyMeasure, x::Any) = -norm(x)^2

@testset "interface" begin
    @test inverse(logdensityof) == logfuncdensity
    @test inverse(logfuncdensity) == logdensityof
    @test inverse(densityof) == funcdensity
    @test inverse(funcdensity) == densityof

    @test @inferred(DensityKind("foo")) === NoDensity()
    @test_throws ArgumentError logdensityof("foo")
    @test_throws ArgumentError densityof("foo")

    for object1 in (MyDensity(), MyMeasure())
        x = [1, 2, 3]

        DensityInterface.test_density_interface(object1, x, -norm(x)^2)

        object2 = logfuncdensity(x -> -norm(x)^2)
        @test DensityKind(object2) === IsDensity()
        DensityInterface.test_density_interface(object2, x, -norm(x)^2)

        object3 = funcdensity(x -> exp(-norm(x)^2))
        @test DensityKind(object3) === IsDensity()
        DensityInterface.test_density_interface(object3, x, -norm(x)^2)
    end
end


text:
@osx_only immutable Cpasswd
   pw_name::Cstring
   pw_passwd::Cstring
   pw_uid::Cint
   pw_gid::Cint
   pw_change::Cint
   pw_class::Cstring
   pw_gecos::Cstring
   pw_dir::Cstring
   pw_shell::Cstring
   pw_expire::Cint
   pw_fields::Cint
end

@linux_only immutable Cpasswd
   pw_name::Cstring
   pw_passwd::Cstring
   pw_uid::Cint
   pw_gid::Cint
   pw_gecos::Cstring
   pw_dir::Cstring
   pw_shell::Cstring
end

@unix_only immutable Cgroup
    gr_name::Cstring
    gr_passwd::Cstring
    gr_gid::Cint
end

immutable User
    name::ASCIIString
    uid::UInt64
    gid::UInt64
    dir::ASCIIString
    shell::ASCIIString

    function User(passwd::Ptr{Cpasswd})
        ps = unsafe_load(passwd)

        new(
            pointer_to_string(ps.pw_name),
            UInt64(ps.pw_uid),
            UInt64(ps.pw_gid),
            pointer_to_string(ps.pw_dir),
            pointer_to_string(ps.pw_shell)
        )
    end
end

function Base.show(io::IO, user::User)
    print(io, "$(user.uid) ($(user.name))")
end

function User(name::ASCIIString)
    ps = ccall((:getpwnam, "libc"), Ptr{Cpasswd}, (Ptr{UInt8},), name)
    User(ps)
end

function User(uid::UInt64)
    ps = ccall((:getpwuid, "libc"), Ptr{Cpasswd}, (UInt64,), uid)
    User(ps)
end

function User()
    uid = ccall((:geteuid, "libc"), Cint, ())
    User(UInt64(uid))
end

immutable Group
    name::ASCIIString
    gid::UInt64

    function Group(group::Ptr{Cgroup})
        gr = unsafe_load(group)

        new(
            pointer_to_string(gr.gr_name),
            UInt64(gr.gr_gid)
        )
    end
end

function Base.show(io::IO, group::Group)
    print(io, "$(group.gid) ($(group.name))")
end

function Group(name::ASCIIString)
    ps = ccall((:getgrnam, "libc"), Ptr{Cgroup}, (Ptr{UInt8},), name)
    Group(ps)
end

function Group(gid::UInt64)
    gr = ccall((:getgrgid, "libc"), Ptr{Cgroup}, (UInt64,), gid)
    Group(gr)
end

function Group()
    gid = ccall((:getegid, "libc"), Cint, ())
    Group(UInt64(gid))
end


text:
"""
    GP{Tm<:MeanFunction, Tk<:Kernel}

A Gaussian Process (GP) with known `mean` and `kernel`. See e.g. [1] for an introduction.

# Zero Mean
If only one argument is provided, assume the mean to be zero everywhere:

```jldoctest
julia> f = GP(Matern32Kernel());

julia> x = randn(5);

julia> mean(f(x)) == zeros(5)
true

julia> cov(f(x)) == kernelmatrix(Matern32Kernel(), x)
true
```

### Constant Mean

If a `Real` is provided as the first argument, assume the mean function is
constant with that value.

```jldoctest
julia> f = GP(5.0, Matern32Kernel());

julia> x = randn(5);

julia> mean(f(x)) == 5.0 .* ones(5)
true

julia> cov(f(x)) == kernelmatrix(Matern32Kernel(), x)
true
```

### Custom Mean

Provide an arbitrary function to compute the mean:

```jldoctest
julia> f = GP(x -> sin(x) + cos(x / 2), Matern32Kernel());

julia> x = randn(5);

julia> mean(f(x)) == sin.(x) .+ cos.(x ./ 2)
true

julia> cov(f(x)) == kernelmatrix(Matern32Kernel(), x)
true
```

[1] - C. E. Rasmussen and C. Williams. "Gaussian processes for machine learning".
MIT Press. 2006.
"""
struct GP{Tm<:MeanFunction,Tk<:Kernel} <: AbstractGP
    mean::Tm
    kernel::Tk
end

GP(mean, kernel::Kernel) = GP(CustomMean(mean), kernel)
GP(mean::Real, kernel::Kernel) = GP(ConstMean(mean), kernel)
GP(kernel::Kernel) = GP(ZeroMean(), kernel)

# AbstractGP interface implementation.

Statistics.mean(f::GP, x::AbstractVector) = _map_meanfunction(f.mean, x)

Statistics.cov(f::GP, x::AbstractVector) = kernelmatrix(f.kernel, x)

Statistics.var(f::GP, x::AbstractVector) = kernelmatrix_diag(f.kernel, x)

Statistics.cov(f::GP, x::AbstractVector, x′::AbstractVector) = kernelmatrix(f.kernel, x, x′)

max_stars_repo_name: JuliaGaussianProcesses/AbstractGP

max_issues_repo_name: JuliaGaussianProcesses/AbstractGP

max_forks_repo_name: JuliaGaussianProcesses/AbstractGP