1# SPDX-FileCopyrightText: 2020 Jonathan Pieper <jpieper@stud.uni-frankfurt.de> 

2# 

3# SPDX-License-Identifier: GPL-3.0-or-later 

4 

5import unittest 

6 

7import logging # System Modules 

8import os 

9from glob import glob 

10 

11# Basic Plotting libraries 

12import matplotlib.pyplot as plt 

13import matplotlib 

14import seaborn as sns 

15from matplotlib import cm 

16from matplotlib.colors import LogNorm 

17from mpl_toolkits.axes_grid1.inset_locator import inset_axes 

18 

19# Math / Science Libraries 

20import pandas as pd 

21import numpy as np 

22import scipy 

23import scipy.stats 

24 

25import ana 

26 

27logging.basicConfig(level=logging.ERROR) 

28 

29 

30def load_data(datapath, focus='Vin', cut=-2): 

31 meas_data = {} 

32 meas_info = {} 

33 all_data = {} 

34 for f in datapath: 

35 f_info = ana.measurement.MeasurementClass().get_info_from_name(f) 

36 sr = f_info[focus] 

37 nr = f_info['nr'] 

38 meas_info[sr] = f_info 

39 meas_data[sr] = pd.read_csv(f, sep='\t') 

40 new_df = meas_data[sr] 

41 new_df[focus] = float(sr[:cut]) 

42 if nr in all_data.keys(): 

43 all_data[nr] = pd.concat([all_data[nr], new_df]) 

44 else: 

45 all_data[nr] = new_df 

46 return meas_data, meas_info, all_data 

47 

48 

49def calc_PSD(meas_data): 

50 meas_obj = {} 

51 for sr, data_df in meas_data.items(): 

52 if len(data_df['Vx']) % 1024: 

53 avg = len(data_df['Vx']) // 1024 

54 d = data_df['Vx'].iloc[:-(len(data_df['Vx']) % 1024)] 

55 else: 

56 d = data_df.Vx 

57 

58 max_len = len(d) 

59 

60 data = { 

61 'data': d, 

62 'info': { 

63 'Nr': 0, 

64 'rate': 1 / data_df.time.diff().mean(), 

65 'length': max_len * data_df.time.diff().mean(), 

66 } 

67 } 

68 

69 meas_obj[sr] = ana.RAW(data, 

70 rate=data['info']['rate'], 

71 nof_first_spectra=32, 

72 calc_first=True, 

73 downsample=False, 

74 ) 

75 return meas_obj 

76 

77 

78def merge_data(meas_obj, cutoff_frequency=.9, cut=-2): 

79 diff_voltages = pd.DataFrame() 

80 for sr, m in meas_obj.items(): 

81 s = m.avg_spec 

82 s = s[s.freq < cutoff_frequency] 

83 if len(s) < 2: 

84 continue 

85 newdf = pd.DataFrame() 

86 newdf['freq'] = s.freq 

87 newdf['S'] = s.S 

88 newdf['Vin'] = float(sr[:cut]) 

89 diff_voltages = pd.concat([diff_voltages, newdf]) 

90 return diff_voltages 

91 

92 

93def plot_PSD_classic(diff_voltages, title, groupby_category='SR', group_name='T/min', 

94 num=10, style=[['science'], {'context': 'talk', 'style': 'white', 'palette': 'bright', }]): 

95 c1 = sns.color_palette("hls", num) 

96 sns.set_palette(c1) 

97 fig, ax = plt.subplots(figsize=(16, 12)) 

98 # g = sns.relplot(x='freq', y='S', hue='Vin', data=diff_voltages, height=5, kind='line') 

99 grouped = diff_voltages.groupby(groupby_category) 

100 for group in grouped.groups.keys(): 

101 if group == 0: 

102 grouped.get_group(group).plot(x='freq', y='S', kind='line', 

103 loglog=True, ax=ax, 

104 label='%.1e %s' % (group, group_name), 

105 color='black', 

106 xlabel='Frequency [Hz]', 

107 ylabel='$S_{V_H}$ [$\\mathrm{V}^2/\\mathrm{Hz}$]', 

108 ) 

109 else: 

110 grouped.get_group(group).plot(x='freq', y='S', kind='line', 

111 loglog=True, ax=ax, 

112 label='%.1e %s' % (group, group_name), 

113 xlabel='Frequency [Hz]', 

114 ylabel='$S_{V_H}$ [$\\mathrm{V}^2/\\mathrm{Hz}$]', 

115 ) 

