text:
struct CPU1Backend{R<:CPU1} <: CPUBackend{R}
    resource::R
end

function backend_init(resource::CPU1, domain, sim)
    CPU1Backend(resource)
end

function backend_update!(state::State, backend::CPU1Backend, sim)
    cpu_kernel!(state, backend.resource, sim)
end

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

struct CPUThreadsBackend{R<:CPUThreads,T} <: CPUBackend{R}
    resource::R
    tiles::T
end

function backend_init(r::CPUThreads{TileSize{N}}, domain::Domain, sim) where N
    dims = gridsize(domain, sim.wave.dx)
    inds = map(s->2:s-1, dims)
    tiles = collect(TileIterator(inds, r.settings.dims))
    CPUThreadsBackend(r, tiles)
end

function backend_init(r::CPUThreads{NTuple{N,Int}}, domain::Domain, sim) where N
    backend_init(CPUThreads(TileSize(r.settings)), domain, sim)
end

function backend_cleanup!(backend::CPUBackend, state, sim)
    nothing
end

move_backend(A::AbstractArray, backend::CPUBackend) = A

function backend_update!(state::State, backend::CPUThreadsBackend, sim)
    tiles = backend.tiles
    Threads.@threads for i = 1:length(tiles)
        inds = tiles[i]
        cpu_kernel!(state, CPUThreads(inds), sim)
    end
end

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

function cpu_kernel_impl(N, indices=:(current_domain))
    quote
        @inbounds @nloops $N I $indices begin
            @nexprs $N w->(
                q₊_w = @nref($N, q, i -> i==w ? I_i+1 : I_i);
                q₋_w = @nref($N, q, i -> i==w ? I_i-1 : I_i)
            )
            # q₊_1 + q₋_1 + q₊_2 + ...
            K = @ncall($N, +, w->(q₊_w + q₋_w))
            # K̄ = $(2*N) - K

            @nexprs $N w->(
                Q_w = q₊_w * @nref($N, Ψₜ, i -> i==w ? I_i+1 : I_i) +
                      q₋_w * @nref($N, Ψₜ, i -> i==w ? I_i-1 : I_i)
            )
            # Q_1 + Q_2 + ...
            Q = @ncall($N, +, Q)

            @nexprs $N w->(γ_w = (1-q₊_w) * γ + (1-q₋_w) * γ)
            # γ_1 + γ_2 + ...
            γ̄ = @ncall($N, +, γ)

            @nref($N, Ψₜ₊₁, I) = 1/(1+λhalf*γ̄) * (
                (2-λsq*K)   * @nref($N, Ψₜ,  I)  +
                (λhalf*γ̄-1) * @nref($N, Ψₜ₋₁,I) +
                λsq*Q)
        end
    end
end

@generated function cpu_kernel!(
        state::BoxState{N,T},
        resource::Union{CPU1,CPUThreads},
        sim::Simulator{N,T,W}) where {N, T, W <: UniformWave{N,T}}
    SETUP  = resource<:CPUThreads ? :(inds = resource.settings) : :()
    REGION = resource<:CPUThreads ? :(i->inds[i]) : :(i->2:size(Ψₜ,i)-1)
    quote
        λ = sim.wave.λ
        λsq = λ^2
        λhalf = λ/2
        γ = T(state.box.γ)
        Ψₜ₋₁ = state.previous
        Ψₜ   = state.current
        Ψₜ₊₁ = state.previous
        q    = state.q
        $(SETUP)
        $(cpu_kernel_impl(N, REGION))
        state
    end
end


text:
using Test
using DataFrames
using Distributions
using LinearAlgebra
using Statistics
using StatsBase

import Random: seed!
import Baysor as B

module DataWrappers

using DataFrames
import Baysor as B

function get_bmm_data(; n_mols::Int=5000, confidence::Bool=false, scale::Float64=0.1, min_molecules_per_cell::Int=10, do_maximize::Bool=false, kwargs...)
    df = DataFrame(:x => rand(n_mols), :y => rand(n_mols), :gene => rand(1:10, n_mols), :confidence => confidence ? ones(n_mols) : rand(n_mols))
    bmd = B.initialize_bmm_data(df, scale=scale, min_molecules_per_cell=min_molecules_per_cell, confidence_nn_id=0, kwargs...);
    if do_maximize
        B.maximize!(bmd)
    end

    return bmd
end

end

@testset "Baysor" begin
    @testset "utils" begin
        @testset "convex_hull" begin
            chull = B.convex_hull([[0, 0], [1, 0], [0, 0], [1, 1], [0, 1]])
            @test all(chull .== [0 0 1 1 0; 0 1 1 0 0])
            @test B.area(chull) ≈ 1.0

            chull = B.convex_hull([[0, 0], [1, 0], [0, 0], [1, 1], [0, 1], [0, 2], [2, 0]])
            @test all(chull .== [0 0 2 0; 0 2 0 0])
            @test B.area(chull) ≈ 2.0

            chull = B.convex_hull([[0, 0], [1, 0], [0, 0], [1, 1], [0, 1], [0, 2], [2, 0], [2, 2]])
            @test all(chull .== [0 0 2 2 0; 0 2 2 0 0])
            @test B.area(chull) ≈ 4.0
        end
    end

    @testset "distributions" begin
        @testset "ShapePrior" begin
            seed!(42)

            means = B.MeanVec{2}([10.0, 20.0])
            std_stds = B.MeanVec{2}([5.0, 10.0])
            stds = hcat([Vector(B.sample_var(B.ShapePrior{2}(means, std_stds, 1000))) for i in 1:100000]...) .^ 0.5
            @test all(stds .> 0)
            @test all(abs.(vec(mean(stds, dims=2)) .- means) .< 1)

            means = B.MeanVec{2}([1000.0, 2000.0])
            stds = hcat([Vector(B.sample_var(B.ShapePrior{2}(means, std_stds, 1000))) for i in 1:100000]...) .^ 0.5
            @test all(stds .> 0)
            @test all(abs.(vec(mean(stds, dims=2)) .- means) .< 1)
            @test all(abs.(vec(std(stds, dims=2)) .- std_stds) .< 1)
        end

        @testset "CategoricalSmoothed" begin
            dist = B.CategoricalSmoothed([0.0, 0.0, 0.0])
            B.maximize!(dist, [missing, 1, 2], [0.5, 0.5, 0.5])
            @test all(dist.counts .≈ [0.5, 0.5, 0.0])
            @test dist.sum_counts ≈ 1.0
        end
    end

    @testset "data_processing" begin
        @testset "initialization" begin
            n_mols = 1000
            df = DataFrame(:x => rand(n_mols), :y => rand(n_mols), :gene => rand(1:10, n_mols))

            df[!, :cluster] = rand(1:5, n_mols)
            df[!, :prior_segmentation] = rand(1:100, n_mols)
            bm_data = B.initialize_bmm_data(df, scale=6.0, min_molecules_per_cell=30);
            @test all(bm_data.cluster_per_molecule .== df.cluster)
            @test all(bm_data.segment_per_molecule .== df.prior_segmentation)

            @test all(length(bm_data.adjacent_points[i]) == length(bm_data.adjacent_weights[i]) for i in 1:length(bm_data.adjacent_points))
            @test !any(i in bm_data.adjacent_points[i] for i in 1:length(bm_data.adjacent_points))

            for i in 5:10:55
                init_params = B.cell_centers_uniformly(df, i; scale=10.0)
                @test size(init_params.centers, 1) == i
                @test size(init_params.centers, 2) == 2
                @test length(init_params.covs) == i
                @test maximum(init_params.assignment) == i
            end
        end

        @testset "parse_parameters" begin
            @test B.parse_scale_std("23.55%", 121.0) ≈ 121.0 * 0.2355
            @test B.parse_scale_std(33.87, 121.0) ≈ 33.87
            @test B.parse_scale_std(nothing, 121.0) ≈ 121.0 * 0.25
        end

        @testset "molecule_graph" begin
            for adj_type in [:triangulation, :knn, :both]
                adj_points, adj_weights = B.build_molecule_graph(DataFrame(rand(1000, 2), [:x, :y]); adjacency_type=adj_type, k_adj=30)[1:2];
                @test all(length.(adj_points) .== length.(adj_weights))
                @test all(length.(adj_points) .== length.(unique.(adj_points)))
            end
        end

        @testset "filter_segmentation" begin
            labels = [1 2 3; 1 2 4; 5 2 7];
            position_data = hcat(vec([[id[1], id[2]] for id in CartesianIndices(labels)])...)
            df_spatial = DataFrame(:x => position_data[1,:], :y => position_data[2,:])
            assignment = B.staining_value_per_transcript(df_spatial, labels)

            expected_results = [
                labels,
                labels,
                [1 2 0; 1 2 0; 0 2 0],
                [0 2 0; 0 2 0; 0 2 0],
                [0 0 0; 0 0 0; 0 0 0]
            ]
            for m in 0:4
                l, ta = B.filter_segmentation_labels!(deepcopy(labels), deepcopy(assignment), min_molecules_per_segment=m)
                @test all(l .== expected_results[m + 1])

                a = B.filter_segmentation_labels!(deepcopy(assignment), min_molecules_per_segment=m)
                @test all(a .== B.staining_value_per_transcript(df_spatial, l))
                @test all(a .== ta)
            end
        end
    end

    @testset "bmm_algorithm" begin
        @testset "noise_composition_density" begin
            for conf in [false, true]
                bm_data = DataWrappers.get_bmm_data(confidence=conf, do_maximize=true);
                dens_exp = mean([mean(c.composition_params.counts[c.composition_params.counts .> 0] ./ c.composition_params.sum_counts) for c in bm_data.components]);
                dens_obs = B.noise_composition_density(bm_data)
                @test abs(dens_exp - dens_obs) < 1e-10
            end
        end

        @testset "distribution_sampling" begin
            bm_data = DataWrappers.get_bmm_data(confidence=true, do_maximize=true);
            @test_nowarn for i in 1:1000
                B.sample_distribution!(bm_data)
            end
        end

        @testset "drop_unused_components" begin
            bm_data = DataWrappers.get_bmm_data(n_mols=1000, confidence=true);
            bm_data.assignment[rand(1:length(bm_data.assignment), 100)] .= 0
            B.maximize!(bm_data)
            B.drop_unused_components!(bm_data)

            @test all(B.num_of_molecules_per_cell(bm_data) .== [c.n_samples for c in bm_data.components])
        end

        @testset "maximize_mvnormal" begin
            is_not_nan = map(1:10000) do i
                n_samps = rand(0:100)
                d = B.MvNormalF([0., 0.], diagm(0 => ones(2)));
                B.maximize!(d, rand(2, n_samps); center_probs=rand(n_samps) .* rand(Binomial(1, 0.5), n_samps))
                !isnan(d.Σ[1])
            end

            @test all(is_not_nan)
        end

        @testset "synthetic_run" begin
            seed!(42)
            for it in 1:6
                n_components = 10;
                n_genes = 20
                scale = 0.2;
                frame_size = (100, 200)
                centers = hcat(rand(n_components) * frame_size[1], rand(n_components) * frame_size[2]);
                sizes = rand(Poisson(50), n_components);
                expressions = rand(Dirichlet(n_genes, 0.1), n_components);

                df_spatial = vcat([DataFrame(
                    :x => rand(Normal(centers[i,1], scale), sizes[i]),
                    :y => rand(Normal(centers[i,2], scale), sizes[i]),
                    :gene => rand(Categorical(expressions[:,i]), sizes[i]),
                    :cell => i
                ) for i in 1:n_components]...)

                noise_size = 100
                df_spatial = vcat(df_spatial, DataFrame(:x => rand(noise_size) * frame_size[1], :y => rand(noise_size) * frame_size[2],
                    :gene => rand(1:n_genes, noise_size), :cell => 0));

                if it % 2 == 0
                    df_spatial[!, :z] = rand(1:3, size(df_spatial, 1)) .* scale ./ 10
                end

                if it > 4
                    df_spatial[!, :cluster] = df_spatial.cell .+ 1
                end
                bm_data_arr = B.initialize_bmm_data(df_spatial; scale=scale, scale_std="5%", min_molecules_per_cell=10);
                bm_data = B.bmm!(bm_data_arr; n_iters=300, new_component_frac=0.3, min_molecules_per_cell=10, assignment_history_depth=50, verbose=false);

                conj_table = counts(B.estimate_assignment_by_history(bm_data)[1], df_spatial.cell)[2:end, 2:end];
                @test all([all(vec(mapslices(maximum, conj_table ./ sum(conj_table, dims=d), dims=d)) .> 0.5) for d in 1:2])

                # Test that it works
                seg_df = B.get_segmentation_df(bm_data)
                gene_names = ["g$i" for i in 1:maximum(df_spatial.gene)]
                seg_df = B.get_segmentation_df(bm_data, gene_names)
                df_res = B.get_cell_stat_df(bm_data, seg_df)
                cm1 = B.convert_segmentation_to_counts(B.composition_data(bm_data), bm_data.assignment)
                cm2 = B.convert_segmentation_to_counts(B.composition_data(bm_data), bm_data.assignment; gene_names=gene_names)
                @test all((size(cm2[:, 2:end]) .== size(cm1)))
            end
        end
    end
