text:
The Grand Locus / Life for statistical sciences

## Are you human?

On the Internet, nobody knows you're a dog.

This is the text of a famous cartoon by Peter Steiner that I reproduced below. This picture marked a turning point in the use of identity on the Internet, when it was realized that you don't have to tell the truth about yourself. The joke in the cartoon pushes it to the limit, as if you do not even have to be human. But is there anything else than humans on the Internet?

Actually yes. The Internet is full of robots or web bots. Those robots are not pieces of metal like Robby the robot. Instead, they are computer scripts that issue network requests and process the response without human intervention. How much of the world traffic those web bots represent is hard to estimate, but sources cited on Wikipedia mention that the vast majority of email is spam (usually sent by spambots), so it might be that humans issue a minority of requests on the Internet.

In my previous post I mentioned that computers do not understand humans. For the same reasons, it is sometimes difficult for a server to determine whether it is processing a request issued by a human or by a robot. In the early days of the Internet, even million dollar sites like Yahoo! and Amazon were not protected against massive poisonous requests and were taken down by a 17 year old school boy. It is around that time that appeared the CAPTCHAs (Completely Automated Public Turing test to tell Computers and Humans Apart) which, as the name says is a reverse Turing test. The gist of reverse Turing tests is to process the request only if the client can perform a task that is very easy for a human, and very difficult for a computer (like recognizing a distorted text, a picture with a house or the meaning of a simple sentence).

In short, telling robots apart from humans is important for security reasons. Yet, most servers do not want to deny access to every robot because search engines like Google use robots called web spiders or web crawlers to index web pages. The current agreement is that servers should indicate their policy in a page called robots.txt (out of curiosity you can check the robots.txt page of the blog, but it only contains the address of the site map). The content tells the robots which pages they should not request, but does not prevent them to do so in any way. Not surprisingly, most spammers or petty hackers do not take the time to read the robots.txt page... perhaps some do not even know it exists. So robots could issue a request to every page anyway, right? Well, let's check that out. In the technical section below I show a very simple Python script to get the reviews of user 2467618 on IMDB.

If Python is installed on your computer (which will be the case if you are using Mac or Linux) you can start a Python session (version 2.x) and use the module urllib as shown.

import urllibcontent = urllib.urlopen(    'http://www.imdb.com/user/ur2467618/comments?order=date&start=0').read()f = open('downloaded_content.html', 'w')f.write(content)f.close()

You can now open the file downloaded_content.html in your home directory with your favorite browser to see what it contains.

In case you are not Python-proficient, you can check out what the script retrives here. Among others, you will notice that it says "Access denied http: 403". Sure enough, the robots.txt file of IMDB says that requests to /user are disallowed. So how do servers protect themselves from unwanted queries issued by robots?

There is no universal answer, but in the case of IMDB like in many others, the answer lies in the HTTP headers. You might notice in the "Access denied" page that it says at the bottom "Browser: Python-urllib/1.17". The issue here is that by default, urllib is honest about the user agent, which is easily intercepted and denied by the server. If we decide to lie about our user agent and claim we issue the request through Chrome, we would do as indicated in the the following technical part instead.

Here is an improved version of the script which masquerades the HTTP headers to pretend we are using Chrome to request the page.

import urllib2import cookielibcookies = cookielib.LWPCookieJar()handlers = [    urllib2.HTTPHandler(),   urllib2.HTTPSHandler(),   urllib2.HTTPCookieProcessor(cookies),]   opener = urllib2.build_opener(*handlers)headers = {   'Accept': 'text/html,application/xhtml+xml,'\             'application/xml;q=0.9,*/*;q=0.8',   'Accept-Charset': 'ISO-8859-1,utf-8;q=0.7,*;q=0.3',   'Accept-Encoding': 'gzip,deflate,sdch',   'Accept-Language': 'en-US,en;q=0.8,fr;q=0.6',   'Connection': 'keep-alive',   'User-Agent': 'Mozilla/5.0 (X11; Linux i686) AppleWebKit/535.19 '\                 '(KHTML, like Gecko) Ubuntu/12.04 '\                 'Chromium/18.0.1025.151 Chrome/18.0.1025.151 '\                 'Safari/535.19',}request = urllib2.Request(    url='http://www.imdb.com/user/ur2467618/comments?order=date&start=0',     headers=headers)connection = opener.open(request)content = connection.read()# The content is gzip-compressed.f = open('downloaded_content.html.gz', 'wb')f.write(content)f.close()

As you can check by decompressing the file downloaded_content.html.gz, we get the same content as if we had issued the request from Chrome.

To download the page, we had to set the HTTP headers to their values when we request the page from Chrome, which makes the code substantially more complicated. We can get the values of those headers very easily in Chrome by clicking on the wrench tool in the top-right corner and then choosing Tools > Developer Tools. This displays a console which has a "Network" item where all requests are analyzed and where you can find the values of all the HTTP headers, data and cookies. You can also use Firebug on Firefox, or a sniffer like Wireshark to analyze the traffic to and from your browser.

By setting those headers to the values they have when the request is issued by a human through a browser, it is much harder to recognize that the request actually comes from a script. One thing still that most servers check is the frequency and the regularity of the requests. I don't know any human who would issue over nine thousand requests with less than a second interval. Neither do system administrators, and this is why they might block the IP address those requests are issued from. Basically, by masquerading HTTP headers, and breaking the regularity patterns in the requests, it becomes very difficult to distinguish humans from robots without CAPTCHAs.

If you ever wondered, the answer is yes: this is what I did to fetch the reviews from IMBD. But then, why would I ever be interested in user 2467618? This is what I will expand on in my next post.


text:
# CONSTRUCTIONS FOR SPARSE ROW-ORTHOGONAL MATRICES WITH A FULL ROW

• Published : 1999.03.01

#### Abstract

In [4], it was shown that an n by n orthogonal matrix which has a row of nonzeros has at least ( log2n + 3)n - log2n +1 nonzero entries. In this paper, the matrices achieving these bounds are constructed. The analogous sparsity problem for m by n row-orthogonal matrices which have a row of nonzeros in conjectured.

#### References

1. Bull. of Korean Math. Soc. Constructions for the sparsest orthogonal matrices G.-S. Cheon;B. L. Shader
2. J. of Combinatorial Theory Series A v.85 How sparse can a matrix with orthogonal rowa be? G.-S. Cheon;B. L. Shader
3. Rocky Mtn. Math. J. A simple proof of Fiedler's conjecture concerning orthogonal matrices B. L. Shader
4. Combinatorial and GraphTheoretical Problems in Linear Algebra Combinatorial Orthogonality L. B. Beasley;R. A. Brualdi;B. L. Shader;R. A. Brualdi(eds.);S. Friedland(eds.);V. Klee(eds.)
5. Discrete Applied Math. J. Sparse orthogonal matrices and the Haar wavelet G.-S. Cheon;B. L. Shader

format: plain


text:
## Project Euler: Problem 79

The statement of problem 79 from Project Euler has been reproduced in its entirety below.

A common security method used for online banking is to ask the user for three random characters from a passcode. For example, if the passcode was 531278, they may ask for the 2nd, 3rd, and 5th characters; the expected reply would be: 317.

The text file, keylog.txt, contains fifty successful login attempts.

Given that the three characters are always asked for in order, analyse the file so as to determine the shortest possible secret passcode of unknown length.

I did not find this problem difficult in terms of either math or programming, but I really enjoyed it because the solution illustrates how a small amount of information about something valuable and secret, a passcode, can compromise its integrity. Although knowing the shortest possible passcode does not grant an attacker an entry, it increases the probability that an attack is successful. Unfortunately, many average users choose the shortest passcode satisfying all of the requirements imposed at the time of creation (e.g. Password1).

As usual, I have implemented my solution in Python.

We first need to read the contents of the text file keylog.txt. We then need to extract the set of passcode samples from the contents of the text file. We also need the set of the digits that appear in the passcode samples; this provides a lower bound on the length of the passcode.

L1 = [] #Initialize a list to hold all integer instances
L2 = [] #Initialize a list to hod all tuple instances

f = open(file_path, 'r') #Create a file object in read-mode

for line in f: #loop through the lines of the text file
digits = map(int, line[0:3]) #Create a list of the integers in the sample
L1.extend(digits) #Append the digits to the list L1
L2.append(tuple(digits)) #Create a tuple out of digits and append to L2

f.close() #Close the file object

S1 = set(L1) #Use the set object to eliminate duplicates
S2 = set(L2)

We initialize two empty lists $L1$ and $L2$. To read the text file, we create a file object $f$. Since a file object is iterable, we can loop over the lines of $f$. Python reads the lines of the text file, not surprisingly, as strings. Since we are interested in integers, we use the function map with arguments int and. The int function will return an error when called if the entire line is passed as an argument due to the trailing newline character ‘\n’, so we only pass the first three characters of the string block. For example, the line ‘123’ will yield [1, 2, 3]. We store the result in the variable digits. We append the three digits that appear in a given sample to the running list $L1$ using the extend method of lists. We choose to represent each sample as a 3-tuple of integers using the function tuple. Since we don’t want to append the elements of the tuple to the list $L2$, but rather the tuple itself, we use the append method of lists. After we have completed the loop, we are done with the file object $f$ and can close it. To eliminate duplicate digits and tuples, we use the set data structure.

d = {1 : set(),
0 : set(),
3 : set(),
2 : set(),
7 : set(),
6 : set(),
9 : set(),
8 : set()}

for x in S2:

We initialize a dictionary $d$ whose keys are the integers represented as integers and whose values are empty sets. The sets will store all of the digits that appear after a given a key. We iterate over the elements of the set $S2$. We look up the dictionary value of the first element of the tuple and add the second and third elements of the tuple to the set. We then look up the value of the second element of the tuple and add the third element of the tuple to the set. When the loop completes, the values of each key will contain all the digits that we know appear after the key based on the 50 samples.
Since there are eight unique digits in the “keylog” sample, the shortest possible passcode must have length greater than or equal to $8$. No digits appear after $0$, so we place it as the last character of the passcode. $0$ appears only after $9$, so we place $9$ as the second-to-last character. Continuing in the fashion, we obtain the candidate passcode
$\displaystyle s_{1}s_{2}\cdots s_{8}=\text{73162890}$
We cannot change any of the characters $s_{i}$ without eliminating a digit, thereby invalidating one of the samples. We cannot add characters while preserving the order without increasing the length. Hence, 73162890 is the shortest possible passcode from which the samples could have come.


text:
1. Mar 30, 2007

Ry122

does cos + sin = 0?

2. Mar 30, 2007

neutrino

cos + sin can never be equal to anything other than cos + sin. sin and cos are just functions. Neither has a "value" as such.

Do you want to if sin(somevalue1) + cos(somevalue2) = 0?

3. Mar 30, 2007

sutupidmath

like neutrino said cos + sin have no meaning. But if we get let's say

sin(2n pi)+ cos[(pi/2)+n pi] =0, and there are more samples like this. if this is what you meant at first place?????

4. Mar 30, 2007

f(x)

did you mean $$cosx+sinx=0$$ ?
$$\mbox{for the above to be true},\\ cosx=-sinx \\ \Rightarrow cotx=-1$$ (or $$sinx=0$$ but this does not hold for that.)

5. Mar 30, 2007

ChaoticLlama

Please only post calculus questions in this forum.

6. Mar 30, 2007

f(x)

LOL,this may well be calculus.....
$$Max(sinx+cosx)= \sqrt{2}$$
$$Min(sinx+cosx)= -\sqrt{2}$$
Since they function is continuous, it must have x such that f(x)=0

7. Mar 30, 2007

Gib Z

$$\cos^2 x + \sin^2 x = 1$$


text:
# Is this a typo in Section 1.8.1 Mathematics for Computer Science?

Am I completely mistake or is it suppose to say $$n^2$$ is a multiple of 2 and therefore $$n$$ must be a multiple of 4?

This is from MIT's Mathematics for Computer Science

• Generally: If $m,n$ are positive integers and $m^{1/n}$ is not an integer then $m^{1/n}$ is irrational. Proved by a similar method. – DanielWainfleet Apr 26 at 1:48
• What makes you think this is a mistake? – Servaes Apr 26 at 7:31

No, it is correct. The point is that $$2d^2=n^2$$ implies $$n^2$$ is even, and only even numbers square to give an even number, so $$n$$ much be even, so $$n^2$$ is then actually a multiple of $$4$$.
No. $$6^2$$ is a multiple of $$2$$ but $$6$$ is not a multiple of $$4$$. If $$n=2k$$ then for sure $$n^2=4k^2$$. So, MIT is right.


text:
380 views

Consider the following grammar:

• $S \rightarrow S$
• $S \rightarrow SS \mid a \mid \epsilon$

Construct the collection of sets of LR (0) items for this grammar and draw its goto graph.

edited | 380 views
0
i think this is ambigous grammer try to gererate string "a" you will gate 2 LMD
0
It is ambiguous grammar and we can't parse it by LR(0) parser. right na?
0

Here  S→  ϵ will take as reduce/ shift action or not effect on LR(0) DFA??

0
So we cannot write LR(0) items if grammar is ambiguous?
+1
afaik, LR(0) items can be constructed if the grammar is ambiguous.

+1 vote

