text:
using GenerativeAD
using FileIO, BSON
using ValueHistories, DistributionsAD
using Flux
using ConditionalDists
using GenerativeModels
using EvalMetrics
using Plots

f = "vae/iris/seed=3/model_activation=tanh_batchsize=32_hdim=256_init_seed=59386742_lr=0.001_nlayers=4_zdim=2.bson"
#######
function get_crit_vae(f)
	model_data = load(f)
	loss = get(model_data["history"], :validation_likelihood)[2][end]
end
#######

#path = "vae/iris/seed=3"
path = "vae/htru2/seed=3"
mfs = filter(x->occursin("model", x), readdir(path, join=true))

crit_vals = map(get_crit_vae, mfs)

seeds = map(x->split(split(x, "init_seed")[2], "_")[1][2:end], mfs)
sfs = filter(x->occursin("reconstruction-sampled", x), readdir(path, join=true))

function get_auc(sfs, seed)
	f = sfs[findfirst(x->occursin(seed,x), sfs)]
	data = load(f)
	auc = EvalMetrics.auc_trapezoidal(EvalMetrics.roccurve(data[:val_labels], vec(data[:val_scores]))...)
end

aucs = map(x->get_auc(sfs, x), seeds)

scatter(crit_vals, aucs, xlabel="validation ELBO", ylabel="validation AUC", title="VAE - iris")
savefig("~/generativead/GenerativeAD.jl/plots/elbo_vs_auc.png")

function get_auc_test(sfs, seed)
	f = sfs[findfirst(x->occursin(seed,x), sfs)]
	data = load(f)
	auc = EvalMetrics.auc_trapezoidal(EvalMetrics.roccurve(data[:tst_labels], vec(data[:tst_scores]))...)
end
aucs = map(x->get_auc_test(sfs, x), seeds)
scatter(crit_vals, aucs, xlabel="validation ELBO", ylabel="test AUC", title="VAE - iris")
savefig("~/generativead/GenerativeAD.jl/plots/elbo_vs_auc_test.png")

sfs = filter(x->occursin("reconstruction-sampled", x), readdir(path, join=true))

function return_rec_auc(f)
	data = load(f)
	auc = EvalMetrics.auc_trapezoidal(EvalMetrics.roccurve(data[:val_labels], vec(data[:val_scores]))...)
	rec = Flux.mean(data[:val_scores][data[:val_labels].==0])
	auc, rec
end

res = map(return_rec_auc, sfs)
aucs = [x[1] for x in res]
recs = [x[2] for x in res]
scatter(recs, aucs, xlabel="validation reconstruction error", ylabel="validation AUC", title="VAE - iris")
savefig("~/generativead/GenerativeAD.jl/plots/rec_sampled_vs_auc.png")

## WAE
############
f = "wae/iris/seed=3/model_activation=tanh_batchsize=64_hdim=512_init_seed=86274140_kernel=gauss_lambda=0.1_lr=0.0001_nlayers=4_sigma=0.001_zdim=2.bson"
model_data = load(f)

#############
function get_crit_wae(f)
	model_data = load(f)
	loss = get(model_data["history"], :validation_likelihood)[2][end]
end
function get_crit_aae(f)
	model_data = load(f)
	loss = get(model_data["history"], :validation_loss)[2][end]
end
#############

path = "wae/iris/seed=3"
mfs = filter(x->occursin("model", x), readdir(path, join=true))
crit_vals = map(get_crit_wae, mfs)
aucs = map(x->get_auc(sfs, x), seeds)
scatter(crit_vals, aucs, xlabel="validation loss", ylabel="validation AUC", title="WAE - iris")

seeds = map(x->split(split(x, "init_seed")[2], "_")[1][2:end], mfs)
sfs = filter(x->occursin("reconstruction-sampled", x), readdir(path, join=true))
res = map(return_rec_auc, sfs)
aucs = [x[1] for x in res]
recs = [x[2] for x in res]
scatter(recs, aucs, xlabel="validation reconstruction error", ylabel="validation AUC", title="WAE - iris")

path = "aae/iris/seed=3"
mfs = filter(x->occursin("model", x), readdir(path, join=true))
crit_vals = map(get_crit_aae, mfs)
aucs = map(x->get_auc(sfs, x), seeds)
scatter(crit_vals, aucs, xlabel="validation loss", ylabel="validation AUC", title="AAE - iris")

seeds = map(x->split(split(x, "init_seed")[2], "_")[1][2:end], mfs)
sfs = filter(x->occursin("reconstruction-sampled", x), readdir(path, join=true))
res = map(return_rec_auc, sfs)
aucs = [x[1] for x in res]
recs = [x[2] for x in res]
scatter(recs, aucs, xlabel="validation reconstruction error", ylabel="validation AUC", title="AAE - iris")

#### FMGAN


text:
using Test
using ReversibleSeismic
using NiLang
using NiLang.AD
using CUDA

CUDA.allowscalar(false)

@testset "identity" begin
    x = CUDA.CuArray([1.0,2.0,3.0])
    y = CUDA.CuArray([5.0,2.0,3.0])
    @instr x += y
    @test Array(x) == [6.0, 4.0, 6.0]
end

@i function i_loss_bennett_gpu!(out!::T, state, param, src, srcv::Vector{T}, c::AbstractMatrix{T}; bennett_k=50, logger=NiLang.BennettLog()) where T
    bennett!((@const bennett_step!), state, bennett_k, 1, (@const param.NSTEP-1), param, src, srcv, c; do_uncomputing=false, logger=logger)
    out! += sum((@skip! abs2), state[param.NSTEP].u)
end

function loss_bennett_gpu(c::AbstractMatrix{T}; nstep, usecuda=true) where T
    nx = size(c, 1) - 2
    ny = size(c, 2) - 2
    param = AcousticPropagatorParams(nx=nx, ny=ny,
         Rcoef=0.2, dx=20.0, dy=20.0, dt=0.05, nstep=nstep)
    src = (nx ÷ 2, ny ÷ 2)
    srcv = Ricker(param, 100.0, 500.0)

    if usecuda
        state = Dict(1=>CuSeismicState(Float64, nx, ny))
        param = togpu(param)
        c = CuArray(c)
    else
        state = Dict(1=>SeismicState(Float64, nx, ny))
        c = copy(c)
    end
    i_loss_bennett_gpu!(0.0, state, param, src, srcv, c)[1]
end

@testset "loss" begin
    T = Float64
    nx = ny = 99
    na = 52
    nb = 41
    param0 = AcousticPropagatorParams(nx=nx, ny=ny,
         Rcoef=0.2, dx=20.0, dy=20.0, dt=0.05, nstep=(na-2)*(nb-1)+2)

    tua = CUDA.zeros(T, nx+2, ny+2, na)
    tφa = CUDA.zeros(T, nx+2, ny+2, na)
    tψa = CUDA.zeros(T, nx+2, ny+2, na)
    tub = CUDA.zeros(T, nx+2, ny+2, 2nb)
    tφb = CUDA.zeros(T, nx+2, ny+2, nb)
    tψb = CUDA.zeros(T, nx+2, ny+2, nb)

    c0 = 1000*CUDA.ones(T, nx+2, ny+2)
    src = (nx ÷ 2, ny ÷ 2)
    srcv0 = Ricker(param0, 100.0, 500.0)
    u = solve(param0, src, copy(srcv0), Array(c0))[:,:,end]
    ls = sum(abs2, u)

    param0 = AcousticPropagatorParams(nx=nx, ny=ny,
         Rcoef=0.2, dx=20.0, dy=20.0, dt=0.05, nstep=(na-2)*(nb-1)+2)
    s0 = CuSeismicState(Float64, nx, ny)
    state = Dict(1=>s0)
    loss_bg = i_loss_bennett_gpu!(0.0, copy(state), togpu(param0), src, copy(srcv0), copy(c0); bennett_k=200)[1]

    s0 = CuSeismicState(Float64, nx, ny)
    state = Dict(1=>s0)
    loss_cpu = i_loss_bennett_gpu!(0.0, Dict(1=>SeismicState(Float64, nx, ny)), deepcopy(param0), src, copy(srcv0), Array(c0))[1]
    @test loss_bg ≈ loss_cpu
    @test loss_bg ≈ ls
end

function getgrad_bennett_gpu(c::AbstractMatrix{T}; nstep) where T
    nx, ny = size(c) .- 2
    param = AcousticPropagatorParams(nx=size(c,1)-2, ny=size(c,2)-2,
          Rcoef=0.2, dx=20.0, dy=20.0, dt=0.05, nstep=nstep) |> togpu
    c = c |> CuArray
    src = size(c) .÷ 2 .- 1
    srcv = Ricker(param, 100.0, 500.0)
    s0 = CuSeismicState(Float64, nx, ny)
    state = Dict(1=>s0)
    NiLang.AD.gradient(Val(1), i_loss_bennett_gpu!, (0.0, state, param, src, srcv, c))[6]
end

function getgrad_treeverse_gpu(c::AbstractMatrix{T}; nstep) where T
    nx, ny = size(c) .- 2
    param = AcousticPropagatorParams(nx=size(c,1)-2, ny=size(c,2)-2,
          Rcoef=0.2, dx=20.0, dy=20.0, dt=0.05, nstep=nstep) |> togpu
    c = c |> CuArray
    src = size(c) .÷ 2 .- 1
    srcv = Ricker(param, 100.0, 500.0)
    s0 = CuSeismicState(Float64, nx, ny)
    function gn(sn)
        g = zero(sn)
        g.u .+= 2 .* sn.u
        return (g, zero(srcv), zero(c))
    end
    treeverse_solve(s0, gn; param=param, src=src, srcv=srcv, c=c)[2][3]
end

"""
obtain gradients numerically, for gradient checking.
"""
function getngrad_gpu(c0::AbstractMatrix{T}, i, j; nstep::Int, δ=1e-4) where T
     c = c0 |> CuArray
     c[ReversibleSeismic.SafeIndex(i,j)] += δ
     fpos = loss_bennett_gpu(c; nstep)
     c = c0 |> CuArray
     c[ReversibleSeismic.SafeIndex(i,j)] -= δ
     fneg = loss_bennett_gpu(c; nstep)
     return (fpos - fneg)/2δ
end

@testset "gradient" begin
    nx = ny = 99
    na = 42
    nb = 51
    nstep=(na-2)*(nb-1)+2
    c = 1000*ones(Float64, nx+2, ny+2)
    g2 = Array(getgrad_bennett_gpu(c; nstep=nstep))
    g3 = Array(getgrad_treeverse_gpu(c; nstep=nstep))
    ng4545 = getngrad_gpu(c, 45, 45; nstep=nstep, δ=1e-3)
    @test isapprox(g2[45,45], ng4545; rtol=1e-2)
    @test isapprox(g3[45,45], ng4545; rtol=1e-2)
end

ext: jl


text:
using skalman
using Plots, LinearAlgebra, FastGaussQuadrature

n = 8
nodes, weights = gausslobatto(n)

# Broadcast and store in a matrix
P = legendrepoly.(nodes,n-1)
P = reduce(hcat,P)
P = transpose(P)

# Broadcast and store in a matrix
PX = derlegendre.(nodes,n-1)
PX = reduce(hcat,PX)
PX = transpose(PX)

M = zeros(n,n)
S = similar(M)
for i = 1:n
    for j = 1:n
    M[i,j] = sum(P[:,i].*P[:,j].*weights)
    S[i,j] = sum(P[:,i].*PX[:,j].*weights)
    end
end

f = sin.(nodes)
b = similar(f)
for i = 1:n
    b[i] = sum(P[:,i].*f.*weights)
end
umodes = M\b
uxmodes = M\(S*umodes)

zz = collect(LinRange(-1.0, 1.0, 101))
# Broadcast and store in a matrix
PP = legendrepoly.(zz,n-1)
PP = reduce(hcat,PP)
PP = transpose(PP)

unodes = PP*umodes
uxnodes = PP*uxmodes
pl = plot(zz,uxnodes)
plot!(pl,zz,cos.(zz))


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

using Plots;