116 ax.set_title(title) 

117 # save_plot('m506', 'png') 

118 

119 return ax 

120 

121 

122def show_PSD_classic(diff_voltages, title, ax=None, groupby_category='Vin', group_name=r'$\mu A$', 

123 num=10, style=[['science'], {'context': 'talk', 'style': 'white', 'palette': 'bright', }]): 

124 if not ax: 

125 fig, ax = plt.subplots(figsize=(12, 8)) 

126 c1 = sns.color_palette("hls", num) 

127 sns.set_palette(c1) 

128 # g = sns.relplot(x='freq', y='S', hue='Vin', data=diff_voltages, height=5, kind='line') 

129 grouped = diff_voltages.groupby(groupby_category) 

130 for group in grouped.groups.keys(): 

131 grouped.get_group(group).plot(x='freq', y='S', kind='line', 

132 loglog=True, ax=ax, 

133 label='%s %s' % (group * 1e-3, group_name), 

134 xlabel='Frequency [Hz]', 

135 ylabel='$S_{V_H}$ [$\\mathrm{V}^2/\\mathrm{Hz}$]', 

136 ) 

137 ax.set_title(title) 

138 return ax 

139 

140 

141def plot_PSD_contour(meas_obj, diff_voltages, title, 

142 cutoff_frequency=.9, 

143 groupby_category='Vin'): 

144 diff_voltages_contour = pd.DataFrame() 

145 for sr, m in meas_obj.items(): 

146 s = m.avg_spec 

147 s = s[s.freq < cutoff_frequency] 

148 if len(s) < 2: 

149 continue 

150 diff_voltages_contour[float(sr[:-2])] = s.S 

151 

152 v = diff_voltages[groupby_category].unique() 

153 v.sort() 

154 frequencies = diff_voltages.freq.unique() 

155 smin, smax = diff_voltages.S.min(), diff_voltages.S.max() 

156 levels = np.logspace(np.log10(smin), 

157 np.log10(smax), 10) 

158 

159 fig, ax = plt.subplots(figsize=(12, 8)) 

160 cs = ax.contourf(v, frequencies, diff_voltages_contour, 

161 norm=LogNorm(vmin=smin, vmax=smax), 

162 levels=levels, 

163 cmap=plt.cm.Blues, 

164 ) 

165 cbar = plt.gcf().colorbar(cs, ax=ax) 

166 cbar.set_label('$S_V^{V_{in}} (f)$') 

167 cbar.set_ticklabels(['%.1e' % _ for _ in levels]) 

168 ax.set_yscale('log') 

169 ax.set_ylabel('$f$ [Hz]') 

170 ax.set_xlabel('$V_{in}$ [$m$V]') 

171 ax.set_title(title) 

172 

173 

174def plot_diff_volt(nr=508, figsize=(16, 12)): 

175 datapath = glob(f'./data/MFN/m{nr}/*') 

176 meas_data, meas_info, all_data = load_data(datapath) 

177 meas_obj = calc_PSD(meas_data) 

178 diff_voltages = merge_data(meas_obj) 

179 

180 sns.set('talk', 'ticks') 

181 plt.rcParams['axes.grid'] = True 

182 plt.rcParams['grid.linestyle'] = '--' 

183 

184 fig, ax = plt.subplots(figsize=figsize) 

185 show_PSD_classic(diff_voltages, f'm{nr}: Compare different input current amplitudes', ax=ax) 

186 

187 inset2 = inset_axes(ax, width='100%', height='100%', 

188 bbox_to_anchor=(.54, .75, .3, .25), 

189 bbox_transform=ax.transAxes) 

190 inset3 = inset_axes(ax, width='100%', height='100%', 

191 bbox_to_anchor=(.1, .1, .3, .25), 

192 bbox_transform=ax.transAxes) 

193 

194 grouped = diff_voltages.groupby('Vin') 

195 for group in grouped.groups.keys(): 

196 g = grouped.get_group(group) 

197 fit_area = g.query('freq > %f and freq < %f' % (2e-2, 7e-1)) 

198 fit_area['lnf'] = np.log10(fit_area.freq) 

199 fit_area['lnS'] = np.log10(fit_area.S) 

200 fit = scipy.stats.linregress(fit_area.lnf, fit_area.lnS) 

201 intercept, slope = fit.intercept, -fit.slope 

