# In-situ tuning of optomechanical crystals with nano-oxidation

Utku Hatipoglu  Institute for Quantum Information and Matter, California
Institute of Technology, Pasadena, CA 91125, USA  Kavli Nanoscience Institute
and Thomas J. Watson, Sr., Laboratory of Applied Physics, California Institute
of Technology, Pasadena, CA 91125, USA  These authors contributed equally to
this work.  Sameer Sonar  Institute for Quantum Information and Matter,
California Institute of Technology, Pasadena, CA 91125, USA  Kavli Nanoscience
Institute and Thomas J. Watson, Sr., Laboratory of Applied Physics, California
Institute of Technology, Pasadena, CA 91125, USA  These authors contributed
equally to this work.  David P. Lake  Institute for Quantum Information and
Matter, California Institute of Technology, Pasadena, CA 91125, USA  Kavli
Nanoscience Institute and Thomas J. Watson, Sr., Laboratory of Applied
Physics, California Institute of Technology, Pasadena, CA 91125, USA  Srujan
Meesala  Institute for Quantum Information and Matter, California Institute of
Technology, Pasadena, CA 91125, USA  Kavli Nanoscience Institute and Thomas J.
Watson, Sr., Laboratory of Applied Physics, California Institute of
Technology, Pasadena, CA 91125, USA  Oskar Painter  Institute for Quantum
Information and Matter, California Institute of Technology, Pasadena, CA
91125, USA  Kavli Nanoscience Institute and Thomas J. Watson, Sr., Laboratory
of Applied Physics, California Institute of Technology, Pasadena, CA 91125,
USA  AWS Center for Quantum Computing, Pasadena, CA 91125, USA
<EMAIL_ADDRESS>
###### Abstract

Optomechanical crystals are a promising device platform for quantum
transduction and sensing. Precise targeting of the optical and acoustic
resonance frequencies of these devices is crucial for future advances on these
fronts. However, fabrication disorder in these wavelength-scale nanoscale
devices typically leads to inhomogeneous resonance frequencies. Here we
achieve in-situ, selective frequency tuning of optical and acoustic resonances
in silicon optomechanical crystals via electric field-induced nano-oxidation
using an atomic-force microscope. Our method can achieve a tuning range $>2$
nm ($0.13\%$) for the optical resonance wavelength in the telecom C-band, and
$>60$ MHz ($1.2\%$) for the acoustic resonance frequency at $5$ GHz. The
tuning resolution of $1.1$ pm for the optical wavelength, and $150$ kHz for
the acoustic frequency allows us to spectrally align multiple optomechanical
crystal resonators using optimal oxidation patterns determined via an inverse
design protocol. Our results establish a method for precise post-fabrication
tuning of optomechancical crystals. This technique can enable coupled
optomechanical resonator arrays, scalable resonant optomechanical circuits,
and frequency matching of microwave-optical quantum transducers.

††journal: opticajournal

##  1 Introduction

Optomechanical crystals (OMCs) provide a coherent interface between optical
photons and acoustic phonons [1]. This capability is now being utilized in
microwave-optical quantum transducers [2, 3, 4, 5, 6] towards connecting
gigahertz frequency superconducting quantum processors [7] via low-loss
optical communication channels. In parallel efforts, OMCs have enabled
coherent control and routing of phonons in chip-scale optomechanical circuits
[8, 9, 10] with demonstrations of nonreciprocal optical transmission [11] and
topologically protected phonon transport [12]. However, scaling to multi-node
quantum networks and more complex optomechanical circuits is challenging due
to variations in optical and acoustic resonance frequencies across devices
caused by fabrication imperfections. In particular, since critical feature
sizes in OMCs are affected by the precision limits of electron beam
lithography and reactive ion etching, the spread in resonance frequencies
across devices can be more than a hundred times the resonance linewidths,
thereby hindering frequency alignment.

Figure 1: Nano-oxidation setup schematic and OMC cavity mode profiles. a,
Simplified schematic of the AFM nano-oxidation setup. The AFM tip is used to
perform nano-oxidation while measuring the optical and acoustic resonances of
the device in real-time with a dimpled optical fiber. The AFM tip is operated
in tapping mode over a grounded silicon sample and biased with a square wave
voltage, allowing for electrochemical formation of an oxide layer on the
silicon surface. Inset illustrates the raster scan used to generate the oxide
layer. b, Electric field intensity profile of the optical mode. The optics-
focused oxidation region is outlined in green. c, Displacement magnitude
profile of the acoustic breathing mode. The acoustics-focused oxidation
regions are outlined in green. Figure 2: Coarse frequency tuning of OMC
cavity. a, Oxidation patterns used for coarse tuning of optical (left) and
acoustic (right) resonances. For the optics-focused pattern, the length is
fixed at $l_{0}=170$ nm and the width, $w$ is varied. For the acoustics-
focused pattern, the width is fixed at $w_{0}=120$ nm and the length, $l$ is
varied. b, Tuning of optical and acoustic resonances as a function of the size
of the oxidation pattern. Optical wavelength shift $\delta\lambda$ with
respect to $w$ (top left) and $l$ (top right). Acoustic frequency shift
$\delta\Omega/2\pi$ with respect to $w$ (bottom left) and $l$ (bottom right).
Solid (dashed) lines represent data taken for aggressive (mild) oxidation.
Dotted lines show simulated values from finite element method simulations.
Error bars represent the standard deviation across five OMC cavities measured
for each dataset.

Post-fabrication tuning of optical resonance frequencies in chip-scale
microcavities has been achieved using a variety of techniques including laser-
assisted thermal oxidation [13], atomic force microscope (AFM) nano-oxidation
[14, 15], strain tuning [16, 17], thermo-optic tuning [18, 19], and gas
condensation [20, 21]. However, since OMCs co-localize optical and acoustic
resonances in a wavelength-scale volume, selective tuning of both optical and
acoustic resonances without compromising the quality factors of either is a
significantly more complex endeavor and an outstanding technical challenge.
Here we use AFM nano-oxidation tuning to demonstrate such control over the
resonance frequencies of OMCs. In our approach, field-induced oxidation of the
silicon device surface with high spatial resolution allows us to tune the
acoustic resonance by modifying the local mass distribution and elasticity,
and the optical resonance by modifying the local refractive index. By using an
inverse design protocol to guide the nano-oxidation sequence in real-time, we
achieve simultaneous alignment of the optical and acoustic resonance
frequencies of multiple OMC cavities.

##  2 AFM nano-oxidation setup

A simplified schematic of the experimental setup is shown in Fig. 1a. We
perform nano-oxidation using a conductive chromium/platinum coated silicon AFM
tip with a radius <25 nm [22], and track the optical and acoustic resonances
in real-time by performing optomechanical spectroscopy via a dimpled optical
fiber coupled to the device. The AFM is operated in tapping mode while a
voltage bias is applied to the conductive tip, and the silicon-on-insulator
(SOI) chip with a silicon device layer resistivity of 5 k$\Omega$-cm is
grounded. When the voltage-biased AFM tip is brought close to the silicon
surface, a strong electric field triggers an electrochemical reaction between
ions in the native water meniscus and the silicon surface, resulting in local
oxidation of silicon. Importantly, this reaction can proceed even when the
silicon surface is covered with native oxide, because the strong local
electric field allows oxyanions (OH- , O-) to diffuse through the native oxide
[23]. In our experiments, we maintain the AFM setup enclosure at ambient
conditions with a relative humidity of 40$\pm$5$\%$ and temperature of
23$\pm$1∘C, and apply a square wave voltage oscillating between $\pm 10$ V at
a frequency of 20 Hz to the AFM tip.

Figure 3: Real-time, in-situ monitoring of the oxidation process. a, Real-time
spectra showing optical resonance tuning during sequential application of the
oxidation pattern. As shown in the inset, the pattern is applied starting at
point A, one pixel at a time, up to point C. AFM tip locations during the
experiment are shown at the top of the plot with letters A, B and C which
correspond to the points labeled on the inset image. Discontinuities along the
horizontal axis are due to proximity effects between the AFM tip and the OMC
cavity during the oxidation sequence. b, Real-time data showing acoustic
frequency shift during pixel-wise oxidation. Oxide pixels are applied on the
surface starting from point A up to point B as shown in the inset.

The OMC cavities in this work were fabricated on an SOI chip with a $220$ nm
thick silicon device layer. The cavities were patterned via electron beam
lithography followed by reactive ion etching, and finally suspended by
removing the buried oxide layer with an anhydrous vapor hydrofluoric acid
etch. The devices were designed to support a fundamental TE-like optical
resonance at a wavelength of $\sim$1550 nm and a breathing acoustic resonance
at a frequency of $\sim$5.1 GHz (Fig. 1b,c) [24]. An on-chip adiabatic
waveguide coupler [25] allows for optical coupling to the OMC cavity using a
dimpled optical fiber. The optical spectrum of the OMC cavity was probed in
reflection mode, and the resonance frequencies were recorded before and after
each step of the nano-oxidation sequence. To measure the acoustic spectrum of
the device, we routed the optical signal reflected from the device to a high-
speed photodetector and measured the power spectral density (PSD) of the
photodetector electrical output on a spectrum analyzer. The result is
proportional to the PSD of thermal displacement fluctuations of the acoustic
mode which are transduced onto the optical signal via the optomechanical
interaction in the device. In these measurements, the laser was blue-detuned
with respect to the optical resonance frequency by a detuning close to the
acoustic resonance frequency, and operated at low power to minimize
optomechanical back-action on the acoustic mode [26].

##  3 Nano-oxidation characterization

From the simulated mode profiles shown in Fig. 1b,c, we identified two
strategic oxidation regions for selective tuning of the optical and acoustic
modes. Nano-oxidation in the region shown with the green rectangle in Fig. 1b
is expected to induce a relatively large change in the optical resonance due
to a high concentration of electric field energy. Here the impact on the
acoustic resonance is expected to be minimal due to two reasons. First, the
displacement amplitude is small in this region. Second, even though the stress
amplitude is large in this region, the Young’s modulus of the oxide layer when
weighted by its thickness is similar to that of the original silicon. In
contrast, nano-oxdiation on the regions shown with the green rectangles in
Fig. 1c leads to a large acoustic resonance shift and a small optical
wavelength shift due to the high concentration of motional mass and low
concentration of electric field energy. In section 5.4 of the supplementary
information, we show the results of finite element method (FEM) simulations
where we investigated the optical and acoustic frequency shifts due to
creation of an oxide pixel at an arbitrary location on the OMC surface. This
allowed us to identify candidate regions for fine tuning of optical and
acoustic resonance frequencies with maximal selectivity. For coarse frequency
tuning, we opted to use larger, rectangular oxide patterns instead of oxide
pixels to maximize the tuning range at the expense of selectivity.

