text:
module consoles

export Console, console, err_console, console_height, console_width

"""
    Console

The `Console` object stores information about the dimensions of the output(::IO)
where objects will be printed
"""
struct Console
    io::IO
    width::Int
    height::Int
end

Console(io::IO) = Console(io, displaysize(io)[2], displaysize(io)[1])
Console() = Console(stdout)

console = Console(stdout)
err_console = Console(stderr)

"""
    console_height()

Get the current console height.
"""
console_height() = displaysize(stdout)[1]
console_height(io::IO) = displaysize(io)[1]

"""
    console_width()

Get the current console width.
"""
console_width() = displaysize(stdout)[2]
console_width(io::IO) = displaysize(io)[2]

end


text:
#Mapping functions
##################
# function index2d(k,nx)
#     j = floor(Int,(k-1)/nx)+1
#     # TODO: add consistency check with ny?
#     return k-(j-1)*nx,j
# end

function index1d(i,j,nx)
    # TODO: add consistency check with ny?
    return i+(j-1)*nx
end

"""
    compute_equilibrium!(...)

* rel_error is the tolerance for global temperature equilibrium (default is 2e-5).
* max_years is the maximum number of annual cycles to be computed when searching for equilibrium
"""
function compute_equilibrium!(discretization; max_years=100, rel_error=2e-5, verbose=true)
    @unpack mesh, model, lu_decomposition, num_steps_year, annual_temperature, rhs, last_rhs  = discretization
    @unpack nx, dof = mesh

    average_temperature = average_temperature_old = area_weighted_average(view(annual_temperature, :, num_steps_year), mesh)

    if verbose
      println("year","  ","Average Temperature")
      println(0,"  ",average_temperature_old)
    end

    for year in 1:max_years
        average_temperature = 0.0
        for time_step in 1:num_steps_year
            old_time_step = (time_step == 1) ? num_steps_year : time_step - 1
            update_rhs!(rhs, mesh, num_steps_year, time_step, view(annual_temperature, :, old_time_step), model, last_rhs)

            # Use in-place operation `ldiv!` instead of `\` to avoid allocations
            ldiv!(view(annual_temperature, :, time_step), lu_decomposition, rhs)

            annual_temperature[1:nx, time_step] .= annual_temperature[1, time_step]
            annual_temperature[dof-nx+1:dof, time_step] .= annual_temperature[dof, time_step]

            average_temperature += area_weighted_average(view(annual_temperature, :, time_step), mesh)
        end
        average_temperature = average_temperature/num_steps_year
        if verbose
          println(year,"  ",average_temperature)
        end

        if (abs(average_temperature-average_temperature_old)<rel_error)
            if verbose
              println("EQUILIBRIUM REACHED!")
            end
            break
        end

        average_temperature_old = average_temperature
    end

    return average_temperature
end

"""
    compute_evolution!(...)

Compute the evolution of the mean temperature with varying CO2 levels.
"""
function compute_evolution!(discretization, co2_concentration_at_step, year_start, year_end; verbose=true)
    @unpack mesh, model, lu_decomposition, num_steps_year, annual_temperature, rhs, last_rhs  = discretization
    @unpack nx, dof = mesh

    if verbose
      println("year","  ","Average Temperature")
    end

    max_years = year_end - year_start + 1

    # Allocate output arrays
    year_at_step = zeros(Float64,max_years*num_steps_year+1)
    mean_temperature_at_step = zeros(Float64,max_years*num_steps_year+1)
    year_array = zeros(Float64,max_years)
    mean_temperature_yearly = zeros(Float64,max_years)

    # Fill year arrays
    year_at_step[1] = year_start + 0.21644  # The simulation starts at the vernal equinox
    for year in 2:size(year_at_step,1)
        year_at_step[year] = year_at_step[year-1] + 1.0/num_steps_year
    end
    year_array[1] = year_start + 0.71644    # The simulation starts at the vernal equinox
    for year in 2:max_years
        year_array[year] = year_array[year-1] + 1.0
    end

    # Initialize the temperature array
    mean_temperature_at_step[1] = area_weighted_average(view(annual_temperature, :, num_steps_year), mesh)
    step = 1

    for year in 1:max_years
        average_temperature = 0.0
        for time_step in 1:num_steps_year
            set_co2_concentration!(model, co2_concentration_at_step[step])
            old_time_step = (time_step == 1) ? num_steps_year : time_step - 1
            update_rhs!(rhs, mesh, num_steps_year, time_step, view(annual_temperature, :, old_time_step), model, last_rhs)

            # Use in-place operation `ldiv!` instead of `\` to avoid allocations
            ldiv!(view(annual_temperature, :, time_step), lu_decomposition, rhs)

            annual_temperature[1:nx, time_step] .= annual_temperature[1, time_step]
            annual_temperature[dof-nx+1:dof, time_step] .= annual_temperature[dof, time_step]

            # Compute mean temperature
            step += 1
            mean_temperature_at_step[step] = area_weighted_average(view(annual_temperature, :, time_step), mesh)
            average_temperature += mean_temperature_at_step[step]
        end
        average_temperature = average_temperature/num_steps_year
        mean_temperature_yearly[year] = average_temperature
        if verbose
          println(year_start+year-1,"  ",average_temperature)
        end
    end

    return (; year_at_step, mean_temperature_at_step, year_array, mean_temperature_yearly)
end

"""
Computes mean temperature in the globe at a specific time
"""
function area_weighted_average(vector,mesh)
    @unpack nx,ny,area,dof = mesh

    average_vector = 0.0

    # Contribution of inner points
    for j=2:ny-1
        for i=1:nx
            row_idx = index1d(i,j,nx)
            average_vector += area[j] * vector[row_idx]
        end
    end

    # Poles
    average_vector += area[1] * (vector[1] + vector[dof])

    return average_vector
end

function compute_matrix(mesh,num_steps_year,model)
    @unpack nx,ny,dof,h,geom,csc2,cot,area = mesh
    @unpack diffusion_coeff,heat_capacity,radiative_cooling_feedback = model
    matrix = spzeros(Float64,dof,dof)
    sh2 = 1 / h^2
    # Inner DOFs (c coefficients are divided by h^2)
    for j=2:ny-1
        for i=1:nx

            # Compute coefficients
            c0 = 2 * sh2 * diffusion_coeff[i,j] * (1 + csc2[j]) + 2 * heat_capacity[i,j] * num_steps_year + radiative_cooling_feedback

            # Get diffusion coefficient
            if (i == 1) # Periodic BC
                d_phi = (diffusion_coeff[2,j] - diffusion_coeff[nx,j]) / (2 * h)
            elseif (i == nx) # Periodic BC
                d_phi = (diffusion_coeff[1,j] - diffusion_coeff[nx-1,j]) / (2 * h)
            else # Inner DOFs
                d_phi = (diffusion_coeff[i+1,j] - diffusion_coeff[i-1,j]) / (2 * h)
            end

            c1 = sh2 * csc2[j] * (diffusion_coeff[i,j] - 0.5 * h * d_phi)
            c3 = sh2 * csc2[j] * (diffusion_coeff[i,j] + 0.5 * h * d_phi)

            d_theta= (diffusion_coeff[i,j+1] - diffusion_coeff[i,j-1]) / (2 * h)
            c2 = sh2 * (diffusion_coeff[i,j] - 0.5 * h * (diffusion_coeff[i,j] * cot[j] + d_theta))
            c4 = sh2 * (diffusion_coeff[i,j] + 0.5 * h * (diffusion_coeff[i,j] * cot[j] + d_theta))

            # Fill matrix A
            row_idx = index1d(i,j,nx)
            col_idx_c0 = row_idx

            if (i==1) # Periodic BC
                col_idx_c1 = index1d(nx,j,nx)
                col_idx_c3 = row_idx + 1
            elseif (i==nx) # Periodic BC
                col_idx_c1 = row_idx - 1
                col_idx_c3 = index1d(1,j,nx)
            else
                col_idx_c1 = row_idx - 1
                col_idx_c3 = row_idx + 1
            end
            col_idx_c2 = row_idx - nx
            col_idx_c4 = row_idx + nx

            matrix[row_idx,col_idx_c0] = c0
            matrix[row_idx,col_idx_c1] = -c1
            matrix[row_idx,col_idx_c2] = -c2
            matrix[row_idx,col_idx_c3] = -c3
            matrix[row_idx,col_idx_c4] = -c4

            # Get rhs
            # rhs[row_idx] = -(c0 * Temp[col_idx_c0] - c1 * Temp[col_idx_c0] - c2 * Temp[col_idx_c2] - c3 * Temp[col_idx_c3] - c4 * Temp[col_idx_c4]) + 4 * heat_capacity[i,j]
        end
    end

    # Poles
    total_area = area[1] + area[2]

    diff_total_np = zeros(Float64,nx)
    diff_total_sp = zeros(Float64,nx)
    for i=1:nx
        diff_total_np[i] = (area[1] * diffusion_coeff[1,1]  + area[2] * diffusion_coeff[i,2]) / total_area
        diff_total_sp[i] = (area[1] * diffusion_coeff[1,ny] + area[2] * diffusion_coeff[i,ny-1]) / total_area
    end
    diff_total_sum_np = sum(diff_total_np)
    diff_total_sum_sp = sum(diff_total_sp)

    gc_np = geom * diff_total_sum_np + radiative_cooling_feedback + 2 * heat_capacity[1,1] * num_steps_year     # north pole
    gc_sp = geom * diff_total_sum_sp + radiative_cooling_feedback + 2 * heat_capacity[1,ny] * num_steps_year    # south pole

    for i=1:nx
        #North pole
        row_idx = i
        matrix[row_idx,row_idx] = gc_np
        for ii=1:nx
            matrix[row_idx,nx+ii] = - geom * diff_total_np[ii]
        end

        #South pole
        row_idx = dof + 1 - i
        matrix[row_idx,row_idx] = gc_sp
        for ii=1:nx
            matrix[row_idx,dof-2*nx+ii] = - geom * diff_total_sp[ii]
        end
    end
    return matrix
end

function update_rhs!(rhs, mesh, num_steps_year, time_step, temperature, model,  last_rhs)
    @unpack nx,ny = mesh
    @unpack heat_capacity, solar_forcing, radiative_cooling_co2 = model
    for j=1:ny
        for i=1:nx
            row_idx = index1d(i,j,nx)

            rhs[row_idx] = 4 * heat_capacity[i,j] * temperature[row_idx] * num_steps_year  - last_rhs[row_idx] + solar_forcing[i,j,time_step] - 2 * radiative_cooling_co2
            if (time_step == 1)
                rhs[row_idx] += solar_forcing[i,j,num_steps_year]
            else
                rhs[row_idx] += solar_forcing[i,j,time_step-1]
            end
        end
    end
    last_rhs .= rhs
end

ext: jl


text:
## Interface for planning heuristics ##
export Heuristic, precompute!, compute
export use_heuristic_cache!, clear_heuristic_cache!

"Global flag as to whether heuristic cache is enabled."
const _use_heuristic_cache = Ref(false)
"Globally enable or disable caching of heuristic values."
use_heuristic_cache!(val::Bool=true) = _use_heuristic_cache[] = val

"Cached heuristic values."
const _heuristic_cache = Dict{Tuple{UInt,Symbol,UInt,UInt}, Real}()
"Clear cache of heuristic values."
clear_heuristic_cache!() = empty!(_heuristic_cache)

"Abstract heuristic type, which defines the interface for planning heuristics."
abstract type Heuristic end

"Precomputes heuristic information given a domain, state, and goal."
precompute!(h::Heuristic, domain::Domain, state::State, spec::Specification) =
    h # Return the heuristic unmodified by default
precompute!(h::Heuristic, domain::Domain, state::State, spec) =
    precompute!(h, domain, state, Specification(spec))
precompute!(h::Heuristic, domain::Domain, state::State) =
    precompute!(h, domain, state, NullGoal())
precompute!(h::Heuristic, domain::Domain) =
    precompute!(h, domain, State(Term[]), NullGoal())

"Returns whether heuristic has been precomputed for a domain, state, and goal."
is_precomputed(h::Heuristic, domain::Domain, state::State, spec::Specification) =
    false

"Computes the heuristic value of state relative to a goal in a given domain."
compute(h::Heuristic, domain::Domain, state::State, spec::Specification) =
    error("Not implemented.")
compute(h::Heuristic, domain::Domain, state::State, spec) =
    compute(h, domain, state, Specification(spec))

"Computes the heuristic value of state relative to a goal in a given domain."
function (h::Heuristic)(domain::Domain, state::State, spec::Specification;
                        cache::Bool=_use_heuristic_cache[])
    if (cache)
        key = (hash(h), domain.name, hash(state), hash(spec))
        if haskey(_heuristic_cache, key) return _heuristic_cache[key] end
    end
    val = compute(h, domain, state, spec)
    if (cache) _heuristic_cache[key] = val end
    return val
end

(h::Heuristic)(domain::Domain, state::State, spec;
               cache::Bool=_use_heuristic_cache[]) =
    h(domain, state, Specification(spec); cache=cache)

include("heuristics/basic.jl")
include("heuristics/hsp.jl")
include("heuristics/ff.jl")


text:
using Test, Fracciones

@testset "Constructores" begin
    @test Fraccion(true, 0x0f) === Fraccion{UInt8}(1, 15) # comprobar tipos
    @test Fraccion(0x00, 5) === Fraccion(0, 1) # ok aunque 0x00 == typemin(UInt8)
    @test Fraccion(2, 0) === Fraccion(1, 0)
    @test Fraccion(-2, 0) === Fraccion(-1, 0)
    # inválidos
    @test_throws ArgumentError Fraccion(0, 0)
    @test_throws ArgumentError Fraccion(Int8(-128), Int8(5))
    @test_throws ArgumentError Fraccion(Int8(5), Int8(-128))
    # otros constructores
    @test Fraccion{Int}(1.0, 2.0) === Fraccion(1,2)
    @test Fraccion{UInt8}(1, 2) === Fraccion(0x01, 0x02)
    @test Fraccion{UInt8}(Fraccion(1,2)) === Fraccion(0x01, 0x02)
    @test Fraccion(0) === Fraccion(0,1)
    @test Fraccion(2) === Fraccion(2,1)
    @test Fraccion{UInt8}(2.0) === Fraccion(0x02,0x01)
    @test Fraccion(2.0) === Fraccion(2,1)
    @test Fraccion(Inf) === Fraccion(1,0)
    @test Fraccion(-Inf) === Fraccion(-1,0)
    @test Fraccion(1,2) === Fraccion(Fraccion(1,2))
end

@testset "Conversiones" begin
    @test Int(Fraccion(4,2)) === 2
    @test float(Fraccion(1,4)) == 0.25
    @test Float64[Fraccion(1,2), Fraccion(1,4)] == [0.5, 0.25]
end

@testset "Signos" begin
    @test Fraccion(5, -2) == Fraccion(-5, 2) == -Fraccion(5,2)
    @test numerador(Fraccion(5,-2)) < 0
    @test denominador(Fraccion(5,-2)) > 0
end

@testset "Propiedades y generadores" begin
    @test numerador(Fraccion(2,4)) == 1
    @test denominador(Fraccion(2,4)) == 2
    @test numerador(Fraccion(3,0)) == 1
    @test reciproco(Fraccion(2, 6)) == 3
    @test reciproco(3) == Fraccion(1, 3)
    @test sign(Fraccion(1,2)) == 1
    @test sign(Fraccion(-1,2)) == -1
    @test sign(Fraccion(0)) == 0
    @test abs(Fraccion(-1,2)) == Fraccion(1,2)
    @test one(Fraccion{UInt8}) === Fraccion(0x01)
    @test one(Fraccion) === Fraccion(1)
    @test zero(Fraccion{UInt8}) === Fraccion(0x00)
    @test zero(Fraccion) === Fraccion(0)
    @test zeros(Fraccion, 2) == Fraccion[0, 0]
    @test ones(Fraccion, 2) == Fraccion.([1, 1])
    @test typemin(Fraccion{Int}) === Fraccion(-1,0)
    @test typemin(Fraccion{UInt8}) === Fraccion(0x00)
    @test typemin(Fraccion{Bool}) === Fraccion(false)
    @test typemax(Fraccion{Int}) === Fraccion(1,0)
    @test typemax(Fraccion{UInt8}) === Fraccion(0x01,0x00)
    @test typemax(Fraccion{Bool}) === Fraccion(true,false)
end

@testset "Operaciones algebraicas" begin
    @test Fraccion(1, 12) + Fraccion(5, 3) == Fraccion(7, 4)
    @test Fraccion(2, 3) - Fraccion(1, 5) == Fraccion(7, 15)
    # con infinitos
    @test Fraccion(1, 0) + Fraccion(1,2) === Fraccion(1, 0)
    @test Fraccion(1, 0) + 3 == Fraccion(1, 0)
    @test Fraccion(1, 3) * Fraccion(3, 4) == Fraccion(1, 4)
    @test Fraccion(1, 2) / Fraccion(3, 4) == Fraccion(2, 3)
    @test_throws ArgumentError Fraccion(0, 1) * Fraccion(-1, 0)

    @test Fraccion(2,3)^-2 == Fraccion(9,4)
    @test Fraccion(2,3)^3 == Fraccion(8,27)
    @test Fraccion(1,3)^0 == 1
    # combinar con otros tipos de números
    @test Fraccion(1,3) + 2 == 2 + Fraccion(1,3) == Fraccion(7,3)
    @test Fraccion(1,3) + 0 == Fraccion(1,3)
    @test 2 - Fraccion(1,0) == Fraccion(-1,0)
    @test Fraccion(1,3) * 2 == 2 * Fraccion(1,3) == Fraccion(2,3)
    @test Fraccion(1,3) * 0 == Fraccion(0,1)
    @test 2 * Fraccion(1,0) == Fraccion(1,0)
    @test Fraccion(1,3) / 2 == Fraccion(1,6)
    @test 2 / Fraccion(1,3) == Fraccion(6)
    @test Fraccion(1,6) / 0 == Fraccion(1,0)
    @test 0 / Fraccion(1,0) == Fraccion(0,1)
end

@testset "Comparaciones" begin
    for (menor, mayor) = (
        (Fraccion(1,4), Fraccion(1,2)), (0.25, Fraccion(1,2)), (Fraccion(1,4), 0.5),
        (Fraccion(0), Fraccion(1,2)), (0, Fraccion(1,2)), (Fraccion(0), 0.5),
        (Fraccion(-1,4), Fraccion(0)), (-0.25, Fraccion(0)), (Fraccion(-1,4), 0),
        (Fraccion(0), Fraccion(1,0)), (0, Fraccion(1,0)), (Fraccion(0), Inf),
        (Fraccion(-1,0), Fraccion(0)), (-Inf, Fraccion(0)), (Fraccion(-1,0), 0),
        (Fraccion(-1,0), Fraccion(1,0)), (-Inf, Fraccion(1,0)), (Fraccion(-1,0), Inf)
    )
        @test menor < mayor
        @test -menor > -mayor
        @test menor ≤ mayor
        @test mayor ≥ menor
        @test menor == menor
        @test mayor == mayor
    end
end

@testset "Función y macro `fraccion`" begin
    @test fraccion(3, Fraccion(5,2)) == Fraccion(6,5)
    @test fraccion(Fraccion(5,2), 3) == Fraccion(5,6)
    @test (@fraccion 3/4) == Fraccion(3, 4)
    @test (@fraccion (3 + (4/3))/5^2) == Fraccion(13, 75) == 13/75
end


text:
################################################################################
#
# Weighting functions for PDE-misfit least-squares norm
#
################################################################################

export src_focusing_fun, src_focusing_mat, src_focusing_array, depth_weighting, depth_weighting_array, compute_srcfcsδ

function src_focusing_fun(q::Array{T, 2}, n::NTuple{2, Int64}, nabc::NTuple{4, Int64}, h::NTuple{2, R}, x_src::Array{R, 1}, z_src::R; δ::R = R(0.001)) where T <: RuC

    wq = Array{T}(undef, size(q))
    Ns = length(x_src)
    next = sz_abc(n, nabc)
    x = reshape(collect((0:next[1]-1)*h[1]).-nabc[1]*h[1], :, 1)
    z = reshape(collect((0:next[2]-1)*h[2]).-nabc[3]*h[2], 1, :)
    for idx_s = 1:Ns
        winv2 = ((x.-x_src[idx_s]).^R(2).+(z.-z_src).^R(2).+δ.^R(2))./δ.^R(2)
        wq[:, idx_s] = q[:, idx_s]./vec(winv2)
    end
    return wq

end

function src_focusing_mat(n::NTuple{2, Int64}, nabc::NTuple{4, Int64}, h::NTuple{2, R}, x_src::R, z_src::R; δ::R = R(0.001))

    next = sz_abc(n, nabc)
    x = reshape(collect((0:next[1]-1)*h[1]).-nabc[1]*h[1], :, 1)
    z = reshape(collect((0:next[2]-1)*h[2]).-nabc[3]*h[2], 1, :)
    return spdiagm(0 => vec(sqrt.(((x.-x_src).^R(2).+(z.-z_src).^R(2).+δ.^R(2))./δ.^R(2))))

end

function src_focusing_array(n::NTuple{2, Int64}, nabc::NTuple{4, Int64}, h::NTuple{2, R}, x_src::Array{R, 1}, z_src::R; δ::R = R(0.001))

    Ns = length(x_src)
    next = sz_abc(n, nabc)
    w = Array{R}(undef, (prod(next), Ns))
    x = reshape(collect((0:next[1]-1)*h[1]).-nabc[1]*h[1], :, 1)
    z = reshape(collect((0:next[2]-1)*h[2]).-nabc[3]*h[2], 1, :)
    for idx_s = 1:Ns
        w[:, idx_s] = vec(((x.-x_src[idx_s]).^R(2).+(z.-z_src).^R(2).+δ.^R(2))./δ.^R(2))
    end
    return w

end

function depth_weighting(q::Array{T, 2}, n::NTuple{2, Int64}, nabc::NTuple{4, Int64}, h::NTuple{2, R}, z_src::R; δ::R = 1f-4) where T <: RuC

    next = sz_abc(n, nabc)
    z = reshape(collect((0:next[2]-1)*h[2]).-nabc[3]*h[2], 1, :)
    winv2 = ((z.-z_src).^R(2).+δ.^R(2))./δ.^R(2)
    return reshape(reshape(q, (next..., size(q, 2)))./winv2, size(q))

end

function depth_weighting_array(n::NTuple{2, Int64}, nabc::NTuple{4, Int64}, h::NTuple{2, R}, z_src::R; δ::R = 1f-4)

    next = sz_abc(n, nabc)
    z = reshape(collect((0:next[2]-1)*h[2]).-nabc[3]*h[2], 1, :)
    return reshape(repeat(sqrt.(((z.-z_src).^2f0.+δ.^2f0)./δ.^2f0), outer = (next[1], 1)), :, 1)

end

function compute_srcfcsδ(m::R, freq::R)
    # Promote focusing on a "basin of attraction" corresponding to one wavelength in radius

    return sqrt(R(3)./m)./freq

end


text:
# LogLinear Transformation

"""
    LogLinearFormula(df::Int64)

Returns GLM formula can be used in `glm` with `FreqTab` data frame.
`df`, is abbreviation for degrees of formula, represents degree of polynomial log-linear method.
This function works very slowly. If a fixed degree formula will be used repeatedly, define the formula as a variable.
"""
function LogLinearFormula(df)
    fml = "@formula(freq ~ 1 +"
    for d in 1:df
        fml *= "scale^$d"
        if d != df
            fml *= " + "
        else
            fml *= ")"
        end
    end
    return eval(Meta.parse(fml))
end

const fml₁ = @formula(freq ~ 1 + scale)
const fml₂ = @formula(freq ~ 1 + scale + scale^2)
const fml₃ = @formula(freq ~ 1 + scale + scale^2 + scale^3)
const fml₄ = @formula(freq ~ 1 + scale + scale^2 + scale^3 + scale^4)
const fml₅ = @formula(freq ~ 1 + scale + scale^2 + scale^3 + scale^4 + scale^5)
const fml₆ = @formula(freq ~ 1 + scale + scale^2 + scale^3 + scale^4 + scale^5 + scale^6)
const fml₇ = @formula(freq ~ 1 + scale + scale^2 + scale^3 + scale^4 + scale^5 + scale^6 + scale^7)
const fml₈ = @formula(freq ~ 1 + scale + scale^2 + scale^3 + scale^4 + scale^5 + scale^6 + scale^7 + scale^8)
const fml₉ = @formula(freq ~ 1 + scale + scale^2 + scale^3 + scale^4 + scale^5 + scale^6 + scale^7 + scale^8 + scale^9)
const fml₁₀ = @formula(freq ~ 1 + scale + scale^2 + scale^3 + scale^4 + scale^5 + scale^6 + scale^7 + scale^8 + scale^9 + scale^10)

mutable struct SmoothedFreqTab <: EG
    table
    raw
    interval
    fit
end
"""
    presmoothing(F::FreqTab, @LogLinearFormula(df::Int64))

Returns presmoothed frequency table as `SmoothedFreqTab` and `glm` fitted object.

Preserving first C moments of original frequency data, passed `LogLinearFormula(C)` to `fml`.
C is a degree of freedom

# Examples
```julia
julia> tab = freqtab(expandtable(ACTmath.scale, ACTmath.xcount));
julia> presmoothing(tab, Equate.fml₄)
Equate.SmoothedFreqTab(40×5 DataFrame
│ Row │ scale   │ freq     │ cumfreq  │ prob        │ cumprob     │
│     │ Float64 │ Float64? │ Float64? │ Float64     │ Float64     │
├─────┼─────────┼──────────┼──────────┼─────────────┼─────────────┤
│ 1   │ 1.0     │ 0.864929 │ 0.864929 │ 0.000199799 │ 0.000199799 │
│ 2   │ 2.0     │ 2.47529  │ 3.34022  │ 0.000571794 │ 0.000771592 │
⋮
│ 38  │ 38.0    │ 31.9351  │ 4291.07  │ 0.00737703  │ 0.991238    │
│ 39  │ 39.0    │ 22.7842  │ 4313.85  │ 0.00526315  │ 0.996501    │
│ 40  │ 40.0    │ 15.1484  │ 4329.0   │ 0.00349928  │ 1.0         │, [1, 2, 3, 3, 3, 4, 4, 4, 4, 4  …  40, 40, 40, 40, 40, 40, 40, 40, 40, 40], 1.0, StatsModels.TableRegressionModel{GLM.GeneralizedLinearModel{GLM.GlmResp{Array{Float64,1},Distributions.Poisson{Float64},GLM.LogLink},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

freq ~ 1 + scale + :(scale ^ 2) + :(scale ^ 3) + :(scale ^ 4)

Coefficients:
─────────────────────────────────────────────────────────────────────────────────
                   Coef.   Std. Error       z  Pr(>|z|)    Lower 95%    Upper 95%
─────────────────────────────────────────────────────────────────────────────────
(Intercept)  -1.35982     0.310243      -4.38    <1e-4   -1.96789     -0.751757
scale         1.30127     0.0728808     17.85    <1e-70   1.15843      1.44412
scale ^ 2    -0.089081    0.00589528   -15.11    <1e-50  -0.100636    -0.0775265
scale ^ 3     0.00254824  0.000195003   13.07    <1e-38   0.00216604   0.00293044
scale ^ 4    -2.67698e-5  2.25114e-6   -11.89    <1e-31  -3.1182e-5   -2.23576e-5
─────────────────────────────────────────────────────────────────────────────────)

```
"""
function presmoothing(F::EG, fml)
    fit1 = glm(fml, F.table, Poisson(), LogLink())
    freq = predict(fit1, DataFrame(scale = F.table.scale))
    tab = DataFrame(scale = F.table.scale, freq = freq, cumfreq = cumsum(freq),
                    prob = freq ./ sum(freq), cumprob = cumsum(freq) ./ sum(freq))
    return SmoothedFreqTab(tab, F.raw, F.interval, fit1)
