- Divided into smaller sub-problems (divide),
- The sub-problems are solved (conquer),
- And the solutions to sub-problems are combined to solve the original and “bigger” problem (combine).

Instead of constructing indefinite number of nested loops destroying the readability of the code and the performance of execution, the “recursive” way utilizes just one block of code which calls itself (hence the term “recursive”) for the smaller problem. The main point is to define a “stop” rule, so that the function does not sink into an infinite recursion depth. While nested loops modify the same object (or address space in the low level sense), recursion moves the “stack pointer”, so each recursion depth uses a different part of the stack (a copy of the objects will be created for each recursion). This illustrates a well-known trade-off in algorithm design: Memory versus performance; recursion enhances performance at the expense of using more memory.

Given the benefits of the recursive divide-and-conquer paradigm, however, the performance of the naive recursion may not still be as good as we may want, since there may be too many smaller problems to navigate through recursion. Sometimes, the same or a similar sub-problem may be visited multiple times and each sub-problem paves way to its sub-problems down to the stop rule – so a naive recursive approach solves the same sub-problems many and many times. The recipe to enhance performance is to “memoize” solved sub-problems: Simply to keep solutions in a look-up table, and first check the table if the value was memoized before commencing with the “divide-and-conquer”. Divide-and-conquer paradigm on steroids with memoization is known as “dynamic programming” or “dp” for short (Cormen, et al. 2009, p. 359).

Now we will illustrate an application of the “dynamic programming” approach using a question from the PuzzlOR site – a site which publishes bi-monthly decision support puzzles for applied mathematicians. Since the answers to expired questions are disclosed by the website itself, there is no problem in describing a solution in detail. The question that appeared in October 2014 is called “Fighters”. The question involves a quadruel: A kind of duel fight that is among four participants instead of the classical two participant case.

The question goes as follows:

Four different fighters are having an all-out battle to determine who among them is the strongest. The image above shows those four fighters: Allan, Barry, Charles, and Dan.Each fighter has varying attack and health abilities. At the start of the battle, they have differing health points: Allan has 10, Barry has 12, Charles has 16, and Dan has 18. Also, each fighter has differing attack points: Allan has 4, Barry has 3, Charles has 2, and Dan has 1.The battle takes place over multiple rounds, each round consisting of a single attack. In each round, one random attacker and one random defender are chosen. When the attacker attacks a defender, the defender loses health points in the amount equivalent to the attacker’s attack points. For example, if Allan is the attacker and Barry is the defender, Barry would lose 4 health points.The fighters continue to randomly attack and defend in subsequent rounds until there is only one fighter left, who is then declared the winner. A fighter is removed from the battle when his life points become zero (or less).Question: Which fighter is most likely to win the battle?

Starting with the initial health levels of 10, 12, 16, and 18, we can divide the main problem of finding the most likely winner into smaller ones, each of which is the ex-post health level configuration after a possible and equally likely attack between any “alive” fighters – “alive” meaning having a non-zero health level. We can divide the problem until we have only one alive fighter – which is the winner. We calculate the actual probability of winning here in this deepest recursion step at the stop rule (only one winner and the probability is zero for other fighters), so this is the “conquer” phase in which the smaller problem is solved.

Each recursion paves way into smaller problems – each one involving a different pair of fighters – number of which is the count of two-member permutations of alive fighters without replacement . We calculate permutations instead of combinations since Fighter X attacking Fighter Y is different than the other way around – the attacker injures the attacked one, so the attack is one-way and not symmetric. We calculate the permutations without replacement, since an attacker cannot attack herself/himself! So in the case of four alive fighters, there are 12 possible attacks. And each attack has a probability of 1/12. The winning probabilities of each fighter calculated from each of these 12 attacks are summed together in a bottom-up manner (from the deepest recursion level involving the last alive fighter, up to the initial question) and this step is the “combine” phase: In each step involving n number of ordered pairs of fighters, the 1/n probability is multiplied with the bottom-up probabilities to allow for carrying the joint probabilities up to the topmost level (the initial question).

