Upto date with main repo

Introductory/mandelbroat/compare_languages/display.py

@@ -0,0 +1,18 @@
+#!/usr/bin/env python
+import sys
+import numpy as np
+import matplotlib.cm as cm
+from matplotlib.colors import LogNorm
+import pylab as pl
+
+filename = 'mand.dat'
+if len(sys.argv)>1:
+    filename = sys.argv[1]
+
+X,Y,Z = np.loadtxt(filename,unpack=True)
+ext=[X[0],X[-1],Y[0],Y[-1]]
+N = int(np.sqrt(len(Z)))
+dat = Z.reshape(N,N)
+print(ext)
+pl.imshow(dat.T,extent=ext,cmap=cm.hot,norm=LogNorm())   # pylab's function for displaying 2D image
+pl.show()

Introductory/mandelbroat/compare_languages/makefile

@@ -1,13 +1,14 @@
-CC = g++-10
+CC = g++-12 
 F90 = gfortran
+CFLAGS = -O3 -fopenmp
 
 all : mandc mandf
 
 mandc : mandc.cc
-	$(CC) -O3 -fopenmp -o $@ $<
+	$(CC) $(CFLAGS) -o $@ $<
 
 mandf : mandf.f90
-	$(F90) -O3 -fopenmp -o $@ $<
+	$(F90) $(CFLAGS) -o $@ $<
 
 clean :
 	rm -f mandc mandf

Introductory/mandelbroat/compare_languages/mandc.cc

@@ -1,14 +1,15 @@
-#include <iostream>
-#include <complex>
-#include <ctime>
-#include <vector>
-#include <omp.h>
-using namespace std;
+#include <iostream> // for printing to stdout
+#include <complex>  // for complex numbers
+#include <ctime>    // for measuring time
+#include <vector>   // for using array/vector
+#include <omp.h>    // openMP and measuing time with openMP
+using namespace std;// everything from standard template library will be in global scope
+// std::cout and std::endl can be called cout, endl, etc...
 
 int Mandelb(const complex<double>& z0, int max_steps)
 {
-  complex<double> z=0;
-  for (int i=0; i<max_steps; i++){
+  complex<double> z=0; // specify type where we start using it. (still need to specify type not like python with auto-typing)
+  for (int i=0; i<max_steps; i++){ // specify int type inside the loop, so it is local to the loop. useful for reducing bugs.
     if (abs(z)>2.) return i;
     z = z*z + z0;
   }
@@ -17,35 +18,38 @@ int Mandelb(const complex<double>& z0, int max_steps)
 
 int main()
 {
-  const int Nx = 1000;
+  const int Nx = 1000; // constants are named cost rather than parameter
   const int Ny = 1000;
-  int max_steps = 1000;
-  double ext[]={-2,1,-1,1};
+  int max_steps = 1000; // allowed to change later.
+  double ext[]={-2,1,-1,1}; // this specifies fixed array of four numbers and initializes it.
 
-  vector<int> mand(Nx*Ny);
-  clock_t startTimec = clock();
-  double start = omp_get_wtime();
+  vector<int> mand(Nx*Ny);  // allocating one dimensional array (vector) of size Nx*Ny
+  // multidimensional dynamic arrays are not standard in C++. One can use pointers to allocate/deallocate memory,
+  // and write one's own class interface for that. Or one has to use extension of C++ (blitz++ is excellent).
+  // Unfortunately the standard C++ still does not support standard class for that.
+  clock_t startTimec = clock();  // cpu time at the start 
+  double start = omp_get_wtime(); // wall time at the start
 
   #pragma omp parallel for
   for (int i=0; i<Nx; i++){
     for (int j=0; j<Ny; j++){
-      double x = ext[0] + (ext[1]-ext[0])*i/(Nx-1.);
-      double y = ext[2] + (ext[3]-ext[2])*j/(Ny-1.);
-      mand[i*Ny+j] = Mandelb(complex<double>(x,y), max_steps);
+      double x = ext[0] + (ext[1]-ext[0])*i/(Nx-1.); // x in the interval ext[0]...ext[1]
+      double y = ext[2] + (ext[3]-ext[2])*j/(Ny-1.); // y in the interval ext[2]...ext[3]
+      mand[i*Ny+j] = Mandelb(complex<double>(x,y), max_steps); // storing values in 2D array using 1D array
     }
   }
 
   clock_t endTimec = clock();
-  double diffc = double(endTimec-startTimec)/CLOCKS_PER_SEC;
-  double diff = omp_get_wtime()-start;
+  double diffc = double(endTimec-startTimec)/CLOCKS_PER_SEC; // how to get seconds from cpu time
+  double diff = omp_get_wtime()-start;                       // openMP time is already in seconds
 
-  clog<<"clock time : "<<diffc<<"s"<<" with wall time="<<diff<<"s "<<endl;
+  clog<<"clock time : "<<diffc<<"s"<<" with wall time="<<diff<<"s "<<endl; // printout of time
 
   for (int i=0; i<Nx; i++){
     for (int j=0; j<Ny; j++){
       double x = ext[0] + (ext[1]-ext[0])*i/(Nx-1.);
       double y = ext[2] + (ext[3]-ext[2])*j/(Ny-1.);
-      cout<<x<<" "<<y<<" "<< 1./mand[i*Ny+j] << endl;
+      cout<<x<<" "<<y<<" "<< 1./mand[i*Ny+j] << endl; // prinout of mandelbrot set
     }
   }
 }

Introductory/mandelbroat/compare_languages/mandf.f90

@@ -7,8 +7,8 @@ INTEGER Function Mandelb(z0, max_steps)
   INTEGER    :: i
   z=0.
   do i=1,max_steps
-     if (abs(z)>2.) then
-        Mandelb = i-1
+     if (abs(z)>2.) then ! |f(z)|>2 we know it will explode for large n, hence not part of the set.
+        Mandelb = i-1      ! return how many iterations it took to know this is not part of the set.
         return 
      end if
      z = z*z + z0
@@ -31,16 +31,16 @@ program mand
   INTEGER, parameter :: max_steps = 1000
   REAL*8  :: ext(4) = (/-2., 1., -1., 1./) ! The limits of plotting
   REAL*8  :: mande(Nx,Ny)
-  REAL    :: start, finish, startw, finishw
+  REAL*8  :: start, finish, startw, finishw
 
   call cpu_time(start)
   startw  = OMP_get_wtime()
 
   !$OMP PARALLEL DO  PRIVATE(j,x,y,z0)
   do i=1,Nx
      do j=1,Ny
-        x = ext(1) + (ext(2)-ext(1))*(i-1.)/(Nx-1.)
-        y = ext(3) + (ext(4)-ext(3))*(j-1.)/(Ny-1.)
+        x = ext(1) + (ext(2)-ext(1))*(i-1.)/(Nx-1.) ! x \in [ext(1)...ext(2)]
+        y = ext(3) + (ext(4)-ext(3))*(j-1.)/(Ny-1.) ! y \in [ext(3)...ext(4)]
         z0 = dcmplx(x,y)
         mande(i,j) = Mandelb(z0, max_steps)
      enddo

Introductory/mandelbroat/compare_languages/manp.py

@@ -1,11 +1,12 @@
 #!/usr/bin/env python
-
-from scipy import *
-from pylab import *
-import time
-
+from scipy import *    # scientific python
+from pylab import *    # plotting library
+import time            # timeing
+from numba import jit  # This is the new line with numba
+
+@jit(nopython=True)    # This is the second new line with numba
 def Mand(z0, max_steps):
-    z = 0j
+    z = 0j  # no need to specify type. To initialize to complex number, just assign 0j==i*0
     for itr in range(max_steps):
         if abs(z)>2.:
             return itr
@@ -18,17 +19,18 @@ def Mand(z0, max_steps):
     Ny = 1000
     max_steps = 1000 #50
 
-    ext = [-2,1,-1,1]
+    ext = [-2,1,-1,1]  # no need to specify this is a python list, just initialize it
     t0 = time.time()
 
-    data = zeros( (Nx,Ny) )
+    data = zeros( (Nx,Ny) ) # initialize a 2D dynamic array (something C++ is still lacking)
 
     for i in range(Nx):
         for j in range(Ny):
             x = ext[0] + (ext[1]-ext[0])*i/(Nx-1.)
             y = ext[2] + (ext[3]-ext[2])*j/(Ny-1.)
-            data[i,j] = Mand(x + y*1j, max_steps)
+            data[i,j] = Mand(x + y*1j, max_steps)  # creating complex number of the fly
+
     print ('clock time: '+str( time.time()-t0) )
-    imshow(transpose(1./data), extent=ext)
-    show()
+    imshow(transpose(1./data), extent=ext)   # pylab's function for displaying 2D image
+    show()                                   # showing the image.