"""
    nrbplot(nurbs::NURBS{I,F},
            npts::Union{I, Vector{I}} [;
            plotargs... ]
           ) where {I<:Integer,F<:AbstractFloat}

Plot a NURBS curve or a surface.

# Arguments:
- `nurbs`: NURBS structure
- `npts`: Number of evaluation points, for a surface a vector with two elements
  for the number of points along the U and V directions respectively.
- `plotargs`: plot arguments to pass to the wrapped plot function

!!! note "Note:"
    For 3D Plots `nrbplot` uses the `Plots.pyplot()` backend. Use plot
    arguments compatible with that backend.

# Examples:
```julia
julia> nrbplot(nrb,subd)
julia> nrbplot(nrb,subd,c=:winter)
```

Plot the test surface with 20 points along the U direction and 30 along the V
direction
```julia
julia> nrbplot(nrbtestsrf,[20;30],c=:copper)
```

!!! note "Note:"
    Contrary to the Matlab version does this ported version supports no
    specialised input arguments just pass named plot arguments which will
    be passed to the wrapped plot functions.

!!! note "Note:"
    It is also possible to overload the Julia `plot()` function with the
    mechanic [Recipes](https://docs.juliaplots.org/latest/recipes).

# Reference
- see also: [`nrbeval`](@ref)
"""
function nrbplot(nurbs::NURBS{I,F},
                 npts::Union{I, Vector{I}};
                 plotargs...
                ) where {I<:Integer,F<:AbstractFloat}
# TODO: this ported version does not do all the tests of the matlab version,
#       instead the type features of julia are used -> check if the matlab
#       errors get caught
# TODO: have a look into Recipes to make the plotting more adaptive and to
#       maybe replace this function

# TODO: make this workaround obsolete with Recipes
# filter out key indicating append to old plot as workaround
# in order not to write a near indetical copy of this function
appendflag = false;
if :nrbplot_append in keys(plotargs)
  appendflag = plotargs[:nrbplot_append];
  plotargs = (; filter(p->first(p) != :nrbplot_append,
                       collect(pairs(plotargs)))...);
end

subd = npts .+ 1;

if typeof(nurbs)<:NURBS2D
  if !(typeof(npts)<:Vector{I}) || length(npts)!=2
    throw(ArgumentError("To plot surfaces a Vector{I} with the length 2 is " *
                        "needed!"));
  end # if
  # default color map
  defaults = (; c=:jet, st=:surface, aspect_ratio=:auto);
  settings = merge( defaults, plotargs );
  # TODO: pyplot 3D does not supports equal axis scale, add that later

  # plot a NURBS surface
  p = nrbeval(nurbs,[collect(range(0,1,length=subd[1])),
                     collect(range(0,1,length=subd[2]))]);

  Plots.pyplot();
  if appendflag
   plot!(p[1,:,:],p[2,:,:],p[3,:,:];settings...);
 else
   plot(p[1,:,:],p[2,:,:],p[3,:,:];settings...);
 end
  # TODO: the julia plotting eco system for 3D plots is still a bit chaotic
  #       come back later and make this solution a bit more robust
else
  if !(typeof(npts)<:I)
    throw(ArgumentError("Valid argument of npts for curves is <:Integer!"));
  end # if

  # plot a NURBS curve
  p = nrbeval(nurbs,collect(range(0,1,length=subd)));

  if all(i->i==0,p[3,:])
    # 2D curve

    # default color map
    defaults = (; aspect_ratio=:equal);
    settings = merge( defaults, plotargs );

    if appendflag
      plot!(p[1,:],p[2,:];settings...);
    else
      plot(p[1,:],p[2,:];settings...);
    end
  else
    # 3D curve

    # default color map
    defaults = (; aspect_ratio=:auto);
    settings = merge( defaults, plotargs );
    # TODO: pyplot 3D does not supports equal axis scale, add that later

    Plots.pyplot();
    if appendflag
      plot!(p[1,:],p[2,:],p[3,:];settings...);
    else
      plot(p[1,:],p[2,:],p[3,:];settings...);
    end
    # TODO: the julia plotting eco system for 3D plots is still a bit chaotic
    #       come back later and make this solution a bit more robust
  end
end # if
end # nrbplot

"""
    nrbplot!(nurbs::NURBS{I,F},
             npts::Union{I, Vector{I}} [;
             plotargs... ]
            ) where {I<:Integer,F<:AbstractFloat}

Plot a NURBS curve or a surface.

Wrapper of [`nrbplot`](@ref) in order to introduce the Julia plot append
mechanic for that plot routine.

!!! note "Note:"
    This function serves only as wrapper for the [`nrbplot`](@ref) function for
    input arguments and usage look at that function.

# Reference
- see [`nrbplot`](@ref)
"""
function nrbplot!(nurbs::NURBS{I,F},
                  npts::Union{I, Vector{I}};
                  plotargs...
                 ) where {I<:Integer,F<:AbstractFloat}
# append nrbplot_append = true to plotargs
defaults = (; nrbplot_append=true);
plotargs = merge( defaults, plotargs );
nplt = nrbplot(nurbs,npts; plotargs...);
end


text:
using Revise, Dates, DataFrames, GLM
#NSymbol = Union{Nothing, Symbol}
#NInt = Union{Nothing, Int} #NOTE: consider deleting this later
import Base: +, -, ==, >, <, ≥, ≤, length, isless, isequal
import StatsModels: implicit_intercept

#include("Finometrics.jl")
using Revise
const FMExpr = Union{Symbol,Expr,Nothing}
#  StatsModels.implicit_intercept(::Type{<:Any}) = true
function getModelMatrix(df::T, f::FormulaTerm)::Matrix{Float64} where
  T <: AbstractDataFrame

  return ModelMatrix(ModelFrame(f, df, model=LinearModel)).m
end

#Same as above but allows for an expression
function getModelMatrix(df::T, exp::V)::Matrix{Float64} where
  {T <: AbstractDataFrame, V <: FMExpr}

  #special case which crashes Formula
  if exp == Symbol("")
    return ones(Float64,size(df,1),1)
  end

  return getModelMatrix(df, get1SidedFormula(exp))
end

#helper function to create a one-sided formula given an expression
#IN: an expression and dataframe
#OUT: A one-sided formula object
function get1SidedFormula(RHS::T)::FormulaTerm where
  {T <: FMExpr}

  return @eval(@formula(identity(1) ~ $RHS)) #WARNING: THIS IS A HACK!!!!
end

function testMM()
  df = DataFrame(x = 1001:2000, y = rand(1000), z=rand(1000),
  f = (i->Symbol(:q,i÷100)).(1:1000))
  #println(df)
  x = getModelMatrix(df, :(x+y))

  println(x)
end

testMM()
#=using Finometrics
using Revise

#use this to test the lag functions
function testlaganddifference(N = 10_000)
  rowindex = collect(1:N)
  periodcol = deepcopy(rowindex)

  #make a test frame
  df = DataFrame(rowindex=rowindex, period = periodcol, group1 = periodcol .÷ 17,
    group2 = periodcol .÷ 29, target1 = rand(N), target2 = rand(N))

  #sort!(df, [:group1, :period, :group2])
  #CASE 1: create lagged field Lgroup1
  lagwithin!(df, :target1, :group1, :period)

  #CASE 2: create lagged fields L2target1 and L2target2
  lagwithin!(df, [:target1, :target2], [:group1, :group2], :period, lags=2, sorted=false)

  #CASE 3: difference :Ltarget1 and create DLtarget1 and LLtarget1
  differencewithin!(df, :Ltarget1, :group1, :period, deletelag=false, sorted=true,
    laggedname=:LLtarget1e, differencedname=:DLtarget1e)

  #CASE 4: difference [:L2target1, :L2target2] creating
  # [:DL2target1e, :DL2target2e], [:LL2target1e, :LL2target2e]
  differencewithin!(df, [:L2target1, :L2target2], [:group1, :group2], :period,
    deletelag=false, sorted=true, differencednames=[:DL2target1e, :DL2target2e],
    laggednames = [:LL2target1e, :LL2target2e])

  #validation for cases 1 and 3
  sort!(df, [:group1, :period])
  for subdf ∈ groupby(df, :group1)
    sN::Int = size(subdf, 1)
    for r ∈ 1:sN
      if r > 1 #CASE 1check
        @assert subdf[r-1, :target1] == subdf[r, :Ltarget1]
      else
        @assert ismissing(subdf[r, :Ltarget1])
      end

      if r > 2 #CASE 3 check
        @assert subdf[r-1, :Ltarget1] == subdf[r, :LLtarget1e]
        @assert subdf[r, :Ltarget1] - subdf[r, :LLtarget1e] == subdf[r, :DLtarget1e]
      else
        @assert ismissing(subdf[r, :LLtarget1e])
        @assert ismissing(subdf[r, :DLtarget1e])

      end
    end
  end

  #validation for cases 2 and 4
  sort!(df, [:group1, :group2, :period])
  for subdf ∈ groupby(df, [:group1, :group2])
    sN::Int = size(subdf, 1)
    for r ∈ 1:sN
      if r > 2 #CASE 2 check
        @assert subdf[r-2, :target1] == subdf[r, :L2target1]
        @assert subdf[r-2, :target2] == subdf[r, :L2target2]
      else
        @assert ismissing(subdf[r, :L2target1])
        @assert ismissing(subdf[r, :L2target2])
      end

      if r > 3 #CASE 4 check
        @assert subdf[r-1, :L2target1] == subdf[r, :LL2target1e]
        @assert subdf[r-1, :L2target2] == subdf[r, :LL2target2e]
        @assert subdf[r, :L2target1] - subdf[r, :LL2target1e] == subdf[r, :DL2target1e]
        @assert subdf[r, :L2target2] - subdf[r, :LL2target2e] == subdf[r, :DL2target2e]
      else
        @assert ismissing(subdf[r, :LL2target1e])
        @assert ismissing(subdf[r, :LL2target2e])
        @assert ismissing(subdf[r, :DL2target1e])
        @assert ismissing(subdf[r, :DL2target2e])

      end
    end
  end

  display(df)
end

#testlaganddifference()

function testYearQuarter()

  #test the first constructor
  yq1::YearQuarter = YearQuarter(1998,2)
  try
    YearQuarter(1985,7)
    @assert false
  catch err
    (typeof(err)<:AssertionError) && error("Error catching invalid values in first constructor")
  end

  #test the second constructor
  yq2::YearQuarter = YearQuarter(2000.4)
  try
    YearQuarter(2000.05)
    @assert false
  catch err
    (typeof(err)<:AssertionError) && error("Error catching invalid values in 2nd constructor")
  end

  #test addition of quarters
  for i ∈ 0:10
    @assert (yq1 + i) == (yq2 - (10-i))
  end
  #@assert yq1 + 10 == yq2
  #@assert yq2 - 10 == yq1

  @assert yq2 - yq1 == 10
  @assert yq1 < yq2
  @assert yq2 > yq1
  @assert yq2 ≥ yq1
  @assert yq1 ≥ yq1

  #println(yq1)
  #println(Float64(yq1))
end=#

abstract type YearPeriod end
##################################Year-Month#####################
struct YearMonth <: YearPeriod
  y::UInt16
  m::UInt8

  function YearMonth(y,m)
    (m≥1 && m≤12) || error("Invalid YearMonth $y $m")
    return new(y,m)
  end
end

#creates a year month from yyyymm date style
function YearMonth(i::Integer)::YearMonth
  local y::UInt16
  local m::UInt8

  (i>1000_00 && i<10000_00) || error("Invalid YearMonth $i")
  (y,m) = (i ÷ 100, i % 100)

  return YearMonth(y,m)
end

#fallback in case we need these
Int(ym::YearMonth)::Int = ym.y*100 + ym.m
Base.show(io::IO, ym::YearMonth) = print(io, "$(ym.y)_$(ym.m)")

#quick way to split the year quarter
Tuple{Int,Int}(ym::YearMonth)::Tuple{Int,Int} = (ym.y, ym.m)

#allows for adding quarters
function +(ym::YearMonth, madded::Integer)::YearMonth
  (y₀::Int, m₀::Int) = Tuple{Int,Int}(ym)

  m::Int = (m₀ + madded) % 12
  (m≤0) && (m+=12) #only allow positive months

  y::Int = (madded - (m-m₀)) ÷ 12 + y₀

  @assert m-m₀ + 12*(y-y₀) == madded #this should never fail

  return YearMonth(y,m)
end

#helper methods to build on the above
+(qadded::Int, yq::YearMonth)::YearMonth = ym + qadded
-(ym::YearMonth, madded::Int)::YearMonth = ym + (-1*madded)

#gets the number of months difference
function -(ym1::YearMonth, ym2::YearMonth)::Int

  Δy::Int = ym1.y - ym2.y
  Δm::Int = ym1.m - ym2.m

  return Δy * 12 + Δm
end

#comparison operators (primary sorting functions, then operator symbols)
isless(ym1::YearMonth, ym2::YearMonth) = ym1.y*100+ym1.m < ym2.y*100 + ym2.m#(ym1.y<ym2.y) || ((ym1.y==ym2.y) && (ym1.m<ym2.m))

