text:
"""
Provides logging functionality for Genie apps.
"""
module Loggers

using Memento, Millboard, Dates
using Genie

import Base.log
export log

"""
    log(message, level = "info"; showst::Bool = true) :: Nothing
    log(message::Any, level::Any = "info"; showst::Bool = false) :: Nothing
    log(message::String, level::Symbol) :: Nothing

Logs `message` to all configured logs (STDOUT, FILE, etc) by delegating to `Lumberjack`.
Supported values for `level` are "info", "warn", "debug", "err" / "error", "critical".
If `level` is `error` or `critical` it will also dump the stacktrace onto STDOUT.

# Examples
```julia
```
"""
function log(message::Union{String,Symbol,Number,Exception}, level::Union{String,Symbol} = "info") :: Nothing
  message = string(message)
  level = string(level)

  if level == "err" || level == "critical"
    level = "warn"
  elseif level == "debug"
    level = "info"
  else
    level = "info"
  end

  root_logger = try
    Memento.config!(level |> string; fmt="[{date}|{level}]: {msg}")
  catch ex
    Memento.config(string(level))
  end

  try
    if isfile(log_path())
      file_logger = getlogger(@__MODULE__)
      setlevel!(file_logger, Genie.config.log_level |> string)
      push!(file_logger, DefaultHandler(log_path(), DefaultFormatter("[{date}|{level}]: {msg}")))

      Base.invoke(Core.eval(@__MODULE__, Meta.parse("Memento.$level")), Tuple{typeof(file_logger),typeof(message)}, file_logger, message)
    else
      Base.invoke(Core.eval(@__MODULE__, Meta.parse("Memento.$level")), Tuple{typeof(root_logger),typeof(message)}, root_logger, message)
    end
  catch ex
    println(string(ex))
  end

  nothing
end

"""
    truncate_logged_output(output::String) :: String

Truncates (shortens) output based on `output_length` settings and appends "..." -- to be used for limiting the output length when logging.

# Examples
```julia
julia> Genie.config.output_length
100

julia> Genie.config.output_length = 10
10

julia> Loggers.truncate_logged_output("abc " ^ 10)
"abc abc ab..."
```
"""
function truncate_logged_output(output::String) :: String
  length(output) > Genie.config.output_length && output[1:Genie.config.output_length] * "..."
end

function log_path(path = Genie.LOG_PATH) :: String
  try
    "$(joinpath(path, haskey(ENV, "GENIE_ENV") ? ENV["GENIE_ENV"] : "dev")).log"
  catch ex
    string(ex) |> println
    # println("...")
    "$(joinpath(path, "dev")).log"
  end
end
function log_path!(path = Genie.LOG_PATH) :: String
  if ! isfile(log_path(path))
    mkpath(path)
    touch(log_path(path))
  end

  log_path(path)
end

"""
    empty_log_queue() :: Vector{Tuple{String,Symbol}}

The Genie log queue is used to push log messages in the early phases of framework bootstrap,
when the logger itself is not available. Once the logger is ready, the queue is emptied and the
messages are logged.
Automatically invoked.
"""
function empty_log_queue() :: Nothing
  for log_message in Genie.GENIE_LOG_QUEUE
    log(log_message...)
  end

  empty!(Genie.GENIE_LOG_QUEUE)

  nothing
end

"""
    macro location()

Provides a macro that injects the FILE and the LINE where the logger was invoked.
"""
macro location()
  :(log(" in $(@__FILE__):$(@__LINE__)", :err))
end

# setup_loggers()
# empty_log_queue()

end

max_stars_repo_name: kiranchowdary/Genie

max_issues_repo_name: kiranchowdary/Genie

max_forks_repo_name: kiranchowdary/Genie


text:
#Framework
include("./storage.jl")
include("./optimizer.jl")
include("./shuffler.jl")

#Problem specific
include("./problems/EuclideanTSP.jl")

module Ma

using Printf
import Base
import ProgressMeter
using ..Storage
using RandomNumbers

function runBlock(code, maxIter, maxTime, momentum=0.9, updateInterval=1)
    if maxIter isa Number || maxTime isa Number
        p = ProgressMeter.Progress(1000)
    else
        p = ProgressMeter.ProgressUnknown()
    end

    startTime = time()
    estimatedThroughput = 1.0
    iterationsDone = 0
    lastProgress = 0
    currentProgress = 0.0
    while true
        if (maxTime isa Number) && time() - startTime > maxTime
            break
        elseif (maxIter isa Number ) && iterationsDone >= maxIter
            break
        end

        iterations = max(1, floor(Int, updateInterval * estimatedThroughput))

        if maxIter isa Number
            iterationsLeft = maxIter - iterationsDone
            iterations = min(iterations, iterationsLeft)
        end

        before = time()
        infos = code(iterations)
        timeTaken = time() - before

        iterationsDone += iterations
        lastThroughput = iterations / timeTaken
        estimatedThroughput = (momentum * estimatedThroughput) + (1 - momentum) * lastThroughput

        ## Reporting
        progresses = Vector{Float64}()

        if maxIter isa Number
            push!(progresses, iterationsDone / maxIter)
        end

        if maxTime isa Number
            push!(progresses, (time() - startTime) / maxTime)
        end

        if length(progresses) > 0
            currentProgress = maximum(progresses) * 1000
            if currentProgress > lastProgress
                improvement = floor(Int, currentProgress - lastProgress)
                lastProgress = currentProgress
                ProgressMeter.next!(p, step=improvement)
            end
        else
            ProgressMeter.next!(p, step=iterations)
        end

        infos = Dict(pairs(infos))
        infos[:t] = @sprintf("%i/s", (iterationsDone) / (time() - startTime))

        strings = map(collect(infos)) do (key, value)
            "$key=$value"
        end
        p.desc = "Optimizing ($(join(strings, ", "))) "
    end

    ProgressMeter.finish!(p)
    println()

    return (iterationsDone, time() - startTime)
end

mutable struct Manager
    problem
    optimizer
    shuffler
    losses::Vector{Float32}
    storage::Dict{Symbol, Any}
    debug::Bool

    function Manager()
        new(missing, missing, missing, Vector{Float32}(), Dict{Symbol, Any}(), false)
    end
end

function allocateStorage!(m::Manager)
    resize!(m.losses, m.shuffler.popSize)
    m.storage[:stateStorage] = Storage.allocateCPU(m.problem.solType,
                                                   (m.shuffler.popSize,
                                                    m.optimizer.solutionStates, ))
    m.storage[:dataStorage] = Storage.allocateCPU(m.problem.dataType, (1, ))[1]
    if !isnothing(m.optimizer.stateType)
        m.storage[:optimStorage] = Storage.allocateCPU(m.optimizer.stateType,
                                                       (m.shuffler.popSize,))
    end
end

function init!(m::Manager)
    Storage.initFromFields(m.storage[:dataStorage], m.problem.data...)

    # We only init the best state,
    for i in 1:m.shuffler.popSize
        m.problem.init(m.storage[:stateStorage][i, 1], m.storage[:dataStorage])
        if !isnothing(m.optimizer.stateType) && !isnothing(m.optimizer.init)
            m.optimizer.init(m.storage[:optimStorage][i])
        end
    end

end

function run!(m::Manager, maxIter=nothing, maxTime=nothing)

    optimStorage = nothing

    runBlock(maxIter, maxTime) do blocksIter
        Base.Threads.@threads for i in 1:m.shuffler.popSize
            rng_xor = RandomNumbers.Xorshifts.Xoroshiro128Star()
            if haskey(m.storage, :optimStorage)
                optimStorage = m.storage[:optimStorage][i]
            end
            state = m.storage[:stateStorage][i]
            loss = m.optimizer.optimize(optimStorage, state,
                                        m.storage[:dataStorage], blocksIter, rng_xor)
            m.losses[i] = loss
        end
        return (
                loss=@sprintf("%.3f - %3f", minimum(m.losses), maximum(m.losses)),
               )
    end
end

function bestSolution(m::Manager)
    bestIx = argmin(m.losses)
    return (m.losses[bestIx], m.storage[:stateStorage][bestIx][1])
end

export Manager

end

import .Ma
import .TSP

import .Storage
import .Optim
import .Shufflers

# Code generation

# data = TSP.readTourFile("./pla85900.tsp")
data = TSP.generateRandomTSPData(10000)
manager = Ma.Manager()
manager.problem = TSP.generateProblemForData(data)
manager.optimizer = Optim.ExhaustiveLocalSearch(manager.problem)
manager.shuffler = Shufflers.Independent(64)
Ma.allocateStorage!(manager)
Ma.init!(manager)
Ma.run!(manager, 10000000, 6)
loss, solution = Ma.bestSolution(manager)
print("LOSS: ", loss)

ext: jl


text:
using ITensors
using ITensorUnicodePlots

import ITensors: state, op

function state(
  ::StateName"iXm",
  ::SiteType"S=1/2")
  return [im,-im]/√2
end

function op(
  ::OpName"iX",
  ::SiteType"S=1/2")
  return [0 im; im 0]
end

# Index with site type
i = Index(2, "S=1/2")

Zp = state("Zp", i)
Zm = state("Zm", i)
Xp = state("Xp", i)
Xm = state("Xm", i)

Z = op("Z", i)
X = op("X", i)

@visualize dag(Zp) * Xp

iXm = state("iXm", i)
inner(Zm, iXm)

Z = ITensor(i', i)
Z[i'=>1, i=>1] = 1
Z[i'=>2, i=>2] = -1

z = [1 0; 0 -1]
Z = ITensor(z, i', i)
Z = op("Z", i)

Zp = state("Zp", i)
Zm = state("Zm", i)

X * Zp == Zm'
noprime(X * Zp) == Zm

(dag(Zm)' * X * Zp)[]
inner(Zm', X, Zp)

apply(X, Zp) == noprime(X * Zp)
inner(Zm, apply(X, Zp))

op("iX", i)
X * im


text:
include("../btg.jl")
include("btgBayesopt.jl")
include("../computation/finitedifference.jl")
using Plots
using Random
Random.seed!(8);
#
# trying out diff random seeds:
# failure: (3)
# success : (4, 5 works with good initialization, 6 works with init, 7 works with and w/o init.)
#use 12 for all prev experiments

#####
##### Define optimization problem, initialize GP with burn-in points, initialize btg parameters
#####

Himmelblau(x) = (x[1]^2 + x[2] -11)^2 + (x[1]+x[2]^2-7)^2 #function to optimize
y, x = sample_points(Himmelblau, [-5, -5], [5, 5]; num = 100) #burn-in points which are used as training data in GP
x = vcat(x, [5 5])
y = vcat(y, 890)
x = vcat(x, [-4.999 4.999])
y = vcat(y, 530)
x = vcat(x, [4.999 -4.999])
y = vcat(y, 610)
x = vcat(x, [-4.999 -4.999])
y = vcat(y, 249.6042059640019)
x = vcat(x, [5 2.55573])
y = vcat(y, 294.57103540062496)
x = vcat(x, [5.0 2.3475172404355718])
y = vcat(y, 279.5672977301108)

@info "y-vals", y[1:5]
### Set BTG parameters and get function handles for pdf, cdf, dpdf, etc.
Fx = linear_polynomial_basis(x)
#train = trainingData(x, Fx, y) #training data triple: location, covariates, labels
train = extensible_trainingData(x, Fx, y)
#rangeθ = reshape(select_single_theta_range(x), 1, 2) #rangetheta is auto-selected
#rangeθ = [0.0001 1000]
rangeθ = [300.0 1500]
@info "rangeθ", rangeθ
rangeλ = [-1.0 1.0]
btg1 = btg(train, rangeθ, rangeλ);
(pdf, cdf, dpdf, cdf_gradient, cdf_hessian) = solve(btg1; derivatives = true) #get function handles for pdf, cdf, derivtives

#####
##### Define acquisition function and derivatives
#####

println("defining acquisition function...")
#lx = [1, -5, -5]; ux = [3000.0, 5, 5] #box-constraints for optimization problem
lx = [1, -5, -5]; ux = [3000.0, 5, 5] #box-constraints for optimization problem
#make input an augmented vector [y, s]: label, location
cdf_fixed(v) = cdf(v[2:end], linear_polynomial_basis(v[2:end]), v[1])
cdf_gradient_fixed(v) = cdf_gradient(v[2:end], linear_polynomial_basis(v[2:end]), v[1])

d = length(lx)
fun(x...) = x[1] #x\in R^d+1
dfun(g, x...) = (g[1] = 1; g[i] = 0 for i = 2:length(g);)
#cdf_wrapper(x) = cdf_fixed(collect(x))
cdf_wrapper(x...) = cdf_fixed(collect(x))
#cdf_gradient_wrapper(g, x...) = (grad = cdf_gradient_fixed(collect(x...)); g[i] = grad[i] for i = 1:length(grad); @info grad; @info g; return nothing)
function cdf_gradient_wrapper(g, x...)
    grad = cdf_gradient_fixed(collect(x))
    @info "grad", grad
    for i = 1:length(grad)
        #@info "i", i
        #@info "grad component i", grad[i]
        global g[i] = grad[i]
    end
    return nothing
end
model = Model(Ipopt.Optimizer)
set_optimizer_attributes(model, "tol" => 1e-2, "max_iter" => 200)
#model = Model(Gurobi.Optimizer)
#NLoptSolver(algorithm=:SLSQP))
register(model, :fun, 3, fun, dfun)
register(model, :cdf_wrapper, 3, cdf_wrapper, cdf_gradient_wrapper)
#####
##### Initial guess for IPopt. Check gradient before starting optimizer
#####
#println("checking gradient at initial point...")
#####
##### Register variables and constraints in model and run optimizer.
#####
#@variable(model, lx[i] <= au[i=1:3] <= ux[i], start = initval[i])
@variable(model, lx[i] <= au[i=1:3] <= ux[i])
#@variable(model, au)
@NLobjective(model, Min, fun(au...))
@NLconstraint(model, cdf_wrapper(au...) == 0.25)

function single_optimization()
    initval = init_constrained_pt(cdf_fixed, lx, ux; quantile = 0.25)
    set_start_value(au[1], initval[1])
    set_start_value(au[2], initval[2])
    set_start_value(au[3], initval[3])
    #println("Running Optimizer...")
    #initval = [200, 4, 4]
    #@info "initval", initval
    #@variable(model, lx[i] <= au[i=1:3] <= ux[i])
    JuMP.optimize!(model)
    println(value.(au))
    vstar = value.(au)
    return (vstar, initval)
end

include("test_script.jl")
function run_func(n)
    x = [];
    vstars = [];
    init_val_cdf = []
    init_vals = []
    for i = 1:n
        Random.seed!(i+200)
        (vstar, initval) = single_optimization();
        push!(x, cdf_fixed(vstar));
        push!(vstars, vstar);
        push!(init_val_cdf, cdf_fixed(initval))
        push!(init_vals, initval)
    end
    return x, vstars, init_val_cdf, init_vals
end

if true
    (res, vstars, init_val_cdf, init_vals) = run_func(4);
    vcdfplot_sequence(vstars; res = res, upper = 2000) #shows cdfs at points we converged to
end

### check derivatives after the fact

#for i = 1:length(vstars)
#    check
#end

#(_, _, plt1, pol) = checkDerivative(cdf_fixed, cdf_gradient_fixed, vstars[1],  nothing, 4, 8)

for i = 1:10 #take 10 BO steps
    #(u_star, s_star) = optimize_acquisition(cdf, cdf_gradient, cdf_hessian)
    #update BTG trainingData and trainingBuffers with new point (xstar, Fxstar, ystar)  #location, covariates, label
    # update
    #update!(btg, s_star, cov_fun(s_star), u_star
    # - Σθ_inv_X, check  compute directly in O(n^2p)
    # - choleskyΣθ, check compute directly, have cholesky factorization anyways
    # - choleskyXΣX, check compute directly
    # - logdetΣθ    compute directly
    # - logdetXΣX compute directly
    #
end
#4.51690700362799, 3.011185163108001

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module MNIST_GPU

using CuArrays

include("getdata.jl")

export trainimages, testimages,
       trainlabels, testlabels,
       traindata,   testdata

function getimage(io::IO)
    img = Array{Int}(undef, NCOLS, NROWS)
    for i in 1:NCOLS, j in 1:NROWS
        img[i, j] = read(io, UInt8)
    end
    return CuArray(img)
end

function getlabel(io::IO)
    i = read(io, UInt8)
    v = zeros(10)
    v[i+1] = 1.
    return CuArray(v)
end

end  # module MNIST_GPU


text:
module LOBMExec

using LatentOrderBookModel

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    parsed_args = parse_commandline()
    rdpp = ReactionDiffusionPricePaths(
        1,
        parsed_args["T"],
        parsed_args["p₀"],
        parsed_args["M"],
        parsed_args["β"],
        parsed_args["L"],
        parsed_args["D"],
        parsed_args["σ"],
        parsed_args["nu"],
        SourceTerm(parsed_args["λ"], parsed_args["μ"]),
    )

    lob_densities,
    price_paths,
    mid_price_bars,
    P⁺s,
    P⁻s = rdpp(parsed_args["SEED"])
    print(mid_price_bars[:, 1])
    return 0
end

end


text:
#############################################################################
#############################################################################
#
# This file defines the Term type with several operations
#
#############################################################################
#############################################################################

##############################
# Term type and construction #
##############################

"""
A term.
"""
struct Term  #structs are immutable by default
    coeff::Int
    degree::Int
    function Term(coeff::Int, degree::Int)
        degree < 0 && error("Degree must be non-negative")
        coeff != 0 ? new(coeff,degree) : new(coeff,0)
    end
end

"""
Creates the zero term.
"""
zero(::Type{Term})::Term = Term(0,0)

"""
Creates the unit term.
"""
one(::Type{Term})::Term = Term(1,0)

###########
# Display #
###########

"""
Show a term.
"""
show(io::IO, t::Term) = print(io, "$(t.coeff)⋅x^$(t.degree)") #\cdot + [TAB]

########################
# Queries about a term #
########################

"""
Check if a term is 0.
"""
iszero(t::Term)::Bool = iszero(t.coeff)

"""
Compare two terms.
"""
isless(t1::Term,t2::Term)::Bool =  t1.degree == t2.degree ? (t1.coeff < t2.coeff) : (t1.degree < t2.degree)

"""
Evaluate a term at a point x.
"""
evaluate(t::Term, x::T) where T <: Number = t.coeff * x^t.degree

##########################
# Operations with a term #
##########################

"""
Add two terms of the same degree.
"""
function +(t1::Term,t2::Term)::Term
    @assert t1.degree == t2.degree
    Term(t1.coeff + t2.coeff, t1.degree)
end

"""
Negate a term.
"""
-(t::Term,) = Term(-t.coeff,t.degree)

"""
Subtract two terms with the same degree.
"""
-(t1::Term, t2::Term)::Term = t1 + (-t2)

"""
Multiply two terms.
"""
*(t1::Term, t2::Term)::Term = Term(t1.coeff * t2.coeff, t1.degree + t2.degree)

"""
Compute the symmetric mod of a term with an integer.
"""
mod(t::Term, p::Int) = Term(mod(t.coeff,p), t.degree)

"""
Compute the derivative of a term.
"""
derivative(t::Term) = Term(t.coeff*t.degree,max(t.degree-1,0))

"""
Divide two terms. Returns a function of an integer.
"""
function ÷(t1::Term,t2::Term) #\div + [TAB]
    @assert t1.degree ≥ t2.degree
    f(p::Int)::Term = Term(mod((t1.coeff * int_inverse_mod(t2.coeff, p)), p), t1.degree - t2.degree)
end

"""
Integer divide a term by an integer.
"""
÷(t::Term, n::Int) = t ÷ Term(n,0)


text:
@doc raw"""
    LogCholeskyMetric <: Metric

The Log-Cholesky metric imposes a metric based on the Cholesky decomposition as
introduced by [^Lin2019].

[^Lin2019]:
    > Lin, Zenhua: "Riemannian Geometry of Symmetric Positive Definite Matrices via
    > Cholesky Decomposition", arXiv: [1908.09326](https://arxiv.org/abs/1908.09326).
"""
struct LogCholeskyMetric <: RiemannianMetric end

cholesky_to_spd(x, W) = (x * x', W * x' + x * W')

tangent_cholesky_to_tangent_spd!(x, W) = (W .= W * x' + x * W')

spd_to_cholesky(p, X) = spd_to_cholesky(p, cholesky(p).L, X)

function spd_to_cholesky(p, x, X)
    w = inv(x) * X * inv(transpose(x))
    # strictly lower triangular plus half diagonal
    return (x, x * (LowerTriangular(w) - Diagonal(w) / 2))
end

@doc raw"""
    distance(M::MetricManifold{SymmetricPositiveDefinite,LogCholeskyMetric}, p, q)

Compute the distance on the manifold of [`SymmetricPositiveDefinite`](@ref)
nmatrices, i.e. between two symmetric positive definite matrices `p` and `q`
with respect to the [`LogCholeskyMetric`](@ref). The formula reads

````math
d_{\mathcal P(n)}(p,q) = \sqrt{
 \lVert ⌊ x ⌋ - ⌊ y ⌋ \rVert_{\mathrm{F}}^2
 + \lVert \log(\operatorname{diag}(x)) - \log(\operatorname{diag}(y))\rVert_{\mathrm{F}}^2 }\ \ ,
````

where $x$ and $y$ are the cholesky factors of $p$ and $q$, respectively,
$⌊\cdot⌋$ denbotes the strictly lower triangular matrix of its argument,
and $\lVert\cdot\rVert_{\mathrm{F}}$ the Frobenius norm.
"""
function distance(
    M::MetricManifold{ℝ,SymmetricPositiveDefinite{N},LogCholeskyMetric},
    p,
    q,
) where {N}
    return distance(CholeskySpace{N}(), cholesky(p).L, cholesky(q).L)
end

@doc raw"""
    exp(M::MetricManifold{SymmetricPositiveDefinite,LogCholeskyMetric}, p, X)

Compute the exponential map on the [`SymmetricPositiveDefinite`](@ref) `M` with
[`LogCholeskyMetric`](@ref) from `p` into direction `X`. The formula reads

````math
\exp_p X = (\exp_y W)(\exp_y W)^\mathrm{T}
````

where $\exp_xW$ is the exponential map on [`CholeskySpace`](@ref), $y$ is the cholesky
decomposition of $p$, $W = y(y^{-1}Xy^{-\mathrm{T}})_\frac{1}{2}$,
and $(\cdot)_\frac{1}{2}$
denotes the lower triangular matrix with the diagonal multiplied by $\frac{1}{2}$.
"""
exp(::MetricManifold{ℝ,SymmetricPositiveDefinite,LogCholeskyMetric}, ::Any...)

function exp!(
    M::MetricManifold{ℝ,SymmetricPositiveDefinite{N},LogCholeskyMetric},
    q,
    p,
    X,
) where {N}
    (y, W) = spd_to_cholesky(p, X)
    z = exp(CholeskySpace{N}(), y, W)
    return copyto!(q, z * z')
end

@doc raw"""
    inner(M::MetricManifold{LogCholeskyMetric,ℝ,SymmetricPositiveDefinite}, p, X, Y)

Compute the inner product of two matrices `X`, `Y` in the tangent space of `p`
on the [`SymmetricPositiveDefinite`](@ref) manifold `M`, as
a [`MetricManifold`](@ref) with [`LogCholeskyMetric`](@ref). The formula reads

````math
    g_p(X,Y) = ⟨a_z(X),a_z(Y)⟩_z,
````

where $⟨\cdot,\cdot⟩_x$ denotes inner product on the [`CholeskySpace`](@ref),
$z$ is the cholesky factor of $p$,
$a_z(W) = z (z^{-1}Wz^{-\mathrm{T}})_{\frac{1}{2}}$, and $(\cdot)_\frac{1}{2}$
denotes the lower triangular matrix with the diagonal multiplied by $\frac{1}{2}$
"""
function inner(
    M::MetricManifold{ℝ,SymmetricPositiveDefinite{N},LogCholeskyMetric},
    p,
    X,
    Y,
) where {N}
    (z, Xz) = spd_to_cholesky(p, X)
    (z, Yz) = spd_to_cholesky(p, z, Y)
    return inner(CholeskySpace{N}(), z, Xz, Yz)
end

@doc raw"""
    log(M::MetricManifold{SymmetricPositiveDefinite,LogCholeskyMetric}, p, q)

Compute the logarithmic map on [`SymmetricPositiveDefinite`](@ref) `M` with
respect to the [`LogCholeskyMetric`](@ref) emanating from `p` to `q`.
The formula can be adapted from the [`CholeskySpace`](@ref) as
````math
\log_p q = xW^{\mathrm{T}} + Wx^{\mathrm{T}},
````
where $x$ is the cholesky factor of $p$ and $W=\log_x y$ for $y$ the cholesky factor
of $q$ and the just mentioned logarithmic map is the one on [`CholeskySpace`](@ref).
"""
log(::MetricManifold{ℝ,SymmetricPositiveDefinite,LogCholeskyMetric}, ::Any...)

function log!(
    M::MetricManifold{ℝ,SymmetricPositiveDefinite{N},LogCholeskyMetric},
    X,
    p,
    q,
) where {N}
    x = cholesky(p).L
    y = cholesky(q).L
    log!(CholeskySpace{N}(), X, x, y)
    return tangent_cholesky_to_tangent_spd!(x, X)
end

@doc raw"""
    vector_transport_to(
        M::MetricManifold{SymmetricPositiveDefinite,LogCholeskyMetric},
        p,
        X,
        q,
        ::ParallelTransport,
    )

Parallel transport the tangent vector `X` at `p` along the geodesic to `q` with respect to
the [`SymmetricPositiveDefinite`](@ref) manifold `M` and [`LogCholeskyMetric`](@ref).
The parallel transport is based on the parallel transport on [`CholeskySpace`](@ref):
Let $x$ and $y$ denote the cholesky factors of `p` and `q`, respectively and
$W = x(x^{-1}Xx^{-\mathrm{T}})_\frac{1}{2}$, where $(\cdot)_\frac{1}{2}$ denotes the lower
triangular matrix with the diagonal multiplied by $\frac{1}{2}$. With $V$ the parallel
transport on [`CholeskySpace`](@ref) from $x$ to $y$. The formula hear reads

````math
\mathcal P_{q←p}X = yV^{\mathrm{T}} + Vy^{\mathrm{T}}.
````
"""
vector_transport_to(
    ::MetricManifold{ℝ,SymmetricPositiveDefinite,LogCholeskyMetric},
    ::Any,
    ::Any,
    ::Any,
    ::ParallelTransport,
)

function vector_transport_to!(
    M::MetricManifold{ℝ,SymmetricPositiveDefinite{N},LogCholeskyMetric},
    Y,
    p,
    X,
    q,
    m::ParallelTransport,
) where {N}
    y = cholesky(q).L
    (x, W) = spd_to_cholesky(p, X)
    vector_transport_to!(CholeskySpace{N}(), Y, x, W, y, m)
    return tangent_cholesky_to_tangent_spd!(y, Y)
end


text:
#OnDieDecap.jl: Package to generate plots for on-die decoupling models
#-------------------------------------------------------------------------------

module OnDieDecap

using Colors
import DelimitedFiles: readdlm
using MDDatasets
import CMDimData
using CMDimData.EasyPlot #set, cons
import CMDimCircuits
import InspectDR
CMDimData.@includepkg EasyPlotInspect

const j = im
include("plotbuilders.jl")
include("plots/supplyimpedance.jl")
include("plots/supplybounce.jl")

run_app() = nothing

end # module

#Last line


text:
export complexity, depth, projectToPositiveSimplex!, initllhvals

"""
Get all nodes in topological order using Tarjan's algoritm.
"""
function getOrderedNodes(root)
  visitedNodes = Vector{SPNNode}()
  visitNode!(root, visitedNodes)
  return visitedNodes
end

function visitNode!(node::Node, visitedNodes)
  # check if we have already visited this node
  if !(node in visitedNodes)

    # visit node
    for child in children(node)
      visitNode!(child, visitedNodes)
    end

    push!(visitedNodes, node)
  end
end
visitNode!(node::Leaf, visitedNodes) = push!(visitedNodes, node)

"""

depth(S)

Compute the depth of the SPN rooted at S.
"""
depth(S::Node) = maximum(ndepth(child, 1) for child in children(S))
depth(S::Leaf) = 0

ndepth(S::Node, d::Int) = maximum(ndepth(child, d+1) for child in children(S))
ndepth(S::Leaf, d::Int) = d

"""

    complexity(S)

Compute the complexity (number of free parameters) of the SPN rooted at S.
"""
complexity(spn::SumProductNetwork) = complexity(spn.root)
complexity(node::SumNode) = mapreduce(c -> complexity(c), +, children(node)) + length(node)
complexity(node::ProductNode) = mapreduce(c -> complexity(c), +, children(node))
complexity(node::IndicatorNode) = 0
complexity(node::UnivariateNode) = length(fieldnames(node.d))
complexity(node::MultivariateNode) = length(fieldnames(node.d))

"""
  sub2ind2(size1, ind1, ind2) -> linear index
"""
sub2ind2(s1, i, j) = i + (j-1)*s1

"""
  sub2ind3(size1, size2, ind1, ind2, ind3) -> linear index
"""
sub2ind3(s1, s2, i, j, k) = i + (j-1)*s1 + (k-1)*s1*s2

"""

    projectToPositiveSimplex!(q::AbstractVector{<:Real}; lowerBound = 0.0, s = 1.0)

Project q to the positive simplex to ensure sum(q) == s.

##### Details
See Algorithm 1 in:
    Duchi, J., Shalev-Shwartz, S., Singer, Y., and Chandra, T.: Efficient projections onto the L 1-ball for learning in high dimensions. In proceeding of ICML 2008

"""
function projectToPositiveSimplex!(q::AbstractVector{<:Real}; lowerBound = 0.0, s = 1.0)

    if sum(q) == 0.0
        q[:] = ones(length(q)) / length(q)
        return q
    end

    if (sum(q) == s) & all(q .> lowerBound)
        return q
    end

    N = length(q)

    U = sort(q, rev=true)
    CSU = cumsum(U)
    CSUU = U .* collect(1:N) .>= (CSU .- s)
    ρ = maximum(findall(CSUU))
    θ = (CSU[ρ] - s) / ρ

    q[:] .-= θ
    q[q .< lowerBound] .= lowerBound
    return q
end

"""
    initllhvals(spn::SumProductNetwork, X::AbstractMatrix)

Construct a log likelihoods data-structure using `spn` and `X`.
"""
function initllhvals(spn::SumProductNetwork, X::AbstractMatrix)
    return initllhvals(spn, size(X,1))
end

function initllhvals(spn::SumProductNetwork, X::AbstractVector)
    return initllhvals(spn, 1)
end

function initllhvals(spn::SumProductNetwork, N::Int)
    idx = Axis{:id}(collect(keys(spn)))
    if N == 1
        return AxisArray(ones(length(idx)) * -Inf, idx)
    else
        return AxisArray(ones(N, length(idx)) * -Inf, 1:N, idx)
    end
end


text:
"""
	module Boundaries
"""
module Boundaries

export Boundary

using ...Defaults
using ..BoundaryConditions
using ..BoundaryLayers
using ..Shapes

"""
	struct Boundary

	Boundary(shape::Shape,bcs,bls) -> bnd

`bcs` is either tuple of boundary conditions, or it is a single boundary condition which applies to all sides.
`bls` is either tuple of boundary layers, or it is a single boundary layer which applies to all sides.

if fewer boundary conditions or layers are supplied than the number of sides of
`shape`, it assumes `noBC` and `noBL`
"""
struct Boundary{TS,TBC,TBL}
    shape::TS
    bcs::TBC
    bls::TBL

    function Boundary(shape::AbstractShape{NDIMS,NSIDES},
        bcs::NTuple{NSIDES,AbstractBC},
        bls::NTuple{NSIDES,AbstractBL}
        ) where {NDIMS,NSIDES}

		bls = map(x->_bls_by_shape(x,shape),bls)
	    return new{typeof(shape),typeof(bcs),typeof(bls)}(shape,bcs,bls)
    end
end

Base.ndims(bnd::Boundary) = ndims(bnd.shape)
Shapes.nsides(bnd::Boundary) = nsides(bnd.shape)

# function Boundary(shape::AbstractShape{NS},bc::AbstractBC=noBC(),bls::AbstractBL=noBL()) where N
	# bcs = fill(bcs,N)
	# bls = fill(bls,N)
	# return Boundary(shape,(bcs...,),(bls...,))
# end
# function Boundary(shape::AbstractShape{N},bcs::AbstractBC,bls) where N
# 	bcs = fill(bcs,N)
# 	return Boundary(shape,(bcs...,),bls)
# end
# function Boundary(shape::AbstractShape{N},bcs,bls::AbstractBL=noBL()) where N
# 	bls = fill(bls,N)
# 	return Boundary(shape,bcs,(bls...,))
# end

function Base.show(io::IO,bnd::Boundary)
	N = nsides(bnd)
	print(io,"Boundary: \n")
	print(io,"\tshape: $(bnd.shape)")
	for i ∈ 1:N
		println(io)
		print(io,"\t", bnd.bcs[i],"  +  ",bnd.bls[i])
	end
end

"""
	conj(::Boundary) -> bnd

Make a new `Boundary` object with conjugated boundary layers.
"""
Base.conj(bnd::Boundary) = Boundary(bnd.shape,conj.(bnd.bcs),conj.(bnd.bls))

_bls_by_shape(bls,shape::AbstractShape) = bls

include("1D/Boundaries.jl")
# include("2D/Boundaries.jl")
# include("3D/Boundaries.jl")

end # module


text:
using Turing, Random, Test
using Turing: Sampler

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@testset "sghmc.jl" begin
    @numerical_testset "sghmc inference" begin
        Random.seed!(125)

        alg = SGHMC(10000, 0.02, 0.5)
        chain = sample(gdemo_default, alg)
        check_gdemo(chain, atol = 0.1)
    end
    @turing_testset "sghmc constructor" begin
        alg = SGHMC(1000, 0.01, 0.1)
        sampler = Sampler(alg)
        @test isa(alg, SGHMC)
        @test isa(sampler, Sampler{<:Turing.SGHMC})

        alg = SGHMC(200, 0.01, 0.1, :m)
        sampler = Sampler(alg)
        @test isa(alg, SGHMC)
        @test isa(sampler, Sampler{<:Turing.SGHMC})

        alg = SGHMC(1000, 0.01, 0.1, :s)
        sampler = Sampler(alg)
        @test isa(alg, SGHMC)
        @test isa(sampler, Sampler{<:Turing.SGHMC})
    end
end

@testset "sgld.jl" begin
    @turing_testset "sgld constructor" begin
        alg = SGLD(1000, 0.25)
        sampler = Sampler(alg)
        @test isa(alg, SGLD)
        @test isa(sampler, Sampler{<:Turing.SGLD})

        alg = SGLD(200, 0.25, :m)
        sampler = Sampler(alg)
        @test isa(alg, SGLD)
        @test isa(sampler, Sampler{<:Turing.SGLD})

        alg = SGLD(1000, 0.25, :s)
        sampler = Sampler(alg)
        @test isa(alg, SGLD)
        @test isa(sampler, Sampler{<:Turing.SGLD})
    end
    @numerical_testset "sgld inference" begin
        Random.seed!(125)

        chain = sample(gdemo_default, SGLD(10000, 0.5))

        # Note: samples are weighted by step sizes cf 4.2 in paper
        v = get(chain, [:lf_eps, :s, :m])
        s_res1weightedMean = sum(v.lf_eps .* v.s) / sum(v.lf_eps)
        m_res1weightedMean = sum(v.lf_eps .* v.m) / sum(v.lf_eps)
        @test s_res1weightedMean ≈ 49/24 atol=0.2
        @test m_res1weightedMean ≈ 7/6 atol=0.2
    end
end


text:
using ExpectationStubs
using ExpectationStubs: ExpectationValueMismatchError
using Test
using Test: @test_logs

@testset "basic stubbing" begin
    @stub foo
    @expect(foo(::Any)=32)
    @test foo(3)==32
    @test foo(4)==32

    @test_logs (:warn, "Expectation already set") @expect(foo(::Any)=34)
    @test foo(5)==34
end

@testset "multitype stubbing" begin
    @stub foo
    @expect(foo(::Int)=20)
    @expect(foo(::Bool)=39)
    @test foo(3)==20
    @test foo(false)==39
end

@testset "mulitple basic keyed stubs" begin
    @stub foo
    @expect(foo(1)=37)
    @expect(foo(2)=38)
    @test foo(2)==38
    @test foo(1)==37
end

@testset "multi-arg stubbing" begin
    @stub foo
    @expect(foo(4, ::Int)=32)
    @test foo(4,3)==32
    @test_throws ExpectationValueMismatchError foo(5, 3)==32
    @test_throws ExpectationValueMismatchError foo(3, 2.0)==32
end

@testset "mixed keyed stubs" begin
    @stub foo
    @expect(foo(1)=30)
    # Throws errors
    #@expect(foo(::Any)=370)
    @test_broken foo(2)==370
    @test foo(1)==30
end

@testset "All expectations used" begin
   @testset "value" begin
        @stub foo
        @expect(foo(1)=30)
        @expect(foo(2)=35)

        foo(1)
        @test !all_expectations_used(foo)

        foo(2)
        @test all_expectations_used(foo)
    end
    @testset "typed" begin
        @stub foo
        @expect(foo(::Int)=30)

        @test !all_expectations_used(foo)
        foo(2)
        @test all_expectations_used(foo)
    end

end

@testset "Expectations use count" begin

    @testset "basic" begin
        @stub foo
        @expect(foo(1)=30)
        @expect(foo(2)=35)
        @expect(foo(3)=40)

        @test @usecount(foo(1)) == 0
        foo(1)
        @test @usecount(foo(1)) == 1
        @test @usecount(foo(2)) == 0
        foo(1)

        @test @usecount(foo(1)) == 2

        @testset "total" begin
            foo(3)
            @test @usecount(foo(::Any)) == 3
            @test @usecount(foo(::Int)) == 3
        end
    end
end

@testset "Expectations used" begin

    @testset "basic" begin
        @stub foo
        @expect(foo(1)=30)
        @expect(foo(2)=35)

        @test foo(1) == 30
        @test @used(foo(1))
        @test !@used(foo(2))
    end

    @testset "mixed" begin
        @stub foo
        @expect(foo(::Int)=310)

        @test foo(105) == 310
        @test @used(foo(105))
        @test @used(foo(::Int))
        @test !@used(foo(2))
        @test !@used(foo(::Bool))
    end

    @testset "using a variable as the value #11" begin
        @stub bar

        b = 11
        @expect bar(b)=51

        @test bar(b) == 51
        @test @used bar(11)
        @test @used bar(b)
    end
end

########################################
include("donotcare.jl")


text:
using Libdl

function find_lib()
    MPI_LIBRARY_PATH = get(ENV, "JULIA_MPI_LIBRARY_PATH") do
        MPI_PATH = get(ENV, "JULIA_MPI_PATH", nothing)
        if MPI_PATH !== nothing
            joinpath(MPI_PATH,"lib")
        else
            nothing
        end
    end

    libmpi = get(ENV, "JULIA_MPI_LIBRARY") do
        find_library(["libmpi", "libmpi_ibm", "msmpi", "libmpich"],
                     MPI_LIBRARY_PATH !== nothing ? [MPI_LIBRARY_PATH] : [])
    end

    if libmpi == ""
        error("No MPI library found.\nEnsure an MPI implementation is loaded, or set the `JULIA_MPI_PATH` variable.")
    end

    return libmpi
end

const libmpi = find_lib()

function find_mpiexec()
    get(ENV, "JULIA_MPIEXEC") do
        MPI_PATH = get(ENV, "JULIA_MPI_PATH", nothing)
        if MPI_PATH !== nothing && Sys.isexecutable(joinpath(MPI_PATH,"bin","mpiexec"))
            joinpath(MPI_PATH,"bin","mpiexec")
        else
            "mpiexec"
        end
    end
end

const mpiexec_path = find_mpiexec()

"""
    mpiexec(fn)

A wrapper function for the MPI launcher executable. Calls `fn(cmd)`, where `cmd` is a `Cmd` object of the MPI launcher.

# Environment Variables

The behaviour of `mpiexec` can be controlled by the following environment variables:

- `JULIA_MPIEXEC`: the name or path of the launcher executable (set at compile time).
- `JULIA_MPIEXEC_ARGS`: additional arguments that are passed to the launcher. These are space seperated, supporting the same quoting rules as
  Julia `Cmd` objects. These can be modified at run time.

# Usage

```jldoctest
julia> mpiexec(cmd -> run(`\$cmd -n 3 echo hello world`));
hello world
hello world
hello world
```
"""
function mpiexec(fn)
    mpiexec_args = Base.shell_split(get(ENV, "JULIA_MPIEXEC_ARGS", ""))
    fn(`$mpiexec_path $mpiexec_args`)
end

const use_stdcall = startswith(basename(libmpi), "msmpi")

macro mpicall(expr)
    @assert expr isa Expr && expr.head == :call && expr.args[1] == :ccall
    # Microsoft MPI uses stdcall calling convention
    # this only affects 32-bit Windows
    # unfortunately we need to use ccall to call Get_library_version
    # so check using library name instead
    if use_stdcall
        insert!(expr.args, 3, :stdcall)
    end
    return esc(expr)
end


text:
function main()
    doubles = [2 * n for n in 1:20]
    push!(doubles, 50)

    checkout_one = Dict()
    for n in doubles
        checkout_one[n] = 1
    end

    score_one = Dict()
    for n in 1:20
        score_one[n] = get(score_one, n, 0) + 1
        score_one[2 * n] = get(score_one, 2 * n, 0) + 1
        score_one[3 * n] = get(score_one, 3 * n, 0) + 1
    end
    score_one[25] = 1
    score_one[50] = 1

    checkout_two = Dict()
    for n in 1:50
        for m in 1:60
            checkout_two[n + m] = get(checkout_two, n + m, 0) + get(checkout_one, n, 0) * get(score_one, m, 0)
        end
    end

    score_two = Dict()
    for n in 1:60
        for m in 1:n
            if n == m
                score_two[n + m] = get(score_two, n + m, 0) + 0.5 * get(score_one, n, 0) * (get(score_one, m, 0) + 1)
            else
                score_two[n + m] = get(score_two, n + m, 0) + get(score_one, n, 0) * get(score_one, m, 0)
            end
        end
    end

    checkout_three = Dict()
    for n in 1:50
        for m in 1:120
            checkout_three[n + m] = get(checkout_three, n + m, 0) + get(checkout_one, n, 0) * get(score_two, m, 0)
        end
    end

    count = 0
    for n in 1:99
        count += get(checkout_one, n, 0)
        count += get(checkout_two, n, 0)
        count += get(checkout_three, n, 0)
    end

    println(Int(count))
end

main()

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
function enforceadhesivity(h::EntropyCone{N, ET}) where {N, ET<:Real}
    lift = h
    last = 1
    for S in indset(h)
        for T in indset(h)
            if setdiff(S, T) == S # no intersection
                I = setdiff(setdiff(ntodim(h.n), S), T)
                cur = copy(h)
                intersect!(cur, submodulareq(cur.n, S, T, I))
                lift *= cur
                last += 1
                equalonsubsetsof!(lift, 1, last, union(S, I))
                equalonsubsetsof!(lift, 1, last, union(T, I))
            end
        end
    end
    unlift(lift)
end

export enforceadhesivity


text:
Base.getindex(a::IndicesArray{T,N}, i::Colon) where {T,N} = a

@propagate_inbounds function Base.getindex(a::IndicesVector, i)
    return _unsafe_getindex(a, parent(a), to_index(axes(a, 1), i))
end

@propagate_inbounds function Base.getindex(a::IndicesArray{T,N}, i...) where {T,N}
    return _unsafe_getindex(a, parent(a), to_indices(a, i))
end

function _unsafe_getindex(::T, a::AbstractArray, inds::Integer) where {T}
    @inbounds(getindex(a, inds))
end

function _unsafe_getindex(::T, a::AbstractArray, inds::AbstractIndex) where {T}
    return rebuild(T, @inbounds(getindex(a, inds)), inds)
end

function _unsafe_getindex(::T, a::AbstractArray, inds::Tuple{Vararg{Int}}) where {T}
    return @inbounds(getindex(a, inds...))
end

function _unsafe_getindex(::T, a::AbstractArray, inds::Tuple) where {T}
    return rebuild(T, @inbounds(getindex(a, values.(inds)...)), _drop_empty(inds))
end

_drop_empty(x::Tuple{Any,Vararg}) = (first(x), _drop_empty(tail(x))...)
_drop_empty(x::Tuple{<:Integer,Vararg}) = _drop_empty(tail(x))
_drop_empty(x::Tuple{}) = ()


text:
# compute models in parallel: multiprocess

addprocs(Sys.CPU_CORES) # add a worker process per core
print_with_color(:white, "Setup:\n")
println("  > Using $(nprocs()-1) worker processes")
n_core = nprocs()-1#2

if nworkers() < n_core
    addprocs(n_core-nworkers(); exeflags="--check-bounds=yes")
end
@assert nprocs() > n_core
@assert nworkers() >= n_core

println(workers())

@everywhere using PBupsModel
@everywhere using GeneralUtils
@everywhere using Base.Test
@everywhere using MAT

# read test data file
mpath = "data"
fname = "testdata.mat"
ratdata, ntrials = LoadData(mpath, fname)

dt = 0.02

# using TrialData load trial data
RightClickTimes, LeftClickTimes, maxT, rat_choice = TrialData(ratdata["rawdata"], 1)

Nsteps = Int(ceil(maxT/dt))

# sigma_a = rand()*4.; sigma_s = rand()*4.; sigma_i = rand()*30.;
# lam = randn()*0.4; B = rand()*12.+5.; bias = randn();
# phi = rand()*1.39+0.01; tau_phi = 0.695*rand()+0.005; lapse = rand();

# known parameter set
sigma_a = 1.; sigma_s = 0.1; sigma_i = 0.2;
lam = -0.5; B = 6.1; bias = 0.1;
phi = 0.3; tau_phi = 0.1; lapse = 0.05*2;

params = [lam, sigma_a, sigma_s, sigma_i, B, phi, tau_phi, bias, lapse]

# LL = LogLikelihood(params, RightClickTimes, LeftClickTimes, Nsteps, rat_choice)

# # write your own tests here
# @test (LL - -0.9972) < 0.0001

# # ###### test2 : Compute Logliklihood of 40 trials
# ntrials = 40
# LLs = SharedArray(Float64, ntrials)
# LL2 = ComputeLL(LLs, params, ratdata["rawdata"], ntrials)

# @test (LL2 - 9.0591) < 0.0001

# # # ###### test3 : Compute Gradients of 40 trials
# LL, LLgrad = ComputeGrad(params, ratdata["rawdata"], ntrials)
# print(LLgrad)

# # ###### test4 : Compute Hessian Matrix of 40 trials
# LL, LLgrad, LLhess = ComputeHess(params, ratdata["rawdata"], ntrials)
# print(LLhess)

# ###### test5 : Model Optimization
# init_params = InitParams()
# result = ModelFitting(init_params, ratdata, ntrials)
# FitSummary(mpath, fname, result)

# ====================== new test ====================== #

# to generalize the model..
# any subset of parameters in 12-param model.
# change the model_likelihood.jl

# 9p / 12p (+3 bias params)
# 8p / 11p sigma_i = 0
# 10p (without adaptation parameters)

# input should be "set of name of parameters", "and their values"
# where shall we handle the sigma_s_L .. sigma_s
# if someone specify sigma_s_L/sigma_s_R,then use them
# if someone does not specify both sigma_s_L and sigma_s_R, then use one of them for sigma_s
# if someone does not specify none of sigma_s_L and sigma_s_R, then use default value for sigma_s

# how do we check that each parameter is set or not.
# compare it with default value... for now

ntrials = 40
LLs = SharedArray(Float64, ntrials)

# ===== 9p ===== #

# # known parameter set
args = ["sigma_a","sigma_s_R","sigma_i","lambda","B","bias","phi","tau_phi","lapse_R"]
x = [1., 0.1, 0.2, -0.5, 6.1, 0.1, 0.3, 0.1, 0.05*2]

LL1 = LogLikelihood(RightClickTimes, LeftClickTimes, Nsteps, rat_choice
                ;make_dict(args, x)...)
print(LL1)

LL2 = ComputeLL(LLs, ratdata["rawdata"], ntrials
    ;make_dict(args, x)...)
# LL2 = ComputeLL(LLs, ratdata["rawdata"], ntrials, args, x)

# @test (LL2 - 9.0591) < 0.0001

# # ###### test3 : Compute Gradients of 40 trials
# LL, LLgrad = ComputeGrad(ratdata["rawdata"], ntrials, args, x)
# print(LLgrad)

# # ###### test4 : Compute Hessian Matrix of 40 trials
# LL, LLgrad, LLhess = ComputeHess(ratdata["rawdata"], ntrials, args, x)
# print(LLhess)

# # ###### test5 : Model Optimization
# args = ["sigma_a","sigma_s_R","sigma_i","lambda","B","bias","phi","tau_phi","lapse_R"]
# init_params = InitParams(args)
# result = ModelFitting(args, init_params, ratdata, ntrials)
# FitSummary(mpath, fname, result)

# ===== 12p ===== #
args_12p = ["sigma_a","sigma_s_R","sigma_s_L","sigma_i","lambda","B","bias","phi","tau_phi","lapse_R","lapse_L","input_gain_weight"]
x_12p = [1., 0.1, 50, 0.2, -0.5, 6.1, 0.1, 0.3, 0.1, 0.05*2, 0.2, 0.4]

ntrials = 40
LLs = SharedArray(Float64, ntrials)
# LL = ComputeLL(LLs, ratdata["rawdata"], ntrials, args_12p, x_12p)
LL = ComputeLL(LLs, ratdata["rawdata"], ntrials
	;make_dict(args_12p, x_12p)...)
print(LL)

# ###### test3 : Compute Gradients of 40 trials
LL, LLgrad = ComputeGrad(ratdata["rawdata"], ntrials, args_12p, x_12p)
print(LLgrad)

###### test4 : Compute Hessian Matrix of 40 trials
LL, LLgrad, LLhess = ComputeHess(ratdata["rawdata"], ntrials, args_12p, x_12p)
print(LLhess)

###### test5 : Model Optimization
init_params = InitParams(args_12p)
result = ModelFitting(args_12p, init_params, ratdata, ntrials)
FitSummary(mpath, fname, result)


text:
const COL_SAMPLER_NAME = :sampler_name
const COL_NUM_METAINFERENCE = :num_metainference
const COL_NUM_PARTICLES = :num_particles
const COL_NUM_REJUVENATION = :num_rejuvenation
const COL_NUM_MCMC_ITERS = :num_mcmc_iters
const COL_NUM_SGD_ITERS = :num_sgd_iters
const COL_AIDE_ESTIMATE = :aide_estimate
const COL_AIDE_STDERR = :aide_estimate_standard_error
const COL_MEDIAN_RUNTIME = :median_runtime


text:
module TestPrettyExpr

using Test
using InitialValues: prettyexpr, def_impl, disambiguate_impl

parseall(code) = Meta.parse("begin $code end")

roundtrip(code::AbstractString) = prettyexpr(parseall(code))
roundtrip(ex::Expr) = parseall(prettyexpr(ex))

isalnn(x) = x isa LineNumberNode

normalize(x) = x
function normalize(x::Expr)
    args = normalize.(filter(!isalnn, x.args))
    if x.head == :block && length(args) == 1
        return args[1]
    end
    return Expr(x.head, args...)
end

eq(x, y) = normalize(x) == normalize(y)

@testset "$label" for (label, ex) in [
    (
        label = "f(x) = x",
        ex = :(f(x) = x),
    ),
    (
        label = "def_impl",
        ex = def_impl(:op, :x, :y),
    ),
    (
        label = "disambiguate_impl",
        ex = disambiguate_impl(:op, :RightType, :x),
    ),
]
    code = prettyexpr(ex)
    @test roundtrip(code) === code
    @test eq(roundtrip(ex), ex)
end

end  # module


text:
using LinearAlgebra
using Polynomials

using SemiLagrangian:
    getbspline, B_SplineLU, get_kl_ku, B_SplineFFT, interpolate!, decal

function test_spline(order, prec)
    setprecision(prec) do
        @time @testset "test spline object of order $order" begin
            p0 = getbspline(order, 0)
            p3 = getbspline(order, 3)
            @test p3 == decal(p0, 3)
            @test isapprox(sum(p0.(range(big"1.0", length = order + 1))), 1, atol = 1e-60)
            @test sum(p0.(range(one(Rational{BigInt}), length = order + 1))) == 1
            for i = 0:order
                x = rand(BigFloat)
                res0 = p0(x + i)
                res3 = p3(x + i + 3)
                @test isapprox(res0, res3, atol = 1e-60)
                @test order == Polynomials.degree(p0[i])
                pol1 = p0[i]
                pol2 = p0[i+1]
                for j = 1:order
                    @test pol1(i + 1) == pol2(i + 1)
                    pol1 = derivative(pol1)
                    pol2 = derivative(pol2)
                end
            end
        end
    end
end
for order = 1:15
    test_spline(order, order < 20 ? 256 : 512)
end
function nb_diff_tol(f1, f2, tol)
    nb = 0
    for i = 1:size(f1, 1)
        if abs(f1[i] - f2[i]) > tol
            #          println("indice=$i diff=$(abs(f1[i]-f2[i]))")
            nb += 1
        end
    end
    return nb
end
"""
    bsplinerec(p, j, x)

Return the value at x in [0,1[ of the B-spline with
integer nodes of degree p with support starting at j.
Implemented recursively using the de Boor's recursion formula
using the [De Boor's Algorithm](https://en.wikipedia.org/wiki/De_Boor%27s_algorithm)

Note that this function has an exponential complexity depending to the order, so it is only there for testing.

```math
B_{i,0}(x) := \\left\\{
\\begin{matrix}
1 & \\mathrm{if}  \\quad t_i ≤ x < t_{i+1} \\\\
0 & \\mathrm{otherwise}
\\end{matrix}
\\right.
```

```math
B_{i,p}(x) := \\frac{x - t_i}{t_{i+p} - t_i} B_{i,p-1}(x)
+ \\frac{t_{i+p+1} - x}{t_{i+p+1} - t_{i+1}} B_{i+1,p-1}(x).
```

"""
function bsplinerec(p, j, x::T) where {T}

    if p == 0
        if j == 0
            return one(T)
        else
            return zero(T)
        end
    else
        w = (x - j) / p
        w1 = (x - j - 1) / p
    end
    (w * bsplinerec(p - 1, j, x) + (1 - w1) * bsplinerec(p - 1, j + 1, x))

end

function test_bspline()

    @time @testset "test verify bspline" begin
        tab_x = [
            1821 // 10000,
            1 // 1234,
            7677 // 8999,
            big"123456789" // big"234567890",
            0 // 1,
        ]
        for ord = 1:10, x in tab_x
            if ord >8
                p1 = Polynomial([big(1//7),2//7,3//7])
                p2 = Polynomial([-big(6//7), 2//7, 3//7])
            else
                p1 = Polynomial([1//7,2//7,3//7])
                p2 = Polynomial([-6//7, 2//7, 3//7])
            end
            for i = 0:(ord-1)
                bsp = getbspline(ord, i)
                b1 = bsp*p1
                b2 = bsp*p2
                bspbis = b1 - b2
                v = bsplinerec(ord, i, x)
                @test  v == bsp(x)
                @test v == bspbis(x)
            end
        end
    end

end

@testset "test get_kl_ku" begin
    kl, ku = get_kl_ku(5)
    @test kl == 2 && ku == 2

    kl, ku = get_kl_ku(6)
    @test kl == 2 && ku == 3
    #    @test kl == 3 && ku == 2
end

test_bspline()


text:
export VerticalCAS_MDP

# Define MDP
mutable struct VerticalCAS_MDP <: MDP{stateType, actType}
    hs::Array{Float64,1}
    vowns::Array{Float64,1}
    vints::Array{Float64,1}
    pras::Array{Int64,1}
    resps::Array{Bool,1}
    taus::Array{Float64,1}
    discount_factor::Float64
    accels::Dict{Int,Tuple{Vector{Float64},Vector{Float64}}}
    velRanges::Dict{Int,Tuple{Float64,Float64}}
    allowedTrans::Dict{Int,Vector{Int}}
    currentTau::Float64
end
VerticalCAS_MDP() = VerticalCAS_MDP(hs,vowns,vints,acts,resps,taus,discount_f,accels,velRanges,allowedTrans,0.0)

# Define necessary functions for VerticalCAS MDP
POMDPs.actionindex(::VerticalCAS_MDP, a::actType) = a + 1
POMDPs.actions(mdp::VerticalCAS_MDP)  = acts
POMDPs.discount(mdp::VerticalCAS_MDP) = mdp.discount_factor
POMDPs.n_actions(::VerticalCAS_MDP)   = length(acts)

function POMDPs.convert_s(::Type{V} where V <: AbstractVector{Float64}, s::stateType, mdp::VerticalCAS_MDP)
    v = [s[1],s[2],s[3],convert(Float64,s[4]),convert(Float64,s[5])] #,s[6]]
    return v
end

function POMDPs.convert_s(::Type{stateType}, v::AbstractVector{Float64}, mdp::VerticalCAS_MDP)
    s = (v[1],v[2],v[3],convert(Int,v[4]), convert(Bool, v[5])) #,v[6])
    return s
end

max_stars_repo_name: cempolat/VerticalCAS

max_issues_repo_name: cempolat/VerticalCAS

max_forks_repo_name: cempolat/VerticalCAS


text:
# This file contains types and structures for technical DEA model
"""
    AbstractTechnicaDEAlModel
An abstract type representing a technical DEA model.
"""
abstract type AbstractTechnicalDEAModel end

"""
    efficiency(model::AbstractTechnicalDEAModel)
Return efficiency scores of a technical DEA model.
# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaio = dea(X, Y);

julia> efficiency(deaio)
11-element Array{Float64,1}:
 1.0
 0.6222896790980051
 0.8198562443845464
 1.0
 0.3103709311127934
 0.5555555555555556
 1.0
 0.7576690895651103
 0.8201058201058201
 0.49056603773584917
 1.0
```
"""
efficiency(model::AbstractTechnicalDEAModel) = model.eff

"""
    nbos(model::AbstractTechnicalDEAModel)
Return number of observations of a technical DEA model.
# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaio = dea(X, Y);

julia> nobs(deaio)
11
```
"""
nobs(model::AbstractTechnicalDEAModel) = model.n

"""
    ninputs(model::AbstractTechnicalDEAModel)
Return number of inputs of a technical DEA model.
## Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaio = dea(X, Y);

julia> ninputs(deaio)
2
```
"""
ninputs(model::AbstractTechnicalDEAModel) = model.m

"""
    noutputs(model::AbstractTechnicalDEAModel)
Return number of outputs of a technical DEA model.
# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaio = dea(X, Y);

julia> noutputs(deaio)
1
```
"""
noutputs(model::AbstractTechnicalDEAModel) = model.s

"""
    peers(model::AbstractTechnicalDEAModel)
Return peers of a technical DEA model.
# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaio = dea(X, Y);

julia> peers(deaio)
11×11 SparseArrays.SparseMatrixCSC{Float64,Int64} with 17 stored entries:
  [1 ,  1]  =  1.0
  [3 ,  1]  =  1.13432
  [2 ,  4]  =  0.424978
  [3 ,  4]  =  0.438005
  [4 ,  4]  =  1.0
  ⋮
  [7 ,  7]  =  1.0
  [8 ,  7]  =  0.114574
  [9 ,  7]  =  1.14815
  [10,  7]  =  0.490566
  [11, 11]  =  1.0
```
"""
peers(model::AbstractTechnicalDEAModel) = model.lambda

"""
    slacks(model::AbstractTechnicalDEAModel, slack::Symbol)
Return slacks a technical DEA model.
# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaio = dea(X, Y);

julia> slacks(deaio, :X)
11×2 Array{Float64,2}:
  0.0          0.0
 -4.41868e-15  0.0
  0.0          8.17926e-15
 -8.03397e-16  0.0
  1.80764e-15  0.0
  4.44444      0.0
  0.0          0.0
  1.60679e-15  0.0
  1.64021      0.0
  9.68683e-15  0.0
  0.0          4.0

julia> slacks(deaio, :Y)
11×1 Array{Float64,2}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
```
"""
function slacks(model::AbstractTechnicalDEAModel, slack::Symbol)::Matrix

    if slack == :X
        return model.slackX
    elseif slack == :Y
        return model.slackY
    end

end


text:
# Included when Blink.jl is loaded

using .Blink

# Create a Blink window
function viewblink(html::AbstractString,
                height::Integer,
                width::Integer,
                debug::Bool)
    w = Window()
    if debug
        opentools(w)
    end
    title(w, "Bio3DView")
    size(w, width, height)
    if Sys.iswindows()
        req_path = replace(path_jquery, "\\" => "\\\\")
    else
        req_path = path_jquery
    end
    loadhtml(w, html)
    return w
end


text:
function MakeRDict(
    R::Union{Array{<:Real,2},SparseArrays.SparseMatrixCSC{<:Real,Int64}},
    model::Model,
    mesh::Mesh,
    )

    ## Make a Dictionary so that the blocks of B are easy to access
    N₋ = sum(model.C.<=0)
    N₊ = sum(model.C.<=0)

    RDict = Dict{Tuple{String,Union{Int,Colon}},SparseArrays.SparseMatrixCSC{Float64,Int64}}()

    ppositions = cumsum(model.C .<= 0)
    qpositions = cumsum(model.C .>= 0)
    for ℓ in ["+", "-"]
        for i = 1:NPhases(model)
            FilBases = repeat(mesh.Fil[(ℓ,i)]', NBases(mesh), 1)[:]
            pitemp = falses(N₋)
            qitemp = falses(N₊)
            if model.C[i] <= 0
                pitemp[ppositions[i]] = mesh.Fil["p"*ℓ,i][1]
            end
            if model.C[i] >= 0
                qitemp[qpositions[i]] = mesh.Fil["q"*ℓ,i][1]
            end
            i_idx = [
                pitemp
                falses((i - 1) * TotalNBases(mesh))
                FilBases
                falses(NPhases(model) * TotalNBases(mesh) - i * TotalNBases(mesh))
                qitemp
            ]
            RDict[(ℓ,i)] = R[i_idx, i_idx]
        end
        FlBases =
            [mesh.Fil["p"*ℓ,:]; repeat(mesh.Fil[ℓ,:]', NBases(mesh), 1)[:]; mesh.Fil["q"*ℓ,:]]
        RDict[ℓ,:] = R[FlBases, FlBases]
    end
    return RDict
end

"""
# Construct the DG approximation to the operator `R`.

    MakeR(
        model::SFFM.Model,
        mesh::DGMesh;
        approxType::String = "projection",
    )

# Arguments
- `model`: a Model object
- `Mmesh`: a Mesh object
- `approxType::String`: (optional) either "interpolation" or
    "projection" (default).

# Output
- a tuple with keys
    - `:R::SparseArrays.SparseMatrixCSC{Float64,Int64}`: an approximation to R
        for the whole space. If ``rᵢ(x)=0`` on any cell, the corresponding
        elements of R are zero.
    - `:RDict::Dict{Tuple{String,Union{Int,Colon}},SparseArrays.SparseMatrixCSC{Float64,Int64}}`: a
        disctionary containing sub-blocks of R. Keys are of the form
        `"PhaseSign"` or just `"Sign"`. i.e. `"1-"` cells in ``Fᵢ⁻``, and
        `"-"` for cells in ``∪ᵢFᵢ⁻``.
"""
function MakeR(
    model::SFFM.Model,
    mesh::Mesh;
    approxType::String = "projection",
    v::Bool = false,
)
    V = SFFM.vandermonde(NBases(mesh))

    EvalR = 1.0 ./ model.r.a(CellNodes(mesh)[:])

    N₋ = sum(model.C .<= 0)
    N₊ = sum(model.C .>= 0)

    R = SparseArrays.spzeros(
        Float64,
        N₋ + N₊ + TotalNBases(mesh) * NPhases(model),
        N₋ + N₊ + TotalNBases(mesh) * NPhases(model),
    )
    # at the boundaries
    R[1:N₋, 1:N₋] = (1.0 ./ model.r.a(model.Bounds[1,1])[model.C .<= 0]).*LinearAlgebra.I(N₋)
    R[(end-N₊+1):end, (end-N₊+1):end] =  (1.0 ./ model.r.a(model.Bounds[1,end])[model.C .>= 0]).* LinearAlgebra.I(N₊)

    # on the interior
    for n = 1:(NIntervals(mesh)*NPhases(model))
        if Basis(mesh) == "legendre"
            if approxType == "interpolation"
                leftM = V.V'
                rightM = V.inv'
            elseif approxType == "projection"
                leftM = V.V' * LinearAlgebra.diagm(V.w)
                rightM = V.V
            end
            temp = leftM*LinearAlgebra.diagm(EvalR[NBases(mesh)*(n-1).+(1:NBases(mesh))])*rightM
        elseif Basis(mesh) == "lagrange"
            if approxType == "interpolation"
                temp = LinearAlgebra.diagm(EvalR[NBases(mesh)*(n-1).+(1:NBases(mesh))])
            elseif approxType == "projection"
                # the first term, LinearAlgebra.diagm(EvalR[NBases(mesh)*(n-1).+(1:NBases(mesh))])
                # is the quadrature approximation of M^r. The quadrature weights to not
                # appear since they cancel when we transform to integral/probability
                # representation. The second term V.V*V.V' is Minv. The last term
                # LinearAlgebra.diagm(V.w)is a result of the conversion to probability
                # / integral representation.
                temp = LinearAlgebra.diagm(EvalR[NBases(mesh)*(n-1).+(1:NBases(mesh))])*V.V*V.V'*LinearAlgebra.diagm(V.w)
            end
        elseif Basis(mesh) == ""
            temp = LinearAlgebra.diagm(EvalR[NBases(mesh)*(n-1).+(1:NBases(mesh))])
        end
        R[NBases(mesh)*(n-1).+(1:NBases(mesh)).+N₋, NBases(mesh)*(n-1).+(1:NBases(mesh)).+N₋] = temp
    end

    # construc the dictionary
    RDict = MakeRDict(R,model,mesh)

    out = (R=R, RDict=RDict)
    v && println("UPDATE: R object created with keys ", keys(out))
    return out
end

"""
Construct the operator `D(s)` from `B, R`.

    MakeD(
        mesh::SFFM.Mesh,
        B::Generator,
        R::NamedTuple{(:R, :RDict)},
    )

# Arguments
- `R`: a tuple as constructed by MakeR
- `B`: a tuple as constructed by MakeB
- `model`: a Model object
- `mesh`: a Mesh object

# Output
- `DDict::Dict{Tuple{String,String},Function(s::Real)}`: a dictionary of functions. Keys are
  of the for `"ℓm"` where `ℓ,m∈{+,-}`. Values are functions with one argument.
  Usage is along the lines of `D["+-"](s=1)`.
"""
function MakeD(
    # mesh::SFFM.Mesh,
    B::Generator,
    R::NamedTuple{(:R, :RDict)},
    v::Bool = false,
)
    DDict = Dict{Tuple{String,String},Any}()
    for ℓ in ["+", "-"], m in ["+", "-"]
        nℓ = sum(B.Fil["p"*ℓ,:]) + sum(B.Fil[ℓ,:]) * size(B.D,1) + sum(B.Fil["q"*ℓ,:])
        Idℓ = SparseArrays.sparse(LinearAlgebra.I,nℓ,nℓ)
        if any(B.Fil["p0",:]) || any(B.Fil["0",:]) || any(B.Fil["q0",:]) # in("0", model.Signs)
            n0 = sum(B.Fil["p0",:]) +
                sum(B.Fil["0",:]) * size(B.D,1) +
                sum(B.Fil["q0",:])
            Id0 = SparseArrays.sparse(LinearAlgebra.I,n0,n0)
            DDict[ℓ,m] = function (; s::Real = 0)
                return if (ℓ == m)
                    R.RDict[ℓ,:] * (
                        B[(ℓ,m),(:,:)] - s * Idℓ +
                        B[(ℓ,"0"),(:,:)] * inv(Matrix(s * Id0 - B[("0","0"),(:,:)])) * B[("0",m),(:,:)]
                    )
                else
                    R.RDict[ℓ,:] * (
                        B[(ℓ,m),(:,:)] +
                        B[(ℓ,"0"),(:,:)] * inv(Matrix(s * Id0 - B[("0","0"),(:,:)])) * B[("0",m),(:,:)]
                    )
                end
            end # end function
        else
            DDict[ℓ,m] = function (; s::Real = 0)
                return if (ℓ == m)
                    R.RDict[ℓ,:] * (B[(ℓ,m),(:,:)] - s * Idℓ)
                else
                    R.RDict[ℓ,:] * B[(ℓ,m),(:,:)]
                end
            end # end function
        end # end if ...
    end # end for ℓ ...
    v && println("UPDATE: D(s) operator created with keys ", keys(DDict))
    return (DDict = DDict)
end
# below to be deprecated
function MakeD(
    mesh::SFFM.Mesh,
    B::Generator,
    R::NamedTuple{(:R, :RDict)},
    v::Bool = false,
)
    DDict = Dict{Tuple{String,String},Any}()
    for ℓ in ["+", "-"], m in ["+", "-"]
        nℓ = sum(mesh.Fil["p"*ℓ,:]) + sum(mesh.Fil[ℓ,:]) * NBases(mesh) + sum(mesh.Fil["q"*ℓ,:])
        Idℓ = SparseArrays.sparse(LinearAlgebra.I,nℓ,nℓ)
        if any(mesh.Fil["p0",:]) || any(mesh.Fil["0",:]) || any(mesh.Fil["q0",:]) # in("0", model.Signs)
            n0 = sum(mesh.Fil["p0",:]) +
                sum(mesh.Fil["0",:]) * NBases(mesh) +
                sum(mesh.Fil["q0",:])
            Id0 = SparseArrays.sparse(LinearAlgebra.I,n0,n0)
            DDict[ℓ,m] = function (; s::Real = 0)
                return if (ℓ == m)
                    R.RDict[ℓ,:] * (
                        B[(ℓ,m),(:,:)] - s * Idℓ +
                        B[(ℓ,"0"),(:,:)] * inv(Matrix(s * Id0 - B[("0","0"),(:,:)])) * B[("0",m),(:,:)]
                    )
                else
                    R.RDict[ℓ,:] * (
                        B[(ℓ,m),(:,:)] +
                        B[(ℓ,"0"),(:,:)] * inv(Matrix(s * Id0 - B[("0","0"),(:,:)])) * B[("0",m),(:,:)]
                    )
                end
            end # end function
        else
            DDict[ℓ,m] = function (; s::Real = 0)
                return if (ℓ == m)
                    R.RDict[ℓ,:] * (B[(ℓ,m),(:,:)] - s * Idℓ)
                else
                    R.RDict[ℓ,:] * B[(ℓ,m),(:,:)]
                end
            end # end function
        end # end if ...
    end # end for ℓ ...
    v && println("UPDATE: D(s) operator created with keys ", keys(DDict))
    return (DDict = DDict)
end

"""
Construct and evaluate ``Ψ(s)``.

Uses newtons method to solve the Ricatti equation
``D⁺⁻(s) + Ψ(s)D⁻⁺(s)Ψ(s) + Ψ(s)D⁻⁻(s) + D⁺⁺(s)Ψ(s) = 0.``

    PsiFun(; s = 0, D, MaxIters = 1000, err = 1e-8)

# Arguments
- `s::Real`: a value to evaluate the LST at
- `D`: a `Dict{String,Function(s::Real)}` as output from MakeD
- `MaxIters::Int`: the maximum number of iterations of newtons method
- `err::Float64`: an error tolerance for terminating newtons method. Terminates
    when `max(Ψ_{n} - Ψ{n-1}) .< eps`.

# Output
- `Ψ(s)::Array{Float64,2}`: a matrix approxiamtion to ``Ψ(s)``.
"""
function PsiFun(
    D::Dict{Tuple{String,String},Any};
    s::Real = 0,
    MaxIters::Int = 1000,
    err::Float64 = 1e-8,
    v::Bool = false,
)

    exitflag = ""

    EvalD = Dict{Tuple{String,String},SparseArrays.SparseMatrixCSC{Float64,Int64}}(("+","-") => D["+","-"](s = s))
    Dimensions = size(EvalD["+","-"])
    for ℓ in [("+","+"), ("-","-"), ("-","+")]
        EvalD[ℓ] = D[ℓ](s = s)
    end
    Psi = zeros(Float64, Dimensions)
    A = EvalD["+","+"]
    B = EvalD["-","-"]
    C = EvalD["+","-"]
    # RA, QA = LinearAlgebra.schur(Matrix(A)) # uncomment for algorithm 1
    # RB, QB = LinearAlgebra.schur(Matrix(B)) # uncomment for algorithm 1
    OldPsi = Psi
    flag = 1
    for n = 1:MaxIters
        ## line below goes with algorithm 4, is quadratically convergent
        Psi = LinearAlgebra.sylvester(Matrix(A), Matrix(B), Matrix(C))
        ## line below goes with algorithm 1
        # Psi = mysyl(RA,QA,RB,QB,C)
        if maximum(abs.(OldPsi - Psi)) < err
            flag = 0
            exitflag = string(
                "Reached err tolerance in ",
                n,
                " iterations with error ",
                string(maximum(abs.(OldPsi - Psi))),
            )
            break
        elseif any(isnan.(Psi))
            flag = 0
            exitflag = string("Produced NaNs at iteration ", n)
            break
        end
        OldPsi = Psi
        ## Algorithm 4 of Bean, O'Reilly, Taylor, 2008,
        ## Algorithms for the Laplace–Stieltjes Transforms of First Return Times for Stochastic Fluid Flows,
        ## Methodol Comput Appl Probab (2008) 10:381–408, DOI 10.1007/s11009-008-9077-3
        A = EvalD["+","+"] + Psi * EvalD["-","+"]
        B = EvalD["-","-"] + EvalD["-","+"] * Psi
        C = EvalD["+","-"] - Psi * EvalD["-","+"] * Psi
        ## Algorithm 1 of the above citation
        # A, B dont change, need only comput schur decomp once
        # C = EvalD["+-"] + Psi * EvalD["-+"] * Psi
    end
    if flag == 1
        exitflag = string(
            "Reached Max Iters ",
            MaxIters,
            " with error ",
            string(maximum(abs.(OldPsi - Psi))),
        )
    end
    v && println("UPDATE: Iterations for Ψ(s=", s,") exited with flag: ", exitflag)
    return Psi
end

function mysyl(RA,QA,RB,QB,C)
    # sylvester but we give it the schur decomps
    # RA, QA = schur(A)
    # RB, QB = schur(B)

    D = -(adjoint(QA) * (C*QB))
    Y, scale = LinearAlgebra.LAPACK.trsyl!('N','N', RA, RB, D)
    LinearAlgebra.rmul!(QA*(Y * adjoint(QB)), inv(scale))
end

"""
Returns the DG approximation to the return probabilities ``ξ`` for the process
``Y(t)``.
NOTE: IMPLEMENTED FOR LAGRANGE BASIS ONLY

    MakeXi(
        B::Generator,
        Ψ::Array;
        mesh=1,
        model=1,
    )

# Arguments
- `B`: an object as returned by `MakeB`
- `Ψ::Array{Float64,2}`: an array as returned by `PsiFun`

# Output
- `ξ::Array{Float64,2}`: a row-vector of first return probabilities
"""
function MakeXi(
    B::Generator,
    Ψ::Array{Float64,2};
    mesh::Mesh = DGMesh(),
    model::Model = Model(),
)
    # BBullet = [B["--"] B["-0"]; B["0-"] B["00"]]
    # invB = inv(Matrix(Bbullet))

    # BBulletPlus = [B["-+"]; B["0+"]]

    # solve the linear system -[ξ 0] [Bmm Bm0; B₀₋ B00]^-1 [B₋₊; B₀₊]Ψ = ξ
    # writing out the system it turns out we only need the index -- and -0
    # blocks of the inverse. Wikipedia tells us that these are

    # the next two lines are a fancy way to do an inverse of a block matrix
    tempMat = inv(Matrix(B[("0","0"),(:,:)]))
    invBmm = inv(B[("-","-"),(:,:)] - B[("-","0"),(:,:)]*tempMat*B[("0","-"),(:,:)])
    invBm0 = -invBmm*B[("-","0"),(:,:)]*tempMat

    A = -(invBmm*B[("-","+"),(:,:)]*Ψ + invBm0*B[("0","+"),(:,:)]*Ψ + LinearAlgebra.I)
    b = zeros(1,size(B[("-","-"),(:,:)],1))

    A[:,1] .= 1.0 # normalisation conditions
    b[1] = 1.0 # normalisation conditions

    ξ = b/A

    return ξ
end

"""
Returns the DG approximation to some quantities regarding the limiting
distribution of a SFFM. See Ouput below
NOTE: IMPLEMENTED FOR LAGRANGE BASIS ONLY

    MakeLimitDistMatrices(
        B::Generator,
        D::Dict{Tuple{String,String},Any},
        R::Dict{Tuple{String,Union{Int,Colon}},SparseArrays.SparseMatrixCSC{Float64,Int64}},
        Ψ::Array{<:Real},
        ξ::Array{<:Real},
        mesh::Mesh,
        model::Model;
    )

# Arguments
- `B`: an object as returned by `MakeB`
- `D`: an object as returned by `MakeD`
- `R`: an object as returned by `MakeR`
- `Ψ::Array{Float64,2}`: an array as returned by `PsiFun`
- `ξ::Array{Float64,2}`: an row-vector as returned by `XiFun`

# Output
marginalX, p, K
- `marginalX::Array{Float64,2}`: a row-vector of the marginal limiting
    distribution of the first buffer ``X(t)``.
- `p::Array{Float64,2}`: a row-vector of the distribution of ``X(t)`` at the
    times of first return of ``Y(t)`` to 0.
- `K::Array{Float64,2}`: the array in the operator exponential of the stationary
    distribution

"""
function MakeLimitDistMatrices(
    B::Generator,
    D::Dict{Tuple{String,String},Any},
    R::Dict{Tuple{String,Union{Int,Colon}},SparseArrays.SparseMatrixCSC{Float64,Int64}},
    Ψ::Array{<:Real},
    ξ::Array{<:Real},
    mesh::SFFM.Mesh;
)
    B00inv = inv(Matrix(B[("0","0"),(:,:)]))
    invBmm = inv(B[("-","-"),(:,:)] - B[("-","0"),(:,:)]*B00inv*B[("0","-"),(:,:)])
    invBm0 = -invBmm*B[("-","0"),(:,:)]*B00inv

    # B00inv = inv(Matrix(B["00"]))
    # invBmm = inv(B["--"] - B["-0"]*B00inv*B["0-"])
    # invBm0 = -invBmm*B["-0"]*B00inv

    αp = ξ * -[invBmm invBm0]

    K = D["+","+"]() + Ψ * D["-","+"]()

    n₊, n₀ = size(B[("+","0"),(:,:)])
    n₋ = size(B[("-","+"),(:,:)],1)

    BBulletPlus = [B[("-","+"),(:,:)]; B[("0","+"),(:,:)]]

    αintegralPibullet = ((αp * BBulletPlus) / -K) * [R["+",:] Ψ*R["-",:]]
    αintegralPi0 = -αintegralPibullet * [B[("+","0"),(:,:)]; B[("-","0"),(:,:)]] * B00inv

    α = sum(αintegralPibullet) + sum(αintegralPi0) + sum(αp)

    p = αp ./ α
    integralPibullet = αintegralPibullet ./ α
    integralPi0 = αintegralPi0 ./ α

    marginalX = zeros(Float64, n₊ + n₋ + n₀)
    idx₊ = [mesh.Fil["p+",:]; repeat(mesh.Fil["+",:]', NBases(mesh), 1)[:]; mesh.Fil["q+",:]]
    marginalX[idx₊] = integralPibullet[1:n₊]
    idx₋ = [mesh.Fil["p-",:]; repeat(mesh.Fil["-",:]', NBases(mesh), 1)[:]; mesh.Fil["q-",:]]
    marginalX[idx₋] = integralPibullet[(n₊+1):end] + p[1:n₋]
    idx₀ = [mesh.Fil["p0",:]; repeat(mesh.Fil["0",:]', NBases(mesh), 1)[:]; mesh.Fil["q0",:]]
    marginalX[idx₀] = integralPi0[:] + p[(n₋+1):end]

    return marginalX, p, K
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
export CuTensor

mutable struct CuTensor{T, N}
    data::CuArray{T, N}
    inds::Vector{Cwchar_t}
    function CuTensor{T, N}(data::CuArray{T, N}, inds::Vector{Cwchar_t}) where {T<:Number, N}
        new(data, inds)
    end
    function CuTensor{T, N}(data::CuArray{N, T}, inds::Vector{<:AbstractChar}) where {T<:Number, N}
        new(data, Cwchar_t.(inds))
    end
end

CuTensor(data::CuArray{T, N}, inds::Vector{<:AbstractChar}) where {T<:Number, N} =
    CuTensor{T, N}(data, convert(Vector{Cwchar_t}, inds))

CuTensor(data::CuArray{T, N}, inds::Vector{Cwchar_t}) where {T<:Number, N} =
    CuTensor{T, N}(data, inds)

Base.size(T::CuTensor) = size(T.data)
Base.size(T::CuTensor, i) = size(T.data, i)
Base.length(T::CuTensor) = length(T.data)
Base.ndims(T::CuTensor) = length(T.inds)
Base.strides(T::CuTensor) = strides(T.data)
Base.eltype(T::CuTensor) = eltype(T.data)
Base.similar(T::CuTensor{Tv, N}) where {Tv, N} = CuTensor{Tv, N}(similar(T.data), copy(T.inds))
Base.collect(T::CuTensor) = (collect(T.data), T.inds)


text:
value(p::ElementPenalty, θ::Number, s::Number)       = s * value(p, θ)
value(p::ElementPenalty, θ::AbstractArray{<:Number})            = sum(x -> value(p, x), θ)
value(p::ElementPenalty, θ::AbstractArray{<:Number}, s::Number) = sum(x -> value(p, x, s), θ)
function value(p::ElementPenalty, θ::AbstractArray{T}, s::AbstractArray{S}) where {T <: Number, S <: Number}
    size(θ) == size(s) || error("lengths of parameter/weights do not match")
    result = zero(value(p, first(θ), first(s)))
    @inbounds for i in eachindex(θ, s)
        result += value(p, θ[i], s[i])
    end
    result
end

prox!(p::ProxableElementPenalty, θ::AbstractArray{<:Number}, s::Number) = map!(θj -> prox(p, θj, s), θ, θ)
function prox!(p::ProxableElementPenalty, θ::AbstractArray{<:Number}, s::AbstractArray{<:Number})
    @assert size(θ) == size(s)
    @inbounds for i in eachindex(θ, s)
        θ[i] = prox(p, θ[i], s[i])
    end
    θ
end
prox(p::ProxableElementPenalty, θ::AbstractArray{<:Number}, s::Number)       = prox!(p, copy(θ), s)
prox(p::ProxableElementPenalty, θ::AbstractArray{<:Number}, s::AbstractArray{<:Number}) = prox!(p, copy(θ), s)

deriv(p::ElementPenalty, θ::Number, s::Number) = s * deriv(p, θ)
grad(p::ElementPenalty, θ::AbstractArray{<:Number}) = grad!(similar(θ), p, θ)
function grad(p::ElementPenalty, θ::AbstractArray{T}, s::S) where {T<:Number, S<:Number}
    grad!(similar(θ, float(promote_type(T, S))), p, θ, s)
end
function grad(p::ElementPenalty, θ::AbstractArray{T}, s::AbstractArray{S}) where {T<:Number, S<:Number}
    grad!(similar(θ, float(promote_type(T, S))), p, θ, s)
end
function grad!(storage::AbstractArray{<:Number}, p::ElementPenalty, θ::AbstractArray{<:Number})
    map!(x -> deriv(p, x), storage, θ)
end
function grad!(storage::AbstractArray{<:Number}, p::ElementPenalty, θ::AbstractArray{<:Number}, s::Number)
    map!(x -> deriv(p, x, s), storage, θ)
end
function grad!(storage::AbstractArray{<:Number}, p::ElementPenalty, θ::AbstractArray{<:Number}, s::AbstractArray{<:Number})
    @assert size(storage) == size(θ) == size(s)
    @inbounds for j in eachindex(θ, s)
        storage[j] = deriv(p, θ[j], s[j])
    end
    storage
end

addgrad(∇j::Number, p::ElementPenalty, θj::Number) = ∇j + deriv(p, θj)
addgrad(∇j::Number, p::ElementPenalty, θj::Number, s::Number) = ∇j + s * deriv(p, θj)
function addgrad!(∇::AbstractArray{<:Number}, p::ElementPenalty, θ::AbstractArray{<:Number})
    @assert size(∇) == size(θ)
    @inbounds for j in eachindex(∇, θ)
        ∇[j] = addgrad(∇[j], p, θ[j])
    end
    ∇
end
function addgrad!(∇::AbstractArray{<:Number}, p::ElementPenalty, θ::AbstractArray{<:Number}, s::Number)
    @assert size(∇) == size(θ)
    @inbounds for j in eachindex(∇, θ)
        ∇[j] = addgrad(∇[j], p, θ[j], s)
    end
    ∇
end
function addgrad!(∇::AbstractArray{<:Number}, p::ElementPenalty, θ::AbstractArray{<:Number}, s::AbstractArray{<:Number})
    @assert size(∇) == size(θ) == size(s)
    @inbounds for j in eachindex(∇, θ, s)
        ∇[j] = addgrad(∇[j], p, θ[j], s[j])
    end
    ∇
end

value(p::ArrayPenalty, A::AbstractArray{<:Number}, λ::Number) = λ * value(p, A)

# --------------------
# AVAILABLE PENALTIES
# --------------------
include("penalties/elementwise.jl")
include("penalties/arraywise.jl")

# common functions
soft_thresh(x::Number, λ::Number) = sign(x) * max(zero(x), abs(x) - λ)

function soft_thresh!(x::AbstractArray{<:Number}, λ::Number)
    for i in eachindex(x)
        @inbounds x[i] = soft_thresh(x[i], λ)
    end
    x
end

# make penalties callable
for T in filter(!isabstracttype, union(subtypes(ElementPenalty),
                                      subtypes(ProxableElementPenalty),
                                      subtypes(ArrayPenalty)))
    @eval (pen::$T)(args...) = value(pen, args...)
end


text:
# ===== VSB Verifications: Vortex Sheet Calculations =====

println("=== BEGINNING Verifications5.jl ===")

# === Imports ===
println("   IMPORTING PACKAGES...")
using VSB
using SimpleVPM
using LinearAlgebra
using Plots
pyplot()
using LaTeXStrings
#using Printf

cd("/Users/Damyn/Documents/BYU/FLOW Lab/VSB/verifications")

# === Boundary Geometry and Plotting Grids ===
println("   GENERATING GEOMETRY...")
NPTS = 50                              # Number of Points on Boundary
del_theta = (2*pi)/(NPTS - 1)           # Step in theta around body
R = 1.0                                 # Cylinder Radius
X_coor = zeros(NPTS)                    # X-coordinates
Y_coor = zeros(NPTS)                    # Y-coordinates
X_tHat = zeros(NPTS)                    # X Tangent Vector
Y_tHat = zeros(NPTS)                    # X Tangent Vector
X_nHat = zeros(NPTS)                    # Y Normal Vector
Y_nHat = zeros(NPTS)                    # Y Normal Vector
theta = zeros(NPTS)                     # Angle from positive X-axis
for i=1:NPTS
    X_coor[i] = R*cos((i-1)*del_theta)
    Y_coor[i] = R*sin((i-1)*del_theta)

    X_tHat[i] = sin((i-1)*del_theta)
    Y_tHat[i] = -cos((i-1)*del_theta)

    X_nHat[i] = cos((i-1)*del_theta)
    Y_nHat[i] = sin((i-1)*del_theta)

    theta[i] = (i-1)*del_theta
end
body_pts = [[X_coor[i], Y_coor[i], 0.0] for i in 1:NPTS]
t_hats = [[X_tHat[i], Y_tHat[i], 0.0] for i in 1:NPTS]
n_hats = [[X_nHat[i], Y_nHat[i], 0.0] for i in 1:NPTS]
body = VSB.Boundary(body_pts,t_hats,n_hats)

# === Problem parameters ===
println("   CALCULATING PROBLEM PARAMETERS...")
magU_inf = 1.0
U_inf(X) = magU_inf .* [1,0,0]
params = VSB.Parameters(body, U_inf, verbose = true)
alphas = params.alphas

# === Analytical solution for gamma ===
println("   CALCULATING GAMMA ANALYTICAL...")
function gamma_analytical_i(i)
    point = body.bodyPTS[i]
    t_hat = body.tHats[i]
    gamma = -dot(U_inf(point), t_hat)
    return gamma
end
gammaA = [gamma_analytical_i(i) for i in 1:NPTS]

# === Numerical solution for gamma ===
println("   CALCULATING GAMMA NUMERICAL...")
function gamma_numerical_i(i)
    point = body.bodyPTS[i]
    gamma = VSB.CalcVS(body, alphas, point)
    return gamma
end
gammaN = [gamma_numerical_i(i) for i in 1:NPTS]

# === Plotting ===
println("   PLOTTING...")
plt1 = plot(theta, gammaA, line = (:red), label = "Analytical")
plot!(theta, gammaN, line = (:blue, :dash), label = "Numerical")
plot!(xlims = (0.0, 2*pi),
      xticks = ([0.0,pi/2,pi,3*pi/2,2*pi],["0",L"\frac{\pi}{2}",L"\pi",L"\frac{3\pi}{2}",L"2\pi"]),
      xlabel = L"\theta",
      ylims = (-1.05, 1.05),
      ylabel = L"\gamma(\theta)",
      grid = false)
# annotate!( [ (1,0.1,text("Analytical",14,:red)) , (1,-0.1,text("Numerical",14,:blue))])
savefig(plt1, "VSB_Verification_VS2.pdf")

println("=== END OF Verifications5.jl ===")


text:
using Test, Suppressor
using Bridge, BridgeSDEInference, StaticArrays, Distributions
using Statistics, Random, LinearAlgebra

const BSI = BridgeSDEInference
using BridgeSDEInference: ℝ

include("test_ODE_solver_change_pt.jl")
include("test_blocking.jl")
include("test_measchange.jl")
include("test_random_walk.jl")
include("test_workspace.jl")
include("test_setup.jl")


text:
using PyPlot

function plot_on_circumference(angularposition, val; ylabel = "", title = "", ylims = [])
    fig, ax = PyPlot.subplots()
    ax.plot(angularposition, val)
    if length(ylims) > 0
        ax.set_ylim(ylims)
    end
    ax.set_ylabel(ylabel)
    ax.set_xlabel("Angular position (deg)")
    ax.grid()
    fig.tight_layout()
    ax.set_title(title)
    return fig
end

function plot_stress_strain(
    angularposition,
    stresscomponent,
    straincomponent;
    ylabels = ["", ""],
    title = "",
)
    fig, ax = PyPlot.subplots(2, 1, sharex = true)
    ax[1].plot(angularposition, stresscomponent)
    ax[1].set_ylabel(ylabels[1])
    ax[1].grid()
    ax[2].plot(angularposition, straincomponent)
    ax[2].set_ylabel(ylabels[2])
    ax[2].grid()
    ax[2].set_xlabel("Angular position (deg)")
    fig.tight_layout()
    fig.suptitle(title)
    return fig
end

function plot_all_stress_strain(angularposition, stress, symmdispgrad, filename)
    fig11 = plot_stress_strain(
        angularposition,
        stress[1, :],
        symmdispgrad[1, :],
        ylabels = [L"\sigma_{11}", L"\epsilon_{11}"],
    )
    fig22 = plot_stress_strain(
        angularposition,
        stress[2, :],
        symmdispgrad[2, :],
        ylabels = [L"\sigma_{22}", L"\epsilon_{22}"],
    )
    fig12 = plot_stress_strain(
        angularposition,
        stress[3, :],
        symmdispgrad[3, :],
        ylabels = [L"\sigma_{12}", L"2\epsilon_{12}"],
    )

    s33 = stress[4, :]
    lowlim = 1.1min(0.0, minimum(s33))
    uplim = 1.1max(0.0, maximum(s33))
    fig33 = plot_on_circumference(
        angularposition,
        s33,
        ylabel = L"\sigma_{33}",
        ylims = [lowlim, uplim],
    )

    fig33.savefig(filename*"33.png")
    fig11.savefig(filename * "11.png")
    fig22.savefig(filename * "22.png")
    fig12.savefig(filename * "12.png")
end

# using PyPlot
# fig,ax = PyPlot.subplots()
# ax.plot(angularposition,productpressure,label="product")
# ax.plot(angularposition,parentpressure,label="parent")
# ax.set_title("Pressure around circumference")
# ax.set_xlabel("Angular position (deg)")
# ax.set_ylabel("Pressure (Pa)")
# ax.grid()
# ax.legend()
# fig.savefig("Pressure-around-circumference.png")


text:
module NDSparseIO

using CSV
using PooledArrays
using NullableArrays
using NDSparseData
using WeakRefStrings
using DataFrames
using NamedTuples

export csv2ndsparse

import NDSparseData: cmpelts, copyelt!, NDSparse

const STR = WeakRefString{UInt8}

# NDSparse optimizations for PooledArrays
@inline cmpelts(a::PooledArray, i, j) = (x=cmp(a.refs[i],a.refs[j]); x)
@inline copyelt!(a::PooledArray, i, j) = (a.refs[i] = a.refs[j])

process_col(a::NullableArray{STR}) = PooledArray{String}(a.values)
process_col(a::NullableArray) = a.values
process_col(a) = a

function csv2ndsparse(f; named::Bool=false, idxcols=nothing, datacol=nothing, datacols=nothing, kw...)
    return NDSparse(CSV.read(f; kw...), named=named, idxcols=idxcols, datacol=datacol, datacols=datacols)
end

function NDSparse(df::DataFrame; named::Bool=false, idxcols=nothing, datacol=nothing, datacols=nothing, kw...)
    cols = map(process_col, df.columns)
    n = length(cols)
    if datacols === nothing
        if datacol === nothing
            datacols = [n]
        elseif isa(datacol, Symbol)
            datacols = [df.colindex.lookup[datacol]]
        else
            datacols = [datacol]
        end
    elseif any(x->isa(x,Symbol), datacols)
        datacols = map(x->df.colindex.lookup[x], datacols)
    end
    if idxcols === nothing
        idxcols = sort(setdiff([1:n;], datacols))
    elseif any(x->isa(x,Symbol), idxcols)
        idxcols = map(x->df.colindex.lookup[x], idxcols)
    end
    if length(datacols) == 1
        data = cols[datacols[1]]
    end
    if named
        idx = Columns(cols[idxcols]..., names=df.colindex.names[idxcols])
        if length(datacols) > 1
            data = Columns(cols[datacols]..., names=df.colindex.names[datacols])
        end
    else
        idx = Columns(cols[idxcols]...)
        if length(datacols) > 1
            data = Columns(cols[datacols]...)
        end
    end
    NDSparse(idx, data)
end

end

max_stars_repo_name: JuliaComputing/NDSparseIO

max_issues_repo_name: JuliaComputing/NDSparseIO


text:
npics     = 6000 # number of pics to use for traning
newdim    = 416 # new size for pic upscaling; if you modify this you need to modify the network
# Load MNIST data:
include(Knet.dir("data","mnist.jl"))
# load raw MNIST: 60_000 pics
_mnist_xtrn,_mnist_ytrn,_mnist_xtst,_mnist_ytst = mnist()
# Select only a few pics
_mnist_xtrn = _mnist_xtrn[:,:,:,1:npics]
_mnist_ytrn = _mnist_ytrn[1:npics]
_mnist_xtst = _mnist_xtst[:,:,:,1:npics]
_mnist_ytst = _mnist_ytst[1:npics]
# Hacky upscale using Images.jl
_mnist_xtrn = imresize(_mnist_xtrn[:,:,1,:],newdim,newdim)
_mnist_xtrn = reshape(_mnist_xtrn,newdim,newdim,1,npics)
_mnist_xtst = imresize(_mnist_xtst[:,:,1,:],newdim,newdim)
_mnist_xtst = reshape(_mnist_xtst,newdim,newdim,1,npics)