text:
struct SVRG_basic_iterable{R<:Real,C<:RealOrComplex{R},Tx<:AbstractArray{C},Tf,Tg}
    F::Array{Tf}            # smooth term
    g::Tg                   # nonsmooth term
    x0::Tx                  # initial point
    N::Int                  # of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}}  # Lipschitz moduli of nabla f_i
    μ::Maybe{Union{Array{R},R}}  # convexity moduli of the gradients
    γ::Maybe{R}             # stepsize
    m::Maybe{Int}           # number of inner loop updates
    plus::Bool              # for SVRG++ variant
end

mutable struct SVRG_basic_state{R<:Real,Tx}
    γ::R                    # stepsize
    m::Int                  # number of inner loop updates
    av::Tx                  # the running average
    z::Tx
    z_full::Tx              # the outer loop argument
    w::Tx                   # the inner loop variable
    ind::Array{Int}         # running idx set
    # some extra placeholders
    ∇f_temp::Tx             # placeholder for gradients
    temp::Tx
end

function SVRG_basic_state(γ::R, m, av::Tx, z::Tx, z_full::Tx, w::Tx, ind) where {R,Tx}
    return SVRG_basic_state{R,Tx}(γ, m, av, z, z_full, w, ind, copy(av), copy(av))
end

function Base.iterate(iter::SVRG_basic_iterable{R}) where {R}
    N = iter.N
    ind = collect(1:N)
    m = iter.m === nothing ? m = 2 * N : m = iter.m
    # updating the stepsize
    if iter.γ === nothing
        if iter.plus
            @warn "provide a stepsize γ"
            return nothing
        else
            if iter.L === nothing || iter.μ === nothing
                @warn "smoothness or convexity parameter absent"
                return nothing
            else
                L_M = maximum(iter.L)
                μ_M = maximum(iter.μ)
                γ = 1 / (10 * L_M)
                # condition Theorem 3.1
                rho = (1 + 4 * L_M * γ^2 * μ_M * (N + 1)) / (μ_M * γ * N * (1 - 4L_M * γ))
                if rho >= 1
                    @warn "convergence condition violated...provide a stepsize!"
                end
            end
        end
    else
        γ = iter.γ # provided γ
    end
    # initializing the vectors
    av = zero(iter.x0)
    for i = 1:N
        ∇f, ~ = gradient(iter.F[i], iter.x0)
        ∇f ./= N
        av .+= ∇f
    end
    z_full = copy(iter.x0)
    z = zero(av)
    w = copy(iter.x0)
    state = SVRG_basic_state(γ, m, av, z, z_full, w, ind)
    print()
    return state, state
end

function Base.iterate(iter::SVRG_basic_iterable{R}, state::SVRG_basic_state{R}) where {R}
    # The inner cycle
    for i in rand(state.ind, state.m)
        gradient!(state.temp, iter.F[i], state.z_full)
        gradient!(state.∇f_temp, iter.F[i], state.w)
        state.temp .-= state.∇f_temp
        state.temp .-= state.av
        state.temp .*= state.γ
        state.temp .+= state.w
        CIAOAlgorithms.prox!(state.w, iter.g, state.temp, state.γ)
        state.z .+= state.w   # keeping track of the sum of w's
    end
    # full update
    state.z_full .= state.w #state.z ./ state.m
    # iter.plus || (state.w .= state.z_full) # only for basic SVRG
    state.z = zero(state.z)  # for next iterate
    state.av .= state.z
    for i = 1:iter.N
        gradient!(state.∇f_temp, iter.F[i], state.z_full)
        state.∇f_temp ./= iter.N
        state.av .+= state.∇f_temp
    end
    # iter.plus && (state.m *= 2) # only for SVRG++

    return state, state
end

solution(state::SVRG_basic_state) = state.z_full


text:
a = "code"
b = "wa.rs"
name = a*b

ext: jl

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using GeometricBase.Config
using GeometricBase.Utils
using GeometricIntegrators.Integrators
using GeometricIntegrators.Integrators.VPRK
using GeometricIntegrators.Tableaus
using GeometricProblems.LotkaVolterra2d
using SimpleSolvers
using Test

SimpleSolvers.set_config(:nls_atol, 8eps())
SimpleSolvers.set_config(:nls_rtol, 2eps())

const Δt = 0.01
const nt = 10
const q₀ = [1.0, 1.0]
const parameters = (a₁=1.0, a₂=1.0, b₁=-1.0, b₂=-2.0)

ode  = lotka_volterra_2d_ode(q₀; params=parameters)
iode = lotka_volterra_2d_iode(q₀; params=parameters)
lode = lotka_volterra_2d_lode(q₀; params=parameters)
ldae = lotka_volterra_2d_ldae(q₀; params=parameters)

int  = IntegratorFIRK(ode, TableauGauss(8), Δt)
sol  = integrate(ode, int, nt)
refx = sol.q[end]