end

"""
    presmoothing(F::EG, degree::Int64)
Examine various degree of freedom for fitting the log liner models to the raw score histogram.

# Examples
```julia
julia> tab = freqtab(expandtable(ACTmath.scale, ACTmath.xcount))
Equate.FreqTab(40×5 DataFrame
│ Row │ scale   │ freq  │ cumfreq │ prob       │ cumprob  │
│     │ Float64 │ Int64 │ Int64   │ Float64    │ Float64  │
├─────┼─────────┼───────┼─────────┼────────────┼──────────┤
│ 1   │ 1.0     │ 1     │ 1       │ 0.000231   │ 0.000231 │
│ 2   │ 2.0     │ 1     │ 2       │ 0.000231   │ 0.000462 │
⋮
│ 38  │ 38.0    │ 38    │ 4291    │ 0.00877801 │ 0.991222 │
│ 39  │ 39.0    │ 23    │ 4314    │ 0.00531301 │ 0.996535 │
│ 40  │ 40.0    │ 15    │ 4329    │ 0.003465   │ 1.0      │, [1, 2, 3, 3, 3, 4, 4, 4, 4, 4  …  40, 40, 40, 40, 40, 40, 40, 40, 40, 40], 1.0)

julia> fit = presmoothing(tab, 10)
Fitting dof = 1 model.
Fitting dof = 2 model.
Fitting dof = 3 model.
Fitting dof = 4 model.
Fitting dof = 5 model.
Fitting dof = 6 model.
Fitting dof = 7 model.
Fitting dof = 8 model.
Fitting dof = 9 model.
Fitting dof = 10 model.
10×7 DataFrame. Omitted printing of 1 columns
│ Row │ degree │ aic     │ aicc    │ bic     │ loglik   │ deviance │
│     │ Int64  │ Float64 │ Float64 │ Float64 │ Float64  │ Float64  │
├─────┼────────┼─────────┼─────────┼─────────┼──────────┼──────────┤
│ 1   │ 1      │ 2235.18 │ 2235.5  │ 2238.56 │ -1115.59 │ 1988.28  │
│ 2   │ 2      │ 661.558 │ 662.225 │ 666.625 │ -327.779 │ 412.654  │
⋮
│ 8   │ 8      │ 290.675 │ 296.675 │ 305.875 │ -136.338 │ 29.7704  │
│ 9   │ 9      │ 292.573 │ 300.159 │ 309.462 │ -136.286 │ 29.6683  │
│ 10  │ 10     │ 294.051 │ 303.48  │ 312.629 │ -136.026 │ 29.1466  │

```
"""
function presmoothing(F::EG, degree::Int64)
    if degree > 11
        throw(ArgumentError("Maximum degree for examining must be less than ore equal to 10."))
    end
    fmls = (fml₁, fml₂, fml₃, fml₄, fml₅, fml₆, fml₇, fml₈, fml₉, fml₁₀)
    AIC, BIC, AICC, DEVIANCE, LL = zeros.(Float64, fill(degree, 5))
    FIT = Any[]
    for d in 1:degree
        println("Fitting dof = $d model.")
        fit = glm(fmls[d], F.table, Poisson(), LogLink())
        push!(FIT, fit)
        AIC[d] = StatsBase.aic(fit)
        BIC[d] = StatsBase.bic(fit)
        AICC[d] = StatsBase.aicc(fit)
        LL[d] = StatsBase.loglikelihood(fit)
        DEVIANCE[d] = StatsBase.deviance(fit)
    end
    DataFrame(degree = [1:1:degree;], aic = AIC, aicc = AICC, bic = BIC, loglik = LL, deviance = DEVIANCE, fit = FIT)
end

mutable struct SmoothedNEATFreqTab <: NEAT
    tableX
    tableV
    rawX
    rawV
    intervalX
    intervalV
    fitX
    fitV
    marginal
end
"""
    presmoothing(F::NEAT, LogLinearFormula(df::Int64), LogLinearFormula(df::Int64))

Returns presmoothed frequency table as `SmoothedNEATFreqTab` and `glm` fitted object,
but `interval` and `marginal` elements, which are returned, are not based on smoothed score.

Preserving first C moments of original frequency data, passed `LogLinearFormula(C)` to `fml`.
C is a degree of freedom
"""
function presmoothing(F::NEAT, fmlX, fmlV)
    # Smoothing X (independent part)
    fitX = glm(fmlX, F.tableX, Poisson(), LogLink())
    freqX = predict(fitX, DataFrame(scale = F.tableX.scale))
    tabX = DataFrame(scale = F.tableX.scale, freq = freqX, cumfreq = cumsum(freqX),
                    prob = freqX ./ sum(freqX), cumprob = cumsum(freqX) ./ sum(freqX))
    # Smoothing V (common part)
    fitV = glm(fmlV, F.tableV, Poisson(), LogLink())
    freqV = predict(fitV, DataFrame(scale = F.tableV.scale))
    tabV = DataFrame(scale = F.tableV.scale, freq = freqV, cumfreq = cumsum(freqV),
                    prob = freqV ./ sum(freqV), cumprob = cumsum(freqV) ./ sum(freqV))
    return SmoothedNEATFreqTab(tabX, tabV, F.rawX, F.rawV, F.intervalX, F.intervalV, fitX, fitV, F.marginal)
end

# Kernel method
function RjX(x, xⱼ, a, μ, hX)
    return (x - a * xⱼ - (1-a)*μ) / (a*hX)
end
struct KernelFreqTab <: EG
    table
    raw
    interval
    Bandwidth
end
"""
    KernelSmoothing(X::SmoothedFreqTab; kernel = :Gaussian, hX = 0.622, scale = X.table.scale)

Returns Kernel smoothed frequency table as `KernelFreqTab`.

# Arguments

- `X` the object which class is `freqtab`

## Optional Arguments

- `kernel`
- `hX` The band width.
- `newint` A new interval to depict the continuized probability line.

# Value

`KernelFreqTab`

In the kernel smoothing, the choice of bandwidth `hX` is an important consideration. When bandwidth is large, the smoothed distribution becomes linear equating fucnction. Contrary, bandwidth is small, it becomes spikye shaped distribution. We recommend to use default value `hX = 6.22` or select it by using `EstBandwidth`.
"""
function KernelSmoothing(X::EG; kernel = :Gaussian, hX = 0.622, newint = X.interval/100)
    scale = (minimum(X.table.scale) - X.interval/2):newint:(maximum(X.table.scale) + X.interval/2)
    # hX = bandwidth of cumulative distribution function
    N = sum(X.table.freq)
    μ = mean(X.raw); σ² = var(X.raw)
    a² = σ² / (σ² + hX^2)
    a = sqrt(a²)
    𝒇hX = zeros(Float64, length(scale))
    FhX = zeros(Float64, length(scale))
    for (i, x) in enumerate(scale), (j, xⱼ) in enumerate(X.table.scale)
        𝒇hX[i] += X.table.prob[j]*pdf.(Normal(0, 1), RjX(x, xⱼ, a, μ, hX)) / (a*hX)
        FhX[i] += X.table.prob[j]*cdf.(Normal(0, 1), RjX(x, xⱼ, a, μ, hX))
    end
    𝒇hX = 𝒇hX ./ sum(𝒇hX) # normalize
    tbl = DataFrame(scale = scale, freq = N .* 𝒇hX, cunfreq = cumsum(N .* 𝒇hX), prob = 𝒇hX, cumprob = cumsum(𝒇hX))
    return KernelFreqTab(tbl, X.raw, X.interval, hX)
end

function BandwidthPenalty(hX, X::SmoothedFreqTab, lprob, rprob, K; kernel = :Gaussian)
    hX = exp(hX[1])
    r = X.table.prob
    μ = mean(X.raw); σ² = var(X.raw)
    a² = σ² / (σ² + hX^2)
    a =sqrt(a²)
    𝒇hX = zeros(Float64, length(X.table.scale))
    r𝒇′hX = zeros(Float64, length(X.table.scale))
    l𝒇′hX = zeros(Float64, length(X.table.scale))
    for (i, x) in enumerate(X.table.scale), (j, xⱼ) in enumerate(X.table.scale)
        R = RjX(x, xⱼ, a, μ, hX)
        𝒇hX[i] += X.table.prob[j]*pdf.(Normal(0, 1), R) / (a*hX)
        r𝒇′hX[i] -= rprob[j]*pdf.(Normal(0, 1), R) / (a*hX)^2 * R
        l𝒇′hX[i] -= lprob[j]*pdf.(Normal(0, 1), R) / (a*hX)^2 * R
    end
    pen1 = sum((r .- 𝒇hX) .^2)
    pen2 = sum((l𝒇′hX .< 0) .* (r𝒇′hX .> 0))
    return pen1 + K * pen2
end

"""
    EstBandwidth(X::EG; kernel = :Gaussian)

Estimate optimal bandwidth `hX` in Gaussian kernel.

# Arguments

- `X` is the object which the class is `freqtab`.

# Value

Optimize information.

"""
function EstBandwidth(X::SmoothedFreqTab; kernel = :Gaussian, K = 1)
    x = X.table.scale
    lscore = x .- X.interval/4
    rscore = x .+ X.interval/4
    lpred = predict(X.fit, DataFrame(scale = lscore))
    rpred = predict(X.fit, DataFrame(scale = rscore))
    opt = optimize(hX -> BandwidthPenalty(hX, X, lpred, rpred, K; kernel = kernel), [0.5], method = BFGS())
    println("Minimizer sould be transformed `exp()` before interpretation. \nMinimizer (after taking `exp()`) = $(exp(opt.minimizer[1]))\n\n")
    return opt
end


text:
"""
$(TYPEDEF)

A homogeneous transformation matrix representing the transformation from one
three-dimensional Cartesian coordinate system to another.
"""
struct Transform3D{T}
    mat::SMatrix{4, 4, T, 16}
    from::CartesianFrame3D
    to::CartesianFrame3D

    @inline function Transform3D(from::CartesianFrame3D, to::CartesianFrame3D, mat::AbstractMatrix{T}) where T
        new{T}(mat, from, to)
    end
end

Base.eltype(::Type{Transform3D{T}}) where {T} = T

@inline function Transform3D(from::CartesianFrame3D, to::CartesianFrame3D, rot::Rotation{3}, trans::SVector{3})
    T = promote_type(eltype(typeof(rot)), eltype(typeof(trans)))
    R = convert(RotMatrix3{T}, rot)
    @inbounds mat = @SMatrix [R[1] R[4] R[7] trans[1];
                              R[2] R[5] R[8] trans[2];
                              R[3] R[6] R[9] trans[3];
                              zero(T) zero(T) zero(T) one(T)]
   Transform3D(from, to, mat)
end

@inline function Transform3D(from::CartesianFrame3D, to::CartesianFrame3D, rot::Rotation{3, T}) where {T}
    R = convert(RotMatrix3{T}, rot)
    @inbounds mat = @SMatrix [R[1] R[4] R[7] zero(T);
                              R[2] R[5] R[8] zero(T);
                              R[3] R[6] R[9] zero(T);
                              zero(T) zero(T) zero(T) one(T)]
   Transform3D(from, to, mat)
end

@inline function Transform3D(from::CartesianFrame3D, to::CartesianFrame3D, trans::SVector{3, T}) where {T}
    @inbounds mat = @SMatrix [one(T) zero(T) zero(T) trans[1];
                              zero(T) one(T) zero(T) trans[2];
                              zero(T) zero(T) one(T) trans[3];
                              zero(T) zero(T) zero(T) one(T)]
    Transform3D(from, to, mat)
end

@inline Base.convert(::Type{Transform3D{T}}, t::Transform3D{T}) where {T} = t
@inline Base.convert(::Type{Transform3D{T}}, t::Transform3D) where {T} = Transform3D(t.from, t.to, similar_type(t.mat, T)(t.mat))

@inline rotation(t::Transform3D) = @inbounds return RotMatrix(t.mat[1], t.mat[2], t.mat[3], t.mat[5], t.mat[6], t.mat[7], t.mat[9], t.mat[10], t.mat[11])
@inline translation(t::Transform3D) = @inbounds return SVector(t.mat[13], t.mat[14], t.mat[15])

function Base.show(io::IO, t::Transform3D)
    println(io, "Transform3D from \"$(string(t.from))\" to \"$(string(t.to))\":")
    angle_axis = AngleAxis(rotation(t))
    angle = rotation_angle(angle_axis)
    axis = rotation_axis(angle_axis)
    print(io, "rotation: $(angle) rad about $(axis), translation: $(translation(t))") # TODO: use fixed Quaternions.jl version once it's updated
end

@inline function Base.:*(t1::Transform3D, t2::Transform3D)
    @framecheck(t1.from, t2.to)
    mat = t1.mat * t2.mat
    Transform3D(t2.from, t1.to, mat)
end

@inline function LinearAlgebra.inv(t::Transform3D)
    rotinv = inv(rotation(t))
    Transform3D(t.to, t.from, rotinv, -(rotinv * translation(t)))
end

Base.one(::Type{Transform3D{T}}, from::CartesianFrame3D, to::CartesianFrame3D) where {T} =
    Transform3D(from, to, one(SMatrix{4, 4, T}))
Base.one(::Type{Transform3D}, from::CartesianFrame3D, to::CartesianFrame3D) = one(Transform3D{Float64}, from, to)
Base.one(::Type{T}, frame::CartesianFrame3D) where {T<:Transform3D} = one(T, frame, frame)

function Random.rand(::Type{Transform3D{T}}, from::CartesianFrame3D, to::CartesianFrame3D) where T
    rot = rand(RotMatrix3{T})
    trans = rand(SVector{3, T})
    Transform3D(from, to, rot, trans)
end

Random.rand(::Type{Transform3D}, from::CartesianFrame3D, to::CartesianFrame3D) = rand(Transform3D{Float64}, from, to)

function Base.isapprox(x::Transform3D, y::Transform3D; atol::Real = 1e-12)
    x.from == y.from && x.to == y.to && isapprox(rotation(x), rotation(y), atol = atol) && isapprox(translation(x), translation(y), atol = atol)
end


text:
module Histograms
using Plots
using JLD
import Base: +
export H1D, H2D, hfill!, hsave, add!

type Histogram1D
    """
    This is a one-dimensional histogram with full outlier handling.
    The convention for numbering the bins is:
    start from 2, underflow bin is bin 1, overflow bin is nBins+2.
    """
    binEdges
    entries
    weights
    weights_squared
    torques
    inertials
end
# the edges must be sorted
# we add two bins, for underflow and overflow
Histogram1D(nBins, min, max) = Histogram1D(
      # the last edge is the overflow
      # linspace gives us nBins _edges_
      # we'll add 1 to get the right number of bins
      linspace(min, max, nBins+1)
    , zeros(Int64, nBins+2)
    , zeros(nBins+2)
    , zeros(nBins+2)
    , zeros(nBins+2)
    , zeros(nBins+2)
)
H1D(nBins, min, max) = Histogram1D(nBins, min, max)

type Histogram2D
    # this is just a list of vectors (x, y)
    binEdges
    entries
    weights
    weights_squared
    torquesX
    torquesY
    inertialsX
    inertialsY
end
H2D(nBinsX, minX, maxX, nBinsY, minY, maxY) = Histogram2D(nBinsX, minX, maxX, nBinsY, minY, maxY)
Histogram2D() = Histogram2D(0,0.0,0.0, 0,0.0,0.0)
Histogram2D(nBinsX, minX, maxX, nBinsY, minY, maxY) = Histogram2D(
      (linspace(minX, maxX, nBinsX+1), linspace(minY, maxY, nBinsY+1))
    , zeros(Int64, (nBinsX+2, nBinsY+2))
    , zeros((nBinsX+2, nBinsY+2))
    , zeros((nBinsX+2, nBinsY+2))
    , zeros((nBinsX+2, nBinsY+2))
    , zeros((nBinsX+2, nBinsY+2))
    , zeros((nBinsX+2, nBinsY+2))
    , zeros((nBinsX+2, nBinsY+2))
)

function find_index(h::Histogram1D, x):
    """
    Returns the index of the bin that this item would be put into
    """
    # underflow
    if x < h.binEdges[1]
        return 1
    end
# we have two extra bins, one for underflow(1) and one for overflow(end+1)
# this means we have one more bin than we have edges
    @inbounds for i in 2:length(h.binEdges)
        if x < h.binEdges[i]
            return i
        end
    end
    # overflow
    return length(h.binEdges)+1
end

function hfill!(h::Histogram1D, x, weight=1.)
    index = find_index(h, x)
    h.entries[index] += 1
    h.weights[index] += weight
    h.weights_squared[index] += weight*weight
    h.torques[index] += x*weight
    h.inertials[index] += x*x*weight
    return
end

stddev(h::Histogram1D) = sum(h.torques)
mean(h::Histogram1D) = sum(h.entries)
function +(h1::Histogram1D, h2::Histogram1D)
    if ! all(h1.binEdges == h2.binEdges)
        println("ERROR! The bin egdes must be the same for both histograms")
        return
    end
    entries = h1.entries + h2.entries
    weights = h1.weights + h2.weights
    weights_squared = h1.weights_squared + h2.weights_squared
    return Histogram1D(edges, entries, weights, weights_squared)
end

function +(h1::Histogram2D, h2::Histogram2D)
    if ! all(h1.binEdges[1] == h2.binEdges[1])
        println("ERROR! The bin egdes must be the same for both histograms")
        return Histogram2D()
    end
    if ! all(h1.binEdges[2] == h2.binEdges[2])
        println("ERROR! The bin egdes must be the same for both histograms")
        return Histogram2D()
    end
    entries = h1.entries + h2.entries
    weights = h1.weights + h2.weights
    weights_squared = h1.weights_squared + h2.weights_squared
    torquesX = h1.torquesX + h2.torquesX
    torquesY = h1.torquesY + h2.torquesY
    inertialsX = h1.inertialsX + h2.inertialsX
    inertialsY = h1.inertialsY + h2.inertialsY
    return Histogram2D(entries, weights, weights_squared, torquesX, torquesY, inertialsX, inertialsY)
end

function add!(h1::Histogram2D, h2::Histogram2D)
    if ! all(h1.binEdges[1] == h2.binEdges[1])
        println("ERROR! The bin egdes must be the same for both histograms")
        return
    end
    if ! all(h1.binEdges[2] == h2.binEdges[2])
        println("ERROR! The bin egdes must be the same for both histograms")
        return
    end
    h1.entries .+= h2.entries
    h1.weights .+= h2.weights
    h1.weights_squared .+= h2.weights_squared
    h1.torquesX .+= h2.torquesX
    h1.torquesY .+= h2.torquesY
    h1.inertialsX .+= h2.inertialsX
    h1.inertialsY .+= h2.inertialsY
    return
end

function find_index(h::Histogram2D, x, y)
    # underflow
    i = 0
    if x < h.binEdges[1][1]
        i = 1
    end
    j = 0
    if y < h.binEdges[2][1]
        j = 1
    end
# we have two extra bins, one for underflow(1) and one for overflow(end+1)
# this means we have one more bin than we have edges
    if i == 0
        @inbounds for i in 2:length(h.binEdges[1])
            if x < h.binEdges[1][i]
                break
            end
        end
    end
    # we shift the index of the binEgdes by starting at 1
    # thats why we return i rather than i-1
    if j == 0
        @inbounds for j in 2:length(h.binEdges[2])
            if y < h.binEdges[2][j]
                break
            end
        end
    end
    if i == 0
        i = length(h.binEdges[1]) + 1
    end
    if j == 0
        j = length(h.binEdges[2]) + 1
    end
    return i,j
end

function hfill!(h::Histogram2D, x, y, weight = 1.0)
    i,j = find_index(h, x, y)
    h.entries[i,j] += 1
    h.weights[i,j] += weight
    h.weights_squared[i,j] += weight^2
    h.torquesX[i,j] += x * weight
    h.torquesY[i,j] += y * weight
    h.inertialsX[i,j] += x^2 * weight
    h.inertialsY[i,j] += y^2 * weight
end

# only use the in-range bins to compute the statistics
function getStatisticSet(h::Histogram2D)
    weights = sum(h.weights[2:end-1])
	torquesX = sum(h.torquesX[2:end-1])
	torquesY = sum(h.torquesY[2:end-1])

	entries = h.entries[2:end-1]
	meanX = torquesX / weights
	meanY = torquesY / weights

    rmsX = sqrt((sum(h.inertialsX[2:end-1]) - torquesX^2 / weights) / weights)
	rmsY = sqrt((sum(h.inertialsY[2:end-1]) - torquesY^2 / weights) / weights)
	return entries, meanX, meanY, rmsX, rmsY
end

function plot(h::Histogram1D)
    prinln("Not yet implemented")
end

function plot(h::Histogram2D)
    xs = h.binEdges[1]
    ys = h.binEdges[2]
    z = h.weights'[1:end-1, 1:end-1]
    heatmap(xs,ys,z)
end

function hsave(h::Histogram1D, histname, filename)
    jldopen(filename, "w") do file
        addrequire(file, Histograms)
        write(file, histname, h)
    end
end

function hsave(h::Histogram2D, histname, filename)
    jldopen(filename, "w") do file
        addrequire(file, Histograms)
        write(file, histname, h)
    end
end

end


text:
"""
type-II discrete cosine transform (DCT)
"""
function dct!(data::Array{T, 3}) where {T <:AbstractFloat}
	dctmat = genmat(size(data, 1))
	dctmat1 = genmat(size(data, 2))'
	@views @inbounds @fastmath for i in 1:size(data, 3)
		data[:, :, i] .= dctmat * data[:, :, i] * dctmat1
	end

	dctmat = genmat(size(data, 3))
	@views @inbounds @fastmath for j in 1:size(data, 2), i in 1:size(data, 1)
		data[i, j, :] .= dctmat * data[i, j, :]
	end
end

"""
inverse discrete cosine transform (DCT)
"""
function idct!(data::Array{T, 3}) where {T <:AbstractFloat}
	dctmat = genmat(size(data, 3))'
	@views @inbounds @fastmath for j in 1:size(data, 2), i in 1:size(data, 1)
		data[i, j, :] .= dctmat * data[i, j, :]
	end

	dctmat = genmat(size(data, 1))'
	dctmat1 = genmat(size(data, 2))
	@views @inbounds @fastmath for i in 1:size(data, 3)
		data[:, :, i] .= dctmat * data[:, :, i] * dctmat1
	end
end

function genmat(h::Int)::Array{Float64, 2}
	mat = Array{Float64, 2}(undef, h, h)
	mat[1, :] .= sqrt(1 / h)
	@inbounds @fastmath for j in 1:h
		@simd for i in 2:h
			mat[i, j] = cos(π * (i - 1) * (j - 0.5) / h) * sqrt(2 / h)
		end
	end
	mat
end


text:
struct FilterMap <: Map end

function return_cache(k::FilterMap,f,a)
  # vals = testitem(a)
  vals = a
  T = eltype(eltype(a))
  r = zeros(T,length(vals))
  c = CachedArray(r)
end

function evaluate!(cache,k::FilterMap,f,a)
  c = cache
  vals = a
  filters = f
  @check size(vals) == size(filters) "Local arrays mismatch"
  setsize!(c,(sum(filters),))
  r = c.array
  i = 0
  for (val,filter) in zip(vals,filters)
    if filter
      i += 1
      r[i] = val
    end
  end
  r
end


text:
# Multi-threading in Julia
# Multi-threading will be native to Julia in the next release (v0.5.x). It is currently under active development.
# This tutorial gives you a quick taste of how to use multi-threading in Julia.

# First load the threading code that's in the Threads module.

using Base.Threads

@inline function cndf2(in::Array{Float64,1})
    out = 0.5 .+ 0.5 .* erf(0.707106781 .* in)
    return out
end

function blackscholes_serial(sptprice::Float64,
                           strike::Array{Float64,1},
                           rate::Float64,
                           volatility::Float64,
                           time::Float64)
    logterm = log10(sptprice ./ strike)
    powterm = .5 .* volatility .* volatility
    den = volatility .* sqrt(time)
    d1 = (((rate .+ powterm) .* time) .+ logterm) ./ den
    d2 = d1 .- den
    NofXd1 = cndf2(d1)
    NofXd2 = cndf2(d2)
    futureValue = strike .* exp(- rate .* time)
    c1 = futureValue .* NofXd2
    call = sptprice .* NofXd1 .- c1
    put  = call .- futureValue .+ sptprice
end

# This next parallel rewrite does two simple things:
# It devectorizes the entire loop. We're now looping through arrays instead of dealing with them in blocks like we did previously. This is devectorized code while the above serial block is vectorized code.

function blackscholes_devec(sptprice::Float64,
				strike::Vector{Float64},
				rate::Float64,
				volatility::Float64,
				time::Float64)
    sqt = sqrt(time)
    put = similar(strike)
    for i = 1:size(strike, 1)
        logterm = log10(sptprice / strike[i])
        powterm = 0.5 * volatility * volatility
        den = volatility * sqt
        d1 = (((rate + powterm) * time) + logterm) / den
        d2 = d1 - den
        NofXd1 = 0.5 + 0.5 * erf(0.707106781 * d1)
        NofXd2 = 0.5 + 0.5 * erf(0.707106781 * d2)
        futureValue = strike[i] * exp(-rate * time)
        c1 = futureValue * NofXd2
        call_ = sptprice * NofXd1 - c1
        put[i] = call_ - futureValue + sptprice
    end
    put
end

# Affixes the `@threads` macro in front of the `for` to tell Julia that this is a multi-threaded block.

