diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 0000000..01a77c0
Binary files /dev/null and b/.DS_Store differ
diff --git a/docs/.DS_Store b/docs/.DS_Store
new file mode 100644
index 0000000..4428674
Binary files /dev/null and b/docs/.DS_Store differ
diff --git a/docs/Project.toml b/docs/Project.toml
index 112b3af..67597e4 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -6,6 +6,7 @@
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 ProxTV = "925ea013-038b-5ab6-a1ab-e0849925e528"
+proxTV_jll = "700117f8-5dbb-54dd-9908-6f3eb0e21f87"
 
 [compat]
 Documenter = "1"
diff --git a/src/.DS_Store b/src/.DS_Store
new file mode 100644
index 0000000..1774f9e
Binary files /dev/null and b/src/.DS_Store differ
diff --git a/src/ProxTV.jl b/src/ProxTV.jl
deleted file mode 100644
index 7102764..0000000
--- a/src/ProxTV.jl
+++ /dev/null
@@ -1,11 +0,0 @@
-module ProxTV
-
-"""
-    hi = hello_world()
-A simple function to return "Hello, World!"
-"""
-function hello_world()
-  return "Hello, World!"
-end
-
-end
diff --git a/src/gen/Project.toml b/src/gen/Project.toml
new file mode 100644
index 0000000..8d6c2a6
--- /dev/null
+++ b/src/gen/Project.toml
@@ -0,0 +1,7 @@
+[deps]
+Clang = "40e3b903-d033-50b4-a0cc-940c62c95e31"
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
+proxTV_jll = "700117f8-5dbb-54dd-9908-6f3eb0e21f87"
+
+[compat]
+julia = "1.6"
diff --git a/src/gen/README.md b/src/gen/README.md
new file mode 100644
index 0000000..ced2cd2
--- /dev/null
+++ b/src/gen/README.md
@@ -0,0 +1,9 @@
+# Wrapping headers
+
+This directory contains a script that can be used to automatically generate wrappers from C headers provided by [proxTV](https://github.com/albarji/proxTV).
+This is done using Clang.jl.
+
+## Usage
+
+Either run `julia wrapper.jl` directly, or include it and call the `main()` function.
+Be sure to activate the project environment in this folder (`julia --project`), which will install `Clang.jl` and `JuliaFormatter.jl`.
diff --git a/src/gen/proxtv.toml b/src/gen/proxtv.toml
new file mode 100644
index 0000000..9cc24d7
--- /dev/null
+++ b/src/gen/proxtv.toml
@@ -0,0 +1,8 @@
+[general]
+use_julia_native_enum_type = true
+print_using_CEnum = false
+library_name = "libproxtv"
+
+[codegen]
+use_julia_bool = true
+use_ccall_macro = true
diff --git a/src/gen/wrapper.jl b/src/gen/wrapper.jl
new file mode 100644
index 0000000..935d34c
--- /dev/null
+++ b/src/gen/wrapper.jl
@@ -0,0 +1,41 @@
+# Script to parse proxTV headers and generate Julia wrappers.
+using proxTV_jll
+using Clang
+using Clang.Generators
+using JuliaFormatter
+
+function main()
+  cd(@__DIR__)
+  include_dir = joinpath(proxTV_jll.artifact_dir, "include")
+  headers = [
+    joinpath(include_dir, header) for
+    header in readdir(include_dir) if endswith(header, ".h")
+  ]
+
+  options = load_options(joinpath(@__DIR__, "proxtv.toml"))
+  options["general"]["output_file_path"] = joinpath("..", "libproxtv.jl")
+  options["general"]["output_ignorelist"] =
+    ["mxGetInf", "sign", "min", "max", "dpttrs_", "dpttrf_"]  # example: "RC_OK"
+
+  args = get_default_args()
+  push!(args, "-I$include_dir")
+  push!(args, "-DNOMATLAB")
+
+  ctx = create_context(headers, args, options)
+  build!(ctx)
+
+  path = options["general"]["output_file_path"]
+  code = read(path, String)
+  code = "using proxTV_jll\n\n" * code
+  code = replace(code, "Cdouble" => "Float64")
+  code = replace(code, "Cint" => "Int32")
+  write(path, code)
+
+  format_file(path, YASStyle())
+  return nothing
+end
+
+# If we want to use the file as a script with `julia wrapper.jl`
+if abspath(PROGRAM_FILE) == @__FILE__
+  main()
+end
diff --git a/src/libproxtv.jl b/src/libproxtv.jl
new file mode 100644
index 0000000..aae68e9
--- /dev/null
+++ b/src/libproxtv.jl
@@ -0,0 +1,623 @@
+using proxTV_jll
+
+function LPnorm(x, n, p)
+  @ccall libproxtv.LPnorm(x::Ptr{Float64}, n::Int32, p::Float64)::Float64
+end
+
+function PN_LP1(y, lambda, x, info, n)
+  @ccall libproxtv.PN_LP1(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+  )::Int32
+end
+
+function PN_LP2(y, lambda, x, info, n)
+  @ccall libproxtv.PN_LP2(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+  )::Int32
+end
+
+struct Workspace
+  n::Int32
+  d::Ptr{Ptr{Float64}}
+  maxd::Int32
+  nd::Int32
+  i::Ptr{Ptr{Int32}}
+  maxi::Int32
+  ni::Int32
+  in::Ptr{Float64}
+  out::Ptr{Float64}
+  warm::Cshort
+  warmDual::Ptr{Float64}
+  warmLambda::Float64
+end
+
+function PN_LPinf(y, lambda, x, info, n, ws)
+  @ccall libproxtv.PN_LPinf(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+    ws::Ptr{Workspace},
+  )::Int32
+end
+
+function PN_LPp(y, lambda, x, info, n, p, ws, positive, objGap)
+  @ccall libproxtv.PN_LPp(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+    p::Float64,
+    ws::Ptr{Workspace},
+    positive::Int32,
+    objGap::Float64,
+  )::Int32
+end
+
+function PN_LPp_v2(y, lambda, x, info, n, p, ws, positive)
+  @ccall libproxtv.PN_LPp_v2(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+    p::Float64,
+    ws::Ptr{Workspace},
+    positive::Int32,
+  )::Int32
+end
+
+function LP1_project(y, lambda, x, n, ws)
+  @ccall libproxtv.LP1_project(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    n::Int32,
+    ws::Ptr{Workspace},
+  )::Int32
+end
+
+function LPp_project(y, lambda, x, info, n, p, ws)
+  @ccall libproxtv.LPp_project(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+    p::Float64,
+    ws::Ptr{Workspace},
+  )::Int32
+end
+
+function solveLinearLP(z, n, p, lambda, s)
+  @ccall libproxtv.solveLinearLP(
+    z::Ptr{Float64},
+    n::Int32,
+    p::Float64,
+    lambda::Float64,
+    s::Ptr{Float64},
+  )::Cvoid
+end
+
+function TV(y, lambda, x, info, n, p, ws)
+  @ccall libproxtv.TV(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+    p::Float64,
+    ws::Ptr{Workspace},
+  )::Int32
+end
+
+function PN_TV1(y, lambda, x, info, n, sigma, ws)
+  @ccall libproxtv.PN_TV1(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+    sigma::Float64,
+    ws::Ptr{Workspace},
+  )::Int32
+end
+
+function linearizedTautString_TV1(y, lambda, x, n)
+  @ccall libproxtv.linearizedTautString_TV1(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    n::Int32,
+  )::Int32
+end
+
+function classicTautString_TV1(signal, n, lam, prox)
+  @ccall libproxtv.classicTautString_TV1(
+    signal::Ptr{Float64},
+    n::Int32,
+    lam::Float64,
+    prox::Ptr{Float64},
+  )::Int32
+end
+
+function hybridTautString_TV1(y, n, lambda, x)
+  @ccall libproxtv.hybridTautString_TV1(
+    y::Ptr{Float64},
+    n::Int32,
+    lambda::Float64,
+    x::Ptr{Float64},
+  )::Cvoid
+end
+
+function hybridTautString_TV1_custom(y, n, lambda, x, backtracksexp)
+  @ccall libproxtv.hybridTautString_TV1_custom(
+    y::Ptr{Float64},
+    n::Int32,
+    lambda::Float64,
+    x::Ptr{Float64},
+    backtracksexp::Float64,
+  )::Cvoid
+end
+
+function classicTautString_TV1_offset(signal, n, lam, prox, offset)
+  @ccall libproxtv.classicTautString_TV1_offset(
+    signal::Ptr{Float64},
+    n::Int32,
+    lam::Float64,
+    prox::Ptr{Float64},
+    offset::Float64,
+  )::Int32
+end
+
+function SolveTVConvexQuadratic_a1_nw(n, b, w, solution)
+  @ccall libproxtv.SolveTVConvexQuadratic_a1_nw(
+    n::Int32,
+    b::Ptr{Float64},
+    w::Float64,
+    solution::Ptr{Float64},
+  )::Cvoid
+end
+
+function PN_TV1_Weighted(Y, W, X, info, n, sigma, ws)
+  @ccall libproxtv.PN_TV1_Weighted(
+    Y::Ptr{Float64},
+    W::Ptr{Float64},
+    X::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+    sigma::Float64,
+    ws::Ptr{Workspace},
+  )::Int32
+end
+
+function tautString_TV1_Weighted(y, lambda, x, n)
+  @ccall libproxtv.tautString_TV1_Weighted(
+    y::Ptr{Float64},
+    lambda::Ptr{Float64},
+    x::Ptr{Float64},
+    n::Int32,
+  )::Int32
+end
+
+function PN_TV1_Trend2_Weighted(Y, W, X, info, n, sigma, ws)
+  @ccall libproxtv.PN_TV1_Trend2_Weighted(
+    Y::Ptr{Float64},
+    W::Ptr{Float64},
+    X::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+    sigma::Float64,
+    ws::Ptr{Workspace},
+  )::Int32
+end
+
+function SolveTVConvexQuadratic_a1(n, b, w, solution)
+  @ccall libproxtv.SolveTVConvexQuadratic_a1(
+    n::Int32,
+    b::Ptr{Float64},
+    w::Ptr{Float64},
+    solution::Ptr{Float64},
+  )::Cvoid
+end
+
+function more_TV2(y, lambda, x, info, n)
+  @ccall libproxtv.more_TV2(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+  )::Int32
+end
+
+function morePG_TV2(y, lambda, x, info, n, ws)
+  @ccall libproxtv.morePG_TV2(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+    ws::Ptr{Workspace},
+  )::Int32
+end
+
+function PG_TV2(y, lambda, x, info, n)
+  @ccall libproxtv.PG_TV2(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+  )::Int32
+end
+
+function DR2L1W_TV(M, N, unary, W1, W2, s, nThreads, maxit, info)
+  @ccall libproxtv.DR2L1W_TV(
+    M::Csize_t,
+    N::Csize_t,
+    unary::Ptr{Float64},
+    W1::Ptr{Float64},
+    W2::Ptr{Float64},
+    s::Ptr{Float64},
+    nThreads::Int32,
+    maxit::Int32,
+    info::Ptr{Float64},
+  )::Int32
+end
+
+function GP_TVp(y, lambda, x, info, n, p, ws)
+  @ccall libproxtv.GP_TVp(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+    p::Float64,
+    ws::Ptr{Workspace},
+  )::Int32
+end
+
+function OGP_TVp(y, lambda, x, info, n, p, ws)
+  @ccall libproxtv.OGP_TVp(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+    p::Float64,
+    ws::Ptr{Workspace},
+  )::Int32
+end
+
+function FISTA_TVp(y, lambda, x, info, n, p, ws)
+  @ccall libproxtv.FISTA_TVp(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+    p::Float64,
+    ws::Ptr{Workspace},
+  )::Int32
+end
+
+function FW_TVp(y, lambda, x, info, n, p, ws)
+  @ccall libproxtv.FW_TVp(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+    p::Float64,
+    ws::Ptr{Workspace},
+  )::Int32
+end
+
+function GPFW_TVp(y, lambda, x, info, n, p, ws)
+  @ccall libproxtv.GPFW_TVp(
+    y::Ptr{Float64},
+    lambda::Float64,
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    n::Int32,
+    p::Float64,
+    ws::Ptr{Workspace},
+  )::Int32
+end
+
+function PD2_TV(y, lambdas, norms, dims, x, info, ns, nds, npen, ncores, maxIters)
+  @ccall libproxtv.PD2_TV(
+    y::Ptr{Float64},
+    lambdas::Ptr{Float64},
+    norms::Ptr{Float64},
+    dims::Ptr{Float64},
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    ns::Ptr{Int32},
+    nds::Int32,
+    npen::Int32,
+    ncores::Int32,
+    maxIters::Int32,
+  )::Int32
+end
+
+function proxDykstraTV2DWeighted(y, W, x, ns, nds, info, ncores, maxIters)
+  @ccall libproxtv.proxDykstraTV2DWeighted(
+    y::Ptr{Float64},
+    W::Ptr{Float64},
+    x::Ptr{Float64},
+    ns::Ptr{Int32},
+    nds::Int32,
+    info::Ptr{Float64},
+    ncores::Int32,
+    maxIters::Int32,
+  )::Int32
+end
+
+function DR2_TV(M, N, unary, W1, W2, norm1, norm2, s, nThreads, maxit, info)
+  @ccall libproxtv.DR2_TV(
+    M::Csize_t,
+    N::Csize_t,
+    unary::Ptr{Float64},
+    W1::Float64,
+    W2::Float64,
+    norm1::Float64,
+    norm2::Float64,
+    s::Ptr{Float64},
+    nThreads::Int32,
+    maxit::Int32,
+    info::Ptr{Float64},
+  )::Int32
+end
+
+function CondatChambollePock2_TV(M, N, Y, lambda, X, alg, maxit, info)
+  @ccall libproxtv.CondatChambollePock2_TV(
+    M::Csize_t,
+    N::Csize_t,
+    Y::Ptr{Float64},
+    lambda::Float64,
+    X::Ptr{Float64},
+    alg::Cshort,
+    maxit::Int32,
+    info::Ptr{Float64},
+  )::Int32
+end
+
+function Yang2_TV(M, N, Y, lambda, X, maxit, info)
+  @ccall libproxtv.Yang2_TV(
+    M::Csize_t,
+    N::Csize_t,
+    Y::Ptr{Float64},
+    lambda::Float64,
+    X::Ptr{Float64},
+    maxit::Int32,
+    info::Ptr{Float64},
+  )::Int32
+end
+
+function Kolmogorov2_TV(M, N, Y, lambda, X, maxit, info)
+  @ccall libproxtv.Kolmogorov2_TV(
+    M::Csize_t,
+    N::Csize_t,
+    Y::Ptr{Float64},
+    lambda::Float64,
+    X::Ptr{Float64},
+    maxit::Int32,
+    info::Ptr{Float64},
+  )::Int32
+end
+
+function Yang3_TV(M, N, O, Y, lambda, X, maxit, info)
+  @ccall libproxtv.Yang3_TV(
+    M::Csize_t,
+    N::Csize_t,
+    O::Csize_t,
+    Y::Ptr{Float64},
+    lambda::Float64,
+    X::Ptr{Float64},
+    maxit::Int32,
+    info::Ptr{Float64},
+  )::Int32
+end
+
+function PD_TV(y, lambdas, norms, dims, x, info, ns, nds, npen, ncores, maxIters)
+  @ccall libproxtv.PD_TV(
+    y::Ptr{Float64},
+    lambdas::Ptr{Float64},
+    norms::Ptr{Float64},
+    dims::Ptr{Float64},
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    ns::Ptr{Int32},
+    nds::Int32,
+    npen::Int32,
+    ncores::Int32,
+    maxIters::Int32,
+  )::Int32
+end
+
+function PDR_TV(y, lambdas, norms, dims, x, info, ns, nds, npen, ncores, maxIters)
+  @ccall libproxtv.PDR_TV(
+    y::Ptr{Float64},
+    lambdas::Ptr{Float64},
+    norms::Ptr{Float64},
+    dims::Ptr{Float64},
+    x::Ptr{Float64},
+    info::Ptr{Float64},
+    ns::Ptr{Int32},
+    nds::Int32,
+    npen::Int32,
+    ncores::Int32,
+    maxIters::Int32,
+  )::Int32
+end
+
+function TV1D_denoise(input, output, width, lambda)
+  @ccall libproxtv.TV1D_denoise(
+    input::Ptr{Float64},
+    output::Ptr{Float64},
+    width::Int32,
+    lambda::Float64,
+  )::Cvoid
+end
+
+function TV1D_denoise_tautstring(input, output, width, lambda)
+  @ccall libproxtv.TV1D_denoise_tautstring(
+    input::Ptr{Float64},
+    output::Ptr{Float64},
+    width::Int32,
+    lambda::Float64,
+  )::Cvoid
+end
+
+function dp(n, y, lam, beta)
+  @ccall libproxtv.dp(n::Int32, y::Ptr{Float64}, lam::Float64, beta::Ptr{Float64})::Cvoid
+end
+
+function radialProjection(x, n, norm, lambda)
+  @ccall libproxtv.radialProjection(
+    x::Ptr{Float64},
+    n::Int32,
+    norm::Float64,
+    lambda::Float64,
+  )::Cvoid
+end
+
+function newWorkspace(n)
+  @ccall libproxtv.newWorkspace(n::Int32)::Ptr{Workspace}
+end
+
+function resetWorkspace(ws)
+  @ccall libproxtv.resetWorkspace(ws::Ptr{Workspace})::Cvoid
+end
+
+function getDoubleWorkspace(ws)
+  @ccall libproxtv.getDoubleWorkspace(ws::Ptr{Workspace})::Ptr{Float64}
+end
+
+function getIntWorkspace(ws)
+  @ccall libproxtv.getIntWorkspace(ws::Ptr{Workspace})::Ptr{Int32}
+end
+
+function freeWorkspace(ws)
+  @ccall libproxtv.freeWorkspace(ws::Ptr{Workspace})::Cvoid
+end
+
+function newWorkspaces(n, p)
+  @ccall libproxtv.newWorkspaces(n::Int32, p::Int32)::Ptr{Ptr{Workspace}}
+end
+
+function freeWorkspaces(wa, p)
+  @ccall libproxtv.freeWorkspaces(wa::Ptr{Ptr{Workspace}}, p::Int32)::Cvoid
+end
+
+function compareDoublesReversed(v1, v2)
+  @ccall libproxtv.compareDoublesReversed(v1::Ptr{Cvoid}, v2::Ptr{Cvoid})::Int32
+end
+
+const STOP_PNLP = 0
+
+const STOP_GAP_PNLP = 1.0e-5
+
+const MAX_ITERS_PNLP = 1000
+
+const SIGMA_PNLP = 0.05
+
+const EPSILON_PNLP = 1.0e-15
+
+const MIN_GRD_HESSIAN_PNLP = 1.0e-15
+
+const LPPROJ_PSMALL = 1.002
+
+const LPPROJ_PLARGE = 100
+
+const MIN_STEP_PNLP = 1.0e-10
+
+const STOP_PN = 1.0e-6
+
+const SIGMA = 0.05
+
+const MAX_ITERS_PN = 100
+
+const STOP_MS = 1.0e-5
+
+const STOP_MSSUB = 1.0e-6
+
+const MAX_ITERS_MS = 100
+
+const STOP_TVLP = 1.0e-5
+
+const MAX_ITERS_TVLP = 10000
+
+const MAX_ITERS_TVLPFW = 1000000
+
+const MAX_ITERS_TVLPGPFW = 1000000
+
+const MAX_NOIMP_TVLP = 10
+
+const OBJGAP_LPPROX_TVLP = 1.0e-15
+
+const LAMBDA_STEPS_TVLP = 1
+
+const LAMBDA_REDUCTION_TVLP = 1000.0
+
+const STOP_STEP_TVLP_FW = 1.0e-15
+
+const FW_CYCLES_TVLP = 10
+
+const MAX_NOIMP_TVLP_GPFW = 10FW_CYCLES_TVLP
+
+const MIN_IMP_TVLP = 1.0e-10
+
+const ENFORCED_CYCLES_TVLP_GPFW = 10
+
+const STOP_PD = 1.0e-6
+
+const MAX_ITERS_PD = 35
+
+const MAX_ITERS_CONDAT = 2500
+
+const STOP_CONDAT = 0
+
+const MAX_ITERS_KOLMOGOROV = 2500
+
+const STOP_KOLMOGOROV = 0
+
+const MAX_ITERS_DR = 35
+
+const MAX_ITERS_YANG = 35
+
+const lapack_int = Int32
+
+const DEBUG_N = 10
+
+const N_INFO = 3
+
+const INFO_ITERS = 0
+
+const INFO_GAP = 1
+
+const INFO_RC = 2
+
+const EPSILON = 1.0e-10
+
+const RC_OK = 0
+
+const RC_ITERS = 1
+
+const RC_STUCK = 2
+
+const RC_ERROR = 3
+
+const WS_MAX_MEMORIES = 100