text:
using Dates

"""
    read_pids(sub_base::String, offset::Integer, N_ids::Integer)

Reads the `PID` block in the subfind output.
"""
function read_pids(sub_base::String, N_ids::Integer, offset::Integer)

    # choose correct file
    sub_file = select_file(sub_base, 0)

    # read the header
    sub_header = read_subfind_header(sub_file)

    # read info for PID datatype
    info = read_info(sub_file)
    pid_info = info[getfield.(info, :block_name) .== "PID"][1]

    # if there are fewer IDs in the file than in the halo they are distributed over multiple files
    if (sub_header.nfof - offset) < N_ids

        # allocate array for IDs
        ids = Array{pid_info.data_type}(undef, N_ids)

        ids_read = 0
        files_read = 0

        # the remainig offset
        offset_remaining = offset

        # loop over sub files until all IDs are read
        while ids_read < N_ids

            sub_file = select_file(sub_base, files_read)

            if !isfile(sub_file)
                error("$sub_file does not exist!")
            end

            # read the header
            sub_header = read_subfind_header(sub_file)

            # if the offset is larger than the number of IDs in this file
            if offset_remaining >= sub_header.nfof
                # read next file
                files_read += 1
                # subtract IDs in this file from offset
                offset_remaining -= sub_header.nfof
                continue
            else
                # find starting position for read-in
                start_pos = offset_remaining + 1

                # number of particles remaining to be read
                if ( sub_header.nfof - offset_remaining < N_ids - ids_read )
                    n_to_read = sub_header.nfof - offset_remaining
                else
                    n_to_read = N_ids - ids_read
                end

                ids[ids_read+1:ids_read+n_to_read] = read_subfind(sub_file, "PID")[start_pos:start_pos+n_to_read-1]

                # update number of IDs alrady read
                ids_read += n_to_read

                # increase for next file read
                files_read += 1
            end

        end # while

        return ids
    else # they can be read from one file
        # position of first ID in array
        start_pos = offset + 1
        return read_subfind(sub_file, "PID")[start_pos:start_pos+N_ids-1]
    end
end

"""
    read_ids_in_halo( sub_base::String, halo::HaloID;
                      halo_type::Integer=1, verbose::Bool=true)

Reads the IDs of all particles contained in a `halo`.
"""
function read_ids_in_halo( sub_base::String, halo::HaloID;
                           halo_type::Integer=1, verbose::Bool=true)

    # select subfind file to read
    sub_file = select_file(sub_base, halo.file)

    # select block with number of particles in fof for
    if halo_type == 1
        # halos
        len_block = "GLEN"
        off_block = "GOFF"

    elseif halo_type == 2
        # subhalos
        len_block = "SLEN"
        off_block = "SOFF"
    end

    if verbose
        @info "Reading number of particles in halo..."
        t1 = Dates.now()
    end

    # read number of IDs in the halo
    N_ids = Int64(read_subfind(sub_file, len_block)[halo.id])

    if verbose
        t2 = Dates.now()
        @info "N_particles to read: $N_ids. Took: $(t2 - t1)"
    end
    # read offset in PID array
    offset = Int64(read_subfind(sub_file, off_block)[halo.id])

    if verbose
        @info "Reading IDs in halo..."
        t1 = Dates.now()
    end
    # read all IDs of the particles contained in a halo
    halo_ids = read_pids(sub_base, N_ids, offset)

    if verbose
        t2 = Dates.now()
        @info "IDs read. Took: $(t2 - t1)"
    end

    return halo_ids
end

"""
    function get_pos_block_name(sub_file, halo_type)

Returns the position block name depending on the halo type.
"""
function get_pos_block_name(sub_file, halo_type)
    if halo_type == 1
        # halos
        return "GPOS"
    elseif halo_type == 2
        # subhalos
        return "SPOS"
    else
        error("No position block for halo type $halo_type")
    end
end

"""
    function get_rad_block_name(sub_file, halo_type)

Returns the radius block name depending on the halo type.
"""
function get_rad_block_name(sub_file, halo_type)
    if halo_type == 1
        # halos
        rad_block = "R200"
        if !block_present(sub_file, rad_block)
            rad_block = "RMEA"

            if !block_present(sub_file, rad_block)
                error("Neither R200 nor RMEA present!")
            end

        end
        return rad_block
    elseif halo_type == 2
        # subhalos
        return "RHMS"
    else
        error("No radius block for halo type $halo_type")
    end
end

"""
    read_particles_in_halo(snap_base::String, blocks::Array{String},
                                sub_base::String, halo::HaloID;
                                radius::Union{Real,Nothing}=nothing,
                                rad_scale::Real=1.0, halo_type::Integer=1,
                                parttype::Integer=0, verbose::Bool=true,
                                use_keys::Bool=true)

Reads all particles of type `parttype` that are contained in a halo defined by its `HaloID`.
If `radius` is given (in simulation units), particles are read within at least this radius
times `rad_scale`. Otherwise, R200, RMEA, or RHMS times `rad_scale` is used depending on
`halo_type` (1, 1, and 2, respectively).
Returns a `Dict` with each of the `blocks` as entries.
"""
function read_particles_in_halo(snap_base::String, blocks::Array{String},
                                sub_base::String, halo::HaloID;
                                radius::Union{Real,Nothing}=nothing,
                                rad_scale::Real=1.0, halo_type::Integer=1,
                                parttype::Integer=0, verbose::Bool=true,
                                use_keys::Bool=true)

    # select subfind file to read
    sub_file = select_file(sub_base, halo.file)
    h = read_header(sub_file)

    if verbose
        @info "Reading IDs in halo..."
        t1 = Dates.now()
    end

    # read all IDs of the particles contained in a halo
    halo_ids = read_ids_in_halo(sub_base, halo, halo_type=halo_type, verbose=verbose)

    if verbose
        t2 = Dates.now()
        @info "IDs read. Took: $(t2 - t1)"

        @info "Reading Data..."
        t1 = Dates.now()
    end

    # position of halo
    pos_block = get_pos_block_name(sub_file, halo_type)
    halo_pos = read_subfind(sub_file, pos_block)[:,halo.id]

    # initial search radius for read-in
    if isnothing(radius)
        rad_block = get_rad_block_name(sub_file, halo_type)
        initial_radius = rad_scale * read_subfind(sub_file, rad_block)[halo.id]
    else
        initial_radius = rad_scale * radius
    end

    # read ids in the halo
    data = read_particles_by_id(snap_base, halo_ids, blocks,
                              parttype=parttype, verbose=verbose,
                              pos0=halo_pos,
                              r0=initial_radius,
                              use_keys=use_keys)

    if verbose
        t2 = Dates.now()
        @info "Data read. Took: $(t2 - t1)"
        println()
        @info "Done!"
    end

    return data

end

"""
    read_particles_in_halo( snap_base::String, block::String,
                            sub_base::String, halo::HaloID;
                            kwargs...)

Reads all particles of type `parttype` that are contained in a halo defined by its `HaloID`.
If `radius` is given (in simulation units), particles are read within at least this radius
times `rad_scale`. Otherwise, R200, RMEA, or RHMS times `rad_scale` is used depending on
`halo_type` (1, 1, and 2, respectively).
Returns an `Array` with the requested `block`.
"""
function read_particles_in_halo(snap_base::String, block::String,
                                sub_base::String, halo::HaloID;
                                kwargs...)

    data = read_particles_in_halo(snap_base, [block], sub_base, halo; kwargs...)

    return data[block]
end


text:
# Julia wrapper for header: /usr/local/include/libavutil/audioconvert.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_get_channel_layout,
    av_get_channel_layout_string,
    av_get_channel_layout_nb_channels

function av_get_channel_layout(name)
    ccall((:av_get_channel_layout,libavutil),UInt64,(Ptr{UInt8},),name)
end

function av_get_channel_layout_string(buf,buf_size::Integer,nb_channels::Integer,channel_layout::UInt64)
    ccall((:av_get_channel_layout_string,libavutil),Void,(Ptr{UInt8},Cint,Cint,UInt64),buf,buf_size,nb_channels,channel_layout)
end

function av_get_channel_layout_nb_channels(channel_layout::UInt64)
    ccall((:av_get_channel_layout_nb_channels,libavutil),Cint,(UInt64,),channel_layout)
end

ext: jl

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
using LinearAlgebra
import FINUFFT

#=
General comment: FINUFFT requires dedicated plans for forward and the adjoint NFFT.
However, it is right now not possible to cache two plans in parallel since the plan
seem to have global state. We therefore only cache one plan and need to create a new
one "on demand" if the other one is needed. This is the purpose of the `isAdjoint` flag.
=#

mutable struct FINUFFTPlan{T,D} <: AbstractNFFTPlan{T,D,1}
  N::NTuple{D,Int64}
  M::Int64
  x::Matrix{T}
  m::Int
  σ::T
  reltol::T
  fftflags::UInt32
  planTrafo::FINUFFT.finufft_plan{T}
  planAdjoint::FINUFFT.finufft_plan{T}
end

################
# constructors
################

function FINUFFTPlan(x::Matrix{T}, N::NTuple{D,Int};
              dims::Union{Integer,UnitRange{Int64}}=1:D,
              fftflags=UInt32(NFFT.FFTW.ESTIMATE),
              kargs...) where {D,T}

  if dims != 1:D
  error("FINUFFT directional plans not yet implemented!")
  end

  M = size(x,2)

  m, σ, reltol = accuracyParams(; kargs...)

  reltol = max(reltol, 1.0e-15)

  x_ = T.(2π) * x
  nodes = ntuple(d->vec(x_[d,:]), D)

  planTrafo = FINUFFT.finufft_makeplan(2, collect(N), -1, 1, reltol;
                          nthreads = Threads.nthreads(), dtype=T,
                          fftw = fftflags, upsampfac=2.0, debug=0)

  FINUFFT.finufft_setpts!(planTrafo, nodes...)

  planAdjoint = FINUFFT.finufft_makeplan(1, collect(N), 1, 1, reltol;
                          nthreads = Threads.nthreads(), dtype=T,
                          fftw = fftflags, upsampfac=2.0, debug=0)

  FINUFFT.finufft_setpts!(planAdjoint, nodes...)

  p = FINUFFTPlan(N, M, x_, m, T(σ), T(reltol), fftflags, planTrafo, planAdjoint)

  finalizer(p -> begin
    #println("Run FINUFFT finalizer")
    FINUFFT.finufft_destroy!(p.planTrafo)
    FINUFFT.finufft_destroy!(p.planAdjoint)
  end, p)

  return p
end

function Base.show(io::IO, p::FINUFFTPlan)
  print(io, "FINUFFTPlan")
end

AbstractNFFTs.size_in(p::FINUFFTPlan) = Int.(p.N)
AbstractNFFTs.size_out(p::FINUFFTPlan) = (Int(p.M),)

function AbstractNFFTs.plan_nfft(::Type{<:Array}, x::Matrix{T}, N::NTuple{D,Int}, rest...;
                   timing::Union{Nothing,TimingStats} = nothing, kargs...) where {T,D}
  t = @elapsed begin
    p = FINUFFTPlan(x, N, rest...; kargs...)
  end
  if timing != nothing
    timing.pre = t
  end
  return p
end

function LinearAlgebra.mul!(fHat::StridedArray, p::FINUFFTPlan{T,D}, f::AbstractArray;
             verbose=false, timing::Union{Nothing,TimingStats} = nothing) where {T,D}

  FINUFFT.finufft_exec!(p.planTrafo, f, fHat)

  return fHat
end

function LinearAlgebra.mul!(f::StridedArray, pl::Adjoint{Complex{T},<:FINUFFTPlan{T,D}}, fHat::AbstractArray;
                     verbose=false, timing::Union{Nothing,TimingStats} = nothing) where {T,D}
  p = pl.parent

  FINUFFT.finufft_exec!(p.planAdjoint, fHat, f)

  return f
end

mutable struct FINNUFFTPlan{T} <: AbstractNNFFTPlan{T,1,1}
  N::Int64
  M::Int64
  x::Matrix{T}
  y::Matrix{T}
  m::Int
  σ::T
  reltol::T
  fftflags::UInt32
end

function FINNUFFTPlan(x::Matrix{T}, y::Matrix{T};
  fftflags=UInt32(NFFT.FFTW.ESTIMATE),
  kargs...) where {T}

  N = size(y,2)
  M = size(x,2)

  m, σ, reltol = accuracyParams(; kargs...)

  reltol = max(reltol, 1.0e-15)

  x_ = 2π * x
  y_ = y

  p = FINNUFFTPlan(N, M, x_, y_, m, T(σ), reltol, fftflags)

  return p
end

function Base.show(io::IO, p::FINNUFFTPlan)
  print(io, "FINNUFFTPlan")
end

AbstractNFFTs.size_in(p::FINNUFFTPlan) = (Int.(p.N),)
AbstractNFFTs.size_out(p::FINNUFFTPlan) = (Int(p.M),)

function LinearAlgebra.mul!(fHat::StridedArray, p::FINNUFFTPlan{T}, f::AbstractArray;
 verbose=false, timing::Union{Nothing,TimingStats} = nothing) where {T}

  D = size(p.x,1)

  forwardPlan = FINUFFT.finufft_makeplan(3,D,-1,1,p.reltol;
                              nthreads = Threads.nthreads(),
                              fftw = p.fftflags)

  nodesX = ntuple(d->vec(p.x[d,:]), D)
  nodesY = ntuple(d->vec(p.y[d,:]), D)

  if D==1
    FINUFFT.finufft_setpts!(forwardPlan, nodesY..., T[], T[], nodesX..., T[], T[])
  elseif D==2
    FINUFFT.finufft_setpts!(forwardPlan, nodesY..., T[], nodesX..., T[])
  else
    FINUFFT.finufft_setpts!(forwardPlan, nodesY..., nodesX...)
  end

  FINUFFT.finufft_exec!(forwardPlan, f, fHat)

  FINUFFT.finufft_destroy!(forwardPlan)

  return fHat
end

function LinearAlgebra.mul!(f::StridedArray, pl::Adjoint{Complex{T},<:FINNUFFTPlan{T}}, fHat::AbstractArray;
         verbose=false, timing::Union{Nothing,TimingStats} = nothing) where {T}
  p = pl.parent

  D = size(p.x,1)

  adjointPlan = FINUFFT.finufft_makeplan(3,D,1,1,p.reltol;
              nthreads = Threads.nthreads(),
              fftw = p.fftflags)

  nodesX = ntuple(d->vec(p.x[d,:]), D)
  nodesY = ntuple(d->vec(p.y[d,:]), D)

  if D==1
    FINUFFT.finufft_setpts!(adjointPlan, nodesX..., T[], T[], nodesY..., T[], T[])
  elseif D==2
    FINUFFT.finufft_setpts!(adjointPlan, nodesX..., T[], nodesY...,T[])
  else
    FINUFFT.finufft_setpts!(adjointPlan, nodesX..., nodesY...)
  end

  FINUFFT.finufft_exec!(adjointPlan, fHat, f)

  FINUFFT.finufft_destroy!(adjointPlan)

  return f
end


text:
include("sysid.jl")

if !isdir("data")
  print("Downloading data...")
  run(`./get-data.sh`)
  run(`mkdir output`)
  println("done.")
end

function runExperiment(filename,tau)
  const nTrain = 300
  const nTest = 700
  const nTotal = nTrain+nTest
  dat = readdlm("data/$filename.dat")
  const u  = vec(dat[:,1])[1:nTotal]
  const y  = vec(dat[:,2])[1:nTotal]

  const uTrain = u[1:nTrain]
  const yTrain = y[1:nTrain]
  const yTest = y[(nTrain+1):end]

  function normalization(r,theta)
    1.0-r*r
  end

  function dndr(r,theta)
    -2*r
  end

  function dndt(r,theta)
    0.0
  end

  sim = LinearSysID(uTrain,100,50,normalization,dndr,dndt)

  target = yTrain

  sim_t = LinearSysID(u,2,2,normalization,dndr,dndt)
  function fullResponse(t,w,u)
    response = phi(sim_t,t,w)
    return response[(nTrain+1):end]
  end
  global scores = zeros(0)

  function printRMSE(old_thetas, t,w,k,o)
    est = phi(sim,t,w)
    rmse = norm(vec(target) - vec(est))
    fR = fullResponse(t,w,u)
    testRMSE = norm(vec(yTest) - vec(fR))
    mean_test = mean(yTest)
    println()
    println(size(t,2), " atoms.")
    response = phi(sim_t,t,w)
    score = 100*(1.0 - testRMSE/norm(yTest.-mean_test))
    global scores
    scores = [scores; score]
    println("Score: ", score)
    println("Tau: ", sum(abs(w)))
    println("thetas: ", t)
    println("weights: ", w)
    println(" ")
    response = phi(sim,t,w)
    return false
  end

  #fwami
  global scores = zeros(0)
  result = ADCG(sim, LSLoss(), target, tau; max_iters = 15, callback = printRMSE)
  writecsv("output/$(filename)_fwami.csv",[["iter" "score"]; [1:length(scores) scores]])

  #fwa-m
  global scores = zeros(0)
  result = ADCG(sim, LSLoss(), target, tau; max_cd_iters=0, max_iters = 15, callback = printRMSE)
  writecsv("output/$(filename)_fwam.csv",[["iter" "score"]; [1:length(scores) scores]])

  #gradient flow
  global scores = zeros(0)
  result = ADCG(sim, LSLoss(), target, tau; max_cd_iters=1, max_iters = 15, callback = printRMSE)
  writecsv("output/$(filename)_gf.csv",[["iter" "score"]; [1:length(scores) scores]])
end

runExperiment("dryer",20.0)
runExperiment("robot_arm",47.0)


text:
using TopOptGridap
using Documenter

DocMeta.setdocmeta!(TopOptGridap, :DocTestSetup, :(using TopOptGridap); recursive=true)

makedocs(;
    modules=[TopOptGridap],
    authors="Mohamed Tarek <mohamed82008@gmail.com> and contributors",
    repo="https://github.com/JuliaTopOpt/TopOptGridap.jl/blob/{commit}{path}#{line}",
    sitename="TopOptGridap.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaTopOpt.github.io/TopOptGridap.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaTopOpt/TopOptGridap.jl",
)


text:
#= make_scenarios_rejection_sample.jl
make timeseries of DR performance with a mean performace near 100%
Give a set of outcomes that each have an associated probability
And information on how often to resample from those outcomes

inputs:
n of total timesteps
prob distribution of performance and outcomes, vdr, p (original)
number of timesteps in which DR performance is constant
should scenarios be randomly selected
how many scenarios are selected?
threshold for rejection

outputs:
vdr[t,o]
p[o]
=#
using StatsBase
function make_scenarios(n_timesteps,v_og,p_og,int_length; randsel = true, nrand = 5, thresh = 0.005)
    # determine total number of possible scenarios
    Tp = convert(Int64,floor(n_timesteps/int_length)) # number of time blocks
    n_omega = length(p_og) # number of possilbe outcomes
    n_new_scenarios = convert(Int64,n_omega^Tp)

    # make vdr placeholder and all_s
    vdr = fill(0.0,(n_timesteps,nrand))
    p = fill(0.0,(nrand))
    all_s = collect(1:n_new_scenarios)

    # rejection sample
    i=1;
    n_sampled = 0
    while i<n_new_scenarios && n_sampled < nrand
        # draw a single sample
        s = sample(all_s,1)
        println("sampled ",s)

        # make vdr data for that sample
        s_vdr = fill(0.0,(n_timesteps))
        s_p = fill(0.0,(n_timesteps))
        for j in 1:Tp #row chunk
            # identify which rows belong to this time period
            if j == Tp
                last_period = n_timesteps
            else
                last_period = int_length*j #t_firsts[i+1]-1
            end
            rows = (int_length*(j-1) + 1):last_period

            # sample correct performance for that time period
            # imagine all possible combos are listed out in an array that is Tp x n_new_scenarios
            ind = mod(fld(s[1]-1,n_omega^(Tp-j)),n_omega)+1
                # fld: largest integer <= x/y
                # n_omega^(Tp-j) is the number of row-wise 'blocks' of the same response in our imagined array
                #           on the last block, there is just one block (j=Tp, this = 1)
                #           on the first block,   (j=1, this = n_omega^(Tp-1), e.g. n of blocks if every n_omega blocks were grouped)
                # divide sample number by above;
                # find remainer of that with number of possible responses
                # add 1 since above is 0-indexed
            s_vdr[rows] = v_og[ind]
            s_p[rows] = p_og[ind]
        end

        # assess if s_vdr is ok
        if abs(1-mean(s_vdr,weights(s_p))) < thresh #mean is weighted by probability
            # add to vdr and p
            vdr[:,n_sampled+1] = s_vdr
            p[n_sampled+1] = prod(s_p)
            # update n_sampled
            n_sampled = n_sampled+1
            println("added sample ",s)
        end

        # remove from possible scenarios to sample
        deleteat!(all_s,find(all_s .== s))

        i=i+1;
    end
    return(vdr,p)
end

max_stars_repo_name: pajalevi/Stoch-UC-model-for-DR

max_issues_repo_name: pajalevi/Stoch-UC-model-for-DR

max_forks_repo_name: pajalevi/Stoch-UC-model-for-DR


text:
using Telperion
using Test

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


text:
module RectModes
import StaticArrays: SVector
using Reexport
@reexport using Luna.Modes
import Luna: Maths
import Luna.PhysData: c, ε_0, μ_0, ref_index_fun, roomtemp
import Luna.Modes: AbstractMode, dimlimits, neff, field, Aeff, N

export RectMode, dimlimits, neff, field, N, Aeff

# core and clad are function-like objects which return the
# (possibly complex) refractive index as a function of freq
# pol is either :x or :y
# a and b are the half widths of the waveguide in each dimension.
struct RectMode{Ta, Tb, Tcore, Tclad} <: AbstractMode
    a::Ta
    b::Tb
    n::Int
    m::Int
    pol::Symbol
    coren::Tcore
    cladn::Tclad
end

function show(io::IO, m::RectMode)
    a = "a(z=0)=$(m.a(0))"
    b = "b(z=0)=$(m.b(0))"
    pol = string(m.pol)
    out = "RectMode{"*join([mode_string(m), a, b, "pol=$pol"], ", ")*"}"
    print(io, out)
end

mode_string(m::RectMode) = string(m.n)*string(m.m)

RectMode(a::Number, args...; kwargs...) = RectMode(z->a, args...; kwargs...)
RectMode(afun, b::Number, args...; kwargs...) = RectMode(afun, z->b, args...; kwargs...)
RectMode(a::Number, b::Number, args...; kwargs...) = RectMode(z->a, z->b, args...; kwargs...)

"convenience constructor assunming single gas filling and specified cladding"
function RectMode(afun, bfun, gas, P, clad; n=1, m=1, pol=:x, T=roomtemp)
    rfg = ref_index_fun(gas, P, T)
    rfs = ref_index_fun(clad)
    coren = (ω; z) -> rfg(2π*c./ω)
    cladn = (ω; z) -> rfs(2π*c./ω)
    RectMode(afun, bfun, n, m, pol, coren, cladn)
end

"convenience constructor for non-constant core index"
function RectMode(afun, bfun, coren, clad; n=1, m=1, pol=:x)
    rfs = ref_index_fun(clad)
    cladn = (ω; z) -> rfs(2π*c./ω)
    RectMode(afun, bfun, n, m, pol, coren, cladn)
end

dimlimits(m::RectMode; z=0) = (:cartesian, (-m.a(z), -m.b(z)), (m.a(z), m.b(z)))

"""
    neff(m::RectMode, ω; z=0)

effective index of rectangular mode with dielectric core and arbitrary
 (metal or dielectric) cladding.

Adapted from
Laakmann, K. D. & Steier, W. H.
Waveguides: characteristic modes of hollow rectangular dielectric waveguides.
Appl. Opt., AO 15, 1334–1340 (1976).

I had to re-derive the result in order to get the complex cladding index contribution
to the real part of neff.
"""
function neff(m::RectMode, ω; z=0)
    εcl = m.cladn(ω, z=z)^2
    εco = m.coren(ω, z=z)^2
    λ = 2π*c./ω
    if m.pol == :x
        ac = εcl/sqrt(complex(εcl - 1))
        bc = 1/sqrt(complex(εcl - 1))
    elseif m.pol == :y
        bc = εcl/sqrt(complex(εcl - 1))
        ac = 1/sqrt(complex(εcl - 1))
    else
        error("RectMode pol must be either :x or :y")
    end
    sqrt(complex(εco - (m.m*λ/(4*m.a(z)))^2*(1 - im*λ/(2π*m.a(z))*ac)^2
                     - (m.n*λ/(4*m.b(z)))^2*(1 - im*λ/(2π*m.b(z))*bc)^2))
end

# here we use cartesian coords, so xs = (x, y)
function field(m::RectMode, xs; z=0)
    if isodd(m.m)
        Ea = cos(m.m*π*xs[1]/(2*m.a(z)))
    else
        Ea = sin(m.m*π*xs[1]/(2*m.a(z)))
    end
    if isodd(m.n)
        Eb = cos(m.n*π*xs[2]/(2*m.b(z)))
    else
        Eb = sin(m.n*π*xs[2]/(2*m.b(z)))
    end
    E = Ea*Eb
    if m.pol == :x
        return SVector(E, 0.0)
    else
        return SVector(0.0, E)
    end
end

N(m::RectMode; z=0) = 0.5*sqrt(ε_0/μ_0)*m.a(z)*m.b(z)

Aeff(m::RectMode; z=0) = 16/9*m.a(z)*m.b(z)

end


text:
"""
    LU solver
"""
struct LUSolver{T} <: LinearSolver
    A::Array{T, 2}
end

function lu_solver(A)
    LUSolver(zero(Array(A)))
end

function linear_solve!(solver::LUSolver{T}, x::Vector{T}, A::SparseMatrixCSC{T,Int}, b::Vector{T}) where T
    solver.A .= A # sparse -> dense
    ldiv!(x, lu!(solver.A, check = false), b) # solve
end

function linear_solve!(solver::LUSolver{T}, x::Matrix{T}, A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T
    solver.A .= A # sparse -> dense
    ldiv!(x, lu!(solver.A, check = false), B) # solve
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NetworkSecurityRulesEvaluationResult <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    protocolMatched::Any # spec type: Union{ Nothing, Bool } # spec name: protocolMatched
    sourceMatched::Any # spec type: Union{ Nothing, Bool } # spec name: sourceMatched
    sourcePortMatched::Any # spec type: Union{ Nothing, Bool } # spec name: sourcePortMatched
    destinationMatched::Any # spec type: Union{ Nothing, Bool } # spec name: destinationMatched
    destinationPortMatched::Any # spec type: Union{ Nothing, Bool } # spec name: destinationPortMatched

    function NetworkSecurityRulesEvaluationResult(;name=nothing, protocolMatched=nothing, sourceMatched=nothing, sourcePortMatched=nothing, destinationMatched=nothing, destinationPortMatched=nothing)
        o = new()
        validate_property(NetworkSecurityRulesEvaluationResult, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(NetworkSecurityRulesEvaluationResult, Symbol("protocolMatched"), protocolMatched)
        setfield!(o, Symbol("protocolMatched"), protocolMatched)
        validate_property(NetworkSecurityRulesEvaluationResult, Symbol("sourceMatched"), sourceMatched)
        setfield!(o, Symbol("sourceMatched"), sourceMatched)
        validate_property(NetworkSecurityRulesEvaluationResult, Symbol("sourcePortMatched"), sourcePortMatched)
        setfield!(o, Symbol("sourcePortMatched"), sourcePortMatched)
        validate_property(NetworkSecurityRulesEvaluationResult, Symbol("destinationMatched"), destinationMatched)
        setfield!(o, Symbol("destinationMatched"), destinationMatched)
        validate_property(NetworkSecurityRulesEvaluationResult, Symbol("destinationPortMatched"), destinationPortMatched)
        setfield!(o, Symbol("destinationPortMatched"), destinationPortMatched)
        o
    end
end # type NetworkSecurityRulesEvaluationResult

const _property_map_NetworkSecurityRulesEvaluationResult = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("protocolMatched")=>Symbol("protocolMatched"), Symbol("sourceMatched")=>Symbol("sourceMatched"), Symbol("sourcePortMatched")=>Symbol("sourcePortMatched"), Symbol("destinationMatched")=>Symbol("destinationMatched"), Symbol("destinationPortMatched")=>Symbol("destinationPortMatched"))
const _property_types_NetworkSecurityRulesEvaluationResult = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("protocolMatched")=>"Bool", Symbol("sourceMatched")=>"Bool", Symbol("sourcePortMatched")=>"Bool", Symbol("destinationMatched")=>"Bool", Symbol("destinationPortMatched")=>"Bool")
Base.propertynames(::Type{ NetworkSecurityRulesEvaluationResult }) = collect(keys(_property_map_NetworkSecurityRulesEvaluationResult))
Swagger.property_type(::Type{ NetworkSecurityRulesEvaluationResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkSecurityRulesEvaluationResult[name]))}
Swagger.field_name(::Type{ NetworkSecurityRulesEvaluationResult }, property_name::Symbol) =  _property_map_NetworkSecurityRulesEvaluationResult[property_name]

function check_required(o::NetworkSecurityRulesEvaluationResult)
    true
end

function validate_property(::Type{ NetworkSecurityRulesEvaluationResult }, name::Symbol, val)
end


text:
"""
    DynaAgent{Tl<:AbstractLearner, Tb<:AbstractTurnBuffer, Tm<:AbstractEnvironmentModel, Tpp<:Function} <: AbstractAgent
    DynaAgent(learner::Tl, buffer::Tb, model::Tm, nsteps::Int=0, preprocessor::Tpp=identity, role=:anonymous) where {Tl, Tb, Tm, Tpp}

See more details at Section (8.2) on Page 162 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*

See also: [`Agent`](@ref), [`AbstractLearner`](@ref), [`AbstractTurnBuffer`](@ref), [`AbstractEnvironmentModel`](@ref)
"""
struct DynaAgent{Tl<:AbstractLearner,
                 Tb<:AbstractTurnBuffer,
                 Tm<:AbstractEnvironmentModel,
                 Tpp<:Function} <: AbstractAgent
    learner::Tl
    buffer::Tb
    model::Tm
    nsteps::Int
    preprocessor::Tpp
    role::Symbol
    function DynaAgent(learner::Tl,
                       buffer::Tb,
                       model::Tm,
                       nsteps::Int=0,
                       preprocessor::Tpp=identity,
                       role=:anonymous) where {Tl, Tb, Tm, Tpp}
        new{Tl, Tb, Tm, Tpp}(learner, buffer, model, nsteps, preprocessor, role)
    end
end

function update!(agent::DynaAgent)
    update!(agent.learner, agent.buffer)  # direct RL
    update!(agent.model, agent.buffer, agent.learner)  # model learning
    plan!(agent.learner, agent.model, agent.nsteps)  # planning
end

"""
    (agent::DynaAgent)(obs)

Take in an `obs` from environment and use `agent.preprocessor` to transform it into an internal `state`.
Then use `agent.learner` to get the `action`.
Return a pair of `state => action`.
"""
function (agent::DynaAgent)(obs)
    s = agent.preprocessor(obs)
    a = agent.learner(s)
    s => a
end

"""
    plan!(learner::QLearner, model::Union{TimeBasedSampleModel, ExperienceSampleModel}, nsteps::Int)

See more details at Section (8.2) on Page 164 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
function plan!(learner::QLearner, model::Union{TimeBasedSampleModel, ExperienceSampleModel}, nsteps::Int)
    for _ in 1:nsteps
        update!(learner, sample(model)...)
    end
end

"""
    plan!(learner::QLearner, model::PrioritizedSweepingSampleModel, nsteps::Int)

See more details at Section (8.2) on Page 170 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
function plan!(learner::QLearner, model::PrioritizedSweepingSampleModel, nsteps::Int)
    for _ in 1:nsteps
        record = sample(model)
        record === nothing && continue
        s, a, r, d, s′ = record
        update!(learner, s, a, r, d, s′)
        for (s̄, ā, r̄, d̄) in model.predecessors[s]
            P = priority(learner, s̄, ā, r̄, d̄, s)
            if P ≥ model.θ
                model.PQueue[(s̄, ā)] = P
            end
        end
    end
end


text:
"""
    div0(a, b)

Compute `a / b`, but return 0 if `b` is 0.
"""
function div0(
    dividend::Number,
    divisor::Number
)

    tmp = dividend / divisor
    return divisor == 0 ? zero(tmp) : tmp

end


text:
include("model.jl")

using Serialization
using Random
using SparseArrays
using Printf
using Statistics
using .Model

println("Loading data...")
X = deserialize("lastfm.jlser")
interactions = Interactions(X)

function select_k(scores::Matrix{<:AbstractFloat}, k::Int)
    n = size(scores, 2)
    y = Array{Int}(undef, k, n)
    ix = Array{Int}(undef, size(scores, 1))
    for i=1:n
        partialsortperm!(ix, view(scores, :, i), k, rev=true)
        y[:, i] = view(ix, 1:k)
    end
    return y
end

function precision_at_k(model::Factorization,
                        interactions::SparseMatrixCSC,
                        users::Vector{Int},
                        k::Int)
    scores = model.item_bias .+ model.item_embeddings' * model.user_embeddings[:, users]
    recommendations = select_k(scores, k)
    pk = zeros(Float32, length(users))
    for i=1:size(recommendations, 2)
        for j=1:size(recommendations, 1)
            user = users[i]
            item = recommendations[j, i]
            if is_positive(interactions, user, item)
                pk[i] += 1.0f0
            end
        end
    end
    return pk / k
end

println("Creating and initializing model")
model = Factorization{Float32}(size(X, 2), size(X, 1), 32)
optimizer = FactorizationAdadelta(model, 1e-6, 0.95)

users_val = abs.(rand(Int, 1000)) .% size(interactions.sparse, 2) .+ 1
println("Training")
init!(interactions, model)
for i=1:5
    t1 = time_ns()
    epoch!(model, optimizer, interactions, 1000)
    t2 = time_ns()
    @printf("Epoch time: %.2fmin ", (t2 - t1)/1e9/60)
    @printf("precision@5: %.3f\n",
            mean(precision_at_k(model, interactions.sparse, users_val, 5)))
end

max_stars_repo_name: plafl/julia-cython-comp

max_issues_repo_name: plafl/julia-cython-comp

max_forks_repo_name: plafl/julia-cython-comp


text:
# 0.2 short presentation for solving a model
expr_RunModel = :(begin;
   include("$(pwd())\\src\\proc_VarsDeclare.jl"); # use pwd() to avoid possible ambiguity of working directory
   include("$(pwd())\\src\\proc_InitPars.jl");
   Guess = ( r = 0.08, L = 0.2 );EasySearch.SteadyState!( 1, Guess, Dt, Dst, Pt, Ps, Pc, env, atol = 1E-8, MaxIter = 1000, PrintMode = "silent", MagicNum = 2.0, StepLen = 0.5 ); Guess = ( r = 0.12, L = 0.75 );
   EasySearch.SteadyState!( env.T, Guess, Dt, Dst, Pt, Ps, Pc, env, atol = 1E-6, MaxIter = 1000, PrintMode = "silent", MagicNum = 2.0, StepLen = 0.5 );
   PerfLog = EasySearch.Transition!( Dt, Dst, Pt, Ps, Pc, env, atol = 1E-3, MaxIter = 500, PrintMode = "silent", MagicNum = 2.0, StepLen = 0.7, ReturnLog = true );
   EasySearch.ProcAfterTransition!( Dt, Dst, Pt, Ps, Pc, env );
   end;
)
# 0.3 short presentation for extracting current demographic data
expr_GetDemogDat = :(
   Dict( "TotalPopu" => Dt[:N],
      "LaborPopu" => sum(Ps[:N][:,1:env.Sr],dims=2)[:],
      "WorkPopuRatio" => sum(Ps[:N][:,1:env.Sr],dims=2)[:] ./ sum(Ps[:N][:,1:env.S],dims=2)[:] ,
      "AgingPopu_65plus" => sum(Ps[:N][:, (65 - env.START_AGE + 1):env.S],dims=2)[:],
      "AgingPopuRatio" => sum(Ps[:N][:, (65 - env.START_AGE + 1):env.S],dims=2)[:] ./ sum(Ps[:N][:,1:env.S],dims=2)[:]
   )
)
# 0.4 a small function to reset env::NamedTuple
f_resetEnv( demogcsvpath::String ) = (
   T = 400, S = 80, Sr = 40, START_AGE = 20, START_YEAR = 1945,
   PATH_DEMOGCSV = demogcsvpath,
   PATH_WAGEPROFILE = "./data/WageProfileCoef.csv", PATH_MA2MB = "./data/MA2MBCoef.csv",
   PATH_M2C = "./data/M2C.csv", PATH_TFPGROW = "./data/tfpGrowthProfile.csv",
)
# 0.5 year range to plot
idx_year2plot = 2010:2110
idx_plot = idx_year2plot .- env.START_YEAR .+ 1
# -----------------------------------
# 0.6 run models under: base, low, high fertility assumptions, and save results
# NOTE: we will use these data to do plotting and other computations
   # 2.0.1 base (re-run it in case of not running SECTION 1)
   env = f_resetEnv( "./data/Demography_base.csv" )
   eval(expr_RunModel)
   dict_Demog_base = eval(expr_GetDemogDat)
   DatPkg_base = ( Dt = Dt, Dst = Dst, Pt = Pt, Ps = Ps, Pc = Pc )
   # 2.0.2 low fertility
   env = f_resetEnv( "./data/Demography_lowfertility.csv" )
   eval(expr_RunModel)
   dict_Demog_low = eval(expr_GetDemogDat)
   DatPkg_low = ( Dt = Dt, Dst = Dst, Pt = Pt, Ps = Ps, Pc = Pc )
   # 2.0.3 high fertility
   env = f_resetEnv( "./data/Demography_highfertility.csv" )
   eval(expr_RunModel)
   dict_Demog_high = eval(expr_GetDemogDat)
   DatPkg_high = ( Dt = Dt, Dst = Dst, Pt = Pt, Ps = Ps, Pc = Pc )

#


text:
#############################################################################
# nuclearnorm.jl
# Handles nuclear norm (the sum of the singular values of a matrix),
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################
export nuclearnorm

### Nuclear norm

struct NuclearNormAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int, Int}

    function NuclearNormAtom(x::AbstractExpr)
        children = (x,)
        return new(:nuclearnorm, hash(children), children, (1,1))
    end
end

function sign(x::NuclearNormAtom)
    return Positive()
end

# The monotonicity
function monotonicity(x::NuclearNormAtom)
    return (NoMonotonicity(),)
end

function curvature(x::NuclearNormAtom)
    return ConvexVexity()
end

function evaluate(x::NuclearNormAtom)
    return sum(svdvals(evaluate(x.children[1])))
end

nuclearnorm(x::AbstractExpr) = NuclearNormAtom(x)

# Create the equivalent conic problem:
#   minimize (tr(U) + tr(V))/2
#   subject to
#            [U A; A' V] ⪰ 0
# see eg Recht, Fazel, Parillo 2008 "Guaranteed Minimum-Rank Solutions of Linear Matrix Equations via Nuclear Norm Minimization"
# http://arxiv.org/pdf/0706.4138v1.pdf
function conic_form!(x::NuclearNormAtom, unique_conic_forms)
    if !has_conic_form(unique_conic_forms, x)
        A = x.children[1]
        m, n = size(A)
        U = Variable(m,m)
        V = Variable(n,n)
        p = minimize(.5*(tr(U) + tr(V)), [U A; A' V] ⪰ 0)
        cache_conic_form!(unique_conic_forms, x, p)
    end
    return get_conic_form(unique_conic_forms, x)
end


text:
module MCDTS

using LinearAlgebra, DataFrames, GLM, Distances, Statistics, StatsBase, Distributions, Neighborhood, DynamicalSystemsBase, DelayEmbeddings
using Random

import Base.show

export mcdts_embedding

include("optim.jl")
include("tree.jl")
include("./Delay preselection statistics/DelaySelectionStatistics.jl")
include("./Cost functions/CCM.jl")
include("./Cost functions/FNN_statistic.jl")
include("./Cost functions/L_statistic.jl")
include("./Cost functions/Prediction_error.jl")
include("tree_computations.jl")
include("utils.jl")

end # module


text:
##
##
## String representations
##
##

function show(io::IO, d::Distribution)
  print(io, @sprintf "%s distribution\n" typeof(d))
  for parameter in typeof(d).names
    if isa(d.(parameter), AbstractArray)
      param = string(ucfirst(string(parameter)), ":\n", d.(parameter), "\n")
    else
      param = string(ucfirst(string(parameter)), ": ", d.(parameter), "\n")
    end
    print(io, param)
  end
  m = mean(d)
  if isa(m, AbstractArray)
    print(io, string("Mean:\n", m, "\n"))
  else
    print(io, string("Mean: ", m, "\n"))
  end
  v = var(d)
  if isa(v, AbstractArray)
    print(io, string("Variance:\n", v))
  else
    print(io, string("Variance: ", v))
  end
end


text:
using Statistics
using LinearAlgebra
using Flux
using Flux: @functor
using LaTeXStrings
using KrylovMethods
using Distributions: Normal
using BSON: @load
using Plots
using CSV, DataFrames
using Dates
using Random

use_gpu = true
if use_gpu == true
    using CUDA
    #CUDA.allowscalar(false)
    cgpu = gpu
else
    cgpu = cpu
    pyplot()
end

pu = cpu # gpu
r_type = Float64
c_type = ComplexF64
u_type = Float32
gmres_type = ComplexF64
a_type = CuArray{gmres_type}
a_type = Array{gmres_type}

include("../src/multigrid/helmholtz_methods.jl")
include("../src/unet/model.jl")
include("../src/unet/data.jl")
include("../src/unet/train.jl")
include("../src/kappa_models.jl")
include("test_utils.jl")

fgmres_func = KrylovMethods.fgmres # gpu_flexible_gmres #

function test_train_unet!(n, f, opt, init_lr, train_size, test_size, batch_size, iterations;
                                    is_save=false, data_augmentetion=false, e_vcycle_input=false,
                                    kappa_type=1, threshold=50, kappa_input=true, kappa_smooth=false, k_kernel=3,
                                    gamma_input=true, kernel=(3,3), smaller_lr=10, v2_iter=10, level=3,
                                    axb=false, norm_input=false, model_type=SUnet, k_type=NaN, resnet_type=SResidualBlock, k_chs=-1, indexes=3, data_path="", full_loss=false, residual_loss=false, gmres_restrt=1, σ=elu, arch=1)

    h = 1.0./n;
    gamma_val = 0.00001
    pad_cells = [10;10]
    kappa = r_type.(ones(r_type,n-1,n-1)|>pu)
    omega = r_type(2*pi*f);
    gamma = gamma_val*2*pi * ones(r_type,size(kappa))
    gamma = r_type.(absorbing_layer!(gamma, pad_cells, omega))|>pu
    test_name = replace("$(Dates.format(now(), "HH_MM_SS")) RADAM ND $(model_type) $(k_type) $(resnet_type) $(k_chs) $(σ) $(indexes) $(k_kernel) g=$(gmres_restrt) t=$(u_type) g=$("$(gamma_input)"[1]) e=$("$(e_vcycle_input)"[1]) r=$("$(residual_loss)"[1]) k=$(kappa_type) $(threshold) n=$(n) f=$(f) m=$(train_size) bs=$(batch_size) lr=$(init_lr) each=$(smaller_lr) i=$(iterations)","."=>"_")
    model = create_model!(e_vcycle_input, kappa_input, gamma_input; kernel=kernel, type=model_type, k_type=k_type, resnet_type=resnet_type, k_chs=k_chs, indexes=indexes, σ=σ, arch=arch)|>cgpu
    model, train_loss, test_loss = train_residual_unet!(model, test_name, n, n, f, kappa, omega, gamma,
                                                        train_size, test_size, batch_size, iterations, init_lr;
                                                        e_vcycle_input=e_vcycle_input, v2_iter=v2_iter, level=level, data_augmentetion=data_augmentetion,
                                                        kappa_type=kappa_type, threshold=threshold, kappa_input=kappa_input, kappa_smooth=kappa_smooth, k_kernel=k_kernel,
                                                        gamma_input=gamma_input, kernel=kernel, smaller_lr=smaller_lr, axb=axb, jac=false, norm_input=norm_input, model_type=model_type, k_type=k_type, k_chs=k_chs, indexes=indexes,
                                                        data_path=data_path, full_loss=full_loss, residual_loss=residual_loss, gmres_restrt=gmres_restrt,σ=σ)

    iter = range(1, length=iterations)
    p = plot(iter, train_loss, label="Train loss")
    plot!(iter, test_loss, label="Test loss")
    yaxis!("Loss", :log10)
    xlabel!("Iterations")
    savefig("test/unet/results/$(test_name) residual graph")

    # New Example Check

    (input,e_true) = generate_random_data!(1, n, n, kappa, omega, gamma; e_vcycle_input=e_vcycle_input, v2_iter=v2_iter, level=level,
                                                          kappa_type=kappa_type, threshold=threshold, kappa_input=kappa_input, kappa_smooth=kappa_smooth, axb=axb, norm_input=norm_input, gmres_restrt=gmres_restrt)[1]

    e_vcycle = input[:,:,1:2,:]
    r = input[:,:,end-1:end,:]
    if gamma_input == true
        input = cat(input, gamma, dims=3)
    end
    input = input|> cgpu
    model_result = model(input)
    model_result = model_result|>cpu

    if is_save == true
        heatmap(e_true[:,:,1,1], color=:grays)
        title!(L"e^{true} = x - \tilde x")
        savefig("test/unet_helmholtz/results/$(test_name) e_true")

        heatmap(r[:,:,1,1], color=:grays)
        title!(L"r = b^{true} - A \tilde x")
        savefig("test/unet_helmholtz/results/$(test_name) r")

        heatmap(model_result[:,:,1,1], color=:grays)
        title!(L"e^{unet} = Unet(r)")
        savefig("test/unet_helmholtz/results/$(test_name) e_unet")
    end

    if e_vcycle_input == true
        if is_save == true
            heatmap(real(e_vcycle[:,:,1,1]), color=:grays)
            title!(L"e^{vcycle} = Vcycle(A,r,e^{(0)}=0)")
            savefig("test/unet_helmholtz/results/$(test_name) e_vcycle")
        end
        @info "$(Dates.format(now(), "HH:MM:SS")) - V-Cycle train error norm = $(norm_diff!(e_vcycle, e_true)), UNet train error norm = $(norm_diff!(model_result, e_true))"
    else
        @info "$(Dates.format(now(), "HH:MM:SS")) - UNet train error norm = $(norm_diff!(model_result, e_true))"
    end

    dataset_size = 10
    restrt = 10
    max_iter = 6
    check_model!(test_name, model, n, n, f, kappa, omega, gamma, e_vcycle_input, 4, kappa_input, gamma_input, kernel, dataset_size, restrt, max_iter; v2_iter=v2_iter, level=level, smooth=kappa_smooth, threshold=threshold, axb=axb, norm_input=norm_input, before_jacobi=false, indexes=indexes, arch=arch)
end

init_lr = 0.0001
opt = RADAM(init_lr)
train_size = 25000
test_size = 100
batch_size = 20
iterations = 120
full_loss = false
gmres_restrt = -1 # 1 -Default, 5 - 5GMRES, -1 Random

test_train_unet!(128, 10.0, opt, init_lr, train_size, test_size, batch_size, iterations;
                    data_augmentetion = false,
                    e_vcycle_input = false,
                    kappa_type = 1,
                    kappa_input = true,
                    threshold = 25,
                    kappa_smooth = true,
                    k_kernel = 5,
                    gamma_input = true,
                    kernel = (3,3),
                    smaller_lr = 48,
                    v2_iter = 10,
                    level = 3,
                    axb = false,
                    norm_input = false,
                    model_type = FFSDNUnet,
                    k_type = TFFKappa,
                    resnet_type = TSResidualBlockI,
                    k_chs = 10,
                    arch = 2,
                    indexes = 3,
                    full_loss = full_loss,
                    residual_loss = false,
                    data_path = "",
                    gmres_restrt = gmres_restrt,
                    σ = elu)


text:
function minimap2consensus(
    ref_fasta_path, query, out_bam_path, out_fasta_path;
    preset, min_ac, min_af, max_dp, threads=3)

    run(pipeline(`minimap2 -ax $preset -t $threads $ref_fasta_path $query`,
                 `samtools sort -O bam -o $out_bam_path -`))

    run(`samtools index $out_bam_path`)

    open(out_fasta_path, "w") do out_fasta_stream
        bam2pileup2consensus(out_bam_path, out_fasta_stream, ref_fasta_path,
                             min_ac, min_af, max_dp)
    end
end

"""
    align_short_reads(ref_fasta_path, fastq_list, out_bam_path, out_fasta_path; <keyword arguments>)

Align short reads with minimap2 & generate consensus FASTA.

Reads are aligned using `minimap2` with preset "sr". To generate a
consensus FASTA, a pileup is generated with `samtools mpileup`, and
the consensus allele is output at each site passing filters (see
keyword arguments); the missing allele "-" is output at sites not
passing filters. Note the consensus FASTA does not contain insertions,
only substitutions and deletions.

# Arguments
- `ref_fasta_path`: filename of reference FASTA. Should be gzipped and
   indexed with `samtools faidx`.
- `fastq_list`: list of fastq files for the sample.
- `out_bam_path`: Output BAM filename.
- `out_fasta_path`: Output consensus FASTA filename.
- `min_ac=10`: Filter for minimum count of the consensus allele.
- `min_af=0.9`: Filter for minimum frequency of the consensus allele.
- `max_dp=typemax(Int)`: Filter for maximum read depth to consider.
- `threads=3`: Number of threads for minimap2 (Note: pileup stage not currently parallelized)
"""
function align_short_reads(
    ref_fasta_path, fastq_list, out_bam_path, out_fasta_path;
    min_ac=10, min_af=0.9, max_dp=typemax(Int), threads=3)

    minimap2consensus(
        ref_fasta_path, fastq_list, out_bam_path, out_fasta_path,
        preset="sr", min_ac=min_ac, min_af=min_af, max_dp=max_dp,
        threads=threads)
end

"""
    align_assembly(ref_fasta_path, asm_fasta_path, out_bam_path, out_fasta_path; <keyword arguments>)

Align assembly with minimap2 & generate re-aligned FASTA.

Reads are aligned using `minimap2` (see keyword arguments for
settings). Then, a pileup is generated with `samtools mpileup`, and
the consensus allele is output at each site passing filters (see
keyword arguments); the missing allele "-" is output at sites not
passing filters. Note the consensus FASTA does not contain insertions,
only substitutions and deletions. Also note that the consensus FASTA
uses the contig names from the reference FASTA, not the new assembly.

# Arguments
- `ref_fasta_path`: filename of reference FASTA. Should be gzipped and
   indexed with `samtools faidx`.
- `asm_fasta_path`: filename of assembly FASTA.
- `out_bam_path`: Output BAM filename.
- `out_fasta_path`: Output consensus FASTA filename.
- `preset="asm5"`: Preset for minimap2. asm5 means <=5% divergence. Other relevant presets are asm10 and asm20.
- `min_ac=1`: Filter for minimum count of the consensus allele.
- `min_af=1`: Filter for minimum frequency of the consensus allele.
- `max_dp=1`: Filter for maximum read depth to consider.
- `threads=3`: Number of threads for minimap2 (Note: pileup stage not currently parallelized)
"""
function align_assembly(
    ref_fasta_path, asm_fasta_path, out_bam_path, out_fasta_path;
    preset="asm5", min_ac=1, min_af=1, max_dp=1, threads=3)

    minimap2consensus(
        ref_fasta_path, asm_fasta_path, out_bam_path, out_fasta_path,
        preset=preset, min_ac=min_ac, min_af=min_af, max_dp=max_dp,
        threads=threads)
end

export align_short_reads
export align_assembly


text:
module Spring #correspond à un namespace
    using Plots

    export SpringParam, InitParam, AnimParam
    export animate_spring

    toto(x) = cos(x)

    function gogo()
        x=[Float64(i) for i ∈ 1:60]
        y=toto.(x)
        plot(x,y)
    end

    Base.@kwdef struct SpringParam
        ls::Float64
        ms::Float64
        ks::Float64
        ns::Int
    end
    #A quoi sert Base.@kwdef? A rentrer des champs nommés
    getvalue(p::SpringParam) = p.ls,p.ms,p.ks,p.ns

    Base.@kwdef struct InitParam
        λ::Float64
        shift::Float64
        pos::Float64
    end
    getvalue(p::InitParam) = p.λ,p.shift,p.pos
    #meme nom de methode mais argument specialise

    Base.@kwdef struct AnimParam
        δt::Float64
        nδt::Int
        nδtperframe::Int
    end
    getvalue(p::AnimParam) = p.δt,p.nδt,p.nδtperframe

    steady_position(i,ls) = (i-1)*ls

    function initial_position(i,ls,λ,shift,pos)
        xs=steady_position(i,ls)
        dx=xs-pos
        xs-λ*dx*shift*exp(-0.5*dx^2/λ^2)
    end

    function update_force!(fx,xc,ks) # ! est une convention de nommage en julia
                                     # pour les fonctions dont on modifie les arguments
                                     # placés en premieres positions
        ns=length(xc)
        for i ∈ 2:ns-1
            fx[i] = -ks*(2xc[i]-xc[i-1]-xc[i+1])
        end
    end

    function update_position!(xt,xc,xp,fx,δt,ms)
        coef=δt^2/ms
        @. xt = 2*xc - xp + fx*coef
    end

    function advance_nδtpf(xc,xp,xt,fx,sp,ap)
        ls,ms,ks,ns=getvalue(sp)
        δt,nδt,nδtperframe=getvalue(ap)

        for _ ∈ 1:nδtperframe #boucle sans indice
            update_force!(fx,xc,ks)
            update_position!(xt,xc,xp,fx,δt,ms)
            xc,xp,xt=xt,xc,xp
        end
        xc,xp,xt
    end

    function animate_spring(sp,ip,ap)
        ls,ms,ks,ns=getvalue(sp)
        λ,shift,pos=getvalue(ip)
        δt,nδt,nδtperframe=getvalue(ap)

        xs=[steady_position(i,ls) for i ∈ 1:ns] #array comprehension
        xc=[initial_position(i,ls,λ,shift,pos) for i ∈ 1:ns]

        dc=zero(xc) #tableau de la meme forme et meme type que xc avec des 0
        #for i ∈ 1:ns
        #    dc[i]=xc[i]-xs[i]
        #end
        #version broadcastée, pas de difference de rapidité en julia
        #dc.=xc.-xs
        #autre notation pour ne pas mettre des . partout
        #@. dc=xc-xs

        fx=zero(xc)
        #update_force!(fx,xc,ks)

        xt=zero(xc)
        xp=copy(xc) #sinon ca pointe juste, on change de nom
        #example tester a=rand(2), b=rand(3), a,b=b,a
        nf=nδt÷nδtperframe
        t=0.0
        anim = @animate for i ∈ 1:nf #creation d'une animation
            xc,xp,xt=advance_nδtpf(xc,xp,xt,fx,sp,ap)
            @. dc = xc - xs
            t+=nδtperframe*δt
            plot(xs,dc,ylims=(-shift,shift),title="t=$t")
        end
        gif(anim,"toto.gif",fps=15)

        #plot(xs,fx)
        #plot(xs,dc)

        #imid=ns÷2 #division entiere \div
        #xs[imid],xc[imid],dc[imid]
        #sp,ip,ap #retour de la fonction main
    end

    # si on laisse main ici il faut l'appeler avec Spring.main() dans le REPL
    # ou l'exporter du module pour s'affranchir du Spring.
    # Mieux faire un fichier main.jl à la racine puis dans le REPL include("main.jl")
    #function main()
    #    sp=Spring.SpringParam(ls=0.1,ms=1,ks=1,ns=400)
    #    ip=InitParam(λ=1.0,shift=1.5,pos=sp.ls*sp.ns/2)
    #    ap=AnimParam(δt=1.0,nδt=1000,nδtperframe=10)

    #end

end # module


text:
using LinearAlgebra: dot, norm

function angular_error(x, y)
    x3 = (1, x...)
    y3 = (1, y...)
    z = dot(x3, y3) / sqrt(dot(x3, x3) * dot(y3, y3))
    # gracefully handle z > 1 in case of rounding errors
    return acos(clamp(z, 0., 1.))
end

function endpoint_error(x, y)
    return norm(x .- y)
end


text:
using LocalHillResponses
using Documenter

makedocs(;
    modules=[LocalHillResponses],
    authors="Wen-Wei Tseng <sosiristseng@gmail.com>",
    repo="https://github.com/NTUMitoLab/LocalHillResponses.jl/blob/{commit}{path}#L{line}",
    sitename="LocalHillResponses.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://NTUMitoLab.github.io/LocalHillResponses.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/NTUMitoLab/LocalHillResponses.jl.git",
)


text:
using Test, Dates

# will export `YearMonth` type
using YearMonths

# constructors
ym = YearMonth(2018, 12)
@test ym == YearMonth("2018-12") # accepts yyyy-mm format
@test ym == YearMonth("201812")  # accepts yyyymm format

# algebra
@test ym - Month(1) == YearMonth(2018, 11)
@test ym - Year(1) == YearMonth(2017, 12)
@test ym + Month(1) == YearMonth(2019, 1)
@test Month(1) + ym == YearMonth(2019, 1)
@test ym + Year(1) == YearMonth(2019, 12)
@test Year(1) + ym == YearMonth(2019, 12)
@test ym + Month(12) == YearMonth(2019, 12)

# comparison
ym1 = YearMonth(2020, 1)
ym2 = YearMonth(2019, 12)
@test ym1 > ym2

# Basic `Dates` API
ym = YearMonth(2018, 5)
@test firstdayofmonth(ym) == Date(2018, 5, 1)
@test lastdayofmonth(ym) == Date(2018, 5, 31)
@test month(ym) == 5
@test year(ym) == 2018
@test yearmonth(ym) == (2018, 5)


text:
using LinearAlgebra
import LinearAlgebra:checksquare

"""
    function sweep_op_internal!(A::AbstractMatrix{Float64}, k::Int64)

    Implement the naive sweep operator described by Goodnight, J. (1979). "A Tutorial on the SWEEP Operator." The American Statistician.
"""
function sweep_op_internal!(A::AbstractMatrix{Float64}, k::Int64)
    p = checksquare(A)
    if k >= p
        throw(ArgumentError("Incorrect k value"))
    end

    @inbounds D = A[k,k] # step 1 D = Aₖₖ

    if D == zero(eltype(A))
        throw(ArgumentError("sweep_op_internal!: the element $k,$k of the matrix is zero. Is the Design Matrix symmetric positive definite?"))
    end

    @inbounds for i in 1:p
        A[k,i] = A[k,i] / D   # step 2 divide row k by D
    end

    @inbounds for i in 1:(k - 1) # step 3: for every other row  i != k
        B = A[i,k]   # let B = Aᵢₖ
        for j in 1:p
            A[i,j] = A[i,j] - B * A[k, j]  # Subtract B * row k from row i
        end
        A[i,k] = - B / D # set Aᵢₖ = -B/D
    end

    @inbounds for i in k + 1:p # step 3: for every other row  i != k
        B = A[i,k]   # let B = Aᵢₖ
        for j in 1:p
            A[i,j] = A[i,j] - B * A[k, j]  # Subtract B * row k from row i
        end
        A[i,k] = - B / D # set Aᵢₖ = -B/D
    end

    @inbounds A[k,k] = 1 / D # step 4 Set Aₖₖ = 1/D

    return A
end

function sweep_op!(A::AbstractMatrix{Float64}, k::Int64)
    sweep_op_internal!(A, k)
end

function sweep_op!(A::AbstractMatrix{Float64}, ks::AbstractVector{Int64})
    for k in ks
        sweep_op_internal!(A, k)
    end
    return A
end

"""
    function sweep_op_full!(A::AbstractMatrix{Float64})

    (internal) Get SSE, error sum of squares, for the full model.
"""
function sweep_op_full!(A::AbstractMatrix{Float64})
    n , p = size(A)

    for k in 1:p-1
        sweep_op_internal!(A, k)
    end
    return A[p,p]
end

"""
    function sweep_op_fullT1SS!(A::AbstractMatrix{Float64})

    (internal) Get SSE, error sum of squares for the full model.
    Also give Type I SS for all independent variables.
"""
function sweep_op_fullT1SS!(A::AbstractMatrix{Float64})
    n , p = size(A)
    TypeISS = Vector{Float64}(undef, p-1)

    for k in 1:p-1
        preSSE = A[p,p]
        sweep_op_internal!(A, k)
        TypeISS[k] = preSSE - A[p,p]
    end
    return A[p,p], TypeISS
end

"""
    function sweep_op_fullT1SS!(A::AbstractMatrix{Float64})

    (internal) Get SSE, error sum of squares for all (p-1) models.
    Also give Type I SS for all independent variables.
"""
function sweep_op_allT1SS!(A::AbstractMatrix{Float64})
    n , p = size(A)
    TypeISS = Vector{Float64}(undef, p-1)
    SSEs = Vector{Float64}(undef, p-1)

    for k in 1:p-1
        preSSE = A[p,p]
        sweep_op_internal!(A, k)
        SSEs[k] = A[p, p]
        TypeISS[k] = preSSE - A[p,p]
    end
    return SSEs, TypeISS
end

"""
    function get_TypeIISS(extended_inverse::AbstractMatrix{Float64})

    (internal) Get Type II SS for all independent variables given an extended inverse (sweep operator already applied)
"""
function get_TypeIISS(extended_inverse::AbstractMatrix{Float64})
    n, p = size(extended_inverse)
    TypeIISS = Vector{Float64}(undef, p-1)

    for k in 1:p-1
        TypeIISS[k] = extended_inverse[k, end]^2 / extended_inverse[k, k]
    end

    return TypeIISS
end


text:
include("AlgAssRelOrd/Order.jl")
include("AlgAssRelOrd/Elem.jl")
include("AlgAssRelOrd/Ideal.jl")
include("AlgAssRelOrd/NEQ.jl")
include("AlgAssRelOrd/Eichler.jl")


text:
module HelloDocs

export GreetUser

#Exported function:

"""
    GreetUser()

Greet User.
"""
function GreetUser()
  print("Hello User!")
end

#NON-exported function:

"""
    greet()

Greet World.
"""
greet() = print("Hello World!")

#NON-exported function:

"""
    GreetYou()

Greet You.
"""
GreetYou() = print("Hello You!")

#

end # module


text:
using Test
using Markdown: MD

@testset "ReturnTypes" begin

    using CpuId

    # Moved upwards temporarily for better diagnostics
    println(cpuinfo())
    println(cpufeaturetable())
    println(hvinfo())
    println()

    # Can't do real testing on results when target machine is unknown.
    # Thus, let's simply check whether the result types are correct,
    # which also fails if a test throws.

    @test isa( CpuId.cpuid(), NTuple{4, UInt32} )
    @test isa( CpuId.cpuid(0x00), NTuple{4, UInt32} )
    @test isa( CpuId.cpuid(0x00, 0x00), NTuple{4, UInt32} )

    # Test the low-level cpufeature querying function
    @test isa( CpuId.cpufeature(CpuId.SSE)      , Bool )
    @test isa( CpuId.cpufeature(:SSE)           , Bool )
    @test isa( CpuId.cpufeaturedesc(:SSE)       , String )
    @test isa( CpuId.cpufeatures()              , Vector{Symbol} )

    @test_throws UndefVarError CpuId.cpufeature(:UNKNOWNFEATURE)
    @test CpuId.cpufeaturedesc(:UNKNOWNFEATURE) ==
                            "Unknown feature, no description available!"

    # LLVM eliminates calls to hasleaf(...) if the executing machine supports
    # that leaf.  Thus test whether the reverse actually throws...
    function test_nonexisting_leaf()
        leaf = 0x8000_008f
        CpuId.hasleaf(leaf) || CpuId._throw_unsupported_leaf(leaf)
        CpuId.cpuid(leaf)
    end
    @test_throws ErrorException test_nonexisting_leaf()

    @test isa( address_size()         , Integer )
    @test isa( physical_address_size(), Integer )
    @test isa( cachelinesize()        , Integer )
    @test isa( cachesize()            , Tuple )
    for i in 0:5
        @test isa( cachesize(i)       , Integer )
    end
    @test isa( cpuarchitecture()      , Symbol )
    @test isa( cpubrand()             , String )
    @test isa( cpumodel()             , Dict )
    @test isa( cpuvendor()            , Symbol )
    @test isa( CpuId.cpuvendorstring(), String )
    @test isa( hypervised()           , Bool )
    @test isa( hvvendor()             , Symbol )
    @test isa( CpuId.hvvendorstring() , String )
    @test isa( hvversion()            , Dict{Symbol,Any} )
    @test isa( hypervised()           , Bool )
    @test isa( simdbits()             , Integer )
    @test isa( simdbytes()            , Integer )
    @test isa( cpucores()             , Integer )
    @test isa( cputhreads()           , Integer )
    @test isa( has_cpu_frequencies()  , Bool )
    @test isa( cpu_base_frequency()   , Integer )
    @test isa( cpu_bus_frequency()    , Integer )
    @test isa( cpu_max_frequency()    , Integer )
    @test isa( cpuinfo()              , MD )
    @test isa( cpufeaturetable()      , MD )
    @test isa( hvinfo()               , MD )

    @test isa( cpucycle()             , UInt64 )
    @test isa( cpucycle_id()          , Tuple{UInt64,UInt64} )

    @test isa( perf_revision()        , Int )
    @test isa( perf_gen_counters()    , Int )
    @test isa( perf_gen_bits()        , Int )
    @test isa( perf_fix_counters()    , Int )
    @test isa( perf_fix_bits()        , Int )

    # Check if trailing null characters are correctly identified
    # as hypervisor vendor KVM
    @test get( CpuId._cpuid_vendor_id, "KVMKVMKVM\0\0\0", :Unknown) === :KVM

    # If we're on Linux, then also dump /proc/cpuinfo for comparison when on a
    # remote CI, but only print data of the first CPU.
    Sys.islinux() && run(`sed -e '/^$/,$d' /proc/cpuinfo`)

end

# Dump the cpuid table of the executing CPU
include("mock.jl")
include("mockdb.jl")

dump_cpuid_table() ; flush(stdout) ; flush(stderr)

print("\n\n-----\nMocking CpuId\n-----\n\n")
flush(stdout) ; flush(stderr)

# Run the known cpuid records
@testset "Mocking" begin
    for i in 1:length(_mockdb)
        # temporarily replace the low-level cpuid function with known records
        eval(quote
            @testset "Mocked #$($i) $(strip(cpubrand()))" begin
                mock_cpuid($i)
                flush(stdout) ; flush(stderr)
                @test isa( cpubrand()       , String )
                @test isa( cpuinfo()        , MD )
                @test isa( cpufeaturetable(), MD )
                @test isa( hvinfo()         , MD )
                println("Tested recorded cpuid table #",$i," for '", strip(cpubrand()), "'")

                for (fns, res) in last(_mockdb[$i])
                    if isa(res, Tuple) && length(res) > 0 && first(res) === :broken
                        @test_broken getfield(CpuId, fns)() == last(res)
                    else
                        @test getfield(CpuId, fns)() == res
                    end
                end

            end
        end)

        #dump_cpuid_table()
        flush(stdout) ; flush(stderr)
    end
end


text:
# Unit Test for by split vector type
using COSMO, Test, Random, LinearAlgebra

rng = Random.MersenneTwister(13131)

convex_set = COSMO.PsdConeTriangle(3)
composite_set = COSMO.CompositeConvexSet([COSMO.ZeroSet(3), COSMO.Nonnegatives(1)])
x1 = rand(3)
x2 = rand(4)

s1 = COSMO.SplitVector(x1, convex_set)
s2 = COSMO.SplitVector(x2, composite_set)
@test size(s1) == (3, )
@test length(s2) == 4
@test getindex(s1, 2) == x1[2]
setindex!(s1, 5, 2)
@test getindex(s1, 2) == 5.
@test firstindex(s2) == 1
@test lastindex(s2) == 4
@test Base.showarg(IOBuffer(), s2, true ) == nothing


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

include("utils.jl")
include("unet3D.jl")


text:
@testset "6.2.7 hyper^m (a+b cosh^n)^p" begin
    (a, b, n, x, ) = @variables a b n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*sinh(c+d*x)^m*(a+b*cosh(c+d*x)^n)^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*sinh(c+d*x)^m*(a+b*cosh(c+d*x)^2)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sinh(c+d*x)^m*(a+b*cosh(c+d*x)^2)^p*when*a+b=0=#

    #= ::Subsubsection:: =#
    #=p>0=#

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [sinh(x)^4/(a - a*cosh(x)^2), x, 3, x/(2*a) - (cosh(x)*sinh(x))/(2*a)]
    @test_int [sinh(x)^3/(a - a*cosh(x)^2), x, 2, -(cosh(x)/a)]
    @test_int [sinh(x)^2/(a - a*cosh(x)^2), x, 2, -x/a]
    @test_int [csch(x)^2/(a - a*cosh(x)^2), x, 3, -(coth(x)/a) + coth(x)^3/(3*a)]
    @test_int [csch(x)^4/(a - a*cosh(x)^2), x, 3, coth(x)/a - (2*coth(x)^3)/(3*a) + coth(x)^5/(5*a)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sinh(c+d*x)^m*(a+b*cosh(c+d*x)^2)^p=#

    #= ::Subsubsection:: =#
    #=p>0=#

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [sinh(x)^7/(a + b*cosh(x)^2), x, 4, -(((a + b)^3*atan((sqrt(b)*cosh(x))/sqrt(a)))/(sqrt(a)*b^(7/2))) + ((a^2 + 3*a*b + 3*b^2)*cosh(x))/b^3 - ((a + 3*b)*cosh(x)^3)/(3*b^2) + cosh(x)^5/(5*b)]
    @test_int [sinh(x)^5/(a + b*cosh(x)^2), x, 4, ((a + b)^2*atan((sqrt(b)*cosh(x))/sqrt(a)))/(sqrt(a)*b^(5/2)) - ((a + 2*b)*cosh(x))/b^2 + cosh(x)^3/(3*b)]
    @test_int [sinh(x)^3/(a + b*cosh(x)^2), x, 3, -(((a + b)*atan((sqrt(b)*cosh(x))/sqrt(a)))/(sqrt(a)*b^(3/2))) + cosh(x)/b]
    @test_int [sinh(x)^1/(a + b*cosh(x)^2), x, 2, atan((sqrt(b)*cosh(x))/sqrt(a))/(sqrt(a)*sqrt(b))]
    @test_int [csch(x)^1/(a + b*cosh(x)^2), x, 4, -((sqrt(b)*atan((sqrt(b)*cosh(x))/sqrt(a)))/(sqrt(a)*(a + b))) - atanh(cosh(x))/(a + b)]
    @test_int [csch(x)^3/(a + b*cosh(x)^2), x, 5, (b^(3/2)*atan((sqrt(b)*cosh(x))/sqrt(a)))/(sqrt(a)*(a + b)^2) + ((a + 3*b)*atanh(cosh(x)))/(2*(a + b)^2) - (coth(x)*csch(x))/(2*(a + b))]
    @test_int [csch(x)^5/(a + b*cosh(x)^2), x, 6, -((b^(5/2)*atan((sqrt(b)*cosh(x))/sqrt(a)))/(sqrt(a)*(a + b)^3)) - ((3*a^2 + 10*a*b + 15*b^2)*atanh(cosh(x)))/(8*(a + b)^3) + ((3*a + 7*b)*coth(x)*csch(x))/(8*(a + b)^2) - (coth(x)*csch(x)^3)/(4*(a + b))]

    @test_int [sinh(x)^6/(a + b*cosh(x)^2), x, 6, ((8*a^2 + 20*a*b + 15*b^2)*x)/(8*b^3) - ((a + b)^(5/2)*atanh((sqrt(a)*tanh(x))/sqrt(a + b)))/(sqrt(a)*b^3) - ((4*a + 7*b)*cosh(x)*sinh(x))/(8*b^2) + (cosh(x)*sinh(x)^3)/(4*b)]
    @test_int [sinh(x)^4/(a + b*cosh(x)^2), x, 5, -(((2*a + 3*b)*x)/(2*b^2)) + ((a + b)^(3/2)*atanh((sqrt(a)*tanh(x))/sqrt(a + b)))/(sqrt(a)*b^2) + (cosh(x)*sinh(x))/(2*b)]
    @test_int [sinh(x)^2/(a + b*cosh(x)^2), x, 4, x/b - (sqrt(a + b)*atanh((sqrt(a)*tanh(x))/sqrt(a + b)))/(sqrt(a)*b)]
    @test_int [sinh(x)^0/(a + b*cosh(x)^2), x, 2, atanh((sqrt(a)*tanh(x))/sqrt(a + b))/(sqrt(a)*sqrt(a + b))]
    @test_int [csch(x)^4/(a + b*cosh(x)^2), x, 4, (b^2*atanh((sqrt(a)*tanh(x))/sqrt(a + b)))/(sqrt(a)*(a + b)^(5/2)) + ((a + 2*b)*coth(x))/(a + b)^2 - coth(x)^3/(3*(a + b))]
    @test_int [csch(x)^6/(a + b*cosh(x)^2), x, 4, -((b^3*atanh((sqrt(a)*tanh(x))/sqrt(a + b)))/(sqrt(a)*(a + b)^(7/2))) - ((a^2 + 3*a*b + 3*b^2)*coth(x))/(a + b)^3 + ((2*a + 3*b)*coth(x)^3)/(3*(a + b)^2) - coth(x)^5/(5*(a + b))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*sinh(c+d*x)^m*(a+b*cosh(c+d*x)^3)^p=#

    @test_int [sinh(x)/(4 - 3*cosh(x)^3), x, 7, atan((1 + 6^(1/3)*cosh(x))/sqrt(3))/(2*2^(1/3)*3^(5/6)) - log(2^(2/3) - 3^(1/3)*cosh(x))/(6*6^(1/3)) + log(2*2^(1/3) + 2^(2/3)*3^(1/3)*cosh(x) + 3^(2/3)*cosh(x)^2)/(12*6^(1/3))]

    #= ::Section:: =#
    #=Integrands*of*the*form*sinh(c+d*x)^m*(a+b*cosh(c+d*x)^4)^p=#

    #= ::Section:: =#
    #=Integrands*of*the*form*sinh(c+d*x)^m*(a+b*cosh(c+d*x)^n)^p=#

    #= ::Title:: =#
    #=Integrands*of*the*form*cosh(c+d*x)^m*(a+b*cosh(c+d*x)^n)^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*cosh(c+d*x)^m*(a+b*cosh(c+d*x)^2)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cosh(c+d*x)^m*(a+b*cosh(c+d*x)^2)^p=#

    #= ::Subsubsection:: =#
    #=p>0=#

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [cosh(x)^7/(a + b*cosh(x)^2), x, 4, -((a^3*atan((sqrt(b)*sinh(x))/sqrt(a + b)))/(b^(7/2)*sqrt(a + b))) + ((a^2 - a*b + b^2)*sinh(x))/b^3 - ((a - 2*b)*sinh(x)^3)/(3*b^2) + sinh(x)^5/(5*b)]
    @test_int [cosh(x)^6/(a + b*cosh(x)^2), x, 6, ((8*a^2 - 4*a*b + 3*b^2)*x)/(8*b^3) - (a^(5/2)*atanh((sqrt(a)*tanh(x))/sqrt(a + b)))/(b^3*sqrt(a + b)) - ((4*a - 3*b)*cosh(x)*sinh(x))/(8*b^2) + (cosh(x)^3*sinh(x))/(4*b)]
    @test_int [cosh(x)^5/(a + b*cosh(x)^2), x, 4, (a^2*atan((sqrt(b)*sinh(x))/sqrt(a + b)))/(b^(5/2)*sqrt(a + b)) - ((a - b)*sinh(x))/b^2 + sinh(x)^3/(3*b)]
    @test_int [cosh(x)^4/(a + b*cosh(x)^2), x, 5, -(((2*a - b)*x)/(2*b^2)) + (a^(3/2)*atanh((sqrt(a)*tanh(x))/sqrt(a + b)))/(b^2*sqrt(a + b)) + (cosh(x)*sinh(x))/(2*b)]
    @test_int [cosh(x)^3/(a + b*cosh(x)^2), x, 3, -((a*atan((sqrt(b)*sinh(x))/sqrt(a + b)))/(b^(3/2)*sqrt(a + b))) + sinh(x)/b]
    @test_int [cosh(x)^2/(a + b*cosh(x)^2), x, 3, x/b - (sqrt(a)*atanh((sqrt(a)*tanh(x))/sqrt(a + b)))/(b*sqrt(a + b))]
    @test_int [cosh(x)^1/(a + b*cosh(x)^2), x, 2, atan((sqrt(b)*sinh(x))/sqrt(a + b))/(sqrt(b)*sqrt(a + b))]
    @test_int [cosh(x)^0/(a + b*cosh(x)^2), x, 2, atanh((sqrt(a)*tanh(x))/sqrt(a + b))/(sqrt(a)*sqrt(a + b))]
    @test_int [sech(x)^1/(a + b*cosh(x)^2), x, 4, atan(sinh(x))/a - (sqrt(b)*atan((sqrt(b)*sinh(x))/sqrt(a + b)))/(a*sqrt(a + b))]
    @test_int [sech(x)^2/(a + b*cosh(x)^2), x, 3, -((b*atanh((sqrt(a)*tanh(x))/sqrt(a + b)))/(a^(3/2)*sqrt(a + b))) + tanh(x)/a]
    @test_int [sech(x)^3/(a + b*cosh(x)^2), x, 5, ((a - 2*b)*atan(sinh(x)))/(2*a^2) + (b^(3/2)*atan((sqrt(b)*sinh(x))/sqrt(a + b)))/(a^2*sqrt(a + b)) + (sech(x)*tanh(x))/(2*a)]
    @test_int [sech(x)^4/(a + b*cosh(x)^2), x, 4, (b^2*atanh((sqrt(a)*tanh(x))/sqrt(a + b)))/(a^(5/2)*sqrt(a + b)) + ((a - b)*tanh(x))/a^2 - tanh(x)^3/(3*a)]
    @test_int [sech(x)^5/(a + b*cosh(x)^2), x, 6, ((3*a^2 - 4*a*b + 8*b^2)*atan(sinh(x)))/(8*a^3) - (b^(5/2)*atan((sqrt(b)*sinh(x))/sqrt(a + b)))/(a^3*sqrt(a + b)) + ((3*a - 4*b)*sech(x)*tanh(x))/(8*a^2) + (sech(x)^3*tanh(x))/(4*a)]

    @test_int [1/(a + b*cosh(x)^2)^2, x, 4, ((2*a + b)*atanh((sqrt(a)*tanh(x))/sqrt(a + b)))/(2*a^(3/2)*(a + b)^(3/2)) - (b*cosh(x)*sinh(x))/(2*a*(a + b)*(a + b*cosh(x)^2))]

    @test_int [1/(a + b*cosh(x)^2)^3, x, 5, ((8*a^2 + 8*a*b + 3*b^2)*atanh((sqrt(a)*tanh(x))/sqrt(a + b)))/(8*a^(5/2)*(a + b)^(5/2)) - (b*cosh(x)*sinh(x))/(4*a*(a + b)*(a + b*cosh(x)^2)^2) - (3*b*(2*a + b)*cosh(x)*sinh(x))/(8*a^2*(a + b)^2*(a + b*cosh(x)^2))]

    @test_int [1/(1 + cosh(x)^2), x, 2, atanh(tanh(x)/sqrt(2))/sqrt(2)]
    @test_int [1/(1 + cosh(x)^2)^2, x, 4, (3*atanh(tanh(x)/sqrt(2)))/(4*sqrt(2)) - (cosh(x)*sinh(x))/(4*(1 + cosh(x)^2))]
    @test_int [1/(1 + cosh(x)^2)^3, x, 5, (19*atanh(tanh(x)/sqrt(2)))/(32*sqrt(2)) - (cosh(x)*sinh(x))/(8*(1 + cosh(x)^2)^2) - (9*cosh(x)*sinh(x))/(32*(1 + cosh(x)^2))]

    @test_int [1/(1 - cosh(x)^2), x, 3, coth(x)]
    @test_int [1/(1 - cosh(x)^2)^2, x, 3, coth(x) - coth(x)^3/3]
    @test_int [1/(1 - cosh(x)^2)^3, x, 3, coth(x) - (2*coth(x)^3)/3 + coth(x)^5/5]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cosh(c+d*x)^m*(a+b*cosh(c+d*x)^2)^(p/2)=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [sqrt(a + b*cosh(x)^2), x, 2, -((im*sqrt(a + b*cosh(x)^2)*Elliptic.E(pi/2 + im*x, -(b/a)))/sqrt(1 + (b*cosh(x)^2)/a))]

    @test_int [sqrt(1 + cosh(x)^2), x, 1, (-im)*Elliptic.E(pi/2 + im*x, -1)]
    @test_int [sqrt(1 - cosh(x)^2), x, 3, coth(x)*sqrt(-sinh(x)^2)]
    @test_int [sqrt(-1 + cosh(x)^2), x, 3, coth(x)*sqrt(sinh(x)^2)]
    @test_int [sqrt(-1 - cosh(x)^2), x, 2, -((im*sqrt(-1 - cosh(x)^2)*Elliptic.E(pi/2 + im*x, -1))/sqrt(1 + cosh(x)^2))]

    @test_int [(a + b*cosh(x)^2)^(3/2), x, 6, -((2*im*(2*a + b)*sqrt(a + b*cosh(x)^2)*Elliptic.E(pi/2 + im*x, -(b/a)))/(3*sqrt(1 + (b*cosh(x)^2)/a))) + (im*a*(a + b)*sqrt(1 + (b*cosh(x)^2)/a)*Elliptic.F(pi/2 + im*x, -(b/a)))/(3*sqrt(a + b*cosh(x)^2)) + (1/3)*b*cosh(x)*sqrt(a + b*cosh(x)^2)*sinh(x)]

    @test_int [(1 + cosh(x)^2)^(3/2), x, 4, -2*im*Elliptic.E(pi/2 + im*x, -1) + (2/3)*im*Elliptic.F(pi/2 + im*x, -1) + (1/3)*cosh(x)*sqrt(1 + cosh(x)^2)*sinh(x)]
    @test_int [(1 - cosh(x)^2)^(3/2), x, 4, (2/3)*coth(x)*sqrt(-sinh(x)^2) + (1/3)*coth(x)*(-sinh(x)^2)^(3/2)]
    @test_int [(-1 + cosh(x)^2)^(3/2), x, 4, (-(2/3))*coth(x)*sqrt(sinh(x)^2) + (1/3)*coth(x)*(sinh(x)^2)^(3/2)]
    @test_int [(-1 - cosh(x)^2)^(3/2), x, 6, (2*im*sqrt(-1 - cosh(x)^2)*Elliptic.E(pi/2 + im*x, -1))/sqrt(1 + cosh(x)^2) + (2*im*sqrt(1 + cosh(x)^2)*Elliptic.F(pi/2 + im*x, -1))/(3*sqrt(-1 - cosh(x)^2)) - (1/3)*cosh(x)*sqrt(-1 - cosh(x)^2)*sinh(x)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [1/sqrt(a + b*cosh(x)^2), x, 2, -((im*sqrt(1 + (b*cosh(x)^2)/a)*Elliptic.F(pi/2 + im*x, -(b/a)))/sqrt(a + b*cosh(x)^2))]

    @test_int [1/sqrt(1 + cosh(x)^2), x, 1, (-im)*Elliptic.F(pi/2 + im*x, -1)]
    @test_int [1/sqrt(1 - cosh(x)^2), x, 3, -((atanh(cosh(x))*sinh(x))/sqrt(-sinh(x)^2))]
    @test_int [1/sqrt(-1 + cosh(x)^2), x, 3, -((atanh(cosh(x))*sinh(x))/sqrt(sinh(x)^2))]
    @test_int [1/sqrt(-1 - cosh(x)^2), x, 2, -((im*sqrt(1 + cosh(x)^2)*Elliptic.F(pi/2 + im*x, -1))/sqrt(-1 - cosh(x)^2))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*cosh(c+d*x)^m*(a+b*cosh(c+d*x)^3)^p=#

    @test_int [1/(a + b*cosh(x)^3), x, 8, (2*atanh((sqrt(a^(1/3) - b^(1/3))*tanh(x/2))/sqrt(a^(1/3) + b^(1/3))))/(3*a^(2/3)*sqrt(a^(1/3) - b^(1/3))*sqrt(a^(1/3) + b^(1/3))) + (2*atanh((sqrt(a^(1/3) + (-1)^(1/3)*b^(1/3))*tanh(x/2))/sqrt(a^(1/3) - (-1)^(1/3)*b^(1/3))))/(3*a^(2/3)*sqrt(a^(1/3) - (-1)^(1/3)*b^(1/3))*sqrt(a^(1/3) + (-1)^(1/3)*b^(1/3))) + (2*atanh((sqrt(a^(1/3) - (-1)^(2/3)*b^(1/3))*tanh(x/2))/sqrt(a^(1/3) + (-1)^(2/3)*b^(1/3))))/(3*a^(2/3)*sqrt(a^(1/3) - (-1)^(2/3)*b^(1/3))*sqrt(a^(1/3) + (-1)^(2/3)*b^(1/3)))]
    @test_int [1/(a - b*cosh(x)^3), x, 8, (2*atanh((sqrt(a^(1/3) + b^(1/3))*tanh(x/2))/sqrt(a^(1/3) - b^(1/3))))/(3*a^(2/3)*sqrt(a^(1/3) - b^(1/3))*sqrt(a^(1/3) + b^(1/3))) + (2*atanh((sqrt(a^(1/3) - (-1)^(1/3)*b^(1/3))*tanh(x/2))/sqrt(a^(1/3) + (-1)^(1/3)*b^(1/3))))/(3*a^(2/3)*sqrt(a^(1/3) - (-1)^(1/3)*b^(1/3))*sqrt(a^(1/3) + (-1)^(1/3)*b^(1/3))) + (2*atanh((sqrt(a^(1/3) + (-1)^(2/3)*b^(1/3))*tanh(x/2))/sqrt(a^(1/3) - (-1)^(2/3)*b^(1/3))))/(3*a^(2/3)*sqrt(a^(1/3) - (-1)^(2/3)*b^(1/3))*sqrt(a^(1/3) + (-1)^(2/3)*b^(1/3)))]
    @test_int [1/(1 + cosh(x)^3), x, 7, -((2*(-(1/3))^(1/4)*atan((-1)^(3/4)*3^(1/4)*tanh(x/2)))/(3*(1 - (-1)^(1/3)))) - (2*(-(1/3))^(1/4)*atanh((-1)^(3/4)*3^(1/4)*tanh(x/2)))/(3*(1 + (-1)^(2/3))) + sinh(x)/(3*(1 + cosh(x)))]
    @test_int [1/(1 - cosh(x)^3), x, 7, -((2*(-1)^(1/4)*atan(((-1)^(3/4)*tanh(x/2))/3^(1/4)))/(3^(3/4)*(1 - (-1)^(2/3)))) - (2*(-1)^(1/4)*atanh(((-1)^(3/4)*tanh(x/2))/3^(1/4)))/(3^(3/4)*(1 + (-1)^(1/3))) - sinh(x)/(3*(1 - cosh(x)))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*cosh(c+d*x)^m*(a+b*cosh(c+d*x)^4)^p=#

    @test_int [1/(a + b*cosh(x)^4), x, 10, (sqrt(sqrt(a) - sqrt(a + b))*atanh((sqrt(sqrt(a) + sqrt(a + b)) - sqrt(2)*a^(1/4)*tanh(x))/sqrt(sqrt(a) - sqrt(a + b))))/(2*sqrt(2)*a^(3/4)*sqrt(a + b)) - (sqrt(sqrt(a) - sqrt(a + b))*atanh((sqrt(sqrt(a) + sqrt(a + b)) + sqrt(2)*a^(1/4)*tanh(x))/sqrt(sqrt(a) - sqrt(a + b))))/(2*sqrt(2)*a^(3/4)*sqrt(a + b)) - (sqrt(sqrt(a) + sqrt(a + b))*log(sqrt(a + b) - sqrt(2)*a^(1/4)*sqrt(sqrt(a) + sqrt(a + b))*tanh(x) + sqrt(a)*tanh(x)^2))/(4*sqrt(2)*a^(3/4)*sqrt(a + b)) + (sqrt(sqrt(a) + sqrt(a + b))*log(sqrt(a + b) + sqrt(2)*a^(1/4)*sqrt(sqrt(a) + sqrt(a + b))*tanh(x) + sqrt(a)*tanh(x)^2))/(4*sqrt(2)*a^(3/4)*sqrt(a + b)), ((sqrt(a) - sqrt(a + b))*atan((a^(1/4)*sqrt(a + b + sqrt(a)*sqrt(a + b)) - sqrt(2)*(a + b)^(3/4)*coth(x))/(a^(1/4)*sqrt(a + b - sqrt(a)*sqrt(a + b)))))/(2*sqrt(2)*a^(3/4)*(a + b)^(1/4)*sqrt(a + b - sqrt(a)*sqrt(a + b))) - ((sqrt(a) - sqrt(a + b))*atan((a^(1/4)*sqrt(a + b + sqrt(a)*sqrt(a + b)) + sqrt(2)*(a + b)^(3/4)*coth(x))/(a^(1/4)*sqrt(a + b - sqrt(a)*sqrt(a + b)))))/(2*sqrt(2)*a^(3/4)*(a + b)^(1/4)*sqrt(a + b - sqrt(a)*sqrt(a + b))) - ((sqrt(a) + sqrt(a + b))*log(sqrt(a)*(a + b)^(1/4) - sqrt(2)*a^(1/4)*sqrt(a + b + sqrt(a)*sqrt(a + b))*coth(x) + (a + b)^(3/4)*coth(x)^2))/(4*sqrt(2)*a^(3/4)*(a + b)^(1/4)*sqrt(a + b + sqrt(a)*sqrt(a + b))) + ((sqrt(a) + sqrt(a + b))*log(sqrt(a)*(a + b)^(1/4) + sqrt(2)*a^(1/4)*sqrt(a + b + sqrt(a)*sqrt(a + b))*coth(x) + (a + b)^(3/4)*coth(x)^2))/(4*sqrt(2)*a^(3/4)*(a + b)^(1/4)*sqrt(a + b + sqrt(a)*sqrt(a + b)))]

    @test_int [1/(a - b*cosh(x)^4), x, 4, atanh((a^(1/4)*tanh(x))/sqrt(sqrt(a) - sqrt(b)))/(2*a^(3/4)*sqrt(sqrt(a) - sqrt(b))) + atanh((a^(1/4)*tanh(x))/sqrt(sqrt(a) + sqrt(b)))/(2*a^(3/4)*sqrt(sqrt(a) + sqrt(b)))]

    @test_int [1/(1 + cosh(x)^4), x, 10, -(atan((sqrt(1 + sqrt(2)) - 2*coth(x))/sqrt(-1 + sqrt(2)))/(4*sqrt(1 + sqrt(2)))) + atan((sqrt(1 + sqrt(2)) + 2*coth(x))/sqrt(-1 + sqrt(2)))/(4*sqrt(1 + sqrt(2))) - (1/8)*sqrt(1 + sqrt(2))*log(sqrt(2) - 2*sqrt(1 + sqrt(2))*coth(x) + 2*coth(x)^2) + (1/8)*sqrt(1 + sqrt(2))*log(1 + sqrt(2*(1 + sqrt(2)))*coth(x) + sqrt(2)*coth(x)^2)]

    @test_int [1/(1 - cosh(x)^4), x, 3, atanh(tanh(x)/sqrt(2))/(2*sqrt(2)) + coth(x)/2]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*cosh(c+d*x)^m*(a+b*cosh(c+d*x)^n)^p=#

    @test_int [1/(a + b*cosh(x)^5), x, 12, (2*atanh((sqrt(a^(1/5) - b^(1/5))*tanh(x/2))/sqrt(a^(1/5) + b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - b^(1/5))*sqrt(a^(1/5) + b^(1/5))) + (2*atanh((sqrt(a^(1/5) + (-1)^(1/5)*b^(1/5))*tanh(x/2))/sqrt(a^(1/5) - (-1)^(1/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(1/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(1/5)*b^(1/5))) + (2*atanh((sqrt(a^(1/5) - (-1)^(2/5)*b^(1/5))*tanh(x/2))/sqrt(a^(1/5) + (-1)^(2/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(2/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(2/5)*b^(1/5))) + (2*atanh((sqrt(a^(1/5) + (-1)^(3/5)*b^(1/5))*tanh(x/2))/sqrt(a^(1/5) - (-1)^(3/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(3/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(3/5)*b^(1/5))) + (2*atanh((sqrt(a^(1/5) - (-1)^(4/5)*b^(1/5))*tanh(x/2))/sqrt(a^(1/5) + (-1)^(4/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(4/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(4/5)*b^(1/5)))]
    @test_int [1/(a + b*cosh(x)^6), x, 7, atanh((a^(1/6)*tanh(x))/sqrt(a^(1/3) + b^(1/3)))/(3*a^(5/6)*sqrt(a^(1/3) + b^(1/3))) + atanh((a^(1/6)*tanh(x))/sqrt(a^(1/3) - (-1)^(1/3)*b^(1/3)))/(3*a^(5/6)*sqrt(a^(1/3) - (-1)^(1/3)*b^(1/3))) + atanh((a^(1/6)*tanh(x))/sqrt(a^(1/3) + (-1)^(2/3)*b^(1/3)))/(3*a^(5/6)*sqrt(a^(1/3) + (-1)^(2/3)*b^(1/3)))]
    @test_int [1/(a + b*cosh(x)^8), x, 9, -(atanh(((-a)^(1/8)*tanh(x))/sqrt((-a)^(1/4) - b^(1/4)))/(4*(-a)^(7/8)*sqrt((-a)^(1/4) - b^(1/4)))) - atanh(((-a)^(1/8)*tanh(x))/sqrt((-a)^(1/4) - im*b^(1/4)))/(4*(-a)^(7/8)*sqrt((-a)^(1/4) - im*b^(1/4))) - atanh(((-a)^(1/8)*tanh(x))/sqrt((-a)^(1/4) + im*b^(1/4)))/(4*(-a)^(7/8)*sqrt((-a)^(1/4) + im*b^(1/4))) - atanh(((-a)^(1/8)*tanh(x))/sqrt((-a)^(1/4) + b^(1/4)))/(4*(-a)^(7/8)*sqrt((-a)^(1/4) + b^(1/4)))]

    @test_int [1/(a - b*cosh(x)^5), x, 12, (2*atanh((sqrt(a^(1/5) + b^(1/5))*tanh(x/2))/sqrt(a^(1/5) - b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - b^(1/5))*sqrt(a^(1/5) + b^(1/5))) + (2*atanh((sqrt(a^(1/5) - (-1)^(1/5)*b^(1/5))*tanh(x/2))/sqrt(a^(1/5) + (-1)^(1/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(1/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(1/5)*b^(1/5))) + (2*atanh((sqrt(a^(1/5) + (-1)^(2/5)*b^(1/5))*tanh(x/2))/sqrt(a^(1/5) - (-1)^(2/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(2/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(2/5)*b^(1/5))) + (2*atanh((sqrt(a^(1/5) - (-1)^(3/5)*b^(1/5))*tanh(x/2))/sqrt(a^(1/5) + (-1)^(3/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(3/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(3/5)*b^(1/5))) + (2*atanh((sqrt(a^(1/5) + (-1)^(4/5)*b^(1/5))*tanh(x/2))/sqrt(a^(1/5) - (-1)^(4/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(4/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(4/5)*b^(1/5)))]
    @test_int [1/(a - b*cosh(x)^6), x, 7, atanh((a^(1/6)*tanh(x))/sqrt(a^(1/3) - b^(1/3)))/(3*a^(5/6)*sqrt(a^(1/3) - b^(1/3))) + atanh((a^(1/6)*tanh(x))/sqrt(a^(1/3) + (-1)^(1/3)*b^(1/3)))/(3*a^(5/6)*sqrt(a^(1/3) + (-1)^(1/3)*b^(1/3))) + atanh((a^(1/6)*tanh(x))/sqrt(a^(1/3) - (-1)^(2/3)*b^(1/3)))/(3*a^(5/6)*sqrt(a^(1/3) - (-1)^(2/3)*b^(1/3)))]
    @test_int [1/(a - b*cosh(x)^8), x, 9, atanh((a^(1/8)*tanh(x))/sqrt(a^(1/4) - b^(1/4)))/(4*a^(7/8)*sqrt(a^(1/4) - b^(1/4))) + atanh((a^(1/8)*tanh(x))/sqrt(a^(1/4) - im*b^(1/4)))/(4*a^(7/8)*sqrt(a^(1/4) - im*b^(1/4))) + atanh((a^(1/8)*tanh(x))/sqrt(a^(1/4) + im*b^(1/4)))/(4*a^(7/8)*sqrt(a^(1/4) + im*b^(1/4))) + atanh((a^(1/8)*tanh(x))/sqrt(a^(1/4) + b^(1/4)))/(4*a^(7/8)*sqrt(a^(1/4) + b^(1/4)))]

    @test_int [1/(1 + cosh(x)^5), x, 11, -((2*atan(tanh(x/2)/sqrt(-((1 - (-1)^(1/5))/(1 + (-1)^(1/5))))))/(5*sqrt(-1 + (-1)^(2/5)))) - (2*sqrt(-((1 + (-1)^(3/5))/(1 - (-1)^(3/5))))*atan(sqrt(-((1 + (-1)^(3/5))/(1 - (-1)^(3/5))))*tanh(x/2)))/(5*(1 + (-1)^(3/5))) + (2*atanh(sqrt((1 - (-1)^(2/5))/(1 + (-1)^(2/5)))*tanh(x/2)))/(5*sqrt(1 - (-1)^(4/5))) + (2*atanh(sqrt((1 - (-1)^(4/5))/(1 + (-1)^(4/5)))*tanh(x/2)))/(5*sqrt(1 + (-1)^(3/5))) + sinh(x)/(5*(1 + cosh(x)))]
    @test_int [1/(1 + cosh(x)^6), x, 7, atanh(tanh(x)/sqrt(2))/(3*sqrt(2)) + atanh(tanh(x)/sqrt(1 - (-1)^(1/3)))/(3*sqrt(1 - (-1)^(1/3))) + atanh(tanh(x)/sqrt(1 + (-1)^(2/3)))/(3*sqrt(1 + (-1)^(2/3)))]
    @test_int [1/(1 + cosh(x)^8), x, 9, atanh(tanh(x)/sqrt(1 - (-1)^(1/4)))/(4*sqrt(1 - (-1)^(1/4))) + atanh(tanh(x)/sqrt(1 + (-1)^(1/4)))/(4*sqrt(1 + (-1)^(1/4))) + atanh(tanh(x)/sqrt(1 - (-1)^(3/4)))/(4*sqrt(1 - (-1)^(3/4))) + atanh(tanh(x)/sqrt(1 + (-1)^(3/4)))/(4*sqrt(1 + (-1)^(3/4)))]

    @test_int [1/(1 - cosh(x)^5), x, 11, -((2*atan(tanh(x/2)/sqrt(-((1 - (-1)^(2/5))/(1 + (-1)^(2/5))))))/(5*sqrt(-1 + (-1)^(4/5)))) + (2*atan(sqrt(-((1 + (-1)^(4/5))/(1 - (-1)^(4/5))))*tanh(x/2)))/(5*sqrt(-1 - (-1)^(3/5))) + (2*atanh(sqrt((1 - (-1)^(1/5))/(1 + (-1)^(1/5)))*tanh(x/2)))/(5*sqrt(1 - (-1)^(2/5))) + (2*atanh(sqrt((1 - (-1)^(3/5))/(1 + (-1)^(3/5)))*tanh(x/2)))/(5*sqrt(1 + (-1)^(1/5))) - sinh(x)/(5*(1 - cosh(x)))]
    @test_int [1/(1 - cosh(x)^6), x, 8, atanh(tanh(x)/sqrt(1 + (-1)^(1/3)))/(3*sqrt(1 + (-1)^(1/3))) + atanh(tanh(x)/sqrt(1 - (-1)^(2/3)))/(3*sqrt(1 - (-1)^(2/3))) + coth(x)/3]
    @test_int [1/(1 - cosh(x)^8), x, 10, atanh(tanh(x)/sqrt(1 - im))/(4*sqrt(1 - im)) + atanh(tanh(x)/sqrt(1 + im))/(4*sqrt(1 + im)) + atanh(tanh(x)/sqrt(2))/(4*sqrt(2)) + coth(x)/4]

    #= ::Title:: =#
    #=Integrands*of*the*form*tanh(c+d*x)^m*(a+b*cosh(c+d*x)^n)^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*tanh(c+d*x)^m*(a+b*cosh(c+d*x)^2)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tanh(c+d*x)^m*(a+b*cosh(c+d*x)^2)^p=#

    @test_int [tanh(x)/(1 + cosh(x)^2), x, 4, log(cosh(x)) - (1/2)*log(1 + cosh(x)^2)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tanh(c+d*x)^m*(a+b*cosh(c+d*x)^2)^(p/2)=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [tanh(x)*sqrt(a + b*cosh(x)^2), x, 4, (-sqrt(a))*atanh(sqrt(a + b*cosh(x)^2)/sqrt(a)) + sqrt(a + b*cosh(x)^2)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [tanh(x)/sqrt(a + b*cosh(x)^2), x, 3, -(atanh(sqrt(a + b*cosh(x)^2)/sqrt(a))/sqrt(a))]
    @test_int [tanh(x)/sqrt(1 + cosh(x)^2), x, 3, -atanh(sqrt(1 + cosh(x)^2))]
    @test_int [tanh(x)/sqrt(1 - cosh(x)^2), x, 4, -atanh(sqrt(-sinh(x)^2))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*tanh(c+d*x)^m*(a+b*cosh(c+d*x)^3)^p=#

    @test_int [tanh(x)^3/(a + b*cosh(x)^3), x, 11, -((b^(2/3)*atan((a^(1/3) - 2*b^(1/3)*cosh(x))/(sqrt(3)*a^(1/3))))/(sqrt(3)*a^(5/3))) + log(cosh(x))/a + (b^(2/3)*log(a^(1/3) + b^(1/3)*cosh(x)))/(3*a^(5/3)) - (b^(2/3)*log(a^(2/3) - a^(1/3)*b^(1/3)*cosh(x) + b^(2/3)*cosh(x)^2))/(6*a^(5/3)) - log(a + b*cosh(x)^3)/(3*a) + sech(x)^2/(2*a)]

    @test_int [tanh(x)/sqrt(a + b*cosh(x)^3), x, 4, -((2*atanh(sqrt(a + b*cosh(x)^3)/sqrt(a)))/(3*sqrt(a)))]
    @test_int [tanh(x)*sqrt(a + b*cosh(x)^3), x, 5, (-(2/3))*sqrt(a)*atanh(sqrt(a + b*cosh(x)^3)/sqrt(a)) + (2/3)*sqrt(a + b*cosh(x)^3)]

    #= ::Section:: =#
    #=Integrands*of*the*form*tanh(c+d*x)^m*(a+b*cosh(c+d*x)^4)^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*tanh(c+d*x)^m*(a+b*cosh(c+d*x)^n)^p=#

    @test_int [tanh(x)/sqrt(a + b*cosh(x)^n), x, 4, -((2*atanh(sqrt(a + b*cosh(x)^n)/sqrt(a)))/(sqrt(a)*n))]

    @test_int [tanh(x)*sqrt(a + b*cosh(x)^n), x, 5, -((2*sqrt(a)*atanh(sqrt(a + b*cosh(x)^n)/sqrt(a)))/n) + (2*sqrt(a + b*cosh(x)^n))/n]
    end

max_stars_repo_path: test/6 Hyperbolic functions/6.2 Hyperbolic cosine/6.2.7 hyper^m (a+b cosh^n)^p.jl


text:
"""
    get_comments(target::String; concatenate::Bool=false) -> Vector{String}

Extract the comments from a passed PDF or all PDFs found recursively in the
passed directory.

# Arguments
- `target::String`: $(TARGET_PDF_ARGUMENT)

# Keywords
- `concatenate::Bool=false`: $(CONCATENATE_KEYWORD)

# Returns
- `Vector{String}`: the comments

# Throws
- Exceptions from: [`get_highlights_comments_pages`](@ref)

# Example
```jldoctest; output = false
using PDFHighlights

path_to_pdf_dir = joinpath(pathof(PDFHighlights) |> dirname |> dirname, "test", "pdf")
path_to_pdf = joinpath(path_to_pdf_dir, "TestPDF.pdf")

get_comments(path_to_pdf_dir; concatenate=true) ==
get_comments(path_to_pdf; concatenate=true) ==
["Comment 1", "Comment 2 Comment 3", "Comment 4", "", "", "", ""]

# output

true
```
"""
function get_comments(target::String; concatenate::Bool=false)::Vector{String}
    _, comments, _ = get_highlights_comments_pages(target; concatenate)
    return comments
end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module AbstractCollisionAvoidanceSystemImpl

export AbstractCollisionAvoidanceSystem

abstract AbstractCollisionAvoidanceSystem

end

module AbstractCollisionAvoidanceSystemInterfaces

export testThreat,
       selectRA

function testThreat()
end

function selectRA()
end

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
##############################################################################
##
## AbstractFixedEffectSolver
##
## this type must defined solve_residuals!, solve_coefficients!
##
##############################################################################
abstract type AbstractFixedEffectSolver{T} end

"""
`solve_residuals!(y, fes, w; method = :cpu, double_precision = true, tol = 1e-8, maxiter = 10000, )`

Returns ``y_i - X_i'\\beta`` where ``\\beta = argmin_{b} \\sum_i y_i - X_i'b``, where `X` denotes the matrix of fixed effects `fes`.

### Arguments
* `y` : A `AbstractVector` or A `AbstractMatrix`
* `fes`: A `Vector{<:FixedEffect}`
* `w`: A vector of weights, i.e. `AbstractWeights`
* `method` : A `Symbol` for the method. Default is :cpu. The option :gpu requires `using CUDA` (in this case, it is recommanded to use the option `double_precision = false`).
* `double_precision::Bool`: Should the demeaning operation use Float64 rather than Float32? Default to true.
* `tol` : Tolerance. Default to 1e-8 if `double_precision = true`, 1e-6 otherwise.
* `maxiter` : Maximum number of iterations

### Returns
* `res` :  Residual of the least square problem
* `iterations`: Number of iterations
* `converged`: Did the algorithm converge?

### Examples
```julia
using  FixedEffects
p1 = repeat(1:5, inner = 2)
p2 = repeat(1:5, outer = 2)
solve_residuals!(rand(10), [FixedEffect(p1), FixedEffect(p2)])
```
"""
function solve_residuals!(y::Union{AbstractVector{<: Number}, AbstractMatrix{<: Number}}, fes::AbstractVector{<: FixedEffect}, w::AbstractWeights = uweights(eltype(y), size(y, 1));
	method::Symbol = :cpu, double_precision::Bool = eltype(y) == Float64,
	tol::Real = double_precision ? 1e-8 : 1e-6, maxiter::Integer = 10000,
	nthreads = method == :cpu ? Threads.nthreads() : 256)
	any((length(fe) != size(y, 1) for fe in fes)) && throw("FixedEffects must have the same length as y")
	any(ismissing.(fes)) && throw("FixedEffects must not have missing values")
	feM = AbstractFixedEffectSolver{double_precision ? Float64 : Float32}(fes, w, Val{method}, nthreads)
	solve_residuals!(y, feM; maxiter = maxiter, tol = tol)
end

"""
Solve a least square problem for a set of FixedEffects

`solve_coefficients!(y, fes, w; method = :cpu, maxiter = 10000, double_precision = true, tol = 1e-8)`

Returns ``\\beta = argmin_{b} \\sum_i w_i(y_i - X_i'b)`` where `X` denotes the matrix of fixed effects `fes`.

### Arguments
* `y` : A `AbstractVector`
* `fes`: A `Vector{<:FixedEffect}`
* `w`: A vector of weights, i.e. `AbstractWeights`
* `method` : A `Symbol` for the method. Default is :cpu. The option :gpu requires `using CUDA` (in this case, it is recommanded to use the option `double_precision = false`).
* `double_precision::Bool`: Should the demeaning operation use Float64 rather than Float32? Default to true.
* `tol` : Tolerance. Default to 1e-8 if `double_precision = true`, 1e-6 otherwise.
* `maxiter` : Maximum number of iterations
* `nthreads` : Number of threads

### Returns
* ``\\beta`` : Solution of the least square problem
* `iterations`: Number of iterations
* `converged`: Did the algorithm converge?
Fixed effects are generally not unique. We standardize the solution
in the following way: the mean of fixed effects within connected components is zero
(except for the first).
This gives the unique solution in the case of two fixed effects.

### Examples
```julia
using  FixedEffects
p1 = repeat(1:5, inner = 2)
p2 = repeat(1:5, outer = 2)
x = rand(10)
solve_coefficients!(rand(10), [FixedEffect(p1), FixedEffect(p2)])
```
"""
function solve_coefficients!(y::AbstractVector{<: Number}, fes::AbstractVector{<: FixedEffect}, w::AbstractWeights = uweights(eltype(y), length(y)); method::Symbol = :cpu, double_precision::Bool = eltype(y) == Float64,
	tol::Real = double_precision ? 1e-8 : 1e-6,  maxiter::Integer = 10000,
	nthreads = method == :cpu ? Threads.nthreads() : 256)
	any(ismissing.(fes)) && throw("Some FixedEffect has a missing value for reference or interaction")
	any((length(fe) != length(y) for fe in fes))  && throw("FixedEffects must have the same length as y")
	feM = AbstractFixedEffectSolver{double_precision ? Float64 : Float32}(fes, w, Val{method}, nthreads)
	solve_coefficients!(y, feM; maxiter = maxiter, tol = tol)
end

function AbstractFixedEffectSolver{T}(fes::Vector{<:FixedEffect}, weights::AbstractWeights, ::Union{Type{Val{:lsmr}}, Type{Val{:lsmr_threads}}, Type{Val{:lsmr_cores}}}) where {T}
	@warn ":lsmr, :lsmr_threads, and :lsmr_cores are deprecated (FixedEffects is now always multi-threaded)"
	AbstractFixedEffectSolver{T}(fes, weights, Val{:cpu})
end

function AbstractFixedEffectSolver{T}(fes::Vector{<:FixedEffect}, weights::AbstractWeights, ::Type{Val{:lsmr_gpu}}) where {T}
	@warn ":lsmr_gpu is deprecated,  use :gpu"
	AbstractFixedEffectSolver{T}(fes, weights, Val{:gpu})
end

max_issues_repo_name: matthieugomez/FixedEffects

max_forks_repo_name: matthieugomez/FixedEffects