function blackscholes_parallel(sptprice::Float64,
				strike::Vector{Float64},
				rate::Float64,
				volatility::Float64,
				time::Float64)
    sqt = sqrt(time)
    put = similar(strike)
    @threads for i = 1:size(strike, 1)
        logterm = log10(sptprice / strike[i])
        powterm = 0.5 * volatility * volatility
        den = volatility * sqt
        d1 = (((rate + powterm) * time) + logterm) / den
        d2 = d1 - den
        NofXd1 = 0.5 + 0.5 * erf(0.707106781 * d1)
        NofXd2 = 0.5 + 0.5 * erf(0.707106781 * d2)
        futureValue = strike[i] * exp(-rate * time)
        c1 = futureValue * NofXd2
        call_ = sptprice * NofXd1 - c1
        put[i] = call_ - futureValue + sptprice
    end
    put
end

function run(iterations)
    sptprice   = 42.0
    initStrike = Float64[ 40.0 + (i / iterations) for i = 1:iterations ]
    rate       = 0.5
    volatility = 0.2
    time       = 0.5

    tic()
    put1 = blackscholes_serial(sptprice, initStrike, rate, volatility, time)
    t1 = toq()
    println("Serial checksum: ", sum(put1))
    tic()
    put2 = blackscholes_devec(sptprice, initStrike, rate, volatility, time)
    t2 = toq()
    println("Parallel checksum: ", sum(put2))
    tic()
    put3 = blackscholes_parallel(sptprice, initStrike, rate, volatility, time)
    t3 = toq()
    println("Parallel checksum: ", sum(put3))
    return t1, t2, t3
end

function driver()
    srand(0)
    tic()
    iterations = 10^6
    blackscholes_serial(0., Float64[], 0., 0., 0.)
    blackscholes_devec(0., Float64[], 0., 0., 0.)
    blackscholes_parallel(0., Float64[], 0., 0., 0.)
    println("SELFPRIMED ", toq())
    tserial, tdevec, tparallel = run(iterations)
    println("Time taken for serial = $tserial")
    println("Time taken for devec = $tdevec")
    println("Time taken for parallel = $tparallel")
    println("Serial rate = ", iterations / tserial, " opts/sec")
    println("Devec rate = ", iterations / tdevec, " opts/sec")
    println("Parallel rate = ", iterations / tparallel, " opts/sec")
end
driver()

max_issues_repo_name: andreasnoack/ParallelWorkshop


text:
#2020 Klaus Stricker - Tests for ReadWriteDlm2
#License is MIT: http://julialang.org/license

# rwd2tests_3.jl

