text:
# Proving correctness of an iterative Fibonacci algorithm

One of the questions in the problem sets that I'm struggling in is this specific number that asks me to prove an iterative Fibonacci algorithm. The algorithm is written below:

function fib(n)
if n = 0 then
return(0)
else
a = 0; b = 1; i = 2;
while i <= n do
c = a + b
a = b
b = c
i = i + 1
return(b)

The way too prove correctness, according to my professor was to make sure that there are these three steps:

1. Initialization - the loop invariant must hold true prior to the first iteration
2. Maintenance - the loop invariant must hold true after an iteration
3. Termination - the loop invariant must hold true when the loop terminates

The loop invariant I've chosen is a <= b since I find this to be true for steps 1 through 3. First of all, I'm not sure if this is a valid loop invariant and this is the only observation I saw since i <= n isn't always true for inputs n that are natural numbers.

Assuming that I've chosen the correct loop invariant, I need to answer the proof by doing three steps so for this number I plan to answer it this way

1. Initialization - before the start of the loop a is assigned a value of 0 while b is assigned a value of 1 which starts the Fibonacci sequence. a <= b holds true prior to the start of the loop
2. Maintenance - during the loop, another variable c is added such that it is equal to the sum of a and b. After which, b is assigned to variable a while c is assigned to b thus making the invariant a <= b true during the iteration.
3. Termination - the loop ends when i > n. Before i is incremented, the procedures in the maintenance step is still done, thus the loop invariant still holds true

For my questions, is my loop invariant a <= b correct? And are the three statements I mentioned above sufficient to prove the correctness of an iterative Fibonacci algorithm?

• Although I'm not a computer scientist, I stongly suspect that your loop invariant must be related to Fibonacci numbers. Something along the lines of b being the (i-1)th Fobinacci number. Couldn't you else "prove" that your code computes anything correctly? Feb 15 '18 at 18:31
• This is the part where I am a bit confused. Because, from what I get from the proof definition, to prove that this is correct, the loop invariant must hold true for the three cases. I am really not sure if it has something to do with correct computation. Feb 15 '18 at 18:38
• What exactly is your question? If it is "Is a<=b a loop invariant in my algorithm?" then the answer is yes. If your question is "Does this loop invariant help me prove the algorithm correct?" then the answer is no. Feb 15 '18 at 18:52

The condition $a \leq b$ is a loop invariant, but it captures very little of the operation of the loop. A better loop invariant is $$a = F_{i-2} \text{ and } b = F_{i-1}.$$ This clearly holds just before the loop. The loop sets $a' = b = F_{i-1}$, $b' = a+b = F_{i-2}+F_{i-1}=F_i$, $i' = i+1$, from which it is easy to prove that it maintains the loop invariant. When the loop has finished, we must have $i = n+1$, and so $a = F_{i-2} = F_{n-1}$ and $b = F_{i-1} = F_n$.

There is a weak point in this argument: How do we know that $i = n+1$ at the end of the loop? We can fix this by strengthening the loop invariant to $$a = F_{i-2} \text{ and } b = F_{i-1} \text{ and } i \leq n+1.$$ This holds just before the loop since $n \geq 1$ (assuming $n \geq 0$, a precondition that needs to be added to the entire function). It is maintained during the loop since at the beginning of each iteration $i \leq n$. Finally, after the loop we have $i \leq n+1$ and $i > n$, which forces $i = n+1$ since $i$ is an integer.

• Now I see the part where the computation of the Fibonacci takes place. I was thinking of something like this but I just couldn't properly sew it in. I guess the part of determining the loop invariant is the most important part for proving iterative algorithms Feb 15 '18 at 19:44


text:
## Biquandle brackets

Sam Nelson (CMC)

### Abstract

Given a finite biquandle X and a commutative ring with identity R, we define an algebraic structure known as a biquandle bracket. Biquandle brackets can be used to define a family of knot and link invariants known as quantum enhancements which include biquandle cocycle invariants and skein polynomials such as the Alexander, Jones and HOMFLYpt polynomials as special cases. As an application we will see a new skein invariant which is not determined by the knot group, the knot quandle or the HOMFLYpt polynomial.

### Where

Millikan 2099, Pomona College

## A tale of two matrices

### Speaker

Stephan Garcia (Pomona College)

### Abstract

Let $A$ and $B$ be $n \times n$ matrices. Although they are typically unequal, $AB$ and $BA$ share many important properties. We survey some known results in the area and discuss a few novel theorems. This talk will be accessible to students. This is joint work with David Sherman (U.~Virginia) and Gary Weiss (U.~Cincinatti).

### Where

Millikan 2099, Pomona College

## Few distinct distances implies no heavy lines

### Abstract

After Guth and Katz's almost tight bound for the distinct distances problem, one of the main open problems is to characterize the structure of planar point sets that determine a small number of distinct distances. We show that if a set of n points determines o(n) distinct distances then no line contains n^{7/8} points of the set and no circle contains n^{5/6} such points. Our analysis combines tools from incidence geometry and additive combinatorics. In the talk, before getting to our result I will spend some time on surveying the distinct distances problem in general. Joint work with Joshua Zahl and Frank de Zeeuw.

### Where

Millikan 2099, Pomona College

## Spherical designs and lattices

### Speaker

Hiren Maharaj (Pomona College)

### Abstract

Let n > 1. A collection of points P on the unit sphere S_{n-2} in R^{n-1} is called a spherical t-design for some positive integer t if the average value of every polynomial in n-1 variables with real coefficients of degree t or less on  S_{n-2} equals the average value of f on the set P. A full-rank lattice in R^{n-1} is called strongly eutactic if its set of normalized minimal vectors forms a spherical 2-design. Given a finite Abelian group G of order n, one can form a corresponding lattice L_G of rank n-1. This talk will be about recent joint work with Albrecht Boettcher, Lenny Fukshansky and Stephan Garcia, in which we show that L_G is strongly eutactic if and only if n is odd or G=(Z/2Z)^k for some positive integer k.

### Where

Millikan 2099, Pomona College

## Belyi maps on elliptic curves and dessin d'enfant on the torus

### Speaker

Edray Goins (Purdue University)

### Abstract

A Belyi map $\beta: P^1(C) \to P^1(C)$ is a rational function with at most three critical values; we may assume these values are ${0, 1, \infty }$. A Dessin d'Enfant is a planar bipartite graph obtained by considering the preimage of a path between two of these critical values, usually taken to be the line segment from 0 to 1. Such graphs can be drawn on the sphere by composing with stereographic projection: $\beta^{-1} ([0,1]) \subseteq P^1(C) \simeq S^2(R)$. Replacing $P^1$ with an elliptic curve $E$, there is a similar definition of a Belyi map $\beta: E(C) \to P^1(C)$. The corresponding Dessin d'Enfant can be drawn on the torus by composing with an elliptic logarithm: $\beta^{-1} ([0,1]) \subseteq E(C) \simeq T^2(R)$. In this talk, we discuss the problems of (1) constructing examples of Belyi maps for elliptic curves and (2) drawing Dessins d'Enfants on the torus. This work is part of PRiME (Purdue Research in Mathematics Experience) with Leonardo Azopardo, Sofia Lyrintzis, Bronz McDaniels, Maxim Millan, Yesid Sanchez Arias, Danny Sweeney, and Sarah Thomaz with assistance by Hongshan Li and Avi Steiner.

### Where

Millikan 2099, Pomona College

## Building a better Bernoulli Factory

Mark Huber (CMC)

### Abstract

Suppose I have a non-fair coin with unknown probability $p$ of heads that I can flip as many times as I'd like. A Bernoulli factory is a way of flipping the coin to construct a single coin flip that has chance $f(p)$ of heads for some function $f(p)$. For example, to get a $p^2(1-p)$-coin, I flip the original coin 3 times, and count it as heads if the sequence is HHT, and tails otherwise. My goal in this talk is to build a fast Bernoulli factory for the innocent looking function $f(p) = Cp$, where $C$ is a given constant. It might look innocent, but it turns out that for a fixed number of coin flips, you cannot build a Bernoulli factory that gives a $Cp$ coin for $C > 1$.

So we're going to use a random number of flips to build our $Cp$-coin. I'll show that when $Cp$ is small, this new method uses (to first order) $C$ flips, which for reasons I'll talk about is probably optimal.

### Where

Millikan 2099, Pomona College

## Hom Quandles

### Speaker

Alissa Crans (Loyola Marymount University)

### Abstract

Analogous to the case for groups, the collection of quandle homomorphisms, Hom(Q, X), has no natural quandle structure. However, if X is an abelian quandle, then the hom set does become a quandle with the obvious pointwise operation. We will consider examples and investigate properties of this hom quandle.

### Where

Millikan 2099, Pomona College

## Solving quadratic equations over Q-bar

### Speaker

Lenny Fukshansky (CMC)

### Abstract

Given a quadratic equation in N variables over a fixed number field K, there exists an algorithm to determine whether it has a non-trivial solution over K and to find such a solution. This matter becomes considerably more complicated for a system of quadratic equations: existence of a general such algorithm would contradict Matijasevich's negative answer to Hilbert's 10th problem. The problem however is more tractable if we allow searching in extensions of K. I will discuss an approach to this problem, which involves height functions, the common tools of Diophantine geometry. Our investigation extends previous results on small-height zeros of quadratic forms, including Cassels' theorem and its various generalizations and contributes to the literature of so-called "absolute" Diophantine results with respect to height.

### Where

Millikan 2099, Pomona College

### Misc. Information

Organizational meeting in the same room preceding the talk at 12:00 noon.

## Applications of combinatorial representation theory to machine learning

### Speaker

Lily Silverstein (CGU)

### Abstract

Some machine learning problems are naturally modeled by probability distributions (or other data) defined over a finite group. In this case the generalized Fourier transform, based on irreducible group representations, is a useful tool for designing efficient algorithms. In the case of the symmetric group, we can use combinatorial objects like Young diagrams to define an analogue to bandlimiting. Finally, I will talk about how certain probabilistic inferences can be performed directly in the Fourier domain, by considering the combinatorial decomposition of tensor products of representations. This talk is expository and based mainly on work done by Jonathan Huang and Risi Kondor.

### Where

Mudd Science Library 126, Pomona College

## Group Frames and Tests for Uniformity

### Speaker

Mike Orrison (HMC)

### Abstract

Group frames are special spanning sets of vectors in representations of finite groups that sometimes behave remarkably like orthogonal bases.  In this talk, I’ll explain why group frames are beginning to play a prominent role in the work I have been doing with Anna Bargagliotti on linear tests of uniformity for probability distributions defined on finite sets with symmetry.  Along the way, I’ll try to convince you that group frames should be a part of the next linear algebra course you teach!

### Where

Mudd Science Library 126, Pomona College

### Misc. Information

Proudly Serving Math Community at the Claremont Colleges Since 2007

format: plain


text:
# Q# factorization Method Program

Is there a method to calculate prime factorization in Q# of Visual Studio of Microsoft Quantum Development Kit?

In the Python language, it is the sympy.factorint method.

Or else is there any Q# sample program?

It is like Shor's algorithm calling Q# from C#.


text:
# Weakly non-completable partial Latin square

An empty cell in a partial Latin square (pLs) is said to be forced if it has a unique admissible entry (compatible with the definition of a Latin square). Attempting to complete a given pLs, one can start by successively filling in these forced entries.

If the pLs is completable, this can lead to the following two situations:

1. The pLs can be completed solely by forcing. Such pLs are usually called strongly completable. Example: $\left(\begin{array}{ccc} 1&-&- \\-&3&- \\ -&-&- \end{array}\right)$; the middle entry in the top row is forced to be 2, and so on.
2. The pLs cannot be completed solely by forcing. At some point one has to make a case analysis (there might still exist a unique completion). Such pLs are usually called weakly completable. Example: $\left(\begin{array}{ccccc}-&-&-&-&5\\-&1&4&-&-\\3&-&5&-&-\\4&-&2&-&-\\-&-&-&2&4\end{array}\right)$; no forced entry, but nevertheless has a unique completion.

However, I'm interested in partial Latin squares that are not completable. If we attempt to fill a non-completable pLS by forcing, we run into either of the following two situations:

1. Forcing leads to a contradiction (i.e. there is a cell with no admissible entry). We might call such a pLs strongly non-completable.
Example: $\left(\begin{array}{ccc}1&-&-\\-&3&-\\-&-&1\end{array}\right)$; both the left and right entry in the middle row are forced to be 2.
2. Forcing does not lead to a contradiction; at some point we have to make a case analysis to show that the pLs is not completable. Let's call such a pLs weakly non-completable.
Example: ???

Question: All non-completable pLs I know of are strongly non-completable. Are there known examples of weakly non-completable partial Latin squares (maybe even a trivial one I'm missing)? Is there a way to systematically construct such examples? More generally, has the notion 4. been studied before?

$$\left(\begin{array}{ccccc}1&2&3&-&-\\2&3&1&-&-\\3&1&2&-&-\\-&-&-&-&-\\-&-&-&-&-\end{array}\right)$$


text:
# Euler's theorem

Euler's theorem

In number theory, Euler's theorem (also known as the Fermat-Euler theorem or Euler's totient theorem) states that if "n" is a positive integer and "a" is coprime to "n", then:$a^\left\{varphi \left(n\right)\right\} equiv 1 pmod\left\{n\right\}$where φ("n") is Euler's totient function and "... ≡ ... (mod "n")" denotes congruence modulo n.

The theorem is a generalization of Fermat's little theorem, and is further generalized by Carmichael's theorem.

The theorem may be used to easily reduce large powers modulo "n". For example, consider finding the last decimal digit of 7222, i.e. 7222 (mod 10). Note that 7 and 10 are coprime, and φ(10) = 4. So Euler's theorem yields 74 ≡ 1 (mod 10), and we get 7222 ≡ 74x55 + 2 ≡ (74)55x72 ≡ 155x72 ≡ 49 ≡ 9 (mod 10).

In general, when reducing a power of "a" modulo "n" (where "a" and "n" are coprime), one needs to work modulo φ("n") in the exponent of "a": :if "x" ≡ "y" (mod φ("n")), then "a""x" ≡ "a""y" (mod "n").

Proofs

Leonhard Euler published a proof in 1736. Using modern terminology, one may prove the theorem as follows: the numbers "a" which are relatively prime to "n" form a group under multiplication mod "n", the group of units of the ring Z/"n"Z. This group has φ("n") elements, and the statement of Euler's theorem follows then from Lagrange's theorem.

Another direct proof: if "a" is coprime to "n", then multiplication by "a"permutes the residue classes mod "n" that are coprime to "n"; in other words (writing "R" forthe set consisting of the φ("n") different such classes) the sets{ "x" : "x" in "R" } and{ "ax" : "x" in "R" } are equal; therefore,their products are equal. Hence, "P" ≡ "a"φ("n")"P" (mod "n") where "P" is the first of those products. Since "P" is coprime to "n", it follows that "a"φ("n") ≡ 1 (mod "n").

The Mizar project has completely formalized and automatically checked a proof of Euler's theorem in the [http://www.mizar.org/JFM/Vol10/euler_2.html EULER_2 file] .

*
* [http://planetmath.org/encyclopedia/EulersTheorem.html Euler's Theorem] at [http://planetmath.org PlanetMath]

Wikimedia Foundation. 2010.

### Look at other dictionaries:

• Euler'sches Theorem — Ausschöpfungstheorem, Adding up Theorem. Bei linear homogenen Produktionsfunktionen gilt: wobei: f´i = partielle Grenzproduktivität des Faktors i, ri = gesamte Einsatzmenge des Faktors i, Q = Output. Bei vollständiger Konkurrenz ist das… …   Lexikon der Economics

• Euler's theorem (differential geometry) — In the mathematical field of differential geometry, Euler s theorem is a result on the curvature of curves on a surface. The theorem establishes the existence of principal curvatures and associated principal directions which give the directions… …   Wikipedia

• Euler's theorem in geometry — In geometry, Euler s theorem, named after Leonhard Euler, states that the distance d between the circumcentre and incentre of a triangle can be expressed as : d^2=R (R 2r) ,where R and r denote the circumradius and inradius respectively (the… …   Wikipedia

• Euler's totient function — For other functions named after Euler, see List of topics named after Leonhard Euler. The first thousand values of φ(n) In number theory, the totient φ(n) of a positive integer n is defined to be the number of positive integers less than or equal …   Wikipedia

• Euler's formula — This article is about Euler s formula in complex analysis. For Euler s formula in algebraic topology and polyhedral combinatorics see Euler characteristic.   Part of a series of articles on The mathematical constant e …   Wikipedia

• Euler's rotation theorem — In kinematics, Euler s rotation theorem states that, in three dimensional space, any displacement of a rigid body such that a point on the rigid body remains fixed, is equivalent to a rotation about a fixed axis through that point. The theorem is …   Wikipedia

• Euler-Bernoulli beam equation — Euler Bernoulli beam theory, or just beam theory, is a simplification of the linear theory of elasticity which provides a means of calculating the load carrying and deflection characteristics of beams. It was first enunciated circa 1750, but was… …   Wikipedia

• Euler's sum of powers conjecture — Euler s conjecture is a disproved conjecture in mathematics related to Fermat s last theorem which was proposed by Leonhard Euler in 1769. It states that for all integers n and k greater than 1, if the sum of n k th powers of positive integers is …   Wikipedia

• Euler's criterion — In mathematics, Euler s criterion is used in determining in number theory whether a given integer is a quadratic residue modulo a prime. DefinitionEuler s criterion states: Let p be an odd prime and a an integer coprime to p . Then a is a… …   Wikipedia

• Euler pseudoprime — An odd composite integer n is called an Euler pseudoprime to base a , if a and n are coprime, and : a^{(n 1)/2} equiv pm 1pmod{n}(where mod refers to the modulo operation).The motivation for this definition is the fact that all prime numbers p… …   Wikipedia


text:
# LabWindows/CVI

cancel
Showing results for
Did you mean:

## How do I compile LW/CVI application using winsock2 SDK from Microsoft?

I'm using CVI 5.0.1, linking with ws2_32.lib and providing the include path to the winsock2 header files. Still getting winsock2.h errors, link errors.
Message 1 of 6
(3,997 Views)

## Re: How do I compile LW/CVI application using winsock2 SDK from Microsoft?

I don't have problem using winsock library, but I am using CVI6.0

In my C file, the order of include files is

winodws.h
winsock2.h
Includes.h (all my CVI inlclude files)

the ....CVI\sdk\include path is added under "Include Paths" option

ws2_32.lib is added to the Project from CVI\sdk\lib folder.

Hope it helps.

vishi

"fermenter" wrote in message
news:506500000008000000B9810000-1042324653000@exchange.ni.com...
> I'm using CVI 5.0.1, linking with ws2_32.lib and providing the include
> path to the winsock2 header files. Still getting winsock2.h errors,
Message 2 of 6
(3,997 Views)

## Re: How do I compile LW/CVI application using winsock2 SDK from Microsoft?

Hello,
Please make sure that you are using the Win32 SDK support that ships with LW/CVI. If you don't have it already installed, you can perform a custom installation and manually add the support for Win32 SDK. Once you have the support installed, you don't have to explicitly specify the directory in the include path. However, you have to add the library file in your CVI project in order to link your files correctly.

Mika Fukuchi
Application Engineer
National Instruments
Message 3 of 6
(3,997 Views)

## Re: How do I compile LW/CVI application using winsock2 SDK from Microsoft?

In my CVI 6.0 installation, i don't have W2_32.lib, from where i can get it?

Message 4 of 6
(3,661 Views)

## Re: How do I compile LW/CVI application using winsock2 SDK from Microsoft?

If you have the Full Development System (FDS) of CVI, then you need to rerun the installer, and when you get to the feature tree, make sure that you enable the LabWindows/CVI>>Win32 SDK feature.

If, on the hand, you have the Base version of CVI, then I'm afraid you're out of luck. These additional SDK files are one of the differences between Base and FDS.

Luis
Message 5 of 6
(3,641 Views)

## Re: How do I compile LW/CVI application using winsock2 SDK from Microsoft?

Thanks LuisG
I think i have to reinstall my CVI and try again.
Vishnu
Message 6 of 6
(3,591 Views)


text:
# Coordinate systems and affines¶

A nibabel (and nipy) image is the association of three things:

• The image data array: a 3D or 4D array of image data
• An affine array that tells you the position of the image array data in a reference space.
• image metadata (data about the data) describing the image, usually in the form of an image header.

This document describes how the affine array describes the position of the image data in a reference space. On the way we will define what we mean by reference space, and the reference spaces that Nibabel uses.

## Introducing Someone¶

We have scanned someone called “Someone”, and we have a two MRI images of their brain, a single EPI volume, and a structural scan. In general we never use the person’s name in the image filenames, but we make an exception in this case:

We can load up the EPI image to get the image data array:

>>> import nibabel as nib
>>> epi_img_data = epi_img.get_fdata()
>>> epi_img_data.shape
(53, 61, 33)

Then we have a look at slices over the first, second and third dimensions of the array.

>>> import matplotlib.pyplot as plt
>>> def show_slices(slices):
...    """ Function to display row of image slices """
...    fig, axes = plt.subplots(1, len(slices))
...    for i, slice in enumerate(slices):
...        axes[i].imshow(slice.T, cmap="gray", origin="lower")
>>>
>>> slice_0 = epi_img_data[26, :, :]
>>> slice_1 = epi_img_data[:, 30, :]
>>> slice_2 = epi_img_data[:, :, 16]
>>> show_slices([slice_0, slice_1, slice_2])
>>> plt.suptitle("Center slices for EPI image")

(png, hires.png, pdf)

We collected an anatomical image in the same session. We can load that image and look at slices in the three axes:

>>> anat_img_data = anat_img.get_fdata()
>>> anat_img_data.shape
(57, 67, 56)
>>> show_slices([anat_img_data[28, :, :],
...              anat_img_data[:, 33, :],
...              anat_img_data[:, :, 28]])
>>> plt.suptitle("Center slices for anatomical image")

As is usually the case, we had a different field of view for the anatomical scan, and so the anatomical image has a different shape, size, and orientation in the magnet.

## Voxel coordinates are coordinates in the image data array¶

As y’all know, a voxel is a pixel with volume.

In the code above, slice_0 from the EPI data is a 2D slice from a 3D image. The plot of the EPI slices displays the slices in grayscale (graded between black for the minimum value, white for the maximum). Each pixel in the slice grayscale image also represents a voxel, because this 2D image represents a slice from the 3D image with a certain thickness.

The 3D array is therefore also a voxel array. As for any array, we can select particular values by indexing. For example, we can get the value for the middle voxel in the EPI data array like this:

>>> n_i, n_j, n_k = epi_img_data.shape
>>> center_i = (n_i - 1) // 2  # // for integer division
>>> center_j = (n_j - 1) // 2
>>> center_k = (n_k - 1) // 2
>>> center_i, center_j, center_k
(26, 30, 16)
>>> center_vox_value = epi_img_data[center_i, center_j, center_k]
>>> center_vox_value
81.5492877960205...

The values (26, 30, 16) are indices into the data array epi_img_data. (26, 30, 16) is therefore a ‘voxel coordinate’ - a coordinate into the voxel array.

A coordinate is a set of numbers giving positions relative to a set of axes. In this case 26 is a position on the first array axis, where the axis is of length epi_img_data.shape[0], and therefore goes from 0 to 52 (epi_img_data.shape == (53, 61, 33)). Similarly 30 gives a position on the second axis (0 to 60) and 16 is the position on the third axis (0 to 32).

## Voxel coordinates and points in space¶

The voxel coordinate tells us almost nothing about where the data came from in terms of position in the scanner. For example, let’s say we have the voxel coordinate (26, 30, 16). Without more information we have no idea whether this voxel position is on the left or right of the brain, or came from the left or right of the scanner.

This is because the scanner allows us to collect voxel data in almost any arbitrary position and orientation within the magnet.

In the case of Someone’s EPI, we took transverse slices at a moderate angle to the floor to ceiling direction. This localizer image from the scanner console has a red box that shows the position of the slice block for someones_epi.nii.gz and a blue box for the slice block of someones_anatomy.nii.gz:

The localizer is oriented to the magnet, so that the left and right borders of the image are parallel to the floor of the scanner room, with the left border being towards the floor and the right border towards the ceiling.

You will see from the labels on the localizer that the center of the EPI voxel data block (at 26, 30, 16 in epi_img_data) is not quite at the center of magnet bore (the magnet isocenter).

We have an anatomical and an EPI scan, and later on we will surely want to be able to relate the data from someones_epi.nii.gz to someones_anatomy.nii.gz. We can’t easily do this at the moment, because we collected the anatomical image with a different field of view and orientation to the EPI image, so the voxel coordinates in the EPI image refer to different locations in the magnet to the voxel coordinates in the anatomical image.

We solve this problem by keeping track of the relationship of voxel coordinates to some reference space. In particular, the affine array stores the relationship between voxel coordinates in the image data array and coordinates in the reference space. We store the relationship of voxel coordinates from someones_epi.nii.gz and the reference space, and also the (different) relationship of voxel coordinates in someones_anatomy.nii.gz to the same reference space. Because we know the relationship of (voxel coordinates to the reference space) for both images, we can use this information to relate voxel coordinates in someones_epi.nii.gz to spatially equivalent voxel coordinates in someones_anatomy.nii.gz.

## The scanner-subject reference space¶

What does “space” mean in the phrase “reference space”? The space is defined by an ordered set of axes. For our 3D spatial world, it is a set of 3 independent axes.

We can decide what space we want to use, by choosing these axes. We need to choose the origin of the axes, their direction and their units.

To start with, we define a set of three orthogonal scanner axes.

### The scanner axes¶

• The origin of the axes is at the magnet isocenter. This is coordinate (0, 0, 0) in our reference space. All three axes pass through the isocenter.
• The units for all three axes are millimeters.
• Imagine an observer standing behind the scanner looking through the magnet bore towards the end of the scanner bed. Imagine a line traveling towards the observer through the center of the magnet bore, parallel to the bed, with the zero point at the magnet isocenter, and positive values closer to the observer. Call this line the scanner-bore axis.
• Draw a line traveling from the scanner room floor up through the magnet isocenter towards the ceiling, at right angles to the scanner bore axis. 0 is at isocenter and positive values are towards the ceiling. Call this the scanner-floor/ceiling axis.
• Draw a line at right angles to the other two lines, traveling from the observer’s left, parallel to the floor, and through the magnet isocenter to the observer’s right. 0 is at isocenter and positive values are to the right. Call this the scanner-left/right.

If we make the axes have order (scanner left-right; scanner floor-ceiling; scanner bore) then we have an ordered set of 3 axes and therefore the definition of a 3D space. Call the first axis the “X” axis, the second “Y” and the third “Z”. A coordinate of $$(x, y, z) = (10, -5, -3)$$ in this space refers to the point in space 10mm to the (fictional observer’s) right of isocenter, 5mm towards the floor from the isocenter, and 3mm towards the foot of the scanner bed. This reference space is sometimes known as “scanner XYZ”. It was the standard reference space for the predecessor to DICOM, called ACR / NEMA 2.0.

### From scanner to subject¶

If the subject is lying in the usual position for a brain scan, face up and head first in the scanner, then scanner-left/right is also the left-right axis of the subject’s head, scanner-floor/ceiling is the anterior-posterior axis of the head and scanner-bore is the inferior-superior axis of the head.

Sometimes the subject is not lying in the standard position. For example, the subject may be lying with their face pointing to the right (in terms of the scanner-left/right axis). In that case “scanner XYZ” will not tell us about the subject’s left and right, but only the scanner left and right. We might prefer to know where we are in terms of the subject’s left and right.

To deal with this problem, most reference spaces use subject- or patient- centered scanner coordinate systems. In these systems, the axes are still the scanner axes above, but the ordering and direction of the axes comes from the position of the subject. The most common subject-centered scanner coordinate system in neuroimaging is called “scanner RAS” (right, anterior, superior). Here the scanner axes are reordered and flipped so that the first axis is the scanner axis that is closest to the left to right axis of the subject, the second is the closest scanner axis to the anterior-posterior axis of the subject, and the third is the closest scanner axis to the inferior-superior axis of the subject. For example, if the subject was lying face to the right in the scanner, then the first (X) axis of the reference system would be scanner-floor/ceiling, but reversed so that positive values are towards the floor. This axis goes from left to right in the subject, with positive values to the right. The second (Y) axis would be scanner-left/right (anterior-posterior in the subject), and the Z axis would be scanner-bore (inferior-superior).

### Naming reference spaces¶

Reading names of reference spaces can be confusing because of different meanings that authors use for the same terms, such as ‘left’ and ‘right’.

We are using the term “RAS” to mean that the axes are (in terms of the subject): left to Right; posterior to Anterior; and inferior to Superior, respectively. Although it is common to call this convention “RAS”, it is not quite universal, because some use “R”, “A” and “S” in “RAS” to mean that the axes starts on the right, anterior, superior of the subject, rather than ending on the right, anterior, superior. In other words, they would use “RAS” to refer to a coordinate system we would call “LPI”. To be safe, we’ll call our interpretation of the RAS convention “RAS+”, meaning that Right, Anterior, Posterior are all positive values on these axes.

Some people also use “right” to mean the right hand side when an observer looks at the front of the scanner, from the foot the scanner bed. Unfortunately, this means that you have to read coordinate system definitions carefully if you are not familiar with a particular convention. We nibabel / nipy folks agree with most of our brain imaging friends and many of our enemies in that we always use “right” to mean the subject’s right.

## Voxel coordinates are in voxel space¶

We have not yet made this explicit, but voxel coordinates are also in a space. In this case the space is defined by the three voxel axes (first axis, second axis, third axis), where 0, 0, 0 is the center of the first voxel in the array and the units on the axes are voxels. Voxel coordinates are therefore defined in a reference space called voxel space.

## The affine matrix as a transformation between spaces¶

We have voxel coordinates (in voxel space). We want to get scanner RAS+ coordinates corresponding to the voxel coordinates. We need a coordinate transform to take us from voxel coordinates to scanner RAS+ coordinates.

In general, we have some voxel space coordinate $$(i, j, k)$$, and we want to generate the reference space coordinate $$(x, y, z)$$.

Imagine we had solved this, and we had a coordinate transform function $$f$$ that accepts a voxel coordinate and returns a coordinate in the reference space:

$(x, y, z) = f(i, j, k)$

$$f$$ accepts a coordinate in the input space and returns a coordinate in the output space. In our case the input space is voxel space and the output space is scanner RAS+.

In theory $$f$$ could be a complicated non-linear function, but in practice, we know that the scanner collects data on a regular grid. This means that the relationship between $$(i, j, k)$$ and $$(x, y, z)$$ is linear (actually affine), and can be encoded with linear (actually affine) transformations comprising translations, rotations and zooms (wikipedia linear transform, wikipedia affine transform).

Scaling (zooming) in three dimensions can be represented by a diagonal 3 by 3 matrix. Here’s how to zoom the first dimension by $$p$$, the second by $$q$$ and the third by $$r$$ units:

$\begin{split}\begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} = \begin{bmatrix} p i\\ q j\\ r k\\ \end{bmatrix} = \begin{bmatrix} p & 0 & 0 \\ 0 & q & 0 \\ 0 & 0 & r \\ \end{bmatrix} \begin{bmatrix} i\\ j\\ k\\ \end{bmatrix}\end{split}$

