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

using JuMP, Gurobi

function TEAVAR(env,
                edges,
                capacity,
                flows,
                demand,
                beta,
                k,
                T,
                Tf,
                scenarios,
                scenario_probs,
				outputfilename;
                explain=false,
                verbose=false,
                utilization=false,
                average=false)

    nedges = length(edges)
    nflows = length(flows)
    ntunnels = length(T)
    nscenarios = length(scenarios)
    p = scenario_probs

	println(Dates.format(now(), "HH:MM:SS"), ": #edges ", nedges, " #flows ", nflows, " #tunnels ", ntunnels, " #scenarios ", nscenarios)

    #CREATE TUNNEL SCENARIO MATRIX
    X  = ones(nscenarios,ntunnels)
    for s in 1:nscenarios
        for t in 1:ntunnels
            if size(T[t],1) == 0
                X[s,t] = 0
            else
                for e in 1:nedges
                    if scenarios[s][e] == 0
                        back_edge = findfirst(x -> x == (edges[e][2],edges[e][1]), edges)
                        if in(e, T[t]) || in(back_edge, T[t])
                        # if in(e, T[t])
                            X[s,t] = 0
                        end
                    end
                end
            end
        end
    end

	println(Dates.format(now(), "HH:MM:SS"), ": created tunnel scenario matrix")

    #CREATE TUNNEL EDGE MATRIX
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    println(Dates.format(now(), "HH:MM:SS"), ": created tunnel edge matrix")

    model = Model(solver=GurobiSolver(env, OutputFlag=1))
    # flow per commodity per path variables
    @variable(model, a[1:nflows, 1:k] >= 0, basename="a", category=:SemiCont)
    # alpha variable
    @variable(model, alpha >= 0, basename="alpha", category=:SemiCont)
    # maximum flow lost in that scenario
    @variable(model, umax[1:nscenarios] >= 0, basename="umax")
    # flow lost per commod per scenario
    @variable(model, u[1:nscenarios, 1:nflows] >= 0, basename="u")

    # capacity constraints for final flow assigned to "a" variables
    for e in 1:nedges
        @constraint(model, sum(a[f,t] * L[Tf[f][t],e] for f in 1:nflows, t in 1:size(Tf[f],1)) <= capacity[e])
    end

      # FLOW LEVEL LOSS
      @expression(model, satisfied[s=1:nscenarios, f=1:nflows], sum(a[f,t] * X[s,Tf[f][t]] for t in 1:size(Tf[f],1)) / demand[f])

      for s in 1:nscenarios
          for f in 1:nflows
              # @constraint(model, (demand[f] - sum(a[f,t] * X[s,Tf[f][t]] for t in 1:size(Tf[f],1))) / demand[f] <= u[s,f])
              @constraint(model, u[s,f] >= 1 - satisfied[s,f])
          end
      end

      for s in 1:nscenarios
          if average
              @constraint(model, umax[s] + alpha >= (sum(u[s,f] for f in 1:nflows)) / nflows)
              # @constraint(model, umax[s] + alpha >= avg_loss[s])
          else
              for f in 1:nflows
                  @constraint(model, umax[s] + alpha >= u[s,f])
              end
          end
      end
      @objective(model, Min, alpha + (1 / (1 - beta)) * sum((p[s] * umax[s] for s in 1:nscenarios)))

    println(Dates.format(now(), "HH:MM:SS"), ": ready to solve")

    solve(model)

    println(Dates.format(now(), "HH:MM:SS"), ": solver finished; explaining")

    if (explain)
        println("Runtime: ", getsolvetime(model))
	println("beta: ", beta)
	println("#flows: ", nflows)
	println("#edges: ", nedges)
	println("#tunnels: ", ntunnels)
	println("#demands: ", nflows, " total demand=", sum(demand[f] for f in 1:nflows))
	println("#scenarios: ", nscenarios, " total_prob= ", sum(p[s] for s in 1:nscenarios))
	println("total_scenario_prob= ", sum(p[s] for s in 1:nscenarios))

	result_u = getvalue(u)
	result_umax = getvalue(umax)
	result_satisfied = getvalue(satisfied)
	result_alpha = getvalue(alpha)
	result_a = getvalue(a)

	# compute allocations and weights
	flow_allocation_1 = Array{Float64}(undef, nflows)
	numtunnel_perflow = 0
	for f in 1:nflows
		numtunnel_perflow = max(numtunnel_perflow, size(Tf[f], 1))
	end
	flow_tunnel_weights = Matrix{Float64}(undef, nflows, numtunnel_perflow)
	all_flow_alloc = 0

	for f in 1:nflows
		flow_totalalloc = 0
		for t in 1:size(Tf[f], 1)
			flow_totalalloc += result_a[f,t]
		end

		all_flow_alloc += flow_totalalloc
		flow_allocation_1[f] = (1 - result_alpha) * min(flow_totalalloc, demand[f])
		for t in 1:size(Tf[f], 1)
			weight = 0
			if flow_totalalloc > 0
				weight = result_a[f, t]/ flow_totalalloc
			end
			flow_tunnel_weights[f, t] = weight
		end
	end
	flow_allocation_1[flow_allocation_1 .< 0] .= 0

	println("totalalloc= ", all_flow_alloc, " netalloc1= ", sum(flow_allocation_1[f] for f in 1:nflows))

	output = open(outputfilename, "w")
	write(output, "---flow allocation with alpha---\n")
	writedlm(output, flow_allocation_1, ',')
	write(output, "---weights per flow tunnel---\n")
	writedlm(output, flow_tunnel_weights, ',')
	flush(output)
	close(output)

        printResults(getobjectivevalue(model), result_alpha, result_a, result_u, result_umax, edges, scenarios, T, Tf, L, capacity, p, demand, verbose=verbose, utilization=utilization)
    end

    return (result_alpha, getobjectivevalue(model), result_a, result_umax, getsolvetime(model))
end


text:
module DEC2D

using StaticArrays

struct Triangulation
    vertices::Vector{SVector{2,Float64}}
    edges::Vector{SVector{2,Int}}
    triangles::Vector{SVector{3,Int}}
end

struct DualTriangulation
    n::Tuple{Int,Int,Int} # #{dual triangles} , #{midpoints}, #{exterior vertices}
    dtriangles::Vector{SVector{2,Float64}} # dual triangles (vertices)
    dedges::Vector{SVector{2,Int}} # dual edges (edges)
    dvertices::Vector{Vector{Int}} # dual vertices (polygons)
end

export DECMesh
struct DECMesh
    tri::Triangulation
    dtri::DualTriangulation
end

include("geometry.jl")
include("triangulation.jl")
include("diffop.jl")

end # module

ext: jl


text:
##
"Compute the square difference between a model spreading and an estimate."
function squarediff(model::FrequencyDirectionSpectra.FDSpectralModel, fdest::Union{FrequencyDirectionSpectra.MEM,FrequencyDirectionSpectra.MLM})
    return squarediff(model, fdest.D, fdest.Ω, fdest.Φ)
end
"Compute the square difference between a model spreading and an estimate seperated into ordinates and frequencies."
function squarediff(model::FrequencyDirectionSpectra.FDSpectralModel, D, Ω, Φ)
    return sum((FrequencyDirectionSpectra.spreading(model, Ω[i], Φ[j]) - D[i, j])^2 for i ∈ 1:size(D, 1), j ∈ 1:size(D, 2))
end
"Compute the square difference between a model sdf and an estimate."
function squarediff(model::FrequencyDirectionSpectra.FDSpectralModel,B::BartlettPeriodogram)
    return squarediff(model, B.Ω, B.ordinates)
end
"Compute the square difference between a model sdf and an estimate seperated into ordinates and frequencies."
function squarediff(model::FrequencyDirectionSpectra.FDSpectralModel, Ω, ordinates)
    return sum((FrequencyDirectionSpectra.sdf(model, ω) - I)^2 for (ω, I) ∈ zip(Ω, ordinates))
end

##
function LSfit(
    model,
    ts,
    Δ,
    θ₀;
    fdest = FrequencyDirectionSpectra.MLM,
    ωₗ = 0,
    ωᵤ = Inf,
    window = 128,
    ndirection = 100,
)
    B = BartlettPeriodogram(ts, Δ, window)[end÷2+1:end]
    Bused = B[ωₗ.<=B.Ω.<=ωᵤ]
    S = fdest(Bused.ordinates, Bused.Ω, range(-π, π, length = ndirection))
    Bz = BartlettPeriodogram(ts[:, 1], Δ, window)[end÷2+1:end][ωₗ.<=B.Ω.<=ωᵤ]
    res1 = optimize(θ -> marginalbounds(θ,ωₗ,ωᵤ) ? squarediff(model(θ), Bz) : Inf, θ₀[1:4])
    if Optim.converged(res1)
        res2 = optimize(θ -> spreadingbounds(θ) ? squarediff(model([res1.minimizer; θ]), S) : Inf, θ₀[5:end])
        if Optim.converged(res2)
            return (true, [res1.minimizer; res2.minimizer])
        end
    end
    return (false, θ₀)
end

##
function marginalbounds(θ,ωₗ,ωᵤ)
    return 4>θ[1]>0 && ωᵤ>θ[2]>ωₗ && 20>θ[3]>1 && 8>θ[4]>3
end
function spreadingbounds(θ)
    return 2π>θ[2]>0 && θ[3]>0 && θ[4]>0 && θ[5]>0
end


text:
module Turing

##############
# Dependency #
########################################################################
# NOTE: when using anything from external packages,                    #
#       let's keep the practice of explictly writing Package.something #
#       to indicate that's not implemented inside Turing.jl            #
########################################################################

using Requires
using Distributions
using ForwardDiff

using LinearAlgebra
using ProgressMeter
using Markdown

#  @init @require Stan="682df890-35be-576f-97d0-3d8c8b33a550" begin
  using Stan
  import Stan: Adapt, Hmc
#  end
import Base: ~, convert, promote_rule, rand, getindex, setindex!
import Distributions: sample
import ForwardDiff: gradient
using Flux: Tracker
import MCMCChain: AbstractChains, Chains

##############################
# Global variables/constants #
##############################

const ADBACKEND = Ref(:reverse_diff)
setadbackend(backend_sym) = begin
  @assert backend_sym == :forward_diff || backend_sym == :reverse_diff
  ADBACKEND[] = backend_sym
end

const ADSAFE = Ref(false)
setadsafe(switch::Bool) = begin
  @info("[Turing]: global ADSAFE is set as $switch")
  ADSAFE[] = switch
end

const CHUNKSIZE = Ref(0) # default chunksize used by AD

setchunksize(chunk_size::Int) = begin
  if ~(CHUNKSIZE[] == chunk_size)
    @info("[Turing]: AD chunk size is set as $chunk_size")
    CHUNKSIZE[] = chunk_size
    global SEEDS = ForwardDiff.construct_seeds(ForwardDiff.Partials{chunk_size,Float64}) # pre-alloced dual parts
  end
end

setchunksize(40)

const PROGRESS = Ref(true)
turnprogress(switch::Bool) = begin
  @info("[Turing]: global PROGRESS is set as $switch")
  PROGRESS[] = switch
end

# Constants for caching
const CACHERESET  = 0b00
const CACHEIDCS   = 0b10
const CACHERANGES = 0b01

#######################
# Sampler abstraction #
#######################

abstract type InferenceAlgorithm end
abstract type Hamiltonian <: InferenceAlgorithm end

"""
    Sampler{T}

Generic interface for implementing inference algorithms.
An implementation of an algorithm should include the following:

1. A type specifying the algorithm and its parameters, derived from InferenceAlgorithm
2. A method of `sample` function that produces results of inference, which is where actual inference happens.

Turing translates models to chunks that call the modelling functions at specified points. The dispatch is based on the value of a `sampler` variable. To include a new inference algorithm implements the requirements mentioned above in a separate file,
then include that file at the end of this one.
"""
mutable struct Sampler{T<:InferenceAlgorithm}
  alg   ::  T
  info  ::  Dict{Symbol, Any}         # sampler infomation
end

include("../deps/deps.jl"); check_deps();
include("helper.jl")
include("transform.jl")
include("core/varinfo.jl")  # core internal variable container
include("trace/trace.jl")   # to run probabilistic programs as tasks

