text:
"Construction of the Kravchuk matrix using three-term recurrence relation"

export kmatrix

begin
    k₁(x, k, j) = ((2*x)/√((2 + 2*j - k)*(k - 1)))
    k₂(k, j) = √(((3 + 2*j - k)*(k - 2))/((2 + 2*j - k)*(k - 1)))
end

function kmatrix(j)
    vk = zeros(Float64, (Integer(2*j + 1), Integer(2*j + 1)))
    Threads.@threads for x in -j:j
        vk[Integer(x + j + 1), 1] = 2.0^(-j) * √binomial(Integer(2*j), Integer(j + x))
        vk[Integer(x + j + 1), 2] = ((2*x)/√(2*j)) * vk[Integer(x + j + 1), 1]
        for k in 2:Integer(2*j+1)
            if k == 2
                vk[Integer(x + j + 1), k] = k₁(x, k, j) * vk[Integer(x + j + 1), k - 1]
            else
                vk[Integer(x + j + 1), k] = k₁(x, k, j) * vk[Integer(x + j + 1), k - 1] - k₂(k, j) * vk[Integer(x + j + 1), k - 2]
            end
        end
    end
    return vk
end

max_stars_repo_name: rurz/UnitRots

max_issues_repo_name: rurz/UnitRots

max_forks_repo_name: rurz/UnitRots


text:
# Various algorithms for computing quantile

function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real, lx::T, rx::T) where {T<:Real}
    rx < lx && throw(ArgumentError("empty bracketing interval [$lx, $rx]"))

    # In some special cases, e.g. #1501, rx == lx`
    # If the distribution is degenerate the check below can fail, hence we skip it
    if rx == lx
        # Returns `lx` of the same type as `(lx + rx) / 2`
        # For specific types such as `Float64` it is more performant than `oftype((lx + rx) / 2, lx)`
        return middle(lx)
    end

    # base tolerance on types to support e.g. `Float32` (avoids an infinite loop)
    # ≈ 3.7e-11 for Float64
    # ≈ 2.4e-5 for Float32
    tol = cbrt(eps(float(T)))^2
    # find quantile using bisect algorithm
    cl = cdf(d, lx)
    cr = cdf(d, rx)
    cl <= p <= cr ||
        throw(ArgumentError("[$lx, $rx] is not a valid bracketing interval for `quantile(d, $p)`"))
    while rx - lx > tol
        m = (lx + rx)/2
        c = cdf(d, m)
        if p > c
            cl = c
            lx = m
        else
            cr = c
            rx = m
        end
    end
    return (lx + rx)/2
end

function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real, lx::Real, rx::Real)
    return quantile_bisect(d, p, promote(lx, rx)...)
end

quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) =
    quantile_bisect(d, p, minimum(d), maximum(d))

# if starting at mode, Newton is convergent for any unimodal continuous distribution, see:
#   Göknur Giner, Gordon K. Smyth (2014)
#   A Monotonically Convergent Newton Iteration for the Quantiles of any Unimodal
#   Distribution, with Application to the Inverse Gaussian Distribution
#   http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf

function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
    x = xs + (p - cdf(d, xs)) / pdf(d, xs)
    T = typeof(x)
    if 0 < p < 1
        x0 = T(xs)
        while abs(x-x0) > max(abs(x),abs(x0)) * tol
            x0 = x
            x = x0 + (p - cdf(d, x0)) / pdf(d, x0)
        end
        return x
    elseif p == 0
        return T(minimum(d))
    elseif p == 1
        return T(maximum(d))
    else
        return T(NaN)
    end
end

function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
    x = xs + (ccdf(d, xs)-p) / pdf(d, xs)
    T = typeof(x)
    if 0 < p < 1
        x0 = T(xs)
        while abs(x-x0) > max(abs(x),abs(x0)) * tol
            x0 = x
            x = x0 + (ccdf(d, x0)-p) / pdf(d, x0)
        end
        return x
    elseif p == 1
        return T(minimum(d))
    elseif p == 0
        return T(maximum(d))
    else
        return T(NaN)
    end
end

function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
    T = typeof(lp - logpdf(d,xs))
    if -Inf < lp < 0
        x0 = T(xs)
        if lp < logcdf(d,x0)
            x = x0 - exp(lp - logpdf(d,x0) + logexpm1(max(logcdf(d,x0)-lp,0)))
            while abs(x-x0) >= max(abs(x),abs(x0)) * tol
                x0 = x
                x = x0 - exp(lp - logpdf(d,x0) + logexpm1(max(logcdf(d,x0)-lp,0)))
            end
        else
            x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0)))
            while abs(x-x0) >= max(abs(x),abs(x0))*tol
                x0 = x
                x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0)))
            end
        end
        return x
    elseif lp == -Inf
        return T(minimum(d))
    elseif lp == 0
        return T(maximum(d))
    else
        return T(NaN)
    end
end

function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
    T = typeof(lp - logpdf(d,xs))
    if -Inf < lp < 0
        x0 = T(xs)
        if lp < logccdf(d,x0)
            x = x0 + exp(lp - logpdf(d,x0) + logexpm1(max(logccdf(d,x0)-lp,0)))
            while abs(x-x0) >= max(abs(x),abs(x0)) * tol
                x0 = x
                x = x0 + exp(lp - logpdf(d,x0) + logexpm1(max(logccdf(d,x0)-lp,0)))
            end
        else
            x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0)))
            while abs(x-x0) >= max(abs(x),abs(x0)) * tol
                x0 = x
                x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0)))
            end
        end
        return x
    elseif lp == -Inf
        return T(maximum(d))
    elseif lp == 0
        return T(minimum(d))
    else
        return T(NaN)
    end
end

# A macro: specify that the quantile (and friends) of distribution D
# is computed using the newton method
macro quantile_newton(D)
    esc(quote
        quantile(d::$D, p::Real) = quantile_newton(d,p)
        cquantile(d::$D, p::Real) = cquantile_newton(d,p)
        invlogcdf(d::$D, lp::Real) = invlogcdf_newton(d,lp)
        invlogccdf(d::$D, lp::Real) = invlogccdf_newton(d,lp)
    end)
end

ext: jl


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("misc.jl")
coefs = Float64.([1 1 1])

thumbstack_stencil = @def_stencil_expression(
            @sum(i,-2,2,(
                @sum(j,-2,2, c[max(abs(i),abs(j))]*D[x+i,y+j,z])
            )) + @sum(i,1,2, c[i]*(D[x,y,z+i]+D[x,y,z-i]))
)

st_inst = NewStencilInstance(thumbstack_stencil, 2, coefs)

radius = 2
nx = 8
ny = 8
nz = 8
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

t_steps = 5
## Actual time steps
data = CreateData(dx,dy,dz,radius)
dout = ApplyStencil(st_inst, data, t_steps)
#heatmap(dout1[:,:,Int(round(dz/2))])
o = Array(dout)


text:
#########################################################################################
# Linear fit with GLM
#########################################################################################
import GLM

function linear_regression_fit_glm(x,y)
    X = hcat(ones(length(x)), x)
    out = GLM.lm(X, y)
    o, s = GLM.coef(out)
    s05, s95 = GLM.confint(out)[2, :]
    return s, s05, s95
end

#########################################################################################
# Linear fit with LsqFit
#########################################################################################
import LsqFit

linearregression(x, p) = @. p[1] + p[2]*x #+ p[3]*log(-x)

function linear_regression_fit_lsqfit(x,y, p0 = [0.0, 0.0])
    linearregression(x, p) = @. p[1] + p[2]*x
    fit = LsqFit.curve_fit(linearregression, x, y, p0)
    s = LsqFit.coef(fit)[2]
    s05, s95 = LsqFit.confidence_interval(fit, 0.05)[2]
    return s, s05, s95
end

#########################################################################################
# Logarithmically corrected fit with LsqFit
#########################################################################################

function logarithmic_corrected_fit_lsqfit(x, y, p0 = [0.0, 0.0, 0.0])
    corrected_fit(x, p) = @. p[1] + p[2]*x + p[3]*log(-x)
    i = findlast(<(0), x)
    if (!isnothing(i) && i < length(x)÷2)
        @info "More than half of the data have log(ε) > 0 for corrected version. Returning standard linear..."
        return linear_regression_fit_glm(x,y)
    elseif !isnothing(i)
        x2 = x[1:i]
        y2 = y[1:i]
    else # last case: all x are negative
        x2 = x
        y2 = y
    end
    if !all(x2 .< 0)
        @info "There still exist positive log(ε) for corrected version. Using standard linear..."
        return linear_regression_fit_glm(x,y)
    end
    try
        fit = LsqFit.curve_fit(corrected_fit, x2, y2, p0)
        s = LsqFit.coef(fit)[2]
        s05, s95 = LsqFit.confidence_interval(fit, 0.05)[2]
        return s, s05, s95
    catch err
        @info "Corrected version errored with $(err). Using standard linear..."
        return linear_regression_fit_glm(x,y)
    end
end

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
#train the neural net for dataset (datax,datay)
#requires that the function nnDCost be defined representing the derivative of the cost function
function nnTrain(nn::neuralNet,datax,datay,epochs::Int64,learnRate::Float64;batchSize=size(datax,2),batchLength=epochs,verbose=false)

    epCnt = 0;
    bCnt  = 1;

    while epCnt < epochs
        #get indices of the data to be used
        (batchSize == size(datax,2)) ? (dIdx=1:batchSize) : (dIdx=sample(1:size(datax,2),batchSize,replace=false));

        #loop over epochs
        for ep=1:min(batchLength,epochs-epCnt)
            #verbose && @printf("starting epoch %3d...\n",ep);

            #these vectors will store the running sums
            c  = 0.0;
            dw = zeros(length(nn.dCdw));
            db = zeros(length(nn.dCdb));

            #loop over training data
            for d=dIdx
                #pick out the dataset
                dx = datax[:,d];
                dy = datay[:,d];

                #forward run
                y = nnForward(nn,dx);

                #add costs to running sum
                c += nnCost(dy,y);

                #back propagate
                nnBackPropagate(nn,dy);

                #add computed gradients to running sum
                dw += nn.dCdw;
                db += nn.dCdb;
            end

            #shift by mean * learning rate
            nn.w -= learnRate*(dw ./ batchSize);
            nn.b -= learnRate*(db ./ batchSize);

            #verbose && println("finished epoch $ep.");
            verbose && @printf("Epoch %5d (iteration %5d of batch %5d): Cost = %10.8f\n",ep+epCnt,ep,bCnt,c/batchSize);
        end
        bCnt  += 1;
        epCnt += batchLength;
    end

end

max_stars_repo_name: krometis/SimpleNeuralNet

max_issues_repo_name: krometis/SimpleNeuralNet

max_forks_repo_name: krometis/SimpleNeuralNet


text:
export  CreateTable,
        addtablecolumn!,
        setunique!,
        setprimary!,
        createtable

"""
```
CreateTable( tablename::AbstractString )
```
This mutable struct represents a `CREATE TABLE` SQL statement where the name of the table is `tablename`.
"""
mutable struct CreateTable <: SQLStatement
    name::String
    columns::Vector{Var}
    types::Vector{String}
    flags::Vector{Vector{String}}
    unique::Vector{Var}
    primary::Vector{Var}

    CreateTable( tablename::AbstractString ) = new( tablename, Var[], String[], String[], Var[], Var[] )
end

function Base.show( io::IO, ct::CreateTable )
    isempty(ct.columns) && return
    print( io, "CREATE TABLE \"", ct.name, "\" (" )

    for ii in eachindex(ct.columns)
        ii > 1 && print( io, "," )
        print( io, "\n    ", ct.columns[ii], " ", ct.types[ii] )
        isempty(ct.flags[ii]) || print( io, " ", join( ct.flags[ii], " " ) )
    end

    isempty(ct.unique) || print( io, ",\n    UNIQUE (", join( ct.unique, ", " ), ")" )
    isempty(ct.primary) || print( io, ",\n    PRIMARY KEY (", join( ct.primary, ", " ), ")" )
    print( io, "\n)" )
end

"""
```
addtablecolumn!(
    sc::Var,
    valtype::AbstractString,
    flags::Vector{T}=String[] )
addtablecolumn!(
    sc::Var,
    valtype::AbstractString,
    flags::Vector{T}=String[] )
```
This function permits the user to add a command to add a column with name `sc` of type `valtype` to a `CreateTable` struct. Any additional flags that this column has, such as `NOT NULL`, etc. can be set in `flags`.

The output of this function is a function that takes a `CreateTable` as argument and returns the updated argument, such that Julia's function chaining syntax can be used.

Example:
```
sq |> addtablecolumn!( "var1", "int" ) |> addtablecolumn!( "var2", "real" )
```
adds two columns to the `CreateTable` struct: a column `var1` of type `int`, and a column `var2` of type `real`.
"""
addtablecolumn!( sc::Var, valtype::AbstractString, flags::Vector{T}=String[] ) where T <: AbstractString = ct::CreateTable -> begin
    push!( ct.columns, sc )
    push!( ct.types, valtype )
    push!( ct.flags, flags )
    ct
end

addtablecolumn!( name::AbstractString, valtype::AbstractString, flags::Vector{T}=String[] ) where T <: AbstractString =
    addtablecolumn!( name |> Var, valtype, flags )

allunique( ct::CreateTable ) = length(ct.columns) == length(unique(ct.columns))

function makeunique!( ct::CreateTable )
    makeunique!( ct.columns )
    ct