#<(ym1::YearMonth, ym2::YearMonth) = isless(ym1, ym2)
#>(ym1::YearMonth, ym2::YearMonth) = isless(ym2, ym1)

#≤(ym1::YearMonth, ym2::YearMonth) = isless(ym1, ym2) || isequal(ym1, ym2)
#≥(ym1::YearMonth, ym2::YearMonth) = isless(ym2, ym1) || isequal(ym1, ym2)

#convert back to a date format
bom(ym::YearMonth)::Date = Date(ym.y, ym.m, 1)
eom(ym::YearMonth)::Date = lastdayofmonth(boq(ym))

Base.length(::YearPeriod) = 1
Base.iterate(yp::YearPeriod) = yp
Base.iterate(ym::YearPeriod, s::Integer) = nothing
Base.Broadcast.broadcastable(yp::YearPeriod) = Ref(yp)

function testYearMonth()

  #test the first constructor
  ym1::YearMonth = YearMonth(1998,2)
  try
    YearMonth(1985,14)
    @assert false
  catch err
    (typeof(err)<:AssertionError) && error("Error catching invalid values in first constructor")
  end

  #test the second constructor
  ym2::YearMonth = YearMonth(2008_04)
  try
    YearMonth(2000_14)
    @assert false
  catch err
    (typeof(err)<:AssertionError) && error("Error catching invalid values in 2nd constructor")
  end

  #test addition of quarters
  for i ∈ 0:122
    @assert (ym1 + i) == (ym2 - (122-i))
  end

  @assert ym2 - ym1 == 122
  @assert ym1 < ym2
  @assert ym2 > ym1
  @assert ym2 ≥ ym1
  @assert ym1 ≥ ym1

end

#@time for i ∈ 1:200_000
#  testYearMonth()
#end

function testymperformance(N::Int)
  ms::Vector{Int} = rand(1:12, N)
  ys::Vector{Int} = rand(0:20, N) .+ 1990
  addedmonths::Vector{Int} = rand(-120:120, N)
  local tot::Int
  local yms::Vector{YearMonth}
  local dts::Vector{Date}

  local addedMonths::Vector{Month} = (Month).(addedmonths)

  benchym::YearMonth = YearMonth(2000,1)
  benchdt::Date = Date(2000,1,1)

  @time begin
    yms = ((y::Int,m::Int)->YearMonth(y,m)).(ys,ms)
    yms .= yms .+ addedmonths
    tot = sum(yms .> benchym)
    print("\nYearMonth: $tot values found in")
  end

  @time begin
    dts = ((y::Int,m::Int)->Date(y,m,1)).(ys,ms)
    dts .= dts .+ addedMonths
    tot = sum(dts .> benchdt)
    print("Date: $tot values found in")
  end

end

#testymperformance(100_000_000)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# Write a function that concatenates two lists. [a,b,c], [1,2,3] → [a,b,c,1,2,3]

list1 = ["a","b","c"]
list2 = [1,2,3]

function concat(x::Vector{String},y::Vector{Int64})
    vcat(x,y)
end

concat(list1,list2)

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _lsf_rotate{T<:AbstractFloat}(Δ::T, vsini::T, ɛ::T)
    n_2 = floor(ceil(2*vsini/Δ)/2)
    vel = collect(-n_2:n_2)*Δ
    e1 = 2*(1.0 - ɛ)
    e2 = pi*ɛ/2.0
    e3 = pi*(1.0 - ɛ/3.0)
    x = vel/vsini
    x1 = abs(1.0 - x.^2)
    return vel, (e1*sqrt(x1) + e2*x1)/e3
end

"""
    lsf_rotate(delta_v, v_sin_i[, epsilon = 0.3]) -> velocity_grid, lsf

### Purpose ###

Create a 1-d convolution kernel to broaden a spectrum from a rotating star.

### Explanation ###

Can be used to derive the broadening effect (LSF, line spread function) due to
rotation on a synthetic stellar spectrum.  Assumes constant limb darkening
across the disk.

### Arguments ###

* `delta_v`: numeric scalar giving the step increment (in km/s) in the output
  rotation kernel
* `v_sin_i`: the rotational velocity projected along the line of sight (km/s)
* `epsilon` (optional numeric argument): numeric scalar giving the
  limb-darkening coefficient, default = 0.6 which is typical for photospheric
  lines.  The specific intensity \$I\$ at any angle \$\\theta\$ from the
  specific intensity \$I_{\\text{cen}}\$ at the center of the disk is given by:

\$\$ I = I_{\\text{cen}}\\cdot(1 - \\varepsilon\\cdot(1 - \\cos(\\theta))) \$\$

### Output ###

The 2-tuple (`velocity_grid`, `lsf`):

* `velocity_grid`: vector of velocity grid with the same number of elements as
  `lsf` (see below)
* `lsf`: the convolution kernel vector for the specified rotational velocity.
  The number of points in `lsf` will be always be odd (the kernel is symmetric)
  and equal to either `ceil(2*v_sin_i/delta_v)` or `ceil(2*v_sin_i/delta_v) + 1`,
  whichever number is odd.  Elements of `lsf` will always be of type
  `AbstractFloat`.  To actually compute the broadening, the spectrum should be
  convolved with the rotational `lsf`

### Example ###

Plot the line spread function for a star rotating at 90 km/s in velocity space
every 3 km/s.  Use [PyPlot.jl](https://github.com/stevengj/PyPlot.jl) for
plotting.

``` julia
using PyPlot
plot(lsf_rotate(3, 90)...)
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
lsf_rotate(Δ::Real, v::Real, ɛ::Real=0.6) =
    _lsf_rotate(promote(float(Δ), float(v), float(ɛ))...)


text:
#=

   Number of walks until no shoes

   From
   Gunnar Blom, Lars Holst, Dennis Sandell:
   "Problems and Snapshots from the World of Probability"
   Page 8f, Problem 1.6 Number of walks until no shoes

   """
   A has a house with one front door and one back door. He places
   two pairs of walking shoes at each door. For each walk, he selects
   one door at random, puts on a pair of shoes, returns after a walk
   to a randomly chosen door, and takes off the shoes at the door.
   We want to determine the expected number of finished walks until
   A discovers that no shoes are available at the door he has selected
   for his next walk.
   """

  For two pair of shoes, the theoretical probability. of winning for each strategy is 12.

   Summary Statistics
   parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

         len   12.1993   10.1602     0.1016    0.1316   6606.4119    0.9999     2942.7225

   Theoretical walks for n: 1:10:
   (1, 4)
   (2, 12)
   (3, 24)
   (4, 40)
   (5, 60)
   (6, 84)
   (7, 112)
   (8, 144)
   (9, 180)
   (10, 220)

  Cf  banachs_match_box_problem.jl
      ~/webppl/number_of_walks_until_no_shoes.wppl

=#

using Turing, StatsPlots
include("jl_utils.jl")

# The theoretical probability
function theoreticalProb2(n,i)
    2*n+4*n*i-2*i*i
end

@model function number_of_walks_until_no_shoes(n,i)

    # In both function we use rand(Bernoulli()), not
    #   flip()
    # since that don't give the proper solution.

    # Take shoes from a random door
    # Note: we only add to a when taking the shoes.
    function takeShoes(a,left,right)
        # Pick a door and check if there are any shoes
        pick = rand(Bernoulli()) == true ? 'l' : 'r'
        if (pick == 'l' && left == 0)  || (pick == 'r' && right == 0)
            return a
        else
            if pick == 'l'
                return leaveShoes(push!(a,pick),left-1,right)
            else
                return leaveShoes(push!(a,pick),left,right-1)
            end
        end
    end

    # Leave shoes at a random door
    function leaveShoes(a,left,right)
        if rand(Bernoulli()) == true
            return takeShoes(a,left+1,right)
        else
            return takeShoes(a,left,right+1)
        end
    end

    a = takeShoes([],n,i)

    len ~ Dirac(length(a))

    return a

end

n = 2
i = 2
# n = 10
# i = 10
println("Theoretical prob:", theoreticalProb2(n,i))
model = number_of_walks_until_no_shoes(n,i)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

# show_var_dist_pct(chns, :len)

# chains_params = Turing.MCMCChains.get_sections(chns, :parameters)
# genq = generated_quantities(model, chains_params)
# show_var_dist_pct(genq)
# println()

println("\nTheoretical walks for n: 1:10:")
for i in 1:10
    println((i,theoreticalProb2(i,i)))
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
struct Trans2DateTimeBoolVariate <: AbstractBoolVariate
    name::String
    basevariate1::AbstractDateTimeVariate
    basevariate2::AbstractDateTimeVariate
    transform::Function
end

Base.length(var::Trans2DateTimeBoolVariate) = length(var.basevariate1)

function slice(boolvar::Trans2DateTimeBoolVariate, fromobs::Integer, toobs::Integer, slicelength::Integer)
    base1 = boolvar.basevariate1
    base2 = boolvar.basevariate2
    f = boolvar.transform
    g = (ms1::Int64, ms2::Int64) -> begin
        ms1 == zero(Int64) || ms2 == zero(Int64) ? false : begin
            dt1 = Dates.epochms2datetime(ms1)
            dt2 = Dates.epochms2datetime(ms2)
            f(dt1, dt2)
        end
    end
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = zip(slice(base1, fromobs, toobs, slicelength), slice(base2, fromobs, toobs, slicelength))
    mapslice2(g, slices, slicelength, Bool)
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
using DataFrames
using DiffDRO
using DiffDRO.TestUtils: get_financial_data
using JuMP
using LinearAlgebra
using LineSearches
using Optim
using Random
using Pkg
using Statistics
using SCS
using Test

Pkg.add(url="https://github.com/andrewrosemberg/PortfolioOpt.jl.git")
using PortfolioOpt

@testset "DiffDRO" begin
    include("dro.jl")
    include("obj.jl")
end

max_stars_repo_name: molet/DiffDRO

max_issues_repo_name: molet/DiffDRO

max_forks_repo_name: molet/DiffDRO


text:
# Interpolation between different time bases

import Grid

function Grid.interp{S<:SecondT}(sig::RegularSignal, ti::AbstractVector{S})
    # Using a Regular Grid is a little awkward, but more performant; it assumes the original basis is 1:length(sig)
    t = sig.time
    # So convert ti to "indices" relative to the signal's time
    tidxs = Array(Float64, length(ti))
    for i=1:length(ti) tidxs[i] = (ti[i] - t[1])/step(t) + 1; end

    g = Grid.InterpGrid(sig.data, Grid.BCnan, Grid.InterpLinear)
    vi = [g[tidx, chan] for tidx in tidxs, chan=1:length(sig)]
    Signal(ti, reshape(vi, tuple(length(ti), size(sig)...)), sig.dims, sig.meta)
end

function Grid.interp{S<:SecondT}(sig::Signal, ti::AbstractVector{S})
    # Grid doesn't support multidimensional irregular interpolation, so work
    # around it for now
    gs = [Grid.InterpIrregular(sig.time, sig.data[:,chan], Grid.BCnan, Grid.InterpLinear) for chan=1:length(sig)]
    vi = eltype(sig.data)[g[t] for t in ti, g in gs]
    Signal(ti, reshape(vi, tuple(length(ti), size(sig)...)), sig.dims, sig.meta)
end


text:
#
# types.jl -
#
# Definitions of common type in `LinearInterpolators`.
#
#------------------------------------------------------------------------------
#
# This file is part of the LinearInterpolators package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2016-2021, Éric Thiébaut.
#

"""
    Limits{T<:AbstractFloat}
Abstract type for defining boundary conditions.
"""
abstract type Limits{T<:AbstractFloat} end # FIXME: add parameter S

"""

`FlatLimits` are boundary conditions that implies that extrapolated positions
yield the same value of the interpolated array at the nearest position.

"""
struct FlatLimits{T<:AbstractFloat} <: Limits{T}
    len::Int # length of dimension in interpolated array
    function FlatLimits{T}(len::Integer) where {T}
        @assert len ≥ 1
        new{T}(len)
    end
end

"""

`SafeFlatLimits` are boundary conditions that implies that extrapolated
positions yield the same value of the interpolated array at the nearest
position.  Compared to `FlatLimits`, the operations are "safe" in the sense
that no `InexactError` get thrown even for very large extrapolation distances.

