In [1]:
import networkx as nx
import numpy as np
from numpy import linalg as la
In [2]:
from sympy import init_printing
init_printing()
In [3]:
import matplotlib.pyplot as plt
%pylab inline
In [4]:
G = nx.DiGraph(selfloops=True)
We use a directed graph with various properties along the nodes and edges. The direction describes the propagation of signals in the system.
There are three kinds of nodes: inputs nodes, internal nodes, and output nodes. There is the same number of input and output nodes (say n). The number of internal nodes may be different. Each internal node has an associated matrix describing its relationship between its incoming and outgoing signals. It suffices for now to take $2 \times 2$ matrices of the form $\begin{pmatrix} t && -r \\ r && t \end{pmatrix}$ corresponding to a beamsplitter, where $r$ and $t$ are the reflectivity and transmissivity of the beamsplitter, respectively. These satisfy $r^2+t^2 = 1$.
In general we may want other matrices, but it's not really necessary.
If the signal along several edges is thought of as a vector, multiplying by the matrix from the left represents the signal traveling through the element. This formalism works only for linear networks.
Let's make an example graph:
In [5]:
rs = np.asarray([0.9,0.5,0.8]) ## some sample values
ts = np.sqrt(1.-rs**2) ## ts are determined from rs
In [6]:
N = 2 ## number of input nodes
In [7]:
for i in range(N): ## make the input and output nodes
G.add_node(i*2,label='x_in_'+str(i))
G.add_node(i*2+1,label='x_out_'+str(i))
for i,(r,t) in enumerate(zip(rs,ts)): ## make the remaining nodes
G.add_node(2*N+i,label='x_'+str(i),M=np.matrix([[t,-r],[r,t]]))
In [8]:
G.nodes(data=True) ## display the nodes
Out[8]:
In [24]:
num_nodes = len(G.nodes(data=True))
Each (directed) edge $j$ has a time delay $\tau_j$. In general a delay line may have an additional phase shift $\exp(i\theta_j)$ which is determined by a number $\theta_j$.
We will also include a pair of indices for each edge. The first index corresponds to the previous node and the second index corresponds to the next node. The indices indicate enumerations of the edges with respect to the input and output nodes, respectively. If the previous or next node is an input or output node of the graph, the index will be $0$.
For now, let's assume that only internal edges have nonzero delays.
For the visualization, it would be nice if for a given node, the incoming and outgoing edges with the same index value would appear as a straight line, since this physically means the signal is being transmitted without reflecting.
In [9]:
## edges to inputs
G.add_edge(0,4,delay=0.,indices=(0,0),theta=0.,edge_type = 'input',edge_num=0)
G.add_edge(2,6,delay=0.,indices=(0,1),theta=0.,edge_type = 'input',edge_num=1)
## edges to outputs
G.add_edge(4,1,delay=0.,indices=(1,0),theta=0.,edge_type = 'output',edge_num=2)
G.add_edge(6,3,delay=0.,indices=(0,0),theta=0.,edge_type = 'output',edge_num=3)
## internal edges
G.add_edge(4,5,delay=1.,indices=(0,0),theta=0.,edge_type = 'internal',edge_num=4)
G.add_edge(5,4,delay=1.,indices=(1,1),theta=0.,edge_type = 'internal',edge_num=5)
G.add_edge(5,6,delay=1.,indices=(0,0),theta=0.,edge_type = 'internal',edge_num=6)
G.add_edge(6,5,delay=1.,indices=(1,1),theta=0.,edge_type = 'internal',edge_num=7)
In [10]:
G.edges(data=True)
Out[10]:
In [13]:
## I can make a diagram for the graph, output to file
A=nx.to_agraph(G)
A.draw('file.ps',prog='neato')
In [21]:
internal_edges = {(edge[0],edge[1]):edge[2] for edge in G.edges(data=True) if edge[2]['edge_type'] == 'internal'}
m = len(internal_edges)
In [25]:
# input_edges = [edge for edge in G.edges(data=True) if edge[2]['edge_type'] == 'input']
# output_edges = [edge for edge in G.edges(data=True) if edge[2]['edge_type'] == 'output']
In [26]:
M1 = np.zeros((m,m))
internal_node_range = range(2*N,num_nodes)
internal_connections = []
for i in internal_node_range: ## internal nodes
outgoing_connections = nx.edges(G,[i])
internal_connections += [connection for connection in outgoing_connections if connection[1] in internal_node_range]
for i in internal_connections:
for j in internal_connections:
if i[1] == j[0]:
matrix_indices = G.edge[i[0]][i[1]]['indices'][0], G.edge[j[0]][j[1]]['indices'][1]
M1[internal_edges[j]['edge_num']-2*N,internal_edges[i]['edge_num']-2*N] = G.node[i[1]]['M'][matrix_indices]
In [27]:
M1
Out[27]:
In [28]:
all_connections = []
for i in range(num_nodes): ## internal nodes
outgoing_connections = nx.edges(G,[i])
all_connections += [connection for connection in outgoing_connections if connection[1] in range(num_nodes)]
In [29]:
all_edges = {(edge[0],edge[1]):edge[2] for edge in G.edges(data=True)}
m_all = len(all_edges)
In [30]:
U = np.zeros((m_all,m_all))
for i in all_connections:
for j in all_connections:
if i[1] == j[0]:
matrix_indices = G.edge[i[0]][i[1]]['indices'][0], G.edge[j[0]][j[1]]['indices'][1]
U[all_edges[j]['edge_num'],all_edges[i]['edge_num']] = G.node[i[1]]['M'][matrix_indices]
In [34]:
## should coincide with M1
M1 = U[4:8,4:8]
In [40]:
M4 = U[:4,:4]
In [ ]:
M3 = U[8:16,4:8]
In [36]:
M4 = U[8:16,8:16]
Using the run_Potapov function of this method generates the variables that will be used for the first part of the visualization. Those are contained in an instance of the Time_Delay_Network. Specifically, the outputs we will want to plot are (1) Time_Delay_Network.roots (2) Time_Delay_Network.spatial_modes.
The roots $r_1,...,r_n$ are a list of complex numbers corresponding to the modes indexed by $1,...,n$. The imaginary part of root $r_k$ correspond to the frequency of mode $k$, and the real part of $r_k$ indicate the decay coefficients of mode $k$.
The spatial_modes are a list $v_1,...,v_n$ of complex-valued vectors. Each vector $v_k$ in the list corresponds to a mode $k$, in the same order as the roots. Each vector has the same length as the number of time delays of the network, $\tau_1,...,\tau_m$. The $l_{th}$ component $v_{k,l}$ of vector $v_k$ indicates the spatially normalized amplitude of mode $k$ along the delay $\tau_l$.
What would be cool is to be able to select one or many modes $1,...,k,...,n$ and to illustrate the spatial component of the signal of the selected modes along the graph. Specifically, the frequency of the root could correspond to a color or a periodic sinusoidal shape (higher frequency would be more blue or a shorter period), or both. The absolute value of the spatial mode component could be indicated by the thickness of the signal along each time delay. A phase shift could be indicated by a shift in the frequency of a sinusoidal signal.
In [11]:
import Potapov_Code
In [225]:
Network = Potapov_Code.Time_Delay_Network.Example3() ## an example network with hardcoded values
In [227]:
Network.run_Potapov(commensurate_roots=True) ## run the analysis
In [238]:
roots = Network.roots ## roots
plt.scatter(map(lambda z: z.real, roots), map(lambda z: z.imag, roots))
Out[238]:
In [237]:
Network.spatial_modes ## the spatial modes
Out[237]: