text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

using FactCheck
using Latte

net = Net(8)
data,  data_value   = MemoryDataLayer(net, :data, (227, 227, 3))
label, label_value = MemoryDataLayer(net, :label, (1,))
data_value[:]  = rand(Float32, size(data_value)...) * 256
label_value[:] = map(floor, rand(Float32, size(label_value)...) * 10)
fc1         = InnerProductLayer(:fc1, net, data, 10)
softmax     = SoftmaxLossLayer(:softmax, net, fc1, label)

init(net)

params = SolverParameters(
    lr_policy    = LRPolicy.Decay(.01f0, 5.0f-7),
    mom_policy   = MomPolicy.Fixed(0.9),
    max_epoch    = 300,
    regu_coef    = .0005)
sgd = SGD(params)

input = get_buffer(net, :fc1value)
label = get_buffer(net, :labelvalue)

expected_grad = zeros(Float32, size(input))

epsilon = 1e-5
forward(net; solver=sgd)
batch_size = size(input)[end]
pred = zeros(input)
expected_loss = 0
for n in 1:batch_size
    pred[:,n] = exp(input[:,n] .- maximum(input[:,n]))
    pred[:,n] /= sum(pred[:,n])
    label_value = round(Int, label[1, n]) + 1
    expected_loss -= log(max(pred[label_value, n], eps(Float32)))
end
expected_loss /= batch_size

facts("Testing SoftmaxLoss layer") do
    context("Forward") do
        @fact get_buffer(net, :softmaxvalue)[1] --> roughly(expected_loss)
    end
    context("Backward") do
        backward(net)
        expected_grad = copy(pred)
        for n in 1:batch_size
            label_value = round(Int, label[1, n]) + 1
            expected_grad[label_value, n] -= 1
        end
        expected_grad /= batch_size
        @fact get_buffer(net, :fc1∇) --> roughly(expected_grad)
    end
end

FactCheck.exitstatus()


text:
# struct Simulation{DT <: Number, ET <: Equation{DT}}
#     Δt::DT
#     equ::ET
#     x::OffsetArray{DT,3,Array{DT,3}}

#     function Simulation(equ::ET, Δt::DT, nt::Int) where {DT, ET <: Equation{DT}}
#         x = OffsetArray(zeros(DT, size(equ.x₀)..., nt+1), axes(equ.x₀)..., 0:nt)
#         x[:,:,0] .= equ.x₀
#         new{DT,ET}(Δt, equ, x)
#     end
# end

struct Simulation{DT <: Number, FT <: Function}
    Δt::DT
    equ::Equation{DT,FT}
    x::OffsetArray{DT,3,Array{DT,3}}

    function Simulation(equ::Equation{DT,FT}, Δt::DT, nt::Int) where {DT, FT}
        x = OffsetArray(zeros(DT, size(equ.x₀)..., nt+1), axes(equ.x₀)..., 0:nt)
        x[:,:,0] .= equ.x₀
        new{DT,FT}(Δt, equ, x)
    end
end

ndims(sim::Simulation) = ndims(sim.equ)
nsamples(sim::Simulation) = nsamples(sim.equ)
ntimesteps(sim::Simulation) = lastindex(sim.x,3)

eachsample(sim::Simulation) = axes(sim.x,2)
eachtimestep(sim::Simulation) = axes(sim.x,3)[1:end]

function run!(sim::Simulation{DT}) where {DT}
    ẋ = zeros(DT, ndims(sim))
    for n in eachtimestep(sim)
        for i in eachsample(sim)
            sim.equ.f(ẋ, sim.x[:,i,n-1])
            sim.x[:,i,n] .= sim.x[:,i,n-1] .+ sim.Δt .* ẋ
        end
    end
    return sim.x
end

ext: jl

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
using Documenter
using BioMetaheuristics, BioMetaheuristics.TestFunctions

# ! Build the full Documentation with Documenter
makedocs(;
    modules=[BioMetaheuristics, BioMetaheuristics.TestFunctions],
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://edwinb-ai.github.io/BioMetaheuristics.jl/stable/",
        assets=String[]
    ),
    pages=[
        "Home" => "index.md",
        "Theory" => "theory.md",
        "Examples" => ["examples.md", "api.md"],
        "Implementations" => "algorithms.md",
        "Benchmark functions" => "benchmarks.md",
        "Reference" => "reference.md",
        "License" => "license.md",
    ],
    repo="https://github.com/edwinb-ai/BioMetaheuristics.jl/blob/{commit}{path}#L{line}",
    sitename="BioMetaheuristics.jl",
    authors="Edwin Bedolla"
)

deploydocs(
    repo="github.com/edwinb-ai/BioMetaheuristics.jl.git",
    push_preview=true
)


text:
using GLVisualize, Colors, FixedPointNumbers
using Base.Test

@test GLVisualize.Intensity{Float32}(Gray(N0f8(0.8))) == GLVisualize.Intensity{Float32}(Float32(N0f8(0.8)))

nothing


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export arcmap
"""
`arcmap(f::Function, fst::WFST, args...; modify_initials=identity, modify_finals=identity, isym=get_isym(fst), osym=get_osym(fst))`

Creates a new WFST from `fst` by applying the function `f` to all its arcs, see [`Arc`](@ref).
The arguments of `f` can be specified in `args`.

The functions `modify_initials` and `modify_finals` operate in the initial and final dictionaries, which keys are the initial/final states and values are the corresponding state weights.

Input and output tables can also be modified using the keywords `isym` and `osym`.

The following example shows how to use `arcmap` to convert the type of weight of a WFST:
```julia
julia> A = WFST(["a","b","c"]); # by default weight is TropicalWeight{Float32}

julia> add_arc!(A,1=>2,"a"=>"a",1);

julia> add_arc!(A,1=>3,"b"=>"c",3);

julia> initial!(A,1); final!(A,2,4); final!(A,3,2)
WFST #states: 3, #arcs: 2, #isym: 3, #osym: 3
|1/0.0f0|
a:a/3.0f0 → (2)
b:c/3.0f0 → (3)
((2/4.0f0))
((3/2.0f0))

julia> trop2prob(x) = ProbabilityWeight{Float64}(exp(-get(x)))
trop2prob (generic function with 1 method)

julia> function trop2prob(arc::Arc)
           ilab = get_ilabel(arc)
           olab = get_olabel(arc)
           w = trop2prob(get_weight(arc))
           n = get_nextstate(arc)
           return Arc(ilab,olab,w,n)
       end
trop2prob (generic function with 2 methods)

julia> trop2prob(initials::Dict) = Dict(i => trop2prob(w) for (i,w) in initials)
trop2prob (generic function with 3 methods)

julia> arcmap(trop2prob,A; modify_initials=trop2prob, modify_finals=trop2prob)
WFST #states: 3, #arcs: 2, #isym: 3, #osym: 3
|1/1.0|
a:a/0.3678794503211975 → (2)
b:c/0.049787066876888275 → (3)
((2/0.018315639346837997))
((3/0.1353352814912796))

```
"""
function arcmap(f::Function, fst::WFST, args...;
                modify_initials = identity,
                modify_finals = identity,
                isym = get_isym(fst),
                osym = get_osym(fst)
               )
  states = [f.(arcs, args...) for arcs in get_states(fst)]
  finals   = modify_finals(  deepcopy(fst.finals)  )
  initials = modify_initials(deepcopy(fst.initials))
  osym = deepcopy(fst.osym)
  isym = deepcopy(fst.isym)
  pfst = WFST(states,initials,finals,isym,osym)
  return pfst
end


text:
export ls, regex_find, safe_isdir, safe_isfile

# safe_isfile, safe_isdir adapted from https://github.com/JuliaPackaging/BinaryProvider.jl/commit/08a314a225206a68665c6f730d7c3feeda1ba615
# Temporary hack around https://github.com/JuliaLang/julia/issues/26685
function safe_isfile(path::String)
  try
    return isfile(path)
  catch err
    return false
  end
end

function safe_isdir(path::String)
  try
    return isdir(path)
  catch err
    return false
  end
end

"""
    find_regex(path::String, r::Regex)

OS-agnostic equivalent to Linux `find`. First argument is a path string, second is a Regex. File strings are postprocessed using Julia's native PCRE Regex engine.

"""
function regex_find(path::String, r::Regex)
  path = realpath(path)
  if Sys.iswindows()
    s = filter(x -> !(isempty(x) || x == path),
      String.(split(read(
          `powershell -Command "(Get-ChildItem -Path $path -File -Force -Recurse).FullName"`,
          String), "\r\n"))
              )
    s = [replace(i, Base.Filesystem.pathsep() => "/") for i in s]
    s2 = s[findall([occursin(r, f) for f in s])]

  else
    s = filter(x -> !(isempty(x) || x == path),
                String.(split(read(
                  `sh -c "find $path -type f"`,
                  String), "\n"))
              )

    s2 = String[]
    m = length(path) + 2
    for (i,f) in enumerate(s)
      s1 = f[m:end]
      if occursin(r, s1)
        push!(s2, f)
      end
    end
  end

  # Julia doesn't seem to handle regex searches in shell
  return sort(s2)
end

@doc """
    ls(str::String)

Similar functionality to Bash ls -1 with OS-agnostic output. Accepts wildcards.
Always returns full path and file name.

    ls()

Return full path and file name of files in current working directory.
""" ls
function ls(s::String)
  safe_isfile(s) && return [realpath(s)]
  safe_isdir(s) && return [joinpath(realpath(s), i) for i in readdir(s)]

  (p,f) = splitdir(s)

  if any([occursin(i, s) for i in regex_chars]) || occursin("*", p) || f == "*"
    # We're actually going to start at the highest-level directory that is
    # uniquely specified, so rather than starting at p from splitdir...
    fpat = String.(split(s, "*"))
    path, ff = splitdir(fpat[1])
    if isempty(ff)
      popfirst!(fpat)

      #= ...but this can leave us with an empty fpat and regex for that
      isn't standardized, so... =#
      if isempty(fpat) || fpat == [""]
        fpat = [".*"]
      end
    else
      fpat[1] = ff
    end

    # In case of empty path ... ?
    if isempty(path)
      path = "."
    end

    # So we're going to check for matches on all but the first m of each string:
    if Sys.iswindows()
        for i = 1:length(fpat)
            fpat[i] = replace(fpat[i], "." => "\\.")
        end
    end
    ff = join(fpat, ".*")
    mpat = Regex(ff * "\$")

    s1 = regex_find(path, mpat)
    if s1 == nothing
      s1 = []
    end
  else
    s1 = try
      glob(f,p)
    catch
      String[]
    end

    # DND DND DND DND DND
    # Control for odd behavior of glob in Linux
    if !(isempty(s1))
      if isempty(p); p = "." ; end
      for (i,s) in enumerate(s1)
        f = splitdir(s)[2]
        s1[i] = joinpath(realpath(p), f)
      end
    end
    # DND DND DND DND DND
  end
  return s1
end

ls() = [joinpath(realpath("."), i) for i in readdir(pwd())]


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

using Markdown

# ╔═╡ be8f2b17-1ffd-44a6-9932-dcf9850b200b
md"""
# Bucles

Temas:
1. bucles `while`
2. bucles `for`

## bucles while

La sintaxis para un `while` es

```julia
while *condición*
    *cuerpo de bucle*
end
```

Por ejemplo, usaremos un `while` para contar o iterar sobre un arreglo.
"""

# ╔═╡ 88e8cf59-9f29-473a-ad15-f9cb06bc8041
n = 0
while n < 10
    n += 1
    println(n)
end

# ╔═╡ 4592214b-7e99-44d2-a42b-bd089458e9e9
myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

i = 1
while i <= length(myfriends)
    friend = myfriends[i]
    println("Hola $friend, ¡gusto en verte!")
    i += 1
end

# ╔═╡ 057fa1b3-74bf-4fc6-8481-9e533fb6ab52
md"""
## bucles for

La sintaxis para un bucle `for` es

```julia
for *var* in *iterable de bucle*
    *cuerpo de bucle*
end
```

Podemos usar un bucle for para generar los datos anteriores
"""

# ╔═╡ fe7047d8-1adc-4530-b7bf-603f369bdff6
for n in 1:10
    println(n)
end

# ╔═╡ c6bed13a-3d1a-4aac-b194-ecaac14f81d8
myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

for friend in myfriends
    println("Hi $friend, it's great to see you!")
end

# ╔═╡ bf6fa9ad-2c42-4b36-a46f-8535b6d9b399
md"""
Nota: se puede reemplazar `in` con `=` ó `∈`.
"""

# ╔═╡ 8969e3d6-0bd5-4c8d-82ab-76601fdcc26f
for n = 1:10
    println(n)
end

# ╔═╡ c3bca3ca-810d-4f77-8d5b-47f81e0faab6
for n ∈ 1:10
    println(n)
end

# ╔═╡ 7e3053c0-c48a-48ca-aa49-d76d4b0ec5ef
md"""
Ahora usemos bucles `for` para crear tablas de sumas donde el valor de cada entrada es la suma de los índices del renglón y la columna.

Primero, inicializamos el arreglo con puros 0s.
"""

# ╔═╡ 392de431-f9ba-4b8f-96ef-a2b07ba40eba
m, n = 5, 5
A = fill(0, (m, n))

# ╔═╡ 7e9cc9af-875d-44c7-a579-0719b935524f
for i in 1:m
    for j in 1:n
        A[i, j] = i + j
    end
end
A

# ╔═╡ 5e68befc-4bba-449f-a584-00ae0686d98d
md"""
Aquí va un poquito de azúcar sintáctica para el bucle `for`
"""

# ╔═╡ 761a72b1-d015-49be-b096-49151298adc0
B = fill(0, (m, n))

# ╔═╡ 48b16d10-1566-4437-8cca-fced7a30b42f
for i in 1:m, j in 1:n
    B[i, j] = i + j
end
B

# ╔═╡ 1d481a45-af10-4ed7-a21c-9eb8a1bf17d6
md"""
La manera mas "Juliana" de crear esta tabla es por medio de un *array comprehension / comprehensión de arreglo*.
"""

# ╔═╡ c92fd9a3-5525-4046-a6bd-08323ad27999
C = [i + j for i in 1:m, j in 1:n]

# ╔═╡ 5a1aa04a-a98c-4292-80f2-b30d0018adfa
md"""
En el próximo ejemplo, embebimos un arreglo de comprehensión en un bucle `for`, generando tablas de adición de tamaño creciente.
"""

# ╔═╡ 66830444-9540-42ed-8fe2-5e6fb87cde62
for n in 1:10
   A = [i + j for i in 1:n, j in 1:n]
   display(A)
end

# ╔═╡ 78391d6a-3f6c-430b-9daf-480f7e12131d
md"""
### Ejercicios

4.1 Crea un diccionario `squares`, que tiene llaves de valores de 1 a 100. El valor asociado a cada llave es el cuadrado de la llave. Guarda los valores asociados a las llaves pares como enteros y las impares como cadenas. Por ejemplo,

```julia
squares[10] == 100
squares[11] == "121"
```

(¡No necesitas condicionales para esto!)
"""

# ╔═╡ bb3eb294-d78d-4e64-868e-84ae0f28c985

# ╔═╡ 44c740ff-80e0-4390-b599-00db419c269c
md"""
4.2 Usa `fill` para crea una matriz de `10x10` de solo `0`'s. Pobla las primeras 10 entradas con el índice de esa entrada. ¿Julia usa el orden de primero columna o primero renglón? (O sea, ¿el "segundo" elemento es el de la primera columna en el primer renglón, ó es el de el primer renglón en la segunda columna?)
"""

# ╔═╡ 074c5f35-f3cc-441f-87c4-b9c407b13efb

# ╔═╡ Cell order:
# ╟─be8f2b17-1ffd-44a6-9932-dcf9850b200b
# ╠═88e8cf59-9f29-473a-ad15-f9cb06bc8041
# ╠═4592214b-7e99-44d2-a42b-bd089458e9e9
# ╟─057fa1b3-74bf-4fc6-8481-9e533fb6ab52
# ╠═fe7047d8-1adc-4530-b7bf-603f369bdff6
# ╠═c6bed13a-3d1a-4aac-b194-ecaac14f81d8
# ╟─bf6fa9ad-2c42-4b36-a46f-8535b6d9b399
# ╠═8969e3d6-0bd5-4c8d-82ab-76601fdcc26f
# ╠═c3bca3ca-810d-4f77-8d5b-47f81e0faab6
# ╟─7e3053c0-c48a-48ca-aa49-d76d4b0ec5ef
# ╠═392de431-f9ba-4b8f-96ef-a2b07ba40eba
# ╠═7e9cc9af-875d-44c7-a579-0719b935524f
# ╟─5e68befc-4bba-449f-a584-00ae0686d98d
# ╠═761a72b1-d015-49be-b096-49151298adc0
# ╠═48b16d10-1566-4437-8cca-fced7a30b42f
# ╟─1d481a45-af10-4ed7-a21c-9eb8a1bf17d6
# ╠═c92fd9a3-5525-4046-a6bd-08323ad27999
# ╟─5a1aa04a-a98c-4292-80f2-b30d0018adfa
# ╠═66830444-9540-42ed-8fe2-5e6fb87cde62
# ╟─78391d6a-3f6c-430b-9daf-480f7e12131d
# ╠═bb3eb294-d78d-4e64-868e-84ae0f28c985
# ╟─44c740ff-80e0-4390-b599-00db419c269c
# ╠═074c5f35-f3cc-441f-87c4-b9c407b13efb

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
export get_fpga_name, get_fpga_state

"""
    get_fpga_name(; device="fpga0" )

Get the name of the FPGA mounted as `device`.
"""
function get_fpga_name(; device="fpga0" )
    is_valid_device( device )

    return readchomp( "$kmod_class_dir/$device/name" )
end

"""
    get_fpga_state(; device="fpga0 )

Get the current state of the FPGA mounted as `device`.
"""
function get_fpga_state(; device="fpga0" )
    is_valid_device( device )

    state = readchomp( "$kmod_class_dir/$device/state" )

    return StatusFromStrings[ state ]
end


text:
function count_bits(chars::AbstractVector{Char})
    c0 = count(chars .== '0')
    c1 = count(chars .== '1')
    c0, c1
end

function most_common_bit(chars::AbstractVector{Char})
    c0, c1 = count_bits(chars)
    c0 > c1 ? '0' : '1'
end

function least_common_bit(chars::AbstractVector{Char})
    c0, c1 = count_bits(chars)
    c0 <= c1 ? '0' : '1'
end

function oxygen_generator(chars::AbstractVector{Char})
    chars .== most_common_bit(chars)
end

function co2_scrubber(chars::AbstractVector{Char})
    chars .== least_common_bit(chars)
end

function rating(matrix::AbstractMatrix{Char}, bit_criteria)
    for bit_position in 1:size(matrix, 2)
        index = bit_criteria(matrix[:, bit_position])
        matrix = matrix[index, :]
        if size(matrix, 1) == 1
            return matrix[1, :]
        end
    end
end

function binary_number(chars::AbstractVector{Char})
    s = String(chars)
    parse(Int, s, base = 2)
end

matrix = vcat((permutedims(Vector{Char}(line)) for line in readlines())...)
gamma = binary_number([most_common_bit(col) for col in eachcol(matrix)])
epsilon = binary_number([least_common_bit(col) for col in eachcol(matrix)])
oxygen = binary_number(rating(matrix, oxygen_generator))
co2 = binary_number(rating(matrix, co2_scrubber))
println(gamma * epsilon)
println(oxygen * co2)


text:
@time using MonteCarlo, MonteCarloObservable, StableDQMC
#  3.567540 seconds (7.22 M allocations: 469.081 MiB, 4.08% gc time, 0.34% compilation time)
using Test
using LinearAlgebra, Random, Dates
using MonteCarlo: @bm, TimerOutputs

struct DummyModel <: MonteCarlo.Model
    lattice
end
DummyModel() = DummyModel(SquareLattice(2))

# check elementwise, not matrix norm
function check(A::Array, B::Array, atol, rtol=atol)
    for (x, y) in zip(A, B)
        if !isapprox(x, y, atol=atol, rtol=rtol)
            @info "$x ≉ $y "
            return false
        end
    end
    true
end
function check(x::Number, y::Number, atol, rtol=atol)
    result = isapprox(x, y, atol=atol, rtol=rtol)
    result || @info "$x ≉ $y"
    result
end

# In case some test failed and left behind a .jld file
for f in readdir()
    if endswith(f, ".jld") || endswith(f, "jld2")
        @warn "Removing $f"
        rm(f)
    end
end

macro benchmark_test(name, code)
    TimerOutputs.timer_expr(MonteCarlo, true, name, code)
end

@testset "All Tests" begin
    println("Utilities")
    @time @testset "Utilities" begin
        @bm function test1(x, y)
            sleep(x+y)
        end
        @bm test2(x, y) = sleep(x+y)
        # @eval MonteCarlo ... will put the body of test3, test4 in the wrong scope
        # eval(timer_expr(...)) doesn't work
        function test3(x, y)
            @benchmark_test "test1(::Any, ::Any)" begin sleep(x+y) end
        end
        test4(x, y) = @benchmark_test "test2(::Any, ::Any)" begin sleep(x+y) end

        x = code_lowered(test1, Tuple{Float64, Float64})[1]
        y = code_lowered(test3, Tuple{Float64, Float64})[1]
        @test x.code == y.code

        x = code_lowered(test2, Tuple{Float64, Float64})[1]
        y = code_lowered(test4, Tuple{Float64, Float64})[1]
        @test x.code == y.code

        @test !MonteCarlo.timeit_debug_enabled()
        enable_benchmarks()
        @test MonteCarlo.timeit_debug_enabled()
        disable_benchmarks()
        @test !MonteCarlo.timeit_debug_enabled()
    end
    # 0.473232 seconds (778.10 k allocations: 47.050 MiB, 2.72% gc time, 99.50% compilation time)

    println("Lattices")
    @time @testset "Lattices" begin
        include("lattices.jl")
    end
    # 26.765963 seconds (50.23 M allocations: 2.407 GiB, 3.35% gc time, 58.25% compilation time)

    println("Model")
    @time @testset "Models" begin
        include("modeltests_IsingModel.jl")
        include("modeltests_HubbardModel.jl")
    end
    # 11.418479 seconds (29.34 M allocations: 1.301 GiB, 4.78% gc time, 7.60% compilation time)

    println("DQMC")
    @time @testset "Flavors" begin
        # include("flavortests_MC.jl")
        include("flavortests_DQMC.jl")
    end
    # 61.075442 seconds (107.31 M allocations: 6.180 GiB, 4.44% gc time, 5.18% compilation time)

    println("Scheduler")
    @time @testset "Scheduler & (DQMC) Updates" begin
        include("updates.jl")
    end
    # 12.962446 seconds (17.87 M allocations: 1017.561 MiB, 2.51% gc time, 7.02% compilation time)

    println("Measurement")
    @time @testset "Measurements" begin
        include("measurements.jl")
    end
    # 13.129428 seconds (21.85 M allocations: 1.240 GiB, 4.31% gc time, 2.01% compilation time)

    println("Integration")
    @time @testset "Integration tests" begin
        include("integration_tests.jl")
    end
    # 20.755311 seconds (39.02 M allocations: 2.359 GiB, 3.11% gc time, 12.65% compilation time)

    println("ED")
    @time @testset "Exact Diagonalization" begin
        include("ED/ED_tests.jl")
    end
    # 61.415264 seconds (49.97 M allocations: 12.824 GiB, 2.79% gc time, 10.31% compilation time)

    println("File IO")
    @time @testset "File IO" begin
        include("FileIO.jl")
    end
    # 81.479393 seconds (120.31 M allocations: 5.953 GiB, 2.20% gc time, 5.36% compilation time)
end


text:
export expandPolygon, insidePolygon

function expandPolygon( x::Vector{Float64}, y::Vector{Float64},  # polygon points
	                     LL )  # distance
# Expand the polygon.  Each line segment goes out a distance of LL
# in the normal direction.

npts = length(x)

if length(y) != npts
   error("length(y) != npts")
end

xs = Array{Float64}(npts*2)
ys = Array{Float64}(npts*2)

idx = 1

for i = 1:npts-1

   theta = atan2(y[i+1] - y[i], x[i+1] - x[i])
   theta += pi/2  # rotate by 90 for normal direction.

   lct = LL * cos(theta)
   lst = LL * sin(theta)

   xs[idx] = x[i] + lct
   ys[idx] = y[i] + lst

   xs[idx+1] = x[i+1] + lct
   ys[idx+1] = y[i+1] + lst

   idx += 2
end # i

# end point
theta = atan2(y[1] - y[npts], x[1] - x[npts])
theta += pi/2

lct = LL * cos(theta)
lst = LL * sin(theta)

xs[idx] = x[npts] + lct
ys[idx] = y[npts] + lst

xs[idx+1] = x[1] + lct
ys[idx+1] = y[1] + lst

return xs, ys
end # function expandPolygon

#---------------------------------------------------------------------

function insidePolygon( x::Vector{Float64}, y::Vector{Float64},  # polygon points
  	                     xp,yp)
#  Is point (xp,yp) inside polygon defined by x,y ?

n = length(x)

if length(y) != n
   error("length(y) != n")
end

inside = false

for i = 1:n
   i1 = i + 1
   if i1 > n
      i1 = 1
   end
   dy = y[i1] - y[i]

   if dy != 0.0
      y1 = (yp <= y[i])
      y2 = (yp >  y[i1])

      if y1 == y2

         rslope = (x[i1] - x[i]) / dy
         d = (xp-x[i]) - (yp-y[i]) * rslope

         if abs(d) < 1.e-4
         	 d = 0.0
         end

         if d <= 0.0
            inside = !inside
         end

      end  # y1 == y2

   end # (dy != 0.0)

end # i

return inside
end # function insidePolygon


text:
"""
`init_subspec!(m::RealBond)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::RealBond)
    if subspec(m) == "ss0"
        return
    elseif subspec(m) == "ss1"
        return ss1!(m)
    elseif subspec(m) == "ss2"
        return ss2!(m)
    else
        error("This subspec should be a 0")
    end
end

"""
```
ss1!(m::RealBond)
```

Initializes subspec 1 of `RealBond`. In this subspecification,
The monetary shock process is shut down with σ_mon = 0., and
ρ_z is fixed so the only parameter that is free is σ_z.
"""
function ss1!(m::RealBond)
    m <= parameter(:ρ_z, 0.0, fixed = true,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")
    m <= parameter(:ρmon, 0.0, fixed = true,
                   description="ρmon: AR(1) coefficient in the monetary policy process.",
                   tex_label="\\rho_z")
    m <= parameter(:σ_z, sqrt(.007), (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_z: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label="\\sigma_{z}")
    m <= parameter(:σ_mon, 0.0, fixed = true,
                   description="σ_mon: The standard deviation of the monetary policy shock.",
                   tex_label="\\sigma_{mon}")
end

"""
```
ss2!(m::RealBond)
```

Initializes subspec 2 of `RealBond`. In this subspecification,
The technology shock process is shut down with σ_z = 0., and
ρ_z is fixed so the only parameter that is free is σ_mon.
"""
function ss2!(m::RealBond)
    m <= parameter(:ρ_z, 0.0, fixed = true,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")
    m <= parameter(:ρmon, 0.0, fixed = true,
                   description="ρmon: AR(1) coefficient in the monetary policy process.",
                   tex_label="\\rho_z")
    m <= parameter(:σ_z, 0.0, fixed = true,
                   description="σ_mon: The standard deviation of the monetary policy shock.",
                   tex_label="\\sigma_{mon}")
    m <= parameter(:σ_mon, 0.2380, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_mon: The standard deviation of the monetary policy shock.",
                   tex_label="\\sigma_{mon}")
end


text:
"""
    orthiter(;input::AbstractArray=[], outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractArray=[], pseudocount::AbstractArray=[], colmeanlist::AbstractArray="", colvarlist::Union{Nothing,AbstractString,AbstractArray}="", dim::Number=3, numepoch::Number=3, lower::Number=0, upper::Number=1.0f+38, evalfreq::Number=5000, offsetStoch::Number=1f-20, initW::Union{Nothing,AbstractArray}=nothing, logdir::Union{Nothing,AbstractString}=nothing, perm::Bool=false)
Online CCA solved by orthogonal iteraction (power method).
Input Arguments
---------
- `input` : Julia Binary file generated by `OnlineCCA.csv2bin` function.
- `outdir` : The directory specified the directory you want to save the result.
- `scale` : {log,ftt,raw}-scaling of the value.
- `pseudocount` : The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated.
- `colmeanlist` : The mean of each column of matrix. The CSV file is generated by `OnlineCCA.sumr` functions.
- `colvarlist` : The variance of each column of matrix. The CSV file is generated by `OnlineCCA.sumr` functions.
- `dim` : The number of dimension of CCA.
- `numepoch` : The number of epoch.
- `lower` : Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)
- `upper` : Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)
- `evalfreq` : Evaluation Frequency of Reconstruction Error
- `offsetStoch` : Off set value for avoding overflow when calculating stochastic gradient
- `initW` : The CSV file saving the initial values of eigenvectors.
- `logdir` : The directory where intermediate files are saved, in every evalfreq (e.g. 5000) iteration.
- `perm` : Whether the data matrix is shuffled at random.
Output Arguments
---------
- `W` : Eigen vectors of covariance matrix (No. columns of the data matrix × dim)
- `λ` : Eigen values (dim × dim)
- `V` : Latent variables of covariance matrix (No. rows of the data matrix × dim)
- `ExpCorVar` : Explained correlation/variance by the eigenvectors
- `TotalCorVar` : Total correlation/variance of the data matrix
- `Objective` : The summation of object function (e.g. w_i' X_i' X_j w_j)
- `AvgLatentCor` : The mean of correlation efficient of latent variables (e.g. cor(X_i w_i, X_j w_j))
- stop : Whether the calculation is converged
"""
function orthiter(;input::AbstractArray=[], outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractArray=[], pseudocount::AbstractArray=[], colmeanlist::AbstractArray="", colvarlist::Union{Nothing,AbstractString,AbstractArray}="", dim::Number=3, numepoch::Number=3, lower::Number=0, upper::Number=1.0f+38, evalfreq::Number=5000, offsetStoch::Number=1f-20, initW::Union{Nothing,AbstractArray}=nothing, logdir::Union{Nothing,AbstractString}=nothing, perm::Bool=false)
    # Initial Setting
    cca = ORTHITER()
    pseudocount, W, v, D, colmeanlist, colmeanvec, colvarlist, colvarvec, N, M, TotalCorVar, lower, upper, evalfreq, offsetStoch = init(input, pseudocount, dim, colmeanlist, colvarlist, initW, logdir, cca, lower, upper, evalfreq, offsetStoch, scale)
    # Perform CCA
    out = orthiter(input, outdir, scale, pseudocount, colmeanlist, colvarlist, dim, numepoch, logdir, cca, W, v, D, colmeanvec, colvarvec, N, M, TotalCorVar, lower, upper, evalfreq, offsetStoch, perm)
    if outdir isa String
        output(outdir, out, cca)
    end
    return out
end

function orthiter(input, outdir, scale, pseudocount, colmeanlist, colvarlist, dim, numepoch, logdir, cca, W, v, D, colmeanvec, colvarvec, N, M, TotalCorVar, lower, upper, evalfreq, offsetStoch, perm)
    N, M = nm(input)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    l = length(input)
    x = Any[]
    normx = Any[]
    for i=1:l
        push!(x, zeros(UInt32, M[i]))
        push!(normx, zeros(Float32, M[i]))
    end
    # If not 0 the calculation is converged
    stop = 0
    s = 1
    n = 1
    stream = Any[]
    # Each epoch s
    progress = Progress(numepoch*N[1])
    while(stop == 0 && s <= numepoch)
        next!(progress)
        if n == 1
            stream = Any[]
            # Stream
            for i=1:l
                push!(stream, ZstdDecompressorStream(open(input[i])))
            end
            # tmpN, tmpN
            for i=1:l
                read!(stream[i], tmpN)
                read!(stream[i], tmpM)
            end
        else
            # x, norm
            for i=1:l
                read!(stream[i], x[i])
                normx[i] = normalizex(x[i], n, stream[i], scale[i], pseudocount[i], colmeanlist[i], colmeanvec[i], colvarlist[i], colvarvec[i])
                if perm
                    normx[i] .= normx[randperm(length(normx[i]))]
                end
            end
            # Update Eigen vector
            for i=1:l
                for k=1:l
                    if i != k
                        W[i] .+= normx[i] * (normx[k]'W[k] * offsetStoch)
                    end
                end
            end
            # close
            if n == N[1]
                for i=1:l
                    close(stream[i])
                end
            end
        end
        n += 1
        if n == N[1] + 1
            # NaN
            checkNaN(W, GD())
            # Retraction
            W .= Retraction(W)
            # save log file
            if logdir isa String
                stop = outputlog(s, input, dim, logdir, W, GD(), TotalCorVar, scale, pseudocount, colmeanlist, colmeanvec, colvarlist, colvarvec, lower, upper, stop)
            end
            s += 1
            n = 1
        end
    end

    # Return, W, λ, V
    out = WλV(W, input, dim, scale, pseudocount, colmeanlist, colmeanvec, colvarlist, colvarvec, TotalCorVar)
    return (out[1], out[2], out[3], out[4], out[5], out[6], out[7], stop)
end


text:
using EngEconomics

# Given
xInit = 120000
yInit = 96000
xBenefits1 = (11000, 10)
xBenefits2 = (9000, 20)
yBenefits = (12000, 20)
xSalvage = 40000
ySalvage = 20000
MARR = 0.08
N = 20

# Find Present Worth Method Better Option
xInitPW = -xInit
xBenefitsPW = xBenefits1[1] * seriesPresentAmountFactor(MARR, xBenefits1[2]) +
				xBenefits2[1] * seriesPresentAmountFactor(MARR, xBenefits2[2])
xSalvagePW = xSalvage * presentWorthFactor(MARR, N)
xPW = xInitPW + xBenefitsPW + xSalvagePW

yInitPW = -yInit
yBenefitsPW = yBenefits[1] * seriesPresentAmountFactor(MARR, yBenefits[2])
ySalvagePW = ySalvage * presentWorthFactor(MARR, N)
yPW = yInitPW + yBenefitsPW + ySalvagePW

if xPW > yPW
	println("X is better")
else
	println("Y is better")
end


text:
using Gadfly

function pf(porder,mesh,a,b,Q)
    basis = phi(porder+1)
    function (x)
        k = min(1+div(x-a,(b-a)/K),K)
        poly = evaluate(Q[1,k,:],basis,𝜒(mesh,k,x))
    end
end

function pf(p::DG1D,Q)
    basis = phi(p)
    function (x)
        k = min(1+div(x-p.a,(p.b-p.a)/p.K),p.K)
        poly = evaluate(Q[1,k,:],basis,𝜒(p.mesh,k,x))
    end
end

plotSolution(p::DG1D,Q) = Gadfly.plot(pf(p,Q),p.a,p.b)


text:
@time begin
@time begin
#    function fall()
    function f(x)
        return x
    end

    function y(h, x)
        return h(x)
    end

    println("f")
    @time f(1.0)
    @time f(1.0)
    println("y")
    @time y(f, 1.0)
    @time y(f, 1.0)

    println("f broad")
    @time f.([1.0, 2.0])
    @time f.([1.0, 2.0])
    println("y broad")
    @time y.(f, [1.0, 2.0])
    @time y.(f, [1.0, 2.0])

    function g(x)
        return 2*x
    end

    println("g g broad")
    @time g(1.0)
    @time g.([1.0, 2.0])

    println("y g broad")
    @time y.(g, [1.0, 2.0])
    @time y.(g, [1.0, 2.0])
    println()
    println("done")
#    end
#
#    @time fall()
#    println("done fall1")
#    @time fall()
#    println("done fall2")

    end

end


text:
using Test
using TestSetExtensions
using SafeTestsets

@time @safetestset "greedy coverage" begin include("test_coverage.jl") end

max_stars_repo_name: adolgert/test-scientific-code

max_issues_repo_name: adolgert/test-scientific-code

max_forks_repo_name: adolgert/test-scientific-code


text:
using TransportBasedInference: evaluate

@testset "Test inverse function" begin

    Ne = 200
    Nx = 3

    idx = [0 0 0; 0 0 1; 0 0 2; 0 1 0; 0 0 1; 0 1 1; 0 2 1; 0 3 1; 1 3 1]

    Nψ = size(idx,1)

    coeff = randn(Nψ)

    B = MultiBasis(CstProHermiteBasis(6), Nx)
    f = ExpandedFunction(B, idx, coeff)
    R = IntegratedFunction(f)

    X = randn(Nx, Ne) .* randn(Nx, Ne)

    F = evaluate(R, X)

    Xmodified = deepcopy(X)

    Xmodified[end,:] .+= cos.(2*π*randn(Ne)) .* exp.(-randn(Ne).^2/2)

    @test norm(Xmodified[end,:]-X[end,:])>1e-6

    S = Storage(f, Xmodified)

    inverse!(Xmodified, F, R, S)

    @test norm(Xmodified[end,:] - X[end,:])<2e-8

    @test norm(Xmodified[1:end-1,:] - X[1:end-1,:])<2e-8
end

@testset "Test inverse HermiteHermiteMapComponent" begin

    Ne = 200
    Nx = 3

    idx = [0 0 0; 0 0 1; 0 0 2; 0 1 0; 0 0 1; 0 1 1; 0 2 1; 0 3 1; 1 3 1]

    Nψ = size(idx,1)

    coeff = randn(Nψ)

    C = HermiteMapComponent(20, Nx, idx, coeff)

    X = randn(Nx, Ne) .* randn(Nx, Ne)

    F = evaluate(C.I, X)

    Xmodified = deepcopy(X)

    Xmodified[end,:] .+= cos.(2*π*randn(Ne)) .* exp.(-randn(Ne).^2/2)

    @test norm(Xmodified[end,:]-X[end,:])>1e-6

    S = Storage(C.I.f, Xmodified)

    inverse!(Xmodified, F, C, S)

    @test norm(Xmodified[end,:] - X[end,:])<1e-5

    @test norm(Xmodified[1:end-1,:] - X[1:end-1,:])<1e-5
end

@testset "Test inverse LinHermiteHermiteMapComponent" begin

    Ne = 200
    Nx = 3

    idx = [0 0 0; 0 0 1; 0 0 2; 0 1 0; 0 0 1; 0 1 1; 0 2 1; 0 3 1; 1 3 1]

    Nψ = size(idx,1)

    coeff = randn(Nψ)

    C = HermiteMapComponent(20, Nx, idx, coeff)
    X = randn(Nx, Ne) .* randn(Nx, Ne)

    X0 = deepcopy(X)
    Xmodified = deepcopy(X)

    L = LinHermiteMapComponent(X, C)
    F = evaluate(L, X)

    @test norm(X -X0)<1e-8

    Xmodified[end,:] .+= cos.(2*π*randn(Ne)) .* exp.(-randn(Ne).^2/2)

    Lmodified = LinHermiteMapComponent(Xmodified, C)

    @test norm(Xmodified[end,:]-X[end,:])>1e-6

    transform!(Lmodified.L, Xmodified)
    S = Storage(Lmodified.C.I.f, Xmodified);
    itransform!(Lmodified.L, Xmodified)

    inverse!(Xmodified, F, L, S)

    @test norm(Xmodified[end,:] - X[end,:])<1e-5

    @test norm(Xmodified[1:end-1,:] - X[1:end-1,:])<1e-5
end


text:
@testset "Primitives" begin
  @testset "Boxes" begin
    b = Box(P2(0, 0), P2(1, 1))
    @test embeddim(b) == 2
    @test paramdim(b) == 2
    @test coordtype(b) == T
    @test minimum(b) == P2(0, 0)
    @test maximum(b) == P2(1, 1)
    @test extrema(b) == (P2(0, 0), P2(1, 1))

    @test measure(b) == T(1)
    @test P2(1, 1) ∈ b
  end

  @testset "Ball" begin
    b = Ball(P3(1,2,3), T(5))
    @test embeddim(b) == 3
    @test paramdim(b) == 3
    @test coordtype(b) == T
    @test center(b) == P3(1,2,3)
    @test radius(b) == T(5)

    b = Ball(P2(0,0), T(2))
    @test measure(b) ≈ π*(2^2)
    b = Ball(P3(0,0,0), T(2))
    @test measure(b) ≈ (4/3)*π*(2^3)

    b = Ball(P2(0,0), T(2))
    @test P2(1,0) ∈ b
    @test P2(0,1) ∈ b
    @test P2(2,0) ∈ b
    @test P2(0,2) ∈ b
    @test P2(3,5) ∉ b

    b = Ball(P3(0,0,0), T(2))
    @test P3(1,0,0) ∈ b
    @test P3(0,0,1) ∈ b
    @test P3(2,0,0) ∈ b
    @test P3(0,0,2) ∈ b
    @test P3(3,5,2) ∉ b
  end

  @testset "Sphere" begin
    s = Sphere(P3(0,0,0), T(1))
    @test embeddim(s) == 3
    @test paramdim(s) == 2
    @test coordtype(s) == T
    @test center(s) == P3(0, 0, 0)
    @test radius(s) == T(1)

    s = Sphere(P2(0,0), T(2))
    @test measure(s) ≈ 2π*2
    s = Sphere(P3(0,0,0), T(2))
    @test measure(s) ≈ 4π*(2^2)

    s = Sphere(P2(0,0), T(2))
    @test P2(1,0) ∉ s
    @test P2(0,1) ∉ s
    @test P2(2,0) ∈ s
    @test P2(0,2) ∈ s
    @test P2(3,5) ∉ s

    s = Sphere(P3(0,0,0), T(2))
    @test P3(1,0,0) ∉ s
    @test P3(0,0,1) ∉ s
    @test P3(2,0,0) ∈ s
    @test P3(0,0,2) ∈ s
    @test P3(3,5,2) ∉ s
  end

  @testset "Cylinder" begin
    c = Cylinder(P3(1,2,3), P3(4,5,6), T(5))
    @test embeddim(c) == 3
    @test paramdim(c) == 3
    @test coordtype(c) == T
    @test radius(c) == T(5)
    @test height(c) ≈ √27

    @test measure(c) ≈ π*5.0^2*√27
  end
end


text:
using ITensors, AutoHOOT, Zygote, OptimKit
using Zygote: @adjoint

const itensorad = AutoHOOT.ITensorsAD
const ad = AutoHOOT.autodiff
const gops = AutoHOOT.graphops

"""Generate an array of networks representing inner products, <p|H_1|p>, ..., <p|H_n|p>, <p|p>
Parameters
----------
peps: a peps network with datatype PEPS
peps_prime: prime of peps used for inner products
peps_prime_ham: prime of peps used for calculating expectation values
Hlocal: An array of MPO operators with datatype LocalMPO
Returns
-------
An array of networks.
"""
function generate_inner_network(
    peps::PEPS,
    peps_prime::PEPS,
    peps_prime_ham::PEPS,
    Hlocal::Array,
)
    network_list = []
    for H_term in Hlocal
        inner = inner_network(
            peps,
            peps_prime,
            peps_prime_ham,
            H_term.mpo,
            [H_term.coord1, H_term.coord2],
        )
        network_list = vcat(network_list, [inner])
    end
    inner = inner_network(peps, peps_prime)
    network_list = vcat(network_list, [inner])
    return network_list
end

# gradient of this function returns nothing.
@adjoint function generate_inner_network(
    peps::PEPS,
    peps_prime::PEPS,
    peps_prime_ham::PEPS,
    Hlocal::Array,
)
    adjoint_pullback(v) = (nothing, nothing, nothing, nothing)
    return generate_inner_network(peps, peps_prime, peps_prime_ham, Hlocal),
    adjoint_pullback
end

function rayleigh_quotient(inners::Array)
    self_inner = itensorad.scalar(inners[length(inners)])
    expectations = itensorad.scalar(sum(inners[1:length(inners)-1]))
    return expectations / self_inner
end

function loss_grad_wrap(peps::PEPS, Hlocal::Array)
    function loss(peps::PEPS)
        peps_prime = prime(peps; ham = false)
        peps_prime_ham = prime(peps; ham = true)
        network_list = generate_inner_network(peps, peps_prime, peps_prime_ham, Hlocal)
        variables = extract_data([peps, peps_prime, peps_prime_ham])
        inners = itensorad.batch_tensor_contraction(network_list, variables...)
        return rayleigh_quotient(inners)
    end
    loss_w_grad(peps::PEPS) = loss(peps), gradient(loss, peps)[1]
    return loss_w_grad
end

"""Update PEPS based on gradient descent
Parameters
----------
peps: a peps network with datatype PEPS
Hlocal: An array of MPO operators with datatype LocalMPO
stepsize: step size used in the gradient descent
num_sweeps: number of gradient descent sweeps/iterations
Returns
-------
An array containing Rayleigh quotient losses after each iteration.
"""
function gradient_descent(peps::PEPS, Hlocal::Array; stepsize::Float64, num_sweeps::Int)
    loss_w_grad = loss_grad_wrap(peps, Hlocal)
    # gradient descent iterations
    losses = []
    for iter = 1:num_sweeps
        l, g = loss_w_grad(peps)
        print("The rayleigh quotient at iteraton $iter is $l\n")
        peps = peps - stepsize * g
        push!(losses, l)
    end
    return losses
end

function optimize(peps::PEPS, Hlocal::Array; num_sweeps::Int, method = "GD")
    @assert(method in ["GD", "LBFGS", "CG"])
    inner(x, peps1, peps2) = peps1 * peps2
    loss_w_grad = loss_grad_wrap(peps, Hlocal)
    scale(peps, alpha) = alpha * peps
    add(peps1, peps2, alpha) = peps1 + alpha * peps2
    linesearch = HagerZhangLineSearch()
    if method == "GD"
        alg = GradientDescent(num_sweeps, 1e-8, linesearch, 2)
    elseif method == "LBFGS"
        alg = LBFGS(
            16;
            maxiter = num_sweeps,
            gradtol = 1e-8,
            linesearch = linesearch,
            verbosity = 2,
        )
    elseif method == "CG"
        alg = ConjugateGradient(;
            maxiter = num_sweeps,
            gradtol = 1e-8,
            linesearch = linesearch,
            verbosity = 2,
        )
    end
    _, _, _, _, history =
        OptimKit.optimize(loss_w_grad, peps, alg; inner = inner, scale! = scale, add! = add)
    return history[:, 1]
end


text:
module TestInfo

using MLJBase
import MLJBase.info_dict
using Test
using OrderedCollections
using LossFunctions

mutable struct DummyProb <: Probabilistic
    an_int::Int
    a_float::Float64
    a_vector::Vector{Float64}
    untyped
end

MLJBase.load_path(::Type{DummyProb})        = "GreatPackage.MLJ.DummyProb"
MLJBase.input_scitype(::Type{DummyProb})    = MLJBase.Table(Finite)
MLJBase.target_scitype(::Type{DummyProb})   = AbstractVector{<:Continuous}
MLJBase.is_pure_julia(::Type{DummyProb})    = true
MLJBase.supports_weights(::Type{DummyProb}) = true
MLJBase.package_name(::Type{DummyProb})     = "GreatPackage"
MLJBase.package_uuid(::Type{DummyProb}) = "6f286f6a-111f-5878-ab1e-185364afe411"
MLJBase.package_url(::Type{DummyProb})     = "https://mickey.mouse.org"
MLJBase.package_license(::Type{DummyProb}) = "MIT"
MLJBase.hyperparameter_ranges(::Type{DummyProb}) =
          (range(Int, :an_int, values=[1,2]),
           range(Float64, :a_float, lower=1, upper=2),
           range(Vector{Float64}, :a_vector, values=[[1.0], [2.0]]),
           nothing)
MLJBase.predict(::DummyProb, fr, X) = nothing

mutable struct DummyDeterm <: Deterministic end

MLJBase.load_path(::Type{DummyDeterm})        = "GreatPackage.MLJ.DummyDeterm"
MLJBase.input_scitype(::Type{DummyDeterm})    = MLJBase.Table(Finite)
MLJBase.target_scitype(::Type{DummyDeterm})   = AbstractVector{<:Continuous}
MLJBase.is_pure_julia(::Type{DummyDeterm})    = true
MLJBase.supports_weights(::Type{DummyDeterm}) = true
MLJBase.package_name(::Type{DummyDeterm})     = "GreatPackage"
MLJBase.package_uuid(::Type{DummyDeterm}) = "6f286f6a-111f-5878-ab1e-185364afe411"
MLJBase.package_url(::Type{DummyDeterm}) = "https://mickey.mouse.org"
MLJBase.package_license(::Type{DummyDeterm}) = "MIT"
MLJBase.predict(::DummyDeterm, fr, X) = nothing

mutable struct DummyInt <: Interval end
MLJBase.load_path(::Type{DummyInt}) = "GreatPackage.MLJ.DummyInt"
MLJBase.input_scitype(::Type{DummyInt}) = MLJBase.Table(Finite)
MLJBase.target_scitype(::Type{DummyInt}) = AbstractVector{<:Continuous}
MLJBase.is_pure_julia(::Type{DummyInt}) = true
MLJBase.supports_weights(::Type{DummyInt}) = true
MLJBase.package_name(::Type{DummyInt}) = "GreatPackage"
MLJBase.package_uuid(::Type{DummyInt}) = "6f286f6a-111f-5878-ab1e-185364afe411"
MLJBase.package_url(::Type{DummyInt}) = "https://mickey.mouse.org"
MLJBase.package_license(::Type{DummyInt}) = "MIT"
MLJBase.predict(::DummyInt, fr, X) = nothing

mutable struct DummyUnsup <: Unsupervised end
MLJBase.load_path(::Type{DummyUnsup}) = "GreatPackage.MLJ.DummyUnsup"
MLJBase.input_scitype(::Type{DummyUnsup}) = MLJBase.Table(Finite)
MLJBase.output_scitype(::Type{DummyUnsup}) = AbstractVector{<:Continuous}
MLJBase.is_pure_julia(::Type{DummyUnsup}) = true
MLJBase.supports_weights(::Type{DummyUnsup}) = true
MLJBase.package_name(::Type{DummyUnsup}) = "GreatPackage"
MLJBase.package_uuid(::Type{DummyUnsup}) = "6f286f6a-111f-5878-ab1e-185364afe411"
MLJBase.package_url(::Type{DummyUnsup}) = "https://mickey.mouse.org"
MLJBase.package_license(::Type{DummyUnsup}) = "MIT"
MLJBase.transform(::DummyUnsup, fr, X) = nothing

@testset "info on probabilistic models" begin
    d = LittleDict{Symbol,Any}(
            :name             => "DummyProb",
            :load_path        => "GreatPackage.MLJ.DummyProb",
            :is_pure_julia    => true,
            :package_uuid     => "6f286f6a-111f-5878-ab1e-185364afe411",
            :package_name     => "GreatPackage",
            :package_license  => "MIT",
            :input_scitype    => MLJBase.Table(Finite),
            :output_scitype   => Unknown,
            :supports_weights => true,
            :supports_online  => false,
            :target_scitype   => MLJBase.AbstractVector{<:Continuous},
            :prediction_type  => :probabilistic,
            :package_url      => "https://mickey.mouse.org",
            :is_supervised    => true,
            :is_wrapper       => false,
            :docstring        => "DummyProb from GreatPackage.jl.\n[Documentation](https://mickey.mouse.org).",
            :implemented_methods  => [:predict, ],
            :hyperparameter_types => ("Int64", "Float64",
                                 "Array{Float64,1}", "Any"),
            :hyperparameters  => (:an_int, :a_float, :a_vector, :untyped),
            :hyperparameter_ranges =>
                (range(Int, :an_int, values=[1,2]),
                 range(Float64, :a_float, lower=1, upper=2),
                 range(Vector{Float64}, :a_vector, values=[[1.0], [2.0]]),
                 nothing))
    @test info_dict(DummyProb) == d
    @test info_dict(DummyProb(42, 3.14, [1.0, 2.0], :cow)) == d
end

@testset "info on deterministic models" begin
    d = LittleDict{Symbol,Any}(
            :name             => "DummyDeterm",
            :load_path        => "GreatPackage.MLJ.DummyDeterm",
            :is_pure_julia    => true,
            :package_uuid     => "6f286f6a-111f-5878-ab1e-185364afe411",
            :package_name     => "GreatPackage",
            :package_license  => "MIT",
            :input_scitype    => MLJBase.Table(Finite),
            :output_scitype   => Unknown,
            :supports_weights => true,
            :supports_online  => false,
            :target_scitype   => MLJBase.AbstractVector{<:Continuous},
            :prediction_type  => :deterministic,
            :package_url      => "https://mickey.mouse.org",
            :is_supervised    => true,
            :is_wrapper       => false,
            :docstring        => "DummyDeterm from GreatPackage.jl.\n[Documentation](https://mickey.mouse.org).",
            :implemented_methods   => [:predict, ],
            :hyperparameter_types  => (),
            :hyperparameters       => (),
            :hyperparameter_ranges => ())

    @test info_dict(DummyDeterm)   == d
    @test info_dict(DummyDeterm()) == d
end

@testset "info on interval models" begin
    d = LittleDict{Symbol,Any}(
            :name => "DummyInt",
            :load_path => "GreatPackage.MLJ.DummyInt",
            :is_pure_julia => true,
            :package_uuid  => "6f286f6a-111f-5878-ab1e-185364afe411",
            :package_name  => "GreatPackage",
            :package_license => "MIT",
            :input_scitype => MLJBase.Table(Finite),
            :output_scitype   => Unknown,
            :supports_weights => true,
            :supports_online => false,
            :target_scitype => MLJBase.AbstractVector{<:Continuous},
            :prediction_type => :interval,
            :package_url   => "https://mickey.mouse.org",
            :is_supervised => true,
            :is_wrapper => false,
            :docstring => "DummyInt from GreatPackage.jl.\n[Documentation](https://mickey.mouse.org).",
            :implemented_methods => [:predict, ],
            :hyperparameter_types  => (),
            :hyperparameters => (),
            :hyperparameter_ranges => ())

    @test info_dict(DummyInt)   == d
    @test info_dict(DummyInt()) == d
end

@testset "info on unsupervised models" begin
    d = LittleDict{Symbol,Any}(
            :name            => "DummyUnsup",
            :load_path       => "GreatPackage.MLJ.DummyUnsup",
            :is_pure_julia   => true,
            :package_uuid    => "6f286f6a-111f-5878-ab1e-185364afe411",
            :package_name    => "GreatPackage",
            :package_license => "MIT",
            :input_scitype   => MLJBase.Table(Finite),
            :target_scitype   => Unknown,
            :supports_weights => true,
            :prediction_type  => :unknown,
            :output_scitype  => MLJBase.AbstractVector{<:Continuous},
            :package_url     => "https://mickey.mouse.org",
            :is_supervised   => false,
            :supports_online => false,
            :is_wrapper      => false,
            :docstring       => "DummyUnsup from GreatPackage.jl.\n[Documentation](https://mickey.mouse.org).",
            :implemented_methods   => [:transform, ],
            :hyperparameter_types  => (),
            :hyperparameters       => (),
            :hyperparameter_ranges => ())

    @test info_dict(DummyUnsup)   == d
    @test info_dict(DummyUnsup()) == d
end

@testset "info for measures" begin
    @test info(rms).name == "RootMeanSquaredError"
end

end
true


text:
using Printf
# ARRAYS ARE INDEXED FROM ONE NOT ZERO IN JULIA!!!!

## I wrote this for displaying arrays, but we don't actually need it. println() works just fine!
#
#function display_array(some_array)
#    for x in some_array
#        print(x, " ")
#    end
#    println()
#end

function array_examples()
    println("--------------------------------------")
    println("Array examples")

    println()

    println("We can create one-dimensional arrays of zeros")
    one_dimensional_array = zeros(Int32, 5)
    println(one_dimensional_array)

    println("We can create two-dimensional arrays of zeros")
    two_dimensional_array = zeros(Int32, 5, 2)
    println(two_dimensional_array)

    println("We can create arrays where we specify the type, the size, and the initial value (in this case undefined)")
    undefined_array = Array{Int32}(undef, 5)
    println(undefined_array)

    println("We can create a int array and set initial values")
    int_array = Int32[1,2,3,4]
    println(int_array)

    println("We can get individual values. REMEMBER JULIA INDEXED FROM ONE NOT ZERO!")
    println(int_array[2]) # The second element really is the second element (not the third as would be the case with zero-indexing)

    println("We can create a float array and set initial values")
    float_array = Float64[1,2,3,4]
    println(float_array)

    println("We can create an array where each element is off a different type")
    array_of_different_types = [99, 3.1415, "foo", true, 'x']
    println(array_of_different_types)

    println("We can create an array from a range (2 to 10, stepping at intervals of 2")
    even_number_array = collect(2:2:10)
    println(even_number_array)

    println("We can create an array using comprehension")
    square_numbers = [n^2 for n in 1:5]
    println(square_numbers)

    println("We can also create multi-dimensional arrays using comprehension")
    another_multi_dim_array = [n * m for n in 1:5, m in 1:5]
    println(another_multi_dim_array)

    println("We can get a specific (one-indexed) element from the array")
    println(array_of_different_types[1])

    println("We can get the last element from the array")
    println(array_of_different_types[end])

    println("We can get a range, in this case items 2-4")
    println(array_of_different_types[2:4])

    println("We can check if an item, e.g. 3.1415 exists in the array")
    println(3.1415 in array_of_different_types)

    println("We can find a specific item in the array, e.g. 'x'")
    println(findfirst(isequal('x'), array_of_different_types))

    int_array = [1,2,3,4,5,6,7,8,10]
    println("We can find all the even numbers in an array using a generic function and findall")
    f(a) = mod(a, 2) == 0 ? true : false
    println(findall(f, int_array))

    println("We can get the number of items that pass the test")
    println(count(f, int_array))

    println("We can get the row and column size of our 5x2 array we created earlier")
    println(size(two_dimensional_array))

    println("We can get the length of an array")
    println(length(int_array))

    println("We can get the sum of all the values in the array")
    println(sum(int_array))

    println("We can insert items at index 2 (second number after colon seems to specify range in target array to overwrite. increase this number to see old items in array go missing!)")
    splice!(int_array, 2:1, [100,200,300,400,500])
    println(int_array)

    println("We can remove items at index 2 to 6")
    splice!(int_array, 2:6)
    println(int_array)

    println("We get the maximum value in the array")
    println(maximum(int_array))

    println("We get the minimum value in the array")
    println(minimum(int_array))

    println("We can perform calculations on an array without looping")
    println(int_array * 2)

    println("We can store functions in an array")
    function_array = [sin, cos, tan]
    println(function_array)

    for n in function_array
        println(n(45))
    end

    two_dimensional_int_array = [1 2 3 4; 5 6 7 8; 9 10 11 12]
    for row = 1:3, col = 1:4
        @printf("row %d, col %d = %d\n", row, col, two_dimensional_int_array[row,col])
    end
    println()

    println("We can print just the items from column 2")
    println(two_dimensional_int_array[:, 2])

    println("We can push items to the end of an array")
    push!(int_array, 999)
    println(int_array)
end

max_stars_repo_name: James-P-D/JuliaDump

max_issues_repo_name: James-P-D/JuliaDump

max_forks_repo_name: James-P-D/JuliaDump


text:
using BinaryBuilder, Pkg

name = "Bonmin"
version = v"1.8.8"

sources = [
    GitSource("https://github.com/coin-or/Bonmin.git",
              "65c56cea1e7c40acd9897a2667c11f91d845bb7b"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/Bonmin

if [[ ${target} == *mingw* ]]; then
    sed -i s/dllimport/dllexport/ /workspace/destdir/include/coin-or/IpoptConfig.h
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/fix_dllexport.patch
fi

# Remove misleading libtool files
rm -f ${libdir}/*.la
update_configure_scripts

# old and custom autoconf
sed -i s/elf64ppc/elf64lppc/ configure

export CPPFLAGS="${CPPFLAGS} -I${includedir} -I${includedir}/coin -I${includedir}/coin-or"
export CXXFLAGS="${CXXFLAGS} -std=c++11"

if [[ ${target} == *mingw* ]]; then
    export LDFLAGS="-L${libdir}"
    export LT_LDFLAGS="-no-undefined"
fi

./configure \
    --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    --disable-pthread-mumps \
    --enable-static \
    --enable-shared \
    lt_cv_deplibs_check_method=pass_all \
    --with-asl-lib="-lipoptamplinterface -lasl"

make
make install
"""

