Snippets

Rene Cejas Bolecek Calorimetry Filtering tools: Average - Outlier

Created by Rene Cejas Bolecek last modified
filterAvOutlierV1.py
Raw 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
#!/usr/bin/python

##################################################################################
# Average C1 values - Outlier v1
# Rel. 10/2015
# Dr. Rene Cejas Bolecek
# Low Temperatures Laboratory, CAB, Argentine.
# email: reneczechdev@gmail.com
##################################################################################

#Libraries

from numpy import zeros, sqrt, where, pi, average, arange, histogram,log,log10 # log is ln
import numpy as np
from PyAstronomy import pyasl 
#http://www.hs.uni-hamburg.de/DE/Ins/Per/Czesla/PyA/PyA/pyaCDoc/installingPyA.html


#User settings
filename = '14_10T_norm_f1'
ColT = 0 #En los archivos de J.G. es 0
ColC1 = 6 #En los archivos de J.G. es 6
#Av. settings
Tcutoff = 0.01

#Outlier settings
polyFit = 1
stdlim = 0.05 
flagDEB = False #True or False

###################################################################
#Output Files
writeIt1 = open(filename+"_Av_V1.dat", "w")
writeIt2 = open(filename+"_Outlier_V1.dat", "w")
###################################################################
#Gnuplot check
#P3(x)=a0+a1*x+a2*x**2+a3*x**3
#a0 = 1.99699e-06-8E-6
#a1 = 5.93816e-07
#a2 = 3.49326e-09
#a3= -3.94088e-11
#plot "14_10T_norm.dat" u 1:7, P3(x), "14_10T_norm_f1.dat" u 1:7, "14_10T_norm_f1_AvV1.dat" u 1:3, "14_10T_norm_f1_OutlierV1.dat" u 1:2
#Var init
Tav = 0
C1av = 0
i0 = 0
flag = True
data   = np.loadtxt(filename+'.dat')
T = data[:,ColT]
C1 = data[:,ColC1]

#Calculation init
for i in range(len(data)-1):
	TCalc = []
	
	TCalc.append(T[i])
	
	C1Calc = []
	
	C1Calc.append(C1[i])
	j = i+1
	
	if (i == i0 and j<=(len(data)-1)):
		flag = True
	if flag:
		counter = 1
		while (True and j <= (len(data)-1)):

			if (abs(T[i0]-T[j])/T[j])> Tcutoff:
				i0 = j
				#http://www.hs.uni-hamburg.de/DE/Ins/Per/Czesla/PyA/PyA/pyaslDoc/aslDoc/outlier.html
				#iin, iout,p , model = pyasl.polyResOutlier(TCalc, C1Calc, deg=3, stdlim=1, controlPlot=True,fullOutput = True)
				iin, iout= pyasl.polyResOutlier(TCalc, C1Calc, deg=polyFit, stdlim=stdlim, controlPlot=False)
				TOutliered = zeros(len(iin))
				C1Outliered = zeros(len(iin))
				for idx in range(len(iin)):
					print idx, iin
					TOutliered[idx] = TCalc[iin[idx]] 
					C1Outliered[idx] = C1Calc[iin[idx]]
					writeIt2.write("%.5f %.10f\n" % (TOutliered[idx],C1Outliered[idx]))
				writeIt1.write("%.5f %.5f %.10f %.10f %d\n" % (np.mean(TCalc),np.std(TCalc),np.mean(C1Calc),np.std(C1Calc),len(C1Calc)))
				break
			
			flag = False
			C1Calc.append(C1[j])
			TCalc.append(T[j])
			counter = counter +1
			if (j == (len(data)-1)):
				if(flagDEB and counter!=0): print "save Tav:", i," ",Tav/(counter*1.0)
				writeIt1.write("%.5f %.5f %.10f %.10f %d\n" % (np.mean(TCalc),np.std(TCalc),np.mean(C1Calc),np.std(C1Calc),len(C1Calc)))
			j = j+1
			if(flagDEB): print "Tav: ", Tav, "j: ", j ,"counter: ",  counter
			

writeIt1.close()	
writeIt2.close()	
print "Process finished"

Comments (0)

HTTPS SSH

You can clone a snippet to your computer for local editing. Learn more.