text:
module VariationalInference

include("../src/Jaynes.jl")
using .Jaynes
using Gen

model = @jaynes (a₀::Float64, b₀::Float64, μ₀::Float64, λ₀::Float64, N::Int) -> begin
    τ ~ Gamma(a₀, b₀)
    μ ~ Normal(μ₀, 1 / (λ₀ * τ))
    x = [{:x => i} ~ Normal(μ, 1 / τ) for i in 1 : N]
    x
end

# Factorized approximation.
variational_family = @jaynes () -> begin
    α, β = learnable(:α), learnable(:β)
    τ ~ Gamma(α^2 + 0.01, β^2 + 0.01)
    μp, logσ = learnable(:μp), learnable(:logσ)
    μ ~ Normal(μp, exp(logσ))
end
init_param!(variational_family, [(:α, ) => 1.0, (:β, ) => 1.0,
                                 (:μp, ) => 0.0, (:logσ, ) => 0.0])

# Black box variational inference.
N = 100
observations = begin
    a₀ = 0.5
    b₀ = 1.0
    μ₀ = 0.0
    λ₀ = 3.0
    tr = simulate(model, (a₀, b₀, μ₀, λ₀, N))
    filter(k -> k[1] isa Pair, tr)
end
display(observations)

update = ParamUpdate(GradientDescent(1e-5, 100),
                     variational_family => [(:α, ), (:β, ), (:μp, ), (:logσ, )])

@time elbo_estimate, traces, elbo_history = black_box_vi!(model, (0.5, 1.0, 0.0, 3.0, N), observations, variational_family, (), update; iters = 1000, verbose = true)
display(get_params(variational_family))

end # module


text:
module PassThrough

using PCLLite
using Cxx
using Base.Test

cxx"""
#include <iostream>
#include <pcl/point_types.h>
#include <pcl/filters/passthrough.h>
"""

cloud = icxx"""
    pcl::PointCloud<pcl::PointXYZ>::Ptr(new pcl::PointCloud<pcl::PointXYZ>);
    """
cloud_filtered = icxx"""
    pcl::PointCloud<pcl::PointXYZ>::Ptr(new pcl::PointCloud<pcl::PointXYZ>);
"""

# fill in the cloud data
icxx"$cloud->width = 5;"
icxx"$cloud->height = 1;"
icxx"$cloud->points.resize($cloud->width * $cloud->height);"

size = convert(Int, icxx"$cloud->points.size();")
for i in 0:size-1
    icxx"""
    $cloud->points[$i].x = 1024 * rand () / (RAND_MAX + 1.0f);
    $cloud->points[$i].y = 1024 * rand () / (RAND_MAX + 1.0f);
    $cloud->points[$i].z = 1024 * rand () / (RAND_MAX + 1.0f);
    """
end

icxx"""
  std::cerr << "Cloud before filtering: " << std::endl;
  for (size_t i = 0; i < $cloud->points.size (); ++i)
    std::cerr << "    " << $cloud->points[i].x << " "
                        << $cloud->points[i].y << " "
                        << $cloud->points[i].z << std::endl;
"""

# Create the filtering object
pass = icxx"pcl::PassThrough<pcl::PointXYZ>();"
icxx"""
$pass.setInputCloud($cloud);
$pass.setFilterFieldName("z");
$pass.setFilterLimits(0.0, 1.0);
$pass.filter(*($cloud_filtered));
"""

icxx"""
  std::cerr << "Cloud after filtering: " << std::endl;
  for (size_t i = 0; i < $cloud_filtered->points.size (); ++i)
    std::cerr << "    " << $cloud_filtered->points[i].x << " "
                        << $cloud_filtered->points[i].y << " "
                        << $cloud_filtered->points[i].z << std::endl;
"""

Base.length{T}(cloud::cxxt"boost::shared_ptr<pcl::PointCloud<$T>>") =
    icxx"$cloud->points.size();"
@test length(cloud_filtered) <= length(cloud)

end # module PassThrough

ext: jl


text:
using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.DeprecatedTest
const MOIU = MOI.Utilities

@testset "Start value attributes and nothing" begin
    model = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
    MOI.DeprecatedTest.start_values_unset_test(
        model,
        MOI.DeprecatedTest.Config(),
    )
end

@testset "Failed copy" begin
    include("../dummy.jl")
    model = DummyModelWithAdd()
    MOIT.failcopytestc(model)
    MOIT.failcopytestia(model)
    MOIT.failcopytestva(model)
    MOIT.failcopytestca(model)
end


text:
# Algorithms for Best Arm Identification in Exponential Family Bandit Models in the Fixed-Confidence Setting

# The nature of the distribution should be precised by choosing a value for typeDistribution before including the current file

# All the algorithms take the following input
# mu : vector of arms means
# delta : risk level
# rate : the exploration rate (a function)

using Distributions
using PyPlot

include("KLfunctions.jl")

if (typeDistribution=="Bernoulli")
   d=dBernoulli
   dup=dupBernoulli
   dlow=dlowBernoulli
   function sample(mu)
       (rand()<mu)
   end
   function bdot(theta)
        exp(theta)/(1+exp(theta))
   end
   function bdotinv(mu)
        log(mu/(1-mu))
   end
elseif (typeDistribution=="Poisson")
   d=dPoisson
   dup=dupPoisson
   dlow=dlowPoisson
   function sample(mu)
       rand(Poisson(mu))
   end
   function bdot(theta)
        exp(theta)
   end
   function bdotinv(mu)
        log(mu)
   end
elseif (typeDistribution=="Exponential")
   d=dExpo
   dup=dupExpo
   dlow=dlowExpo
   function sample(mu)
       -mu*log(rand())
   end
   function bdot(theta)
        -log(-theta)
   end
   function bdotinv(mu)
        -exp(-mu)
   end
elseif (typeDistribution=="Gaussian")
   # sigma (std) must be defined !
   d=dGaussian
   dup=dupGaussian
   dlow=dlowGaussian
   function sample(mu)
       mu+sigma*randn()
   end
   function bdot(theta)
        sigma^2*theta
   end
   function bdotinv(mu)
        mu/sigma^2
   end
end

# COMPUTING THE OPTIMAL WEIGHTS

function dicoSolve(f, xMin, xMax, delta=1e-11)
  # find m such that f(m)=0 using dichotomix search
  l = xMin
  u = xMax
  sgn = f(xMin)
  while u-l>delta
    m = (u+l)/2
    if f(m)*sgn>0
      l = m
    else
      u = m
    end
  end
  m = (u+l)/2
  return m
end

function I(alpha,mu1,mu2)
    if (alpha==0)|(alpha==1)
       return 0
    else
        mid=alpha*mu1 + (1-alpha)*mu2
        return alpha*d(mu1,mid)+(1-alpha)*d(mu2,mid)
    end
end

muddle(mu1, mu2, nu1, nu2) = (nu1*mu1 + nu2*mu2)/(nu1+nu2)

function cost(mu1, mu2, nu1, nu2)
  if (nu1==0)&(nu2==0)
     return 0
  else
     alpha=nu1/(nu1+nu2)
     return((nu1 + nu2)*I(alpha,mu1,mu2))
  end
