results.py 25 KB
Newer Older
1
2
3
4
5
# results.py - Collate results from MRS fits
#
# Author: Will Clarke <william.clarke@ndcn.ox.ac.uk>
#         Saad Jbabdi <saad@fmrib.ox.ac.uk>
#
6
# Copyright (C) 2020 University of Oxford
7
8
9
10
11
12
# SHBASECOPYRIGHT

import pandas as pd
import fsl_mrs.utils.models as models
import fsl_mrs.utils.quantify as quant
import fsl_mrs.utils.qc as qc
13
from fsl_mrs.utils.misc import FIDToSpec, SpecToFID, calculate_lap_cov
14

15
import numpy as np
16
from copy import deepcopy
17

18

19
20
21
22
class FitRes(object):
    """
       Collects fitting results
    """
23
24
25

    def __init__(self, model, method, names, metab_groups, baseline_order, B, ppmlim):
        """ Short class initilisation """
26
        # Initilise some basic parameters - includes fitting options
William Clarke's avatar
William Clarke committed
27
        # Populate parameter names
28
29
30
31
        self.model = model
        self.fill_names(names, baseline_order=baseline_order, metab_groups=metab_groups)
        self.method = method
        self.ppmlim = ppmlim
32
        self.baseline_order = baseline_order
33
        self.base_poly = B
34

35
        self.concScalings = {'internal': None, 'internalRef': None, 'molarity': None, 'molality': None, 'info': None}
36

37
    def loadResults(self, mrs, fitResults):
38
39
        "Load fitting results and calculate some metrics"
        # Populate data frame
40
41
        if fitResults.ndim == 1:
            self.fitResults = pd.DataFrame(data=fitResults[np.newaxis, :], columns=self.params_names)
42
43
44
        else:
            self.fitResults = pd.DataFrame(data=fitResults, columns=self.params_names)
        self.params = self.fitResults.mean().values
45
46
47
48

        # Store prediction, baseline, residual
        self.pred = self.predictedFID(mrs, mode='Full')
        self.baseline = self.predictedFID(mrs, mode='Baseline')
49
50
51
        self.residuals = mrs.FID - self.pred

        # Calculate single point crlb and cov
52
53
54
55
56
57
58
59
60
61
        first, last = mrs.ppmlim_to_range(self.ppmlim)
        _, _, forward, _, _ = models.getModelFunctions(self.model)

        def forward_lim(p):
            return forward(p, mrs.frequencyAxis,
                           mrs.timeAxis,
                           mrs.basis,
                           self.base_poly,
                           self.metab_groups,
                           self.g)[first:last]
62
        data = mrs.get_spec(ppmlim=self.ppmlim)
63
64
65
66
67
68
        # self.crlb      = calculate_crlb(self.params,forward_lim,data)
        self.cov = calculate_lap_cov(self.params, forward_lim, data)
        self.crlb = np.diagonal(self.cov)
        std = np.sqrt(self.crlb)
        self.corr = self.cov / (std[:, np.newaxis] * std[np.newaxis, :])
        self.mse = np.mean(np.abs(FIDToSpec(self.residuals)[first:last])**2)
69
70

        with np.errstate(divide='ignore', invalid='ignore'):
71
72
            self.perc_SD = np.sqrt(self.crlb) / self.params * 100
        self.perc_SD[self.perc_SD > 999] = 999   # Like LCModel :)
73
74
75
76
77
78
79
        self.perc_SD[np.isnan(self.perc_SD)] = 999

        # Calculate mcmc metrics
        if self.method == 'MH':
            self.mcmc_cov = self.fitResults.cov().values
            self.mcmc_cor = self.fitResults.corr().values
            self.mcmc_var = self.fitResults.var().values
Saad Jbabdi's avatar
Saad Jbabdi committed
80
            self.mcmc_samples = self.fitResults.values
81

Saad Jbabdi's avatar
Saad Jbabdi committed
82
83
84
85
        # VB metrics
        if self.method == 'VB':
            self.vb_cov = self.optim_out.cov
            self.vb_var = self.optim_out.var
86
            std = np.sqrt(self.vb_var)
87
            self.vb_corr = self.vb_cov / (std[:, np.newaxis] * std[np.newaxis, :])
Saad Jbabdi's avatar
Saad Jbabdi committed
88

89
        self.hzperppm = mrs.centralFrequency / 1E6
90
91

        # Calculate QC metrics
92
        self.FWHM, self.SNR = qc.calcQC(mrs, self, ppmlim=(0.2, 4.2))
93

94
95
        # Run relative concentration scaling to tCr in 'default' 1H MRS case.
        # Create combined metab at same time to avoid later errors.
96
        if (('Cr' in self.metabs) and ('PCr' in self.metabs)):
97
            self.combine([['Cr', 'PCr']])
98
            self.calculateConcScaling(mrs)
99

100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
    def calculateConcScaling(self,
                             mrs,
                             quant_info=None,
                             internal_reference=['Cr', 'PCr'],
                             verbose=False):
        """Run calculation of internal and (if possible) water concentration scaling.

        :param mrs: MRS object
        :type mrs: fsl_mrs.core.mrs.MRS
        :param quant_info: Quantification information for water scaling, defaults to None
        :type quant_info: fsl_mrs.utils.quantify.QuantificationInfo, optional
        :param internal_reference: Internal referencing emtabolite, defaults to ['Cr', 'PCr'] i.e. tCr
        :type internal_reference: list, optional
        :param verbose: Enable for verbose output, defaults to False
        :type verbose: bool, optional
        """

        self.intrefstr = '+'.join(internal_reference)
        self.referenceMetab = internal_reference

        internalRefScaling = quant.quantifyInternal(internal_reference, self.getConc(), self.metabs)

        if mrs.H2O is not None and quant_info is not None:
            molalityScaling, molarityScaling, ref_info = quant.quantifyWater(mrs,
                                                                             self,
                                                                             quant_info,
                                                                             verbose=verbose)
127
128
129
130
131
132

            self.concScalings = {
                'internal': internalRefScaling,
                'internalRef': self.intrefstr,
                'molarity': molarityScaling,
                'molality': molalityScaling,
133
134
                'quant_info': quant_info,
                'ref_info': ref_info}
135
        else:
136
137
138
139
140
            self.concScalings = {
                'internal': internalRefScaling,
                'internalRef': self.intrefstr,
                'molarity': None,
                'molality': None,
141
142
                'quant_info': None,
                'ref_info': None}
143
144

    def combine(self, combineList):
145
146
147
        """Combine two or more basis into single result"""
        # Create extra entries in the fitResults DF , add to param_names and recalculate
        for toComb in combineList:
148
149
            newstr = '+'.join(toComb)
            if newstr in self.metabs:
150
151
152
153
                continue
            ds = pd.Series(np.zeros(self.fitResults.shape[0]), index=self.fitResults.index)
            jac = np.zeros(self.cov.shape[0])
            for metab in toComb:
154
155
156
157
                if metab not in self.metabs:
                    raise ValueError(f'Metabolites to combine must be in res.metabs. {metab} not found.')
                ds = ds.add(self.fitResults[metab])
                index = self.metabs.index(metab)
158
                jac[index] = 1.0
159
160

            self.fitResults[newstr] = pd.Series(ds, index=self.fitResults.index)
William Clarke's avatar
William Clarke committed
161
            self.params_names_inc_comb.append(newstr)
162
163
164
            self.metabs.append(newstr)
            newCRLB = jac @ self.cov @ jac
            self.crlb = np.concatenate((self.crlb, newCRLB[np.newaxis]))
165

166
        self.numMetabs = len(self.metabs)
167
168
169
        # self.params = self.fitResults.mean().values
        with np.errstate(divide='ignore', invalid='ignore'):
            params = self.fitResults.mean().values
170
171
            self.perc_SD = np.sqrt(self.crlb) / params * 100
        self.perc_SD[self.perc_SD > 999] = 999   # Like LCModel :)
172
        self.perc_SD[np.isnan(self.perc_SD)] = 999
173

174
175
176
177
        if self.method == 'MH':
            self.mcmc_cov = self.fitResults.cov().values
            self.mcmc_cor = self.fitResults.corr().values
            self.mcmc_var = self.fitResults.var().values
