for_alex.jl

import DSP: conv
using LinearAlgebra

function boundaries( ::Val{:zeros}, x::Vector, half_width::Int )
    return (padded = true, z = vcat( zeros(half_width), x, zeros(half_width) ) )
end
function boundaries( ::Val{:repeating}, x::Vector, half_width::Int )
    return (padded = true, z = vcat( repeat( [ first(x) ], half_width), x, repeat( [ last(x) ], half_width) ) )
end
function boundaries( ::Val{:mirroring}, x::Vector, half_width::Int )
    return (padded = true, z = vcat( reverse(x[1:half_width]), x, reverse(x[(end-half_width):end]) ))
end
boundaries( anything, x::Vector, half_width::Int ) = (padded = false, z = x)
boundary_mapper(s::Symbol, x::Vector, half_width::Int) = boundaries( Val(s), x, half_width )

function SavitskyGolay( x::Vector, window::Int, ∂_order::Int, polynomial_order::Int;
                        Δt::Float64 = 1., boundarycondition = :repeating)
    len = length( x )
    half_width = ( ( window - 1 ) / 2 ) |> Int
    @assert ( len > window > polynomial_order ) "Invalid window, or polynomial order."
    @assert ( window > 1 ) "Invalid window size. Must be positive."
    @assert ( window % 2 == 1 ) "Window size must be an odd number."
    @assert ( ∂_order < polynomial_order ) "Polynomial order must be greater then the derivative order"
    A = LinearAlgebra.pinv( (-half_width:half_width) .^ transpose( 0:polynomial_order ) )
    true_ν = Matrix{Float64}(undef, window, window)

    true_ν = A[ 1+∂_order,:] .*= factorial(∂_order) / ( (-Δt) ^ Float64( ∂_order) )
    padded, newx = boundary_mapper( boundarycondition, x, half_width)
    lft_pad,rgt_pad = 1:half_width, (len - half_width + 1) : len
    if padded
        rng = ((2half_width+1):(len + 2half_width))
        return conv(newx, true_ν)[rng]
    else
        return conv(newx, true_ν)
    end
end

using Plots

#IR test
y = zeros(100)
y[50] = 1.0
plot( y, label = "Impulse");
plot!( SavitskyGolay(y, 33, 0, 6; boundarycondition = :zeros),
        label = "SG 0th derivative")

ω = 11

#remove Boundary effects
θ = (1:3000) .* (pi/1000)
y = cos.( θ ) .+ randn( length( θ ) ) / 20

plot(cos.( θ ), label = "analytic sine");
plot!(SavitskyGolay(cos.( θ ), ω, 0, 7; boundarycondition = :mirroring), label = "SG 0th derivative")

plot(-sin.( θ ), label = "analytic 1st derivative");
plot!(SavitskyGolay( cos.( θ ), ω, 1, 7;  Δt = θ.step |> Float64,
                    boundarycondition = :mirroring),
        label = "SG 1st derivative")

plot(-cos.( θ ), label = "analytic 2nd derivative");
plot!(SavitskyGolay(cos.( θ ), ω, 2, 5; Δt = θ.step |> Float64,
                    boundarycondition = :mirroring),
            label = "SG 2nd derivative")

BC_none = SavitskyGolay(y, ω, 0, 3; boundarycondition = :none)
BC_zeros = SavitskyGolay(y, ω, 0, 5; boundarycondition = :zeros)
BC_repeating = SavitskyGolay(y, ω, 0, 3; boundarycondition = :repeating)
BC_mirroring = SavitskyGolay(y, ω, 0, 3; boundarycondition = :mirroring)

plot(y, legend = false);
plot!(BC_none);
plot!(BC_zeros);
plot!(BC_repeating);
plot!(BC_mirroring)

function ⊥(operator::Function, a::AbstractArray, b::AbstractArray)
    a_dim_lens = length( size( a ) )
    operator.(a, reshape( b, ( ones(Int, a_dim_lens)..., size(b)... ) ) )
end

⊥( -, ⊥( ^, 1:3, 1:3 ), 1:3)