text:
export fista

mutable struct FISTA{matT,vecT} <: AbstractLinearSolver
  A::matT
  reg::Regularization
  x::vecT
  xᵒˡᵈ::vecT
  res::vecT
  res_norm::Float64
  res_norm_old::Float64
  ρ::Float64
  t::Float64
  tᵒˡᵈ::Float64
  iterations::Int64
  relTol::Float64
  normalizeReg::Bool
  regFac::Float64
end

"""
    FISTA(A::matT, x::vecT=zeros(eltype(A),size(A,2))
          ; reg=nothing, regName=["L1"], λ=[0.0], kargs...) where {matT,vecT}

creates a `FISTA` object for the system matrix `A`.

# Arguments
* `A`                       - system matrix
* `x::vecT`                 - Array with the same type and size as the solution
* (`reg=nothing`)           - Regularization object
* (`regName=["L1"]`)        - name of the Regularization to use (if reg==nothing)
* (`λ=[0.0]`)               - Regularization paramter
* (`ρ::Float64=1`)          - step size for gradient step
* (`t::Float64=1.0`)        - parameter for predictor-corrector step
* (`relTol::Float64=1.e-5`) - tolerance for stopping criterion
* (`iterations::Int64=50`)  - maximum number of iterations
"""
function FISTA(A::matT, x::vecT=zeros(eltype(A),size(A,2)); reg=nothing, regName=["L1"]
              , λ=[0.0]
              , ρ::Float64=1.0
              , t::Float64=1.0
              , relTol::Float64=eps()
              , iterations::Int64=50
              , normalizeReg::Bool=false
              , kargs...) where {matT,vecT}

  if reg == nothing
    reg = Regularization(regName, λ, kargs...)
  end

  xᵒˡᵈ = similar(x)
  res = similar(x,size(A,1))

  return FISTA(A,vec(reg)[1],x,xᵒˡᵈ,res,0.0,0.0,ρ,t,t,iterations,relTol,normalizeReg,1.0)
end

"""
    init!(it::FISTA{matT,T}, b::vecT
              ; A::matT=solver.A
              , x::vecT=similar(b,0)
              , t::Float64=1.0) where T

(re-) initializes the FISTA iterator
"""
function init!(solver::FISTA{matT,vecT}, b::vecT
              ; A::matT=solver.A
              , x::vecT=similar(b,0)
              , t::Float64=1.0
              ) where {matT,vecT}

  solver.A = A
  if isempty(x)
    if !isempty(b) #!iszero(b)
      solver.x[:] .= adjoint(A) * b
    else
      solver.x[:] .= zeros(T,size(A,2))
    end
  else
    solver.x[:] .= x
  end
  solver.xᵒˡᵈ[:] .= solver.x  # this could also be zero
  solver.res[:] .= A*solver.x-b
  solver.res_norm_old = 0.0
  solver.res_norm = norm(solver.res)
  solver.t = t
  solver.tᵒˡᵈ = t
  # normalization of regularization parameters
  if solver.normalizeReg
    solver.regFac = norm(b,1)/length(b)
  else
    solver.regFac = 1.0
  end
end

"""
    solve(solver::FISTA, b::Vector)

solves an inverse problem using FISTA.

# Arguments
* `solver::FISTA`                 - the solver containing both system matrix and regularizer
* `b::vecT`                     - data vector
* `A::matT=solver.A`            - operator for the data-term of the problem
* (`startVector::vecT=similar(b,0)`)  - initial guess for the solution
* (`solverInfo=nothing`)          - solverInfo object

when a `SolverInfo` objects is passed, the residuals are stored in `solverInfo.convMeas`.
"""
function solve(solver::FISTA, b::vecT; A::matT=solver.A, startVector::vecT=similar(b,0), solverInfo=nothing, kargs...) where {matT,vecT}
  # initialize solver parameters
  init!(solver, b; A=A, x=startVector)

  # log solver information
  solverInfo != nothing && storeInfo(solverInfo,solver.x,solver.res_norm)

  # perform FISTA iterations
  for (iteration, item) = enumerate(solver)
    solverInfo != nothing && storeInfo(solverInfo,solver.x,solver.res_norm)
  end

  return solver.x
end