"""
struct SafeFlatLimits{T<:AbstractFloat} <: Limits{T}
    inf::T   # bound for all neighbors before range
    sup::T   # bound for all neighbors after range
    len::Int # length of dimension in interpolated array
    function SafeFlatLimits{T}(inf, sup, len::Integer) where {T}
        @assert inf < sup
        @assert len ≥ 1
        new{T}(inf, sup, len)
    end
end


text:
using .MCMCChains

function convert_a3d(a3d_array, cnames, ::Val{:mcmcchains}; start=1)
  pi = filter(p -> length(p) > 2 && p[end-1:end] == "__", cnames)
  p = filter(p -> !(p in  pi), cnames)

  MCMCChains.Chains(a3d_array,
    cnames,
    Dict(
      :parameters => p,
      :internals => pi
    );
    start=start
  )
end


text:
struct DirichletMultinomial{T <: Real} <: DiscreteMultivariateDistribution
    n::Int
    α::Vector{T}
    α0::T

    function DirichletMultinomial{T}(n::Integer, α::Vector{T}) where T
        α0 = sum(abs, α)
        sum(α) == α0 || throw(ArgumentError("alpha must be a positive vector."))
        n > 0 || throw(ArgumentError("n must be a positive integer."))
        new{T}(Int(n), α, α0)
    end
end
DirichletMultinomial(n::Integer, α::Vector{T}) where {T <: Real} = DirichletMultinomial{T}(n, α)
DirichletMultinomial(n::Integer, α::Vector{T}) where {T <: Integer} = DirichletMultinomial(n, float(α))
DirichletMultinomial(n::Integer, k::Integer) = DirichletMultinomial(n, ones(k))

Base.show(io::IO, d::DirichletMultinomial) = show(io, d, (:n, :α,))

# Parameters
ncategories(d::DirichletMultinomial) = length(d.α)
length(d::DirichletMultinomial) = ncategories(d)
ntrials(d::DirichletMultinomial) = d.n
params(d::DirichletMultinomial) = (d.n, d.α)
@inline partype(d::DirichletMultinomial{T}) where {T<:Real} = T

# Statistics
mean(d::DirichletMultinomial) = d.α .* (d.n / d.α0)
function var(d::DirichletMultinomial{T}) where T <: Real
    v = fill(d.n * (d.n + d.α0) / (1 + d.α0), length(d))
    p = d.α / d.α0
    for i in eachindex(v)
        @inbounds v[i] *= p[i] * (1 - p[i])
    end
    v
end
function cov(d::DirichletMultinomial{T}) where T <: Real
    v = var(d)
    c = d.α * d.α'
    multiply!(c, -d.n * (d.n + d.α0) / (d.α0^2 * (1 + d.α0)))
    for i in 1:length(d)
        @inbounds c[i, i] = v[i]
    end
    c
end

# Evaluation
function insupport(d::DirichletMultinomial, x::AbstractVector{T}) where T<:Real
    k = length(d)
    length(x) == k || return false
    for xi in x
        (isinteger(xi) && xi >= 0) || return false
    end
    return sum(x) == ntrials(d)
end
function _logpdf(d::DirichletMultinomial{S}, x::AbstractVector{T}) where {T<:Real, S<:Real}
    c = lgamma(S(d.n + 1)) + lgamma(d.α0) - lgamma(d.n + d.α0)
    for j in eachindex(x)
        @inbounds xj, αj = x[j], d.α[j]
        c += lgamma(xj + αj) - lgamma(xj + 1) - lgamma(αj)
    end
    c
end

# Sampling
function _rand!(d::DirichletMultinomial, x::AbstractVector{T}) where T<:Real
    multinom_rand!(ntrials(d), rand(Dirichlet(d.α)), x)
end

# Fit Model
# Using https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2945396/pdf/nihms205488.pdf
struct DirichletMultinomialStats <: SufficientStats
    n::Int
    s::Matrix{Float64}  # s_{jk} = ∑_i x_{ij} ≥ (k - 1),  k = 1,...,(n - 1)
    tw::Float64
    DirichletMultinomialStats(n::Int, s::Matrix{Float64}, tw::Real) = new(n, s, Float64(tw))
end
function suffstats(::Type{DirichletMultinomial}, x::Matrix{T}) where T<:Real
    ns = sum(x, dims=1)  # get ntrials for each observation
    n = ns[1]       # use ntrails from first ob., then check all equal
    all(ns .== n) || error("Each sample in X should sum to the same value.")
    d, m = size(x)
    s = zeros(d, n)
    @inbounds for k in 1:n, i in 1:m, j in 1:d
        if x[j, i] >= k
            s[j, k] += 1.0
        end
    end
    DirichletMultinomialStats(n, s, m)
end
function suffstats(::Type{DirichletMultinomial}, x::Matrix{T}, w::Array{Float64}) where T<:Real
    length(w) == size(x, 2) || throw(DimensionMismatch("Inconsistent argument dimensions."))
    ns = sum(x, dims=1)
    n = ns[1]
    all(ns .== n) || error("Each sample in X should sum to the same value.")
    d, m = size(x)
    s = zeros(d, n)
    @inbounds for k in 1:n, i in 1:m, j in 1:d
        if x[j, i] >= k
            s[j, k] += w[i]
        end
    end
    DirichletMultinomialStats(n, s, sum(w))
end
function fit_mle(::Type{DirichletMultinomial}, ss::DirichletMultinomialStats;
                 tol::Float64 = 1e-8, maxiter::Int = 1000)
    k = size(ss.s, 2)
    α = ones(size(ss.s, 1))
    rng = 0.0:(k - 1)
    @inbounds for iter in 1:maxiter
        α_old = copy(α)
        αsum = sum(α)
        denom = ss.tw * sum(inv, αsum .+ rng)
        for j in eachindex(α)
            αj = α[j]
            num = αj * sum(vec(ss.s[j, :]) ./ (αj .+ rng))  # vec needed for 0.4
            α[j] = num / denom
        end
        maximum(abs, α_old - α) < tol && break
    end
    DirichletMultinomial(ss.n, α)
end


text:
include("rob_lib.jl")

function mark_it_all(r::Robot)
    num_verticale = digit_movements!(r, Sud)
    num_horizont = digit_movements!(r, West)

    side = Nord

    # while we dont reach the corner
    while !(isborder(r, Ost) && (isborder(r, Ost) || isborder(r, Sud)))
        putmarker!(r)
        putmarkers!(r, side)
        move!(r, Ost)
        side = inverse(side)
    end

    putmarker!(r)
    putmarkers!(r, side)

    digit_movements!(r, Sud)
    digit_movements!(r, West)

    Back!(r, Nord, num_verticale)
    Back!(r, Ost, num_horizont)
end

max_stars_repo_name: NikitaDestroyer/Basics

max_issues_repo_name: NikitaDestroyer/Basics

max_forks_repo_name: NikitaDestroyer/Basics


text:
# core

Base.length(dict::CardinalDict{V,K}) where {K,V} = sum(dict.present)
Base.endof(dict::CardinalDict{V,K}) where {K,V} = findlast(dict.present, true)
Base.isempty(dict::CardinalDict{V,K}) where {K,V} = sum(dict.present) == 0
isfull(dict::CardinalDict{V,K}) where {K,V} = all(dict.present)

Base.eltype(dict::CardinalDict{V,K}) where {K,V} = Pair{K,V}

Base.:(==)(a_dict::D, b_dict::D) where D<:CardinalDict{V,K} where {K,V} =
    keymax(a_dict) == keymax(b_dict) && keys(a_dict) == keys(b_dict) && content(a_dict) == content(b_dict)
Base.:(==)(a_dict::CardinalDict{V,K}, b_dict::CardinalDict{J,W}) where {J,W,K,V} =
    false

Base.:(!=)(a_dict::D, b_dict::D) where D<:CardinalDict{V,K} where {K,V} =
    keymax(a_dict) != keymax(b_dict)  || keys(a_dict) != keys(b_dict) || content(a_dict) != content(b_dict)
Base.:(!=)(a_dict::CardinalDict{V,K}, b_dict::CardinalDict{J,W}) where {J,W,K,V} =
    true

# keys, content, in (∈, ∋, ∉, ∌)

function tally_active(dict::CardinalDict{V,K}) where {K,V}
    return sum(dict.present)
end
function tally_passive(dict::CardinalDict{V,K}) where {K,V}
    return length(dict.present) - sum(dict.present)
end

function Base.keys(dict::CardinalDict{V,K}) where {K,V}
    result = Vector{K}(tally_active(dict))
    allkeys = one(K):keymax(dict)
    result[:] = allkeys[dict.present]
    return result
end

function Base.values(dict::CardinalDict{V,K}) where {K,V}
    result = Vector{V}(tally_active(dict))
    allkeys = keys(dict)
    for (i,ky) in enumerate(allkeys)
        result[i] = getindex(dict, ky)
    end
    return result
end

function Base.in(value::V, dict::CardinalDict{V,K}) where {K,V}
    return length(dict)!=0 && in(value, content(dict))
end

#@inline Base.:(∈)(value::V, dict::CardinalDict{V,K}) where {K,V} = in(value, dict) # not needed
@inline Base.:(∋)(dict::CardinalDict{V,K}, value::V) where {K,V} = in(value, dict)
@inline Base.:(∉)(value::V, dict::CardinalDict{V,K}) where {K,V} = !in(value, dict)
@inline Base.:(∌)(dict::CardinalDict{V,K}, value::V) where {K,V} = !in(value, dict)

# delete!, clearindex!, empty!

function Base.delete!(dict::CardinalDict{V,K}, key::K) where {K,V}
    check_index(dict, key)
    setindex!(dict.present, false, key)
    return dict
end
@inline Base.delete!(dict::CardinalDict{V,K}, key::J) where {J,K,V} =
    delete!(dict, key%K)

function Base.empty!(dict::CardinalDict{V,K}) where {K,V}
    0 < key <= keymax(dict) || throw(ErrorException("Key (index) $(key) is outside of the domain 1:$(keymax(dict))."))
    dict.present[:] = falses(length(dict.present))
    return dict
end

# iteration

function Base.start(dict::CardinalDict{V,K}) where {K,V}
    (1, keys(dict))
end

function Base.next(dict::CardinalDict{V,K}, state) where {K,V}
    index, ks = state
    (ks[index], dict.present[index]), (index+1, ks)
end

function Base.done(dict::CardinalDict{V,K}, state) where {K,V}
    state[1] > length(dict)
end

# nonstandard

keymax(dict::CardinalDict{V,K}) where {K,V} = K(length(dict.content))

function clearindex!(dict::CardinalDict{V,K}, key::K) where K where V
    check_index(dict, key)
     @inbounds setindex!(dict.present, false, key)
    return nothing
end
@inline clearindex!(dict::CardinalDict{V,K}, key::J) where J where K where V =
    clearindex!(dict, key%K)


text:
# verify that solve_tr_subproblem! finds the minimum
using NearlyNewton, Test
using LinearAlgebra: norm, dot

trsolver = NWI()

n = 2
gr = [-0.74637,0.52388]
# H is not posdef!
H = [0.945787 -3.07884; -3.07884 -1.27762]

s = zeros(n)
Δ = 1.0
solution = trsolver(gr, H, Δ, s ;maxiter=100)
println("Trust region : $Δ")
println("Norm of step : $(norm(solution.p,2))")

for j in 1:100000
    bad_s = rand(n)
    bad_s ./= norm(bad_s)  # boundary
    model(s2) = dot(gr, s2) + (dot(s2, H * s2))/2
    @test model(s) <= model(bad_s) + 1e-8
end

# H is now posdef!
Hp = H*H
s = zeros(n)
Δ = 1.0
solution = trsolver(gr, Hp, Δ, s ;maxiter=100)
println("Trust region : $Δ")
println("Norm of step : $(norm(solution.p,2))")
println("Is it interior? $(solution.interior)")

for j in 1:100000
    bad_s = rand(n)
    if norm(bad_s) > Δ
        bad_s ./= norm(bad_s)  # boundary
    end
    model(s2) = dot(gr, s2) + (dot(s2, Hp * s2))/2
    @test model(s) <= model(bad_s) + 1e-8
end

Hhard = [-0.945787 0.0; 0.0 -1.27762]
Δ = 1.0
solution = trsolver(gr, Hhard, Δ, s ;maxiter=100)

or

# From tR book p180
s2 = zeros(4)
gr2 = [0.0, 1.0, 1.0, 1.0]
Hhard2 = [-2.0 0.0 0.0 0.0;0.0 -1.0 0.0 0.0;0.0 0.0 0.0 0.0;0.0 0.0 0.0 1.0]
solution = trsolver(gr2,  Hhard2, Δ, s2 ;maxiter=100)


text:
struct Integrator

end


text:
using TensorMatrices_lemon
using Base.Test

# creating a tenmat
#=
		"""
			matrix = rand(6,8)
			rowindex = [3,1]
			colindex = [2,4]
			tensorsize = (2,4,3,2)
			tenmat = Tenmat(matrix, rowindex, colindex, tensorsize)
		"""
		""" [test]
			tensorSize(tenmat)  == (2,4,3,2)

