Checkexp

Project.toml

@@ -1,7 +1,7 @@
 name = "NonlinearNormalForm"
 uuid = "05e19671-dec8-4f15-984f-54eaa6ca64be"
 authors = ["Matt Signorelli"]
-version = "0.3.3"
+version = "0.4.0"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
@@ -26,7 +26,7 @@ Reexport = "1"
 ReferenceFrameRotations = "3"
 SkewLinearAlgebra = "1"
 StaticArrays = "1"
-TPSAInterface = "0.2.1"
+TPSAInterface = "0.2.2"
 julia = "1.9"
 
 [extras]

src/NonlinearNormalForm.jl

@@ -43,7 +43,6 @@ export        TaylorMap,
               VectorField,
 
               norm,
-              to_SO3,
 
               compose,
 
@@ -60,27 +59,18 @@ export        TaylorMap,
               normalize_eigenmode!,
               locate_modes!,
               moveback_unstable!,
-
-              normal!,
+
               normal,
-              gofix!,
-              gofix,
-              linear_a!,
-              linear_a,
               c_map,
               ci_map,
               c_jacobian,
               ci_jacobian,
 
               log!,
-              pb,
-
-
 
-              inv_with_log,
               equilibrium_moments,
               factorize,
-              fast_canonize,
+              canonize,
               compute_de_moivre,
               compute_sagan_rubin
 

src/de_moivre.jl

@@ -1,45 +1,45 @@
-function compute_de_moivre(a1::DAMap{V0}, ::Val{linear}=Val{false}()) where {V0<:StaticArray,linear}
+function compute_de_moivre(a1::DAMap{V0}, ::Val{linear}=Val{true}()) where {V0<:StaticArray,linear}
   # jp_mat[i] in FPP is J matrix restricted to i-th plane
   # ip_mat[i] in FPP is identity matrix restricted to i-th plane
   # jt_mat in FPP is symplectic s matrix
-  nv = isodd(length(V0)) ? length(V0)+1 : length(V0)
+  nhv = nhvars(a1)
   if linear
-    let a1_mat = jacobian(a1, VARS_CPARAM), a1i_mat = inv(a1_mat)
-      B = StaticArrays.sacollect(SVector{Int(nv/2),typeof(a1_mat)}, a1_mat*jp_mat(a1, i)*a1i_mat for i in 1:Int(nv/2))
-      K = StaticArrays.sacollect(SVector{Int(nv/2),typeof(a1_mat)}, -j_mat(a1)*B[i] for i in 1:Int(nv/2))
-      E = StaticArrays.sacollect(SVector{Int(nv/2),typeof(a1_mat)}, -B[i]*j_mat(a1) for i in 1:Int(nv/2))
-      H = StaticArrays.sacollect(SVector{Int(nv/2),typeof(a1_mat)}, a1_mat*ip_mat(a1, i)*a1i_mat for i in 1:Int(nv/2))
+    let a1_mat = jacobian(a1, HVARS), a1i_mat = inv(a1_mat)
+      B = StaticArrays.sacollect(SVector{Int(nhv/2),typeof(a1_mat)}, a1_mat*jp_mat(a1, i)*a1i_mat for i in 1:Int(nhv/2))
+      K = StaticArrays.sacollect(SVector{Int(nhv/2),typeof(a1_mat)}, -j_mat(a1)*B[i] for i in 1:Int(nhv/2))
+      E = StaticArrays.sacollect(SVector{Int(nhv/2),typeof(a1_mat)}, -B[i]*j_mat(a1) for i in 1:Int(nhv/2))
+      H = StaticArrays.sacollect(SVector{Int(nhv/2),typeof(a1_mat)}, a1_mat*ip_mat(a1, i)*a1i_mat for i in 1:Int(nhv/2))
       return (; H=H, B=B, E=E, K=K)
     end
   else
     let
       tmp = zero(a1)
       a1i = inv(a1)