end

function xkofy(y, k, mu, delta = 1e-11)
  # return x_k(y), i.e. finds x such that g_k(x)=y
  g(x)=(1+x)*cost(mu[1], mu[k], 1/(1+x), x/(1+x))-y
  xMax=1
  while g(xMax)<0
       xMax=2*xMax
  end
  return dicoSolve(x->g(x), 0, xMax, 1e-11)
end

function aux(y,mu)
  # returns F_mu(y) - 1
  K = length(mu)
  x = [xkofy(y, k, mu) for k in 2:K]
  m = [muddle(mu[1], mu[k], 1, x[k-1]) for k in 2:K]
  return (sum([d(mu[1],m[k-1])/(d(mu[k], m[k-1])) for k in 2:K])-1)
end

function oneStepOpt(mu, delta = 1e-11)
  yMax=0.5
  if d(mu[1], mu[2])==Inf
     # find yMax such that aux(yMax,mu)>0
     while aux(yMax,mu)<0
          yMax=yMax*2
     end
  else
     yMax=d(mu[1],mu[2])
  end
  y = dicoSolve(y->aux(y, mu), 0, yMax, delta)
  x =[xkofy(y, k, mu, delta) for k in 2:length(mu)]
  unshift!(x, 1)
  nuOpt = x/sum(x)
  return nuOpt[1]*y, nuOpt
end

function OptimalWeights(mu, delta=1e-11)
  # returns T*(mu) and w*(mu)
  K=length(mu)
  IndMax=find(mu.==maximum(mu))
  L=length(IndMax)
  if (L>1)
     # multiple optimal arms
     vOpt=zeros(1,K)
     vOpt[IndMax]=1/L
     return 0,vOpt
  else
     mu=vec(mu)
     index=sortperm(mu,rev=true)
     mu=mu[index]
     unsorted=vec(collect(1:K))
     invindex=zeros(Int,K)
     invindex[index]=unsorted
     # one-step optim
     vOpt,NuOpt=oneStepOpt(mu,delta)
     # back to good ordering
     nuOpt=NuOpt[invindex]
     NuOpt=zeros(1,K)
     NuOpt[1,:]=nuOpt
     return vOpt,NuOpt
  end
end

# OPTIMAL ALGORITHMS

function TrackAndStop(mu,delta,rate)
  condition = true
  K=length(mu)
  N = zeros(1,K)
  S = zeros(1,K)
  # initialization
  for a in 1:K
      N[a]=1
      S[a]=sample(mu[a])
  end
  t=K
  Best=1
  while (condition)
       Mu=S./N
       # Empirical best arm
       IndMax=find(Mu.==maximum(Mu))
       Best=IndMax[floor(Int,length(IndMax)*rand())+1]
       I=1
       if (length(IndMax)>1)
          # if multiple maxima, draw one them at random
          I = Best
       else
       	  # compute the stopping statistic
       	  NB=N[Best]
          SB=S[Best]
          muB=SB/NB
          MuMid=(SB+S)./(NB+N)
          Index=collect(1:K)
          splice!(Index,Best)
          Score=minimum([NB*d(muB,MuMid[i])+N[i]*d(Mu[i],MuMid[i]) for i in Index])
          if (Score > rate(t,0,delta))
            # stop
            condition=false
          elseif (t >10000000)
            # stop and outputs (0,0)
            condition=false
            Best=0
            print(N)
            print(S)
            N=zeros(1,K)
          else
            if (minimum(N) <= max(sqrt(t) - K/2,0))
               # forced exploration
               I=indmin(N)
            else
               # continue and sample an arm
	       val,Dist=OptimalWeights(Mu,1e-11)
               # choice of the arm
               I=indmax(Dist-N/t)
            end
	 end
       end
       # draw the arm
       t+=1
       S[I]+=sample(mu[I])
       N[I]+=1
   end
   recommendation=Best
   return (recommendation,N)
end

function TrackAndStop2(mu,delta,rate)
  # Uses a Tracking of the cummulated sum
  condition = true
  K=length(mu)
  N = zeros(1,K)
  S = zeros(1,K)
  # initialization
  for a in 1:K
      N[a]=1
      S[a]=sample(mu[a])
  end
  t=K
  Best=1
  SumWeights=ones(1,K)/K
  while (condition)
       Mu=S./N
       # Empirical best arm
       IndMax=find(Mu.==maximum(Mu))
       Best=IndMax[floor(Int,length(IndMax)*rand())+1]
       I=1
       # Compute the stopping statistic
       NB=N[Best]
       SB=S[Best]
       muB=SB/NB
       MuMid=(SB+S)./(NB+N)
       Index=collect(1:K)
       splice!(Index,Best)
       Score=minimum([NB*d(muB,MuMid[i])+N[i]*d(Mu[i],MuMid[i]) for i in Index])
       if (Score > rate(t,0,delta))
            # stop
            condition=false
          elseif (t >1000000)
            # stop and output (0,0)
            condition=false
            Best=0
            print(N)
            print(S)
            N=zeros(1,K)
          else
            # continue and sample an arm
	    val,Dist=OptimalWeights(Mu,1e-11)
            SumWeights=SumWeights+Dist
	    # choice of the arm
            if (minimum(N) <= max(sqrt(t) - K/2,0))
               # forced exploration
               I=indmin(N)
            else
               I=indmax(SumWeights-N)
            end
       end
       # draw the arm
       t+=1
       S[I]+=sample(mu[I])
       N[I]+=1
   end
   recommendation=Best
   return (recommendation,N)
end

# Chernoff stopping rule coupled with different sampling rules

function ChernoffTarget(mu,delta,rate,Target=ones(1,length(mu))/length(mu))
  # sampling rule : choose arm maximizing (Target - empirical proportion)
  condition = true
  K=length(mu)
  N = zeros(1,K)
  S = zeros(1,K)
  # initialization
  for a in 1:K
      N[a]=1
      S[a]=sample(mu[a])
  end
  t=K
  Best=1
  while (condition)
       Mu=S./N
       Ind=find(Mu.==maximum(Mu))
       # Empirical best arm
       Best=Ind[floor(Int,length(Ind)*rand())+1]
       # Compute the stopping statistic
       NB=N[Best]
       SB=S[Best]
       muB=SB/NB
       MuMid=(SB+S)./(NB+N)
       Index=collect(1:K)
       splice!(Index,Best)
       Score=minimum([NB*d(muB,MuMid[i])+N[i]*d(Mu[i],MuMid[i]) for i in Index])
       if (Score > rate(t,0,delta))
            # stop
            condition=false
       elseif (t >1000000)
            # stop and return (0,0)
            condition=false
            Best=0
            print(N)
            print(S)
            N=zeros(1,K)
       else
	    I=indmax(Target-N/t)
	    t+=1
	    S[I]+=sample(mu[I])
	    N[I]+=1
	end
   end
   recommendation=Best
   return (recommendation,N)
end