using Turing.Traces
using Turing.VarReplay

###########
# Exports #
###########

# Turing essentials - modelling macros and inference algorithms
export @model, @~, @VarName                   # modelling
export MH, Gibbs                              # classic sampling
export HMC, SGLD, SGHMC, HMCDA, NUTS          # Hamiltonian-like sampling
export IS, SMC, CSMC, PG, PIMH, PMMH, IPMCMC  # particle-based sampling
export sample, setchunksize, resume           # inference
export auto_tune_chunk_size!, setadbackend, setadsafe # helper
export turnprogress  # debugging

# Turing-safe data structures and associated functions
export TArray, tzeros, localcopy, IArray

export @sym_str

export UnivariateGMM2, Flat, FlatPos

##################
# Inference code #
##################

include("core/util.jl")         # utility functions
include("core/compiler.jl")     # compiler
include("core/container.jl")    # particle container
include("core/io.jl")           # I/O
include("samplers/sampler.jl")  # samplers
include("core/ad.jl")           # Automatic Differentiation

end


text:
export ispercolation!, clustersize!, clustersizefreq!, clusternumber!, average_clustersize!, strength!

"""
ispercolation!(latt::TwoDLattice)

check whether the system is percolated or not.
"""
function ispercolation!(latt::TwoDLattice)
    if latt.lattice_properties.ispercolationcheck
        return latt.lattice_properties.ispercolation
    end
    latt.lattice_properties.ispercolationcheck = true

    if !latt.lattice_properties.islabeled
        label_components!(latt)
        latt.lattice_properties.islabeled = true
    end

    labeled_site = latt.lattice_properties.labeled_lattice_sites

    tmpset = (Set(labeled_site[1,:]) ∩ Set(labeled_site[end,:])) ∩ (Set(labeled_site[:,1]) ∩ Set(labeled_site[:,end]))
    setdiff!(tmpset, 0)
    latt.lattice_properties.ispercolation = ifelse(isempty(tmpset), false, true)
    if latt.lattice_properties.ispercolation
        percoclusterlabel = pop!(tmpset)
        latt.lattice_properties.percoclusterlabel = percoclusterlabel
        latt.lattice_properties.percoclustersize = sum(latt.lattice_properties.labeled_lattice_sites .== percoclusterlabel)
    end

    return latt.lattice_properties.ispercolation
end

"""
clustersize!(latt::Lattice)

check each cluster size.
"""
function clustersize!(latt::Lattice)
    if latt.lattice_properties.isclustersizecheck
        return latt.lattice_properties.clustersize
    end
    latt.lattice_properties.isclustersizecheck = true

    if !latt.lattice_properties.islabeled
        label_components!(latt)
    end

    clustersize = zeros(Int, latt.lattice_properties.nclusters)
    for i in 1:latt.lattice_config.linear_size^2
        sitei = latt.lattice_properties.labeled_lattice_sites[i]
        if sitei != 0 && sitei != latt.lattice_properties.percoclusterlabel
            clustersize[latt.lattice_properties.labeled_lattice_sites[i]] += 1
        end
    end

    latt.lattice_properties.clustersize = clustersize

    return latt.lattice_properties.clustersize
end

"""
clustersizefreq!(latt::Lattice)

check the cluster size distribution.
"""
function clustersizefreq!(latt::Lattice)
    if latt.lattice_properties.isclustersizefreqcheck
        return latt.lattice_properties.clustersizefreq
    end

    if !latt.lattice_properties.isclustersizecheck
        clustersize!(latt)
    end

    latt.lattice_properties.clustersizefreq = delete!(countmap(latt.lattice_properties.clustersize), 0)
    latt.lattice_properties.isclustersizefreqcheck = true

    return latt.lattice_properties.clustersizefreq
end

"""
clusternumber!(latt::Lattice)

Cluster Number \$n_s\$ is the probability that a occupied site belonging to s-cluster.
Formal definition of cluster number is defined by exclunding the percolation cluster.
But in this calculation, it is included.
"""
function clusternumber!(latt::Lattice)
    if latt.lattice_properties.isclusternumbercheck
        return latt.lattice_properties.clusternumber
    end

    clustersizefreq!(latt)
    clusternumber = Dict(key => latt.lattice_properties.clustersizefreq[key] / latt.lattice_config.linear_size^2
                         for key in keys(latt.lattice_properties.clustersizefreq))
    latt.lattice_properties.clusternumber = clusternumber
    latt.lattice_properties.isclusternumbercheck = true

    return latt.lattice_properties.clusternumber
end

"""
average_clustersize!(latt::Lattice)

Average cluster size \$S\$ is defined by
\$S = \\frac{\\sum_{s=1} s^2 n_s}{p}\$.
"""
function average_clustersize!(latt::Lattice)
    if latt.lattice_properties.isaverage_clustersizecheck
        return latt.lattice_properties.average_clustersize
    end

    if !latt.lattice_properties.isclusternumbercheck
        clusternumber!(latt)
    end

    clusternumberdict = latt.lattice_properties.clusternumber
    average_clustersize = sum(collect(keys(clusternumberdict)).^2 .* collect(values(clusternumberdict))
                             ) / latt.lattice_config.p
    latt.lattice_properties.average_clustersize = average_clustersize

    return latt.lattice_properties.average_clustersize
end

"""
strength!(latt::Lattice)

Strength \$P(p)\$ is the probability that an arbitrary site belongs to the percolating cluster.

\$P(p) = \\frac{\\text{# of sites in percolating cluster with occupied probability}~ p}{L^2}\$  where \$L\$ is the system linear size
"""
function strength!(latt::Lattice)
    if latt.lattice_properties.isstrengthcheck
        return latt.lattice_properties.strength
    end

    if ispercolation!(latt)
       latt.lattice_properties.strength = latt.lattice_properties.percoclustersize / latt.lattice_config.linear_size^2
    end

    return latt.lattice_properties.strength
end


text:
# ---
# title: 413. Arithmetic Slices
# id: problem413
# author: Indigo
# date: 2021-06-20
# difficulty: Medium
# categories: Math, Dynamic Programming
# link: <https://leetcode.com/problems/arithmetic-slices/description/>
# hidden: true
# ---
#
# A sequence of numbers is called arithmetic if it consists of at least three
# elements and if the difference between any two consecutive elements is the
# same.
#
# For example, these are arithmetic sequences:
#
#
#
#     1, 3, 5, 7, 9
#     7, 7, 7, 7
#     3, -1, -5, -9
#
# The following sequence is not arithmetic.
#
#
#
#     1, 1, 2, 5, 7
#
#
#
# A zero-indexed array A consisting of N numbers is given. A slice of that array
# is any pair of integers (P, Q) such that 0 <= P < Q < N.
#
# A slice (P, Q) of the array A is called arithmetic if the sequence:
# A[P], A[P \+ 1], ..., A[Q - 1], A[Q] is arithmetic. In particular, this means
# that P + 1 < Q.
#
# The function should return the number of arithmetic slices in the array A.
#
#
#
# **Example:**
#
#
#
#     A = [1, 2, 3, 4]
#
#     return: 3, for 3 arithmetic slices in A: [1, 2, 3], [2, 3, 4] and [1, 2, 3, 4] itself.
#
#
#
## @lc code=start
using LeetCode

function number_of_arithmetic_slices(nums::Vector{Int})
    dp, res = 0, 0
    for i in 3:length(nums)
        if nums[i] + nums[i - 2] == (nums[i - 1] << 1)
            res += (dp += 1)
        else
            dp = 0
        end
    end
    return res
end
## @lc code=end


text:
abstract type AbstractCacheStrategy end

struct DefaultCache <: AbstractCacheStrategy end
struct SolutionCache <: AbstractCacheStrategy end
struct EndpointCache <: AbstractCacheStrategy end

abstract type AbstractRenderCache{M,T} end

struct SolutionRenderCache{M,T,G} <: AbstractRenderCache{M,T}
    # metric
    m::M

    # max time
    max_time::T

    # size information
    height::Int
    width::Int

    # geodesics themselves in 2d array
    geodesics::AbstractMatrix{G}

    function SolutionRenderCache(
        m::AbstractMetricParams{T},
        max_time::T,
        height,
        width,
        cache::AbstractVector{SciMLBase.EnsembleSolution{T,N,Vector{O}}},
    ) where {T,N,O}

        geodesics = Matrix{O}(undef, (height, width))

        # populate store
        for (col, simsol) in enumerate(cache)
            for (row, sol) in enumerate(simsol)
                geodesics[col, row] = sol
            end
        end

        # return instance
        new{typeof(m),T,O}(m, max_time, height, width, geodesics)
    end
end

struct EndpointRenderCache{M,T,P} <: AbstractRenderCache{M,T}
    # metric
    m::M

    # max time
    max_time::T

    # size information
    height::Int
    width::Int

    points::Matrix{P}
end

function Base.show(io::IO, cache::AbstractRenderCache{M}) where {M}
    repr = "$(Base.typename(typeof(cache)).name){$M} (dimensions $(cache.height)x$(cache.width))"
    write(io, repr)
end


text:
function inner_product(a1, a2)
  somme = 0
  for (word1, count1) in a1
    for (word2, count2) in a2
      if word2 == word1
        somme += count1 * count2
        break
      end
    end
  end
  return somme
end

function len(a)
  return √(innter_product(a, a))
end

function θ(a1, a2)
  cosθ = inner_product(a1, a2) / √(inner_product(a1, a1) * inner_product(a2, a2))
  return acos(cosθ)
end