@testset "$(rpad("VPRK integrators",80))" begin

    sol = integrate(iode, TableauVPRK(:pglrk, 2, TableauGauss(1), -1), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-6

    sol = integrate(iode, TableauVPRK(:pglrk, 4, TableauGauss(2), +1), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 8E-7

    sol = integrate(iode, TableauVPRK(:pglrk, 6, TableauGauss(3), -1), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-12

    sol = integrate(iode, TableauVPLobattoIIIA(2), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    sol = integrate(iode, TableauVPLobattoIIIA(3), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 8E-7

    sol = integrate(iode, TableauVPLobattoIIIA(4), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 3E-11

    sol = integrate(iode, TableauVPLobattoIIIB(2), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-6

    sol = integrate(iode, TableauVPLobattoIIIB(3), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 8E-7

    sol = integrate(iode, TableauVPLobattoIIIB(4), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-11

end

@testset "$(rpad("VPRK integrators with standard projection",80))" begin

    int = IntegratorVPRKpStandard(iode, TableauVPGLRK(1), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = IntegratorVPRKpStandard(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = IntegratorVPRKpStandard(iode, TableauVPGLRK(3), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

end

@testset "$(rpad("VPRK integrators with symplectic projection",80))" begin

    int = IntegratorVPRKpSymplectic(iode, TableauVPGLRK(1), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-6

    int = IntegratorVPRKpSymplectic(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = IntegratorVPRKpSymplectic(iode, TableauVPGLRK(3), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 8E-12

end

@testset "$(rpad("VPRK integrators with symmetric projection",80))" begin

    int = IntegratorVPRKpSymmetric(iode, TableauVPGLRK(1), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = IntegratorVPRKpSymmetric(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = IntegratorVPRKpSymmetric(iode, TableauVPGLRK(3), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-16

end

@testset "$(rpad("VPRK integrators with midpoint projection",80))" begin

    int = IntegratorVPRKpMidpoint(iode, TableauVPGLRK(1), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = IntegratorVPRKpMidpoint(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = IntegratorVPRKpMidpoint(iode, TableauVPGLRK(3), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-16

end

@testset "$(rpad("VPRK integrators with internal projection",80))" begin

    int = IntegratorVPRKpInternal(iode, TableauVPGLRK(1), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-6

    int = IntegratorVPRKpInternal(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = IntegratorVPRKpInternal(iode, TableauVPGLRK(3), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-12

    int = IntegratorVPRKpInternal(iode, TableauVPGLRK(4), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 8E-16

end

@testset "$(rpad("VPRK integrators with projection on secondary constraint",80))" begin

    # TODO: reactivate

    # int = IntegratorVPRKpSecondary(ldae, TableauVPGLRK(1), Δt)
    # sol = integrate(ldae, int, nt)
    # @test relative_maximum_error(sol.q, refx) < 2E-6

    # int = IntegratorVPRKpSecondary(ldae, TableauVPGLRK(2), Δt)
    # sol = integrate(ldae, int, nt)
    # @test relative_maximum_error(sol.q, refx) < 8E-7

    # int = IntegratorVPRKpSecondary(ldae, TableauVPGLRK(3), Δt)
    # sol = integrate(ldae, int, nt)
    # @test relative_maximum_error(sol.q, refx) < 4E-12

end

@testset "$(rpad("VPRK integrators with variational projection",80))" begin

    intV1 = IntegratorVPRKpVariational(iode, TableauVPGLRK(1), Δt)
    solV1 = integrate(iode, intV1, nt)
    @test relative_maximum_error(solV1.q, refx) < 8E-7

    intV2 = IntegratorVPRKpVariational(iode, TableauVPGLRK(2), Δt)
    solV2 = integrate(iode, intV2, nt)
    @test relative_maximum_error(solV2.q, refx) < 8E-8

    intV3 = IntegratorVPRKpVariational(iode, TableauVPGLRK(3), Δt)
    solV3 = integrate(iode, intV3, nt)
    @test relative_maximum_error(solV3.q, refx) < 1E-11

    intQ1 = IntegratorVPRKpVariationalQ(iode, TableauVPGLRK(1), Δt)
    solQ1 = integrate(iode, intQ1, nt)
    @test relative_maximum_error(solQ1.q, refx) < 4E-5

    intQ2 = IntegratorVPRKpVariationalQ(iode, TableauVPGLRK(2), Δt)
    solQ2 = integrate(iode, intQ2, nt)
    @test relative_maximum_error(solQ2.q, refx) < 2E-4

    intQ3 = IntegratorVPRKpVariationalQ(iode, TableauVPGLRK(3), Δt)
    solQ3 = integrate(iode, intQ3, nt)
    @test relative_maximum_error(solQ3.q, refx) < 1E-8

    intP1 = IntegratorVPRKpVariationalP(iode, TableauVPGLRK(1), Δt)
    solP1 = integrate(iode, intP1, nt)
    @test relative_maximum_error(solP1.q, refx) < 8E-7

    intP2 = IntegratorVPRKpVariationalP(iode, TableauVPGLRK(2), Δt)
    solP2 = integrate(iode, intP2, nt)
    @test relative_maximum_error(solP2.q, refx) < 8E-8

    intP3 = IntegratorVPRKpVariationalP(iode, TableauVPGLRK(3), Δt)
    solP3 = integrate(iode, intP3, nt)
    @test relative_maximum_error(solP3.q, refx) < 1E-11

    @test relative_maximum_error(solV1.q, solP1.q[end]) == 0
    @test relative_maximum_error(solV2.q, solP2.q[end]) == 0
    @test relative_maximum_error(solV3.q, solP3.q[end]) == 0

end

@testset "$(rpad("VPRK integrators with projection on Runge-Kutta tableau",80))" begin

    int = IntegratorVPRKpTableau(iode, CoefficientsPGLRK(5), Δt*5)
    sol = integrate(iode, int, div(nt,5))
    @test relative_maximum_error(sol.q, refx) < 8E-12

    int = IntegratorVPRKpTableau(iode, CoefficientsPGLRK(6), Δt*5)
    sol = integrate(iode, int, div(nt,5))
    @test relative_maximum_error(sol.q, refx) < 4E-14

end

@testset "$(rpad("Degenerate symplectic partitioned Runge-Kutta methods",80))" begin

    int = IntegratorVPRKdegenerate(iode, TableauVPGLRK(1), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-5

    int = IntegratorVPRKdegenerate(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 4E-7

    int = IntegratorVPRKdegenerate(iode, TableauVPGLRK(3), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 2E-10

end

@testset "$(rpad("VSPRK integrators with Legendre projection",80))" begin

    int = IntegratorVPRKpLegendre(iode, TableauVPGLRK(1), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-6

    int = IntegratorVPRKpLegendre(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = IntegratorVPRKpLegendre(iode, TableauVPGLRK(3), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 8E-16

end


text:
# Define a function/method
function test1(a::Int64)
    r = a * 53
    return r
end

# Cause specialization of the method
r1 = test1(Int64(1))
println("JIT-compiled version 1: ", r1)

# Collect specialization for compiling
mis = []
m = methods(test1)
mi = m.ms[1].specializations[1]
push!(mis, mi)

# Native compile specialization
native_code = ccall(:jl_simple_create_native, Ptr{Cvoid}, (Vector{Core.MethodInstance},), mis)

# Remember names of native functions
functionObject = mi.cache.functionObject[1:end-1]
specFunctionObject = mi.cache.specFunctionObject[1:end-1]

# Produce a shared library containing the native code
ccall(:jl_simple_dump_native, Cvoid, (Ptr{Cvoid}, Cstring,), native_code, "test.a")
#run(`ar x test.a`)
#run(`ld -shared -fPIC -o test.so data.so text.so`)
run(`ld -shared -fPIC -o test.so --whole-archive test.a`)

# Load the shared library
using Libdl
lib = Libdl.dlopen("./test.so", RTLD_GLOBAL)

# Get pointers to the native functions
jl_fun = Libdl.dlsym(lib, functionObject)
julia_fun = Libdl.dlsym(lib, specFunctionObject)

# Redefine the function
function test1(a::Int64)
    a * 3
end

# Test the redefined function
r2 = test1(Int64(1))
println("JIT-compiled version 2: ", r2)

# Link the native code into the specialization
m = methods(test1)
mi = m.ms[1].specializations[1]
mi.cache.invoke = jl_fun
mi.cache.specptr = julia_fun

# Test the function with its native code
r3 = test1(Int64(1))
println("Native compiled version 1: ", r3)

# Close the library
Libdl.dlclose(lib)
return nothing


text:
is_valid(x, encoding::Symbol = :none) = is_valid(x, Val(encoding))

is_valid(x, ::Val) = true

function binarize(x::Vector{T}, d = length(x); binarization = :one_hot) where {T <: Number}
    return collect(Iterators.flatten(map(i -> binarize(i, d; binarization), x)))
end

function binarize(x::T, domain_size::Int; binarization = :one_hot) where {T <: Number}
    return binarize(x, domain(0:(domain_size - 1)), Val(binarization))
end

function binarize(x::T, d::D; binarization=:one_hot
) where {T <: Number, D <: DiscreteDomain{T}}
    return binarize(x, d, Val(binarization))
end

function binarize(x::Vector{T}, d::Vector{D}; binarization=:one_hot
) where {T <: Number, D <: DiscreteDomain{T}}
    return collect(Iterators.flatten(map(elt -> binarize(elt[1], elt[2]; binarization), zip(x,d))))
end

function debinarize(x; binarization = :one_hot)
    ds::Int = if binarization == :domain_wall
        (-1 + sqrt(1 + 4 * length(x))) ÷ 2 + 1
    else
        sqrt(length(x))
    end
    return debinarize(x, ds; binarization)
end

function debinarize(x, domain_size::Int; binarization = :one_hot)
    return debinarize(x, domain(0:(domain_size - 1)); binarization)
end

function debinarize(x, domains_sizes::Vector{Int}; binarization = :one_hot)
    domains = map(ds -> domain(0:(ds - 1)), domains_sizes)
    return debinarize(x, domains; binarization)
end

function debinarize(x, d::D; binarization = :one_hot) where {D <: DiscreteDomain}
    k::Int = length(x) / length(d)
    if binarization == :domain_wall
        typeof(d) <: RangeDomain && first(get_domain(d)) == 0 && (k += 1)
    end
    domains = fill(d, k)
    return debinarize(x, domains; binarization)
end

function debinarize(x, domains::Vector{D}; binarization = :one_hot
) where {D <: DiscreteDomain}
    return debinarize(x, domains, Val(binarization))
end


text:
include("Include.jl")

# setup trade logic -
function my_technical_agent_trade_logic(agent_model, asset_price_array, iteration_index)

	# initialize -
	order_model_array = Array{VLOrderModel,1}()
	agent_position_array = agent_model.agent_position_array

	# get the parameters from this agent -
	parameters_dict = agent_model.agent_parameters

	# get the c-parameter array -
	c_parameter_array = parameters_dict["c_parameter_array"]
	price_lag_array = parameters_dict["price_lag_array"]

	# compute the setpoint array -
	(number_of_iterations, number_of_assets) = size(asset_price_array)
	for asset_index = 1:number_of_assets

		# what is the lag?
		asset_lag_index = price_lag_array[asset_index]

		# what is the lagged price for this asset?
		price_setpoint_value = 1.0
		if (asset_lag_index < iteration_index)
			price_setpoint_value = asset_price_array[(iteration_index - asset_lag_index),asset_index]
		end

		# build the order model -
		c_value = c_parameter_array[asset_index]
		current_price = asset_price_array[asset_index]
		desired_position_size = c_value * sign(current_price - price_setpoint_value)
		Δω = agent_position_array[iteration_index, asset_index] - desired_position_size

		# build the order -
		# an order has this data -
		local_asset_index = asset_index
    	action = (Δω < 0) ? -1 : 1
    	order_type = 1
    	# quantity = abs(Δω)
        quantity = 2
    	price = current_price
    	agent_id = agent_model.agent_id
		order_model = VLOrderModel(local_asset_index, action, order_type, quantity, price, agent_id);
		push!(order_model_array, order_model)
	end

	# return -
	return order_model_array
end

# build the game world -
number_of_iterations = 20
number_of_agents = 1
number_of_assets = 2
game_world = build(VLGameWorld; number_of_iterations=number_of_iterations, number_of_assets=number_of_assets);

# build agent array -
my_agent_array = Array{VLAgentModel,1}()
for agent_index = 1:number_of_agents

    # let's use the defaults for the various agent settings -
    agent_model = build(VLAgentModel; number_of_iterations=number_of_iterations, number_of_assets=number_of_assets)

    # setup parameters for this agent -
    c_array = rand(1:number_of_assets, number_of_assets)
    lag_array = rand(1:1, number_of_assets)
    p = Dict("c_parameter_array" => c_array, "price_lag_array" => lag_array)
    agent_model.agent_parameters = p

    # let's have all my agents be technical traders -
    agent_model.agent_trade_logic = my_technical_agent_trade_logic

    # grab -
    push!(my_agent_array, agent_model)
end

# attach the agents to the world -
game_world.agent_array = my_agent_array

# run the model -
R = run(game_world)


text:
module BioCookbook

import BioSequences
import BioSequences
import BioAlignments
import Statistics

module msa
import FASTX
include("concatenate.jl")
include("simulation.jl")
export concatenate, simAlignments
end

using .msa

module mapping
import XAM
include("extractReads.jl")
export extractReads

module viz
import XAM
import DataFrames
include("clipanalysis.jl")
end

end

using .mapping

end

max_stars_repo_name: compbiocore/BioJulia-cookbook

max_issues_repo_name: compbiocore/BioJulia-cookbook

max_forks_repo_name: compbiocore/BioJulia-cookbook


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

function rhs!(du, u, t,
              mesh::ParallelP4estMesh{3}, equations,
              initial_condition, boundary_conditions, source_terms,
              dg::DG, cache)
  # Start to receive MPI data
  @trixi_timeit timer() "start MPI receive" start_mpi_receive!(cache.mpi_cache)

  # Prolong solution to MPI interfaces
  @trixi_timeit timer() "prolong2mpiinterfaces" prolong2mpiinterfaces!(
    cache, u, mesh, equations, dg.surface_integral, dg)

  # Start to send MPI data
  @trixi_timeit timer() "start MPI send" start_mpi_send!(
    cache.mpi_cache, mesh, equations, dg, cache)

  # Reset du
  @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache)

  # Calculate volume integral
  @trixi_timeit timer() "volume integral" calc_volume_integral!(
    du, u, mesh,
    have_nonconservative_terms(equations), equations,
    dg.volume_integral, dg, cache)

  # Prolong solution to interfaces
  @trixi_timeit timer() "prolong2interfaces" prolong2interfaces!(
    cache, u, mesh, equations, dg.surface_integral, dg)

  # Calculate interface fluxes
  @trixi_timeit timer() "interface flux" calc_interface_flux!(
    cache.elements.surface_flux_values, mesh,
    have_nonconservative_terms(equations), equations,
    dg.surface_integral, dg, cache)

  # Prolong solution to boundaries
  @trixi_timeit timer() "prolong2boundaries" prolong2boundaries!(
    cache, u, mesh, equations, dg.surface_integral, dg)

  # Calculate boundary fluxes
  @trixi_timeit timer() "boundary flux" calc_boundary_flux!(
    cache, t, boundary_conditions, mesh, equations, dg.surface_integral, dg)

  # Prolong solution to mortars
  @trixi_timeit timer() "prolong2mortars" prolong2mortars!(
    cache, u, mesh, equations, dg.mortar, dg.surface_integral, dg)

  # Calculate mortar fluxes
  @trixi_timeit timer() "mortar flux" calc_mortar_flux!(
    cache.elements.surface_flux_values, mesh,
    have_nonconservative_terms(equations), equations,
    dg.mortar, dg.surface_integral, dg, cache)

  # Finish to receive MPI data
  @trixi_timeit timer() "finish MPI receive" finish_mpi_receive!(
    cache.mpi_cache, mesh, equations, dg, cache)

  # Calculate MPI interface fluxes
  @trixi_timeit timer() "MPI interface flux" calc_mpi_interface_flux!(
    cache.elements.surface_flux_values, mesh,
    have_nonconservative_terms(equations), equations,
    dg.surface_integral, dg, cache)

  # Calculate surface integrals
  @trixi_timeit timer() "surface integral" calc_surface_integral!(
    du, u, mesh, equations, dg.surface_integral, dg, cache)

  # Apply Jacobian from mapping to reference element
  @trixi_timeit timer() "Jacobian" apply_jacobian!(
    du, mesh, equations, dg, cache)

  # Calculate source terms
  @trixi_timeit timer() "source terms" calc_sources!(
    du, u, t, source_terms, equations, dg, cache)

  # Finish to send MPI data
  @trixi_timeit timer() "finish MPI send" finish_mpi_send!(cache.mpi_cache)

  return nothing
end

function prolong2mpiinterfaces!(cache, u,
                                mesh::ParallelP4estMesh{3},
                                equations, surface_integral, dg::DG)
  @unpack mpi_interfaces = cache
  index_range = eachnode(dg)

  @threaded for interface in eachmpiinterface(dg, cache)
    # Copy solution data from the local element using "delayed indexing" with
    # a start value and a step size to get the correct face and orientation.
    # Note that in the current implementation, the interface will be
    # "aligned at the primary element", i.e., the index of the primary side
    # will always run forwards.
    local_side = mpi_interfaces.local_sides[interface]
    local_element = mpi_interfaces.local_element_ids[interface]
    local_indices = mpi_interfaces.node_indices[interface]

    i_element_start, i_element_step_i, i_element_step_j = index_to_start_step_3d(local_indices[1], index_range)
    j_element_start, j_element_step_i, j_element_step_j = index_to_start_step_3d(local_indices[2], index_range)
    k_element_start, k_element_step_i, k_element_step_j = index_to_start_step_3d(local_indices[3], index_range)

    i_element = i_element_start
    j_element = j_element_start
    k_element = k_element_start
    for j in eachnode(dg)
      for i in eachnode(dg)
        for v in eachvariable(equations)
          mpi_interfaces.u[local_side, v, i, j, interface] = u[v, i_element, j_element, k_element, local_element]
        end
        i_element += i_element_step_i
        j_element += j_element_step_i
        k_element += k_element_step_i
      end
      i_element += i_element_step_j
      j_element += j_element_step_j
      k_element += k_element_step_j
    end
  end

  return nothing
end

function calc_mpi_interface_flux!(surface_flux_values,
                                  mesh::ParallelP4estMesh{3},
                                  nonconservative_terms,
                                  equations, surface_integral, dg::DG, cache)
  @unpack local_element_ids, node_indices, local_sides = cache.mpi_interfaces
  @unpack contravariant_vectors = cache.elements
  index_range = eachnode(dg)

  @threaded for interface in eachmpiinterface(dg, cache)
    # Get element and side index information on the local element
    local_element = local_element_ids[interface]
    local_indices = node_indices[interface]
    local_direction = indices2direction(local_indices)
    local_side = local_sides[interface]

    # Create the local i,j,k indexing on the local element used to pull normal direction information
    i_element_start, i_element_step_i, i_element_step_j = index_to_start_step_3d(local_indices[1], index_range)
    j_element_start, j_element_step_i, j_element_step_j = index_to_start_step_3d(local_indices[2], index_range)
    k_element_start, k_element_step_i, k_element_step_j = index_to_start_step_3d(local_indices[3], index_range)

    i_element = i_element_start
    j_element = j_element_start
    k_element = k_element_start

    # Initiate the node indices to be used in the surface for loop,
    # the surface flux storage must be indexed in alignment with the local element indexing
    local_surface_indices = surface_indices(local_indices)
    i_surface_start, i_surface_step_i, i_surface_step_j = index_to_start_step_3d(local_surface_indices[1], index_range)
    j_surface_start, j_surface_step_i, j_surface_step_j = index_to_start_step_3d(local_surface_indices[2], index_range)
    i_surface = i_surface_start
    j_surface = j_surface_start

    for j in eachnode(dg)
      for i in eachnode(dg)
        # Get the normal direction on the local element
        # Contravariant vectors at interfaces in negative coordinate direction
        # are pointing inwards. This is handled by `get_normal_direction`.
        normal_direction = get_normal_direction(local_direction, contravariant_vectors,
                                                i_element, j_element, k_element,
                                                local_element)

        calc_mpi_interface_flux!(surface_flux_values, mesh, nonconservative_terms, equations,
                                 surface_integral, dg, cache,
                                 interface, normal_direction,
                                 i, j, local_side,
                                 i_surface, j_surface, local_direction, local_element)

        # Increment local element indices to pull the normal direction
        i_element += i_element_step_i
        j_element += j_element_step_i
        k_element += k_element_step_i
        # Increment the surface node indices along the local element
        i_surface += i_surface_step_i
        j_surface += j_surface_step_i
      end
      # Increment local element indices to pull the normal direction
      i_element += i_element_step_j
      j_element += j_element_step_j
      k_element += k_element_step_j
      # Increment the surface node indices along the local element
      i_surface += i_surface_step_j
      j_surface += j_surface_step_j
    end
  end

  return nothing
end

# Inlined version of the interface flux computation for conservation laws
@inline function calc_mpi_interface_flux!(surface_flux_values,
                                          mesh::P4estMesh{3},
                                          nonconservative_terms::Val{false}, equations,
                                          surface_integral, dg::DG, cache,
                                          interface_index, normal_direction,
                                          interface_i_node_index, interface_j_node_index, local_side,
                                          surface_i_node_index, surface_j_node_index,
                                          local_direction_index, local_element_index)
  @unpack u = cache.mpi_interfaces
  @unpack surface_flux = surface_integral

  u_ll, u_rr = get_surface_node_vars(u, equations, dg,
                                     interface_i_node_index, interface_j_node_index, interface_index)

  if local_side == 1
    flux_ = surface_flux(u_ll, u_rr, normal_direction, equations)
  else # local_side == 2
    flux_ = -surface_flux(u_ll, u_rr, -normal_direction, equations)
  end

  for v in eachvariable(equations)
    surface_flux_values[v, surface_i_node_index, surface_j_node_index,
                        local_direction_index, local_element_index] = flux_[v]
  end
end

end # muladd


text:
module TestSectorPerformance
using AlphaVantage
using Test

TEST_SLEEP_TIME =  parse(Float64, get(ENV, "TEST_SLEEP_TIME", "15"))
MAX_TESTS = parse(Int64, get(ENV, "MAX_TESTS", "1"))

@testset "Sector Performance" begin
    data = sector_performance()
    @test isa(data, AlphaVantageResponse)
    sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
end

end # module


text:
export elementproductmatmul!

"""
    elementproductmatmul!(d, A, B, c)

Compute in-place the product `(A ∘ B)*c`, where `∘` denotes the element-wise  product of two matrices.
"""
function elementproductmatmul!(d::Array{Float64,1}, A::Array{Float64,2}, B::Array{Float64,2}, c::Array{Float64,1})
    nx, ny = size(A)
    @avx for i = 1:nx
        di = zero(eltype(d))
        for j = 1:ny
            di += (A[i, j] * B[i, j]) * c[j]
        end
        d[i] = di
    end
end


text:
##############################################################################
##
## Combination behaves like [A B C ...] without forming it
##
##############################################################################

struct Combination{T} <: AbstractMatrix{T}
    A::Tuple
    cumlength::Vector{Int}
end

function Combination(A::Union{AbstractVector{T}, AbstractMatrix{T}}...) where {T}
    Combination{T}(A, cumsum([size(x, 2) for x in A]))
end

Base.size(c::Combination) = (size(c.A[1], 1), c.cumlength[end])
Base.size(c::Combination, i::Integer) = size(c)[i]

function Base.view(c::Combination, ::Colon, j)
    index = searchsortedfirst(c.cumlength, j)
    newj = index == 1 ? j : j - c.cumlength[index-1]
    view(c.A[index], :, newj)
end

##############################################################################
##
## Returns base of [A B C ...]
## Important: it must be the case that it returns the in order, that is [A B A] returns [true true false] not [false true true]
##
##
##############################################################################
function basis(@nospecialize(xs::AbstractVector...))
    invXX = invsym!(crossprod(collect(xs)))
    return diag(invXX) .> 0
end

function crossprod(xs::Vector{<:AbstractVector})
    XX = zeros(length(xs), length(xs))
    for i in 1:length(xs)
        for j in 1:i
            XX[i, j] = xs[i]' * xs[j]
        end
    end
    for i in 1:length(xs)
        for j in (i+1):length(xs)
            XX[i, j] = XX[j, i]
        end
    end
    return XX
end

# generalized 2inverse
function invsym!(X::AbstractMatrix)
    # The C value adjusts the check to the relative scale of the variable. The C value is equal to the corrected sum of squares for the variable, unless the corrected sum of squares is 0, in which case C is 1. If you specify the NOINT option but not the ABSORB statement, PROC GLM uses the uncorrected sum of squares instead. The default value of the SINGULAR= option, 107, might be too small, but this value is necessary in order to handle the high-degree polynomials used in the literature to compare regression routin
    tols = max.(diag(X), 1)
    for j in 1:size(X, 1)
        d = X[j,j]
        if abs(d) < tols[j] * sqrt(eps())
            X[j,:] .= 0
            X[:,j] .= 0
        else
            X[j,:] = X[j,:] ./ d
            for i in 1:size(X, 1)
                if i != j
                    X[i,:] .= X[i,:] .- X[i,j] .* X[j,:]
                    X[i,j] = -X[i,j] / d
                end
            end
            X[j,j] = 1 / d
        end
    end
    return X
end

function getcols(X::AbstractMatrix,  basecolX::AbstractVector)
    sum(basecolX) == size(X, 2) ? X : X[:, basecolX]
end


text:
struct DistributedVisualizationData{A<:AbstractPData}
  visdata::A
end

function Base.getproperty(x::DistributedVisualizationData, sym::Symbol)
  if sym == :grid
    map_parts(i->i.grid,x.visdata)
  elseif sym == :filebase
    get_part(x.visdata).filebase
  elseif sym == :celldata
    map_parts(i->i.celldata,x.visdata)
  elseif sym == :nodaldata
    map_parts(i->i.nodaldata,x.visdata)
  else
    getfield(x, sym)
  end
end

function Base.propertynames(x::DistributedVisualizationData, private::Bool=false)
  (fieldnames(typeof(x))...,:grid,:filebase,:celldata,:nodaldata)
end

# Define how each object is visualized

function Visualization.visualization_data(
  model::DistributedDiscreteModel{Dc},
  filebase::AbstractString;
  labels=get_face_labeling(model)) where Dc

  parts = get_part_ids(model.models)
  nparts = length(parts)
  vd = map_parts(
    parts,model.models,model.gids.partition,labels.labels) do part,model,gids,labels

    vd = visualization_data(model,filebase;labels=labels)
    vd_cells = vd[end]
    push!(vd_cells.celldata, "gid" => gids.lid_to_gid)
    push!(vd_cells.celldata, "part" => gids.lid_to_part)
    vd
  end
  r = []
  for i in 0:Dc
    push!(r,DistributedVisualizationData(map_parts(x->x[i+1],vd)))
  end
  r
end

function Visualization.visualization_data(
  trian::DistributedTriangulation,
  filebase::AbstractString;
  order=-1,
  nsubcells=-1,
  celldata=nothing,
  cellfields=nothing)

  trians = trian.trians
  parts = get_part_ids(trians)
  nparts = length(trians)

  cdat = _prepare_cdata(trians,celldata)
  fdat = _prepare_fdata(trians,cellfields)

  vd = map_parts(
    parts,trians,cdat,fdat) do part,trian,celldata,cellfields
    _celldata = Dict{Any,Any}(celldata)
    # we do not use "part" since it is likely to be used by the user
    if haskey(_celldata,"piece")
      @unreachable "piece is a reserved cell data name"
    end
    _celldata["piece"] = fill(part,num_cells(trian))
    vd = visualization_data(
      trian,filebase;
      order=order,nsubcells=nsubcells,
      celldata=_celldata,cellfields=cellfields)
    @assert length(vd) == 1
    vd[1]
  end
  [DistributedVisualizationData(vd)]
end

function _prepare_cdata(trians,a::Nothing)
  map_parts(trians) do t
    Dict()
  end
end

function _prepare_cdata(trians,a)
  if length(a) == 0
    return map_parts(trians) do t
      Dict()
    end
  end
  ks = []
  vs = []
  for (k,v) in a
    push!(ks,k)
    push!(vs,v)
  end
  map_parts(vs...) do vs...
    b = []
    for i in 1:length(vs)
      push!(b,ks[i]=>vs[i])
    end
    b
  end
end

function _prepare_fdata(trians,a::Nothing)
  map_parts(trians) do t
    Dict()
  end
end

function _prepare_fdata(trians,a)
  if length(a) == 0
    return map_parts(trians) do t
      Dict()
    end
  end
  ks = []
  vs = []
  for (k,v) in a
    push!(ks,k)
    push!(vs,v.fields)
  end
  map_parts(vs...) do vs...
    b = []
    for i in 1:length(vs)
      push!(b,ks[i]=>vs[i])
    end
    b
  end
end

# Vtk related

function Visualization.write_vtk_file(
  grid::AbstractPData{<:Grid}, filebase; celldata, nodaldata)
  pvtk = Visualization.create_vtk_file(grid,filebase;celldata=celldata,nodaldata=nodaldata)
  map_parts(vtk_save,pvtk)
end

function Visualization.create_vtk_file(
  grid::AbstractPData{<:Grid}, filebase; celldata, nodaldata)
  parts = get_part_ids(grid)
  nparts = length(parts)
  map_parts(parts,grid,celldata,nodaldata) do part,g,c,n
    Visualization.create_pvtk_file(
      g,filebase;
      part=part,nparts=nparts,
      celldata=c,nodaldata=n)
  end
end


text:
# Test of Dante
using Dante, LinearAlgebra, Test

@testset "Riemann Problem" begin
   param, state_GV = solve("PARAM_test.toml")
   # Obtain the initial states
   Rho, U, P, tEnd = setInitRiemann(param.RiemannProblemType, false)
   # Exact solution
   xe, re, ue, pe, ee, te, Me, se =
      getEulerExactSol(Rho[1], U[1], P[1], Rho[end], U[end], P[end], tEnd, 3,
      param.x[3:end-2])
   ϵ = norm(state_GV[3:end-2,3,3,8] .- pe)
   @show ϵ
   @test ϵ < 0.25
end


text:
@testset verbose=true "Iterators" begin
    ℓₘₐₓ = 20

    @testset "D" begin
        Drange = mapslices(
            ℓm′m -> tuple(ℓm′m...),
            SphericalFunctions.WignerDrange(ℓₘₐₓ);
            dims=[2]
        )[:, 1]
        𝔇 = Diterator(Drange, ℓₘₐₓ)
        for (ℓ, 𝔇ˡ) in enumerate(𝔇)
            ℓ -= 1
            Dˡ = [
                (ℓ, m′, m)
                for m in -ℓ:ℓ, m′ in -ℓ:ℓ
            ]
            @test 𝔇ˡ == Dˡ
        end
        @test Base.IteratorSize(typeof(𝔇)) == Base.HasShape{1}()
        @test Base.IteratorEltype(typeof(𝔇)) == Base.HasEltype()
        collection = collect(𝔇)
        @test length(collection) == length(𝔇)
        @test size(collection, 1) == size(𝔇, 1)
        @test size(collection) == size(𝔇)
        @test eltype(collection) == eltype(𝔇)

        for ℓₘᵢₙ in 0:ℓₘₐₓ
            𝔇 = Diterator(Drange, ℓₘₐₓ, ℓₘᵢₙ)
            for (ℓ, 𝔇ˡ) in enumerate(𝔇)
                ℓ += ℓₘᵢₙ - 1
                Dˡ = [
                    (ℓ, m′, m)
                    for m in -ℓ:ℓ, m′ in -ℓ:ℓ
                ]
                @test 𝔇ˡ == Dˡ
            end
        end
    end

    @testset "d" begin
        drange = mapslices(
            ℓm′m -> tuple(ℓm′m...),
            SphericalFunctions.WignerDrange(ℓₘₐₓ);
            dims=[2]
        )[:, 1]
        𝔡 = diterator(drange, ℓₘₐₓ)
        for (ℓ, 𝔡ˡ) in enumerate(𝔡)
            ℓ -= 1
            dˡ = [
                (ℓ, m′, m)
                for m in -ℓ:ℓ, m′ in -ℓ:ℓ
            ]
            @test 𝔡ˡ == dˡ
        end
        @test Base.IteratorSize(typeof(𝔡)) == Base.HasShape{1}()
        @test Base.IteratorEltype(typeof(𝔡)) == Base.HasEltype()
        collection = collect(𝔡)
        @test length(collection) == length(𝔡)
        @test size(collection, 1) == size(𝔡, 1)
        @test size(collection) == size(𝔡)
        @test eltype(collection) == eltype(𝔡)

        for ℓₘᵢₙ in 0:ℓₘₐₓ
            𝔡 = diterator(drange, ℓₘₐₓ, ℓₘᵢₙ)
            for (ℓ, 𝔡ˡ) in enumerate(𝔡)
                ℓ += ℓₘᵢₙ - 1
                dˡ = [
                    (ℓ, m′, m)
                    for m in -ℓ:ℓ, m′ in -ℓ:ℓ
                ]
                @test 𝔡ˡ == dˡ
            end
        end
    end

    @testset "Y" begin
        Yrange = mapslices(
            ℓm -> tuple(ℓm...),
            SphericalFunctions.Yrange(ℓₘₐₓ);
            dims=[2]
        )[:, 1]
        𝔜 = Yiterator(Yrange, ℓₘₐₓ)
        for (ℓ, 𝔜ˡ) in enumerate(𝔜)
            ℓ -= 1
            Yˡ = [
                (ℓ, m)
                for m in -ℓ:ℓ
            ]
            @test 𝔜ˡ == Yˡ
        end
        @test Base.IteratorSize(typeof(𝔜)) == Base.HasShape{1}()
        @test Base.IteratorEltype(typeof(𝔜)) == Base.HasEltype()
        collection = collect(𝔜)
        @test length(collection) == length(𝔜)
        @test size(collection, 1) == size(𝔜, 1)
        @test size(collection) == size(𝔜)
        @test eltype(collection) == eltype(𝔜)

        for ℓₘᵢₙ in 0:ℓₘₐₓ
            𝔜 = Yiterator(Yrange, ℓₘₐₓ, ℓₘᵢₙ)
            for (ℓ, 𝔜ˡ) in enumerate(𝔜)
                ℓ += ℓₘᵢₙ - 1
                Yˡ = [
                    (ℓ, m)
                    for m in -ℓ:ℓ
                ]
                @test 𝔜ˡ == Yˡ
            end
            iₘᵢₙ = Ysize(ℓₘᵢₙ-1)+1
            @test Yrange[iₘᵢₙ] == (ℓₘᵢₙ, -ℓₘᵢₙ)
            𝔜 = Yiterator(Yrange[iₘᵢₙ:end], ℓₘₐₓ, ℓₘᵢₙ, 1)
            for (ℓ, 𝔜ˡ) in enumerate(𝔜)
                ℓ += ℓₘᵢₙ - 1
                Yˡ = [
                    (ℓ, m)
                    for m in -ℓ:ℓ
                ]
                @test 𝔜ˡ == Yˡ
            end
        end
    end

end


text:
module TestTypedTables

include("preamble.jl")
using BangBang
using BangBang: SingletonVector
using ConstructionBase: constructorof
using Test
using TypedTables

@testset "push!!" begin
    tints = Table(a = [1], b = [2])
    @test push!!(tints, (a = 3.0, b = 4.0))::Table == Table(a = [1.0, 3.0], b = [2.0, 4.0])
    @test tints == Table(a = [1], b = [2])

    tints = Table(a = [1], b = [2])
    @test push!!(tints, (a = 3, b = 4)) === tints
    @test tints == Table(a = [1, 3], b = [2, 4])
end

@testset "append!!(::Table, ::$(constructorof(typeof(src)))) [widen]" for src in Any[
    [(a = 3.0, b = 4.0)],
    SingletonVector((a = 3.0, b = 4.0)),
]
    tints = Table(a = [1], b = [2])
    @test append!!(tints, src)::Table == Table(
        a = [1.0, 3.0],
        b = [2.0, 4.0],
    )
    @test tints == Table(a = [1], b = [2])
end

@testset "append!!(::Table, ::$(constructorof(typeof(src)))) [mutation]" for src in Any[
    # Broken:
    # [(a = 3, b = 4)],
    SingletonVector((a = 3, b = 4)),
]
    tints = Table(a = [1], b = [2])
    @test append!!(tints, src) === tints == Table(
        a = [1, 3],
        b = [2, 4],
    )
end

@testset "append!!" begin
    if versionof(TypedTables) >= v"1.2.4" && VERSION >= v"1.4" # could be older
        tints = Table(a = [1], b = [2])
        @test append!!(tints, [(a = 3, b = 4)]) === tints
        @test tints == Table(a = [1, 3], b = [2, 4])
    end
end

@testset "type inference" begin
    @test (@inferred singletonof(Table, (a=1,))) == Table((a=[1]))
end

end  # module


text:
devices = Dict{Symbol, DeviceModel}(:Generators => DeviceModel(PSY.ThermalStandard, ThermalDispatch),
                                    :Loads =>  DeviceModel(PSY.PowerLoad, StaticPowerLoad))
branches = Dict{Symbol, DeviceModel}(:L => DeviceModel(PSY.Line, StaticLineUnbounded))
services = Dict{Symbol, ServiceModel}()

@testset "Operation set ref models" begin
    model_ref = ModelReference(CopperPlatePowerModel, devices, branches, services);
    op_model = OperationModel(TestOptModel, model_ref, c_sys5)
    set_transmission_ref!(op_model, DCPLLPowerModel)
    @test op_model.model_ref.transmission == DCPLLPowerModel

    model_ref = ModelReference(CopperPlatePowerModel, devices, branches, services);
    op_model = OperationModel(TestOptModel, model_ref, c_sys5)
    new_devices = Dict{Symbol, DeviceModel}(:Generators => DeviceModel(PSY.ThermalStandard, ThermalBasicUnitCommitment),
                                            :Loads =>  DeviceModel(PSY.PowerLoad, StaticPowerLoad))
    set_devices_ref!(op_model, new_devices)
    @test op_model.model_ref.devices[:Generators].formulation == ThermalBasicUnitCommitment
    jump_model = op_model.canonical.JuMPmodel
    @test ((JuMP.VariableRef, MOI.ZeroOne) in JuMP.list_of_constraint_types(jump_model)) == true

    model_ref = ModelReference(DCPPowerModel, devices, branches, services);
    op_model = OperationModel(TestOptModel, model_ref, c_sys5)
    new_branches = Dict{Symbol, DeviceModel}(:L => DeviceModel(PSY.Line, StaticLine))
    set_branches_ref!(op_model, new_branches)
    @test op_model.model_ref.branches[:L].formulation == StaticLine
end

@testset "Operation set models" begin
    model_ref = ModelReference(CopperPlatePowerModel, devices, branches, services);
    op_model = OperationModel(TestOptModel, model_ref, c_sys5)
    set_device_model!(op_model, :Generators, DeviceModel(PSY.ThermalStandard, ThermalBasicUnitCommitment))
    @test op_model.model_ref.devices[:Generators].formulation == ThermalBasicUnitCommitment
    jump_model = op_model.canonical.JuMPmodel
    @test ((JuMP.VariableRef, MOI.ZeroOne) in JuMP.list_of_constraint_types(jump_model)) == true

    model_ref = ModelReference(DCPPowerModel, devices, branches, services);
    op_model = OperationModel(TestOptModel, model_ref, c_sys5)
    set_branch_model!(op_model, :L, DeviceModel(PSY.Line, StaticLine))
    @test op_model.model_ref.branches[:L].formulation == StaticLine
end


text:
#using TensorToolbox, Test, LinearAlgebra

println("\n\n****Testing TTtensor.jl")

N=3
G=CoreCell(undef,N)
G[1]=rand(1,4,3)
G[2]=rand(3,6,4)
G[3]=rand(4,3,1)

println("\n...Test core tensors G of sizes: ")
[println(size(G[n])) for n=1:N]

T=TTtensor(G)

println("\n...Test size of TTtensor T: ", size(T))
println("\n...Test ndims of TTtensor T: ", ndims(T))
println("\n...Test TTsvd of TTtensor T: ", TTrank(T))

X=randTTtensor([6,5,4,3],[4,3,2]);
Y=randTTtensor([6,5,4,3],[3,2,2]);
Z=X+Y
err=norm(full(Z)-(full(X)+full(Y)))
println("\n...Test addition: ")
@test err ≈ 0 atol=1e-10

z=innerprod(X,Y)
err=norm(z-innerprod(full(X),full(Y)))
println("\n...Test dot product: ")
@test err ≈ 0 atol=1e-10

X=rand(5,4,3,2)

println("\n...Test tensor X of size: ", size(X))
println("\n...Testing TTsvd.")

T=TTsvd(X)

println("Results:")
println("Type of output T: ", typeof(T))
@test isa(T,TTtensor)
err=norm(full(T) - X)
println("\n...Testing function full, i.e. contracted product (conprod): norm(full(T)-X) = ", err)
@test err ≈ 0 atol=1e-10

println("\n...Testing TTsvd with requested rank.")

G=CoreCell(undef,4);
I=[6,5,6,5]
R=[4,4,4]
G[1]=reshape(rand(I[1],R[1]),(1,I[1],R[1]))
G[2]=reshape(rand(R[1],I[2],R[2]),(R[1],I[2],R[2]))
G[3]=reshape(rand(R[2],I[3],R[3]),(R[2],I[3],R[3]))
G[4]=reshape(rand(R[3],I[4]),(R[3],I[4],1))
T=TTtensor(G)
X=full(T)
println("Created tensor x of size ",size(X),", with TTrank = ",TTrank(T), ".")
println("T1=TTsvd(X,reqrank=$R)")
T1=TTsvd(X,reqrank=R)
err=norm(X-full(T1))
println("norm(X-full(T1)) = ",err)
@test err ≈ 0 atol=1e-8

println("\n...Testing recompression.")

S=TTsvd(T)
err=norm(full(S) - X)
println("S=TTsvd(T): norm(full(S)-X) = ", err)
@test err ≈ 0 atol=1e-10

println("\n...Testing reorthogonalization - function reorth.")

I=[6,5,4,3]
R=[4,3,2]
N=length(R)
T=randTTtensor(I,R)
println("Test TT-tensor T of size: ", size(T)," with flags T.lorth = ",T.lorth," and T.rorth = ",T.rorth,".")
println("\nPerforming (default) left orthogonalization: Tl=reorth(T).")
Tl=reorth(T)
println("Flags: Tl.lorth = ",Tl.lorth,", Tl.rorth = ",Tl.rorth,".")
@test Tl.lorth == true
println("\nPerforming right orthogonalization: Tr=reorth(T,\"right\").")
Tr=reorth(T,"right")
println("Flags: Tr.lorth = ",Tr.lorth,", Tr.rorth = ",Tr.rorth,".")

Gr=Tr.cores;
for n=1:N-1
    Gn=tenmat(Gr[n+1],1);
    @show norm(Diagonal(ones(size(Gn,1)))- Gn*Gn')
end

@test Tr.rorth == true

println("\n...Testing reorthogonalization with overwriting - function reorth!.")
println("\nPerforming (default) left orthogonalization: reorth!(T)")
reorth!(T)
println("Flags: T.lorth = ",T.lorth,", T.rorth = ",T.rorth,".")
@test T.lorth == true
println("\nPerforming right orthogonalization: reorth!(T,\"right\")")
reorth!(T,"right")
println("Flags: T.lorth = ",T.lorth,", T.rorth = ",T.rorth,".")
@test T.lorth == true && T.rorth==true

println("\n...Testing TTsvd on TT-tensors.")
epsilon=1e-8
println("\nFixed precision problem with ϵ = ",epsilon,".")
G=CoreCell(undef,4);
I=[6,5,6,5]
R=[4,4,4]
r=2
G[1]=reshape(rand(I[1],r)*rand(r,R[1]),(1,I[1],R[1]))
G[2]=reshape(rand(R[1]*I[2],r)*rand(r,R[2]),(R[1],I[2],R[2]))
G[3]=reshape(rand(R[2]*I[3],r)*rand(r,R[3]),(R[2],I[3],R[3]))
G[4]=reshape(rand(R[3],r)*rand(r,I[4]),(R[3],I[4],1))
T=TTtensor(G)
println("Created TT-tensor T of size ",size(T),", with TTrank = ",TTrank(T), ", but its actual TT-rank is (2,2,2)")
println("T1=TTsvd(T,ϵ)")
T1=TTsvd(T)
rnk=TTrank(T1)
println("TTrank(T1) = ",rnk)
@test rnk==(2,2,2)
println("norm(full(T)-full(T1)) = ",norm(full(T)-full(T1)))
@test norm(full(T)-full(T1)) ≈ 0 atol=epsilon

R=[2,2,2]
println("\nFixed rank problem with reqrank = ",R,".")
T=TTsvd(rand(6,5,4,3))
println("Created TT-tensor T of size ",size(T),", with TTrank = ",TTrank(T),".")
println("T1=TTsvd(T,reqrank=$R)")
T1=TTsvd(T,reqrank=R)
rnk=TTrank(T1)
println("TTrank(T1) = ",rnk)
@test rnk==(2,2,2)

println("\n... Testing contraction with vectors - function TTtv.")
X=randTTtensor([6,5,4,3,2],[3,3,3,3])
G=X.cores
N=ndims(X)
Isz=[size(X)...]
J=1
L=N-J+1
u=VectorCell(undef,L)
for i=1:L
    u[i]=randn(Isz[J+i-1])
end
Xf=full(X)
let
sz=copy(Isz)
for i=1:L
    sz[i]=1
    global Xf=reshape(tenmat(Xf,J+i-1)'*u[i],tuple(sz...))
end
end
@test TTtv(G,u) ≈ dropdims(Xf) atol=1e-12
J=2;
L=N-J+1;
v=VectorCell(undef,L)
[v[i]=randn(Isz[J+i-1]) for i=1:L]
@test norm(TTtv(X,v,J) - dropdims(contract(G[1],TTtv(G[2:end],v)))) ≈ 0 atol=1e-12
for J=3:N
    L=N-J+1;
    v=VectorCell(undef,L)
    [v[i]=randn(Isz[J+i-1]) for i=1:L]
    T=TTtv(X,v,J)
    Gtest=CoreCell(undef,J-1)
    [Gtest[n]=copy(G[n]) for n=1:J-2]
    Gtest[J-1]=contract(G[J-1],TTtv(G[J:end],v))
    @test norm(full(T)-contract(Gtest)) ≈ 0 atol=1e-12
end

println("\n... Testing transformation of ktensor to TTtensor - function kten2TT.")
Isz=[5,4,3,2]
R=3
T=randktensor(Isz,R)
println("Test ktensor T of size ",size(T),", and k-rank = $R, transofrming to TTtensor Ttt.")
Ttt=kten2TT(T)
err=norm(full(T)-full(Ttt))
println("norm(full(T)-full(Ttt)) = ", err)
@test err ≈ 0 atol=1e-8

println("\n... Testing contraction of full tensor to TTtensor - function contract.")
N=5
Isz=repeat([5],N,1)
R=repeat([3],N-1,1)
T=randTTtensor(Isz,R)
Xsz=[Isz[1:N-1]...;2]
X=rand(Xsz...)
Tfull=full(T)
T1=contract(Tfull,collect(1:N-1),X,collect(1:N-1),[2,1])
Y=contract(T,X,1,4)
@test norm(T1-full(Y)) ≈ 0 atol=1e-12

Xsz=[Isz[1:N-2]...;2]
X=rand(Xsz...)
T2=contract(Tfull,collect(2:N-1),X,collect(1:N-2),[1,3,2])
Y=contract(T,X,2,3)
@test norm(T2-full(Y)) ≈ 0 atol=1e-12

println("\n... Testing element-wise product of two TTtensor - function ewprod.")
N=4;
Isz=[8,7,6,5];
Rx=[7,6,5];
Ry=[6,5,4];
X=randTTtensor(Isz,Rx);
Y=randTTtensor(Isz,Ry);
Zfull=full(X).*full(Y);
Z=ewprod(X,Y)
@test norm(Zfull-full(Z)) ≈ 0 atol=1e-8


text:
using ConcreteAbstractions
using Test

@base struct AbstractFoo{T}
    a
    b::Int
    c::T
    d::Vector{T}
end

@extend struct Foo <: AbstractFoo
    e::T
end

@testset begin
    @test Foo <: AbstractFoo
    foo = Foo(10,10,10,[10],5)
    @test isa(foo, Foo{Int})
    @test fieldnames(typeof(foo)) == (:a,:b,:c,:d,:e)
    @test isa(foo.a, Int)
    @test isa(foo.b, Int)
    @test isa(foo.c, Int)
    @test isa(foo.d, Vector{Int})
    @test isa(foo.e, Int)
end


text:
#Multi-threaded:
#JULIA_NUM_THREADS=<number-of-threads> julia --project=@. integrate_ephemeris.jl
#Single-threaded:
#julia --project=@. integrate_ephemeris.jl

@show Threads.nthreads()

using PlanetaryEphemeris
using Dates, JLD

#script parameters (TODO: use ArgParse.jl instead)
maxsteps = 1000000
# jd0 = datetime2julian(DateTime(1969,6,28,0,0,0)) #starting time of integration
jd0 = datetime2julian(DateTime(2000,1,1,12)) #starting time of integration
nyears = 2031.0 - year(julian2datetime(jd0))
@show jd0, J2000, jd0-J2000
dense = true #false
dynamics = DE430!
@show dynamics
nast = 343 #16 # number of asteroid perturbers
quadmath = false #true # use quadruple precision
###bodyind = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 17, 18, 19, 20, 21, 22, 23, 25, 27, 28, 30, 40, 41, 46, 55, 62, 73, 113, 115, 277, 322] # SS + 25 ast perturbers
bodyind = 1:(11+16) #1:(11+nast) # body indices in output

# integration parameters
order = 25
abstol = 1.0E-20

#integrator warmup
PlanetaryEphemeris.propagate(1, jd0, nyears, output=false, dense=dense, dynamics=dynamics, nast=nast, quadmath=quadmath, bodyind=bodyind, order=order, abstol=abstol)
println("*** Finished warmup")

# perform full integration
PlanetaryEphemeris.propagate(maxsteps, jd0, nyears, dense=dense, dynamics=dynamics, nast=nast, quadmath=quadmath, bodyind=bodyind, order=order, abstol=abstol)
println("*** Finished full integration")


text:
InputFolder = './Images/CytopacqCells/Movie1/';
OutputFolder = './Results/Images/CytopacqCells/Movie1/';

Fill = 0;
Lbl = 1;

@iA = '*.tif';

@fxg_mLocThr [iA] > [M];
params.Sigmas = [4 4];
params.MeanBox = [50 50];
params.AddThr = 2;
params.MinArea = 0;
params.IgnoreZero = 0;
/endf

@fxm_mModBinWat [M] > [M2];
params.SmallHolesArea = 0;
params.DistStdRad = 5;
params.MinDistLocVar = 9;
/endf

/keep M2 > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
using Dates

@testset "Date.TimeTypes" begin

x, code, vpos, vlen, tlen = Parsers.xparse(Date, "")
@test x === Date(0)
@test code == INVALID | EOF
x, code, vpos, vlen, tlen = Parsers.xparse(Date, "2018-01-01")
@test x === Date(2018, 1, 01)
@test code == OK | EOF
x, code, vpos, vlen, tlen = Parsers.xparse(DateTime, "2018-01-01")
@test x === DateTime(2018, 1, 01)
@test code == OK | EOF
x, code, vpos, vlen, tlen = Parsers.xparse(Time, "01:02:03")
@test x === Time(1, 2, 3)
@test code == OK | EOF

x, code, vpos, vlen, tlen = Parsers.xparse(Date, "")
@test x === Date(0)
@test code == INVALID | EOF
x, code, vpos, vlen, tlen = Parsers.xparse(Date, "\"\"")
@test x === Date(0)
@test code == QUOTED | INVALID | EOF
x, code, vpos, vlen, tlen = Parsers.xparse(Date, "2018-01-01")
@test x === Date(2018, 1, 01)
@test code == OK | EOF
x, code, vpos, vlen, tlen = Parsers.xparse(Date, "\"2018-01-01\"")
@test x === Date(2018, 1, 01)
@test code == QUOTED | OK | EOF
x, code, vpos, vlen, tlen = Parsers.xparse(DateTime, "\"2018-01-01")
@test x === DateTime(2018, 1, 1)
@test code == OK | QUOTED | EOF | INVALID_QUOTED_FIELD
x, code, vpos, vlen, tlen = Parsers.xparse(Date, "\"abcd\"")
@test x === Date(0)
@test code == QUOTED | INVALID | EOF

x, code, vpos, vlen, tlen = Parsers.xparse(Date, "NA", sentinel=["NA"])
@test x === Date(0)
@test code === SENTINEL | EOF
x, code, vpos, vlen, tlen = Parsers.xparse(Date, "\\N", sentinel=["\\N"])
@test x === Date(0)
@test code === SENTINEL | EOF
x, code, vpos, vlen, tlen = Parsers.xparse(Date, "NA2", sentinel=["NA"])
@test x === Date(0)
@test code === SENTINEL | INVALID_DELIMITER | EOF
x, code, vpos, vlen, tlen = Parsers.xparse(Date, "-", sentinel=["-"])
@test x === Date(0)
@test code === SENTINEL | EOF
x, code, vpos, vlen, tlen = Parsers.xparse(Date, "£", sentinel=["£"])
@test x === Date(0)
@test code === SENTINEL | EOF
x, code, vpos, vlen, tlen = Parsers.xparse(Date, "null")
@test x === Date(0)
@test code === INVALID_DELIMITER | EOF
x, code, vpos, vlen, tlen = Parsers.xparse(Date, ",")
@test x === Date(0)
@test code === INVALID | DELIMITED
x, code, vpos, vlen, tlen = Parsers.xparse(Date, "1,")
@test x === Date(1)
@test code === OK | DELIMITED

@test Parsers.parse(DateTime, "1996/Feb/15", Parsers.Options(dateformat="yy/uuu/dd")) === DateTime(1996, 2, 15)
@test Parsers.parse(DateTime, "1996, Jan, 15", Parsers.Options(dateformat="yyyy, uuu, dd")) === DateTime(1996, 1, 15)

end


text:
## Simple random graph generators ##

"""
    rand_geom(p::Real)

Sample a value from geometric distribution with probability of success `p`.
This code is adapted from the [Distributions.jl](https://github.com/JuliaStats/Distributions.jl)
package.
"""
rand_geom(p::Real) = floor(Int, -randexp() / log1p(-p))

"""
    gilbert_fill!(g::Graph, p::Real)

Fill graph by choosing each edge randomly with probability `p`.
"""
function gilbert_fill!(g::UndirectedGraph, p::Real)
    i, j = 1, 1
    n = nodecount(g)
    while true
        j += 1 + rand_geom(p)
        while j > n
            i += 1
            i > n && return
            j = j - n + i
        end
        addedge!(g, i, j)
    end
end

function gilbert_fill!(g::DirectedGraph, p::Real)
    i, j = 1, 1
    n = nodecount(g)
    while true
        j_ = j + 1 + rand_geom(p)
        j = (j < i && j_ >= i) ? j_ + 1 : j_
        while j > n
            i += 1
            i > n && return
            j -= n
            j += j >= i
        end
        addedge!(g, i, j)
    end
end

"""
    graph_gilbert(n::Integer, p::Real)

Generate a random graph according to Gilbert model with `n` nodes
and expected density `p`.

An independent Bernoulli trial is simulated for each possible edge in the graph
and it is added with probability `p`.
Note that in this case the number of edges is different each time, but the
matiematical expectation for the graph's density is `p`.

# References

Gilbert, E.N. (1959). "Random Graphs". Annals of Mathematical Statistics. 30 (4): 1141–1144.
doi:10.1214/aoms/1177706098.

Batagelj, V. and Brandes, U. (2005) "Efficient generation of large random networks".
Physical Review E. 71 (3).
"""
function graph_gilbert(n::Integer, p::Real; kwargs...)
    @assert(0 <= p <= 1, "edge probability (p) should be in [0, 1] interval")
    g = Graph(n; kwargs...)
    p > eps() && gilbert_fill!(g, p)
    return g
end

"""
    graph_erdos_renyi(n::Integer, m::Integer)

Generate an [Erdős–Rényi random graph](https://en.wikipedia.org/wiki/Erdős–Rényi_model)
with `n` nodes and `m` edges.

From mathematical viewpoint, the graph is picked at random with uniform probability
from a set of all possible graphs having `n` nodes and `m` edges.
In practice, an empty graph with `n` nodes is created first and
`m` edges get added between the randomly chosen pairs of its nodes.

# References

Erdős, P.; Rényi, A. (1959). "On Random Graphs. I". Publicationes Mathematicae. 6: 290–297.
"""
function graph_erdos_renyi(n::Integer, m::Integer; kwargs...)
    g = Graph(n; kwargs...)
    if m < maxedgecount(g)
        while edgecount(g) < m
            x, y = rand(1:n), rand(1:n)
            if x == y || hasedge(g, x, y)
                continue
            else
                addedge!(g, x, y)
            end
        end
    elseif m == maxedgecount(g)
        gilbert_fill!(g, 1.)
    else
        error("too many edges (m) for Erdős–Rényi graph with n nodes")
    end
    return g
end

@deprecate(erdosrenyi, graph_erdos_renyi)
@deprecate(graph_erdosrenyi, graph_erdos_renyi)

"""
    graph_random(n::Integer, m::Integer)
    graph_random(n::Integer, p::AbstractFloat)

A shorthand for creating random networks. There are two similar but distinct models:

* `(n, m)` model, also known as Erdős–Rényi model. It randomly selects a
  graph from a set of all possible graphs which have `n` nodes and `m` edges.

* `(n, p)` model, introduced by Gilbert. For a graph on `n` nodes,
  it includes each edge with probability `p` independently of the others.
  This gives a graph with mathematical expectation of density equal to `p`.

This method is just a shorthand for [`graph_erdos_renyi`](@ref) and
[`graph_gilbert`](@ref), which implement these models and are chosen automatically
based on the argument types.

If you need a graph of high density, use `(n, p)` model as it is much faster.
If you want an exact number of edges with lower density, use `(n, m)` model.

# Examples

```julia-repl
julia> g = graph_random(4, 1.)
4-node 12-edge directed multigraph

julia> g = graph_random(4, 1)
4-node 1-edge directed multigraph
```

# References

[Wikipedia](https://en.wikipedia.org/wiki/Erdős–Rényi_model)
"""
graph_random(n::Integer, m::Integer; kwargs...) =
    graph_erdos_renyi(n::Integer, m::Integer; kwargs...)
graph_random(n::Integer, p::AbstractFloat; kwargs...) =
    graph_gilbert(n::Integer, p::AbstractFloat; kwargs...)


text:
function get_domain_problem_objects(domain_path, problem_path)
    domain = load_domain(domain_path)
    problem = load_problem(problem_path)

    return domain, problem
end


text:
using CompScienceMeshes
using BEAST

using Test

m = meshrectangle(1.0, 1.0, 0.5, 3)
nodes = skeleton(m,0)
srt_bnd_nodes = sort.(skeleton(boundary(m),0))
int_nodes = submesh(nodes) do node
    !(sort(node) in srt_bnd_nodes)
end

@test length(int_nodes) == 1

L0 = BEAST.lagrangec0d1(m, int_nodes)
@test numfunctions(L0) == 1
@test length(L0.fns[1]) == 6
@test length(m) == 8

Lx = BEAST.lagrangecxd0(m)
@test numfunctions(Lx) == 8

Id = BEAST.Identity()
Q1, st1 = BEAST.allocatestorage(Id, L0, Lx, Val{:bandedstorage}, BEAST.LongDelays{:ignore})
BEAST.assemble_local_matched!(Id, L0, Lx, st1)

Q2, st2 = BEAST.allocatestorage(Id, L0, Lx, Val{:bandedstorage}, BEAST.LongDelays{:ignore})
BEAST.assemble_local_mixed!(Id, L0, Lx, st2)

@test isapprox(Q1, Q2, atol=1e-8)

RT = raviartthomas(m)
BC = buffachristiansen(m)

Q1, st1 = BEAST.allocatestorage(Id, BC, RT, Val{:bandedstorage}, BEAST.LongDelays{:ignore})
BEAST.assemble_local_refines!(Id, BC, RT, st1)

Q2, st2 = BEAST.allocatestorage(Id, BC, RT, Val{:bandedstorage}, BEAST.LongDelays{:ignore})
BEAST.assemble_local_mixed!(Id, BC, RT, st2)
@test isapprox(Q1, Q2, atol=1e-8)


text:
using Test
using Transit
using JSON

import DataStructures: list

# Convert the value given to transit and then read
# back in the resulting JSON. Not a round trip, since
# we are just reading back in the raw JSON.
function round_trip(x, verbose=false)
    io = IOBuffer()
    Transit.write(io, x)
    seekstart(io)
    Transit.parse(io)
end

function test_round_trip(value, verbose=false)
    x = round_trip(value, verbose)
    x == value
end

@test test_round_trip(1)
@test test_round_trip(1.0)
@test test_round_trip("hello")
@test test_round_trip("~hello")
@test test_round_trip(true)
@test test_round_trip(false)
@test test_round_trip(2//3)
@test test_round_trip([1,2,3])
@test test_round_trip(list(1,2,"hello"))

@test test_round_trip([:hello, :hello, :hello])
@test test_round_trip([:aaaa, :bbbb, :aaaa, :bbbb])

symbols = [Transit.TSymbol("hello"), Transit.TSymbol("hello"), Transit.TSymbol("hello")]

@test test_round_trip(symbols)
@test test_round_trip([:aaaa, :bbbb, :aaaa, :bbbb])

@test test_round_trip([:aaaa, :bbbb, :aaaa, :bbbb], true)

@test test_round_trip(Dict([("aaa", 11), ("bbb", 22)]), true)

@test test_round_trip(Transit.TSet([1,2,3]))

# @test test_round_trip(Transit.TSet(Any[Transit.TSet([1,2]), Transit.TSet([2,3])]))


text:
plusnumber(t::Int) = t>0 ? "+$t" : "$t"

#treat thsi table as a one-off since it has no common elements w/ other tables
function maketable7(bb::DataFrame)
  T7_PANEL_A_COLKEYS = [:StdDev, :Mean, :Range, :Return, :AbsReturn]
  T7_PANEL_A_NAME::String = "t7panela-$(REPLICATION_TYPE[])-$(DIST_TYPE[])-$(OUT_SUFFIX[])"
  T7_PANEL_A_KEYROW_RANGES::Dict = Dict(-3=>-6:-3, 0=>-2:2, 3=>3:6)
  local days::Int = SEO_DAYS
  local Fret::Symbol = SEO_RET
  local minallowed::Int = SEO_MIN_ALLOWED

  requirefield!(bb, Fret, minallowed)
  deleterows!(bb, (!).(bb.oneex))

  table7panela(bb,
      days=days,
      tablename= T7_PANEL_A_NAME,
      colkeys = T7_PANEL_A_COLKEYS,
      keyrowranges = T7_PANEL_A_KEYROW_RANGES,
      Fret=Fret
      )
end

function requirefield!(df::DataFrame, Fret::Symbol, minallowed::Int; Fgroup=:permno)
  df.tokeep = trues(size(df,1))
  for sdf ∈ groupby(df, :permno)
    tokeep::Bool = sum((!ismissing).(sdf[!,Fret])) ≥ minallowed
    sdf.tokeep .= tokeep
  end

  deleterows!(df, (!).(df.tokeep))
  select!(df, Not(:tokeep))

  return nothing
end

function table7panela(bb::AbstractDataFrame;
    days = SEO_DAYS,
    tablename::String = error("tablename is required"),
    tablepath::String = TABLE_PATH,
    decimals::Int = DEFAULT_DECIMALS,
    colkeys::Vector{Symbol} = error("colkeys is required"),
    colnames::Vector{Vector{String}} = [["", "Daily Cross-", "", ""],
      ["Event", "Sectional Returns", "", "Average"],
      ["Day"; (string).(colkeys);]],
    widthcolnames::Vector{Vector{Int}} = [[1,2,1,2], [1,2,1,2], ones(Int, length(colnames[3]))],
    descrownames::Vector{String} = (plusnumber).(-days:days),
    keyrowranges::Dict = error("keyrow ranges is required"),
    Fret::Symbol = error("Fret is required"),
    )

  local parallel::Bool = PARALLEL[] && (!DISABLE_FM_PARALLEL)
  #bb = bb[completecases[]]

  #################################
  #basic lookups
  i2day::Vector{Int} = collect(-days:days)
  sbbs::GroupedDataFrame = groupby(bb, :day)
  bbidx::Dict = Dict((i::Int->sbbs[i][1,:day]=>sbbs[i]).(1:length(sbbs)))

  #preallocation and housekeeping

  Ncols::Int = sum(widthcolnames[1])-1
  Nrows::Int = 2*days+1
  desccontent::Vector{Vector{String}} = (i->Vector{String}(undef, Ncols)).(1:Nrows)

  #create an index to make this part of the code more readable
  colidx::Dict = Dict(colkey=>i for (i, colkey) ∈ enumerate(colkeys))

  #now populate the the table
  for (i,rowkey) ∈ enumerate(descrownames)
    day::Int = i2day[i] #get the period offset
    sbb::SubDataFrame = bbidx[day]

    #compute the following stats for each row
    daystd::Float64 = std(skipmissing(sbb[!, Fret]))*100
    daymean::Float64 = mean(skipmissing(sbb[!, Fret]))*100

    desccontent[i][colidx[:StdDev]] = num2str(daystd, decimals)
    desccontent[i][colidx[:Mean]] = num2str(daymean, decimals)

    if haskey(keyrowranges, day)

      dayₗ::Int = minimum(keyrowranges[day])
      dayₕ::Int = maximum(keyrowranges[day])
      desccontent[i][colidx[:Range]] = "\\text{$(plusnumber(dayₗ)) to $(plusnumber(dayₕ))}"

      ssbb = view(bb, (r::DataFrameRow->r.day ∈ keyrowranges[day]).(eachrow(bb)), [:permno, Fret])
      #=aggbb::DataFrame = aggregate(groupby(ssbb, :permno),
        v::AbstractVector{MFloat64}->prod(1 .+ skipmissing(v)))
      rets::Vector{Float64} = aggbb[!,Symbol(Fret, :_function)] .- 1

      retmean::Float64 = mean(skipmissing(rets))*100
      retabsmean::Float64 = mean(skipmissing((abs).(rets)))*100=#

      retmean::Float64 = mean(skipmissing(ssbb[!,Fret]))*100
      retabsmean::Float64 = mean(skipmissing((abs).(ssbb[!,Fret])))*100

      desccontent[i][colidx[:Return]] = num2str(retmean, decimals)
      desccontent[i][colidx[:AbsReturn]] = num2str(retabsmean, decimals)
    else
      desccontent[i][colidx[:Range]] = ""
      desccontent[i][colidx[:Return]] = ""
      desccontent[i][colidx[:AbsReturn]] = ""
    end
  end

  alignmentstring::String = "rrr c rr"

  #println(desccontent)

  #display(desccontent)
  #construct the tex table
  tt::String = textable(colnames=colnames,
    descrownames=descrownames,
    desccontent=desccontent,
    alignmentstring=alignmentstring,
    widthcolnames=widthcolnames,
    rowlabelheader=true)

  writetextable(tt, path=tablepath, outname="$tablename.tex")
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Distributed
@everywhere using Mamba

## Data
line = Dict{Symbol, Any}(
  :x => [1, 2, 3, 4, 5],
  :y => [1, 3, 3, 3, 5]
)

line[:xmat] = [ones(5) line[:x]]

## Model Specification
model = Model(
  y = Stochastic(1,
    (xmat, beta, s2) -> MvNormal(xmat * beta, sqrt(s2)),
    false
  ),
  beta = Stochastic(1, () -> MvNormal(2, sqrt(100))),
  s2 = Stochastic(() -> InverseGamma(0.01, 0.01))
)

## Initial Values
inits = [
  Dict{Symbol, Any}(
    :y => line[:y],
    :beta => rand(Normal(0, 1), 2),
    :s2 => rand(Gamma(1, 1))
  )
  for i in 1:3
]

## Tuning Parameters
scale1 = [0.5, 0.25]
summary1 = identity
eps1 = 0.5

scale2 = 0.5
summary2 = x -> [mean(x); sqrt(var(x))]
eps2 = 0.1

## User-Defined Sampling Scheme
scheme = [
  ABC(:beta, scale1, summary1, eps1, kernel=Normal, maxdraw=100, nsim=3),
  ABC(:s2,   scale2, summary2, eps2, kernel=Epanechnikov, maxdraw=100, nsim=3)
]
setsamplers!(model, scheme)

## MCMC Simulation with Approximate Bayesian Computation
sim = mcmc(model, line, inits, 10000, burnin=1000, chains=3)
describe(sim)


text:
export OceanBC,
    VelocityBC,
    VelocityDragBC,
    TemperatureBC,
    Impenetrable,
    Penetrable,
    NoSlip,
    FreeSlip,
    KinematicStress,
    Insulating,
    TemperatureFlux

using ..BalanceLaws
using ..DGMethods.NumericalFluxes

"""
    OceanBC(velocity    = Impenetrable(NoSlip())
            temperature = Insulating())

The standard boundary condition for OceanModel. The default options imply a "no flux" boundary condition.
"""
Base.@kwdef struct OceanBC{M, T}
    velocity::M = Impenetrable(NoSlip())
    temperature::T = Insulating()
end

abstract type VelocityBC end
abstract type VelocityDragBC end
abstract type TemperatureBC end

"""
    Impenetrable(drag::VelocityDragBC) :: VelocityBC

Defines an impenetrable wall model for velocity. This implies:
  - no flow in the direction normal to the boundary, and
  - flow parallel to the boundary is subject to the `drag` condition.
"""
struct Impenetrable{D <: VelocityDragBC} <: VelocityBC
    drag::D
end

"""
    Penetrable(drag::VelocityDragBC) :: VelocityBC

Defines an penetrable wall model for velocity. This implies:
  - no constraint on flow in the direction normal to the boundary, and
  - flow parallel to the boundary is subject to the `drag` condition.
"""
struct Penetrable{D <: VelocityDragBC} <: VelocityBC
    drag::D
end

"""
    NoSlip() :: VelocityDragBC

Zero velocity at the boundary.
"""
struct NoSlip <: VelocityDragBC end

"""
    FreeSlip() :: VelocityDragBC

No surface drag on velocity parallel to the boundary.
"""
struct FreeSlip <: VelocityDragBC end

"""
    KinematicStress(stress) :: VelocityDragBC

Applies the specified kinematic stress on velocity normal to the boundary.
Prescribe the net inward kinematic stress across the boundary by `stress`,
a function with signature `stress(problem, state, aux, t)`, returning the flux (in m²/s²).
"""
struct KinematicStress <: VelocityDragBC end

"""
    Insulating() :: TemperatureBC

No temperature flux across the boundary
"""
struct Insulating <: TemperatureBC end

"""
    TemperatureFlux(flux) :: TemperatureBC

Prescribe the net inward temperature flux across the boundary by `flux`,
a function with signature `flux(problem, state, aux, t)`, returning the flux (in m⋅K/s).
"""
struct TemperatureFlux <: TemperatureBC end

"""
    ocean_boundary_state!(nf, boundaries::Tuple, ::HBModel,
                          Q⁺, A⁺, n, Q⁻, A⁻, bctype)
applies boundary conditions for the first-order and gradient fluxes
dispatches to a function in OceanBoundaryConditions.jl based on bytype defined by a problem such as SimpleBoxProblem.jl
"""
@generated function ocean_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    boundaries::Tuple,
    ocean,
    Q⁺,
    A⁺,
    n,
    Q⁻,
    A⁻,
    bctype,
    t,
    args...,
)
    N = fieldcount(boundaries)
    return quote
        Base.Cartesian.@nif(
            $(N + 1),
            i -> bctype == i, # conditionexpr
            i -> ocean_boundary_state!(
                nf,
                boundaries[i],
                ocean,
                Q⁺,
                A⁺,
                n,
                Q⁻,
                A⁻,
                t,
            ), # expr
            i -> error("Invalid boundary tag")
        ) # elseexpr
        return nothing
    end
end

"""
    ocean_boundary_state!(nf, boundaries::Tuple, ::HBModel,
                          Q⁺, A⁺, D⁺, n, Q⁻, A⁻, D⁻, bctype)
applies boundary conditions for the second-order fluxes
dispatches to a function in OceanBoundaryConditions.jl based on bytype defined by a problem such as SimpleBoxProblem.jl
"""
@generated function ocean_boundary_state!(
    nf::NumericalFluxSecondOrder,
    boundaries::Tuple,
    ocean,
    Q⁺,
    D⁺,
    A⁺,
    n,
    Q⁻,
    D⁻,
    A⁻,
    bctype,
    t,
    args...,
)
    N = fieldcount(boundaries)
    return quote
        Base.Cartesian.@nif(
            $(N + 1),
            i -> bctype == i, # conditionexpr
            i -> ocean_boundary_state!(
                nf,
                boundaries[i],
                ocean,
                Q⁺,
                D⁺,
                A⁺,
                n,
                Q⁻,
                D⁻,
                A⁻,
                t,
            ), # expr
            i -> error("Invalid boundary tag")
        ) # elseexpr
        return nothing
    end
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using CVChannel
using Test

"""
This script shows the multiplicativity of the Siddhu channel with itself.
"""

println("In this script we look at the parallel communication value of the Siddhu channel with itself.")
println("Since we know cv(N) is 2 always, all we need to do is check that cvPPT(N ⊗ N) = 4 always.")
@testset "Verify that cv is 2 and cv(N ⊗ N) is multiplicative" begin
    scan_range =[0:0.01:0.5;]
    is_mult = true
    ctr = 1
    println("Now starting to scan over s.")
    for s_id in scan_range
        ctr % 10 == 0 ? println("Now on ", ctr, " of ", length(scan_range), " points.") : nothing
        #To save time I don't check the single copy cv
        # parallel channel uses
        sidchan(X) = siddhuChannel(X,s_id)
        sid_chan= Choi(sidchan,3,3)
        par_cv, = pptCV(parChoi(sid_chan, sid_chan), :dual)
        !isapprox(par_cv,4,atol=3e-6) ? is_mult = false : nothing
        ctr += 1
    end
    @test is_mult
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
module EDClust

using Distributions
using SpecialFunctions
using StatsBase
using ProgressMeter

include("EMPolya.jl")
include("MMPolya.jl")
include("fitPolya.jl")

export fitPolya

end


text:
import Base: *, ==, +, isless, iszero, zero, iterate, eltype, size, length, getindex, setindex!,
push!, empty!, show

# Helper functions
include("utils.jl")

# Coefficients
include("barecoefficient.jl")
include("coefficient.jl")
include("coefficientsum.jl")

# Operators
include("majoranaproduct.jl")

# Terms
include("term.jl")
include("termsum.jl")

max_stars_repo_name: aormorningstar/SymbolicMajoranas

max_issues_repo_name: aormorningstar/SymbolicMajoranas

max_forks_repo_name: aormorningstar/SymbolicMajoranas


text:
# pseudospectral radius and abscissa
using Pseudospectra, Test

@testset "PS radius & abscissa" begin
    A = Pseudospectra.grcar(100)

    opts = Dict{Symbol,Any}()

    # also check that this gets the eigenvalues
    ps_data = new_matrix(A,opts)
    topeig = maximum(imag(ps_data.ps_dict[:ews]))
    @test isapprox(topeig, 2.2617668, atol=1e-6)

    # reference values from Guglielmi & Overton (2012)
    r,zr = psa_radius(A,1e-4,ps_data.ps_dict[:ews])
    @test isapprox(r, 2.85216, atol=1e-5)
    α,za = psa_abscissa(A,1e-4,ps_data.ps_dict[:ews])
    @test isapprox(α, 2.41276, atol=1e-5)

    # Note: if we ever interface to a Hamiltonian eigensolver, we could use
    # the ϵ=1e-8 case in S&P.
end


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# begin enum bk_termtype
const bk_termtype = UInt32
const SOURCE = (UInt32)(0)
const SINK = (UInt32)(1)
# end enum bk_termtype

error = Cvoid

const bk_error = Ptr{Cvoid}

graph_int = Cvoid

const bk_graph_int = Ptr{Cvoid}

graph_short = Cvoid

const bk_graph_short = Ptr{Cvoid}

graph_float = Cvoid

const bk_graph_float = Ptr{Cvoid}

graph_double = Cvoid

const bk_graph_double = Ptr{Cvoid}

max_stars_repo_name: Gnimuc/BKmaxflow

max_issues_repo_name: Gnimuc/BKmaxflow

max_forks_repo_name: Gnimuc/BKmaxflow


text:
function code_typed(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = mcjob(func, types; kwargs...)
    GPUCompiler.code_typed(job; kwargs...)
end

function code_warntype(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = mcjob(func, types; kwargs...)
    GPUCompiler.code_warntype(io, job; kwargs...)
end

function code_llvm(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = mcjob(func, types; kwargs...)
    GPUCompiler.code_llvm(io, job; kwargs...)
end

function code_native(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = mcjob(func, types; kwargs...)
    GPUCompiler.code_native(io, job; kwargs...)
end

# aliases without ::IO argument
for method in (:code_warntype, :code_llvm, :code_native)
    @eval begin
        $method(@nospecialize(func), @nospecialize(types); kwargs...) =
            $method(stdout, func, types; kwargs...)
    end
end


text:
function value(obj::TricubicFunction, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(obj, "value", jdouble, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end


text:
export find_index

import Base: show

"""
    find_index(conditionFn::F) where { F <: Function }

Creates a find operator, which emits only the index of the first value emitted by the source Observable that meets some condition.
Indices are 1-based.

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Arguments
- `conditionFn::F`: condition function with `(data::T) -> Bool` signature

# Examples
```jldoctest
using Rocket

source = from([ 0, 1, 2, 3, 4, 5, 6, 7, 8 ])
subscribe!(source |> find_index((d) -> d !== 0 && d % 2 == 0), logger())
;

# output

[LogActor] Data: 3
[LogActor] Completed
```

See also: [`find`](@ref), [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
find_index(conditionFn::F) where { F <: Function } = enumerate() + find(d -> conditionFn(d[2])) + map(Int, d -> d[1])


text:
@testset "Header" begin
    @testset "Header Enums" begin
        let
            ids = vcat(-12345, 1:50, 52:97, 103)
            enums = collect(instances(HeaderEnum))
            @test enums == map(HeaderEnum, ids)
        end
    end
end


text:
module ChainMap
import MacroTools
import DataStructures
import NullableArrays

include("nonstandard.jl")
include("annotate.jl")
include("chain.jl")
include("lambda.jl")
include("LazyCall.jl")
include("recursive.jl")
include("unweave.jl")
include("all.jl")

end


text:
using Knet
include("ImageEncoder.jl")
include("Generator.jl")
include("Evaluations.jl")
# Sentence encoder is not yet complete so we are using a random vector
# I haven't implemented a loop yet since everything is not complete, but the following code implements a forward pass

x = param(rand(128,128,3,5)) # input image batch
imgFeatures = ImageEncoder(x) # 16 x 16 x 512 x 5
pooledFeatures = pool_features(imgFeatures) # 1 x 1 x 512 x 5
objDetections = obj_detector(pooledFeatures) # 58 x 5

y = param(rand(rnn_out_dim,batch_size)) # Encoded sentence, we start with random since it is not fully implemented
z = param(randn(noise_dim,batch_size)) # Noise to input to the the generator

generatedImage = forward(y,z,imgFeatures) # 128 x 128 x 3 x 5

generatorLoss = AdversarialGenerator(generatedImage) # Evaluation of generator


text:
import Pkg
println(@__DIR__)
if !occursin("DFQExamples.jl", Base.active_project())
    Pkg.add(url="https://github.com/angus-lewis/DiscretisedFluidQueues")
    # Pkg.rm("DiscretisedFluidQueues")
    Pkg.activate(@__DIR__)
    #Pkg.develop(url=(@__DIR__)*"/../DiscretisedFluidQueues.jl")
    #CPkg.develop(url=(@__DIR__)*"/../FluidQueues.jl")
    Pkg.develop(url=(@__DIR__)*"/../FluidFluidQueues.jl")
    Pkg.instantiate()
    display(!occursin("DFQExamples.jl", Base.active_project()))
end

using CSV
using DiscretisedFluidQueues
using DataFrames
using Distributed
import Distributions
using FluidQueues
using FluidFluidQueues
using Plots
using LinearAlgebra
using JLD2
using JSON
using Random
using SparseArrays
using StableRNGs
using Statistics

include((@__DIR__)*"/error_metrics.jl")
include((@__DIR__)*"/helper_functions.jl")

plotlyjs()


text:
include("../../common/numder/src/numder_module.jl")

module nonlineq1d
using ..numder

struct NonlineqTolerance
    max_iter::Int
    abs_step_size_tol::Float64
    rel_step_size_tol::Float64
end

function NonlineqTolerance(;
    max_iter = 10_000,
    abs_step_size_tol = 1e-15,
    rel_step_size_tol = 1e-10)

    return NonlineqTolerance(
        max_iter,
        abs_step_size_tol,
        rel_step_size_tol)
end

function __check_bracketing(fcn, x0, x1)
    f0    = fcn(x0)
    f1    = fcn(x1)
    if f0 * f1 > 0
        error("Initial bracket doesnt bracket a root with opposite sings")
    end

    return (f0, f1)
end

function __check_step_size(x0, x1, options)

    break_flag = false

    abs_tol = options.abs_step_size_tol
    rel_tol = options.rel_step_size_tol * max(abs(x0), abs(x1))

    dx = x1 - x0
    if abs(dx) <= max(abs_tol, rel_tol)
        break_flag = true
    end

    return break_flag
end

struct SolverStatistics
    iter::Int
end

function bisection(fcn, x0, x1; options = NonlineqTolerance() )

    (f0, f1) = __check_bracketing(fcn, x0, x1)

    x_sol = NaN64
    iter = 0
    for outer iter = 1:options.max_iter

        x_mid = 0.5 * x0 + 0.5 * x1

        break_flag = __check_step_size(x0, x1, options)
        if break_flag
            x_sol = x_mid
            break
        end

        f_mid = fcn(x_mid)

        if f_mid * f0 < 0.0
            x0, x1 = x0, x_mid
            f0, f1 = f0, f_mid
        elseif f_mid * f1 < 0.0
            x0, x1 = x_mid, x1
            f0, f1 = f_mid, f1
        else
            # exactly got to 0
            x_sol = x_mid
            break
        end

    end

    stat = SolverStatistics(iter)
    return (x_sol, stat)
end

function regula_falsi(fcn, x0, x1; options = NonlineqTolerance())

    (f0, f1) = __check_bracketing(fcn, x0, x1)

    x_sol = NaN64

    iter = 0
    for outer iter = 1:options.max_iter

        break_flag = __check_step_size(x0, x1, options)
        if break_flag
            x_sol = 0.5 * x0 + 0.5 * x1
            break
        end

        x_lin = x0 - f0 * (x1 - x0) / (f1 - f0)
        f_lin = fcn(x_lin)

        x_mid = 0.5 * x0 + 0.5 * x1
        f_mid = fcn(x_mid)

        if f_lin * f_mid < 0.0
            # take this interval
            x0, x1 = x_lin, x_mid
            f0, f1 = f_lin, f_mid

            continue
        end

        if f_mid * f0 < 0.0

            if x0 <= x_lin <= x_mid
                # f_lin has the same sign as f_mid
                x0, x1 = x0, x_lin
                f0, f1 = f0, f_lin
            else # not in the interval, no improvement achievable
                x0, x1 = x0, x_mid
                f0, f1 = f0, f_mid
            end

        elseif f_mid * f1 < 0.0
            x0, x1 = x_mid, x1
            f0, f1 = f_mid, f1

            if x_mid <= x_lin <= x1
                # f_lin has the same sign as f_mid
                x0, x1 = x_lin, x1
                f0, f1 = f_lin, f1
            else # not in the interval, no improvement achievable
                x0, x1 = x_mid, x1
                f0, f1 = f_mid, f1
            end

        else
            # exactly got to 0
            x_sol = x_mid
            break
        end

    end # for

    stat = SolverStatistics(iter)
    return (x_sol, stat)
end

function regula_falsi_ilinois(fcn, x0, x1; options = NonlineqTolerance())

    (f0, f1) = __check_bracketing(fcn, x0, x1)
    x_prev   = NaN64
    x_pprev  = NaN64
    x_sol    = NaN64

    iter = 0
    for outer iter = 1:options.max_iter

        break_flag = __check_step_size(x0, x1, options)
        if break_flag
            x_sol = 0.5 * x0 + 0.5 * x1
            break
        end

        if x_prev == x_pprev
            if x0 == x_prev
                y0, y1 = f0 * 0.5, f1
            else
                y0, y1 = f0, f1 * 0.5
            end
        else
            y0, y1 = f0, f1
        end

        x_lin = x0 - y0 * (x1 - x0) / (y1 - y0)
        f_lin = fcn(x_lin)

        if f_lin * f0 < 0.0
            x0, x1  = x0, x_lin
            f0, f1  = f0, f_lin
            x_pprev = iter >= 2 ? x_prev : NaN64
            x_prev  = x0
        elseif f_lin * f1 < 0.0
            x0, x1  = x_lin, x1
            f0, f1  = f_lin, f1
            x_pprev = iter >= 2 ? x_prev : NaN64
            x_prev  = x1
        else
            # exactly got to 0
            x_sol = x_lin
            break
        end

    end

    stat = SolverStatistics(iter)
    return (x_sol, stat)
end

function newton(fcn, x0, x1; options = NonlineqTolerance())

    xk    = 0.5 * x0 + 0.5 * x1
    x_sol = NaN64
    iter  = 0
    for outer iter = 1:options.max_iter

        fk = fcn(xk)
        df = numder.numdiff_fw(fcn, xk, sqrt(eps(xk)), fk)

        dx  = - fk / (df + eps(df))
        xk1 = xk + dx

        break_flag = __check_step_size(xk, xk1, options)
        if break_flag
            x_sol = xk1
            break
        end

        xk = xk1
    end

    stat = SolverStatistics(iter)
    return (x_sol, stat)
end

end # module

module ntest
using ..nonlineq1d

function test_1()
    fcn(x) = x^2 - 4

    x0 = 0.1
    x1 = 5

    x_bis      = nonlineq1d.bisection(fcn, x0, x1)
    x_regfal   = nonlineq1d.regula_falsi(fcn, x0, x1)
    x_regfalil = nonlineq1d.regula_falsi_ilinois(fcn, x0, x1)
    x_new      = nonlineq1d.newton(fcn, x0, x1)

end

fcn(x) = cos(x)

x0 = pi/2 - pi/4
x1 = pi/2 + pi/4

x_bis      = nonlineq1d.bisection(fcn, x0, x1)
x_regfal   = nonlineq1d.regula_falsi(fcn, x0, x1)
x_regfalil = nonlineq1d.regula_falsi_ilinois(fcn, x0, x1)
x_new      = nonlineq1d.newton(fcn, x0, x1)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# graddesc.jl
module GradDesc

using InfDecomp_Base

using CPUTime
export my_gradient_descent

struct Solution_Stats
    optimum     :: Float64
    q           :: Vector{Float64}
    nm_pr∇      :: Float64
    max_η       :: Float64
    status      :: Symbol  # can be one of   :grad0 :step0 :iter
    time        :: Float64
end

# Compute the projection operator onto the image of the matrix B.
# Uses SVD, where singular values less than `eps` are considered 0
# Requires that B be matrix with AT LEAST AS MANY ROWS AS COLUMNS !!!!
function compute_projector( B :: AbstractMatrix{Float64}, eps::Float64=1.e-10 ) :: Matrix{Float64}
    (m,n) = size(B)
    # @assert m >= n "compute_projector(B) requires that B has at least as many rows as columns"

    # SVD
    B_full = full(B)
    F = LinAlg.svdfact!(B_full, thin=false)
    NZ = [ j    for j in 1:min(m,n) if F[:S][j] > eps ]

    # partial isometry:
    PI = view(F[:U], 1:m, NZ)

    P = PI*PI'

    @assert P*P≈P "compute_projector(B): Something went wrong -- P^2 != P :("
    @assert P'≈P  "compute_projector(B): Something went wrong -- P^T != P :("

    return P
end #^ compute_projector()

function initial_interior_point(e::My_Eval, q::Vector{Float64}) :: Void
    for x = 1:e.n_x
        p_x = 0.
        for y = 1:e.n_y
            for z = 1:e.n_z
                p_x += e.prb_xyz[x,y,z]
            end #^ for z
        end #^ for y

        for y = 1:e.n_y
            for z = 1:e.n_z
                i = e.varidx[x,y,z]
                if i>0
                    q[i] = e.marg_xy[x,y] * e.marg_xz[x,z] / p_x
                end #^ if ∃i
            end #^ for z
        end #^ for y
    end #^ for x
    ;
end

function my_gradient_descent(e::My_Eval;
                             max_iter        :: Int64    =1000,
                             eps_grad        :: Float64  =1.e-20,
                             eps_steplength  :: Float64  =1.e-20,
                             stepfactor      :: Float64  =.1        )   :: Solution_Stats
    # Function for steplength
    steplength(q::Float64, g::Float64) = (  g > 0 ?   q/g   :  Inf  )

    CPUtic()

    # Definitions of vectors
    q_0  :: Vector{Float64} = zeros(e.n)  # standard interior feasible point (also initial q)
    q    :: Vector{Float64} = zeros(e.n)  # iterate interior feasible point
    ∇    :: Vector{Float64} = zeros(e.n)  # gradient
    pr∇  :: Vector{Float64} = zeros(e.n)  # projected gradient (onto tangent space)
    P    :: Matrix{Float64} = eye(e.n) - compute_projector(e.Gt) # projection operator (onto lin space)

    # initial solution
    initial_interior_point(e,q_0)

    # main loop
    local nm_pr∇ :: Float64
    local max_η  :: Float64
    local status = :iter
    local iter :: Int64
    local best_obj_val :: Float64 = Inf # stores the best objective value found
    local q_best :: Vector{Float64} = zeros(e.n) # stores the best objective value corresponding input
    local nm_pr∇_best :: Float64
    local iter_best :: Int64
    local max_η_best :: Float64
    q .= q_0
    for iter = 1:max_iter
        # compute gradient
        InfDecomp_Base.∇f(e,∇,q,Float64(0.))

        # project gradient onto tangent space
        pr∇ .=  P*∇

        max_η  = -1.
        nm_pr∇ = norm(pr∇)
        if nm_pr∇ ≤ eps_grad
            status = :grad0
            break
        end

        # max steplength which retains feasibility
        max_η = Inf
        for i in 1:e.n
            max_η = min(  max_η ,  steplength(q[i],pr∇[i])  )
        end

        # max eigenvalue of the Hessian
        if max_η ≤ eps_steplength
            status = :step0
            break
        end

        # check if the objective function is better
        # if so, copy the best feasible q into q_best
        obj_val = -condEntropy(e,q,Float64(0.))
        if obj_val <= best_obj_val
            best_obj_val = obj_val
            q_best .= q
            nm_pr∇_best = nm_pr∇
            iter_best = iter
            max_η_best = max_η
        end

        if iter%10==1
            @show iter nm_pr∇ max_η
            @show q
            @show pr∇
            @show obj_val
        end

        q .-= (stepfactor*min(1.,max_η)) .* pr∇

    end #^ for --- main loop

    tm = CPUtoc()

    println("Terminated with")
    @show iter nm_pr∇ max_η best_obj_val q_best nm_pr∇_best iter_best

    return Solution_Stats(best_obj_val, q_best, nm_pr∇_best, max_η_best, status, tm)
end #^ my_gradient_descent()

end # module

max_stars_repo_name: dot-at/newpid

max_issues_repo_name: dot-at/newpid

max_forks_repo_name: dot-at/newpid


text:
@testset "Woodbury" begin
seed!(123)
n = 5

d = 1 .+ rand(n)
dl = -rand(n-1)
du = -rand(n-1)
v = randn(n)
B = randn(n,2)

U = randn(n,2)
V = randn(2,n)
C = randn(2,2)

for elty in (Float32, Float64, ComplexF32, ComplexF64, Int)
    if elty == Int
        seed!(61516384)
        d = rand(1:100, n)
        dl = -rand(0:10, n-1)
        du = -rand(0:10, n-1)
        v = rand(1:100, n)
        B = rand(1:100, n, 2)

        # Woodbury
        U = rand(1:100, n, 2)
        V = rand(1:100, 2, n)
        C = rand(1:100, 2, 2)
    else
        d = convert(Vector{elty}, d)
        dl = convert(Vector{elty}, dl)
        du = convert(Vector{elty}, du)
        v = convert(Vector{elty}, v)
        B = convert(Matrix{elty}, B)
        U = convert(Matrix{elty}, U)
        V = convert(Matrix{elty}, V)
        C = convert(Matrix{elty}, C)
    end
    ε = eps(abs2(float(one(elty))))
    T = Tridiagonal(dl, d, du)

    # Matrix for A
    W = Woodbury(T, U, C, V)
    @test size(W, 1) == n
    @test size(W) == (n, n)
    @test axes(W) === (Base.OneTo(n), Base.OneTo(n))
    @test axes(W, 1) === Base.OneTo(n)
    F = Matrix(W)
    @test F ≈ T + U*C*V
    @test Array(W) == Matrix(W)
    if elty <: Real
        @test eltype(Matrix{Float64}(W)) === Float64
        @test eltype(Array{Float64}(W)) === Float64
    end
    @test W*v ≈ F*v
    iFv = F\v
    @test norm(W\v - iFv)/norm(iFv) <= n*cond(F)*ε # Revisit. Condition number is wrong
    @test norm(inv(W)*v - iFv)/norm(iFv) <= n*cond(F)*ε
    @test abs((det(W) - det(F))/det(F)) <= n*cond(F)*ε # Revisit. Condition number is wrong
    iWv = similar(iFv)
    if elty != Int
        iWv = ldiv!(W, copy(v))
        @test iWv ≈ iFv
    end
    @test W'*v ≈ F'*v
    @test transpose(W)*v ≈ transpose(F)*v
    @test (W + W) \ v ≈ (2*Matrix(W)) \ v

    # Diagonal matrix for A (lu(A) fails)
    D = Diagonal(d)
    W = Woodbury(D, U, C, V)
    F = Matrix(W)
    @test F ≈ D + U*C*V
    @test W*v ≈ F*v
    iFv = F\v
    @test norm(W\v - iFv)/norm(iFv) <= n*cond(F)*ε # Revisit. Condition number is wrong

    # Check division with diagonal matrix
    iFvmat = F\Diagonal(v)
    @test norm(W\Diagonal(v) - iFvmat)/norm(iFvmat) <= n*cond(F)*ε # Revisit. Condition number is wrong

    # Factorization for A
    W = Woodbury(lu(T), U, C, V)
    F = Matrix(W)
    @test_broken W*v ≈ F*v   # needs mul!(dst, ::LU, v, ...)
    iFv = F\v
    @test norm(W\v - iFv)/norm(iFv) <= n*cond(F)*ε # Revisit. Condition number is wrong
    @test abs((det(W) - det(F))/det(F)) <= n*cond(F)*ε # Revisit. Condition number is wrong
    iWv = similar(iFv)
    if elty != Int
        iWv = ldiv!(W, copy(v))
        @test iWv ≈ iFv
    end
end

# Sparse U and V
n = 5
d = 1 .+ rand(n)
dl = -rand(n-1)
du = -rand(n-1)
T = Tridiagonal(dl, d, du)
U = sprandn(n,2,0.2)
V = sprandn(2,n,0.2)
C = randn(2,2)
v = randn(n)
ε = eps()

W = Woodbury(T, U, C, V)
F = Matrix(W)
iFv = F\v
@test norm(W\v - iFv)/norm(iFv) <= n*cond(F)*ε

W = Woodbury(lu(T), U, C, V)
F = Matrix(W)
iFv = F\v
@test norm(W\v - iFv)/norm(iFv) <= n*cond(F)*ε

# Display
iob = IOBuffer()
show(iob, W)
str = String(take!(iob))
@test occursin("C:", str)

# Vector-valued U and scalar-valued C
U = rand(n)
C = 2.3
V = rand(1,n)
W = Woodbury(T, U, C, V)
F = Matrix(W)

v = randn(n)
ε = eps()
iFv = F\v
@test norm(W\v - iFv)/norm(iFv) <= n*cond(F)*ε

@test Matrix(2*W) ≈ 2*Matrix(W)
@test Matrix(W*2) ≈ 2*Matrix(W)
R = rand(size(W)...)
@test Matrix(W + R) ≈ Matrix(W) + R
@test Matrix(R + W) ≈ Matrix(W) + R

show(iob, W)

Wc = copy(W)

# Mismatched sizes
@test_throws DimensionMismatch Woodbury(rand(5,5),rand(5,2),rand(2,3),rand(3,5))
@test_throws DimensionMismatch Woodbury(rand(5,5),rand(5,2),rand(3,3),rand(2,5))

# spec builder
n = 10
num_nz = 3
for i in 1:5  # repeat 5 times
    args = [(rand(1:n), rand(1:n), rand()) for j in 1:num_nz]
    r, v, c = WoodburyMatrices.sparse_factors(Float64, n, args...)
    out = r*v*c

    by_hand = zeros(n, n)
    for (i, j, v) in args
        by_hand[i, j] = v
    end

    @test maximum(abs,out - by_hand) == 0.0
end

end  # @testset "Woodbury"