diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ada4a1f3..16ffdd7b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -15,7 +15,7 @@ jobs: fail-fast: false matrix: version: - - '1.9' + - 'lts' - '1' #- 'nightly' os: diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 00000000..2532ce7c --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,48 @@ +# Changelog + +## 0.6.0 + +### Most Important Breaking Changes +- **Breaking**: There have been multiple changes to arguments in the Scanner.toml related to both devices and the scanner, please find a detailed migration guide [https://magneticparticleimaging.github.io/MPIMeasurements.jl/dev/config/upgrade.html#v0.5-to-v0.6](here) +- **Breaking**: the receive transfer function is now defined per receive channel instead of per scanner, this allows the sequence to flexibly select receive channels and assemble the correct TF + +### Improved support for arbitrary waveform components +- define an `ArbitraryElectricalComponent` by using an amplitude, phase and base waveform (`values`) +- `values` can either be a vector or the filename to an .h5 file with field "/values" located in the new `Waveforms` folder of the Scanner +- added TX controller for arbitrary waveform components (see next section) + +### Updated TxDAQController and Feedback +- completely reworked internals of `TxDAQController` device +- added new `ControlSequence` type structure to implement a tx controller to control arbitrary waveform and DC enabled channels, the type of ControlSequence to be used is automatically detected based on the requirements of the sequence that should be controlled, the old behaviour is implemented as `CrossCouplingControlSequence` +- split `amplitudeAccuracy` and `fieldToVolDeviation` settings into relative and absolute values to improve flexibility, the two conditions are combined as OR +- `phaseAccuracy` has a unit now +- added caching of last control values to increase control speed of repeating measurements +- feedback calibration is now handled as a complex valued, optionally frequency dependent transfer function +- forward calibration of tx channels can now be a complex number to include a phase shift +- removed `correctCrossCoupling` setting from TxDAQControllerParams, if any field sets decouple=true the controller will try to decouple it + + +### New MPS Measurement Protocol +Updated the `MPSMeasurementProtocol` used to measure hybrid system matrix measurements in an MPS. +New features include: +- offsets are now defined as a `ProtocolOffsetElectricalChannel` directly in the sequence instead of the protocol +- added a wait time per offset channel, to account for slow DC sources. Any data recorded during this settling time will be discarded +- added functionality to RedPitayaDAQ channels to use H-bridges for switching the polarity of DC offsets +- new algorithm ordering the channels to reduce total wait time +- save data in proper system matrix format for hybrid reconstruction + +### New MultiSequenceSystemMatrixProtocol +Measures a (hybrid) system matrix that is defined by one `Sequence` per position, this can be used to flexibly vary any component of the field sequence like amplitudes, phases and offsets. The individual measurements will be saved together as frames in a joint MDF file. Between the individual measurements the system can be instructed to wait until the value of a temperature sensor is below a configurable threshhold. + +### General Updates +- the phase of signal components can now also be one of {"cosine", "cos", "sine", "sin", "-cosine", "-cos", "-sine", "-sin"} instead of giving the phase directly +- a magnetic field that needs to be decoupled will also require control +- frequency dividers can now be rational, the trajectory length will still be an integer using the lcm of the numerators +- all field amplitudes can now be given as a voltage, circumventing field control +- add `block` keyword argument to `startProtocol` of the ConsoleProtocolHandler to only return from the function when the protocol is finished +- devices now have a config file parameter containing the file from which they were initialized +- re-included implementation for fiber optical temperature sensor (FOTemp.jl) +- small fixes regarding different Isel robot versions +- renamed the `saveAsSystemMatrix` parameter to `saveInCalibFolder` +- removed `defaultSequence` parameter from scanner as the sequence is always defined in the protocol +- small updates to documentation \ No newline at end of file diff --git a/Project.toml b/Project.toml index 9897b3e9..5d571f00 100644 --- a/Project.toml +++ b/Project.toml @@ -1,12 +1,13 @@ name = "MPIMeasurements" uuid = "a874a27a-e9a7-503d-98b6-bf61df4bb725" authors = ["Tobias Knopp "] -version = "0.5.0" +version = "0.6.0" [deps] DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" LibSerialPort = "a05a14c7-6e3b-5ba9-90a2-45558833e1df" @@ -14,17 +15,20 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MPIFiles = "371237a9-e6c1-5201-9adb-3d8cfa78fa9f" MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" Mmap = "a63ad114-7e13-5084-954f-fe012c677804" +NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" RedPitayaDAQServer = "c544963a-496b-56d4-a5fe-f99a3f174c8f" -Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +RelocatableFolders = "05181044-ff0b-4ac5-8273-598c1e38db00" ReplMaker = "b873ce64-0db9-51f5-a568-4457d8e49576" Scratch = "6c6a2e73-6563-6170-7368-637461726353" Sockets = "6462fe0b-24de-5631-8697-dd941f90decc" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StringEncodings = "69024149-9ee7-55f6-a4c4-859efe599b68" TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" ThreadPools = "b189fb0b-2eb5-4ed4-bc0c-d34c51242431" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" +UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] @@ -32,19 +36,21 @@ Aqua = "0.8" DataStructures = "0.18, 0.19" Dates = "1" DocStringExtensions = "0.8, 0.9" +FFTW = "1.8" HDF5 = "0.16, 0.17" InteractiveUtils = "1" LibSerialPort = "0.5.2" LinearAlgebra = "1" -MPIFiles = "0.14.2, 0.15, 0.16, 0.17" +MPIFiles = "0.15,0.16,0.17" MacroTools = "0.5.6" Mmap = "1" +NaturalSort = "1" Pkg = "1" Plots = "1.39" ProgressMeter = "1.5" REPL = "1" -RedPitayaDAQServer = "0.6, 0.7, 0.8, 0.10" -Reexport = "0.2, 1.0" +RedPitayaDAQServer = "0.6, 0.7, 0.8, 0.9, 0.10, 0.11" +RelocatableFolders = "1" ReplMaker = "0.2.7" Scratch = "1.1.0" Sockets = "1" @@ -54,8 +60,9 @@ TOML = "1" Test = "1" ThreadPools = "2.1.1" UUIDs = "1" +UnicodePlots = "3" Unitful = "1.13, 1.14, 1.15, 1.16, 1.17" -julia = "1.9" +julia = "1.10" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" diff --git a/config/SimpleRedPitayaScanner/Scanner.toml b/config/SimpleRedPitayaScanner/Scanner.toml index 19464ec9..71491fd4 100644 --- a/config/SimpleRedPitayaScanner/Scanner.toml +++ b/config/SimpleRedPitayaScanner/Scanner.toml @@ -7,7 +7,6 @@ topology = "FFL" gradient = "42T/m" #datasetStore = "/opt/data/HeadScanner" datasetStore = "~/.mpi/Data" -defaultSequence = "1DSequence" defaultProtocol = "MPIMeasurement" producerThreadID = 2 consumerThreadID = 3 diff --git a/config/SimpleSimulatedScanner/Scanner.toml b/config/SimpleSimulatedScanner/Scanner.toml index b3d7505c..3d683896 100644 --- a/config/SimpleSimulatedScanner/Scanner.toml +++ b/config/SimpleSimulatedScanner/Scanner.toml @@ -6,7 +6,6 @@ name = "SimpleSimulatedScanner" topology = "FFL" # Of course ;) gradient = "42T/m" # Note: This is later parsed with Unitful datasetStore = "/home/foerger/Documents/DatasetStores/store1" -defaultSequence = "1DSequence" defaultProtocol = "MPIMeasurement" producerThreadID = 2 consumerThreadID = 3 diff --git a/docs/make.jl b/docs/make.jl index 4525aafe..add0db6e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,18 +7,25 @@ using Unitful makedocs( sitename = "MPIMeasurements", authors = "Tobias Knopp et al.", - format = Documenter.HTML(prettyurls = false), + format = Documenter.HTML(prettyurls = false, size_threshold = 500000,), modules = [MPIMeasurements], pages = [ "Home" => "index.md", "Installation" => "installation.md", - "Framework" => Any[ + "Framework Explanations" => Any[ "Scanner" => "framework/scanner.md", "Devices" => "framework/devices.md", "Sequences" => "framework/sequences.md", "Protocols" => "framework/protocols.md", "Examples" => "framework/examples.md", ], + "Configuration Files / Parameters" => Any[ + "Upgrade Guide" => "config/upgrade.md", + "Scanner" => "config/scanner.md", + "Devices" => "config/devices.md", + "Sequences" => "config/sequence.md", + "Protocols" => "config/protocols.md", + ], "Library" => Any[ "Framework" => Any[ "Scanner" => "lib/framework/scanner.md", @@ -26,20 +33,16 @@ makedocs( "Sequence" => "lib/framework/sequence.md", "Protocol" => "lib/framework/protocol.md", ], - "Base" => Any[ - "Devices" => Any[ - "Robots" => Any[ - "Interface" => "lib/base/devices/robots/interface.md", - "Isel" => "lib/base/devices/robots/isel.md" - ], - "Virtual" => Any[ - "Serial Port Pool" => "lib/base/devices/virtual/serialportpool.md", - ] + + "Devices" => Any[ + "Robots" => Any[ + "Interface" => "lib/base/devices/robots/interface.md", + "Isel" => "lib/base/devices/robots/isel.md" ], - "Protocols" => Any[ - "MPIMeasurement" => "lib/base/protocols/mpimeasurement.md" + "Virtual" => Any[ + "Serial Port Pool" => "lib/base/devices/virtual/serialportpool.md", ] - ] + ], ], ], warnonly = [:docs_block, :autodocs_block, :cross_references], diff --git a/docs/src/config/devices.md b/docs/src/config/devices.md new file mode 100644 index 00000000..a4f4d0c2 --- /dev/null +++ b/docs/src/config/devices.md @@ -0,0 +1,13 @@ +# Device Parameters +## DAQ +```@docs +MPIMeasurements.RedPitayaDAQParams +MPIMeasurements.RedPitayaLUTChannelParams +MPIMeasurements.DAQRxChannelParams +MPIMeasurements.DAQTxChannelParams +``` + +## TxController +```@docs +MPIMeasurements.TxDAQControllerParams +``` \ No newline at end of file diff --git a/docs/src/config/protocols.md b/docs/src/config/protocols.md new file mode 100644 index 00000000..e4b990cd --- /dev/null +++ b/docs/src/config/protocols.md @@ -0,0 +1,6 @@ +# Protocols + +```@autodocs +Modules = [MPIMeasurements] +Filter = t -> typeof(t) === DataType && t <: MPIMeasurements.ProtocolParams +``` \ No newline at end of file diff --git a/docs/src/config/scanner.md b/docs/src/config/scanner.md new file mode 100644 index 00000000..f066b832 --- /dev/null +++ b/docs/src/config/scanner.md @@ -0,0 +1,4 @@ +# Scanner Parameters +```@docs +MPIMeasurements.MPIScannerGeneral +``` \ No newline at end of file diff --git a/docs/src/config/sequence.md b/docs/src/config/sequence.md new file mode 100644 index 00000000..47347807 --- /dev/null +++ b/docs/src/config/sequence.md @@ -0,0 +1,17 @@ +# Sequence Parameters +## Magnetic Field +```@docs +MagneticField +``` + +## Channels +```@autodocs +Modules = [MPIMeasurements] +Filter = t -> typeof(t) === DataType && t <: MPIMeasurements.TxChannel +``` + +## Components +```@autodocs +Modules = [MPIMeasurements] +Filter = t -> typeof(t) === DataType && t <: MPIMeasurements.ElectricalComponent +``` \ No newline at end of file diff --git a/docs/src/config/upgrade.md b/docs/src/config/upgrade.md new file mode 100644 index 00000000..ac784e48 --- /dev/null +++ b/docs/src/config/upgrade.md @@ -0,0 +1,37 @@ +# Upgrade Guide +This page should give you hints on what to change in your configuration files when upgrading between breaking releases of this package. This guide focusses on removed or changed parameters not on new features that have a sensible default value, for that please check the release notes. + +## v0.5 to v0.6 +### Scanner and Devices +- the `transferFunction` parameter is no longer in the [General] section, but is now a per-channel parameter attached to each `DAQRxChannelParams`. Use `transferFunction = "tfFile.h5:2"` to select channel 2 out of the file containing multilpe channels +- the `feedback` group is removed from the tx channels and replaced by two other parameters: `feedback.channelID` is replaced by `feedbackChannelID` in the tx channel directly and `feedback.calibration` is replaced by the `transferFunction` parameter in the corresponding receive channel. The `transferFunction` parameter can correctly parse a single (complex) value with units + +### TxDAQController + +**Old parameter:** \ +`amplitudeAccuracy` \ +**Replaced by:** \ +`absoluteAmplitudeAccuracy`, absolute control accuracy threshhold with units of magnetic field, e.g. "50µT" \ +`relativeAmplitudeAcccuracy`, defined as the allowable deviation of the amplitude as a ratio of the desired amplitude, e.g. 0.001 + +**Old parameter:** \ +`fieldToVoltDeviation` \ +**Replaced by:** \ +`fieldToVoltAbsDeviation`, absolute threshhold for how much the actual field amplitude is allowed to vary from the expected value in units of magnetic field, to still accept the system as safe, e.g. "5mT" +`fieldToVoltRelDeviation`, relative threshhold for allowed deviation +*Values will be used as rtol and atol of `isapprox`* + +**Old parameter:** \ +`controlPause` \ +**Replaced by:** \ +`timeUntilStable`, time in s to wait before evaluating the feedback signals after ramping + +**Changed parameters:** \ +`phaseAccuracy`, is now a Unitful value, specify e.g. "0.1°" + +**Removed parameters:** \ +`correctCrossCoupling`, if a field has `decouple = true` the controller will correct the cross coupling + + +#### Calibrations +Because the feedback calibration (or transfer function) is now a complex value it is possible to include the phase shift between feedback and field into this number. This allows the field phase in the Sequence.toml to be set to the desired/nominal value ("0.0rad" or "sin" for sine excitation and "pi/2*rad" or "cos" for cosine excitation) instead of correcting for the feedback phase shift using that phase value. \ No newline at end of file diff --git a/docs/src/framework/devices.md b/docs/src/framework/devices.md index 521d121d..96652792 100644 --- a/docs/src/framework/devices.md +++ b/docs/src/framework/devices.md @@ -20,7 +20,7 @@ The fields of each `Device` can be grouped into three parts. ### Common Device Fields Every `Device` must have the fields `deviceID, params, present` and `dependencies`, as these are the fields used during automatic instantiation of a `Scanner`. -The `deviceID` is the name of a specific `Device` instance and corresponds to the name/key used in the `Scanner.toml`. The `params` field contains the user provided configuration parameters (see [Device Parameter](@ref)). +The `deviceID` is the name of a specific `Device` instance and corresponds to the name/key used in the `Scanner.toml`. The `params` field contains the user provided configuration parameters (see [Device Parameter Field](@ref)). Lastly the `present` field denotes if a `Device` was succesfully initialized and the `dependencies` field contains a `Dict` containing all the dependent `Devices` of the current `Device`. diff --git a/docs/src/framework/scanner.md b/docs/src/framework/scanner.md index 8d726f3b..edeb5189 100644 --- a/docs/src/framework/scanner.md +++ b/docs/src/framework/scanner.md @@ -34,7 +34,7 @@ The `Scanner.toml` contains the configuration parameters of the `Scanner` and is ### General Section The general section contains the details of the scanner. The fields correspond to the `scanner` group of the [MPI data format (MDF)](https://github.com/MagneticParticleImaging/MDF) and are used when writing a measurement to disk. -All fields that have units will be parsed with [Unitful](https://github.com/PainterQubits/Unitful.jl) and should therefore be denoted as strings with the unit attached without a space. This also applies to the [Devices section](@ref). +All fields that have units will be parsed with [Unitful](https://github.com/PainterQubits/Unitful.jl) and should therefore be denoted as strings with the unit attached without a space. This also applies to the [Device Section](@ref). ```toml [General] @@ -46,20 +46,17 @@ topology = "" gradient = "XXT/m" ``` -### Runtime Section - -The protocol section contains hints for `Protocol`, scripts and GUI implementations, such as on which Julia threads to run certain threads or which default `Sequences` to display. +In addition, the General section contains hints for `Protocol`, scripts and GUI implementations, such as on which Julia threads to run certain threads or which default `Sequences` to display. ```toml -[Runtime] defaultProtocol = "" -transferFunction = "" -dataStorage = " +datasetStore = " # Which Julia threads to run common tasks on +producerThreadID = 1 protocolThreadID = 2 -producerThreadID = 3 -consumerThreadID = 4 -serialThreadID = 5 + +consumerThreadID = 3 +serialThreadID = 4 ``` ### Device Section diff --git a/docs/src/framework/sequences.md b/docs/src/framework/sequences.md index 15a849f5..1c6f39b3 100644 --- a/docs/src/framework/sequences.md +++ b/docs/src/framework/sequences.md @@ -14,6 +14,7 @@ The general section contains a description string for the sequence, as well as t ```toml [General] +name = "" description = "" targetScanner = "" baseFrequency = "125MHz" diff --git a/docs/src/lib/base/protocols/mpimeasurement.md b/docs/src/lib/base/protocols/mpimeasurement.md deleted file mode 100644 index fc036e26..00000000 --- a/docs/src/lib/base/protocols/mpimeasurement.md +++ /dev/null @@ -1,5 +0,0 @@ -# MPIMeasurementProtocol - -```@docs -MPIMeasurements.MPIMeasurementProtocolParams -``` \ No newline at end of file diff --git a/docs/src/lib/framework/scanner.md b/docs/src/lib/framework/scanner.md index 01f8135c..29fa52f1 100644 --- a/docs/src/lib/framework/scanner.md +++ b/docs/src/lib/framework/scanner.md @@ -11,7 +11,6 @@ MPIMeasurements.scannerName MPIMeasurements.scannerTopology MPIMeasurements.scannerGradient MPIMeasurements.scannerDatasetStore -MPIMeasurements.defaultSequence MPIMeasurements.defaultProtocol ``` diff --git a/src/Devices/DAQ/DAQ.jl b/src/Devices/DAQ/DAQ.jl index 46859ebb..e9ab1e7e 100644 --- a/src/Devices/DAQ/DAQ.jl +++ b/src/Devices/DAQ/DAQ.jl @@ -1,11 +1,11 @@ import Base: setindex!, getindex -export AbstractDAQ, DAQParams, SinkImpedance, SINK_FIFTY_OHM, SINK_HIGH, DAQTxChannelSettings, DAQChannelParams, DAQFeedback, DAQTxChannelParams, DAQRxChannelParams, +export AbstractDAQ, DAQParams, SinkImpedance, SINK_FIFTY_OHM, SINK_HIGH, DAQTxChannelSettings, DAQChannelParams, DAQTxChannelParams, DAQRxChannelParams, createDAQChannels, createDAQParams, startTx, stopTx, setTxParams, readData, numRxChannelsTotal, numTxChannelsTotal, numRxChannelsActive, numTxChannelsActive, currentPeriod, getDAQ, getDAQs, channelIdx, limitPeak, sinkImpedance, allowedWaveforms, isWaveformAllowed, - feedbackChannelID, feedbackCalibration, calibration + feedbackChannelID, feedbackTransferFunction, transferFunction, hasTransferFunction, calibration abstract type AbstractDAQ <: Device end abstract type DAQParams <: DeviceParams end @@ -15,6 +15,13 @@ abstract type DAQParams <: DeviceParams end SINK_HIGH end +@enum TxValueRange begin + #POSITIVE + #NEGATIVE + BOTH + HBRIDGE +end + struct DAQTxChannelSettings "Applied channel voltage. Dimensions are (components, channels, periods)." amplitudes::Array{typeof(1.0u"V"), 3} @@ -50,55 +57,111 @@ abstract type DAQChannelParams end abstract type TxChannelParams <: DAQChannelParams end abstract type RxChannelParams <: DAQChannelParams end -Base.@kwdef struct DAQFeedback - channelID::AbstractString - calibration::Union{typeof(1.0u"T/V"), Nothing} = nothing +Base.@kwdef struct DAQHBridge{N} + channelID::Union{String, Vector{String}} + manual::Bool = false + deadTime::typeof(1.0u"s") = 0.0u"s" + level::Matrix{typeof(1.0u"V")} # Can this be simplified by a convention for h-bridges? end +negativeLevel(bridge::DAQHBridge) = bridge.level[:, 1] +positiveLevel(bridge::DAQHBridge) = bridge.level[:, 2] +level(bridge::DAQHBridge, x::Number) = signbit(x) ? negativeLevel(bridge) : positiveLevel(bridge) +manual(bridge::DAQHBridge) = bridge.manual +deadTime(bridge::DAQHBridge) = bridge.deadTime +id(bridge::DAQHBridge) = bridge.channelID + +function createDAQChannels(::Type{DAQHBridge}, dict::Dict{String, Any}) + splattingDict = Dict{Symbol, Any}() + splattingDict[:channelID] = dict["channelID"] isa Vector ? dict["channelID"] : [dict["channelID"]] + N = length(splattingDict[:channelID]) + splattingDict[:level] = reshape(uparse.(dict["level"]), N, :) + + if haskey(dict, "manual") + splattingDict[:manual] = dict["manual"] + end -Base.@kwdef struct DAQTxChannelParams <: TxChannelParams + if haskey(dict, "deadTime") + splattingDict[:deadTime] = uparse(dict["deadTime"]) + end + return DAQHBridge{N}(;splattingDict...) +end + +""" +Parameters for a general DAQTxChannel + +$FIELDS +""" +Base.@kwdef mutable struct DAQTxChannelParams <: TxChannelParams channelIdx::Int64 limitPeak::typeof(1.0u"V") + limitSlewRate::typeof(1.0u"V/s") = 1000.0u"V/µs" # default is basically no limit sinkImpedance::SinkImpedance = SINK_HIGH allowedWaveforms::Vector{Waveform} = [WAVEFORM_SINE] - feedback::Union{DAQFeedback, Nothing} = nothing - calibration::Union{typeof(1.0u"V/T"), Nothing} = nothing + feedbackChannelID::Union{String, Nothing} = nothing + calibration::Union{TransferFunction, String, Nothing} = nothing end -Base.@kwdef struct DAQRxChannelParams <: RxChannelParams +""" +Parameters for a general DAQRxChannel + +$FIELDS +""" +Base.@kwdef mutable struct DAQRxChannelParams <: RxChannelParams + "required, Integer, channel idx of the corresponding DAQ" channelIdx::Int64 + "optional, can be either a filename to a .h5 file in the TransferFunctions folder or a single (complex) number" + transferFunction::Union{TransferFunction, String, Nothing} = nothing end -"Create DAQ channel description from device dict part." -function createDAQChannels(dict::Dict{String, Any}) - channels = Dict{String, DAQChannelParams}() - for (key, value) in dict - splattingDict = Dict{Symbol, Any}() - if value["type"] == "tx" - splattingDict[:channelIdx] = value["channel"] - splattingDict[:limitPeak] = uparse(value["limitPeak"]) +function createDAQChannel(::Type{DAQTxChannelParams}, dict::Dict{String,Any}) + splattingDict = Dict{Symbol, Any}() + splattingDict[:channelIdx] = dict["channel"] + splattingDict[:limitPeak] = uparse(dict["limitPeak"]) - if haskey(value, "sinkImpedance") - splattingDict[:sinkImpedance] = value["sinkImpedance"] == "FIFTY_OHM" ? SINK_FIFTY_OHM : SINK_HIGH - end + if haskey(dict, "limitSlewRate") + splattingDict[:limitSlewRate] = uparse(dict["limitSlewRate"]) + end - if haskey(value, "allowedWaveforms") - splattingDict[:allowedWaveforms] = toWaveform.(value["allowedWaveforms"]) - end + if haskey(dict, "sinkImpedance") + splattingDict[:sinkImpedance] = dict["sinkImpedance"] == "FIFTY_OHM" ? SINK_FIFTY_OHM : SINK_HIGH + end + + if haskey(dict, "allowedWaveforms") + splattingDict[:allowedWaveforms] = toWaveform.(dict["allowedWaveforms"]) + end - if haskey(value, "feedback") - channelID=value["feedback"]["channelID"] - calibration=uparse(value["feedback"]["calibration"]) + if haskey(dict, "feedbackChannelID") + splattingDict[:feedbackChannelID] = dict["feedbackChannelID"] + end - splattingDict[:feedback] = DAQFeedback(channelID=channelID, calibration=calibration) - end + if haskey(dict, "calibration") + splattingDict[:calibration] = parse_into_tf(dict["calibration"]) + end - if haskey(value, "calibration") - splattingDict[:calibration] = uparse.(value["calibration"]) - end + return DAQTxChannelParams(;splattingDict...) +end + +function createDAQChannel(::Type{DAQRxChannelParams}, dict::Dict{String,Any}) + splattingDict = Dict{Symbol, Any}() + splattingDict[:channelIdx] = dict["channel"] + + if haskey(dict, "transferFunction") + splattingDict[:transferFunction] = parse_into_tf(dict["transferFunction"]) + end - channels[key] = DAQTxChannelParams(;splattingDict...) + return DAQRxChannelParams(;splattingDict...) +end + +"Create DAQ channel description from device dict part." +function createDAQChannels(dict::Dict{String, Any}) + channels = Dict{String, DAQChannelParams}() + for (key, value) in dict + if value["type"] == "tx" + channels[key] = createDAQChannel(DAQTxChannelParams, value) elseif value["type"] == "rx" - channels[key] = DAQRxChannelParams(channelIdx=value["channel"]) + channels[key] = createDAQChannel(DAQRxChannelParams, value) + elseif value["type"] == "txSlow" + channels[key] = createDAQChannel(RedPitayaLUTChannelParams, value) end end @@ -128,7 +191,7 @@ function createDAQParams(DAQType::Type{T}, dict::Dict{String, Any}) where {T <: splattingDict = dict_to_splatting(mainDict) splattingDict[:channels] = createDAQChannels(DAQType, channelDict) - + # TODO: check if DAQ type can actually support all types of channels try return DAQType(;splattingDict...) catch e @@ -173,6 +236,7 @@ end @mustimplement numComponentsMax(daq::AbstractDAQ) @mustimplement canPostpone(daq::AbstractDAQ) @mustimplement canConvolute(daq::AbstractDAQ) +@mustimplement channel(daq::AbstractDAQ, channelID::AbstractString) getDAQs(scanner::MPIScanner) = getDevices(scanner, AbstractDAQ) getDAQ(scanner::MPIScanner) = getDevice(scanner, AbstractDAQ) @@ -184,17 +248,131 @@ getDAQ(scanner::MPIScanner) = getDevice(scanner, AbstractDAQ) # return factor # end -channel(daq::AbstractDAQ, channelID::AbstractString) = daq.params.channels[channelID] -channelIdx(daq::AbstractDAQ, channelID::AbstractString) = channel(daq, channelID).channelIdx -channelIdx(daq::AbstractDAQ, channelIDs::Vector{<:AbstractString}) = [channel(daq, channelID).channelIdx for channelID in channelIDs] -limitPeak(daq::AbstractDAQ, channelID::AbstractString) = channel(daq, channelID).limitPeak -sinkImpedance(daq::AbstractDAQ, channelID::AbstractString) = channel(daq, channelID).sinkImpedance -allowedWaveforms(daq::AbstractDAQ, channelID::AbstractString) = channel(daq, channelID).allowedWaveforms -isWaveformAllowed(daq::AbstractDAQ, channelID::AbstractString, waveform::Waveform) = waveform in allowedWaveforms(daq, channelID) -feedback(daq::AbstractDAQ, channelID::AbstractString) = channel(daq, channelID).feedback -feedbackChannelID(daq::AbstractDAQ, channelID::AbstractString) = feedback(daq, channelID).channelID -feedbackCalibration(daq::AbstractDAQ, channelID::AbstractString) = feedback(daq, channelID).calibration -calibration(daq::AbstractDAQ, channelID::AbstractString) = channel(daq, channelID).calibration +channelIdx(channel::DAQChannelParams) = channel.channelIdx +channelIdx(daq::AbstractDAQ, channelID::AbstractString) = channelIdx(channel(daq, channelID)) +channelIdx(daq::AbstractDAQ, channelIDs::Vector{<:AbstractString}) = [channelIdx(channel(daq, channelID)) for channelID in channelIDs] + +limitPeak(channel::DAQTxChannelParams) = channel.limitPeak +limitPeak(daq::AbstractDAQ, channelID::AbstractString) = limitPeak(channel(daq, channelID)) + +sinkImpedance(channel::DAQTxChannelParams) = channel.sinkImpedance +sinkImpedance(daq::AbstractDAQ, channelID::AbstractString) = sinkImpedance(channel(daq, channelID)) + +allowedWaveforms(channel::DAQTxChannelParams) = channel.allowedWaveforms +allowedWaveforms(daq::AbstractDAQ, channelID::AbstractString) = allowedWaveforms(channel(daq, channelID)) + +isWaveformAllowed(channel::DAQTxChannelParams, waveform::Waveform) = waveform in allowedWaveforms(channel) +isWaveformAllowed(daq::AbstractDAQ, channelID::AbstractString, waveform::Waveform) = isWaveformAllowed(channel(daq, channelID), waveform) + +feedbackChannelID(channel::DAQTxChannelParams) = channel.feedbackChannelID +feedbackChannelID(daq::AbstractDAQ, channelID::AbstractString) = feedbackChannelID(channel(daq, channelID)) + +feedbackChannel(daq::AbstractDAQ, channel_::DAQTxChannelParams) = channel(daq, feedbackChannelID(channel_)) +feedbackChannel(daq::AbstractDAQ, channelID::AbstractString) = feedbackChannel(daq, channel(daq,channelID)) + +feedbackTransferFunction(daq::AbstractDAQ, channel::DAQTxChannelParams) = transferFunction(daq, feedbackChannel(daq, channel)) +feedbackTransferFunction(daq::AbstractDAQ, channelID::AbstractString) = feedbackTransferFunction(daq, channel(daq, channelID)) + +hasTransferFunction(channel::DAQRxChannelParams) = !isnothing(channel.transferFunction) +hasTransferFunction(daq::AbstractDAQ, channelID::AbstractString) = hasTransferFunction(channel(daq,channelID)) + +function splitTFName(name::AbstractString) + if endswith(name, r"\.h5:\d+") + tmp = rsplit(name, ":"; limit=2) + return tmp[1], parse(Int, tmp[2]) + else + return name, nothing + end +end +transferFunction(::AbstractDAQ, ::Nothing) = nothing +transferFunction(daq::AbstractDAQ, channelID::AbstractString) = transferFunction(daq, channel(daq, channelID)) +function transferFunction(dev::Union{MPIScanner, AbstractDAQ}, channel::DAQRxChannelParams) + if isa(channel.transferFunction, String) + filename, tfChannel = splitTFName(channel.transferFunction) + channel.transferFunction = TransferFunction(joinpath(configDir(dev), "TransferFunctions", filename), channels = tfChannel) + else + channel.transferFunction + end +end + +calibration(daq::AbstractDAQ, channelID::AbstractString) = calibration(daq, channel(daq,channelID)) +function calibration(dev::Union{MPIScanner, AbstractDAQ}, channel::DAQTxChannelParams) + if isa(channel.calibration, String) + filename, tfChannel = splitTFName(channel.calibration) + channel.calibration = TransferFunction(joinpath(configDir(dev), "TransferFunctions", filename), channels = tfChannel) + else + channel.calibration + end +end + +calibration(dev::Device, channelID::AbstractString, frequencies) = calibration.([dev], [channelID], frequencies) +calibration(daq::AbstractDAQ, channelID::AbstractString, frequency::Real) = calibration(daq, channel(daq, channelID), frequency) +calibration(dev::Union{MPIScanner, AbstractDAQ}, channel::DAQTxChannelParams, frequency::Real) = calibration(dev, channel)(frequency) + +export applyForwardCalibration!, applyForwardCalibration + +function applyForwardCalibration(seq::Sequence, device::Device) + seqCopy = deepcopy(seq) + applyForwardCalibration!(seqCopy, device) + return seqCopy +end + +function applyForwardCalibration!(seq::Sequence, device::Device) + + for channel in periodicElectricalTxChannels(seq) + off = offset(channel) + if dimension(off) != dimension(1.0u"V") + isnothing(calibration(device, id(channel))) && throw(ScannerConfigurationError("An offset value in channel $(id(channel)) requires calibration but no calibration is configured on the DAQ channel!")) + offsetVolts = off*abs(calibration(device, id(channel), 0)) # use DC value for offsets + offset!(channel, uconvert(u"V",offsetVolts)) + end + + for comp in components(channel) + amp = amplitude(comp) + pha = phase(comp) + if dimension(amp) != dimension(1.0u"V") + isnothing(calibration(device, id(channel))) && throw(ScannerConfigurationError("An amplitude value in channel $(id(channel)) requires calibration but no calibration is configured on the DAQ channel!")) + f_comp = ustrip(u"Hz", txBaseFrequency(seq)) / divider(comp) + complex_comp = (amp*exp(im*pha)) * calibration(device, id(channel), f_comp) + amplitude!(comp, uconvert(u"V",abs(complex_comp))) + phase!(comp, angle(complex_comp)u"rad") + if comp isa ArbitraryElectricalComponent + N = length(values(comp)) + f_awg = rfftfreq(N, f_comp*N) + calib = calibration(device, id(channel), f_awg) ./ (abs.(calibration(device, id(channel), f_comp))*exp.(im*2*pi*range(0,length(f_awg)-1).*angle(calibration(device, id(channel),f_comp)))) # since amplitude and phase are already calibrated for the base frequency, here we need to remove that factor + calib = ustrip.(NoUnits, calib) + values!(comp, irfft(rfft(values(comp)).*calib, N)) + end + end + end + end + + for lutChannel in acyclicElectricalTxChannels(seq) + if lutChannel isa StepwiseElectricalChannel + values = lutChannel.values + if dimension(values[1]) != dimension(1.0u"V") + isnothing(calibration(device, id(lutChannel))) && throw(ScannerConfigurationError("A value in channel $(id(lutChannel)) requires calibration but no calibration is configured on the DAQ channel!")) + values = values.*calibration(device, id(lutChannel)) + lutChannel.values = values + end + elseif lutChannel isa ContinuousElectricalChannel + amp = lutChannel.amplitude + off = lutChannel.offset + if dimension(amp) != dimension(1.0u"V") + isnothing(calibration(device, id(lutChannel))) && throw(ScannerConfigurationError("An amplitude value in channel $(id(lutChannel)) requires calibration but no calibration is configured on the DAQ channel!")) + amp = amp*calibration(device, id(lutChannel)) + lutChannel.amplitude = amp + end + if dimension(off) != dimension(1.0u"V") + isnothing(calibration(device, id(lutChannel))) && throw(ScannerConfigurationError("An offset value in channel $(id(lutChannel)) requires calibration but no calibration is configured on the DAQ channel!")) + off = off*calibration(device, id(lutChannel)) + lutChannel.offset = off + end + end + end + + nothing +end diff --git a/src/Devices/DAQ/RedPitayaDAQ.jl b/src/Devices/DAQ/RedPitayaDAQ.jl index fca3fe17..548dd3ec 100644 --- a/src/Devices/DAQ/RedPitayaDAQ.jl +++ b/src/Devices/DAQ/RedPitayaDAQ.jl @@ -2,6 +2,11 @@ export RampingMode, NONE, HOLD, STARTUP @enum RampingMode NONE HOLD STARTUP export RedPitayaDAQParams +""" +Parameters for a DAQ of type `RedPitayaDAQ` + +$(FIELDS) +""" Base.@kwdef mutable struct RedPitayaDAQParams <: DAQParams "All configured channels of this DAQ device." channels::Dict{String, DAQChannelParams} @@ -21,74 +26,59 @@ Base.@kwdef mutable struct RedPitayaDAQParams <: DAQParams counterTriggerSourceChannel::DIOPins = DIO7_P end -Base.@kwdef struct RedPitayaTxChannelParams <: TxChannelParams - channelIdx::Int64 - limitPeak::typeof(1.0u"V") - sinkImpedance::SinkImpedance = SINK_HIGH - allowedWaveforms::Vector{Waveform} = [WAVEFORM_SINE] - feedback::Union{DAQFeedback, Nothing} = nothing - calibration::Union{typeof(1.0u"V/T"), typeof(1.0u"V/A"), Nothing} = nothing -end +""" +Parameters for a TxChannel of type RedPitayaLUTChannel +$(FIELDS) +""" Base.@kwdef struct RedPitayaLUTChannelParams <: TxChannelParams channelIdx::Int64 calibration::Union{typeof(1.0u"V/T"), typeof(1.0u"V/A"), Nothing} = nothing + range::TxValueRange = BOTH + hbridge::Union{DAQHBridge, Nothing} = nothing + switchTime::typeof(1.0u"s") = 0.0u"s" + switchEnable::Bool = true end -"Create the params struct from a dict. Typically called during scanner instantiation." -function RedPitayaDAQParams(dict::Dict{String, Any}) - return createDAQParams(RedPitayaDAQParams, dict) -end +calibration(dev::Device, channel::RedPitayaLUTChannelParams) = channel.calibration -function createDAQChannels(::Type{RedPitayaDAQParams}, dict::Dict{String, Any}) - # TODO This is mostly copied from createDAQChannels, maybe manage to get rid of the duplication - channels = Dict{String, DAQChannelParams}() - for (key, value) in dict - splattingDict = Dict{Symbol, Any}() - if value["type"] == "tx" - splattingDict[:channelIdx] = value["channel"] - splattingDict[:limitPeak] = uparse(value["limitPeak"]) +function createDAQChannel(::Type{RedPitayaLUTChannelParams}, dict::Dict{String, Any}) + splattingDict = Dict{Symbol, Any}() + splattingDict[:channelIdx] = dict["channel"] + if haskey(dict, "calibration") + calib = uparse.(dict["calibration"]) - if haskey(value, "sinkImpedance") - splattingDict[:sinkImpedance] = value["sinkImpedance"] == "FIFTY_OHM" ? SINK_FIFTY_OHM : SINK_HIGH - end - - if haskey(value, "allowedWaveforms") - splattingDict[:allowedWaveforms] = toWaveform.(value["allowedWaveforms"]) - end + if unit(upreferred(calib)) == upreferred(u"V/T") + calib = calib .|> u"V/T" + elseif unit(upreferred(calib)) == upreferred(u"V/A") + calib = calib .|> u"V/A" + else + error("The values have to be either given as a V/T or in V/A. You supplied the type `$(eltype(calib))`.") + end + splattingDict[:calibration] = calib + end - if haskey(value, "feedback") - channelID=value["feedback"]["channelID"] - calibration=uparse(value["feedback"]["calibration"]) + if haskey(dict, "hbridge") + splattingDict[:hbridge] = createDAQChannels(DAQHBridge, dict["hbridge"]) + end - splattingDict[:feedback] = DAQFeedback(channelID=channelID, calibration=calibration) - end + if haskey(dict, "range") + splattingDict[:range] = dict["range"] + end - if haskey(value, "calibration") - splattingDict[:calibration] = uparse.(value["calibration"]) - end + if haskey(dict, "switchTime") + splattingDict[:switchTime] = uparse(dict["switchTime"]) + end - channels[key] = RedPitayaTxChannelParams(;splattingDict...) - elseif value["type"] == "rx" - channels[key] = DAQRxChannelParams(channelIdx=value["channel"]) - elseif value["type"] == "txSlow" - calib = nothing - if haskey(value, "calibration") - calib = uparse.(value["calibration"]) - - if unit(upreferred(calib)) == upreferred(u"V/T") - calib = calib .|> u"V/T" - elseif unit(upreferred(calib)) == upreferred(u"V/A") - calib = calib .|> u"V/A" - else - error("The values have to be either given as a V/t or in V/A. You supplied the type `$(eltype(calib))`.") - end - end - channels[key] = RedPitayaLUTChannelParams(channelIdx=value["channel"], calibration = calib) - end + if haskey(dict, "switchEnable") + splattingDict[:switchEnable] = dict["switchEnable"] end + return RedPitayaLUTChannelParams(; splattingDict...) +end - return channels +"Create the params struct from a dict. Typically called during scanner instantiation." +function RedPitayaDAQParams(dict::Dict{String, Any}) + return createDAQParams(RedPitayaDAQParams, dict) end export RedPitayaDAQ @@ -120,6 +110,15 @@ Base.@kwdef mutable struct RedPitayaDAQ <: AbstractDAQ end function _init(daq::RedPitayaDAQ) + # force all calibration transfer functions to be read from disk + for channelID in keys(daq.params.channels) + if isa(channel(daq, channelID), DAQTxChannelParams) + calibration(daq, channelID) + elseif isa(channel(daq, channelID), DAQRxChannelParams) + transferFunction(daq, channelID) + end + end + # Restart the DAQ if necessary try daq.rpc = RedPitayaCluster(daq.params.ips, triggerMode = daq.params.triggerMode) @@ -162,6 +161,8 @@ optionalDependencies(::RedPitayaDAQ) = [TxDAQController, SurveillanceUnit] Base.close(daq::RedPitayaDAQ) = daq.rpc +channel(daq::RedPitayaDAQ, channelID::AbstractString) = daq.params.channels[channelID] + export usesCounterTrigger usesCounterTrigger(daq::RedPitayaDAQ) = daq.params.useCounterTrigger @@ -200,7 +201,7 @@ function setRampingParams(daq::RedPitayaDAQ, sequence::Sequence) for channel in txChannels m = nothing idx = channel[2].channelIdx - if channel[2] isa RedPitayaTxChannelParams + if channel[2] isa DAQTxChannelParams m = idx elseif channel[2] isa RedPitayaLUTChannelParams # Map to fast DAC @@ -308,9 +309,8 @@ function setSequenceParams(daq::RedPitayaDAQ, luts::Vector{Union{Nothing, Array{ @debug "There are no LUTs to set." end - @info "Set sequence params" - stepsPerRepetition = div(periodsPerFrame(daq.rpc), daq.acqPeriodsPerPatch) + @debug "Set sequence params" samplesPerStep=div(samplesPerPeriod(daq.rpc) * periodsPerFrame(daq.rpc), stepsPerRepetition) result = execute!(daq.rpc) do batch @add_batch batch samplesPerStep!(daq.rpc, div(samplesPerPeriod(daq.rpc) * periodsPerFrame(daq.rpc), stepsPerRepetition)) @add_batch batch clearSequence!(daq.rpc) @@ -380,12 +380,7 @@ function createLUT!(lut::Array{Float32}, start, channelMapping) lutValues = [] lutIdx = [] for (lutChannel, seqChannel) in channelMapping - tempValues = values(seqChannel) - if !isnothing(lutChannel.calibration) - tempValues = tempValues.*lutChannel.calibration - end - tempValues = ustrip.(u"V", tempValues) - push!(lutValues, tempValues) + push!(lutValues, ustrip.(u"V", values(seqChannel))) push!(lutIdx, lutChannel.channelIdx) end @@ -444,6 +439,408 @@ function getTiming(daq::RedPitayaDAQ) return sampleTiming end + +function prepareProtocolSequences(base::Sequence, daq::RedPitayaDAQ; numPeriodsPerOffset::Int64 = 1) + cpy = deepcopy(base) + + # Prepare offset sequences + offsetVector = ProtocolOffsetElectricalChannel[] + fieldMap = Dict{ProtocolOffsetElectricalChannel, MagneticField}() + calibsize = Int64[] + for field in cpy + for channel in channels(field, ProtocolOffsetElectricalChannel) + push!(offsetVector, channel) + fieldMap[channel] = field + push!(calibsize, length(values(channel))) + end + end + + stepfreq = 125e6/lcm(dfDivider(cpy)) + if stepfreq > 10e3 + if stepfreq/numPeriodsPerOffset > 10e3 + error("One period of the configured drivefield is only $(round(1/stepfreq*1e6,digits=4)) us long. To ensure that the RP can keep up with the offsets, please use a number of periods per offset that is a multiple of 50 and at least $(ceil(stepfreq/10e3/50)*50)!") + else + if (numPeriodsPerOffset/50 != numPeriodsPerOffset÷50) + error("One period of the configured drivefield is only $(round(1/stepfreq*1e6,digits=4)) us long. To ensure that the RP can keep up with the offsets, please use a number of periods per offset that is a multiple of 50!") + else + # if the DF period is shorter than 0.1ms (10kHz) we use 50 DF periods for every step, this allows us to use DF periods of up to 10kHz*50=500kHz and results in a worst case wait time inefficiency of 0.1ms*50 = 5ms + # it would be poossible to use other numbers for periodsPerStep, but this will probably work best for both extremes + periodsPerStep = 50 + stepfreq = stepfreq/periodsPerStep + end + end + else + # if a DF period is longer than 0.1ms we can just use a single step for every period, ensuring most efficient use of waittimes + periodsPerStep = 1 + end + + stepsPerOffset = numPeriodsPerOffset÷periodsPerStep + + # Compute step duration + stepduration = (1/stepfreq) + + @debug "prepareProtocolSequences: Stepduration: $stepduration, Stepfreq: $stepfreq" + allSteps = nothing + enables = nothing + hbridges = nothing + if any(x -> requireHBridge(x, daq), offsetVector) + offsets, allSteps, enables, hbridges, permutation = prepareHSwitchedOffsets(offsetVector, daq, cpy, stepduration) + else + offsets, allSteps, enables, hbridges, permutation = prepareOffsets(offsetVector, daq, cpy, stepduration) + end + + numOffsets = length(first(allSteps)[2]) + df = lcm(dfDivider(cpy)) + + numMeasOffsets = length(permutation) + numWaitOffsets = numOffsets-numMeasOffsets + numPeriodsPerFrame = (numMeasOffsets * stepsPerOffset + numWaitOffsets) * periodsPerStep + + divider = df * numPeriodsPerFrame + + # offsets and permutation refer to offsets without multiple df per period and include every offset only once + # We need to perform two operations: + # 1. build an index mask that handles the repetition of offsets, such that offsets[mask] contains the desired repetitions of offsets (not the wait times) + stepRepetition = collect(1:numOffsets) + if numPeriodsPerOffset > 1 + + stepRepetition = Int[] + for i in 1:numOffsets + if i in permutation # if offset i is a offset that should be saved + append!(stepRepetition, repeat([i],stepsPerOffset)) # we have to repeat it by the amount of steps per Offset + else + push!(stepRepetition, i) # otherwise we just output it once + end + end + + # apply step repetition to permutation, every valid step is repeated stepsPerOffset times + updatedPermutation = repeat(permutation, inner=stepsPerOffset) + # to remove the repeated values, add 1 to all indizes starting from each value that is identical to its previous neighbor + sortedIdx = sortperm(updatedPermutation) + for i=2:length(updatedPermutation) + if updatedPermutation[sortedIdx[i]] == updatedPermutation[sortedIdx[i-1]] + updatedPermutation[sortedIdx[i:end]] .+= 1 + end + end + + offsets = offsets[stepRepetition,:] + + permutation = Int64[] + for patch in updatedPermutation + idx = (patch-1)*periodsPerStep+1:patch*periodsPerStep + append!(permutation, idx) + end + + end + @debug "prepareProtocolSequences: Permutation after update" permutation offsets + + # Generate actual StepwiseElectricalChannel + for (ch, steps) in allSteps + stepwise = StepwiseElectricalChannel(id = id(ch), divider = divider, values = identity.(steps)[stepRepetition], enable = enables[ch][stepRepetition]) + fieldMap[ch][id(stepwise)] = stepwise + + # Generate H-Bridge channels + if haskey(hbridges, ch) + offsetChannel = channel(daq, id(ch)) + hbridgeChannels = id(offsetChannel.hbridge) + hbridgeValues = hbridges[ch] + for (i, hbridgeChannel) in enumerate(hbridgeChannels) + hsteps = map(x -> x[i], hbridgeValues) + hbridgeStepwise = StepwiseElectricalChannel(id = hbridgeChannel, divider = divider, values = hsteps[stepRepetition], enable = fill(true, length(stepRepetition))) + fieldMap[ch][hbridgeChannel] = hbridgeStepwise + end + end + end + + return cpy, permutation, offsets, calibsize, numPeriodsPerFrame +end + +function prepareHSwitchedOffsets(offsetVector::Vector{ProtocolOffsetElectricalChannel}, daq::RedPitayaDAQ, seq::Sequence, stepduration) + switchingIndices = findall(x->requireHBridge(x,daq), offsetVector) + switchingChannel = offsetVector[switchingIndices] + regularChannel = isempty(switchingChannel) ? offsetVector : offsetVector[filter(!in(switchingIndices), 1:length(offsetVector))] + allSteps = Dict{ProtocolOffsetElectricalChannel, Vector{Any}}() + allEnables = Dict{ProtocolOffsetElectricalChannel, Vector{Union{Bool, Missing}}}() + validSteps = Union{Bool, Missing}[] + + # Prepare values with missing for h-bridge switch + # Permute patches by considering quadrants. Use lower bits of number to encode permutations + # for example 011 means channel 1 is disabled, channel 2 and 3 h-bridge is enabled + for comb in 0:2^length(switchingChannel)-1 + quadrant = map(Bool, digits(comb, base = 2, pad=length(switchingChannel))) + + offsets = Vector{Vector{Any}}() + + # Same order here as in Iterators.flatten below + for ch in regularChannel + push!(offsets, values(ch)) + end + for (i, isPositive) in enumerate(quadrant) + push!(offsets, values(switchingChannel[i], isPositive)) + end + + orderedChannel = collect(Iterators.flatten((regularChannel, switchingChannel))) + + # For each quadrant we can generate the offsets individually + quadrantSteps, couldPause, anySwitching, perm = prepareOffsetSwitches(offsets, orderedChannel, daq, stepduration) + + for (i, ch) in enumerate(orderedChannel[perm]) + steps = get(allSteps, ch, []) + enables = get(allEnables, ch, Union{Bool, Missing}[]) + push!(steps, quadrantSteps[i]...) + push!(enables, map(x-> channel(daq, id(ch)).switchEnable ? true : !x, couldPause[i])...) + + # If not at the end then add missing, here H-bridge switches will be inserted later + if comb != 2^length(switchingChannel)-1 + push!(steps, missing) + push!(enables, missing) + end + + allSteps[ch] = steps + allEnables[ch] = enables + end + + # Steps are valid if no channel has a switch pause + push!(validSteps, map(!, anySwitching)...) + if comb != 2^length(switchingChannel)-1 + push!(validSteps, missing) + end + end + + # Compute switch timing + deadTimes = map(x-> x.hbridge.deadTime, filter(x-> x.range == HBRIDGE, map(x->channel(daq, id(x)), offsetVector))) + maxTime = maximum(map(ustrip, deadTimes)) + numSwitchPeriods = Int64(ceil(maxTime/stepduration)) + + hbridges = prepareHBridgeLevels(allSteps, daq, numSwitchPeriods) + + # Set enable to false during hbridge switching + for (ch, enables) in allEnables + temp = Bool[] + foreach(x-> ismissing(x) ? push!(temp, fill(false, numSwitchPeriods)...) : push!(temp, x), enables) + allEnables[ch] = temp + end + + # Likewise set valid to false during hbridge switching + tempValid = Bool[] + foreach(x-> ismissing(x) ? push!(tempValid, fill(false, numSwitchPeriods)...) : push!(tempValid, x), validSteps) + validSteps = tempValid + + # Prepare H-Bridge pauses + for (ch, steps) in allSteps + temp = [] + # Add numSwitchPeriods 0 for every missing value + foreach(x-> ismissing(x) ? push!(temp, fill(zero(eltype(steps[1])), numSwitchPeriods)...) : push!(temp, x), steps) + allSteps[ch] = temp + end + + # Generate offsets without permutation + offsets = map(values, offsetVector) + offsets = hcat(prepareOffsets(offsets, daq)...) + # Map each offset combination to its original index + offsetDict = Dict{Vector, Int64}() + for (i, row) in enumerate(eachrow(offsets)) + offsetDict[row] = i + end + + # Permute offsets + # Sort patches according to their value in the index dictionary + permoffsets = hcat(map(x-> allSteps[x], offsetVector)...) + patchPermutation = sortperm(map(pair -> validSteps[pair[1]] ? offsetDict[identity(pair[2])] : typemax(Int64), enumerate(eachrow(permoffsets)))) + # Pause/switching patches are sorted to the end and need to be removed + patchPermutation = patchPermutation[1:end-length(filter(!identity, validSteps))] + + # Remove (negative) sign from all channels with an h-bridge + for (ch, steps) in filter(x->haskey(hbridges, x[1]), allSteps) + allSteps[ch] = abs.(steps) + end + + return permoffsets, allSteps, allEnables, hbridges, patchPermutation +end + +function prepareOffsets(offsetVector::Vector{ProtocolOffsetElectricalChannel}, daq::RedPitayaDAQ, seq::Sequence, stepduration) + allSteps = Dict{ProtocolOffsetElectricalChannel, Vector{Any}}() + enables = Dict{ProtocolOffsetElectricalChannel, Vector{Bool}}() + + offsets = map(values, offsetVector) + steps, couldPause, anySwitching, perm = prepareOffsetSwitches(offsets, offsetVector, daq, stepduration) + validSteps = map(!, anySwitching) + for (i, ch) in enumerate(offsetVector[perm]) + allSteps[ch] = steps[i] + enables[ch] = map(x-> channel(daq, id(ch)).switchEnable ? true : !x, couldPause[i]) + end + + # Map each offset combination to its original index + offsets = hcat(prepareOffsets(offsets, daq)...) + offsetDict = Dict{Vector, Int64}() + for (i, row) in enumerate(eachrow(offsets)) + offsetDict[row] = i + end + + # Permute offsets + # Sort patches according to their value in the index dictionary + permoffsets = hcat(map(x-> allSteps[x], offsetVector)...) + patchPermutation = sortperm(map(pair -> validSteps[pair[1]] ? offsetDict[identity(pair[2])] : typemax(Int64), enumerate(eachrow(permoffsets)))) + # Switching patches are sorted to the end and need to be removed + patchPermutation = patchPermutation[1:end-length(filter(!identity, validSteps))] + + hbridges = prepareHBridgeLevels(allSteps, daq) + # Remove (negative) sign from all channels with an h-bridge + for (ch, steps) in filter(x->haskey(hbridges, x[1]), allSteps) + allSteps[ch] = abs.(steps) + end + + return permoffsets, allSteps, enables, hbridges, patchPermutation +end + +function prepareHBridgeLevels(allSteps, daq::RedPitayaDAQ, switchSteps::Int64 = 0) + hbridges = Dict{ProtocolOffsetElectricalChannel, Vector{Any}}() + for ch in filter(ch -> channel(daq, id(ch)).range == HBRIDGE , keys(allSteps)) + offsetChannel = channel(daq, id(ch)) + steps = allSteps[ch] + hsteps = [] + for (i, step) in enumerate(steps) + if ismissing(step) # Assumes only one missing per switch and no missing at last index + push!(hsteps, fill(level(offsetChannel.hbridge, steps[i + 1]), switchSteps)...) + else + push!(hsteps, level(offsetChannel.hbridge, step)) + end + end + hbridges[ch] = hsteps + end + return hbridges +end + +function requireHBridge(offsetChannel::ProtocolOffsetElectricalChannel{T}, daq::RedPitayaDAQ) where T + offsets = values(offsetChannel) + ch = channel(daq, id(offsetChannel)) + if ch.range == BOTH + return false + elseif ch.range == HBRIDGE + # Assumption protocol offset channel can only produce one sign change + # Consider 0 to be "positive" + first = offsets[1] + last = offsets[end] + return signbit(first) != signbit(last) && (!iszero(first) && !iszero(last)) + else + error("Unexpected channel range") + end +end + +function prepareOffsetSwitches(offsets::Vector{Vector{T}}, channels::Vector{ProtocolOffsetElectricalChannel}, daq::RedPitayaDAQ, stepduration::Float64) where T + switchSteps = Dict{ProtocolOffsetElectricalChannel, Int64}() + for ch in channels + daqChannel = channel(daq, id(ch)) + if !iszero(daqChannel.switchTime) + switchSteps[ch] = Int64(ceil(ustrip(u"s", daqChannel.switchTime/stepduration))) + end + end + + if isempty(switchSteps) + offsets = prepareOffsets(offsets, daq) + return offsets, [fill(false, length(first(offsets))) for ch in channels], fill(false, length(first(offsets))), 1:length(channels) + end + + # Sorting is not optimal, but allows us in the next step to consider pauses individually and not have to take max of all possible pauses at a point + perm = sortperm(map(x-> haskey(switchSteps, x) ? switchSteps[x] : 0, channels)) + + + # The idea is to realize the following pattern for n-channel (shown with three: x, y, z) + # A switch S_x means the next "proper" offset value is to be repeated S_x times, the amount of switch steps x requires + # x_i means the i'th value of the x offsets + # ()_i means loop the pattern in the bracket i times + + # i = number of x offsets, j = # y offsets , k = # z offsets + # x: (S_z - S_y (S_y - S_ x (S_x, x_i)_i)_j)_k + # x: (S_z - S_y (S_y - S_ x (S_x, y_j)_i)_j)_k + # x: (S_z - S_y (S_y - S_ x (S_x, z_k)_i)_j)_k + + sortedOffsets = offsets[perm] + sortedChannels = channels[perm] + sortedSwitchSteps = map(ch -> get(switchSteps, ch, 0), sortedChannels) + + sortedOffsetsWithSwitches = [] + sortedCouldPause = [] + anyChannelSwitching = Bool[] + + # Each channel will add the same switch "pauses" + # This array contains the precomputed S_x, S_y - S_x, ..., + sortedLoopSteps = zeros(Int64, size(sortedSwitchSteps)) + previous = 0 + for (i, pause) in enumerate(sortedSwitchSteps) + sortedLoopSteps[i] = pause - previous + previous = pause + end + + # The pattern will be constructed with a "nested" for loop. + # However, the nesting depth depends on the number of channels. To make this generic to n-channel we will use CartesianIndices + iterations = map(length, sortedOffsets) + loop = CartesianIndices(Tuple(iterations)) + + # Here we want to have a loop as follows (but generic for n-channels): + # for ch in channels + # for k = 1:length(z) + # for j = 1:length(y) + # for i = 1:length(x) + # # ... create pattern + for (channelIdx, channel) in enumerate(sortedChannels) + offsetVec = sortedOffsets[channelIdx] + resultOffsets = eltype(offsetVec)[] # Store the offsets of the current channel + resultCouldPause = Bool[] # Store when the current channel doesn't need to send (i.e others switch and he doesn't yet) + resultAnySwitching = Bool[] # Store when any channel is switching + + previous = CartesianIndex(Tuple(fill(-1, length(iterations)))) # This will be used to detect which index is currently "running"/changing + for indices in loop + value = offsetVec[indices[channelIdx]] + + # Apply (loop) switch steps for each changing index, i.e. an index that is not equal to its previous iteration value + requiredSwitches = map(!, iszero.(Tuple(indices - previous))) + + # Compute how many switching "pauses" are to be inserted + switches = 0 + switchIndices = findall(requiredSwitches) + if !isempty(switchIndices) + switches = sum(sortedLoopSteps[findall(requiredSwitches)]) + end + # If the current channel is switching, we need to subtract that from couldPause + ownSwitches = in(channelIdx, switchIndices) ? sortedSwitchSteps[channelIdx] : 0 + + # + 1 is for the actual offset + push!(resultOffsets, fill(value, 1 + switches)...) + push!(resultCouldPause, vcat(fill(true, switches - ownSwitches), fill(false, 1 + ownSwitches))...) + + # anySwitching only needs to be computed once + if isempty(anyChannelSwitching) + push!(resultAnySwitching, push!(fill(true, switches), false)...) + end + + previous = indices + end + + push!(sortedOffsetsWithSwitches, resultOffsets) + push!(sortedCouldPause, resultCouldPause) + if isempty(anyChannelSwitching) + anyChannelSwitching = resultAnySwitching + end + end + + return sortedOffsetsWithSwitches, sortedCouldPause, anyChannelSwitching, perm +end + +function prepareOffsets(offsetVectors::Vector{Vector{T}}, daq::RedPitayaDAQ) where T + result = Vector{Vector{Any}}() + inner = 1 + numOffsets = prod(length.(offsetVectors)) + for offsets in offsetVectors + outer = div(numOffsets, inner * length(offsets)) + steps = repeat(offsets, inner = inner, outer = outer) + inner*=length(offsets) + push!(result, steps) + end + return map(x -> identity.(x), result) # Improve typing from Vector{Vector{Any}} +end + + #### Producer/Consumer #### mutable struct RedPitayaAsyncBuffer <: AsyncBuffer samples::Union{Matrix{Int16}, Nothing} @@ -506,7 +903,7 @@ function startProducer(channel::Channel, daq::RedPitayaDAQ, numFrames; isControl # Start pipeline @debug "Pipeline started" try - @debug currentWP(daq.rpc) + @debug "Producer:" currentWP(daq.rpc) readSamples(rpu, startSample, samplesToRead, channel, chunkSize = chunkSize) catch e @info "Attempting reconnect to reset pipeline" @@ -573,18 +970,22 @@ end function setup(daq::RedPitayaDAQ, sequence::Sequence) stopTx(daq) setupRx(daq, sequence) - setupTx(daq, sequence) - prepareTx(daq, sequence) - setSequenceParams(daq, sequence) + sequenceVolt = applyForwardCalibration(sequence, daq) + setupTx(daq, sequenceVolt) + prepareTx(daq, sequenceVolt) + setSequenceParams(daq, sequenceVolt) end function setupTx(daq::RedPitayaDAQ, sequence::Sequence) @debug "Setup tx" + periodicChannels = periodicElectricalTxChannels(sequence) if any([length(component.amplitude) > 1 for channel in periodicChannels for component in periodicElectricalComponents(channel)]) error("The Red Pitaya DAQ cannot work with more than one period in a frame or frequency sweeps yet.") end + + @debug "SetupTx: Outputting the following offsets" offsets=[offset(chan) for chan in periodicChannels] # Iterate over sequence(!) channels execute!(daq.rpc) do batch @@ -601,11 +1002,19 @@ function setupTxChannel(daq::RedPitayaDAQ, channel::PeriodicElectricalChannel, b end function setupTxChannel!(batch::ScpiBatch, daq::RedPitayaDAQ, channel::PeriodicElectricalChannel, baseFreq) channelIdx_ = channelIdx(daq, id(channel)) # Get index from scanner(!) channel - offsetVolts = offset(channel)*calibration(daq, id(channel)) - @add_batch batch offsetDAC!(daq.rpc, channelIdx_, ustrip(u"V", offsetVolts)) - for (idx, component) in enumerate(components(channel)) + @add_batch batch offsetDAC!(daq.rpc, channelIdx_, ustrip(u"V", offset(channel))) + + for (idx, component) in enumerate(periodicElectricalComponents(channel)) # in the future add the SweepElectricalComponents as well setupTxComponent!(batch, daq, channel, component, idx, baseFreq) end + + awgComps = arbitraryElectricalComponents(channel) + if length(awgComps) > 1 + throw(SequenceConfigurationError("The channel with the ID $(id(channel)) defines more than one arbitrary electrical component, which is not supported.")) + elseif length(awgComps) == 1 + setupTxComponent!(batch, daq, channel, awgComps[1], 0, baseFreq) + end + end function setupTxComponent(daq::RedPitayaDAQ, channel::PeriodicElectricalChannel, component, componentIdx, baseFreq) @@ -615,7 +1024,7 @@ function setupTxComponent(daq::RedPitayaDAQ, channel::PeriodicElectricalChannel, end function setupTxComponent!(batch::ScpiBatch, daq::RedPitayaDAQ, channel::PeriodicElectricalChannel, component::Union{PeriodicElectricalComponent, SweepElectricalComponent}, componentIdx, baseFreq) if componentIdx > 3 - throw(SequenceConfigurationError("The channel with the ID `$(id(channel))` defines more than three fixed waveform components, this is more than what the RedPitaya offers. Use an arbitrary waveform is a fourth is necessary.")) + throw(SequenceConfigurationError("The channel with the ID `$(id(channel))` defines more than three fixed waveform components, this is more than what the RedPitaya offers. Use an arbitrary waveform if a fourth is necessary.")) end channelIdx_ = channelIdx(daq, id(channel)) # Get index from scanner(!) channel freq = ustrip(u"Hz", baseFreq) / divider(component) @@ -623,7 +1032,7 @@ function setupTxComponent!(batch::ScpiBatch, daq::RedPitayaDAQ, channel::Periodi waveform_ = uppercase(fromWaveform(waveform(component))) if !isWaveformAllowed(daq, id(channel), waveform(component)) throw(SequenceConfigurationError("The channel with the ID `$(id(channel))` "* - "defines a waveforms of $waveform_, but the scanner channel does not allow this.")) + "defines a waveform of $waveform_, but the scanner channel does not allow this.")) end @add_batch batch signalTypeDAC!(daq.rpc, channelIdx_, componentIdx, waveform_) end @@ -635,7 +1044,7 @@ function setupTxComponent!(batch::ScpiBatch, daq::RedPitayaDAQ, channel::Periodi waveform_ = uppercase(fromWaveform(waveform(component))) if !isWaveformAllowed(daq, id(channel), waveform(component)) throw(SequenceConfigurationError("The channel with the ID `$(id(channel))` "* - "defines a waveforms of $waveform_, but the scanner channel does not allow this.")) + "defines a waveform of $waveform_, but the scanner channel does not allow this.")) end # Waveform is set together with amplitudes for arbitrary waveforms end @@ -651,9 +1060,9 @@ function setupRx(daq::RedPitayaDAQ, sequence::Sequence) @assert txBaseFrequency(sequence) == 125.0u"MHz" "The base frequency is fixed for the Red Pitaya "* "and must thus be 125 MHz and not $(txBaseFrequency(sequence))." - # The decimation can only be a power of 2 beginning with 8 + # The decimation has to be divisible by 2 and must be 8 or larger decimation_ = upreferred(txBaseFrequency(sequence)/rxSamplingRate(sequence)) - if decimation_ in [2^n for n in 3:8] + if iseven(decimation_) && decimation_ >= 8 daq.decimation = decimation_ else throw(ScannerConfigurationError("The decimation derived from the rx bandwidth of $(rxBandwidth(sequence)) and "* @@ -672,14 +1081,14 @@ function setupRx(daq::RedPitayaDAQ, sequence::Sequence) # TODO possibly move some of this into abstract daq daq.refChanIDs = [] - txChannels = [channel[2] for channel in daq.params.channels if channel[2] isa RedPitayaTxChannelParams] - daq.refChanIDs = unique([tx.feedback.channelID for tx in txChannels if !isnothing(tx.feedback)]) + txChannels = [channel[2] for channel in daq.params.channels if channel[2] isa DAQTxChannelParams] + daq.refChanIDs = unique(filter!(!isnothing, feedbackChannelID.(txChannels))) # Construct view to save bandwidth rxIDs = sort(union(channelIdx(daq, daq.rxChanIDs), channelIdx(daq, daq.refChanIDs))) selection = [false for i = 1:length(daq.rpc)] for i in map(x->div(x -1, 2) + 1, rxIDs) - @debug i + @debug "setupRx: Adding RP id $i to cluster view" selection[i] = true end daq.rpv = RedPitayaClusterView(daq.rpc, selection) @@ -717,7 +1126,9 @@ function startTx(daq::RedPitayaDAQ; isControlStep=false) end function stopTx(daq::RedPitayaDAQ) - masterTrigger!(daq.rpc, false) + if masterTrigger(daq.rpc) + masterTrigger!(daq.rpc, false) # only deactivate trigger if it is currently active + end execute!(daq.rpc) do batch @add_batch batch serverMode!(daq.rpc, CONFIGURATION) for channel in 1:2*length(daq.rpc) @@ -738,6 +1149,7 @@ function clearTx!(daq::RedPitayaDAQ) for comp = 1:4 @add_batch batch amplitudeDAC!(daq.rpc, channel, comp, 0.0) end + @add_batch batch offsetDAC!(daq.rpc, channel, 0.0) end end end @@ -749,38 +1161,31 @@ end function prepareTx(daq::RedPitayaDAQ, sequence::Sequence) stopTx(daq) - @debug "Preparing amplitude and phase" allAmps = Dict{String, Vector{typeof(1.0u"V")}}() allPhases = Dict{String, Vector{typeof(1.0u"rad")}}() allAwg = Dict{String, Vector{typeof(1.0u"V")}}() - for channel in periodicElectricalTxChannels(sequence) + for channel in periodicElectricalTxChannels(sequence) name = id(channel) - amps = [] - phases = [] - wave = nothing - for comp in periodicElectricalComponents(channel) - # Lengths check == 1 happens in setupTx already - amp = amplitude(comp) - if dimension(amp) == dimension(1.0u"T") - amp = (amp * calibration(daq, name)) - end - push!(amps, amp) - push!(phases, phase(comp)) + # This should set all unused components to zero, but what about unused channels? + amps = zeros(typeof(1.0u"V"), 3) + phases = zeros(typeof(1.0u"rad"), 3) + wave = zeros(typeof(1.0u"V"), 2^14) + for (i,comp) in enumerate(periodicElectricalComponents(channel)) + # Length check < 3 happens in setupTx already + amps[i] = amplitude(comp) + phases[i] = phase(comp) end + # Length check < 1 happens in setupTx already comps = arbitraryElectricalComponents(channel) - if length(comps) > 2 - throw(ScannerConfigurationError("Channel $(id(channel)) defines more than one arbitrary electrical component, which is not supported.")) - elseif length(comps) == 1 - wave = values(comps[1]) - if dimension(wave[1]) != dimension(1.0u"V") - wave = wave.*calibration(daq, name) - end - allAwg[name] = wave + if length(comps)==1 + wave = scaledValues(comps[1]) end + allAwg[name] = wave allAmps[name] = amps allPhases[name] = phases end + @debug "prepareTx: Outputting the following amplitudes and phases:" allAmps allPhases awg=lineplot(1:2^14,ustrip.(u"V",hcat(collect(Base.values(allAwg))...)),ylabel="AWG / V", name=collect(Base.keys(allAwg)), canvas=DotCanvas, border=:ascii, height=10, xlim=(1,2^14)) setTxParams(daq, allAmps, allPhases, allAwg) end @@ -909,7 +1314,7 @@ numTxChannelsActive(daq::RedPitayaDAQ) = numChan(daq.rpc) #TODO: Currently, all numRxChannelsActive(daq::RedPitayaDAQ) = numRxChannelsReference(daq)+numRxChannelsMeasurement(daq) numRxChannelsReference(daq::RedPitayaDAQ) = length(daq.refChanIDs) numRxChannelsMeasurement(daq::RedPitayaDAQ) = length(daq.rxChanIDs) -numComponentsMax(daq::RedPitayaDAQ) = 3 +numComponentsMax(daq::RedPitayaDAQ) = 4 canPostpone(daq::RedPitayaDAQ) = true canConvolute(daq::RedPitayaDAQ) = false diff --git a/src/Devices/DAQ/SimpleSimulatedDAQ.jl b/src/Devices/DAQ/SimpleSimulatedDAQ.jl index 942b1fff..e81f1eef 100644 --- a/src/Devices/DAQ/SimpleSimulatedDAQ.jl +++ b/src/Devices/DAQ/SimpleSimulatedDAQ.jl @@ -54,6 +54,7 @@ end neededDependencies(::SimpleSimulatedDAQ) = [SimulationController] optionalDependencies(::SimpleSimulatedDAQ) = [TxDAQController, SurveillanceUnit] +channel(daq::SimpleSimulatedDAQ, channelID::AbstractString) = daq.params.channels[channelID] Base.close(daq::SimpleSimulatedDAQ) = nothing @@ -90,8 +91,8 @@ function setupTx(daq::SimpleSimulatedDAQ, sequence::Sequence) channelMapping = scannerChannel.channelIdx # Switch to getter? # Activate corresponding receive channels - if !isnothing(scannerChannel.feedback) - feedbackChannelID = scannerChannel.feedback.channelID + if !isnothing(scannerChannel.feedbackChannelID) + feedbackChannelID = scannerChannel.feedbackChannelID scannerFeedbackChannel = daq.params.channels[feedbackChannelID] feedbackChannelIdx = scannerFeedbackChannel.channelIdx # Switch to getter? push!(daq.refChanIDs, feedbackChannelID) @@ -184,7 +185,7 @@ function readDataPeriods(daq::SimpleSimulatedDAQ, startPeriod::Integer, numPerio if dimension(daq.amplitude[1]) == dimension(u"T") factor = 1.0u"T/T" else - factor = 1/scannerChannel.calibration + factor = 1/scannerChannel.calibration # TODO/JA: figure out how to include change here end temperatureRise = daq.params.temperatureRise[sendChannelID] @@ -215,12 +216,12 @@ function readDataPeriods(daq::SimpleSimulatedDAQ, startPeriod::Integer, numPerio Bᵢ = Bₘₐₓ.*sin.(2π*f*t.+ϕ) # Desired, ideal field without drift Bᵣ = (Bₘₐₓ.+ΔB).*sin.(2π*f*t.+ϕ.+Δϕ) # Field with drift of phase and amplitude - uₜₓ .+= Bᵢ.*sendChannel.calibration + uₜₓ .+= Bᵢ.*sendChannel.calibration # TODO/JA: figure out how to include change here uᵣₓ .+= simulateLangevinInduced(t, Bᵣ, f, ϕ.+Δϕ) # f is not completely correct due to the phase change, but this is only a rough approximation anyways # Assumes the same induced voltage from the field as given out with uₜₓ, # just with a slight change in phase and amplitude - uᵣₑ .+= Bᵣ.*sendChannel.calibration + uᵣₑ .+= Bᵣ.*sendChannel.calibration # TODO/JA: figure out how to include change here end # Assumes one reference and one measurement channel for each send channel diff --git a/src/Devices/Device.jl b/src/Devices/Device.jl index 533f6282..635c6c0b 100644 --- a/src/Devices/Device.jl +++ b/src/Devices/Device.jl @@ -13,7 +13,7 @@ abstract type DeviceParams end Base.close(device::Device) = @warn "The device type `$(typeof(device))` has no `close` function defined." function validateDeviceStruct(device::Type{<:Device}) - requiredFields = [:deviceID, :params, :optional, :present, :dependencies] + requiredFields = [:deviceID, :params, :optional, :present, :dependencies, :configFile] missingFields = [x for x in requiredFields if !in(x, fieldnames(device))] if length(missingFields) > 0 msg = "Device struct $(string(device)) is missing the required fields " * join(missingFields, ", ", " and ") @@ -33,6 +33,8 @@ macro add_device_fields(paramType) present::Bool = false #"Vector of dependencies for this device." dependencies::Dict{String,Union{Device,Missing}} + #"Path of the config used to create the device (either Scanner.toml or Device.toml)" + configFile::Union{String, Nothing} = nothing end) end @@ -51,6 +53,22 @@ isPresent(device::Device) = device.present "Retrieve the dependencies of a device." dependencies(device::Device) = device.dependencies +"Retrieve the path to the configuration file of a device." +configFile(device::Device) = device.configFile + +"Retrieve the configuration directory of the scanner that the device is contructed from." +function configDir(device::Device) + if !isnothing(device.configFile) + if basename(device.configFile) == "Scanner.toml" + return dirname(device.configFile) # if the device is read from a Scanner.toml the directory is the scanner directory + else + return joinpath(splitpath(device.configFile)[1:end-2]...) # otherwise the Device.toml is in the Devices folder + end + else + return nothing + end +end + "Retrieve all dependencies of a certain type." dependencies(device::Device, type::DataType) = [dependency for dependency in Base.values(dependencies(device)) if dependency isa type] diff --git a/src/Devices/Devices.jl b/src/Devices/Devices.jl index 262d5ccd..0892c7bd 100644 --- a/src/Devices/Devices.jl +++ b/src/Devices/Devices.jl @@ -15,7 +15,7 @@ include("Virtual/Virtual.jl") # List our own enums to avoid accidentally converting a different enum # Did not list enums like LakeShoreF71GaussMeterConnectionModes atm, because their convert function uses specific strings # and not the enum name -for enum in [RedPitayaDAQServer.TriggerMode, RampingMode, TemperatureControlMode] +for enum in [RedPitayaDAQServer.TriggerMode, RampingMode, TemperatureControlMode, TxValueRange] @eval begin function Base.convert(::Type{$enum}, x::String) try diff --git a/src/Devices/Robots/CollisionModule/SimpleBoreCollisionModule.jl b/src/Devices/Robots/CollisionModule/SimpleBoreCollisionModule.jl index eb6b1b1d..da6abda6 100644 --- a/src/Devices/Robots/CollisionModule/SimpleBoreCollisionModule.jl +++ b/src/Devices/Robots/CollisionModule/SimpleBoreCollisionModule.jl @@ -1,5 +1,10 @@ export SimpleBoreCollisionModule, SimpleBoreCollisionModuleParams +""" +Parameters for a `SimpleBoreCollisionModule`` + +$FIELDS +""" Base.@kwdef struct SimpleBoreCollisionModuleParams <: DeviceParams "Diameter of scanner in the y-z plane" scannerDiameter::typeof(1.0u"mm") @@ -15,16 +20,7 @@ end SimpleBoreCollisionModuleParams(dict::Dict) = params_from_dict(SimpleBoreCollisionModuleParams, dict) Base.@kwdef mutable struct SimpleBoreCollisionModule <: AbstractCollisionModule3D - "Unique device ID for this device as defined in the configuration." - deviceID::String - "Parameter struct for this devices read from the configuration." - params::SimpleBoreCollisionModuleParams - "Flag if the device is optional." - optional::Bool = false - "Flag if the device is present." - present::Bool = false - "Vector of dependencies for this device." - dependencies::Dict{String,Union{Device,Missing}} + @add_device_fields SimpleBoreCollisionModuleParams end collisionModuleType(cm::SimpleBoreCollisionModule) = PositionCollisionType() diff --git a/src/Devices/Robots/IgusRobot.jl b/src/Devices/Robots/IgusRobot.jl index 348585b7..acc1836c 100644 --- a/src/Devices/Robots/IgusRobot.jl +++ b/src/Devices/Robots/IgusRobot.jl @@ -31,6 +31,11 @@ const MODES = Dict( "PROFILE_POSITION" => 1, "HOMING" => 6) +""" +Parameters for a Robot of type `IgusRobot` + +$FIELDS +""" Base.@kwdef struct IgusRobotParams <: DeviceParams defaultVelocity::Vector{typeof(1.0u"mm/s")} = [10.0u"mm/s"] axisRange::Vector{Vector{typeof(1.0u"mm")}} = [[0,500.0]]u"mm" diff --git a/src/Devices/Robots/IselRobot.jl b/src/Devices/Robots/IselRobot.jl index ad975454..44f42dd9 100644 --- a/src/Devices/Robots/IselRobot.jl +++ b/src/Devices/Robots/IselRobot.jl @@ -31,6 +31,11 @@ const iselErrorCodes = Dict( abstract type IselRobotParams <: DeviceParams end +""" +Parameters for a Robot of type `IselRobot` that is connected directly via a serial port + +$FIELDS +""" Base.@kwdef struct IselRobotPortParams <: IselRobotParams axisRange::Vector{Vector{typeof(1.0u"mm")}} = [[0,420],[0,420],[0,420]]u"mm" defaultVel::Vector{typeof(1.0u"mm/s")} = [10,10,10]u"mm/s" @@ -54,6 +59,11 @@ Base.@kwdef struct IselRobotPortParams <: IselRobotParams end IselRobotPortParams(dict::Dict) = params_from_dict(IselRobotPortParams, prepareRobotDict(dict)) +""" +Parameters for a Robot of type `IselRobot` that is connected via a serial port pool + +$FIELDS +""" Base.@kwdef struct IselRobotPoolParams <: IselRobotParams axisRange::Vector{Vector{typeof(1.0u"mm")}} = [[0,420],[0,420],[0,420]]u"mm" defaultVel::Vector{typeof(1.0u"mm/s")} = [10,10,10]u"mm/s" @@ -151,22 +161,19 @@ function _enable(robot::IselRobot) robot.enableTime = time() _enable(robot, controllerVersion(robot)) end -function _enable(robot::IselRobot, version::IseliMCS8) - # NOP +function _enable(robot::IselRobot, version) + _setMotorCurrent(robot, version, true) end -function _enable(robot::IselRobot, version::IselC142) - _setMotorCurrent(robot, true) +function _enable(robot::IselRobot, version::IseliMCS8) + _setMotorCurrent(robot, version, true) end function _disable(robot::IselRobot) writeIOOutput(robot, zeros(Bool, 8)) _disable(robot, controllerVersion(robot)) end -function _disable(robot::IselRobot, version::IseliMCS8) - # NOP -end -function _disable(robot::IselRobot, version::IselC142) - _setMotorCurrent(robot, false) +function _disable(robot::IselRobot, version) + _setMotorCurrent(robot, version, false) end @@ -304,7 +311,12 @@ macro minimumISELversion(version::Int) end) end -function _setMotorCurrent(rob::IselRobot, power::Bool) +_setMotorCurrent(rob::IselRobot, power::Bool) = _setMotorCurrent(rob, controllerVersion(rob), power) +function _setMotorCurrent(rob::IselRobot, ::IseliMCS8, power::Bool) + ret = queryIsel(rob, string("@0B3,", Int(power))) + checkIselError(ret) +end +function _setMotorCurrent(rob::IselRobot, ::IselC142, power::Bool) ret = queryIsel(rob, string("@0B65529,", Int(!power))) checkIselError(ret) end @@ -327,7 +339,7 @@ function mm2steps(rob::IselRobot, len::Vector{<:Unitful.Velocity}) return round.(Int64, temp .* rob.params.stepsPermm) end -steps2mm(rob::IselRobot, steps::Vector{Integer}) = steps2mm(rob, Int64.(steps)) +steps2mm(rob::IselRobot, steps::Vector{<:Integer}) = steps2mm(rob, Int64.(steps)) steps2mm(rob::IselRobot, steps::Vector{Int64}) = steps ./ rob.params.stepsPermm * u"mm" @@ -342,7 +354,7 @@ function setRefVelocity(rob::IselRobot, vel::Vector{<:Unitful.Velocity}) cmd = string("@0d", vel[1], ",", vel[2], ",", vel[3], ",", vel[3]) else # TODO: check if this distinction is really necessary, both should be able to use @0d - cmd = string("@0Id", " ", vel[1], ",", vel[2], ",", vel[3], ",", vel[3]) + cmd = string("@0Id", vel[1], ",", vel[2], ",", vel[3], ",", vel[3]) end ret = queryIsel(rob, cmd) checkIselError(ret) diff --git a/src/Devices/Sensors/Temperature/ArduinoTemperatureSensor.jl b/src/Devices/Sensors/Temperature/ArduinoTemperatureSensor.jl index db346565..cab27709 100644 --- a/src/Devices/Sensors/Temperature/ArduinoTemperatureSensor.jl +++ b/src/Devices/Sensors/Temperature/ArduinoTemperatureSensor.jl @@ -44,7 +44,11 @@ function _init(sensor::ArduinoTemperatureSensor) @info "Connection to ArduinoTempBox established." ard = SimpleArduino(;commandStart = params.commandStart, commandEnd = params.commandEnd, sd = sd) sensor.ard = ard - setMaximumTemps(sensor, params.maxTemps) + try + setMaximumTemps(sensor, params.maxTemps) + catch e + @warn "Temperature Sensor does not support setMaximumTemps!" error=e + end end function initSerialDevice(sensor::ArduinoTemperatureSensor, params::ArduinoTemperatureSensorPortParams) diff --git a/src/Devices/Sensors/Temperature/FOTemp.jl b/src/Devices/Sensors/Temperature/FOTemp.jl index b7b5333c..84c88f9e 100644 --- a/src/Devices/Sensors/Temperature/FOTemp.jl +++ b/src/Devices/Sensors/Temperature/FOTemp.jl @@ -1,3 +1,3 @@ -export FOTemp, gettemperature, gettemperatures, logtemperature, logtemperatures +export FOTemp include("FOTempLowLevel.jl") diff --git a/src/Devices/Sensors/Temperature/FOTempLowLevel.jl b/src/Devices/Sensors/Temperature/FOTempLowLevel.jl index ea23d2ce..296b0bd4 100644 --- a/src/Devices/Sensors/Temperature/FOTempLowLevel.jl +++ b/src/Devices/Sensors/Temperature/FOTempLowLevel.jl @@ -65,8 +65,8 @@ Returns the averaged temperature of a channel. """ function getAveragedTemperature(ft::FOTemp,channel::Char) temp = query(ft.sd, "?01 $channel") - acq = readuntil(ft.sd.sp, '\n', ft.sd.timeout_ms) - return parse(Float64,split(temp," ")[2]) / 10 + #acq = readuntil(ft.sd.sp, '\n', ft.sd.timeout_ms) + return parse(Float64,split(temp," ")[2]) / 10 end """ @@ -74,8 +74,8 @@ Returns the averaged temperature of all channel. """ function getAveragedTemperature(ft::FOTemp) temp = query(ft.sd, "?02") - acq = readuntil(ft.sd.sp, '\n', ft.sd.timeout_ms) - return parse(Float64,split(temp," ")[2]) / 10 + #acq = readuntil(ft.sd.sp, '\n', ft.sd.timeout_ms) + return parse(Float64,split(temp," ")[2]) / 10 end """ diff --git a/src/Devices/Sensors/Temperature/TemperatureSensor.jl b/src/Devices/Sensors/Temperature/TemperatureSensor.jl index e77dce78..766b98ef 100644 --- a/src/Devices/Sensors/Temperature/TemperatureSensor.jl +++ b/src/Devices/Sensors/Temperature/TemperatureSensor.jl @@ -3,7 +3,7 @@ abstract type TemperatureSensor <: Device end include("ArduinoTemperatureSensor.jl") include("DummyTemperatureSensor.jl") -#include("FOTemp.jl") +include("FOTemp.jl") include("SimulatedTemperatureSensor.jl") Base.close(t::TemperatureSensor) = nothing diff --git a/src/Devices/Virtual/SerialPortPool.jl b/src/Devices/Virtual/SerialPortPool.jl index 5b722ed9..8e269f45 100644 --- a/src/Devices/Virtual/SerialPortPool.jl +++ b/src/Devices/Virtual/SerialPortPool.jl @@ -76,6 +76,8 @@ function getSerialDevice(pool::SerialPortPool, description::String; kwargs...) @error e end end + else + throw(ScannerConfigurationError("No suitable SerialPort for `$description` found in SerialPortPool!")) end return nothing end @@ -106,7 +108,7 @@ function initSerialDevice(device::Device, query::String, response::String) pool = dependency(device, SerialPortPool) sd = getSerialDevice(pool, query, response; serial_device_splatting(device.params)...) if isnothing(sd) - throw(ScannerConfigurationError("Device $(deviceID(device)) found no fitting serial port.")) + throw(ScannerConfigurationError("Device $(deviceID(device)) failed to connect to serial port.")) end return sd else @@ -120,7 +122,7 @@ function initSerialDevice(device::Device, description::String) pool = dependency(device, SerialPortPool) sd = getSerialDevice(pool, description; serial_device_splatting(device.params)...) if isnothing(sd) - throw(ScannerConfigurationError("Device $(deviceID(device)) found no fitting serial port with description $description.")) + throw(ScannerConfigurationError("Device $(deviceID(device)) failed to connect to serial port with description $description.")) end return sd else diff --git a/src/Devices/Virtual/TxDAQController.jl b/src/Devices/Virtual/TxDAQController.jl index 85b07e65..8ec7c5a0 100644 --- a/src/Devices/Virtual/TxDAQController.jl +++ b/src/Devices/Virtual/TxDAQController.jl @@ -1,33 +1,65 @@ export TxDAQControllerParams, TxDAQController, controlTx +""" +Parameters for a `TxDAQController`` + +$FIELDS +""" Base.@kwdef mutable struct TxDAQControllerParams <: DeviceParams - phaseAccuracy::Float64 - amplitudeAccuracy::Float64 - controlPause::Float64 + "Angle, required, allowed deviation of the excitation phase" + phaseAccuracy::typeof(1.0u"rad") + "Number, required, allowed relative deviation of the excitation amplitude" + relativeAmplitudeAccuracy::Float64 + "Magnetic field, default: 50µT, allowed absolute deviation of the excitation amplitude" + absoluteAmplitudeAccuracy::typeof(1.0u"T") = 50.0u"µT" + "Integer, default: 20, maximum number of steps to try to control the system" maxControlSteps::Int64 = 20 - fieldToVoltDeviation::Float64 = 0.2 - correctCrossCoupling::Bool = false + "Bool, default: false, control the DC value of the excitation field (only posible for DC enabled DF amplifiers)" + controlDC::Bool = false + "Float, default: 0.0, time in seconds to wait before the DF is stable after ramping" + timeUntilStable::Float64 = 0.0 + "Float, default: 0.002, time in seconds that the DF should be averaged during the control measurement" minimumStepDuration::Float64 = 0.002 + "Float, default: 0.2, relative deviation allowed between forward calibration and actual system state" + fieldToVoltRelDeviation::Float64 = 0.2 + "Magnetic field, default: 5.0mT, absolute deviation allowed between forward calibration and actual system state" + fieldToVoltAbsDeviation::typeof(1.0u"T") = 5.0u"mT" + "Magnetic field, default: 40mT, maximum field amplitude that the controller should allow" + maxField::typeof(1.0u"T") = 40.0u"mT" end TxDAQControllerParams(dict::Dict) = params_from_dict(TxDAQControllerParams, dict) struct SortedRef end struct UnsortedRef end -struct ControlledChannel - seqChannel::PeriodicElectricalChannel - daqChannel::TxChannelParams -end +abstract type ControlSequence end -mutable struct ControlSequence +macro ControlSequence_fields() + return esc(quote + #"Goal sequence of the controller, amplitudes are in T" targetSequence::Sequence + #"Current best guess of a sequence to reach the goal, amplitudes are in V" currSequence::Sequence - # Periodic Electric Components - simpleChannel::OrderedDict{PeriodicElectricalChannel, TxChannelParams} + #"Dictionary containing the sequence channels needing control and their respective DAQ channels" + controlledChannelsDict::OrderedDict{PeriodicElectricalChannel, TxChannelParams} + #"Vector of indices into the receive channels selecting the correct feedback channels for the controlledChannels" + refIndices::Vector{Int64} + #"Vector of TransferFunctions to be applied to the feedback channels" + refTFs::Vector{TransferFunction} + end) +end + +mutable struct CrossCouplingControlSequence <: ControlSequence + @ControlSequence_fields sinLUT::Union{Matrix{Float64}, Nothing} cosLUT::Union{Matrix{Float64}, Nothing} - refIndices::Vector{Int64} - # Arbitrary Waveform +end + +mutable struct AWControlSequence <: ControlSequence + @ControlSequence_fields + rfftIndices::BitArray{3} # Matrix of size length(controlledChannelsDict) x 4 (max. num of components) x len(rfft) + dcSearch::Vector{@NamedTuple{V::Vector{Float64}, B::Vector{Float64}}} + #lutMap::Dict{String, Dict{AcyclicElectricalTxChannel, Int}} # for every field ID contain a dict of removed LUTChannels together with the corresponding channel end @enum ControlResult UNCHANGED UPDATED INVALID @@ -35,133 +67,316 @@ end Base.@kwdef mutable struct TxDAQController <: VirtualDevice @add_device_fields TxDAQControllerParams - ref::Union{Array{Float32, 4}, Nothing} = nothing # TODO remove when done - cont::Union{Nothing, ControlSequence} = nothing # TODO remove when done + ref::Union{Array{Float32, 4}, Nothing} = nothing + cont::Union{Nothing, ControlSequence} = nothing + startFrame::Int64 = 1 + controlResults::OrderedDict{String, Union{typeof(1.0im*u"V/T"), Dict{Float64,typeof(1.0im*u"V/T")}}} = Dict{String, Union{typeof(1.0im*u"V/T"), Dict{Float64,typeof(1.0im*u"V/T")}}}() + lastDCResults::Union{Vector{@NamedTuple{V::Vector{Float64}, B::Vector{Float64}}},Nothing} = nothing + lastChannelIDs::Vector{String} = String[] +end + +function calibration(txCont::TxDAQController, channelID::AbstractString) + if haskey(txCont.controlResults, channelID) && txCont.controlResults[channelID] isa typeof(1.0u"V/T") + return txCont.controlResults[channelID] + else + return calibration(dependency(txCont, AbstractDAQ), channelID) + end +end + +function calibration(txCont::TxDAQController, channelID::AbstractString, frequency::Real) + if haskey(txCont.controlResults, channelID) && txCont.controlResults[channelID] isa Dict && haskey(txCont.controlResults[channelID],frequency) + return txCont.controlResults[channelID][frequency] + else + return calibration(dependency(txCont, AbstractDAQ), channelID, frequency) + end end function _init(tx::TxDAQController) # NOP end +function close(txCont::TxDAQController) + # NOP +end + neededDependencies(::TxDAQController) = [AbstractDAQ] optionalDependencies(::TxDAQController) = [SurveillanceUnit, Amplifier, TemperatureController] -function ControlSequence(txCont::TxDAQController, target::Sequence, daq::AbstractDAQ) - # Prepare Periodic Electrical Components - seqControlledChannel = getControlledChannel(txCont, target) - simpleChannel = createPeriodicElectricalComponentDict(seqControlledChannel, target, daq) - sinLUT, cosLUT = createLUTs(seqControlledChannel, target) +function checkIfControlPossible(txCont::TxDAQController, target::Sequence) + daq = dependency(txCont, AbstractDAQ) + + if any(.!isa.(getControlledChannels(target), ElectricalTxChannel)) + error("A field that requires control can only have ElectricalTxChannels.") + end + + if !all(unitIsTesla.(getControlledChannels(target))) + error("All values corresponding to a field that should be controlled need to be given in T, instead of V or A!") + end + + periodicChannels = [channel for channel in getControlledChannels(target) if typeof(channel) <: PeriodicElectricalChannel] + lutChannels = [channel for channel in getControlledChannels(target) if typeof(channel) <: AcyclicElectricalTxChannel] + + if !isempty(lutChannels) #&& !txCont.params.controlDC + error("A field that requires control can only have PeriodicElectricalChannels if controlDC is set to false for the TxDAQController!") + end + + outputsPeriodic = channelIdx(daq, id.(periodicChannels)) + if !allunique(outputsPeriodic) + error("Multiple periodic field channels are output on the same DAC output. This can not be controlled!") + end + + mapToFastDAC(i) = begin x = mod(i,6); if x==1 || x==2; return 2*(i÷6)+x end end + outputsLUT = channelIdx(daq, id.(lutChannels)) + if !all([x ∈ outputsPeriodic for x = mapToFastDAC.(outputsLUT)]) + error("If AcyclicElectricalTxChannels should be controlled they must map to the same output as a controlled PeriodicElectricalChannel!") + end + +end - currSeq = target +function ControlSequence(txCont::TxDAQController, target::Sequence) + daq = dependency(txCont, AbstractDAQ) - _name = "Control Sequence for target $(name(target))" - description = "" - _targetScanner = targetScanner(target) - _baseFrequency = baseFrequency(target) - general = GeneralSettings(;name=_name, description = description, targetScanner = _targetScanner, baseFrequency = _baseFrequency) - acq = AcquisitionSettings(;channels = RxChannel[], bandwidth = rxBandwidth(target)) + currSeq = prepareSequenceForControl(txCont, target) - _fields = MagneticField[] - for field in fields(target) - if control(field) - _id = id(field) - safeStart = safeStartInterval(field) - safeTrans = safeTransitionInterval(field) - safeEnd = safeEndInterval(field) - safeError = safeErrorInterval(field) - # Init purely periodic electrical component channel - periodicChannel = [deepcopy(channel) for channel in periodicElectricalTxChannels(field) if length(arbitraryElectricalComponents(channel)) == 0] - periodicComponents = [comp for channel in periodicChannel for comp in periodicElectricalComponents(channel)] - for channel in periodicChannel - otherComp = filter(!in(periodicElectricalComponents(channel)), periodicComponents) - for comp in periodicElectricalComponents(channel) - # TODO smarter init value, maybe min between 1/10 and target (<- target might not be V) - amplitude!(comp, simpleChannel[channel].limitPeak/10) - end - if txCont.params.correctCrossCoupling - for comp in otherComp - copy = deepcopy(comp) - amplitude!(copy, 0.0u"V") - push!(channel, copy) + measuredFrames = max(cld(txCont.params.minimumStepDuration, ustrip(u"s", dfCycle(currSeq))),1) + discardedFrames = cld(txCont.params.timeUntilStable, ustrip(u"s", dfCycle(currSeq))) + txCont.startFrame = discardedFrames + 1 + acqNumFrames(currSeq, discardedFrames+measuredFrames) + + applyForwardCalibration!(currSeq, txCont) # uses the forward calibration to convert the values for the field from T to V + + seqControlledChannels = getControlledChannels(currSeq) + + @debug "ControlSequence" seqControlledChannels + + # Dict(PeriodicElectricalChannel => TxChannelParams) + controlledChannelsDict = createControlledChannelsDict(seqControlledChannels, daq) # should this be changed to components instead of channels? + refIndices, refTFs = createReferenceIndexMapping(controlledChannelsDict, daq) + + controlSequenceType = decideControlSequenceType(target, txCont.params.controlDC) + @debug "ControlSequence: Decided on using ControlSequence of type $controlSequenceType" + + if controlSequenceType==CrossCouplingControlSequence + + # temporarily remove first (and only) component from each channel + comps = [popfirst!(channel.components) for channel in periodicElectricalTxChannels(fields(currSeq)[1])] + + # insert all components into each channel in the same order, all comps from the other channels are set to 0T + for (i, channel) in enumerate(periodicElectricalTxChannels(fields(currSeq)[1])) + for (j, comp) in enumerate(comps) + copy_ = deepcopy(comp) + if i!=j + amplitude!(copy_, 0.0u"V") end + push!(channel.components, copy_) end end - # Init arbitrary waveform - # TODO Implement AWG - contField = MagneticField(;id = _id, channels = periodicChannel, safeStartInterval = safeStart, safeTransitionInterval = safeTrans, - safeEndInterval = safeEnd, safeErrorInterval = safeError, control = true) - push!(_fields, contField) + + sinLUT, cosLUT = createLUTs(seqControlledChannels, currSeq) # TODO/JA: check if this makes a difference between target and currSeq, though it should not + + return CrossCouplingControlSequence(target, currSeq, controlledChannelsDict, refIndices, refTFs, sinLUT, cosLUT) + + elseif controlSequenceType == AWControlSequence # use the new controller + + # AW offset will get ignored -> should be zero + if any(x->any(mean.(scaledValues.(arbitraryElectricalComponents(x))).>10u"µT"), getControlledChannels(target)) + error("The DC-component of arbitrary waveform components cannot be handled during control! Please remove any DC-offset from your waveform and use the offset parameter of the corresponding channel!") + end + + rfftIndices = createRFFTindices(controlledChannelsDict, target, daq) + + cont = AWControlSequence(target, currSeq, controlledChannelsDict, refIndices, refTFs, rfftIndices, []) + + ## Apply last DC result + if !isnothing(txCont.lastDCResults) && (txCont.lastChannelIDs == id.(seqControlledChannels)) + cont.dcSearch = txCont.lastDCResults[end-1:end] + Ω = calcDesiredField(cont) + initTx = calcControlMatrix(cont) + last = cont.dcSearch[end] + previous = cont.dcSearch[end-1] + initTx[:,1] .= previous.V .- ((previous.B.-Ω[:,1]).*(last.V.-previous.V))./(last.B.-previous.B) + @info "Would have reused last DC Results" initTx[:,1] + #updateControlSequence!(cont, initTx) end + + return cont end +end - # Create Ref Indexing - mapping = Dict( b => a for (a,b) in enumerate(channelIdx(daq, daq.refChanIDs))) - controlOrderChannelIndices = [channelIdx(daq, ch.feedback.channelID) for ch in collect(Base.values(simpleChannel))] - refIndices = [mapping[x] for x in controlOrderChannelIndices] +function decideControlSequenceType(target::Sequence, controlDC::Bool=false) + + hasAWComponent = any(isa.([component for channel in getControlledChannels(target) for component in channel.components], ArbitraryElectricalComponent)) + moreThanOneComponent = any(x -> length(x.components) > 1, getControlledChannels(target)) + moreThanThreeChannels = length(getControlledChannels(target)) > 3 + moreThanOneField = length(getControlledFields(target)) > 1 + needsDecoupling_ = needsDecoupling(target) + @debug "decideControlSequenceType:" hasAWComponent moreThanOneComponent moreThanThreeChannels moreThanOneField needsDecoupling_ controlDC + + if needsDecoupling_ && !hasAWComponent && !moreThanOneField && !moreThanThreeChannels && !moreThanOneComponent && !controlDC + return CrossCouplingControlSequence + elseif needsDecoupling_ + throw(SequenceConfigurationError("The given sequence can not be controlled! To control a field with decoupling it cannot have an AW component ($hasAWComponent), more than one field ($moreThanOneField), more than three channels ($moreThanThreeChannels) nor more than one component per channel ($moreThanOneComponent). DC control ($controlDC) is also not possible")) + elseif !hasAWComponent && !moreThanOneField && !moreThanThreeChannels && !moreThanOneComponent && !controlDC + return CrossCouplingControlSequence + else + return AWControlSequence + end +end + +""" +This function creates the correct indices into a channel-wise rfft of the reference channels for each channel in the controlledChannelsDict +""" +function createRFFTindices(controlledChannelsDict::OrderedDict{PeriodicElectricalChannel, TxChannelParams}, seq::Sequence, daq::AbstractDAQ) + # Goal: create a Vector of index masks for each component + # N for every single-frequency component that has N periods in the sequence + # every Mth sample for all arbitraryWaveform components + numControlledChannels = length(keys(controlledChannelsDict)) + rfftSize = Int(div(rxNumSamplingPoints(seq),2)+1) + index_mask = falses(numControlledChannels, numComponentsMax(daq)+1, rfftSize) + + refChannelIdx = [channelIdx(daq, ch.feedbackChannelID) for ch in collect(Base.values(controlledChannelsDict))] + + + for (i, channel) in enumerate(keys(controlledChannelsDict)) + for (j, comp) in enumerate(components(channel)) + dfCyclesPerPeriod = Int(dfSamplesPerCycle(seq)/divider(comp)) + if isa(comp, PeriodicElectricalComponent) + index_mask[i,j,dfCyclesPerPeriod+1] = true + elseif isa(comp, ArbitraryElectricalComponent) + # the frequency samples have a spacing dfCyclesPerPeriod in the spectrum, only use a maximum number of 2^13+1 points, since the waveform has a buffer length of 2^14 (=> rfft 2^13+1) + N_harmonics = findlast(x->x>1e-8, abs.(rfft(values(comp)/0.5length(values(comp))))) + index_mask[i,j,dfCyclesPerPeriod+1:dfCyclesPerPeriod:min(rfftSize, N_harmonics*dfCyclesPerPeriod+1)] .= true + @debug "createRFFTindices: AWG component" divider(comp) sum(index_mask[i,j,:]) N_harmonics + end + end + index_mask[i,end,1] = ~any(index_mask[findall(x->x==refChannelIdx[i], refChannelIdx),end,1]) && channel.dcEnabled # use the first DC enabled channel going to each ref channel to control the DC value with + end + + # TODO/JA: test if this works + # Do a or on the mask of all different components + allComponentMask = .|(eachslice(index_mask, dims=2)...) + uniqueRefChIdx = unique(refChannelIdx) + combinedRefChannelMask = falses(length(uniqueRefChIdx), size(allComponentMask, 2)) + # if two controlled channels are on the same ref channel they should be combined, maybe there is a more elegant way to do this, but I cant figure it out + for i in 1:numControlledChannels + combinedRefChannelMask[findfirst(x->x==refChannelIdx[i], uniqueRefChIdx),:] .|= allComponentMask[i,:] + end + # if the total number of trues is smaller than the full array, there is an overlap of different components in the same FFT bin somewhere + if sum(combinedRefChannelMask) < sum(index_mask) + throw(SequenceConfigurationError("The controller can not control two different components, that have the same frequency on the same reference channel! This might also include arbitrary waveform components overlapping with normal components in frequency space or multiple TX channels using the same reference channel!")) + end + allOffsets = offset.(keys(controlledChannelsDict)) + for refCh in uniqueRefChIdx + if length(unique(allOffsets[findall(x->x==refCh, refChannelIdx)]))>1 + throw(SequenceConfigurationError("The offsets of all channels going to the same ref channel need to be identical!")) + end + end + + return index_mask +end - currSeq = Sequence(;general = general, acquisition = acq, fields = _fields) +allComponentMask(cont::AWControlSequence) = .|(eachslice(cont.rfftIndices, dims=2)...) - duration = div(txCont.params.minimumStepDuration, ustrip(u"s", dfCycle(currSeq)), RoundUp) - acqNumFrames(currSeq, max(duration, 1)) +function createLUTs(seqChannel::Vector{PeriodicElectricalChannel}, seq::Sequence) + N = rxNumSamplingPoints(seq) + D = length(seqChannel) - return ControlSequence(target, currSeq, simpleChannel, sinLUT, cosLUT, refIndices) + dfCyclesPerPeriod = Int[dfSamplesPerCycle(seq)/divider(components(chan)[i]) for (i,chan) in enumerate(seqChannel)] + + sinLUT = zeros(D,N) + cosLUT = zeros(D,N) + for d=1:D + for n=1:N + sinLUT[d,n] = sin(2 * pi * (n-1) * dfCyclesPerPeriod[d] / N) + cosLUT[d,n] = cos(2 * pi * (n-1) * dfCyclesPerPeriod[d] / N) + end + end + return sinLUT, cosLUT end -acyclicElectricalTxChannels(cont::ControlSequence) = acyclicElectricalTxChannels(cont.targetSequence) -periodicElectricalTxChannels(cont::ControlSequence) = periodicElectricalTxChannels(cont.targetSequence) -acqNumFrames(cont::ControlSequence) = acqNumFrames(cont.targetSequence) -acqNumFrameAverages(cont::ControlSequence) = acqNumFrameAverages(cont.targetSequence) -acqNumFrames(cont::ControlSequence, x) = acqNumFrames(cont.targetSequence, x) -acqNumFrameAverages(cont::ControlSequence, x) = acqNumFrameAverages(cont.targetSequence, x) +""" +Returns a vector of indices into the the DAQ data after the channels have already been selected by the daq.refChanIDs to select the correct reference channels in the order of the controlledChannels +""" +function createReferenceIndexMapping(controlledChannelsDict::OrderedDict{PeriodicElectricalChannel, TxChannelParams}, daq::AbstractDAQ) + # Dict(RedPitaya ChannelIndex => Index in daq.refChanIDs) + mapping = Dict( b => a for (a,b) in enumerate(channelIdx(daq, daq.refChanIDs))) + # RedPitaya ChannelIndex der Feedback-Kanäle in Reihenfolge der Kanäle im controlledChannelsDict + controlOrderChannelIndices = [channelIdx(daq, ch.feedbackChannelID) for ch in collect(Base.values(controlledChannelsDict))] + # Index in daq.refChanIDs in der Reihenfolge der Kanäle im controlledChannelsDict + refIndices = [mapping[x] for x in controlOrderChannelIndices] + # Feedback TransferFunction for the controlled Channels + refTFs = [feedbackTransferFunction(daq, ch) for ch in collect(Base.values(controlledChannelsDict))] + return refIndices, refTFs +end + +function prepareSequenceForControl(txCont::TxDAQController, seq::Sequence) + # Ausgangslage: Eine Sequenz enthält eine Reihe an Feldern, von denen ggf nicht alle geregelt werden sollen + # Jedes dieser Felder enthält eine Reihe an Channels, die nicht alle geregelt werden können (nur periodicElectricalChannel) + + # Ziel: Eine Sequenz, die, wenn Sie abgespielt wird, alle Informationen beinhaltet, die benötigt werden, um die Regelung durchzuführen, Sendefelder in V + checkIfControlPossible(txCont, seq) + + _name = "Control Sequence for target $(name(seq))" + description = "" + _targetScanner = targetScanner(seq) + _baseFrequency = baseFrequency(seq) + general = GeneralSettings(;name=_name, description = description, targetScanner = _targetScanner, baseFrequency = _baseFrequency) + acq = AcquisitionSettings(;channels = RxChannel[], bandwidth = rxBandwidth(seq)) # uses the default values of 1 for numPeriodsPerFrame, numFrames, numAverages, numFrameAverages + + return Sequence(;general = general, acquisition = acq, fields = [deepcopy(f) for f in fields(seq) if control(f)]) +end -function createPeriodicElectricalComponentDict(seqControlledChannel::Vector{PeriodicElectricalChannel}, seq::Sequence, daq::AbstractDAQ) +function createControlledChannelsDict(seqControlledChannels::Vector{PeriodicElectricalChannel}, daq::AbstractDAQ) missingControlDef = [] dict = OrderedDict{PeriodicElectricalChannel, TxChannelParams}() - for seqChannel in seqControlledChannel + for seqChannel in seqControlledChannels name = id(seqChannel) daqChannel = get(daq.params.channels, name, nothing) - if isnothing(daqChannel) || isnothing(daqChannel.feedback) || !in(daqChannel.feedback.channelID, daq.refChanIDs) + if isnothing(daqChannel) || isnothing(daqChannel.feedbackChannelID) || !in(daqChannel.feedbackChannelID, daq.refChanIDs) + @debug "Found missing control def: " name isnothing(daqChannel) isnothing(daqChannel.feedbackChannelID) !in(daqChannel.feedbackChannelID, daq.refChanIDs) push!(missingControlDef, name) else dict[seqChannel] = daqChannel end end - + if length(missingControlDef) > 0 message = "The sequence requires control for the following channel " * join(string.(missingControlDef), ", ", " and ") * ", but either the channel was not defined or had no defined feedback channel." throw(IllegalStateException(message)) end - # Check that we only control three channels, as our RedPitayaDAQs only have 3 signal components atm - if length(dict) > 3 - throw(IllegalStateException("Sequence requires controlling of more than four channels, which is currently not implemented.")) - end - - # Check that channels only have one component - if any(x -> length(x.components) > 1, seqControlledChannel) - throw(IllegalStateException("Sequence has channel with more than one component. Such a channel cannot be controlled by this controller")) - end return dict end -function controlTx(txCont::TxDAQController, seq::Sequence, ::Nothing = nothing) - daq = dependency(txCont, AbstractDAQ) - setupRx(daq, seq) - control = ControlSequence(txCont, seq, daq) - return controlTx(txCont, seq, control) -end +############################################################################### +############## Top-Level functions for interacting with the controller +############################################################################### +function controlTx(txCont::TxDAQController, seq::Sequence) + if needsControlOrDecoupling(seq) + daq = dependency(txCont, AbstractDAQ) + setupRx(daq, seq) + control = ControlSequence(txCont, seq) # depending on the controlled channels and settings this will select the appropiate type of ControlSequence + return controlTx(txCont, control) + else + @warn "The sequence you selected does not need control, even though the protocol wanted to control!" + return seq + end +end -function controlTx(txCont::TxDAQController, seq::Sequence, control::ControlSequence) +function controlTx(txCont::TxDAQController, control::ControlSequence) # Prepare and check channel under control daq = dependency(txCont, AbstractDAQ) Ω = calcDesiredField(control) txCont.cont = control + if !checkVoltLimits(calcControlMatrix(control), control) + error("Initial guess for the controller already exceeds the possibilities of the DAQ!") + end # Start Tx su = nothing @@ -186,7 +401,7 @@ function controlTx(txCont::TxDAQController, seq::Sequence, control::ControlSeque end if !isempty(amps) # Only enable amps that amplify a channel of the current sequence - txChannelIds = id.(vcat(acyclicElectricalTxChannels(seq), periodicElectricalTxChannels(seq))) + txChannelIds = id.(vcat(acyclicElectricalTxChannels(control.targetSequence), periodicElectricalTxChannels(control.targetSequence))) amps = filter(amp -> in(channelId(amp), txChannelIds), amps) @sync for amp in amps @async turnOn(amp) @@ -195,16 +410,30 @@ function controlTx(txCont::TxDAQController, seq::Sequence, control::ControlSeque # Hacky solution controlPhaseDone = false + controlPhaseError = nothing i = 1 - len = length(keys(control.simpleChannel)) + try + while !controlPhaseDone && i <= txCont.params.maxControlSteps @info "CONTROL STEP $i" # Prepare control measurement setup(daq, control.currSequence) + + if haskey(ENV, "JULIA_DEBUG") && contains(ENV["JULIA_DEBUG"],"checkEachControlStep") + menu = REPL.TerminalMenus.RadioMenu(["Continue", "Abort"], pagesize=2) + choice = REPL.TerminalMenus.request("Please confirm the current values for control:", menu) + if choice == 1 + println("Continuing...") + else + println("Control cancelled") + error("Control cancelled!") + end + end + channel = Channel{channelType(daq)}(32) - buffer = AsyncBuffer(FrameSplitterBuffer(daq, StorageBuffer[DriveFieldBuffer(1, zeros(ComplexF64,len, len, 1, acqNumFrames(control.currSequence)), control)]), daq) - @info "Control measurement started" + buffer = AsyncBuffer(FrameSplitterBuffer(daq, StorageBuffer[DriveFieldBuffer(1, zeros(ComplexF64, controlMatrixShape(control)..., 1, acqNumFrames(control.currSequence)), control)]), daq) + @debug "Control measurement started" producer = @async begin @debug "Starting control producer" endSample = asyncProducer(channel, daq, control.currSequence, isControlStep=true) @@ -221,14 +450,18 @@ function controlTx(txCont::TxDAQController, seq::Sequence, control::ControlSeque end end wait(consumer) - @info "Control measurement finished" - - @info "Evaluating control step" - uRef = read(sink(buffer, DriveFieldBuffer))[:, :, 1, end] - if !isnothing(uRef) - controlPhaseDone = controlStep!(control, txCont, uRef, Ω) == UNCHANGED + @debug "Control measurement finished" + + @debug "Evaluating control step" + tmp = read(sink(buffer, DriveFieldBuffer)) + @debug "Size of calc fields from ref" size(tmp) + + Γ = mean(tmp[:, :, 1, txCont.startFrame:end],dims=3)[:,:,1] # calcFieldsFromRef happened here already + if !isnothing(Γ) + controlPhaseDone = controlStep!(control, txCont, Γ, Ω) == UNCHANGED if controlPhaseDone @info "Could control" + updateCachedCalibration(txCont, control) else @info "Could not control" end @@ -239,6 +472,8 @@ function controlTx(txCont::TxDAQController, seq::Sequence, control::ControlSeque end catch ex @error "Exception during control loop" exception=(ex, catch_backtrace()) + resetCachedCalibration(txCont) + controlPhaseError = ex finally try stopTx(daq) @@ -273,47 +508,135 @@ function controlTx(txCont::TxDAQController, seq::Sequence, control::ControlSeque end if !controlPhaseDone - error("TxDAQController $(deviceID(txCont)) could not control.") + if isnothing(controlPhaseError) + error("TxDAQController $(deviceID(txCont)) could not reach a stable field after $(txCont.params.maxControlSteps) steps.") + else + error("TxDAQController $(deviceID(txCont)) failed control with the following message:\n$(sprint(showerror, controlPhaseError))") + end end return control end +""" +Returns a Sequence that is merged from the control result and all uncontrolled field of the given ControlSequence +""" +function getControlResult(cont::ControlSequence)::Sequence + + # Use the magnetic field that are controlled from currSeq and all uncontrolled fields and general settings from target -function setup(daq::RedPitayaDAQ, sequence::ControlSequence) - stopTx(daq) - setupRx(daq, sequence.targetSequence) - setupTx(daq, sequence.currSequence) - prepareTx(daq, sequence.currSequence) - setSequenceParams(daq, sequence.targetSequence) -end + _name = "Control Result for target $(name(cont.targetSequence))" + general = GeneralSettings(;name=_name, description = description(cont.targetSequence), targetScanner = targetScanner(cont.targetSequence), baseFrequency = baseFrequency(cont.targetSequence)) + acq = cont.targetSequence.acquisition -getControlledChannel(::TxDAQController, seq::Sequence) = [channel for field in seq.fields if field.control for channel in field.channels if typeof(channel) <: PeriodicElectricalChannel && length(arbitraryElectricalComponents(channel)) == 0] -getUncontrolledChannel(::TxDAQController, seq::Sequence) = [channel for field in seq.fields if !field.control for channel in field.channels if typeof(channel) <: PeriodicElectricalChannel] + _fields = MagneticField[] + for field in fields(cont.currSequence) + _id = id(field) + safeStart = safeStartInterval(field) + safeTrans = safeTransitionInterval(field) + safeEnd = safeEndInterval(field) + safeError = safeErrorInterval(field) + contField = MagneticField(;id = _id, channels = deepcopy(channels(field)), safeStartInterval = safeStart, safeTransitionInterval = safeTrans, + safeEndInterval = safeEnd, safeErrorInterval = safeError, decouple = false, control = false) + push!(_fields, contField) + end + for field in fields(cont.targetSequence) + if !control(field) + push!(_fields, field) + # TODO/JA: if there are LUT channels sharing a channel with the controlled fields, we should be able to use the DC calibration that has been found + # and insert it into the corresponding LUT channels as a calibration value + end + end -function createLUTs(seqChannel::Vector{PeriodicElectricalChannel}, seq::Sequence) - N = rxNumSamplingPoints(seq) - D = length(seqChannel) - cycle = ustrip(dfCycle(seq)) - base = ustrip(dfBaseFrequency(seq)) - dfFreq = [base/x.components[1].divider for x in seqChannel] + return Sequence(;general = general, acquisition = acq, fields = _fields) +end + +setup(daq::AbstractDAQ, sequence::ControlSequence) = setup(daq, getControlResult(sequence)) + +function updateCachedCalibration(txCont::TxDAQController, cont::ControlSequence) + finalCalibration = diag(calcControlMatrix(cont) ./ calcDesiredField(cont)) + channelIds = id.(getControlledChannels(cont)) + dividers = divider.(MPIMeasurements.getPrimaryComponents(cont)) + frequencies = ustrip(u"Hz", txBaseFrequency(cont.currSequence)) ./ dividers + for i in axes(finalCalibration, 1) + if !isnan(finalCalibration[i]) + if !haskey(txCont.controlResults, channelIds[i]) + txCont.controlResults[channelIds[i]] = Dict{Float64,typeof(1.0im*u"V/T")}() + end + txCont.controlResults[channelIds[i]][frequencies[i]] = finalCalibration[i]*u"V/T" + @debug "Cached control result: $(round(finalCalibration[i], digits=2)*u"V/T") for channel $(channelIds[i]) at $(round(frequencies[i],digits=3)) Hz" + end + end - sinLUT = zeros(D,N) - cosLUT = zeros(D,N) - for d=1:D - Y = round(Int64, cycle*dfFreq[d] ) - for n=1:N - sinLUT[d,n] = sin(2 * pi * (n-1) * Y / N) - cosLUT[d,n] = cos(2 * pi * (n-1) * Y / N) + nothing +end + +function updateCachedCalibration(txCont::TxDAQController, cont::AWControlSequence) + finalCalibration = calcControlMatrix(cont) ./ calcDesiredField(cont) + + calibrationResults = findall(x->!isnan(x), finalCalibration) + channelIDs = id.(getControlledChannels(cont)) + freqAxis = rfftfreq(rxNumSamplingPoints(cont.currSequence),ustrip(u"Hz",2*rxBandwidth(cont.currSequence))) + + for res in calibrationResults + chId = channelIDs[res[1]] + f = freqAxis[res[2]] + if !haskey(txCont.controlResults, chId) + txCont.controlResults[chId] = Dict{Float64,typeof(1.0im*u"V/T")}() end + txCont.controlResults[chId][f] = finalCalibration[res]*u"V/T" + @debug "Cached calibration result:" chId f finalCalibration[res] end - return sinLUT, cosLUT + + if length(cont.dcSearch) >= 2 + txCont.lastDCResults = cont.dcSearch[end-1:end] + else + txCont.lastDCResults = nothing + end + txCont.lastChannelIDs = channelIDs + + @debug "Cached DC result" txCont.lastDCResults +end + +function resetCachedCalibration(txCont::TxDAQController) + txCont.controlResults = Dict{String, Union{typeof(1.0im*u"V/T"), Dict{Float64,typeof(1.0im*u"V/T")}}}() + txCont.lastDCResults = nothing + txCont.lastChannelIDs = String[] + @debug "Reset cached calibration" + nothing end + +# TODO/JA: check if changes here needs changes somewhere else +getControlledFields(seq::Sequence) = [field for field in seq.fields if field.control] +getControlledChannels(seq::Sequence) = [channel for field in seq.fields if field.control for channel in field.channels] +# The elements of collect(getControlledChannels(cont)) are always identical (===) to getControlledChannels(cont.currSequence) +getControlledChannels(cont::ControlSequence) = keys(cont.controlledChannelsDict) # maybe add collect here as well? for most uses this not needed +getControlledDAQChannels(cont::ControlSequence) = collect(Base.values(cont.controlledChannelsDict)) +getPrimaryComponents(cont::CrossCouplingControlSequence) = [components(channel)[i] for (i,channel) in enumerate(getControlledChannels(cont))] + +acyclicElectricalTxChannels(cont::ControlSequence) = acyclicElectricalTxChannels(cont.targetSequence) +periodicElectricalTxChannels(cont::ControlSequence) = periodicElectricalTxChannels(cont.targetSequence) +acqNumFrames(cont::ControlSequence) = acqNumFrames(cont.targetSequence) +acqNumFrameAverages(cont::ControlSequence) = acqNumFrameAverages(cont.targetSequence) +acqNumFrames(cont::ControlSequence, x) = acqNumFrames(cont.targetSequence, x) +acqNumFrameAverages(cont::ControlSequence, x) = acqNumFrameAverages(cont.targetSequence, x) +numControlledChannels(cont::ControlSequence) = length(getControlledChannels(cont)) + +controlMatrixShape(cont::AWControlSequence) = (numControlledChannels(cont), size(cont.rfftIndices,3)) +controlMatrixShape(cont::CrossCouplingControlSequence) = (numControlledChannels(cont), numControlledChannels(cont)) + + +################################################################################# +########## Functions for the control steps +################################################################################# + + + controlStep!(cont::ControlSequence, txCont::TxDAQController, uRef) = controlStep!(cont, txCont, uRef, calcDesiredField(cont)) -controlStep!(cont::ControlSequence, txCont::TxDAQController, uRef, Ω::Matrix) = controlStep!(cont, txCont, calcFieldFromRef(cont, uRef), Ω) -function controlStep!(cont::ControlSequence, txCont::TxDAQController, Γ::Matrix, Ω::Matrix) - if checkFieldDeviation(Γ, Ω, txCont) +controlStep!(cont::ControlSequence, txCont::TxDAQController, uRef, Ω::Matrix{<:Complex}) = controlStep!(cont, txCont, calcFieldsFromRef(cont, uRef), Ω) +function controlStep!(cont::ControlSequence, txCont::TxDAQController, Γ::Matrix{<:Complex}, Ω::Matrix{<:Complex}) + if fieldAccuracyReached(cont, txCont, Γ, Ω) return UNCHANGED elseif updateControl!(cont, txCont, Γ, Ω) return UPDATED @@ -322,30 +645,139 @@ function controlStep!(cont::ControlSequence, txCont::TxDAQController, Γ::Matrix end end -calcFieldFromRef(cont::ControlSequence, uRef; frame::Int64 = 1, period::Int64 = 1) = calcFieldFromRef(cont, uRef, UnsortedRef(), frame = frame, period = period) -function calcFieldFromRef(cont::ControlSequence, uRef::Array{Float32, 4}, ::UnsortedRef; frame::Int64 = 1, period::Int64 = 1) - return calcFieldFromRef(cont, uRef[:, :, :, frame], UnsortedRef(), period = period) +#fieldAccuracyReached(cont::ControlSequence, txCont::TxDAQController, uRef) = fieldAccuracyReached(cont, txCont, calcFieldFromRef(cont, uRef)) +fieldAccuracyReached(cont::ControlSequence, txCont::TxDAQController, Γ::Matrix{<:Complex}) = fieldAccuracyReached(cont, txCont, Γ, calcDesiredField(cont)) +function fieldAccuracyReached(cont::CrossCouplingControlSequence, txCont::TxDAQController, Γ::Matrix{<:Complex}, Ω::Matrix{<:Complex}) + + abs_deviation = abs.(Ω) .- abs.(Γ) + rel_deviation = abs_deviation ./ abs.(Ω) + rel_deviation[abs.(Ω).<1e-15] .= 0 # relative deviation does not make sense for a zero goal + phase_deviation = angle.(Ω).-angle.(Γ) + phase_deviation[abs.(Ω).<1e-15] .= 0 # phase deviation does not make sense for a zero goal + + if !needsDecoupling(cont.targetSequence) + abs_deviation = diag(abs_deviation) + rel_deviation = diag(rel_deviation) + phase_deviation = diag(phase_deviation) + # elseif isa(cont, AWControlSequence) + # abs_deviation = abs_deviation[allComponentMask(cont)]' # TODO/JA: keep the distinction between the channels (maybe as Vector{Vector{}}), instead of putting everything into a long vector with unknown order + # rel_deviation = rel_deviation[allComponentMask(cont)]' + # phase_deviation = phase_deviation[allComponentMask(cont)]' + # Γt = checkVoltLimits(Γ,cont,return_time_signal=true)' + # Ωt = checkVoltLimits(Ω,cont,return_time_signal=true)' + + # diff = (Ωt .- Γt) + # @debug "fieldAccuracyReached" diff=lineplot(1:rxNumSamplingPoints(cont.currSequence),diff, canvas=DotCanvas, border=:ascii) + # @info "fieldAccuracyReached2" max_diff = maximum(abs.(diff)) + end + @debug "Check field deviation [T]" Ω Γ + @debug "Ω - Γ = " abs_deviation rel_deviation phase_deviation + @info "Observed field deviation:\nabs:\t$(abs_deviation*1000) mT\nrel:\t$(rel_deviation*100) %\nphi:\t$(phase_deviation/pi*180)°\n allowed: $(txCont.params.absoluteAmplitudeAccuracy|>u"mT"), $(txCont.params.relativeAmplitudeAccuracy*100) %, $(uconvert(u"°",txCont.params.phaseAccuracy))" + phase_ok = abs.(phase_deviation) .< ustrip(u"rad", txCont.params.phaseAccuracy) + amplitude_ok = (abs.(abs_deviation) .< ustrip(u"T", txCont.params.absoluteAmplitudeAccuracy)) .| (abs.(rel_deviation) .< txCont.params.relativeAmplitudeAccuracy) + @debug "Field deviation:" amplitude_ok phase_ok + return all(phase_ok) && all(amplitude_ok) end -function calcFieldFromRef(cont::ControlSequence, uRef::Array{Float32, 3}, ::UnsortedRef; period::Int64 = 1) - return calcFieldFromRef(cont, view(uRef[:, cont.refIndices, :], :, :, period), SortedRef()) + +function fieldAccuracyReached(cont::AWControlSequence, txCont::TxDAQController, Γ::Matrix{<:Complex}, Ω::Matrix{<:Complex}) + + Γt = transpose(checkVoltLimits(Γ,cont,return_time_signal=true)) + Ωt = transpose(checkVoltLimits(Ω,cont,return_time_signal=true)) + @debug "fieldAccuracyReached" transpose(Γ[allComponentMask(cont)]) abs.(Γ[allComponentMask(cont)])' angle.(Γ[allComponentMask(cont)])'# abs.(Ω[allComponentMask(cont)])' angle.(Ω[allComponentMask(cont)])' + diff = (Ωt .- Γt) + zero_mean_diff = diff .- mean(diff, dims=1) + @debug "fieldAccuracyReached" diff=lineplot(1:rxNumSamplingPoints(cont.currSequence),diff*1000, canvas=DotCanvas, border=:ascii, ylabel="mT", name=dependency(txCont, AbstractDAQ).refChanIDs[cont.refIndices]) + @info "Observed field deviation (time-domain):\nmax_diff:\t$(maximum(abs.(diff))*1000) mT\nmax_diff (w/o DC): \t$(maximum(abs.(zero_mean_diff))*1000)" + amplitude_ok = abs.(diff).< ustrip(u"T", txCont.params.absoluteAmplitudeAccuracy) + return all(amplitude_ok) end -function calcFieldsFromRef(cont::ControlSequence, uRef::Array{Float32, 4}) - len = length(keys(cont.simpleChannel)) + +#updateControl!(cont::ControlSequence, txCont::TxDAQController, uRef) = updateControl!(cont, txCont, calcFieldFromRef(cont, uRef), calcDesiredField(cont)) +function updateControl!(cont::ControlSequence, txCont::TxDAQController, Γ::Matrix{<:Complex}, Ω::Matrix{<:Complex}) + @debug "Updating control values" + oldTx = calcControlMatrix(cont) + + if !validateAgainstForwardCalibrationAndSafetyLimit(oldTx, Γ, cont, txCont) + error("Last control step produced unexpected results! Either your forward calibration is inaccurate or the system is not in the expected state (e.g. amp not on)!" ) + end + newTx = updateControlMatrix(cont, txCont, Γ, Ω, oldTx) + + if validateAgainstForwardCalibrationAndSafetyLimit(newTx, Ω, cont, txCont) && checkVoltLimits(newTx, cont) + updateControlSequence!(cont, newTx) + return true + else + error("The new tx values are not allowed! Either your forward calibration is inaccurate or the system can not produce the requested field strength!") + return false + end +end + +# Γ: Matrix from Ref +# Ω: Desired Matrix +# oldTx: Last Set Matrix +function updateControlMatrix(cont::CrossCouplingControlSequence, txCont::TxDAQController, Γ::Matrix{<:Complex}, Ω::Matrix{<:Complex}, oldTx::Matrix{<:Complex}) + if needsDecoupling(cont.targetSequence) + β = Γ*inv(oldTx) + else + β = diagm(diag(Γ))*inv(diagm(diag(oldTx))) + end + newTx = inv(β)*Ω + @debug "Last TX matrix [V]:" oldTx + @debug "Ref matrix [T]:" Γ + @debug "Desired matrix [T]:" Ω + @debug "New TX matrix [V]:" newTx + return newTx +end + +function updateControlMatrix(cont::AWControlSequence, txCont::TxDAQController, Γ::Matrix{<:Complex}, Ω::Matrix{<:Complex}, oldTx::Matrix{<:Complex}) + # For now we completely ignore coupling and hope that it can find good values anyways + # The problem is, that to achieve 0 we will always output zero, but we would need a much more sophisticated method to solve this + newTx = oldTx./Γ.*Ω + + # handle DC separately: + if txCont.params.controlDC + push!(cont.dcSearch, (V=oldTx[:,1], B=Γ[:,1])) + if length(cont.dcSearch)==1 + my_sign(x) = if x<0; -1 else 1 end + testOffset = real.(Ω[:,1])*u"T" .- 2u"mT"*my_sign.(real.(Ω[:,1])) + newTx[:,1] = ustrip.(u"V", testOffset.*[abs(calibration(dependency(txCont, AbstractDAQ), id(channel))(0)) for channel in getControlledChannels(cont)]) + else + @debug "History of dcSearch" cont.dcSearch + last = cont.dcSearch[end] + previous = cont.dcSearch[end-1] + newTx[:,1] .= previous.V .- ((previous.B.-Ω[:,1]).*(last.V.-previous.V))./(last.B.-previous.B) + end + end + + #@debug "Last TX matrix [V]:" oldTx=lineplot(1:rxNumSamplingPoints(cont.currSequence),checkVoltLimits(oldTx,cont,return_time_signal=true)') + #@debug "Ref matrix [T]:" Γ=lineplot(1:rxNumSamplingPoints(cont.currSequence),checkVoltLimits(Γ,cont,return_time_signal=true)') + #@debug "Desired matrix [V]:" Ω=lineplot(1:rxNumSamplingPoints(cont.currSequence),checkVoltLimits(Ω,cont,return_time_signal=true)') + #@debug "New TX matrix [T]:" newTx=lineplot(1:rxNumSamplingPoints(cont.currSequence),checkVoltLimits(newTx,cont,return_time_signal=true)') + + return newTx +end + +################################################################################# +########## Functions for calculating the field matrix in T from the reference channels +################################################################################# + +function calcFieldsFromRef(cont::CrossCouplingControlSequence, uRef::Array{Float32, 4}) + len = numControlledChannels(cont) N = rxNumSamplingPoints(cont.currSequence) - dividers = Int64[divider(components(channel)[1]) for channel in keys(cont.simpleChannel)] + dividers = divider.(getPrimaryComponents(cont)) + frequencies = ustrip(u"Hz", txBaseFrequency(cont.currSequence)) ./ dividers Γ = zeros(ComplexF64, len, len, size(uRef, 3), size(uRef, 4)) sorted = uRef[:, cont.refIndices, :, :] for i = 1:size(Γ, 4) for j = 1:size(Γ, 3) - calcFieldFromRef!(view(Γ, :, :, j, i), cont, view(sorted, :, :, j, i), SortedRef()) + _calcFieldFromRef!(view(Γ, :, :, j, i), cont, view(sorted, :, :, j, i), SortedRef()) end end for d =1:len - c = ustrip(u"T/V", collect(Base.values(cont.simpleChannel))[d].feedback.calibration) + c = ustrip(u"T/V", 1 ./cont.refTFs[d](frequencies[d])) for e=1:len - correction = c * -im * dividers[e]/dividers[d] * 2/N + correction = c * dividers[e]/dividers[d] * 2/N for j = 1:size(Γ, 3) for i = 1:size(Γ, 4) Γ[d, e, j, i] = correction * Γ[d, e, j, i] @@ -356,23 +788,18 @@ function calcFieldsFromRef(cont::ControlSequence, uRef::Array{Float32, 4}) return Γ end -function calcFieldFromRef(cont::ControlSequence, uRef, ::SortedRef) - len = length(keys(cont.simpleChannel)) +function calcFieldsFromRef(cont::AWControlSequence, uRef::Array{Float32,4}) N = rxNumSamplingPoints(cont.currSequence) - dividers = Int64[divider(components(channel)[1]) for channel in keys(cont.simpleChannel)] - Γ = zeros(ComplexF64, len, len) - calcFieldFromRef!(Γ, cont, uRef, SortedRef()) - for d =1:len - c = ustrip(u"T/V", collect(Base.values(cont.simpleChannel))[d].feedback.calibration) - for e=1:len - correction = c * -im * dividers[e]/dividers[d] * 2/N - Γ[d,e] = correction * Γ[d,e] - end - end - return Γ + # do rfft channel wise and correct with the transfer function, return as (num control channels x len rfft x periods x frames) Matrix, the selection of [:,:,1,1] is done in controlTx + spectrum = rfft(uRef, 1)/0.5N + spectrum[1,:,:,:] ./= 2 + sortedSpectrum = permutedims(spectrum[:, cont.refIndices, :, :], (2,1,3,4)) + frequencies = ustrip.(u"Hz",rfftfreq(N, rxSamplingRate(cont.currSequence))) + fbTF = reduce(vcat, transpose([ustrip.(u"V/T", tf(frequencies)) for tf in cont.refTFs])) + return sortedSpectrum./fbTF end -function calcFieldFromRef!(Γ::AbstractArray{ComplexF64, 2}, cont::ControlSequence, uRef, ::SortedRef) +function _calcFieldFromRef!(Γ::AbstractArray{ComplexF64, 2}, cont::CrossCouplingControlSequence, uRef, ::SortedRef) len = size(Γ, 1) for d=1:len for e=1:len @@ -384,141 +811,181 @@ function calcFieldFromRef!(Γ::AbstractArray{ComplexF64, 2}, cont::ControlSequen return Γ end -function calcDesiredField(cont::ControlSequence) - seqChannel = keys(cont.simpleChannel) - temp = [ustrip(amplitude(components(ch)[1])) * exp(im*ustrip(phase(components(ch)[1]))) for ch in seqChannel] - return convert(Matrix{ComplexF64}, diagm(temp)) -end -checkFieldDeviation(uRef, cont::ControlSequence, txCont::TxDAQController) = checkFieldDeviation(calcFieldFromRef(cont, uRef), cont, txCont) -checkFieldDeviation(Γ::Matrix, cont::ControlSequence, txCont::TxDAQController) = checkFieldDeviation(Γ, calcDesiredField(cont), txCont) -function checkFieldDeviation(Γ::Matrix, Ω::Matrix, txCont::TxDAQController) - if txCont.params.correctCrossCoupling - diff = Ω - Γ - else - diff = diagm(diag(Ω)) - diagm(diag(Γ)) +################################################################################# +########## Functions for creating and applying the matrix representation +########## of the controller to the sequence(s) +################################################################################# + + +# Convert Target Sequence to Matrix in T + +function calcDesiredField(cont::CrossCouplingControlSequence) + desiredField = zeros(ComplexF64, controlMatrixShape(cont)) + for (i, channel) in enumerate(getControlledChannels(cont.targetSequence)) + comp = components(channel)[1] + desiredField[i, i] = ustrip(u"T", amplitude(comp)) * exp(im*ustrip(u"rad", phase(comp))) end - deviation = maximum(abs.(diff)) / maximum(abs.(Ω)) # TODO maximum(abs(diff)./abs(omega)) - @debug "Check field deviation" Ω Γ - @debug "Ω - Γ = " diff - @info "deviation = $(deviation) allowed= $(txCont.params.amplitudeAccuracy)" - return deviation < txCont.params.amplitudeAccuracy + return desiredField end -updateControl!(cont::ControlSequence, txCont::TxDAQController, uRef) = updateControl!(cont, txCont, calcFieldFromRef(cont, uRef), calcDesiredField(cont)) -function updateControl!(cont::ControlSequence, txCont::TxDAQController, Γ::Matrix, Ω::Matrix) - @debug "Updating control values" - κ = calcControlMatrix(cont) - newTx = updateControlMatrix(Γ, Ω, κ, correct_coupling = txCont.params.correctCrossCoupling) - if checkFieldToVolt(κ, Γ, cont, txCont) && checkVoltLimits(newTx, cont, txCont) - updateControlSequence!(cont, newTx) - return true - else - @warn "New control values are not allowed" - return false +function calcDesiredField(cont::AWControlSequence) + # generate desired spectrum per control channel that can be compared to the result of calcFieldsFromRef in the controlStep later + # size: (rfft of sequence x num control channels) + # the separation of individual components is done by the index masks + + desiredField = zeros(ComplexF64, controlMatrixShape(cont)) + + for (i, channel) in enumerate(getControlledChannels(cont.targetSequence)) + + if cont.rfftIndices[i,end,1] + desiredField[i,1] = ustrip(u"T", offset(channel)) + end + + for (j, comp) in enumerate(components(channel)) + if isa(comp, PeriodicElectricalComponent) + if ustrip(u"T",amplitude(comp)) == 0 + @warn "You tried to control a field to 0 T, this will just output 0 V on that channel, since this controller can not correct cross coupling" + end + desiredField[i, cont.rfftIndices[i,j,:]] .= ustrip(u"T",amplitude(comp)) * exp(im*ustrip(u"rad",phase(comp)-pi/2)) # The phase given in the component is for a sine, but the FFT-phase uses a cosine + elseif isa(comp, ArbitraryElectricalComponent) + desiredField[i, cont.rfftIndices[i,j,:]] .= rfft(ustrip.(u"T",scaledValues(comp)))[2:sum(cont.rfftIndices[i,j,:])+1]/(0.5*2^14) # the buffer length should always be 2^14 currently + end + end end + + return desiredField end -# Γ: Matrix from Ref -# Ω: Desired Matrix -# κ: Last Set Matrix -function updateControlMatrix(Γ::Matrix, Ω::Matrix, κ::Matrix; correct_coupling::Bool = false) - if correct_coupling - β = Γ*inv(κ) - else - β = diagm(diag(Γ))*inv(diagm(diag(κ))) + +# Convert Last Tx (currSequence) to Matrix in V +function calcControlMatrix(cont::CrossCouplingControlSequence) + oldTx = zeros(ComplexF64, controlMatrixShape(cont)) + for (i, channel) in enumerate(getControlledChannels(cont)) + for (j, comp) in enumerate(periodicElectricalComponents(channel)) + oldTx[i, j] = ustrip(u"V", amplitude(comp)) * exp(im*ustrip(u"rad", phase(comp))) + end end - - newTx = inv(β)*Ω - @debug "Last matrix:" κ - @debug "Ref matrix" Γ - @debug "Desired matrix" Ω - @debug "New matrix" newTx - return newTx + return oldTx end -function calcControlMatrix(cont::ControlSequence) - # TxDAQController only works on one field atm (possible future TODO: update to multiple fields, matrix per field) - field = fields(cont.currSequence)[1] - len = length(keys(cont.simpleChannel)) - κ = zeros(ComplexF64, len, len) - # In each channel the first component is the channels "own" component, the following are the ordered correction components of the other channel - # -> For Channel 2 its components in the matrix row should be c2 c1 c3 for a 3x3 matrix - for (i, channel) in enumerate([channel for channel in periodicElectricalTxChannels(field) if length(arbitraryElectricalComponents(channel)) == 0]) - next = 2 - comps = periodicElectricalComponents(channel) - for j = 1:len - comp = nothing - if (i == j) - comp = comps[1] - elseif next <= length(comps) - comp = comps[next] - next+=1 +function calcControlMatrix(cont::AWControlSequence) + oldTx = zeros(ComplexF64, controlMatrixShape(cont)) + for (i, channel) in enumerate(getControlledChannels(cont)) + if cont.rfftIndices[i,end,1] + oldTx[i,1] = ustrip(u"V", offset(channel)) + end + for (j, comp) in enumerate(components(channel)) + if isa(comp, PeriodicElectricalComponent) + oldTx[i, cont.rfftIndices[i,j,:]] .= ustrip(u"V",amplitude(comp)) * exp(im*ustrip(u"rad",phase(comp)-pi/2)) # The phase given in the component is for a sine, but the FFT-phase uses a cosine + elseif isa(comp, ArbitraryElectricalComponent) + oldTx[i, cont.rfftIndices[i,j,:]] .= rfft(ustrip.(u"V",scaledValues(comp)))[2:sum(cont.rfftIndices[i,j,:])+1]/(0.5*2^14) # the buffer length should always be 2^14 currently end + end + end + return oldTx +end - val = 0.0 - if !isnothing(comp) - r = ustrip(u"V", amplitude(comp)) - angle = ustrip(u"rad", phase(comp)) - val = r*cos(angle) + r*sin(angle)*im - end - κ[i, j] = val + +# Convert New Tx from matrix in V to currSequence +function updateControlSequence!(cont::CrossCouplingControlSequence, newTx::Matrix) + for (i, channel) in enumerate(periodicElectricalTxChannels(cont.currSequence)) + for (j, comp) in enumerate(components(channel)) + amplitude!(comp, abs(newTx[i, j])*1.0u"V") + phase!(comp, angle(newTx[i, j])*1.0u"rad") end end - return κ -end - -function updateControlSequence!(cont::ControlSequence, newTx::Matrix) - # TxDAQController only works on one field atm (possible future TODO: update to multiple fields, matrix per field) - field = fields(cont.currSequence)[1] - for (i, channel) in enumerate([channel for channel in periodicElectricalTxChannels(field) if length(arbitraryElectricalComponents(channel)) == 0]) - comps = periodicElectricalComponents(channel) - j = 1 - for (k, comp) in enumerate(comps) - val = 0.0 - # First component is a diagonal entry from the matrix - if k == 1 - val = newTx[i, i] - # All other components are "in order" and skip the diagonal entry - else - if j == i - j+=1 - end - val = newTx[i, j] - j+=1 +end + +function updateControlSequence!(cont::AWControlSequence, newTx::Matrix) + for (i, channel) in enumerate(periodicElectricalTxChannels(cont.currSequence)) + if cont.rfftIndices[i,end,1] + offset!(channel, real(newTx[i,1])*1.0u"V") + end + for (j, comp) in enumerate(components(channel)) + if isa(comp, PeriodicElectricalComponent) + amplitude!(comp, abs.(newTx[i, cont.rfftIndices[i,j,:]])[]*1.0u"V") + phase!(comp, angle.(newTx[i, cont.rfftIndices[i,j,:]])[]*1.0u"rad"+(pi/2)u"rad") + elseif isa(comp, ArbitraryElectricalComponent) + spectrum = zeros(ComplexF64, 2^13+1) + spectrum[2:sum(cont.rfftIndices[i,j,:])+1] .= newTx[i, cont.rfftIndices[i,j,:]] + amplitude!(comp, 1.0u"V") + phase!(comp, 0.0u"rad") + values!(comp, irfft(spectrum, 2^14)*(0.5*2^14)) end - amplitude!(comp, abs(val)*1.0u"V") - phase!(comp, angle(val)*1.0u"rad") end end +end + +################################################################################# +########## Functions for checking the matrix representation for safety and plausibility +################################################################################# + +function calcExpectedField(tx::Matrix{<:Complex}, cont::CrossCouplingControlSequence) + dividers = divider.(getPrimaryComponents(cont)) + frequencies = ustrip(u"Hz", txBaseFrequency(cont.currSequence)) ./ dividers + calibFieldToVoltEstimate = [ustrip(u"V/T", chan.calibration(frequencies[i])) for (i,chan) in enumerate(getControlledDAQChannels(cont))] + B_fw = tx ./ calibFieldToVoltEstimate + return B_fw end -function checkFieldToVolt(oldTx, Γ, cont::ControlSequence, txCont::TxDAQController) - calibFieldToVoltEstimate = [ustrip(u"V/T", ch.calibration) for ch in collect(Base.values(cont.simpleChannel))] - calibFieldToVoltMeasured = abs.(diag(oldTx) ./ diag(Γ)) - deviation = abs.(1.0 .- calibFieldToVoltMeasured./calibFieldToVoltEstimate) - @debug "We expected $(calibFieldToVoltEstimate) and got $(calibFieldToVoltMeasured), deviation: $deviation" - valid = maximum( deviation ) < txCont.params.fieldToVoltDeviation - if !valid - @warn "Measured field to volt deviates by $deviation from estimate, exceeding allowed deviation" - end - return valid +function calcExpectedField(tx::Matrix{<:Complex}, cont::AWControlSequence) + N = rxNumSamplingPoints(cont.currSequence) + frequencies = ustrip.(u"Hz",rfftfreq(N, rxSamplingRate(cont.currSequence))) + calibFieldToVoltEstimate = reduce(vcat,transpose([ustrip.(u"V/T", chan.calibration(frequencies)) for chan in getControlledDAQChannels(cont)])) + B_fw = tx ./ calibFieldToVoltEstimate + return B_fw end -function checkVoltLimits(newTx, cont::ControlSequence, txCont::TxDAQController) - validChannel = zeros(Bool, size(newTx, 2)) - for i = 1:size(newTx, 2) +function validateAgainstForwardCalibrationAndSafetyLimit(tx::Matrix{<:Complex}, B::Matrix{<:Complex}, cont::ControlSequence, txCont::TxDAQController) + # step 1 apply forward calibration to tx -> B_fw + B_fw = calcExpectedField(tx, cont) + + # step 2 check B_fw against B (rel. and abs. Accuracy) + forwardCalibrationAgrees = isapprox.(abs.(B_fw), abs.(B), rtol = txCont.params.fieldToVoltRelDeviation, atol=ustrip(u"T",txCont.params.fieldToVoltAbsDeviation)) + + # step 3 check if B_fw and B are both below safety limit + isSafe(Btest) = abs.(Btest). u"V/µs" + validSlew = maximum(abs.(slew_rate), dims=2) .< getproperty.(getControlledDAQChannels(cont),:limitSlewRate) + validPeak = maximum(abs.(testSignalTime), dims=2) .< ustrip.(u"V", getproperty.(getControlledDAQChannels(cont),:limitPeak)) + + valid = all(validSlew) && all(validPeak) + @debug "Check Volt Limit" p=lineplot(1:N,testSignalTime', canvas=DotCanvas, border=:ascii) maximum(abs.(testSignalTime), dims=2) maximum(abs.(slew_rate), dims=2) + if !valid + @debug "Valid Tx Channel" validSlew validPeak + @warn "New control sequence exceeds voltage limits (slew rate or peak) of tx channel" + end + return valid + end +end diff --git a/src/MPIMeasurements.jl b/src/MPIMeasurements.jl index 4f72c1fd..3ec53470 100644 --- a/src/MPIMeasurements.jl +++ b/src/MPIMeasurements.jl @@ -15,17 +15,23 @@ using InteractiveUtils using Mmap using Scratch using StringEncodings +using Statistics using DocStringExtensions using MacroTools using LibSerialPort +using UnicodePlots using LinearAlgebra +using HDF5 +using FFTW +using NaturalSort +using RelocatableFolders using ReplMaker import REPL import REPL: LineEdit, REPLCompletions import REPL: TerminalMenus import Base.write, Base.take!, Base.put!, Base.isready, Base.isopen, Base.eltype, Base.close, Base.wait, Base.length, Base.push! -import Base: ==, isequal, hash, isfile, push!, pop!, empty!, getindex, setindex!, firstindex, lastindex, length, iterate, delete!, deleteat!, keys, haskey +import Base: ==, isequal, hash, isfile, push!, pop!, empty!, getindex, setindex!, firstindex, lastindex, length, iterate, delete!, deleteat!, keys, haskey, values # Reexporting MPIFiles is disliked by Aqua since there are undefined exports. Therefore, I disabled reexporting here. #using Reexport @@ -39,11 +45,11 @@ import MPIFiles: hasKeyAndValue, rxBandwidth, rxNumChannels, rxNumSamplingPoints using RedPitayaDAQServer - -const scannerConfigurationPath = [normpath(string(@__DIR__), "../config")] # Push custom configuration directories here +const DEFAULT_SCANNER_CONFIG_PATH = @path joinpath(@__DIR__, "..", "config") +const scannerConfigurationPath = AbstractString[DEFAULT_SCANNER_CONFIG_PATH] # Push custom configuration directories here export addConfigurationPath -addConfigurationPath(path::String) = !(path in scannerConfigurationPath) ? pushfirst!(scannerConfigurationPath, path) : nothing +addConfigurationPath(path) = !(path in scannerConfigurationPath) ? pushfirst!(scannerConfigurationPath, path) : nothing # Circular reference between Scanner.jl and Protocol.jl. Thus we predefine the protocol """ diff --git a/src/Protocols/ContinousMeasurementProtocol.jl b/src/Protocols/ContinousMeasurementProtocol.jl index 38cac750..49318d3e 100644 --- a/src/Protocols/ContinousMeasurementProtocol.jl +++ b/src/Protocols/ContinousMeasurementProtocol.jl @@ -1,6 +1,8 @@ export ContinousMeasurementProtocol, ContinousMeasurementProtocolParams """ -Parameters for the MPIMeasurementProtocol +Parameters for the `ContinousMeasurementProtocol`` + +$FIELDS """ Base.@kwdef mutable struct ContinousMeasurementProtocolParams <: ProtocolParams "Foreground frames to measure. Overwrites sequence frames" @@ -62,12 +64,7 @@ function _init(protocol::ContinousMeasurementProtocol) protocol.cancelled = false protocol.finishAcknowledged = false if protocol.params.controlTx - controllers = getDevices(protocol.scanner, TxDAQController) - if length(controllers) > 1 - throw(IllegalStateException("Cannot unambiguously find a TxDAQController as the scanner has $(length(controllers)) of them")) - end - protocol.txCont = controllers[1] - protocol.txCont.currTx = nothing + protocol.txCont = getDevice(protocol.scanner, TxDAQController) else protocol.txCont = nothing end @@ -122,7 +119,7 @@ function _execute(protocol::ContinousMeasurementProtocol) while !measPauseOver || protocol.stopped handleEvents(protocol) if !notifiedStop && protocol.stopped - put!(protocol.biChannel, OperationSuccessfulEvent(StopEvent())) + put!(protocol.biChannel, OperationSuccessfulEvent(PauseEvent())) notifiedStop = true end if notifiedStop && !protocol.stopped @@ -198,11 +195,7 @@ function asyncMeasurement(protocol::ContinousMeasurementProtocol) return protocol.seqMeasState end -function cleanup(protocol::ContinousMeasurementProtocol) - # NOP -end - -function stop(protocol::ContinousMeasurementProtocol) +function pause(protocol::ContinousMeasurementProtocol) protocol.stopped = true end diff --git a/src/Protocols/MPIForceProtocol.jl b/src/Protocols/MPIForceProtocol.jl index 4aa0ad94..9954bad2 100644 --- a/src/Protocols/MPIForceProtocol.jl +++ b/src/Protocols/MPIForceProtocol.jl @@ -1,6 +1,8 @@ export MPIForceProtocol, MPIForceProtocolParams """ -Parameters for the MPIForceProtocol +Parameters for the `MPIForceProtocol` + +$FIELDS """ Base.@kwdef mutable struct MPIForceProtocolParams <: ProtocolParams "If set the tx amplitude and phase will be set with control steps" @@ -166,19 +168,10 @@ function performExperiment(protocol::MPIForceProtocol) end end - -function cleanup(protocol::MPIForceProtocol) - # NOP -end - -function stop(protocol::MPIForceProtocol) +function pause(protocol::MPIForceProtocol) protocol.stopped = true end -function resume(protocol::MPIForceProtocol) - put!(protocol.biChannel, OperationNotSupportedEvent(ResumeEvent())) -end - function cancel(protocol::MPIForceProtocol) protocol.cancelled = true end diff --git a/src/Protocols/MPIMeasurementProtocol.jl b/src/Protocols/MPIMeasurementProtocol.jl index cad1647c..ad8d1e99 100644 --- a/src/Protocols/MPIMeasurementProtocol.jl +++ b/src/Protocols/MPIMeasurementProtocol.jl @@ -1,6 +1,8 @@ export MPIMeasurementProtocol, MPIMeasurementProtocolParams """ -Parameters for the MPIMeasurementProtocol +Parameters for the `MPIMeasurementProtocol` + +$FIELDS """ Base.@kwdef mutable struct MPIMeasurementProtocolParams <: ProtocolParams "Foreground frames to measure. Overwrites sequence frames" @@ -15,8 +17,8 @@ Base.@kwdef mutable struct MPIMeasurementProtocolParams <: ProtocolParams saveTemperatureData::Bool = false "Sequence to measure" sequence::Union{Sequence, Nothing} = nothing - "Remember background measurement" - rememberBGMeas::Bool = false + "Do not measure background but reuse the last BG measurement if suitable" + reuseLastBGMeas::Bool = false end function MPIMeasurementProtocolParams(dict::Dict, scanner::MPIScanner) sequence = nothing @@ -78,12 +80,11 @@ function timeEstimate(protocol::MPIMeasurementProtocol) if !isnothing(protocol.params.sequence) params = protocol.params seq = params.sequence - totalFrames = (params.fgFrames + params.bgFrames) * acqNumFrameAverages(seq) + totalFrames = (params.fgFrames + params.bgFrames*params.measureBackground*!params.reuseLastBGMeas) * acqNumFrameAverages(seq) samplesPerFrame = rxNumSamplingPoints(seq) * acqNumAverages(seq) * acqNumPeriodsPerFrame(seq) totalTime = (samplesPerFrame * totalFrames) / (125e6/(txBaseFrequency(seq)/rxSamplingRate(seq))) time = totalTime * 1u"s" est = string(time) - @show est end return est end @@ -114,7 +115,7 @@ function _execute(protocol::MPIMeasurementProtocol) end function performMeasurement(protocol::MPIMeasurementProtocol) - if (length(protocol.bgMeas) == 0 || !protocol.params.rememberBGMeas) && protocol.params.measureBackground + if (length(protocol.bgMeas) == 0 || !protocol.params.reuseLastBGMeas) && protocol.params.measureBackground if askChoices(protocol, "Press continue when background measurement can be taken", ["Cancel", "Continue"]) == 1 throw(CancelException()) end @@ -123,9 +124,9 @@ function performMeasurement(protocol::MPIMeasurementProtocol) @debug "Taking background measurement." protocol.unit = "BG Frames" measurement(protocol) - protocol.bgMeas = read(sink(protocol.seqMeasState.sequenceBuffer, MeasurementBuffer)) deviceBuffers = protocol.seqMeasState.deviceBuffers push!(protocol.protocolMeasState, vcat(sinks(protocol.seqMeasState.sequenceBuffer), isnothing(deviceBuffers) ? SinkBuffer[] : deviceBuffers), isBGMeas = true) + protocol.bgMeas = read(protocol.protocolMeasState, MeasurementBuffer) if askChoices(protocol, "Press continue when foreground measurement can be taken", ["Cancel", "Continue"]) == 1 throw(CancelException()) end @@ -205,19 +206,6 @@ function asyncMeasurement(protocol::MPIMeasurementProtocol) return protocol.seqMeasState end - -function cleanup(protocol::MPIMeasurementProtocol) - # NOP -end - -function stop(protocol::MPIMeasurementProtocol) - put!(protocol.biChannel, OperationNotSupportedEvent(StopEvent())) -end - -function resume(protocol::MPIMeasurementProtocol) - put!(protocol.biChannel, OperationNotSupportedEvent(ResumeEvent())) -end - function cancel(protocol::MPIMeasurementProtocol) protocol.cancelled = true #put!(protocol.biChannel, OperationNotSupportedEvent(CancelEvent())) @@ -269,6 +257,14 @@ function handleEvent(protocol::MPIMeasurementProtocol, event::DatasetStoreStorag mdf = event.mdf data = read(protocol.protocolMeasState, MeasurementBuffer) isBGFrame = measIsBGFrame(protocol.protocolMeasState) + if protocol.params.reuseLastBGMeas + try + data = cat(data,protocol.bgMeas,dims=4) + isBGFrame = cat(isBGFrame,trues(size(protocol.bgMeas,4)), dims=1) + catch + @error "The last BG measurement is not compatible with your current measurement, cant write BG into file! Data is size $(size(data)), BG is size $(size(protocol.bgMeas))" + end + end drivefield = read(protocol.protocolMeasState, DriveFieldBuffer) appliedField = read(protocol.protocolMeasState, TxDAQControllerBuffer) temperature = read(protocol.protocolMeasState, TemperatureBuffer) diff --git a/src/Protocols/MPSMeasurementProtocol.jl b/src/Protocols/MPSMeasurementProtocol.jl index bbf93fd4..dde61073 100644 --- a/src/Protocols/MPSMeasurementProtocol.jl +++ b/src/Protocols/MPSMeasurementProtocol.jl @@ -1,6 +1,8 @@ export MPSMeasurementProtocol, MPSMeasurementProtocolParams """ -Parameters for the MPSMeasurementProtocol +Parameters for the `MPSMeasurementProtocol` + + $FIELDS """ Base.@kwdef mutable struct MPSMeasurementProtocolParams <: ProtocolParams "Foreground frames to measure. Overwrites sequence frames" @@ -10,27 +12,23 @@ Base.@kwdef mutable struct MPSMeasurementProtocolParams <: ProtocolParams "If set the tx amplitude and phase will be set with control steps" controlTx::Bool = false "If unset no background measurement will be taken" - measureBackground::Bool = true + measureBackground::Bool = false "Remember background measurement" rememberBGMeas::Bool = false - "Tracer that is being used for the measurement" + #"Tracer that is being used for the measurement" #tracer::Union{Tracer, Nothing} = nothing - "If the temperature should be safed or not" + "If the temperature should be saved or not" saveTemperatureData::Bool = false "Sequence to measure" sequence::Union{Sequence, Nothing} = nothing - - # TODO: This is only for 1D MPS systems for now - "Start value of the MPS offset measurement. Overwrites parts of the sequence definition." - offsetStart::typeof(1.0u"T") = -0.012u"T" - "Stop value of the MPS offset measurement. Overwrites parts of the sequence definition." - offsetStop::typeof(1.0u"T") = 0.012u"T" - "Number of values of the MPS offset measurement. Overwrites parts of the sequence definition." - offsetNum::Integer = 101 + "Sort patches" + sortPatches::Bool = true + "Flag if the measurement should be saved as a system matrix or not" + saveInCalibFolder::Bool = false "Number of periods per offset of the MPS offset measurement. Overwrites parts of the sequence definition." - dfPeriodsPerOffset::Integer = 100 - "Number of periods per offset which should be deleted. Acquired total number of periods is `dfPeriodsPerOffset + deletedDfPeriodsPerOffset`. Overwrites parts of the sequence definition." - deletedDfPeriodsPerOffset::Integer = 1 + dfPeriodsPerOffset::Integer = 50 + "If true all periods per offset are averaged" + averagePeriodsPerOffset::Bool = true end function MPSMeasurementProtocolParams(dict::Dict, scanner::MPIScanner) sequence_ = nothing @@ -62,6 +60,11 @@ Base.@kwdef mutable struct MPSMeasurementProtocol <: Protocol seqMeasState::Union{SequenceMeasState, Nothing} = nothing protocolMeasState::Union{ProtocolMeasState, Nothing} = nothing + sequence::Union{Sequence, Nothing} = nothing + offsetfields::Union{Matrix{Float64}, Nothing} = nothing + patchPermutation::Vector{Union{Int64, Nothing}} = Union{Int64, Nothing}[] + calibsize::Vector{Int64} = Int64[] + bgMeas::Array{Float32, 4} = zeros(Float32,0,0,0,0) done::Bool = false cancelled::Bool = false @@ -80,9 +83,9 @@ function requiredDevices(protocol::MPSMeasurementProtocol) end function _init(protocol::MPSMeasurementProtocol) - if isnothing(sequence(protocol)) - throw(IllegalStateException("Protocol requires a sequence")) - end + #if isnothing(sequence(protocol)) + # throw(IllegalStateException("Protocol requires a sequence")) + #end protocol.done = false protocol.cancelled = false protocol.finishAcknowledged = false @@ -94,7 +97,32 @@ function _init(protocol::MPSMeasurementProtocol) protocol.bgMeas = zeros(Float32,0,0,0,0) protocol.protocolMeasState = ProtocolMeasState() - setupSequence(protocol) + #try + seq, perm, offsets, calibsize, numPeriodsPerFrame = prepareProtocolSequences(protocol.params.sequence, getDAQ(scanner(protocol)); numPeriodsPerOffset = protocol.params.dfPeriodsPerOffset) + + # For each patch assign nothing if invalid or otherwise index in "proper" frame + patchPerm = Vector{Union{Int64, Nothing}}(nothing, numPeriodsPerFrame) + if !protocol.params.sortPatches + # perm is arranged in a way that the first offset dimension switches the fastest + # if the patches should be saved in the order they were measured, we need to sort perm + perm = sort(perm) + end + + # patchPerm contains the "target" for every patch that is measured, nothing for discarded patches, different indizes for non-averaged patches, identical indizes for averaged-patches + # Same target for all frames to be averaged + part = protocol.params.averagePeriodsPerOffset ? protocol.params.dfPeriodsPerOffset : 1 + for (i, patches) in enumerate(Iterators.partition(perm, part)) + patchPerm[patches] .= i + end + + protocol.sequence = seq + protocol.patchPermutation = patchPerm + protocol.offsetfields = ustrip.(u"T", offsets) # TODO make robust + protocol.calibsize = calibsize + @debug "Prepared Protocol Sequence: $(length(patchPerm)) measured and $(length(perm)) valid patches in Permutation" + #catch e + # throw(e) + #end return nothing end @@ -103,13 +131,15 @@ function timeEstimate(protocol::MPSMeasurementProtocol) est = "Unknown" if !isnothing(sequence(protocol)) params = protocol.params - seq = params.sequence - totalFrames = (params.fgFrames + params.bgFrames) * acqNumFrameAverages(seq) - samplesPerFrame = rxNumSamplingPoints(seq) * acqNumAverages(seq) * acqNumPeriodsPerFrame(seq) - totalTime = (samplesPerFrame * totalFrames) / (125e6/(txBaseFrequency(seq)/rxSamplingRate(seq))) - time = totalTime * 1u"s" - est = string(time) - @info "The estimated duration is $est s." + seq = sequence(protocol) + fgFrames = params.fgFrames * acqNumFrameAverages(seq) + bgFrames = params.measureBackground*params.bgFrames * acqNumFrameAverages(seq) + txSamplesPerFrame = lcm(dfDivider(seq)) * size(protocol.patchPermutation, 1) + fgTime = (txSamplesPerFrame * fgFrames) / txBaseFrequency(seq) |> u"s" + bgTime = (txSamplesPerFrame * bgFrames) / txBaseFrequency(seq) |> u"s" + perc_wait = round(Int,sum(isnothing.(protocol.patchPermutation))/size(protocol.patchPermutation,1)*100) + est = "FG: $(timeFormat(fgTime)) ($(perc_wait)% waiting)\nBG: $(timeFormat(bgTime))" + @info "The estimated duration is FG: $fgTime ($(perc_wait)% waiting), BG: $bgTime." end return est end @@ -120,68 +150,9 @@ function enterExecute(protocol::MPSMeasurementProtocol) protocol.finishAcknowledged = false end -function getRelevantDfDivider(protocol::MPSMeasurementProtocol) - sequence_ = sequence(protocol) - - divider = nothing - offsetFieldIdx = nothing - offsetChannelIdx = nothing - for (fieldIdx, field) in enumerate(fields(sequence_)) - for (channelIdx, channel) in enumerate(channels(field)) - if channel isa PeriodicElectricalChannel - @warn "The protocol currently always uses the first component of $(id(channel)) for determining the divider." - divider = channel.components[1].divider - elseif channel isa ContinuousElectricalChannel - offsetFieldIdx = fieldIdx - offsetChannelIdx = channelIdx - end - end - end - - if isnothing(divider) - throw(ProtocolConfigurationError("The sequence `$(name(sequence_))` for the protocol `$(name(protocol))` does not define a PeriodicElectricalChannel and thus a divider.")) - end - - if isnothing(offsetFieldIdx) || isnothing(offsetChannelIdx) - throw(ProtocolConfigurationError("The sequence `$(name(sequence_))` for the protocol `$(name(protocol))` does not define a ContinuousElectricalChannel and thus an offset field description.")) - end - - return divider, offsetFieldIdx, offsetChannelIdx -end - -function createOffsetChannel(protocol::MPSMeasurementProtocol; deletedPeriodsPerOffset=protocol.params.deletedDfPeriodsPerOffset) - sequence_ = sequence(protocol) - - divider, offsetFieldIdx, offsetChannelIdx = getRelevantDfDivider(protocol) - - stepDivider = divider * (protocol.params.dfPeriodsPerOffset + deletedPeriodsPerOffset) - offsetDivider = stepDivider * protocol.params.offsetNum - - chanOffset = (protocol.params.offsetStop + protocol.params.offsetStart) / 2 - amplitude = abs(protocol.params.offsetStop - protocol.params.offsetStart) / 2 - - oldChannel = sequence_.fields[offsetFieldIdx].channels[offsetChannelIdx] - newChannel = ContinuousElectricalChannel(id=oldChannel.id, dividerSteps=stepDivider, divider=offsetDivider, amplitude=amplitude, phase=oldChannel.phase, offset=chanOffset, waveform=WAVEFORM_SAWTOOTH_RISING) - - @info "Values used for the creation of the offset channel" divider stepDivider offsetDivider chanOffset amplitude - - return newChannel, offsetFieldIdx, offsetChannelIdx -end - -function setupSequence(protocol::MPSMeasurementProtocol; deletedPeriodsPerOffset=protocol.params.deletedDfPeriodsPerOffset) - sequence_ = sequence(protocol) - newChannel, offsetFieldIdx, offsetChannelIdx = createOffsetChannel(protocol, deletedPeriodsPerOffset=deletedPeriodsPerOffset) - - @info "Offset values used for the measurement" MPIMeasurements.values(newChannel) - - sequence_.fields[offsetFieldIdx].channels[offsetChannelIdx] = newChannel - protocol.params.sequence = sequence_ -end - function _execute(protocol::MPSMeasurementProtocol) @debug "Measurement protocol started" - setupSequence(protocol) performMeasurement(protocol) put!(protocol.biChannel, FinishedNotificationEvent()) @@ -217,7 +188,7 @@ function performMeasurement(protocol::MPSMeasurementProtocol) acqNumFrames(sequence(protocol), protocol.params.fgFrames) @debug "Starting foreground measurement." - protocol.unit = "Frames" + protocol.unit = "Offsets" measurement(protocol) deviceBuffers = protocol.seqMeasState.deviceBuffers push!(protocol.protocolMeasState, vcat(sinks(protocol.seqMeasState.sequenceBuffer), isnothing(deviceBuffers) ? SinkBuffer[] : deviceBuffers), isBGMeas = false) @@ -267,38 +238,92 @@ function measurement(protocol::MPSMeasurementProtocol) end function asyncMeasurement(protocol::MPSMeasurementProtocol) - scanner_ = scanner(protocol) - sequence = protocol.params.sequence - daq = getDAQ(scanner_) + scanner_ = scanner(protocol) + sequence, protocol.seqMeasState = SequenceMeasState(protocol) + protocol.seqMeasState.producer = @tspawnat scanner_.generalParams.producerThreadID asyncProducer(protocol.seqMeasState.channel, protocol, sequence) + bind(protocol.seqMeasState.channel, protocol.seqMeasState.producer) + protocol.seqMeasState.consumer = @tspawnat scanner_.generalParams.consumerThreadID asyncConsumer(protocol.seqMeasState) + return protocol.seqMeasState +end + +function SequenceMeasState(protocol::MPSMeasurementProtocol) + sequence = protocol.sequence + daq = getDAQ(scanner(protocol)) deviceBuffer = DeviceBuffer[] + + + # Setup everything as defined per sequence if protocol.params.controlTx sequence = controlTx(protocol.txCont, sequence) push!(deviceBuffer, TxDAQControllerBuffer(protocol.txCont, sequence)) end setup(daq, sequence) - protocol.seqMeasState = SequenceMeasState(daq, sequence) + + # Now for the buffer chain we want to reinterpret periods to frames + # This has to happen after the RedPitaya sequence is set, as that code repeats the full sequence for each frame + acqNumFrames(protocol.sequence, acqNumFrames(protocol.sequence) * periodsPerFrame(daq.rpc)) + setupRx(daq, daq.decimation, samplesPerPeriod(daq.rpc), 1) + + numFrames = acqNumFrames(protocol.sequence) + + # Prepare buffering structures: + # RedPitaya-> MPSBuffer -> Splitter --> FrameBuffer{Mmap} + # |-> DriveFieldBuffer + numValidPatches = length(filter(x->!isnothing(x), protocol.patchPermutation)) + averages = protocol.params.averagePeriodsPerOffset ? protocol.params.dfPeriodsPerOffset : 1 + numFrames = div(numValidPatches, averages) + bufferSize = (rxNumSamplingPoints(protocol.sequence), length(rxChannels(protocol.sequence)), 1, numFrames) + buffer = FrameBuffer(protocol, "meas.bin", Float32, bufferSize) + + # buffers = StorageBuffer[buffer] + # TODO: reenable DriveFieldBuffer + #if protocol.params.controlTx + #len = length(keys(sequence.refIndices)) + # push!(buffers, DriveFieldBuffer(1, zeros(ComplexF64, controlMatrixShape(sequence)..., 1, numFrames), sequence)) + #end + + # buffer = FrameSplitterBuffer(daq, StorageBuffer[buffer]) + buffer = MPSBuffer(buffer, protocol.patchPermutation, averages, 1, acqNumPeriodsPerFrame(protocol.sequence)) + + channel = Channel{channelType(daq)}(32) + deviceBuffer = DeviceBuffer[] if protocol.params.saveTemperatureData - push!(deviceBuffer, TemperatureBuffer(getTemperatureSensor(scanner_), acqNumFrames(protocol.params.sequence))) + push!(deviceBuffer, TemperatureBuffer(getTemperatureSensor(scanner(protocol)), numFrames)) end - protocol.seqMeasState.deviceBuffers = deviceBuffer - protocol.seqMeasState.producer = @tspawnat scanner_.generalParams.producerThreadID asyncProducer(protocol.seqMeasState.channel, protocol, sequence) - bind(protocol.seqMeasState.channel, protocol.seqMeasState.producer) - protocol.seqMeasState.consumer = @tspawnat scanner_.generalParams.consumerThreadID asyncConsumer(protocol.seqMeasState) - return protocol.seqMeasState -end - -function cleanup(protocol::MPSMeasurementProtocol) - # NOP + return sequence, SequenceMeasState(numFrames, channel, nothing, nothing, AsyncBuffer(buffer, daq), deviceBuffer, asyncMeasType(protocol.sequence)) end -function stop(protocol::MPSMeasurementProtocol) - put!(protocol.biChannel, OperationNotSupportedEvent(StopEvent())) +mutable struct MPSBuffer <: IntermediateBuffer + target::StorageBuffer + permutation::Vector{Union{Int64, Nothing}} + average::Int64 + counter::Int64 + total::Int64 end - -function resume(protocol::MPSMeasurementProtocol) - put!(protocol.biChannel, OperationNotSupportedEvent(ResumeEvent())) +function push!(mpsBuffer::MPSBuffer, frames::Array{T,4}) where T + from = nothing + to = nothing + for i = 1:size(frames, 4) + frameIdx = div(mpsBuffer.counter - 1, mpsBuffer.total) + patchCounter = mod1(mpsBuffer.counter, mpsBuffer.total) + patchIdx = mpsBuffer.permutation[patchCounter] + if !isnothing(patchIdx) + if mpsBuffer.average > 1 + to = insert!(+, mpsBuffer.target, frameIdx * mpsBuffer.total + patchIdx, view(frames, :, :, :, i:i)./mpsBuffer.average) + else + to = insert!(mpsBuffer.target, frameIdx * mpsBuffer.total + patchIdx, view(frames, :, :, :, i:i)) + end + end + mpsBuffer.counter += 1 + end + if !isnothing(from) && !isnothing(to) + return (start = from, stop = to) + else + return nothing + end end +sinks!(buffer::MPSBuffer, sinks::Vector{SinkBuffer}) = sinks!(buffer.target, sinks) function cancel(protocol::MPSMeasurementProtocol) protocol.cancelled = true @@ -335,7 +360,8 @@ function handleEvent(protocol::MPSMeasurementProtocol, event::ProgressQueryEvent reply = nothing if !isnothing(protocol.seqMeasState) framesTotal = protocol.seqMeasState.numFrames - framesDone = min(index(sink(protocol.seqMeasState.sequenceBuffer, MeasurementBuffer)) - 1, framesTotal) + measFramesDone = protocol.seqMeasState.sequenceBuffer.target.counter-1 + framesDone = length(unique(protocol.patchPermutation[1:measFramesDone]))-1 reply = ProgressEvent(framesDone, framesTotal, protocol.unit, event) else reply = ProgressEvent(0, 0, "N/A", event) @@ -348,30 +374,47 @@ handleEvent(protocol::MPSMeasurementProtocol, event::FinishedAckEvent) = protoco function handleEvent(protocol::MPSMeasurementProtocol, event::DatasetStoreStorageRequestEvent) store = event.datastore scanner = protocol.scanner + sequence = protocol.sequence mdf = event.mdf - data = read(protocol.protocolMeasState, MeasurementBuffer) - - if protocol.params.deletedDfPeriodsPerOffset > 0 - numSamples_ = size(data, 1) - numChannels_ = size(data, 2) - numPeriods_ = size(data, 3) - numFrames_ = size(data, 4) - numPeriodsPerOffset_ = div(numPeriods_, protocol.params.offsetNum) - - data = reshape(data, (numSamples_, numChannels_, numPeriodsPerOffset_, protocol.params.offsetNum, numFrames_)) - data = data[:, :, protocol.params.deletedDfPeriodsPerOffset+1:end, :, :] # Kick out first N periods - data = reshape(data, (numSamples_, numChannels_, :, numFrames_)) + data = read(protocol.protocolMeasState, MeasurementBuffer) + data = reshape(data, rxNumSamplingPoints(sequence), length(rxChannels(sequence)), :, protocol.params.fgFrames + protocol.params.bgFrames * protocol.params.measureBackground) + acqNumFrames(sequence, size(data, 4)) + + periodsPerOffset = size(protocol.patchPermutation,1)÷size(protocol.offsetfields,1) - # Reset sequence since the info is used for the MDF - setupSequence(protocol, deletedPeriodsPerOffset=0) + offsetPerm = zeros(Int64, size(data, 3)) + for (index, patch) in enumerate(protocol.patchPermutation) + if !isnothing(patch) + offsetPerm[patch] = (index-1)÷periodsPerOffset + 1 + end end + offsets = protocol.offsetfields[offsetPerm, :] - isBGFrame = measIsBGFrame(protocol.protocolMeasState) + + isBGFrame = reduce(&, reshape(measIsBGFrame(protocol.protocolMeasState), size(data, 3), :), dims = 1)[1, :] drivefield = read(protocol.protocolMeasState, DriveFieldBuffer) appliedField = read(protocol.protocolMeasState, TxDAQControllerBuffer) temperature = read(protocol.protocolMeasState, TemperatureBuffer) - filename = saveasMDF(store, scanner, protocol.params.sequence, data, isBGFrame, mdf, drivefield = drivefield, temperatures = temperature, applied = appliedField) + + + filename = nothing + + periodsPerOffset = protocol.params.averagePeriodsPerOffset ? 1 : protocol.params.dfPeriodsPerOffset + isBGFrame = repeat(isBGFrame, inner = div(size(data, 3), periodsPerOffset)) + data = reshape(data, size(data, 1), size(data, 2), periodsPerOffset, :) + + # All periods in one frame (should) have same offset + offsets = reshape(offsets, (periodsPerOffset, :, size(offsets, 2)))[1, :, :] + + # Pad zeros in z or z and y direction if the protocol uses less than 3 offset dimensions. + @assert size(offsets, 2) <= 3 "The number of offset dimensions must be below or equal to 3 since the MDF specifies Ox3 with directions in x, y and z." + offsets = hcat(offsets, zeros(eltype(offsets), (size(offsets, 1), 3 - size(offsets, 2)))) + + calibsize_ = vcat(protocol.calibsize, ones(eltype(protocol.calibsize), 3 - length(protocol.calibsize))) + offsets = reshape(offsets, (calibsize_..., :)) # make calib size "visible" to storing function + filename = saveasMDF(store, scanner, sequence, data, offsets, isBGFrame, mdf, storeAsSystemMatrix=protocol.params.saveInCalibFolder, drivefield = drivefield, temperatures = temperature, applied = appliedField) + @info "The measurement was saved at `$filename`." put!(protocol.biChannel, StorageSuccessEvent(filename)) end @@ -379,4 +422,4 @@ end protocolInteractivity(protocol::MPSMeasurementProtocol) = Interactive() protocolMDFStudyUse(protocol::MPSMeasurementProtocol) = UsingMDFStudy() -sequence(protocol::MPSMeasurementProtocol) = protocol.params.sequence \ No newline at end of file +sequence(protocol::MPSMeasurementProtocol) = protocol.sequence \ No newline at end of file diff --git a/src/Protocols/MechanicalMPIMeasurementProtocol.jl b/src/Protocols/MechanicalMPIMeasurementProtocol.jl index b9015edc..87cb7eb0 100644 --- a/src/Protocols/MechanicalMPIMeasurementProtocol.jl +++ b/src/Protocols/MechanicalMPIMeasurementProtocol.jl @@ -1,6 +1,8 @@ export MechanicalMPIMeasurementProtocol, MechanicalMPIMeasurementProtocolParams """ -Parameters for the MechanicalMPIMeasurementProtocol +Parameters for the `MechanicalMPIMeasurementProtocol` + +$FIELDS """ Base.@kwdef mutable struct MechanicalMPIMeasurementProtocolParams <: ProtocolParams "Foreground frames to measure. Overwrites sequence frames" @@ -215,18 +217,6 @@ function asyncMeasurement(protocol::MechanicalMPIMeasurementProtocol) return protocol.seqMeasState end -function cleanup(protocol::MechanicalMPIMeasurementProtocol) - # NOP -end - -function stop(protocol::MechanicalMPIMeasurementProtocol) - put!(protocol.biChannel, OperationNotSupportedEvent(StopEvent())) -end - -function resume(protocol::MechanicalMPIMeasurementProtocol) - put!(protocol.biChannel, OperationNotSupportedEvent(ResumeEvent())) -end - function cancel(protocol::MechanicalMPIMeasurementProtocol) protocol.cancelled = true #put!(protocol.biChannel, OperationNotSupportedEvent(CancelEvent())) diff --git a/src/Protocols/MultiSequenceSystemMatrixProtocol.jl b/src/Protocols/MultiSequenceSystemMatrixProtocol.jl new file mode 100644 index 00000000..16b63b5b --- /dev/null +++ b/src/Protocols/MultiSequenceSystemMatrixProtocol.jl @@ -0,0 +1,581 @@ +export MultiSequenceSystemMatrixProtocol, MultiSequenceSystemMatrixProtocolParams +""" +Parameters for the `MultiSequenceSystemMatrixProtocol` + +$FIELDS +""" +Base.@kwdef mutable struct MultiSequenceSystemMatrixProtocolParams <: ProtocolParams + "If set the tx amplitude and phase will be set with control steps" + controlTx::Bool = false + "If the temperature should be safed or not" + saveTemperatureData::Bool = false + "Sequences to measure" + sequences::Union{Vector{Sequence},Nothing} = nothing + "SM Positions mapped to the natural sorting of sequence tomls" + positions::Union{Positions, Nothing} = nothing + "Flag if the calibration should be saved as a system matrix or not" + saveInCalibFolder::Bool = false + "Seconds to wait between measurements" + waitTime::Float64 = 0.0 + "Maximum temperature the coils can be at to start the next measurement" + coolDownTo::Float64 = 100 + "Sequence to use as a BG measurement" + bgSequence::Union{Sequence, Nothing} = nothing + "Number of BG measurements before and after the measurement" + numBGMeas::Int = 1 +end + +function MultiSequenceSystemMatrixProtocolParams(dict::Dict, scanner::MPIScanner) + positions = nothing + if haskey(dict, "Positions") + posDict = dict["Positions"] + + positions = Positions(posDict) + delete!(dict, "Positions") + end + + if haskey(dict, "sequences") + sequences = Sequences(scanner, dict["sequences"]) + pop!(dict, "sequences") + end + + if haskey(dict, "bgSequence") + bgSequence = Sequence(scanner, dict["bgSequence"]) + pop!(dict, "bgSequence") + else + bgSequence = nothing + end + params = params_from_dict(MultiSequenceSystemMatrixProtocolParams, dict) + params.sequences = sequences + params.positions = positions + params.bgSequence = bgSequence + return params +end +MultiSequenceSystemMatrixProtocolParams(dict::Dict) = params_from_dict(MultiSequenceSystemMatrixProtocolParams, dict) + +Base.@kwdef mutable struct MultiSequenceSystemMatrixProtocol <: Protocol + @add_protocol_fields MultiSequenceSystemMatrixProtocolParams + + systemMeasState::Union{SystemMatrixMeasState,Nothing} = nothing + allSequences::Union{Vector{Sequence}, Nothing} = nothing + + done::Bool = false + cancelled::Bool = false + finishAcknowledged::Bool = false + paused::Bool = false + restored::Bool = false + measuring::Bool = false + txCont::Union{TxDAQController,Nothing} = nothing +end + +function requiredDevices(protocol::MultiSequenceSystemMatrixProtocol) + result = [AbstractDAQ] + if protocol.params.controlTx + push!(result, TxDAQController) + end + if protocol.params.saveTemperatureData + push!(result, TemperatureSensor) + end + return result +end + +function _init(protocol::MultiSequenceSystemMatrixProtocol) + if isnothing(protocol.params.sequences) + throw(IllegalStateException("Protocol requires sequences")) + end + + if any(acqNumFrameAverages.(protocol.params.sequences) .!= 1) + throw(ProtocolConfigurationError("The sequences for a MultiSequenceSystemMatrixProtocol currently do not support numFrameAverages != 1")) + end + + if length(protocol.params.sequences) != length(protocol.params.positions) + @warn "The MultiSequenceSystemMatrixProtocol has $(length(protocol.params.sequences)) sequences but you configured $(length(protocol.params.positions)) positions." + end + + if isnothing(getTemperatureSensor(protocol.scanner)) && protocol.params.coolDownTo < 100 + throw(ProtocolConfigurationError("The Protocol requires cooling down to $(protocol.params.coolDownTo) °C between measurements, but no temperature sensor is present in the scanner!)")) + end + + protocol.systemMeasState = SystemMatrixMeasState() + numBGMeas = protocol.params.numBGMeas + if !isnothing(protocol.params.bgSequence) + protocol.allSequences = [repeat([protocol.params.bgSequence], numBGMeas); protocol.params.sequences; repeat([protocol.params.bgSequence], numBGMeas)] + else + protocol.allSequences = protocol.params.sequences + end + numPos = length(protocol.allSequences) + + if numPos == 0 + throw(ProtocolConfigurationError("Protocol has no sequences configured!")) + end + measIsBGPos = zeros(Bool, numPos) + if !isnothing(protocol.params.bgSequence) + measIsBGPos[1:numBGMeas] .= true + measIsBGPos[end-numBGMeas+1:end] .= true + end + + framesPerPos = zeros(Int64, numPos) + posToIdx = zeros(Int64, numPos) + for (i, seq) in enumerate(protocol.allSequences) + framesPerPos[i] = acqNumFrames(seq) + end + numTotalFrames = sum(framesPerPos) + posToIdx[1] = 1 + posToIdx[2:end] = cumsum(framesPerPos)[1:end-1] .+ 1 + measIsBGFrame = zeros(Bool, numTotalFrames) + + protocol.systemMeasState.measIsBGPos = measIsBGPos + protocol.systemMeasState.posToIdx = posToIdx + protocol.systemMeasState.measIsBGFrame = measIsBGFrame + protocol.systemMeasState.currPos = 1 + protocol.systemMeasState.positions = protocol.params.positions + + + numRxChannels = length(rxChannels(protocol.params.sequences[1])) # kind of hacky, but actual rxChannels for RedPitaya are only set when setupRx is called + rxNumSamplingPoints = rxNumSamplesPerPeriod(protocol.params.sequences[1]) + numPeriods = acqNumPeriodsPerFrame(protocol.params.sequences[1]) + + #= Initialization of signals happens in execute, to handle restoring the protocol=# + + protocol.systemMeasState.currentSignal = zeros(Float32, rxNumSamplingPoints, numRxChannels, numPeriods, 1) + + + + if protocol.params.saveTemperatureData + sensor = getTemperatureSensor(protocol.scanner) + protocol.systemMeasState.temperatures = zeros(numChannels(sensor), numTotalFrames) + end + if protocol.params.controlTx + protocol.txCont = getDevice(protocol.scanner, TxDAQController) + else + protocol.txCont = nothing + end + return nothing +end + +function timeEstimate(protocol::MultiSequenceSystemMatrixProtocol) + if protocol.params.waitTime >= 5 + t = protocol.params.waitTime*length(protocol.params.sequences) * 1u"s" + est = timeFormat(t) + else + est = "Unknown" + end + #if !isnothing(protocol.params.sequence) + # params = protocol.params + # seq = params.sequence + # totalFrames = (params.fgFrames + params.bgFrames) * acqNumFrameAverages(seq) + # samplesPerFrame = rxNumSamplingPoints(seq) * acqNumAverages(seq) * acqNumPeriodsPerFrame(seq) + # totalTime = (samplesPerFrame * totalFrames) / (125e6 / (txBaseFrequency(seq) / rxSamplingRate(seq))) + # time = totalTime * 1u"s" + # est = string(time) + # @show est + #end + return est +end + +function enterExecute(protocol::MultiSequenceSystemMatrixProtocol) + protocol.paused = false + protocol.cancelled = false + protocol.finishAcknowledged = false + protocol.restored = false + protocol.systemMeasState.currPos = 1 +end + +function initMeasData(protocol::MultiSequenceSystemMatrixProtocol) + if isfile(protocol, "meta.toml") + message = """Found existing calibration file! \n + Should it be resumed?""" + if askConfirmation(protocol, message) + restore(protocol) + end + end + # Initialize Signals + if !protocol.restored + numRxChannels = length(rxChannels(protocol.params.sequences[1])) # kind of hacky, but actual rxChannels for RedPitaya are only set when setupRx is called + rxNumSamplingPoints = rxNumSamplesPerPeriod(protocol.params.sequences[1]) + numPeriods = acqNumPeriodsPerFrame(protocol.params.sequences[1]) + numTotalFrames = length(protocol.systemMeasState.measIsBGFrame) + rm(file(protocol, "signals.bin"), force=true) + signals = mmap!(protocol, "signals.bin", Float32, (rxNumSamplingPoints, numRxChannels, numPeriods, numTotalFrames)) + protocol.systemMeasState.signals = signals + protocol.systemMeasState.signals[:] .= 0.0 + end +end + + +function _execute(protocol::MultiSequenceSystemMatrixProtocol) + @debug "Measurement protocol started" + + initMeasData(protocol) + + finished = false + notifiedStop = false + while !finished + finished = performMeasurements(protocol) + + # paused + notifiedStop = false + while protocol.paused + handleEvents(protocol) + protocol.cancelled && throw(CancelException()) + if !notifiedStop + put!(protocol.biChannel, OperationSuccessfulEvent(PauseEvent())) + notifiedStop = true + end + if !protocol.paused + put!(protocol.biChannel, OperationSuccessfulEvent(ResumeEvent())) + end + sleep(0.05) + end + end + + + put!(protocol.biChannel, FinishedNotificationEvent()) + while !(protocol.finishAcknowledged) + handleEvents(protocol) + protocol.cancelled && throw(CancelException()) + end + + @info "Protocol finished." + close(protocol.biChannel) + @debug "Protocol channel closed after execution." +end + +function performMeasurements(protocol::MultiSequenceSystemMatrixProtocol) + finished = false + calib = protocol.systemMeasState + + while !finished + + wasRestored = protocol.restored + protocol.restored = false + + timeWaited = @elapsed begin + wait(calib.producer) + wait(calib.consumer) + end + diffTime = protocol.params.waitTime - timeWaited + if diffTime > 0.0 && !wasRestored && protocol.systemMeasState.currPos > 1 + @info "Wait $diffTime s for next measurement" + for _ in 1:diffTime/0.1 + handleEvents(protocol) + if protocol.paused || protocol.cancelled + break + end + sleep(0.1) + end + end + if protocol.paused || protocol.cancelled + enterPause(protocol) + finished = false + break + end + + tempSensor = getTemperatureSensor(protocol.scanner) + if !isnothing(tempSensor) + for i=1:60 + T = getTemperature(tempSensor,1) + if T > protocol.params.coolDownTo + @info "Coils still cooling down: $T °C (goal: $(protocol.params.coolDownTo) °C)" + sleep(1) + else + @info "Coil temperature of $T °C is below goal of $(protocol.params.coolDownTo) °C. Starting next measurement..." + break + end + if i==60 + error("Timeout while waiting for coils to cool down!") + end + end + end + + performMeasurement(protocol) + if protocol.systemMeasState.currPos > length(protocol.allSequences) + calib = protocol.systemMeasState + daq = getDAQ(protocol.scanner) + wait(calib.consumer) + wait(calib.producer) + stopTx(daq) + finished = true + end + end + + return finished +end + +function enterPause(protocol::MultiSequenceSystemMatrixProtocol) + calib = protocol.systemMeasState + wait(calib.consumer) + wait(calib.producer) +end + + +function performMeasurement(protocol::MultiSequenceSystemMatrixProtocol) + # Prepare + calib = protocol.systemMeasState + index = calib.currPos + @info "Measurement $index of $(length(protocol.allSequences))" + if index == 1 || calib.measIsBGPos[index] != calib.measIsBGPos[index-1] + text = if calib.measIsBGPos[index]; "background" else "foreground" end + if askChoices(protocol, "Press continue when $text measurement can be taken", ["Cancel", "Continue"]) == 1 + throw(CancelException()) + end + end + + daq = getDAQ(protocol.scanner) + + sequence = protocol.allSequences[index] + if protocol.params.controlTx + sequence = controlTx(protocol.txCont, sequence) + end + setup(daq, sequence) + + channel = Channel{channelType(daq)}(32) + calib.producer = @tspawnat protocol.scanner.generalParams.producerThreadID asyncProducer(channel, daq, sequence) + bind(channel, calib.producer) + calib.consumer = @tspawnat protocol.scanner.generalParams.consumerThreadID asyncConsumer(channel, protocol, index) + while !istaskdone(calib.producer) + handleEvents(protocol) + # Dont want to throw cancel here + sleep(0.05) + end + + # Increment measured positions + calib.currPos += 1 +end + +# function prepareDAQ(protocol::MultiSequenceSystemMatrixProtocol) +# calib = protocol.systemMeasState +# daq = getDAQ(protocol.scanner) + +# sequence = protocol.allSequences[calib.currPos] +# if protocol.params.controlTx +# if isnothing(protocol.contSequence) || protocol.restored || (calib.currPos == 1) +# sequence = controlTx(protocol.txCont, sequence) +# end +# #if isempty(protocol.systemMeasState.drivefield) +# # len = length(keys(protocol.contSequence.simpleChannel)) +# # drivefield = zeros(ComplexF64, len, len, size(calib.signals, 3), size(calib.signals, 4)) +# # calib.drivefield = mmap!(protocol, "observedField.bin", drivefield) +# # applied = zeros(ComplexF64, len, len, size(calib.signals, 3), size(calib.signals, 4)) +# # calib.applied = mmap!(protocol, "appliedFiled.bin", applied) +# #end +# #sequence = protocol.contSequence +# end +# setup(daq, sequence) +# if protocol.restored +# protocol.restored = false +# end +# end + +function asyncConsumer(channel::Channel, protocol::MultiSequenceSystemMatrixProtocol, index) + calib = protocol.systemMeasState + @info "readData" + daq = getDAQ(protocol.scanner) + numFrames = acqNumFrames(protocol.allSequences[index]) + startIdx = calib.posToIdx[index] + stopIdx = startIdx + numFrames - 1 + + # Prepare Buffer + deviceBuffer = DeviceBuffer[] + if protocol.params.saveTemperatureData + tempSensor = getTemperatureSensor(protocol.scanner) + push!(deviceBuffer, TemperatureBuffer(view(calib.temperatures, :, startIdx:stopIdx), tempSensor)) + end + + sinks = StorageBuffer[] + push!(sinks, FrameBuffer(1, view(calib.signals, :, :, :, startIdx:stopIdx))) + #sequence = protocol.params.sequence + #if protocol.params.controlTx + # sequence = protocol.contSequence + # push!(sinks, DriveFieldBuffer(1, view(calib.drivefield, :, :, :, startIdx:stopIdx), sequence)) + # push!(deviceBuffer, TxDAQControllerBuffer(1, view(calib.applied, :, :, :, startIdx:stopIdx), protocol.txCont)) + #end + + sequenceBuffer = AsyncBuffer(FrameSplitterBuffer(daq, sinks), daq) + asyncConsumer(channel, sequenceBuffer, deviceBuffer) + + calib.measIsBGFrame[startIdx:stopIdx] .= calib.measIsBGPos[index] + + calib.currentSignal = calib.signals[:, :, :, stopIdx:stopIdx] + + @info "store" + timeStore = @elapsed store(protocol, index) + @info "done after $timeStore" +end + +function store(protocol::MultiSequenceSystemMatrixProtocol, index) + + sysObj = protocol.systemMeasState + params = MPIFiles.toDict(sysObj.positions) + params["currPos"] = index + 1 # Safely stored up to and including index + #params["paused"] = protocol.paused + #params["currentSignal"] = sysObj.currentSignal + params["waitTime"] = protocol.params.waitTime + params["measIsBGPos"] = sysObj.measIsBGPos + params["posToIdx"] = sysObj.posToIdx + params["measIsBGFrame"] = sysObj.measIsBGFrame + #params["temperatures"] = vec(sysObj.temperatures) + params["sequences"] = toDict.(protocol.params.sequences) + if !isnothing(protocol.params.bgSequence) + params["bgSequence"] = toDict(protocol.params.bgSequence) + end + + filename = file(protocol, "meta.toml") + filename_backup = file(protocol, "meta.toml.backup") + if isfile(filename) + mv(filename, filename_backup, force=true) + end + + open(filename, "w") do f + TOML.print(f, params) + end + + Mmap.sync!(sysObj.signals) + #Mmap.sync!(sysObj.drivefield) + #Mmap.sync!(sysObj.applied) + rm(filename_backup, force=true) + return +end + +function restore(protocol::MultiSequenceSystemMatrixProtocol) + + sysObj = protocol.systemMeasState + params = MPIFiles.toDict(sysObj.positions) + if isfile(protocol, "meta.toml") + params = TOML.parsefile(file(protocol, "meta.toml")) + sysObj.currPos = params["currPos"] + protocol.paused = false + protocol.params.waitTime = params["waitTime"] + sysObj.measIsBGPos = params["measIsBGPos"] + sysObj.posToIdx = params["posToIdx"] + sysObj.measIsBGFrame = params["measIsBGFrame"] + #temp = params["temperatures"] + #if !isempty(temp) && (length(sysObj.temperatures) == length(temp)) + # sysObj.temperatures[:] .= temp + #end + + sysObj.positions = Positions(params) + seq = protocol.params.sequences + + storedSeq = sequenceFromDict.(params["sequences"]) + if storedSeq != seq + message = "Stored sequences do not match initialized sequence. Use stored sequence instead?" + if askChoices(protocol, message, ["Cancel", "Use"]) == 1 + throw(CancelException()) + end + seq = storedSeq + protocol.params.sequence # What is going on here? + end + + if !isnothing(protocol.params.bgSequence) + protocol.allSequences = [protocol.params.bgSequence; protocol.params.sequences; protocol.params.bgSequence] + else + protocol.allSequences = protocol.params.sequences + end + + # Drive Field + #if isfile(protocol, "observedField.bin") # sysObj.drivefield is still empty at point of (length(sysObj.drivefield) == length(drivefield)) + # sysObj.drivefield = mmap(protocol, "observedField.bin", ComplexF64) + #end + #if isfile(protocol, "appliedField.bin") + # sysObj.applied = mmap(protocol, "appliedField.bin", ComplexF64) + #end + + + sysObj.signals = mmap(protocol, "signals.bin", Float32) + + numTotalFrames = sum(acqNumFrames, protocol.allSequences) + numRxChannels = length(rxChannels(seq[1])) # kind of hacky, but actual rxChannels for RedPitaya are only set when setupRx is called + rxNumSamplingPoints = rxNumSamplesPerPeriod(seq[1]) + numPeriods = acqNumPeriodsPerFrame(seq[1]) + paramSize = (rxNumSamplingPoints, numRxChannels, numPeriods, numTotalFrames) + if size(sysObj.signals) != paramSize + throw(DimensionMismatch("Dimensions of stored signals $(size(sysObj.signals)) does not match initialized signals $paramSize")) + end + + protocol.restored = true + @info "Restored system matrix measurement" + end +end + +function cleanup(protocol::MultiSequenceSystemMatrixProtocol) + rm(dir(protocol), force = true, recursive = true) +end + +function pause(protocol::MultiSequenceSystemMatrixProtocol) + calib = protocol.systemMeasState + if calib.currPos <= length(calib.positions) + # OperationSuccessfulEvent is put when it actually is in the stop loop + protocol.paused = true + else + # paused has no concept once all measurements are done + put!(protocol.biChannel, OperationUnsuccessfulEvent(PauseEvent())) + end +end + +function resume(protocol::MultiSequenceSystemMatrixProtocol) + protocol.paused = false + protocol.restored = true + # OperationSuccessfulEvent is put when it actually leaves the stop loop +end + +function cancel(protocol::MultiSequenceSystemMatrixProtocol) + protocol.cancelled = true # Set cancel s.t. exception can be thrown when appropiate + protocol.paused = true # Set stop to reach a known/save state +end + +function handleEvent(protocol::MultiSequenceSystemMatrixProtocol, event::ProgressQueryEvent) + put!(protocol.biChannel, ProgressEvent(protocol.systemMeasState.currPos, length(protocol.allSequences), "Position", event)) +end + +function handleEvent(protocol::MultiSequenceSystemMatrixProtocol, event::DataQueryEvent) + data = nothing + if event.message == "CURR" + data = protocol.systemMeasState.currentSignal + elseif event.message == "BG" + sysObj = protocol.systemMeasState + index = sysObj.currPos + while index > 1 && !sysObj.measIsBGPos[index] + index = index - 1 + end + startIdx = sysObj.posToIdx[index] + data = copy(sysObj.signals[:, :, :, startIdx:startIdx]) + else + put!(protocol.biChannel, UnknownDataQueryEvent(event)) + end + put!(protocol.biChannel, DataAnswerEvent(data, event)) +end + +function handleEvent(protocol::MultiSequenceSystemMatrixProtocol, event::DatasetStoreStorageRequestEvent) + if false + # TODO this should be some sort of storage failure event + put!(protocol.biChannel, IllegaleStateEvent("Calibration measurement is not done yet. Cannot save!")) + else + store = event.datastore + scanner = protocol.scanner + mdf = event.mdf + data = protocol.systemMeasState.signals + positions = protocol.systemMeasState.positions + isBackgroundFrame = protocol.systemMeasState.measIsBGFrame + temperatures = nothing + if protocol.params.saveTemperatureData + temperatures = protocol.systemMeasState.temperatures + end + drivefield = nothing + if !isempty(protocol.systemMeasState.drivefield) + drivefield = protocol.systemMeasState.drivefield + end + applied = nothing + if !isempty(protocol.systemMeasState.applied) + applied = protocol.systemMeasState.applied + end + filename = saveasMDF(store, scanner, protocol.params.sequences[1], data, positions, isBackgroundFrame, mdf; storeAsSystemMatrix=protocol.params.saveInCalibFolder, temperatures=temperatures, drivefield=drivefield, applied=applied) + @show filename + put!(protocol.biChannel, StorageSuccessEvent(filename)) + end +end + +handleEvent(protocol::MultiSequenceSystemMatrixProtocol, event::FinishedAckEvent) = protocol.finishAcknowledged = true + +protocolInteractivity(protocol::MultiSequenceSystemMatrixProtocol) = Interactive() +protocolMDFStudyUse(protocol::MultiSequenceSystemMatrixProtocol) = UsingMDFStudy() diff --git a/src/Protocols/Protocol.jl b/src/Protocols/Protocol.jl index 5b4a17af..8f774eb6 100644 --- a/src/Protocols/Protocol.jl +++ b/src/Protocols/Protocol.jl @@ -7,8 +7,8 @@ export Protocol, ProtocolParams, name, description, scanner, params, abstract type ProtocolParams end -export ProtocolState, PS_UNDEFINED, PS_INIT, PS_RUNNING, PS_PAUSED, PS_FINISHED, PS_FAILED -@enum ProtocolState PS_UNDEFINED PS_INIT PS_RUNNING PS_PAUSED PS_FINISHED PS_FAILED +export ProtocolState, PS_UNDEFINED, PS_INIT, PS_RUNNING, PS_PAUSED, PS_FINISHED, PS_FAILED, PS_STOPPED +@enum ProtocolState PS_UNDEFINED PS_INIT PS_RUNNING PS_PAUSED PS_FINISHED PS_FAILED PS_STOPPED export @add_protocol_fields macro add_protocol_fields(paramType) @@ -95,10 +95,29 @@ abstract type ProtocolEvent end @mustimplement _init(protocol::Protocol) # Prepare necessary data structures, ex. buffer for samples + background meas @mustimplement _execute(protocol::Protocol) @mustimplement enterExecute(protocol::Protocol) # Prepare execution tracking flags, ex. execution done or cancelled flags -@mustimplement cleanup(protocol::Protocol) -@mustimplement stop(protocol::Protocol) -@mustimplement resume(protocol::Protocol) -@mustimplement cancel(protocol::Protocol) + + +### all Protocols must implement these events +function cleanup(protocol::Protocol) + # NOP +end + +function pause(protocol::Protocol) + put!(protocol.biChannel, OperationNotSupportedEvent(PauseEvent())) +end + +function resume(protocol::Protocol) + put!(protocol.biChannel, OperationNotSupportedEvent(ResumeEvent())) +end + +function stop(protocol::Protocol) + put!(protocol.biChannel, OperationNotSupportedEvent(StopEvent())) +end + +function cancel(protocol::Protocol) + put!(protocol.biChannel, OperationNotSupportedEvent(CancelEvent())) +end +### function init(protocol::Protocol) @debug "Initializing protocol $(name(protocol)) with inner initializer." @@ -276,13 +295,15 @@ function askChoices(protocol::Protocol, message::AbstractString, choices::Vector end end -handleEvent(protocol::Protocol, event::PauseEvent) = stop(protocol) -handleEvent(protocol::Protocol, event::StopEvent) = stop(protocol) # TODO Differentiate stop and pause +handleEvent(protocol::Protocol, event::PauseEvent) = pause(protocol) +handleEvent(protocol::Protocol, event::StopEvent) = stop(protocol) handleEvent(protocol::Protocol, event::ResumeEvent) = resume(protocol) handleEvent(protocol::Protocol, event::CancelEvent) = cancel(protocol) handleEvent(protocol::Protocol, event::ProtocolEvent) = put!(biChannel(protocol), UndefinedEvent(event)) #handleEvent(protocol::Protocol, event::InfoQueryEvent) = +# by default all event are not supported + function handleEvents(protocol::Protocol) while isready(protocol.biChannel) event = take!(protocol.biChannel) @@ -312,4 +333,5 @@ include("MPIForceProtocol.jl") include("RobotMPIMeasurementProtocol.jl") include("RobotBasedProtocol.jl") include("ContinousMeasurementProtocol.jl") +include("MultiSequenceSystemMatrixProtocol.jl") #include("TransferFunctionProtocol.jl") \ No newline at end of file diff --git a/src/Protocols/RobotBasedMagneticFieldStaticProtocol.jl b/src/Protocols/RobotBasedMagneticFieldStaticProtocol.jl index 091a28ba..2aadecb3 100644 --- a/src/Protocols/RobotBasedMagneticFieldStaticProtocol.jl +++ b/src/Protocols/RobotBasedMagneticFieldStaticProtocol.jl @@ -1,5 +1,10 @@ export RobotBasedMagneticFieldStaticProtocolParams, RobotBasedMagneticFieldStaticProtocol, measurement +""" +Parameter for the `RobotBasedMagneticFieldStaticProtocol` + +$FIELDS +""" Base.@kwdef mutable struct RobotBasedMagneticFieldStaticProtocolParams <: RobotBasedProtocolParams sequence::Union{Sequence, Nothing} = nothing positions::Union{GridPositions, Nothing} = nothing @@ -169,13 +174,13 @@ function performMeasurement(protocol::RobotBasedMagneticFieldStaticProtocol) addMeasuredPosition(protocol.measurement, nextPosition(protocol).data, field=field_) # fieldError=fieldError_, fieldFrequency=fieldFrequency_, timestamp=timestamp_, temperature=temperature_) end -function stop(protocol::RobotBasedMagneticFieldStaticProtocol) +function pause(protocol::RobotBasedMagneticFieldStaticProtocol) if protocol.currPos <= length(protocol.params.positions) # OperationSuccessfulEvent is put when it actually is in the stop loop protocol.stopped = true else # Stopped has no concept once all measurements are done - put!(protocol.biChannel, OperationUnsuccessfulEvent(StopEvent())) + put!(protocol.biChannel, OperationUnsuccessfulEvent(PauseEvent())) end end diff --git a/src/Protocols/RobotBasedProtocol.jl b/src/Protocols/RobotBasedProtocol.jl index a04aeb38..b2e8e4df 100644 --- a/src/Protocols/RobotBasedProtocol.jl +++ b/src/Protocols/RobotBasedProtocol.jl @@ -31,7 +31,7 @@ function _execute(protocol::RobotBasedProtocol) handleEvents(protocol) protocol.cancelled && throw(CancelException()) if !notifiedStop - put!(protocol.biChannel, OperationSuccessfulEvent(StopEvent())) + put!(protocol.biChannel, OperationSuccessfulEvent(PauseEvent())) notifiedStop = true end if !protocol.stopped diff --git a/src/Protocols/RobotBasedSystemMatrixProtocol.jl b/src/Protocols/RobotBasedSystemMatrixProtocol.jl index b81905aa..f69a8af5 100644 --- a/src/Protocols/RobotBasedSystemMatrixProtocol.jl +++ b/src/Protocols/RobotBasedSystemMatrixProtocol.jl @@ -1,6 +1,8 @@ export RobotBasedSystemMatrixProtocol, RobotBasedSystemMatrixProtocolParams """ -Parameters for the RobotBasedSystemMatrixProtocol +Parameters for the `RobotBasedSystemMatrixProtocol` + +$FIELDS """ Base.@kwdef mutable struct RobotBasedSystemMatrixProtocolParams <: RobotBasedProtocolParams "Minimum wait time between robot movements" @@ -82,11 +84,11 @@ function SystemMatrixMeasState() RegularGridPositions([1,1,1],[0.0,0.0,0.0],[0.0,0.0,0.0]), 1, Array{Float32,4}(undef,0,0,0,0), Array{Float32,4}(undef,0,0,0,0), Vector{Bool}(undef,0), Vector{Int64}(undef,0), Vector{Bool}(undef,0), - Matrix{Float64}(undef,0,0), Array{ComplexF64,4}(undef,0,0,0,0), Array{ComplexF64,4}(undef,0,0,0,0)) + Matrix{Float32}(undef,0,0), Array{ComplexF64,4}(undef,0,0,0,0), Array{ComplexF64,4}(undef,0,0,0,0)) end function requiredDevices(protocol::RobotBasedSystemMatrixProtocol) - result = [AbstractDAQ, Robot, SurveillanceUnit] + result = [AbstractDAQ, Robot] if protocol.params.controlTx push!(result, TxDAQController) end @@ -260,10 +262,10 @@ function prepareDAQ(protocol::RobotBasedSystemMatrixProtocol) protocol.restored = false end if isempty(protocol.systemMeasState.drivefield) - len = length(keys(protocol.contSequence.simpleChannel)) - drivefield = zeros(ComplexF64, len, len, size(calib.signals, 3), size(calib.signals, 4)) + bufferShape = controlMatrixShape(protocol.contSequence) + drivefield = zeros(ComplexF64, bufferShape[1], bufferShape[2], size(calib.signals, 3), size(calib.signals, 4)) calib.drivefield = mmap!(protocol, "observedField.bin", drivefield) - applied = zeros(ComplexF64, len, len, size(calib.signals, 3), size(calib.signals, 4)) + applied = zeros(ComplexF64, bufferShape[1], bufferShape[2], size(calib.signals, 3), size(calib.signals, 4)) calib.applied = mmap!(protocol, "appliedFiled.bin", applied) end sequence = protocol.contSequence @@ -271,7 +273,7 @@ function prepareDAQ(protocol::RobotBasedSystemMatrixProtocol) acqNumFrames(sequence, calib.measIsBGPos[calib.currPos] ? protocol.params.bgFrames : protocol.params.fgFrames) #acqNumFrameAverages(sequence, calib.measIsBGPos[calib.currPos] ? 1 : protocol.params.fgFrames) - acqNumFrameAverages(sequence, 1) + #acqNumFrameAverages(sequence, 1) setup(daq, sequence) end @@ -289,7 +291,9 @@ function postMovement(protocol::RobotBasedSystemMatrixProtocol) channelIdx = id.(vcat(acyclicElectricalTxChannels(protocol.params.sequence), periodicElectricalTxChannels(protocol.params.sequence))) amps = filter(amp -> in(channelId(amp), channelIdx), amps) end - enableACPower(su) + if !isnothing(su) + enableACPower(su) + end if tempControl != nothing disableControl(tempControl) end @@ -326,7 +330,9 @@ function postMovement(protocol::RobotBasedSystemMatrixProtocol) if tempControl != nothing enableControl(tempControl) end - disableACPower(su) + if !isnothing(su) + disableACPower(su) + end end end @@ -339,19 +345,30 @@ function asyncConsumer(channel::Channel, protocol::RobotBasedSystemMatrixProtoco stopIdx = startIdx + numFrames - 1 # Prepare Buffer - deviceBuffer = DeviceBuffer[] + deviceBuffer = StorageBuffer[] if protocol.params.saveTemperatureData tempSensor = getTemperatureSensor(protocol.scanner) push!(deviceBuffer, TemperatureBuffer(view(calib.temperatures, :, startIdx:stopIdx), tempSensor)) end sinks = StorageBuffer[] - push!(sinks, SimpleFrameBuffer(1, view(calib.signals, :, :, :, startIdx:stopIdx))) + buffer = FrameBuffer(1, view(calib.signals, :, :, :, startIdx:stopIdx)) + if acqNumFrameAverages(protocol.params.sequence) > 1 + buffer = AverageBuffer(buffer, protocol.params.sequence) + end + push!(sinks, buffer) + sequence = protocol.params.sequence if protocol.params.controlTx - sequence = protocol.contSequence - push!(sinks, DriveFieldBuffer(1, view(calib.drivefield, :, :, :, startIdx:stopIdx), sequence)) - push!(deviceBuffer, TxDAQControllerBuffer(1, view(calib.applied, :, :, :, startIdx:stopIdx), protocol.txCont)) + contSeq = protocol.contSequence + dfBuffer = DriveFieldBuffer(1, view(calib.drivefield, :, :, :, startIdx:stopIdx), contSeq) + #controlBuffer = TxDAQControllerBuffer(1, view(calib.applied, :, :, :, startIdx:stopIdx), protocol.txCont) + if acqNumFrameAverages(protocol.params.sequence) > 1 + dfBuffer = AverageBuffer(dfBuffer, rxNumSamplesPerPeriod(sequence), 2, acqNumPeriodsPerFrame(sequence), acqNumFrameAverages(sequence)) + #controlBuffer = AverageBuffer(controlBuffer, rxNumSamplesPerPeriod(sequence), 2, acqNumPeriodsPerFrame(sequence), acqNumFrameAverages(sequence)) + end + push!(sinks, dfBuffer) + #push!(deviceBuffer, controlBuffer) end sequenceBuffer = AsyncBuffer(FrameSplitterBuffer(daq, sinks), daq) @@ -473,14 +490,14 @@ function restore(protocol::RobotBasedSystemMatrixProtocol) end -function stop(protocol::RobotBasedSystemMatrixProtocol) +function pause(protocol::RobotBasedSystemMatrixProtocol) calib = protocol.systemMeasState if calib.currPos <= length(calib.positions) # OperationSuccessfulEvent is put when it actually is in the stop loop protocol.stopped = true else # Stopped has no concept once all measurements are done - put!(protocol.biChannel, OperationUnsuccessfulEvent(StopEvent())) + put!(protocol.biChannel, OperationUnsuccessfulEvent(PauseEvent())) end end @@ -502,7 +519,7 @@ end function handleEvent(protocol::RobotBasedSystemMatrixProtocol, event::DataQueryEvent) data = nothing - if event.message == "CURR" + if event.message == "SIGNAL" data = protocol.systemMeasState.currentSignal elseif event.message == "BG" sysObj = protocol.systemMeasState diff --git a/src/Protocols/RobotBasedTDesignFieldProtocol.jl b/src/Protocols/RobotBasedTDesignFieldProtocol.jl index f94a928c..8448b525 100644 --- a/src/Protocols/RobotBasedTDesignFieldProtocol.jl +++ b/src/Protocols/RobotBasedTDesignFieldProtocol.jl @@ -1,5 +1,10 @@ export RobotBasedTDesignFieldProtocolParams, RobotBasedTDesignFieldProtocol, measurement +""" +Parameters for the `RobotBasedTDesignFieldProtocol` + +$FIELDS +""" Base.@kwdef mutable struct RobotBasedTDesignFieldProtocolParams <: RobotBasedProtocolParams sequence::Union{Sequence, Nothing} = nothing radius::typeof(1.0u"mm") = 0.0u"mm" @@ -211,13 +216,13 @@ function finishMeasurement(protocol::RobotBasedTDesignFieldProtocol, gauss::Lake end end -function stop(protocol::RobotBasedTDesignFieldProtocol) +function pause(protocol::RobotBasedTDesignFieldProtocol) if protocol.currPos <= length(protocol.positions) # OperationSuccessfulEvent is put when it actually is in the stop loop protocol.stopped = true else # Stopped has no concept once all measurements are done - put!(protocol.biChannel, OperationUnsuccessfulEvent(StopEvent())) + put!(protocol.biChannel, OperationUnsuccessfulEvent(PauseEvent())) end end @@ -252,8 +257,5 @@ function handleEvent(protocol::RobotBasedTDesignFieldProtocol, event::FileStorag put!(protocol.biChannel, StorageSuccessEvent(filename)) end -function cleanup(protocol::RobotBasedTDesignFieldProtocol) - # NOP -end protocolInteractivity(protocol::RobotBasedTDesignFieldProtocol) = Interactive() diff --git a/src/Protocols/RobotMPIMeasurementProtocol.jl b/src/Protocols/RobotMPIMeasurementProtocol.jl index 9d082336..2d65c36b 100644 --- a/src/Protocols/RobotMPIMeasurementProtocol.jl +++ b/src/Protocols/RobotMPIMeasurementProtocol.jl @@ -1,6 +1,8 @@ export RobotMPIMeasurementProtocol, RobotMPIMeasurementProtocolParams """ -Parameters for the RobotMPIMeasurementProtocol +Parameters for the `RobotMPIMeasurementProtocol` + +$FIELDS """ Base.@kwdef mutable struct RobotMPIMeasurementProtocolParams <: ProtocolParams "Foreground position" @@ -236,19 +238,6 @@ function asyncMeasurement(protocol::RobotMPIMeasurementProtocol) return protocol.seqMeasState end - -function cleanup(protocol::RobotMPIMeasurementProtocol) - # NOP -end - -function stop(protocol::RobotMPIMeasurementProtocol) - put!(protocol.biChannel, OperationNotSupportedEvent(StopEvent())) -end - -function resume(protocol::RobotMPIMeasurementProtocol) - put!(protocol.biChannel, OperationNotSupportedEvent(ResumeEvent())) -end - function cancel(protocol::RobotMPIMeasurementProtocol) protocol.cancelled = true #put!(protocol.biChannel, OperationNotSupportedEvent(CancelEvent())) diff --git a/src/Protocols/Storage/ChainableBuffer.jl b/src/Protocols/Storage/ChainableBuffer.jl index 9a671756..4f9e34a8 100644 --- a/src/Protocols/Storage/ChainableBuffer.jl +++ b/src/Protocols/Storage/ChainableBuffer.jl @@ -5,7 +5,7 @@ mutable struct AverageBuffer{T} <: IntermediateBuffer where {T<:Number} end AverageBuffer(buffer::StorageBuffer, samples, channels, periods, avgFrames) = AverageBuffer{Float32}(buffer, zeros(Float32, samples, channels, periods, avgFrames), 1) AverageBuffer(buffer::StorageBuffer, sequence::Sequence) = AverageBuffer(buffer, rxNumSamplesPerPeriod(sequence), length(rxChannels(sequence)), acqNumPeriodsPerFrame(sequence), acqNumFrameAverages(sequence)) -function push!(avgBuffer::AverageBuffer{T}, frames::Array{T,4}) where {T<:Number} +function push!(avgBuffer::AverageBuffer{T}, frames::AbstractArray{T,4}) where {T<:Number} #setIndex - 1 = how many frames were written to the buffer # Compute how many frames there will be @@ -26,7 +26,7 @@ function push!(avgBuffer::AverageBuffer{T}, frames::Array{T,4}) where {T<:Number # Insert into buffer toFrames = fr + fit toAvg = avgBuffer.setIndex + fit - avgBuffer.buffer[:, :, :, avgBuffer.setIndex:toAvg] = frames[:, :, :, fr:toFrames] + avgBuffer.buffer[:, :, :, avgBuffer.setIndex:toAvg] = @view frames[:, :, :, fr:toFrames] avgBuffer.setIndex += length(avgBuffer.setIndex:toAvg) fr = toFrames + 1 @@ -49,31 +49,52 @@ sinks!(buffer::AverageBuffer, sinks::Vector{SinkBuffer}) = sinks!(buffer.target, abstract type MeasurementBuffer <: SequenceBuffer end # TODO Error handling? Throw own error or crash with index error -mutable struct SimpleFrameBuffer{A<: AbstractArray{Float32, 4}} <: MeasurementBuffer +mutable struct FrameBuffer{A<: AbstractArray{Float32, 4}} <: MeasurementBuffer nextFrame::Integer data::A end -function SimpleFrameBuffer(sequence::Sequence) +function FrameBuffer(sequence::Sequence) numFrames = acqNumFrames(sequence) rxNumSamplingPoints = rxNumSamplesPerPeriod(sequence) numPeriods = acqNumPeriodsPerFrame(sequence) numChannel = length(rxChannels(sequence)) + @debug "Creating FrameBuffer with size $rxNumSamplingPoints x $numChannel x $numPeriods x $numFrames" buffer = zeros(Float32, rxNumSamplingPoints, numChannel, numPeriods, numFrames) - return SimpleFrameBuffer(1, buffer) + return FrameBuffer(1, buffer) end -function insert!(buffer::SimpleFrameBuffer, from::Integer, frames::Array{Float32,4}) +function FrameBuffer(protocol::Protocol, file::String, sequence::Sequence) + numFrames = acqNumFrames(sequence) + rxNumSamplingPoints = rxNumSamplesPerPeriod(sequence) + numPeriods = acqNumPeriodsPerFrame(sequence) + numChannel = length(rxChannels(sequence)) + return FrameBuffer(protocol, file, Float32, (rxNumSamplingPoints, numChannel, numPeriods, numFrames)) +end +function FrameBuffer(protocol::Protocol, f::String, args...) + @debug "Creating memory-mapped FrameBuffer with size $(args[2])" + rm(file(protocol, f), force=true) + mapped = mmap!(protocol, f, args...) + return FrameBuffer(1, mapped) +end + +function insert!(op, buffer::FrameBuffer, from::Integer, frames::AbstractArray{Float32, 4}) + to = from + size(frames, 4) - 1 + frames = op(frames, view(buffer.data, :, :, :, from:to)) + insert!(buffer, from, frames) +end +function insert!(buffer::FrameBuffer, from::Integer, frames::AbstractArray{Float32,4}) to = from + size(frames, 4) - 1 buffer.data[:, :, :, from:to] = frames + buffer.nextFrame = to return to end -function push!(buffer::SimpleFrameBuffer, frames::Array{Float32,4}) +function push!(buffer::FrameBuffer, frames::AbstractArray{Float32,4}) from = buffer.nextFrame to = insert!(buffer, from, frames) buffer.nextFrame = to + 1 return (start = from, stop = to) end -read(buffer::SimpleFrameBuffer) = buffer.data -index(buffer::SimpleFrameBuffer) = buffer.nextFrame +read(buffer::FrameBuffer) = buffer.data +index(buffer::FrameBuffer) = buffer.nextFrame abstract type FieldBuffer <: SequenceBuffer end mutable struct DriveFieldBuffer{A <: AbstractArray{ComplexF64, 4}} <: FieldBuffer @@ -81,10 +102,16 @@ mutable struct DriveFieldBuffer{A <: AbstractArray{ComplexF64, 4}} <: FieldBuffe data::A cont::ControlSequence end +function insert!(op, buffer::DriveFieldBuffer, from::Integer, frames::AbstractArray{ComplexF64, 4}) + to = from + size(frames, 4) - 1 + frames = op(frames, view(buffer.data, :, :, :, from:to)) + insert!(buffer, from, frames) +end function insert!(buffer::DriveFieldBuffer, from::Integer, frames::Array{ComplexF64,4}) - # TODO duplicate to SimpleFrameBuffer + # TODO duplicate to FrameBuffer to = from + size(frames, 4) - 1 buffer.data[:, :, :, from:to] = frames + buffer.nextFrame = to return to end function push!(buffer::DriveFieldBuffer, frames::Array{ComplexF64,4}) @@ -93,6 +120,8 @@ function push!(buffer::DriveFieldBuffer, frames::Array{ComplexF64,4}) buffer.nextFrame = to + 1 return (start = from, stop = to) end +insert!(op, buffer::DriveFieldBuffer, from::Integer, frames::Array{Float32,4}) = insert!(op, buffer, from, calcFieldsFromRef(buffer.cont, frames)) +insert!(buffer::DriveFieldBuffer, from::Integer, frames::Array{Float32,4}) = insert!(buffer, from, calcFieldsFromRef(buffer.cont, frames)) push!(buffer::DriveFieldBuffer, frames::Array{Float32,4}) = push!(buffer, calcFieldsFromRef(buffer.cont, frames)) read(buffer::DriveFieldBuffer) = buffer.data index(buffer::DriveFieldBuffer) = buffer.nextFrame @@ -120,6 +149,44 @@ function push!(buffer::FrameSplitterBuffer, frames) end return result end +function insert!(buffer::FrameSplitterBuffer, from, frames) + uMeas, uRef = retrieveMeasAndRef!(frames, buffer.daq) + result = nothing + if !isnothing(uMeas) + for buf in buffer.targets + measSinks = length(sinks(buf, MeasurementBuffer)) + fieldSinks = length(sinks(buf, DriveFieldBuffer)) + if measSinks > 0 && fieldSinks == 0 + # Return latest measurement result + result = insert!(buf, from, uMeas) + elseif measSinks == 0 && fieldSinks > 0 + insert!(buf, from, uRef) + else + @warn "Unexpected sink combination $(typeof.(sinks(buf)))" + end + end + end + return result +end +function insert!(op, buffer::FrameSplitterBuffer, from, frames) + uMeas, uRef = retrieveMeasAndRef!(frames, buffer.daq) + result = nothing + if !isnothing(uMeas) + for buf in buffer.targets + measSinks = length(sinks(buf, MeasurementBuffer)) + fieldSinks = length(sinks(buf, DriveFieldBuffer)) + if measSinks > 0 && fieldSinks == 0 + # Return latest measurement result + result = insert!(op, buf, from, uMeas) + elseif measSinks == 0 && fieldSinks > 0 + insert!(op, buf, from, uRef) + else + @warn "Unexpected sink combination $(typeof.(sinks(buf)))" + end + end + end + return result +end function sinks!(buffer::FrameSplitterBuffer, sinks::Vector{SinkBuffer}) for buf in buffer.targets sinks!(buf, sinks) @@ -148,9 +215,8 @@ end function TxDAQControllerBuffer(tx::TxDAQController, sequence::ControlSequence) numFrames = acqNumFrames(sequence.targetSequence) numPeriods = acqNumPeriodsPerFrame(sequence.targetSequence) - # TODO function for length(keys(simpleChannel)) - len = length(keys(sequence.simpleChannel)) - buffer = zeros(ComplexF64, len, len, numPeriods, numFrames) + bufferShape = controlMatrixShape(sequence) + buffer = zeros(ComplexF64, bufferShape[1], bufferShape[2], numPeriods, numFrames) return TxDAQControllerBuffer(1, buffer, tx) end update!(buffer::TxDAQControllerBuffer, start, stop) = insert!(buffer, calcControlMatrix(buffer.tx.cont), start, stop) diff --git a/src/Protocols/Storage/MDF.jl b/src/Protocols/Storage/MDF.jl index 3aaf380c..314b68ac 100644 --- a/src/Protocols/Storage/MDF.jl +++ b/src/Protocols/Storage/MDF.jl @@ -58,7 +58,7 @@ end function MPIFiles.saveasMDF(store::DatasetStore, scanner::MPIScanner, sequence::Sequence, data::Array{Float32,4}, - positions::Positions, isBackgroundFrame::Vector{Bool}, mdf::MDFv2InMemory; storeAsSystemMatrix::Bool = false, deltaSampleSize::Union{Vector{typeof(1.0u"m")}, Nothing} = nothing, temperatures::Union{Array{Float32}, Nothing}=nothing, drivefield::Union{Array{ComplexF64}, Nothing}=nothing, applied::Union{Array{ComplexF64}, Nothing}=nothing) + positions::Union{Positions, AbstractArray}, isBackgroundFrame::Vector{Bool}, mdf::MDFv2InMemory; storeAsSystemMatrix::Bool = false, deltaSampleSize::Union{Vector{typeof(1.0u"m")}, Nothing} = nothing, temperatures::Union{Array{Float32}, Nothing}=nothing, drivefield::Union{Array{ComplexF64}, Nothing}=nothing, applied::Union{Array{ComplexF64}, Nothing}=nothing) if storeAsSystemMatrix study = MPIFiles.getCalibStudy(store) @@ -135,6 +135,30 @@ function fillMDFCalibration(mdf::MDFv2InMemory, positions::GridPositions; deltaS return end +function fillMDFCalibration(mdf::MDFv2InMemory, offsetFields::Union{AbstractArray, Nothing}; deltaSampleSize::Union{Vector{typeof(1.0u"m")}, Nothing} = nothing) + + # /calibration/ subgroup + + if !isnothing(deltaSampleSize) + deltaSampleSize = Float64.(ustrip.(uconvert.(Unitful.m, deltaSampleSize))) : nothing + end + + method = "hybrid" + order = "xyz" + + calibsize = size(offsetFields)[1:end-1] + offsetFields = reshape(offsetFields, prod(calibsize), :) + + mdf.calibration = MDFv2Calibration(; + deltaSampleSize = deltaSampleSize, + method = method, + offsetFields = permutedims(offsetFields), # Switch dimensions since the MDF specifies Ox3 but Julia is column major + order = order, + size = collect(calibsize) + ) + + return +end fillMDFCalibration(mdf::MDFv2InMemory, positions::Positions; deltaSampleSize::Union{Vector{typeof(1.0u"m")}, Nothing} = nothing) = @warn "Storing positions of type $(typeof(positions)) in MDF is not implemented" @@ -275,9 +299,9 @@ function fillMDFAcquisition(mdf::MDFv2InMemory, scanner::MPIScanner, sequence::S #acqGradient(mdf, acqGradient(sequence)) # TODO: Impelent in Sequence #MPIFiles.acqGradient(mdf, ustrip.(u"T/m", scannerGradient(scanner))) - MPIFiles.acqNumAverages(mdf, acqNumAverages(sequence)) + MPIFiles.acqNumAverages(mdf, acqNumAverages(sequence)*acqNumFrameAverages(sequence)) MPIFiles.acqNumFrames(mdf, length(measIsBackgroundFrame(mdf))) # important since we might have added BG frames - MPIFiles.acqNumPeriodsPerFrame(mdf, acqNumPeriodsPerFrame(sequence)) + MPIFiles.acqNumPeriodsPerFrame(mdf, size(measData(mdf), 3)) # Hacky to support cases where periods of data != periods of sequence offsetField_ = acqOffsetField(sequence) MPIFiles.acqOffsetField(mdf, isnothing(offsetField_) || !all(x-> x isa Unitful.MagneticFlux, offsetField_) ? nothing : ustrip.(u"T", offsetField_)) MPIFiles.acqStartTime(mdf, Dates.unix2datetime(time())) #seqCont.startTime) # TODO as parameter, start time from protocol @@ -285,7 +309,9 @@ function fillMDFAcquisition(mdf::MDFv2InMemory, scanner::MPIScanner, sequence::S # /acquisition/drivefield/ subgroup MPIFiles.dfBaseFrequency(mdf, ustrip(u"Hz", dfBaseFrequency(sequence))) MPIFiles.dfCycle(mdf, ustrip(u"s", dfCycle(sequence))) - MPIFiles.dfDivider(mdf, dfDivider(sequence)) + dividers = dfDivider(sequence) + dividers[.!isinteger.(dividers)] .= 0 + MPIFiles.dfDivider(mdf, Int.(dividers)) MPIFiles.dfNumChannels(mdf, dfNumChannels(sequence)) MPIFiles.dfPhase(mdf, ustrip.(u"rad", dfPhase(sequence))) MPIFiles.dfStrength(mdf, ustrip.(u"T", dfStrength(sequence))) @@ -301,13 +327,12 @@ function fillMDFAcquisition(mdf::MDFv2InMemory, scanner::MPIScanner, sequence::S MPIFiles.rxUnit(mdf, "V") # transferFunction - if hasTransferFunction(scanner) - numFreq = div(numSamplingPoints_,2)+1 - freq = collect(0:(numFreq-1))./(numFreq-1).*ustrip(u"Hz", rxBandwidth(sequence)) - tf_ = TransferFunction(scanner) - tf = tf_(freq,1:numRxChannels_) - MPIFiles.rxTransferFunction(mdf, tf) - MPIFiles.rxInductionFactor(mdf, tf_.inductionFactor) + if any(hasTransferFunction.([getDAQ(scanner)], id.(rxChannels(sequence)))) + tfs = transferFunction.([getDAQ(scanner)], id.(rxChannels(sequence))) + tf = reduce((x,y)->MPIFiles.combine(x,y,interpolate=true), tfs) + sampledTF = ustrip.(sampleTF(tf, mdf)) + MPIFiles.rxTransferFunction(mdf, sampledTF) + MPIFiles.rxInductionFactor(mdf, tf.inductionFactor) end end diff --git a/src/Protocols/Storage/ProducerConsumer.jl b/src/Protocols/Storage/ProducerConsumer.jl index d0c6ea3c..df129526 100644 --- a/src/Protocols/Storage/ProducerConsumer.jl +++ b/src/Protocols/Storage/ProducerConsumer.jl @@ -1,17 +1,15 @@ -SequenceMeasState(x, sequence::ControlSequence, sequenceBuffer::Nothing = nothing) = SequenceMeasState(x, sequence, StorageBuffer[]) -function SequenceMeasState(x, sequence::ControlSequence, sequenceBuffer::Vector{StorageBuffer}) +SequenceMeasState(daq::RedPitayaDAQ, sequence::ControlSequence, sequenceBuffer::Nothing = nothing) = SequenceMeasState(daq, sequence, StorageBuffer[]) +function SequenceMeasState(daq::RedPitayaDAQ, sequence::ControlSequence, sequenceBuffer::Vector{StorageBuffer}) numFrames = acqNumFrames(sequence.targetSequence) numPeriods = acqNumPeriodsPerFrame(sequence.targetSequence) - # TODO function for length(keys(simpleChannel)) - len = length(keys(sequence.simpleChannel)) - buffer = DriveFieldBuffer(1, zeros(ComplexF64, len, len, numPeriods, numFrames), sequence) - avgFrames = acqNumFrameAverages(sequence.targetSequence) - if avgFrames > 1 - samples = rxNumSamplesPerPeriod(sequence.targetSequence) - periods = acqNumPeriodsPerFrame(sequence.targetSequence) - buffer = AverageBuffer(buffer, samples, len, periods, avgFrames) - end - return SequenceMeasState(x, sequence.targetSequence, push!(sequenceBuffer, buffer)) + bufferShape = controlMatrixShape(sequence) + buffer = DriveFieldBuffer(1, zeros(ComplexF64, bufferShape[1], bufferShape[2], numPeriods, numFrames), sequence) + frameAvgs = acqNumFrameAverages(sequence.targetSequence) + if frameAvgs > 1 + samples = rxNumSamplesPerPeriod(sequence.targetSequence) + buffer = AverageBuffer(buffer, samples, length(daq.refChanIDs), numPeriods, frameAvgs) + end + return SequenceMeasState(daq, sequence.targetSequence, push!(sequenceBuffer, buffer)) end SequenceMeasState(protocol::Protocol, x, sequenceBuffer::Union{Nothing, Vector{StorageBuffer}} = nothing) = SequenceMeasState(getDAQ(scanner(protocol)), x, sequenceBuffer) function SequenceMeasState(daq::RedPitayaDAQ, sequence::Sequence, sequenceBuffer::Union{Nothing, Vector{StorageBuffer}} = nothing) @@ -19,7 +17,7 @@ function SequenceMeasState(daq::RedPitayaDAQ, sequence::Sequence, sequenceBuffer # Prepare buffering structures @debug "Allocating buffer for $numFrames frames" - buffer = SimpleFrameBuffer(sequence) + buffer = FrameBuffer(sequence) if acqNumFrameAverages(sequence) > 1 buffer = AverageBuffer(buffer, sequence) end @@ -115,7 +113,7 @@ function asyncProducer(channel::Channel, protocol::Protocol, sequence::Sequence) end asyncConsumer(measState::SequenceMeasState) = asyncConsumer(measState.channel, measState.sequenceBuffer, measState.deviceBuffers) -function asyncConsumer(channel::Channel, sequenceBuffer::StorageBuffer, deviceBuffers::Union{Vector{DeviceBuffer}, Nothing} = nothing) +function asyncConsumer(channel::Channel, sequenceBuffer::StorageBuffer, deviceBuffers = nothing) @debug "Consumer start" while isopen(channel) || isready(channel) while isready(channel) diff --git a/src/Scanner.jl b/src/Scanner.jl index e0f08d9d..8478a5c7 100644 --- a/src/Scanner.jl +++ b/src/Scanner.jl @@ -11,7 +11,7 @@ export MPIScanner, MPIScannerGeneral, scannerBoreSize, scannerFacility, Recursively find all concrete types of the given type. """ -function deepsubtypes(type::DataType) +function deepsubtypes(type::Type) subtypes_ = subtypes(type) allSubtypes = subtypes_ for subtype in subtypes_ @@ -26,7 +26,11 @@ end Retrieve the concrete type of a given supertype corresponding to a given string. """ -function getConcreteType(supertype_::DataType, type::String) +concreteTypesCache = Dict{String, Type}() +function getConcreteType(supertype_::Type, type::String) + if haskey(concreteTypesCache, type) + return concreteTypesCache[type] + end knownTypes = deepsubtypes(supertype_) foundImplementation = nothing for Implementation in knownTypes @@ -34,6 +38,7 @@ function getConcreteType(supertype_::DataType, type::String) foundImplementation = Implementation end end + push!(concreteTypesCache, type=>foundImplementation) return foundImplementation end @@ -52,21 +57,26 @@ function initiateDevices(configDir::AbstractString, devicesParams::Dict{String, # Get implementations for all devices in the specified order for deviceID in order params = nothing + configFile = nothing if haskey(devicesParams, deviceID) params = devicesParams[deviceID] + configFile = joinpath(configDir, "Scanner.toml") else params = deviceParams(configDir, deviceID) + configFile = joinpath(configDir, "Devices", deviceID*".toml") end if !isnothing(params) dependencies = filter(kv -> in(kv[1], get(params, "dependencies", String[])), devices) - device = Device(deviceID, params; dependencies = dependencies, robust = robust) + device = Device(deviceID, params; dependencies = dependencies, robust = robust, configFile = configFile) if !isOptional(device) && !isPresent(device) @error "The device with ID `$deviceID` should be present but isn't." end devices[deviceID] = device + else + throw(ScannerConfigurationError("The device ID `$deviceID` was not found in the configuration. Please check your configuration.")) end end @@ -74,7 +84,7 @@ function initiateDevices(configDir::AbstractString, devicesParams::Dict{String, return devices end -function Device(deviceID::String, deviceParams::Dict{String, Any}; dependencies::Dict{String, Device}, robust::Bool = false) +function Device(deviceID::String, deviceParams::Dict{String, Any}; dependencies::Dict{String, Device}, robust::Bool = false, configFile = nothing) device = nothing params = copy(deviceParams) @@ -103,7 +113,7 @@ function Device(deviceID::String, deviceParams::Dict{String, Any}; dependencies: end # Construct device - device = DeviceImpl(deviceID=deviceID, params=paramsInst, dependencies=dependencies_) # All other fields must have default values! + device = DeviceImpl(deviceID=deviceID, params=paramsInst, dependencies=dependencies_, configFile = configFile) # All other fields must have default values! # Initialize device checkDependencies(device) @@ -128,19 +138,24 @@ function getFittingDeviceParamsType(params::Dict{String, Any}, deviceType::Strin fittingDeviceParams = [] lastException = nothing + lastBacktrace = nothing for (i, paramType) in enumerate(tempDeviceParams) try tempParams = paramType(copy(params)) push!(fittingDeviceParams, tempParams) catch ex lastException = ex + lastBacktrace = Base.catch_backtrace() end end if length(fittingDeviceParams) == 1 return fittingDeviceParams[1] elseif length(fittingDeviceParams) == 0 && !isnothing(lastException) - throw(lastException) + Base.printstyled("ERROR: "; color=:red, bold=true) + Base.showerror(stdout, lastException) + Base.show_backtrace(stdout, lastBacktrace) + throw("The above error occured during device creation!") else return nothing end @@ -168,6 +183,8 @@ end General description of the scanner. Note: The fields correspond to the root section of an MDF file. + +$FIELDS """ Base.@kwdef struct MPIScannerGeneral "Bore size of the scanner." @@ -184,12 +201,8 @@ Base.@kwdef struct MPIScannerGeneral gradient::Union{typeof(1u"T/m"), Nothing} = nothing "Path of the dataset store." datasetStore::String = "" - "Default sequence of the scanner." - defaultSequence::String = "" "Default protocol of the scanner." defaultProtocol::String = "" - "Location of the scanner's transfer function." - transferFunction::String = "" "Thread ID of the producer thread." producerThreadID::Int32 = 2 "Thread ID of the consumer thread." @@ -204,12 +217,14 @@ end $(SIGNATURES) Basic description of a scanner. + +$(FIELDS) """ mutable struct MPIScanner "Name of the scanner" name::String - "Path to the used configuration directory." - configDir::String + "Path to the used configuration file." + configFile::String "General parameters of the scanner like its bore size or gradient." generalParams::MPIScannerGeneral "Device instances instantiated by the scanner from its configuration." @@ -226,12 +241,12 @@ mutable struct MPIScanner configDir = findConfigDir(name) params = getScannerParams(configDir) - @debug "Instantiating scanner `$name` from configuration file at `$filename`." + @info "Instantiating scanner `$name` from configuration file at `$(joinpath(configDir, "Scanner.toml"))`." generalParams = params_from_dict(MPIScannerGeneral, params["General"]) @assert generalParams.name == name "The folder name and the scanner name in the configuration do not match." devices = initiateDevices(configDir, params["Devices"], robust = robust) - scanner = new(name, configDir, generalParams, devices) + scanner = new(name, joinpath(configDir, "Scanner.toml"), generalParams, devices) return scanner end @@ -268,18 +283,15 @@ end "Name of the scanner" name(scanner::MPIScanner) = scanner.name +"Path to the used configuration file" +configFile(scanner::MPIScanner) = scanner.configFile + "Path to the used configuration directory." -configDir(scanner::MPIScanner) = scanner.configDir +configDir(scanner::MPIScanner) = dirname(scanner.configFile) "General parameters of the scanner like its bore size or gradient." generalParams(scanner::MPIScanner) = scanner.generalParams -"Location of the scanner's transfer function." -transferFunction(scanner::MPIScanner) = scanner.generalParams.transferFunction - -"Check, whether the scanner has a transfer function defined." -hasTransferFunction(scanner::MPIScanner) = transferFunction(scanner) != "" - """ $(SIGNATURES) @@ -401,15 +413,14 @@ init(scanner::MPIScanner, devices::Vector{<:Device}; kwargs...) = init(scanner, This does not initialize devices that themselves depend on the given devices. """ function init(scanner::MPIScanner, deviceIDs::Vector{String} = getDeviceIDs(scanner); kwargs...) - configDir = scanner.configDir - params = getScannerParams(configDir)["Devices"] + params = getScannerParams(configDir(scanner))["Devices"] devices = dependenciesDFS(deviceIDs, params) order = filter(x -> in(x, devices), params["initializationOrder"]) getDevices(close, scanner, order) - deviceDict = initiateDevices(configDir, params; kwargs..., order = order) + deviceDict = initiateDevices(configDir(scanner), params; kwargs..., order = order) for (id, device) in deviceDict scanner.devices[id] = device @@ -439,9 +450,6 @@ scannerGradient(scanner::MPIScanner) = scanner.generalParams.gradient "Path of the dataset store." scannerDatasetStore(scanner::MPIScanner) = scanner.generalParams.datasetStore -"Default sequence of the scanner." -defaultSequence(scanner::MPIScanner) = scanner.generalParams.defaultSequence - "Default protocol of the scanner." defaultProtocol(scanner::MPIScanner) = scanner.generalParams.defaultProtocol @@ -481,13 +489,21 @@ function Sequence(configdir::AbstractString, name::AbstractString) end return sequenceFromTOML(path) end - +function Sequences(configdir::AbstractString, name::AbstractString) + path = joinpath(configdir, "Sequences", name) + if !isdir(path) + error("Sequence-Directory $(path) not available!") + end + paths = sort(collect(readdir(path, join = true)), lt=natural) + return [sequenceFromTOML(p) for p in paths] +end """ $(SIGNATURES) Constructor for a sequence of `name` from the configuration directory specified for the scanner. """ Sequence(scanner::MPIScanner, name::AbstractString) = Sequence(configDir(scanner), name) +Sequences(scanner::MPIScanner, name::AbstractString) = Sequences(configDir(scanner), name) function Sequence(scanner::MPIScanner, dict::Dict) sequence = sequenceFromDict(dict) @@ -526,8 +542,6 @@ function MPIFiles.TransferFunction(configdir::AbstractString, name::AbstractStri return TransferFunction(path) end -MPIFiles.TransferFunction(scanner::MPIScanner) = TransferFunction(configDir(scanner),transferFunction(scanner)) - #### Protocol #### function getProtocolList(scanner::MPIScanner) path = joinpath(configDir(scanner), "Protocols/") diff --git a/src/Sequences/ContinuousElectricalChannel.jl b/src/Sequences/ContinuousElectricalChannel.jl index e7f0ad36..a1de9874 100644 --- a/src/Sequences/ContinuousElectricalChannel.jl +++ b/src/Sequences/ContinuousElectricalChannel.jl @@ -1,7 +1,7 @@ export ContinuousElectricalChannel "Electrical channel with a stepwise definition of values." -Base.@kwdef struct ContinuousElectricalChannel <: AcyclicElectricalTxChannel # TODO: Why is this named continuous? +Base.@kwdef mutable struct ContinuousElectricalChannel <: AcyclicElectricalTxChannel # TODO: Why is this named continuous? "ID corresponding to the channel configured in the scanner." id::AbstractString "Divider of sampling frequency." @@ -13,11 +13,13 @@ Base.@kwdef struct ContinuousElectricalChannel <: AcyclicElectricalTxChannel # T "Phase of the component for each period of the field." phase::typeof(1.0u"rad") "Offset of the channel. If defined in Tesla, the calibration configured in the scanner will be used." - offset::Union{typeof(1.0u"T"), typeof(1.0u"A")} = 0.0u"T" + offset::Union{typeof(1.0u"T"), typeof(1.0u"V"), typeof(1.0u"A")} = 0.0u"T" "Waveform of the component." waveform::Waveform = WAVEFORM_SINE end +unitIsTesla(chan::ContinuousElectricalChannel) = (dimension(chan.offset) == dimension(u"T")) && (dimension(chan.amplitude)==dimension(u"T")) + channeltype(::Type{<:ContinuousElectricalChannel}) = StepwiseTxChannel() function createFieldChannel(channelID::AbstractString, channelType::Type{ContinuousElectricalChannel}, channelDict::Dict{String, Any}) @@ -51,7 +53,17 @@ function createFieldChannel(channelID::AbstractString, channelType::Type{Continu end if haskey(channelDict, "phase") - phase = uparse.(channelDict["phase"]) + phaseDict = Dict("cosine"=>0.0u"rad", "cos"=>0.0u"rad","sine"=>pi/2u"rad", "sin"=>pi/2u"rad","-cosine"=>pi*u"rad", "-cos"=>pi*u"rad","-sine"=>-pi/2u"rad", "-sin"=>-pi/2u"rad") + + try + phase = uparse(channelDict["phase"]) + catch + if haskey(phaseDict, channelDict["phase"]) + phase = phaseDict[channelDict["phase"]] + else + error("The value $(channelDict["phase"]) for the phase could not be parsed. Use either a unitful value, or one of the predefined keywords ($(keys(phaseDict)))") + end + end else phase = 0.0u"rad" # Default phase end diff --git a/src/Sequences/MagneticField.jl b/src/Sequences/MagneticField.jl index 7991e328..8dfbb2a9 100644 --- a/src/Sequences/MagneticField.jl +++ b/src/Sequences/MagneticField.jl @@ -5,6 +5,8 @@ Description of a magnetic field. The field can either be electromagnetically or mechanically changed. The mechanical movement of e.g. an iron yoke would be defined within two channels, one electrical and one mechanical. + +$FIELDS """ Base.@kwdef struct MagneticField "Unique ID of the field description." @@ -50,6 +52,14 @@ function getindex(field::MagneticField, index::String) throw(KeyError(index)) end setindex!(field::MagneticField, txChannel::TxChannel, i::Integer) = channels(field)[i] = txChannel +function setindex!(field::MagneticField, txChannel::TxChannel, i::String) + for (index, channel) in enumerate(channels(field)) + if isequal(id(channel), i) + return setindex!(field, txChannel, index) + end + end + push!(field, txChannel) +end firstindex(field::MagneticField) = start_(field) lastindex(field::MagneticField) = length(field) keys(field::MagneticField) = map(id, field) @@ -86,7 +96,7 @@ safeEndInterval(field::MagneticField) = field.safeEndInterval export safeErrorInterval safeErrorInterval(field::MagneticField) = field.safeErrorInterval -control(field::MagneticField) = field.control +control(field::MagneticField) = field.control || field.decouple decouple(field::MagneticField) = field.decouple export electricalTxChannels @@ -101,6 +111,9 @@ periodicElectricalTxChannels(field::MagneticField) = channels(field, PeriodicEle export acyclicElectricalTxChannels acyclicElectricalTxChannels(field::MagneticField) = channels(field, AcyclicElectricalTxChannel) +export protocolTxChannels +protocolTxChannels(field::MagneticField) = channels(field, ProtocolTxChannel) + function toDict!(dict, field::MagneticField) for structField in [x for x in fieldnames(typeof(field)) if !in(x, [:id, :channels])] dict[String(structField)] = toDictValue(getproperty(field, structField)) diff --git a/src/Sequences/PeriodicElectricalChannel.jl b/src/Sequences/PeriodicElectricalChannel.jl index 6ae57dc2..cc9b1c59 100644 --- a/src/Sequences/PeriodicElectricalChannel.jl +++ b/src/Sequences/PeriodicElectricalChannel.jl @@ -2,11 +2,15 @@ export PeriodicElectricalComponent, SweepElectricalComponent, PeriodicElectricalChannel, ArbitraryElectricalComponent -"Component of an electrical channel with periodic base function." +""" +Component of an electrical channel with periodic base function. + +$FIELDS +""" Base.@kwdef mutable struct PeriodicElectricalComponent <: ElectricalComponent id::AbstractString "Divider of the component." - divider::Integer + divider::Rational{Int} "Amplitude (peak) of the component for each period of the field." amplitude::Union{Vector{typeof(1.0u"T")}, Vector{typeof(1.0u"V")}, Vector{typeof(1.0u"A")}} # Is it really the right choice to have the periods here? Or should it be moved to the MagneticField? "Phase of the component for each period of the field." @@ -15,9 +19,13 @@ Base.@kwdef mutable struct PeriodicElectricalComponent <: ElectricalComponent waveform::Waveform = WAVEFORM_SINE end -"Sweepable component of an electrical channel with periodic base function. +""" +Sweepable component of an electrical channel with periodic base function. Note: Does not allow for changes in phase since this would make the switch -between frequencies difficult." +between frequencies difficult. + +$FIELDS +""" Base.@kwdef mutable struct SweepElectricalComponent <: ElectricalComponent "Divider of the component." divider::Vector{Integer} @@ -27,25 +35,42 @@ Base.@kwdef mutable struct SweepElectricalComponent <: ElectricalComponent waveform::Waveform = WAVEFORM_SINE end +""" +Arbitrary waveform component of an electrical channel with periodic base function defined by a vector of $(RedPitayaDAQServer._awgBufferSize) values. + +$FIELDS +""" Base.@kwdef mutable struct ArbitraryElectricalComponent <: ElectricalComponent id::AbstractString "Divider of the component." - divider::Integer - values::Union{Vector{typeof(1.0u"T")}, Vector{typeof(1.0u"A")}, Vector{typeof(1.0u"V")}} + divider::Rational{Int} + "Amplitude scale of the base waveform for each period of the field" + amplitude::Union{Vector{typeof(1.0u"T")}, Vector{typeof(1.0u"A")}, Vector{typeof(1.0u"V")}} + "Phase of the component for each period of the field." + phase::Vector{typeof(1.0u"rad")} + "Values for the base waveform of the component, will be multiplied by `amplitude`" + values::Vector{Float64} end -"""Electrical channel based on based on periodic base functions. Only the -PeriodicElectricalChannel counts for the cycle length calculation""" -Base.@kwdef struct PeriodicElectricalChannel <: ElectricalTxChannel +""" +Electrical channel based on based on periodic base functions. Only the +PeriodicElectricalChannel counts for the cycle length calculation + +$FIELDS +""" +Base.@kwdef mutable struct PeriodicElectricalChannel <: ElectricalTxChannel "ID corresponding to the channel configured in the scanner." id::AbstractString "Components added for this channel." components::Vector{ElectricalComponent} "Offset of the channel. If defined in Tesla, the calibration configured in the scanner will be used." - offset::Union{typeof(1.0u"T"), typeof(1.0u"V")} = 0.0u"T" + offset::Union{typeof(1.0u"T"), typeof(1.0u"A"), typeof(1.0u"V")} = 0.0u"T" isDfChannel::Bool = true + dcEnabled::Bool = true end +unitIsTesla(chan::PeriodicElectricalChannel) = (dimension(offset(chan)) == dimension(u"T")) && all([dimension(amplitude(comp))==dimension(u"T") for comp in components(chan)]) + # Indexing Interface length(ch::PeriodicElectricalChannel) = length(components(ch)) function getindex(ch::PeriodicElectricalChannel, index::Integer) @@ -61,6 +86,14 @@ function getindex(ch::PeriodicElectricalChannel, index::String) throw(KeyError(index)) end setindex!(ch::PeriodicElectricalChannel, comp::ElectricalComponent, i::Integer) = components(ch)[i] = comp +function setindex!(ch::PeriodicElectricalChannel, comp::ElectricalComponent, i::String) + for (index, cmp) in enumerate(components(ch)) + if isequal(id(cmp), i) + return setindex!(field, comp, index) + end + end + push!(ch, comp) +end firstindex(ch::PeriodicElectricalChannel) = start_(ch) lastindex(ch::PeriodicElectricalChannel) = length(ch) keys(ch::PeriodicElectricalChannel) = map(id, ch) @@ -92,10 +125,12 @@ function createFieldChannel(channelID::AbstractString, ::Type{PeriodicElectrical tmp = uparse.(channelDict["offset"]) if eltype(tmp) <: Unitful.Current tmp = tmp .|> u"A" + elseif eltype(tmp) <: Unitful.Voltage + tmp = tmp .|> u"V" elseif eltype(tmp) <: Unitful.BField tmp = tmp .|> u"T" else - error("The value for an offset has to be either given as a current or in tesla. You supplied the type `$(eltype(tmp))`.") + error("The value for an offset has to be either given as a voltage, current or in tesla. You supplied the type `$(eltype(tmp))`.") end splattingDict[:offset] = tmp end @@ -103,6 +138,9 @@ function createFieldChannel(channelID::AbstractString, ::Type{PeriodicElectrical if haskey(channelDict, "isDfChannel") splattingDict[:isDfChannel] = channelDict["isDfChannel"] end + if haskey(channelDict, "dcEnabled") + splattingDict[:dcEnabled] = channelDict["dcEnabled"] + end components = Vector{ElectricalComponent}() componentsDict = [(k, v) for (k, v) in channelDict if v isa Dict] @@ -117,10 +155,9 @@ end function createChannelComponent(componentID::AbstractString, componentDict::Dict{String, Any}) if haskey(componentDict, "type") type = pop!(componentDict, "type") - knownComponents = MPIFiles.concreteSubtypes(ElectricalComponent) - index = findfirst(x -> x == type, string.(knownComponents)) - if !isnothing(index) - createChannelComponent(componentID, knownComponents[index], componentDict) + concreteType = getConcreteType(ElectricalComponent, type) + if !isnothing(concreteType) + createChannelComponent(componentID, concreteType, componentDict) else error("Component $componentID has an unknown channel type `$type`.") end @@ -130,8 +167,44 @@ function createChannelComponent(componentID::AbstractString, componentDict::Dict end function createChannelComponent(componentID::AbstractString, ::Type{PeriodicElectricalComponent}, componentDict::Dict{String, Any}) - divider = componentDict["divider"] + + divider, amplitude, phase = extractBasicComponentProperties(componentDict) + + if haskey(componentDict, "waveform") + waveform = toWaveform(componentDict["waveform"]) + else + waveform = WAVEFORM_SINE # Default to sine + end + + return PeriodicElectricalComponent(id=componentID, divider=divider, amplitude=amplitude, phase=phase, waveform=waveform) +end + +function createChannelComponent(componentID::AbstractString, ::Type{ArbitraryElectricalComponent}, componentDict::Dict{String, Any}) + divider, amplitude, phase = extractBasicComponentProperties(componentDict) + + if componentDict["values"] isa AbstractString # case 1: filename to waveform + filename = joinpath(homedir(), ".mpi", "Waveforms", componentDict["values"]) + try + values = reshape(h5read(filename, "/values"),:) + catch + throw(SequenceConfigurationError("Could not load the waveform $(componentDict["values"]), either the file does not exist at $filename or the file structure is wrong")) + end + else # case 2: vector of real numbers + values = componentDict["values"] + end + + if abs(values[1]-values[end])>0.01 # more than 1% of max value in 1 of 2^14 samples -> slew > 160 + @warn "The first and last value of your selected waveform are producing a jump of size $(abs(values[1]-values[end]))! Please check your waveform, if this is intended!" + end + + return ArbitraryElectricalComponent(id=componentID, divider=divider,amplitude=amplitude, phase=phase, values=values) +end + +function extractBasicComponentProperties(componentDict::Dict{String, Any}) + divider = eval(Meta.parse("$(componentDict["divider"])")) + if divider isa Dict; divider = divider["num"]//divider["den"] end + if !(divider isa Rational); divider=rationalize(divider) end amplitude = uparse.(componentDict["amplitude"]) if eltype(amplitude) <: Unitful.Current amplitude = amplitude .|> u"A" @@ -144,36 +217,28 @@ function createChannelComponent(componentID::AbstractString, ::Type{PeriodicElec end if haskey(componentDict, "phase") - phase = uparse.(componentDict["phase"]) + phaseDict = Dict("sine"=>0.0u"rad", "sin"=>0.0u"rad","cosine"=>pi/2u"rad", "cos"=>pi/2u"rad","-sine"=>pi*u"rad", "-sin"=>pi*u"rad","-cosine"=>-pi/2u"rad", "-cos"=>-pi/2u"rad") + phase = [] + for x in componentDict["phase"] + try + push!(phase, uparse.(x)) + catch + if haskey(phaseDict, x) + push!(phase, phaseDict[x]) + else + error("The value $x for the phase could not be parsed. Use either a unitful value, or one of the predefined keywords ($(keys(phaseDict)))") + end + end + end else phase = fill(0.0u"rad", length(divider)) # Default phase end - - if haskey(componentDict, "waveform") - waveform = toWaveform(componentDict["waveform"]) - else - waveform = WAVEFORM_SINE # Default to sine - end - return PeriodicElectricalComponent(id=componentID, divider=divider, amplitude=amplitude, phase=phase, waveform=waveform) -end - -function createChannelComponent(componentID::AbstractString, ::Type{ArbitraryElectricalComponent}, componentDict::Dict{String, Any}) - divider = componentDict["divider"] - values = uparse.(componentDict["values"]) - if eltype(values) <: Unitful.Current - values = values .|> u"A" - elseif eltype(values) <: Unitful.Voltage - values = values .|> u"V" - elseif eltype(values) <: Unitful.BField - values = values .|> u"T" - else - error("The values have to be either given as a current or in tesla. You supplied the type `$(eltype(values))`.") - end - return ArbitraryElectricalComponent(id=componentID, divider=divider, values=values) + return divider, amplitude, phase end -export offset +export offset, offset! offset(channel::PeriodicElectricalChannel) = channel.offset +offset!(channel::PeriodicElectricalChannel, offset::Union{typeof(1.0u"T"),typeof(1.0u"V")}) = channel.offset = offset export components components(channel::PeriodicElectricalChannel) = channel.components @@ -183,16 +248,27 @@ periodicElectricalComponents(channel::PeriodicElectricalChannel) = components(ch export arbitraryElectricalComponents arbitraryElectricalComponents(channel::PeriodicElectricalChannel) = components(channel, ArbitraryElectricalComponent) -cycleDuration(channel::PeriodicElectricalChannel, baseFrequency::typeof(1.0u"Hz")) = lcm([comp.divider for comp in components(channel)])/baseFrequency +cycleDuration(channel::PeriodicElectricalChannel, baseFrequency::typeof(1.0u"Hz")) = lcm([numerator(comp.divider) for comp in components(channel)])/baseFrequency isDfChannel(channel::PeriodicElectricalChannel) = channel.isDfChannel -export divider +# TODO/JA: check if this can automatically be implemented for all setters (and getters) +function waveform!(channel::PeriodicElectricalChannel, componentId::AbstractString, value) + index = findfirst(x -> id(x) == componentId, channel.components) + if !isnothing(index) + waveform!(channel.components[index], value) + else + throw(ArgumentError("Channel $(id(channel)) has no component with id $componentid")) + end +end + +export divider, divider! divider(component::ElectricalComponent, trigger::Integer=1) = length(component.divider) == 1 ? component.divider[1] : component.divider[trigger] +divider!(component::PeriodicElectricalComponent,value::Union{Rational,Integer}) = component.divider = value export amplitude, amplitude! -amplitude(component::PeriodicElectricalComponent; period::Integer=1) = component.amplitude[period] -function amplitude!(component::PeriodicElectricalComponent, value::Union{typeof(1.0u"T"),typeof(1.0u"V")}; period::Integer=1) +amplitude(component::Union{PeriodicElectricalComponent,ArbitraryElectricalComponent}; period::Integer=1) = component.amplitude[period] +function amplitude!(component::Union{PeriodicElectricalComponent,ArbitraryElectricalComponent}, value::Union{typeof(1.0u"T"),typeof(1.0u"V"),typeof(1.0u"A")}; period::Integer=1) if eltype(component.amplitude) != typeof(value) && length(component.amplitude) == 1 component.amplitude = typeof(value)[value] else @@ -200,31 +276,23 @@ function amplitude!(component::PeriodicElectricalComponent, value::Union{typeof( end end amplitude(component::SweepElectricalComponent; trigger::Integer=1) = component.amplitude[period] -amplitude(component::ArbitraryElectricalComponent) = maximum(abs.(component.values)) export phase, phase! -phase(component::PeriodicElectricalComponent, trigger::Integer=1) = component.phase[trigger] -phase!(component::PeriodicElectricalComponent, value::typeof(1.0u"rad"); period::Integer=1) = component.phase[period] = value +phase(component::Union{PeriodicElectricalComponent,ArbitraryElectricalComponent}, trigger::Integer=1) = component.phase[trigger] +phase!(component::Union{PeriodicElectricalComponent,ArbitraryElectricalComponent}, value::typeof(1.0u"rad"); period::Integer=1) = component.phase[period] = value phase(component::SweepElectricalComponent, trigger::Integer=1) = 0.0u"rad" -phase!(component::ArbitraryElectricalComponent, value::typeof(1.0u"rad"); period::Integer=1) = component.phase[period] = value -phase(component::ArbitraryElectricalComponent, trigger::Integer=1) = 0.0u"rad" #component.phase[period] + +export values, values! +values(component::ArbitraryElectricalComponent) = component.values +values!(component::ArbitraryElectricalComponent, values::Vector{Float64}) = component.values = values +scaledValues(component::ArbitraryElectricalComponent) = amplitude(component) .* circshift(values(component), round(Int,phase(component)/(2π*u"rad")*length(values(component)))) export waveform, waveform! waveform(component::ElectricalComponent) = component.waveform waveform!(component::ElectricalComponent, value) = component.waveform = value waveform(::ArbitraryElectricalComponent) = WAVEFORM_ARBITRARY - -values(component::ArbitraryElectricalComponent) = component.values - -function waveform!(channel::PeriodicElectricalChannel, componentId::AbstractString, value) - index = findfirst(x -> id(x) == componentId, channel.components) - if !isnothing(index) - waveform!(channel.components[index], value) - else - throw(ArgumentError("Channel $(id(channel)) has no component with id $componentid")) - end -end +waveform!(::ArbitraryElectricalComponent, value) = error("Can not change the waveform type of an ArbitraryElectricalComponent. Use values!() to change the waveform.") export id id(component::PeriodicElectricalComponent) = component.id diff --git a/src/Sequences/ProtocolOffsetElectricalChannel.jl b/src/Sequences/ProtocolOffsetElectricalChannel.jl new file mode 100644 index 00000000..cae9d58f --- /dev/null +++ b/src/Sequences/ProtocolOffsetElectricalChannel.jl @@ -0,0 +1,55 @@ +export ProtocolOffsetElectricalChannel + +""" +Offset channnel that gets translated into the correct representation by the `MPSMeasurementProtocol` + +$FIELDS +""" +Base.@kwdef mutable struct ProtocolOffsetElectricalChannel{T<:Union{typeof(1.0u"T"),typeof(1.0u"A"),typeof(1.0u"V")}} <: ProtocolTxChannel + "ID corresponding to the channel configured in the scanner." + id::AbstractString + offsetStart::T + offsetStop::T + numOffsets::Int64 +end + +channeltype(::Type{<:ProtocolOffsetElectricalChannel}) = StepwiseTxChannel() + +function createFieldChannel(channelID::AbstractString, ::Type{<:ProtocolOffsetElectricalChannel}, channelDict::Dict{String, Any}) + offsetStart = uparse.(channelDict["offsetStart"]) + if eltype(offsetStart) <: Unitful.Current + offsetStart = offsetStart .|> u"A" + elseif eltype(offsetStart) <: Unitful.Voltage + offsetStart = offsetStart .|> u"V" + elseif eltype(offsetStart) <: Unitful.BField + offsetStart = offsetStart .|> u"T" + else + error("The value for an offsetStart has to be either given as a current or in tesla. You supplied the type `$(eltype(tmp))`.") + end + + offsetStop = uparse.(channelDict["offsetStop"]) + if eltype(offsetStop) <: Unitful.Current + offsetStop = offsetStop .|> u"A" + elseif eltype(offsetStop) <: Unitful.Voltage + offsetStop = offsetStop .|> u"V" + elseif eltype(offsetStop) <: Unitful.BField + offsetStop = offsetStop .|> u"T" + else + error("The value for an offsetStop has to be either given as a current or in tesla. You supplied the type `$(eltype(tmp))`.") + end + + numOffsets = channelDict["numOffsets"] + + return ProtocolOffsetElectricalChannel(;id = channelID, offsetStart = offsetStart, offsetStop = offsetStop, numOffsets = numOffsets) +end + +values(channel::ProtocolOffsetElectricalChannel{T}) where T = collect(range(channel.offsetStart, channel.offsetStop, length = channel.numOffsets)) +values(channel::ProtocolOffsetElectricalChannel, isPositive::Bool) = filter(x-> signbit(x) != isPositive, values(channel)) + +function toDict!(dict, channel::ProtocolOffsetElectricalChannel) + dict["type"] = string(typeof(channel)) + for field in [x for x in fieldnames(typeof(channel)) if !in(x, [:id])] + dict[String(field)] = toDictValue(getproperty(channel, field)) + end + return dict +end \ No newline at end of file diff --git a/src/Sequences/Sequence.jl b/src/Sequences/Sequence.jl index 2711d4dc..5a262276 100644 --- a/src/Sequences/Sequence.jl +++ b/src/Sequences/Sequence.jl @@ -53,6 +53,14 @@ function getindex(seq::Sequence, index::String) throw(KeyError(index)) end setindex!(seq::Sequence, field::MagneticField, i::Integer) = fields(seq)[i] = field +function setindex!(seq::Sequence, field::MagneticField, i::String) + for (index, f) in enumerate(fields(seq)) + if id(f) == i + return setindex!(seq, field, index) + end + end + push!(seq, field) +end firstindex(seq::Sequence) = start_(seq) lastindex(seq::Sequence) = length(seq) keys(seq::Sequence) = map(id, seq) @@ -161,6 +169,12 @@ function fieldDictToFields(fieldsDict::Dict{String, Any}) if haskey(fieldDict, "safeErrorInterval") splattingDict[:safeErrorInterval] = uparse(fieldDict["safeErrorInterval"]) end + if haskey(splattingDict, :decouple) && splattingDict[:decouple] + if haskey(splattingDict, :control) && !splattingDict[:control] + throw(SequenceConfigurationError("A field that should be decoupled always needs to be controlled, please fix field \"$(fieldID)\"")) + end + splattingDict[:control] = true + end field = MagneticField(;splattingDict...) push!(fields, field) @@ -173,10 +187,9 @@ end function createFieldChannel(channelID::AbstractString, channelDict::Dict{String, Any}) if haskey(channelDict, "type") type = pop!(channelDict, "type") - knownChannels = MPIFiles.concreteSubtypes(TxChannel) - index = findfirst(x -> x == type, string.(knownChannels)) - if !isnothing(index) - createFieldChannel(channelID, knownChannels[index], channelDict) + concreteType = getConcreteType(TxChannel, type) + if !isnothing(concreteType) + createFieldChannel(channelID, concreteType, channelDict) else error("Channel $channelID has an unknown channel type `$type`.") end @@ -277,6 +290,17 @@ function phase!(channel::PeriodicElectricalChannel, componentId::AbstractString, end end +export divider! +function divider!(channel::PeriodicElectricalChannel, componentId::AbstractString, value::Union{Rational,Integer}) + index = findfirst(x -> id(x) == componentId, channel.components) + if !isnothing(index) + divider!(channel.components[index], value) + else + throw(ArgumentError("Channel $(id(channel)) has no component with id $componentid")) + end +end + + export acqGradient acqGradient(sequence::Sequence) = nothing # TODO: Implement @@ -285,7 +309,7 @@ function acqNumPeriodsPerFrame(sequence::Sequence) #TODO: We can't limit this to acyclic channels. What is the correct number of periods per frame with mechanical channels? if hasAcyclicElectricalTxChannels(sequence) channels = acyclicElectricalTxChannels(sequence) - samplesPerCycle = lcm(dfDivider(sequence)) + samplesPerCycle = dfSamplesPerCycle(sequence) numPeriods = [div(c.divider, samplesPerCycle) for c in channels ] if minimum(numPeriods) != maximum(numPeriods) @@ -303,7 +327,7 @@ function acqNumPeriodsPerPatch(sequence::Sequence) #TODO: We can't limit this to acyclic channels. What is the correct number of periods per frame with mechanical channels? if hasAcyclicElectricalTxChannels(sequence) channels = acyclicElectricalTxChannels(sequence) - samplesPerCycle = lcm(dfDivider(sequence)) + samplesPerCycle = dfSamplesPerCycle(sequence) stepsPerCycle = [ typeof(c) <: StepwiseElectricalChannel ? div(c.divider,length(c.values)*samplesPerCycle) : div(c.dividerSteps,samplesPerCycle) for c in channels ] @@ -323,13 +347,14 @@ export acqOffsetField function acqOffsetField(sequence::Sequence) # TODO: This is a hack for getting the required information for the MPSMeasurementProtocol. Can we find a generalized solution? if hasAcyclicElectricalTxChannels(sequence) - @warn "This is a hack for the MPSMeasurementProtocol. It might result in wrong MDF settings in other cases." - channels = acyclicElectricalTxChannels(sequence) - offsetChannel = first([channel for channel in channels if channel isa ContinuousElectricalChannel]) - values_ = MPIMeasurements.values(offsetChannel) - values3D = reshape([values_ fill(0.0u"T", length(values_)) fill(0.0u"T", length(values_))], (length(values_), 1, 3)) - - return values3D + #@warn "This is a hack for the MPSMeasurementProtocol. It might result in wrong MDF settings in other cases." + #channels = acyclicElectricalTxChannels(sequence) + #offsetChannel = first([channel for channel in channels if channel isa ContinuousElectricalChannel]) + #values_ = MPIMeasurements.values(offsetChannel) + #values3D = reshape([values_ fill(0.0u"T", length(values_)) fill(0.0u"T", length(values_))], (length(values_), 1, 3)) + # + #return values3D + return nothing else return nothing end @@ -341,8 +366,11 @@ dfBaseFrequency(sequence::Sequence) = baseFrequency(sequence) export txBaseFrequency txBaseFrequency(sequence::Sequence) = dfBaseFrequency(sequence) # Alias, since this might not only concern the drivefield +export dfSamplesPerCycle +dfSamplesPerCycle(sequence::Sequence) = lcm(numerator.(dfDivider(sequence))) + export dfCycle -dfCycle(sequence::Sequence) = lcm(dfDivider(sequence))/dfBaseFrequency(sequence) |> u"s" +dfCycle(sequence::Sequence) = dfSamplesPerCycle(sequence)/dfBaseFrequency(sequence) |> u"s" export txCycle txCycle(sequence::Sequence) = dfCycle(sequence) # Alias, since this might not only concern the drivefield @@ -351,7 +379,7 @@ export dfDivider function dfDivider(sequence::Sequence) # TODO: How do we integrate the mechanical channels and non-periodic channels and sweeps? channels = dfChannels(sequence) maxComponents = maximum([length(channel.components) for channel in channels]) - result = zeros(Int64, (dfNumChannels(sequence), maxComponents)) + result = ones(Rational{Int64}, (dfNumChannels(sequence), maxComponents)) # Otherwise the dfCycle is miscalculated in the case of a different amount of components for (channelIdx, channel) in enumerate(channels) for (componentIdx, component) in enumerate(channel.components) @@ -395,6 +423,11 @@ function dfStrength(sequence::Sequence) # TODO: How do we integrate the mechanic for (channelIdx, channel) in enumerate(channels) for (componentIdx, component) in enumerate(channel.components) for (periodIdx, strength) in enumerate(amplitude(component)) # TODO: What do we do if this is in volt? The conversion factor is with the scanner... Remove the volt version? + if strength isa Unitful.Voltage + # HACKY way to save Voltage DFs into the MDF -> use Tesla as mV (everything > 1 will be V) + @warn "Your DF is defined as a voltage! Will save $(uconvert(u"mV",strength)) as $(ustrip(u"mV", strength).*u"T")" + strength = ustrip(u"mV", strength).*u"T" + end result[periodIdx, channelIdx, componentIdx] = strength end end @@ -456,11 +489,12 @@ rxNumChannels(sequence::Sequence) = length(rxChannels(sequence)) export rxNumSamplingPoints function rxNumSamplingPoints(sequence::Sequence) result = upreferred(rxSamplingRate(sequence)*dfCycle(sequence)) - if !isinteger(result) - throw(ScannerConfigurationError("The selected combination of divider and decimation results in non-integer sampling points.")) + if !(result≈round(result)) + @debug "rxNumSamplingPoints" result≈round(result) rxSamplingRate(sequence) dfCycle(sequence) + throw(ScannerConfigurationError("The selected combination of divider ($(dfSamplesPerCycle(sequence))) and decimation ($(Int(baseFrequency(sequence)/rxSamplingRate(sequence)))) results in non-integer sampling points ($result).")) end - return Int64(result) + return round(Int64,result) end export rxNumSamplesPerPeriod @@ -471,7 +505,7 @@ rxChannels(sequence::Sequence) = rxChannels(sequence.acquisition) for T in [Sequence, GeneralSettings, AcquisitionSettings, MagneticField, TxChannel, ContinuousElectricalChannel, ContinuousMechanicalRotationChannel, ContinuousMechanicalTranslationChannel, PeriodicElectricalChannel, PeriodicElectricalComponent, SweepElectricalComponent, StepwiseElectricalChannel, - StepwiseMechanicalRotationChannel, StepwiseMechanicalTranslationChannel, ArbitraryElectricalComponent] + StepwiseMechanicalRotationChannel, StepwiseMechanicalTranslationChannel, ArbitraryElectricalComponent, ProtocolOffsetElectricalChannel] @eval begin @generated function ==(x::$T, y::$T) fieldEqualities = [:(x.$field == y.$field) for field in fieldnames($T)] diff --git a/src/Sequences/StepwiseElectricalChannel.jl b/src/Sequences/StepwiseElectricalChannel.jl index 556bee8e..83043811 100644 --- a/src/Sequences/StepwiseElectricalChannel.jl +++ b/src/Sequences/StepwiseElectricalChannel.jl @@ -1,7 +1,11 @@ export StepwiseElectricalChannel -"Electrical channel with a stepwise definition of values." -Base.@kwdef struct StepwiseElectricalChannel <: AcyclicElectricalTxChannel +""" +Electrical channel with a stepwise definition of values. + +$FIELDS +""" +Base.@kwdef mutable struct StepwiseElectricalChannel <: AcyclicElectricalTxChannel "ID corresponding to the channel configured in the scanner." id::AbstractString "Divider of the component." @@ -9,14 +13,16 @@ Base.@kwdef struct StepwiseElectricalChannel <: AcyclicElectricalTxChannel "Values corresponding to the individual steps." values::Union{Vector{typeof(1.0u"T")}, Vector{typeof(1.0u"A")}, Vector{typeof(1.0u"V")}} "TBD" - enable::Vector{Bool} + enable::Vector{Bool} = Bool[] end +unitIsTesla(chan::StepwiseElectricalChannel) = all(dimension(chan.values) .== [dimension(u"T")]) + channeltype(::Type{<:StepwiseElectricalChannel}) = StepwiseTxChannel() function createFieldChannel(channelID::AbstractString, channelType::Type{StepwiseElectricalChannel}, channelDict::Dict{String, Any}) divider = channelDict["divider"] - enable = haskey(channelDict, "enable") ? parsePossibleURange(channelDict["enable"]) : Bool[] + enable = haskey(channelDict, "enable") ? channelDict["enable"] : Bool[] values = parsePossibleURange(channelDict["values"]) if eltype(values) <: Unitful.Current values = values .|> u"A" @@ -36,6 +42,7 @@ function createFieldChannel(channelID::AbstractString, channelType::Type{Stepwis end values(channel::StepwiseElectricalChannel) = channel.values +values!(channel::StepwiseElectricalChannel, val::Union{Vector{typeof(1.0u"T")}, Vector{typeof(1.0u"A")}, Vector{typeof(1.0u"V")}}) = channel.values = val cycleDuration(channel::StepwiseElectricalChannel, baseFrequency::typeof(1.0u"Hz")) = upreferred(channel.divider/baseFrequency) stepsPerCycle(channel::StepwiseElectricalChannel) = length(channel.values) diff --git a/src/Sequences/TxChannel.jl b/src/Sequences/TxChannel.jl index 41958a0b..3966f4c5 100644 --- a/src/Sequences/TxChannel.jl +++ b/src/Sequences/TxChannel.jl @@ -8,11 +8,12 @@ abstract type TxMechanicalMovementType end struct RotationTxChannel <: TxMechanicalMovementType end struct TranslationTxChannel <: TxMechanicalMovementType end -export TxChannel, ElectricalTxChannel, AcyclicElectricalTxChannel, MechanicalTxChannel, ElectricalComponent +export TxChannel, ElectricalTxChannel, AcyclicElectricalTxChannel, MechanicalTxChannel, ProtocolTxChannel, ElectricalComponent abstract type TxChannel end abstract type ElectricalTxChannel <: TxChannel end abstract type AcyclicElectricalTxChannel <: ElectricalTxChannel end abstract type MechanicalTxChannel <: TxChannel end +abstract type ProtocolTxChannel <: TxChannel end abstract type ElectricalComponent end @@ -51,6 +52,9 @@ cycleDuration(::T, var) where T <: TxChannel = error("The method has not been im export id id(channel::TxChannel) = channel.id +export unitIsTesla +unitIsTesla(channel::ElectricalTxChannel) = false #default fallback + function toDict!(dict, channel::TxChannel) for field in [x for x in fieldnames((typeof(channel))) if x != :id] dict[String(field)] = toDictValue(getproperty(channel, field)) @@ -65,4 +69,5 @@ include("ContinuousElectricalChannel.jl") include("ContinuousMechanicalTranslationChannel.jl") include("StepwiseMechanicalTranslationChannel.jl") include("StepwiseMechanicalRotationChannel.jl") -include("ContinuousMechanicalRotationChannel.jl") \ No newline at end of file +include("ContinuousMechanicalRotationChannel.jl") +include("ProtocolOffsetElectricalChannel.jl") \ No newline at end of file diff --git a/src/Utils/Console/ConsoleProtocolHandler.jl b/src/Utils/Console/ConsoleProtocolHandler.jl index 88cc2a02..17dcb7a8 100644 --- a/src/Utils/Console/ConsoleProtocolHandler.jl +++ b/src/Utils/Console/ConsoleProtocolHandler.jl @@ -71,7 +71,7 @@ function initProtocol(cph::ConsoleProtocolHandler) end export startProtocol -function startProtocol(cph::ConsoleProtocolHandler) +function startProtocol(cph::ConsoleProtocolHandler; block=false) try @info "Execute protocol with name `$(name(cph.protocol))`" @@ -101,6 +101,17 @@ function startProtocol(cph::ConsoleProtocolHandler) cph.protocolState = PS_INIT @debug "Start event handler" cph.eventHandler = Timer(timer -> eventHandler(cph, timer), 0.0, interval=0.05) + if block + while (cph.protocolState != PS_FINISHED || isopen(cph.eventHandler)) + #@info cph.protocolState isopen(cph.eventHandler) + if cph.protocolState==PS_FAILED + @error "Protocol failed!" + return false + end + sleep(0.5) + end + #@info cph.protocolState isopen(cph.eventHandler) + end return true end catch e @@ -157,7 +168,7 @@ function eventHandler(cph::ConsoleProtocolHandler, timer::Timer) end catch ex - @error "The eventhandler catched an exception." + @error "The eventhandler caught an exception." confirmFinishedProtocol(cph) close(timer) #showError(ex) diff --git a/src/Utils/DictToStruct.jl b/src/Utils/DictToStruct.jl index cad537f1..063c762b 100644 --- a/src/Utils/DictToStruct.jl +++ b/src/Utils/DictToStruct.jl @@ -46,6 +46,16 @@ function parsePossibleURange(val::String) end parsePossibleURange(val::Vector{String}) = parsePossibleURange.(val) +function parse_into_tf(value::String) + if occursin(".h5", value) # case 1: filename to transfer function, the TF will be read the first time calibration() is called, (done in _init(), to prevent delays while using the device) + calibration_tf = value + else # case 2: single value, extended into transfer function with no frequency dependency + calibration_value = uparse(value) + calibration_tf = TransferFunction([0,Inf],ComplexF64[ustrip(calibration_value), ustrip(calibration_value)], units=[unit(calibration_value)]) + end + return calibration_tf +end + function dict_to_splatting(dict::Dict) splattingDict = Dict{Symbol, Any}() for (key, value) in dict diff --git a/src/Utils/MmapFiles.jl b/src/Utils/MmapFiles.jl index da83d9b3..58483e8b 100644 --- a/src/Utils/MmapFiles.jl +++ b/src/Utils/MmapFiles.jl @@ -3,23 +3,33 @@ dir(protocol::Protocol) = joinpath(dir(scanner(protocol)), name(protocol)) file(protocol::Protocol, file::String) = joinpath(mkpath(dir(protocol)), file) isfile(protocol::Protocol, name::String) = isfile(file(protocol, name)) -function mmap!(protocol::Protocol, f::String, array::Array{T,N}) where {T, N} - open(file(protocol, f), "w+") do io - write(io, N) - for s in size(array) - write(io, s) - end - write(io, array) +function mmap!(protocol::Protocol, f::String, array::Array{T,N}) where {T,N} + open(file(protocol, f), "w+") do io + write(io, N) + for s in size(array) + write(io, s) end - return mmap(protocol, f, T) + write(io, array) + end + return mmap(protocol, f, T) end -function mmap(protocol::Protocol, f::String, eltype::Type{T}) where {T} - io = open(file(protocol, f), "r+") - N = Base.read(io, Int64) - dims = [] - for i = 1:N - push!(dims, Base.read(io, Int64)) +function mmap!(protocol::Protocol, f::String, eltype::DataType, size::Union{Tuple{Vararg{Int64}}, Vector{Int64}}) + open(file(protocol, f), "w+") do io + write(io, length(size)) + for s in size + write(io, s) end - return Mmap.mmap(io, Array{eltype,N}, tuple(dims...)) + end + return mmap(protocol, f, eltype) +end + +function mmap(protocol::Protocol, f::String, eltype::Type{T}) where {T} + io = open(file(protocol, f), "r+") + N = Base.read(io, Int64) + dims = [] + for i = 1:N + push!(dims, Base.read(io, Int64)) + end + return Mmap.mmap(io, Array{eltype,N}, tuple(dims...)) end \ No newline at end of file diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl index 6800f982..425dd4d4 100644 --- a/src/Utils/Utils.jl +++ b/src/Utils/Utils.jl @@ -14,4 +14,22 @@ end Base.convert(::Type{IPAddr}, str::AbstractString) = parse(IPAddr, str) # TODO: Remove this type piracy -Base.convert(::Type{ClusterTriggerSetup}, str::AbstractString) = stringToEnum(str, ClusterTriggerSetup) \ No newline at end of file +Base.convert(::Type{ClusterTriggerSetup}, str::AbstractString) = stringToEnum(str, ClusterTriggerSetup) + +# should be available in Unitful but isnt https://github.com/PainterQubits/Unitful.jl/issues/240 +function timeFormat(t::Unitful.Time) + v = ustrip(u"s",t) + if v>3600 + x = Int((v%3600)÷60) + return "$(Int(v÷3600)):$(if x<10; "0" else "" end)$(x) h" + elseif v>60 + x = round(v%60,digits=1) + return "$(Int(v÷60)):$(if x<10; "0" else "" end)$(x) min" + elseif v>0.5 + return "$(round(v,digits=2)) s" + elseif v>0.5e-3 + return "$(round(v*1e3,digits=2)) ms" + else + return "$(round(v*1e6,digits=2)) µs" + end +end \ No newline at end of file diff --git a/test/Scanner/ExtensibleScannerTest.jl b/test/Scanner/ExtensibleScannerTest.jl index 0ff1b967..51e4e401 100644 --- a/test/Scanner/ExtensibleScannerTest.jl +++ b/test/Scanner/ExtensibleScannerTest.jl @@ -16,16 +16,7 @@ MPIMeasurements, but also from the outside. FlexibleDAQParams(dict::Dict) = params_from_dict(FlexibleDAQParams, dict) Base.@kwdef mutable struct FlexibleDAQ <: MPIMeasurements.AbstractDAQ - "Unique device ID for this device as defined in the configuration." - deviceID::String - "Parameter struct for this devices read from the configuration." - params::FlexibleDAQParams - "Flag if the device is optional." - optional::Bool = false - "Flag if the device is present." - present::Bool = false - "Vector of dependencies for this device." - dependencies::Dict{String, Union{Device, Missing}} + MPIMeasurements.@add_device_fields FlexibleDAQParams end function MPIMeasurements.init(daq::FlexibleDAQ) diff --git a/test/Scanner/ScannerConstructionTest.jl b/test/Scanner/ScannerConstructionTest.jl index 7304b91c..5433de5d 100644 --- a/test/Scanner/ScannerConstructionTest.jl +++ b/test/Scanner/ScannerConstructionTest.jl @@ -6,6 +6,7 @@ @test_throws MPIMeasurements.ScannerConfigurationError MPIScanner("TestBrokenDeviceMissingOptional") @test_throws MPIMeasurements.ScannerConfigurationError MPIScanner("TestBrokenDeviceMissingPresent") @test_throws MPIMeasurements.ScannerConfigurationError MPIScanner("TestBrokenDeviceMissingParams") + @test_throws MPIMeasurements.ScannerConfigurationError MPIScanner("TestBrokenDeviceMissingConfigFile") end # Dependencies are checked correctly diff --git a/test/TestConfigs/TestBrokenDeviceMissingConfigFile/Scanner.toml b/test/TestConfigs/TestBrokenDeviceMissingConfigFile/Scanner.toml new file mode 100644 index 00000000..d4a2a3f1 --- /dev/null +++ b/test/TestConfigs/TestBrokenDeviceMissingConfigFile/Scanner.toml @@ -0,0 +1,17 @@ +[General] +boreSize = "1337mm" +facility = "My awesome institute" +manufacturer = "Me, Myself and I" +name = "TestBrokenDeviceMissingConfigFile" +topology = "FFL" +gradient = "42T/m" +datasetStore = "./tmp/TestSimpleSimulatedScannerStore" + +[Devices] +initializationOrder = [ + "testDevice", +] + +[Devices.testDevice] +deviceType = "TestMissingConfigFileDevice" +dependencies = [] diff --git a/test/TestConfigs/TestBrokenDeviceMissingDependencies/Scanner.toml b/test/TestConfigs/TestBrokenDeviceMissingDependencies/Scanner.toml index 43085bea..baf24a3d 100644 --- a/test/TestConfigs/TestBrokenDeviceMissingDependencies/Scanner.toml +++ b/test/TestConfigs/TestBrokenDeviceMissingDependencies/Scanner.toml @@ -6,7 +6,6 @@ name = "TestBrokenDeviceMissingDependencies" topology = "FFL" gradient = "42T/m" datasetStore = "./tmp/TestSimpleSimulatedScannerStore" -defaultSequence = "SimpleSimulatedSequence" [Devices] initializationOrder = [ diff --git a/test/TestConfigs/TestBrokenDeviceMissingID/Scanner.toml b/test/TestConfigs/TestBrokenDeviceMissingID/Scanner.toml index bf1d68ec..88096b60 100644 --- a/test/TestConfigs/TestBrokenDeviceMissingID/Scanner.toml +++ b/test/TestConfigs/TestBrokenDeviceMissingID/Scanner.toml @@ -6,7 +6,6 @@ name = "TestBrokenDeviceMissingID" topology = "FFL" gradient = "42T/m" datasetStore = "./tmp/TestSimpleSimulatedScannerStore" -defaultSequence = "SimpleSimulatedSequence" [Devices] initializationOrder = [ diff --git a/test/TestConfigs/TestBrokenDeviceMissingOptional/Scanner.toml b/test/TestConfigs/TestBrokenDeviceMissingOptional/Scanner.toml index db827780..aa209378 100644 --- a/test/TestConfigs/TestBrokenDeviceMissingOptional/Scanner.toml +++ b/test/TestConfigs/TestBrokenDeviceMissingOptional/Scanner.toml @@ -6,7 +6,6 @@ name = "TestBrokenDeviceMissingOptional" topology = "FFL" gradient = "42T/m" datasetStore = "./tmp/TestSimpleSimulatedScannerStore" -defaultSequence = "SimpleSimulatedSequence" [Devices] initializationOrder = [ diff --git a/test/TestConfigs/TestBrokenDeviceMissingParams/Scanner.toml b/test/TestConfigs/TestBrokenDeviceMissingParams/Scanner.toml index 5c802273..83e0bd1e 100644 --- a/test/TestConfigs/TestBrokenDeviceMissingParams/Scanner.toml +++ b/test/TestConfigs/TestBrokenDeviceMissingParams/Scanner.toml @@ -6,7 +6,6 @@ name = "TestBrokenDeviceMissingParams" topology = "FFL" gradient = "42T/m" datasetStore = "./tmp/TestSimpleSimulatedScannerStore" -defaultSequence = "SimpleSimulatedSequence" [Devices] initializationOrder = [ diff --git a/test/TestConfigs/TestBrokenDeviceMissingPresent/Scanner.toml b/test/TestConfigs/TestBrokenDeviceMissingPresent/Scanner.toml index ed117d40..da19da1b 100644 --- a/test/TestConfigs/TestBrokenDeviceMissingPresent/Scanner.toml +++ b/test/TestConfigs/TestBrokenDeviceMissingPresent/Scanner.toml @@ -6,7 +6,6 @@ name = "TestBrokenDeviceMissingPresent" topology = "FFL" gradient = "42T/m" datasetStore = "./tmp/TestSimpleSimulatedScannerStore" -defaultSequence = "SimpleSimulatedSequence" [Devices] initializationOrder = [ diff --git a/test/TestConfigs/TestDeviceMissingDependencies/Scanner.toml b/test/TestConfigs/TestDeviceMissingDependencies/Scanner.toml index c22e9fa7..e9e8a499 100644 --- a/test/TestConfigs/TestDeviceMissingDependencies/Scanner.toml +++ b/test/TestConfigs/TestDeviceMissingDependencies/Scanner.toml @@ -6,7 +6,6 @@ name = "TestDeviceMissingDependencies" topology = "FFL" gradient = "42T/m" datasetStore = "./tmp/TestSimpleSimulatedScannerStore" -defaultSequence = "SimpleSimulatedSequence" [Devices] initializationOrder = [ diff --git a/test/TestConfigs/TestDeviceUninitDependencies/Scanner.toml b/test/TestConfigs/TestDeviceUninitDependencies/Scanner.toml index 3a4db612..b029f9bc 100644 --- a/test/TestConfigs/TestDeviceUninitDependencies/Scanner.toml +++ b/test/TestConfigs/TestDeviceUninitDependencies/Scanner.toml @@ -6,7 +6,6 @@ name = "TestDeviceUninitDependencies" topology = "FFL" gradient = "42T/m" datasetStore = "./tmp/TestSimpleSimulatedScannerStore" -defaultSequence = "SimpleSimulatedSequence" [Devices] initializationOrder = [ diff --git a/test/TestConfigs/TestDeviceWorkingScanner/Scanner.toml b/test/TestConfigs/TestDeviceWorkingScanner/Scanner.toml index f6866ac9..29ff545b 100644 --- a/test/TestConfigs/TestDeviceWorkingScanner/Scanner.toml +++ b/test/TestConfigs/TestDeviceWorkingScanner/Scanner.toml @@ -6,7 +6,6 @@ name = "TestDeviceWorkingScanner" topology = "FFL" gradient = "42T/m" datasetStore = "./tmp/TestSimpleSimulatedScannerStore" -defaultSequence = "SimpleSimulatedSequence" [Devices] initializationOrder = [ diff --git a/test/TestConfigs/TestDeviceWrongDependencies/Scanner.toml b/test/TestConfigs/TestDeviceWrongDependencies/Scanner.toml index d2074249..eb1b4ecf 100644 --- a/test/TestConfigs/TestDeviceWrongDependencies/Scanner.toml +++ b/test/TestConfigs/TestDeviceWrongDependencies/Scanner.toml @@ -6,7 +6,6 @@ name = "TestDeviceWrongDependencies" topology = "FFL" gradient = "42T/m" datasetStore = "./tmp/TestSimpleSimulatedScannerStore" -defaultSequence = "SimpleSimulatedSequence" [Devices] initializationOrder = [ diff --git a/test/TestConfigs/TestFlexibleScanner/Scanner.toml b/test/TestConfigs/TestFlexibleScanner/Scanner.toml index 3f375828..f8c59283 100644 --- a/test/TestConfigs/TestFlexibleScanner/Scanner.toml +++ b/test/TestConfigs/TestFlexibleScanner/Scanner.toml @@ -6,7 +6,6 @@ name = "TestFlexibleScanner" topology = "FFL" gradient = "42T/m" datasetStore = "./tmp/TestFlexibleScannerStore" -defaultSequence = "" [Devices] initializationOrder = [ diff --git a/test/TestConfigs/TestRedPitayaScanner/Scanner.toml b/test/TestConfigs/TestRedPitayaScanner/Scanner.toml index 66b02096..0127c259 100644 --- a/test/TestConfigs/TestRedPitayaScanner/Scanner.toml +++ b/test/TestConfigs/TestRedPitayaScanner/Scanner.toml @@ -6,7 +6,6 @@ name = "TestRedPitayaScanner" topology = "FFL" gradient = "42T/m" datasetStore = "./tmp/TestRedPitayaScannerStore" -defaultSequence = "" [Devices] initializationOrder = [ diff --git a/test/TestConfigs/TestSimpleSimulatedScanner/Scanner.toml b/test/TestConfigs/TestSimpleSimulatedScanner/Scanner.toml index 8cb74a3e..699147ee 100644 --- a/test/TestConfigs/TestSimpleSimulatedScanner/Scanner.toml +++ b/test/TestConfigs/TestSimpleSimulatedScanner/Scanner.toml @@ -6,7 +6,6 @@ name = "TestSimpleSimulatedScanner" topology = "FFL" gradient = "42T/m" datasetStore = "./tmp/TestSimpleSimulatedScannerStore" -defaultSequence = "SimpleSimulatedSequence" [Devices] initializationOrder = [ diff --git a/test/TestDevices.jl b/test/TestDevices.jl index e22cefb9..b6533c32 100644 --- a/test/TestDevices.jl +++ b/test/TestDevices.jl @@ -46,6 +46,8 @@ Base.@kwdef mutable struct TestDevice <: Device present::Bool = false "Vector of dependencies for this device." dependencies::Dict{String, Union{Device, Missing}} + "Path of the config used to create the device (either Scanner.toml or Device.toml)" + configFile::Union{String, Nothing} = nothing initRan::Bool = false end @@ -73,6 +75,8 @@ Base.@kwdef mutable struct TestDependencyDevice <: Device present::Bool = false "Vector of dependencies for this device." dependencies::Dict{String, Union{Device, Missing}} + "Path of the config used to create the device (either Scanner.toml or Device.toml)" + configFile::Union{String, Nothing} = nothing initRan::Bool = false end @@ -100,6 +104,8 @@ Base.@kwdef mutable struct TestDeviceB <: Device present::Bool = false "Vector of dependencies for this device." dependencies::Dict{String, Union{Device, Missing}} + "Path of the config used to create the device (either Scanner.toml or Device.toml)" + configFile::Union{String, Nothing} = nothing initRan::Bool = false end @@ -123,6 +129,8 @@ Base.@kwdef mutable struct TestMissingIDDevice <: Device present::Bool = false "Vector of dependencies for this device." dependencies::Dict{String, Union{Device, Missing}} + "Path of the config used to create the device (either Scanner.toml or Device.toml)" + configFile::Union{String, Nothing} = nothing initRan::Bool = false end @@ -138,6 +146,8 @@ Base.@kwdef mutable struct TestMissingParamsDevice <: Device present::Bool = false "Vector of dependencies for this device." dependencies::Dict{String, Union{Device, Missing}} + "Path of the config used to create the device (either Scanner.toml or Device.toml)" + configFile::Union{String, Nothing} = nothing initRan::Bool = false end @@ -153,6 +163,8 @@ Base.@kwdef mutable struct TestMissingOptionalDevice <: Device present::Bool = false "Vector of dependencies for this device." dependencies::Dict{String, Union{Device, Missing}} + "Path of the config used to create the device (either Scanner.toml or Device.toml)" + configFile::Union{String, Nothing} = nothing initRan::Bool = false end @@ -168,6 +180,8 @@ Base.@kwdef mutable struct TestMissingPresentDevice <: Device #present::Bool = false "Vector of dependencies for this device." dependencies::Dict{String, Union{Device, Missing}} + "Path of the config used to create the device (either Scanner.toml or Device.toml)" + configFile::Union{String, Nothing} = nothing initRan::Bool = false end @@ -183,6 +197,25 @@ Base.@kwdef mutable struct TestMissingDependencyDevice <: Device present::Bool = false #"Vector of dependencies for this device." #dependencies::Dict{String, Union{Device, Missing}} + "Path of the config used to create the device (either Scanner.toml or Device.toml)" + configFile::Union{String, Nothing} = nothing + + initRan::Bool = false +end + +Base.@kwdef mutable struct TestMissingConfigFileDevice <: Device + "Unique device ID for this device as defined in the configuration." + deviceID::String + "Parameter struct for this devices read from the configuration." + params::TestDeviceParams + "Flag if the device is optional." + optional::Bool = false + "Flag if the device is present." + present::Bool = false + "Vector of dependencies for this device." + dependencies::Dict{String, Union{Device, Missing}} + #"Path of the config used to create the device (either Scanner.toml or Device.toml)" + #configFile::Union{String, Nothing} = nothing initRan::Bool = false end \ No newline at end of file