178

179
180
    def predictedSpec(self, mrs, mode='Full', noBaseline=False, ppmlim=None, shift=True):
        if ppmlim is None:
William Clarke's avatar
William Clarke committed
181
            ppmlim = self.ppmlim
182

William Clarke's avatar
William Clarke committed
183
        if mode.lower() == 'full':
184
185
186
187
188
189
190
            out = models.getFittedModel(
                self.model,
                self.params,
                self.base_poly,
                self.metab_groups,
                mrs,
                noBaseline=noBaseline)
William Clarke's avatar
William Clarke committed
191
        elif mode.lower() == 'baseline':
192
193
            out = models.getFittedModel(self.model, self.params, self.base_poly,
                                        self.metab_groups, mrs, baselineOnly=True)
William Clarke's avatar
William Clarke committed
194
        elif mode in self.metabs:
195
196
197
198
199
200
201
202
203
            out = models.getFittedModel(
                self.model,
                self.params,
                self.base_poly,
                self.metab_groups,
                mrs,
                basisSelect=mode,
                baselineOnly=False,
                noBaseline=noBaseline)
William Clarke's avatar
William Clarke committed
204
205
206
        else:
            raise ValueError('Unknown mode, must be one of: Full, baseline or a metabolite name.')

207
        first, last = mrs.ppmlim_to_range(ppmlim=ppmlim, shift=shift)
William Clarke's avatar
William Clarke committed
208

209
210
        return out[first:last]

211
    def predictedFID(self, mrs, mode='Full', noBaseline=False, no_phase=False):
212
        if mode.lower() == 'full':
213
214
215
216
217
218
            out = models.getFittedModel(
                self.model,
                self.params,
                self.base_poly,
                self.metab_groups,
                mrs,
219
220
                noBaseline=noBaseline,
                no_phase=no_phase)
221
        elif mode.lower() == 'baseline':
222
223
224
225
226
227
228
229
            out = models.getFittedModel(
                self.model,
                self.params,
                self.base_poly,
                self.metab_groups,
                mrs,
                baselineOnly=True,
                no_phase=no_phase)
230
        elif mode in self.metabs:
231
232
233
234
235
236
237
238
            out = models.getFittedModel(
                self.model,
                self.params,
                self.base_poly,
                self.metab_groups,
                mrs,
                basisSelect=mode,
                baselineOnly=False,
239
240
                noBaseline=noBaseline,
                no_phase=no_phase)
241
242
243
244
245
        else:
            raise ValueError('Unknown mode, must be one of: Full, baseline or a metabolite name.')
        return SpecToFID(out)

    def __str__(self):
246
        out = "----- Fitting Results ----\n"
247
248
249
250
        out += " names          = {}\n".format(self.params_names)
        out += " params         = {}\n".format(self.params)
        out += " CRLB           = {}\n".format(self.crlb)
        out += " MSE            = {}\n".format(self.mse)
251
        # out += " cov            = {}\n".format(self.cov)
252
253
254
255
256
        out += " phi0 (deg)     = {}\n".format(self.getPhaseParams()[0])
        out += " phi1 (deg/ppm) = {}\n".format(self.getPhaseParams(phi1='deg_per_ppm')[1])
        out += " gamma (Hz)     = {}\n".format(self.getLineShapeParams())
        out += " eps (ppm)      = {}\n".format(self.getShiftParams())
        out += " b_norm         = {}\n".format(self.getBaselineParams())
257
258

        return out
259
260

    def fill_names(self, names, baseline_order=0, metab_groups=None):
261
262
263
264
265
        """
        mrs            : MRS Object
        baseline_order : int
        metab_groups   : list (by default assumes single metab group)
        """
266
267
268
        self.metabs = deepcopy(names)
        self.original_metabs = deepcopy(names)
        self.numMetabs = len(self.metabs)
269

270
271
272
273
274
275
        self.params_names = []
        self.params_names.extend(names)

        if metab_groups is None:
            g = 1
        else:
276
            g = max(metab_groups) + 1
277
278
279

        self.g = g
        self.metab_groups = metab_groups
280

281
282
        for i in range(g):
            self.params_names.append(f'gamma_{i}')
