code: import string
  
  import numpy as np
  import matplotlib.pyplot as plt
  
  from VertexTissue import globals as const
  from VertexTissue.Analysis import *
  
  from VertexTissue.funcs import euclidean_distance
  
  from periodic_2D import forces, omegas
  
  
  if __name__ == '__main__':
  
      square_length = lambda  G, t  : euclidean_distance(G.nodes[0]['pos'], G.nodes[2]['pos'])
  
      fig, axs = plt.subplots(2,int(np.ceil(len(omegas)*len(forces)/2)))
      fig.set_size_inches(12, 8)
      axs = axs.flatten()
  
      patterns=[]
      i=0
  
      ke=const.mu_apical
      kv = ke*60
      eta = const.eta
      kappa = eta + 2*kv
  
      alphabet_string = string.ascii_lowercase
      alphabet_list = list(alphabet_string)
  
  
      linewidth=3
      for f in forces:
          for omega in omegas:
  
  
  
  
  
              ax=axs[i]
              lbl = alphabet_list[i]
              i+=1
  
              arg = (kappa*ke/(omega*kv**2)+omega*eta/ke)/2
              delta = -np.arctan(arg)
  
              num = ke**2+(omega*kv)**2
              denom = (kappa*omega*ke)**2+(kv*eta*omega**2)**2
  
  
              denom2 = (kappa*ke)**3+kappa*ke*(omega*kv*eta)**2
      
              res = analyze_network_evolution(path='./data/viscoelastic/',
                                      pattern=f'periodic_force_{f}_freq_{omega}_*.pickle',
                                      func=square_length)
  
              res=np.array(res)
              t=res[:,0]
              ax.plot(t, res[:,1],linewidth=linewidth, label='numerical')
              A=f
              lam = 2*ke/eta + 1/60
              gamma = ke/eta * (2*A)
              B=(2.0/(lam*eta))*(0+gamma/lam)
              sol = (3.4+B)+gamma*(1/const.mu_apical-2.0/(const.eta*lam))*t - B*np.exp(-lam*t)
  
              
  
              num2 = -kv*2*A*omega*ke*eta*kv**2
              l_final = const.l_apical + 2*A/(2*omega*kv+omega*eta)
              l_trans = -2*np.exp(-lam*t)*(num2)/denom2
              amp = 2*A*np.sqrt(num/denom)
              sol = l_final+amp*np.sin(omega*t+delta) +l_trans
  
  
              ax.plot(t, sol, '--',linewidth=linewidth, label='theoretical')
              
              ax.set_xlabel('t (seconds)', fontsize=14)
              ax.set_ylabel('length', fontsize=14)
              ax.title.set_text(f'({lbl}) Force={f}, max error = {np.max(np.abs(sol-res[:,1])):.3e}')
              ax.title.set_fontsize(16)
              ax.legend(loc='right')
  
      
      
     
  
  
  
  
  
  
      
      plt.tight_layout()
      plt.show() import os
  
  import networkx as nx
  import numpy as np
  
  
  
  
  
  
  
  def validate(n,attr='pos'): 
      try:
          G=nx.read_gpickle(f't{n}.pickle')
          G1=nx.read_gpickle(f't_fast{n}.pickle')
  
          vals = nx.get_node_attributes(G,attr)
          vals1 = nx.get_node_attributes(G1,attr)
  
          max_diff  = np.max([np.max(np.abs(v-v1)) for v,v1 in zip(vals.values(),vals1.values()) ])
          print(f'{n}: {max_diff}')
  
  
          return True
      except:
          return False
  
  
  for i in range(300):
      cont= validate(i)
      if not cont:
          break
apis: numpy.sqrt matplotlib.pyplot.tight_layout numpy.arctan numpy.abs numpy.exp matplotlib.pyplot.show numpy.array numpy.sin
apis: numpy.abs