<jupyter_start><jupyter_text>![Matplotlib](images/matplotlib_logo-s.png) Matplotlib Let's take a quick look at the matplotlib essentials.First let's import the neccesary modules:<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np<jupyter_output><empty_output><jupyter_text>And now let's go with some examples.  1D plot A tangent between 0 and 4$\pi$.$$y = \tan(x), \; x\in[0., 4\pi]$$<jupyter_code>x = np.linspace( 0., 4.*np.pi, num=200 )
y = np.tan(x)

plt.plot(x, y)

plt.xlabel(r'$x$')        #using LaTeX syntax
plt.ylabel(r'$\tan(x)$')
plt.grid(True)
#plt.savefig('tan.png')<jupyter_output><empty_output><jupyter_text>1D polar plot  $\theta(\rho) = 2\pi \rho, \; \rho \in [0, 2]$<jupyter_code>r = np.arange(0, 2., 0.01)
theta = 2 * np.pi * r

plt.polar(theta, r, color='r', linewidth=3)
plt.grid(True)

plt.title("A line plot on a polar axis")<jupyter_output><empty_output><jupyter_text>Exercise Plot the function $$y = \frac{\sin(x)}{x}, \; x\in[0.1, 6\pi]$$and try to * set the figure *legend** change the *line style*hint: checkout the [matplotlib examples](http://matplotlib.org/examples/index.html)<jupyter_code># %load files/sol_sinc.py<jupyter_output><empty_output><jupyter_text>![SciPy](images/scipy_logo.png) The SciPy framework builds on top NumPy and provides a large number of numerical algorithms for working with data. Some of the topics that SciPy covers are:* Special functions ([scipy.special](http://docs.scipy.org/doc/scipy/reference/special.html))* Integration/ODEs ([scipy.integrate](http://docs.scipy.org/doc/scipy/reference/integrate.html))* Optimization ([scipy.optimize](http://docs.scipy.org/doc/scipy/reference/optimize.html))* Interpolation ([scipy.interpolate](http://docs.scipy.org/doc/scipy/reference/interpolate.html))* Fourier Transforms ([scipy.fftpack](http://docs.scipy.org/doc/scipy/reference/fftpack.html))* Signal Processing ([scipy.signal](http://docs.scipy.org/doc/scipy/reference/signal.html))* Linear Algebra ([scipy.linalg](http://docs.scipy.org/doc/scipy/reference/linalg.html))* Sparse Eigenvalue Problems ([scipy.sparse](http://docs.scipy.org/doc/scipy/reference/sparse.html))* Statistics ([scipy.stats](http://docs.scipy.org/doc/scipy/reference/stats.html))* Multi-dimensional image processing ([scipy.ndimage](http://docs.scipy.org/doc/scipy/reference/ndimage.html))* File IO ([scipy.io](http://docs.scipy.org/doc/scipy/reference/io.html))In many cases, you will want to import specific names from `scipy` subpackages. However, as a start, it is helpful to do the following import:<jupyter_code>import scipy as sp<jupyter_output><empty_output><jupyter_text>Approach One of the most important skills in data science is to be able to find Python functions and classes in a module and learn how to use them yourself. Here are some recommended steps on how to go about this:* Find the online documentation for the package you are using.* Try to find the subpackage or even the function that looks like will do the job.* Import the module, function or class and use tab completion and `?` to explore it.* Try using the function or class for an extremely simple case where you know the answer.* Then try using for your real problem.<jupyter_code>%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Data FittingLet's randomize a exponential decay function and try to do a linear fit.We have$$D(t) = N e^{(-t/\tau)}$$ * Define the function<jupyter_code>def D(t, N, tau):
    return N * np.exp(-t/tau)<jupyter_output><empty_output><jupyter_text>* randomize data<jupyter_code>ndata = 50
tdata = np.linspace(0., 5., ndata)
data_pure = D(tdata, N=7., tau=3.)
data_measure = data_pure + 0.3 * np.random.normal(size=ndata)<jupyter_output><empty_output><jupyter_text>* plot<jupyter_code>plt.plot(tdata, data_pure)
plt.plot(tdata, data_measure)<jupyter_output><empty_output><jupyter_text>* Look for fit scipy module and fit<jupyter_code>from scipy.optimize import curve_fit

popt, pcov = curve_fit(D, tdata, data_measure)

print (popt)
N_fit, tau_fit = popt
data_fit = D(tdata, N=N_fit, tau=tau_fit)<jupyter_output>[ 7.14801215  2.89312668]<jupyter_text>* plot result<jupyter_code>plt.plot(tdata, data_pure)
plt.plot(tdata, data_measure)
plt.plot(tdata, data_fit, label="fit")

plt.legend()<jupyter_output><empty_output><jupyter_text>Exercise: simple integration  Calculate $$I = \int_0^\infty e^{-x} dx$$hint: [GIYF](https://www.google.es/search?q=scipy%20integrate%20function&rct=j)<jupyter_code># %load files/sol_integ.py
#help(integrate)<jupyter_output><empty_output><jupyter_text>Real life example: solar activity period Adapted from [Python for scientific use, Part II: Data analysis](http://linuxgazette.net/115/andreasen.html) by <NAME>. See also a more modern approach to the subject using [Pandas](http://pandas.pydata.org/),  the Python Data Analysis Library, on [this Josh Hemann notebook](http://nbviewer.ipython.org/gist/jhemann/4569783)In the next example we will use the Fast Fourier Transform (FFT) in order to transform time-dependent data into the frequency domain. By doing so, it is possible to analyse if any predominant frequencies exists - i.e. if there is any periodicity in the data. The example data we will use is the sunspot activity measurements from year 1700 to year 2014 provided by the Sunspot index and Long-term Solar Observation. The data set is the yearly mean total sunspot observations available via HTTP.The datafile is also included in `files/SN_y_tot_V2.0.txt`. The sunspot data have been used for illustrating the power of FFT with respect to finding a periodicity in sunspot activity in various computer languages, e.g. Matlab and BASIC.   Retrieve dataSolar spot data are available online from the [datafiles section](http://www.sidc.be/silso/datafiles) of the Sunspot index and Long-term Solar Observation webpage. We will use the yearly mean total sunspot number [plain text version.]( versionhttp://www.sidc.be/silso/DATA/SN_y_tot_V2.0.txt)A quick web search point us to the `urllib` python module and to a very simple retrival procedure.<jupyter_code>import urllib
url = "http://www.sidc.be/silso/DATA/SN_y_tot_V2.0.txt"
raw_data = urllib.request.urlopen(url)
solarspots = np.loadtxt(raw_data)
solarspots
years = solarspots[:,0]
spots = solarspots[:,1]
plt.figure(figsize=(10,6))

plt.plot(years, spots)
plt.plot(years, spots, "ro", label="spots" )

plt.xlabel("year")
plt.ylabel("yearly sunspots")
plt.legend()
plt.grid()

#plt.savefig("yearly_sunsports.png", dpi=300)<jupyter_output><empty_output><jupyter_text>FFTLet's take a look at the [SciPy Fast Fourier Transform module manual](http://docs.scipy.org/doc/scipy/reference/tutorial/fftpack.html), and straight from there...<jupyter_code>from scipy.fftpack import fft, fftfreq
spots_fft = fft(spots)
plt.scatter(spots_fft.real, spots_fft.imag)

#plt.xlim(-5000,5000)
plt.xlabel("real(FFT)")
plt.ylabel("imag(FFT)")
plt.grid()<jupyter_output><empty_output><jupyter_text>In order to construct a periodogram, i.e. a graph of power vs. frequency, we first compute the power of the FFT signal which is simply the **FFT signal squared**. We only need the part of the signal ranging from zero to a frequency equal to the Nyquist frequency, which is equal to half the maximum frequency, since frequencies above the Nyquist frequency correspond to negative frequencies. The frequency range is calculated from 0-N/2 as N/(2T) where N is the number of samples and T is the sampling time.<jupyter_code>N = spots_fft.size
timestep = 1 # year
freq = fftfreq(N, d=timestep)   # 1/year
power = np.abs(spots_fft)   # proportional to this

freq = freq[1:N/2]
power = power[1:N/2]
plt.plot(freq, power)

plt.xlabel("Frequency [1/year]")
plt.ylabel(r"$|FFT|^2$")
plt.grid()<jupyter_output><empty_output><jupyter_text>Which in years becomes<jupyter_code>period = 1./freq
plt.plot(period, power)

plt.xlabel("Period [years]")
plt.ylabel(r"$|FFT|^2$")
plt.grid()

#plt.xlim( (0,30) )<jupyter_output><empty_output><jupyter_text>Solar spots period<jupyter_code>max_power_index = np.argmax(power)
max_power_period = period[max_power_index]
print ('Maximum solar spots period is %1.2f years' %(max_power_period) )<jupyter_output>Maximum solar spots period is 10.86 years<jupyter_text>Scipy direct method After doing the exercise, while looking for information, I found that scipy allows us to directly compute the periodgram with [`scipy.signal.periodgram`](http://scipy.github.io/devdocs/generated/scipy.signal.periodogram.html)<jupyter_code>from scipy import signal
frec=1.  # sampling frequency: 1 year
f, Pxx_spectrum = signal.periodogram(spots, frec, scaling='spectrum')
plt.plot(f, Pxx_spectrum)
plt.xlim(.0,.5)
max_power_index = np.argmax(Pxx_spectrum)
max_power_period = 1./f[max_power_index]
print ('Maximum solar spots period is %1.2f years' %(max_power_period) )<jupyter_output>Maximum solar spots period is 10.86 years
{'fraction_non_alphanumeric': 0.12359048788656916, 'fraction_numerical': 0.014290499051021548, 'mean_word_length': 5.843391902215432, 'pattern_counts': {'":': 0, '<': 71, '<?xml version=': 0, '>': 71, 'https://': 1, 'lorem ipsum': 0, 'www.': 4, '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': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '21197', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 3289, 'n_tokens_neox': 3088, 'n_words': 860, 'path': 'scipy-tutorial.ipynb', 'repo_name': 'iamc/Pyss15-scientific-python-tutorial', 'size': '8967'}