text:
#include("./src/SBMLCases.jl")
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using SBMLCases, DataStructures, JSON, Dates

cases_dict = JSON.parsefile("./cases.json"; dicttype = OrderedDict)
build_dict = JSON.parsefile("./build.log"; dicttype = OrderedDict)

### run all cases

required_time = @elapsed run_and_update_status!(
    cases_dict;
    build_dict = build_dict,
    range = 1:1780
)

#heta_version = ENV["heta_version"]
date = Dates.format(Dates.now(), "yyyy-mm-dd HH:MM:SS")
heta_version = length(ARGS) > 0 ? ARGS[1] : "unknown"
pkg = Pkg.installed()
report = Dict(
    "cases" => cases_dict,
    "required_time" => required_time,
    "date" => date,
    "heta_version" => heta_version,
    "solver" => "SimSolver"
)
if haskey(pkg, "SimSolver")
    report["solver_version"] = string(pkg["SimSolver"])
end

open("./results.json", "w") do f
    JSON.print(f, report, 4)
end

### Draft
#=
# get report of builds
report = case_build_errors(
    cases_dict["01112"]; # 01112 00025
    build_dict = build_dict
)

# get report of simulation
report = case_sim_result(
    cases_dict["01001"];
    cases_path = "./cases/semantic", # path to models
    output_path = "./cases/output" # path to meta files
)
=#

max_stars_repo_name: insysbio/sbml-heta-cases

max_issues_repo_name: insysbio/sbml-heta-cases

max_forks_repo_name: insysbio/sbml-heta-cases


text:
function locinfo(frame::Frame)
    if frame.framecode.scope isa Method
        meth = frame.framecode.scope
        ret = JuliaInterpreter.whereis(meth)
        ret === nothing && return nothing
        deffile, _ = ret
        ret = JuliaInterpreter.whereis(frame)
        ret === nothing && return nothing
        current_file, current_line = ret
        local body, defline
        try # https://github.com/timholy/CodeTracking.jl/issues/31
            body, defline = CodeTracking.definition(String, meth)
        catch
            return nothing
        end
        if deffile != current_file || defline > current_line
            isfile(current_file) || return nothing
            body = read(current_file, String)
            defline = 0 # We are not sure where the context start in cases like these, could be improved?
        end
        return defline, current_line, body
    else
        println("not yet implemented")
    end
end

# Used for the tests
const _print_full_path = Ref(true)

function locdesc(frame::Frame; current_line=false)
    sprint() do io
        if frame.framecode.scope isa Method
            locdesc(io, frame; current_line=current_line)
        else
            println(io, "not yet implemented")
        end
    end
end

function locdesc(io, frame::Frame; current_line=false)
    framecode = frame.framecode
    meth = framecode.scope
    @assert meth isa Method
    argnames = framecode.src.slotnames[2:meth.nargs]
    spectypes = Any[Any for i=1:length(argnames)]
    print(io, meth.name,'(')
    first = true
    for (argname, argT) in zip(argnames, spectypes)
        first || print(io, ", ")
        first = false
        print(io, argname)
        !(argT === Any) && print(io, "::", argT)
    end
    line = current_line ? JuliaInterpreter.linenumber(frame) : meth.line
    path = string(_print_full_path[] ? meth.file : basename(String(meth.file)), ":", line)
    path = CodeTracking.replace_buildbot_stdlibpath(String(path))
    print(io, ") at ", path)
end

ext: jl


text:
using Random,Distributions,StatsBase,Plots,HypothesisTests,Measures; pyplot()
Random.seed!(3)

dist = Gamma(2, 2.5)
distH0 = Exponential(5)
n = 100
data = rand(dist,n)

Fhat = ecdf(data)
diffF(dist, x) = sqrt(n)*(Fhat(x) - cdf.(dist,x))
xGrid = 0:0.001:30
ksStat = maximum(abs.(diffF(distH0, xGrid)))

M = 10^5
KScdf(x) = sqrt(2pi)/x*sum([exp(-(2k-1)^2*pi^2 ./(8x.^2)) for k in 1:M])

println("p-value calculated via series: ",
	1-KScdf(ksStat))
println("p-value calculated via Kolmogorov distribution: ",
	1-cdf(Kolmogorov(),ksStat),"\n")

println(ApproximateOneSampleKSTest(data,distH0))

p1 = plot(xGrid, Fhat(xGrid),
	c=:black, lw=1, label="ECDF from data")
p1 = plot!(xGrid, cdf.(dist,xGrid),
	c=:blue, ls=:dot, label="CDF under \n actual distribution")
p1 = plot!(xGrid, cdf.(distH0,xGrid),
	c=:red, ls=:dot, label="CDF under \n postulated H0",
	xlims=(0,20), ylims=(0,1), xlabel = "x", ylabel = "Probability")

p2= plot(cdf.(dist,xGrid), diffF(dist, xGrid),lw=0.5,
	c=:blue,	label="KS Process under \n actual distribution")
p2 = plot!(cdf.(distH0,xGrid), diffF(distH0, xGrid), lw=0.5,
	c=:red, xlims=(0,1), label="KS Process under \n postulated H0",
    xlabel = "t", ylabel = "K-S Process")

plot(p1, p2, legend=:bottomright, size=(800, 400), margin = 5mm)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
@testset "Empty anim" begin
    anim = @animate for i in []
    end
    @test_throws ArgumentError gif(anim)
end

@userplot CirclePlot
@recipe function f(cp::CirclePlot)
    x, y, i = cp.args
    n = length(x)
    inds = circshift(1:n, 1 - i)
    linewidth --> range(0, 10, length = n)
    seriesalpha --> range(0, 1, length = n)
    aspect_ratio --> 1
    label --> false
    x[inds], y[inds]
end

@testset "Circle plot" begin
    n = 10
    t = range(0, 2π, length = n)
    x = sin.(t)
    y = cos.(t)

    anim = @animate for i in 1:n
        circleplot(x, y, i)
    end
    @test filesize(gif(anim).filename) > 10_000
    @test filesize(mov(anim).filename) > 10_000
    @test filesize(mp4(anim).filename) > 10_000
    @test filesize(webm(anim).filename) > 10_000

    @gif for i in 1:n
        circleplot(x, y, i, line_z = 1:n, cbar = false, framestyle = :zerolines)
    end every 5
end

@testset "html" begin
    p = plot([sin, cos], zeros(0), leg = false, xlims = (0, 2π), ylims = (-1, 1))
    anim = Animation()
    for x in range(0, stop = 2π, length = 10)
        push!(p, x, Float64[sin(x), cos(x)])
        frame(anim)
    end

    agif = gif(anim)
    html = tempname() * ".html"
    open(html, "w") do io
        show(io, MIME("text/html"), agif)
    end
    @test filesize(html) > 10_000
    @test showable(MIME("image/gif"), agif)

    agif = mp4(anim)
    html = tempname() * ".html"
    open(html, "w") do io
        show(io, MIME("text/html"), agif)
    end
    @test filesize(html) > 10_000
end


text:
module LapSolve
##### PLEASE READ THE NOTICE ON THE BOTTOM, to know the original authors of this code
using SparseArrays
export solve_lap, solve_stiff_lap
include("StiffWrapper.jl")
# Zero markers used in hungarian algorithm
# 0 => NON   => Non-zero
# 1 => Z     => ordinary zero
# 2 => STAR  => starred zero
# 3 => PRIME => primed zero
const NON = Int8(0)
const Z = Int8(1)
const STAR = Int8(2)
const PRIME = Int8(3)

"""
solve_lap(cost_matrix) -> (assignment, cost)

# Example
```julia
julia> using LapSolve

julia> M=rand(1:100,3,3)
3×3 Array{Int64,2}:
 66  83  92
 17  92  65
 30  84  18

julia> solve_lap(M)
([2, 1, 3], 118)
```
"""
function solve_lap(costMat::AbstractMatrix)
    rowNum, colNum = size(costMat)
    # currently, the function `solve_lap` automatically transposes `cost matrix` when there are more workers than jobs.
    costMatrix = rowNum <= colNum ? costMat : transpose(costMat)
    matching = build_matching(costMatrix)
    assignment = zeros(Int, rowNum)
    rows = rowvals(matching)
    for c = 1:size(matching,2), i in nzrange(matching, c)
        r = rows[i]
        if matching[r,c] == STAR
            if rowNum ≤ colNum
                assignment[r] = c
            else
                assignment[c] = r
            end
        end
    end
    # calculate minimum cost
    cost = sum(costMat[i...] for i in zip(1:rowNum, assignment) if i[2] != 0)
    return assignment,cost
end

"""
solve_stiff_lap(cost_matrix) -> assignment in cartesian coords

# Example
```julia
julia> using LapSolve

julia> M=rand([1,2,Inf],5,5)
5×5 Array{Float64,2}:
 Inf    Inf      1.0    1.0    1.0
 Inf    Inf    Inf    Inf    Inf
   2.0    1.0  Inf      1.0    2.0
 Inf      2.0    1.0    2.0  Inf
   2.0    2.0    2.0  Inf      1.0

julia> solve_stiff_lap(M)
6-element Array{Tuple{Int64,Int64},1}:
 (1, 4)        ← row 1 is assigned to column 4
 (2, -1)       ← -1 for missing assignment, row 2 is assigned to nothing
 (3, 2)
 (4, 3)
 (5, 5)
 (-1, 1)       ← column 1 is assigned to nothing
```
"""
function solve_stiff_lap(C::AbstractMatrix{T}, penalty=1.05) where T <: Real
    m,n=size(C)
    assignment,cost=solve_lap(StiffWrapper(C,penalty))
    c=[(i,assignment[i]) for i in 1:(m+n)]
    for k in eachindex(c)
        i,j=c[k]
        i<=m||(i=-1)
        j<=n||(j=-1)
        c[k]=(i,j)
    end
    [a for a in c if a[1]+a[2]>-2]
end
const suggeststiffsolver="in the cost matrix contains only Inf, unsolvable problem via solve_lap, try solve_stiff_lap"
function build_matching(costMat::AbstractMatrix{T}) where T <: Real
    rowNum, colNum = size(costMat)
    colNum >= rowNum || throw(ArgumentError("Non-square matrix should have more columns than rows."))

    # preliminaries:
    # "no lines are covered;"
    rowCovered = falses(rowNum)
    colCovered = falses(colNum)

    # for tracking indices
    rowCoveredIdx = Int[]
    colCoveredIdx = Int[]
    rowUncoveredIdx = Int[]
    colUncoveredIdx = Int[]
    sizehint!(rowCoveredIdx, rowNum)
    sizehint!(colCoveredIdx, colNum)
    sizehint!(rowUncoveredIdx, rowNum)
    sizehint!(colUncoveredIdx, colNum)

    # "no zeros are starred or primed."
    # use a sparse matrix Zs to store these three kinds of zeros:
    # 0 => NON   => Non-zero
    # 1 => Z     => ordinary zero
    # 2 => STAR  => starred zero
    # 3 => PRIME => primed zero
    Zs = spzeros(Int8, rowNum, colNum)

    # for tracking changes per row/col of A
    Δrow = zeros(T,rowNum)
    Δcol = zeros(T,colNum)

    # "subtract" minimum from each row
    @inbounds for i in 1:rowNum
        mw=costMat[i,1]
        @inbounds for j in 2:colNum
            cost=costMat[i,j]
            if cost<mw
                mw=cost
            end
        end
        Δrow[i] = -mw
    end
    for i in Δrow
        isinf(i) && error("row $i $suggeststiffsolver")
    end
    # "subtract" minimum from each column
    if rowNum == colNum
        @inbounds for j in 1:colNum
            mw=costMat[1,j]+Δrow[1]
            @inbounds for i in 2:rowNum
                cost=costMat[i,j]+Δrow[i]
                if cost<mw
                    mw=cost
                end
            end
            Δcol[j] = -mw
        end
    end
    for i in Δcol
        isinf(i) && error("column $i $suggeststiffsolver")
    end

    # for tracking those starred zero
    rowSTAR = falses(rowNum)
    colSTAR = falses(colNum)
    # since looping through a row in a SparseMatrixCSC is costy(not cache-friendly),
    # a row to column index mapping of starred zeros is also tracked here.
    row2colSTAR = Dict{Int,Int}()
    for ii in CartesianIndices(costMat)
        r,c=Tuple(ii)
        cost=costMat[r,c]+Δrow[r]+Δcol[c]
        # "consider a zero Z of the matrix;"
        if iszero(cost)
            Zs[r,c] = Z
            # "if there is no starred zero in its row and none in its column, star Z.
            #  repeat, considering each zero in the matrix in turn;"
            if !colSTAR[c] && !rowSTAR[r]
                Zs[r,c] = STAR
                rowSTAR[r] = true
                colSTAR[c] = true
                row2colSTAR[r] = c
                # "then cover every column containing a starred zero."
                colCovered[c] = true
            end
        end
    end

    # preliminaries done, go to step 1
    stepNum = 1

    # if the assignment is already found, exit
    stepNum = exit_criteria(colCovered, size(Zs))

    # pre-allocation
    sequence = Tuple{Int,Int}[]       # used in step 2
    minLocations = Tuple{Int,Int}[]   # used in step 3

    # these three steps are parallel with those in the paper:
    # J. Munkres, "Algorithms for the Assignment and Transportation Problems",
    # Journal of the Society for Industrial and Applied Mathematics, 5(1):32–38, 1957 March.
    while stepNum != 0
        if stepNum == 1
            stepNum = step1!(Zs, rowCovered, colCovered, rowSTAR, row2colSTAR)
        elseif stepNum == 2
            empty!(sequence)
            stepNum = step2!(Zs, sequence, rowCovered, colCovered, rowSTAR, row2colSTAR)
        elseif stepNum == 3
            empty!(rowCoveredIdx)
            empty!(colCoveredIdx)
            empty!(rowUncoveredIdx)
            empty!(colUncoveredIdx)
            empty!(minLocations)
            stepNum = step3!(costMat, Zs, minLocations, rowCovered, colCovered, rowSTAR, row2colSTAR, Δrow, Δcol, rowCoveredIdx, colCoveredIdx, rowUncoveredIdx, colUncoveredIdx)
        end
    end

    return Zs
end

"""
    exit_criteria(colCovered, ZsDims) -> stepNum
We adjust Munkres's algorithm in order to deal with rectangular matrices,
so only K columns are counted here, where K = min(size(Zs))
"""
function exit_criteria(colCovered, ZsDims)
    count = 0
    @inbounds for i in eachindex(colCovered)
        colCovered[i] && (count += 1;)
    end
    if count == ZsDims[1]
        # algorithm exits
        return 0
    else
        # "otherwise, return to step 1."
        return 1
    end
end

"""
Step 1 of the original Munkres' Assignment Algorithm
"""
function step1!(Zs, rowCovered, colCovered, rowSTAR, row2colSTAR)
    colLen = size(Zs,2)
    rows = rowvals(Zs)
    # step 1:
    zeroCoveredNum = 0
    # "repeat until all zeros are covered."
    while zeroCoveredNum < nnz(Zs)
        zeroCoveredNum = 0
        for c = 1:colLen, i in nzrange(Zs, c)
            r = rows[i]
            # "choose a non-covered zero and prime it"
            if colCovered[c] == false && rowCovered[r] == false
                Zs[r,c] = PRIME
                # "consider the row containing it."
                # "if there is a starred zero Z in this row"
                if rowSTAR[r]
                    # "cover this row and uncover the column of Z"
                    rowCovered[r] = true
                    colCovered[row2colSTAR[r]] = false
                else
                    # "if there is no starred zero in this row,
                    #  go at once to step 2."
                    return 2
                end
            else
                # otherwise, this zero is covered
                zeroCoveredNum += 1
            end
        end
    end
    # "go to step 3."
    return 3
end

"""
Step 2 of the original Munkres' Assignment Algorithm
"""
function step2!(Zs, sequence, rowCovered, colCovered, rowSTAR, row2colSTAR)
    ZsDims = size(Zs)
    rows = rowvals(Zs)
    # step 2:
    flag = false
    # "there is a sequence of alternating primed and starred zeros, constructed
    #  as follows:"
    # "let Z₀ denote the uncovered 0′.[there is only one.]"
    for c = 1:ZsDims[2], i in nzrange(Zs, c)
        r = rows[i]
        # note that Z₀ is an **uncovered** 0′
        if Zs[r,c] == PRIME && colCovered[c] == false && rowCovered[r] == false
            # push Z₀(uncovered 0′) into the sequence
            push!(sequence, (r, c))
            # "let Z₁ denote the 0* in the Z₀'s column.(if any)"
            # find 0* in Z₀'s column
            for j in nzrange(Zs, c)
                Z₁r = rows[j]
                if Zs[Z₁r, c] == STAR
                    # push Z₁(0*) into the sequence
                    push!(sequence, (Z₁r, c))
                    # set sequence continue flag => true
                    flag = true
                    break
                end
            end
            break
        end
    end
    # "let Z₂ denote the 0′ in Z₁'s row(there will always be one)."
    # "let Z₃ denote the 0* in the Z₂'s column."
    # "continue until the sequence stops at a 0′ Z₂ⱼ, which has no 0* in its column."
    while flag
        flag = false
        r = sequence[end][1]
        # find Z₂ in Z₃'s row (always exits)
        for c = 1:ZsDims[2]
            if Zs[r,c] == PRIME
                # push Z₂ into the sequence
                push!(sequence, (r, c))
                # find 0* in Z₂'s column
                for j in nzrange(Zs, c)
                    Z₃r = rows[j]
                    if Zs[Z₃r, c] == STAR
                        # push Z₃(0*) into the sequence
                        push!(sequence, (Z₃r, c))
                        # set sequence continue flag => true
                        flag = true
                        break
                    end
                end
                break
            end
        end
    end

    # "unstar each starred zero of the sequence;"
    for i in 2:2:length(sequence)-1
        Zs[sequence[i]...] = Z
    end

    # clean up
    fill!(rowSTAR, false)
    empty!(row2colSTAR)

    # "and star each primed zero of the sequence."
    for i in 1:2:length(sequence)
        Zs[sequence[i]...] = STAR
    end

    for c = 1:ZsDims[2], i in nzrange(Zs, c)
        r = rows[i]
        # "erase all primes;"
        if Zs[r,c] == PRIME
            Zs[r,c] = Z
        end
        # "and cover every column containing a 0*."
        if Zs[r,c] == STAR
            colCovered[c] = true
            rowSTAR[r] = true
            row2colSTAR[r] = c
        end
    end

    # "uncover every row"
    fill!(rowCovered, false)

    # "if all columns are covered, the starred zeros form the desired independent set."
    return exit_criteria(colCovered, ZsDims)
end

"""
Step 3 of the original Munkres' Assignment Algorithm
"""
function step3!(costMat::AbstractMatrix{T}, Zs, minLocations, rowCovered, colCovered, rowSTAR, row2colSTAR,
                Δrow, Δcol, rowCoveredIdx, colCoveredIdx, rowUncoveredIdx, colUncoveredIdx) where T <: Real

    @inbounds for i in eachindex(rowCovered)
        rowCovered[i] ? push!(rowCoveredIdx, i) : push!(rowUncoveredIdx, i)
    end

    @inbounds for i in eachindex(colCovered)
        colCovered[i] ? push!(colCoveredIdx, i) : push!(colUncoveredIdx, i)
    end

    h = typemax(T)
    @inbounds for j in colUncoveredIdx, i in rowUncoveredIdx
        cost = costMat[i,j] + Δcol[j] + Δrow[i]
        if cost <= h
            if cost != h
                h = cost
                empty!(minLocations)
            end
            push!(minLocations, (i,j))
        end
    end

    for i in rowCoveredIdx
        Δrow[i] += h
    end

    for i in colUncoveredIdx
        Δcol[i] -= h
    end

    # mark new zeros in Zs and remove those elements that will no longer be zero
    # produce new zeros
    for loc in minLocations
        Zs[loc...] = Z
    end
    # reduce old zeros
    rows = rowvals(Zs)
    for c in colCoveredIdx, i in nzrange(Zs, c)
        r = rows[i]
        if rowCovered[r]
            if Zs[r,c] == STAR
                rowSTAR[r] = false
                delete!(row2colSTAR, r)
            end
            Zs[r,c] = NON
        end
    end

    dropzeros!(Zs)

    # "return to step 1."
    return stepNum = 1
end

#### NOTICE!!!
#### most functions in this file are taken with almost no modifications from https://github.com/Gnimuc/Hungarian.jl
#### credit for those functions go to the original contributors: Yupei Qi <qiyupei@gmail.com>, et al
#### LICENSE is the Unlicense @ https://github.com/Gnimuc/Hungarian.jl/blob/0af029c5bc25bdff8d75296f9ae510e3f236f0aa/LICENSE
#### --------------------------------------------------------------------------------------------------------

end # module

max_stars_repo_name: francescoalemanno/SolveLAP


text:
immutable AdamsBashforth5 <: TimeIntergral_onestep
  rf :: Tuple{Vararg{Array{Float64,1}}}

  function AdamsBashforth5()
    rf1 = [1.0]
    rf2 = [3/2;-1/2]
    rf3 = [23/12; -4/3; 5/12]
    rf4 = [55/24; -59/24; 37/24; -3/8]
    rf5 = [1901/720; -1387/360; 109/30; -637/360; 251/720]
    rf = (rf1, rf2,rf3,rf4,rf5)
    new(rf)
  end
end

function eval!(tf :: AdamsBashforth5, p :: ModelParams,
  data :: Tuple, trends :: Tuple,
  t :: Float64, dt :: Float64)

  new_trend = fTrend(data, p,t)
  all_trends = (new_trend, trends...)
  ab_order = length(all_trends)
  for ii = 1:ab_order
    delta_ii = map(all_trends[ii]) do f
      f* (p.rf[ii] * dt)
    end
    data = map(+, data, delta_ii)
  end
  TimeForward!(p,t,dt)
  return (data,new_trend)
end


text:
"""
FunctionTransform
```
    f(x) = abs.(x)
    tr = FunctionTransform(f)
```
Take a function or object `f` as an argument which is going to act on each vector individually.
Make sure that `f` is supposed to act on a vector by eventually using broadcasting
For example `f(x)=sin(x)` -> `f(x)=sin.(x)`
"""
struct FunctionTransform{F} <: Transform
    f::F
end

apply(t::FunctionTransform, X::T; obsdim::Int = defaultobs) where {T} = mapslices(t.f, X, dims = feature_dim(obsdim))

duplicate(t::FunctionTransform,f) = FunctionTransform(f)


text:
# # Fitting a polynomial with kernels
#
# First, let us load the relevant packages and define our problem.
# We want to interpolate a simple polynomial: `(2x-1)² + 2y + xy - 3`.

using KernelMachines, Statistics, Plots

f(x, y) = (2x-1)^2 + 2y + x * y - 3
f((x, y)) = f(x, y)
rg = 0:0.01:1
N = length(rg)
flat = hcat(repeat(rg, inner=N), repeat(rg, outer=N))
truth = map(f, eachrow(flat))
theme(:wong)
surface(rg, rg, reshape(truth, N, N), clims=(-3, 1))

# Let us generate a `6 x 6` trainig grid.

N_train = 6
rg_train = range(0, 1, length=N_train)
X = hcat(repeat(rg_train, inner=N_train), repeat(rg_train, outer=N_train))
Y = map(f, eachrow(X));

# Now, let us train a Kernel Machine on the problem.

krm = KernelMachineRegression(X, Y;
    dims=(3, 3, 3), kernel=multiplicativegaussiankernel)
fit!(krm)

pred_krm = predict(krm, flat)
surface(rg, rg, reshape(pred_krm, N, N), clims=(-3, 1))

# In this problem, a simple Kernel Ridge regression also performs very well.

kr = KernelRegression(X, Y)
fit!(kr)
pred_kr = predict(kr, flat)
surface(rg, rg, reshape(pred_kr, N, N), clims=(-3, 1))


text:
module Qt5

__precompile__(false)

import Base: +,-,*,/

# using ImmutableArrays

# typealias P2 Vector2{Float64}
# typealias P3 Vector3{Float64}
# typealias Point Union(P2,P3)

# const ORIGIN = P3(0,0,0)

# P2(p::P3) = P2(p[1], p[2])
# P3(p::P2, z::Real) = P3(p[1], p[2], z)
# P3(p::P2) = P3(p, 0.0)

include("measures.jl")
include("core.jl")
include("cppdefs.jl")
include("widgets.jl")

end # module

q = Qt5


text:
#Moore et al 2009 mean vectors:
M=Array{Any,1}(undef,8)
M[1]=vec([0.0234 0.0192 0.0129 0.0075 0.0031 0.0002])
M[2]=vec([0.0162 0.0141 0.0112 0.0073 0.0034 0.0002])
M[3]=vec([0.0107 0.0098 0.0092 0.0070 0.0039 0.0003])
M[4]=vec([0.0065 0.0064 0.0070 0.0064 0.0048 0.0006])
M[5]=vec([0.0033 0.0034 0.0042 0.0042 0.0043 0.0009])
M[6]=vec([0.0064 0.0074 0.0105 0.0116 0.0140 0.0041])
M[7]=vec([0.0121 0.0140 0.0192 0.0204 0.0231 0.0084])
M[8]=vec([0.0184 0.0230 0.0333 0.0359 0.0409 0.0137])

#Moore et al 2009 covariance matrices:
S=Array{Any,1}(undef,8)
S[1]=[0.00000959 0.00000556 0.00000138 -0.00000034 -0.00000024 -0.00000003;
        0.00000556 0.00000493 0.00000193 0.00000060 0.00000023 0.00000001;
        0.00000138 0.00000193 0.00000282 0.00000223 0.00000119 0.00000007;
        -0.00000034 0.00000060 0.00000223 0.00000232 0.00000119 0.00000007;
        -0.00000024 0.00000023 0.00000119 0.00000119 0.00000071 0.00000005;
        -0.00000003 0.00000001 0.00000007 0.00000007 0.00000005 0.00000001]
S[2]=[0.00000346 0.00000186 -0.00000011 -0.00000060 -0.00000062 -0.00000005;
        0.00000186 0.00000228 0.00000086 0.00000033 -0.00000007 0.00000003;
        -0.00000011 0.00000086 0.00000231 0.00000221 0.00000145 0.00000023;
        -0.00000060 0.00000033 0.00000221 0.00000266 0.00000191 0.00000034;
        -0.00000062 -0.00000007 0.00000145 0.00000191 0.00000175 0.00000041;
        -0.00000005 0.00000003 0.00000023 0.00000034 0.00000041 0.00000021]
S[3]=[0.00000241 0.00000144 0.00000035 -0.00000031 -0.00000063 -0.00000006;
        0.00000144 0.00000138 0.00000076 0.00000015 -0.00000021 -0.00000001;
        0.00000035 0.00000076 0.00000161 0.00000156 0.00000121 0.00000016;
        -0.00000031 0.00000015 0.00000156 0.00000227 0.00000209 0.00000031;
        -0.00000063 -0.00000021 0.00000121 0.00000209 0.00000225 0.00000037;
        -0.00000006 -0.00000001 0.00000016 0.00000031 0.00000037 0.00000013]
S[4]=[0.00000166 0.00000091 0.00000034 -0.00000009 -0.00000080 -0.00000025;
        0.00000091 0.00000097 0.00000071 0.00000025 -0.00000041 -0.00000015;
        0.00000034 0.00000071 0.00000118 0.00000103 0.00000072 0.00000003;
        -0.00000009 0.00000025 0.00000103 0.00000137 0.00000162 0.00000025;
        -0.00000080 -0.00000041 0.00000072 0.00000162 0.00000290 0.00000065;
        -0.00000025 -0.00000015 0.00000003 0.00000025 0.00000065 0.00000050]
S[5]=[0.00000178 0.00000132 0.00000104 0.00000081 0.00000018 -0.00000014;
        0.00000132 0.00000127 0.00000121 0.00000099 0.00000034 -0.00000010;
        0.00000104 0.00000121 0.00000150 0.00000142 0.00000110 0.00000013;
        0.00000081 0.00000099 0.00000142 0.00000158 0.00000177 0.00000042;
        0.00000018 0.00000034 0.00000110 0.00000177 0.00000351 0.00000131;
        -0.00000014 -0.00000010 0.00000013 0.00000042 0.00000131 0.00000081]
S[6]=[0.00000715 0.00000586 0.00000409 0.00000292 0.00000005 -0.00000075;
        0.00000586 0.00000589 0.00000520 0.00000398 0.00000027 -0.00000114;
        0.00000409 0.00000520 0.00000634 0.00000541 0.00000188 -0.00000097;
        0.00000292 0.00000398 0.00000541 0.00000528 0.00000392 0.00000070;
        0.00000005 0.00000027 0.00000188 0.00000392 0.00000995 0.00000657;
        -0.00000075 -0.00000114 -0.00000097 0.00000070 0.00000657 0.00000819]
S[7]=[0.00002625 0.00001981 0.00001058 0.00000544 -0.00000654 -0.00001122;
        0.00001981 0.00001745 0.00001314 0.00000822 -0.00000431 -0.00001228;
        0.00001058 0.00001314 0.00001629 0.00001226 0.00000035 -0.00001311;
        0.00000544 0.00000822 0.00001226 0.00001170 0.00000742 -0.00000500;
        -0.00000654 -0.00000431 0.00000035 0.00000742 0.00002241 0.00001782;
        -0.00001122 -0.00001228 -0.00001311 -0.00000500 0.00001782 0.00003987]
S[8]=[0.00001186 0.00001134 0.00001139 0.00000919 0.00000395 -0.00000186;
        0.00001134 0.00001484 0.00002034 0.00001907 0.00001531 0.00000087;
        0.00001139 0.00002034 0.00003467 0.00003546 0.00003555 0.00000708;
        0.00000919 0.00001907 0.00003546 0.00003907 0.00004604 0.00001733;
        0.00000395 0.00001531 0.00003555 0.00004604 0.00007306 0.00004953;
        -0.00000186 0.00000087 0.00000708 0.00001733 0.00004953 0.00006542];

max_stars_repo_name: gaelforget/MarineEcosystemNotebooks

max_forks_repo_name: gaelforget/MarineEcosystemNotebooks


text:
import TensorDecompositions
import NTFk

csize = (2, 3, 4)
tsize = (5, 10, 15)
tucker_orig = NTFk.rand_tucker(csize, tsize; factors_nonneg=true, core_nonneg=true)

T_orig = TensorDecompositions.compose(tucker_orig)
T_orig .*= 1000

sizes = [csize, (1,3,4), (3,3,4), (2,2,4), (2,4,4), (2,3,3), (2,3,5)]
tucker_estimated, csize_estimated = NTFk.analysis(T_orig, sizes, 3; eigmethod=[false,false,false], progressbar=false, tol=1e-16, max_iter=100000, lambda=0.);

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
"""
    HarmonicBond(i, j, b0, kb)

A harmonic bond between two atoms.
"""
struct HarmonicBond{T} <: SpecificInteraction
    i::Int
    j::Int
    b0::T
    kb::T
end

@inline @inbounds function force(b::HarmonicBond,
                coords,
                s::Simulation)
    ab = vector(coords[b.i], coords[b.j], s.box_size)
    c = b.kb * (norm(ab) - b.b0)
    f = c * normalize(ab)
    return [b.i, b.j], [f, -f]
end

@inline @inbounds function potential_energy(b::HarmonicBond,
                                            s::Simulation)
    dr = vector(s.coords[b.i], s.coords[b.j], s.box_size)
    r = norm(dr)
    b.kb / 2 * (r - b.b0) ^ 2
end


text:
using BEAST

function hassemble(
    operator::BEAST.AbstractOperator,
    test_functions,
    trial_functions;
    compressor=:aca,
    tol=1e-4,
    treeoptions=BoxTreeOptions(nmin=100),
    maxrank=100,
    threading=:single,
    quadstrat=BEAST.defaultquadstrat(operator, test_functions, trial_functions),
    verbose=false,
    svdrecompress=true
)

    @views blkasm = BEAST.blockassembler(operator, test_functions, trial_functions)

    @views function assembler(Z, tdata, sdata)
        @views store(v,m,n) = (Z[m,n] += v)
        blkasm(tdata,sdata,store)
    end

    @views farblkasm = BEAST.blockassembler(
        operator,
        test_functions,
        trial_functions,
        quadstrat=quadstrat
    )

    @views function farassembler(Z, tdata, sdata)
        @views store(v,m,n) = (Z[m,n] += v)
        farblkasm(tdata,sdata,store)
    end

    test_tree = create_tree(test_functions.pos, treeoptions)
    trial_tree = create_tree(trial_functions.pos, treeoptions)

    @time hmat = HMatrix(
        assembler,
        test_tree,
        trial_tree,
        compressor=compressor,
        T=scalartype(operator),
        tol=tol,
        maxrank=maxrank,
        threading=threading,
        farmatrixassembler=farassembler,
        verbose=verbose,
        svdrecompress=svdrecompress
    )
    return hmat
end

# The following to function ensure that no dynamic dispatching is
# performed since we know already that all triangles are well-separate

# Copied from BEAST/examples/quadstrat.jl
function BEAST.quaddata(op, tref, bref,
    tels, bels, qs::BEAST.DoubleNumQStrat)

    qs = BEAST.DoubleNumWiltonSauterQStrat(qs.outer_rule, qs.inner_rule, 1, 1, 1, 1, 1, 1)
    BEAST.quaddata(op, tref, bref, tels, bels, qs)
end

# Copied from BEAST/examples/quadstrat.jl
function BEAST.quadrule(op, tref, bref,
    i ,τ, j, σ, qd, qs::BEAST.DoubleNumQStrat)

    return BEAST.DoubleQuadRule(
        qd.tpoints[1,i],
        qd.bpoints[1,j])
end

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: sbadrian/FastBEAST

max_forks_repo_name: sbadrian/FastBEAST


text:
"""
    SemidirectProductOperation(action::AbstractGroupAction)

Group operation of a semidirect product group. The operation consists of the operation
`opN` on a normal subgroup `N`, the operation `opH` on a subgroup `H`, and an automorphism
`action` of elements of `H` on `N`. Only the action is stored.
"""
struct SemidirectProductOperation{A<:AbstractGroupAction} <: AbstractGroupOperation
    action::A
end

function show(io::IO, op::SemidirectProductOperation)
    print(io, "SemidirectProductOperation($(op.action))")
end

const SemidirectProductGroup{N,H,A} =
    GroupManifold{ProductManifold{Tuple{N,H}},SemidirectProductOperation{A}}

@doc doc"""
    SemidirectProductGroup(N::GroupManifold, H::GroupManifold, A::AbstractGroupAction)

A group that is the semidirect product of a normal group $N$ and a subgroup $H$, written
$G = N ⋊_θ H$, where $θ: H × N \to N$ is an automorphism action
of $H$ on $N$. The group $G$ has the composition rule

````math
g \circ g' = (n, h) \circ (n', h') = (n \circ θ_h(n'), h \circ h')
````

and the inverse

````math
g^{-1} = (n, h)^{-1} = (θ_{h^{-1}}(n^{-1}), h^{-1}).
````
"""
function SemidirectProductGroup(N::GroupManifold, H::GroupManifold, A::AbstractGroupAction)
    N === g_manifold(A) || error("Subgroup $(N) must be the G-manifold of action $(A)")
    H === base_group(A) || error("Subgroup $(H) must be the base group of action $(A)")
    op = SemidirectProductOperation(A)
    M = ProductManifold(N, H)
    return GroupManifold(M, op)
end

function show(io::IO, G::SemidirectProductGroup)
    M = base_manifold(G)
    N, H = M.manifolds
    A = G.op.action
    print(io, "SemidirectProductGroup($(N), $(H), $(A))")
end

_padpoint!(G::SemidirectProductGroup, y) = y

_padvector!(G::SemidirectProductGroup, v) = v

inv(G::GT, e::Identity{GT}) where {GT<:SemidirectProductGroup} = e

function inv!(G::SemidirectProductGroup, y, x)
    M = base_manifold(G)
    N, H = M.manifolds
    A = G.op.action
    nx, hx = submanifold_components(G, x)
    ny, hy = submanifold_components(G, y)
    inv!(H, hy, hx)
    ninv = inv(N, nx)
    apply!(A, ny, hy, ninv)
    @inbounds _padpoint!(G, y)
    return y
end
inv!(G::AG, y, e::Identity{AG}) where {AG<:SemidirectProductGroup} = identity!(G, y, e)

identity(G::GT, e::Identity{GT}) where {GT<:SemidirectProductGroup} = e

function identity!(G::SemidirectProductGroup, y, x)
    M = base_manifold(G)
    N, H = M.manifolds
    nx, hx = submanifold_components(G, x)
    ny, hy = submanifold_components(G, y)
    identity!(N, ny, nx)
    identity!(H, hy, hx)
    @inbounds _padpoint!(G, y)
    return y
end
identity!(G::GT, e::E, ::E) where {GT<:SemidirectProductGroup,E<:Identity{GT}} = e

compose(G::GT, x, e::Identity{GT}) where {GT<:SemidirectProductGroup} = x
compose(G::GT, e::Identity{GT}, x) where {GT<:SemidirectProductGroup} = x
compose(G::GT, e::E, ::E) where {GT<:SemidirectProductGroup,E<:Identity{GT}} = e

function compose!(G::SemidirectProductGroup, z, x, y)
    M = base_manifold(G)
    N, H = M.manifolds
    A = G.op.action
    nx, hx = submanifold_components(G, x)
    ny, hy = submanifold_components(G, y)
    nz, hz = submanifold_components(G, z)
    compose!(H, hz, hx, hy)
    zₙtmp = apply(A, hx, ny)
    compose!(N, nz, nx, zₙtmp)
    @inbounds _padpoint!(G, z)
    return z
end
compose!(G::GT, z, ::Identity{GT}, y) where {GT<:SemidirectProductGroup} = copyto!(z, y)
compose!(G::GT, z, x, ::Identity{GT}) where {GT<:SemidirectProductGroup} = copyto!(z, x)
function compose!(G::GT, z, e::E, ::E) where {GT<:SemidirectProductGroup,E<:Identity{GT}}
    return identity!(G, z, e)
end

function translate_diff!(G::SemidirectProductGroup, vout, x, y, v, conv::LeftAction)
    M = base_manifold(G)
    N, H = M.manifolds
    A = G.op.action
    nx, hx = submanifold_components(G, x)
    ny, hy = submanifold_components(G, y)
    nv, hv = submanifold_components(G, v)
    nvout, hvout = submanifold_components(G, vout)
    translate_diff!(H, hvout, hx, hy, hv, conv)
    nw = apply_diff(A, hx, ny, nv)
    nz = apply(A, hx, ny)
    translate_diff!(N, nvout, nx, nz, nw, conv)
    @inbounds _padvector!(G, vout)
    return vout
end

function hat!(G::SemidirectProductGroup, V, x, v)
    M = base_manifold(G)
    N, H = M.manifolds
    dimN = manifold_dimension(N)
    dimH = manifold_dimension(H)
    @assert length(v) == dimN + dimH
    nx, hx = submanifold_components(G, x)
    nV, hV = submanifold_components(G, V)
    hat!(N, nV, nx, view(v, 1:dimN))
    hat!(H, hV, hx, view(v, dimN+1:dimN+dimH))
    @inbounds _padvector!(G, V)
    return V
end

function vee!(G::SemidirectProductGroup, v, x, V)
    M = base_manifold(G)
    N, H = M.manifolds
    dimN = manifold_dimension(N)
    dimH = manifold_dimension(H)
    @assert length(v) == dimN + dimH
    nx, hx = submanifold_components(G, x)
    nV, hV = submanifold_components(G, V)
    vee!(N, view(v, 1:dimN), nx, nV)
    vee!(H, view(v, dimN+1:dimN+dimH), hx, hV)
    return v
end

function zero_tangent_vector(G::SemidirectProductGroup, x)
    v = allocate_result(G, zero_tangent_vector, x)
    zero_tangent_vector!(G, v, x)
    return v
end

function zero_tangent_vector!(G::SemidirectProductGroup, v, x)
    M = base_manifold(G)
    N, H = M.manifolds
    nx, hx = submanifold_components(G, x)
    nv, hv = submanifold_components(G, v)
    zero_tangent_vector!(N, nv, nx)
    zero_tangent_vector!(H, hv, hx)
    return v
end

function isapprox(G::SemidirectProductGroup, x, y; kwargs...)
    M = base_manifold(G)
    N, H = M.manifolds
    nx, hx = submanifold_components(G, x)
    ny, hy = submanifold_components(G, y)
    return isapprox(N, nx, ny; kwargs...) && isapprox(H, hx, hy; kwargs...)