A rotation in three dimensions can be represented as a 3 by 3 rotation matrix (wikipedia rotation matrix). For example, here is a rotation by $$\theta$$ radians around the third array axis:

$\begin{split}\begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} = \begin{bmatrix} \cos(\theta) & -\sin(\theta) & 0 \\ \sin(\theta) & \cos(\theta) & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} i\\ j\\ k\\ \end{bmatrix}\end{split}$

This is a rotation by $$\phi$$ radians around the second array axis:

$\begin{split}\begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} = \begin{bmatrix} \cos(\phi) & 0 & \sin(\phi) \\ 0 & 1 & 0 \\ -\sin(\phi) & 0 & \cos(\phi) \\ \end{bmatrix} \begin{bmatrix} i\\ j\\ k\\ \end{bmatrix}\end{split}$

A rotation of $$\gamma$$ radians around the first array axis:

$\begin{split}\begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\gamma) & -\sin(\gamma) \\ 0 & \sin(\gamma) & \cos(\gamma) \\ \end{bmatrix} \begin{bmatrix} i\\ j\\ k\\ \end{bmatrix}\end{split}$

Zoom and rotation matrices can be combined by matrix multiplication.

Here’s a scaling of $$p, q, r$$ units followed by a rotation of $$\theta$$ radians around the third axis followed by a rotation of $$\phi$$ radians around the second axis:

$\begin{split}\begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} = \begin{bmatrix} \cos(\phi) & 0 & \sin(\phi) \\ 0 & 1 & 0 \\ -\sin(\phi) & 0 & \cos(\phi) \\ \end{bmatrix} \begin{bmatrix} \cos(\theta) & -\sin(\theta) & 0 \\ \sin(\theta) & \cos(\theta) & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} p & 0 & 0 \\ 0 & q & 0 \\ 0 & 0 & r \\ \end{bmatrix} \begin{bmatrix} i\\ j\\ k\\ \end{bmatrix}\end{split}$

This can also be written:

\begin{align}\begin{aligned}\begin{split}M = \begin{bmatrix} \cos(\phi) & 0 & \sin(\phi) \\ 0 & 1 & 0 \\ -\sin(\phi) & 0 & \cos(\phi) \\ \end{bmatrix} \begin{bmatrix} \cos(\theta) & -\sin(\theta) & 0 \\ \sin(\theta) & \cos(\theta) & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} p & 0 & 0 \\ 0 & q & 0 \\ 0 & 0 & r \\ \end{bmatrix}\end{split}\\\begin{split}\begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} = M \begin{bmatrix} i\\ j\\ k\\ \end{bmatrix}\end{split}\end{aligned}\end{align}

This might be obvious because the matrix multiplication is the result of applying each transformation in turn on the coordinates output from the previous transformation. Combining the transformations into a single matrix $$M$$ works because matrix multiplication is associative – $$ABCD = (ABC)D$$.

A translation in three dimensions can be represented as a length 3 vector to be added to the length 3 coordinate. For example, a translation of $$a$$ units on the first axis, $$b$$ on the second and $$c$$ on the third might be written as:

$\begin{split}\begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} = \begin{bmatrix} i\\ j\\ k\\ \end{bmatrix} + \begin{bmatrix} a \\ b \\ c \\ \end{bmatrix}\end{split}$

We can write our function $$f$$ as a combination of matrix multiplication by some 3 by 3 rotation / zoom matrix $$M$$ followed by addition of a 3 by 1 translation vector $$(a, b, c)$$

$\begin{split}\begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} = M \begin{bmatrix} i\\ j\\ k\\ \end{bmatrix} + \begin{bmatrix} a\\ b\\ c\\ \end{bmatrix}\end{split}$

