text:
## support / pmf for Base types ##############################################

"""
    support(distribution)

Return a collection with unique elements containing at least all
the possible outcomes from `distribution`. For numbers and some other types,
it is sorted.

# Examples
```jldoctest
julia> support(1:3)
1:3

julia> support([4, 2, 4])
2:4

julia> support(([3], [1]))
2-element Array{Array{Int64,1},1}:
 [1]
 [3]
```
"""
function support end

issortable(::Type) = false

for T in (Real,AbstractChar,String)
    @eval begin
        issortable(::Type{<:$T}) = true
        issortable(::Type{<:Tuple{Vararg{$T}}}) = true
        issortable(::Type{<:AbstractArray{<:$T}}) = true
    end
end

support(x::Union{AbstractArray{<:Integer},Tuple{Vararg{Integer}}}) =
    UnitRange(extrema(x)...)

support(x::Union{AbstractArray,Tuple})  =
    issortable(eltype(x)) ?
        unique!(sort!(vec(collect(x)))) :
        unique!(vec(collect(x)))

"""
    pmf(distribution, x)

Return the evaluation of the probability mass function of a (discrete)
`distribution` at `x`.

# Examples
```jldoctest
julia> [pmf([1, 2, 4, 1], x) for x=0:4]
5-element Array{Float64,1}:
 0.0
 0.5
 0.25
 0.0
 0.25
```
"""
pmf(distribution, x)

pmf_count(distribution) = 1.0

function pmf(A::Union{AbstractArray,Tuple}, x)
    x ∈ A || return 0.0
    c = count(isequal(x), A)
    c * 1.0 / length(A)
end

pmf(A::AbstractRange, x) = Float64(x ∈ A) / length(A)

function pmf(A::Union{AbstractArray,Tuple}; normalized::Bool=true)
    d = Dict{eltype(A),Float64}()
    r = normalized ? 1/length(A) : 1.0
    for x in A
        p = get(d, x, 0.0)
        d[x] = p + r
    end
    PMF(A, d, normalized=normalized, count=length(A))
end

## PMF #######################################################################

"""
    pmf(distribution; normalized::Bool=true)

Return the probability mass function of a (discrete) `distribution`
as a `PMF` object, which caches evaluations.
If `normalized` is `false`, the evaluations of the function must be
interpreted not as probabilities but as "weigths", which are computed
according to the specificities of the `distribution`
(the sum of the weights somehow correspond to the size of the sample
space with "multiplicities"). Not all distributions support the
`normalized` keyword.

# Examples
```jldoctest
julia> v = [1, 2, 3, 1]; f = pmf(v);

julia> f(1) # result computed from pmf(v, 1) and cached in f
0.5

julia> f # when displayed, all values are computed and cached
pmf for [1, 2, 3, 1] with support of length 3:
  1 => 0.5
  2 => 0.25
  3 => 0.25

julia> pmf(v, normalized=false)
non-normalized pmf for [1, 2, 3, 1] with support of length 3:
  1 => 2.0
  2 => 1.0
  3 => 1.0
```
"""
pmf(d) = PMF(d)

# a struct to cache values of pmf
# function, and map interface mostly for printing purposes
mutable struct PMF{T,D} <: AbstractDict{T,Float64}
    d::D
    pmf::Dict{T,Float64}
    cached::Bool                      # all values cached
    count::Float64                    # sample count with multiplicities
    support::Union{Nothing,Vector{T}} # !== nothing when keys is sorted
end

function PMF(d)
    T = gentype(d)
    PMF{T,typeof(d)}(d, Dict{T,Float64}(), false, -pmf_count(d), nothing)
end

function PMF(d, probas::Dict{T,Float64};
             normalized::Bool=true, count=nothing) where T
    gentype(d) == T ||
        throw(ArgumentError("distribution and dictionary are incompatibles"))

    if count === nothing
        if normalized
            count = 1.0
        else
            count = sum(values(probas))
        end
    else
        count > 0 || throw(ArgumentError("count must be > 0"))
    end

    resort!(PMF(d, probas, true, Float64(normalized ? -count : count),
                nothing))
end

pmf(f::PMF) = f

isnormalized(f::PMF) = f.count < 0.0 || f.count == 1.0 # use isapprox ?

"""
    normalize!(f::PMF)
    denormalize!(f::PMF)

Make the evaluations of `f` normalized (their sum is approximately 1.0) or
denormalized.

# Examples
```julia-repl
julia> f = pmf(rand(1:3, 3), normalized=false)
Non-normalized pmf for [3, 3, 1] with support of length 2:
  1 => 1.0
  3 => 2.0

julia> normalize!(f)
pmf for [3, 3, 1] with support of length 2:
  1 => 0.333333
  3 => 0.666667

julia> denormalize!(f)
Non-normalized pmf for [3, 3, 1] with support of length 2:
  1 => 1.0
  3 => 2.0
```
"""
normalize!, denormalize!

function normalize!(f::PMF)
    cacheall!(f)
    isnormalized(f) && return f
    replace!(f.pmf) do x
        first(x) => last(x) / f.count
    end
    f.count = -f.count
    f
end

function denormalize!(f::PMF)
    cacheall!(f)
    f.count > 0 && return f
    f.count = -f.count
    replace!(f.pmf) do x
        first(x) => last(x) * f.count
    end
    f
end

function Base.filter!(g, f::PMF)
    cacheall!(f)
    # TODO: optimize (avoid two passes over f.pmf)
    s = 0.0
    replace!(f.pmf) do kv
        if g(kv[1])
            kv
        else
            s += kv[2]
            kv[1] => 0.0
        end
    end
    filter!(x -> x[2] != 0.0, f.pmf)
    s == 0.0 && return f # nothing happened
    if f.count > 0 # non-normalized
        f.count -= s
    else # normalized
        f.count = 1 - s # make non-normalized
        # this allows to keep probability-like values, but which don't
        # sum up to 1.0, so f can't be considered normalized anymore
        # old value of count is totally discarded
    end
    resort!(f)
end

function cacheall!(f::PMF)
    if !f.cached
        for x in support(f.d)
            f(x)
        end
        f.cached = true
        resort!(f)
    end
    f
end

function resort!(f::PMF)
    @assert f.cached
    if issortable(keytype(f))
        f.support = sort!(collect(keys(f.pmf)))
    end
    f
end

function (f::PMF)(x)
    if haskey(f.pmf, x)
        f.pmf[x]
    elseif f.cached
        0.0
    else
        p = pmf(f.d, x)
        if p != 0.0
            f.pmf[x] = p
        end
        p
    end
end

support(f::PMF) = keys(f)

### Dict interface (for printing)

function Base.keys(f::PMF)
    cacheall!(f)
    f.support !== nothing ?
        f.support :
        keys(f.pmf)
end

function Base.values(f::PMF)
    cacheall!(f)
    f.support !== nothing ?
        (f.pmf[x] for x in f.support) :
        values(f.pmf)
end

Base.getindex(f::PMF, x) = cacheall!(f).pmf[x]
Base.length(f::PMF) = length(cacheall!(f).pmf)

function Base.summary(io::IO, f::PMF)
    isnormalized(f) ||
        print(io, "non-normalized ")
    print(io, "pmf for ", f.d, " with support of length ", length(f))
end

function Base.iterate(f::PMF, iter=iterate(keys(cacheall!(f))))
    iter === nothing && return nothing
    (iter[1] => f.pmf[iter[1]]), iterate(keys(f), iter[2])
end

## sampling

# hack: PMF <: AbstractDict for printing purposes, but we want it to be
# also a distribution, so we need to specialize gentype

Sampler(::Type{RNG}, f::PMF, n::Repetition) where {RNG<:AbstractRNG} =
    Sampler(RNG, f.d, n)

Random.gentype(::Type{<:PMF{T}}) where {T} = T


text:
using UnPack: @unpack
using Zygote: gradient
using ChainRulesCore: NoTangent

struct H11{I,A,V} <: JuLIP.AbstractCalculator
    image::I
    AuN::A
    potential::V
end

function H11(jatoms)
    image = ImagePotential(jatoms)
    AuN = AuNCoupling()
    potential = H11DiabaticElement()
    H11(image, AuN, potential)
end

function H11DiabaticElement()
    V11AuO = repulsive(A₁, α₁, rcutoff)
    V11NO = morse(F₁, γ₁, r₁NO)

    potentials = DefaultDict{Tuple, JuLIP.AnalyticFunction}(zero_potential)
    potentials[1,3] = V11AuO
    potentials[2,3] = V11NO

    DiabaticElement(potentials)
end

function JuLIP.energy(calc::H11, at::JuLIP.AbstractAtoms; domain=1:length(at))
    e1 = JuLIP.energy(calc.image, at)
    e2 = JuLIP.energy(calc.AuN, at; domain=domain)
    e3 = JuLIP.energy(calc.potential, at; domain=domain)
    return e1 + e2 + e3
end

function JuLIP.forces(calc::H11, at::JuLIP.AbstractAtoms; kwargs...)
    f1 = JuLIP.forces(calc.image, at)
    f2 = JuLIP.forces(calc.AuN, at; kwargs...)
    f3 = JuLIP.forces(calc.potential, at; kwargs...)
    return f1 + f2 + f3
end

struct ImagePotential{V,T} <: JuLIP.AbstractCalculator
    V::V
    mO::T
    mN::T
    mtotal::T
    iN::Int
    iO::Int
    slab_plane::T
end

function ImagePotential(jatoms)
    iN = findfirst(isequal(7), jatoms.Z)
    iO = findfirst(isequal(8), jatoms.Z)
    mN = jatoms.M[iN]
    mO = jatoms.M[iO]
    mtotal = mO + mN

    gold_positions = jatoms.X[jatoms.Z .== 79]
    slab_plane = maximum(x->x[3], gold_positions)

    V = image(D, C, zimage)
    ImagePotential(V, mO, mN, mtotal, iN, iO, slab_plane)
end

JuLIP.energy(calc::ImagePotential, at::JuLIP.AbstractAtoms) = JuLIP.evaluate(calc, at.X)
JuLIP.forces(calc::ImagePotential, at::JuLIP.AbstractAtoms) = -JuLIP.evaluate_d(calc, at.X)

function JuLIP.evaluate(calc::ImagePotential, R)
    @unpack mO, mN, mtotal, iN, iO = calc
    zcom = (R[iO][3]*mO + R[iN][3]*mN) / mtotal - calc.slab_plane
    return JuLIP.evaluate(calc.V, zcom)
end

function JuLIP.evaluate_d(calc::ImagePotential, R)
    @unpack mO, mN, mtotal, iN, iO = calc
    zcom = (R[iO][3]*mO + R[iN][3]*mN) / mtotal - calc.slab_plane
    image = JuLIP.evaluate_d(calc.V, zcom)
    out = zero(R)
    out[iN] = SVector{3}(0, 0, image * mN / mtotal)
    out[iO] = SVector{3}(0, 0, image * mO / mtotal)
    return out
end

struct AuNCoupling{V,Z} <: JuLIP.SitePotential
    V::V
    zlist::Z
end

function AuNCoupling()
    V(cosθ, r) = B₁*(exp(-2β₁*(r-r₁AuN)) - exp(-2β₁*(rcutoff-r₁AuN))
                - 2cosθ^2*(exp(-β₁*(r-r₁AuN)) - exp(-β₁*(rcutoff-r₁AuN)))
            )
    zlist = JuLIP.Potentials.ZList([:Au, :N, :O])
    AuNCoupling(V, zlist)
end

function JuLIP.evaluate!(tmp, V::AuNCoupling, Rs, Zs, z0)
    Es = 0.0
    i0 = JuLIP.z2i(V, z0)
    if i0 == 2 # Evaluated only at N site
        cosθ = calculate_cosθ(Rs, Zs)
        for (R, Z) in zip(Rs, Zs)
            i = JuLIP.z2i(V, Z)
            if i == 1
                r = norm(R)
                Es += V.V(cosθ, r)
            end
        end
    end
    return Es
end

function calculate_cosθ(Rs, Zs)
    closest = find_closest(Rs, Zs, 8)
    NO_distance = Rs[closest]
    cosθ = NO_distance[3] / norm(NO_distance)
    return cosθ
end

"""
Find the nearest atom with `z == ztarget`.
This is used to make sure the NO molecule interacts only with itself,
not the periodic images.
"""
function find_closest(Rs, Zs, ztarget)
    target_indices = findall(isequal(ztarget), Zs)
    max_dist = 100
    closest = target_indices[1]
    for i in target_indices
        dist = norm(Rs[i])
        if dist < max_dist
            closest = i
            max_dist = dist
        end
    end
    return closest
end

"""
This is the slowest part of the code, a hand-coded derivative would be faster,
or possibly a different autodiff package.
"""
function JuLIP.evaluate_d!(dEs, tmp, V::AuNCoupling, Rs, Zs, z0)
    grad = gradient(x -> JuLIP.evaluate!(tmp, V, x, Zs, z0), Rs)[1]
    if grad === nothing
        for i=1:length(dEs)
            dEs[i] = zero(dEs[i])
        end
    else
        for i=1:length(grad)
            if (grad[i] === nothing) || (grad[i] == NoTangent())
                dEs[i] = zero(dEs[i])
            else
                dEs[i] = grad[i]
            end
        end
    end
    return dEs
end

JuLIP.cutoff(::AuNCoupling) = rcutoff

ext: jl


text:
## R code 9.12
m8.3 <- quap(
    alist(

        log_gdp_std ~ dnorm( mu , sigma ) ,
        mu <- a[cid] + b[cid]*( rugged_std - 0.215 ) ,
        a[cid] ~ dnorm( 1 , 0.1 ) ,
        b[cid] ~ dnorm( 0 , 0.3 ) ,
        sigma ~ dexp( 1 )
    ) , data=dd )
precis( m8.3 , depth=2 )

## R code 9.16
m9.1 <- ulam(
    alist(
        log_gdp_std ~ dnorm( mu , sigma ) ,
        mu <- a[cid] + b[cid]*( rugged_std - 0.215 ) ,
        a[cid] ~ dnorm( 1 , 0.1 ) ,
        b[cid] ~ dnorm( 0 , 0.3 ) ,
        sigma ~ dexp( 1 )
    ) , data=dat_slim , chains=4 , cores=4 )

## R code 9.22
y <- c(-1,1)
set.seed(11)
m9.2 <- ulam(
    alist(
        y ~ dnorm( mu , sigma ) ,
        mu <- alpha ,
        alpha ~ dnorm( 0 , 1000 ) ,
        sigma ~ dexp( 0.0001 )
    ) , data=list(y=y) , chains=3 )

## R code 9.23
precis( m9.2 )

## R code 9.24
set.seed(11)
m9.3 <- ulam(
    alist(
        y ~ dnorm( mu , sigma ) ,
        mu <- alpha ,
        alpha ~ dnorm( 1 , 10 ) ,
        sigma ~ dexp( 1 )
    ) , data=list(y=y) , chains=3 )

## R code 9.26
set.seed(384)
m9.4 <- ulam(
    alist(
        y ~ dnorm( mu , sigma ) ,
        mu <- a1 + a2 ,
        a1 ~ dnorm( 0 , 1000 ),
        a2 ~ dnorm( 0 , 1000 ),
        sigma ~ dexp( 1 )
    ) , data=list(y=y) , chains=3 )
precis( m9.4 )

## R code 9.27
m9.5 <- ulam(
    alist(
        y ~ dnorm( mu , sigma ) ,
        mu <- a1 + a2 ,
        a1 ~ dnorm( 0 , 10 ),
        a2 ~ dnorm( 0 , 10 ),
        sigma ~ dexp( 1 )
    ) , data=list(y=y) , chains=3 )

## R code 9.30
m5.8s <- ulam(
    alist(
        height ~ dnorm( mu , sigma ) ,
        mu <- a + bl*leg_left + br*leg_right ,
        a ~ dnorm( 10 , 100 ) ,
        bl ~ dnorm( 2 , 10 ) ,
        br ~ dnorm( 2 , 10 ) ,
        sigma ~ dexp( 1 )
    ) ,
    data=d, chains=4,
    start=list(a=10,bl=0,br=0.1,sigma=1) )

## R code 9.31
m5.8s2 <- ulam(
    alist(
        height ~ dnorm( mu , sigma ) ,
        mu <- a + bl*leg_left + br*leg_right ,
        a ~ dnorm( 10 , 100 ) ,
        bl ~ dnorm( 2 , 10 ) ,
        br ~ dnorm( 2 , 10 ) ,
        sigma ~ dexp( 1 )
    ) ,
    data=d, chains=4,
    constraints=list(br="lower=0"),
    start=list(a=10,bl=0,br=0.1,sigma=1) )


text:
function get_twister(;
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=0.8,
    contact=true,
    contact_type=:nonlinear,
    spring=0.0,
    damper=0.0,
    num_bodies=5,
    joint_type=:Prismatic,
    height=1.0,
    radius=0.05,
    T=Float64)

    # Parameters
    ex = [1.0; 0.0; 0.0]
    ey = [0.0; 1.0; 0.0]
    ez = [0.0; 0.0; 1.0]
    axes = [ex, ey, ez]

    vert11 = [0.0; 0.0; height / 2.0]
    vert12 = -vert11
    vert = height / 2.0

    # Links
    origin = Origin{T}()
    bodies = [Box(3 * radius, 2 * radius, height, height, color=RGBA(1.0, 0.0, 0.0)) for i = 1:num_bodies]

    # Constraints
    jointb1 = JointConstraint(Floating(origin, bodies[1],
        spring=0.0,
        damper=0.0)) # TODO remove the spring and damper from floating base

    if num_bodies > 1
        joints = [JointConstraint(Prototype(joint_type, bodies[i - 1], bodies[i], axes[i%3+1];
            parent_vertex=vert12,
            child_vertex=vert11,
            spring=spring,
            damper=damper)) for i = 2:num_bodies]
        joints = [jointb1; joints]
    else
        joints = [jointb1]
    end

    if contact
        n = num_bodies
        normal = [[0.0; 0.0; 1.0] for i = 1:n]
        friction_coefficient = friction_coefficient * ones(n)
        contacts1 = contact_constraint(bodies[1], normal[1],
            friction_coefficient=friction_coefficient[1],
            contact_point=vert11,
            contact_type=contact_type) # to avoid duplicating the contact points
        contacts2 = contact_constraint(bodies, normal,
            friction_coefficient=friction_coefficient,
            contact_points=fill(vert12, n),
            contact_type=contact_type)
        mech = Mechanism(origin, bodies, joints, [contacts1; contacts2],
            gravity=gravity,
            timestep=timestep)
    else
        mech = Mechanism(origin, bodies, joints,
            gravity=gravity,
            timestep=timestep,
            spring=spring,
            damper=damper)
    end
    return mech
end

function initialize_twister!(mechanism::Mechanism{T};
    base_position=[0.0, -1.0, 0.0],
    base_orientation=RotX(0.6 * π),
    base_linear_velocity=zeros(3),
    base_angular_velocity=zeros(3),
    relative_linear_velocity=zeros(3),
    relative_angular_velocity=zeros(3)) where T

    bodies = mechanism.bodies
    pbody = bodies[1]
    h = 1.0
    vert11 = [0.0; 0.0; h / 2.0]
    vert12 = -vert11

    # set position and velocities
    set_maximal_configurations!(mechanism.origin, pbody,
        child_vertex=base_position,
        Δq=base_orientation)
    set_maximal_velocities!(pbody,
        v=base_linear_velocity,
        ω=base_angular_velocity)

    previd = pbody.id
    for body in mechanism.bodies[2:end]
        set_maximal_configurations!(get_body(mechanism, previd), body,
            parent_vertex=vert12,
            child_vertex=vert11)
        set_maximal_velocities!(get_body(mechanism, previd), body,
            parent_vertex=vert12,
            child_vertex=vert11,
            Δv=relative_linear_velocity,
            Δω=relative_angular_velocity)
        previd = body.id
    end
end


text:
using ParallelAccelerator

@acc function mttkrp(I::Int, J::Int,
                     K::Int, L::Int,
                     B::Array{Float64, 3}, C::Matrix{Float64},
                     D::Matrix{Float64})

    A = Matrix{Float64}(I, J)
    for i = 1:I
        for j = 1:J
            s = 0.0
            for k = 1:K
                t = 0.0
                for l = 1:L
                    t = t + B[i, k, l] * D[l, j]
                end
                s = s + t * C[k, j]
            end
            A[i, j] = s
        end
    end
    return A
end

I, J, K, L = (200, 200, 200, 200)

B = randn(I, K, L)
C = randn(K, J)
D = randn(L, J)

res1 = mttkrp(I, J, K, L, B, C, D)

@inbounds function mttkrp_o(I::Int, J::Int,
                     K::Int, L::Int,
                     B::Array{Float64, 3}, C::Matrix{Float64},
                     D::Matrix{Float64})

    tic()
    A = Matrix{Float64}(I, J)
    for i = 1:I
        for j = 1:J
            s = 0.0
            for k = 1:K
                t = 0.0
                for l = 1:L
                    t = t + B[i, k, l] * D[l, j]
                end
                s = s + t * C[k, j]
            end
            A[i, j] = s
        end
    end
    println("Julia ", toq())
    return A
 end