202 voltage = group * 1e-3 

203 

204 inset2.plot(voltage, 10 ** intercept, 'o') 

205 inset3.plot(voltage, slope, 'o') 

206 

207 inset2.set_xlabel(r'Current [$\mu\mathrm{A}$]') 

208 inset2.set_ylabel('$S_{V_H} (f=1\\;$Hz$)$') 

209 inset2.set_yscale('log') 

210 

211 inset3.set_xlabel(r'Current [$\mu\mathrm{A}$]') 

212 inset3.set_ylabel('$\\alpha$') 

213 

214 

215def plot_diff_sweeprates(nr=507, figsize=(16, 12)): 

216 datapath = glob(f'./data/MFN/m{nr}/*') 

217 meas_data, meas_info, all_data = load_data(datapath, focus='SR', cut=None) 

218 meas_obj = calc_PSD(meas_data) 

219 f_max = (8 / (2 * np.pi)) 

220 diff_voltages = merge_data(meas_obj, cutoff_frequency=f_max, cut=None) 

221 title = 'm%s: Different Sweeprates ($\\tau = 100~\\mathrm{ms}$; $f_{Ref} = 17~\\mathrm{Hz}$)' % nr 

222 

223 sns.set('talk', 'ticks') 

224 plt.rcParams['axes.grid'] = True 

225 plt.rcParams['grid.linestyle'] = '--' 

226 

227 ax = plot_PSD_classic(diff_voltages, title, groupby_category='Vin') 

228 

229 inset2 = inset_axes(ax, width='100%', height='100%', 

230 bbox_to_anchor=(.15, .5, .3, .25), 

231 bbox_transform=ax.transAxes) 

232 inset3 = inset_axes(ax, width='100%', height='100%', 

233 bbox_to_anchor=(.08, .09, .3, .24), 

234 bbox_transform=ax.transAxes) 

235 

236 grouped = diff_voltages.groupby('Vin') 

237 for group in grouped.groups.keys(): 

238 g = grouped.get_group(group) 

239 fit_area = g.query('freq > %f and freq < %f' % (8e-2, 7e-1)) 

240 fit_area['lnf'] = np.log10(fit_area.freq) 

241 fit_area['lnS'] = np.log10(fit_area.S) 

242 fit = scipy.stats.linregress(fit_area.lnf, fit_area.lnS) 

243 intercept, slope = fit.intercept, -fit.slope 

244 

245 if group == 0: 

246 continue 

247 

248 voltage = group * 1e3 

249 

250 inset2.plot(voltage, 10 ** intercept, 'o') 

251 inset3.plot(voltage, slope, 'o') 

252 

253 inset2.set_xlabel('Sweeprate [$\\mathrm{mT}/\\mathrm{min}$]') 

254 inset2.set_ylabel('$S_{V_H} (f=1\\;$Hz$)$') 

255 inset2.set_yscale('log') 

256 

257 inset3.set_xlabel('Sweeprate [$\\mathrm{mT}/\\mathrm{min}$]') 

258 inset3.set_ylabel('$\\alpha$') 

259 

260 

261class PlotMasterCase(unittest.TestCase): 

262 def setUp(self): 

263 self.meas = {} 

264 

265 filename = 'output/' 

266 if not os.path.exists(filename): 

267 os.makedirs(filename) 

268 

269 # figures? 

270 self.save_figures = True 

271 self.ext = 'pdf' # or png, pgf 

272 self.figsize = 16, 12 

273 self.style = dict( 

274 default=True, 

275 figsize=self.figsize, 

276 latex=True, 

277 grid=True, 

278 set_mpl_params=True, 

279 ) 

280 

281 self.m = ana.Hloop(57) 

282 self.m.style.set_style(**self.style) 

283 self.plot = ana.plot.Plot() 

284 

285 def fit_data(self, mfn): 

286 df_fit = pd.DataFrame() 

287 for field in np.sort(mfn.data.Field.unique()): 

288 signal = mfn.data.query('Field == @field') 

289 mean = signal.Vx.mean() 

290 signal.loc[:, 'fittedVx'] = signal.Vx - mean 

291 

292 df_fit = pd.concat([df_fit, signal]) 

293 return df_fit 

294 

295 def split_data(self, mfn, df_fit): 

296 neg_saturation = mfn.data.Field.min() 

297 field_range = mfn.data.Field.max() - mfn.data.Field.min() 

298 steps = 4 