We could record the parameters necessary for $$f$$ as the 3 by 3 matrix, $$M$$ and the 3 by 1 vector $$(a, b, c)$$.

In fact, the 4 by 4 image affine array does include exactly this information. If $$m_{i,j}$$ is the value in row $$i$$ column $$j$$ of matrix $$M$$, then the image affine matrix $$A$$ is:

$\begin{split}A = \begin{bmatrix} m_{1,1} & m_{1,2} & m_{1,3} & a \\ m_{2,1} & m_{2,2} & m_{2,3} & b \\ m_{3,1} & m_{3,2} & m_{3,3} & c \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}\end{split}$

Why the extra row of $$[0, 0, 0, 1]$$? We need this row because we have rephrased the combination of rotations / zooms and translations as a transformation in homogenous coordinates (see wikipedia homogenous coordinates). This is a trick that allows us to put the translation part into the same matrix as the rotations / zooms, so that both translations and rotations / zooms can be applied by matrix multiplication. In order to make this work, we have to add an extra 1 to our input and output coordinate vectors:

$\begin{split}\begin{bmatrix} x\\ y\\ z\\ 1\\ \end{bmatrix} = \begin{bmatrix} m_{1,1} & m_{1,2} & m_{1,3} & a \\ m_{2,1} & m_{2,2} & m_{2,3} & b \\ m_{3,1} & m_{3,2} & m_{3,3} & c \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} i\\ j\\ k\\ 1\\ \end{bmatrix}\end{split}$

This results in the same transformation as applying $$M$$ and $$(a, b, c)$$ separately. One advantage of encoding transformations this way is that we can combine two sets of [rotations, zooms, translations] by matrix multiplication of the two corresponding affine matrices.

In practice, although it is common to combine 3D transformations using 4 by 4 affine matrices, we usually apply the transformations by breaking up the affine matrix into its component $$M$$ matrix and $$(a, b, c)$$ vector and doing:

$\begin{split}\begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} = M \begin{bmatrix} i\\ j\\ k\\ \end{bmatrix} + \begin{bmatrix} a\\ b\\ c\\ \end{bmatrix}\end{split}$

As long as the last row of the 4 by 4 is $$[0, 0, 0, 1]$$, applying the transformations in this way is mathematically the same as using the full 4 by 4 form, without the inconvenience of adding the extra 1 to our input and output vectors.

## The inverse of the affine gives the mapping from scanner to voxel¶

The affine arrays we have described so far have another pleasant property — they are usually invertible. As y’all know, the inverse of a matrix $$A$$ is the matrix $$A^{-1}$$ such that $$I = A^{-1} A$$, where $$I$$ is the identity matrix. Put another way:

\begin{align}\begin{aligned}\begin{split}\begin{bmatrix} x\\ y\\ z\\ 1\\ \end{bmatrix} = A \begin{bmatrix} i\\ j\\ k\\ 1\\ \end{bmatrix}\end{split}\\\begin{split}A^{-1}\begin{bmatrix} x\\ y\\ z\\ 1\\ \end{bmatrix} = A^{-1} A \begin{bmatrix} i\\ j\\ k\\ 1\\ \end{bmatrix}\end{split}\\\begin{split}\begin{bmatrix} i\\ j\\ k\\ 1\\ \end{bmatrix} = A^{-1} \begin{bmatrix} x\\ y\\ z\\ 1\\ \end{bmatrix}\end{split}\end{aligned}\end{align}

That means that the inverse of the affine matrix gives the transformation from scanner RAS+ coordinates to voxel coordinates in the image data.

Now imagine we have affine array $$A$$ for someones_epi.nii.gz, and affine array $$B$$ for someones_anatomy.nii.gz. $$A$$ gives the mapping from voxels in the image data array of someones_epi.nii.gz to millimeters in scanner RAS+. $$B$$ gives the mapping from voxels in image data array of someones_anatomy.nii.gz to the same scanner RAS+. Now let’s say we have a particular voxel coordinate $$(i, j, k)$$ in the data array of someones_epi.nii.gz, and we want to find the voxel in someones_anatomy.nii.gz that is in the same spatial position. Call this matching voxel coordinate $$(i', j', k')$$ . We first apply the transform from someones_epi.nii.gz voxels to scanner RAS+ ($$A$$) and then apply the transform from scanner RAS+ to voxels in someones_anatomy.nii.gz ($$B^{-1}$$):

$\begin{split}\begin{bmatrix} i'\\ j'\\ k'\\ 1\\ \end{bmatrix} = B^{-1} A \begin{bmatrix} i\\ j\\ k\\ 1\\ \end{bmatrix}\end{split}$

## The affine by example¶

We can get the affine from the nibabel image object. Here is the affine for the EPI scan:

>>> # Set numpy to print 3 decimal points and suppress small values
>>> import numpy as np
>>> np.set_printoptions(precision=3, suppress=True)
>>> # Print the affine
>>> epi_img.affine
array([[  3.   ,   0.   ,   0.   , -78.   ],
[  0.   ,   2.866,  -0.887, -76.   ],
[  0.   ,   0.887,   2.866, -64.   ],
[  0.   ,   0.   ,   0.   ,   1.   ]])

As you see, the last row is $$[0, 0, 0, 1]$$

### Applying the affine¶

To make the affine simpler to apply, we split it into $$M$$ and $$(a, b, c)$$:

>>> M = epi_img.affine[:3, :3]
>>> abc = epi_img.affine[:3, 3]

Then we can define our function $$f$$:

>>> def f(i, j, k):
...    """ Return X, Y, Z coordinates for i, j, k """
...    return M.dot([i, j, k]) + abc

The labels on the localizer image give the impression that the center voxel of someones_epi.nii.gz was a little above the magnet isocenter. Now we can check:

>>> epi_vox_center = (np.array(epi_img_data.shape) - 1) / 2.
>>> f(epi_vox_center[0], epi_vox_center[1], epi_vox_center[2])
array([ 0.   , -4.205,  8.453])

That means the center of the image field of view is at the isocenter of the magnet on the left to right axis, and is around 4.2mm posterior to the isocenter and ~8.5 mm above the isocenter.

The parameters in the affine array can therefore give the position of any voxel coordinate, relative to the scanner RAS+ reference space.

We get the same result from applying the affine directly instead of using $$M$$ and $$(a, b, c)$$ in our function. As above, we need to add a 1 to the end of the vector to apply the 4 by 4 affine matrix.

>>> epi_img.affine.dot(list(epi_vox_center) + [1])
array([ 0.   , -4.205,  8.453,  1.   ])

In fact nibabel has a function apply_affine that applies an affine to an $$(i, j, k)$$ point by splitting the affine into $$M$$ and $$abc$$ then multiplying and adding as above:

>>> from nibabel.affines import apply_affine
>>> apply_affine(epi_img.affine, epi_vox_center)
array([ 0.   , -4.205,  8.453])

Now we can apply the affine, we can use matrix inversion on the anatomical affine to map between voxels in the EPI image and voxels in the anatomical image.

>>> import numpy.linalg as npl
>>> epi_vox2anat_vox = npl.inv(anat_img.affine).dot(epi_img.affine)

What is the voxel coordinate in the anatomical corresponding to the voxel center of the EPI image?

>>> apply_affine(epi_vox2anat_vox, epi_vox_center)
array([28.364, 31.562, 36.165])

The voxel coordinate of the center voxel of the anatomical image is:

>>> anat_vox_center = (np.array(anat_img_data.shape) - 1) / 2.
>>> anat_vox_center
array([28. , 33. , 27.5])

The voxel location in the anatomical image that matches the center voxel of the EPI image is nearly exactly half way across the first axis, a voxel or two back from the anatomical voxel center on the second axis, and about 9 voxels above the anatomical voxel center. We can check the localizer image by eye to see whether this makes sense, by seeing how the red EPI field of view center relates to the blue anatomical field of view center and the blue anatomical image field of view.

### The affine as a series of transformations¶

You can think of the image affine as a combination of a series of transformations to go from voxel coordinates to mm coordinates in terms of the magnet isocenter. Here is the EPI affine broken down into a series of transformations, with the results shown on the localizer image:

We start by putting the voxel grid onto the isocenter coordinate system, so a translation of one voxel equates to a translation of one millimeter in the isocenter coordinate system. Our EPI image would then have the black bounding box in the image above. Next we scale the voxels to millimeters by scaling by the voxel size (green bounding box). We could do this with an affine:

>>> scaling_affine = np.array([[3, 0, 0, 0],
...                            [0, 3, 0, 0],
...                            [0, 0, 3, 0],
...                            [0, 0, 0, 1]])

After applying this affine, when we move one voxel in any direction, we are moving 3 millimeters in that direction:

>>> one_vox_axis_0 = [1, 0, 0]
>>> apply_affine(scaling_affine, one_vox_axis_0)
array([3, 0, 0])

Next we rotate the scaled voxels around the first axis by 0.3 radians (see rotate around first axis):

>>> cos_gamma = np.cos(0.3)
>>> sin_gamma = np.sin(0.3)
>>> rotation_affine = np.array([[1, 0, 0, 0],
...                            [0, cos_gamma, -sin_gamma, 0],
...                            [0, sin_gamma, cos_gamma, 0],
...                            [0, 0, 0, 1]])
>>> affine_so_far = rotation_affine.dot(scaling_affine)
>>> affine_so_far
array([[ 3.   ,  0.   ,  0.   ,  0.   ],
[ 0.   ,  2.866, -0.887,  0.   ],
[ 0.   ,  0.887,  2.866,  0.   ],
[ 0.   ,  0.   ,  0.   ,  1.   ]])

The EPI voxel block coordinates transformed by affine_so_far are at the position of the yellow box on the figure.

Finally we translate the 0, 0, 0 coordinate at the bottom, posterior, left corner of our array to be at its final position relative to the isocenter, which is -78, -76, -64:

>>> translation_affine = np.array([[1, 0, 0, -78],
...                                [0, 1, 0, -76],
...                                [0, 0, 1, -64],
...                                [0, 0, 0, 1]])
>>> whole_affine = translation_affine.dot(affine_so_far)
>>> whole_affine
array([[  3.   ,   0.   ,   0.   , -78.   ],
[  0.   ,   2.866,  -0.887, -76.   ],
[  0.   ,   0.887,   2.866, -64.   ],
[  0.   ,   0.   ,   0.   ,   1.   ]])

This brings the affine-transformed voxel coordinates to the red box on the figure, matching the position on the localizer.

## Other reference spaces¶

The scanner RAS+ reference space is a “real-world” space, in the sense that a coordinate in this space refers to a position in the real world, in a particular scanner in a particular room.

Imagine that we used some fancy software to register someones_epi.nii.gz to a template image, such as the Montreal Neurological Institute (MNI) template brain. The registration has moved the voxels around in complicated ways — the image has changed shape to match the template brain. We probably do not want to know how the voxel locations relate to the original scanner, but how they relate to the template brain. So, what reference space should we use?

In this case we use a space defined in terms of the template brain — the MNI reference space.

• The origin (0, 0, 0) point is defined to be the point that the anterior commissure of the MNI template brain crosses the midline (the AC point).
• Axis units are millimeters.
• The Y axis follows the midline of the MNI brain between the left and right hemispheres, going from posterior (negative) to anterior (positive), passing through the AC point. The template defines this line.
• The Z axis is at right angles to the Y axis, going from inferior (negative) to superior (positive), with the superior part of the line passing between the two hemispheres.
• The X axis is a line going from the left side of the brain (negative) to right side of the brain (positive), passing through the AC point, and at right angles to the Y and Z axes.

These axes are defined with reference to the template. The exact position of the Y axis, for example, is somewhat arbitrary, as is the definition of the origin. Left and right are left and right as defined by the template. These are the axes and the space that MNI defines for its template.

A coordinate in this reference system gives a position relative to the particular brain template. It is not a real-world space because it does not refer to any particular place but to a position relative to a template.

The axes are still left to right, posterior to anterior and inferior to superior in terms of the template subject. This is still an RAS+ space — the MNI RAS+ space.

An image aligned to this template will therefore have an affine giving the relationship between voxels in the aligned image and the MNI RAS+ space.

There are other reference spaces. For example, we might align an image to the Talairach atlas brain. This brain has a different shape and size than the MNI brain. The origin is the AC point, but the Y axis passes through the point that the posterior commissure crosses the midline (the PC point), giving a slightly different trajectory from the MNI Y axis. Like the MNI RAS+ space, the Talairach axes also run left to right, posterior to anterior and inferior superior, so this is the Talairach RAS+ space.

There are conventions other than RAS+ for the reference space. For example, DICOM files map input voxel coordinates to coordinates in scanner LPS+ space. Scanner LPS+ space uses the same scanner axes and isocenter as scanner RAS+, but the X axis goes from right to the subject’s Left, the Y axis goes from anterior to Posterior, and the Z axis goes from inferior to Superior. A positive X coordinate in this space would mean the point was to the subject’s left compared to the magnet isocenter.

## Nibabel always uses an RAS+ output space¶

Nibabel images always use RAS+ output coordinates, regardless of the preferred output coordinates of the underlying format. For example, we convert affines for DICOM images to output RAS+ coordinates instead of LPS+ coordinates. We chose this convention because it is the most popular in neuroimaging; for example, it is the standard used by NIfTI and MINC formats.

Nibabel does not enforce a particular RAS+ space. For example, NIfTI images contain codes that specify whether the affine maps to scanner or MNI or Talairach RAS+ space. For the moment, you have to consult the specifics of each format to find which RAS+ space the affine maps to.


text:
Mathematics

# Jensen Polynomials, the Riemann-zeta Function, and SYK

A new paper by Michael Griffin, Ken Ono, Larry Rolen, and Don Zagier appears to have a made some intriguing steps when it comes to the Riemann Hypothesis (RH). The paper is titled, ‘Jensen polynomials for the Riemann zeta function and other sequences’. The preprint originally appeared in arXiv [arXiv:1902.07321 [math.NT]] in February 2019. It was one of a long list of papers that I wanted to read over the summer. And with the final version now in the Proceedings of the National Academy of Sciences (PNAS), I would like to discuss a bit about the author’s work and one way in which it relates to my own research.

First, the regular reader will recall that in a past post on string mass and the Riemann-zeta function, we discussed the RH very briefly, including the late Sir Michael Atiyah’s claim to have solved it, and finally the separate idea of a stringy proof. The status of Atiyah’s claim still seems unclear, though I mentioned previously that it doesn’t look like it will hold. The idea of a stringy proof also remains a distant dream. But we may at least recall from this earlier post some basic properties of the RH.

What is very interesting about the Griffin et al paper is that it returns to a rather old approach to the RH, based on George Pólya’s research in 1927. The authors also build on the work of Johan Jensen. The connection is as follows. It was the former, Pólya, a Hungarian mathematician, who proved that, for the Riemann-zeta function $\zeta{s}$ at its point of symmetry, the RH is equivalent to the hyperbolicity of Jensen polynomials. For the inquisitive reader, as an entry I recommend this 1990 article in the Journal of Mathematical Analysis and Applications by George Csordas, Rirchard S. Varga, and Istvan Vincze titled, ‘Jensen polynomials with applications to the Riemann zeta-function’.

Pólya’s work is generally very interesting, something I have been familiarising myself with in relation to the Sachdev-Ye-Kitaev model (more on this later) and quantum gravity. When it comes to the RH, his approach was left mostly abandoned for decades. But Griffin et al formulate what is basically a new general framework, leveraging Pólya’s insights, and in the process proving a few new theorems and even proving criterion pertaining to the RH.

1. Hyperbolicity of Polynomials

I won’t discuss their paper in length, instead focusing on a particular section of the work. But as a short entry to their study, Griffin et al pick up from the work of Pólya, summarising his result about how the RH is equivalent to the hyperbolicity of all Jensen polynomials associated with a particular sequence of Taylor coefficients,

$\displaystyle (-1 + 4z^{2}) \Lambda(\frac{1}{2} + z) = \sum_{n=0}^{\infty} \frac{\gamma (n)}{n!} \cdot z^{2n} \ \ (1)$

Where ${\Lambda(s) = \pi^{-s/2} \Gamma (s/2)\zeta{s} = \Lambda (1 - s)}$, as stated in the paper. Now, if I am not mistaken, the sequence of Taylor coefficients belongs to what is called the Laguerre-Pólya class, in which case if there is some function ${f(x)}$ that belongs to this class, the function satisfies the Laguerre inequalities.

Additionally,  the Jensen polynomial can be seen in (1). Written generally, a Jensen polynomial is of the form ${g_{n}(t) := \sum_{k = 0}^{n} {n \choose k} \gamma_{k}t^{k}}$, where ${\gamma_{k}}$‘s are positive and they satisfy the Turán inequalities ${\gamma_{k}^{2} - \gamma_{k - 1} \gamma_{k + 1} \geq 0}$.

Now, given that a polynomial with real coefficients is hyperbolic if all of its zeros are real, where read in Griffin et al how the Jensen polynomial of degree ${d}$ and shift ${n}$ in the arbitrary sequence of real numbers ${\{ \alpha (0), \alpha (1), ... \}}$ is the following polynomial,

$\displaystyle J_{\alpha}^{d,n} (X) := \sum_{j = 0}^{d} {d \choose j} \alpha (n + j)X^{j} \ \ (2)$

Where ${n}$ and ${d}$ are the non-negative integers and where, I think, ${J_{\alpha}^{d,n} (X)}$ is the hyperbolicity of polynomials. Now, recall that we have our previous Taylor coefficients ${\gamma}$. From the above result, the following statement is given that the RH is equivalent to ${J_{\gamma}^{d,n}(X)}$ – the hyperbolicity of polynomials – for all non-negative integers. What is very curious, and what I would like to look into a bit more, is how this conditions holds under differentiation. In any case, as the authors point out, one can prove the RH by showing hyperbolicity for ${J_{\alpha}^{d,n} (X)}$; but proving the RH is of course notoriously difficult!

Alternatively, another path may be chosen. My understanding is that Griffin-Ono-Rolen-Zagier use shifts in ${n}$ for small ${d}$, because, from what I understand about hyperbolic polynomials, one wants to limit the hyperbolicity in the ${d}$ direction. Then the idea, should I not be corrected, is to study the asymptotic behaviour of ${\gamma(n)}$.

This is the general entry, from which the authors then go on to consider a number of theorems. I won’t go through all of the theorems. One can just as well read the paper and the proofs. What I want to do is focus particularly on Theorem 3.

2. Theorem 3

Aside from the more general considerations and potential breakthroughs with respect to the RH, one of my interests triggered in the Griffin-Ono-Rolen-Zagier paper has to do with my ongoing studies concerning Gaussian Unitary Ensembles (GUE) and Random Matrix Theory (RMT) in the context of the Sachdev-Ye-Kitaev (SYK) model (plus similar models) and quantum gravity. Moreover, RMT has become an interest in relation to chaos and complexity, not least because in SYK and similar models we consider late-time behaviour of quantum black holes in relation to theories of quantum chaos and random matrices.

But for now, one thing that is quite fascinating about Jensen polynomials for the Riemann-zeta function is the proof in Griffin et al of the GUE random matrix model prediction. That is, the derivative aspect GUE random matrix model prediction for the zeros of Jensen polynomials. One of the claims here is that the GUE and the RH are satisfied by the symmetric version of the zeta function. To quote in length,

‘To make this precise, recall that Dyson, Montgomery, and Odlyzko [9, 10, 11] conjecture that the nontrivial zeros of the Riemann zeta function are distributed like the eigenvalues of random Hermitian matrices. These eigenvalues satisfy Wigner’s Semicircular Law, as do the roots of the Hermite polynomials ${H_{d}(X)}$, when suitably normalized, as ${d \rightarrow +\infty}$ (see Chapter 3 of [12]). The roots of ${J){\gamma}^{d,0} (X)}$, as ${d \rightarrow +\infty}$ approximate the zeros of ${\Lambda (\frac{1}{2} + z)}$ (see [1] or Lemma 2.2 of [13]), and so GUE predicts that these roots also obey the Semicircular Law. Since the derivatives of ${\Lambda (\frac{1}{2} + z)}$ are also predicted to satisfy GUE, it is natural to consider the limiting behavior of ${J_{\gamma}^{d,n}(X)}$ as ${n \rightarrow +\infty}$. The work here proves that these derivative aspect limits are the Hermite polynomials ${H_{d}(X)}$, which, as mentioned above, satisfy GUE in degree aspect.’

I think Theorem 3 raises some very interesting, albeit searching questions. I also think it possibly raises or inspires (even if naively) some course of thought about the connection of insights being made in SYK and SYK-like models, RMT more generally, and even studies of the zeros of the Riemann-zeta function in relation to quantum black holes. In my own mind, I also immediately think of the Hilbert-Polya hypothesis and the Jensen polynomials in this context, as well as ideas pertaining to the eigenvalues of Hamiltonians in different random matrix models of quantum chaos. There is connection and certainly also an interesting analogy here. To what degree? It is not entirely clear, from my current vantage. There are also some differences that need to be considered in all of these areas. But it may not be naive to ask, in relation to some developing inclinations in SYK and other tensor models, about how GUE random matrices and local Riemann zeros are or may be connected.

Perhaps I should save such considerations for a separate article.

Standard


text:
## Simple questions about Lp spaces

1. The problem statement, all variables and given/known data

My question is just on the definition of L.

Is L=Lp where p=∞, i.e.,

is a measurable function in L if ∫Alf(x)l<∞?

2. Relevant equations

*L: The space of all bounded measurable functions on [0,1] (bounded except for possibly on a set of measure zero)

*A measurable function is said to belong to Lp if ∫Alf(x)lp<∞.

3. The attempt at a solution

Looks like it would be true based on the definition of Lp but I am really not sure since Royden only gives the one definition on L.
 PhysOrg.com science news on PhysOrg.com >> City-life changes blackbird personalities, study shows>> Origins of 'The Hoff' crab revealed (w/ Video)>> Older males make better fathers: Mature male beetles work harder, care less about female infidelity
 Sorry, this should have have been posted in the calculus and beyond section.
 The infinity norm just take the maximum value of f. For general p norms, as p increases, it puts more weight on the larger terms. If you want an intuitive clue to the idea, notice that i you are measuring the hypotenuse based off of two values (the two sides), then if one of the sides is much larger than the other, than the hypotenuse comes pretty close to the length of the longer side. The infinity norm goes further, and simply gives you the length of the longer side, completely ignoring other lengths. So does it make sense why L^infinity is all the bounded functions etc? The problem with your expression with the integral of magnitude to the infinity, is when that is p, we really mean to take the pth root of the integral. So in the infinity norm, you would take the infinitieth root. So instead (they don't do that), they take the limit in p of the pth root, and prove that it is the supremum of |f|, but I might be forgetting a detail about sets of nonzero measure or something.

## Simple questions about Lp spaces

So you need to take the limit of the p-norm as p goes to infty and show that this gives the "essential" supremum (i.e. the L infty norm). You prove this using the squeeze lemma.

Lets assume the measure of the entire space is finite for simplicity. I think the same arguments could be adapted, but suppose the measure of the whole space is M. According to Jensen's inequality applied to the convex function phi(x)=x^p, and using the original measure divided by M (giving a probability measure), we can prove that
$\|f\|_p \leq M^{1/p}\|f\|_\infty$
As long as M is not zero, then we find that the limsup of the p-norms is less than or equal to the L infinity norm.

For the other direction, let epsilon be positive. We know that there is a measurable subset E of positive measure on which the value of f is at least equal to $\|f\|_\infty - \epsilon$.

The Lp norm of f is bounded below by its Lp norm calculated over E.
$\|f\|_p\geq \left( \int_E |f|^p d\mu\right)^{1/p}\geq \left( (\|f\|_\infty-\epsilon)^p\mu(E)\right)^{1/p}=(\|f\|_\infty-\epsilon)(\mu(E))^{1/p}$

As long as the measure of E is positive, then we may conclude that the limsup of the Lp norm as p gets large is the L infinity norm.

 Similar discussions for: Simple questions about Lp spaces Thread Forum Replies Calculus & Beyond Homework 6 Calculus & Beyond Homework 7 Calculus 1 Calculus & Beyond Homework 3 Linear & Abstract Algebra 8


text:
# Harmonic oscillator's energy levels

1. May 18, 2013

### Rorshach

1. The problem statement, all variables and given/known data
Okay, this one confuses me a bit:

A particle is in a one-dimensional harmonic oscillator. At time t = 0 is given by its wave function

ψ(x)=Nx3exp(-mωx2/2hbarred)

a) At this point you measure the particle's energy. What measurement values are available? Also determine the corresponding probabilities!
b) After the power supply that gave outcome E = 3hbarredω/2 measure the particle's position immediately. What is the probability of finding the particle in the classically forbidden region? (The classically forbidden region is defined by the condition that V (x)> = E_total)