end

ext: jl

max_stars_repo_name: repushko/Baysor

max_issues_repo_name: repushko/Baysor

max_forks_repo_name: repushko/Baysor


text: Dict{String,Any}("best-day-ever" => Dict{String,Any}("value" => "1987-07-05T17:45:00Z","type" => "datetime"),"numtheory" => Dict{String,Any}("perfection" => Dict{String,Any}("value" => Any[Dict{String,Any}("value" => "6","type" => "integer"), Dict{String,Any}("value" => "28","type" => "integer"), Dict{String,Any}("value" => "496","type" => "integer")],"type" => "array"),"boring" => Dict{String,Any}("value" => "false","type" => "bool")))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
# Support for examining mingw-style pseudo relocations

@struct immutable runtime_pseudo_reloc_item_v1
  addend::UInt32
  target::UInt32
end

@struct immutable runtime_pseudo_reloc_item_v2
  sym::UInt32
  target::UInt32
  flags::UInt32
end

@struct immutable runtime_pseudo_reloc_v2
  magic1::UInt32
  magic2::UInt32
  version::UInt32
end

let findsymbol(h, name) =
    first(filter(x->ObjFileBase.symname(x)==name, COFF.Symbols(h)))

  global dump_pseudo_relocs
  function dump_pseudo_relocs(h)
    start = findsymbol(h, "__RUNTIME_PSEUDO_RELOC_LIST__")
    last = findsymbol(h, "__RUNTIME_PSEUDO_RELOC_LIST_END__")
    seekstart(start)
    startpos = position(h)
    magic_struct = unpack(h, runtime_pseudo_reloc_v2)
    @assert (magic_struct.magic1 == 0 && magic_struct.magic2 == 0 &&
      magic_struct.version == 1) || error("Unsupported pseudo relocation format")
    while position(h) < startpos + (last.entry.Value - start.entry.Value)
      item = unpack(h, runtime_pseudo_reloc_item_v2)
      @show item
    end
  end

end


text:
#################################################################
#
#    Definition of Likelihood-type model
#
#   (Basic MCMC model type based on evaluating the log-target)
#
#   Examples of other possible models: MCMCHierarchicalModel,
#      MCMCGPModel, MCMCKernelModel
#
#################################################################

export MCMCLikModel

println("Loading MCMCLikModel model")

# The likelihood Model type
type MCMCLikelihoodModel <: MCMCModel
	eval::Function              # log-likelihood evaluation function
	evalg::FunctionOrNothing 			# gradient vector evaluation function
	evalt::FunctionOrNothing 			# tensor evaluation function
	evaldt::FunctionOrNothing 			# tensor derivative evaluation function
  evalallg::FunctionOrNothing 		# 2-tuple (log-lik, gradient vector) evaluation function
  evalallt::FunctionOrNothing 		# 3-tuple (log-lik, gradient vector, tensor) evaluation function
  evalalldt::FunctionOrNothing 		# 4-tuple (log-lik, gradient vector, tensor, tensor derivative) evaluation function
	pmap::Dict                  # map to/from parameter vector from/to user-friendly variables
	size::Int                   # parameter vector size
	init::Vector{Float64}       # parameter vector initial values
	scale::Vector{Float64}      # scaling hint on parameters

	MCMCLikelihoodModel(f::Function,
						g::FunctionOrNothing, ag::FunctionOrNothing,
						t::FunctionOrNothing, at::FunctionOrNothing,
						dt::FunctionOrNothing, adt::FunctionOrNothing,
						init::Vector{Float64},
						sc::Vector{Float64},
						pmap::Dict) = begin

		s = size(init, 1)

		@assert ispartition(pmap, s) "param map is not a partition of parameter vector"
		@assert size(sc,1) == s "scale parameter size ($(size(sc,1))) different from initial values ($s)"

    fin = [f, g, ag, t, at, dt, adt]
    fout = Array(FunctionOrNothing, 7)

		# check that all functions can be called with a vector of Float64 as argument
		for i = 1:7
			if isgeneric(fin[i]) && !method_exists(fin[i], (Vector{Float64},))
				if method_exists(fin[i], (Float64,))
					fout[i] = x::Vector{Float64} -> fin[i](x[1])
				else
				  error("one of the supplied functions cannot be called with Vector{Float64}")
          #TODO : make error message print which function is problematic
        end
			else
			    fout[i] = fin[i]
			end
		end

		# check that initial values are in the support of likelihood function
		@assert isfinite(fout[1](init)) "Initial values out of model support, try other values"

    new(fout[1], fout[2], fout[4], fout[6], fout[3], fout[5], fout[7], pmap, s, init, sc)
	end
end

function show(io::IO, res::MCMCLikelihoodModel)
  print(io, "LikelihoodModel, with $(length(res.pmap)) parameter(s)")
  hasgradient(res) && print(io, ", with gradient")
  hastensor(res) && print(io, "/tensor")
  hasdtensor(res) && print(io, "/dtensor")
  println(io)
end

typealias MCMCLikModel MCMCLikelihoodModel

# Model creation using expression parsing and autodiff
function MCMCLikelihoodModel(	m::Expr;
								gradient::Bool=false,
								init=nothing,
								pmap=nothing,
								scale::F64OrVectorF64 = 1.0,
								args...)
	# when using expressions, initial values are passed in keyword args
	#  with one arg by parameter, therefore there is not need for an init arg
	@assert init == nothing "'init' kwargs not allowed for model as expression\n"

	# same thing with 'pmap'
	@assert pmap == nothing "'pmap' kwargs not allowed for model as expression\n"

	# generate lik function
	f, s, p, i = generateModelFunction(m; gradient=false, args...) # loglik only function

	# generate gradient function if requested
	if gradient
		g, s, p, i = generateModelFunction(m; gradient=true, args...) # loglik and gradient function
	else
		g = nothing
	end

	MCMCLikelihoodModel(f, allgrad=g, init=i, pmap=p, scale=scale)
end

# Model creation : with user supplied functions
function MCMCLikelihoodModel(	lik::Function;
								grad::FunctionOrNothing = nothing,
								tensor::FunctionOrNothing = nothing,
								dtensor::FunctionOrNothing = nothing,
								allgrad::FunctionOrNothing = nothing,
								alltensor::FunctionOrNothing = nothing,
								alldtensor::FunctionOrNothing = nothing,
								init::F64OrVectorF64 = [1.0],
								scale::F64OrVectorF64 = 1.0,
								pmap::Union(Nothing, Dict) = nothing)

	# convert init to vector if needed
	init = isa(init, Float64) ? [init] : init

	# expand scale to parameter vector size if needed
	scale = isa(scale, Float64) ? scale * ones(length(init)) : scale

	# all parameters named "pars" by default
	if pmap == nothing ; pmap = Dict([:pars], [(1, size(init))]) ; end

	# now build missing functions, if any
	fmat = Array(FunctionOrNothing, 3, 2)
	for (i, f1, f2) in [(1, grad, allgrad),
						(2, tensor, alltensor),
						(3, dtensor, alldtensor)]
		fmat[i,:] = [f1 f2]
		if f1==nothing && f2!=nothing # only the tuple version is supplied
			fmat[i,1] = (v) -> f2(v)[end]
		elseif f1!=nothing && f2==nothing # only the single version is supplied
			if i == 1
				fmat[i,2] = (v) -> (lik(v), f1(v))
			else
				@assert isa(fmat[i-1,2], Function) "missing function !"
				fmat[i,2] = (v) -> tuple(fmat[i-1,2](v)..., f1(v))
			end
		end
	end

	MCMCLikelihoodModel(lik,
						fmat[1,1], fmat[1,2],
						fmat[2,1], fmat[2,2],
						fmat[3,1], fmat[3,2],
						init, scale, pmap)
end

# Model creation with multivariate Distribution as input

function MCMCLikelihoodModel(d::MultivariateDistribution;
  grad::FunctionOrNothing = nothing,
  tensor::FunctionOrNothing = nothing,
  dtensor::FunctionOrNothing = nothing,
  allgrad::FunctionOrNothing = nothing,
  alltensor::FunctionOrNothing = nothing,
  alldtensor::FunctionOrNothing = nothing,
	init::F64OrVectorF64 = [1.0], scale::F64OrVectorF64 = 1.0, pmap::Union(Nothing, Dict) = nothing)
  @assert method_exists(logpdf, (typeof(d), Vector{Float64})) "logpdf function not defined for $d"

  fout = Array(FunctionOrNothing, 2)
  fout[1] = x::Vector{Float64} -> logpdf(d, x)

  if grad == nothing
  	if method_exists(gradlogpdf, (typeof(d), Vector{Float64}))
  		fout[2] = x::Vector{Float64} -> gradlogpdf(d, x)
  	end
  end

	MCMCLikelihoodModel(fout[1]; grad=fout[2], tensor=tensor, dtensor=dtensor,
    allgrad=allgrad, alltensor=alltensor, alldtensor=alldtensor,
    init=init, scale=scale, pmap=pmap)
end

# Model creation with univariate Distribution as input

function MCMCLikelihoodModel(d::UnivariateDistribution;
  grad::FunctionOrNothing = nothing,
  tensor::FunctionOrNothing = nothing,
  dtensor::FunctionOrNothing = nothing,
  allgrad::FunctionOrNothing = nothing,
  alltensor::FunctionOrNothing = nothing,
  alldtensor::FunctionOrNothing = nothing,
	init::F64OrVectorF64 = [1.0], scale::F64OrVectorF64 = 1.0, pmap::Union(Nothing, Dict) = nothing)
  @assert method_exists(logpdf, (typeof(d), Float64)) "logpdf function not defined for $d"

  fout = Array(FunctionOrNothing, 2)
  fout[1] = x::Vector{Float64} -> logpdf(d, x[1])

  if grad == nothing
  	if method_exists(gradlogpdf, (typeof(d), Float64))
  		fout[2] = x::Vector{Float64} -> [gradlogpdf(d, x[1])]
  	end
  end

	MCMCLikelihoodModel(fout[1]; grad=fout[2], tensor=tensor, dtensor=dtensor,
    allgrad=allgrad, alltensor=alltensor, alldtensor=alldtensor,
    init=init, scale=scale, pmap=pmap)
end


text:
"""
    FVector(type::VectorSpaceType, data, basis::AbstractBasis)

Decorator indicating that the vector `data` contains coordinates of a vector from a fiber
of a vector bundle of type `type`. `basis` is an object describing the basis of that space
in which the coordinates are given.

Conversion between `FVector` representation and the default representation of an object
(for example a tangent vector) for a manifold should be done using [`get_coordinates`](@ref)
and [`get_vector`](@ref).

# Examples

```julia-repl
julia> using Manifolds

julia> M = Sphere(2)
Sphere(2, ℝ)

julia> p = [1.0, 0.0, 0.0]
3-element Vector{Float64}:
 1.0
 0.0
 0.0

julia> X = [0.0, 2.0, -1.0]
3-element Vector{Float64}:
  0.0
  2.0
 -1.0

julia> B = DefaultOrthonormalBasis()
DefaultOrthonormalBasis(ℝ)

julia> fX = TFVector(get_coordinates(M, p, X, B), B)
TFVector([2.0, -1.0], DefaultOrthonormalBasis(ℝ))

julia> X_back = get_vector(M, p, fX.data, fX.basis)
3-element Vector{Float64}:
 -0.0
  2.0
 -1.0
```
"""
struct FVector{TType<:VectorSpaceType,TData,TBasis<:AbstractBasis}
    type::TType
    data::TData
    basis::TBasis
end

const TFVector = FVector{TangentSpaceType}
const CoTFVector = FVector{CotangentSpaceType}

function TFVector(data, basis::AbstractBasis)
    return TFVector{typeof(data),typeof(basis)}(TangentSpace, data, basis)
end
function CoTFVector(data, basis::AbstractBasis)
    return CoTFVector{typeof(data),typeof(basis)}(CotangentSpace, data, basis)
end

function Base.show(io::IO, fX::TFVector)
    return print(io, "TFVector(", fX.data, ", ", fX.basis, ")")
end
function Base.show(io::IO, fX::CoTFVector)
    return print(io, "CoTFVector(", fX.data, ", ", fX.basis, ")")
end