299 df = pd.DataFrame() 

300 for i in range(steps): 

301 stepsize = (field_range / steps) 

302 start = neg_saturation + i * stepsize 

303 end = start + stepsize 

304 df_tmp = df_fit.query('Field >= @start and Field < @end') 

305 df_tmp.reset_index(inplace=True) 

306 df_tmp['Group'] = i 

307 for j, f in enumerate(np.sort(df_tmp.Field.unique())): 

308 df_tmp.loc[df_tmp.query('Field == @f').index, 'SubGroup'] = j 

309 df_tmp.loc[df_tmp.query('Field == @f').index, 'x'] = df_tmp.loc[ 

310 df_tmp.query('Field == @f').index, 'fittedVx'] + i * .2 

311 df = pd.concat([df, df_tmp]) 

312 return df 

313 

314 def plot_hist(self, mfn, 

315 filename = '', 

316 factor = 1e3, 

317 steps = 4): 

318 # plt.style.use(['science']) 

319 sns.set('paper', 'white', rc={"axes.facecolor": (0, 0, 0, 0), 

320 'figure.figsize': (6, 4), 

321 }) 

322 neg_saturation = mfn.data.Field.min() 

323 field_range = mfn.data.Field.max() - mfn.data.Field.min() 

324 

325 df = self.split_data(mfn, self.fit_data(mfn)) 

326 

327 for i in range(steps): 

328 stepsize = (field_range / steps) 

329 start = neg_saturation + i * stepsize 

330 end = start + stepsize 

331 df_tmp = df.query('Field >= @start and Field < @end').copy() 

332 df_tmp.loc[:, 'fittedVx'] *= factor 

333 

334 # Initialize the FacetGrid object 

335 pal = sns.dark_palette(color="blue", n_colors=10, input='huls') 

336 g = sns.FacetGrid(df_tmp, row="Field", hue="Field", aspect=5, 

337 height=.8, palette=pal) 

338 g.map(sns.kdeplot, "fittedVx", clip_on=True, shade=True, alpha=1, 

339 lw=1.5, bw_adjust=.2) 

340 g.map(sns.kdeplot, "fittedVx", clip_on=True, color="w", lw=2, 

341 bw_adjust=.2) 

342 # g.map(sns.distplot, "fittedVx", hist=True, norm_hist=True, 

343 # hist_kws={'alpha': .5}) 

344 g.map(plt.axhline, y=0, lw=2, clip_on=True) # .set(xlim=(-1.3,1.3)) 

345 

346 def label(x, color, label): 

347 ax = plt.gca() 

348 ax.text(0, .3, label, fontweight="bold", color=color, 

349 ha="left", va="center", fontsize=16, 

350 transform=ax.transAxes) 

351 

352 g.map(label, "fittedVx") 

353 

354 # Set the subplots to overlap 

355 g.fig.subplots_adjust(hspace=-.3) 

356 

357 # Remove axes details that don't play well with overlap 

358 g.set_titles("") 

359 g.set(yticks=[], xlim=(-2.5, 2)) # , xticks=[-100 + 50*_ for _ in range(5)]) 

360 g.despine(bottom=True, left=True) 

361 unit = 'mV' if factor == 1e3 else '\\muV' if factor == 1e6 else 'V' 

362 plt.gca().set_xlabel('$V_H$ [$\\mathrm{%s}$]' % unit) 

363 if filename: 

364 plt.savefig('output/%s' % (filename.format(i))) 

365 

366 def test_plot_hist(self): 

367 mfn = ana.MFN(446) 

368 self.plot_hist(mfn, 'hist446_{}.pdf') 

369 

370 def test_plot_hist2(self): 

371 mfn = ana.MFN(447) 

372 self.plot_hist(mfn, 'hist447_{}.pdf') 

373 

374 def test_plot_compare_hyst(self): 

375 self.plot.plot_hloops(to_show=[23, 22, 119, 120], 

376 xlim=[-150, 150], 

377 filename='output/hysteresis_45') 

378 

379 def test_plot_compare_hyst_par(self): 

380 self.plot.plot_hloops_90(mirror=True, 

381 filename='output/hysteresis_90') 

382 

383 def test_plot_compare_multi_hyst_plusses(self): 

384 insets = [ 

385 ((.05, .05, .35, .3), 

386 ((-225, 0), (.1, .7), 'blue', .1)), 

387 ((.05, .64, .3, .35), 

388 ((-520, -200), (.05, .35), 'green', .1)), 

389 ((.59, .64, .25, .25), 

390 ((0, 200), (-.33, .05), 'red', .1)), 

391 ] 