text:
function ContinueOnEvent()
    return ContinueOnEvent(())
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function event_occurred(obj::ContinueOnEvent, arg0::SpacecraftState, arg1::EventDetector, arg2::jboolean)
    return jcall(obj, "eventOccurred", Action, (SpacecraftState, EventDetector, jboolean), arg0, arg1, arg2)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function init(obj::EventHandler, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function reset_state(obj::EventHandler, arg0::EventDetector, arg1::SpacecraftState)
    return jcall(obj, "resetState", SpacecraftState, (EventDetector, SpacecraftState), arg0, arg1)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
#
# Various data containers.
#

const OneTo = Base.OneTo

include("data/entity.jl")
include("data/record.jl")
include("data/column.jl")
include("data/datavector.jl")
include("data/dictvector.jl")

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 LunoAPI

using HTTP, Base64, JSON

include("api.jl")

export set_auth_token, get_auth_token, get_price, get_balance, stop_all_orders,
    execute_market_order, post_limit_order

end # module


text:
export WhithamGreenNaghdi

"""
    WhithamGreenNaghdi(param;kwargs)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
the modified Green-Naghdi model proposed by V. Duchêne, S. Israwi and R. Talhouk.

# Argument
`param` is of type `NamedTuple` and must contain
- dimensionless parameters `ϵ` (nonlinearity) and `μ` (dispersion);
- numerical parameters to construct the mesh of collocation points as `mesh = Mesh(param)`.

## Keywords
- `SGN`: if `true` computes the Serre-Green-Naghdi (SGN) instead of Whitham-Green-Naghdi (WGN) system (default is `false`);
- `iterative`: solve the elliptic problem through GMRES if `true`, LU decomposition if `false` (default is `true`);
- `precond`: use a (left) preconditioner for GMRES if `true` (default), choose `precond` as the preconditioner if provided;
- `gtol`: relative tolerance of the GMRES algorithm (default is `1e-14`);
- `restart`: the corresponding option of the GMRES algorithm (default is `100`);
- `maxiter`: the corresponding option of GMRES (default is `nothing`);
- `ktol`: tolerance of the Krasny filter (default is `0`, i.e. no filtering);
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `verbose`: prints information if `true` (default is `true`).

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!` and in particular
1. a function `WhithamGreenNaghdi.f!` to be called in the time-integration solver;
2. a function `WhithamGreenNaghdi.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `WhithamGreenNaghdi.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v)`, where

    - `η` is the surface deformation;
    - `v` is the derivative of the trace of the velocity potential;
4. additionally, a handy function `WhithamGreenNaghdi.mapfrofull` which from data matrix returns the Tuple of real vectors `(η,v,u)`, where

	- `u` corresponds to the layer-averaged velocity.

"""
mutable struct WhithamGreenNaghdi <: AbstractModel

	label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	mapfrofull	:: Function
	param	:: NamedTuple
	kwargs  :: NamedTuple

    function WhithamGreenNaghdi(param::NamedTuple;SGN=false,dealias=0,ktol=0,iterate=true,gtol=1e-14,precond=true,restart=nothing,maxiter=nothing,verbose=true)
		if SGN == true
			label = string("Serre-Green-Naghdi")
		else
			label = string("Whitham-Green-Naghdi")
		end
		if verbose @info string("model ",label) end

		kwargs = (iterate=iterate,SGN=SGN,dealias=dealias,ktol=ktol,gtol=gtol,precond=precond,verbose=verbose)
		μ 	= param.μ
		ϵ 	= param.ϵ
		mesh = Mesh(param)
		param = ( ϵ = ϵ, μ = μ, xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N )

		k = mesh.k
		x 	= mesh.x
		x₀ = mesh.x[1]

		∂ₓ	=  1im * mesh.k
        if SGN == true
			F₁ = 1 ./(1 .+ μ/3*k.^2)
			F₀ = sqrt(μ)*∂ₓ
	    else
			F₁ 	= tanh.(sqrt(μ)*abs.(k))./(sqrt(μ)*abs.(k))
			F₁[1] 	= 1
			F₀ = 1im * sqrt.(3*(1 ./F₁ .- 1)).*sign.(k)
		end
		if precond == true
			#Precond = Diagonal( 1 .+ μ/3*k.^2 )
			Precond = Diagonal( 1 ./  F₁ )
		elseif precond == false
			Precond = Diagonal( ones(size(k)) )
		else
			Precond = precond
		end
		K = mesh.kmax * (1-dealias/(2+dealias))
		Π⅔ 	= abs.(mesh.k) .<= K # Dealiasing low-pass filter
		if dealias == 0
			if verbose @info "no dealiasing" end
			Π⅔ 	= ones(size(mesh.k))
		elseif verbose
			@info string("dealiasing : spectral scheme for power ", dealias + 1," nonlinearity ")
		end
		if iterate == true && verbose
			@info "elliptic problem solved with GMRES method"
		elseif verbose
			@info "elliptic problem solved with LU decomposition"
		end
		FFT = exp.(-1im*k*(x.-x₀)');
        IFFT = exp.(1im*k*(x.-x₀)')/length(x);
		M₀ = IFFT * Diagonal( F₀ .* Π⅔)* FFT
		IFFTF₀ = IFFT * Diagonal( F₀ .* Π⅔)
        Id = Diagonal(ones(size(x)));
		h = zeros(Complex{Float64}, mesh.N)
		u, fftv, fftη, fftu, hdu = (similar(h),).*ones(5)
		L = similar(FFT)
		if maxiter == nothing maxiter = mesh.N end
		if restart == nothing restart = min(20,mesh.N) end

		# Evolution equations are ∂t U = f!(U)
		function f!(U)
			fftη .= U[:,1]
			h .= 1 .+ ϵ*ifft(fftη)
			fftv .= U[:,2]
			if iterate == false
				L .= Id - 1/3 * Diagonal(Π⅔) * FFT * Diagonal( 1 ./h ) * M₀ * Diagonal( h.^3 ) * IFFTF₀
				fftu .= L \ fftv
			elseif iterate == true
		        function LL(hatu)
		            hatu- 1/3 *Π⅔.*fft( 1 ./h .* ifft( F₀ .* Π⅔.*fft( h.^3 .* ifft( F₀ .* Π⅔.* hatu ) ) ) )
				end
				fftu .= gmres( LinearMap(LL, length(h); issymmetric=false, ismutating=false) , fftv ;
						restart = restart, maxiter = maxiter, Pl = Precond, reltol = gtol )
			end
			u .= ifft(fftu)
			hdu .= h .* ifft(Π⅔.*F₀.*fftu)
		   	U[:,1] .= -∂ₓ.*Π⅔.*(fftu .+ ϵ * fft(ifft(fftη) .* u))
			U[:,2] .= -∂ₓ.*Π⅔.*(fftη .+ ϵ * fft( u.*ifft(fftv)
							.- 1/2 * u.^2 .- 1/2 * hdu.^2 ) )
			U[abs.(U).< ktol ].=0
		end

		# Discrete Fourier transform with, possibly, dealiasing and Krasny filter.
		function mapto(data::InitialData)
			U = [Π⅔ .* fft(data.η(x)) Π⅔ .*fft(data.v(x))]
			U[abs.(U).< ktol ].=0
			return U
		end

		# Return `(η,v)`, where
		# - `η` is the surface deformation;
		# - `v` is the derivative of the trace of the velocity potential.
		# Inverse Fourier transform and takes the real part.
		function mapfro(U)
			real(ifft(U[:,1])),real(ifft(U[:,2]))
		end
		# Returns `(η,v,u)`, where
		# - `η` is the surface deformation;
		# - `v` is the derivative of the trace of the velocity potential;
		# - `u` corresponds to the layer-averaged velocity.
		# Inverse Fourier transform and take the real part, plus solves the costly elliptic problem for `u`.
		function mapfrofull(U)
				fftη .= U[:,1]
			   	h .= 1 .+ ϵ*ifft(fftη)
				L .= Id - 1/3 * Diagonal(Π⅔) * FFT * Diagonal( 1 ./h ) * M₀ * Diagonal( h.^3 ) * IFFTF₀

				   real(ifft(U[:,1])),real(ifft(U[:,2])),real(ifft(L \ U[:,2]))
		end

        new(label, f!, mapto, mapfro, mapfrofull, param, kwargs)
    end
end


text:
#Pedro F Felzenszwalb and Daniel P Huttenlocher. Efficient graph-based image segmentation.International Journal of Computer Vision, 59(2):167–181, 2004
# Adapted to use RGB
module CDNGraph
export config
import RGBDSegmentation:
    CDNImage,
    clusterize
import ImageSegmentation:
    felzenszwalb,
    ImageEdge,
    meantype,
    SegmentedImage,
    _colon
using LinearAlgebra
using DataStructures

"""
    edge = ImageEdge(index1, index2, weight)
Construct an edge in a Region Adjacency Graph. `index1` and `index2` are the integers corresponding to individual pixels/voxels (in the sense of linear indexing via `sub2ind`), and `weight` is the edge weight (measures the dissimilarity between pixels/voxels).
"""
struct ImageEdgeCDN
    index1::Int
    index2::Int
    weightRGB::Float64
    weightXYZ::Float64
    weightN::Float64
end
function felzenszwalb(edges::Array{ImageEdgeCDN}, num_vertices::Int, k::Real, min_size::Int = 0)

    num_edges = length(edges)
    G = IntDisjointSets(num_vertices)
    set_size = ones(num_vertices)
    thresholdRGB = fill(convert(Float64,k), num_vertices)
    thresholdXYZ = fill(convert(Float64,k), num_vertices)
    thresholdN   = fill(convert(Float64,k), num_vertices)

    sort!(edges, lt = (x,y)->(x.weightRGB<y.weightRGB))

    for edge in edges
        (wRGB, wXYZ, wN) = (edge.weightRGB, edge.weightXYZ, edge.weightN)
        a = find_root(G, edge.index1)
        b = find_root(G, edge.index2)
        if a!=b
            if (wRGB <= min(thresholdRGB[a], thresholdRGB[b])) && (wXYZ <= min(thresholdXYZ[a], thresholdXYZ[b]))  && (wN <= min(thresholdN[a], thresholdN[b]))
                merged_root = union!(G, a, b)
                set_size[merged_root] = set_size[a] + set_size[b]
                thresholdRGB[merged_root] = wRGB + k/set_size[merged_root]
                thresholdXYZ[merged_root] = wXYZ + k/set_size[merged_root]
                thresholdN[merged_root] = wN + k/set_size[merged_root]
            end
        end
    end

    #merge small segments
    for edge in edges
        a = find_root(G, edge.index1)
        b = find_root(G, edge.index2)
        if a!=b && (set_size[a] < min_size || set_size[b] < min_size)
            union!(G, a, b)
        end
    end

    segments = OrderedSet()
    for i in 1:num_vertices
        push!(segments, find_root(G, i))
    end

    num_sets = length(segments)
    segments2index = Dict{Int, Int}()
    for (i, s) in enumerate(segments)
        segments2index[s]=i
    end

    index_map = Array{Integer}(undef, num_vertices)
    for i in 1:num_vertices
        index_map[i] = segments2index[find_root(G, i)]
    end

    return index_map, num_sets
end
function felzenszwalb(img::CDNImage{T}, k::Real, edge_weight::Function, min_size::Int = 0) where T<:Real
    _, rows, cols = size(img)
    num_vertices = rows*cols
    num_edges = 4*rows*cols - 3*rows - 3*cols + 2
    edges = Array{ImageEdgeCDN}(undef, num_edges)

    R = CartesianIndices((size(img,2),size(img,3)))
    I1, Iend = first(R), last(R)
    num = 1
    for I in R
        for J in CartesianIndices(_colon(max(I1, I-I1), min(Iend, I+I1)))
            if I >= J
                continue
            end
            edges[num] = ImageEdgeCDN(
                (I[2]-1)*rows+I[1],
                (J[2]-1)*rows+J[1],
                edge_weight(img[:,I[1],I[2]],img[:,J[1],J[2]])...)
            num += 1
        end
    end

    index_map, num_segments = felzenszwalb(edges, num_vertices, k, min_size)

    labels = zeros(Integer, (size(img,2),size(img,3)))
    for j=1:size(img, 3), i=1:size(img, 2)
       labels[i, j] = index_map[(j-1)*rows+i]
    end
    return SegmentedImage(img, labels)
end
function weight(a::Vector, b::Vector)
    return (norm(a[1:3]-b[1:3]),norm(a[4:6]-b[4:6]),acos(clamp(dot(a[7:9], b[7:9]),0,1)))
end
function felzenszwalb(img::CDNImage, k::Real, min_size::Int=0)
    return felzenszwalb(img,k, weight, min_size)
end
struct Config
   min_size::Integer
end
function Config(;min_size::Integer)
    return Config(min_size)
end
function clusterize(cfg::Config, img::CDNImage, k::Integer)
    return felzenszwalb(img, k, cfg.min_size)
end
end


text:
function displacement(r::AbstractMatrix, δ::Int)
    sqrt((r[1+δ, 2] - r[1, 2])^2 + (r[1+δ, 1] - r[1, 1])^2)
end

function displacement(r::AbstractMatrix, t::Int, δ::Int)
    sqrt((r[t+δ, 2] - r[t, 2])^2 + (r[t+δ, 1] - r[t, 1])^2)
end

squared_displacement(r::AbstractMatrix, t::Int, δ::Int) = displacement(r, t, δ)^2

function extract(df::DataFrame, n::Int, id::Symbol, x::Symbol)
    Vector(@where(df, cols(id) .== n)[!, x])
    # Vector{Float64}(df[df[:, id] .== n, x])
end

function extract(df::DataFrame, n::Int, id::Symbol, label::Vector{Symbol})
    Matrix(@where(df, cols(id) .== n)[!, label])
    # df[df[!, id] .== n, label]
end

max_stars_repo_name: SosUts/DiffusionMoments

max_issues_repo_name: SosUts/DiffusionMoments

max_forks_repo_name: SosUts/DiffusionMoments


text:
function uniquerows(a)

ind = findall(abs.(a[:]).<1e-10);
a[ind] .= 0;

n = size(a,1);
ia = unique(i -> a[i,:], 1:n);
c = a[ia,:];
ic = [findfirst(i->all(j->abs(a[k,j] - c[i,j])<1e-8,1:size(c,2)),1:size(c,1)) for k = 1:size(a,1)];

# d,tm = sortrows(c);
# e = d[2:end,:]-d[1:end-1,:];
# e = sum(e.*e,dims=2);
# i = argmin(e[:]);
# display(d[i-1:i+1,:])
# display(minimum(e[:]))

# display(size(c))
# display(size(a))
# display(size(ic))
#
# ent = unique(ic[:]);
# ndof = length(ent);
# display([minimum(ic) maximum(ic) ndof])

# ii = setdiff(1:n,ia);
# nc = length(ia);
# ic = zeros(Int,n);
# ic[ia] = 1:nc;
#
# #ic[ii] = xiny(a[ii,:],c,1);
#
# # nb = length(ii);
# # for k = 1:nb
# #     for q = 1:nc
# #         if maximum(abs.(a[ii[k],:]-c[q,:]))<1e-7
# #             ic[ii[k]] = q;
# #             break;
# #         end
# #     end
# # end
#
# while (true)
#     tm = a[ii,:];
#     display(size(tm))
#     if length(ii)==1
#         ij = 1;
#         b = reshape(tm[1,:],1,length(tm));
#     else
#         ij = unique(ix -> tm[ix,:], 1:size(tm,1));
#         b = tm[ij,:];
#     end
#     # dsiplay(ii)
#     # display(size(a))
#     # display(size(tm))
#
#     nb = size(b,1);
#     ind = zeros(Int,nb);
#     for i = 1:nc
#         if maximum(abs.(b[1,:]-c[i,:]))<1e-7
#             ind[1] = i;
#             break;
#         end
#     end
#     m = ind[1];
#     for k = 2:nb
#         found = 0;
#         for q = m:nc
#             if maximum(abs.(b[k,:]-c[q,:]))<1e-7
#                 ind[k] = q;
#                 m = q;
#                 found = 1;
#                 break;
#             end
#         end
#         if found==0
#             for q = 1:m-1
#                 if maximum(abs.(b[k,:]-c[q,:]))<1e-7
#                     ind[k] = q;
#                     m = q;
#                     break;
#                 end
#             end
#         end
#     end
#     if length(ii)==1
#         ic[ii] = ind;
#     else
#         ic[ii[ij]] = ind;
#     end
#
#     i1 = setdiff(1:length(ii),ij);
#     if length(i1)==0
#         break;
#     else
#         ii = ii[i1];
#     end
# end

# ii = setdiff(1:n,ia);
# nc = length(ia);
# ic = zeros(Int,n);
# ic[ia] = 1:nc;
#
# b,ij = sortrows(a[ii,:]);
# m = 0;
# for i = 1:nc
#     if maximum(abs.(b[1,:]-c[i,:]))<1e-7
#         m = i;
#         break;
#     end
# end
# ic[ii[ij[1]]] = m;
#
# for k = 2:length(ii)
#     if maximum(abs.(b[k,:]-b[k-1,:]))<1e-7
#         ic[ii[ij[k]]] = ic[ii[ij[k-1]]];
#     else
#         for q = 1:nc
#             if maximum(abs.(b[k,:]-c[q,:]))<1e-7
#                 ic[ii[ij[k]]] = q;
#                 m = q;
#                 break;
#             end
#         end
#     end
# end

return c,ia,ic

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
include("../src/hypergraph-affinity-functions.jl")

names = ["Clothing, Shoes & Accessories"
"Electronics and Office"
"Home, Furniture & Appliances"
"Home Improvement & Patio"
"Baby"
"Toys, Games, and Video Games"
"Food, Household & Pets"
"Pharmacy, Health & Beauty"
"Sports, Fitness & Outdoors"
"Auto, Tires & Industrial"
"Other"]

dataset = "walmart-trips"
r = 25
H, labels = read_hypergraph_data(dataset,r)

## Categories 1 and 7 make up over 50% of the dataset.
# At a high level, this suggests that a lot of shopping can be broadly
# categorized into grocery shopping or clothes shopping. We take the subset
# of the data corresponding to these two (almost balanced) categories, to obtain
# a two-class hypergraph

CF = findall(x->in(x,[1,7]),labels)
labels = labels[CF]
classes = round.(Int64,(labels .- 1) ./6)
n = length(CF)

H2 = H[:,CF]
order = vec(sum(H2,dims = 2))
nonempty = findall(x->x>0,order)
H = H2[nonempty,:]
# class 1 is food

# Turn labels to two classes
alpha, B1, B2, H1, H2, R1, R2, N = Hypergraph_to_Scores(H,classes,25)
hr, hd = clique_expansion_homophily(N,r)
@show alpha

## Plot
lw = 2.5
ms = 9
gfs = 15
tfs = 15
titlesize = 14
s1 = 350
s2 = 300
k = 14
color1 = :blue
color2 = :gray
title = "$k-size hyperedge, $namelabel"
title = ""
leg = false
xlab = "Affinity type t"
ylab = "Affinity"
p = Plots.plot(fontfamily = "helvetica",linewidth = 2,legend = false,grid = false, size = (s1,s2),color = :gray)
title!(p,title,titlefontsize = titlesize)
xaxis!(p,xlab, xticks = k/5:k/5:k,tickfontsize=tfs,guidefontsize = gfs)
yaxis!(p,ylab,tickfontsize=tfs,guidefontsize = gfs)

plot!(p,1:k,H2[k,1:k],color = color2, linewidth = lw, markersize = ms,
    markershape = :circle, markerstrokecolor = color2, label = class2)
plot!(p,1:k,H1[k,1:k],color = color1, linewidth = lw, markersize = ms,
    markershape = :circle,markerstrokecolor = color1, label = class1)

savefig("Figures/ClothesFood/Affinities_$(k)_Label_$class1.pdf")

## Next
for k = 4:4:12
s1 = 400
s2 = 300
title = "$k-size hyperedge"
title = ""
leg = false
xlab = "Affinity type t"
ylab = "Affinity / Baseline"
p = Plots.plot([1; k],[1; 1],yscale = :log10,fontfamily = "helvetica",linewidth = 2,legend = false,grid = false, size = (s1,s2),color = :gray)
title!(p,title,titlefontsize = titlesize)
xaxis!(p,xlab, xticks = k/5:k/5:k,tickfontsize=tfs,guidefontsize = gfs)
yaxis!(p,ylab,tickfontsize=tfs,guidefontsize = gfs)
plot!(p,1:k,R1[k,1:k],color = color1, linewidth = lw, markersize = ms,
    markershape = :circle,markerstrokecolor = color1, label = class1)
plot!(p,1:k,R2[k,1:k],color = color2, linewidth = lw, markersize = ms,
    markershape = :circle, markerstrokecolor = color2, label = class2)

savefig("Figures/ClothesFood/Ratios_$(k)_Label_$class1.pdf")
end

## Group Homophily Index, print usable LaTeX code for building table
include("../src/hypergraph-affinity-functions.jl")
S1,S2  = GHI(H,classes,14)

max_stars_repo_name: nveldt/HypergraphHomophily

max_issues_repo_name: nveldt/HypergraphHomophily

max_forks_repo_name: nveldt/HypergraphHomophily


text:
# Flags
# =====
#
# Bitwise flags (or FLAG).

"the read is paired in sequencing, no matter whether it is mapped in a pair"
const SAM_FLAG_PAIRED        = UInt16(0x1)

"the read is mapped in a proper pair"
const SAM_FLAG_PROPER_PAIR   = UInt16(0x2)

"the read itself is unmapped; conflictive with SAM_FLAG_PROPER_PAIR"
const SAM_FLAG_UNMAP         = UInt16(0x4)

"the mate is unmapped"
const SAM_FLAG_MUNMAP        = UInt16(0x8)

"the read is mapped to the reverse strand"
const SAM_FLAG_REVERSE       = UInt16(0x10)

"the mate is mapped to the reverse strand"
const SAM_FLAG_MREVERSE      = UInt16(0x20)

"this is read1"
const SAM_FLAG_READ1         = UInt16(0x40)

"this is read2"
const SAM_FLAG_READ2         = UInt16(0x80)

"not primary alignment"
const SAM_FLAG_SECONDARY     = UInt16(0x100)

"QC failure"
const SAM_FLAG_QCFAIL        = UInt16(0x200)

"optical or PCR duplicate"
const SAM_FLAG_DUP           = UInt16(0x400)

"supplementary alignment"
const SAM_FLAG_SUPPLEMENTARY = UInt16(0x800)


text:
using Revise
using AdFem
using PyPlot
using LinearAlgebra
using Statistics
using MAT
function f(x, y)
        1/(1+x^2) + x * y + y^2
    end

mmesh = Mesh(50, 50, 1/50, degree=2)

left = bcnode((x,y)->x<1e-5, mmesh)
right = bcedge((x1,y1,x2,y2)->(x1>0.049-1e-5) && (x2>0.049-1e-5), mmesh)

t1 = eval_f_on_boundary_edge((x,y)->1.0e-4, right, mmesh)
t2 = eval_f_on_boundary_edge((x,y)->0.0, right, mmesh)
rhs = compute_fem_traction_term(t1, t2, right, mmesh)

nu = 0.3
x = gauss_nodes(mmesh)
E = Variable(ones(get_ngauss(mmesh)))
# E = constant(eval_f_on_gauss_pts(f, mmesh))

D = compute_plane_stress_matrix(E, nu*ones(get_ngauss(mmesh)))
K = compute_fem_stiffness_matrix(D, mmesh)

bdval = [eval_f_on_boundary_node((x,y)->0.0, left, mmesh);
        eval_f_on_boundary_node((x,y)->0.0, left, mmesh)]
DOF = [left;left .+ mmesh.ndof]
K, rhs = impose_Dirichlet_boundary_conditions(K, rhs, DOF, bdval)
u = K\rhs
U = matread("fenics/data2.mat")["u"]

using Random; Random.seed!(233)
# idx = rand(1:mmesh.ndof, 1000)
idx = 1:mmesh.ndof
idx = [idx; idx .+ mmesh.ndof]

loss = sum((U[idx] - u[idx])^2) * 1e10
sess = Session(); init(sess)

_loss = Float64[]
cb = (vs, iter, loss)->begin
    global _loss
    push!(_loss, loss)
    printstyled("[#iter $iter] loss=$loss\n", color=:green)
    nu = vs[1]
    close("all")
    visualize_scalar_on_gauss_points(nu, mmesh)
    if mod(iter, 10)==1
        matwrite("fenics/data3-inverse$iter.mat", Dict("iter"=>iter,"loss"=>_loss, "E"=>nu))
        savefig("fenics/inverse3_pixel$iter.png")
    end
end

# run(sess, loss)
loss_ = BFGS!(sess, loss, vars = [E], callback = cb)


text:
module Module1

export splitrange, gapdict_from_bitarray, primegaps_from_range, combine_primegaps, primegapsdict

using Primes
using StatsBase
# using ThreadsX
using Distributed
using Plots

# greet() = print("Hello World!")

# splits up range into smaller chunks, avoids crashes
function splitrange(N, maxN=10^8)
  i=3
  RangeArray = []
  while i+maxN  <= N
    RangeArray = vcat(RangeArray, (i, i+maxN))
    i += maxN
  end
  RangeArray = vcat(RangeArray, (i, N))
end

# goes through a bitarray and records all gaps between trues
function gapdict_from_bitarray(Primesbitarray)
  i = 0
  i1 = i
  gapsarray = []
  while true
    i1 = findnext(Primesbitarray, i+1)
    try
      append!(gapsarray, [(i1 - i)])
    catch
      break
    end
    i = i1
  end
  prevgap = popfirst!(gapsarray)
  endgap = length(Primesbitarray) -i
  return prevgap, countmap(gapsarray), endgap
end

function primegaps_from_range(Range)
  gapdict_from_bitarray(primesmask(Range[1], Range[2]))
end

function combine_primegaps(GapsStructArray)
  EndGaps = Dict(2=>0)
  for (i,GapsStruct) in pairs(IndexLinear(), GapsStructArray[1:end-1])
    addcounts!(EndGaps, [(GapsStruct[end] +GapsStructArray[i+1][1] -1)])
  end

  IntGaps = Dict(2=>0)
  for GapsStruct in GapsStructArray
    IntGaps = merge(+, IntGaps, GapsStruct[2])
  end

  return sort(delete!(merge(+, IntGaps, EndGaps), 0))
end

function primegapsdict(N)
  # GapsStructArray = map(primegaps_from_range, splitrange(N))
  # GapsStructArray = ThreadsX.map(primegaps_from_range, splitrange(N))
  GapsStructArray = pmap(primegaps_from_range, splitrange(N))
  combine_primegaps(GapsStructArray)
end

end # module

max_stars_repo_name: chortbauer/primegaps

max_issues_repo_name: chortbauer/primegaps

max_forks_repo_name: chortbauer/primegaps


text:
export BrickletIndustrialDual020mAV2CurrentCallbackConfiguration
struct BrickletIndustrialDual020mAV2CurrentCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletIndustrialDual020mAV2ChannelLEDStatusConfig
struct BrickletIndustrialDual020mAV2ChannelLEDStatusConfig
    min::Integer
    max::Integer
    config::Integer
end

export BrickletIndustrialDual020mAV2SPITFPErrorCount
struct BrickletIndustrialDual020mAV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletIndustrialDual020mAV2Identity
struct BrickletIndustrialDual020mAV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletIndustrialDual020mAV2
"""
Measures two DC currents between 0mA and 20mA (IEC 60381-1)
"""
mutable struct BrickletIndustrialDual020mAV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletIndustrialDual020mAV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_industrial_dual_0_20ma_v2")
        deviceInternal = package.BrickletIndustrialDual020mAV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_current
"""
    $(SIGNATURES)

Returns the current of the specified channel.

It is possible to detect if an IEC 60381-1 compatible sensor is connected
and if it works probably.

If the returned current is below 4mA, there is likely no sensor connected
or the connected sensor is defective. If the returned current is over 20mA,
there might be a short circuit or the sensor is defective.

If you want to get the value periodically, it is recommended to use the
:cb:`Current` callback. You can set the callback configuration
with :func:`Set Current Callback Configuration`.
"""
function get_current(device::BrickletIndustrialDual020mAV2, channel)
    return device.deviceInternal.get_current(channel)
end

export set_current_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Current` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Current` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_current_callback_configuration(device::BrickletIndustrialDual020mAV2, channel, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_current_callback_configuration(channel, period, value_has_to_change, option, min, max)
end

export get_current_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Current Callback Configuration`.
"""
function get_current_callback_configuration(device::BrickletIndustrialDual020mAV2, channel)
    return device.deviceInternal.get_current_callback_configuration(channel)
end

export set_sample_rate
"""
    $(SIGNATURES)

Sets the sample rate to either 240, 60, 15 or 4 samples per second.
The resolution for the rates is 12, 14, 16 and 18 bit respectively.

.. csv-table::
 :header: "Value", "Description"
 :widths: 10, 100

 "0",    "240 samples per second, 12 bit resolution"
 "1",    "60 samples per second, 14 bit resolution"
 "2",    "15 samples per second, 16 bit resolution"
 "3",    "4 samples per second, 18 bit resolution"
"""
function set_sample_rate(device::BrickletIndustrialDual020mAV2, rate)
    device.deviceInternal.set_sample_rate(rate)
end

export get_sample_rate
"""
    $(SIGNATURES)

Returns the sample rate as set by :func:`Set Sample Rate`.
"""
function get_sample_rate(device::BrickletIndustrialDual020mAV2)
    return device.deviceInternal.get_sample_rate()
end

export set_gain
"""
    $(SIGNATURES)

Sets a gain between 1x and 8x. If you want to measure a very small current,
you can increase the gain to get some more resolution.

Example: If you measure 0.5mA with a gain of 8x the return value will be
4mA.
"""
function set_gain(device::BrickletIndustrialDual020mAV2, gain)
    device.deviceInternal.set_gain(gain)
end

export get_gain
"""
    $(SIGNATURES)

Returns the gain as set by :func:`Set Gain`.
"""
function get_gain(device::BrickletIndustrialDual020mAV2)
    return device.deviceInternal.get_gain()
end

export set_channel_led_config
"""
    $(SIGNATURES)

Each channel has a corresponding LED. You can turn the LED off, on or show a
heartbeat. You can also set the LED to "Channel Status". In this mode the
LED can either be turned on with a pre-defined threshold or the intensity
of the LED can change with the measured value.

You can configure the channel status behavior with :func:`Set Channel LED Status Config`.
"""
function set_channel_led_config(device::BrickletIndustrialDual020mAV2, channel, config)
    device.deviceInternal.set_channel_led_config(channel, config)
end

export get_channel_led_config
"""
    $(SIGNATURES)

Returns the channel LED configuration as set by :func:`Set Channel LED Config`
"""
function get_channel_led_config(device::BrickletIndustrialDual020mAV2, channel)
    return device.deviceInternal.get_channel_led_config(channel)
end

export set_channel_led_status_config
"""
    $(SIGNATURES)

Sets the channel LED status config. This config is used if the channel LED is
configured as "Channel Status", see :func:`Set Channel LED Config`.

For each channel you can choose between threshold and intensity mode.

In threshold mode you can define a positive or a negative threshold.
For a positive threshold set the "min" parameter to the threshold value in nA
above which the LED should turn on and set the "max" parameter to 0. Example:
If you set a positive threshold of 10mA, the LED will turn on as soon as the
current exceeds 10mA and turn off again if it goes below 10mA.
For a negative threshold set the "max" parameter to the threshold value in nA
below which the LED should turn on and set the "min" parameter to 0. Example:
If you set a negative threshold of 10mA, the LED will turn on as soon as the
current goes below 10mA and the LED will turn off when the current exceeds 10mA.

In intensity mode you can define a range in nA that is used to scale the brightness
of the LED. Example with min=4mA and max=20mA: The LED is off at 4mA and below,
on at 20mA and above and the brightness is linearly scaled between the values
4mA and 20mA. If the min value is greater than the max value, the LED brightness
is scaled the other way around.
"""
function set_channel_led_status_config(device::BrickletIndustrialDual020mAV2, channel, min, max, config)
    device.deviceInternal.set_channel_led_status_config(channel, min, max, config)
end

export get_channel_led_status_config
"""
    $(SIGNATURES)

Returns the channel LED status configuration as set by
:func:`Set Channel LED Status Config`.
"""
function get_channel_led_status_config(device::BrickletIndustrialDual020mAV2, channel)
    return device.deviceInternal.get_channel_led_status_config(channel)
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletIndustrialDual020mAV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletIndustrialDual020mAV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletIndustrialDual020mAV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletIndustrialDual020mAV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletIndustrialDual020mAV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletIndustrialDual020mAV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletIndustrialDual020mAV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletIndustrialDual020mAV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletIndustrialDual020mAV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletIndustrialDual020mAV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletIndustrialDual020mAV2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletIndustrialDual020mAV2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletIndustrialDual020mAV2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using WaterLily
using StaticArrays

function block(L=2^5;Re=250,U=1,amp=0,ϵ=0.5,thk=2ϵ+√2)
    # Set viscosity
    ν=U*L/Re

    # Create dynamic block geometry
    function sdf(x,t)
        y = x .- SVector(0.,clamp(x[2],-L/2,L/2))
        √sum(abs2,y)-thk/2
    end
    function map(x,t)
        α = amp*cos(t*U/L); R = @SMatrix [cos(α) sin(α); -sin(α) cos(α)]
        R * (x.-SVector(3L+L*sin(t*U/L)+0.01,4L))
    end
    body = AutoBody(sdf,map)

    Simulation((6L+2,6L+2),zeros(2),L;U,ν,body,ϵ)
end

using BenchmarkTools
test() = @benchmark sim_step!(sim,π/4,remeasure=true) setup=(sim=block())

# include("TwoD_plots.jl")
# sim_gif!(block();duration=4π,step=π/16,remeasure=true)
# sim_gif!(block(amp=π/4);duration=8π,step=π/16,remeasure=true,μbody=true,cfill=:Blues,legend=false,border=:none)


text:
import ..Tagging: hastag, traithastag, tag, mergetag
import ..Util: mergef
import ..Traits: traits
export AbstractΩ, defΩ, defω, resolve, idtype, replacetags, like

# # Sample Space
# A sample space represents a set of possible values.
# Sample spaces are structured in the sense that that are composed of parts.
# Different parts are indicated by an id (see IDS)

"Abstract Sample Space"
abstract type AbstractΩ end

# # Tags
"""
`tag(ω::AbstractΩ, tags)`

tag `ω` with `tags`.
"""
function tag(ω::AbstractΩ, tags, mergefunc=mergetag)
  replacetags(ω, mergef(mergefunc, ω.tags, tags))
end

rmtag(ω::AbstractΩ, tag) =
  replacetags(ω, rmkey(ω.tags, tag))

updatetag(ω::AbstractΩ, tag, val) =
  replacetags(ω, update(ω.tags, tag, val))

traithastag(t::AbstractΩ, tag) = traithastag(t.tags, tag)
hastag(ω::AbstractΩ, tag) = hastag(ω.tags, tag)

traits(ω::AbstractΩ) = traits(ω.tags) # FIXME: Do this at the type level

# # Defaults
"Default sample space"
function defΩ end

"Default sample space object"
function defω end

# FIXME MOve this somewhere (shouldnt really be in AbstractΩ)

function replacetags end

function resolve end

"`ids(ω)` Collection of ids in `ω`, i.e. domain of ω"
function ids end

function idtype end

function like end


text:
# define the special PCNone type for use with direct solvers

"""
  Preconditioner type for direct solve.
  Do not use with PetscLinearOperator.

  **Public Fields**

   * none
"""
mutable struct PCNone <: AbstractPC
  is_shared::Bool
end

"""
  Outer constructor

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
"""
function PCNone(mesh::AbstractMesh, sbp::AbstractOperator,
                eqn::AbstractSolutionData, opts::Dict)

  is_shared = false
  return PCNone(is_shared)
end

function free(pc::PCNone)

  return nothing
end

function calcPC(pc::PCNone, mesh::AbstractMesh, sbp::AbstractOperator,
                eqn::AbstractSolutionData, opts::Dict, ctx_residual, t)

  return nothing
end

function applyPC(pc::PCNone, mesh::AbstractMesh, sbp::AbstractOperator,
                 eqn::AbstractSolutionData, opts::Dict, t, b::AbstractVector,
                 x::AbstractVector)

  return nothing
end

function applyPCTranspose(pc::PCNone, mesh::AbstractMesh, sbp::AbstractOperator,
                 eqn::AbstractSolutionData, opts::Dict, t, b::AbstractVector,
                 x::AbstractVector)

  return nothing
end

function getBasePC(pc::PCNone)

  # this is the bottom of the recursion tree
  return pc
end


text:
"""
**`Helium.csv2he`** -*Function*.

    Helium.csv2he(source, outputFile, DataType; hasColNames::Bool=true,
                   hasRowNames::Bool=false, strMiss::String="na",
                   sep::String =",", skipCol::Int64=0) => file

Convert a CSV file that contains a matrix into the Helium binary file format.
The function first evaluate the size of the matrix, then read the data and save
it in Helium format. If no output path is provided then a Helium file will be
created with the same CSV file name.
If there are not column names, set `hasColNames = false`. If
there are row names, set `hasRowNames = true`.

"""
function csv2he(csvFile::String, heFile::String, matType::DataType;
                hasColNames::Bool=true, hasRowNames::Bool=false,
                strMiss::String="na", sep::String=",", skipCol::Int64=0)

    # HeInfo constructor
    he = HeInfo(abspath(csvFile), matType, 1, 1,
                hasColNames, hasRowNames,
                strMiss, sep, skipCol)

    # Convert CSV to matrix
    mat, cNames, rNames, supp = csv2mat(he)

    writehe(mat, heFile, colNames = cNames, rowNames = rNames, supplement = supp)

end

"""
**`Helium.csv2he`** -*Function*.

    `Helium.csv2he(source, DataType; hasColNames::Bool=true,
                   hasRowNames::Bool=false, strMiss::String = "na",
                   sep::String =",", skipCol::Int64=0)` => `file`

Convert a CSV file that contains a matrix into the Helium binary file format.
The function first evaluate the size of the matrix, then read the data and save
it in Helium format. If no output path is provided then a Helium file will be
created with the same CSV file name.
If there are not column names, set `hasColNames = false`. If
there are row names, set `hasRowNames = true`.

"""
function csv2he(csvFile::String, matType::DataType;
                hasColNames::Bool=true, hasRowNames::Bool=false,
                strMiss::String="na", sep::String=",", skipCol::Int64 =0)

    heFile = string(csvFile[1:end-3], "he")

    csv2he(csvFile, heFile, matType,
                    hasColNames = hasColNames, hasRowNames = hasRowNames,
                    strMiss = strMiss, sep = sep, skipCol = skipCol)

    println("Location of the Helium file:")
    println(realpath(heFile))
end

"""
**`Helium.csv2mat`** -*Function*.

    `Helium.csv2mat(HeAttributes)` => `matrix, column names, row names`

Extract data matrix, column and row names from a CSV file.
The function first evaluate the size of the matrix, then read the data.
If there are not column names, set `hasColNames = false`. If
there are row names, set `hasRowNames = true`.

"""
function csv2mat(he::HeAttributes)

    # Get the number of columns
    frstRow = readline(he.fileName)
    he.numCols = length(split(frstRow, he.sep))

    # Get the number of rows
    global nRows = 0
    open(he.fileName) do io
        while !eof(io)
            readline(io)
            nRows += 1
        end
    end
    he.numRows = nRows

    # Get matrix
    mat, colNames, rowNames, matXtra = Dict(true=>getmat,
                                            false=>getskipmat)[he.skipCol == 0](he)
    return mat, colNames, rowNames, matXtra

end


text:
@macromethod d(x -> x) esc(:($x -> 1))
@macromethod d(x -> y) esc(:($x -> 0))

@macromethod d(x -> x^(:T{n, Number})) esc(:($x -> $n*($x^($n-1))))

@macromethod d(x -> a+b) esc(:($x -> @d($x -> $a)($x) + @d($x -> $b)($x)))
@macromethod d(x -> a-b) esc(:($x -> @d($x -> $a)($x) - @d($x -> $b)))

@macromethod d(x -> a*b) esc(:($x -> $a * @d($x -> $b)($x) +
                                     $b * @d($x -> $a)($x)))

@macromethod d(x -> a/b) esc(:($x -> ($b * @d($x -> $a)($x) -
                                      $a * @d($x -> $b)($x))/
                                     ($b^2)))

@d(x -> 4x)(1) == 4

@d(x -> 10x^2 + 5x)(10) == (x-> 20x + 5)(10)

@d(x -> 20x^2/(x+2))(12) == (x-> (40x*(x+2) - (20x^2))/((x+2)^2))(12)


text:
using Spotify
using Spotify.Library
if !Spotify.has_ig_access_token()
    Spotify.get_implicit_grant()
    Spotify.wait_for_ig_access()
end
offset = 0
st, repeat_after_sec = library_get_saved_tracks(1, offset, "NO")
@info "Listing $(st.total) tracks. Pauses may occur due the the current API rate limit."

repeat_after_sec = 0
while true
    while true
        repeat_after_sec !=0 && sleep(repeat_after_sec + 0.5)
        st, repeat_after_sec = library_get_saved_tracks(50, offset, "NO")
        repeat_after_sec == 0 && break
    end
    for (i, t) in enumerate(st.items)
        println(lpad(i + offset, 5), t.track.name)
    end
    offset += length(st.items)
    offset >= st.total && break
end


text:
import TPT
using Base.Test

using DataFrames
using Optim
using Plots; pyplot()

println("--- TM Optim ---")

# Load elemental parameters
p = readtable(joinpath("data", "parameters", "tm_optim.csv"))
N = size(p, 1) # number of elements

#
# Prepare a directory to store the results
#
!isdir("results") && mkdir("results")

resdir = joinpath("results", "tm_optim")
!isdir(resdir) && mkdir(resdir)

#
# Load experimental g(r) from csv files
#
g_exp = Vector{Vector{Tuple{Float64,Float64}}}(N)
for i in 1:N
    elem = p[:X][i]
    T = convert(Int, p[:T][i])

    data = readtable(joinpath("data", "rdf", "$(elem)$(T).csv"))

    M = size(data, 1)

    g_exp[i] = Vector{Tuple{Float64,Float64}}(M)
    for j in 1:M
        g_exp[i][j] = (data[:r][j] / 0.5291, data[:g][j])
    end
end

p[:T] = convert(DataArrays.DataArray{Float64,1}, p[:T])

#
# Prepare TPTSystem for each system
#
ahs = [ TPT.AHS(σ = p[:σ][i], ρ = p[:ρ][i]) for i in 1:N ]
wca = [ TPT.LWCA(ahs[i], p[:T][i], struct=:full) for i in 1:N ]
tb = [ TPT.HHTB(p[:zd][i], p[:Ed][i], p[:Wd][i], p[:R0][i]) for i in 1:N ]

a₀ = p[:a]
a = zeros(N)
res = zeros(N)
sys = Vector{TPT.TPTSystem}(N)
F = Vector{Vector{Float64}}(N)

#
# Performe optimization of rc
#
Threads.@threads for i in 1:N
  function fopt(a::Float64)::Float64
    pse = TPT.BretonnetSilbert(p[:zs][i], p[:rc][i], a)
    nfe = TPT.NFE(wca[i], pse)
    nfetb = TPT.NFETB(nfe, tb[i])
    sys[i] = TPT.TPTSystem(wca[i], nfetb)

    g = TPT.paircorrelation(sys[i])[1,1]

    residue::Float64 = 0

    for j in 1:length(g_exp[i])
      residue += abs(g_exp[i][j][2] - g(g_exp[i][j][1]))
    end

    return residue
  end

  # fnlopt(x::Vector{Float64}, g)::Float64 = fopt(x[1])
  #
  # opt = NLopt.Opt(:GN_DIRECT, 1)
  # NLopt.min_objective!(opt, fnlopt)
  # NLopt.stopval!(opt, 1e-0)
  # NLopt.xtol_abs!(opt, 1e-3)
  # NLopt.ftol_abs!(opt, 1e-4)
  # NLopt.lower_bounds!(opt, [0.5rc₀[i]])
  # NLopt.upper_bounds!(opt, [1.5rc₀[i]])
  #
  # (fmin, xmin, result) = NLopt.optimize(opt, [rc₀[i]])
  #
  # res[i] = result
  # rc[i] = xmin[1]
  # fopt(xmin)

  result = Optim.optimize(fopt, 0.25, 0.35, rel_tol = 1e-3)
  a[i] = Optim.minimizer(result)
  res[i] = Optim.minimum(result)

  fopt(a[i])

  #
  # Density dependency of free-energy
  #
  # F[i] = Vector{Float64}(5)
  # for k in 1:5
  #   ρ = (0.8 + (k-1)*0.1) * p[:ρ][i]
  #
  #   ahs = TPT.AHS(σ = p[:σ][i], ρ = ρ)
  #   wca = TPT.LWCA(ahs, p[:T][i], struct=:full)
  #   tb = TPT.WHTB(p[:zd][i], p[:rd][i], version=:original)
  #   pse = TPT.BretonnetSilbert(p[:zs][i], p[:rc][i], a[i])
  #   nfe = TPT.NFE(wca, pse)
  #   nfetb = TPT.NFETB(nfe, tb)
  #   sys = TPT.TPTSystem(wca, nfetb, m = p[:m][i])
  #
  #   F[i][k] = TPT.helmholtz(sys)
  # end
end

#
# Print the numerical results to console
#
println("The optimized values of a:")
for i in 1:N
  @printf "  %2s: %1.3f\n" p[:X][i] a[i]
end

rmin = zeros(N)
println("The positions of minimum of pairpotentials:")
for i in 1:N
  rmin[i] = sys[i].ref.rmin[1,1]
  @printf "  %2s: %1.3f\n" p[:X][i] rmin[i]
end

σ_hs = zeros(N)
println("The effective hard-sphere diameters:")
for i in 1:N
  σ_hs[i] = sys[i].ref.trial.σ[1]
  @printf "  %2s: %1.3f\n" p[:X][i] σ_hs[i]
end

#
# Save the numerical results as a csv file
#
df = DataFrame(X = p[:X], a = round(a, 3), res = round(res, 1), rmin = round(rmin, 3), σ_hs = round(σ_hs, 3))
writetable(joinpath(resdir, "results.csv"), df)

# Overwrite the parameters for the subsequent tests
p[:a] = round(a, 3)
writetable(joinpath("data", "parameters", "tm_optim.csv"), p)

#
# Save the graphical results as png files
#

# 1. Pairpotentials
for i in 1:N
  u_nfe = TPT.pairpotential(sys[i].pert.nfe)[1,1]
  u_tb = TPT.pairpotential(sys[i].pert.tb)[1,1]
  u_tot = TPT.pairpotential(sys[i].pert)[1,1]

  plot([u_nfe, u_tb, u_tot], 2, 20, ylims=(-0.05, 0.05), labels = ["NFE" "TB" "total"], xlabel="r (a.u.)", ylabel="u(r) (a.u.)")
  vline!([σ_hs[i]], label="HS")
  png(joinpath(resdir, "$(i)-$(p[:X][i])_u"))
end

# 2. Pair correlation functions
for i in 1:N
  g_hs = TPT.paircorrelation(ahs[i])[1,1]
  g_wca = TPT.paircorrelation(sys[i])[1,1]

  plot([g_hs, g_wca], 2, 20, labels=["HS" "WCA"], xlabel="r (a.u.)", ylabel="g(r)", ylims=:auto)
  plot!(g_exp[i], label="exp", xlims=(2,20), ylims=:auto)
  png(joinpath(resdir, "$(i)-$(p[:X][i])_g"))
end

# 3. Density dependency of free-energy
# for i in 1:N
#   plot(0.8:0.1:1.2, 2625.5*F[i], xlabel="Relative density", ylabel="F (kJ/mol)", ylims=:auto)
#   png(joinpath(resdir, "$(i)-$(p[:X][i])_F"))
# end

#
# Check the results
#
@testset "TM Optim" begin
  # for i in 1:N
  #   @test isapprox(a[i], a₀[i], atol=1e-3)
  # end
end # testset


text:
export Diag

# Diag can have either Vector storage, in which case
# it is a general Diag tensor, or scalar storage,
# in which case it has a uniform value
mutable struct Diag{T} <: TensorStorage
  data::T
  Diag{T}(data::Vector) where {T<:AbstractVector} = new{T}(data)
  Diag{T}(size::Integer) where {T<:AbstractVector{S}} where S = new{T}(Vector{S}(undef,size))
  Diag{T}(x::Number,size::Integer) where {T<:AbstractVector} = new{T}(fill(convert(eltype(T),x),size))
  #Diag{T}() where {T} = new{T}(Vector{T}())
  # Make a uniform Diag storage
  Diag{T}(x::Number) where {T<:Number} = new{T}(convert(T,x))
  # Determine the storage type parameter from the input
  Diag(data::T) where T = new{T}(data)
end

data(D::Diag) = D.data
eltype(D::Diag) = eltype(data(D))
getindex(D::Diag{T},i::Int) where {T<:AbstractVector}= data(D)[i]
# Version of getindex for uniform (scalar) storage
getindex(D::Diag{T},i::Int) where {T<:Number} = data(D)
*(D::Diag{T},x::S) where {T<:AbstractVector,S<:Number} = Dense{promote_type(eltype(D),S)}(x*data(D))
*(x::Number,D::Diag) = D*x

#
# Type promotions involving Diag
# Useful for knowing how conversions should work when adding and contracting
#

Base.promote_type(::Type{Diag{T1}},::Type{Diag{T2}}) where {T1<:AbstractVector,T2<:AbstractVector} = Diag{promote_type(T1,T2)}

Base.promote_type(::Type{Diag{T1}},::Type{Diag{T2}}) where {T1<:Number,T2<:Number} = Diag{promote_type(T1,T2)}

Base.promote_type(::Type{Diag{T1}},::Type{Diag{T2}}) where {T1<:Number,T2<:Vector{S2}} where S2 = Diag{Vector{promote_type(T1,S2)}}

Base.promote_type(::Type{Diag{T1}},::Type{Diag{T2}}) where {T1<:Vector{S1},T2<:Number} where S1 = promote_type(Diag{T2},Diag{T1})

Base.promote_type(::Type{Dense{T1}},::Type{Diag{T2}}) where {T1<:Number,T2<:Vector{S2}} where S2 = Dense{promote_type(T1,S2)}

Base.promote_type(::Type{Dense{T1}},::Type{Diag{T2}}) where {T1<:Number,T2<:Number} = Dense{promote_type(T1,T2)}

Base.promote_type(::Type{Diag{T1}},::Type{Dense{T2}}) where {T1,T2} = promote_type(Dense{T2},Diag{T1})

# TODO: define length for Diag, should this be
# the length of the diagonal
# or the product of the tensor dimensions?
# How is the length defined for scalar storage?
length(D::Diag{T}) where {T<:AbstractVector} = length(data(D))

# convert to Dense
function storage_dense(D::Diag,is::IndexSet)
  return Dense{eltype(D)}(vec(storage_convert(Array,D,is)))
end

# convert to complex
storage_complex(D::Diag) = Diag(complex(data(D)))

copy(D::Diag{T}) where {T} = Diag{T}(copy(data(D)))

# TODO: implement this in a sparse way
# For now, we will just make them dense since the output is dense anyway
function storage_outer(D1::Diag{T1},is1::IndexSet,D2::Diag{T2},is2::IndexSet) where {T1,T2}
  A1 = storage_dense(D1,is1)
  A2 = storage_dense(D2,is2)
  return storage_outer(A1,is1,A2,is2)
end

# Convert to an array by setting the diagonal elements
function storage_convert(::Type{Array},D::Diag,is::IndexSet)
  A = zeros(eltype(D),dims(is))
  for i = 1:minDim(is)
    A[fill(i,length(is))...] = D[i]
  end
  return A
end

Base.convert(::Type{Diag{T}},D::Diag) where T = Diag{T}(data(D))

function storage_convert(::Type{Dense{T}},D::Diag,is::IndexSet) where T
  return Dense{T}(vec(storage_convert(Array,D,is)))
end

storage_fill!(D::Diag,x::Number) = fill!(data(D),x)

function diag_getindex(Tstore::Diag{<:AbstractVector},
                       val::Int)
  return getindex(data(Tstore),val)
end

# Uniform case
function diag_getindex(Tstore::Diag{<:Number},
                       val::Int)
  return data(Tstore)
end

# Get diagonal elements
# Gives zero for off-diagonal elements
function storage_getindex(Tstore::Diag{T},
                          Tis::IndexSet,
                          vals::Union{Int, AbstractVector{Int}}...) where {T}
  if all(==(vals[1]),vals)
    return diag_getindex(Tstore,vals[1])
  else
    return zero(eltype(T))
  end
end

# Set diagonal elements
# Throw error for off-diagonal
function storage_setindex!(Tstore::Diag{<:AbstractVector},
                           Tis::IndexSet,
                           x::Union{<:Number, AbstractArray{<:Number}},
                           vals::Union{Int, AbstractVector{Int}}...)
  all(y->y==vals[1],vals) || error("Cannot set off-diagonal element of Diag storage")
  return setindex!(data(Tstore),x,vals[1])
end

function storage_setindex!(Tstore::Diag{<:Number},
                           Tis::IndexSet,
                           x::Union{<:Number, AbstractArray{<:Number}},
                           vals::Union{Int, AbstractVector{Int}}...)
  error("Cannot set elements of a uniform Diag storage")
end

# Add generic Diag's in-place
function _add!(Bstore::Diag,
               Bis::IndexSet,
               Astore::Diag,
               Ais::IndexSet,
               x::Number = 1)
  # This is just used to check if the index sets
  # are permutations of each other, maybe
  # this should be at the ITensor level
  p = calculate_permutation(Bis,Ais)
  Adata = data(Astore)
  Bdata = data(Bstore)
  if x == 1
    Bdata .+= Adata
  else
    Bdata .+= x .* Adata
  end
end

# Uniform Diag case
function _add!(Bstore::Diag{BT},
               Bis::IndexSet,
               Astore::Diag{AT},
               Ais::IndexSet,
               x::Number = 1) where {BT<:Number,AT<:Number}
  # This is just used to check if the index sets
  # are permutations of each other, maybe
  # this should be at the ITensor level
  p = calculate_permutation(Bis,Ais)
  if x == 1
    Bstore.data += data(Astore)
  else
    Bstore.data += x * data(Astore)
  end
end

function _add!(Bstore::Dense,
               Bis::IndexSet,
               Astore::Diag,
               Ais::IndexSet,
               x::Number = 1)
  # This is just used to check if the index sets
  # are permutations of each other
  p = calculate_permutation(Bis,Ais)
  Bdata = data(Bstore)
  mindim = minDim(Bis)
  reshapeBdata = reshape(Bdata,dims(Bis))
  if x == 1
    for ii = 1:mindim
      # TODO: this should be optimized, maybe use
      # strides instead of reshape?
      reshapeBdata[fill(ii,order(Bis))...] += Astore[ii]
    end
  else
    for ii = 1:mindim
      reshapeBdata[fill(ii,order(Bis))...] += x*Astore[ii]
    end
  end
end

# TODO: implement this
# This should not require any permutation
function storage_copyto!(Bstore::Diag,
                         Bis::IndexSet,
                         Astore::Diag,
                         Ais::IndexSet,
                         x::Number = 1)
  p = calculate_permutation(Bis,Ais)
  Adata = data(Astore)
  Bdata = data(Bstore)
  if is_trivial_permutation(p)
    if x == 1
      Bdata .= Adata
    else
      Bdata .= x .* Adata
    end
  else
    reshapeBdata = reshape(Bdata,dims(Bis))
    reshapeAdata = reshape(Adata,dims(Ais))
    if x == 1
      _add!(reshapeBdata,reshapeAdata,p,(a,b)->b)
    else
      _add!(reshapeBdata,reshapeAdata,p,(a,b)->x*b)
    end
  end
end

# TODO: Move to tensorstorage.jl?
function storage_mult!(Astore::Diag,
                       x::Number)
  Adata = data(Astore)
  rmul!(Adata, x)
end

# TODO: Move to tensorstorage.jl?
function storage_mult(Astore::Diag{T},
                      x::Number) where T
  Bdata = x*data(Astore)
  return Diag(Bdata)
end

# TODO: make this a special version of storage_add!()
# This shouldn't do anything to the data
function storage_permute!(Bstore::Diag,
                          Bis::IndexSet,
                          Astore::Diag,
                          Ais::IndexSet)
  p = calculate_permutation(Bis,Ais)
  Adata = data(Astore)
  Bdata = data(Bstore)
  if is_trivial_permutation(p)
    Bdata .= Adata
  else
    reshapeBdata = reshape(Bdata,dims(Bis))
    permutedims!(reshapeBdata,reshape(Adata,dims(Ais)),p)
  end
end

# TODO: move this to tensorstorage.jl?
function storage_dag(Astore::Diag,Ais::IndexSet)
  return dag(Ais),storage_conj(Astore)
end

# TODO: move this to tensorstorage.jl?
function storage_scalar(D::Diag{T}) where {T<:AbstractVector}
  length(D)==1 && return D[1]
  throw(ErrorException("Cannot convert Diag -> Number for length of data greater than 1"))
end

function storage_scalar(D::Diag{T}) where {T<:Number}
  return data(T)
end

# Contract function for Diag*Dense
function _contract(Cinds::IndexSet,
                   Clabels::Vector{Int},
                   Astore::Diag{Vector{SA}},
                   Ainds::IndexSet,
                   Alabels::Vector{Int},
                   Bstore::Dense{SB},
                   Binds::IndexSet,
                   Blabels::Vector{Int}) where {SA<:Number,SB<:Number}
  SC = promote_type(SA,SB)

  # Convert the arrays to a common type
  # since we will call BLAS
  Astore = convert(Diag{Vector{SC}},Astore)
  Bstore = convert(Dense{SC},Bstore)

  Adims = dims(Ainds)
  Bdims = dims(Binds)
  Cdims = dims(Cinds)

  # Create storage for output tensor
  # This needs to be filled with zeros, since in general
  # the output of Diag*Dense will be sparse
  Cstore = Dense{SC}(zero(SC),prod(Cdims))

  # This seems unnecessary, and makes this function
  # non-generic. I think Dense should already store
  # the properly-shaped data (and be parametrized by the
  # order, and maybe store the size?)
  Adata = data(Astore)
  Bdata = reshape(data(Bstore),Bdims)
  Cdata = reshape(data(Cstore),Cdims)

  # Functions to do the contraction
  if(length(Alabels)==0)
    _contract_scalar!(Cdata,Clabels,Bdata,Blabels,Adata[1])
  elseif(length(Blabels)==0)
    _contract_scalar!(Cdata,Clabels,Adata,Alabels,Bdata[1])
  else
    _contract_diag_dense!(Cdata,Clabels,Adata,Alabels,Bdata,Blabels)
  end

  return Cstore
end

# Contract function for Diag uniform * Dense
function _contract(Cinds::IndexSet,
                   Clabels::Vector{Int},
                   Astore::Diag{SA},
                   Ainds::IndexSet,
                   Alabels::Vector{Int},
                   Bstore::Dense{SB},
                   Binds::IndexSet,
                   Blabels::Vector{Int}) where {SA<:Number,SB<:Number}
  SC = promote_type(SA,SB)

  # Convert the arrays to a common type
  # since we will call BLAS
  #Astore = convert(Diag{SC},Astore)
  #Bstore = convert(Dense{SC},Bstore)

  Bdims = dims(Binds)
  Cdims = dims(Cinds)

  # Create storage for output tensor
  # This needs to be filled with zeros, since in general
  # the output of Diag*Dense will be sparse
  Cstore = Dense{SC}(zero(SC),prod(Cdims))

  # This seems unnecessary, and makes this function
  # non-generic. I think Dense should already store
  # the properly-shaped data (and be parametrized by the
  # order, and maybe store the size?)
  Adata = data(Astore)
  Bdata = reshape(data(Bstore),Bdims)
  Cdata = reshape(data(Cstore),Cdims)

  # Functions to do the contraction
  if(length(Alabels)==0)
    _contract_scalar!(Cdata,Clabels,Bdata,Blabels,Adata[1])
  elseif(length(Blabels)==0)
    _contract_scalar!(Cdata,Clabels,Adata,Alabels,Bdata[1])
  elseif(length(Alabels)==2 && length(Clabels)==length(Blabels))
    # This is just a replacement of the indices
    # TODO: This logic should be higher up
    for i = 1:length(Clabels)
      if Blabels[i] > 0
        Cinds[i] = Binds[i]
      else
        if Alabels[1] > 0
          Cinds[i] = Ainds[1]
        else
          Cinds[i] = Ainds[2]
        end
      end
    end
    Cdata .= Adata .* Bdata
  else
    _contract_diag_dense!(Cstore,Cinds,Clabels,
                          Astore,Ainds,Alabels,
                          Bstore,Binds,Blabels)
  end

  return Cstore
end

# Contract function for Dense*Diag (calls Diag*Dense)
function _contract(Cinds::IndexSet,
                   Clabels::Vector{Int},
                   Astore::Dense,
                   Ainds::IndexSet,
                   Alabels::Vector{Int},
                   Bstore::Diag,
                   Binds::IndexSet,
                   Blabels::Vector{Int})
  return _contract(Cinds,Clabels,
                   Bstore,Binds,Blabels,
                   Astore,Ainds,Alabels)
end

# Contract function for Diag*Diag
function _contract(Cinds::IndexSet,
                   Clabels::Vector{Int},
                   Astore::Diag{Vector{SA}},
                   Ainds::IndexSet,
                   Alabels::Vector{Int},
                   Bstore::Diag{Vector{SB}},
                   Binds::IndexSet,
                   Blabels::Vector{Int}) where {SA<:Number,SB<:Number}
  SC = promote_type(SA,SB)

  # Convert the arrays to a common type
  # since we will call BLAS
  Astore = convert(Diag{Vector{SC}},Astore)
  Bstore = convert(Diag{Vector{SC}},Bstore)
  Cstore = Diag{Vector{SC}}(minDim(Cinds))

  Adata = data(Astore)
  Bdata = data(Bstore)
  Cdata = data(Cstore)

  if(length(Alabels)==0)
    _contract_scalar!(Cdata,Clabels,Bdata,Blabels,Adata[1])
  elseif(length(Blabels)==0)
    _contract_scalar!(Cdata,Clabels,Adata,Alabels,Bdata[1])
  else
    _contract_diag_diag!(Cdata,Clabels,Adata,Alabels,Bdata,Blabels)
  end

  return Cstore
end

# Contract function for Diag uniform * Diag uniform
function _contract(Cinds::IndexSet,
                   Clabels::Vector{Int},
                   Astore::Diag{SA},
                   Ainds::IndexSet,
                   Alabels::Vector{Int},
                   Bstore::Diag{SB},
                   Binds::IndexSet,
                   Blabels::Vector{Int}) where {SA<:Number,SB<:Number}
  SC = promote_type(SA,SB)

  # Convert the arrays to a common type
  # since we will call BLAS
  Astore = convert(Diag{SC},Astore)
  Bstore = convert(Diag{SC},Bstore)
  Cstore = Diag{SC}(minDim(Cinds))

  Adata = data(Astore)
  Bdata = data(Bstore)
  Cdata = data(Cstore)

  if length(Clabels) == 0  # If all indices of A and B are contracted
    # all indices are summed over, just add the product of the diagonal
    # elements of A and B
    min_dim = minimum(dims(Ainds)) # == length(Bdata)
    # Need to set to zero since
    # Cdata is originally uninitialized memory
    # (so causes random output)
    Cdata = zero(SC)
    for i = 1:min_dim
      Cdata += Adata*Bdata
    end
  else
    # not all indices are summed over, set the diagonals of the result
    # to the product of the diagonals of A and B
    Cdata = Adata*Bdata
  end

  Cstore.data = Cdata
  return Cstore
end

function _contract_diag_dense!(Cdata::Array{T,NC},Clabels::Vector{Int},
                               Adata::Vector{T},Alabels::Vector{Int},
                               Bdata::Array{T,NB},Blabels::Vector{Int}) where {T,NB,NC}
  if all(i -> i < 0, Blabels)  # If all of B is contracted
    min_dim = minimum(size(Bdata))
    if length(Clabels) == 0
      # all indices are summed over, just add the product of the diagonal
      # elements of A and B
      for i = 1:min_dim
        Cdata[1] += Adata[i]*Bdata[ntuple(_->i,Val(NB))...]
      end
    else
      # not all indices are summed over, set the diagonals of the result
      # to the product of the diagonals of A and B
      # TODO: should we make this return a Diag storage?
      for i = 1:min_dim
        Cdata[ntuple(_->i,Val(NC))...] = Adata[i]*Bdata[ntuple(_->i,Val(NB))...]
      end
    end
  else
    astarts = zeros(Int,length(Alabels))
    bstart = 0
    cstart = 0
    b_cstride = 0
    nbu = 0
    for ib = 1:length(Blabels)
      ia = findfirst(==(Blabels[ib]),Alabels)
      if(!isnothing(ia))
        b_cstride += stride(Bdata,ib)
        bstart += astarts[ia]*stride(Bdata,ib)
      else
        nbu += 1
      end
    end

    c_cstride = 0
    for ic = 1:length(Clabels)
      ia = findfirst(==(Clabels[ic]),Alabels)
      if(!isnothing(ia))
        c_cstride += stride(Cdata,ic)
        cstart += astarts[ia]*stride(Cdata,ic)
      end
    end

    # strides of the uncontracted dimensions of
    # Bdata
    bustride = zeros(Int,nbu)
    custride = zeros(Int,nbu)
    # size of the uncontracted dimensions of
    # Bdata, to be used in CartesianIndices
    busize = zeros(Int,nbu)
    n = 1
    for ib = 1:length(Blabels)
      if Blabels[ib] > 0
        bustride[n] = stride(Bdata,ib)
        busize[n] = size(Bdata,ib)
        ic = findfirst(==(Blabels[ib]),Clabels)
        custride[n] = stride(Cdata,ic)
        n += 1
      end
    end

    boffset_orig = 1-sum(strides(Bdata))
    coffset_orig = 1-sum(strides(Cdata))
    cartesian_inds = CartesianIndices(Tuple(busize))
    for inds in cartesian_inds
      boffset = boffset_orig
      coffset = coffset_orig
      for i in 1:nbu
        ii = inds[i]
        boffset += ii*bustride[i]
        coffset += ii*custride[i]
      end
      for j in 1:length(Adata)
        Cdata[cstart+j*c_cstride+coffset] += Adata[j]*Bdata[bstart+j*b_cstride+boffset]
      end
    end
  end
end

function _contract_diag_dense!(Cstore::Dense,Cinds,Clabels::Vector{Int},
                               Astore::Diag{TA},Ainds,Alabels::Vector{Int},
                               Bstore::Dense,Binds,Blabels::Vector{Int}) where {TA<:Number}
  if all(i -> i < 0, Blabels)  # If all of B is contracted
    Bdims = dims(Binds)
    min_dim = minimum(Bdims)
    rB = reshape(data(Bstore),Bdims)
    if length(Clabels) == 0
      # all indices are summed over, just add the product of the diagonal
      # elements of A and B
      # TODO: replace this with manual strides
      for i = 1:min_dim
        Cstore[1] += Astore[i]*rB[ntuple(_->i,length(Binds))...]
      end
    else
      # not all indices are summed over, set the diagonals of the result
      # to the product of the diagonals of A and B
      # TODO: should we make this return a Diag storage?
      for i = 1:min_dim
        Cstore[ntuple(_->i,length(Cinds))...] = Astore[i]*rB[ntuple(_->i,length(Binds))...]
      end
    end
  else
    astarts = zeros(Int,length(Alabels))
    bstart = 0
    cstart = 0
    b_cstride = 0
    nbu = 0
    for ib = 1:length(Blabels)
      ia = findfirst(==(Blabels[ib]),Alabels)
      if(!isnothing(ia))
        b_cstride += stride(Binds,ib)
        bstart += astarts[ia]*stride(Binds,ib)
      else
        nbu += 1
      end
    end

    c_cstride = 0
    for ic = 1:length(Clabels)
      ia = findfirst(==(Clabels[ic]),Alabels)
      if(!isnothing(ia))
        c_cstride += stride(Cinds,ic)
        cstart += astarts[ia]*stride(Cinds,ic)
      end
    end

    # strides of the uncontracted dimensions of
    # Bdata
    bustride = zeros(Int,nbu)
    custride = zeros(Int,nbu)
    # size of the uncontracted dimensions of
    # Bdata, to be used in CartesianIndices
    busize = zeros(Int,nbu)
    custride = zeros(Int,nbu)
    # size of the uncontracted dimensions of
    # Bdata, to be used in CartesianIndices
    busize = zeros(Int,nbu)
    n = 1
    for ib = 1:length(Blabels)
      if Blabels[ib] > 0
        bustride[n] = stride(Binds,ib)
        busize[n] = dim(Binds,ib) #size(Bstore,ib)
        ic = findfirst(==(Blabels[ib]),Clabels)
        custride[n] = stride(Cinds,ic)
        n += 1
      end
    end

    min_dim = minimum(dims(Binds))
    boffset_orig = 1-sum(strides(Binds))
    coffset_orig = 1-sum(strides(Cinds))
    cartesian_inds = CartesianIndices(Tuple(busize))
    for inds in cartesian_inds
      boffset = boffset_orig
      coffset = coffset_orig
      for i in 1:nbu
        ii = inds[i]
        boffset += ii*bustride[i]
        coffset += ii*custride[i]
      end
      for j in 1:min_dim
        Cstore[cstart+j*c_cstride+coffset] += Astore[j]*Bstore[bstart+j*b_cstride+boffset]
      end
    end
  end
end

# Maybe this works for uniform storage as well
function _contract_diag_diag!(Cdata::Vector{T},Clabels::Vector{Int},
                              Adata::Vector{T},Alabels::Vector{Int},
                              Bdata::Vector{T},Blabels::Vector{Int}) where T
  if length(Clabels) == 0  # If all indices of A and B are contracted
    # all indices are summed over, just add the product of the diagonal
    # elements of A and B
    Adim = length(Adata)  # == length(Bdata)
    # Need to set to zero since
    # Cdata is originally uninitialized memory
    # (so causes random output)
    Cdata[1] = zero(T)
    for i = 1:Adim
      Cdata[1] += Adata[i]*Bdata[i]
    end
  else
    min_dim = min(length(Adata),length(Bdata))
    # not all indices are summed over, set the diagonals of the result
    # to the product of the diagonals of A and B
    for i = 1:min_dim
      Cdata[i] = Adata[i]*Bdata[i]
    end
  end
end


text:
using Optim
using Parameters
using LinearAlgebra
using GeometryTypes

### If C=0 we have original Zinchenko 1997 stabilization scheme
@with_kw struct KineticStabilizer
    C::AbstractFloat = 0.4
    ftol::AbstractFloat = 1e-6
end

function F(v,points,faces,zc::KineticStabilizer)

    Cp = zc.C

    s = 0
    for ti in 1:length(faces)
        v1,v2,v3 = faces[ti]

        a = norm(points[v2] - points[v1])
        b = norm(points[v3] - points[v2])
        c = norm(points[v1] - points[v3])

        xava = dot(points[v2]-points[v1],v[v2]-v[v1])
        xbvb = dot(points[v3]-points[v2],v[v3]-v[v2])
        xcvc = dot(points[v1]-points[v3],v[v1]-v[v3])

        ### This part is responsible for first part in the sum
        s += v2>v1 ? 4 * xava^2 : 0
        s += v3>v2 ? 4 * xbvb^2 : 0
        s += v1>v3 ? 4 * xcvc^2 : 0

        ### This part is needed for ordinary things
        Cdelta = 1/4*sqrt(1 - 2* (a^4 + b^4 + c^4)/(a^2 + b^2 + c^2)^2)
        A = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( a^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        B = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( b^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        C = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( c^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        DCdelta = -A*xava - B*xbvb - C*xcvc

        s += Cp*DCdelta^2/Cdelta^2
    end

    return s
end

function gradF!(v,points,faces,storage,zc::KineticStabilizer)

    Cp = zc.C

    for i in 1:length(storage)
        storage[i] = Point(0,0,0)
    end

    for i in 1:length(faces)
        v1,v2,v3 = faces[i]
        a = norm(points[v2] - points[v1])
        b = norm(points[v3] - points[v2])
        c = norm(points[v1] - points[v3])

        Cdelta = 1/4*sqrt(1 - 2 * (a^4 + b^4 + c^4)/(a^2 + b^2 + c^2)^2)
        A = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( a^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        B = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( b^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        C = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( c^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)

        xa = points[v2] - points[v1]
        xb = points[v3] - points[v2]
        xc = points[v1] - points[v3]

        xava = dot(points[v2]-points[v1],v[v2]-v[v1])
        xbvb = dot(points[v3]-points[v2],v[v3]-v[v2])
        xcvc = dot(points[v1]-points[v3],v[v1]-v[v3])

        DCdelta = -A*xava - B*xbvb - C*xcvc

        A_ = 2*Cp*DCdelta*A/Cdelta^2
        B_ = 2*Cp*DCdelta*B/Cdelta^2
        C_ = 2*Cp*DCdelta*C/Cdelta^2

        storage[v1] += (A_ - 8*xava)*xa - C_*xc
        storage[v2] += (B_ - 8*xbvb)*xb - A_*xa
        storage[v3] += (C_ - 8*xcvc)*xc - B_*xb
    end
end

function stabilize!(v,points,faces,n,zc::KineticStabilizer)
    #h2 = hvec(points,faces,n)

    ### Let's try to make it simpler

    function f(x::Vector)
        # vv = reshape(x,3,div(length(x),3))
        return F(x,points,faces,zc)
    end

    function g!(storage::Vector,x::Vector)

        # vv = reshape(x,size(points)...)
        # gradf = reshape(storage,size(points)...)
        gradF!(x,points,faces,storage,zc)
        for i in 1:size(v,2)
            P = I - n[i]*n[i]'
            storage[i] = P*storage[i]
        end
    end

    res = optimize(f,g!,v,method=Optim.ConjugateGradient(),ftol=zc.ftol)

    ### One needs to be a little bit carefull about what types would Julia return
    return res
    #v[:,:] = reshape(res.minimum, size(v)...)[:,:]
end