<jupyter_start><jupyter_text>Ξεκινάμε με μια απλη επανάληψη του πως μπορούμε να κατασκευάσουμε ένα γραμμικό χώρο και να σχεδιάσουμε τη συνάρτηση ημίτονο απο 0 έως 2π.<jupyter_code>xx=np.linspace(0.,2.*np.pi,200) #Δημιουργούμε ένα ομοιογενή χώρο 200 στοιχείων απο 0 έως 2π
plt.plot(xx,np.sin(xx)) #Εκτυπώνουμε τα σημεία (x,sin(x)) χρησιμοποιώντας την συνάρτηση sin() της Numpy
plt.xlabel('X')
plt.ylabel('Y')<jupyter_output><empty_output><jupyter_text>Παρακαάτω θα προσπαθήσουμε να εκτιμήσουμε τη περιοδικότητα της ηλιακής δραστηριότητας μέσω του αριθμού του κηλίδων $R$.Ο αριθμός των κηλίδων δίνεται από τη σχέση $ R= N_s + 10  N_g$ όπου $N_s$ ο αριθμός των κηλίδων και $N_g$ ο αριθμών των ομάδων (groups) που διακρίνονται πάνω στον ηλιακό δίσκο.Τα δεδομένα όπως φαίνονται και παρακάτω μέσα από την εντολή **head** του linux bash (μπορούμε να εκτελέσουμε μια εντολή του bash (για linux) ή του command line (για windows) μέσω του ! (θαυμαστικό)) εκτείνονται από το 1818 σε 8 στήλες.Στήλη 1-3 έτος, μήνας, μέραΣτήλη 4 ημερομηνία σε έτη (δεκαδική μορφή)Στήλη 5: Ημερήσιος αριθμός κηλίδων. Στη περίπτωση που δεν έχουμε διαθέσιμες μετρήσεις έχουμε τον αριθμό -1Στήλη 6: Τυπική διακύμανση του αριθμού των κηλίδωνΣτήλη 7: Αριθμός παρατηρήσεων για τη συγκεκριμένη ήμερα.Στήλη 8: Κενό αν η παρατήρηση θεωρείται έγκυρη, '*' αν υπάρχει περίπτωση επανεξέτασης.<jupyter_code>!head SN_d_tot_V2.0.csv
A=np.loadtxt('SN_d_tot_V2.0.csv',delimiter=';') #Εισάγουμε τα δεδομένα μέσα σε ένα δισδιάστατο πίνακα (A)
t=A[:,3] #χρονος (δεκαδικη μορφη σε ετη)
R=A[:,4] #αριθμος κηλιδων
N=R.shape[0] #αριθμος ολων των παρατηρησεων<jupyter_output><empty_output><jupyter_text>Για λογους ευκολιας θα θεωρησουμε οτι τις ημερες οπου δεν ειχαμε παρατηρησεις δεν ειχαμε κηλιδες, αρα αναζητουμε τα σημεια του πινακα των παρατηρησεων R οπου υπαρχει η τιμη -1. Στα σημεια αυτα θα βαλουμε τη τιμη 0.<jupyter_code>R[R==-1]=0.
plt.plot(t,R)
#plt.xlim(1860,1861)<jupyter_output><empty_output><jupyter_text>Χρήση του αλγορίθμου FFT(https://en.wikipedia.org/wiki/Fast_Fourier_transform) Μέσω του αλγορίθμού FFT που εμπεριέχεται στη βιβλιοθήκη numpy.fft [1] θα υπολογίσουμε το διακριτό μετασχηματισμό fourier της χρονοσειράς $R$ (αριθμός κηλίδων) $${\displaystyle X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-i2\pi kn/N}\qquad k=0,\dots ,N-1.}$$[1] Cooley, <NAME>., and <NAME>, 1965, "An algorithm for the        machine calculation of complex Fourier series," *Math. Comput.*        19: 297-301.        https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm<jupyter_code>import numpy.fft as ff
dt=t[400]-t[399] #Χρονική διαφορά ανάμεσα σε δύο παρατηρήσεις (δεκαδική μορφή)
print('Χρονική διαφορά ανάμεσα στη παρατήρηση 400 και 399: {:.4f}'.format(dt)) #Εκτύπωσε τη (μεχρι 4 δεκαδικα)
freq=ff.fftfreq(N,d=dt) #Χώρος συχνοτήτων για τις συγκεκριμένες παρατηρήσεις μέσω της βιβλιοθηκης fft της numpy
fR=np.abs(ff.fft(R))/N  #Πλατος του σήματος για τις παραπάνω συχνότητες<jupyter_output>Χρονική διαφορά ανάμεσα στη παρατήρηση 400 και 399: 0.0020<jupyter_text>θα κρατήσουμε μόνο τις θετικές συχνότητες (ο μετασχηματισμός είναι συμμετρικός ως προς τη συχνότητα 0)<jupyter_code>F=freq[freq>0]
T=1/F #περιοδος
FR=fR[freq>0]
plt.plot(T,FR) #Διάγραμμα του πλάτους των συχνοτήτων για κάθε περίοδο
plt.xlabel('Period (yrs)')<jupyter_output><empty_output><jupyter_text>Αναζητούμε τη περίοδο όπου το μετασχηματισμένο σήμα είναι μέγιστο.<jupyter_code>T[FR.argmax()]<jupyter_output><empty_output><jupyter_text>Παρατηρούμε ότι η παραπάνω τιμή είναι διαφορετική από την γνωστή η οποία είναι περίπου 11 έτη. Ο λόγος για τον οποίο παρατηρούμε αυτή τη διαφορά ενδέχεται να βρίσκεται στην εκτίμηση των χρονικών διαφορών dt. Για να το εξετάσουμε δημιουργούμε μια λίστα με $N-1$ στοιχεία ή οποία θα περιλαμβάνει όλες τις χρονικές διαφορές $dt=t_i-t_{i-1}$<jupyter_code>dts=np.array([t[i]-t[i-1] for i  in range(1,N)])
plt.hist(dts)<jupyter_output><empty_output><jupyter_text>Παρατηρούμε ότι οντως οι χρονικές μεταβολές διαφέρουν. Αυτο μπορεί να οφείλεται στη δεκαδική μορφή των ημερομηνιων και να μπορει να διoρθωθει εκμεταλλευόμενοι ολόκληρη την ημερομηνια όπως δίνεται από τις στήλες 1-3.Για λόγους εξάσκησης όμως θα υποθέσουμε οτι δεν έχουμε αυτη τη δυνατότητα και απλά έχουμε δεδομένα όπου οι παρατηρήσεις παρθηκαν σε μη τακτά χρονικά διαστήματα. Για να λύσουμε αυτό το πρόβλημα θα αντικαταστήσουμε τα πραγματικά δεδομένα με ένα γραμμμικό συνεχές (linear interpolation), ενώνοντας τα σημεία μεταξύ τους με ευθείες.<jupyter_code>from scipy.interpolate import interp1d
Rp=interp1d(t,R) #Δημιουργία ενος interpolated "συνεχούς" σήματος με βάση τα δεδομένα (t,s)
tt=np.linspace(t[0],t[-1],N) #Δημιουργία ενός "συνεχούς" ισοχρονικών διαστημάτων
dt=tt[4]-tt[3]
print(dt)
freq=ff.fftfreq(N,d=dt)
fRp=np.abs(ff.fft(Rp(t)))/N
plt.plot(1./freq[freq>0],fRp[freq>0])
plt.xlim(0,20)
print(1./freq[fRp[freq>0].argmax()])<jupyter_output>0.00273788182881
11.7153963455<jupyter_text>Μέσω τη παραπάνω μεθόδου παρατηρούμε ότι η τιμή είναι η αναμένομενη, στα 11.7 έτη. Στη συνέχεια θα παρουσιάσουμε πως μπορούμε να παράξουμε τα ίδια αποτελεσματα μέσω της βιβλιοθήκης ανάλυσης δεδομένων ```pandas``` όπου και θα χρησιμοιήσουμε την ακριβή ημερομηνία κάθε παρατήρησης.<jupyter_code>import pandas as pd
df=pd.read_csv('SN_d_tot_V2.0.csv',delimiter=';',
               names=['Year','Month','Day','Date_dec','R','dR','n','valid'])
df['Date']=pd.to_datetime(df[['Year','Month','Day']])
df.set_index('Date')
#df
N=df.shape[0]
freq=ff.fftfreq(df.shape[0],d=1.) #dt= 1 day
fs=np.abs(ff.fft(df['R']))/N
period=(1./freq[freq>0])/365.
MaxFs=period[fs[freq>0].argmax()]
plt.plot(period,fs[freq>0]) 
plt.vlines(MaxFs,0,fs[freq>0].max(),linestyle='--',linewidth=0.75,label='{:.2f} Years'.format(MaxFs))
plt.xlim(0,20)
plt.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1127127474817645, 'fraction_numerical': 0.03004515456755818, 'mean_word_length': 6.355044699872286, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 46, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '277660', 'license': 'MIT', 'n_tokens_mistral': 4539, 'n_tokens_neox': 3035, 'n_words': 606, 'path': 'Sun Periodicity.ipynb', 'repo_name': 'UOAPythonWorkshop/WorkshopNotebooks', 'size': '5766'}