-      B = StaticArrays.sacollect(SVector{Int(nv/2)}, begin
+      B = StaticArrays.sacollect(SVector{Int(nhv/2)}, begin
         setray!(tmp.v, v_matrix=jp_mat(a1, i))
         a1∘tmp∘a1i
-      end for i in 1:Int(nv/2)
+      end for i in 1:Int(nhv/2)
       )
       setray!(tmp.v, v_matrix=-j_mat(a1))
-      K = StaticArrays.sacollect(SVector{Int(nv/2)}, begin
+      K = StaticArrays.sacollect(SVector{Int(nhv/2)}, begin
         tmp∘B[i]
-      end for i in 1:Int(nv/2)
+      end for i in 1:Int(nhv/2)
       )
-      E = StaticArrays.sacollect(SVector{Int(nv/2)}, begin
+      E = StaticArrays.sacollect(SVector{Int(nhv/2)}, begin
         B[i]∘tmp
-      end for i in 1:Int(nv/2)
+      end for i in 1:Int(nhv/2)
       )
-      H = StaticArrays.sacollect(SVector{Int(nv/2)}, begin
+      H = StaticArrays.sacollect(SVector{Int(nhv/2)}, begin
         setray!(tmp.v, v_matrix=ip_mat(a1, i))
         a1∘tmp∘a1i
-      end for i in 1:Int(nv/2)
+      end for i in 1:Int(nhv/2)
       )
       return (; H=H, B=B, E=E, K=K)
     end
   end
 end
 
-function compute_de_moivre(a1::DAMap, ::Val{linear}=Val{false}()) where {linear}
+function compute_de_moivre(a1::DAMap, ::Val{linear}=Val{true}()) where {linear}
   # jp_mat[i] in FPP is J matrix restricted to i-th plane
   # ip_mat[i] in FPP is identity matrix restricted to i-th plane
   # jt_mat in FPP is symplectic s matrix
@@ -83,8 +83,8 @@ end
 # =================================================================================== #
 # Construct special matrices for lattice functions
 function ip_mat(m::DAMap{V0}, i) where {V0<:StaticArray}
-  nv = isodd(length(V0)) ? length(V0)+1 : length(V0)
-  return StaticArrays.sacollect(SMatrix{nv,nv,real(eltype(V0))}, 
+  nhv = nhvars(m) #isodd(length(V0)) ? length(V0)+1 : length(V0)
+  return StaticArrays.sacollect(SMatrix{nhv,nhv,real(eltype(V0))}, 
   begin
     if col == 2*i-1 && row == 2*i-1
       1
@@ -93,50 +93,42 @@ function ip_mat(m::DAMap{V0}, i) where {V0<:StaticArray}
     else
       0
     end
-  end for col in 1:nv for row in 1:nv)
+  end for col in 1:nhv for row in 1:nhv)
 end
 
 function ip_mat(m::DAMap, i)
-  nv = nvars(m)
-  if isodd(nv) 
-    nv += 1
-  end
+  nhv = nhvars(m)
   ip = zeros(real(eltype(m.v0)), nhv, nhv)
   ip[2*i-1, 2*i-1] = 1
   ip[2*i, 2*i] = 1
   return ip
 end
 
 function jp_mat(m::DAMap{V0}, i) where {V0<:StaticArray}
-  nv = isodd(length(V0)) ? length(V0)+1 : length(V0)
-  coast = isodd(length(V0))
-  coast_plane = coast && i == Int(nv/2)
-  return StaticArrays.sacollect(SMatrix{nv,nv,real(eltype(V0))}, 
+  nhv = nhvars(m) 
+  return StaticArrays.sacollect(SMatrix{nhv,nhv,real(eltype(V0))}, 
   begin
     if col == 2*i && row == 2*i-1
       1
-    elseif col == 2*i-1 && row == 2*i #&& !coast_plane
+    elseif col == 2*i-1 && row == 2*i
       -1
     else
       0
     end
-  end for col in 1:nv for row in 1:nv)
+  end for col in 1:nhv for row in 1:nhv)
 end
 
 function jp_mat(m::DAMap, i)
-  nv = nvars(m)
-  if isodd(nv) 
-    nv += 1
-  end
-  jp = zeros(real(eltype(m.v0)), nv, nv)
+  nhv = nhvars(m)
+  jp = zeros(real(eltype(m.v0)), nhv, nhv)
   jp[2*i-1, 2*i] = 1
   jp[2*i, 2*i-1] = -1
   return jp
 end
 
 function j_mat(m::DAMap{V0}) where {V0<:StaticArray}
-  nv = isodd(length(V0)) ? length(V0)+1 : length(V0)
-  return StaticArrays.sacollect(SMatrix{nv,nv,real(eltype(V0))}, 
+  nhv = nhvars(m) 
+  return StaticArrays.sacollect(SMatrix{nhv,nhv,real(eltype(V0))}, 
   begin
     if fld1(col,2) != fld1(row,2) 
       0
@@ -149,16 +141,13 @@ function j_mat(m::DAMap{V0}) where {V0<:StaticArray}
         0
       end
     end
-  end for col in 1:nv for row in 1:nv)
+  end for col in 1:nhv for row in 1:nhv)
 end
 
 function j_mat(m::DAMap)
-  nv = nvars(m)
-  if isodd(nv) 
-    nv += 1
-  end
-  j = zeros(real(eltype(m.v0)), nv, nv)
-  for i in 1:Int(nv/2)
+  nhv = nhvars(m)
+  j = zeros(real(eltype(m.v0)), nhv, nhv)
+  for i in 1:Int(nhv/2)
     j[2*i-1, 2*i] = 1
     j[2*i, 2*i-1] = -1
   end

src/map.jl

@@ -265,10 +265,11 @@ end
 
 # zero but new type potentially
 function zero(::Type{$t{V0,V,Q,S}}, m::TaylorMap) where {V0,V,Q,S}
-  return _zero($t{V0,V,Q,S}, getinit(eltype(V)), nvars(m), m)
+  return _zero($t{V0,V,Q,S}, getinit(m), nvars(m), m)
 end
 
-zero(::Type{$t}, m::TaylorMap{V0,V,Q,S}) where {V0,V,Q,S} = zero($t{V0,V,Q,S}, m)
+# Keep same initializer if not specified
+zero(::Type{$t}, m::TaylorMap{V0,V,Q,S}) where {V0,V,Q,S} = _zero($t{V0,V,Q,S}, getinit(m), nvars(m), m)
 
 # Copy ctor including optional TPSA init change
 function $t(m::TaylorMap; init::AbstractTPSAInit=getinit(m))
@@ -347,11 +348,6 @@ function $t(;
     np = ndiffs(init) - nv
   end
 
-
-  #println(nv)
-  #println(np)
-
-  #println(np)
   nn = nv+np
   # Check if nv+np agrees with ndiffs in init
   nn == ndiffs(init) || error("Number of variables + parameters does not agree with the number of differentials in the TPSA")