11import jax
2-
3- jax .config .update ("jax_enable_x64" , True )
4- from time import time
5-
62import jax .numpy as jnp
73import matplotlib .pyplot as plt
4+ import numpy as np
85from scipy .constants import speed_of_light
96
107from simphony .simulation .jax_tools import python_based_while_loop
118
9+ jax .config .update ("jax_enable_x64" , True )
10+
1211
1312# @jax.jit
1413def _initial_poles (model_order , frequency , alpha ):
@@ -42,18 +41,13 @@ def _lstsq_matrices(model_order, transfer_function, phi0, phi1):
4241 M = jnp .zeros (((num_ports ** 2 ) * (model_order ), (model_order )), dtype = complex )
4342 B = jnp .zeros (((num_ports ** 2 ) * (model_order )), dtype = complex )
4443
45- A1 = phi0
46- Q1 , R11 = jnp .linalg .qr (A1 )
47-
4844 iter = 0
4945 for i in range (num_ports ):
5046 for j in range (num_ports ):
5147 D = jnp .diag (transfer_function [:, i , j ])
5248 A_block = jnp .hstack ([phi0 , - D @ phi1 ]) # never build the big matrix
5349 Q , R = jnp .linalg .qr (A_block , mode = "reduced" )
5450
55- R11 = R [: model_order + 1 , : model_order + 1 ]
56- R12 = R [: model_order + 1 , model_order + 1 :]
5751 R22 = R [model_order + 1 :, model_order + 1 :]
5852 Q2 = Q [:, model_order + 1 :]
5953
@@ -80,9 +74,6 @@ def _lstsq_matrices(model_order, transfer_function, phi0, phi1):
8074 return M , B
8175
8276
83- import numpy as np
84-
85-
8677def _full_lstsq_matrices (transfer_function , phi0 , phi1 ):
8778 D = []
8879 V = []
@@ -395,9 +386,9 @@ def optimize_order(bias_fn, min_order, max_order):
395386 C_max_minus_1 , * _ = bias_fn (max_order - 1 )
396387 lambda_lower = jnp .abs (C_max_minus_1 - C_max )
397388 lambda_upper = C_min - C_max
398- l = jnp .log10 (lambda_lower )
399- u = jnp .log10 (lambda_upper )
400- complexity_penalty = 10 ** (0.5 * (u + l ))
389+ lower_log = jnp .log10 (lambda_lower )
390+ upper_log = jnp .log10 (lambda_upper )
391+ complexity_penalty = 10 ** (0.5 * (upper_log + lower_log ))
401392
402393 # TODO: implement Golden Section Search
403394 # to minimize C - complexity_penalty * order
@@ -490,10 +481,7 @@ def main():
490481 )
491482 residues = jnp .reshape (residues [1 :], (10 , 1 , 1 ))
492483 feedthrough = jnp .zeros ((1 , 1 ), dtype = complex )
493- N = len (poles )
494-
495484 f = jnp .linspace (0.001 , 10 / (2 * jnp .pi ), 100 )
496- aortic_response = pole_residue_response (f , poles , residues , feedthrough )
497485
498486 _mzi , info = sax .circuit (
499487 netlist = {
@@ -531,7 +519,6 @@ def mzi(wl=1.55):
531519 f_max = speed_of_light / 1.50e-6
532520 # f_min = speed_of_light / 1.565e-6
533521 # f_max = speed_of_light / 1.5350e-6
534- f_center = 0.5 * (f_min + f_max )
535522 frequency = jnp .linspace (f_min , f_max , 1000 )
536523
537524 plt .plot (
@@ -598,15 +585,12 @@ def main2():
598585 )
599586 residues = jnp .reshape (residues [1 :], (10 , 1 , 1 ))
600587 feedthrough = jnp .zeros ((1 , 1 ), dtype = complex )
601- N = len (poles )
602-
603588 f = jnp .linspace (0.001 , 10 / (2 * jnp .pi ), 100 )
604589 response = pole_residue_response (f , poles , residues , feedthrough )
605590
606591 poles1 , residues1 , feedthrough1 , error = vector_fitting (
607592 10 , response , f , max_iterations = 5
608593 )
609- toc = time ()
610594 H = pole_residue_response (f , poles1 , residues1 , feedthrough1 )
611595 plt .plot (f , jnp .abs (H [:, 0 , 0 ]) ** 2 )
612596 plt .plot (f , jnp .abs (response [:, 0 , 0 ]) ** 2 , "r--" )
0 commit comments