"""
    AbstractFibreVector{TType<:VectorSpaceType}

Type for a vector from a vector space (fibre of a vector bundle) of type `TType` of a manifold.
While a [`AbstractManifold`](@ref) does not necessarily require this type, for example when it is
implemented for `Vector`s or `Matrix` type elements, this type can be used for more
complicated representations, semantic verification, or even dispatch for different
representations of tangent vectors and their types on a manifold.
"""
abstract type AbstractFibreVector{TType<:VectorSpaceType} end

"""
    TVector = AbstractFibreVector{TangentSpaceType}

Type for a tangent vector of a manifold. While a [`AbstractManifold`](@ref) does not necessarily
require this type, for example when it is implemented for `Vector`s or `Matrix` type
elements, this type can be used for more complicated representations, semantic verification,
or even dispatch for different representations of tangent vectors and their types on a
manifold.
"""
const TVector = AbstractFibreVector{TangentSpaceType}

"""
    CoTVector = AbstractFibreVector{CotangentSpaceType}

Type for a cotangent vector of a manifold. While a [`AbstractManifold`](@ref) does not necessarily
require this type, for example when it is implemented for `Vector`s or `Matrix` type
elements, this type can be used for more complicated representations, semantic verification,
or even dispatch for different representations of cotangent vectors and their types on a
manifold.
"""
const CoTVector = AbstractFibreVector{CotangentSpaceType}

Base.:+(X::FVector, Y::FVector) = FVector(X.type, X.data + Y.data, X.basis)

Base.:-(X::FVector, Y::FVector) = FVector(X.type, X.data - Y.data, X.basis)
Base.:-(X::FVector) = FVector(X.type, -X.data, X.basis)

Base.:*(a::Number, X::FVector) = FVector(X.type, a * X.data, X.basis)

allocate(x::FVector) = FVector(x.type, allocate(x.data), x.basis)
allocate(x::FVector, ::Type{T}) where {T} = FVector(x.type, allocate(x.data, T), x.basis)

function Base.copyto!(X::FVector, Y::FVector)
    copyto!(X.data, Y.data)
    return X
end

function number_eltype(
    ::Type{FVector{TType,TData,TBasis}},
) where {TType<:VectorSpaceType,TData,TBasis}
    return number_eltype(TData)
end
number_eltype(v::FVector) = number_eltype(v.data)

Base.show(io::IO, ::TangentSpaceType) = print(io, "TangentSpace")
Base.show(io::IO, ::CotangentSpaceType) = print(io, "CotangentSpace")

"""
    vector_space_dimension(M::AbstractManifold, V::VectorSpaceType)

Dimension of the vector space of type `V` on manifold `M`.
"""
vector_space_dimension(::AbstractManifold, ::VectorSpaceType)

function vector_space_dimension(M::AbstractManifold, ::TCoTSpaceType)
    return manifold_dimension(M)
end


