PatternVectors.jl

About this Package

• ZeroPattern: convenient representation for vectors of the form:

[0,0,0,0,0,0...]

• FillPattern: convenient representation for vectors of the form:

[a,a,a,a,a,a...]

• InitialValuePattern: convenient representation for vectors of the form:

[a,b,b,b,b,b...]

• FinalValuePattern: convenient representation for vectors of the form:

[a,a,a,a,a,a...,b]

• PaddedFillPattern: convenient representation for vectors of the form:

[x,a,a,a,a,a...,y]

• EvenOddPattern: convenient representation for vectors of the form:

[a,b,a,b,a,b...]

• PaddedEvenOddPattern: convenient representation for vectors of the form:

[x,a,b,a,b,a,b...,y]

Defining a new pattern is easy and fully supported by the library. The module is standalone.

Performances Comparison

Simple multiplication

using PatternVectors, BenchmarkTools
n=10_000
x=PatternVector(n,EvenOddPattern(0.2,2.3))
y=randn(n)
x_c=collect(x)
@btime @. $x*$y;
@btime @. $x_c*$y;

  2.632 μs (3 allocations: 78.20 KiB)
  3.345 μs (3 allocations: 78.20 KiB)

Flipping sign based on index and sum

using PatternVectors, BenchmarkTools, LinearAlgebra
n=10_000
function f_std_scalar(f_x)
	N=length(f_x)
	sum_z=zero(eltype(f_x))
	for i in 1:N
		w=ifelse(isodd(i),1,-1)
		sum_z+=@views @inbounds w*f_x[i]
	end
	return sum_z
end

function f_std_scalar_2(f_x)
	N=length(f_x)
	sum_z=zero(eltype(f_x))
	for i in 1:N
		if(isodd(i))
			sum_z+=@views @inbounds f_x[i]
		else
			sum_z-=@views @inbounds f_x[i]
		end
	end
	return sum_z
end

function f_std_vec(f_x)
	N=length(f_x)
	idx=1:N
	W=@. ifelse(isodd(idx),1,-1)
	return sum(W.*f_x)
end

function f_std_vec_linear_algebra(f_x)
	N=length(f_x)
	idx=1:N
	W=@. ifelse(isodd(idx),1,-1)
	return dot(W,f_x)
end

function f_apv(f_x)
	N=length(f_x)
	W=PatternVector(n,EvenOddPattern(1,-1))
	return sum(W.*f_x)
end
x=randn(n)
f_x=@. sin(x)+x*cos(x)
@btime f_std_scalar($f_x);
@btime f_std_scalar_2($f_x);
@btime f_std_vec($f_x);
@btime f_std_vec_linear_algebra($f_x);
@btime f_apv($f_x);

  9.267 μs (0 allocations: 0 bytes)
  9.267 μs (0 allocations: 0 bytes)
  6.870 μs (6 allocations: 156.39 KiB)
  12.985 μs (2 allocations: 78.16 KiB)
  3.608 μs (3 allocations: 78.20 KiB)

Simpson Integration

using PatternVectors, BenchmarkTools,LinearAlgebra
n2=10_000
function f_simpson_std_scalar(f_x)
	N=length(f_x)
	sum_z=zero(eltype(f_x))
	for i in 1:N
		w=ifelse(i==1,1/3,ifelse(i==N,1/3,4/3))
		sum_z+=@views @inbounds w*f_x[i]
	end
	return sum_z
end

function f_simpson_std_vec(f_x)
	N=length(f_x)
	idx=1:N
	W=@. ifelse(idx==1,1/3,ifelse(idx==N,1/3,4/3))
	return sum(W.*f_x)
end

function f_simpson_std_vec_linear_algebra(f_x)
	N=length(f_x)
	idx=1:N
	W=@. ifelse(idx==1,1/3,ifelse(idx==N,1/3,4/3))
	return dot(W,f_x)
end

function f_simpson_apv(f_x)
	N=length(f_x)
	W=PatternVector(N,PaddedEvenOddPattern(1/3,4/3,4/3,1/3))
	return sum(W.*f_x)
end

function f_simpson_apv_linear_algebra(f_x)
	N=length(f_x)
	W=PatternVector(N,PaddedEvenOddPattern(1/3,4/3,4/3,1/3))
	return dot(W,f_x)
end

x2=randn(n2)
f_x2=@. sin(x2)+x2*cos(x2)
@btime f_simpson_std_scalar($f_x2);
@btime f_simpson_std_vec($f_x2);
@btime f_simpson_std_vec_linear_algebra($f_x2);
@btime f_simpson_apv($f_x2);
@btime f_simpson_apv_linear_algebra($f_x2);

  9.267 μs (0 allocations: 0 bytes)
  7.051 μs (6 allocations: 156.39 KiB)
  4.224 μs (3 allocations: 78.20 KiB)
  4.065 μs (3 allocations: 78.20 KiB)
  9.277 μs (0 allocations: 0 bytes)

To be noticed the performance improvements thanks to the usage of PatternVector, and to be noticed that the first function proposed is not compatible with the CUDA.jl stack.

New Patterns Definition

• Define a container (namely MyNewPattern{T}) inheriting from AbstractPattern{T} where T is the data type you are going to store in the array.

• Define a constructor for it.

• Implement the interface:

  • pattern_minimum_size: define the minimum size for your type

  • getindex_pattern: define the logic to extract scalars

  • getindex_pattern_range: define the logic to extract ranges

  • materialize_pattern: define the way to compute broadcasted functions with the new pattern

In case you need to mix your pattern with other existing patterns, you will have to specify we is the more general pattern able to store both:

• determine_mixed_pattern(::Type{T}, ::Type{V}) where {T <: MyNewPattern{M}, V <: MyOldPattern{N}} where {M, N} = MyWinningPattern{promote_type(M, N)}

And you are ready to go!

In case you want to use your new pattern in AD applications you will have to provide two additional implementations:

• ChainRulesCore.rrule(::Type{MyNewPattern}, args...): the rrule for the constructor of your newly defined pattern.

• pattern_to_vector_pullback: the rrule to convert from pattern to array.

For more details have a look at the already implemented types.