end

"""
```
setunique!( scs::Vector{Var} )
setunique!( scs::Vector{T} ) where T <: AbstractString
setunique!( scs::Union{Var, AbstractString}... )
```
This function permits the user to add a command that sets the unique clause to enforce uniqueness on the combination of the variables in `scs` of a `CreateTable` struct.

The output of this function is a function that takes a `CreateTable` as argument and returns the updated argument, such that Julia's function chaining syntax can be used.

Example:
```
sq |> setunique!( "var1", "var2" )
```
sets the uniqueness constraint that all combinations of values of `var1` and `var2` in the table must be unique.
"""
setunique!( scs::Vector{Var} ) = ct::CreateTable -> begin
    ct.unique = deepcopy(scs)
    ct
end

setunique!( scs::Vector{T} ) where T <: AbstractString = Var.(scs) |> setunique!
setunique!( scs::SVar... ) = Var.( scs |> collect ) |> setunique!

"""
```
setprimary!( scs::Vector{Var} )
setprimary!( scs::Vector{T} ) where T <: AbstractString
setprimary!( scs::Union{Var, AbstractString}... )
```
This function permits the user to add a command that sets the primary key clause to the variables in `scs` in a `CreateTable` struct. A primary key cluase is like the unique clause, with the added constraint that none of the variables can have a null value.

The output of this function is a function that takes a `CreateTable` as argument and returns the updated argument, such that Julia's function chaining syntax can be used.

Example:
```
sq |> setprimary!( "var1", "var2" )
```
sets the primary key clause to the combination of `var1` and `var2`.
"""
setprimary!( scs::Vector{Var} ) = ct::CreateTable -> begin
    ct.primary = deepcopy(scs)
    ct
end

setprimary!( colnames::Vector{T} ) where T <: AbstractString = Var.(colnames) |> setprimary!
setprimary!( scs::SVar... ) = Var.( scs |> collect ) |> setprimary!

function sanitise!( ct::CreateTable )
    ct |> makeunique!
    filter!( var -> var ∈ ct.columns, ct.unique ) |> unique!
    filter!( var -> var ∈ ct.columns, ct.primary ) |> unique!
    ct
end

"""
```
createtable( conn::Conn )
```
This function permits the user to execute the SQL equivalent of a `CreateTable` struct on the database that `conn` connects to.

The output of this function is a function that takes a `CreateTable` as argument and returns the result of the executed SQL statement.

Example:
```
sq |> createtable(conn)
```
"""
createtable( conn::Conn ) = ct::CreateTable -> LibPQ.execute( conn, string( ct |> sanitise!, ";" ) )


text:
# This code presents an example to illustrate how SmolyakApprox can be used

using SmolyakApprox

function test_smolyak_derivative()

  d  = 5  # Set the number of dimensions
  mu = 3  # Set the level

  grid, multi_ind = smolyak_grid(chebyshev_gauss_lobatto,d,mu)  # Construct the Smolyak grid and the multi index

  # An arbitrary test function

  function test(grid)

    y_value = (grid[:,1].+1).^0.1.*exp.(grid[:,2]).*log.(grid[:,3].+2).^0.2.*(grid[:,4].+2).^0.8.*(grid[:,5].+7).^0.1

    return y_value

  end

  y = test(grid)  # Evaluate the test function on the Smolyak grid

  point = [0.75, 0.45, 0.82, -0.15, -0.95]  # Choose a point to evaluate the approximated function

  # One way of computing the weights and evaluating the approximated function

  weights = smolyak_weights(y,grid,multi_ind)        # Compute the Smolyak weights
  y_hat = smolyak_evaluate(weights,point,multi_ind)  # Evaluate the approximated function

  #= A second way of computing the weights and evaluating the approximated function that
  computes the interpolation matrix just once. =#

  interp_mat = smolyak_inverse_interpolation_matrix(grid,multi_ind)  # Compute the interpolation matrix
  w = smolyak_weights(y,interp_mat)             # Compute the Smolyak weights
  y_hatt = smolyak_evaluate(w,point,multi_ind)  # Evaluate the approximated function

  # Evaluate the exact function at point

  y_actual = test(point')

  derivatives_1 = smolyak_derivative_finite_difference(weights,point,multi_ind)
  derivatives_2 = smolyak_gradient(weights,point,multi_ind)
  derivatives_3 = smolyak_derivative(weights,point,multi_ind,1)
  derivatives_4 = smolyak_derivative(weights,point,multi_ind,3)

  grid_full, multi_ind_full = smolyak_grid_full(chebyshev_gauss_lobatto,d,mu)  # Construct the Smolyak grid and the multi index
  y_full = test(grid_full)                                                     # Evaluate the test function on the Smolyak grid
  weights_full = smolyak_weights_full(y_full,grid_full,multi_ind_full)         # Compute the Smolyak weights

  derivatives_3_full = smolyak_derivative_full(weights_full,point,multi_ind_full,1)         # Evaluate the approximated function
  derivatives_4_full = smolyak_gradient_full(weights_full,point,multi_ind_full) # Evaluate the approximated function

  return derivatives_1, derivatives_2, derivatives_3, derivatives_4, derivatives_3_full, derivatives_4_full

end

test_smolyak_derivative()


text:
using REPL
using REPL.LineEdit
using REPL.REPLCompletions
using JuliaInterpreter: moduleof, locals, eval_code
import ..Atom: @msg

const normal_prefix = @static Sys.iswindows() ? "\e[33m" : "\e[38;5;166m"
const compiled_prefix = "\e[96m"

function debugprompt()
  try
    panel = REPL.LineEdit.Prompt("debug> ";
              prompt_prefix = isCompileMode() ? compiled_prefix : normal_prefix,
              prompt_suffix = Base.text_colors[:normal],
              complete = JunoDebuggerRPELCompletionProvider(),
              on_enter = REPL.return_callback)

    panel.hist = REPL.REPLHistoryProvider(Dict{Symbol,Any}(:junodebug => panel))
    REPL.history_reset_state(panel.hist)

    search_prompt, skeymap = LineEdit.setup_search_keymap(panel.hist)
    search_prompt.complete = REPL.LatexCompletions()

    panel.on_done = (s, buf, ok) -> begin
      if !ok
        LineEdit.transition(s, :abort)
        REPL.LineEdit.reset_state(s)
        return false
      end
      @msg working()

      line = String(take!(buf))

      if isempty(line)
        @msg doneWorking()
        return true
      end

      try
        r = interpret(line)
        r ≠ nothing && display(r)
      catch err
        display_error(stderr, err, stacktrace(catch_backtrace()))
      end
      println()
      LineEdit.reset_state(s)

      @msg doneWorking()
      @msg updateWorkspace()

      return true
    end

    panel.keymap_dict = LineEdit.keymap(Dict{Any,Any}[skeymap, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults])

    REPL.run_interface(Base.active_repl.t, REPL.LineEdit.ModalInterface([panel, search_prompt]))
  catch e
    @msg doneWorking()
    @msg updateWorkspace()
    e isa InterruptException || rethrow(e)
  end
end

# completions

struct JunoDebuggerRPELCompletionProvider <: REPL.CompletionProvider end

function LineEdit.complete_line(c::JunoDebuggerRPELCompletionProvider, s, state::DebuggerState = STATE)
  partial = REPL.beforecursor(s.input_buffer)
  full = LineEdit.input_string(s)

  frame = active_frame(state)

  # module-aware repl backend completions
  comps, range, should_complete = REPLCompletions.completions(full, lastindex(partial), moduleof(frame))
  ret = map(REPLCompletions.completion_text, comps)

  # local completions -- should be shown first
  @>> locals(frame) filter!(v -> begin
    # ref: https://github.com/JuliaDebug/JuliaInterpreter.jl/blob/master/src/utils.jl#L365-L370
    if v.name == Symbol("#self") && (v.value isa Type || sizeof(v.value) == 0)
      return false
    else
      return startswith(string(v.name), partial)
    end
  end) map(v -> string(v.name)) prepend!(ret) unique!

  return ret, partial[range], should_complete
end

# Evaluation
function interpret(code::AbstractString, state::DebuggerState = STATE)
  state.frame === nothing && return
  eval_code(active_frame(state), code)
end


text:
export
    CLG

"""
    CLG(paramdict::Dict)

Constructs an*sfunc representing a Conditional linear Gaussian.  These sfuncs may have
both discrete and continuous parents.  For each combination of discrete parents, there is
a `LinearGaussian` that depends on the continuous parents.

`CLG`s are implemented as a `Table` with a `LinearGaussian`.

The `paramdict` parameter defines the discrete and continuous parents, and the linear
gaussean values where the length of a key is the count of the discrete inputs, the
length of the tuple in a value is the count of continuous inputs, and the rest of the
values are used to build the parameters for `CLG` itself.  For example,

```
    Dict((:x,1) => ((-1.0, 1.0, 2.0), 3.0, 1.0),
         (:x,2) => ((-2.0, 4.0, 2.0), 3.0, 1.0),
         (:x,3) => ((-3.0, 2.0, 2.0), 3.0, 1.0),
         (:y,1) => ((-4.0, 5.0, 2.0), 3.0, 1.0),
         (:y,2) => ((-5.0, 3.0, 2.0), 3.0, 1.0),
         (:y,3) => ((-6.0, 6.0, 2.0), 3.0, 1.0))
```

    - the keys define two(2) discrete parents, with values `[:x,:y]` and `[1,2,3]`
    - in the values, the first tuple defines three(3) continuous parents for each
      underlying `LinearGausian`, with values `-1.0:-6.0`, `1.0:6.0`, and `2.0`
    - the values `3.0` and `1.0` are mean/stddev of the underlying `LinearGaussian`

See also: [`Table`](@ref), [`LinearGaussian`](@ref)
"""
function CLG(paramdict::Dict)
    numdiscreteinputs = length(collect(keys(paramdict))[1])
    numcontinuousinputs = length(collect(values(paramdict))[1][1])
    sfmaker(weights) = LinearGaussian(weights[1], weights[2], weights[3])
    return Table(NTuple{numcontinuousinputs, Float64}, Float64, numdiscreteinputs, paramdict, sfmaker)
end

#=

This representation is deprecated, but there may be useful implementation details.

# TODO (MRH): Convert CLG constructor into a call to a method that builds a Tables as in make_CLG above

struct CLG{NumDiscreteInputs, NumContinuousInputs, O} <:
    SFunc{Tuple{NTuple{NumDiscreteInputs, O},
                NTuple{NumContinuousInputs, Float64}},
          Float64,
          Tuple{}} # TODO (MRH): Params type
    # each segment corresponds to one setting of the discrete parents
    # and includes a linear weight for each continuous parent and a bias term
    segments :: Dict{<: Array{O, N} where N, <: Tuple{Array{Float64, N} where N, Float64}}
    variance :: Float64

    function CLG(num_discrete_inputs, num_continous_inputs, seg, var)
        T = typeof(seg).parameters[1].parameters[1]
        new{num_discrete_inputs, num_continous_inputs, T}(seg, var)
    end
end

#############################
#                           #
# Helper functions for CLGs #
#                           #
#############################

# Return all 2-tuples of discrete parent combinations and continuous
# parent combinations
function get_parent_combos(:: CLG{M,N}, parent_ranges) where {M,N}
    discrete_ranges = Array{Array{Any, 1}, 1}(parent_ranges[1:M])
    continuous_ranges = Array{Array{Any, 1}, 1}(parent_ranges[M+1:M+N])
    if isempty(discrete_ranges)
        discrete_combos = [[]]
    else
        discrete_combos = cartesian_product(discrete_ranges)
    end
    if isempty(continuous_ranges)
        continuous_combos = [[]]
    else
        continuous_combos = cartesian_product(continuous_ranges)
    end
    l = Array{Array{Any,1},1}([discrete_combos, continuous_combos])
    return map(a -> Tuple(a), cartesian_product(l))
end

function overlaps(int1, int2) :: Bool
    (l1, u1) = int1
    (l2, u2) = int2
    return l1 <= l2 && u1 >= l2 || l1 <= u2 && u1 >= u2
end

# Return the minimum and maximum probabilities of a normal distribution
# over an interval, when the mean is bounded and the variance is given.
function minmaxprob(lower, upper, lower_mean, upper_mean, variance)
    if lower == -Inf || upper == Inf return (0,1) end
    # The point can be anywhere between lower and upper
    # The mean can be anywhere between lower_mean and upper_mean
    # We need to find the minimum and maximum possible density
    # To find this, we compute the minimum and maximum possible distances
    # from the point to the mean.
    # If the point interval overlaps the mean interval, the minimum possible
    # distance is zero.
    # Otherwise the minimum possible distance is the min distance between
    # mean endpoints and interval endpoints.
    # The maximum possible distance is always the max distance between mean
    # endpoints and interval endpoints.
    d1 = abs(lower_mean - lower)
    d2 = abs(upper_mean - lower)
    d3 = abs(lower_mean - upper)
    d4 = abs(upper_mean - upper)
    dmin = overlaps((lower, upper), (lower_mean, upper_mean)) ? 0 :
             max(min(d1, d2, d3, d4), 0)
    dmax = max(d1, d2, d3, d4)
    densmin = normal_density(dmax, 0, variance)
    densmax = normal_density(dmin, 0, variance)
    diff = upper - lower
    pmin = min(max(densmin * diff, 0), 1)
    pmax = min(max(densmax * diff, 0), 1)
    return (pmin, pmax)
end

# Compute numerical lower and upper bounds on the density by partitioning
# the interval into num_partitions and computing bounds in each partition
function numerical_bounds(lower, upper, lmean, umean,
                          num_partitions, variance)
    if lower == -Inf || upper == Inf
        return (0,1)
    end
    start = lower
    step = (upper - lower) / num_partitions
    l = 0.0
    u = 0.0
    for i = 1:num_partitions
        (x,y) = minmaxprob(start, start + step, lmean, umean, variance)
        l += x
        u += y
        start += step
    end
    return (l,u)
end

#= ============================================================
   Operators
=============================================================== =#

function make_range(sf :: CLG, parranges, size :: Int)
    sd = sqrt(sf.variance)
    parent_combos = get_parent_combos(sf, parranges)

    # Prepare for accumulating the list of candidates by determining
    # whether the mean should be included for each segment and the number
    # of candidates around the mean. We choose the number of values per
    # segment so that the total is approximately the number of values desired.
    values_per_segment = max(div(size, length(parent_combos)), 1)
    if mod(values_per_segment, 2) == 1
        use_mean = true
        values_per_segment -= 1
    else
        use_mean = false
    end
    pairs_per_segment = values_per_segment / 2

    # Get candidates for including in the range
    # For each combinatiuon of parent values, we get the segment according
    # to the discrete values and the mean according to the continuous values.
    # We then spread candidates around the mean according to the standard
    # deviation.
    candidates = []
    for (discrete_combo, continuous_combo) in parent_combos
        (weights, bias) = sf.segments[discrete_combo]
        mean = linear_value(weights, bias, continuous_combo)
        if use_mean
            push!(candidates, mean)
            for i = 1:pairs_per_segment
                push!(candidates, mean - i * sd)
                push!(candidates, mean + i * sd)
            end
        else
            for i = 1:pairs_per_segment
                push!(candidates, mean - (i - 0.5) * sd)
                push!(candidates, mean + (i - 0.5) * sd)
            end
        end
    end
    sort!(candidates)

    # Select size of the candidates evenly spread
    skip = div(length(candidates), size)
    position = div(mod(length(candidates), size), 2) + 1
    result = []
    for i = 1:size
        push!(result, candidates[position])
        position += skip
    end

    return result
end

@op_impl begin
    struct CLGSupport{NumDiscreteInputs, NumContinuousInputs, O} <: Support{CLG{NumDiscreteInputs, NumContinuousInputs, O}} end
    function support((parranges, size, curr))
        if isempty(parranges)
            return nothing
        end

        old_size = length(curr)
        if old_size >= size
            return curr
        end
        result = make_range(sf, parranges, size)
        # We must make sure that any value in the current range is kept
        # so we replace the closest element in the new range with an element
        # from the current range
        i = 1
        j = 1
        is_current = fill(false, size)
        while i <= old_size && j <= size
            if curr[i] < result[j]
                if j == 1 || is_current[j-1] ||
                        result[j] - curr[i] < curr[i] - result[j-1]
                    result[j] = curr[i]
                    is_current[j] = true
                    i += 1
                    j += 1
                else
                    result[j-1] = curr[i]
                    is_current[j-1] = true
                    i += 1
                end
            else
                j += 1
            end
        end
        if i <= old_size
            tail = old_size - i
            result[size-tail:new_size] = curr[old_size-tail:old_size]
        end
        return result
    end
end

# We bound the integral under the density by dividing each interval of the
# range into num_partitions partitions and computing the lower and upper
# bound in each partition. However, this is complicated by the fact that
# we don't know the density function because we only have ranges for the
# parents.
@op_impl begin
    mutable struct CLGBoundedProbs <: BoundedProbs{CLG}
        numpartitions::Int64 = 100
    end

    function bounded_probs((range, parrange))
        vec = typeof(sf).parameters
        M = vec[1]
        N = vec[2]

        intervals = make_intervals(range)
        prs = Array{Array{Union{Symbol, Tuple{Float64, Float64}}, 1}, 1}(undef, M+N)
        for i = 1:M
            prs[i] = parrange[i]
        end
        for i = M+1:M+N
            prs[i] = make_intervals(parrange[i])
        end

        parent_combos = get_parent_combos(sf, prs)
        # Normally, the length of the range and the number of intervals is the same
        # However, for an empty range, there is still one unbounded interval
        if isempty(range)
            lower = Array{Float64}(undef, length(parent_combos))
            upper = Array{Float64}(undef, length(parent_combos))
        else
            lower = Array{Float64}(undef, length(parent_combos) * length(range))
            upper = Array{Float64}(undef, length(parent_combos) * length(range))
        end
        pos = 1
        for (discrete_combo, continuous_combo) in parent_combos
            (weights, bias) = sf.segments[discrete_combo]
            (lmean, umean) = bounded_linear_value(weights, bias, continuous_combo)
            ls = []
            us = []
            for interval in intervals
                (il, iu) = interval
                (l,u) = numerical_bounds(il, iu, lmean, umean,
                                op_impl.numpartitions, sf.variance)
                push!(ls, l)
                push!(us, u)
            end
            # We get better bounds by considering the bounds on other intervals.
            # This is especially important for intervals with -Inf or Inf as an
            # endpoint.
            for i = 1:length(intervals)
                otherls = 0.0
                otherus = 0.0
                for j = 1:i-1
                    otherls += ls[j]
                    otherus += us[j]
                end
                for j = i+1:length(intervals)
                    otherls += ls[j]
                    otherus += us[j]
                end
                l = max(ls[i], 1 - otherus)
                u = min(us[i], 1 - otherls)
                lower[pos] = l
                upper[pos] = u
                pos += 1
            end
        end
        return (lower, upper)
    end
end
=#


text:
function adapt_randn(rng::AbstractRNG, x::AbstractArray{<:ForwardDiff.Dual}, dims...)
    return adapt_randn(rng, ForwardDiff.valtype(eltype(x)), x, dims...)
end

## Binomial ##

function binomlogpdf(n::Int, p::ForwardDiff.Dual{T}, x::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val = ForwardDiff.value(p)
    Δ = ForwardDiff.partials(p)
    return FD(binomlogpdf(n, val, x),  Δ * (x / val - (n - x) / (1 - val)))
end

## Integer dual ##

function BetaBinomial(n::ForwardDiff.Dual{<:Any, <:Integer}, α::Real, β::Real; check_args = true)
    return BetaBinomial(ForwardDiff.value(n), α, β; check_args = check_args)
end
Binomial(n::ForwardDiff.Dual{<:Any, <:Integer}, p::Real) = Binomial(ForwardDiff.value(n), p)
function Erlang(α::ForwardDiff.Dual{<:Any, <:Integer}, θ::Real; check_args = true)
    return Erlang(ForwardDiff.value(α), θ, check_args = check_args)
end

## Negative binomial ##

function nbinomlogpdf(r::ForwardDiff.Dual{T}, p::ForwardDiff.Dual{T}, k::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val_p = ForwardDiff.value(p)
    val_r = ForwardDiff.value(r)
    Δ_r = ForwardDiff.partials(r)
    dr = _nbinomlogpdf_grad_1(val_r, val_p, k)
    Δ_p = ForwardDiff.partials(p)
    dp = _nbinomlogpdf_grad_2(val_r, val_p, k)
    Δ = ForwardDiff._mul_partials(Δ_r, Δ_p, dr, dp)
    return FD(nbinomlogpdf(val_r, val_p, k),  Δ)
end
function nbinomlogpdf(r::Real, p::ForwardDiff.Dual{T}, k::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val_p = ForwardDiff.value(p)
    Δ_p = ForwardDiff.partials(p) * _nbinomlogpdf_grad_2(r, val_p, k)
    return FD(nbinomlogpdf(r, val_p, k),  Δ_p)
end
function nbinomlogpdf(r::ForwardDiff.Dual{T}, p::Real, k::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val_r = ForwardDiff.value(r)
    Δ_r = ForwardDiff.partials(r) * _nbinomlogpdf_grad_1(val_r, p, k)
    return FD(nbinomlogpdf(val_r, p, k),  Δ_r)
end

## ForwardDiff broadcasting support ##
# If we use Distributions >= 0.24, then `DISTRIBUTIONS_HAS_GENERIC_UNIVARIATE_PDF` is `true`.
# In Distributions 0.24 `logpdf` is defined for inputs of type `Real` which are then
# converted to the support of the distributions (such as integers) in their concrete implementations.
# Thus it is no needed to have a special function for dual numbers that performs the conversion
# (and actually this method leads to method ambiguity errors since even discrete distributions now
# define logpdf(::MyDistribution, ::Real), see, e.g.,
# JuliaStats/Distributions.jl@ae2d6c5/src/univariate/discrete/binomial.jl#L119).
if !DISTRIBUTIONS_HAS_GENERIC_UNIVARIATE_PDF
    @eval begin
        function Distributions.logpdf(d::DiscreteUnivariateDistribution, k::ForwardDiff.Dual)
            return logpdf(d, convert(Integer, ForwardDiff.value(k)))
        end
    end
end


text:
using DoubleExponentialFormulas
using Plots, LaTeXStrings
using Dierckx
using DelimitedFiles, Printf
using TwoFAST

# %% Specify a redshift
redshift = [0, 0.3, 0.5, 0.7, 1, 1.3, 1.5]
ired = 2
z = redshift[ired]

# %% Read in linear power spectra computed by CLASS and spline interpolate in k
filename = @sprintf("data/matterpower05_z%1d_pk.dat", ired)
d = readdlm(filename, comments = true)
pk = Spline1D(d[:, 1], d[:, 2])

# %% Compute correlation function with `quaddeo`
kmin, kmax = 1e-4, Inf
R1 = Float64.(10:2:200)
ξ1 = zeros(length(R1))
for i = 1:length(R1)
   func1(k) = k^2 * pk(k) * sin(k * R1[i]) / (k * R1[i]) / 2 / π^2
   ξ1[i], err = quaddeo(func1, R1[i], 0, kmin, kmax)
end

# %% Compute correlation function with `xicalc` in TwoFAST
R2, ξ2 = xicalc(pk, 0, 0)

# %% Plot results and save to "compare_xi_twofast.pdf"
p = plot(
   R1,
   R1 .^ 2 .* ξ1,
   m = 2,
   c = :red,
   lab = "Computed with quaddeo",
   ylab = L"R^2\xi(R)",
   xlab = L"R~[Mpc/h]",
   title = "Redshift = $z",
   xlims = [0,200],
)
p = plot!(R2, R2 .^ 2 .* ξ2, ls = :dash, c = :red, lab = "with TwoFAST")
savefig("compare_xi_twofast.pdf")
display(p)


text:
module tanH
    using LoopVectorization

    function func!(output_matrix::Array{Float32}, value_matrix::Array{Float32})
        @avxt for i in eachindex(value_matrix)
            value_matrix[i] = ifelse(value_matrix[i]>3.0f38, 3.0f38, value_matrix[i])
            value_matrix[i] = ifelse(value_matrix[i]<-3.0f38, -3.0f38, value_matrix[i])
            output_matrix[i] = tanh(value_matrix[i])
        end
    end

    function get_∇biases!(∇biases::Array{Float32}, value_matrix::Array{Float32}, δ::Array{Float32})
        @avxt for i in eachindex(value_matrix)
            ∇biases[i] = (1-tanh(value_matrix[i])^2)*δ[i]
        end
    end

    function get_name()
        return "tanH"
    end
end


text:
#Test code
#-------------------------------------------------------------------------------

using GracePlot

@warn("No real test code yet... just ensuring that \"using\" works.")

@info("Run sample/runsamples.jl for sample usage.")

:Test_Complete


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, DataFrames, Plots
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.15"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
binomialRV(n,p) = sum(rand(n) .< p)

# ╔═╡ 3c3fe880-2162-11eb-2efd-b53935597f04
begin
	p, n, N = 0.5, 10, 10^6
end;

# ╔═╡ 3c400b1c-2162-11eb-2196-09023e52b6af
begin
	bDist = Binomial(n,p)
	xGrid = 0:n
	bPmf = [pdf(bDist,i) for i in xGrid]
	data = [binomialRV(n,p) for _ in 1:N]
	pmfEst = counts(data,0:n)/N
end

# ╔═╡ 3c412d30-2162-11eb-3931-e9ae8415ba38
begin
	plot( xGrid, pmfEst,
		line=:stem, marker=:circle,
		c=:blue, ms=10, msw=0, lw=4, label="MC estimate")
	plot!( xGrid, bPmf,
		 line=:stem, marker=:xcross, c=:red,
		 ms=6, msw=0, lw=2, label="PMF", xticks=(0:1:10),
		 ylims=(0,0.3), xlabel="x", ylabel="Probability")
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.15"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═3c3fe880-2162-11eb-2efd-b53935597f04
# ╠═3c400b1c-2162-11eb-2196-09023e52b6af
# ╠═3c412d30-2162-11eb-3931-e9ae8415ba38
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
using DelimitedFiles:readdlm
using Random:shuffle!
using Flux
using Plots
using CUDA

using DrugRepositioning

CUDA.allowscalar(false)

"原始数据存储路径"
data_dir=joinpath(@__DIR__,"data","origin")

"对比方法测试结果"
compare_dir=joinpath(@__DIR__,"data","all_result")
ms=compared_scores(compare_dir)

"高玲的方法的测试结果"
gao_dir=joinpath(@__DIR__,"data","all_result","gao")
P=readdlm(joinpath(gao_dir,"2","0mean5_P.txt"))[:]
FPR=readdlm(joinpath(gao_dir,"2","0mean5_FPR.txt"))[:]
TPR=readdlm(joinpath(gao_dir,"2","0mean5_TPR.txt"))[:]
R=TPR
pushfirst!(ms,MetricsScore("GPred_2",R,P,FPR,TPR,[],[]))
P=readdlm(joinpath(gao_dir,"1","0mean5_P.txt"))[:]
FPR=readdlm(joinpath(gao_dir,"1","0mean5_FPR.txt"))[:]
TPR=readdlm(joinpath(gao_dir,"1","0mean5_TPR.txt"))[:]
R=TPR
pushfirst!(ms,MetricsScore("GPred_1",R,P,FPR,TPR,[],[]))

"读取CSV"
function read_csv(file_path;d_type::DataType=Float64,i_type::DataType=String,h_dtype::DataType=String)
    f=readdlm(file_path,',')
    header=f[1,2:end]|>Array{String}
    indexer=f[2:end,1]|>Array{String}
    data=f[2:end,2:end]|>Array{d_type}
    indexer,header,data
end

"药物疾病关联"
function drug_disease()
    return read_csv(joinpath(data_dir,"drug_disease_mat.csv");d_type=Bool)
end

"疾病相似性"
function disease_similarity()
    return read_csv(joinpath(data_dir,"disease_similarity_mat.csv");d_type=Float32)
end

"药物靶标域"
drug_domain()=read_csv(joinpath(data_dir,"drug_target_domain_mat.csv");d_type=Bool)
"药物靶标GO"
drug_go()=read_csv(joinpath(data_dir,"drug_target_go_mat.csv");d_type=Bool)
"药物公共化学子结构"
drug_pubchem()=read_csv(joinpath(data_dir,"drug_pubchem_mat.csv");d_type=Bool)

drugs,diseases,RD=drug_disease()
_,_,D1=disease_similarity()
_,_,RP1=drug_domain()
_,_,RP2=drug_go()
_,_,RP3=drug_pubchem()
_,_,RP4=drug_disease()

props=Dict{Symbol,Array{Matrix}}(
    :drug=>[RP1',RP2',RP3',RP4'],
    :disease=>[]
)

dataset=Dataset(drugs,diseases,RD,props)
add_similarity!(dataset,:disease,D1)
model=CombineModel(dataset;feature_left=128,feature_right=512)|>gpu
# trainer=Trainer(model,Flux.Losses.logitcrossentropy,dataset;epochs=80)
# training(trainer,Flux.Optimise.ADAM(3e-4))
res=ValidateResult("result")
df=analyze_by_drugs(res)
R,P,FPR,TPR=pr_roc(res)
my_score=MetricsScore("MSRD",R,P,FPR,TPR,[],[])
pushfirst!(ms,my_score)
pr_roc_curve(ms)|>display
plot(res)

max_stars_repo_name: yeruoforever/DrugRepositioning

max_issues_repo_name: yeruoforever/DrugRepositioning

max_forks_repo_name: yeruoforever/DrugRepositioning


text:
"Abstract barrier function type"
abstract type BarrierFunction <: CertificateFunction end

"""
    ControlBarrierFunction

Control barrier function h for a control affine system.
# Fields
- `h`: function h(x) that represents the CBF
- `∇h`: function ∇h(x) defining the gradient of the CBF
- `α`: extended class K function α(h(x))
"""
struct ControlBarrierFunction <: BarrierFunction
    h
    ∇h
    α
end

"""
	ControlBarrierFunction(h, α)

Construct a CBF given h(x) and extended class K function α(h(x)).
"""
function ControlBarrierFunction(h, α)
    ∇h(x) = gradient(h, x)[1]'
    return ControlBarrierFunction(h, ∇h, α)
end

"""
	(h::ControlBarrierFunction)(x)

Evaluate CBF at state x.
"""
function (h::ControlBarrierFunction)(x)
    return h.h(x)
end


text:
using CompScienceMeshes
using BEAST

m = CompScienceMeshes.tetmeshsphere(1.0, 0.45)
bnd_m = boundary(m)

m1 = skeleton(m,1)
X = BEAST.nedelecc3d(m,m1)

Y = BEAST.ttrace(X,bnd_m)
@assert length(geometry(Y)) == length(bnd_m)


text:
@testset "VM 1D1V with spin" begin
    σ, μ = 0.02, 0.0
    kx, α = 1.004355, 0.001
    xmin, xmax = 0, 2π / kx
    domain = [xmin, xmax, xmax - xmin]
    nx = 8
    mesh = OneDGrid(xmin, xmax, nx)
    n_particles = 2
    mesh = OneDGrid(xmin, xmax, nx)
    spline_degree = 3

    df = CosSumGaussianSpin([[kx]], [α], [[σ]], [[μ]])

    particle_group2 = ParticleGroup{1,1}(n_particles; n_spin=3)
    sampler = ParticleSampler{1,1}(:sobol, n_particles)

    sample!(particle_group2, sampler, df, mesh)

    ref = [
        3.127970342747129 4.691955514120693
        0.01734694403902491 -0.018034876317136343
        6.249684744808763 6.255940685494258
        0.0 0.0
        0.0 0.0
        1.0 1.0
    ]

    @test particle_group2.array ≈ ref

    particle_group = ParticleGroup{1,1}(n_particles; n_spin=3)

    for i_part in 1:n_particles
        x = zeros(1)
        v = zeros(1)
        s = zeros(3)
        w = zeros(1)
        x = GEMPIC.get_x(particle_group2, i_part)
        v = GEMPIC.get_v(particle_group2, i_part)

        s1 = GEMPIC.get_spin(particle_group2, i_part, 1)
        s2 = GEMPIC.get_spin(particle_group2, i_part, 2)
        s3 = GEMPIC.get_spin(particle_group2, i_part, 3)

        w = GEMPIC.get_weights(particle_group2, i_part)

        GEMPIC.set_x!(particle_group, i_part, x[1])
        GEMPIC.set_v!(particle_group, i_part, v[1])
        GEMPIC.set_spin!(particle_group, i_part, 1, s1)
        GEMPIC.set_spin!(particle_group, i_part, 2, s2)
        GEMPIC.set_spin!(particle_group, i_part, 3, s3)
        GEMPIC.set_weights!(particle_group, i_part, w[1])
    end

    kernel_smoother0 = ParticleMeshCoupling1D(mesh, n_particles, spline_degree, :galerkin)
    kernel_smoother1 = ParticleMeshCoupling1D(
        mesh, n_particles, spline_degree - 1, :galerkin
    )
    kernel_smoother2 = ParticleMeshCoupling1D(
        mesh, n_particles, spline_degree - 2, :galerkin
    )

    rho = zeros(Float64, nx)
    efield_poisson = zeros(Float64, nx)

    # Init!ialize the field solver
    maxwell_solver = Maxwell1DFEM(mesh, spline_degree)
    # efield by Poisson
    solve_poisson!(efield_poisson, particle_group, kernel_smoother0, maxwell_solver, rho)

    ref = [
        -0.5811006688338523,
        -0.8336082614158025,
        -2.7872165528999013,
        1.3710914975949415,
        -1.3738807953831014,
        2.7885897256327588,
        0.8341900701471621,
        0.5819349851577953,
    ]

    @test efield_poisson ≈ ref

    # Initialize the arrays for the spline coefficients of the fields
    k0 = 12.0523
    E0 = 10
    ww = 12.104827940833333

    Ey(x) = E0 * cos(k0 * x)
    Ez(x) = E0 * sin(k0 * x)
    Ay(x) = -E0 / ww * sin(k0 * x)
    Az(x) = E0 / ww * cos(k0 * x)

    efield_dofs = [efield_poisson, zeros(Float64, nx), zeros(Float64, nx)]
    afield_dofs = [zeros(Float64, nx), zeros(Float64, nx)]

    l2projection!(efield_dofs[2], maxwell_solver, Ey, spline_degree)
    l2projection!(efield_dofs[3], maxwell_solver, Ez, spline_degree)
    l2projection!(afield_dofs[1], maxwell_solver, Ay, spline_degree)
    l2projection!(afield_dofs[2], maxwell_solver, Az, spline_degree)

    @test efield_dofs[2] ≈ [
        -0.17074650510231487,
        0.17074650468466593,
        -0.17074650409994616,
        0.1707465033481691,
        -0.17074650242932413,
        0.17074650134342606,
        -0.1707465000904745,
        0.1707464986704385,
    ]
    @test efield_dofs[3] ≈ [
        -1.0681800094815476e-5,
        1.6022700129154877e-5,
        -2.136360014868184e-5,
        2.670450014853168e-5,
        -3.204540010485847e-5,
        3.738630007371513e-5,
        -4.272719996293515e-5,
        4.806809980886038e-5,
    ]
    @test afield_dofs[1] ≈ [
        8.824412992125462e-7,
        -1.3236619477267236e-6,
        1.7648825950343306e-6,
        -2.2061032407091485e-6,
        2.6473238828043315e-6,
        -3.088544525898473e-6,
        3.5297651624492664e-6,
        -3.970985795411606e-6,
    ]
    @test afield_dofs[2] ≈ [
        -0.014105653210181833,
        0.01410565317567925,
        -0.014105653127374468,
        0.014105653065269032,
        -0.01410565298936165,
        0.014105652899653832,
        -0.014105652796145382,
        0.014105652678833834,
    ]

    propagator = HamiltonianSplittingSpin(
        maxwell_solver,
        kernel_smoother0,
        kernel_smoother1,
        kernel_smoother2,
        particle_group,
        efield_dofs,
        afield_dofs,
    )

    efield_dofs_n = propagator.e_dofs

    thdiag = TimeHistoryDiagnosticsSpin(
        particle_group, maxwell_solver, kernel_smoother0, kernel_smoother1
    )

    steps, Δt = 2, 0.1

    # Strang splitting
    strang_splitting!(propagator, Δt, 1)

    solve_poisson!(efield_poisson, particle_group, kernel_smoother0, maxwell_solver, rho)

    @test efield_poisson ≈ [-0.579, -0.836, -2.782, 1.379, -1.381, 2.784, 0.836, 0.580] atol =
        1e-2

    strang_splitting!(propagator, Δt, 1)

    solve_poisson!(efield_poisson, particle_group, kernel_smoother0, maxwell_solver, rho)

    @test efield_poisson ≈ [-0.572, -0.848, -2.759, 1.417, -1.420, 2.761, 0.848, 0.573] atol =
        1e-2

    ref = [
        3.117 4.702
        -0.125 0.124
        6.249 6.255
        -6.256e-5 6.112e-5
        3.774e-9 3.532e-9
        1.000 1.000
    ]
    @test propagator.particle_group.array ≈ ref atol = 1e-2

    write_step!(
        thdiag,
        Δt,
        spline_degree,
        efield_dofs,
        afield_dofs,
        efield_dofs_n,
        efield_poisson,
        propagator,
    )

    GEMPIC.save("particles", 1, particle_group)
end


text:
module test_bukdu_multipart

using Test
using Bukdu
using HTTP: HTTP, Form, Multipart

function preview(multipart::Multipart)
    mark(multipart)
    data = read(multipart)
    reset(multipart)
    [multipart.filename, String(data)]
end

routes() do
    post("/upload") do conn::Conn
        pre = preview.((conn.params.user_file1, conn.params.user_file2))
        render(Text, pre)
    end
end

Bukdu.start(8191, host="127.0.0.1")

part1 = Multipart("mycoolfile1.txt", IOBuffer("hello"))
part2 = Multipart("mycoolfile2.txt", IOBuffer("world"))
r = HTTP.post("http://127.0.0.1:8191/upload", Form(["user_file1" => part1, "user_file2" => part2]))
@test r.body == Vector{UInt8}("""(["mycoolfile1.txt", "hello"], ["mycoolfile2.txt", "world"])""")

seekstart(part2.data)
r = HTTP.post("http://127.0.0.1:8191/upload", Form(["user_file1" => Multipart("", IOBuffer()), "user_file2" => part2]))
@test r.body == Vector{UInt8}("""(["", ""], ["mycoolfile2.txt", "world"])""")

Bukdu.stop()

Routing.empty!()

end # module test_bukdu_multipart


text:
fl = ILt(s -> 1/s^3, talbot)
@test isapprox( fl(1.0), 0.5; atol =  1e-9)
@test isapprox( fl.([1.0,2.0]), [0.5, 2.0]; atol =  1e-9)

fl = ILt(s -> 1/s^3, gwr,8)
@test isapprox( fl(1.0), 0.5; atol = 1e-5)

p = ILtPair(Weeks(s -> s/(1+s^2)), cos)

## Optim is broken
# e1 = abserr(p, 10.0)
# optimize(p,10.0)
# e2 = abserr(p, 10.0)
# @test e2 < e1

@test typeof(Weeks(iltpair_power(5))) <: ILtPair

# Default method is talbot
fl1 = ILt( s -> 1/s)
@test fl1.iltfunc == talbot


text:
using QuantLib

mutable struct Problem{F <: CostFunction, C <: Constraint, T}
  costFunction::F
  constraint::C
  initialValue::Vector{T}
  currentValue::Vector{T}
  functionValue::Float64
  squaredNorm::Float64
  functionEvaluation::Int
  gradientEvaluation::Int
end

function Problem(costFunction::F, constraint::C, initialValue::Vector{T}) where {F <: CostFunction, C <: Constraint, T}
  currentValue = initialValue
  functionValue = 0.0
  squaredNorm = 0.0
  functionEvaluation = 0
  gradientEvaluation = 0

  return Problem{F, C, T}(costFunction, constraint, initialValue, currentValue, functionValue, squaredNorm, functionEvaluation, gradientEvaluation)
end

## Problem methods ##
function reset!(p::Problem)
  p.functionEvaluation = p.gradientEvaluation = 0
  p.functionValue = p.squaredNorm = 0.0

  return p
end

function value!(p::Problem, x::Vector{T}) where {T}
  p.functionEvaluation += 1
  return QuantLib.value(p.costFunction, x)
end

function values!(p::Problem, x::Vector{Float64})
  p.functionEvaluation += 1
  return QuantLib.func_values(p.costFunction, x)
end

function extrapolate!(p::Problem, i_highest::Int, factor::Float64, values::Vector{T}, sum_array::Vector{T},
                    vertices::Vector{Vector{T}}) where {T}
  pTry = zeros(length(sum_array))
  while true
    dimensions = length(values) - 1
    factor1 = (1.0 - factor) / dimensions
    factor2 = factor1 - factor
    @inbounds pTry = sum_array * factor1 - vertices[i_highest] * factor2
    factor *= 0.5

    if test(p.constraint, pTry) || abs(factor) <= EPS_VAL
      break
    end
  end

  if abs(factor) <= EPS_VAL
    return values[i_highest]
  end
  factor *= 2
  vTry = value!(p, pTry) # TODO check this
  if vTry < values[i_highest]
    @inbounds values[i_highest] = vTry
    @inbounds sum_array[:] = sum_array + pTry - vertices[i_highest]
    @inbounds vertices[i_highest] = pTry
  end
  # if vTry == 0.02226748007400235
  #   println("sum_array: $sum_array")
  #   println("values: $values")
  #   println("pTry: $pTry")
  #   println("vertices: $vertices")
  #   error("Ooops")
  # end
  return vTry
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1LoadBalancerStatus <: SwaggerModel
    ingress::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1LoadBalancerIngress} } # spec name: ingress

    function IoK8sApiCoreV1LoadBalancerStatus(;ingress=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1LoadBalancerStatus, Symbol("ingress"), ingress)
        setfield!(o, Symbol("ingress"), ingress)
        o
    end
end # type IoK8sApiCoreV1LoadBalancerStatus

const _property_map_IoK8sApiCoreV1LoadBalancerStatus = Dict{Symbol,Symbol}(Symbol("ingress")=>Symbol("ingress"))
const _property_types_IoK8sApiCoreV1LoadBalancerStatus = Dict{Symbol,String}(Symbol("ingress")=>"Vector{IoK8sApiCoreV1LoadBalancerIngress}")
Base.propertynames(::Type{ IoK8sApiCoreV1LoadBalancerStatus }) = collect(keys(_property_map_IoK8sApiCoreV1LoadBalancerStatus))
Swagger.property_type(::Type{ IoK8sApiCoreV1LoadBalancerStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1LoadBalancerStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1LoadBalancerStatus }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1LoadBalancerStatus[property_name]

function check_required(o::IoK8sApiCoreV1LoadBalancerStatus)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1LoadBalancerStatus }, name::Symbol, val)
end


text:
using RandomWalker
using Test

@testset "RandomWalker.jl" begin
    # Write your tests here.
    # Write here
end


text:
module EclipseREPL
function equivalent(command)
  if(strip(command)=="?")
    return helpcode("?")
  end
  hcode = replace(command, r"^\s*\?", "")
  if hcode != command
    return helpcode(hcode)
  end
  return command
end
function validate(statement)
  expression = parse(statement,1; greedy=true, raise=false)
  if !isa(expression[1],Expr)
    return (:complete,statement)
  end
  if expression[1].head==:error
    return (expression[1].head, expression[1].args[1])
  end
  return (expression[1].head,nothing)
end
function execute(command)
  status = nothing
  mimeType=nothing
  try
    statement = equivalent(command)
    state,message = validate(statement)
    if(state==:error)
      status="error"
      println(message)
      return
    end
    if(state==:incomplete)
      status="incomplete"
      return
    end
    status="complete"
    result=include_string(statement)
    if result==nothing
      return
    end
    response=stringmime("text/plain",result)
    resultType = string(typeof(result))
    if resultType=="Gadfly.Plot"
      mimeType = "text/html"
      println(stringmime(mimeType,result))
      return
    end
    if contains(response,"matplotlib")
      mimeType = "text/html"
      eval(:($(Expr(:using, symbol("PyPlot")))))
      figure = PyPlot.gcf()
      pygui(false) # after the fact - review
      println(stringmime("image/svg+xml",figure))
      return
    end
    if rstrip(statement)[end]!=';'
       mimeType="text/plain"
        println(response)
    end
  catch e
    showerror(STDOUT, e); println()
    status = "error"
    println()
  finally
    flush_all()
    println()
    println("<<<<$status>>>>")
    println("<<<<$mimeType>>>>")
    println("<<<<eox>>>>")
  end
end
function flush_all()
    Libc.flush_cstdio()
    flush(STDOUT)
    flush(STDERR)
end
using Base
function Base.workspace()
    last = Core.Main
    b = last.Base
    e = EclipseREPL
    ccall(:jl_new_main_module, Any, ())
    m = Core.Main
    ccall(:jl_add_standard_imports, Void, (Any,), m)
    eval(m,
         Expr(:toplevel,
              :(const Base = $(Expr(:quote, b))),
              :(const LastMain = $(Expr(:quote, last))),
              :(const EclipseREPL = $(Expr(:quote, e)))))
    empty!(Base.package_locks)
    nothing
end
function helpcode(code::AbstractString)
    if VERSION < v"0.4.0-dev+2891" # old Base.@help macro
        return "Base.@help " * code
    else # new Base.Docs.@repl macro from julia@08663d4bb05c5b8805a57f46f4feacb07c7f2564
        code_ = strip(code)
        # as in base/REPL.jl, special-case keywords so that they parse
        if(haskey(Docs.keywords, symbol(code_)))
          return "eval(:(Base.Docs.@repl \$(symbol(\"$code_\"))))"
        else
          return "Base.Docs.@repl $code_"
        end
    end
end
end

max_stars_repo_name: JuliaEditorSupport/JuliaDT

max_issues_repo_name: JuliaEditorSupport/JuliaDT

max_forks_repo_name: JuliaComputing/JuliaDT


text:
@testset "HyperSphere" begin
    @test Sphere(Point(0,0,0), 4) == HyperSphere{3,Int}(Point{3,Int}((0,0,0)),4)
    c = Circle(Point(0,0), 1)
    @test c == HyperSphere(Point(0,0), 1)
    @test origin(c) == Point(0,0)
    @test !( isinside(c, 1, 2) )
    @test isinside(c, 0.5, 0.5)

    centered_rect = centered(HyperSphere)
    @test centered_rect == HyperSphere{3,Float32}(Point3f0(0), 0.5f0)
    centered_rect = centered(HyperSphere{2})
    @test centered_rect == HyperSphere{2,Float32}(Point2f0(0),0.5f0)

    centered_rect = centered(HyperSphere{2, Float64})
    @test centered_rect == HyperSphere{2,Float64}(Point(0.,0.), 0.5)

    centered_rect = centered(HyperSphere{3, Float32})
    @test centered_rect == HyperSphere{3,Float32}(Point3f0(0), 0.5f0)

    @test widths(centered_rect) == Vec3f0(1)
    @test radius(centered_rect) == 0.5f0
    @test maximum(centered_rect) == Vec3f0(0.5f0)
    @test minimum(centered_rect) == Vec3f0(-0.5f0)
    s = Sphere(Point3f0(0), 1f0)
    f = decompose(Face{2, Int}, s, 3)
    # TODO 54 linesegments for 3 facets is too much.
    @test length(f) == 54
end


text:
import Base: @deprecate

## deprecated in v0.16

@deprecate(Bcast!(buffer, count::Integer, root::Integer, comm::Comm),
           Bcast!(view(buffer, 1:count), root, comm), false)

@deprecate(Scatter!(sendbuf, recvbuf, count::Integer, root::Integer, comm::Comm),
           Scatter!(UBuffer(sendbuf, count), recvbuf, root, comm), false)
@deprecate(Scatter!(sendbuf, recvbuf::AbstractArray, count::Integer, root::Integer, comm::Comm),
           Scatter!(UBuffer(sendbuf, count), view(recvbuf, 1:count), root, comm), false)
@deprecate(Scatter!(sendbuf, recvbuf::Nothing, count::Integer, root::Integer, comm::Comm),
           Scatter!(UBuffer(sendbuf, count), MPI.IN_PLACE, root, comm), false)
@deprecate(Scatter!(sendbuf::Nothing, recvbuf, count::Integer, root::Integer, comm::Comm),
           Scatter!(nothing, recvbuf, root, comm), false)
@deprecate(Scatter!(sendbuf::Nothing, recvbuf::AbstractArray, count::Integer, root::Integer, comm::Comm),
           Scatter!(nothing, view(recvbuf, 1:count), root, comm), false)

@deprecate(Scatter(sendbuf, count::Integer, root::Integer, comm::Comm),
           MPI.Comm_rank(comm) == root ?
           Scatter!(UBuffer(sendbuf, count), similar(sendbuf, count), root, comm) :
           Scatter!(nothing, similar(sendbuf, count), root, comm), false)

@deprecate(Scatterv!(sendbuf, recvbuf, counts::Vector, root::Integer, comm::Comm),
           MPI.Comm_rank(comm) == root ?
           Scatterv!(VBuffer(sendbuf, counts), recvbuf, root, comm) :
           Scatterv!(nothing, recvbuf, root, comm), false)
@deprecate(Scatterv!(sendbuf, recvbuf::AbstractArray, counts::Vector, root::Integer, comm::Comm),
           MPI.Comm_rank(comm) == root ?
           Scatterv!(VBuffer(sendbuf, counts), view(recvbuf,1:counts[MPI.Comm_rank(comm)+1]), root, comm) :
           Scatterv!(nothing, view(recvbuf,1:counts[MPI.Comm_rank(comm)+1]), root, comm), false)
@deprecate(Scatterv!(sendbuf, recvbuf::Nothing, counts::Vector, root::Integer, comm::Comm),
           MPI.Comm_rank(comm) == root ?
           Scatterv!(VBuffer(sendbuf, counts), MPI.IN_PLACE, root, comm) :
           Scatterv!(nothing, MPI.IN_PLACE, root, comm), false)
@deprecate(Scatterv!(sendbuf::Nothing, recvbuf, counts::Vector, root::Integer, comm::Comm),
           Scatterv!(nothing, recvbuf, root, comm), false)
@deprecate(Scatterv!(sendbuf::Nothing, recvbuf::AbstractArray, counts::Vector, root::Integer, comm::Comm),
           Scatterv!(nothing, view(recvbuf,1:counts[MPI.Comm_rank(comm)+1]), root, comm), false)

@deprecate(Scatterv(sendbuf, counts::Vector, root::Integer, comm::Comm),
           MPI.Comm_rank(comm) == root ?
           Scatterv!(VBuffer(sendbuf, counts), similar(sendbuf, counts[Comm_rank(comm) + 1]), root, comm) :
           Scatterv!(nothing, similar(sendbuf, counts[Comm_rank(comm) + 1]), root, comm), false)

@deprecate(Gather!(sendbuf, recvbuf, count::Integer, root::Integer, comm::Comm),
           Gather!(sendbuf, MPI.Comm_rank(comm) == root ? UBuffer(recvbuf, count) : nothing, root, comm), false)
@deprecate(Gather!(sendbuf::AbstractArray, recvbuf, count::Integer, root::Integer, comm::Comm),
           Gather!(view(sendbuf, 1:count), MPI.Comm_rank(comm) == root ? UBuffer(recvbuf, count) : nothing, root, comm), false)
@deprecate(Gather!(sendbuf::Nothing, recvbuf, count::Integer, root::Integer, comm::Comm),
           Gather!(MPI.IN_PLACE, UBuffer(recvbuf, count), root, comm), false)
@deprecate(Gather!(sendbuf, recvbuf::Nothing, count::Integer, root::Integer, comm::Comm),
           Gather!(sendbuf, nothing, root, comm), false)
@deprecate(Gather!(sendbuf::AbstractArray, recvbuf::Nothing, count::Integer, root::Integer, comm::Comm),
           Gather!(view(sendbuf, 1:count), nothing, root, comm), false)

@deprecate(Gather(sendbuf, count::Integer, root::Integer, comm::Comm),
           Gather(sendbuf, root, comm), false)
@deprecate(Gather(sendbuf::AbstractArray, count::Integer, root::Integer, comm::Comm),
           Gather(view(sendbuf, 1:count), root, comm), false)

@deprecate(Gatherv!(sendbuf, recvbuf, counts::Vector{Cint}, root::Integer, comm::Comm),
           Gatherv!(sendbuf, MPI.Comm_rank(comm) == root ? VBuffer(recvbuf, counts) : nothing, root, comm), false)
@deprecate(Gatherv!(sendbuf::AbstractArray, recvbuf, counts::Vector{Cint}, root::Integer, comm::Comm),
           Gatherv!(view(sendbuf,1:counts[MPI.Comm_rank(comm)+1]), MPI.Comm_rank(comm) == root ? VBuffer(recvbuf, counts) : nothing, root, comm), false)
@deprecate(Gatherv!(sendbuf::Nothing, recvbuf, counts::Vector{Cint}, root::Integer, comm::Comm),
           Gatherv!(MPI.IN_PLACE, VBuffer(recvbuf, counts), root, comm), false)
@deprecate(Gatherv!(sendbuf, recvbuf::Nothing, counts::Vector{Cint}, root::Integer, comm::Comm),
           Gatherv!(sendbuf, nothing, root, comm), false)
@deprecate(Gatherv!(sendbuf::AbstractArray, recvbuf::Nothing, counts::Vector{Cint}, root::Integer, comm::Comm),
           Gatherv!(view(sendbuf,1:counts[MPI.Comm_rank(comm)+1]), nothing, root, comm), false)

@deprecate(Gatherv(sendbuf, counts::Vector{Cint}, root::Integer, comm::Comm),
           Gatherv!(sendbuf, Comm_rank(comm) == root ? VBuffer(similar(sendbuf, sum(counts)), counts) : nothing, root, comm), false)
@deprecate(Gatherv(sendbuf::AbstractArray, counts::Vector{Cint}, root::Integer, comm::Comm),
           Gatherv!(view(sendbuf,1:counts[MPI.Comm_rank(comm)+1]), Comm_rank(comm) == root ? VBuffer(similar(sendbuf, sum(counts)), counts) : nothing, root, comm), false)

@deprecate(Allgather!(sendbuf, recvbuf, count::Integer, comm::Comm),
           Allgather!(sendbuf, UBuffer(recvbuf, count), comm::Comm), false)
@deprecate(Allgather!(sendbuf::AbstractArray, recvbuf, count::Integer, comm::Comm),
           Allgather!(view(sendbuf,1:count), UBuffer(recvbuf, count), comm::Comm), false)
@deprecate(Allgather!(sendbuf::Nothing, recvbuf, count::Integer, comm::Comm),
           Allgather!(MPI.IN_PLACE, UBuffer(recvbuf, count), comm::Comm), false)
@deprecate(Allgather!(sendrecvbuf, count::Integer, comm::Comm),
           Allgather!(UBuffer(sendrecvbuf, count), comm), false)

@deprecate(Allgather(sendbuf, count::Integer, comm::Comm),
           Allgather(view(sendbuf, 1:count), comm), false)

@deprecate(Allgatherv!(sendbuf, recvbuf, counts::Vector{Cint}, comm::Comm),
           Allgatherv!(sendbuf, VBuffer(recvbuf, counts), comm), false)
@deprecate(Allgatherv!(sendbuf::AbstractArray, recvbuf, counts::Vector{Cint}, comm::Comm),
           Allgatherv!(view(sendbuf,1:counts[MPI.Comm_rank(comm)+1]), VBuffer(recvbuf, counts), comm), false)
@deprecate(Allgatherv!(sendrecvbuf, counts::Vector{Cint}, comm::Comm),
           Allgatherv!(VBuffer(sendrecvbuf, counts), comm), false)

@deprecate(Allgatherv(sendbuf, counts::Vector{Cint}, comm::Comm),
           Allgatherv!(sendbuf, VBuffer(similar(sendbuf, sum(counts)), counts), comm), false)

@deprecate(Alltoall!(sendbuf, recvbuf, count::Integer, comm::Comm),
           Alltoall!(UBuffer(sendbuf, count), UBuffer(recvbuf, count), comm), false)
@deprecate(Alltoall!(sendbuf::InPlace, recvbuf, count::Integer, comm::Comm),
           Alltoall!(sendbuf, UBuffer(recvbuf, count), comm), false)
@deprecate(Alltoall!(sendrecvbuf, count::Integer, comm::Comm),
           Alltoall!(UBuffer(sendrecvbuf, count), comm), false)

@deprecate(Alltoall(sendbuf, count::Integer, comm::Comm),
           Alltoall(UBuffer(sendbuf, count), comm), false)

@deprecate(Alltoallv!(sendbuf, recvbuf, sendcounts, recvcounts, comm::Comm),
           Alltoallv!(VBuffer(sendbuf, sendcounts), VBuffer(recvbuf, recvcounts), comm), false)

@deprecate(Alltoallv(buf, sendcounts, recvcounts, comm::Comm),
           Alltoallv!(VBuffer(buf, sendcounts), VBuffer(similar(buf,sum(recvcounts)), recvcounts), comm), false)

@deprecate(Reduce!(send::AbstractArray, recv::AbstractArray, count::Integer, op, root::Integer, comm::Comm),
           Reduce!(view(send, 1:count), view(recv, 1:count), op, root, comm), false)
@deprecate(Reduce!(send, recv::AbstractArray, count::Integer, op, root::Integer, comm::Comm),
           Reduce!(send, view(recv, 1:count), op, root, comm), false)
@deprecate(Reduce!(send::AbstractArray, recv, count::Integer, op, root::Integer, comm::Comm),
           Reduce!(view(send, 1:count), recv, op, root, comm), false)
@deprecate(Reduce!(send, recv, count::Integer, op, root::Integer, comm::Comm),
           Reduce!(send, recv, op, root, comm), false)
@deprecate(Reduce!(sendrecvbuf::AbstractArray, count::Integer, op, root::Integer, comm::Comm),
           Reduce!(view(sendrecvbuf, 1:count), op, root, comm), false)
@deprecate(Reduce!(sendrecvbuf, count::Integer, op, root::Integer, comm::Comm),
           Reduce!(sendrecvbuf, op, root, comm), false)

@deprecate(Allreduce!(send::AbstractArray, recv::AbstractArray, count::Integer, op, comm::Comm),
           Allreduce!(view(send, 1:count), view(recv, 1:count), op, comm), false)
@deprecate(Allreduce!(send, recv::AbstractArray, count::Integer, op, comm::Comm),
           Allreduce!(send, view(recv, 1:count), op, comm), false)
@deprecate(Allreduce!(send::AbstractArray, recv, count::Integer, op, comm::Comm),
           Allreduce!(view(send, 1:count), recv, op, comm), false)
@deprecate(Allreduce!(send, recv, count::Integer, op, comm::Comm),
           Allreduce!(send, recv, op, comm), false)
@deprecate(Allreduce!(sendrecvbuf::AbstractArray, count::Integer, op, comm::Comm),
           Allreduce!(view(sendrecvbuf, 1:count), op, comm), false)
@deprecate(Allreduce!(sendrecvbuf, count::Integer, op, comm::Comm),
           Allreduce!(sendrecvbuf, op, comm), false)

@deprecate(Scan!(sendbuf::AbstractArray, recvbuf::AbstractArray, count::Integer, op, comm::Comm),
           Scan!(view(sendbuf, 1:count), view(recvbuf, 1:count), op, comm), false)
@deprecate(Scan!(sendbuf, recvbuf::AbstractArray, count::Integer, op, comm::Comm),
           Scan!(sendbuf, view(recvbuf, 1:count), op, comm), false)
@deprecate(Scan!(sendbuf::AbstractArray, recvbuf, count::Integer, op, comm::Comm),
           Scan!(view(sendbuf, 1:count), recvbuf, op, comm), false)
@deprecate(Scan!(sendbuf, recvbuf, count::Integer, op, comm::Comm),
           Scan!(sendbuf, recvbuf, op, comm), false)
@deprecate(Scan!(sendrecvbuf::AbstractArray, count::Integer, op, comm::Comm),
           Scan!(view(sendrecvbuf, 1:count), op, comm), false)
@deprecate(Scan!(sendrecvbuf, count::Integer, op, comm::Comm),
           Scan!(sendrecvbuf, op, comm), false)

@deprecate(Exscan!(sendbuf::AbstractArray, recvbuf::AbstractArray, count::Integer, op, comm::Comm),
           Exscan!(view(sendbuf, 1:count), view(recvbuf, 1:count), op, comm), false)
@deprecate(Exscan!(sendbuf, recvbuf::AbstractArray, count::Integer, op, comm::Comm),
           Exscan!(sendbuf, view(recvbuf, 1:count), op, comm), false)
@deprecate(Exscan!(sendbuf::AbstractArray, recvbuf, count::Integer, op, comm::Comm),
           Exscan!(view(sendbuf, 1:count), recvbuf, op, comm), false)
@deprecate(Exscan!(sendbuf, recvbuf, count::Integer, op, comm::Comm),
           Exscan!(sendbuf, recvbuf, op, comm), false)
@deprecate(Exscan!(sendrecvbuf::AbstractArray, count::Integer, op, comm::Comm),
           Exscan!(view(sendrecvbuf, 1:count), op, comm), false)
@deprecate(Exscan!(sendrecvbuf, count::Integer, op, comm::Comm),
           Exscan!(sendrecvbuf, op, comm), false)

@deprecate(Get(origin_buffer, count::Integer, target_rank::Integer, target_disp::Integer, win::Win),
           Get(view(origin_buffer, 1:count), target_rank, target_disp, win), false)
@deprecate(Get(origin_buffer::Ref, count::Integer, target_rank::Integer, target_disp::Integer, win::Win),
           Get(origin_buffer, target_rank, target_disp, win), false)
@deprecate(Put(origin_buffer, count::Integer, target_rank::Integer, target_disp::Integer, win::Win),
           Put(view(origin_buffer, 1:count), target_rank, target_disp, win), false)
@deprecate(Put(origin_buffer::Ref, count::Integer, target_rank::Integer, target_disp::Integer, win::Win),
           Put(origin_buffer, target_rank, target_disp, win), false)
@deprecate(Accumulate(origin_buffer, count::Integer, target_rank::Integer, target_disp::Integer, op::Op, win::Win),
           Accumulate(view(origin_buffer, 1:count), target_rank, target_disp, op, win), false)
@deprecate(Accumulate(origin_buffer::Ref, count::Integer, target_rank::Integer, target_disp::Integer, op::Op, win::Win),
           Accumulate(origin_buffer, target_rank, target_disp, op, win), false)
@deprecate(Get_accumulate(origin_buffer, result_buffer, count::Integer, target_rank::Integer, target_disp::Integer, op::Op, win::Win),
           Get_accumulate(view(origin_buffer,1:count), view(result_buffer,1:count), target_rank, target_disp, op, win), false)
@deprecate(Get_accumulate(origin_buffer::Ref, result_buffer::Ref, count::Integer, target_rank::Integer, target_disp::Integer, op::Op, win::Win),
           Get_accumulate(origin_buffer, result_buffer, target_rank, target_disp, op, win), false)


text:
# 1:2*length(layerconfig) -> encoder lstm (2k-1 for weight and 2k for bias in k'th encoder layer)
# 2*length(layerconfig)+1:4*length(layerconfig) -> decoder lstm (2k-1+2length for weight and 2k+2length for bias in kth decoder layer)
function initparams(atype, layerconfig, encembed, decembed, invocab, outvocab, winit)
    len = length(layerconfig)
    parameters = Array(Any, 4*length(layerconfig) + 4)

    # embedding and prediction layers
    parameters[end] = winit * randn(outvocab, decembed) # decoder embedding
    parameters[end-1] = winit * randn(invocab, encembed) # encoder embedding
    parameters[end-2] = winit * randn(layerconfig[end], outvocab) # final layer weight
    parameters[end-3] = zeros(1, outvocab) # final layer bias

    # encoder lstm
    input = encembed
    for k=1:len
        parameters[2k-1] = winit * randn(input+layerconfig[k], 4layerconfig[k])
        parameters[2k] = zeros(1, 4layerconfig[k])
        parameters[2k][1:layerconfig[k]] = 1 # forget gate bias
        input = layerconfig[k]
    end

    # decoder lstm
    input2 = decembed
    for k=1:len
        parameters[2k-1+2len] = winit * randn(input2+layerconfig[k], 4*layerconfig[k])
        parameters[2k+2len] = zeros(1, 4layerconfig[k])
        parameters[2k+2len][1:layerconfig[k]] = 1 # forget gate bias
        input2 = layerconfig[k]
    end
    return map(p->convert(atype, p), parameters)
end

# 1:2length(layerconfig) -> encoder states (k is cell and k-1 is hidden for kth layer)
# 2length(layerconfig)+1:4length(layerconfig) -> decoder states
# first hidden layer of the decoder left nothing intentionally
function initstates(atype, layerconfig, batchsize)
    len = length(layerconfig)
    states = Array(Any, 4len)

    # initialize encoder states
    for k=1:len
        states[2k-1] = zeros(batchsize, layerconfig[k])
        states[2k] = zeros(batchsize, layerconfig[k])
    end

    # initialize decoder states
    for k=1:len
        states[2k-1+2len] = zeros(batchsize, layerconfig[k])
        states[2k + 2len] = zeros(batchsize, layerconfig[k])
    end
    states = map(s->convert(atype, s), states)
    states[2len+1] = nothing
    return states
end

function lstm(weight, bias, hidden, cell, input)
    gates   = hcat(input,hidden) * weight .+ bias
    hsize   = size(hidden,2)
    forget  = sigm(gates[:,1:hsize])
    ingate  = sigm(gates[:,1+hsize:2hsize])
    outgate = sigm(gates[:,1+2hsize:3hsize])
    change  = tanh(gates[:,1+3hsize:end])
    cell    = cell .* forget + ingate .* change
    hidden  = outgate .* tanh(cell)
    return (hidden,cell)
end

# multi layer lstm forward function
function forward(parameters, states, input)
    x = input
    for i=1:2:length(states)
        (states[i], states[i+1]) = lstm(parameters[i], parameters[i+1], states[i], states[i+1], x)
        x = states[i]
    end
    return x
end

function loss(parameters, states, inputseq, outseq)
    total = 0.0
    count = 0.0
    atype = typeof(AutoGrad.getval(parameters[1]))
    hlayers = length(states)/4
    hlayers = convert(Int, hlayers)

    # go through the encoder
    enstates = states[1:2hlayers]
    final_hidden = nothing
    for i=1:length(inputseq)
        input = convert(atype, inputseq[i])
        x = input * parameters[end-1]
        final_hidden = forward(parameters[1:2hlayers], enstates, x)
    end

    # initialize the hidden state of the decoder, decstates[1] is the initial hidden of decoder lstm
    decstates = states[2hlayers+1:4hlayers]
    decstates[1] = final_hidden

    # go through the decoder
    for i=1:length(outseq)-1
        input = convert(atype, outseq[i])
        x = input * parameters[end]
        hidden = forward(parameters[2hlayers+1:4hlayers], decstates, x)
        ypred = hidden * parameters[end-2] .+ parameters[end-3]
        ynorm = logp(ypred, 2)
        ygold = convert(atype, outseq[i+1])
        count += size(ygold, 1)
        total += sum(ygold .* ynorm)
    end
    return - total / count
end

lossgradient = grad(loss)


text:
module RocketCatchErrorOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: catch_error()" begin

    println("Testing: operator catch_error()")

    run_proxyshowcheck("CatchError", catch_error((err, obs) -> never()), args = (check_subscription = true, ))

    run_testset([
        (
            source = from(1:5) |> safe() |> map(Int, (d) -> d == 4 ? error(4) : d) |> catch_error((err, obs) -> from(1:3)),
            values = @ts([ 1, 2, 3, 1, 2, 3, c ])
        ),
        (
            source = from(1:5) |> async(0) |> safe() |> map(Int, (d) -> d == 4 ? error(4) : d) |> catch_error((err, obs) -> obs) |> take(10),
            values = @ts([ 1 ] ~ [ 2 ] ~ [ 3 ] ~ [ 1 ] ~ [ 2 ] ~ [ 3 ] ~ [ 1 ] ~ [ 2 ] ~ [ 3 ] ~ [ 1, c ])
        ),
        (
            source = from(1:5) |> switch_map(Int, (d) -> d == 4 ? faulted(Int, 4) : of(d)) |> catch_error((err, obs) -> of(5)),
            values = @ts([ 1, 2, 3, 5, c ])
        ),
        (
            source = from(1:5) |> safe() |> map(Int, (d) -> d == 4 ? error(4) : d) |> catch_error((err, obs) -> completed(Int)),
            values = @ts([ 1, 2, 3, c ])
        ),
        (
            source = from(1:5) |> safe() |> map(Int, (d) -> d == 4 ? error(4) : d) |> catch_error((err, obs) -> never(Int)),
            values = @ts([ 1, 2, 3 ])
        ),
        (
            source = from(1:5) |> safe() |> map(Int, (d) -> d == 4 ? error(4) : d) |> catch_error((err, obs) -> throw("e")),
            values = @ts([ 1, 2, 3, e("e") ])
        )
    ])
end

end


text:
using Base.Test
import ArchGDAL; const AG = ArchGDAL

function read(f, filename)
    return AG.registerdrivers() do
        AG.read(filename) do dataset
            f(dataset)
end end end

@testset "Testing Displays for different objects" begin
    read("data/point.geojson") do dataset
        print(dataset)
    end;

    read("data/point.geojson") do dataset
        print(AG.getlayer(dataset, 0))
    end;

    read("data/point.geojson") do dataset
        AG.getfeature(AG.getlayer(dataset, 0), 2) do feature
           print(feature)
        end
    end;

    read("gdalworkshop/world.tif") do dataset
        print(dataset)
    end;

    read("gdalworkshop/world.tif") do dataset
        print(AG.getband(dataset, 1))
    end;
end


text:
#Code written by Skipper af Rosenborg and Esben Bager
function createDataSplits(XArrays, YArrays, t, trainingSize)
    trainingXArrays = Array{Array{Float64, 2}}(industriesTotal)
    trainingYArrays = Array{Array{Float64, 2}}(industriesTotal)

    validationXRows = Array{Array{Float64, 2}}(industriesTotal)
    validationY     = Array{Array{Float64, 2}}(industriesTotal)

    OOSXArrays      = Array{Array{Float64, 2}}(industriesTotal)
    OOSYArrays      = Array{Array{Float64, 2}}(industriesTotal)

    OOSRow          = Array{Array{Float64, 2}}(industriesTotal)
    OOSY            = Array{Array{Float64, 2}}(industriesTotal)
    for i = 1:industriesTotal
        trainingXArrays[i] = XArrays[i][t:(t-1+trainingSize), :]
        trainingYArrays[i] = YArrays[i][t:(t-1+trainingSize), :]

        validationXRows[i] = XArrays[i][(t+trainingSize):(t+trainingSize), :]
        validationY[i]     = YArrays[i][(t+trainingSize):(t+trainingSize), :]

        OOSXArrays[i]      = XArrays[i][(t+1):(t+trainingSize), :]
        OOSYArrays[i]      = YArrays[i][(t+1):(t+trainingSize), :]

        OOSRow[i]          = XArrays[i][(t+trainingSize+1):(t+trainingSize+1), :]
        OOSY[i]            = YArrays[i][(t+trainingSize+1):(t+trainingSize+1), :]
    end
    return trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY
end

function calculatePvalue(return1N, returnPerfect, returnCEO)
    nominator = returnPerfect-returnCEO
    if returnCEO < 0
        nominator = abs(returnCEO-returnPerfect)
    end

    denominator = returnPerfect-return1N
    if return1N < 0
        denominator = abs(return1N-returnPerfect)
    end

    PValue = 1 - nominator/denominator
    return PValue
end

function findPerfectResults(trainingXArrays, Xrow, Yvalues, gamma)
    indexes = 10
    Sigma =  cov(trainingXArrays[1][:,1:10])

    #A=U^(T)U where U is upper triangular with real positive diagonal entries
    F = lufact(Sigma)
    U = F[:U]  #Cholesky factorization of Sigma

    M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0, Threads=1))
    @variables M begin
            w[1:indexes]
            u[1:indexes]
            z
            y
    end

    @objective(M,Min, gamma*y - Yvalues'*w)
    @constraint(M, 0 .<= w)
    @constraint(M, sum(w[i] for i=1:indexes) == 1)
    @constraint(M, norm([2*U'*w;y-1]) <= y+1)
    solve(M)
    wPerfect = getvalue(w)
    forecastRow = (exp.(Xrow')-1)*100
    periodPerfectReturn = forecastRow*wPerfect

    return wPerfect, periodPerfectReturn
end

function findPerfectResultsRFR(trainingXArrays, Xrow, Yvalues, gamma)
    indexes = 11
    Sigma =  cov(trainingXArrays[1][:,1:11])

    #A=U^(T)U where U is upper triangular with real positive diagonal entries
    F = lufact(Sigma)
    U = F[:U]  #Cholesky factorization of Sigma

    M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0, Threads=1))
    @variables M begin
            w[1:indexes]
            u[1:indexes]
            z
            y
    end

    @objective(M,Min, gamma*y - Yvalues'*w)
    @constraint(M, 0 .<= w)
    @constraint(M, sum(w[i] for i=1:indexes) == 1)
    @constraint(M, norm([2*U'*w;y-1]) <= y+1)
    solve(M)
    wPerfect = getvalue(w)
    forecastRow = zeros(11)
	forecastRow[1:10] = (exp.(Xrow[1:10])-1)*100
	forecastRow[11]   = (exp.(Xrow[11])-1)
    periodPerfectReturn = forecastRow*wPerfect

    return wPerfect, periodPerfectReturn
end

function performMVOptimization(expectedReturns, U, gamma, Xrow, Yvalues)
    indexes = 10
    M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0, Threads=1))
    @variables M begin
            w[1:indexes]
            u[1:indexes]
            z
            y
    end

    @objective(M,Min, gamma*y - expectedReturns'*w)
    @constraint(M, 0 .<= w)
    @constraint(M, sum(w[i] for i=1:indexes) == 1)
    @constraint(M, norm([2*U'*w;y-1]) <= y+1)
    solve(M)
    wStar = getvalue(w)

    forecastRow = (exp.(Xrow)-1)*100

    periodReturn = forecastRow'*wStar
    period1NReturn = forecastRow'*w1N

    return period1NReturn, periodReturn, wStar, forecastRow
end

function performMVOptimizationRISK(expectedReturns, U, gamma, Xrow, Yvalues)
    indexes = 11
    M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0))
    @variables M begin
            w[1:indexes]
            y
    end

    @objective(M,Min, gamma*y - expectedReturns'*w)
    @constraint(M, 0 .<= w)
    @constraint(M, sum(w[i] for i=1:indexes) == 1)
    @constraint(M, norm([2*U'*w;y-1]) <= y+1)
    solve(M)
    wStar = getvalue(w)
    u = getvalue(y)

    forecastRow = zeros(11)
	forecastRow[1:10] = (exp.(Xrow[1:10])-1)*100
	forecastRow[11]   = (exp.(Xrow[11])-1)

    periodReturn = forecastRow'*wStar
    period1NReturn = forecastRow'*w1N

    return period1NReturn, periodReturn, wStar, forecastRow, u
end

function generateExpectedReturns(betaArray, trainingXArrays, trainingYArrays, validationXRows)
    industriesTotal = 10

    expectedReturns = zeros(10)
    bRows = 240
    bCols = floor.(Int,zeros(industriesTotal))
    for i = 1:industriesTotal
        bCols[i] = size(trainingXArrays[i])[2]
    end

    for i = 1:10
        errorArray = (trainingXArrays[i]*betaArray[i,1:bCols[i]]-trainingYArrays[i])*(1/bRows)
        errorSum = sum(errorArray)
        expectedReturns[i] = validationXRows[i][:]'*betaArray[i,1:bCols[i]]-errorSum
    end
    return expectedReturns
end

function runCEO(trainingXArrays, trainingYArrays, modelConfigRow, gamma)
    industriesTotal = 10
    periodMean = zeros(10)

    l1, l2, l3, l4 = modelConfigRow

    maxPredictors = 1412
    bCols = floor.(Int,zeros(industriesTotal))
    Sigma =  cov(trainingXArrays[1][:,1:10])

    periodMean = zeros(10)
    for i = 1:10
        periodMean[i] = mean(trainingYArrays[i][:,1])
    end

    periodCap = zeros(10)
    for i = 1:10
        periodCap[i] = maximum(trainingYArrays[i][:])
    end

    #A=U^(T)U where U is upper triangular with real positive diagonal entries
    F = lufact(Sigma)
    U = F[:U]  #Cholesky factorization of Sigma

    for i = 1:industriesTotal
        bCols[i] = size(trainingXArrays[i])[2]
    end
    gamma = gamma
    bRows = size(trainingXArrays[1])[1]
    M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0, Threads=1))

    @variables M begin
            b[1:industriesTotal, 1:maxPredictors] #Beta in LASSO; industry i and predictors 1:1304
            z[1:industriesTotal, 1:maxPredictors] #auxilliary for 1norm
            q[1:industriesTotal] #auxilliary for norm
            f[1:industriesTotal, 1:bRows] #expected return for industry i at time t
            w[1:industriesTotal, 1:bRows] #weights for industry i at time t
            y[1:bRows] #auxilliary variables for MV optimization
    end

    #@objective(M,Min,sum(0.5*q[i] + l2*ones(maxPredictors)'*z[i,:] for i=1:industriesTotal)+ l3*(sum(gamma*y[t]-f[:,t]'*w[:,t] for t=1:bRows)) - l4*(sum(w[i,t]*trainingYArrays[i][t] for i=1:10, t=1:bRows)))

    @objective(M,Min,sum(l1*q[i] + l2*ones(maxPredictors)'*z[i,:] for i=1:industriesTotal)+ l3*(sum(gamma*y[t]-sum(f[i,t] for i=1:10) for t=1:bRows)) - l4*(sum(w[i,t]*trainingYArrays[i][t] for i=1:10, t=1:bRows)))

    for i = 1:industriesTotal
        for t = 1:bRows
            @constraint(M, (periodMean[i]+periodCap[i])*w[i,t] >= f[i,t])
        end
    end
    for i = 1:industriesTotal
        @constraint(M, norm( [1-q[i];2*(trainingXArrays[i]*b[i,1:bCols[i]]-trainingYArrays[i])] ) <= 1+q[i]) #second order cone constraint SOC
    end
    @constraint(M,  b .<= z)
    @constraint(M, -z .<= b)

    for i = 1:10
        errorArray = (trainingXArrays[i]*b[i,1:bCols[i]]-trainingYArrays[i])*(1/bRows)
        errorSum = sum(errorArray)
        for t = 1:bRows
            prediction = trainingXArrays[i][t, :]'*b[i,1:bCols[i]]
            @constraint(M, prediction - errorSum >= f[i, t]) ##== originally
        end
    end

    @constraint(M, 0 .<= w)
    for t = 1:bRows
        @constraint(M, sum(w[i,t] for i=1:industriesTotal) == 1)
        @constraint(M, norm([2*U'*w[:,t];y[t]-1]) <= y[t]+1)
    end

    @time(solve(M))
    betaArray = getvalue(b)
    fVars = getvalue(f)
    yvars = getvalue(y)
    #Insert shrinking here
    return betaArray, U
end

function runCEORFR(trainingXArrays, trainingYArrays, modelConfigRow, gamma, rfRatesVec)
    industriesTotal = 10

    l1, l2, l3, l4 = modelConfigRow

    maxPredictors = 1412
    bCols = floor.(Int,zeros(industriesTotal))
    trainX = hcat(trainingXArrays[1][:,1:10], rfRatesVec)
    Sigma =  cov(trainX)

    periodMean = zeros(11)
    for i = 1:10
        periodMean[i] = mean(trainingYArrays[i][:,1])
    end

    periodCap = zeros(11)
    for i = 1:10
        periodCap[i] = maximum(trainingYArrays[i][:])
    end
    periodCap[11] = maximum(rfRatesVec)

    periodMean[11] = mean(rfRatesVec)

    #A=U^(T)U where U is upper triangular with real positive diagonal entries
    F = lufact(Sigma)
    U = F[:U]  #Cholesky factorization of Sigma

    for i = 1:industriesTotal
        bCols[i] = size(trainingXArrays[i])[2]
    end
    gamma = gamma
    bRows = size(trainingXArrays[1])[1]
    M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0, Threads=1))

    @variables M begin
            b[1:industriesTotal, 1:maxPredictors] #Beta in LASSO; industry i and predictors 1:1412
            z[1:industriesTotal, 1:maxPredictors] #auxilliary for 1norm
            q[1:industriesTotal] #auxilliary for norm
            f[1:(industriesTotal+1), 1:bRows] #expected return for industry i at time t
            w[1:(industriesTotal+1), 1:bRows] #weights for industry i at time t
            y[1:bRows] #auxilliary variables for MV optimization
    end

    #@objective(M,Min,sum(0.5*q[i] + l2*ones(maxPredictors)'*z[i,:] for i=1:industriesTotal)+ l3*(sum(gamma*y[t]-f[:,t]'*w[:,t] for t=1:bRows)) - l4*(sum(w[i,t]*trainingYArrays[i][t] for i=1:10, t=1:bRows)))

    @objective(M,Min,sum(l1*q[i] + l2*ones(maxPredictors)'*z[i,:] for i=1:industriesTotal)+ l3*(sum(gamma*y[t]-sum(f[i,t] for i=1:11) for t=1:bRows)) - l4*(sum(w[11,t]*rfRatesVec[t] for t=1:bRows) + sum(w[i,t]*trainingYArrays[i][t] for i=1:10, t=1:bRows)))

    for i = 1:(industriesTotal+1)
        for t = 1:bRows
            @constraint(M, (periodMean[i]+periodCap[i])*w[i,t] >= f[i,t])
        end
    end
    for i = 1:industriesTotal
        @constraint(M, norm( [1-q[i];2*(trainingXArrays[i]*b[i,1:bCols[i]]-trainingYArrays[i])] ) <= 1+q[i]) #second order cone constraint SOC
    end
    @constraint(M,  b .<= z)
    @constraint(M, -z .<= b)

    for i = 1:10
        errorArray = (trainingXArrays[i]*b[i,1:bCols[i]]-trainingYArrays[i])*(1/bRows)
        errorSum = sum(errorArray)
        for t = 1:bRows
            prediction = trainingXArrays[i][t, :]'*b[i,1:bCols[i]]
            @constraint(M, prediction - errorSum >= f[i, t]) ##== originally
        end
    end

    #RFR
    for t = 1:bRows
        @constraint(M, rfRatesVec[t] >= f[11, t])
    end

    @constraint(M, 0 .<= w)
    for t = 1:bRows
        @constraint(M, sum(w[i,t] for i=1:(industriesTotal+1)) == 1)
        @constraint(M, norm([2*U'*w[:,t];y[t]-1]) <= y[t]+1)
    end

    @time(solve(M))
    betaArray = getvalue(b)
    #Insert shrinking here
    return betaArray, U
end

function generateXandYs(industries, modelMatrix)
    industriesTotal = length(industries)
    for i = 1:industriesTotal
        industry = industries[i]
        mainData = loadIndexDataLOGReturn(industry, path)

        #path = "/zhome/9f/d/88706/SpecialeCode/"
        #path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/"
        dateAndRecession = Array(mainData[:,end-1:end])
        mainDataArr = Array(mainData[:,1:end-2])

        colNames = names(mainData)

        nRows = size(mainDataArr)[1]
        nCols = size(mainDataArr)[2]

        fileName = path*"Results/IndexData/CEOTest/"*industry*"/"*string(trainingSize)*"-1/"*string(trainingSize)*"_"

        mainYarr = mainDataArr[:, nCols:nCols]

        mainXarr = mainDataArr[:,1:nCols-1]

        VIX = modelMatrix[i, 1]
        raw = modelMatrix[i, 2]
        expTrans = modelMatrix[i, 3]
        timeTrans = modelMatrix[i, 4]
        TA = modelMatrix[i,5]

        if raw == 1
            fileName = fileName*"Raw"
            mainXarr = mainXarr[:,1:10]
        else
            if VIX == 1
                fileName = fileName*"VIX_Macro"
                mainXarr = mainXarr[:,1:end-1]
            else
                fileName = fileName*"Macro"
            end
        end

        # Transform with time elements
        if timeTrans == 1
            fileName = fileName*"Time"
            mainXarr = expandWithTime3612(mainXarr)
        end

        # Transform #
        if expTrans == 1
            fileName = fileName*"Exp"
            mainXarr = expandWithTransformations(mainXarr)
        end

        if TA == 1
            fileName = fileName*"TA"
            mainXarr = expandWithMAandMomentum(mainXarr, mainYarr, (nCols-1))
        end

        # Standardize
        standX = zScoreByColumn(mainXarr)
        standY = mainYarr

        XArrays[i] = standX
        YArrays[i] = standY
    end
    return XArrays, YArrays
end

max_stars_repo_name: skipperafrosenborg/Thesis

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

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
module TestFinSets
using Test

using Catlab.Theories
using Catlab.CategoricalAlgebra.ShapeDiagrams, Catlab.CategoricalAlgebra.FinSets

# Category of finite ordinals
#############################

f = FinOrdFunction([1,3,4], 5)
g = FinOrdFunction([1,1,2,2,3], 3)
h = FinOrdFunction([3,1,2], 3)

# Evaluation.
@test map(f, 1:3) == [1,3,4]
@test map(id(FinOrd(3)), 1:3) == [1,2,3]
@test map(FinOrdFunction(x -> (x % 3) + 1, 3, 3), 1:3) == [2,3,1]

# Composition and identities.
@test dom(f) == FinOrd(3)
@test codom(f) == FinOrd(5)
@test compose(f,g) == FinOrdFunction([1,2,2], 3)
@test compose(g,h) == FinOrdFunction([3,3,1,1,2], 3)
@test compose(compose(f,g),h) == compose(f,compose(g,h))
@test force(compose(id(dom(f)),f)) == f
@test force(compose(f,id(codom(f)))) == f

# Limits
########

# Terminal object.
@test terminal(FinOrd) == FinOrd(1)

# Binary Product.
span = product(FinOrd(2), FinOrd(3))
@test apex(span) == FinOrd(6)
@test force(left(span)) == FinOrdFunction([1,2,1,2,1,2])
@test force(right(span)) == FinOrdFunction([1,1,2,2,3,3])

# N-ary Product.
cone = product([FinOrd(2), FinOrd(3)])
@test apex(cone) == FinOrd(6)
@test force(leg(cone,1)) == FinOrdFunction([1,2,1,2,1,2])
@test force(leg(cone,2)) == FinOrdFunction([1,1,2,2,3,3])
@test apex(product(FinOrd[])) == FinOrd(1)

# Equalizer.
f, g = FinOrdFunction([1,2,3]), FinOrdFunction([3,2,1])
@test equalizer(f,g) == FinOrdFunction([2], 3)
@test equalizer([f,g]) == FinOrdFunction([2], 3)

# Equalizer in case of identical functions.
f = FinOrdFunction([4,2,3,1], 5)
@test equalizer(f,f) == force(id(FinOrd(4)))
@test equalizer([f,f]) == force(id(FinOrd(4)))

# Equalizer matching nothing.
f, g = id(FinOrd(5)), FinOrdFunction([2,3,4,5,1], 5)
@test equalizer(f,g) == FinOrdFunction(Int[], 5)
@test equalizer([f,g]) == FinOrdFunction(Int[], 5)

# Pullback.
span = pullback(Cospan(FinOrdFunction([1,1,3,2],4), FinOrdFunction([1,1,4,2],4)))
@test apex(span) == FinOrd(5)
@test force(left(span)) == FinOrdFunction([1,2,1,2,4], 4)
@test force(right(span)) == FinOrdFunction([1,1,2,2,4], 4)

# Pullback from a singleton set: the degenerate case of a product.
span = pullback(Cospan(FinOrdFunction([1,1]), FinOrdFunction([1,1,1])))
@test apex(span) == FinOrd(6)
@test force(left(span)) == FinOrdFunction([1,2,1,2,1,2])
@test force(right(span)) == FinOrdFunction([1,1,2,2,3,3])

# Pullback using generic limit interface
f,g = FinOrdFunction([1,1,3,2],4), FinOrdFunction([1,1,4,2],4)
cone = limit(Diagram([FinOrd(4),FinOrd(4),FinOrd(4)], [(1,3,f),(2,3,g)]))
@test apex(cone) == FinOrd(5)
@test force(leg(cone,1)) == FinOrdFunction([1,2,1,2,4],4)
@test force(leg(cone,2)) == FinOrdFunction([1,1,2,2,4],4)

# Colimits
##########

# Initial object.
@test initial(FinOrd) == FinOrd(0)

# Binary Coproduct.
cospan = coproduct(FinOrd(2), FinOrd(3))
@test base(cospan) == FinOrd(5)
@test left(cospan) == FinOrdFunction([1,2], 5)
@test right(cospan) == FinOrdFunction([3,4,5], 5)

# N-ary Coproduct.
cocone = coproduct([FinOrd(2), FinOrd(3)])
@test base(cocone) == FinOrd(5)
@test leg(cocone,1) == FinOrdFunction([1,2], 5)
@test leg(cocone,2) == FinOrdFunction([3,4,5], 5)
@test base(coproduct(FinOrd[])) == FinOrd(0)

# Coequalizer from a singleton set.
f, g = FinOrdFunction([1], 3), FinOrdFunction([3], 3)
@test coequalizer(f,g) == FinOrdFunction([1,2,1], 2)
@test coequalizer([f,g]) == FinOrdFunction([1,2,1], 2)

# Coequalizer in case of identical functions.
f = FinOrdFunction([4,2,3,1], 5)
@test coequalizer(f,f) == force(id(FinOrd(5)))
@test coequalizer([f,f]) == force(id(FinOrd(5)))

# Coequalizer identifying everything.
f, g = id(FinOrd(5)), FinOrdFunction([2,3,4,5,1], 5)
@test coequalizer(f,g) == FinOrdFunction(repeat([1],5))
@test coequalizer([f,g]) == FinOrdFunction(repeat([1],5))

# Pushout from the empty set: the degenerate case of the coproduct.
f, g = FinOrdFunction(Int[], 2), FinOrdFunction(Int[], 3)
cospan = pushout(Span(f,g))
@test base(cospan) == FinOrd(5)
@test left(cospan) == FinOrdFunction([1,2], 5)
@test right(cospan) == FinOrdFunction([3,4,5], 5)

# Pushout from a singleton set.
f, g = FinOrdFunction([1], 2), FinOrdFunction([2], 3)
cospan = pushout(Span(f,g))
@test base(cospan) == FinOrd(4)
h, k = left(cospan), right(cospan)
@test compose(f,h) == compose(g,k)
@test h == FinOrdFunction([1,2], 4)
@test k == FinOrdFunction([3,1,4], 4)

# Same thing with generic colimit interface
diag = Diagram([FinOrd(1),FinOrd(2),FinOrd(3)],[(1,2,f), (1,3,g)])
cocone = colimit(diag)
@test base(cocone) == FinOrd(4)
h, k = leg(cocone,2), leg(cocone,3)
@test compose(f,h) == compose(g,k)
@test h == FinOrdFunction([1,2], 4)
@test k == FinOrdFunction([3,1,4], 4)

# Pushout from a two-element set, with non-injective legs.
f, g = FinOrdFunction([1,1], 2), FinOrdFunction([1,2], 2)
cospan = pushout(Span(f,g))
@test base(cospan) == FinOrd(2)
h, k = left(cospan), right(cospan)
@test compose(f,h) == compose(g,k)
@test h == FinOrdFunction([1,2], 2)
@test k == FinOrdFunction([1,1], 2)

# Same thing with generic colimit interface
diag = Diagram([FinOrd(2),FinOrd(2),FinOrd(2)],[(1,2,f),(1,3,g)])
cocone = colimit(diag)
@test base(cocone) == FinOrd(2)
h,k = leg(cocone,2), leg(cocone,3)
@test compose(f,h) == compose(g,k)
@test h == FinOrdFunction([1,2], 2)
@test k == FinOrdFunction([1,1], 2)

end