diff --git a/Manifest.toml b/Manifest.toml index 3fcfa38..ebf07d9 100644 --- a/Manifest.toml +++ b/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" diff --git a/Project.toml b/Project.toml index 61a728c..00fb2d9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "BLUEs" uuid = "b3a7f272-e305-45d1-bcf3-14d22bb67726" authors = ["G Jake Gebbie "] -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" diff --git a/src/convolutions.jl b/src/convolutions.jl new file mode 100644 index 0000000..d69cd54 --- /dev/null +++ b/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 diff --git a/src/dimensional_data.jl b/src/dimensional_data.jl index 50d4e50..c81cad2 100644 --- a/src/dimensional_data.jl +++ b/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")