'''This program serves for estimating the fractal dimension of the point set F of Example 9.6.
'''

# from __future__ import print_function
import numpy as np

n = 10

N = np.zeros(n, dtype='int')
number_of_boxes = 4**np.arange(1,n+1)
box_size = 4.**-np.arange(1,n+1)

for j in range(n):
    grid = np.linspace(0,1,number_of_boxes[j]+1)   
    F = 1. / np.arange(1, number_of_boxes[j] + 2)
    
    '''We count how many boxes intersect with F = {0, 1, 1/2, 1/3, 1/4, ...}'''
    hist, bin_edges = np.histogram(F, bins = grid + 1e-15)
    N[j] = sum(hist>0)
    
'''Calculate box dimension d'''            
d = - (np.log(N[1:]) - np.log(N[:-1])) / (
                    np.log(box_size[1:]) - np.log(box_size[:-1]))

'''Approximate box dimension d'''
d_approximation = - np.log(N) / np.log(box_size)

print('Boxes \t\t Approximation of d \t d')
print('{0:7}'.format(number_of_boxes[0]), '\t', '{0:.12f}'.format(d_approximation[0]))
for j in range(1, n):
    print('{0:7}'.format(number_of_boxes[j]), '\t', 
          '{0:.12f}'.format(d_approximation[j]), '\t', 
          '{0:.12f}'.format(d[j-1]))