NYCJUG/2022-11-08
Beginner's regatta
Here we examine the usefulness of right-to-left evaluation for non-commutative functions like subtraction and division. This order of evaluation allows us to easily evaluate an alternating sum, in the case of "-/", or, less usefully, an alternating product by "%/".
Apparently it's not easy to produce an alternating sum in other languages, this discussion from StackOverflow offered as evidence: https://stackoverflow.com/questions/52880477/alternating-sums-of-a-list
Write a program that accepts 9 integers from the user and stores them in a list. Next, compute the alternating sum of all of the elements in the list. For example, if the user enters 1 4 9 16 9 7 4 9 11 then it computes 1 – 4 + 9 – 16 + 9 – 7 + 4 – 9 + 11 = –2
myList = [] value = None count = 0 while count != 9: value = int(input("Please enter a number: ")) myList.append(value) count = count + 1 if count == 9: break print(myList) def newList(mylist): return myList[0] - myList[1] + myList[2] - myList[3] + myList[4] - myList[5] + myList[6] - myList[7] + myList[8] x = newList(myList) print(x)My code returns the correct answer, but I need it to print out the actual alternating sums as in the example. I have been stuck on this for a while. I am having a mental block on this and havent been able to find anything similar to this online. I appreciate any help or tips. Also, this is python 3. Thank you.
In J or APL, the answer to this is -/. Since a string of subtractions is evaluated right to left, it becomes an alternating sum, as we see below.
Reduction with Non-commutative Functions
To illustrate what how J does this, consider this simple example:
-/1 2 3 4 _2
First we evaluate 3-4, then subtract this result from 2: 2-(3-4). Since J is evaluating right to left, the above expression essentially becomes "(1-2)+(3-4)".
Because evaluation is right to left, something as simple as -/ gives us an alternating sum, which shows up frequently in certain kinds of arithmetic progressions. For example, an expression for 1/e is given like this: [from https://math.hmc.edu/funfacts/e-is-irrational/]
1/e = 1 – 1/1! + 1/2! – 1/3! + …
Looking back to the expression to calculate the reciprocal of e, in J we could write this alternating sum for y terms using extended precision as:
13 : '-/%!2+i.x: y' [: -/ [: % [: ! 2 + [: i. x:
We could use this to evaluate the base of the natural logarithm e by taking the reciprocal of the result to an arbitrary precision using extended precision integers in J:
calce=: 13 : '%-/%!2+i.x: y' calce 10 3991680r1468457 20j18":calce 10 NB. Format rational to 20-digit decimal with 18 places past the decimal 2.718281842777827338
How it Works
Let's take a look at our series generator in detail. Before we apply minus reduction to get our alternating sum, let's look more closely at how we generate this series which is the reciprocal of successive factorials.
So, 10 successive factorials would be !i.10 and the reciprocal of each would be %!i.10. Once we have the underlying series of terms, we simply -/ them. Here are the numbers, their factorials followed by their reciprocals on the bottom. Notice how the middle row corresponds to the denominators in the above equation for 1/e.
(i. 10),(!i.10),:%i. 10 0 1 2 3 4 5 6 7 8 9 1 1 2 6 24 120 720 5040 40320 362880 _ 1 0.5 0.333333 0.25 0.2 0.166667 0.142857 0.125 0.111111
So the full expression for e to 10 terms is
%-/%!i.10 2.71828
So it looks like we are on the right track. Now all we have to do is ensure that the argument is an extended integer by applying x: and ensure we work on scalars, giving us:
calce=: [:%-/@(%@!@i.)@x:"0
Formatting the result to display the extended precision number to 18 digits past the decimal point:
20j18":calce 10 2.718281842777827338
Comparing this to the known value (which is accurate to less than 16 digits):
20j18":^1 2.718281828459045091
Once we exceed the number of digits J provides us in floating point, how do we evaluate our precision? This is an interesting problem to be dealt with later. For now, let's see how many iterations we need to calculate e to the default floating point precision and slightly beyond that.
20j18":^1 2.718281828459045091 20j18":calce 20 2.718281828459045235 20j18":calce 30 2.718281828459045235
We see that there is no difference between arguments of 20 and 30, indicating we have greater precision than the 18 digits past the decimal displayed here. Let's compare to 28 digits:
30j28":calce 20 2.7182818284590452353665869897 30j28":calce 30 2.7182818284590452353602874714 30j28":^1 2.7182818284590450907955982984
In fact, our basis of comparison does not give us 18 digits, so the last two or three are suspect at least, and in fact wrong as shown by reference to other sources, like this fun site: https://apod.nasa.gov/htmltest/gifcity/e.1mil . The first 400 or so digits they display are:
#'2.71828182845904523536028747135266249775724709369995' 52
Comparing this to the first 50 digits we calculate:
2.71828182845904523536028747135266249775724709369995
52j50":calce 20
2.71828182845904523536658698967595342428154276021909
52j50":calce 30
2.71828182845904523536028747135266252501193123422785
52j50":calce 40
2.71828182845904523536028747135266249775724709369996
So we get to almost 50 digits with only 40 iterations.
Alternating Product
The alternating product does not seem to be as widely used as the alternating sum. There appears to be this exception but I do not understand what it means:
The alternating product has applications throughout mathematics. In differential geometry, one takes the maximal alternating product of the tangent space to get a canonical line bundle over a manifold. Intermediate alternating products give rise to differential forms (sections of these products over the manifold).
However, we can use the same behavior we saw with minus reduce but with division to easily give us another useful set of mathematical series called continued fractions. Looking at this essay, we see that these series are useful for calculating many extended precision values using the fork +%.
The continued fraction for the square root of 2 is
√2= 1+1/(2+1/(2+1/(2+1/(2+...))))
or
So, to calculate this in J, we can do this:
(+%)/ 1,10$2 1.41421
Looking at more digits:
20j18":(+%)/ 1,10$2 1.414213551646054778
Comparing this with a known value to a limited number of digits:
20j18":%:2 1.414213562373095145
Using more terms to get better precision:
20j18":(+%)/ 1,20$2 1.414213562373094923 20j18":(+%)/ 1,30$2 1.414213562373095145
As with the calculation of e above, we see that we need about 30 terms to match the floating point precision.
Show-and-tell
How Much Precision?
In the previous section, we saw how to calculate a number of quantities to an arbitrary number of digits. This is fine for quantities that are easy to check, like e or the square root of 2. However, how do we know what our precision is for a more obscure quantity?
Assuming we have an accurate converging series, we can easily figure out how many digits each additional so many terms buys us. Let's start with our continued fraction version for the square root of 2.
calcSqrt2=: 13 : '(+%)/ 1,y$2x'"0
We force extended precision by using the extended version "2x". The argument to this function is the number of terms (+1) we want to use.
So, comparing our approximations with the known floating point value:
20j18":&> (%:2),calcSqrt2 20 30 1.414213562373095145 1.414213562373094923 1.414213562373095145
Interestingly, we can use this to show us the internal precision of the J calculation using the built-in square root function %::
102j100":&> (%:2),calcSqrt2 30 40 1.4142135623730951454746218587388284504413604736328125000000000000000000000000000000000000000000000000 1.4142135623730951454746218587388284504413604736328125000000000000000000000000000000000000000000000000 1.4142135623730951454746218587388284504413604736328125000000000000000000000000000000000000000000000000
Here we see that concatenating the floating point number to the extended precision ones forces the entire array to floating point, degrading the precision, and that we seem to have about 52 digits of precision in floating point:
#'1.4142135623730951454746218587388284504413604736328125' 54 <pre> This seems plausible if we are using extended precision floating point in the interpreter. Alas, this is too good to be true. Subsequent comparison to our continued fraction calculation shows us that the apparent precision past the first 16 digits or so are spurious. <pre> (102j100":%:2),102j100":&>calcSqrt2 30 40 1.4142135623730951454746218587388284504413604736328125000000000000000000000000000000000000000000000000 1.4142135623730950488016834827469619979965358813986413927534448659660525668915837365375964417201725072 1.4142135623730950488016887242095822171633115114204027766433406990747923271591003220117048108247971040
Comparing the floating-point to the 30-term version:
0 i.~=/(102j100":%:2),:102j100":&>calcSqrt2 30 NB. Where does floating-point first differ from 30-term version? 17
So %:2 gives us only 14 (16-#'1.') accurate digits past the decimal. We can test that the 30-term evaluation has more digits than this by comparing it to the 40-term one:
50j48":&> calcSqrt2 30 40 1.414213562373095048801683482746961997996535881399 1.414213562373095048801688724209582217163311511420
Looking up where they first differ:
=/50j48":&> calcSqrt2 30 40 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 i.~=/102j100":&>calcSqrt2 30 40 25
How Quickly Does Precision Increase?
Now that we have a way to figure out how many digits are good, we can use this to figure out how many digits each additional term adds to our precision.
0 i.~=/50j48":&> calcSqrt2 30 40 25 0 i.~=/50j48":&> calcSqrt2 40 50 32 0 i.~=/50j48":&> calcSqrt2 50 60 40 0 i.~=/50j48":&> calcSqrt2 60 70 48
This looks like each additional ten terms adds seven or eight digits of decimal precision.
We can check if this seems true for many more terms. Note how we have to display these to more digits to compare them properly.
0 i.~=/1000j998":&> calcSqrt2 100 200 77 0 i.~=/1000j998":&> calcSqrt2 200 300 155 0 i.~=/1000j998":&> calcSqrt2 300 400 231 0 i.~=/1000j998":&> calcSqrt2 400 500 308
Taking the successive differences for each increase of 100 terms:
2-~/\77 155 231 308 78 76 77
It looks like each additional term gives us about 0.77 more decimal digits.
Another interesting statistic to consider is how much more time it takes to add 100 more terms.
tms=. 6!:2 'calcSqrt2 100' tms=. tms,6!:2 'calcSqrt2 200' tms=. tms,6!:2 'calcSqrt2 300' tms=. tms,6!:2 'calcSqrt2 400' tms=. tms,6!:2 'calcSqrt2 500' tms 0.0031639 0.0181793 0.0452892 0.0880701 0.150514 2%~/\tms 5.74585 2.49125 1.94462 1.70902
Generalizing this:
#tms=. tms,,6!:2 &> (<'calcSqrt2 '),&.>":&.>600+100*i.10 15 2%~/\tms 5.74585 2.49125 1.94462 1.70902 1.56226 1.45142 1.39857 1.34634 1.31036 1.27864 1.26327 1.22964 1.22426 1.21314
This tells us that each additional 100 terms takes 1.2 times as long to compute.
Advanced topics
Work with Multi-threading
Here's a summary of some code I put together over the past few months to use multi-threading to parallelize a time-consuming process, which is code to examine photos and transpose, or "flip", them to be upright.
NB.* startFlipping: cover function for multi-threaded photo flipper "flipN".
startFlipping=: 3 : 0
12 startFlipping y
:
y=. '\' (],[ #~ [~:[:{:]) y rplc '/\'
subds=. 0{"1 dir y,'*.' NB. Process all subdirs' photos
if. 0 e. $subds do. subds=. ,<'' end.
svdir=. 1!:43'' NB. Keep original directory to return to.
if. -.nameExists 'MUTEX' do. MUTEX=: 10 T. 0 end. NB. Fast mutex
if. x<nt=. 1 T. '' do. NB. Have as many threads as requested?
maxt=. 1{8 T. '' NB. Maximum # threads supported
createThreads (maxt-nt)<.0>.x-nt NB. As many more as we need, up to max.
end.
tmst=. tmsd=. '' NB. Timings for threads, subdirs
for_sd. subds do. sd=. ;sd NB. Process each sub-dir
1!:44 ] y,sd
if. 0~:#FLS=: 0{"1 dir '*.jpg' do. tmsd=. tmsd,6!:1''
tmst=. tmst,<flipN x NB. Save thread times, run x of them
tmsd=. tmsd,~6!:1''
end.
end.
1!:44 svdir
tmst;tmsd
NB.EG startFlipping PICDRVguess,':/DCIM/' NB. PICDRVguess is drive letter of photo chip.
)
This uses the function below to create the threads if they have not already been created. We count the number of existing threads with 1 T. '' and restrict the number we create according to the maximum allowed using 1{8 T. .
createThreads=: 3 : '{{ 0 T. '''' }} ^:y]'''''
Now we look at the cover for the core routine where we spin off the threads in a for loop:
flipN=: 3 : 0
tms=. i.0
while. y<#FLS do. pyx=. i.0
for_ix. i. y<.#FLS do.
if. 0<#FLS do. pyx=. pyx,flip1Photo t. '']ix end.
tms=. tms,-/&>_2{.&.>pyx
end.
end.
NB.EG tms=. flipN 10 NB. Run on 10 threads
)
Finally, the deepest function that concerns multi-threading:
NB.* flip1Photo: flip a single photo.
flip1Photo=: 3 : 0
assert. *./nameExists &> 'MUTEX';'FLS' NB. Need these 2 globals: FLS=: 0{"1 dir '*.jpg' [ MUTEX=: 10 T. 0 NB. Fast mutex
tntxt=. ' '-.~":tn=. 3 T. '' [ tm=. 6!:1'' NB. Current thread number is "tn". "tm" is session start time.
if. 0<#FLS do.
myFl=. 'ThisFL',tntxt [ 11 T. MUTEX NB. Lock mutex
(nms)=: split FLS [ nms=. myFl,' FLS' NB. Take 1 filename from list
13 T. MUTEX NB. Unlock mutex
if. 0~:#>".myFl do. flipPhotoByOrientation >".myFl end.
end.
tn,(6!:1''),tm NB. Thread #, end, start session time.
NB.EG flip1Photo t. '']0
)
The above uses this to give each thread a different random seed. This is important for simulations as we want to run numerous different simulations, not the same one repeatedly.
randomizeSeed=: 3 : '9!:1](|-/7!:1'''')+(3 T. '''')+<.13#.|.(6!:1''''),(6!:4''''),6!:0'''']y' NB. Add thread number...
This innermost threading function is the only place we use the mutex primitive 13 T. MUTEX based on the previously-created mutex variable MUTEX. It would seem prudent to have multiple mutexes if there are multiple places in the code where we might need this but I have not tested how well this works.
Some Experiments and Results
This code has been used daily for about three months without error on a couple of different machines. At this point, it seems reliable.
The only set of experiments done recently with this code was one where we compared flipping the photos directly on the SD card versus first moving them to a hard-drive and flipping them there. The result of these experiments were that flipping on the SD card seemed as reliable as flipping on a hard-drive but it was quite a bit slower. This is probably due to the relative speed of reading and writing from an SD card versus a hard-disk with an SSD cache.
Learning and Teaching J
Locklin's Comparison of Q versus J
Scott Locklin recently published an essay comparing J to Q. The TLDR section starts with this
J is a larger and more powerful programming language than Q in many ways (as a language, not as a tool). Q is a lot easier to use, especially if your problem set falls into the broad category of timeseries/trading problems. J is more complete in that there are a lot of ideas baked into J; it is a huge language in comparison to Q, but Q solves specific problems.
Locklin proposes that J's handling of higher-dimensional arrays is overkill. He asserts that "[m]ost sane people work with 2-d arrays and tables work just fine for this purpose....' Those of us who work with photos might take issue with this. Q essentially works only with tables expressed as dictionaries with equal-length entries. He continues this exposition in the context of J's rank-centric array handling, something useful for working with higher-dimensional arrays but not so much for the "very few people [who] will ever in their lives use rank>2 arrays."
Another difference between the two languages is Q's profusion of data types. This perhaps overlooks J's many numeric types but these difference reflect the different orientations of the languages. J was developed from the perspective of mathematics and, perhaps, number theory. Q is dominated by efficiency considerations which is why it has data types which sacrifice generality for specificity. In J, we don't worry about integers being promoted to floating point as they exceed the integer domain but this does lead to performance overhead and variability.
He makes a good point about J namespaces, that they are flexible and useful, but says that the name of the current namespace is less obvious in J than in Q which apparently displays the current locale in the command-line interface whereas J requires you to enter "coname ".
In the discussion of "point-free coding", which we know as "tacit" in J, he mentions using J's "13 :" to create a tacit expression written with explicit variables. I do this frequently myself, as we saw above in the "Beginner's Regatta", as I find the use of explicit names more readable, so I think this is a good point. In fact, it relates to a discussion in the APL community about direct definitions in Dyalog APL. These are pieces of code in curly braces with α (lower-case alpha) for the left argument and ω (lower-case omega) for the right argument, like this nested expression:
','{1↓¨(1,⍺=⍵)⊂⍺,⍵}'comma,separated,list'
comma separated list
The consensus in the APL community seems to be that having named arguments is more readable than inferring them.
Code Comparison
The example code he provides seems a bit strange. He compares this J verb
histogram =: <: @ (#/.~) @ (i. @#@[ , I.)
to this Q function:
histogram:{@[(1+max x)#0;x;+;1]}
Of the J, he says "[t]his is the kind of thing I would write (I didn’t) as it uses lots of parens @’s and tildes" but that "Q is more direct; it basically acts like there are J-@’s everywhere, except you don’t have to write the bloody things...."
He continues
I dunno maybe most people find these equivalently obtuse (the [;;;] stuff ends up doing a lot of interesting things that are not obvious any more than [, is in J), but the Q version is a lot more straightforward to me; and this sort of parenthesis @ soup makes me wish hooks and forks were a sublanguage in J rather than the normal way to do things.
I agree that most people would probably find both expressions hard to comprehend but that speaks more to familiarity with each language than to some context-free notion of readability. Additionally, I find the J example to be a bit odd and not terribly useful as a histogram creator or as an example. I can't speak to the Q version, but the J version takes an interesting approach that's worth looking at. However, after seeing its limitations, I could not resist the urge to come up with a better example, one that handles
]rr=. 10?@$10 6 3 9 0 3 3 0 6 2 9 histogram rr 6 3 9 0 3 3 0 6 2 9 NB. Not what I was expecting 0 5 10 histogram rr NB. Supposed to be dyadic but what is the left argument? 2 4 4 (i.#rr) histogram rr NB. This is more useful 2 0 1 3 0 0 2 0 0 2 /:~rr 0 0 2 3 3 3 6 6 9 9
Here's something I put together in an attempt to make a more J-like and more useful version.
histogram1=: 13 : '<:2-~/\I. (i.#x) e.~/:x,y' histogram1 [: <: 2 -~/\ [: I. ([: i. [: # [) e.~ [: /: , 0 5 10 histogram1 rr NB. Buckets as left argument 6 4 2 5 10 histogram1 rr 4 4
This histogram1 does a better job than histogram for a number of argument types.
histogram 'abbcccdd'
|domain error: histogram
|histogram[0]
dbr 1 [ dbr 0
histogram1 'abbcccdd'
0 0 0 0 0 0 0
'abcd' histogram1 'abbcccdd'
1 2 3
'abcde' histogram1 'abbcccdd'
1 2 3 2
0 0.2 0.4 0.6 0.8 1 histogram 0.1*rr
2 1 3 0 2 2
0 0.2 0.4 0.6 0.8 1 histogram1 0.1*rr
2 4 0 2 2
/:~0.1*rr NB. To help judge the correctness of the result
0 0 0.2 0.3 0.3 0.3 0.6 0.6 0.9 0.9
'abcd' histogram 'abbcccdd'
1 2 3 2
APL Histogram Example
Just for something completely different, here is some (lightly edited) APL code I recently found for drawing a histogram. It looks very much like something I wrote myself back when large sheets of ice covered much of North America.
X←6 3 9 0 3 3 0 6 2 9
' ⎕'[¯1 0↓(⌽⍳(~⎕IO)+⌈/A)∘.≤A←+/((⌊/X)+⍳1+(⌈/X)-⌊/X)∘.=X]
⎕
⎕ ⎕ ⎕ ⎕
⎕ ⎕⎕ ⎕ ⎕
X[⍋X] ⍝ To check result by inspection
0 0 2 3 3 3 6 6 9 9
This handles numeric vectors and assumes the buckets are unit lengths but also does display a useful result.
The Problem with Hooks and Forks
Locklin makes the valid point that "...in J, you literally have to count verbs to figure out how the composite train works; then you have to look to the left and do it all over again at least once, because hooks and forks behave differently in dyadic and monadic use." This applies when we are reading a fragment of J code but not so much when we are writing it. He goes on to point out that "...we now have the dissect toolkit to figure out what’s going on with these snarls, but, well, you need something called “dissect” to read the code."
He avers that he is "...definitely of the Art [Arthur Whitney] camp on this one; everything should be atop; hooks and forks should be some kind of sublanguage for mentats. I can do it, I just mostly wish I didn’t have to."
This may be a case of a beautiful theory - function composition - slain by an ugly fact: the difficulty of interpreting sequences of more than about three symbols at a time.
Lack of Advanced Documentation
Locklin has good things to say about the value of NuVoc but makes the point that "[t]here are no advanced books on how to use J; they are all unfortunately beginners books." This contrasts with some of the very good advanced books on q like the ones by Nick Psaris, like this one on machine learning and this one on tips for building maintainable Kdb code.
In a way, some of what we do here with NYCJUG attempts to remedy this by providing details on worked-out problems but he does have a good point that there is little guidance above the level of basic code management.
In conclusion, it seems like Locklin's main argument is that the differences between the languages stem from J's greater generality as opposed to Q's focus on a certain class of problems. He does raise some good points about the problem with reading tacit expressions or long hooks or forks and the lack of advanced documentation for J. These are things we should work on to improve J.
Teaching Practices Inventory
Here is a short summary of some teaching practices advocated at the University of British Columbia. These are the sort of things to keep in mind as we develop better pedagogical materials for J.
Adding open-ended questions
Depending on the goals for using the Teaching Practices Inventory, open-ended questions can be added at the end of the inventory. For example, at UBC we added 3 questions:
- What do you see as the biggest barrier to achieving more effective student learning in the courses you teach?
- What changes could be made at UBC to help you teach more effectively?
- What changes could be made at UBC to increase your satisfaction/enjoyment of teaching?
Teaching Practices Inventory categories
- Course information provided (including learning goals or outcomes)
- Supporting materials provided
- In class features and activities
- Assignments
- Feedback and testing
- Other (diagnostics, pre-post testing, new methods with measures, …)
- Training and guidance of teaching assistants
- Collaboration or sharing in teaching