text:
## evaluation metrics  for Gibbs sampler

#include("ntl_gibbs.jl")

function mean_arrival_time_Lp(T_inferred::Vector{Int},T_truth::Vector{Int},p::Real)
  """
  Returns mean (normalized by number of elements in `T_inferred`) L^`p` norm
    of difference between `T_inferred` and `T_truth`.
  """
  if size(T_inferred,1) != size(T_truth,1)
    error("T_inferred and T_truth must be the same length.")
  end
  K = size(T_inferred,1)
  if isfinite(p)
    d = (1/K)*sum( abs.((T_inferred .- T_truth).^p) ).^(1/p)
  else # L^∞-norm
    d = (1/K)*maximum( abs.(T_inferred .- T_truth) )
  end
  return d
end

function mean_arrival_time_Lp(T_inferred::Array{Int,2},T_truth::Vector{Int},p::Real)
  """
  Returns mean (normalized by number of columns of `T_inferred`) of
    mean (normalized by number of rows of `T_inferred`) L^`p` norm of
    differences between columns of `T_inferred` and `T_truth`.
  """
  s = zero(Float64)
  for j in 1:size(T_inferred,2)
    s += mean_arrival_time_Lp(T_inferred[:,j],T_truth,p)
  end
  return s./size(T_inferred,2)
end

function total_variation_distance(p::Vector{Float64},q::Vector{Float64})
  """
  Returns total variation distance between two discrete probabiliy distributions
    `p` and `q`. Elements of `p` and `q` are assumed to corresponding support points.
    If `size(p,1) != size(q,1)` then the shorter of the two is padded with zeros.
  """
  if abs(1.0 - sum(p)) > eps() || abs(1.0 - sum(q)) > eps()
    error("Elements of `p` or `q` do not sum to 1.")
  end
  np = size(p,1)
  nq = size(q,1)
  if np==nq
    d = 0.5*sum( abs.(p .- q) )
  elseif np > nq
    d = total_variation_distance(p,[q; zeros(Float64,np-nq)])
  else
    d = total_variation_distance([p; zeros(Float64,nq-np)],q)
  end
  return d
end

function total_variation_distance(p::Vector{Int64},q::Vector{Int64})
  """
  Returns total variation distance between two discrete probabiliy distributions
    `p` and `q`. Elements of `p` and `q` are assumed to corresponding support points.
    If `size(p,1) != size(q,1)` then the shorter of the two is padded with zeros.
  """
  # if abs(1.0 - sum(p)) > eps() || abs(1.0 - sum(q)) > eps()
  #   error("Elements of `p` or `q` do not sum to 1.")
  # end
  np = size(p,1)
  nq = size(q,1)
  if np==nq
    d = 0.5*sum( abs.(p .- q) )
  elseif np > nq
    d = total_variation_distance(p,[q; zeros(Int64,np-nq)])
  else
    d = total_variation_distance([p; zeros(Int64,nq-np)],q)
  end
  return d
end

function deviance_information_criterion(PP::Vector{Int},T_gibbs::Array{Float64,2},alpha_gibbs::Vector{Float64})
  """
  Returns DIC based on arrival times and α (marginalizes Ψ_j's)
  """
  T_mean = mean(T_gibbs,2)
  alpha_mean = mean(alpha_gibbs)
  D_param_bar = log_CPPF(PP,T_mean,alpha_mean)

  M = size(T_gibbs,2)
  D_bar = zero(Float64)
  for m in 1:M
    D_bar += log_CPPF(PP,T[:,m],alpha[m])
  end

  return 2*D_bar/M - D_param_bar

end

function predictive_logprob(PP_train::Vector{Int},T_train::Vector{Int},
  PP_train_test::Vector{Int},T_train_test::Vector{Int},
  predictive_ia_dist::DiscreteDistribution,alpha::Float64)

  f = (x,y) -> predictive_ia_dist
  predictive_logprob(PP_train,T_train,PP_train_test,T_train_test,f,alpha)
end

function predictive_logprob(PP_train::Vector{Int},T_train::Vector{Int},
  PP_train_test::Vector{Int},T_train_test::Vector{Int},
  predictive_ia_dist::Function,alpha::Float64)
  """
  Returns predictive log-probability of `PP_train_test` and `T_train_test` given `PP_train` and
    `predictive_ia_dist`, which is a function that evaluates interarrivals in `T_train_test`
    at a specific parameter setting (e.g., a sample from the posterior given the
    training data)

  `PP_train_test` and `T_train_test` should be the stochastic extension (growth) of `PP_train` and `T_train`,
    i.e., when `PP_train` has `K_train` elements, then the first `K_train` elements
    of `PP_train_test` correspond to those same elements and `PP_train_test[k] >= PP_train[k]` should
    be true for all `k <= K_train`. Similarly, `T_train_test` shoud be an extension of `T_train`.
  """
  logp_train = logp_partition(PP_train,T_train,alpha,predictive_ia_dist,false)
  logp_test_train = logp_partition(PP_train_test,T_train_test,alpha,predictive_ia_dist,false)
  return logp_test_train - logp_train
end

function predictive_logprob(logp_train::Float64,
  PP_train_test::Vector{Int},T_train_test::Vector{Int},
  predictive_ia_dist::DiscreteDistribution,alpha::Float64)
  """
  Returns predictive log-probability when log-probability of training data
    (`logp_train`) has already been calculated.
  """

  f = (x,y) -> predictive_ia_dist
  predictive_logprob(logp_train,PP_train_test,T_train_test,f,alpha)
end

function predictive_logprob(logp_train::Float64,
  PP_train_test::Vector{Int},T_train_test::Vector{Int},
  predictive_ia_dist::Function,alpha::Float64)
  """
  Returns predictive log-probability when log-probability of training data
    (`logp_train`) has already been calculated.
  """

  logp_test_train = logp_partition(PP_train_test,T_train_test,alpha,predictive_ia_dist,true)
  return logp_test_train - logp_train
end

function sample_predicted_arrival_times(ia_dist::Function,T_end::Int,
  K_end::Int,n_end::Int,n_preds::Int)
  """
  Samples predicted interarrival times from `ia_dist` and adds them to `T_end`,
    which corresponds to the `K_end`-th arrival time, until `n_end + n_preds` is exceeded.
    The first sampled arrival time is forced to be greater than `n_end`.
  """
  zero_shift = 1 - minimum(ia_dist(1,1))
  K_start = K_end
  # sample first arrival
  new_arr = T_end + rand(ia_dist(T_end,K_start)) + zero_shift
  while new_arr < n_end
    new_arr = T_end + rand(ia_dist(T_end,K_start)) + zero_shift
  end

  if new_arr > (n_end + n_preds)
    return []
  else
    T = [new_arr]
    j = 1
    K = K_start + j
    while T[end] < (n_end + n_preds)
      j += 1
      append!(T,rand(ia_dist(T[j-1],K_start+j-1)) + T[j-1] + zero_shift)
    end
    if T[end] > (n_end + n_preds)
      pop!(T)
    end
    return T
  end

end

function sample_predicted_sequence(PP_train::Vector{Int},T_end::Int,
  ia_dist::DiscreteDistribution,alpha::Float64,n_preds::Int)

  f = (x,y) -> ia_dist
  sample_predicted_sequence(PP_train,T_end,f,alpha,n_preds)

end

function sample_predicted_sequence(PP_train::Vector{Int},T_end::Int,
  ia_dist::Function,alpha::Float64,n_preds::Int)

  K_train = size(PP_train,1)
  n_train = sum(PP_train)
  T_pred = sample_predicted_arrival_times(ia_dist,T_end,K_train,n_train,n_preds)
  K_test = size(T_pred,1)
  K_end = K_train + K_test

  PP = vcat(PP_train,zeros(Int64,size(T_pred,1)))
  Z = zeros(Int64,n_preds)

  if K_test==0
    for n in 1:n_preds
      Z[n] = wsample(1:K_train, PP .- alpha)
      PP[Z[n]] += 1
    end
    return Z,PP,T_pred
  else

    k = K_train
    for n in 1:n_preds
      if n+n_train <= T_pred[end] && n+n_train == T_pred[k+1-K_train]
        k += 1
        PP[k] = 1
        Z[n] = k
        k > K_end ? k = K_end : nothing
      else
        Z[n] = wsample(1:k, PP[1:k] .- alpha) # discounted size-biased sample
        PP[Z[n]] += 1
      end
    end
    return Z,PP,T_pred
  end
end

function sample_predicted_partition(PP_train::Vector{Int},T_end::Int,
  ia_dist::DiscreteDistribution,alpha::Float64,n_preds::Int)

  f = (x,y) -> ia_dist
  return sample_predicted_partition(PP_train,T_end,f,alpha,n_preds)

end

function sample_predicted_partition(PP_train::Vector{Int},T_end::Int,
  ia_dist::Function,alpha::Float64,n_preds::Int)

  Z_pred,PP,T_pred = sample_predicted_sequence(PP_train,T_end,ia_dist,alpha,n_preds)
  return PP,T_pred

end


text:
# add Bukdu     # in pkg mode

using Bukdu

struct WelcomeController <: ApplicationController     # 1
  conn::Conn
end

struct RestController <: ApplicationController
  conn::Conn
end

function index(c::WelcomeController)                   # 2
  render(JSON, "Hello World from Bukdu!")
end

function init(c::RestController)
  render(JSON, (:init, c.params.region, c.params.site_id, c.params.channel_id))
end

routes() do
  get("/", WelcomeController, index)                        # 3A
  get("/init/region/:region/site/:site_id/channel/:channel_id/",
    RestController, init, :site_id=>Int, :channel_id=>Int)  # 3B
end

Bukdu.start(8080)                                       # 4
# Bukdu Listening on 127.0.0.1:8080
# Task (runnable) @0x00000000735447d0

# INFO: GET     WelcomeController   index           200 /

ext: jl


text:
using DelimitedFiles
using MPI

function forcecondition(sys::AbstractSystem,t::Float64)
  if sys.t >= t
    return true
  else
    return false
  end
end

function runtrajs!(ntraj::Int64,nsteps::Int64,dt::Float64,t0::Float64,
                   x0::Vector{Float64},system::AbstractThermostattedSystem,
                   integrator::StochasticEuler,mcvb::MCVBTCallback)

  # run trajectories for evaluating gradients
  for traj in 1:ntraj
    system.x = copy(x0)
    system.t = t0
    #integrator.model.potentials[end].condition(sys::AbstractSystem) = forcecondition(sys,dt*rand(0:nsteps))
    trand::Float64 = dt*rand(0:nsteps)
#    println("traj $traj ",trand)
    f(sys::AbstractSystem) = forcecondition(sys,trand)
    integrator.model.potentials[end].condition = f
    initialize!(mcvb)
    runtraj!(nsteps,dt,integrator,mcvb)
  end

end

"""
Train forces using stochastic gradient descent method with parallelization
"""
# TODO add IO capabilities
function mcvbttrainsgdpar!(nepochs::Int64,totntraj::Int64,nsteps::Int64,dt::Float64,
                           t0::Float64,x0::Vector{Float64},lrf::Float64,lrv::Float64,
                           system::AbstractThermostattedSystem,model::MixedModel,
                           integrator::StochasticEuler,mcvb::MCVBTCallback;
                           restart=true,printevery::Int64=1,Ffile::String="coeffF",
                           Vfile::String="coeffV")

  comm::MPI.Comm = MPI.COMM_WORLD
  rank::Int64 = MPI.Comm_rank(comm)
  nprocs::Int64 = MPI.Comm_size(comm)

  # calculate number of trajectories for this processor
  ntraj::Int64 = floor(Int64,totntraj/nprocs)
  extra::Int64 = totntraj % nprocs
  if (nprocs - rank - 1) < extra
    ntraj += 1
  end

  if rank == 0
    gradf_buf::Vector{Float64} = zeros(Float64, size(mcvb.gradf))
    gradv_buf::Vector{Float64} = zeros(Float64, size(mcvb.gradv))
  end

  # start optimization
  for epoch in 1:nepochs

    # print update
    if rank == 0
      println("# Epoch: $epoch")
      flush(stdout)
    end

    # zero out the averages
    initializeavgs!(mcvb)

    # run trajectories and average gradients
    runtrajs!(ntraj,nsteps,dt,t0,x0,system,integrator,mcvb)

    # average the gradients
    if rank == 0

      # receive information from each sub process and add to average
      for src in 1:(nprocs-1)
        MPI.Recv!(gradf_buf,src,100*src+0,comm)
        mcvb.gradf += gradf_buf
        MPI.Recv!(gradv_buf,src,100*src+1,comm)
        mcvb.gradv += gradv_buf
        mcvb.dkl += MPI.recv(src,100*src+2,comm)[1]
        mcvb.aval += MPI.recv(src,100*src+3,comm)[1]
      end

      # average them over the number of trajectories
      average!(totntraj,mcvb)

      # TODO make nicer with IO
      # write average quantities to stdout
      println(mcvb.dkl," ",mcvb.aval)
      flush(stdout)

      # update coefficients for forces
      updateparams!(lrf,mcvb.gradf,model.potentials[end])

      # update coefficients for valuebaseline
      updateparams!(lrv,mcvb.gradv,mcvb.vbl)

      # TODO file system handling and outpout
      if restart
        writedlm("running_coeffF.txt",model.potentials[end].theta)
        writedlm("running_coeffV.txt",mcvb.vbl.theta)
      end

      if epoch % printevery == 0
        # write updated coefficients
        writedlm(Ffile*"_$epoch.txt",model.potentials[end].theta)
        writedlm(Vfile*"_$epoch.txt",mcvb.vbl.theta)
      end

    else

      # send gradients to main processor
      MPI.Send(mcvb.gradf,0,100*rank+0,comm)
      MPI.Send(mcvb.gradv,0,100*rank+1,comm)
      MPI.send(mcvb.dkl,0,100*rank+2,comm)
      MPI.send(mcvb.aval,0,100*rank+3,comm)

    end

    # broadcast the updated parameters
    MPI.Barrier(comm)
    MPI.Bcast!(model.potentials[end].theta,0,comm)
    MPI.Bcast!(mcvb.vbl.theta,0,comm)
    MPI.Barrier(comm)

    # change the parameters for the model on the other processors
    if rank != 0
      setparams!(model.potentials[end])
      setparams!(mcvb.vbl)
    end

    # synchronize
    MPI.Barrier(comm)

  end

end


text:
using Bridge, Distributions, StaticArrays

PLOT = :pyplot
include("plot.jl")

#diag2(x, y) = SDiagonal(@SVector [x,y])
diag2(x, y) = @SMatrix [ x 0. ; 0. y]

struct FitzHughNagumo  <: ContinuousTimeProcess{SVector{2,Float64}}
    α::Float64
    β::Float64
    γ1::Float64
    γ2::Float64
    ϵ::Float64
    s::Float64
