From 04778ab314747ae08dfca0cb093c4bf2d1a07489 Mon Sep 17 00:00:00 2001 From: Justin Ackers Date: Fri, 12 May 2023 08:29:07 +0200 Subject: [PATCH 01/23] Updated interface of TransferFunction --- src/TransferFunction.jl | 44 ++++++++++++++++++---------------------- test/TransferFunction.jl | 6 +++--- 2 files changed, 23 insertions(+), 27 deletions(-) diff --git a/src/TransferFunction.jl b/src/TransferFunction.jl index 2abed86f..a628ff4a 100644 --- a/src/TransferFunction.jl +++ b/src/TransferFunction.jl @@ -3,15 +3,19 @@ export TransferFunction, sampleTF, setTF mutable struct TransferFunction freq::Vector{Float64} data::Matrix{ComplexF64} + interpolator::Vector{AbstractInterpolation} inductionFactor::Vector{Float64} function TransferFunction(freq_, datain::Array{T}, inductionFactor=ones(size(datain, 2))) where {T<:Complex} - freq = freq_[1]:(freq_[2]-freq_[1]):freq_[end] + if length(freq_) != size(datain,1); error("Length of frequency axis ($(length(freq_))) does not match the number of samples in the data ($(size(datain,1)))!") end data = reshape(deepcopy(datain), size(datain,1), size(datain, 2)) - return new(freq_, data, inductionFactor) + interpolator = [extrapolate(interpolate((freq_,), channel, Gridded(Linear())), Interpolations.Flat()) for channel in eachcol(data)] + return new(freq_, data, interpolator, inductionFactor) end end +Base.show(io::IO, ::MIME"text/plain", tf::TransferFunction) = print(io, "MPIFiles.TransferFunction: \n\t$(size(tf.data,2)) channels\n\t$(size(tf.data,1)) frequency samples from $(tf.freq[1]) Hz to $(tf.freq[end]) Hz") + function TransferFunction(freq_, ampdata, phasedata, args...) data = ampdata.*exp.(im.*phasedata) return TransferFunction(freq_, data, args...) @@ -19,7 +23,7 @@ end function TransferFunction(filename::String; kargs...) filenamebase, ext = splitext(filename) - if ext == ".h5" + if ext == ".h5" tf = load_tf(filename) else #if ext == "s1p" || ext == "s2p" tf = load_tf_fromVNA(filename; kargs...) @@ -27,29 +31,21 @@ function TransferFunction(filename::String; kargs...) return tf end -function getindex(tmf::TransferFunction, x::UnitRange, chan::Integer=1) - a = tmf.data[x,chan] - return a +function getindex(tmf::TransferFunction, args...) + try + return getindex(tmf.data, args...) + catch e + @warn "The indexing using square brackets on TransferFunction objects now always operates on the integer indizes of the underlying transfer function data. To use frequency interpolation, use tf(freq, channel) instead of tf[[freq],channel]." + rethrow(e) + end end -#function getindex(tmf::TransferFunction, x::Real, chan::Integer=1) -# a = tmf.interp[chan](x) -# return a -#end - -function getindex(tmf::TransferFunction, X::Union{Vector,AbstractRange}, chan::Integer=1) - I = extrapolate(interpolate((tmf.freq,), tmf.data[:,chan], Gridded(Linear())), Interpolations.Flat()) - - return [I(x) for x in X] +function (tmf::TransferFunction)(x, chan::Integer=1) + if chan>length(tmf.interpolator); error("The TransferFunction only has $(length(tmf.interpolator)) channel(s), unable to access channel $(chan)") end + return tmf.interpolator[chan](x) end -function getindex(tmf::TransferFunction, X::Union{Vector,AbstractRange}, chan::AbstractRange) - out = zeros(ComplexF64, length(X), length(chan)) - for d=1:length(chan) - out[:,d] = tmf[X,d] - end - return out -end +(tmf::TransferFunction)(x, chan::AbstractArray) = hcat([tmf(x,c) for c in chan]...) function load_tf(filename::String) tf = h5read(filename,"/transferFunction") @@ -59,6 +55,7 @@ function load_tf(filename::String) end function combine(tf1, tf2) + if tf1.freq != tf2.freq; error("The frequency axes of the transfer functions do not match. Can not combine!") end freq = tf1.freq data = cat(tf1.data,tf2.data, dims=2) inductionFactor = cat(tf1.inductionFactor, tf2.inductionFactor, dims=1) @@ -126,8 +123,7 @@ end function sampleTF(tmf::TransferFunction, f::MPIFile) freq = rxFrequencies(f) numChan = rxNumChannels(f) - numFreq = length(freq) - return tmf[freq,1:numChan] + return tmf(freq,1:numChan) end diff --git a/test/TransferFunction.jl b/test/TransferFunction.jl index e3c2fab8..1d4c2d98 100644 --- a/test/TransferFunction.jl +++ b/test/TransferFunction.jl @@ -15,13 +15,13 @@ tfh5path = joinpath(tmpdir,"transferFunction","example.h5") @test a.freq == b.freq @test a.data == b.data - @test a[[1],1] == a[[0.0],1] + @test a[1,1] == a(0.0,1) c = MPIFiles.combine(MPIFiles.combine(a,a),a) rm(tfh5path) MPIFiles.save(tfh5path, c) - @test a[[1],1] == c[[1],2] - @test a[[1],1] == c[[1],3] + @test a[1,1] == c[1,2] + @test a[1,1] == c[1,3] measBruker = MPIFile(fnMeasBruker) tf = sampleTF(c, measBruker) From 659c8e9e578e50e2641f9633afbd6ee702a3f1e3 Mon Sep 17 00:00:00 2001 From: Justin Ackers Date: Tue, 16 May 2023 11:52:19 +0200 Subject: [PATCH 02/23] add units to TransferFunction --- Project.toml | 7 ++++--- src/MPIFiles.jl | 1 + src/TransferFunction.jl | 39 ++++++++++++++++++++++++++++++--------- 3 files changed, 35 insertions(+), 12 deletions(-) diff --git a/Project.toml b/Project.toml index ec7bf2ec..09e8f554 100644 --- a/Project.toml +++ b/Project.toml @@ -12,8 +12,8 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" Graphics = "a2bd30eb-e257-5431-a919-1863eab51364" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" -ImageMetadata = "bc367c6b-8a6b-528e-b4bd-a4b897500b49" ImageAxes = "2803e5a7-5153-5ecf-9a86-9b4c37f5f5ac" +ImageMetadata = "bc367c6b-8a6b-528e-b4bd-a4b897500b49" Inflate = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" @@ -33,15 +33,17 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" UnitfulAngles = "6fb2a4bd-7999-5318-a3b2-8ad61056cd98" +UnitfulParsableString = "06c00241-927a-4d5b-bb5e-6b5a2ada3567" [compat] AxisArrays = "0.3, 0.4" CodecZlib = "0.7" +DocStringExtensions = "0.8, 0.9" FileIO = "1.0" Graphics = "0.4, 1.0" HDF5 = "0.14, 0.15, 0.16" -ImageMetadata = "0.9" ImageAxes = "0.6" +ImageMetadata = "0.9" Inflate = "0.1.2" Interpolations = "0.12, 0.13, 0.14" Reexport = "1.0" @@ -49,7 +51,6 @@ SparsityOperators = "0.3, 0.4" StableRNGs = "1.0.0" Unitful = "0.17, 1.2" UnitfulAngles = "0.6.1" -DocStringExtensions = "0.8, 0.9" julia = "1.5" [extras] diff --git a/src/MPIFiles.jl b/src/MPIFiles.jl index e3b7351b..0350060c 100644 --- a/src/MPIFiles.jl +++ b/src/MPIFiles.jl @@ -26,6 +26,7 @@ using Pkg.Artifacts using Unitful using InteractiveUtils using UnitfulAngles +using UnitfulParsableString using Inflate, SHA using StableRNGs using REPL: fielddoc diff --git a/src/TransferFunction.jl b/src/TransferFunction.jl index a628ff4a..bc2205b1 100644 --- a/src/TransferFunction.jl +++ b/src/TransferFunction.jl @@ -5,20 +5,33 @@ mutable struct TransferFunction data::Matrix{ComplexF64} interpolator::Vector{AbstractInterpolation} inductionFactor::Vector{Float64} + units::Vector{Unitful.FreeUnits} - function TransferFunction(freq_, datain::Array{T}, inductionFactor=ones(size(datain, 2))) where {T<:Complex} + function TransferFunction(freq_::Vector{<:Real}, datain::Array{<:Complex}; inductionFactor::Vector{<:Real}=ones(size(datain, 2)), units::Vector=Unitful.FreeUnits[Unitful.NoUnits for i in 1:size(datain, 2)]) + + parsed_units = Unitful.FreeUnits[] + for tmp in units + if isa(tmp, String); tmp = uparse(tmp) end # get correct unit from unit strings + if !isa(tmp, Unitful.FreeUnits); tmp = Unitful.unit(tmp) end # get correct unit for numbers and quantities, e.g. unit(1) or unit(1u"V/T"), includes the case that the string parsing returned one of these types + push!(parsed_units, tmp) + end + if length(freq_) != size(datain,1); error("Length of frequency axis ($(length(freq_))) does not match the number of samples in the data ($(size(datain,1)))!") end + if size(datain, 2) != length(inductionFactor); error("Number of channels in data ($(size(datain, 2))) does not match the number of channels in the given inductionFactor ($(size(inductionFactor,1)))!") end + if size(datain, 2) != length(units); error("Number of channels in data ($(size(datain, 2))) does not match the number of channels in the given units ($(size(units,1)))!") end + data = reshape(deepcopy(datain), size(datain,1), size(datain, 2)) interpolator = [extrapolate(interpolate((freq_,), channel, Gridded(Linear())), Interpolations.Flat()) for channel in eachcol(data)] - return new(freq_, data, interpolator, inductionFactor) + return new(freq_, data, interpolator, inductionFactor, parsed_units) end end -Base.show(io::IO, ::MIME"text/plain", tf::TransferFunction) = print(io, "MPIFiles.TransferFunction: \n\t$(size(tf.data,2)) channels\n\t$(size(tf.data,1)) frequency samples from $(tf.freq[1]) Hz to $(tf.freq[end]) Hz") +Base.show(io::IO, ::MIME"text/plain", tf::TransferFunction) = print(io, "MPIFiles.TransferFunction: \n\t$(size(tf.data,2)) channels, units of $(string.(tf.units))\n\t$(size(tf.data,1)) frequency samples from $(tf.freq[1]) Hz to $(tf.freq[end]) Hz") -function TransferFunction(freq_, ampdata, phasedata, args...) +function TransferFunction(freq_::Vector{<:Real}, ampdata::Array{<:Real,N}, phasedata::Array{<:Real,N}; kwargs...) where N + if size(ampdata) != size(phasedata); error("The size of ampdata and phasedata must match!") end data = ampdata.*exp.(im.*phasedata) - return TransferFunction(freq_, data, args...) + return TransferFunction(freq_, data; kwargs...) end function TransferFunction(filename::String; kargs...) @@ -42,7 +55,7 @@ end function (tmf::TransferFunction)(x, chan::Integer=1) if chan>length(tmf.interpolator); error("The TransferFunction only has $(length(tmf.interpolator)) channel(s), unable to access channel $(chan)") end - return tmf.interpolator[chan](x) + return tmf.interpolator[chan](x) .* tmf.units[chan] end (tmf::TransferFunction)(x, chan::AbstractArray) = hcat([tmf(x,c) for c in chan]...) @@ -51,21 +64,29 @@ function load_tf(filename::String) tf = h5read(filename,"/transferFunction") freq = h5read(filename,"/frequencies") inductionFactor = h5read(filename,"/inductionFactor") - return TransferFunction(freq,tf,inductionFactor) + unit = [] + try + unit = h5read(filename, "/unit") + catch # if h5read fails, it should mean that there is no units, maybe do something better here + return TransferFunction(freq, tf, inductionFactor=inductionFactor) + end + return TransferFunction(freq, tf, inductionFactor=inductionFactor, units=uparse.(unit)) end -function combine(tf1, tf2) +function combine(tf1::TransferFunction, tf2::TransferFunction) if tf1.freq != tf2.freq; error("The frequency axes of the transfer functions do not match. Can not combine!") end freq = tf1.freq data = cat(tf1.data,tf2.data, dims=2) inductionFactor = cat(tf1.inductionFactor, tf2.inductionFactor, dims=1) - return TransferFunction(freq, data, inductionFactor) + units = cat(tf1.units, tf2.units, dims=1) + return TransferFunction(freq, data, inductionFactor=inductionFactor, units=units) end function save(filename::String, tf::TransferFunction) h5write(filename, "/transferFunction", tf.data) h5write(filename, "/frequencies", tf.freq) h5write(filename, "/inductionFactor", tf.inductionFactor) + h5write(filename, "/unit", string.(tf.units)) return nothing end From bc00873dc167b09971ed66193963fd09bf388364 Mon Sep 17 00:00:00 2001 From: "Justin Ackers @ mpi05" Date: Wed, 5 Jul 2023 17:29:24 +0200 Subject: [PATCH 03/23] nicer print --- src/TransferFunction.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TransferFunction.jl b/src/TransferFunction.jl index bc2205b1..19fc501e 100644 --- a/src/TransferFunction.jl +++ b/src/TransferFunction.jl @@ -26,7 +26,7 @@ mutable struct TransferFunction end end -Base.show(io::IO, ::MIME"text/plain", tf::TransferFunction) = print(io, "MPIFiles.TransferFunction: \n\t$(size(tf.data,2)) channels, units of $(string.(tf.units))\n\t$(size(tf.data,1)) frequency samples from $(tf.freq[1]) Hz to $(tf.freq[end]) Hz") +Base.show(io::IO, ::MIME"text/plain", tf::TransferFunction) = print(io, "MPIFiles.TransferFunction: \n\t$(size(tf.data,2)) channel(s), units of $(string.(tf.units))\n\t$(size(tf.data,1)) frequency samples from $(tf.freq[1]) Hz to $(tf.freq[end]) Hz") function TransferFunction(freq_::Vector{<:Real}, ampdata::Array{<:Real,N}, phasedata::Array{<:Real,N}; kwargs...) where N if size(ampdata) != size(phasedata); error("The size of ampdata and phasedata must match!") end From 8bd96ca19859fb0257b041361637e29e58256bb0 Mon Sep 17 00:00:00 2001 From: nHackel Date: Tue, 26 Sep 2023 16:13:40 +0200 Subject: [PATCH 04/23] Allow concrete subtypes to work with concrete parametric types --- src/Utils.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Utils.jl b/src/Utils.jl index 992f6b3b..76afe41c 100644 --- a/src/Utils.jl +++ b/src/Utils.jl @@ -101,10 +101,10 @@ macro keyoptional(expr) end export concreteSubtypes -function concreteSubtypes(type::DataType) +function concreteSubtypes(type::Union{DataType, UnionAll}) subtypes_ = subtypes(type) # Only add filtered ones but check all subtypes - allSubtypes = filter(x -> !isabstracttype(x), subtypes_) + allSubtypes = filter(x -> isconcretetype(x) || isstructtype(x), subtypes_) for subtype in subtypes_ subsubtypes_ = concreteSubtypes(subtype) allSubtypes = vcat(allSubtypes, subsubtypes_) From 0c7e3e150b16b3eaa0b9bb4d0b9ba5c9d7f1b067 Mon Sep 17 00:00:00 2001 From: Tobias Knopp Date: Tue, 26 Sep 2023 23:52:34 +0200 Subject: [PATCH 05/23] remove reexport altogether and in turn bump version --- Project.toml | 4 +--- src/MPIFiles.jl | 27 +++++++++++---------------- 2 files changed, 12 insertions(+), 19 deletions(-) diff --git a/Project.toml b/Project.toml index 163dab2d..f61eba96 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MPIFiles" uuid = "371237a9-e6c1-5201-9adb-3d8cfa78fa9f" authors = ["Tobias Knopp "] -version = "0.13.0" +version = "0.14.0" [deps] AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9" @@ -23,7 +23,6 @@ Mmap = "a63ad114-7e13-5084-954f-fe012c677804" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" LinearOperatorCollection = "a4a2c56f-fead-462a-a3ab-85921a5f2575" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" @@ -46,7 +45,6 @@ ImageMetadata = "0.9" ImageAxes = "0.6" Inflate = "0.1.2" Interpolations = "0.12, 0.13, 0.14" -Reexport = "1.0" LinearOperatorCollection = "1.1" StableRNGs = "1.0.0" Unitful = "0.17, 1.2" diff --git a/src/MPIFiles.jl b/src/MPIFiles.jl index 9b0df446..5456cd59 100644 --- a/src/MPIFiles.jl +++ b/src/MPIFiles.jl @@ -3,22 +3,21 @@ module MPIFiles using UUIDs using Graphics: @mustimplement -using Reexport using LinearOperatorCollection using FFTW -@reexport using AxisArrays +using AxisArrays const axes = Base.axes -@reexport using Interpolations -@reexport using HDF5 -@reexport using Dates -@reexport using DelimitedFiles -@reexport using ImageMetadata +using Interpolations +using HDF5 +using Dates +using DelimitedFiles +using ImageMetadata using ImageAxes -@reexport using LinearAlgebra -@reexport using Random -@reexport using Mmap -@reexport using Statistics -@reexport using Unitful +using LinearAlgebra +using Random +using Mmap +using Statistics +using Unitful using CodecZlib using Tar using Pkg.PlatformEngines @@ -32,10 +31,6 @@ using StableRNGs using REPL: fielddoc using DocStringExtensions -if VERSION < v"1.1" - isnothing(x) = x == nothing -end - ### global import list ### import Base: convert, get, getindex, haskey, iterate, length, ndims, range, read, show, time, write, close From 3bd3f3937a37d028fa31570fefb0a94d5db88df8 Mon Sep 17 00:00:00 2001 From: Jonas Schumacher Date: Wed, 27 Sep 2023 08:44:43 +0200 Subject: [PATCH 06/23] Add Aqua and fix arising issues (#67) * Add Aqua and fix arising issues * Re-add re-export for Unitful * Comment in all tests * Update Project.toml --------- Co-authored-by: Tobias Knopp --- Project.toml | 19 +++++++++++-------- README.md | 3 ++- docs/src/lowlevel.md | 2 +- src/Conversion.jl | 2 +- src/DatasetStore/BrukerDatasetStore.jl | 5 +++-- src/DatasetStore/MDFDatasetStore.jl | 5 +++-- src/MPIFiles.jl | 2 +- src/Measurements.jl | 2 +- src/Positions/Positions.jl | 2 +- test/runtests.jl | 5 +++++ 10 files changed, 29 insertions(+), 18 deletions(-) diff --git a/Project.toml b/Project.toml index f61eba96..cd7854a0 100644 --- a/Project.toml +++ b/Project.toml @@ -9,22 +9,22 @@ CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" Graphics = "a2bd30eb-e257-5431-a919-1863eab51364" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" -ImageMetadata = "bc367c6b-8a6b-528e-b4bd-a4b897500b49" ImageAxes = "2803e5a7-5153-5ecf-9a86-9b4c37f5f5ac" +ImageMetadata = "bc367c6b-8a6b-528e-b4bd-a4b897500b49" Inflate = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LinearOperatorCollection = "a4a2c56f-fead-462a-a3ab-85921a5f2575" Mmap = "a63ad114-7e13-5084-954f-fe012c677804" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" -LinearOperatorCollection = "a4a2c56f-fead-462a-a3ab-85921a5f2575" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" @@ -35,28 +35,31 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" UnitfulAngles = "6fb2a4bd-7999-5318-a3b2-8ad61056cd98" [compat] +Aqua = "0.6" AxisArrays = "0.3, 0.4" CodecZlib = "0.7" -FileIO = "1.0" +DelimitedFiles = "1.9.0" +DocStringExtensions = "0.8, 0.9" FFTW = "1.3" +FileIO = "1.0" Graphics = "0.4, 1.0" HDF5 = "0.14, 0.15, 0.16" -ImageMetadata = "0.9" ImageAxes = "0.6" +ImageMetadata = "0.9" Inflate = "0.1.2" Interpolations = "0.12, 0.13, 0.14" LinearOperatorCollection = "1.1" StableRNGs = "1.0.0" -Unitful = "0.17, 1.2" +Unitful = "1.13, 1.14, 1.15, 1.16, 1.17" UnitfulAngles = "0.6.1" -DocStringExtensions = "0.8, 0.9" julia = "1.5" [extras] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3" LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3" Scratch = "6c6a2e73-6563-6170-7368-637461726353" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "HTTP", "LazyArtifacts", "Scratch"] +test = ["Aqua", "Test", "HTTP", "LazyArtifacts", "Scratch"] diff --git a/README.md b/README.md index 4e3756d9..0615574e 100644 --- a/README.md +++ b/README.md @@ -6,4 +6,5 @@ MPIFiles.jl is a Julia package for loading and storing magnetic particle imaging [![DOI](http://joss.theoj.org/papers/10.21105/joss.01331/status.svg)](https://doi.org/10.21105/joss.01331) [![Build status](https://github.com/MagneticParticleImaging/MPIFiles.jl/workflows/CI/badge.svg)](https://github.com/MagneticParticleImaging/MPIFiles.jl/actions) [![codecov.io](http://codecov.io/github/MagneticParticleImaging/MPIFiles.jl/coverage.svg?branch=master)](http://codecov.io/github/MagneticParticleImaging/MPIFiles.jl?branch=master) -[![License](https://img.shields.io/github/license/MagneticParticleImaging/MPIFiles.jl?color=green&style=flat)](https://github.com/MagneticParticleImaging/MPIFiles.jl/blob/master/LICENSE) \ No newline at end of file +[![License](https://img.shields.io/github/license/MagneticParticleImaging/MPIFiles.jl?color=green&style=flat)](https://github.com/MagneticParticleImaging/MPIFiles.jl/blob/master/LICENSE) +[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) \ No newline at end of file diff --git a/docs/src/lowlevel.md b/docs/src/lowlevel.md index 9af862f4..7c0d31d6 100644 --- a/docs/src/lowlevel.md +++ b/docs/src/lowlevel.md @@ -23,7 +23,7 @@ scannerFacility, scannerOperator, scannerManufacturer, scannerName, scannerTopol # acquisition parameters acqStartTime, acqNumFrames, acqNumAverages, acqGradient, acqOffsetField, -acqNumPeriodsPerFrame, acqSize +acqNumPeriodsPerFrame # drive-field parameters dfNumChannels, dfStrength, dfPhase, dfBaseFrequency, dfCustomWaveform, dfDivider, diff --git a/src/Conversion.jl b/src/Conversion.jl index bee5fee3..154ce335 100644 --- a/src/Conversion.jl +++ b/src/Conversion.jl @@ -1,5 +1,5 @@ # This file contains routines to generate MDF files -export saveasMDF, loadDataset, loadMetadata, loadMetadataOnline, setparam!, compressCalibMDF +export saveasMDF, loadDataset, loadMetadata, setparam!, compressCalibMDF export getFFdataPerPos, prepareAsMDFSingleMeasurement, convertCustomSF, blockAverage diff --git a/src/DatasetStore/BrukerDatasetStore.jl b/src/DatasetStore/BrukerDatasetStore.jl index c89b7e31..31a98974 100644 --- a/src/DatasetStore/BrukerDatasetStore.jl +++ b/src/DatasetStore/BrukerDatasetStore.jl @@ -1,10 +1,11 @@ -export BrukerDatasetStore, BrukerStore, findBrukerFiles, getNewCalibNum +export BrukerDatasetStore, findBrukerFiles, getNewCalibNum struct BrukerDatasetStore <: DatasetStore path::String end -if ispath("/opt/mpidata") +@static if ispath("/opt/mpidata") + export BrukerStore const BrukerStore = BrukerDatasetStore("/opt/mpidata") end diff --git a/src/DatasetStore/MDFDatasetStore.jl b/src/DatasetStore/MDFDatasetStore.jl index 9b9064bf..790af305 100644 --- a/src/DatasetStore/MDFDatasetStore.jl +++ b/src/DatasetStore/MDFDatasetStore.jl @@ -1,4 +1,4 @@ -export MDFDatasetStore, MDFStore, addStudy, getMDFStudyFolderName, calibdir, getMDFStore, getCalibStudy +export MDFDatasetStore, addStudy, getMDFStudyFolderName, calibdir, getMDFStore, getCalibStudy struct MDFDatasetStore <: DatasetStore path::String @@ -24,7 +24,8 @@ struct MDFDatasetStore <: DatasetStore end end -if ispath("/opt/mpidata/Bruker") +@static if ispath("/opt/mpidata/Bruker") + export MDFStore const MDFStore = MDFDatasetStore("/opt/data/Bruker") end path(e::Experiment{MDFDatasetStore}) = joinpath( path(e.study), string(e.num)*".mdf" ) diff --git a/src/MPIFiles.jl b/src/MPIFiles.jl index 5456cd59..b43cdecf 100644 --- a/src/MPIFiles.jl +++ b/src/MPIFiles.jl @@ -63,7 +63,7 @@ export scannerFacility, scannerOperator, scannerManufacturer, scannerName, # acquisition parameters export acqStartTime, acqNumFrames, acqNumAverages, - acqGradient, acqOffsetField, acqNumPeriodsPerFrame, acqSize + acqGradient, acqOffsetField, acqNumPeriodsPerFrame # drive-field parameters export dfNumChannels, dfStrength, dfPhase, dfBaseFrequency, dfCustomWaveform, diff --git a/src/Measurements.jl b/src/Measurements.jl index cd63dae4..f901d3c2 100644 --- a/src/Measurements.jl +++ b/src/Measurements.jl @@ -1,4 +1,4 @@ -export getMeasurements, getMeasurementsFD, getMeasurementsLowLevel +export getMeasurements, getMeasurementsFD function measDataConv(f::MPIFile, args...) data = measDataTD(f, args...) diff --git a/src/Positions/Positions.jl b/src/Positions/Positions.jl index f8a53664..a7d7765a 100644 --- a/src/Positions/Positions.jl +++ b/src/Positions/Positions.jl @@ -4,7 +4,7 @@ export Positions, GridPositions, RegularGridPositions, ChebyshevGridPositions, export SpatialDomain, AxisAlignedBox, Ball export loadTDesign, getPermutation export fieldOfView, fieldOfViewCenter, shape -export idxToPos, posToIdx, posToLinIdx, spacing, isSubgrid, deriveSubgrid, toDict +export posToIdx, posToLinIdx, spacing, isSubgrid, deriveSubgrid, toDict abstract type Positions end abstract type GridPositions<:Positions end diff --git a/test/runtests.jl b/test/runtests.jl index 994af6cb..6080e3d6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,6 +9,7 @@ using UUIDs using Unitful using Scratch using LazyArtifacts +using Aqua const datadir = joinpath(artifact"data", "data") @info "The test data is located at $datadir." @@ -21,6 +22,10 @@ mkpath(joinpath(tmpdir,"mdfim")) mkpath(joinpath(tmpdir,"positions")) mkpath(joinpath(tmpdir,"transferFunction")) +@testset "Aqua" begin + Aqua.test_all(MPIFiles) +end + include("General.jl") include("Cartesian.jl") include("Positions.jl") From 5a2f6e311615e404cc3e65d5077f91a55dcf77fe Mon Sep 17 00:00:00 2001 From: Jonas Schumacher Date: Wed, 27 Sep 2023 09:45:38 +0200 Subject: [PATCH 07/23] Fix CI for docs (#65) * Update Positions.jl * Try fix for Documenter update --- docs/Project.toml | 3 +++ docs/make.jl | 3 ++- src/Positions/Positions.jl | 2 +- 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index dfa65cd1..d31eb643 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,2 +1,5 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" + +[compat] +Documenter = "1" \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index db13bce3..7eaab185 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -15,7 +15,8 @@ makedocs( "Frequency Filter" => "frequencyFilter.md", "Reconstruction Results" => "reconstruction.md" # "Positions" => "positions.md" - ] + ], + warnonly = [:missing_docs] ) deploydocs( diff --git a/src/Positions/Positions.jl b/src/Positions/Positions.jl index a7d7765a..0c2560de 100644 --- a/src/Positions/Positions.jl +++ b/src/Positions/Positions.jl @@ -220,7 +220,7 @@ end function posToIdxFloat(grid::RegularGridPositions, pos) idx = 0.5 .* (shape(grid) .* ((pos .- fieldOfViewCenter(grid)) ./ ( 0.5 .* fieldOfView(grid) ) .+ 1) .+ 1) - idx = [isnan(val) ? one(eltype(idx)) : val for val in idx] + #idx = [isnan(val) ? one(eltype(idx)) : val for val in idx] return idx end From 107ffaf80174584de50be67d501c64cef3e84a46 Mon Sep 17 00:00:00 2001 From: Jonas Schumacher Date: Wed, 27 Sep 2023 14:07:27 +0200 Subject: [PATCH 08/23] Change `getindex` to return KeyError (#62) * Change `getindex` to return KeyErrror * Fix errors with missing induction factor * Some checks * Fix stupid logic mistake * Remove warning * Fix tests --- src/Conversion.jl | 20 ++++++++------------ src/MDF.jl | 12 +++--------- src/MDFInMemory.jl | 7 +++++-- src/Measurements.jl | 9 +++++++++ 4 files changed, 25 insertions(+), 23 deletions(-) diff --git a/src/Conversion.jl b/src/Conversion.jl index 154ce335..90eebc23 100644 --- a/src/Conversion.jl +++ b/src/Conversion.jl @@ -4,18 +4,18 @@ export getFFdataPerPos, prepareAsMDFSingleMeasurement, convertCustomSF, blockAve function setparam!(params::Dict, parameter, value) - if value != nothing + if !(isnothing(value) || ismissing(value)) params[parameter] = value end end # we do not support all conversion possibilities function loadDataset(f::MPIFile; frames=1:acqNumFrames(f), applyCalibPostprocessing=false, - numPeriodAverages=1, numPeriodGrouping=1, experimentNumber=nothing, + numPeriodAverages=1, numPeriodGrouping=1, experimentNumber=missing, fixDistortions=false, kargs...) params = loadMetadata(f) - if experimentNumber != nothing + if !ismissing(experimentNumber) params[:experimentNumber] = experimentNumber end @@ -63,7 +63,7 @@ function loadDataset(f::MPIFile; frames=1:acqNumFrames(f), applyCalibPostprocess cat(zeros(Bool,acqNumFGFrames(f)),ones(Bool,acqNumBGFrames(f)), dims=1)) snr = calibSNR(f) - if snr == nothing + if isnothing(snr) @info "calculate SNR" snr = calculateSystemMatrixSNR(f, data, numPeriodAverages=numPeriodAverages, numPeriodGrouping=numPeriodGrouping) end @@ -151,8 +151,6 @@ function loadMeasParams(f, params = Dict{Symbol,Any}(); skipMeasData = false) return params end - - function appendBGDataset(params::Dict, filenameBG::String; kargs...) params = MPIFile(filenameBG) do fBG appendBGDataset(params, fBG; kargs...) @@ -171,7 +169,6 @@ function appendBGDataset(params::Dict, fBG::MPIFile; frames=1:acqNumFrames(fBG)) return params end - isConvertibleToMDF(f::MPIFile) = true function isConvertibleToMDF(f::BrukerFile) # check if raw data is consistent @@ -196,16 +193,15 @@ function saveasMDF(filenameOut::String, filenameIn::String; kargs...) end function saveasMDF(filenameOut::String, f::MPIFile; filenameBG = nothing, enforceConversion=false, kargs...) - # This is a hack. Needs to be fixed properly - if(haskey(kargs, :SNRThresh) || haskey(kargs, :sparsityTrafoRedFactor)) && calibSNR(f) != nothing + if(haskey(kargs, :SNRThresh) || haskey(kargs, :sparsityTrafoRedFactor)) && !isnothing(calibSNR(f)) compressCalibMDF(filenameOut, f; kargs...) return end if enforceConversion || isConvertibleToMDF(f) - params = loadDataset(f;kargs...) - if filenameBG != nothing + params = loadDataset(f; kargs...) + if !isnothing(filenameBG) appendBGDataset(params, filenameBG) end saveasMDF(filenameOut, params) @@ -746,7 +742,7 @@ function saveasMDF(file::HDF5.File, params::Dict{Symbol,Any}) write(file, "/scanner/topology", get(params,:scannerTopology,"FFP")) # acquisition parameters - write(file, "/acquisition/numAverages", params[:acqNumAverages]) + write(file, "/acquisition/numAverages", params[:acqNumAverages]) write(file, "/acquisition/numFrames", get(params,:acqNumFrames,1)) write(file, "/acquisition/numPeriodsPerFrame", get(params,:acqNumPeriodsPerFrame,1)) write(file, "/acquisition/startTime", "$( get(params,:acqStartTime, Dates.unix2datetime(time())) )") diff --git a/src/MDF.jl b/src/MDF.jl index f4a22856..899042a6 100644 --- a/src/MDF.jl +++ b/src/MDF.jl @@ -69,13 +69,7 @@ function h5haskey(filename, parameter) end end -function getindex(f::MDFFile, parameter) - if haskey(f.file, parameter) - return read(f.file, parameter) - else - return nothing - end -end +getindex(f::MDFFile, parameter) = read(f.file, parameter) function getindex(f::MDFFile, parameter, default) #if !haskey(f.param_cache,parameter) @@ -106,8 +100,8 @@ studyUuid(f::MDFFileV2) = @keyrequired UUID(f["/study/uuid"]) studyDescription(f::MDFFileV1)::Union{String, Missing} = "n.a." studyDescription(f::MDFFileV2)::Union{String, Missing} = @keyrequired f["/study/description"] function studyTime(f::MDFFile) - t = f["/study/time"] - if typeof(t)==String + t = @keyoptional f["/study/time"] + if typeof(t) == String return DateTime(t) else return nothing diff --git a/src/MDFInMemory.jl b/src/MDFInMemory.jl index 8eb7adb0..7e87f810 100644 --- a/src/MDFInMemory.jl +++ b/src/MDFInMemory.jl @@ -1407,7 +1407,11 @@ function inMemoryMDFFromMDFFileV2(mdfFile::MDFFileV2)::MDFv2InMemory end # Add measurements data - measDataRaw(mdf, mdfFile.mmap_measData) + if !isnothing(mdfFile.mmap_measData) + measDataRaw(mdf, mdfFile.mmap_measData) + else + @warn "The measurement data could not be read. Please check closely." + end return mdf end @@ -1477,7 +1481,6 @@ measObservedDriveField(mdf::MDFv2InMemory, measDriveFields) = mdf.custom["measOb measAppliedDriveField(mdf::MDFv2InMemory) = @keyoptional mdf.custom["measAppliedDriveField"] measAppliedDriveField(mdf::MDFv2InMemory, measTransmit) = mdf.custom["measAppliedDriveField"] = measTransmit - auxiliaryData(mdf::MDFv2InMemory) = @keyoptional mdf.custom["auxiliaryData"] auxiliaryData(mdf::MDFv2InMemory, auxiliaryData) = mdf.custom["auxiliaryData"] = auxiliaryData diff --git a/src/Measurements.jl b/src/Measurements.jl index f901d3c2..d57eda45 100644 --- a/src/Measurements.jl +++ b/src/Measurements.jl @@ -268,6 +268,10 @@ function getMeasurements(f::MPIFile, neglectBGFrames=true; if tfCorrection && !measIsTFCorrected(f) tf = rxTransferFunction(f) inductionFactor = rxInductionFactor(f) + if ismissing(inductionFactor) + @warn "The file is missing the induction factor. The induction factor will be set to 1." + inductionFactor = ones(Float64, rxNumChannels(f)) + end J = size(data,1) dataF = rfft(data, 1) @@ -323,6 +327,11 @@ function getMeasurementsFD(f::MPIFile, args...; if tfCorrection && !measIsTFCorrected(f) tf = rxTransferFunction(f) inductionFactor = rxInductionFactor(f) + if ismissing(inductionFactor) + @warn "The file is missing the induction factor. The induction factor will be set to 1." + inductionFactor = ones(Float64, rxNumChannels(f)) + end + data[2:end,:,:,:] ./= tf[2:end,:,:,:] if all(tf[1,:,:,:] .!= 0) && !any(isnan.(tf[1,:,:,:])) From e3c9557a59baacfc8412d40a7046f91c627e9b5d Mon Sep 17 00:00:00 2001 From: nHackel Date: Mon, 9 Oct 2023 15:27:47 +0000 Subject: [PATCH 09/23] Fix bug in system matrix loading for brukerfile where original row vector was mistakenly mutated --- src/Brukerfile.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/Brukerfile.jl b/src/Brukerfile.jl index d01eda03..9b11ff19 100644 --- a/src/Brukerfile.jl +++ b/src/Brukerfile.jl @@ -540,17 +540,20 @@ function systemMatrix(b::BrukerFileCalib, rows, bgCorrection=true) nFreq = div(rxNumSamplingPoints(b)*numSubPeriods(b),2)+1 if numSubPeriods(b) > 1 + rows_ = collect(rows) stepsize = numSubPeriods(b) for k=1:length(rows) freq = rows[k][1] rec = rows[k][2] - rows[k] = CartesianIndex((freq-1)*stepsize+1, rec) + rows_[k] = CartesianIndex((freq-1)*stepsize+1, rec) end + else + rows_ = rows end s = open(sfFilename) data = Mmap.mmap(s, Array{ComplexF64,3}, (prod(calibSize(b)), nFreq, rxNumChannels(b))) - S = data[:,rows] + S = data[:,rows_] close(s) rmul!(S, 1.0/acqNumAverages(b)) return S From 5b023427be6359a78cab599f75ab8ba270a1bcd4 Mon Sep 17 00:00:00 2001 From: Jonas Schumacher Date: Wed, 18 Oct 2023 12:31:34 +0200 Subject: [PATCH 10/23] Fix error with missing conversion factor --- src/Measurements.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Measurements.jl b/src/Measurements.jl index d57eda45..1afde70c 100644 --- a/src/Measurements.jl +++ b/src/Measurements.jl @@ -7,7 +7,7 @@ function measDataConv(f::MPIFile, args...) data = map(Float32, data) end a = rxDataConversionFactor(f) - if a!=nothing + if !ismissing(a) for d=1:size(data,2) slice = view(data,:,d,:,:) rmul!(slice, a[1,d]) From 5d706bfa9a26fe48ffbf8f1f1b456a084c0fe74c Mon Sep 17 00:00:00 2001 From: "Justin Ackers @ mpi05" Date: Wed, 22 Nov 2023 14:48:46 +0100 Subject: [PATCH 11/23] add compat for aqua --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 70ad17a2..cfcd4691 100644 --- a/Project.toml +++ b/Project.toml @@ -53,6 +53,7 @@ LinearOperatorCollection = "1.1" StableRNGs = "1.0.0" Unitful = "1.13, 1.14, 1.15, 1.16, 1.17" UnitfulAngles = "0.6.1" +UnitfulParsableString = "0.1.6" julia = "1.5" [extras] From 61a711220d257141f9c4259012f90fc448656e7f Mon Sep 17 00:00:00 2001 From: Jonas Schumacher Date: Thu, 23 Nov 2023 15:55:11 +0100 Subject: [PATCH 12/23] `dataConversionFactor` and `inductionFactor` are optional --- src/Measurements.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Measurements.jl b/src/Measurements.jl index 1afde70c..0a512acd 100644 --- a/src/Measurements.jl +++ b/src/Measurements.jl @@ -7,7 +7,7 @@ function measDataConv(f::MPIFile, args...) data = map(Float32, data) end a = rxDataConversionFactor(f) - if !ismissing(a) + if !isnothing(a) for d=1:size(data,2) slice = view(data,:,d,:,:) rmul!(slice, a[1,d]) @@ -268,7 +268,7 @@ function getMeasurements(f::MPIFile, neglectBGFrames=true; if tfCorrection && !measIsTFCorrected(f) tf = rxTransferFunction(f) inductionFactor = rxInductionFactor(f) - if ismissing(inductionFactor) + if isnothing(inductionFactor) @warn "The file is missing the induction factor. The induction factor will be set to 1." inductionFactor = ones(Float64, rxNumChannels(f)) end @@ -327,7 +327,7 @@ function getMeasurementsFD(f::MPIFile, args...; if tfCorrection && !measIsTFCorrected(f) tf = rxTransferFunction(f) inductionFactor = rxInductionFactor(f) - if ismissing(inductionFactor) + if isnothing(inductionFactor) @warn "The file is missing the induction factor. The induction factor will be set to 1." inductionFactor = ones(Float64, rxNumChannels(f)) end From 4efd8566fc62820307aa5e455c83a18316fa87dd Mon Sep 17 00:00:00 2001 From: Jonas Schumacher Date: Thu, 23 Nov 2023 16:03:56 +0100 Subject: [PATCH 13/23] Change macro to @keyoptional --- src/MDF.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MDF.jl b/src/MDF.jl index 899042a6..91fe4023 100644 --- a/src/MDF.jl +++ b/src/MDF.jl @@ -263,12 +263,12 @@ function rxHasTransferFunction(f::MDFFile) haskey(f.file, "/acquisition/receiver/transferFunction") end rxInductionFactor(f::MDFFileV1) = nothing -rxInductionFactor(f::MDFFileV2) = @keyrequired f["/acquisition/receiver/inductionFactor"] +rxInductionFactor(f::MDFFileV2) = @keyoptional f["/acquisition/receiver/inductionFactor"] rxUnit(f::MDFFileV1)::Union{String, Missing} = "a.u." rxUnit(f::MDFFileV2)::Union{String, Missing} = @keyrequired f["/acquisition/receiver/unit"] rxDataConversionFactor(f::MDFFileV1) = repeat([1.0, 0.0], outer=(1,rxNumChannels(f))) -rxDataConversionFactor(f::MDFFileV2) = @keyrequired f["/acquisition/receiver/dataConversionFactor"] +rxDataConversionFactor(f::MDFFileV2) = @keyoptional f["/acquisition/receiver/dataConversionFactor"] # measurements function measData(f::MDFFileV1, frames=1:acqNumFrames(f), periods=1:acqNumPeriodsPerFrame(f), From 1d924456b53bb2910a71947f6e24f070bba9caa4 Mon Sep 17 00:00:00 2001 From: Jonas Schumacher Date: Tue, 5 Dec 2023 17:12:15 +0100 Subject: [PATCH 14/23] Fix getMeasurements for files without bg frames when setting bgCorrection = true --- src/Measurements.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/Measurements.jl b/src/Measurements.jl index 0a512acd..dd134f0d 100644 --- a/src/Measurements.jl +++ b/src/Measurements.jl @@ -196,10 +196,14 @@ function getMeasurements(f::MPIFile, neglectBGFrames=true; data = getAveragedMeasurements(f; frames=idx[frames], numAverages=numAverages, kargs...) - - if bgCorrection + + idxBG = measBGFrameIdx(f) + hasBGFrames = length(idxBG) > 0 + if bgCorrection && !hasBGFrames + @warn "Background correction was selected but there are no background frames in the file." + elseif bgCorrection && hasBGFrames @debug "Applying bg correction ..." - idxBG = measBGFrameIdx(f) + dataBG = getAveragedMeasurements(f; frames=idxBG, kargs...) if interpolateBG blockLen = measBGFrameBlockLengths(measIsBGFrame(f)) @@ -245,8 +249,7 @@ function getMeasurements(f::MPIFile, neglectBGFrames=true; data = getAveragedMeasurements(f; frames=frames, numAverages=numAverages, kargs...) end - if bgCorrection - idxBG = measBGFrameIdx(f) + if bgCorrection && hasBGFrames dataBG = getAveragedMeasurements(f; frames=idxBG, kargs...) data[:,:,:,:] .-= mean(dataBG, dims=4) From e110a272e8e57442ce0575e0e63b2c5b5258f836 Mon Sep 17 00:00:00 2001 From: Jonas Schumacher Date: Wed, 6 Dec 2023 11:03:02 +0100 Subject: [PATCH 15/23] Bump version for patches --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index cd7854a0..2280aa45 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MPIFiles" uuid = "371237a9-e6c1-5201-9adb-3d8cfa78fa9f" authors = ["Tobias Knopp "] -version = "0.14.0" +version = "0.14.1" [deps] AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9" From 2df65176f87c3a6bf4d130aa8d6483d6988e801d Mon Sep 17 00:00:00 2001 From: Artyom Tsanda Date: Fri, 8 Dec 2023 13:52:37 +0100 Subject: [PATCH 16/23] fix for the data type used in MDFv2Reconstruction --- src/MDFInMemory.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MDFInMemory.jl b/src/MDFInMemory.jl index 7e87f810..32af6a26 100644 --- a/src/MDFInMemory.jl +++ b/src/MDFInMemory.jl @@ -520,7 +520,7 @@ Reconstruction group of an in-memory MDF """ mutable struct MDFv2Reconstruction <: MDFv2InMemoryPart "Reconstructed data" - data::Union{Array{Number, 3}, Missing} + data::Union{Array{Float32, 3}, Missing} "Field of view of reconstructed data; optional" fieldOfView::Union{Vector{Float64}, Nothing} "Center of the reconstructed data (relative to scanner origin/center); optional" From e1e9190e5fd9ece35cbd320225832ead01a70374 Mon Sep 17 00:00:00 2001 From: nHackel Date: Fri, 15 Dec 2023 10:36:25 +0100 Subject: [PATCH 17/23] Make tracer injection time optional --- src/MDF.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/MDF.jl b/src/MDF.jl index 91fe4023..93aafd8f 100644 --- a/src/MDF.jl +++ b/src/MDF.jl @@ -143,14 +143,15 @@ tracerSolute(f::MDFFileV2)::Union{Vector{String}, Missing} = @keyrequired _makeS tracerSolute(f::MDFFileV1)::Union{Vector{String}, Missing} = ["Fe"] function tracerInjectionTime(f::MDFFile)::Union{Vector{DateTime}, Nothing} p = typeof(f) <: MDFFileV1 ? "/tracer/time" : "/tracer/injectionTime" - if isnothing(f[p]) + time = @keyoptional f[p] + if isnothing(nothing) return nothing end - if typeof(f[p]) == String - return [DateTime(f[p])] + if typeof(time) == String + return [DateTime(time)] else - return [DateTime(y) for y in f[p]] + return [DateTime(y) for y in time] end end #tracerInjectionTime(f::MDFFileV2) = DateTime( f["/tracer/injectionTime"] ) From f500e355de645a26894e22ba9919d6a1ce11d7f4 Mon Sep 17 00:00:00 2001 From: Tobias Knopp Date: Tue, 2 Jan 2024 14:23:11 +0100 Subject: [PATCH 18/23] fix major bug --- src/FrequencyFilter.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/FrequencyFilter.jl b/src/FrequencyFilter.jl index ba6a0a05..cf947b08 100644 --- a/src/FrequencyFilter.jl +++ b/src/FrequencyFilter.jl @@ -45,7 +45,7 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0, filterFrequenciesByMinFreq!(freqIndices, f, minFreq) end - if maxFreq < nFreq + if maxFreq < rxBandwidth(f) filterFrequenciesByMaxFreq!(freqIndices, f, maxFreq) end From e771bc5637ed838680c8dbf7aaff1c3aa3d91c0a Mon Sep 17 00:00:00 2001 From: "Justin Ackers @ mpi05" Date: Mon, 8 Jan 2024 16:07:14 +0100 Subject: [PATCH 19/23] improved and exported "combine" function --- src/TransferFunction.jl | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/TransferFunction.jl b/src/TransferFunction.jl index b5afd6d5..56ccbcdf 100644 --- a/src/TransferFunction.jl +++ b/src/TransferFunction.jl @@ -1,4 +1,4 @@ -export TransferFunction, sampleTF, setTF +export TransferFunction, sampleTF, setTF, combine mutable struct TransferFunction freq::Vector{Float64} @@ -80,13 +80,23 @@ function load_tf(filename::String) return TransferFunction(freq, tf, inductionFactor=inductionFactor, units=uparse.(unit)) end -function combine(tf1::TransferFunction, tf2::TransferFunction) - if tf1.freq != tf2.freq; error("The frequency axes of the transfer functions do not match. Can not combine!") end - freq = tf1.freq - data = cat(tf1.data,tf2.data, dims=2) - inductionFactor = cat(tf1.inductionFactor, tf2.inductionFactor, dims=1) - units = cat(tf1.units, tf2.units, dims=1) - return TransferFunction(freq, data, inductionFactor=inductionFactor, units=units) +""" +$(TYPEDSIGNATURES) + +Combine two `TransferFunctions` along their channel dimension. If interpolate=false, will only work if the frequency samples are identical. +""" +function combine(tf1::TransferFunction, tf2::TransferFunction; interpolate=false) + if !interpolate + if tf1.freq != tf2.freq; error("The frequency axes of the transfer functions do not match. Can not combine!") end + freq = tf1.freq + data = cat(tf1.data,tf2.data, dims=2) + else + freq = unique(sort(cat(tf1.freq, tf2.freq, dims=1))) + data = cat(tf1(freq, :), tf2(freq, :), dims=2) + end + inductionFactor = cat(tf1.inductionFactor, tf2.inductionFactor, dims=1) + units = cat(tf1.units, tf2.units, dims=1) + return TransferFunction(freq, data, inductionFactor=inductionFactor, units=units) end function save(filename::String, tf::TransferFunction) From a5d02abc59bcd8d03ff748d44f898afff4898221 Mon Sep 17 00:00:00 2001 From: "Justin Ackers @ mpi05" Date: Mon, 8 Jan 2024 16:07:46 +0100 Subject: [PATCH 20/23] added documentation and small fixes --- Project.toml | 2 +- docs/make.jl | 3 +- docs/src/transferfunction.md | 69 +++++++++++++++++++++++++++++++ src/TransferFunction.jl | 78 +++++++++++++++++++++++++++++++----- 4 files changed, 139 insertions(+), 13 deletions(-) create mode 100644 docs/src/transferfunction.md diff --git a/Project.toml b/Project.toml index cfcd4691..3d923d0a 100644 --- a/Project.toml +++ b/Project.toml @@ -39,7 +39,7 @@ UnitfulParsableString = "06c00241-927a-4d5b-bb5e-6b5a2ada3567" Aqua = "0.6" AxisArrays = "0.3, 0.4" CodecZlib = "0.7" -DelimitedFiles = "1.9.0" +DelimitedFiles = "1" DocStringExtensions = "0.8, 0.9" FFTW = "1.3" FileIO = "1.0" diff --git a/docs/make.jl b/docs/make.jl index 7eaab185..b38018b9 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -13,7 +13,8 @@ makedocs( "Measurements" => "measurements.md", "System Matrices" => "systemmatrix.md", "Frequency Filter" => "frequencyFilter.md", - "Reconstruction Results" => "reconstruction.md" + "Reconstruction Results" => "reconstruction.md", + "Transfer Functions" => "transferfunction.md" # "Positions" => "positions.md" ], warnonly = [:missing_docs] diff --git a/docs/src/transferfunction.md b/docs/src/transferfunction.md new file mode 100644 index 00000000..0c9dc5b2 --- /dev/null +++ b/docs/src/transferfunction.md @@ -0,0 +1,69 @@ +# Transfer functions + +MPIFiles defines the type `TransferFunction` which represents the system-properties of a linear time-invariant system in frequency space ([see also](https://en.wikipedia.org/wiki/Transfer_function)). In the context of MPIFiles this usually includes the properties of the receive chain of an MPI system, but can also hold information about other parts of the signal chain. + +## Basic construction and data access +The `TransferFunction` object is constructed from samples in frequency space and offers two ways to access the underlying data. The first uses indexing to directly access the underlying data: + +```@setup tf +using MPIFiles +``` +```@repl tf + f = collect(range(0,1e6,step=1e3)); + # TransferFunction of a simple lowpass filter + tf = TransferFunction(f, 1 ./ (1 .+ im*f/1e4 )) + tf[1] # Value of tf at first frequency sample (f = 0 Hz) +``` + +The second method has to be called on the object itself and uses linear interpolation to access the tf at any frequency: +```@repl tf + tf(0) # Value of tf at f = 0 Hz + tf(1e4) # Value of tf at f = 10000 Hz +``` + +A `TransferFunction` can have multiple channels, which adds a second index to both data access functions. Directly accessing multiple channels is also possible. The complex array used to construct the `TransferFunction` needs to have the shape [number of frequency samples, channels]. + +```@repl tf +tf = TransferFunction(f, [1 ./(1 .+im*f/1e4) 1 ./(1 .+im*f/1e3)]) +tf[11,1] +tf[11,2] +tf(1e4,1) +tf(1e4,2) +tf(1e4, [1,2]) +tf(1e4,:) +``` +## Units +To attach units to the `TransferFunction` the keyword-parameter `units` can to be used to give a `Unitful` unit to every channel of the tf. This can be useful if the transfer function is not dimensionless but relates two physical quantities, e.g. voltage and current in the form of an impedance. All **interpolated** accesses to tf data then return a `Unitful.Quantity`. + +```@repl tf +R = 1; # Ohm +L = 10e-6; # Henry +RL = TransferFunction(f, R .+ im*2pi*f*L, units=["V/A"]) +RL([0,100e3]) +``` + +## Saving and loading +A `TransferFunction` object can be saved to and loaded from a .h5 file. + +```@docs +MPIFiles.save(filename::String, tf::TransferFunction) +MPIFiles.TransferFunction(::String) +``` + +## Additional constructors + +In addition to the constructor taking a single (complex) array it is also possible to give two arrays representing amplitude and phase. + +It is also possible to construct a `TransferFunction` from the transfer function data included in an `MPIFile`. + +```@docs + MPIFiles.TransferFunction(freq::Vector{<:Real}, ampdata::Array{<:Real,N}, phasedata::Array{<:Real,N}) where N + MPIFiles.TransferFunction(::MPIFile) +``` + +## Other interesting functions +```@docs + TransferFunction + combine(::TransferFunction, ::TransferFunction) +``` + diff --git a/src/TransferFunction.jl b/src/TransferFunction.jl index 56ccbcdf..7349ca90 100644 --- a/src/TransferFunction.jl +++ b/src/TransferFunction.jl @@ -1,5 +1,19 @@ export TransferFunction, sampleTF, setTF, combine +""" +$(TYPEDEF) +$(TYPEDFIELDS) + + + TransferFunction(freq_::Vector{<:Real}, datain::Array{<:Complex}; inductionFactor::Vector{<:Real}=ones(size(datain, 2)), units::Vector=Unitful.FreeUnits[Unitful.NoUnits for i in 1:size(datain, 2)]) + +Create a `TransferFunction` from a complex data array at frequencies `freq_`. + +# Optional Keyword-Arguments: +- `inductionFactor::Vector{<:Real}`: induction factor for each channel +- `units::Vector`: units for each channel, can be either Unitful.FreeUnits or a string that can be parsed as a Unitful unit + +""" mutable struct TransferFunction freq::Vector{Float64} data::Matrix{ComplexF64} @@ -7,8 +21,8 @@ mutable struct TransferFunction inductionFactor::Vector{Float64} units::Vector{Unitful.FreeUnits} - function TransferFunction(freq_::Vector{<:Real}, datain::Array{<:Complex}; inductionFactor::Vector{<:Real}=ones(size(datain, 2)), units::Vector=Unitful.FreeUnits[Unitful.NoUnits for i in 1:size(datain, 2)]) + function TransferFunction(freq_::Vector{<:Real}, datain::Array{<:Number}; inductionFactor::Vector{<:Real}=ones(size(datain, 2)), units::Vector=Unitful.FreeUnits[Unitful.NoUnits for i in 1:size(datain, 2)]) parsed_units = Unitful.FreeUnits[] for tmp in units if isa(tmp, String); tmp = uparse(tmp) end # get correct unit from unit strings @@ -28,12 +42,26 @@ end Base.show(io::IO, ::MIME"text/plain", tf::TransferFunction) = print(io, "MPIFiles.TransferFunction: \n\t$(size(tf.data,2)) channel(s), units of $(string.(tf.units))\n\t$(size(tf.data,1)) frequency samples from $(tf.freq[1]) Hz to $(tf.freq[end]) Hz") -function TransferFunction(freq_::Vector{<:Real}, ampdata::Array{<:Real,N}, phasedata::Array{<:Real,N}; kwargs...) where N +""" +$(TYPEDSIGNATURES) + +Create a `TransferFunction` from separate amplitude and phase arrays at frequencies `freq`. + +`ampdata` and `phasedata` should have the following shape: [frequencies, channels] +""" +function TransferFunction(freq::Vector{<:Real}, ampdata::Array{<:Real,N}, phasedata::Array{<:Real,N}; kwargs...) where N if size(ampdata) != size(phasedata); error("The size of ampdata and phasedata must match!") end data = ampdata.*exp.(im.*phasedata) - return TransferFunction(freq_, data; kwargs...) + return TransferFunction(freq, data; kwargs...) end +""" +$(TYPEDSIGNATURES) + +Create a `TransferFunction` from a data file at `filename`. + +The file can be either a h5-File created with this package. Keyword arguments will be passed to `load_tf_fromVNA` +""" function TransferFunction(filename::String; kargs...) filenamebase, ext = splitext(filename) if ext == ".h5" @@ -44,6 +72,11 @@ function TransferFunction(filename::String; kargs...) return tf end +""" +$(TYPEDSIGNATURES) + +Create a `TransferFunction` from the tf data saved in a MPIFile (see `rxTransferFunction`) +""" function TransferFunction(file::MPIFile) tf_file = rxTransferFunction(file) inductionFactor = rxInductionFactor(file) @@ -51,21 +84,38 @@ function TransferFunction(file::MPIFile) return TransferFunction(f, abs.(tf_file), angle.(tf_file), inductionFactor) end -function getindex(tmf::TransferFunction, args...) +""" + tf[i,j] + +Directly access the underlying data of a `TransferFunction` +""" +function getindex(tf::TransferFunction, args...) try - return getindex(tmf.data, args...) + return getindex(tf.data, args...) catch e @warn "The indexing using square brackets on TransferFunction objects now always operates on the integer indizes of the underlying transfer function data. To use frequency interpolation, use tf(freq, channel) instead of tf[[freq],channel]." rethrow(e) end end -function (tmf::TransferFunction)(x, chan::Integer=1) - if chan>length(tmf.interpolator); error("The TransferFunction only has $(length(tmf.interpolator)) channel(s), unable to access channel $(chan)") end - return tmf.interpolator[chan](x) .* tmf.units[chan] +""" + tf(f, chan::Integer) + +Interpolated access to a `TransferFunction` at frequencies `f` and single channel `chan` +""" +function (tf::TransferFunction)(f, chan::Integer=1) + if chan>length(tf.interpolator); error("The TransferFunction only has $(length(tf.interpolator)) channel(s), unable to access channel $(chan)") end + return tf.interpolator[chan](f) .* tf.units[chan] end -(tmf::TransferFunction)(x, chan::AbstractArray) = hcat([tmf(x,c) for c in chan]...) +""" + tf(f, chan::AbstractVector{<:Integer}) + +Interpolated access to a `TransferFunction` at frequencies `f` and channels `chan` +""" +(tf::TransferFunction)(f, chan::AbstractVector{<:Integer}) = hcat([tf(f,c) for c in chan]...) + +(tf::TransferFunction)(f, ::Colon) = tf(f, axes(tf.data,2)) function load_tf(filename::String) tf = h5read(filename,"/transferFunction") @@ -99,6 +149,11 @@ function combine(tf1::TransferFunction, tf2::TransferFunction; interpolate=false return TransferFunction(freq, data, inductionFactor=inductionFactor, units=units) end +""" +$(TYPEDSIGNATURES) + +Save `tf` as a h5 file to `filename` +""" function save(filename::String, tf::TransferFunction) h5write(filename, "/transferFunction", tf.data) h5write(filename, "/frequencies", tf.freq) @@ -107,6 +162,7 @@ function save(filename::String, tf::TransferFunction) return nothing end +"""TODO: fix function and write docs""" function load_tf_fromVNA(filename::String; frequencyWeighting=false, R = 50.0, #Ω @@ -158,10 +214,10 @@ function load_tf_fromVNA(filename::String; end -function sampleTF(tmf::TransferFunction, f::MPIFile) +function sampleTF(tf::TransferFunction, f::MPIFile) freq = rxFrequencies(f) numChan = rxNumChannels(f) - return tmf(freq,1:numChan) + return tf(freq,1:numChan) end From a4bb0a17bd9cb1d79184aa1423af3d3918ade3ed Mon Sep 17 00:00:00 2001 From: "Justin Ackers @ mpi05" Date: Mon, 8 Jan 2024 16:07:51 +0100 Subject: [PATCH 21/23] added tests --- test/TransferFunction.jl | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/test/TransferFunction.jl b/test/TransferFunction.jl index 1d4c2d98..7d2b8a16 100644 --- a/test/TransferFunction.jl +++ b/test/TransferFunction.jl @@ -26,4 +26,40 @@ tfh5path = joinpath(tmpdir,"transferFunction","example.h5") measBruker = MPIFile(fnMeasBruker) tf = sampleTF(c, measBruker) @test size(tf) == (rxNumFrequencies(measBruker), rxNumChannels(measBruker)) + + f = collect(range(0,1e6,step=1e3)); + data = 1 ./ (1 .+ im*f/1e4 ); + tf = TransferFunction(f, data) # TransferFunction of a simple lowpass filter + @test tf[1] == 1.0 + @test tf[1:6,1] == data[1:6] + @test tf(0) == tf(0,1) + + data = [1 ./(1 .+im*f/1e4) 1 ./(1 .+im*f/1e3)] + tf = TransferFunction(f, data) + @test tf[11,1] == data[11,1] + @test tf[11,2] == data[11,2] + @test tf[11,:] == data[11,:] + @test tf(1e4, [1,2]) == [data[11,1] 1 ./(1 .+im*1e4/1e3)] + @test tf(1e4,:) == tf(1e4, [1,2]) + @test tf(1e4,:) == tf(1e4, 1:2) + @test tf(0) == tf(0,1) + + tf = TransferFunction(f, data, units=["V/V", "A/V"]) + @test tf(0) == 1u"V/V" + @test tf(0,2) == 1u"A/V" + + + f1 = collect(range(0,1e6,step=1e3)); + data1 = 1 ./ (1 .+ im*f1/1e4 ); + + f2 = collect(range(0,1e7,step=1e4)); + data2 = 1 ./ (1 .+ im*f2/1e4 ); + + tf1 = TransferFunction(f1, data1) + tf2 = TransferFunction(f2, data2) + + @test_throws ErrorException tf_c = combine(tf1,tf2) + tf_c = combine(tf1,tf2,interpolate=true) + @test tf_c(101.5e3, :) ≈ [tf1(101.5e3) tf2(101.5e3)] atol=1e-12 + end From b57e6d65bb3d7882f36ed4591d9cda6458e98d8d Mon Sep 17 00:00:00 2001 From: "Justin Ackers @ mpi05" Date: Mon, 8 Jan 2024 16:14:58 +0100 Subject: [PATCH 22/23] bugfix for tracerInjectionTime --- src/MDF.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MDF.jl b/src/MDF.jl index 93aafd8f..ea84a2a8 100644 --- a/src/MDF.jl +++ b/src/MDF.jl @@ -144,7 +144,7 @@ tracerSolute(f::MDFFileV1)::Union{Vector{String}, Missing} = ["Fe"] function tracerInjectionTime(f::MDFFile)::Union{Vector{DateTime}, Nothing} p = typeof(f) <: MDFFileV1 ? "/tracer/time" : "/tracer/injectionTime" time = @keyoptional f[p] - if isnothing(nothing) + if isnothing(time) return nothing end From f6a217daec0878294397d934b3123e772ae0b69c Mon Sep 17 00:00:00 2001 From: "Justin Ackers @ mpi05" Date: Mon, 8 Jan 2024 16:21:46 +0100 Subject: [PATCH 23/23] remove Test from main deps (make aqua happy) --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index 2bbe24ca..6794e1e9 100644 --- a/Project.toml +++ b/Project.toml @@ -29,7 +29,6 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" Tar = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" UnitfulAngles = "6fb2a4bd-7999-5318-a3b2-8ad61056cd98"