Spiral of Theodorus…

Settembre 10, 2021·Biagio

The spiral of Theodorus is a spiral composed of right triangles. Hundreds of years ago, Theodorus of Cyrene constructed continuous right triangles and got a beautiful spiral. He used that spiral to prove that all non-square integers from 3–17 are irrational.

How would you plot this spiral? At each step, you need to draw a segment of length 1, perpendicular to the hypotenuse of the previous triangle. There are two perpendicular directions, and you want to choose the one that moves counterclockwise.

If we step outside the xy plane, we can compute the cross product of the unit vector in the z direction with the vector (x, y). The cross product will be perpendicular to both, and by the right-hand rule, it will point in the counterclockwise direction.

The cross product of (0, 0, 1) and (x, y, 0) is (-y, x, 0), so the direction we want to go in the xy plane is (-y, x). We divide this vector by its length to get a vector of length 1, then add it to our previous point.

Here is a code written in Python to plot the spiral

import matplotlib.pyplot as plt
    
def vertex(x, y):
    h = (x**2 + y**2)**0.5
    return (x - y/h, y + x/h)
    
plt.axes().set_aspect(1)
plt.axis('off')
    
# base of the first triangle
plt.plot([0, 1], [0, 0])
    
N = 17
x_old, y_old = 1, 0

for n in range(1, N):
    x_new, y_new = vertex(x_old, y_old)
    # draw short side
    plt.plot([x_old, x_new], [y_old, y_new])
    # draw hypotenuse
    plt.plot([0, x_new], [0, y_new])
    x_old, y_old = x_new, y_new

plt.show()

How to generate and read QR Code in Python

Settembre 4, 2021·Biagio

QR code is a type of matrix barcode that is machine readable optical label which contains information about the item to which it is attached. In practice, QR codes often contain data for a locator, identifier, or tracker that points to a website or application, etc.

Problem Statement :

Generate and read QR codes in Python using qrcode and OpenCV libraries

Installing required dependencies:

pyqrcode module is a QR code generator. The module automates most of the building process for creating QR codes. This module attempts to follow the QR code standard as closely as possible. The terminology and the encoding used in pyqrcode come directly from the standard.

pip install pyqrcode

Install an additional module pypng to save image in png format:

pip install pypng

Import Libraries

import pyqrcode
import png
from pyqrcode import QRCode
import cv2
import numpy as np

Create QR Code:

# OUTPUT SECTION

# String which represents the QR code
s = "http://www.raucci.net"

# output file name

filename = "qrcode.png"

# Generate QR Code

img = pyqrcode.create (s)

# Create and save the svg file naming "brqr.svg"

img.svg("brqr.svg", scale=8)

# Create and save the svg file naming "brqr.png"

img.png("brqr.png", scale=6)

Read QR Code

Here we will be using OpenCV for that, as it is popular and easy to integrate with the webcam or any video.

# INPUT SECTION

# read the QRCODE image

img = cv2.imread('brqr.png')

# initialize the cv2 QRCode detector
detector = cv2.QRCodeDetector()

val, pts, st_code=detector.detectAndDecode(img)

print(val)

Here the output

Estimating Pi…

Settembre 2, 2021·Biagio

We will be exploring a very interesting, and simple for that matter, application of statistics to help us estimate the value of Pi.

For this method, we will be imagining a simple scenario. Imagine for a moment, that we have the unit circle, inside a square.

By having the unit circle, we immediately figure out that the area of the square will be four since the radius of the circle is defined at one, which means that our square will have sides with a value of two. Now here’s where things get interesting, if we take the ratio of both of areas, we end up getting the following:

$\frac{{{A_c}}}{{{A_s}}} = \frac{{\pi {r^2}}}{{{{\left( {2r} \right)}^2}}} = \frac{\pi }{4}$

Both of these geometric figures end up having a ratio of pi over four between them, which is an important value for our next step in which we use a bit of imagination.

from “Ordine e disordine” by Luciano De Crescenzo

For a moment, imagine that you have a circle inside a square on the ground; suppose it starts raining. Some drops will most likely fall inside the circle and others will likely fall inside the square but outside the circle. Using this concept is how we will code our estimator pi, throwing some random numbers into the unit circle equation, as shown below:

$x^2+y^2=1$

Furthermore, taking a ratio of throws that landed inside our circle and the total number of throws, we can then formulate the following:

$\frac{{N{\text{of throws inside circle}}}}{{N {\text{of total throws}}}} = \frac{{{N_C}}}{{{N_T}}}$

And by combining our ratio between the unit circle with the square with this new equation, we can assume the next equation

$\frac{{{N_C}}}{{{N_T}}} \simeq \frac{\pi }{4} \Leftrightarrow \frac{{4{N_C}}}{{{N_T}}} \simeq \pi$

With this equation, we can finally start coding our estimator and see how close we can get to pi’s actual value.

The code

import matplotlib.pyplot as plt
import numpy as np
import time as t
from progress.bar import Bar


tin = t.time()


# numeri di punti della simulazione
n=80000

# Vettore coordinate  x e y dei punti casuali
x = np.random.rand(n)
y = np.random.rand(n)

Pi_Greco = np.zeros(n)
#Vettore distanza

d= (x**2 + y**2)**(1/2)

Ps = 0
Pq = 0

bar = Bar('Processing', max=n)

for i in range(n):
	Pq = Pq + 1
	if d[i] < 1:
		Ps = Ps+1
	Pi_Greco [i] = 4*(Ps/Pq)
	bar.next()

bar.finish()

Pi_Greco_Reale = np.ones(n)*np.pi

tfin = t.time()

print('Valore di Pi Greco: ', 4*Ps/Pq)
print('elapsed time: ', round(tfin-tin, 3))

plt.figure(1)
plt.plot(Pi_Greco, 'red', label='estimate of Pi')
plt.plot(Pi_Greco_Reale, 'green', label='Exact value')

plt.xlabel('throws')
plt.ylabel('value')
plt.title('Monte Carlo simulation')


plt.legend()
plt.show()

Results for our pi estimation. Pi is represented by the green line, while red represents our estimations.

It’s quite interesting to see our estimation start with low accuracy but as we increase our attempts, we start to get a convergence on the the value of pi, as shown by our green line. Statistical methods like this, and other more complex versions, are nice tools to understand and experiment within the world of physics. I highly recommend taking a look into some Statistical Mechanic concepts to see the beauty behind the application of statistics and probability in physics, and maybe take some time play with these concepts in Python!

Python animations with Matplotlib

Settembre 1, 2021·Biagio

In Python, plotting graphs is straightforward — you can use powerful libraries like Matplotlib. It happens, however, that you need to visualize the trend over time of some variables – that is, you need to animate the graphs.

Luckily, it’s just as easy to create animations as it is to create plots with Matplotlib.

Matplotlib

Matplotlib – as you can read on the official site – is a comprehensive library for creating static, animated, and interactive visualizations in Python. You can plot interactive graphs, histograms, bar charts, and so on.

How to Install Matplotlib

Installing Matplotlib is simple. Just open up your terminal and run:

pip install matplotlib

Numpy

Also, if you don’t have numpy, please install it so you can follow the examples in this tutorial:

pip install numpy

How to Plot with Matplotlib

Even though this tutorial is about animations in Matplotlib, first let’s create a simple static graph of a sine wave:

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(0, 10, 0.1)
y = np.sin(x)

fig = plt.figure()
ax = plt.axes(xlim=(0, 10), ylim=(-1.1, 1.1))
diagram = plt.plot(x, y)
 
plt.show()

How to Animate with Matplotlib

To create an animation with Matplotlib you need to use Matplotlib’s animation framework’s FuncAnimation class.

For instance, let’s create an animation of a sine wave:

import numpy as np
from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation

fig = plt.figure()
ax = plt.axes(xlim=(0, 5), ylim=(-1.5, 1.5))
line, = ax.plot([], [], lw=2)

def init():
    line.set_data([], [])
    return line,

def animate(i):
    x = np.linspace(0, 5, 500)
    y = np.sin(2 * np.pi * (x + 0.02 * i))
    line.set_data(x, y)
    return line,

sinewawe_animation = FuncAnimation(fig, animate, init_func=init, frames=200, interval=20, blit=True)
sinewawe_animation.save("Animation.gif")

plt.show()

We have:

Let’s then go through the code above in a bit more detail to better understand how animations work with Matplotlib.

Lines 1–3

Here you add the required libraries. In particular, we add the FuncAnimation class that can be used to create an animation for you.

Lines 5–7

fig = plt.figure()
ax = plt.axes(xlim=(0, 5), ylim=(-1.5, 1.5))
line, = ax.plot([], [], lw=2)

Here you first create an empty window for the animation figure. Then you create an empty line object. This line is later modified to form the animation.

Lines 9–11

def init():
    line.set_data([], [])
    return line,

Here you create an init() function that sets the initial state for the animation.

Lines 13–17

You then create an animate() function. This is the function that gives rise to the sine wave. It takes the frame number i as its argument, then it creates a sine wave that is shifted according to the frame number (the bigger it is, the more the wave is shifted). Finally, it returns the updated line object. Now the animation framework updates the graph based on how the line has changed.

def animate(i):
    x = np.linspace(0, 5, 500)
    y = np.sin(2 * np.pi * (x + 0.02 * i))
    line.set_data(x, y)
    return line,

Line 19

sinewawe_animation = FuncAnimation(fig, animate, init_func=init, frames=200, interval=20, blit=True)

This line of code puts it all together and creates the actual animation object. It simply:

Renders an animation to the figure fig by repeatedly calling the animate() function starting from the initial state defined by init()
The number of frames rendered to “one round of animation” is 200.
A delay between two frames is 20 milliseconds (1000ms / 20ms = 50 FPS).
(The blit=True makes sure only changed pieces of the plot are drawn, to improve the efficiency)

Line 21

sinewawe_animation.save("Animation.gif")

This piece of code is used to generate an animated gif (the same one I used in this tutorial to show you the animation)

Line 22

You guessed it — this line shows the animation.

A (apparently) simple problem …

Agosto 26, 2021·Biagio

So when I say a simple mathematical problem most would think that I am kidding but I am not, there are many unsolved mathematical problems in the world but this is so simple but yet unsolved.

The problem is called “3n + 1 problem” or “Collatz conjecture”. To understand the problem first we need to understand what it is, so basically just pick natural number if the number is odd then we mutiple the number with 3 and add 1, if the number is even we divide it by 2. We apply these conditions to resultant value. in other words:

In modular arithmetic notation, define the function f as follows:

$f(n)={\begin{cases}{\frac {n}{2}}&{\text{if }}n\equiv 0{\pmod {2}}\\[4px]3n+1&{\text{if }}n\equiv 1{\pmod {2}}.\end{cases}}$

Now form a sequence by performing this operation repeatedly, beginning with any positive integer, and taking the result at each step as the input at the next.

The Collatz conjecture is: This process will eventually reach the number 1, regardless of which positive integer is chosen initially.

To demonstrate the problem let’s consider a number 5; since it is odd we apply 3n+1:

$3 \cdot 5 + 1 = 16\left( {{\text{even}}} \right);16/2 = 8;8/2 = 4;4/2 = 2;2/2 = 1$

So we get a value of one but we apply conditions fruther we will be stuck in loop which is 4, 2, 1.

This video sums up the problem well

If the conjecture is false, it can only be because there is some starting number which gives rise to a sequence that does not contain 1. Such a sequence would either enter a repeating cycle that excludes 1, or increase without bound. No such sequence has been found.

Since I am from engineering background here is the Python code for Collatz conjecture:

def collatz(n):
    while n > 1:
        print(n, end=' ')
        if (n % 2):
            # n is odd
            n = 3*n + 1
        else:
            # n is even
            n = n//2
    print(1, end='')
 
 
n = int(input('Enter n: '))
print('Sequence: ', end='')
collatz(n)

The above code is the demonstration of Collatz conjecture…

Punto fisso

Novembre 26, 2020·Biagio

Un punto fisso per una funzione definita da un insieme in sé è un elemento coincidente con la sua immagine. Siano $f:A \to A$ e $x \in A$ . $x$ è un punto fisso per $f$ se $$x=f(x)$$ Si tratta di un punto che la funzione mappa in sé stesso.

Con l’ausilio di Python, valutiamo il punto fisso della funzione $y=\sqrt{x}$ .


import matplotlib.pyplot as plt
from numpy import array,linspace,sqrt,sin
from numpy.linalg import norm

def fixedp(f,x0,tol=10e-5,maxiter=100):
 """ Fixed point algorithm """
 e = 1
 itr = 0
 xp = []
 while(e > tol and itr < maxiter):
    x = f(x0)      # fixed point equation
    e = norm(x0-x) # error at the current step
    x0 = x
    xp.append(x0)  # save the solution of the current step
    itr = itr + 1
 return x,xp

f = lambda x : sqrt(x)
x_start = .5
xf,xp = fixedp(f,x_start)
x = linspace(0,2,1000)
y = f(x)
plt.plot(x,y,xp,f(xp),'bo',
     x_start,f(x_start),'ro',xf,f(xf),'go',x,x,'k')
stringa = "Fixed Point: " + str(round(xf, 4))

Exponential sums make interesting pictures

Novembre 25, 2020·Biagio

Exponential sums are a specialized area of math that studies series with terms that are complex exponentials.

$\sum\limits_{n=1}^{N}{e{}^{2\pi if(n)}}$

Exponential sums also make pretty pictures. If you make a scatter plot of the sequence of partial sums you can get surprising shapes. This is related to the trickiness of estimating such sums: the partial sums don’t simply monotonically converge to a limit. By the plot of an exponential sum we mean the sequence of partial sums, plotted in the complex plane, with successive points joined by straight line segments. That is, we start at the origin; draw a line interval corresponding to the first term of the sum; from the end of this interval draw another, corresponding to the second term of the sum; and so on.

Following this article, we playing around with polynomials with dates in the denominator. If we take that suggestion, with

$f(n)=\frac{n}{dd}+\frac{{n}^{2}}{mm}+\frac{{n}^{3}}{yy}$

and with today’s date, we get the curve below:

Here’s the python code that produced the image.

import matplotlib.pyplot as plt 
from numpy import array, pi, exp, log 
N = 20000 
def f(n): 
    return n/25 + n**2/11+ n**3/20 
z = array( [exp( 2*pi*1j*f(n) ) for n in range(0, N)] ) 
z = z.cumsum() 

plt.plot(z.real, z.imag, color='#333399') 
plt.axes().set_aspect(1) 
plt.show()

An interesting picture is the exponential sum with f(n)=(log n)⁴ and N=5000. The graph was dubbed “the Loch Ness monster” by John Loxton in a 1981 article.