Merge pull request #13 from terasakisatoshi/add-example-visualize-tan…

…gent-sp

Add example visualize tangent sp

Project.toml

@@ -6,6 +6,7 @@ version = "0.2.1"
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Example = "7876af07-990d-54b4-ab0e-23690620f79a"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"

experiments/notebook/tangent_space.jl

@@ -0,0 +1,91 @@
+# -*- coding: utf-8 -*-
+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .jl
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.3.3
+#   kernelspec:
+#     display_name: Julia 1.3.1
+#     language: julia
+#     name: julia-1.3
+# ---
+
+using Plots
+using Distributions
+using LinearAlgebra
+using ForwardDiff
+
+# +
+struct Point{T<:Real}
+    x::T
+    y::T
+end
+
+
+Point(xy::Vector{T}) where {T} = Point{T}(T.(xy)...)
+
+# +
+struct Tp
+    zfunc::Function
+    v::Vector{Float64}
+    w::Vector{Float64}
+    p::Point{Float64}
+    invA::Matrix{Float64}
+end
+
+function Tp(zfunc, p_xy)
+    zfuncwrap(xy) = zfunc(xy...)
+    g = xy -> ForwardDiff.gradient(zfuncwrap, xy)
+    ∇g_x, ∇g_y = g(p_xy)
+
+    v = [
+        1.0
+        0.0
+        ∇g_x
+    ]
+
+    w = [
+        0.0
+        1.0
+        ∇g_x
+    ]
+
+    v = normalize(v)
+    w = normalize(w)
+
+    A = [
+        v[1] w[1]
+        v[2] w[2]
+    ]
+
+    Tp(zfunc, v, w, Point(p_xy), inv(A))
+end
+
+function tangent_plane(tp::Tp, x, y)
+    z_vec = [tp.v[3] tp.w[3]] * tp.invA * [x, y]
+    z = z_vec[1] + tp.zfunc(tp.p.x, tp.p.y)
+    return z
+end
+# -
+
+d = MvNormal([0.5, 0.5], I(2))
+zfunc(x, y) = pdf(d, [x, y])
+tp = Tp(zfunc, [0.0, 0.0])
+
+function visualize(tp::Tp; camera = (0, 30))
+    tp_xs = range(tp.p.x - 0.5, tp.p.x + 0.5, length = 10)
+    tp_ys = range(tp.p.y - 0.5, tp.p.y + 0.5, length = 10)
+    xs = range(tp.p.x - 3.0, tp.p.x + 3.0, length = 20)
+    ys = range(tp.p.y - 3.0, tp.p.y + 3.0, length = 20)
+    p = plot()
+    wireframe!(p, tp_xs, tp_ys, (x, y) -> tangent_plane(tp, x, y), camera = camera)
+    surface!(p, xs, ys, (x, y) -> tp.zfunc(x, y), camera = camera)
+    return p
+end
+
+@gif for θ in 1:180
+    visualize(tp, camera = (θ, 30))
+end