function ChernoffBC(mu,delta,rate)
  # Chernoff stopping rule,  sampling based on the "best challenger"
  # described in experimental section of [Garivier and Kaufmann 2016]
  condition = true
  K=length(mu)
  N = zeros(1,K)
  S = zeros(1,K)
  # initialization
  for a in 1:K
      N[a]=1
      S[a]=sample(mu[a])
  end
  t=K
  Best=1
  while (condition)
       Mu=S./N
       # Empirical best arm
       IndMax=find(Mu.==maximum(Mu))
       Best=IndMax[floor(Int,length(IndMax)*rand())+1]
       I=1
       # Compute the stopping statistic
       NB=N[Best]
       SB=S[Best]
       MuB=SB/NB
       MuMid=(SB+S)./(NB+N)
       Challenger=1
       Score=Inf
       for i=1:K
	  if i!=Best
             score=NB*d(MuB,MuMid[i])+N[i]*d(Mu[i],MuMid[i])
	     if (score<Score)
		 Challenger=i
		 Score=score
	     end
	  end
       end
       if (Score > rate(t,0,delta))
            # stop
            condition=false
          elseif (t >1000000)
            # stop and return (0,0)
            condition=false
            Best=0
            print(N)
            print(S)
            N=zeros(1,K)
          else
            # continue and sample an arm
	    val,Dist=OptimalWeights(Mu,1e-11)
            if (minimum(N) <= max(sqrt(t) - K/2,0))
               # forced exploration
               I=indmin(N)
             else
               # choose between the arm and its Challenger
               I=(NB/(NB+N[Challenger]) < Dist[Best]/(Dist[Best]+Dist[Challenger]))?Best:Challenger
               #I=(d(MuB,MuMid[Challenger])>d(Mu[Challenger],MuMid[Challenger]))?Best:Challenger
               #I=(N[Best]<N[Challenger])?Best:Challenger
            end
       end
       # draw the arm
       t+=1
       S[I]+=sample(mu[I])
       N[I]+=1
   end
   recommendation=Best
   return (recommendation,N)
end

function ChernoffBC2(mu,delta,rate)
  # Chernoff stopping rule + alternative choice between the empirical best and its "challenger"
  # Faster, requires no computation of Optimal Weights
  condition = true
  K=length(mu)
  N = zeros(1,K)
  S = zeros(1,K)
  # initialization
  for a in 1:K
      N[a]=1
      S[a]=sample(mu[a])
  end
  t=K
  Best=1
  while (condition)
       Mu=S./N
       # Empirical best arm
       IndMax=find(Mu.==maximum(Mu))
       Best=IndMax[floor(Int,length(IndMax)*rand())+1]
       I=1
       # Compute the stopping statistic
       NB=N[Best]
       SB=S[Best]
       MuB=SB/NB
       MuMid=(SB+S)./(NB+N)
       Challenger=1
       Score=Inf
       for i=1:K
	  if i!=Best
             score=NB*d(MuB,MuMid[i])+N[i]*d(Mu[i],MuMid[i])
	     if (score<Score)
		 Challenger=i
		 Score=score
	     end
	  end
       end
       if (Score > rate(t,0,delta))
            # stop
            condition=false
          elseif (t >1000000)
            # stop and return (0,0)
            condition=false
            Best=0
            print(N)
            print(S)
            N=zeros(1,K)
          else
            # continue and sample an arm
	    if (minimum(N) <= max(sqrt(t) - K/2,0))
               # forced exploration
               I=indmin(N)
             else
               # choose between the arm and its Challenger
               I=(N[Best]<N[Challenger])?Best:Challenger
               #I=(d(MuB,MuMid[Challenger])>d(Mu[Challenger],MuMid[Challenger]))?Best:Challenger
            end
       end
       # draw the arm
       t+=1
       S[I]+=sample(mu[I])
       N[I]+=1
   end
   recommendation=Best
   return (recommendation,N)
end

function ChernoffKLLUCB(mu,delta,rate)
  # Chernoff stopping rule, KL-LUCB sampling rule
  condition = true
  K=length(mu)
  N = zeros(1,K)
  S = zeros(1,K)
  # initialization
  for a in 1:K
      N[a]=1
      S[a]=sample(mu[a])
  end
  t=K
  Best=1
  while (condition)
       Mu=S./N
       Ind=find(Mu.==maximum(Mu))
       # Empirical best arm
       Best=Ind[round(Int,floor(length(Ind)*rand())+1)]
       # Compute the stopping statistic
       NB=N[Best]
       SB=S[Best]
       muB=SB/NB
       MuMid=(SB+S)./(NB+N)
       Index=collect(1:K)
       splice!(Index,Best)
       Score=minimum([NB*d(muB,MuMid[i])+N[i]*d(Mu[i],MuMid[i]) for i in Index])
       # Find the challenger
       UCB=zeros(1,K)
       LCB=dlow(Mu[Best],rate(t,0,delta)/N[Best])
       for a in 1:K
	  if a!=Best
	     UCB[a]=dup(Mu[a],rate(t,0,delta)/N[a])
          end
       end
       Ind=find(UCB.==maximum(UCB))
       Challenger=Ind[round(Int,floor(length(Ind)*rand())+1)]
       # draw both arms
       t=t+2
       S[Best]+=sample(mu[Best])
       N[Best]+=1
       S[Challenger]+=sample(mu[Challenger])
       N[Challenger]+=1
       # check stopping condition
       condition=(Score <= rate(t,0,delta))
       if (t>1000000)
	  condition=false
          Best=0
          N=zeros(1,K)
       end
   end
   recommendation=Best
   return (recommendation,N)
end

# KL-LUCB [Kaufmann and Kalyanakrishnan 2013]

function KLLUCB(mu,delta,rate)
  condition = true
  K=length(mu)
  N = zeros(1,K)
  S = zeros(1,K)
  # initialization
  for a in 1:K
      N[a]=1
      S[a]=sample(mu[a])
  end
  t=K
  Best=1
  while (condition)
       Mu=S./N
       Ind=find(Mu.==maximum(Mu))
       # Empirical best arm
       Best=Ind[round(Int,floor(length(Ind)*rand())+1)]
       # Find the challenger
       UCB=zeros(1,K)
       LCB=dlow(Mu[Best],rate(t,N[Best],delta)/N[Best])
       for a in 1:K
	  if a!=Best
	     UCB[a]=dup(Mu[a],rate(t,N[a],delta)/N[a])
          end
       end
       Ind=find(UCB.==maximum(UCB))
       Challenger=Ind[round(Int,floor(length(Ind)*rand())+1)]
       # draw both arms
       t=t+2
       S[Best]+=sample(mu[Best])
       N[Best]+=1
       S[Challenger]+=sample(mu[Challenger])
       N[Challenger]+=1
       # check stopping condition
       condition=(LCB < UCB[Challenger])
       if (t>1000000)
	  condition=false
          Best=0
          N=zeros(1,K)
       end
   end
   recommendation=Best
   return (recommendation,N)
end

# Racing algorithms

function ChernoffRacing(mu,delta,rate)
  condition = true
  K=length(mu)
  N = zeros(1,K)
  S = zeros(1,K)
  # initialization
  for a in 1:K
      N[a]=1
      S[a]=sample(mu[a])
  end
  round=1
  t=K
  Best=1
  Remaining=collect(1:K)
  while (length(Remaining)>1)
       # Drawn all remaining arms
       for a in Remaining
	  S[a]+=sample(mu[a])
	  N[a]+=1
       end
       round+=1
       t+=length(Remaining)
       # Check whether the worst should be removed
       Mu=S./N
       MuR=Mu[Remaining]
       MuBest=maximum(MuR)
       IndBest=find(MuR.==MuBest)[1]
       IndBest=IndBest[floor(Int,rand()*length(IndBest))+1]
       Best=Remaining[IndBest]
       MuWorst=minimum(MuR)
       IndWorst=find(MuR.==MuWorst)[1]
       IndWorst=IndWorst[floor(Int,rand()*length(IndWorst))+1]
       if (round*(d(MuBest, (MuBest+MuWorst)/2)+d(MuWorst,(MuBest+MuWorst)/2)) > rate(t,0,delta))
          # remove Worst arm
          splice!(Remaining,IndWorst)
       end
       if (t>1000000)
	  Remaining=[]
          Best=0
          N=zeros(1,K)
       end
   end
   recommendation=Best
   return (recommendation,N)
end

function KLRacing(mu,delta,rate)
  condition = true
  K=length(mu)
  N = zeros(1,K)
  S = zeros(1,K)
  # initialization
  for a in 1:K
      N[a]=1
      S[a]=sample(mu[a])
  end
  round=1
  t=K
  Best=1
  Remaining=collect(1:K)
  while (length(Remaining)>1)
       # Drawn all remaining arms
       for a in Remaining
	  S[a]+=sample(mu[a])
	  N[a]+=1
       end
       round+=1
       t+=length(Remaining)
       # Check whether the worst should be removed
       Mu=S./N
       MuR=Mu[Remaining]
       MuBest=maximum(MuR)
       IndBest=find(MuR.==MuBest)[1]
       Best=IndBest[floor(Int,rand()*length(IndBest))+1]
       Best=Remaining[Best]
       MuWorst=minimum(MuR)
       IndWorst=find(MuR.==MuWorst)[1]
       IndWorst=IndWorst[floor(Int,rand()*length(IndWorst))+1]
       if (dlow(MuBest,rate(t,round,delta)/round) > dup(MuWorst,rate(t,round,delta)/round))
          # remove Worst arm
          splice!(Remaining,IndWorst)
       end
       if (t>1000000)
	  Remaining=[]
          Best=0
          N=zeros(1,K)
       end
   end
   recommendation=Best
   return (recommendation,N)
end

# UGapE [Gabillon et al., 2012]

function UGapE(mu,delta,rate)
  condition = true
  K=length(mu)
  N = zeros(1,K)
  S = zeros(1,K)
  # initialization
  for a in 1:K
      N[a]=1
      S[a]=sample(mu[a])
  end
  t=K
  Best=1
  while (condition)
       Mu=S./N
       Ind=find(Mu.==maximum(Mu))
       # Empirical best arm
       Best=Ind[round(Int,floor(length(Ind)*rand())+1)]
       # Find the challenger
       UCB=zeros(1,K)
       LCB=zeros(1,K)
       for a in 1:K
          UCB[a]=dup(Mu[a],rate(t,N[a],delta)/N[a])
          LCB[a]=dlow(Mu[a],rate(t,N[a],delta)/N[a])
       end
       B=zeros(1,K)
       for a in 1:K
          Index=collect(1:K)
          splice!(Index,a)
          B[a] = maximum(UCB[Index])-LCB[a]
       end
       Value=minimum(B)
       Best=indmin(B)
       UCB[Best]=0
       Challenger=indmax(UCB)
       # choose which arm to draw
       t=t+1
       I=(N[Best]<N[Challenger])?Best:Challenger
       S[I]+=sample(mu[I])
       N[I]+=1
       # check stopping condition
       condition=(Value > 0)
       if (t>1000000)
	  condition=false
          Best=0
          N=zeros(1,K)
       end
   end
   recommendation=Best
   return (recommendation,N)
end

# Pure-Exploration Thompson Sampling [Russo, 2016] + Chernoff Stopping rule

function ChernoffPTS(mu,delta,rate,frac,alpha=1,beta=1)
  # Chernoff stopping rule combined with the PTS sampling rule
  condition = true
  K=length(mu)
  N = zeros(1,K)
  S = zeros(1,K)
  # initialization
  for a in 1:K
      N[a]=1
      S[a]=sample(mu[a])
  end
  t=K
  Best=1
  while (condition)
       Mu=S./N
       Ind=find(Mu.==maximum(Mu))
       # Empirical best arm
       Best=Ind[floor(Int,length(Ind)*rand())+1]
       # Compute the stopping statistic
       NB=N[Best]
       SB=S[Best]
       muB=SB/NB
       MuMid=(SB+S)./(NB+N)
       Index=collect(1:K)
       splice!(Index,Best)
       Score=minimum([NB*d(muB,MuMid[i])+N[i]*d(Mu[i],MuMid[i]) for i in Index])
       if (Score > rate(t,0,delta))
            # stop
            condition=false
       elseif (t >1000000)
            condition=false
            Best=0
            print(N)
            print(S)
            N=zeros(1,K)
       else
            TS=zeros(K)
            for a=1:K
	       TS[a]=rand(Beta(alpha+S[a], beta+N[a]-S[a]), 1)[1]
            end
            I = indmax(TS)
            if (rand()>frac)
               J=I
               condition=true
               while (I==J)
                  TS=zeros(K)
                  for a=1:K
	             TS[a]=rand(Beta(alpha+S[a], beta+N[a]-S[a]), 1)[1]
                  end
                  J= indmax(TS)
               end
               I=J
            end
            # draw arm I
	    t+=1
	    S[I]+=sample(mu[I])
	    N[I]+=1
	end
   end
   recommendation=Best
   return (recommendation,N)
end

max_stars_repo_name: jsfunc/best-arm-identification

max_issues_repo_name: ClaireVern/best-arm-identification

max_forks_repo_name: ClaireVern/best-arm-identification


text:
module Day8

export solvePart1, solvePart2

import DelimitedFiles: readdlm
import DataStructures: deque, popfirst!

function getInput()
    open("day8input.txt") do f
        return readdlm(f, Int64)[:]
    end
end

struct Node
    children
    metadata
end

function extractTreeRecursive(q)
    numChildren = popfirst!(q)
    numMetadata = popfirst!(q)
    children = []
    for i in 1:numChildren
        push!(children, extractTreeRecursive(q))
    end
    metadata = []
    for i in 1:numMetadata
        push!(metadata, popfirst!(q))
    end
    return Node(children, metadata)
end

function extractTree(numbers)
    q = deque(Int64)
    push!(q, numbers...)
    tree = extractTreeRecursive(q)
    return tree
end