392 fig, (ax, iax) = self.plot.plot_hloops_compare_90( 

393 to_show=[414 + _ for _ in range(4)], 

394 insets=insets, 

395 nograd=True 

396 ) 

397 plt.savefig('output/hloops_plusses.pdf') 

398 

399 def test_plot_compare_multi_hyst_crosses(self): 

400 alpha = .08 

401 insets = [ 

402 ((.05, .05, .33, .33), 

403 ((-200, -100), (-.9, -.7), 'red', alpha)), 

404 ((.05, .64, .31, .35), 

405 ((-70, -0), (-1.1, -.7), 'blue', alpha)), 

406 ((.42, .74, .31, .25), 

407 ((40, 140), (-.9, -.7), 'green', alpha)), 

408 ((.7, .06, .3, .3), 

409 ((160, 210), (-.77, -.6), 'yellow', alpha)), 

410 ] 

411 fig, (ax, iax) = self.plot.plot_hloops_compare_90( 

412 to_show=[414, 415, 416, 417], structure="crosses", 

413 ylim=(-1.26, -.3), xlim=(-300, 250), 

414 insets=insets, 

415 nograd=True, 

416 nodiff=True) 

417 plt.savefig('output/hloops_crosses.pdf') 

418 

419 def test_sa_static_fields(self): 

420 self.plot.plot_static_fields(filename='output/sr785-static') 

421 

422 def test_sa_first_sweeps(self): 

423 self.plot.first_sweeps(filename='output/sr785-first') 

424 

425 def test_sa_fspans(self): 

426 self.plot.multiple_fspans(filename='output/sr785-fspans') 

427 

428 def test_sa_sweeprates(self): 

429 self.plot.compare_sweeprates(filename='output/sr785-sweeprates') 

430 

431 def test_sa_voltages(self): 

432 self.plot.compare_voltages(filename='output/sr785-voltages') 

433 

434 def test_sa_temp(self): 

435 self.plot.sv_temp_effect(filename='output/sr785-temp') 

436 

437 def test_sa_fast(self): 

438 self.plot.fast_sweeprate(show_numbers=[320, 321, 325, 326, 

439 323, 324, 329, 322, 

440 327, 328, 336, 333, 

441 332, 331], 

442 filename='output/sr785-fast') 

443 

444 def plot_compare_time(self, m, fields, add=''): 

445 fig, axes, tax = m.plot_compare_timesignals(fields) 

446 for i, ax in enumerate(axes): 

447 ax.set_xlabel('Time [s]') 

448 tax[i].set_xlabel(r'\Large \# Counts') 

449 if i == 0: 

450 ax.set_ylabel('$V_x$ [V]') 

451 

452 plt.suptitle(r'Timesignals inside the Hysteresis %s' % add) 

453 

454 def test_time(self): 

455 for nr, fields in [(446, [-.04, .01]), 

456 (447, [.04, .015])]: 

457 mfn = ana.MFN(nr) 

458 mfn.style.set_style(default=True, grid=False, 

459 size='talk', style='ticks', 

460 palette='deep', latex=True) 

461 self.plot_compare_time(mfn, 

462 fields, 

463 '(m{})'.format(nr)) 

464 plt.savefig('output/daq-time-{}.pdf'.format(nr)) 

465 

466 def test_info446(self): 

467 m446 = ana.MFN(446) 

468 m446.style.set_style(figsize=(12,8), size='notebook') 

469 m446.plot_info(show_field=.1, gridspec_kw=dict(wspace=.3, 

470 hspace=.45)) 

471 

472 def test_info447(self): 

473 m447 = ana.MFN(447) 

474 m447.style.set_style(figsize=(12,8), size='notebook') 

475 m447.plot_info(show_field=.15, gridspec_kw=dict(wspace=.3, 

476 hspace=.45)) 

477 

478 def test_daq_sweeprates(self): 

479 plot_diff_sweeprates(nr=507, figsize=(16, 12)) 

480 plt.savefig('output/daq-sweeprate.pdf') 

481 

482 def test_daq_volt(self): 

483 plot_diff_volt(nr=506, figsize=(16, 12)) 

484 plt.savefig('output/daq-volt.pdf') 

485 

486if __name__ == '__main__': 

487 unittest.main()