Merge remote-tracking branch 'origin/master' into JS/fix-rowindex

.buildkite/pipeline.yml

@@ -0,0 +1,37 @@
+steps:
+  - label: "Nvidia GPUs -- RegularizedLeastSquares.jl"
+    plugins:
+      - JuliaCI/julia#v1:
+          version: "1.10"
+    agents:
+      queue: "juliagpu"
+      cuda: "*"
+    command: |
+      julia --color=yes --project -e '
+      using Pkg
+      Pkg.add("TestEnv")
+      using TestEnv
+      TestEnv.activate();
+      Pkg.add("CUDA")
+      Pkg.instantiate()
+      include("test/gpu/cuda.jl")'
+    timeout_in_minutes: 30
+
+  - label: "AMD GPUs -- RegularizedLeastSquares.jl"
+    plugins:
+      - JuliaCI/julia#v1:
+          version: "1.10"
+    agents:
+      queue: "juliagpu"
+      rocm: "*"
+      rocmgpu: "*"
+    command: |
+      julia --color=yes --project -e '
+      using Pkg
+      Pkg.add("TestEnv")
+      using TestEnv
+      TestEnv.activate();
+      Pkg.add("AMDGPU")
+      Pkg.instantiate()
+      include("test/gpu/rocm.jl")'
+    timeout_in_minutes: 30

.gitignore

@@ -3,5 +3,6 @@
 *.jl.mem
 docs/build/
 docs/site/
+docs/src/generated/
 
 Manifest.toml

Project.toml

@@ -1,7 +1,7 @@
 name = "RegularizedLeastSquares"
 uuid = "1e9c538a-f78c-5de5-8ffb-0b6dbe892d23"
 authors = ["Tobias Knopp <[email protected]>"]
-version = "0.14.0"
+version = "0.16.5"
 
 [deps]
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -16,20 +16,33 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 VectorizationBase = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
+[weakdeps]
+GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+
+
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
 
 [compat]
 IterativeSolvers = "0.9"
 julia = "1.9"
 StatsBase = "0.33, 0.34"
-FLoops = "0.2"
 VectorizationBase = "0.19, 0.21"
-LinearOperatorCollection = "1.2"
-LinearOperators = "2.3.3"
 FFTW = "1.0"
+FLoops = "0.2"
+GPUArrays = "8, 9, 10"
+CUDA = "4, 5"
+JLArrays = "0.1.2"
+LinearOperatorCollection = "2"
+LinearOperators = "2.3.3"
 
 [targets]
-test = ["Test", "Random", "FFTW"]
+test = ["Test", "Random", "FFTW", "JLArrays"]
+
+[extensions]
+RegularizedLeastSquaresGPUArraysExt = "GPUArrays"
+RegularizedLeastSquaresCUDAExt = "CUDA"

docs/Project.toml

@@ -1,10 +1,15 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ImageGeoms = "9ee76f2b-840d-4475-b6d6-e485c9297852"
+ImagePhantoms = "71a99df6-f52c-4da1-bd2a-69d6f37f3252"
+JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
 LinearOperatorCollection = "a4a2c56f-fead-462a-a3ab-85921a5f2575"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+RadonKA = "86de8297-835b-47df-b249-c04e8db91db5"
 RegularizedLeastSquares = "1e9c538a-f78c-5de5-8ffb-0b6dbe892d23"
-UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a"
 
 [compat]
-Documenter = "1.0, 1.1, 1.2"
+Documenter = "1"

docs/make.jl