2. Relevant equations

En=(n+1/2)hbarredω

3. The attempt at a solution
I don't really have an idea how to come up with energy levels that are available, or how to calculate the probabilities of those levels.

2. May 18, 2013

### BruceW

"one-dimensional harmonic oscillator" should imply the form of the Hamiltonian. You can almost guess what it should look like. If it is not in your notes, then you could google it if you are having trouble.

3. May 19, 2013

### Rorshach

I got notice on violation of rules of the forum by posting a problem that I did not try to solve on my own- while it is not entirely true, I understand it could be perceived as such. I am sincerely sorry, won't happen again. I tried to come with idea how to find out how many possible energy levels are available for my function, but only thing that comes to my mind is this basic formula for energy E=(n+1/2)hbarredω. Hamiltonian for my function is of course
h=-hbarred2/2m(d/dx)+1/2mω2x2 , but how to connect it to my search for available energy levels?

4. May 19, 2013

### fluidistic

I don't really understand why would the Hamiltonian help to solve for a).
At first glance the psi function they give you is a supperposition of only 2 eigenstates (can you find them?) that have a definite energy. I guess that these energies are the one you could possibly measure.

As for the energy level formula $E_n =\hbar \omega (n+1/2)$, it comes from solving the Schrödinger's equation and assuming that the solution behaves as $f(x)v(x)$ where v(x) is a decaying exponential that solves the Schrödinger's equation when x tends to infinity and the function f(x) are the Hermite polynomials that you can determine when solving the original Schrödinger's equation using a power series "guess function" for f(x). You get that the solution doesn't diverge if and only if the power series is truncated for some n. This condition gives precisely that $E_n =\hbar \omega (n+1/2)$.

5. May 20, 2013

### Rorshach

I don't know how to find those two eigenstates, how can You tell from first glance that there will be two levels?

6. May 20, 2013

### fluidistic

Hint: Look at the first few Hermite polynomials and at the psi(x) they give you.

7. May 20, 2013

### Rorshach

So it is connected to power to which x is taken right? if x was to the power of 10 there would be 6 levels?

8. May 20, 2013

### fluidistic

If you only had x to the power of 10, I think so. But if you had an x^10 as well as say an x^9 in your psi(x), then no, not only necessarily 6 levels.
The idea is that the $\psi _n (x)$ form a basis for any function $\Psi (x,t)$ for a fixed t. In your case t is fixed at $t=0$.
So that you can write your $\Psi (x,0)$ as a linear combination of the $\psi _n(x)$. One writes $\Psi (x,0)=\Psi(x)=\sum _{n=0}^\infty c_n \psi _n (x)$.
We're lucky enough that the from the $\Psi (x)$ they give you, it's not hard to determine the coefficients $c_n$ and the $\psi _n (x)$ involved.
I'll also tell you a "trick". The associated probabilities with each energy levels are simply the values of $|c_n|^2$. But in your exercise, I think you're meant to find them by integration (I've never done it, so I'll have to do it myself too!). Though it's good to know 2 methods so you can compare the results afterwards.

9. May 20, 2013

### Rorshach

okay, but how do I determine the cn and the ψn(x)? It sounds similar to distributing the function into the Fourier series, but I have never had a clue how to perform Fourier series distribution on a function... if I am right and this is the thing to do, could you show me on example how to perform Fourier series distribution?

10. May 20, 2013

### fluidistic

I'm not sure it's a "Fourier expansion" but the idea is pretty similar. We must exploit the orthogonality property of the eigenstates; or the fact that they form a basis for any psi(x).
We have $\Psi(x)=\sum _{n=0}^\infty c_n \psi _n (x)$. Multiply both sides by $\psi ^* _m (x)$ and integrate from - infinity to positive infinity. After a bit of algebra you'll reach that $c_m = \int _{-\infty}^\infty \psi ^* _ m (x) \Psi (x) dx$.

11. May 20, 2013

### Rorshach

I'm sorry, but I still don't know how to get ψn. It seems to me that without it I cannot do anything.

12. May 20, 2013

### fluidistic

You get them when you solve the -normalized- Schrödinger's equation ($\psi '' + (2\varepsilon - y^2)\psi =0$ where $y=x(m \omega / \hbar )^{1/2}$ and $\varepsilon = E/(\hbar \omega )$) and assume that the solutions you're looking for are of the form $\psi (y)=v(y)\exp (-y^2/2)$. And in order to calculate v(y) you propose a power series of the form $v(y)=\sum _{k=0}^\infty C_k y^k$.
In order for $\psi (y)$ to converge you'll get that the power series for v(y) must be truncated at some value for k=n.
You'll also get the recurrence relation $C_{k+2}=C_k \frac{2k+1-2\varepsilon }{(k+2)(k+1)}$. So for any value of n, you get that v(y) is a Hermite polynomial.
And thus that the $\psi (y)$ turns out to be of the form $\psi _n (y)= H_ n (y)e^{-y^2/2}$.
After replacing the variables back and normalizing, one reaches that $\psi _n (x)=\left ( \frac{m\omega }{\pi \hbar 2^{2n}n!^2} \right ) ^{1/4} \exp \left ( - \frac{m\omega x^2}{2 \hbar} \right ) H_n \left [ \left ( \frac{m\omega }{\hbar} \right ) ^{1/2 }x \right ]$.

So yes, in your exercise I don't know how you could solve it without knowing the eigenfunctions $\psi _n (x)$. Either you look at them in a book/Internet, either you find them yourself (which is a pain but maybe necessary).

13. May 22, 2013

### Rorshach

Ok, but how do I know into which polynomials I am supposed to distribute my wave function? I only know that the hermite polynomial of the third power of x is equal to 8x3-12x. Am I supposed to look for some coefficients that will convert that polynomial into my wave function?

14. May 23, 2013

### vela

Staff Emeritus
What are the first few eigenfunctions for the one-dimensional harmonic oscillator?

15. May 23, 2013

### Rorshach

Those are the functions that follow the formula $$ψ_n=C_nexp(\frac{-\xi^2}{2})H_n(\xi)$$, right?

16. May 23, 2013

### vela

Staff Emeritus
Right, assuming you have the correct expressions for the normalization constant and $\xi$. I'm going to use $\phi_n$ to denote the eigenfunctions, so you don't confuse it with the given wavefunction $\psi$.

You want to find the appropriate coefficients $a_0, a_1, a_2, \dots$ such that
$$\psi(x) = a_0 \phi_0(x) + a_1 \phi_1(x) + a_2 \phi_2(x) + a_3\phi_3(x) + \cdots$$ Substitute in for $\psi$, $\phi_0$, etc., and see if you can figure out how to solve for the coefficients. You should easily be able to see why $a_n=0$ for $n\ge 4$.

17. May 24, 2013

### Rorshach

$a_n=0$, because the highest power in the formula is 3, right? I got to the point where I am supposed to match the coefficients, but I am somewhat confused by hermitian $H_n(\xi)=\sum_{n=0}^\infty c_n\xi^n$, namely: I know that hermitian with the power of $\xi$ being equal to 3 is $H_3(\xi)=8\xi^3-12\xi$ am I supposed to fix the constants for the elements of the hermitian sum so it would match the formula for the hermitian of third power $exp(-\frac{1}{2}\xi^2)C_3(8\xi^3-12\xi)$, or I am supposed to fix the constants so it would only leave $\xi^3$?

Last edited: May 24, 2013
18. May 24, 2013

### vela

Staff Emeritus
I don't understand what you're asking. Write it out explicitly in math.

19. May 24, 2013

### Rorshach

ok, so my wavefunction is $\psi(x)=Nx^3exp(-\frac{m\omega x^2}{2\hbar})$. Now we substitute certain values: $\omega=\sqrt{\frac{k}{m}}, \alpha=\sqrt{\frac{m\omega}{\hbar}}, \lambda=\frac{2E}{\hbar\omega}, \xi=x\alpha$. Then we can write the formula for the eigenfunctions: $\psi(\xi)=C_nH(\xi)exp(-\frac{\xi^2}{2})$. We can write that $H_n(\xi)\sum_{n=0}^\infty\ c_n\xi^n$, right? Now I know that hermit polynomial $H(\xi)$for $\xi^3$ is equal to $8\xi^3-12\xi$. And now I don't know what coefficients I am supposed to be looking for: the ones that will give me just $\xi^3$, or the ones that will give me $8\xi^3-12\xi$.

Last edited: May 24, 2013
20. May 24, 2013

### vela

Staff Emeritus


text:
# Acclerating plane

this is another post, but i just want to post it in this room and see if ANYONE knows how to solve this:
I've tried everything, now i'm getting depererate

An inclined plane that makes an angle of 28° to the horizontal is mounted on wheels. A small block of mass m = 1.2 kg rests on the plane, held there by a coefficient of static friction µ = 0.73.

The plane is accelerating to the right. What is the minimum acceleration in order that the block slides down the plane?

Any help would be great !!!!!

## Answers and Replies

Related Introductory Physics Homework Help News on Phys.org
Pyrrhus
Homework Helper
Apply Newton's 2nd Law and solve for a.

The 0 degrees with the horizontal system of coordinates will be the easier one.

i tried literally everything, would anyone please help me solve it??? lol...

i came up with the following, after lots of equation, i used it to solve for a, didn't work either

g sin (theta) - mu (a sin (theta) + g cos (theta) ) = a cos (theta)

Pyrrhus
Homework Helper
Well, with the coordinate system set as i said above.

$$\sum F_{x} = ma = - n \sin \theta + F_{f} \cos \theta$$

$$\sum F_{y} = 0 = F_{friction} \sin \theta + n \cos \theta - mg$$

#### Attachments

• 15.2 KB Views: 320
Last edited:


text:
How to delay visible light

I want to delay visible light (~450nm-600nm) by 10ns. One way would just be to have it travel about 10ft since the speed of light is about 1ft/ns. Could I reduce that length by sending it through some high index of refraction material? Dispersion is a bit of a concern in that case I guess. Any recommendations for what that material should be?

This may not be a complete solution, but you can use mirrors! It's not hard to design a mirror-based optical system which will have the property you are looking for. This reminds me of an interesting Project-Euler question.

Or a less fancy but more practical way would be to use two simple planar mirrors, carefully aligned.

• I agree with using mirrors. A more practical system is to use just two ordinary planar mirrors, that are parallel (or almost-parallel) to each other. The light can bounce back and forth 10 times no problem, and it's easy to align and set up. – Steve Byrnes Jul 4 '13 at 19:11
• @SteveB Indeed, I should add that to my answer! – Ali Jul 4 '13 at 19:13

