"""
`x, y = noisy_function(fun, x; noise = 0.01, f_rand = randn)`

Generates a noisy response `y` for the given function `fun`
by adding `noise .* f_randn(length(x))` to the result of `fun(x)`.
"""
function noisy_function(fun::Function, x::AbstractVector{T}; noise::Real = 0.01, f_rand::Function = randn) where T<:Real
    x_vec = collect(x)
    y = fun.(x_vec) .+ noise .* f_rand.()
    x_vec, y
end

"""
`x, y = noisy_sin(n, start, stop; noise = 0.3, f_rand = randn)`

Generates `n` noisy equally spaces samples of a sinus from `start` to `stop`
by adding `noise .* f_randn(length(x))` to the result of `fun(x)`.
"""
function noisy_sin(n::Int = 50, start::Real = 0, stop::Real = 2π; noise::Real = 0.3, f_rand::Function = randn)
    noisy_function(sin, range(start, stop=stop, length=n); noise = noise, f_rand = f_rand)
end

"""
`x, y = noisy_poly(coef, x; noise = 0.01, f_rand = randn)`

Generates a noisy response for a polynomial of degree `length(coef)`
using the vector `x` as input and adding `noise .* f_randn(length(x))` to the result.
The vector `coef` contains the coefficients for the terms of the polynome.
The first element of `coef` denotes the coefficient for the term with
the highest degree, while the last element of `coef` denotes the intercept.
"""
function noisy_poly(coef::AbstractVector{R}, x::AbstractVector{T}; noise::Real = 0.1, f_rand::Function = randn) where {T<:Real,R<:Real}
    n = length(x)
    m = length(coef)
    x_vec = collect(x)
    y = zeros(n)
    @inbounds for i = 1:n
        for k = 1:m
            y[i] += coef[k] * x_vec[i]^(m-k)
        end
    end
    y .+= noise .* f_rand.()
    x_vec, y
end


"""
`x, y = noisy_spiral(n, a, theta, b; noise = 0.01, f_rand = randn)`

Generates `n` noisy responses for a spiral with two labels. Uses the radius, angle
and scaling arguments to space the points in 2D space and adding `noise .* f_randn(n)`
to the response.
"""
function noisy_spiral(n::Int = 97, a::Real = 6.5, theta::Real = 16.0, b::Real=104.0; noise::Real = 0.1, f_rand::Function = randn)
    x = zeros(Float64, (2, 2*n))
    y = zeros(Int, 2*n)
    index = 0:1.0:(n-1)
    for i = 1:n
        _angle = index[i]*pi/theta
    	_radius = a * (b-index[i]) / b
    	x_coord = _radius * sin(_angle)
    	y_coord = _radius * cos(_angle)
    	x[1, i] = x_coord
    	x[2, i] = y_coord
        x[1, n+i] = -(x_coord)
    	x[2, n+i] = -(y_coord)
        y[i] = 1
        y[n+i] = 0
    end
    x[1, :] .+= noise .* f_rand.()
    x[2, :] .+= noise .* f_rand.()
    x, y
end