Prior to experiments on OMC devices, we characterized the effect of the AFM
tapping amplitude and scanning velocity on nano-oxidation using a bare silicon
chip. We observed that lower tapping amplitude and slower scanning velocities
result in wider and thicker oxide lines. We measured the three-dimensional
profiles of field-induced oxide lines and pixels using AFM measurements as
described in section 5.2 of the supplementary information. Based on the amount
of frequency tuning required, we operated the AFM in two distinct modes which
we refer to as ‘mild’ and ‘aggressive’ tuning modes. In the mild tuning mode,
the tapping amplitude and scan speed were set to 15 nm and 100 nm/s,
respectively, resulting in oxide thickness of approximately 1.2 nm. In the
aggressive tuning mode, the tapping amplitude and scan speed were set to 5 nm
and 50 nm/s, respectively, leading to oxide thickness of approximately 2.5 nm.
In the mild tuning mode, nano-oxidation was found to generate single pixels
with a lateral size of approximately 25 nm. In order to obtain thicker and
more uniform oxide patterns with a raster scan, we used a line spacing of 13
nm with significant overlap between neighboring lines as illustrated in the
inset of Fig. 1a. This procedure allowed us to generate oxide layers with a
thickness of 1.6 nm and 3.2 nm in the mild and aggressive tuning modes,
respectively.

After establishing these nano-oxidation parameters, we performed optical and
acoustic resonance frequency tuning experiments on a single OMC cavity. We
specifically studied the impact of two oxidation patterns: (i) an optics-
focused pattern (Fig. 2a left) composed of four rectangles located between the
five central ellipses in the beam; here the lengths of the rectangles were
kept constant at $l_{0}$ = 170 nm, while their width, $w$ was varied, and (ii)
an acoustics-focused pattern (Fig. 2a, right) composed of two rectangles
situated at the edge of the OMC cavity; here the widths of the rectangles were
kept fixed at $w_{0}=$ 120 nm, while their length, $l$ was varied. As shown in
Fig. 2b, the experimentally measured tuning of the resonance frequencies
(solid and dashed lines) agrees with simulation results (dotted lines). For
the optics-focused pattern, increasing $w$ leads to a large reduction in the
optical wavelength and a smaller reduction in acoustic frequency. For the
acoustics-focused pattern, increasing $l$ leads to a decrease in acoustic
frequency due to increasing motional mass ($\Omega_{m}\propto\sqrt{1/m}$) and
a decrease in the optical wavelength due to a decrease in the effective
dielectric constant. Aggressive oxidation at $w=$ 320 nm resulted in more than
1.2 nm blue shift of the optical wavelength while aggressive oxidation at $l$
= 1.6 $\mu$m resulted in more than 50 MHz reduction in the acoustic frequency.
In Fig. 2b, each data point represents the average value and the standard
deviation obtained from five devices. For each device, the patterns were
applied cumulatively, starting from the center and adding more oxide to
increase either $w$ or $l$ at each step.

Figure 4: Simultaneous optical and acoustic resonance frequency alignment of
three OMC cavities. Initial a, optical and b, acoustic resonance spectra. c,
Joint patterns used for coarse alignment of the three OMC cavities. The
patterns are shown with a color code matched to the OMC spectra, and are
applied in the order shown from left to right. Resonance frequencies were
recorded after every step of the nano-oxidation protocol. d, Optical and e,
acoustic spectra measured after completion of the coarse alignment of
sequence. f, Optical and g, acoustic spectra after completion of the real-time
fine-tuning sequence. Insets in panels d and f show a zoomed-in view of the
optical resonances. The satellite acoustic resonance in magenta in panels e
and g is due to a neighboring OMC cavity within the evanescent field of the
optical coupling waveguide.

##  4 Simultaneous optical and acoustic frequency tuning

To achieve simultaneous tuning of optical and acoustic resonances, we define a
joint pattern consisting of a union of the previously introduced optics-
focused and acoustics-focused patterns. This is parameterized by geometric
parameters $w$ and $l$, and binary variables $M_{\text{O}}$ ($M_{\text{A}}$)
which denote mild or aggressive oxidation in the optical (acoustic) components
of the joint pattern. Using the experimental data shown in Fig. 2, we
interpolate $\delta\lambda$ and $\delta\Omega_{\text{m}}$ as a function of the
parameters $\\{w,l,M_{\text{O}},M_{\text{A}}\\}$. To achieve the desired
frequency shifts, we developed an inverse design algorithm to find an optimal
combination of the acoustic and optics focused patterns. The algorithm takes
the desired frequency and wavelength shift as the input and generates the
required joint oxidation pattern by searching for the most suitable
combination of $l$ and $w$. The algorithm starts the search from milder
oxidation modes to prevent a reduction in optical Q-factor as explained in
supplementary section 5.5.

While the inverse design algorithm is sufficient for coarse tuning of the OMC
wavelength and frequency, fine frequency alignment to within the linewidth of
the optical and acoustic resonances requires pixel-by-pixel nano-oxidation
with real-time feedback. This method allows for $\sim$1 pm precision for
optical wavelength tuning and $\sim$150 kHz precision for acoustic resonance
tuning. To quantify the sensitivity of the optical resonance wavelength to
single oxide pixels, we applied a pixel-by-pixel linear oxide pattern
perpendicular to the OMC long axis with a pitch of 13 nm shown in the inset of
Fig 3a. The resonance spectra recorded in Fig. 3a for each incremental pixel
in the oxidation sequence reveal a decrease in the per-pixel wavelength shift
towards the edges of the OMC, with an observed shift of 1.1 pm/px at the
extreme points A and C. This spatial dependence is in agreement with
predictions from FEM simulations, as discussed in supplementary section 5.4.
We perform a similar sensitivity analysis for the acoustic resonance via a
pixel-by-pixel linear oxide pattern applied with a pitch of 40 nm along the
edge of the OMC, shown in the inset of Fig 3b. The acoustic resonance spectrum
is recorded after the addition of every oxide pixel, and the results are shown
in Fig. 3b. We measure a maximum frequency shift of -1 MHz/px at point A which
is in line with the center of the OMC cavity. As we move away from the center,
the sensitivity eventually decreases to a value of -150 kHz/px as measured at
the extreme point B.

After characterizing coarse and fine tuning of optical and acoustic resonance
frequencies of a single OMC cavity, we used these techniques to perform
frequency alignment of both the optical and acoustic resonances of three OMC
cavities. Due to constraints on the maximum achievable tuning range, we pre-
select cavities with initial optical wavelength spread <1 nm and acoustic
frequency spread <20 MHz, shown on Fig. 4a, b. We begin with coarse tuning
steps using oxidation patterns generated by the inverse design algorithm as
shown in Fig. 4c. This aligns the optical resonance wavelengths across all
cavities to within a tolerance of <200 pm and acoustic resonance frequencies
to withina tolerance of 5 MHz. The results of the coarse tuning are shown on
Fig. 4d,e. We then proceed to real-time, pixel-by-pixel oxidation with
feedback to achieve fine tuning of the resonance frequencies. After four
rounds of fine-tuning, the optical and acoustic resonance frequencies were
successfully tuned to within 2 pm and 500 kHz respectively (see Fig 4f, g). A
step-by-step description of the frequency alignment procedure is described in
detail in supplementary section 5.6.

##  5 Outlook

We have shown simultaneous, real-time tuning of optical and acoustic
resonances in OMC cavities via AFM nano-oxidation, demonstrating optical
wavelength and acoustic frequency shifts >2 nm (0.13%) and >60 MHz (1.2%),
respectively. Such in-situ tuning methods could enable the realization of
coupled optomechanical resonator arrays to study topological phases of photons
and phonons [27, 28, 29, 12]. In the context of OMC-based microwave-optical
quantum transducers [2, 3, 4, 5, 6], the nano-oxidation tuning demonstrated
here will enable post-fabrication frequency alignment of transducers in
multiple nodes of a quantum network. Moreover, these techniques are applicable
to other material systems such as silicon nitride [30], where tuning of
optomechanical devices may be required [31]. Finally, AFM nano-oxidation could
be a valuable technique for fundamental studies of two-level systems (TLS) in
amorphous materials [32] by enabling spatially targeted creation of TLS in
acoustic [1] and superconducting [33] devices.

