Remove more code! (#134)

Author: amontoison

Project.toml

@@ -4,19 +4,18 @@ version = "0.1.0"
 authors = ["Caleb Derrickson (@CalebDerrickson), Alexis Montoison (@amontoison)", "Chris Geoga (@cgeoga)"]
 
 [deps]
-KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [weakdeps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Vecchia = "8d73829f-f4b0-474a-9580-cecc8e084068"
 
 [extensions]
-VecchiaMLECUDAExt = "CUDA"
+VecchiaMLECUDAExt = ["CUDA", "KernelAbstractions"]
 VecchiaMLEVecchiaExt = ["StaticArrays", "Vecchia"]
 
 [compat]
@@ -33,7 +32,6 @@ NLPModels = "0.21.5"
 NLPModelsIpopt = "0.10.4, 0.11"
 NLPModelsJuMP = "0.13.2"
 NLPModelsTest = "0.10.2"
-Random = "1.10"
 SparseArrays = "1.10"
 Test = "1.10"
 UnoSolver = "0.1.9"

README.md

@@ -22,36 +22,3 @@ julia> ]
 pkg> add https://github.com/exanauts/VecchiaMLE.jl.git
 pkg> test VecchiaMLE
 ```
-
-## Configuration
-
-The struct `VecchiaMLEInput` is available for the user to specify their analysis. At the minimum, we require the following:
-
-```
-* k::Int                      # Number of conditioning points per point for the Vecchia Approximation.
-* samples::Matrix{Float64}    # Matrix of samples (each row is a sample).
-```
-Other options may be passed as keyword arguments. Such options are:
-
-```
-* ptset::AbstractVector     # The locations of the analysis. May be passed as a matrix or vector of vectors.
-* rowsL::AbstractVector     # The sparsity pattern rows of L if the user gives one. MUST BE IN CSC FORMAT! 
-* colptrL::AbstractVector   # The column pointer of L if the user gives one. MUST BE IN CSC FORMAT!
-* sparsityGeneration        # The method by which to generate a sparsity pattern. See SPARSITY_GEN.
-* metric::Distances.metric  # The metric by which nearest neighbors are determined. Defaults to Euclidean.
-```
-
-## Usage
-
-Once the user has initialized an instance of `VecchiaMLEInput`, they may pass it to `sparsity_pattern` to recover the inverse cholesky.
-
-## Getting samples from a covariance matrix
-
-```
-n = 100                                                            # Dimension of Covariance matrix
-number_of_samples = 100                                            # Analysis samples number
-params = [5.0, 0.2, 2.25]                                          # Parameters for Matern matrix. See BesselK.jl. 
-ptset = VecchiaMLE.generate_safe_xyGrid(n)                         # Generates 100 equiparitioned points on [0, 10] x [0, 10].
-MatCov = VecchiaMLE.generate_MatCov(params, ptset)                 # size n x n. 
-samples = VecchiaMLE.generate_samples(MatCov, number_of_samples)   # Generate samples. 
-```

docs/Project.toml

@@ -1,7 +1,10 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 VecchiaMLE = "499233a1-c46f-47cc-ab67-5a37788e9870"
+Vecchia = "8d73829f-f4b0-474a-9580-cecc8e084068"
 
 [compat]
 Documenter = "1.0"

docs/make.jl

@@ -8,10 +8,7 @@ makedocs(
                            prettyurls = get(ENV, "CI", nothing) == "true",
                            collapselevel = 1),
   sitename = "VecchiaMLE.jl",
-  pages = ["Home" => "index.md",
-    "API" => "api.md",
-    "Vecchia model" => "vecchia_model.md",
-    "Tutorials" => "how_to_run.md"]
+  pages = ["Home" => "index.md"],  # "Tutorials" => "vecchia_model.md"
 )
 
 deploydocs(

docs/src/api.md

@@ -6,11 +6,7 @@ using NLPModelsIpopt
 
 n = 400
 number_of_samples = 100
-
-params = [5.0, 0.2, 2.25, 0.25]
-ptset = VecchiaMLE.generate_safe_xyGrid(n)
-MatCov = generate_MatCov(params, ptset)
-samples = generate_samples(MatCov, number_of_samples)
+samples = ...
 
 P = ones(n, n)
 P = tril(P)
@@ -29,11 +25,7 @@ using NLPModelsIpopt
 
 n = 400
 number_of_samples = 100
-
-params = [5.0, 0.2, 2.25, 0.25]
-ptset = VecchiaMLE.generate_safe_xyGrid(n)
-MatCov = generate_MatCov(params, ptset)
-samples = generate_samples(MatCov, number_of_samples)
+samples = ...
 
 P = ones(n, n)
 P = triu(P)
@@ -43,3 +35,50 @@ nlp_U = VecchiaModel(I, J, samples; format=:coo, uplo=:U)
 output = ipopt(nlp_U)
 U = recover_factor(nlp_U, output.solution)
 ```
+
+```@example Vecchia
+using VecchiaMLE
+using Vecchia
+using StaticArrays
+using LinearAlgebra
+
+# Generate some fake data with an exponential covariance function.
+# Note that each _column_ of sim is an iid replicate, in keeping with formatting
+# standards of the many R packages for GPs. In this demonstration, n is small
+# and the number of replicates is large to demonstrate asymptotic correctness.
+pts = rand(SVector{2,Float64}, 100)
+z   = randn(length(pts), 1000)
+K   = Symmetric([exp(-norm(x-y)) for x in pts, y in pts])
+sim = cholesky(K).L * z
+
+# Create a VecchiaModel using the options specified by Vecchia.jl, in
+# particular choosing an ordering (in this case RandomOrdering()) and a
+# conditioning set design (in this case KNNConditioning(10)). This also returns
+# the permutation for you to use with subsequent data operations. See the docs
+# and myriad extensions of Vecchia.jl for more information on ordering and
+# conditioning set design options.
+(perm, nlp) = VecchiaModel(pts, sim, RandomOrdering(), KNNConditioning(10);
+                           lvar_diag=fill(inv(sqrt(1.0)),  length(pts)),
+                           uvar_diag=fill(inv(sqrt(1e-2)), length(pts)),
+                           lambda=1e-3)
+
+# Now bring in some optimizers and fit the nonparametric model. This gives a U
+# such that Σ^{-1} ≈ U*U', where Σ is the covariance matrix for each column of sim.
+using MadNLP
+result = madnlp(nlp; tol=1e-10)
+U      = UpperTriangular(recover_factor(nlp, result.solution))
+
+# KL divergence from the true covariance:
+K_perm = K[perm, perm]
+kl     = (tr(U'*K_perm*U) - length(pts)) + (-2*logdet(U) - logdet(K_perm))
+
+# compare that with what you get from a parametric Vecchia model using the
+# correct kernel and the same permutation.
+para_vecchia = VecchiaApproximation(pts[perm], (x,y,p)->exp(-norm(x-y)), sim[perm,:];
+                                    ordering=NoPermutation())
+para_U  = rchol(para_vecchia, Float64[]).U
+para_kl = (tr(para_U'*K_perm*para_U) - length(pts)) + (-2*logdet(para_U) - logdet(K_perm))
+
+println("KL divergence with true parametric kernel:  ", para_kl)
+println("KL divergence with nonparametric estimator: ", kl)
+```