@testset "3_  other types" begin

    # Test Complex and Rational parsing
    a = Any[complex(-1,-2) complex(1.2,-2) complex(1e9,3e19) 1//3 -1//5 -2//-4 1//-0 -0//1]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "-1-2im;1,2-2,0im;1,0e9+3,0e19im;1//3;-1//5;1//2;1//0;0//1\n"
    b = readdlm2("test.csv", Any)
    rm("test.csv")
    @test isequaldlm(a, b, Any)

    # Complex and Rational - tolerance with blancs, i/j and different signes
    write("test.csv", "    -1-2j; 1,2 - 2,0i ;1.0E9+3.0E19im;   -1//-3;1//-5;  1//2 ;1//-0;-0//1 \n")
    b = readdlm2("test.csv", Any)
    rm("test.csv")
    @test b == a

    # Test Complex and Rational parsing - decimal = '.', delimiter = \t
    a = Any[complex(-1,-2) complex(1.2,-2) complex(1e9,3e19) 1//3 -1//5 -2//-4 1//-0 -0//1]
    writedlm2("test.csv", a, '\t', decimal='.')
    @test read("test.csv", String) == "-1-2im\t1.2-2.0im\t1.0e9+3.0e19im\t1//3\t-1//5\t1//2\t1//0\t0//1\n"
    b = readdlm2("test.csv",'\t', Any, decimal='.')
    rm("test.csv")
    @test isequaldlm(a, b, Any)

    #  Test different types with header and Any Array - decimal comma
    a = Any["Nr" "Wert";1 Date(2017);2 DateTime(2018);3 Dates.Time(23,54,45,123,456,78);4 1.5e10+5im;5 1500//5;6 1.5e10]
    writedlm2("test.csv", a)
    @test read("test.csv", String) ==
    """
    Nr;Wert
    1;2017-01-01
    2;2018-01-01T00:00:00.0
    3;23:54:45,123456078
    4;1,5e10+5,0im
    5;300//1
    6;1,5e10
    """
    b = readdlm2("test.csv", Any, header=true)
    rm("test.csv")
    @test a[2:end,:] == b[1]
    @test a[1:1,:] == b[2]

    #  Test different types with header and Any Array - english decimal
    a = Any["Nr" "Value";1 Date(2017);2 DateTime(2018);3 Dates.Time(23,54,45,123,456,78);4 1.5e10+5im;5 1500//5;6 1.5e10]
    writedlm2("test.csv", a, '\t', decimal='.')
    @test read("test.csv", String) ==
    """
    Nr\tValue
    1\t2017-01-01
    2\t2018-01-01T00:00:00.0
    3\t23:54:45.123456078
    4\t1.5e10+5.0im
    5\t300//1
    6\t1.5e10
    """
    b = readdlm2("test.csv", '\t', Any, decimal='.', header=true)
    rm("test.csv")
    @test a[2:end,:] == b[1]
    @test a[1:1,:] == b[2]

    # Test Complex Array read and write
    a = Complex[complex(-1,-2) complex(1.2,-2) complex(-1e9,3e-19)]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "-1-2im;1,2-2,0im;-1,0e9+3,0e-19im\n"
    b = readdlm2("test.csv", Complex)
    rm("test.csv")
    @test a == b
    @test typeof(b) == Array{Complex,2}

    a = Complex[complex(-1,-2) complex(1.2,-2) complex(-1e9,3e-19)]
    writedlm2("test.csv", a, imsuffix="i")
    @test read("test.csv", String) == "-1-2i;1,2-2,0i;-1,0e9+3,0e-19i\n"
    b = readdlm2("test.csv", Complex)
    rm("test.csv")
    @test a == b
    @test typeof(b) == Array{Complex,2}

    a = Any["test" "test2";complex(-1,-2) complex(1.2,-2);complex(-1e9,3e-19) complex(1,115)]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "test;test2\n-1-2im;1,2-2,0im\n-1,0e9+3,0e-19im;1+115im\n"
    b = readdlm2("test.csv", Any)
    rm("test.csv")
    @test a == b
    @test typeof(b) == Array{Any,2}

    a = Any["test1" "test2";complex(-1,-2) complex(1.2,-2);complex(-1e9,3e-19) complex(1,115)]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "test1;test2\n-1-2im;1,2-2,0im\n-1,0e9+3,0e-19im;1+115im\n"
    b, h = readdlm2("test.csv", Complex, header=true)
    rm("test.csv")
    @test b == Complex[complex(-1,-2) complex(1.2,-2);complex(-1e9,3e-19) complex(1,115)]
    @test h == AbstractString["test1" "test2"]

    a = Complex[complex(-1//3,-7//5) complex(1,-1//3) complex(-1//2,3e-19)]
    writedlm2("test.csv", a, imsuffix="i")
    @test read("test.csv", String) == "-1//3-7//5*i;1//1-1//3*i;-0,5+3,0e-19i\n"
    b = readdlm2("test.csv", Complex)
    rm("test.csv")
    @test a == b
    @test typeof(b) == Array{Complex,2}

    a = Complex[complex(-1//3,-7//5) complex(1,-1//3) complex(-1//2,3e-19)]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "-1//3-7//5*im;1//1-1//3*im;-0,5+3,0e-19im\n"
    b = readdlm2("test.csv", Complex)
    rm("test.csv")
    @test a == b
    @test typeof(b) == Array{Complex,2}

    # Test Time read and write
    a = Time[Time(12,54,43,123,456,789) Time(13,45);Time(1,45,58,0,0,1) Time(23,59,59)]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "12:54:43,123456789;13:45:00\n01:45:58,000000001;23:59:59\n"
    b = readdlm2("test.csv", Time)
    rm("test.csv")
    @test b == a
    @test typeof(b) == Array{Time,2}

    a = Any["test1" "test2";Time(12,54,43,123,456,789) Time(13,45);Time(1,45,58,0,0,1) Time(23,59,59)]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "test1;test2\n12:54:43,123456789;13:45:00\n01:45:58,000000001;23:59:59\n"
    b, h = readdlm2("test.csv", Time, header=true)
    rm("test.csv")
    @test b == Time[Time(12,54,43,123,456,789) Time(13,45);Time(1,45,58,0,0,1) Time(23,59,59)]
    @test h == AbstractString["test1" "test2"]

    # Test Rational read and write
    a = Rational[456//123 123//45;456//23 1203//45]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "152//41;41//15\n456//23;401//15\n"
    b = readdlm2("test.csv", Rational)
    rm("test.csv")
    @test b == a
    @test typeof(b) == Array{Rational,2}

    a = Any["test1" "test2";456//123 123//45;456//23 1203//45]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "test1;test2\n152//41;41//15\n456//23;401//15\n"
    b, h = readdlm2("test.csv", Rational, header=true)
    rm("test.csv")
    @test b == Rational[456//123 123//45;456//23 1203//45]
    @test h == AbstractString["test1" "test2"]

    # Test DateTime read and write
    a = DateTime[DateTime(2017) DateTime(2016);DateTime(2015) DateTime(2014)]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "2017-01-01T00:00:00.0;2016-01-01T00:00:00.0\n2015-01-01T00:00:00.0;2014-01-01T00:00:00.0\n"
    b = readdlm2("test.csv", DateTime)
    c = readdlm2("test.csv", DateTime, tables=true) # test tables interface
    @test ReadWriteDlm2.mttoarray(c) == a
    rm("test.csv")
    @test b == a
    @test typeof(b) == Array{DateTime,2}

    a = Any["test1" "test2";DateTime(2017) DateTime(2016);DateTime(2015) DateTime(2014)]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "test1;test2\n2017-01-01T00:00:00.0;2016-01-01T00:00:00.0\n2015-01-01T00:00:00.0;2014-01-01T00:00:00.0\n"
    b, h = readdlm2("test.csv", DateTime, header=true)
    rm("test.csv")
    @test b == DateTime[DateTime(2017) DateTime(2016);DateTime(2015) DateTime(2014)]
    @test h == AbstractString["test1" "test2"]

    # Test Date read and write
    a = Date[Date(2017) Date(2016);Date(2015) Date(2014)]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "2017-01-01;2016-01-01\n2015-01-01;2014-01-01\n"
    b = readdlm2("test.csv", Date)
    rm("test.csv")
    @test b == a
    @test typeof(b) == Array{Date,2}

    a = Any["test1" "test2";Date(2017) Date(2016);Date(2015) Date(2014)]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "test1;test2\n2017-01-01;2016-01-01\n2015-01-01;2014-01-01\n"
    b, h = readdlm2("test.csv", Date, header=true)
    rm("test.csv")
    @test b == Date[Date(2017) Date(2016);Date(2015) Date(2014)]
    @test h == AbstractString["test1" "test2"]

    # Test alternative rs Tuple - is decimal ignored?
    a = Float64[1.1 1.2;2.1 2.2]
    writedlm2("test.csv", a, decimal='€')
    @test read("test.csv", String) == "1€1;1€2\n2€1;2€2\n"
    b = readdlm2("test.csv", Any, rs=(r"(\d)€(\d)", s"\1.\2"), decimal='n')
    rm("test.csv")
    @test a == b
end


text:
# This file contains algorithms needed for the integratiorn of
# transcendental functions from chapter 5 of the book
#
#    Manuel Bronstein, Symbolic Integration I: Transcendental Functions,
#    2nd ed., Springer 2005.
#

"""
    HermiteReduce(f, D) -> (g, h, r)

Hermite reduction - quadratic version.

Given a field `k`, a derivation `D` on `k(t)` and `f` in `k(t)`, return `g`, `h`, `r` in `k(t)`
such that `f=D(g)+h+r`, `h` is simple and `r` is reduced.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 5.4, p. 139.
"""
function HermiteReduce(f::FracElem{P}, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}}
    iscompatible(f, D) || error("rational function f must be in the domain of derivation D")
    fp, fs, fn = CanonicalRepresentation(f, D)
    a = numerator(fn)
    d = denominator(fn)
    ds = Squarefree(d)
    g = zero(f)
    for i=2:length(ds)
        if degree(ds[i])>0
            v = ds[i]
            u = divexact(d, v^i)
            for j=i-1:-1:1
                b, c = gcdx(u*D(v), v, -1//j*a)
                g += b//v^j
                a = -j*c - u*D(b)
            end
            d = u*v
        end
    end
    q, r = divrem(a, d)
    g, r//d, q + fp + fs
end

"""
    PolynomialReduce(p, D) -> (q, r)

Polynomial reduction.

Given a field `k`, a derivation `D` on `k(t)` and `p` in `k[t]` where `t` is a nonlinear
monomial over `k`, return `q`, `r` in `k[t]`
such that `p=D(q)+r` and `degree(r)<degree(D(t))`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 5.4, p. 141.
"""
function PolynomialReduce(p::P, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}}
    iscompatible(p, D) || error("polynomial p must be in the domain of derivation D")
    isnonlinear(D) || error("monomial of derivation D must be nonlinear")
    δ = degree(D)
    λ = leading_coefficient(D)
    t = gen(domain(D))
    if degree(p)<δ
        return (zero(p), p)
    end
    m = degree(p) - δ +1
    q0 = (leading_coefficient(p)//(m*λ))*t^m
    q, r = PolynomialReduce(p - D(q0), D)
    q0 + q, r
end

"""
    ResidueReduce(f, D) -> (ss, Ss, ρ)

Polynomial reduction.

Given a field `k`, a derivation `D` on `k(t)` and `f` in `k(t)`, return either `ρ=1`
and `ss=[s₁,...,sₘ]`, `sᵢ` in `k[z]`, `Ss=[S₁,...Sₘ]`, `Sᵢ` in `k[z,t]` where
`z` is a new indeterminate over `k` such that `f-D(g)` is in `k[t]` where
`g=∑ᵢ∑_{α:sᵢ(α)=0}α*log(Sᵢ(α,t))`, or `ρ=0` and such `ss` and `Ss` such that
`f+h` and `f+h-D(g)` do not have an elementary integral over `k(t)` for any
`h` in `k⟨t⟩`.

(Here, `k⟨t⟩` denotes the elements of `k(t)` which are reduced w.r.t. `D`.)

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 5.6, p. 153.
"""
function ResidueReduce(f::F, D::Derivation; symbol=:α) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    iscompatible(f, D) || error("rational function f must be in the domain of derivation D")
    issimple(f, D) || error("rational function f must be simple with respect to derivation D")
    d = denominator(f)
    if isone(d)
        return PolyElem[], PolyElem{PolyElem}[], 1
    end
    p, a = divrem(numerator(f), d)
    # For SubResultant with respect to t we have to construct the
    # polynomial ring k[z][t] with z, t in this order (!)
    kz, z = PolynomialRing(base_ring(d), symbol)
    kzt, t = PolynomialRing(kz, var(parent(d)))
    if degree(D(d))<=degree(d)
        r, Rs = SubResultant(d(t), a(t) - z*D(d)(t))
    else
        r, Rs = SubResultant(a(t) - z*D(d)(t), d(t))
    end
    κD = CoefficientLiftingDerivation(kz, BaseDerivation(D))
    ns, ss = SplitSquarefreeFactor(r, κD)
    ds = degree.(ss)
    ss = typeof(d)[ss[i] for i=1:length(ss) if ds[i]>0]
    Ss = typeof(t)[]
    ii = 1
    for i=1:length(ds)
        if ds[i]>0
            if i==degree(d)
                S = d(t)
            else
                m = findfirst(z->z==i, degree.(Rs))
                S =  Rs[m]
                As = Squarefree(leading_coefficient(S))
                for j=1:length(As)
                    S = divexact(S, gcd(As[j],ss[ii])^j)
                end
            end
            push!(Ss, S)
            ii+=1
        end
    end
    ρ = degree(prod(ns))<=0 ? 1 : 0
    ss, Ss, ρ
end

"""
    ConstantPart(ss, Ss, D) -> (αs, gs, ss1, Ss1)

Given the output `ss=[s₁,...,sₘ]`, `sᵢ` in `k[z]`, `Ss=[S₁,...Sₘ]`, `Sᵢ` in `k[z,t]`
of `ResidueReduce` and the derivation D on k(t) which was used by `ResidueReduce`
return  `αs=[α₁,...,αᵣ]` consisting of the roots of those `sᵢ` that have all its roots in `Const(k)` and the corresponding
`gs=[g₁,...,gᵣ]` where `gⱼ` in `k[t]` and `gⱼ(t) = Sᵢ(αⱼ,t)`. The remaining `sᵢ` that have at least one root not in `Const(k)`
and the corresponding `Sᵢ` are returned in `ss1` and `Ss1`.
"""
function ConstantPart(ss::Vector{P}, Ss::Vector{PP}, D::Derivation) where  {P<:PolyElem, PP<:PolyElem{P}}
    #TODO: better name for this function...
    length(ss)==length(Ss) || error("lengths must match")
    if isempty(ss)
        return Term[], 0,  ss, Ss
    end
    Ss1 = eltype(Ss)[]
    ss1 = eltype(ss)[]
    gs = Term[]
    Dg = 0
    D1 = CoefficientLiftingDerivation(parent(ss[1]), BaseDerivation(D))
    for i=1:length(ss)
        αs = constant_roots(ss[i], D1)
        if length(αs)==degree(ss[i]) # all roots found
            for α in αs
                g = map_coefficients(c->c(α), Ss[i])
                if !isconstant(g, D)
                    push!(gs, FunctionTerm(log, α, g))
                    Dg += α*D(g)//g
                end
            end
        else
            RT = LogToReal(SumOfLogTerms(ss[i], Ss[i]))
            αs = constant_roots(ss[i], D1, useQQBar=true)
            if all([isrational(real(α)) && isrational(imag(α)) for α in αs])
                for α in αs
                    u = rationalize(real(α))
                    v = rationalize(imag(α))
                    if iszero(v)
                        g = map_coefficients(c->c(u), Ss[i])
                        if !isconstant(g, D)
                            push!(gs, FunctionTerm(log, u, g))
                            Dg += u*D(g)//g
                        end
                    elseif v > 0
                        var = string(symbols(parent(Ss[i]))[1])
                        F = base_ring(ss[i])
                        if !iszero(u)
                            # Ignore log terms with contstant arguments. In some cases the denominator of the constant argument
                            # might be zero after substitutiong (u,v). So this avoids divison by zero in these cases.
                            if degree(RT.LT.arg)>0 || (!isconstant(numerator(constant_coefficient(RT.LT.arg))(u,v), BaseDerivation(D)) &&
                                                       !isconstant(denominator(constant_coefficient(RT.LT.arg))(u,v), BaseDerivation(D)))
                                # TODO: Think about avoiding division by zero like below for atan term.
                                # But maybe this case cannot occur anyway, at least I know of no example where it does.
                                g = polynomial(F, [numerator(c)(u, v)//denominator(c)(u, v) for c in coefficients(RT.LT.arg)], var)
                                push!(gs, FunctionTerm(log, RT.LT.coeff*u, g))
                                Dg += RT.LT.coeff*u*D(g)//g
                            end
                        end
                        for AT in RT.ATs
                            # Ignore atan terms with contstant arguments. In some cases the denominator of the constant argument
                            # might be zero after substitutiong (u,v). So this avoids divison by zero in these cases.
                            if degree(AT.arg)>0 || (!isconstant(numerator(constant_coefficient(AT.arg))(u,v), BaseDerivation(D)) &&
                                                    !isconstant(denominator(constant_coefficient(AT.arg))(u,v), BaseDerivation(D)))
                                if all([!iszero(denominator(c)(u, v)) for c in coefficients(AT.arg)])
                                    # Ignore atan term if substitution of (u,v) in argument would cause divion by zero.
                                    # This requires more thought, but it seems to work...
                                    g = polynomial(F, [numerator(c)(u, v)//denominator(c)(u, v) for c in coefficients(AT.arg)], var)
                                    push!(gs, FunctionTerm(atan, AT.coeff*v, g))
                                    Dg += AT.coeff*v*D(g)//(1 + g^2)
                                end
                            end
                        end
                    end
                end
            else
                push!(Ss1, Ss[i])
                push!(ss1, ss[i])
            end
        end
    end
    gs, Dg, ss1, Ss1
end

"""
    IntegratePrimitivePolynomial(p, D) -> (q, ρ)

Integration of polynomials in a primitive extension.

Given a field `k`, a derivation `D` on `k(t)` such that `t` is a primitive monomial over `k`
and `p` in `k[t]`, return either `ρ=1` and `q` in `k[t]` such that `p-D(q)` is in `k`, or
`ρ=0` and  `q` in `k[t]` such that `p-D(q)` does not have an elementary integral over
`k(t)`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 5.8, p. 158.
"""
function IntegratePrimitivePolynomial(p::P, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}}
    iscompatible(p, D) || error("rational function p must be in the domain of derivation D")
    isprimitive(D) || error("monomial of derivation D must be primitive")
    if degree(p)<=0
        return zero(p), 1
    end
    a = leading_coefficient(p)
    b, c, β = LimitedIntegrate(a, leading_coefficient(D), BaseDerivation(D))
    if β<=0
        return zero(p), 0
    end
    m = degree(p)
    t = gen(parent(p))
    q0 = (c//(m+1))*t^(m+1) + b*t^m
    q, β = IntegratePrimitivePolynomial(p-D(q0), D)
    return q+q0, β
end

"""
    IntegrateHyperexponentialPolynomial(p, D) -> (q, ρ)

Integration of hyperexponential polynomials.

Given a field `k`, a derivation `D` on `k(t)` such that `t` is a hyperexponential monomial over `k`
and `p` in `k⟨t⟩=k[t,t⁻¹]`, return either `ρ=1` and `q` in `k[t,t⁻¹]` such that `p-D(q)` is in `k`, or
`ρ=0` and  `q` in `k[t,t⁻¹]` such that `p-D(q)` does not have an elementary integral over
`k(t)`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 5.9, p. 162.
"""
function IntegrateHyperexponentialPolynomial(p::F, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    iscompatible(p, D) || error("rational function p must be in the domain of derivation D")
    ishyperexponential(D) || error("monomial of derivation D must be hyperexponential")
    isreduced(p, D) || error("raional function p must be reduced.")
    q = zero(p)
    ρ = 1
    if iszero(p)
        # Note: with the conventions in Bronstein's book the for loop below would
        # run from i=+infinity to i=-infinity for p=0, i.e., not at all.
        # But it's cleaner to handle the case p=0 separately.
        return q, ρ
    end
    t = gen(base_ring(parent(p)))
    w = coeff(MonomialDerivative(D), 1) # Dt/t
    for i=valuation(p, t):-valuation_infinity(p)
        if i!=0
            a = coeff(numerator(p), i + degree(denominator(p)))
            v, ρ1 = RischDE(i*w, a, BaseDerivation(D))
            if ρ1==0
                ρ = 0
            else
                q += v*(t//1)^i
            end
        end
    end
    q, ρ
end

"""
    IntegrateHypertangentPolynomial(p, D) -> (q, c)

Integration of hypertangent polynomials.

Given a field `k` not conataining `√-1` a derivation `D` on `k(t)` such that `t` is a hypertangent monomial over `k`
and `p` in `k[t]`, return `q` in `k[t]` and `c` in `k` such that `p-D(q)-c*D(t²+1)/(t²+1)` is in `k` and
`p-D(q)` does not have an elementary integral over `k(t)` if `D(c)≠0`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 5.10, p. 167.
"""
function IntegrateHypertangentPolynomial(p::P, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}}
    iscompatible(p, D) || error("rational function p must be in the domain of derivation D")
    ishypertangent(D) || error("monomial of derivation D must be hypertangent")
    t = gen(parent(p))
    q, r = PolynomialReduce(p, D)
    α = constant_coefficient(divexact(MonomialDerivative(D), t^2+1))
    c = coeff(r, 1)//(2*α)
    q, c
end

"""
    IntegrateHypertangentReduced(p, D) -> (q, ρ)

Integration of hypertangent reduced elements.

Given a field `k` not conataining `√-1` a derivation `D` on `k(t)` such that `t` is a hypertangent monomial over `k`
and `p` in `k⟨t⟩`, return either `ρ=1` and `q` in `k⟨t⟩` such that `p-D(q)` is in `k[t]`, or
`ρ=0` and  `q` in `k⟨t⟩` such that  `p-D(q)` does not have an elementary integral over `k(t).

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 5.10, p. 169.
"""
function IntegrateHypertangentReduced(p::F, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    iscompatible(p, D) || error("rational function p must be in the domain of derivation D")
    ishypertangent(D) || error("monomial of derivation D must be hypertangent")
    isreduced(p, D) || error("rational function p must be reduced.")
    Z = zero(parent(p))
    if iszero(p)
        return Z, 1
    end
    t = gen(base_ring(parent(p)))
    Q = t^2+1
    m = -(iszero(p) ? typemax(Int) : valuation(p, Q))
    if m<=0
        return Z, 1
    end
    h = numerator(Q^m*p)
    q, r = divrem(h, Q)
    a = coeff(r, 1)
    b = coeff(r, 0)
    η = constant_coefficient(divexact(MonomialDerivative(D), Q))
    c, d, ρ = CoupledDESystem(zero(parent(a)), 2*m*η, a, b, BaseDerivation(D))
    if ρ<=0
        return Z, 0
    end
    q0 = (c*t + d)//Q^m
    q, ρ = IntegrateHypertangentReduced(p - D(q0), D)
    q + q0, ρ
end

"""
    Integrate(f, D) -> (αs, lgs, g1, r, ρ)

Given a field `k` (with `√-1` not in `k` in the hypertangent case), a derivation
`D` on `k(t)`, and `f` in `k(t)`,
return either `ρ=1`, `αs=[α₁,...αᵣ]` with `αᵢ` in `Const(k)`, `lgs=[lg₁,...lgᵣ]`
with `lgᵢ` in `k[t]`,
and `g1` in `k(t)`  such that `g=g1+∑αᵢ*log(lgᵢ)` is elementary over `k(t)` and `r=f-D(g)` is in `k`, or
`ρ=0` and such `αs`, `lgs` and `g1` such that `r=f-D(g)` does not have an elementary integral over
`k(t)`.

This function corresponds to `IntegreatePrimitive`, `IntegrateHyperexponential` and
`IntegrateTangent` joined together, see [Bronstein's book](https://link.springer.com/book/10.1007/b138171),
Section 5.8 p. 160, Section 5.9 p. 163, and Section 5.10 p. 172, respectively.
"""
function Integrate(f:: F, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    if isbasic(D)
        return IntegrateRationalFunction(f), zero(base_ring(parent(f))), 1
    end
    g1, h, r = HermiteReduce(f, D)
    ss, Ss, ρ = ResidueReduce(h, D)
    g2, Dg2, ss1, Ss1 = ConstantPart(ss, Ss, D)
    if !isempty(ss1)
        throw(AlgebraicNumbersInvolved())
        # From this point on all computations have to be performed in extension fields
        # over field of algebraic (instead of merely rational) numbers, i.e., replace
        # Nemo.QQ by Nemo.QQBar, see http://nemocas.github.io/Nemo.jl/latest/algebraic/
        # For the time being, we start all over again, now with QQBar enabled
        # and discard all results that have been obtained until now.
        # TODO: transform the so far obtained results to fields over QQBar and proceed.
    end
    if ρ<=0
        g = vcat(IdTerm(g1), g2)
        return g, f - D(g1) - Dg2, 0
    end
    p = h - Dg2 + r
    if isprimitive(D)
        @assert isone(denominator(p))
        q, ρ = IntegratePrimitivePolynomial(numerator(p), D)
        g = vcat(IdTerm(g1 + q), g2)
        f1 = f - D(g1 + q) - Dg2
    elseif ishyperexponential(D)
        q, ρ = IntegrateHyperexponentialPolynomial(p, D)
        g = vcat(IdTerm(g1 + q), g2)
        f1 = f - D(g1 + q) - Dg2
    elseif ishypertangent(D)
        q1, ρ = IntegrateHypertangentReduced(p, D)
        if ρ<=0
            g = vcat(IdTerm(g1 + q1), g2)
            f1 = f - D(g1 + q1) - Dg2
        else
            @assert isone(denominator(p - D(q1)))
            q2, c = IntegrateHypertangentPolynomial(numerator(p - D(q1)), D)
            if iszero(BaseDerivation(D)(c))
                t = gen(base_ring(parent(f)))
                η = constant_coefficient(divexact(MonomialDerivative(D), 1 + t^2))
                g = vcat(IdTerm(g1 + q1 + q2), g2, FunctionTerm(log, c, 1 + t^2))
                f1 = f - D(g1 + q1 + q2) - Dg2 - 2*c*η*t
            else
                g = vcat(IdTerm(g1 + q1 + q2), g2)
                f1 = f - D(g1 + q1 + q2) - Dg2
                ρ = 0
            end
        end
    else
        H = MonomialDerivative(D)
        throw(NotImplementedError("Integrate: monomial derivative = $H\n@ $(@__FILE__):$(@__LINE__)"))
    end
    if ρ<=0
        return g, f1, ρ
    end
    @assert isone(denominator(f1)) && degree(numerator(f1)) <= 0
    f1 = constant_coefficient(numerator(f1))
    h, f2,  ρ = Integrate(f1, BaseDerivation(D))
    @assert ρ<=0 || iszero(f2)
    vcat(h, g), f2, ρ
end

"""
    InFieldDerivative(f, D) -> (u, ρ)

Given a field `K`, a derivation `D` on `K` and `f` in `K`, return either `ρ=0`, in which case
the equation `D(u)=f` does not have a solution `u` in `K`, or `ρ=1` and a solution `u` in `K`
of that equation.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 5.12, p. 175.
"""
function InFieldDerivative(f::F, D::Derivation) where F<:FieldElement
    # base case f \in constant field, D must be the null derivation
    iscompatible(f, D) || error("field element f must be in the domain of derivation D")
    iszero(D) || error("base case only for null derivations")
    if iszero(f)
        return zero(f), 1
    else
        return zero(f), 0
    end
end

function InFieldDerivative(f::F, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    # See Bronstein's book, Section 5.12, p. 175
    iscompatible(f, D) || error("rational function f must be in the domain of derivation D")
    Z = zero(f)
    no_solution = (Z, 0)
    g, h, r = HermiteReduce(f, D)
    if !isone(denominator(h))
        # no solution, h returned by HermiteReduce not a polynomial
        return no_solution
    end
    p = h + r # reduced
    @assert isreduced(p, D)
    if isprimitive(D)
        # p reduced => polynomial in this case
        p = numerator(p)
        q, ρ = IntegratePrimitivePolynomial(p, D)
        if ρ<=0
            return no_solution
        end
        a0 = p-D(q)
        @assert degree(a0)<=0 # p-D(q) ∈ k
        a = constant_coefficient(a0)
        v, c, ρ = LimitedIntegrate(a, leading_coefficient(D), BaseDerivation(D))
        if ρ<=0
            return no_solution
        end
        t = gen(parent(p))
        return g + q + v + c*t, 1
    else
        if ishyperexponential(D)
            q, ρ = IntegrateHyperexponentialPolynomial(p, D)
            if ρ<=0
                return no_solution
            end
            a0 = p-D(q)
        elseif ishypertangent(D)
            q1, ρ = IntegrateHypertangentReduced(p, D)
            if ρ<=0
                return no_solution
            end
            p1 = p - D(q1)
            @assert isone(denominator(p1))
            q2, c = IntegrateHypertangentPolynomial(numerator(p1), D)
            if !iszero(c)
                # This is my (maybe too naive) interpretation of "prevent
                # introducing a new logarithm". See parenthetical comment
                # on top of p.176.
                return no_solution
            end
            q = q1 + q2
            a0 = p1 - D(q2)
        else
            H = MonomialDerivative(D)
            throw(NotImplementedError("InFieldDerivative: monomial deivative =$H\n@ $(@__FILE__):$(@__LINE__)"))
        end
        @assert isone(denominator(a0)) && degree(numerator(a0))<=0 # p-D(q) ∈ k
        a = constant_coefficient(numerator(a0))
        v, ρ = InFieldDerivative(a, BaseDerivation(D))
        if ρ<=0
            return no_solution
        end
        return g + q + v, 1
    end
end

"""
    InFieldLogarithmicDerivative(f, D) -> (u, ρ)

Given a field `K`, a derivation `D` on `K` and `f` in `K`, return either `ρ=0`, in which case
the equation `D(u)/u=f` does not have a nonzero  solution `u` in `K`, or `ρ=1` and a nonzero
solution `u` in `K` of that equation.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 5.12, p. 176.
"""
function InFieldLogarithmicDerivative(f::F, D::Derivation) where F<:FieldElement
    n, v, β = InFieldLogarithmicDerivativeOfRadical(f, D, expect_one=true)
    @assert β<=0 || n==1
    v, β
end

"""
    InFieldLogarithmicDerivativeOfRadical(f, D) -> (n, u, ρ)

Given a field `K`, a derivation `D` on `K` and `f` in `K`, return either `ρ=0`, in which case
the equation `D(u)/u=n*f` does not have a solution `u≠0` in `K` and integer `n≠0`, or `ρ=1` and a solution
`u` and `n` of that equation.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 5.12, p. 177.
"""
function InFieldLogarithmicDerivativeOfRadical(f::F, D::Derivation; expect_one::Bool=false) where F<:FieldElement
    # base case f \in constant field, D must be the null derivation
    iscompatible(f, D) || error("field element f must be in the domain of derivation D")
    iszero(D) || error("base case only for null derivations")
    if iszero(f)
        return 1, one(f), 1
    else
        return 0, zero(f), 0 # no solution (n must be !=0)
    end
end

function InFieldLogarithmicDerivativeOfRadical(f::F, D::Derivation; expect_one::Bool=false) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    iscompatible(f, D) || error("rational function f must be in the domain of derivation D")
    if iszero(f)
        return 1, one(f), 1
    end
    Z = zero(f)
    no_solution = (0, Z, 0)
    if !issimple(f, D)
        return no_solution
    end
    if isone(denominator(f))
        m = 1
        Dg = Z
        v = one(f)
    else
        Rs, Gs, ρ = ResidueReduce(f, D)
        @assert ρ>=1
        if length(Rs)==0
            m = 1
            Dg = Z
            v = one(f)
        else
            if !all([all(isrational.(coefficients(R))) for R in Rs])
                # no solution, Rothstein-Trager resultant has non-rational coefficients
                return no_solution
            end
            Rs1 = [map_coefficients(c->convert(fmpq, rationalize(c)), R) for R in Rs]
            As1 = [roots(R) for R in Rs1]
            if !(all([length(As1[i])==degree(Rs1[i]) for i=1:length(As1)]) &&
                all([all(isrational.(As1[i])) for i=1:length(As1)]) )
                # no solution, not all roots of the Rothstein-Trager resultant are rational
                return no_solution
            end
            As = [[rationalize_over_Int(a) for a in as] for as in As1]
            if expect_one
                if !all([all([isone(denominator(a)) for a in as]) for as in As])
                    # no solution, not all roots of the Rothstein-Trager resultant are integers"
                    return no_solution
                end
                m = 1
            else
                m = gcd([gcd(denominator.(as)) for as in As])
            end
            Dg = Z
            v = one(f)
            for i=1:length(As)
                for a in As[i]
                    lc = leading_coefficient(Gs[i])(a) # make monic (!!!)
                    g = map_coefficients(c->c(a)//lc, Gs[i])
                    Dg += a*D(g)//g
                    # Make g = g+Z an element of the field k(t),
                    # otherwise exponentiation with negative numbers would not work.
                    v *= (g+Z)^numerator(m*a)
                end
            end
        end
    end
    p = f - Dg
    if iszero(p)
        if expect_one && m!=1
            # no solution, m not == 1 as expected"
            return no_solution
        end
        return m, v, 1
    end
    if !(isone(denominator(p)))
        # no solution, p=f-Dg not a polynomial
        return no_solution
    end
    if !(degree(numerator(p))<max(1,degree(D)))
        return no_solution
    end
    H = MonomialDerivative(D)
    if isprimitive(D)
        p0 = constant_coefficient(numerator(p))
        n, u, ρ = InFieldLogarithmicDerivativeOfRadical(p0, BaseDerivation(D))
        if ρ<=0
            return no_solution
        end
        if expect_one && n!=1
            # no solution, n not == 1 as expected
            return no_solution
        end
        N = lcm(n, m)
        # Make u = u + Z an element of the field k(t),
        # otherwise exponentiation with negative numbers would not work.
        U = v^div(N, m)*(u+Z)^div(N, n)
        return N, U, 1
    elseif ishyperexponential(D)
        p0 = constant_coefficient(numerator(p))
        w = coeff(H, 1)
        n, e, u, ρ = ParametricLogarithmicDerivative(p0, w, BaseDerivation(D))
        if ρ<=0
            return no_solution
        end
        if expect_one && n!=1
            # no solution, m not == 1 as expected"
            return no_solution
        end
        N = lcm(n, m)
        t = gen(parent(numerator(f)))
        # Make u = u + Z, t = t + Z elements of the field k(t),
        # otherwise exponentiation with negative numbers would not work.
        U = v^div(N, m)*(u+Z)^div(N, n)*(t+Z)^div(e*N, n)
        return N, U, 1
    elseif ishypertangent(D)
        p0 = numerator(p)
        @assert degree(p0)<=1
        a = coeff(p0, 0)
        b = coeff(p0, 1)
        t = gen(parent(H))
        η = constant_coefficient(divexact(H, t^2+1))
        if !contains_I(parent(f))
            c = (1//2)*b//η
            if !isrational(c)
                return no_solution
            end
            n, u, ρ = InFieldLogarithmicDerivativeOfRadical(a, BaseDerivation(D))
            if ρ<=0
                return no_solution
            end
            e0 = n*rationalize_over_Int(c)
            if !isone(denominator(e0))
                return no_solution
            end
            e = numerator(e0)
            N = lcm(n, m)
            U = v^div(N, m)*(u+Z)^div(N, n)*(t^2+1+Z)^div(e*N, n)
            return N, U, 1
        elseif iszero(real(a))
            # Note: This case is not treated in 5.12 of Bronsteins's book, altough it seems
            # to be the one relevant for checking the condition of Theorem 5.10.1.
            # I did not prove that in the case of Theorem 5.12 real(a)=0 always holds true.
            ai = imag(a)
            c1 = (1//2)*(b+ai)//η
            c2 = (1//2)*(b-ai)//η
            if !isrational(c1) || !isrational(c2)
                return no_solution
            end
            # implicitely set u = 1 => D(u)//u = 0
            c1 = rationalize_over_Int(c1)
            c2 = rationalize_over_Int(c2)
            n = lcm(denominator(c1), denominator(c2))
            e1 = numerator(c1*n)
            e2 = numerator(c2*n)
            N = lcm(n, m)
            I = get_I(parent(Z))
            U = v^div(N, m)*(t-I)^div(e1*N, n)*(t+I)^div(e2*N, n)
            return N, U, 1
        else
            # TODO...
            throw(NotImplementedError("InFieldLogarithmicDerivativeOfRadical: hypertangent case with I in K and real(a)!=0\n@ $(@__FILE__):$(@__LINE__)"))
        end
    else
        throw(NotImplementedError("InFieldLogarithmicDerivativeOfRadical: monomial derivative $H\n@ $(@__FILE__):$(@__LINE__)"))
    end
end


text:
############################ STOKES ############################3
struct Stokes{N,T} <: AbstractPDE{N}
    μ::T
end
Stokes(;μ,dim=3)             = Stokes{dim}(μ)
Stokes{N}(μ::T) where {N,T}  = Stokes{N,T}(μ)

getname(::Stokes) = "Stokes"

default_kernel_eltype(::Stokes{N}) where {N} = Mat{N,N,Float64,N*N}
default_density_eltype(::Stokes{N}) where {N} = Vec{N,Float64}

# Single Layer
function (SL::SingleLayerKernel{T,S})(x,y)::T  where {T,S<:Stokes}
    N = ambient_dim(S)
    x==y && return zero(T)
    μ = SL.op.μ
    r = x .- y
    d = norm(r)
    RRT = r*transpose(r)
    if N==2
        ID = Mat{2,2,Float64,4}(1,0,0,1)
        γ  = -log(d)
    elseif N==3
        ID = Mat3{Float64}(1,0,0,0,1,0,0,0,1)
        γ  = 1/d
   end
    return 1/(4π*(N-1)*μ) * (γ*ID + RRT/d^N)
end

# Double Layer Kernel
function (DL::DoubleLayerKernel{T,S})(x,y,ny)::T where {T,S<:Stokes}
    N = ambient_dim(S)
    μ = DL.op.μ
    r = x .- y
    d = norm(r)
    if N==2
        x==y && return zero(T)
        return 1/π * dot(r,ny)/d^4 .* r*transpose(r)
    elseif N==3
        x==y && return zero(T)
        ID = Mat3{Float64}(1,0,0,0,1,0,0,0,1)
        RRT = r*transpose(r) # r ⊗ rᵗ
        return 3/(4π) * dot(r,ny)/d^5 .* RRT
    end
end

# TODO: code the hypersingular and adjoint double layer operators for Stokes in 2 and 3d
# Adjoint Double Layer Kernel
# function (HS::AdjointDoubleLayerKernel{N,T,Op})(x,y,ny)::T where {N,T,Op<:Stokes}
#     μ = HS.op.μ
#     r = x .- y
#     d = norm(r)
#     if N==2
#         x==y && return zero(T)
#         return -1/π * dot(r,ny)/d^4 .* r*transpose(r)
#     elseif N==3
#         x==y && return zero(T)
#         ID = Mat3{Float64}(1,0,0,0,1,0,0,0,1)
#         RRT = r*transpose(r) # r ⊗ rᵗ
#         return -3/(4π) * dot(r,ny)/d^5 .* RRT
#     end
# end
# AdjointDoubleLayerKernel{N}(op::Op,args...) where {N,Op<:Stokes} = AdjointDoubleLayerKernel{N,Mat{N,N,Float64,N*N},Op}(op,args...)

# # Hypersingular
# function (DL::HyperSingularKernel{N,T,Op})(x,y,nx,ny)::T where {N,T,Op<:Stokes}
#     μ = DL.op.μ
#     r = x .- y
#     d = norm(r)
#     if N==2
#         x==y && return zero(T)
#         return -1/π * dot(r,ny)/d^4 .* r*transpose(r)
#     elseif N==3
#         x==y && return zero(T)
#         ID = Mat3{Float64}(1,0,0,0,1,0,0,0,1)
#         RRT = r*transpose(r) # r ⊗ rᵗ
#         return -3/(4π) * dot(r,ny)/d^5 .* RRT
#     end
# end
# HyperSingularKernel{N}(op::Op,args...) where {N,Op<:Stokes} = HyperSingularKernel{N,Mat{N,N,Float64,N*N},Op}(op,args...)

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# SPDX-License-Identifier: MIT

#======================================================================
This file wrapped a few basic linear algebra operations in order to
have more fine-grained control between CPU/GPU code and different
semifield algebra operations.
======================================================================#

# Element-wise division.
eldiv!(out, x::AbstractArray, y::AbstractArray) = broadcast!(/, out, x, y)

# Element-wise multiplication.
elmul!(out, x::AbstractArray, y::AbstractArray) =  broadcast!(*, out, x, y)
function elmul!(out, x::AbstractSparseArray{T}, y::AbstractArray{T}) where T
    fill!(out, zero(T))
    broadcast!(*, out, x, y)
end
elmul!(out, x::CuSparseVector, y::AbstractVector) = elmul_svdv!(out, x, y)
elmul!(out, x::CuSparseVector, Y::AbstractMatrix) = elmul_svdm!(out, x, Y)

# Matrix multiplication.
matmul!(out, X::AbstractMatrix{T}, y::AbstractArray{T}) where T =
    mul!(out, X, y, one(T), zero(T))
function matmul!(out, X::AbstractSparseMatrix{T}, y::AbstractArray{T}) where T
    fill!(out, zero(T))
    mul!(out, X, y, one(T), zero(T))
end
matmul!(out, X::CuMatrix{T}, Y::CuMatrix{T}) where T =
    mul_dmdm!(out, X, Y)
matmul!(out, X::CuSparseMatrixCSR{T}, y::AbstractVector{T}) where T =
    mul_smdv!(out, X, y)
matmul!(out, X::CuSparseMatrixCSR{T}, Y::AbstractMatrix{T}) where T =
    mul_smdm!(out, X, Y)

max_stars_repo_name: MartinKocour/MarkovModels

max_issues_repo_name: MartinKocour/MarkovModels

max_forks_repo_name: MartinKocour/MarkovModels


text:
centered{N,T}(C::Type{HyperCube{N,T}}) = C(Vec{N,T}(-0.5), T(1))
centered{T<:HyperCube}(::Type{T}) = centered(HyperCube{ndims_or(T, 3), eltype_or(T, Float32)})

centered{N,T}(R::Type{HyperRectangle{N,T}}) = R(Vec{N,T}(-0.5), Vec{N,T}(1))
centered{T<:HyperRectangle}(::Type{T}) = centered(HyperRectangle{ndims_or(T, 3), eltype_or(T, Float32)})

centered{N,T}(S::Type{HyperSphere{N,T}}) = S(Vec{N,T}(0), T(0.5))
centered{T<:HyperSphere}(::Type{T}) = centered(HyperSphere{ndims_or(T, 3), eltype_or(T, Float32)})


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

function selection_sort!(A)
    for i in 1:(length(A)-1)
        min_index = i

        for j in i+1:length(A)
            if A[j] < A[min_index]
                min_index = j
            end
        end

        exchange!(A, i, min_index)
    end
end

max_stars_repo_name: blallen/CodePractice

max_issues_repo_name: blallen/CodePractice

max_forks_repo_name: blallen/CodePractice


text:
#   Author: V. Vitale
#   Feb 2022

include("mps.jl")
include("mpo.jl")

function save_rhoi(A::MPS)
    right_normalize!(A)
    p=zeros(A.N)
    @tensor temp[:]:=A.data[1][1,-1,2]*conj(A.data[1][1,-2,2])
    p[1]=real(temp[1,1])
    for i in 1:A.N
        sA = size(A.data[i])
        U,S,V = svd(reshape(A.data[i],(sA[1]*sA[2],sA[3])),full=false)
        V=V'
        A.data[i] = reshape( U,( sA[1], sA[2], :))
        if i<A.N
            S=diagm(S)
            @tensor A.data[i+1][:] := S[-1,1 ] * V[ 1,2 ] * A.data[i+1][2,-2,-3]
            @tensor temp[:]:=A.data[i+1][1,-1,2]*conj(A.data[i+1][1,-2,2])
            p[i+1]=real(temp[1,1])
        else
            @tensor temp[:]:=A.data[i][1,-1,2]*conj(A.data[i][1,-2,2])
            p[i]=real(temp[1,1])
        end

    end
    return p
end

function sample_shadows(A::MPS,nu::Int,nm::Int,lA::Int)
    if nm==1
        shadows=[]
        ket=[[1 ; 0],[0 ; 1]]
        for r in 1:nu
            WlA=MPO()
            Initialize!("Local_Haar",WlA,lA)

            if lA<A.N
                W=MPO()
                Initialize!("Id",W,A.N)
                for i in 1:lA
                    W.data[i]=Base.copy(WlA.data[i])
                end
            else
                W=copy(WlA)
            end

            state=W*A
            p_state=save_rhoi(state)[1:lA]
            for m in 1:nm
                shadow=MPO()
                Initialize!("Zeros",shadow,2,1,lA)
                k=zeros(Int,lA)
                p_extract=rand(lA)
                k[p_extract.>p_state].=1
                for (qubit,spin) in enumerate(k)
                    shadow.data[qubit][1,1,:,:]=3*kron(ket[spin+1],ket[spin+1]')-1.0*I(2)
                end
                push!(shadows,(WlA*shadow)*adjoint(WlA))
            end
        end
    else
        shadows=[]
        ket=[[1 ; 0],[0 ; 1]]
        for r in 1:nu
            WlA=MPO()
            Initialize!("Local_Haar",WlA,lA)

            if lA<A.N
                W=MPO()
                Initialize!("Id",W,A.N)
                for i in 1:lA
                    W.data[i]=Base.copy(WlA.data[i])
                end
            else
                W=copy(WlA)
            end

            state=W*A
            p_state=save_rhoi(state)[1:lA]
            for m in 1:nm
                shadow=MPO()
                Initialize!("Zeros",shadow,2,1,lA)
                k=zeros(Int,lA)
                p_extract=rand(lA)
                k[p_extract.>p_state].=1
                for (qubit,spin) in enumerate(k)
                    shadow.data[qubit][1,1,:,:]=3*kron(ket[spin+1],ket[spin+1]')-1.0*I(2)
                end
                push!(shadows,(WlA*shadow)*adjoint(WlA))
            end
        end
    end
    return shadows
end

function calculate_shadows_MPO(nu::Int,nm::Int,lU::Array,lstr::Array,qubit_set::Array)
    ket=[[1 ; 0],[0 ; 1]]
    chi=1
    d=2
    rho=Dict()
    for r in 1:nu
        counts=lstr[r]
        rho[r]=MPO()
        Initialize!("Zeros",rho[r],d,chi,length(qubit_set))
        for (step_m,m) in enumerate(keys(counts))
            str=string(m-1,base=2,pad=length(qubit_set))
            for (step,j) in enumerate(qubit_set)
                kj = parse(Int32,str[j])+1
                temp=Array{Complex{Float64}}(undef,chi,chi,d,d)
                temp[chi,chi,:,:]=3*(lU[r][j]')*kron(ket[kj],ket[kj]')*lU[r][j]-1.0*I(2)

                if step_m==1
                    rho[r].data[step]=temp*((counts[m]/nm)^(1/length(qubit_set)))
                else
                    rho[r].data[step]=++(rho[r].data[step],temp*((counts[m]/nm)^(1/length(qubit_set))))
                end
            end
        end
    end
    return rho
end

function calc_dist(rhos::Array,nu::Int,nm::Int)

    distrhos=zeros(nu*nm,nu*nm)
    idx=zeros(Int,nu*nm,nu*nm)

    for i in 1:nu*nm
        for j in i:nu*nm
            distrhos[i,j]=real(tr(rhos[i]*adjoint(rhos[i]))-
                tr(rhos[i]*adjoint(rhos[j]))-
                tr(rhos[j]*adjoint(rhos[i]))+
                tr(rhos[j]*adjoint(rhos[j])))^(1/2)
        end
    end
    distrhos=distrhos+distrhos'
    for i in 1:length(rhos)
        idx[i,:]=sortperm(distrhos[i,:])
        distrhos[i,:]=distrhos[i,idx[i,:]]
    end

    return distrhos
end

max_stars_repo_name: aanniikkeettr/Tensor-Networks

max_issues_repo_name: aanniikkeettr/Tensor-Networks

max_forks_repo_name: aanniikkeettr/Tensor-Networks


text:
function _uview_destride(vec::Strided.StridedView{T,3}, i) where {T}
    stride_last = vec.strides[end]

    start = vec.offset + (i-1) * stride_last
    len   = prod(vec.size[1:end-1])

    return reshape(uview(vec.parent, start:start+len), vec.size[1], vec.size[2])
end

function _uview_vec_destride(vec::Strided.StridedView{T,3}, i) where {T}
    stride_last = vec.strides[end]

    start = vec.offset + (i-1) * stride_last
    len   = prod(vec.size[1:end-1])

    return uview(vec.parent, start+1:start+len)
end


text:
################ Define Asymmetric Symbolic Error-based Intervals ##############

struct AsymESIP{F<:AbstractPolytope}
    equation
    factors
    errors
    domain::F
end

function init_asym_esip(input::AbstractHyperrectangle)
    n = dim(input)
    equation = [I zeros(n)]
    factors = zeros(n, 0)
    errors = zeros(0, n + 1)  # extra dimension for constants
    return AsymESIP(equation, factors, errors, input)
end

function sym_bounds_matrix(a::AsymESIP)
    eq = a.equation
    sym_lo = eq + min.(0, a.factors) * a.errors
    sym_up = eq + max.(0, a.factors) * a.errors
    return sym_lo, sym_up
end

function bounds_matrix(a::AsymESIP)
    sym_lo, sym_up = sym_bounds_matrix(a)

    los = low(a.domain)
    his = high(a.domain)

    lb = max.(0, sym_lo[:, 1:end-1]) * los + min.(0, sym_lo[:, 1:end-1]) * his + sym_lo[:, end]
    ub = max.(0, sym_up[:, 1:end-1]) * his + min.(0, sym_up[:, 1:end-1]) * los + sym_up[:, end]

    return lb, ub
end


text:
module AutomotiveSensors

using AutomotiveSimulator
using AutomotiveVisualization
using Parameters
using StaticArrays
using Random
using LinearAlgebra
using Distributions

export
    AbstractSensor,
    measure

abstract type AbstractSensor end
abstract type AbstractNoiseModel end

export
    occlusion_checker,
    OcclusionOverlay

include("occlusion_checker.jl")

export
    AbstractNoiseModel,
    LinearNoise

include("noise_model.jl")

export
    GaussianSensor,
    GaussianSensorOverlay

include("gaussian_sensor.jl")
include("false_positive.jl")

export
    PerfectSensor

include("perfect_sensor.jl")

end