#    σ::SDiagonal{2,Float64}
#    a::SDiagonal{2,Float64}
    σ::SMatrix{2,2,Float64}
    a::SMatrix{2,2,Float64}

    FitzHughNagumo(σ, α = 1/3,  β = 0.08*0.7, γ1 = 0.08, γ2 = 0.08*0.8,  ϵ = 1.,  s = 1.) = new(α, β, γ1, γ2, ϵ, s, σ, σ*σ') #'
end

#Bridge.b(t, x, P::FitzHughNagumo) = @SVector [-P.α*x[1]^3+ P.ϵ*(x[1]-x[2]) + P.s, P.γ1*x[1]- P.γ2*x[2] + P.β]
#phi234(t, x, P::FitzHughNagumo) = @SMatrix [0. -x[1]^3 + x[1]-x[2] 1. ; x[1] 0. 0. ] #γ, ϵ, s
function param(θ, σ)
    β, γ, ϵ = θ
    σ1, σ2 = σ
#   σ,    α, β, γ1,γ2, ϵ,  s
    diag2(σ1,σ2), ϵ, β, γ, 1., ϵ, 0.5*ϵ  # α == ϵ, s = 0.5*ϵ
end

Bridge.bderiv(t,x, P::FitzHughNagumo) = @SMatrix [(P.ϵ - 3*P.α*(x[1]*x[1])) P.γ1 ; -P.ϵ -P.γ2]
Bridge.b(t, x, P::FitzHughNagumo) = @SVector [-P.α*x[1]*x[1]*x[1] + P.ϵ*(x[1]-x[2]) + P.s, P.γ1*x[1]- P.γ2*x[2] + P.β]
#phi1234(t, x, P::FitzHughNagumo) = @SMatrix [0. 0. -x[1]*x[1]*x[1] + x[1]-x[2] 1. ; 1. x[1] 0. 0. ]#β, γ, ϵ, s
#intercept1234(t, x, P::FitzHughNagumo) = @SVector [0, - P.γ2*x[2]]
phi123(t, x, P::FitzHughNagumo) = @SMatrix [0. 0. 1.-x[1]*x[1]*x[1] + x[1]-x[2] ; 1. x[1] 0. ] #β, γ, ϵ
intercept123(t, x, P::FitzHughNagumo) = @SVector [0, - P.γ2*x[2]]

Bridge.σ(t, x, P::FitzHughNagumo) = P.σ
Bridge.a(t, x, P::FitzHughNagumo) = P.a
Bridge.Γ(t, x, P::FitzHughNagumo) = inv(P.a)
Bridge.constdiff(::FitzHughNagumo) = true

function conjugateb(YY, th, xi, P, phif, intc)
        n = length(th)
        G = zeros(n,n)
        mu = zeros(th)
        for Y in YY
            for i in 1:length(Y)-1
                phi = phif(Y[i]..., P)
                Gphi = Bridge.Γ(Y[i]..., P)*phi
                zi = phi'*Gphi
                dy = Y.yy[i+1] - Y.yy[i]
                ds = Y.tt[i+1] - Y.tt[i]
                mu = mu + Vector(Gphi'*(dy - intc(Y[i]..., P)*ds))
                G[:] = G +  Matrix(zi * (Y.tt[i+1]-Y.tt[i]))
            end
        end #for m
        WW = G .+ diagm(xi)
        WL = chol(WW,Val{:L})
        th° = WL'\(randn(n)+WL\mu)
end

mcstart(yy) = (zeros(yy), zeros(yy), 0)
function mcnext(mc, x)
    m, m2, n = mc
    delta = x - m
    n = n + 1
    m = m + delta*(1/n)
    m2 = m2 + map((x,y)->x.*y, delta, x - m)
    m, m2, n
end
function mcbandste(mc)
    m, m2, k = mc
    Q = sqrt(2.)*erfinv(0.95)

    ste = eltype(m)[sqrt(v) for v in m2 * (1/(k - 1))]*sqrt(1/k)
    m -Q*ste, m +Q*ste
end
function mcband(mc)
    m, m2, k = mc
    Q = sqrt(2.)*erfinv(0.95)

    std = eltype(m)[sqrt(v) for v in m2 * (1/(k - 1))]
    m-Q*std, m+Q*std
end

function MyProp(u, v, P, proptype=:mdb)
    if proptype == :aff
        cs = Bridge.CSpline(u[1], v[1],  Bridge.b(u..., P),   Bridge.b(v..., P))
        #println(Bridge.b(u..., P),  Bridge.b(v..., P))
        return BridgeProp(P, u..., v..., P.a, cs)
    elseif proptype == :lin
        error("proptype $proptype not implemented")
        return GuidedProp(P, u..., v..., Pt)
    elseif proptype == :dh
        return DHBridgeProp(P, u..., v...)
    else
        error("proptype $proptype not implemented")
    end
end

function main(propid, m, budget)

############## Configuration ###################################
Random.seed!(10)
K = 1000000 #100000
tim = time()
budget *= 60

simid = 1
#propid = 1
proptype = [:dh, :aff, :aff, :aff, :aff][propid]
eulertype = [:mdb, :eul, :tcs, :theta, :utheta][propid]
simname =["ex", "nex"][simid] * "$proptype$eulertype$m"
println(simname)

INNOS = false  # compute innovations from bridge?
STIME = false # use a grid with smaller steps towards the end-point of a bridge

#θ =  [[0.6, 1.4][simid], 1.5, 10., 0.5*10] # [β, γ, ϵ, s]
θ =  [[0.6, 1.4][simid], 1.5, 10.] # [β, γ, ϵ]
θtrue=copy(θ)
#scaleθ = 0.08*[1., 1., 5., 5.]
scaleθ = 0.08*[1., 1., 5.]
σ = [0.25, 0.2]
#σ = [0.20, 0.15]
scaleσ = [0.02,0.02] #([0.01, 0.01],[0.01, 0.01],[0.01, 0.01])[propid]

σtrue = copy(σ)
Ptrue = FitzHughNagumo(param(θtrue, σtrue)...)

if !INNOS && scaleθ == 0
    error("θ does not get updated")
end

n = 400 # number of segments #evil case: n200 T200; nice case n400, T300

#m = 80 # number of euler steps per segments
mextra = div(2000, m) #factor of extra steps for truth
TT = 300.
dt = TT/n/m
println("Δ $(TT/n) dt $dt dt-sim $(TT/n/m/mextra)")
tt = range(0., stop=TT, length=n*m+1)
tttrue = range(0., stop=TT, length=n*mextra*m+1)
ttf = tt[1:m:end]

uu = @SVector [0., 1.]

r = [(:xrange,(-2,2)), (:yrange,(-1,3))]

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

global Y = euler(uu, sample(tttrue, Wiener{SVector{2,Float64}}()), Ptrue)
global Yobs = Y[1:mextra*m:end] #subsample
assert(endof(Yobs) == n+1)

Y2 = SamplePath(tt, copy(Y.yy[1:mextra:end])) #subsample
normalize(tt) = (tt - tt[1])/(tt[end]-tt[1])
tau(t, tmin, tmax) = tmin + (t-tmin).*(2-(t-tmin)/(tmax-tmin))

θs = Float64[]
iter = 0
if STIME
    tts = [tau(collect(tt[m*(i-1)+1:m*(i)+1]), tt[m*(i-1)+1],tt[m*(i)+1]) for i in 1:n]
else
    tts = [(collect(tt[m*(i-1)+1:m*(i)+1])) for i in 1:n]
end
sss = deepcopy(tts)

BB = SamplePath[Y2[m*(i-1)+1:m*i+1] for i in 1:n]
ZZ = SamplePath[Y2[m*(i-1)+1:m*i+1] for i in 1:n]
for i in 1:n
    BB[i].yy[:] =  map(x -> BB[i].yy[1] + x*(BB[i].yy[end] - BB[i].yy[1]),normalize(tts[i]))
end

BBnew = copy(BB)
BBall = vcat([BB[i][1:end-1] for i in 1:n]...)

#
if PLOT == :winston
        xr = (5,7)
        plot2(Yobs, "o", "o"; xrange=xr, yrange=(-3,3),linewidth=0.1)
        display(oplot2(BBall, "b", "b"; xrange=xr, yrange=(-3,3), linewidth=0.7))
elseif PLOT == :pyplot
#    PyPlot.clf();plot(Y[1:3000], linewidth=0.2, color="b"); plot(Yobs[1:5], linewidth=0, marker="o")
end

################### Prior ###################################
conjθs = [1,2,3] # set of conjugate thetas [β, γ, ϵ=alpha, s]
phi = phi123
intercept = intercept123
#estθ = [false, false, true, false] #params to estimate
estθ = [true, true, true] #params to estimate
estσ = [true, true]
# arbitrary starting values

σ[1] = 0.3
σ[2] = 0.3
#σ = copy(σtrue)
θ = 0.5 + 1.0*rand(length(θ))
#θ = copy(θtrue)

for i in 1:length(θ) # start with truth for parameters not to be estimated
    if !estθ[i]
        θ[i] = θtrue[i]
        scaleθ[i] = 0.
    end
end
for i in 1:length(σ) # start with truth for parameters not to be estimated
    if !estσ[i]
        σ[i] = σtrue[i]
        scaleσ[i] = 0.
    end
end

P = FitzHughNagumo(param(θ, σ)...)

#if PLOT
#    plot2(Y, "r","r" ; yrange=(-3,3),linewidth=0.5)
#    display(oplot2(Y2, "b", "b"; yrange=(-3,3), linewidth=0.7))
#end

# reserve some space
i = 1
ww = Array{SVector{2,Float64},1}(length(m*(i-1)+1:m*(i)+1)) ## checkme!
yy = copy(ww)

# Prior
Alpha = 0.002
Beta= 0.002

piσ²(s2) = pdf(InverseGamma(Alpha,Beta), s2[1])*pdf(InverseGamma(Alpha,Beta), s2[2])
lq(x, si) = Bridge.logpdfnormal(x, si^2*I)
PiError = InverseGamma(Alpha,Beta)

#xi = 1./[50., 50., 50., 50.]
xi = 1./[50., 50., 50.]

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

# Bookkeeping
mkpath(joinpath("output",simname))
try # save cp of this file as documentation
    cp(@__FILE__(), joinpath("output",simname,"$simname.jl"); force=true)
catch
end

open(joinpath("output",simname,"truth.txt"), "w") do f
    println(f, "beta gamma eps s sigma eta")
    println(f, join(round.([θtrue ; σtrue],3)," "))
end

open(joinpath("output",simname,"params.txt"), "w") do f
    println(f, "n beta gamma eps s sigma eta")
end

# initialize
bacc = zeros(Int, n-1)
siacc = 0
thacc = 0

mc = mcstart(vcat([BB[i][1:end-1] for i in 1:n]...).yy )
mcparams = mcstart([θ ; σ])

perf = @timed while true

    #σ[:] = σtrue
    #θ[:] = θtrue

    P = FitzHughNagumo(param(θ, σ)...)
    iter += 1

    for i in 1:n-1
        P° = MyProp(Yobs[i], Yobs[i+1], P, proptype)
        Z° = sample!(SamplePath(tts[i],ww), Wiener{SVector{2,Float64}}())
        if eulertype == :tcs
            B = ubridge!(SamplePath(sss[i], yy), Z°, P°)
        elseif eulertype == :mdb
            B = bridge!(SamplePath(sss[i], yy), Z°, P°, Bridge.mdb!)
        elseif eulertype == :eul
            B = bridge!(SamplePath(sss[i], yy), Z°, P°)
        elseif eulertype == :theta
            B = bridge!(SamplePath(sss[i], yy), Z°, P°, Bridge.thetamethod!)
        elseif eulertype == :utheta
            B = ubridgetheta!(SamplePath(sss[i], yy), Z°, P°)
        elseif eulertype == :heun
            B = bridge!(SamplePath(sss[i], yy), Z°, P°, Bridge.heun!)
        else
            error("$eulertype not implemented")
        end
        if iter == 1
             BB[i] = B
             if !INNOS
                ZZ[i] = copy(Z°)
             end
        end
        if eulertype == :tcs
            llold = ullikelihoodtrapez(BB[i], P°)
            llnew = ullikelihoodtrapez(B, P°)
        else
            llold = llikelihood(BB[i], P°)
            llnew = llikelihood(B, P°)
        end
        if rand() < exp(llnew - llold)
            bacc[i] += 1
            BB[i] = copy(B)
            if !INNOS
                ZZ[i] = copy(Z°)
            end
        end
    end

    if INNOS
        # update theta using X

        if iter % 1 == 1
            BBall = vcat([BB[i][1:end-1] for i in 1:n]...)
            θ° = θ + (2rand(length(θ)) .- 1).*scaleθ
            Pθ = FitzHughNagumo(param(θ, σ)...)
            Pθ° = FitzHughNagumo(param(θ°, σ)...)
            ll = girsanov(BBall, Pθ°, Pθ)
            if rand() < exp(ll)
                θ = θ°
                thacc += 1
            end
        end

        #update conjugate theta
     #   if any(estθ[conjθs])
     #       θ[conjθs[estθ]] = conjugateb(BB, θ[conjθs],xi[conjθs], P, phi, intercept)[estθ[conjθs]]
     #   end
    end
    # update sigma (and theta)
    if iter % 1 == 0
        σ° = σ .* exp.(scaleσ .* randn(length(σ)))
        θ° = θ + (2rand(length(θ)) .- 1).*scaleθ/3
        Pσ = FitzHughNagumo(param(θ, σ)...)
        Pσ° = FitzHughNagumo(param(θ°, σ°)...)
        ll = 0.
        for i in 1:n # reparametrization
            P° = MyProp(Yobs[i], Yobs[i+1], Pσ, proptype)
            P°° = MyProp(Yobs[i], Yobs[i+1], Pσ°, proptype)
            if eulertype == :mdb
                if !INNOS
                    Z = ZZ[i]
                else
                    Z = Bridge.mdbinnovations(BB[i], P°)
                end
                #Z2 = sample(Z.tt, Wiener{SVector{2,Float64}}())
                #Z.yy[:] = sqrt(.9)*Z.yy + sqrt(0.1)*Z2.yy
                BBnew[i] = bridge(Z, P°°, Bridge.mdb!)
                ll += lptilde(P°°) - lptilde(P°) + llikelihood(BBnew[i], P°°) - llikelihood(BB[i], P°)
            elseif eulertype == :tcs
                if !INNOS
                    Z = ZZ[i]
                else
                    Z = Bridge.uinnovations(BB[i], P°)
                end
                BBnew[i] = Bridge.ubridge(Z, P°°)
                ll += lptilde(P°°) - lptilde(P°) + ullikelihoodtrapez(BBnew[i], P°°) - ullikelihoodtrapez(BB[i], P°)
            elseif eulertype == :eul
                if !INNOS
                    Z = ZZ[i]
                else
                    Z = Bridge.innovations(BB[i], P°)
                end
                BBnew[i] = bridge(Z, P°°)
                ll += lptilde(P°°) - lptilde(P°) + llikelihood(BBnew[i], P°°) - llikelihood(BB[i], P°)
            elseif eulertype == :theta
                if !INNOS
                    Z = ZZ[i]
                else
                    Z = thetainnovations(BB[i], P°)
                end
                BBnew[i] = bridge(Z, P°°, Bridge.thetamethod!)
                ll += lptilde(P°°) - lptilde(P°) + llikelihood(BBnew[i], P°°) - llikelihood(BB[i], P°)
             elseif eulertype == :utheta
                Z = Bridge.uinnovations(BB[i], P°)
                BBnew[i] = Bridge.ubridge(Z, P°°)
                ll += lptilde(P°°) - lptilde(P°) + ullikelihood(BBnew[i], P°°) - ullikelihood(BB[i], P°)
            elseif eulertype == :heun
                if !INNOS
                    Z = ZZ[i]
                else
                    Z = heuninnovations(BB[i], P°)
                end
                BBnew[i] = bridge(Z, P°°, Bridge.heun!)
                ll += lptilde(P°°) - lptilde(P°) + llikelihood(BBnew[i], P°°) - llikelihood(BB[i], P°)
            else
                error("$eulertype not implemented")
            end

        end

        if rand() < exp(ll) #* (piσ²(σ°.^2)/piσ²(σ.^2)) * (prod(σ°.^2)/prod(σ.^2))

            σ = σ°
            θ = θ°
            siacc += 1
            for i in 1:n
                BB[i] = BBnew[i]
            end
         #  print("acc")
        end
    end
    open(joinpath("output",simname,"params.txt"), "a") do f; println(f, iter, " ", join(round.([θ ; σ],8)," ")) end
    println(iter, "\t", join(round.([θ./θtrue; σ./σtrue; Inf; 100mean(bacc)/iter; 100minimum(bacc)/iter; Inf; 100thacc/iter; 100siacc/iter; ll  ],3),"\t"))

    BBall = vcat([BB[i][1:end-1] for i in 1:n]...)

    mc = mcnext(mc, BBall.yy)
    mcparams = mcnext(mcparams, [θ ; σ])
    if PLOT == :winston && iter % 10 == 0
        xr = (5,7)

        plot2(Yobs, "o", "o"; xrange=xr, yrange=(-3,3),linewidth=0.1)
        display(oplot2(BBall, "b", "g"; xrange=xr, yrange=(-3,3), linewidth=0.7))
    elseif PLOT == :pyplot && iter % 10 == 1
        #xr = (5,7)

        if iter % 1000 == 1
            PyPlot.clf()
            display(plot(Yobs[1:min(10,n)]; linewidth=0.0, marker="o", color="b"))
        end

        display(plot(BBall[m*1:m*min(9, n-1)], linewidth=0.1, color="b"))
    end
    if iter >= K || abs(time() - tim) > budget
        break
    end
end

open(joinpath("output",simname,"info.txt"), "w") do f
    println(f, "n $n m $m T $TT A $Alpha B $Beta")
    println(f, "Y0 = $uu")
    println(f, "xi = $xi")
    println(f, "acc = ", 100bacc/iter, " (br) ",  100siacc/iter, " (si)",  100thacc/iter, " (th)")
    println(f, perf)
end

if PLOT == :winston

    plot2(Y, "r","b" ; yrange=(-3,3),linewidth=0.5)
    oplot2(Yobs,"+r", "+b")
    savefig(joinpath("output",simname,"truth.pdf"))

    plot2(Y, "r","b" ; yrange=(-3,3),linewidth=0.5)
    oplot2(Yobs,"or","ob";symbolsize=0.3)
    mcb = mcband(mc);
    oplot2(SamplePath(tt, mcb[1]),"r","b";linewidth=0.5)
    oplot2(SamplePath(tt, mcb[2]),"r","b";linewidth=0.5)
    hcat(mcbandste(mcparams)..., [θtrue; σtrue])

    savefig(joinpath("output",simname,"band.pdf"))

    xr = (6,10)
    plot2(Y, "r","b" ; xrange=xr,yrange=(-3,3),linewidth=0.5)
    oplot2(Yobs,"or","ob";xrange=xr,symbolsize=0.3)
    mcb = mcband(mc);
    oplot2(SamplePath(tt, mcb[1]),"r","b";xrange=xr,linewidth=0.5)
    oplot2(SamplePath(tt, mcb[2]),"r","b";xrange=xr,linewidth=0.5)
    savefig(joinpath("output",simname,"bandpart.pdf"))

    plot2(vcat([BB[i] for i in 1:n]...), "r", "b"; yrange=(-3,3), linewidth=0.5)
    oplot2(Yobs,"+r", "+b")
    savefig(joinpath("output",simname,"sample.pdf"))
end
mc, mcparams
end
#include("fitzhugh_nagumo_full.jl"); for i in 1:3, m in [10,20,50]; main(i, m, 0.3); end


text:
using jInv.Mesh
using jInv.Utils
using jInv.ForwardShare
using ShapeReconstructionPaLS.Utils
using ParamLevelSet
using ShapeReconstructionPaLS.ShapeFromSilhouette;
using Statistics
using Distributed
using LinearAlgebra
using SparseArrays
using Test

println("Test SoftMax")
n = 20;
m = 10;
u = rand(20).-0.1;
At = sprand(n,m,4.0/m);
At.nzval[:] .= 1.0;

# d0,Jt = softMaxProjWithSensMat(At,u);
d0,Jt = softMaxProjWithSigmoid(At,u);

du0 = 0.1*randn(size(u));

U = spdiagm(0 => u);
dataMax = maximum(At'*U,dims=2)[:];
dataMax = SigmoidFunc(maximum(At'*U,dims=2)[:]);

r = norm(dataMax .- d0)./norm(d0);
@test r < 0.025
for ii = 1:8
	du = (0.5^ii)*du0;
	ut = u + du;
	# dt, = softMaxProjWithSensMat(At,ut);
	dt, = softMaxProjWithSigmoid(At,ut);
	println("norm(dt-d0): ",norm(dt-d0),", norm(dt - d0 - J0*dm): ",norm(dt - d0 - Jt'*du));
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
export Descent, ADAM, Momentum, Nesterov, RMSProp,
ADAGrad, AdaMax, ADADelta, AMSGrad, NADAM, ADAMW,RADAM, OADAM, AdaBelief

include("ObjectiveType.jl")
include("Gradient.jl")
include("Optimizers.jl")
include("Update.jl")
include("Constraint.jl")


text:
module Plotly

using ..Figures, ..Pages, ..JSON

export Scatter, Layout, Config, newPlot
export Font, Title, ColorBar, Gradient, Marker, Line, Error, Margin, Axis

include("utils.jl"); using .Utilities
include("traceutils.jl"); using .TraceUtilities
include("layoututils.jl"); using .LayoutUtilities
include("config.jl"); using .Configs

include("traces.jl");  using .Traces

include("layouts.jl"); using .Layouts

function newPlot(io::IO,id,traces;layout=default["layout"],config=default["config"])
    isa(traces,Vector) || (traces = [traces])

    data = similar(traces,Dict)
    figure(io,id, Style(
        "width" => string(layout.width.value,"px"),
        "height" => string(layout.height.value+30,"px")))
    for (i,trace) in enumerate(traces)
        data[i] = Utilities.diff(default[trace.type.value],trace)
    end
    layout = Utilities.diff(default["layout"],layout)
    config = Utilities.diff(default["config"],config)
    script = "Plotly.newPlot('$(Figures.current[])',$(json(data)),$(json(layout)),$(json(config)));"
    # println("***********************************************")
    # println(script)
    # println("***********************************************")
    print(io,script)
    return io
end

function react(id,traces)
    typeof(traces) <: Vector || (traces = [traces])

    data = similar(traces,Dict)
    figure(id)
    for (i,trace) in enumerate(traces)
        data[i] = diff(default[trace.type.value],trace)
    end
    Pages.broadcast("script",
        """Plotly.react("$(Figures.current[])",$(json(data)));""")
end

function plot(id,traces)
    typeof(traces) <: Vector || (traces = [traces])

    data = similar(traces,Dict)
    figure(id)
    for (i,trace) in enumerate(traces)
        data[i] = diff(default[trace.type.value],trace)
    end
    Pages.broadcast("script",
        """Plotly.plot("$(Figures.current[])",$(json(data)));""")
end

# function Base.display(::Figures.Display, trace::Trace)
#     println("Try to plot the trace.")
#     # id = isempty(Figures.current[]) ? p.plot.divid : Figures.current[]
#     # Pages.broadcast("script","""Plotly.newPlot("$(id)",$(p.plot.data),$(p.plot.layout),$(PlotlyJS.JSON.json(p.options)));""")
#     # println("""Plotly.newPlot("$(id)",$(p.plot.data),$(p.plot.layout),$(PlotlyJS.JSON.json(p.options)));""")
# end

end


text:
using Test
using IntervalArithmetic
using PolynomialBasis
using CartesianMesh
# using Revise
using HDGElasticity

function allequal(v1,v2)
    return all(v1 .== v2)
end

function allapprox(v1,v2)
    return all(v1 .≈ v2)
end

function distance_function(coords,xc)
    return coords[1,:] .- xc
end

function Base.isequal(c1::HDGElasticity.CellMap,c2::HDGElasticity.CellMap)
    flag = true
    flag = flag && allapprox(c1.xiL,c2.xiL)
    flag = flag && allapprox(c1.xiR,c2.xiR)
    flag = flag && allapprox(c1.xL,c2.xL)
    flag = flag && allapprox(c1.xR,c2.xR)
end

x0 = [0.,0.]
widths = [2.,1.]
nelements = [2,1]
mesh = UniformMesh(x0,widths,nelements)

domain = HDGElasticity.cell_domain(mesh)
testdomain = [IntervalBox(0..1,0..1),IntervalBox(1..2,0..1)]
@test allequal(domain,testdomain)

connectivity = HDGElasticity.cell_connectivity(mesh)
testconn = [[(0,0),(2,4),(0,0),(0,0)],
            [(0,0),(0,0),(0,0),(1,2)]]
@test all([allequal(connectivity[i],testconn[i]) for i = 1:2])

@test HDGElasticity.is_interior_cell(connectivity[1]) == false
@test allequal(HDGElasticity.interior_cells(connectivity),[false,false])

basis = TensorProductBasis(2,1)
poly = InterpolatingPolynomial(1,basis)
NF = HDGElasticity.number_of_basis_functions(basis)
coords = HDGElasticity.nodal_coordinates(mesh,basis)
xc = 0.75
coeffs = reshape(distance_function(coords,xc),NF,:)
cellsign = HDGElasticity.cell_signatures(coeffs,poly)
testcellsign = [0,1]
@test allequal(cellsign,testcellsign)

facemaps = HDGElasticity.reference_cell_facemaps(2)
cellmap = HDGElasticity.CellMap(domain[1])
facescale = HDGElasticity.face_determinant_jacobian(cellmap)
isinteriorcell = HDGElasticity.interior_cells(connectivity)
dgmesh = HDGElasticity.DGMesh(domain,connectivity,isinteriorcell,
    cellsign,facemaps,facescale,cellmap)
@test allequal(dgmesh.domain,domain)
@test all([allequal(dgmesh.connectivity[i],connectivity[i]) for i = 1:2])
@test allequal(dgmesh.cellsign,cellsign)
@test allequal(dgmesh.facescale,facescale)
@test isequal(dgmesh.cellmap,cellmap)

dgmesh = HDGElasticity.DGMesh(mesh,coeffs,poly)
@test allequal(dgmesh.domain,domain)
@test all([allequal(dgmesh.connectivity[i],connectivity[i]) for i = 1:2])
@test allequal(dgmesh.cellsign,cellsign)
@test allequal(dgmesh.facescale,facescale)
@test HDGElasticity.dimension(dgmesh) == 2
cellmap = HDGElasticity.CellMap(dgmesh.domain[1])
@test isequal(cellmap,dgmesh.cellmap)


text:
using SGDOptim
using Base.Test

n = 10
k = 8

# SampleSeq

# case 1: X is a matrix, and y is a vector

X = rand(1:100, (5, n))
y = rand(1:100, n)
ord = randperm(n)[1:k]
ss = collect(sample_seq(X, y, ord))

@test length(ss) == k
for i = 1:k
    j = ord[i]
    @test ss[i] == (X[:,j], y[j])
end

# case 2: both X and Y are matrices

X = rand(1:100, (5, n))
Y = rand(1:100, (3, n))
ord = randperm(n)[1:k]
ss = collect(sample_seq(X, Y, ord))

@test length(ss) == k
for i = 1:k
    j = ord[i]
    @test ss[i] == (X[:, j], Y[:, j])
end

# case 3: mini-batches in natural order

@test SGDOptim.batches(10, 3) == UnitRange{Int}[1:3, 4:6, 7:9, 10:10]
@test SGDOptim.batches(12, 3) == UnitRange{Int}[1:3, 4:6, 7:9, 10:12]

X = rand(1:100, (5, n))
y = rand(1:100, n)

bs = UnitRange{Int}[1:3, 4:6, 7:9, 10:10]
ss = collect(minibatch_seq(X, y, 3))

@test length(ss) == length(bs)
for i = 1:length(bs)
    b = bs[i]
    @test ss[i] == (X[:, b], y[b])
end

ord = [4, 2, 3, 1, 2]
ss = collect(minibatch_seq(X, y, 3, ord))

@test length(ss) == length(ord)
for i = 1:length(ord)
    b = bs[ord[i]]
    @test ss[i] == (X[:, b], y[b])
end


text:
module Birtwistle2007

using PyPlot

export runSim

include("model/simulation.jl")
using .Sim

include("plotFunc.jl");

function runSim()
    plotFunc_timecourse(Sim)
end

end  # module


text:
@testset "PermuteDims" begin
    @test (PermuteDims <: Augmentor.AffineOperation) == false
    @test (PermuteDims <: Augmentor.ImageOperation) == false
    @test (PermuteDims <: Augmentor.Operation) == true

    @testset "constructor" begin
        @test_throws MethodError PermuteDims()
        @test_throws MethodError PermuteDims(())
        @test_throws MethodError PermuteDims(Float64)
        @test_throws ArgumentError PermuteDims((2,2))
        @test_throws ArgumentError PermuteDims(2,2)
        @test typeof(PermuteDims(1,2)) <: PermuteDims{2} <: Augmentor.Operation
        @test typeof(PermuteDims((1,2))) <: PermuteDims{2} <: Augmentor.Operation
        @test typeof(PermuteDims(1)) <: PermuteDims{1} <: Augmentor.Operation
        @test typeof(PermuteDims((1,))) <: PermuteDims{1} <: Augmentor.Operation
        @test typeof(PermuteDims((1,2))) <: PermuteDims{2} <: Augmentor.Operation
        @test typeof(PermuteDims((3,1,2))) <: PermuteDims{3} <: Augmentor.Operation
        @test str_show(PermuteDims((1,))) == "Augmentor.PermuteDims(1)"
        @test str_show(PermuteDims((1,2))) == "Augmentor.PermuteDims(1, 2)"
        @test str_showconst(PermuteDims((1,3,2))) == "PermuteDims(1, 3, 2)"
        @test str_showcompact(PermuteDims((3,2,1))) == "Permute dimension order to (3, 2, 1)"
    end
    @testset "eager" begin
        @test Augmentor.supports_eager(PermuteDims) === true
        @test @inferred(Augmentor.supports_eager(typeof(PermuteDims(2,1)))) === true
        f = (img) -> permutedims(img, (2,1))
        imgs = [
            (rect),
            (Augmentor.prepareaffine(rect)),
            (OffsetArray(rect, -2, -1)),
            (view(rect, IdentityRange(1:2), IdentityRange(1:3))),
        ]
        @testset "single image" begin
            for img_in in imgs
                img_out = f(img_in)
                @test_throws MethodError Augmentor.applyeager(PermuteDims(3,2,1), img_in)
                res = @inferred(Augmentor.applyeager(PermuteDims(2,1), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
        @testset "multiple images" begin
            for img_in1 in imgs, img_in2 in imgs
                img_out1 = f(img_in1)
                img_out2 = f(img_in2)
                img_in = (img_in1, img_in2)
                img_out = (img_out1, img_out2)
                @test_throws MethodError Augmentor.applyeager(PermuteDims(3,2,1), img_in)
                res = @inferred(Augmentor.applyeager(PermuteDims(2,1), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(PermuteDims) === false
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(PermuteDims) === false
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(PermuteDims) === true
        @test @inferred(Augmentor.supports_lazy(typeof(PermuteDims(2,1)))) === true
        f = (img) -> permuteddimsview(img, (2,1))
        imgs = [
            (rect),
            (Augmentor.prepareaffine(rect)),
            (OffsetArray(rect, -2, -1)),
            (view(rect, IdentityRange(1:2), IdentityRange(1:3))),
        ]
        @testset "single image" begin
            for img_in in imgs
                img_out = f(img_in)
                @test_throws MethodError Augmentor.applylazy(PermuteDims(3,2,1), img_in)
                res = @inferred(Augmentor.applylazy(PermuteDims(2,1), img_in))
                @test res == img_out
                @test res == Augmentor.applyeager(PermuteDims(2,1), img_in)
                @test typeof(res) == typeof(img_out)
            end
        end
        @testset "multiple images" begin
            for img_in1 in imgs, img_in2 in imgs
                img_out1 = f(img_in1)
                img_out2 = f(img_in2)
                img_in = (img_in1, img_in2)
                img_out = (img_out1, img_out2)
                @test_throws MethodError Augmentor.applylazy(PermuteDims(3,2,1), img_in)
                res = @inferred(Augmentor.applylazy(PermuteDims(2,1), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
    end
    @testset "view" begin
        @test Augmentor.supports_view(PermuteDims) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(PermuteDims) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(PermuteDims) === false
    end
end

# --------------------------------------------------------------------

@testset "Reshape" begin
    @test (Reshape <: Augmentor.AffineOperation) == false
    @test (Reshape <: Augmentor.ImageOperation) == false
    @test (Reshape <: Augmentor.Operation) == true

    @testset "constructor" begin
        @test_throws MethodError Reshape()
        @test_throws MethodError Reshape(())
        @test_throws MethodError Reshape(Float64)
        @test_throws MethodError Reshape(1.,2.)
        @test typeof(Reshape(2,2)) <: Reshape{2} <: Augmentor.Operation
        @test typeof(Reshape(1,2)) <: Reshape{2} <: Augmentor.Operation
        @test typeof(Reshape((1,2))) <: Reshape{2} <: Augmentor.Operation
        @test typeof(Reshape(1)) <: Reshape{1} <: Augmentor.Operation
        @test typeof(Reshape((1,))) <: Reshape{1} <: Augmentor.Operation
        @test typeof(Reshape((1,2))) <: Reshape{2} <: Augmentor.Operation
        @test typeof(Reshape((3,1,2))) <: Reshape{3} <: Augmentor.Operation
        @test str_show(Reshape((1,))) == "Augmentor.Reshape{1}(1)"
        @test str_show(Reshape((1,2))) == "Augmentor.Reshape{2}(1, 2)"
        @test str_showconst(Reshape((1,3,2))) == "Reshape(1, 3, 2)"
        @test str_showcompact(Reshape((3,2,1))) == "Reshape array to 3×2×1"
        @test str_showcompact(Reshape(10)) == "Reshape array to 10-element vector"
    end
    @testset "eager" begin
        @test Augmentor.supports_eager(Reshape) === false
        @test @inferred(Augmentor.supports_eager(typeof(Reshape(2,1)))) === false
        imgo = Augmentor.plain_array(reshape(rect, (3,2,1)))
        imgs = [
            (rect),
            (Augmentor.prepareaffine(rect)),
            (OffsetArray(rect, -2, -1)),
            (view(rect, IdentityRange(1:2), IdentityRange(1:3))),
        ]
        @testset "single image" begin
            for img_in in imgs
                img_out = imgo
                res = @inferred(Augmentor.applyeager(Reshape(3,2,1), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
        @testset "multiple images" begin
            for img_in1 in imgs, img_in2 in imgs
                img_out1 = imgo
                img_out2 = imgo
                img_in = (img_in1, img_in2)
                img_out = (img_out1, img_out2)
                res = @inferred(Augmentor.applyeager(Reshape(3,2,1), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(Reshape) === false
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(Reshape) === false
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(Reshape) === true
        @test @inferred(Augmentor.supports_lazy(typeof(Reshape(2,1)))) === true
        imgs = [
            (rect),
            (Augmentor.prepareaffine(rect)),
            (OffsetArray(rect, -2, -1)),
            (view(rect, IdentityRange(1:2), IdentityRange(1:3))),
        ]
        @testset "single image" begin
            for img_in in imgs
                img_out = reshape(Augmentor.plain_axes(img_in), (3,2,1))
                res = @inferred(Augmentor.applylazy(Reshape(3,2,1), img_in))
                @test res == Augmentor.applyeager(Reshape(3,2,1), img_in)
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
        @testset "multiple images" begin
            for img_in1 in imgs, img_in2 in imgs
                img_out1 = reshape(Augmentor.plain_axes(img_in1), (3,2,1))
                img_out2 = reshape(Augmentor.plain_axes(img_in2), (3,2,1))
                img_in = (img_in1, img_in2)
                img_out = (img_out1, img_out2)
                res = @inferred(Augmentor.applylazy(Reshape(3,2,1), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
    end
    @testset "view" begin
        @test Augmentor.supports_view(Reshape) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(Reshape) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(Reshape) === false
    end
end


text:
file_dir = Pkg.dir("D4M")*"/examples/2Apps/1EntityAnalysis/Entity.jld";

using JLD

E = load(file_dir)["E"]

E = logical(E);

A = sqIn(E);

d = diag(Adj(A))

A = putAdj(A,Adj(A)-diagm(d));

####

i,j,v = findnz(Adj(A));

An = putAdj(A, sparse(i,j,v ./ min(d[i],d[j])));

####

x = "LOCATION/new york,";

p = StartsWith("PERSON/,");

print( (A[p,x] > 4) & (An[p,x] > 0.3));

####

p0 = "PERSON/john kennedy,";

p1 = Row(A[p,p0] + A[p0,p]);

using PyPlot

spy(A[p1,p1])

p2 = Row( A[p1,p1] - (A[p,p0]+ A[p0,p]));

A[p2,p2] > 1


text:
const VALIDATION_DEBUG_SIGNAL = Ref{Bool}(false)

function debug_callback(args...)
    VALIDATION_DEBUG_SIGNAL[] = true
    default_debug_callback(args...)
end

const debug_callback_c = @cfunction(debug_callback, UInt32, (VkDebugUtilsMessageSeverityFlagBitsEXT, VkDebugUtilsMessageTypeFlagBitsEXT, Ptr{vk.VkDebugUtilsMessengerCallbackDataEXT}, Ptr{Cvoid}))
const API_VERSION = v"1.2"
const VALIDATION_LAYER = "VK_LAYER_KHRONOS_validation"
const WITH_VALIDATION = Ref(false)
const INSTANCE_LAYERS = [
]
const INSTANCE_EXTENSIONS = [
]
const DEVICE_EXTENSIONS = [
]
const ENABLED_FEATURES = PhysicalDeviceFeatures(
)

let available_layers = unwrap(enumerate_instance_layer_properties())
    if VALIDATION_LAYER ∈ getproperty.(available_layers, :layer_name)
        push!(INSTANCE_LAYERS, VALIDATION_LAYER)
        push!(INSTANCE_EXTENSIONS, "VK_EXT_debug_utils")
        WITH_VALIDATION[] = true
    else
        @warn "Validation layer not found."
    end
end

@testset "Vulkan tests" begin
    include("init.jl")

    @testset "Utilities" begin
        @testset "Function pointers" begin
            @test function_pointer(instance, "vkEnumeratePhysicalDevices") == get_instance_proc_addr("vkEnumeratePhysicalDevices"; instance)
            @test function_pointer(device, "vkQueueWaitIdle") == get_device_proc_addr(device, "vkQueueWaitIdle")
        end
    end

    @testset "Debugging" begin
        @testset "Validation" begin
            if WITH_VALIDATION[]
                @test VALIDATION_DEBUG_SIGNAL[]
            end
        end
    end

    @testset "Bitmask flags" begin
        buffer = Buffer(device, 24, BUFFER_USAGE_INDIRECT_BUFFER_BIT, VK_SHARING_MODE_EXCLUSIVE, [0]; flags=BUFFER_CREATE_PROTECTED_BIT & BUFFER_CREATE_SPARSE_ALIASED_BIT)
        buffer_2 = Buffer(device, 24, BUFFER_USAGE_SHADER_DEVICE_ADDRESS_BIT, VK_SHARING_MODE_EXCLUSIVE, [0]; flags=0)
        @test buffer isa Buffer
        @test buffer_2 isa Buffer
    end

    @testset "Buffers" begin
        include("buffers.jl")
    end
end

GC.gc()


text:
nt = 1000
seed = 99
ssm = acqr(0.8,1.5,2,3,0.5,2)
m = acqr_simula(ssm,nt,seed)
y = m["y"]

#######################
print("Testing acqr_em_starting: ")
ssm = acqr_em_starting(y)
if sum(ssm.c^2) > 1e-6
    println("OK")
end


text:
using Distributed
using DifferentialEquations
using Plots

addprocs()
@everywhere using DifferentialEquations

function f(du,u,p,t)
  du[1] = p[1] * u[1] - p[2] * u[1]*u[2]
  du[2] = -3 * u[2] + u[1]*u[2]
end

function g(du,u,p,t)
    du[1] = p[3]*u[1]
    du[2] = p[4]*u[2]
  end

  p = [1.5,1.0,0.1,0.1]
prob = SDEProblem(f,g,[1.0,1.0],(0.0,10.0),p)

function prob_func(prob,i,repeat)
    x = 0.3rand(2)
    # remake(prob,p=[p[1:2];x])
    remake(prob)
  end

  ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
sim = solve(ensemble_prob,SRIW1(),trajectories=10)
using Plots; plotly()
using Plots; plot(sim,linealpha=0.6,color=:blue,vars=(0,1),title="Phase Space Plot")
plot!(sim,linealpha=0.6,color=:red,vars=(0,2),title="Phase Space Plot")

summ = EnsembleSummary(sim,0:0.1:10)
pyplot() # Note that plotly does not support ribbon plots
plot(summ,fillalpha=0.5)


text: include("./TikzPlotScatter.jl")


text:
"""
Clamped cylinder, transient vibration.

The cylinder is banged at the initial time (i.e. it is given an initial velocity).
"""
module test_loops

using LinearAlgebra
using SparseArrays
using Arpack
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: default_layout_3d, plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using PlotlyJS
using Gnuplot
using FinEtools.MeshExportModule.VTKWrite: vtkwritecollection

using BenchmarkTools
using ProfileView

E = 200e9;
nu = 0.3;
rho = 7800.0
R = 0.1
L = 0.8
ksi = 0.0
omegad = 2*pi*1000
color = "blue"
omegaf = 2*pi*10^5
tend = 15 * (2*pi/omegaf)

function loop!(M, K, ksi, U0, V0, tend, dt, force!, peek)
    U1 = deepcopy(U0)
    V1 = deepcopy(U0)
    A0 = deepcopy(U0)
    A1 = deepcopy(U0)
    F = deepcopy(U0)
    E = deepcopy(U0)
    Vp = deepcopy(U0)
    C = deepcopy(U0)
    C .= (ksi*2*omegad) .* vec(diag(M))
    invMC = deepcopy(U0)
    invMC .= 1.0 ./ (vec(diag(M)) .+ (dt/2) .* C)
    nsteps = Int64(round(tend/dt))
    if nsteps*dt < tend
        dt = tend / (nsteps+1)
    end
    t = 0.0

    # @btime @. $U1 = $U0 + $dt*$V0 + (($dt^2)/2)*$A0; #1
    # @btime $force!($F, $t); #2
    # @btime $F .-=  mul!($E, $K, $U1) #3
    # @btime @. $Vp = $V0 + ($dt/2) * $A0 #4
    # @btime @. $F -= $C * $Vp #5
    # @btime $A1 .= $invMC .* $F #6
    # @btime @. $A1 = $invMC * $F #6
    # @btime @inbounds for i in eachindex($A1)
    #     $A1[i] = $invMC[i] * $F[i] #6
    # end
    # @btime @. $V1 = $V0 + ($dt/2)* ($A0+$A1); #7
    # @btime $peek(1, $U0, $t) #8

    A0 .= invMC .* force!(F, t);
    peek(0, U0, t)
    @time for step in 1:nsteps
        # Displacement update
        @. U1 = U0 + dt*V0 + ((dt^2)/2)*A0;
        # External loading
        force!(F, t);
        # Add elastic restoring forces
        F .-=  mul!(E, K, U1)
        # Add damping forces
        @. Vp = V0 + (dt/2) * A0
        @. F -= C * Vp
        # Compute the new acceleration.
        A1 .= invMC .* F
        # # Update the velocity
        @. V1 = V0 + (dt/2)* (A0+A1);
        t = t + dt
        U0, U1 = U1, U0
        V0, V1 = V1, V0
        A0, A1 = A1, A0
        peek(step, U0, t)
    end
end

cylindrical!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) = begin
    r = vec(XYZ); r[2] = 0.0;
    csmatout[:, 3] .= vec(r)/norm(vec(r))
    csmatout[:, 2] .= (0.0, 1.0, 0.0) #  this is along the axis
    cross3!(view(csmatout, :, 1), view(csmatout, :, 2), view(csmatout, :, 3))
    return csmatout
end

struct ThreadBuffer
    # colrange::UnitRange{Int64}
    kcolumns::SparseMatrixCSC{Float64, Int64}
    uv::SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}
    result::Vector{Float64}
end

struct ThreadBuff
    colrange::UnitRange{Int64}
    result::Vector{Float64}
end

function _execute(n = 8, thickness = 0.01, visualize = true)

    # @info "Mesh: $n elements per side"
    # Mesh
    tolerance = min(R, L) / n  / 100

    tolerance = R/n/1000
    fens, fes = T3block(360.0, L, n, 4*n);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a/180*pi), y-L/2, R*cos(a/180*pi))
    end
    bfes = meshboundary(fes)
    candidates = connectednodes(bfes)
    fens, fes = mergenodes(fens, fes, tolerance, candidates)
    bfes = meshboundary(fes)
    @info "Mesh $(count(fens)) nodes, $(count(fes)) elements"

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
    ocsys = CSys(3, 3, cylindrical!)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = FEMMShellT3FFModule.make(IntegDomain(fes, TriRule(1), thickness), ocsys, mater)
    # Set up
    femm.drilling_mass_scale = 1.0
    femm.drilling_stiffness_scale = 1.0

# Construct the requisite fields, geometry and displacement
# Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

# Apply EBC's
    l1 = connectednodes(meshboundary(fes))
    for i in 1:6
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

# Assemble the system matrix
    FEMMShellT3FFModule.associategeometry!(femm, geom0)
    K = FEMMShellT3FFModule.stiffness(femm, geom0, u0, Rfield0, dchi);
    M = FEMMShellT3FFModule.mass(femm, SysmatAssemblerSparseDiag(), geom0, dchi);

    U0 = gathersysvec(dchi)
    R0 = deepcopy(U0)

    @show nth = Base.Threads.nthreads()
    chunk = Int(floor(dchi.nfreedofs / nth))

    threadbuffers = ThreadBuffer[];
    for th in 1:nth
        # This thread will work on a subset of the columns.
        colrange = th < nth ? (chunk*(th-1)+1:chunk*(th)+1-1) : (chunk*(th-1)+1:dchi.nfreedofs)
        push!(threadbuffers, ThreadBuffer(K[:, colrange], view(U0, colrange), deepcopy(R0)));
    end

    function parmul!(R, tbs)
        tasks = [];
        for th in 1:length(tbs)
            tb = tbs[th]
            push!(tasks, Threads.@spawn begin
                mul!(tb.result, tb.kcolumns, tb.uv)
            end);
        end
        Threads.wait(tasks[1]);
        R .= tbs[1].result
        for th in 2:length(tbs)
            Threads.wait(tasks[th]);
            R .+= tbs[th].result
        end
        R
    end

    threadbuffs = ThreadBuff[];
    for th in 1:nth
        colrange = th < nth ? (chunk*(th-1)+1:chunk*(th)+1-1) : (chunk*(th-1)+1:dchi.nfreedofs)
        push!(threadbuffs, ThreadBuff(colrange, deepcopy(R0)));
    end

    function parcoomul!(R, IB, JB, VB, tbs, U)
        function innercoomul!(A, IB, JB, VB, C)
            nz = length(IB)
            @assert length(JB) == nz
            @assert length(VB) == nz
            @inbounds @simd for i in 1:nz
                ir = IB[i]
                jc = JB[i]
                A[ir] += VB[i] * C[jc]
            end
            nothing
        end
        tasks = [];
        for th in 1:length(tbs)
            tb = tbs[th]
            push!(tasks, Threads.@spawn begin
                tb.result .= zero(eltype(tb.result))
                innercoomul!(tb.result, view(IB, tb.colrange), view(JB, tb.colrange), view(VB, tb.colrange), U)
            end);
        end
        Threads.wait(tasks[1]);
        R .= tbs[1].result
        for th in 2:length(tbs)
            Threads.wait(tasks[th]);
            R .+= tbs[th].result
        end
        R
    end

    function coomul!(A, IB, JB, VB, C)
        nz = length(IB)
        @assert length(JB) == nz
        @assert length(VB) == nz
        A .= zero(eltype(A))
        @inbounds @simd for i in 1:nz
            ir = IB[i]
            jc = JB[i]
            A[ir] += VB[i] * C[jc]
        end
    end

    IK, JK, VK = findnz(K)

    oU0 = deepcopy(U0)
    oU0 .= rand()

    U0 .= oU0
    mul!(R0, K, U0)
    R01 = deepcopy(R0)

    U0 .= oU0
    parmul!(R0, threadbuffers)
    @show norm(R0 - R01)/norm(R01)

    U0 .= oU0
    parcoomul!(R0, IK, JK, VK, threadbuffs, U0)
    @show norm(R0 - R01)/norm(R01)

    U0 .= oU0
    coomul!(R0, IK, JK, VK, U0)
    @show norm(R0 - R01)/norm(R01)

    @btime mul!($R0, $K, $U0)
    @btime $coomul!($R0, $IK, $JK, $VK, $U0)
    @btime $parmul!($R0, $threadbuffers)
    @btime $parcoomul!($R0, $IK, $JK, $VK, $threadbuffs, $U0)

    # oR0 = deepcopy(R0)
    # oR0 .= rand()
    # oU0 = deepcopy(U0)
    # oU0 .= rand()

    # R0 .= oR0; U0 .= oU0
    # @. R0 = U0 + 0.1*R0
    # R01 = deepcopy(R0)

    # R0 .= oR0; U0 .= oU0
    # @inbounds @simd for i in eachindex(U0)
    #     R0[i] = U0[i] + 0.1*R0[i]
    # end
    # @show norm(R0 - R01)

    # R0 .= oR0; U0 .= oU0
    # Threads.@threads for i in eachindex(U0)
    #     @inbounds R0[i] = U0[i] + 0.1*R0[i]
    # end
    # @show norm(R0 - R01)

    # ranges = []
    # for th in 1:nth
    #     # This thread will work on a subset of the columns.
    #     colrange = th < nth ? (chunk*(th-1)+1:chunk*(th)+1-1) : (chunk*(th-1)+1:dchi.nfreedofs)
    #     push!(ranges, colrange)
    # end
    # R0 .= oR0; U0 .= oU0
    # Threads.@threads for r in ranges
    #     @inbounds for j in r
    #         R0[j] = U0[j] + 0.1*R0[j]
    #     end
    # end
    # @show norm(R0 - R01)

    # @btime @. $R0 = $U0 + 0.1*$R0
    # @btime for i in eachindex($U0)
    #     $R0[i] = $U0[i] + 0.1*$R0[i]
    # end
    # @btime @inbounds @simd for i in eachindex($U0)
    #     $R0[i] = $U0[i] + 0.1*$R0[i]
    # end
    # @btime Threads.@threads for i in eachindex($U0)
    #     @inbounds $R0[i] = $U0[i] + 0.1*$R0[i]
    # end
    # ranges = []
    # for th in 1:2
    #     # This thread will work on a subset of the columns.
    #     colrange = th < nth ? (chunk*(th-1)+1:chunk*(th)+1-1) : (chunk*(th-1)+1:dchi.nfreedofs)
    #     push!(ranges, colrange)
    # end
    # @btime Threads.@threads for r in $ranges
    #     @inbounds for j in r
    #         $R0[j] = $U0[j] + 0.1*$R0[j]
    #     end
    # end

    nothing
end

function test(ns = [64])
    @info "Clamped cylinder vibration"
    for n in ns
        _execute(n, 0.01, !false)
    end
    return true
end

function allrun()
    println("#####################################################")
    println("# test_convergence ")
    test()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing

 using .Main.test_loops; Main.test_loops.test(64*[2, 3, 4, 5])