symmetric problems post processing routines finished

tri2full() and buildRES() functions were added as a complement in vectfit3.py module. These functions are useful when processing fitting results of symmetric matrix functions.

Author: SebastianL18

README.md

@@ -41,7 +41,7 @@ The vector fitting function has many options that can be modified via _opts_ dic
 Default options of vector fitting are already defined into vectfit3.py module. There, _opts_ dictionary contains the following configuration:
 
     opts={
-        "lowert_mat" : False, # F(s) samples belong to a full matrix
+        "symm_mat"   : False, # F(s) samples belong to a full matrix
         "relax"      : True,  # Use vector fitting with relaxed non triviality
         "stable"     : True,  # Enforce stable poles
         "asymp"      : 2,     # Include only D in fitting (not E).
@@ -135,21 +135,24 @@ In the following figures vector fitting results for some test cases are shown:
  * Differences between _vectfit3.py_ and the original MATLAB implementation are listed below:
 
     - All options for _vectfit3_ configuration are defined as boolean variables, except asymp which has 3 possible states.
-    - The new option "lowert_mat" for _vectfit3_ configuration is added. This indicates when $F(s)$ samples belong to a lower triangular matrix function, that reduces the number of elements to fit for a symmetric matrix function.
+    - The new option "symm_mat" for _vectfit3_ configuration is added. This indicates when $F(s)$ samples belong to a symmetric matrix function in CMO. This is a common practice because it reduces the number of elements to fit. See test 4
     - A new method to sort the poles computed during the identification process is implemented.
     - Real and complex data is meticulously treated and differentiated through the entire process.
     - General code organization.
     - A new method to compute error plots is implemented. The new error is $log_{10}(error_{relative})$
     - Error graphs are now plotted outside the magnitude axis as a subplot in the same figure.
+    - Now "cmplx_ss" and "symm_mat" flags are also members of SER "dictionary". This is helpful when using the space-state model in external post-processing routines.
+    - ss2pr() subroutine is renamed as buildRES() and is modified to return just the residues matrixes, because poles are already given by vectfit() main function.
+    - tri2full subroutine, which is used to transform compressed results of a symmetric problem, is now internally called from buildRES() prior to compute residue matrixes. However can be imported from the module if it is needed.
 
 ### In development
 Nowadays, some final details in _vectfit3.py_ are still in progress:
-* The function _"buildPOLRES()"_, to build the poles and residues model from the state-space model generated by vector fitting needs to be finished.
-* Cases for symmetric problems that are represented with lower triangular matrixes need to be reconstructed to obtain a full state-space o pole-residue model.
+* _"buildRES()"_ function that computes residues matrixes from SER, only considers symmetric data reduction. Therefore a subroutine to map the results from element-wise representation "vectfit default" to a full matrix representation, when original data belong to a non-symmetric matrix function must be created. This function should take into account the mapping method used, for instance Column Major Order (CMO) and Row Major Order (RMO).
 
 To contribute, give suggestions or report any bug please contact me:
  * Sebastian Loaiza Elejalde
     - _Dsc Student in Power Systems_
+    - [_CINVESTAV - GDL_](https://unidad.gdl.cinvestav.mx/)
     - sebloel18@gmail.com 📬
     - sebastian.loaiza@cinvestav.mx 📬
 

vectfit3.py

@@ -64,7 +64,7 @@
 # Dictionary which contains the default settings fot vectfit. 
 # Any key can be modified to change som options included in vectfit3
 opts={
-    "lowert_mat" : False, # Indicates when F(s) samples belong to a lower triangular matrix (symmetric problem)
+    "symm_mat"   : False, # Indicates when F(s) samples belong to a lower triangular matrix (symmetric problem)
     "relax"      : True,  # Use vector fitting with relaxed non triviality
     "stable"     : True,  # Enforce stable poles
     "asymp"      : 2,     # Include only D in fitting (not E). See [4]
@@ -280,7 +280,7 @@ def vectfitPlot(F,fit,s,opts,initialState=False):
     plt.show()
 
 # vectfit3() subroutine.
-def buildSER(Ac,Br,Cc,Dr,Er,cmplx_ss,lowert_mat=False):
+def buildSER(Ac,Br,Cc,Dr,Er,cmplx_ss,symm_mat):
     """Function to build the state-space model of the fitted function such as:
         F(s) = C * (sI-A)^-1 * B + D + sE
     
@@ -292,13 +292,13 @@ def buildSER(Ac,Br,Cc,Dr,Er,cmplx_ss,lowert_mat=False):
          - D: Constant terms matrix. Real matrix of shape [Nc x 1]
          - E: Proportional terms matrix. Real matrix of shape [Nc x 1]
          - cmplx_ss: Boolean option from opts{} which indicates if a complex system is needed
-         - lower_mat: Boolean option from opts{} which indicates if data comes from a symmetric problem
+         - symm_mat: Boolean option from opts{} which indicates if data comes from a symmetric matrix problem
         
         Results.
         
          - SER: Dictionary that storages A,B,C,D,E system matrixes adapted as indicated by cmplx_ss 
     """
-    SER=dict(A=Ac,B=Br,C=Cc,D=Dr,E=Er,cmplxType=True,symmData=lowert_mat) #complex state-space system
+    SER=dict(A=Ac,B=Br,C=Cc,D=Dr,E=Er,cmplx_type=True,symm_mat=symm_mat) #complex state-space system
     if not(cmplx_ss):
         # Real state-space system is required so the matrixes are modified
         Ar=np.real(Ac) #importing real poles to Ar
@@ -325,7 +325,7 @@ def buildSER(Ac,Br,Cc,Dr,Er,cmplx_ss,lowert_mat=False):
         SER["A"]=Ar
         SER["B"]=Br
         SER["C"]=Cr
-        SER["cmplxType"]=False
+        SER["cmplx_type"]=False
     return SER
 
 # build_RES() subroutine.
@@ -344,24 +344,22 @@ def tri2full(SER, real2cmplx=False):
     C=SER["C"]
     D=SER["D"]
     E=SER["E"]
-    if real2cmplx and A.dtype==np.float64: # Real to complex transformation is required
+    if real2cmplx and not(SER["cmplx_type"]): # Real to complex transformation is required
         # Complex versions of A and C matrixes:
         Ac=np.zeros(A.shape, dtype=np.complex128)
         Cc=np.zeros(C.shape, dtype=np.complex128)
-        # For complex poles B is a vector of ones:
-        B=np.ones(B.shape, dtype=np.float64)
-        for m in range(A.shape[0]-1):
-            if A[m,m+1]!=0: #case for complex poles
-                Ac[m,m]=A[m,m]+1j*A[m,m+1]          #reference pole 
-                Ac[m+1,m+1]=A[m+1,m+1]-1j*A[m,m+1]  #conjugated pole
-                #ERROR BUILDING C complex
-                Cc[:,m]=C[:,m]+1j*C[:,m+1]          #Reference C value
-                Cc[:,m+1]=np.conj(Cc[:,m])          #Conjugated C value
-            else: #case pure real poles 
+        for m in range(A.shape[0]):
+            if B[m,0]==1: #case for pure real poles
                 Ac[m,m]=A[m,m]
                 Cc[:,m]=C[:,m]
-        #Updating A and C with their complex forms:
+            elif B[m,0]==2: #case for complex poles 
+                Ac[m,m]=A[m,m]+1j*A[m,m+1]          #reference pole 
+                Ac[m+1,m+1]=A[m+1,m+1]+1j*A[m+1,m]  #conjugated pair
+                Cc[:,m]=C[:,m]+1j*C[:,m+1]
+                Cc[:,m+1]=np.conj(Cc[:,m])
+        #Updating arrays with their complex forms:
         A=Ac
+        B=np.ones(B.shape, dtype=np.float64) #for complex poles B is a vector of ones
         C=Cc
     # Unzip process parameters:
     n=B.shape[0]    # Order of approximation
@@ -400,20 +398,34 @@ def tri2full(SER, real2cmplx=False):
     SER["C"]=Cf
     SER["D"]=Df
     SER["E"]=Ef
-    SER["symmData"]=False
+    SER["symm_mat"]=False
     return SER
 
 def buildRES(SER):
     """Function to generate the residues matrix of the fitted function computed by vectfit
-        *Returns a tuple with:
-            - R: residue matrixes stacked as a 3D array of shape [Ny x Ny x n]. Ny is the matrix function size and n the aproximation order
+        *Returns residues matrixes stacked in Res:
+            - Res: a 3D array of shape [Ny x Ny x n]. where Ny is the matrix function size and n the aproximation order
     """
     # state-space model matrixes in SER:
     n=SER["A"].shape[0]   # order of aproximation
-    if SER["symmData"] and not(SER["cmplxType"]):
+    if SER["symm_mat"]:
         # Data needs to be resized to a full matrix representation instead of lower trinagular and element-wise representation
         #and also the space-state model needs to be converte to a complex
         SER=tri2full(SER, real2cmplx=True)
+        C=SER["C"]
+        B=SER["B"]
+        Ny=C.shape[0]
+        Res=np.zeros((Ny,Ny,n), dtype=np.complex128) #Residues matrixes
+        Rk=np.zeros((Ny,Ny), dtype=np.complex128) #Instant Res values
+        i=0 # auxiliar index
+        for k in range(n):
+            Rk[:,:]=0
+            for m in range(Ny):
+                i=m*n+k
+                #in order to perform the correct matrix operation C and B values must be reshaped 
+                Rk+=np.reshape(C[:,i],(Ny,1))@np.reshape(B[i,:],(1,Ny))   
+            Res[:,:,k]=Rk
+    return Res
 
 # * ----------------------------------------------------------  main vectfit3 function ---------------------------------------------------------- *
 
@@ -833,7 +845,7 @@ def vectfit(F,s,poles,weights,opts=opts):
     C=SERC
     D=SERD
     E=SERE
-    SER=buildSER(A,B,C,D,E,opts["cmplx_ss"],opts["lowert_mat"])
+    SER=buildSER(A,B,C,D,E,opts["cmplx_ss"],opts["symm_mat"])
 
     # Vector fitting process finished.
     return (SER,poles,rmserr,fit)

vectfit_testing.py

@@ -47,6 +47,7 @@
     print(" v/ Fitting process completed. Aproximation error achieved = ",rmserr)
     print("\nFinal poles computed:\n",poles)
 
+    
 elif test==2:
     print("Test 2: 18th order frequency response F(s) of two dimentions") # -------------------------------------------------------------------- #
 
@@ -252,12 +253,13 @@
 
     # vector fitting configuration
     from vectfit3 import opts
-    opts["asymp"]=3        # Modified to include D and E in fitting
-    opts["logx"]=False     # Modified to use linear axis for x
-    opts["spy2"]=False     # Modified to omit graphs generation into the iterative application of vectfit
-    opts["phaseplot"]=True # Modified to include the phase angle graph in the results
-    opts["skip_res"]=True  # Modified to skip residue computation during the iterative execution of vector fitting
-    opts["cmplx_ss"]=False 
+    opts["asymp"]=3         # Modified to include D and E in fitting
+    opts["logx"]=False      # Modified to use linear axis for x
+    opts["spy2"]=False      # Modified to omit graphs generation into the iterative application of vectfit
+    opts["phaseplot"]=True  # Modified to include the phase angle graph in the results
+    opts["skip_res"]=True   # Modified to skip residue computation during the iterative execution of vector fitting
+    opts["symm_mat"]=True   # Modified to indicate that F(s) samples belong to the symmetric matrix Y(s)
+    opts["cmplx_ss"]=False  # Modified to create a real only space-state model
     # Remaining options by default
 
     print("\n * Applying 5 iterations of vector fitting...")
@@ -284,6 +286,9 @@
         print("     ...",itr+1," iterations applied")
     print(" v/ Fitting process completed. Aproximation error achieved = ",rmserr)
     print("\nFinal poles computed:\n",poles)
-    from vectfit3 import tri2full
-    SER=tri2full(SER, real2cmplx=True)
-    print("\nExpanded SER matrix = \n",SER["C"])
+
+    from vectfit3 import buildRES
+    Res=buildRES(SER)
+    print("\nResidues matrixes computed:\n")
+    for k in range(n):
+        print("\n",Res[:,:,k],"\n")