"""

# initialization
"""
	A = rand(4,7,2,9,6)
	rowindex = [5,2]
	colindex = [3,4,1]
	tenmat = Tenmat(A,rowindex,colindex)

	A = rand(4,7,6);
	tenmat = Tenmat(A,[3,1],[2])
	size(tenmat) == (24,7) || println("wrong")
	size(tenmat,2) == 7 || println("wrong")
	getindex(tenmat, 2, 6) == getindex(tenmat.matrix, 2, 6) || println("wrong")
	setindex!(tenmat, 100.0 , 2, 6)
	getindex(tenmat, 2, 6) == 100.0 || println("wrong")
	similar(tenmat, Int)
	similar(tenmat, Int, (5,2,7))
	similar(tenmat, Int, [1,3],[2])
	similar(tenmat, Int, [4,2,5], [1,3], (2,4,3,1,6))
	similar(tenmat)
	similar(tenmat, (5,2,7))
	similar(tenmat, [1,3],[2])
	similar(tenmat, [4,2,5], [1,3], (2,4,3,1,6))

"""
=#


text:
include(joinpath(@__DIR__, "config.jl"))
using JuMP, Gurobi, MLKernels

solver = with_optimizer(Gurobi.Optimizer; OutputFlag=0, Threads=1)

### data settings ###
data_sets = ["Annthyroid", "Arrhythmia", "Arrhythmia", "Cardiotocography", "Glass",
             "HeartDisease", "Hepatitis", "InternetAds", "Ionosphere", "PageBlocks",
             "Parkinson", "Pima", "SpamBase", "Stamps", "WPBC", "Wilt"]
data_dirs = [joinpath("dami", x) for x in data_sets]
data_info = "dami_outlier_detection"

### learning scenario ###
initial_pool_strategy = [("Pnin", Dict(:n => 25))]
split_strategy = [("Sl", Dict())]
initial_sample_strat = :nn_biased_sample

#### models ####
models = [Dict(:type => :SVDDneg, :param => Dict{Symbol, Any}())]
init_strategies = [WangCombinedInitializationStrategy(solver, 2.0.^range(-4, stop=4, step=1.0), FixedCStrategy(1.0))]
classify_precision = SVDD.OPT_PRECISION

#### oracle ####
oracle_param = [
    Dict{Symbol, Any}(
        :type => QuerySynthesisCVWrapperOracle,
        :param => Dict{Symbol, Any}(
            :subtype => QuerySynthesisSVMOracle,
        )
    ),
]

#### query strategies ####
basline_query_strategies = [
    Dict(:type => :RandomQss, :param => Dict{Symbol, Any}(:epsilon => 0.1)),
    Dict(:type => :RandomOutlierQss, :param => Dict{Symbol, Any}(:epsilon => 0.1))
]

fancy_query_strategies = [
    Dict(:type => :DecisionBoundaryQss, :param => Dict{Symbol, Any}()),
    Dict(:type => :ExplorativeMarginQss, :param => Dict{Symbol, Any}(:solver => solver, :lambda => 1.0, :use_penalty => false)),
    Dict(:type => :ExplorativeMarginQss, :param => Dict{Symbol, Any}(:solver => solver, :lambda => 1.0, :use_penalty => true)),
]

qs_optimizers = [
    ParticleSwarmOptimization(eps=1.0, swarmsize=1_000, maxiter=1_000),
    BlackBoxOptimization(:dxnes, eps=1.0),
]

query_strategies = basline_query_strategies

for qs in fancy_query_strategies
    for o in qs_optimizers
        qs_parametrized = deepcopy(qs)
        qs_parametrized[:param][:optimizer] = o
        push!(query_strategies, qs_parametrized)
    end
end

num_al_iterations = 100
num_resamples_initial_pool = 10
exp_name = "eval_part_2_qss"

max_stars_repo_name: englhardt/des-evaluation

max_issues_repo_name: englhardt/des-evaluation

max_forks_repo_name: englhardt/des-evaluation


text:
"""
	function build_mlp(idim::Int, hdim::Int, odim::Int, nlayers::Int; activation::String="relu", lastlayer::String="")
	Creates a chain with `nlayers` layers of `hdim` neurons with transfer function `activation`.
	input and output dimension is `idim` / `odim`
	If lastlayer is no specified, all layers use the same function.
	If lastlayer is "linear", then the last layer is forced to be Dense.
	It is also possible to specify dimensions in a vector.

```juliadoctest
julia> build_mlp(4, 11, 1, 3, activation="relu")
Chain(Dense(4, 11, relu), Dense(11, 11, relu), Dense(11, 1, relu))

julia> build_mlp([4, 11, 11, 1], activation="relu")
Chain(Dense(4, 11, relu), Dense(11, 11, relu), Dense(11, 1, relu))

julia> build_mlp(4, 11, 1, 3, activation="relu", lastlayer="tanh")
Chain(Dense(4, 11, relu), Dense(11, 11, relu), Dense(11, 1, tanh))
```
"""
build_mlp(ks::Vector{Int}, fs::Vector) = Flux.Chain(map(i -> Dense(i[2],i[3],i[1]), zip(fs,ks[1:end-1],ks[2:end]))...)

build_mlp(idim::Int, hdim::Int, odim::Int, nlayers::Int; kwargs...) =
	build_mlp(vcat(idim, fill(hdim, nlayers-1)..., odim); kwargs...)

function build_mlp(ks::Vector{Int}; activation::String = "relu", lastlayer::String = "")
	activation = (activation == "linear") ? "identity" : activation
	fs = Array{Any}(fill(eval(:($(Symbol(activation)))), length(ks) - 1))
	if !isempty(lastlayer)
		fs[end] = (lastlayer == "linear") ? identity : eval(:($(Symbol(lastlayer))))
	end
	build_mlp(ks, fs)
end

"""
	vec_size(conv_dims, scalings, channel)

Computes the size of vectorized ouput, given size of input, convolutional scalings and last channel.
"""
function vec_size(conv_dims, scalings, channel)
	ho = conv_dims[1]/(reduce(*, scalings)) # height before reshaping
	wo = conv_dims[2]/(reduce(*, scalings)) # width before reshaping
	# this ensures size compatibility, ho/wo are Ints
	(ho == floor(Int, ho)) ? ho = floor(Int, ho) : error("your input size and scaling is not compatible")
	(wo == floor(Int, wo)) ? wo = floor(Int, wo) : error("your input size and scaling is not compatible")
	ho*wo*channel, ho, wo
end

# these two functions are taken from GMExtensions.jl, but I wanted to avoid the dependency since
# it is an unregistered package
# TODO: maybe use this in GANomaly constructor?
"""
	conv_encoder(idim, zdim, kernelsizes, channels, scalings[; activation, densedims])

Constructs a convolutional encoder.

# Arguments
	- `idim`: size of input - (h,w,c)
	- `zdim`: latent space dimension
	- `kernelsizes`: kernelsizes for conv. layers (only odd numbers are allowed)
	- `channels`: channel numbers
	- `scalings`: scalings vector
	- `activation`: default relu
	- `densedims`: if set, more than one dense layers are used

# Example
```julia-repl
julia> encoder = conv_encoder((64, 48, 1), 2, (3, 5, 5), (2, 4, 8), (2, 2, 2), densedims = (256,))
Chain(Conv((3, 3), 1=>2, relu), MaxPool((2, 2), pad = (0, 0, 0, 0), stride = (2, 2)), Conv((5, 5), 2=>4, relu), MaxPool((2, 2), pad = (0, 0, 0, 0), stride = (2, 2)), Conv((5, 5), 4=>8, relu), MaxPool((2, 2), pad = (0, 0, 0, 0), stride = (2, 2)), #44, Dense(384, 2))
julia> encoder(randn(Float32, 64, 48, 1, 2))
2×2 Array{Float32,2}:
  0.247844   0.133781
 -0.605763  -0.494911
```
"""
function conv_encoder(idim::Union{Tuple, Vector}, zdim::Int, kernelsizes::Union{Tuple, Vector},
	channels::Union{Tuple, Vector}, scalings::Union{Tuple, Vector};
	activation = "relu", densedims::Union{Tuple, Vector} = [], batchnorm=false)
	nconv = length(kernelsizes)
	(nconv == length(channels) == length(scalings)) ? nothing : error("incompatible input dimensions")
	(length(idim) == 3) ? nothing : error("idim must be (h, w, c)")
	# also check that kernelsizes are all odd numbers
	(all(kernelsizes .% 2 .== 1)) ? nothing : error("not implemented for odd kernelsizes")

	# string to function
	actf = eval(:($(Symbol(activation))))

	# initialize dimensions
	cins = vcat(idim[3], channels[1:end-1]...) # channels in
	couts = channels # channels out
	din, _, _ = vec_size(idim, scalings, channels[end]) # the size of vectorized output after all convolutions

	# now build a vector of layers to be used later
	layers = Array{Any,1}()
	# first add the convolutional and maxpooling layers
	for (k, ci, co, s) in zip(kernelsizes, cins, couts, scalings)
		pad = Int((k-1)/2)
		# paddding so that input and output size are same
		push!(layers, Conv((k,k), ci=>co, actf; pad = (pad, pad)))
		push!(layers, MaxPool((s,s)))
		# generally, bn is not recommended for training vaes
		batchnorm ? push!(layers, BatchNorm(co)) : nothing
	end

	# reshape
	push!(layers, x -> reshape(x, din, :))

	# and dense layers
	ndense = length(densedims)
	dins = vcat(din, densedims...)
	douts = vcat(densedims..., zdim)
	dacts = vcat([actf for _ in 1:ndense]..., identity)
	for (_di, _do, _da) in zip(dins, douts, dacts)
		push!(layers, Dense(_di, _do, _da))
	end

	Flux.Chain(layers...)
end

"""
	conv_decoder(odim, zdim, kernelsizes, channels, scalings[; activation, densedims, vec_output,
		vec_output_dim])

Constructs a convolutional encoder.

# Arguments
- `odim`: size of outnput - (h,w,c)
- `zdim`: latent space dimension
- `kernelsizes`: kernelsizes for conv. layers (only odd numbers are allowed)
- `channels`: channel numbers
- `scalings`: scalings vector
- `activation`: default relu
- `densedims`: if set, more than one dense layers are used
- `vec_output`: output is vectorized (default false)
- `vec_output_dim`: determine what the final size of the vectorized input should be (e.g. add extra dimensions for variance estimation)

# Example
```julia-repl
julia> decoder = conv_decoder((64, 48, 1), 2, (5, 5, 3), (8, 4, 2), (2, 2, 2); densedims = (256,))
Chain(Dense(2, 256, relu), Dense(256, 384, relu), #19, ConvTranspose((2, 2), 8=>8, relu), Conv((5, 5), 8=>4, relu), ConvTranspose((2, 2), 4=>4, relu), Conv((5, 5), 4=>2, relu), ConvTranspose((2, 2), 2=>2, relu), Conv((3, 3), 2=>1))
julia> y = decoder(randn(Float32, 2, 2));
julia> size(y)
(64, 48, 1, 2)
```
"""
function conv_decoder(odim::Union{Tuple, Vector}, zdim::Int, kernelsizes::Union{Tuple, Vector},
	channels::Union{Tuple, Vector}, scalings::Union{Tuple, Vector};
	activation = relu, densedims::Union{Tuple, Vector} = [],
	vec_output = false, vec_output_dim = nothing, batchnorm=false)
	nconv = length(kernelsizes)
	(nconv == length(channels) == length(scalings)) ? nothing : error("incompatible input dimensions")
	(length(odim) == 3) ? nothing : error("odim must be (h, w, c)")
	# also check that kernelsizes are all odd numbers
	(all(kernelsizes .% 2 .== 1)) ? nothing : error("not implemented for odd kernelsizes")

	# string to function
	actf = eval(:($(Symbol(activation))))

	# initialize some stuff
	cins = channels # channels in
	couts = vcat(channels[2:end]..., odim[3]) # channels out
	dout, ho, wo = vec_size(odim, scalings, channels[1]) # size of vectorized output after convolutions

	# now build a vector of layers to be used later
	layers = Array{Any,1}()

	# start with dense layers
	ndense = length(densedims)
	dins = vcat(zdim, densedims...)
	douts = vcat(densedims..., dout)
	for (_di, _do) in zip(dins, douts)
		push!(layers, Dense(_di, _do, actf))
	end

	# reshape
	push!(layers, x -> reshape(x, ho, wo, channels[1], :))

	# add the transpose and convolutional layers
	acts = vcat([actf for _ in 1:nconv-1]..., identity)
	for (k, ci, co, s, act) in zip(kernelsizes, cins, couts, scalings, acts)
		pad = Int((k-1)/2)
		# use convtranspose for upscaling - there are other posibilities, however this seems to be ok
		push!(layers, ConvTranspose((s,s), ci=>ci, actf, stride = (s,s)))
		push!(layers, Conv((k,k), ci=>co, act; pad = (pad, pad)))
		# generally, bn is not recommended for training vaes
		batchnorm ? push!(layers, BatchNorm(co)) : nothing
	end

	# if you want vectorized output
	if vec_output
		push!(layers, x -> reshape(x, reduce(*, odim), :))
		# if you want the final vector to be a different size than the output
		# (e.g. you want extra dimensions for variance estimation)
		(vec_output_dim == nothing) ? nothing : push!(layers, Dense(reduce(*, odim), vec_output_dim))
	end

	Flux.Chain(layers...)
end


text:
Color = Tuple{Float32, Float32, Float32, Float32}

mutable struct Renderer <: AbstractRenderer
	camera::Camera
	resources::Dict{Symbol, WeakRef}
	renderState::RenderStateHolder
	toRender::Vector{Renderable}
	sortFunc::Function
	clearColor::Union{Color, Nothing}
	clearStencil::Union{Int, Nothing}
	clearDepth::Union{Float64, Nothing}
	shaderPreamble::Array{UInt8}

	Renderer() = new(Camera(), Dict{Symbol, WeakRef}(), RenderStateHolder(), Vector{Renderable}(undef, 0), identity, (0f0, 0f0, 0f0, 1f0), 0, 1.0, UInt8[])
end

function init(renderer::Renderer)
end

function done(renderer::Renderer)
	while !isempty(renderer.resources)
		id, res = first(renderer.resources)
		if res.value != nothing
		  done(res.value)
		else
		  delete!(renderer.resources, id)
		end
	end
end

function add_renderer_resource(resource::Resource)
	renderer = resource.renderer
	if haskey(renderer.resources, resource.id)
		resource.id = Symbol("$(resource.id)_$(object_id(resource))")
		@assert !haskey(renderer.resources, resource.id)
	end
	renderer.resources[resource.id] = WeakRef(resource)
end

function remove_renderer_resource(resource::Resource)
	delete!(resource.renderer.resources, resource.id)
end

function has_resource(renderer::Renderer, id::Symbol)
	!haskey(renderer.resources, id) && return false
	renderer.resources[id].value != nothing && return true
	delete!(renderer.resources, id)
	return false
end

get_resource(renderer::Renderer, id::Symbol) = renderer.resources[id].value

function add(renderer::Renderer, renderable::Renderable)
	push!(renderer.toRender, renderable)
end

set_clear_color(renderer::Renderer, c) = renderer.clearColor = c
set_clear_stencil(renderer::Renderer, s) = renderer.clearStencil = s
set_clear_depth(renderer::Renderer, d) = renderer.clearDepth = d

function set_shader_preamble(renderer::Renderer, preamble::AbstractArray{UInt8})
	renderer.shaderPreamble = preamble[:]
end

function render_frame(renderer::Renderer)
	gl_clear_buffers(renderer.clearColor, renderer.clearDepth, renderer.clearStencil)
	# cull
	frustum = getfrustum(renderer.camera)
	filter!(r->!Shapes.outside(frustum, getbound(r)), renderer.toRender)
	#filter!(r->Shapes.intersect(frustum, getbound(r)), renderer.toRender)
	# sort
	if renderer.sortFunc != identity
		sort!(renderer.toRender, lt = renderer.sortFunc)
	end
	# render
	for r in renderer.toRender
		render(r, renderer)
	end

	empty!(renderer.toRender)
end

function apply(holder::RenderStateHolder, renderer::Renderer)
	for state in values(holder.states)
		set_and_apply(renderer.renderState, state)
	end
end


text:
# coding: utf-8

# In[1]:

#details of insertion sort can be found in the following link
# https://github.com/je-suis-tm/search-and-sort/blob/master/bubble%2C%20selection%20and%20insertion%20sort.jl
function insertion_sort(arr)

    for i in 2:length(arr)

        target=arr[i]

        j=i

        while (j!=1) && (target<arr[j-1])

            arr[j]=arr[j-1]
            j-=1
        end

        arr[j]=target

    end

    return arr

end

# In[2]:

#shell sort is a variation of insertion sort
#slicing method [i:group:length(arr)] is the key
#we gotta cut the original list into a few small groups

#assume we have a list of n elements
#the first step, we get a few sublists by [i:group:length(arr)]
#we apply insertion sort on those sublists and concatenate
#next, we do group/default
#and we obtain a few new sublists by [i:group:length(arr)]
#we apply insertion sort on new sublists and concatenate
#we keep using /default method and do insertion sort and concatenate
#til group reaches zero
#we concatenate sublists and do a final insertion sort
#we shall end up with a sorted list
#the time complexity is quite difficult to calculate
#it falls somewhere between linear and quadratic

# In[3]:

function shell_sort(arr,default=3)

    #the first step is to initialize default value
    #we will use this variable to divide the list and do slicing
    #in this case, default equals to 3
    #you can change the default number to any digit
    #bear in mind that this variable keeps the size of each sublist reasonable
    #we keep group divided by 3 until it reaches zero
    group=Int32(floor(length(arr)/default))

    while group>0

        for i in 1:group

            #equivalent to [i::group] in python
            subset=arr[i:group:length(arr)]

            #instead of reinventing the wheel
            #i directly call it in shell sort
            subset=insertion_sort(subset)

            arr[i:group:length(arr)]=subset

        end

        group=Int32(floor(group/default))

    end

    return insertion_sort(arr)

end

# In[4]:

for _ in 1:100

    test_arr=rand(1000)

    if !(shell_sort(test_arr)==sort(test_arr))
        printstyled("Erreur",color=:red)
    end

end

max_stars_repo_path: shell sort.jl

max_stars_repo_name: je-suis-tm/search-and-sort

max_issues_repo_path: shell sort.jl

max_issues_repo_name: HafizSabia/search-and-sort

max_forks_repo_path: shell sort.jl

max_forks_repo_name: HafizSabia/search-and-sort


text:
module StatefulIterators

export StatefulIterator, ArrayIterator, peek, available, reset!

import Base: start, next, done, read, read!, readbytes,
             position, seek, seekstart, seekend, skip,
             eof, copy, eltype,
             readuntil, readline

# --- the type + basic operations

type StatefulIterator{T,S}
    iter::T
    state::S
end

StatefulIterator(iter) = StatefulIterator(iter, start(iter))
StatefulIterator(s::StatefulIterator) = StatefulIterator(s.iter, s.state)

copy(s::StatefulIterator) = StatefulIterator(s)

eltype{T,S}(::StatefulIterator{T,S}) = eltype(T)

macro unlimited_loop(s, i, x, block)
    quote
        $(esc(i)) = 0
        for $(esc(x)) in $(esc(s))
            $(esc(i)) += 1
            $(esc(block))
        end
    end
end

macro limited_loop(s, i, n, x, block)
    quote
        $(esc(i)) = 0
        if $(esc(n)) > 0
            for $(esc(x)) in $(esc(s))
                $(esc(i)) += 1
                $(esc(block))
                $(esc(i)) >= $(esc(n)) && break
            end
            $(esc(i)) < $(esc(n)) && throw(EOFError())
        end
    end
end

macro preserving_state(s, block)
    quote
        save = $(esc(s)).state
        try
            $(esc(block))
        finally
            $(esc(s)).state = save
        end
    end
end

# --- iter protocol

start(s::StatefulIterator) = s

next{T,S}(s::StatefulIterator{T,S}, ::StatefulIterator{T,S}) = next(s)
function next(s::StatefulIterator)
    item, s.state = next(s.iter, s.state)
    item, s
end

done(s::StatefulIterator, ::StatefulIterator) = done(s)
done(s::StatefulIterator) = done(s.iter, s.state)

# --- basic stream-like extensions

assertNotDone(s::StatefulIterator) = done(s) && throw(EOFError())

"""read is like next() but returns only the value.  similar to next(),
it should be used with done() or you will get an error when reading
past the end of the collection.  the error is EOFError to be
consistent with streams."""
read(s::StatefulIterator) = (assertNotDone(s); next(s)[1])

compatible(U, T) = U == Any || eltype(T) <: U

function read!{T,S,U}(s::StatefulIterator{T,S}, a::Array{U})
    n = length(a)
    if compatible(U, T)
        @limited_loop s i n x begin
            @inbounds a[i] = x
        end
    else
        b = reinterpret(eltype(T), a)
        m = length(b)
        @limited_loop s i m x begin
            @inbounds b[i] = x
        end
        t, u = sizeof(eltype(T)), sizeof(U)
        if t > u
            # need to fill in the Us that didn't cover an entire eltype(T)
            nn = fld(m * t, u)
            if n > nn
                c = read(s, U, n-nn)
                for i in 1:(n-nn)
                    a[nn+i] = c[i]
                end
            end
        end
    end
end

function readbytes(s::StatefulIterator, nb=typemax(Int))
    b = Array(UInt8, 0)
    readbytes!(s, b, nb)
    b
end

function readbytes!{T,S}(s::StatefulIterator{T,S}, b::Array{UInt8}, nb=tyemax(Int))
    try
        n = min(nb, available(s, UInt8))
        if length(b) < n
            resize!(b, n)
        end
        a = reinterpret(eltype(T), a)
        m = length(a)
        @limited_loop s i m x begin
            @inbounds a[i] = x
        end
    catch  # available failed
        t = sizeof(eltype(T))
        i, a = 0, Array(UInt8, t)
        while !done(s) && i < nb
            read!(s, a)
            j = 0
            while i < nb && j < t
                i += 1
                j += 1
                if i > length(b)
                    push!(b, a[j])
                else
                    b[i] = a[j]
                end
            end
        end
    end
end

"""the next() value, without advancing the state.  similar to next(),
it should be used with done() or you will get an error when reading
past the end of the collection."""
peek(s::StatefulIterator) = (assertNotDone(s); next(s.iter, s.state)[1])

position(s::StatefulIterator) = s.state

seek{T,S}(s::StatefulIterator{T,S}, p::S) = s.state = p

seekstart(s::StatefulIterator) = s.state = start(s.iter)

function seekend(s::StatefulIterator)
    @unlimited_loop s i _ begin end
end

"""the number of read() values remaining."""
function available(s::StatefulIterator)
    @preserving_state s begin
        @unlimited_loop s i _ begin end
        i
    end
end

function skip(s::StatefulIterator, offset)
    @limited_loop s i offset x begin end
end

# we do not support mark / unmark / reset

eof(s::StatefulIterator) = done(s)

# TODO - support serialize / deserialize

# --- arbitrary types and arrays of data

function read(s::StatefulIterator, dims::Int...)
    a = Array(eltype(s), dims)
    n = prod(dims)
    @limited_loop s i n x begin
        @inbounds a[i] = x
    end
    a
end

function read{T,S,U}(s::StatefulIterator{T,S}, ::Type{U})
    if compatible(U, T)
        read(s)
    else
        n = cld(sizeof(U), sizeof(eltype(T)))
        reinterpret(U, read(s, n))[1]
    end
end

function read{T,S,U}(s::StatefulIterator{T,S}, ::Type{U}, dims...)
    if compatible(U, T)
        read(s, dims...)
    else
        m = prod(dims)
        n = cld(m * sizeof(U), sizeof(eltype(T)))
        reshape(resize!(reinterpret(U, read(s, n)), m), dims)
    end
end

function peek(s::StatefulIterator, dims::Int...)
    @preserving_state s begin
        read(s, dims...)
    end
end

function peek(s::StatefulIterator, t::Type)
    @preserving_state s begin
        read(s, t)
    end
end

function peek(s::StatefulIterator, t::Type, dims...)
    @preserving_state s begin
        read(s, t, dims...)
    end
end

function available{T,S,U}(s::StatefulIterator{T,S}, ::Type{U})
    if compatible(U, T)
        available(s)
    else
        t, u = sizeof(eltype(T)), sizeof(U)
        fld(available(s) * t, u)
    end
end

function skip{T,S,U}(s::StatefulIterator{T,S}, ::Type{U}, offset)
    if compatible(U, T)
        skip(s, offset)
    else
        t, u = sizeof(eltype(T)), sizeof(U)
        skip(s, offset * cld(u, t))
    end
end

# --- optimisations for known types

typealias LinearIndexed Union{Array, LinSpace, UnitRange, ASCIIString}

seekend{T<:LinearIndexed,S<:Integer}(s::StatefulIterator{T,S}) = s.state = length(s.iter) + 1

seekend{T<:UnitRange,S<:Integer}(s::StatefulIterator{T,S}) = s.state = s.iter.stop + 1

endof(r::StepRange) = r.start + r.step * (fld(r.stop - r.start, r.step) + 1)

seekend{T<:StepRange,S<:Integer}(s::StatefulIterator{T,S}) = s.state = endof(s.iter)

available{T<:LinearIndexed,S<:Integer}(s::StatefulIterator{T,S}) = length(s.iter) - s.state + 1

available{T<:UnitRange,S<:Integer}(s::StatefulIterator{T,S}) = s.iter.stop - s.state + 1

available{T<:StepRange,S<:Integer}(s::StatefulIterator{T,S}) = fld(s.iter.stop - s.state, s.iter.step) + 1

limit(n, lim) = n > lim ? throw(EOFError()) : n

skip{T<:LinearIndexed,S<:Integer}(s::StatefulIterator{T,S}, offset) = s.state = limit(s.state + offset, length(s.iter) + 1)

skip{T<:UnitRange,S<:Integer}(s::StatefulIterator{T,S}, offset) = s.state = limit(s.state + offset, s.iter.stop + 1)

skip{T<:StepRange,S<:Integer}(s::StatefulIterator{T,S}, offset) = s.state = limit(s.state + offset * s.iter.step, endof(s.iter))

# we don't need to add the versions with an explicit type - like
# skip(s, type, offset) - because the implementations earlier will
# delegate to these optimized versions.

# --- avoid known bad types

# we cannot peek iters that are themselves stateful (because the state
# is not actually a separate state).
for X in (Task, StatefulIterator)
    msg = "$X lacks explicit state"
    @eval copy{T<:$X,S}(s::StatefulIterator{T,S}) = error($msg)
    @eval peek{T<:$X,S}(s::StatefulIterator{T,S}, args...) = error($msg)
    @eval position{T<:$X,S}(s::StatefulIterator{T,S}) = error($msg)
    @eval seek{T<:$X,S}(s::StatefulIterator{T,S}, p) = error($msg)
    @eval seekstart{T<:$X,S}(s::StatefulIterator{T,S}) = error($msg)
    @eval available{T<:$X,S}(s::StatefulIterator{T,S}, args...) = error($msg)
end

# --- text-specific routines

# julia-0.4/base/io.jl
function readuntil{T<:AbstractString,S}(s::StatefulIterator{T,S}, t::AbstractString)
    l = length(t)
    if l == 0
        return ""
    elseif l > 40
        warn("readuntil(SteatfulIterator,AbstractString) will perform poorly with a long string")
    end
    out = IOBuffer()
    m = Array(Char, l)  # last part of stream to match
    t = collect(t)
    i = 0
    while !eof(s)
        i += 1
        c = read(s, Char)
        write(out, c)
        if i <= l
            m[i] = c
        else
            m[1:l-1] = m[2:l]
            m[l] = c
        end
        if i >= l && m == t
            break
        end
    end
    return takebuf_string(out)
end

readline{T<:AbstractString,S}(s::StatefulIterator{T,S}) = readuntil(s, "\n")

end


text:
@testset "Empirical" begin
  # homogeneous field has zero variogram
  γ = EmpiricalVariogram(Matrix(1.0I, 3, 3), ones(3), nlags=2, maxlag=2.)
  x, y, n = values(γ)
  @test x ≈ [1/2, 3/2]
  @test isnan(y[1]) && y[2] == 0.
  @test n == [0, 3]

  # test spatial data interface
  sdata = PointSetData(Dict(:z => [1.,0.,1.]), [25. 50. 75.; 25. 75. 50.])
  γ = EmpiricalVariogram(sdata, :z, nlags=20, maxlag=1.)
  x, y, n = values(γ)
  @test length(x) == 20
  @test length(y) == 20
  @test length(n) == 20

  # empirical variogram on integer coordinates
  γ = EmpiricalVariogram(Matrix(1I, 3, 3), ones(3), nlags=2, maxlag=2)
  x, y, n = values(γ)
  @test x ≈ [1/2, 3/2]
  @test isnan(y[1]) && y[2] == 0.
  @test n == [0, 3]

  # empirical variogram with only missing data
  X = rand(2,3); z = [NaN, NaN, NaN]
  γ = EmpiricalVariogram(X, z, maxlag=1., nlags=5)
  x, y, n = values(γ)
  @test x == [.1, .3, .5, .7, .9]
  @test all(isnan.(y))
  @test all(iszero.(n))

  # invalid number of coordinates for observations
  X = rand(3,2); z = [1, 2, 3]
  @test_throws AssertionError EmpiricalVariogram(X, z)

  # directional variogram and known anisotropy ratio
  imgdata = readdlm(joinpath(datadir,"anisotropic.tsv"))
  geodata = RegularGridData{Float64}(Dict(:z => imgdata))
  γhor = DirectionalVariogram(geodata, (1.,0.), :z, maxlag=50.)
  γver = DirectionalVariogram(geodata, (0.,1.), :z, maxlag=50.)
  γₕ = Variography.fit(GaussianVariogram, γhor)
  γᵥ = Variography.fit(GaussianVariogram, γver)
  @test range(γₕ) / range(γᵥ) ≈ 3. atol=.1

  if visualtests
    TI = training_image("WalkerLake")[1:20,1:20,1]
    xwalker = Float64[i for i=1:20 for j=1:20]
    ywalker = Float64[j for i=1:20 for j=1:20]
    zwalker = Float64[TI[i,j] for i=1:20 for j=1:20]
    γwalker = EmpiricalVariogram(hcat(xwalker,ywalker)', zwalker, maxlag=15.)
    @plottest plot(γwalker) joinpath(datadir,"EmpiricalVariograms.png") !istravis

    @plottest begin
      p1 = plot(γhor, showbins=false, label="horizontal")
      plot!(γver, showbins=false, label="vertical")
      p2 = plot(γₕ, maxlag=50., label="horizontal")
      plot!(γᵥ, maxlag=50., label="vertical")
      plot(p1, p2, layout=(2,1))
    end joinpath(datadir,"DirectionalVariograms.png") !istravis
  end
end


text:
function getD(k::Int,n::Int,x::AbstractArray{T}=[1.:1.:n;]) where T<:Real
    D = spdiagm(n-1,n,0=>fill(-1.,n-1),1=>fill(1.,n-1))
    if k>= 1
        for i in 1:k
            D1 = spdiagm(n-i-1,n-i,0=>fill(-1.,n-i-1),1=>fill(1.,n-i-1))
            diffx = x[(i+1):n]-x[1:(n-i)]
            xdiag = spdiagm(0=> i ./ diffx)
            D = D1*xdiag*D
        end
    end
    D
end

function getTM(k::Int,n::Int,x::Vector{T}=[1.:1.:n;]) where T<:Real
    if (k==1) & (n<=200)
        D = getD(k,n,x)
        return [spdiagm(k+1,n,0=>fill(1.,k+1));D]
    else
        D = getD(k-1,n,x)
        diffx = x[(k+1):n]-x[1:(n-k)]
        xdiag = spdiagm(0=> k ./ diffx)
        return [spdiagm(k,n,0=> fill(1.,k));xdiag*D]
    end
end

function forward_solve!(β::AbstractArray{T},bw::Int,TM::SparseMatrixCSC{T,Int},θ::AbstractArray{T}) where T<:Real
    n = length(θ)
    copyto!(β,θ)
    @inbounds for i in 1:n
        for j in 1:min(i-1,bw-1)
            β[i] -= TM[i,i-j]*β[i-j]
        end
        β[i] /= TM[i,i]
    end
    return zero(T)
end

function back_solve!(∇β::AbstractArray{T},bw::Int,TMᵀ::SparseMatrixCSC{T,Int},x::AbstractArray{T}) where T<:Real
    n = length(x)
    copyto!(∇β,x)
    @inbounds for i in n:-1:1
        for j in 1:min(n-i,bw-1)
            ∇β[i] -= TMᵀ[i,i+j]*∇β[i+j]
        end
        ∇β[i] /= TMᵀ[i,i]
    end
    return zero(T)
end

function thin(xgrid::AbstractArray{T},w::AbstractArray{T},ybar::AbstractArray{T},sse::T,nbins::Int) where T<:Real
    xmin = minimum(xgrid)
    range = maximum(xgrid)-xmin
    gap = range/nbins
    w_new = zeros(T,nbins)
    xgrid_new = zeros(T,nbins)
    ybar_new = zeros(T,nbins)
    sse_new = 0.0
    for i in 1:nbins
        subset = nothing
        if i==1
            subset = xmin .<= xgrid .<= xmin+gap
        else
            subset = xmin+gap*(i-1) .< xgrid .<= xmin+i*gap
        end
        w_new[i] = sum(w[subset])
        if w_new[i]>0
            xgrid_new[i] = sum(w[subset].*xgrid[subset])/w_new[i]
            ybar_new[i] = sum(w[subset].*ybar[subset])/w_new[i]
            sse_new += sum(w[subset].* abs2.(ybar[subset].-ybar_new[i]))
        else
            xgrid_new[i] = xmin+gap*(i-0.5)
            ybar_new[i] = mean(y)
        end
    end
    sse_new += sse
    return xgrid_new,w_new,ybar_new,sse_new
end


text:
using ResampleAndFit, Test, DataFrames, Dates
import FileTools
using Test

tests = ["filtdata.jl",
		 "resampledata.jl",
		 "interpdata.jl",
		 "fitdata.jl",
		 "fillnans.jl",
		 "correctdata.jl",
		 "mergedata.jl"];
# Run all tests in the list
for i in tests
	include(i)
end


text:
dimension(sig) = length(sig)

basis_labels(sig::Tuple) = 1:dimension(sig)
basis_labels(sig::NamedTuple{names}) where names = names

basis_blade_label(sig::Tuple, indices, labels=basis_labels(sig)) = "v$(join(labels[indices]))"
basis_blade_label(sig::NamedTuple, indices, labels=nothing) = join(basis_labels(sig)[indices])

abstract type MetricSignature end

"""
	OffsetSignature(sig, indices)