function sumMetadata(tree::Node)
    summ = 0
    for child in tree.children
        summ += sumMetadata(child)
    end
    summ += sum(tree.metadata)
    return summ
end

function solvePart1()
    numbers = getInput()
    tree = extractTree(numbers)
    return sumMetadata(tree)
end

function getValue(node::Node)
    n = length(node.children)
    if n == 0
        return sum(node.metadata)
    end
    value = 0
    for entry in node.metadata
        if entry <= n
            value += getValue(node.children[entry])
        end
    end
    return value
end

function solvePart2()
    numbers = getInput()
    tree = extractTree(numbers)
    return getValue(tree)
end

end


text:
#!/usr/bin/env julia
using JuMP

using Cbc

type clause
	x::Int
	y::Int
end

function menu(clauses::Array{clause}, x_max::Int)
	println("Which algorithm would you like to run?")
	println("1. Fast")
	println("2. Exact")
	print("3. Quit\n> ")
	choice = strip(readline(STDIN))
	println()
	if choice == "3"
		println("now exiting max2sat solver ...")
		return true
 	elseif choice == "1"
		best, ans =approx(clauses, x_max)
		println("Clauses satisfied = ", int(best), " out of ", length(clauses))
		println("Truth assignment = ", bitpack(int(getValue(ans).innerArray)))
	elseif choice == "2"
		best, ans = exact(clauses, x_max)
		println("Clauses satisfied = ", int(best), " out of ", length(clauses))
		println("Truth assignment = ", ans)
	else
		println("Bad input, please try again (1-3)")
	end
	return false
end

function prep(filename::String)
	println("Reading file from ", filename)
	f = open(filename)
	lines=readlines(f)
	len = length(lines)
	close(f)

	x_max= 0
	clauses = Array(clause, len)

	for index = 1:len
		items=int(split(lines[index]))
		clauses[index]=clause(items[1],items[2])
		if(abs(items[1]) > x_max)
			x_max = abs(items[1])
		end
		if(abs(items[2]) > x_max)
			x_max = abs(items[2])
		end
	end
	return x_max, clauses
end

function sat_solver(filename::String)
	x_max, clauses = prep(filename)

	while !menu(clauses, x_max) == true
		println()
	end

end

function exact(clauses::Array{clause}, x_max::Int)
	best = 0
	x = BitArray(x_max + 1)
	ans = deepcopy(x)
	while !x[x_max+1]
		temp = check(x, clauses, best)
		if best < temp
			best = temp
			ans = deepcopy(x)
		end
		increment(x)
	end
	return best, ans[1:length(ans)-1]
end

function check_clause(c::clause, bs::BitArray)
	idx = abs(c.x)
	idy = abs(c.y)
	if(c.x < 0)
		if(c.y<0)
			!bs[idx] || !bs[idy]
		else
			!bs[idx] || bs[idy]
		end
	else
		if(c.y<0)
			bs[idx] || !bs[idy]
		else
			bs[idx] || bs[idy]
		end
	end
end

function check(bs::BitArray, clauses::Array{clause}, best::Int)
	count = 0
	len = length(clauses)
	for i = 1: len
		if(check_clause(clauses[i], bs))
			count +=1
		end
		if count + len -i < best
			break
		end
	end
	return count
end

function increment(bs::BitArray)
	bslen = length(bs)
	assert(bslen != 0)
	for loop = 1:bslen
		if ((bs[loop] $= 0x1)== 0x1)
			break
		end
	end
end

function approx(clauses::Array{clause}, x_max::Int)
	m = Model(solver=CbcSolver(threads=6))

	len = length(clauses)
	@defVar(m, 0<= x[1:x_max] <= 1, Int)
	#@defVar(m, 0<= x[1:x_max] <= 1)
	@defVar(m, 0<= z[1:len] <= 1)

	for i = 1:len
		idx = abs(clauses[i].x)
		idy = abs(clauses[i].y)
		if(clauses[i].x < 0)
			if(clauses[i].y<0)
				@addConstraint(m, (1-x[idx]) + (1-x[idy]) >= z[i])
			else
				@addConstraint(m, (1-x[idx]) + x[idy] >= z[i])
			end
		else
			if(clauses[i].y<0)
				@addConstraint(m, x[idx] + (1-x[idy]) >= z[i])
			else
				@addConstraint(m, x[idx] + x[idy] >= z[i])
			end
		end
	end
	@setObjective(m, :Max, sum{z[i], i=1:len} )

	status = solve(m)
	#println(m)

	satisfied = getObjectiveValue(m)

	println(status)

	return satisfied, x
end

function main()
	sat_solver("instance.txt")

end

main()

max_stars_repo_path: Julia Implemenation/Max2Sat.jl

max_issues_repo_path: Julia Implemenation/Max2Sat.jl

max_forks_repo_path: Julia Implemenation/Max2Sat.jl


text:
using Combinatorics

function part1(input)
    ingredients = reduce(hcat, parse_line.(eachline(input)))
    return foldl(max, (score1(ingredients[1:(end - 1), :], proportions)
                       for proportions1 in partitions(100, size(ingredients, 2))
                       for proportions in permutations(proportions1)))
end

function part2(input)
    ingredients = reduce(hcat, parse_line.(eachline(input)))
    return foldl(max, (score2(ingredients, proportions)
                       for proportions1 in partitions(100, size(ingredients, 2))
                       for proportions in permutations(proportions1)))
end

parse_line(line) = parse.(Int, rstrip.(split(line)[3:2:11], ','))

score1(ingredients, proportions) = prod(max.(0, ingredients * proportions))

function score2(ingredients, proportions)
    n = max.(0, ingredients * proportions)
    return prod(n[1:(end - 1)]) * (n[end] == 500)
end


text:
#############################
##  Numerical experiments associated with section 5.4.1
#############################
## Nonlinear equations are treated by first defining an objective function, then using standard iterative methods
## such as Newton's method via NLsolve.jl. Using more sophisticated methods, linesearches etc. can sometimes result in significant speed boosts.
##

using NLsolve, ApproxFun, MultivariateOrthogonalPolynomials, BandedMatrices, BlockBandedMatrices, Plots
using SparseVolterraExamples

#####################################################
## Problem in Equation (25), u_1(x) = exp(x)
#####################################################
####
## Define g(x), the nonlinearity and the kernel
g(x) = exp(x)+1/3*x*(1-exp(3*x))
    gFun = Fun(g,Jacobi(0,1,0..1))
    f(u) = u^3
    K(x,y) = x
    n = 15 # Desired polynomial order of approximation for accuracy
####
## Define the objective function and compute the solution using a standard Newton method
objective(x) = triNonLinearVolterraObjective(x,f,K,gFun,n)
sol = nlsolve(objective,zeros(n),method=:newton,ftol=1e-8)
####
## Plot function aganinst analytic solution
plot(Fun(Jacobi(0,1,0..1),sol.zero),grid=false,legend=:topleft,xlabel="x",ylabel="u(x)",label="sparse method")
plot!(exp,0,1,grid=false,legend=:topleft,xlabel="x",ylabel="u(x)",label="analytic exp(x)")
####
## Or plot error. A simple loop over polynomial degree can then produce Figure 7(a).
plot(Fun(Jacobi(0,1,0..1),sol.zero)-Fun(x->exp(x),Jacobi(0,1,0..1)),grid=false,legend=:bottomleft,xlabel="x",ylabel="error",label=false)

