Simple Fortran kernel for Jupyter

Switching compilers

%comipler: compiler-name .file-extention

default compiler is intel fortran

In [1]:
use, intrinsic :: iso_fortran_env
print *, compiler_version()
print *, compiler_options()
end
 GCC version 9.1.0
 -mtune=generic -march=x86-64 -Wall

gfortran

In [2]:
%compiler: gfortran
use, intrinsic :: iso_fortran_env
print *, compiler_version()
print *, compiler_options()
end
 GCC version 9.1.0
 -mtune=generic -march=x86-64 -Wall

Change file extention

default file extention is .f90 (free format)

for FIXED FORMAT source use .f/.for

In [4]:
%compiler: gfortran .f
      PARAMETER(N = 100)
      OPEN(9, FILE = 'SINCOS.XY')
      PI = 4.0 * ATAN(1.0)
      DO 10 I = 0, N
        X = 2.0 * I * PI / N  
        WRITE(9, *) X, SIN(X), COS(X)
   10 CONTINUE
      STOP 'NORMAL TERMINATION'
      END
STOP NORMAL TERMINATION

Write text file

  • %writefile fn [-a]
  • %%writefile fn [-a]
  • -a append
In [5]:
%writefile: Laurence_Sterne.txt
O Cupid! Cupid! prince of Gods and men.”
write to file:Laurence_Sterne.txt
In [6]:
%%writefile: Laurence_Sterne.txt -a
Laurence Sterne. (17131768).  A Sentimental Journey through France and Italy.
append to file:Laurence_Sterne.txt

Compile C program

  • %compiler: gcc .c
In [9]:
%compiler: gcc .c
#include <stdio.h>
int main(void) {
    printf("Hello world\n");
    return 0;  
}
Hello world

Compile C++ program

  • %compiler: g++ .c++
In [7]:
%compiler: g++ .c++
#include <iostream>

using namespace std;

int main() {
    std::cout << "Hello World!" << std::endl;
    return 0;
}
Hello World!

Separate module compilation

%module: object-file-name

In [8]:
%module: m
%fcflags: -O2
!
! The Sieve of Eratosthenes
!
module m_eratos
    implicit none
contains
    integer function nprimes(n) 
        integer, intent(in) :: n
        logical, allocatable :: tab(:)
        integer :: i
        allocate(tab(2:n), source = .true.)
        do i = 2, int(sqrt(real(n)))
            if (tab(i)) tab(i**2::i) = .false.
        end do
        nprimes = count(tab)
    end function nprimes    
end module m_eratos
[gfort kernel] module objects created successfully: m.o

Compiler options

%fcflags: arg1 arg2 ... %ldflags: arg1 arg2 ...

gfortran fcflags source.f90 ldflags

to use module, put object-file-names to fcflags

In [10]:
%fcflags: m.o
program test
    use m_eratos
    implicit none
    integer :: i 
    do i = 1, 8
       write(*, *) i, 10**i, nprimes(10**i), nint(10**i / log(10.0**i))
       write(9, *) i, 10**i, nprimes(10**i), nint(10**i / log(10.0**i))
    end do
end program test
           1          10           4           4
           2         100          25          22
           3        1000         168         145
           4       10000        1229        1086
           5      100000        9592        8686
           6     1000000       78498       72382
           7    10000000      664579      620421
           8   100000000     5761455     5428681

Matplotlib Ready

%fig: args

import matplotlib.pyplot as plt
import numpy as np
_fig = plt.figure( args )

...

( user inputs are executed by Python's exec command )

...

plt.show()

note

I don't know why, but plt.figure() doesn't work properlly with exec/eval commands. So it is hard coded in the kernel.

plot number of primes $\pi(n)$

$$ \large\pi(n)\simeq{n \over {\log(n)}} $$
In [12]:
%fcflags: m.o -O2  
program prime
    use m_eratos
    implicit none
    integer :: i 
    open(9, file = 'primes.csv')
    write(9, '(5g0)') 'n', ',', 'primes', ',', 'n/log(n)'
    do i = 2, 501
      ! write(*, *) i, ',', nprimes(i)
       write(9, '(5g0)') i, ',', nprimes(i), ',', i / log(real(i)) 
    end do
    stop 'normal end'
end program prime
STOP normal end
In [13]:
%fig:
x, y1, y2 = np.loadtxt('primes.csv', unpack = True, delimiter=',', skiprows=1) 
plt.grid(which='both')
plt.plot(x, y1)
plt.plot(x, y2)
plt.title('Prime numbers', fontsize = 25)
plt.xlabel('n', fontsize = 20)
plt.ylabel('number of primes', fontsize = 20)
In [14]:
%fig:
x, y1, y2, y3 = np.loadtxt('fort.9', unpack = True)
plt.yscale('log')
plt.grid(which='both')
plt.plot(x, y2)
plt.plot(x, y3)
plt.title('Prime number theorem', fontsize = 25)
plt.xlabel('10^x', fontsize = 20)
plt.ylabel('number of primes', fontsize = 20)

subplot

sin cos curves

In [15]:
%fig:figsize=(6, 2)
x, y1, y2 = np.loadtxt('SINCOS.XY', unpack = True)
plt.subplot(1, 2, 1)
plt.grid(which='both')
plt.plot(x, y1)
plt.subplot(1, 2, 2)
plt.grid(which='both')
plt.plot(x, y2)

Pandas examples

DataFrame is AxesSubplot object.

Internal variable _fig is accessible and assignable (e.g. _fig = rhs). However multiple assignment is not possble.

In [16]:
%fig:
import pandas as pd
ax = _fig.add_subplot(1, 1, 1)

x = np.arange(100)
df = pd.DataFrame({'x': x, 'sin(x)': np.sin(x / 5)})
df.plot(ax=ax)
In [17]:
%fig:
from pandas import *
fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(9, 6))
_fig = fig
# Multiple assignment is not allowed for _fig.

df = DataFrame(np.random.randn(1000, 4))
df.plot(ax=axes[0], legend=False)
df.plot(ax=axes[1], legend=False, kind='hist')

BMP example

In [18]:
%module: BMP
    module m_bmp
      use, intrinsic :: iso_fortran_env
      implicit none
      type :: t_bmp_file_header
        sequence  
        integer(int16) :: bfType = transfer('BM', 0_int16) ! BitMap
        integer(int32) :: bfSize          ! file size in bytes
        integer(int16) :: bfReserved1 = 0 ! always 0
        integer(int16) :: bfReserved2 = 0 ! always 0
        integer(int32) :: bfOffBits
      end type t_bmp_file_header
      !
      type :: t_bmp_info_header
        sequence
        integer(int32) :: biSize          = Z'28' ! size of bmp_info_header ; 40bytes 
        integer(int32) :: biWidth
        integer(int32) :: biHeight
        integer(int16) :: biPlanes        = 1 ! always 1
        integer(int16) :: biBitCount
        integer(int32) :: biCompression   = 0 ! 0:nocompression, 1:8bitRLE, 2:4bitRLE, 3:bitfield
        integer(int32) :: biSizeImage
        integer(int32) :: biXPelsPerMeter = 3780 ! 96dpi
        integer(int32) :: biYPelsPerMeter = 3780 ! 96dpi 
        integer(int32) :: biClrUsed       = 0
        integer(int32) :: biClrImportant  = 0 
      end type t_bmp_info_header
      !
      type :: t_rgb
        sequence
        integer(int8) :: ib, ig, ir
      end type t_rgb  
    contains   
      subroutine wr(bmp, fn)
        type(t_rgb), intent(in) :: bmp(:, :)
        character(len = *), intent(in) :: fn
        type(t_bmp_file_header) :: bmp_file_header
        type(t_bmp_info_header) :: bmp_info_header
        integer :: i, nx, ny
        nx = size(bmp, 1)
        ny = size(bmp, 2)
        bmp_file_header%bfSize      = 14 + 40 + 0 + (3 * nx + mod(ny, 4)) * ny
        bmp_file_header%bfOffBits   = 14 + 40
        bmp_info_header%biWidth     = nx
        bmp_info_header%biHeight    = ny
        bmp_info_header%biBitCount  = 24 
        bmp_info_header%biSizeImage = (3 * nx + mod(ny, 4)) * ny
        open(9, file = fn // '.bmp', access = 'stream', status = 'unknown')
        write(9) bmp_file_header
        write(9) bmp_info_header
        write(9) (bmp(:, i), repeat(achar(0), mod(nx, 4)), i = 1, ny)
        close(9)
      end subroutine wr  
    end module m_bmp
[gfort kernel] module objects created successfully: BMP.o
In [20]:
%fcflags: BMP.o
    module m_mandel
      use, intrinsic :: iso_fortran_env
      implicit none
      integer, parameter :: kd = real64
      integer, parameter :: maxiter = 255
    contains
      pure elemental integer function mandel(c)
        complex(kd), intent(in) :: c
        complex(kd) :: z
        z = c
        do mandel = maxiter, 1, -1
          if (abs(z) > 2.0_kd) exit
          z = z**2 + c
        end do 
      end function mandel
    end module m_mandel


    program mandelbrot
      use m_bmp 
      use m_mandel
      implicit none
      complex(kd), allocatable :: c(:, :)
      type(t_rgb), allocatable :: bmp(:, :)
      integer  :: ix, iy, nx, ny
      real(kd) :: xmin, ymin, xmax, ymax
      xmin = 0.15_kd
      xmax = 0.25_kd
      ymin = 0.50_kd
      ymax = 0.60_kd
      nx = 320
      ny = 320
      allocate( c(nx, ny) )
      forall (ix = 1:nx, iy = 1:ny) c(ix, iy) = cmplx(pos(ix, nx, xmin, xmax), pos(iy, ny, ymin, ymax)) 
      bmp = to_rgb( mandel(c) ) 
      call wr(bmp, 'mandel')
      stop 'normal end'
    contains
      pure real(kd) function pos(i, n, rmin, rmax) 
        integer , intent(in) :: i, n
        real(kd), intent(in) :: rmin, rmax
        pos = (rmax - rmin) / (n - 1) * i + rmin
      end function pos
      !
      pure elemental type(t_rgb) function to_rgb(k)
        integer, intent(in) :: k
        to_rgb = t_rgb(10 * k,  k, 5 * k)
      end function to_rgb
    end program mandelbrot
STOP normal end

PIL

In [21]:
%fig:
from PIL import Image
im1 = Image.open('mandel.bmp')
plt.imshow(im1)

Display Images (MIME: image/xxx)

%image: fn1 fn2 ...

In [22]:
%image: mandel.bmp

Python ; limited functionality (exec line by line)

  • %py:
  • print text -> eval(text) ; not pyhon2 nor python3 print
In [23]:
%py:
from sympy import *
x = Symbol('x')
bernoulli_gen = x / (exp(x) - 1)
print bernoulli_gen

bernoulli_ser = series(bernoulli_gen , x, 0, 10).removeO()
btmp = Poly(bernoulli_ser, x).coeffs()
bernoulli = [0, 0, 0, 0, 0, 0]

indx = [0, 1, 2, 4, 6, 8]
for i in range(6): bernoulli[i] = factorial(indx[i]) * btmp[5 - i]
print bernoulli
x/(exp(x) - 1)
[1, -1/2, 1/6, -1/30, 1/42, -1/30]

Bernoulli numbers

$$ B_0=1, B_1=-{1\over2}, B_2={1\over6}, B_4=-{1\over30}, B_6={1\over42}, B_8=-{1\over30} $$