added zeros and ones for allocated TPSs

docs/src/man/o_all.md

@@ -4,11 +4,13 @@ The following functions from Base have been overloaded for operations with `TPS`
 +, -, *, /, ^, ∘, inv, atan, hypot, abs, sqrt, exp, log, 
 sin, cos, tan, csc, sec, cot, sinc, sinh, cosh, tanh, csch, 
 sech, coth, asin, acos, atan, acsc, asec, acot, asinh, acosh, 
-atanh, acsch, asech, acoth, zero, one, real, imag, conj, angle, 
-complex, promote_rule, getindex, setindex!, ==, <, >, <=, >=, 
-!=, isequal, show, copy!
+atanh, acsch, asech, acoth, zero, zeros, one, ones, real, imag, 
+conj, angle, complex, promote_rule, getindex, setindex!, ==, <, 
+>, <=, >=, !=, isequal, show, copy!
 ```
 
+`zeros` and `ones` are overloaded from Base so that allocated `TPS`/`ComplexTPS`s are placed in each element. If we didn't explicity overload these functions, every element would correspond to the exact same heap-allocated TPS, which is problematic when setting individual monomial coefficients of the same TPS.
+
 `GTPSA.jl` overloads (and exports) the following functions from the corresponding packages:
  **`LinearAlgebra`**: `norm` 
 **`SpecialFunctions`**: `erf`, `erfc`

src/GTPSA.jl

@@ -40,7 +40,9 @@ import Base:  +,
               asech ,
               acoth ,
               zero  ,
+              zeros ,
               one   ,
+              ones  ,
               real  ,
               imag  ,
               conj  ,

src/getset.jl

@@ -356,7 +356,7 @@ function gradient(t::Union{TPS,ComplexTPS}; include_params=false)
 end
 
 """
-    jacobian!(result::Matrix{<:Union{Float64,ComplexF64}}, m::Vector{<:Union{TPS,ComplexTPS}}; include_params=false)
+    jacobian!(result::Matrix{<:Union{Float64,ComplexF64}}, m::Union{Vector{<:Union{TPS,ComplexTPS}}, SubArray{<:Union{TPS,ComplexTPS},N,P,I,true}}; include_params=false)  where {N,P,I}
 
 Extracts the first-order partial derivatives (evaluated at 0) from the Vector of TPSs. 
 and fills the `result` matrix in-place. The partial derivatives wrt the parameters will 
@@ -371,7 +371,7 @@ in the TPSs.
 ### Output
 - `result`         -- Preallocated matrix to fill with the Jacobian of `m`
 """
-function jacobian!(result::Matrix{<:Union{Float64,ComplexF64}}, m::Vector{<:Union{TPS,ComplexTPS}}; include_params=false)
+function jacobian!(result::Matrix{<:Union{Float64,ComplexF64}}, m::Union{Vector{<:Union{TPS,ComplexTPS}}, SubArray{<:Union{TPS,ComplexTPS},N,P,I,true}}; include_params=false) where {N,P,I}
   n = numvars(first(m))
   if include_params
     n += numparams(first(m))
@@ -390,7 +390,7 @@ function jacobian!(result::Matrix{<:Union{Float64,ComplexF64}}, m::Vector{<:Unio
 end
 
 """
-    jacobian(m::Vector{<:Union{TPS,ComplexTPS}}; include_params=false)
+    jacobian(m::Union{Vector{<:Union{TPS,ComplexTPS}}, SubArray{<:Union{TPS,ComplexTPS},N,P,I,true}}; include_params=false) where {N,P,I}
 
 Extracts the first-order partial derivatives (evaluated at 0) from the Vector of TPSs. 
 The partial derivatives wrt the parameters will also be extracted when the `include_params` 
@@ -404,7 +404,7 @@ the first-order monomial coefficients already in the TPSs.
 ### Output
 - `J`              -- Jacobian of `m`
 """
-function jacobian(m::Vector{<:Union{TPS,ComplexTPS}}; include_params=false)
+function jacobian(m::Union{Vector{<:Union{TPS,ComplexTPS}}, SubArray{<:Union{TPS,ComplexTPS},N,P,I,true}}; include_params=false) where {N,P,I}
   n = numvars(first(m))
   if include_params
     n += numparams(first(m))
@@ -419,7 +419,7 @@ function jacobian(m::Vector{<:Union{TPS,ComplexTPS}}; include_params=false)
 end
 
 """
-    jacobiant!(result::Matrix{<:Union{Float64,ComplexF64}}, m::Vector{<:Union{TPS,ComplexTPS}}; include_params=false)
+    jacobiant!(result::Matrix{<:Union{Float64,ComplexF64}}, m::Union{Vector{<:Union{TPS,ComplexTPS}}, SubArray{<:Union{TPS,ComplexTPS},N,P,I,true}}; include_params=false) where {N,P,I}
 
 Extracts the first-order partial derivatives (evaluated at 0) from the Vector of TPSs, 
 as the transpose of the Jacobian. Because of Julia's column-major indexing vs. C's row-major, 
@@ -435,7 +435,7 @@ first-order monomial coefficients already in the TPSs and filling `result`.
 ### Output
 - `result`         -- Preallocated matrix to fill with the transpose of the Jacobian of `m`
 """
-function jacobiant!(result::Matrix{<:Union{Float64,ComplexF64}}, m::Vector{<:Union{TPS,ComplexTPS}}; include_params=false)
+function jacobiant!(result::Matrix{<:Union{Float64,ComplexF64}}, m::Union{Vector{<:Union{TPS,ComplexTPS}}, SubArray{<:Union{TPS,ComplexTPS},N,P,I,true}}; include_params=false) where {N,P,I}
   n = numvars(first(m))
   if include_params
     n += numparams(first(m))
@@ -449,7 +449,7 @@ function jacobiant!(result::Matrix{<:Union{Float64,ComplexF64}}, m::Vector{<:Uni
 end
 
 """
-    jacobiant(m::Vector{<:Union{TPS,ComplexTPS}}; include_params=false)
+    jacobiant(m::Union{Vector{<:Union{TPS,ComplexTPS}}, SubArray{<:Union{TPS,ComplexTPS},N,P,I,true}}; include_params=false) where {N,P,I}
 
 Extracts the first-order partial derivatives (evaluated at 0) from the Vector of TPSs, 
 as the transpose of the Jacobian. Because of Julia's column-major indexing vs. C's row-major, 
@@ -465,7 +465,7 @@ first-order monomial coefficients already in the TPSs.
 ### Output
 - `Jt`             -- Transpose of the Jacobian of `m`
 """
-function jacobiant(m::Vector{<:Union{TPS,ComplexTPS}}; include_params=false)
+function jacobiant(m::Union{Vector{<:Union{TPS,ComplexTPS}}, SubArray{<:Union{TPS,ComplexTPS},N,P,I,true}}; include_params=false) where {N,P,I}
   n = numvars(first(m))
   if include_params
     n += numparams(first(m))

src/operators.jl

@@ -29,6 +29,29 @@ function one(ct1::ComplexTPS)::ComplexTPS
   return ct
 end
 
+# --- zeros and ones (taken from Base.array.jl) --- 
+# We overload this because we want each element of the array to be a separate allocated TPS
+for (fname, felt) in ((:zeros, :zero), (:ones, :one))
+  @eval begin
+      $fname(::Type{T}, dims::Base.DimOrInd...) where {T<:Union{TPS,ComplexTPS}} = $fname(T, dims)
+      $fname(::Type{T}, dims::NTuple{N, Union{Integer, Base.OneTo}}) where {T<:Union{TPS,ComplexTPS},N} = $fname(T, map(to_dim, dims))
+      function $fname(::Type{T}, dims::NTuple{N, Integer}) where {T<:Union{TPS,ComplexTPS},N}
+          a = Array{T,N}(undef, dims)
+          for idx in eachindex(a)
+            a[idx] = $felt(T)
+          end
+          return a
+      end
+      function $fname(::Type{T}, dims::Tuple{}) where {T<:Union{TPS,ComplexTPS}}
+          a = Array{T}(undef)
+          for idx in eachindex(a)
+            a[idx] = $felt(T)
+          end
+          return a
+      end
+  end
+end
+
 # --- Unary ---
 # TPS:
 function +(t1::TPS)::TPS