#####################################################
## Problem in Equation (26), u_2(x) = sin(x)
#####################################################
####
## Define g(x), the nonlinearity and the kernel
g(x) = sin(x)+1/4*sin(x)^2-1/4*x^2
    gFun = Fun(g,Jacobi(0,1,0..1))
    f(u) = u^2
    K(x,y) = x-y
    n = 14 # Desired polynomial order of approximation for accuracy
####
## Define the objective function and compute the solution using a standard Newton method
objective(x) = triNonLinearVolterraObjective(x,f,K,gFun,n)
sol = nlsolve(objective,zeros(n),method=:newton,ftol=1e-15)
####
## Plot function aganinst analytic solution
plot(Fun(Jacobi(0,1,0..1),sol.zero),grid=false,legend=:topleft,xlabel="x",ylabel="u(x)",label="sparse method")
plot!(sin,0,1,grid=false,legend=:topleft,xlabel="x",ylabel="u(x)",label="analytic sin(x)")
####
## Or plot error. A simple loop over polynomial degree can then produce Figure 7(b).
plot(Fun(Jacobi(0,1,0..1),sol.zero)-Fun(x->sin(x),Jacobi(0,1,0..1)),grid=false,legend=:bottomleft,xlabel="x",ylabel="error",label=false)

max_stars_repo_path: examples/Section 5-4-1 - Nonlinear Volterra integral equations.jl


text:
import JSON
import Requests
using Requests.get
using Requests.Response

using Mongo
"""
global definition of dictionary to store metadata per database
"""
metaDictionary = Dict{UTF8String,Dict{UTF8String,Any}}()
"""
global definition of mongodb client
"""
client = MongoClient()
securityCollection = MongoCollection(client,"aimsqaunt","security")

dataCollection = MongoCollection(client , "aimsquant","data")
"""
function to return the GET arguements for downloading the meta data for
the datasets in the particular database defined by database_code
"""
function getqueryargs(database_code:: AbstractString; per_page::Int = 100 ,
 sort_by="id" , page::Int = 1 )

    queryArgs = Dict{Any , Any}("database_code" => database_code,
                                "per_page" => per_page,
                                "sort_by" => sort_by,
                                "page" =>page,
                                "api_key" =>getapikey())

    return queryArgs
end

"""
function to insert metadata per database
"""
function insertmetadata(database_code::AbstractString , data :: Dict{UTF8String,Any})
    metaDictionary[database_code] = data

    #To check if insertion is successful
    temp = get(metaDictionary , database_code , 0)

    if temp == 0
        error("Error in insertion")
        return false
    else
        return true

    end
end

"""
function to get the metadata from the metadata database
"""
function getmetadata(database_code::AbstractString)

   temp = get(metaDictionary , database_code , 0)

    if temp == 0
        error("Value doesn't exist in database")
        return false
    else
        return temp

    end
end

"""
function to get the Quandl API key
"""
function getapikey()
    if !ispath(joinpath(pwd(),"token/"))
        error("Api Key is not initialized")
    end

    api_key = readall(joinpath(pwd(),"token/auth_token"))

    if api_key == ""
        println("Empty API Key")
    else
        println("Using API key " , api_key)
    end

    return api_key

end

"""
function to get the basic url for the Quandl
"""
function getbaseurl()
    path = "https://www.quandl.com/api/"
    return path
end

"""
function to download the metadata for all the datasets present in the database ("NSE" , "WIKI")
"""
function getlistofdatasets(code::AbstractString)

    if length(code) == 0
        error("Please pass a valid code!")
    end
    final_path = getbaseurl() * "v3/datasets.json/"

    resp = get(final_path , query = getqueryargs(code))

    if resp.status != 200
        println("Error in processing the request")
    end

    data = Requests.json(resp)

    insertmetadata(code , data["meta"])

    #println(data["meta"])

    #println ("Reached here")
    storealldatasetscode(code)
    return nothing
end
"""
function to check if data exists in mongodb document
"""
function checkifexistsinmongodbdocument(collection:: MongoCollection , data :: Dict{UTF8String,Any} )
    dataCount = count(collection , data)

    if dataCount == 0
        return false
    else
        return true
    end
end

"""
function to insert security data
The data is not available in Quandl Dataset for the started trading day and ended trading day , So if another data source is added then its needed to modified
"""
function insertintosecuritydocument(data :: Dict{UTF8String,Any})

    securityID = get(data , "id" , "NULL")
    ISIN = get(data , "ISIN" , "NULL")
    ticker = get(data , "dataset_code" , "NULL")
    exchange = get(data , "database_code" , "NULL")
    frequency = get(data , "frequency" , "NULL")
    dataType = get(data , "type" , "NULL")
    sourceName  = "Quandl"
    QuandlDict = Dict{Any,Any}("name" => sourceName,
                               "id" => get(data , "id" , "NULL"),
                               "newest_available_date" => get(data , "newest_available_date" , "NULL"),
                               "oldest_available_date" => get(data , "oldest_available_date" , "NULL"),
                               "refreshed_at" => get(data ,"refreshed_at"  , "NULL"),
                               "description" => get(data , "description" , "NULL"),
                               "dataset_code" => get(data , "dataset_code" , "NULL"),
                               "database_code" => get(data , "database_code" , "NULL"))

    dataDict = Dict{Any , Any}("dataType" => dataType,
                                "frequency" => frequency,
                                "source" => QuandlDict
                              )
    securityData = Dict{Any , Any}("securityID" => securityID,
        "ISIN" => ISIN,
        "ticker" => ticker,
        "exchange" => exchange,
        "name" => sourceName,
        "dataSources" => dataDict)
    insert(securityCollection , securityData)

end

"""
function to extract all the datasets code from metadata
"""
function storealldatasetscode(code::AbstractString)

    if length(code) == 0
        error("Please pass a valid code!")
    end
    metaDict = getmetadata(code)
    final_path = getbaseurl() * "v3/datasets.json/"

    if metaDict == false
        error(code * "data doesn't exist in database !!!")
    else
        total_pages = metaDict["total_pages"]
        for i = 1:total_pages
            resp = get(final_path , query = getqueryargs(code , page = i))

            if resp.status != 200
                error("Error in processing the query")
            else
                respJSON = Requests.json(resp)
                dataArray = respJSON["datasets"]
                len = length(dataArray)
                for j = 1:len
                    dataset  = dataArray[j]
                    println(dataset)
                    insertintosecuritydocument(dataset)
                    database_code = dataset["database_code"]
                    dataset_code = dataset["dataset_code"]
                    getDataURL = getbaseurl() *  "v3/datasets/" *database_code * "/" * dataset_code *".json"
                    queryArgs = Dict{Any , Any}("api_key" => getapikey())
                    dataResp = get(getDataURL , query = queryArgs)
                    #println(dataResp)
                    #println ("Reached here in loop")
                    if dataResp.status != 200
                        error("Error in processing the query dataset_code query")
                    else
                        dataRespJSON = Requests.json(dataResp)
                        handlesecuritydata(dataRespJSON["dataset"])

                    end
                end
            end
        end
    end
end

"""
function to insert in data document
"""
function insertindatadocument(securityID::AbstractString,column, data :: Dict{Int , Array{Any , 1}})

    for k in keys(data)
        data = get(data , k, "NULL")
        if data == "NULL"
            println("Something Wrong While inserting Data")
        else
            year  = k
            quandlData = Dict{Any ,Any}("priority" => 1,
                                        "source" => Dict{Any,Any}("name" => "Quandl","id" => securityID),
                                        "column" =>column,
                                        "data" => data
                                        )
            toInsertDict = Dict{Any ,Any}("securityID" => securityID , "Year" => year , "data" => quandlData)
            insert(dataCollection , toInsertDict)
        end

    end
end

"""
function to handle Data for each security
"""
function handlesecuritydata(data :: Dict{UTF8String,Any})
    securityID = get(data , "id" , "NULL")
    column = get(data ,"column", "NULL")
    source = Dict{Any ,Any}("name" => "Quandl" ,
                            "id" => securityID)
    dataToInsert = filterdata(data["data"])

    #insertindatadocument(securityID,column,dataToInsert)

end

"""
function to filter data base on year
"""
function filterdata(data :: Array{Any , 1})
    length  = length(data)
    dataPerYear = Dict{Int, Array{Any , 1}}
    previousDate = 0
    arr = Array{Any , 1}[]
    for i = 1 : length
        dataRow = data[i]
        date = parse(Int , (split(dataRow[1] , "-")[1]) )
        if date == previousDate
           push!(arr , dataRow)
        else
            if previousDate == 0
                nothing
            else
                insert!(dataPerYear ,previousDate , arr)
            end
            arr = Array{Any , 1}[]
            previousDate = date
            push!(arr , dataRow)
        end
    end
    return dataPerYear
end
"""
function to set the Quandl API key
"""

function setauthtoken(token::AbstractString)

    if length(token) != 20 && length(token) != 0
        error("Invalid Token : must be 20 characters long or be an empty")
    end
    println(pwd())

    a = joinpath(pwd(),"token")
    println("Printing the path")
    println(a)
    if !ispath(joinpath(pwd(),"token/"))
        println("Creating new directory")
        mkdir(joinpath(pwd(),"token"))
    end

    open(joinpath(pwd(),"token/auth_token"),"w") do token_file
        write(token_file , token)
    end

    return nothing
end

#println("Starting the quandl")
#setauthtoken("JHKaDwdS-RtM26RxPauV")
getlistofdatasets("NSE")

max_stars_repo_name: AnkitAggarwalPEC/Quandl-and-Julia

max_issues_repo_name: AnkitAggarwalPEC/Quandl-and-Julia

max_forks_repo_name: AnkitAggarwalPEC/Quandl-and-Julia


text:
#! format: off

using SIIPExamples #hide
using PowerSimulationsDynamics
using PowerSystems
using Sundials
using Plots
gr()
PSD = PowerSimulationsDynamics

file_dir = joinpath(
    pkgdir(SIIPExamples),
    "script",
    "4_PowerSimulationsDynamics_examples",
    "Data",
)
omib_sys = System(joinpath(file_dir, "omib_sys.json"))

time_span = (0.0, 30.0)
perturbation_trip = BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1")
sim = PSD.Simulation(ResidualModel, omib_sys, pwd(), time_span, perturbation_trip)

show_states_initial_value(sim)

x0_init = PSD.get_initial_conditions(sim)

PSD.execute!(
    sim, #simulation structure
    IDA(), #Sundials DAE Solver
    dtmax = 0.02,
); #Arguments: Maximum timestep allowed

results = read_results(sim)

angle = get_state_series(results, ("generator-102-1", :δ));
plot(angle, xlabel = "time", ylabel = "rotor angle [rad]", label = "rotor angle")

volt = get_voltage_magnitude_series(results, 102);
plot(volt, xlabel = "time", ylabel = "Voltage [pu]", label = "V_2")

sim2 = PSD.Simulation(ResidualModel, omib_sys, pwd(), time_span)

small_sig = small_signal_analysis(sim2)

small_sig.reduced_jacobian

small_sig.eigenvalues

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: NREL-SIIP/Examples

max_issues_repo_name: NREL-SIIP/Examples

max_forks_repo_name: NREL-SIIP/Examples


text:
@time using DamagedShearBand
const DSB = DamagedShearBand
#using Plots # add Plots before using it
#pyplot(size=(1000,400)) # slow but interactive backend for Plots, add PlotlyJS package first

# imposed values
ϵ̇₁₁ = -1e-6
σ₂₂ = -1e6 # negative in compression

# initial higher compressive stress multiplier (S = σ₁₁/σ₂₂)
S = 1.0

# Create a Rheology type instance containing elastic moduli and damage parameters. Change default values by  supplying keywords arguments.
r = DSB.Rheology(D₀=0.4,n=34.0)

# initial damage
D_i = r.D₀

# time integration parameters :
ϵ₁₁_goal = -0.002
tspan = (0.0,ϵ₁₁_goal/ϵ̇₁₁) # initial and final simulation time
Δt_i = 0.0001 # initial timestep

# Build stress and strain tensors
σᵢⱼ_i = DSB.build_principal_stress_tensor(r,S,σ₂₂,D_i ; abstol=1e-15) # takes care of the plane strain constraint by solving non linear out of plane strain wrt σ₃₃ using Newton algorithm
ϵᵢⱼ_i = DSB.compute_ϵij(r,D_i,σᵢⱼ_i)

n_iter_saved = 5000
n_iter_printed = 10
op = DSB.OutputParams(save_frequency = n_iter_saved/(tspan[2]-tspan[1]),
                      save_period = nothing,
                      print_frequency = n_iter_printed/(tspan[2]-tspan[1]),
                      print_period = nothing)

sp = DSB.SolverParams(newton_abstol = 1e-12,
                      newton_maxiter = 100,
                      time_maxiter = nothing,
                      e₀ = (D=1e-1, σ=100.0, ϵ=1e-3))

p = DSB.Params(op,sp)
# Integrate over time, by solving for the unknowns [σ₁₁next, σ₃₃next, ϵ₂₂next] at each timestep
# -- e₀ is the absolute tolerance used to adapt time stepping, may be a scalar or a NamedTuple with keys (D,σ,ϵ)
# -- abstol is the absolute tolerance on the unknowns in the Newton solver
t_vec, σᵢⱼ_vec, ϵᵢⱼ_vec, D_vec = DSB.adaptative_time_integration(r,p,σᵢⱼ_i,ϵᵢⱼ_i,D_i,ϵ̇₁₁,Δt_i,tspan)

