text:
using Statistics,Serialization,StatsPlots
path_to_data = "data//baseline_GER//"
filename_prefix = "figures//"

include("$(path_to_data)covid_par_ini.jl")
worker_results = deserialize(open("$(path_to_data)batchdata.dat"))

datapoint = fld(T,datat)

gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
gdpstd = std(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
nn = size(gdpmean)[1]
pl1 = plot(datat*collect(2:nn),[gdpmean[2:nn], gdpmean[2:nn].-gdpstd[2:nn],gdpmean[2:nn].+gdpstd[2:nn]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black], label = ["GDP" "" ""])
savefig(pl1,"$(filename_prefix)gdpdyn.pdf")

unempmean = mean(hcat((results -> results[:unempltraj]).(worker_results)...), dims=2)
unempstd = std(hcat((results -> results[:unempltraj]).(worker_results)...), dims=2)
pl2 = plot(datat*collect(2:nn),[unempmean[2:nn], unempmean[2:nn].-unempstd[2:nn],unempmean[2:nn].+unempstd[2:nn]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black],label = ["unemp" "" ""])
savefig(pl2,"$(filename_prefix)unempdyn.pdf")

savmean = mean(cat((results -> results[:f_av_savings]).(worker_results)...,dims=3), dims=3)
savstd = std(cat((results -> results[:f_av_savings]).(worker_results)...,dims=3), dims=3)
pl4 = plot(datat*collect(2:nn),[savmean[2:nn,1], savmean[2:nn,1].-savstd[2:nn,1],savmean[2:nn,1].+savstd[2:nn,1]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black], label = ["man" "" ""])
plot!(datat*collect(2:nn),[savmean[2:nn,2], savmean[2:nn,2].-savstd[2:nn,2],savmean[2:nn,2].+savstd[2:nn,2]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:red :black :black], label = ["ser" "" ""])
plot!(datat*collect(2:nn),[savmean[2:nn,3], savmean[2:nn,3].-savstd[2:nn,3],savmean[2:nn,3].+savstd[2:nn,3]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:green :black :black], label = ["food" "" ""])
savefig(pl4,"$(filename_prefix)savdyn.pdf")

RKIR0mean = mean(hcat((results -> results[:RKIR0traj]).(worker_results)...), dims=2)
RKIR0std = std(hcat((results -> results[:RKIR0traj]).(worker_results)...), dims=2)
nn1 = size(RKIR0mean)[1]
pl3 = plot(collect(1:nn1),[RKIR0mean, RKIR0mean.-RKIR0std,RKIR0mean.+RKIR0std], linestyle= [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black],label = ["RKIR0" "" ""])
savefig(pl3,"$(filename_prefix)RKIR0dyn.pdf")

totinfmean = mean(hcat((results -> results[:totinftraj]).(worker_results)...), dims=2)
totinfstd = std(hcat((results -> results[:totinftraj]).(worker_results)...), dims=2)
include("$(path_to_data)emp_traj_100k.jl")
nzer= (Int(ceil(virustime/datat))-1)*datat+2*corlatent
emptotinftraj = vcat(zeros(nzer),emptotinf)
nne = size(emptotinf)[1]
pl5 = plot(collect(1:nn1),[totinfmean, totinfmean.-totinfstd,totinfmean.+totinfstd], linestyle= [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black],label = ["totinf" "" ""])
plot!(emptotinftraj, linecolor = [:green], label="GER", linewidth = [2])
savefig(pl5,"$(filename_prefix)totinfdyn.pdf")

pl5a = plot(collect(1:nne),[totinfmean[nzer+1:nzer+nne], totinfmean[nzer+1:nzer+nne].-totinfstd[nzer+1:nzer+nne],totinfmean[nzer+1:nzer+nne].+totinfstd[nzer+1:nzer+nne]], linestyle= [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black],label = ["totinf" "" ""])
plot!(emptotinftraj[nzer+1:nzer+nne], linecolor = [:green], label="GER", linewidth = [2])
savefig(pl5a,"$(filename_prefix)totinfemp.pdf")

hcapfrac = zeros(datapoint+1)
for th = 1:datapoint+1
    hcapfrac[th] = hcap / (nhh * icufrac)
end
infmean = mean(cat((results -> results[:inftraj]).(worker_results)...,dims=3), dims=3)
infstd = std(cat((results -> results[:inftraj]).(worker_results)...,dims=3), dims=3)
pl6 = plot(datat*collect(2:nn),[infmean[2:nn,1], infmean[2:nn,1].-infstd[2:nn,1],infmean[2:nn,1].+infstd[2:nn,1],hcapfrac[2:nn]], linestyle = [:solid :dot :dot :solid], linewidth = [2 1 1 1], linecolor = [:blue :black :black :black],label = ["y" "" "" "cap"])
plot!(datat*collect(2:nn),[infmean[2:nn,2], infmean[2:nn,2].-infstd[2:nn,2],infmean[2:nn,2].+infstd[2:nn,2]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:red :black :black],label = ["o" "" ""])
plot!(datat*collect(2:nn),[infmean[2:nn,3], infmean[2:nn,3].-infstd[2:nn,3],infmean[2:nn,3].+infstd[2:nn,3]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:green :black :black],label = ["tot" "" ""])
savefig(pl6,"$(filename_prefix)infdyn.pdf")

concountmean = mean(hcat((results -> results[:contact_count_traj]).(worker_results)...), dims=2)
concountstd = std(hcat((results -> results[:contact_count_traj]).(worker_results)...), dims=2)
nnhc = nn1 - size(concountmean)[1]
concountmean = vcat(zeros(nnhc),concountmean)
concountstd = vcat(zeros(nnhc),concountstd)
pl17 = plot(collect(1:nn1),[concountmean, concountmean.-concountstd,concountmean.+concountstd], linestyle= [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black],label = ["totcont" "" ""])
savefig(pl17,"$(filename_prefix)concountdyn.pdf")

totcasmean = mean(cat((results -> results[:totcastraj]).(worker_results)...,dims=3), dims=3)
totcasstd = std(cat((results -> results[:totcastraj]).(worker_results)...,dims=3), dims=3)
pl7 = plot(datat*collect(2:nn),[totcasmean[2:nn,1], totcasmean[2:nn,1].-totcasstd[2:nn,1],totcasmean[2:nn,1].+totcasstd[2:nn,1]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black], label = ["y" "" ""])
plot!(datat*collect(2:nn),[totcasmean[2:nn,2], totcasmean[2:nn,2].-totcasstd[2:nn,2],totcasmean[2:nn,2].+totcasstd[2:nn,2]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:red :black :black],label = ["o" "" ""])
plot!(datat*collect(2:nn),[totcasmean[2:nn,3], totcasmean[2:nn,3].-totcasstd[2:nn,3],totcasmean[2:nn,3].+totcasstd[2:nn,3]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:green :black :black], label = ["tot" "" ""])
savefig(pl7,"$(filename_prefix)casdyn.pdf")

pubh = (datat/(365*gdpmean[2]))*hcat((results -> results[:pubacctraj]).(worker_results)...)
pubmean = mean(pubh, dims=2)
pubstd = std(pubh, dims=2)
pl8 = plot(datat*collect(1:nn),[pubmean, pubmean.-pubstd,pubmean.+pubstd], linestyle= [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black],label = ["pubacc" "" ""])
savefig(pl8,"$(filename_prefix)pubaccdyn.pdf")

inactmean = mean(hcat((results -> results[:inactivefirmstraj]).(worker_results)...), dims=2)
inactstd = std(hcat((results -> results[:inactivefirmstraj]).(worker_results)...), dims=2)
pl9 = plot(datat*collect(1:nn),[inactmean, inactmean.-inactstd,inactmean.+inactstd], linestyle= [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black],label = ["inact firms" "" ""])
savefig(pl9,"$(filename_prefix)inactf.pdf")

gdplossmean = mean(hcat((results -> results[:gdplosstraj]).(worker_results)...), dims=2)
gdplossstd = std(hcat((results -> results[:gdplosstraj]).(worker_results)...), dims=2)
pl10 = plot(datat*collect(1:nn),[gdplossmean, gdplossmean.-gdplossstd,gdplossmean.+gdplossstd], linestyle= [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black],label = ["gdp loss" "" ""])
savefig(pl10,"$(filename_prefix)gdplossdyn.pdf")

nzer=(Int(ceil(virustime/datat))-1)*datat+2*corlatent
empR0smtraj = vcat(zeros(nzer),empR0sm)
RKIR0smmean = mean(hcat((results -> results[:RKIR0smtraj]).(worker_results)...), dims=2)
RKIR0smstd = std(hcat((results -> results[:RKIR0smtraj]).(worker_results)...), dims=2)
nnhc = nn1 - size(RKIR0smmean)[1]
RKIR0smmean = vcat(zeros(nnhc),RKIR0smmean)
RKIR0smstd = vcat(zeros(nnhc),RKIR0smstd)
pl16 = plot(collect(1:nn1),[RKIR0smmean, RKIR0smmean.-RKIR0smstd,RKIR0smmean.+RKIR0smstd], linestyle= [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black],label = ["RKIR0" "" ""])
plot!(empR0smtraj, linecolor = [:green], label="GER", linewidth = [2])
savefig(pl16,"$(filename_prefix)RKIR0smdyn.pdf")

pl16a = plot(collect(1:nne),[RKIR0smmean[nzer+1:nzer+nne], RKIR0smmean[nzer+1:nzer+nne].-RKIR0smstd[nzer+1:nzer+nne],RKIR0smmean[nzer+1:nzer+nne].+RKIR0smstd[nzer+1:nzer+nne]], linestyle= [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black],label = ["RKIR0" "" ""])
plot!(empR0smtraj[nzer+1:nzer+nne], linecolor = [:green], label="GER", linewidth = [2])
savefig(pl16a,"$(filename_prefix)RKIR0smemp.pdf")

bankmean = mean(cat((results -> results[:bankrupttraj]).(worker_results)...,dims=3), dims=3)
bankstd = std(cat((results -> results[:bankrupttraj]).(worker_results)...,dims=3), dims=3)
pl8 = plot(datat*collect(2:nn),[bankmean[2:nn,1], bankmean[2:nn,1].-bankstd[2:nn,1],bankmean[2:nn,1].+bankstd[2:nn,1]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black], label = ["man" "" ""])
plot!(datat*collect(2:nn),[bankmean[2:nn,2], bankmean[2:nn,2].-bankstd[2:nn,2],bankmean[2:nn,2].+bankstd[2:nn,2]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:red :black :black], label = ["ser" "" ""])
plot!(datat*collect(2:nn),[bankmean[2:nn,3], bankmean[2:nn,3].-bankstd[2:nn,3],bankmean[2:nn,3].+bankstd[2:nn,3]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:green :black :black], label = ["food" "" ""])
savefig(pl8,"$(filename_prefix)bankruptdyn.pdf")

unempsecmean = mean(cat((results -> results[:unempsectraj]).(worker_results)...,dims=3), dims=3)
unempsecstd = std(cat((results -> results[:unempsectraj]).(worker_results)...,dims=3), dims=3)
pl22 = plot(datat*collect(2:nn),[unempsecmean[2:nn,1], unempsecmean[2:nn,1].-unempsecstd[2:nn,1],unempsecmean[2:nn,1].+unempsecstd[2:nn,1]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black], label = ["man" "" ""])
plot!(datat*collect(2:nn),[unempsecmean[2:nn,2], unempsecmean[2:nn,2].-unempsecstd[2:nn,2],unempsecmean[2:nn,2].+unempsecstd[2:nn,2]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:red :black :black], label = ["ser" "" ""])
plot!(datat*collect(2:nn),[unempsecmean[2:nn,3], unempsecmean[2:nn,3].-unempsecstd[2:nn,3],unempsecmean[2:nn,3].+unempsecstd[2:nn,3]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:green :black :black], label = ["food" "" ""])
savefig(pl22,"$(filename_prefix)unempsecdyn.pdf")

gdplossh = (results -> results[:togdploss]).(worker_results)
totcash = (results -> results[:totcas]).(worker_results)
pubacch = (results -> results[:pubacctraj]).(worker_results)
pubacc = []
for i = 1:size(pubacch)[1]
    push!(pubacc,pubacch[i][datapoint+1]*datat/(365*gdpmean[2]))
end

pl19 =  boxplot(gdplossh, label = "Av. GDP Loss")
savefig(pl19,"$(filename_prefix)totgdploss.pdf")

pl20 =  boxplot(totcash, label = "Tot. Cas.")
savefig(pl20,"$(filename_prefix)totcas.pdf")

pl21 =  boxplot(pubacc, label = "Pub. Acc.")
savefig(pl21,"$(filename_prefix)pubacc.pdf")

polswcounth = (results -> results[:polswitchcount]).(worker_results)
badpoltimeh = (results -> results[:badpoltime]).(worker_results)

pl23 =  boxplot(polswcounth, label = "Pol. Sw.")
savefig(pl23,"$(filename_prefix)polswitch.pdf")

pl24 =  boxplot(badpoltimeh, label = "Restr.")
savefig(pl24,"$(filename_prefix)badpol.pdf")


text:
using DelimitedFiles
using LinearAlgebra
using SparseArrays

import CPLEX

using CLaP

# Read model from file
cbf = MOI.FileFormats.CBF.Model()
fname = ARGS[1]
fdir = dirname(fname)
finst = basename(fname)[1:end-4]

MOI.read_from_file(cbf, fname)

# Create standard form problem
micp_optimizer = MOI.OptimizerWithAttributes(CPLEX.Optimizer,
    "CPX_PARAM_THREADS" => 1,                       # Single thread
    "CPXPARAM_MIP_Strategy_MIQCPStrat" => 2,        # 1: NL-B&B, 2:OA
)

sf = build_standard_form(micp_optimizer, cbf, bridge_type=Float64)

# Solve root node with aggressive L-a-P cuts
micp = sf.model

# Time limit
# TODO: provide as input
MOI.set(micp, MOI.TimeLimitSec(), 3600.0)  # 1-hour time limit

MOI.optimize!(micp)

@show MOI.get(micp, MOI.TerminationStatus())    # Termination status
@show MOI.get(micp, MOI.ObjectiveBound())       # Lower bound

# Export solution
# TODO: check for solution status first
x_micp = MOI.get(micp, MOI.VariablePrimal(), sf.var_indices)

writedlm(joinpath(@__DIR__, "res/$finst.sol"), x_micp)

ext: jl

max_stars_repo_name: mtanneau/CLaP

max_issues_repo_name: mtanneau/CLaP

max_forks_repo_name: mtanneau/CLaP


text:
module ConstructionTools

export Subdomains,
add_subdomain,
build_domain,
build_pc_domain,
line_defect,
site_defect,
add_halfspace_waveguide,
add_halfspace_waveguides,
add_planar_waveguide,
add_planar_waveguides,
add_pc_waveguide,
add_pc_waveguides

include("ConstructionToolsBase.jl")
using .ConstructionToolsBase

include("PhotonicCrystal.jl")
using .PhotonicCrystal

include("Waveguides.jl")
using .Waveguides

end #module


text:
function movingaverage(::Type{T} = Float64; lim::Integer = -1) where T<:Real
	buffer = Vector{T}(0)
	if lim == -1
		# unlimited buffer
		return (y::T) -> begin
			push!(buffer, y)
			return mean(buffer)
		end
	else
		# limited size buffer
		return (y) -> begin
			push!(buffer, y)
			if length(buffer) > lim shift!(buffer) end
			return mean(buffer)
		end
	end
end

test = movingaverage()
@show test(1.0) # mean([1])
@show test(2.0) # mean([1, 2])
@show test(3.0) # mean([1, 2, 3])

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# FUNCTION recover_primal
"""Computes the Ridge regressor
INPUT
  ℓ           - LossFunction to use
  Y           - Vector of observed responses
  Z           - Matrix of observed features
  γ           - Regularization parameter
OUTPUT
  w           - Optimal regressor"""
function recover_primal(ℓ::Regression, Y, Z, γ)
  CM = Matrix(I, size(Z,2), size(Z,2))/γ + Z'*Z      # The capacitance matrix
  α = -Y + Z*(CM\(Z'*Y))            # Matrix Inversion Lemma
  return -γ*Z'*α                    # Regressor
end

using LIBLINEAR

function recover_primal(ℓ::Classification, Y, Z, γ)
  solverNumber = LibLinearSolver(ℓ)
  if isa(ℓ, SubsetSelection.Classification)
    model = LIBLINEAR.linear_train(Y, Z'; verbose=false, C=γ, solver_type=Cint(solverNumber))
    return Y[1]*model.w
  # else
  #   model = linear_train(Y, Z'; verbose=false, C=γ, solver_type=Cint(solverNumber), eps = ℓ.ε)
  #   return model.w
  end
end

function LibLinearSolver(ℓ::L1SVR)
  return 13
end
function LibLinearSolver(ℓ::L2SVR)
  return 12
end

function LibLinearSolver(ℓ::LogReg)
  return 7
end
function LibLinearSolver(ℓ::L1SVM)
  return 3
end
function LibLinearSolver(ℓ::L2SVM)
  return 2
end


text:
__precompile__(true)

module Stochastic

using Distributions
using Gadfly

import Distributions: rand
import Gadfly: plot

export

  # random process types
  RandomProcess,
  BrownianMotion,

  # methods
  params

  # source files
  include("Plots.jl")
  include("BrownianMotion.jl")

end


text:
using Test
using ScottishTaxBenefitModel.HistoricBenefits
using ScottishTaxBenefitModel
using .Definitions, .Utils, .Randoms
using .ExampleHelpers

## NOTE this test has the 2019 OBR data and 2019Q4 as a target jammed on - will need
## changing with update versions

@testset "historic tests" begin

    be_99 = 66.75
    pe_99 = 66.75
    @test benefit_ratio(
        1999,
        be_99,
        bereavement_allowance_or_widowed_parents_allowance_or_bereavement
    ) ≈ 1

    @test benefit_ratio(
        2020,
        be_99,
        bereavement_allowance_or_widowed_parents_allowance_or_bereavement
    ) ≈ be_99/121.95

    @test benefit_ratio(
        1999,
        pe_99,
        state_pension
    ) ≈ 1

    d = Date( 2019, 1, 1)
    @test DLA_RECEIPTS[nearest( d, DLA_RECEIPTS ),:Scotland] ≈ 182_154
    @test PIP_RECEIPTS[nearest( d, PIP_RECEIPTS ),:Scotland] ≈ 220_043

    @test should_switch_dla_to_pip( "000000000000", 2016, 1 )
    @test ! should_switch_dla_to_pip("99999999999", 2016, 1 )

    names = ExampleHouseholdGetter.initialise()
    scot = ExampleHouseholdGetter.get_household( "mel_c2_scot" ) # scots are a married couple
    head = scot.people[SCOT_HEAD]
    spouse = scot.people[SCOT_SPOUSE]

    disable_seriously!( head ) # don't really need this ...
    head.onerand = String(fill( '9', 100 )) # 0.99999.. switch will never be made
    @test randchunk( head.onerand, 3, 3 ) ≈ 0.999
    head.dla_self_care_type = high
    head.dla_mobility_type = high
    switch_dla_to_pip!( head, 2016, 3 )
    @test head.dla_self_care_type == high
    @test head.dla_mobility_type == high
    @test head.pip_mobility_type == no_pip
    @test head.pip_daily_living_type == no_pip

    head.onerand = String(fill( '0', 100 )) # switch always made
    @test randchunk( head.onerand, 3, 3 ) == 0
    head.dla_self_care_type = high
    head.dla_mobility_type = high
    head.pip_mobility_type = no_pip
    head.pip_daily_living_type = no_pip
    switch_dla_to_pip!( head, 2016, 3 )
    @test head.dla_self_care_type == missing_lmh
    @test head.dla_mobility_type == missing_lmh
    @test head.pip_mobility_type == enhanced_pip
    @test head.pip_daily_living_type == enhanced_pip

end # testset

@testset "Run on actual Data" begin
    # just see if switcher roughly works on actual data
    nhhs,npeople = init_data( reset=true )
    nyears = 6
    pips_dlas = DataFrame(
        year = zeros( Int, nyears ),
        orig_dla_self_care = zeros( Int, nyears ),
        new_dla_self_care = zeros( Int, nyears ),
        orig_dla_mobility = zeros( Int, nyears ),
        new_dla_mobility = zeros( Int, nyears ),
        orig_pip_daily_living = zeros( Int, nyears ),
        new_pip_daily_living = zeros( Int, nyears ),
        orig_pip_mobility = zeros( Int, nyears ),
        new_pip_mobility = zeros( Int, nyears ),
    )
    for hno in 1:nhhs
        hh = get_household(hno)
        year = hh.interview_year
        month = hh.interview_month
        yp = year - 2014
        # println( "on year $year month $month yp=$yp hno=$hno")
        for (pid,pers) in hh.people # just unweighted counts
            dr = pips_dlas[yp,:]
            pips_dlas[yp,:year] = year
            if pers.dla_self_care_type != missing_lmh
                dr.orig_dla_self_care += 1
            end
            if pers.dla_mobility_type != missing_lmh
                dr.orig_dla_mobility += 1
            end
            if pers.pip_daily_living_type != no_pip
                dr.orig_pip_daily_living += 1
            end
            if pers.pip_mobility_type != no_pip
                dr.orig_pip_mobility += 1
            end
            println( pers.onerand )
            switch_dla_to_pip!( pers, year, month )
            # ss = should_switch_dla_to_pip( pers.onerand, year, month )
            # println( "should switch $ss ")
            if pers.dla_self_care_type != missing_lmh
                dr.new_dla_self_care += 1
            end
            if pers.dla_mobility_type != missing_lmh
                dr.new_dla_mobility += 1
            end
            if pers.pip_daily_living_type != no_pip
                dr.new_pip_daily_living += 1
            end
            if pers.pip_mobility_type != no_pip
                dr.new_pip_mobility += 1
            end

        end # people
    end # hhs
    println( pips_dlas )
end # testset


text:
function toh01()
#
#  generalization of K.C. Toh and L. N. Trefethen
#
    k = [(-10:9)...]
    z = 2.0 .* (k .+ 0.5) ./ 19 + im .* sin.((2*pi) .* (k .+ 0.5) ./ 19)
    p = poly(z)
    z = [z ones(20)]

    p, PolyZeros(z)
end


text:
using ValidatedNumerics, ValidatedNumerics.RootFinding
using ForwardDiff

using FactCheck

const Dual = ForwardDiff.Dual

facts("Promotion between Dual and Interval") do
    @fact promote(ForwardDiff.Dual(2, 1), Interval(1, 2)) -->
        (Dual(2..2, 1..1), Dual(1..2, 0..0))

    @fact promote(Interval(2, 3), Dual(2, 1)) --> (Dual(Interval(2, 3), Interval(0)),
            Dual(Interval(2.0), Interval(1.0)))

end


text:
using GPMLj, Test, Random, PyCall

Random.seed!(123)

@testset "Priors" begin

    @testset "Gaussian" begin
        prior = gpflow.priors.Gaussian(0.0, 1.0)
        typeof(prior)<:gpflow.Prior
        @test typeof(prior)<:gpflow.priors.Gaussian
        @test typeof(prior.o)<:PyObject
        temp = prior.o
        @test instantiate!(prior) == temp
    end

    @testset "LogNormal" begin
        prior = gpflow.priors.LogNormal(0.0, 1.0)
        typeof(prior)<:gpflow.Prior
        @test typeof(prior)<:gpflow.priors.LogNormal
        @test typeof(prior.o)<:PyObject
        temp = prior.o
        @test instantiate!(prior) == temp
    end

    @testset "Gamma" begin
        prior = gpflow.priors.Gamma((1,), 1.0)
        typeof(prior)<:gpflow.Prior
        @test typeof(prior)<:gpflow.priors.Gamma
        @test typeof(prior.o)<:PyObject
        temp = prior.o
        @test instantiate!(prior) == temp
    end

    @testset "Laplace" begin
        prior = gpflow.priors.Laplace(0.0, 1.0)
        typeof(prior)<:gpflow.Prior
        @test typeof(prior)<:gpflow.priors.Laplace
        @test typeof(prior.o)<:PyObject
        temp = prior.o
        @test instantiate!(prior) == temp
    end

    @testset "Uniform" begin
        prior = gpflow.priors.Uniform()
        typeof(prior)<:gpflow.Prior
        @test typeof(prior)<:gpflow.priors.Uniform
        @test typeof(prior.o)<:PyObject
        temp = prior.o
        @test instantiate!(prior) == temp
    end

end


text:
# ---
# title: 350. Intersection of Two Arrays II
# id: problem350
# author: zhwang
# date: 2022-01-19
# difficulty: Easy
# categories: Hash Table, Two Pointers, Binary Search, Sort
# link: <https://leetcode.com/problems/intersection-of-two-arrays-ii/description/>
# hidden: true
# ---
#
# Given two arrays, write a function to compute their intersection.
#
# **Example 1:**
#
#
#
#     Input: nums1 = [1,2,2,1], nums2 = [2,2]
#     Output: [2,2]
#
#
# **Example 2:**
#
#
#
#     Input: nums1 = [4,9,5], nums2 = [9,4,9,8,4]
#     Output: [4,9]
#
# **Note:**
#
#   * Each element in the result should appear as many times as it shows in both arrays.
#   * The result can be in any order.
#
# **Follow up:**
#
#   * What if the given array is already sorted? How would you optimize your algorithm?
#   * What if _nums1_ 's size is small compared to _nums2_ 's size? Which algorithm is better?
#   * What if elements of _nums2_ are stored on disk, and the memory is limited such that you cannot load all elements into the memory at once?
#
#
## @lc code=start
using LeetCode

function intersection_of_arrays(nums1::Vector{Int}, nums2::Vector{Int})::Vector{Int}
    nums1, nums2 = sort(nums1), sort(nums2)
    n1, n2 = length(nums1), length(nums2)
    i1, i2, res = 1, 1, Int[]
    while i1 <= n1 && i2 <= n2
        if nums1[i1] > nums2[i2]
            i2 += 1
        elseif nums1[i1] < nums2[i2]
            i1 += 1
        else
            push!(res, nums1[i1])
            i1 += 1
            i2 += 1
        end
    end
    return res
end

## @lc code=end


text:
"""
Function `condition` is similar to `assume`, however it does not generate new values for `vname`
from `distr`. Instead, it condition on an existing value of `vname`. This existing value can be
mutated externally by some inference method like SMC2, HMC, PMMH etc.
"""
immutable Conditional
  vals   :: Array{Distributions.VariateForm, 1}
  distrs :: Array{Distributions.Distribution,1}
  pos  :: TArray{Int,1}  # N.B. TArray is only meaningful for ParticleSamplers, since each particle needs a counter for its passed parameters.
end

# Condition on a variable during inference, e.g. in SMC2, we sample x | theta in the inner SMC.
#  A new value is sampled from the prior, if no values for theta is found in `Sampler.conditionals`.
#  Otherwise re-use or replay the value in `Sampler.conditionals`.

# theta = condition( spl :: Sampler, distr :: Distribution )

# Utility functions

# flatten( c :: Conditional )
# logpdf (c :: Conditional ) is directly added to `Trace.logweight` in each call of `condition`

#Params does not contain a counter, since it is Xparticle specific
#the length of the arrays should be equal to the number of elements stored in them
type Params
  values::Array{Any} # values of the random variables in a run
  dists::Array{Any} # the prior on the parameters
end
Params() = Params([],[])

# copies the entire structure
## TODO: #is shallow copy enough here?
clone(ps::Params) = Params(deepcopy(ps.values), deepcopy(ps.dists))

# TODO: rename prior to logpdf
# function prior(params :: Params) #prior logdensity for the current values
#   result = 0
#   for i=1:length(params.values)
#     result += logpdf(params.dists[i], params.values[i], true)
#   end
#   return result
# end
prior(params :: Params) = reduce(+,map((d,v)->logpdf(d,v,true), params.dists, params.values))
logpdf(params :: Params) = prior(params)

#Flatten should return a flat array with all the elements and a function that can rebuild the
#original data structure from such an array.
#It takes a distribution as the second argument to compute the constraints.
#The returned values are transformed to cover the whole real line.
function flatten(ps :: Params)
  a = map(flatten, ps.values, ps.dists)
  rebuilders = map(x -> x[2], a)
  vectors = map(x -> x[1], a)
  sizes = map(length, vectors)
  flat = Float64[]
  foldl(append!, flat, vectors)
  function reb(flat)
    consumed = 0
    objects = Any[]
    while(consumed < length(flat))
      which = length(objects) + 1
      n = sizes[which]
      push!(objects, rebuilders[which](flat[(consumed + 1) : (consumed + n)]))
      consumed += n
    end
    return Params(objects, deepcopy(ps.dists))
  end
  return (flat, reb)
end


text:
module Mechanical1D
using HoJMTK

@parameters t
const D = Differential(t)

include("mech1d-core.jl")
include("mech1d-sources.jl")
include("mech1d-components.jl")

export
    Fixed,
    connect_series,
    connect_parallel,
    connect_star,
    ConstantVelocity,
    Force,
    Damper,
    Mass,
    Spring
end


text:
function add_variables!(m::JuMP.AbstractModel, ps::AbstractVector{REoptInputs{T}}) where T <: AbstractScenario

	dvs_idx_on_techs = String[
		"dvSize",
		"dvPurchaseSize",
	]
	dvs_idx_on_techs_time_steps = String[
        "dvCurtail",
		"dvRatedProduction",
	]
	dvs_idx_on_storagetypes = String[
		"dvStoragePower",
		"dvStorageEnergy",
	]
	dvs_idx_on_storagetypes_time_steps = String[
		"dvDischargeFromStorage"
	]
	for p in ps
		_n = string("_", p.s.site.node)
		for dv in dvs_idx_on_techs
			x = dv*_n
			m[Symbol(x)] = @variable(m, [p.techs.all], base_name=x, lower_bound=0)
		end

		for dv in dvs_idx_on_techs_time_steps
			x = dv*_n
			m[Symbol(x)] = @variable(m, [p.techs.all, p.time_steps], base_name=x, lower_bound=0)
		end

		for dv in dvs_idx_on_storagetypes
			x = dv*_n
			m[Symbol(x)] = @variable(m, [p.s.storage.types.elec], base_name=x, lower_bound=0)
		end

		for dv in dvs_idx_on_storagetypes_time_steps
			x = dv*_n
			m[Symbol(x)] = @variable(m, [p.s.storage.types.all, p.time_steps], base_name=x, lower_bound=0)
		end

		dv = "dvGridToStorage"*_n
		m[Symbol(dv)] = @variable(m, [p.s.storage.types.elec, p.time_steps], base_name=dv, lower_bound=0)

		dv = "dvGridPurchase"*_n
		m[Symbol(dv)] = @variable(m, [p.time_steps], base_name=dv, lower_bound=0)

		dv = "dvPeakDemandTOU"*_n
		m[Symbol(dv)] = @variable(m, [p.ratchets, 1], base_name=dv, lower_bound=0)

		dv = "dvPeakDemandMonth"*_n
		m[Symbol(dv)] = @variable(m, [p.months, 1], base_name=dv, lower_bound=0)

		dv = "dvProductionToStorage"*_n
		m[Symbol(dv)] = @variable(m, [p.s.storage.types.all, p.techs.all, p.time_steps], base_name=dv, lower_bound=0)

		dv = "dvStoredEnergy"*_n
		m[Symbol(dv)] = @variable(m, [p.s.storage.types.all, 0:p.time_steps[end]], base_name=dv, lower_bound=0)

		dv = "MinChargeAdder"*_n
		m[Symbol(dv)] = @variable(m, base_name=dv, lower_bound=0)

        if !isempty(p.s.electric_tariff.export_bins)
            dv = "dvProductionToGrid"*_n
            m[Symbol(dv)] = @variable(m, [p.techs.elec, p.s.electric_tariff.export_bins, p.time_steps], base_name=dv, lower_bound=0)
        end

		ex_name = "TotalTechCapCosts"*_n
		m[Symbol(ex_name)] = @expression(m, p.third_party_factor *
			sum( p.cap_cost_slope[t] * m[Symbol("dvPurchaseSize"*_n)][t] for t in p.techs.all )
		)

		ex_name = "TotalStorageCapCosts"*_n
		m[Symbol(ex_name)] = @expression(m, p.third_party_factor *
			sum(p.s.storage.attr[b].net_present_cost_per_kw * m[Symbol("dvStoragePower"*_n)][b] for b in p.s.storage.types.elec)
			+ sum(p.s.storage.attr[b].net_present_cost_per_kwh * m[Symbol("dvStorageEnergy"*_n)][b] for b in p.s.storage.types.all)
		)

		ex_name = "TotalPerUnitSizeOMCosts"*_n
		m[Symbol(ex_name)] = @expression(m, p.third_party_factor * p.pwf_om *
			sum( p.om_cost_per_kw[t] * m[Symbol("dvSize"*_n)][t] for t in p.techs.all )
		)

        ex_name = "TotalPerUnitProdOMCosts"*_n
		m[Symbol(ex_name)] = 0

		add_elec_utility_expressions(m, p; _n=_n)

		#################################  Objective Function   ########################################
		m[Symbol("Costs"*_n)] = @expression(m,
			# Capital Costs
			m[Symbol("TotalTechCapCosts"*_n)] + m[Symbol("TotalStorageCapCosts"*_n)] +

			## Fixed O&M, tax deductible for owner
			m[Symbol("TotalPerUnitSizeOMCosts"*_n)] * (1 - p.s.financial.owner_tax_pct) +

			# Utility Bill, tax deductible for offtaker, including export benefit
			m[Symbol("TotalElecBill"*_n)] * (1 - p.s.financial.offtaker_tax_pct)
		);
    end
    add_bounds(m, ps)
end

"""
add non-negative bounds to decision variables
"""
function add_bounds(m::JuMP.AbstractModel, ps::AbstractVector{REoptInputs{T}}) where T <: AbstractScenario

	dvs_idx_on_techs = String[
		"dvSize",
		"dvPurchaseSize",
	]
	dvs_idx_on_techs_time_steps = String[
        "dvCurtail",
		"dvRatedProduction",
	]
	dvs_idx_on_storagetypes = String[
		"dvStoragePower",
		"dvStorageEnergy",
	]
	dvs_idx_on_storagetypes_time_steps = String[
		"dvDischargeFromStorage",
		"dvGridToStorage",
	]
	for p in ps
        _n = string("_", p.s.site.node)

		for dv in dvs_idx_on_techs
			x = dv*_n
			@constraint(m, [tech in p.techs.all], -m[Symbol(x)][tech] ≤ 0 )
		end

		for dv in dvs_idx_on_techs_time_steps
			x = dv*_n
            @constraint(m, [tech in p.techs.all, ts in p.time_steps],
                -m[Symbol(x)][tech, ts] ≤ 0
            )
		end

		for dv in dvs_idx_on_storagetypes
			x = dv*_n
            @constraint(m, [b in p.s.storage.types.elec],
                -m[Symbol(x)][b] ≤ 0
            )
		end

		for dv in dvs_idx_on_storagetypes_time_steps
			x = dv*_n
            @constraint(m, [b in p.s.storage.types.elec, ts in p.time_steps],
                -m[Symbol(x)][b, ts] ≤ 0
            )
		end

		dv = "dvGridPurchase"*_n
		@constraint(m, [ts in p.time_steps], -m[Symbol(dv)][ts] ≤ 0)

		dv = "dvPeakDemandTOU"*_n
		@constraint(m, [r in p.ratchets], -m[Symbol(dv)][r, 1] ≤ 0)

		dv = "dvPeakDemandMonth"*_n
		@constraint(m, [mth in p.months], -m[Symbol(dv)][mth, 1] ≤ 0)

		dv = "dvProductionToStorage"*_n
        @constraint(m, [b in p.s.storage.types.elec, tech in p.techs.all, ts in p.time_steps],
            -m[Symbol(dv)][b, tech, ts] ≤ 0
        )

		dv = "dvStoredEnergy"*_n
        @constraint(m, [b in p.s.storage.types.elec, ts in 0:p.time_steps[end]],
            -m[Symbol(dv)][b, ts] ≤ 0
        )

		dv = "MinChargeAdder"*_n
		@constraint(m, -m[Symbol(dv)] ≤ 0)
    end
end

function build_reopt!(m::JuMP.AbstractModel, ps::AbstractVector{REoptInputs{T}}) where T <: AbstractScenario
    add_variables!(m, ps)
    @warn "Outages are not currently modeled in multinode mode."
    @warn "Diesel generators are not currently modeled in multinode mode."
    for p in ps
        _n = string("_", p.s.site.node)

        for b in p.s.storage.types.all
            if p.s.storage.attr[b].max_kw == 0 || p.s.storage.attr[b].max_kwh == 0
                @constraint(m, [ts in p.time_steps], m[Symbol("dvStoredEnergy"*_n)][b, ts] == 0)
                @constraint(m, m[Symbol("dvStorageEnergy"*_n)][b] == 0)
                @constraint(m, m[Symbol("dvStoragePower"*_n)][b] == 0)
                @constraint(m, [t in p.techs.elec, ts in p.time_steps_with_grid],
                            m[Symbol("dvProductionToStorage"*_n)][b, t, ts] == 0)
                @constraint(m, [ts in p.time_steps], m[Symbol("dvDischargeFromStorage"*_n)][b, ts] == 0)
                @constraint(m, [ts in p.time_steps], m[Symbol("dvGridToStorage"*_n)][b, ts] == 0)
            else
                add_storage_size_constraints(m, p, b; _n=_n)
                add_general_storage_dispatch_constraints(m, p, b; _n=_n)
				if b in p.s.storage.types.elec
					add_elec_storage_dispatch_constraints(m, p, b; _n=_n)
				elseif b in p.s.storage.types.hot
					add_hot_thermal_storage_dispatch_constraints(m, p, b; _n=_n)
				elseif b in p.s.storage.types.cold
					add_cold_thermal_storage_dispatch_constraints(m, p, b; _n=_n)
				end
            end
        end

        if any(max_kw->max_kw > 0, (p.s.storage.attr[b].max_kw for b in p.s.storage.types.elec))
            add_storage_sum_constraints(m, p; _n=_n)
        end

        add_production_constraints(m, p; _n=_n)

        if !isempty(p.techs.all)
            add_tech_size_constraints(m, p; _n=_n)
            if !isempty(p.techs.no_curtail)
                add_no_curtail_constraints(m, p; _n=_n)
            end
        end

        add_elec_load_balance_constraints(m, p; _n=_n)

        if !isempty(p.s.electric_tariff.export_bins)
            add_export_constraints(m, p; _n=_n)
        end

        if !isempty(p.s.electric_tariff.monthly_demand_rates)
            add_monthly_peak_constraint(m, p; _n=_n)
        end

        if !isempty(p.s.electric_tariff.tou_demand_ratchet_time_steps)
            add_tou_peak_constraint(m, p; _n=_n)
        end

		if !(p.s.electric_utility.allow_simultaneous_export_import) & !isempty(p.s.electric_tariff.export_bins)
			add_simultaneous_export_import_constraint(m, p; _n=_n)
		end

        if p.s.electric_tariff.demand_lookback_percent > 0
            add_demand_lookback_constraints(m, p; _n=_n)
        end

        if !isempty(p.s.electric_tariff.coincpeak_periods)
            add_coincident_peak_charge_constraints(m, p; _n=_n)
        end

    end
end

function add_objective!(m::JuMP.AbstractModel, ps::AbstractVector{REoptInputs{T}}) where T <: AbstractScenario
	if !(any(p.s.settings.add_soc_incentive for p in ps))
		@objective(m, Min, sum(m[Symbol(string("Costs_", p.s.site.node))] for p in ps))
	else # Keep SOC high
		@objective(m, Min, sum(m[Symbol(string("Costs_", p.s.site.node))] for p in ps)
        - sum(sum(sum(m[Symbol(string("dvStoredEnergy_", p.s.site.node))][b, ts]
            for ts in p.time_steps) for b in p.s.storage.types.elec) for p in ps) / (8760. / ps[1].hours_per_time_step))
	end  # TODO need to handle different hours_per_time_step?
	nothing
end

function run_reopt(m::JuMP.AbstractModel, ps::AbstractVector{REoptInputs{T}}) where T <: AbstractScenario

	build_reopt!(m, ps)

	add_objective!(m, ps)

	@info "Model built. Optimizing..."
	tstart = time()
	optimize!(m)
	opt_time = round(time() - tstart, digits=3)
	if termination_status(m) == MOI.TIME_LIMIT
		status = "timed-out"
    elseif termination_status(m) == MOI.OPTIMAL
        status = "optimal"
    else
        status = "not optimal"
        @warn "REopt solved with " termination_status(m), ", returning the model."
        return m
	end
	@info "REopt solved with " termination_status(m)
	@info "Solving took $(opt_time) seconds."

	tstart = time()
	results = reopt_results(m, ps)
	time_elapsed = time() - tstart
	@info "Results processing took $(round(time_elapsed, digits=3)) seconds."
	results["status"] = status
	results["solver_seconds"] = opt_time
	return results
end

function reopt_results(m::JuMP.AbstractModel, ps::AbstractVector{REoptInputs{T}}) where T <: AbstractScenario
	# TODO address Warning: The addition operator has been used on JuMP expressions a large number of times.
	results = Dict{Union{Int, String}, Any}()
	for p in ps
		results[p.s.site.node] = reopt_results(m, p; _n=string("_", p.s.site.node))
	end
	return results
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
export reconstruction

include("RecoParameters.jl")
include("DirectReconstruction.jl")
include("IterativeReconstruction.jl")

"""
    reconstruction(acqData::AcquisitionData, recoParams::Dict)

Performs image reconstruction of an AcquisitionData object.
Parameters are specified in a dictionary.

Reconstruction types are specified by the symbol `:reco`.
Valid reconstruction names are:
* :direct - direct Fourier reconstruction
* :standard           - iterative reconstruction for all contrasts, coils & slices independently
* :multiEcho          - iterative joint reconstruction of all echo images
* :multiCoil          - SENSE-type iterative reconstruction
* :multiCoilMultiEcho - SENSE-type iterative reconstruction of all echo images
"""
function reconstruction(acqData::AcquisitionData, recoParams::Dict)
  # check dimensionality of encoding
  encodingDims = dims(trajectory(acqData))
  if encodingDims==3 && numSlices(acqData)>1
    @error "reconstruction of multiple 3d-encoded volumina is not yet supported"
  end

  # load reconstruction parameters
  recoParams = merge(defaultRecoParams(), recoParams)

  # iterative reco
  reconSize, weights, sparseTrafo, reg, normalize, encOps, solvername, senseMaps = setupIterativeReco(acqData, recoParams)
  if recoParams[:reco] == "standard"
    return reconstruction_simple(acqData, reconSize[1:encodingDims], reg, sparseTrafo, weights, solvername, normalize, encOps, recoParams)
  elseif recoParams[:reco] == "multiEcho"
    return reconstruction_multiEcho(acqData, reconSize[1:encodingDims], reg, sparseTrafo, weights, solvername, normalize, encOps, recoParams)
  elseif recoParams[:reco] == "multiCoil"
    return reconstruction_multiCoil(acqData, reconSize[1:encodingDims], reg, sparseTrafo, weights, solvername, senseMaps, normalize, encOps, recoParams)
  elseif recoParams[:reco] == "multiCoilMultiEcho"
    return reconstruction_multiCoilMultiEcho(acqData, reconSize[1:encodingDims], reg, sparseTrafo, weights, solvername, senseMaps, normalize, encOps, recoParams)
  end

  # direct reco
  if recoParams[:reco] != "direct"
    @error "reco modell $(recoParams[:reco]) not found"
  end
  reconSize, weights, cmap = setupDirectReco(acqData, recoParams)
  return reconstruction_direct(acqData, reconSize[1:encodingDims], weights, cmap)
end

"""
    reconstruction(acqData::AcquisitionData, recoParams::Dict,filename::String; force=false)

performs the same image reconstrucion as `reconstruction(acqData::AcquisitionData, recoParams::Dict)`
and saves the image in a file with name `filename`.
If `force=false`, the reconstructed image is loaded from the the file `filename` if the latter is
present.
"""
function reconstruction(acqData::AcquisitionData, recoParams::Dict, filename::String;
                        force=false)
  if !force && isfile(filename)
    return recoImage( RecoFileIBI(filename) )
  else
    I = reconstruction(acqData, recoParams)
    saveasRecoFile(filename, I, recoParams)
    return I
  end
end


text:
export softmax_regression

import ..SoftmaxRegression

using mlpack._Internal.io

import mlpack_jll
const softmax_regressionLibrary = mlpack_jll.libmlpack_julia_softmax_regression

# Call the C binding of the mlpack softmax_regression binding.
function softmax_regression_mlpackMain()
  success = ccall((:softmax_regression, softmax_regressionLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module softmax_regression_internal
  import ..softmax_regressionLibrary

import ...SoftmaxRegression

# Get the value of a model pointer parameter of type SoftmaxRegression.
function IOGetParamSoftmaxRegression(paramName::String)::SoftmaxRegression
  SoftmaxRegression(ccall((:IO_GetParamSoftmaxRegressionPtr, softmax_regressionLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type SoftmaxRegression.
function IOSetParamSoftmaxRegression(paramName::String, model::SoftmaxRegression)
  ccall((:IO_SetParamSoftmaxRegressionPtr, softmax_regressionLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeSoftmaxRegression(stream::IO, model::SoftmaxRegression)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeSoftmaxRegressionPtr, softmax_regressionLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeSoftmaxRegression(stream::IO)::SoftmaxRegression
  buffer = read(stream)
  SoftmaxRegression(ccall((:DeserializeSoftmaxRegressionPtr, softmax_regressionLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    softmax_regression(; [input_model, labels, lambda, max_iterations, no_intercept, number_of_classes, test, test_labels, training, verbose])

This program performs softmax regression, a generalization of logistic
regression to the multiclass case, and has support for L2 regularization.  The
program is able to train a model, load  an existing model, and give predictions
(and optionally their accuracy) for test data.

Training a softmax regression model is done by giving a file of training points
with the `training` parameter and their corresponding labels with the `labels`
parameter. The number of classes can be manually specified with the
`number_of_classes` parameter, and the maximum number of iterations of the
L-BFGS optimizer can be specified with the `max_iterations` parameter.  The L2
regularization constant can be specified with the `lambda` parameter and if an
intercept term is not desired in the model, the `no_intercept` parameter can be
specified.

The trained model can be saved with the `output_model` output parameter. If
training is not desired, but only testing is, a model can be loaded with the
`input_model` parameter.  At the current time, a loaded model cannot be trained
further, so specifying both `input_model` and `training` is not allowed.

The program is also able to evaluate a model on test data.  A test dataset can
be specified with the `test` parameter. Class predictions can be saved with the
`predictions` output parameter.  If labels are specified for the test data with
the `test_labels` parameter, then the program will print the accuracy of the
predictions on the given test set and its corresponding labels.

For example, to train a softmax regression model on the data `dataset` with
labels `labels` with a maximum of 1000 iterations for training, saving the
trained model to `sr_model`, the following command can be used:

```julia
julia> using CSV
julia> dataset = CSV.read("dataset.csv")
julia> labels = CSV.read("labels.csv"; type=Int)
julia> sr_model, _ = softmax_regression(labels=labels,
            training=dataset)
```

Then, to use `sr_model` to classify the test points in `test_points`, saving the
output predictions to `predictions`, the following command can be used:

```julia
julia> using CSV
julia> test_points = CSV.read("test_points.csv")
julia> _, predictions = softmax_regression(input_model=sr_model,
            test=test_points)
```

# Arguments

 - `input_model::SoftmaxRegression`: File containing existing model
      (parameters).
 - `labels::Array{Int, 1}`: A matrix containing labels (0 or 1) for the
      points in the training set (y). The labels must order as a row.
 - `lambda::Float64`: L2-regularization constant  Default value `0.0001`.

 - `max_iterations::Int`: Maximum number of iterations before termination.
       Default value `400`.

 - `no_intercept::Bool`: Do not add the intercept term to the model.
      Default value `false`.

 - `number_of_classes::Int`: Number of classes for classification; if
      unspecified (or 0), the number of classes found in the labels will be
      used.  Default value `0`.

 - `test::Array{Float64, 2}`: Matrix containing test dataset.
 - `test_labels::Array{Int, 1}`: Matrix containing test labels.
 - `training::Array{Float64, 2}`: A matrix containing the training set
      (the matrix of predictors, X).
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output_model::SoftmaxRegression`: File to save trained softmax
      regression model to.
 - `predictions::Array{Int, 1}`: Matrix to save predictions for test
      dataset into.

"""
function softmax_regression(;
                            input_model::Union{SoftmaxRegression, Missing} = missing,
                            labels = missing,
                            lambda::Union{Float64, Missing} = missing,
                            max_iterations::Union{Int, Missing} = missing,
                            no_intercept::Union{Bool, Missing} = missing,
                            number_of_classes::Union{Int, Missing} = missing,
                            test = missing,
                            test_labels = missing,
                            training = missing,
                            verbose::Union{Bool, Missing} = missing,
                            points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, softmax_regressionLibrary), Nothing, ());

  IORestoreSettings("Softmax Regression")

  # Process each input argument before calling mlpackMain().
  if !ismissing(input_model)
    softmax_regression_internal.IOSetParamSoftmaxRegression("input_model", convert(SoftmaxRegression, input_model))
  end
  if !ismissing(labels)
    IOSetParamURow("labels", labels)
  end
  if !ismissing(lambda)
    IOSetParam("lambda", convert(Float64, lambda))
  end
  if !ismissing(max_iterations)
    IOSetParam("max_iterations", convert(Int, max_iterations))
  end
  if !ismissing(no_intercept)
    IOSetParam("no_intercept", convert(Bool, no_intercept))
  end
  if !ismissing(number_of_classes)
    IOSetParam("number_of_classes", convert(Int, number_of_classes))
  end
  if !ismissing(test)
    IOSetParamMat("test", test, points_are_rows)
  end
  if !ismissing(test_labels)
    IOSetParamURow("test_labels", test_labels)
  end
  if !ismissing(training)
    IOSetParamMat("training", training, points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output_model")
  IOSetPassed("predictions")
  # Call the program.
  softmax_regression_mlpackMain()

  return softmax_regression_internal.IOGetParamSoftmaxRegression("output_model"),
         IOGetParamURow("predictions")
end


text:
# This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license

using SpecialFunctions
using Base.Test

const SF = SpecialFunctions

@testset "error functions" begin
    @test SF.erf(Float16(1)) ≈ 0.84270079294971486934
    @test SF.erf(1) ≈ 0.84270079294971486934
    @test SF.erfc(1) ≈ 0.15729920705028513066
    @test SF.erfc(Float16(1)) ≈ 0.15729920705028513066
    @test SF.erfcx(1) ≈ 0.42758357615580700442
    @test SF.erfcx(Float32(1)) ≈ 0.42758357615580700442
    @test SF.erfcx(Complex64(1)) ≈ 0.42758357615580700442
    @test SF.erfi(1) ≈ 1.6504257587975428760
    @test SF.erfinv(0.84270079294971486934) ≈ 1
    @test SF.erfcinv(0.15729920705028513066) ≈ 1
    @test SF.dawson(1) ≈ 0.53807950691276841914

    @test SF.erf(1+2im) ≈ -0.53664356577856503399-5.0491437034470346695im
    @test SF.erfc(1+2im) ≈ 1.5366435657785650340+5.0491437034470346695im
    @test SF.erfcx(1+2im) ≈ 0.14023958136627794370-0.22221344017989910261im
    @test SF.erfi(1+2im) ≈ -0.011259006028815025076+1.0036063427256517509im
    @test SF.dawson(1+2im) ≈ -13.388927316482919244-11.828715103889593303im

    for elty in [Float32,Float64]
        for x in logspace(-200, -0.01)
            @test isapprox(SF.erf(SF.erfinv(x)), x, atol=1e-12*x)
            @test isapprox(SF.erf(SF.erfinv(-x)), -x, atol=1e-12*x)
            @test isapprox(SF.erfc(SF.erfcinv(2*x)), 2*x, atol=1e-12*x)
            if x > 1e-20
                xf = Float32(x)
                @test isapprox(SF.erf(SF.erfinv(xf)), xf, atol=1e-5*xf)
                @test isapprox(SF.erf(SF.erfinv(-xf)), -xf, atol=1e-5*xf)
                @test isapprox(SF.erfc(SF.erfcinv(2xf)), 2xf, atol=1e-5*xf)
            end
        end
        @test SF.erfinv(one(elty)) == Inf
        @test SF.erfinv(-one(elty)) == -Inf
        @test_throws DomainError SF.erfinv(convert(elty,2.0))

        @test SF.erfcinv(zero(elty)) == Inf
        @test_throws DomainError SF.erfcinv(-one(elty))
    end

    @test SF.erfinv(one(Int)) == SF.erfinv(1.0)
    @test SF.erfcinv(one(Int)) == SF.erfcinv(1.0)
end

@testset "airy" begin
    @test_throws SF.AmosException SF.airyai(200im)
    @test_throws SF.AmosException SF.airybi(200)

    for T in [Float32, Float64, Complex64,Complex128]
        @test SF.airyai(T(1.8)) ≈ 0.0470362168668458052247
        @test SF.airyaiprime(T(1.8)) ≈ -0.0685247801186109345638
        @test SF.airybi(T(1.8)) ≈ 2.595869356743906290060
        @test SF.airybiprime(T(1.8)) ≈ 2.98554005084659907283
    end
    for T in [Complex64, Complex128]
        z = convert(T,1.8 + 1.0im)
        @test SF.airyaix(z) ≈ SF.airyai(z) * exp(2/3 * z * sqrt(z))
        @test SF.airyaiprimex(z) ≈ SF.airyaiprime(z) * exp(2/3 * z * sqrt(z))
        @test SF.airybix(z) ≈ SF.airybi(z) * exp(-abs(real(2/3 * z * sqrt(z))))
        @test SF.airybiprimex(z) ≈ SF.airybiprime(z) * exp(-abs(real(2/3 * z * sqrt(z))))
    end
    @test_throws MethodError SF.airyai(complex(big(1.0)))

    for x = -3:3
        @test SF.airyai(x) ≈ SF.airyai(complex(x))
        @test SF.airyaiprime(x) ≈ SF.airyaiprime(complex(x))
        @test SF.airybi(x) ≈ SF.airybi(complex(x))
        @test SF.airybiprime(x) ≈ SF.airybiprime(complex(x))
        if x >= 0
            @test SF.airyaix(x) ≈ SF.airyaix(complex(x))
            @test SF.airyaiprimex(x) ≈ SF.airyaiprimex(complex(x))
        else
            @test_throws DomainError SF.airyaix(x)
            @test_throws DomainError SF.airyaiprimex(x)
        end
        @test SF.airybix(x) ≈ SF.airybix(complex(x))
        @test SF.airybiprimex(x) ≈ SF.airybiprimex(complex(x))
    end
end

@testset "bessel functions" begin
    bessel_funcs = [(SF.bessely0, SF.bessely1, SF.bessely), (SF.besselj0, SF.besselj1, SF.besselj)]
    @testset "$z, $o" for (z, o, f) in bessel_funcs
        @test z(Float32(2.0)) ≈ z(Float64(2.0))
        @test o(Float32(2.0)) ≈ o(Float64(2.0))
        @test z(2) ≈ z(2.0)
        @test o(2) ≈ o(2.0)
        @test z(2.0 + im) ≈ f(0, 2.0 + im)
        @test o(2.0 + im) ≈ f(1, 2.0 + im)
    end
    @testset "besselj error throwing" begin
        @test_throws MethodError SF.besselj(1.2,big(1.0))
        @test_throws MethodError SF.besselj(1,complex(big(1.0)))
        @test_throws MethodError SF.besseljx(1,big(1.0))
        @test_throws MethodError SF.besseljx(1,complex(big(1.0)))
    end
    @testset "besselh" begin
        true_h133 = 0.30906272225525164362 - 0.53854161610503161800im
        @test SF.besselh(3,1,3) ≈ true_h133
        @test SF.besselh(-3,1,3) ≈ -true_h133
        @test SF.besselh(3,2,3) ≈ conj(true_h133)
        @test SF.besselh(-3,2,3) ≈ -conj(true_h133)
        @testset "Error throwing" begin
            @test_throws SF.AmosException SF.besselh(1,0)
            @test_throws MethodError SF.besselh(1,big(1.0))
            @test_throws MethodError SF.besselh(1,complex(big(1.0)))
            @test_throws MethodError SF.besselhx(1,big(1.0))
            @test_throws MethodError SF.besselhx(1,complex(big(1.0)))
        end
    end
    @testset "besseli" begin
        true_i33 = 0.95975362949600785698
        @test SF.besseli(3,3) ≈ true_i33
        @test SF.besseli(-3,3) ≈ true_i33
        @test SF.besseli(3,-3) ≈ -true_i33
        @test SF.besseli(-3,-3) ≈ -true_i33
        @test SF.besseli(Float32(-3),Complex64(-3,0)) ≈ -true_i33
        @testset "Error throwing" begin
            @test_throws SF.AmosException SF.besseli(1,1000)
            @test_throws DomainError SF.besseli(0.4,-1.0)
            @test_throws MethodError SF.besseli(1,big(1.0))
            @test_throws MethodError SF.besseli(1,complex(big(1.0)))
            @test_throws MethodError SF.besselix(1,big(1.0))
            @test_throws MethodError SF.besselix(1,complex(big(1.0)))
        end
    end
    @testset "besselj" begin
        @test SF.besselj(0,0) == 1
        for i = 1:5
            @test SF.besselj(i,0) == 0
            @test SF.besselj(-i,0) == 0
            @test SF.besselj(-i,Float32(0)) == 0
            @test SF.besselj(-i,Float32(0)) == 0
        end

        j33 = SF.besselj(3,3.)
        @test SF.besselj(3,3) == j33
        @test SF.besselj(-3,-3) == j33
        @test SF.besselj(-3,3) == -j33
        @test SF.besselj(3,-3) == -j33
        @test SF.besselj(3,3f0) ≈ j33
        @test SF.besselj(3,complex(3.)) ≈ j33
        @test SF.besselj(3,complex(3f0)) ≈ j33
        @test SF.besselj(3,complex(3)) ≈ j33

        j43 = SF.besselj(4,3.)
        @test SF.besselj(4,3) == j43
        @test SF.besselj(-4,-3) == j43
        @test SF.besselj(-4,3) == j43
        @test SF.besselj(4,-3) == j43
        @test SF.besselj(4,3f0) ≈ j43
        @test SF.besselj(4,complex(3.)) ≈ j43
        @test SF.besselj(4,complex(3f0)) ≈ j43
        @test SF.besselj(4,complex(3)) ≈ j43

        @test j33 ≈ 0.30906272225525164362
        @test j43 ≈ 0.13203418392461221033
        @test SF.besselj(0.1, complex(-0.4)) ≈ 0.820421842809028916 + 0.266571215948350899im
        @test SF.besselj(3.2, 1.3+0.6im) ≈ 0.01135309305831220201 + 0.03927719044393515275im
        @test SF.besselj(1, 3im) ≈ 3.953370217402609396im
        @test SF.besselj(1.0,3im) ≈ SF.besselj(1,3im)
        @testset "Error throwing" begin
            @test_throws DomainError    SF.besselj(0.1, -0.4)
            @test_throws SF.AmosException SF.besselj(20,1000im)
            @test_throws MethodError SF.besselj(big(1.0),3im)
        end
    end

    @testset "besselk" begin
        true_k33 = 0.12217037575718356792
        @test SF.besselk(3,3) ≈ true_k33
        @test SF.besselk(-3,3) ≈ true_k33
        true_k3m3 = -0.1221703757571835679 - 3.0151549516807985776im
        @test SF.besselk(3,complex(-3)) ≈ true_k3m3
        @test SF.besselk(-3,complex(-3)) ≈ true_k3m3
        # issue #6564
        @test SF.besselk(1.0,0.0) == Inf
        @testset "Error throwing" begin
            @test_throws SF.AmosException SF.besselk(200,0.01)
            @test_throws DomainError SF.besselk(3,-3)
            @test_throws MethodError SF.besselk(1,big(1.0))
            @test_throws MethodError SF.besselk(1,complex(big(1.0)))
            @test_throws MethodError SF.besselkx(1,big(1.0))
            @test_throws MethodError SF.besselkx(1,complex(big(1.0)))
        end
    end

    @testset "bessely" begin
        y33 = SF.bessely(3,3.)
        @test SF.bessely(3,3) == y33
        @test SF.bessely(3.,3.) == y33
        @test SF.bessely(3,Float32(3.)) ≈ y33
        @test SF.bessely(-3,3) ≈ -y33
        @test y33 ≈ -0.53854161610503161800
        @test SF.bessely(3,complex(-3)) ≈ 0.53854161610503161800 - 0.61812544451050328724im
        @testset "Error throwing" begin
            @test_throws SF.AmosException SF.bessely(200.5,0.1)
            @test_throws DomainError SF.bessely(3,-3)
            @test_throws DomainError SF.bessely(0.4,-1.0)
            @test_throws DomainError SF.bessely(0.4,Float32(-1.0))
            @test_throws DomainError SF.bessely(1,Float32(-1.0))
            @test_throws DomainError SF.bessely(Cint(3),Float32(-3.))
            @test_throws DomainError SF.bessely(Cint(3),Float64(-3.))

            @test_throws MethodError SF.bessely(1.2,big(1.0))
            @test_throws MethodError SF.bessely(1,complex(big(1.0)))
            @test_throws MethodError SF.besselyx(1,big(1.0))
            @test_throws MethodError SF.besselyx(1,complex(big(1.0)))
        end
    end

    @testset "besselhx" begin
        for elty in [Complex64,Complex128]
            z = convert(elty, 1.0 + 1.9im)
            @test SF.besselhx(1.0, 1, z) ≈ convert(elty,-0.5949634147786144 - 0.18451272807835967im)
            @test SF.besselhx(Float32(1.0), 1, z) ≈ convert(elty,-0.5949634147786144 - 0.18451272807835967im)
        end
        @testset "Error throwing" begin
            @test_throws MethodError SF.besselh(1,1,big(1.0))
            @test_throws MethodError SF.besselh(1,1,complex(big(1.0)))
            @test_throws MethodError SF.besselhx(1,1,big(1.0))
            @test_throws MethodError SF.besselhx(1,1,complex(big(1.0)))
        end
    end
    @testset "scaled bessel[ijky] and hankelh[12]" begin
        for x in (1.0, 0.0, -1.0), y in (1.0, 0.0, -1.0), nu in (1.0, 0.0, -1.0)
            z = Complex128(x + y * im)
            z == zero(z) || @test SF.hankelh1x(nu, z) ≈ SF.hankelh1(nu, z) * exp(-z * im)
            z == zero(z) || @test SF.hankelh2x(nu, z) ≈ SF.hankelh2(nu, z) * exp(z * im)
            (nu < 0 && z == zero(z)) || @test SF.besselix(nu, z) ≈ SF.besseli(nu, z) * exp(-abs(real(z)))
            (nu < 0 && z == zero(z)) || @test SF.besseljx(nu, z) ≈ SF.besselj(nu, z) * exp(-abs(imag(z)))
            z == zero(z) || @test SF.besselkx(nu, z) ≈ SF.besselk(nu, z) * exp(z)
            z == zero(z) || @test SF.besselyx(nu, z) ≈ SF.bessely(nu, z) * exp(-abs(imag(z)))
        end
        @test SF.besselkx(1, 0) == Inf
        @testset "Error throwing" begin
            @test_throws SF.AmosException SF.hankelh1x(1, 0)
            @test_throws SF.AmosException SF.hankelh2x(1, 0)
            @test_throws SF.AmosException SF.besselix(-1, 0)
            @test_throws SF.AmosException SF.besseljx(-1, 0)
            @test_throws SF.AmosException SF.besselyx(1, 0)
            @test_throws DomainError SF.besselix(0.4,-1.0)
            @test_throws DomainError SF.besseljx(0.4, -1.0)
            @test_throws DomainError SF.besselkx(0.4,-1.0)
            @test_throws DomainError SF.besselyx(0.4,-1.0)
        end
    end
    @testset "issue #6653" begin
        @testset "$f" for f in (SF.besselj,SF.bessely,SF.besseli,SF.besselk,SF.hankelh1,SF.hankelh2)
            @test f(0,1) ≈ f(0,Complex128(1))
            @test f(0,1) ≈ f(0,Complex64(1))
        end
    end
end

# useful test functions for relative error, which differ from isapprox
# in that relerrc separately looks at the real and imaginary parts
relerr(z, x) = z == x ? 0.0 : abs(z - x) / abs(x)
relerrc(z, x) = max(relerr(real(z),real(x)), relerr(imag(z),imag(x)))
≅(a,b) = relerrc(a,b) ≤ 1e-13

@testset "gamma and friends" begin
    @testset "digamma" begin
        @testset "$elty" for elty in (Float32, Float64)
            @test SF.digamma(convert(elty, 9)) ≈ convert(elty, 2.140641477955609996536345)
            @test SF.digamma(convert(elty, 2.5)) ≈ convert(elty, 0.7031566406452431872257)
            @test SF.digamma(convert(elty, 0.1)) ≈ convert(elty, -10.42375494041107679516822)
            @test SF.digamma(convert(elty, 7e-4)) ≈ convert(elty, -1429.147493371120205005198)
            @test SF.digamma(convert(elty, 7e-5)) ≈ convert(elty, -14286.29138623969227538398)
            @test SF.digamma(convert(elty, 7e-6)) ≈ convert(elty, -142857.7200612932791081972)
            @test SF.digamma(convert(elty, 2e-6)) ≈ convert(elty, -500000.5772123750382073831)
            @test SF.digamma(convert(elty, 1e-6)) ≈ convert(elty, -1000000.577214019968668068)
            @test SF.digamma(convert(elty, 7e-7)) ≈ convert(elty, -1428572.005785942019703646)
            @test SF.digamma(convert(elty, -0.5)) ≈ convert(elty, .03648997397857652055902367)
            @test SF.digamma(convert(elty, -1.1)) ≈ convert(elty,  10.15416395914385769902271)

            @test SF.digamma(convert(elty, 0.1)) ≈ convert(elty, -10.42375494041108)
            @test SF.digamma(convert(elty, 1/2)) ≈ convert(elty, -γ - log(4))
            @test SF.digamma(convert(elty, 1)) ≈ convert(elty, -γ)
            @test SF.digamma(convert(elty, 2)) ≈ convert(elty, 1 - γ)
            @test SF.digamma(convert(elty, 3)) ≈ convert(elty, 3/2 - γ)
            @test SF.digamma(convert(elty, 4)) ≈ convert(elty, 11/6 - γ)
            @test SF.digamma(convert(elty, 5)) ≈ convert(elty, 25/12 - γ)
            @test SF.digamma(convert(elty, 10)) ≈ convert(elty, 7129/2520 - γ)
        end
    end

    @testset "trigamma" begin
        @testset "$elty" for elty in (Float32, Float64)
            @test SF.trigamma(convert(elty, 0.1)) ≈ convert(elty, 101.433299150792758817)
            @test SF.trigamma(convert(elty, 0.1)) ≈ convert(elty, 101.433299150792758817)
            @test SF.trigamma(convert(elty, 1/2)) ≈ convert(elty, π^2/2)
            @test SF.trigamma(convert(elty, 1)) ≈ convert(elty, π^2/6)
            @test SF.trigamma(convert(elty, 2)) ≈ convert(elty, π^2/6 - 1)
            @test SF.trigamma(convert(elty, 3)) ≈ convert(elty, π^2/6 - 5/4)
            @test SF.trigamma(convert(elty, 4)) ≈ convert(elty, π^2/6 - 49/36)
            @test SF.trigamma(convert(elty, 5)) ≈ convert(elty, π^2/6 - 205/144)
            @test SF.trigamma(convert(elty, 10)) ≈ convert(elty, π^2/6 - 9778141/6350400)
        end
    end

    @testset "invdigamma" begin
        @testset "$elty" for elty in (Float32, Float64)
            for val in [0.001, 0.01, 0.1, 1.0, 10.0]
                @test abs(SF.invdigamma(SF.digamma(convert(elty, val))) - convert(elty, val)) < 1e-8
            end
        end
        @test abs(SF.invdigamma(2)) == abs(SF.invdigamma(2.))
    end

    @testset "polygamma" begin
        @test SF.polygamma(20, 7.) ≈ -4.644616027240543262561198814998587152547
        @test SF.polygamma(20, Float16(7.)) ≈ -4.644616027240543262561198814998587152547
    end

    @testset "eta" begin
        @test SF.eta(1) ≈ log(2)
        @test SF.eta(2) ≈ pi^2/12
        @test SF.eta(Float32(2)) ≈ SF.eta(2)
        @test SF.eta(Complex64(2)) ≈ SF.eta(2)
    end
end

@testset "zeta" begin
    @test SF.zeta(0) ≈ -0.5
    @test SF.zeta(2) ≈ pi^2/6
    @test SF.zeta(Complex64(2)) ≈ SF.zeta(2)
    @test SF.zeta(4) ≈ pi^4/90
    @test SF.zeta(1,Float16(2.)) ≈ SF.zeta(1,2.)
    @test SF.zeta(1.,Float16(2.)) ≈ SF.zeta(1,2.)
    @test SF.zeta(Float16(1.),Float16(2.)) ≈ SF.zeta(1,2.)
    @test isnan(SF.zeta(NaN))
    @test isnan(SF.zeta(1.0e0))
    @test isnan(SF.zeta(1.0f0))
    @test isnan(SF.zeta(complex(0,Inf)))
    @test isnan(SF.zeta(complex(-Inf,0)))
end

#(compared to Wolfram Alpha)
@testset "digamma, trigamma, polygamma & zeta" begin
    for x in -10.2:0.3456:50
        @test 1e-12 > relerr(SF.digamma(x+0im), SF.digamma(x))
    end
    @test SF.digamma(7+0im) ≅ 1.872784335098467139393487909917597568957840664060076401194232
    @test SF.digamma(7im) ≅ 1.94761433458434866917623737015561385331974500663251349960124 + 1.642224898223468048051567761191050945700191089100087841536im
    @test SF.digamma(-3.2+0.1im) ≅ 4.65022505497781398615943030397508454861261537905047116427511+2.32676364843128349629415011622322040021960602904363963042380im
    @test SF.trigamma(8+0im) ≅ 0.133137014694031425134546685920401606452509991909746283540546
    @test SF.trigamma(8im) ≅ -0.0078125000000000000029194973110119898029284994355721719150 - 0.12467345030312762782439017882063360876391046513966063947im
    @test SF.trigamma(-3.2+0.1im) ≅ 15.2073506449733631753218003030676132587307964766963426965699+15.7081038855113567966903832015076316497656334265029416039199im
    @test SF.polygamma(2, 8.1+0im) ≅ -0.01723882695611191078960494454602091934457319791968308929600
    @test SF.polygamma(30, 8.1+2im) ≅ -2722.8895150799704384107961215752996280795801958784600407589+6935.8508929338093162407666304759101854270641674671634631058im
    @test SF.polygamma(3, 2.1+1im) ≅ 0.00083328137020421819513475400319288216246978855356531898998-0.27776110819632285785222411186352713789967528250214937861im
    @test 1e-11 > relerr(SF.polygamma(3, -4.2 + 2im),-0.0037752884324358856340054736472407163991189965406070325067-0.018937868838708874282432870292420046797798431078848805822im)
    @test SF.polygamma(13, 5.2 - 2im) ≅ 0.08087519202975913804697004241042171828113370070289754772448-0.2300264043021038366901951197725318713469156789541415899307im
    @test 1e-11 > relerr(SF.polygamma(123, -47.2 + 0im), 5.7111648667225422758966364116222590509254011308116701029e291)
    @test SF.zeta(4.1+0.3im, -3.2+0.1im) ≅ -281.34474134962502296077659347175501181994490498591796647 + 286.55601240093672668066037366170168712249413003222992205im
    @test SF.zeta(4.1+0.3im, 3.2+0.1im) ≅ 0.0121197525131633219465301571139288562254218365173899270675-0.00687228692565614267981577154948499247518236888933925740902im
    @test SF.zeta(4.1, 3.2+0.1im) ≅ 0.0137637451187986846516125754047084829556100290057521276517-0.00152194599531628234517456529686769063828217532350810111482im
    @test 1e-12 > relerrc(SF.zeta(1.0001, -4.5e2+3.2im), 10003.765660925877260544923069342257387254716966134385170 - 0.31956240712464746491659767831985629577542514145649468090im)
    @test SF.zeta(3.1,-4.2) ≅ SF.zeta(3.1,-4.2+0im) ≅ 149.7591329008219102939352965761913107060718455168339040295
    @test 1e-15 > relerrc(SF.zeta(3.1+0im,-4.2), SF.zeta(3.1,-4.2+0im))
    @test SF.zeta(3.1,4.2) ≅ 0.029938344862645948405021260567725078588893266227472565010234
    @test SF.zeta(27, 3.1) ≅ 5.413318813037879056337862215066960774064332961282599376e-14
    @test SF.zeta(27, 2) ≅ 7.4507117898354294919810041706041194547190318825658299932e-9
    @test 1e-12 > relerr(SF.zeta(27, -105.3), 1.3113726525492708826840989036205762823329453315093955e14)
    @test SF.polygamma(4, -3.1+Inf*im) == SF.polygamma(4, 3.1+Inf*im) == 0
    @test SF.polygamma(4, -0.0) == Inf == -SF.polygamma(4, +0.0)
    @test SF.zeta(4, +0.0) == SF.zeta(4, -0.0) ≅ pi^4 / 90
    @test SF.zeta(5, +0.0) == SF.zeta(5, -0.0) ≅ 1.036927755143369926331365486457034168057080919501912811974
    @test SF.zeta(Inf, 1.) == 1
    @test SF.zeta(Inf, 2.) == 0
    @test isnan(SF.zeta(NaN, 1.))
    @test isa([SF.digamma(x) for x in [1.0]], Vector{Float64})
    @test isa([SF.trigamma(x) for x in [1.0]], Vector{Float64})
    @test isa([SF.polygamma(3,x) for x in [1.0]], Vector{Float64})
    @test SF.zeta(2 + 1im, -1.1) ≅ SF.zeta(2 + 1im, -1.1+0im) ≅ -64.580137707692178058665068045847533319237536295165484548 + 73.992688148809018073371913557697318846844796582012921247im
    @test SF.polygamma(3,5) ≈ SF.polygamma(3,5.)

    @test SF.zeta(-3.0, 7.0) ≅ -52919/120
    @test SF.zeta(-3.0, -7.0) ≅ 94081/120
    @test SF.zeta(-3.1, 7.2) ≅ -587.457736596403704429103489502662574345388906240906317350719
    @test SF.zeta(-3.1, -7.2) ≅ 1042.167459863862249173444363794330893294733001752715542569576
    @test SF.zeta(-3.1, 7.0) ≅ -518.431785723446831868686653718848680989961581500352503093748
    @test SF.zeta(-3.1, -7.0) ≅ 935.1284612957581823462429983411337864448020149908884596048161
    @test SF.zeta(-3.1-0.1im, 7.2) ≅ -579.29752287650299181119859268736614017824853865655709516268 - 96.551907752211554484321948972741033127192063648337407683877im
    @test SF.zeta(-3.1-0.1im, -7.2) ≅ 1025.17607931184231774568797674684390615895201417983173984531 + 185.732454778663400767583204948796029540252923367115805842138im
    @test SF.zeta(-3.1-0.1im, 7.2 + 0.1im) ≅ -571.66133526455569807299410569274606007165253039948889085762 - 131.86744836357808604785415199791875369679879576524477540653im
    @test SF.zeta(-3.1-0.1im, -7.2 + 0.1im) ≅ 1035.35760409421020754141207226034191979220047873089445768189 + 130.905870774271320475424492384335798304480814695778053731045im
    @test SF.zeta(-3.1-0.1im, -7.0 + 0.1im) ≅ 929.546530292101383210555114424269079830017210969572819344670 + 113.646687807533854478778193456684618838875194573742062527301im
    @test SF.zeta(-3.1, 7.2 + 0.1im) ≅ -586.61801005507638387063781112254388285799318636946559637115 - 36.148831292706044180986261734913443701649622026758378669700im
    @test SF.zeta(-3.1, -7.2 + 0.1im) ≅ 1041.04241628770682295952302478199641560576378326778432301623 - 55.7154858634145071137760301929537184886497572057171143541058im
    @test SF.zeta(-13.4, 4.1) ≅ -3.860040842156185186414774125656116135638705768861917e6
    @test SF.zeta(3.2, -4) ≅ 2.317164896026427640718298719837102378116771112128525719078
    @test SF.zeta(3.2, 0) ≅ 1.166773370984467020452550350896512026772734054324169010977
    @test SF.zeta(-3.2+0.1im, 0.0) ≅ SF.zeta(-3.2+0.1im, 0.0+0im) ≅ 0.0070547946138977701155565365569392198424378109226519905493 + 0.00076891821792430587745535285452496914239014050562476729610im
    @test SF.zeta(-3.2, 0.0) ≅ SF.zeta(-3.2, 0.0+0im) ≅ 0.007011972077091051091698102914884052994997144629191121056378

    @test 1e-14 > relerr(SF.eta(1+1e-9), 0.693147180719814213126976796937244130533478392539154928250926)
    @test 1e-14 > relerr(SF.eta(1+5e-3), 0.693945708117842473436705502427198307157819636785324430166786)
    @test 1e-13 > relerr(SF.eta(1+7.1e-3), 0.694280602623782381522315484518617968911346216413679911124758)
    @test 1e-13 > relerr(SF.eta(1+8.1e-3), 0.694439974969407464789106040237272613286958025383030083792151)
    @test 1e-13 > relerr(SF.eta(1 - 2.1e-3 + 2e-3 * im), 0.69281144248566007063525513903467244218447562492555491581+0.00032001240133205689782368277733081683574922990400416791019im)
    @test 1e-13 > relerr(SF.eta(1 + 5e-3 + 5e-3 * im), 0.69394652468453741050544512825906295778565788963009705146+0.00079771059614865948716292388790427833787298296229354721960im)
    @test 1e-12 > relerrc(SF.zeta(1e-3+1e-3im), -0.5009189365276307665899456585255302329444338284981610162-0.0009209468912269622649423786878087494828441941303691216750im)
    @test 1e-13 > relerrc(SF.zeta(1e-4 + 2e-4im), -0.5000918637469642920007659467492165281457662206388959645-0.0001838278317660822408234942825686513084009527096442173056im)

    # Issue #7169: (TODO: better accuracy should be possible?)
    @test 1e-9  > relerrc(SF.zeta(0 + 99.69im), 4.67192766128949471267133846066040655597942700322077493021802+3.89448062985266025394674304029984849370377607524207984092848im)
    @test 1e-12 > relerrc(SF.zeta(3 + 99.69im), 1.09996958148566565003471336713642736202442134876588828500-0.00948220959478852115901654819402390826992494044787958181148im)
    @test 1e-9  > relerrc(SF.zeta(-3 + 99.69im), 10332.6267578711852982128675093428012860119184786399673520976+13212.8641740351391796168658602382583730208014957452167440726im)
    @test 1e-13 > relerrc(SF.zeta(2 + 99.69im, 1.3), 0.41617652544777996034143623540420694985469543821307918291931-0.74199610821536326325073784018327392143031681111201859489991im)
end

@testset "vectorization of 2-arg functions" begin
    binary_math_functions = [
        SF.besselh, SF.hankelh1, SF.hankelh2, SF.hankelh1x, SF.hankelh2x,
        SF.besseli, SF.besselix, SF.besselj, SF.besseljx, SF.besselk, SF.besselkx, SF.bessely, SF.besselyx,
        SF.polygamma, SF.zeta
    ]
    @testset "$f" for f in binary_math_functions
        x = y = 2
        v = [f(x,y)]
        @test f.([x],y) == v
        @test f.(x,[y]) == v
        @test f.([x],[y]) == v
    end
end

@testset "MPFR" begin
    @testset "bessel functions" begin
        setprecision(53) do
            @test SF.besselj(4, BigFloat(2)) ≈ SF.besselj(4, 2.)
            @test SF.besselj0(BigFloat(2)) ≈ SF.besselj0(2.)
            @test SF.besselj1(BigFloat(2)) ≈ SF.besselj1(2.)
            @test SF.bessely(4, BigFloat(2)) ≈ SF.bessely(4, 2.)
            @test SF.bessely0(BigFloat(2)) ≈ SF.bessely0(2.)
            @test SF.bessely1(BigFloat(2)) ≈ SF.bessely1(2.)
        end
    end

    let err(z, x) = abs(z - x) / abs(x)
        @test 1e-60 > err(SF.eta(parse(BigFloat,"1.005")), parse(BigFloat,"0.693945708117842473436705502427198307157819636785324430166786"))
        @test 1e-60 > err(exp(SF.eta(big(1.0))), 2.0)
    end

    let a = parse(BigInt, "315135")
        @test typeof(SF.erf(a)) == BigFloat
        @test typeof(SF.erfc(a)) == BigFloat
    end
end

@testset "Base Julia issue #17474" begin
    @test SF.f64(complex(1f0,1f0)) === complex(1.0, 1.0)
    @test SF.f64(1f0) === 1.0

    @test typeof(SF.eta(big"2")) == BigFloat
    @test typeof(SF.zeta(big"2")) == BigFloat
    @test typeof(SF.digamma(big"2")) == BigFloat

    @test_throws MethodError SF.trigamma(big"2")
    @test_throws MethodError SF.trigamma(big"2.0")
    @test_throws MethodError SF.invdigamma(big"2")
    @test_throws MethodError SF.invdigamma(big"2.0")

    @test_throws MethodError SF.eta(Complex(big"2"))
    @test_throws MethodError SF.eta(Complex(big"2.0"))
    @test_throws MethodError SF.zeta(Complex(big"2"))
    @test_throws MethodError SF.zeta(Complex(big"2.0"))
    @test_throws MethodError SF.zeta(1.0,big"2")
    @test_throws MethodError SF.zeta(1.0,big"2.0")
    @test_throws MethodError SF.zeta(big"1.0",2.0)
    @test_throws MethodError SF.zeta(big"1",2.0)

    @test typeof(SF.polygamma(3, 0x2)) == Float64
    @test typeof(SF.polygamma(big"3", 2f0)) == Float32
    @test typeof(SF.zeta(1, 2.0)) == Float64
    @test typeof(SF.zeta(1, 2f0)) == Float64 # BitIntegers result in Float64 returns
    @test typeof(SF.zeta(2f0, complex(2f0,0f0))) == Complex{Float32}
    @test typeof(SF.zeta(complex(1,1), 2f0)) == Complex{Float64}
    @test typeof(SF.zeta(complex(1), 2.0)) == Complex{Float64}
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests of titles.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Titles" begin
    title = "This is a very very long title that will be displayed above the table."

    expected = """

This is a very very long title that will be displayed above the table.

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data,
                          standalone = false,
                          tf = tf_html_default,
                          title = title)

    @test result == expected

    expected = """

This is a very very long title that will be displayed above the table.

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data,
                          standalone = false,
                          tf = tf_html_default,
                          title = title,
                          title_alignment = :c)

    @test result == expected

    expected = """

This is a very very long title that will be displayed above the table.

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1.0
    1

    2
    true
    2.0
    2

    3
    false
    3.0
    3

    4
    true
    4.0
    4

    5
    false
    5.0
    5

    6
    true
    6.0
    6

"""

    result = pretty_table(String, data,
                          standalone = false,
                          tf = tf_html_default,
                          title = title,
                          title_alignment = :r)

    @test result == expected
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Auxiliary functions for MGEO.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export conf_design_vars, conf_mgeo, convert, pf_to_csv, sort_pareto!

################################################################################
#                                  Overloads
################################################################################

"""
    convert(::Type{Matrix}, pareto::Vector{Pareto_Point{Nv,Nf}}) where {Nv,Nf}

Convert the structure that stores the Pareto frontier `pareto` to an array.

# Returns

A matrix of `Float64` in which the i-th line of the array is

    [var[1] var[2] ... var[n] f[1] f[2] ... f[nf]],

where `n` is the number of design variables and `nf` is the number of
objective functions.

"""
function convert(::Type{Matrix},
                 pareto::Vector{Pareto_Point{Nv,Nf}}) where {Nv,Nf}

    # Get the number of points.
    np = length(pareto)

    # Allocate the array.
    pf_array = zeros(np, Nv+Nf)

    # Fill the array.
    for i=1:np
        pf_array[i,:] = [pareto[i].vars' pareto[i].f']
    end

    pf_array
end

"""
    pf_to_csv(pareto::Vector{Pareto_Point{Nv,Nf}}, filename::String, design_vars::SVector{Nv,T}) where {Nv,Nf,T}

Save the Pareto frontier `pareto` to a CSV file `filename` considering the
configuration of the design variables in `design_vars`.

"""
function pf_to_csv(pareto::Vector{Pareto_Point{Nv,Nf}},
                   filename::String,
                   design_vars::SVector{Nv,T}) where {Nv,Nf,T}

    # Convert Pareto frontier.
    pf_array = convert(Matrix, pareto)

    # Array of variable names.
    var_names = Matrix{String}(undef, 1, Nv)

    for i = 1:Nv
        var_names[1,i] = design_vars[i].name
    end

    # Array of function names.
    func_names = Matrix{String}(undef, 1, Nf)

    for i = 1:Nf
        func_names[1,i] = "Obj. Func. $i"
    end

    # Write the information to the file.
    writedlm(filename, [ [ var_names func_names ]; pf_array ], ",")

    nothing
end

"""
    sort_pareto!(pareto::Vector{Pareto_Point{Nv,Nf}}, fobj::Int64) where {Nv,Nf}

Sort the points in Pareto frontier `pareto` using the objective function index
`fobj`. Notice that the sort will be performed in place, modifying the vector
`pareto`.

"""
function sort_pareto!(pareto::Vector{Pareto_Point{Nv,Nf}},
                      fobj::Int64) where {Nv,Nf}

    # Check the  input.
    nf = length(pareto[1].f)

    ( (fobj < 1) || (fobj > nf) ) && @error("fobj must be between 1 and the number of objective functions.")

    sort!(pareto, lt=(a,b)->begin
              # Check if data is valid.
              return (a.f[fobj] < b.f[fobj])
          end
         )

    nothing
end

################################################################################
#                              Private Functions
################################################################################

"""
    bitstrtonum(design_vars::SVector{Nv, T}, string::BitArray, factors::AbstractVector) where {Nv,T}

Convert the bit string `string` to real values using the information of the
design variables in `design_vars` and the integer factors `factors`.

# Returns

An array of `Float64` with the real values of the design variables.

"""
function bitstrtonum(design_vars::SVector{Nv, T},
                     string::BitArray,
                     factors::AbstractVector) where {Nv,T}

    # Array of variables.
    vars = Vector{Float64}(undef, Nv)

    # Loop for each design variable.
    @inbounds for i = 1:Nv
        # Convert the bits to integer.
        string_index_i = design_vars[i].index
        string_index_f = design_vars[i].index+design_vars[i].bits-1
        var_int = dot(factors[1:design_vars[i].bits],
                      string[string_index_i:string_index_f])

        # Convert the integer to real number.
        Δ       = (design_vars[i].max - design_vars[i].min)*var_int
        vars[i] = design_vars[i].min + Δ/(design_vars[i].full_scale)
    end

    SVector{Nv, Float64}(vars)
end

"""
    check_dominance!(pareto::Vector{Pareto_Point{Nv,Nf}}, candidate::Pareto_Point{Nv,Nf}, mgeod::MGEO_Structure{Nv,Nf}) where {Nv, Nf}

Check the dominance of the point `candidate` in the Pareto frontier `pareto`
using the MGEO configuration in `mgeod`. Notice that the variable `pareto` can
be modified depending on the dominance of the `candidate`.

# Returns

`true` if the point was added, `false` otherwise.

"""
function check_dominance!(pareto::Vector{Pareto_Point{Nv,Nf}},
                          candidate::Pareto_Point{Nv,Nf},
                          mgeod::MGEO_Structure{Nv,Nf}) where {Nv, Nf}
    # If Pareto frontier is empty, then add the point and return true.
    if length(pareto) == 0
        push!(pareto, candidate)
        return true
    end

    remove = Int64[]

    # Loop through the entire list of Pareto points.
    add_point = true

    @inbounds for p=1:length(pareto)
        # Variable to check if the point in the Pareto frontier list is
        # dominated by the candidate.
        point_dominated     = true
        candidate_dominated = true

        # Loop through objective functions.
        for i=1:Nf
            # If both objective functions are equal (given `mgeo_eps`), then
            # nothing can be stated about the dominance of both points.
            if abs(candidate.f[i] - pareto[p].f[i]) < mgeod.mgeo_eps
                continue
            end

            if candidate.f[i] < pareto[p].f[i]
                candidate_dominated = false
            elseif candidate.f[i] > pareto[p].f[i]
                point_dominated = false
            end
        end

        # If the point is dominated by the candidate and dominates the candidate
        # at the same time, then the candidate and the point is the same. Thus,
        # mark to do not add the candidate and exit the loop.

        if (point_dominated && candidate_dominated)
            add_point = false
            break
        end

        # If the point in the Pareto frontier is dominated, then mark it to be
        # excluded.
        if point_dominated
            push!(remove, p)
        # If the candidate is dominated by any point in the list, stop the
        # search and do not add it.
        elseif candidate_dominated
            add_point = false
            break
        end
    end

    # Remove the points dominated by the candidate point.
    @inbounds deleteat!(pareto, remove)

    # Check if the point must be added.
    add_point && push!(pareto, candidate)

    # Return if a point was added.
    add_point
end


text:
@testset "fsm_stream_reassembler_single.jl" begin
    @testset begin
        sr = StreamReassembler(65000)
        @test sr.cur_idx == 0
        @test read!(sr.output) == ""
        @test sr.output |> !eof
    end

    @testset begin
        sr = StreamReassembler(65000)
        push_substring!(sr, "a", 0)
        @test sr.cur_idx == 1
        @test read!(sr.output) == "a"
        @test sr.output |> !eof
    end

    @testset begin
        sr = StreamReassembler(65000)
        push_substring!(sr, "a", 0, true)
        @test sr.cur_idx == 1
        @test read!(sr.output) == "a"
        @test sr.output |> eof
    end

    @testset begin
        sr = StreamReassembler(65000)
        push_substring!(sr, "", 0, true)
        @test sr.cur_idx == 0
        @test read!(sr.output) == ""
        @test sr.output |> eof
    end

    @testset begin
        sr = StreamReassembler(65000)
        push_substring!(sr, "b", 0, true)
        @test sr.cur_idx == 1
        @test read!(sr.output) == "b"
        @test sr.output |> eof
    end

    @testset begin
        sr = StreamReassembler(65000)
        push_substring!(sr, "", 0)
        @test sr.cur_idx == 0
        @test read!(sr.output) == ""
        @test sr.output |> !eof
    end

    @testset begin
        sr = StreamReassembler(8)
        push_substring!(sr, "abcdefgh", 0)
        @test sr.cur_idx == 8
        @test read!(sr.output) == "abcdefgh"
        @test sr.output |> !eof
    end

    @testset begin
        sr = StreamReassembler(8)
        push_substring!(sr, "abcdefgh", 0, true)
        @test sr.cur_idx == 8
        @test read!(sr.output) == "abcdefgh"
        @test sr.output |> eof
    end

    @testset begin
        sr = StreamReassembler(8)
        push_substring!(sr, "abc", 0)
        @test sr.cur_idx == 3
        push_substring!(sr, "bcdefgh", 1, true)

        @test read!(sr.output) == "abcdefgh"
        @test sr.output |> eof
    end

    @testset begin
        sr = StreamReassembler(8)
        push_substring!(sr, "abc", 0)
        @test sr.cur_idx == 3
        @test sr.output |> !eof

        push_substring!(sr, "ghX", 6, true)
        @test sr.cur_idx == 3
        @test sr.output |> !eof

        push_substring!(sr, "cdefg", 2)
        @test sr.cur_idx == 8
        @test sr.output |> !eof

        @test read!(sr.output) == "abcdefgh"
        @test sr.output |> !eof
    end

    @testset begin
        sr = StreamReassembler(8)
        push_substring!(sr, "abc", 0)
        @test sr.cur_idx == 3
        @test sr.output |> !eof

        push_substring!(sr, "", 6)
        @test sr.cur_idx == 3
        @test sr.output |> !eof

        push_substring!(sr, "de", 3, true)
        @test sr.cur_idx == 5

        @test read!(sr.output) == "abcde"
        @test sr.output |> eof
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
mutable struct MWSingleLayerTDIO{T} <: RetardedPotential{T}
    "speed of light in medium"
    speed_of_light::T
    "weight of the weakly singular term"
    ws_weight::T
    "weight of the hypersingular term"
    hs_weight::T
    "number of temporal differentiations in the weakly singular term"
    ws_diffs::Int
    "number of temporal differentiations in the hyper singular term"
    hs_diffs::Int
end

function Base.:*(a::Number, op::MWSingleLayerTDIO)
	@info "scalar product a * op (SL)"
	MWSingleLayerTDIO(
		op.speed_of_light,
		a * op.ws_weight,
		a * op.hs_weight,
		op.ws_diffs,
		op.hs_diffs)
end

mutable struct MWDoubleLayerTDIO{T} <: RetardedPotential{T}
    speed_of_light::T
    weight::T
    num_diffs::Int
end

function Base.:*(a::Number, op::MWDoubleLayerTDIO)
	@info "scalar product a * op (DL)"
	MWDoubleLayerTDIO(
		op.speed_of_light,
		a * op.weight,
		op.num_diffs)
end

mutable struct MWDoubleLayerTransposedTDIO{T} <: RetardedPotential{T}
	speed_of_light::T
    weight::T
    num_diffs::Int
end

function Base.:*(a::Number, op::MWDoubleLayerTransposedTDIO)
	@info "scalar product a * op (DL)"
	MWDoubleLayerTransposedTDIO(
		op.speed_of_light,
		a * op.weight,
		op.num_diffs)
end

MWSingleLayerTDIO(;speedoflight) = MWSingleLayerTDIO(speedoflight,-1/speedoflight,-speedoflight,2,0)
MWDoubleLayerTDIO(;speedoflight) = MWDoubleLayerTDIO(speedoflight, one(speedoflight), 0)

module TDMaxwell3D
import ...BEAST

function singlelayer(;speedoflight, numdiffs=0)
	@assert numdiffs >= 0
	numdiffs == 0 && return BEAST.integrate(BEAST.MWSingleLayerTDIO(speedoflight,-1/speedoflight,-speedoflight,2,0))
	return BEAST.MWSingleLayerTDIO(speedoflight,-1/speedoflight,-speedoflight,2+numdiffs-1,numdiffs-1)
end

function doublelayer(;speedoflight, numdiffs=0)
	@assert numdiffs >= -1
	numdiffs == -1 && BEAST.integrate(BEAST.MWDoubleLayerTDIO(speedoflight,1.0,0))
	return BEAST.MWDoubleLayerTDIO(speedoflight,1.0,numdiffs)
end

end # module TDMaxwell3D

export TDMaxwell3D

function quaddata(op::MWSingleLayerTDIO, testrefs, trialrefs, timerefs,
        testels, trialels, timeels)

    dmax = numfunctions(timerefs)-1
    bn = binomial.((0:dmax),(0:dmax)')

    V = eltype(testels[1].vertices)
    ws = WiltonInts84.workspace(V)
    # quadpoints(testrefs, testels, (3,)), bn, ws
    quadpoints(testrefs, testels, (3,)), bn, ws
end

quadrule(op::MWSingleLayerTDIO, testrefs, trialrefs, timerefs,
        p, testel, q, trialel, r, timeel, qd) = WiltonInts84Strat(qd[1][1,p],qd[2],qd[3])

struct TransposedStorage{F}
	store::F
end

@inline (f::TransposedStorage)(v,m,n,k) = f.store(v,n,m,k)

function allocatestorage(op::MWDoubleLayerTDIO, testST, basisST,
	::Type{Val{:bandedstorage}},
	::Type{LongDelays{:ignore}})

    # tfs = spatialbasis(testST)
    # bfs = spatialbasis(basisST)
	X, T = spatialbasis(testST), temporalbasis(testST)
	Y, U = spatialbasis(basisST), temporalbasis(basisST)

	if CompScienceMeshes.refines(geometry(Y), geometry(X))
		testST = Y⊗T
		basisST = X⊗U
	end

    M = numfunctions(X)
    N = numfunctions(Y)

    K0 = fill(typemax(Int), M, N)
    K1 = zeros(Int, M, N)

    function store(v,m,n,k)
        K0[m,n] = min(K0[m,n],k)
        K1[m,n] = max(K1[m,n],k)
    end

    aux = EmptyRP(op.speed_of_light)
    print("Allocating memory for convolution operator: ")
    assemble!(aux, testST, basisST, store)
    println("\nAllocated memory for convolution operator.")

	maxk1 = maximum(K1)
	bandwidth = maximum(K1 .- K0 .+ 1)
	data = zeros(eltype(op), bandwidth, M, N)
	Z = SparseND.Banded3D(K0, data, maxk1)
    store1(v,m,n,k) = (Z[m,n,k] += v)
    return ()->Z, store1
end

function assemble!(dl::MWDoubleLayerTDIO, W::SpaceTimeBasis, V::SpaceTimeBasis, store,
    threading=Threading{:multi})

	X, T = spatialbasis(W), temporalbasis(W)
	Y, U = spatialbasis(V), temporalbasis(V)
	if CompScienceMeshes.refines(geometry(Y), geometry(X))
		@assert !CompScienceMeshes.refines(geometry(X), geometry(Y))
		W = Y⊗T
		V = X⊗U
		op = MWDoubleLayerTransposedTDIO(dl.speed_of_light, dl.weight, dl.num_diffs)
		assemble!(op, W, V, store)
		return
	end

	P = Threads.nthreads()
	Y, S = spatialbasis(W), temporalbasis(W)
	splits = [round(Int,s) for s in range(0, stop=numfunctions(Y), length=P+1)]

	@info "Starting assembly with $P threads:"
	Threads.@threads for i in 1:P
		lo, hi = splits[i]+1, splits[i+1]
		lo <= hi || continue
		Y_p = subset(Y, lo:hi)
		store2 = (v,m,n,k) -> store(v,lo+m-1,n,k)
		assemble_chunk!(dl, Y_p ⊗ S, V, store2)
	end

	# return assemble_chunk!(dl, W, V, store1)
end

function quaddata(op::MWDoubleLayerTDIO, testrefs, trialrefs, timerefs,
        testels, trialels, timeels)

    dmax = numfunctions(timerefs)-1
    bn = binomial.((0:dmax),(0:dmax)')

    V = eltype(testels[1].vertices)
    ws = WiltonInts84.workspace(V)
    # quadpoints(testrefs, testels, (3,)), bn, ws
    quadpoints(testrefs, testels, (3,)), bn, ws
end

quadrule(op::MWDoubleLayerTDIO, testrefs, trialrefs, timerefs,
        p, testel, q, trialel, r, timeel, qd) = WiltonInts84Strat(qd[1][1,p],qd[2],qd[3])

function quaddata(op::MWDoubleLayerTransposedTDIO,
		testrefs, trialrefs, timerefs,
        testels, trialels, timeels)

    dmax = numfunctions(timerefs)-1
    bn = binomial.((0:dmax),(0:dmax)')

    V = eltype(testels[1].vertices)
    ws = WiltonInts84.workspace(V)
    quadpoints(testrefs, testels, (3,)), bn, ws
end

quadrule(op::MWDoubleLayerTransposedTDIO, testrefs, trialrefs, timerefs,
        p, testel, q, trialel, r, timeel, qd) = WiltonInts84Strat(qd[1][1,p],qd[2],qd[3])

function momintegrals!(z, op::MWDoubleLayerTransposedTDIO,
	g, f, T, τ, σ, ι, qr::WiltonInts84Strat)

	op1 = MWDoubleLayerTDIO(op.speed_of_light, op.weight, op.num_diffs)
	momintegrals!(z, op1, g, f, T, τ, σ, ι, qr::WiltonInts84Strat)
	w = similar(z)
	permutedims!(w, z, [2,1,3])

end

@inline function tmRoR(d, iG, bns)
    sgn = isodd(d) ? -1 : 1
    r = sgn * iG[d+2]
end

# build
# ``\int (D-R)^d/R (y-b) dy`` from
# ``(ξ-b) \int R^k dy`` and
# ``\int R^k (y-ξ) dy``
@inline function tmRoRf(d, iG, iGξy, bξ, bns)
    sgn = isodd(d) ? -1 : 1
    iGf = iGξy[d+2] + bξ * iG[d+2]
    r = sgn * iGf
end

"""
    Q = qd(T,dh,::Val{N})

Q[k] is the factor in front resulting from differentiating t^(k-1) dh times.
"""
@generated function qh(::Type{T},dh,n::Type{Val{N}}) where {N,T}
    xp = quote end
    for k in 1:N
        qk = Symbol(:q,k)
        d = k-1
        xp1 = quote
            $(qk) = one($T)
            for v in 0 : dh-1
                $(qk) *= ($d-v)
            end
        end
        append!(xp.args, xp1.args)
    end
    xp1 = :(())
    for k in 1:N
        qk = Symbol(:q,k)
        push!(xp1.args, :($qk))
    end
    push!(xp.args, xp1)
    return xp
end

function innerintegrals!(zl, op::MWSingleLayerTDIO,
        p, # test_point, test_time
        U, V, W, # local_test_space, local_trial_space, local_temporal_space
        τ, σ, ι, # test_element, trial_element, spherial_shell
        qr, w)   # inner_quadrature_rule, outer_quadrature_weight

	T = typeof(w)

    sol = op.speed_of_light
    #Rmax = sol * tmax

    dx = w
    x = cartesian(p)
    n = cross(σ[1]-σ[3],σ[2]-σ[3])
    n /= norm(n)
    ξ = x - ((x-σ[1]) ⋅ n) * n

    r = ι[1]
    R = ι[2]

    @assert r < R
    @assert degree(W) <= 3

    ∫G, ∫Gξy, = WiltonInts84.wiltonints(σ[1],σ[2],σ[3],x,r,R,Val{2},qr.workspace)

    αg = 1 / volume(τ) / 2
	αf = 1 / volume(σ) / 2
	αG = 1 / 4π
	α = αg * αf * αG * op.ws_weight * dx
	β = 4 * αg * αf * αG * op.hs_weight * dx

	ds = op.ws_diffs
	dh = op.hs_diffs

    qhs = qh(T,dh,Val{4})
    qss = qh(T,ds,Val{4})

    bn = qr.binomials

    #solpowers = collect(sol^p for p ∈ 0:numfunctions(W)-1)
    sol2 = sol*sol
    sol3 = sol2*sol
    sol4 = sol3*sol
    sol5 = sol4*sol
    solpowers = (one(sol), sol, sol2, sol3, sol4, sol5)

    for i in 1 : numfunctions(U)
        a = τ[i]
        g = (x-a)
        for j in 1 : numfunctions(V)
            b = σ[j]; bξ = ξ-b
            for k in 1 : numfunctions(W)
				d = k-1 # ranges from 0 to numfunctions(W)-1
				sgn = isodd(d) ? -1 : 1
				# hyper singular contribution
				if d >= dh
                    @assert dh == 0
                    q = qhs[k]
                    Ih = tmRoR(d-dh, ∫G, bn) # \int (cTmax-R)^(d-dh)/R dy
                    #zl[i,j,k] += β * q * Ih / sol^(d-dh)
                    zl[i,j,k] += β * q * Ih / solpowers[d-dh+1]
				end
				# weakly singular contribution
				if d >= ds
                    q = qss[k]
                    Is = tmRoRf(d-ds, ∫G, ∫Gξy, bξ, bn) # \int (cTmax-R)^(d-ds)/R (y-b) dy
                    #zl[i,j,k] += α * q * (g ⋅ Is) / sol^(d-ds)
                    zl[i,j,k] += α * q * (g ⋅ Is) / solpowers[d-ds+1]
				end
            end
        end
    end

end

function innerintegrals!(z, op::MWDoubleLayerTDIO,
    p,
    U, V, W,
    τ, σ, ι,
    qr, w)

	T = typeof(w)

    sol = op.speed_of_light

    dx = w
    x = cartesian(p)
    n = cross(σ[1]-σ[3],σ[2]-σ[3])
    n /= norm(n)
    ξ = x - ((x-σ[1]) ⋅ n) * n

    r = ι[1]
    R = ι[2]

    @assert r < R
    @assert degree(W) <= 3

    #N = max(degree(W), 0)
    ∫G, ∫Gξy, ∫∇G = WiltonInts84.wiltonints(σ[1],σ[2],σ[3],x,r,R,Val{2},qr.workspace)
	@assert isapprox(∫∇G[2] , point(0,0,0), atol=1e-8)

    αg = 1 / volume(τ) / 2
	αf = 1 / volume(σ) / 2
	αG = 1 / 4 / π
	α = αG * op.weight * dx # * αg * αf

	ds = op.num_diffs

    @inline function tmGRoR(d, iGG)
        sgn = isodd(d) ? -1 : 1
        r = sgn * iGG[d+1]
    end

    Ux = U(p)
    Vx = αf * @SVector[(x-σ[1]), (x-σ[2]), (x-σ[3])]

    for i in 1 : numfunctions(U)
        # a = τ[i]
        # g = αg * (x-τ[i])
        g = Ux[i].value
        for j in 1 : numfunctions(V)
            # b = σ[j]
            # f = αf * (x-σ[j])
            # f = Vx[j].value
            f = Vx[j]
            fxg = f × g
            for k in 1 : numfunctions(W)
				d = k-1
				sgn = isodd(d) ? -1 : 1
				if d >= ds
					q = one(T)
					for p in 0 : ds-1
						q *= (d-p)
					end
                    # @assert q == 1
                    z[i,j,k] += -α * q * ( fxg ⋅ tmGRoR(d-ds, ∫∇G) ) / sol^(d-ds)
				end
            end
        end
    end

end