## References

  * [1] G. S. MacCabe, H. Ren, J. Luo, J. D. Cohen, H. Zhou, A. Sipahigil, M. Mirhosseini, and O. Painter, “Nano-acoustic resonator with ultralong phonon lifetime,” Science 370, 840–843 (2020).
  * [2] M. Mirhosseini, A. Sipahigil, M. Kalaee, and O. Painter, “Superconducting qubit to optical photon transduction,” Nature 588, 599–603 (2020).
  * [3] S. Meesala, S. Wood, D. Lake, P. Chiappina, C. Zhong, A. D. Beyer, M. D. Shaw, L. Jiang, and O. Painter, “Non-classical microwave-optical photon pair generation with a chip-scale transducer,” arXiv preprint arXiv:2303.17684 (2023).
  * [4] W. Jiang, F. M. Mayor, S. Malik, R. Van Laer, T. P. McKenna, R. N. Patel, J. D. Witmer, and A. H. Safavi-Naeini, “Optically heralded microwave photon addition,” Nature Physics 19, 1423–1428 (2023).
  * [5] M. J. Weaver, P. Duivestein, A. C. Bernasconi, S. Scharmer, M. Lemang, T. C. v. Thiel, F. Hijazi, B. Hensen, S. Gröblacher, and R. Stockill, “An integrated microwave-to-optics interface for scalable quantum computing,” Nature Nanotechnology pp. 1–7 (2023).
  * [6] H. Zhao, A. Bozkurt, and M. Mirhosseini, “Electro-optic transduction in silicon via gigahertz-frequency nanomechanics,” Optica 10, 790–796 (2023).
  * [7] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver, “Superconducting qubits: Current state of play,” Annual Review of Condensed Matter Physics 11, 369–395 (2020).
  * [8] K. Fang, M. H. Matheny, X. Luan, and O. Painter, “Optical transduction and routing of microwave phonons in cavity-optomechanical circuits,” Nature Photonics 10, 489–496 (2016).
  * [9] R. N. Patel, Z. Wang, W. Jiang, C. J. Sarabalis, J. T. Hill, and A. H. Safavi-Naeini, “Single-mode phononic wire,” Phys. Rev. Lett. 121, 040501 (2018).
  * [10] A. Zivari, N. Fiaschi, R. Burgwal, E. Verhagen, R. Stockill, and S. Gröblacher, “On-chip distribution of quantum information using traveling phonons,” Science Advances 8, eadd2811 (2022).
  * [11] K. Fang, J. Luo, A. Metelmann, M. H. Matheny, F. Marquardt, A. A. Clerk, and O. Painter, “Generalized non-reciprocity in an optomechanical circuit via synthetic magnetism and reservoir engineering,” Nature Physics 13, 465–471 (2017).
  * [12] H. Ren, T. Shah, H. Pfeifer, C. Brendel, V. Peano, F. Marquardt, and O. Painter, “Topological phonon transport in an optomechanical system,” Nature Communications 13, 3476 (2022).
  * [13] C. L. Panuski, I. Christen, M. Minkov, C. J. Brabec, S. Trajtenberg-Mills, A. D. Griffiths, J. J. McKendry, G. L. Leake, D. J. Coleman, C. Tran _et al._ , “A full degree-of-freedom spatiotemporal light modulator,” Nature Photonics 16, 834–842 (2022).
  * [14] K. Hennessy, C. Högerle, E. Hu, A. Badolato, and A. Imamoğlu, “Tuning photonic nanocavities by atomic force microscope nano-oxidation,” Applied physics letters 89, 041118 (2006).
  * [15] A. Yokoo, T. Tanabe, E. Kuramochi, and M. Notomi, “Ultrahigh-q nanocavities written with a nanoprobe,” Nano letters 11, 3634–3642 (2011).
  * [16] C. W. Wong, P. T. Rakich, S. G. Johnson, M. Qi, H. I. Smith, E. P. Ippen, L. C. Kimerling, Y. Jeon, G. Barbastathis, and S.-G. Kim, “Strain-tunable silicon photonic band gap microcavities in optical waveguides,” Applied physics letters 84, 1242–1244 (2004).
  * [17] I. Luxmoore, E. D. Ahmadi, B. Luxmoore, N. Wasley, A. Tartakovskii, M. Hugues, M. Skolnick, and A. Fox, “Restoring mode degeneracy in h1 photonic crystal cavities by uniaxial strain tuning,” Applied Physics Letters 100, 121116 (2012).
  * [18] M. S. Nawrocka, T. Liu, X. Wang, and R. R. Panepucci, “Tunable silicon microring resonator with wide free spectral range,” Applied physics letters 89, 071110 (2006).
  * [19] A. Faraon, D. Englund, I. Fushman, J. Vučković, N. Stoltz, and P. Petroff, “Local quantum dot tuning on photonic crystal chips,” Applied Physics Letters 90, 213110 (2007).
  * [20] S. Mosor, J. Hendrickson, B. Richards, J. Sweet, G. Khitrova, H. Gibbs, T. Yoshie, A. Scherer, O. Shchekin, and D. Deppe, “Scanning a photonic crystal slab nanocavity by condensation of xenon,” Applied Physics Letters 87, 141105 (2005).
  * [21] S. Strauf, M. T. Rakher, I. Carmeli, K. Hennessy, C. Meier, A. Badolato, M. J. DeDood, P. M. Petroff, E. L. Hu, E. Gwinn _et al._ , “Frequency control of photonic crystal membrane resonators by monolayer deposition,” Applied Physics Letters 88, 043116 (2006).
  * [22] “Atomic force microscopy: Park systems,” .
  * [23] R. Garcia, R. V. Martinez, and J. Martinez, “Nano-chemistry and scanning probe nanolithographies,” Chemical Society Reviews 35, 29–38 (2006).
  * [24] J. Chan, A. H. Safavi-Naeini, J. T. Hill, S. Meenehan, and O. Painter, “Optimized optomechanical crystal cavity with acoustic radiation shield,” Applied Physics Letters 101, 081115 (2012).
  * [25] S. Gröblacher, J. T. Hill, A. H. Safavi-Naeini, J. Chan, and O. Painter, “Highly efficient coupling from an optical fiber to a nanoscale silicon optomechanical cavity,” Applied Physics Letters 103, 181104 (2013).
  * [26] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, “Cavity optomechanics,” Reviews of Modern Physics 86, 1391 (2014).
  * [27] V. Peano, C. Brendel, M. Schmidt, and F. Marquardt, “Topological phases of sound and light,” Physical Review X 5, 031011 (2015).
  * [28] V. Peano, M. Houde, C. Brendel, F. Marquardt, and A. A. Clerk, “Topological phase transitions and chiral inelastic transport induced by the squeezing of light,” Nature communications 7, 1–8 (2016).
  * [29] C. Brendel, V. Peano, O. J. Painter, and F. Marquardt, “Pseudomagnetic fields for sound at the nanoscale,” Proceedings of the National Academy of Sciences 114, E3390–E3395 (2017).
  * [30] F.-S. Chien, Y. Chou, T. Chen, W.-F. Hsieh, T.-S. Chao, and S. Gwo, “Nano-oxidation of silicon nitride films with an atomic force microscope: Chemical mapping, kinetics, and applications,” Journal of Applied Physics 89, 2465–2472 (2001).
  * [31] R. D. Delaney, M. D. Urmey, S. Mittal, B. M. Brubaker, J. M. Kindem, P. S. Burns, C. A. Regal, and K. W. Lehnert, “Superconducting-qubit readout via low-backaction electro-optic transduction,” Nature 606, 489–493 (2022).
  * [32] J. Lisenfeld, A. Bilmes, S. Matityahu, S. Zanker, M. Marthaler, M. Schechter, G. Schön, A. Shnirman, G. Weiss, and A. V. Ustinov, “Decoherence spectroscopy with individual two-level tunneling defects,” Scientific Reports 6, 23786 (2016).
  * [33] M. Chen, J. C. Owens, H. Putterman, M. Schäfer, and O. Painter, “Phonon engineering of atomic-scale defects in superconducting quantum circuits,” arXiv preprint arXiv:2310.03929 (2023).

Funding
This work was supported by the ARO/LPS Cross Quantum Technology Systems
program (grant W911NF-18-1-0103), the U.S. Department of Energy Office of
Science National Quantum Information Science Research Centers (Q-NEXT, award
DE-AC02-06CH11357), the Institute for Quantum Information and Matter, an NSF
Physics Frontiers Center (grant PHY-1125565) with support of the Gordon and
Betty Moore Foundation, the Kavli Nanoscience Institute at Caltech, and the
AWS Center for Quantum Computing. S.M. acknowledges support from the IQIM
Postdoctoral Fellowship.

Disclosures. The authors declare no conflicts of interest.

Supplementary information. is available in the online version of the paper.

Competing interests. The authors declare no competing interests.

Data Availability Statement. Correspondence and requests for materials should
be sent to OP (opainter@caltech.edu).

## Supplementary Section

###  5.1 Experimental setup

Figure 5: AFM integrated with optical test setup. A dimpled fiber positioning
stage is integrated into an AFM, allowing for real-time oxidation and optical
testing. Image on the right shows a close-up view of the AFM cantilever, the
dimpled fiber and the OMC chip.

To enable in-situ nano-oxidation and real-time tuning of OMC cavities, an
optical test setup has been integrated into the enclosure of an AFM as shown
in Fig. 5. The system has three sub-parts that are controlled independently:
1) the optical fiber positioning system, comprising a three-axes linear
motorized positioning stage and two rotational positioning stages, 2) the
sample mount of the AFM with x and y motional degrees of freedom, and 3) the
AFM cantilever with z motional degree of freedom. Additionally, the AFM
enclosure is equipped with a temperature and humidity sensor to ensure that
environmental conditions are within the desired regime. The entire assembly
sits on an optical table to minimize vibrations during the oxidation process.

###  5.2 Characterization of nano-oxidation

Figure 6: Characterization of nano-oxidation on a silicon chip. a, Height
profile of the oxide slab created by nano-oxidation (blue trace), and the
depth profile obtained after removal of the oxide slab using a hydrofluoric
acid etch (red trace). Insets show the AFM scan of the oxide slab (top) and
the etched trench (bottom) along with horizontal line cuts used to sample the
traces on the main plot. b, Height profile of an oxide pixel along x and y
axes indicating a lateral size of $\sim$25 nm. Inset shows the AFM scan of the
oxide pixel along with x and y line cuts used to sample the traces on the main
plot.

The nano-oxidation process generates oxide beneath and above the surface of
silicon. Here we describe the procedure used to assess the ratio of oxide
height and depth above and below the silicon surface, respectively. We first
measure the post-oxidation height by performing AFM scans in the tapping mode.
Subsequently, we etch the oxide in an anhydrous hydrofluoric acid vapor
etching tool (SPTS PRImaxx uEtch), and perform AFM scans of the resulting
trenches in silicon. The data obtained, as depicted in Figure 6a, indicates
that the oxide thickness beneath the surface is approximately 70$\%$ of the
total measured thickness above the surface. Figure 6b shows the
characterization of a single oxide pixel generated in the mild oxidation mode.
The size of the pixel is approximately 25nm along both x and y axes. We
utilize such oxide pixels for fine-tuning of both optical and acoustic
resonances in real-time as discussed in the main text.

###  5.3 Effect of nano-oxidation on optical Q factor

Figure 7: Effect of nano-oxidation on the optical Q factor during coarse
tuning. Change in optical intrinsic Q-factor with respect to the a, width $w$
of the optics-focused oxidation pattern, and b, length $l$ of the acoustics-
focused pattern. Data is obtained in both mild and aggressive oxidation modes.
Error bars represent the standard deviation measured across five OMCs.

In addition to altering the mechanical and optical characteristics of the
local silicon surface, the nano-oxidation process can introduce roughness and
additional scattering sites. Consequently, we observe a reduction in the
optical Q-factor during our coarse tuning experiments. As shown in Figure 7,
the reduction in Q-factor is more pronounced upon increasing the size of the
oxidation patterns and the thickness of the oxide layer. For optics-focused
patterns generated in regions with high electric field intensity, we observe a
substantial impact on the optical Q-factor. A maximum reduction in $Q_{i}$ of
approximately 35% is noted for the aggressive oxidation case when $w<$ 300 nm.
For the mild oxidation case, the reduction in $Q_{i}$ is less than 20%. In the
case of the acoustics-focused pattern, the reduction in $Q_{i}$ is less than
35% for aggressive oxidation case, and less than 20% for mild oxidation case
when $l<$ 1.5 $\mu$m.

###  5.4 FEM simulations for pixel-by-pixel oxidation and selectivity analysis

Figure 8: FEM simulations for pixel-by-pixel oxidation. a, Optical resonance
wavelength shift, $\delta\lambda$ and b, mechanical resonance frequency shift,
$\delta\Omega/2\pi$ simulated as a function of the location of the oxide
pixel. Each oxide pixel is 25x25 nm in size with a thickness of 3.2 nm. c,
Interpolated normalized optical selectivity, $\delta\lambda/\delta\Omega$ and
d, interpolated normalized acoustic selectivity, $\delta\Omega/\delta\lambda$.
The normalization is performed with respect to the maximum shifts obtained in
the pixel-by-pixel oxidation profiles in panels a and b.

Figure 8a and b show the results of FEM simulations where we mapped out how
the location of a single oxidation pixel affects the optical wavelength and
acoustic frequency, respectively. These profiles closely match the optical
intensity and acoustic displacement profiles shown in Fig. 1. This is expected
because oxidation at points with larger optical field concentration leads to
larger change in optical energy, and thereby in the optical resonance
wavelength. Similarly, oxidation at the points of maximum displacement leads
to a larger fractional change in effective mass of the acoustic mode, and
thereby in the acoustic resonance frequency. The exact material properties of
the thin film oxide such as the density, Young’s modulus and relative
permittivity are specific to the oxidation process. In our approach, we use
these material properties as variables and then we find what scaling allows us
to match the coarse tuning experiment (Fig. 2) and simulations. This way, we
found that the density, Young’s modulus and relative permittivity of oxide are
3500 kg/$\text{m}^{\text{3}}$, 70 GPa and 2, respectively. In the acoustic
simulations, we use the anisotropic elasticity tensor of silicon with
($C_{11},C_{12},C_{44}$) = (166, 64, 80) GPa and assume a [110]
crystallographic orientation for the x-axis.