William Clarke's avatar
William Clarke committed
283
284
285
286
        if self.model.lower() == 'voigt':
            for i in range(g):
                self.params_names.append(f'sigma_{i}')

287
288
289
        for i in range(g):
            self.params_names.append(f"eps_{i}")

290
291
292
        self.params_names.extend(['Phi0', 'Phi1'])

        for i in range(0, baseline_order + 1):
293
294
295
            self.params_names.append(f"B_real_{i}")
            self.params_names.append(f"B_imag_{i}")

William Clarke's avatar
William Clarke committed
296
        self.params_names_inc_comb = deepcopy(self.params_names)
297

298
    def to_file(self, filename, what='concentrations'):
299
        """
300
        Save results to csv file
301
302
303
304

        Parameters:
        -----------
        filename : str
305
        what     : one of 'summary', 'concentrations, 'qc', 'parameters', 'concentrations-mh','parameters-mh'
306
307
        """

308
        if what == 'summary':
309
310
            df = pd.DataFrame()
            df['Metab'] = self.metabs
311

312
313
            if self.concScalings['internal'] is not None:
                concstr = f'/{self.concScalings["internalRef"]}'
314
                df[concstr] = self.getConc(scaling='internal')
315
                df[concstr + ' CRLB'] = self.getUncertainties(type='internal')
316
            else:
317
                df['Raw conc'] = self.getConc()
318
                df['Raw CRLB'] = self.getUncertainties(type='raw')
319
            if self.concScalings['molality'] is not None:
320
                df['mMol/kg'] = self.getConc(scaling='molality')
321
                df['mMol/kg CRLB'] = self.getUncertainties(type='molality')
322
            if self.concScalings['molarity'] is not None:
323
                df['mM'] = self.getConc(scaling='molarity')
324
                df['mM CRLB'] = self.getUncertainties(type='molarity')
325
326
            df['%CRLB'] = self.getUncertainties()

327
            SNR = self.getQCParams()[0]
328
            SNR.index = SNR.index.str.replace('SNR_', '')
329
330
            SNR.index.name = 'Metab'
            SNR = SNR.reset_index(name='SNR')
331
            df = df.merge(SNR, how='outer')
332
333

            FWHM = self.getQCParams()[1]
334
            FWHM.index = FWHM.index.str.replace('fwhm_', '')
335
336
            FWHM.index.name = 'Metab'
            FWHM = FWHM.reset_index(name='FWHM')
337
            df = df.merge(FWHM, how='outer')
338
339

        elif what == 'concentrations':
340
            scaling_type = ['raw']
341
342
343
344
345
346
347
            if self.concScalings['internal'] is not None:
                scaling_type.append('internal')
            if self.concScalings['molality'] is not None:
                scaling_type.append('molality')
            if self.concScalings['molarity'] is not None:
                scaling_type.append('molarity')

348
349
            std = [self.getUncertainties(type=st) for st in scaling_type]
            mean = [self.getConc(scaling=st, function='mean').T for st in scaling_type]
Saad Jbabdi's avatar
Saad Jbabdi committed
350
351

            d = {}
352
353
            d['mean'] = {key: value for key, value in zip(scaling_type, mean)}
            d['std'] = {key: value for key, value in zip(scaling_type, std)}
Saad Jbabdi's avatar
Saad Jbabdi committed
354

355
            dict_of_df = {k: pd.DataFrame(v) for k, v in d.items()}
Saad Jbabdi's avatar
Saad Jbabdi committed
356
            df = pd.concat(dict_of_df, axis=1)
357
            df.insert(0, 'Metabs', self.metabs)
358

359
        elif what == 'concentrations-mh':
360
            scaling_type = ['raw']
361
362
363
364
365
366
367
368
369
            if self.concScalings['internal'] is not None:
                scaling_type.append('internal')
            if self.concScalings['molality'] is not None:
                scaling_type.append('molality')
            if self.concScalings['molarity'] is not None:
                scaling_type.append('molarity')

            all_df = []
            for st in scaling_type:
370
371
                df = self.getConc(scaling=st, function=None).T
                df.insert(0, 'scaling', st)
372
373
374
375
376
377
                df.index.name = 'metabolite'
                all_df.append(df)

            df = pd.concat(all_df)
            df.reset_index(inplace=True)

Saad Jbabdi's avatar
Saad Jbabdi committed
378
379
380
        elif what == 'qc':

            SNR = self.getQCParams()[0]
381
            SNR.index = SNR.index.str.replace('SNR_', '')
Saad Jbabdi's avatar
Saad Jbabdi committed
382
383
            SNR.index.name = 'Metab'
            SNR = SNR.reset_index(name='SNR')
384

Saad Jbabdi's avatar
Saad Jbabdi committed
385
            FWHM = self.getQCParams()[1]
386
            FWHM.index = FWHM.index.str.replace('fwhm_', '')
Saad Jbabdi's avatar
Saad Jbabdi committed
387
388
            FWHM.index.name = 'Metab'
            FWHM = FWHM.reset_index(name='FWHM')
389

390
            df = SNR.merge(FWHM, how='outer')
391

392
        elif what == 'parameters':
Saad Jbabdi's avatar
Saad Jbabdi committed
393
            if self.method == 'MH':
394
                mean = self.fitResults.T.mean(axis=1)
395
396
397
                std = self.fitResults.T.std(axis=1)
                mean = mean[~mean.index.str.contains('\\+')].to_numpy()
                std = std[~std.index.str.contains('\\+')].to_numpy()
Saad Jbabdi's avatar
Saad Jbabdi committed
398
            else:
399
                mean = self.params
400
                std = np.sqrt(np.diagonal(self.cov))
Saad Jbabdi's avatar
Saad Jbabdi committed
401
402
            df = pd.DataFrame()
            df['parameter'] = self.params_names
403
404
            df['mean'] = mean
            df['std'] = std
405
406
407
408
409

        elif what == 'parameters-mh':
            df = self.fitResults.T
            df.index.name = 'parameter'
            df.reset_index(inplace=True)
410

411
        df.to_csv(filename, index=False, header=True)
412
413

    # Functions to return physically meaningful units from the fitting results
414
    def getConc(self, scaling='raw', metab=None, function='mean'):
415
        if function is None:
416
417
            def dfFunc(m):
                return self.fitResults[m]
418
        elif metab is None:
419
420
            def dfFunc(m):
                return self.fitResults[m].apply(function).to_numpy()
421
        else:
422
423
424
            def dfFunc(m):
                return self.fitResults[m].apply(function)

425
426
        # Extract concentrations from parameters.
        if metab is not None:
427
            if isinstance(metab, list):
428
429
430
431
432
433
                for mm in metab:
                    if mm not in self.metabs:
                        raise ValueError(f'{mm} is not a recognised metabolite.')
            else:
                if metab not in self.metabs:
                    raise ValueError(f'{metab} is not a recognised metabolite.')
William Clarke's avatar
William Clarke committed
434
            rawConc = dfFunc(metab)
435
        else:
William Clarke's avatar
William Clarke committed
436
            rawConc = dfFunc(self.metabs)
437

438
        if scaling == 'raw':
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
            return rawConc
        elif scaling == 'internal':
            if self.concScalings['internal'] is None:
                raise ValueError('Internal concetration scaling not calculated, run calculateConcScaling method.')
            return rawConc * self.concScalings['internal']

        elif scaling == 'molality':
            if self.concScalings['molality'] is None:
                raise ValueError('Molality concetration scaling not calculated, run calculateConcScaling method.')
            return rawConc * self.concScalings['molality']

        elif scaling == 'molarity':
            if self.concScalings['molarity'] is None:
                raise ValueError('Molarity concetration scaling not calculated, run calculateConcScaling method.')
            return rawConc * self.concScalings['molarity']
        else:
            raise ValueError(f'Unrecognised scaling value {scaling}.')

457
    def getPhaseParams(self, phi0='degrees', phi1='seconds', function='mean'):
458
459
460
461
462
463
        """Return the two phase parameters in specified units"""
        if function is None:
            p0 = self.fitResults['Phi0'].values
            p1 = self.fitResults['Phi1'].values
        else:
            p0 = self.fitResults['Phi0'].apply(function)
464
            p1 = self.fitResults['Phi1'].apply(function)
465
466

        if phi0.lower() == 'degrees':
467
            p0 *= np.pi / 180.0
468
        elif (phi0.lower() == 'radians') or (phi0.lower() == 'raw'):
469
            pass
470
471
        else:
            raise ValueError('phi0 must be degrees or radians')
472

473
        if phi1.lower() == 'seconds':
474
            p1 *= -1.0 / (2 * np.pi)
475
        elif phi1.lower() == 'deg_per_ppm':
476
            p1 *= -180.0 / np.pi * self.hzperppm
477
        elif phi1.lower() == 'deg_per_hz':
478
            p1 *= -180.0 / np.pi * 1.0
479
480
        elif phi1.lower() == 'raw':
            pass
481
        else:
482
            raise ValueError('phi1 must be seconds, deg_per_ppm, deg_per_hz or raw.')
483

484
        return p0, p1
485

486
    def getShiftParams(self, units='ppm', function='mean'):
487
488
        """ Return shift parameters (eps) in specified units - default = ppm."""
        if function is None:
489
            shiftParams = np.zeros([self.fitResults.shape[0], self.g])
490
            for g in range(self.g):
491
                shiftParams[:, g] = self.fitResults[f'eps_{g}'].values
492
493
494
495
        else:
            shiftParams = []
            for g in range(self.g):
                shiftParams.append(self.fitResults[f'eps_{g}'].apply(function))
496
            shiftParams = np.asarray(shiftParams)
497
498

        if units.lower() == 'ppm':
499
            shiftParams *= 1 / (2 * np.pi * self.hzperppm)
500
        elif units.lower() == 'hz':
501
            shiftParams *= 1 / (2 * np.pi)
502
503
        elif units.lower() == 'raw':
            pass
504
        else:
505
            raise ValueError('Units must be Hz, ppm or raw.')
506
507
508

        return shiftParams

509
    def getLineShapeParams(self, units='Hz', function='mean'):
510
        """ Return line broadening parameters (gamma and sigma) in specified units - default = Hz."""
511
        if self.model == 'lorentzian':
512
            if function is None:
513
                gamma = np.zeros([self.fitResults.shape[0], self.g])
514
                for g in range(self.g):
515
                    gamma[:, g] = self.fitResults[f'gamma_{g}'].values
516
517
518
519
            else:
                gamma = []
                for g in range(self.g):
                    gamma.append(self.fitResults[f'gamma_{g}'].apply(function))
520
                gamma = np.asarray(gamma)
521

522
523
            if units.lower() == 'hz':
                gamma *= 1 / (np.pi)
524
            elif units.lower() == 'ppm':
525
                gamma *= 1 / (np.pi * self.hzperppm)
526
527
            elif units.lower() == 'raw':
                pass
528
            else:
529
                raise ValueError('Units must be Hz, ppm or raw.')
William Clarke's avatar
William Clarke committed
530
            combined = gamma
531
532
            return combined, gamma
        elif self.model == 'voigt':
William Clarke's avatar
William Clarke committed
533
            if function is None:
534
535
                gamma = np.zeros([self.fitResults.shape[0], self.g])
                sigma = np.zeros([self.fitResults.shape[0], self.g])
William Clarke's avatar
William Clarke committed
536
                for g in range(self.g):
537
538
                    gamma[:, g] = self.fitResults[f'gamma_{g}'].values
                    sigma[:, g] = self.fitResults[f'sigma_{g}'].values
William Clarke's avatar
William Clarke committed
539
540
541
542
543
544
            else:
                gamma = []
                sigma = []
                for g in range(self.g):
                    gamma.append(self.fitResults[f'gamma_{g}'].apply(function))
                    sigma.append(self.fitResults[f'sigma_{g}'].apply(function))
545
                gamma = np.asarray(gamma)
William Clarke's avatar
William Clarke committed
546
547
                sigma = np.asarray(sigma)

548
549
550
            if units.lower() == 'hz':
                gamma *= 1 / (np.pi)
                sigma *= 1 / (np.pi)
William Clarke's avatar
William Clarke committed
551
            elif units.lower() == 'ppm':
552
553
                gamma *= 1 / (np.pi * self.hzperppm)
                sigma *= 1 / (np.pi * self.hzperppm)
554
555
            elif units.lower() == 'raw':
                pass
William Clarke's avatar
William Clarke committed
556
            else:
557
                raise ValueError('Units must be Hz, ppm or raw.')
558

559
560
561
562
            combined = gamma / 2 + np.sqrt((gamma**2 / 4.0) + (sigma * 2 * np.log(2))**2)
            return combined, gamma, sigma

    def getBaselineParams(self, complex=True, normalise=True):
563
564
        """ Return normalised complex baseline parameters."""
        bParams = []
565
        for b in range(self.baseline_order + 1):
566
567
568
            breal = self.fitResults[f'B_real_{b}'].mean()
            bimag = self.fitResults[f'B_imag_{b}'].mean()
            if complex:
569
                bParams.append(breal + 1j * bimag)
570
            else:
571
572
                bParams.extend([breal, bimag])

573
        bParams = np.asarray(bParams)
574
575
        if normalise:
            with np.errstate(divide='ignore', invalid='ignore'):
576
                return bParams / np.abs(bParams[0])
577
578
        else:
            return bParams
579

580
    def getQCParams(self, metab=None):
581
        """Returns peak wise SNR and FWHM (in Hz)"""
582
        if metab is None:
583
            return self.SNR.peaks.mean(), self.FWHM.mean()
584
        else:
585
            return self.SNR.peaks['SNR_' + metab].mean(), self.FWHM['fwhm_' + metab].mean()
586

587
588
589
    def getUncertainties(self, type='percentage', metab=None):
        """ Return the uncertainties (SD) on concentrations.
        Can either be in raw, molarity or molality or percentage uncertainties.
William Clarke's avatar
William Clarke committed
590
591
592

        """
        abs_std = []
593
594
        if metab is None:
            metab = self.metabs
595
596
        elif isinstance(metab, str):
            metab = [metab, ]
597
        for m in metab:
William Clarke's avatar
William Clarke committed
598
            if self.method == 'Newton':
599
                index = self.params_names_inc_comb.index(m)
William Clarke's avatar
William Clarke committed
600
601
                abs_std.append(np.sqrt(self.crlb[index]))
            elif self.method == 'MH':
602
                abs_std.append(self.fitResults[m].std())
William Clarke's avatar
William Clarke committed
603
            elif self.method == 'VB':
604
                index = self.params_names_inc_comb.index(m)
William Clarke's avatar
William Clarke committed
605
606
607
                abs_std.append(np.sqrt(self.vb_var[index]))
        abs_std = np.asarray(abs_std)
        if type.lower() == 'raw':
608
            return abs_std
William Clarke's avatar
William Clarke committed
609
610
611
612
613
        elif type.lower() == 'molarity':
            return abs_std * self.concScalings['molarity']
        elif type.lower() == 'molality':
            return abs_std * self.concScalings['molality']
        elif type.lower() == 'internal':
614
615
            internal_ref = self.concScalings['internalRef']
            if self.method == 'Newton':
616
                internalRefIndex = self.params_names_inc_comb.index(internal_ref)
617
618
619
620
621
622
                internalRefSD = np.sqrt(self.crlb[internalRefIndex])
            elif self.method == 'MH':
                internalRefSD = self.fitResults[internal_ref].std()
            elif self.method == 'VB':
                index = self.params_names_inc_comb.index(internal_ref)
                internalRefSD = np.sqrt(self.vb_var[index])
623
            abs_std = np.sqrt(abs_std**2 + internalRefSD**2)
William Clarke's avatar
William Clarke committed
624
625
            return abs_std * self.concScalings['internal']
        elif type.lower() == 'percentage':
626
            vals = self.fitResults[metab].mean().to_numpy()
627
628
            perc_SD = abs_std / vals * 100
            perc_SD[perc_SD > 999] = 999   # Like LCModel :)
William Clarke's avatar
William Clarke committed
629
630
631
632
            perc_SD[np.isnan(perc_SD)] = 999
            return perc_SD
        else:
            raise ValueError('type must either be "absolute" or "percentage".')