end
function isapprox(G::SemidirectProductGroup, x, v, w; kwargs...)
    M = base_manifold(G)
    N, H = M.manifolds
    nx, hx = submanifold_components(G, x)
    nv, hv = submanifold_components(G, v)
    nw, hw = submanifold_components(G, w)
    return isapprox(N, nx, nv, nw; kwargs...) && isapprox(H, hx, hv, hw; kwargs...)
end
function isapprox(G::GT, x, e::Identity{GT}; kwargs...) where {GT<:SemidirectProductGroup}
    return isapprox(G, e, x; kwargs...)
end
function isapprox(G::GT, e::Identity{GT}, x; kwargs...) where {GT<:SemidirectProductGroup}
    return isapprox(G, identity(G, x), x; kwargs...)
end
function isapprox(
    ::GT,
    ::E,
    ::E;
    kwargs...,
) where {GT<:SemidirectProductGroup,E<:Identity{GT}}
    return true
end


text:
using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct
using ..Units
using Parameters, Unitful

# Computational routines needed for computations with the MaterialParams structure

# with tuple & vector - apply for all phases in MatParam
function compute_param!(
    fn::F,
    rho::AbstractVector,
    MatParam::NTuple{N,AbstractMaterialParamsStruct},
    args,
) where {F,N}
    return rho .= map(x -> fn(x, P, T), MatParam)
end

# each individual calcuation
function compute_param(
    fn::F, MatParam::NTuple{N,AbstractMaterialParamsStruct}, args
) where {F,N}
    return map(x -> fn(x, args), MatParam)
end

#---------------------------------------------------------------------------------------------------------------------------#
#Computational routines for Phases

# performs computation given a single Phase
@inline @generated function compute_param(
    fn::F,MatParam::NTuple{N,AbstractMaterialParamsStruct}, Phase::Int64, args
) where {F, N}
    quote
        T = isempty(args) ? 0.0 : zero(typeof(args).types[1])
        out = T
        Base.Cartesian.@nexprs $N i-> out += MatParam[i].Phase == Phase ? fn(MatParam[i], args) : T
    end
end

function compute_param(
    fn::F, MatParam::AbstractVector{AbstractMaterialParamsStruct}, Phase::Int64, args
) where F
    return compute_param(fn, Tuple(MatParam), Phase, args)
end

@inline function compute_param!(
    fn::F,
    rho::AbstractArray,
    MatParam::NTuple{N,AbstractMaterialParamsStruct},
    Phases::AbstractArray{<:Integer,ndim},
    args,
) where {F,ndim,N}
    @inbounds for I in eachindex(Phases)
        k = keys(args)
        v = getindex.(values(args), I)
        argsi = (; zip(k, v)...)
        rho[I] = compute_param(fn, MatParam, Phases[I], argsi)
    end
end

function compute_param!(
    fn::F,
    rho::AbstractArray,
    MatParam::AbstractVector{AbstractMaterialParamsStruct},
    Phases::AbstractArray{<:Integer,ndim},
    args,
) where {F,ndim}
    return compute_param!(fn, rho, Tuple(MatParam), Phases, args)
end

#Phase PhaseRatios

function compute_param!(
    fn::F,
    rho::AbstractArray{_T,N},
    MatParam::NTuple{K,AbstractMaterialParamsStruct},
    PhaseRatio::AbstractArray{_T,M},
    args,
) where {F,_T,N,M,K}
    if M != (N + 1)
        error("The PhaseRatios array should have one dimension more than the other arrays")
    end

    @inbounds for I in CartesianIndices(rho)
        k = keys(args)
        v = getindex.(values(args), Tuple(I)...)
        argsi = (; zip(k, v)...)
        # unroll the dot product
        val = Ref(zero(_T))
        ntuple(Val(K)) do i
            Base.@_inline_meta
            val[] += PhaseRatio[I.I..., i] * fn(MatParam[i], argsi)
        end
        rho[I] = val[]
    end
end

#Multiplies parameter with the fraction of a phase
function compute_param_times_frac(
    fn::F,
    PhaseRatios::NTuple{N,T},
    MatParam::NTuple{N,AbstractMaterialParamsStruct},
    argsi,
) where {F,N,T}
    # Unrolled dot product
    val = Ref(zero(T))
    ntuple(Val(N)) do i
        Base.@_inline_meta
        val[] += PhaseRatios[i] * fn(MatParam[i], argsi)
    end
    return val[]
end


text:
using LoopVectorization, Test#, OffsetArrays

@testset "Array Wrappers" begin
	@show @__LINE__

    function addone!(y, x)
        @turbo for i ∈ eachindex(x,y)
            y[i] = x[i] + 1
        end
        y
    end

    x5 = PermutedDimsArray(rand(10,10), (2,1));
    y5 = PermutedDimsArray(rand(10,10), (2,1));
    y5t = x5 .+ 1;
    @test LoopVectorization.check_args(x5)
    @test y5t == addone!(y5, x5)

    x6 = view(x5, 1:3, 1:3);
    y6 = view(y5, 1:3, 1:3);
    fill!(y6, NaN);
    @test LoopVectorization.check_args(x6)
    addone!(y6, x6); # also
    @test y5t == y5 # Test that `NaN`s replaced with correct answers, and that other values were untouched.

    A = rand(12,13,14,15);
    pA = PermutedDimsArray(A, (3,1,4,2));
    @test addone!(similar(pA), pA) == pA .+ 1
    ppA = PermutedDimsArray(pA, (4,2,3,1));
    @test LoopVectorization.check_args(ppA)
    @test addone!(similar(ppA), ppA) == ppA .+ 1

    x = rand(10,10); xc = copy(x);
    xv = view(x, 1:2:10, 1:3);
    xcv = view(xc, 1:2:10, 1:3);
    @test LoopVectorization.check_args(xv)
    addone!(xv, xv);
    xcv .+= 1;
    @test x == xc
end


text:
const julia_version = v"1.4.2"

include("../common.jl")

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
export fwi_objective

function fwi_objective(model_full::Model, source::judiVector, dObs::judiVector, srcnum::Integer; options=Options(), frequencies=[])
# Setup time-domain linear or nonlinear foward and adjoint modeling and interface to OPESCI/devito
    pm = load_pymodel()
    ac = load_acoustic_codegen()

    # Load full geometry for out-of-core geometry containers
    typeof(dObs.geometry) == GeometryOOC && (dObs.geometry = Geometry(dObs.geometry))
    typeof(source.geometry) == GeometryOOC && (source.geometry = Geometry(source.geometry))
    length(model_full.n) == 3 ? dims = [3,2,1] : dims = [2,1]   # model dimensions for Python are (z,y,x) and (z,x)

    # for 3D modeling, limit model to area with sources/receivers
    if options.limit_m == true && length(model_full.n) == 3 # only supported for 3D
        model = deepcopy(model_full)
        model = limit_model_to_receiver_area(source.geometry,dObs.geometry,model,options.buffer_size)
    else
        model = model_full
    end

    # Source/receiver parameters
    tmaxSrc = source.geometry.t[1]
    tmaxRec = dObs.geometry.t[1]

    # Set up Python model structure (force origin to be zero due to current devito bug)
    modelPy = pm."Model"(origin=(0.,0.,0.), spacing=model.d, shape=model.n, vp=process_physical_parameter(sqrt.(1f0./model.m), dims), nbpml=model.nb,
                       rho=process_physical_parameter(model.rho, dims), space_order=options.space_order)
    dtComp = modelPy.critical_dt

    # Extrapolate input data to computational grid
    qIn = time_resample(source.data[1],source.geometry,dtComp)[1]
    if typeof(dObs.data[1]) == SegyIO.SeisCon
        data = convert(Array{Float32,2},dObs.data[1][1].data)
        dObs = judiVector(dObs.geometry,data)
    end
    dObserved = time_resample(dObs.data[1],dObs.geometry,dtComp)[1]
    ntComp = size(dObserved,2)
    ntSrc = Int(trunc(tmaxSrc/dtComp+1))
    ntRec = Int(trunc(tmaxRec/dtComp+1))

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

    # Forward modeling to generate synthetic data and background wavefields
    if options.optimal_checkpointing == true
        op_F = pycall(ac."forward_modeling", PyObject, modelPy, PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(qIn))), PyReverseDims(copy(transpose(rec_coords))), op_return=true)
        argout1, argout2 = pycall(ac."adjoint_born", PyObject, modelPy, PyReverseDims(copy(transpose(rec_coords))), PyReverseDims(copy(transpose(dObserved))),
                                  op_forward=op_F, is_residual=false, n_checkpoints=options.num_checkpoints, maxmem=options.checkpoints_maxmem)
    elseif ~isempty(options.frequencies)
                typeof(options.frequencies) == Array{Any,1} && (options.frequencies = options.frequencies[srcnum])
                dPredicted, uf_real, uf_imag = pycall(ac."forward_freq_modeling", PyObject, modelPy, PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(qIn))), PyReverseDims(copy(transpose(rec_coords))),
                                                      options.frequencies, space_order=options.space_order, nb=model.nb, factor=options.dft_subsampling_factor)
                argout1 = .5f0*dot(vec(dPredicted) - vec(dObserved), vec(dPredicted) - vec(dObserved))*dtComp  # FWI objective function value
                argout2 = pycall(ac."adjoint_freq_born", Array{Float32, length(model.n)}, modelPy, PyReverseDims(copy(transpose(rec_coords))), PyReverseDims(copy(transpose((dPredicted - dObserved)))),
                                 options.frequencies, factor=options.dft_subsampling_factor, uf_real, uf_imag, space_order=options.space_order, nb=model.nb)
    else
        dPredicted, u0 = pycall(ac."forward_modeling", PyObject, modelPy, PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(qIn))), PyReverseDims(copy(transpose(rec_coords))), save=true, tsub_factor=options.subsampling_factor)
        argout1 = .5f0*dot(vec(dPredicted) - vec(dObserved), vec(dPredicted) - vec(dObserved))*dtComp # FWI objective function value
        argout2 = pycall(ac."adjoint_born", Array{Float32}, modelPy, PyReverseDims(copy(transpose(rec_coords))), PyReverseDims(copy(transpose((dPredicted  - dObserved)))), tsub_factor=options.subsampling_factor,
                         u=u0, is_residual=true)
    end
    argout2 = remove_padding(argout2, model.nb, true_adjoint=options.sum_padding)
    if options.limit_m==true && length(model_full.n) == 3
        argout2 = extend_gradient(model_full,model,argout2)
    end
    return [argout1; vec(argout2)]
end


text:
@doc raw"""
# Summary

    abstract type PaddableLayer <: Layer

Abstract Type to hold all Paddable Layers (i.e.  `ConvLayer` & `PoolLayer`)

# Subtypes

    ConvLayer
    PoolLayer

# Supertype Hierarchy

    PaddableLayer <: Layer <: Any
"""
abstract type PaddableLayer <: Layer end

@doc raw"""
# Summary

    abstract type ConvLayer <: PaddableLayer

Abstract Type to hold all ConvLayer

# Subtypes

    Conv1D
    Conv2D
    Conv3D

# Supertype Hierarchy

    ConvLayer <: PaddableLayer <: Layer <: Any
"""
abstract type ConvLayer <: PaddableLayer end

export ConvLayer

### Convolution layers

@doc raw"""
# Summary

    mutable struct Conv2D <: ConvLayer

# Fields

    channels    :: Integer
    f           :: Tuple{Integer,Integer}
    s           :: Tuple{Integer,Integer}
    inputS      :: Tuple
    outputS     :: Tuple
    padding     :: Symbol
    W           :: Array{F,4} where F
    dW          :: Array{F,4} where F
    K           :: Array{F,2} where F
    dK          :: Array{F,2} where F
    B           :: Array{F,4} where F
    dB          :: Array{F,4} where F
    actFun      :: Symbol
    keepProb    :: AbstractFloat
    V           :: Dict{Symbol,Array{F,4} where F}
    S           :: Dict{Symbol,Array{F,4} where F}
    V̂dk         :: Array{F,2} where F
    Ŝdk         :: Array{F,2} where F
    forwCount   :: Integer
    backCount   :: Integer
    updateCount :: Integer
    prevLayer   :: Union{Nothing, Layer}
    nextLayers  :: Array{Layer,1}

# Supertype Hierarchy

    Conv2D <: ConvLayer <: PaddableLayer <: Layer <: Any
"""
mutable struct Conv2D <: ConvLayer
    channels::Integer

    """
        filter size
    """
    f::Tuple{Integer, Integer}

    """
        stides size
    """
    s::Tuple{Integer, Integer}

    inputS::Tuple
    outputS::Tuple

    """
        Padding mode `:valid` or `:same`
    """
    padding::Symbol

    W::Array{F, 4} where {F}
    dW::Array{F,4} where {F}

    """
        Unrolled filters
    """
    K::Array{F, 2} where {F}
    dK::Array{F, 2} where {F}

    B::Array{F, 4} where {F}
    dB::Array{F, 4} where {F}

    actFun::Symbol
    keepProb::AbstractFloat

    # Z::Array{T,4} where {T}
    #
    # dA::Array{T,4} where {T}
    # A::Array{T,4} where {T}

    V::Dict{Symbol, Array{F, 4}  where {F}}
    S::Dict{Symbol, Array{F, 4}  where {F}}

    V̂dk::Array{F, 2} where {F}
    Ŝdk::Array{F, 2} where {F}

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    prevLayer::L where {L<:Union{Layer,Nothing}}
    nextLayers::Array{Layer,1}

    function Conv2D(
        c,
        f::Tuple{Integer, Integer};
        prevLayer=nothing,
        strides::Tuple{Integer, Integer}=(1, 1),
        padding::Symbol=:valid,
        activation::Symbol=:noAct,
        keepProb=1.0
    )

        if prevLayer == nothing
            T = Any
        elseif prevLayer isa Array
            T = eltype(prevLayer)
        else
            T = eltype(prevLayer.W)
        end

        new(c,
            f,
            strides,
            (0,), #inputS
            (0,), #outputS
            padding,
            Array{T,4}(undef,0,0,0,0), #W
            Array{T,4}(undef,0,0,0,0), #dW
            Matrix{T}(undef,0,0), #K
            Matrix{T}(undef,0,0), #dK
            Array{T,4}(undef,0,0,0,0), #B
            Array{T,4}(undef,0,0,0,0), #dB
            activation,
            keepProb,
            # Array{T,4}(undef, 0,0,0,0), #Z
            # Array{T,4}(undef, 0,0,0,0), #dA
            # Array{T,4}(undef, 0,0,0,0), #A
            Dict(:dw=>Array{T,4}(undef,0,0,0,0),
                 :db=>Array{T,4}(undef,0,0,0,0)), #V
            Dict(:dw=>Array{T,4}(undef,0,0,0,0),
                 :db=>Array{T,4}(undef,0,0,0,0)), #S
            Matrix{T}(undef,0,0), #V̂dk
            Matrix{T}(undef,0,0), #Ŝdk
            0, #forwCount
            0, #backCount
            0, #updateCount
            prevLayer,
            Array{Layer,1}(undef, 0)  #nextLayer
            )

    end #function Conv2D