text:
int_rules_3_2_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 3.2.2*(f+g*x)^m*(h+i*x)^q*(A+B*log(e*((a+b*x)*over*(c+d*x))^n))^p =#
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(h') + ~(i') * ~x) * (~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))), ~x) => ((~f + ~g * ~x) ^ (~m + 1) * (~h + ~i * ~x) * (~A + ~B * log(~e * ((~a + ~b * ~x) / (~c + ~d * ~x)) ^ ~n))) / (~g * (~m + 2)) + ((~i * (~b * ~c - ~a * ~d)) / (~b * ~d * (~m + 2))) * Antiderivative((~f + ~g * ~x) ^ ~m * ((~A - ~B * ~n) + ~B * log(~e * ((~a + ~b * ~x) / (~c + ~d * ~x)) ^ ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~A, ~B, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ(~b * ~f - ~a * ~g, 0) && (EqQ(~d * ~h - ~c * ~i, 0) && IGtQ(~m, -2))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(h') + ~(i') * ~x) * (~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)), ~x) => ((~f + ~g * ~x) ^ (~m + 1) * (~h + ~i * ~x) * (~A + ~B * log((~e * (~a + ~b * ~x) ^ ~n) / (~c + ~d * ~x) ^ ~n))) / (~g * (~m + 2)) + ((~i * (~b * ~c - ~a * ~d)) / (~b * ~d * (~m + 2))) * Antiderivative((~f + ~g * ~x) ^ ~m * ((~A - ~B * ~n) + ~B * log((~e * (~a + ~b * ~x) ^ ~n) / (~c + ~d * ~x) ^ ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~A, ~B, ~m, ~n], ~x) && (EqQ(~n + ~mn, 0) && (IGtQ(~n, 0) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ(~b * ~f - ~a * ~g, 0) && (EqQ(~d * ~h - ~c * ~i, 0) && IGtQ(~m, -2))))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(h') + ~(i') * ~x) ^ ~(q') * (~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) ^ ~(p'), ~x) => (~b * ~c - ~a * ~d) ^ (~m + ~q + 1) * (~g / ~b) ^ ~m * (~i / ~d) ^ ~q * Subst(Antiderivative(((~x) ^ ~m * (~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p) / (~b - ~d * ~x) ^ (~m + ~q + 2), ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~A, ~B, ~n, ~p], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ(~b * ~f - ~a * ~g, 0) && (EqQ(~d * ~h - ~c * ~i, 0) && IntegersQ(~m, ~q))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(h') + ~(i') * ~x) ^ ~(q') * (~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) ^ ~(p'), ~x) => (~b * ~c - ~a * ~d) ^ (~m + ~q + 1) * (~g / ~b) ^ ~m * (~i / ~d) ^ ~q * Subst(Antiderivative(((~x) ^ ~m * (~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p) / (~b - ~d * ~x) ^ (~m + ~q + 2), ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~A, ~B, ~n, ~p], ~x) && (EqQ(~n + ~mn, 0) && (IGtQ(~n, 0) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ(~b * ~f - ~a * ~g, 0) && (EqQ(~d * ~h - ~c * ~i, 0) && IntegersQ(~m, ~q))))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(h') + ~(i') * ~x) ^ ~(q') * (~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) ^ ~(p'), ~x) => (((~d) ^ 2 * ((~g * (~a + ~b * ~x)) / ~b) ^ ~m) / ((~i) ^ 2 * (~b * ~c - ~a * ~d) * ((~i * (~c + ~d * ~x)) / ~d) ^ ~m * ((~a + ~b * ~x) / (~c + ~d * ~x)) ^ ~m)) * Subst(Antiderivative((~x) ^ ~m * (~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p, ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~A, ~B, ~m, ~n, ~p, ~q], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ(~b * ~f - ~a * ~g, 0) && (EqQ(~d * ~h - ~c * ~i, 0) && EqQ(~m + ~q + 2, 0))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(h') + ~(i') * ~x) ^ ~(q') * (~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) ^ ~(p'), ~x) => (((~d) ^ 2 * ((~g * (~a + ~b * ~x)) / ~b) ^ ~m) / ((~i) ^ 2 * (~b * ~c - ~a * ~d) * ((~i * (~c + ~d * ~x)) / ~d) ^ ~m * ((~a + ~b * ~x) / (~c + ~d * ~x)) ^ ~m)) * Subst(Antiderivative((~x) ^ ~m * (~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p, ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~A, ~B, ~m, ~n, ~p, ~q], ~x) && (EqQ(~n + ~mn, 0) && (IGtQ(~n, 0) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ(~b * ~f - ~a * ~g, 0) && (EqQ(~d * ~h - ~c * ~i, 0) && EqQ(~m + ~q + 2, 0))))))
    #= Antiderivative(((~f')+(~g')*(~x))^(~m')*((~h')+(~i')*(~x))^(~q')*((~A')+(~B')*log((~e')*((~a')+(~b')*(~x)) ^(~n')*((~c')+(~d')*(~x))^(~mn)))^(~p'),~x) :=  b*d*(f+g*x)^(m+1)/(g*i*(b*c-a*d)*(h+i*x)^(m+1)*((a+b*x)/(c+d*x))^(m+1) )* Subst(Antiderivative(x^m*(A+B*log(e*x^n))^p,x),x,(a+b*x)/(c+d*x)) <--  FreeQ([a,b,c,d,e,f,g,h,i,A,B,m,n,p,q],x) && EqQ(n+mn,0) &&  NeQ(b*c-a*d,0) && EqQ(b*f-a*g,0) && EqQ(d*h-c*i,0) && EqQ(m+q+2,0) =#
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(h') + ~(i') * ~x) ^ ~(q') * (~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) ^ ~(p'), ~x) => (~b * ~c - ~a * ~d) ^ (~q + 1) * (~i / ~d) ^ ~q * Subst(Antiderivative((((~b * ~f - ~a * ~g) - (~d * ~f - ~c * ~g) * ~x) ^ ~m * (~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p) / (~b - ~d * ~x) ^ (~m + ~q + 2), ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegersQ(~m, ~q) && (IGtQ(~p, 0) && EqQ(~d * ~h - ~c * ~i, 0))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(h') + ~(i') * ~x) ^ ~(q') * (~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) ^ ~(p'), ~x) => (~b * ~c - ~a * ~d) ^ (~q + 1) * (~i / ~d) ^ ~q * Subst(Antiderivative((((~b * ~f - ~a * ~g) - (~d * ~f - ~c * ~g) * ~x) ^ ~m * (~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p) / (~b - ~d * ~x) ^ (~m + ~q + 2), ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~A, ~B, ~n], ~x) && (EqQ(~n + ~mn, 0) && (IGtQ(~n, 0) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegersQ(~m, ~q) && (IGtQ(~p, 0) && EqQ(~d * ~h - ~c * ~i, 0))))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(h') + ~(i') * ~x) ^ ~(q') * (~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) ^ ~(p'), ~x) => (~b * ~c - ~a * ~d) * Subst(Antiderivative((((~b * ~f - ~a * ~g) - (~d * ~f - ~c * ~g) * ~x) ^ ~m * ((~b * ~h - ~a * ~i) - (~d * ~h - ~c * ~i) * ~x) ^ ~q * (~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p) / (~b - ~d * ~x) ^ (~m + ~q + 2), ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegersQ(~m, ~q) && IGtQ(~p, 0)))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(h') + ~(i') * ~x) ^ ~(q') * (~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) ^ ~(p'), ~x) => (~b * ~c - ~a * ~d) * Subst(Antiderivative((((~b * ~f - ~a * ~g) - (~d * ~f - ~c * ~g) * ~x) ^ ~m * ((~b * ~h - ~a * ~i) - (~d * ~h - ~c * ~i) * ~x) ^ ~q * (~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p) / (~b - ~d * ~x) ^ (~m + ~q + 2), ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~A, ~B, ~n], ~x) && (EqQ(~n + ~mn, 0) && (IGtQ(~n, 0) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegersQ(~m, ~q) && IGtQ(~p, 0)))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(h') + ~(i') * ~x) ^ ~(q') * (~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) ^ ~(p'), ~x) => Unintegrable((~f + ~g * ~x) ^ ~m * (~h + ~i * ~x) ^ ~q * (~A + ~B * log(~e * ((~a + ~b * ~x) / (~c + ~d * ~x)) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~A, ~B, ~m, ~n, ~p, ~q], ~x)
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(h') + ~(i') * ~x) ^ ~(q') * (~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) ^ ~(p'), ~x) => Unintegrable((~f + ~g * ~x) ^ ~m * (~h + ~i * ~x) ^ ~q * (~A + ~B * log((~e * (~a + ~b * ~x) ^ ~n) / (~c + ~d * ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~A, ~B, ~m, ~n, ~p, ~q], ~x) && (EqQ(~n + ~mn, 0) && IntegerQ(~n))
    @apply_utils Antiderivative((~w) ^ ~(m') * (~y) ^ ~(q') * (~(A') + ~(B') * log(~(e') * (~u / ~v) ^ ~(n'))) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~w, ~x) ^ ~m * ExpandToSum(~y, ~x) ^ ~q * (~A + ~B * log(~e * (ExpandToSum(~u, ~x) / ExpandToSum(~v, ~x)) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~e, ~A, ~B, ~m, ~n, ~p, ~q], ~x) && (LinearQ([~u, ~v, ~w, ~y], ~x) && Not(LinearMatchQ([~u, ~v, ~w, ~y], ~x)))
    @apply_utils Antiderivative((~w) ^ ~(m') * (~y) ^ ~(q') * (~(A') + ~(B') * log(~(e') * (~u) ^ ~(n') * (~v) ^ ~mn)) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~w, ~x) ^ ~m * ExpandToSum(~y, ~x) ^ ~q * (~A + ~B * log((~e * ExpandToSum(~u, ~x) ^ ~n) / ExpandToSum(~v, ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~e, ~A, ~B, ~m, ~n, ~p, ~q], ~x) && (EqQ(~n + ~mn, 0) && (IGtQ(~n, 0) && (LinearQ([~u, ~v, ~w, ~y], ~x) && Not(LinearMatchQ([~u, ~v, ~w, ~y], ~x)))))
    @apply_utils Antiderivative(~(w') * (~(A') + ~(B') * log(~(e') * (~u) ^ ~(n') * (~v) ^ ~mn)) ^ ~(p'), ~x) => Subst(Antiderivative(~w * (~A + ~B * log(~e * (~u / ~v) ^ ~n)) ^ ~p, ~x), ~e * (~u / ~v) ^ ~n, (~e * (~u) ^ ~n) / (~v) ^ ~n) <-- FreeQ([~e, ~A, ~B, ~n, ~p], ~x) && (EqQ(~n + ~mn, 0) && (LinearQ([~u, ~v], ~x) && Not(IntegerQ(~n))))
    #= Antiderivative((~w')*((~A')+(~B')*log((~e')*((~f')*(~u)^(~q')*(~v)^(~mq))^(~n')))^(~p'),~x) :=   Subst(Antiderivative(w*(A+B*log(e*f^n*(u/v)^(n*q)))^p,x),e*f^n*(u/v)^(n*q),e*( f*(u^q/v^q))^n) <-- FreeQ([e,f,A,B,n,p,q],x) && EqQ(q+mq,0) &&  LinearQ([u,v],x) && Not(IntegerQ(n)) =#
    @apply_utils Antiderivative((~(f') + ~(g') * ~x + ~(h') * (~x) ^ 2) ^ ~(m') * (~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) ^ ~(p'), ~x) => ((~h) ^ ~m / ((~b) ^ ~m * (~d) ^ ~m)) * Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~m * (~A + ~B * log(~e * ((~a + ~b * ~x) / (~c + ~d * ~x)) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~A, ~B, ~n, ~p], ~x) && (EqQ(~b * ~d * ~f - ~a * ~c * ~h, 0) && (EqQ(~b * ~d * ~g - ~h * (~b * ~c + ~a * ~d), 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x + ~(h') * (~x) ^ 2) ^ ~(m') * (~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) ^ ~(p'), ~x) => ((~h) ^ ~m / ((~b) ^ ~m * (~d) ^ ~m)) * Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~m * (~A + ~B * log((~e * (~a + ~b * ~x) ^ ~n) / (~c + ~d * ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~A, ~B, ~n, ~p], ~x) && (EqQ(~n + ~mn, 0) && (IGtQ(~n, 0) && (EqQ(~b * ~d * ~f - ~a * ~c * ~h, 0) && (EqQ(~b * ~d * ~g - ~h * (~b * ~c + ~a * ~d), 0) && IntegerQ(~m)))))
    @apply_utils Antiderivative((~P2x) ^ ~(m') * (~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) ^ ~(p'), ~x) => With([f = Coeff(~P2x, ~x, 0), g = Coeff(~P2x, ~x, 1), h = Coeff(~P2x, ~x, 2)], (~b * ~c - ~a * ~d) * Subst(Antiderivative(((((((~b) ^ 2 * f - ~a * ~b * g) + (~a) ^ 2 * h) - (((2 * ~b * ~d * f - ~b * ~c * g) - ~a * ~d * g) + 2 * ~a * ~c * h) * ~x) + (((~d) ^ 2 * f - ~c * ~d * g) + (~c) ^ 2 * h) * (~x) ^ 2) ^ ~m * (~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p) / (~b - ~d * ~x) ^ (2 * (~m + 1)), ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B, ~n], ~x) && (PolyQ(~P2x, ~x, 2) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegerQ(~m) && IGtQ(~p, 0))))
    @apply_utils Antiderivative((~P2x) ^ ~(m') * (~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) ^ ~(p'), ~x) => With([f = Coeff(~P2x, ~x, 0), g = Coeff(~P2x, ~x, 1), h = Coeff(~P2x, ~x, 2)], (~b * ~c - ~a * ~d) * Subst(Antiderivative(((((((~b) ^ 2 * f - ~a * ~b * g) + (~a) ^ 2 * h) - (((2 * ~b * ~d * f - ~b * ~c * g) - ~a * ~d * g) + 2 * ~a * ~c * h) * ~x) + (((~d) ^ 2 * f - ~c * ~d * g) + (~c) ^ 2 * h) * (~x) ^ 2) ^ ~m * (~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p) / (~b - ~d * ~x) ^ (2 * (~m + 1)), ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B, ~n], ~x) && (PolyQ(~P2x, ~x, 2) && (EqQ(~n + ~mn, 0) && (IGtQ(~n, 0) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegerQ(~m) && IGtQ(~p, 0))))))
end

max_stars_repo_path: src/rules/3 Logarithms/3.2.2 (f+g x)^m (h+i x)^q (A+B log(e ((a+b x) over (c+d x))^n))^p.jl


text:
export PSO, start!, step!, release!

immutable PSO{C, A} <: PopulationBasedOptAlgo{C}
    ctx::C
    populationSize::Int
    dnaSize::Int
    vMax::Float32
    w::Float32
    c1::Float32
    c2::Float32
    mutationChance::Float32
    mutationStrength::Float32
    comparer::AbstractComparer
    popMan::PopulationManager
    pBests::Vector{Entity}
    velocities::Vector{A}
end

PSO(
    ctx,
    dnaSize,
    comparer::AbstractComparer,
    decode::Function;
    populationSize = 100,
    vMax = 0.1f0,
    w = 0.8f0,
    c1 = 2f0,
    c2 = 2f0,
    mutationChance = 0.02f0,
    mutationStrength = 0.05f0) =
PSO(
    ctx,
    populationSize,
    dnaSize,
    vMax,
    w,
    c1,
    c2,
    mutationChance,
    mutationStrength,
    comparer,
    PopulationManager(
        ctx,
        populationSize,
        dnaSize,
        comparer,
        dna -> begin
            dna2 = normalized(dna)
            local result
            try
                result = decode(dna2)
            finally
                if result != dna2
                    release!(dna2)
                end
            end
            result
        end),
    Vector{Entity}(),
    [(x -> array(ctx))(x) for x in 1:populationSize])

function release!(pso::PSO)
    release!(pso.popMan)
    for best in pso.pBests
        release!(best)
    end
    empty!(pso.pBests)
    for vel in pso.velocities
        release!(vel)
    end
    empty!(pso.velocities)
end

function start!(pso::PSO)
    reset!(pso.comparer)
    randomize!(pso.popMan)
    for idx in 1:pso.populationSize
        push!(pso.pBests, copy(pso.popMan.population[idx]))
        scope!(pso.ctx) do this
            pso.velocities[idx][] = randu(pso.ctx, Float32, pso.dnaSize) .* (pso.vMax * 2) .- pso.vMax
        end
    end
end

function step!(pso::PSO)
    reset!(pso.comparer)

    candidate = createCandidate(pso.popMan)

    population = pso.popMan.population
    gBest = get(pso.popMan.best).dna.array
    for idx in 1:pso.populationSize
        pBest = pso.pBests[idx]
        dna = scope!(pso.ctx) do this
            rnd1 = randu(pso.ctx, Float32, pso.dnaSize)
            rnd2 = randu(pso.ctx, Float32, pso.dnaSize)
            loc = population[idx].dna.array
            vel = pso.velocities[idx]
            newVel = (pso.w .* vel) .+ (rnd1 .* pso.c1) .* (pBest.dna.array .- loc) .+ (rnd2 .* pso.c2) .* (gBest .- loc)
            newVel = maxOf(newVel, -pso.vMax)
            newVel = minOf(newVel, pso.vMax)
            vel[] = newVel
            eval!(vel)
            loc = loc .+ vel
            childDNA = DNA(pso.ctx, loc)
            this.result(childDNA.array)
            mutate!(childDNA, pso.mutationChance, pso.mutationStrength, false)
            eval!(childDNA)
        end
        entity = push!(candidate, dna)
        if pso.comparer(entity.body, pBest.body)
            release!(pBest)
            pso.pBests[idx] = copy(entity)
        end
    end

    set!(pso.popMan, candidate, true)
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
struct DetectionCallback
    detection_times::Vector{OptTimePoint}
    detection_types::Vector{UInt8}

    tracking_times::Vector{OptTimePoint}
    tracking_sources::Vector{UInt32}
    tracking_types::Vector{UInt8}

    max_num_infected::UInt32
end
DetectionCallback(sz::Integer, max_num_infected::Integer=10^8) = DetectionCallback(
    Vector{OptTimePoint}(missing, sz),
    fill(UInt8(0), sz),
    Vector{OptTimePoint}(missing, sz),
    fill(UInt32(0), sz),
    fill(UInt8(0), sz),
    max_num_infected
)

function saveparams(dict, cb::DetectionCallback, prefix::AbstractString="")
  dict[prefix*"detection_times"] = optreal2float32.(cb.detection_times)
  dict[prefix*"detection_types"] = cb.detection_types

  dict[prefix*"tracking_times"] = optreal2float32.(cb.tracking_times)
  dict[prefix*"tracking_sources"] = cb.tracking_sources
  dict[prefix*"tracking_types"] = cb.tracking_types
end

function reset!(cb::DetectionCallback)
  fill!(cb.detection_times, missing)
  fill!(cb.detection_types, 0)
  fill!(cb.tracking_sources, 0)
  fill!(cb.tracking_types, 0)
end

function (cb::DetectionCallback)(event::MocosSim.Event, state::MocosSim.SimState, params::MocosSim.SimParams)
  eventkind = MocosSim.kind(event)
  contactkind = MocosSim.contactkind(event)
  subject = MocosSim.subject(event)
  if MocosSim.isdetection(eventkind)
    cb.detection_times[subject] = MocosSim.time(event)
    cb.detection_types[subject] = MocosSim.detectionkind(event) |> UInt8
  elseif MocosSim.istracking(eventkind)
    cb.tracking_times[subject] = MocosSim.time(event)
    cb.tracking_sources[subject] = MocosSim.source(event)
    cb.tracking_types[subject] = MocosSim.trackingkind(event) |> UInt8
  end
  return MocosSim.numinfected(state.stats) < cb.max_num_infected
end

function save_infections_and_detections(path::AbstractString, simstate::MocosSim.SimState, callback::DetectionCallback)
  f = jldopen(path, "w", compress=true)
  try
    MocosSim.saveparams(f, simstate)
    saveparams(f, callback)
  finally
    close(f)
  end
  nothing
end

max_stars_repo_name: Novirmor/MocosSimLauncher

max_issues_repo_name: Novirmor/MocosSimLauncher

max_forks_repo_name: Novirmor/MocosSimLauncher


text:
module Addition

# include and implement the Addition interface
include("../interface/addition/IAddition.jl")
function doAdd(a, b)
    a + b
end
end

max_stars_repo_name: sonolac/JuliaInterface

max_issues_repo_name: sonolac/JuliaInterface

max_forks_repo_name: sonolac/JuliaInterface


text:
## ICD10 #####
"""
    ICD10

An ICD-10 code (standard WHO version).

Fields are `level`, indicating if it is a 3, 4, or 5 digit code, and `data` containing the code.
"""
struct ICD10 <: AbstractICD10
  level::Int8
  data::NTuple{5,UInt8}
end

"""
    ICD10(str::String, validateinput = false)

Create an ICD10 code from a string.

ICD10 codes have the form `ANN[.][N[N]]` where `A` is a letter A-Z and
`N` is a digit 0-9. Parts in brackets are optional.

If `validateinput=true` then the input string is checked for valid format using
a regex.
"""
function ICD10(str::AbstractString, validateinput = false)
  punct = occursin(".", str) ? true : false
  if validateinput
    validICD10input(str) || throw(
      DomainError(
        str,
        "ICD-10 codes should have format `ANN[.][N[N]]` where `A` is letter A-Z, `N` is decimal digit, and parts in brackets are optional",
      ),
    )
  end
  level = punct ? length(str) - 1 : length(str)
  ch4 = punct ? 5 : 4

  if level == 3
    return ICD10(level, (NTuple{3,UInt8}.(str)..., UInt8(0), UInt8(0)))
  elseif level == 4
    return ICD10(level, (NTuple{3,UInt8}.(str)..., UInt8(str[ch4]), UInt8(0)))
  elseif level == 5
    return ICD10(level, (NTuple{3,UInt8}.(str)..., UInt8(str[ch4]), UInt8(str[ch4+1])))
  else
    throw(DomainError(str, "ICD-10 codes should be 3-5 characters excluding the period"))
  end
end

"""
    ICD10(icd::T) where T <: AbstractICD10

Create an ICD-10 code from another type of ICD-10 code.

Note this just changes the type. It does not translate concepts between versions.
"""
ICD10(icd::T) where {T<:AbstractICD10} = ICD10(icd)

function validICD10input(str)
  icdfmt = r"^[[:upper:]][[:digit:]]{2}(\.){0,1}([[:digit:]]{1,2}){0,1}$"
  return occursin(icdfmt, str)
end


text:
module MesoscaleML

"""
    func(x)

Returns double the number `x` plus `1`.
"""
func(x) = 2x + 1
export func

end # module


text:
# ---
# title: 1593. Split a String Into the Max Number of Unique Substrings
# id: problem1593
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Backtracking
# link: <https://leetcode.com/problems/split-a-string-into-the-max-number-of-unique-substrings/description/>
# hidden: true
# ---
#
# Given a string `s`, return _the maximum  number of unique substrings that the
# given string can be split into_.
#
# You can split string `s` into any list of  **non-empty substrings** , where
# the concatenation of the substrings forms the original string. However, you
# must split the substrings such that all of them are **unique**.
#
# A **substring** is a contiguous sequence of characters within a string.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "ababccc"
#     Output: 5
#     **Explanation** : One way to split maximally is ['a', 'b', 'ab', 'c', 'cc']. Splitting like ['a', 'b', 'a', 'b', 'c', 'cc'] is not valid as you have 'a' and 'b' multiple times.
#
#
# **Example 2:**
#
#
#
#     Input: s = "aba"
#     Output: 2
#     **Explanation** : One way to split maximally is ['a', 'ba'].
#
#
# **Example 3:**
#
#
#
#     Input: s = "aa"
#     Output: 1
#     **Explanation** : It is impossible to split the string any further.
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 16`
#
#   * `s` contains only lower case English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "Testing Coloring Submodule" begin

@test ColoringParams(0, 1, 0) == ColoringParams(0, 1, "gray")

Random.seed!(1)

N = 64
I = Images.shepp_logan(N,N)

IC = colorize(I, ColoringParams(0,1,"gray"), 0.0, 1.0)
@test typeof(IC) == Array{RGBA{Normed{UInt8,8}},2}

exportImage("img/coloring1.png", IC)
@testImg("coloring1.png")

IC = colorize(I, ColoringParams(0,1,"viridis"), 0.0, 1.0)
exportImage("img/coloring2.png", IC)
@testImg("coloring2.png")

exportImage("img/coloring3.png", I, vmin=0.0, vmax=0.5, colormap="blue", normalize=true)
@testImg("coloring3.png")

# Movie

Itemp = Array(reshape(range(0,1,length=90000),30,30,100))
exportMovie("img/movie.gif", Itemp, vmin=0.0, vmax=1.0, colormap="viridis", normalize=true)
@testImg("movie.gif")

exportMovies("img/movie.gif", [Itemp,Itemp,Itemp])
@testImg("movie_xy.gif")
@testImg("movie_xz.gif")
@testImg("movie_yz.gif")

# Overlay

O = zeros(N,N)
O[10:30,10:30] .= 1.0

OC = colorize(O, ColoringParams(0,1,"red"), 0.0, 1.0)

over = overlay(IC, OC, cmap("red")[1])
exportImage("img/coloring4.png", over)
@testImg("coloring4.png")

over = blend(IC, OC)
exportImage("img/coloring5.png", over)
@testImg("coloring5.png")

I = reshape(Images.shepp_logan(N,N),1,N,N)
O = reshape(O,1,N,N)
I = cat(I,O,dims=1)

OC = colorize(I, [ColoringParams(0,1,"gray"), ColoringParams(0,1,"red"),],
     [0.0,0.0], [1.0,1.0])

exportImage("img/coloring6.png", OC)
#@testImg("coloring6.png")

end


text:
module BasicStructs
using Parameters
export ImageConstants, ConfigurtationStruct, ResultMetrics
"""
constants associated with image over which we will evaluate segmentations
"""
@with_kw struct ImageConstants
    mvspx ::Float64 # Voxelspacing x
    mvspy::Float64 # Voxelspacing y
    mvspz::Float64 #mean Voxelspacing z
    isZConst::Bool # true if slices thickness is the same in all image
    ZPositions::Vector{Float64} # array of true physical positions (in mm) of slices relative to the begining - used in case we have variable thickness of slices
    numberOfVox::Int64 # number of voxels in image
end #ImageConstants
"""
configuration struct that when passed will marks what kind of metrics we are intrested in

    """
@with_kw struct ConfigurtationStruct
    anyFuzzy::Bool = false# is any of the metric calculated fuzzy

end #ConfigurtationStruct

"""
Struct holding all resulting metrics - if some metric was not calculated its value is just -1
"""
@with_kw struct ResultMetrics
    dice::Float64 = -1.0 #dice coefficient
    jaccard::Float64 =  -1.0 #jaccard coefficient
    gce::Float64 =  -1.0 #global consistency error
    vol::Float64 =  -1.0 # Volume metric
    randInd::Float64 = -1.0 # Rand Index
    mi::Float64 = -1.0 # mutual information
    ic::Float64 = -1.0 # interclass correlation
    Kc::Float64 = -1.0 # Kohen Cappa
    Mahalanobis::Float64 = -1.0 # Mahalanobis distance
    Hausdorff::Float64 = -1.0 # Hausdorff distance

end #ResultMetrics

end#BasicStructs


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/NodeNumbering.jl/blob/master/LICENSE

"""
    adjacency_visualization(RCM_adjacency::Dict{Int, Vector{Int}})

Visualizes adjacency graph Dicts as matrices. Result can be plotted using `PyPlot` function `matshow()`.
A part of the Reverse Cuthill-McKee Algorithm.

# Examples

```jldoctest
julia> RCM_adjacency = Dict{Int, Vector{Int}}(
                       7 => [5, 6, 9],
                       9 => [7, 8],
                       4 => [2, 6],
                       2 => [4, 3, 1],
                       3 => [5, 6, 2, 1],
                       8 => [6, 9],
                       5 => [7, 3],
                       6 => [7, 3, 8, 4],
                       1 => [3, 2]);

julia> adjacency_visualization(RCM_adjacency)
9×9 Array{Float64,2}:
 1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0
 1.0  1.0  1.0  0.0  1.0  1.0  0.0  0.0  0.0
 0.0  1.0  0.0  1.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  1.0  0.0  1.0  0.0  0.0
 0.0  0.0  1.0  1.0  0.0  1.0  1.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0  1.0  1.0  0.0  1.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  1.0  1.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0

```
"""
function adjacency_visualization(RCM_adjacency::Dict{Int, Vector{Int}})
    m=eye(length(RCM_adjacency))
    for (k, v) in RCM_adjacency
        for i in v
          m[k, i]=1
        end
    end
    return m
end


text:
# test ModelSetup/Medium/initwf.jl
include("typedefine.jl")
include("initwf.jl")

nDZ = 100
nHX = 100
ext = 10
iflag = 2

elwf = initelwf(nDZ,nHX,ext,iflag)

acwf = initacwf(nDZ,nHX,ext,iflag)

nHY = 100

elwf3d = initelwf(nDZ,nHX,nHY,ext,iflag)

acwf3d = initacwf(nDZ,nHX,nHY,ext,iflag)


text:
@testset "ESS" begin
    @model function demo(x)
        m ~ Normal()
        x ~ Normal(m, 0.5)
    end
    demo_default = demo(1.0)

    @model function demodot(x)
        m = Vector{Float64}(undef, 2)
        @. m ~ Normal()
        x ~ Normal(m[2], 0.5)
    end
    demodot_default = demodot(1.0)

    @turing_testset "ESS constructor" begin
        Random.seed!(0)
        N = 500

        s1 = ESS()
        s2 = ESS(:m)
        for s in (s1, s2)
            @test DynamicPPL.alg_str(Sampler(s, demo_default)) == "ESS"
        end

        c1 = sample(demo_default, s1, N)
        c2 = sample(demo_default, s2, N)
        c3 = sample(demodot_default, s1, N)
        c4 = sample(demodot_default, s2, N)

        s3 = Gibbs(ESS(:m), MH(:s))
        c5 = sample(gdemo_default, s3, N)
    end

    @numerical_testset "ESS inference" begin
        Random.seed!(1)
        chain = sample(demo_default, ESS(), 5_000)
        check_numerical(chain, [:m], [0.8], atol = 0.1)

        Random.seed!(1)
        chain = sample(demodot_default, ESS(), 5_000)
        check_numerical(chain, ["m[1]", "m[2]"], [0.0, 0.8], atol = 0.1)

        Random.seed!(100)
        alg = Gibbs(
            CSMC(15, :s),
            ESS(:m))
        chain = sample(gdemo(1.5, 2.0), alg, 10_000)
        check_numerical(chain, [:s, :m], [49/24, 7/6], atol=0.1)

        # MoGtest
        Random.seed!(125)
        alg = Gibbs(
            CSMC(15, :z1, :z2, :z3, :z4),
            ESS(:mu1), ESS(:mu2))
        chain = sample(MoGtest_default, alg, 6000)
        check_MoGtest_default(chain, atol = 0.1)

        # Different "equivalent" models.
        Random.seed!(125)
        check_gdemo_models(ESS(), 1_000)
    end
end


text:
# Author: Gabrio Rizzuti, grizzuti3@gatech.edu
# Date: October 2020

using InvertibleNetworks, LinearAlgebra, Test, Statistics

###############################################################################
# Initialization and invertibility

# Input
nx = 28
ny = 28
nc = 4
batchsize = 10
X = rand(Float32, nx, ny, nc, batchsize)

# Layer and initialization
AL = AffineLayer(nx, ny, nc; logdet=false)
Y = AL.forward(X)

# Test invertibility
@test isapprox(norm(X - AL.inverse(AL.forward(X)))/norm(X), 0f0, atol=1f-6)
@test isapprox(norm(X - AL.forward(AL.inverse(X)))/norm(X), 0f0, atol=1f-6)

# Test with logdet enabled
AL = AffineLayer(nx, ny, nc; logdet=true)
Y, lgdt = AL.forward(X)

# Test invertibility
@test isapprox(norm(X - AL.inverse(AL.forward(X)[1]))/norm(X), 0f0, atol=1f-6)
@test isapprox(norm(X - AL.forward(AL.inverse(X))[1])/norm(X), 0f0, atol=1f-6)

###############################################################################
# Gradient Test

AL = AffineLayer(nx, ny, nc; logdet=true)
X = randn(Float32, nx, ny, nc, batchsize)
X0 = randn(Float32, nx, ny, nc, batchsize)
dX = X - X0; dX *= norm(X0)/norm(dX)

# Forward pass
Y = AL.forward(X)[1]

function loss(AL, X, Y)

    # Forward pass
    if AL.logdet
        Y_, lgdet = AL.forward(X)
    else
        Y_ = AL.forward(X)
    end

    # Residual and function value
    ΔY = Y_ - Y
    f = .5f0/batchsize*norm(ΔY)^2
    AL.logdet == true && (f -= lgdet)

    # Back propagation
    ΔX, X_ = AL.backward(ΔY./batchsize, Y_)

    # Check invertibility
    isapprox(norm(X - X_)/norm(X), 0f0, atol=1f-6)

    return f, ΔX, get_grads(AL)
end

# Gradient test for X
maxiter = 5
print("\nGradient test affine layer\n")
f0, ΔX = loss(AL, X0, Y)[1:2]
h = .1f0
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)
for j=1:maxiter
    f = loss(AL, X0 + h*dX, Y)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dX, ΔX))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test for parameters
logdet=false
AL = AffineLayer(nx, ny, nc; logdet=logdet)
θ = deepcopy(get_params(AL))
θ[2].data = randn(Float32, size(θ[2].data))
AL0 = AffineLayer(nx, ny, nc; logdet=logdet)
θ0 = get_params(AL0)
θ0[2].data = randn(Float32, size(θ0[2].data))
set_params!(AL0, θ0)
θ0 = deepcopy(θ0)
dθ = θ-θ0; dθ .*= norm.(θ0)./(norm.(dθ).+1f-10)

maxiter = 5
print("\nGradient test affine layer\n")
f0, ΔX, Δθ = loss(AL0, X, Y)
h = 0.1f0
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)
for j=1:maxiter
    set_params!(AL0, θ0+h*dθ)
    f = loss(AL0, X, Y)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(dθ, Δθ))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)

###################################################################################################
# Jacobian-related tests

# Gradient test

# Initialization
logdet=true
AL = AffineLayer(nx, ny, nc; logdet=logdet)
θ = deepcopy(get_params(AL))
AL0 = AffineLayer(nx, ny, nc; logdet=logdet)
θ0 = deepcopy(get_params(AL0))
X = randn(Float32, nx, ny, nc, batchsize)

# Perturbation (normalized)
dθ = θ-θ0
dθ[2].data = randn(Float32, size(dθ[2].data))
dθ .*= norm.(θ0)./(norm.(dθ).+1f-10)
dX = randn(Float32, nx, ny, nc, batchsize); dX *= norm(X)/norm(dX)

# Jacobian eval
logdet ? ((dY, Y, lgdet, GNdθ) = AL.jacobian(dX, dθ, X)) : ((dY, Y) = AL.jacobian(dX, dθ, X))

# Test
print("\nJacobian test\n")
h = 0.1f0
maxiter = 5
err9 = zeros(Float32, maxiter)
err10 = zeros(Float32, maxiter)
for j=1:maxiter
    set_params!(AL, θ+h*dθ)
    logdet ? ((Y_loc, _) = AL.forward(X+h*dX)) : (Y_loc = AL.forward(X+h*dX))
    err9[j] = norm(Y_loc - Y)
    err10[j] = norm(Y_loc - Y - h*dY)
    print(err9[j], "; ", err10[j], "\n")
    global h = h/2f0
end

@test isapprox(err9[end] / (err9[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err10[end] / (err10[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Adjoint test

set_params!(AL, θ)
logdet ? ((dY, Y, _, _) = AL.jacobian(dX, dθ, X)) : ((dY, Y) = AL.jacobian(dX, dθ, X))
dY_ = randn(Float32, size(dY))
logdet ? ((dX_, dθ_, _, _) = AL.adjointJacobian(dY_, Y)) : ((dX_, dθ_, _) = AL.adjointJacobian(dY_, Y))
a = dot(dY, dY_)
b = dot(dX, dX_)+dot(dθ, dθ_)
@test isapprox(a, b; rtol=1f-3)


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

export
    av_buffersink_read,
    av_buffersink_read_samples

function av_buffersink_read(ctx,buf)
    ccall((:av_buffersink_read,libavfilter),Cint,(Ptr{AVFilterContext},Ptr{Ptr{AVFilterBufferRef}}),ctx,buf)
end

function av_buffersink_read_samples(ctx,buf,nb_samples::Integer)
    ccall((:av_buffersink_read_samples,libavfilter),Cint,(Ptr{AVFilterContext},Ptr{Ptr{AVFilterBufferRef}},Cint),ctx,buf,nb_samples)
end


text:
"""
    glsol(sysg, sysf; poles = missing, sdeg = missing, mindeg = false, solgen = false, minreal = true, fast = true,
          atol = 0, atol1 = atol, atol2 = atol, rtol, ) -> (sysx, info, sysgen)

Determine for the descriptor systems `sysg = (Ag-λEg,Bg,Cg,Dg)` and
`sysf = (Af-λEf,Bf,Cf,Df)` with the transfer function matrices `G(λ)` and `F(λ)`,
respectively, the descriptor system `sysx` with the transfer function matrix `X(λ)`
such that `X(λ)` is the solution of the linear rational equation

    X(λ)G(λ) = F(λ) .      (1)

If `solgen = true`, the descriptor system `sysgen` is determined representing a generator of
all solutions of (1). Its transfer function matrix has the form `GEN(λ) = [ X0(λ); XN(λ) ]`,
such that any `X(λ)` can be generated as

    X(λ) = X0(λ) + Z(λ)*XN(λ) ,

where `X0(λ)` is a particular solution satisfying `X0(λ)G(λ) = F(λ)`,
`XN(λ)` is a proper rational left nullspace basis of `G(λ)` satisfying `XN(λ)G(λ) = 0`, and
`Z(λ)` is an arbitrary rational matrix with suitable dimensions. If `solgen = false`, `sysgen` is
set to `nothing`.

The call with

    glsol(sysgf, pf; poles = missing, sdeg = missing, mindeg = false, solgen = false, minreal = true, fast = true,
          atol = 0, atol1 = atol, atol2 = atol, rtol, ) -> (sysx, info, sysgen)

uses the compound descriptor system `sysgf = (A-λE,B,[Cg; Cf],[Dg; Df])`,
where `Cf` has `pf` rows, to define
the descriptor systems `sysg = (A-λE,B,Cg,Dg)` and `sysf = (A-λE,B,Cf,Df)`
(i.e., `Ag-λEg = Af-λEf = A-λE` and `Bg = Bf = B`).

The generator `sysgen` has a descriptor system realization `sysgen = (A0-λE0,B0, [C0; CN],[D0; DN])`,
which is usually not minimal  (unobservable and/or non-dynamic modes present), with

                 ( Ai-λEi    *       *    )
       A0-λE0  = (   0     Af-λEf    *    ) ,
                 (   0       0     Al-λEl )

                 ( *  )
           B0  = ( *  ),   ( C0 )  = ( C1 C2 C3 )
                 ( Bl )    ( CN )    ( 0  0  Cl )

with `El`, `Ef` and `Ai` invertible and upper triangular, `Ei` nillpotent
and upper triangular, and `DN` full column rank. The dimensions of the diagonal blocks of `A0-λE0` are
returned in the named tuple `info` as the components `info.nf`, `info.ninf` and `info.nl`, respectively.

A minimal order descriptor system realization of the proper basis `XN(λ)` is `(Al-λEl,Bl,Cl,DN)`,
where `Cl` and `DN` have `pr` columns (returned in `info.pr`), representing the dimension of the
left nullspace basis. The normal rank `nrank` of  `G(λ)` is returned in `info.nrank`.

If `mindeg = false`, the solution `sysx` is determined in the form `sysx = (A0+F*CN-λE0,B0+F*DN,C0,D0)`,
where the matrix `F = 0`, unless a nonzero stabilizing gain is used such that `Al+F*Bl-λEl` has stable eigenvalues.
The vector `poles` specified as a keyword argument, can be used to specify the desired eigenvalues
alternatively to or jointly with enforcing a desired stability degree `sdeg` of eigenvalues.
The dimension `nl` of `Al` is the number of freely assignable poles of the solution `X(λ)` and is returned in `info.nl`.
The eigenvalues of `Af-λEf` contain the finite zeros of `G(λ)`, while
the zeros  of `Ai-λEi` contain the infinite zeros of `G(λ)`.
The norm of the employed gain `F` is returned in `info.fnorm`. If `G(λ)` has infinite zeros, then the solution
`X(λ)` may have infinite poles. The integer vector `info.rdeg` contains the relative row degrees of
`X(λ)` (i.e., the numbers of integrators/delays needed to make each row of `X(λ)` proper).

If `mindeg = true`, a minimum degree solution is determined as `X(λ) = X0(λ) + Z(λ)XN(λ)`, where `Z(λ)` is
determined using order reduction based on a Type 2 minimum dynamic cover. This computation involves using
non-orthogonal transformations whose worst condition number is returned in `info.tcond`, in conjunction with
using feedback and feedforward gains, whose norms are returned in `info.fnorm`. High values of these quantities
indicate a potential loss of numerical stability of computations.

If `minreal = true`, the computed realization `sysx` is minimal.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively,
the absolute tolerance for the nonzero elements of `Ag`, `Bg`, `Cg`, `Dg`, `Af`, `Bf`, `Cf`, `Df`,
the absolute tolerance for the nonzero elements of `Eg` and `Ef`,
and the relative tolerance for the nonzero elements of
   `Ag`, `Bg`, `Cg`, `Dg`, `Af`, `Bf`, `Cf`, `Df`, `Eg` and `Ef`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the working machine epsilon
and `n` is the maximal order of the systems `sysg` and `sysf`. The keyword argument `atol` can be used
to simultaneously set `atol1 = atol`, `atol2 = atol`.

The rank determinations in the performed reductions
 are based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

_Method:_  The dual of method of [1] to solve rational systems is used.

_References:_

[1] A. Varga, "Computation of least order solutions of linear rational equations",
Proc. MTNS'04, Leuven, Belgium, 2004.
"""
function glsol(sysg::DescriptorStateSpace, sysf::DescriptorStateSpace; kwargs...)
   size(sysg,2) == size(sysf,2) || throw(DimensionMismatch("sysg and sysf must have the same number of inputs"))
   sysx, info1, sysgen = grsol(gdual(sysg), gdual(sysf); kwargs...)
   info = (nrank = info1.nrank, rdeg = info1.rdeg, tcond = info1.tcond, fnorm = info1.fnorm, pl = info1.mr, nl = info1.nr, nf = info1.nf, ninf = info1.ninf)
   return gdual(sysx,rev = true), info, sysgen === nothing ? sysgen : gdual(sysgen,rev = true)
end
function glsol(sysgf::DescriptorStateSpace, pf::Int; kwargs...)
   p = size(sysgf,1);
   (pf <= p && pf >= 0) || throw(DimensionMismatch("pf must be at most $p, got $pf"))
   sysx, info1, sysgen = grsol(gdual(sysgf,rev = true), pf; kwargs...)
   info = (nrank = info1.nrank, rdeg = info1.rdeg, tcond = info1.tcond, fnorm = info1.fnorm, pl = info1.mr, nl = info1.nr, nf = info1.nf, ninf = info1.ninf)
   return gdual(sysx,rev = true), info, sysgen === nothing ? sysgen : gdual(sysgen,rev = true)
end
"""
    grsol(sysg, sysf; poles = missing, sdeg = missing, mindeg = false, solgen = false, minreal = true, fast = true,
          atol = 0, atol1 = atol, atol2 = atol, rtol, ) -> (sysx, info, sysgen)

Determine for the descriptor systems `sysg = (Ag-λEg,Bg,Cg,Dg)` and
`sysf = (Af-λEf,Bf,Cf,Df)` with the transfer function matrices `G(λ)` and `F(λ)`,
respectively, the descriptor system `sysx` with the transfer function matrix `X(λ)`
such that `X(λ)` is the solution of the linear rational equation

    G(λ)X(λ) = F(λ) .      (1)

If `solgen = true`, the descriptor system `sysgen` is determined representing a generator of
all solutions of (1). Its transfer function matrix has the form `GEN(λ) = [ X0(λ) XN(λ) ]`,
such that any `X(λ)` can be generated as

    X(λ) = X0(λ) + XN(λ)*Z(λ) ,

where `X0(λ)` is a particular solution satisfying `G(λ)X0(λ) = F(λ)`,
`XN(λ)` is a proper rational right nullspace basis of `G(λ)` satisfying `G(λ)XN(λ) = 0`, and
`Z(λ)` is an arbitrary rational matrix with suitable dimensions. If `solgen = false`, `sysgen` is
set to `nothing`.

The call with

    grsol(sysgf, mf; poles = missing, sdeg = missing, mindeg = false, solgen = false, minreal = true, fast = true,
          atol = 0, atol1 = atol, atol2 = atol, rtol, ) -> (sysx, info, sysgen)

uses the compound descriptor system `sysgf = (A-λE,[Bg Bf],C,[Dg Df])`,
where `Bf` has `mf` columns, to define
the descriptor systems `sysg = (A-λE,Bg,C,Dg)` and `sysf = (A-λE,Bf,C,Df)`
(i.e., `Ag-λEg = Af-λEf = A-λE` and `Cg = Cf = C`).

The generator `sysgen` has a descriptor system realization `sysgen = (A0-λE0,[B0 BN],C0,[D0 DN])`,
which is usually not minimal  (uncontrollable and/or non-dynamic modes present), with

                   ( Ar-λEr    *       *    )
       A0-λE0    = (   0     Af-λEf    *    ) ,
                   (   0       0     Ai-λEi )

                   ( B1 | Br )
       [B0 | BN] = ( B2 | 0  ),  Cg  =   ( Cr   *    *  ) ,
                   ( B3 | 0  )

with `Er`, `Ef` and `Ai` invertible and upper triangular, `Ei` nillpotent
and upper triangular, and `DN` full row rank. The dimensions of the diagonal blocks of `A0-λE0` are
returned in the named tuple `info` as the components `info.nr`, `info.nf`, and `info.ninf`, respectively.

A minimal order descriptor system realization of the proper basis `XN(λ)` is `(Ar-λEr,Br,Cr,DN)`,
where `Br` and `DN` have `mr` columns (returned in `info.mr`), representing the dimension of the
right nullspace basis. The normal rank `nrank` of  `G(λ)` is returned in `info.nrank`.

If `mindeg = false`, the solution `sysx` is determined in the form `sysx = (A0+BN*F-λE0,B0,C0+DN*F,D0)`,
where the matrix `F = 0`, unless a nonzero stabilizing gain is used such that `Ar+Br*F-λEr` has stable eigenvalues.
The vector `poles` specified as a keyword argument, can be used to specify the desired eigenvalues
alternatively to or jointly with enforcing a desired stability degree `sdeg` of eigenvalues.
The dimension `nr` of `Ar` is the number of freely assignable poles of the solution `X(λ)` and is returned in `info.nr`.
The eigenvalues of `Af-λEf` contain the finite zeros of `G(λ)`, while
the zeros  of `Ai-λEi` contain the infinite zeros of `G(λ)`.
The norm of the employed gain `F` is returned in `info.fnorm`. If `G(λ)` has infinite zeros, then the solution
`X(λ)` may have infinite poles. The integer vector `info.rdeg` contains the relative column degrees of
`X(λ)` (i.e., the numbers of integrators/delays needed to make each column of `X(λ)` proper).

If `mindeg = true`, a minimum degree solution is determined as `X(λ) = X0(λ) + XN(λ)*Z(λ)`, where `Z(λ)` is
determined using order reduction based on a Type 2 minimum dynamic cover. This computation involves using
non-orthogonal transformations whose worst condition number is returned in `info.tcond`, in conjunction with
using feedback and feedforward gains, whose norms are returned in `info.fnorm`. High values of these quantities
indicate a potential loss of numerical stability of computations.

If `minreal = true`, the computed realization `sysx` is minimal.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively,
the absolute tolerance for the nonzero elements of `Ag`, `Bg`, `Cg`, `Dg`, `Af`, `Bf`, `Cf`, `Df`,
the absolute tolerance for the nonzero elements of `Eg` and `Ef`,
and the relative tolerance for the nonzero elements of
   `Ag`, `Bg`, `Cg`, `Dg`, `Af`, `Bf`, `Cf`, `Df`, `Eg` and `Ef`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the working machine epsilon
and `n` is the maximal order of the systems `sysg` and `sysf`. The keyword argument `atol` can be used
to simultaneously set `atol1 = atol`, `atol2 = atol`.

The rank determinations in the performed reductions
 are based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

_Method:_  The method of [1] to solve rational systems is used.

_References:_

[1] A. Varga, "Computation of least order solutions of linear rational equations",
Proc. MTNS'04, Leuven, Belgium, 2004.
"""
function grsol(sysg::DescriptorStateSpace{T1}, sysf::DescriptorStateSpace{T2};
               poles::Union{AbstractVector,Missing} = missing, sdeg::Union{Real,Missing} = missing,
               atol::Real = zero(real(T1)), atol1::Real = atol, atol2::Real = atol,
               rtol::Real = ((size(sysg.A,1)+1)*eps(real(float(one(T1)))))*iszero(max(atol1,atol2)),
               fast::Bool = true, mindeg::Bool = false, solgen::Bool = false, minreal::Bool = true) where {T1,T2}
   size(sysg,1) == size(sysf,1) || throw(DimensionMismatch("sysg and sysf must have the same number of outputs"))
   return grsol(gir([sysg sysf]; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol), size(sysf,2);
   poles = poles, sdeg = sdeg, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, mindeg = mindeg, solgen = solgen, minreal = minreal)
end
function grsol(sysgf::DescriptorStateSpace{T}, mf::Int;
               poles::Union{AbstractVector,Missing} = missing, sdeg::Union{Real,Missing} = missing,
               atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
               rtol::Real = ((size(sysgf.A,1)+1)*eps(real(float(one(T)))))*iszero(max(atol1,atol2)),
               fast::Bool = true, mindeg::Bool = false, solgen::Bool = false, minreal::Bool = true) where T

   T1 = T <: BlasFloat ? T : promote_type(Float64,T)

   p, mgf = size(sysgf);
   (mf <= mgf && mf >= 0) || throw(DimensionMismatch("mf must be at most $mgf, got $mf"))
   m = mgf-mf
   Ts = sysgf.Ts;

   ONE = one(T1)

   disc = !iszero(Ts)
   sdeg_nomissing = !ismissing(sdeg)
   poles_nomissing = !ismissing(poles)

   stabilize = sdeg_nomissing || poles_nomissing

   # set default values of sdeg if poles = missing
   # sdeg_nomissing || (sdeg = disc ? real(T1)(0.95) : real(T1)(-0.05))

   # sort desired eigenvalues
   if poles_nomissing
      tempc = poles[imag.(poles) .> 0]
      if !isempty(tempc)
         tempc1 = conj(poles[imag.(poles) .< 0])
         isequal(tempc[sortperm(real(tempc))],tempc1[sortperm(real(tempc1))]) ||
                 error("poles must be a self-conjugated complex vector")
      end
      # check that all eigenvalues are inside of the stability region
      sdeg_nomissing && ( ((disc && any(abs.(poles) .> sdeg) )  || (!disc && any(real.(poles) .> sdeg)))  &&
            error("The elements of poles must lie in the stability region of interest") )
   end

   At, Et, Bt, Ct, Dt = dssdata(T1,sysgf);

   n = size(At,1)

   Bf = view(Bt,:,m+1:mgf)
   Df = view(Dt,:,m+1:mgf)
   Bg = view(Bt,:,1:m)
   Dg = view(Dt,:,1:m)

   # compute the Kronecker-like form of the system matrix of G to obtain
   #                 ( Br  Ar - s Er      *             *           *     )
   #   Q'*(A-sE)*Z = (  0     0      Afin - s Efin      *           *     )
   #                 (  0     0           0        Ainf - s Einf    *     )
   #                 (  0     0           0            0         Al - s El)

   M, N, Q, Z, νr, μr, νi, nf, νl, μl = sklf(At, Et, Bg, Ct, Dg,
                                             finite_infinite = true, fast = fast, ut = true,
                                             atol1 = atol1, atol2 = atol2, rtol = rtol)

   #                        ( B1 )
   # compute Q'*[-Bf;-Df] = ( B2 )
   #                        ( B3 )
   #                        ( B4 )
   f = Q'*[-Bf;-Df];

   # determine the orders of Ar, Ai and Al
   nr = sum(νr); ninf = sum(νi); nl = sum(μl);
   # determine the orders of regular part and of the maximal invertible part
   nreg = nf+ninf; ninv = nr+nreg+nl;
   mr = n+m-ninv  # number of columns of Br

   # check compatibility condition, i.e., B4 = 0.
   tola = max(atol1,atol2)
   iszero(tola) &&
      (tola = n*eps(real(T1))*max(opnorm(Et,1), opnorm(At,1), opnorm(Bt,1), opnorm(Ct,Inf)))
      # println("p = $p n = $n m = $m mf = $mf")
      # println("νl = $νl  nr = $nr  nreg = $nreg n = $n  p = $p")
      # println("sum(νl) = $(sum(νl))")
      # println("sum(νl) = $(sum(νl))")
      # println("f = $f")
      # println("f[nr+nreg+1:end,:] = $(f[nr+nreg+1:end,:])")
   sum(νl) > 0 && mf > 0 && maximum(abs.(f[nr+nreg+1:end,:])) >= tola &&
         error("System not compatible")

   # form X0 = (A0-sE0,B0,C0,D0) and XN = (Ar-sEr,Br,Cr,DN), where
   #
   #                ( Ar - s*Er     *           *        )        ( B1 )        ( Br )
   #       A0-sE0 = (    0       Af - s*Ef      *        ) , B0 = ( B2 ) , BN = ( 0  )
   #                (    0          0        Ai - s*Einf )        ( B3 )        ( 0  )
   #
   #           C0 = ( Cr  Cf Ci );  D0 = 0, with
   #
   #             [ 0  Im ]*Z := [ DN Cr  Cf Ci Cl ];
   n0 = nr+nreg;
   i1 = 1:nr;
   A0 = M[1:n0,mr+1:mr+n0]; E0 = N[1:n0,mr+1:mr+n0];
   B0 = f[1:n0,:];    C0 = Z[n+1:n+m,mr+1:mr+n0]; D0 = zeros(T1,m,mf)
   Br = view(M,i1,1:mr); BN = [Br; zeros(T1,nreg,mr)]; DN = view(Z,n+1:n+m,1:mr); # CN = C0
   Ar = view(A0,i1,i1)
   Er = view(E0,i1,i1)
   Cr = view(C0,:,i1)

   rdeg = zeros(Int,mf)
   if length(νi) > 1
      i2 = nr+nf+1:n0
      # compute relative column degrees as the number of infinite poles
      # (in excess with 1 with respect to the controllable infinite eigenvalues)
      for i = 1:mf
          _, _, _, ni, = sklf_right!(copy(E0[i2,i2]), copy(A0[i2,i2]), copy(B0[i2,i:i]), copy(C0[:,i2]);
                                    fast = fast, atol1 = atol1, atol2 = atol2, atol3 = atol1, rtol = rtol,
                                    withQ = false, withZ = false)
          ni > 1  && (rdeg[i] = ni-1)
      end
   end

   if mindeg
      # form the pair of generators if requested
      # copy A0, E0 and C0, to ensure are not later modified (due to bindings)
      solgen ? sysgen = dss(copy(A0),E0,[B0 BN],copy(C0),[D0 DN],Ts = Ts) : sysgen = nothing

      # block-diagonalize E0 to allow working on the proper part
      # exploit Er is upper triangular
      i2 = nr+1:n0;
      Y = -E0[i1,i2]
      ldiv!(UpperTriangular(Er),Y)
      mul!(view(A0,i1,i2), Ar, Y, ONE, ONE)
      mul!(view(C0,:,i2), Cr, Y, ONE, ONE)
      # Y = -E0[i1,i1]\E0[i1,i2];
      # A0[i1,i2] = A0[i1,i2]+A0[i1,i1]*Y;
      # C0[:,i2] = C0[:,i2]+C0[:,i1]*Y;

      # form a pair of generators for the proper part
      sysr = dss(Ar, Er, [B0[i1,:] A0[i1,i2] Br], Cr, [D0 C0[:,i2] DN], Ts = Ts);

      # compute minimum order for the proper part
      sysr, _, info2 = grmcover2(sysr,mf+nreg, fast = fast, atol1 = atol1, atol2 = atol2)
      aa, ee, bb, cc, dd = dssdata(sysr); na = size(aa,1); ee == I && (ee = Matrix{T1}(I,na,na))
      i3 = mf+1:mf+nreg;
      A0 = [aa bb[:,i3]; zeros(T1,nreg,na) A0[i2,i2]]; E0 = blockdiag(ee,E0[i2,i2])
      B0 = [ bb[:,1:mf]; B0[i2,:]]; C0 = [cc dd[:,i3]]; D0 = dd[:,1:mf];
      tcond = info2.tcond;
      fnorm = max(info2.fnorm,info2.gnorm);
   else
      fnorm = 0
      tcond = 1
      if stabilize
         # make spurious poles stable
         if nr > 0
            i1 = 1:nr;
            F = saloc(Ar, Er, Br; evals = poles, sdeg = sdeg, disc = disc, fast = fast, atol1 = atol1, atol2 = atol2, atol3 = atol1)[1]
            #F = saloc(A0[i1,i1], E0[i1,i1], Br; evals = poles, sdeg = sdeg, disc = disc, fast = fast, atol1 = atol1, atol2 = atol2, atol3 = atol1)[1]
            #A0[i1,i1] = A0[i1,i1]+Br*F; C0[:,i1] = C0[:,i1]+DN*F
            mul!(Ar, Br, F, ONE, ONE)
            mul!(Cr, DN, F, ONE, ONE)
         end
         fnorm = norm(F)
      end
      # form a pair of generators
      solgen ? sysgen = dss(A0, E0, [B0 BN], C0, [D0 DN], Ts = Ts) : sysgen = nothing
   end

   # eliminate possible uncontrollable eigenvalues and non-dynamic modes
   minreal ? sysx = gminreal(dss(A0, E0, B0, C0, D0, Ts = Ts), obs = false, fast = fast, atol1 = atol1, atol2 = atol2) :
             sysx = dss(A0, E0, B0, C0, D0, Ts = Ts)
   # sysx = gss2ss(sysx,tol,'triu');
   info = (nrank = ninv-n, rdeg = rdeg, tcond = tcond, fnorm = fnorm, mr = mr, nr = nr, nf = nf, ninf = ninf)
   return sysx, info, sysgen
# end GRSOL
end


text:
"""
MulAddAccumulate(μ, A, B)

represents the vector satisfying v[k+1] == A[k]*v[k]+B[k] with v[1] == μ
"""
mutable struct MulAddAccumulate{T} <: AbstractCachedVector{T}
    data::Vector{T}
    A::AbstractVector
    B::AbstractVector
    datasize::Tuple{Int}
    function MulAddAccumulate{T}(data::Vector{T}, A::AbstractVector, B::AbstractVector, datasize::Tuple{Int}) where T
        length(A) == length(B) || throw(ArgumentError("lengths must match"))
        new{T}(data, A, B, datasize)
    end
end

size(M::MulAddAccumulate) = size(M.A)

MulAddAccumulate(data::Vector{T}, A::AbstractVector, B::AbstractVector, datasize::Tuple{Int}) where T =
    MulAddAccumulate{T}(data, A, B, datasize)
MulAddAccumulate(μ, A, B) = MulAddAccumulate([μ], A, B, (1,))
MulAddAccumulate(A, B) = MulAddAccumulate(A[1]+B[1], A, B)

function LazyArrays.cache_filldata!(K::MulAddAccumulate, inds)
    A,B = K.A,K.B
    @inbounds for k in inds
        K.data[k] = muladd(A[k], K.data[k-1], B[k])
    end
end

@simplify function *(D::Derivative, P::SemiclassicalJacobi)
    Q = SemiclassicalJacobi(P.t, P.a+1,P.b+1,P.c+1,P)
    A,B,_ = recurrencecoefficients(P)
    α,β,_ = recurrencecoefficients(Q)

    d = AccumulateAbstractVector(*, A ./ Vcat(1,α))
    v1 = AccumulateAbstractVector(+, B ./ A)
    v2 = MulAddAccumulate(Vcat(0,0,α[2:∞]) ./ α, Vcat(0,β ./ α) ./ α);
    v3 = AccumulateAbstractVector(*, Vcat(A[1]A[2], A[3:∞] ./ α))
    Q * (_BandedMatrix(Vcat(((1:∞) .* d)', (((1:∞) .* (v1 .+ B[2:end]./A[2:end]) .- (2:∞) .* (α .* v2 .+ β ./ α)) .* v3)'), ∞, 2,-1))'
end

@simplify function *(D::Derivative, wQ::Weighted{<:Any,<:SemiclassicalJacobi})
    Q = wQ.P
    P = SemiclassicalJacobi(Q.t, Q.a-1,Q.b-1,Q.c-1)
    Weighted(P) * ((-sum(orthogonalityweight(Q))/sum(orthogonalityweight(P))) * (Q \ (D * P))')
end


text:
using TimeSeriesInterface
using Test
using Dates
using Statistics

include("timeseries.jl")
include("forecast.jl")
include("models.jl")


text:
function get_free_state_parameters(obj::AbstractGradientConverter)
    return jcall(obj, "getFreeStateParameters", jint, ())
end


text:
function reshop_ovf(mdl, ovf)
    # 0-based vs 1-based
    ovf_vidx = ovf.vidx-1
    argsC = Vector{Cint}(ovf.args .- 1)

    ovf_def = emp_ovf(mdl, ovf.name, ovf_vidx, argsC)

    for (k, v) in ovf.params
        emp_ovf_param(ovf_def, k, v)
    end

    emp_ovf_check(ovf_def)
end


text:
module VoxelwiseConstantMapping

using AbbreviatedTypes
using HCubature
using Statistics: mean

include("vxlwise_const.jl")
include("integral.jl")

end # module


text:
type Person
    firstname::String
    lastname::String
    sex::Char
    age::Float64
    children::Array{String, 1}
end

p1 = Person("Alan", "Bates", 'M', 45.5, ["Jeff", "Stephan"])
#> Person("Alan", "Bates",'M',45.5,String["Jeff","Stephan"])

people = Person[] #> 0-element Array{Person,1}
push!(people, p1)
push!(people, Person("Julia", "Smith", 'F', 27, ["Viral"]))

show(people)
#> [Person("Alan", "Bates", 'M',45.5,String["Jeff","Stephan"]),
# Person("Julia", "Smith",'F',27.0,String["Viral"])]

fullname(p::Person) = "$(p.firstname) $(p.lastname)"
fullname(p::Person) = string(p.firstname, " ", p.lastname)
print(fullname(p1)) #> Alan Bates

type Family
     name::String
     members::Array{String, 1}
     big::Bool
     Family(name::String) = new(name, String[], false)
     Family(name::String, members) = new(name, members, length(members) > 4)
end

fam = Family("Bates-Smith", ["Alan", "Julia", "Jeff", "Stephan", "Viral"])
#> Family("Bates-Smith",String["Alan","Julia","Jeff","Stephan","Viral"],true)


text:
baremodule ShowCode

export @sc_intel,
    @sc_ircode,
    @sc_llvm,
    @sc_native,
    post_godbolt,
    sc_intel,
    sc_ircode,
    sc_llvm,
    sc_native

#=
macro sc_sexpr end
macro sc_lowered end
macro sc_typed end
=#
macro sc_ircode end
macro sc_llvm end
macro sc_native end
macro sc_intel end

function sc_ircode end
function sc_llvm end
function sc_native end
function sc_intel end

function from_llvm end
function from_native end

function post_godbolt end

module Implementations

using ..ShowCode: ShowCode
import ..ShowCode:
    @sc_intel,
    @sc_ircode,
    @sc_llvm,
    @sc_native,
    post_godbolt,
    sc_intel,
    sc_ircode,
    sc_llvm,
    sc_native

import AbstractTrees
import DefaultApplication
import HTTP
import JSON
using Accessors: @set
using Base64: base64encode
using InteractiveUtils:
    InteractiveUtils,
    code_llvm,
    code_lowered,
    code_native,
    code_typed,
    gen_call_with_extracted_types_and_kwargs,
    print_llvm,
    print_native
using Requires: @require
using UnPack: @unpack

include("core.jl")
include("code_ircode.jl")
include("ircode.jl")
include("llvm.jl")
include("native.jl")
include("godbolt.jl")

function __init__()
    @require Cthulhu="f68482b8-f384-11e8-15f7-abe071a5a75f" include("cthulhu.jl")
end

end

const CONFIG = Implementations.ShowCodeConfig()

Implementations.finalize_module()

end


text:
using INMET
using Test

@testset "INMET.jl" begin
  @testset "stations" begin
    # automatic stations
    df = INMET.stations()
    @test all(isequal("Automatica"), df.TP_ESTACAO)

    # manual stations
    df = INMET.stations(:manual)
    @test all(isequal("Convencional"), df.TP_ESTACAO)
  end

  @testset "series" begin
    df = INMET.series(:A301, Date(2021,1,1), Date(2021,7,31))
    @test all(isequal("A301"), df.CD_ESTACAO)
  end

  @testset "on" begin
    dfday  = INMET.on(Date(2021,7,1))
    dfhour = INMET.on(DateTime(2021,7,1,1))
    @test size(dfday, 1) > size(dfhour, 1)
    @test all(isequal("2021-07-01"), dfday.DT_MEDICAO)
    @test all(isequal("2021-07-01"), dfhour.DT_MEDICAO)
    @test all(isequal("0100"), dfhour.HR_MEDICAO)
  end
end


text:
# This should stay as the first method because it's used in a test
# (or change the test)
function checkname(fdef::Expr, name)
    fproto = fdef.args[1]
    (fdef.head === :where || fdef.head == :(::)) && return checkname(fproto, name)
    fdef.head === :call || return false
    if fproto isa Expr
        # A metaprogramming-generated function
        fproto.head === :$ && return true   # uncheckable, let's assume all is well
        # Is the check below redundant?
        fproto.head === :. || return false
        # E.g. `function Mod.bar.foo(a, b)`
        return checkname(fproto.args[end], name)
    end
    isa(fproto, Symbol) || isa(fproto, QuoteNode) || isa(fproto, Expr) || return false
    return checkname(fproto, name)
end
checkname(fname::Symbol, name::Symbol) = begin
    fname === name && return true
    startswith(string(name), string('#', fname, '#')) && return true
    return false
end
checkname(fname::Symbol, ::Nothing) = true
checkname(fname::QuoteNode, name) = checkname(fname.value, name)

function isfuncexpr(ex, name=nothing)
    # Strip any macros that wrap the method definition
    while ex isa Expr && ex.head === :macrocall && length(ex.args) == 3
        ex = ex.args[3]
    end
    isa(ex, Expr) || return false
    if ex.head === :function || ex.head === :(=)
        return checkname(ex.args[1], name)
    end
    return false
end

function linerange(def::Expr)
    start, haslinestart = findline(def, identity)
    stop, haslinestop  = findline(def, Iterators.reverse)
    (haslinestart & haslinestop) && return start:stop
    return nothing
end
linerange(arg) = linerange(convert(Expr, arg))  # Handle Revise's RelocatableExpr

function findline(ex, order)
    ex.head === :line && return ex.args[1], true
    for a in order(ex.args)
        a isa LineNumberNode && return a.line, true
        if a isa Expr
            ln, hasline = findline(a, order)
            hasline && return ln, true
        end
    end
    return 0, false
end

fileline(lin::LineInfoNode)   = String(lin.file), lin.line
fileline(lnn::LineNumberNode) = String(lnn.file), lnn.line

# This is piracy, but it's not ambiguous in terms of what it should do
Base.convert(::Type{LineNumberNode}, lin::LineInfoNode) = LineNumberNode(lin.line, lin.file)

# This regex matches the pseudo-file name of a REPL history entry.
const rREPL = r"^REPL\[(\d+)\]$"

"""
    src = src_from_file_or_REPL(origin::AbstractString, repl = Base.active_repl)

Read the source for a function from `origin`, which is either the name of a file
or "REPL[\$i]", where `i` is an integer specifying the particular history entry.
Methods defined at the REPL use strings of this form in their `file` field.

If you happen to have a file where the name matches `REPL[\$i]`, first pass it through
`abspath`.
"""
function src_from_file_or_REPL(origin::AbstractString, args...)
    # This Varargs design prevents an unnecessary error when Base.active_repl is undefined
    # and `origin` does not match "REPL[$i]"
    m = match(rREPL, origin)
    if m !== nothing
        return src_from_REPL(m.captures[1], args...)
    end
    isfile(origin) || return nothing
    return read(origin, String)
end

function src_from_REPL(origin::AbstractString, repl = Base.active_repl)
    hist_idx = parse(Int, origin)
    hp = repl.interface.modes[1].hist
    return hp.history[hp.start_idx+hist_idx]
end

function basepath(id::PkgId)
    id.name ∈ ("Main", "Base", "Core") && return ""
    loc = Base.locate_package(id)
    loc === nothing && return ""
    return dirname(dirname(loc))
end

"""
    path = maybe_fix_path(path)

Return a normalized, absolute path for a source file `path`.
"""
function maybe_fix_path(file)
    if !isabspath(file)
        # This may be a Base or Core method
        newfile = Base.find_source_file(file)
        if isa(newfile, AbstractString)
            file = normpath(newfile)
        end
    end
    return maybe_fixup_stdlib_path(file)
end

safe_isfile(x) = try isfile(x); catch; false end
const BUILDBOT_STDLIB_PATH = dirname(abspath(joinpath(String((@which uuid1()).file), "..", "..", "..")))
replace_buildbot_stdlibpath(str::String) = replace(str, BUILDBOT_STDLIB_PATH => Sys.STDLIB)
"""
    path = maybe_fixup_stdlib_path(path::String)

Return `path` corrected for julia issue [#26314](https://github.com/JuliaLang/julia/issues/26314) if applicable.
Otherwise, return the input `path` unchanged.

Due to the issue mentioned above, location info for methods defined one of Julia's standard libraries
are, for non source Julia builds, given as absolute paths on the worker that built the `julia` executable.
This function corrects such a path to instead refer to the local path on the users drive.
"""
function maybe_fixup_stdlib_path(path)
    if !safe_isfile(path)
        maybe_stdlib_path = replace_buildbot_stdlibpath(path)
        safe_isfile(maybe_stdlib_path) && return maybe_stdlib_path
    end
    return path
end

function postpath(filename, pre)
    idx = findfirst(pre, filename)
    idx === nothing && error(pre, " not found in ", filename)
    post = filename[first(idx) + length(pre) : end]
    post[1:1] == Base.Filesystem.path_separator && return post[2:end]
    return post
end

if Base.VERSION < v"1.1"
    function splitpath(p::String)
        # splitpath became available with Julia 1.1
        # Implementation copied from Base except doesn't handle the drive
        out = String[]
        isempty(p) && (pushfirst!(out,p))  # "" means the current directory.
        while !isempty(p)
            dir, base = _splitdir_nodrive(p)
            dir == p && (pushfirst!(out, dir); break)  # Reached root node.
            if !isempty(base)  # Skip trailing '/' in basename
                pushfirst!(out, base)
            end
            p = dir
        end
        return out
    end
    splitpath(p::AbstractString) = splitpath(String(p))

    _splitdir_nodrive(path::String) = _splitdir_nodrive("", path)
    function _splitdir_nodrive(a::String, b::String)
        m = match(Base.Filesystem.path_dir_splitter,b)
        m === nothing && return (a,b)
        a = string(a, isempty(m.captures[1]) ? m.captures[2][1] : m.captures[1])
        a, String(m.captures[3])
    end
end

# Robust across Julia versions
getpkgid(project::AbstractString, libname) = getpkgid(Base.project_deps_get(project, libname), libname)
getpkgid(id::PkgId, libname) = id
getpkgid(uuid::UUID, libname) = PkgId(uuid, libname)