The augmented production is $S^{'} \rightarrow S$.

$\textbf{GOTO Graph:}$

Here, each of $I_0$, $I_1$, $I_2$, $I_3$ is a set of $LR(0)$ items. And hence $I_0$, $I_1$, $I_2$, $I_3$ are the collection of sets of $LR(0)$ items.

by Boss (15.2k points)
selected by
+1 vote

Augmented production S'-->S

by Boss (11.8k points)
reshown
0
Sir we can write LR(0) items..

But there will be conflicts in dfa (of canonical collection of LR(0) items).. so the grammar will not be LR(0)
+2
Thats all is asked in the question rt?
0
Yes..


text:
Math in Focus Grade 5 Chapter 9 Practice 7 Answer Key Real-World Problems: Decimals

Practice the problems of Math in Focus Grade 5 Workbook Answer Key Chapter 9 Practice 7 Real-World Problems: Decimals to score better marks in the exam.

Math in Focus Grade 5 Chapter 9 Practice 7 Answer Key Real-World Problems: Decimals

Question 1.
Mrs. Lee uses 0.025 kilogram of wax to make a candle. On Monday, she made 50 candles. On Tuesday, she made 4 times as many candles as on Monday. How much wax did she use to make the candles on Tuesday?
Explanation:
given information,
on Monday 0.025 km of wax use to make 50 candles.
on Tuesday x kg of wax used to make 4 x 50 = 200 candles.
how much wax ‘x ‘ is to be find out.
0.025 ÷ 50 = 0.0005
0.025 ÷ 50
=0.0005 kilograms of wax is used for one candle.
for 4 x 50 = 200 candles.
0.0005 x 200 = 0.1 kg of wax is used.

Question 2.
One lap of a race track measures 4.68 kilometers. During a race of 56 laps, a driver stops to refuel after completing 48 laps. How many more kilometers does he have to drive to finish the race?

37.2 kilometers
Explanation:
one lap = 4.65 kilometers
56 x 4.65 = 260.4 kilometers
48 x 4.65 = 223.2 kilometers
37.2 more kilometers does he have to drive to finish the race.
260.4 – 223.2 = 37.2 kilometers

Question 3.
Mrs. Rahlee bought 300 yards of ribbon to make flowers. She used 1.22 yards to make one large flower. She made 200 such large flowers She used all of the remaining ribbon to make 100 small flowers. What was the length of ribbon Mrs. Rahlee used to make one small flower?
0.56 yards length of ribbon Mrs. Rahlee used to make one small flower.
Explanation:
total 300 yards of ribbon to make flowers.
1.22 yards for one flower.
200 x 1.22 = 244 yards used to make one large flower.
300 – 244 = 56 yards
56 ÷ 100 = 0.56
0.56 yards length of ribbon Mrs. Rahlee used to make one small flower.

Question 4.
Britta bought some carrots and apples for $24.80. A carrot and an apple cost$0.90 altogether. She bought more carrots than apples. The cost of the extra number of carrots was $6.80. How many apples did Britta buy? Answer: 20 apples Explanation: c carrot a apple x is number of carrots y is number of apples c + a = 0.90 x c + y a = 24.80 6.80 + y a = 24.80 y a = 24.80 – 6.8 = 18 But a+c=.90 y=18÷ 0.9 = 20 y is the number of apples Solve. Show your work. Question 5. A plastic tub has a capacity of 1 3.5 quarts. It can hold 3 times as much liquid as a pail. The pail can hold twice as much liquid as a can. Find the capacity of the pail and that of the can in quarts. Answer: 2.25 quarts Explanation: capacity of tub = 13.5 quarts. Tub can hold three times as much as a pail. capacity of pail = 13.5 ÷ 3 = 4.5 quarts. The pail can hold twice as much liquid as a can. capacity of can =4.5 ÷ 2 = 2.25 quarts. Question 6. Marcy paid$35 for 10 kilograms of raisins. She divided the raisins equally into two containers. Then she sold the raisins in the first container at $4.50 per kilogram and those in the second container at$5.50 per kilogram. How much money did Marcy earn after selling all the raisins?
Answer: $15 Explanation: Since she divided raisins equally both containers have 5kg each. Total earnings from selling container one = 4.5 * 5 =$22.5
Total earnings from selling container two = 5.50 * 5 = $27.5 Total earnings from both containers = 22.5 + 27.5 =$50
Total cost of raisins = $35 Total Profit =$50 – $35 =$15


text:
# Subsonic and supersonic diffusers have the following geometry

This question was previously asked in
ISRO Refrigeration and Air Conditioning 2020 Official
View all ISRO Scientist ME Papers >
1. Divergent and convergent respectively
2. Both divergent
3. Both convergent
4. Convergent and divergent respectively

Option 1 : Divergent and convergent respectively
Free
ST 1: General Awareness
6847
15 Questions 15 Marks 15 Mins

## Detailed Solution

Concept:

All flows must satisfy the continuity and momentum relations as well as the energy and state equations.

Application of the continuity and momentum equations to a differential flow yield:

$$\frac{{dA}}{A} = \left( {M{a^2} - 1} \right)\frac{{dV}}{V}$$

This equation governs the shape of a nozzle or diffuser in subsonic or supersonic isentropic flow.

Diffusers:

Diffusers are designed to decelerate fluids, therefore dV is negative, dV < 0

From the above equation,

dV < 0, Ma > 1 (supersonic) ⇒ dA < 0, convergent

dV < 0, Ma < 1 (subsonic) ⇒ dA > 0, divergent

Therefore, for diverging diffuser, Flow should be subsonic.

Nozzles:

Nozzles are designed to accelerate fluids, therefore dV is positive, dV > 0

From the above equation,

dV > 0, Ma > 1 (supersonic) ⇒ dA > 0, divergent

dV > 0, Ma < 1 (subsonic) ⇒ dA < 0, convergent


text:
# How to construct this collision vector?

I have a 3d game, but for movement I use only two axes (X and Z) - X is vertical, Z is horizontal

Here's what I'm trying to achieve:

Player is the circle, InputDirection is Blue, NewInputDirection is Green. I'm struggling on getting the desired green direction.

Image 1 The player moves up (he's pushing onto the wall)

InputDirection: (-1f, 0f, 0f)

NewInputDirection: (0f, 0f, 0f) because we should stop him.

Image 2 The player moves towards the wall and the bottom. When he reaches the wall, his horizontal axis should be blocked - set to 0, so he slides)

InputDirection: (0.5f, 0f, -0.5f)

NewInputDirection: (0.5f, 0f, 0f)

Image 3 Same thing as above, but the wall is rotated.

Right now, my code looks like this:

RaycastHit hit;
bool collides = Physics.Raycast(this.playerRigidbody.position, inputDirection, out hit, 1f); // ray in front of the player

if(collides){ // if there's a collision
// calculate the direction from the player to the collision contact point
Vector3 collisionNormal = (hit.point - playerRigidbody.position).normalized;
// idk...
}

I don't know how to calculate it. I've tried various ways. Should I also take the wall's rotation into account?

The simple answer here is vector projection. If n is a unit vector (length 1), then...

float projection = Vector3.Dot(n, someVector);

is the signed length of the component of someVector that's parallel to n. We can subtract that component to leave only the perpendicular component:

Vector3 perpendicular = someVector - n * projection;

It looks like you have more going on here than just the vector math question, so I'll expand on that below.

You seem to be using this to let the player character move up to and along a wall, without penetrating through it.

For that purpose, I'd recommend using a spherecast or capsulecast instead. The reason is that the ray can squeeze into spaces that your player character couldn't, so if you just position yourself at the end of the ray, you might clip into an obstacle the ray barely missed to one side or the other.

RaycastHit hit;
bool collides = Physics.SphereCast(
playerRigidbody.position,
inputDirection,
out hit,
1f);

if(collides){ // if there's a collision
// I've found that adding a slight extra margin here helps avoid unwanted intersections.
float uninterruptedDistance = hit.distance - 0.001f;
// Compute our blue arrow.
Vector3 directTravel = inputDirection * uninterruptedDistance;

// Subtract out the component of the direction that goes into the wall.
// Leaving the direction that's parallel.
float projection = Vector3.Dot(inputDirection, hit.normal);
// Compute our green arrow direction.
Vector3 slideDirection = (inputDirecton - hit.normal * projection).normalized;

// If you want to move as far as you can along the blue, then start moving along the green...
float remainingDistance = moveDistance - uninterruptedDistance;
Vector3 slideTravel = slideDirection * remainingDistance;

Vector3 totalTravel = directTravel + slideTravel;

// Do what you want with it.
}

Note that we scanned for obstacles along only the blue arrow - so it's possible moving to the end of the green arrow will put us in collision with something we didn't scan for. So if you want to be really exacting here, you can recurse and fire a new ray/sphere/capsule cast along the new direction. I show an example of doing that in this answer.

• Thank you for the answer, really helped me to understand what's going on. I'm just wondering, about the provided link (detecting other objects, eg. when you're standing in the corner between 2 walls). Wouldn't it be better to implement SphereCastAll` in the above example, instead of shooting multiple raycasts as in the link? – Jacob Aug 21 '19 at 16:10
• The whole point of my question was that I don't want to "correct" the player's movement, but the direction he had chosen. In the function, I'm given an input direction and a speed. I need to use these two to construct a new direction vector which will not collide with the environment. Wouldn't it be better to just iterate through all RaycastHit's in the SphereCastAll and subtract the colliding direction of each of the collisions from the input direction? Theoretically, it should work. – Jacob Aug 21 '19 at 16:41
• No, SphereCastAll is redundant for that purpose. It keeps scanning past the first collision, finding all possible collisions along that line. But you want to stop at the first collision - any obstacles past that, your player never hits. After you pick a direction that avoids that collision, you're now looking at moving in a different direction. You have no guarantees about how far the object can move in this new direction unless you do a new query pointing in this direction. – DMGregory Aug 21 '19 at 18:19


text:
## Introduction

Since Soddy first suggested the existence of isotopes in 19131, isotopes have drawn a great deal of attention due to its application in the fields of physics, chemistry, biomedicine, and geology. Generally, isotopes with different nuclear masses could change the energy levels within isotopic atoms and molecules, thus lead to a frequency shift in the atomic or molecular spectrum2, which has been widely used to identify the species of the isotopes and to investigate the static structure of the isotopologues. Moreover, for isotopic molecules, the vibrational modes depend sensitively on the masses of its constituent isotopic atoms, which provides an important method to determine the mechanism of a chemical reaction via the kinetic isotope effect3,4, namely, the fact that heavier isotopes tend to react more slowly than lighter ones. However, a real-time measurement of the motions of the isotopic atoms in molecule (molecular vibration) is a long-standing challenge over the last century, due to the awesome rapidity of the molecular vibration.

Recent advances in strong-field physics have provided efficient approaches to probe both the molecular structure and dynamics using the table-top laser. These new methods rely on the recollision of an electron, removed from the molecule by a strong laser field, with its parent ion. The molecular structure and dynamics are encoded in the amplitude and phase of the emitted high-order harmonics. It stimulates the development of high-order harmonic spectroscopy (HHS)5,6,7,8,9 as an emerging tool for ultrafast detection with femtosecond to attosecond time resolutions. Apart from HHS, some other techniques based on strong-field ionization, such as photoelectron holography and photoelectron diffraction10,11,12,13,14,15,16,17, and so on, have also been demonstrated to image the molecular structure and dynamics. Up to now, many works have been carried out to investigate the effects of nuclear motion in strong-field ionization18,19,20,21 and molecular high-order harmonic generation (MHOHG)22,23,24,25,26. In 2005, Lein showed theoretically that the laser-driven nuclear motion will introduce an amplitude modulation (AM) (see Fig. 1) in harmonic signals via the nuclear autocorrelation function27, which denotes the overlap between the initial and time-dependent nuclear wave function that evolves from the moment of ionization until the recollision. By analyzing the AMs in high-order harmonic generation (HHG) from isotopic molecules (H2 and D2), the intracycle nuclear dynamics has been theoretically predicted27 and experimentally detected28,29. Nevertheless, this method is restricted because the propagation and other inherent physical factors, such as two-center interference27,30 and energy-dependent rescattering cross sections31, may affect the harmonic intensity. Moreover, in the presence of intense lasers, the nuclear motion will lead to larger internuclear distances R and a decrease in the ionization potential Ip18,19,20, which can result in an increase in the ionization rate and thus a strong AM in MHOHG. These factors complicate the retrieval of nuclear dynamics by AM. Apart from AM, frequency modulation (FM) is an alternative way commonly used in various applications, e.g., signal processing and telecommunications. By considering the frequency shift in the atomic spectrum of isotopes, it stimulates us to ask whether the nuclear motion in intense laser fields can induce a frequency shift in the MHOHG spectrum (see Fig. 1). Compared to AM, FM is more stable and insensitive to the laser parameters provided that the ionization saturation is avoided and the pulse length is properly adopted. It thus can provide an alternative powerful way to identify the nuclear dynamics. After the prediction by Bian and Bandrauk in ref. 32, FM has received a lot of attention in theoretical studies31,33,34,35. However, the FM in isotopic MHOHG has never been observed in experiment and the measurement of nuclear motion based on the FM is not addressed.

In the present work, we report the experimental observation of FM in MHOHG from isotopic molecules H2 and D2. High-order harmonics generated from isotopic molecules show obvious red shift with respect to those from Ar atom. The red shift is demonstrated to originate from the laser-induced nuclear motion of isotopic molecules, which strengthens harmonic emission at the trailing edge of the laser pulse. From the observed frequency shift, the nuclear motions of H2 and D2 are successfully retrieved, which agree well with the calculations from non-Born–Oppenheimer time-dependent Schrödinger equation (NBO–TDSE).

## Results

### Experimental observation of FM

The experiment is carried out by adopting a Ti:sapphire laser, and H2 and D2 molecules (see the methods). These isotopes have attracted extensive interest as a prototype. Figure 2a–c displays the spatially resolved harmonic spectra generated from atomic gas Ar and the hydrogen isotopes H2 and D2, respectively. Their ionization potentials are very close. The spatially integrated HHG signals are presented by the dash-dotted (Ar), solid (H2), and dashed (D2) lines in Fig. 2d. One can see that the harmonic intensities from D2 are higher than those from H2, which is in consistent with previous studies27,28,29,36,37. More importantly, the measured harmonics from H2 and D2 present obvious frequency shift with respect to those from Ar. As shown in Fig. 2e–h, each harmonic from both H2 and D2 shows a red shift relative to that from Ar. While for D2, the frequency shift is larger than that of H2.

### Gas pressure dependence of FM

HHG in gas medium includes the individual response, as well as the copropagation of laser and harmonic fields. The propagation effect can possibly induce a frequency shift in HHG38,39,40, which depends sensitively on the gas pressure. However, in our experiment the ionizations of the three gases are weak (below 4%), and also the gas pressure is low. Then the frequency shift induced by the propagation effect will be inappreciable. To check this effect, we measured the harmonic spectra generated from Ar, H2, and D2 at different gas pressures. With the gas pressure changing from 15 to 35 torr, the intensity of each harmonic from these three gases exhibits a quadratic increase, which indicates a good phase matching in our experiment. More than that, the central wavelengths of each harmonic from the three gases are nearly unchanged as shown in Fig. 3a–c. For a clear insight, in Fig. 3d, we present the central wavelength of H17 for Ar (diamonds), H2 (squares), and D2 (circles) as a function of the gas pressure. The frequency shift of H2 and D2 relative to that of Ar keeps almost constant as the gas pressure varies. These results indicate that the influence of propagation effect on the harmonic frequency shift is negligible in our experiment. Besides, the experimental conditions used for HHG from H2 and D2 are exactly the same, the differences in the harmonic spectra can be mainly attributed to the individual response of isotope molecules in the driving laser field.

### Theoretical simulation of FM

It has been reported that the nonadiabatic effect of the time-dependent laser intensity can induce a blue or red shift when HHG is dominant at the leading or trailing edge of the laser pulse41,42,43. For H2 (or D2), the ionization rate depends sensitively on the internuclear distance R18,19. Due to the laser-driven nuclear motion, the average internuclear distance at the trailing edge can be larger than that at the leading edge of the laser pulse, which makes the ionization, as well as the HHG signals stronger at the trailing edge, and therefore induces a red shift in the harmonic spectrum. In contrast, since the laser intensity used in our experiment is far smaller than the ionization saturation threshold of Ar, the ionization and HHG of Ar atom mainly occurs at central part of the laser pulse and is symmetric with respect to the pulse center (t = 0). Then no obvious net shift exists in the harmonics from Ar. Since the nuclear dynamics is avoided for Ar, it can serve as a benchmark to evaluate the frequency shift of harmonics from the two isotopic molecules. For a given harmonic order, the frequency shift caused by the nonadiabatic effect can be obtained via the time derivative of the laser pulse, namely, $$\Delta \omega = \alpha _q\frac{{\partial I(t)}}{{\partial t}}|_{t = t_{\mathrm{i}}}$$44,45,46. Here, $$I(t) = I_0{\mathrm{exp}}\left( {\frac{{ - 4{\mathrm{ln}}(2)t^2}}{{\tau ^2}}} \right)$$ with I0 = 1.5 × 1014 W cm−2, τ = 30 fs is the envelope of the laser pulse, ti is the ionization moment of the given harmonic, and α q is its phase coefficient, which can be evaluated according to the strong-field approximation (SFA) model47. The time derivative with a positive (negative) sign means a blue (red) shift of this harmonic. Owing to the slower nuclear motion of heavier nuclei, the dominant harmonic emission of D2 occurs later than that of H2 (see Supplementary Note 1 and Supplementary Figure 1). As a result, the HHG from D2 experiences a more rapid change of the effective laser intensity (namely, a larger value of $$\left| {\frac{{\partial I(t)}}{{\partial t}}} \right|$$), which therefore gives rise to a larger red shift in the harmonic spectrum as observed in our experiment. Besides the nonadiabatic effect, the nuclear motion can lead to the variation of the ionization potential and the complex recombination dipole, which may affect the harmonic phase accumulated during the electron excursion and influence the MHOHG48,49. To evaluate these influences, we have performed simulations with the modified SFA model48, which indicates that the frequency shift induced by these two effects is far smaller than our experimental observations (see Supplementary Note 2 and Supplementary Figures 2 and 3). Moreover, the laser-driven nonadiabatic alignment may also lead to a red shift in MHOHG. We have evaluated this influence by considering the time-dependence of the laser-driven alignment under our experiment condition (see Supplementary Note 3 and Supplementary Figure 4). Our calculations show that the red shift induced by the molecular alignment is about one order of magnitude smaller than our experimental observations. Note also that the fluctuation of the laser carrier-envelope phase, which is not fixed in our experiment, will not affect the measured frequency shift of MHOHG because a multi-cycle laser pulse is used in our experiment. Therefore, the main contribution to the frequency shift shown in Fig. 2 is attributed to the nonadiabatic effect induced by the nuclear motion32.

In Fig. 4a, we present the relative frequency shift of H15-H23 for H2 (squares) and D2 (circles). The relative frequency shifts gradually decrease as the harmonic order increases. The experiment is also simulated by solving the NBO–TDSE (see the methods). The calculated frequency shifts of H2 and D2 are presented by the dashed lines in Fig. 4a, which are in agreement with the experimental observations. Some difference in quantity may arise from the uncertainties of experimental parameters. In Fig. 4b, we have calculated the asymmetry coefficients of HHG signals for H2 (squares) and D2 (circles). Here, the HHG asymmetry is defined as η(ω) = (P+(ω) − P(ω))/(P+(ω) + P(ω))31, where $$P_ + (w) = {\int}_0^{ + \infty } g(\omega ,t)\mathrm{d}t$$ and $$P_ - (w) = {\int}_{ - \infty }^0 g(\omega ,t)\mathrm{d}t$$ are the amount of harmonic ω generated at the trailing and leading edges of the laser pulse, respectively. g(ω, t) is the time-frequency spectrogram calculated with the Gabor transform. In our calculation, the width of the time window used in the Gabor transform is 0.1 fs, which corresponds to a filter with the width of 10ω0 (ω0 is laser frequency) in the frequency domain. One can see that for each harmonic order, the harmonic emission is more pronounced at the trailing edge of the laser pulse (namely, η(ω) > 0). Therefore, all the harmonics exhibit a red shift in the spectrum as shown in Fig. 4a. Moreover, the trend of the red shift agrees qualitatively with the asymmetry coefficients. This agreement suggests that the observed red shift indeed results from the delayed emission of HHG with respect to the center of the laser pulse.

### Monitoring the nuclear dynamics by FM

As mentioned above, the frequency shift mainly arises from the asymmetry of the ionization (and so HHG) with respect to the center of laser pulse (t = 0 fs)32. Previous studies have shown that the ionization rate of H2 (D2) is approximately linearly dependent on internuclear distance R before it reaches 2 a.u.50. Thus the relative frequency shift of the q-th harmonic with respect to 0 can be estimated as:

$$\frac{{\Delta \omega }}{{\omega _0}} = \frac{{\mathop {\sum}\limits_{t_{\mathrm{i}} < 0} R(t_{\mathrm{i}}) - \mathop {\sum}\limits_{t_{\mathrm{i}} > 0} R(t_{\mathrm{i}})}}{{\mathop {\sum}\limits_{t_{\mathrm{i}}} R(t_{\mathrm{i}})}},$$
(1)

where ω0 is the frequency of the driving laser, ti is the ionization moment of the electron (contributing to the q-th harmonic generation) in each half optical cycle. ti < 0 and ti > 0 mean the ionization occurs at the leading and trailing edges of laser pulse, respectively. For a given harmonic, ti can be calculated according to the three-step model47,51. To retrieve the nuclear motion, we consider to employ the commonly used linear harmonic oscillator model52 to describe the two-body vibrations of H2 and D2. The simulations with the NBO–TDSE show that the harmonic oscillator model works well in a low-ionization case. In our experiment, the ionization is below 4%, therefore the harmonic oscillator model is applicable. In the harmonic oscillator model, the potential V(r) of H2 (or D2) can be approximatively expressed as $$V(r) = V_0 + \frac{k}{2}(r - R_e)^2$$, where V0 and k are constants and Re is the equilibrium internuclear distance of H2 and D2. Then the laser-driven nuclear motion can be derived in the form of (for details, see Supplementary Note 4)

$$R(t) = A{\mathrm{sin}}({\mathrm{\Omega }}t + \phi ) + BI(t) + R_{\mathrm{e}}.$$
(2)

On the right side of Eq. (2), the first term denotes the inherent harmonic vibration, A, Ω, and ϕ are the corresponding amplitude, frequency, and phase of the vibration. The second term represents the laser-nucleus interaction. Inserting Eq. (2) into Eq. (1), the frequency shift of a specific harmonic can be expressed as a function of A, B, Ω, and ϕ. By fitting the observed frequency shifts of H15-H23 to Eq. (1) with the least square method, the four parameters can be determined. Then the nuclear motion R(t) can be retrieved. Figure 4c–d shows the retrieved nuclear vibrations (dash-dotted line) of H2 and D2, respectively. As shown in this figure, the maximum of the retrieved R(t) of H2 is about 1.505 a.u., which is slightly larger than that of D2 (1.485 a.u.). Moreover, the retrieved R(t) of H2 oscillates with a period of 8.2 fs. In contrast, it is 11.4 fs for D2. The retrieved oscillation periods of H2 and D2 are very close to the vibrational periods of H2 and D2 in their ground electronic state (7.5 and 10.6 fs). The ratio of these two retrieved periods is also very close to the expected mass ratio of $$\sqrt 2$$. The results calculated from the NBO–TDSE are also presented as the solid lines. From Fig. 4c, d, one can see that due to the inherent harmonic vibration of the harmonic oscillator model, the retrieved nuclear motion shows much deeper modulation at the beginning when compared to the simulated one. While with the increase of the laser intensity, the simulated nuclear motion also turns to oscillate after t = − 20 fs due to vibrational excitation. Despite the initial oscillation, the main structures of the retrieved nuclear motion R(t), such as the dynamic range and the overall trend, can agree well with the theoretical predictions in the range of [−20, 20] fs where most of the HHG signals are generated. It should be explained that to compare with the NBO–TDSE calculations, the initial internuclear distance used in the experimental fitting is obtained from the NBO–TDSE simulation [namely, the initial values of the solid lines in Fig. 4c, d]. It is given by the expectation value of R with the ground state wavefunction Ψ0, namely, 〈Ψ0|R0〉. Note that the expectation values 〈Ψ0|R0〉 are different for H2 and D2 due to their different field-free Hamiltonians (depending on the nuclear mass). Moreover, the so-called equilibrium internuclear distance Re is defined as the minimum of the BO potential VBO(R), namely, where dVBO(R)/dR = 053. Since the BO potential VBO(R) is slightly asymmetric with respect to Re, the expectation value of R obtained from the NBO–TDSE simulation (1.44 a.u. for H2 and 1.43 a.u. for D2) is slightly different from the so-called equilibrium internuclear distance Re (1.4  a.u. for both H2 and D254). To study the stability of the retrievals, we have also performed the fitting with different harmonic orders or using the known values of the oscillation frequency Ω. The obtained results are all in good agreement with the NBO–TSDE simulations (see Supplementary Note 4 and Supplementary Figures 5 and 6). Considering the simplicity of the harmonic oscillator model and the uncertainty of experimental parameters, the agreement of the retrieved nuclear motions with the TDSE predictions is very satisfying.

## Discussion

In summary, we experimentally observed the red shift in HHG from isotopic molecules H2 and D2. The red shift is primarily attributed to the laser-driven nuclear motion in H2 and D2, which strengthens the ionization rate and harmonic emission due to larger internuclear distance R and lower Ip at the trailing edge of the laser pulse. By using a linear harmonic oscillator model, the nuclear vibrations of H2 and D2 are successfully retrieved from the observed frequency shift. The FM effect in MHOHG is universal, which can be directly applied to other light molecules if the ionization rate is sensitive to nuclear motion. Moreover, in our experiment the molecules are not pre-aligned, the alignment effect is negligible. In principle, the FM technology can be extended to aligned molecules with any alignment angles with respect to the laser polarization. The alignment-angle-dependent FM can not only be used to extract the ultrafast electron-nuclear dynamics, but also be possible to image molecular structure.

In previous studies of AM, the intensity ratios of HHG from isotopic molecules reveal the nuclear dynamics of $${\mathrm{H}}_2^ +$$ and $${\mathrm{D}}_2^ +$$ within the time window from ionization to recombination in one laser cycle, namely, intracycle dynamics. In contrast, in the present work, the observed frequency shift provides a monitoring of the nuclear vibrations of H2 and D2 at each ionization moment in the laser pulse, namely, intercycle dynamics. Therefore, FM in MHOHG reveals a different physical process and is complementary with the method of AM27,28,29 for probing the nuclear dynamics. These findings may provide a deep insight into some of the most fundamental events in chemistry and facilitate the development of HHS.

## Methods

### Experimental methods

The experiment is performed by using a commercial Ti:sapphire laser system (Legend Elite-Duo, Coherent, Inc.), which delivers the 30 fs, 800 nm pulses at a repetition rate of 1 kHz. The output laser pulse is focused to a 2-mm-long gas cell by a 600-mm focal-length lens. In Fig. 2, the stagnation pressure of the gases is 30 torr and the gas cell is placed 2 mm after the laser focus to ensure the phase matching of the short quantum path. The laser energy used in our experiment is maintained at 1.5 mJ and the corresponding intensity is estimated to be 1.5 × 1014 W cm−2. The generated harmonic spectrum is detected by a homemade flat-field soft x-ray spectrometer consisting of a flat-field grating (1200 grooves mm−1) and a slit with a width of about 0.1 mm and height of 15 mm. High-order harmonics are dispersed by the grating and imaged onto the microchannel plate (MCP) fitted with a phosphor screen. The image on the screen is read out by a CCD camera.

To accurately evaluate the frequency shift in MHOHG, we have calibrated the spectrometer by using the atomic lines of carbon in terms of a procedure similar to that in refs. 55,56. The atomic lines are produced by focusing several millijoules of the driving laser pulse to interact with a 0.5-mm-thick graphite sheet placed at the position where HHG occurs. We record the generated atomic lines and read their coordinates on the phosphor screen. By assigning the observed atomic lines to the known literature data of carbon, we can then achieve the calibration of the spectrometer. Details of the calibration are provided in Supplementary Methods.

### Theoretical methods

To simulate the HHG process and nuclear dynamics of H2 and D2, we numerically solve the NBO–TDSE with one active electron27,57. Since the electron and nuclear motions follow the linearly polarized laser field, we adopt the one-dimensional model,

$$i\frac{{\partial \Psi (z,R,t)}}{{\partial t}} = [H_{\mathrm{e}}(t) + H_{\mathrm{n}}(t) - E(t)z]\Psi (z,R,t),$$
$$H_{\mathrm{n}} = - \frac{1}{{2\mu }}\frac{{\partial ^2}}{{\partial R^2}} + \frac{1}{R},$$
$$H_{\mathrm{e}} = - \frac{1}{2}\frac{{\partial ^2}}{{\partial z^2}} + V_{{\mathrm{en}}}(R,z),$$
$$V_{{\mathrm{en}}}(R,z) = - \frac{{Z(R,|z + R/2|)}}{{\sqrt {(z + R/2)^2 + 0.5} }} - \frac{{Z(R,|z - R/2|)}}{{\sqrt {(z - R/2)^2 + 0.5} }} + V_{{\mathrm{BO}}}^ + (R).$$
(3)

Here, z is the electron coordinate, R is the internuclear distance, and E(t) is the driving laser field. Hn and He are the Hamiltonians for the nuclei and electron, respectively. Ven(R, z) is the Coulomb potential of the electron–nucleus interaction. μ is the reduced mass of two nuclei and $$V_{{\mathrm{BO}}}^ + (R)$$ is the lowest BO potential of H$$_2^ +$$. In order to faithfully mimic the nuclear dynamics of H2 and D2, we have adopted an effective nuclear charge $$Z(R,\xi ) = \left[ {1 + e^{ - \xi ^2/\sigma ^2(R)}} \right]/2$$, where σ(R) is an R-dependent screening parameter. By adjusting the parameter σ(R) at each internuclear distance, the energy of the ground electronic state of hydrogen and hydrogen ion can be well reproduced.

The TDSE is solved using the B-spline method with Crank–Nicolson scheme58. The initial wave function is taken as the ground state of the hydrogen molecule calculated using the imaginary-time propagation method. The wave function is numerically discretized in the grid with a size of −50 ≤ z ≤ 50 a.u. and R ≤ 18 a.u. The discrete steps for the electron coordinate and nuclear coordinate are δz = 0.1 and δR = 0.05 a.u., respectively. The temporal step is δt = 0.05 a.u. Then the harmonic spectrum can be obtained by performing the Fourier transform of the laser-induced electron dipole moment d(t) = 〈Ψ(z,R,t)|z|Ψ(z,R,t)〉. The expectation of the internuclear distance is calculated by $$R(t) = \frac{{\left\langle {\Psi \left( {z,R,t} \right)\left| R \right|\Psi \left( {z,R,t} \right)} \right\rangle }}{{\left\langle {\Psi \left( {z,R,t} \right)\Psi \left( {z,R,t} \right)} \right\rangle }}$$59.

To confirm the results of the one-electron model, we have also solved the NBO–TDSE with two electrons22,23, which reads

$$i\frac{{\partial \psi (z_1,z_2,R,t)}}{{\partial t}} = \left[ {H_{\mathrm{e}} + H_{\mathrm{n}} + E(t)(z_1 + z_2)} \right]\psi (z_1,z_2,R,t),$$
$$H_{\mathrm{n}} = - \frac{1}{{2\mu }}\frac{{\partial ^2}}{{\partial R^2}} + \frac{1}{R},$$
$$H_{\mathrm{e}} = - \frac{1}{2}\frac{{\partial ^2}}{{\partial z_1^2}} - V_{{\mathrm{en}}}(z_1) - \frac{1}{2}\frac{{\partial ^2}}{{\partial z_2^2}} - V_{{\mathrm{en}}}(z_2) + V_{{\mathrm{ee}}}(z_1 - z_2),$$
$$V_{{\mathrm{en}}}(z_i) = \frac{1}{{\sqrt {(z_i + R/2)^2 + \beta (R)^2} }} + \frac{1}{{\sqrt {(z_i - R/2)^2 + \beta (R)^2} }},$$
$$V_{{\mathrm{ee}}} = \frac{1}{{\sqrt {(z_1 - z_2)^2 + \gamma (R)^2} }}.$$
(4)

Hn and He are the Hamiltonians for the nuclei and electrons. Ven(z i ) and Vee are the Coulomb potential for the electron–nucleus and electron–electron interactions. β(R), γ(R) are the R-dependent softening parameters. Due the huge computation of the two-electron model, we calculate only the harmonic spectrum for H2. The obtained red shift is consistent with the one-electron model, as well as the experimental observations.

### Data availability

All the data that support the findings of this study are available from the corresponding author upon reasonable request.


text:
# Stripping out the Business Logic

published on August 7th, 2013

What I am writing day-in and day-out can probably best be described as business programs. To me business programs are programs whose essential purpose it is to keep track of some state, control processes, make decisions, and communicate with external services. One question that intrigues me is how one can seperate out the business logic from the rest of the program.

Business logic is, as the name implies, about business first and logic second, and since business needs can change on a whim, the business logic needs to change too, usually much more often than the rest of the application. Furthermore, the business logic is the raison d'être of the whole application, so it makes sense to have this logic in one place instead of intermingling it with implementation specific logic. Finally, when I write such programs it always seems I have to do way too much work to achieve what I want. To me, seperating out the business logic into a seperate domain specific language (DSL), module, workflow pattern, or whatever, holds the promise that to create additional functionality I only need to implement that functionality in the business logic DSL or module, reducing the amount of work needed and thus speeding up development considerably.

However, while there is surely some interest in this topic, I have yet to find the best way to achieve this. As far as I can tell there may be a lot of activity, but this activity doesn't result in well communicated solutions that are easy to comprehend. Maybe I just don't get it, but I have yet to find somebody to tell me: this is how you do it. It seems that quite often business logic is simply mixed with the rest of the application. In this post I want to explore several possibilities of how business logic can be seperated from the rest of the application and described concisely.

## Goals

In an ideal world, the representation of the business logic should fulfill the following criteria:

Succinct

: Since one of the objectives is to have the business logic in one place, the representation of it should be as small as possible. Boilerplate and verboseness should be avoided.

Unambiguous

: Each instance of a representation should have only one meaning. This may seem obvious, but as far as I understand this is not necessarily true for BPMN. This means that a rountrip translation (representation $\rightarrow$ program $\rightarrow$ representation) should at least be theoretically possible and lead to equivalent results.

Easy to understand

: The representation of the business logic should convey the logic to a reasonably intelligent reader. The reader should not have to have intimate knowledge with how this representation is translated to a program in order to understand what will happen.

Testable

: Since the representation will contain the entire business logic, which will change often and might be complex, it is important that this logic is easily testable to make sure we got everything right.

Representable

: The representation should be easily representable in the language one is working with. If it is it becomes easier to serialize and deserialize a representation along with its state and to display it to a (support) user.

Complete

: We should be able to represent everything we want to represent in this representation, otherwise we'll have to work around it at which point we would have probably better done without it.

## Some Possible Solutions

### Finite-State Machines

Finite-state machines are sometimes[^1]^,^[^2] suggested for this task. Finite-state machines are well understood, easy to draw and are able to model all processes which only have a, well, finite number of states. On a high level, this should include all sensible business processes. However, once we include storing and modifying arbitrary data in our application, we are clearly out of finite-state machine territory. Therefore, to describe all of the business logic we will clearly need some extension to finite state machines. If that extension is ad-hoc then we might lose the benefits that finite-state machines provide in the first place.

### Workflow Languages

There are also a large number of workflow languages, many of which are found in some enterprise software vendors packages.[^5]^,^[^6] To me the most interesting are the Business Process Modell and Notation^[Business Process Modell and Notation website.] and the Yet Another Workflow Language[^3] which was created by Arthur ter Hofstede, an academic how studies the field and has written an influential review paper.[^4] These languages seem interesting and have clearly some industry weight behind them, it seems fairly difficult to find good introductory material on them that is not (a) a 1000-page tome of a standard or (b) some software vendors sales pitch documents. Some of these languages also have problems, for example there seems to be no general unambiguous way to translate BPMN into a program.[^7]

Overall, I am a little disappointed with the quality of the introductory information that I have found so far. The most promising document seems to be the review paper by ter Hofstede and the YAWL language.

### A Process DSL

What we have done at work is that we have basically defined an expression language that is a C++ DSL. With this language you can write expressions such as

if_(customerHasValidEmailAddress())
.else_(sendVerificationEmail() and showEmailVerificationPage())

and so on. Each of the mentioned steps has to be implemented in the host language to perform the desired action. This means of course that the implementation burden is still there, but at least steps that are already defined can be reused in a different context.

One problem I see with this approach is that there is no way to store state inside the process description, in essence this is a language that has boolean expressions and stores everything else in global state. This also means that for every new piece of state, however temporary, requires a new thing that can contain it. I have found that this can limit reuse and sometimes make these processes confusing. But this DSL seems essentially like an implementation of the previously described workflow languages.

### A More Powerful DSL

This is where I would like to go. A DSL that can hold some temporary state, where variables can be bound and reused at a later state. Ideally, this DSL would also permit a way to manage the data model directly, so that it would essentially be able to implement the entire business logic without having to fall back to implementations in the host language all the time. Alas, right now, this is only an idea and I am not sure how it will work.

## Outlook

This post has barely scratched the surface of the task of abstracting out the business language into a seperate representation. I hope I have transmitted the idea of why this is desireable and outlined a few of the possible approaches that I have found during my inital research of the topic. In future posts I hope that I'll be able to explore some of these approaches in more depth.

[^1]: My question on stackoverflow about a year ago on the same topic.

[^2]: Shopify technical blog: Why developers should be force-fed finite state machines.

[^4]: Workflow Patterns review paper by of workflow pattern systems by Arthur ter Hofstede.

[^6]: OS Workflow.


text:
# How is mass determined from an empirical formula?

Well, normally we determine $\text{empirical formula}$ from a mass...
...or at least from $\text{a percentage composition}$.


text:
(Note: This tutorial is a fork of the official dask tutorial, which you can find here)

In this tutorial, we will use dask.dataframe to do parallel operations on dask dataframes look and feel like Pandas dataframes but they run on the same infrastructure that powers dask.delayed.

Let’s start by installing dask with:

conda install -c conda-forge dask

dask has many backends, but the one that’s most versitile (and often most performant) is the distributed backend, so let’s use that here. Run:

from dask.distributed import Client
client = Client()
client

to create a set of workers on your computer. This cluster is running only on your own computer, but it operates exactly the way it would work on a cloud cluster where these workers live on other computers.

When you type client in a jupyter notebook, you should see the clusters status pop up like this:

This status tells me that I have four processes running on my computer, each of which running 4 threads (for 16 cores total). In addition, it tells me that the cluster sees that I have about 34GB of memory for it to play with.

One cool thing you’ll see is that this panel includes a link to a Dashboard. Click this link, and a new tab will open in your browser that shows you, in real time, what your dask cluster is doing! We’ll come back to this later.

Note: these workers communication via networking protocols, so your firewall / anti-virus may ask if it’s ok to let Python accept connections (it is!). If you can’t get your system running, it may be due to your firewall / anti-virus. In that case, try client = Client(processes=False). This will create four workers hiding in the same process (obviating the need for network protocols for communication between processes). Don’t worry if that last bit doesn’t make sense – those are implementation details you don’t ahve to worry about! :)

### Real Data¶

Lets try this with an extract of flights in the USA across several years. This data is specific to flights out of the three airports in the New York City area. To begin, download NYC Flight data here, and unzip it somewhere you can find it.

(1) Look at the data in the folder. As you will see it is in many files. Unlike pandas, however, dask can load collections of files all at once. Run import dask.dataframe as dd, then using dd.read_csv with an asterix wildcard to load all the files at once with something like:

import dask.dataframe as dd
parse_dates={'Date': [0, 1, 2]})

(The parse_dates option just tells dask that columns 1, 2, and 3 together form a date (they’re year, month, and day of the month respectively), and that dask should convert them to a single column called Date that is of type datetime. This isn’t special to dask, it’s just a pandas read_csv option.)

(2) Look at your dataframe by just calling it (df). Notice that the respresentation of the dataframe object contains no data - Dask has just done enough to read the start of the first file, and infer the column names and types.

(3) now use df.head() to ask dask to actually load a few rows.

We can view the start and end of the data

(4) Now try and run df.tail().

### What just happened?¶

Unlike pandas.read_csv which reads in the entire file before inferring datatypes, dask.dataframe.read_csv only reads in a sample from the beginning of the file (or first file if using a glob). These inferred datatypes are then enforced when reading all partitions.

In this case, the datatypes inferred in the sample are incorrect. The first n rows have no value for CRSElapsedTime (which pandas infers as a float), and later on turn out to be strings (object dtype). Note that Dask gives an informative error message about the mismatch. When this happens you have a few options:

• Specify dtypes directly using the dtype keyword. This is the recommended solution, as it’s the least error prone (better to be explicit than implicit) and also the most performant.

• Increase the size of the sample keyword (in bytes)

• Use assume_missing to make dask assume that columns inferred to be int (which don’t allow missing values) are actually floats (which do allow missing values). In our particular case this doesn’t apply.

In our case we’ll use the first option and directly specify the dtypes of the offending columns: set TailNum to str, CRSElapsedTime to float and Cancelled to bool.

## Computations with dask.dataframe¶

Suppose we want to compute the maximum of the DepDelay (short for Departure Delays) column. With just pandas, we would loop over each file to find the individual maximums, then find the final maximum over all the individual maximums

maxes = []
for fn in filenames:
maxes.append(df.DepDelay.max())

final_max = max(maxes)

But this is a far cry from what we do with smaller pandas dataset, which is just

df = pd.read_csv(filename, dtype=dtype)
df.DepDelay.max()

To solve this, we can use dask! dask.dataframe lets us write pandas-like code, that operates on larger than memory datasets in parallel. To get the maximium value of the column, we can run:

[ ]:

%time df.DepDelay.max().compute()

CPU times: user 127 ms, sys: 29.9 ms, total: 157 ms
Wall time: 1.98 s

1435.0

This writes the delayed computation for us and then runs it.

Some things to note:

1. Dask uses lazy evaluation, which means it doesn’t actually do any work till we call .compute(). If you just ran df.DepDelay.max(), you’d just get back a placeholder:

> df.DepDelay.max()

dd.Scalar<series-..., dtype=float64>

1. Dask will delete intermediate results (like the full pandas dataframe for each file) as soon as possible.

• This lets us handle datasets that are larger than memory

• This means that repeated computations will have to load all of the data in each time (run the code above again, is it faster or slower than you would expect?)

(5) Want to know more about what dask is doing? First, just store df.DepDelay.max() as a new variable and look at it. Can you see the placeholder? It tells you what it will eventually be (a scalar (single value) float64). But you don’t see the actual maximum value.

(6) Now we can use .visualize() on that variable to see the plan dask has for getting you the maximimum value if you were to run your code right now. First, however, note you’ll need to install the package graphviz by running pip install graphviz (for some reason the conda install doesn’t seem to be working).

## Let’s do some manipulations!¶

In this section we do a few dask.dataframe computations. If you are comfortable with Pandas then these should be familiar. You will have to think about when to call .compute().

When you run .compute(), you can also move over to your Dashboard (remember that link in the panel that openned when you created client?) to see your cluster working!

(7) In total, what is the average delay (DepDelay) for non-canceled flights?

(8) In total, how many non-cancelled flights were taken from each airport?

Hint: use df.groupby and size.

(9) What was the average departure delay from each airport? (DepDelay)

(10) How many rows are in our dataset?

Note you’ve probably just found what I think is one of the trickiest things to figure out with a distributed computing system like dask: some functions don’t follow the rules of delayed evaluation, and when you run them, dask executes immediately. len is one of those functions.

As we’ll discuss below, a good rule of thumb is that anything that isn’t either (a) a method on the dask dataframe (i.e. something that isn’t of the form df.METHOD()), or an explicit dask function (something you call from the dask library, like dd.to_parquet), it will probably execute immediately.

len is pretty safe since all it returns is a single number, so it’s unlikely to blow up your computer if it runs before you want, but you have to be careful with other functions, because if you were to us a pandas function like pd.to_numeric(df), it will execute immediately, and it will do so by pulling all the data into your client session as a single giant dataframe, and with a larger dataset, it is likely to crash you computer when you do so.

An alternative that won’t execute until you run compute is df.shape[0] (see how it takes the form of a method on our dask dataframe df?) Try it to see!

## Sharing Intermediate Results¶

When computing all of the above, we sometimes did the same operation more than once. For most operations, dask.dataframe hashes the arguments, allowing duplicate computations to be shared, and only computed once.

For example, lets compute the mean and standard deviation for departure delay of all non-canceled flights. Since dask operations are lazy, those values aren’t the final results yet. They’re just the recipe required to get the result.

If we compute them with two calls to compute, there is no sharing of intermediate computations.

[ ]:

non_cancelled = df[~df.Cancelled]
mean_delay = non_cancelled.DepDelay.mean()
std_delay = non_cancelled.DepDelay.std()

[ ]:

%%time

mean_delay_res = mean_delay.compute()
std_delay_res = std_delay.compute()

CPU times: user 376 ms, sys: 43.4 ms, total: 419 ms
Wall time: 3.88 s

But lets try by passing both to a single compute call.

[ ]:

import dask

[ ]:

%%time

CPU times: user 237 ms, sys: 26.5 ms, total: 264 ms
Wall time: 2 s

Using dask.compute takes roughly 1/2 the time. This is because the task graphs for both results are merged when calling dask.compute, allowing shared operations to only be done once instead of twice. In particular, using dask.compute only does the following once:

• the calls to read_csv

• the filter (df[~df.Cancelled])

• some of the necessary reductions (sum, count)

To see what the merged task graphs between multiple results look like (and what’s shared), you can use the dask.visualize function (we might want to use filename='graph.pdf' to zoom in on the graph better):

[ ]:

dask.visualize(mean_delay, std_delay)

But lets try by passing both to a single compute call.

(11) Now you try – for each airport, get both the average delay and the maximum delay in the data for flights that weren’t cancelled.

(As another example of places that distributed computing gets tricky, while it’s easy to do means and maxes, you can’t easily get the median of a group (e.g.df.groupby('Origin)['DepDelay'].median()) in a distributed manner. There are workarounds, but these are the kinds of “fragile” things that make distributed computing a bit of a pain).

## Converting Scheduled Departure Time (CRSDepTime) to a timestamp¶

This dataset stores the scheduled departure times of flights in the format HHMM, which are read in as integers in read_csv:

[ ]:

crs_dep_time = df.CRSDepTime.head(10)
crs_dep_time

0    1540
1    1540
2    1540
3    1540
4    1540
5    1540
6    1540
7    1540
8    1540
9    1540
Name: CRSDepTime, dtype: int64

To combine the data we have on flight dates with these departure times into a single column that includes all that information in a single datetime format, we have to convert the time of day into a Timedelta type, then add those Timedeltas to our Date column to get a single column that knows the date AND time of day of departures. (Implicitly, dates think of themselves as living at midnight, so adding a timedelta of 10 hours moves them to 10am).

To convert these to timestamps of scheduled departure time, we need to convert these integers into pd.Timedelta objects, and then combine them with the Date column.

In pandas we’d do this using the pd.to_timedelta function, and a bit of arithmetic:

[ ]:

import pandas as pd

# Get the first 10 dates to complement our crs_dep_time

# Get hours as an integer, convert to a timedelta
hours = crs_dep_time // 100
hours_timedelta = pd.to_timedelta(hours, unit='h')

# Get minutes as an integer, convert to a timedelta
minutes = crs_dep_time % 100
minutes_timedelta = pd.to_timedelta(minutes, unit='m')

# Apply the timedeltas to offset the dates by the departure time
departure_timestamp = date + hours_timedelta + minutes_timedelta
departure_timestamp

0   1990-01-01 15:40:00
1   1990-01-02 15:40:00
2   1990-01-03 15:40:00
3   1990-01-04 15:40:00
4   1990-01-05 15:40:00
5   1990-01-06 15:40:00
6   1990-01-07 15:40:00
7   1990-01-08 15:40:00
8   1990-01-09 15:40:00
9   1990-01-10 15:40:00
dtype: datetime64[ns]

## Other Things to Know¶

### Persist (caching)¶

Sometimes you want your computers to keep intermediate results in memory (when it’s possible). To “cache” data in this manner, you can to use the persist() method to tell dask what results you want to keep around. You can read more here.

### Using Parquet¶

Dask works really efficiently with data stored in parquet format, as noted here, so keeping data in parquet (as I’ve suggested you do anyway!) will help you with dask.

## Limitations¶

### What doesn’t work?¶

Dask.dataframe only covers a small but well-used portion of the Pandas API. This limitation is for two reasons:

1. The Pandas API is huge

2. Some operations are genuinely hard to do in parallel (e.g. sort)

Additionally, some important operations like set_index work, but are slower than in Pandas because they include substantial shuffling of data, and may write out to disk.

### What definitely works?¶

• Trivially parallelizable operations (fast):

• Elementwise operations: df.x + df.y

• Row-wise selections: df[df.x > 0]

• Loc: df.loc[4.0:10.5]

• Common aggregations: df.x.max()

• Is in: df[df.x.isin([1, 2, 3])]

• Datetime/string accessors: df.timestamp.month

• Cleverly parallelizable operations (also fast):

• groupby-aggregate (with common aggregations): df.groupby(df.x).y.max()

• value_counts: df.x.value_counts

• Drop duplicates: df.x.drop_duplicates()

• Join on index: dd.merge(df1, df2, left_index=True, right_index=True)

• Operations requiring a shuffle (slow-ish, unless on index)

• Set index: df.set_index(df.x)

• groupby-apply (with anything): df.groupby(df.x).apply(myfunc)

• Join not on the index: pd.merge(df1, df2, on='name')

• Ingest operations

• Files: dd.read_csv, dd.read_parquet, dd.read_json, dd.read_orc, etc.

• Pandas: dd.from_pandas

• Anything supporting numpy slicing: dd.from_array

• From any set of functions creating sub dataframes via dd.from_delayed.

• Dask.bag: mybag.to_dataframe(columns=[...])