text:
using JuMP
using PiecewiseLinearOpt

m = Model()
@variable(m, x)
d = [0, 1, 2, 3, 4]
f(x) = x^2
fd = [f(x) for x in d]
println(all_variables(m))

PiecewiseLinearOpt.piecewiselinear(m, x, d, fd; method = :DLog)

println(all_variables(m))

ext: jl


text:
using Pkg;

Pkg.add("DifferentialEquations")
Pkg.add("ModelingToolkit")
Pkg.add("DynamicalSystems")
Pkg.add("ForwardDiff")
Pkg.add("Graphs")
Pkg.add("Plots")
Pkg.add("Pluto")
Pkg.add("PlutoUI")
Pkg.add("LaTeXStrings")

max_stars_repo_name: csimal/Julia-for-Dynamical-Systems

max_issues_repo_name: csimal/Julia-for-Dynamical-Systems

max_forks_repo_name: csimal/Julia-for-Dynamical-Systems


text:
using PowerModels
using Ipopt
using JuMP
using .PP2PM

function run_powermodels(json_path)
    pm = PP2PM.load_pm_from_json(json_path)
    ipopt_solver = JuMP.with_optimizer(Ipopt.Optimizer, print_level=0)

    result = PowerModels.run_ac_opf(pm, ipopt_solver,
                                    setting = Dict("output" => Dict("branch_flows" => true)))
    return result
end

max_stars_repo_name: Zamwell/pandapower

max_issues_repo_name: Zamwell/pandapower

max_forks_repo_name: Zamwell/pandapower


text:
using StrideArrays, LinearAlgebra, Aqua, ArrayInterface
using Test

import InteractiveUtils
InteractiveUtils.versioninfo(stdout; verbose = true)

@show StrideArrays.VectorizationBase.register_count()
const START_TIME = time()

@time @testset "StrideArrays.jl" begin
    Aqua.test_all(StrideArrays)
    # @test isempty(detect_unbound_args(StrideArrays))
    @time include("matmul_tests.jl")
    @time include("misc.jl")
    @time include("broadcast_tests.jl")
end

const ELAPSED_MINUTES = (time() - START_TIME)/60
@test ELAPSED_MINUTES < 180


text:
module BounWrite
export setboun, timelimit
include("PathFile.jl");
##############################

"""

    setboun(times, heads, location, FEHMprefix)

Sets boundary conditions in the `FEHMprefix.boun` file for FEHM simulations.

Arguments:
- `times`: vector of years at which pumping rate from well is changed
- `hears`: pumping rate.
- `location`: location (FEHM code) of pumping well.
- `FEHMprefix` : FEHM files prefix in `./contaminantion-model/

"""

function setboun(times, heads, location, FEHMprefix)

