text:
@everywhere using Distributions
include("../utils.jl")

# compute As given Basis B and Parameters θ and store the results in As
@everywhere function computeAs(
    B::Matrix{Float64},
    θ::Array{Float64,1},
    As::Array{Float64,3},
    Avees::Array{Float64,1} = Array(Float64, size(θ,1))
  )

  Avees = B*θ
  n::Int64 = length(Avees)
  nCells::Int64 = round(Int, n/6)
  for i = 1:nCells
    As[i,1,1] = Avees[6*(i-1)+1]
    As[i,1,2] = Avees[6*(i-1)+2]
    As[i,1,3] = Avees[6*(i-1)+3]
    As[i,2,1] = Avees[6*(i-1)+4]
    As[i,2,2] = Avees[6*(i-1)+5]
    As[i,2,3] = Avees[6*(i-1)+6]
  end
end

@everywhere function computeAndEvaluateTransformationSingleThread(
    B::Matrix{Float64},
    θ::Array{Float64,1},
    As::Array{Float64,3},
    X::Array{Float64,2},
    Y::Array{Float64,2},
    Tx::Array{Float64,2},
    t::Float64,
    nCols::Int64,
    nRows::Int64,
    incX::Float64,
    incY::Float64,
    Nsteps::Int64,
    Δt::Float64,
    nSteps::Int64,
    ΔtAs::Array{Float64,3},
    expT::Array{Float64,3} = Array(Float64, size(As,1), 3, 3),
    Avees::Array{Float64,1} = Array(Float64, size(θ,1)),
    bdryConstraints::Bool = false,
  )
  computeAs(B,θ,As,Avees)
  solver2d(X,Tx,t,As,nCols,nRows,incX,incY,Nsteps,Δt,nSteps,ΔtAs,bdryConstraints,expT)

  nPts = size(X,1)
  leastSquares = 0
  for i=1:nPts
    leastSquares += (Tx[i,1]-Y[i,1])^2 + (Tx[i,2]-Y[i,2])^2
  end
  leastSquares
end

@everywhere function computeAndEvaluateTransformationParallel(
    B::Matrix{Float64},
    θ::Array{Float64,1},
    As::Array{Float64,3},
    X::Array{Float64,2},
    Y::Array{Float64,2},
    Tx::DArray{Float64,2},
    t::Float64,
    nCols::Int64,
    nRows::Int64,
    incX::Float64,
    incY::Float64,
    Nsteps::Int64,
    Δt::Float64,
    nSteps::Int64,
    ΔtAs::Array{Float64,3},
    expT::Array{Float64,3} = Array(Float64, size(As,1), 3, 3),
    Avees::Array{Float64,1} = Array(Float64, size(θ,1)),
    bdryConstraints::Bool = false,
  )
  computeAs(B,θ,As,Avees)
  solver2dParallel(X,Tx,t,As,nCols,nRows,incX,incY,Nsteps,Δt,nSteps,ΔtAs,bdryConstraints,expT)

  threads::Array{Int64,1} = workers()
  nThreads::Int64 = length(threads)
  refs::Array{RemoteRef,1} = Array(RemoteRef, nThreads)
  partialLeastSquares::Array{Float64,1} = Array(Float64, nThreads)
  for index::Int64=1:nThreads
    p = threads[index]
    refs[index] = @spawnat(p, (
      indicesOnThread::Array{Int64,1} = localindexes(Tx)[1];
      localArray = localpart(Tx);
      currLeastSquares::Float64 = 0.0;
      for i::Int64 = 1:length(indicesOnThread);
        j = indicesOnThread[i];
        currLeastSquares += (localArray[i,1]-Y[j,1])^2 + (localArray[i,2]-Y[j,2])^2;
      end;
      currLeastSquares
      ))
  end

  # wait for the workers to finish and sum
  map(wait, refs)
  leastSquares::Float64 = 0.0
  for index::Int64=1:nThreads
    leastSquares += fetch(refs[index])
  end
  leastSquares
end

@everywhere function computeAndEvaluateTransformationParallelShared(
    B::Matrix{Float64},
    θ::Array{Float64,1},
    As::Array{Float64,3},
    X::Array{Float64,2},
    Y::Array{Float64,2},
    Tx::SharedArray{Float64,2},
    t::Float64,
    nCols::Int64,
    nRows::Int64,
    incX::Float64,
    incY::Float64,
    Nsteps::Int64,
    Δt::Float64,
    nSteps::Int64,
    ΔtAs::Array{Float64,3},
    expT::Array{Float64,3} = Array(Float64, size(As,1), 3, 3),
    Avees::Array{Float64,1} = Array(Float64, size(θ,1)),
    bdryConstraints::Bool = false,
  )
  computeAs(B,θ,As,Avees)
  solver2dParallelShared(X,Tx,t,As,nCols,nRows,incX,incY,Nsteps,Δt,nSteps,ΔtAs,bdryConstraints,expT)
  nPts = size(X,1)
  leastSquares = 0
  for i=1:nPts
    leastSquares += (Tx[i,1]-Y[i,1])^2 + (Tx[i,2]-Y[i,2])^2
  end
  leastSquares
end

@everywhere function calculateLogProbGaussian(
    θ::Array{Float64,1},
    θmean::Array{Float64,1},
    invCov::Array{Float64,2},
  )
  -0.5*(θ-θmean)'*invCov*(θ-θmean)
end

function metropolis(
    f::Function,
    θ::Array{Float64,1},
    iterations::Int64,
    σ::Float64,
    priorCov::Array{Float64,2},
    jumpCov::Array{Float64,2},
    storeTrace::Bool = false,
  )
  dim::Int64 = size(θ,1)
  invPriorCov::Array{Float64,2} = inv(priorCov)
  priorDist::Distributions.MvNormal = MvNormal(zeros(dim), priorCov)
  jumpDist::Distributions.MvNormal = MvNormal(zeros(dim), jumpCov)
  prevθ::Array{Float64,1} = copy(θ)
  newθ::Array{Float64,1} = Array(Float64, dim)
  prevVal::Float64 = -0.5*f(θ)/(σ^2)
  currVal::Float64 = 0.0
  prevPrior::Float64 = calculateLogProbGaussian(θ, zeros(dim), invPriorCov)[1]  #log(pdf(priorDist,θ))
  currPrior::Float64 = prevPrior
  maxVal::Float64 = prevVal
  bestθ::Array{Float64,1} = copy(θ)
  accepted::Int64 = 0

  if storeTrace
    llTrace::Array{Float64,1} = Array(Float64, iterations)
    acceptanceTrace::Array{Float64,1} = Array(Float64, iterations)
    lpTrace::Array{Float64,1} = Array(Float64, iterations)
    lpllTrace::Array{Float64,1} = Array(Float64, iterations)
  end

  for i=1:iterations
    if i%1000 == 0
      println("Running iteration $i out of $iterations")
    end
    acceptanceRate = accepted/i
    newθ = rand(jumpDist) + prevθ
    currVal = -0.5*f(newθ)/(σ^2)
    currPrior = calculateLogProbGaussian(newθ, zeros(dim), invPriorCov)[1] # log(pdf(priorDist,newθ))
    logRatio = (currVal-prevVal) + (currPrior-prevPrior)
    if logRatio>0 || logRatio < log(rand())
      accepted += 1
      prevθ = newθ
      prevVal = currVal
      prevPrior = currPrior
      if currVal > maxVal
        maxVal = currVal
        bestθ = newθ
      end
    end
    if storeTrace
      llTrace[i] = currVal
      lpTrace[i] = currPrior
      acceptanceTrace[i] = acceptanceRate
      lpllTrace[i] = currVal+currPrior
    end
  end

  if storeTrace
    return bestθ, acceptanceTrace, llTrace, lpTrace, lpllTrace
  end
  bestθ
end

function particleFilter(
    f::Function, # likelihood
    θ::Array{Float64,1},
    nParticles::Int64,
    iterations::Int64,
    σ::Float64,
    priorCov::Array{Float64,2},
    jumpCov::Array{Float64,2},
    plotWeights::Bool = false,
    storeTrace::Bool = false,
  )
  dim::Int64 = size(θ,1)
  invPriorCov::Array{Float64,2} = inv(priorCov)
  priorDist::Distributions.MvNormal = MvNormal(zeros(dim), priorCov)
  jumpDist::Distributions.MvNormal = MvNormal(zeros(dim), jumpCov)
  prevSamples::Array{Float64, 2} = rand(priorDist, nParticles)'
  currSamples::Array{Float64, 2} = Array(Float64, nParticles, dim)
  newθ::Array{Float64, 1} = Array(Float64, dim)
  weights::Array{Float64,1} = Array(Float64, nParticles)
  cumSumWeights::Array{Float64,1} = Array(Float64, nParticles)
  perturbation::Array{Float64,1} = Array(Float64, dim)
  sumWeights::Float64 = 0.0
  maxVal::Float64 = -Inf
  bestθ::Array{Float64,1} = θ

  if storeTrace
    llTrace::Array{Float64,1} = zeros(iterations)
    lpTrace::Array{Float64,1} = zeros(iterations)
    lpllTrace::Array{Float64,1} = zeros(iterations)
  end

  for i = 1:iterations
    if i%10 == 0
      println("Running iteration $i out of $iterations")
    end
    # update weights
    for j = 1:nParticles
      for k = 1:dim
        newθ[k] = prevSamples[j,k] # avoiding squeeze
      end
      logLikelihood::Float64 = -0.5*f(newθ)/(σ^2)
      logPrior::Float64 = calculateLogProbGaussian(newθ, zeros(dim), invPriorCov)[1]
      weights[j] = logLikelihood + logPrior

      if weights[j] > maxVal
        maxVal = weights[j]
        bestθ = newθ
        if storeTrace
          llTrace[i] = logLikelihood
          lpTrace[i] = logPrior
          lpllTrace[i] = logLikelihood + logPrior
        end
      end
    end

    if storeTrace
      if lpllTrace[i] == 0 && i > 1
        llTrace[i] = llTrace[i-1]
        lpTrace[i] = lpTrace[i-1]
        lpllTrace[i] = lpllTrace[i-1]
      end
    end

    maxWeight = weights[1]
    for j = 1:nParticles
      if weights[j] != Inf && weights[j] > maxWeight
        maxWeight = weights[j]
      end
    end

    weights -= maxWeight
    sumWeights = 0.0
    for j = 1:nParticles
      weights[j] = e^weights[j]
      sumWeights += weights[j]
    end

    weights /= sumWeights
    if plotWeights
      PyPlot.clf()
      PyPlot.plot(1:nParticles, cumsum(sort(weights)),label="dist", color="blue")
      savefig(string("plot", i, ".svg"))
    end

    # resample
    cumSumWeights = cumsum(weights)
    for j = 1:nParticles
      cumWeightSample::Float64 = rand()
      for k = 1:nParticles
        if cumSumWeights[k] >= cumWeightSample
          currSamples[j,:] = prevSamples[k,:]
          break
        end
      end
    end

    # perturb
    for j = 1:nParticles
      perturbation = rand(jumpDist)
      for k = 1:dim
        prevSamples[j,k] = perturbation[k] + currSamples[j,k]
      end
    end
  end

  if storeTrace
    return bestθ, llTrace, lpTrace, lpllTrace
  end
  bestθ
end

