Author Archives: Picaud Vincent

Julia dispatching, @enum versus Type

Re-posted from: https://pixorblog.wordpress.com/2018/02/23/julia-dispatch-enum-vs-type-comparison/

1 Context

Imagine the following scenario: a function func must efficiently calls a (short) subroutine belonging to predefined set of specialized functions. The chosen specialization is done thanks to a tag argument. Something like:

specialization(x,tag_1) = specialization 1
specialization(x,tag_2) = specialization 2
etc.

function func(tag) 
  # some code
  for i in range
     specialization(x,tag)
  end 
  # some code
end

Like specialization is called often, an efficient solution to perform the dispatching is desired.

Here I check two Julia solutions, one using @enum, the other using Type.

2 `@enum` based approach

Coming from the C++ world, I will start from the following C++ code to propose a Julia solution:

#include <iostream>
#include <type_traits>

enum Tags { A, B, C }; // I do not used enum class on purpose, 
                       // to have A,B,C... as global identifiers

auto specialized(int x, std::integral_constant<Tags, A>) 
{ return 1 * x; }
auto specialized(int x, std::integral_constant<Tags, B>) 
{ return 2 * x; }
auto specialized(int x, std::integral_constant<Tags, C>) 
{ return 3 * x; }

template <Tags TAG>
constexpr auto Tags_v = std::integral_constant<Tags,TAG>();

template <Tags TAG>
auto func(std::integral_constant<Tags, TAG> tag)
{
  return specialized(10, tag); // no run-time penalty
}

int main() 
{ 
  std::cout << func(Tags_v<B>); 
}

In C++, there is no run-time penalty because the right specialization is chosen at compile-time. In basic situations like this one, the compiler even inline everything.

Julia provides an @enum macro. To mimic the C++ std::integral_constant we use the Type{Val{.}} trick.

@enum Tags A B C

specialized(x::Int,::Type{Val{A}}) = 1*x
specialized(x::Int,::Type{Val{B}}) = 2*x
specialized(x::Int,::Type{Val{C}}) = 3*x

func(::Type{Val{e}}) where {e} = specialized(10,Val{e})

This work as exepected:

func(Val{B})

The generated assembly code confirms that everything is inlined as in C++.

@code_native func(Val{B})

	.text
Filename: none
	pushq	%rbp
	movq	%rsp, %rbp
Source line: 1
	movl	$20, %eax
	popq	%rbp
	retq
	nopl	(%rax,%rax)

Comparison with C++:

The C++ code is 13 lines long whereas the Julia one is only 5 lines long. The call site syntax is nearly the same:

func(Tags_v<A>);

versus

func(Val{B})

However, concerning the Julia solution, I have a minor regret. AFAIK in Julia we can not filter argument e with something like

func(::Type{Val{e}}) where {e::Tags} = specialized(10,Val{e}) # illegal code

This has several disadvantages. The code is less readable, not clearly showing our intent and argument restriction. Also, in case of bad use of the function, the error message is delayed until no specialization is found. For instance you can write:

func(Val{1})

and you will get the following error message:

julia> ERROR: MethodError: no method matching specialized(::Int64, ::Type{Val{1}})
Closest candidates are:
  specialized(::Int64, !Matched::Type{Val{A::Tags = 0}}) at none:1
  specialized(::Int64, !Matched::Type{Val{B::Tags = 1}}) at none:1
  specialized(::Int64, !Matched::Type{Val{C::Tags = 2}}) at none:1
Stacktrace:
 [1] func(::Type{Val{1}}) at ./none:1

The problem is that you have to search through the Stacktrace to find the origin of the error. In this simple example the answer is immediate but for deeper stracktrace maybe this can be the source of painful debugging.

In comparison, a compilation attempt of this C++ code:

func(Tags_v<1>)

prints an error message pointing directly at the call site:

test.cpp:22:21: error: invalid conversion from ‘int’ to ‘Tags’ [-fpermissive]
   func(Tags_v<1>);
               ^~~~~~~~~

3 `Type` based solution

In C++ a possible implementation is:

#include <iostream>
#include <type_traits>

struct Tags {};
struct A : Tags {};
struct B : Tags {};
struct C : Tags {};

auto specialized(int x, A) { return 1 * x; }
auto specialized(int x, B) { return 2 * x; }
auto specialized(int x, C) { return 3 * x; }

template <typename TAG,typename ENABLED = std::enable_if_t<std::is_base_of<Tags,TAG>::value>>
auto func(TAG tag)
{
  return specialized(10, tag); // no run-time penalty
}

int main() 
{ 
  std::cout << func(B()); 
}

A Julia equivalent can be:

abstract type Tags end

struct A <: Tags end
struct B <: Tags end
struct C <: Tags end

specialized(x::Int,::Type{A}) = 1*x
specialized(x::Int,::Type{B}) = 2*x
specialized(x::Int,::Type{C}) = 3*x

func(::Type{T}) where {T<:Tags} = specialized(10,T)

func(B)

A look at the generated assembly code also confirms that everything is inlined:

@code_native func(B)

	.text
Filename: none
	pushq	%rbp
	movq	%rsp, %rbp
Source line: 1
	movl	$20, %eax
	popq	%rbp
	retq
	nopl	(%rax,%rax)

Comparison with C++:

The Julia code is still shorter. The two syntax

func(B)

versus

func(B()); // (or func<A>() if we choose another implementation)

is in favor of Julia (as Julia directly supports DataType, no instantiation of B is required).

With the Type approach Julia does not suffer anymore from the “argument filtering” problem we had with the @enum approach.

For instance:

func(1)

prints the following error message:

  : ERROR: MethodError: no method matching func(::Int64)
  : Closest candidates are:
  :   func(!Matched::Type{T<:Tags}) where T<:Tags at none:1

which points directly to the call site.

4 Default value and keyword argument

4.1 Default value

In C++ or Julia we can modify the func function to support a default value:

In C++:

enum approach

template <Tags TAG=Tags::A>
auto func(std::integral_constant<Tags, TAG> tag=Tags_v<A>)
{
  return specialized(10, tag); 
}

std::cout << func()           // prints 10
std::cout << func(Tags_v<B>)  // prints 20

Type approach

template <typename TAG=A,typename ENABLED = std::enable_if_t<std::is_base_of_v<Tags,TAG>>>
auto func(TAG tag=A())
{
  return specialized(10, tag); // no run-time penalty
}

std::cout << func()     // prints 10
std::cout << func(B())  // prints 20

In Julia:

@enum approach

func(::Type{Val{e}}=Val{A}) where {e} = specialized(10,Val{e})

func()       # prints 10
func(Val{B}) # prints 20

Type approach

func(::Type{T}=A) where {T<:Tags} = specialized(10,T)
func()  # prints 10
func(B) # prints 20

We can also check with @code_native that everything is inlined as before.

Comparison with C++:

Here there is a clear advantage in favor of Julia. In the C++ code you have to set the default value at two different places.

4.2 Keyword argument

There is no direct solution in C++. In Julia you can write:

@enum approach

func_kwa(;tag::Type{Val{e}}=Val{A}) where {e} = specialized(10,Val{e})

Type approach

func_kwa(;tag::Type{T}=A) where {T<:Tags} = specialized(10,T)

However there is a bad surprise if you look at the generated code (my Julia version is v0.6):

@code_native func_kwa(tag=B)

The code is so long that I put it at the end of this post, a lot of work is done at run-time! In our context this is not an acceptable solution due to the induced performance penalty.

5 Conclusions

We have compared two solutions to select a specialized subroutine given a tag passed as function argument. One solution is based on @enum, the other one is based on Type.

For the considered basic example we have verified that there is not run-time penalty and that like in C++, Julia can inline the function.

IMHO the better solution seems to be one based on Type. The reason are twofold:

better code readability,
AFAIK with the @enum/Type{Val{.}} approach, one can not filter argument directly at the call site.

We also have checked that:

we can use default argument without run-time penalty,
we can not use keyword argument because it seems that a lot of work is done at run-time (at least with my Julia version).

6 Annex: asm code for `func_kwa(tag=B)`

Julia version:

 versioninfo()

  Julia Version 0.6.2
  Commit d386e40c17 (2017-12-13 18:08 UTC)
  Platform Info:
    OS: Linux (x86_64-pc-linux-gnu)
    CPU: Intel(R) Xeon(R) CPU E5-2603 v3 @ 1.60GHz
    WORD_SIZE: 64
    BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
    LAPACK: libopenblas64_
    LIBM: libopenlibm
    LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

Note: same conclusion with the @enum approach, func_kwa(tag=Val{B}).

@code_native func_kwa(tag=B)

	   .text
   Filename: <missing>
	   pushq	%rbp
	   movq	%rsp, %rbp
	   pushq	%r15
	   pushq	%r14
	   pushq	%r13
	   pushq	%r12
	   pushq	%rbx
	   subq	$120, %rsp
	   movq	%rdi, %r14
	   movabsq	$140389732029024, %r13  # imm = 0x7FAF081B9260
	   movq	%fs:0, %r15
	   addq	$-10888, %r15           # imm = 0xD578
	   leaq	-64(%rbp), %r12
	   xorps	%xmm0, %xmm0
	   movups	%xmm0, -64(%rbp)
	   movq	$0, -48(%rbp)
	   movups	%xmm0, -96(%rbp)
	   movups	%xmm0, -112(%rbp)
	   movups	%xmm0, -128(%rbp)
	   movups	%xmm0, -144(%rbp)
	   movq	$0, -80(%rbp)
	   movq	$26, -160(%rbp)
	   movq	(%r15), %rax
	   movq	%rax, -152(%rbp)
	   leaq	-160(%rbp), %rax
	   movq	%rax, (%r15)
	   movq	$0, -72(%rbp)
   Source line: 0
	   movq	$0, -120(%rbp)
	   movq	8(%r14), %rcx
	   sarq	%rcx
	   testq	%rcx, %rcx
	   jle	L271
	   movq	(%r14), %rdx
	   movq	24(%r14), %rsi
	   movl	$1, %eax
	   leaq	-24387768(%r13), %rdi
	   nopw	(%rax,%rax)
   L176:
	   leaq	-1(%rax), %rbx
	   cmpq	%rsi, %rbx
	   jae	L431
	   movq	-8(%rdx,%rax,8), %rbx
	   testq	%rbx, %rbx
	   je	L364
	   movq	%rbx, -144(%rbp)
	   movq	%rbx, -136(%rbp)
	   cmpq	%rdi, %rbx
	   jne	L465
	   cmpq	%rsi, %rax
	   jae	L379
	   movq	(%rdx,%rax,8), %rbx
	   testq	%rbx, %rbx
	   je	L416
	   movq	%rbx, -128(%rbp)
	   movq	%rbx, -120(%rbp)
	   addq	$2, %rax
	   decq	%rcx
	   jne	L176
	   movq	%rbx, -80(%rbp)
	   jmp	L286
   L271:
	   leaq	43425776(%r13), %rbx
	   movq	%rbx, -120(%rbp)
	   movq	%rbx, -80(%rbp)
   L286:
	   leaq	-22171552(%r13), %rax
	   movq	%rax, -64(%rbp)
	   movq	%rbx, -56(%rbp)
	   addq	$-22171848, %r13        # imm = 0xFEADAF38
	   movq	%r13, -48(%rbp)
	   movabsq	$jl_apply_generic, %rax
	   movl	$3, %esi
	   movq	%r12, %rdi
	   callq	*%rax
	   movq	%rax, -72(%rbp)
	   movq	(%rax), %rax
	   movq	-152(%rbp), %rcx
	   movq	%rcx, (%r15)
	   leaq	-40(%rbp), %rsp
	   popq	%rbx
	   popq	%r12
	   popq	%r13
	   popq	%r14
	   popq	%r15
	   popq	%rbp
	   retq
   L364:
	   movabsq	$jl_throw, %rax
	   movq	%r13, %rdi
	   callq	*%rax
   L379:
	   movq	%rsp, %rcx
	   leaq	-16(%rcx), %rsi
	   movq	%rsi, %rsp
	   incq	%rax
	   movq	%rax, -16(%rcx)
	   movabsq	$jl_bounds_error_ints, %rax
	   movl	$1, %edx
	   movq	%r14, %rdi
	   callq	*%rax
   L416:
	   movabsq	$jl_throw, %rax
	   movq	%r13, %rdi
	   callq	*%rax
   L431:
	   movq	%rsp, %rcx
	   leaq	-16(%rcx), %rsi
	   movq	%rsi, %rsp
	   movq	%rax, -16(%rcx)
	   movabsq	$jl_bounds_error_ints, %rax
	   movl	$1, %edx
	   movq	%r14, %rdi
	   callq	*%rax
   L465:
	   movabsq	$jl_gc_pool_alloc, %rax
	   movl	$1456, %esi             # imm = 0x5B0
	   movl	$32, %edx
	   movq	%r15, %rdi
	   callq	*%rax
	   movq	%rax, %rbx
	   leaq	-21540240(%r13), %rax
	   movq	%rax, -8(%rbx)
	   movq	%rbx, -112(%rbp)
	   xorps	%xmm0, %xmm0
	   movups	%xmm0, (%rbx)
	   movq	-648819312(%r13), %rax
	   movq	56(%rax), %rax
	   testq	%rax, %rax
	   jne	L545
	   movabsq	$jl_throw, %rax
	   movq	%r13, %rdi
	   callq	*%rax
   L545:
	   movq	%rax, -104(%rbp)
	   movq	%rax, -64(%rbp)
	   leaq	-25545896(%r13), %rax
	   movq	%rax, -56(%rbp)
	   movabsq	$jl_f_getfield, %rax
	   xorl	%edi, %edi
	   movl	$2, %edx
	   movq	%r12, %rsi
	   callq	*%rax
	   movq	%rax, -96(%rbp)
	   movq	%rax, (%rbx)
	   testq	%rax, %rax
	   je	L632
	   movq	-8(%rbx), %rcx
	   andl	$3, %ecx
	   cmpq	$3, %rcx
	   jne	L632
	   testb	$1, -8(%rax)
	   jne	L632
	   movabsq	$jl_gc_queue_root, %rax
	   movq	%rbx, %rdi
	   callq	*%rax
   L632:
	   movl	$1432, %esi             # imm = 0x598
	   movl	$16, %edx
	   movq	%r15, %rdi
	   movabsq	$jl_gc_pool_alloc, %rax
	   callq	*%rax
	   addq	$-647092848, %r13       # imm = 0xD96E2590
	   movq	%r13, -8(%rax)
	   movq	%rax, -88(%rbp)
	   movq	%r14, (%rax)
	   movq	%rax, 8(%rbx)
	   testq	%rax, %rax
	   je	L712
	   movq	-8(%rbx), %rax
	   andl	$3, %eax
	   cmpq	$3, %rax
	   jne	L712
	   movabsq	$jl_gc_queue_root, %rax
	   movq	%rbx, %rdi
	   callq	*%rax
   L712:
	   movq	$-1, 16(%rbx)
	   movabsq	$jl_throw, %rax
	   movq	%rbx, %rdi
	   callq	*%rax
	   nop

Direct convolution

By: Picaud Vincent

Re-posted from: https://pixorblog.wordpress.com/2016/07/17/direct-convolution/

For small kernels, direct convolution beats FFT based one. I present here a basic implementation. This implementation allows to compute

$\gamma[k]=\sum\limits_{i\in\Omega^\alpha}\alpha[i]\beta[k+\lambda i],\text{ with }\lambda\in\mathbb{Z}^*\text{\ \ \ \ \ \ \ \ \ \ \ \ (1)}$

From time to time we will use the notation $\gamma=\alpha\bigodot\limits_\lambda\beta$ .

An arbitrary stride $\lambda$ has been introduced to define:

$\lambda=-1$ convolution
$\lambda=+1$ cross-correlation
$\lambda=\pm 2^n$ the stationary wavelet transform (the so called “à trous” algorithm)

Also note that with proper boundary extension (periodic and zero padding essentially), changing the sign of $\lambda$ gives the adjoint operator:

$\langle \delta, \alpha\bigodot\limits_\lambda\beta \rangle = \langle \alpha\bigodot\limits_{-\lambda}\delta,\beta \rangle$

Disclaimer

Maybe the following is overwhelmingly detailed for a simple task like Eq. (1), but I have found some interests in writing this once for all. Maybe it can be useful for someone else.

Some notations

We note $\Omega$ our vector domain (or support), for instance

$\Omega^\alpha =\llbracket i_{\min} ,i_{\max} \rrbracket$

means that $\alpha[i]$ is defined for

$i\in \llbracket i_{\min}, i_{\max} \rrbracket$

To get interval lower/upper bounds we use the notation

$i_{\min}=l(\Omega^\alpha)\text{ and }i_{\max}=u(\Omega^\alpha)$

We denote by $\lambda\Omega$ the scaled domain $\Omega$ defined by:

$\lambda\Omega=\llbracket \lambda^+\,i_{min}+\lambda^-\,i_{max}, \lambda^+\,i_{max}+\lambda^-\,i_{min} \rrbracket$

where $\lambda^+=\max{(0,\lambda)}$ and $\lambda^-=\min{(0,\lambda)}$

Finally we use $A\setminus B$ the relative complement of $B$ with respect to the set $A$ defined by

$A\setminus B = \{ i\ /\ (i\in A) \wedge (i\notin B) \}$

This set is not necessary connex, however like we are working in $\mathbb{Z}$ , it is sufficient to introduce the left and right parts (that can be empty)

$(A\setminus B)_{\text{Left}}=\llbracket l(A), \min{(u(A),l(B)-1)} \rrbracket$

$(A\setminus B)_{\text{Right}}=\llbracket \max{(l(A),u(B)+1)}, u(A) \rrbracket$

Goal

Given two vectors $\alpha$ , $\beta$ defined on $\Omega^\alpha$ , $\Omega^\beta$ we want to define and implement an algorithm that computes $\gamma[k]$ for $k\in\Omega^\gamma$ .

First step, no boundary extension

We need to define the $\Omega^\gamma_1$ the domain that does not violate $\beta$ domain of definition. This can be expressed as

$\Omega^\gamma_1=\{k\in\mathbb{Z}\ /\ \forall i \in \Omega^\alpha \Rightarrow k+\lambda i \in \Omega^\beta \}$

Let’s write the details,

$(\forall i \in \Omega^\alpha \Rightarrow k+\lambda i \in \Omega^\beta)\Leftrightarrow (\forall i \in \Omega^\alpha \Rightarrow l(\Omega^\beta)-\lambda i \le k \le u(\Omega^\beta)-\lambda i)$

$\Leftrightarrow \max\limits_{i\in \Omega^\alpha} l(\Omega^\beta)-\lambda i \le k \le \min\limits_{i\in \Omega^\alpha} u(\Omega^\beta)-\lambda i$

$\Leftrightarrow l(\Omega^\beta)-l(\lambda \Omega^\alpha) \le k \le u(\Omega^\beta)-u(\lambda \Omega^\alpha)$

hence we have

$\boxed{ \Omega^\gamma_1=\llbracket l(\Omega^\beta)-l(\lambda \Omega^\alpha) , u(\Omega^\beta)-u(\lambda \Omega^\alpha) \rrbracket }$

Thus the computation of $\gamma[k],\ k\in\Omega^\gamma$ (Eq. 1) is splitted into two parts:

one part $\Omega^\gamma \cap \Omega^\gamma_1$ free of boundary effect,
one part $\Omega^\gamma \setminus \Omega^\gamma_1$ that requires boundary extension $\tilde{\beta}=\Phi(\beta,k)$

The algorithm takes the following form:

$latex-test.png$

Second step, boundary extensions

Usually we define some classical boundary extensions. These extensions are computed from $\beta[.]$ and are sometimes entailed by a validity condition. For a better clarity I give explicit lower/upper bounds:

$\Omega^\beta = \llbracket j_{\min} , j_{\max} \rrbracket \neq \emptyset$

Left boundary $(j<j_{min})$	$\tilde{\beta}_j = \Phi_{left}(\beta,j)$	validity condition
Mirror	$\tilde{\beta}_j = \beta[2\,j_{min}-j]$	$2\,j_{min}-j_{max} \le j$
Periodic (or cyclic)	$\tilde{\beta}_j = \beta[j_{max}-j_{min}+j+1]$	$2\,j_{min}-j_{max}-1 \le j$
Constant	$\tilde{\beta}_j = \beta[j_{min}]$	none
Zero padding	$\tilde{\beta}_j = 0$	none

Right boundary $(j>j_{max})$	$\tilde{\beta}_j = \Phi_{right}(\beta,j)$	validity condition
Mirror	$\tilde{\beta}_j = \beta[2\,j_{max}-j]$	$j\le 2\,j_{max}-j_{min}$
Periodic (or cyclic)	$\tilde{\beta}_j = \beta[-j_{max}+j_{min}+j-1]$	$j\le 2\,j_{max}-j_{min}+1$
Constant	$\tilde{\beta}_j = \beta[j_{max}]$	none
Zero padding	$\tilde{\beta}_j = 0$	none

As we want something general we want to get rid of these validity conditions.

Periodic case

Starting from a vector $\beta$ defined on $\llbracket L=0, U \rrbracket$ we want to define a periodic function $\tilde{\beta}$ of period $T=U+1$ . This function must fulfills the $\tilde{\beta}[j+T]=\tilde{\beta}[j]$ relation.

We can do that by considering $\tilde{\beta}=\beta \circ \phi^P_U(j)$ where

$\phi^P_U(j)=\bmod_F(j,U+1)$

and $\bmod_F$ is the modulus function associated to a floored division.

For a vector defined on an arbitrary domain $\llbracket j_{\min}, j_{\max} \rrbracket$ , we first translate the indices

$\tau_{j_{\min}}(j)=j-j_{\min}$

and then translate them back using $\tau^{(-1)}_{j_{\min}}=\tau_{-j_{\min}}$

Putting all together, we build a periodized vector

$\boxed{\tilde{\beta} = \beta \circ \phi^P_{j_{\min},j_{\max}}}$

where

$\phi^P_{j_{\min},j_{\max}} = \tau^{(-1)}_{j_{\min}} \circ \phi^P_{j_{\max}- j_{\min}} \circ \tau_{j_{\min}}$

$\boxed{\phi^P_{j_{\min},j_{\max}} = j_{\min} + \bmod_F(j-j_{\min},j_{\max}- j_{\min}+1)}$

Mirror Symmetry case

Starting from a vector $\beta$ defined on $\llbracket L=0, U \rrbracket$ we can extend it by mirror symmetry on $\llbracket U+1, 2U \rrbracket$ using $\tilde{\beta}=\beta\circ \phi^M_U$ with

$\phi^M_U(j)=U-|U-j|$

The resulting vector $\tilde{\beta}=\beta\circ \phi^M_U$ fulfills the $\tilde{\beta}[U-j]=\tilde{\beta}[U+j]$ relation for $j\in \llbracket 0, U \rrbracket$ .

To get a “global” definition we then periodize it on $\llbracket 0, 2U-1 \rrbracket$ using $\phi^P_{2U-1}$ (attention $2U-1$ and not $2U$ , otherwise the component $0$ is duplicated!).

For an arbitrary domain $\llbracket j_{\min}, j_{\max} \rrbracket$ we use index translation as for the periodic case. Putting everything together we get:

$\boxed{\tilde{\beta} = \beta \circ \phi^M_{j_{\min},j_{\max}}}$

where

$\phi^M_{j_{\min},j_{\max}} = \tau^{(-1)}_{j_{\min}} \circ \phi^M_{j_{\max}- j_{\min}} \circ \phi^P_{2(j_{\max}- j_{\min})-1} \circ \tau_{j_{\min}}$

$\boxed{ \phi^M_{j_{\min},j_{\max}} =j_{\max}-|j_{\max}-j_{\min}-\bmod_F(j-j_{\min},2(j_{\max}-j_{\min}))| }$

Boundary extensions

To use the algorithm with boundary extensions, you only have to define:

$\tilde{\beta}=\Phi(\beta,k+\lambda i)=\beta[\phi^X[k+\lambda i]]$

where $X$ is the boundary extension you have chosen (periodic, constant). You do not have to take care of any validity condition, these formula are general.

Implementation

This is a straightforward implementation following as close as possible the presented formula. We did not try to optimize it, this would have obscured the presentation. Some ideas: reverse $\alpha$ for $\lambda<0$ (access memory in the right order), use simd, or C++ meta-programming with loop unrolling for fixed $\alpha$ size, specialize regarding to Vector/StridedVector or $\lambda=\pm 1$

Preamble

Index translation / domain definition

There is however one last thing we have to explain. In languages like Julia, C we are manipulating arrays having a common starting index: $1$ in Julia, Fortran or $0$ in C, C++

For this reason we do not manipulate $\alpha$ on $\Omega^\alpha$ but an another translated array $\tilde{\alpha}$ defined on $\llbracket 1, N^\alpha \rrbracket$ (Julia) or $\llbracket 0, N^\alpha-1 \rrbracket$ (C++).

To cover all cases, I assume that the starting index is denoted by $\tilde{i}_0$ .

The array $\tilde{\alpha}$ is defined by:

$\alpha[i] = \tilde{\alpha}[\tilde{i}] = \tilde{\alpha}[i-l(\Omega^\alpha)+\tilde{i}_0]$

Hence we must modify the initiale Eq. (1) to use $\tilde{\alpha}$ instead of $\alpha$

$\gamma[k]=\sum\limits_{i\in\Omega^\alpha}\alpha[i]\beta[k+\lambda i] = \sum\limits_{i\in\Omega^\alpha}\tilde{\alpha}[i-l(\Omega^\alpha)+\tilde{i}_0]\beta[k+\lambda i]$

With $\tilde{i}=i-l(\Omega^\alpha)+\tilde{i}_0$ we have

$i\in\Omega^\alpha \Leftrightarrow \tilde{i}\in\llbracket \tilde{i}_0,u(\Omega^\alpha)-l(\Omega^\alpha)+\tilde{i}_0 \rrbracket$

and

$k+\lambda i = k+ \lambda \tilde{i} + \underbrace{\lambda (l(\Omega^\alpha) - \tilde{i}_0)}_{\beta\_\text{offset}}$

Thus, Eq (1) becomes:

$\boxed{ \gamma[k]=\sum\limits_{\tilde{i}=\tilde{i}_0}^{u(\Omega^\alpha)-l(\Omega^\alpha)+\tilde{i}_0}\tilde{\alpha}[\tilde{i}]\beta[k+ \lambda \tilde{i} + \lambda (l(\Omega^\alpha) - \tilde{i}_0)]}$

The $2$ other arrays are less problematic:

For $\beta$ array, which is our input array, we implicitly use $\Omega^\beta = \llbracket \tilde{i}_0, \tilde{i}_0 + \text{length}(\beta) - 1 \rrbracket$ . This does not reduce the generality of the subroutine.
For $\gamma$ which is the output array, as for $\beta$ we assume it is defined on $\llbracket \tilde{i}_0, \tilde{i}_0 + \text{length}(\gamma) - 1 \rrbracket$ , but we provide $\Omega^\gamma\subset \llbracket \tilde{i}_0, \tilde{i}_0 + \text{length}(\gamma) - 1 \rrbracket$ to define the components we want to compute. The other components, $\llbracket \tilde{i}_0, \tilde{i}_0 + \text{length}(\gamma) - 1 \rrbracket \setminus \Omega^\gamma$ , will remain unmodified by the subroutine.

Definition of $\alpha\_\text{offset}$

As we have seen before, the convolution subroutine will have $\tilde{\alpha}$ as argument, but we also need $\Omega^\alpha$ . For the driver subroutine we do not directly provide this interval because its length is redundant with $\tilde{\alpha}$ length. Instead we provide an $\alpha\_\text{offset}$ offset. $\Omega^\alpha$ is deduced from:

$\Omega^\alpha = \llbracket -\alpha\_\text{offset}, -\alpha\_\text{offset} + \text{length}(\tilde{\alpha}) -1 \rrbracket$

Note: this definition does not depend on $\tilde{i}_0$ .

With $\alpha\_\text{offset}=0$ you are in the “usual situation”. If you have a window size of $2n+1$ , taking $\alpha\_\text{offset}=n$ returns the middle of the window. Here, in the Fig. below, the graphical representation of an arbitrary case: a filter if size $4$ , with $\alpha\_\text{offset}=2$ and $\lambda=3$ .

Julia

Auxiliary subroutines

We start by defining the basic operations on sets:

function scale(λ::Int64,Ω::UnitRange)
    ifelse(λ>0,
           UnitRange(λ*start(Ω),λ*last(Ω)),
           UnitRange(λ*last(Ω),λ*start(Ω)))
end

function compute_Ωγ1(Ωα::UnitRange,
                     λ::Int64,
                     Ωβ::UnitRange)
    
    λΩα = scale(λ,Ωα)

    UnitRange(start(Ωβ)-start(λΩα),
              last(Ωβ)-last(λΩα))
end

# Left & Right relative complements A\B
#
function relelativeComplement_left(A::UnitRange,
                                   B::UnitRange)
    UnitRange(start(A),
              min(last(A),start(B)-1))
end

function relelativeComplement_right(A::UnitRange,
                                    B::UnitRange)
    UnitRange(max(start(A),last(B)+1),
              last(A))
end

Boundary extensions

We then define the boundary extensions. Nothing special there, we only had to check that the Julia mod(x,y) function is the floored division version (by opposition to the rem(x,y) function which is the rounded toward zero division version).

const tilde_i0 = Int64(1)

function boundaryExtension_zeroPadding{T}(β::StridedVector{T},
                                          k::Int64)
    kmin = tilde_i0
    kmax = length(β) + kmin - 1
    
    if (k>=kmin)&&(k<=kmax)
        β[k]
    else
        T(0)
    end
end

function boundaryExtension_constant{T}(β::StridedVector{T},
                                       k::Int64)
    kmin = tilde_i0
    kmax = length(β) + kmin - 1

    if k<kmin
        β[kmin]
    elseif k<=kmax
        β[k]
    else
        β[kmax]
    end
end

function boundaryExtension_periodic{T}(β::StridedVector{T},
                                       k::Int64)
    kmin = tilde_i0
    kmax = length(β) + kmin - 1

    β[kmin+mod(k-kmin,1+kmax-kmin)]
end

function boundaryExtension_mirror{T}(β::StridedVector{T},
                                     k::Int64)
    kmin = tilde_i0
    kmax = length(β) + kmin - 1

    β[kmax-abs(kmax-kmin-mod(k-kmin,2*(kmax-kmin)))]
end

# For the user interface
#
boundaryExtension = 
    Dict(:ZeroPadding=>boundaryExtension_zeroPadding,
         :Constant=>boundaryExtension_constant,
         :Periodic=>boundaryExtension_periodic,
         :Mirror=>boundaryExtension_mirror)

Main subroutine

Finally we define the main subroutine. Its arguments have been defined in the preamble part. I just added one @simd & @inbounds because this has a significant impact concerning perfomance (see end of this post).

function direct_conv!{T}(tilde_α::StridedVector{T},
                         Ωα::UnitRange,
                         λ::Int64,
                         β::StridedVector{T},
                         γ::StridedVector{T},
                         Ωγ::UnitRange,
                         LeftBoundary::Symbol,
                         RightBoundary::Symbol)
    # Sanity check
    @assert λ!=0
    @assert length(tilde_α)==length(Ωα)
    @assert (start(Ωγ)>=1)&&(last(Ωγ)<=length(γ))

    # Initialization
    Ωβ = UnitRange(1,length(β))
    tilde_Ωα = 1:length(Ωα)
    
    for k in Ωγ
        γ[k]=0 
    end

    rΩγ1=intersect(Ωγ,compute_Ωγ1(Ωα,λ,Ωβ))
    
    # rΩγ1 part: no boundary effect
    #
    β_offset = λ*(start(Ωα)-tilde_i0)
    @simd for k in rΩγ1
        for i in tilde_Ωα
            @inbounds γ[k]+=tilde_α[i]*β[k+λ*i+β_offset]
        end
    end

    # Left part
    #
    rΩγ1_left = relelativeComplement_left(Ωγ,rΩγ1)
    Φ_left = boundaryExtension[LeftBoundary]
    
    for k in rΩγ1_left
        for i in tilde_Ωα
            γ[k]+=tilde_α[i]*Φ_left(β,k+λ*i+β_offset)
        end
    end

    # Right part
    #
    rΩγ1_right = relelativeComplement_right(Ωγ,rΩγ1)
    Φ_right = boundaryExtension[RightBoundary]
    
    for k in rΩγ1_right
        for i in tilde_Ωα
            γ[k]+=tilde_α[i]*Φ_right(β,k+λ*i+β_offset)
        end
    end
end

# Some UI functions, γ inplace modification 
#
function direct_conv!{T}(tilde_α::StridedVector{T},
                         α_offset::Int64,
                         λ::Int64,

                         β::StridedVector{T},

                         γ::StridedVector{T},
                         Ωγ::UnitRange,
                         
                         LeftBoundary::Symbol,
                         RightBoundary::Symbol)

    Ωα = UnitRange(-α_offset,
                   length(tilde_α)-α_offset-1)
    
    direct_conv!(tilde_α,
                 Ωα,
                 λ,
                 
                 β,

                 γ,
                 Ωγ,

                 LeftBoundary,
                 RightBoundary)
end

# Some UI functions, allocates γ 
#
function direct_conv{T}(tilde_α::StridedVector{T},
                        α_offset::Int64,
                        λ::Int64,

                        β::StridedVector{T},

                        LeftBoundary::Symbol,
                        RightBoundary::Symbol)

    γ = Array{T,1}(length(β))
    
    direct_conv!(tilde_α,
                 α_offset,
                 λ,

                 β,

                 γ,
                 UnitRange(1,length(γ)),

                 LeftBoundary,
                 RightBoundary)

    γ
end

In C/C++

As this post is already long I will not provide a complete code here. The only trap is to use the right mod function.

C/C++ modulus operator % is not standardized. Only the D%d=D-d*(D/d) relation is invariant allowing to define the Euclidean division. On the other side a lot of CPU x86 idiv, truncate toward zero, as a consequence C/C++ generally uses this direction.

To be sure, we have to explicitly use our F-mod function:

// Floored mod
int modF(int D, int d)
{
    int r = std::fmod(D,d);
    if((r > 0 && d < 0) || (r < 0 && d > 0)) r = r + d;
    return r;
}

You can read:

Usages examples

Basic usages

Beware that due to the asymmetric role of $\alpha$ and $\beta$ the proposed approach does preserve all the mathematical properties of the $\alpha\bigodot\limits_\lambda\beta$ operator.

Commutativity:

$\alpha\bigodot\limits_{\lambda=-1}\beta=\beta\bigodot\limits_{\lambda=-1}\alpha$

only for ZeroPadding

Adjoint operator:

$\forall \lambda\in\mathbb{Z}^*,\ \langle \alpha\bigodot\limits_{\lambda}v ,w \rangle_E = \langle v , \alpha\bigodot\limits_{-\lambda} w \rangle_F$

only for ZeroPadding and Periodic

I have assumed $\mathbb{R}$ arrays (not $\mathbb{C}$ ones): some conjugation are missing
Not considered here, but extension to n-dimensional & separable filters is immediate

push!(LOAD_PATH,"./")
using DirectConv

α=rand(4);
β=rand(10);

# Check adjoint operator
# -> restricted to ZeroPadding & Periodic
#    (asymmetric role of α and β)
#    
vβ=rand(length(β))
d1=dot(direct_conv(α,2,-3,vβ,:ZeroPadding,:ZeroPadding),β)
d2=dot(direct_conv(α,2,+3,β,:ZeroPadding,:ZeroPadding),vβ)

@assert abs(d1-d2)<sqrt(eps())

d1=dot(direct_conv(α,-1,-3,vβ,:Periodic,:Periodic),β)
d2=dot(direct_conv(α,-1,+3,β,:Periodic,:Periodic),vβ)

@assert abs(d1-d2)<sqrt(eps())

# Check commutativity 
# -> λ = -1 (convolution) and
#    restricted to ZeroPadding
#    (asymmetric role of α and β)
v1=zeros(20)
v2=zeros(20)
direct_conv!(α,0,-1,
             β,v1,UnitRange(1,20),:ZeroPadding,:ZeroPadding)
direct_conv!(β,0,-1,
             α,v2,UnitRange(1,20),:ZeroPadding,:ZeroPadding)

@assert (norm(v1-v2)<sqrt(eps()))

# Check Interval splitting
# (should work for any boundary extension type)
#
γ=direct_conv(α,3,2,β,:Mirror,:Periodic) # global computation
Γ=zeros(length(γ))
Ω1=UnitRange(1:3)
Ω2=UnitRange(4:length(γ))
direct_conv!(α,3,2,β,Γ,Ω1,:Mirror,:Periodic) # compute on Ω1
direct_conv!(α,3,2,β,Γ,Ω2,:Mirror,:Periodic) # compute on Ω2

@assert (norm(γ-Γ)<sqrt(eps()))

Performance?

In a previous post I gave a short derivation of the Savitzky-Golay filters. I used a FFT based convolution to apply the filters. It is interesting to compare the performance of the presented direct approach vs the FFT one.

push!(LOAD_PATH,"./")
using DirectConv

function apply_filter{T}(filter::StridedVector{T},signal::StridedVector{T})

    @assert isodd(length(filter))

    halfWindow = round(Int,(length(filter)-1)/2)
    
    padded_signal = 
        [signal[1]*ones(halfWindow);
         signal;
         signal[end]*ones(halfWindow)]

    filter_cross_signal = conv(filter[end:-1:1],
                               padded_signal)

    filter_cross_signal[2*halfWindow+1:end-2*halfWindow]
end

# Now we can create a (very) rough benchmark
M=Array(Float64,0,3)
β=rand(1000000);
for halfWidth in 1:2:40
    α=rand(2*halfWidth+1);

    fft_t0 = time()
    fft_v = apply_filter(α,β)
    fft_t1 = time()

    direct_t0 = time()
    direct_v = direct_conv(α,halfWidth,1,β,
                           :Constant,:Constant)
    direct_t1 = time()

    @assert (norm(fft_v -direct_v)<sqrt(eps()))

    M=vcat(M,
           Float64[length(α)
                   (fft_t1-fft_t0)*1e3
                   (direct_t1-direct_t0)*1e3])
end
M

We see that for small filters direct method can easily be 10 time faster than the FFT approach!

Conclusion: for small filters, use a direct approach!

Discussion

Optimization/performance

If I have time I will try to benchmark two basic implementations, a Julia one vs a C/C++ one. I’m a beginner in Julia language, with C++, I’m more at home.

I would be curious to see the difference between a basic implementation and an optimized one in Julia. Just to see how optimization can obfuscate (or not) the initial code and the performance gain. In C++ you generally have a lot of boiler-plate code (meta-programming).

Applications

The basic Eq. (1) is common tool that can be used for:

