optimize computing speed of array processing

There is one bottleneck with computation of topographs.

First I read a bunch of .tif or .h5 files (sequential topographs at different positions on the crystal rocking curve). Then I combine all these in 3D array, data[k,j,i] where k index represents the snapshot number (angle on the rocking curve) and [i,j] is the XY location on the topograph. Then I do something straightforward: run a double "for" loop to calculate rocking curve statistics for each [i,j] and this takes time.

It becomes quite time consuming with ~1000x1000 pixels.
Is there a way to speed up the python code?

This main loop in topography processing could be optimized for execution time, perhaps by using optimized routines from numpy, scipy, or other numerical analysis packages. To optimize, we may need to divide this long code block so as to identify routines already implemented in numpy, scipy, or other.

dtxrd/rctopo

Lines 448 to 658 in 1b2708c

    
             for j in yrange: 
        
               for i in xrange: 
        
                 rcurve=data[:,j,i] 
        
                 #------------------------------------------------- 
        
                 # Background determination 
        
                 #------------------------------------------------- 
        
                 if opts.bkg==-1:     
        
                   bkg=0.25*(rcurve[0]+rcurve[1]+rcurve[len(rcurve)-2]+rcurve[len(rcurve)-1]) #; print 'bkg = ', bkg   
        
                   b1=0.5*(rcurve[0]+rcurve[1]) 
        
                   b2=0.5*(rcurve[len(rcurve)-2]+rcurve[len(rcurve)-1]) 
        
                   t1=th_beg 
        
                   t2=th_end 
        
                   A=(b2-b1)/(t2-t1) 
        
                   B=b2-A*t2 
        
                   bkg_lin=A*th+B 
        
                 #------------------------------------------------ 
        
                 # Deglitching procedure 
        
                 #------------------------------------------------                   
        
                 if df>1.0: 
        
                   index=range(len(rcurve)) 
        
                   index=index[2:len(rcurve)-2] 
        
                   for k in index: 
        
                     y=rcurve[k] 
        
                     mean=0.25*(rcurve[k-2]+rcurve[k-1]+rcurve[k+1]+rcurve[k+2]) 
        
                     if y > df*mean or y < (2.0-df)*mean: 
        
                       rcurve[k]=mean 
        
                       print("X= ", dx*i, "Y= ",dy*j, " glitch found")    
        
                 #--------------------------------------------------------------------------------------- 
        
                 # subtracting bkg since it can vary from pixel to pixel             
        
                 if opts.conduct == 1: 
        
                     #rcurve_tot = rcurve_tot + [-1.0*(rcurve-bkg_lin)]  # TC         
        
                     rcurve = -1.0*(rcurve-bkg_lin)  # TC 
        
                     bkgs = 0.0    # TC       
        
                     rcmax=max(rcurve) 
        
                     rmax=rcmax       # TC 
        
                     #rmax=bkg        # TC 
        
                 if opts.conduct == 0: 
        
                     #rcurve_tot = rcurve_tot + [rcurve-bkg]        # RC                                 
        
                     bkgs = bkg  #RC 
        
                     rcmax=max(rcurve) 
        
                     rmax=rcmax-bkg # RC      
        
                 ###debug:   
        
                 #print 'rmax= ', rmax 
        
                 #--------------------------------------------------------------------------------------- 
        
           #     tot=sum(rcurve)/len(rcurve)      #; print 'tot = ', tot 
        
                 rcurve_l=list(rcurve) 
        
                 max_i=rcurve_l.index(rcmax) 
        
                 thi0=th[max_i]-th0 
        
                 #---------------------------------------------------------       
        
                 # define central pixel                                      
        
                 indx0=int(0.5*(xrange[0]+xrange[len(xrange)-1])) 
        
                 indy0=int(0.5*(yrange[0]+yrange[len(yrange)-1])) 
        
                 #-------------------------------------------------------- 
        
                 if (i==indx0) and (j==indy0): 
        
                       print("thi0 = ", thi0) 
        
                       #ugol=th-th0 #th_deg     
        
                       p0=[bkgs,rcmax,thi0,fwhm_0] #bkg max center 
        
                       p, unc, dte, stat, cov = fit1d(residuals,p0,ugol,rcurve)  
        
                       com0=sum((rcurve-p[0])*(ugol))/sum(rcurve-p[0]); print("com = ", com0) 
        
                       print('--------------------------------------------------------------------------') 
        
                       print('Gaussian fit at the central pixel') 
        
                       print('X0 [mm] = ', float(xrange[0] + indx0)*dx) 
        
                       print('Y0 [mm] = ', float(yrange[0] + indy0)*dy) 
        
                       print('bkg  = ', bkg) 
        
                       print('bkgs = ', bkgs) 
        
                       print('p = ', p) 
        
                       print('FWHM =', 2.0*abs(p[3])*sqrt(log(2)), ' urad') 
        
                       #-------------------------------------------------- 
        
                       # assign rcurve from pixel for output 
        
                       # subtract background  
        
                       prcurve = rcurve - bkgs                
        
                       ##################################################################################               
        
                       # start plotting things ---------------------------------- 
        
                       mpl.rcParams.update({'font.size': 12}) 
        
                       f1=plt.figure(figsize=(8.0,7.0)) 
        
                       plt.plot(ugol,prcurve, 'bo',label="local RC")  
        
                       plt.plot(ugol,prcurve, 'b-') 
        
                       ugol_g=arange(ugol[0],ugol[len(ugol)-1]+0.1,0.1) 
        
                       plt.plot(ugol_g,gauss(p,ugol_g)-p[0],'b--',label="Gauss. fit for local RC")             
        
                       plt.title ('pixel location [mm]: X ='+str(round((nx-i)*dx,3))+', Y = '+str(round((ny-j)*dy,3)),size=16)  
        
                       plt.xlabel('$\\theta - \\theta_0$, [$\mu$rad]', size=16) 
        
                       plt.ylabel('Counts', size=16)                                         
        
           ##################################################################################################################  
        
           ####### PIXEL REJECTION CRITERIA #################################################################################  
        
           ##################################################################################################################   
        
                 if opts.conduct == 0: #RC 
        
                     crit1 = bkg<dyn_range*bkg0 and rcmax>tr*bkg and rcmax<tot_range*bkg0             
        
                 elif opts.conduct == 1: #TC 
        
                     crit1 = bkg<dyn_range*bkg0 and rcmax>tr*bkg and rcmax<tot_range*bkg0 and bkg>bkg0 
        
                 # explicitly:           
        
                 #if bkg<dyn_range*bkg0 and rcmax>tr*bkg and rcmax<tot_range*bkg0:  #RC 
        
                 #if bkg<dyn_range*bkg0 and rcmax>tr*bkg and rcmax<tot_range*bkg0 and bkg>bkg0:            #TC  
        
                 if crit1 == True: 
        
                     # add rcurve to rcurve_tot - afer deglitching if it is activated 
        
                     # AND after pixel rejection criteria 
        
                     # subtracting bkg since it can vary from pixel to pixel       
        
                     if opts.conduct == 1: 
        
                         rcurve_tot = rcurve_tot + [-1.0*(rcurve-bkg_lin)]  # TC         
        
                     #rmax=bkg        # TC 
        
                     if opts.conduct == 0: 
        
                         rcurve_tot = rcurve_tot + [rcurve-bkg]        # RC                                 
        
                     ################################################################################## 
        
                     # extracting curve parameters if the pixel complies with the criterion 
        
                     ##################################################################################                         
        
                     if opts.stat==1: 
        
                           p0=[bkgs,rcmax,thi0,fwhm_0] 
        
                           p, unc, dte, stat, cov = fit1d(residuals,p0,th-th0,rcurve) 
        
                           if opts.conduct == 0: #RC 
        
                             crit2 = abs(p[1]) < dyn_range*bkg0 and abs(p[2]) < 0.5*abs(th_end-th_beg) and abs(p[3]) < abs(th_end-th_beg)  # and cov !='NA' # as an option 
        
                           elif opts.conduct == 1: #TC   
        
                             crit2 = abs(p[1]) < dyn_range*bkg0 and abs(p[1])>tr*bkg and abs(p[2]) < 0.5*abs(th_end-th_beg) and abs(p[3]) < abs(th_end-th_beg) 
        
                           #if abs(p[0]) < dyn_range*bkg0 and abs(p[2]) < 0.5*abs(th_end-th_beg) and abs(p[3]) < abs(th_end-th_beg):  # and cov !='NA' # as an option 
        
                           #if abs(p[0]) < dyn_range*bkg0 and abs(p[1])>tr*bkg and abs(p[2]) < 0.5*abs(th_end-th_beg) and abs(p[3]) < abs(th_end-th_beg):  # TC  
        
                           if crit2 == True: 
        
                             # 
        
                             fwhm_g=2.0*abs(p[3])*sqrt(log(2.0)) 
        
                             fwhm[i,j]=fwhm_g;   fwhm1=fwhm1+[fwhm[i,j]]                   
        
                             if zf == 0:  
        
                                 if opts.conduct == 0:    peak[i,j] = p[1]      #RC 
        
                                 elif opts.conduct == 1:  peak[i,j] = p[1]/bkg  #TC                  
        
                             elif zf == -1: 
        
                                 if opts.conduct == 0: peak[i,j] = abs(p[1])*abs(p[3])*sqrt(pi)/abs(zf) 
        
                                 if opts.conduct == 1: peak[i,j] = abs(p[1])*abs(p[3])*sqrt(pi)/abs(zf)/bkg 
        
                             else: 
        
                                 peak[i,j] = p[1]/abs(zf) 
        
                                 #if opts.conduct == 0:    peak[i,j] = p[1]      #RC 
        
                                 #elif opts.conduct == 1:  peak[i,j] = p[1]/bkg  #TC   
        
                             ##---------------------------------------------------------     
        
                             peak1=peak1+[peak[i,j]]   # p[1] - in gauss function background already subtracted                
        
                             thmid[i,j]=p[2];       thmid1=thmid1+[thmid[i,j]] 
        
                             #com[i,j]=sum((rcurve-p[0])*(th-th0))/sum(rcurve-p[0]); com1=com1+[com[i,j]]                                     
        
                             com[i,j]=p[2];         com1=com1+[com[i,j]]                
        
                             stdev[i,j]=abs(p[3])/sqrt(2.0); stdev1=stdev1+[stdev[i,j]] 
        
                             thneg[i,j]=p[2]-fwhm_g/2.0; thneg1=thneg1+[thneg[i,j]] 
        
                             thpos[i,j]=p[2]+fwhm_g/2.0; thpos1=thpos1+[thpos[i,j]] 
        
                           else: 
        
                             fwhm[i,j]=1e9 
        
                             peak[i,j]=1e9 
        
                             thmid[i,j]=1e9 
        
                             com[i,j]=1e9 
        
                             thneg[i,j]=1e9 
        
                             thpos[i,j]=1e9   
        
                             stdev[i,j]=1e9 
        
                     else: 
        
                         try:  
        
                           stat=curvestat(th-th0,rcurve,bkgs) 
        
                           #print 'stat4 =' , stat[4] 
        
                           #print 'stat5 = ', stat[5] 
        
                           #print 'stat6 = ', stat[6] 
        
                           #print 'stat7 = ', stat[7] 
        
                           if opts.conduct == 0: 
        
                             crit2 = abs(stat[1]-bkgs) < dyn_range*bkg0 and abs(stat[6]) < 0.5*abs(th_end-th_beg) and abs(stat[5]) < abs(th_end-th_beg) 
        
                           if opts.conduct == 1: 
        
                             crit2 = abs(stat[1]-bkgs) < dyn_range*bkg0 and abs(stat[1])>tr*bkg and abs(stat[6]) < 0.5*abs(th_end-th_beg) and abs(stat[5]) < abs(th_end-th_beg) 
        
                           # 
        
                           if crit2 == True: 
        
                             fwhm[i,j]=stat[5];         fwhm1=fwhm1+[fwhm[i,j]] 
        
                             if zf == 0: 
        
                                 if opts.conduct == 0:   peak[i,j]=stat[1]-bkgs   #RC 
        
                                 elif opts.conduct == 1: peak[i,j]=stat[1]/bkg    #TC 
        
                             elif zf == -1: 
        
                                 if opts.conduct == 0:   peak[i,j]=stat[8]/abs(zf) 
        
                                 elif opts.conduct == 1: peak[i,j]=stat[8]/abs(zf)/bkg    
        
                             else:  
        
                                 peak[i,j]=stat[1]/abs(zf) 
        
                                 #if opts.conduct == 0:   peak[i,j]=stat[1]/abs(zf)   #RC 
        
                                 #elif opts.conduct == 1: peak[i,j]=stat[1]/abs(zf)   #TC 
        
                                 #                                             
        
                             peak1=peak1+[peak[i,j]]    # peak intensity                  
        
                             thmid[i,j]=stat[4];        thmid1=thmid1+[thmid[i,j]] 
        
                             com[i,j]=stat[6];          com1=com1+[com[i,j]] 
        
                             stdev[i,j]=sqrt(abs(stat[7]));  stdev1=stdev1+[stdev[i,j]]  # variance to standard deviation 
        
                             thneg[i,j]=stat[2];        thneg1=thneg1+[thneg[i,j]] 
        
                             thpos[i,j]=stat[3];        thpos1=thpos1+[thpos[i,j]] 
        
                           else: 
        
                             fwhm[i,j]=1e9 
        
                             peak[i,j]=1e9 
        
                             thmid[i,j]=1e9 
        
                             com[i,j]=1e9 
        
                             thneg[i,j]=1e9 
        
                             thpos[i,j]=1e9 
        
                             stdev[i,j]=1e9                                                             
        
           #################################################################### 
        
           ## Testing 
        
           #################################################################### 
        
           #                if com[i,j] > 0: 
        
           #                  print "flag!", i, j 
        
                             #com[i,j]=-1e3                   
        
           #                if (i==xi) and (j==yi): 
        
           #                    print "fwhm check", fwhm[i,j] 
        
           #                    print "com check", com[i,j] 
        
           ####################################################################             
        
                         except UnboundLocalError: 
        
                           fwhm[i,j]=1e9 
        
                           peak[i,j]=1e9 
        
                           thmid[i,j]=1e9 
        
                           com[i,j]=1e9 
        
                           thneg[i,j]=1e9 
        
                           thpos[i,j]=1e9 
        
                           stdev[i,j]=1e9 
        
                 else:   
        
                   fwhm[i,j]=1e9 
        
                   peak[i,j]=1e9 
        
                   thmid[i,j]=1e9 
        
                   com[i,j]=1e9 
        
                   thneg[i,j]=1e9 
        
                   thpos[i,j]=1e9 
        
                   stdev[i,j]=1e9

The comments in the code suggest a set of subdivisions:

Background determination
Deglitching procedure
subtracting bkg since it can vary from pixel to pixel
define central pixel
printing and visualization
PIXEL REJECTION CRITERIA
extracting curve parameters if the pixel complies with the criterion
Testing (commented out)

As a first optimization, the printing and visualization should be separated from the computational code.

Background determination and deglitching probably both have library routines that could be used, even using the topo-specific math.

This code, for example, looks as if it may be a weighted average:

dtxrd/rctopo

Lines 455 to 462 in 1b2708c

    
           bkg=0.25*(rcurve[0]+rcurve[1]+rcurve[len(rcurve)-2]+rcurve[len(rcurve)-1]) #; print 'bkg = ', bkg   
        
           b1=0.5*(rcurve[0]+rcurve[1]) 
        
           b2=0.5*(rcurve[len(rcurve)-2]+rcurve[len(rcurve)-1]) 
        
           t1=th_beg 
        
           t2=th_end 
        
           A=(b2-b1)/(t2-t1) 
        
           B=b2-A*t2 
        
           bkg_lin=A*th+B

@sstoupin emailed:

Hi Pete,

Thanks for the quick response and suggestions.

I saw your edits [here] - makes sense. I will go though LGTM at some point - there are many redundancies for sure.

The deglitching procedure is something rarely used (its for noisy data) it gets enabled with an option, otherwise not functioning. Needs more testing anyway.

The bottleneck I think originates from the the double "for" loop where either curvestat or fit1d functions are used (both are my functions, perhaps there are better developments). I think I found now one that could be better. https://lmfit.github.io/lmfit-py/model.html

this can perhaps replace my fit1d (for gaussian profile fitting)

Not sure about curvestat - this one finds left/right slope positions (which correspond to FWHM).

I think some sort of vectorized form for these functions would be ideal, don't know yet if these exist.

I agree that vector functions will replace the need for the double loop. So far, my code changes are not for speed. One easy thing to take out of the loop is this part:

      if opts.bkg==-1:    
        bkg=0.25*(rcurve[0]+rcurve[1]+rcurve[len(rcurve)-2]+rcurve[len(rcurve)-1]) #; print 'bkg = ', bkg  
        b1=0.5*(rcurve[0]+rcurve[1])
        b2=0.5*(rcurve[len(rcurve)-2]+rcurve[len(rcurve)-1])
        t1=th_beg
        t2=th_end
        A=(b2-b1)/(t2-t1)
        B=b2-A*t2
        bkg_lin=A*th+B

At minimum, the if decision is executed 10^6 times. A vector replacement before the double for loop would execute zero or once depending on opts.bkg value. And, easy to demonstrate.

We should add some timing collection first, so we can identify if any change is an actual improvement of execution time.

This section is a crude evaluation of the rocking curve background:
bkg - constant, average of 4 points on tails
bkg_lin is a line
the background can be unique for each pixel but generally very similar (unless it is a bad/hot pixel of area detector).
What lines in the code would you suggest for time tracking?

For starters, just the overall time, as shown in changeset 2393444

	for j in yrange:
	for i in xrange:
	rcurve=data[:,j,i]
	#-------------------------------------------------
	# Background determination
	#-------------------------------------------------
	if opts.bkg==-1:
	bkg=0.25*(rcurve[0]+rcurve[1]+rcurve[len(rcurve)-2]+rcurve[len(rcurve)-1]) #; print 'bkg = ', bkg
	b1=0.5*(rcurve[0]+rcurve[1])
	b2=0.5*(rcurve[len(rcurve)-2]+rcurve[len(rcurve)-1])
	t1=th_beg
	t2=th_end
	A=(b2-b1)/(t2-t1)
	B=b2-A*t2
	bkg_lin=A*th+B
	#------------------------------------------------
	# Deglitching procedure
	#------------------------------------------------
	if df>1.0:
	index=range(len(rcurve))
	index=index[2:len(rcurve)-2]
	for k in index:
	y=rcurve[k]
	mean=0.25*(rcurve[k-2]+rcurve[k-1]+rcurve[k+1]+rcurve[k+2])
	if y > dfmean or y < (2.0-df)mean:
	rcurve[k]=mean
	print("X= ", dxi, "Y= ",dyj, " glitch found")
	#---------------------------------------------------------------------------------------
	# subtracting bkg since it can vary from pixel to pixel
	if opts.conduct == 1:
	#rcurve_tot = rcurve_tot + [-1.0*(rcurve-bkg_lin)] # TC
	rcurve = -1.0*(rcurve-bkg_lin) # TC
	bkgs = 0.0 # TC
	rcmax=max(rcurve)
	rmax=rcmax # TC
	#rmax=bkg # TC
	if opts.conduct == 0:
	#rcurve_tot = rcurve_tot + [rcurve-bkg] # RC
	bkgs = bkg #RC
	rcmax=max(rcurve)
	rmax=rcmax-bkg # RC
	###debug:
	#print 'rmax= ', rmax
	#---------------------------------------------------------------------------------------
	# tot=sum(rcurve)/len(rcurve) #; print 'tot = ', tot
	rcurve_l=list(rcurve)
	max_i=rcurve_l.index(rcmax)
	thi0=th[max_i]-th0
	#---------------------------------------------------------
	# define central pixel
	indx0=int(0.5*(xrange[0]+xrange[len(xrange)-1]))
	indy0=int(0.5*(yrange[0]+yrange[len(yrange)-1]))
	#--------------------------------------------------------
	if (i==indx0) and (j==indy0):
	print("thi0 = ", thi0)
	#ugol=th-th0 #th_deg
	p0=[bkgs,rcmax,thi0,fwhm_0] #bkg max center
	p, unc, dte, stat, cov = fit1d(residuals,p0,ugol,rcurve)
	com0=sum((rcurve-p[0])*(ugol))/sum(rcurve-p[0]); print("com = ", com0)
	print('--------------------------------------------------------------------------')
	print('Gaussian fit at the central pixel')
	print('X0 [mm] = ', float(xrange[0] + indx0)*dx)
	print('Y0 [mm] = ', float(yrange[0] + indy0)*dy)
	print('bkg = ', bkg)
	print('bkgs = ', bkgs)
	print('p = ', p)
	print('FWHM =', 2.0abs(p[3])sqrt(log(2)), ' urad')
	#--------------------------------------------------
	# assign rcurve from pixel for output
	# subtract background
	prcurve = rcurve - bkgs
	##################################################################################
	# start plotting things ----------------------------------
	mpl.rcParams.update({'font.size': 12})
	f1=plt.figure(figsize=(8.0,7.0))
	plt.plot(ugol,prcurve, 'bo',label="local RC")
	plt.plot(ugol,prcurve, 'b-')
	ugol_g=arange(ugol[0],ugol[len(ugol)-1]+0.1,0.1)
	plt.plot(ugol_g,gauss(p,ugol_g)-p[0],'b--',label="Gauss. fit for local RC")
	plt.title ('pixel location [mm]: X ='+str(round((nx-i)dx,3))+', Y = '+str(round((ny-j)dy,3)),size=16)
	plt.xlabel('$\\theta - \\theta_0$, [$\mu$rad]', size=16)
	plt.ylabel('Counts', size=16)
	##################################################################################################################
	####### PIXEL REJECTION CRITERIA #################################################################################
	##################################################################################################################
	if opts.conduct == 0: #RC
	crit1 = bkg<dyn_rangebkg0 and rcmax>trbkg and rcmax<tot_range*bkg0
	elif opts.conduct == 1: #TC
	crit1 = bkg<dyn_rangebkg0 and rcmax>trbkg and rcmax<tot_range*bkg0 and bkg>bkg0
	# explicitly:
	#if bkg<dyn_rangebkg0 and rcmax>trbkg and rcmax<tot_range*bkg0: #RC
	#if bkg<dyn_rangebkg0 and rcmax>trbkg and rcmax<tot_range*bkg0 and bkg>bkg0: #TC
	if crit1 == True:
	# add rcurve to rcurve_tot - afer deglitching if it is activated
	# AND after pixel rejection criteria
	# subtracting bkg since it can vary from pixel to pixel
	if opts.conduct == 1:
	rcurve_tot = rcurve_tot + [-1.0*(rcurve-bkg_lin)] # TC
	#rmax=bkg # TC
	if opts.conduct == 0:
	rcurve_tot = rcurve_tot + [rcurve-bkg] # RC
	##################################################################################
	# extracting curve parameters if the pixel complies with the criterion
	##################################################################################
	if opts.stat==1:
	p0=[bkgs,rcmax,thi0,fwhm_0]
	p, unc, dte, stat, cov = fit1d(residuals,p0,th-th0,rcurve)
	if opts.conduct == 0: #RC
	crit2 = abs(p[1]) < dyn_rangebkg0 and abs(p[2]) < 0.5abs(th_end-th_beg) and abs(p[3]) < abs(th_end-th_beg) # and cov !='NA' # as an option
	elif opts.conduct == 1: #TC
	crit2 = abs(p[1]) < dyn_rangebkg0 and abs(p[1])>trbkg and abs(p[2]) < 0.5*abs(th_end-th_beg) and abs(p[3]) < abs(th_end-th_beg)
	#if abs(p[0]) < dyn_rangebkg0 and abs(p[2]) < 0.5abs(th_end-th_beg) and abs(p[3]) < abs(th_end-th_beg): # and cov !='NA' # as an option
	#if abs(p[0]) < dyn_rangebkg0 and abs(p[1])>trbkg and abs(p[2]) < 0.5*abs(th_end-th_beg) and abs(p[3]) < abs(th_end-th_beg): # TC
	if crit2 == True:
	#
	fwhm_g=2.0abs(p[3])sqrt(log(2.0))
	fwhm[i,j]=fwhm_g; fwhm1=fwhm1+[fwhm[i,j]]
	if zf == 0:
	if opts.conduct == 0: peak[i,j] = p[1] #RC
	elif opts.conduct == 1: peak[i,j] = p[1]/bkg #TC
	elif zf == -1:
	if opts.conduct == 0: peak[i,j] = abs(p[1])abs(p[3])sqrt(pi)/abs(zf)
	if opts.conduct == 1: peak[i,j] = abs(p[1])abs(p[3])sqrt(pi)/abs(zf)/bkg
	else:
	peak[i,j] = p[1]/abs(zf)
	#if opts.conduct == 0: peak[i,j] = p[1] #RC
	#elif opts.conduct == 1: peak[i,j] = p[1]/bkg #TC
	##---------------------------------------------------------
	peak1=peak1+[peak[i,j]] # p[1] - in gauss function background already subtracted
	thmid[i,j]=p[2]; thmid1=thmid1+[thmid[i,j]]
	#com[i,j]=sum((rcurve-p[0])*(th-th0))/sum(rcurve-p[0]); com1=com1+[com[i,j]]
	com[i,j]=p[2]; com1=com1+[com[i,j]]
	stdev[i,j]=abs(p[3])/sqrt(2.0); stdev1=stdev1+[stdev[i,j]]
	thneg[i,j]=p[2]-fwhm_g/2.0; thneg1=thneg1+[thneg[i,j]]
	thpos[i,j]=p[2]+fwhm_g/2.0; thpos1=thpos1+[thpos[i,j]]
	else:
	fwhm[i,j]=1e9
	peak[i,j]=1e9
	thmid[i,j]=1e9
	com[i,j]=1e9
	thneg[i,j]=1e9
	thpos[i,j]=1e9
	stdev[i,j]=1e9
	else:
	try:
	stat=curvestat(th-th0,rcurve,bkgs)
	#print 'stat4 =' , stat[4]
	#print 'stat5 = ', stat[5]
	#print 'stat6 = ', stat[6]
	#print 'stat7 = ', stat[7]
	if opts.conduct == 0:
	crit2 = abs(stat[1]-bkgs) < dyn_rangebkg0 and abs(stat[6]) < 0.5abs(th_end-th_beg) and abs(stat[5]) < abs(th_end-th_beg)
	if opts.conduct == 1:
	crit2 = abs(stat[1]-bkgs) < dyn_rangebkg0 and abs(stat[1])>trbkg and abs(stat[6]) < 0.5*abs(th_end-th_beg) and abs(stat[5]) < abs(th_end-th_beg)
	#
	if crit2 == True:
	fwhm[i,j]=stat[5]; fwhm1=fwhm1+[fwhm[i,j]]

	if zf == 0:
	if opts.conduct == 0: peak[i,j]=stat[1]-bkgs #RC
	elif opts.conduct == 1: peak[i,j]=stat[1]/bkg #TC
	elif zf == -1:
	if opts.conduct == 0: peak[i,j]=stat[8]/abs(zf)
	elif opts.conduct == 1: peak[i,j]=stat[8]/abs(zf)/bkg
	else:
	peak[i,j]=stat[1]/abs(zf)
	#if opts.conduct == 0: peak[i,j]=stat[1]/abs(zf) #RC
	#elif opts.conduct == 1: peak[i,j]=stat[1]/abs(zf) #TC
	#
	peak1=peak1+[peak[i,j]] # peak intensity
	thmid[i,j]=stat[4]; thmid1=thmid1+[thmid[i,j]]
	com[i,j]=stat[6]; com1=com1+[com[i,j]]
	stdev[i,j]=sqrt(abs(stat[7])); stdev1=stdev1+[stdev[i,j]] # variance to standard deviation
	thneg[i,j]=stat[2]; thneg1=thneg1+[thneg[i,j]]
	thpos[i,j]=stat[3]; thpos1=thpos1+[thpos[i,j]]
	else:
	fwhm[i,j]=1e9
	peak[i,j]=1e9
	thmid[i,j]=1e9
	com[i,j]=1e9
	thneg[i,j]=1e9
	thpos[i,j]=1e9
	stdev[i,j]=1e9
	####################################################################
	## Testing
	####################################################################
	# if com[i,j] > 0:
	# print "flag!", i, j
	#com[i,j]=-1e3
	# if (i==xi) and (j==yi):
	# print "fwhm check", fwhm[i,j]
	# print "com check", com[i,j]
	####################################################################

	except UnboundLocalError:
	fwhm[i,j]=1e9
	peak[i,j]=1e9
	thmid[i,j]=1e9
	com[i,j]=1e9
	thneg[i,j]=1e9
	thpos[i,j]=1e9
	stdev[i,j]=1e9

	else:
	fwhm[i,j]=1e9
	peak[i,j]=1e9
	thmid[i,j]=1e9
	com[i,j]=1e9
	thneg[i,j]=1e9
	thpos[i,j]=1e9
	stdev[i,j]=1e9