function particleFilterParallel(
    f::Function, # likelihood
    θ::Array{Float64,1},
    nParticles::Int64,
    iterations::Int64,
    σ::Float64,
    priorCov::Array{Float64,2},
    jumpCov::Array{Float64,2},
    plotWeights::Bool = false,
    storeTrace::Bool = false,
  )
  dim::Int64 = size(θ,1)
  invPriorCov::Array{Float64,2} = inv(priorCov)
  priorDist::Distributions.MvNormal = MvNormal(zeros(dim), priorCov)
  jumpDist::Distributions.MvNormal = MvNormal(zeros(dim), jumpCov)
  prevSamples::SharedArray{Float64, 2} = convert(SharedArray, rand(priorDist, nParticles)')
  currSamples::SharedArray{Float64, 2} = SharedArray(Float64, nParticles, dim)
  weights::SharedArray{Float64,1} = SharedArray(Float64, nParticles)
  cumSumWeights::Array{Float64,1} = Array(Float64, nParticles)
  perturbation::Array{Float64,1} = Array(Float64, dim)
  sumWeights::Float64 = 0.0
  maxVal::SharedArray{Float64,1} = convert(SharedArray, [-Inf])
  bestθ::SharedArray{Float64,1} = convert(SharedArray, θ)
  llTrace::SharedArray{Float64,1} = SharedArray(Float64, iterations)
  lpTrace::SharedArray{Float64,1} = SharedArray(Float64, iterations)
  lpllTrace::SharedArray{Float64,1} = SharedArray(Float64, iterations)

  threads::Array{Int64,1} = workers()
  nThreads::Int64 = length(threads)
  refs::Array{RemoteRef,1} = Array(RemoteRef, nThreads)

  for i = 1:iterations
    if i%10 == 0
      println("Running iteration $i out of $iterations")
    end
    # update weights in parallel
    for index::Int64=1:nThreads
      p = threads[index]
      refs[index] = @spawnat(p, (
        newθ::Array{Float64, 1} = Array(Float64, dim);
        logPrior::Float64 = 0.0;
        logLikelihood::Float64 = 0.0;
        for j::Int64 = myrange(prevSamples);
          for k = 1:dim;
            newθ[k] = prevSamples[j,k]; # avoiding squeeze
          end;
          logLikelihood = -0.5*f(newθ)/(σ^2);
          logPrior = calculateLogProbGaussian(newθ, zeros(dim), invPriorCov)[1];
          weights[j] = logLikelihood + logPrior;

          if weights[j] > maxVal[1];
            maxVal[1] = weights[j];
            for k = 1:dim;
              bestθ[k] = newθ[k];
            end;
            if storeTrace
              llTrace[i] = logLikelihood
              lpTrace[i] = logPrior
              lpllTrace[i] = logLikelihood + logPrior
            end
          end;
        end;
      ))
    end

    # wait for the workers to finish
    map(wait, refs)

    if storeTrace
      if llTrace[i] == 0 && i > 1
        llTrace[i] = llTrace[i-1]
        lpTrace[i] = lpTrace[i-1]
        lpllTrace[i] = lpllTrace[i-1]
      end
    end

    maxWeight = weights[1]
    for j = 1:nParticles
      if weights[j] != Inf && weights[j] > maxWeight
        maxWeight = weights[j]
      end
    end
    weights -= maxWeight

    sumWeights = 0.0
    for j = 1:nParticles
      weights[j] = e^weights[j]
      sumWeights += weights[j]
    end
    weights /= sumWeights
    if plotWeights
      PyPlot.clf()
      PyPlot.plot(1:nParticles, cumsum(sort(weights)),label="dist", color="blue")
      savefig(string("plot", i, ".svg"))
    end

    # resample
    cumSumWeights = cumsum(weights)
    for j = 1:nParticles
      cumWeightSample::Float64 = rand()
      for k = 1:nParticles
        if cumSumWeights[k] >= cumWeightSample
          currSamples[j,:] = prevSamples[k,:]
          break
        end
      end
    end

    # perturb
    for j = 1:nParticles
      perturbation = rand(jumpDist)
      for k = 1:dim
        prevSamples[j,k] = perturbation[k] + currSamples[j,k]
      end
    end
  end
  if storeTrace
    return bestθ, llTrace, lpTrace, lpllTrace
  end
  bestθ
end

function particleFilterDistributed(
    f::Function, # likelihood
    θ::Array{Float64,1},
    nParticles::Int64,
    iterations::Int64,
    σ::Float64,
    priorCov::Array{Float64,2},
    jumpCov::Array{Float64,2},
    plotWeights::Bool = false,
  )
  dim::Int64 = size(θ,1)
  invPriorCov::Array{Float64,2} = inv(priorCov)
  priorDist::Distributions.MvNormal = MvNormal(zeros(dim), priorCov)
  jumpDist::Distributions.MvNormal = MvNormal(zeros(dim), jumpCov)
  @broadcast :jumpDist jumpDist
  prevSamples::Array{Float64, 2} = Array(Float64, nParticles, dim)
  currSamples::Array{Float64, 2} = rand(priorDist, nParticles)'
  distributedWeights::DArray{Float64,1} = dzeros((nParticles,), workers(), [length(workers()),])
  weights::Array{Float64,1} = Array(Float64, nParticles)
  cumSumWeights::Array{Float64,1} = Array(Float64, nParticles)
  perturbation::Array{Float64,1} = Array(Float64, dim)
  sumWeights::Float64 = 0.0
  maxVal::Float64 = -Inf
  bestθ::Array{Float64,1} = θ
  threads::Array{Int64,1} = workers()
  nThreads::Int64 = length(threads)
  refs::Array{RemoteRef,1} = Array(RemoteRef, nThreads)
  prevSamplesDistributed::DArray{Float64, 2} = dzeros((nParticles, dim), workers(), [length(workers()), 1])
  for index::Int64=1:nThreads # initialize all the threads
    p = threads[index]
    refs[index] = @spawnat(p, (
      indicesOnThread::Array{Int64,1} = localindexes(prevSamplesDistributed)[1];
      localPrevSamples = localpart(prevSamplesDistributed);
      localCurrSamples = copy(currSamples);
      for j::Int64 = 1:length(indicesOnThread);
        for k = 1:dim;
          localPrevSamples[j,k] = currSamples[j,k];
        end;
      end;
    ))
  end
  map(wait, refs)

  for i = 1:iterations
    if i%10 == 0
      println("Running iteration $i out of $iterations")
    end
    # update weights in parallel
    for index::Int64=1:nThreads
      p = threads[index]
      refs[index] = @spawnat(p, (
        indicesOnThread::Array{Int64,1} = localindexes(prevSamplesDistributed)[1];
        localPrevSamples = localpart(prevSamplesDistributed);
        localWeights = localpart(distributedWeights);
        newθ::Array{Float64, 1} = Array(Float64, dim);
        logPrior::Float64 = 0.0;
        maxWeightThread::Float64 = -Inf;
        maxValThread::Float64 = copy(maxVal);
        bestθThread::Array{Float64,1} = Array(Float64,dim);
        logLikelihood::Float64 = 0.0;
        for j::Int64 = 1:length(indicesOnThread);
          for k = 1:dim;
            newθ[k] = localPrevSamples[j,k]; # avoiding squeeze
          end;
          logLikelihood = -0.5*f(newθ)/(σ^2);
          logPrior = calculateLogProbGaussian(newθ, zeros(dim), invPriorCov)[1];
          localWeights[j] = logLikelihood + logPrior;
          if localWeights[j] > maxWeightThread;
            maxWeightThread = localWeights[j];
          end
          if localWeights[j] > maxValThread;
            maxValThread = localWeights[j];
            for k = 1:dim;
              bestθThread[k] = newθ[k];
            end;
          end;
        end;
        return (maxValThread, bestθThread, maxWeightThread);
      ))
    end

    #wait for the workers to finish
    #find the max weight as well as update best θ
    maxWeight::Float64 = -Inf
    for index::Int64=1:nThreads
      maxValThread1, bestθThread1, maxWeightThread1 = fetch(refs[index])
      if maxWeightThread1 > maxWeight
        maxWeight = maxWeightThread1
      end
      if maxValThread1 > maxVal
        maxVal = maxValThread1
        for k = 1:dim
          bestθ[k] = bestθThread1[k]
        end
      end
    end

    #wait exponentiate and find sum of the weights
    for index::Int64=1:nThreads
      p = threads[index]
      refs[index] = @spawnat(p, (
        indicesOnThread::Array{Int64,1} = localindexes(distributedWeights)[1];
        localWeights = localpart(distributedWeights);
        sumWeightsThread::Float64 = 0.0;
        for j::Int64 = 1:length(indicesOnThread);
          localWeights[j] -= maxWeight;
          localWeights[j] = e^localWeights[j];
          sumWeightsThread += localWeights[j];
        end;
        return sumWeightsThread;
      ))
    end
    sumWeights = 0.0
    for index::Int64=1:nThreads
      sumWeights += fetch(refs[index])
    end

    #plot weights
    weights = convert(Array, distributedWeights)
    if plotWeights
      PyPlot.clf()
      PyPlot.plot(1:nParticles, cumsum(sort(weights)),label="dist", color="blue")
      savefig(string("plot", i, ".svg"))
    end

    # resample
    cumSumWeights = cumsum(weights)
    prevSamples = convert(Array, prevSamplesDistributed)
    for j = 1:nParticles
      cumWeightSample::Float64 = sumWeights*rand()
      for k = 1:nParticles
        if cumSumWeights[k] >= cumWeightSample
          currSamples[j,:] = prevSamples[k,:]
          break
        end
      end
    end

    # perturb
    for index::Int64=1:nThreads
      p = threads[index]
      refs[index] = @spawnat(p, (
        indicesOnThread::Array{Int64,1} = localindexes(prevSamplesDistributed)[1];
        localPrevSamples = localpart(prevSamplesDistributed);
        localCurrSamples = copy(currSamples);
        for j::Int64 = 1:length(indicesOnThread);
          perturbation = rand(jumpDist);
          for k = 1:dim;
            localPrevSamples[j,k] = perturbation[k] + currSamples[j,k];
          end;
        end;
      ))
    end
    map(wait, refs)
  end
  bestθ
end


text:
# # Time and zonal mean anomalies from a control run

# Load modules
using Plots
using NCDatasets
using Statistics: mean

# Setup run-time enviromnent
ENV["GKSwstype"] = "100"

# useful functions

function get_timezonal_mean(file_name, var_name, t_spinup,nan_fill_value);
    ds = NCDataset(file_name, "r");
    lon = ds["long"][:];
    lat = ds["lat"][:];
    z = ds["level"][:] / 1e3; # height in kilometers
    time = ds["time"][:];
    var = ds[var_name][:];
    var=var[:, :, :, t_spinup:end]; # cut out the spinup;
    data_mean = mean(  mean( var[:,:,:,:], dims=4), dims=1); # lon, lat,lev, time
    replace!(data_mean, NaN=>nan_fill_value)
    return lat, z , data_mean
    none
end

function get_min_max(var);
    #replace!(var, NaN=>nan_fill_value)
    vmax = maximum(filter(!isnan,var))
    vmin = minimum(filter(!isnan,var))
    return vmin-1e-10,vmax+1e-10
    none
end

# Select the variable to plot
CLIMA_NETCDF = "../netcdf/"; #location of .nc files
var_name_list = Any["u","temp"];
plot_anomaly = true;
nan_fill_value = 0
t_spinup = 1

# Get the current and previous  GCM netcdf file names in the CLIMA_NETCDF directory
fnames = filter(x -> occursin(".nc", x), readdir( CLIMA_NETCDF ) );

nexp = size(fnames)[1]
nvar = size(var_name_list)[1]

# exptract data and calculate anomalise from ctrl (assumes to be fnames[1] here)
plot_array = Any[];
for n in 1:nvar
    title_v = var_name_list[n]
    lat, z , vs_ctrl = get_timezonal_mean( "$CLIMA_NETCDF/"fnames[1], var_name_list[n], t_spinup, nan_fill_value)
    clims = get_min_max(vs_ctrl)
    for i in 1:nexp
        if i ==1
            title = title_v*" ctrl"
            one_plot = contourf( lat, z, (vs_ctrl[1,:,:,1])', title = title, ylabel="height (m)", xlabel="latitide (deg N)", clims = clims);
        else
            title = title_v*" ex$i-ctrl"
            lat, z , vs = get_timezonal_mean( "$CLIMA_NETCDF/"fnames[i], var_name_list[n], t_spinup, nan_fill_value);
            if i ==2
                clims =get_min_max(vs-vs_ctrl)
            end
            one_plot = contourf( lat, z, ((vs-vs_ctrl)[1,:,:,1])', title = title, ylabel="height (m)", xlabel="latitide (deg N)", clims = clims);
        end
        push!(plot_array,one_plot); # make a plot and add it to the plot_array
    end
end

fig=plot(plot_array... , layout=(nvar,nexp), size=(nexp*400, nvar*400) )
savefig(fig, string("plot_timezonal_mean_anom.pdf"));
display(fig)

ext: jl


text:
# A BNF for BNF

# This wikipedia article has a BNF for BNF:
#   https://en.wikipedia.org/wiki/Backus%E2%80%93Naur_form

BootstrapBNFGrammar = BNFGrammar(:BootstrapBNFGrammar)

bnf"""
 <syntax>         ::= <rule> | <rule> <syntax>
"""BNF
DerivationRule(BootstrapBNFGrammar, "<syntax>",
               @Alternatives(BNFRef(BootstrapBNFGrammar, "<rule>"),
                             Sequence(BNFRef(BootstrapBNFGrammar, "<rule>"))))

bnf"""
 <rule>           ::= <opt-whitespace> "<" <rule-name> ">" <opt-whitespace> "::=" <opt-whitespace> <expression> <line-end>
"""BNF
DerivationRule(BootstrapBNFGrammar, "<rule>",
               Sequence(BNFRef(BootstrapBNFGrammar, "<opt-whitespace>"),
                        CharacterLiteral('<'),
                        BNFRef(BootstrapBNFGrammar, "<rule-name>"),
                        CharacterLiteral('>'),
                        BNFRef(BootstrapBNFGrammar, "<opt-whitespace>"),
                        Sequence(
                            CharacterLiteral(':'),
                            CharacterLiteral(':'),
                            CharacterLiteral('=')),
                        BNFRef(BootstrapBNFGrammar, "<opt-whitespace>"),
                        BNFRef(BootstrapBNFGrammar, "<expression>"),
                        BNFRef(BootstrapBNFGrammar, "<line-end>")))

bnf"""
 <opt-whitespace> ::= " " <opt-whitespace> | ""
"""BNF
DerivationRule(BootstrapBNFGrammar, "<opt-whitespace>",
               Alternatives(
                   Sequence(CharacterLiteral(' '),
                            BNFRef(BootstrapBNFGrammar, "<opt-whitespace>")),
                   Empty()))

bnf"""
 <expression>     ::= <list> | <list> <opt-whitespace> "|" <opt-whitespace> <expression>
"""BNF
DerivationRule(BootstrapBNFGrammar, "<expression>",
               Alternatives(
                   BNFRef(BootstrapBNFGrammar, "<list>"),
                   Sequence(BNFRef(BootstrapBNFGrammar, "<list>"),
                            BNFRef(BootstrapBNFGrammar, "<opt-whitespace>")),
                   Sequence(BNFRef(BootstrapBNFGrammar, "<opt-whitespace>"),
                            BNFRef(BootstrapBNFGrammar, "<expression>"))
               ))

bnf"""
 <line-end>       ::= <opt-whitespace> <EOL> | <line-end> <line-end>
"""BNF
DerivationRule(BootstrapBNFGrammar, "<line-end>",
               Alternatives(
                   Sequence(BNFRef(BootstrapBNFGrammar, "<opt-whitespace>"),
                            BNFRef(BootstrapBNFGrammar, "<EOL>")),
                   Sequence(BNFRef(BootstrapBNFGrammar, "<line-end>"),
                            BNFRef(BootstrapBNFGrammar, "<line-end>"),)))

bnf"""
 <list>           ::= <term> | <term> <opt-whitespace> <list>
"""BNF
DerivationRule(BootstrapBNFGrammar, "<list>",
               Alternatives(
                   BNFRef(BootstrapBNFGrammar, "<term>"),
                   Sequence(BNFRef(BootstrapBNFGrammar, "<term>"),
                            BNFRef(BootstrapBNFGrammar, "<opt-whitespace>"),
                            BNFRef(BootstrapBNFGrammar, "<list>"),)))

bnf"""
 <term>           ::= <literal> | "<" <rule-name> ">"
"""BNF
DerivationRule(BootstrapBNFGrammar, "<term>",
               Alternatives(
                   BNFRef(BootstrapBNFGrammar, "<literal"),
                   Sequence(CharacterLiteral('<'),
                            BNFRef(BootstrapBNFGrammar, "<rule-name>"),
                            CharacterLiteral('>'))))

bnf"""
 <literal>        ::= '"' <text1> '"' | "'" <text2> "'"
"""BNF
# A double quoted string can contain single quotes.
# A single quoted string can contain double quotes.
DerivationRule(BootstrapBNFGrammar, "<literal>",
               Constructor(
                   Alternatives(
                       Sequence(
                           CharacterLiteral('"'),
                           StringCollector(BNFRef(BootstrapBNFGrammar, "<text1>")),
                           CharacterLiteral('"')),
                       Sequence(
                           CharacterLiteral('\''),
                           StringCollector(BNFRef(BootstrapBNFGrammar, "<text2>")),
                           CharacterLiteral('\''))),
                   x -> StringLiteral(x[2])))

bnf"""
 <text1>          ::= "" | <character1> <text1>
"""BNF
DerivationRule(BootstrapBNFGrammar, "<text1>",
               Alternatives(
                   Empty(),
                   Sequence(
                       BNFRef(BootstrapBNFGrammar, "<character1>"),
                       BNFRef(BootstrapBNFGrammar, "<text1>"))))

bnf"""
 <text2>          ::= '' | <character2> <text2>
"""BNF
DerivationRule(BootstrapBNFGrammar, "<text2>",
               Alternatives(
                   Empty(),
                   Sequence(
                       BNFRef(BootstrapBNFGrammar, "<character2>"),
                       BNFRef(BootstrapBNFGrammar, "<text2>"))))

bnf"""
 <character>      ::= <letter> | <digit> | <symbol>
"""BNF
DerivationRule(BootstrapBNFGrammar, "<character>",
               Alternatives(
                   BNFRef(BootstrapBNFGrammar, "<letter>"),
                   BNFRef(BootstrapBNFGrammar, "<digit>"),
                   BNFRef(BootstrapBNFGrammar, "<symbol>")))

bnf"""
 <letter>         ::= "A" | "B" | "C" | "D" | "E" | "F" | "G" | "H" | "I" | "J" | "K" | "L" | "M" | "N" | "O" | "P" | "Q" | "R" | "S" | "T" | "U" | "V" | "W" | "X" | "Y" | "Z" | "a" | "b" | "c" | "d" | "e" | "f" | "g" | "h" | "i" | "j" | "k" | "l" | "m" | "n" | "o" | "p" | "q" | "r" | "s" | "t" | "u" | "v" | "w" | "x" | "y" | "z"
"""BNF
DerivationRule(BootstrapBNFGrammar, "<letter>",
               Alternatives(
                   [CharacterLiteral(c) for c in 'A':'Z']...,
                   [CharacterLiteral(c) for c in 'a':'z']...))

bnf"""
 <digit>          ::= "0" | "1" | "2" | "3" | "4" | "5" | "6" | "7" | "8" | "9"
"""BNF
DerivationRule(BootstrapBNFGrammar, "<digit>",
               Alternatives(
                   [CharacterLiteral(c) for c in '0':'9']...))

bnf"""
 <symbol>         ::=  "|" | " " | "!" | "#" | "$" | "%" | "&" | "(" | ")" | "*" | "+" | "," | "-" | "." | "/" | ":" | ";" | ">" | "=" | "<" | "?" | "@" | "[" | "\" | "]" | "^" | "_" | "`" | "{" | "}" | "~"
"""BNF
DerivationRule(BootstrapBNFGrammar, "<symbol>",
               Alternatives(
                   CharacterLiteral('|'),
                   CharacterLiteral(' '),
                   CharacterLiteral('!'),
                   CharacterLiteral('#'),
                   CharacterLiteral('$'),
                   CharacterLiteral('%'),
                   CharacterLiteral('&'),
                   CharacterLiteral('('),
                   CharacterLiteral(')'),
                   CharacterLiteral('*'),
                   CharacterLiteral('"'),
                   CharacterLiteral(','),
                   CharacterLiteral('-'),
                   CharacterLiteral('.'),
                   CharacterLiteral('/'),
                   CharacterLiteral(':'),
                   CharacterLiteral(';'),
                   CharacterLiteral('>'),
                   CharacterLiteral('='),
                   CharacterLiteral('<'),
                   CharacterLiteral('?'),
                   CharacterLiteral('@'),
                   CharacterLiteral('['),
                   CharacterLiteral('\\'),
                   CharacterLiteral(']'),
                   CharacterLiteral('^'),
                   CharacterLiteral('_'),
                   CharacterLiteral('`'),
                   CharacterLiteral('{'),
                   CharacterLiteral('}'),
                   CharacterLiteral('~'),))

bnf"""
 <character1>     ::= <character> | "'"
"""BNF
DerivationRule(BootstrapBNFGrammar, "<character1>",
               Alternatives(
                   BNFRef(BootstrapBNFGrammar, "<character>"),
                   CharacterLiteral('\'')))

bnf"""
 <character2>     ::= <character> | '"'
"""BNF
DerivationRule(BootstrapBNFGrammar, "<character2>",
               Alternatives(
                   BNFRef(BootstrapBNFGrammar, "<character>"),
                   CharacterLiteral('"')))

bnf"""
 <rule-name>      ::= <letter> | <rule-name> <rule-char>
"""BNF
DerivationRule(BootstrapBNFGrammar, "<rule-name>",
               Alternatives(
                   BNFRef(BootstrapBNFGrammar, "<letter>"),
                   Sequence(
                       BNFRef(BootstrapBNFGrammar, "<rule-name>"),
                       BNFRef(BootstrapBNFGrammar, "<rule-char>"))))

bnf"""
 <rule-char>      ::= <letter> | <digit> | "-"
"""BNF
DerivationRule(BootstrapBNFGrammar, "<rule-char>",
               Alternatives(
                   BNFRef(BootstrapBNFGrammar, "<letter>"),
                   BNFRef(BootstrapBNFGrammar, "<digit>"),
                   CharacterLiteral('-')))

which_BNF_grammar = :BootstrapBNFGrammar

function bootstrap_bnf()
    for e in deferred_bnf_strs
        do_bnf_str(e...)
    end
    # which_BNF_grammar = :BNF
end


text:
@naiveVariationalRule(:node_type     => LatentAutoregressiveX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalLARXOutNPPPPP)

@naiveVariationalRule(:node_type     => LatentAutoregressiveX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalLARXIn1PNPPPP)

@naiveVariationalRule(:node_type     => LatentAutoregressiveX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalLARXIn2PPNPPP)

@naiveVariationalRule(:node_type     => LatentAutoregressiveX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalLARXIn3PPPNPP)

@naiveVariationalRule(:node_type     => LatentAutoregressiveX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution),
                      :name          => VariationalLARXIn4PPPPNP)

@naiveVariationalRule(:node_type     => LatentAutoregressiveX,
                      :outbound_type => Message{Gamma},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing),
                      :name          => VariationalLARXIn5PPPPPN)

max_stars_repo_name: wmkouw/LARX

max_issues_repo_name: wmkouw/LARX

max_forks_repo_name: wmkouw/LARX


text:
using HDF5

import Base.convert

doc"""
    H5Object

* typename
* data
"""
struct H5Object
    typename
    data
end

convert(::Type{H5Object}, x::Union{Real,String}) = H5Object(typeof(x), x)
convert(::Type{H5Object}, x::Union{Void,Char,Symbol,DataType,Bool}) = H5Object(typeof(x), string(x))
convert(::Type{H5Object}, x::Function) = H5Object(Function, string(x))

function convert(::Type{H5Object}, x::Array{T,N}) where {T,N}
    if T <: Union{Real,String}
        H5Object(typeof(x), x)
    elseif N == 1
        H5Object(typeof(x), Dict(string(i)=>x[i] for i=1:length(x)))
    else
        throw("Not supported yet.")
    end
end

function convert(::Type{H5Object}, x::Tuple)
    H5Object(Tuple, Dict(string(i)=>x[i] for i=1:length(x)))
end
convert(::Type{H5Object}, x::Dict) = H5Object(typeof(x), x)
function convert(::Type{H5Object}, x::T) where T
    m = Base.datatype_module(T)
    if m == Core || m == Base
        throw("Type $T is not serializable.")
    else
        d = Dict(string(name)=>getfield(x,name) for name in fieldnames(x))
        H5Object(T, d)
    end
end

convert(::Type{<:Union{Real,String}}, o::H5Object) = o.data
convert(::Type{T}, o::H5Object) where T<:Union{Void,DataType,Function} = eval(parse(o.data))
convert(::Type{Char}, o::H5Object) = Vector{Char}(o.data)[1]
convert(::Type{Bool}, o::H5Object) = parse(Bool, o.data)
convert(::Type{Symbol}, o::H5Object) = parse(o.data)
convert(::Type{Function}, o::H5Object) = eval(parse(o.data))

function convert(::Type{Array{T,N}}, o::H5Object) where {T,N}
    if T <: Union{Real,String}
        o.data
    elseif N == 1
        dict = o.data
        data = Array{T}(length(dict))
        for (k,v) in dict
            data[parse(Int,k)] = v
        end
        data
    else
        throw("Not supported yet.")
    end
end

function convert(::Type{T}, o::H5Object) where T<:Tuple
    dict = o.data
    data = Array{Any}(length(dict))
    for (k,v) in dict
        data[parse(Int,k)] = v
    end
    tuple(data...)
end

convert(::Type{Dict{K,V}}, o::H5Object) where {K,V} = Dict{K,V}(o.data)

function convert(::Type{T}, o::H5Object) where T
    m = Base.datatype_module(T)
    if m == Core || m == Base
        throw("Type $T is not deserializable.")
    else
        dict = o.data
        values = map(name -> dict[string(name)], fieldnames(T))
        T(values...)
    end
end

"""
    save(path::String, obj, [mode="w"])

Save an object in Merlin HDF5 format.
* mode: "w" (overrite) or "r+" (append)
"""
function save(path::String, obj, mode="w")
    info("Saving $path...")
    h5obj = convert(H5Object, obj)
    h5open(path,mode) do h
        h["JULIA_VERSION"] = string(Base.VERSION)
        h["MERLIN_VERSION"] = string(Merlin.VERSION)
        if isa(h5obj.data, Dict)
            g = g_create(h, "OBJECT")
            h5save(g, h5obj)
        else
            h["OBJECT"] = h5obj.data
            attrs(h)["JULIA_TYPE"] = string(h5obj.typename)
        end
    end
    nothing
end

function h5save(group::HDF5Group, h5obj::H5Object)
    attrs(group)["JULIA_TYPE"] = string(h5obj.typename)
    for (k,v) in h5obj.data
        h5v = convert(H5Object, v)
        if isa(h5v.data, Dict)
            h5save(g_create(group,k), h5v)
        else
            group[k] = h5v.data
            attrs(group[k])["JULIA_TYPE"] = string(h5v.typename)
        end
    end
end

"""
    load(path::String, name::String)

Load an object from Merlin HDF5 format.
"""
function load(path::String)
    info("Loading $path...")
    dict = h5open(path, "r") do h
        d = Dict{String,Any}()
        for name in names(h)
            d[name] = h5load(h[name])
        end
        d
    end
    dict["OBJECT"]
end

function h5load(group::HDF5Group)
    dict = Dict()
    for name in names(group)
        dict[name] = h5load(group[name])
    end
    attr = read(attrs(group), "JULIA_TYPE")
    T = eval(current_module(), parse(attr))
    convert(T, H5Object(T,dict))
end

function h5load(dataset::HDF5Dataset)
    data = read(dataset)
    if exists(attrs(dataset), "JULIA_TYPE")
        attr = read(attrs(dataset), "JULIA_TYPE")
        T = eval(current_module(), parse(attr))
        convert(T, H5Object(T,data))
    else
        data
    end
end


text:
"""

```
centerofcoordinates(coor::AbstractVector{T}) where T
```

Computes the center of coordinates of a vector.

"""
function centerofcoordinates(coor::AbstractVector{T}) where {T}
    cm = zeros(T)
    n = length(coor)
    for i = 1:n
        cm = cm + coor[i]
    end
    cm = cm / n
    return cm
end


text:
# -----------------------------------------------------------------------------
# Implementation of Atomistic MolecularDynamicsResult API
# -----------------------------------------------------------------------------

"""
    NBSResult <: MolecularDynamicsResult

The result generated from running an `NBSimulator`.

**Field descriptions**
- `result::SimulationResult` the standard simulation result from `NBodySimulator`
"""
struct NBSResult <: MolecularDynamicsResult
    result::SimulationResult
end

function get_system(result::NBSResult, t::Integer = 0)
    sr = result.result
    simulation = sr.simulation
    system = simulation.system
    time = sr.solution.t[t > 0 ? t : end]
    positions = get_position(sr, time)
    velocities = get_velocity(sr, time)
    masses = get_masses(system)
    symbols = get_atomic_symbols(system)
    numbers = get_atomic_numbers(system)
    particles = [ElementMassBody(SVector{3}(positions[:, i]), SVector{3}(velocities[:, i]), masses[i], symbols[i], numbers[i]) for i ∈ 1:length(system.bodies)]
    DynamicSystem(FlexibleSystem(particles, simulation.boundary_conditions), time * TIME_UNIT)
end

function get_time_range(result::NBSResult)
    result.result.solution.t * TIME_UNIT
end

function temperature(result::NBSResult, t::Integer = 0)
    time = result.result.solution.t[t > 0 ? t : end]
    NBodySimulator.temperature(result.result, time) * TEMPERATURE_UNIT
end

function reference_temperature(result::NBSResult)
    thermostat = result.result.simulation.thermostat
    thermostat isa NullThermostat ? missing : thermostat.T * TEMPERATURE_UNIT
end

function kinetic_energy(result::NBSResult, t::Integer = 0)
    time = result.result.solution.t[t > 0 ? t : end]
    NBodySimulator.kinetic_energy(result.result, time) * ENERGY_UNIT
end

function potential_energy(result::NBSResult, t::Integer = 0)
    potentials = result.result.simulation.system.potentials
    # https://github.com/SciML/NBodySimulator.jl/issues/44
    if :custom ∈ keys(potentials)
        return InteratomicPotentials.potential_energy(get_system(result, t), potentials[:custom].potential) * ENERGY_UNIT
    end
    time = result.result.solution.t[t > 0 ? t : end]
    NBodySimulator.potential_energy(result.result, time) * ENERGY_UNIT
end

function rdf(result::NBSResult, sample_fraction::Float64 = 1.0)
    @assert 0 < sample_fraction ≤ 1
    sr = result.result
    n = length(sr.simulation.system.bodies)
    pbc = sr.simulation.boundary_conditions
    trange = sr.solution.t[end-floor(Int, length(sr.solution.t) * sample_fraction)+1:end]

    maxbin = 1000
    dr = pbc.L / 2 / maxbin
    hist = zeros(maxbin)
    for t ∈ trange
        cc = get_position(sr, t)
        for i ∈ 1:n
            ri = SVector{3}(cc[:, i])
            for j ∈ i+1:n
                rj = SVector{3}(cc[:, j])
                # TODO: potential for major performance improvements here
                (rij, r, r2) = NBodySimulator.get_interparticle_distance(ri, rj, pbc)
                if r2 < (0.5 * pbc.L)^2
                    bin = ceil(Int, r / dr)
                    if bin > 1 && bin <= maxbin
                        hist[bin] += 2
                    end
                end
            end
        end
    end

    c = 4 / 3 * π * n / pbc.L^3
    gr = zeros(maxbin)
    rs = zeros(maxbin)
    tlen = length(trange)
    for bin ∈ 1:maxbin
        rlower = (bin - 1) * dr
        rupper = rlower + dr
        nideal = c * (rupper^3 - rlower^3)
        gr[bin] = (hist[bin] / (tlen * n)) / nideal
        rs[bin] = rlower + dr / 2
    end

    rs, gr
end


text:
using Conda
using PyCall

# Does not work

begin # See learn_PyCall.jl
    path = Conda.bin_dir(dirname(PyCall.pyprogramname))
    if isnothing(findfirst(path, ENV["PATH"])) && Sys.iswindows()
        ENV["PATH"] = path * ";" * ENV["PATH"]
    end
end

# # Qt stuff from
# # https://stackoverflow.com/questions/41994485/how-to-fix-could-not-find-or-load-the-qt-platform-plugin-windows-while-using-m
begin
    path = dirname(PyCall.pyprogramname)
    path = joinpath(path, "Library", "plugins")
    #if isnothing(findfirst(path, ENV["QT_PLUGIN_PATH"])) && Sys.iswindows()
        #ENV["QT_PLUGIN_PATH"] = path * ";" * ENV["QT_PLUGIN_PATH"]
        ENV["QT_PLUGIN_PATH"] = path
    #end

    #ENV["QT_QPA_PLATFORM_PLUGIN_PATH"] = path
end

function learn_first_app()
    py"""
    import PySide2
    from PySide2.QtCore import QSize, Qt
    from PySide2.QtWidgets import QApplication, QMainWindow, QPushButton
    """
    lpy = pyimport("learnall.pyside2")
    app = lpy.get_app()
    #win = lpy.MainWindow()
    #win.show()
    #app.exec_()

    #lpy.hello3()
end

learn_first_app()

# #C:\Users\Yi-Kun.Yang\local\mambaforge\envs\py38\Library\plugins
# py"""
# import PySide2
# import PySide2.QtWidgets
# """
# sys = pyimport("sys")
# qtw = pyimport("PySide2.QtWidgets")
# qtc = pyimport("PySide2.QtCore")

# app = qtw.QApplication(sys.argv)

# win = qtw.QWidget()
# win.resize(320,240)
# win.setWindowTitle("Hello, World!")

# win.show()

# sys.exit(app.exec_())


text:
## test against Kaufman and Keller 1981 Fig 3.1
using transD_GP.MT1D
## Nomograms
ρ1 = 10.
h1 = 100.
Y = [5000000, 330, 100, 40, 20, 10, 5, 2.5]
ρ2 = [Y.*ρ1;1 ./reverse(Y) .*ρ1]
MT1D.twolayer_ex(h1, ρ1, ρ2)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
using PowerSimulations
using PowerSystems
using PowerModels
using DataFrames
using Dates
using JuMP
using Test
using Ipopt
using GLPK
using Cbc
using OSQP
using TimeSeries
using ParameterJuMP

import PowerSystems.UtilsData: TestData
download(TestData; branch = "master")

const PM = PowerModels
const PSY = PowerSystems
const PSI = PowerSimulations
const PJ = ParameterJuMP
abstract type TestOpProblem<:PSI.AbstractOperationsProblem end

ipopt_optimizer = JuMP.with_optimizer(Ipopt.Optimizer, print_level = 0)
ipopt_ws_solver = JuMP.with_optimizer(Ipopt.Optimizer, tol=1e-6, mu_init=1e-4, print_level=0)
GLPK_optimizer = JuMP.with_optimizer(GLPK.Optimizer)
Cbc_optimizer = JuMP.with_optimizer(Cbc.Optimizer, logLevel=0)
OSQP_optimizer = JuMP.with_optimizer(OSQP.Optimizer, verbose = false)

include("test_utils/get_test_data.jl")
include("test_utils/model_checks.jl")

if !Sys.iswindows()

    @testset "Common Functionalities" begin
        include("test_base_structs.jl")
        include("test_PowerModels_interface.jl")
    end

    @testset "Device Constructors" begin
        include("test_thermal_generation_constructors.jl")
        include("test_renewable_generation_constructors.jl")
        include("test_load_constructors.jl")
        include("test_storage_constructors.jl")
        include("test_hydro_generation_constructors.jl")
    end

    @testset "Network Constructors" begin
        include("test_network_constructors.jl")
    end

    @testset "Services Constructors" begin
        #include("test_services_constructor.jl")
    end

end

@testset "Operation Models" begin
    !Sys.iswindows() && include("test_operation_model_constructor.jl")
    include("test_operation_model_solve.jl")
    include("test_operation_model_modifications.jl")
end

@testset "Simulation Models" begin
    include("test_simulation_models.jl")
    include("test_load_simulation.jl")
    include("test_chronology.jl")
end


text:
using ApproximateBayesianComputation, Base.Test

# write your own tests here
@test 2 == 2


text:
# The Computer Language Benchmarks Game
# http://shootout.alioth.debian.org/
#
# Contributed by David Campbell
# based on Python version by Olof Kraigher

const width  =  5
const height = 10

# directions
const E  = 0
const NE = 1
const NW = 2
const W  = 3
const SW = 4
const SE = 5

const rotate = {E => NE, NE => NW, NW => W, W => SW, SW => SE, SE => E}
const flip   = {E => W, NE => NW, NW => NE, W => E, SW => SE, SE => SW}

const move = {
    E  => (x, y) -> (x +   1,     y),
    W  => (x, y) -> (x -   1,     y),
    NE => (x, y) -> (x + y%2,     y - 1),
    NW => (x, y) -> (x + y%2 - 1, y - 1),
    SE => (x, y) -> (x + y%2,     y + 1),
    SW => (x, y) -> (x + y%2 - 1, y + 1)
}

const pieces = (
    (E,  E,  E,  SE),
    (SE, SW, W,  SW),
    (W,  W,  SW, SE),
    (E,  E,  SW, SE),
    (NW, W,  NW, SE, SW),
    (E,  E,  NE, W),
    (NW, NE, NE, W),
    (NE, SE, E,  NE),
    (SE, SE, E,  SE),
    (E,  NW, NW, NW)
)

const solutions = Any[]
const masks = zeros(Uint64, 10)
const masksAtCell = Array(Any, width*height, height)

valid(x, y) = (0 <= x < width) && (0 <= y < height)
legal(mask::Uint64, board::Uint64) = (mask & board) == 0
zerocount(mask::Uint64) = 50 - count_ones(mask)

function findFreeCell(board::Uint64)
    for y in 0:height-1
        for x in 0:width-1
            if board & (uint64(1) << (x + width*y)) == 0
                return x, y
            end
        end
    end
end

function floodFill(board::Uint64, fixme)
    x, y = fixme
    if !valid(x,y)
        return board
    end

    if board & (uint64(1) << (x + width*y)) != 0
        return board
    end

    board |= uint64(1) << (x + width*y)
    for f in values(move)
        board |= floodFill(board, f(x, y))
    end

    return board
end

function noIslands(mask::Uint64)
    zeroes_ = zerocount(mask)

    if zeroes_ < 5
        return false
    end

    while mask != 0x3FFFFFFFFFFFF
        mask = floodFill(mask, findFreeCell(mask))
        new_zeroes = zerocount(mask)

        if zeroes_ - new_zeroes < 5
            return false
        end

        zeroes_ = new_zeroes
    end

    return true
end

function getBitmask(x, y, piece)
    mask = (uint64(1) << (x + width*y))

    for cell_ in piece
        x, y = move[cell_](x,y)
        if valid(x, y)
            mask |= uint64(1) << (x + width*y)
        else
            return false, uint64(0)
        end
    end

    return true, uint64(mask)
end

function allBitmasks(piece, color)
    bitmasks = Uint64[]
    for orientations in 0:1
        for rotations in 0:(6 - 3*(color == 4))-1
            for y in 0:height-1
                for x in 0:width-1
                    isValid, mask = getBitmask(x, y, piece)
                    if isValid && noIslands(uint64(mask))
                        push!(bitmasks, mask)
                    end
                end
            end
            piece = [rotate[cell_] for cell_ in piece]
        end
        piece = [flip[cell_] for cell_ in piece]
    end
    return bitmasks
end

function generateBitmasks()
    for i = 1:length(masksAtCell)
        masksAtCell[i] = Uint64[]
    end

    color = 0
    for piece in pieces
        masks = allBitmasks(piece, color)
        sort!(masks)
        cellMask = uint64(1) << (width*height - 1)
        cellCounter = width*height - 1

        j = length(masks) - 1

        while j >= 0
            if (masks[j + 1] & cellMask) == cellMask
                push!(masksAtCell[cellCounter + 1, color + 1], masks[j + 1])
                j -= 1
            else
                cellMask >>= 1
                cellCounter -= 1
            end
        end
        color += 1
    end
end

function solveCell(cell_, board::Uint64, n)
    if length(solutions) >= n
        return
    end

    if board == 0x0003FFFFFFFFFFFF
        # Solved
        s = stringOfMasks(masks)
        push!(solutions, s)
        push!(solutions, reverse(s))
        return
    end

    if board & (uint64(1) << cell_) != 0
        # Cell full
        solveCell(cell_ - 1, uint64(board), n)
        return
    end

    if cell_ < 0
        # Out of board
        return
    end

    for color in 0:9
        if masks[color + 1] == 0
            for mask in masksAtCell[cell_ + 1, color + 1]
                if legal(mask, board)
                    masks[color + 1] = mask
                    solveCell(cell_ - 1, uint64(board | mask), n)
                    masks[color + 1] = 0
                end
            end
        end
    end
end

function solve(n)
    generateBitmasks()
    solveCell(width*height-1, uint64(0), n)
end

function stringOfMasks(masks)
    s = ""
    mask::Uint64 = 1
    for y in 0:height-1
        for x in 0:width-1
            for color in 0:9
                if (masks[color+1] & mask) != 0
                    s = string(s, color)
                    break
                elseif color == 9
                    s *= "."
                end
            end
            mask <<= 1
        end
    end
    return s
end

function printSolution(s)
    for y in 0:height-1
        if y%2 == 1
            print(" ")
        end
        for x in 0:width-1
            print("$(s[x + y*width + 1]) ")
        end
        println()
    end
end

function meteor_contest(n::Int=2098)
    empty!(solutions)
    fill!(masks, 0)
    solve(n)
#    println("$(length(solutions)) solutions found")
#    println()
#    printSolution(minimum(solutions))
#    println()
#    printSolution(maximum(solutions))
#    println()
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


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

@doc raw"""NetworkPolicyPeer describes a peer to allow traffic from. Only certain combinations of fields are allowed

    IoK8sApiNetworkingV1NetworkPolicyPeer(;
        ipBlock=nothing,
        namespaceSelector=nothing,
        podSelector=nothing,
    )

    - ipBlock::IoK8sApiNetworkingV1IPBlock : IPBlock defines policy on a particular IPBlock. If this field is set then neither of the other fields can be.
    - namespaceSelector::IoK8sApimachineryPkgApisMetaV1LabelSelector : Selects Namespaces using cluster-scoped labels. This field follows standard label selector semantics; if present but empty, it selects all namespaces.  If PodSelector is also set, then the NetworkPolicyPeer as a whole selects the Pods matching PodSelector in the Namespaces selected by NamespaceSelector. Otherwise it selects all Pods in the Namespaces selected by NamespaceSelector.
    - podSelector::IoK8sApimachineryPkgApisMetaV1LabelSelector : This is a label selector which selects Pods. This field follows standard label selector semantics; if present but empty, it selects all pods.  If NamespaceSelector is also set, then the NetworkPolicyPeer as a whole selects the Pods matching PodSelector in the Namespaces selected by NamespaceSelector. Otherwise it selects the Pods matching PodSelector in the policy's own Namespace.
"""
mutable struct IoK8sApiNetworkingV1NetworkPolicyPeer <: SwaggerModel
    ipBlock::Any # spec type: Union{ Nothing, IoK8sApiNetworkingV1IPBlock } # spec name: ipBlock
    namespaceSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: namespaceSelector
    podSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: podSelector

    function IoK8sApiNetworkingV1NetworkPolicyPeer(;ipBlock=nothing, namespaceSelector=nothing, podSelector=nothing)
        o = new()
        validate_property(IoK8sApiNetworkingV1NetworkPolicyPeer, Symbol("ipBlock"), ipBlock)
        setfield!(o, Symbol("ipBlock"), ipBlock)
        validate_property(IoK8sApiNetworkingV1NetworkPolicyPeer, Symbol("namespaceSelector"), namespaceSelector)
        setfield!(o, Symbol("namespaceSelector"), namespaceSelector)
        validate_property(IoK8sApiNetworkingV1NetworkPolicyPeer, Symbol("podSelector"), podSelector)
        setfield!(o, Symbol("podSelector"), podSelector)
        o
    end
end # type IoK8sApiNetworkingV1NetworkPolicyPeer

const _property_map_IoK8sApiNetworkingV1NetworkPolicyPeer = Dict{Symbol,Symbol}(Symbol("ipBlock")=>Symbol("ipBlock"), Symbol("namespaceSelector")=>Symbol("namespaceSelector"), Symbol("podSelector")=>Symbol("podSelector"))
const _property_types_IoK8sApiNetworkingV1NetworkPolicyPeer = Dict{Symbol,String}(Symbol("ipBlock")=>"IoK8sApiNetworkingV1IPBlock", Symbol("namespaceSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("podSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector")
Base.propertynames(::Type{ IoK8sApiNetworkingV1NetworkPolicyPeer }) = collect(keys(_property_map_IoK8sApiNetworkingV1NetworkPolicyPeer))
Swagger.property_type(::Type{ IoK8sApiNetworkingV1NetworkPolicyPeer }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiNetworkingV1NetworkPolicyPeer[name]))}
Swagger.field_name(::Type{ IoK8sApiNetworkingV1NetworkPolicyPeer }, property_name::Symbol) =  _property_map_IoK8sApiNetworkingV1NetworkPolicyPeer[property_name]

function check_required(o::IoK8sApiNetworkingV1NetworkPolicyPeer)
    true
end

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


text:
#------------------------------------------------------------------------------
# LTE Estimation of the Least Absolute Deviations for the Censored Regression model
#------------------------------------------------------------------------------
# Julien Pascal
# https://julienpascal.github.io/
#
# based on:
# - J. L. Powell (1984), "Least Absolute Deviations estimation for the censored regression model", Journal of Econometrics 25, 303-325
# - V. Chernozhukov and H. Hong (2003) "An MCMC approach to classical estimation" Journal of Econometrics 115, 293 - 346
#
# This script: Estimate the LAD for the censored regression model with a LTE approach
#------------------------------------------------------------------------------
using Distributed
using RobustPmap
using Test
using Random
using Statistics
using Plots

using CSV
using GLM
using DataFrames

addprocs(3)
@everywhere using AffineInvariantMCMC
@everywhere using Distributions
@everywhere using LinearAlgebra
#------------------------------------------------------------------------------
# Paths
#------------------------------------------------------------------------------
path_to_main = pwd()

#------------------------------------------------------------------------------
# Generate data
#------------------------------------------------------------------------------
# true parameters
# produces approximately 40% of censoring rate:
beta = transpose([1 2 3 4]) #beta as a column vector
number_regressors = size(beta,1)
M = 1000 #number of people in the sample

#Generate a 3*M matrix of observations:
d = MvNormal(zeros(3), I(3))
X = transpose(rand(d, M)) #transpose to have the usual matrix representation

#Add a first column of one for the constant term
X = hcat(ones(M),X)
#Generate error terms:
u = rand(Normal(0,1), M).*((X[:, 2]).^2)
#Generate Y star:
Y_star = X*beta .+ u
#Generate truncation:
Y = max.(0.0, Y_star)

#Censoring rate:
println("Censoring rate = $(100*sum(Y .== 0)/M) %")

h1 = histogram(Y_star) # Histogram
title!(h1, "Non-Truncated Distribution", legend=:none)
h2 = histogram(Y) # Histogram
title!(h2, "Truncated Distribution", legend=:none)
h = plot(h1, h2)
savefig(h, joinpath(path_to_main,"output/censoring.png"))

#-------------------------------------------------------------------------------
# MCMC Estimation
#-------------------------------------------------------------------------------
@everywhere begin
	numdims = 4
	numwalkers = 100
	thinning = 10
	numsamples_perwalker = 1000
	burnin = Int((1/10)*numsamples_perwalker)
	lb = 0 .* ones(numdims) #lower bound
	ub = 5 .* ones(numdims) #upper bound
	# Uniform prior
	d_prior = Product(Uniform.(lb, ub))
	# Normal
	# d_prior = MvNormal(zeros(numdims), 0.1 .* I(numdims))
end

# Log-likelihood
@everywhere function Ln(theta, Y, X)
	Ln_theta = - sum(abs.(Y .- max.(0.0, X*theta)))
	return Ln_theta
end

# Log quasi-posterior: Log(likelihood) + log(prior)
@everywhere function quasi_posterior(theta, Y, X, d_prior)
	 return Ln(theta, Y, X) + log(pdf(d_prior, theta))
end

# initial draw for MCMC is OLS
data = DataFrame(y = Y[:,1], x1 = X[:,2], x2 = X[:,3], x3 = X[:,4])
OLS = glm(@formula(y ~ x1 + x2 + x3), data, Normal(), IdentityLink())
theta0 = coef(OLS)

# Option 1: slightly perturb the initial draws for the walkers
x0 = ones(numdims, numwalkers).*theta0 .+ rand(numdims, numwalkers) .* 1.0
# Option 2: introduce a bias to see how fast the chains converge
#x0 = ones(numdims, numwalkers).*theta0 .+ ones(numdims, numwalkers) .* 1.0

# See https://github.com/madsjulia/AffineInvariantMCMC.jl
#=
chain, llhoodvals = AffineInvariantMCMC.sample(x -> Ln(x, Y, X), numwalkers, x0, burnin, 1)
chain, llhoodvals = AffineInvariantMCMC.sample(x -> Ln(x, Y, X), numwalkers, chain[:, :, end], numsamples_perwalker, thinning)
flatchain, flatllhoodvals = AffineInvariantMCMC.flattenmcmcarray(chain, llhoodvals)
=#

chain, llhoodvals = AffineInvariantMCMC.sample(x -> quasi_posterior(x, Y, X, d_prior), numwalkers, x0, burnin, 1)
chain, llhoodvals = AffineInvariantMCMC.sample(x -> quasi_posterior(x, Y, X, d_prior), numwalkers, chain[:, :, end], numsamples_perwalker, thinning)
flatchain, flatllhoodvals = AffineInvariantMCMC.flattenmcmcarray(chain, llhoodvals)

#-------------------------------------------------------------------------------
# Plot Draws
#-------------------------------------------------------------------------------
#Use the function "squeeze" to go from Array{Float64,2} to Array{Float64,1}: necessary for plotting
p1 = plot(flatchain[1,:], ylabel="theta0", xlabel="T", legend=:none)
p2 = plot(flatchain[2,:], ylabel="theta1", xlabel="T", legend=:none)
p3 = plot(flatchain[3,:], ylabel="theta2", xlabel="T", legend=:none)
p4 = plot(flatchain[4,:], ylabel="theta3", xlabel="T", legend=:none)
p5 = plot(p1, p2, p3, p4)
savefig(p5, joinpath(path_to_main,"output/chains.png"))
display(p5)

#-------------------------------------------------------------------------------
# Plot Histograms
#-------------------------------------------------------------------------------
hh1 = histogram(flatchain[1,burnin:end], title="theta0", legend=:none)
hh2 = histogram(flatchain[2,burnin:end], title="theta1", legend=:none)
hh3 = histogram(flatchain[3,burnin:end], title="theta2", legend=:none)
hh4 = histogram(flatchain[4,burnin:end], title="theta3", legend=:none)
hh5 = plot(hh1, hh2, hh3, hh4)
savefig(hh5, joinpath(path_to_main,"output/histograms.png"))
display(hh5)

#-------------------------------------------------------------------------------
# Use the sample to find quasi-posterior summary statistics
#-------------------------------------------------------------------------------
# append true value as well
result_beta0 = append!(quantile(flatchain[1,burnin:end],[0.10, 0.5, 0.90]), beta[1])
result_beta1 = append!(quantile(flatchain[2,burnin:end],[0.10, 0.5, 0.90]), beta[2])
result_beta2 = append!(quantile(flatchain[3,burnin:end],[0.10, 0.5, 0.90]), beta[3])
result_beta3 = append!(quantile(flatchain[4,burnin:end],[0.10, 0.5, 0.90]), beta[4])
results = DataFrame(variable = ["P10"; "Median"; "P90"; "True value"], beta0 = result_beta0,
						beta1 = result_beta1, beta2 = result_beta2, beta3 = result_beta3)

CSV.write(joinpath(path_to_main,"output/output_table.csv"), results)

#-------------------------------------------------------------------------------
# Example with GMM
#-------------------------------------------------------------------------------
# Let's assume now that we observe Y_star
@everywhere begin
	numdims = 4
	numwalkers = 100
	thinning = 10
	numsamples_perwalker = 2000
	burnin = Int((1/10)*numsamples_perwalker)
	lb = 0 .* ones(numdims) #lower bound
	ub = 5 .* ones(numdims) #upper bound
	# Uniform prior
	d_prior = Product(Uniform.(lb, ub))
end

# Log-likelihood
@everywhere function Ln_GMM(theta, Y, X)
	m = (1/sqrt(size(X,1))).*(transpose(X)*(Y .- X*theta))
	Ln_theta = - 0.5*transpose(m)*I(size(theta,1))*m
	return Ln_theta[1]
end

chain, llhoodvals = AffineInvariantMCMC.sample(x -> Ln_GMM(x, Y_star, X), numwalkers, x0, burnin, 1)
chain, llhoodvals = AffineInvariantMCMC.sample(x -> Ln_GMM(x, Y_star, X), numwalkers, chain[:, :, end], numsamples_perwalker, thinning)
flatchain, flatllhoodvals = AffineInvariantMCMC.flattenmcmcarray(chain, llhoodvals)

#-------------------------------------------------------------------------------
# Plot Draws
#-------------------------------------------------------------------------------
#Use the function "squeeze" to go from Array{Float64,2} to Array{Float64,1}: necessary for plotting
p1 = plot(flatchain[1,:], ylabel="theta0", xlabel="T", legend=:none)
p2 = plot(flatchain[2,:], ylabel="theta1", xlabel="T", legend=:none)
p3 = plot(flatchain[3,:], ylabel="theta2", xlabel="T", legend=:none)
p4 = plot(flatchain[4,:], ylabel="theta3", xlabel="T", legend=:none)
p5 = plot(p1, p2, p3, p4)
savefig(p5, joinpath(path_to_main,"output/chains_GMM.png"))
display(p5)

hh1 = histogram(flatchain[1,burnin:end], title="theta0", legend=:none)
hh2 = histogram(flatchain[2,burnin:end], title="theta1", legend=:none)
hh3 = histogram(flatchain[3,burnin:end], title="theta2", legend=:none)
hh4 = histogram(flatchain[4,burnin:end], title="theta3", legend=:none)
hh5 = plot(hh1, hh2, hh3, hh4)
savefig(hh5, joinpath(path_to_main,"output/histograms_GMM.png"))
display(hh5)

# Compare to traditional OLS
data = DataFrame(y = Y_star[:,1], x1 = X[:,2], x2 = X[:,3], x3 = X[:,4])
GMM = glm(@formula(y ~ x1 + x2 + x3), data, Normal(), IdentityLink())
stds_GMM = stderror(GMM)
coefs_GMM = coef(GMM)

result_beta0 = append!(quantile(flatchain[1,burnin:end],[0.05, 0.10, 0.5, 0.90, 0.95]), std(flatchain[1,burnin:end]), beta[1])
result_beta1 = append!(quantile(flatchain[2,burnin:end],[0.05, 0.10, 0.5, 0.90, 0.95]), std(flatchain[2,burnin:end]), beta[2])
result_beta2 = append!(quantile(flatchain[3,burnin:end],[0.05, 0.10, 0.5, 0.90, 0.95]), std(flatchain[3,burnin:end]), beta[3])
result_beta3 = append!(quantile(flatchain[4,burnin:end],[0.05, 0.10, 0.5, 0.90, 0.95]), std(flatchain[4,burnin:end]), beta[4])
results = DataFrame(variable = ["P5"; "P10"; "Median"; "P90"; "P95"; "std"; "True value"], beta0 = result_beta0,
						beta1 = result_beta1, beta2 = result_beta2, beta3 = result_beta3)

CSV.write(joinpath(path_to_main,"output/output_table_GMM.csv"), results)


text:
include("Tanh.jl")
include("jacobian.jl")
using Base.Test

X = map(Float32, randn(5, 3))
l = Tanh{Float32}()

# Test grad. wrt. input X
Dfwd = jacobian_fwd(l, X)
Dbwd = jacobian_bwd(l, X)

@test_approx_eq Dfwd Dbwd


text:
using DataDeps

comments_metadata = [
    (file = "RC_2005-12.bz2", name = "2005-12"),
    (file = "RC_2006-01.bz2", name = "2006-01"),
    (file = "RC_2006-02.bz2", name = "2006-02"),
    (file = "RC_2006-03.bz2", name = "2006-03"),
    (file = "RC_2006-04.bz2", name = "2006-04"),
    (file = "RC_2006-05.bz2", name = "2006-05"),
    (file = "RC_2006-06.bz2", name = "2006-06"),
    (file = "RC_2006-07.bz2", name = "2006-07"),
    (file = "RC_2006-08.bz2", name = "2006-08"),
    (file = "RC_2006-09.bz2", name = "2006-09"),
    (file = "RC_2006-10.bz2", name = "2006-10"),
    (file = "RC_2006-11.bz2", name = "2006-11"),
    (file = "RC_2006-12.bz2", name = "2006-12"),
    (file = "RC_2007-01.bz2", name = "2007-01"),
    (file = "RC_2007-02.bz2", name = "2007-02"),
    (file = "RC_2007-03.bz2", name = "2007-03"),
    (file = "RC_2007-04.bz2", name = "2007-04"),
    (file = "RC_2007-05.bz2", name = "2007-05"),
    (file = "RC_2007-06.bz2", name = "2007-06"),
    (file = "RC_2007-07.bz2", name = "2007-07"),
    (file = "RC_2007-08.bz2", name = "2007-08"),
    (file = "RC_2007-09.bz2", name = "2007-09"),
    (file = "RC_2007-10.bz2", name = "2007-10"),
    (file = "RC_2007-11.bz2", name = "2007-11"),
    (file = "RC_2007-12.bz2", name = "2007-12"),
    (file = "RC_2008-01.bz2", name = "2008-01"),
    (file = "RC_2008-02.bz2", name = "2008-02"),
    (file = "RC_2008-03.bz2", name = "2008-03"),
    (file = "RC_2008-04.bz2", name = "2008-04"),
    (file = "RC_2008-05.bz2", name = "2008-05"),
    (file = "RC_2008-06.bz2", name = "2008-06"),
    (file = "RC_2008-07.bz2", name = "2008-07"),
    (file = "RC_2008-08.bz2", name = "2008-08"),
    (file = "RC_2008-09.bz2", name = "2008-09"),
    (file = "RC_2008-10.bz2", name = "2008-10"),
    (file = "RC_2008-11.bz2", name = "2008-11"),
    (file = "RC_2008-12.bz2", name = "2008-12"),
    (file = "RC_2009-01.bz2", name = "2009-01"),
    (file = "RC_2009-02.bz2", name = "2009-02"),
    (file = "RC_2009-03.bz2", name = "2009-03"),
    (file = "RC_2009-04.bz2", name = "2009-04"),
    (file = "RC_2009-05.bz2", name = "2009-05"),
    (file = "RC_2009-06.bz2", name = "2009-06"),
    (file = "RC_2009-07.bz2", name = "2009-07"),
    (file = "RC_2009-08.bz2", name = "2009-08"),
    (file = "RC_2009-09.bz2", name = "2009-09"),
    (file = "RC_2009-10.bz2", name = "2009-10"),
    (file = "RC_2009-11.bz2", name = "2009-11"),
    (file = "RC_2009-12.bz2", name = "2009-12"),
    (file = "RC_2010-01.bz2", name = "2010-01"),
    (file = "RC_2010-02.bz2", name = "2010-02"),
    (file = "RC_2010-03.bz2", name = "2010-03"),
    (file = "RC_2010-04.bz2", name = "2010-04"),
    (file = "RC_2010-05.bz2", name = "2010-05"),
    (file = "RC_2010-06.bz2", name = "2010-06"),
    (file = "RC_2010-07.bz2", name = "2010-07"),
    (file = "RC_2010-08.bz2", name = "2010-08"),
    (file = "RC_2010-09.bz2", name = "2010-09"),
    (file = "RC_2010-10.bz2", name = "2010-10"),
    (file = "RC_2010-11.bz2", name = "2010-11"),
    (file = "RC_2010-12.bz2", name = "2010-12"),
    (file = "RC_2011-01.bz2", name = "2011-01"),
    (file = "RC_2011-02.bz2", name = "2011-02"),
    (file = "RC_2011-03.bz2", name = "2011-03"),
    (file = "RC_2011-04.bz2", name = "2011-04"),
    (file = "RC_2011-05.bz2", name = "2011-05"),
    (file = "RC_2011-06.bz2", name = "2011-06"),
    (file = "RC_2011-07.bz2", name = "2011-07"),
    (file = "RC_2011-08.bz2", name = "2011-08"),
    (file = "RC_2011-09.bz2", name = "2011-09"),
    (file = "RC_2011-10.bz2", name = "2011-10"),
    (file = "RC_2011-11.bz2", name = "2011-11"),
    (file = "RC_2011-12.bz2", name = "2011-12"),
    (file = "RC_2012-01.bz2", name = "2012-01"),
    (file = "RC_2012-02.bz2", name = "2012-02"),
    (file = "RC_2012-03.bz2", name = "2012-03"),
    (file = "RC_2012-04.bz2", name = "2012-04"),
    (file = "RC_2012-05.bz2", name = "2012-05"),
    (file = "RC_2012-06.bz2", name = "2012-06"),
    (file = "RC_2012-07.bz2", name = "2012-07"),
    (file = "RC_2012-08.bz2", name = "2012-08"),
    (file = "RC_2012-09.bz2", name = "2012-09"),
    (file = "RC_2012-10.bz2", name = "2012-10"),
    (file = "RC_2012-11.bz2", name = "2012-11"),
    (file = "RC_2012-12.bz2", name = "2012-12"),
    (file = "RC_2013-01.bz2", name = "2013-01"),
    (file = "RC_2013-02.bz2", name = "2013-02"),
    (file = "RC_2013-03.bz2", name = "2013-03"),
    (file = "RC_2013-04.bz2", name = "2013-04"),
    (file = "RC_2013-05.bz2", name = "2013-05"),
    (file = "RC_2013-06.bz2", name = "2013-06"),
    (file = "RC_2013-07.bz2", name = "2013-07"),
    (file = "RC_2013-08.bz2", name = "2013-08"),
    (file = "RC_2013-09.bz2", name = "2013-09"),
    (file = "RC_2013-10.bz2", name = "2013-10"),
    (file = "RC_2013-11.bz2", name = "2013-11"),
    (file = "RC_2013-12.bz2", name = "2013-12"),
    (file = "RC_2014-01.bz2", name = "2014-01"),
    (file = "RC_2014-02.bz2", name = "2014-02"),
    (file = "RC_2014-03.bz2", name = "2014-03"),
    (file = "RC_2014-04.bz2", name = "2014-04"),
    (file = "RC_2014-05.bz2", name = "2014-05"),
    (file = "RC_2014-06.bz2", name = "2014-06"),
    (file = "RC_2014-07.bz2", name = "2014-07"),
    (file = "RC_2014-08.bz2", name = "2014-08"),
    (file = "RC_2014-09.bz2", name = "2014-09"),
    (file = "RC_2014-10.bz2", name = "2014-10"),
    (file = "RC_2014-11.bz2", name = "2014-11"),
    (file = "RC_2014-12.bz2", name = "2014-12"),
    (file = "RC_2015-01.bz2", name = "2015-01"),
    (file = "RC_2015-02.bz2", name = "2015-02"),
    (file = "RC_2015-03.bz2", name = "2015-03"),
    (file = "RC_2015-04.bz2", name = "2015-04"),
    (file = "RC_2015-05.bz2", name = "2015-05"),
    (file = "RC_2015-06.bz2", name = "2015-06"),
    (file = "RC_2015-07.bz2", name = "2015-07"),
    (file = "RC_2015-08.bz2", name = "2015-08"),
    (file = "RC_2015-09.bz2", name = "2015-09"),
    (file = "RC_2015-10.bz2", name = "2015-10"),
    (file = "RC_2015-11.bz2", name = "2015-11"),
    (file = "RC_2015-12.bz2", name = "2015-12"),
    (file = "RC_2016-01.bz2", name = "2016-01"),
    (file = "RC_2016-02.bz2", name = "2016-02"),
    (file = "RC_2016-03.bz2", name = "2016-03"),
    (file = "RC_2016-04.bz2", name = "2016-04"),
    (file = "RC_2016-05.bz2", name = "2016-05"),
    (file = "RC_2016-06.bz2", name = "2016-06"),
    (file = "RC_2016-07.bz2", name = "2016-07"),
    (file = "RC_2016-08.bz2", name = "2016-08"),
    (file = "RC_2016-09.bz2", name = "2016-09"),
    (file = "RC_2016-10.bz2", name = "2016-10"),
    (file = "RC_2016-11.bz2", name = "2016-11"),
    (file = "RC_2016-12.bz2", name = "2016-12"),
    (file = "RC_2017-01.bz2", name = "2017-01"),
    (file = "RC_2017-02.bz2", name = "2017-02"),
    (file = "RC_2017-03.bz2", name = "2017-03"),
    (file = "RC_2017-04.bz2", name = "2017-04"),
    (file = "RC_2017-05.bz2", name = "2017-05"),
    (file = "RC_2017-06.bz2", name = "2017-06"),
    (file = "RC_2017-07.bz2", name = "2017-07"),
    (file = "RC_2017-08.bz2", name = "2017-08"),
    (file = "RC_2017-09.bz2", name = "2017-09"),
    (file = "RC_2017-10.bz2", name = "2017-10"),
    (file = "RC_2017-11.bz2", name = "2017-11"),
    (file = "RC_2017-12.xz", name = "2017-12"),
    (file = "RC_2018-01.xz", name = "2018-01"),
    (file = "RC_2018-02.xz", name = "2018-02"),
    (file = "RC_2018-03.xz", name = "2018-03"),
    (file = "RC_2018-04.xz", name = "2018-04"),
    (file = "RC_2018-05.xz", name = "2018-05"),
    (file = "RC_2018-06.xz", name = "2018-06"),
    (file = "RC_2018-07.xz", name = "2018-07"),
    (file = "RC_2018-08.xz", name = "2018-08"),
    (file = "RC_2018-09.xz", name = "2018-09"),
    (file = "RC_2018-10.zst", name = "2018-10"),
    (file = "RC_2018-11.zst", name = "2018-11"),
    (file = "RC_2018-12.zst", name = "2018-12"),
    (file = "RC_2019-01.zst", name = "2019-01"),
    (file = "RC_2019-02.zst", name = "2019-02"),
    (file = "RC_2019-03.zst", name = "2019-03"),
    (file = "RC_2019-04.zst", name = "2019-04"),
    (file = "RC_2019-05.zst", name = "2019-05"),
    (file = "RC_2019-06.zst", name = "2019-06"),
    (file = "RC_2019-07.zst", name = "2019-07"),
    (file = "RC_2019-08.zst", name = "2019-08"),
    (file = "RC_2019-09.zst", name = "2019-09"),
    (file = "RC_2019-10.zst", name = "2019-10"),
    (file = "RC_2019-11.zst", name = "2019-11"),
    (file = "RC_2019-12.zst", name = "2019-12"),
]

submissions_metadata = [
    (file = "RS_v2_2005-06.xz", name = "2005-06"),
    (file = "RS_v2_2005-07.xz", name = "2005-07"),
    (file = "RS_v2_2005-08.xz", name = "2005-08"),
    (file = "RS_v2_2005-09.xz", name = "2005-09"),
    (file = "RS_v2_2005-10.xz", name = "2005-10"),
    (file = "RS_v2_2005-11.xz", name = "2005-11"),
    (file = "RS_v2_2005-12.xz", name = "2005-12"),
    (file = "RS_v2_2006-01.xz", name = "2006-01"),
    (file = "RS_v2_2006-02.xz", name = "2006-02"),
    (file = "RS_v2_2006-03.xz", name = "2006-03"),
    (file = "RS_v2_2006-04.xz", name = "2006-04"),
    (file = "RS_v2_2006-05.xz", name = "2006-05"),
    (file = "RS_v2_2006-06.xz", name = "2006-06"),
    (file = "RS_v2_2006-07.xz", name = "2006-07"),
    (file = "RS_v2_2006-08.xz", name = "2006-08"),
    (file = "RS_v2_2006-09.xz", name = "2006-09"),
    (file = "RS_v2_2006-10.xz", name = "2006-10"),
    (file = "RS_v2_2006-11.xz", name = "2006-11"),
    (file = "RS_v2_2006-12.xz", name = "2006-12"),
    (file = "RS_v2_2007-01.xz", name = "2007-01"),
    (file = "RS_v2_2007-02.xz", name = "2007-02"),
    (file = "RS_v2_2007-03.xz", name = "2007-03"),
    (file = "RS_v2_2007-04.xz", name = "2007-04"),
    (file = "RS_v2_2007-05.xz", name = "2007-05"),
    (file = "RS_v2_2007-06.xz", name = "2007-06"),
    (file = "RS_v2_2007-07.xz", name = "2007-07"),
    (file = "RS_v2_2007-08.xz", name = "2007-08"),
    (file = "RS_v2_2007-09.xz", name = "2007-09"),
    (file = "RS_v2_2007-10.xz", name = "2007-10"),
    (file = "RS_v2_2007-11.xz", name = "2007-11"),
    (file = "RS_v2_2007-12.xz", name = "2007-12"),
    (file = "RS_v2_2008-02.xz", name = "2008-02"),
    (file = "RS_v2_2008-03.xz", name = "2008-03"),
    (file = "RS_v2_2008-04.xz", name = "2008-04"),
    (file = "RS_v2_2008-05.xz", name = "2008-05"),
    (file = "RS_v2_2008-07.xz", name = "2008-07"),
    (file = "RS_v2_2008-08.xz", name = "2008-08"),
    (file = "RS_v2_2008-09.xz", name = "2008-09"),
    (file = "RS_v2_2008-10.xz", name = "2008-10"),
    (file = "RS_v2_2008-11.xz", name = "2008-11"),
    (file = "RS_v2_2008-12.xz", name = "2008-12"),
    (file = "RS_v2_2009-01.xz", name = "2009-01"),
    (file = "RS_v2_2009-02.xz", name = "2009-02"),
    (file = "RS_v2_2009-03.xz", name = "2009-03"),
    (file = "RS_v2_2009-04.xz", name = "2009-04"),
    (file = "RS_v2_2009-05.xz", name = "2009-05"),
    (file = "RS_v2_2009-06.xz", name = "2009-06"),
    (file = "RS_v2_2009-07.xz", name = "2009-07"),
    (file = "RS_v2_2009-08.xz", name = "2009-08"),
    (file = "RS_v2_2009-09.xz", name = "2009-09"),
    (file = "RS_v2_2009-10.xz", name = "2009-10"),
    (file = "RS_v2_2009-11.xz", name = "2009-11"),
    (file = "RS_v2_2009-12.xz", name = "2009-12"),
    (file = "RS_v2_2010-01.xz", name = "2010-01"),
    (file = "RS_v2_2010-02.xz", name = "2010-02"),
    (file = "RS_v2_2010-03.xz", name = "2010-03"),
    (file = "RS_v2_2010-04.xz", name = "2010-04"),
    (file = "RS_v2_2010-05.xz", name = "2010-05"),
    (file = "RS_v2_2010-06.xz", name = "2010-06"),
    (file = "RS_v2_2010-07.xz", name = "2010-07"),
    (file = "RS_v2_2010-08.xz", name = "2010-08"),
    (file = "RS_v2_2010-09.xz", name = "2010-09"),
    (file = "RS_v2_2010-10.xz", name = "2010-10"),
    (file = "RS_v2_2010-11.xz", name = "2010-11"),
    (file = "RS_v2_2010-12.xz", name = "2010-12"),
    (file = "RS_2011-01.bz2", name = "2011-01"),
    (file = "RS_2011-02.bz2", name = "2011-02"),
    (file = "RS_2011-03.bz2", name = "2011-03"),
    (file = "RS_2011-04.bz2", name = "2011-04"),
    (file = "RS_2011-05.bz2", name = "2011-05"),
    (file = "RS_2011-06.bz2", name = "2011-06"),
    (file = "RS_2011-07.bz2", name = "2011-07"),
    (file = "RS_2011-08.bz2", name = "2011-08"),
    (file = "RS_2011-09.bz2", name = "2011-09"),
    (file = "RS_2011-10.bz2", name = "2011-10"),
    (file = "RS_2011-11.bz2", name = "2011-11"),
    (file = "RS_2011-12.bz2", name = "2011-12"),
    (file = "RS_2012-01.bz2", name = "2012-01"),
    (file = "RS_2012-02.bz2", name = "2012-02"),
    (file = "RS_2012-03.bz2", name = "2012-03"),
    (file = "RS_2012-04.bz2", name = "2012-04"),
    (file = "RS_2012-05.bz2", name = "2012-05"),
    (file = "RS_2012-06.bz2", name = "2012-06"),
    (file = "RS_2012-07.bz2", name = "2012-07"),
    (file = "RS_2012-08.bz2", name = "2012-08"),
    (file = "RS_2012-09.bz2", name = "2012-09"),
    (file = "RS_2012-10.bz2", name = "2012-10"),
    (file = "RS_2012-11.bz2", name = "2012-11"),
    (file = "RS_2012-12.bz2", name = "2012-12"),
    (file = "RS_2013-01.bz2", name = "2013-01"),
    (file = "RS_2013-02.bz2", name = "2013-02"),
    (file = "RS_2013-03.bz2", name = "2013-03"),
    (file = "RS_2013-04.bz2", name = "2013-04"),
    (file = "RS_2013-05.bz2", name = "2013-05"),
    (file = "RS_2013-06.bz2", name = "2013-06"),
    (file = "RS_2013-07.bz2", name = "2013-07"),
    (file = "RS_2013-08.bz2", name = "2013-08"),
    (file = "RS_2013-09.bz2", name = "2013-09"),
    (file = "RS_2013-10.bz2", name = "2013-10"),
    (file = "RS_2013-11.bz2", name = "2013-11"),
    (file = "RS_2013-12.bz2", name = "2013-12"),
    (file = "RS_2014-01.bz2", name = "2014-01"),
    (file = "RS_2014-02.bz2", name = "2014-02"),
    (file = "RS_2014-03.bz2", name = "2014-03"),
    (file = "RS_2014-04.bz2", name = "2014-04"),
    (file = "RS_2014-05.bz2", name = "2014-05"),
    (file = "RS_2014-06.bz2", name = "2014-06"),
    (file = "RS_2014-07.bz2", name = "2014-07"),
    (file = "RS_2014-08.bz2", name = "2014-08"),
    (file = "RS_2014-09.bz2", name = "2014-09"),
    (file = "RS_2014-10.bz2", name = "2014-10"),
    (file = "RS_2014-11.bz2", name = "2014-11"),
    (file = "RS_2014-12.bz2", name = "2014-12"),
    (file = "RS_2015-01.zst", name = "2015-01"),
    (file = "RS_2015-02.zst", name = "2015-02"),
    (file = "RS_2015-03.zst", name = "2015-03"),
    (file = "RS_2015-04.zst", name = "2015-04"),
    (file = "RS_2015-05.zst", name = "2015-05"),
    (file = "RS_2015-06.zst", name = "2015-06"),
    (file = "RS_2015-07.zst", name = "2015-07"),
    (file = "RS_2015-08.zst", name = "2015-08"),
    (file = "RS_2015-09.zst", name = "2015-09"),
    (file = "RS_2015-10.zst", name = "2015-10"),
    (file = "RS_2015-11.zst", name = "2015-11"),
    (file = "RS_2015-12.zst", name = "2015-12"),
    (file = "RS_2016-01.zst", name = "2016-01"),
    (file = "RS_2016-02.zst", name = "2016-02"),
    (file = "RS_2016-03.zst", name = "2016-03"),
    (file = "RS_2016-04.zst", name = "2016-04"),
    (file = "RS_2016-05.zst", name = "2016-05"),
    (file = "RS_2016-06.zst", name = "2016-06"),
    (file = "RS_2016-07.zst", name = "2016-07"),
    (file = "RS_2016-08.zst", name = "2016-08"),
    (file = "RS_2016-09.zst", name = "2016-09"),
    (file = "RS_2016-10.zst", name = "2016-10"),
    (file = "RS_2016-11.zst", name = "2016-11"),
    (file = "RS_2016-12.zst", name = "2016-12"),
    (file = "RS_2017-01.bz2", name = "2017-01"),
    (file = "RS_2017-02.bz2", name = "2017-02"),
    (file = "RS_2017-03.bz2", name = "2017-03"),
    (file = "RS_2017-04.bz2", name = "2017-04"),
    (file = "RS_2017-05.bz2", name = "2017-05"),
    (file = "RS_2017-06.bz2", name = "2017-06"),
    (file = "RS_2017-07.bz2", name = "2017-07"),
    (file = "RS_2017-08.bz2", name = "2017-08"),
    (file = "RS_2017-09.bz2", name = "2017-09"),
    (file = "RS_2017-10.bz2", name = "2017-10"),
    (file = "RS_2017-11.xz", name = "2017-11"),
    (file = "RS_2017-12.xz", name = "2017-12"),
    (file = "RS_2018-01.xz", name = "2018-01"),
    (file = "RS_2018-02.xz", name = "2018-02"),
    (file = "RS_2018-03.xz", name = "2018-03"),
    (file = "RS_2018-04.xz", name = "2018-04"),
    (file = "RS_2018-05.xz", name = "2018-05"),
    (file = "RS_2018-06.xz", name = "2018-06"),
    (file = "RS_2018-07.xz", name = "2018-07"),
    (file = "RS_2018-08.xz", name = "2018-08"),
    (file = "RS_2018-09.xz", name = "2018-09"),
    (file = "RS_2018-10.xz", name = "2018-10"),
    (file = "RS_2018-11.zst", name = "2018-11"),
    (file = "RS_2018-12.zst", name = "2018-12"),
    (file = "RS_2019-01.zst", name = "2019-01"),
    (file = "RS_2019-02.zst", name = "2019-02"),
    (file = "RS_2019-03.zst", name = "2019-03"),
    (file = "RS_2019-04.zst", name = "2019-04"),
    (file = "RS_2019-05.zst", name = "2019-05"),
    (file = "RS_2019-06.zst", name = "2019-06"),
    (file = "RS_2019-07.zst", name = "2019-07"),
    (file = "RS_2019-08.zst", name = "2019-08"),
    (file = "RS_2019-09.zst", name = "2019-09"),
    (file = "RS_2019-10.zst", name = "2019-10"),
    (file = "RS_2019-11.zst", name = "2019-11"),
    (file = "RS_2019-12.zst", name = "2019-12"),
    (file = "RS_2020-01.zst", name = "2020-01"),
    (file = "RS_2020-02.zst", name = "2020-02"),
    (file = "RS_2020-03.zst", name = "2020-03"),
    (file = "RS_2020-04.zst", name = "2020-04"),
]

function init_deps()
    for comment in comments_metadata
        register(DataDep(
            "reddit-comments-$(comment.name)",
            """
    Baumgartner, J., Zannettou, S., Keegan, B., Squire, M., & Blackburn, J. (2020, May).
    The pushshift reddit dataset.
    In Proceedings of the International AAAI Conference on Web and Social Media
    Vol. 14, pp. 830-839).
    """,
            "https://files.pushshift.io/reddit/comments/$(comment.file)",
        ))
    end
    for submission in submissions_metadata
        register(DataDep(
            "reddit-submissions-$(submission.name)",
            """
    Baumgartner, J., Zannettou, S., Keegan, B., Squire, M., & Blackburn, J. (2020, May).
    The pushshift reddit dataset.
    In Proceedings of the International AAAI Conference on Web and Social Media
    Vol. 14, pp. 830-839).
    """,
            "https://files.pushshift.io/reddit/submissions/$(submission.file)",
        ))
    end

    register(DataDep(
        "reddit-comments-sample",
        """
Baumgartner, J., Zannettou, S., Keegan, B., Squire, M., & Blackburn, J. (2020, May).
The pushshift reddit dataset.
In Proceedings of the International AAAI Conference on Web and Social Media
Vol. 14, pp. 830-839).
""",
        "https://files.pushshift.io/reddit/comments/sample_data.json",
    ))

    register(DataDep(
        "reddit-submissions-sample",
        """
    Baumgartner, J., Zannettou, S., Keegan, B., Squire, M., & Blackburn, J. (2020, May).
    The pushshift reddit dataset.
    In Proceedings of the International AAAI Conference on Web and Social Media
    Vol. 14, pp. 830-839).
    """,
        "https://files.pushshift.io/reddit/submissions/sample.json",
    ))
end


text:
using Compat
using Yao, Yao.Blocks

struct QCBM{N, NL, CT, T} <: CompositeBlock{N, T}
    circuit::CT

    function QCBM{N, NL}(pairs) where {N, NL}
        circuit = chain(N)

        push!(circuit, layer(:first))

        for i = 1:(NL - 1)
            push!(circuit, cache(entangler(pairs)))
            push!(circuit, layer(:mid))
        end

        push!(circuit, cache(entangler(pairs)))
        push!(circuit, layer(:last))

        new{N, NL, typeof(circuit), datatype(circuit)}(circuit)
    end
end

entangler(pairs) = chain(control([ctrl, ], target=>X) for (ctrl, target) in pairs)

layer(x::Symbol) = layer(Val(x))
layer(::Val{:first}) = rollrepeat(chain(Rx(0.0), Rz(0.0)))
layer(::Val{:last}) = rollrepeat(chain(Rz(0.0), Rx(0.0)))
layer(::Val{:mid}) = rollrepeat(chain(Rz(0.0), Rx(0.0), Rz(0.0)))

import Yao.Blocks: apply!, blocks
apply!(r::AbstractRegister, blk::QCBM) = apply!(r, blk.circuit)
blocks(blk::QCBM) = blocks(blk.circuit)

function print_block(io::IO, qcbm::QCBM)
    printstyled(io, "QCBM"; bold=true, color=:red)
end

function initialize!(qcbm::QCBM)
    params = 2pi * rand(nparameters(qcbm))
    dispatch!(qcbm, params)
end

function (x::QCBM{N})(nbatch::Int=1) where N
    apply!(zero_state(N, nbatch), x.circuit)
end

function gradient(qcbm::QCBM{N, NL}, kernel, ptrain) where {N, NL}
    prob = abs2.(statevec(qcbm()))
    grad = zeros(real(datatype(qcbm)), nparameters(qcbm))
    idx = 0
    for ilayer = 1:2:(2 * NL + 1)
        idx = grad_layer!(grad, idx, prob, qcbm, qcbm.circuit[ilayer], kernel, ptrain)
    end
    grad
end

function grad_layer!(grad, idx, prob, qcbm, layer, kernel, ptrain)
    count = idx
    for each_line in blocks(layer)
        for each in blocks(each_line)
            gradient!(grad, count+1, prob, qcbm, each, kernel, ptrain)
            count += 1
        end
    end
    count
end

function gradient!(grad, idx, prob, qcbm, gate, kernel, ptrain)
    dispatch!(+, gate, pi / 2)
    prob_pos = abs2.(statevec(qcbm()))

    dispatch!(-, gate, pi)
    prob_neg = abs2.(statevec(qcbm()))

    dispatch!(+, gate, pi / 2) # set back

    grad_pos = Kernels.expect(kernel, prob, prob_pos) - Kernels.expect(kernel, prob, prob_neg)
    grad_neg = Kernels.expect(kernel, ptrain, prob_pos) - Kernels.expect(kernel, ptrain, prob_neg)
    grad[idx] = grad_pos - grad_neg
    grad
end

qcbm = QCBM{4, 5}([1=>2, 2=>3, 3=>4, 4=>1])
nparameters(qcbm)
initialize!(qcbm)
dispatch!(-, qcbm, 1:64)

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
export gradtest, gradtest_assembleStiffAndForce

@doc raw"""
    gradtest(f::Function, x0::Array{Float64}; scale::Float64 = 1.0)

Testing the gradients of a vector function `f`.
Here `x0` is $n$-dimensional vector, `f` takes an $n$-dimensional vector as inputs, and
outputs a $m$ dimensional vector. The function tests the second order convergence of `f`
```math
\|f(x_0+\gamma c) - f(x_0) - \gamma \nabla f(x_0) c\|_2 = \mathcal{O}(\gamma^2)
```
"""
function gradtest(f::Function, x0::Array{Float64}; scale::Float64 = 1.0)
    v0 = rand(Float64,size(x0))
    γs = scale ./10 .^(1:5)
    err2 = []
    err1 = []
    f0, J = f(x0)
    for i = 1:5
        f1, _ = f(x0+γs[i]*v0)
        push!(err1, norm(f1-f0))
        # @show f1, f0, 2γs[i]*J*v0
        push!(err2, norm(f1-f0-γs[i]*J*v0))
        # push!(err2, norm((f1-f2)/(2γs[i])-J*v0))
        # #@show "test ", f1, f2, f1-f2
    end
    close("all")
    loglog(γs, err2, label="Automatic Differentiation")
    loglog(γs, err1, label="Finite Difference")
    loglog(γs, γs.^2 * 0.5*abs(err2[1])/γs[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
    loglog(γs, γs * 0.5*abs(err1[1])/γs[1], "--",label="\$\\mathcal{O}(\\gamma)\$")
    plt.gca().invert_xaxis()
    legend()
    println("Finite difference: $err1")
    println("Automatic differentiation: $err2")
    return err1, err2
end

function gradtest_assembleStiffAndForce(domain::Domain; scale::Float64 = 1.0)
    function f(s)
        domain.state[domain.dof_to_eq] = s
        rhs, K = assembleStiffAndForce(domain, 0.0)
        return rhs, K
    end
    gradtest(f, rand(domain.neqs), scale=scale)
end


text:
using Test, LinearAlgebra

using XPORTA: write_ieq, IEQ

using BellScenario

@testset "N-(N-4)-N scenario" begin

    @time for N = 6:6
        d = N - 4

        filepath = "./data/quick_adjacency_decomposition/"
        filename = "$N-$d-$N"

        println("$filename scenario")

        PM = LocalSignaling(N, N, d)

        M = zeros(Int64, N, N)
        M[[1, 2, 3], 1] .= 1

        for i = 1:(N-3)
            M[i+3, i] = 1
        end

        BG_seed = BellGame(M, d)

        vertices = LocalPolytope.vertices(PM)
        canonical_facets = LocalPolytope.adjacency_decomposition(
            vertices,
            BG_seed,
            PM,
            max_vertices = 36,
        )
        games = collect(keys(canonical_facets))

        pretty_print_txt(games, filepath * filename)
        write_ieq(filepath * "ieq/" * filename, convert(IEQ, games))
    end

end


text:
# Test real-to-complex FFT using a backwards plan (BRFFT).

using PencilFFTs
using MPI
using LinearAlgebra
using FFTW

using Test

MPI.Init()
comm = MPI.COMM_WORLD

let dev_null = @static Sys.iswindows() ? "nul" : "/dev/null"
    MPI.Comm_rank(comm) == 0 || redirect_stdout(open(dev_null, "w"))
end

@testset "BRFFT: odd = $odd" for odd ∈ (false, true)
    dims_real = (12, 13, 16 + odd)  # dimensions in physical (real) space
    dims_coef = (12, 13, 9)         # dimensions in coefficient (complex) space
    pen = Pencil(dims_coef, comm)
    uc = PencilArray{ComplexF64}(undef, pen)
    uc .= 0
    uc[2, 4, 3] = 1 + 2im

    plan_c2r = PencilFFTPlan(uc, Transforms.BRFFT(dims_real))
    @test size(plan_c2r) == dims_coef  # = size of input

    ur = plan_c2r * uc
    @test size_global(ur, LogicalOrder()) == dims_real

    # Equivalent using FFTW.
    # Note that by default FFTW performs the c2r transform along the first
    # dimension, while PencilFFTs does it along the last one.
    uc_fftw = gather(uc)
    ur_full = gather(ur)
    if uc_fftw !== nothing
        bfft!(uc_fftw, (1, 2))
        ur_fftw = brfft(uc_fftw, dims_real[end], 3)
        @test ur_full ≈ ur_fftw
    end

    # Check normalisation
    uc_back = plan_c2r \ ur
    @test isapprox(uc_back, uc; atol = 1e-8)
end


text:
using SimplexSplitting, InvariantDistribution
include("invariantize_embedding.jl")

function m1(t::Triangulation; n_pts = 100, sample_randomly = false)
	mm_dd(t, n_randpts = n_pts, sample_randomly = sample_randomly)
end

function m2(t::Triangulation; n_pts = 100, sample_randomly = false)
	mm_dd2(t, n_randpts = n_pts, sample_randomly = sample_randomly)
end

function m3(t::Triangulation; n_pts = 100, sample_randomly = false)
	dim = size(t.points, 2)
	convex_coeffs = get_coeffs(dim, n_pts, sample_randomly)
	simplices, imsimplices = prepare_mm_dd(t)
	mm_dd3(t, simplices, convex_coeffs)
end

function m4(t::Triangulation; n_pts = 100, sample_randomly = false)
	mm_dd4(t, n_randpts = n_pts, sample_randomly = sample_randomly)
end

@testset "Markov matrix approx: version comparison" begin
	using BenchmarkTools
	emb = invariant_gaussian_embedding(npts = 30)
	t = triang_from_embedding(Embedding(emb))

	@testset "Uniform point cloud"
		@time M1 = m1(t, n_randpts = 100, sample_randomly = false)
		@time M2 = m2(t, n_randpts = 100, sample_randomly = false)
		@time M3 = m2(t, n_randpts = 100, sample_randomly = false)
		@time M4 = m2(t, n_randpts = 100, sample_randomly = false)
		@time MEX = mm_p(t)

		@test all(sum(M1, 2) .≈ 1.0)
		@test all(sum(M2, 2) .≈ 1.0)
		@test all(sum(M3, 2) .≈ 1.0)
		@test all(sum(M4, 2) .≈ 1.0)
		@test all(sum(MEX, 2) .≈ 1.0)

		@show "Maximum discrepancy between exact and approx estimators:"
		@show maximum(abs(M4 .- MEX))
	end

	@testset "Random point cloud"
		@time M2 = mm_dd2(t, n_randpts = 100, sample_randomly = true)
		@time M4 = mm_dd4(t, n_randpts = 50, sample_randomly = true)
	end
end


text:
using DelayDiffEq
using Test

# Test parameterized delayed logistic equation

# delayed logistic equation
f_inplace(du, u, h, p, t) = (du[1] = p[1] * u[1] * (1 - h(p, t-1; idxs = 1)))
f_scalar(u, h, p, t) = p[1] * u * (1 - h(p, t-1))

# simple history function
h(p, t; idxs=nothing) = 0.1

@testset for inplace in (true, false)
  # define problem
  # we specify order_discontinuity_t0 = 1 to indicate that the discontinuity at
  # t = 0 is of first order
  prob = DDEProblem(inplace ? f_inplace : f_scalar,
                    inplace ? [0.1] : 0.1,
                    h, (0.0, 50.0), [0.3];
                    constant_lags = [1],
                    order_discontinuity_t0 = 1)

  # solve problem with initial parameter:
  sol1 = solve(prob, MethodOfSteps(Tsit5()))
  @test length(sol1) == 21
  @test first(sol1(12)) ≈ 0.884 atol = 1e-4
  @test first(sol1[end]) ≈ 1 atol = 1e-5

  # solve problem with updated parameter
  prob.p[1] = 1.4
  sol2 = solve(prob, MethodOfSteps(Tsit5()))
  @test length(sol2) == 47
  @test first(sol2(12)) ≈ 1.125 atol = 5e-4
  @test first(sol2[end]) ≈ 0.994 atol = 2e-5
end

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
#=
Geometry in 2d euclidean space
=#

"""
deltaangle(a::Real, b::Real)

Return the minimum δ such that
    a + δ = mod(b, 2π)
"""
deltaangle(a::Real, b::Real) = atan2(sin(b-a), cos(b-a))

# distance between two angles
angledist( a::Real, b::Real ) = abs(deltaangle(a,b))

# linear interpolation between angles
lerp_angle(a::Real, b::Real, t::AbstractFloat) = a + deltaangle(a, b)*t

# true if the values are collinear within a tolerance
function are_collinear(a::AbstractVec, b::AbstractVec, c::AbstractVec, tol::Float64=1e-8)
    # http://mathworld.wolfram.com/Collinear.html
    # if val = 0 then they are collinear
    val = a.x*(b.y-c.y) + b.x*(c.y-a.y)+c.x*(a.y-b.y)
    abs(val) < tol
end

function inertial2body(point::VecE2, reference::VecSE2)

    #=
    Convert a point in an inertial cartesian coordinate frame
    to be relative to a body's coordinate frame

    The body's position is given relative to the same inertial coordinate frame
    =#

    s, c = sin(reference.θ), cos(reference.θ)
    Δx = point.x - reference.x
    Δy = point.y - reference.y
    VecE2(c*Δx + s*Δy, c*Δy - s*Δx)
end
function inertial2body(point::VecSE2, reference::VecSE2)

    #=
    Convert a point in an inertial cartesian coordinate frame
    to be relative to a body's coordinate frame

    The body's position is given relative to the same inertial coordinate frame
    =#

    s, c = sin(reference.θ), cos(reference.θ)
    Δx = point.x - reference.x
    Δy = point.y - reference.y
    VecSE2(c*Δx + s*Δy, c*Δy - s*Δx, point.θ - reference.θ)
end
function body2inertial(point::VecE2, reference::VecSE2)

    #=
    Convert a point in a body-relative cartesian coordinate frame
    to be relative to a the inertial coordinate frame the body is described by
    =#

    c, s = cos(reference.θ), sin(reference.θ)
    VecE2(c*point.x -s*point.y + reference.x, s*point.x +c*point.y + reference.y)
end
function body2inertial(point::VecSE2, reference::VecSE2)

    #=
    Convert a point in a body-relative cartesian coordinate frame
    to be relative to a the inertial coordinate frame the body is described by
    =#

    c, s = cos(reference.θ), sin(reference.θ)
    VecSE2(c*point.x -s*point.y + reference.x, s*point.x +c*point.y + reference.y, reference.θ + point.θ)
end

include("1d.jl")
include("lines.jl")
include("line_segments.jl")
include("rays.jl")
include("projectiles.jl")
include("solids.jl")


text:
module RunningVectors
abstract AbstractRunningVector

type RunningVector{T} <: AbstractRunningVector
	d::Vector{T}
	offset::Int
end
type RunningBitVector <: AbstractRunningVector
	d::BitVector
	offset::Int

end
abstract AbstractRunningSumVector <: AbstractRunningVector
type RunningSumVector{T} <: AbstractRunningSumVector
	d::Vector{T}
	offset::Int
	sum::Int
end
type RunningSumBitVector <: AbstractRunningSumVector
	d::BitVector
	offset::Int
	sum::Int
end
RunningSumBitVector() = RunningSumBitVector(BitVector(0),0,0)
RunningVector{T}(::Type{T}) = RunningVector(Array(T,0),0)
RunningBitVector() = RunningBitVector(BitVector(0),0)

Base.getindex(r::AbstractRunningVector, inds::Colon) = getindex(r.d, 1:endof(r))
Base.getindex(r::AbstractRunningVector, inds) = getindex(r.d, inds-r.offset)
Base.setindex!(r::AbstractRunningVector, v, inds) = setindex!(r.d, v, inds-r.offset)
Base.length(r::AbstractRunningVector) = r.offset + length(r.d)
Base.endof(r::AbstractRunningVector) = length(r)
Base.size(r::AbstractRunningVector) = (length(r),)
Base.resize!(r::AbstractRunningVector,n) = resize!(r.d, n-r.offset)
Base.eltype(r::AbstractRunningVector) = eltype(r.d)
runningendof(x) = typemax(Int)
runningendof(r::AbstractRunningVector) = endof(AbstractRunningVector)
available(r::AbstractRunningVector) = r.offset+1:length(r)

function freeuntil!(r::AbstractRunningVector, i::Int)
	r.offset >= i && return r
	todelete = min(length(r.d),i-r.offset)
	deleteat!(r.d, 1:todelete)
	r.offset = i
	r
end

Base.push!(r::AbstractRunningVector, items...) = push!(r.d, items...)
Base.append!(r::AbstractRunningVector, items) = append!(r.d, items)
function Base.append!(r::AbstractRunningSumVector, items)
	r.sum+=sum(items)
	append!(r.d,items)
end
function Base.push!(r::AbstractRunningSumVector, items...)
	r.sum+=sum(items)
	push!(r.d,items...)
end
function Base.setindex!(r::AbstractRunningSumVector, v, inds)
	oldv = getindex(r.d, inds-r.offset)
	r.sum-=sum(oldv)
	r.sum+=sum(v)
	setindex!(r.d,v, inds-r.offset)
end
Base.sum(r::AbstractRunningSumVector) = r.sum

function Base.show(io::IO, r::AbstractRunningVector)
	print(io, string(typeof(r)))
	print(io, " length=$(length(r)), offset=$(r.offset)")
end

export RunningVector, RunningBitVector, RunningSumVector, RunningSumBitVector, freeuntil!, runningendof, available,
		AbstractRunningVector, AbstractRunningSumVector
end


text:
function calculate_losses!(result, data)
    for (l, branch) in result["solution"]["branch"]
        branch["ploss"] = branch["pf"] + branch["pt"]
    end
    result["totalloss"] = sum(branch["ploss"] for (l, branch) in result["solution"]["branch"])

    result["totalload"] = sum(load["pd"] for (l, load) in data["load"])

    result["totalgen"] = sum(gen["pg"] for (l, gen) in result["solution"]["gen"])

    return result["totalloss"]
end

"This function accepts a data dictionary, and returns a new data dictionary
with perturbation applied"
function create_impedance_perturbation(data_input, α, ϵ, λ)
    # Create a copy of the input so that a new dataset is returned
    data = deepcopy(data_input)

    # First apply the Laplace noise to each branch
    distribution = Distributions.Laplace(0, 3 * α / ϵ)
    sum_g = 0
    sum_b = 0
    sum_g_shunt = 0
    sum_b_shunt = 0
    for (l, branch) in data["branch"]
        noise = Random.rand(distribution, 1)[1]
        noise_shunt = Random.rand(distribution, 1)[1]
        z = branch["br_r"] + im*branch["br_x"]
        y = 1 / z
        g = real(y)
        b = imag(y)
        g_shunt = branch["g_to"] #assumes it is equal to branch["g_fr"]
        b_shunt = branch["b_to"] #assumes it is equal to branch["b_fr"]
        # println("g is: ", g)
        # println("b is: ", b)

        # Handle the case where resistance is 0 to avoid undefined behaviour
        if branch["br_r"] != 0
            r = b / g # Algorithm 1 and eq. (15)
            pert_b = b + noise
            pert_g = pert_b / r
        else
            pert_g = 0
            pert_b = b + noise # Just apply perturbation to b
        end

        pert_b_shunt = b_shunt + noise_shunt
        pert_g_shunt = g_shunt # Removed the perturbation of g_shunt

        branch["g_obj"] = pert_g
        branch["b_obj"] = pert_b
        branch["g_shunt_obj"] = pert_g_shunt
        branch["b_shunt_obj"] = pert_b_shunt
        # println("perturbed g is: ", pert_g)
        # println("perturbed b is: ", pert_b)
        sum_g += g
        sum_b += b
        sum_g_shunt += g_shunt
        sum_b_shunt += b_shunt
    end

    # Apply the noisy mean limit values for s4 and s5
    n = length(data["branch"])
    lap_μ = Distributions.Laplace(0, 3 * α / (n * ϵ))

    μ_g = (1 / n) * sum_g + Random.rand(lap_μ, 1)[1] # Eq 16
    μ_b = (1 / n) * sum_b + Random.rand(lap_μ, 1)[1] # Eq 17
    μ_g_shunt = (1 / n) * sum_g_shunt + Random.rand(distribution, 1)[1]
    μ_b_shunt = (1 / n) * sum_b_shunt + Random.rand(distribution, 1)[1]
    # In the case that random perturbations create negative values, floor these at 0
    if μ_g_shunt < 0
        μ_g_shunt = 0
    end
    if μ_b_shunt < 0
        μ_b_shunt = 0
    end

    # Add noisy mean limit values to our data dictionary
    # Note that we use min and max to handle negative parameters
    data["g_lb"] = min(μ_g / λ, μ_g * λ)
    data["g_ub"] = max(μ_g / λ, μ_g * λ)
    data["b_lb"] = min(μ_b / λ, μ_b * λ)
    data["b_ub"] = max(μ_b / λ, μ_b * λ)
    data["g_lb_shunt"] = min(μ_g_shunt / λ, μ_g_shunt * λ)
    data["g_ub_shunt"] = max(μ_g_shunt / λ, μ_g_shunt * λ)
    data["b_lb_shunt"] = min(μ_b_shunt / λ, μ_b_shunt * λ)
    data["b_ub_shunt"] = max(μ_b_shunt / λ, μ_b_shunt * λ)
    @show data["b_lb_shunt"], data["b_ub_shunt"]

    return data
end

function overwrite_impedances_in_data!(result, data)
    for (l, branch) in result["solution"]["branch"]
        y = branch["g"] + im* branch["b"]
        z = 1/y
        r = real(z)
        x = imag(z)
        b = branch["b_shunt"]
        data["branch"][l]["br_r"] = r
        data["branch"][l]["br_x"] = x
        data["branch"][l]["b_fr"] = b
        data["branch"][l]["b_to"] = b
    end
end