В этой статье будет приведен очень короткий код поиска периодов тропического года, звездного года и звездных суток в данных радиоактивного распада плутония.
Вот здесь на 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 math
import random
from statistics import mean
koeff=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]/average1
for jk in range(int(len(matrixd)/koeff)):
averaged=mean(harrad[jk])
for mk in range(koeff): harrad[jk][mk]=harrad[jk][mk]/averaged
Далее привожу функцию сравнения:
import numpy
import numpy as np
from numpy.fft import fft, fftfreq
import matplotlib.pyplot as plt
from numpy import linspace, loadtxt, ones, convolve
from numpy import array, zeros, argmin, inf
from numpy.linalg import norm
def _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 numpy
numpy.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])
v1000