platforms = supported_platforms()
filter!(!Sys.isfreebsd, platforms)
platforms = expand_cxxstring_abis(platforms)
filter!(x -> cxxstring_abi(x) != "cxx03", platforms)
platforms = expand_gfortran_versions(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libbonmin", :libbonmin),
    ExecutableProduct("bonmin", :amplexe),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Cbc_jll", v"2.10.5"),
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("Ipopt_jll", v"3.13.4"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: pazner/Yggdrasil

max_issues_repo_name: pazner/Yggdrasil

max_forks_repo_name: pazner/Yggdrasil


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

@test integrate(x^3*arcsech(a+b*x), x) == :(1//4*x^4*arcsech(a+b*x)+-1//4*a^4*b^-4*arcsech(a+b*x)+(1/2)*a*b^-4*(1+2*a^2)*arctan(((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(a+b*x)^-1*(1+a+b*x))+-1//12*b^-4*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(2+17*a^2)*(1+a+b*x)+-1//12*b^-2*x^2*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(1+a+b*x)+1//3*a*b^-4*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(a+b*x)*(1+a+b*x))
@test integrate(x^2*arcsech(a+b*x), x) == :(1//3*x^3*arcsech(a+b*x)+-1//6*b^-3*(1+6*a^2)*arctan(((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(a+b*x)^-1*(1+a+b*x))+1//3*a^3*b^-3*arcsech(a+b*x)+-1//6*x*b^-2*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(1+a+b*x)+5//6*a*b^-3*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(1+a+b*x))
@test integrate(x*arcsech(a+b*x), x) == :((1/2)*x^2*arcsech(a+b*x)+a*b^-2*arctan(((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(a+b*x)^-1*(1+a+b*x))+-1//2*a^2*b^-2*arcsech(a+b*x)+-1//2*b^-2*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(1+a+b*x))
@test integrate(arcsech(a+b*x), x) == :(-2*b^-1*arctan(((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2))+b^-1*(a+b*x)*arcsech(a+b*x))
@test integrate(x^-1*arcsech(a+b*x), x) == :(-1//2*Polylogarithms.polylog(2,-1*exp(2*arcsech(a+b*x)))+arcsech(a+b*x)*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+arcsech(a+b*x)*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+-1*arcsech(a+b*x)*log(1+exp(2*arcsech(a+b*x)))+Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x))))
@test integrate(x^-2*arcsech(a+b*x), x) == :(-1*x^-1*arcsech(a+b*x)+-1*b*a^-1*arcsech(a+b*x)+2*b*a^-1*(1+-1*a^2)^-1//2*arctanh((1+a)^(1/2)*(1+-1a)^-1//2*tanh((1/2)*arcsech(a+b*x))))
@test integrate(x^-3*arcsech(a+b*x), x) == :(-1//2*x^-2*arcsech(a+b*x)+(1/2)*a^-2*b^2*arcsech(a+b*x)+-1*a^-2*b^2*(1+-1*a^2)^-3//2*(1+-2*a^2)*arctanh((1+a)^(1/2)*(1+-1a)^-1//2*tanh((1/2)*arcsech(a+b*x)))+(1/2)*b*a^-1*x^-1*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(1+-1*a^2)^-1*(1+a+b*x))
@test integrate(x^-4*arcsech(a+b*x), x) == :(-1//3*x^-3*arcsech(a+b*x)+-1//3*a^-3*b^3*arcsech(a+b*x)+1//3*a^-3*b^3*(1+-1*a^2)^-5//2*(2+-5*a^2+6*a^4)*arctanh((1+a)^(1/2)*(1+-1a)^-1//2*tanh((1/2)*arcsech(a+b*x)))+1//6*b*a^-1*x^-2*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(1+-1*a^2)^-1*(1+a+b*x)+-1//6*a^-2*b^2*x^-1*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(1+-1*a^2)^-2*(2+-5*a^2)*(1+a+b*x))
@test integrate(x^2*arcsech(a+b*x)^2, x) == :(-1//3*x*b^-2+1//3*x^3*arcsech(a+b*x)^2+2*a*b^-3*log(a+b*x)+-2//3*b^-3*arcsech(a+b*x)*arctan(exp(arcsech(a+b*x)))+-1//3*im*b^-3*Polylogarithms.polylog(2,im*exp(arcsech(a+b*x)))+1//3*im*b^-3*Polylogarithms.polylog(2,-1*im*exp(arcsech(a+b*x)))+1//3*a^3*b^-3*arcsech(a+b*x)^2+-4*a^2*b^-3*arcsech(a+b*x)*arctan(exp(arcsech(a+b*x)))+-2*im*a^2*b^-3*Polylogarithms.polylog(2,im*exp(arcsech(a+b*x)))+2*im*a^2*b^-3*Polylogarithms.polylog(2,-1*im*exp(arcsech(a+b*x)))+2*a*b^-3*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(1+a+b*x)*arcsech(a+b*x)+-1//3*b^-3*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(a+b*x)*(1+a+b*x)*arcsech(a+b*x))
@test integrate(x*arcsech(a+b*x)^2, x) == :((1/2)*x^2*arcsech(a+b*x)^2+-1*b^-2*log(a+b*x)+-1//2*a^2*b^-2*arcsech(a+b*x)^2+-1*b^-2*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(1+a+b*x)*arcsech(a+b*x)+-2*im*a*b^-2*Polylogarithms.polylog(2,-1*im*exp(arcsech(a+b*x)))+2*im*a*b^-2*Polylogarithms.polylog(2,im*exp(arcsech(a+b*x)))+4*a*b^-2*arcsech(a+b*x)*arctan(exp(arcsech(a+b*x))))
@test integrate(arcsech(a+b*x)^2, x) == :(b^-1*arcsech(a+b*x)^2*(a+b*x)+-4*b^-1*arcsech(a+b*x)*arctan(exp(arcsech(a+b*x)))+-2*im*b^-1*Polylogarithms.polylog(2,im*exp(arcsech(a+b*x)))+2*im*b^-1*Polylogarithms.polylog(2,-1*im*exp(arcsech(a+b*x))))
@test integrate(x^-1*arcsech(a+b*x)^2, x) == :((1/2)*Polylogarithms.polylog(3,-1*exp(2*arcsech(a+b*x)))+-2*Polylogarithms.polylog(3,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+-2*Polylogarithms.polylog(3,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+arcsech(a+b*x)^2*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+arcsech(a+b*x)^2*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+-1*arcsech(a+b*x)^2*log(1+exp(2*arcsech(a+b*x)))+-1*Polylogarithms.polylog(2,-1*exp(2*arcsech(a+b*x)))*arcsech(a+b*x)+2*Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))*arcsech(a+b*x)+2*Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))*arcsech(a+b*x))
@test integrate(x^-2*arcsech(a+b*x)^2, x) == :(-1*x^-1*arcsech(a+b*x)^2+-1*b*a^-1*arcsech(a+b*x)^2+-2*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+2*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+-2*b*a^-1*(1+-1*a^2)^-1//2*arcsech(a+b*x)*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+2*b*a^-1*(1+-1*a^2)^-1//2*arcsech(a+b*x)*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x))))
@test integrate(x^-3*arcsech(a+b*x)^2, x) == :(-1//2*x^-2*arcsech(a+b*x)^2+(1/2)*a^-2*b^2*arcsech(a+b*x)^2+a^-2*b^2*(1+-1*a^2)^-1*log(x*(a+b*x)^-1)+a^-2*b^2*(1+-1*a^2)^-3//2*Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+-1*a^-2*b^2*(1+-1*a^2)^-3//2*Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+-2*a^-2*b^2*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+2*a^-2*b^2*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+a^-2*b^2*(1+-1*a^2)^-3//2*arcsech(a+b*x)*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+-1*a^-2*b^2*(1+-1*a^2)^-3//2*arcsech(a+b*x)*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+-2*a^-2*b^2*(1+-1*a^2)^-1//2*arcsech(a+b*x)*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+2*a^-2*b^2*(1+-1*a^2)^-1//2*arcsech(a+b*x)*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+a^-1*b^2*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(1+-1*a^2)^-1*(1+-1*a*(a+b*x)^-1)^-1*(a+b*x)^-1*(1+a+b*x)*arcsech(a+b*x))
@test integrate(x*arcsech(a+b*x)^3, x) == :((1/2)*x^2*arcsech(a+b*x)^3+-3//2*b^-2*arcsech(a+b*x)^2+3//2*b^-2*Polylogarithms.polylog(2,-1*exp(2*arcsech(a+b*x)))+3*b^-2*arcsech(a+b*x)*log(1+exp(2*arcsech(a+b*x)))+-1//2*a^2*b^-2*arcsech(a+b*x)^3+-6*im*a*b^-2*Polylogarithms.polylog(3,im*exp(arcsech(a+b*x)))+6*im*a*b^-2*Polylogarithms.polylog(3,-1*im*exp(arcsech(a+b*x)))+6*a*b^-2*arcsech(a+b*x)^2*arctan(exp(arcsech(a+b*x)))+-3//2*b^-2*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*arcsech(a+b*x)^2*(1+a+b*x)+-6*im*a*b^-2*Polylogarithms.polylog(2,-1*im*exp(arcsech(a+b*x)))*arcsech(a+b*x)+6*im*a*b^-2*Polylogarithms.polylog(2,im*exp(arcsech(a+b*x)))*arcsech(a+b*x))
@test integrate(arcsech(a+b*x)^3, x) == :(b^-1*arcsech(a+b*x)^3*(a+b*x)+-6*im*b^-1*Polylogarithms.polylog(3,-1*im*exp(arcsech(a+b*x)))+-6*b^-1*arcsech(a+b*x)^2*arctan(exp(arcsech(a+b*x)))+6*im*b^-1*Polylogarithms.polylog(3,im*exp(arcsech(a+b*x)))+-6*im*b^-1*Polylogarithms.polylog(2,im*exp(arcsech(a+b*x)))*arcsech(a+b*x)+6*im*b^-1*Polylogarithms.polylog(2,-1*im*exp(arcsech(a+b*x)))*arcsech(a+b*x))
@test integrate(x^-1*arcsech(a+b*x)^3, x) == :(6*Polylogarithms.polylog(4,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+6*Polylogarithms.polylog(4,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+-3//4*Polylogarithms.polylog(4,-1*exp(2*arcsech(a+b*x)))+arcsech(a+b*x)^3*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+arcsech(a+b*x)^3*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+-1*arcsech(a+b*x)^3*log(1+exp(2*arcsech(a+b*x)))+-6*Polylogarithms.polylog(3,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))*arcsech(a+b*x)+-6*Polylogarithms.polylog(3,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))*arcsech(a+b*x)+3*arcsech(a+b*x)^2*Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+3*arcsech(a+b*x)^2*Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+-3//2*arcsech(a+b*x)^2*Polylogarithms.polylog(2,-1*exp(2*arcsech(a+b*x)))+3//2*Polylogarithms.polylog(3,-1*exp(2*arcsech(a+b*x)))*arcsech(a+b*x))
@test integrate(x^-2*arcsech(a+b*x)^3, x) == :(-1*x^-1*arcsech(a+b*x)^3+-1*b*a^-1*arcsech(a+b*x)^3+-6*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(3,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+6*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(3,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+-6*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))*arcsech(a+b*x)+-3*b*a^-1*(1+-1*a^2)^-1//2*arcsech(a+b*x)^2*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+3*b*a^-1*(1+-1*a^2)^-1//2*arcsech(a+b*x)^2*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+6*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))*arcsech(a+b*x))
@test integrate(x^-3*arcsech(a+b*x)^3, x) == :(-1//2*x^-2*arcsech(a+b*x)^3+(1/2)*a^-2*b^2*arcsech(a+b*x)^3+-6*a^-2*b^2*(1+-1*a^2)^-1//2*Polylogarithms.polylog(3,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+-3*a^-2*b^2*(1+-1*a^2)^-3//2*Polylogarithms.polylog(3,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+3*a^-2*b^2*(1+-1*a^2)^-1*Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+3*a^-2*b^2*(1+-1*a^2)^-1*Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+3*a^-2*b^2*(1+-1*a^2)^-3//2*Polylogarithms.polylog(3,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+6*a^-2*b^2*(1+-1*a^2)^-1//2*Polylogarithms.polylog(3,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+-3//2*a^-2*b^2*(1+-1*a^2)^-1*arcsech(a+b*x)^2+-6*a^-2*b^2*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))*arcsech(a+b*x)+-3*a^-2*b^2*(1+-1*a^2)^-3//2*Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))*arcsech(a+b*x)+-3*a^-2*b^2*(1+-1*a^2)^-1//2*arcsech(a+b*x)^2*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+3*a^-2*b^2*(1+-1*a^2)^-1*arcsech(a+b*x)*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+3*a^-2*b^2*(1+-1*a^2)^-1*arcsech(a+b*x)*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+3*a^-2*b^2*(1+-1*a^2)^-3//2*Polylogarithms.polylog(2,a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))*arcsech(a+b*x)+3*a^-2*b^2*(1+-1*a^2)^-1//2*arcsech(a+b*x)^2*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+6*a^-2*b^2*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))*arcsech(a+b*x)+-3//2*a^-2*b^2*(1+-1*a^2)^-3//2*arcsech(a+b*x)^2*log(1+-1*a*(1+(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+3//2*a^-2*b^2*(1+-1*a^2)^-3//2*arcsech(a+b*x)^2*log(1+-1*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(arcsech(a+b*x)))+3//2*a^-1*b^2*((1+a+b*x)^-1*(1+-1a+-1*b*x))^(1/2)*(1+-1*a^2)^-1*(1+-1*a*(a+b*x)^-1)^-1*(a+b*x)^-1*arcsech(a+b*x)^2*(1+a+b*x))
@test integrate(x^3*arcsech(x^(1/2)), x) == :(1//4*x^4*arcsech(x^(1/2))+-3//20*x^-1//2*(1+x^-1//2)^-1//2*(1+-1x)^3*(-1+x^-1//2)^-1//2+-1//4*x^-1//2*(1+x^-1//2)^-1//2*(-1+x^-1//2)^-1//2*(1+-1x)+1//4*x^-1//2*(1+x^-1//2)^-1//2*(1+-1x)^2*(-1+x^-1//2)^-1//2+1//28*x^-1//2*(1+x^-1//2)^-1//2*(1+-1x)^4*(-1+x^-1//2)^-1//2)
@test integrate(x^2*arcsech(x^(1/2)), x) == :(1//3*x^3*arcsech(x^(1/2))+-1//3*x^-1//2*(1+x^-1//2)^-1//2*(-1+x^-1//2)^-1//2*(1+-1x)+-1//15*x^-1//2*(1+x^-1//2)^-1//2*(1+-1x)^3*(-1+x^-1//2)^-1//2+2//9*x^-1//2*(1+x^-1//2)^-1//2*(1+-1x)^2*(-1+x^-1//2)^-1//2)
@test integrate(x*arcsech(x^(1/2)), x) == :((1/2)*x^2*arcsech(x^(1/2))+-1//2*x^-1//2*(1+x^-1//2)^-1//2*(-1+x^-1//2)^-1//2*(1+-1x)+1//6*x^-1//2*(1+x^-1//2)^-1//2*(1+-1x)^2*(-1+x^-1//2)^-1//2)
@test integrate(arcsech(x^(1/2)), x) == :(x*arcsech(x^(1/2))+-1*x^-1//2*(1+x^-1//2)^-1//2*(-1+x^-1//2)^-1//2*(1+-1x))
@test integrate(x^-1*arcsech(x^(1/2)), x) == :(arcsech(x^(1/2))^2+-1*Polylogarithms.polylog(2,-1*exp(2*arcsech(x^(1/2))))+-2*arcsech(x^(1/2))*log(1+exp(2*arcsech(x^(1/2)))))
@test integrate(x^-2*arcsech(x^(1/2)), x) == :(-1*x^-1*arcsech(x^(1/2))+(1/2)*x^-3//2*(1+x^-1//2)^-1//2*(-1+x^-1//2)^-1//2*(1+-1x)+(1/2)*x^-1//2*(1+x^-1//2)^-1//2*(1+-1x)^(1/2)*(-1+x^-1//2)^-1//2*arctanh((1+-1x)^(1/2)))
@test integrate(x^-3*arcsech(x^(1/2)), x) == :(-1//2*x^-2*arcsech(x^(1/2))+1//8*x^-5//2*(1+x^-1//2)^-1//2*(-1+x^-1//2)^-1//2*(1+-1x)+1//16*x^-3//2*(1+x^-1//2)^-1//2*(-1+x^-1//2)^-1//2*(3+-3x)+3//16*x^-1//2*(1+x^-1//2)^-1//2*(1+-1x)^(1/2)*(-1+x^-1//2)^-1//2*arctanh((1+-1x)^(1/2)))
@test integrate(x^-4*arcsech(x^(1/2)), x) == :(-1//3*x^-3*arcsech(x^(1/2))+1//18*x^-7//2*(1+x^-1//2)^-1//2*(-1+x^-1//2)^-1//2*(1+-1x)+1//48*x^-3//2*(1+x^-1//2)^-1//2*(-1+x^-1//2)^-1//2*(5+-5x)+1//72*x^-5//2*(1+x^-1//2)^-1//2*(-1+x^-1//2)^-1//2*(5+-5x)+5//48*x^-1//2*(1+x^-1//2)^-1//2*(1+-1x)^(1/2)*(-1+x^-1//2)^-1//2*arctanh((1+-1x)^(1/2)))
@test integrate(arcsech(x^-1), x) == :(x*arccosh(x)+-1*(1+x)^(1/2)*(-1+x)^(1/2))
@test integrate(x^-1*arcsech(a*x^n), x) == :((1/2)*n^-1*arcsech(a*x^n)^2+-1//2*n^-1*Polylogarithms.polylog(2,-1*exp(2*arcsech(a*x^n)))+-1*n^-1*arcsech(a*x^n)*log(1+exp(2*arcsech(a*x^n))))
@test integrate(x^-1*arcsech(a*x^5), x) == :(-1//10*Polylogarithms.polylog(2,-1*exp(2*arcsech(a*x^5)))+1//10*arcsech(a*x^5)^2+-1//5*arcsech(a*x^5)*log(1+exp(2*arcsech(a*x^5))))
@test integrate(arcsech(c*exp(a+b*x)), x) == :((1/2)*b^-1*arcsech(c*exp(a+b*x))^2+-1//2*b^-1*Polylogarithms.polylog(2,-1*exp(2*arcsech(c*exp(a+b*x))))+-1*b^-1*arcsech(c*exp(a+b*x))*log(1+exp(2*arcsech(c*exp(a+b*x)))))
@test integrate(x^3*exp(arcsech(a*x)), x) == :(1//4*x^4*exp(arcsech(a*x))+1//12*a^-1*x^3+-1//8*x*a^-3*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2)+1//8*a^-4*((1+a*x)^-1)^(1/2)*(1+a*x)^(1/2)*arcsin(a*x))
@test integrate(x*exp(arcsech(a*x)), x) == :((1/2)*x*a^-1+(1/2)*x^2*exp(arcsech(a*x))+(1/2)*a^-2*((1+a*x)^-1)^(1/2)*(1+a*x)^(1/2)*arcsin(a*x))
@test integrate(x^-4*exp(arcsech(a*x)), x) == :(-1//3*x^-3*exp(arcsech(a*x))+1//12*a^-1*x^-4+1//8*a*x^-2*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2)+1//8*a^3*((1+a*x)^-1)^(1/2)*(1+a*x)^(1/2)*arctanh((1+a*x)^(1/2)*(1+-1*a*x)^(1/2))+1//12*a^-1*x^-4*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2))
@test integrate(x^-5*exp(arcsech(a*x)), x) == :(-1//4*x^-4*exp(arcsech(a*x))+1//20*a^-1*x^-5+1//15*a*x^-3*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2)+1//20*a^-1*x^-5*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2)+2//15*a^3*x^-1*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2))
@test integrate(x^-6*exp(arcsech(a*x)), x) == :(-1//5*x^-5*exp(arcsech(a*x))+1//30*a^-1*x^-6+1//16*a^3*x^-2*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2)+1//16*a^5*((1+a*x)^-1)^(1/2)*(1+a*x)^(1/2)*arctanh((1+a*x)^(1/2)*(1+-1*a*x)^(1/2))+1//24*a*x^-4*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2)+1//30*a^-1*x^-6*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2))
@test integrate(x^-7*exp(arcsech(a*x)), x) == :(-1//6*x^-6*exp(arcsech(a*x))+1//42*a^-1*x^-7+1//35*a*x^-5*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2)+1//42*a^-1*x^-7*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2)+4//105*a^3*x^-3*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2)+8//105*a^5*x^-1*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2))
@test integrate(x^-8*exp(arcsech(a*x)), x) == :(-1//7*x^-7*exp(arcsech(a*x))+1//56*a^-1*x^-8+1//48*a*x^-6*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2)+1//56*a^-1*x^-8*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2)+5//128*a^5*x^-2*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2)+5//128*a^7*((1+a*x)^-1)^(1/2)*(1+a*x)^(1/2)*arctanh((1+a*x)^(1/2)*(1+-1*a*x)^(1/2))+5//192*a^3*x^-4*((1+a*x)^-1)^-1//2*(1+-1*a*x)^(1/2))
@test integrate(x^7*exp(arcsech(a*x^2)), x) == :(1//8*x^8*exp(arcsech(a*x^2))+1//24*a^-1*x^6+1//16*a^-4*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*arcsin(a*x^2)+-1//16*a^-3*x^2*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*(1+-1*a^2*x^4)^(1/2))
@test integrate(x^6*exp(arcsech(a*x^2)), x) == :(1//7*x^7*exp(arcsech(a*x^2))+2//35*a^-1*x^5+2//21*a^-7//2*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*Elliptic.F(arcsin(x*a^(1/2)),-1)+-2//21*x*a^-3*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*(1+-1*a^2*x^4)^(1/2))
@test integrate(x^4*exp(arcsech(a*x^2)), x) == :(1//5*x^5*exp(arcsech(a*x^2))+2//15*a^-1*x^3+-2//5*a^-5//2*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*Elliptic.F(arcsin(x*a^(1/2)),-1)+2//5*a^-5//2*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*Elliptic.E(arcsin(x*a^(1/2)),-1))
@test integrate(x^3*exp(arcsech(a*x^2)), x) == :(1//4*a^-1*x^2+1//4*x^4*exp(arcsech(a*x^2))+1//4*a^-2*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*arcsin(a*x^2))
@test integrate(x^2*exp(arcsech(a*x^2)), x) == :(1//3*x^3*exp(arcsech(a*x^2))+2//3*x*a^-1+2//3*a^-3//2*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*Elliptic.F(arcsin(x*a^(1/2)),-1))
@test integrate(x*exp(arcsech(a*x^2)), x) == :(a^-1*log(x)+(1/2)*x^2*exp(arcsech(a*x^2))+-1//2*a^-1*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*arctanh((1+-1*a^2*x^4)^(1/2)))
@test integrate(exp(arcsech(a*x^2)), x) == :(x*exp(arcsech(a*x^2))+-2*a^-1*x^-1+-2*a^-1//2*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*Elliptic.E(arcsin(x*a^(1/2)),-1)+2*a^-1//2*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*Elliptic.F(arcsin(x*a^(1/2)),-1)+-2*a^-1*x^-1*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*(1+-1*a^2*x^4)^(1/2))
@test integrate(x^-2*exp(arcsech(a*x^2)), x) == :(-1*x^-1*exp(arcsech(a*x^2))+2//3*a^-1*x^-3+-2//3*a^(1/2)*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*Elliptic.F(arcsin(x*a^(1/2)),-1)+2//3*a^-1*x^-3*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*(1+-1*a^2*x^4)^(1/2))
@test integrate(x^-3*exp(arcsech(a*x^2)), x) == :(-1//2*x^-2*exp(arcsech(a*x^2))+1//4*a^-1*x^-4+1//4*a*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*arctanh((1+-1*a^2*x^4)^(1/2))+1//4*a^-1*x^-4*((1+a*x^2)^-1)^(1/2)*(1+a*x^2)^(1/2)*(1+-1*a^2*x^4)^(1/2))
@test integrate(x^4*exp(2*arcsech(a*x)), x) == :(-1//2*a^-5*arctan(((1+a*x)^-1*(1+-1*a*x))^(1/2))+-1//30*a^-5*(1+a*x)^3*(4+45*((1+a*x)^-1*(1+-1*a*x))^(1/2))+1//4*a^-5*(1+a*x)*(4+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))+1//5*a^-5*(1+a*x)^4*(1+-1*a*x)+5//4*a^-5*((1+a*x)^-1*(1+-1*a*x))^(1/2)*(1+a*x)^2+1//10*a^-5*((1+a*x)^-1*(1+-1*a*x))^(1/2)*(1+a*x)^4*(5+-6*((1+a*x)^-1*(1+-1*a*x))^(1/2)))
@test integrate(x^3*exp(2*arcsech(a*x)), x) == :(-1*x*a^-3+1//4*a^-4*(1+a*x)^3*(1+-1*a*x)+1//6*a^-4*(1+a*x)^2*(3+-8*((1+a*x)^-1*(1+-1*a*x))^(1/2))+1//6*a^-4*((1+a*x)^-1*(1+-1*a*x))^(1/2)*(1+a*x)^3*(4+-3*((1+a*x)^-1*(1+-1*a*x))^(1/2)))
@test integrate(x^2*exp(2*arcsech(a*x)), x) == :(-2*a^-3*arctan(((1+a*x)^-1*(1+-1*a*x))^(1/2))+1//12*a^-3*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^4*(1+a*x)^3+(1/2)*a^-3*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))*(1+a*x)+-1//6*a^-3*((1+a*x)^-1*(1+-1*a*x))^(1/2)*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^3*(1+a*x)^2)
@test integrate(x*exp(2*arcsech(a*x)), x) == :(2*a^-2*log(1+a*x)+4*a^-2*log(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))+-1//2*a^-2*(1+a*x)^2+a^-2*(1+2*((1+a*x)^-1*(1+-1*a*x))^(1/2))*(1+a*x))
@test integrate(exp(2*arcsech(a*x)), x) == :(-1x+-4*a^-1*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+4*a^-1*arctan(((1+a*x)^-1*(1+-1*a*x))^(1/2)))
@test integrate(x^-1*exp(2*arcsech(a*x)), x) == :(-1*log(1+a*x)+-2*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2+-2*log(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))+2*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1)
@test integrate(x^-2*exp(2*arcsech(a*x)), x) == :(2*a*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2+-4//3*a*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-3)
@test integrate(x^-3*exp(2*arcsech(a*x)), x) == :(a^2*(2+-2*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+(1/2)*a^2*arctanh(((1+a*x)^-1*(1+-1*a*x))^(1/2))+-1*a^2*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-4+2*a^2*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-3+-3//2*a^2*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2)
@test integrate(x^-4*exp(2*arcsech(a*x)), x) == :(-1*a^3*(4+-4*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-1*a^3*(4+4*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+2*a^3*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-4+-7//3*a^3*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-3+-4//5*a^3*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-5+3//2*a^3*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2)
@test integrate(x^-5*exp(2*arcsech(a*x)), x) == :(a^4*(8+8*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-3*a^4*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-4+2*a^4*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-5+3*a^4*(8+-8*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-11//8*a^4*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2+-2//3*a^4*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-6+-1//8*a^4*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2+1//4*a^4*arctanh(((1+a*x)^-1*(1+-1*a*x))^(1/2))+8//3*a^4*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-3)
@test integrate(x^-6*exp(2*arcsech(a*x)), x) == :(-1*a^5*(4+-4*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-1*a^5*(4+4*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+2*a^5*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-6+4*a^5*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-4+-35//12*a^5*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-3+-18//5*a^5*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-5+-4//7*a^5*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-7+-1//12*a^5*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-3+1//8*a^5*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2+11//8*a^5*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2)
@test integrate(x^4*exp(-1*arcsech(a*x)), x) == :(-1*x*a^-4+-1//5*a^-5*((1+a*x)^-1*(1+-1*a*x))^(1/2)*(1+a*x)^5+-1//15*a^-5*(1+a*x)^3*(15+17*((1+a*x)^-1*(1+-1*a*x))^(1/2))+1//6*a^-5*(1+a*x)^2*(9+4*((1+a*x)^-1*(1+-1*a*x))^(1/2))+1//20*a^-5*(1+a*x)^4*(5+16*((1+a*x)^-1*(1+-1*a*x))^(1/2)))
@test integrate(x^3*exp(-1*arcsech(a*x)), x) == :(1//4*a^-4*arctan(((1+a*x)^-1*(1+-1*a*x))^(1/2))+-1//4*a^-4*((1+a*x)^-1*(1+-1*a*x))^(1/2)*(1+a*x)^4+-1//8*a^-4*(1+a*x)^2*(8+5*((1+a*x)^-1*(1+-1*a*x))^(1/2))+1//8*a^-4*(1+a*x)*(8+((1+a*x)^-1*(1+-1*a*x))^(1/2))+1//12*a^-4*(1+a*x)^3*(4+9*((1+a*x)^-1*(1+-1*a*x))^(1/2)))
@test integrate(x^2*exp(-1*arcsech(a*x)), x) == :(-1*x*a^-2+-1//3*a^-3*((1+a*x)^-1*(1+-1*a*x))^(1/2)*(1+a*x)^3+1//6*a^-3*(1+a*x)^2*(3+4*((1+a*x)^-1*(1+-1*a*x))^(1/2)))
@test integrate(x*exp(-1*arcsech(a*x)), x) == :(a^-2*arctan(((1+a*x)^-1*(1+-1*a*x))^(1/2))+(1/2)*a^-2*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))*(1+a*x)+1//4*a^-2*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^2*(1+a*x)^2)
@test integrate(exp(-1*arcsech(a*x)), x) == :(a^-1*log(1+a*x)+2*a^-1*log(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))+-1*a^-1*((1+a*x)^-1*(1+-1*a*x))^(1/2)*(1+a*x))
@test integrate(x^-1*exp(-1*arcsech(a*x)), x) == :(-2*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-2*arctan(((1+a*x)^-1*(1+-1*a*x))^(1/2)))
@test integrate(x^-2*exp(-1*arcsech(a*x)), x) == :(a*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-1*a*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2+-1*a*arctanh(((1+a*x)^-1*(1+-1*a*x))^(1/2)))
@test integrate(x^-3*exp(-1*arcsech(a*x)), x) == :(a^2*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2+-1*a^2*(2+-2*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-1*a^2*(2+2*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-2//3*a^2*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-3)
@test integrate(x^-4*exp(-1*arcsech(a*x)), x) == :(a^3*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-3+a^3*(2+2*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+a^3*(4+-4*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-1*a^3*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2+-1//2*a^3*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-4+-1//4*a^3*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2+-1//4*a^3*arctanh(((1+a*x)^-1*(1+-1*a*x))^(1/2)))
@test integrate(x^-5*exp(-1*arcsech(a*x)), x) == :(a^4*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-4+a^4*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2+-3*a^4*(8+-8*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-3*a^4*(8+8*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-4//3*a^4*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-3+-2//5*a^4*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-5+-1//6*a^4*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-3+1//4*a^4*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2)
@test integrate(x^-6*exp(-1*arcsech(a*x)), x) == :(a^5*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-5+a^5*(4+-4*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-1*a^5*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2+3*a^5*(8+8*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-13//8*a^5*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-4+-3//8*a^5*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2+-1//3*a^5*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-6+-1//8*a^5*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-4+-1//8*a^5*arctanh(((1+a*x)^-1*(1+-1*a*x))^(1/2))+1//4*a^5*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-3+19//12*a^5*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-3)
@test integrate(x^-7*exp(-1*arcsech(a*x)), x) == :(a^6*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-6+a^6*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2+-5*a^6*(16+-16*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-5*a^6*(16+16*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-1+-19//10*a^6*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-5+-11//6*a^6*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-3+-5//12*a^6*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-3+-2//7*a^6*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-7+-1//10*a^6*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-5+1//4*a^6*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-4+3//8*a^6*(1+-1*((1+a*x)^-1*(1+-1*a*x))^(1/2))^-2+9//4*a^6*(1+((1+a*x)^-1*(1+-1*a*x))^(1/2))^-4)
@test integrate(x^4*(1+-1*c^2*x^2)^-1*exp(arcsech(c*x)), x) == :(-1//2*c^-5*log(1+-1*c^2*x^2)+-1//2*c^-3*x^2+-2//3*c^-5*((1+c*x)^-1)^-1//2*(1+-1*c*x)^(1/2)+-1//3*c^-3*x^2*((1+c*x)^-1)^-1//2*(1+-1*c*x)^(1/2))
@test integrate(x^3*(1+-1*c^2*x^2)^-1*exp(arcsech(c*x)), x) == :(c^-4*arctanh(c*x)+-1*x*c^-3+(1/2)*c^-4*((1+c*x)^-1)^(1/2)*(1+c*x)^(1/2)*arcsin(c*x)+-1//2*x*c^-3*((1+c*x)^-1)^-1//2*(1+-1*c*x)^(1/2))
@test integrate(x^2*(1+-1*c^2*x^2)^-1*exp(arcsech(c*x)), x) == :(-1//2*c^-3*log(1+-1*c^2*x^2)+-1*c^-3*((1+c*x)^-1)^-1//2*(1+-1*c*x)^(1/2))
@test integrate(x*(1+-1*c^2*x^2)^-1*exp(arcsech(c*x)), x) == :(c^-2*arctanh(c*x)+c^-2*((1+c*x)^-1)^(1/2)*(1+c*x)^(1/2)*arcsin(c*x))
@test integrate((1+-1*c^2*x^2)^-1*exp(arcsech(c*x)), x) == :(c^-1*log(x)+-1//2*c^-1*log(1+-1*c^2*x^2)+-1*c^-1*((1+c*x)^-1)^(1/2)*(1+c*x)^(1/2)*arctanh((1+c*x)^(1/2)*(1+-1*c*x)^(1/2)))
@test integrate(x^-1*(1+-1*c^2*x^2)^-1*exp(arcsech(c*x)), x) == :(-1*c^-1*x^-1+-1*c^-1*x^-1*((1+c*x)^-1)^-1//2*(1+-1*c*x)^(1/2)+arctanh(c*x))
@test integrate(x^-2*(1+-1*c^2*x^2)^-1*exp(arcsech(c*x)), x) == :(c*log(x)+-1//2*c*log(1+-1*c^2*x^2)+-1//2*c^-1*x^-2+-1//2*c*((1+c*x)^-1)^(1/2)*(1+c*x)^(1/2)*arctanh((1+c*x)^(1/2)*(1+-1*c*x)^(1/2))+-1//2*c^-1*x^-2*((1+c*x)^-1)^-1//2*(1+-1*c*x)^(1/2))
@test integrate(x^-3*(1+-1*c^2*x^2)^-1*exp(arcsech(c*x)), x) == :(c^2*arctanh(c*x)+-1*c*x^-1+-1//3*c^-1*x^-3+-2//3*c*x^-1*((1+c*x)^-1)^-1//2*(1+-1*c*x)^(1/2)+-1//3*c^-1*x^-3*((1+c*x)^-1)^-1//2*(1+-1*c*x)^(1/2))
@test integrate((d*x+a*d*b^-1)^-1*arcsech(a+b*x), x) == :((1/2)*d^-1*arcsech(a+b*x)^2+-1//2*d^-1*Polylogarithms.polylog(2,-1*exp(2*arcsech(a+b*x)))+-1*d^-1*arcsech(a+b*x)*log(1+exp(2*arcsech(a+b*x))))
@test integrate(x^3*arcsech(a+b*x^4), x) == :(-1//2*b^-1*arctan(((1+a+b*x^4)^-1*(1+-1a+-1*b*x^4))^(1/2))+1//4*b^-1*(a+b*x^4)*arcsech(a+b*x^4))
@test integrate(x^(-1+n)*arcsech(a+b*x^n), x) == :(-2*b^-1*n^-1*arctan(((1+a+b*x^n)^-1*(1+-1a+-1*b*x^n))^(1/2))+b^-1*n^-1*(a+b*x^n)*arcsech(a+b*x^n))

max_stars_repo_path: test/7 Inverse hyperbolic functions/7.5 Inverse hyperbolic secant/7.5.2 Inverse hyperbolic secant functions.jl


text:
module IdealizedExperiments

using Oceananigans.Units

include("three_layer_constant_fluxes.jl")

two_day_suite_parameters = Dict{Symbol, Any}(
    :free_convection          => Dict{Symbol, Any}(:momentum_flux => 0.0,     :buoyancy_flux => 1.2e-7, :f => 1e-4),
    :strong_wind              => Dict{Symbol, Any}(:momentum_flux => -1e-3,   :buoyancy_flux => 0.0,    :f => 1e-4),
    :strong_wind_weak_cooling => Dict{Symbol, Any}(:momentum_flux => -7e-4,   :buoyancy_flux => 6e-8,   :f => 1e-4),
    :weak_wind_strong_cooling => Dict{Symbol, Any}(:momentum_flux => -3.3e-4, :buoyancy_flux => 1.1e-7, :f => 1e-4),
    :strong_wind_weak_heating => Dict{Symbol, Any}(:momentum_flux => -1e-3,   :buoyancy_flux => -4e-8,  :f => 1e-4),
    :strong_wind_no_rotation  => Dict{Symbol, Any}(:momentum_flux => -2e-4,   :buoyancy_flux => 0.0,    :f => 0.0),
)

for (name, set) in two_day_suite_parameters
    set[:name] = string(name)
    set[:stop_time] = 2days
end

four_day_suite_parameters = Dict{Symbol, Any}(
    :free_convection          => Dict{Symbol, Any}(:momentum_flux => 0.0,     :buoyancy_flux => 7.0e-8, :f => 1e-4),
    :strong_wind              => Dict{Symbol, Any}(:momentum_flux => -8e-4,   :buoyancy_flux => 0.0,    :f => 1e-4),
    :strong_wind_weak_cooling => Dict{Symbol, Any}(:momentum_flux => -6.5e-4, :buoyancy_flux => 4e-8,   :f => 1e-4),
    :weak_wind_strong_cooling => Dict{Symbol, Any}(:momentum_flux => -3e-4,   :buoyancy_flux => 7e-8,   :f => 1e-4),
    :strong_wind_no_rotation  => Dict{Symbol, Any}(:momentum_flux => -1e-4,   :buoyancy_flux => 0.0,    :f => 0.0),
)

for (name, set) in four_day_suite_parameters
    set[:name] = string(name)
    set[:stop_time] = 4days
end

six_day_suite_parameters = Dict{Symbol, Any}(
    :free_convection          => Dict{Symbol, Any}(:momentum_flux => 0.0,     :buoyancy_flux => 5e-8, :f => 1e-4),
    :strong_wind              => Dict{Symbol, Any}(:momentum_flux => -7e-4,   :buoyancy_flux => 0.0,  :f => 1e-4),
    :strong_wind_weak_cooling => Dict{Symbol, Any}(:momentum_flux => -5.5e-4, :buoyancy_flux => 3e-8, :f => 1e-4),
    :weak_wind_strong_cooling => Dict{Symbol, Any}(:momentum_flux => -2.2e-4, :buoyancy_flux => 5e-8, :f => 1e-4),
    :strong_wind_no_rotation  => Dict{Symbol, Any}(:momentum_flux => -7e-5,   :buoyancy_flux => 0.0,  :f => 0.0),
)

for (name, set) in six_day_suite_parameters
    set[:name] = string(name)
    set[:stop_time] = 6days
end

end # module

max_stars_repo_name: ali-ramadhan/LESbrary

max_issues_repo_name: ali-ramadhan/LESbrary

max_forks_repo_name: ali-ramadhan/LESbrary


text:
export SCS_init, SCS_solve, SCS_finish, SCS_version

function SCS_version()
    return unsafe_string(ccall((:scs_version, SCS.scs), Cstring, ()))
end

for (T, lib) in zip([SCS.Direct, SCS.Indirect], [SCS.scs, SCS.scsindir])
    @eval begin
        function SCS_init(::Type{$T}, data::SCSData, cone::SCSCone)
            # Initialize the info struct
            info = SCSInfo(0, convert(Int128, 0), convert(Int128, 0), 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

            p_work = ccall((:scs_init, $lib), Ptr{Void},
                (Ptr{SCSData}, Ptr{SCSCone}, Ptr{SCSInfo}),
                &data, &cone, &info)

            return p_work, info
        end

        function SCS_solve(::Type{$T}, p_work::Ptr{Void}, data::SCSData, cone::SCSCone, info::SCSInfo, solution::SCSSolution)
            solution_ptr = pointer([solution])

            info_ptr = pointer([info])
            status = ccall((:scs_solve, $lib), Int,
                (Ptr{Void}, Ptr{SCSData}, Ptr{SCSCone}, Ptr{SCSSolution}, Ptr{SCSInfo}),
                p_work, &data, &cone, solution_ptr, info_ptr)

            solution = unsafe_load(solution_ptr)
            info = unsafe_load(info_ptr)

            return status, solution, info, p_work
        end

        function SCS_solve(::Type{$T}, p_work::Ptr{Void}, data::SCSData, cone::SCSCone, info::SCSInfo)
            solution = SCSSolution(pointer(zeros(data.n)), pointer(zeros(data.m)), pointer(zeros(data.m)))
            return SCS_solve($T, p_work, data, cone, info, solution)
        end

        function SCS_solve(::Type{$T}, data::SCSData, cone::SCSCone, solution::SCSSolution)
            p_work, info = SCS_init($T, data, cone)
            return SCS_solve($T, p_work, data, cone, info, solution)
        end

        function SCS_finish(::Type{$T}, p_work::Ptr{Void})
            ccall((:scs_finish, $lib), Void,
                (Ptr{Void}, ),
                p_work)
        end
    end
end

SCS_init(data::SCSData, cone::SCSCone) = SCS_init(SCS.Direct, data, cone)

function SCS_solve(p_work::Ptr{Void}, data::SCSData, cone::SCSCone,
        info::SCSInfo, solution::SCSSolution)
    return SCS_solve(SCS.Direct, p_work, data, cone, info, solution)
end

function SCS_solve(p_work::Ptr{Void}, data::SCSData, cone::SCSCone,
        info::SCSInfo)
    return SCS_solve(SCS.Direct, p_work, data, cone, info)
end

function SCS_solve(data::SCSData, cone::SCSCone, solution::SCSSolution)
    return SCS_solve(SCS.Direct, data, cone, solution)
end

SCS_finish(p_work::Ptr{Void}) = SCS_finish(SCS.Direct, p_work)


text:
function solve()
    n = parse(Int, readline())
    a = [parse(Int, x) for x in split(readline())]

    z = [(a[i], i) for i=1:n]
    sort!(z, rev=true)

    dp = zeros(Int, (n+1, n+1))
    for s=1:n, i=0:s
        li = i
        ri = s-i
        l = if li > 0
            dp[li, ri+1] + z[s][1] * abs(li - z[s][2])
        else
            0
        end
        r = if ri > 0
            dp[li+1, ri] + z[s][1] * abs((n-ri+1) - z[s][2])
        else
            0
        end
        dp[li+1, ri+1] = max(l, r)
    end
    maximum([dp[i+1, n-i+1] for i=0:n])
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
module OpticsPolynomials

include("jacobi.jl")
include("cheby.jl")
include("legendre.jl")
include("zernike.jl")

end


text:
# @everywhere using RegularizedLeastSquares, LinearOperators, TensorDecompositions, LinearAlgebra
#@everywhere using RegularizedLeastSquares, LinearOperators, LinearAlgebra
export smRecovery

###################################################################
# System Matrix recovery using DCT-sparsity and low-rank constraing
###################################################################
function smRecovery(y::Matrix{T}, samplingIdx::Array{Int64}, params::Dict) where T
  # constrain number of threads for distributed computation
  BLAS.set_num_threads(1)

  shape = params[:shape]

  # sampling operator
  P = SamplingOp(samplingIdx, shape)

  # sparsifying transformation and regularization
  if !haskey(params, :sparseTrafo)
    params[:sparseTrafo] = DCTOp(ComplexF64,params[:shape],2)
  end
  reg1Name = get(params, :reg1, "L1")
  reg1 = Regularization(reg1Name, params[:λ1]; params...)

  # low rank regularization
  lrProx = get(params, :reg2, "Nothing")
  λ_lr = get(params, :λ_lr, 1.0)
  reg2 = lrReg(lrProx, λ_lr, params)

  # precalculate inverse for split Bregman iteration
  precon = invPrecon(samplingIdx, prod(shape); params...)
  params[:precon] = precon
  params[:iterationsCG] = 1

  # normalized measurement
  y2, y_norm = getNormalizedMeasurement(y)

  # # reconstruction
  # sf = map(x->splitBregman(P,x,[reg1,reg2]; params...),y2)
  #
  # sfMat = zeros(ComplexF64,prod(shape),size(y,2))
  # for i=1:size(y,2)
  #   # undo normalization
  #   sfMat[:,i] .= y_norm[i]*sf[i]
  # end

  # reconstruction
  sfMat = zeros(ComplexF64,prod(shape),size(y,2))
  solver = SplitBregman(P;reg=[reg1,reg2], params...)
  for k=1:size(y,2)
    sfMat[:,k] = solve(solver, y2[k])
    # undo normalization
    sfMat[:,k] *= y_norm[k]
  end

  return sfMat
end

#########################
# low rank regularization
#########################
function lrReg(lrProx::AbstractString, λ::AbstractFloat, params::Dict{Symbol,Any})
  if lrProx == "LR"
    proxMap = (x,y;kargs...)->proxLR!(x,y,params[:shape];kargs...)
    reg = Regularization(prox! = proxMap, λ=λ, params=params)
  elseif lrProx == "FR"
    proxMap = (x,y;kargs...)->proxFR!(x,y,params[:shape];kargs...)
    reg = Regularization(prox! = proxMap, λ=λ, params=params)
  elseif lrProx == "Nothing"
    reg = Regularization(prox! = (x,y;kargs...)->x, λ=λ, params=params)
  else
    error("proximal map $(lrProx) is not supported")
  end
  return reg
end

# low rank regularization for 3d (using HOSVD)
function proxLR!(x, λ::Real, shape::NTuple{3,Int64}; kargs...)
  guvw_r = hosvd(reshape(real.(x), shape),shape)
  guvw_i = hosvd(reshape(imag.(x), shape),shape)
  proxL1!(guvw_r.core,λ)
  proxL1!(guvw_i.core,λ)
  x[:] .= vec(compose(guvw_r)+1im*compose(guvw_i))
end

# low rank regularization for 2d (using SVD)
function proxLR!(x, λ::Real, shape::NTuple{2,Int64}; kargs...)
  U,S,V = svd(reshape(x,shape))
  proxL1!(S,λ)
  x[:] .= vec(U*Matrix(Diagonal(S))*V')
end

# fixed rank constraing for 3d (using HOSVD)
function proxFR!(x, λ::Real, shape::NTuple{3,Int64}; rank::NTuple{3,Int64}=(1,1,1), kargs...)
  guvw_r = hosvd(reshape(real.(x), shape),rank)
  guvw_i = hosvd(reshape(imag.(x), shape),rank)
  x[:] = vec(compose(guvw_r)+1im*compose(guvw_i))
end

# fixed rank constraing for 3d (using SVD)
function proxFR!(x, λ::Real, shape::NTuple{2,Int64}; rank::Int64=1, kargs...)
  U,S,V = svd(reshape(x, shape))
  S[rank+1:end] .= 0.0
  x[:] .= vec(U*Matrix(Diagonal(S))*V')
end

##########################
# inverse of normal matrix
##########################
mutable struct diagPrecon{T}
  diag::Vector{T}
end

function ldiv!(y::Vector{T},P::diagPrecon{T},x::Vector{T}) where T
  y[:] .= P.diag .* x
end

function ldiv!(P::diagPrecon{T},x::Vector{T}) where T
  x[:] .= P.diag .* x
end

\(P::diagPrecon{T}, x::Vector{T}) where T = P.diag .* x

function invPrecon(samplingIdx, ncol; μ::Float64 = 1.e-2, λ::Float64 = 1.e2, ρ::Float64=1.0, kargs...)
  normalP = normalPDiag(samplingIdx, ncol)
  diag = [μ*normalP[i]+λ+ρ for i=1:length(normalP)]
  diag = diag.^-1
  return diagPrecon(diag)
end

function normalPDiag(samplingIdx, ncol; T=ComplexF64)
  normalP = zeros(T,ncol)
  for i=1:length(samplingIdx)
    normalP[samplingIdx[i]] = 1.0
  end
  return normalP
end

function getNormalizedMeasurement(y::Matrix{T}) where T
  y_norm = [norm(y[:,k]) for k=1:size(y,2)]
  y2 = @DArray [ComplexF64.(y[:,k]/y_norm[k]) for k=1:size(y,2)]
  return y2, y_norm
end


text:
# s_t = TT*s_{t-1} + RR*ε_t, ε_t ~ iidN(0, QQ)
# y_t = DD+ ZZ*s_t + η_t, η_t ~ iidN(0, HH)
#y is nobs x ny with each row being time period

function chand_recursion(y::Matrix{S},
                         T::Matrix{S}, R::Matrix{S}, C::Vector{S},
                         Q::Matrix{S}, Z::Matrix{S}, D::Vector{S}, H::Matrix{S},
                         s_pred::Vector{S} = Vector{S}(0), P_pred::Matrix{S} = Matrix{S}(0,0);
                         allout::Bool = false, Nt0::Int = 0) where {S<:AbstractFloat}
    # Dimensions
    Ns = size(T, 1) # number of states
    Ny, Nt = size(y) # number of periods of data

    # Initialize s_pred and P_pred
    if isempty(s_pred) || isempty(P_pred)
        s_pred, P_pred = init_stationary_states(T, R, C, Q)
    end

    # V_{t|t-1} = Var(y_t|y_{1:t-1})
    V_pred = Z*P_pred*Z' + H
    V_pred = 0.5*(V_pred+V_pred')
    invV_pred = inv(V_pred)

    invV_t1 = invV_pred

    # We write ΔP in terms of W_t and M_t           ΔP = W_t * M_t * W_t'
    W_t = T*P_pred*Z' # (Eq 12)
    M_t = -invV_pred    # (Eq 13)
    kal_gain = W_t*invV_pred

    # Initialize loglikelihoods to zeros so that the ones we don't update (those in presample) contribute zero towards sum of loglikelihoods.
    loglh = Vector{Float64}(Nt)
    zero_vec = zeros(Ny)
    ν_t = zero_vec
    P =  P_pred

    for t in 1:Nt
        # Step 1: Compute forecast error, ν_t and evaluate likelihoood
        yhat = Z*s_pred + D
        ν_t = y[:, t] - yhat

        # This is a symptom of a bad parameter draw: tells callee to simply redraw
        if !isposdef(V_pred)
              return NaN
        else
            loglh[t] = logpdf(MvNormal(zero_vec, V_pred), ν_t)
        end

        # Step 2: Compute s_{t+1} using Eq. 5
        if t < Nt
            s_pred = T*s_pred + kal_gain*ν_t
            if allout
                P = P + W_t*M_t*W_t'
            end
        end

        # Intermediate calculations to re-use
        ZW_t = Z*W_t
        MWpZp = M_t*(ZW_t')
        TW_t = T*W_t

        # Step 3: Update forecast error variance F_{t+1} (Eq 19)
        V_t1 = V_pred
        invV_t1 = invV_pred
        V_pred  = V_pred + ZW_t*MWpZp        # F_{t+1}
        V_pred  = 0.5*(V_pred+V_pred')
        invV_pred = inv(V_pred)

        # Step 4: Update Kalman Gain (Eq 20). Recall that kalgain = K_t * V_t⁻¹
        # Kalgain_{t+1} = (Kalgain_t*V_{t-1} + T*W_t*M_t*W_t'*Z')V_t⁻¹
        kal_gain = (kal_gain*V_t1 + TW_t*MWpZp)*invV_pred

        # Step 5: Update W
        W_t = TW_t - kal_gain*ZW_t    # W_{t+1}

        # Step 6: Update M
        M_t = M_t + MWpZp*invV_t1*MWpZp'     # M_{t+1}
        M_t = 0.5*(M_t + M_t')
    end
    loglh = remove_presample!(Nt0, loglh)
    if allout
        s_TT = s_pred + P'*Z'*invV_t1*(ν_t)
        P_TT = P - P'*Z'*invV_t1*Z*P
        return sum(loglh), loglh, s_TT, P_TT
    else
        return sum(loglh)
    end
end

function remove_presample!(Nt0::Int, loglh::Vector{S}) where {S<:AbstractFloat}
    out = remove_presample!(Nt0, loglh, Array{Float64, 2}(0,0), Array{Float64,3}(0,0,0), Array{Float64, 2}(0,0), Array{Float64, 3}(0,0,0), outputs = [:loglh])
    return out[1]
end

#=function chand_recursion(regime_indices::Vector{Range{Int}}, y::Matrix{S},
                         Ts::Vector{Matrix{S}}, Rs::Vector{Matrix{S}}, Cs::Vector{Vector{S}},
                         Qs::Vector{Matrix{S}}, Zs::Vector{Matrix{S}}, Ds::Vector{Vector{S}},
                         Es::Vector{Matrix{S}},
                         s_0::Vector{S} = Vector{S}(0), P_0::Matrix{S} = Matrix{S}(0,0);
                         outputs::Vector{Symbol} = [:loglh, :pred, :filt],
                         Nt0::Int = 0) where {S<:AbstractFloat}
    Ns = size(Ts[1], 1)
    Nt = size(y, 2)
    @assert first(regime_indices[1]) == 1
    @assert last(regime_indices[end]) == Nt

    loglh = Vector{S}(0)

    for i=1:length(regime_indices)
        ts = regime_indices[i]
        loglh, loglh_i = chand_recursion(y[:ts], Ts[i], Rs[i], Cs[i], Qs[i], Zs[i], Ds[i], Es[i],
                                         s_t, P_t; allout = true, Nt0 = 0)
        loglh[ts] = loglh_i
    end
    loglh, _, _, _, _ = remove_presample!(Nt0, loglh, zeros(0), zeros(0), zeros(0), zeros(0))
    return sum(loglh), loglh
end =#


text:
mutable struct PushApplier{T}
    path::TreePath{RigidBody{T}, Joint{T}}
    application_point::Point3D{SVector{3, T}}
    jacobian::PointJacobian{Matrix{T}}
    new_push::Bool
    t0::T
    Δt::T
    force::FreeVector3D{SVector{3, T}}
end

function PushApplier(mechanism::Mechanism{T}, application_point::Point3D;
        Δt::Number=zero(T), force::AbstractVector=zero(SVector{3, T})) where T
    world = root_body(mechanism)
    worldframe = root_frame(mechanism)
    frame = application_point.frame
    treepath = RigidBodyDynamics.path(mechanism, world, body_fixed_frame_to_body(mechanism, frame))
    jacobian = PointJacobian(worldframe, zeros(3, num_velocities(mechanism)))
    PushApplier(treepath, application_point, jacobian, false, typemax(T), T(Δt), FreeVector3D(worldframe, SVector{3, T}(force)))
end

function (controller::PushApplier)(τ, t, state::MechanismState)
    if controller.new_push
        controller.t0 = t
        controller.new_push = false
    end
    if controller.t0 <= t <= controller.t0 + controller.Δt
        jacobian = controller.jacobian
        point = transform(state, controller.application_point, jacobian.frame)
        point_jacobian!(jacobian, state, controller.path, point)
        mul!(τ, transpose(jacobian), controller.force)
    else
        τ .= 0
    end
    τ
end