update scripts

benchmark1/latex/reco_images.tex

@@ -0,0 +1,95 @@
+\documentclass{standalone} %[varwidth]
+
+%\def\size{0.15\textwidth}
+%\def\sizespace{0.18\textwidth}
+
+\usepackage{tikz}
+\usetikzlibrary{calc}
+\usepackage{tikzscale}
+\usetikzlibrary{positioning,shapes,shapes.multipart}
+\usepackage{pgfplots}
+\pgfplotsset{compat=1.5.1}
+
+\tikzset{
+  ../pics/square/.default={1},
+  ../pics/square/.style = {
+    code = {
+    \draw[pic actions] (0,0) rectangle (#1,0.6#1);
+    }
+  }
+} 
+
+\pgfplotsset{width=4.5cm,compat=1.5.1}
+\newcommand{\zb}[1]{\mbox{\boldmath{${#1}$}}}
+\newcommand{\zbs}[1]{\mbox{\boldmath\scriptsize{${#1}$}}}
+
+\begin{document}
+
+\begin{tikzpicture}
+\tikzstyle{nodered} = [inner sep = 1pt, draw=red, line width=1pt, fill=red]
+\tikzstyle{nodewhite} = [inner sep = 1pt, draw=white, line width=1pt, fill=white]
+
+%& Conc Series &$\lambda^1$ & $\lambda^1$ windowed & $\lambda^2$ windowed & mask & new reco \\
+\node[inner sep=1pt] (conc1) at (0,0) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_bart_rf1.png}};
+\node[right = 1pt of conc1, nodewhite] (conc2) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_bart_rf2.png}};
+\node[right = 1pt of conc2, nodewhite] (conc3) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_bart_rf3.png}};
+\node[right = 1pt of conc3, nodewhite] (conc4) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_bart_rf4.png}};
+
+% Labels above figures
+\node[above = -0.2mm of conc1]{$R=1$};
+\node[above = -0.2mm of conc2]{$R=2$};
+\node[above = -0.2mm of conc3]{$R=3$};
+\node[above = -0.2mm of conc4]{$R=4$};
+
+\node[rotate=90, left = 2mm of conc1, xshift=5mm]{BART};
+
+\node[inner sep=1pt, below = 2pt of conc1] (conc1) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_diff_bart_rf1.png}};
+\node[right = 1pt of conc1, nodewhite] (conc2) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_diff_bart_rf2.png}};
+\node[right = 1pt of conc2, nodewhite] (conc3) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_diff_bart_rf3.png}};
+\node[right = 1pt of conc3, nodewhite] (conc4) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_diff_bart_rf4.png}};
+\node[rotate=90, left = 2mm of conc1, xshift=9mm]{Diff BART};
+
+\node[inner sep=1pt, below=2pt of conc1] (conc1) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_mrireco_toeplitz1_oversamp2.0_rf1.png}};
+\node[right = 1pt of conc1, nodewhite] (conc2) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_mrireco_toeplitz1_oversamp2.0_rf2.png}};
+\node[right = 1pt of conc2, nodewhite] (conc3) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_mrireco_toeplitz1_oversamp2.0_rf3.png}};
+\node[right = 1pt of conc3, nodewhite] (conc4) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_mrireco_toeplitz1_oversamp2.0_rf4.png}};
+\node[rotate=90, left = 2mm of conc1, xshift=8mm]{MRIReco};
+
+
+\node[inner sep=1pt, below = 2pt of conc1] (conc1) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_diff_mrireco_rf1.png}};
+\node[right = 1pt of conc1, nodewhite] (conc2) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_diff_mrireco_rf2.png}};
+\node[right = 1pt of conc2, nodewhite] (conc3) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_diff_mrireco_rf3.png}};
+\node[right = 1pt of conc3, nodewhite] (conc4) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_diff_mrireco_rf4.png}};
+\node[rotate=90, left = 2mm of conc1, xshift=12mm]{Diff MRIReco};
+
+\node[inner sep=1pt, below = 2pt of conc1] (conc1) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_diff_rf1.png}};
+\node[right = 1pt of conc1, nodewhite] (conc2) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_diff_rf2.png}};
+\node[right = 1pt of conc2, nodewhite] (conc3) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_diff_rf3.png}};
+\node[right = 1pt of conc3, nodewhite] (conc4) 
+{\includegraphics[width=0.2\textwidth]{../reco/imgCG_diff_rf4.png}};
+\node[rotate=90, left = 2mm of conc1, xshift=9mm]{Diff Both};
+
+
+\centering
+\end{tikzpicture}
+
+\end{document}

benchmark1/latex/timings.tex

@@ -0,0 +1,109 @@
+\documentclass{standalone} %[varwidth]
+\usepackage{pgfplots}
+\usepackage{pgfplotstable}
+\usetikzlibrary{calc}
+%s\usepackage{siunitx}
+%\def\size{0.15\textwidth}
+%\def\sizespace{0.18\textwidth}
+
+\usepackage{tikz}
+\usepackage{tikzscale}
+\usetikzlibrary{positioning,shapes,shapes.multipart}
+\usepackage{pgfplots}
+\pgfplotsset{compat=1.5.1}
+
+\tikzset{
+  ../../../pics/square/.default={1},
+  ../../../pics/square/.style = {
+    code = {
+    \draw[pic actions] (0,0) rectangle (#1,0.6#1);
+    }
+  }
+} 
+
+\pgfplotsset{width=4.5cm,compat=1.5.1}
+\newcommand{\zb}[1]{\mbox{\boldmath{${#1}$}}}
+\newcommand{\zbs}[1]{\mbox{\boldmath\scriptsize{${#1}$}}}
+
+\definecolor{ibilight}{RGB}{193,216,237}
+\definecolor{ibidark}{RGB}{0,73,146}
+\definecolor{uke2}{RGB}{170,156,143}
+\definecolor{uke3}{RGB}{87,87,86}
+\definecolor{ukesec1}{RGB}{255,223,0}
+\definecolor{ukesec2}{RGB}{239,123,5}
+\definecolor{ukesec3}{RGB}{104,195,205}
+\definecolor{ukesec4}{RGB}{138,189,36}
+\definecolor{tuhh}{RGB}{45,198,214}
+
+\begin{document}
+
+
+
+\begin{tikzpicture}[node distance=-0.05cm and -0.05cm]%[scale=1, every node/.style={scale=0.8}]
+\def\sizex{6cm}
+\def\sizey{5cm}
+
+\pgfplotsset{footnotesize,
+             grid=both,
+             width=\sizex,
+             height=\sizey,
+             extra y ticks={1},
+             axis background/.style={fill=white},
+             ylabel=Time [s],
+             legend style={%at={(0.99,0.98)},anchor=north east,
+                           at={(0.35,0.7)},anchor=west,
+						   cells={anchor=west},}},
+
+\tikzset{every mark/.append style={scale=0.5}}
+
+
+
+\begin{axis}[
+    xtick={1,4,8,12},
+    extra y tick labels={},
+    title={$R=1$},
+    xlabel=\# Threads
+    ]
+
+    \pgfplotstableread[col sep=comma]{../reco/recoTimes.csv}\DataBart
+    \pgfplotstableread[col sep=comma]{../reco/recoTimes_mrireco_toeplitz1_oversamp2.0.csv}\DataMRIReco
+    \pgfplotstableread[col sep=comma]{../reco/recoTimes_mrireco_toeplitz0_oversamp1.25.csv}\DataMRIRecoOver
+
+    \addplot [ukesec2, mark=*,line width=2] table [x expr=\thisrowno{0}, y expr=\thisrowno{1}] {\DataBart};
+    \addlegendentry{BART };
+
+      \addplot [ibidark, mark=*,line width=2] table [x expr=\thisrowno{0}, y expr=\thisrowno{1}] {\DataMRIReco};
+    \addlegendentry{MRIReco};
+
+    \addplot [ukesec3, mark=*,line width=2] table [x expr=\thisrowno{0}, y expr=\thisrowno{1}] {\DataMRIRecoOver};
+     \addlegendentry{MRIReco*};
+
+\end{axis}
+
+
+\begin{axis}[
+    xtick={1,4,8,12},
+    extra y tick labels={},
+    xlabel=\# Threads,
+    title={$R=4$}, xshift=\sizex
+    ]
+
+    \pgfplotstableread[col sep=comma]{../reco/recoTimes.csv}\DataBart
+    \pgfplotstableread[col sep=comma]{../reco/recoTimes_mrireco_toeplitz1_oversamp2.csv}\DataMRIReco
+    \pgfplotstableread[col sep=comma]{../reco/recoTimes_mrireco_toeplitz0_oversamp1.25.csv}\DataMRIRecoOver
+
+    \addplot [ukesec2, mark=*,line width=2] table [x expr=\thisrowno{0}, y expr=\thisrowno{4}] {\DataBart};
+    %\addlegendentry{BART };
+
+      \addplot [ibidark, mark=*,line width=2] table [x expr=\thisrowno{0}, y expr=\thisrowno{4}] {\DataMRIReco};
+    %\addlegendentry{MRIReco};
+
+    \addplot [ukesec3, mark=*,line width=2] table [x expr=\thisrowno{0}, y expr=\thisrowno{4}] {\DataMRIRecoOver};
+
+
+\end{axis}
+
+\end{tikzpicture}
+
+
+\end{document}

benchmark1/makeImages.jl

@@ -0,0 +1,95 @@
+using ImageUtils, HDF5, LinearAlgebra, RegularizedLeastSquares
+
+
+imagesBart = Any[]
+imagesMRIReco = Any[]
+imagesMRIReco2 = Any[]
+
+rf = [1,2,3,4]
+
+f_img  = "./reco/imgCG_bart.h5"
+for i = 1:length(rf)
+  h5open(f_img, "r") do file
+    A = read(file, "/rf$(rf[i])_re")
+    B = read(file, "/rf$(rf[i])_im")
+    im_ = A.+im.*B
+    push!(imagesBart, im_)
+    filename = "./reco/imgCG_bart_rf$(rf[i]).png"
+    exportImage(filename, abs.(im_), colormap="viridis")
+  end
+end
+
+
+f_img  = "./reco/imgCG_mrireco_toeplitz1_oversamp2.0.h5"
+for i = 1:length(rf)
+  h5open(f_img, "r") do file
+    A = read(file, "/rf$(rf[i])_re")
+    B = read(file, "/rf$(rf[i])_im")
+    im_ = A.+im.*B
+    push!(imagesMRIReco, im_)
+    filename = "./reco/imgCG_mrireco_toeplitz1_oversamp2.0_rf$(rf[i]).png"
+    exportImage(filename, abs.(im_),colormap="viridis" )
+  end
+end
+
+f_img  = "./reco/imgCG_mrireco_toeplitz0_oversamp1.25.h5"
+for i = 1:length(rf)
+  h5open(f_img, "r") do file
+    A = read(file, "/rf$(rf[i])_re")
+    B = read(file, "/rf$(rf[i])_im")
+    im_ = A.+im.*B
+    push!(imagesMRIReco2, im_)
+    filename = "./reco/imgCG_mrireco_toeplitz0_oversamp1.25_rf$(rf[i]).png"
+    exportImage(filename, abs.(im_), colormap="viridis" )
+  end
+end
+
+function optimalScaling(I,Ireco)
+  N = length(I)
+
+  # This is a little trick. We usually are not interested in simple scalings
+  # and therefore "calibrate" them away
+  alpha = norm(Ireco)>0 ? (dot(vec(I),vec(Ireco))+dot(vec(Ireco),vec(I))) /
+          (2*dot(vec(Ireco),vec(Ireco))) : 1.0
+  I2 = Ireco.*alpha
+
+  return I2
+end
+
+A1 = abs.(imagesMRIReco[1])
+A1 .= A1 ./maximum(A1)
+
+B1 = abs.(imagesBart[1])
+B1 .= B1 ./maximum(B1)
+
+for i = 1:length(rf)
+  A = abs.(imagesMRIReco[i])
+  A .= A ./maximum(A)
+
+  B = abs.(imagesBart[i])
+  B .= optimalScaling(A,B)
+
+  D = A.-B 
+  filename = "./reco/imgCG_diff_rf$(rf[i]).png"
+  exportImage(filename, D, colormap="viridis", vmin=-0.1, vmax=0.1, normalize=false )
+
+  ###
+  A .= optimalScaling(A1,A)
+  D = A.-A1 
+  filename = "./reco/imgCG_diff_mrireco_rf$(rf[i]).png"
+  exportImage(filename, D, colormap="viridis", vmin=-0.1, vmax=0.1, normalize=false )  
+
+  B .= optimalScaling(B1,B)
+  D = B.-B1 
+  filename = "./reco/imgCG_diff_bart_rf$(rf[i]).png"
+  exportImage(filename, D, colormap="viridis", vmin=-0.1, vmax=0.1, normalize=false )  
+
+  @info nrmsd(A,B)   nrmsd(abs.(imagesMRIReco[1]),A)  nrmsd(abs.(imagesBart[1]),B)
+
+end
+
+
+
+
+
+

benchmark1/recoBrainBart.m

@@ -38,9 +38,10 @@
 end
 
 %% write output images and reco times to files
-f_times = './reco/recoTimes_bart.csv';
-f_img = './reco/imgCG_bart.h5';
-dlmwrite(f_times, times, '-append');
+f_times = ['./reco/recoTimes_bart.csv'];
+f_img = ['./reco/imgCG_bart.h5'];
+nthreads = str2double(getenv('OMP_NUM_THREADS'));
+dlmwrite(f_times, horzcat(nthreads,times), '-append');
 if ~isfile(f_img)
     for i=1:length(rf)
         h5create(f_img,['/rf' num2str(rf(i)) '_re'], size(img_cg{i}));

benchmark1/recoBrainMRIReco.jl

@@ -1,11 +1,14 @@
-using HDF5, MRIReco, LinearAlgebra, DelimitedFiles, BenchmarkTools
-FFTW.set_num_threads(1);BLAS.set_num_threads(1)
+using HDF5, MRIReco, DelimitedFiles, BenchmarkTools
 
 filename = "./data/rawdata_brain_radial_96proj_12ch.h5"
 data = permutedims(h5read(filename, "rawdata"),[3,2,1,4])
 traj = permutedims(h5read(filename, "trajectory"),[3,2,1])
 N = 300
 Nc = 12
+toeplitz = parse(Int,get(ENV,"TOEPLITZ","0"))
+oversamplingFactor = parse(Float64,get(ENV,"OVERSAMPLING","1.25"))
+
+@info "Threads = $(Threads.nthreads()) Toeplitz=$(toeplitz)  OversamplingFactor=$(oversamplingFactor)"   
 
 #############################################
 # load data and form Acquisition data object
@@ -33,6 +36,8 @@ params[:regularization] = "L2"
 params[:λ] = 1.e-2
 params[:iterations] = 100
 params[:solver] = "cgnr"
+params[:toeplitz] = toeplitz == 1
+params[:oversamplingFactor] = oversamplingFactor
 params[:senseMaps] = reshape(sensitivity, N, N, 1, Nc)
 
 # @time img_ref = reconstruction(acqData, params).data
@@ -49,20 +54,21 @@ params[:relTol] = 0.0
 for i = 1:length(rf)
   @info "r=$(rf[i])"
   # undersample profiles
-  acqDataSub = convertUndersampledData(sample_kspace(acqData, Float64.(rf[i]), "regular"))
+  global acqDataSub = convertUndersampledData(sample_kspace(acqData, Float64.(rf[i]), "regular"))
   # SENSE reconstruction while monitoring error
   # run twice to take factor out precompilation effects
   img_cg[i] = reconstruction(acqDataSub, params).data
-  times[i] = @elapsed img_cg[i] = reconstruction(acqDataSub, params).data
+  times[i] = @belapsed reconstruction(acqDataSub, params).data
 end
 
 ##############
 # write output
 ##############
-f_times = "./reco/recoTimes_mrireco.csv"
-f_img  = "./reco/imgCG_mrireco.h5"
+f_times = "./reco/recoTimes_mrireco_toeplitz$(toeplitz)_oversamp$(oversamplingFactor).csv"
+f_img  = "./reco/imgCG_mrireco_toeplitz$(toeplitz)_oversamp$(oversamplingFactor).h5"
+
 open(f_times,"a") do file
-  writedlm(file, transpose(times))
+  writedlm(file, hcat(Threads.nthreads(), transpose(times)), ',')
 end
 
 if !isfile(f_img)
@@ -74,4 +80,4 @@ if !isfile(f_img)
   end
 end
 
-exit()
+exit()