@@ -1,27 +1,51 @@
-using Documenter, RegularizedLeastSquares, LinearOperatorCollection, Wavelets
+using Documenter, Literate, RegularizedLeastSquares
+
+# Generate examples
+OUTPUT_BASE = joinpath(@__DIR__(), "src/generated")
+INPUT_BASE = joinpath(@__DIR__(), "src/literate")
+for (_, dirs, _) in walkdir(INPUT_BASE)
+    for dir in dirs
+        OUTPUT = joinpath(OUTPUT_BASE, dir)
+        INPUT = joinpath(INPUT_BASE, dir)
+        for file in filter(f -> endswith(f, ".jl"), readdir(INPUT))
+            Literate.markdown(joinpath(INPUT, file), OUTPUT)
+        end
+    end
+end
 
 makedocs(
     format = Documenter.HTML(;
         prettyurls=get(ENV, "CI", "false") == "true",
         canonical="https://github.com/JuliaImageRecon/RegularizedLeastSquares.jl",
         assets=String[],
+        collapselevel=1,
     ),
     repo="https://github.com/JuliaImageRecon/RegularizedLeastSquares.jl/blob/{commit}{path}#{line}",
     modules = [RegularizedLeastSquares],
     sitename = "RegularizedLeastSquares.jl",
     authors = "Tobias Knopp, Mirco Grosser, Martin Möddel, Niklas Hackelberg, Andrew Mao, Jakob Assländer",
     pages = [
         "Home" => "index.md",
-        "Getting Started" => "gettingStarted.md",
+        "Getting Started" => "generated/examples/getting_started.md",
+        "Examples" => Any["Compressed Sensing" => "generated/examples/compressed_sensing.md", "Computed Tomography" => "generated/examples/computed_tomography.md"],
         "Solvers" => "solvers.md",
-        "Regularization" => "regularization.md",
-        "API" => Any["Solvers" => "API/solvers.md",
+        "Regularization" => "generated/explanations/regularization.md",
+        "How to" => Any[
+            "Weighting" => "generated/howto/weighting.md",
+            "Normal Operator" => "generated/howto/normal_operator.md",
+            "Multi-Threading" => "generated/howto/multi_threading.md",
+            "GPU Acceleration" => "generated/howto/gpu_acceleration.md",
+            "Efficient Kaczmarz" => "generated/howto/efficient_kaczmarz.md",
+            "Callbacks" => "generated/howto/callbacks.md",
+            "Plug-and-Play Regularization" => "generated/howto/plug-and-play.md"
+        ],
+        "API Reference" => Any["Solvers" => "API/solvers.md",
         "Regularization Terms" => "API/regularization.md"],
 
     ],
     pagesonly = true,
     checkdocs = :none,
-    doctest   = true,
+    doctest   = false,
     doctestfilters = [r"(\d*)\.(\d{4})\d+"]
     )
 

docs/src/API/solvers.md

@@ -5,6 +5,9 @@ REPL one can access this documentation by entering the help mode with `?`
 ## solve!
 ```@docs
 RegularizedLeastSquares.solve!(::AbstractLinearSolver, ::Any)
+RegularizedLeastSquares.init!(::AbstractLinearSolver, ::Any)
+RegularizedLeastSquares.init!(::AbstractLinearSolver, ::AbstractSolverState, ::AbstractMatrix)
+
 ```
 
 ## ADMM
@@ -42,8 +45,11 @@ RegularizedLeastSquares.POGM
 RegularizedLeastSquares.SplitBregman
 ```
 
-## Miscellaneous Functions
+## Miscellaneous
 ```@docs
+RegularizedLeastSquares.solverstate
+RegularizedLeastSquares.solversolution
+RegularizedLeastSquares.solverconvergence
 RegularizedLeastSquares.StoreSolutionCallback
 RegularizedLeastSquares.StoreConvergenceCallback
 RegularizedLeastSquares.CompareSolutionCallback

docs/src/index.md

@@ -4,25 +4,34 @@
 
 ## Introduction
 
-RegularizedLeastSquares.jl is a Julia package for solving large scale linear systems
-using different types of algorithms. Ill-conditioned problems arise in many areas of practical interest. To solve these problems, one often resorts to regularization techniques and non-linear problem formulations. This packages provides implementations for a variety of solvers, which are used in fields such as MPI and MRI.
+RegularizedLeastSquares.jl is a Julia package for solving large linear systems using various types of algorithms. Ill-conditioned problems arise in many areas of practical interest. Regularisation techniques and nonlinear problem formulations are often used to solve these problems. This package provides implementations for a variety of solvers used in areas such as MPI and MRI. In particular, this package serves as the optimization backend of the Julia packages [MPIReco.jl](https://github.com/MagneticParticleImaging/MPIReco.jl) and [MRIReco.jl](https://github.com/MagneticResonanceImaging/MRIReco.jl).
 
-The implemented methods range from the $l_2$-regularized CGNR method to more general optimizers such as the Alternating Direction of Multipliers Method (ADMM) or the Split-Bregman method.
+The implemented methods range from the $l^2_2$-regularized CGNR method to more general optimizers such as the Alternating Direction of Multipliers Method (ADMM) or the Split-Bregman method.
 
-For convenience, implementations of popular regularizers, such as $l_1$-regularization and TV regularization, are provided. On the other hand, hand-crafted regularizers can be used quite easily. For this purpose, a `Regularization` object needs to be build. The latter mainly contains the regularization parameter and a function to calculate the proximal map of a given input.
+For convenience, implementations of popular regularizers, such as $l_1$-regularization and TV regularization, are provided. On the other hand, hand-crafted regularizers can be used quite easily.
 
 Depending on the problem, it becomes unfeasible to store the full system matrix at hand. For this purpose, RegularizedLeastSquares.jl allows for the use of matrix-free operators. Such operators can be realized using the interface provided by the package [LinearOperators.jl](https://github.com/JuliaSmoothOptimizers/LinearOperators.jl).
 Other interfaces can be used as well, as long as the product `*(A,x)` and the adjoint `adjoint(A)` are provided. A number of common matrix-free operators are provided by the package [LinearOperatorColection.jl](https://github.com/JuliaImageRecon/LinearOperatorCollection.jl).
 
-## Installation
+## Features
 
-Within Julia, use the package manager:
-```julia
-using Pkg
-Pkg.add("RegularizedLeastSquares")
-```
-This adds the latest release of the package is added. To install a different version, please consult the [Pkg documentation](https://pkgdocs.julialang.org/dev/managing-packages/#Adding-packages).
+* Variety of optimization algorithms optimized for least squares problems
+* Support for matrix-free operators
+* GPU support
 
 ## Usage
 
   * See [Getting Started](@ref) for an introduction to using the package
+
+## See also
+
+Packages:
+* [ProximalAlgorithms.jl](https://github.com/JuliaFirstOrder/ProximalAlgorithms.jl)
+* [ProximalOperators.jl](https://github.com/JuliaFirstOrder/ProximalOperators.jl)
+* [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl)
+* [RegularizedOptimization.jl](https://github.com/JuliaSmoothOptimizers/RegularizedOptimization.jl)
+
+Organizations:
+* [JuliaNLSolvers](https://github.com/JuliaNLSolvers)
+* [JuliaSmoothOptimizers](https://github.com/JuliaSmoothOptimizers)
+* [JuliaFirstOrder](https://github.com/JuliaFirstOrder)

docs/src/literate/examples/compressed_sensing.jl

@@ -0,0 +1,69 @@
+# # Compressed Sensing Example
+# In this example we will go through a simple example from the field of Compressed Sensing.
+# In addtion to RegularizedLeastSquares.jl, we will need the packages [LinearOperatorCollection.jl](https://github.com/JuliaImageRecon/LinearOperatorCollection.jl), [ImagePhantoms.jl](https://github.com/JuliaImageRecon/ImagePhantoms.jl), [ImageGeoms.jl](https://github.com/JuliaImageRecon/ImageGeoms.jl) and Random.jl, as well as [CairoMakie.jl](https://docs.makie.org/stable/) for visualization. We can install them the same way we did RegularizedLeastSquares.jl.
+
+# ## Preparing the Inverse Problem
+# To get started, let us generate a simple phantom using the ImagePhantoms and ImageGeom packages:
+
+using ImagePhantoms, ImageGeoms
+N = 256
+image = shepp_logan(N, SheppLoganToft())
+size(image)
+
+# This produces a 256x256 image of a Shepp-Logan phantom. 
+
+# In this example, we consider a problem in which we randomly sample a third of the pixels in the image. Such a problem and the corresponding measurement can be constructed with the packages LinearOperatorCollection and Random:
+
+# We first randomly shuffle the indices of the image and then select the first third of the indices to sample.
+using Random, LinearOperatorCollection
+randomIndices = shuffle(eachindex(image))
+sampledIndices = sort(randomIndices[1:div(end, 3)])
+
+# Afterwards we build a sampling operator which samples the image at the selected indices.
+A = SamplingOp(eltype(image), pattern = sampledIndices , shape = size(image));
+
+# Then we apply the sampling operator to the vectorized image to obtain the sampled measurement vector
+b = A*vec(image);
+
+# To visualize our image we can use CairoMakie:
+using CairoMakie
+function plot_image(figPos, img; title = "", width = 150, height = 150)
+  ax = CairoMakie.Axis(figPos; yreversed=true, title, width, height)
+  hidedecorations!(ax)
+  heatmap!(ax, img)
+end
+fig = Figure()
+plot_image(fig[1,1], image, title = "Image")
+samplingMask = fill(false, size(image))
+samplingMask[sampledIndices] .= true
+plot_image(fig[1,2], image .* samplingMask, title = "Sampled Image")
+resize_to_layout!(fig)
+fig
+
+# As we can see in the right image, only a third of the pixels are sampled. The goal of the inverse problem is to recover the original image from this measurement vector.
+
+# ## Solving the Inverse Problem
+# To recover the image from the measurement vector, we solve the TV-regularized least squares problem:
+# ```math
+# \begin{equation}
+#   \underset{\mathbf{x}}{argmin} \frac{1}{2}\vert\vert \mathbf{A}\mathbf{x}-\mathbf{b} \vert\vert_2^2 + \lambda\vert\vert\mathbf{x}\vert\vert_{\text{TV}} .
+# \end{equation}
+# ```
+
+# For this purpose we build a TV regularizer with regularization parameter $λ=0.01$
+
+using RegularizedLeastSquares
+reg = TVRegularization(0.01; shape=size(image));
+
+# We will use the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) to solve our inverse problem. Thus, we build the corresponding solver
+
+solver = createLinearSolver(FISTA, A; reg=reg, iterations=20);
+
+# and apply it to our measurement vector
+img_approx = solve!(solver,b)
+
+# To visualize the reconstructed image, we need to reshape the result vector to the correct shape. Afterwards we can use CairoMakie again:
+img_approx = reshape(img_approx,size(image));
+plot_image(fig[1,3], img_approx, title = "Reconstructed Image")
+resize_to_layout!(fig)
+fig