ff = open(string(fehmdir, FEHMprefix, ".boun"),"w+");
write(ff,"boun
model 3 w21
tran
time\n");
write(ff,string(length(times)));
write(ff,"\n");
times = times./8;
for i = 1:length(times)
	tt = times[i];
	write(ff,"$tt \n");
end
write(ff,"year \ndsw \n");
for i = 1:length(heads)
	hh = heads[i];
	write(ff, "$hh \n");
end
write(ff,"end\n");
write(ff, "-");
write(ff, string(location));
write(ff, " 0 0 1\n ");
write(ff," ");
close(ff);

end

##################################################

"""

    timelimit(days, location, FEHMprefix)

Sets simulation horizing in the `FEHMprefix.data` FEHM file.

Argument:

- `days` : number of days limit
- `location` : location (FEHM code) of pumping well.
- `FEHMprefix` : FEHM files prefix in `./contaminantion-model/

"""

function timelimit(days, location, FEHMprefix)

gg = open(string(fehmdir,FEHMprefix,".data"),"w+");
write(gg,"** Cr model
cont
avs 1 1e20
geo
liquid
concentration
head
velocity
material
endavs
head 1000
nobr
bous
1
airw
3
30.0 0.1
zone
file
gs01_outside.zone
flow
-3 0 0 1798.613472817     1 1e10
-5 0 0 1772.054959807     1 1e10

pres
1 0 0 1776 1 2

zone
110
4.971500000E+05 501750          501750          4.971500000E+05 4.971500000E+05 501750          501750          4.971500000E+05
5.373500000E+05 5.373500000E+05 5.406500000E+05 5.406500000E+05 5.373500000E+05 5.373500000E+05 5.406500000E+05 5.406500000E+05
-1 -1 -1 -1 1 1 1 1

flow
-110 0 0 0.000717864089902394     1 0

zonn
file
hole_coord-20120702-recharge-600p.zonn
flow
-9420 0 0 -0.0003413616876478    1 0
-9620 0 0 -0.00007255554857619   1 0

hyco
1 0 0 -2.51795 -2.51795 -2.51795

rock
1 0 0 2500. 1010. 1.0E-02

rlp
1 0. 0. 1. 1. 0. 1.

1 0 0 1

ppor
-1
1 0 0 0.000215968

stea
shea 1e-6
stim 1e20
smul 1e10
end
time
365.25 $days 10000 001 2006 1 0
0.0  -1.2 1.0 365.25

ctrl
-5 1e-6 1 350 bcgs
1 0 0 2

1.0 0 1.0
10 1000 1e-8 1e20
0 2
sol
1 -1
iter
1e-5 1e-5 0.1 1e-4 1.1
13 0 0 5 1200.
node
1
5471
zone
file
gs01_outside.zone
zonn
file
hole_coord-20140523-screen-change-100.zonn
flxz
1\n")
write(gg,"$location")
write(gg,"\nboun
file\n")
write(gg,FEHMprefix)
write(gg,".boun
zonn
file\n")
write(gg,FEHMprefix)
write(gg,".source.zonn
trac
file\n")
write(gg,FEHMprefix)
write(gg,".trac
stop");
close(gg);
end
end

max_stars_repo_name: josephbakarji/socio-hydro-abm

max_issues_repo_name: josephbakarji/socio-hydro-abm

max_forks_repo_name: josephbakarji/socio-hydro-abm


text:
#Function for calculating Rydberg state lifetimes
#23/07/2017

function Radiative_lifetimes(atom,nn,ll,jj)

  temp = 300 #Approximate room temperature in K

  if atom=="87Rb"
    ZZ, spin, eneigval, params, delta_nlj = Rb87numbers(nn,ll,jj)
  end
  dict = load(string("Data/",atom,"levels.jld"))

  #Radiative lifetimes
  if atom=="87Rb"
    loadparam, stateenvec, Rdecayratevec, qvec1, qvec4 = Rb87DecayPaths(ll,jj)
  end

  for qq = qvec1:qvec4
    #Get matrix element
    matrixelement, matrixelementSI = matrix_elements(atom,nn,ll,jj,loadparam[qq,1],loadparam[qq,2],loadparam[qq,3])
    #Compute energies
    enlow = stateenvec[qq]-dict["ionlim"] #Energy of lower state in eV
    eneigvaleV = eneigval*atomenergy #Energy of state |n,l,j⟩ in eV
    #Calculate lifetime excluding black body radiation
    omega = abs((eneigvaleV-enlow)/rpceV) #Angular frequency of transimission
    term1 = (omega^3)/(3*pi*vacpmtvty*rpcJ*(lightc^3)) #splitting terms to make reading easier
    term2 = (2*jj+1)/(2*loadparam[qq,3]+1)

    Rdecayratevec[qq] = term1*term2*matrixelementSI^2 #Radiative decay rate
  end

  Rdecayrate = sum(Rdecayratevec)

  #Account for black body radiation following Beterov et al. (2009)

  neff = nn - delta_nlj #Effective principal quantum number

  if atom=="87Rb"
    BBdecayrate = Rb87blackbody(neff,temp,ll,jj)
  end

  decayrate = Rdecayrate+BBdecayrate #total decay rate
  lifetime = 1/decayrate

  return lifetime
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
function halo_regions_initalized_correctly(arch, FT, Nx, Ny, Nz)
    # Just choose something anisotropic to catch Δx/Δy type errors.
    Lx, Ly, Lz = 10, 20, 30

    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (Lx, Ly, Lz))
    field = CellField(FT, arch, grid)

    # Fill the interior with random numbers.
    data(field) .= rand(FT, Nx, Ny, Nz)

    Hx, Hy, Hz = grid.Hx, grid.Hy, grid.Hz

    # The halo regions should still just contain zeros.
    (all(field.data[1-Hx:0,          :,          :] .== 0) &&
     all(field.data[Nx+1:Nx+Hx,      :,          :] .== 0) &&
     all(field.data[:,          1-Hy:0,          :] .== 0) &&
     all(field.data[:,      Ny+1:Ny+Hy,          :] .== 0) &&
     all(field.data[:,               :,     1-Hz:0] .== 0) &&
     all(field.data[:,               :, Nz+1:Nz+Hz] .== 0))
end

function halo_regions_correctly_filled(arch, FT, Nx, Ny, Nz)
    # Just choose something anisotropic to catch Δx/Δy type errors.
    Lx, Ly, Lz = 100, 200, 300

    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (Lx, Ly, Lz))
    field = CellField(FT, arch, grid)
    fbcs = HorizontallyPeriodicBCs()

    data(field) .= rand(FT, Nx, Ny, Nz)
    fill_halo_regions!(field.data, fbcs, arch, grid)

    Hx, Hy, Hz = grid.Hx, grid.Hy, grid.Hz

    (all(field.data[1-Hx:0,   1:Ny,   1:Nz] .== field.data[Nx-Hx+1:Nx, 1:Ny,       1:Nz]) &&
     all(field.data[1:Nx,   1-Hy:0,   1:Nz] .== field.data[1:Nx,      Ny-Hy+1:Ny,  1:Nz]) &&
     all(field.data[1:Nx,     1:Ny,     1-Hz:0] .== field.data[1:Nx,        1:Ny,   1:1]) &&
     all(field.data[1:Nx,     1:Ny, Nz+1:Nz+Hz] .== field.data[1:Nx,        1:Ny, Nz:Nz]))
end

@testset "Halo regions" begin
    println("Testing halo regions...")

    Ns = [(8, 8, 8), (8, 8, 4), (10, 7, 5),
          (1, 8, 8), (1, 9, 5),
          (8, 1, 8), (5, 1, 9),
          (8, 8, 1), (5, 9, 1),
          (1, 1, 8)]

    @testset "Initializing halo regions" begin
        println("  Testing initializing halo regions...")
        for arch in archs, FT in float_types, N in Ns
            @test halo_regions_initalized_correctly(arch, FT, N...)
        end
    end

    @testset "Filling halo regions" begin
        println("  Testing filling halo regions...")
        for arch in archs, FT in float_types, N in Ns
            @test halo_regions_correctly_filled(arch, FT, N...)
        end
    end
end


text:
#=
julia> fib(n) = n < 3 ? 1 : fib(n-1) + fib(n-2)

julia> fib.(1:10)
10-element Array{Int64,1}:
  1
  1
  2
  3
  5
  8
 13
 21
 34
 55
 =#

function fib(n)
    if n < 3
        return (result = 1, counter = 1)
    else
        result1, counter1 = fib(n - 1)
        result2, counter2 = fib(n - 2)
        return (result = result1 + result2, counter = 1 + counter1 + counter2)
    end
end

# How many times does the function get called?
#=
julia> fib(6)
(result = 8, counter = 15)

julia> fib(10)
(result = 55, counter = 109)

julia> fib(20)
(result = 6765, counter = 13529)
=#

# How well does it perform?
#=
julia> @btime fib(40)
  569.906 ms (0 allocations: 0 bytes)
(result = 102334155, counter = 204668309)
=#

const fib_cache = Dict()

_fib(n) = n < 3 ? 1 : fib(n-1) + fib(n-2)

function fib(n)
    if haskey(fib_cache, n)
        return fib_cache[n]
    else
        value = _fib(n)
        fib_cache[n] = value
        return value
    end
end

# Benchmark performance. Make sure cache is emptied at each
#=
julia> @btime fib(40) setup=(empty!(fib_cache))
  38.731 ns (0 allocations: 0 bytes)
102334155
=#


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ a6c0b440-5688-11eb-0635-39912417c1f9
begin
	using NeXLDatabase, Gadfly, DataFrames, PlutoUI
md"""
## Plot k-ratios
This script makes it easy to plot some or all of the k-ratios associated with a single element over different measurement sets.
  1) Select an element
  2) Specify a minimum k-ratio to consider (eliminates materials with trace quantities of the element)
  3) Select one or more sets of measurements to display.  Measurements are as identified by material, standard material and beam energy.
"""
end

# ╔═╡ 3742ed50-1852-11eb-01e3-5b27d6c85250
db=openNeXLDatabase("C:\\Users\\nritchie\\Documents\\DrWatson\\K-Ratio Project\\data\\exp_pro\\kratio.db")

# ╔═╡ e793f370-1852-11eb-24b6-05fa57c3a3c6
md"""
Select an element or elements: $(@bind els MultiSelect(map(el->symbol(el),elements[1:94]),default=["Al"]))
"""

# ╔═╡ 78885a50-1859-11eb-0b6c-3368b717d267
md"""
Specify the minimum k-ratio: $(@bind mink NumberField(0.00 : 0.01 : 1.0; default=0.1))
"""

# ╔═╡ a51b0a40-1921-11eb-2c1c-ab31e48e7846
allkrs = mapreduce(el->findall(db, DBKRatio, elm=parse(Element,el), mink=mink), append!, els)

# ╔═╡ 22f61d10-185a-11eb-2c53-85f7ef0fb78a
md"""
Select which set of measurements to display:

$(@bind krsnames MultiSelect(sort(unique(map(kr->NeXLDatabase.krname(kr), allkrs)))))
"""

# ╔═╡ 2a8bd5f0-185b-11eb-0b96-8592774730a0
begin
	krs = filter(kr->NeXLDatabase.krname(kr) in krsnames, allkrs)
	if length(krs)>0
		plot(krs) |> SVG(joinpath(homedir(),"Desktop","k-ratio plot.svg"),6inch,4inch)
		plot(krs)
	else
		md"No k-ratios have been selected"
	end
end

# ╔═╡ cb991c00-1888-11eb-02fb-b33a153473dc
asa(DataFrame, DBCampaign[ read(db, DBCampaign, kr.campaign) for kr in krs])

# ╔═╡ 6185fa20-19d9-11eb-04d3-3132221ca798
begin
	if length(krs)>0
		plot(krs, style=:XY)
	else
		md"No k-ratios have been selected"
	end
end

# ╔═╡ Cell order:
# ╟─a6c0b440-5688-11eb-0635-39912417c1f9
# ╟─3742ed50-1852-11eb-01e3-5b27d6c85250
# ╟─e793f370-1852-11eb-24b6-05fa57c3a3c6
# ╟─78885a50-1859-11eb-0b6c-3368b717d267
# ╟─a51b0a40-1921-11eb-2c1c-ab31e48e7846
# ╟─22f61d10-185a-11eb-2c53-85f7ef0fb78a
# ╠═2a8bd5f0-185b-11eb-0b96-8592774730a0
# ╟─cb991c00-1888-11eb-02fb-b33a153473dc
# ╟─6185fa20-19d9-11eb-04d3-3132221ca798


text:
module OldED

using LinearAlgebra
using SparseArrays
import Base: view, copy, eltype, size, Array
import Base: +, -, *, /, sum
import LinearAlgebra: mul!

include("Basis.jl")
include("Operation.jl")
include("Operation_construct.jl")
include("Measure.jl")
include("Spin.jl")
#-----------------------------------------------------------------------------------------------------
# Krylov space
#-----------------------------------------------------------------------------------------------------
function reducespace(
    vector_space::AbstractMatrix;
    rtol::Real=1e-3
)
    mat = vector_space' * vector_space
    vals, vecs = eigen(Hermitian(mat))
    pos = vals .> rtol
    vals_pos = vals[pos]
    vecs_pos = vecs[:, pos]
    null_vector_space = vector_space * vecs_pos
    for i = 1:length(vals_pos)
        null_vector_space[:, i] ./= sqrt(vals_pos[i])
    end
    return null_vector_space
end
# Special case for single vector
reducespace(v::AbstractVector) = normalize(v)
#-----------------------------------------------------------------------------------------------------
export krylovspace
function krylovspace(
    H,
    vector_space::AbstractMatrix;
    rtol::Real=1e-3
)
    vs_rank = rank(vector_space)
    cat_space = hcat(vector_space, H * vector_space)
    new_space = reducespace(cat_space, rtol=rtol)
    new_rank = size(new_space, 2)
    while new_rank > vs_rank
        vs_rank = new_rank
        cat_space = hcat(new_space, H * new_space)
        new_space = reducespace(cat_space, rtol=rtol)
        new_rank = size(new_space, 2)
    end
    return new_space
end
# Special case for vector
krylovspace(H, v::AbstractVector; rtol::Real=1e-3) = krylovspace(H, reshape(v, :, 1), rtol=rtol)

end # module ExactDiagonalization


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    API implementation for J2 osculating orbit propagator.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

get_epoch(orbp::OrbitPropagatorJ2osc) = orbp.j2oscd.j2d.epoch

function get_mean_elements(orbp::OrbitPropagatorJ2osc)
    j2d = orbp.j2oscd.j2d
    orb = KeplerianElements(
        j2d.epoch + j2d.Δt / 86400,
        j2d.al_k * j2d.j2_gc.R0,
        j2d.e_k,
        j2d.i_k,
        j2d.Ω_k,
        j2d.ω_k,
        j2d.f_k
    )

    return orb
end

"""
    init_orbit_propagator(Val(:J2osc), epoch::Number, a_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, f_0::Number, dn_o2::Number = 0, ddn_o6::Number = 0; j2_gc::J2_GravCte{T} = j2_gc_egm08) where T
    init_orbit_propagator(Val(:J2osc), orb_0::Orbit, dn_o2::Number = 0, ddn_o6::Number = 0; j2_gc::J2_GravCte = j2_gc_egm08)

Initialize the J2 osculating orbit propagator.

# Args

- `epoch::Number`: Initial orbit epoch [Julian Day].
- `a_0::Number`: Initial mean semi-major axis [m].
- `e_0::Number`: Initial mean eccentricity.
- `i_0::Number`: Initial mean inclination [rad].
- `Ω_0::Number`: Initial mean right ascension of the ascending node [rad].
- `ω_0::Number`: Initial mean argument of perigee [rad].
- `f_0::Number`: Initial mean true anomaly [rad].
- `dn_o2::Number`: (OPTIONAL) First time derivative of mean motion divided by 2
    \\[rad/s²] (**Default** = 0).
- `ddn_o6::Number`: (OPTIONAL) Second time derivative of mean motion divided by 6
    \\[rad/s³] (**Default** = 0).
- `orb_0::Orbit`: Object of type [`Orbit`](@ref) with the initial mean orbital
    elements [SI].

# Keywords

- `j2_gc::J2_GravCte`: (OPTIONAL) J2 orbit propagator gravitational constants.
    (**Default** = `j2_gc_egm08`).
"""
function init_orbit_propagator(
    ::Val{:J2osc},
    epoch::Number,
    a_0::Number,
    e_0::Number,
    i_0::Number,
    Ω_0::Number,
    ω_0::Number,
    f_0::Number,
    dn_o2::Number = 0,
    ddn_o6::Number = 0;
    j2_gc::J2_GravCte{T} = j2_gc_egm08
) where T
    # Create the new Two Body propagator structure.
    j2oscd = j2osc_init(
        epoch,
        a_0,
        e_0,
        i_0,
        Ω_0,
        ω_0,
        f_0,
        dn_o2,
        ddn_o6;
        j2_gc = j2_gc
    )

    # Create and return the orbit propagator structure.
    return OrbitPropagatorJ2osc(j2oscd)
end

function init_orbit_propagator(
    ::Val{:J2osc},
    orb_0::Orbit,
    dn_o2::Number = 0,
    ddn_o6::Number = 0;
    j2_gc::J2_GravCte = j2_gc_egm08
)
    # Convert the orbit representation to Keplerian elements.
    k_0 = convert(KeplerianElements, orb_0)

    return init_orbit_propagator(
        Val(:J2osc),
        k_0.t,
        k_0.a,
        k_0.e,
        k_0.i,
        k_0.Ω,
        k_0.ω,
        k_0.f,
        dn_o2,
        ddn_o6;
        j2_gc = j2_gc
    )
end

function propagate!(orbp::OrbitPropagatorJ2osc, t::Number)
    # Auxiliary variables.
    j2oscd = orbp.j2oscd

    # Propagate the orbit.
    r_i, v_i = j2osc!(j2oscd, t)

    # Return.
    return r_i, v_i
end

function step!(orbp::OrbitPropagatorJ2osc, Δt::Number)
    # Auxiliary variables.
    j2oscd = orbp.j2oscd

    # Propagate the orbit.
    r_i, v_i = j2osc!(j2oscd, j2oscd.Δt + Δt)

    # Return the information about the step.
    return r_i, v_i
end


text:
using DisplayStructure
const DS = DisplayStructure
using Terming
const T = Terming

##### Views #####

function gen_views(h, w, title, str)
    display_stack =  DS.DisplayStack()
    push!(display_stack, :form => DS.Panel(title, [h, w], [1, 1]))
    push!(display_stack, :str => DS.Label(str, [5, 5]))

    return display_stack
end

##### Models #####

abstract type Model end

struct MovableStr <: Model
    reachable_size::Tuple{UnitRange{Int64}, UnitRange{Int64}}
    pos::Vector{Int}
end

function MovableStr(
    y_range::UnitRange{Int64},
    x_range::UnitRange{Int64}
)
    return MovableStr((y_range, x_range),[5, 5])
end

function move(str::MovableStr, direction::Symbol)
    if direction === :up
        new_y = str.pos[1] - 1
        (new_y in str.reachable_size[1]) && (str.pos[1] = new_y; return true)
    elseif direction === :down
        new_y = str.pos[1] + 1
        (new_y in str.reachable_size[1]) && (str.pos[1] = new_y; return true)
    elseif direction === :right
        new_x = str.pos[2] + 1
        (new_x in str.reachable_size[2]) && (str.pos[2] = new_x; return true)
    elseif direction === :left
        new_x = str.pos[2] - 1
        (new_x in str.reachable_size[2]) && (str.pos[2] = new_x; return true)
    end

    return false
end

##### Controls #####

function init_term()
    T.raw!(true)
    T.alt_screen(true)
    T.cshow(false)
    T.clear()
end

function reset_term()
    T.raw!(false)
    T.alt_screen(false)
    T.cshow(true)
end

struct App
    views::DisplayStack
    model::Model
end

function App()
    h, w = T.displaysize()
    title = " Title "
    str = "會動的字串"

    return App(
        gen_views(h, w, title, str),
        MovableStr(2:(h-1), 2:(w-1-textwidth(str)+1))
    )
end

render(app::App) = DS.render(app.views)

function update(app::App)
    app.views[:str].pos .= app.model.pos
    render(app)
end

function handle_quit(::App)
    keep_running = false
    T.cmove_line_last()
    T.println("\nShutting down")
    return keep_running
end

function handle_event(app::App)
    is_running = true
    while is_running
        sequence = T.read_stream()
        if sequence == "\e"
            is_running = handle_quit(app)
        elseif sequence == "w"
            move(app.model, :up) && update(app)
        elseif sequence == "s"
            move(app.model, :down) && update(app)
        elseif sequence == "d"
            move(app.model, :right) && update(app)
        elseif sequence == "a"
            move(app.model, :left) && update(app)
        end
    end
end

function Base.run(app::App)
    init_term()
    render(app)
    handle_event(app)
    reset_term()
    return
end

##### Main #####

run(App())


text:
include("lib.jl")
include("analyze.jl")
include("execute.jl")
include("docs.jl")

using .lib
using .analyze
using .execute
using .docs

d=Data()

argAnalyze(d)

if d.mode == CMMain
	exec(d)
elseif d.mode == CMHelp
	help()
elseif d.mode == CMVersion
	version()
end

max_stars_repo_path: Time Measurement/Julia/main.jl

max_issues_repo_path: Time Measurement/Julia/main.jl

max_forks_repo_path: Time Measurement/Julia/main.jl


text:
##
# Configurable constants
######################################################
## constants while working in local file system mode
const fs_pfx = "/tmp/cc"

const part_idx_location     = joinpath(fs_pfx, "part_idx")
const doc_to_id_location    = joinpath(fs_pfx, "doc_to_id")
const id_to_doc_location    = joinpath(fs_pfx, "id_to_doc")
const stat_location         = joinpath(fs_pfx, "stats")
#const mime_location         = joinpath(fs_pfx, "mime_counts")
#const domain_location       = joinpath(fs_pfx, "domain_counts")
#const charset_location       = joinpath(fs_pfx, "charset_counts")
const docs_location         = joinpath(fs_pfx, "docs")


text:
# ============================================================================
# exported from MatrixProcessing on 10.09.2020
"""
    get_pairwise_correlation_matrix(vectorized_video, tau_max=25)

Computes pairwise correlation of the input signals accordingly to the formula
presented in paper "Clique topology reveals intrinsic geometric structure
in neural correlations" by Chad Giusti et al.

The Computations are done only for upper half of the matrix, the lower half is
a copy of upper half. Computation-wise the difference is at level of 1e-16, but
this causes that inverse is not the same as non-inverse matrix.

"""
function get_pairwise_correlation_matrix(vectorized_video, tau_max=25)
    number_of_signals = size(vectorized_video,1)
    T = size(vectorized_video,2)

    C_ij = zeros(number_of_signals,number_of_signals);
    # log_C_ij = zeros(number_of_signals,number_of_signals);

     # this is given in frames
    lags = -tau_max:1:tau_max

    for row=1:number_of_signals
        for column=row:number_of_signals
            signal_ij = vectorized_video[row,:];
            signal_ji = vectorized_video[column,:];

            # cross_corelation
            ccg_ij = crosscov(signal_ij, signal_ji, lags);
            ccg_ij = ccg_ij ./ T;

            A = sum(ccg_ij[tau_max+1:end]);
            B = sum(ccg_ij[1:tau_max+1]);
            r_i_r_j = 1;
            C_ij[row, column] = max(A, B)/(tau_max*r_i_r_j);
            C_ij[column, row] = C_ij[row, column]
            # log_C_i_j[row, column] = log10(abs(C_ij[row, column]));
        end
    end

    return C_ij
end

"""
    get_subimg_correlations(video_array, centers, sub_img_size, shift)

Computes the correlation between the subimages and subimages shifted by values
from range -@shift:@shift and returns array with frames of size
length(@centers) x length(@centers) with the number of frames equal to the
number of rames in @video_array.

Each of the subimage is center around values stored in  @centers
"""
# TODO Check if this is the same as some of the functions from the ImageProcessing
function get_subimg_correlations(video_array, centers, sub_img_size, shift)
    half_size = ceil(Int,(sub_img_size-1)/2)
    half_range = half_size + shift
    h, w, len = get_video_dimension(video_array)
    extracted_pixels = zeros(sub_img_size, sub_img_size, len)

    for frame = 1:len
        img = video_array[frame]
        for index_x = 1:size(centers,2)
            c_x = centers[2, index_x]
            for index_y = 1:size(centers,2)
                c_y = centers[1, index_y]
                subimage = img[(c_x-half_range):(c_x+half_range),
                                (c_y-half_range):(c_y+half_range)]
                center = img[(c_x-half_size):(c_x+half_size), (c_y-half_size):(c_y+half_size)]

                for left_boundary = 1:(2*shift+1)
                    for lower_boundary = 1:(2*shift+1)
                        corelation = center .* subimage[left_boundary:left_boundary+sub_img_size-1, lower_boundary:lower_boundary+sub_img_size-1]
                        corelation = sum(corelation)
                        extracted_pixels[index_x, index_y, frame] += corelation
                    end
                end
                extracted_pixels[index_x, index_y, frame] /= 256*(sub_img_size^2)*(shift*2)^2
            end
        end
    end
    return extracted_pixels
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    TriRefinement()

Refinement of polygonal meshes into triangles.
A n-gon is subdivided into n-2 triangles.
"""
struct TriRefinement <: RefinementMethod end

function refine(mesh, method::TriRefinement)
  if eltype(mesh) <: Triangle
    # go ahead and refine
    _refine(mesh, method)
  else
    # triangulate non-triangle elements
    triangulate(mesh)
  end
end

function _refine(mesh, ::TriRefinement)
  # retrieve geometry and topology
  points = vertices(mesh)
  connec = topology(mesh)

  # convert to half-edge structure
  t = convert(HalfEdgeTopology, connec)

  # add centroids of elements
  ∂₂₀ = Boundary{2,0}(t)
  epts = map(1:nelements(t)) do elem
    ps = view(points, ∂₂₀(elem))
    cₒ = sum(coordinates, ps) / length(ps)
    Point(cₒ)
  end

  # original vertices
  vpts = points

  # new points in refined mesh
  newpoints = [vpts; epts]

  offset = length(vpts)

  # connect vertices into new triangles
  newconnec = Connectivity{Triangle,3}[]
  for elem in 1:nelements(t)
    verts = CircularVector(∂₂₀(elem))
    for i in 1:length(verts)
      u = elem + offset
      v = verts[i]
      w = verts[i+1]
      tri = connect((u, v, w))
      push!(newconnec, tri)
    end
  end

  SimpleMesh(newpoints, newconnec)
end


text:
include("builder.jl")
include("dielectric.jl")
include("dftk.jl")
include("quantum_espresso.jl")
using FFTW
using LinearAlgebra
using DataFrames
import JSON
import TimerOutputs

function setup_threading()
    FFTW.set_num_threads(4)
    BLAS.set_num_threads(4)
end

function random_logfile(prefix)
    logfile = "$(prefix)_$(@sprintf "%05d" rand(UInt16)).log"
    while isfile(logfile)
        logfile = "$(prefix)_$(@sprintf "%05d" rand(UInt16)).log"
    end
    logfile
end

function run_step_on_cases(payload, cases)
    spay = string(payload)
    header_printed = false

    # TODO Stopfile to stop iterations

    anymissing = false
    for dir in sort(collect(keys(cases)))
        config = cases[dir]

        if isfile(joinpath(config.directory, string(payload) * ".done"))
            continue
        end

        for i in 1:10  # To avoid infinite loop
            if !header_printed
                println("#####" * "#"^length(spay) * "#####")
                println("#--  $(string(payload))  --#")
                println("#####" * "#"^length(spay) * "#####")
                header_printed = true
            end

            logfile = random_logfile(joinpath(config.directory, string(payload)))
            println(basename(dir), "  ->  ", basename(dir) * "/" * basename(logfile))
            res = open(logfile, "w") do fp
                redirect_stdout(fp) do
                    payload(config)
                end
            end
            println("    ... $res\n")
            if res == :exhausted
                rm(logfile)
                anymissing || touch(joinpath(config.directory, string(payload) * ".done"))
                break
            elseif res == :missing
                anymissing = true
                rm(logfile)
            end
        end
    end
end

function firstrun(config)
    lockfile = joinpath(config.directory, "firstrun.running")
    firstprefix = joinpath(config.directory, "firstrun")
    (isfile(lockfile) || has_results(firstprefix)) && return :exhausted
    write(lockfile, read(`hostname`, String))

    lattice, atoms = load_system(config.case, config.n_repeat)
    firstα = get(config, :firstα, 0.5)
    localiser = config.construct_localiser(lattice, atoms)
    firstmixing = config.mixings[config.firstmixkey](firstα, localiser)

    try
        println("#\n#-- Firstrun $(config.case) n_repeat=$(config.n_repeat) -- " *
                "mixing=$firstmixing\n#")
        TimerOutputs.reset_timer!(DFTK.timer)
        run_dftk(lattice, atoms; fileprefix=firstprefix, mixing=firstmixing,
                 maxiter=100, config.kwargs_scf...)
        println(DFTK.timer)
    catch
        rm(joinpath(config.directory, "firstrun.done"), force=true)
        rethrow()
    finally
        rm(lockfile)
    end

    :done
end

function optimal_damping(config)
    lattice, atoms = load_system(config.case, config.n_repeat)
    localiser = config.construct_localiser(lattice, atoms)

    # Load firstrun results
    firstprefix = joinpath(config.directory, "firstrun")
    has_results(firstprefix) || return :missing
    firstmixing = config.mixings[config.firstmixkey](0.5, localiser)
    scfres = run_dftk(lattice, atoms; fileprefix=firstprefix, mixing=firstmixing,
                      config.kwargs_scf...)

    mixkey = nothing
    lockfile = nothing
    resultfile = nothing
    for mix in sort(collect(keys(config.mixings)))
        lockfile = joinpath(config.directory, mix * ".running")
        resultfile = joinpath(config.directory, mix * ".optimal_damping.json")
        (isfile(lockfile) || isfile(resultfile)) && continue
        write(lockfile, read(`hostname`, String))
        mixkey = mix
        break
    end
    mixkey === nothing && return :exhausted

    try
        mixing = config.mixings[mixkey](1.0, localiser)
        println("#\n#-- optimal_damping $(config.case) n_repeat=$(config.n_repeat) " *
                "-- mixing=$mixing\n#")
        fileprefix = joinpath(config.directory, "$(mixkey)_optimal_damping")

        TimerOutputs.reset_timer!(DFTK.timer)
        arnoldires = run_arnoldi_dielectric(scfres, mixing=mixing, force_lda=true,
                                            fileprefix=fileprefix, iterations=60)
        println(DFTK.timer)

        # ε = 1 - χ0 (vc + fxc)
        # In linear response regime: F(ρₙ) = (1 - ε) ρₙ
        # SCF Damped fixed-point scheme is
        #     ρₙ₊₁ = ρₙ + α (F(ρₙ) - ρₙ)
        #          = ρₙ + α (1 - ε - 1)ρₙ
        #          = (1 - αε) ρₙ
        # which converges iff
        #    σ(1 - αε) ⊂ (-1, 1)
        # but where things converge better the further the eigenvalues are from 1 and -1.
        # Optimising properly the problem
        #    min_α max(|1 - αε|, |1 + αε|))
        # shows the best α to give
        #    min σ(1 - αε) = max σ(1 - αε) ∈ (-1, 1)
        # or
        #    α_opt = 2 / (λ_max(ε) + λ_min(ε))

        eigenvalues  = real(arnoldires.λ)
        best_damping = 2 / (maximum(abs, eigenvalues) + minimum(abs, eigenvalues))

        # Since we are doing LDA here but our calculation is PBE there is really no
        # point in keeping more than 3 significant digits
        best_damping = round(best_damping, sigdigits=3)

        open(resultfile, "w") do fp
            save_dict = Dict(
                "eigenvalues" => real(arnoldires.λ),
                "residual_norms" => real(arnoldires.residual_norms),
                "condition_number" => maximum(abs, eigenvalues) / minimum(abs, eigenvalues),
                "optimal_damping" => best_damping,
            )
            JSON.print(fp, save_dict)
        end
    catch
        rm(joinpath(config.directory, "optimal_damping.done"), force=true)
        rethrow()
    finally
        rm(lockfile)
    end

    :done
end

function scf(config)
    lattice, atoms = load_system(config.case, config.n_repeat)
    localiser = config.construct_localiser(lattice, atoms)

    mixkey = nothing
    lockfile = nothing
    resultfile = nothing
    fileprefix = nothing
    for mix in sort(collect(keys(config.mixings)))
        lockfile = joinpath(config.directory, mix * ".running")
        fileprefix = joinpath(config.directory, mix)
        (isfile(lockfile) || has_results(fileprefix))  && continue
        write(lockfile, read(`hostname`, String))
        mixkey = mix
        break
    end
    mixkey === nothing && return :exhausted

    try
        # Load optimal damping results
        resultfile = joinpath(config.directory, mixkey * ".optimal_damping.json")
        isfile(resultfile) || return :missing
        αopt = open(JSON.parse, resultfile, "r")["optimal_damping"]

        mixing = config.mixings[mixkey](αopt, localiser)
        println("#\n#-- SCF $(config.case) n_repeat=$(config.n_repeat) -- " *
                "mixing=$mixing\n#")

        TimerOutputs.reset_timer!(DFTK.timer)
        run_dftk(lattice, atoms; fileprefix=fileprefix, mixing=mixing, config.kwargs_scf...)
        println(DFTK.timer)
    catch
        rm(joinpath(config.directory, "scf.done"), force=true)
        rethrow()
    finally
        rm(lockfile)
    end
    plot_convergence_summary(config)

    :done
end

function collect_minimal_total_energy(config::NamedTuple)
    collect_minimal_total_energy(config.directory, config.mixings)
end
function collect_minimal_total_energy(directory::AbstractString, mixings::Dict)
    # First pass: For DFTK runs collect mimimal total energy
    #       (The implementation in QE is slightly different
    #        resulting in slightly deviating energies)
    mintotal = 0
    for mixkey in sort(collect(keys(mixings)))
        resultfile = joinpath(directory, mixkey * ".json")
        isfile(resultfile) || continue

        energies = open(JSON.parse, resultfile, "r")["energies"]
        energies === nothing && continue
        isempty(energies) && continue

        mintotal = min(mintotal, minimum(energies[max(end-4, 1):end]))
    end

    mintotal
end

function plot_convergence_summary(config;
                                  iteration=30, tol=1e-10,
                                  outputfile=joinpath(config.directory, "summary.svg"))
    style = Dict(
        "None"          => (color=4, linestyle=:dot,   mark=:cross),
        "TF"            => (color=5, linestyle=:dot,   mark=:xcross),
        "Kerker"        => (color=2, linestyle=:solid, mark=:xcross),
        "RestaSiO2"     => (color=0, linestyle=:solid, mark=:xcross),
        "RestaGaAs"     => (color=1, linestyle=:dash,  mark=:xcross),
        "HybridSiO2"    => (color=0, linestyle=:solid, mark=:dtriangle),
        "HybridGaAs"    => (color=1, linestyle=:dash,  mark=:star4),
        "HybridVac"     => (color=4, linestyle=:dot,   mark=:star6),
        "Custom"        => (color=5, linestyle=:dash,  mark=:star8),
    )
    errors = Dict{String, Vector{Float64}}()

    # First pass: For DFTK runs collect mimimal total energy
    #       (The implementation in QE is slightly different
    #        resulting in slightly deviating energies)
    mintotal = collect_minimal_total_energy(config)

    # Add DFTK errors
    for mixkey in sort(collect(keys(config.mixings)))
        resultfile = joinpath(config.directory, mixkey * ".json")
        isfile(resultfile) || continue

        energies = open(JSON.parse, resultfile, "r")["energies"]
        energies === nothing && continue
        values = abs.(energies .- mintotal)
        idx = findfirst(values .< 1e-16)
        isnothing(idx) || (values = values[1:idx-1])
        isempty(values) && continue

        errors[mixkey] = values
    end

    # Add Quantum-Espresso errors
    qe_fn = joinpath(config.directory, "TF/espresso.pwo")
    if isfile(qe_fn)
        errors["TF"] = [p.accuracy for p in parse_pwo(qe_fn)]
    end

    # If no data return
    isempty(errors) && return plot()

    ymax = 0
    maxiter = 0
    p = plot(legend=:topright, title=config.case)
    for mixkey in sort(collect(keys(errors)))
        dampfile = joinpath(config.directory, mixkey * ".optimal_damping.json")
        αopt = isfile(dampfile) ? open(JSON.parse, dampfile, "r")["optimal_damping"] : NaN
        mixkey == "TF" && (αopt = 0.7)  # QE calculations all use α=0.7

        maxiter = max(maxiter, length(errors[mixkey]))
        ymax = max(ymax, maximum(errors[mixkey]))
        plot!(p, errors[mixkey]; label=mixkey * " (α=$αopt)", style[mixkey]...)
    end
    xlims!(p, (0, min(iteration, maxiter)))
    ylims!(p, (tol, 2 * ymax))
    xaxis!(p, "Iteration")
    yaxis!(p, "Total energy absolute error", :log10)

    if isnothing(outputfile)
        p
    else
        savefig(p, outputfile)
    end
end

function dump_ldos(config, temperatures=[1e-3, 2e-3, 4e-3, 8e-3, 16e-3, 32e-3, 64e-3, 128e-3])
    lattice, atoms = load_system(config.case, config.n_repeat)
    localiser = config.construct_localiser(lattice, atoms)

    # Load firstrun results
    firstprefix = joinpath(config.directory, "firstrun")
    has_results(firstprefix) || return nothing
    firstmixing = config.mixings[config.firstmixkey](0.5, localiser)
    scfres = run_dftk(lattice, atoms; fileprefix=firstprefix, mixing=firstmixing,
                      config.kwargs_scf...)

    isnothing(temperatures) && (temperatures = [scfres.ham.basis.model.temperature])
    h5open(joinpath(config.directory, "ldos.hdf5"), "w") do h5f
        h5f["temperatures"] = temperatures
        h5f["temperature_scf"] = scfres.ham.basis.model.temperature
        h5f["ρ_real", "shuffle", (), "compress", 3] = real(scfres.ρ.real)
        for temp in temperatures
            stemp = @sprintf "%.3f" temp
            ldos = LDOS(scfres.εF, scfres.ham.basis, scfres.eigenvalues, scfres.ψ,
                        temperature=temp)
            h5f["ldos_$stemp", "shuffle", (), "compress", 3] = ldos
        end
    end
end

function collect_results(config; tol=1.6e-10)
    df = DataFrame(case=String[],
                   repeat=Int[],
                   size=Float64[],
                   mixing=String[],
                   damping=Union{Float64,Missing}[],
                   iterations=Union{Missing,Int}[],
                   condition_number=Union{Float64,Missing}[],
                   λmin=Union{Float64,Missing}[],
                   λmax=Union{Float64,Missing}[],
                   resmin=Union{Float64,Missing}[],
                   resmax=Union{Float64,Missing}[],
                  )

    for cfgkey in sort(collect(keys(config)))
        cfg        = config[cfgkey]
        minenergy  = collect_minimal_total_energy(cfg)
        lattice, _ = load_system(cfg.case, cfg.n_repeat)

        # Some made-up size function for now:
        size = norm(lattice * [0.0, 0.0, 1.0])

        for mixkey in sort(collect(keys(cfg.mixings)))
            scfdata  = Dict{String, Any}()
            dampdata = Dict{String, Any}()
            dampfile = joinpath(cfg.directory, mixkey * ".optimal_damping.json")
            scffile  = joinpath(cfg.directory, mixkey * ".json")

            isfile(dampfile) && (dampdata = open(JSON.parse, dampfile, "r"))
            isfile(scffile) && (scfdata  = open(JSON.parse, scffile, "r"))

            iterations = missing
            if "energies" in keys(scfdata)
                itfirst = findfirst(abs.(scfdata["energies"] .- minenergy) .≤ tol)
                iterations = isnothing(itfirst) ? typemax(Int) : itfirst
            end
            push!(df, (cfg.case,
                       cfg.n_repeat,
                       size,
                       mixkey,
                       get(dampdata, "optimal_damping", missing),
                       iterations,
                       get(dampdata, "condition_number", missing),
                       get(dampdata, "eigenvalues", [missing])[1],
                       get(dampdata, "eigenvalues", [missing])[end],
                       get(dampdata, "residual_norms", [missing])[1],
                       get(dampdata, "residual_norms", [missing])[end],
                 ))
        end
    end

    df
end


text:
# tests

# data generation
srand(0)
Np = 50
θ0 = randn(Np)
θ0[ rand(Np) .> 0.1 ] = 0.

Nobs = 10000
x = randn(Np, Nobs)
u = vec(θ0' * x) + 0.1randn(Nobs)
y = ( 1 ./ (1. + exp(-u)) .> 0.5 ) + 0.

# convergence comparisons
model = m.riskmodel(m.LinearPred(Np), m.LogisticLoss())

θ1 = m.fbs(model, zeros(Np), x, y,
           reg=m.SqrL2Reg(10.0), λ=2e-3, maxsteps=100, maxtime=10)
θ2 = m.fista(model, zeros(Np), x, y,
             reg=m.SqrL2Reg(10.0), λ=2e-3, maxsteps=100, maxtime=10)
θ3 = m.fasta(model, zeros(Np), x, y,
             reg=m.SqrL2Reg(10.0), λ=1e-3, maxsteps=100, maxtime=10)

@test sumabs2(θ1-θ2) < 1e-2
@test sumabs2(θ2-θ3) < 1e-2

θ1 = m.fbs(model, zeros(Np), x, y,
           reg=m.L1Reg(10.0), λ=2e-3, maxsteps=1000, maxtime=10, cbinterval=1000)
θ2 = m.fista(model, zeros(Np), x, y,
             reg=m.L1Reg(10.0), λ=2e-3, maxsteps=1000, maxtime=10, cbinterval=1000)
θ3 = m.fasta(model, zeros(Np), x, y,
             reg=m.L1Reg(10.0), λ=1e-3, maxsteps=100, maxtime=10, cbinterval=100)

@test sumabs2(θ1-θ2) < 1e-2
@test sumabs2(θ2-θ3) < 1e-2

# fbs!(riskmodel(LinearPred(Np), MyLoss()), zeros(Np), x, y,
#      reg=SqrL2Reg(1.0), λ=1e-3, maxsteps=20, maxtime=10, cbinterval = 2)

# fista!(riskmodel(LinearPred(Np), MyLoss()), zeros(Np), x, y,
#      reg=SqrL2Reg(1.0), λ=1e-3, maxsteps=20, maxtime=10, cbinterval = 2)

# fasta!(riskmodel(LinearPred(Np), MyLoss()), zeros(Np), x, y,
#      reg=SqrL2Reg(1.0), λ=1e-3, maxsteps=20, maxtime=10, cbinterval = 1)
# fasta2!(riskmodel(LinearPred(Np), MyLoss()), zeros(Np), x, y,
#      reg=SqrL2Reg(1.0), λ=1e-3, maxsteps=20, maxtime=10, cbinterval = 1)

# fbs!(riskmodel(LinearPred(Np), MyLoss()), zeros(Np), x, y,
#      reg=L1Reg(10.0), λ=1e-3, maxsteps=20, maxtime=10, cbinterval = 2)

# fista!(riskmodel(LinearPred(Np), MyLoss()), zeros(Np), x, y,
#      reg=L1Reg(10.0), λ=1e-3, maxsteps=20, maxtime=10, cbinterval = 2)

# fasta!(riskmodel(LinearPred(Np), MyLoss()), zeros(Np), x, y,
#      reg=L1Reg(10.0), λ=1e-3, maxsteps=20, maxtime=10, cbinterval = 1)
# fasta2!(riskmodel(LinearPred(Np), MyLoss()), zeros(Np), x, y,
#      reg=L1Reg(10.0), λ=1e-3, maxsteps=20, maxtime=10, cbinterval = 1)

# fasta!(riskmodel(LinearPred(Np), MyLoss()), zeros(Np), x, y,
#      reg=L1Reg(10.0), λ=10., maxsteps=50, maxtime=10, cbinterval = 2)
# fasta2!(riskmodel(LinearPred(Np), MyLoss()), ones(Np), x, y,
#      reg=L1Reg(10.0), λ=10., maxsteps=50, maxtime=10, cbinterval = 2)
# fasta2!(riskmodel(LinearPred(Np), MyLoss()), ones(Np), x, y,
#      reg=L1Reg(1.0), λ=10., maxsteps=50, maxtime=10, cbinterval = 2)


text:
include("pedigree_run.jl")
include("pedigree_utils.jl")

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
######################### STaylor1

abstract type AbstractSeries{T<:Number} <: Number end

"""
    STaylor1{N,T<:Number} <: AbstractSeries{T}

DataType for polynomial expansions in one independent variable.

**Fields:**

- `coeffs :: NTuple{N,T}` Expansion coefficients; the ``i``-th
    component is the coefficient of degree ``i-1`` of the expansion.

Note that `STaylor1` variables are callable. For more information, see
[`evaluate`](@ref).
"""
struct STaylor1{N,T<:Number} <: AbstractSeries{T}
    coeffs::NTuple{N,T}
    function STaylor1{N,T}(coeffs::NTuple{N,T}) where {N, T <: Number}
        new(coeffs)
    end
end
function STaylor1(coeffs::NTuple{N,T}) where {N, T <: Number}
    STaylor1{N,T}(coeffs)
end

## Outer constructors ##

"""
    STaylor1(x::T, v::Val{N})

Shortcut to define the independent variable of a `STaylor1{N,T}` polynomial of
given `N` with constant term equal to `x`.
"""
@generated function STaylor1(x::T, v::Val{N}) where {N,T<:Number}
    y = Any[:(zero($T)) for i=1:N]
    tup = :((x,))
    push!(tup.args, y...)
    return quote
        Base.@_inline_meta
        STaylor1{(N+1),T}($tup)
    end
end
function STaylor1(coeffs::Vector{T}, l::Val{L}, v::Val{N}) where {L,N,T<:Number}
    STaylor1{(N+1),T}(ntuple(i -> (i < L+1) ? coeffs[i] : zero(T),  N+1))
end
@inline function STaylor1(coeffs::Vector{T}) where {T<:Number}
    STaylor1{length(coeffs),T}(tuple(coeffs...))
end
@inline STaylor1(x::STaylor1{N,T}) where {N,T<:Number} = x


text:
module Tools

include("blast/BLAST.jl")

end  # module Tools


text:
using GLPKMathProgInterface

const mathprogbase_test = joinpath(Pkg.dir("MathProgBase"), "test")

include(joinpath(mathprogbase_test, "linprog.jl"))
linprogtest(GLPKSolverLP())

include(joinpath(mathprogbase_test, "linproginterface.jl"))
linprogsolvertest(GLPKSolverLP())

include(joinpath(mathprogbase_test, "mixintprog.jl"))
mixintprogtest(GLPKSolverMIP())

include(joinpath(mathprogbase_test,"conicinterface.jl"))
coniclineartest(GLPKSolverLP())
coniclineartest(GLPKSolverMIP())

lpm = GLPKInterfaceLP.LinearQuadraticModel(GLPKSolverLP(it_lim=9513, tol_bnd=4.149))
@test lpm.param.it_lim == 9513
@test lpm.param.tol_bnd == 4.149
mipm = GLPKInterfaceMIP.LinearQuadraticModel(GLPKSolverMIP(it_lim=5910, tol_obj=1.52e-3))
@test mipm.smplxparam.it_lim == 5910
@test mipm.param.tol_obj == 1.52e-3


text:
### Data structures ###
abstract type InputData end
abstract type TSData <:InputData end
abstract type OptData <: InputData end
abstract type ClustResult end

"FullInputData"
struct FullInputData <: TSData
 region::String
 years::Array{Int64}
 N::Int
 data::Dict{String,Array}
end

"ClustData \n weights: this is the absolute weight. E.g. for a year of 365 days, sum(weights)=365"
struct ClustData <: TSData
 region::String
 years::Array{Int64}
 K::Int
 T::Int
 data::Dict{String,Array}
 weights::Array{Float64}
 mean::Dict{String,Array}
 sdv::Dict{String,Array}
 deltas::Array{Float64,2}
 k_ids::Array{Int64}
end

"ClustDataMerged"
struct ClustDataMerged <: TSData
 region::String
 years::Array{Int64}
 K::Int
 T::Int
 data::Array
 data_type::Array{String}
 weights::Array{Float64}
 mean::Dict{String,Array}
 sdv::Dict{String,Array}
 deltas::Array{Float64,2}
 k_ids::Array{Int64}
end

"ClustResultAll"
struct ClustResultAll <: ClustResult
 best_results::ClustData
 best_ids::Array{Int,1}
 best_cost::Float64
 data_type::Array{String}
 clust_config::Dict{String,Any}
 centers::Array{Array{Float64},1}
 weights::Array{Array{Float64},1}
 clustids::Array{Array{Int,1},1}
 cost::Array{Float64,1}
 iter::Array{Int,1}
end

# TODO: not used yet, but maybe best to implement this one later for users who just want to use clustering but do not care about the locally converged solutions
"ClustResultBest"
struct ClustResultBest <: ClustResult
 best_results::ClustData
 best_ids::Array{Int,1}
 best_cost::Float64
 data_type::Array{String}
 clust_config::Dict{String,Any}
end

"ClustResultSimple"
struct ClustResultSimple <: ClustResult
 best_results::ClustData
 #TODO: clust_data::ClustData
 clust_config::Dict{String,Any}
end

"SimpleExtremeValueDescr"
struct SimpleExtremeValueDescr
   data_type::String
   extremum::String
   peak_def::String
   periods::Int64
   "Replace default constructor to only allow certain entries"
   function SimpleExtremeValueDescr(data_type::String,
                                    extremum::String,
                                    peak_def::String,
                                    periods::Int64)
       # only allow certain entries
       if !(extremum in ["min","max"])
         @error("extremum - "*extremum*" - not defined")
       elseif !(peak_def in ["absolute","integral"])
         @error("peak_def - "*peak_def*" - not defined")
       end
       new(data_type,extremum,peak_def,periods)
   end
end

"""
    SimpleExtremeValueDescr(data_type::String, extremum::String, peak_def::String)
"""
function SimpleExtremeValueDescr(data_type::String,
                                 extremum::String,
                                 peak_def::String)
   return SimpleExtremeValueDescr(data_type, extremum, peak_def, 1)
end

#### Constructors for data structures###

"""
    FullInputData(region::String,
                        N::Int;
                        el_price::Array=[],
                        el_demand::Array=[],
                        solar::Array=[],
                        wind::Array=[]
                        )
Constructor for FullInputData with optional data input
"""
function FullInputData(region::String,
                      N::Int;
                      el_price::Array=[],
                      el_demand::Array=[],
                      solar::Array=[],
                      wind::Array=[]
                      )
 dt = Dict{String,Array}()
 !isempty(el_price) && (dt["el_price"]=el_price)
 !isempty(el_demand) &&  (dt["el_demand"]=el_demand)
 !isempty(wind) && (dt["wind"]=wind)
 !isempty(solar) && (dt["solar"]=solar)
 # TODO: Check dimensionality of N and supplied input data streams Nx1
 isempty(dt) && @error("Need to provide at least one input data stream")
 FullInputData(region,N,dt)
end

"""
  ClustData(region::String,
                         years::Array{Int64,1},
                         K::Int,
                         T::Int;
                         el_price::Array=[],
                         el_demand::Array=[],
                         solar::Array=[],
                         wind::Array=[],
                         weights::Array{Float64}=ones(K),
                         mean::Dict{String,Array}=Dict{String,Array}(),
                         sdv::Dict{String,Array}=Dict{String,Array}(),
                         deltas::Array{Float64,2}=ones(T,K),
                         k_ids::Array{Int64,1}=collect(1:K)
                         )
constructor 1 for ClustData: provide data individually
"""
function ClustData(region::String,
                         years::Array{Int64,1},
                         K::Int,
                         T::Int;
                         el_price::Array=[],
                         el_demand::Array=[],
                         solar::Array=[],
                         wind::Array=[],
                         weights::Array{Float64}=ones(K),
                         mean::Dict{String,Array}=Dict{String,Array}(),
                         sdv::Dict{String,Array}=Dict{String,Array}(),
                         deltas::Array{Float64,2}=ones(T,K),
                         k_ids::Array{Int64,1}=collect(1:K)
                         )
   dt = Dict{String,Array}()
   mean_sdv_provided = ( !isempty(mean) && !isempty(sdv))
   if !isempty(el_price)
     dt["el_price"]=el_price
     if !mean_sdv_provided
       mean["el_price"]=zeros(T)
       sdv["el_price"]=ones(T)
     end
   end
   if !isempty(el_demand)
     dt["el_demand"]=el_demand
     if !mean_sdv_provided
       mean["el_demand"]=zeros(T)
       sdv["el_demand"]=ones(T)
     end
   end
   if !isempty(wind)
     dt["wind"]=wind
     if !mean_sdv_provided
       mean["wind"]=zeros(T)
       sdv["wind"]=ones(T)
     end
   end
   if !isempty(solar)
     dt["solar"]=solar
     if !mean_sdv_provided
       mean["solar"]=zeros(T)
       sdv["solar"]=ones(T)
     end
   end
   isempty(dt) && @error("Need to provide at least one input data stream")
   # TODO: Check dimensionality of K T and supplied input data streams KxT
   ClustData(region,years,K,T,dt,weights,mean,sdv,deltas,k_ids)
end

"""
    ClustData(region::String,
                      years::Array{Int64,1},
                      K::Int,
                      T::Int,
                      data::Dict{String,Array},
                      weights::Array{Float64},
                      deltas::Array{Float64,2},
                      k_ids::Array{Int64,1};
                      mean::Dict{String,Array}=Dict{String,Array}(),
                      sdv::Dict{String,Array}=Dict{String,Array}()
                      )
constructor 2 for ClustData: provide data as dict
"""
function ClustData(region::String,
                       years::Array{Int64,1},
                       K::Int,
                       T::Int,
                       data::Dict{String,Array},
                       weights::Array{Float64},
                       deltas::Array{Float64,2},
                       k_ids::Array{Int64,1};
                       mean::Dict{String,Array}=Dict{String,Array}(),
                       sdv::Dict{String,Array}=Dict{String,Array}()
                       )
 isempty(data) && @error("Need to provide at least one input data stream")
 mean_sdv_provided = ( !isempty(mean) && !isempty(sdv))
 if !mean_sdv_provided
   for (k,v) in data
     mean[k]=zeros(T)
     sdv[k]=ones(T)
   end
 end
 # TODO check if right keywords are used
 ClustData(region,years,K,T,data,weights,mean,sdv,deltas,k_ids)
end

"""
    ClustData(data::ClustDataMerged)
constructor 3: Convert ClustDataMerged to ClustData
"""
function ClustData(data::ClustDataMerged)
 data_dict=Dict{String,Array}()
 i=0
 for k in data.data_type
   i+=1
   data_dict[k] = data.data[(1+data.T*(i-1)):(data.T*i),:]
 end
 ClustData(data.region,data.years,data.K,data.T,data_dict,data.weights,data.mean,data.sdv,data.deltas,data.k_ids)
end

"""
    ClustData(data::FullInputData,K,T)
constructor 4: Convert FullInputData to ClustData
"""
function ClustData(data::FullInputData,
                                 K::Int,
                                 T::Int)
  data_reshape = Dict{String,Array}()
  for (k,v) in data.data
     data_reshape[k] =  reshape(v,T,K)
  end
  return ClustData(data.region,data.years,K,T,data_reshape,ones(K),ones(T,K),collect(1:K))
end

"""
    ClustDataMerged(region::String,
                        years::Array{Int64,1},
                        K::Int,
                        T::Int,
                        data::Array,
                        data_type::Array{String},
                        weights::Array{Float64},
                        k_ids::Array{Int64,1};
                        mean::Dict{String,Array}=Dict{String,Array}(),
                        sdv::Dict{String,Array}=Dict{String,Array}()
                        )
constructor 1: construct ClustDataMerged
"""
function ClustDataMerged(region::String,
                       years::Array{Int64,1},
                       K::Int,
                       T::Int,
                       data::Array,
                       data_type::Array{String},
                       weights::Array{Float64},
                       k_ids::Array{Int64,1};
                       deltas::Array{Float64}=ones(T,K),
                       mean::Dict{String,Array}=Dict{String,Array}(),
                       sdv::Dict{String,Array}=Dict{String,Array}()
                       )
 mean_sdv_provided = ( !isempty(mean) && !isempty(sdv))
 if !mean_sdv_provided
   for dt in data_type
     mean[dt]=zeros(T)
     sdv[dt]=ones(T)
   end
 end
 ClustDataMerged(region,years,K,T,data,data_type,weights,mean,sdv,deltas,k_ids)
end

"""
    ClustDataMerged(data::ClustData)
constructor 2: convert ClustData into merged format
"""
function ClustDataMerged(data::ClustData)
 n_datasets = length(keys(data.data))
 data_merged= zeros(data.T*n_datasets,data.K)
 data_type=String[]
 i=0
 for (k,v) in data.data
   i+=1
   data_merged[(1+data.T*(i-1)):(data.T*i),:] = v
   push!(data_type,k)
 end
 if maximum(data.deltas)!=1
   throw(@error "You cannot recluster data with different Δt")
 end
 ClustDataMerged(data.region,data.years,data.K,data.T,data_merged,data_type,data.weights,data.mean,data.sdv,data.deltas,data.k_ids)
end


text:
using ASTInterface
using Test

@testset "ASTInterface.jl" begin
    # Write your tests here.
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(x^3*arctanh(a+b*x)^2, x) == :(1//4*x^4*arctanh(a+b*x)^2+1//12*b^-4*(a+b*x)^2+1//12*b^-4*log(1+-1*(a+b*x)^2)+a*b^-4*arctanh(a+b*x)+-1*a*x*b^-3+-1//4*b^-4*arctanh(a+b*x)^2*(1+a^4+6*a^2)+1//4*b^-4*(1+6*a^2)*log(1+-1*(a+b*x)^2)+1//6*b^-4*(a+b*x)^3*arctanh(a+b*x)+a*b^-4*(1+a^2)*Polylogarithms.polylog(2,-1*(1+-1a+-1*b*x)^-1*(1+a+b*x))+(1/2)*b^-4*(1+6*a^2)*(a+b*x)*arctanh(a+b*x)+-1*a*b^-4*(a+b*x)^2*arctanh(a+b*x)+-1*a*b^-4*arctanh(a+b*x)^2*(1+a^2)+2*a*b^-4*(1+a^2)*arctanh(a+b*x)*log(2*(1+-1a+-1*b*x)^-1))
@test integrate(x^2*arctanh(a+b*x)^2, x) == :(-1//3*b^-3*arctanh(a+b*x)+1//3*x*b^-2+1//3*x^3*arctanh(a+b*x)^2+-1*a*b^-3*log(1+-1*(a+b*x)^2)+-1//3*b^-3*(1+3*a^2)*Polylogarithms.polylog(2,-1*(1+-1a+-1*b*x)^-1*(1+a+b*x))+1//3*b^-3*(a+b*x)^2*arctanh(a+b*x)+1//3*b^-3*arctanh(a+b*x)^2*(1+3*a^2)+-2*a*b^-3*(a+b*x)*arctanh(a+b*x)+-1//3*b^-3*(2+6*a^2)*arctanh(a+b*x)*log(2*(1+-1a+-1*b*x)^-1)+1//3*a*b^-3*arctanh(a+b*x)^2*(3+a^2))
@test integrate(x*arctanh(a+b*x)^2, x) == :((1/2)*b^-2*log(1+-1*(a+b*x)^2)+(1/2)*x^2*arctanh(a+b*x)^2+a*b^-2*Polylogarithms.polylog(2,-1*(1+-1a+-1*b*x)^-1*(1+a+b*x))+b^-2*(a+b*x)*arctanh(a+b*x)+-1*a*b^-2*arctanh(a+b*x)^2+-1//2*b^-2*arctanh(a+b*x)^2*(1+a^2)+2*a*b^-2*arctanh(a+b*x)*log(2*(1+-1a+-1*b*x)^-1))
@test integrate(arctanh(a+b*x)^2, x) == :(b^-1*arctanh(a+b*x)^2+-1*b^-1*Polylogarithms.polylog(2,-1*(1+-1a+-1*b*x)^-1*(1+a+b*x))+b^-1*arctanh(a+b*x)^2*(a+b*x)+-2*b^-1*arctanh(a+b*x)*log(2*(1+-1a+-1*b*x)^-1))
@test integrate(x^-1*arctanh(a+b*x)^2, x) == :((1/2)*Polylogarithms.polylog(3,1+-2*(1+a+b*x)^-1)+-1//2*Polylogarithms.polylog(3,1+-2*b*x*(1+-1a)^-1*(1+a+b*x)^-1)+arctanh(a+b*x)^2*log(2*b*x*(1+-1a)^-1*(1+a+b*x)^-1)+Polylogarithms.polylog(2,1+-2*(1+a+b*x)^-1)*arctanh(a+b*x)+-1*arctanh(a+b*x)^2*log(2*(1+a+b*x)^-1)+-1*Polylogarithms.polylog(2,1+-2*b*x*(1+-1a)^-1*(1+a+b*x)^-1)*arctanh(a+b*x))
@test integrate(x^-2*arctanh(a+b*x)^2, x) == :(-1*x^-1*arctanh(a+b*x)^2+b*(1+-1*a^2)^-1*Polylogarithms.polylog(2,1+-2*(1+a+b*x)^-1)+b*(2+-2a)^-1*Polylogarithms.polylog(2,-1*(1+-1a+-1*b*x)^-1*(1+a+b*x))+-1*b*(1+-1*a^2)^-1*Polylogarithms.polylog(2,1+-2*b*x*(1+-1a)^-1*(1+a+b*x)^-1)+-1*b*(2+2a)^-1*Polylogarithms.polylog(2,1+-2*(1+a+b*x)^-1)+b*(1+a)^-1*arctanh(a+b*x)*log(2*(1+a+b*x)^-1)+b*(1+-1a)^-1*arctanh(a+b*x)*log(2*(1+-1a+-1*b*x)^-1)+-2*b*(1+-1*a^2)^-1*arctanh(a+b*x)*log(2*(1+a+b*x)^-1)+2*b*(1+-1*a^2)^-1*arctanh(a+b*x)*log(2*b*x*(1+-1a)^-1*(1+a+b*x)^-1))
@test integrate(x^-3*arctanh(a+b*x)^2, x) == :(-1//2*x^-2*arctanh(a+b*x)^2+b^2*(1+-1*a^2)^-2*log(x)+1//4*b^2*(1+a)^-2*Polylogarithms.polylog(2,1+-2*(1+a+b*x)^-1)+1//4*b^2*(1+-1a)^-2*Polylogarithms.polylog(2,-1*(1+-1a+-1*b*x)^-1*(1+a+b*x))+a*b^2*(1+-1*a^2)^-2*Polylogarithms.polylog(2,1+-2*(1+a+b*x)^-1)+(1/2)*b^2*(1+-1a)^-2*arctanh(a+b*x)*log(2*(1+-1a+-1*b*x)^-1)+-1*a*b^2*(1+-1*a^2)^-2*Polylogarithms.polylog(2,1+-2*b*x*(1+-1a)^-1*(1+a+b*x)^-1)+-1*b*x^-1*(1+-1*a^2)^-1*arctanh(a+b*x)+-1*b^2*(1+a)^-2*(2+-2a)^-1*log(1+a+b*x)+-1//2*b^2*(1+a)^-1*(1+-1a)^-2*log(1+-1a+-1*b*x)+-1//2*b^2*(1+a)^-2*arctanh(a+b*x)*log(2*(1+a+b*x)^-1)+-2*a*b^2*(1+-1*a^2)^-2*arctanh(a+b*x)*log(2*(1+a+b*x)^-1)+2*a*b^2*(1+-1*a^2)^-2*arctanh(a+b*x)*log(2*b*x*(1+-1a)^-1*(1+a+b*x)^-1))
@test integrate(x^-1*arctanh(1+b*x)^2, x) == :((1/2)*Polylogarithms.polylog(3,1+2*b^-1*x^-1)+-1*arctanh(1+b*x)^2*log(-2*b^-1*x^-1)+-1*Polylogarithms.polylog(2,1+2*b^-1*x^-1)*arctanh(1+b*x))
@test integrate((c*e+d*e*x)^3*(a+b*arctanh(c+d*x)), x) == :(1//4*b*x*e^3+-1//4*b*d^-1*e^3*arctanh(c+d*x)+1//4*d^-1*e^3*(c+d*x)^4*(a+b*arctanh(c+d*x))+1//12*b*d^-1*e^3*(c+d*x)^3)
@test integrate((c*e+d*e*x)^2*(a+b*arctanh(c+d*x)), x) == :(1//3*d^-1*e^2*(c+d*x)^3*(a+b*arctanh(c+d*x))+1//6*b*d^-1*e^2*(c+d*x)^2+1//6*b*d^-1*e^2*log(1+-1*(c+d*x)^2))
@test integrate((a+b*arctanh(c+d*x))*(c*e+d*e*x), x) == :((1/2)*b*e*x+(1/2)*e*d^-1*(c+d*x)^2*(a+b*arctanh(c+d*x))+-1//2*b*e*d^-1*arctanh(c+d*x))
@test integrate((c*e+d*e*x)^-1*(a+b*arctanh(c+d*x)), x) == :(a*d^-1*e^-1*log(c+d*x)+(1/2)*b*d^-1*e^-1*Polylogarithms.polylog(2,c+d*x)+-1//2*b*d^-1*e^-1*Polylogarithms.polylog(2,-1c+-1*d*x))
@test integrate((c*e+d*e*x)^-2*(a+b*arctanh(c+d*x)), x) == :(b*d^-1*e^-2*log(c+d*x)+-1*d^-1*e^-2*(c+d*x)^-1*(a+b*arctanh(c+d*x))+-1//2*b*d^-1*e^-2*log(1+-1*(c+d*x)^2))
@test integrate((c*e+d*e*x)^-3*(a+b*arctanh(c+d*x)), x) == :((1/2)*b*d^-1*e^-3*arctanh(c+d*x)+-1//2*b*d^-1*e^-3*(c+d*x)^-1+-1//2*d^-1*e^-3*(c+d*x)^-2*(a+b*arctanh(c+d*x)))
@test integrate((a+b*arctanh(c+d*x))^2*(c*e+d*e*x)^3, x) == :(-1//4*d^-1*e^3*(a+b*arctanh(c+d*x))^2+(1/2)*a*b*x*e^3+1//3*b^2*d^-1*e^3*log(1+-1*(c+d*x)^2)+1//4*d^-1*e^3*(a+b*arctanh(c+d*x))^2*(c+d*x)^4+1//12*b^2*d^-1*e^3*(c+d*x)^2+(1/2)*b^2*d^-1*e^3*(c+d*x)*arctanh(c+d*x)+1//6*b*d^-1*e^3*(c+d*x)^3*(a+b*arctanh(c+d*x)))
@test integrate((a+b*arctanh(c+d*x))^2*(c*e+d*e*x)^2, x) == :(1//3*x*b^2*e^2+1//3*d^-1*e^2*(a+b*arctanh(c+d*x))^2+-1//3*b^2*d^-1*e^2*Polylogarithms.polylog(2,-1*(1+-1c+-1*d*x)^-1*(1+c+d*x))+-1//3*b^2*d^-1*e^2*arctanh(c+d*x)+1//3*d^-1*e^2*(a+b*arctanh(c+d*x))^2*(c+d*x)^3+-2//3*b*d^-1*e^2*(a+b*arctanh(c+d*x))*log(2*(1+-1c+-1*d*x)^-1)+1//3*b*d^-1*e^2*(c+d*x)^2*(a+b*arctanh(c+d*x)))
@test integrate((a+b*arctanh(c+d*x))^2*(c*e+d*e*x), x) == :(-1//2*e*d^-1*(a+b*arctanh(c+d*x))^2+a*b*e*x+(1/2)*e*b^2*d^-1*log(1+-1*(c+d*x)^2)+(1/2)*e*d^-1*(a+b*arctanh(c+d*x))^2*(c+d*x)^2+e*b^2*d^-1*(c+d*x)*arctanh(c+d*x))
@test integrate((a+b*arctanh(c+d*x))^2*(c*e+d*e*x)^-1, x) == :((1/2)*b^2*d^-1*e^-1*Polylogarithms.polylog(3,1+-2*(1+-1c+-1*d*x)^-1)+2*d^-1*e^-1*(a+b*arctanh(c+d*x))^2*arctanh(1+-2*(1+-1c+-1*d*x)^-1)+-1//2*b^2*d^-1*e^-1*Polylogarithms.polylog(3,-1+2*(1+-1c+-1*d*x)^-1)+b*d^-1*e^-1*(a+b*arctanh(c+d*x))*Polylogarithms.polylog(2,-1+2*(1+-1c+-1*d*x)^-1)+-1*b*d^-1*e^-1*(a+b*arctanh(c+d*x))*Polylogarithms.polylog(2,1+-2*(1+-1c+-1*d*x)^-1))
@test integrate((a+b*arctanh(c+d*x))^2*(c*e+d*e*x)^-2, x) == :(d^-1*e^-2*(a+b*arctanh(c+d*x))^2+-1*b^2*d^-1*e^-2*Polylogarithms.polylog(2,-1+2*(1+c+d*x)^-1)+-1*d^-1*e^-2*(a+b*arctanh(c+d*x))^2*(c+d*x)^-1+2*b*d^-1*e^-2*(a+b*arctanh(c+d*x))*log(2+-2*(1+c+d*x)^-1))
@test integrate((a+b*arctanh(c+d*x))^2*(c*e+d*e*x)^-3, x) == :((1/2)*d^-1*e^-3*(a+b*arctanh(c+d*x))^2+b^2*d^-1*e^-3*log(c+d*x)+-1//2*b^2*d^-1*e^-3*log(1+-1*(c+d*x)^2)+-1//2*d^-1*e^-3*(a+b*arctanh(c+d*x))^2*(c+d*x)^-2+-1*b*d^-1*e^-3*(c+d*x)^-1*(a+b*arctanh(c+d*x)))
@test integrate((a+b*arctanh(c+d*x))^2*(c*e+d*e*x)^-4, x) == :(1//3*d^-1*e^-4*(a+b*arctanh(c+d*x))^2+-1//3*b^2*d^-1*e^-4*(c+d*x)^-1+-1//3*b^2*d^-1*e^-4*Polylogarithms.polylog(2,-1+2*(1+c+d*x)^-1)+-1//3*d^-1*e^-4*(a+b*arctanh(c+d*x))^2*(c+d*x)^-3+1//3*b^2*d^-1*e^-4*arctanh(c+d*x)+-1//3*b*d^-1*e^-4*(c+d*x)^-2*(a+b*arctanh(c+d*x))+2//3*b*d^-1*e^-4*(a+b*arctanh(c+d*x))*log(2+-2*(1+c+d*x)^-1))
@test integrate((a+b*arctanh(c+d*x))^2*(c*e+d*e*x)^-5, x) == :(1//4*d^-1*e^-5*(a+b*arctanh(c+d*x))^2+-1//3*b^2*d^-1*e^-5*log(1+-1*(c+d*x)^2)+-1//4*d^-1*e^-5*(a+b*arctanh(c+d*x))^2*(c+d*x)^-4+-1//12*b^2*d^-1*e^-5*(c+d*x)^-2+2//3*b^2*d^-1*e^-5*log(c+d*x)+-1//2*b*d^-1*e^-5*(c+d*x)^-1*(a+b*arctanh(c+d*x))+-1//6*b*d^-1*e^-5*(c+d*x)^-3*(a+b*arctanh(c+d*x)))
@test integrate((a+b*arctanh(c+d*x))^3*(c*e+d*e*x)^2, x) == :(1//3*d^-1*e^2*(a+b*arctanh(c+d*x))^3+a*x*b^2*e^2+(1/2)*b^3*d^-1*e^2*Polylogarithms.polylog(3,1+-2*(1+-1c+-1*d*x)^-1)+(1/2)*b^3*d^-1*e^2*log(1+-1*(c+d*x)^2)+-1//2*b*d^-1*e^2*(a+b*arctanh(c+d*x))^2+1//3*d^-1*e^2*(a+b*arctanh(c+d*x))^3*(c+d*x)^3+b^3*d^-1*e^2*(c+d*x)*arctanh(c+d*x)+(1/2)*b*d^-1*e^2*(a+b*arctanh(c+d*x))^2*(c+d*x)^2+-1*b*d^-1*e^2*(a+b*arctanh(c+d*x))^2*log(2*(1+-1c+-1*d*x)^-1)+-1*b^2*d^-1*e^2*(a+b*arctanh(c+d*x))*Polylogarithms.polylog(2,1+-2*(1+-1c+-1*d*x)^-1))
@test integrate((a+b*arctanh(c+d*x))^3*(c*e+d*e*x), x) == :(-1//2*e*d^-1*(a+b*arctanh(c+d*x))^3+(1/2)*e*d^-1*(a+b*arctanh(c+d*x))^3*(c+d*x)^2+-3//2*e*b^3*d^-1*Polylogarithms.polylog(2,-1*(1+-1c+-1*d*x)^-1*(1+c+d*x))+3//2*b*e*d^-1*(a+b*arctanh(c+d*x))^2+-3*e*b^2*d^-1*(a+b*arctanh(c+d*x))*log(2*(1+-1c+-1*d*x)^-1)+3//2*b*e*d^-1*(a+b*arctanh(c+d*x))^2*(c+d*x))
@test integrate((a+b*arctanh(c+d*x))^3*(c*e+d*e*x)^-1, x) == :(2*d^-1*e^-1*(a+b*arctanh(c+d*x))^3*arctanh(1+-2*(1+-1c+-1*d*x)^-1)+-3//4*b^3*d^-1*e^-1*Polylogarithms.polylog(4,1+-2*(1+-1c+-1*d*x)^-1)+3//4*b^3*d^-1*e^-1*Polylogarithms.polylog(4,-1+2*(1+-1c+-1*d*x)^-1)+-3//2*b*d^-1*e^-1*(a+b*arctanh(c+d*x))^2*Polylogarithms.polylog(2,1+-2*(1+-1c+-1*d*x)^-1)+-3//2*b^2*d^-1*e^-1*(a+b*arctanh(c+d*x))*Polylogarithms.polylog(3,-1+2*(1+-1c+-1*d*x)^-1)+3//2*b*d^-1*e^-1*(a+b*arctanh(c+d*x))^2*Polylogarithms.polylog(2,-1+2*(1+-1c+-1*d*x)^-1)+3//2*b^2*d^-1*e^-1*(a+b*arctanh(c+d*x))*Polylogarithms.polylog(3,1+-2*(1+-1c+-1*d*x)^-1))
@test integrate((a+b*arctanh(c+d*x))^3*(c*e+d*e*x)^-2, x) == :(d^-1*e^-2*(a+b*arctanh(c+d*x))^3+-1*d^-1*e^-2*(a+b*arctanh(c+d*x))^3*(c+d*x)^-1+-3//2*b^3*d^-1*e^-2*Polylogarithms.polylog(3,-1+2*(1+c+d*x)^-1)+-3*b^2*d^-1*e^-2*(a+b*arctanh(c+d*x))*Polylogarithms.polylog(2,-1+2*(1+c+d*x)^-1)+3*b*d^-1*e^-2*(a+b*arctanh(c+d*x))^2*log(2+-2*(1+c+d*x)^-1))
@test integrate((a+b*arctanh(c+d*x))^3*(c*e+d*e*x)^-3, x) == :((1/2)*d^-1*e^-3*(a+b*arctanh(c+d*x))^3+-3//2*b^3*d^-1*e^-3*Polylogarithms.polylog(2,-1+2*(1+c+d*x)^-1)+-1//2*d^-1*e^-3*(a+b*arctanh(c+d*x))^3*(c+d*x)^-2+3//2*b*d^-1*e^-3*(a+b*arctanh(c+d*x))^2+3*b^2*d^-1*e^-3*(a+b*arctanh(c+d*x))*log(2+-2*(1+c+d*x)^-1)+-3//2*b*d^-1*e^-3*(a+b*arctanh(c+d*x))^2*(c+d*x)^-1)
@test integrate((a+b*arctanh(c+d*x))^3*(c*e+d*e*x)^-4, x) == :(1//3*d^-1*e^-4*(a+b*arctanh(c+d*x))^3+b^3*d^-1*e^-4*log(c+d*x)+(1/2)*b*d^-1*e^-4*(a+b*arctanh(c+d*x))^2+-1//2*b^3*d^-1*e^-4*Polylogarithms.polylog(3,-1+2*(1+c+d*x)^-1)+-1//2*b^3*d^-1*e^-4*log(1+-1*(c+d*x)^2)+-1//3*d^-1*e^-4*(a+b*arctanh(c+d*x))^3*(c+d*x)^-3+b*d^-1*e^-4*(a+b*arctanh(c+d*x))^2*log(2+-2*(1+c+d*x)^-1)+-1*b^2*d^-1*e^-4*(c+d*x)^-1*(a+b*arctanh(c+d*x))+-1*b^2*d^-1*e^-4*(a+b*arctanh(c+d*x))*Polylogarithms.polylog(2,-1+2*(1+c+d*x)^-1)+-1//2*b*d^-1*e^-4*(a+b*arctanh(c+d*x))^2*(c+d*x)^-2)
@test integrate((2+2x)^-1*arctanh(1+x), x) == :(-1//4*Polylogarithms.polylog(2,-1+-1x)+1//4*Polylogarithms.polylog(2,1+x))
@test integrate((d*x+a*d*b^-1)^-1*arctanh(a+b*x), x) == :((1/2)*d^-1*Polylogarithms.polylog(2,a+b*x)+-1//2*d^-1*Polylogarithms.polylog(2,-1a+-1*b*x))
@test integrate((e+f*x)^3*(a+b*arctanh(c+d*x)), x) == :(1//4*f^-1*(e+f*x)^4*(a+b*arctanh(c+d*x))+1//12*b*d^-4*f^3*(c+d*x)^3+(1/2)*b*d^-4*f^2*(c+d*x)^2*(d*e+-1*c*f)+-1//8*b*d^-4*f^-1*(-1*f+d*e+-1*c*f)^4*log(1+c+d*x)+1//4*b*f*x*d^-3*(f^2*(1+6*c^2)+6*d^2*e^2+-12*c*d*e*f)+1//8*b*d^-4*f^-1*(f+d*e+-1*c*f)^4*log(1+-1c+-1*d*x))
@test integrate((e+f*x)^2*(a+b*arctanh(c+d*x)), x) == :(1//3*f^-1*(e+f*x)^3*(a+b*arctanh(c+d*x))+1//6*b*d^-3*f^2*(c+d*x)^2+b*f*x*d^-2*(d*e+-1*c*f)+-1//6*b*d^-3*f^-1*(d*e+-1*f*(1+c))^3*log(1+c+d*x)+1//6*b*d^-3*f^-1*(f+d*e+-1*c*f)^3*log(1+-1c+-1*d*x))
@test integrate((a+b*arctanh(c+d*x))*(e+f*x), x) == :((1/2)*f^-1*(e+f*x)^2*(a+b*arctanh(c+d*x))+(1/2)*b*f*x*d^-1+-1//4*b*d^-2*f^-1*(d*e+-1*f*(1+c))^2*log(1+c+d*x)+1//4*b*d^-2*f^-1*(f+d*e+-1*c*f)^2*log(1+-1c+-1*d*x))
@test integrate(a+b*arctanh(c+d*x), x) == :(a*x+(1/2)*b*d^-1*log(1+-1*(c+d*x)^2)+b*d^-1*(c+d*x)*arctanh(c+d*x))
@test integrate((e+f*x)^-1*(a+b*arctanh(c+d*x)), x) == :(f^-1*(a+b*arctanh(c+d*x))*log(2*d*(1+c+d*x)^-1*(f+d*e+-1*c*f)^-1*(e+f*x))+(1/2)*b*f^-1*Polylogarithms.polylog(2,1+-2*(1+c+d*x)^-1)+-1*f^-1*(a+b*arctanh(c+d*x))*log(2*(1+c+d*x)^-1)+-1//2*b*f^-1*Polylogarithms.polylog(2,1+-2*d*(1+c+d*x)^-1*(f+d*e+-1*c*f)^-1*(e+f*x)))
@test integrate((a+b*arctanh(c+d*x))^2*(e+f*x)^3, x) == :(1//4*f^-1*(a+b*arctanh(c+d*x))^2*(e+f*x)^4+d^-4*(a+b*arctanh(c+d*x))^2*(d*e+-1*c*f)*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)+-1//4*d^-4*f^-1*(a+b*arctanh(c+d*x))^2*(d^4*e^4+f^4*(1+c^4+6*c^2)+d^2*e^2*f^2*(6+6*c^2)+-4*c*f*d^3*e^3+-4*c*d*e*f^3*(3+c^2))+1//12*b^2*d^-4*f^3*(c+d*x)^2+1//12*b^2*d^-4*f^3*log(1+-1*(c+d*x)^2)+x*b^2*d^-3*f^2*(d*e+-1*c*f)+-1*b^2*d^-4*f^2*(d*e+-1*c*f)*arctanh(c+d*x)+-1*b^2*d^-4*(d*e+-1*c*f)*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)*Polylogarithms.polylog(2,-1*(1+-1c+-1*d*x)^-1*(1+c+d*x))+1//4*f*b^2*d^-4*(f^2*(1+6*c^2)+6*d^2*e^2+-12*c*d*e*f)*log(1+-1*(c+d*x)^2)+1//6*b*d^-4*f^3*(c+d*x)^3*(a+b*arctanh(c+d*x))+b*d^-4*f^2*(c+d*x)^2*(a+b*arctanh(c+d*x))*(d*e+-1*c*f)+(1/2)*a*b*f*x*d^-3*(f^2*(1+6*c^2)+6*d^2*e^2+-12*c*d*e*f)+(1/2)*f*b^2*d^-4*(c+d*x)*(f^2*(1+6*c^2)+6*d^2*e^2+-12*c*d*e*f)*arctanh(c+d*x)+-2*b*d^-4*(a+b*arctanh(c+d*x))*(d*e+-1*c*f)*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)*log(2*(1+-1c+-1*d*x)^-1))
@test integrate((a+b*arctanh(c+d*x))^2*(e+f*x)^2, x) == :(1//3*d^-3*(a+b*arctanh(c+d*x))^2*(f^2*(1+3*c^2)+3*d^2*e^2+-6*c*d*e*f)+1//3*f^-1*(a+b*arctanh(c+d*x))^2*(e+f*x)^3+-1//3*b^2*d^-3*f^2*arctanh(c+d*x)+-1//3*b^2*d^-3*(f^2*(1+3*c^2)+3*d^2*e^2+-6*c*d*e*f)*Polylogarithms.polylog(2,-1*(1+-1c+-1*d*x)^-1*(1+c+d*x))+1//3*x*b^2*d^-2*f^2+f*b^2*d^-3*(d*e+-1*c*f)*log(1+-1*(c+d*x)^2)+-2//3*b*d^-3*(a+b*arctanh(c+d*x))*(f^2*(1+3*c^2)+3*d^2*e^2+-6*c*d*e*f)*log(2*(1+-1c+-1*d*x)^-1)+-1//3*d^-3*f^-1*(a+b*arctanh(c+d*x))^2*(d*e+-1*c*f)*(d^2*e^2+f^2*(3+c^2)+-2*c*d*e*f)+1//3*b*d^-3*f^2*(c+d*x)^2*(a+b*arctanh(c+d*x))+2*a*b*f*x*d^-2*(d*e+-1*c*f)+2*f*b^2*d^-3*(c+d*x)*(d*e+-1*c*f)*arctanh(c+d*x))
@test integrate((a+b*arctanh(c+d*x))^2*(e+f*x), x) == :(d^-2*(a+b*arctanh(c+d*x))^2*(d*e+-1*c*f)+(1/2)*f^-1*(a+b*arctanh(c+d*x))^2*(e+f*x)^2+(1/2)*f*b^2*d^-2*log(1+-1*(c+d*x)^2)+-1*b^2*d^-2*(d*e+-1*c*f)*Polylogarithms.polylog(2,-1*(1+-1c+-1*d*x)^-1*(1+c+d*x))+-1//2*d^-2*f^-1*(a+b*arctanh(c+d*x))^2*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)+a*b*f*x*d^-1+f*b^2*d^-2*(c+d*x)*arctanh(c+d*x)+-2*b*d^-2*(a+b*arctanh(c+d*x))*(d*e+-1*c*f)*log(2*(1+-1c+-1*d*x)^-1))
@test integrate((a+b*arctanh(c+d*x))^2, x) == :(d^-1*(a+b*arctanh(c+d*x))^2+d^-1*(a+b*arctanh(c+d*x))^2*(c+d*x)+-1*b^2*d^-1*Polylogarithms.polylog(2,-1*(1+-1c+-1*d*x)^-1*(1+c+d*x))+-2*b*d^-1*(a+b*arctanh(c+d*x))*log(2*(1+-1c+-1*d*x)^-1))
@test integrate((a+b*arctanh(c+d*x))^2*(e+f*x)^-1, x) == :(f^-1*(a+b*arctanh(c+d*x))^2*log(2*d*(1+c+d*x)^-1*(f+d*e+-1*c*f)^-1*(e+f*x))+(1/2)*b^2*f^-1*Polylogarithms.polylog(3,1+-2*(1+c+d*x)^-1)+-1*f^-1*(a+b*arctanh(c+d*x))^2*log(2*(1+c+d*x)^-1)+-1//2*b^2*f^-1*Polylogarithms.polylog(3,1+-2*d*(1+c+d*x)^-1*(f+d*e+-1*c*f)^-1*(e+f*x))+b*f^-1*(a+b*arctanh(c+d*x))*Polylogarithms.polylog(2,1+-2*(1+c+d*x)^-1)+-1*b*f^-1*(a+b*arctanh(c+d*x))*Polylogarithms.polylog(2,1+-2*d*(1+c+d*x)^-1*(f+d*e+-1*c*f)^-1*(e+f*x)))
@test integrate((a+b*arctanh(c+d*x))^3*(e+f*x)^2, x) == :(1//3*d^-3*(a+b*arctanh(c+d*x))^3*(f^2*(1+3*c^2)+3*d^2*e^2+-6*c*d*e*f)+1//3*f^-1*(a+b*arctanh(c+d*x))^3*(e+f*x)^3+(1/2)*b^3*d^-3*f^2*log(1+-1*(c+d*x)^2)+(1/2)*b^3*d^-3*(f^2*(1+3*c^2)+3*d^2*e^2+-6*c*d*e*f)*Polylogarithms.polylog(3,1+-2*(1+-1c+-1*d*x)^-1)+-1//2*b*d^-3*f^2*(a+b*arctanh(c+d*x))^2+a*x*b^2*d^-2*f^2+b^3*d^-3*f^2*(c+d*x)*arctanh(c+d*x)+(1/2)*b*d^-3*f^2*(a+b*arctanh(c+d*x))^2*(c+d*x)^2+-1*b*d^-3*(a+b*arctanh(c+d*x))^2*(f^2*(1+3*c^2)+3*d^2*e^2+-6*c*d*e*f)*log(2*(1+-1c+-1*d*x)^-1)+-1*b^2*d^-3*(a+b*arctanh(c+d*x))*(f^2*(1+3*c^2)+3*d^2*e^2+-6*c*d*e*f)*Polylogarithms.polylog(2,1+-2*(1+-1c+-1*d*x)^-1)+-3*f*b^3*d^-3*(d*e+-1*c*f)*Polylogarithms.polylog(2,-1*(1+-1c+-1*d*x)^-1*(1+c+d*x))+3*b*f*d^-3*(a+b*arctanh(c+d*x))^2*(d*e+-1*c*f)+-1//3*d^-3*f^-1*(a+b*arctanh(c+d*x))^3*(d*e+-1*c*f)*(d^2*e^2+f^2*(3+c^2)+-2*c*d*e*f)+-6*f*b^2*d^-3*(a+b*arctanh(c+d*x))*(d*e+-1*c*f)*log(2*(1+-1c+-1*d*x)^-1)+3*b*f*d^-3*(a+b*arctanh(c+d*x))^2*(c+d*x)*(d*e+-1*c*f))
@test integrate((a+b*arctanh(c+d*x))^3*(e+f*x), x) == :(d^-2*(a+b*arctanh(c+d*x))^3*(d*e+-1*c*f)+(1/2)*f^-1*(a+b*arctanh(c+d*x))^3*(e+f*x)^2+-3//2*f*b^3*d^-2*Polylogarithms.polylog(2,-1*(1+-1c+-1*d*x)^-1*(1+c+d*x))+-1//2*d^-2*f^-1*(a+b*arctanh(c+d*x))^3*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)+3//2*b*f*d^-2*(a+b*arctanh(c+d*x))^2+3//2*b^3*d^-2*(d*e+-1*c*f)*Polylogarithms.polylog(3,1+-2*(1+-1c+-1*d*x)^-1)+-3*b*d^-2*(a+b*arctanh(c+d*x))^2*(d*e+-1*c*f)*log(2*(1+-1c+-1*d*x)^-1)+-3*f*b^2*d^-2*(a+b*arctanh(c+d*x))*log(2*(1+-1c+-1*d*x)^-1)+-3*b^2*d^-2*(a+b*arctanh(c+d*x))*(d*e+-1*c*f)*Polylogarithms.polylog(2,1+-2*(1+-1c+-1*d*x)^-1)+3//2*b*f*d^-2*(a+b*arctanh(c+d*x))^2*(c+d*x))
@test integrate((a+b*arctanh(c+d*x))^3, x) == :(d^-1*(a+b*arctanh(c+d*x))^3+d^-1*(a+b*arctanh(c+d*x))^3*(c+d*x)+3//2*b^3*d^-1*Polylogarithms.polylog(3,1+-2*(1+-1c+-1*d*x)^-1)+-3*b*d^-1*(a+b*arctanh(c+d*x))^2*log(2*(1+-1c+-1*d*x)^-1)+-3*b^2*d^-1*(a+b*arctanh(c+d*x))*Polylogarithms.polylog(2,1+-2*(1+-1c+-1*d*x)^-1))
@test integrate((a+b*arctanh(c+d*x))^3*(e+f*x)^-1, x) == :(f^-1*(a+b*arctanh(c+d*x))^3*log(2*d*(1+c+d*x)^-1*(f+d*e+-1*c*f)^-1*(e+f*x))+-1*f^-1*(a+b*arctanh(c+d*x))^3*log(2*(1+c+d*x)^-1)+-3//4*b^3*f^-1*Polylogarithms.polylog(4,1+-2*d*(1+c+d*x)^-1*(f+d*e+-1*c*f)^-1*(e+f*x))+3//4*b^3*f^-1*Polylogarithms.polylog(4,1+-2*(1+c+d*x)^-1)+-3//2*b*f^-1*(a+b*arctanh(c+d*x))^2*Polylogarithms.polylog(2,1+-2*d*(1+c+d*x)^-1*(f+d*e+-1*c*f)^-1*(e+f*x))+-3//2*b^2*f^-1*(a+b*arctanh(c+d*x))*Polylogarithms.polylog(3,1+-2*d*(1+c+d*x)^-1*(f+d*e+-1*c*f)^-1*(e+f*x))+3//2*b*f^-1*(a+b*arctanh(c+d*x))^2*Polylogarithms.polylog(2,1+-2*(1+c+d*x)^-1)+3//2*b^2*f^-1*(a+b*arctanh(c+d*x))*Polylogarithms.polylog(3,1+-2*(1+c+d*x)^-1))
@test integrate((c+d*x^3)^-1*arctanh(a+b*x), x) == :(-1//6*c^-2//3*d^-1//3*Polylogarithms.polylog(2,d^1//3*(b*c^1//3+d^1//3*(1+-1a))^-1*(1+-1a+-1*b*x))+1//6*c^-2//3*d^-1//3*Polylogarithms.polylog(2,-1*d^1//3*(b*c^1//3+-1*d^1//3*(1+a))^-1*(1+a+b*x))+-1//6*(-1)^1//3*c^-2//3*d^-1//3*Polylogarithms.polylog(2,-1*(-1)^2//3*d^1//3*(b*c^1//3+-1*(-1)^2//3*d^1//3*(1+a))^-1*(1+a+b*x))+-1//6*(-1)^2//3*c^-2//3*d^-1//3*Polylogarithms.polylog(2,-1*(-1)^1//3*d^1//3*(b*c^1//3+-1*(-1)^1//3*d^1//3*(1+-1a))^-1*(1+-1a+-1*b*x))+-1//6*c^-2//3*d^-1//3*log(b*(b*c^1//3+d^1//3*(1+-1a))^-1*(c^1//3+x*d^1//3))*log(1+-1a+-1*b*x)+1//6*(-1)^1//3*c^-2//3*d^-1//3*Polylogarithms.polylog(2,(-1)^2//3*d^1//3*(b*c^1//3+(-1)^2//3*d^1//3*(1+-1a))^-1*(1+-1a+-1*b*x))+1//6*(-1)^2//3*c^-2//3*d^-1//3*Polylogarithms.polylog(2,(-1)^1//3*d^1//3*(b*c^1//3+(-1)^1//3*d^1//3*(1+a))^-1*(1+a+b*x))+1//6*c^-2//3*d^-1//3*log(b*(b*c^1//3+-1*d^1//3*(1+a))^-1*(c^1//3+x*d^1//3))*log(1+a+b*x)+-1//6*(-1)^1//3*c^-2//3*d^-1//3*log(b*(b*c^1//3+-1*(-1)^2//3*d^1//3*(1+a))^-1*(c^1//3+x*(-1)^2//3*d^1//3))*log(1+a+b*x)+-1//6*(-1)^2//3*c^-2//3*d^-1//3*log(b*(b*c^1//3+-1*(-1)^1//3*d^1//3*(1+-1a))^-1*(c^1//3+-1*x*(-1)^1//3*d^1//3))*log(1+-1a+-1*b*x)+1//6*(-1)^1//3*c^-2//3*d^-1//3*log(b*(b*c^1//3+(-1)^2//3*d^1//3*(1+-1a))^-1*(c^1//3+x*(-1)^2//3*d^1//3))*log(1+-1a+-1*b*x)+1//6*(-1)^2//3*c^-2//3*d^-1//3*log(b*(b*c^1//3+(-1)^1//3*d^1//3*(1+a))^-1*(c^1//3+-1*x*(-1)^1//3*d^1//3))*log(1+a+b*x))
@test integrate((c+d*x^2)^-1*arctanh(a+b*x), x) == :(-1//4*d^-1//2*(-1c)^-1//2*Polylogarithms.polylog(2,-1*d^(1/2)*(b*(-1c)^(1/2)+-1*d^(1/2)*(1+a))^-1*(1+a+b*x))+-1//4*d^-1//2*(-1c)^-1//2*Polylogarithms.polylog(2,-1*d^(1/2)*(b*(-1c)^(1/2)+-1*d^(1/2)*(1+-1a))^-1*(1+-1a+-1*b*x))+1//4*d^-1//2*(-1c)^-1//2*Polylogarithms.polylog(2,d^(1/2)*(b*(-1c)^(1/2)+d^(1/2)*(1+a))^-1*(1+a+b*x))+1//4*d^-1//2*(-1c)^-1//2*Polylogarithms.polylog(2,d^(1/2)*(b*(-1c)^(1/2)+d^(1/2)*(1+-1a))^-1*(1+-1a+-1*b*x))+-1//4*d^-1//2*(-1c)^-1//2*log(b*(b*(-1c)^(1/2)+-1*d^(1/2)*(1+a))^-1*((-1c)^(1/2)+x*d^(1/2)))*log(1+a+b*x)+-1//4*d^-1//2*(-1c)^-1//2*log(b*(b*(-1c)^(1/2)+-1*d^(1/2)*(1+-1a))^-1*((-1c)^(1/2)+-1*x*d^(1/2)))*log(1+-1a+-1*b*x)+1//4*d^-1//2*(-1c)^-1//2*log(b*(b*(-1c)^(1/2)+d^(1/2)*(1+a))^-1*((-1c)^(1/2)+-1*x*d^(1/2)))*log(1+a+b*x)+1//4*d^-1//2*(-1c)^-1//2*log(b*(b*(-1c)^(1/2)+d^(1/2)*(1+-1a))^-1*((-1c)^(1/2)+x*d^(1/2)))*log(1+-1a+-1*b*x))
@test integrate((c+d*x)^-1*arctanh(a+b*x), x) == :((1/2)*d^-1*Polylogarithms.polylog(2,1+-2*(1+a+b*x)^-1)+-1//2*d^-1*Polylogarithms.polylog(2,1+-2*b*(1+a+b*x)^-1*(d+b*c+-1*a*d)^-1*(c+d*x))+d^-1*arctanh(a+b*x)*log(2*b*(1+a+b*x)^-1*(d+b*c+-1*a*d)^-1*(c+d*x))+-1*d^-1*arctanh(a+b*x)*log(2*(1+a+b*x)^-1))
@test integrate((c+d*x^-1)^-1*arctanh(a+b*x), x) == :((1/2)*d*c^-2*Polylogarithms.polylog(2,c*(c+b*d+-1*a*c)^-1*(1+-1a+-1*b*x))+-1//2*d*c^-2*Polylogarithms.polylog(2,c*(c+a*c+-1*b*d)^-1*(1+a+b*x))+(1/2)*d*c^-2*log(b*(c+b*d+-1*a*c)^-1*(d+c*x))*log(1+-1a+-1*b*x)+(1/2)*b^-1*c^-1*(1+a+b*x)*log(1+a+b*x)+(1/2)*b^-1*c^-1*(1+-1a+-1*b*x)*log(1+-1a+-1*b*x)+-1//2*d*c^-2*log(-1*b*(c+a*c+-1*b*d)^-1*(d+c*x))*log(1+a+b*x))
@test integrate((c+d*x^-2)^-1*arctanh(a+b*x), x) == :(-1//4*d^(1/2)*(-1c)^-3//2*Polylogarithms.polylog(2,(-1c)^(1/2)*(b*d^(1/2)+(-1c)^(1/2)*(1+a))^-1*(1+a+b*x))+-1//4*d^(1/2)*(-1c)^-3//2*Polylogarithms.polylog(2,(-1c)^(1/2)*(b*d^(1/2)+(-1c)^(1/2)*(1+-1a))^-1*(1+-1a+-1*b*x))+1//4*d^(1/2)*(-1c)^-3//2*Polylogarithms.polylog(2,(-1c)^(1/2)*((-1c)^(1/2)*(1+a)+-1*b*d^(1/2))^-1*(1+a+b*x))+1//4*d^(1/2)*(-1c)^-3//2*Polylogarithms.polylog(2,(-1c)^(1/2)*((-1c)^(1/2)+-1*a*(-1c)^(1/2)+-1*b*d^(1/2))^-1*(1+-1a+-1*b*x))+(1/2)*b^-1*c^-1*(1+a+b*x)*log(1+a+b*x)+(1/2)*b^-1*c^-1*(1+-1a+-1*b*x)*log(1+-1a+-1*b*x)+-1//4*d^(1/2)*(-1c)^-3//2*log(b*(b*d^(1/2)+(-1c)^(1/2)*(1+a))^-1*(d^(1/2)+-1*x*(-1c)^(1/2)))*log(1+a+b*x)+-1//4*d^(1/2)*(-1c)^-3//2*log(b*(b*d^(1/2)+(-1c)^(1/2)*(1+-1a))^-1*(d^(1/2)+x*(-1c)^(1/2)))*log(1+-1a+-1*b*x)+1//4*d^(1/2)*(-1c)^-3//2*log(-1*b*((-1c)^(1/2)*(1+a)+-1*b*d^(1/2))^-1*(d^(1/2)+x*(-1c)^(1/2)))*log(1+a+b*x)+1//4*d^(1/2)*(-1c)^-3//2*log(-1*b*((-1c)^(1/2)*(1+-1a)+-1*b*d^(1/2))^-1*(d^(1/2)+-1*x*(-1c)^(1/2)))*log(1+-1a+-1*b*x))
@test integrate((c+d*x^-3)^-1*arctanh(a+b*x), x) == :(-1//6*c^-4//3*d^1//3*Polylogarithms.polylog(2,c^1//3*(c^1//3*(1+a)+-1*b*d^1//3)^-1*(1+a+b*x))+1//6*c^-4//3*d^1//3*Polylogarithms.polylog(2,c^1//3*(b*d^1//3+c^1//3*(1+-1a))^-1*(1+-1a+-1*b*x))+(1/2)*b^-1*c^-1*(1+a+b*x)*log(1+a+b*x)+(1/2)*b^-1*c^-1*(1+-1a+-1*b*x)*log(1+-1a+-1*b*x)+-1//6*(-1)^1//3*c^-4//3*d^1//3*Polylogarithms.polylog(2,(-1)^2//3*c^1//3*(b*d^1//3+(-1)^2//3*c^1//3*(1+-1a))^-1*(1+-1a+-1*b*x))+-1//6*(-1)^2//3*c^-4//3*d^1//3*Polylogarithms.polylog(2,(-1)^1//3*c^1//3*(b*d^1//3+(-1)^1//3*c^1//3*(1+a))^-1*(1+a+b*x))+-1//6*c^-4//3*d^1//3*log(-1*b*(c^1//3*(1+a)+-1*b*d^1//3)^-1*(d^1//3+x*c^1//3))*log(1+a+b*x)+1//6*(-1)^1//3*c^-4//3*d^1//3*Polylogarithms.polylog(2,(-1)^2//3*c^1//3*(-1*b*d^1//3+(-1)^2//3*c^1//3*(1+a))^-1*(1+a+b*x))+1//6*(-1)^2//3*c^-4//3*d^1//3*Polylogarithms.polylog(2,(-1)^1//3*c^1//3*(-1*b*d^1//3+(-1)^1//3*c^1//3*(1+-1a))^-1*(1+-1a+-1*b*x))+1//6*c^-4//3*d^1//3*log(b*(b*d^1//3+c^1//3*(1+-1a))^-1*(d^1//3+x*c^1//3))*log(1+-1a+-1*b*x)+-1//6*(-1)^1//3*c^-4//3*d^1//3*log(b*(b*d^1//3+(-1)^2//3*c^1//3*(1+-1a))^-1*(d^1//3+x*(-1)^2//3*c^1//3))*log(1+-1a+-1*b*x)+-1//6*(-1)^2//3*c^-4//3*d^1//3*log(b*(b*d^1//3+(-1)^1//3*c^1//3*(1+a))^-1*(d^1//3+-1*x*(-1)^1//3*c^1//3))*log(1+a+b*x)+1//6*(-1)^1//3*c^-4//3*d^1//3*log(-1*b*(-1*b*d^1//3+(-1)^2//3*c^1//3*(1+a))^-1*(d^1//3+x*(-1)^2//3*c^1//3))*log(1+a+b*x)+1//6*(-1)^2//3*c^-4//3*d^1//3*log(-1*b*(-1*b*d^1//3+(-1)^1//3*c^1//3*(1+-1a))^-1*(d^1//3+-1*x*(-1)^1//3*c^1//3))*log(1+-1a+-1*b*x))
@test integrate((c+d*x^(1/2))^-1*arctanh(a+b*x), x) == :(c*d^-2*Polylogarithms.polylog(2,b^(1/2)*(c*b^(1/2)+d*(-1+-1a)^(1/2))^-1*(c+d*x^(1/2)))+c*d^-2*Polylogarithms.polylog(2,b^(1/2)*(c*b^(1/2)+-1*d*(-1+-1a)^(1/2))^-1*(c+d*x^(1/2)))+d^-1*x^(1/2)*log(1+a+b*x)+-1*c*d^-2*Polylogarithms.polylog(2,b^(1/2)*(c*b^(1/2)+d*(1+-1a)^(1/2))^-1*(c+d*x^(1/2)))+-1*c*d^-2*Polylogarithms.polylog(2,b^(1/2)*(c*b^(1/2)+-1*d*(1+-1a)^(1/2))^-1*(c+d*x^(1/2)))+-1*d^-1*x^(1/2)*log(1+-1a+-1*b*x)+c*d^-2*log(d*(c*b^(1/2)+d*(-1+-1a)^(1/2))^-1*((-1+-1a)^(1/2)+-1*b^(1/2)*x^(1/2)))*log(c+d*x^(1/2))+c*d^-2*log(-1*d*(c*b^(1/2)+-1*d*(-1+-1a)^(1/2))^-1*((-1+-1a)^(1/2)+b^(1/2)*x^(1/2)))*log(c+d*x^(1/2))+c*d^-2*log(c+d*x^(1/2))*log(1+-1a+-1*b*x)+-1*c*d^-2*log(d*(c*b^(1/2)+d*(1+-1a)^(1/2))^-1*((1+-1a)^(1/2)+-1*b^(1/2)*x^(1/2)))*log(c+d*x^(1/2))+-1*c*d^-2*log(-1*d*(c*b^(1/2)+-1*d*(1+-1a)^(1/2))^-1*((1+-1a)^(1/2)+b^(1/2)*x^(1/2)))*log(c+d*x^(1/2))+-1*c*d^-2*log(c+d*x^(1/2))*log(1+a+b*x)+-2*b^-1//2*d^-1*(1+-1a)^(1/2)*arctanh(b^(1/2)*x^(1/2)*(1+-1a)^-1//2)+2*b^-1//2*d^-1*(1+a)^(1/2)*arctan(b^(1/2)*x^(1/2)*(1+a)^-1//2))
@test integrate((c+d*x^-1//2)^-1*arctanh(a+b*x), x) == :(c^-3*d^2*Polylogarithms.polylog(2,b^(1/2)*(c*(1+-1a)^(1/2)+d*b^(1/2))^-1*(d+c*x^(1/2)))+c^-3*d^2*Polylogarithms.polylog(2,-1*b^(1/2)*(c*(1+-1a)^(1/2)+-1*d*b^(1/2))^-1*(d+c*x^(1/2)))+-1*c^-3*d^2*Polylogarithms.polylog(2,b^(1/2)*(c*(-1+-1a)^(1/2)+d*b^(1/2))^-1*(d+c*x^(1/2)))+-1*c^-3*d^2*Polylogarithms.polylog(2,-1*b^(1/2)*(c*(-1+-1a)^(1/2)+-1*d*b^(1/2))^-1*(d+c*x^(1/2)))+d*c^-2*x^(1/2)*log(1+-1a+-1*b*x)+c^-3*d^2*log(c*(c*(1+-1a)^(1/2)+d*b^(1/2))^-1*((1+-1a)^(1/2)+-1*b^(1/2)*x^(1/2)))*log(d+c*x^(1/2))+c^-3*d^2*log(c*(c*(1+-1a)^(1/2)+-1*d*b^(1/2))^-1*((1+-1a)^(1/2)+b^(1/2)*x^(1/2)))*log(d+c*x^(1/2))+c^-3*d^2*log(d+c*x^(1/2))*log(1+a+b*x)+(1/2)*b^-1*c^-1*(1+a+b*x)*log(1+a+b*x)+(1/2)*b^-1*c^-1*(1+-1a+-1*b*x)*log(1+-1a+-1*b*x)+-1*d*c^-2*x^(1/2)*log(1+a+b*x)+-1*c^-3*d^2*log(c*(c*(-1+-1a)^(1/2)+d*b^(1/2))^-1*((-1+-1a)^(1/2)+-1*b^(1/2)*x^(1/2)))*log(d+c*x^(1/2))+-1*c^-3*d^2*log(c*(c*(-1+-1a)^(1/2)+-1*d*b^(1/2))^-1*((-1+-1a)^(1/2)+b^(1/2)*x^(1/2)))*log(d+c*x^(1/2))+-1*c^-3*d^2*log(d+c*x^(1/2))*log(1+-1a+-1*b*x)+-2*d*b^-1//2*c^-2*(1+a)^(1/2)*arctan(b^(1/2)*x^(1/2)*(1+a)^-1//2)+2*d*b^-1//2*c^-2*(1+-1a)^(1/2)*arctanh(b^(1/2)*x^(1/2)*(1+-1a)^-1//2))
@test integrate((a+b*x+c*x^2)^-1*arctanh(d+e*x), x) == :((1/2)*(b^2+-4*a*c)^-1//2*Polylogarithms.polylog(2,1+(e*(b+(b^2+-4*a*c)^(1/2))+2*c*(1+-1d))^-1*(1+d+e*x)^-1*(-4*c*(d+e*x)+-2*e*(b+(b^2+-4*a*c)^(1/2))+4*c*d))+-1//2*(b^2+-4*a*c)^-1//2*Polylogarithms.polylog(2,1+(1+d+e*x)^-1*(2c+b*e+-1*e*(b^2+-4*a*c)^(1/2)+-2*c*d)^-1*(-4*c*(d+e*x)+-2*e*(b+-1*(b^2+-4*a*c)^(1/2))+4*c*d))+(b^2+-4*a*c)^-1//2*arctanh(d+e*x)*log(2*e*(e*(b+-1*(b^2+-4*a*c)^(1/2))+2*c*(1+-1d))^-1*(1+d+e*x)^-1*(b+-1*(b^2+-4*a*c)^(1/2)+2*c*x))+-1*(b^2+-4*a*c)^-1//2*arctanh(d+e*x)*log(2*e*(e*(b+(b^2+-4*a*c)^(1/2))+2*c*(1+-1d))^-1*(1+d+e*x)^-1*(b+(b^2+-4*a*c)^(1/2)+2*c*x)))
@test integrate((1+-1*(c+d*x)^2)^-1*(a+b*arctanh(c+d*x))*(c*e+d*e*x), x) == :(e*d^-1*(a+b*arctanh(c+d*x))*log(2*(1+-1c+-1*d*x)^-1)+(1/2)*b*e*d^-1*Polylogarithms.polylog(2,-1*(1+-1c+-1*d*x)^-1*(1+c+d*x))+-1//2*e*b^-1*d^-1*(a+b*arctanh(c+d*x))^2)

max_stars_repo_path: test/7 Inverse hyperbolic functions/7.3 Inverse hyperbolic tangent/7.3.5 u (a+b arctanh(c+d x))^p.jl


text:
# Digit factorials

function solve()
    factorials = [factorial(i) for i in 0:9]
    return sum(n for n in 3:2_540_160
               if sum(factorials[parse(Int, i) + 1] for i in string(n)) == n)
end

print(Int(solve()))

max_stars_repo_name: wephy/project-euler

max_issues_repo_name: wephy/project-euler

max_forks_repo_name: wephy/project-euler


text:
using LTC
using Plots
gr()
using BenchmarkTools
using DiffEqSensitivity
using OrdinaryDiffEq
using DiffEqFlux
using GalacticOptim
# using BlackBoxOptim
using ModelingToolkit
using IterTools: ncycle

include("half_cheetah_data_loader.jl")

function train_cheetah(n, solver=VCABM(), sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true));)

  cbg = function (p,l,pred,y;doplot=true)
    display(l)
    if doplot
      fig = plot([ŷ[1,1] for ŷ in pred], label="ŷ1")
      plot!(fig, [ŷ[2,1] for ŷ in pred], label="ŷ2")
      plot!(fig, [yi[1,1] for yi in y], label="y1")
      plot!(fig, [yi[2,1] for yi in y], label="y2")
      display(fig)
    end
    return false
  end

  batchsize=10
  seq_len=20
  train_dl, test_dl, valid_dl = get_dl(batchsize=batchsize, seq_len=seq_len)

  wiring = LTC.NCPWiring(17,17;
    n_sensory=17, n_inter=5, n_command=4, n_motor=17,
    sensory_in=-1, rec_sensory=4, sensory_inter=2, sensory_command=1, sensory_motor=0,
    inter_in=2, rec_inter=2, inter_command=3, inter_motor=1,                       # inter_in = sensory_out
    command_in=0, rec_command=4, command_motor=4,                   # command_in = inter_out
    motor_in=0, rec_motor=2)

  net = LTC.Net(wiring, name=:net)
  sys = ModelingToolkit.structural_simplify(net)

  model = DiffEqFlux.FastChain(LTC.Mapper(wiring.n_in),
                               LTC.RecurMTK(LTC.MTKCell(wiring.n_in, wiring.n_out, net, sys, solver, sensealg)),
                               LTC.Mapper(wiring.n_out),
                               )

  opt = Flux.Optimiser(ClipValue(1.00), ExpDecay(1, 0.1, 200, 1e-4), ADAM())
  # opt = Optim.LBFGS()
  # opt = BBO()
  # opt = ParticleSwarm(;lower=lb, upper=ub)
  # opt = Fminbox(GradientDescent())
  AD = GalacticOptim.AutoZygote()
  # AD = GalacticOptim.AutoModelingToolkit()
  LTC.optimize(model, LTC.loss_seq, cbg, opt, AD, ncycle(train_dl,n)), model
end

@time res1,model = train_cheetah(100)
# @time res1,model = train_cheetah(500, AutoTsit5(Rosenbrock23()))
# @time res1,model = train_cheetah(5, VCABM(), InterpolatingAdjoint(checkpointing=true, autojacvec=ReverseDiffVJP(true)))
# @time res1,model = train_cheetah(5, Tsit5(), InterpolatingAdjoint(checkpointing=true))
# @time res1,model = train_cheetah(5, AutoTsit5(Rosenbrock23()), InterpolatingAdjoint(checkpointing=true, autojacvec=ReverseDiffVJP(true)))
# @time res1,model = train_cheetah(5, Tsit5(), InterpolatingAdjoint(checkpointing=true, autojacvec=ReverseDiffVJP(true)))
# @time res1,model = train_cheetah(1, VCABM(), BacksolveAdjoint())


text:
mutable struct S a::Int64 end
function f(s::S) s.a = "foo" end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia