--- title: "Spacekime Analytics (Time Complexity and Inferential Uncertainty)" subtitle: "Kime Dynamics (kynamics)" author: "SOCR Team (Ivo Dinov)" date: "`r format(Sys.time(),'%m/%d/%Y')`" output: html_document: theme: spacelab highlight: tango includes: before_body: TCIU_header.html toc: true number_sections: true toc_depth: 2 toc_float: collapsed: false smooth_scroll: true code_folding: hide --- ```{r setup, include=FALSE, echo=FALSE} knitr::opts_chunk$set(echo = TRUE, warings = FALSE, error = TRUE) ``` # Foundations of Dynamical System in Spacekime ## Kime-based dynamical system A spacekime dynamical (*kynamical*) system describes the evolution of a process or a phenomenon over a specific trajectory (path curve) in the kime plane. Describing a kynamical system requires (1) a specific simple (non-intersecting) parametric curve (could be closed) $\gamma :\mathbb{R} → \mathbb{C}$; (2) formulation of the process that is evolving over $\gamma$; and (3) a rule indicating how the process evolves over kime $κ=\gamma(s)=(κ_1 (s),κ_2 (s))\in \mathbb{C}$, where the curve parameter $s\in \mathbb{R}$. The resulting kynamical system represents a mathematical model describing the dynamic kime evolution of the system. Let’s examine each of the three components of a kynamical system. ## The kime trajectory The kime curve describes the path along which we are tracking, exploring, following, or modeling the dynamical changes of the process of interest. Planar curves in $\mathbb{C}$ can be represented in different coordinates, including Cartesian and Polar coordinates. For instance, implicit paths may be defined by *implicit* analytical equations such as $f(x,y)=0$, for some specific function $f:\mathbb{C}→\mathbb{R}$. Most of the time, implicit equations are not separable and do not permit explicit solutions for one of the variables. Occasionally, such equations may allow explicit solutions such as $y=g(x)$, for some specific function $g$. An alternative *parametric* representation of a plane curve in Cartesian coordinates uses $\gamma(s)=(κ_1 (s),κ_2 (s)): \mathbb{R}→\mathbb{R}^2$ for some specific functions $κ_1,κ_2: \mathbb{R}→\mathbb{R}$. Polar coordinate representation of plane curves is another alternative that expresses the location of each point in terms of a phase-angle and a radial-distance from the origin, $\gamma(A,ϕ)=Ae^{iϕ},A∈\mathbb{R}^+,ϕ∈[-π,π)$. ## The state space Kynamical systems also require the special formulation of the process that evolves over kime trajectories. This formulation is based on a set of variables that give a complete description of the system at any particular kime point in C. These descriptor variables describe the mathematical state of the process system. A kynamical system evaluated at a specific numerical vector representing an instance observation of these variables at a particular kime, would provide a complete picture of the state of the entire system at that kime. To model a specific kynamical system, the investigator has to identify the number, type, and interactions between all variables to form the complete description for the mathematical kynamics model. These variables that completely describe the state of the kynamical system are called the state variables, and the corresponding space covering the set of all the possible values of the state variables is the state space. The state space can be discrete, consisting of isolated points, e.g., when a state variable could only take on integer values, continuous, consisting of a contiguous set of points, e.g., when a state variable could take on any real value, or mixed, when some state variables are continuous and some are discrete. The special case of continuous-and-finite-dimensional state spaces are called the phase spaces, and the finite number of state variables represents the dimension of the kynamical system. Note that state spaces can also be infinite-dimensional. ## The kime evolution rule The last component of a kynamical system requires specific rules for the kime evolution of the process. The rule must be defined to ensure that the state variables represent a holistic description of the complete state space. In other words, the state variable vector at a particular kime must completely determine the system evolution to all future states. For instance, if the kime evolution depends on some extraneous variable not included in the state space, then the rule combined with the state space would not represent the complete process kynamics. This situation would necessitate a change in the rule, reformulation of the state space, or augmentation of the necessary state variables describing the kynamical system. Typically, kime evolution rules involve either *discrete* or *continuous* kime increments (not mixed discrete-and-continuous changes). In the *case of discrete* kime increments and using parametric kime trajectory representation, then the system evolves in kime steps that can be indexed by the natural numbers $s=0,1,2,\cdots \in \mathbb{N}$, and the state of the system at kime $κ(s)=(κ_1 (s),κ_2 (s))$ can be expressed as $X_κ$. In other cases, the kime evolution rule may be based on a function $κ(s+1)=f(κ(s))$, i.e., a function transforming the state of the system at one kime, $κ(s)$, into the state of the system at the next kime, $κ(s+1)$. Initial constraints or boundary conditions, such as $κ_o=κ(s=0)$ provide a recipe for determining the state at subsequent kime points, i.e., $κ_1=κ(s=1)=f(κ(0))=f(κ_o)$, $κ_2=κ(s=2)=f(κ(1))=f(f(κ_o))$, and so forth for all future states. In such discrete kynamical systems, the states along the curve $\gamma(s)=(κ_1 (s),κ_2 (s))$ are determined by the state-transition function $f(∙)$ and the initial state $κ_o=κ(s=0)$. In *continuous kynamical systems*, states evolve continuously, as a smooth flow, over the kime trajectory path $\gamma(s)=(κ_1 (s),κ_2 (s))$. As the curve parameter $s\in\mathbb{R}$ increases (think of $s$ as the curve-length parameter), we traverse the simple curve in kime and the state $X_κ=X_{κ(s)}$ evolution at kime $κ=k(s)$ can be thought of as a point navigating through the state space. In continuous kynamical systems, the evolution rule needs to specify the motion of the point $X_κ$ along $\gamma$ in terms of the motion velocity, e.g., $v(κ)=F(X_κ)$ at $κ=k(s)$. Assuming that the initial starting state at kime $κ_o=κ(s=0)$ is $X_o$ and we have specified a kime-trajectory $κ=\gamma(s)=(κ_1 (s),κ_2 (s))\in\mathbb{C}$, then the kynamics along $\gamma$ would represent a path through state space, $X_{κ(s)}$ starting with $X_o$. # Examples of kynamical systems ## Simulation example 1 Let’s define $$κ(s)=\left (κ_1 (s)=\frac{1}{2} s^2+1,\ \ κ_2 (s)=2+6(s-\frac{1}{3})(s-\frac{1}{2})(s-1) \right )$$ and the state $$X_κ^{(1)}= \begin{pmatrix} x_1 (κ) \\ x_2 (κ) \\ x_3 (κ) \end{pmatrix} = \begin{pmatrix} κ_2^2 - κ_1 κ_2^2 \\ κ_1 - κ_2 \\ κ_1 + κ_2 \end{pmatrix},$$ where $0≤s≤1, \ 1<κ_1<3/2, \ 1<κ_2<2.1$. The graphs below show a kime plane curve $\gamma(s)=(κ_1 (s),κ_2 (s))$ and the kime-surface $X_κ^{(1)}$. This 1D Parameterization of the kime-dynamics (kynamics) simulation example illustrates a parametric complex-plane curve and a corresponding kime-surface defined in terms of kime. Note that visualization of kime-surfaces is generally difficult and can be confused with the graphs of functions that have 2D input (domain) and a 1D output (range), which can be drawn as surfaces embedded in $\mathbb{R}^3$. In reality, parameterized kime-functions have a different flavor as they have a 2D input $κ\in \mathbb{C}$ and a 3D output $X_κ^{(1)}$, which implies that plotting them requires 5D coordinates. Note that the thick blue tube curling through the kime-surface represents the kynamics along $\gamma$ as a path through the state space, $X_κ^{(1)}$. The three thinner curves embedded in the three cardinal projection planes illustrate the silhouettes of the dynamical system. ```{r eval=T, message=FALSE, warning=FALSE, error=FALSE} library(plotly) # Define curve parameter s <- seq(0, 1, length.out=400) # define kime functions # k1 <- 2*s^2 + 1 # k2 <- 2 + 12*(s-1/3)*(s-1/2)*(s-1) k1 <- (1/2)*s^2 + 1 k2 <- 2 + 6*(s-1/3)*(s-1/2)*(s-1) # Define state space X_{\kappa} x <- k1^2 - k1*k2^2 y <- k1 - k2 z <- k1 + k2 data1 <- data.frame(x, y, z) # define axes labels f <- list(family = "Courier New, monospace", size = 18, color = "black") xlab <- list(title = "k1", titlefont = f) ylab <- list(title = "k2", titlefont = f) zlab <- list(title = "Kime-surface", titlefont = f) # define kime-surface xsurf <- k1 ysurf <- k2 z_fun <- function(x,y) { x + y } zsurf <- outer(xsurf, ysurf, FUN="z_fun") # 2D Gamma Kime-space plot # k1 <- (1/2)*s^2 + 1 # k2 <- 2 + 6*(s-1/3)*(s-1/2)*(s-1) # Define state space X_{\kappa} # x <- k1^2 - k1*k2^2 # y <- k1 - k2 plot_ly(x = ~k1, y = ~k2, type = 'scatter', mode = 'lines', line = list(width = 4)) %>% layout(xaxis = list(title = "k1", range = c(0.95, 1.55)), yaxis = list(title = "k2", range = c(0.95, 2.05))) # 3D Kime-surface Plot # plot_ly(x = ~k1, y = ~k2, z = ~t(zsurf), type="surface", opacity=0.7) p <- plot_ly(data1, x = ~x, y = ~y, z = ~z, type = 'scatter3d', mode = 'lines', # name="Kime-space curve projection on Kime-surface", line = list(width = 22, color = ~c, colorscale = list(c(0,'red'), c(1,'blue')))) %>% # add kime-surface add_trace(z = ~zsurf, type = 'surface', opacity=0.4) %>% # trace the z-axis add_trace(data1, x = 0, y =0, z = ~z, type="scatter3d", mode="lines", line = list(width = 10, color = ~c , colorscale = list(c(0,'red'), c(1,'blue'))), name="Z", hoverinfo="none") %>% # add projection on plane x=0.2 add_trace(data1, x = 0.2, y = ~y, z = ~z, type="scatter3d", mode="lines", line = list(width = 2, dash="solid", color = ~c, colorscale = list(c(0,'gray'), c(1,'black'))), name="Z", hoverinfo="none") %>% # add projection on plane y=0.2 add_trace(data1, x = ~x, y = 0.2, z = ~z, type="scatter3d", mode="lines", line = list(width = 2, dash="solid", color = ~c, colorscale = list(c(0,'gray'), c(1,'black'))), name="Z", hoverinfo="none") %>% # add projection on plane z=1.8 add_trace(data1, x = ~x, y = ~y, z = 1.8, type="scatter3d", mode="lines", line = list(width = 10, dash="solid", color = ~c, colorscale = list(c(0,'gray'), c(1,'black'))), name="Z", hoverinfo="none") %>% layout(title = "Kynamical System Simulation Example", showlegend = FALSE, scene = list(aspectmode = "manual", xaxis = xlab, yaxis = ylab, zaxis = zlab, aspectratio = list(x=1, y=1, z=1))) # 1:1:1 aspect ratio p ``` ## Simulation example 2 Let’s define $κ=(κ_1,κ_2)$ and the state $$X_κ^{(2)}=\begin{pmatrix} x_1 (κ) \\ x_2 (κ) \\ x_3 (κ) \end{pmatrix} = \begin{pmatrix} κ_2^3 - κ_1 κ_2 \\ κ_1-κ_2 \\ κ_1+κ_2 \end{pmatrix},$$ where $-2≤κ_1≤2; -1≤κ_2≤1$. Note that the pair of arguments in this example are independent. - *Input*: 2D kime-space parameterization $\gamma (κ)=(κ_1,κ_2 )$ of a simple closed curve; rectangular parametric boundary $∂D$, where $D=\{-2≤κ_1≤2 ∩ -1≤κ_2≤1\}$. - *Output*: $X_κ^{(2)}=(κ_2^3-κ_1 κ_2,\ κ_1+κ_2,\ κ_1-κ_2 )'$. Left 3D display shows a checkerboard surface texture-mapping providing cues of surface curvature. The right 3D rendering depicts the embedding within the 3D manifold of the curvy-linear representation of the initial kime-region boundary as a simple-closed curve representation $\partial D\subset \mathbb{C}$. The graphs below show a generic 2D Parameterization of the kime-dynamics over a two-parameter complex-plane curve representing the boundary $∂D$ of a rectangle $D⊂\mathbb{C}$ and a corresponding kime-surface. These show the corresponding 3D surface kime-parameterized by $(κ_1,κ_2 )$ and a checkerboard shading as a proxy of the surface geodesic metric. The second plot depicts the corresponding transformed boundary embedded in the manifold. Just like in the earlier single-parameter case, the thick blue tube curling through the kime-manifold represents the kynamics along $\gamma=\gamma(κ_1,κ_2 )$, as a path through the state space, $X_κ^{(2)}$. The three thinner curves embedded in the three cardinal projection planes illustrate the silhouettes of the dynamical system. ```{r eval=T, message=FALSE, warning=FALSE, error=FALSE} # Kime domain rectangle kime_grid_length <- 100 kime_k1_range <- seq(-2,2, length=kime_grid_length) kime_k2_range <- seq(-1,1,length=kime_grid_length) kime_grid <- expand.grid(kime_k1_range,kime_k2_range) # dim(kime_grid) # [1] 900 2 # head(kime_grid) # Var1 Var2 # 1 -2.000000 -1 # 2 -1.862069 -1 # 3 -1.724138 -1 # 4 -1.586207 -1 # 5 -1.448276 -1 # 6 -1.310345 -1 # Plot Input Grid plot_ly(type="scatter", mode = 'lines', xaxis ="k1", yaxis ="k2") %>% layout(xaxis = list(gridcolor = "black", range = c(-3.1, 3.1)), yaxis = list(gridcolor = "black", range = c(-2.1, 2.1)), shapes = list( list(type = "rect", fillcolor = "lightblue", line = list(width = 4, color = "blue"), opacity = 0.3, x0 = -2, x1 = 2, xref = "k1", y0 = -1, y1 = 1, yref = "k2")) ) # ax <- list(zeroline = TRUE, showline = TRUE, mirror = "ticks", gridcolor = toRGB("gray50"), # gridwidth = 2, zerolinecolor = toRGB("red"), zerolinewidth = 4, # linecolor = toRGB("black"),linewidth = 2) # s <- seq(-1, 4, length=10) # plot_ly() %>% # layout(xaxis = ax, yaxis = ax) # Define 3D Manifold x <- kime_grid$Var1^3 - kime_grid$Var1*kime_grid$Var2 y <- kime_grid$Var1 - kime_grid$Var2 z <- kime_grid$Var1 + kime_grid$Var2 # Define manifold colors # Create an empty array of strings with the same shape as the meshgrid # storing alternating two colors in a checkerboard pattern. colorPair <- c(rep('white', 5), rep('black',5)) #colorPair <- c('white','black') surf_colors <- array("", dim = c(kime_grid_length,kime_grid_length)) # Check checkerboard pattern for (j in 1:kime_grid_length) { for (i in 1:kime_grid_length) { surf_colors[j,i] <- colorPair[1+ (i+j) %% length(colorPair)] # surf_colors[j,i] <- colorPair[1+ ((-1+i+(j)/2) %% 10)] } } # Check checkerboard pattern # image(matrix(unclass(as.factor(surf_colors)), nrow=kime_grid_length)) plot_ly(z= data.matrix(data.frame(unclass(surf_colors))) , type="heatmap") %>% hide_colorbar() # %>% hide_colorbar() kime_plot <- plot_ly() kime_plot <- add_trace(kime_plot, x=x, y=y, z=z, type="mesh3d", name="3D Surface (kime)", showlegend=FALSE, showscale = FALSE, intensity = unclass(as.factor(surf_colors)), #y colors = colorRamp(rainbow(8)), opacity=0.95) %>% # colors = colorRamp(surf_colors), opacity=0.9) %>% layout(scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1, z=1))) kime_plot dim(matrix(unclass(as.factor(surf_colors)), nrow=kime_grid_length)) s <- seq(0, 1, length.out=400) # define kime BOUNDARY (\partial D) CURVE (rectangular region D) k1 <- k2 <- s for (i in 1:4) { if (i==1) { # 0<=s<1/4 k1[1:100] <- 16*s[1:100] -2 # (-2,-1) --> (2,-1) k2[1:100] <- -1 } else if (i==2) { # 1/4<=s<1/2 k1[101:200] <- 2 k2[101:200] <- 8*s[101:200] - 3 # (2,-1) --> (2,1) } else if (i==3) { # 1/2<=s<3/4 k1[201:300] <- 10-16*s[201:300] # (2,1) --> (-2,1) k2[201:300] <- 1 } else { # if (i==4) # 3/4<=s<1 k1[301:400] <- -2 k2[301:400] <- 7 - 8*s[301:400] # (-2,1) --> (-2,-1) } } # Plot Input kime ROI boundary Grid plot_ly(x=~k1, y=~k2, type="scatter", mode = 'lines', xaxis ="k1", yaxis ="k2", line = list(width = 10, color = "black")) %>% layout(xaxis = list(gridcolor = "black", range = c(-3.1, 3.1)), yaxis = list(gridcolor = "black", range = c(-2.1, 2.1)), shapes = list( list(type = "rect", fillcolor = "lightblue", line = list(width = 4, color = "blue"), opacity = 0.3, x0 = -2, x1 = 2, xref = "k1", y0 = -1, y1 = 1, yref = "k2")) ) # Define state space X_{\kappa} x <- k1^3 - k1*k2 y <- k1 - k2 z <- k1 + k2 data1 <- data.frame(x, y, z) # define axes labels f <- list(family = "Courier New, monospace", size = 18, color = "black") xlab <- list(title = "k1", titlefont = f) ylab <- list(title = "k2", titlefont = f) zlab <- list(title = "Kime-surface", titlefont = f) # define kime-surface xsurf <- k1 ysurf <- k2 z_fun <- function(x,y) { x + y } zsurf <- outer(xsurf, ysurf, FUN="z_fun") # 3D Kime-surface Plot p <- plot_ly(x = ~x, y = ~y, z = ~z, type = 'scatter3d', mode = 'lines', # name="Kime-space curve projection on Kime-surface", line = list(width = 22, color = ~c, colorscale = list(c(0,'red'), c(1,'blue')))) %>% # add tube over boundary of kime-region add_trace(z = ~zsurf, type = 'surface', opacity=0.4) %>% # # trace the z-axis # add_trace(x = 0, y =0, z = ~z, type="scatter3d", mode="lines", # line = list(width = 10, color = ~c , colorscale = list(c(0,'red'), c(1,'blue'))), # name="Z", hoverinfo="none") %>% # add projection on plane x=10.1 add_trace(data1, x = 10.1, y = ~y, z = ~z, type="scatter3d", mode="lines", line = list(width = 2, dash="solid", color = ~c, colorscale = list(c(0,'gray'), c(1,'black'))), name="Z", hoverinfo="none") %>% # add projection on plane y=3.1 add_trace(data1, x = ~x, y = 3.1, z = ~z, type="scatter3d", mode="lines", line = list(width = 2, dash="solid", color = ~c, colorscale = list(c(0,'gray'), c(1,'black'))), name="Z", hoverinfo="none") %>% # add projection on plane z=3.1 add_trace(data1, x = ~x, y = ~y, z = 3.1, type="scatter3d", mode="lines", line = list(width = 10, dash="solid", color = ~c, colorscale = list(c(0,'gray'), c(1,'black'))), name="Z", hoverinfo="none") %>% layout(title = "Kynamical System Simulation Example", showlegend = FALSE, scene = list(aspectmode = "manual", xaxis = xlab, yaxis = ylab, zaxis = zlab, aspectratio = list(x=1, y=1, z=1))) # 1:1:1 aspect ratio p ``` ## Torus example 3 Let’s define $κ=(κ_1,κ_2 )$, $-π≤κ_1<π$, $-π≤κ_2<π$ and the state $$X_κ^{(3)}=\begin{pmatrix} x_1 (κ) \\ x_2 (κ) \\ x_3 (κ) \end{pmatrix} = \begin{pmatrix} (6+2 \cos(κ_2))\cos(κ_1) \\ (6+2 \cos(κ_2))\sin(κ_1) \\ \sin(κ_2) \end{pmatrix}.$$ - *Input*: 2D kime-space parameterization $\gamma (κ)=(κ_1,κ_2 )$, $κ_1=s-π$, $κ_2=(1/10)((κ_1-1)^3+5(κ_1-1)^2+2(κ_1-5)- 8)$ representing a simple kime-curve. - *Output*: $X_κ^{(3)}$. 3D display shows the torus surface rendered in 3D along with the curvy-linear representation of the initial kime-curve representation. The figures below show the 2D Torus Parametrization and the kime-dynamics over a two-parameter complex-plane curve and the corresponding 2D kime-surface. The second graph illustrates the corresponding 3D torus surface parameterized by $(κ_1,κ_2 )$ the corresponding transformed kime-curve embedded in the torus manifold. The thick red tube curling through the torus surface represents the kynamics along $\gamma=\gamma(κ_1,κ_2 )$, as a path through the state space, $X_κ^{(3)}$. The three thinner curves embedded in the three cardinal projection planes illustrate the silhouettes of the dynamical system. ```{r eval=T, message=FALSE, warning=FALSE, error=FALSE} # Kime domain rectangle library(plotly) kime_grid_length <- 100 kime_k1_range <- seq(-pi, pi, length=kime_grid_length) kime_k2_range <- seq(-pi, pi,length=kime_grid_length) kime_grid <- expand.grid(kime_k1_range,kime_k2_range) # Plot Input Grid plot_ly(type="scatter", mode = 'lines', xaxis ="k1", yaxis ="k2") %>% layout(xaxis = list(gridcolor = "black", range = c(-4.1, 4.1), dtick = 0.2), yaxis = list(gridcolor = "black", range = c(-4.1, 4.1), dtick = 0.2), #shapes = lapply(x0=-3, x1=3, y=(-9:9)/5, hline), shapes = list( list(type = "rect", fillcolor = "lightblue", line = list(width = 4, color = "blue"), opacity = 0.3, x0 = -pi, x1 = pi, xref = "k1", y0 = -pi, y1 = pi, yref = "k2")) ) # Define 3D Manifold x <- 5 + 2*cos(kime_grid$Var2)*cos(kime_grid$Var1) y <- 5 + 2*cos(kime_grid$Var2)*sin(kime_grid$Var1) z <- 2*sin(kime_grid$Var2) # Define manifold colors # Create an empty array of strings with the same shape as the meshgrid # storing alternating two colors in a checkerboard pattern. #colorPair <- c(rep('white', 5), rep('black',5)) colorPair <- c('white','black') surf_colors <- array("", dim = c(kime_grid_length,kime_grid_length)) # Check checkerboard pattern for (j in 1:kime_grid_length) { for (i in 1:kime_grid_length) { surf_colors[j,i] <- colorPair[1+ (i+j) %% length(colorPair)] # surf_colors[j,i] <- colorPair[1+ ((-1+i+(j)/2) %% 10)] } } # Check checkerboard pattern # image(data.matrix(data.frame(unclass(surf_colors)))) plot_ly(z= data.matrix(data.frame(unclass(surf_colors))) , type="heatmap") %>% hide_colorbar() kime_plot <- plot_ly() kime_plot <- add_trace(kime_plot, x=x, y=y, z=z, type="mesh3d", name="3D Surface (kime)", showlegend=FALSE, showscale = FALSE, intensity = unclass(as.factor(surf_colors)), #y colors = colorRamp(rainbow(8)), opacity=0.95) %>% # colors = colorRamp(surf_colors), opacity=0.9) %>% layout(scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1, z=1))) kime_plot dim(matrix(unclass(as.factor(surf_colors)), nrow=kime_grid_length)) # s <- seq(0, 1, length.out=400) # # # define kime BOUNDARY (\partial D) CURVE (rectangular region D) # k1 <- k2 <- s # for (i in 1:4) { # if (i==1) { # 0<=s<1/4 # k1[1:100] <- 16*s[1:100] -2 # (-2,-1) --> (2,-1) # k2[1:100] <- -1 # } else if (i==2) { # 1/4<=s<1/2 # k1[101:200] <- 2 # k2[101:200] <- 8*s[101:200] - 3 # (2,-1) --> (2,1) # } else if (i==3) { # 1/2<=s<3/4 # k1[201:300] <- 10-16*s[201:300] # (2,1) --> (-2,1) # k2[201:300] <- 1 # } else { # if (i==4) # 3/4<=s<1 # k1[301:400] <- -2 # k2[301:400] <- 7 - 8*s[301:400] # (-2,1) --> (-2,-1) # } # } # # # Plot Input kime ROI boundary Grid # plot_ly(x=~k1, y=~k2, type="scatter", mode = 'lines', xaxis ="k1", yaxis ="k2", # line = list(width = 10, color = "black")) %>% # layout(xaxis = list(gridcolor = "black", range = c(-3.1, 3.1)), # yaxis = list(gridcolor = "black", range = c(-2.1, 2.1)), # shapes = list( # list(type = "rect", # fillcolor = "lightblue", line = list(width = 4, color = "blue"), opacity = 0.3, # x0 = -2, x1 = 2, xref = "k1", # y0 = -1, y1 = 1, yref = "k2")) # ) # # # Define state space X_{\kappa} # x <- k1^3 - k1*k2 # y <- k1 - k2 # z <- k1 + k2 # # data1 <- data.frame(x, y, z) # # # define axes labels # f <- list(family = "Courier New, monospace", size = 18, color = "black") # xlab <- list(title = "k1", titlefont = f) # ylab <- list(title = "k2", titlefont = f) # zlab <- list(title = "Kime-surface", titlefont = f) # # # define kime-surface # xsurf <- k1 # ysurf <- k2 # z_fun <- function(x,y) { x + y } # zsurf <- outer(xsurf, ysurf, FUN="z_fun") # # # 3D Kime-surface Plot # p <- plot_ly(x = ~x, y = ~y, z = ~z, type = 'scatter3d', # mode = 'lines', showlegend = F, # name="Kime-space curve projection on Kime-surface", # line = list(width = 22, color = ~c, colorscale = list(c(0,'red'), c(1,'blue')))) %>% # # add tube over boundary of kime-region # add_trace(z = ~zsurf, type = 'surface', opacity=0.4) %>% # # # trace the z-axis # # add_trace(x = 0, y =0, z = ~z, type="scatter3d", mode="lines", # # line = list(width = 10, color = ~c , colorscale = list(c(0,'red'), c(1,'blue'))), # # name="Z", hoverinfo="none") %>% # # add projection on plane x=10.1 # add_trace(data1, x = 10.1, y = ~y, z = ~z, type="scatter3d", mode="lines", # line = list(width = 2, dash="solid", color = ~c, # colorscale = list(c(0,'gray'), c(1,'black'))), # name="Z", hoverinfo="none") %>% # # add projection on plane y=3.1 # add_trace(data1, x = ~x, y = 3.1, z = ~z, type="scatter3d", mode="lines", # line = list(width = 2, dash="solid", color = ~c, # colorscale = list(c(0,'gray'), c(1,'black'))), # name="Z", hoverinfo="none") %>% # # add projection on plane z=3.1 # add_trace(data1, x = ~x, y = ~y, z = 3.1, type="scatter3d", mode="lines", # line = list(width = 10, dash="solid", color = ~c, # colorscale = list(c(0,'gray'), c(1,'black'))), # name="Z", hoverinfo="none") %>% # layout(title = "Kynamical System Simulation Example", showlegend = FALSE, # scene = list(aspectmode = "manual", # xaxis = xlab, yaxis = ylab, zaxis = zlab, # aspectratio = list(x=1, y=1, z=1))) # 1:1:1 aspect ratio # p ############ Torus a <- 6 # Torus radius r <- 2 # Tube radius ###### Delaunay Triangulation library(geometry) # Parameter defining how many points to take. grid.df <- expand.grid( u = seq(0, 2*pi, length.out = 20), v = seq(0, 2*pi, length.out = 20) ) # Define a Torus parametrization (see Wikipedia) mat <- matrix( c(x<- (a + r*cos(grid.df$u))*cos(grid.df$v), # x y <- (a + r*cos(grid.df$u))*sin(grid.df$v), # y z <- r*sin(grid.df$u) # z ), ncol = 3, dimnames = list(NULL, c("x", "y", "z")) ) Delaunay.mat <- delaunayn(grid.df) #, options="Qcc Qc Qt Qz", output.options=TRUE) # Delaunay.mat.t <- t(Delaunay.mat) # Use mat to plot (not the 2D grid grid.df) # Plotly layout axs <- list( backgroundcolor="rgb(200,200,200)", # gray gridcolor="rgb(255,255,255)", # white showbackground=TRUE, zerolinecolor="rgb(255,255,255)" # white ) # Apply the colormap # Compute the mean of z for each row of the Delaunay vertices zmean <- apply(Delaunay.mat, MARGIN=1, function(row){mean(mat[row,3])}) library(scales) # Determine the 2-cell face's colors # plotted color result will be slightly different # since colour_ramp uses CIELAB instead of RGB # could use colorRamp for exact replication facecolor = colour_ramp( # brewer_pal(palette="RdBu")(10); # brewer_pal("div")(10) colorRampPalette(c("pink", "purple"))(10)) # (rescale(x=zmean)) scene = list(camera = list(eye = list(x = 2, y = 1, z = 1)), aspectmode = "manual", aspectratio = list(x=2, y=2, z=1)) plot_ly( # x = x, y = y, z = z, # vertex (0-cell) coordinates # i = i, j = j, k = k, # indices to the vertices, which describe a 2-cell (face or a simplex) of dimension dim # see docs: https://plot.ly/r/3d-mesh-plots/ x = mat[, 1], y = mat[, 2], z = mat[, 3], # JavaScript is 0 based index so subtract 1 i=Delaunay.mat[, 1]-1, j=Delaunay.mat[, 2]-1, k=Delaunay.mat[, 3]-1, # facecolor = facecolor, # mat[, 3], type = "mesh3d", opacity = 0.5, contour=list(show=TRUE, color="#000", width=15) ) %>% # add_trace(x=mat[, 1], y=mat[, 2], z=0, type="scatter3d", mode="lines", line = list(color="rgb(0, 255, 0)")) %>% layout(title="Torus Triangulation", scene=scene) #### Parametric curve on torus k1 <- grid.df$u -pi k2 <- ((k1-1)^3 + 5*(k1-1)^2 + 2*(k1-5) - 8)/10 matCurve <- matrix( c(xCurve<- grid.df$u -pi, # x yCurve <- ((xCurve-1)^3 + 5*(xCurve-1)^2 + 2*(xCurve-5) - 8)/10, # y zCurve <- r*sin(grid.df$u) # z ), ncol = 3 ) # 3D Torus with Curve-Tube Plot p <- plot_ly(x = ~x, y = ~y, z = ~z, type = 'scatter3d', mode = 'lines', showlegend = F, # name="Kime-space curve projection on Kime-surface", line = list(width = 22, color = ~c, colorscale = list(c(0,'red'), c(1,'blue')))) %>% # add tube over boundary of kime-region add_trace(z = ~zsurf, type = 'surface', opacity=0.4) # %>% # # trace the z-axis # add_trace(x = 0, y =0, z = ~z, type="scatter3d", mode="lines", # line = list(width = 10, color = ~c , colorscale = list(c(0,'red'), c(1,'blue'))), # name="Z", hoverinfo="none") %>% plot_ly( # x = x, y = y, z = z, # vertex (0-cell) coordinates # i = i, j = j, k = k, # indices to the vertices, which describe a 2-cell (face or a simplex) of dimension dim # see docs: https://plot.ly/r/3d-mesh-plots/ x = mat[, 1], y = mat[, 2], z = mat[, 3], # JavaScript is 0 based index so subtract 1 i=Delaunay.mat[, 1]-1, j=Delaunay.mat[, 2]-1, k=Delaunay.mat[, 3]-1, facecolor = facecolor, type = "mesh3d", opacity = 0.5, contour=list(show=TRUE, color="#000", width=15)) %>% add_trace(x=matCurve[, 1], y=matCurve[, 2], z=matCurve, type="scatter3d", mode="lines", line = list(color="rgb(0, 255, 0)")) %>% layout(title="Torus Triangulation", scene=scene) kime_grid_length <- 20 kime_k1_range <- seq(-pi,pi, length=kime_grid_length) kime_k2_range <- seq(-pi,pi, length=kime_grid_length) kime_grid <- expand.grid(kime_k1_range,kime_k2_range) # Plot 2D Input Grid + kime curve plot_ly(type="scatter", mode = 'lines', xaxis ="k1", yaxis ="k2") %>% add_lines(x=matCurve[, 1], y=matCurve[, 2], mode="lines", line = list(color="rgb(255, 0, 0)", width=5)) %>% layout(showlegend=FALSE, xaxis = list(gridcolor = "black", range = c(-4, 4), dtick = 0.2), yaxis = list(gridcolor = "black", range = c(-4, 4), dtick = 0.2), shapes = list( list(type = "rect", fillcolor = "lightblue", line = list(width = 4, color = "blue"), opacity = 0.3, x0 = -pi, x1 = pi, xref = "k1", y0 = -pi, y1 = pi, yref = "k2")) ) # Checkerboard Torus + kime-curve x<- (a + r*cos(grid.df$u))*cos(grid.df$v) # x y <- (a + r*cos(grid.df$u))*sin(grid.df$v) # y z <- r*sin(grid.df$u) # z # Define manifold colors # Create an empty array of strings with the same shape as the meshgrid # storing alternating two colors in a checkerboard pattern. #colorPair <- c(rep('white', 5), rep('black',5)) colorPair <- c('white','black') surf_colors <- array("", dim = c(length(grid.df$u), length(grid.df$v))) # Check checkerboard pattern for (j in 1:length(grid.df$u)) { for (i in 1:length(grid.df$v)) { surf_colors[j,i] <- colorPair[1+ (i+j) %% length(colorPair)] # surf_colors[j,i] <- colorPair[1+ ((-1+i+(j)/2) %% 10)] } } # Check checkerboard pattern # image(matrix(unclass(as.factor(surf_colors)), nrow=kime_grid_length)) plot_ly(z= data.matrix(data.frame(unclass(surf_colors))) , type="heatmap") %>% hide_colorbar() kime_plot <- plot_ly() kime_plot <- add_trace(kime_plot, x=x, y=y, z=z, type="mesh3d", name="3D Surface (kime)", showlegend=FALSE, intensity = unclass(as.factor(surf_colors)), #y colors = colorRamp(rainbow(8)), opacity=0.95) %>% # colors = colorRamp(surf_colors), opacity=0.9) %>% layout(scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1, z=1))) kime_plot ###################################### TEST # shape=="complex") # phi <- seq(from = 0, to = 2*pi, by = ((2*pi - 0)/(200 - 1))) # psi <- seq(from = 0, to = 2*pi, by = ((2*pi - 0)/(200 - 1))) # # rendering (u,v) parametric surfaces requires x,y,z arguments to be 2D arrays # # In out case, the three coordinates have to be 200*200 parameterized tensors/arrays # r1 = 2 + sin(3 * phi + 5 * psi) # r = 2 + sin(7phi+5psi) # x1 = (r1 * cos(phi)) %o% sin(psi) # x = r*cos(phi)*sin(psi) # y1 = (r1 * sin(phi)) %o% sin(psi) # y = r*sin(phi)*sin(psi) # z1 = r1 %o% cos(psi) # z = r*cos(psi) # # p <- plot_ly(hoverinfo="none", showscale = FALSE) %>% # add_trace(x = ~x1, y = ~y1, z = ~z1, type = 'surface', opacity=1, visible=T) %>% # layout(title = "Torus", showlegend = FALSE) # p # Checkerboard TORUS phi <- grid.df$u # seq(from = 0, to = 2*pi, by = ((2*pi - 0)/(200 - 1))) psi <- grid.df$v # seq(from = 0, to = 2*pi, by = ((2*pi - 0)/(200 - 1))) xOuter <- outer(phi, psi, function(phi, psi) { cos(phi)*(a+r*cos(psi)) }) yOuter <- outer(phi, psi, function(phi, psi) { sin(phi)*(a+r*cos(psi)) }) zOuter <- outer(phi, psi, function(phi, psi) { r*sin(psi) }) plot_ly(hoverinfo="none", showscale = FALSE) %>% add_trace(x=~xOuter, y=~yOuter, z=~zOuter, type='surface', opacity=1.0, visible=T) %>% layout(title = "Torus", showlegend = FALSE, scene = list(aspectmode = "manual", aspectratio = list(x=2, y=2, z=1))) ################ checkerboard color pattern a <- 6 # Torus radius (from torus gravitational center) r <- 2 # Tube radius (from torus circular core) # phi <- seq(from = 0, to = 2*pi, by = ((2*pi - 0)/(20 - 1))) # psi <- seq(from = 0, to = 2*pi, by = ((2*pi - 0)/(20 - 1))) xOuter <- outer(phi, psi, function(phi, psi) { cos(phi)*(a+r*cos(psi)) }) yOuter <- outer(phi, psi, function(phi, psi) { sin(phi)*(a+r*cos(psi)) }) zOuter <- outer(phi, psi, function(phi, psi) { r*sin(psi) }) # checkerboard color pattern colorPair <- c('white','black') surf_colors <- array("", dim = c(length(phi), length(phi))) # Check checkerboard pattern for (j in 1:length(phi)) { for (i in 1:length(phi)) { surf_colors[j,i] <- colorPair[1+ (i+j) %% length(colorPair)] # surf_colors[j,i] <- colorPair[1+ ((-1+i+(j)/2) %% 10)] } } # Check checkerboard pattern # image(matrix(unclass(as.factor(surf_colors)), nrow=length(phi))) plot_ly(z= data.matrix(data.frame(unclass(surf_colors))) , type="heatmap") %>% hide_colorbar() ### Kime tube on torus #### Parametric curve on torus # k1 <- phi -pi # k2 <- ((phi-1)^3 + 5*(phi-1)^2 + 2*(phi-5) - 8)/10 xCurveOuter <- outer(phi, psi, function(phi, psi) { cos(phi -pi)*(a+r*cos(((phi -pi-1)^3 + 5*(phi -pi-1)^2 + 2*(phi -pi-5) - 8)/10)) }) yCurveOuter <- outer(phi, psi, function(phi, psi) { sin(phi -pi)*(a+r*cos(((phi -pi-1)^3 + 5*(phi -pi-1)^2 + 2*(phi -pi-5) - 8)/10)) }) zCurveOuter <- outer(phi, psi, function(phi, psi) { r*sin(((phi -pi-1)^3 + 5*(phi -pi-1)^2 + 2*(phi -pi-5) - 8)/10) }) plot_ly(x=~as.vector(xCurveOuter), y=~as.vector(yCurveOuter), z=~as.vector(zCurveOuter), type="scatter3d", mode="markers+lines", line=list(color='red', width=15)) plot_ly(hoverinfo="none") %>% add_trace(x=~xOuter, y=~yOuter, z=~zOuter, type='surface', opacity=0.01, visible=TRUE, showlegend= FALSE, showscale = FALSE) %>% add_trace(x=~as.vector(xCurveOuter), y=~as.vector(yCurveOuter), z=~as.vector(zCurveOuter), name="Tube", type='scatter3d', mode='markers+lines', showscale = FALSE, hoverinfo = 'Tube', line=list(color='red', width=15), showlegend = FALSE) %>% # add projection on plane x=8.5 add_trace(x = 8.5, y = ~as.vector(yCurveOuter), z=~as.vector(zCurveOuter), type="scatter3d", mode="lines", line = list(width = 2, dash="solid", color = ~c, colorscale = list(c(0,'gray'), c(1,'black'))), name="Z", hoverinfo="none") %>% # add projection on plane y=-8.5 add_trace(x = ~as.vector(xCurveOuter), y = -8.5, z = ~as.vector(xCurveOuter), type="scatter3d", mode="lines", line = list(width = 2, dash="solid", color = ~c, colorscale = list(c(0,'gray'), c(1,'black'))), name="Z", hoverinfo="none") %>% # add projection on plane z=-2.5 add_trace(x = ~as.vector(xCurveOuter), y = ~as.vector(yCurveOuter), z = ~-2.5, type="scatter3d", mode="lines", line = list(width = 2, dash="solid", color = ~c, colorscale = list(c(0,'gray'), c(1,'black'))), name="Z", hoverinfo="none") %>% layout(title = "Torus", showlegend = FALSE, scene = list(aspectmode = "manual", aspectratio = list(x=2, y=2, z=1))) ``` ## fMRI example [See this TCIU Section (Kime-series/kime-surfaces reconstruction, visualization, and predictive analytics)](https://www.socr.umich.edu/TCIU/HTMLs/Chapter4_TCIU_Predictive_Analytics.html). ## Kynamics example of bacteria growth and decay We will demonstrate a kime dynamical system modeling the growth/decay of a specific population. Denote the bacteria population size at time $t$ by $f(t)$. As each bacterium divides into two bacteria, the entire population grows. Generally speaking not all bacteria would divide at the same time and some will die; hence the process is stochastic. For simplicity, assume that the kime radial growth of the bacteria colony is modeled by an exponential growth model, $f(t)=A_o a^t$, where $a > 1$ is the (growth) parameter and the constant $A_o$ represents the initial starting population at $t=0$. The overall kynamic growth or decay will be modeled by a bivariate joint distribution that depends on time $t$ and the random kime-phase $\phi \sim \Phi[-\pi, \pi)$, for some symmetric finitely supported kime-phase distribution $\Phi$. For a fixed time $t$, the bacteria population size is controlled by the marginal phase function $$g(\phi)=\left ( \frac{1}{2a}\right )^{\frac{\phi}{\pi}}:[-\pi, \pi) \to \left [ \frac{1}{2a}, 2a \right ].$$ Now that we have the individual marginal distributions for the kime-magnitude (time) and the kime-phase, we can formulate the joint distribution for bacteria growth with respect to complex-time (kime). ### Multivariate Joint Probability Distribution Copula Modeling In special cases where the components $\{X_i\}$ are independent, ${{\bf{X}}=\{X_1,X_2,\cdots,X_n\}}$, computing the joint multivariate distribution involves simply the multiplication of the individual marginal probability distributions, $f_{X_i}(x_i)$, i.e., $f_{{\bf{X}}}({\bf{x}})=\prod_{i=1}^n{f_{X_i}(x_i)}$. However, in general, there are intricate inter-dependencies between variables that need to be accounted for when computing the joint probability distribution $f_{{\bf{X}}}({\bf{x}})$ from the corresponding marginal univariate distribution models, $f_{X_i}(x_i)$. There are alternative strategies to express the joint probability distribution in terms of the marginals. A common approach requires *linking* the marginal distributions via connecting **copula function models**. Copula functions enable the formulation of multivariate probability distributions in terms of a set of given univariate marginal densities subject to some interdependence. For a pair of density functions, $f_X (x)$ and $f_Y (y)$, and their corresponding cumulative probability distribution functions, $F_X (x)$ and $F_Y (y)$, the joint bivariate cumulative distribution function model, $F_{XY} (x,y)$ can be expressed via a distribution copula $C(∙,∙)$ function. Similarly, the bivariate copula density function model, $f_{XY} (x,y)$ can be explicated in terms of some appropriate probability density copula function, $c(∙,∙)$: $$F_{XY} (x,y)=F_X (x)\times F_Y (y)\times C(F_X (x),F_Y (y)),$$ $$f_{XY} (x,y)=f_X (x)\times f_Y (y)\times c(f_X (x),f_Y (y)).$$ Specifying different copula function models yields certain types of joint probability distribution reconstructions. Some copula methods are only applicable to the continuous distributions; others are only limited to certain ranges. Let's consider three alternative copula models. * The [Farlie–Gumbel–Morgenstern (FGM) copula](https://en.wikipedia.org/wiki/Copula_(probability_theory)#List_of_copula_density_functions_and_applications). For a given constant $K$, the FGM model is commonly used for the bivariate normal cases: $$C_{XY}^{FGM} (x,y)= 1+K[1-2\times F_X(x)][1-2\times F_Y(y)].$$ * The [Exponential copula](https://en.wikipedia.org/wiki/Copula_(probability_theory)#List_of_copula_density_functions_and_applications) represents an alternative copula model type that is more stable with departures from normality of the marginal distributions: $$C_{XY}^{Exp} (x,y)=\frac{1}{1-ρ} \times \exp \left ( \frac{ρ[\ln⁡(1-F_X (x))+ \ln(1-F_Y(y))]}{1-ρ} \right ).$$ * The (bivariate) [Gaussian copula](https://en.wikipedia.org/wiki/Copula_(probability_theory)#List_of_copula_density_functions_and_applications) is an alternative which provides stable, consistent and reasonable estimates of the joint distribution model. $$C_{XY}^{Gauss} (x,y)=\frac{1}{\sqrt{1-ρ^2}} \times \exp \left (-\frac{ρ^2 (a^2+b^2 )-2abρ}{2(1-ρ^2 )} \right ),$$ where the pairwise correlation $-1 < \rho < 1$ and $$a(x)= \sqrt{2}\times erf^{-1} (2\times F_X (x)-1),$$ $$b(y)= \sqrt{2}\times erf^{-1} (2\times F_Y (y)-1),$$ and the *error function* is $$erf(s)= \frac{2}{\sqrt{\pi}}\int_{0}^{\infty}{e^{-t^2}dt}.$$ Similarly for the [copula blend of the pair of marginal density functions](https://en.wikipedia.org/wiki/Copula_(probability_theory)#Mathematical_derivation_of_copula_density_function): $$f_{XY}(x,y) = {\partial^2 F_{XY}(x,y) \over\partial x\,\partial y }= {\partial^2 C(F_X(x),F_Y(y)) \over\partial x\,\partial y} =\\ {\partial^2 C(u,v) \over\partial u\,\partial v} \times {\partial F_X(x) \over\partial x} \times {\partial F_Y(y) \over\partial y} = c(u,v) f_X(x) f_Y(y) \\ \Longrightarrow \frac{f_{XY}(x,y)}{f_X(x) f_Y(y) } = c(u,v).$$ ### Kynamic Modeling of Bacterial Growth using Gaussian Copula Let's now return to the kime dynamical system of bacteria growth/decay of a specific bio-population. Recall that for some (growth) parameter $a > 1$, the kime radial growth of the bacteria colony is modeled by an exponential function $f(t)=A_o a^t$. The kime-phase population change can be modeled as a stochastic process, $\phi \sim \Phi[-\pi, \pi)$, for some symmetric kime-phase distribution, $\Phi$. Let's assume for a fixed time point, $t$, the population size kynamics in phase-space is controlled by the function $$g(\phi)=\left ( \frac{1}{2a}\right )^{\frac{\phi}{\pi}}:[-\pi, \pi) \to \left [ \frac{1}{2a}, 2a \right ].$$ Note that the effect of this phase-model function is random ranging from shrinking by half to doubling the population. Of course, this phase-space dynamics corresponds to a fixed kime-magnitude (time) and depends on the phase distribution $\Phi[-\pi, \pi)$. Therefore, the kynamics of the bacteria population can be expressed via a (Gaussian) copula: $$F(t,\phi)= f(t)\ g(\phi)\ C_{T \Phi}^{Gauss} (t,\phi) \Longrightarrow$$ $$F(t,\phi)= \underbrace{A_o a^t}_{f(t)} \times \underbrace{\left ( \frac{1}{2a}\right )^{\frac{\phi}{\pi}}}_{g(\phi)} \times \underbrace{\frac{1}{\sqrt{1-ρ^2}} \times \exp \left (-\frac{ρ^2 (a^2+b^2 )-2abρ)}{2(1-ρ^2 )} \right )}_{\text{Gaussian copula}}.$$ We can now illustrate alternative trajectories of the dynamics of the bacteria population: - Classical Time Dynamics: $\phi=0$. - Path Kynamics over curves $\gamma:\mathbb{R}\to \mathbb{C}$. - Full kime-dynamics as kime-surfaces over the entire complex time domain. Let's examine each of these three progressively more complex models, starting with the classical time dynamics: $\phi=0$. The interactive plot below shows the shapes of the main model functions: $f(s)$, $g(\phi)$, $C_{S \Phi}^{Gauss} (s,\phi)$, and $F(s,\phi)$, over time $s\sim S$ and kime-phase $\phi\sim\Phi$. As the radial growth is exponential, we display the logarithm of the joint distribution, $\ln(F(s,\phi))$, which appears linear in this plot. ```{r eval=T, message=FALSE, warning=FALSE, error=FALSE} library(extraDistr) library(pracma) # set parameters a <- 2 # base of population exponential growth/decay w.r.t. time, a>0 Ao <- 100 # initial bacterial colony size at time t=0, Ao >0 rho <- 0.3 # initial dependence (correlation), -1 < rho = Corr(t,phi) < 1 len <- 100 s <- seq(0, 10, length=len) # Time parameter grid phi <- seq(-pi, pi, length=len) # Phase parameter grid # Base function definitions # f(t) f <- function (t=0) { # input is argument-value (s[t]), not argument-index (t) # if (t < s[1]) { val = 0 } # else if (t > s[len]) { val = 0 } # else { val=Ao * a^t } return ( Ao * a^t ) } # f(-1); f(s[1]); f(s[len]); f(s[len]+10) # Note: since "integral" expects a vectorized smooth function, which our densities # are not, we need a work around, piece-wise integration? F1 <- function (x) { # inputs are arguments-indices # CDF as an abstract function # message(paste('input to F1 =',x)) if (x < 1) { val = 0 } else if (x > len) { val = 1 } else { # normalize the CDF to integral(f) val= integral(f, xmin=s[1], xmax=s[x], reltol=1e-12, method="Kronrod")/147588 } return(val) } # F1(-1); F1(1); F1(len); F1(len+10) # g(phi) g <- function (ph=0) { # input is argument-value (phi[t]) # Any other symmetric distribution over a finite domain can be used, see # https://en.wikipedia.org/wiki/Symmetric_probability_distribution#Partial_list_of_examples # draw random Phases: phiNew ~ Phi[-pi, pi) # l = rlaplace(1000, mu=0, sigma=0.5) # hist(l) # Laplace Distribution phiNew = rlaplace(length(phi), mu=0, sigma=0.5) phi <- phiNew return ( (1/(2*a))^(ph/pi) ) } # g(-1); f(phi[1]); f(phi[len]); f(phi[len]+10) gMax <- max(abs(g(phi))) fMax <- max(abs(f(s))) G1 <- function (x) { # inputs are arguments-indices # CDF as an abstract function if (x < 1) { val = 0 } else if (x > len) { val = 1 } else { # normalize the CDF to integral(f) val= integral(g, xmin=phi[1], xmax=phi[x], reltol=1e-12, method="Kronrod")/903678 } return(val) } # G1(-1); G1(1); G1(len); G1(len+10) # Gauss-copula # C_{XY}^{Gauss} (x,y)=\frac{1}{\sqrt{1-ρ^2}} \times \exp \left (-\frac{ρ^2 (a^2+b^2 )-2abρ)}{2(1-ρ^2 )} \right ), # where the pairwise correlation $-1 < \rho < 1$ and # a(x)= \sqrt{2}\times erf^{-1} (2\times F_X (x)-1), # b(y)= \sqrt{2}\times erf^{-1} (2\times F_Y (y)-1), # and the error function is # $$erf^{-1}(s)= \frac{2}{\sqrt{\pi}}\int_{0}^{\infty}{e^{-t^2}dt}.$$ CGauss <- function (t=len/2, s=len/2) { # arguments are indices! if (t < 1) { A= -1e+7 } else if (t >= len) { A= 1e+7 } else { A = sqrt(2) * erfinv(2*(F1(t))-1) } if (s < 1) { B= -1e+7 } else if (s > len) { B= 1e+7 } else { B = sqrt(2) * erfinv(2*(G1(s))-1) } # message(paste('A=', A, 'B=', B)) val = (1/(sqrt(1-rho^2))) * exp(-(rho^2*(A^2+B^2) - 2*A*B*rho)/(2*(1-rho^2))) if (is.finite(val)) return (val) else return (0) } # CGauss(50, 50) # F(t) fJoint <- function (t=len/2, ph=len/2, rho = rho) { return ( f(s[t]) * g(phi[ph]) * CGauss(t, ph) ) } # fJoint(50, 50) y1 <- f(s) # f: first marginal y2 <- g(phi) # g: second marginal y3 <- c() # Gaussian Copula (. , s=sMax/2) for (i in 1:len) y3[i]=CGauss(i, 55) y4 <- c() # joint density f(t,t) for (i in 1:len) y4[i] <- fJoint(t=i, ph=55, rho) plot_ly() %>% add_trace(x=~s, y=~log(y1), type='scatter', mode='markers+lines', name="log(f): log(Longitudinal marginal)") %>% add_trace(x=~s, y=~y2, type='scatter', mode='markers+lines', name="g: Kime-phase marginal", opacity=1) %>% add_trace(x=~s, y=~y3, type='scatter', mode='markers+lines', name="C_Gauss(t,s): Gaussian Copula", opacity=0.5) %>% add_trace(x=~s, y=~log(y4), type='scatter', mode='markers+lines', name="log(f(t,s, rho)): Joint Kime-time and Kime-phase Effects") %>% layout(title = "Copula-based Modeling of Joint Distribution", xaxis = list(title = "time"), yaxis = list(title = "Functional Value"), # showlegend = FALSE, legend = list(orientation = 'h'), scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1))) ``` Next we can display the raw and smoothed versions of the joint distribution of bacteria growth in kime as flat 2D images. ```{r errorTesting, message=FALSE, warning=FALSE, error=FALSE} #### Path Kynamics over curves \gamma:\mathbb{R}\to \mathbb{C} ### Part 2: Path Kynamics over curves \gamma:\mathbb{R}\to \mathbb{C} ### Kime tube gamma(k1,k2) on torus ########## CARTESIAN #### Parametric curve on torus # define kime-surface s <- seq(0, 10, length=100) # Time parameter grid phi <- seq(-pi, pi, length=100) # Phase parameter grid xsurf <- s # longitudinal trajectory ysurf <- phi # phase trajectory zsurf <- mat<-matrix(0, nrow=len, ncol=len) # joint density f(t,s) for (i in 1:len) { for (j in 1:len) { zsurf[i,j] <- (fJoint(t=i, ph=j, rho))/(10^4) } if (i %% 30 ==0) message(paste0('Progress: ... ', round((100*i)/len, 0), '%')) } # image(zsurf) plot_ly(z=~zsurf, type="heatmap") # %>% hide_colorbar() # image(zsurf) ``` The next plot may be somewhat counter-intuitive, as it intentionally renders the joint distribution in *Cartesian coordinates*, $\ln(F(s,\phi))$, instead of poloar coordinates. However, it shows the individual marginal density dynamics over time $s$ (steady growth) and over kime-phase $\phi$ (stochastic). The surface display also includes overlays of a pair of canonical curves corresponding to time-dynamics for a fixed phase $\phi=0$ (orange), and phase-dynamics for a fixed time $s=s(t=50)=5$ (green). ```{r errorTesting2, message=FALSE, warning=FALSE, error=FALSE} # 3D Kime-surface Plot # 3D Kime-surface Plot fnt <- list(family = "Courier New, monospace", size = 18, color = "black") xlab <- list(title = "time", titlefont = fnt) ylab <- list(title = "phase", titlefont = fnt) zlab <- list(title = "Kime-surface", titlefont = fnt) plot_ly(x=~xsurf, y=~ysurf, z=~t(zsurf), type="surface", opacity=0.95, showlegend=FALSE, name="Joint Kime-Density-Surface") %>% # add vital days after surgery projection curve, for a fixed BMI=22 add_trace(x=~xsurf, y=ysurf[22], z=~zsurf[, 22], type="scatter3d", mode="lines", line=list(width=20), name="Phi=-1.8") %>% # add BMI projection curve vital days after surgery=1000 add_trace(x=xsurf[70], y=~ysurf, z=~zsurf[70,], type="scatter3d", mode="lines", line=list(width=20), name="Time=7") %>% hide_colorbar() %>% layout(title = "Joint Kynamical System Simulation\n (Gaussian copula, Rho=0.3)", scene = list(aspectmode = "manual", xaxis = xlab, yaxis = ylab, zaxis = zlab, aspectratio = list(x=1, y=1, z=1))) # 1:1:1 aspect ratio ``` 2D image representations of the full kynamics are not particularly useful, other than illustrating the natural polar representation of complex-time in terms of radial time and angular phase. ```{r errorTesting3, message=FALSE, warning=FALSE, error=FALSE} #### Full kime-dynamics as kime-surfaces over the entire complex time domain. ###### POLAR Coordinates ... # Change basis from from polar coordinates to Cartesian coordinates matrix_ON <- matrix(0, nrow = 101, ncol = 101) for (t in 1:length(s)) { for (p in 1:length(phi)) { x = as.integer(51+5*s[t]*cos(phi[p])) y = as.integer(51+5*s[t]*sin(phi[p])) if (x>100) x=100 if (x<1) x=1 if (y>100) y=100 if (y<1) y=1 matrix_ON[x,y] <- zsurf[t,p] # zsurfSmooth3$z # zsurfSmooth2[t,p] # zsurfSmooth # zsurf } } dim(matrix_ON) # 101 * 101 # image(t(matrix_ON), col = hcl.colors(100, "terrain"), axes = FALSE) plot_ly(z= matrix_ON, type="heatmap") # Inspect visually matrix_ON # library(gplots) # heatmap.2(matrix_ON, Rowv=FALSE, Colv=FALSE, dendrogram='none', cellnote=matrix_ON, notecol="black", trace='none', key=FALSE, lwid = c(.01, .99), lhei = c(.01, .99), margins = c(5, 15 )) # smooth/blur the matrices library (spatstat) matrix_ON_smooth <- as.matrix(blur(as.im(matrix_ON), sigma=1.5)) # heatmap.2( matrix_ON_smooth, Rowv=FALSE, Colv=FALSE, dendrogram='none', cellnote=matrix_ON, notecol="black", trace='none', key=FALSE, lwid = c(.01, .99), lhei = c(.01, .99), margins = c(5, 15 )) ``` The ultimate kime-surface representation of the bacteria kynamics in kime can be visualized as a raw (not very useful) or a (heavily) smoothed (more intuitive) surface. ```{r message=FALSE, warning=FALSE, error=FALSE} #library(plotly) plot_ly(z = ~matrix_ON_smooth, type = "surface", showscale=FALSE) %>% hide_colorbar() %>% layout(title = "Smoothed Surface") plot_ly(z = ~matrix_ON, type = "surface") %>% hide_colorbar() %>% layout(title = "Raw Surface") ``` Now let's superimpose the canonical *time* and *phase* curves. The *canonical time curve* represents one typical instantiation of the observed temporal process (for a fixed phase, e.g., $\phi=0$). The *canonical phase curve* depicts repeated IID random sampling, in a controlled experiment, for one fixed time point, e.g., $time=7$, where different phase arguments correspond to independent sampling instances of the bacteria population at that time. ```{r eval=T, message=FALSE, warning=FALSE, error=FALSE} # Plot matrix_ON_smooth matToPlot <- matrix_ON_smooth # matrix_ON # Generate Time-curve (phi=0) t <- 1:length(s) curveS_x <- as.integer(51+5*s[t]*cos(phi[1])) curveS_y <- as.integer(51+5*s[t]*sin(phi[1])) curveS_z <- (matToPlot[curveS_x, curveS_y])[ , 1] # plot_ly( x=~curveS_x, y=~curveS_y, z=~curveS_z, # type="scatter3d", mode="lines") # Generate phase trajectory (t=7) curveP_x <- as.integer(51+5*s[70]*cos(phi)) curveP_y <- as.integer(51+5*s[70]*sin(phi)) curveP_z <- (matToPlot[curveP_x, curveP_y])[ , 1] # plot_ly( x=~curveP_x, y=~curveP_y, z=~curveP_z, # type="scatter3d", mode="lines") # Project Curves/paths onto 3D Kime-surface Plot # add projection curve phi=0. Plot matrix_ON_smooth? , # opacity=0.1, showscale=FALSE plot_ly(z = ~t(matToPlot), type = "surface", opacity=0.2) %>% # Add time-series (phi=0) add_trace(x=~curveS_x, y=~curveS_y, z=~curveS_z+0.2, type="scatter3d", mode="lines", line=list(width=20), name="Phi=0") %>% # Add phase trajectory (time=5) add_trace(x=~curveP_x, y=~curveP_y, z=~curveP_z+0.2, type="scatter3d", mode="lines", line=list(width=20, color="blue"), name="Time=5") %>% hide_colorbar() %>% layout(title = "Canonical Kime Curves (Smooth Surface)") # Plot matrix_ON? matToPlot <- matrix_ON # matrix_ON_smooth # Generate Time-curve (phi=0) t <- 1:length(s) curveS_x <- as.integer(51+5*s[t]*cos(phi[1])) curveS_y <- as.integer(51+5*s[t]*sin(phi[1])) curveS_z <- (matToPlot[curveS_x, curveS_y])[ , 1] # plot_ly( x=~curveS_x, y=~curveS_y, z=~curveS_z, # type="scatter3d", mode="lines") # Generate phase trajectory (t=7) curveP_x <- as.integer(51+5*s[70]*cos(phi)) curveP_y <- as.integer(51+5*s[70]*sin(phi)) curveP_z <- (matToPlot[curveP_x, curveP_y])[ , 1] # plot_ly( x=~curveP_x, y=~curveP_y, z=~curveP_z, # type="scatter3d", mode="lines") # Project Curves/paths onto 3D Kime-surface Plot # add projection curve phi=0. Plot matrix_ON_smooth? , # opacity=0.1, showscale=FALSE plot_ly(z = ~matToPlot, type = "surface", opacity=0.2) %>% # Add time-series (phi=0) add_trace(x=~curveS_x, y=~curveS_y, z=~curveS_z+0.2, type="scatter3d", mode="lines", line=list(width=20), name="Phi=0") %>% # Add phase trajectory (time=5) add_trace(x=~curveP_x, y=~curveP_y, z=~curveP_z+0.2, type="scatter3d", mode="lines", line=list(width=20, color="blue"), name="Time=7") %>% hide_colorbar() %>% layout(title = "Canonical Kime Curves (Raw Surface)") # # fix the Plot_Ly Text Labels # x <- vector() # y <- vector() # i <- 1 # for (t in 1:length(s)) { # for (p in 1:length(phi)) { # x = as.integer(51+5*t*cos(phi[p])) # y = as.integer(51+5*t*sin(phi[p])) # i <- i+1 # } # } # # hoverText <- cbind(x=1:101, y=1:101, height=as.vector(t(matrix_ON_smooth))) # tail(mytext) # custom_txt <- matrix(NA, nrow=101, ncol=101) # hoverTextON <- cbind(x=1:101, y=1:101, height=as.vector(t(matrix_ON_smooth))) # tail(mytext) # custom_txtON <- matrix(NA, nrow=101, ncol=101) # # for (x in 1:101) { # for (y in 1:101) { # t = sqrt((x-51)^2 + (y-51)^2) # p = atan2(y-51, x-51) # custom_txt[x,y]<-paste(' fMRI: ', round(hoverText[(x-1)*101+y, 3], 3), # '\n time: ', round(t, 0), # '\n phi: ', round(p, 2)) # custom_txtON[x,y]<-paste(' fMRI: ', # round(hoverTextON[(x-1)*101+y,3],3), # '\n time: ', round(t, 0), # '\n phi: ', round(p, 2)) # } # } # dim(custom_txt); custom_txt[10:12,10:12]; custom_txtON[10:12,10:12] # # xx2 <- 51 + c(-50:50) %o% cos(seq(-pi, pi, 2*pi/20)) # yy2 <- 51 + c(-50:50) %o% sin(seq(-pi, pi, 2*pi/20)) # #zz2 <- as.vector(matrix_ON_smooth) %o% rep(1, 21*21) # zz2 <- matrix_ON_smooth # # #plot 2D into 3D and make the text of the diameter (time), height (r), and phase (phi) # f <- list(family = "Courier New, monospace", size = 18, color = "black") # x <- list(title = "k1", titlefont = f) # y <- list(title = "k2", titlefont = f) # z <- list(title = "fMRI Kime-series", titlefont = f) # # # Plot ON kime-surface # plot_ly(x = ~xx2, y = ~yy2, z = ~zz2, type = "surface", # scatterpolar # text = custom_txt, hoverinfo = "text", showlegend = FALSE) %>% # #add_trace(x=~xx2, y=~yy2, z=~ww2, colors = c("blue", "yellow"), # # type="surface", text = custom_txtOFF, hoverinfo = "text", # # opacity=0.3, showscale = FALSE, showlegend = FALSE) %>% # # trace the main Z-axis # # add_trace(x=51, y=51, z=0:5, type="scatter3d", mode="lines", # # line = list(width = 10, color="red"), name="Space(x)", # # hoverinfo="none", showlegend = FALSE) %>% # layout(dragmode = "turntable", title = "ON Kime-Surface/Kime-Series at a fixed voxel location", # scene = list(xaxis = x, yaxis = y, zaxis = z), showlegend = FALSE) ``` # Symmetries, Generators and Operators ## Symmetry Symmetry relates to the action of various transformations (e.g., linear operators), which plays a key role in physics and quantum mechanics. The identity transformation $I$ on a space $V$ is a mapping $I:V\to V$ that changes nothing at all, $I(v)=v,\ \forall v\in V$. Heuristically, symmetry can be explained as a transformation that leaves the object unchanged, e.g., you can't tell that a circle in the plane, or a perfect ball in 3D, is rotated as these have angular symmetries. On the other hand, a square in 2D, or a cube in 3D, have very specific rotational symmetries (related to multiples of $90^o$ spherical coordinate rotations). ### Groups А symmetry transformation leaves the object of interest unchanged. Group theory is important for representing, modeling, and computing symmetries. А group, $G$, is а set of elements (e.g., transformations) along with an group-operation (element product, e.g., rotation composition) subject to some group axioms: - Closure: for each pair of elements $a,b\in G$, $a\cdot b \in G$, - Associativity: $\forall a,b,c\in G$, $(a\cdot b)\cdot c=a\cdot (b\cdot c)$, - Identity element: $\exists! I \in G$, such that $\forall a \in G$, $I\cdot a=a \cdot I = a$, - Inverse element: $\forall a \in G$, $\exists a^{−1}\in G$, such that $a^{-1}\cdot a= a\cdot a^{-1} =I$ (identity element). Note that there are *continuous* symmetries, such as rotations of а 3D ball, and a *discrete* symmetries, such mirror-imaging and integer rotations by $90^o$ for a 3D cube. A special property of continuous symmetries makes them especially important - they have elements that approximate arbitrarily close to the identity element $I$. For instance, all rotational symmetries of а perfect 3D Ball. We can get as close and precise to an identity rotation, $I=0^o$, of a 3D ball, by small but non-trivial rotations $I\approx \epsilon 1^o$, where $\epsilon\to 0^+$. This is not the case for discrete symmetries, e.g., thе set rotational transformations that leave а 2D square unchanged (*invariant*) include the 4 rotations by $0^o,90^0, 180^o, 270^o$, along with mirror-imaging. ### Infinitesimal transformations and generators A group element $g\in G$ is close to the identity $I$ when $g(\epsilon) = I + \epsilon \mathcal{G}$, where $\epsilon >0$ and $\mathcal{G}$ is a *generator*. These group elements represent *infinitesimal transformations* and in isolation, they change very little the objects they act on. However, repeating the application of an infinitesimal transformation multiple times would represent another element in the group $$g^n(\epsilon)=\underbrace{g(\epsilon)\cdot g(\epsilon)\cdot \ ...\ \cdot g(\epsilon)}_{n} = (I + \epsilon \mathcal{G})^n,$$ whose effect may be much more pronounced. For instance, repeating a tiny rotation several times (in the same direction) results in a *large* rotation. SUppose we are given some fixed finite transformation parameter, $\theta$, e.g., $\theta= 55^o$, then we can choose a large $n\in \mathbb{N}$ so that $\epsilon = \frac{\theta}{n}\to 0^+$. Recalling the definition of the Euler number $e$ and the corresponding exponential function, $e^x$, we have: $$h(\theta) =\lim_{n\to\infty} {\left (I+\underbrace{\frac{\theta}{n}}_{\epsilon}\mathcal{G}\right )^n} \equiv e^{\theta\mathcal{G}}.$$ Тhе *generator* $\mathcal{G}$ effectively creates in the limit the new finite transformation $h\in G$. We will examine 3 types of generators: - *Spatial translation*: $\mathcal{G}_x^{trans}=\partial_x \equiv \frac{\partial}{\partial x}$, - *Temporal translation* (kime-magnitude): $\mathcal{G}_t^{trans}=\partial_t \equiv \frac{\partial}{\partial t}$, and - *Phase translation* (kime-phase): $\mathcal{G}_{\phi}^{trans}=\partial_{\phi} \equiv \frac{\partial}{\partial \phi}$. Let's start with а function $$f=f\left (\overbrace{x}^{space}, \overbrace{\underbrace{\kappa}_{(t,\phi)}}^{kime}\right )=f(x,t,\phi).$$ ### Spatial translation We aim to prove that to generate а *spatial translation* $T$, such that $Т f(x,t,\phi) = f(x+a,t,\phi)$, we need to use the *spatial translation generator* $G_x^{trans}=\partial_x$. Given any translation parameter $\theta\equiv a\in \mathbb{R}$, we can use the Taylor series expansion of the exponential function, $e^x = \sum_{n=0}^{\infty}{\frac{x^n}{n!}}$, to expand right hand size for $\mathcal{G}\equiv \mathcal{G}_x^{trans}$: $$e^{\theta\mathcal{G}_x^{trans}} f(x,t,\phi) \equiv e^{a\mathcal{G}_x^{trans}} f(x,t,\phi)=$$ $$ \left ( 1 + a\mathcal{G}_x^{trans} +\frac{a^2}{2!}(\mathcal{G}_x^{trans})^2 + \cdots \frac{a^k}{k!}(\mathcal{G}_x^{trans})^k + \cdots \right ) f(x,t,\phi)= $$ $$\left ( 1 + a\partial_x +\frac{a^2}{2!}(\partial_x)^2 + \cdots \frac{a^k}{k!}(\partial_x)^k + \cdots \right ) f(x,t,\phi)\equiv f(x + a,t,\phi).$$ The last equation is just the Taylor expansion of $f$ at $(x+a,t,\phi)$. Note that $(\partial_x)^k f(x,t,\phi)\equiv \frac{\partial^k}{\partial x^k}f$ is the *rate of change* of $f(x,t,\phi)$ in the $x$-direction. Hence, multiplying this rate of change by a constant distance parameter $a$ propagates the function forward in the $x$-direction. To a first (linear) approximation, the total corresponding change of $f(x,t,\phi)$ is $a \partial_x f(x,t,\phi)$ and the value of $f$ at the location $(x+a,t,\phi)$ is (approximately) $f(x,t,\phi)+a \partial_x f(x,t,\phi)$. ### Temporal translation Similarly, to generate а *temporal translation* $Т f(x,t,\phi) = f(x,t+b,\phi)$, we need to use the *temporal translation generator* $\mathcal{G}_t^{trans}=\partial_t$. For a translation parameter $b\in \mathbb{R}^+$ we expand right hand size for $\mathcal{G}\equiv \mathcal{G}_t^{trans}$: $$e^{b\mathcal{G}_t^{trans}}= \left ( 1 + b\mathcal{G}_t^{trans} +\frac{b^2}{2!}(\mathcal{G}_t^{trans})^2 + \cdots \frac{b^k}{k!}(\mathcal{G}_t^{trans})^k + \cdots \right ) f(x,t,\phi)= $$ $$\left ( 1 + b\partial_t +\frac{b^2}{2!}(\partial_t)^2 + \cdots \frac{b^k}{k!}(\partial_t)^k + \cdots \right ) f(x,t,\phi)\equiv f(x ,t+b,\phi). $$ The last equation is the Taylor expansion of $f$ at $(x,t+b,\phi)$. Note that $(\partial_t)^k f(x,t,\phi)\equiv \frac{\partial^k}{\partial t^k}f$ is the *rate of change* of $f(x,t,\phi)$ in *time* (kime-magnitude). Multiplying this rate of change by a constant time increment $b$ propagates the function forward in time and an approximation to the total corresponding change of $f(x,t,\phi)$ is $b \partial_t f(x,t,\phi)$. ### Kime-Phase translation Finally, we generate а *kime (phase) translation* $Т f(x,t,\phi) = f(x,t,\phi +c)$ using the *phase translation generator* $\mathcal{G}_{\phi}^{trans}=\partial_{\phi}$. For a fixed phase-translation parameter $c\in [\pi,\pi)$, we expand right hand size for $\mathcal{G}\equiv \mathcal{G}_{\phi}^{trans}$: $$e^{c\mathcal{G}_{\phi}^{trans}}= \left ( 1 + c\mathcal{G}_{\phi}^{trans} +\frac{c^2}{2!}(\mathcal{G}_{\phi}^{trans})^2 + \cdots \frac{c^k}{k!}(\mathcal{G}_{\phi}^{trans})^k + \cdots \right ) f(x,t,\phi)= $$ $$\left ( 1 + c\partial_{\phi} +\frac{c^2}{2!}(\partial_{\phi})^2 + \cdots \frac{c^k}{k!}(\partial_{\phi})^k + \cdots \right ) f(x,t,\phi)\equiv f(x ,t,\phi +c). $$ Again, the last equation is just the Taylor expansion of $f$ at $(x,t,\phi +c)$ and $(\partial_{\phi})^k f(x,t,\phi)\equiv \frac{\partial^k}{\partial {\phi}^k}f$ is the *rate of change* of $f(x,t,\phi)$ in kime-phase direction. Summarizing, multiplying this rate of change by a constant phase increment $c$ propagates the function forward in the kime-phase direction and an approximation to the total corresponding change of $f(x,t,\phi)$ is $c \partial_{\phi} f(x,t,\phi)$. In essence, we describe all continuous symmetries by acting iteratively, with infinitesimal transformations corresponding to a specific generator, $\mathcal{G}_{(\cdot)}^{trans}$, on the wavefunction, $f(x,t,\phi)$. ## Operators ### Bra-Kets The *state vector*, **ket**, has the mathematical notation $|\Psi\rangle$, which represents all information (tractable, measurable quantities we can observe) about a system. For each observable quantity (e.g., energy, momentum, position, spin, etc.), every different possible state of the system corresponds to а different ket $\{\ |\Psi_{\alpha}\rangle\ \}_{\alpha}$. As an example, an electron rotating around an atomic nucleus can have different observable *energies*; from the ground state with the lowest energy, $|\Psi_1^e\rangle$, to excited higher-energy states, $|\Psi_2^e\rangle$, $|\Psi_3^e\rangle$, etc. The same electron has observable *positions* which can be tracked by $\{\ |\Psi_{\alpha}^x\rangle\ \}$; from initial position at $x= О$, $|\Psi_{0}^x\rangle$, to other positions such as $x= 0.123$, $|\Psi_{0.123}^x\rangle$ and $x= 5$, $|\Psi_{5}^x\rangle$. Note that some observable quantities may have discrete, finite, countably, or uncountably many possible states. The *bra* is the *Hermitian conjugate* (equivalent to transposed, $(\cdot)^t$, complex conjugate, $(\cdot)^*$) of the *ket* $\langle \Phi | \equiv |\Phi\rangle^{\dagger} \equiv (|\Phi\rangle^{*})^{t}$. The *ket* notation allows us to extract/compute/represent specific information about the state of the system. As an example, consider the system *momentum* (product of mass and velocity, $p=m\ v$) of our system. The momentum states can be expressed as eigenvalues of the momentum operator and the momentum eigenvector kets $$\overbrace{\underbrace{\hat{p}}_{\text{operator}}}^{\text{momentum}}|\Psi^{p}\rangle = \underbrace{p}_{\mathbb{R}\ \text{eigenvalue}} \overbrace{|\Psi^{p}\rangle}^{\text{eigenvector}}$$ As the problem context identifies the observable quantity we are interested in, going forward, we will be suppressing the superindex and writing these as $\hat{p}|\Psi\rangle = {p}|\Psi\rangle$. The corresponding momenta at momentum state configurations $|\Psi_1\rangle$ and $|\Psi_2\rangle$ are $p_1,p_2\in \mathbb{R}$, i.e., the two possible observable momenta, $\hat{p}|\Psi_1\rangle = {p_1}|\Psi_1\rangle$ and $\hat{p}|\Psi_2\rangle = {p_2}|\Psi_2\rangle$. In practice, due to Heisenberg uncertainty, the system is rarely in a single *base state*. Rather it's mostly in superposition (linear combinations) of *base states*. The intrinsic system uncertainty suggests that the momentum changes every time we measure position and vice-versa, the position changes each time we measure the momentum. This implies that the order of measuring certain quantities matters; just like washing and drying our hands come in a natural order (the reverse order yields a different outcome). Specifically, $$\hat{p}\hat{x} |\Psi \rangle - \hat{x}\hat{p} |\Psi \rangle =\overbrace{\underbrace{(\hat{p}\hat{x}-\hat{x}\hat{p})}_{[\hat{x},\hat{p}]}}^{\text {commutator}}|\Psi \rangle =$$ $$-i\hbar\partial_x \hat{x}|\Psi \rangle+ \hat{x}i\hbar\partial_x |\Psi \rangle = \underbrace{-(i\hbar\partial_x \hat{x})|\Psi \rangle -\hat{x}(i\hbar\partial_x |\Psi \rangle )}_{\text{product diff rule}} + \hat{x}i\hbar\partial_x |\Psi \rangle =$$ $$-i\hbar|\Psi \rangle \not=0.$$ Note that in principle, $\hat{p}_i\hat{x}_j |\Psi \rangle =-i\hbar\delta_{i,j}|\Psi \rangle \not=0$ , as $i,j\in\{1,2,3\}$. Also, the rate of change of the system with respect to space $x$, time $t$ or kime-phase $\phi$ is defined as $\partial_{\alpha}|\Psi \rangle$, where $\alpha\in\{x,t,\phi\}$. ### Schrodinger equation As we showed above, $\partial_t$ is a multiple of the quantum energy operator $\hat{E}\equiv -i\hbar\partial_t$, i.e., $\hat{E}|\Psi \rangle \equiv -i\hbar\partial_t |\Psi \rangle =E|\Psi \rangle.$ Total energy in classical mechanics is the total sum of *kinetic energy* $T$ and *potential energy* $V$, $E = T + V$. The kinetic energy is $T= \frac{1}{2}mv^2 =\frac{p^2}{2m}$, where the momentum $р = mv$, $m$ is the mass, and $v$ is the velocity. By replacing the momentum $p$ by the momentum operator, $\hat{p}=-i\hbar\partial_t$,we transform the classical kinetic energy into а quantum operator $\hat{T}=\frac{\hat{p}^2}{2m}\equiv \frac{(-i\hbar\partial_t)^2}{2m}=-\frac{\hbar^2 \partial_t^2}{2m}$. Similarly, we can express the classical potential energy, which usually only depends on the spatial position $x$ (not time or kime-phase), $V = V(x)$, as an operator $\hat{V}=V(\hat{x})$. To derive the [Shrodinger equation](https://en.wikipedia.org/wiki/Schr%C3%B6dinger_equation), we transform the total energy equation, $E=T+V$, to an operator equation, $\hat{E} = \hat{T} + \hat{V}$. Starting with $\hat{E}\equiv -i\hbar\partial_t$, we can expand the terms as follows: $$\hat{E}|\Psi \rangle \equiv -i\hbar\partial_t |\Psi \rangle =E|\Psi \rangle =(\hat{T} + \hat{V})|\Psi \rangle=$$ $$\left (\overbrace{\underbrace{-\frac{\hbar^2 \partial_t^2}{2m} +V(\hat{x})}_{\hat{H}}}^{\text{Hamiltonian operator}} \right )|\Psi \rangle= -\frac{\hbar^2 \partial_t^2}{2m}|\Psi \rangle +V(\hat{x})|\Psi \rangle .$$ Many problems require solving the Schrodinger equation for a specific potential $V(x)$ under some given boundary conditions. ### Operator moments Designing a controlled experiment involving a system of a pair of particles we may sometimes measure a momentum value $p_1$ and other times a value $p_2$. The superposition of the 2 base states is expressed as $$|\Psi\rangle = \underbrace{a_1}_{\in\mathbb{C}}|\Psi_1\rangle + \underbrace{a_2}_{\in\mathbb{C}}|\Psi_2\rangle,$$ where $|a_1|^2,|a_2|^2\in\mathbb{R}$ represent the probability distribution (odds) of measuring momenta $p_1$ and $p_2$, respectively. In the general *discrete case*, $$|\Psi\rangle = \sum_{i} {\left ( \underbrace{\langle \Psi_i | \Psi\rangle}_{\text{coeff}\ \in\mathbb{C}} \overbrace{|\Psi_i\rangle}^{\text{base}} \right )}.$$ Similarly, in the *continuous case*, the summation transforms to integration $$|\Psi\rangle = \int_{\mathbb{R}} {\left ( \underbrace{\langle p | \Psi\rangle}_{\text{coeff}\ \in\mathbb{C}} |p\rangle dp \right )}\equiv \int_{\mathbb{R}} {\Psi(p) |p\rangle dp }.$$ In our simple Bernoulli experiment, the momentum operator $\hat{p}$ acting on the momentum state $|\Psi\rangle$ will be expressed as a linear operator acting on this superposition: $$\hat{p}|\Psi\rangle =\hat{p}\left ( a_1|\Psi_1\rangle + a_2|\Psi_2\rangle \right )= a_1\hat{p}|\Psi_1\rangle + a_2\hat{p}|\Psi_2\rangle= a_1p_1|\Psi_1\rangle + a_2p_2|\Psi_2\rangle.$$ This representation makes sense in a statistical framework. Once an observation is made, the system yields a stochastic measure of the momentum - in this case one of two possible momentum values $p_1$ or $p_2$. In general, there may be infinitely many, or even continuous, state outcomes. Extending the [notions of statistical moments](https://wiki.socr.umich.edu/index.php/AP_Statistics_Curriculum_2007_Distrib_MeanVar) (e.g., first moment is the *mean* and the second moment is the *variance*), we can define *operator moments*: - *Operator mean*: Weighted average of the possible observable momenta. For our simple momentum Bernoulli experiment, $\langle \hat{p}\rangle \equiv \langle p\rangle =\underbrace{|a_1|^2}_{\text {probability}} \overbrace{p_1}^{\text {momentum}} + |a_2|^2 p_2$. - *Operator variance*: Weighted average of the squared deviations of the possible observable momenta from the mean. For our simple momentum Bernoulli experiment, $$(\Delta{p})^2\equiv \langle \hat{p}^2\rangle -\langle \hat{p}\rangle^2 \equiv \left ( \underbrace{|a_1|^2}_{\text {probability}} \overbrace{p_1^2}^{\text {momentum}} + |a_2|^2 p_2^2 \right ) -\left (|a_1|^2 p_1 + |a_2|^2 p_2 \right ).$$ Based on *orthonormality* assumptions, $\langle \Psi_i | \Psi_j\rangle =\delta_{i,j}=\left\{ \begin{array}{ll} 1, i=j\\ 0, {\text{otherwise}} \end{array} \right .$, given a concrete basis $|\Psi\rangle$, e.g., momentum $|\Psi^p\rangle$ or position $|\Psi^x\rangle$, the operator moments are computed as follows: $$\underbrace{\langle\Psi |}_{\text{bra}} \overbrace{\hat{p}}^{\text{operator}} \underbrace{| \Psi \rangle}_{\text{ket}} \equiv (a_1^*\langle \Psi_1 | +a_2^*\langle \Psi_2 |) \ \hat{p}\ (a_1 | \Psi_1\rangle +a_2 | \Psi_2 \rangle)=$$ $$(a_1^*\langle \Psi_1 | +a_2^*\langle \Psi_2 |) (a_1 \hat{p}| \Psi_1\rangle +a_2 \hat{p}| \Psi_2 \rangle)=$$ $$(a_1^*\langle \Psi_1 | +a_2^*\langle \Psi_2 |) (a_1 p_1| \Psi_1\rangle +a_2 p_2| \Psi_2 \rangle)=$$ $$|a_1|^2 \underbrace{\langle \Psi_1 | \Psi_1\rangle}_{\text{1}} + a_1^*a_2\overbrace{\langle \Psi_1 | \Psi_2\rangle}^{0} + a_2^*a_1\langle \Psi_2 | \Psi_1\rangle + |a_2|^2\langle \Psi_2 | \Psi_2\rangle=$$ $$|a_1|^2 p_1 + |a_2|^2 p_2.$$ # Curvature One approach to capture the change of a function, such as a real or complex valued surface, that depends on two or more parameters is to track the curvature of the surface, which depends on the metric defining the manifold. There are [different types of curvatures](https://en.wikipedia.org/wiki/Curvature) that capture the various properties of manifolds. Here, we will focus on surfaces (2-manifolds) and illustrate the relation between function changes with respect to kime and the [Gaussian curvature](https://en.wikipedia.org/wiki/Gaussian_curvature), which the product of the two *principal curvatures* $Κ = \kappa_1\times \kappa_2$. Given a point on the surface $\mathcal{S}=\{(t,\phi,f(t,\phi))\ | \ t>0, -\pi\leq\phi<\pi \}$, the surface normal vector is perpendicular (orthogonal) to the surface. The intersections of the surface and any normal plane containing the normal vector generates normal-section curves on the surface whose curvature represents the *normal curvature*. The extrema of all normal curvatures (min and max) corresponding to different normal sections lead to the principal curvatures ($\kappa_1$ for min and $\kappa_2$ for max). For a fixed point $\bf{p}\in\mathcal{S}$, the two surface principal curvatures $\kappa_1$ and $\kappa_2$ represent the eigenvalues of the [shape operator](https://en.wikipedia.org/wiki/Shape_operator). By tracking the min/max bending of the surface along different directions the *curvature* captures the "change" of the surface parameterized by kime. The Gaussian curvature of the surface at $\bf{p}$ can also be expressed as the determinant of the [Hessian matrix](https://en.wikipedia.org/wiki/Hessian_matrix) of the second order derivatives of the implicit function $f(t,\phi)$, or the product of the eigenvalues of the Hessian. ## First and second fundamental forms Given a surface metric $g$, the *first fundamental form* is $$\mathbb{I}\equiv (g_{{ij}})={\begin{pmatrix}g_{{11}}&g_{{12}}\\g_{{21}}&g_{{22}}\end{pmatrix}}={\begin{pmatrix}E&F\\F&G\end{pmatrix}},$$ where the tensor components are calculated as scalar product of tangent vectors $X_1$ and $X_2$: $$g_{ij}=\langle X_{i} | X_{j}\rangle \equiv X_{i}\cdot X_{j}, \ i,j\in\{1,2\}.$$ In Cartesian $(u,v)$ or polar $(t,\phi)$ coordinates, the (linear distance) *line element* $ds$ may be expressed in terms of the first fundamental form coefficients: $$ds^{2}=E du^{2}+2Fdu\ dv+Gdv^{2},$$ $$ds^2 = dt^2 + G(t,\phi) d\phi^2.$$ The (quadratic) *area element* is $dA = |Xu \times Xv| du dv$ is: $$dA=|X_{u}\times X_{v}|\ du\ dv={\sqrt {\langle X_{u}|X_{u}\rangle \langle X_{v}|X_{v}\rangle -\left\langle X_{u}|X_{v}\right\rangle ^{2}}}\,du\,dv={\sqrt {EG-F^{2}}}\,du\,dv.$$ For example, consider a parametric curve on a unitary 2D sphere parameterized (in spherical coordinates) by: $$X(u,v)={\begin{bmatrix} \cos u\sin v\\\sin u\sin v\\\cos v\end{bmatrix}}\equiv \left (\cos u\sin v\ {\bf{i}} + \sin u\sin v\ {\bf{j}} + \cos v \ {\bf{k}} \right ) ,\ 0\leq u \lt 2\pi,\ 0\leq v\lt \pi\ ,$$ where the unitary base vectors are: $$\bf{i}={\begin{bmatrix} 1\\ 0 \\ 0 \end{bmatrix}},\ \bf{j}={\begin{bmatrix} 0\\ 1 \\ 0 \end{bmatrix}},\ \bf{k}={\begin{bmatrix} 0\\ 0 \\ 1 \end{bmatrix}}\ .$$ The plot below shows the spherical surface along with several curves, including geodesics (shortest paths between two points on the manifold). ```{r eval=T, message=FALSE, warning=FALSE, error=FALSE} u <- seq(from = 0, to = 2*pi, by = ((2*pi - 0)/(200 - 1))) v <- seq(from = 0, to = pi, by = ((pi - 0)/(200 - 1))) # rendering (u,v) parametric surfaces requires x,y,z arguments to be 2D arrays # In out case, the three coordinates have to be 200*200 parameterized tensors/arrays x1 = cos(u) %o% sin(v) # x = \cos u\sin v y1 = sin(u) %o% sin(v) # y = \sin u\sin v z1 = rep(1,length(u)) %o% cos(v) # z = 1*cos(v) # p <- plot_ly(hoverinfo="none", showscale = FALSE) %>% # add_trace(x = ~x1, y = ~y1, z = ~z1, type = 'surface', opacity=0.5, visible=T) %>% # layout(title = "Spherical Surface Example", showlegend = FALSE) # p t <- 1:length(u) # Curve 1 curveS_x <- cos(u[t]) * sin(v[t]) curveS_y <- sin(u[t]) * sin(v[t]) curveS_z <- cos(v[t]) # plot_ly( x=~curveS_x, y=~curveS_y, z=~curveS_z, # type="scatter3d", mode="lines") # Curve 2 curveP_x <- tanh(t/5) * cos(t/5) curveP_y <- tanh(t/5) * sin(t/5) curveP_z <- 1/cosh(t/5) ## Meridian curveMer_x = cos(0) * sin(v) # x = \cos u\sin v curveMer_y = sin(0) * sin(v) # y = \sin u\sin v curveMer_z = cos(v) # z = 1*cos(v) ## Parallel curvePar_x = cos(u) * sin(1) # x = \cos u\sin v curvePar_y = sin(u) * sin(1) # y = \sin u\sin v curvePar_z = cos(1) # z = 1*cos(v) # Project Curves/paths onto 3D Kime-surface Plot plot_ly(x = ~x1, y = ~y1, z = ~z1, type = 'surface', opacity=0.5, visible=T) %>% # Add Curve 1 add_trace(x=~curveS_x, y=~curveS_y, z=~curveS_z, type="scatter3d", mode="lines", line=list(width=20), name="Curve 1") %>% # Add Curve 2 add_trace(x=~curveP_x, y=~curveP_y, z=~curveP_z, type="scatter3d", mode="lines", line=list(width=20), name="Curve 2") %>% # Add Meridian add_trace(x=~curveMer_x, y=~curveMer_y, z=~curveMer_z, type="scatter3d", mode="lines", line=list(width=20), name="Meridian") %>% # Add Parallel add_trace(x=~curvePar_x, y=~curvePar_y, z=~curvePar_z, type="scatter3d", mode="lines", line=list(width=20), name="Meridian") %>% layout(title = "Spherical Surface Example", showlegend = FALSE) ``` The coefficients of the first fundamental form are computed as inner-products of the first-order partial derivatives with respect to $u$ and $v$: $${\begin{aligned} \frac{\partial X}{\partial u}=X_{u}&={\begin{bmatrix}-\sin u\sin v\\\cos u\sin v\\0\end{bmatrix}},\\ \frac{\partial X}{\partial v}=X_{v}&={\begin{bmatrix}\cos u\cos v\\\sin u\cos v\\-\sin v\end{bmatrix}}.\end{aligned}}$$ Hence, $${\begin{aligned} E&=\langle X_{u}| X_{u}\rangle=\sin ^{2}v\\F&=\langle X_{u} | X_{v}\rangle=0\\G&=\langle X_{v} | X_{v}\rangle=1\end{aligned}}.$$ And the first fundamental form is: $$\mathbb{I}\equiv {\begin{bmatrix}E&F\\F&G\end{bmatrix}}={\begin{bmatrix}\sin ^{2}v&0\\0&1\end{bmatrix}}.$$ For the same general parametric surface $X(u,v)$, the *second fundamental form* $\mathbb{II}$ is defined in terms of the linearly independent first order partial derivatives, $X_u$ and $X_v$, and the field of surface unitary normal vectors $n\equiv \frac{X_u \times X_v}{||X_u \times X_v ||}$: $$\mathbb{II} =L\ du^{2}+2M\ du\ dv + N\ dv^{2}\ = {\begin{bmatrix}du & dv \end{bmatrix}}\ {\begin{bmatrix}L&M\\M&N\end{bmatrix}}{\begin{bmatrix}du\\ dv\end{bmatrix}}\ ,$$ where the three coefficients $L, M, N$ are defined at a given point represent the projections of the second partial derivatives of $X$ at that point onto the normal line to $S$: $$L=\langle {X}_{uu} | \mathbf {n} \rangle \equiv -\langle {X}_{u} | {\mathbf {n}}_u \rangle \quad M=\langle X_{uv} | \mathbf {n} \rangle \equiv -\langle {X}_{u} | {\mathbf {n}}_v \rangle \quad N=\langle X_{vv} | \mathbf {n} \rangle\ \equiv -\langle {X}_{v} | {\mathbf {n}}_v \rangle.$$ The equivalences above are true since the normal $\mathbf {n}$ is perpendicular to the surface tangent plane containing the tangent vectors ${X}_{u},{X}_{v}$. Thus, $\langle {X}_{u} | \mathbf {n} \rangle = 0 =\langle {X}_{v} | \mathbf {n} \rangle$, and $$\langle {X}_{uu} | \mathbf {n} \rangle + \langle {X}_{u} | {\mathbf {n}}_u \rangle=0,$$ $$\langle {X}_{vv} | \mathbf {n} \rangle + \langle {X}_{v} | {\mathbf {n}}_v \rangle=0,$$ and $$\langle {X}_{uv} | \mathbf {n} \rangle + \langle {X}_{u} | {\mathbf {n}}_v \rangle=0.$$ In the spherical example above, $X(u,v)=(\cos u\sin v,\sin u\sin v,\cos v)'$, $R=1$, the second order derivatives are: $${\begin{aligned} \frac{\partial^2 X}{\partial u^2}=X_{uu}&={\begin{bmatrix}-\cos u\sin v\\ -\sin u\sin v\\0\end{bmatrix}},\\ \frac{\partial^2 X}{\partial u \partial v}=X_{uv}& ={\begin{bmatrix}-\sin u\cos v\\ -\cos u\cos v\\0\end{bmatrix}},\\ \frac{\partial^2 X}{\partial v^2}=X_{vv}&={\begin{bmatrix}-\cos u\sin v\\ -\sin u\sin v\\-\cos v\end{bmatrix}}.\end{aligned}}$$ For a sphere of radius $R$, the second fundamental form is the product of $R$ times the first fundamental form, $\mathbb{II} = R\times \mathbb{I}$. This is because the normal ${\bf{n}}=\frac{1}{R}X$, and thus, $\frac{\partial \bf{n}}{\partial u}\equiv n_u=\frac{1}{R} X_u$ and $\frac{\partial \bf{n}}{\partial v}\equiv n_v=\frac{1}{R} X_v$. This yields that $L=\langle {X}_{uu} | \mathbf {n} \rangle$ Thus, we can compute the second fundamental form components: $$\mathbb{II}\equiv {\begin{bmatrix}L&M\\M&N\end{bmatrix}}= {\begin{bmatrix} \overbrace{\langle {X}_{uu} | \mathbf {n} \rangle}^{-\langle {X}_{u} | {\mathbf {n}}_u \rangle} & \overbrace{\langle {X}_{uv} | \mathbf {n} \rangle}^{-\langle {X}_{u} | {\mathbf {n}}_v \rangle} \\ \\ \underbrace{\langle {X}_{uv} | \mathbf {n} \rangle}_{-\langle {X}_{u} | {\mathbf {n}}_v \rangle} & \underbrace{\langle {X}_{vv} | \mathbf {n} \rangle}_{-\langle {X}_{v} | {\mathbf {n}}_v \rangle} \end{bmatrix}} = \frac{1}{R} \mathbb{I}\equiv \frac{1}{R}{\begin{bmatrix}E&F\\F&G\end{bmatrix}}= \frac{1}{R}{\begin{bmatrix}\sin ^{2}v&0\\0&1\end{bmatrix}} .$$ The Gaussian curvature of a surface embedded in $\mathbb{R}^3$ is the ratio of the determinants of the second and first fundamental forms: $$\frac {\det(\mathbb {II} )}{\det(\mathbb {I} )}={\frac {LN-M^{2}}{EG-F^{2}}}.$$ In the spherical example, as the radius $R=1$, the surface Gaussian curvature is constant, $\frac {\det(\mathbb {II} )}{\det(\mathbb {I} )}= 1$. The [Brioschi formula](https://en.wikipedia.org/wiki/Gaussian_curvature#Alternative_formulas) expresses the Gaussian curvature purely in terms of the first fundamental form: $$K={\frac {{\begin{vmatrix}-{\frac {1}{2}}E_{vv}+F_{uv}-{\frac {1}{2}}G_{uu}&{\frac {1}{2}}E_{u}&F_{u}-{\frac {1}{2}}E_{v}\\F_{v}-{\frac {1}{2}}G_{u}&E&F\\{\frac {1}{2}}G_{v}&F&G\end{vmatrix}}-{\begin{vmatrix}0&{\frac {1}{2}}E_{v}&{\frac {1}{2}}G_{u}\\{\frac {1}{2}}E_{v}&E&F\\{\frac {1}{2}}G_{u}&F&G\end{vmatrix}}}{\left(EG-F^{2}\right)^{2}}}.$$ For any orthogonal parametrization, including polar coordinates in complex-time, the mixed term $F = 0$ and the Gaussian curvature is simplified to: $$K=-{\frac {1}{2{\sqrt {EG}}}}\left({\frac {\partial }{\partial u}}{\frac {G_{u}}{\sqrt {EG}}}+{\frac {\partial }{\partial v}}{\frac {E_{v}}{\sqrt {EG}}}\right) .$$ In the case of explicit surface definition via a function $z = J(x,y)$, the Gaussian curvature is: $$K={\frac {\langle J_{xx} | J_{yy}\rangle -J_{xy}^{2}}{\left(1+J_{x}^{2}+J_{y}^{2}\right)^{2}}}.$$ This differential geometric framework allows us to represent, model, understand, track, and analyze system dynamics with respect to complex time (kime). The basic idea is that the kime-surface metric tensor and different types of curvature measures provide mechanisms to capture *change with respect to kime magnitude and kime phase*. Also [Gauss' Theorema Egregium](https://en.wikipedia.org/wiki/Theorema_egregium) indicates that the surface Gaussian curvature can be computed by measuring length along surface geodesics via the first fundamental form and its first and second order partial derivatives. Furthermore, the Gaussian curvature is determined by the intrinsic surface metric tensor irrespective of the ambient space the surface is embedded in. Hence, the Gaussian curvature is an *intrinsic invariant* that is preserved under isometric deformations of the surface.