add: astropy;change: bind fit

.gitignore

@@ -17,5 +17,8 @@ __pycache__
 .vscode
 .vs
 
+# dev
+*.ipynb
+
 *.code-workspace
 *.exe

main.py

@@ -10,21 +10,24 @@ BF = BindFilter()
 path = 'result/bind'
 file_list = os.listdir(path)
 
-reg = re.compile(r'(([0-9]{4})-([0-9])-([0.9])).txt')
+reg = re.compile(r'(([0-9]{4})-([0-9])-([0-9])).txt')
 
 pbar = tqdm(desc="Gaussian Mixture Bind Filter", total=len(file_list))
-
 for file in file_list:
     name, E, n, m = reg.match(file).groups()
     file = os.path.join(path, file)
 
-    BF(file)
-    BF.draw('result/GMM/' + name + '.png')
-    np.savetxt('result/bind-GMM/' + name + '.txt', BF.fit_data, fmt='%d')
+    BF.filter(file)
+    # BF.draw('result/GMM/' + name + '.png')
+    # np.savetxt('result/bind-GMM/' + name + '.txt', BF.fit_data, fmt='%d')
     binds.append(Bind(int(E), int(n), int(m), BF.fit_data))
 
     pbar.update(1)
-
-    break
-
+pbar.close()
+
+pbar = tqdm(desc="Bind Linear Fit", total=len(binds))
+for bind in binds:
+    bind.fit()
+    bind.draw('result/L-FIT/' + bind.name + '.png')
+    pbar.update(1)
 pbar.close()

qdx/Bind.py

@@ -1,28 +1,107 @@
-from cProfile import label
 import numpy as np
 import matplotlib.pyplot as plt
+from astropy.modeling import models, fitting
+from sklearn.mixture import GaussianMixture
 from sklearn.linear_model import LinearRegression
 
+def get_hist(data, maxN=50):
+    step = 1
+    edge = np.arange(data.min(), data.max() + 1, step)
+    count, _ = np.histogram(data, bins=edge)
+    while count.max() <= maxN:
+        step += 1
+        edge = np.arange(data.min(), data.max() + 1, step)
+        count, _ = np.histogram(data, bins=edge)
+    return count, (edge[1:] + edge[:-1]) / 2
 
 class Bind(object):
-    k, b = 0, 0
+    flag = 0
+    L, C = 40, 0
+    k1, k2, b = 0, 0, 0
+    K1, C1, K2, C2 = 0, 0, 0, 0
 
-    def __init__(self, E, n, m, data):
-        self.E = E
+    def __init__(self, n, m):
         self.n, self.m = n, m
-        self.x, self.y = data[:, 0], data[:, 1]
+
+    def add_data(self, E, data):
+        if self.flag == 0:
+            self.flag = 1
+            self.E1 = E
+            self.x1, self.y1 = data[:, 0].reshape(-1, 1), data[:, 1].reshape(-1, 1)
+        else:
+            self.E2 = E
+            self.x2, self.y2 = data[:, 0].reshape(-1, 1), data[:, 1].reshape(-1, 1)
 
     def fit(self):
-        self.reg = LinearRegression()
-        self.reg.fit(self.x, self.y)
-        self.k = self.reg.coef_[0][0]
-        self.b = self.reg.intercept_[0]
+        self.reg1 = LinearRegression()
+        self.reg1.fit(self.x1, self.y1)
+        self.K1 = self.reg1.coef_[0][0]
+        self.C1 = self.reg1.intercept_[0]
+
+        self.reg2 = LinearRegression()
+        self.reg2.fit(self.x2, self.y2)
+        self.K2 = self.reg2.coef_[0][0]
+        self.C2 = self.reg2.intercept_[0]
+
+    def solve(self):
+        self.k1 = (self.E1 - self.E2) / (self.C1 - self.C2)
+        self.k2 = -(self.K1 + self.K2) / (self.k1 * 2)
+        self.b = (self.E1 - self.C1 * self.k1 + self.E2 - self.C2 * self.k1) / 2
+
+    def split(self, n=7):
+        self.cluster = []
+        data = np.concatenate((self.x1, self.y1), axis=1)
+
+        model = GaussianMixture(n_components=n)
+        model.fit(data)
+
+        ny = model.predict(data)
+        for i in np.unique(ny):
+            idx = np.where(ny == i)[0]
+            self.cluster.append(data[idx])
+
+        self.cluster = np.array(self.cluster)
+
+    def fit2(self):
+        self.fx = []
+        self.fy = []
+        fitter = fitting.LevMarLSQFitter()
+
+        for data in self.cluster:
+            x1, x2 = data[:, 0], data[:, 1]
+            c1, e1 = get_hist(x1)
+            c2, e2 = get_hist(x2)
+            modelx1 = models.Gaussian1D(amplitude=c1.max(), mean=x1.mean(), stddev=x1.std())
+            modelx2 = models.Gaussian1D(amplitude=c2.max(), mean=x2.mean(), stddev=x2.std())
+            fitted_model1 = fitter(modelx1, e1, c1)
+            fitted_model2 = fitter(modelx2, e2, c2)
+            self.fx.append(fitted_model1.mean.value)
+            self.fy.append(fitted_model2.mean.value)
+
+    def pX(self, x1=None, x2=None):
+        x1 = x1 if x1 else self.x1
+        x2 = x2 if x2 else self.y1
+        self.px = (self.k2 * x2 - self.k1 * x1) * self.L / self.E1 + self.C
 
     def draw(self, title):
         fig = plt.figure(figsize=(8, 8))
         ax = fig.add_subplot(1, 1, 1)
-        ax.scatter(self.x, self.y, s=0.1, c='k', label='raw')
-        ax.plot(self.x, self.reg.predict(self.x), c='r', label='fit')
+        ax.scatter(self.x1, self.y1, s=0.1, c='black', label=r'$E_1$')
+        ax.scatter(self.x2, self.y2, s=0.1, c='dimgray', label=r'$E_2$')
+        ax.plot(self.x1, self.reg1.predict(self.x1), c='red', label=r'$x_2={:.4f}x_1+{:.4f},\ R^2={:.5f}$'.format(self.K1, self.C1, self.RSquare1))
+        ax.plot(self.x2, self.reg2.predict(self.x2), c='orangered', label=r'$x_2={:.4f}x_1+{:.4f},\ R^2={:.5f}$'.format(self.K2, self.C2, self.RSquare2))
         ax.legend()
-        fig.savefig(title)
+        fig.savefig(title, facecolor='w', transparent=False)
         plt.close()
+
+    @property
+    def name(self):
+        return '{:d}-{:d}-{:d}-{:d}'.format(self.E1, self.E2, self.n, self.m)
+
+    @property
+    def RSquare1(self):
+        return self.reg1.score(self.x1, self.y1)
+
+    @property
+    def RSquare2(self):
+        return self.reg2.score(self.x2, self.y2)

qdx/BindFilter.py

@@ -2,12 +2,11 @@ import numpy as np
 from matplotlib import pyplot as plt
 from sklearn.mixture import GaussianMixture
 
-
 class BindFilter(object):
     def __init__(self):
         pass
 
-    def fit(self, file):
+    def filter(self, file):
         self.clusters = []
         self.fit_data = np.array([])
 
@@ -22,6 +21,7 @@ class BindFilter(object):
             self.fit_data = idx if len(idx) > len(self.fit_data) else self.fit_data
             self.clusters.append(data[idx])
 
+        self.data = data
         self.fit_data = data[self.fit_data]
 
     def draw(self, title):

qdx/Block.py

@@ -1,6 +1,5 @@
 from .Bind import Bind
 
-
 class Block(object):
     binds = []