deconvolution procedures,
decimated and undecimated wavelet transforms,

For wavelet transform especially the undecimated one, AFAIK Eq. (1) is really the good choice. I will certainly write some posts on these stuff.

Some extra reading:

The FFT way: Algorithms for Efficient Computation of Convolution, K. Pavel
The Winograd’s minimal filtering algorithms way: Fast Algorithms for Convolutional Neural Networks, A. Lavin, S. Gray
The OpenCL/GUPU way: Case study: High performance convolution using OpenCL __local memory

Code

The code is on github.

Complement: more domains

The $\Omega^\gamma_2$ domain

We have introduced $\Omega^\gamma_1$ the domain that does not violate $\beta$ domain of definition (given $\Omega^\alpha$ and $\Omega^\beta$ ).

To be exhaustive we can introduce $\Omega^\gamma_2$ the domain that use at least one $(i,k+\lambda i)\in \Omega^\alpha \times \Omega^\beta$ .

This domain is:

$\Omega^\gamma_2=\{ k\in\mathbb{Z}\ /\ \exists i \in \Omega^\alpha \Rightarrow k+\lambda i \in \Omega^\beta \}$

following arguments similar to those used for $\Omega^\gamma_1$ we get:

$\boxed{ \Omega^\gamma_2=\llbracket l(\Omega^\beta)-u(\lambda \Omega^\alpha) , u(\Omega^\beta)-l(\lambda \Omega^\alpha) \rrbracket }$

The $\Omega^\beta_{2'}$ domain

We can also ask for the “dual” question: given $\Omega^\alpha$ and $\Omega^\gamma$ what is the domain of $\beta$ , $\Omega^\beta_{2'}$ , involved in the computation of $\gamma$

By definition, this domain must fulfill the following relation:

$\Omega^\gamma_2(\Omega^\beta_{2'})=\Omega^\gamma$

hence, using the previous result

$\llbracket l(\Omega^\beta_{2'})-u(\lambda \Omega^\alpha) , u(\Omega^\beta_{2'})-l(\lambda \Omega^\alpha) \rrbracket = \llbracket l(\Omega^\gamma),u(\Omega^\gamma) \rrbracket$

which gives:

$\boxed{ \Omega^\beta_{2'} = \llbracket l(\Omega^\gamma)+u(\lambda \Omega^\alpha),u(\Omega^\gamma)+l(\lambda \Omega^\alpha) \rrbracket }$

Savitzky-Golay filters & Julia

By: Picaud Vincent

Re-posted from: https://pixorblog.wordpress.com/2016/07/13/savitzky-golay-filters-julia/

I always have found that presentations of the Savitzky-Golay filters were over tricky. I present here a simple derivation of these formula and a possible implementation in Julia.

Derivation of the formula

Given a polynomial of degree $d$ with values

$p(x_i)=\sum\limits_{j=0}^d c_j x_i^j$

defined for a window of size $2n+1$ , $\{x_i,\ i\in\llbracket -n,n \rrbracket \}$

We want to find the value of its k-order derivative $p^{(k)}(x_0)$ in the middle of the window assuming that the $c_j$ are founded solving a least-squares problem:

$\min\limits_{\mathbf{c}} \frac{1}{2} \| \mathbf{V} \mathbf{c} - \mathbf{y} \|_2^2$

where $\{y_i, i \in\llbracket -n,n \rrbracket \}$ is our signal values and $\mathbf{V}$ is the Vandermonde matrix:

$\mathbf{V}= \left( \begin{array}{c|c|c} \vdots & \vdots & \vdots \\ 1 & x_i^{(j-1)} & x_i^d \\ \vdots & \vdots & \vdots \end{array} \right)$

using the normal equation

$\mathbf{c}=(\mathbf{V}^t.\mathbf{V})^{-1}.\mathbf{V}^t.\mathbf{y}$

and a QR decomposition, $\mathbf{V}=\mathbf{Q}.\mathbf{R}$ we get

$\mathbf{c}=\mathbf{R}^{-1}.\mathbf{Q}^t.\mathbf{y}$

now we can express all the polynomial values $p(x_i)$ in a vector $\mathbf{p}=\mathbf{V}.\mathbf{c}$ . Lets rewrite this in matrix form:

$\underbrace{\left( \begin{array}{c} p(x_{-n}) \\ \vdots \\ p(x_{0}) \\ \vdots \\ p(x_{+n}) \end{array} \right)}\limits_{\mathbf{p}}=\underbrace{ \left( \begin{array}{c|c|c} \vdots & \vdots & \vdots \\ 1 & x_i^{(j-1)} & x_i^d \\ \vdots & \vdots & \vdots \end{array} \right)}\limits_{\mathbf{V}}.\underbrace{\left( \begin{array}{c} c_0 \\ \vdots \\ c_n \end{array} \right)}\limits_{\mathbf{c}}$

Now the “big trick” is to write the Taylor series and to remember that this formula is exact for polynomial functions:

$\forall i,\ P(x_i) = \sum\limits_{j=0}^d \frac{x_i^j}{j!} P^{(j)}(x_0)$

Lets rewrite this in matrix form:

$\underbrace{ \left( \begin{array}{c} p(x_{-n}) \\ \vdots \\ p(x_{0}) \\ \vdots \\ p(x_{n}) \\ \end{array} \right) }_{\mathbf{p}} = \underbrace{ \left( \begin{array}{c|c|c} \vdots & \vdots & \vdots \\ 1 & \frac{x_i^{(j-1)}}{(j-1)!} & \frac{x_i^{d}}{d!} \\ \vdots & \vdots & \vdots \end{array} \right) }_{\mathbf{T}} \underbrace{ \left( \begin{array}{c} P^{(0)}(x_0) \\ \vdots \\ P^{(k)}(x_0) \\ \vdots \\ P^{(d)}(x_0) \\ \end{array} \right) }_{\mathbf{p^\delta}}$

With a good eye we see that $\mathbf{V}=\mathbf{T}.\mathbf{D}$ where $\mathbf{D}$ is a diagonal matrix:

$\underbrace{ \left( \begin{array}{c|c|c} \vdots & \vdots & \vdots \\ 1 & x_i^{(j-1)} & x_i^d \\ \vdots & \vdots & \vdots \end{array} \right)}\limits_{\mathbf{V}} = \underbrace{ \left( \begin{array}{c|c|c} \vdots & \vdots & \vdots \\ 1 & \frac{x_i^{(j-1)}}{(j-1)!} & \frac{x_i^{d}}{d!} \\ \vdots & \vdots & \vdots \end{array} \right) }_{\mathbf{T}}.\underbrace{\left( \begin{array}{ccc} 1 & & \\ & (j-1)! & \\ & & d! \end{array} \right)}\limits_{\mathbf{D}}$

That’s all, we only have to group pieces:

$\mathbf{V}.\mathbf{c}=\mathbf{P}=\mathbf{T}.\mathbf{p^\delta}=\mathbf{V}.\mathbf{D}^{-1}.\mathbf{p^\delta}$

With the QR decomposition $\mathbf{V}=\mathbf{Q}.\mathbf{R}$ and $\mathbf{c}=\mathbf{R}^{-1}.\mathbf{Q}^t.\mathbf{y}$ we get:

$\mathbf{Q}.\mathbf{Q}^t.\mathbf{y}=\mathbf{Q}.\mathbf{R}.\mathbf{D}^{-1}.\mathbf{p^\delta}$

using the fact that $\mathbf{Q}^t.\mathbf{Q}=\mathbf{I}$ we get:

$\mathbf{Q}^t.\mathbf{y}=\mathbf{R}.\mathbf{D}^{-1}.\mathbf{p^\delta}$

hence we have:

$\boxed{ \mathbf{p^\delta} = \mathbf{D}.\mathbf{R}^{-1}.\mathbf{Q}^t.\mathbf{y} }$

which is our final formula.

Symbolic computation to check that it works

We can use mathematica to do a symbolic computation using $\mathbf{p^\delta} = \mathbf{D}.\mathbf{R}^{-1}.\mathbf{Q}^t.\mathbf{y}$ .

For a window width of $2n+1=7$ points and a polynomial of degree $d=2$ we get:

n = 3; d = 2;
V = Table[If[j != 0, i^j, 1], {i, -n, n}, {j, 0, d}];
{Qt, R} = QRDecomposition[V];
DD = DiagonalMatrix[Table[Factorial[i], {i, 0, d}]];
DD.Inverse[R].Qt // TeXForm

$\left( \begin{array}{ccccccc} -\frac{2}{21} & \frac{1}{7} & \frac{2}{7} & \frac{1}{3} & \frac{2}{7} & \frac{1}{7} & -\frac{2}{21} \\ -\frac{3}{28} & -\frac{1}{14} & -\frac{1}{28} & 0 & \frac{1}{28} & \frac{1}{14} & \frac{3}{28} \\ \frac{5}{42} & 0 & -\frac{1}{14} & -\frac{2}{21} & -\frac{1}{14} & 0 & \frac{5}{42} \end{array} \right)$

The first row is the smoothing $P(0)$ filter, the second one the smoothed first order derivative $P'(0)$ filter and the last one the smoothed second order derivative $P''(0)$ filter.

A Julia implementation

Here I present a direct implementation in julia.

We first initialize a Vandermonde matrix:

function vandermonde(halfWindow::Int, polyDeg::Int,T::Type=Float64)
    
    @assert halfWindow>=0
    @assert polyDeg>=0
    
    x=T[i for i in -halfWindow:halfWindow]

    n = polyDeg+1
    m = length(x)
    
    V = Array{T}(m, n)
    
    for i = 1:m
        V[i,1] = T(1)
    end
    for j = 2:n
        for i = 1:m
            V[i,j] = x[i] * V[i,j-1]
        end
    end

    return V
end

We then compute the filter coefficients with the presented formula.

Attention we return the transposed matrix because in Julia it is more convenient to use SG[:,1] which is a julia vector. SG[1,:] would be a (1,n) julia matrix.

function SG(halfWindow::Int, polyDeg::Int,T::Type=Float64)

    @assert 2*halfWindow>polyDeg
    
    V=vandermonde(halfWindow,polyDeg,T)
    Q,R=qr(V)
    SG=R\Q'

    for i in 1:size(SG,1)
        SG[i,:]*=factorial(i-1)
    end
    
# CAVEAT: returns the transposed matrix

    return SG'
end

The final step to use the filter is to provide function to do the cross-correlation.

I do not want to talk too much about this subroutine because in a future post I will show how to efficiently compute such kind of convolution. Here we use a FFT, but with a short filter it is much more efficient (around 10 times faster) to use a direct computation. I will show how to implement

$\gamma[k]=\sum\limits_i\alpha[i]\beta[k+\lambda i],\ with\ \lambda\in\mathbb{Z}^*$

which can be used to compute discrete and stationary wavelet transform for instance.

One last thing, here we use constant padding to manage the boundaries.

function apply_filter{T}(filter::StridedVector{T},signal::StridedVector{T})

    @assert isodd(length(filter))

    halfWindow = round(Int,(length(filter)-1)/2)
    
    padded_signal = 
        [signal[1]*ones(halfWindow);
         signal;
         signal[end]*ones(halfWindow)]

    filter_cross_signal = conv(filter[end:-1:1], padded_signal)

    return filter_cross_signal[2*halfWindow+1:end-2*halfWindow]
end

Finally I have included a small usage example. To see well the effect of Savitzky-Golay filters, I have over smoothed with a large window width $2.n+1$ , $n=20$

using Winston

s=readdlm("signal.txt")[:,1]

sg=SG(20,3) # halt-window, polynomal degree

#________________

smoothed=apply_filter(sg[:,1],s)

plot(s,"r")
oplot(smoothed)
title("Smoothed")
savefig("smoothed.png")

#________________

smoothed_d1=apply_filter(sg[:,2],s)

plot(smoothed_d1)
title("Smoothed derivative")
savefig("smoothed_d1.png")

#________________

smoothed_d2=apply_filter(sg[:,3],s)

plot(smoothed_d2)
title("Smoothed 2-derivative")
savefig("smoothed_d2.png")

Here is the resulting plots:

Final word

To reproduce these figures you can find the complete code on github.

Table of Contents

1 Context

2 @enum based approach

3 Type based solution

4 Default value and keyword argument

4.1 Default value

4.2 Keyword argument

5 Conclusions

6 Annex: asm code for func_kwa(tag=B)

Disclaimer

Some notations

Goal

First step, no boundary extension

Second step, boundary extensions

Periodic case

Mirror Symmetry case

Boundary extensions

Implementation

Preamble

Index translation / domain definition

Definition of

Julia

Auxiliary subroutines

Boundary extensions

Main subroutine

In C/C++

Usages examples

Basic usages

Performance?

Discussion

Optimization/performance

Applications

Code

Complement: more domains

The domain

The domain

Derivation of the formula

Symbolic computation to check that it works

A Julia implementation

Final word

2 `@enum` based approach

3 `Type` based solution

6 Annex: asm code for `func_kwa(tag=B)`

Definition of $\alpha\_\text{offset}$

The $\Omega^\gamma_2$ domain

The $\Omega^\beta_{2'}$ domain