In Figure 8c and d, we plot the normalized optical selectivity,
$\delta\lambda/\delta\Omega$ and normalized acoustic selectivity,
$\delta\Omega/\delta\lambda$. Here, $\delta\lambda$ and $\delta\Omega$ are
normalized with respect to the maximum value calculated in the pixel-by-pixel
oxidation profiles in Figure 8a, b. In the case of fine tuning experiments
where we seek small frequency shifts, we concentrate on the high selectivity
regions. However, the locations with maximum optical selectivity,
$\delta\lambda/\delta\Omega$ do not coincide with those that produce maximum
optical tuning $\delta\lambda$. Likewise, the regions with maximum acoustic
selectivity are relatively small and localized at the edge of the silicon-air
boundary. As a result, for coarse tuning experiments, it is more practical to
use large oxide patches as opposed to individual pixels to achieve large
tuning at the expense of selectivity. Details of the pattern generation
algorithm used to generate the large oxide patches are given in Appendix E.

###  5.5 Nano-oxidation inverse pattern generation algorithm.

Figure 9: Pattern generation matrices for coarse tuning. An example pair of
experimentally determined pattern generation matrices used for aggressive
optical - aggressive acoustic oxidation. a, Optical wavelength shift,
$\delta\lambda$ and b, mechanical frequency shift, $\delta\Omega_{m}/2\pi$ as
a function of the width of the optics-focused component, $w$, and length of
the acoustics-focused component, $l$ of a joint pattern. All other dimensions
and parameters are kept constant. The matrices are generated by spline
interpolation of experimentally measured frequency shift data. c, Pattern
dimension contours corresponding to $\delta\lambda=-1.4$ nm and
$\delta\Omega_{m}/2\pi=-45$ MHz. The intersection of both curves is used to
determine $l$ and $w$ for the joint pattern used for coarse tuning. Inset
shows the joint pattern corresponding to the calculated dimensions.

The results from coarse tuning experiments in Fig. 2 of the main text show
frequency shifts generated by either the optics-focused or the acoustic-
focused pattern type. In practice, we require a union of both types of
patterns to tune both optical and acoustic resonances. To enable the
generation of such joint patterns, we first interpolate results from coarse
oxidation experiments to create an oxidation map with frequency shifts
measured with just the optics-focused and acoustic-focused patterns. Our
inverse pattern generation algorithm uses this map to find a suitable joint
pattern to target the desired simultaneous shift in optical and acoustic
resonances. The pattern search begins with mild oxidation patterns and
progresses to aggressive oxidation patterns if larger frequency shifts are
desired. A pattern is deemed suitable if the calculated resonance shifts can
target the desired optical wavelength to within 200 pm and the desired
acoustic resonance to within 5 MHz. The pattern generation maps presented in
Fig. 9a,b correspond to the aggressive optical - aggressive acoustic oxidation
scenario. There are three other pairs of pattern generation maps corresponding
using mild optical - mild acoustic, mild optical - aggressive acoustic, and
mild optical - aggressive acoustic oxidation scenarios. To illustrate the
pattern search process, Figs. 9a, b show contours corresponding to an example
scenario where shifts, $\delta\lambda=-1.4$ nm and $\delta\Omega_{m}/2\pi=-45$
MHz, are desired. As shown in Fig. 9c, the intersection of these two contours
allows us to estimate the required joint pattern with a width, $w\sim$ 120nm
for the optics-focused component and length, $l\sim$ 860nm for the acoustics-
focused component.

###  5.6 Simultaneous tuning of three OMCs: step-by-step process

Figure 10: Simultaneous optical and acoustic frequency alignment of three
OMCs. Step 0 shows the initial condition of three OMCs. In steps 1-3, we
generate coarse tuning patterns via the inverse pattern algorithm and apply
them to the corresponding OMCs. These patterns are color coded for each OMC
and shown in the right column of the figure. In steps 4-7, we perform real-
time fine tuning of both optical and acoustic resonances. Optical and acoustic
spectra monitored during the nano-oxidation sequence in real-time are shown in
the right column of the figure for each step. In step 7, the final states of e
optical and acoustic resonances on the three devices are shown. The inset
shows a zoomed-in plot of the three optical resonances. This seven step
process allows for optical and acoustic resonance matching to within 2 pm and
500 kHz, respectively.

In this section, we describe in detail the step-by-step tuning protocol used
for frequency alignment of three OMC cavities. This involves a series of
coarse tuning steps followed by real-time fine tuning steps. The experimental
results for each step are shown in Fig. 10.

  1. 1.

Step 0: We started with three OMC cavities (red, magenta, blue) whose optical
wavelengths and acoustic frequencies are $\lambda$ = (1553.891, 1553.381,
1553.123) pm and $\Omega_{m}/2\pi$ = (5.057, 5.047, 5.053) MHz respectively.

  2. 2.

Step 1-3: We feed $\lambda$ and $\Omega_{m}/2\pi$ to the inverse pattern
generation algorithm which generates patterns to tune each OMC to the target
frequencies.

  3. 3.

Step 4-7: At the end of step 3, all 3 OMC cavities are within a suitable range
for fine tuning (typically 200 pm and 5 MHz for optics and acoustic
resonances). We then employ pixel-by-pixel oxidation to fine tune the
resonances while monitoring their values in real-time. The corresponding real-
time wavelengths and frequencies are shown in the right most column of Fig.
10.

At the end of three coarse and four fine tuning steps, we achieve optical
wavelengths of 1553.031 nm, 1553.029 nm and 1553.031 nm and acoustic
frequencies of 5.0261 GHz, 5.0264 GHz and 5.0266 GHz for the three OMCs. This
indicates simultaneous optical wavelength and acoustic frequency matching
within 2 pm and 500 kHz, respectively.

###  5.7 Frequency distribution of fabricated devices

The OMC cavity design in this study incorporates geometric features with
critical dimensions smaller than 100 nm. In our fabrication process, electron
beam lithography and reactive ion etching offer a feature size precision of
approximately 5 nm [1]. In Figure 11a, a scanning electron microscope (SEM)
image of a randomly selected OMC highlights geometric imperfections in the
fabricated device. To assess how these imperfections manifest in the measured
optical and acoustic resonances, we fabricated 20 devices with nominally
identical geometry positioned adjacent to each other on a chip.

Figures 11b and c show the measured optical wavelengths with a standard
deviation of 2.58 nm and acoustic resonance frequencies with a standard
deviation of 19.8 MHz, respectively. The highlighted area in Figures 11b (c),
with a tuning range of 2 nm (60 MHz), signifies the number of devices that can
be brought into resonance solely through optical (acoustic) tuning. The data
indicates that approximately 50% (90%) of OMCs fabricated with the identical
geometric design can be frequency-matched optically (acoustically) after
fabrication.

Figure 11: Resonance frequency distribution of 20 OMCs fabricated with
nominally identical device geometry. a, SEM image of the defect unit cells of
a randomly chosen OMC cavity. Inset on the right shows the magnified view of
an ellipse with irregular boundaries. b, Optical resonance wavelengths and c,
acoustic resonance frequencies of 20 OMC cavities. Shaded areas highlight
cavities that can be independently tuned to the same optical wavelength (red)
or acoustic frequency (blue) with maximum tuning ranges set at 2 nm and 60
MHz, respectively. The data indicates that approximately 50% (90%) of the
devices can be frequency-matched for only optics (acoustics).


# Logic-based AI for Interpretable Board Game Winner Prediction with Tsetlin
Machine

Charul Giri  Centre for AI Research
University of Agder
Grimstad, Norway
<EMAIL_ADDRESS>Ole-Christoffer Granmo  Centre for AI Research
University of Agder
Grimstad, Norway
<EMAIL_ADDRESS>Herke van Hoof  Informatics Institute
University of Amsterdam
Amsterdam, the Netherlands
<EMAIL_ADDRESS>Christian D. Blakely  AI and Real-Time Analytics
PwC Switzerland
Zurich, Switzerland
<EMAIL_ADDRESS>
###### Abstract

Hex is a turn-based two-player connection game with a high branching factor,
making the game arbitrarily complex with increasing board sizes. As such, top-
performing algorithms for playing Hex rely on accurate evaluation of board
positions using neural networks. However, the limited interpretability of
neural networks is problematic when the user wants to understand the reasoning
behind the predictions made. In this paper, we propose to use propositional
logic expressions to describe winning and losing board game positions,
facilitating precise visual interpretation. We employ a Tsetlin Machine (TM)
to learn these expressions from previously played games, describing where
pieces must be located or not located for a board position to be strong.
Extensive experiments on $6\times 6$ boards compare our TM-based solution with
popular machine learning algorithms like XGBoost, InterpretML, decision trees,
and neural networks, considering various board configurations with $2$ to $22$
moves played. On average, the TM testing accuracy is $92.1\%$, outperforming
all the other evaluated algorithms. We further demonstrate the global
interpretation of the logical expressions, and map them down to particular
board game configurations to investigate local interpretability. We believe
the resulting interpretability establishes building blocks for accurate
assistive AI and human-AI collaboration, also for more complex prediction
tasks.

###### Index Terms:

Tsetlin Machine, Winner Prediction, Interpretable AI, Hex

##  I Introduction

Board game winner prediction is a critical part of state-of-the-art AIs for
board game playing, such as AlphaZero [22]. To guide a Monte-Carlo tree search
towards a winning move, for instance, one needs to assess intermediate board
configurations. Lately, state-of-the-art solutions use deep neural network
architectures to evaluate board configurations and to propose which move
sequences to explore. Such black-box algorithms achieve unsurpassed accuracy
in board game position evaluation. However, their complex nature makes them
hard to interpret as they do not reveal their internal decision mechanism.
There are approximate post-hoc methods for explaining specific outputs such as
LIME, yet, being approximate, they can be arbitrarily erroneous [18].

As AI integrates into our daily lives, interpretability becomes increasingly
crucial. Whereas black-box methods can have significant ethical, legal, and
social implications by slipping through biased and wrong decisions [14],
algorithms that explain their output facilitate societal acceptance by being
transparent [15]. An AI-powered recruitment tool developed by Amazon was, for
instance, found to be gender-biased and therefore scrapped [7].

Lack of interpretation is problematic because it hinders quality assurance and
human-computer collaboration. In contrast, an interpretable system enables
human assessment of how the AI reasons, for overseeing the AI decision
process. If the AI, in addition, can provide transparent decision
alternatives, it becomes possible to augment AI-based decision-making with
human intelligence, through collaboration. Indeed, humans and AI collaborating
may be more capable than each alone. As exemplified by Hipp et al. [12], the
combination of the tactical superiority of computers and the strategical
superiority of humans led to the victory of amateur players assisted by AI in
Cyborg Chess111Cyborg Chess, also called Advanced Chess, is a form of chess
where the human player is assisted by a computer chess program to explore
potential moves more effectively., facing grandmasters and supercomputers.

The above findings are in contrast to the traditional view that complex black-
box models are necessary to maximize accuracy, creating a trade-off between
interpretability and accuracy. Further, Rudin suggests that interpretable
models can perform competitively if the representation of the problem is
adequate [18].

Finding competitive interpretable representations for complex board games like
Hex is an unresolved challenge. Contemporary interpretable techniques, like
decision trees and XGBoost, employ greedy learning strategies or rely on model
simplifications. They, therefore, struggle with learning complex patterns such
as board configurations that lead to a win or loss.

In this paper, we propose a competitive interpretable representation for the
board game Hex, learned by a Tsetlin Machine (TM) [9]. TM is a recent
technique for pattern recognition that uses a team of Tsetlin automata [23] to
learn patterns expressed in propositional logic. TMs have obtained competitive
accuracy, memory footprint, and learning speed on several benchmark datasets
[1]. They have been particularly successful in natural language processing,
including explainable aspect-based sentiment analysis [24]. Being based on
finite-state automata, they further support Markov chain-based convergence
analysis [26].

The contributions of the paper can be summarized as follows. We propose an
approach to representing board game configurations in Hex using propositional
logic. The expressions describe the presence and absence of pieces with
conjunctive clauses (AND-rules). This representation allows us to capture
board configurations of arbitrary complexity, inspired by human logical
reasoning [17]. Our approach thus falls in the domain of rule-based
classifiers, which provide humans with easy-to-understand interpretation of
the classifiers’ decision making. We further demonstrate how our
representation can be aggregated and visualized for global and local
interpretation.

##  II Board Game Winner Prediction in Hex

###  II-A Problem Overview

Hex is an abstract two-player strategy board game played on an $n\times n$
rhombus board. In brief, the players alternate with placing pieces on empty
board slots seeking to form a connection from one edge to the opposite edge.
The goal of the game is to form an unbroken chain of pieces from one end of
the board to the other. The first player (Black) needs to connect the top edge
with the bottom edge, while the second player (White) must form a connection
between the left and the right edges.

Predicting the conclusion of the game from any intermediate board state is a
challenging task. For example, a 11$\times 11$ Hex board has 2.4 $\times$
$10^{56}$ legal positions, and a typical board configuration provides almost
three times as many move options compared to Chess[5]. The resulting large
branching factor makes it only possible to do relatively shallow game tree
searches. As the game approaches the end, predicting the winner through search
becomes more feasible, however, from the opening and intermediate board
states, only state-of-the-art neural networks are capable of sufficiently
accurate winner prediction. Note that unlike Go, Hex cannot end in a draw,
producing a two-class prediction problem.

###  II-B Hex Game Strategies

To facilitate the interpretability study later in the paper, we here provide a
brief overview of basic playing strategies for Hex. In Hex, a good offense is
equal to a good defense. If you have completed your connection, it will
automatically mean that you have prevented the opponent from completing her
connection. And, if you have prevented your opponent’s connection, you have,
as a consequence, simultaneously completed your own connection. This is why a
game of Hex never results in a draw [19].

Figure 1: Completed Hex game on an $11\times 11$ board where Black is the
winner

Fig. 1 shows a completed Hex game played on an $11\times 11$ board. As seen,
Black wins the game by forming an unbroken chain from the top edge to the
bottom edge. Notice how the continuous chain from top to bottom prevents White
from forming a continuous chain from left edge to right edge.

In Hex, the first player always has the advantage. This is because in a finite
game, there exists a strategy whereby a first mover cannot lose [20]. The
following swap rule balances this bias.

Swap Rule: After the first player (Black) has made its move, the second player
can choose either to swap sides, taking the first move as its own, or to
respond to the first player, continuing playing as White.

Virtual connections and bridge patterns form the basis for playing
competitively, illustrated in Fig. 2.

Virtual Connection: A virtual connection secures a chain without putting the
pieces directly adjacent to each other, lining out a strong chain. The most
common type of virtual connection is a bridge pattern.

Bridge Pattern: A bridge pattern is created when two pieces of the same color
are placed so that they are separated only by two neighbouring empty hex
cells. In Fig. 3, for example, the link from b2 to c3 is secured by Black,
despite the pieces not being adjacent to each other. If White tries to break
the link by putting a pieces in b3, for instance, Black can still connect
through c2, and vice-versa [4].

Figure 2: Bridge pattern Figure 3: Expansion through bridge pattern

###  II-C Current Approaches

In 2000, the Hexy [2] program won the first Computer Olympiad for Hex. It
exploited alpha-beta search and H-search to find and combine virtual
connections to form more complex patterns. Later, MoHex[3] became champion in
the 2009 Hex Computer Olympiad. MoHex is a Monte Carlo tree search (MCTS)
based approach. It combines MCTS with H-search and ICA to prune the game tree
and speed up the tree search. It further exploits Depth-first Proof Number
search to find the best moves in parallel. MoHex 2.0 [13] is an improvement
over MoHex, using the MiniMax algorithm to learn patterns for weighing the
MCTS simulations.

From 2009 to 2016, different versions of MoHex won all of the Computer
Olympiads. In 2016, however, inspired by the success of AlphaGo [21],
researchers started applying neural networks to the game of Hex. NeuroHex [25]
applied Deep Q-Learning with experience replay. Another program called
MoHexNet [11] used depth-1 trees found by NeuroHEX, and combined the depth-1
trees with MCTS. As a result, MoHexNet won the 2016 Computer Olympiad. Since
then, recent work introduced various ways of improving MoHexNet, leading to
MoHex-CNN [10] and MoHex-3HNN[8], which won the Computer Olympiads in the
following years.

Extensive research has been carried out on the policy evaluation methods for
predicting next best moves, playing out the entire game with MCTS. While a
MCTS search by itself is interpretable, the neural network-based evaluation of
board positions are difficult to interpret due to the lack of transparency of
neural networks. In our present paper, we focus on predicting the winner
directly, from any board position, without playing out the game. Our goal is
to address the lacking transparency of the neural network-based approaches,
hindering for instance human-AI collaboration.

Our paper sheds light on this matter by introducing an interpretable method
for predicting the winner for $6\times 6$ Hex boards. Also, Hex, just like Go,
is a pattern oriented game. So this research can be regarded as a step in the
direction of developing AI assistance for the game of Go as well.

##  III Interpretable Winner Prediction with Logical Rules

In this section, we present the basics of the Tsetlin Machine, to lay the
foundation for our board game winner prediction approach. We further describe
our method for representing winning and losing board game configurations.

Figure 4: A two-action Tsetlin Automaton with $2N$ states.

Figure 5: TM learning dynamics for an XOR-gate training sample, with input
($x_{1}=0,x_{2}=1$) and output target $y=1$.

###  III-A Pattern Representation

A TM in its simplest form takes a feature vector
$\mathbf{x}=[x_{1},x_{2},\ldots,x_{o}]\in\\{0,1\\}^{o}$ of $o$ propositional
values as input and assigns the vector a class $\hat{y}\in\\{0,1\\}$. To
minimize classification error, the TM produces $n$ self-contained patterns.
The input vector $\mathbf{x}$ first provides the literal set
$L=\\{l_{1},l_{2},\ldots,l_{2o}\\}=\\{x_{1},x_{2},\ldots,x_{o},\lnot
x_{1},\lnot x_{2},\ldots,\lnot x_{o}\\}$, consisting of the input features and
their negations. By selecting subsets $L_{j}\subseteq L$ of the literals, the
TM can build arbitrarily complex patterns through ANDing the literal selection
to form conjunctive clauses:

$C_{j}(\mathbf{x})=\bigwedge_{l_{k}\in L_{j}}l_{k}.$ (1)

Above, $j\in\\{1,2,\ldots,n\\}$ refers to a particular clause $C_{j}$ and
$k\in\\{1,2,\ldots,2o\\}$ refers to a particular literal $l_{k}$. As an
example, the clause $C_{j}(\mathbf{x})=x_{1}\land\lnot x_{2}$ consists of the
literals $L_{j}=\\{x_{1},\lnot x_{2}\\}$ and evaluates to $1$ when $x_{1}=1$
and $x_{2}=0$.

###  III-B Representing Board Game Configurations

We represent a board game configuration as a clause, with each literal
referring to whether a piece of a certain color is present or absent in a
specific board location. The following is an example of such a clause,
capturing the $6\times 6$ Hex board pattern visualized in Fig. 6:

$x_{10}\land x_{21}\land\neg x_{6}\land\neg x_{14}\land\neg x_{17}\land\neg
x_{33}\land\neg x_{38}\land\neg x_{51}\land\neg x_{59}\land\neg
x_{60}\land\neg x_{66}\land\neg x_{71}\land\neg x_{72}$.

Features $x_{1},x_{2},\ldots,x_{36}$ refer to the Black pieces, index starting
from the top left of the board and increasing along the rows. For instance,
literal $x_{10}$ in the clause specifies that there must be a Black piece in
position _d2_ , while literal $\neg x_{6}$ signify that a Black piece cannot
be located in location _f1_. Correspondingly, features
$x_{31},x_{32},\ldots,x_{72}$ represent the White pieces. This particular
example pattern provides evidence for Black winning because the positions _c4_
and _d2_ are forming a bridge connection. However, the final prediction is not
based on a single individual pattern but by sharing the knowledge of multiple
learned patterns.

Figure 6: Pattern structure for a sample clause

###  III-C Tsetlin Automata Teams

The fundamental learning unit of TM is a Tsetlin Automaton(TA). The TM assigns
one TA per literal $l_{k}$ per clause $C_{j}$ to build the clauses. The TA
assigned to literal $l_{k}$ of clause $C_{j}$ decides whether $l_{k}$ is
_Excluded_ or _Included_ in $C_{j}$. Fig. 4 depicts a two-action TA with $2N$
states. For states $1$ to $N$, the TA performs action _Exclude_ (Action 1),
while for states $N+1$ to $2N$ it performs action _Include_ (Action 2). As
feedback to the action performed, the environment responds with either a
Reward or a Penalty. If the TA receives a Reward, it moves deeper into the
side of the action. If it receives a Penalty, it moves towards the middle and
eventually switches action.

With $n$ clauses and $2o$ literals, we get $n\times 2o$ TAs. We organize the
states of these in a $n\times 2o$ matrix
$A=[a_{k}^{j}]\in\\{1,2,\ldots,2N\\}^{n\times 2o}$. We will use the function
$g(\cdot)$ to map the automaton state $a_{k}^{j}$ to Action $0$ (_Exclude_)
for states $1$ to $N$ and to Action $1$ (_Include_) for states $N+1$ to $2N$:
$g(a_{k}^{j})=a_{k}^{j}>N$.

We can connect the states $a_{k}^{j}$ of the TAs assigned to clause $C_{j}$
with its composition as follows:

$C_{j}(\mathbf{x})=\bigwedge_{l_{k}\in
L_{j}}l_{k}=\bigwedge_{k=1}^{2o}\left[g(a_{k}^{j})\Rightarrow l_{k}\right].$
(2)

Here, $l_{k}$ is one of the literals and $a_{k}^{j}$ is the state of its TA in
clause $C_{j}$. The logical _imply_ operator $\Rightarrow$ implements the
_Exclude_ /_Include_ action. That is, the _imply_ operator is always $1$ if
$g(a_{k}^{j})=0$ (_Exclude_), while if $g(a_{k}^{j})=1$ (_Include_) the truth
value is decided by the truth value of the literal.

###  III-D Classification

The odd-numbered half of the clauses vote for class $\hat{y}=0$ and the even-
numbered half vote for $\hat{y}=1$. Classification is performed as a majority
vote:

$\hat{y}=0\leq\sum_{j=1,3,\ldots}^{n-1}\bigwedge_{k=1}^{2o}\left[g(a_{k}^{j})\Rightarrow
l_{k}\right]-\sum_{j=2,4,\ldots}^{n}\bigwedge_{k=1}^{2o}\left[g(a_{k}^{j})\Rightarrow
l_{k}\right].$ (3)

As such, the odd-numbered clauses have positive polarity, while the even-
numbered ones have negative polarity. As an example, consider the input vector
$\mathbf{x}=[0,1]$ in the lower part of Fig. 5. The figure depicts two clauses
of positive polarity, $C_{1}(\mathbf{x})=x_{1}\land\lnot x_{2}$ and
$C_{3}(\mathbf{x})=\lnot x_{1}\land\lnot x_{2}$ (the negative polarity clauses
are not shown). Both of the clauses evaluate to zero, leading to class
prediction $\hat{y}=0$.

Input | Clause |  1 |  0
---|---|---|---
Literal | 1 0 | 1 0
Include Literal | P(Reward) |  $\frac{s-1}{s}$ NA | 0 0
P(Inaction) |  $\ \ \frac{1}{s}$ NA |  $\frac{s-1}{s}$ $\frac{s-1}{s}$
| P(Penalty) | 0 NA |  $\ \ \frac{1}{s}$ $\frac{1}{s}$
Exclude Literal | P(Reward) |  0 $\frac{1}{s}$ |  $\ \ \frac{1}{s}$ $\frac{1}{s}$
P(Inaction) |  $\ \ \frac{1}{s}$ $\frac{s-1}{s}$ |  $\frac{s-1}{s}$ $\frac{s-1}{s}$
| P(Penalty) |  $\frac{s-1}{s}$ 0 | 0 0
TABLE I: Type I Feedback
Input | Clause |  1 |  0
---|---|---|---
Literal | 1 0 | 1 0
Include Literal | P(Reward) | 0 NA | 0 0
P(Inaction) | 1.0 NA | 1.0 1.0
| P(Penalty) | 0 NA | 0 0
Exclude Literal | P(Reward) | 0 0 | 0 0
P(Inaction) | 1.0 0 | 1.0 1.0
| P(Penalty) | 0 1.0 | 0 0
TABLE II: Type II Feedback

###  III-E Learning

The upper part of Fig. 5 illustrates learning. A TM learns online, processing
one training example $(\mathbf{x},y)$ at a time. Based on $(\mathbf{x},y)$,
the TM rewards and penalizes its TAs, which amounts to incrementing and
decrementing their states. There are two kinds of feedback: Type I Feedback
encourages clauses to encode frequently-encountered patterns and Type II
Feedback increases the discrimination power of the patterns.

Figure 7: TM learning dynamics for three Hex board positions

Type I feedback is given stochastically to clauses with positive polarity when
$y=1$ and to clauses with negative polarity when $y=0$. Conversely, Type II
Feedback is given stochastically to clauses with positive polarity when $y=0$
and to clauses with negative polarity when $y=1$. The probability of a clause
being updated is based on the vote sum $v$:
$v=\sum_{j=1,3,\ldots}^{n-1}\bigwedge_{k=1}^{2o}\left[g(a_{k}^{j})\Rightarrow
l_{k}\right]-\sum_{j=2,4,\ldots}^{n}\bigwedge_{k=1}^{2o}\left[g(a_{k}^{j})\Rightarrow
l_{k}\right]$. The voting error is calculated as:

$\epsilon=\begin{cases}T-v&y=1\\\ T+v&y=0.\end{cases}$ (4)

Here, $T$ is a user-configurable voting margin yielding an ensemble effect.
The probability of updating each clause is
$P(\mathrm{Feedback})=\frac{\epsilon}{2T}$.

A random sampling from $P(\mathrm{Feedback})$ decides which clauses to update.
The following TA state updates of the chosen clauses’ can be formulated as
matrix additions, subdividing Type I Feedback into feedback Type Ia and Type
Ib:

$A^{*}_{t+1}=A_{t}+F^{\mathit{II}}+F^{Ia}-F^{Ib}.$ (5)

Here, $A_{t}=[a^{j}_{k}]\in\\{1,2,\ldots,2N\\}^{n\times 2o}$ contains the
states of the TAs at time step $t$ and $A^{*}_{t+1}$ contains the updated
state for time step $t+1$ (before clipping). The matrices
$F^{\mathit{Ia}}\in\\{0,1\\}^{n\times 2o}$ and
$F^{\mathit{Ib}}\in\\{0,1\\}^{n\times 2o}$ contains Type I Feedback. A zero-
element means no feedback and a one-element means feedback. As shown in Table
I, two rules govern Type I feedback:

  * •

Type Ia Feedback is given with probability $\frac{s-1}{s}$ whenever both
clause and literal are $1$-valued.222Note that the probability $\frac{s-1}{s}$
is replaced by $1$ when boosting true positives. It penalizes _Exclude_
actions and rewards _Include_ actions. The purpose is to remember and refine
the patterns manifested in the current input $\mathbf{x}$. This is achieved by
increasing selected TA states. The user-configurable parameter $s$ controls
pattern frequency, i.e., a higher $s$ produces less frequent patterns.

  * •

Type Ib Feedback is given with probability $\frac{1}{s}$ whenever either
clause or literal is $0$-valued. This feedback rewards _Exclude_ actions and
penalizes _Include_ actions to coarsen patterns, combating overfitting. Thus,
the selected TA states are decreased.

The matrix $F^{\mathit{II}}\in\\{0,1\\}^{n\times 2o}$ contains Type II
Feedback to the TAs, given per Table III.

  * •

Type II Feedback penalizes _Exclude_ actions to make the clauses more
discriminative, combating false positives. That is, if the literal is
$0$-valued and the clause is $1$-valued, TA states below $N+1$ are increased.
Eventually the clause becomes $0$-valued for that particular input, upon
inclusion of the $0$-valued literal.

The final updating step for training example $(\mathbf{x},y)$ is to clip the
state values to make sure that they stay within value $1$ and $2N$:

$A_{t+1}=\mathit{clip}\left(A^{*}_{t+1},1,2N\right).$ (6)

For example, both of the clauses in Fig. 5 receives Type I Feedback over
several training examples, making them resemble the input associated with
$y=1$.

###  III-F Learning Board Configuration Patterns

((a)) Type Ia Feedback

((b)) Type Ib Feedback

((c)) Type II Feedback

Figure 8: Visualization of Tsetlin Automata state updates

Fig. 7 demonstrates how Type Ia, Type Ib, and Type II Feedback operates on
board game configurations to update a clause. We use $3\times 3$ boards to
visualize learning. The starting point is an empty clause, which is updated
based on an input board that led to a Black win (upper left). Since the clause
is empty, it evaluates to True. Type Ia Feedback then makes the clause
resemble the input board by introducing literals from the board, building a
pattern for winning positions. The next new winning board configuration given
as input does not match the clause, leading to Type Ib Feedback that removes
literals from the clause. Finally, a losing board position is given as input.
Since the clause evaluates to True, Feedback Type II is given, refining the
clause to become more discriminative (upper left). It is the states of the
underlying Tsetlin Automata that drive these changes, illustrated by Fig. 8).

##  IV Empirical Evaluation

In this section we evaluate our proposed method and demonstrate how to
interpret the board configuration patterns produced.

###  IV-A Dataset

The dataset we use for evaluation consists of $287,000$ game board
configurations, which were created through self-play using a modified version
of MoHex2.0, with MCTS and data augmentation. The games lasted for up to 28
moves, with the majority of games ending before 20 moves (Fig.
LABEL:fig:datastat).

Since TM requires input features to be binarized, each board is represented
with a $72$ bit vector, as described in the previous section.

Class | Game Boards
---|---
Black |  175968
White |  111826
Total |  287794
TABLE III: Dataset distribution Figure 9: Distribution of games per number of
moves played

###  IV-B Results

Our experimental setup is as follows. As previously explained in Section III,
a TM has two hyper-parameters $T$ and $s$, along with the number of clauses
$n$. For our experiments we initialized our TM with $10,000$ clauses and set
the other hyperparameter $T$ and $s$ to $8000$ and $100$, respectively. As a
general rule, TMs perform well with a $T$ value equal to $80\%$ of the total
number of clauses per class, leaving only $s$ to be tuned. We found an
appropriate $s$ value through a limited manual trial and error search, so
performance may possibly improve trough a more thorough automated grid search.
We used $67\%$ of the data to train the TM for $200$ epochs, and the remaining
$33$% of the data for evaluating prediction accuracy.

Table IV contains the results of our experiments, including a comparison with
six popular machine learning algorithms. For a fair comparison, we used the
same training and testing datasets for all of the algorithms. Each algorithm
was first trained using default sklearn parameters, and the top-performing
ones were then fine-tuned to maximize performance. As seen in the table, we
obtained an average testing accuracy of $92.1\%$ with our approach,
outperforming the other interpretable methods like Decision Trees, XGBoost[6],
and InterpretML[16] by a wide margin.

Also, our technique performs better than a fine-tuned black-box neural network
(NN) configuration. The NN consisted of five layers with Sigmoid activation
functions, trained with the Adam optimizer. We further utilized binary
Crossentropy loss function, and exploited $L_{2}$ regularization and dropout
to reduce overfitting. Still, the maximum accuracy we obtained with the NN was
about 1.5% lower than our TM-based approach.

To investigate whether accuracy varies depending on the stage of a game (e.g.,
beginning, middle, or end), we grouped the data according to the number of
moves played (Fig. 9). Fig. 10 summarises how training and testing accuracy
varies by the number of moves played. Observe how the performance of TM
remains more or less unaffected by number of moves played, indicating a strong
ability to predict the winner at all stages of a game. That is, our method is
not only able to predict the winner accurately later in a game, when more
information is available (intermediate and ending games), it can also predict
the winner very early in the game, signifying the importance of the first
moves played.

Earlier in this paper we discussed that TM learns patterns by forming clauses,
so we performed a precision analysis of the clauses. Fig. 11 illustrates the
precision histograms for positive polarity and negative polarity clauses. It
can be observed that a majority of the clauses fall in a fairly precise group,
which signify that the predictions are made with a higher confidence.