Metric signature `sig` with offset `indices`, enabling non-standard indexing for multivectors.

Example
===
```jldoctest
julia> lorentzian = OffsetSignature((-1,1,1,1), 0:3) # zero-based indexing
⟨-+++⟩ with indices 0:3
(pretty-printed OffsetSignature{(-1, 1, 1, 1), 0:3}())

julia> (1:4)'basis(lorentzian) # construct spacetime 4-vector
Grade-1 Multivector{⟨-+++⟩ with indices 0:3, 1, Vector{Int64}}:
 1 v0
 2 v1
 3 v2
 4 v3

julia> ans[0] # get first ("time") component
1
```
"""
struct OffsetSignature{sig,indices} <: MetricSignature end
OffsetSignature(sig,indices) = OffsetSignature{sig,indices}()
Base.getindex(::OffsetSignature{sig}, args...) where sig = getindex(sig, args...)
Base.length(::OffsetSignature{sig}) where sig = length(sig)
basis_labels(::OffsetSignature{sig}) where sig = basis_labels(sig)
basis_blade_label(::OffsetSignature{sig,labels}, indices) where {sig,labels} = basis_blade_label(sig, indices, labels)

"""
Convert human-readable basis vector index or symbol into 1-based integer index
"""
normalize_bv_index(sig, i) = i
normalize_bv_index(sig, i::Symbol) = findfirst(==(i), basis_labels(sig))
function normalize_bv_index(sig::Union{Tuple,NamedTuple}, i::Integer)
	r = 1:dimension(sig)
	i ∈ r ? i : error("index $i is outside range $r for signature $sig")
end
function normalize_bv_index(osig::OffsetSignature{sig,indices}, ioffset::Integer) where {sig,indices}
	i = findfirst(==(ioffset), indices)
	isnothing(i) && error("index $ioffset is outside range $indices for $osig")
	i
end

"""
Pretty-print metric signature in short, non-parseable form such as ⟨-+++⟩
"""
show_signature(sig) = repr(sig) # fallback
show_signature(sig::Tuple) = "⟨$(join(map(s -> get(Dict(+1=>"+", -1=>"-"), s, s), sig)))⟩"
function show_signature(sig::NamedTuple)
	items = map(zip(keys(sig), sig)) do (label, square)
		s = get(Dict(+1=>"+", -1=>"-"), square, square)
		"$label$s"
	end
	"⟨$(join(items, ","))⟩"
end
function show_signature(::OffsetSignature{sig,indices}) where {sig,indices}
	"$(show_signature(sig)) with indices $indices"
end


text:
module ITensorChainRules

using ChainRulesCore
using ..ITensors

include("zygoterules.jl")

ITensors.dag(z::AbstractZero) = z

function ChainRulesCore.rrule(::typeof(getindex), x::ITensor, I...)
  y = getindex(x, I...)
  function getindex_pullback(ȳ)
    # TODO: add definition `ITensor(::Tuple{}) = ITensor()`
    # to ITensors.jl so no splatting is needed here.
    x̄ = ITensor(inds(x)...)
    x̄[I...] = unthunk(ȳ)
    Ī = broadcast_notangent(I)
    return (NoTangent(), x̄, Ī...)
  end
  return y, getindex_pullback
end

# Specialized version in order to avoid call to `setindex!`
# within the pullback, should be better for taking higher order
# derivatives in Zygote.
function ChainRulesCore.rrule(::typeof(getindex), x::ITensor)
  y = x[]
  function getindex_pullback(ȳ)
    x̄ = ITensor(unthunk(ȳ))
    return (NoTangent(), x̄)
  end
  return y, getindex_pullback
end

function setinds_pullback(ȳ, x, a...)
  x̄ = ITensors.setinds(ȳ, inds(x))
  ā = broadcast_notangent(a)
  return (NoTangent(), x̄, ā...)
end

function inv_op(f::Function, args...; kwargs...)
  return error(
    "Trying to differentiate `$f` but the inverse of the operation (`inv_op`) `$f` with arguments $args and keyword arguments $kwargs is not defined.",
  )
end

function inv_op(::typeof(prime), x, n::Integer=1; kwargs...)
  return prime(x, -n; kwargs...)
end

function inv_op(::typeof(replaceprime), x, n1n2::Pair; kwargs...)
  return replaceprime(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(swapprime), x, n1n2::Pair; kwargs...)
  return swapprime(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(addtags), x, args...; kwargs...)
  return removetags(x, args...; kwargs...)
end

function inv_op(::typeof(removetags), x, args...; kwargs...)
  return addtags(x, args...; kwargs...)
end

function inv_op(::typeof(replacetags), x, n1n2::Pair; kwargs...)
  return replacetags(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(swaptags), x, n1n2::Pair; kwargs...)
  return swaptags(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(replaceind), x, n1n2::Pair; kwargs...)
  return replaceind(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(replaceinds), x, n1n2::Pair; kwargs...)
  return replaceinds(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(swapind), x, args...; kwargs...)
  return swapind(x, reverse(args)...; kwargs...)
end

function inv_op(::typeof(swapinds), x, args...; kwargs...)
  return swapinds(x, reverse(args)...; kwargs...)
end

for fname in (
  :prime,
  :setprime,
  :noprime,
  :replaceprime,
  :swapprime,
  :addtags,
  :removetags,
  :replacetags,
  :settags,
  :swaptags,
  :replaceind,
  :replaceinds,
  :swapind,
  :swapinds,
)
  @eval begin
    function ChainRulesCore.rrule(
      f::typeof($fname), x::Union{ITensor,MPS,MPO}, a...; kwargs...
    )
      y = f(x, a...; kwargs...)
      function f_pullback(ȳ)
        x̄ = inv_op(f, unthunk(ȳ), a...; kwargs...)
        if !hassameinds(x, x̄)
          error(
            "Trying to differentiate function `$f` with arguments $a and keyword arguments $kwargs. The forward pass indices $(inds(x)) do not match the reverse pass indices $(inds(x̄)). Likely this is because the priming/tagging operation you tried to perform is not invertible. Please write your code in a way where the index manipulation operation you are performing is invertible. For example, `prime(A::ITensor)` is invertible, with an inverse `prime(A, -1)`. However, `noprime(A)` is in general not invertible since the information about the prime levels of the original tensor are lost. Instead, you might try `prime(A, -1)` or `replaceprime(A, 1 => 0)` which are invertible.",
          )
        end
        ā = broadcast_notangent(a)
        return (NoTangent(), x̄, ā...)
      end
      return y, f_pullback
    end
  end
end

# TODO: This is not being called by Zygote for some reason,
# using a Zygote overload directly instead. Figure out
# why, maybe raise an issue.
#function ChainRulesCore.rrule(::typeof(adjoint), x::ITensor)
#  y = prime(x)
#  function adjoint_pullback(ȳ)
#    return setinds_pullback(ȳ, x)
#  end
#  return y, adjoint_pullback
#end

# Special case for contracting a pair of ITensors
function ChainRulesCore.rrule(::typeof(*), x1::ITensor, x2::ITensor)
  y = x1 * x2
  function contract_pullback(ȳ)
    x̄1 = ȳ * dag(x2)
    x̄2 = dag(x1) * ȳ
    return (NoTangent(), x̄1, x̄2)
  end
  return y, contract_pullback
end

function ChainRulesCore.rrule(::typeof(*), x1::Number, x2::ITensor)
  y = x1 * x2
  function contract_pullback(ȳ)
    x̄1 = ȳ * dag(x2)
    x̄2 = dag(x1) * ȳ
    return (NoTangent(), x̄1[], x̄2)
  end
  return y, contract_pullback
end

function ChainRulesCore.rrule(::typeof(*), x1::ITensor, x2::Number)
  y = x1 * x2
  function contract_pullback(ȳ)
    x̄1 = ȳ * dag(x2)
    x̄2 = dag(x1) * ȳ
    return (NoTangent(), x̄1, x̄2[])
  end
  return y, contract_pullback
end

# TODO: use some contraction sequence optimization here
function ChainRulesCore.rrule(::typeof(*), x1::ITensor, x2::ITensor, xs::ITensor...)
  y = *(x1, x2, xs...)
  function contract_pullback(ȳ)
    tn = [x1, x2, xs...]
    N = length(tn)
    env_contracted = Vector{ITensor}(undef, N)
    for n in 1:length(tn)
      tn_left = tn[1:(n - 1)]
      # TODO: define contract([]) = ITensor(1.0)
      env_left = isempty(tn_left) ? ITensor(1.0) : contract(tn_left)
      tn_right = tn[reverse((n + 1):end)]
      env_right = isempty(tn_right) ? ITensor(1.0) : contract(tn_right)
      env_contracted[n] = dag(env_left) * ȳ * dag(env_right)
    end
    return (NoTangent(), env_contracted...)
  end
  return y, contract_pullback
end

function ChainRulesCore.rrule(::typeof(+), x1::ITensor, x2::ITensor)
  y = x1 + x2
  function add_pullback(ȳ)
    return (NoTangent(), ȳ, ȳ)
  end
  return y, add_pullback
end

function ChainRulesCore.rrule(::typeof(itensor), x::Array, a...)
  y = itensor(x, a...)
  function itensor_pullback(ȳ)
    uȳ = permute(unthunk(ȳ), a...)
    x̄ = reshape(array(uȳ), size(x))
    ā = broadcast_notangent(a)
    return (NoTangent(), x̄, ā...)
  end
  return y, itensor_pullback
end

function ChainRulesCore.rrule(::Type{ITensor}, x::Array{<:Number}, a...)
  y = ITensor(x, a...)
  function ITensor_pullback(ȳ)
    # TODO: define `Array(::ITensor)` directly
    uȳ = Array(unthunk(ȳ), a...)
    x̄ = reshape(uȳ, size(x))
    ā = broadcast_notangent(a)
    return (NoTangent(), x̄, ā...)
  end
  return y, ITensor_pullback
end

function ChainRulesCore.rrule(::Type{ITensor}, x::Number)
  y = ITensor(x)
  function ITensor_pullback(ȳ)
    x̄ = ȳ[]
    return (NoTangent(), x̄)
  end
  return y, ITensor_pullback
end

function ChainRulesCore.rrule(::typeof(dag), x)
  y = dag(x)
  function dag_pullback(ȳ)
    x̄ = dag(unthunk(ȳ))
    return (NoTangent(), x̄)
  end
  return y, dag_pullback
end

function ChainRulesCore.rrule(::typeof(permute), x::ITensor, a...)
  y = permute(x, a...)
  function permute_pullback(ȳ)
    x̄ = permute(unthunk(ȳ), inds(x))
    ā = broadcast_notangent(a)
    return (NoTangent(), x̄, ā...)
  end
  return y, permute_pullback
end

broadcast_notangent(a) = broadcast(_ -> NoTangent(), a)

@non_differentiable broadcast_notangent(::Any)
@non_differentiable Index(::Any...)
@non_differentiable delta(::Any...)
@non_differentiable dag(::Index)
@non_differentiable inds(::Any...)
@non_differentiable commoninds(::Any...)
@non_differentiable noncommoninds(::Any...)
@non_differentiable uniqueinds(::Any...)
@non_differentiable SiteType(::Any)
@non_differentiable ITensors._sitetypes(::Any)
@non_differentiable addtags(::TagSet, ::Any)

#
# MPO/MPS
#

# TODO: Define a more general version in ITensors.jl
function _contract(::Type{ITensor}, ψ::MPS, ϕ::MPS; kwargs...)
  T = ITensor(1)
  for n in 1:length(ψ)
    T = T * ψ[n] * ϕ[n]
  end
  return T
end

function _contract(::Type{MPO}, ψ::MPS, ϕ::MPS; kwargs...)
  ψmat = convert(MPO, ψ)
  ϕmat = convert(MPO, ϕ)
  return contract(ψmat, ϕmat; kwargs...)
end

function ChainRulesCore.rrule(::typeof(apply), x1::Vector{ITensor}, x2::MPS; kwargs...)
  y = apply(x1, x2; kwargs...)
  function apply_pullback(ȳ)
    N = length(x1) + 1

    # Apply circuit and store intermediates in the forward direction
    x1x2 = Vector{MPS}(undef, N)
    x1x2[1] = x2
    for n in 2:N
      x1x2[n] = apply(x1[n - 1], x1x2[n - 1]; move_sites_back=false)
    end
    x1x2dag = dag.(x1x2)

    # Apply circuit and store intermediates in the reverse direction

    # XXX: Which one is correct?
    # This works to optimize "Ry" but not "Rx"
    #x1dag = [swapprime(x, 0 => 1) for x in x1]

    # This fails to optimize "Ry" and "Rx"
    #x1dag = [dag(x) for x in x1]

    x1dag = [swapprime(dag(x), 0 => 1) for x in x1]

    x1dag_ȳ = Vector{MPS}(undef, N)
    x1dag_ȳ[end] = ȳ
    for n in (N - 1):-1:1
      x1dag_ȳ[n] = apply(x1dag[n], x1dag_ȳ[n + 1]; kwargs...)
    end

    x̄1 = similar(x1)
    for n in 1:length(x1)
      x1dag_ȳ′ = prime(x1dag_ȳ[n + 1], inds(x1[n]; plev=0))
      x̄1[n] = _contract(ITensor, x1dag_ȳ′, x1x2dag[n]; kwargs...)
    end
    x̄2 = x1dag_ȳ[end]

    return (NoTangent(), x̄1, x̄2)
  end
  return y, apply_pullback
end

function ChainRulesCore.rrule(::typeof(inner), x1::MPS, x2::MPO, x3::MPS; kwargs...)
  y = inner(x1, x2, x3; kwargs...)
  function inner_pullback(ȳ)
    x1dag = dag(x1)
    x̄1 = ȳ * contract(x2, x3; kwargs...)
    x̄2 = ȳ * dag(_contract(MPO, x1dag, x3; kwargs...))
    x̄3 = ȳ * dag(contract(x2, x1dag; kwargs...))

    @assert siteinds(x1) == siteinds(x̄1)
    @assert hassameinds(siteinds, x2, x̄2)
    @assert siteinds(x3) == siteinds(x̄3)

    return (NoTangent(), x̄1, x̄2, x̄3)
  end
  return y, inner_pullback
end

function ChainRulesCore.rrule(::typeof(inner), x1::MPS, x2::MPS; kwargs...)
  y = inner(x1, x2)
  function inner_pullback(ȳ)
    x̄1 = ȳ * dag(x2)
    # `dag` of `x1` gets reversed by `inner`
    x̄2 = x1 * ȳ
    return (NoTangent(), x̄1, x̄2)
  end
  return y, inner_pullback
end

end