Dyalog APL Problem Solving Competition 2020 — Phase I
+ +Annotated Solutions
+ + + + +Table of Contents
-
+
- Introduction +
- Problem 1 – Take a Dive +
- Problem 2 – Another Step in the Proper Direction +
- Problem 3 – Past Tasks Blast +
- Problem 4 – Bioinformatics +
- Problem 5 – Future and Present Value +
- Problem 6 – Merge +
- Problem 7 – UPC +
- Problem 8 – Balancing the Scales +
- Problem 9 – Upwardly Mobile +
Introduction
+After Phase I, here are my solutions to Phase II problems. The full code is included in the post, but everything is also available on GitHub.
+A PDF of the problems descriptions is available on the competition website, or directly from my GitHub repo.
+The submission guidelines gave a template where everything is defined in a Contest2020.Problems Namespace. I kept the default values for ⎕IO and ⎕ML because the problems were not particularly easier with ⎕IO←0.
++This post is still a work in progress! I will try to write explanations for every problem below.
+
Problem 1 – Take a Dive
+∇ score←dd DiveScore scores
+ :If 7=≢scores
+ scores←scores[¯2↓2↓⍋scores]
+ :ElseIf 5=≢scores
+ scores←scores[¯1↓1↓⍋scores]
+ :Else
+ scores←scores
+ :EndIf
+ score←2(⍎⍕)dd×+/scores
+∇Problem 2 – Another Step in the Proper Direction
+∇ steps←{p}Steps fromTo;segments;width
+ width←|-/fromTo
+ :If 0=⎕NC'p' ⍝ No left argument: same as Problem 5 of Phase I
+ segments←0,⍳width
+ :ElseIf p<0 ⍝ -⌊p is the number of equally-sized steps to take
+ segments←(-⌊p){0,⍵×⍺÷⍨⍳⍺}width
+ :ElseIf p>0 ⍝ p is the step size
+ segments←p{⍵⌊⍺×0,⍳⌈⍵÷⍺}width
+ :ElseIf p=0 ⍝ As if we took zero step
+ segments←0
+ :EndIf
+ ⍝ Take into account the start point and the direction.
+ steps←fromTo{(⊃⍺)+(-×-/⍺)×⍵}segments
+∇Problem 3 – Past Tasks Blast
+∇ urls←PastTasks url;r;paths
+ r←HttpCommand.Get url
+ paths←('[a-zA-Z0-9_/]+\.pdf'⎕S'&')r.Data
+ urls←('https://www.dyalog.com/'∘,)¨paths
+∇Problem 4 – Bioinformatics
+⍝ Test if a DNA string is a reverse palindrome.
+isrevp←{⍵≡⌽'TAGC'['ATCG'⍳⍵]}
+
+⍝ Generate all subarrays (position, length) pairs, for
+⍝ 4 ≤ length ≤ 12.
+subarrays←{⊃,/(⍳⍵),¨¨3↓¨⍳¨12⌊1+⍵-⍳⍵}
+
+∇ r←revp dna;positions
+ positions←subarrays⍴dna
+ ⍝ Filter subarrays which are reverse palindromes.
+ r←↑({isrevp dna[¯1+⍵[1]+⍳⍵[2]]}¨positions)/positions
+∇Problem 5 – Future and Present Value
+⍝ First solution: ((1+⊢)⊥⊣) computes the total return
+⍝ for a vector of amounts ⍺ and a vector of rates
+⍝ ⍵. It is applied to every prefix subarray of amounts
+⍝ and rates to get all intermediate values. However,
+⍝ this has quadratic complexity.
+⍝ rr←(,\⊣)((1+⊢)⊥⊣)¨(,\⊢)
+
+⍝ Second solution: We want to be able to use the
+⍝ recurrence relation (recur) and scan through the
+⍝ vectors of amounts and rates, accumulating the total
+⍝ value at every time step. However, APL evaluation is
+⍝ right-associative, so a simple Scan
+⍝ (recur\amounts,¨values) would not give the correct
+⍝ result, since recur is not associative and we need
+⍝ to evaluate it left-to-right. (In any case, in this
+⍝ case, Scan would have quadratic complexity, so would
+⍝ not bring any benefit over the previous solution.)
+⍝ What we need is something akin to Haskell's scanl
+⍝ function, which would evaluate left to right in O(n)
+⍝ time. This is what we do here, accumulating values
+⍝ from left to right. (This is inspired from
+⍝ dfns.ascan, although heavily simplified.)
+rr←{recur←{⍵[1]+⍺×1+⍵[2]} ⋄ 1↓⌽⊃{(⊂(⊃⍵)recur⍺),⍵}/⌽⍺,¨⍵}Problem 6 – Merge
+∇ val←ns getval var
+ :If ''≡var ⍝ literal '@'
+ val←'@'
+ :ElseIf (⊂var)∊ns.⎕NL ¯2
+ val←⍕ns⍎var
+ :Else
+ val←'???'
+ :EndIf
+∇∇ text←templateFile Merge jsonFile;template;ns
+ template←⊃⎕NGET templateFile 1
+ ns←⎕JSON⊃⎕NGET jsonFile
+ ⍝ We use a simple regex search and replace on the
+ ⍝ template.
+ text←↑('@[a-zA-Z]*@'⎕R{ns getval ¯1↓1↓⍵.Match})template
+∇Problem 7 – UPC
+ +⍝ Left and right representations of digits. Decoding
+⍝ the binary representation from decimal is more
+⍝ compact than writing everything explicitly.
+lrepr←⍉(7⍴2)⊤13 25 19 61 35 49 47 59 55 11
+rrepr←~¨lrepr∇ bits←WriteUPC digits;left;right
+ :If (11=≢digits)∧∧/digits∊0,⍳9
+ left←,lrepr[1+6↑digits;]
+ right←,rrepr[1+6↓digits,CheckDigit digits;]
+ bits←1 0 1,left,0 1 0 1 0,right,1 0 1
+ :Else
+ bits←¯1
+ :EndIf
+∇∇ digits←ReadUPC bits
+ :If 95≠⍴bits ⍝ incorrect number of bits
+ digits←¯1
+ :Else
+ ⍝ Test if the barcode was scanned right-to-left.
+ :If 0=2|+/bits[3+⍳7]
+ bits←⌽bits
+ :EndIf
+ digits←({¯1+lrepr⍳⍵}¨(7/⍳6)⊆42↑3↓bits),{¯1+rrepr⍳⍵}¨(7/⍳6)⊆¯42↑¯3↓bits
+ :If ~∧/digits∊0,⍳9 ⍝ incorrect parity
+ digits←¯1
+ :ElseIf (⊃⌽digits)≠CheckDigit ¯1↓digits ⍝ incorrect check digit
+ digits←¯1
+ :EndIf
+ :EndIf
+∇Problem 8 – Balancing the Scales
+∇ parts←Balance nums;subsets;partitions
+ ⍝ This is a brute force solution, running in
+ ⍝ exponential time. We generate all the possible
+ ⍝ partitions, filter out those which are not
+ ⍝ balanced, and return the first matching one. There
+ ⍝ are more advanced approach running in
+ ⍝ pseudo-polynomial time (based on dynamic
+ ⍝ programming, see the "Partition problem" Wikipedia
+ ⍝ page), but they are not warranted here, as the
+ ⍝ input size remains fairly small.
+
+ ⍝ Generate all partitions of a vector of a given
+ ⍝ size, as binary mask vectors.
+ subsets←{1↓2⊥⍣¯1⍳2*⍵}
+ ⍝ Keep only the subsets whose sum is exactly
+ ⍝ (+/nums)÷2.
+ partitions←nums{((2÷⍨+/⍺)=⍺+.×⍵)/⍵}subsets⍴nums
+ :If 0=≢,partitions
+ ⍝ If no partition satisfy the above
+ ⍝ criterion, we return ⍬.
+ parts←⍬
+ :Else
+ ⍝ Otherwise, we return the first possible
+ ⍝ partition.
+ parts←nums{((⊂,(⊂~))⊃↓⍉⍵)/¨2⍴⊂⍺}partitions
+ :EndIf
+∇Problem 9 – Upwardly Mobile
+∇ weights←Weights filename;mobile;branches;mat
+ ⍝ Put your code and comments below here
+
+ ⍝ Parse the mobile input file.
+ mobile←↑⊃⎕NGET filename 1
+ branches←⍸mobile∊'┌┴┐'
+ ⍝ TODO: Build the matrix of coefficients mat.
+
+ ⍝ Solve the system of equations (arbitrarily setting
+ ⍝ the first variable at 1 because the system is
+ ⍝ overdetermined), then multiply the coefficients by
+ ⍝ their least common multiple to get the smallest
+ ⍝ integer weights.
+ weights←((1∘,)×(∧/÷))mat[;1]⌹1↓[2]mat
+∇Table of Contents
-
+
- Introduction +
- 1. Let’s Split! +
- 2. Character Building +
- 3. Excel-lent Columns +
- 4. Take a Leap +
- 5. Stepping in the Proper Direction +
- 6. Please Move to the Front +
- 7. See You in a Bit +
- 8. Zigzag Numbers +
- 9. Rise and Fall +
- 10. Stacking It Up +
Introduction
+I’ve always been quite fond of APL and its “array-oriented” approach of programmingSee my previous post on simulating the Ising model with APL. It also contains more background on APL.
+
+. Every year, Dyalog (the company behind probably the most popular APL implementation) organises a competition with various challenges in APL.
The Dyalog APL Problem Solving Competition consists of two phases:
+-
+
- Phase I consists of 10 short puzzles (similar to what one can find on Project Euler or similar), that can be solved by a one-line APL function. +
- Phase II is a collection of larger problems, that may require longer solutions and a larger context (e.g. reading and writing to files), often in a more applied setting. Problems are often inspired by existing domains, such as AI, bioinformatics, and so on. +
In 2018, I participated in the competition, entering only Phase ISince I was a student at the time, I was eligible for a prize, and I won $100 for a 10-line submission, which is quite good!
+
+ (my solutions are on GitHub). This year, I entered in both phases. I explain my solutions to Phase I in this post. Another post will contain annotated solutions for Phase II problems.
The full code for my submission is on GitHub at dlozeve/apl-competition-2020, but everything is reproduced in this post.
+1. Let’s Split!
+++Write a function that, given a right argument
+Ywhich is a scalar or a non-empty vector and a left argumentXwhich is a single non-zero integer so that its absolute value is less or equal to≢Y, splitsYinto a vector of two vectors according toX, as follows:If
+X>0, the first vector contains the firstXelements ofYand the second vector contains the remaining elements.If
+X<0, the second vector contains the last|Xelements ofYand the first vector contains the remaining elements.
Solution: (0>⊣)⌽((⊂↑),(⊂↓))
There are three nested trains hereTrains are nice to read (even if they are easy to abuse), and generally make for shorter dfns, which is better for Phase I.
+
+. The first one, ((⊂↑),(⊂↓)), uses the two functions Take (↑) and Drop (↓) to build a nested array consisting of the two outputs we need. (Take and Drop already have the behaviour needed regarding negative arguments.) However, if the left argument is positive, the two arrays will not be in the correct order. So we need a way to reverse them if X<0.
The second train (0>⊣) will return 1 if its left argument is positive. From this, we can use Rotate (⌽) to correctly order the nested array, in the last train.
2. Character Building
+++UTF-8 encodes Unicode characters using 1-4 integers for each character. Dyalog APL includes a system function,
+⎕UCS, that can convert characters into integers and integers into characters. The expression'UTF-8'∘⎕UCSconverts between characters and UTF-8.Consider the following:
++'UTF-8'∘⎕UCS 'D¥⍺⌊○9' +68 194 165 226 141 186 226 140 138 226 151 139 57 + 'UTF-8'∘⎕UCS 68 194 165 226 141 186 226 140 138 226 151 139 57 +D¥⍺⌊○9How many integers does each character use?
++'UTF-8'∘⎕UCS¨ 'D¥⍺⌊○9' ⍝ using ]Boxing on +┌──┬───────┬───────────┬───────────┬───────────┬──┐ +│68│194 165│226 141 186│226 140 138│226 151 139│57│ +└──┴───────┴───────────┴───────────┴───────────┴──┘The rule is that an integer in the range 128 to 191 (inclusive) continues the character of the previous integer (which may itself be a continuation). With that in mind, write a function that, given a right argument which is a simple integer vector representing valid UTF-8 text, encloses each sequence of integers that represent a single character, like the result of
+'UTF-8'∘⎕UCS¨'UTF-8'∘⎕UCSbut does not use any system functions (names beginning with⎕)
Solution: {(~⍵∊127+⍳64)⊂⍵}
First, we build a binary array from the string, encoding each continuation character as 0, and all the others as 1. Next, we can use this binary array with Partitioned Enclose (⊂) to return the correct output.
3. Excel-lent Columns
+++A Microsoft Excel spreadsheet numbers its rows counting up from 1. However, Excel’s columns are labelled alphabetically — beginning with A–Z, then AA–AZ, BA–BZ, up to ZA–ZZ, then AAA–AAZ and so on.
+Write a function that, given a right argument which is a character scalar or non-empty vector representing a valid character Excel column identifier between A and XFD, returns the corresponding column number
+
Solution: 26⊥⎕A∘⍳
We use the alphabet ⎕A and Index Of (⍳) to compute the index in the alphabet of every character. As a train, this can be done by (⎕A∘⍳). We then obtain an array of numbers, each representing a letter from 1 to 26. The Decode (⊥) function can then turn this base-26 number into the expected result.
4. Take a Leap
+++Write a function that, given a right argument which is an integer array of year numbers greater than or equal to 1752 and less than 4000, returns a result of the same shape as the right argument where 1 indicates that the corresponding year is a leap year (0 otherwise).
+A leap year algorithm can be found here.
+
Solution: 1 3∊⍨(0+.=400 100 4∘.|⊢)
According to the algorithm, a year is a leap year in two situations:
+-
+
- if it is divisible by 4, but not 100 (and therefore not 400), +
- if it is divisible by 400 (and therefore 4 and 100 as well). +
The train (400 100 4∘.|⊢) will test if each year in the right argument is divisible by 400, 100, and 4, using an Outer Product. We then use an Inner Product to count how many times each year is divisible by one of these numbers. If the count is 1 or 3, it is a leap year. Note that we use Commute (⍨) to keep the dfn as a train, and to preserve the natural right-to-left reading of the algorithm.
5. Stepping in the Proper Direction
+++Write a function that, given a right argument of 2 integers, returns a vector of the integers from the first element of the right argument to the second, inclusively.
+
Solution: {(⊃⍵)+(-×-/⍵)×0,⍳|-/⍵}
First, we have to compute the range of the output, which is the absolute value of the difference between the two integers |-/⍵. From this, we compute the actual sequence, including zeroIf we had ⎕IO←0, we could have written ⍳|1+-/⍵, but this is the same number of characters.
+
+: 0,⍳|-/⍵.
This sequence will always be nondecreasing, but we have to make it decreasing if needed, so we multiply it by the opposite of the sign of -/⍵. Finally, we just have to start the sequence at the first element of ⍵.
6. Please Move to the Front
+++Write a function that, given a right argument which is an integer vector and a left argument which is an integer scalar, reorders the right argument so any elements equal to the left argument come first while all other elements keep their order.
+
Solution: {⍵[⍋⍺≠⍵]}
⍺≠⍵ will return a binary vector marking as 0 all elements equal to the left argument. Using this index to sort in the usual way with Grade Up will return the expected result.
7. See You in a Bit
+++A common technique for encoding a set of on/off states is to use a value of \(2^n\) for the state in position \(n\) (origin 0), 1 if the state is “on” or 0 for “off” and then add the values. Dyalog APL’s component file permission codes are an example of this. For example, if you wanted to grant permissions for read (access code 1), append (access code 8) and rename (access code 128) then the resulting code would be 137 because that’s 1 + 8 + 128.
+Write a function that, given a non-negative right argument which is an integer scalar representing the encoded state and a left argument which is an integer scalar representing the encoded state settings that you want to query, returns 1 if all of the codes in the left argument are found in the right argument (0 otherwise).
+
Solution: {f←⍸∘⌽(2∘⊥⍣¯1)⋄∧/(f⍺)∊f⍵}
The difficult part is to find the set of states for an integer. We need a function that will return 1 8 128 (or an equivalent representation) for an input of 137. To do this, we need the base-2 representations of \(137 = 1 + 8 + 128 = 2^0 + 2^3 + 2^7 =
+10010001_2\). The function (2∘⊥⍣¯1) will return the base-2 representation of its argument, and by reversing and finding where the non-zero elements are, we find the correct exponents (1 3 7 in this case). That is what the function f does.
Next, we just need to check that all elements of f⍺ are also in f⍵.
8. Zigzag Numbers
+++A zigzag number is an integer in which the difference in magnitude of each pair of consecutive digits alternates from positive to negative or negative to positive.
+Write a function that takes a single integer greater than or equal to 100 and less than 1015 as its right argument and returns a 1 if the integer is a zigzag number, 0 otherwise.
+
Solution: ∧/2=∘|2-/∘×2-/(10∘⊥⍣¯1)
First, we decompose a number into an array of digits, using (10∘⊥⍣¯1) (Decode (⊥) in base 10). Then, we Reduce N Wise to compute the difference between each pair of digits, take the sign, and ensure that the signs are indeed alternating.
9. Rise and Fall
+++Write a function that, given a right argument which is an integer scalar or vector, returns a 1 if the values of the right argument conform to the following pattern (0 otherwise):
++
+- The elements increase or stay the same until the “apex” (the highest value) is reached
+- After the apex, any remaining values decrease or remain the same
+
Solution: {∧/(⍳∘≢≡⍋)¨(⊂((⊢⍳⌈/)↑⊢),⍵),⊂⌽((⊢⍳⌈/)↓⊢),⍵}
How do we approach this? First we have to split the vector at the “apex”. The train (⊢⍳⌈/) will return the index of (⍳) the maximum element.
Combined with Take (↑) and Drop (↓), we build a two-element vector containing both parts, in ascending order (we Reverse (⌽) one of them). Note that we have to Ravel (,) the argument to avoid rank errors in Index Of.
Next, (⍳∘≢≡⍋) on each of the two vectors will test if they are non-decreasing (i.e. if the ranks of all the elements correspond to a simple range from 1 to the size of the vector).
10. Stacking It Up
+++Write a function that takes as its right argument a vector of simple arrays of rank 2 or less (scalar, vector, or matrix). Each simple array will consist of either non-negative integers or printable ASCII characters. The function must return a simple character array that displays identically to what
+{⎕←⍵}¨displays when applied to the right argument.
Solution: {↑⊃,/↓¨⍕¨⍵}
The first step is to Format (⍕) everything to get strings. A lot of trial-and-error is always necessary when dealing with nested arrays, and this being about formatting exacerbates the problem.
+
+ The next step would be to “stack everything vertically”, so we will need Mix (↑) at some point. However, if we do it immediately we don’t get the correct result:
Mix is padding with spaces both horizontally (necessary as we want the output to be a simple array of characters) and vertically (not what we want). We will have to decompose everything line by line, and then mix all the lines together. This is exactly what SplitSplit is the dual of Mix.
+
+ (↓) does:
{↓¨⍕¨⍵}(3 3⍴⍳9)(↑'Adam' 'Michael')(⍳10) '*'(5 5⍴⍳25)
+┌───────────────────┬─────────────────┬──────────────────────┬─┬───────────────
+│┌─────┬─────┬─────┐│┌───────┬───────┐│┌────────────────────┐│*│┌──────────────
+││1 2 3│4 5 6│7 8 9│││Adam │Michael│││1 2 3 4 5 6 7 8 9 10││ ││ 1 2 3 4 5
+│└─────┴─────┴─────┘│└───────┴───────┘│└────────────────────┘│ │└──────────────
+└───────────────────┴─────────────────┴──────────────────────┴─┴───────────────
+
+ ─────────────────────────────────────────────────────────────┐
+ ┬──────────────┬──────────────┬──────────────┬──────────────┐│
+ │ 6 7 8 9 10│11 12 13 14 15│16 17 18 19 20│21 22 23 24 25││
+ ┴──────────────┴──────────────┴──────────────┴──────────────┘│
+ ─────────────────────────────────────────────────────────────┘Next, we clean this up with Ravel (,) and we can Mix to obtain the final result.
Table of Contents
-The Ising model is a model used to represent magnetic dipole moments in statistical physics. Physical details are on the Wikipedia page, but what is interesting is that it follows a complex probability distribution on a lattice, where each site can take the value +1 or -1.
-
Mathematical definition
-We have a lattice \(\Lambda\) consisting of sites \(k\). For each site, there is a moment \(\sigma_k \in \{ -1, +1 \}\). \(\sigma = -(\sigma_k)_{k\in\Lambda}\) is called the configuration of the lattice.
-The total energy of the configuration is given by the Hamiltonian \[ -H(\sigma) = -\sum_{i\sim j} J_{ij}\, \sigma_i\, \sigma_j, -\] where \(i\sim j\) denotes neighbours, and \(J\) is the interaction matrix.
-The configuration probability is given by: \[ -\pi_\beta(\sigma) = \frac{e^{-\beta H(\sigma)}}{Z_\beta} -\] where \(\beta = (k_B T)^{-1}\) is the inverse temperature, and \(Z_\beta\) the normalisation constant.
-For our simulation, we will use a constant interaction term \(J > 0\). If \(\sigma_i = \sigma_j\), the probability will be proportional to \(\exp(\beta J)\), otherwise it would be \(\exp(\beta J)\). Thus, adjacent spins will try to align themselves.
-Simulation
-The Ising model is generally simulated using Markov Chain Monte Carlo (MCMC), with the Metropolis-Hastings algorithm.
-The algorithm starts from a random configuration and runs as follows:
--
-
- Select a site \(i\) at random and reverse its spin: \(\sigma'_i = -\sigma_i\) -
- Compute the variation in energy (hamiltonian) \(\Delta E = H(\sigma') - H(\sigma)\) -
- If the energy is lower, accept the new configuration -
- Otherwise, draw a uniform random number \(u \in ]0,1[\) and accept the new configuration if \(u < \min(1, e^{-\beta \Delta E})\). -
Implementation
-The simulation is in Clojure, using the Quil library (a Processing library for Clojure) to display the state of the system.
-This post is “literate Clojure”, and contains core.clj. The complete project can be found on GitHub.
The application works with Quil’s functional mode, with each function taking a state and returning an updated state at each time step.
-The setup function generates the initial state, with random initial spins. It also sets the frame rate. The matrix is a single vector in row-major mode. The state also holds relevant parameters for the simulation: \(\beta\), \(J\), and the iteration step.
(defn setup [size]
- "Setup the display parameters and the initial state"
- (q/frame-rate 300)
- (q/color-mode :hsb)
- (let [matrix (vec (repeatedly (* size size) #(- (* 2 (rand-int 2)) 1)))]
- {:grid-size size
- :matrix matrix
- :beta 10
- :intensity 10
- :iteration 0}))Given a site \(i\), we reverse its spin to generate a new configuration state.
-(defn toggle-state [state i]
- "Compute the new state when we toggle a cell's value"
- (let [matrix (:matrix state)]
- (assoc state :matrix (assoc matrix i (* -1 (matrix i))))))In order to decide whether to accept this new state, we compute the difference in energy introduced by reversing site \(i\): \[ \Delta E = -J\sigma_i \sum_{j\sim i} \sigma_j. \]
-The filter some? is required to eliminate sites outside of the boundaries of the lattice.
(defn get-neighbours [state idx]
- "Return the values of a cell's neighbours"
- [(get (:matrix state) (- idx (:grid-size state)))
- (get (:matrix state) (dec idx))
- (get (:matrix state) (inc idx))
- (get (:matrix state) (+ (:grid-size state) idx))])
-
-(defn delta-e [state i]
- "Compute the energy difference introduced by a particular cell"
- (* (:intensity state) ((:matrix state) i)
- (reduce + (filter some? (get-neighbours state i)))))We also add a function to compute directly the hamiltonian for the entire configuration state. We can use it later to log its values across iterations.
-(defn hamiltonian [state]
- "Compute the Hamiltonian of a configuration state"
- (- (reduce + (for [i (range (count (:matrix state)))
- j (filter some? (get-neighbours state i))]
- (* (:intensity state) ((:matrix state) i) j)))))Finally, we put everything together in the update-state function, which will decide whether to accept or reject the new configuration.
(defn update-state [state]
- "Accept or reject a new state based on energy
- difference (Metropolis-Hastings)"
- (let [i (rand-int (count (:matrix state)))
- new-state (toggle-state state i)
- alpha (q/exp (- (* (:beta state) (delta-e state i))))]
- ;;(println (hamiltonian new-state))
- (update (if (< (rand) alpha) new-state state)
- :iteration inc)))The last thing to do is to draw the new configuration:
-(defn draw-state [state]
- "Draw a configuration state as a grid"
- (q/background 255)
- (let [cell-size (quot (q/width) (:grid-size state))]
- (doseq [[i v] (map-indexed vector (:matrix state))]
- (let [x (* cell-size (rem i (:grid-size state)))
- y (* cell-size (quot i (:grid-size state)))]
- (q/no-stroke)
- (q/fill
- (if (= 1 v) 0 255))
- (q/rect x y cell-size cell-size))))
- ;;(when (zero? (mod (:iteration state) 50)) (q/save-frame "img/ising-######.jpg"))
- )And to reset the simulation when the user clicks anywhere on the screen:
-(defn mouse-clicked [state event]
- "When the mouse is clicked, reset the configuration to a random one"
- (setup 100))(q/defsketch ising-model
- :title "Ising model"
- :size [300 300]
- :setup #(setup 100)
- :update update-state
- :draw draw-state
- :mouse-clicked mouse-clicked
- :features [:keep-on-top :no-bind-output]
- :middleware [m/fun-mode])Conclusion
-The Ising model is a really easy (and common) example use of MCMC and Metropolis-Hastings. It allows to easily and intuitively understand how the algorithm works, and to make nice visualizations!
-Table of Contents
-
-
- What is an L-system? -
- Implementation details -
- Common file format for L-systems -
- Variations on L-systems -
- Usage notes -
- References -
L-systems are a formal way to make interesting visualisations. You can use them to model a wide variety of objects: space-filling curves, fractals, biological systems, tilings, etc.
-See the Github repo: https://github.com/dlozeve/lsystems
-What is an L-system?
-A few examples to get started
-
Definition
-An L-system is a set of rewriting rules generating sequences of symbols. Formally, an L-system is a triplet of:
--
-
- an alphabet \(V\) (an arbitrary set of symbols) -
- an axiom \(\omega\), which is a non-empty word of the alphabet (\(\omega \in V^+\)) -
- a set of rewriting rules (or productions) \(P\), each mapping a symbol to a word: \(P \subset V \times V^*\). Symbols not present in \(P\) are assumed to be mapped to themselves. -
During an iteration, the algorithm takes each symbol in the current word and replaces it by the value in its rewriting rule. Not that the output of the rewriting rule can be absolutely anything in \(V^*\), including the empty word! (So yes, you can generate symbols just to delete them afterwards.)
-At this point, an L-system is nothing more than a way to generate very long strings of characters. In order to get something useful out of this, we have to give them meaning.
-Drawing instructions and representation
-Our objective is to draw the output of the L-system in order to visually inspect the output. The most common way is to interpret the output as a sequence of instruction for a LOGO-like drawing turtle. For instance, a simple alphabet consisting only in the symbols \(F\), \(+\), and \(-\) could represent the instructions “move forward”, “turn right by 90°”, and “turn left by 90°” respectively.
-Thus, we add new components to our definition of L-systems:
--
-
- a set of instructions, \(I\). These are limited by the capabilities of our imagined turtle, so we can assume that they are the same for every L-system we will consider:
-
-
-
Forwardmakes the turtle draw a straight segment.
-TurnLeftandTurnRightmakes the turtle turn on itself by a given angle.
-PushandPopallow the turtle to store and retrieve its position on a stack. This will allow for branching in the turtle’s path.
-Stay, which orders the turtle to do nothing.
-
- - a distance \(d \in \mathbb{R_+}\), i.e. how long should each forward segment should be. -
- an angle \(\theta\) used for rotation. -
- a set of representation rules \(R \subset V \times I\). As before, they will match a symbol to an instruction. Symbols not matched by any rule will be associated to
Stay.
-
Finally, our complete L-system, representable by a turtle with capabilities \(I\), can be defined as \[ L = (V, \omega, P, d, \theta, -R). \]
-One could argue that the representation is not part of the L-system, and that the same L-system could be represented differently by changing the representation rules. However, in our setting, we won’t observe the L-system other than by displaying it, so we might as well consider that two systems differing only by their representation rules are different systems altogether.
-Implementation details
-The LSystem data type
-The mathematical definition above translate almost immediately in a Haskell data type:
--- | L-system data type
-data LSystem a = LSystem
- { name :: String
- , alphabet :: [a] -- ^ variables and constants used by the system
- , axiom :: [a] -- ^ initial state of the system
- , rules :: [(a, [a])] -- ^ production rules defining how each
- -- variable can be replaced by a sequence of
- -- variables and constants
- , angle :: Float -- ^ angle used for the representation
- , distance :: Float -- ^ distance of each segment in the representation
- , representation :: [(a, Instruction)] -- ^ representation rules
- -- defining how each variable
- -- and constant should be
- -- represented
- } deriving (Eq, Show, Generic)Here, a is the type of the literal in the alphabet. For all practical purposes, it will almost always be Char.
Instruction is just a sum type over all possible instructions listed above.
Iterating and representing
-From here, generating L-systems and iterating is straightforward. We iterate recursively by looking up each symbol in rules and replacing it by its expansion. We then transform the result to a list of Instruction.
Drawing
-The only remaining thing is to implement the virtual turtle which will actually execute the instructions. It goes through the list of instructions, building a sequence of points and maintaining an internal state (position, angle, stack). The stack is used when Push and Pop operations are met. In this case, the turtle builds a separate line starting from its current position.
The final output is a set of lines, each being a simple sequence of points. All relevant data types are provided by the Gloss library, along with the function that can display the resulting Picture.
Common file format for L-systems
-In order to define new L-systems quickly and easily, it is necessary to encode them in some form. We chose to represent them as JSON values.
-Here is an example for the Gosper curve:
-{
- "name": "gosper",
- "alphabet": "AB+-",
- "axiom": "A",
- "rules": [
- ["A", "A-B--B+A++AA+B-"],
- ["B", "+A-BB--B-A++A+B"]
- ],
- "angle": 60.0,
- "distance": 10.0,
- "representation": [
- ["A", "Forward"],
- ["B", "Forward"],
- ["+", "TurnRight"],
- ["-", "TurnLeft"]
- ]
-}Using this format, it is easy to define new L-systems (along with how they should be represented). This is translated nearly automatically to the LSystem data type using Aeson.
Variations on L-systems
-We can widen the possibilities of L-systems in various ways. L-systems are in effect deterministic context-free grammars.
-By allowing multiple rewriting rules for each symbol with probabilities, we can extend the model to probabilistic context-free grammars.
-We can also have replacement rules not for a single symbol, but for a subsequence of them, thus effectively taking into account their neighbours (context-sensitive grammars). This seems very close to 1D cellular automata.
-Finally, L-systems could also have a 3D representation (for instance space-filling curves in 3 dimensions).
-Usage notes
--
-
- Clone the repository:
git clone [[https://github.com/dlozeve/lsystems]]
- - Build:
stack build
- - Execute
stack exec lsystems-exe -- examples/penroseP3.jsonto see the list of options
- - (Optional) Run tests and build documentation:
stack test --haddock
-
Usage: stack exec lsystems-exe -- --help
lsystems -- Generate L-systems
-
-Usage: lsystems-exe FILENAME [-n|--iterations N] [-c|--color R,G,B]
- [-w|--white-background]
- Generate and draw an L-system
-
-Available options:
- FILENAME JSON file specifying an L-system
- -n,--iterations N Number of iterations (default: 5)
- -c,--color R,G,B Foreground color RGBA
- (0-255) (default: RGBA 1.0 1.0 1.0 1.0)
- -w,--white-background Use a white background
- -h,--help Show this help text
-Apart from the selection of the input JSON file, you can adjust the number of iterations and the colors.
-stack exec lsystems-exe -- examples/levyC.json -n 12 -c 0,255,255
References
--
-
- Prusinkiewicz, Przemyslaw; Lindenmayer, Aristid (1990). The Algorithmic Beauty of Plants. Springer-Verlag. ISBN 978-0-387-97297-8. http://algorithmicbotany.org/papers/#abop -
- Weisstein, Eric W. “Lindenmayer System.” From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/LindenmayerSystem.html -
- Corte, Leo. “L-systems and Penrose P3 in Inkscape.” The Brick in the Sky. https://thebrickinthesky.wordpress.com/2013/03/17/l-systems-and-penrose-p3-in-inkscape/ -
-
+
- + Dyalog APL Problem Solving Competition 2020 — Phase II - August 2, 2020 + + +
- + Dyalog APL Problem Solving Competition 2020 — Phase I - August 2, 2020 + +
- Operations Research and Optimization: where to start? - May 27, 2020 @@ -105,14 +113,6 @@ public key: RWQ6uexORp8f7USHA7nX9lFfltaCA9x6aBV06MvgiGjUt6BVf6McyD26 Ising model simulation in APL - March 5, 2018 -
- - Ising model simulation - February 5, 2018 - - -
- - Generating and representing L-systems - January 18, 2018 - -
Table of Contents
-
+
- Introduction +
- 1. Let’s Split! +
- 2. Character Building +
- 3. Excel-lent Columns +
- 4. Take a Leap +
- 5. Stepping in the Proper Direction +
- 6. Please Move to the Front +
- 7. See You in a Bit +
- 8. Zigzag Numbers +
- 9. Rise and Fall +
- 10. Stacking It Up +
Introduction
+I’ve always been quite fond of APL and its “array-oriented” approach of programmingSee my previous post on simulating the Ising model with APL. It also contains more background on APL.
+
+. Every year, Dyalog (the company behind probably the most popular APL implementation) organises a competition with various challenges in APL.
The Dyalog APL Problem Solving Competition consists of two phases:
+-
+
- Phase I consists of 10 short puzzles (similar to what one can find on Project Euler or similar), that can be solved by a one-line APL function. +
- Phase II is a collection of larger problems, that may require longer solutions and a larger context (e.g. reading and writing to files), often in a more applied setting. Problems are often inspired by existing domains, such as AI, bioinformatics, and so on. +
In 2018, I participated in the competition, entering only Phase ISince I was a student at the time, I was eligible for a prize, and I won $100 for a 10-line submission, which is quite good!
+
+ (my solutions are on GitHub). This year, I entered in both phases. I explain my solutions to Phase I in this post. Another post will contain annotated solutions for Phase II problems.
The full code for my submission is on GitHub at dlozeve/apl-competition-2020, but everything is reproduced in this post.
+1. Let’s Split!
+++Write a function that, given a right argument
+Ywhich is a scalar or a non-empty vector and a left argumentXwhich is a single non-zero integer so that its absolute value is less or equal to≢Y, splitsYinto a vector of two vectors according toX, as follows:If
+X>0, the first vector contains the firstXelements ofYand the second vector contains the remaining elements.If
+X<0, the second vector contains the last|Xelements ofYand the first vector contains the remaining elements.
Solution: (0>⊣)⌽((⊂↑),(⊂↓))
There are three nested trains hereTrains are nice to read (even if they are easy to abuse), and generally make for shorter dfns, which is better for Phase I.
+
+. The first one, ((⊂↑),(⊂↓)), uses the two functions Take (↑) and Drop (↓) to build a nested array consisting of the two outputs we need. (Take and Drop already have the behaviour needed regarding negative arguments.) However, if the left argument is positive, the two arrays will not be in the correct order. So we need a way to reverse them if X<0.
The second train (0>⊣) will return 1 if its left argument is positive. From this, we can use Rotate (⌽) to correctly order the nested array, in the last train.
2. Character Building
+++UTF-8 encodes Unicode characters using 1-4 integers for each character. Dyalog APL includes a system function,
+⎕UCS, that can convert characters into integers and integers into characters. The expression'UTF-8'∘⎕UCSconverts between characters and UTF-8.Consider the following:
++'UTF-8'∘⎕UCS 'D¥⍺⌊○9' +68 194 165 226 141 186 226 140 138 226 151 139 57 + 'UTF-8'∘⎕UCS 68 194 165 226 141 186 226 140 138 226 151 139 57 +D¥⍺⌊○9How many integers does each character use?
++'UTF-8'∘⎕UCS¨ 'D¥⍺⌊○9' ⍝ using ]Boxing on +┌──┬───────┬───────────┬───────────┬───────────┬──┐ +│68│194 165│226 141 186│226 140 138│226 151 139│57│ +└──┴───────┴───────────┴───────────┴───────────┴──┘The rule is that an integer in the range 128 to 191 (inclusive) continues the character of the previous integer (which may itself be a continuation). With that in mind, write a function that, given a right argument which is a simple integer vector representing valid UTF-8 text, encloses each sequence of integers that represent a single character, like the result of
+'UTF-8'∘⎕UCS¨'UTF-8'∘⎕UCSbut does not use any system functions (names beginning with⎕)
Solution: {(~⍵∊127+⍳64)⊂⍵}
First, we build a binary array from the string, encoding each continuation character as 0, and all the others as 1. Next, we can use this binary array with Partitioned Enclose (⊂) to return the correct output.
3. Excel-lent Columns
+++A Microsoft Excel spreadsheet numbers its rows counting up from 1. However, Excel’s columns are labelled alphabetically — beginning with A–Z, then AA–AZ, BA–BZ, up to ZA–ZZ, then AAA–AAZ and so on.
+Write a function that, given a right argument which is a character scalar or non-empty vector representing a valid character Excel column identifier between A and XFD, returns the corresponding column number
+
Solution: 26⊥⎕A∘⍳
We use the alphabet ⎕A and Index Of (⍳) to compute the index in the alphabet of every character. As a train, this can be done by (⎕A∘⍳). We then obtain an array of numbers, each representing a letter from 1 to 26. The Decode (⊥) function can then turn this base-26 number into the expected result.
4. Take a Leap
+++Write a function that, given a right argument which is an integer array of year numbers greater than or equal to 1752 and less than 4000, returns a result of the same shape as the right argument where 1 indicates that the corresponding year is a leap year (0 otherwise).
+A leap year algorithm can be found here.
+
Solution: 1 3∊⍨(0+.=400 100 4∘.|⊢)
According to the algorithm, a year is a leap year in two situations:
+-
+
- if it is divisible by 4, but not 100 (and therefore not 400), +
- if it is divisible by 400 (and therefore 4 and 100 as well). +
The train (400 100 4∘.|⊢) will test if each year in the right argument is divisible by 400, 100, and 4, using an Outer Product. We then use an Inner Product to count how many times each year is divisible by one of these numbers. If the count is 1 or 3, it is a leap year. Note that we use Commute (⍨) to keep the dfn as a train, and to preserve the natural right-to-left reading of the algorithm.
5. Stepping in the Proper Direction
+++Write a function that, given a right argument of 2 integers, returns a vector of the integers from the first element of the right argument to the second, inclusively.
+
Solution: {(⊃⍵)+(-×-/⍵)×0,⍳|-/⍵}
First, we have to compute the range of the output, which is the absolute value of the difference between the two integers |-/⍵. From this, we compute the actual sequence, including zeroIf we had ⎕IO←0, we could have written ⍳|1+-/⍵, but this is the same number of characters.
+
+: 0,⍳|-/⍵.
This sequence will always be nondecreasing, but we have to make it decreasing if needed, so we multiply it by the opposite of the sign of -/⍵. Finally, we just have to start the sequence at the first element of ⍵.
6. Please Move to the Front
+++Write a function that, given a right argument which is an integer vector and a left argument which is an integer scalar, reorders the right argument so any elements equal to the left argument come first while all other elements keep their order.
+
Solution: {⍵[⍋⍺≠⍵]}
⍺≠⍵ will return a binary vector marking as 0 all elements equal to the left argument. Using this index to sort in the usual way with Grade Up will return the expected result.
7. See You in a Bit
+++A common technique for encoding a set of on/off states is to use a value of \(2^n\) for the state in position \(n\) (origin 0), 1 if the state is “on” or 0 for “off” and then add the values. Dyalog APL’s component file permission codes are an example of this. For example, if you wanted to grant permissions for read (access code 1), append (access code 8) and rename (access code 128) then the resulting code would be 137 because that’s 1 + 8 + 128.
+Write a function that, given a non-negative right argument which is an integer scalar representing the encoded state and a left argument which is an integer scalar representing the encoded state settings that you want to query, returns 1 if all of the codes in the left argument are found in the right argument (0 otherwise).
+
Solution: {f←⍸∘⌽(2∘⊥⍣¯1)⋄∧/(f⍺)∊f⍵}
The difficult part is to find the set of states for an integer. We need a function that will return 1 8 128 (or an equivalent representation) for an input of 137. To do this, we need the base-2 representations of \(137 = 1 + 8 + 128 = 2^0 + 2^3 + 2^7 =
+10010001_2\). The function (2∘⊥⍣¯1) will return the base-2 representation of its argument, and by reversing and finding where the non-zero elements are, we find the correct exponents (1 3 7 in this case). That is what the function f does.
Next, we just need to check that all elements of f⍺ are also in f⍵.
8. Zigzag Numbers
+++A zigzag number is an integer in which the difference in magnitude of each pair of consecutive digits alternates from positive to negative or negative to positive.
+Write a function that takes a single integer greater than or equal to 100 and less than 1015 as its right argument and returns a 1 if the integer is a zigzag number, 0 otherwise.
+
Solution: ∧/2=∘|2-/∘×2-/(10∘⊥⍣¯1)
First, we decompose a number into an array of digits, using (10∘⊥⍣¯1) (Decode (⊥) in base 10). Then, we Reduce N Wise to compute the difference between each pair of digits, take the sign, and ensure that the signs are indeed alternating.
9. Rise and Fall
+++Write a function that, given a right argument which is an integer scalar or vector, returns a 1 if the values of the right argument conform to the following pattern (0 otherwise):
++
+- The elements increase or stay the same until the “apex” (the highest value) is reached
+- After the apex, any remaining values decrease or remain the same
+
Solution: {∧/(⍳∘≢≡⍋)¨(⊂((⊢⍳⌈/)↑⊢),⍵),⊂⌽((⊢⍳⌈/)↓⊢),⍵}
How do we approach this? First we have to split the vector at the “apex”. The train (⊢⍳⌈/) will return the index of (⍳) the maximum element.
Combined with Take (↑) and Drop (↓), we build a two-element vector containing both parts, in ascending order (we Reverse (⌽) one of them). Note that we have to Ravel (,) the argument to avoid rank errors in Index Of.
Next, (⍳∘≢≡⍋) on each of the two vectors will test if they are non-decreasing (i.e. if the ranks of all the elements correspond to a simple range from 1 to the size of the vector).
10. Stacking It Up
+++Write a function that takes as its right argument a vector of simple arrays of rank 2 or less (scalar, vector, or matrix). Each simple array will consist of either non-negative integers or printable ASCII characters. The function must return a simple character array that displays identically to what
+{⎕←⍵}¨displays when applied to the right argument.
Solution: {↑⊃,/↓¨⍕¨⍵}
The first step is to Format (⍕) everything to get strings. A lot of trial-and-error is always necessary when dealing with nested arrays, and this being about formatting exacerbates the problem.
+
+ The next step would be to “stack everything vertically”, so we will need Mix (↑) at some point. However, if we do it immediately we don’t get the correct result:
Mix is padding with spaces both horizontally (necessary as we want the output to be a simple array of characters) and vertically (not what we want). We will have to decompose everything line by line, and then mix all the lines together. This is exactly what SplitSplit is the dual of Mix.
+
+ (↓) does:
{↓¨⍕¨⍵}(3 3⍴⍳9)(↑'Adam' 'Michael')(⍳10) '*'(5 5⍴⍳25)
+┌───────────────────┬─────────────────┬──────────────────────┬─┬───────────────
+│┌─────┬─────┬─────┐│┌───────┬───────┐│┌────────────────────┐│*│┌──────────────
+││1 2 3│4 5 6│7 8 9│││Adam │Michael│││1 2 3 4 5 6 7 8 9 10││ ││ 1 2 3 4 5
+│└─────┴─────┴─────┘│└───────┴───────┘│└────────────────────┘│ │└──────────────
+└───────────────────┴─────────────────┴──────────────────────┴─┴───────────────
+
+ ─────────────────────────────────────────────────────────────┐
+ ┬──────────────┬──────────────┬──────────────┬──────────────┐│
+ │ 6 7 8 9 10│11 12 13 14 15│16 17 18 19 20│21 22 23 24 25││
+ ┴──────────────┴──────────────┴──────────────┴──────────────┘│
+ ─────────────────────────────────────────────────────────────┘Next, we clean this up with Ravel (,) and we can Mix to obtain the final result.
Dyalog APL Problem Solving Competition 2020 — Phase II
+ +Annotated Solutions
+ + + + +Table of Contents
-
+
- Introduction +
- Problem 1 – Take a Dive +
- Problem 2 – Another Step in the Proper Direction +
- Problem 3 – Past Tasks Blast +
- Problem 4 – Bioinformatics +
- Problem 5 – Future and Present Value +
- Problem 6 – Merge +
- Problem 7 – UPC +
- Problem 8 – Balancing the Scales +
- Problem 9 – Upwardly Mobile +
Introduction
+After Phase I, here are my solutions to Phase II problems. The full code is included in the post, but everything is also available on GitHub.
+A PDF of the problems descriptions is available on the competition website, or directly from my GitHub repo.
+The submission guidelines gave a template where everything is defined in a Contest2020.Problems Namespace. I kept the default values for ⎕IO and ⎕ML because the problems were not particularly easier with ⎕IO←0.
++This post is still a work in progress! I will try to write explanations for every problem below.
+
Problem 1 – Take a Dive
+∇ score←dd DiveScore scores
+ :If 7=≢scores
+ scores←scores[¯2↓2↓⍋scores]
+ :ElseIf 5=≢scores
+ scores←scores[¯1↓1↓⍋scores]
+ :Else
+ scores←scores
+ :EndIf
+ score←2(⍎⍕)dd×+/scores
+∇Problem 2 – Another Step in the Proper Direction
+∇ steps←{p}Steps fromTo;segments;width
+ width←|-/fromTo
+ :If 0=⎕NC'p' ⍝ No left argument: same as Problem 5 of Phase I
+ segments←0,⍳width
+ :ElseIf p<0 ⍝ -⌊p is the number of equally-sized steps to take
+ segments←(-⌊p){0,⍵×⍺÷⍨⍳⍺}width
+ :ElseIf p>0 ⍝ p is the step size
+ segments←p{⍵⌊⍺×0,⍳⌈⍵÷⍺}width
+ :ElseIf p=0 ⍝ As if we took zero step
+ segments←0
+ :EndIf
+ ⍝ Take into account the start point and the direction.
+ steps←fromTo{(⊃⍺)+(-×-/⍺)×⍵}segments
+∇Problem 3 – Past Tasks Blast
+∇ urls←PastTasks url;r;paths
+ r←HttpCommand.Get url
+ paths←('[a-zA-Z0-9_/]+\.pdf'⎕S'&')r.Data
+ urls←('https://www.dyalog.com/'∘,)¨paths
+∇Problem 4 – Bioinformatics
+⍝ Test if a DNA string is a reverse palindrome.
+isrevp←{⍵≡⌽'TAGC'['ATCG'⍳⍵]}
+
+⍝ Generate all subarrays (position, length) pairs, for
+⍝ 4 ≤ length ≤ 12.
+subarrays←{⊃,/(⍳⍵),¨¨3↓¨⍳¨12⌊1+⍵-⍳⍵}
+
+∇ r←revp dna;positions
+ positions←subarrays⍴dna
+ ⍝ Filter subarrays which are reverse palindromes.
+ r←↑({isrevp dna[¯1+⍵[1]+⍳⍵[2]]}¨positions)/positions
+∇Problem 5 – Future and Present Value
+⍝ First solution: ((1+⊢)⊥⊣) computes the total return
+⍝ for a vector of amounts ⍺ and a vector of rates
+⍝ ⍵. It is applied to every prefix subarray of amounts
+⍝ and rates to get all intermediate values. However,
+⍝ this has quadratic complexity.
+⍝ rr←(,\⊣)((1+⊢)⊥⊣)¨(,\⊢)
+
+⍝ Second solution: We want to be able to use the
+⍝ recurrence relation (recur) and scan through the
+⍝ vectors of amounts and rates, accumulating the total
+⍝ value at every time step. However, APL evaluation is
+⍝ right-associative, so a simple Scan
+⍝ (recur\amounts,¨values) would not give the correct
+⍝ result, since recur is not associative and we need
+⍝ to evaluate it left-to-right. (In any case, in this
+⍝ case, Scan would have quadratic complexity, so would
+⍝ not bring any benefit over the previous solution.)
+⍝ What we need is something akin to Haskell's scanl
+⍝ function, which would evaluate left to right in O(n)
+⍝ time. This is what we do here, accumulating values
+⍝ from left to right. (This is inspired from
+⍝ dfns.ascan, although heavily simplified.)
+rr←{recur←{⍵[1]+⍺×1+⍵[2]} ⋄ 1↓⌽⊃{(⊂(⊃⍵)recur⍺),⍵}/⌽⍺,¨⍵}Problem 6 – Merge
+∇ val←ns getval var
+ :If ''≡var ⍝ literal '@'
+ val←'@'
+ :ElseIf (⊂var)∊ns.⎕NL ¯2
+ val←⍕ns⍎var
+ :Else
+ val←'???'
+ :EndIf
+∇∇ text←templateFile Merge jsonFile;template;ns
+ template←⊃⎕NGET templateFile 1
+ ns←⎕JSON⊃⎕NGET jsonFile
+ ⍝ We use a simple regex search and replace on the
+ ⍝ template.
+ text←↑('@[a-zA-Z]*@'⎕R{ns getval ¯1↓1↓⍵.Match})template
+∇Problem 7 – UPC
+ +⍝ Left and right representations of digits. Decoding
+⍝ the binary representation from decimal is more
+⍝ compact than writing everything explicitly.
+lrepr←⍉(7⍴2)⊤13 25 19 61 35 49 47 59 55 11
+rrepr←~¨lrepr∇ bits←WriteUPC digits;left;right
+ :If (11=≢digits)∧∧/digits∊0,⍳9
+ left←,lrepr[1+6↑digits;]
+ right←,rrepr[1+6↓digits,CheckDigit digits;]
+ bits←1 0 1,left,0 1 0 1 0,right,1 0 1
+ :Else
+ bits←¯1
+ :EndIf
+∇∇ digits←ReadUPC bits
+ :If 95≠⍴bits ⍝ incorrect number of bits
+ digits←¯1
+ :Else
+ ⍝ Test if the barcode was scanned right-to-left.
+ :If 0=2|+/bits[3+⍳7]
+ bits←⌽bits
+ :EndIf
+ digits←({¯1+lrepr⍳⍵}¨(7/⍳6)⊆42↑3↓bits),{¯1+rrepr⍳⍵}¨(7/⍳6)⊆¯42↑¯3↓bits
+ :If ~∧/digits∊0,⍳9 ⍝ incorrect parity
+ digits←¯1
+ :ElseIf (⊃⌽digits)≠CheckDigit ¯1↓digits ⍝ incorrect check digit
+ digits←¯1
+ :EndIf
+ :EndIf
+∇Problem 8 – Balancing the Scales
+∇ parts←Balance nums;subsets;partitions
+ ⍝ This is a brute force solution, running in
+ ⍝ exponential time. We generate all the possible
+ ⍝ partitions, filter out those which are not
+ ⍝ balanced, and return the first matching one. There
+ ⍝ are more advanced approach running in
+ ⍝ pseudo-polynomial time (based on dynamic
+ ⍝ programming, see the "Partition problem" Wikipedia
+ ⍝ page), but they are not warranted here, as the
+ ⍝ input size remains fairly small.
+
+ ⍝ Generate all partitions of a vector of a given
+ ⍝ size, as binary mask vectors.
+ subsets←{1↓2⊥⍣¯1⍳2*⍵}
+ ⍝ Keep only the subsets whose sum is exactly
+ ⍝ (+/nums)÷2.
+ partitions←nums{((2÷⍨+/⍺)=⍺+.×⍵)/⍵}subsets⍴nums
+ :If 0=≢,partitions
+ ⍝ If no partition satisfy the above
+ ⍝ criterion, we return ⍬.
+ parts←⍬
+ :Else
+ ⍝ Otherwise, we return the first possible
+ ⍝ partition.
+ parts←nums{((⊂,(⊂~))⊃↓⍉⍵)/¨2⍴⊂⍺}partitions
+ :EndIf
+∇Problem 9 – Upwardly Mobile
+∇ weights←Weights filename;mobile;branches;mat
+ ⍝ Put your code and comments below here
+
+ ⍝ Parse the mobile input file.
+ mobile←↑⊃⎕NGET filename 1
+ branches←⍸mobile∊'┌┴┐'
+ ⍝ TODO: Build the matrix of coefficients mat.
+
+ ⍝ Solve the system of equations (arbitrarily setting
+ ⍝ the first variable at 1 because the system is
+ ⍝ overdetermined), then multiply the coefficients by
+ ⍝ their least common multiple to get the smallest
+ ⍝ integer weights.
+ weights←((1∘,)×(∧/÷))mat[;1]⌹1↓[2]mat
+∇Table of Contents
-
+
- Introduction +
- Problem 1 – Take a Dive +
- Problem 2 – Another Step in the Proper Direction +
- Problem 3 – Past Tasks Blast +
- Problem 4 – Bioinformatics +
- Problem 5 – Future and Present Value +
- Problem 6 – Merge +
- Problem 7 – UPC +
- Problem 8 – Balancing the Scales +
- Problem 9 – Upwardly Mobile +
Introduction
+After Phase I, here are my solutions to Phase II problems. The full code is included in the post, but everything is also available on GitHub.
+A PDF of the problems descriptions is available on the competition website, or directly from my GitHub repo.
+The submission guidelines gave a template where everything is defined in a Contest2020.Problems Namespace. I kept the default values for ⎕IO and ⎕ML because the problems were not particularly easier with ⎕IO←0.
++This post is still a work in progress! I will try to write explanations for every problem below.
+
Problem 1 – Take a Dive
+∇ score←dd DiveScore scores
+ :If 7=≢scores
+ scores←scores[¯2↓2↓⍋scores]
+ :ElseIf 5=≢scores
+ scores←scores[¯1↓1↓⍋scores]
+ :Else
+ scores←scores
+ :EndIf
+ score←2(⍎⍕)dd×+/scores
+∇Problem 2 – Another Step in the Proper Direction
+∇ steps←{p}Steps fromTo;segments;width
+ width←|-/fromTo
+ :If 0=⎕NC'p' ⍝ No left argument: same as Problem 5 of Phase I
+ segments←0,⍳width
+ :ElseIf p<0 ⍝ -⌊p is the number of equally-sized steps to take
+ segments←(-⌊p){0,⍵×⍺÷⍨⍳⍺}width
+ :ElseIf p>0 ⍝ p is the step size
+ segments←p{⍵⌊⍺×0,⍳⌈⍵÷⍺}width
+ :ElseIf p=0 ⍝ As if we took zero step
+ segments←0
+ :EndIf
+ ⍝ Take into account the start point and the direction.
+ steps←fromTo{(⊃⍺)+(-×-/⍺)×⍵}segments
+∇Problem 3 – Past Tasks Blast
+∇ urls←PastTasks url;r;paths
+ r←HttpCommand.Get url
+ paths←('[a-zA-Z0-9_/]+\.pdf'⎕S'&')r.Data
+ urls←('https://www.dyalog.com/'∘,)¨paths
+∇Problem 4 – Bioinformatics
+⍝ Test if a DNA string is a reverse palindrome.
+isrevp←{⍵≡⌽'TAGC'['ATCG'⍳⍵]}
+
+⍝ Generate all subarrays (position, length) pairs, for
+⍝ 4 ≤ length ≤ 12.
+subarrays←{⊃,/(⍳⍵),¨¨3↓¨⍳¨12⌊1+⍵-⍳⍵}
+
+∇ r←revp dna;positions
+ positions←subarrays⍴dna
+ ⍝ Filter subarrays which are reverse palindromes.
+ r←↑({isrevp dna[¯1+⍵[1]+⍳⍵[2]]}¨positions)/positions
+∇Problem 5 – Future and Present Value
+⍝ First solution: ((1+⊢)⊥⊣) computes the total return
+⍝ for a vector of amounts ⍺ and a vector of rates
+⍝ ⍵. It is applied to every prefix subarray of amounts
+⍝ and rates to get all intermediate values. However,
+⍝ this has quadratic complexity.
+⍝ rr←(,\⊣)((1+⊢)⊥⊣)¨(,\⊢)
+
+⍝ Second solution: We want to be able to use the
+⍝ recurrence relation (recur) and scan through the
+⍝ vectors of amounts and rates, accumulating the total
+⍝ value at every time step. However, APL evaluation is
+⍝ right-associative, so a simple Scan
+⍝ (recur\amounts,¨values) would not give the correct
+⍝ result, since recur is not associative and we need
+⍝ to evaluate it left-to-right. (In any case, in this
+⍝ case, Scan would have quadratic complexity, so would
+⍝ not bring any benefit over the previous solution.)
+⍝ What we need is something akin to Haskell's scanl
+⍝ function, which would evaluate left to right in O(n)
+⍝ time. This is what we do here, accumulating values
+⍝ from left to right. (This is inspired from
+⍝ dfns.ascan, although heavily simplified.)
+rr←{recur←{⍵[1]+⍺×1+⍵[2]} ⋄ 1↓⌽⊃{(⊂(⊃⍵)recur⍺),⍵}/⌽⍺,¨⍵}Problem 6 – Merge
+∇ val←ns getval var
+ :If ''≡var ⍝ literal '@'
+ val←'@'
+ :ElseIf (⊂var)∊ns.⎕NL ¯2
+ val←⍕ns⍎var
+ :Else
+ val←'???'
+ :EndIf
+∇∇ text←templateFile Merge jsonFile;template;ns
+ template←⊃⎕NGET templateFile 1
+ ns←⎕JSON⊃⎕NGET jsonFile
+ ⍝ We use a simple regex search and replace on the
+ ⍝ template.
+ text←↑('@[a-zA-Z]*@'⎕R{ns getval ¯1↓1↓⍵.Match})template
+∇Problem 7 – UPC
+ +⍝ Left and right representations of digits. Decoding
+⍝ the binary representation from decimal is more
+⍝ compact than writing everything explicitly.
+lrepr←⍉(7⍴2)⊤13 25 19 61 35 49 47 59 55 11
+rrepr←~¨lrepr∇ bits←WriteUPC digits;left;right
+ :If (11=≢digits)∧∧/digits∊0,⍳9
+ left←,lrepr[1+6↑digits;]
+ right←,rrepr[1+6↓digits,CheckDigit digits;]
+ bits←1 0 1,left,0 1 0 1 0,right,1 0 1
+ :Else
+ bits←¯1
+ :EndIf
+∇∇ digits←ReadUPC bits
+ :If 95≠⍴bits ⍝ incorrect number of bits
+ digits←¯1
+ :Else
+ ⍝ Test if the barcode was scanned right-to-left.
+ :If 0=2|+/bits[3+⍳7]
+ bits←⌽bits
+ :EndIf
+ digits←({¯1+lrepr⍳⍵}¨(7/⍳6)⊆42↑3↓bits),{¯1+rrepr⍳⍵}¨(7/⍳6)⊆¯42↑¯3↓bits
+ :If ~∧/digits∊0,⍳9 ⍝ incorrect parity
+ digits←¯1
+ :ElseIf (⊃⌽digits)≠CheckDigit ¯1↓digits ⍝ incorrect check digit
+ digits←¯1
+ :EndIf
+ :EndIf
+∇Problem 8 – Balancing the Scales
+∇ parts←Balance nums;subsets;partitions
+ ⍝ This is a brute force solution, running in
+ ⍝ exponential time. We generate all the possible
+ ⍝ partitions, filter out those which are not
+ ⍝ balanced, and return the first matching one. There
+ ⍝ are more advanced approach running in
+ ⍝ pseudo-polynomial time (based on dynamic
+ ⍝ programming, see the "Partition problem" Wikipedia
+ ⍝ page), but they are not warranted here, as the
+ ⍝ input size remains fairly small.
+
+ ⍝ Generate all partitions of a vector of a given
+ ⍝ size, as binary mask vectors.
+ subsets←{1↓2⊥⍣¯1⍳2*⍵}
+ ⍝ Keep only the subsets whose sum is exactly
+ ⍝ (+/nums)÷2.
+ partitions←nums{((2÷⍨+/⍺)=⍺+.×⍵)/⍵}subsets⍴nums
+ :If 0=≢,partitions
+ ⍝ If no partition satisfy the above
+ ⍝ criterion, we return ⍬.
+ parts←⍬
+ :Else
+ ⍝ Otherwise, we return the first possible
+ ⍝ partition.
+ parts←nums{((⊂,(⊂~))⊃↓⍉⍵)/¨2⍴⊂⍺}partitions
+ :EndIf
+∇Problem 9 – Upwardly Mobile
+∇ weights←Weights filename;mobile;branches;mat
+ ⍝ Put your code and comments below here
+
+ ⍝ Parse the mobile input file.
+ mobile←↑⊃⎕NGET filename 1
+ branches←⍸mobile∊'┌┴┐'
+ ⍝ TODO: Build the matrix of coefficients mat.
+
+ ⍝ Solve the system of equations (arbitrarily setting
+ ⍝ the first variable at 1 because the system is
+ ⍝ overdetermined), then multiply the coefficients by
+ ⍝ their least common multiple to get the smallest
+ ⍝ integer weights.
+ weights←((1∘,)×(∧/÷))mat[;1]⌹1↓[2]mat
+∇Table of Contents
-
+
- Introduction +
- 1. Let’s Split! +
- 2. Character Building +
- 3. Excel-lent Columns +
- 4. Take a Leap +
- 5. Stepping in the Proper Direction +
- 6. Please Move to the Front +
- 7. See You in a Bit +
- 8. Zigzag Numbers +
- 9. Rise and Fall +
- 10. Stacking It Up +
Introduction
+I’ve always been quite fond of APL and its “array-oriented” approach of programmingSee my previous post on simulating the Ising model with APL. It also contains more background on APL.
+
+. Every year, Dyalog (the company behind probably the most popular APL implementation) organises a competition with various challenges in APL.
The Dyalog APL Problem Solving Competition consists of two phases:
+-
+
- Phase I consists of 10 short puzzles (similar to what one can find on Project Euler or similar), that can be solved by a one-line APL function. +
- Phase II is a collection of larger problems, that may require longer solutions and a larger context (e.g. reading and writing to files), often in a more applied setting. Problems are often inspired by existing domains, such as AI, bioinformatics, and so on. +
In 2018, I participated in the competition, entering only Phase ISince I was a student at the time, I was eligible for a prize, and I won $100 for a 10-line submission, which is quite good!
+
+ (my solutions are on GitHub). This year, I entered in both phases. I explain my solutions to Phase I in this post. Another post will contain annotated solutions for Phase II problems.
The full code for my submission is on GitHub at dlozeve/apl-competition-2020, but everything is reproduced in this post.
+1. Let’s Split!
+++Write a function that, given a right argument
+Ywhich is a scalar or a non-empty vector and a left argumentXwhich is a single non-zero integer so that its absolute value is less or equal to≢Y, splitsYinto a vector of two vectors according toX, as follows:If
+X>0, the first vector contains the firstXelements ofYand the second vector contains the remaining elements.If
+X<0, the second vector contains the last|Xelements ofYand the first vector contains the remaining elements.
Solution: (0>⊣)⌽((⊂↑),(⊂↓))
There are three nested trains hereTrains are nice to read (even if they are easy to abuse), and generally make for shorter dfns, which is better for Phase I.
+
+. The first one, ((⊂↑),(⊂↓)), uses the two functions Take (↑) and Drop (↓) to build a nested array consisting of the two outputs we need. (Take and Drop already have the behaviour needed regarding negative arguments.) However, if the left argument is positive, the two arrays will not be in the correct order. So we need a way to reverse them if X<0.
The second train (0>⊣) will return 1 if its left argument is positive. From this, we can use Rotate (⌽) to correctly order the nested array, in the last train.
2. Character Building
+++UTF-8 encodes Unicode characters using 1-4 integers for each character. Dyalog APL includes a system function,
+⎕UCS, that can convert characters into integers and integers into characters. The expression'UTF-8'∘⎕UCSconverts between characters and UTF-8.Consider the following:
++'UTF-8'∘⎕UCS 'D¥⍺⌊○9' +68 194 165 226 141 186 226 140 138 226 151 139 57 + 'UTF-8'∘⎕UCS 68 194 165 226 141 186 226 140 138 226 151 139 57 +D¥⍺⌊○9How many integers does each character use?
++'UTF-8'∘⎕UCS¨ 'D¥⍺⌊○9' ⍝ using ]Boxing on +┌──┬───────┬───────────┬───────────┬───────────┬──┐ +│68│194 165│226 141 186│226 140 138│226 151 139│57│ +└──┴───────┴───────────┴───────────┴───────────┴──┘The rule is that an integer in the range 128 to 191 (inclusive) continues the character of the previous integer (which may itself be a continuation). With that in mind, write a function that, given a right argument which is a simple integer vector representing valid UTF-8 text, encloses each sequence of integers that represent a single character, like the result of
+'UTF-8'∘⎕UCS¨'UTF-8'∘⎕UCSbut does not use any system functions (names beginning with⎕)
Solution: {(~⍵∊127+⍳64)⊂⍵}
First, we build a binary array from the string, encoding each continuation character as 0, and all the others as 1. Next, we can use this binary array with Partitioned Enclose (⊂) to return the correct output.
3. Excel-lent Columns
+++A Microsoft Excel spreadsheet numbers its rows counting up from 1. However, Excel’s columns are labelled alphabetically — beginning with A–Z, then AA–AZ, BA–BZ, up to ZA–ZZ, then AAA–AAZ and so on.
+Write a function that, given a right argument which is a character scalar or non-empty vector representing a valid character Excel column identifier between A and XFD, returns the corresponding column number
+
Solution: 26⊥⎕A∘⍳
We use the alphabet ⎕A and Index Of (⍳) to compute the index in the alphabet of every character. As a train, this can be done by (⎕A∘⍳). We then obtain an array of numbers, each representing a letter from 1 to 26. The Decode (⊥) function can then turn this base-26 number into the expected result.
4. Take a Leap
+++Write a function that, given a right argument which is an integer array of year numbers greater than or equal to 1752 and less than 4000, returns a result of the same shape as the right argument where 1 indicates that the corresponding year is a leap year (0 otherwise).
+A leap year algorithm can be found here.
+
Solution: 1 3∊⍨(0+.=400 100 4∘.|⊢)
According to the algorithm, a year is a leap year in two situations:
+-
+
- if it is divisible by 4, but not 100 (and therefore not 400), +
- if it is divisible by 400 (and therefore 4 and 100 as well). +
The train (400 100 4∘.|⊢) will test if each year in the right argument is divisible by 400, 100, and 4, using an Outer Product. We then use an Inner Product to count how many times each year is divisible by one of these numbers. If the count is 1 or 3, it is a leap year. Note that we use Commute (⍨) to keep the dfn as a train, and to preserve the natural right-to-left reading of the algorithm.
5. Stepping in the Proper Direction
+++Write a function that, given a right argument of 2 integers, returns a vector of the integers from the first element of the right argument to the second, inclusively.
+
Solution: {(⊃⍵)+(-×-/⍵)×0,⍳|-/⍵}
First, we have to compute the range of the output, which is the absolute value of the difference between the two integers |-/⍵. From this, we compute the actual sequence, including zeroIf we had ⎕IO←0, we could have written ⍳|1+-/⍵, but this is the same number of characters.
+
+: 0,⍳|-/⍵.
This sequence will always be nondecreasing, but we have to make it decreasing if needed, so we multiply it by the opposite of the sign of -/⍵. Finally, we just have to start the sequence at the first element of ⍵.
6. Please Move to the Front
+++Write a function that, given a right argument which is an integer vector and a left argument which is an integer scalar, reorders the right argument so any elements equal to the left argument come first while all other elements keep their order.
+
Solution: {⍵[⍋⍺≠⍵]}
⍺≠⍵ will return a binary vector marking as 0 all elements equal to the left argument. Using this index to sort in the usual way with Grade Up will return the expected result.
7. See You in a Bit
+++A common technique for encoding a set of on/off states is to use a value of \(2^n\) for the state in position \(n\) (origin 0), 1 if the state is “on” or 0 for “off” and then add the values. Dyalog APL’s component file permission codes are an example of this. For example, if you wanted to grant permissions for read (access code 1), append (access code 8) and rename (access code 128) then the resulting code would be 137 because that’s 1 + 8 + 128.
+Write a function that, given a non-negative right argument which is an integer scalar representing the encoded state and a left argument which is an integer scalar representing the encoded state settings that you want to query, returns 1 if all of the codes in the left argument are found in the right argument (0 otherwise).
+
Solution: {f←⍸∘⌽(2∘⊥⍣¯1)⋄∧/(f⍺)∊f⍵}
The difficult part is to find the set of states for an integer. We need a function that will return 1 8 128 (or an equivalent representation) for an input of 137. To do this, we need the base-2 representations of \(137 = 1 + 8 + 128 = 2^0 + 2^3 + 2^7 =
+10010001_2\). The function (2∘⊥⍣¯1) will return the base-2 representation of its argument, and by reversing and finding where the non-zero elements are, we find the correct exponents (1 3 7 in this case). That is what the function f does.
Next, we just need to check that all elements of f⍺ are also in f⍵.
8. Zigzag Numbers
+++A zigzag number is an integer in which the difference in magnitude of each pair of consecutive digits alternates from positive to negative or negative to positive.
+Write a function that takes a single integer greater than or equal to 100 and less than 1015 as its right argument and returns a 1 if the integer is a zigzag number, 0 otherwise.
+
Solution: ∧/2=∘|2-/∘×2-/(10∘⊥⍣¯1)
First, we decompose a number into an array of digits, using (10∘⊥⍣¯1) (Decode (⊥) in base 10). Then, we Reduce N Wise to compute the difference between each pair of digits, take the sign, and ensure that the signs are indeed alternating.
9. Rise and Fall
+++Write a function that, given a right argument which is an integer scalar or vector, returns a 1 if the values of the right argument conform to the following pattern (0 otherwise):
++
+- The elements increase or stay the same until the “apex” (the highest value) is reached
+- After the apex, any remaining values decrease or remain the same
+
Solution: {∧/(⍳∘≢≡⍋)¨(⊂((⊢⍳⌈/)↑⊢),⍵),⊂⌽((⊢⍳⌈/)↓⊢),⍵}
How do we approach this? First we have to split the vector at the “apex”. The train (⊢⍳⌈/) will return the index of (⍳) the maximum element.
Combined with Take (↑) and Drop (↓), we build a two-element vector containing both parts, in ascending order (we Reverse (⌽) one of them). Note that we have to Ravel (,) the argument to avoid rank errors in Index Of.
Next, (⍳∘≢≡⍋) on each of the two vectors will test if they are non-decreasing (i.e. if the ranks of all the elements correspond to a simple range from 1 to the size of the vector).
10. Stacking It Up
+++Write a function that takes as its right argument a vector of simple arrays of rank 2 or less (scalar, vector, or matrix). Each simple array will consist of either non-negative integers or printable ASCII characters. The function must return a simple character array that displays identically to what
+{⎕←⍵}¨displays when applied to the right argument.
Solution: {↑⊃,/↓¨⍕¨⍵}
The first step is to Format (⍕) everything to get strings. A lot of trial-and-error is always necessary when dealing with nested arrays, and this being about formatting exacerbates the problem.
+
+ The next step would be to “stack everything vertically”, so we will need Mix (↑) at some point. However, if we do it immediately we don’t get the correct result:
Mix is padding with spaces both horizontally (necessary as we want the output to be a simple array of characters) and vertically (not what we want). We will have to decompose everything line by line, and then mix all the lines together. This is exactly what SplitSplit is the dual of Mix.
+
+ (↓) does:
{↓¨⍕¨⍵}(3 3⍴⍳9)(↑'Adam' 'Michael')(⍳10) '*'(5 5⍴⍳25)
+┌───────────────────┬─────────────────┬──────────────────────┬─┬───────────────
+│┌─────┬─────┬─────┐│┌───────┬───────┐│┌────────────────────┐│*│┌──────────────
+││1 2 3│4 5 6│7 8 9│││Adam │Michael│││1 2 3 4 5 6 7 8 9 10││ ││ 1 2 3 4 5
+│└─────┴─────┴─────┘│└───────┴───────┘│└────────────────────┘│ │└──────────────
+└───────────────────┴─────────────────┴──────────────────────┴─┴───────────────
+
+ ─────────────────────────────────────────────────────────────┐
+ ┬──────────────┬──────────────┬──────────────┬──────────────┐│
+ │ 6 7 8 9 10│11 12 13 14 15│16 17 18 19 20│21 22 23 24 25││
+ ┴──────────────┴──────────────┴──────────────┴──────────────┘│
+ ─────────────────────────────────────────────────────────────┘Next, we clean this up with Ravel (,) and we can Mix to obtain the final result.
Table of Contents
-The Ising model is a model used to represent magnetic dipole moments in statistical physics. Physical details are on the Wikipedia page, but what is interesting is that it follows a complex probability distribution on a lattice, where each site can take the value +1 or -1.
-
Mathematical definition
-We have a lattice \(\Lambda\) consisting of sites \(k\). For each site, there is a moment \(\sigma_k \in \{ -1, +1 \}\). \(\sigma = -(\sigma_k)_{k\in\Lambda}\) is called the configuration of the lattice.
-The total energy of the configuration is given by the Hamiltonian \[ -H(\sigma) = -\sum_{i\sim j} J_{ij}\, \sigma_i\, \sigma_j, -\] where \(i\sim j\) denotes neighbours, and \(J\) is the interaction matrix.
-The configuration probability is given by: \[ -\pi_\beta(\sigma) = \frac{e^{-\beta H(\sigma)}}{Z_\beta} -\] where \(\beta = (k_B T)^{-1}\) is the inverse temperature, and \(Z_\beta\) the normalisation constant.
-For our simulation, we will use a constant interaction term \(J > 0\). If \(\sigma_i = \sigma_j\), the probability will be proportional to \(\exp(\beta J)\), otherwise it would be \(\exp(\beta J)\). Thus, adjacent spins will try to align themselves.
-Simulation
-The Ising model is generally simulated using Markov Chain Monte Carlo (MCMC), with the Metropolis-Hastings algorithm.
-The algorithm starts from a random configuration and runs as follows:
--
-
- Select a site \(i\) at random and reverse its spin: \(\sigma'_i = -\sigma_i\) -
- Compute the variation in energy (hamiltonian) \(\Delta E = H(\sigma') - H(\sigma)\) -
- If the energy is lower, accept the new configuration -
- Otherwise, draw a uniform random number \(u \in ]0,1[\) and accept the new configuration if \(u < \min(1, e^{-\beta \Delta E})\). -
Implementation
-The simulation is in Clojure, using the Quil library (a Processing library for Clojure) to display the state of the system.
-This post is “literate Clojure”, and contains core.clj. The complete project can be found on GitHub.
The application works with Quil’s functional mode, with each function taking a state and returning an updated state at each time step.
-The setup function generates the initial state, with random initial spins. It also sets the frame rate. The matrix is a single vector in row-major mode. The state also holds relevant parameters for the simulation: \(\beta\), \(J\), and the iteration step.
(defn setup [size]
- "Setup the display parameters and the initial state"
- (q/frame-rate 300)
- (q/color-mode :hsb)
- (let [matrix (vec (repeatedly (* size size) #(- (* 2 (rand-int 2)) 1)))]
- {:grid-size size
- :matrix matrix
- :beta 10
- :intensity 10
- :iteration 0}))Given a site \(i\), we reverse its spin to generate a new configuration state.
-(defn toggle-state [state i]
- "Compute the new state when we toggle a cell's value"
- (let [matrix (:matrix state)]
- (assoc state :matrix (assoc matrix i (* -1 (matrix i))))))In order to decide whether to accept this new state, we compute the difference in energy introduced by reversing site \(i\): \[ \Delta E = -J\sigma_i \sum_{j\sim i} \sigma_j. \]
-The filter some? is required to eliminate sites outside of the boundaries of the lattice.
(defn get-neighbours [state idx]
- "Return the values of a cell's neighbours"
- [(get (:matrix state) (- idx (:grid-size state)))
- (get (:matrix state) (dec idx))
- (get (:matrix state) (inc idx))
- (get (:matrix state) (+ (:grid-size state) idx))])
-
-(defn delta-e [state i]
- "Compute the energy difference introduced by a particular cell"
- (* (:intensity state) ((:matrix state) i)
- (reduce + (filter some? (get-neighbours state i)))))We also add a function to compute directly the hamiltonian for the entire configuration state. We can use it later to log its values across iterations.
-(defn hamiltonian [state]
- "Compute the Hamiltonian of a configuration state"
- (- (reduce + (for [i (range (count (:matrix state)))
- j (filter some? (get-neighbours state i))]
- (* (:intensity state) ((:matrix state) i) j)))))Finally, we put everything together in the update-state function, which will decide whether to accept or reject the new configuration.
(defn update-state [state]
- "Accept or reject a new state based on energy
- difference (Metropolis-Hastings)"
- (let [i (rand-int (count (:matrix state)))
- new-state (toggle-state state i)
- alpha (q/exp (- (* (:beta state) (delta-e state i))))]
- ;;(println (hamiltonian new-state))
- (update (if (< (rand) alpha) new-state state)
- :iteration inc)))The last thing to do is to draw the new configuration:
-(defn draw-state [state]
- "Draw a configuration state as a grid"
- (q/background 255)
- (let [cell-size (quot (q/width) (:grid-size state))]
- (doseq [[i v] (map-indexed vector (:matrix state))]
- (let [x (* cell-size (rem i (:grid-size state)))
- y (* cell-size (quot i (:grid-size state)))]
- (q/no-stroke)
- (q/fill
- (if (= 1 v) 0 255))
- (q/rect x y cell-size cell-size))))
- ;;(when (zero? (mod (:iteration state) 50)) (q/save-frame "img/ising-######.jpg"))
- )And to reset the simulation when the user clicks anywhere on the screen:
-(defn mouse-clicked [state event]
- "When the mouse is clicked, reset the configuration to a random one"
- (setup 100))(q/defsketch ising-model
- :title "Ising model"
- :size [300 300]
- :setup #(setup 100)
- :update update-state
- :draw draw-state
- :mouse-clicked mouse-clicked
- :features [:keep-on-top :no-bind-output]
- :middleware [m/fun-mode])Conclusion
-The Ising model is a really easy (and common) example use of MCMC and Metropolis-Hastings. It allows to easily and intuitively understand how the algorithm works, and to make nice visualizations!
-Table of Contents
-
-
- What is an L-system? -
- Implementation details -
- Common file format for L-systems -
- Variations on L-systems -
- Usage notes -
- References -
L-systems are a formal way to make interesting visualisations. You can use them to model a wide variety of objects: space-filling curves, fractals, biological systems, tilings, etc.
-See the Github repo: https://github.com/dlozeve/lsystems
-What is an L-system?
-A few examples to get started
-
Definition
-An L-system is a set of rewriting rules generating sequences of symbols. Formally, an L-system is a triplet of:
--
-
- an alphabet \(V\) (an arbitrary set of symbols) -
- an axiom \(\omega\), which is a non-empty word of the alphabet (\(\omega \in V^+\)) -
- a set of rewriting rules (or productions) \(P\), each mapping a symbol to a word: \(P \subset V \times V^*\). Symbols not present in \(P\) are assumed to be mapped to themselves. -
During an iteration, the algorithm takes each symbol in the current word and replaces it by the value in its rewriting rule. Not that the output of the rewriting rule can be absolutely anything in \(V^*\), including the empty word! (So yes, you can generate symbols just to delete them afterwards.)
-At this point, an L-system is nothing more than a way to generate very long strings of characters. In order to get something useful out of this, we have to give them meaning.
-Drawing instructions and representation
-Our objective is to draw the output of the L-system in order to visually inspect the output. The most common way is to interpret the output as a sequence of instruction for a LOGO-like drawing turtle. For instance, a simple alphabet consisting only in the symbols \(F\), \(+\), and \(-\) could represent the instructions “move forward”, “turn right by 90°”, and “turn left by 90°” respectively.
-Thus, we add new components to our definition of L-systems:
--
-
- a set of instructions, \(I\). These are limited by the capabilities of our imagined turtle, so we can assume that they are the same for every L-system we will consider:
-
-
-
Forwardmakes the turtle draw a straight segment.
-TurnLeftandTurnRightmakes the turtle turn on itself by a given angle.
-PushandPopallow the turtle to store and retrieve its position on a stack. This will allow for branching in the turtle’s path.
-Stay, which orders the turtle to do nothing.
-
- - a distance \(d \in \mathbb{R_+}\), i.e. how long should each forward segment should be. -
- an angle \(\theta\) used for rotation. -
- a set of representation rules \(R \subset V \times I\). As before, they will match a symbol to an instruction. Symbols not matched by any rule will be associated to
Stay.
-
Finally, our complete L-system, representable by a turtle with capabilities \(I\), can be defined as \[ L = (V, \omega, P, d, \theta, -R). \]
-One could argue that the representation is not part of the L-system, and that the same L-system could be represented differently by changing the representation rules. However, in our setting, we won’t observe the L-system other than by displaying it, so we might as well consider that two systems differing only by their representation rules are different systems altogether.
-Implementation details
-The LSystem data type
-The mathematical definition above translate almost immediately in a Haskell data type:
--- | L-system data type
-data LSystem a = LSystem
- { name :: String
- , alphabet :: [a] -- ^ variables and constants used by the system
- , axiom :: [a] -- ^ initial state of the system
- , rules :: [(a, [a])] -- ^ production rules defining how each
- -- variable can be replaced by a sequence of
- -- variables and constants
- , angle :: Float -- ^ angle used for the representation
- , distance :: Float -- ^ distance of each segment in the representation
- , representation :: [(a, Instruction)] -- ^ representation rules
- -- defining how each variable
- -- and constant should be
- -- represented
- } deriving (Eq, Show, Generic)Here, a is the type of the literal in the alphabet. For all practical purposes, it will almost always be Char.
Instruction is just a sum type over all possible instructions listed above.
Iterating and representing
-From here, generating L-systems and iterating is straightforward. We iterate recursively by looking up each symbol in rules and replacing it by its expansion. We then transform the result to a list of Instruction.
Drawing
-The only remaining thing is to implement the virtual turtle which will actually execute the instructions. It goes through the list of instructions, building a sequence of points and maintaining an internal state (position, angle, stack). The stack is used when Push and Pop operations are met. In this case, the turtle builds a separate line starting from its current position.
The final output is a set of lines, each being a simple sequence of points. All relevant data types are provided by the Gloss library, along with the function that can display the resulting Picture.
Common file format for L-systems
-In order to define new L-systems quickly and easily, it is necessary to encode them in some form. We chose to represent them as JSON values.
-Here is an example for the Gosper curve:
-{
- "name": "gosper",
- "alphabet": "AB+-",
- "axiom": "A",
- "rules": [
- ["A", "A-B--B+A++AA+B-"],
- ["B", "+A-BB--B-A++A+B"]
- ],
- "angle": 60.0,
- "distance": 10.0,
- "representation": [
- ["A", "Forward"],
- ["B", "Forward"],
- ["+", "TurnRight"],
- ["-", "TurnLeft"]
- ]
-}Using this format, it is easy to define new L-systems (along with how they should be represented). This is translated nearly automatically to the LSystem data type using Aeson.
Variations on L-systems
-We can widen the possibilities of L-systems in various ways. L-systems are in effect deterministic context-free grammars.
-By allowing multiple rewriting rules for each symbol with probabilities, we can extend the model to probabilistic context-free grammars.
-We can also have replacement rules not for a single symbol, but for a subsequence of them, thus effectively taking into account their neighbours (context-sensitive grammars). This seems very close to 1D cellular automata.
-Finally, L-systems could also have a 3D representation (for instance space-filling curves in 3 dimensions).
-Usage notes
--
-
- Clone the repository:
git clone [[https://github.com/dlozeve/lsystems]]
- - Build:
stack build
- - Execute
stack exec lsystems-exe -- examples/penroseP3.jsonto see the list of options
- - (Optional) Run tests and build documentation:
stack test --haddock
-
Usage: stack exec lsystems-exe -- --help
lsystems -- Generate L-systems
-
-Usage: lsystems-exe FILENAME [-n|--iterations N] [-c|--color R,G,B]
- [-w|--white-background]
- Generate and draw an L-system
-
-Available options:
- FILENAME JSON file specifying an L-system
- -n,--iterations N Number of iterations (default: 5)
- -c,--color R,G,B Foreground color RGBA
- (0-255) (default: RGBA 1.0 1.0 1.0 1.0)
- -w,--white-background Use a white background
- -h,--help Show this help text
-Apart from the selection of the input JSON file, you can adjust the number of iterations and the colors.
-stack exec lsystems-exe -- examples/levyC.json -n 12 -c 0,255,255
References
--
-
- Prusinkiewicz, Przemyslaw; Lindenmayer, Aristid (1990). The Algorithmic Beauty of Plants. Springer-Verlag. ISBN 978-0-387-97297-8. http://algorithmicbotany.org/papers/#abop -
- Weisstein, Eric W. “Lindenmayer System.” From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/LindenmayerSystem.html -
- Corte, Leo. “L-systems and Penrose P3 in Inkscape.” The Brick in the Sky. https://thebrickinthesky.wordpress.com/2013/03/17/l-systems-and-penrose-p3-in-inkscape/ -