diff --git a/Project.toml b/Project.toml index a5ae1e8..d44a356 100644 --- a/Project.toml +++ b/Project.toml @@ -47,9 +47,9 @@ HDF5 = "0.15, 0.16" ImageUtils = "0.2" Images = "0.23, 0.24, 0.25" LoggingExtras = ">= 0.4.2" -MPIFiles = "0.15" +MPIFiles = "0.15, 0.16" MPIMeasurements = "0.5" -MPIReco = "0.5" +MPIReco = "0.4, 0.5" MPISphericalHarmonics = "0.0.10" Reexport = "0.2,1.0" CairoMakie = "0.10.2" diff --git a/src/Devices/SurveillanceWidget.jl b/src/Devices/SurveillanceWidget.jl index c1a1134..6a4ae44 100644 --- a/src/Devices/SurveillanceWidget.jl +++ b/src/Devices/SurveillanceWidget.jl @@ -15,7 +15,7 @@ function TemperatureLog(filename::String) p = TOML.parsefile(filename) return TemperatureLog(p["temperatures"], p["times"], p["numChan"]) elseif ext == ".mdf" - temps = h5read(filename, "/measurement/_temperatures") + temps = h5read(filename, "/measurement/_monitoring/temperature/observed") times = [DateTime(0)+Dates.Second(1)*j for j=1:size(temps,2)] return TemperatureLog(vec(temps), times, size(temps,1)) else diff --git a/src/Devices/TemperatureControlWidget.jl b/src/Devices/TemperatureControlWidget.jl index ffc92f5..67e40b7 100644 --- a/src/Devices/TemperatureControlWidget.jl +++ b/src/Devices/TemperatureControlWidget.jl @@ -32,7 +32,7 @@ function TemperatureControllerWidget(tempCont::TemperatureController) #push!(m, m.canvas) #set_gtk_property!(m,:expand, m.canvas, true) - showall(m) + show(m) updateTarget(m) updateMaximum(m) diff --git a/src/Devices/TemperatureSensorWidget.jl b/src/Devices/TemperatureSensorWidget.jl index a928e19..1e05c36 100644 --- a/src/Devices/TemperatureSensorWidget.jl +++ b/src/Devices/TemperatureSensorWidget.jl @@ -71,7 +71,7 @@ function initCallbacks(m::TemperatureSensorWidget) end signal_connect(m["btnLoadTemp"], :clicked) do w - filter = Gtk.GtkFileFilter(pattern=String("*.toml, *.mdf"), mimetype=String("application/toml")) + filter = Gtk4.GtkFileFilter(pattern=String("*.toml, *.mdf"), mimetype=String("application/toml")) diag = open_dialog("Select Temperature File", mpilab[]["mainWindow"], (filter, )) do filename m.updating = true if filename != "" diff --git a/src/MPIUI.jl b/src/MPIUI.jl index 852623a..a98bea8 100644 --- a/src/MPIUI.jl +++ b/src/MPIUI.jl @@ -20,6 +20,7 @@ using MPISphericalHarmonics, SphericalHarmonicExpansions # for MagneticFieldView using NLsolve # for MagneticFieldViewer: findFFP() using DataFrames, CSV # for MagneticFieldViewer: export as csv using Unitful +using MPIFiles.HDF5 #open MDF files import CairoMakie ENV["MPILIB_UI"] = "Nothing" diff --git a/src/Protocol/DataHandler/OnlineRecoHandler.jl b/src/Protocol/DataHandler/OnlineRecoHandler.jl index 0aec068..b0a76d2 100644 --- a/src/Protocol/DataHandler/OnlineRecoHandler.jl +++ b/src/Protocol/DataHandler/OnlineRecoHandler.jl @@ -190,7 +190,7 @@ end if m.freq != nothing # here we merge frequencies and channels - data = reshape(data, size(data,1)*size(data,2), size(data,3), size(data,4)) + # data = reshape(data, size(data,1)*size(data,2), size(data,3), size(data,4)) data = data[m.freq, :, :] end diff --git a/src/Reconstruction/ReconstructionParameter.jl b/src/Reconstruction/ReconstructionParameter.jl index a89c474..92e95c2 100644 --- a/src/Reconstruction/ReconstructionParameter.jl +++ b/src/Reconstruction/ReconstructionParameter.jl @@ -179,15 +179,29 @@ function getParams(m::ReconstructionParameter) params[:regName] = ["L1", "TV"] end=# - params[:solver] = Kaczmarz + #params[:solver] = Kaczmarz +# + #if params[:solver] == Kaczmarz + # if params[:lambdaL1] == 0.0 + # params[:reg] = AbstractRegularization[L2Regularization(Float32(params[:lambd]))] + # else + # params[:reg] = AbstractRegularization[L2Regularization(Float32(params[:lambd])), L1Regularization(Float32(params[:lambdaL1]))] + # end + # append!(params[:reg], [PositiveRegularization(), RealRegularization()]) + #end + + params[:solver] = linearSolverList()[max(get_gtk_property(m["cbSolver"],:active, Int64) + 1,1)] + # Small hack + if params[:solver] == "fusedlasso" + params[:loadasreal] = true + params[:lambd] = [params[:lambdaL1], params[:lambdaTV]] + params[:regName] = ["L1", "TV"] + end - if params[:solver] == Kaczmarz - if params[:lambdaL1] == 0.0 - params[:reg] = AbstractRegularization[L2Regularization(Float32(params[:lambd]))] - else - params[:reg] = AbstractRegularization[L2Regularization(Float32(params[:lambd])), L1Regularization(Float32(params[:lambdaL1]))] - end - append!(params[:reg], [PositiveRegularization(), RealRegularization()]) + if params[:solver] == "kaczmarz" + #params[:loadasreal] = true + params[:lambd] = [params[:lambd], params[:lambdaL1] ] #params[:lambdaTV], params[:lambdaL1]] + params[:regName] = ["L2", "L1"] #"TV", "L1"] end firstFrame = get_gtk_property(m["adjFrame"], :value, Int64) diff --git a/src/Viewer/RawDataViewer.jl b/src/Viewer/RawDataViewer.jl index 2383a47..2be3597 100644 --- a/src/Viewer/RawDataViewer.jl +++ b/src/Viewer/RawDataViewer.jl @@ -244,8 +244,21 @@ end # TODO: Ensure that the measurements fit together (num samples / patches) # otherwise -> error - numFGFrames = minimum(acqNumFGFrames.(fs)) - numBGFrames = minimum(acqNumBGFrames.(fs)) + numFGFrames = unique(acqNumFGFrames.(fs)) + if length(numFGFrames) > 1 + numFGFrames = minimum(numFGFrames) + @warn "Different number of foreground frames in data, limit frames to $numFGFrames" + else + numFGFrames = first(numFGFrames) + end + + numBGFrames = unique(acqNumBGFrames.(fs)) + if length(numBGFrames) > 1 + numBGFrames = minimum(numBGFrames) + @warn "Different number of background frames in data, limit frames to $numBGFrames" + else + numBGFrames = first(numBGFrames) + end dataFGVec = Any[] dataBGVec = Any[] @@ -271,7 +284,7 @@ end tfCorrection=get_gtk_property(m["cbCorrTF"], :active, Bool)) push!(dataFGVec, data) - if acqNumBGFrames(f) > 0 + if numBGFrames > 0 dataBG = getMeasurements(f, false, frames=measBGFrameIdx(f), bgCorrection=false, spectralLeakageCorrection = get_gtk_property(m["cbSLCorr"], :active, Bool), tfCorrection=get_gtk_property(m["cbCorrTF"], :active, Bool)) @@ -381,42 +394,30 @@ end numFreq = floor(Int, size(data,1) ./ 2 .+ 1) maxPoints = 5000 - sp = length(minTP:maxTP) > maxPoints ? round(Int,length(minTP:maxTP) / maxPoints) : 1 - - steps = minTP:sp:maxTP - dataCompressed = zeros(length(steps), size(data,2)) - if sp > 1 - for j=1:size(data,2) - for l=1:length(steps) - st = steps[l] - en = min(st+sp,steps[end]) - med_ = median(data[st:en,j]) - max_ = maximum(data[st:en,j]) - min_ = minimum(data[st:en,j]) - dataCompressed[l,j] = rand(Bool) ? max_ : min_ #abs(max_ - med_) > abs(med_ - min_) ? max_ : min_ - end - end - else - dataCompressed = data[steps,:] - end + # Reduce each channel to either min or max value if intervals of more than maxpoint points + reduceOp = (vals -> map(slice -> rand(Bool) ? maximum(slice) : minimum(slice), eachslice(vals, dims=2))) + timeCompressed, dataCompressed = prepareTimeDataForPlotting(data, timePoints, indexInterval=(minTP,maxTP), + maxpoints=maxPoints, reduceOp = reduceOp) - fTD, axTD, lTD1 = CairoMakie.lines(timePoints[steps], dataCompressed[:,1], + fTD, axTD, lTD1 = CairoMakie.lines(timeCompressed, dataCompressed[:,1], figure = (; figure_padding=4, resolution = (1000, 800), fontsize = 11), axis = (; title = "Time Domain"), color = CairoMakie.RGBf(colors[1]...), label = labels_[1]) for j=2:size(data,2) - CairoMakie.lines!(axTD, timePoints[steps],dataCompressed[:,j], + CairoMakie.lines!(axTD, timeCompressed, dataCompressed[:,j], color = CairoMakie.RGBf(colors[j]...), #linewidth=3) label = labels_[j]) end if length(m.dataBG) > 0 && get_gtk_property(m["cbShowBG"], :active, Bool) - CairoMakie.lines!(axTD, timePoints[minTP:sp:maxTP],dataBG[minTP:sp:maxTP,1], color=:black, + _, dataBGCompressed = prepareTimeDataForPlotting(dataBG, timePoints, indexInterval=(minTP,maxTP), + maxpoints=maxPoints, reduceOp = reduceOp) + CairoMakie.lines!(axTD, timeCompressed, dataBGCompressed[:, 1], color=:black, label="BG", linestyle = :dash) end CairoMakie.autolimits!(axTD) - if timePoints[steps[end]] > timePoints[steps[1]] - CairoMakie.xlims!(axTD, timePoints[steps[1]], timePoints[steps[end]]) + if timeCompressed[end] > timeCompressed[1] + CairoMakie.xlims!(axTD, timeCompressed[1], timeCompressed[end]) end if !autoRangingTD && maxValTD > minValTD CairoMakie.ylims!(axTD, minValTD, maxValTD) @@ -436,33 +437,26 @@ end spFr = length(minFr:maxFr) > maxPoints ? round(Int,length(minFr:maxFr) / maxPoints) : 1 stepsFr = minFr:spFr:maxFr - freqDataCompressed = zeros(length(stepsFr), size(freqdata,2)) - if spFr > 1 - for j=1:size(freqdata,2) - for l=1:length(stepsFr) - st = stepsFr[l] - en = min(st+spFr,stepsFr[end]) - freqDataCompressed[l,j] = maximum(freqdata[st:en,j]) - end - end - else - freqDataCompressed = freqdata[stepsFr,:] - end + freqsCompressed, freqDataCompressed = prepareTimeDataForPlotting(freqdata, freq, indexInterval=(minFr,maxFr), + maxpoints=maxPoints, reduceOp = vals -> map(maximum, eachslice(vals, dims=2))) + - fFD, axFD, lFD1 = CairoMakie.lines(freq[stepsFr],freqDataCompressed[:,1], + fFD, axFD, lFD1 = CairoMakie.lines(freqsCompressed,freqDataCompressed[:,1], figure = (; figure_padding=4, resolution = (1000, 800), fontsize = 11), axis = (; title = "Frequency Domain", yscale=log10), color = CairoMakie.RGBf(colors[1]...), label=labels_[1]) for j=2:size(data,2) - CairoMakie.lines!(axFD, freq[stepsFr], freqDataCompressed[:,j], + CairoMakie.lines!(axFD, freqsCompressed, freqDataCompressed[:,j], color = CairoMakie.RGBf(colors[j]...), label=labels_[j]) end if length(m.dataBG) > 0 && get_gtk_property(m["cbShowBG"], :active, Bool) - CairoMakie.lines!(axTD, timePoints[minTP:sp:maxTP],dataBG[minTP:sp:maxTP,1], color=:black, - label="BG", linestyle = :dash) + #CairoMakie.lines!(axTD, timePoints[minTP:sp:maxTP],dataBG[minTP:sp:maxTP,1], color=:black, + # label="BG", linestyle = :dash) # isn't this duplicate? if showFD - CairoMakie.lines!(axFD, freq[minFr:spFr:maxFr],abs.(rfft(dataBG,1)[minFr:spFr:maxFr,1]) / size(dataBG,1), + _, freqDataBGCompressed = prepareTimeDataForPlotting(abs.(rfft(dataBG,1)), freq, indexInterval=(minFr,maxFr), + maxpoints=maxPoints, reduceOp = vals -> map(maximum, eachslice(vals, dims=2))) + CairoMakie.lines!(axFD, freqsCompressed, freqDataBGCompressed[:, 1] / size(dataBG,1), color=:black, label="BG", linestyle = :dash) end end diff --git a/src/Viewer/Viewer.jl b/src/Viewer/Viewer.jl index 9dacb40..dedb0af 100644 --- a/src/Viewer/Viewer.jl +++ b/src/Viewer/Viewer.jl @@ -1,3 +1,63 @@ +""" + prepareTimeDataForPlotting(data::AbstractMatrix, times = 1:length(data); timeInterval = nothing, maxpoints = length(data), reduceOp = nothing) + +Reduces time-series `data` to at most `maxpoints` by either resampling or reducing the data with `reduceOp`. The time-series is defined by `times` and `timeInterval`. +If `timeInterval` is not given, the whole time-series is used. If `reduceOp` is not given, the data is resampled. +The function returns the reduced time-series `t` and the reduced data `d`. +""" +function prepareTimeDataForPlotting(data::AbstractMatrix, times = 1:length(data); indexInterval = nothing, timeInterval = nothing, maxpoints = length(data), reduceOp = nothing) + idx1, idx2 = timeSlice(times, indexInterval, timeInterval) + if length(idx1:idx2) > maxpoints + t, d = reducePoints(reduceOp, data, times, idx1, idx2, maxpoints) + else + t = times[idx1:idx2] + d = data[idx1:idx2, :] + end + # return [[Point2f(t[l], d[l, k]) for l in 1:size(d,1)] for k in 1:size(d, 2)] + return t, d +end +# Main variant is defined for matrices, so we reshape into Nx1 matrix +function prepareTimeDataForPlotting(data::AbstractVector, args...; kwargs...) + # return first(prepareTimeDataForPlotting(reshape(data, :, 1), args...; kwargs...)) + t, d = prepareTimeDataForPlotting(reshape(data, :, 1), args...; kwargs...) + return t, vec(d) +end + + +timeSlice(times, indexInterval, timeInterval) = error("Cannot specify both indexInterval and timeInterval") +timeSlice(times, indexInterval, ::Nothing) = first(indexInterval), last(indexInterval) +function timeSlice(times, ::Nothing, timeInterval) + t1, t2 = timeInterval + idx1 = findfirst(t -> (t>t1), times) + idx2 = findlast(t -> (t min(i + stepsize -1, idx2), steps))) + d[i, :] = reduceOp(@view data[st:en, :]) + #@info "Reducing data from $st to $en with $(d[i, :])" + end + return t, d +end +function reducePoints(::Nothing, data, times, idx1, idx2, maxpoints) + d = DSP.resample(data[idx1:idx2], maxpoints / length(idx1:idx2)) + t = range(times[idx1], times[idx2], length=length(d)) + return t, d +end + include("BaseViewer.jl") include("SimpleDataViewer.jl") include("DataViewer/DataViewer.jl") diff --git a/src/builder/temperatureControllerWidget.ui b/src/builder/temperatureControllerWidget.ui index dbcb450..6de370d 100644 --- a/src/builder/temperatureControllerWidget.ui +++ b/src/builder/temperatureControllerWidget.ui @@ -1,7 +1,7 @@ - + True False @@ -18,11 +18,6 @@ True True - - False - True - 0 - @@ -31,18 +26,8 @@ True True - - False - True - 1 - - - False - True - 0 - @@ -54,24 +39,24 @@ True False + + 1 + 0 + 2 + - - 1 - 0 - 2 - True False Target Temps + + 0 + 0 + 2 + - - 0 - 0 - 2 - @@ -79,11 +64,11 @@ True True True + + 2 + 0 + - - 2 - 0 - @@ -91,34 +76,34 @@ True True True + + 2 + 1 + - - 2 - 1 - True False Maximum Temps + + 0 + 2 + 2 + - - 0 - 2 - 2 - True False + + 1 + 2 + 2 + - - 1 - 2 - 2 - @@ -126,11 +111,11 @@ True True True + + 2 + 2 + - - 2 - 2 - @@ -138,11 +123,11 @@ True True True + + 2 + 3 + - - 2 - 3 - @@ -150,22 +135,22 @@ True True False + + 3 + 0 + - - 3 - 0 - True True True + + 3 + 1 + - - 3 - 1 - @@ -173,29 +158,24 @@ True True False + + 3 + 2 + - - 3 - 2 - True True False + + 3 + 3 + - - 3 - 3 - - - False - True - 1 - @@ -209,11 +189,6 @@ True True - - False - True - 0 - @@ -222,11 +197,6 @@ True True - - False - True - 1 - @@ -235,11 +205,6 @@ True True - - False - True - 2 - @@ -248,18 +213,8 @@ True True - - False - True - 3 - - - False - True - 2 - @@ -268,11 +223,6 @@ True 3 - - False - True - 3 -