diff --git a/.gitignore b/.gitignore
index 80516d7..5ac67ce 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,3 @@
 Manifest.toml
 .DS_Store
+docs/build/
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index ad24f56..0049dd1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,6 +3,7 @@ uuid = "3630a16b-0f2f-4d88-afbf-c7d59eccf553"
 version = "0.1.0"
 
 [deps]
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
@@ -12,9 +13,9 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
+JuliaFormatter = "1.0.33"
 Manifolds = "0.11.20"
 ManifoldsBase = "2.3.5"
 Manopt = "0.5.37"
@@ -22,3 +23,9 @@ ProgressMeter = "1.11.0"
 RecursiveArrayTools = "4.3"
 TensorOperations = "5.6"
 julia = "1.10"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]
diff --git a/src/api.jl b/src/api.jl
index de87f6e..6535640 100644
--- a/src/api.jl
+++ b/src/api.jl
@@ -1,4 +1,4 @@
-# User-facing convenience API: tucker, save_result, load_result
+# User-facing convenience API: save_result, load_result
 export save_result, load_result
 
 
@@ -8,10 +8,36 @@ include("api/nncpd.jl")
 include("api/btd.jl")
 include("api/tucker.jl")
 
+
 """
     save_result(path::AbstractString, result)
 
-Save a result to a file.
+Save a TensorKitchen result object to `path` for later use.
+
+* This uses Julia's built-in `Serialization` format, so the file is a Julia-native binary file. 
+* It is intended for saving results during local experiments, benchmarks, and development workflows.
+* Common file extensions are `.jls`, `.julia`, `.bin`, or `.tkresult`; the extension is only a convention.
+
+# Examples
+
+```julia
+A = randn(20, 15, 10)
+res = cpd(A, 35)
+
+save_result("cpd_rank35.jls", res)
+```
+
+* You can also save richer records by serializing a `NamedTuple` of the result and additional metadata.
+
+```julia
+record = (
+    method = :cpd,
+    rank = 35,
+    input_size = size(A),
+    result = res,
+)
+save_result("experiment_cpd_rank35.jls", record)
+``` 
 """
 function save_result(path::AbstractString, result)
     open(path, "w") do io
@@ -19,10 +45,36 @@ function save_result(path::AbstractString, result)
     end
     return path
 end
+
 """
     load_result(path::AbstractString)
 
-Load a result from a file
+Load a previously saved TensorKitchen result or experiment record from `path` for later use.
+
+* The file must be written with `save_result`, or otherwise created using Julia's `Serialization.serialize`.
+
+# Examples
+
+```julia
+res = load_result("cpd_rank35.jls")
+
+weights(res)
+factors(res)
+reconstruct(res)
+```
+If the saved object was an experiment record, you can access the stored fields:
+
+```julia
+record = load_result("experiment_cpd_rank35.jls")
+
+record.method
+record.rank
+record.input_size
+record.result
+record.result.weights
+record.result.factors
+reconstruct(record.result)
+```
 """
 function load_result(path::AbstractString)
     open(path, "r") do io
diff --git a/src/api/approx.jl b/src/api/approx.jl
index 81aa067..b48b729 100644
--- a/src/api/approx.jl
+++ b/src/api/approx.jl
@@ -130,7 +130,7 @@ For the generic join path:
 
 ##Notes##
 * `:als` is not a solver option for `approx(...)`. However, if `approx(...)` auto-routes to `cpd(...)` or `btd(...)`, then those specialized pipelines may support ALS separately.
-* `warm_steps` and `warm_init` are not part of the generic `approx(...)` path. Generic joins start from the selected initializer and then use manifold solvers for refinement.
+* `warm_steps` and `warm_init` are not part of the generic `approx(...)` path. Generic joins start from random initial point and then use manifold solvers for refinement.
 * For generic mixed joins, use manifold solvers such as `:rgd`, `:rcg`, or `:lbfgs`.
 """
 function approx(
diff --git a/src/manifolds/join.jl b/src/manifolds/join.jl
index 2a4b1f7..6531ccb 100644
--- a/src/manifolds/join.jl
+++ b/src/manifolds/join.jl
@@ -19,18 +19,18 @@ join_product(base, r) constructs the product-domain manifold used to
 parameterize sums of r structured components. It does not construct the
 image join/secant variety itself.
 """
-
 function _segre_flat_factors(M::Manifolds.Segre)
     dims = factor_dims(M)
     return (Euclidean(1), (Sphere(n - 1) for n in dims)...)
 end
 
 """
-    join_product(base, r) constructs a ProductManifold
+    join_product(base, r) 
 
-Construct rank-`r` join parameter manifold as a plain `ProductManifold`.
-For `Manifolds.Segre`, uses flattened `(Euclidean(1), Sphere, ..., Sphere)` factors
-per component to preserve current CP parameter layout.
+Decides how to expand base into r components:
+* Manifolds.Segre → flattened (Euclidean(1), Sphere, ...) × r (one λ + spheres per rank-1).
+* Manifolds.Tucker / ProductManifold → repeat each factor r times.
+* Generic manifold → ProductManifold(base, base, ..., base).
 """
 function join_product(base::T, r::Int) where {T<:AbstractManifold}
     r >= 2 || throw(ArgumentError("Join rank must be at least 2, got r=$r."))
diff --git a/src/manifolds/secant.jl b/src/manifolds/secant.jl
index 3fe4e34..2e34248 100644
--- a/src/manifolds/secant.jl
+++ b/src/manifolds/secant.jl
@@ -1,6 +1,5 @@
 # manifolds/secant.jl — ProductManifold join constructors for Segre and Tucker families.
 #
-# Pipeline: Manopt uses manifold M for retraction (exp map) and project (egrad→Riemannian grad).
 # - ProductManifold(Euclidean(1), Sphere×...) is the placeholder for Segre (CPD rank-1).
 # - CanonicalCP: Euclidean(r) × Π_m ProductManifold(Sphere(n_m-1), ..., Sphere(n_m-1)) for rank-r CPD.
 # - `join_product(Manifolds.Segre(...), r)` gives flattened CP layout.
@@ -13,8 +12,9 @@ export CPJoin, TuckerJoin, SegreProduct, CanonicalCP
 """
     SegreProduct(dims, r)
 
-Product manifold `Manifolds.Segre(dims) × ... × Manifolds.Segre(dims)` (r factors).
-Each component point uses the `Manifolds.Segre` layout `[[λ], x₁, …, x_d]`.
+* Product manifold `Manifolds.Segre(dims) × ... × Manifolds.Segre(dims)` (r factors).
+* Each component point uses the `Manifolds.Segre` layout `[[λ], x₁, …, x_d]`.
+* `SegreProduct` is used to build a join model for CPD.
 """
 function SegreProduct(dims::NTuple{N,Int}, r::Int) where {N}
     r >= 2 || throw(
@@ -25,15 +25,19 @@ end
 
 SegreProduct(dims::Vector{T}, r::Int) where {T<:Int} = SegreProduct(Tuple(dims), r)
 
-# Canonical rank-r CP: ℝ^r × ∏_m (S^(n_m-1))^r
+
+"""
+    CanonicalCP(dims, r)
+
+* CanonicalCP: Euclidean(r) × Π_m ProductManifold(Sphere(n_m-1), ..., Sphere(n_m-1)) 
+* Used for rank-r CPD.
+"""
 function CanonicalCP(dims::NTuple{N,Int}, r::Int) where {N}
     r >= 2 || throw(ArgumentError("Canonical rank-r CP manifold needs r>=2, got r=$r"))
     mode_factors = ntuple(m -> ProductManifold(ntuple(_ -> Sphere(dims[m] - 1), r)...), N)
     return ProductManifold(Euclidean(r), mode_factors...)
 end
 
-CanonicalCP(dims::Vector{T}, r::Int) where {T<:Int} = CanonicalCP(Tuple(dims), r)
-
 CPJoin(M::Manifolds.Segre, r::Int) = join_product(M, r)
 CPJoin(dims::NTuple{N,T}, r::Int) where {N,T<:Int} = CPJoin(Manifolds.Segre(dims), r)
 CPJoin(dims::Vector{T}, r::Int) where {T<:Int} = CPJoin(Tuple(dims), r)