Method | Hyperparameter | Training Accuracy | Testing Accuracy
---|---|---|---
Naïve bayes |  default |  72.02 |  72.02
Decision Tree |  max depth=100 |  100 |  87.1
K-nearest neighbour |  k=1 |  100 |  84
xgboost |  |  84.88 |  88.90
interpretML(EBC) |  default |  74.56 |  74.33
Neural Network |  |  97.6 |  90.7
Tsetlin Machine |  Clauses = 10000 T = 8000 S = 100 max weight=255 |  95.6 |  92.1

TABLE IV: Comparison with other machine learning techniques

((a)) Training accuracy per #moves played

((b)) Testing accuracy per #moves played

Figure 10: Accuracy at different stages of a game

((a)) Positive Polarity Clauses

((b)) Negative Polarity Clauses

Figure 11: Clause precision

###  IV-C Interpretation

We now investigate how to visualize the clauses for global and local
interpretation.

####  IV-C1 Global Interpretability

Global interpretability [15] describes how a model behaves in general. It
takes a holistic view of the features, providing insight into what features
are informative overall and how they interact in the model.

In a game of Hex, analysing the interaction of features is crucial for
understanding what constitutes a winning pattern. Each clause a TM learns
captures one particular feature interaction, as displayed in Fig. 12. The
figure contains the ten most impactful positive polarity clauses (Fig. 12(a))
and the ten most impactful negative ones (Fig. 12(b)), visualized as board
patterns. The positive polarity clauses capture winning patterns for Black,
while the negative polarity clauses represent losing patterns.

((a)) Positive polarity clauses

((b)) Negative polarity clauses

Figure 12: Ten most impactful clauses per polarity

The clauses $C_{j}$ have been selected by combining precision and data
coverage:

$Score(j)=Precision(j)^{\alpha}\times Coverage(j)$ (7)

calculated from the True Positive (TP), False Negative(FN), and False Positive
(FP) rates:

$\displaystyle Precision=$ $\displaystyle\frac{TP}{TP+FP}$ (8) $\displaystyle
Coverage=$ $\displaystyle\frac{TP}{TP+FN}.$ (9)

The $\alpha$ parameter is for weighing precision against coverage. That is, a
higher $\alpha$ puts more emphasis on clause precision. Since a majority of
the clauses provides high coverage, they tend to have relatively lower
precision. Therefore, we use $\alpha=10$ to single out clauses with high
precision, yet providing reasonable coverage.

((a)) Positive polarity clause

((b)) Negative polarity clause

Figure 13: Example of clause interpretation

Consider, for example, the negative polarity clause in Fig. 13(b). This clause
predicts loss for Black due to the string of negated Black pieces from _a3_ to
_f2_. Hence, the clause states that Black at this point not yet has formed a
continuous chain from top to bottom.

Now, observe the positive polarity clause in Fig. 13(a). This clause predicts
a win for Black. As can be seen from the board, White lacks pieces from _a4_
to _c6_. Accordingly, the pattern opens up for a potential chain for Black.
The missing white pieces also provide evidence for a possible bridge between
_a4_ and _b2_ , increasing the probability of a Black chain forming before the
White one.

####  IV-C2 Local Interpretability

Local interpretability explains specific data points, as opposed to global
interpretation which characterises the entire model at once [15]. In our case,
global interpretation also forms the basis for local interpretation. That is,
by combining the pattern information from the individual clauses that evaluate
to True for a specific board configuration, more insights about the board
configuration is obtained.

Since TM clauses work together to form an accurate prediction, it can in some
cases be insufficient to interpret a single clause alone. Instead, we
aggregate the information from all the clauses that took part in classifying
the input. Then we capture exactly which literals the TM relied on for the
given classification. Consider the example input board in Fig. 14(a). White
went on and won the game from this position, as predicted by the TM. Fig.
14(b) interprets this TM prediction through aggregating the True clauses. That
is, each Hex cell reports the number of times a Black or White piece is
assigned to the cell by one of the clauses that are True for the given board
position (only considering non-negated features).

We can see that a higher count is given to Hex cell _b5_. The high count
captures that _b5_ forms a bridge to _c3_ as well as a secure edge link to
_a5_ and _a6_ , which completes the continuous chain for White. It is
impossible for Black to break this chain. If Black tries to block White by
playing _b4_ , White can continue the chain by playing _c4_ , or the other the
way around. In the next step, if Black tries to block White by playing _a6_ ,
White can complete the chain by playing _a5_ , and vice-versa.

Also notice the pattern forming in the lower right of Fig. 14(b), with a lower
count. Even though White can build a chain here as an alternative, it is a
weaker option since Black can block White’s path from _f1_ to _f6_ simply by
playing _f6_.

In conclusion, this example demonstrate that the TM not only recognizes
individual patterns but is also able to aggregate the patterns to assess
chance of winning.

((a)) Sample test input board

((b)) Local interpretation of clauses for board in Fig. 14(a)

Figure 14: Local aggregated interpretation

##  V Conclusion

In this paper, we presented an interpretable winner prediction approach for
the game of Hex. By using simple AND rules in propositional logic, we can
specify strong and weak board configurations. Our empirical results dispute
the widespread assumption that one cannot escape trading off interpretability
against predictive performance. Despite being inherently interpretable, our
approach outperforms several other popular ML methods, achieving an accuracy
of $92$% on test data. Straightforward board configuration patterns enable
global and local interpretation, through visualization. By aggregating the
active clauses, we can further assess the importance of specific piece
formations, revealing the reasons behind our model’s predictions.

In our further work we plan to go from winner prediction to predicting the
winning move, to solve the game of Hex on larger boards in an interpretable
manner.

## References

  * [1] Kuruge Darshana Abeyrathna, Bimal Bhattarai, Morten Goodwin, Saeed Rahimi Gorji, Ole-Christoffer Granmo, Lei Jiao, Rupsa Saha, and Rohan K Yadav.  Massively parallel and asynchronous tsetlin machine architecture supporting almost constant-time scaling.  In International Conference on Machine Learning, pages 10–20. PMLR, 2021.
  * [2] Vadim V Anshelevich.  The game of hex: An automatic theorem proving approach to game programming.  In AAAI/IAAI, pages 189–194, 2000.
  * [3] Broderick Arneson, Ryan Hayward, and Philip Henderson.  Mohex wins hex tournament.  ICGA journal, 32(2):114, 2009.
  * [4] Cameron Browne.  Hex Strategy: Making the right connections.  CRC Press, 2009.
  * [5] Stephan Chalup, Drew Mellor, and Frances Rosamond.  The machine intelligence hex project.  Computer Science Education, 15, 07 2005.
  * [6] Tianqi Chen and Carlos Guestrin.  Xgboost: A scalable tree boosting system.  In Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, pages 785–794, 2016.
  * [7] Jeffrey Dastin.  Amazon scraps secret ai recruiting tool that showed bias against women, 2018.
  * [8] Chao Gao, Kei Takada, and Ryan Hayward.  Hex 2018: Mohex3hnn over deepezo.  J. Int. Comput. Games Assoc., 41(1):39–42, 2019.
  * [9] Ole-Christoffer Granmo.  The tsetlin machine - A game theoretic bandit driven approach to optimal pattern recognition with propositional logic.  CoRR, abs/1804.01508, 2018.
  * [10] Ryan Hayward and Noah Weninger.  Hex 2017: Mohex wins the 11x11 and 13x13 tournaments.  ICGA Journal, 39(3-4):222–227, 2017.
  * [11] Ryan Hayward, Noah Weninger, Kenny Young, Kei Takada, and Tianli Zhang.  Mohex wins 2016 hex 11x11 and 13x13 tournaments.
  * [12] Jason Hipp, Thomas Flotte, James Monaco, Jerome Cheng, Anant Madabhushi, Yukako Yagi, Jaime Rodriguez-Canales, Michael Emmert-Buck, Michael C Dugan, Stephen Hewitt, et al.  Computer aided diagnostic tools aim to empower rather than replace pathologists: Lessons learned from computational chess.  Journal of pathology informatics, 2, 2011.
  * [13] Shih-Chieh Huang, Broderick Arneson, Ryan B Hayward, Martin Müller, and Jakub Pawlewicz.  Mohex 2.0: a pattern-based mcts hex player.  In International Conference on Computers and Games, pages 60–71. Springer, 2013.
  * [14] Davinder Kaur, Suleyman Uslu, Arjan Durresi, Sunil Badve, and Murat Dundar.  Trustworthy explainability acceptance: A new metric to measure the trustworthiness of interpretable ai medical diagnostic systems.  In Conference on Complex, Intelligent, and Software Intensive Systems, pages 35–46. Springer, 2021.
  * [15] Christoph Molnar.  Interpretable machine learning.  Lulu. com, 2020.
  * [16] Harsha Nori, Samuel Jenkins, Paul Koch, and Rich Caruana.  Interpretml: A unified framework for machine learning interpretability.  arXiv preprint arXiv:1909.09223, 2019.
  * [17] I. Noveck, R. B. Lea, George M. Davidson, and D. O’brien.  Human reasoning is both logical and pragmatic.  Intellectica, 11:81–109, 1991.
  * [18] Cynthia Rudin.  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  Nature Machine Intelligence, 1(5):206–215, 2019.
  * [19] Mark Schachner.  The game of hex: a study in graph theory and algebraic topology.  Mathematics REU, 2019.
  * [20] Ulrich Schwalbe and Paul Walker.  Zermelo and the early history of game theory.  Games and economic behavior, 34(1):123–137, 2001.
  * [21] David Silver, Aja Huang, Chris J Maddison, Arthur Guez, Laurent Sifre, George Van Den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Veda Panneershelvam, Marc Lanctot, et al.  Mastering the game of go with deep neural networks and tree search.  nature, 529(7587):484–489, 2016.
  * [22] David Silver, Thomas Hubert, Julian Schrittwieser, Ioannis Antonoglou, Matthew Lai, Arthur Guez, Marc Lanctot, Laurent Sifre, Dharshan Kumaran, Thore Graepel, et al.  A general reinforcement learning algorithm that masters chess, shogi, and go through self-play.  Science, 362(6419):1140–1144, 2018.
  * [23] Michael Lvovitch Tsetlin.  On behaviour of finite automata in random medium.  Avtomat. i Telemekh, 22(10):1345–1354, 1961.
  * [24] Rohan Kumar Yadav, Lei Jiao, Ole-Christoffer Granmo, and Morten Goodwin.  Human-level interpretable learning for aspect-based sentiment analysis.  In The Thirty-Fifth AAAI Conference on Artificial Intelligence (AAAI-21). AAAI, 2021.
  * [25] Kenny Young, Gautham Vasan, and Ryan Hayward.  Neurohex: A deep q-learning hex agent.  In Computer Games, pages 3–18. Springer, 2016.
  * [26] Xuan Zhang, Lei Jiao, Ole-Christoffer Granmo, and Morten Goodwin.  On the Convergence of Tsetlin Machines for the IDENTITY- and NOT Operators.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021\.


<EMAIL_ADDRESS>
Programing Using High Level Design With Python and FORTRAN: A Study Case in
Astrophysics

Eduardo S<EMAIL_ADDRESS>Oswaldo D. Miranda
<EMAIL_ADDRESS>
Instituto Nacional de Pesquisas Espaciais - Divisão de Astrofísica,
Av. dos Astronautas 1758, São José dos Campos, 12227-010 SP, Brazil