end #mutable struct Conv2D

Base.show(io::IO, l::Conv2D) = print(io, "Conv2D($(l.channels), $(l.f), $(l.s), $(l.padding), $(l.actFun), $(l.inputS), $(l.outputS), ...)")

export Conv2D

@doc raw"""
# Summary

    mutable struct Conv1D <: ConvLayer

# Fields

    channels    :: Integer
    f           :: Integer
    s           :: Integer
    inputS      :: Tuple
    outputS     :: Tuple
    padding     :: Symbol
    W           :: Array{F,3} where F
    dW          :: Array{F,3} where F
    K           :: Array{F,2} where F
    dK          :: Array{F,2} where F
    B           :: Array{F,3} where F
    dB          :: Array{F,3} where F
    actFun      :: Symbol
    keepProb    :: AbstractFloat
    V           :: Dict{Symbol,Array{F,3} where F}
    S           :: Dict{Symbol,Array{F,3} where F}
    V̂dk         :: Array{F,2} where F
    Ŝdk         :: Array{F,2} where F
    forwCount   :: Integer
    backCount   :: Integer
    updateCount :: Integer
    prevLayer   :: Union{Nothing, Layer}
    nextLayers  :: Array{Layer,1}

# Supertype Hierarchy

    Conv1D <: ConvLayer <: PaddableLayer <: Layer <: Any
"""
mutable struct Conv1D <: ConvLayer
    channels::Integer

    """
        filter size
    """
    f::Integer

    """
        stides size
    """
    s::Integer

    inputS::Tuple
    outputS::Tuple

    padding::Symbol

    W::Array{F, 3}  where {F}
    dW::Array{F, 3}  where {F}

    """
        Unrolled filters
    """
    K::Array{F, 2} where {F}
    dK::Array{F, 2} where {F}

    B::Array{F, 3}  where {F}
    dB::Array{F, 3}  where {F}

    actFun::Symbol
    keepProb::AbstractFloat

    # Z::Array{T,3} where {T}
    # dA::Array{T,3} where {T}
    # A::Array{T,3} where {T}

    V::Dict{Symbol, Array{F, 3}  where {F}}
    S::Dict{Symbol, Array{F, 3}  where {F}}

    V̂dk::Array{F, 2} where {F}
    Ŝdk::Array{F, 2} where {F}

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    prevLayer::L where {L<:Union{Layer,Nothing}}
    nextLayers::Array{Layer,1}

    function Conv1D(c,
                    f::Integer;
                    prevLayer=nothing,
                    strides::Integer=1,
                    padding::Symbol=:valid,
                    activation::Symbol=:noAct,
                    keepProb=1.0)

        if prevLayer == nothing
            T = Any
        elseif prevLayer isa Array
            T = eltype(prevLayer)
        else
            T = eltype(prevLayer.W)
        end

        new(c,
            f,
            strides,
            (0,), #inputS
            (0,), #outputS
            padding,
            Array{T,3}(undef,0,0,0), #W
            Array{T,3}(undef,0,0,0), #dW
            Matrix{T}(undef,0,0), #K
            Matrix{T}(undef,0,0), #dK
            Array{T,3}(undef,0,0,0), #B
            Array{T,3}(undef,0,0,0), #dB
            activation,
            keepProb,
            # Array{T,3}(undef, 0,0,0), #Z
            # Array{T,3}(undef, 0,0,0), #dA
            # Array{T,3}(undef, 0,0,0), #A
            Dict(:dw=>Array{T,3}(undef,0,0,0),
                 :db=>Array{T,3}(undef,0,0,0)), #V
            Dict(:dw=>Array{T,3}(undef,0,0,0),
                 :db=>Array{T,3}(undef,0,0,0)), #S
            Matrix{T}(undef,0,0), #V̂dk
            Matrix{T}(undef,0,0), #Ŝdk
            0, #forwCount
            0, #backCount
            0, #updateCount
            prevLayer,
            Array{Layer,1}(undef, 0)  #nextLayer
            )

    end #function Conv1D

end #mutable struct Conv1D

Base.show(io::IO, l::Conv1D) = print(io, "Conv1D($(l.channels), $(l.f), $(l.s), $(l.padding), $(l.actFun), $(l.inputS), $(l.outputS), ...)")

export Conv1D

@doc raw"""
# Summary

    mutable struct Conv3D <: ConvLayer

# Fields

    channels    :: Integer
    f           :: Tuple{Integer,Integer,Integer}
    s           :: Tuple{Integer,Integer,Integer}
    inputS      :: Tuple
    outputS     :: Tuple
    padding     :: Symbol
    W           :: Array{F,5} where F
    dW          :: Array{F,5} where F
    K           :: Array{F,2} where F
    dK          :: Array{F,2} where F
    B           :: Array{F,5} where F
    dB          :: Array{F,5} where F
    actFun      :: Symbol
    keepProb    :: AbstractFloat
    V           :: Dict{Symbol,Array{F,5} where F}
    S           :: Dict{Symbol,Array{F,5} where F}
    V̂dk        :: Array{F,2} where F
    Ŝdk         :: Array{F,2} where F
    forwCount   :: Integer
    backCount   :: Integer
    updateCount :: Integer
    prevLayer   :: Union{Nothing, Layer}
    nextLayers  :: Array{Layer,1}

# Supertype Hierarchy

    Conv3D <: ConvLayer <: PaddableLayer <: Layer <: Any
"""
mutable struct Conv3D <: ConvLayer
    channels::Integer

    """
        filter size
    """
    f::Tuple{Integer, Integer, Integer}

    """
        stides size
    """
    s::Tuple{Integer, Integer, Integer}

    inputS::Tuple
    outputS::Tuple

    padding::Symbol

    W::Array{F, 5}  where {F}
    dW::Array{F, 5}  where {F}

    """
        Unrolled filters
    """
    K::Array{F, 2} where {F}
    dK::Array{F, 2} where {F}

    B::Array{F, 5}  where {F}
    dB::Array{F, 5}  where {F}

    actFun::Symbol
    keepProb::AbstractFloat

    # Z::Array{T,5} where {T}
    # dA::Array{T,5} where {T}
    # A::Array{T,5} where {T}

    V::Dict{Symbol, Array{F, 5}  where {F}}
    S::Dict{Symbol, Array{F, 5}  where {F}}

    V̂dk::Array{F, 2} where {F}
    Ŝdk::Array{F, 2} where {F}

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    prevLayer::L where {L<:Union{Layer,Nothing}}
    nextLayers::Array{Layer,1}

    function Conv3D(c,
                    f::Tuple{Integer, Integer, Integer};
                    prevLayer=nothing,
                    strides::Tuple{Integer, Integer, Integer}=(1, 1, 1),
                    padding::Symbol=:valid,
                    activation::Symbol=:noAct,
                    keepProb=1.0)

        if prevLayer == nothing
            T = Any
        elseif prevLayer isa Array
            T = eltype(prevLayer)
        else
            T = eltype(prevLayer.W)
        end

        new(c,
            f,
            strides,
            (0,), #inputS
            (0,), #outputS
            padding,
            Array{T,5}(undef,0,0,0,0,0), #W
            Array{T,5}(undef,0,0,0,0,0), #dW
            Matrix{T}(undef,0,0), #K
            Matrix{T}(undef,0,0), #dK
            Array{T,5}(undef,0,0,0,0,0), #B
            Array{T,5}(undef,0,0,0,0,0), #dB
            activation,
            keepProb,
            # Array{T,5}(undef, 0,0,0,0,0), #Z
            # Array{T,5}(undef, 0,0,0,0,0), #dA
            # Array{T,5}(undef, 0,0,0,0,0), #A
            Dict(:dw=>Array{T,5}(undef,0,0,0,0,0),
                 :db=>Array{T,5}(undef,0,0,0,0,0)), #V
            Dict(:dw=>Array{T,5}(undef,0,0,0,0,0),
                 :db=>Array{T,5}(undef,0,0,0,0,0)), #S
            Matrix{T}(undef,0,0), #V̂dk
            Matrix{T}(undef,0,0), #Ŝdk
            0, #forwCount
            0, #backCount
            0, #updateCount
            prevLayer,
            Array{Layer,1}(undef, 0)  #nextLayer
            )

    end #function Conv3D

end #mutable struct Conv3D

Base.show(io::IO, l::Conv3D) = print(io, "Conv3D($(l.channels), $(l.f), $(l.s), $(l.padding), $(l.actFun), $(l.inputS), $(l.outputS), ...)")

export Conv3D

###  Pooling layers

@doc raw"""
# Summary

    abstract type PoolLayer <: PaddableLayer

Abstract Type to hold all the `PoolLayer`s

# Subtypes

    AveragePoolLayer
    MaxPoolLayer

# Supertype Hierarchy

    PoolLayer <: PaddableLayer <: Layer <: Any
"""
abstract type PoolLayer <: PaddableLayer end

export PoolLayer

@doc raw"""
# Summary

    abstract type MaxPoolLayer <: PoolLayer

Abstract Type to hold all the `MaxPoolLayer`s

# Subtypes

    MaxPool1D
    MaxPool2D
    MaxPool3D

# Supertype Hierarchy

    MaxPoolLayer <: PoolLayer <: PaddableLayer <: Layer <: Any
"""
abstract type MaxPoolLayer <: PoolLayer end

export MaxPoolLayer

@doc """
    MaxPool2D(
        f::Tuple{Integer,Integer}=(2,2);
        prevLayer=nothing,
        strides::Tuple{Integer,Integer}=f,
        padding::Symbol=:valid,
    )

# Summary

    mutable struct MaxPool2D <: MaxPoolLayer

# Fields

    channels    :: Integer
    f           :: Tuple{Integer,Integer}
    s           :: Tuple{Integer,Integer}
    inputS      :: Tuple
    outputS     :: Tuple
    padding     :: Symbol
    forwCount   :: Integer
    backCount   :: Integer
    updateCount :: Integer
    prevLayer   :: Union{Nothing, Layer}
    nextLayers  :: Array{Layer,1}

# Supertype Hierarchy

    MaxPool2D <: MaxPoolLayer <: PoolLayer <: PaddableLayer <: Layer <: Any
"""
mutable struct MaxPool2D <: MaxPoolLayer
    channels::Integer

    """
        filter size
    """
    f::Tuple{Integer, Integer}

    """
        stides size
    """
    s::Tuple{Integer, Integer}

    inputS::Tuple
    outputS::Tuple

    padding::Symbol

    # dA::Array{T,4} where {T}
    # A::Array{T,4} where {T}

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    prevLayer::L where {L<:Union{Layer,Nothing}}
    nextLayers::Array{Layer,1}

    function MaxPool2D(
        f::Tuple{Integer,Integer}=(2,2);
        prevLayer=nothing,
        strides::Tuple{Integer,Integer}=f,
        padding::Symbol=:valid,
    )

        if prevLayer == nothing
           T = Any
        elseif prevLayer isa Array
           T = eltype(prevLayer)
        else
           T = eltype(prevLayer.W)
        end

        return new(
                   0, #channels
                   f,
                   strides,
                   (0,), #inputS
                   (0,), #outputS
                   padding,
                   # Array{T,4}(undef,0,0,0,0), #dA
                   # Array{T,4}(undef,0,0,0,0), #A
                   0, #forwCount
                   0, #backCount
                   0, #updateCount
                   prevLayer,
                   Array{Layer,1}(undef,0),
                   )
    end #function MaxPool2D

end #mutable struct MaxPool2D

Base.show(io::IO, l::MaxPool2D) = print(io, "MaxPool2D($(l.channels), $(l.f), $(l.s), $(l.padding), $(l.inputS), $(l.outputS), ...)")

export MaxPool2D

@doc """
    MaxPool1D(
        f::Integer=2;
        prevLayer=nothing,
        strides::Integer=f,
        padding::Symbol=:valid,
    )

# Summary

    mutable struct MaxPool1D <: MaxPoolLayer

# Fields

    channels    :: Integer
    f           :: Integer
    s           :: Integer
    inputS      :: Tuple
    outputS     :: Tuple
    padding     :: Symbol
    forwCount   :: Integer
    backCount   :: Integer
    updateCount :: Integer
    prevLayer   :: Union{Nothing, Layer}
    nextLayers  :: Array{Layer,1}

# Supertype Hierarchy

    MaxPool1D <: MaxPoolLayer <: PoolLayer <: PaddableLayer <: Layer <: Any
"""
mutable struct MaxPool1D <: MaxPoolLayer
    channels::Integer

    """
        filter size
    """
    f::Integer

    """
        stides size
    """
    s::Integer

    inputS::Tuple
    outputS::Tuple

    padding::Symbol

    # dA::Array{T,3} where {T}
    # A::Array{T,3} where {T}

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    prevLayer::L where {L<:Union{Layer,Nothing}}
    nextLayers::Array{Layer,1}

    function MaxPool1D(
        f::Integer=2;
        prevLayer=nothing,
        strides::Integer=f,
        padding::Symbol=:valid,
    )

        if prevLayer == nothing
           T = Any
        elseif prevLayer isa Array
           T = eltype(prevLayer)
        else
           T = eltype(prevLayer.W)
        end

        return new(
                   0, #channels
                   f,
                   strides,
                   (0,), #inputS
                   (0,), #outputS
                   padding,
                   # Array{T,3}(undef,0,0,0), #dA
                   # Array{T,3}(undef,0,0,0), #A
                   0, #forwCount
                   0, #backCount
                   0, #updateCount
                   prevLayer,
                   Array{Layer,1}(undef,0),
                   )
    end #function MaxPool1D

end #mutable struct MaxPool1D

Base.show(io::IO, l::MaxPool1D) = print(io, "MaxPool1D($(l.channels), $(l.f), $(l.s), $(l.padding), $(l.inputS), $(l.outputS), ...)")

export MaxPool1D

@doc """
    MaxPool3D(
        f::Tuple{Integer,Integer,Integer}=(2,2,2);
        prevLayer=nothing,
        strides::Tuple{Integer,Integer,Integer}=f,
        padding::Symbol=:valid,
    )

# Summary

    mutable struct MaxPool3D <: MaxPoolLayer

# Fields

    channels    :: Integer
    f           :: Tuple{Integer,Integer,Integer}
    s           :: Tuple{Integer,Integer,Integer}
    inputS      :: Tuple
    outputS     :: Tuple
    padding     :: Symbol
    forwCount   :: Integer
    backCount   :: Integer
    updateCount :: Integer
    prevLayer   :: Union{Nothing, Layer}
    nextLayers  :: Array{Layer,1}

# Supertype Hierarchy

    MaxPool3D <: MaxPoolLayer <: PoolLayer <: PaddableLayer <: Layer <: Any
"""
mutable struct MaxPool3D <: MaxPoolLayer
    channels::Integer

    """
        filter size
    """
    f::Tuple{Integer, Integer, Integer}

    """
        stides size
    """
    s::Tuple{Integer, Integer, Integer}

    inputS::Tuple
    outputS::Tuple

    padding::Symbol

    # dA::Array{T,5} where {T}
    # A::Array{T,5} where {T}

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    prevLayer::L where {L<:Union{Layer,Nothing}}
    nextLayers::Array{Layer,1}

    function MaxPool3D(
        f::Tuple{Integer,Integer,Integer}=(2,2,2);
        prevLayer=nothing,
        strides::Tuple{Integer,Integer,Integer}=f,
        padding::Symbol=:valid,
    )

        if prevLayer == nothing
           T = Any
        elseif prevLayer isa Array
           T = eltype(prevLayer)
        else
           T = eltype(prevLayer.W)
        end

        return new(
                   0, #channels
                   f,
                   strides,
                   (0,), #inputS
                   (0,), #outputS
                   padding,
                   # Array{T,5}(undef,0,0,0,0,0), #dA
                   # Array{T,5}(undef,0,0,0,0,0), #A
                   0, #forwCount
                   0, #backCount
                   0, #updateCount
                   prevLayer,
                   Array{Layer,1}(undef,0),
                   )
    end #function MaxPool3D

end #mutable struct MaxPool3D

Base.show(io::IO, l::MaxPool3D) = print(io, "MaxPool3D($(l.channels), $(l.f), $(l.s), $(l.padding), $(l.inputS), $(l.outputS), ...)")

export MaxPool3D

#### AveragePoolLayer

@doc raw"""
# Summary

    abstract type AveragePoolLayer <: PoolLayer

# Subtypes

    AveragePool1D
    AveragePool2D
    AveragePool3D

# Supertype Hierarchy

    AveragePoolLayer <: PoolLayer <: PaddableLayer <: Layer <: Any
"""
abstract type AveragePoolLayer <: PoolLayer end

export AveragePoolLayer

@doc raw"""
    AveragePool2D(
        f::Tuple{Integer,Integer}=(2,2);
        prevLayer=nothing,
        strides::Tuple{Integer,Integer}=f,
        padding::Symbol=:valid,
    )

# Summary

    mutable struct AveragePool2D <: AveragePoolLayer

# Fields

    channels    :: Integer
    f           :: Tuple{Integer,Integer}
    s           :: Tuple{Integer,Integer}
    inputS      :: Tuple
    outputS     :: Tuple
    padding     :: Symbol
    forwCount   :: Integer
    backCount   :: Integer
    updateCount :: Integer
    prevLayer   :: Union{Nothing, Layer}
    nextLayers  :: Array{Layer,1}

# Supertype Hierarchy

    AveragePool2D <: AveragePoolLayer <: PoolLayer <: PaddableLayer <: Layer <: Any
"""
mutable struct AveragePool2D <: AveragePoolLayer
    channels::Integer

    """
        filter size
    """
    f::Tuple{Integer, Integer}

    """
        stides size
    """
    s::Tuple{Integer, Integer}

    inputS::Tuple
    outputS::Tuple

    padding::Symbol

    # dA::Array{T,4} where {T}
    # A::Array{T,4} where {T}

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    prevLayer::L where {L<:Union{Layer,Nothing}}
    nextLayers::Array{Layer,1}

    function AveragePool2D(
        f::Tuple{Integer,Integer}=(2,2);
        prevLayer=nothing,
        strides::Tuple{Integer,Integer}=f,
        padding::Symbol=:valid,
    )

        if prevLayer == nothing
           T = Any
        elseif prevLayer isa Array
           T = eltype(prevLayer)
        else
           T = eltype(prevLayer.W)
        end

        return new(
                   0, #channels
                   f,
                   strides,
                   (0,), #inputS
                   (0,), #outputS
                   padding,
                   # Array{T,4}(undef,0,0,0,0), #dA
                   # Array{T,4}(undef,0,0,0,0), #A
                   0, #forwCount
                   0, #backCount
                   0, #updateCount
                   prevLayer,
                   Array{Layer,1}(undef,0),
                   )
    end #function AveragePool2D

end #mutable struct AveragePool2D

Base.show(io::IO, l::AveragePool2D) = print(io, "AveragePool2D($(l.channels), $(l.f), $(l.s), $(l.padding), $(l.inputS), $(l.outputS), ...)")

export AveragePool2D

@doc raw"""
    AveragePool1D(
        f::Integer=2;
        prevLayer=nothing,
        strides::Integer=f,
        padding::Symbol=:valid,
    )

# Summary

    mutable struct AveragePool1D <: AveragePoolLayer

# Fields

    channels    :: Integer
    f           :: Integer
    s           :: Integer
    inputS      :: Tuple
    outputS     :: Tuple
    padding     :: Symbol
    forwCount   :: Integer
    backCount   :: Integer
    updateCount :: Integer
    prevLayer   :: Union{Nothing, Layer}
    nextLayers  :: Array{Layer,1}

# Supertype Hierarchy

    AveragePool1D <: AveragePoolLayer <: PoolLayer <: PaddableLayer <: Layer <: Any
"""
mutable struct AveragePool1D <: AveragePoolLayer
    channels::Integer

    """
        filter size
    """
    f::Integer

    """
        stides size
    """
    s::Integer

    inputS::Tuple
    outputS::Tuple

    padding::Symbol

    # dA::Array{T,3} where {T}
    # A::Array{T,3} where {T}

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    prevLayer::L where {L<:Union{Layer,Nothing}}
    nextLayers::Array{Layer,1}

    function AveragePool1D(
        f::Integer=2;
        prevLayer=nothing,
        strides::Integer=f,
        padding::Symbol=:valid,
    )

        if prevLayer == nothing
           T = Any
        elseif prevLayer isa Array
           T = eltype(prevLayer)
        else
           T = eltype(prevLayer.W)
        end

        return new(
                   0, #channels
                   f,
                   strides,
                   (0,), #inputS
                   (0,), #outputS
                   padding,
                   # Array{T,3}(undef,0,0,0), #dA
                   # Array{T,3}(undef,0,0,0), #A
                   0, #forwCount
                   0, #backCount
                   0, #updateCount
                   prevLayer,
                   Array{Layer,1}(undef,0),
                   )
    end #function AveragePool1D

end #mutable struct AveragePool1D

Base.show(io::IO, l::AveragePool1D) = print(io, "AveragePool1D($(l.channels), $(l.f), $(l.s), $(l.padding), $(l.inputS), $(l.outputS), ...)")

export AveragePool1D

@doc raw"""
    AveragePool3D(
        f::Tuple{Integer,Integer,Integer}=(2,2,2);
        prevLayer=nothing,
        strides::Tuple{Integer,Integer,Integer}=f,
        padding::Symbol=:valid,
    )

# Summary

    mutable struct AveragePool3D <: AveragePoolLayer

# Fields

    channels    :: Integer
    f           :: Tuple{Integer,Integer,Integer}
    s           :: Tuple{Integer,Integer,Integer}
    inputS      :: Tuple
    outputS     :: Tuple
    padding     :: Symbol
    forwCount   :: Integer
    backCount   :: Integer
    updateCount :: Integer
    prevLayer   :: Union{Nothing, Layer}
    nextLayers  :: Array{Layer,1}

# Supertype Hierarchy

    AveragePool3D <: AveragePoolLayer <: PoolLayer <: PaddableLayer <: Layer <: Any
"""
mutable struct AveragePool3D <: AveragePoolLayer
    channels::Integer

    """
        filter size
    """
    f::Tuple{Integer, Integer, Integer}

    """
        stides size
    """
    s::Tuple{Integer, Integer, Integer}

    inputS::Tuple
    outputS::Tuple

    padding::Symbol

    # dA::Array{T,5} where {T}
    # A::Array{T,5} where {T}

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    prevLayer::L where {L<:Union{Layer,Nothing}}
    nextLayers::Array{Layer,1}

    function AveragePool3D(
        f::Tuple{Integer,Integer,Integer}=(2,2,2);
        prevLayer=nothing,
        strides::Tuple{Integer,Integer,Integer}=f,
        padding::Symbol=:valid,
    )

        if prevLayer == nothing
           T = Any
        elseif prevLayer isa Array
           T = eltype(prevLayer)
        else
           T = eltype(prevLayer.W)
        end

        return new(
                   0, #channels
                   f,
                   strides,
                   (0,), #inputS
                   (0,), #outputS
                   padding,
                   # Array{T,5}(undef,0,0,0,0,0), #dA
                   # Array{T,5}(undef,0,0,0,0,0), #A
                   0, #forwCount
                   0, #backCount
                   0, #updateCount
                   prevLayer,
                   Array{Layer,1}(undef,0),
                   )
    end #function AveragePool3D

end #mutable struct AveragePool3D

Base.show(io::IO, l::AveragePool3D) = print(io, "AveragePool3D($(l.channels), $(l.f), $(l.s), $(l.padding), $(l.inputS), $(l.outputS), ...)")

export AveragePool3D


text:
mutable struct Travis <: Remote
    repo_name::String
    username::String
    token::String
    user_code::Int
    repo_code::Int
end

Travis(repo_name) =
    Travis(repo_name, settings("username"), settings("travis_token"), 0, 0)

base_url(t::Travis) = "https://api.travis-ci.org"

headers(t::Travis) = Dict(
    "Travis-API-Version" => "3",
    "Authorization" => "token $(t.token)",
    "User-Agent" => user_agent)

function user(t::Travis)
    @info "getting travis user"
    talk_to(HTTP.get, t, "/user") |> json_parse
end

function sync(t::Travis)
    @info "syncing travis"
    talk_to(HTTP.post, t, "/user/$(t.user_code)/sync")
end

function repos(t::Travis)
    @info "getting travis repos"
    talk_to(HTTP.get, t, "/repos") |> json_parse
end

function repo(t::Travis)
    @info "getting travis repo"
    talk_to(HTTP.get, t, "/repo/$(t.username)%2F$(t.repo_name)") |> json_parse
end

function activate(t::Travis)
    @info "creating travis"
    talk_to(HTTP.post, t, "/repo/$(t.repo_code)/activate")
end

function get_keys(t::Travis)
    @info "getting travis keys"
    talk_to(HTTP.get, t, "/repo/$(t.repo_code)/env_vars") |> json_parse
end

function delete_key(t::Travis, key_id)
    @info "deleting travis key"
    talk_to(HTTP.delete, t, "/repo/$(t.repo_code)/env_var/$key_id")
end

function delete_keys(t::Travis, name)
    foreach(
        travis_key ->
            if travis_key["name"] == name
                delete_key(t, travis_key["id"])
            end,
        get_keys(t)["env_vars"]
    )
end

function add_key(t::Travis, name, value; public = false)
    @info "creating travis key"
    talk_to(HTTP.post, t, "/repo/$(t.repo_code)/env_vars", Dict(
        "env_var.name" => name,
        "env_var.value" => value,
        "env_var.public" => public
    ))
end

function exists(t::Travis)
    any(repos(t)["repositories"]) do repo
        repo["slug"] == user_repo(t)
    end
end

function user!(t::Travis)
    t.user_code = user(t)["id"]
end

function repo!(t::Travis)
    t.repo_code = repo(t)["id"]
end


text:
ENV["TMPDIR"] = pwd()
(fname, tio) = mktemp()
println(fname, " created as a \"temporary\" file.")
println(tio, msg)
close(tio)
println("\"", msg, "\" written to ", fname)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#
# The given combinator.
#

function GivenQuery(base::Query, param::Query)
    tag = decoration(output(param), :tag, Symbol, Symbol(""))
    @assert tag != Symbol("")
    remparams = filter(p -> p.first != tag, slots(input(base)))
    q = RecordQuery(istub(base), param)
    q >>
    Query(
        GivenSig(),
        [FieldQuery(output(q), 1) >> base],
        Input(domain(output(q)), mode(input(base)) |> setslots(remparams)),
        output(base))
end

function GivenQuery(base::Query)
    oty = fields(domain(input(base)))[2]
    tag = decoration(oty, :tag, Symbol, Symbol(""))
    @assert tag != Symbol("")
    remparams = filter(p -> p.first != tag, slots(input(base)))
    Query(
        GivenSig(),
        [base],
        Input(domain(input(base)), mode(input(base)) |> setslots(remparams)),
        output(base))
end

GivenQuery(base::Query, param::Pair{Symbol, Query}) =
    GivenQuery(base, param.second |> decorate(:tag => param.first))

GivenQuery(base::Query, params::Vector{Query}) =
    isempty(params) ? base : GivenQuery(GivenQuery(base, params[1:end-1]), params[end])

GivenQuery(base::Query, param, params...) =
    GivenQuery(GivenQuery(base, params...), param)

immutable GivenSig <: AbstractSignature
end

ev(::GivenSig, args::Vector{Query}, ::Input, oty::Output, iflow::InputFlow) =
    OutputFlow(
        oty,
        given_impl(args..., iflow.ctx, iflow.frameoffs, iflow.slotflows, values(iflow), domain(iflow)))

function given_impl(
        arg::Query, ctx::InputContext, frameoffs::InputFrame, slotflows::InputSlotFlows, dv::DataVector,
        dom::Domain)
    fs = fields(dom)
    tag = decoration(fs[2], :tag, Symbol, Symbol(""))
    pmap = Dict{Symbol,OutputFlow}(slotflows)
    pmap[tag] = OutputFlow(fs[2], column(dv, 2))
    pkeys = collect(keys(pmap))
    sort!(pkeys)
    pflows = InputSlotFlow[pkey => pmap[pkey] for pkey in pkeys]
    iflow′ = InputFlow(
        ctx,
        dom,
        dv,
        frameoffs,
        pflows)
    return column(ev(arg, iflow′))
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
module BlockBootstrap

using Random

include("samplingFunctions.jl")
include("bootstrapSampling.jl")
export MBBsampler, NBBsampler, CBBsampler

include("bootstrapStatistic.jl")
export bootstrapStatistic

end


text:
function format(x)
    io = IOBuffer()
    show(IOContext(io, :compact => true, :typeinfo => typeof(x)), x)
    String(take!(io))
end

function row(r, i; format = TableWidgets.format)
	fields = propertynames(r)

    node("tr",
        node("th", format(i)),
        (node("td", format(getproperty(r, field))) for field in fields)...)
end

rendertable(t; kwargs...) = render_row_iterator(Tables.rows(t); kwargs...)
rendertable(t, r::Integer; kwargs...) = render_row_iterator(Iterators.take(Tables.rows(t), r); kwargs...)

function render_row_iterator(t;
    format = TableWidgets.format, className = "is-striped is-hoverable", row = TableWidgets.row)

    fr, lr = Iterators.peel(t)

    names = propertynames(fr)
    headers = node("tr", node("th", ""), (node("th", string(n)) for n in names)...) |> node("thead")

    first_row = row(fr, 1; format = format)
    body = node("tbody", first_row, (row(r, i+1; format = format) for (i, r) in enumerate(lr))...)
    className = "table interact-widget $className"
    n = slap_design!(node("table", headers, body, className = className))
end

"""
`head(t, r=6)`

Show first `r` rows of table `t` as HTML table.
"""
function head(t, r = 6; kwargs...)
    t isa AbstractObservable || (t = Observable{Any}(t))
    r isa AbstractObservable || (r = Observable{Int}(r))
    h = @map rendertable(&t, &r; kwargs...)

    Widget{:head}([:rows => r, :head => h], output = t, layout = i -> i[:head])
end

function toggled(wdg::AbstractWidget; readout = true, label = "Show")
    toggled_wdg = togglecontent(wdg, label = label, value = readout)
    Widget{:toggled}([:toggle => toggled_wdg], output = observe(wdg), layout = i -> i[:toggle])
end

"""
`dataeditor(t, rows; label = "Show table")`

Create a textbox to preprocess a table: displays the result using `toggled(head(t, rows))`.
"""
function dataeditor(t, args...; readout = true, label = "Show table", kwargs...)
    (t isa AbstractObservable) || (t = Observable{Any}(t))
    output = Observable{Any}(t[])
    wdg = Widget{:dataeditor}(output = output)
    wdg[:input] = t
    wdg[:display] = toggled(head(observe(wdg), args...; kwargs...); readout = readout, label = label)
    wdg[:text] = textarea(placeholder = "Write transformation to apply to the table")
    parsetext!(wdg; text = observe(wdg, :text), on = (observe(wdg, :text)), default = identity)
    wdg[:apply] = button("Apply")
    wdg[:reset] = button("Reset", className = "is-danger")
    @map! observe(wdg) begin
        &wdg[:apply]
        (wdg[:function][])(&t)
    end
    @on begin
        &wdg[:reset]
        wdg[:text][] = ""
        observe(wdg)[] = t[]
    end
    @layout! wdg Widgets.div(:text, hbox(:apply, hskip(1em), :reset), :display)
end


text:
# For `AbstractArray`, we do not differentiate between left and right indices:

"""
    checked_similar_from_indices(C, T, indleft, indright, A, conjA = :N)

Returns an object similar to `A` which has an `eltype` given by `T` and whose left indices
correspond to the indices `indleft` from `op(A)`, and its right indices correspond to the
indices `indright` from `op(A)`, where `op` is `conj` if `conjA == :C` or does nothing if
`conjA == :N` (default). Here, `C` is a potential candidate for the similar object. If
`C === nothing`, or its its `eltype` or shape does not match, a new object is allocated and
returned. Otherwise, `C` is returned.
"""
checked_similar_from_indices(C, T::Type, p1::IndexTuple, p2::IndexTuple, A::AbstractArray,
    CA::Symbol = :N) = checked_similar_from_indices(C, T, (p1..., p2...), A, CA)

"""
    checked_similar_from_indices(C, T, indoA, indoB, indleft, indright,
                                    A, B, conjA = :N, conjB= :N)

Returns an object similar to `A` which has an `eltype` given by `T` and dimensions/sizes
corresponding to a selection of those of `opA(A)` and `opB(B)` concatenated. Out of the
collection of indices in `indoA` of `opA(A)` and `indoB` of `opB(B)`, we construct an
object whose left (right) indices correspond to indices `indleft` (`indright`) from that
collection. Furthermore, `C` is a potential candidate for the similar object. If
`C === nothing`, or its its `eltype` or shape does not match, a new object is allocated and
returned. Otherwise, `C` is returned.
"""
checked_similar_from_indices(C, T::Type, poA::IndexTuple, poB::IndexTuple, p1::IndexTuple,
    p2::IndexTuple, A::AbstractArray, B::AbstractArray, CA::Symbol = :N, CB::Symbol = :N) =
    checked_similar_from_indices(C, T, poA, poB, (p1..., p2...), A, B, CA, CB)

"""
    scalar(C)

Returns the single element of a tensor-like object with zero indices or dimensions.
"""
function scalar end

"""
    add!(α, A, conjA, β, C, indleft, indright)

Implements `C = β*C+α*permute(op(A))` where `A` is permuted such that the left (right)
indices of `C` correspond to the indices `indleft` (`indright`) of `A`, and `op` is `conj`
if `conjA == :C` or the identity map if `conjA == :N` (default). Together,
`(indleft..., indright...)` is a permutation of 1 to the number of indices (dimensions) of
`A`.
"""
add!(α, A::AbstractArray, CA::Symbol, β, C::AbstractArray, indleft::IndexTuple,
        indright::IndexTuple) = add!(α, A, CA, β, C, (indleft..., indright...))

"""
    trace!(α, A, conjA, β, C, indleft, indright, cind1, cind2)

Implements `C = β*C+α*partialtrace(op(A))` where `A` is permuted and partially traced,
such that the left (right) indices of `C` correspond to the indices `indleft` (`indright`)
of `A`, and indices `cindA1` are contracted with indices `cindA2`. Furthermore, `op` is
`conj` if `conjA == :C` or the identity map if `conjA=:N` (default). Together,
`(indleft..., indright..., cind1, cind2)` is a permutation of 1 to the number of indices
(dimensions) of `A`.
"""
trace!(α, A::AbstractArray, CA::Symbol, β, C::AbstractArray, indleft::IndexTuple,
        indright::IndexTuple, cind1::IndexTuple, cind2::IndexTuple) =
    trace!(α, A, CA, β, C, (indleft..., indright...), cind1, cind2)

"""
    contract!(α, A, conjA, B, conjB, β, C, oindA, cindA, oindB, cindB, indleft, indright, syms = nothing)

Implements `C = β*C+α*contract(opA(A),opB(B))` where `A` and `B` are contracted, such that
the indices `cindA` of `A` are contracted with indices `cindB` of `B`. The open indices
`oindA` of `A` and `oindB` of `B` are permuted such that `C` has left (right) indices
corresponding to indices `indleft` (`indright`) out of `(oindA..., oindB...)`. The
operation `opA` (`opB`) acts as `conj` if `conjA` (`conjB`) equal `:C` or as the identity
map if `conjA` (`conjB`) equal `:N`. Together, `(oindA..., cindA...)` is a permutation of
1 to the number of indices of `A` and `(oindB..., cindB...)` is a permutation of 1 to the
number of indices of `C`. Furthermore, `length(cindA) == length(cindB)`,
`length(oindA)+length(oindB)` equals the number of indices of `C` and `(indleft...,
indright...)` is a permutation of `1` ot the number of indices of `C`.

The final argument `syms` is optional and can be either `nothing`, or a tuple of three
symbols, which are used to identify temporary objects in the cache to be used for permuting
`A`, `B` and `C` so as to perform the contraction as a matrix multiplication.
"""
contract!(α, A::AbstractArray, CA::Symbol, B::AbstractArray, CB::Symbol,
        β, C::AbstractArray, oindA::IndexTuple, cindA::IndexTuple, oindB::IndexTuple,
        cindB::IndexTuple, indleft::IndexTuple, indright::IndexTuple, syms = nothing) =
    contract!(α, A, CA, B, CB, β, C,
                oindA, cindA, oindB, cindB, (indleft..., indright...), syms)

# actual implementations for AbstractArray with ind = (indleft..., indright...)
Base.@pure function similartype(A, T, sz)
    Core.Compiler.return_type(similar, Tuple{typeof(A), Type{T}, typeof(sz)})
end

function checked_similar_from_indices(C, ::Type{T}, ind::IndexTuple{N}, A::AbstractArray,
        CA::Symbol) where {T,N}

    sz = map(n->size(A, n), ind)
    if C !== nothing && C isa AbstractArray && sz == size(C) && T == eltype(C)
        CT = similartype(A, T, sz)
        return C::CT
    else
        return similar(A, T, sz)
    end
end
function checked_similar_from_indices(C, ::Type{T}, poA::IndexTuple, poB::IndexTuple,
        ind::IndexTuple{N}, A::AbstractArray, B::AbstractArray,
        CA::Symbol, CB::Symbol) where {T,N}

    oszA = map(n->size(A,n), poA)
    oszB = map(n->size(B,n), poB)
    sz = let osz = (oszA..., oszB...)
        map(n->osz[n], ind)
    end
    if C !== nothing && C isa AbstractArray && sz == size(C) && T == eltype(C)
        CT = similartype(A, T, sz)
        return C::CT
    else
        return similar(A, T, sz)
    end
end

scalar(C::AbstractArray) = ndims(C)==0 ? C[1] : throw(DimensionMismatch())

function add!(α, A::AbstractArray{<:Any, N}, CA::Symbol,
        β, C::AbstractArray{<:Any, N}, indCinA) where {N}

    N == length(indCinA) || throw(IndexError("Invalid permutation of length $N: $indCinA"))
    if CA == :N
        if isbitstype(eltype(A)) && isbitstype(eltype(C))
            @unsafe_strided A C _add!(α, A, β, C, (indCinA...,))
        else
            _add!(α, StridedView(A), β, StridedView(C), (indCinA...,))
        end
    elseif CA == :C
        if isbitstype(eltype(A)) && isbitstype(eltype(C))
            @unsafe_strided A C _add!(α, conj(A), β, C, (indCinA...,))
        else
            _add!(α, conj(StridedView(A)), β, StridedView(C), (indCinA...,))
        end
    elseif CA == :A
        if isbitstype(eltype(A)) && isbitstype(eltype(C))
            @unsafe_strided A C _add!(α, map(adjoint, A), β, C, (indCinA...,))
        else
            _add!(α, map(adjoint, StridedView(A)), β, StridedView(C), (indCinA...,))
        end
    else
        throw(ArgumentError("Unknown conjugation flag: $CA"))
    end
    return C
end
_add!(α, A::AbstractStridedView{<:Any,N},
        β, C::AbstractStridedView{<:Any,N}, indCinA::IndexTuple{N}) where N =
    LinearAlgebra.axpby!(α, permutedims(A, indCinA), β, C)

function trace!(α, A::AbstractArray{<:Any, NA}, CA::Symbol, β, C::AbstractArray{<:Any, NC},
        indCinA, cindA1, cindA2) where {NA,NC}

    NC == length(indCinA) ||
        throw(IndexError("Invalid selection of $NC out of $NA: $indCinA"))
    NA-NC == 2*length(cindA1) == 2*length(cindA2) ||
        throw(IndexError("invalid number of trace dimension"))
    if CA == :N
        if isbitstype(eltype(A)) && isbitstype(eltype(C))
            @unsafe_strided A C _trace!(α, A, β, C,
                (indCinA...,), (cindA1...,), (cindA2...,))
        else
            _trace!(α, StridedView(A), β, StridedView(C),
                (indCinA...,), (cindA1...,), (cindA2...,))
        end
    elseif CA == :C
        if isbitstype(eltype(A)) && isbitstype(eltype(C))
            @unsafe_strided A C _trace!(α, conj(A), β, C,
                (indCinA...,), (cindA1...,), (cindA2...,))
        else
            _trace!(α, conj(StridedView(A)), β, StridedView(C),
                (indCinA...,), (cindA1...,), (cindA2...,))
        end
    elseif CA == :A
        if isbitstype(eltype(A)) && isbitstype(eltype(C))
            @unsafe_strided A C _trace!(α, map(adjoint, A), β, C,
                (indCinA...,), (cindA1...,), (cindA2...,))
        else
            _trace!(α, map(adjoint, StridedView(A)), β, StridedView(C),
                (indCinA...,), (cindA1...,), (cindA2...,))
        end
    else
        throw(ArgumentError("Unknown conjugation flag: $CA"))
    end
    return C
end

function _trace!(α, A::AbstractStridedView,
        β, C::AbstractStridedView, indCinA::IndexTuple{NC},
        cindA1::IndexTuple{NT}, cindA2::IndexTuple{NT}) where {NC,NT}

    sizeA = i->size(A, i)
    strideA = i->stride(A, i)
    tracesize = sizeA.(cindA1)
    tracesize == sizeA.(cindA2) || throw(DimensionMismatch("non-matching trace sizes"))
    size(C) == sizeA.(indCinA) || throw(DimensionMismatch("non-matching sizes"))

    newstrides = (strideA.(indCinA)..., (strideA.(cindA1) .+ strideA.(cindA2))...)
    newsize = (size(C)..., tracesize...)
    if A isa UnsafeStridedView
        A2 = UnsafeStridedView(A.ptr, newsize, newstrides, A.offset, A.op)
    else
        A2 = StridedView(A.parent, newsize, newstrides, A.offset, A.op)
    end

    if α != 1
        if β == 0
            Strided._mapreducedim!(x->α*x, +, zero, newsize, (C, A2))
        elseif β == 1
            Strided._mapreducedim!(x->α*x, +, nothing, newsize, (C, A2))
        else
            Strided._mapreducedim!(x->α*x, +, y->β*y, newsize, (C, A2))
        end
    else
        if β == 0
            return Strided._mapreducedim!(identity, +, zero, newsize, (C, A2))
        elseif β == 1
            Strided._mapreducedim!(identity, +, nothing, newsize, (C, A2))
        else
            Strided._mapreducedim!(identity, +, y->β*y, newsize, (C, A2))
        end
    end
    return C
end

function contract!(α, A::AbstractArray, CA::Symbol, B::AbstractArray, CB::Symbol,
        β, C::AbstractArray,
        oindA::IndexTuple, cindA::IndexTuple, oindB::IndexTuple, cindB::IndexTuple,
        indCinoAB::IndexTuple, syms::Union{Nothing, NTuple{3,Symbol}} = nothing)

    trivCind = _trivtuple(indCinoAB)
    if  indCinoAB !== trivCind
        # try reverse order
        indCinoBA = getindex.(Ref(indCinoAB),
                        ((length(oindA) .+ _trivtuple(oindB))..., _trivtuple(oindA)...))
        if indCinoBA == trivCind
            return contract!(α, B, CB, A, CA, β, C,
                                oindB, cindB, oindA, cindA, indCinoBA, syms)
        end
    end

    pA = (oindA...,cindA...)
    (length(pA) == ndims(A) && TupleTools.isperm(pA)) ||
        throw(IndexError("invalid permutation of length $(ndims(A)): $pA"))
    pB = (oindB...,cindB...)
    (length(pB) == ndims(B) && TupleTools.isperm(pB)) ||
        throw(IndexError("invalid permutation of length $(ndims(B)): $pB"))
    (length(oindA) + length(oindB) == ndims(C)) ||
        throw(IndexError("non-matching output indices in contraction"))
    (ndims(C) == length(indCinoAB) && isperm(indCinoAB)) ||
        throw(IndexError("invalid permutation of length $(ndims(C)): $indCinoAB"))

    sizeA = i->size(A, i)
    sizeB = i->size(B, i)
    sizeC = i->size(C, i)

    csizeA = sizeA.(cindA)
    csizeB = sizeB.(cindB)
    osizeA = sizeA.(oindA)
    osizeB = sizeB.(oindB)

    csizeA == csizeB ||
        throw(DimensionMismatch("non-matching sizes in contracted dimensions"))
    sizeAB = let osize = (osizeA..., osizeB...)
        i->osize[i]
    end
    sizeAB.(indCinoAB) == size(C) ||
        throw(DimensionMismatch("non-matching sizes in uncontracted dimensions"))

    TC = eltype(C)
    if use_blas() && TC <: BlasFloat
        if isblascontractable(A, oindA, cindA, CA) && eltype(A) == TC
            A2 = A
            CA2 = CA
        else
            if syms === nothing
                A2 = similar_from_indices(TC, oindA, cindA, A, CA)
            else
                A2 = cached_similar_from_indices(syms[1], TC, oindA, cindA, A, CA)
            end
            add!(1, A, CA, 0, A2, oindA, cindA)
            CA2 = :N
            oindA = _trivtuple(oindA)
            cindA = _trivtuple(cindA) .+ length(oindA)
        end
        if isblascontractable(B, cindB, oindB, CB) && eltype(B) == TC
            B2 = B
            CB2 = CB
        else
            if syms === nothing
                B2 = similar_from_indices(TC, cindB, oindB, B, CB)
            else
                B2 = cached_similar_from_indices(syms[2], TC, cindB, oindB, B, CB)
            end
            add!(1, B, CB, 0, B2, cindB, oindB)
            CB2 = :N
            cindB = _trivtuple(cindB)
            oindB = _trivtuple(oindB) .+ length(cindB)
        end
        ipC = TupleTools.invperm(indCinoAB)
        oindAinC = TupleTools.getindices(ipC, _trivtuple(oindA))
        oindBinC = TupleTools.getindices(ipC, length(oindA) .+ _trivtuple(oindB))
        if isblascontractable(C, oindAinC, oindBinC, :D)
            C2 = C
            _blas_contract!(α, A2, CA2, B2, CB2, β, C2,
                                oindA, cindA, oindB, cindB, oindAinC, oindBinC)
        else
            if syms === nothing
                C2 = similar_from_indices(TC, oindAinC, oindBinC, C, :N)
            else
                C2 = cached_similar_from_indices(syms[3], TC, oindAinC, oindBinC, C, :N)
            end
            _blas_contract!(1, A2, CA2, B2, CB2, 0, C2,
                                oindA, cindA, oindB, cindB,
                                _trivtuple(oindA), length(oindA) .+ _trivtuple(oindB))
            add!(α, C2, :N, β, C, indCinoAB, ())
        end
    else
        _native_contract!(α, A, CA, B, CB, β, C, oindA, cindA, oindB, cindB, indCinoAB)
    end
    return C
end

function isblascontractable(A::AbstractArray{T,N}, p1::IndexTuple, p2::IndexTuple,
        C::Symbol) where {T,N}

    T <: LinearAlgebra.BlasFloat || return false
    @unsafe_strided A isblascontractable(A, p1, p2, C)
end

function isblascontractable(A::AbstractStridedView{T,N}, p1::IndexTuple, p2::IndexTuple,
        C::Symbol) where {T,N}

    T <: LinearAlgebra.BlasFloat || return false
    strideA = i->stride(A, i)
    sizeA = i->size(A,i)

    canfuse1, d1, s1 = _canfuse(sizeA.(p1), strideA.(p1))
    canfuse2, d2, s2 = _canfuse(sizeA.(p2), strideA.(p2))

    if C == :D # destination
        return A.op == identity && canfuse1 && canfuse2 && s1 == 1
    elseif (C == :C && A.op == identity) || (C == :N && A.op == conj)# conjugated
        return canfuse1 && canfuse2 && s2 == 1
    else
        return canfuse1 && canfuse2 && (s1 == 1 || s2 == 1)
    end
end

_canfuse(::Dims{0}, ::Dims{0}) = true, 1, 1
_canfuse(dims::Dims{1}, strides::Dims{1}) = true, dims[1], strides[1]
function _canfuse(dims::Dims{N}, strides::Dims{N}) where {N}
    if dims[1] == 0
        return true, 0, 1
    elseif dims[1] == 1
        return _canfuse(Base.tail(dims), Base.tail(strides))
    else
        b, d, s = _canfuse(Base.tail(dims), Base.tail(strides))
        if b && (s == dims[1]*strides[1] || d == 1)
            dnew = dims[1]*d
            return true, dnew, (dnew == 0 || dnew == 1) ? 1 : strides[1]
        else
            return false, dims[1]*d, strides[1]
        end
    end
end
_trivtuple(t::NTuple{N}) where {N} = ntuple(identity, Val(N))

function _blas_contract!(α, A::AbstractArray{T}, CA, B::AbstractArray{T}, CB,
        β, C::AbstractArray{T}, oindA, cindA, oindB, cindB, oindAinC, oindBinC) where
        {T<:LinearAlgebra.BlasFloat}

    sizeA = i->size(A, i)
    sizeB = i->size(B, i)
    sizeC = i->size(C, i)

    csizeA = sizeA.(cindA)
    csizeB = sizeB.(cindB)
    osizeA = sizeA.(oindA)
    osizeB = sizeB.(oindB)

    @unsafe_strided A B C begin
        A2 = sreshape(permutedims(A, (oindA..., cindA...)), (prod(osizeA), prod(csizeA)))
        B2 = sreshape(permutedims(B, (cindB..., oindB...)), (prod(csizeB), prod(osizeB)))
        C2 = sreshape(permutedims(C, (oindAinC..., oindBinC...)),
                (prod(osizeA), prod(osizeB)))
        if CA == :N && CB == :N
            mul!(C2, A2, B2, α, β)
        elseif (CA == :C || CA == :A) && CB == :N
            mul!(C2, conj(A2), B2, α, β)
        elseif CA == :N && (CB == :C || CB == :A)
            mul!(C2, A2, conj(B2), α, β)
        elseif (CA == :C || CA == :A) && (CB == :C || CB == :A)
            mul!(C2, conj(A2), conj(B2), α, β)
        else
            throw(ArgumentError("unknown conjugation flag $CA and $CB"))
        end
    end
    return C
end

function _native_contract!(α, A::AbstractArray, CA::Symbol, B::AbstractArray, CB::Symbol,
        β, C::AbstractArray, oindA, cindA, oindB, cindB, indCinoAB)

    sizeA = i->size(A, i)
    sizeB = i->size(B, i)
    sizeC = i->size(C, i)

    csizeA = sizeA.(cindA)
    csizeB = sizeB.(cindB)
    osizeA = sizeA.(oindA)
    osizeB = sizeB.(oindB)

    ipC = TupleTools.invperm(indCinoAB)
    if CA == :N
        opA = identity
    elseif CA == :C
        opA = conj
    elseif CA == :A
        opA = adjoint
    else
        throw(ArgumentError("unknown conjugation flag $CA"))
    end
    if CB == :N
        opB = identity
    elseif CB == :C
        opB = conj
    elseif CB == :A
        opB = adjoint
    else
        throw(ArgumentError("unknown conjugation flag $CB"))
    end

    let opA = opA, opB = opB, α = α
        AS = sreshape(permutedims(StridedView(A), (oindA..., cindA...)),
            (osizeA..., one.(osizeB)..., csizeA...))
        BS = sreshape(permutedims(StridedView(B), (oindB..., cindB...)),
            (one.(osizeA)..., osizeB..., csizeB...))
        CS = sreshape(permutedims(StridedView(C), ipC),
            (osizeA..., osizeB..., one.(csizeA)...))
        tsize = (osizeA..., osizeB..., csizeA...)
        if α != 1
            op1 = (x,y) -> α*opA(x)*opB(y)
            if β == 0
                Strided._mapreducedim!(op1, +, zero, tsize, (CS, AS, BS))
            elseif β == 1
                Strided._mapreducedim!(op1, +, nothing, tsize, (CS, AS, BS))
            else
                Strided._mapreducedim!(op1, +, y->β*y, tsize, (CS, AS, BS))
            end
        else
            op2 = (x,y) -> opA(x)*opB(y)
            if β == 0
                Strided._mapreducedim!(op2, +, zero, tsize, (CS, AS, BS))
            elseif β == 1
                Strided._mapreducedim!(op2, +, nothing, tsize, (CS, AS, BS))
            else
                Strided._mapreducedim!(op2, +, y->β*y, tsize, (CS, AS, BS))
            end
        end
    end
    return C
end