"""
  iterate(it::FISTA{matT,vecT}, iteration::Int=0) where {matT,vecT}

performs one fista iteration.
"""
function iterate(solver::FISTA{matT,vecT}, iteration::Int=0) where {matT,vecT}
  if done(solver, iteration) return nothing end

  # gradient step
  solver.xᵒˡᵈ[:] .= solver.x
  solver.x[:] .= solver.x - solver.ρ* (solver.A' * solver.res)

  # proximal map
  solver.reg.prox!(solver.x, solver.regFac*solver.ρ*solver.reg.λ; solver.reg.params...)

  # predictor-corrector update
  solver.tᵒˡᵈ = solver.t
  solver.t = (1. + sqrt(1. + 4. * solver.tᵒˡᵈ^2)) / 2.
  solver.x[:] .= solver.x + (solver.tᵒˡᵈ-1)/solver.t*(solver.x-solver.xᵒˡᵈ)

  # update residual
  # solver.res .= solver.A*solver.x-solver.b
  solver.res .+= solver.A*(solver.x-solver.xᵒˡᵈ)  # this relies on a proper initialization in init()
  solver.res_norm_old = solver.res_norm
  solver.res_norm = norm(solver.res)

  # return the residual-norm as item and iteration number as state
  return solver.res_norm, iteration+1
end

@inline converged(solver::FISTA) = ( abs(solver.res_norm-solver.res_norm_old)/solver.res_norm_old < solver.relTol )

@inline done(solver::FISTA,iteration::Int) = converged(solver) || iteration>=solver.iterations


text:
# # STOPPING CRIITERIA, AS CONTROLS

# `StoppingCriterion`objects are defined in EarlyStopping.jl

# non-wrapping stopping criteria that are nonsensical to apply if
# `IterationControl.training_losses(model)` is not overloaded:
const ATOMIC_CRIITERIA_NEEDING_TRAINING_LOSSES = [:PQ, ]
const ATOMIC_CRIITERIA_NEEDING_LOSS = [:Threshold, :GL, :PQ, :Patience]

# stopping criterion that wrap a single stopping criterion (must have
# `:criterion` as a field):
const EARLY_STOPPING_WRAPPERS = [:Warmup, ]

for ex in ATOMIC_CRIITERIA_NEEDING_LOSS
    quote
        needs_loss(::$ex) = true
    end |> eval
end

for ex in ATOMIC_CRIITERIA_NEEDING_TRAINING_LOSSES
    quote
        needs_training_losses(::$ex) = true
    end |> eval
end

for ex in EARLY_STOPPING_WRAPPERS
    quote
        needs_loss(wrapper::$ex) =
            needs_loss(wrapper.criterion)
        needs_training_losses(wrapper::$ex) =
            needs_training(wrapper.criterion)
    end |> eval
end

function update!(c::StoppingCriterion,
                 model,
                 verbosity,
                 n, state=nothing)
    _loss = loss(model)
    _training_losses = training_losses(model)
    if _training_losses !== nothing && !isempty(_training_losses)
        for tloss in _training_losses
            state = ES.update_training(c, tloss, state)
        end
    end
    state = ES.update(c, _loss, state)
    return state
end

function takedown(c::StoppingCriterion, verbosity, state)
    if done(c, state)
        message =  ES.message(c, state)
        verbosity > 0 && @info message
        return (done = true, log = message)
    else
        return (done = false, log = "")
    end
end

ext: jl


text:
# Compute windowed versions of the Single Layer potential applied to the known density function of the solution of an incident plane wave on a Dirichlet sphere.
using Cubature;
using GSL;

function field(k,the,r)
	res = exp(1im*k*cos(the)*r); # Incident field
	prevTerm = 20*k;
	term = 0;
	for n =0:round(1.4*k+40)
		n = Int64(n);
		term = (2*n+1)*1im^n*sf_bessel_jl(n,k)*(sf_bessel_jl(n,k*r)+1im*sf_bessel_yl(n,k*r))/(sf_bessel_jl(n,k)+1im*sf_bessel_yl(n,k))*sf_legendre_Plm(n,0,cos(the))
		if((abs(term) > abs(prevTerm) ) & (n > 10))
			return res;
		end
		res = res - term;
		prevTerm = term;
	end
	return res
end

function qDens(k,the)
	res = 1im*k*cos(the)*exp(1im*k*cos(the)); # Normal derivative of the incident field
	prevTerm = 20*k;
	term = 0;
	for n =0:round(1.4*k+40)
		n = Int64(n);
		if(n==0)
			derivHank = exp(1im*k)/k*(1+1im/k)
		else
			derivHank = (sf_bessel_jl(n-1,k)+1im*sf_bessel_yl(n-1,k) -(n+1)/k*(sf_bessel_jl(n,k)+1im*sf_bessel_yl(n,k)))
		end
		term = (2*n+1)*1im^n*sf_bessel_jl(n,k)*k*derivHank/(sf_bessel_jl(n,k)+1im*sf_bessel_yl(n,k))*sf_legendre_Plm(n,0,cos(the))
		if((abs(term) > abs(prevTerm) ) & (n > 10))
			return res;
		end
		res = res - term;
		prevTerm = term;
	end
	return res
end

function dis(x,y)
	phis = x[1];
	thetas=x[2];
	rstar = x[3];

	phi = y[1];
	theta=y[2];
	sqrt((cos(theta)-rstar*cos(thetas))^2 + (sin(theta)*sin(phi)-rstar*sin(thetas)*sin(phis))^2 + (sin(theta)*cos(phi)-rstar*sin(thetas)*cos(phis))^2);
end

function inte(x,y,k)
	theta=y[2]; # Jacobian is sin(theta)
	exp(k*1im*dis(x,y) )/4/pi/dis(x,y)*qDens(k,theta)*sin(theta)
end

function windInte(x,y,k,a,b)
	phis = x[1];
	thetas=x[2];
	phi = y[1];
	theta=y[2];
	r = sqrt( (sin(theta)*sin(phi)-sin(thetas)*sin(phis))^2 + (sin(theta)*cos(phi)-sin(thetas)*cos(phis))^2)
	if ((r <= a) || (abs(r-1) < 1e-8))
		return	exp(k*1im*dis(x,y) )/4/pi/dis(x,y)*qDens(k,theta)*sin(theta);
	elseif r < b
		return 	exp(k*1im*dis(x,y) )/4/pi/dis(x,y)*qDens(k,theta)*sin(theta)*exp(2*exp((a-b)/(r-a))/( (r-a)/(b-a)-1));
	else
		return 0;
	end
end

function pCmplx(rea,ima)
	if ima > 0
		return string(rea, " + ", ima, "im ");
	else
		return string(rea, " ",   ima, "im ");
	end
end

function intePol(x,y,k)
	rPol = y[1]
	thetaPol = y[2];
	if((rPol > 2) || (rPol < 0) || (x[3]^2*sin(thetaPol)^2/x[1]^2 + cos(thetaPol)^2 + (sin(thetaPol)-x[2]*cos(thetaPol)/x[3])^2 < 0) )
		print("Error: r=", rPol, ", theta=", thetaPol, ", under sqrt= " , x[3]^2*sin(thetaPol)^2/x[1]^2 + cos(thetaPol)^2 + (sin(thetaPol)-x[2]*cos(thetaPol)/x[3])^2, "\n")
	end
	xCo = -(1/2)*rPol*sqrt(-rPol^2+4)*x[3]*sin(thetaPol)/(sqrt(x[3]^2*sin(thetaPol)^2/x[1]^2 + cos(thetaPol)^2 + (sin(thetaPol)-x[2]*cos(thetaPol)/x[3])^2)*x[1])+(1-(1/2)*rPol^2)*x[1];
	if(abs(xCo) > 1)
		print("Error: r=", rPol, ", theta=", thetaPol, ", xCo= " , xCo, ", x = ", x, "\n")
	end
	exp(k*1im*rPol)/4/pi*qDens(k,acos(xCo))
end

xt = [0.8, pi/4, 1.0];
xtCar = xt[3]*[cos(xt[2]), sin(xt[2])*sin(xt[1]), sin(xt[2])*cos(xt[1])];
xo = [0.8, pi/4, 1.4];
xi = [2*pi/3, -0.1, 0.5];
percDecay = 0.5;

ks = 2.^(3:5);
Ts = linspace(0.9,1.6,2);
tol = 1e-8;

print(field(ks[1],xo[2],xo[3] ), "= total field ext, BC outside = ", exp(1im*ks[1]*cos(xo[2])*xo[3]), "\n")

# Test to check whether the integral without and with window satisfies the boundary condition:
(valr,err) = hcubature(y -> begin real(inte(xt,y,kt)); end, [0,0],[2*pi,pi])
(vali,eri) = hcubature(y -> begin imag(inte(xt,y,kt)); end, [0,0],[2*pi,pi])
print(kt , " =k, integral at x=", xt, " is ", pCmplx(valr,vali), ", b(x) = ", exp(1im*kt*cos(xt[2])), "\n")

bt = 1.9;
at = 1.2;
(valr,err) = hcubature(y -> begin real(windInte(xt,y,kt,at,bt)); end, [0,0],[2*pi,pi])
(vali,eri) = hcubature(y -> begin imag(windInte(xt,y,kt,at,bt)); end, [0,0],[2*pi,pi])
print(kt , " =k, wind integral at x=", xt, " is ", pCmplx(valr,vali), ", b(x) = ", exp(1im*kt*cos(xt[2])), ", a=", at, ", b=", bt, "\n")
# Alternatives to hcubature: hcubature_v, hquadrature, hquadrature_v: the _v should execute in parallel

# Test the integral using polar coordinates
if false
	# This integral should be small as there is no critical point in the range
	(valr,err) = hcubature(y -> begin real(intePol(xtCar,y,ks[ki])); end, [0.67,1.77],[0.7,1.8])
	(valr,err) = hcubature(y -> begin imag(intePol(xtCar,y,ks[ki])); end, [0.67,1.77],[0.7,1.8])
	print(valr, " + (", vali, ") polar coords \n");
elseif false
	(valr,err) = hcubature(y -> begin real(intePol(xtCar,y,ks[ki])); end, [0,0],[2,2*pi])
	(vali,eri) = hcubature(y -> begin imag(intePol(xtCar,y,ks[ki])); end, [0,0],[2,2*pi])
	print(valr, " + (", vali, ") polar coords \n");
end

if false
	# Integrate along some path of steepest descent
	function integr(p) exp(50*im*(-0.72764585274627917e-1*p^2+(.59216178402453730*1im)*p+(.40824829046386303*(-0.72764585274627917e-1*p^2+(.59216178402453730*1im)*p))*sqrt(-(-0.72764585274627917e-1*p^2+(.59216178402453730*1im)*p)^2+4)+(1/3*(1-(1/2)*(-0.72764585274627917e-1*p^2+(.59216178402453730*1im)*p)^2))*sqrt(3)))*(0.45572809000084120e-1-0.72764585274627917e-1*p^2+(.59216178402453730*1im)*p+(-0.72764585274627917e-1*p^2+(.59216178402453730*1im)*p)^2)*(-.14552917054925583*p+.59216178402453730*1im) end # thi=1,sdi=1
	maxP = (17/50)*log(10);
	(valr,err) = hcubature(x -> begin real(integr(x[1])); end, 0, maxP)
	(vali,eri) = hcubature(x -> begin imag(integr(x[1])); end, 0, maxP)
	print(valr+vali*1im, " =val, err = ", err+eri*1im, " , maxP =",maxP ," , integr at maxP = ", integr(maxP), '\n')
end

intVals = (0+0im)*zeros(length(ks),2)
intWindVals = (0+0im)*zeros(length(ks),length(Ts));
extVals = (0+0im)*zeros(length(ks),2)
extWindVals = (0+0im)*zeros(length(ks),length(Ts));
realVals = (0+0im)*zeros(length(ks),2)
windVals = (0+0im)*zeros(length(ks),length(Ts));
realErrs = (0+0im)*zeros(length(ks),1)
windErrs = (0+0im)*zeros(length(ks),length(Ts));

for ki = 1:length(ks)
	print("\n k = ", ks[ki], " ", now());
	valr = 0.0; err = 0.0; vali = 0.0; eri = 0.0;
	(valr,err) = hcubature(y -> begin real(inte(xt,y,ks[ki])); end, [0,0],[2*pi,pi], reltol = tol)
	(vali,eri) = hcubature(y -> begin imag(inte(xt,y,ks[ki])); end, [0,0],[2*pi,pi], reltol = tol)
	realVals[ki,2] = exp(1im*ks[ki]*cos(xt[2]));
	realVals[ki,1] = valr+1im*vali - realVals[ki,2];
	realErrs[ki,1] = err+1im*eri;

	(valr,err) = hcubature(y -> begin real(inte(xo,y,ks[ki])); end, [0,0],[2*pi,pi], reltol = tol)
	(vali,eri) = hcubature(y -> begin imag(inte(xo,y,ks[ki])); end, [0,0],[2*pi,pi], reltol = tol)
	extVals[ki,2] = field(ks[ki],xo[2], xo[3]);
	extVals[ki,1] = valr+1im*vali + extVals[ki,2]; # Maybe also try valr+1im*vali + exp(1im*ks[ki]*cos(xo[2])*xo[3]);
	print("Did BC and exterior, starting interior at ", now());

	(valr,err) = hcubature(y -> begin real(inte(xi,y,ks[ki])); end, [0,0],[2*pi,pi], reltol = tol)
	(vali,eri) = hcubature(y -> begin imag(inte(xi,y,ks[ki])); end, [0,0],[2*pi,pi], reltol = tol)
	intVals[ki,2] = exp(1im*ks[ki]*cos(xi[2])*xi[3]);
	intVals[ki,1] = valr+1im*vali - intVals[ki,2];

	for ti = 1:length(Ts)
		print("Starting windowed version with T = ", Ts[ti], " at ", now());
		valr = 0.0; err = 0.0; vali = 0.0; eri = 0.0;
		(valr,err) = hcubature(y -> begin real(windInte(xt,y,ks[ki],(1-percDecay)*Ts[ti],Ts[ti])); end, [0,0],[2*pi,pi], reltol = tol)
		(vali,eri) = hcubature(y -> begin imag(windInte(xt,y,ks[ki],(1-percDecay)*Ts[ti],Ts[ti])); end, [0,0],[2*pi,pi], reltol = tol)
		windVals[ki,ti] = valr+1im*vali - realVals[ki,2];
		windErrs[ki,ti] = err+1im*eri;

		(valr,err) = hcubature(y -> begin real(windInte(xo,y,ks[ki],(1-percDecay)*Ts[ti],Ts[ti])); end, [0,0],[2*pi,pi], reltol = tol)
		(vali,eri) = hcubature(y -> begin imag(windInte(xo,y,ks[ki],(1-percDecay)*Ts[ti],Ts[ti])); end, [0,0],[2*pi,pi], reltol = tol)
		extWindVals[ki,ti] = valr+1im*vali + extVals[ki,2];
		print(", T ext/int, ", now());

		(valr,err) = hcubature(y -> begin real(windInte(xi,y,ks[ki],(1-percDecay)*Ts[ti],Ts[ti])); end, [0,0],[2*pi,pi], reltol = tol)
		(vali,eri) = hcubature(y -> begin imag(windInte(xi,y,ks[ki],(1-percDecay)*Ts[ti],Ts[ti])); end, [0,0],[2*pi,pi], reltol = tol)
		intWindVals[ki,ti] = valr+1im*vali - intVals[ki,2];
	end
	f = open("simRes.txt","w");
	write(f,"xboundary=", string(xt), "xouter=", string(xo), "xinterior=", string(xi), "\n");
	write(f,"percDecay=", string(percDecay), ", tol=", tol, "\n")
	write(f,"ks=", string(ks), "\n Ts=", string(Ts), "\n")
	write(f,"intVals [scat-inc, inc] ", string(intVals), "\n")
	write(f,"intWindVals [scat-inc] ", string(intWindVals), "\n")
	write(f,"extVals [scat+exact field, exact field] ", string(extVals), "\n")
	write(f,"extWindVals [scat+exact field] ", string(extWindVals), "\n")
	write(f,"realVals [scat-inc, inc] ", string(realVals), "\n")
	write(f,"windVals (only [scattered - incident]) ", string(windVals), "\n")
	write(f,"realErrs", string(realErrs), "\n")
	write(f,"windErrs", string(windErrs), "\n")
	close(f);

	writedlm("intVals", intVals)
	writedlm("intWindVals", intWindVals)
	writedlm("extVals", extVals)
	writedlm("extWindVals", extWindVals)
	writedlm("realVals", realVals)
	writedlm("realErrs", realErrs)
	writedlm("windVals", windVals)
	writedlm("windErrs", windErrs)
end
print("\n\n")
run(`cat simRes.txt`)
print("\n\n")

max_stars_repo_name: popsomer/asyBEM

max_issues_repo_name: popsomer/asyBEM

max_forks_repo_name: popsomer/asyBEM


text:
using LIBSVM

# function input D is a distance matrix, K is a Kernel matrix
# count number of recurrences per time point i
"""
    REC(D::AbstractArray, rec_threshold::Float64, temp_excl::Int = 0)

Count the number of observations (recurrences) which fall into a radius `rec_threshold` of a distance matrix `D`. Exclude steps which are closer than `temp_excl` to be count as recurrences (default: `temp_excl = 5`)

Marwan, N., Carmen Romano, M., Thiel, M., & Kurths, J. (2007). Recurrence plots for the analysis of complex systems. Physics Reports, 438(5-6), 237–329. http://doi.org/10.1016/j.physrep.2006.11.001
"""
function REC(D::AbstractArray, rec_threshold::Float64, temp_excl::Int = 0)
  rec_out = init_REC(D)
  REC!(rec_out, D, rec_threshold, temp_excl)
  return(rec_out)
end

"""
    init_REC(D::Array{Float64, 2})
    init_REC(T::Int)

get object for memory efficient `REC!()` versions. Input can be a distance matrix `D` or the number of timesteps (observations) `T`.
"""
function init_REC(T::Int)
  rec_out = zeros(Float64, T)
end

function init_REC(D::Array{Float64, 2})
  rec_out = zeros(Float64, size(D, 1))
end

"""
    REC!(rec_out::AbstractArray, D::AbstractArray, rec_threshold::Float64, temp_excl::Int = 0)

Memory efficient version of `REC()` for use within a loop. `rec_out` is preallocated output, should be initialised with `init_REC()`.
"""
function REC!(rec_out::AbstractArray, D::AbstractArray, rec_threshold::Float64, temp_excl::Int = 0)
    N = size(D, 1)
    @assert temp_excl < N - 1
    @inbounds for i = 1:N
    rec_out[i] = 0.0
    if(i-temp_excl-1 >= 1)
      for j = 1:(i-temp_excl-1)
        if D[i, j] < rec_threshold && i != j
          rec_out[i] = rec_out[i] + 1
        end
      end
    end
    if(i+temp_excl+1 <= N)
      for j = (i+temp_excl+1):N
        if D[i, j] < rec_threshold && i != j
          rec_out[i] = rec_out[i] + 1
        end
      end
    end
  end
  broadcast!(/, rec_out, rec_out, N)
  broadcast!(-, rec_out, 1.0, rec_out)
  return(rec_out)
end

"""
    KDE(K)

Compute a Kernel Density Estimation (the Parzen sum), given a Kernel matrix `K`.

Parzen, E. (1962). On Estimation of a Probability Density Function and Mode. The Annals of Mathematical Statistics, 33, 1–1065–1076.
"""
function KDE(K::AbstractArray)
    KDEout = init_KDE(size(K, 1))
    KDE!(KDEout, K)
    return(KDEout)
end

"""
    init_KDE(K::Array{Float64, 2})
    init_KDE(T::Int)

Returns `KDE_out` object for usage in `KDE!()`. Use either a Kernel matrix `K` or the number of time steps/observations `T` as argument.
"""
function init_KDE(T::Int)
  KDE_out = zeros(Float64, T)
end

function init_KDE(K::Array{Float64, 2})
  KDE_out = zeros(Float64, size(K, 1))
end

"""
    KDE!(KDE_out, K)

Memory efficient version of `KDE()`. Additionally uses preallocated `KDE_out` object for writing the results. Initialize `KDE_out` with `init_KDE()`.
"""
function KDE!(KDE_out, K::AbstractArray)
  @assert size(K, 1) == size(KDE_out, 1)
  mean!(KDE_out, K)
  broadcast!(-, KDE_out, 1.0, KDE_out)
  return(KDE_out)
end

"""
    init_T2(VAR::Int, T::Int)
    init_T2{tp}(data::AbstractArray{tp,2})

initialize `t2_out` object for `T2!` either with number of variables `VAR` and observations/time steps `T` or with a two dimensional `data` matrix (time * variables)
"""
function init_T2(VAR::Int, T::Int)
  diagS =  zeros(Float64, VAR, VAR);
  Qinv =  zeros(Float64, VAR, VAR);
  data_norm = zeros(Float64, T, VAR);
  cdata = zeros(Float64, T, VAR);
  maha = zeros(Float64, T);
  t2_out = (maha, diagS, cdata, Qinv, data_norm)
  return(t2_out)
end

function init_T2(data::AbstractArray{tp,N}) where {tp, N}
  VAR = size(data, N)
  T = size(data, 1)
  diagS =  zeros(tp, VAR, VAR);
  Qinv =  zeros(tp, VAR, VAR);
  data_norm = zeros(tp, T, VAR);
  cdata = zeros(tp, T, VAR);
  maha = zeros(tp, T);
  t2_out = (maha, diagS, cdata, Qinv, data_norm)
  return(t2_out)
end

"""
    T2!(t2_out, data, Q[, mv])

Memory efficient version of `T2()`, for usage within a loop etc. Initialize the `t2_out` object with `init_T2()`.
`t2_out[1]` contains the squred Mahalanobis distance after computation.
"""
function T2!(t2_out::Tuple{Array{tp,1},Array{tp,2},Array{tp,2},Array{tp,2},Array{tp,2}}
             , data::AbstractArray{tp,2}, Q::AbstractArray{tp,2}, mv = 0) where {tp}
  (maha, diagS, cdata, Qinv, data_norm) = t2_out
  if(mv == 0)
    copyto!(cdata, data)
  elseif(size(mv, 1) == 1)
    copyto!(cdata, data .- mv)
  elseif(size(mv, 1) != 1)
    copyto!(cdata, data .- mv')
  end
  USVt = svd(Q)
  copyto!(view(diagS, diagind(diagS)), (USVt.S .+ 1e-10)  .^ (-0.5))
  transpose!(Qinv, USVt.U * diagS * USVt.Vt)
  copyto!(data_norm, cdata * Qinv)
  broadcast!(*, data_norm, data_norm, data_norm)
  sum!(maha, data_norm)
  return(t2_out[1])
end

# Hotelling's T^2 (Mahalanobis distance to the data mean)
# input is time * var data matrix
"""
    T2{tp}(data::AbstractArray{tp,2}, Q::AbstractArray[, mv])

Compute Hotelling's T2 control chart (the squared Mahalanobis distance to the data's mean vector (`mv`), given the covariance matrix `Q`).
Input data is a two dimensional data matrix (observations * variables).

Lowry, C. A., & Woodall, W. H. (1992). A Multivariate Exponentially Weighted Moving Average Control Chart. Technometrics, 34, 46–53.
"""
function T2(data::AbstractArray{tp,2}, Q::AbstractArray{tp, 2}, mv = 0) where {tp}
  t2_out = init_T2(data)
  T2!(t2_out, data, Q, mv)
  return(t2_out[1])
end

# mean of k nearest neighbor distances
"""
    KNN_Gamma(knn_dists_out)

This function computes the mean distance of the K nearest neighbors given a `knn_dists_out` object from `knn_dists()` as input argument.

Harmeling, S., Dornhege, G., Tax, D., Meinecke, F., & Müller, K.-R. (2006). From outliers to prototypes: Ordering data. Neurocomputing, 69(13-15), 1608–1618. http://doi.org/10.1016/j.neucom.2005.05.015
"""
function KNN_Gamma(knn_dists_out::Tuple{Int64,Array{Int64,1},Array{Float64,1},Array{Int64,2},Array{Float64,2}})
  NNdists = knn_dists_out[5]
  N = size(NNdists,1)
  KNN_Gamma_out = zeros(Float64, N)
  mean!(KNN_Gamma_out, NNdists)
  return(KNN_Gamma_out)
end

#T: number of timesteps in the datacube
"""
    init_KNN_Gamma(T::Int)
    init_KNN_Gamma(knn_dists_out)

initialize a `KNN_Gamma_out` object for `KNN_Gamma!` either with `T`, the number of observations/time steps or with a `knn_dists_out` object.
"""
function init_KNN_Gamma(T::Int)
  KNN_Gamma_out = zeros(Float64, T)
end

function init_KNN_Gamma(knn_dists_out::Tuple{Int64,Array{Int64,1},Array{Float64,1},Array{Int64,2},Array{Float64,2}})
  KNN_Gamma_out = zeros(Float64, size(knn_dists_out[2], 1))
end

"""
    KNN_Gamma!(KNN_Gamma_out, knn_dists_out)

Memory efficient version of `KNN_Gamma`, to be used in a loop. Initialize `KNN_Gamma_out` with `init_KNN_Gamma()`.
"""
function KNN_Gamma!(KNN_Gamma_out, knn_dists_out::Tuple{Int64,Array{Int64,1},Array{Float64,1},Array{Int64,2},Array{Float64,2}})
  NNdists = knn_dists_out[5]
  @assert size(NNdists,1) ==  size(KNN_Gamma_out,1) || error("input size KNN_Gamma_out and NNdists not equal")
  N = size(NNdists,1)
  mean!(KNN_Gamma_out, NNdists)
  return(KNN_Gamma_out)
end

"""
    init_KNN_Delta(T, VAR, k)

return a `KNN_Delta_out` object to be used for `KNN_Delta!`. Input: time steps/observations `T`, variables `VAR`, number of K nearest neighbors `k`.
"""
function init_KNN_Delta(T::Int, VAR::Int, k::Int)
  r = Array{Float64}(undef, T)
  x_i = Array{Float64}(undef, 1, VAR)
  d_x = Array{Float64}(undef, k, VAR)
  KNN_Delta_out = (r, x_i, d_x)
  return(KNN_Delta_out)
end

function init_KNN_Delta(knn_dists_out::Tuple{Int64,Array{Int64,1},Array{Float64,1},Array{Int64,2},Array{Float64,2}}, VAR::Int)
  T = size(knn_dists_out[2], 1)
  K = knn_dists_out[1]
  r = Array{Float64}(undef, T)
  x_i = Array{Float64}(undef, 1, VAR)
  d_x = Array{Float64}(undef, K, VAR)
  KNN_Delta_out = (r, x_i, d_x)
  return(KNN_Delta_out)
end

"""
    KNN_Delta!(KNN_Delta_out, knn_dists_out, data)

Memory Efficient Version of `KNN_Delta()`. `KNN_Delta_out[1]` is the vector difference of the k-nearest neighbors.
"""
function KNN_Delta!(KNN_Delta_out::Tuple{Array{Float64,1},Array{Float64,2},Array{Float64,2}}
                    , knn_dists_out::Tuple{Int64,Array{Int64,1},Array{Float64,1},Array{Int64,2},Array{Float64,2}}
                    , data::AbstractArray{Float64,2})
    (r, x_i, d_x) = KNN_Delta_out
    indices = knn_dists_out[4]
    K = size(indices, 2)
    T = size(data,1) #dimensions
    VAR = size(data, 2)
    @assert size(x_i, 2) == size(d_x, 2) == VAR || error("size of d_X and x_i not size of size(datacube, 4)")
    @assert size(r, 1) == T || error("size of r not of size(datacube ,1)")
    dists = 0.0
    for i = 1:T
        dists = 0.0
        inds=indices[i,:]
        for k=1:length(inds), j=1:size(data,2)
          d_x[k,j] = data[i,j] - data[inds[k],j]
        end
        sum!(x_i, d_x)
        for j = 1:VAR
            dists += (x_i[1,j] / K)^2
        end
        r[i] = sqrt.(dists)
    end
    return(r)
end

"""
    KNN_Delta(knn_dists_out, data)

Compute Delta as vector difference of the k-nearest neighbors. Arguments are a `knn_dists()` object (`knn_dists_out`) and a `data` matrix (observations * variables)

Harmeling, S., Dornhege, G., Tax, D., Meinecke, F., & Müller, K.-R. (2006). From outliers to prototypes: Ordering data. Neurocomputing, 69(13-15), 1608–1618. http://doi.org/10.1016/j.neucom.2005.05.015
"""
function KNN_Delta(knn_dists_out::Tuple{Int64,Array{Int64,1},Array{Float64,1},Array{Int64,2},Array{Float64,2}}
                    , data::AbstractArray{Float64,2})
  KNN_Delta_out = init_KNN_Delta(knn_dists_out, size(data, 2))
  KNN_Delta!(KNN_Delta_out, knn_dists_out, data)
  return(KNN_Delta_out[1])
end

"""
    init_UNIV(T::Int, VAR::Int)
    init_UNIV{tp}(data::AbstractArray{tp, 2})

initialize a `univ_out` object to be used in `UNIV!()` either with number of time steps/observations `T` and variables `VAR` or with a `data` matrix observations * variables.
"""
function init_UNIV(T::Int, VAR::Int)
    var_dat = zeros(tp, T)
    varquants = zeros(Float64,T, VAR)
    univ_out = (var_dat, varquants)
    return univ_out
end

function init_UNIV(data::AbstractArray{tp, N}) where {tp,N}
     T = size(data, 1)
     VAR = size(data, N)
     var_dat = zeros(tp, T)
     varquants = zeros(Float64,T, VAR)
     univ_out = (var_dat, varquants)
     return univ_out
end

"""
    UNIV!(univ_out, data)

Memory efficient version of `UNIV()`, input an `univ_out` object from `init_UNIV()` and some `data` matrix observations * variables
"""
function UNIV!(univ_out::Tuple{Array{tp,1},Array{Float64,2}}, data::AbstractArray{tp, N}) where {tp, N}
  (var_dat, varquants) = univ_out
  @assert size(var_dat, 1) == size(varquants, 1) == size(data, 1)
  @assert size(varquants, 2) == size(data, 2)
  for variable = 1:size(data, 2)
    # copy with little allocation
    copyto!(var_dat, view(data, :, variable))
    get_quantile_scores!(view(varquants, :, variable), var_dat)
  end
  maximum!(var_dat, varquants)
  return var_dat
end

"""
    UNIV(data)

order the values in each varaible and return their maximum, i.e. any of the variables in `data` (observations * variables) is above a given quantile,
the highest quantile will be returned.
"""
function UNIV(data::AbstractArray{tp, N}) where {tp, N}
  univ_out = init_UNIV(data)
  return UNIV!(univ_out, data)
end

"""
    SVDD_train(K, nu)

train a one class support vecort machine model (i.e. support vector data description), given a kernel matrix K and and the highest possible percentage of outliers `nu`.
Returns the model object (`svdd_model`). Requires LIBSVM.

Tax, D. M. J., & Duin, R. P. W. (1999). Support vector domain description. Pattern Recognition Letters, 20, 1191–1199.
Schölkopf, B., Williamson, R. C., & Bartlett, P. L. (2000). New Support Vector Algorithms. Neural Computation, 12, 1207–1245.
"""
function SVDD_train(K::AbstractArray, nu::Float64)
# function in LIBSVM.jl for optional parameter settings.
  #svdd_model = svmtrain(fill(1, size(K, 1)), K
  #                  , kernel_type = Int32(4), svm_type  = Int32(2), nu = nu
  #                  , probability_estimates = false);
  svdd_model = svmtrain(K,  svmtype = OneClassSVM, nu = nu, kernel = Kernel.Precomputed);
  return(svdd_model)
end

"""
    SVDD_predict(svdd_model, K)

predict the outlierness of an object given the testing Kernel matrix `K` and the `svdd_model` from SVDD_train(). Requires LIBSVM.

Tax, D. M. J., & Duin, R. P. W. (1999). Support vector domain description. Pattern Recognition Letters, 20, 1191–1199.
Schölkopf, B., Williamson, R. C., & Bartlett, P. L. (2000). New Support Vector Algorithms. Neural Computation, 12, 1207–1245.
"""
function SVDD_predict(svdd_model::LIBSVM.SVM, K::AbstractArray)
    (predicted_labels, decision_values) = svmpredict(svdd_model, K)
end

"""
    KNFST_predict(model, K)

predict the outlierness of some data (represented by the kernel matrix `K`), given some KNFST `model` from `KNFST_train(K)`. Compute `K`with `kernel_matrix()`.

Paul Bodesheim and Alexander Freytag and Erik Rodner and Michael Kemmler and Joachim Denzler:
"Kernel Null Space Methods for Novelty Detection". Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2013.

"""
function KNFST_predict(KNFST_mod, K)
  knfst_out = init_KNFST(size(K,2),KNFST_mod)
  KNFST_predict!(knfst_out, KNFST_mod, K)
  return(knfst_out)
end

"""
    init_KNFST(T, KNFST_mod)

initialize a `KNFST_out`object for the use with `KNFST_predict!`, given `T`, the number of observations and the model output `KNFST_train(K)`.
"""
function init_KNFST(T::Int, KNFST_mod::Tuple{Array{Float64,2},Array{Float64,2}})
  diffs = zeros(Float64, T, size(KNFST_mod[2], 2));
  scores = zeros(Float64, T);
  Ktransposed = zeros(Float64, T, size(KNFST_mod[1],1));
  KNFST_out = (scores, diffs, Ktransposed)
  return(KNFST_out)
end

"""
    KNFST_predict!(KNFST_out, KNFST_mod, K)

predict the outlierness of some data (represented by the kernel matrix `K`), given a `KNFST_out` object (`init_KNFST()`), some KNFST model (`KNFST_mod = KNFST_train(K)`)
and the testing kernel matrix K.

Paul Bodesheim and Alexander Freytag and Erik Rodner and Michael Kemmler and Joachim Denzler:
"Kernel Null Space Methods for Novelty Detection". Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2013.
"""
function KNFST_predict!(scores, diffs, Ktransposed, proj, targetValue, K)
  @assert size(scores, 1) == size(K, 2) == size(diffs, 1)
  @assert size(Ktransposed, 1) == size(K, 2) && size(Ktransposed, 2) == size(K, 1)
  transpose!(Ktransposed, K)
  # projected test samples: Ktransposed * model["proj"]
  mul!(diffs, Ktransposed, proj)
  # differences to the target value:
  broadcast!(-,diffs,diffs, targetValue)
  broadcast!(*,diffs, diffs, diffs)
  sum!(scores, diffs)
  broadcast!(sqrt, scores, scores)
  return(scores)
end

function KNFST_predict!(KNFST_out::Tuple{Array{Float64,1},Array{Float64,2},Array{Float64,2}}
                                     , KNFST_mod::Tuple{Array{Float64,2},Array{Float64,2}}, K::Array{Float64,2})
  (scores, diffs, Ktransposed) = KNFST_out
  (proj, targetValue) = KNFST_mod
  KNFST_predict!(scores, diffs, Ktransposed, proj, targetValue, K)
  return(scores)
end

# Learning method for novelty detection with KNFST according to the work:
#
# Paul Bodesheim and Alexander Freytag and Erik Rodner and Michael Kemmler and Joachim Denzler:
# "Kernel Null Space Methods for Novelty Detection".
# Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2013.
#
# Please cite that paper if you are using this code!
#
# proj = calculateKNFST(K, labels)
#
# calculates projection matrix of KNFST
#
# INPUT:
#   K -- (n x n) kernel matrix containing similarities of n training samples
#   labels -- (n x 1) vector containing (multi-class) labels of the n training samples
#
# OUTPUT:
#   proj -- projection matrix for data points (project x via kx*proj,
#           where kx is row vector containing kernel values of x and
#           training data)
#
#
function calculateKNFST(K, labels)

    classes = unique(labels);

    ### check labels
    if length(unique(labels)) == 1
      error("calculateKNFST.jl: not able to calculate a nullspace from data of a single class using KNFST (input variable 'labels' only contains a single value)");
    end

    ### check kernel matrix
    (n,m) = size(K);
    if n != m
        error("calculateKNFST.jl: kernel matrix must be quadratic");
    end

    ### calculate weights of orthonormal basis in kernel space
    centeredK = copy(K); # because we need original K later on again
    centerKernelMatrix(centeredK);
    (basisvecsValues,basisvecs) = eigen(centeredK);
    basisvecs = basisvecs[:,basisvecsValues .> 1e-12];
    basisvecsValues = basisvecsValues[basisvecsValues .> 1e-12];
    basisvecsValues = Diagonal(1 ./ sqrt.(basisvecsValues));
    basisvecs = basisvecs*basisvecsValues;

    ### calculate transformation T of within class scatter Sw:
    ### T= B'*Sw*B = H*H'  and H = B'*K*(I-L) and L a block matrix
    L = zeros(n,n);
    for i=1:length(classes)

       L[labels .== classes[i], labels .== classes[i]] .= 1.0 / sum(labels .== classes[i]);

    end

    ### need Matrix M with all entries 1/m to modify basisvecs which allows usage of
    ### uncentered kernel values:  (eye(size(M))-M)*basisvecs
    M = ones(m,m)./m;

    ### compute helper matrix H
    H = ((I - M)*basisvecs)' * K * (I - L);

    ### T = H*H' = B'*Sw*B with B=basisvecs
    T = H*H';

    ### calculate weights for null space
    eigenvecs = nullspace(T);

    if size(eigenvecs,2) < 1

      (eigenvals,eigenvecs) = eigen(T);
      (min_val,min_ID) = findmin(eigenvals);
      eigenvecs = eigenvecs[:,min_ID];

    end

    ### calculate null space projection and return it
    proj = ((I - M) * basisvecs) * eigenvecs;

end

#############################################################################################################

function centerKernelMatrix(kernelMatrix)
# centering the data in the feature space only using the (uncentered) Kernel-Matrix
#
# INPUT:
#       kernelMatrix -- uncentered kernel matrix
# OUTPUT:
#       centeredKernelMatrix -- centered kernel matrix

  ### get size of kernelMatrix
  n = size(kernelMatrix, 1);

  ### get mean values of each row/column
  columnMeans = mean(kernelMatrix, dims = 2); ### NOTE: columnMeans = rowMeans because kernelMatrix is symmetric
  matrixMean = mean(columnMeans);

  centeredKernelMatrix = kernelMatrix;

  for k=1:n
    for j=1:n
      centeredKernelMatrix[k,j] -= columnMeans[j];
      centeredKernelMatrix[j,k] -= columnMeans[j];
    end
  end
  #This line will not have any effect
  centeredKernelMatrix .+= matrixMean;

end

"""
    KNFST_train(K)

train a one class novelty KNFST model on a Kernel matrix `K` according to
Paul Bodesheim and Alexander Freytag and Erik Rodner and Michael Kemmler and Joachim Denzler:
"Kernel Null Space Methods for Novelty Detection". Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2013.

# Output
`(proj, targetValue)`
`proj` 	-- projection vector for data points (project x via kx*proj, where kx is row vector containing kernel values of x and training data)
`targetValue` -- value of all training samples in the null space
"""
function KNFST_train(K)

    # get number of training samples
    n = size(K,1);

    # include dot products of training samples and the origin in feature space (these dot products are always zero!)
    K_ext = [K zeros(n,1); zeros(1,n) 0];

    # create one-class labels + a different label for the origin
    labels = push!(ones(n),0);

    # get model parameters
    proj = calculateKNFST(K_ext,labels);
    targetValue = mean(K_ext[labels.==1,:]*proj, dims = 1);
    proj = proj[1:n,:];

    # return both variables
    return proj, targetValue

end

"""
    Dist2Centers(centers::AbstractArray{tp, 2}) where {tp}

Compute the distance to the nearest centers of i.e. a K-means clustering output.
Large Distances to the nearest center are anomalies. data: Observations * Variables.
# Example
`(proj, targetValue)`

"""
function Dist2Centers(data::AbstractArray{tp,2}, centers::AbstractArray{tp,2}) where {tp}
  (minD, minDIdx) = findmin(pairwise(Euclidean(), data', centers'), 2)
  clustassgin = zeros(Int, size(data, 1))
  for i = 1:size(data, 1)
    (u, clustassgin[i]) = ind2sub((size(data, 1), size(centers, 1)), minDIdx[i])
  end
  return minD, clustassgin
end

###################################
# end


text:
using LinearAlgebra
function loss(w, s=4.0)
	return (w[1]*w[2] - s)^2/2
end
function grad_loss(w, s=4.0)
	term = w[1]*w[2] - s
	return term*[w[2], w[1]]
end
function hess_loss(w, s=4.0)
	term = w[1]*w[2] - s
	dterm = [w[2] w[1]]
	return term*[0 1;1 0] .+
	[w[2]*dterm; w[1]*dterm]
end
next(w, η, s=4.0) = w - η*grad_loss(w, s)
dnext(w, η, s=4.0) = I(2) - η*hess_loss(w, s)
function d2next(w, η, s=4.0)
    term = w[1]*w[2] - s
    dterm = [w[2] w[1]]
    d2term = [0 1;1 0]
    dfirst = [0 0;dterm; dterm; 0 0]
    dsecond = [0 2*w[2]; w[2] w[1];
		   w[2] w[1]; 2*w[1] 0]
    return -η*(dfirst .+ dsecond)
end

max_stars_repo_name: nishaChandramoorthy/neuralOMET

max_issues_repo_name: nishaChandramoorthy/neuralOMET

max_forks_repo_name: nishaChandramoorthy/neuralOMET


text:
using YAXArrays, Documenter

exampledir = joinpath(@__DIR__,"src","examples")
allex = map(readdir(exampledir)) do fname
  n = splitext(fname)[1]
  n => joinpath("examples",fname)
end

makedocs(
    modules = [YAXArrays],
    clean   = true,
    format   = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "YAXArrays.jl",
    authors = "Fabian Gans",
    pages    = [ # Compat: `Any` for 0.4 compat
        "Home" => "index.md",
        "Tutorial" => "man/tutorial.md",
        "Examples" => allex,
        "Manual" => Any[
            "man/datasets.md",
            "man/yaxarrays.md",
            "man/applying functions.md",
            "man/iterators.md",
        ],
        "Docstring Reference" => "api.md"
        ]
)

deploydocs(
    #deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/JuliaDataCubes/YAXArrays.jl.git",
    push_preview = true
)

max_stars_repo_name: meggart/YAXArrays

max_forks_repo_name: meggart/YAXArrays


text:
using Pkg

handle("ping") do
  "pong"
end

handle("args") do args
  if length(filter!(!isempty, args)) !== 0
    append!(Base.ARGS, args)
  end
end

handle("echo") do data
  data
end

handle("cd") do path
  hideprompt() do
    cd(path)
    @info "Changed working directory to `$path`"
  end
end

handle("activateProject") do dir
  hideprompt() do
    Pkg.activate(dir)
  end
  update_project()
end

handle("activateParentProject") do dir
  hideprompt() do
    if (path = search_up_file("Project.toml", dir)) === nothing
      @warn "No project file found for `$dir`"
      return
    end
    Pkg.activate(path)
  end
  update_project()
end

handle("activateDefaultProject") do
  hideprompt() do
    Pkg.activate()
  end
  update_project()
end

function evalsimple(code, mod = "Main")
  mod = mod == "Main" ? Main : Core.eval(Main, :(module $(gensym(mod)) end))
  return with_logger(JunoProgressLogger()) do
    include_string(mod, code)
  end
end
handle(evalsimple, "evalsimple")

handle("exit") do
  exit()
end

handle("packages") do
  finddevpackages()
end

const PlotPaneEnabled = Ref(true)

handle("enableplotpane") do enable
  if enable
    Media.setdisplay(Media.Graphical, PlotPane())
    PlotPaneEnabled[] = true
  else
    Media.unsetdisplay(Media.Graphical)
    PlotPaneEnabled[] = false
  end
  nothing
end

handle("cancelCallback") do args...
  # TODO: Do something sensible here.
  # Until then it's better to silently fail than spam the REPL with "unrecognised
  # message" warnings.
end

handle("reportinfo") do
  io = IOBuffer()
  versioninfo(io)
  println(io)
  println(io, project_status())
  String(take!(io))
end


text:
src = "../src"
if !in(src,LOAD_PATH)
    push!(LOAD_PATH, src)
end
using Tel
using Base.Test


text:
using PMDPs
using Test
using POMDPSimulators

using StaticArrays
using POMDPs
using MCTS, DiscreteValueIteration
using StatsBase
using Plots
using POMDPSimulators
using D3Trees
using POMDPPolicies
using POMDPLinter
using Random
using DataFrames
using POMDPSimulators

using DrWatson
using Distributions
using BeliefUpdaters

"""
n_edges = 2,
c_init = 2,
selling_horizon_end = [10,10],
demand = Float64[2,2],
user_budgets = BudgetPerUnit(Distributions.Uniform(5,30)),
actions =  Action[0,15,30,45,1000]
"""
# params_mcts = Dict(pairs((solver= DPWSolver, n_iterations=1, depth=100, exploration_constant=40.0, keep_tree=true, show_progress=false)))
# params_mdp = Dict(pairs( (n_edges = 10, c_init = 30,  demand = 5*Float64[1,1,1,1,1,1,1,1,1,1], selling_horizon_end = collect(155:5:200), actions= 15:2:300, objective=:revenue)))
# mdp_mc = PMDPs.create_PMDP(PMDPg; params_mdp...)
# planner = PMDPs.get_MCTS_planner(mdp_mc, params_mcts)
# h_mc = run_sim(mdp_mc, planner; max_steps = 201, rng_seed = 1)

for d in [7500,10000,12500]
    println("Processng $d")
    # mdp_params = Dict(pairs( (n_edges = 3, c_init = 2, demand = 5*Float64[1,1,1], selling_horizon_end = [40,45,50], actions = 15:5:90, objective=:utilization)))
    mdp_params = Dict(pairs( (n_edges = 10, c_init = 30,  demand = 5*Float64[1,1,1,1,1,1,1,1,1,1], selling_horizon_end = collect(155:5:200), actions= 15:2:300, objective=:revenue)))
    # mcts_params = Dict(pairs((solver= MCTSSolver, n_iterations=5000, depth=500, exploration_constant=40.0, reuse_tree=true)))
    exp_params = Dict(pairs((n_runs = 10, vi=false, save=:stats)))
    mcts_params = Dict(pairs((solver= DPWSolver, n_iterations=d, depth=100, exploration_constant=40.0, keep_tree=true, show_progress=false)))
    params = Dict(:mdp=>mdp_params, :mcts=>mcts_params, :exp=>exp_params)

    result, filepath =  makesim(params);

    r = result[:r]
    u = result[:u]
    result[:t]

    println()
    display("rewards: $r")
    display("utilization: $u")
end

# for d in 5:1:6
#     println("Processing $d")
#     mdp_params = Dict(pairs( (n_edges = 3, c_init = 2, demand = 5*Float64[1,1,1], selling_horizon_end = [40,45,50], actions = 15:5:90, objective=:utilization)))
#     #mdp_params = Dict(pairs( (n_edges = 2, c_init = 1, demand = Float64[1,1], selling_horizon_end = [20,25], actions = 15:5:30)))
#     # mcts_params = Dict(solver= MCTSSolver, n_iterations=1000, depth=30, exploration_constant=40.0, reuse_tree=true)
#     mcts_params = Dict(pairs( (solver= DPWSolver, n_iterations=d, depth=30, exploration_constant=40.0, enable_state_pw = true, keep_tree=true, show_progress=false)))
#     exp_params = Dict(pairs((n_runs = 10, vi=true, save=:stats)))
#     params = Dict(:mdp=>mdp_params, :mcts=>mcts_params, :exp=>exp_params)

#     result, filepath =  makesim(params);

#     r = result[:r]
#     u = result[:u]
#     result[:t]
#     println()
#     display("rewards: $r")
#     display("utilization: $u")
# end

# map(length, values(mdp_params))
# sum(mdp_params[:demand])

# # analyze results
mdp_dir = readdir(datadir("sims"))[2]
exp = readdir(datadir("sims", mdp_dir))[1]
# res = wload(datadir("sims", mdp_dir, exp))

df = collect_results(
    datadir("sims", mdp_dir);
    white_list = [:r, :u, :t, :params]
)
df


text:
include("sparse/abstractsparse.jl")

module SparseMatrix

importall Base
import Base.NonTupleType, Base.float, Base.Order, Base.Sort.Forward

export SparseMatrixCSC,
       blkdiag, dense, diag, diagm, droptol!, dropzeros!, etree, full,
       getindex, ishermitian, issparse, issym, istril, istriu, nnz,
       setindex!, sparse, sparsevec, spdiagm, speye, spones,
       sprand, sprandbool, sprandn, spzeros, symperm, trace, tril, tril!,
       triu, triu!, nonzeros

include("sparse/sparsematrix.jl")
include("sparse/csparse.jl")

end # module SparseMatrix

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
include("CSD.jl")

"""
Epochs of Channel Sample.
"""
function epochsample(x,fs,epochs,chs;fun=nothing)
    nepoch = size(epochs,1)
    epochis = round.(Int,epochs.*SecondPerUnit.*fs)
    minepochlength = minimum(diff(epochis,dims=2))
    ys=Array{Float64}(undef,length(chs),minepochlength,nepoch)
    for i in 1:nepoch
        y = x[chs,range(max(1,epochis[i,1]),length=minepochlength)]
        if !isnothing(fun)
            y=fun(y)
        end
        ys[:,:,i] = y
    end
    return nepoch==1 ? dropdims(ys,dims=3) : ys
end

"reshape data stream in the same shape of chmask, where masked channels are replaced with local average"
function reshape2mask(ys, chmask; replacemask = true)
    nrow, ncol = size(chmask)
    yss = Array{Float64}(undef, nrow, ncol, size(ys)[2:end]...)
    for c = 1:ncol
        yss[:, c, :, :] = ys[c:ncol:end, :, :] # channel counting in chmask is from cols(left -> right), then rows(up -> down)
    end
    if replacemask # Local unmasked channels averaging replacement for masked channels
        rchmask = copy(chmask)
        for (r, c) in Tuple.(findall(chmask))
            rd = r - 1;ru = r + 1
            while 1 <= rd && rchmask[rd, c]
                rd -= 1
            end
            while ru <= nrow && rchmask[ru, c]
                ru += 1
            end
            if rd < 1
                if ru > nrow
                    @error "No Unmasked Channels in Data to Replace Masked Channel[$r, $c]."
                else
                    yss[r, c, :, :] = yss[ru, c, :, :]
                end
            else
                if ru > nrow
                    yss[r, c, :, :] = yss[rd, c, :, :]
                else
                    yss[r, c, :, :] = (yss[rd, c, :, :] .+ yss[ru, c, :, :]) / 2
                end
            end
            rchmask[r, c] = false
        end
    end
    return yss
end

"Spatial and Temporal filter of data stream"
function stfilter(rm;spatialtype=:none,ir=1,or=2,temporaltype=:none,ti=1:size(rm,2),hbordervalue=nothing)
    nd = ndims(rm)
    if nd==3
        nr,nc,n = size(rm)
        y = Array{Float64}(undef,nr,nc,n)
        for i in 1:n
            y[:,:,i] = stfilter(rm[:,:,i],spatialtype=spatialtype,ir=ir,or=or,temporaltype=temporaltype,ti=ti,hbordervalue=hbordervalue)
        end
        return y
    else
        if spatialtype==:all
            rm .-= mean(rm,dims=1)
        elseif spatialtype==:annulus
            foreach(i->rm[i,:] .-= dropdims(mean(rm[filter(j-> 1<=j<=size(rm,1),union(i-or:i-ir,i+ir:i+or)),:],dims=1),dims=1),1:size(rm,1))
        end
        if temporaltype == :all
            rm .-= mean(rm,dims=2)
        elseif temporaltype == :sub
            rm .-= mean(rm[:,ti],dims=2)
        end
        if !isnothing(hbordervalue)
            rm[[1,end],:].=hbordervalue
        end
        return rm
    end
end

"Remove line noise and its harmonics by notch filter"
function rmline!(y,fs;freq=60,nh=3,bw=3)
    for i=1:nh
        f = iirnotch(freq*i,bw;fs=fs)
        for j=1:size(y,1)
            y[j,:]=filtfilt(f,y[j,:])
        end
    end
    return y
end

"High pass and low pass filtering"
function hlpass(y,fs;high=0,low=Inf)
    fy=copy(y)
    if high>0
        f = digitalfilter(Highpass(high;fs=fs), Butterworth(4))
        for j=1:size(fy,1)
            fy[j,:]=filtfilt(f,fy[j,:])
        end
    end
    if low<Inf
        f = digitalfilter(Lowpass(low;fs=fs), Butterworth(4))
        for j=1:size(fy,1)
            fy[j,:]=filtfilt(f,fy[j,:])
        end
    end
    return fy
end

"time to sample index"
time2sampleindex(t,fs;secondperunit=1,minsampleindex::Integer=0,maxsampleindex::Integer=typemax(Int64)) = clamp(round(Int,t*secondperunit*fs),minsampleindex,maxsampleindex)
"sample index to time"
sampleindex2time(s,fs;secondperunit=1,mintime::Real=0,maxtime::Real=typemax(Float64)) = clamp(s/fs/secondperunit,mintime,maxtime)

"Convert epochs in time to sample index"
function epoch2sampleindex(epochs,fs;minsampleindex::Integer=1,maxsampleindex::Integer=typemax(Int64))
    nepoch = size(epochs,1)
    epochis = time2sampleindex.(epochs,fs;secondperunit=SecondPerUnit,minsampleindex,maxsampleindex)
    minepochlength = minimum(diff(epochis,dims=2))
    si = [range(epochis[i,1],length=minepochlength) for i in 1:nepoch]
    return nepoch==1 ? si[1] : si
end

"""
Get digital rising and falling edges in a analog stream
return di: edges index
       dv: edges active/inactive state
"""
function parsedigitalinanalog(x,ht,lt=ht,n=length(x);activehigh=true)
    di=Int[];dv=Bool[]
    if activehigh
        cs=x[1]>ht
        for i in 2:n
            t = x[i]
            if t > ht && !cs
                cs=true;push!(di,i);push!(dv,cs)
            elseif t < lt && cs
                cs=false;push!(di,i);push!(dv,cs)
            end
        end
    else
        cs=x[1]<lt
        for i in 2:n
            t = x[i]
            if t < lt && !cs
                cs=true;push!(di,i);push!(dv,cs)
            elseif t > ht && cs
                cs=false;push!(di,i);push!(dv,cs)
            end
        end
    end
    return di,dv
end

"Multi-Taper power spectrum estimation"
function powerspectrum(x,fs;freqrange=[0,100],nw=4)
    nd=ndims(x)
    if nd==3
        n=size(x,1);ne=size(x,3)
        ps = mt_pgram(x[1,:,1],fs=fs,nw=nw)
        fi = map(f->freqrange[1]<=f<=freqrange[2],freq(ps))
        freqs = freq(ps)[fi]
        p = Array{Float64}(undef,n,length(freqs),ne)
        for i in 1:n,j in 1:ne
            p[i,:,j] = power(mt_pgram(x[i,:,j],fs=fs,nw=nw))[fi]
        end
    elseif nd == 2
        ps = mt_pgram(x[1,:],fs=fs,nw=nw)
        fi = map(f->freqrange[1]<=f<=freqrange[2],freq(ps))
        freqs = freq(ps)[fi];n = size(x,1)
        p = Array{Float64}(undef,n,length(freqs))
        p[1,:]=power(ps)[fi]
        for i in 2:n
            p[i,:]=power(mt_pgram(x[i,:],fs=fs,nw=nw))[fi]
        end
    else
        ps = mt_pgram(x,fs=fs,nw=nw)
        fi = map(f->freqrange[1]<=f<=freqrange[2],freq(ps))
        freqs = freq(ps)[fi];p = power(ps)[fi]
    end
    return p,freqs
end

@doc raw"""
Discrete Fourier Transform at frequencies

```math
DFT[k] = \sum_{n=0}^{N-1} x[n] e^{\frac{-2\Pi ink}{N}}, k = 0:N-1, N = length(x), fₛ sampling x, fₖ = fₛ/N sampling DTFT
```

1. signal
2. simpling frequency of signal
3. at which frequencies DFT are directly evaluated
"""
function dft(x,fs,f...)
    N = length(x)
    fₖ = fs/N
    ks = round.(Int,f./fₖ)
    Fs = [zero(ComplexF64) for _ in f]
    Ω = [exp(-im*2π*n/N) for n in 0:(N-1)]
    @inbounds for n in 0:(N-1)
        @inbounds for i in eachindex(f)
            Fs[i] += x[n+1] * Ω[((n*ks[i])%N)+1]
        end
    end
    return Fs
end


text:
using ElasticWave

nz = 173; nx = 267; # model size
ext= 20; # thickness of PML layers
iflag = 2; #when iflag=1, four side obsorbing boundary, iflag=2, free surface
dx = 10. ; dz = 10. ; dt = 1e-3; # grid size
tmax = 6; f0=10.; # recording time and dominant frequency of source wavelet
vp=3000.*ones(nz,nx);  vs=3000./sqrt(3)*ones(nz,nx); rho=2.5*ones(nz,nx);
# vp[150:end,:]=4000.;  vs[150:end,:]=4000./sqrt(3); #physical model
fidMtx = CreateFidMtx(nz, nx, ext, iflag, vp, vs, rho, dz, dx, dt, f0);

(lambda, mu) = vel2lm(vp, vs, rho);
spt1 = InitSnapShot(nz, nx, ext, iflag, dt, 1);
spt2 = InitSnapShot(nz, nx, ext, iflag, dt, 2);
wfd1 = InitWfd(nz, nx, ext, iflag, dt, 1);
wfd2 = InitWfd(nz, nx, ext, iflag, dt, 2);

# give random number to spt1
n = length(spt1.Vzz);
spt1.Vzz = randn(n);  spt1.Vzx = randn(n);  spt1.Vxz = randn(n);  spt1.Vxx = randn(n);
spt1.Tzzz = randn(n); spt1.Tzzx = randn(n); spt1.Txxz = randn(n); spt1.Txxx = randn(n);
spt1.Txzz = randn(n); spt1.Txzx = randn(n);
OneStepForward!(spt2, spt1, fidMtx);

#convert spt to wfd
if iflag == 1
   Nz = nz + 2*ext
   zu = ext
elseif iflag == 2
   Nz = nz +   ext
   zu = 0
end
Nx = nx + 2*ext
tmp = reshape(spt1.Vzz+spt1.Vzx, Nz, Nx);
tmp = tmp[zu+1:zu+nz, ext+1:ext+nx]; wfd1.Vz = vec(tmp);
tmp = reshape(spt1.Vxz+spt1.Vxx, Nz, Nx);
tmp = tmp[zu+1:zu+nz, ext+1:ext+nx]; wfd1.Vx = vec(tmp);
tmp = reshape(spt1.Tzzz+spt1.Tzzx, Nz, Nx);
tmp = tmp[zu+1:zu+nz, ext+1:ext+nx]; wfd1.Tzz = vec(tmp);
tmp = reshape(spt1.Txxz+spt1.Txxx, Nz, Nx);
tmp = tmp[zu+1:zu+nz, ext+1:ext+nx]; wfd1.Txx = vec(tmp);
tmp = reshape(spt1.Txzz+spt1.Txzx, Nz, Nx);
tmp = tmp[zu+1:zu+nz, ext+1:ext+nx]; wfd1.Txz = vec(tmp);

# save the boundary of spt1
lb = (4*nz+(nx-4)*4)*5;
bd1 = zeros(lb);
tmp=reshape(wfd1.Vz,nz,nx);
il=   1; iu=   2*nz    ; bd1[il:iu]=vec(tmp[:,1:2]);
il=iu+1; iu=iu+2*nz    ; bd1[il:iu]=vec(tmp[:,end-1:end]);
il=iu+1; iu=iu+2*(nx-4); bd1[il:iu]=vec((tmp[1:2,3:nx-2])');
il=iu+1; iu=iu+2*(nx-4); bd1[il:iu]=vec((tmp[end-1:end,3:nx-2])');

tmp=reshape(wfd1.Vx,nz,nx);
il=iu+1; iu=iu+2*nz    ; bd1[il:iu]=vec(tmp[:,1:2]);
il=iu+1; iu=iu+2*nz    ; bd1[il:iu]=vec(tmp[:,end-1:end]);
il=iu+1; iu=iu+2*(nx-4); bd1[il:iu]=vec((tmp[1:2,3:nx-2])');
il=iu+1; iu=iu+2*(nx-4); bd1[il:iu]=vec((tmp[end-1:end,3:nx-2])');

tmp=reshape(wfd1.Tzz,nz,nx);
il=iu+1; iu=iu+2*nz    ; bd1[il:iu]=vec(tmp[:,1:2]);
il=iu+1; iu=iu+2*nz    ; bd1[il:iu]=vec(tmp[:,end-1:end]);
il=iu+1; iu=iu+2*(nx-4); bd1[il:iu]=vec((tmp[1:2,3:nx-2])');
il=iu+1; iu=iu+2*(nx-4); bd1[il:iu]=vec((tmp[end-1:end,3:nx-2])');

tmp=reshape(wfd1.Txx,nz,nx);
il=iu+1; iu=iu+2*nz    ; bd1[il:iu]=vec(tmp[:,1:2]);
il=iu+1; iu=iu+2*nz    ; bd1[il:iu]=vec(tmp[:,end-1:end]);
il=iu+1; iu=iu+2*(nx-4); bd1[il:iu]=vec((tmp[1:2,3:nx-2])');
il=iu+1; iu=iu+2*(nx-4); bd1[il:iu]=vec((tmp[end-1:end,3:nx-2])');

tmp=reshape(wfd1.Txz,nz,nx);
il=iu+1; iu=iu+2*nz    ; bd1[il:iu]=vec(tmp[:,1:2]);
il=iu+1; iu=iu+2*nz    ; bd1[il:iu]=vec(tmp[:,end-1:end]);
il=iu+1; iu=iu+2*(nx-4); bd1[il:iu]=vec((tmp[1:2,3:nx-2])');
il=iu+1; iu=iu+2*(nx-4); bd1[il:iu]=vec((tmp[end-1:end,3:nx-2])');

bv = zeros(bd1);
extractBoundary!(bv, spt1);

tmp = reshape(spt2.Vzz+spt2.Vzx, Nz, Nx);
tmp = tmp[zu+1:zu+nz, ext+1:ext+nx]; wfd2.Vz = vec(tmp);
tmp = reshape(spt2.Vxz+spt2.Vxx, Nz, Nx);
tmp = tmp[zu+1:zu+nz, ext+1:ext+nx]; wfd2.Vx = vec(tmp);
tmp = reshape(spt2.Tzzz+spt2.Tzzx, Nz, Nx);
tmp = tmp[zu+1:zu+nz, ext+1:ext+nx]; wfd2.Tzz = vec(tmp);
tmp = reshape(spt2.Txxz+spt2.Txxx, Nz, Nx);
tmp = tmp[zu+1:zu+nz, ext+1:ext+nx]; wfd2.Txx = vec(tmp);
tmp = reshape(spt2.Txzz+spt2.Txzx, Nz, Nx);
tmp = tmp[zu+1:zu+nz, ext+1:ext+nx]; wfd2.Txz = vec(tmp);

vz2 = reshape(wfd2.Vz, nz, nx);   vz1 = zeros(vz2);
vx2 = reshape(wfd2.Vx, nz, nx);   vx1 = zeros(vx2);
Tzz2 = reshape(wfd2.Tzz, nz, nx); Tzz1 = zeros(Tzz2);
Txx2 = reshape(wfd2.Txx, nz, nx); Txx1 = zeros(Txx2);
Txz2 = reshape(wfd2.Txz, nz, nx); Txz1 = zeros(Txz2);

vz = reshape(wfd1.Vz, nz, nx);
vx = reshape(wfd1.Vx, nz, nx);
Tzz = reshape(wfd1.Tzz, nz, nx);
Txx = reshape(wfd1.Txx, nz, nx);
Txz = reshape(wfd1.Txz, nz, nx);

a1 =  9/8; a2 = -1/24;
for j = 3:nx-2
    for i = 3 : nz-2
        c  = dt*(mu[i,j]+mu[i+1,j]+mu[i,j+1]+mu[i+1,j+1])/4
        c1 = c*a1/dz; c2 = c*a2/dz;
        t1 = c1*(vx2[i+1,j]-vx2[i  ,j])
        t2 = c2*(vx2[i+2,j]-vx2[i-1,j])
        c1 = c*a1/dx; c2 = c*a2/dx;
        t3 = c1*(vz2[i,j+1]-vz2[i,j  ])
        t4 = c2*(vz2[i,j+2]-vz2[i,j-1])
        Txz1[i,j] = Txz2[i,j]-(t1+t2+t3+t4)
    end
end

il = (4*nz+(nx-4)*4)*4;
Txz1[:,1:2] = reshape(bd1[il+1:il+2*nz], nz, 2); il = il+2*nz;
Txz1[:,end-1:end] = reshape(bd1[il+1:il+2*nz], nz, 2); il = il+2*nz;
Txz1[1:2,3:end-2] = (reshape(bd1[il+1:il+2*(nx-4)],(nx-4),2))'; il = il+2*(nx-4);
Txz1[end-1:end,3:end-2] = (reshape(bd1[il+1:il+2*(nx-4)],(nx-4),2))';

for j = 3:nx-2
    for i = 3 : nz-2
        c  = dt*(lambda[i,j]+2*mu[i,j])
        c1 = c*a1/dx; c2 = c*a2/dx;
        t1 = c1*(vx2[i,j  ]-vx2[i,j-1])
        t2 = c2*(vx2[i,j+1]-vx2[i,j-2])
        c  = dt*lambda[i,j]
        c1 = c*a1/dz; c2 = c*a2/dz;
        t3 = c1*(vz2[i  ,j]-vz2[i-1,j])
        t4 = c2*(vz2[i+1,j]-vz2[i-2,j])
        Txx1[i,j] = Txx2[i,j]-(t1+t2+t3+t4)
    end
end
il = (4*nz+(nx-4)*4)*3;
Txx1[:,1:2] = reshape(bd1[il+1:il+2*nz], nz, 2); il = il+2*nz;
Txx1[:,end-1:end] = reshape(bd1[il+1:il+2*nz], nz, 2); il = il+2*nz;
Txx1[1:2,3:end-2] = (reshape(bd1[il+1:il+2*(nx-4)],(nx-4),2))'; il = il+2*(nx-4);
Txx1[end-1:end,3:end-2] = (reshape(bd1[il+1:il+2*(nx-4)],(nx-4),2))';

for j = 3:nx-2
    for i = 3 : nz-2
        c  = dt*(lambda[i,j]+2*mu[i,j])
        c1 = c*a1/dz; c2 = c*a2/dz;
        t1 = c1*(vz2[i  ,j]-vz2[i-1,j])
        t2 = c2*(vz2[i+1,j]-vz2[i-2,j])
        c  = dt*lambda[i,j]
        c1 = c*a1/dx; c2 = c*a2/dx;
        t3 = c1*(vx2[i,j  ]-vx2[i,j-1])
        t4 = c2*(vx2[i,j+1]-vx2[i,j-2])
        Tzz1[i,j] = Tzz2[i,j]-(t1+t2+t3+t4)
    end
end
il = (4*nz+(nx-4)*4)*2;
Tzz1[:,1:2] = reshape(bd1[il+1:il+2*nz], nz, 2); il = il+2*nz;
Tzz1[:,end-1:end] = reshape(bd1[il+1:il+2*nz], nz, 2); il = il+2*nz;
Tzz1[1:2,3:end-2] = (reshape(bd1[il+1:il+2*(nx-4)],(nx-4),2))'; il = il+2*(nx-4);
Tzz1[end-1:end,3:end-2] = (reshape(bd1[il+1:il+2*(nx-4)],(nx-4),2))';

# Vx
for j = 3:nx-2
    for i = 3 : nz-2
        c  = 2*dt/(rho[i,j]+rho[i,j+1])
        c1 = c*a1/dx; c2 = c*a2/dx;
        t1 = c1*(Txx1[i,j+1]-Txx1[i,j  ])
        t2 = c2*(Txx1[i,j+2]-Txx1[i,j-1])
        c1 = c*a1/dz; c2 = c*a2/dz;
        t3 = c1*(Txz1[i  ,j]-Txz1[i-1,j])
        t4 = c2*(Txz1[i+1,j]-Txz1[i-2,j])
        vx1[i,j] = vx2[i,j]-(t1+t2+t3+t4)
    end
end
il = (4*nz+(nx-4)*4);
vx1[:,1:2] = reshape(bd1[il+1:il+2*nz], nz, 2); il = il+2*nz;
vx1[:,end-1:end] = reshape(bd1[il+1:il+2*nz], nz, 2); il = il+2*nz;
vx1[1:2,3:end-2] = (reshape(bd1[il+1:il+2*(nx-4)],(nx-4),2))'; il = il+2*(nx-4);
vx1[end-1:end,3:end-2] = (reshape(bd1[il+1:il+2*(nx-4)],(nx-4),2))';

# vz
for j = 3:nx-2
    for i = 3 : nz-2
        c  = 2*dt/(rho[i+1,j]+rho[i,j])
        c1 = c*a1/dz; c2 = c*a2/dz;
        t1 = c1*(Tzz1[i+1,j]-Tzz1[i  ,j])
        t2 = c2*(Tzz1[i+2,j]-Tzz1[i-1,j])
        c1 = c*a1/dx; c2 = c*a2/dx;
        t3 = c1*(Txz1[i,j  ]-Txz1[i,j-1])
        t4 = c2*(Txz1[i,j+1]-Txz1[i,j-2])
        vz1[i,j] = vz2[i,j]-(t1+t2+t3+t4)
    end
end
il = 0;
vz1[:,1:2] = reshape(bd1[il+1:il+2*nz], nz, 2); il = il+2*nz;
vz1[:,end-1:end] = reshape(bd1[il+1:il+2*nz], nz, 2); il = il+2*nz;
vz1[1:2,3:end-2] = (reshape(bd1[il+1:il+2*(nx-4)],(nx-4),2))'; il = il+2*(nx-4);
vz1[end-1:end,3:end-2] = (reshape(bd1[il+1:il+2*(nx-4)],(nx-4),2))';

# =====test write boundary to a binary file =====
using ElasticWave
nz = 100; nx = 100; # model size
ext= 20; # thickness of PML layers
iflag = 2; #when iflag=1, four side obsorbing boundary, iflag=2, free surface
dx = 10. ; dz = 10. ; dt = 1e-3; # grid size
tmax = 3.; f0=10.; # recording time and dominant frequency of source wavelet
vp=3000.*ones(nz,nx);  vs=3000./sqrt(3)*ones(nz,nx); rho=2.5e-4*ones(nz,nx);
vp[150:end,:]=4000.;  vs[150:end,:]=4000./sqrt(3);
fidMtx = CreateFidMtx(nz, nx, ext, iflag, vp, vs, rho, dz, dx, dt, f0);

isx = 50; isz = 1; ot=0.
flags = vec([false false true true false]); # add source to normal stress, simulate explosive sources
src = InitSource(isz, isx, nz, nx, ext, iflag, f0, ot, dt, flags);

root = homedir();
path = joinpath(root, "Desktop/bv.bin")
MultiStepForward(path, src, fidMtx, tmax)
bv = ReadWfdBoundary(path)

max_stars_repo_name: WneleiGao/ElasticWave


text:
#!/usr/bin/julia

function compose(f::Function, g::Function)
  return x -> f(g(x))
end

value = 0.5
for pair in [(sin, asin), (cos, acos), (x -> x^3, x -> x^(1/3))]
  func, inverse = pair
  println(compose(func, inverse)(value))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using BinaryBuilder

name = "GLFW"
version = v"3.3.2"

# Collection of sources required to build glfw
sources = [
    ArchiveSource("https://github.com/glfw/glfw/releases/download/$(version)/glfw-$(version).zip",
    "08a33a512f29d7dbf78eab39bd7858576adcc95228c9efe8e4bc5f0f3261efc7")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/glfw-*/
mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
    -DBUILD_SHARED_LIBS=ON \
    -DCMAKE_BUILD_TYPE=Release \
    -DGLFW_BUILD_EXAMPLES=false \
    -DGLFW_BUILD_TESTS=false \
    -DGLFW_BUILD_DOCS=OFF
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["libglfw", "glfw3"], :libglfw)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"),
    Dependency("Libglvnd_jll"),
    Dependency("Xorg_libXcursor_jll"),
    Dependency("Xorg_libXi_jll"),
    Dependency("Xorg_libXinerama_jll"),
    Dependency("Xorg_libXrandr_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
const empty_nb = "{\n \"cells\": [\n  {\n   \"outputs\": [],\n   \"cell_type\": \"markdown\",\n   \"source\": [\n    \"---\\n\",\n    \"\\n\",\n    \"*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*\"\n   ],\n   \"metadata\": {}\n  }\n ],\n \"nbformat_minor\": 3,\n \"metadata\": {\n  \"language_info\": {\n   \"file_extension\": \".jl\",\n   \"mimetype\": \"application/julia\",\n   \"name\": \"julia\",\n   \"version\": \"1.3.1\"\n  },\n  \"kernelspec\": {\n   \"name\": \"julia-1.3\",\n   \"display_name\": \"Julia 1.3.1\",\n   \"language\": \"julia\"\n  }\n },\n \"nbformat\": 4\n}\n"
const empty_script = "# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl\n\n"

function rm_if_empty(filepath::String)
    s = open(filepath) do file
        read(file, String)
    end
    if s == empty_nb || s == empty_script
        @info "Removing empty file: $filepath"
        rm(filepath)
    end
end

function postprocess_notebook(nb)
    txt = read(nb, String)
    open(nb, "w") do f
        write(
            f,
            replace(
                txt,
                r"\"outputs\":\ \[\]\,\n\ \ \ \"cell_type\":\ \"markdown\"" =>
                    "\"cell_type\": \"markdown\"",
            ),
        )
    end
end

"""
`literate_file(folder::AbstractString, file::AbstractString)`

Checks if the file has been modified since the last Weave and updates the notebook and tests accordingly.
* `folder` = Name of the folder the tutorial is in
* `file` = Name of the tutorial file
* `force` = force weave irrespective of file changes
"""
function literate_file(folder, file; force = false, kwargs...)
    filename = splitext(file)[1]
    srcpath = joinpath(SCRIPT_DIR, folder, file)
    testpath = joinpath(TEST_DIR, folder)
    notebookpath = joinpath(NB_DIR, folder)
    notebookfilepath = joinpath(notebookpath, join([filename, ".ipynb"]))
    configpath = joinpath(SCRIPT_DIR, folder, filename * "_config.json")

    config = get(kwargs, :config, Dict())
    if isfile(configpath)
        @info "found config file for $filename"
        config = read_json(configpath)
    end

    literate = get(config, "literate", true)

    if literate
        make_test = get(config, "test", true)
        make_notebook = get(config, "notebook", true)

        if make_test &&
           (mtime(srcpath) > mtime(testpath) || mtime(testpath) == 0.0 || force)
            @warn "Updating tests for $filename."
            fn = Literate.script(srcpath, testpath; config = config, kwargs...)
            rm_if_empty(fn)
        else
            @warn "Skipping tests for $filename."
        end
        if make_notebook && (
            mtime(srcpath) > mtime(notebookfilepath) ||
            mtime(notebookfilepath) == 0.0 ||
            force
        )
            @warn "Converting $filename to Jupyter Notebook."
            fn = Literate.notebook(srcpath, notebookpath; config = config, kwargs...)
            postprocess_notebook(fn)
            rm_if_empty(fn)
        else
            @warn "Skipping Jupyter Notebook for $filename."
        end
    else
        @warn "Skipping literate for $filename per config"
    end
end

"""
`literate_folder(folder::AbstractString)`

Checks the files present in the specified folder for modifications and updates the notebook and tests accordingly.

* `folder` = Name of the folder to check
"""
function literate_folder(folder; force = false, kwargs...)
    for file in readdir(joinpath(SCRIPT_DIR, folder))
        if splitext(file)[end] == ".jl"
            println("")
            println("Building $(joinpath(folder, file))")
            try
                literate_file(folder, file, force = force; kwargs...)
            catch
                @error "failed to build $(joinpath(folder, file))"
            end
            println("")
        end
    end
end

function literate_folder(example::Type{<:Examples}; kwargs...)
    folder = get_dir(example)
    literate_folder(folder; kwargs...)
end

"""
`literate_all()`

Checks every tutorial for modifications and updates the notebook accordingly.
"""
function literate_all(; force = false, kwargs...)
    for folder in readdir(SCRIPT_DIR)
        literate_folder(folder; force = force, kwargs...)
    end
end

function notebook(example::Type{<:Examples})
    literate_folder(example; execute = false)
    IJulia.notebook(dir = joinpath(NB_DIR, get_dir(example)))
end