Package pysls :: Module lib
[frames] | no frames]

Source Code for Module pysls.lib

  1  """ 
  2  Various aux functions. 
  3  There should be more bla bla at places like this ... 
  4   
  5  """ 
  6   
  7  import scipy.signal as spsig 
  8  import numpy as np 
  9  import math 
 10  import sp, li, pl, tem, ca 
 11  import copy as pythoncopy 
 12  import matplotlib.pyplot as plt 
 13  import cPickle 
 14   
 15   
 16   
17 -def peak_detect(array,n):
18 """detects the n highest peaks in the array and returns their wavelength and their height in a list of two lists""" 19 20 kernel = 1 21 arr = spsig.medfilt(array, kernel)#permits to eventually smooth the signal, the second role of this part is to avoid the modification of the input array 22 size = len(arr) 23 min = np.min(arr) 24 25 m = [] #it will contain the indices of the peaks 26 peaksheight = [] 27 i = 1 #counts the number of peaks detected 28 while i <= n: 29 30 index = np.argmax(arr) 31 peakheight = max(arr) 32 if index + 1 < size and index - 1 >= 0 and arr[index+1] != min and arr[index-1] != min: #the two first conditions avoid errors with the two last ; 33 i += 1 34 m.append(index+1) 35 peaksheight.append(peakheight) 36 delta1 = - 1 37 delta2 = - 1 38 j = 0 39 while delta2 < 1 and delta1 < 1:#the choice of '1' is arbitrary, it avoids that a very small increase is detected as a peak, be conscious that this value makes sense to detect peaks in the cross correlation, for a spectra for instance it should be for example '0.05' 40 j += 1 # it represents somehow the width of the peak 41 42 if index - j > 0: 43 delta1 = arr[index-j-1]-arr[index-j] 44 else: 45 delta1 = 1000 #it is to go out of the 'while' structure 46 arr[index] = min 47 if index + j < size -1: 48 delta2 = arr[index+j+1]-arr[index+j] 49 else: 50 delta2 = 1000 51 arr[index] = min 52 53 arr[index-j:index+j+1] = min 54 if np.all(arr == min): #it means there is no more peaks to find 55 i = n + 1 56 57 m = np.array(m) 58 peaksheight = np.array(peaksheight) 59 return np.row_stack((m, peaksheight))#eventuellement donner la largeur, the peaks are automatically sorted from to highest to the smallest
60 61 62
63 -def mean_spectra(spectralist):
64 """calculates the average flux (with and without continum) of a spectra listand returns a spectra""" 65 66 m = [] 67 specnumber = len(spectralist) 68 for j in range(specnumber): 69 m.append(spectralist[j].ztopix(spectralist[j].zsdss)) 70 71 flux_cs, flux, noise = [], [], [] 72 for i in range(len(spectralist[0].fluxes_cs)): 73 a, b, c = [], [], [] 74 for j in range(len(spectralist)): 75 index = i+m[j] 76 if index < len(spectralist[j].fluxes_cs): 77 a.append(spectralist[j].fluxes_cs[index]) 78 b.append(spectralist[j].fluxes[index]) 79 c.append(spectralist[j].noise[i]) 80 flux_cs.append(np.mean(a)) 81 flux.append(np.mean(b)) 82 noise.append(np.mean(c)) 83 84 spec = sp.spectra(np.array(flux_cs),np.array( flux), np.array(noise)) 85 name = 'mean spectra from ' + str(specnumber) + ' files' 86 spec.name = name 87 spec.setangst() 88 spec.calc_csc() 89 return spec
90 91
92 -def candidat_search(directory_in, directory_out, file, tuplelist, cand_id, minpeaks, minqual, minborder, maxborder, noisetuples = [], masktuples = [], fluxchoice = 'fluxes_csc', n = 10, ratiolimit = 0.5):
93 """from a fit file, it writes candidats files (if they satisfy the selection criteria). It returns the number of candidats created. This function is quite complicated, it calls many other functions !""" 94 cand_number = 0 95 spec = sp.spectra() 96 spec = spec.read(directory_in, file) 97 98 spec.initmask() 99 spec.mask_borders(minborder, maxborder) 100 spec.mask_angst(noisetuples) 101 spec.mask_angst(masktuples, spec.zsdss) 102 103 template = tem.factory(tuplelist) 104 105 ccsp = sp.crosscorrspec(spec, template) 106 ccsp.crosscorr(fluxchoice) 107 108 peaks = peak_detect(ccsp.corr, n) 109 print len(peaks[0]), 110 if peaks != []: 111 mainpeaks = mainpeaksdet(peaks, spec) 112 print len(mainpeaks[0]), 113 z = mainpeaks[0] 114 115 for i in range(len(mainpeaks[0])): 116 temp = template.copy() 117 temp.setglobalz(z[i]) 118 temp.optimize(ccsp.spectra, fluxchoice) 119 120 cand_number += create_cand(spec, temp, z[i], directory_out, ratiolimit, cand_id, cand_number, minpeaks, minqual, minborder, maxborder) 121 return cand_number
122 123
124 -def create_cand(spec, temp, z, directory_out, ratiolimit, cand_id, cand_number, minpeaks, minqual, minborder, maxborder):
125 """This is more precisely the function that creats the candidats files, it applies the criteria of selection""" 126 ratio = sigtonoise(spec, temp) 127 number = 0 128 129 lines = [] 130 quality = 0 131 132 for j in range(len(ratio)): 133 if ratio[j] > ratiolimit: 134 temp.linelist[j].strength = ratio[j] 135 lines.append(temp.linelist[j]) 136 quality += ratio[j] 137 138 peaksnumber = len(lines) 139 if peaksnumber >= minpeaks and quality >= minqual: #ou 0 140 141 number += 1 142 143 name = str(quality)[0:5] + '_' + str(peaksnumber).zfill(2) + '_' + str(cand_id+cand_number+number).zfill(5) + '.ca' 144 cand = ca.candidat(name, spec, lines, z, quality) 145 cand.write(directory_out) 146 147 return number
148 149
150 -def sigtonoise(spectra, tem):
151 """it returns a list containing the ratio between the signal and the noise for the lines in the template""" 152 ratiolist = [] 153 154 for i in range(len(tem.linelist)): 155 156 width = tem.linelist[i].width 157 centerangst = tem.linelist[i].angst 158 159 n = spectra.angsttopix(centerangst) 160 i1 = spectra.angsttopix(centerangst)-spectra.angsttopix(centerangst-width) 161 i2 = spectra.angsttopix(centerangst+width) - spectra.angsttopix(centerangst) 162 163 av_noise = np.mean(spectra.noise[n - i1:n + i2 + 1]) 164 if av_noise != 0: 165 ratiolist.append(tem.linelist[i].height/av_noise) 166 else: 167 ratiolist.append(0) 168 169 return ratiolist
170
171 -def mainpeaksdet(peaks, sp):
172 """returns a list of the 'real' peaks in a signal, it checks that two peaks very close to each other are not considered as two different ones ; the minimum distance bewteen to peaks is arbitrary""" 173 z = sp.pixtoz(peaks[0]) 174 heights = peaks[1] 175 mainpeaks = [[z[0]],[heights[0]]] 176 177 for i in range(len(z)-2): 178 deltax1 = np.absolute(z[i+1] - z[i]) 179 if deltax1 > 0.03: 180 mainpeaks[0].append(z[i+1]) 181 mainpeaks[1].append(heights[i+1]) 182 return mainpeaks
183
184 -def sortlist(directory, filelist):
185 """returns a list of candidats sorted by quality""" 186 quality = [] 187 for file in filelist: 188 qual = file[len(directory):len(directory)+5] #this is not a very safe way to read the quality from the candidat 189 quality.append(float(qual)) 190 index = np.argsort(quality) 191 index = index[::-1] 192 filelist = np.array(filelist) 193 return filelist[index]
194 195 196
197 -def meanfilt(arr, N):
198 """filters a signal using a mean fliter ; N must be odd [not used]""" 199 n = (N-1)/2 200 size = len(arr) 201 a = [] 202 for i in range(size): 203 if i < n : 204 a.append(np.mean(arr[0:i+n+1])) 205 elif i > size - n - 1 : 206 a.append(np.mean(arr[i-n:size])) 207 else: 208 a.append(np.mean(arr[i-n:i+n+1])) 209 return np.array(a) 210 211
212 -def select_lines(reftuples, z, minangst, maxangst):
213 """it is to check that the wavelengths contained in the tuples are in the range of the spectrometer""" 214 useful_tuples = [] 215 216 for tuple in reftuples: 217 if tuple[1] * (1 + z) > minangst and tuple[1] * (1 + z) < maxangst: 218 useful_tuples.append(tuple) 219 220 return useful_tuples
221