organize and update convolutions

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.11.1"
 manifest_format = "2.0"
-project_hash = "7b213f9da7628bbc0245f8ee42e7bb7f64e5ca6c"
+project_hash = "23ef8814b2d710e68611e0eb07933c2a7150191a"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
@@ -20,9 +20,9 @@ version = "1.5.0"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "d80af0733c99ea80575f612813fa6aa71022d33a"
+git-tree-sha1 = "50c3c56a52972d78e8be9fd135bfb91c9574c140"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "4.1.0"
+version = "4.1.1"
 
     [deps.Adapt.extensions]
     AdaptStaticArraysExt = "StaticArrays"
@@ -49,15 +49,16 @@ version = "1.1.2"
 
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra"]
-git-tree-sha1 = "3640d077b6dafd64ceb8fd5c1ec76f7ca53bcf76"
+git-tree-sha1 = "d60a1922358aa203019b7857a2c8c37329b8736c"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.16.0"
+version = "7.17.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
     ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
     ArrayInterfaceCUDAExt = "CUDA"
     ArrayInterfaceCUDSSExt = "CUDSS"
+    ArrayInterfaceChainRulesCoreExt = "ChainRulesCore"
     ArrayInterfaceChainRulesExt = "ChainRules"
     ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
     ArrayInterfaceReverseDiffExt = "ReverseDiff"
@@ -71,6 +72,7 @@ version = "7.16.0"
     CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
     CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e"
     ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
     GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
     SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

Project.toml

@@ -1,7 +1,7 @@
 name = "BLUEs"
 uuid = "b3a7f272-e305-45d1-bcf3-14d22bb67726"
 authors = ["G Jake Gebbie <ggebbie@whoi.edu>"]
-version = "0.2.3"
+version = "0.2.4"
 
 [deps]
 AlgebraicArrays = "8af735f6-f3e5-4048-bdaa-40a2355e9eea"
@@ -16,7 +16,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulLinearAlgebra = "c14bd059-d406-4571-8f61-9bd20e53c30b"
 
 [compat]
-AlgebraicArrays = "1.0.3"
+#AlgebraicArrays = "1.0.6-DEV"
 DimensionalData = "0.29"
 Measurements = "2"
 Statistics = "1"

src/convolutions.jl

@@ -0,0 +1,146 @@
+
+"""
+function convolve(x::DimArray{T},E::AbstractDimArray) where T <: Number
+
+Take the convolution of E and x
+Account for proper overlap of dimension.
+Sum and take into account units.
+Return an `AbstractDimArray`
+"""
+function convolve(x::VectorArray,E::AbstractDimArray) 
+    tnow = last(first(rangedims(x)))
+    lags = first(dims(E))
+    vals = sum([E[ii,:] ⋅ x[Near(tnow-ll),:] for (ii,ll) in enumerate(lags)])
+    #(vals isa Number) ? (return DimArray([vals],first(dims(x)))) : (return DimArray(vals,first(dims(x))))
+    #(vals isa Number) ? (return vals) : (return DimArray(vals,first(dims(x))))
+    (vals isa Number) ? (return VectorArray(DimArray([vals],first(rangedims(x))))) : (return VectorArray(AlgebraicArray(vals,first(rangedims(x)))))
+end
+"""
+function convolve(x::AbstractDimArray, E::AbstractDimArray, coeffs::UnitfulMatrix}
+    the `coeffs` argument signifies that x is a 3D array (i.e. >1 state variables)
+
+    this function both convolves, and linearly combines the propagated state variables
+"""
+# function convolve(x::AbstractDimArray,E::AbstractDimArray, coeffs::UnitfulMatrix)
+#     statevars = x.dims[3]
+#     mat = UnitfulMatrix(transpose([convolve(x[:,:,At(s)], E) for s in statevars])) * coeffs
+#     return getindexqty(mat, 1,1)
+# end
+
+function convolve(x::VectorArray, M::AbstractDimArray, t::Number)
+    lags = first(dims(M))
+    return sum([M[ii,:] ⋅ x[Near(t-ll),:] for (ii,ll) in enumerate(lags)])
+end
+
+#coeffs signifies that x is 3D 
+function convolve(x::AbstractDimArray, E::AbstractDimArray, t::Number, coeffs::UnitfulMatrix)
+    statevars = x.dims[3]
+    mat = UnitfulMatrix(transpose([convolve(x[:,:,At(s)], E, t) for s in statevars]))*coeffs
+    return getindexqty(mat, 1,1) 
+end
+
+#don't handle the ndims(M) == 3 case here but I'll get back to it
+function convolve(x::VectorArray, M::AbstractDimArray, Tx::Union{Ti, Vector}, coeffs::UnitfulMatrix)
+    if ndims(M) == 2
+        return DimArray([convolve(x,M,Tx[tt], coeffs) for (tt, yy) in enumerate(Tx)], Tx)
+    elseif ndims(M) == 3
+        error("some code should go here")
+    else
+        error("M has wrong number of dims") 
+    end
+end
+
+#function convolve(x::AbstractDimArray,M::AbstractDimArray,Tx::Union{Ti,Vector})
+function convolve(x::VectorArray,M::AbstractDimArray,Tx::Union{Ti,Vector})
+    if ndims(M) == 2 
+        return VectorArray(DimArray([convolve(x,M,Tx[tt]) for (tt,yy) in enumerate(Tx)],Tx))
+    elseif ndims(M) == 3
+
+        # do a sample calculation to get units.
+        Msmall = M[:,:,1]
+        yunit = unit.(vec(convolve(x,Msmall,Tx))[1]) # assume everything has the same units
+
+        y = DimArray(zeros(length(Tx),last(size(M)))yunit,(Tx,last(dims(M))))
+        for (ii,vv) in enumerate(last(dims(M)))
+
+            Msmall = M[:,:,ii]
+            y[:,ii] = convolve(x,Msmall,Tx)
+
+        end
+        return y
+    else
+        error("M has wrong number of dims")
+    end
+end
+# basically repeats previous function: any way to simplify?
+function convolve(P::MatrixArray, M::AbstractDimArray, Tx::Union{Ti,Vector}) 
+    T2 = typeof(parent(convolve(first(P),M,Tx)))
+    Pyx = Array{T2}(undef,size(P))
+    for i in eachindex(P)
+        #Pyx[i] = parent(parent(convolve(P[i],M,Tx,coeffs)))
+        Pyx[i] = parent(convolve(P[i],M,Tx))
+    end
+    return MatrixArray(DimArray(Pyx,domaindims(P)))
+end
+
+function convolve(P::MatrixArray,M) 
+    #function convolve(P::DimArray{T},M) where T<: AbstractDimArray
+    # became more complicated when returning a scalar was not allowed
+    T2 = typeof(first(parent(convolve(first(P),M))))
+    Pyx = Array{T2}(undef,size(P))
+    for i in eachindex(P)
+        #Pyx[i] = convolve(P[i],M)
+        Pyx[i] = first(parent(convolve(P[i],M)))
+    end
+    return AlgebraicArray(Pyx,RowVector(["1"]),rangedims(P))
+end
+
+# basically repeats previous function: any way to simplify?
+function convolve(P::MatrixArray,M::AbstractDimArray,coeffs::DimVector)
+    T2 = typeof(first(parent(convolve(first(P),M,coeffs))))
+    #T2 = typeof(convolve(first(P),M,coeffs))
+    Pyx = Array{T2}(undef,size(P))
+    for i in eachindex(P)
+        #        Pyx[i] = convolve(P[i],M,coeffs)
+        Pyx[i] = first(parent(convolve(P[i],M,coeffs)))
+    end
+    #return DimArray(Pyx,dims(P))
+    return transpose(VectorArray(DimArray(Pyx,rangedims(P))))
+    #return AlgebraicArray(Pyx,RowVector(["1"]),rangedims(P))
+    #return MatrixArray(DimArray(Pyx,(RowVector(["1"]),rangedims(P))))
+end
+
+function convolve(x::VectorArray, M::AbstractDimArray, coeffs::DimVector) 
+    statevars = dims(x,3) # equal to rangedims(x)[3]
+    vals = sum([convolve(x[:,:,At(s)], M)  * coeffs[At(s)] for s in statevars])
+    #return sum([convolve(x[:,:,At(s)], M)  * coeffs[At(s)] for s in statevars])
+    (vals isa Number) ? (return VectorArray(DimArray([vals],first(rangedims(x))))) : (return VectorArray(AlgebraicArray(vals,first(rangedims(x)))))
+end
+
+function convolve(x::VectorArray, M::AbstractDimArray, Tx::Ti, coeffs::DimVector) # where T <: Number
+    if ndims(M) == 2
+        return VectorArray(DimArray([convolve(x, M, Tx[tt], coeffs) for tt in eachindex(Tx)], Tx))
+    elseif ndims(M) == 3
+        error("some code should go here")
+    else
+        error("M has wrong number of dims") 
+    end
+end
+# basically repeats previous function: any way to simplify?
+function convolve(P::MatrixArray, M::AbstractDimArray, Tx::Ti, coeffs::DimVector) 
+    #T2 = typeof(convolve(first(P),M,Tx,coeffs))
+    #T2 = typeof(first(parent(convolve(first(P),M,Tx,coeffs))))
+    #T2 = typeof(parent(parent(convolve(first(P),M,Tx,coeffs))))
+    T2 = typeof(parent(convolve(first(P),M,Tx,coeffs)))
+    Pyx = Array{T2}(undef,size(P))
+    for i in eachindex(P)
+        #Pyx[i] = parent(parent(convolve(P[i],M,Tx,coeffs)))
+        Pyx[i] = parent(convolve(P[i],M,Tx,coeffs))
+    end
+    return MatrixArray(DimArray(Pyx,domaindims(P)))
+end
+
+function convolve(x::VectorArray, M::AbstractDimArray, t::Number, coeffs::DimVector)
+    statevars = dims(x,3)
+    return sum([convolve(x[:,:,At(s)], M, t)  * coeffs[At(s)] for s in statevars])
+end

src/dimensional_data.jl

@@ -1,3 +1,4 @@
+
 ext = Base.get_extension(AlgebraicArrays, :AlgebraicArraysDimensionalDataExt)
 if !isnothing(ext)
     RowVector = ext.RowVector
@@ -50,141 +51,6 @@ end
 
 #standard_error(P::AbstractDimArray{T,2}) where T <: Number = DimArray(.√diag(P),first(dims(P)))
 
-"""
-function convolve(x::DimArray{T},E::AbstractDimArray) where T <: Number
-
-Take the convolution of E and x
-Account for proper overlap of dimension.
-Sum and take into account units.
-Return an `AbstractDimArray`
-"""
-function convolve(x::VectorArray,E::AbstractDimArray) 
-    tnow = last(first(rangedims(x)))
-    lags = first(dims(E))
-    vals = sum([E[ii,:] ⋅ x[Near(tnow-ll),:] for (ii,ll) in enumerate(lags)])
-    #(vals isa Number) ? (return DimArray([vals],first(dims(x)))) : (return DimArray(vals,first(dims(x))))
-    #(vals isa Number) ? (return vals) : (return DimArray(vals,first(dims(x))))
-    (vals isa Number) ? (return VectorArray(DimArray([vals],first(rangedims(x))))) : (return VectorArray(AlgebraicArray(vals,first(rangedims(x)))))
-end
-"""
-function convolve(x::AbstractDimArray, E::AbstractDimArray, coeffs::UnitfulMatrix}
-    the `coeffs` argument signifies that x is a 3D array (i.e. >1 state variables)
-
-    this function both convolves, and linearly combines the propagated state variables
-"""
-# function convolve(x::AbstractDimArray,E::AbstractDimArray, coeffs::UnitfulMatrix)
-#     statevars = x.dims[3]
-#     mat = UnitfulMatrix(transpose([convolve(x[:,:,At(s)], E) for s in statevars])) * coeffs
-#     return getindexqty(mat, 1,1)
-# end
-
-function convolve(x::VectorArray, M::AbstractDimArray, t::Number)
-    lags = first(dims(M))
-    return sum([M[ii,:] ⋅ x[Near(t-ll),:] for (ii,ll) in enumerate(lags)])
-end
-
-#coeffs signifies that x is 3D 
-function convolve(x::AbstractDimArray, E::AbstractDimArray, t::Number, coeffs::UnitfulMatrix)
-    statevars = x.dims[3]
-    mat = UnitfulMatrix(transpose([convolve(x[:,:,At(s)], E, t) for s in statevars]))*coeffs
-    return getindexqty(mat, 1,1) 
-end
-
-#don't handle the ndims(M) == 3 case here but I'll get back to it
-function convolve(x::VectorArray, M::AbstractDimArray, Tx::Union{Ti, Vector}, coeffs::UnitfulMatrix)
-    if ndims(M) == 2
-        return DimArray([convolve(x,M,Tx[tt], coeffs) for (tt, yy) in enumerate(Tx)], Tx)
-    elseif ndims(M) == 3
-        error("some code should go here")
-    else
-        error("M has wrong number of dims") 
-    end
-end
-
-function convolve(x::AbstractDimArray,M::AbstractDimArray,Tx::Union{Ti,Vector})
-    if ndims(M) == 2 
-        return DimArray([convolve(x,M,Tx[tt]) for (tt,yy) in enumerate(Tx)],Tx)
-    elseif ndims(M) == 3
-
-        # do a sample calculation to get units.
-        Msmall = M[:,:,1]
-        yunit = unit.(vec(convolve(x,Msmall,Tx))[1]) # assume everything has the same units
-
-        y = DimArray(zeros(length(Tx),last(size(M)))yunit,(Tx,last(dims(M))))
-        for (ii,vv) in enumerate(last(dims(M)))
-
-            Msmall = M[:,:,ii]
-            y[:,ii] = convolve(x,Msmall,Tx)
-
-        end
-        return y
-    else
-        error("M has wrong number of dims")
-    end
-end
-
-function convolve(P::MatrixArray,M) 
-    #function convolve(P::DimArray{T},M) where T<: AbstractDimArray
-    # became more complicated when returning a scalar was not allowed
-    T2 = typeof(first(parent(convolve(first(P),M))))
-    Pyx = Array{T2}(undef,size(P))
-    for i in eachindex(P)
-        #Pyx[i] = convolve(P[i],M)
-        Pyx[i] = first(parent(convolve(P[i],M)))
-    end
-    return AlgebraicArray(Pyx,RowVector(["1"]),rangedims(P))
-end
-
-# basically repeats previous function: any way to simplify?
-function convolve(P::MatrixArray,M::AbstractDimArray,coeffs::DimVector)
-    T2 = typeof(first(parent(convolve(first(P),M,coeffs))))
-    #T2 = typeof(convolve(first(P),M,coeffs))
-    Pyx = Array{T2}(undef,size(P))
-    for i in eachindex(P)
-        #        Pyx[i] = convolve(P[i],M,coeffs)
-        Pyx[i] = first(parent(convolve(P[i],M,coeffs)))
-    end
-    #return DimArray(Pyx,dims(P))
-    return transpose(VectorArray(DimArray(Pyx,rangedims(P))))
-    #return AlgebraicArray(Pyx,RowVector(["1"]),rangedims(P))
-    #return MatrixArray(DimArray(Pyx,(RowVector(["1"]),rangedims(P))))
-end
-
-function convolve(x::VectorArray, M::AbstractDimArray, coeffs::DimVector) 
-    statevars = dims(x,3) # equal to rangedims(x)[3]
-    vals = sum([convolve(x[:,:,At(s)], M)  * coeffs[At(s)] for s in statevars])
-    #return sum([convolve(x[:,:,At(s)], M)  * coeffs[At(s)] for s in statevars])
-    (vals isa Number) ? (return VectorArray(DimArray([vals],first(rangedims(x))))) : (return VectorArray(AlgebraicArray(vals,first(rangedims(x)))))
-end
-
-function convolve(x::VectorArray, M::AbstractDimArray, Tx::Ti, coeffs::DimVector) # where T <: Number
-    if ndims(M) == 2
-        return VectorArray(DimArray([convolve(x, M, Tx[tt], coeffs) for tt in eachindex(Tx)], Tx))
-    elseif ndims(M) == 3
-        error("some code should go here")
-    else
-        error("M has wrong number of dims") 
-    end
-end
-# basically repeats previous function: any way to simplify?
-function convolve(P::MatrixArray, M::AbstractDimArray, Tx::Ti, coeffs::DimVector) 
-    #T2 = typeof(convolve(first(P),M,Tx,coeffs))
-    #T2 = typeof(first(parent(convolve(first(P),M,Tx,coeffs))))
-    #T2 = typeof(parent(parent(convolve(first(P),M,Tx,coeffs))))
-    T2 = typeof(parent(convolve(first(P),M,Tx,coeffs)))
-    Pyx = Array{T2}(undef,size(P))
-    for i in eachindex(P)
-        #Pyx[i] = parent(parent(convolve(P[i],M,Tx,coeffs)))
-        Pyx[i] = parent(convolve(P[i],M,Tx,coeffs))
-    end
-    return MatrixArray(DimArray(Pyx,domaindims(P)))
-end
-
-function convolve(x::VectorArray, M::AbstractDimArray, t::Number, coeffs::DimVector)
-    statevars = dims(x,3)
-    return sum([convolve(x[:,:,At(s)], M, t)  * coeffs[At(s)] for s in statevars])
-end
-
 # function uncertainty_units(x::DimArray{T}) where T<: Number
 
 #     unitlist = unit.(x)
@@ -247,3 +113,5 @@ function combine(x0::Estimate,y1::Estimate,E1::Function)
     P = x0.P - dP
     return Estimate(v,P)
 end
+
+include("convolutions.jl")