We first start with the divide-and-conquer paradigm without memoization and see the performance:

## for non-repeated permutations library(gtools) ### non memoized version # get the winning probabilities of each player # win means only one player remains with non-zero health puzzlor.fighters.nm <- function() { health <- c(10, 12, 16, 18) # initial health points of players a:d attack <- c(4, 3, 2, 1) # attack points of players a:d # permutations of players. first column attackers, second column attacked. # self attack not allowed (no repetitions) at.combs <- gtools::permutations(4, 2, 1:4) # get matrix of results through ply results <- t(apply(at.combs, 1, fight.rec.nm, health, attack, 1/12)) # sum all probabilities of each player result <- colSums(results) # update player names names(result) <- c("Allan", "Barry", "Charles", "Dan") return(result) } # recursive function to recalculate the win probabilities after an attack, # incorporating bottom-up cumulative joint probabilities fight.rec.nm <- function(at, health, attack, prob) { # at: vector of attacker and attacked no's # at1: attacker no # at2: attacked no # health: vector for current health levels of all players # attack: vector for attack points of all players # prob: carried joint probability at1 <- at[1] # attacker no at2 <- at[2] # attacked no # update the health of at2 after attack health[at2] <- max(health[at2] - attack[at1], 0) players <- which(health > 0) # index of alive players # if only one player survives, so is the winner if (length(players) == 1) { result <- rep(0,4) # vector for results result[players] <- prob # update the probability of the winner return(result) } else { # permutations of indices of attackers and attackeds. # self attack not allowed (no repetitions in permutations) at.combs <- gtools::permutations(length(players), 2, players) probs <- nrow(at.combs) # total count of attack permutations prob.new <- 1 / probs # probability of each attack permutation # collect matrix of results recursing the function results <- t(apply(at.combs, 1, fight.rec.nm, health, attack, prob.new)) result <- colSums(results) # sum the probabilities for each player result <- prob * result # multiply with the joint probability return(result) } }

The key points here in the above code are:

- In each recursion, the sub-problem is defined by the attacker, attacked, current health points, and the carried probability of that sub-problem (1 / the number of sub-problems generated one level up). The attack points of fighters do not change (line 35).
- In each sub-problem the health points are updated by subtracting the attack point of the attacker from the health point of the attacked. Note that, the health level can be as low as 0 (line 48).
- New sub-problems are defined for the non-repeated ordered-pairs of remaining alive fighters – provided that more than one fighter is alive. The permutations are generated by the use of “gtools” package, a handy package to generate combinations and permutations with different options. Note that the base function “expand.grid” generates only the repeated permutations which is not useful for our purpose, since a fighter is not supposed to attack herself/himself. The carried probability is 1 / the number of new sub-problems. The same updated health levels are carried to the new recursion level in each new sub-problem (lines 62-68). If only one fighter with non-zero health level remains, the recursion stops and the carried probability is recorded for the alive fighter (0 probability for the others) (lines 53-56).
- Each sub-problem returns a 4-tuple vector of winning probabilities for each fighter. The probabilities account for the carried probability (1 / number of ordered pairs). All results from sub-problems at a certain level return a matrix and the columns are summed up and multiplied by the carried probability to yield a vector of combined probabilities to be carried one level up until the original problem (lines 69-71).

We construct two functions: The first one (“puzzlor.fighters.nm”) (line 9) is a wrapper to initiate necessary objects and to start the initial problem calling the recursive function. The second function (“fight.rec.nm”) (line 35) is the recursive function which does all the steps described in detail above.

The performance of this non-memoized version is … well it is still running as at the time I am writing these sentences.

The reason for the performance bottleneck is the fact that too many sub-problems are generated until only one winner remains and the number of “tree branches” grows geometrically with each recursion (multiplied by the number of possible ordered attack pairs in each step). And those sub-problems are revisited many times since the fights can yield the same health level configuration at different instances . We can enhance the performance without redefining the problem but by augmenting it with an additional data structure:

- Each sub-problem, after the attack is realized, is defined by the 4-tuple of health levels (line 48).
- And each of these sub-problems yields a 4-tuple of the winning probabilities of each fighter (line 56, line 71).
- So for each sub-problem defined by the 4-tuple configuration of health levels, we can memoize the 4-tuple winning probabilities of each fighter.
- We must create an object to record the winning probabilities for each health level configuration. But what should be the R object type for this task? We can subset a matrix with at most two indices: row and column, but we need four indices for each health level. And in order to subset a vector (the 4-tuple probabilities), we can at most use one index in a matrix. The solution is to use a higher dimension object called “array” in R. Arrays can have a dimension higher than two. In our example we need to have a 5D array to index for each health configuration and still yield a 1D vector object of probabilities

The new code which combines recursion and memoization and hence is an application of “dynamic programming” (dp) is as follows:

## for non-repeated permutations library(gtools) # get the winning probabilities of each player # win means only one player remains with non-zero health puzzlor.fighters.5D <- function() { health <- c(10, 12, 16, 18) # initial health points of players a:d attack <- c(4, 3, 2, 1) # attack points of players a:d ## global 5D object for memoization result.ar <<- array(dim = c(health + 1, 4), dimnames = NULL) # permutations of players. first column attackers, second column attacked. # self attack not allowed (no repetitions). at.combs <- gtools::permutations(4, 2, 1:4) # get matrix of results through ply results <- t(apply(at.combs, 1, fight.rec.5D, health, attack, 1/12)) result <- colSums(results) # sum all probabilities of each player names(result) <- c("Allan", "Barry", "Charles", "Dan") # update player names return(result) } # recursive function to recalculate the win probabilities after an attack, # incorporating bottom-up cumulative joint probabilities fight.rec.5D <- function(at, health, attack, prob) { # at: vector of attacker and attacked no's # at1: attacker no # at2: attacked no # health: vector for current health levels of all players # attack: vector for attack points of all players # prob: carried joint probability at1 <- at[1] # attacker no at2 <- at[2] # attacked no # update the health of at2 after attack health[at2] <- max(health[at2] - attack[at1], 0) # create a matrix for indexing the 5D array mat.index <- cbind(matrix(rep(health + 1, 4), nrow = 4, byrow = T), 1:4) mat.row <- "["(result.ar, mat.index) # get the memoized health configuration # if the probabilities for the health conf. are already memoized if (!is.na(mat.row[1])) { # get the result and multiply with the prob result <- prob * mat.row return(result) } else { players <- which(health > 0) # index of alive players # if only one player survives, so is the winner if (length(players) == 1) { # result <- rep(0,4) # vector for results result[players] <- prob # update the probability of the winner return(result) } else { # permutations of indices of attackers and attackeds. # self attack not allowed (no repetitions in permutations) at.combs <- gtools::permutations(length(players), 2, players) probs <- nrow(at.combs) # total count of attack permutations prob.new <- 1 / probs # probability of each attack permutation # collect matrix of results recursing the function results <- t(apply(at.combs, 1, fight.rec.5D, health, attack, prob.new)) result <- colSums(results) # sum the probabilities for each player # memoize the health configuration "["(result.ar, mat.index) <<- result # multiply with the joint probability result <- prob * result return(result) } } }

The main points of the modified code are as follows:

- We create a 5D array, to hold the memoized values. The first four dimensions are for the health values and the last dimension is for the memoized probabilities. Note that, each dimension is 1 + the initial health in order to account for 0 health level. When subsetting the array, the health values will be incremented (line 13).
- And also note that, the array is superassigned (“<<-“) so it is defined in the global environment. This is important, since whatever the recursion level is, the memoization should be done into one global object, not a copy of that object at the current depth – otherwise we cannot collect all results and look them up from the same collection.
- We can subset a 5D object with 4 indices (leaving the last index so that a vector is returned) just like we subset a matrix for rows and/or columns: array.object[x,y,z,a,]. However, we should first subset each index from the health level configuration vector – which does not look very pretty considering higher dimension objects. A more elegant and generalized way to subset the 5D array is to use the subsetting operator “[” as a function: Note that every operator in R is a function following the convention of LISP (one of the languages that inspired S and R) (line 48). We should first create a matrix containing the index values and then use the “[” operator as a function (line 46). Note that we can follow the same subsetting approach for modifying the array for memoization purposes (line 81).
- In each sub-problem, the array will be indexed first to get the memoized probabilities – if any. If the probabilities for the given health level configuration is already memoized, those values will be returned (after multiplying with the probability carried from one higher level) (lines 51-55).
- For the last level of recursion (where there is only one fighter left), there is no need to memoize the value, since just returning a vector of three 0’s and one value for carried probability without further recursion should be as efficient as indexing the 5D array (lines 62-66).
- For the cases in which the configuration is not yet memoized and there are more than one alive fighters, the recursion will be followed as in the first code (lines 72-85), however, this time the results will be memoized into the 5D array before multiplying with the carried probability (line 81).

Now let’s execute this new function and compare its results with the first one. Please note that the first code is still running after several hours on an octacore machine at 2.7 GHz current and 3.5 GHz max clock speed!

> microbenchmark(puzzlor.fighters.5D(), times = 1) Unit: seconds expr min lq mean median uq max puzzlor.fighters.5D() 7.983507 7.983507 7.983507 7.983507 7.983507 7.983507 neval 1

The code executes in less than 8 seconds. You can check the correct solutions from the PuzzlOR web page.

A downside to the dp approach is that, it cannot be utilized with the “mcmapply” and similar R functions of the “parallel” package for multicore parallelism, since these functions create copies of the global objects as well as the local objects for each thread. So all values cannot be memoized into the same object. However the gains from the memoization outweigh the disadvantage of using a single core.

You can view the whole code (with an additional option with a 4D array and a matrix, which slightly underperforms the 5D case) from my Github page.

Notes:

- In cases where R returns an error for reaching the maximum recursion depth, the limit for recursion depth up to 5e5 levels can be defined with:
> options(expressions= 500000)

- Sometimes, despite an enough level of maximum recursion depth, the usage of stack may reach the maximum available size. To tackle with this issue answers to a stackoverflow question summarize the solution. From the unix shell, to check the current stack size:
$ ulimit -s # print default 8192

Now you can increase the stack size by:

$ ulimit -s 16384 # enlarge stack limit to 16 megs

And now, you can check the new stack size available to R by calling a new R session from the shell that you increased the stack size with and calling the following:

> Cstack_info()["size"] size 15938355

The methods to enhance a brute force algorithm may be to restrict the set of values to be searched, to rewrite the code in order to choose more efficient ways that run with fewer instructions at the hardware level or to incorporate multicore or GPU parallelism.

In this post, I will illustrate three “brute force” methods to solve a problem, enhancing in each successive method so that we experience a performance gain of more than 100x times in the third method as compared to the first one. A knowledge of the low level mechanics of the R interpreter is important in selecting the optimal method. The problem I choose is the 124th problem in Project Euler. Since Project Euler requests participants not to disclose the solutions publicly, I will just share portions of the algorithm in order to emphasize the efficiency gains from choosing the right tools and methods:

## Ordered radicals

## Problem 124

The radical of

n, rad(n), is the product of the distinct prime factors ofn. For example, 504 = 2^{3}× 3^{2}× 7, so rad(504) = 2 × 3 × 7 = 42.If we calculate rad(

n) for1≤n≤ 10, then sort them on rad(n), and sorting onnif the radical values are equal, we get:

UnsortedSorted

n

rad(n)

n

rad(n)

k 1 1 1 1 1 2 2 2 2 2 3 3 4 2 3 4 2 8 2 4 5 5 3 3 5 6 6 9 3 6 7 7 5 5 7 8 2 6 6 8 9 3 7 7 9 10 10 10 10 10Let E(

k) be thekth element in the sortedncolumn; for example, E(4) = 8 and E(6) = 9.If rad(

n) is sorted for 1 ≤n≤ 100000, find E(10000).

So the main point is to calculate the product of the distinct or unique prime factors of each number up to 1e5. Prime sieving is the trivial part, since it is a common task in many Project Euler problems. Although I have my own sieve algorithms, I personally prefer to use the “Primes” function in “numbers” package. Another method is making a system call to open-source “primesieve” tool. The last part (ordering and getting the number at certain rank) is also trivial.

The common method in all three versions is to check the divisibility of all numbers <= 1e5 with sieved primes. In each version, we track a vector of boolean values with a length of 1e5 to record divisibility and a vector of 1’s with a length of 1e5 to track multiplied distinct factors (in order to yield radicals).

In the naivest first version, we check the divisibility of each number <= 1e5 with a certain prime by the modulo operator “%%”. Note that the numbers vector will be probed against all primes below 1e5:

# whether numbers are divisible by the prime factor, by modulo bool.new <- nums %% prime == 0 # update the radicals vector with the prime factor radicals[bool.new] <<- radicals[bool.new] * prime

In this code snippet, “nums” is the global vector holding the sequence of numbers 1 to 1e5. Although this option is the “brutest” of the three, it still enjoys the efficiency benefit of indexing a vector with another boolean vector of the same length – which takes linear time – instead of searching for the numbers. As long as the vector is not too large and the numbers to be indexed are not too sparse, we may not need to have an efficient hash function without collision.

The performance of the first version that runs the above snippet for radical updating is 12.8 seconds:

> microbenchmark(ordered.radicals.3(), times = 10) Unit: seconds expr min lq mean median uq max ordered.radicals.3() 12.61865 12.64935 12.80029 12.75633 12.85513 13.39739 neval 10

However, we do not need to probe the divisibility of the numbers in the sequence with the prime using the modulo operator “%%”. We know that the numbers divisible by a prime are the multiples of that prime! Using the vector indexing approach, we can get a sequence of the multiples of the prime upto 1e5, and toggle those indices in a clean boolean vector of FALSE values to TRUE, yielding the same result in a much faster manner:

# empty F vector bool.new <- bool.new1 # the indices that are divisible by the prime, as a sequence, without modulo divisible.inds <- (1:floor(lim / prime)) * prime # whether numbers are divisible by the prime factor bool.new[divisible.inds] <- T # update the radicals vector with the prime factor radicals[bool.new] <<- radicals[bool.new] * prime

In the above snippet, first of all an empty boolean vector of 1e5 FALSE values (bool.new) is copied from a template vector (bool.new1) in line 2. In the most critical part, a sequence of the multiples of the prime below 1e5 is formed (line 5). And then the boolean vector of FALSE values is indexed with the multiples and those values are toggled to TRUE (line 8). The last line for updating the radicals is the same line in the first option (line 11).

The performance of this second option is 3 seconds, recording more than 4x gain on top of the first option:

> microbenchmark(ordered.radicals.2(), times = 10) Unit: seconds expr min lq mean median uq max ordered.radicals.2() 2.933848 2.982204 3.022639 3.016654 3.066471 3.118391 neval 10

However, this performance is not good enough. There are 9592 primes below 1e5. And both codes check for all those primes. To enhance the performance we should first review sieving algorithms and reflect on the problem at hand:

A number which is not prime is a composite number: It has at least two prime factors. A composite number must have at least one prime factor less than or equal to its square root. We can arrive at this conclusion by induction: If the smallest prime factor is larger than the square root, the other factor(s) must be larger further. However, this is a contradiction since the product of those factors would be larger than the number itself then. In prime sieving up to some limit, the prime factors <= the square root of the limit are probed. Following this logic, we can deduce that a composite number can have at most one distinct prime factor larger than its square root. Using the same reasoning, we can probe the primes <= sqrt(1e5). If we divide all numbers <= 1e5 with all powers of the primes <= sqrt(1e5) less than 1e5 and are divisible by the numbers, we will have the last remaining prime factors >= sqrt(1e5) for all numbers (and for the case of prime numbers, we will have the primes themselves). Square root of 1e5 is 316 and there are only 65 primes less than 316. So instead of 9592 primes, we will probe for only 65 primes. The maximum power of a prime that is below 1e5 is 16 for the case of prime number 2. So the gains from probing 147 times fewer primes justify the extra work for probing powers of those prime factors.

# empty F vector bool.new <- bool.new1 # empty vector of 1's radicals.new <<- radicals.new1 # maximum power of the prime to be checked max.power <- floor(log(lim, prime)) # the numbers whic are divisible by the prime divisible.inds <- (1:floor(lim / prime)) * prime # whether numbers are divisible by the prime factor bool.new[divisible.inds] <- T # update the radicals vector with the prime factor radicals[bool.new] <<- radicals[bool.new] * prime # find the maximum power of the prime that all numbers are divisible with sapply(1:max.power, update.radicals, prime, lim) # divide all numbers with the maximum power of the prime if it is divisible nums <<- nums / radicals.new

In this last version, we have an additional data structure: an interim vector of 1’s (radicals.new) to collect the max powers of a prime that numbers are divisible with (line 5). We get the maximum power of the the prime that is below 1e5 using logarithm (line 8). After updating the radicals as before, we repeat the procedure for each power of the prime to populate the interim vector of prime powers (radicals.new) using a ply formula (line 20). After the vector is populated, numbers vector is divided by this prime powers vector (line 23).

After all 65 primes are probed and numbers vector is updated for all primes, the remaining numbers vector is the last prime factors to be multiplied with the vector for collecting distinct prime factors (radicals) to get the resulting radicals.

The performance is 109 milliseconds for this last version:

> microbenchmark(ordered.radicals.1(), times = 10) Unit: milliseconds expr min lq mean median uq max ordered.radicals.1() 106.1128 107.258 110.9804 109.0543 110.4271 129.5838 neval 10

We sped up the code nearly 30 times as compared with the second version and nearly 120 times as compared with the first “raw” version. Though we are still in the realm of “brute forcing”, by incorporating a little number theory and a better way to check divisibility, our force is now much more elegant and intelligent and manipulates the low level mechanics of the higher level R language better.

Note: The code snippets given are just critical excerpts from the original code given in order to illustrate the progression in algorithm enhancement.

]]>There are several ways to edit and run Julia code. The Jupyter notebook “IJulia” is the easiest one to deploy:

```
julia> Pkg.add("IJulia")
julia> using IJulia
julia> notebook(detached=true)
```

These commands will install the IJulia package and start a Julia kernel inside the default browser, and will persist even if Julia shell is exited.

The notebook can even be invoked from the UNIX shell with the following one-liner command:

$ julia -e "using IJulia; notebook(detached=true)" &

The notebook is good for prototyping short code snippets and executing them as block “cells”. However for larger and more complex coding tasks, an integrated development environment (IDE) is preferred. “Getting Started with Julia”, a comprehensive beginner’s resource for the language by Ivo Balbaert, mentions the IDE options available for Julia as Julia Studio, Sublime-IJulia and Juno.

GUI based IDE’s usually absorb a large portion of the valuable CPU cycles and RAM that would otherwise be available for resource hungry data analytics tasks. For those coders who devote scarce computer resources to data tasks, lightweight solutions are preferred. Vim is among the most popular lightweight text editors and is highly versatile in a variety of languages due to its large base of plugins/extensions. For a coder who gets used to the practical key bindings, Vim acts as a natural extension of the body.

Although there are some extensions to be installed for Julia, for those who are used to work with Vim editor, there is not yet an integrated extension such as Nvim-R that combines Vim and R in a seamless manner. However several open-source extensions and applications can be combined to create a working environment for Julia through which you can edit Julia code with syntax highlighting and execute selected portions of the code inside the interactive shell.

As also mentioned in “Getting Started with Julia”, julia-vim plugin is a good starting point, but we have to combine it with some other resources to create an efficient environment. The prerequisites for the recipe are:

- Vim editor
- GNU Screen or Tmux terminal multiplexers
- julia-vim plugin for Vim
- vim-slime plugin for Vim
- and of course Julia source code or binaries

Most of the prerequsites are available in repositories of mainstream Linux distributions or can be downloaded and installed from the source code repositories following the links. Apart from the lightweight nature of Vim and the power and speed it derives from key bindings, another advantage of such a recipe is that, it can even be run on command-line (shell only) run levels and headless systems accessed through SSH.

We need the terminal multiplexer so that, both the code editor and interactive shell can be viewed at the same time , julia-vim plugin for Julia syntax highlighting in Vim and vim-slime plugin to send the Julia code from the Vim editor to interactive shell to be interpreted/compiled (not only for Julia, but any language with a REPL shell).

The default terminal multiplexer for vim-slime is GNU Screen and the below scripts will utilize it. However the option can be changed to Tmux easily by:

```
let g:slime_target = "tmux"
```

command from the Vim editor or appending to .vimrc file.

After the installations, we first create the steps that GNU screen will follow when invoked from a shell script:

layout new split -s screen 0 exec vim focus next screen 1 exec julia focus next layout save default

We will save these steps into a file (named “screenjulia” in my example) that will be invoked by the -c argument of screen command in order to override the default configuration and create a new layout for the session. The file creates a new layout, splits the screen session horizontally, names the upper region as “0” and executes vim inside (the julia filename will be appended by the script), names the lower region as “1” and executes Julia inside.

Due to difficulties to pass arguments to the custom configuration file, the optional filename argument will be inserted by the script, if provided:

#!/bin/bash # vim-julia-slime binder sed -i "s/^exec vim.*$/exec vim $1/" /path/to/screenjulia screen -S 1 -c /path/to/screenjulia

The sed command inserts the Julia filename if provided by the optional $1 argument to the end of “exec vim” line in the configuration file “screenjulia” (replace “/path/to/screenjulia” with your actual path). The default session name for screen is “1”.

In order to enable the Julia syntax highlighting feature, the filename must have a .jl extension. To enable these features, it is useful to add the following lines into the .vimrc file:

filetype plugin indent on syntax on

Also enabling line numbers will simplify the line selection in visual mode, so it should also be added to the .vimrc file:

set nu

To use this script as a custom command, it should be made executable and a symlink pointing to the file should be created inside a location that is declared in the $PATH environment variable of your .bashrc or /etc/environment (such as /usr/bin) (replace /path/to/script with your actual path to the script file):

$ chmod +x /path/to/juliascriptname

$ sudo ln -s /path/to/juliascriptname /usr/bin/juliascriptname

Now, invoking the scriptname (replace it with any name you’d like for your command) command from anywhere will initiate the following layout:

$ juliascriptname juliafilename.jl

The code lines are selected in visual mode in Vim and sent by C-c, C-c command (insert mode should be exited) – simply holding CTRL and double tapping c. In each session, the first time this sequence is entered, Vim will prompt for screen session and window name. The default values as defined in our script are 1 and 1. Anytime you can change these values with C-c, v or :SlimeConfig.

The visual mode works by v for character-wise and V for line-wise selections. Without any selection, C-c, C-c binding will send the current block. vip command also selects the current block or paragraph. To be more precise in the line range of code, you should enter VXgg where X is the number of last line to be selected (set nu option in .vimrc file is important in this sense)

With a small code sample for demonstration purposes, selecting first 3 lines with V3gg from Vim and sending to Julia window by C-c, C-c result in:

To switch between panes you should use C-a [tab] and the complete layout can be exited through screen by C-a \ binding.

With an Nvim-R like plugin, vim + Julia can be integrated better with autocompletion and an object browser.

]]>