res2 = mttkrp_o(I, J, K, L, B, C, D)

for i = 1:I
    exit = false
    for j = 1:J
        if res2[i, j] != res1[i, j]
            println("different $(res2[i, j]) $(res1[i,j]) $(res1[j, i])")
            exit = true
            break
        end
    end
    if exit
        break
    end
end

function compare(A, B)
    for i in 1:I
        for j in 1:J
            if A[i, j] != B[i, j]
                println("different $((i, j)) : $(A[i, j]) != $(B[i, j])")
                return false
            end
        end
    end
    return true
end

println(compare(res1, res2))

max_stars_repo_name: malekbr/julia-to-tiramisu

max_issues_repo_name: malekbr/julia-to-tiramisu

max_forks_repo_name: malekbr/julia-to-tiramisu


text:
module StatLearnTest
using OnlineStats, Distributions, Base.Test

function linearmodeldata(n, p, corr = 0)
    # linear model data with correlated predictors
    V = zeros(p, p)
    for j in 1:p, i in 1:p
        V[i, j] = corr^abs(i - j)
    end
    x = rand(MvNormal(ones(p), V), n)'
    β = vcat(1.:5, zeros(p-5))
    y = x*β + randn(n)
    (β, x, y)
end

@testset "StatLearn" begin
    n, p = 500, 5
    x = randn(n, p)
    β = collect(linspace(-1, 1, p))
    β_with_intercept = vcat(0.0, β)
    xβ = x*β

    generate(::LinearRegression, xβ) = xβ + randn(size(xβ, 1))
    generate(::L1Regression, xβ) = xβ + randn(size(xβ, 1))
    generate(::LogisticRegression, xβ) = [rand(Bernoulli(1 / (1 + exp(-η)))) for η in xβ]
    generate(::PoissonRegression, xβ) = [rand(Poisson(exp(η))) for η in xβ]
    generate(::QuantileRegression, xβ) = xβ + randn(size(xβ, 1))
    generate(::SVMLike, xβ) = [rand(Bernoulli(1 / (1 + exp(-η)))) for η in xβ]
    generate(::HuberRegression, xβ) = xβ + randn(size(xβ, 1))

    @testset "methods" begin
        y = x*β + randn(n)
        o = StatLearn(x, y)
        @test predict(o, x) ≈ x * o.β + o.β0
        @test coef(o) == vcat(o.β0, o.β)
    end
    @testset "loss" begin
        y = generate(LinearRegression(), xβ)
        o = StatLearn(x, y, LinearRegression())
        @test loss(o, x, y) ≈ .5 * mean(abs2(y - predict(o, x)))

        y = generate(L1Regression(), xβ)
        o = StatLearn(x, y, L1Regression())
        @test loss(o, x, y) ≈ mean(abs(y - predict(o, x)))

        y = generate(LogisticRegression(), xβ)
        o = StatLearn(x, y, LogisticRegression())
        η = o.β0 + x * o.β
        l = mean([-y[i] * η[i] + log(1.0 + exp(η[i])) for i in 1:length(η)])
        @test loss(o, x, y) ≈ l

        y = generate(PoissonRegression(), xβ)
        o = StatLearn(x, y, PoissonRegression())
        η = o.β0 + x * o.β
        # TODO
        # @test loss(o, x, y) ≈ mean(-y .* η + exp(η))

        y = generate(QuantileRegression(.7), xβ)
        o = StatLearn(x, y, QuantileRegression(.7))
        r = y - o.β0 - x * o.β
        @test loss(o, x, y) ≈ mean([r[i] * (o.model.τ - (r[i]<0)) for i in 1:n])

        y = generate(SVMLike(), xβ)
        o = StatLearn(x, y, SVMLike())
        η = o.β0 + x * o.β
        @test loss(o, x, y) ≈ mean([max(0.0, 1.0 - y[i] * η[i]) for i in 1:n])

        y = generate(HuberRegression(4.), xβ)
        o = StatLearn(x, y, HuberRegression(4.))
        δ = o.model.δ
        r = y - o.β0 - x * o.β
        v = [abs(r[i]) < δ? 0.5 * r[i]^2 : δ * (abs(r[i]) - 0.5 * δ) for i in 1:n]
        @test loss(o, x, y) ≈ mean(v)
    end
end

end #module