# extract σ₁₁ and ϵ₁₁
σ₁₁_vec = [σᵢⱼ[1,1] for σᵢⱼ in σᵢⱼ_vec]
ϵ₁₁_vec = [ϵᵢⱼ[1,1] for ϵᵢⱼ in ϵᵢⱼ_vec]
ϵkk_vec = [tr(ϵᵢⱼ) for ϵᵢⱼ in ϵᵢⱼ_vec]
KI_vec = [DSB.compute_KI(r,σᵢⱼ_vec[i], D_vec[i]) for i in eachindex(σᵢⱼ_vec)]
free_energy_vec = [DSB.compute_free_energy(r,D_vec[i],ϵᵢⱼ_vec[i]) for i in eachindex(t_vec)]
## plot time series
pp1 = plot(t_vec,ϵkk_vec)
pp2 = plot(t_vec,σ₁₁_vec)
plot!(pp1,pp2)
plot(t_vec,ϵ₁₁_vec)
plot(t_vec,KI_vec./r.K₁c)
plot!(t_vec,D_vec)
plot(t_vec,strain_energy_vec)

## plot stress-strain along x with options on plot
layout = @layout [a ; b]

p1 = plot(-ϵ₁₁_vec,-σ₁₁_vec)
#plot(p1,-ϵ₁₁_vec,-σ₁₁_vec)
title!(p1,"D0 = $(r.D₀) & epsdot = $(ϵ̇₁₁) & p=$(-σ₂₂)")
ylabel!("-σ₁₁")

p2 = plot(-ϵ₁₁_vec,D_vec,c=:red)
#plot(p2,-ϵ₁₁_vec,D_vec,c=:red)
ylabel!("D")
xlabel!("-ϵ₁₁")

pl = plot(p1,p2,layout=layout)
#xlims!((0,0.005))

#savefig(p1,"./examples/singlepoint_D0=$(r.D₀)_var_epsdot.pdf")


text:
using Test, Dates, Temporal

@testset "Input/Output" begin
    @testset "Flat Files" begin
        rows, cols = 5, 3
        filepath = "tmp.csv"
        a = TS(randn(rows, cols))
        tswrite(a, filepath)
        @test filepath in readdir()
        b = tsread(filepath)
        @test a == b
    end
    @testset "Web Downloads" begin
        @testset "Quandl" begin
            crude = TS()
            try
                database = "CHRIS"
                table = "CME_CL1"
                crude = quandl("$database/$table", from="2016-01-01", thru=string(today()))
                @test size(crude, 2) == 8
                @test op(crude).fields == [:Open]
                meta = quandl_meta("CHRIS", "CME_CL1")
                @test meta["database_code"] == database
                perpage = 10
                search = quandl_search(db="CHRIS", qry="corn", perpage=perpage)
                @test length(search["datasets"]) == perpage
            catch
                @test_skip !isempty(crude)
            end
        end
        @testset "Yahoo" begin
            apple = TS()
            try
                apple = yahoo("AAPL", from="2016-01-01")
                @test size(apple, 2) == 6
                @test op(apple).fields == [:Open]
            catch
                @test_skip !isempty(apple)
            end
        end
        @testset "Google" begin
            apple = TS()
            try
                apple = google("AAPL", from="2016-01-01")
                @test op(apple).fields == [:Open]
            catch
                @test_skip !isempty(apple)
            end
        end
    end
end


text:
export fwi_objective

fwi_objective(model_full::Model, source::judiVector, dObs::judiVector) = fwi_objective(model_full, source, dObs; options=Options())

function fwi_objective(model_full::Model, source::judiVector, dObs::judiVector, options::Options)
# Setup time-domain linear or nonlinear foward and adjoint modeling and interface to OPESCI/devito

    # assert this is for single source LSRTM
    @assert source.nsrc == 1 "Multiple sources are used in a single-source fwi_objective"
    @assert dObs.nsrc == 1 "Multiple-source data is used in a single-source fwi_objective"

    # Load full geometry for out-of-core geometry containers
    dObs.geometry = Geometry(dObs.geometry)
    source.geometry = Geometry(source.geometry)

    # Limit model to area with sources/receivers
    if options.limit_m == true
        model = deepcopy(model_full)
        model, _ = limit_model_to_receiver_area(source.geometry,dObs.geometry,model,options.buffer_size)
    else
        model = model_full
    end

    # Set up Python model
    modelPy = devito_model(model, options)
    dtComp = get_dt(model; dt=options.dt_comp)

    # Extrapolate input data to computational grid
    qIn = time_resample(source.data[1], source.geometry, dtComp)[1]
    obsd = (typeof(dObs.data[1]) == SegyIO.SeisCon || typeof(dObs.data[1]) == H5SeisCon) ? convert(Array{Float32,2}, dObs.data[1][1].data) : dObs.data[1]
    dObserved = time_resample(obsd, dObs.geometry, dtComp)[1]

    # Set up coordinates
    src_coords = setup_grid(source.geometry, model.n)  # shifts source coordinates by origin
    rec_coords = setup_grid(dObs.geometry, model.n)    # shifts rec coordinates by origin

    if options.optimal_checkpointing == true
        argout1, argout2 = pycall(ac."J_adjoint_checkpointing", Tuple{Float32, Array{Float32, modelPy.dim}},
                                  modelPy, src_coords, qIn,
                                  rec_coords, dObserved, is_residual=false, return_obj=true, isic=options.isic,
                                  t_sub=options.subsampling_factor, space_order=options.space_order)
    elseif ~isempty(options.frequencies)
        argout1, argout2 = pycall(ac."J_adjoint_freq", Tuple{Float32,  Array{Float32, modelPy.dim}},
                                  modelPy, src_coords, qIn,
                                  rec_coords, dObserved, is_residual=false, return_obj=true, isic=options.isic,
                                  freq_list=options.frequencies, t_sub=options.subsampling_factor,
                                  space_order=options.space_order)
    else
        argout1, argout2 = pycall(ac."J_adjoint_standard", Tuple{Float32,  Array{Float32, modelPy.dim}},
                                  modelPy, src_coords, qIn,
                                  rec_coords, dObserved, is_residual=false, return_obj=true,
                                  t_sub=options.subsampling_factor, space_order=options.space_order,
                                  isic=options.isic)
    end
    argout2 = remove_padding(argout2, modelPy.padsizes; true_adjoint=options.sum_padding)
    if options.limit_m==true
        argout2 = extend_gradient(model_full, model, argout2)
    end
    return Ref{Float32}(argout1), PhysicalParameter(argout2, model_full.d, model_full.o)
end


text:
using WeibullCount

begin
    """
    Tests the Van Wijngaarden transformation on the Leibniz formula for π:

        1 - (1 / 3) + (1 / 5) - (1 / 7) + ... = π / 4

    """

    a = 1 ./ (1:2:100)
    coeffs = ones(length(a))
    coeffs[2:2:length(a)] .= -1  # alternating minus-sign

    true_sum = π / 4
    partial = sum(a .* coeffs)
    trans = WeibullCount.van_wijngaarden_transformation(a)

    @test abs(true_sum - partial) .> abs(true_sum - trans)
end