В этой статье будет приведен очень короткий код поиска периодов тропического года, звездного года и звездных суток в данных радиоактивного распада плутония.
Вот здесь на GitHub.com можно скачать все данные и код
Всего данные по почти 19 дней одного года days.dat
И 2 дня в мае 2005 года 0505.txt
И 2 дня в мае 2006 года 0506.txt
Средняя продолжительность тропического года по весеннему равноденствию это на Января 1, 2000 365.2421897 или 365 дней, 5 часов, 48 минут, 45.19 секунд.
Звездный год это 365 дней и 369 минут
Звездные сутки 1436 минут движение Земли в течении которого планета совершает полный оборот относительно сферы неподвижных звезд.
Загружаем данные и расфасовываем по рядам по 60 значений
import mathimport randomfrom statistics import meankoeff=60;matrix = [line.strip() for line in open('/users/andrejeremcuk/downloads/0505.txt')];matrix1 = [line.strip() for line in open('/users/andrejeremcuk/downloads/0506.txt')];matrixd = [line.strip() for line in open('/users/andrejeremcuk/downloads/days.dat')];arra=[[ 0 for e in range(koeff)] for t in range(int(len(matrix)/koeff))];harra=[0 for t in range(int(len(matrix)/koeff))];arra1=[[ 0 for e in range(koeff)] for t in range(int(len(matrix)/koeff))];harra1=[0 for t in range(int(len(matrix)/koeff))];arrad=[[ 0 for e in range(koeff)] for t in range(int(len(matrixd)/koeff))];harrad=[0 for t in range(int(len(matrixd)/koeff))];for i in range(len(matrix)): matrix[i]=int(matrix[i]);matrix1[i]=int(matrix1[i]);for i in range(len(matrixd)): matrixd[i]=int(matrixd[i]);
Сортируем функцией sorted по возрастанию
z=0;for jk in range(int(len(matrix)/koeff)): #int(len(matrix)/koeff) for mk in range(koeff): arra[jk][mk]=float(matrix[z]);arra1[jk][mk]=float(matrix1[z]);z=z+1;for jk in range(int(len(matrix)/koeff)): harra[jk]=sorted(arra[jk]);harra1[jk]=sorted(arra1[jk])z=0;for jk in range(int(len(matrixd)/koeff)): for mk in range(koeff): arrad[jk][mk]=float(matrixd[z]);z=z+1;for jk in range(int(len(matrixd)/koeff)): harrad[jk]=sorted(arrad[jk]);
Делим на среднее значение
for jk in range(int(len(matrix)/koeff)): average=mean(harra[jk]) for mk in range(koeff): harra[jk][mk]=harra[jk][mk]/average average1=mean(harra1[jk]) for mk in range(koeff): harra1[jk][mk]=harra1[jk][mk]/average1for jk in range(int(len(matrixd)/koeff)): averaged=mean(harrad[jk]) for mk in range(koeff): harrad[jk][mk]=harrad[jk][mk]/averaged
Далее привожу функцию сравнения:
import numpyimport numpy as npfrom numpy.fft import fft, fftfreqimport matplotlib.pyplot as pltfrom numpy import linspace, loadtxt, ones, convolvefrom numpy import array, zeros, argmin, inffrom numpy.linalg import normdef _trackeback(D): i, j = array(D.shape) - 1 p, q = [i], [j] while (i > 0 and j > 0): tb = argmin((D[i-1, j-1], D[i-1, j], D[i, j-1])) if (tb == 0): i = i - 1 j = j - 1 elif (tb == 1): i = i - 1 elif (tb == 2): j = j - 1 p.insert(0, i) q.insert(0, j) p.insert(0, 0) q.insert(0, 0) return (array(p), array(q))def dtw(dr): r, c = len(dr), len(dr[0]); D = zeros((r + 1, c + 1)); D[0, 1:] = inf D[1:, 0] = inf for i in range(r): for j in range(c): D[i+1, j+1] = dr[i][j]; for i in range(r): for j in range(c): D[i+1, j+1] += min(D[i, j], D[i, j+1], D[i+1, j]) D = D[1:, 1:]; dist = D[-1, -1] / sum(D.shape); return dist, D, _trackeback(D)
def stat(x,y): pikk=[[float("inf") for m in range(len(y))] for t in range(len(x))]; for e in range(len(x)): for l in range(len(y)): pikk[e][l]=(((abs(x[e]-y[l])))*(abs(e-l)+1)); dt=dtw(pikk) p1=dt[2][0];p=dt[2][1]; res=0;uj=0;rmu=0; for e in range(len(p)): for l in range(len(p1)): pl=p[l];pe=p[e];ml=p1[l];me=p1[e]; #pl=l;pe=e;ml=l;me=e; ret=0;re=1;rea=0;red=0;ab=abs(y[pl]-y[pe]);ac=(pe-pl);rmt=abs(1-x[ml]/y[pl]);rmr=abs(1-x[me]/y[pe]); au=abs(x[ml]-x[me]);ao=(me-ml);hip=math.sqrt(ab**2+ac**2);hipe=math.sqrt(ao**2+au**2); if (1<ml<58)&(1<pl<58)&(1<me<58)&(1<pe<58): rea=abs(1-abs(abs(y[pl-1]*y[pl]*y[pl+1]*y[pl+2])/y[pl]**4))+abs(1-abs((y[pe-1]*y[pe]*y[pe+1]*y[pe+2])/y[pe]**4)) rea+=abs(1-abs(abs(x[ml-1]*x[ml]*x[ml+1]*x[ml+2])/x[ml]**4))+abs(1-abs((x[me-1]*x[me]*x[me+1]*x[me+2])/x[me]**4)) else: rea=0; if (hip!=0)and(hipe!=0): if (math.asin(abs(au)/hipe)!=0)&(math.asin((ab)/hip)!=0): red=0+abs(1-abs(math.asin((ab)/hip)/math.asin((au)/hipe))); red+=abs(1-abs(math.asin((au)/hipe))/math.asin((ab)/hip)); if (x[me]!=0)and(y[pe]!=0)and(x[ml]!=0)and(y[pl]!=0): ret=0+abs(1-abs((0.00+abs(float(y[pl])/float(y[pe])))/(0.00+abs(float(x[ml])/float(x[me]))))); ret+=abs(1-abs((float(x[ml])/float(x[me]))/(float(y[pl])/float(y[pe])))); re=1+abs(1-abs(float(abs(pe-pl)+1))/float(abs(me-ml)+1));rme=(y[pe]+0)+(y[pl]+0)+(x[me]+0)+(x[ml]+0);rmu+=rme; #print(rme); #print(rea,red,ret,re,(abs(e-l)+1),rme); if rea==0: rea=1; res+=abs(ret*red*red*rea)*(abs(e-l)+1)*rme*re;uj=uj+1; if uj==0: uj=1; return res/uj#/rmu#*rmu *(abs(x[ml]-y[pe])+0.01) *(abs(x[ml]-y[pe])+abs(x[me]-y[pl]))
Далее пример кода сравнения 1440 раз всех графиков с отступом от 0 до +1440 минут 2005 и 2006 года:
def dist(rast): statistic=0;statis=0; for i in range(1440): jk=i+000;x=harra1[jk];y=harra[jk+rast]; rt=stat(x,y); statistic+=rt; jk=i+000;x=harra[jk];y=harra1[jk+rast]; rt=stat(x,y); statis+=rt; return statis,statistic;d=[];for i in range(1440): print(i); d.append((dist(i),i)); print(d[-1]);d004=[];for i in range(len(d)): t4=int(d[i][0][0]+d[i][0][1])/2 d004.append(t4)import numpynumpy.argsort(d004)
Для данных за 2 суток я исследовал окрестности +1440 минут
>>> for i in range(len(d1)):
... print(d01[i],1440-i)
...
1529.0 1440
1499.0 1439
1510.0 1438
1530.5 1437
1439.5 1436
1534.5 1435
1516.0 1434
1473.0 1433
1503.0 1432
1486.0 1431
1514.5 1430
1507.5 1429
1525.0 1428
1460.0 1427
1558.5 1426
1562.5 1425
1495.0 1424
1511.0 1423
1527.5 1422
1507.0 1421
1506.0 1420
1498.0 1419
1517.0 1418
1471.0 1417
1478.0 1416
1588.0 1415
1513.5 1414
>>> import numpy
>>> numpy.argsort(d01)
array([ 38, 4, 78, 71, 59, 76, 13, 101, 74, 83, 41, 90,
23,
7, 24, 9, 33, 48, 100, 70, 29, 54, 44, 91, 63, 103,
93, 87, 42, 32, 16, 97, 57, 30, 21, 61, 39, 1, 64,
77, 99, 49, 52, 8, 79, 95, 66, 51, 45, 20, 40, 27,
60, 19, 11, 86, 34, 92, 2, 88, 17, 85, 26, 10, 6,
81, 22, 36, 73, 58, 56, 50, 12, 65, 69, 18, 46, 0,
68, 84, 3, 67, 96, 5, 28, 102, 31, 55, 47, 43, 75,
98, 94, 35, 89, 14, 37, 15, 72, 53, 82, 62, 25, 80])
>>> d01[38]
1419.0
>>> d01[4]
1439.5
И вот результат второй минимум выпадает на +1436 минуту d01[4] -4 значение
>>> for i in range(70):
... print(d04[340+i])
...
(1444, 340)
(1451, 341)
(1531, 342)
(1433, 343)
(1495, 344)
(1460, 345)
(1494, 346)
(1489, 347)
(1460, 348)
(1492, 349)
(1491, 350)
(1499, 351)
(1477, 352)
(1499, 353)
(1431, 354)
(1520, 355)
(1481, 356)
(1454, 357)
(1530, 358)
(1496, 359)
(1443, 360)
(1434, 361)
(1450, 362)
(1470, 363)
(1493, 364)
(1485, 365)
(1428, 366)
(1491, 367)
(1516, 368)
(1414, 369)
(1521, 370)
(1452, 371)
(1455, 372)
(1384, 373)
(1525, 374)
(1395, 375)
(1436, 376)
(1430, 377)
(1478, 378)
(1457, 379)
(1444, 380)
(1429, 381)
(1476, 382)
(1608, 383)
(1499, 384)
(1439, 385)
(1467, 386)
(1466, 387)
(1516, 388)
(1471, 389)
(1471, 390)
(1456, 391)
(1472, 392)
(1442, 393)
(1435, 394)
(1495, 395)
(1476, 396)
(1484, 397)
(1474, 398)
(1489, 399)
(1422, 400)
(1440, 401)
(1477, 402)
(1489, 403)
(1494, 404)
(1472, 405)
(1541, 406)
(1469, 407)
(1522, 408)
(1444, 409)
Далее для звездного года. Тут надо определить что алгоритм сравнивает картинки коллективного взаимодействия, поэтому гистограмный метод находит 368 как абсолютный минимум. Здесь же 373 минута третья от конца 369 минута тоже локальный минимум и 348 минута.
Всего рассчитано до +1440 минут двух дней в 2005 и 2006 годах.
>>> sor=numpy.argsort(d004)
>>> sor[:90] # всего до sor[:1440]
array([ 137, 140, 373, 583, 2, 47, 375, 266, 709, 40,
439,
29, 751, 226, 245, 433, 217, 789, 437, 524, 659, 369,
124, 919, 954, 544, 76, 256, 103, 991, 419, 106, 323,
400, 606, 626, 120, 99, 259, 321, 11, 1369, 235, 366,
158, 1258, 240, 381, 176, 118, 491, 377, 207, 354, 123,
825, 415, 167, 897, 588, 865, 46, 1129, 343, 60, 97,
842, 361, 211, 394, 628, 435, 1218, 878, 268, 376, 153,
165, 880, 494, 643, 609, 1439, 301, 85, 107, 252, 385,
173, 455])