I don't know if the following would suit you, but you can run light through a coiled optical fiber.

• @Jim: This is not obvious to me. Care to explain? – akhmeteli Jul 5 '13 at 1:01
• This is in fact exactly analogous to the way nuclear and particle physics signals are delayed for in some cases and was exactly the answer that occurred to be for signals that are point-like in origin (i.e. not images) or will tolerate being pixelated at the input head of a bundle of fibers. The speed of light in typical fibers is around $2/3\,c$. – dmckee --- ex-moderator kitten Jul 5 '13 at 3:25

Wikipedia's list of refractive indices doesn't list any above 4, and I doubt you want to pass your light through 2.5 feet of something. You're probably better off just having the light bounce back and forth between some mirrors or passing it through a coiled optical fiber (like Ali and akhmeteli suggested.


text:
SG++
optimization.cpp

On this page, we look at an example application of the sgpp::optimization module.

Versions of the example are given in all languages currently supported by SG++: C++, Python, Java, and MATLAB.

For instructions on how to compile and run the example, please see Installation and Usage.

The example interpolates a bivariate test function like the tutorial.cpp (Start Here) example. However, we use B-splines here instead to obtain a smoother interpolant. The resulting sparse grid function is then minimized with the method of steepest descent. For comparison, we also minimize the objective function with Nelder-Mead's method.

First, we include all the necessary headers, including those of the sgpp::base and sgpp::optimization module.

#include <sgpp_base.hpp>
#include <algorithm>
#include <iostream>
#include <iterator>

The function $$f\colon [0, 1]^d \to \mathbb{R}$$ to be minimized is called objective function and has to derive from sgpp::optimization::ScalarFunction. In the constructor, we give the dimensionality of the domain (in this case $$d = 2$$). The eval method evaluates the objective function and returns the function value $$f(\vec{x})$$ for a given point $$\vec{x} \in [0, 1]^d$$. The clone method returns a std::unique_ptr to a clone of the object and is used for parallelization (in case eval is not thread-safe).

class ExampleFunction : public sgpp::optimization::ScalarFunction {
public:
ExampleFunction() : sgpp::optimization::ScalarFunction(2) {}
double eval(const sgpp::base::DataVector& x) {
// minimum is f(x) = -2 for x[0] = 3*pi/16, x[1] = 3*pi/14
return std::sin(8.0 * x[0]) + std::sin(7.0 * x[1]);
}
virtual void clone(std::unique_ptr<sgpp::optimization::ScalarFunction>& clone) const {
clone = std::unique_ptr<sgpp::optimization::ScalarFunction>(new ExampleFunction(*this));
}
};

Now, we can start with the main function.

void printLine() {
std::cout << "----------------------------------------"
"----------------------------------------\n";
}
int main(int argc, const char* argv[]) {
(void)argc;
(void)argv;
std::cout << "sgpp::optimization example program started.\n\n";
// increase verbosity of the output

Here, we define some parameters: objective function, dimensionality, B-spline degree, maximal number of grid points, and adaptivity.

// objective function
ExampleFunction f;
// dimension of domain
const size_t d = f.getNumberOfParameters();
// B-spline degree
const size_t p = 3;
// maximal number of grid points
const size_t N = 30;
const double gamma = 0.95;

First, we define a grid with modified B-spline basis functions and an iterative grid generator, which can generate the grid adaptively.

With the iterative grid generator, we generate adaptively a sparse grid.

printLine();
std::cout << "Generating grid...\n\n";
if (!gridGen.generate()) {
std::cout << "Grid generation failed, exiting.\n";
return 1;
}

Then, we hierarchize the function values to get hierarchical B-spline coefficients of the B-spline sparse grid interpolant $$\tilde{f}\colon [0, 1]^d \to \mathbb{R}$$.

printLine();
std::cout << "Hierarchizing...\n\n";
sgpp::base::DataVector functionValues(gridGen.getFunctionValues());
sgpp::base::DataVector coeffs(functionValues.getSize());
// solve linear system
if (!sleSolver.solve(hierSLE, functionValues, coeffs)) {
std::cout << "Solving failed, exiting.\n";
return 1;
}

We define the interpolant $$\tilde{f}$$ and its gradient $$\nabla\tilde{f}$$ for use with the gradient method (steepest descent). Of course, one can also use other optimization algorithms from sgpp::optimization::optimizer.

printLine();
std::cout << "Optimizing smooth interpolant...\n\n";
double fX0;
double ftX0;

The gradient method needs a starting point. We use a point of our adaptively generated sparse grid as starting point. More specifically, we use the point with the smallest (most promising) function value and save it in x0.

{
sgpp::base::GridStorage& gridStorage = grid.getStorage();
// index of grid point with minimal function value
size_t x0Index =
std::distance(functionValues.getPointer(),
std::min_element(functionValues.getPointer(),
functionValues.getPointer() + functionValues.getSize()));
x0 = gridStorage.getCoordinates(gridStorage[x0Index]);
fX0 = functionValues[x0Index];
ftX0 = ft.eval(x0);
}
std::cout << "x0 = " << x0.toString() << "\n";
std::cout << "f(x0) = " << fX0 << ", ft(x0) = " << ftX0 << "\n\n";

We apply the gradient method and print the results.

const double fXOpt = f.eval(xOpt);
std::cout << "\nxOpt = " << xOpt.toString() << "\n";
std::cout << "f(xOpt) = " << fXOpt << ", ft(xOpt) = " << ftXOpt << "\n\n";

For comparison, we apply the classical gradient-free Nelder-Mead method directly to the objective function $$f$$.

printLine();
std::cout << "Optimizing objective function (for comparison)...\n\n";
const double ftXOptNM = ft.eval(xOptNM);
std::cout << "\nnxOptNM = " << xOptNM.toString() << "\n";
std::cout << "f(xOptNM) = " << fXOptNM << ", ft(xOptNM) = " << ftXOptNM << "\n\n";
printLine();
std::cout << "\nsgpp::optimization example program terminated.\n";
return 0;
}

The example program outputs the following results:

sgpp::optimization example program started.

--------------------------------------------------------------------------------
Generating grid...

100.0% (N = 29, k = 3)
Done in 3ms.
--------------------------------------------------------------------------------
Hierarchizing...

Solving linear system (automatic method)...
estimated nnz ratio: 59.8%
constructing matrix (100.0%)
nnz ratio: 58.0%
Done in 0ms.
Done in 1ms.
--------------------------------------------------------------------------------
Optimizing smooth interpolant...

x0 = [0.625, 0.75]
f(x0) = -1.81786, ft(x0) = -1.81786

9 steps, f(x) = -2.000780
Done in 1ms.

xOpt = [0.589526, 0.673268]
f(xOpt) = -1.99999, ft(xOpt) = -2.00078

--------------------------------------------------------------------------------
Optimizing objective function (for comparison)...

We see that both the gradient-based optimization of the smooth sparse grid interpolant and the gradient-free optimization of the objective function find reasonable approximations of the minimum, which lies at $$(3\pi/16, 3\pi/14) \approx (0.58904862, 0.67319843)$$.


text:
# Two Algorithms to Find Primes in Patterns

@article{Sorenson2020TwoAT,
title={Two Algorithms to Find Primes in Patterns},
author={J. Sorenson and Jonathan Webster},
journal={Math. Comput.},
year={2020},
volume={89},
pages={1953-1968}
}
• Published 23 July 2018
• Mathematics, Computer Science
• Math. Comput.
Let $k\ge 1$ be an integer, and let $P= (f_1(x), \ldots, f_k(x) )$ be $k$ admissible linear polynomials over the integers, or \textit{the pattern}. We present two algorithms that find all integers $x$ where $\max{ \{f_i(x) \} } \le n$ and all the $f_i(x)$ are prime. Our first algorithm takes at most $O_P(n/(\log\log n)^k)$ arithmetic operations using $O(k\sqrt{n})$ space. Our second algorithm takes slightly more time, $O_P(n/(\log \log n)^{k-1})$ arithmetic operations, but uses only $n^{1/c… 1 Citations An Algorithm and Estimates for the Erdős-Selfridge Function (work in progress) • Mathematics Open Book Series • 2020 A new algorithm to compute the value of g(k) is presented, and computational evidence is provided to support the claim that$\hat{g}(K)$estimates reasonably well in practice, and it is proved that for large$x, $G(x,k)$ is asymptotic to $x/\hat{ g)$.

## References

SHOWING 1-10 OF 42 REFERENCES
• Mathematics
ANTS
• 1998
The behavior of ρ *(x), in particular, the point at which ρ (x) first exceeds π (x), and its asymptotic growth is examined.
Strong pseudoprimes to twelve prime bases
• Mathematics
Math. Comput.
• 2017
An algorithm to find all integers $n\le B$ that are strong pseudoprimes to the first $m$ prime bases is presented; with a reasonable heuristic assumption it is shown that it takes at most $B^{2/3+o(1)}$ time.
The Pseudosquares Prime Sieve
The pseudosquares prime sieve is presented, which finds all primes up to n in sublinear time using very little space and the primes generated by the algorithm are proven prime unconditionally.
Dissecting a Sieve to Cut Its Need for Space
• W. Galway
• Computer Science, Mathematics
ANTS
• 2000
A “dissected” sieving algorithm which enumerates primes in the interval [x 1, x 2], using $$O(x_{2}^{1/3})$$ bits of memory and using arithmetic operations on numbers of $$O(\rm ln \it x_{2}$$ bits.
Integers free of small prime factors in arithmetic progressions *
• T. Xuan
• Mathematics
Nagoya Mathematical Journal
• 2000
For real x ≥ y ≥ 2 and positive integers a, q, let Φ(x, y; a, q) denote the number of positive integers ≤ x, free of prime factors ≤ y and satisfying n ≡ a (mod q). By the fundamental lemma of sieve,
An improved sieve of Eratosthenes
The sieve of Eratosthenes will be able to use it to factor integers, and not just to produce lists of consecutive primes, and also has close ties to Voronoi's work on the Dirichlet divisor problem.
Statistical evidence for small generating sets
• Computer Science, Mathematics
• 1993
This work gives additional evidence, independent of the ERH, that primality testing can be done in deterministic polynomial time; if the bound on G(n) is correct, there is a deterministic primality test using O(log n)2 multiplications modulo n.
A heuristic asymptotic formula concerning the distribution of prime numbers
• Mathematics
• 1962
Suppose fi, f2, -*, fk are polynomials in one variable with all coefficients integral and leading coefficients positive, their degrees being hi, h2, **. , hk respectively. Suppose each of these
Prime numbers and computer methods for factorization
1. The Number of Primes Below a Given Limit.- What Is a Prime Number?.- The Fundamental Theorem of Arithmetic.- Which Numbers Are Primes? The Sieve of Eratosthenes.- General Remarks Concerning
Doubly focused enumeration of locally square polynomial values
Let f be a nonconstant squarefree polynomial. Which of the values f(c + 1), f(c + 2), . . . , f(c + H) are locally square at all small primes? This paper presents an algorithm that answers this


text:
3160
Q:

# Find the missing number in the series?4, 18, ?, 100, 180, 294, 448

 A) 48 B) 50 C) 58 D) 60

Answer:   A) 48

Explanation:

Take a Look at this series :

$23$ - $22$ =  8 - 4 = 4

$33$ - $32$ = 27 - 9 = 18

$43$ - $42$ = 64 - 16 = 48

$53$ - $52$ = 125 - 25 = 100

$63$ - $62$ = 216 - 36 = 180

$73$ - $72$ = 343 - 49 = 294

$83$ - $82$ = 512 - 64 = 448

Q:

Select the set in which the numbers are related in the same way as are the numbers of the following set.
(5, 7, 11)

 A) (15, 17, 23) B) (14, 16, 18) C) (23, 29, 31) D) (17, 19, 21)

