#+TITLE: Programming Language Interstices #+SUBTITLE: Bridging Fortran and Python #+AUTHOR: Rohit Goswami # I need the footnotes to be inlined #+STARTUP: fninline #+EXCLUDE_TAGS: noexport #+BEGIN_SRC emacs-lisp :exports none :eval always (require 'ox-extra) (ox-extras-activate '(ignore-headlines)) (eval unpackaged/org-export-html-with-useful-ids-mode) ;; Stop using citeproc-org by default (setq org-export-before-parsing-hook '(org-ref-acronyms-before-parsing org-ref-glossary-before-parsing org-attach-expand-links)) #+END_SRC #+RESULTS: | org-ref-acronyms-before-parsing | org-ref-glossary-before-parsing | org-attach-expand-links | * Configuration :ignoreheading:ignore: :PROPERTIES: :VISIBILITY: folded :END: # Kanged from https://gitlab.com/oer/oer-reveal/blob/master/org/config.org # Also https://gitlab.com/oer/emacs-reveal-howto/-/blob/master/howto.org ** General Options :ignoreheading:ignore: # No Table of contents, no section numbers #+OPTIONS: toc:nil num:nil # Enable smart quotes #+OPTIONS: ':t ** RevealJS Options :ignoreheading:ignore: # Enable: browser history, fragment IDs in URLs, mouse wheel, links between presentations #+OPTIONS: reveal_history:t reveal_fragmentinurl:t #+OPTIONS: reveal_mousewheel:t reveal_inter_presentation_links:t # Disable separate PDF pages for each fragment. Just use one per slide. #+OPTIONS: reveal_pdfseparatefragments:nil # Display notes on separate page for PDF export. #+REVEAL_EXPORT_NOTES_TO_PDF: separate-page # Transition styles: none/fade/slide/convex/concave/zoom/cube #+REVEAL_TRANS: fade # Set a base theme, then override #+REVEAL_THEME: robot-lung #+REVEAL_EXTRA_CSS: /Users/rohitgoswami/.config/doom/reveal/extras/rlExtras.css #+REVEAL_EXTRA_CSS: /Users/rohitgoswami/.config/doom/reveal/extras/oerFragments.css #+REVEAL_EXTRA_CSS: /Users/rohitgoswami/.config/doom/reveal/extras/noImgBoxes.css #+REVEAL_EXTRA_CSS: /Users/rohitgoswami/.config/doom/reveal/extras/betterCite.css #+REVEAL_EXTRA_CSS: /Users/rohitgoswami/.config/doom/reveal/extras/moreCode.css #+REVEAL_MARGIN: 0.2 #+REVEAL_PREAMBLE:
#+REVEAL_PLUGINS: (notes search zoom) # The following variables are non-standard. # Do not display TOC-progress on title slide. #+REVEAL_TITLE_SLIDE_STATE: no-toc-progress # Do not display TOC-progress on TOC slide. #+REVEAL_TOC_SLIDE_STATE: no-toc-progress # Do not include TOC slide in TOC-progress. #+REVEAL_TOC_SLIDE_CLASS: no-toc-progress # Use different heading for TOC. #+REVEAL_TOC_SLIDE_TITLE: Agenda ** External Resources :ignoreheading:ignore: # Note that doom-emacs sets this variable # https://github.com/hlissner/doom-emacs/blob/develop/modules/lang/org/contrib/present.el #+REVEAL_EXTRA_CSS: /Users/rohitgoswami/.config/doom/reveal/rjs/plugin/accessibility/helper.css #+REVEAL_EXTRA_CSS: /Users/rohitgoswami/.config/doom/reveal/rjs/plugin/toc-progress/toc-progress.css #+REVEAL_EXTRA_CSS: /Users/rohitgoswami/.config/doom/reveal/rjs/dist/theme/toc-style.css #+REVEAL_EXTRA_CSS: /Users/rohitgoswami/.config/doom/reveal/rjs/dist/theme/fonts/source-sans-pro/source-sans-pro.css # Allow to selectively hide links. # #+REVEAL_EXTRA_SCRIPTS: ("/Users/rohitgoswami/.config/doom/reveal/rjs/dist/theme/hidelinks.js") #+REVEAL_EXTRA_SCRIPTS: ("/Users/rohitgoswami/.config/doom/reveal/rjs/dist/theme/hidelinks.js" "/Users/rohitgoswami/.config/doom/reveal/sfeir-school-theme/dist/js/sfeir-theme.js") # The following creates an empty footer, for which the css style defines # a height that agrees with the TOC-progress footer’s height. # In this way, the footer’s height is taken into account by reveal.js’s # size calculations. #+REVEAL_SLIDE_FOOTER:%s
@@ @@latex:\\verb|%s|@@" (org-html-encode-plain-text Modern Documentation across languages) Modern Documentation across languages))
#+MACRO: redcode (eval (format "@@html:%s
@@ @@latex:\\rverb|%s|@@" (org-html-encode-plain-text Modern Documentation across languages) Modern Documentation across languages))
#+MACRO: greencode (eval (format "@@html:%s
@@ @@latex:\\gverb|%s|@@" (org-html-encode-plain-text Modern Documentation across languages) Modern Documentation across languages))
#+MACRO: bluecode (eval (format "@@html:%s
@@ @@latex:\\bverb|%s|@@" (org-html-encode-plain-text Modern Documentation across languages) Modern Documentation across languages))
** References :ignoreheading:ignore:
bibliographystyle:unsrt
#+LATEX_HEADER: \addbibresource{./refs.bib}
** LaTeX Options :ignoreheading:ignore:
# Setup for PDF generation via LaTeX export.
#+LATEX_CLASS_OPTIONS: [a4paper]
#+LATEX_HEADER: \usepackage[backend=biber,style=alphabetic]{biblatex}
#+LATEX_HEADER: \newenvironment{notes}{\par\footnotesize}{\par}
#+LATEX_HEADER: \newenvironment{NOTES}{\par\footnotesize}{\par}
#+LATEX_HEADER: \newenvironment{leftcol}{\begin{minipage}{.49\textwidth}}{\end{minipage}}
#+LATEX_HEADER: \newenvironment{rightcol}{\begin{minipage}{.49\textwidth}}{\end{minipage}}
#+LATEX_HEADER: \newenvironment{leftcol30}{\begin{minipage}{.29\textwidth}}{\end{minipage}}
#+LATEX_HEADER: \newenvironment{leftcol40}{\begin{minipage}{.39\textwidth}}{\end{minipage}}
#+LATEX_HEADER: \newenvironment{leftcol60}{\begin{minipage}{.59\textwidth}}{\end{minipage}}
#+LATEX_HEADER: \newenvironment{leftcol70}{\begin{minipage}{.69\textwidth}}{\end{minipage}}
#+LATEX_HEADER: \newenvironment{rightcol30}{\begin{minipage}{.29\textwidth}}{\end{minipage}}
#+LATEX_HEADER: \newenvironment{rightcol40}{\begin{minipage}{.39\textwidth}}{\end{minipage}}
#+LATEX_HEADER: \newenvironment{rightcol60}{\begin{minipage}{.59\textwidth}}{\end{minipage}}
#+LATEX_HEADER: \newenvironment{rightcol70}{\begin{minipage}{.69\textwidth}}{\end{minipage}}
#+LATEX_HEADER: \usepackage{newunicodechar}
#+LATEX_HEADER: \newunicodechar{≈}{$\approx$}
#+LATEX_HEADER: \newunicodechar{⋮}{\vdots}
#+LATEX_HEADER: \newunicodechar{ }{~}
#+LATEX_HEADER: \usepackage{xcolor}
#+LATEX_HEADER: \definecolor{darkred}{rgb}{0.3, 0.0, 0.0}
#+LATEX_HEADER: \definecolor{darkgreen}{rgb}{0.0, 0.3, 0.1}
#+LATEX_HEADER: \definecolor{darkblue}{rgb}{0.0, 0.1, 0.3}
#+LATEX_HEADER: \definecolor{darkorange}{rgb}{1.0, 0.55, 0.0}
#+LATEX_HEADER: \definecolor{sienna}{rgb}{0.53, 0.18, 0.09}
#+LATEX_HEADER: \hypersetup{colorlinks,linkcolor=darkblue,citecolor=darkblue,urlcolor=darkgreen}
#+LATEX_HEADER: \usepackage{newverbs}
#+LATEX_HEADER: \newverbcommand{\rverb}{\color{darkred}}{}
#+LATEX_HEADER: \newverbcommand{\gverb}{\color{darkgreen}}{}
#+LATEX_HEADER: \newverbcommand{\bverb}{\color{darkblue}}{}
* Start Here :ignoreheading:ignore:
* Brief Introduction
** Hello!
- Find me here: https://rgoswami.me
- Who?
+ *Rohit Goswami* MInstP
- Doctoral Researcher, University of Iceland, Faculty of Physical Sciences
#+begin_leftcol
[[file:logos/physUoI.png]]
#+ATTR_HTML: :width 50% :height 50%
file:logos/rannisLogo.png
#+ATTR_HTML: :width 40% :height 40%
[[file:logos/ccLogo.png]]
#+end_leftcol
#+begin_rightcol
#+ATTR_HTML: :width 40% :height 40%
file:logos/EPH_logo.jpg
#+ATTR_HTML: :width 50% :height 40%
[[file:logos/quansightlabs.jpeg]]
#+end_rightcol
** Logistics
#+ATTR_REVEAL: :frag appear
- All contents are [[https://github.com/HaoZeke/haozeke.github.io][hosted on GitHub]]
+ Slides are in ~presentations/lesHouches21~
#+ATTR_REVEAL: :frag appear
- Questions are welcome *anytime*
* Programming Languages
** Motivation
#+begin_quote
“If a program or package (the words are used interchangeably) is to *have a long life* and to be of *wide application* in its field, it is essential for it to be *easily moved* from one machine to another.
It used to be common to dismiss such movement with the statement, *‘There is no such thing as a machine-independent program.’*
Nonetheless, a great many packages *do now move* from one machine to another”cite:lyonUsingAnsFortran1980
#+end_quote
--> Through the magic of *automated coding* and *standards*
** Language Standards
#+begin_quote
“The standard is the contract between the compiler writer and the application developer.”cite:clermanModernFortranStyle2012
#+end_quote
#+BEGIN_SRC ditaa :file images/hello-program.png :cmdline -r -s 2.5 :cache yes
+------+ +----------+ assembly +-----------+
| Code | --> | Compiler | ----------> | Assembler | ---+
+------+ +----------+ +-----------+ |
relocatable machine code |
+----------------------------<-------------------------+
|
|
| +--------+ executable +--------+ +--------+
+---> | Linker | -----------> | Loader | --> | Memory |
+--------+ +--------+ +--------+
#+END_SRC
#+RESULTS[f1fc83b64fe81184a3a817828e6dd0aec3714f25]:
[[file:images/hello-program.png]]
#+ATTR_REVEAL: :frag appear
- Sadly interpreters are not really Code -> Memory
#+ATTR_REVEAL: :frag appear
+ Python itself is an interpreter which is commonly implemented in C
** Changing Standards
#+begin_leftcol
#+begin_src fortran
character(10) BLAH*8
character*8 :: BLAH_ONE(10)
character(8) :: BLAH_ONE(10)
#+end_src
#+begin_src python
#!/usr/bin/env python
print("Hello World")
print "Hello World"
#+end_src
#+end_leftcol
#+begin_rightcol
#+DOWNLOADED: screenshot @ 2021-09-08 23:12:16
[[file:images/Why_Care_About_New_Standards/2021-09-08_23-12-16_screenshot.png]]
#+end_rightcol
** F77 ∉ F90 always
#+ATTR_HTML: :width 70% :height 70%
[[file:images/Why_Care_About_New_Standards/2021-09-08_23-14-26_screenshot.png]]
#+ATTR_HTML: :width 70% :height 70%
[[file:images/Why_Care_About_New_Standards/2021-09-08_23-14-38_screenshot.png]]
* LFortran
** Introduction
#+DOWNLOADED: screenshot @ 2021-09-09 00:50:56
[[file:images/Introduction/2021-09-09_00-50-56_screenshot.png]]
| *Language* | *Files* | *Lines* | *Code* | *Comments* | *Blanks* |
|------------+---------+---------+---------+------------+----------|
| C | 3 | 1023 | 694 | 191 | 138 |
| C Header | 57 | 14237 | 11416 | 1041 | 1780 |
| CMake | 11 | 430 | 361 | 16 | 53 |
| C++ | 54 | 30745 | 26911 | 1560 | 2274 |
| C++ Header | 1 | 8938 | 8297 | 348 | 293 |
| FORTRAN | 20 | 1738 | 1344 | 174 | 220 |
| Python | 2 | 224 | 191 | 4 | 29 |
|------------+---------+---------+---------+------------+----------|
| *Total* | *148* | *57335* | *49214* | *3334* | *4787* |
** Structure
#+BEGIN_SRC ditaa :file images/hello-lfortran.png :cmdline -r -s 3.5 :cache yes
+------+ +-----------+ +----------------+
| Code | --> | Abstract | --> | (ASR) Abstract | ---+
| | | Syntax | | Semantic | |
| | | Tree (AST)| | Representation | |
+------+ +-----------+ +----------------+ |
|
+--------+ |
+---------| Passes |--------------<------------------+
| +--------+
|
| +--------+
| +-->| Python |
| | +--------+
| |
| +----------+ | +------+
+---> | Backends | --+-->| LLVM |
+----------+ | +------+
|
| +-----+
+-->| C++ |
+-----+
#+END_SRC
#+RESULTS[e593bcddd9ca0862ed60f7f70716ce54612ef366]:
[[file:images/hello-lfortran.png]]
** Features
- Runtime Libraries :: Pure Fortran | Impure
- ASR :: Guarantees compilation --> Wrappers
- Parser :: Currently almost all of F2018; including F77
- LLVM :: Canonical backend, includes compile time evaluated expressions
- Jupyter :: Native execution as a kernel
* Practical Matters
** Numerical Methods 101
$$Ax=b \\
x=A^{-1}b
$$
- Numerical stability issues
- Floating Point Numbers IEEE 754 cite:goldbergWhatEveryComputer1991
- Convergence can be tied to machine epsilon....
+ Or a *magic number* $(1e-5,1e-16)$
+ Or *Chemical/Spectroscopic* accuracy
#+begin_src fortran :exports both
program main
print*, tiny(0.d0)
end program
#+end_src
#+RESULTS:
: 2.2250738585072014e-308
** Yawn
:PROPERTIES:
:reveal_background: #99e6ff
:END:
#+begin_quote
“I don't care because *I use libraries.* ”
#+end_quote
#+begin_quote
“I don't care because *I write my own algorithms*.”
#+end_quote
** Implementing Sine
#+DOWNLOADED: screenshot @ 2021-09-09 11:10:27
#+ATTR_HTML: :width 80% :height 40%
[[file:images/Implementing_Sine/2021-09-09_11-10-27_screenshot.png]]
#+begin_src fortran
real(dp) function dsin(x) result(r)
real(dp), parameter :: pi = 3.1415926535897932384626433832795_dp
real(dp), intent(in) :: x ! Assume folded to [-2𝜋,2𝜋]
real(dp) :: y, xnew, twoPi, invTwoPi
if (abs(x) < pi/2) then
r = kernel_dsin(x)
else ! fold to pi/2 from https://realtimecollisiondetection.net/blog/?p=9
y = modulo(xnew, 2*pi)
y = min(y, pi - y)
y = max(y, -pi - y)
y = min(y, pi - y)
r = kernel_dsin(y)
end if
end function
#+end_src
*** Folding I
- Symmetric about $𝜋/2$ --> $x-𝜋/2$
[[file:images/Implementing_Sine/2021-09-09_12-43-34_screenshot.png]]
- $x=-abs(x)$
#+DOWNLOADED: screenshot @ 2021-09-09 12:44:54
[[file:images/Implementing_Sine/2021-09-09_12-44-54_screenshot.png]]
*** Folding II
- $x = -abs(x – 𝜋/2) + 𝜋/2$ -> $x=min(x,𝜋-x)$
#+DOWNLOADED: screenshot @ 2021-09-09 12:46:39
[[file:images/Implementing_Sine/2021-09-09_12-46-39_screenshot.png]]
*** Folding III
- Repeat rightward to get $x=max(x,-𝜋-x)$
#+DOWNLOADED: screenshot @ 2021-09-09 12:47:42
[[file:images/Implementing_Sine/2021-09-09_12-47-42_screenshot.png]]
- Complete with $x=min(x,𝜋-x)$
#+DOWNLOADED: screenshot @ 2021-09-09 12:48:04
[[file:images/Implementing_Sine/2021-09-09_12-48-04_screenshot.png]]
** Kernel Sine
#+begin_leftcol
#+DOWNLOADED: screenshot @ 2021-09-09 11:22:47
#+ATTR_HTML: :width 50% :height 50%
[[file:images/Kernel_Sine/2021-09-09_11-22-47_screenshot.png]]
cite:ChevillardJoldesLauter2010
#+begin_src bash
nix-shell -p sollya rlwrap
rlwrap -A sollya
prec=120;
f=sin(x);
d=[-pi/2;pi/2];
# Use min/max poly
p=remez(f,16,d,
default,1e-16,1e-30);
p; # Returns coefficients
# Zero out even terms
#+end_src
#+end_leftcol
#+begin_rightcol
#+begin_src fortran
! Accurate on [-pi/2,pi/2] to about 1e-16
elemental real(dp) function kernel_dsin2(x) result(res)
real(dp), intent(in) :: x
real(dp), parameter :: S1 = 0.9999999999999990771_dp
real(dp), parameter :: S2 = -0.16666666666664811048_dp
real(dp), parameter :: S3 = 8.333333333226519387e-3_dp
real(dp), parameter :: S4 = -1.9841269813888534497e-4_dp
real(dp), parameter :: S5 = 2.7557315514280769795e-6_dp
real(dp), parameter :: S6 = -2.5051823583393710429e-8_dp
real(dp), parameter :: S7 = 1.6046585911173017112e-10_dp
real(dp), parameter :: S8 = -7.3572396558796051923e-13_dp
real(dp) :: z
z = x*x
res = x * (S1+z*(S2+z*(S3+z*(S4+z*(S5+z*
(S6+z*(S7+z*S8)))))))
end function
#+end_src
#+end_rightcol
** Auxiliary Functions
#+begin_src fortran
elemental real(dp) function dmodulo(x, y) result(r)
real(dp), intent(in) :: x, y
r = x-floor(x/y)*y
end function
elemental integer function dfloor(x) result(r)
real(dp), intent(in) :: x
if (x >= 0) then
r = x
else
r = x-1
end if
end function
elemental real(dp) function dabs(x) result(r)
real(dp), intent(in) :: x
if (x >= 0) then
r = x
else
r = -x
end if
end function
#+end_src
** Accuracy
With an evenly spaced grid of ~0.01~:
#+ATTR_REVEAL: :frag appear
- $1E-16$ over kernel range $-\pi/2, \pi/2$
#+ATTR_REVEAL: :frag appear
- $1E-16$ over $[-20,20]$
#+ATTR_REVEAL: :frag appear
- $1E-12$ from $1E5$
#+ATTR_REVEAL: :frag appear
- $1E-7$ from $1E10$
*** Post Mortem
- The problem turned out to be ~modulo~
+ IEEE defines a remainder for trignometric functions ~ieee754_rem_pio2~
#+begin_quote
Do not *naively assume* "handcrafted" algorithms are IEEE 754 accurate
#+end_quote
* Fortran, C, Python
- F2003 :: Introduced the ~ISO_C_BINDING~
- Described in great detail on [[https://www.fortran90.org/src/best-practices.html#python-interface][fortran90.org]]:
#+begin_src ditaa :file images/hello-stdfcpy.png :cmdline -r -s 1.5 :cache yes
+--------------+ iso_c_binding +------------+
| Fortran Code | --------------> | Call in C |
+--------------+ compiler type +------------+
checking |
v
+--------+ +--------+
| Python | | Cython |
| Code |<--------------| cffi |
+--------+ +--------+
#+end_src
#+RESULTS[c649e6871e7d2b61cff7d210af3db78f85b0fcc5]:
[[file:images/hello-stdfcpy.png]]
* F2PY
** History
- Developed by Pearu Peterson cite:petersonF2PYToolConnecting2009
+ July 9, 1999 :: ~f2py.py~ --> Fortran to Python Interface Generator (FPIG)
+ January 22, 2000 :: ~f2py2e~ --> Fortran to Python Interface Generator, 2nd edition.
+ July 19, 2007 :: ~numpy.f2py~ --> f2py2e moved to NumPy project. This is current stable code of f2py.
- Used extensively for F77
+ [[https://numpy.org/][NumPy]] cite:waltNumPyArrayStructure2011, [[https://scipy.org/][SciPy]] cite:virtanenSciPyFundamentalAlgorithms2020
+ [[https://msspec.cnrs.fr/][MsSpec]] cite:sebilleauMsSpec1MultipleScattering2011 :)
** Design
#+BEGIN_SRC ditaa :file images/hello-f2py.png :cmdline -r -s 2.5 :cache yes
+------+ +--------------+ uses C syntax in pyf
| Code | --> | crackfortran | ---------->----------+
+------+ +--------------+ |
match rules, generate wrappers, build library |
+----------------------------<---------------------+
| +---------+ callbacks +--------+
| | CPython | <---------| Python |
+---> | Library | --------->| Code |
+---------+ +--------+
#+END_SRC
#+RESULTS[7f80fadd5cfd09443054d3693a2d688d9a0b639d]:
[[file:images/hello-f2py.png]]
- A *best effort* wrapper
+ Specifications via ~.pyf~ or inline comments
+ *Not* a compiler
- Can rewrite code :)
** Roadmap
- Moves about as fast as spectroscopy codes :)
+ Also picking up speed in 2021
- Implementing newer standards (90, 95, 2003, 2008, 2018, 2020Y)
+ Automating guarantees
- Documentation
** Relevance
:PROPERTIES:
:reveal_background: #99e6ff
:END:
#+begin_quote
Writing *efficient* wrappers without being a language lawyer
#+end_quote
* Bonus: Reading Code
** Reading Code I
#+begin_leftcol
#+ATTR_REVEAL: :frag appear
#+begin_src asm
main:
push rbp
mov rbp, rsp
mov DWORD PTR [rbp-4], 3
mov eax, 0
pop rbp
ret
__static_initialization_
and_destruction_0(int, int):
push rbp
mov rbp, rsp
sub rsp, 16
mov DWORD PTR [rbp-4], edi
mov DWORD PTR [rbp-8], esi
cmp DWORD PTR [rbp-4], 1
jne .L5
cmp DWORD PTR [rbp-8], 65535
jne .L5
mov edi, OFFSET FLAT:_ZStL8
__ioinit
#+end_src
#+end_leftcol
#+begin_rightcol
#+ATTR_REVEAL: :frag appear
#+begin_src asm
call std::ios_base::Init::Init()
[complete object constructor]
mov edx, OFFSET FLAT:__dso_handle
mov esi, OFFSET FLAT:_ZStL8__ioinit
mov edi, OFFSET FLAT:_ZNSt8ios_base4InitD1Ev
call __cxa_atexit
.L5:
nop
leave
ret
_GLOBAL__sub_I_main:
push rbp
mov rbp, rsp
mov esi, 65535
mov edi, 1
call __static_initialization_
and_destruction_0(int, int)
pop rbp
ret
#+end_src
#+end_rightcol
#+ATTR_REVEAL: :frag appear
- But who *writes _assembly_* anyway?
** Reading Code II
#+begin_leftcol
#+ATTR_REVEAL: :frag appear
#+begin_src cpp
int main ()
{
int D.48918;
{
int a;
a = 3;
D.48918 = 0;
return D.48918;
}
D.48918 = 0;
return D.48918;
}
void _GLOBAL__sub_I_main.cpp ()
{
__static_initialization_
and_destruction_0 (1, 65535);
}
#+end_src
#+end_leftcol
#+begin_rightcol
#+ATTR_REVEAL: :frag appear
#+begin_src cpp
void __static_initialization_
and_destruction_0 (int __initialize_p,
int __priority)
{
if (__initialize_p == 1) goto