Received on September , 2011 / accepted on *****, 2010

###### Abstract

In this work, we present a short review about the high level design
methodology (HLDM), that is based on the use of very high level (VHL)
programing language as main, and the use of the intermediate level (IL)
language only for the critical processing time. The languages used are Python
(VHL) and FORTRAN (IL). Moreover, this methodology, making use of the oriented
object programing (OOP), permits to produce a readable, portable and reusable
code. Also is presented the concept of computational framework, that naturally
appears from the OOP paradigm. As an example, we present the framework called
PYGRAWC (Python framework for Gravitational Waves from Cosmological origin).
Even more, we show that the use of HLDM with Python and FORTRAN produces a
powerful tool for solving astrophysical problems.

Keywords: Computational Physics, Cosmology, Programming Methodology, HLDM.

1\. INTRODUCTION

Python [1, 2] is a very-high level dynamic language programming. The Python
interpreter is available for different operational systems (OS). This means
that is possible to write a code which can run in different OS without
requiring any modification. However, in order to have the best performance,
the critical computational part of the code should be written in a compiled
language like c/c++ or FORTRAN. This way of writing codes is called of high
level design (HLD) [3, 4]. Another important fact is that, in general, $80$%
of the runtime is spent in $20$% of the code (Pareto Principle) [5]. Thus, the
use of a VHL for the principal part of the code permits a more agile
processing. Although the HLD goes beyond of a mixing of different languages,
an important feature of HLD is related to the use of oriented object
programming paradigm, OOP, working together an Unified Modeling Language (UML)
class diagram.

2\. The High Level Design

The first step to write an efficient code consists in dividing the problem in
classes. In this section we present programs that make this task easier, as
for example, dia [6] and dia2code [7]. The first one permits us to write class
diagrams, and with the second program we can generate a frame code. That is,
the class in the diagram image which is converted to a class in code
structure. Through this section these concepts will become clearer.

2.1 Planning Before Programing

We start this section showing an example of class diagram333All documentation
about how to install and use the dia software can be found in [6].. Through
this paper, we will use an example derived from cosmology. In particular, the
main characteristics of a cosmological model are: the age of the Universe, the
scale factor which describes how the radius of the Universe evolves with time,
and the density of matter/energy.

The class diagram that represents this cosmological model is showed in the
figure 1. The attributes of this class are the cosmological parameters, at the
present, for total matter (self.omegam - $\Omega_{m}$), barionic matter
(self.omegab -$\Omega_{b}$), dark energy (self.omegal - $\Omega_{\Lambda}$)
and Hubble parameter (self.h- $h$ )444The Hubble constant at the present time
is written in terms of $h$ by $H_{0}=100\,h\,{\rm km}\,{\rm s}^{-1}\,{\rm
Mpc}^{-1}$ (where $1\,{\rm Mpc}=3.086\times 10^{24}\,{\rm cm}$)..

The file is saved with the name cosmo.dia. Now, it is possible to generate a
structured code with dia2code (see [7] for details), using the command
dia2code cosmo.dia -t Python. The code generated by this example, and all
examples used in this article, can be downloaded from [8]. It is possible to
generate c++ and java structured code choosing the name of the equivalent
language from the dia2code command.

However, this way to structure a code class is only a start point. It is
necessary to do a better organization and fill the methods with the equivalent
operations.

Figure 1: cosmo.dia, an example of class diagram for the basic characteristics
of a cosmological model.

2.2 Mixing Python and FORTRAN

There is a very useful tool, called f2py [9], which permits to do a wrapper of
a FORTRAN 77 code to Python. That is, it compiles the FORTRAN subroutine in a
format which can be used by Python module. The f2py is contained in the
package numpy [10]. Below, we present a simple example:

C FILE hiword.f

subroutine hiword(a,b)

real*8 a,b

cf2py intent(in) a

cf2py intent(out) b

b = a*a

write(6,*) ’b = ’,b,’, a = ’,a

return

end

The comment cf2py allows the f2py wrapper can be identified with both the
input and output variables in the function hiword. Giving the name hiword.f to
the file contained in the above code, we can compile it from the following
command:

f2py -c hiword.f -m hiword

In this case, the -c means compile, and -m generate a Python module with name
hiword. Below, we present an example how to call the function hiword in a
Python code.

$[1]>>>$import hiword

$[2]>>>$print hiword.__doc__

$[3]>>>$ This module ’hiword’ is auto-generated with f2py (version:2).

$[4]>>>$Functions:

$[5]>>>$ b = hiword(a)

$[6]>>>$hiword.hiword(5)

$[7]>>>$ b = 25.000000000000000 , a = 5.0000000000000000

$[8]>>>$ 25.0

The text in front of $>>$ represents what is printed in the display. For more
details and examples see [9].

2.2 Optimizing the Code for Multi-Core Machines

Another interesting fact about Python is that it has a lot of modules. One of
this is the multiprocessing that permits to write a parallel code in an easy
way. As an example of using this module in scientific computing, consider the
following equation:

$f(x)=\int_{a}^{b}{g(x,k)}dk,$ (1)

where $a\leq k\leq b$.

In many cases $g(x,k)$ can not be written in a separated form. In this case,
the integral equation must be evaluated for each $x$ in a given range
$[x_{0},x_{f}]$. However, we can divide the range $[x_{0},x_{f}]$ by the
number of central processor units (CPU) of a cluster compute (or multi-core
machine), and so we can calculate $f(x)$ in parallel mode. In the figure 2, it
is showed the class diagram of ppvector, that is a class we developed to do
this type of operation in parallel model, for multi-core machine, based on the
module multiprocessing. The source code can be downloaded from [11].

Figure 2: ppvector, a Python module for construction of parallel scientific
code in a multi-core machine.

The code below shows the use of ppvector:

import multiprocessing as mpg

from ppvector import ppvector

from scipy.integrate import romberg

np=10000; zmax=20.0; deltaz=zmax/np

g= mpg.Array(’d’,[0 for i in range(np)]) # The d indicate duble precision

z= mpg.Array(’d’,[zmax-i*deltaz for i in range(np)])

#Define a function that will be calculate the integral in parallel

#k is the starter point of the sub-range

#E is the lenght of the range

#n is the number of CPU’s of machine

def f(x):

def f2(k):

return (x+k)**(-2.0)

return romberg(f2,5.0,20.0)

def fun(k,E,n):

k2=k+E

for i in range(k,k2+1):

zloc=z[k]

g[k]=f(zloc)

C1= ppvector(np,fun) # Star the ppvector class

C1.runProcess() # Executing the parallel calculus.

The function Array, of multiprocessing module, allocates a matrix in a global
memory which can be accessed by all CPU’s. In line $25$ is passed on the
length of the vector and the function that divides the job in sub-ranges. In
line $26$, the parallel code is called and executed.

3\. Python Framework for Cosmological Gravitational Waves - PYGRAWC

A framework is a set of classes, interfaces and patterns to solve a group of
problems. It is like a little application with statical and dynamical
structures to solve a set of restrict problems. So, a framework is more than a
simple library (we refer the reader to [12, 13, 14]).

In figure 3 is presented the class diagram of the core of PyGraWC. It is a
framework that we are developing to study gravitational waves from
cosmological origin. Here, it is only showed the class name and the relation
among their several components.

Figure 3: The PYGRAWC framework class diagram.

The class cosmo describes the background cosmology. The class PressSchechter
is based on a Press-Schechter-like formalism [15] and it describes both the
evolution of dark matter halos and the infall of barionic matter in these
halos. The class csfr describes the evolution of the cosmic star formation
rate. The class smbh describes the evolution of supermassive black holes in
the centers of galaxies. The classes bhestelar and bhmassivo calculate the
stochastic background of gravitational waves generated by: the collapse of
stars to form black holes [15] and the growth of supermassive black holes (in
progress). All details about the astrophysical model and the results obtained
from this framework can be seen in [15, 16, 17].

4\. Final Considerations

In this work is presented a High Level Design methodology (HLD), that consists
in the mixing of a very-hight level interpreted language (VHL) with an
intermediated compiled language (IL). Using tools of software engineering,
like UML, and also framework concept, we can write efficient scientific codes
without spending a lot of time in the development phase. Here, it was used
Python (VHL) and Fortran (IL) and it was showed that this combination can be
easily done giving excellent results, as can be seen by the presentation of
Python framework for Gravitational Waves from Cosmological origin (PyGraWC).

Acknowledgments: E.S. Pereira would like to thank the Brazilian Agency CAPES
for support. O.D. Miranda would like to thank the Brazilian Agency CNPq for
partial support (grant 300713/2009-6)

## References

  * [1] http://www.Python.org/ $<$accessed in 09 September 2011$>$
  * [2] GUPTA, R, 2002. Making Use of Python. Wiley Publishing.
  * [3] HINSEN, K, LANGTANGEN, HP, SKAVHAUG, O, ØDEGÅRD, Å, 2006. Using BSP and Python to simplify parallel programming. Future Generation Computer Systems. 22: 123-157.
  * [4] HINSEN, K, 2007. Parallel scripting wiht Python. Computing in science & engineering. Nov./Dec. 82-89.
  * [5] BEHNEL, S, BRADSHAW, R, CITRO, C, DALCIN, L, SELJEBOTN, D, SMITH, K,.Cython: The Best of Both Worlds. CISE, 13, 2, 31-39.
  * [6] http://projects.gnome.org/dia/ $<$accessed in 09 September 2011$>$
  * [7] http://dia2code.sourceforge.net/ $<$accessed in 09 September 2011$>$
  * [8]<EMAIL_ADDRESS>$<$accessed in 10 September 2011$>$
  * [9] PETERSON,P, 2009. F2PY: a tool for connecting Fortran and Python programs IJCSE, 4, 4, 296-305.
  * [10] http://www.scipy.org/ $<$accessed in 09 September 2011$>$
  * [11]<EMAIL_ADDRESS>$<$accessed in 09 September 2011$>$
  * [12] FAYAD, M. E. 2000. Introduction to the computing surveys electronic symposium on object-oriented application frameworks. ACM Comput. Surv., 32,1,1-9.
  * [13] Fayad, M. E., SCHIMIDT, D. C. 1997. Object-oriented application frameworks. Commun.ACM, 40, 10, 32-38.
  * [14] GOVONI, D. 1999. Java Application Frameworks. John Wiley & Sons.
  * [15] PEREIRA, ES, MIRANDA, OD, 2011. Supermassive Black Holes: Connecting the Growth to the Cosmic Star Formation Rate. accept in Mon. Not. R. Astron. Soc. Letter.
  * [16] PEREIRA, ES, MIRANDA, OD, 2010. Stochastic background of gravitational waves generated by pre-galactic black holes. Mon. Not. R. Astron. Soc. 401: 1924-1932.
  * [17] PEREIRA, ES, MIRANDA, OD, 2010. Massive Black Hole Binary Systems in Hierarchical Scenario of Structure Formation. International Journal of Modern Physics D. 19: 1271-1274