Answer & Explanation Answer: C) (23, 29, 31)

Explanation:
Report Error

0 87
Q:

Which number will replace the question mark (?) in the following series?

20, 10, 10, 15, 30, ?

 A) 90 B) 75 C) 60 D) 45

Answer & Explanation Answer: B) 75

Explanation:
Report Error

0 248
Q:

Which of the following sequence of signs will correctly solve the given equation by replacing the question marks ?

1105 ? 65 ? 835 ? 25 ? 5 = 727

 A) ÷, +, −, × B) +, ÷, −, × C) +, −, ×, ÷ D) ×, −, +, ÷

Answer & Explanation Answer: A) ÷, +, −, ×

Explanation:
Report Error

4 526
Q:

Which number will replace the question mark (?) in the following series?

25, 42, 67, 101, 145, ?

 A) 200 B) 202 C) 199 D) 201

Answer & Explanation Answer: A) 200

Explanation:
Report Error

8 871
Q:

Select the set in which the numbers are related in the same way as are the numbers of the following set.
(8, 112, 28)

 A) (14, 364, 52) B) (12, 216, 34) C) (16, 282, 36) D) (18, 216, 22)

Answer & Explanation Answer: A) (14, 364, 52)

Explanation:
Report Error

4 638
Q:

Which number will replace the question mark (?) in the following series?
77, 92, 125, 178, 253, ?

 A) 350 B) 351 C) 352 D) 353

Answer & Explanation Answer: C) 352

Explanation:
Report Error

6 2791
Q:

Select the set in which the numbers are related in the same way as are the numbers of the following set.
(64, 7, 27)

 A) (127, 15, 161) B) (216, 15, 64) C) (343, 18, 125) D) (216, 11, 125)

Answer & Explanation Answer: D) (216, 11, 125)

Explanation:
Report Error

5 377
Q:

Which number will replace the question mark (?) in the following series ?

12, 22, 42, 72, ?

 A) 102 B) 112 C) 92 D) 122

Answer & Explanation Answer: B) 112

Explanation:
Report Error

3 1104


text:
Verifying Genetic Algorithm on Simple Functions¶

• Goal: find the minimum of $f(\boldsymbol{x}) = \boldsymbol{x}_1^2 + \boldsymbol{x}_2^2 + \boldsymbol{x}_3^2 + \boldsymbol{x}_4^2 + \boldsymbol{x}_5^2$

More hints on Jupyter notebook/programming norms¶

• Typesetting math equations: \$f(x) = x_1^2 + x_2^2 + x_3^2 + x_4^2 + x_5^2\\$
• Rules to name variables/functions: camel case, snake case, PEP 8

Naming conventions according to PEP 8:

Type Naming Convention Examples
Function Use a lowercase word or words. Separate words by underscores to improve readability. function, my_function
Variable Use a lowercase single letter, word, or words. Separate words with underscores to improve readability. x, var, my_variable
Class Start each word with a capital letter. Do not separate words with underscores. This style is called camel case. Model, MyClass
Method Use a lowercase word or words. Separate words with underscores to improve readability. class_method, method
Constant Use an uppercase single letter, word, or words. Separate words with underscores to improve readability. CONSTANT, MY_CONSTANT, MY_LONG_CONSTANT
Module Use a short, lowercase word or words. Separate words with underscores to improve readability. module.py, my_module.py
Package Use a short, lowercase word or words. Do not separate words with underscores. package, mypackage
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from itertools import product

# make sure the result is reproducible
np.random.seed(0x1a2b3c4d)

In [2]:
INIT_SIZE = 50
POP_SIZE = 100
CROSS_RATE = 0.95
MUTATION_RATE = 0.005
N_CROSS = 500
N_GEN = 200
N_DIM = 5

Generating initial population¶

An important trick of Python programming is to avoid using loops (because they are slow!). Check out how I generate the initial population without explicit loops.

• We want to avoid negative numbers as they are encoded by two's complement, but we also want to include negative numbers in our genetic algorithm. What do we do? We limit the range of a variable $x$ between 0 and 255, and then when we put $x - 127$ into the equation instead. In this way, we are able to generate both positive and negative integers smoothly.
• Notice how I specify the data type with dtype=np.uint8. This will be useful afterwards.
In [3]:
init_pop = np.random.randint(0, 100, size=(INIT_SIZE, N_DIM), dtype=np.uint8)
print('shape:', init_pop.shape)
print('excerpt:\n', init_pop[:10, :])

shape: (50, 5)
excerpt:
[[86 52 51 37 62]
[43  0 76 34 77]
[73 41 84 82  8]
[15 80 14 87 62]
[86  2 17 95  3]
[19 39 44 67 69]
[ 1  8  4 29 65]
[ 0 16 43 91 57]
[30 20 10  6 25]
[55 48 25 72 23]]

Generating chromosome¶

I am using np.uint8 because it corresponds to an unsigned 8-bit integer. We can use np.unpackbits to convert such arrays to bits.

In [4]:
test_uint8_arr = np.array([0, 1, 2], dtype=np.uint8)
test_uint8_arr_bits = np.unpackbits(test_uint8_arr)
print(test_uint8_arr_bits)
# use np.packbits to reverse the operation
test_uint8_arr_repack = np.packbits(test_uint8_arr_bits)
print(test_uint8_arr_repack)

[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0]
[0 1 2]

In [5]:
# chromosome can be acquired rather easily
init_pop_ch = np.unpackbits(init_pop, axis=1)
print('shape:', init_pop_ch.shape)
print('excerpt:\n', init_pop_ch[:10, :])

shape: (50, 40)
excerpt:
[[0 1 0 1 0 1 1 0 0 0 1 1 0 1 0 0 0 0 1 1 0 0 1 1 0 0 1 0 0 1 0 1 0 0 1 1
1 1 1 0]
[0 0 1 0 1 0 1 1 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 0 1 0 0 1 0 0
1 1 0 1]
[0 1 0 0 1 0 0 1 0 0 1 0 1 0 0 1 0 1 0 1 0 1 0 0 0 1 0 1 0 0 1 0 0 0 0 0
1 0 0 0]
[0 0 0 0 1 1 1 1 0 1 0 1 0 0 0 0 0 0 0 0 1 1 1 0 0 1 0 1 0 1 1 1 0 0 1 1
1 1 1 0]
[0 1 0 1 0 1 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 1 0 1 1 1 1 1 0 0 0 0
0 0 1 1]
[0 0 0 1 0 0 1 1 0 0 1 0 0 1 1 1 0 0 1 0 1 1 0 0 0 1 0 0 0 0 1 1 0 1 0 0
0 1 0 1]
[0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 1 1 0 1 0 1 0 0
0 0 0 1]
[0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1 0 1 1 0 1 0 1 1 0 1 1 0 0 1 1
1 0 0 1]
[0 0 0 1 1 1 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 1 0 0 0 0 1
1 0 0 1]
[0 0 1 1 0 1 1 1 0 0 1 1 0 0 0 0 0 0 0 1 1 0 0 1 0 1 0 0 1 0 0 0 0 0 0 1
0 1 1 1]]

Computing fitness score¶

We can let the fitness score be $\frac{1}{f(\boldsymbol{x}) + 1}$. When $f(\boldsymbol{x})$ increases, the fitness score decreases, which is exactly what we want.

In [6]:
def fitness_score(bitarr):
# get numerical array, convert to np.int for calculation
arr = np.packbits(bitarr).astype(np.int)
# allow both positive and negative values
arr -= 127
# computing f(x) without loops!
fx = np.sum(np.square(arr))
return 1.0 / (fx + 1.0)

Selection¶

In [7]:
def selection(pop_ch):
#if pop_ch.shape[0] < POP_SIZE: return pop_ch

# apply fitness score function on each sample
fitnesses = np.apply_along_axis(fitness_score, axis=1, arr=pop_ch)
total_fitness = np.sum(fitnesses)

# compute fitness proportions
fitness_prop = fitnesses / total_fitness
#print(fitness_prop)

# generate partial sums with np.cumsum
# no need to compute m repetitively
prop_cumsum = np.cumsum(fitness_prop)
#print(prop_cumsum)

samples = []
# begin sampling
for _ in range(POP_SIZE):
t = np.random.uniform(0.0, 1.0)
for i in range(prop_cumsum.size):
if t < prop_cumsum[i]:
samples.append(pop_ch[i, :])
break

return np.asarray(samples)

Crossover¶

In [8]:
def crossover(pop_ch):
# partition the population into two groups
males, females = np.array_split(pop_ch, 2, axis=0)

# get indices of each group
male_ind = np.arange(males.shape[0])
female_ind = np.arange(females.shape[0])
assert male_ind.size * female_ind.size >= POP_SIZE

# keep track of combinations that are used
used_ind = set()

result = []

# here we assume the size of Cartesian product is far greater than
# our desired sample size. Therefore, the probability of used index
# collision is very low
for _ in range(POP_SIZE):

# generate an unused combination
comb = None
while True:
x = np.random.randint(0, male_ind.size - 1)
y = np.random.randint(0, female_ind.size - 1)
if (x, y) not in used_ind:
comb = (x, y)
break

male = males[comb[0], :]
female = females[comb[1], :]

t = np.random.uniform(0.0, 1.0)
if t < CROSS_RATE:
# partition male and female sample (at the middle)
male1, male2 = np.array_split(male, 2)
female1, female2 = np.array_split(female, 2)
# crossover two samples
new1 = np.concatenate([male1, female2])
new2 = np.concatenate([female1, male2])
result.append(new1)
result.append(new2)
else:
result.append(male)
result.append(female)

return np.asarray(result)

Mutation¶

In [9]:
def flip_bit(bit):
return 1 if bit == 0 else 0

def mutation(pop_ch):
# get a copy of pop_ch to avoid modifying the original array
result = pop_ch.copy()
# get total number of bits in the population
total_bits = pop_ch.size
# generate a random number for each bit
random_probs = np.random.uniform(0.0, 1.0, size=total_bits)
# find out indices that are less than mutation probability
mutation_indices = np.where(random_probs < MUTATION_RATE)[0]
# convert flat index to 2D index
pop_indices = np.unravel_index(mutation_indices, pop_ch.shape)
# flip the bits
for i in range(mutation_indices.size):
idx = (pop_indices[0][i], pop_indices[1][i])
#print(idx)
result[idx] = flip_bit(result[idx])
return result

Putting things together¶

In [10]:
# a function to generate statistics of each generation
def get_gen_stat(pop_ch):
fitnesses = np.apply_along_axis(fitness_score, axis=1, arr=pop_ch)
return fitnesses

In [11]:
pop_ch = init_pop_ch
stats = []
for _ in range(N_GEN):
pop_ch = selection(pop_ch)
pop_ch = crossover(pop_ch)
pop_ch = mutation(pop_ch)
stats.append(get_gen_stat(pop_ch))

# add one final selection to limit the output size
pop_ch = selection(pop_ch)

In [12]:
final_pop = np.packbits(pop_ch, axis=1).astype(np.int)
print(final_pop - 127)

[[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0 -1]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[-1  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[-1  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[-1  0  0  0  0]
[ 0  0  0  0  0]
[-1  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]
[ 0  0  0  0  0]]

In [13]:
# graph generation stats
max_fitnesses = [np.max(arr) for arr in stats]
min_fitnesses = [np.min(arr) for arr in stats]
avg_fitnesses = [np.mean(arr) for arr in stats]

xs = np.arange(len(max_fitnesses))

plt.plot(xs, max_fitnesses, 'b', label='max fitness')
plt.plot(xs, min_fitnesses, 'r', label='min fitness')
plt.plot(xs, avg_fitnesses, 'g', label='average fitness')
plt.legend()
plt.xlabel('generation')
plt.ylabel('fitness')
plt.show()

In [ ]: