.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/nddata_NNLS_Tikhonov_190104.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_nddata_NNLS_Tikhonov_190104.py: 2D ILT test =========== 2D regularization test .. GENERATED FROM PYTHON SOURCE LINES 7-343 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/images/sphx_glr_nddata_NNLS_Tikhonov_190104_001.png :alt: True $T_{1}$,$T_{2}$ :srcset: /auto_examples/images/sphx_glr_nddata_NNLS_Tikhonov_190104_001.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_nddata_NNLS_Tikhonov_190104_002.png :alt: 2D test data :srcset: /auto_examples/images/sphx_glr_nddata_NNLS_Tikhonov_190104_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_nddata_NNLS_Tikhonov_190104_003.png :alt: Solution :srcset: /auto_examples/images/sphx_glr_nddata_NNLS_Tikhonov_190104_003.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none ---------- logging output to /home/jmfranck/pyspecdata.0.log ---------- not in notebook, calling show | .. code-block:: Python try: get_ipython().magic('load_ext pyspecdata.ipy') in_notebook = True except: from pyspecdata import * in_notebook = False from pyspecdata import nnls_regularized from numpy import * import time fl=figlist_var() init_logging(level='debug') # got the following from here: # https://medium.com/pythonhive/python-decorator-to-measure-the-execution-time-of-methods-fa04cb6bb36d #l_line = '' #def timeit(method,n_times=5): # def timed(*args, **kw): # timing = zeros(n_times+1) # timing[0] = time.time() # for j in range(n_times): # result = method(*args, **kw) # timing[j+1] = time.time() # time_diff = (timing[-1]-timing[0])/float(n_times) # if 'log_time' in kw: # name = kw.get('log_name', method.__name__.upper()) # kw['log_time'][name] = int(time_diff * 1000) # else: # print('%r %2.2f ms (average of %d runs)' %(method.__name__, time_diff * 1000, n_times) + l_line) # return result # return timed # # ## # # #R = r_[1.:100:500j] # distribution of T2 relaxation rates #peaks = [(80,4,1),(20,0.5,0.5),(30,0.5,0.25)] #calcd = False #for mu,sigma,A in peaks: # if not calcd: # P = A*exp(-(R-mu)**2/(2*sigma**2)) # calcd = True # else: # P += A*exp(-(R-mu)**2/(2*sigma**2)) #P = nddata(P,'R') #print("your labels are",P.dimlabels) #P.setaxis('R',R.ravel()) # # ## Vary R as we move along the rows # #fl.next('distribution function') #fl.plot(P) # # ## # # #time_endpoint = 0.2 #t = nddata(r_[1e-3:time_endpoint:2048j],'t') # column vectors give functions of time #R = P.fromaxis('R') #test_data = exp(-R*t).dot(P) #logger.debug(strm('when constructing test_data, shape of the data is',ndshape(test_data),"len of axis_coords_error",len(test_data.axis_coords_error))) # ##test_data += random.normal(scale = 0.01,size=(2048,1)) #test_data.add_noise(0.01) #fl.next('test data function') #fl.plot(test_data) #xlim(-time_endpoint/10,time_endpoint) # # ## Do the basic NNLS fit # # # #logger.debug(strm('before nnls, shape of the data is',ndshape(test_data),"len of axis_coords_error",len(test_data.axis_coords_error))) #test_fit = test_data.C.nnls('t',R,lambda x,y: exp(-y*x)) #fl.next('fit an exponential distribution',legend=True) #fl.plot(test_data, alpha=0.5, label='test signal') #K = test_fit.get_prop('nnls_kernel').C ##note that order doesn't matter for the following dot (done by dimension name) #fl.plot(test_fit.C.dot(K), alpha=0.5, label='fit') #fl.next('what does the fit distribution look like?') #fl.plot(test_fit) # # ## Now add regularization # # #def L_curve(l,r_norm,x_norm, show_l=None, s=1e-3, # **kwargs): # """plot L-curve using # # Parameters # ========== # l: double # lambda values # r_norm: double # norm of the residual # x_norm: double # norm of solution vector""" # print(l.shape, r_norm.shape, x_norm.shape) # r_to_x = UnivariateSpline(log10(r_norm),log10(x_norm), s=s) # r_fine = linspace(*tuple(log10(r_norm[r_[0,-1]]).tolist()+[1000])) # plot(log10(r_norm),log10(x_norm),'o',**kwargs) # a = gca().get_ylim() # plot(r_fine,r_to_x(r_fine), alpha=0.3) # gca().set_ylim(a) # if show_l is not None: # l_to_r = UnivariateSpline(l,log10(r_norm),s=1e-3) # this_r = l_to_r(show_l) # this_x = r_to_x(this_r) # OLDplot(this_r,this_x,'ro') # annotate('chosen $\lambda=%3g$'%show_l, (this_r,this_x)) # for j,this_l in enumerate(l): # annotate('%5g'%this_l, (log10(r_norm[j]),log10(x_norm[j])), # ha='left',va='bottom',rotation=45) # ylabel('$\log_{10}(x$ norm$)$') # xlabel('$\log_{10}($ residual $)$') # # ## # # ##l = sqrt(logspace(-8,4,10)) # I do this because it gives me a fairly even spacing of points ###@timeit ##def nonvec_lcurve(l): ## x_norm = empty_like(l) ## r_norm = empty_like(l) ## for j,this_l in enumerate(l): ## x = test_data.C.nnls('t', ## R,lambda x,y: exp(-y*x), l=this_l) ## r_norm[j] = x.get_prop('nnls_residual') ## x_norm[j] = linalg.norm(x.data) ## return x,x_norm,r_norm ##x,x_norm,r_norm = nonvec_lcurve(l) ###x_norm = map(linalg.norm,x) # to be fair, this calculation is done outside the timing, below ## ## ### ## ## ##fl.next('L-curve', legend=True);clf() ##L_curve(l, r_norm, x_norm, markersize=10, alpha=0.5, label='manual loop') # # ## ## Vectorized version of lambda curve # ## # # #l = sqrt(logspace(-8,4,10)) # I do this because it gives me a fairly even spacing of points ##@timeit #def vec_lcurve(l): # return test_data.C.nnls('t', # R,lambda x,y: exp(-y*x), l=l) #x = vec_lcurve(l) # # ## # # #fl.next('L-curve') #logger.debug(strm("check dtype of residual:",x.get_prop('nnls_residual').data.dtype)) #L_curve(l, x.get_prop('nnls_residual').data, x.C.run(linalg.norm,'R').data, # markersize=5, alpha=0.5, label='compiled loop') # # ## ## 1.5D example -- very simplistic dataset ## ## (where the data is duplicated along the $\Omega$ dimension) ## ## unlike for the numpy version, I skip straight to the ## vectorized/parallel version. # # #l = sqrt(logspace(-8,4,10)) # I do this because it gives me a fairly even spacing of points #test_data_2d = test_data * nddata(r_[1,1,1],r'\Omega') ##@timeit #def multifreq_lcurve(l): # return test_data_2d.C.nnls('t', # R,lambda x,y: exp(-y*x), l=l) #x = multifreq_lcurve(l) # ## # #fl.next('L-curve') #chosen_l = 0.025 #L_curve(l, x.get_prop('nnls_residual')[r'\Omega',0].data, # x.C.run(linalg.norm,'R')[r'\Omega',0].data, markersize=5, alpha=0.5, # show_l=chosen_l, # label='compiled loop') # ## and show the final result ## here, I omit the SVD (allows negative) result # # #fl.next(r'show result where $\lambda$ set to knee') #test_data.nnls('t', R, lambda x,y: exp(-x*y), l=0.025) #fl.plot(test_data) #fl.plot(P) # # ## ## 1.5D code -- a more complicated/realistic example ## ## First, generate a test distribution, where we assume we have two dimensions: ## ## * A $t_2$ dimension, which is actually typically a frequency dimension ## * A $t_{indirect}$ dimension -- representing, *e.g.*, the delay of an inversion recovery experiment. # # #time_endpoint = 0.2 #R = nddata(r_[1.:100:100j],'R') # distribution of T2 relaxation rates #t2 = nddata(r_[1e-3:time_endpoint:2048j],'t2') # column vectors give functions of time #peaks_2D = [(20,1,0.8,0.025,5e-3,1), # (30,1,0.8,0.075,1e-2,1), # (80,5,1,0.175,1e-2,1)] #data_dist = 0 #for x_mu,x_sigma,x_amp,y_mu,y_sigma,y_amp in peaks_2D: # data_dist += x_amp*y_amp*exp(-(R-x_mu)**2/2/x_sigma**2 # -(t2-y_mu)**2/2/y_sigma**2) #fl.next('Test 1.5D true distribution') #fl.image(data_dist) # # ## Convert to the time domain, where the $R$ dimension is replaced by a ## $t_{indirect}$ dimension, and add noise # # #n_indirect = 128 #t_indirect = nddata((r_[0:n_indirect]+1)/double(n_indirect), # 't_indirect') #kernel = exp(-t_indirect*R) #test_data_2d = kernel.C.dot(data_dist) #fl.next('Test 1.5D dataset') #test_data_2d.add_noise(0.1) #fl.image(test_data_2d) # # ## Generate the L-curve # #l = sqrt(logspace(-8,4,10)) ##@timeit #def multifreq_lcurve(l): # return test_data_2d.C.nnls('t_indirect', # R,lambda x,y: exp(-y*x), l=l) #x = multifreq_lcurve(l) # # ## # # #print(ndshape(test_data_2d)) #print(ndshape(x)) # # ## # # #fl.next('L-curve') #heel_lambda = 0.215 #L_curve(l, x.get_prop('nnls_residual').C.sum('t2').data, # x.C.run(linalg.norm,'R').sum('t2').data, markersize=5, alpha=0.5, # label='compiled loop', # show_l=heel_lambda, s=1e-2) # # ## Final result for 1.5D test data # #fl.next(r'show result where $\lambda$ set to knee') #result = test_data_2d.C.nnls('t_indirect',R, # lambda x,y: exp(-y*x), l=heel_lambda) #fl.image(result) # # ## ## Testing 2D BRD extension def Gaussian_2D(x_axis,y_axis,mu_x,mu_y,sigma_x,sigma_y): this_Gaussian = exp(-(x_axis-mu_x)**2/2/sigma_x**2 -(y_axis-mu_y)**2/2/sigma_y**2) return this_Gaussian def Gaussian_2D_corr(theta,x_axis,y_axis,mu_x,mu_y,sigma_x,sigma_y): diff_xy = cos(theta)*x_axis + sin(theta)*y_axis sum_xy = -sin(theta)*x_axis + cos(theta)*y_axis diff_mu = cos(theta)*mu_x + sin(theta)*mu_y sum_mu = -sin(theta)*mu_x + cos(theta)*mu_y this_Gaussian = exp(-(diff_xy-diff_mu)**2/2/sigma_x**2 -(sum_xy-sum_mu)**2/2/sigma_y**2) return this_Gaussian Nx = 100 Ny = 50 x_min = 3e-3; x_max = 3; y_min = 3e-3; y_max = 2; x_axis_log = nddata(linspace(log10(x_min),log10(x_max),Nx),r'log(T1)') y_axis_log = nddata(linspace(log10(y_min),log10(y_max),Ny),r'log(T2)') fl.next(r'True $T_{1}$,$T_{2}$') dist = Gaussian_2D(x_axis_log,y_axis_log,-1.6,-0.5,0.15,0.15) dist += Gaussian_2D_corr(-195*pi/180,x_axis_log,y_axis_log,-0.1,0.01,0.35,0.15) dist += Gaussian_2D_corr(-115*pi/180,x_axis_log,y_axis_log,-1.8,-1.5,0.09,0.15) fl.image(dist) N_tau1 = 30 N_tau2 = 1000 tau1_min = 5e-4; tau1_max = 4; tau2_min = 5e-4; tau2_max = 1.2 tau1_axis = nddata(logspace(log10(tau1_min),log10(tau1_max),N_tau1),'tau1') tau2_axis = nddata(linspace(tau2_min,tau2_max,N_tau2),'tau2') s = dist*exp(-tau2_axis/(10**y_axis_log))*(1.-2*exp(-tau1_axis/(10**x_axis_log))) s.sum('log(T1)').sum('log(T2)') s /= amax(s.data) s.add_noise(0.03) s.reorder('tau1') fl.next('2D test data') fl.image(s) x_axis = nddata(10**x_axis_log.data.real,'T1') y_axis = nddata(10**y_axis_log.data.real,'T2') solution = s.C.nnls(('tau1','tau2'), (x_axis,y_axis), (lambda x1,x2: 1.-2*exp(-x1/x2), lambda y1,y2: exp(-y1/y2)), l='BRD') fl.next('Solution') fl.image(solution) # if in_notebook: print("in notebook, not calling show") else: print("not in notebook, calling show") fl.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (31 minutes 21.596 seconds) .. _sphx_glr_download_auto_examples_nddata_NNLS_Tikhonov_190104.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: nddata_NNLS_Tikhonov_190104.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: nddata_NNLS_Tikhonov_190104.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_