text:
## Friday, January 20, 2017

### An Efficient Algorithm to Determine the Order of a Permutation.

Consider at set $S = \{1,2,3\ldots n\}$ and recall that $S_n$ a set of all permutations of $S$ form a group under composition (<$S_n, *$>). E.g. $S=\{1,2,3,4\}$ , $\pi_1 = [3\,4\,2\,1]$ and $\pi_2 = [4\,1\,3\,2]$. To understand composition operation it may be easier to look at a permutation $\pi \in S_n$ as bijective function $\pi : S \rightarrow S$. The composition $\pi_1 * \pi_2$ would be equivalent to $\pi_1(\pi_2(i)), \forall i\in S$ which results in $\pi_1 * \pi_2 = [1\,3\,2\,4]$. We use the notation $\pi^k = \pi*\pi*\pi\ldots *\pi$ to denote a composition applied $k$ times. Since <$S_n,*$> is a group it must have an identity element $e$ which is a permutation such that $e(i) = i, \forall i\in S$.

Now lets focus on the following question:
Given any permutation $\pi \in S_n$ is there an integer $k \geq 0$ such that $\pi^k = e$ ?

Lets try few examples with $\pi_1 = [3\, 4\, 2\, 1]$, $\pi_1^2 = [2\, 1\, 4\, 3]$, $\pi_1^3 = [4\, 3\, 1\, 2]$, $\pi_1^4 = [1\, 2\, 3\, 4]$. So for this example $\pi_1$ we were able to find $k=4$ such that $\pi_1^4 = e$. Later we will show that finding such a $k$ is always possible for any permutation $\pi$ and such a $k$ is called the order of the permutation.

Recall that every permutation can be decomposed into one or more disjoint cycles. We call a permutation cyclic if it has exactly one cycle in it. E.g. $\pi_c = [3\, 1\, 4\, 2]$ is a cyclic permutation in $S_4$ with cycle being $1 \rightarrow 3 \rightarrow 3 \rightarrow 4 \rightarrow 2 \rightarrow 1$.

Lemma: If $\sigma \in S_n$ is a cyclic permutation then $\sigma^n = e$

Proof: Since $\sigma$ is cyclic all elements $S$ form a directed cycle with $n$ nodes. Starting at any node $i$ we need to traverse $n$ edges to get back to $i$. This can be algebraically expressed as $\underbrace{\sigma(\sigma(\ldots \sigma(i))\ldots)}_{n\text{-times}} =i$
$\Box$.

The above lemma states that if a permutation is cyclic on $n$ elements then composing it $n$ times on itself results in an identity permutation. Now we are ready to prove our main theorem:

Theorem: $\forall \pi \in S_n$ there exists a $k \geq 0$ such that $\pi^k = e$.

Proof: Let $\pi = \sigma_1\sigma_2\ldots \sigma_m$ be the cycle decomposition of $\pi$. Let $|\sigma_i|$ denote size of cycle. Since each of $\sigma_i$ is cyclic permutation in $S_{|\sigma_i|}$ by the previous lemma $\sigma_i^{|\sigma_i|} = e_{|\sigma_i|}$. Where $e_{|\sigma_i|}$ is an identity permutation in $S_{|\sigma_i|}$. Also notice that any multiple of $|\sigma_i|$ is also going to result in an identity permutation (i.e. for any integer $q \geq 1$ , $\sigma_i^{q|\sigma_i|} = e_{|\sigma_i|}$). Let $t = |\sigma_1|\times |\sigma_2|\times \ldots |\sigma_m|$ then $\pi^t = \sigma_1^t\sigma_2^t\ldots \sigma_m^t$ since $t$ is a multiple of $|\sigma_i|, 1\leq i \leq m$ we have $\pi^t = e_{|\sigma_1|}e_{|\sigma_2|}\ldots e_{|\sigma_m|} = e$
$\Box$.

From the above proof its clear that choosing $k = |\sigma_1|\times |\sigma_2|\ldots |\sigma_m|$ which is a multiple of all cycle lengths guarantees that $\pi^k = e$. However from a computational efficiency purpose we are much better if we let that multiple be the least common multiple of $|\sigma_i|,\, 1\leq i\leq m$. Following is an efficient algorithm -- linear time and takes at most $n$ l.c.m computations -- to find the order of any permutation:

// Finds the order of a valid permutation ( 1 \leq perm[i] \leq n).
size_t FindPermutationOrder(const std::vector< unsigned int > &perm) {
size_t n = perm.size();
size_t perm_order = 1; // This overflows if the l.c.m does not fit in 64-bits.
std::vector< bool > cycle_mask(n, true);

// Find the disjoint cycle lengths in the permutation.
size_t cycle_index = 0;
while (cycle_index < n) {
// Find the start of the next disjoint cycle.
if (cycle_mask[cycle_index]) {
size_t i = cycle_index;
size_t cycle_length = 0;
do {
cycle_mask[i] = false;
i = perm[i] - 1;
cycle_length++;
} while (i != cycle_index);
perm_order = boost::math::lcm(perm_order, cycle_length);
}
cycle_index++;
}
return perm_order;
}

Following are few examples on the algorithm:
Order of [  5  3  1  2  4  ] is 5
Order of [  3  4  2  1  ] is 4
Order of [  1  2  3  4  5  ] is 1
Order of [  7  8  9  10  1  2  4  3  6  5  ] is 5

## Thursday, December 22, 2016

### Simpler Expected Analysis of Randomized Selection Algorithm.

In this post I present a simpler proof to derive expected asymptotic bounds for the Randomized Selection algorithm -- Las Vegas version to be more precise. Given a set $S = \{ a_1\,a_2\,a_3\ldots a_n \}$ of keys and an integer $k \in [1,n]$ the selection problem asks to report the $k^{th}$ largest key in the set $S$. Recall that this problem can be solved optimally using $O(n)$ comparison operations using the BFPRT (Blum, Floyd, Pratt, Rivest and Tarjan 1973) algorithm. However a correct implementation of BFPRT algorithm may take a lot of code. On the other hand randomized algorithms for many deterministic algorithms are often simpler to implement. The algorithm chooses a random pivot for partitioning the input and recursively picks on of the partition until the partition (sub problem) size is just one. Following is an implementation of RandSelect in around 25 lines.

Our goal here is the analyze the expected run-time and expected number of steps for the problem size converge to unity -- more precisely the expected value of the variable rand_walk_steps in the following code. I have purposefully kept the analysis very rigorous to make it accessible to everyone.

TL;DR "Expected number of comparisons in RandSelect is $\Theta(n)$ and expected number of rand_walk_steps is $\Theta(\log(n))$".

// Randomized selection algorithm to pick k^th largest element
// in a mutable array.
template < typename T >
T RandSelect(T *array, size_t n, size_t k) {
assert(k > 0 && k <= n);
size_t problem_size = n, select = k;

// What is the expected value of rand_walk_steps ??
size_t rand_walk_steps = 0;

while (problem_size > 1) {
size_t i = 0, j = problem_size - 1;
size_t ridx = RandArrayIndex(problem_size);
const T pivot = array[ridx];

// Partition the sub problem with pivot.
while (i <= j) {
if (array[i] <= pivot) {
i++;
} else if (array[j] > pivot) {
j--;
} else {
std::swap(array[i], array[j]);
}
}

problem_size = (i >= select) ? i : problem_size - i;
array = (i >= select) ? array : &array[i];
select = (i >= select) ? select : select - i;
++rand_walk_steps;
}
return array[0];
}

Let $X_i$ be the random variable corresponding to problem size in $i^{th}$ step ($1 \geq i \geq n-1$) in the algorithm -- each step corresponds to one iteration of the outer most while loop. If the algorithm terminates in $j$ steps then all $X_k$ with $k>j$ are set to $0$.
$$\begin{array}{lll} X_i &=& \text{Random variable for problem size in }i^{th}\text{ iteration of the outer most}\textsf{ while}\text{ loop} \\ X_i && \left\{ \begin{array}{lr} \in [1,X_{i-1}-1] & \textsf{If } X_{i-1} > 1 \\ = 0 & \textsf{Otherwise} \end{array} \right. \\ X_0 &=& n \,\,\,\,\,\,\, \text{Initial problem size and }\,\,\, E[X_0] = n \\ && \\ \end{array}$$

Lemma-1: If $X_{i-1} >1$ Then $E[X_i | X_{i-1}] = \frac{X_{i-1}}{2}$

Proof: For a given $X_{i-1} > 1$ the random variable $X_{i}$ is uniformly distributed in the range $[1,X_{i-1}-1]$. Hence the expected value is simply the average $\frac{1 + (X_{i-1}-1)}{2} = \frac{X_{i-1}}{2}$ $\Box$.

Lemma-2: $E[X_i] = \frac{n}{2^i}$

Proof:
$$\begin{array}{llll} &E[X_i|X_{i-1}] &=& \frac{X_{i-1}}{2}\,\,\,\,\,\text{From Lemma-1}\\ &E[\,E[X_i|X_{i-1}]\,] &=& E[\frac{X_{i-1}}{2}] \\ \Rightarrow & E[X_i] &=& E[\frac{X_{i-1}}{2}] = \frac{E[X_{i-1}]}{2}\,\,\,\,\,\text{for any two random variables } X,Y\,\,E[E[X|Y]] = E[X]\\ &&=&\frac{1}{2}\times E[X_{i-1}] = \frac{1}{2^2}\times E[X_{i-2}] \ldots \,\,\,\,\, \text{Apply recursively} \\ &&\ldots& \\ &&=& \frac{1}{2^i}\times E[X_{i-i}] = \frac{E[X_0]}{2^i} = \frac{n}{2^i} \,\, \Box\\ \end{array}$$

Theorem-1: Expected number of comparisons done by RandSelect is $\Theta(n)$.

Proof: The number comparisons done by the algorithm in each iteration is $\leq 3\times X_i$. We at most $2$ comparisons in the inner if..else branch and $1$ comparison to test $i <= j$, so a total of $3$ comparisons in the each iteration of the inner most while loop. The loop runs at most $X_i$ times in the $i^{th}$ iteration of the outer most while loop. Finally we need $1$ comparison to test the termination of the while loop. Let $T$ be the random variable corresponding the total number of comparisons in the algorithm. So we need to find $E[T]$.
$$\begin{array}{ll} &T \leq 3(X_0 + X_1 + X_2 + \ldots X_{n-1}) + 1 & \text{Total number of comparisons in the Algorithm. }\\ &&\\ \Rightarrow& \sum_{i=0}^{n-1} X_i \leq T \leq 3\sum_{i=0}^{n-1} X_{i} + 1 & \\ & \sum_{i=0}^{n-1} E[X_i] \leq E[T] \leq 3\sum_{i=0}^{n-1} E[X_{i}] + 1 & \text{Linearity of expectation.} \\ &&\\ & \sum_{i=0}^{n-1} \frac{n}{2^i} \leq E[T] \leq 3\sum_{i=0}^{n-1} \frac{n}{2^i} + n & \text{Apply Lemma-2.} \\ & 2n\left(1-\frac{1}{2^n}\right) \leq E[T] \leq 6n\left(1-\frac{1}{2^n}\right) +n & \text{Use: } \sum_{i=0}^{n-1}\frac{1}{2^i} = \frac{1-1/2^n}{1-1/2}.\\ &&\\ \Rightarrow &n \leq E[T] < 6n + 1 = 6n +1 & \text{Use: for } n\geq 1\,\,\,,\frac{1}{2} \leq \left(1-\frac{1}{2^n}\right) < 1 \\ &&\\ \Rightarrow & E[T] = \Theta(n)\,\,\Box & \\ \end{array}$$

We now return the question of the expected number of steps for the problem size to reach unity (i.e. expected number of iterations of the outer most while loop). This is also referred as the expected number of recursive calls in Motwani and Raghavan's book. However the following analysis is much simpler than the application of a calculus based lemma based on random walks in the book. Let $Y_{i}$ be a random variable corresponding to the number of iterations starting with a problem size of $X_{i} > 1$.
$$\begin{array}{lll} Y_{i} &:=&\text{random variable for the number of steps with problem size }\,\, X_{i} >1 \\ Y_i &=& 1 + Y_{i+1} \,\,\,\,\,\, \text{(Recursive formulation)}\\ Y = Y_0 &=& \underbrace{1 + 1 + 1 \dots + 1}_{k \text{-times}} + Y_{k} \,\,\,\,\,\, \text{(Expand } k \text{ times)}\\ && \\ \end{array}$$

Theorem-2: Expected number of steps for the problem size to reach unity, $E[Y] = \Theta(\log(n))$.

Proof:
$$\begin{array}{ll} &Y_i \leq X_{i}\,\,\,\, \text{ (Number of steps are at most problem size at } i^{th} \text{ step which is } X_i)\\ & \\ \Rightarrow &E[Y_i] \leq E[X_i] \,\,\,\, \text{( Let } Z = X_i-Y_i\,\,\,\,,\,\,\,\, E[Z] = E[X_i]-E[Y_i] \geq 0 \text{)}\\ &E[Y_i] \leq E[X_i] = \frac{n}{2^i} \,\,\,\, \text{ (Apply Lemma-2)} \\ &Y = Y_0 = \underbrace{1 + 1 + 1 \dots + 1}_{t \text{-times}} + Y_{t} \\ \Rightarrow & E[Y] = \underbrace{E[1] + E[1] \dots + E[1]}_{t \text{-times}} + E[Y_{t}] = t + E[Y_{t}] \\ &E[Y] = t + E[Y_{t}] \leq t + E[X_{t}] \leq t + \frac{n}{2^t} \\ &\\ \Rightarrow &t \leq E[Y] \leq t + \frac{n}{2^t} \,\,\,\, \text{After } t \text{ steps }\\ & \\ &\text{When } t = \log_2(n) \text{ the above inequality becomes the following } \\ & \log_2(n) \leq E[Y] \leq \log_2(n) + \frac{n}{2^{log_2(n)}} = \log_2(n) +1 \\ \Rightarrow & E[Y] =\Theta(log(n)) \,\, \Box \\ \end{array}$$

This completes our expected analysis of the RandSelect algorithm. The analysis just uses elementary probability theory without relying on a Calculus based random walk theorem in the "Randomized Algorithms" book by Motwani and Raghavan -- (Theorem-1.3 on page-15).

## Thursday, November 24, 2016

### Integer Division Algorithm with no Division or Multiplication Operators.

I was working through a number theory problem today and I needed an efficient algorithm to find a largest integer $n$ s.t $a\times n \leq b$ for any given integers $a, b$ (w.l.o.g $a < b$). Clearly $n$ is the quotient when $b$ is divided by $a$. Following is an algorithm which can compute the quotient and reminder in $\Theta(\log(a) + \log(b))$ operations without using multiplication or division operators. In fact if the number of bits in $a$ and $b$ are known then the algorithm runs in $\Theta(\log(a) - \log(b))$ operations.

typedef unsigned int uint;

// Integer division of b with a without using division
// or multiplication operators.
div_t integer_div(uint a, uint b) {
div_t result;
uint bit_diff = bit_count(a) - bit_count(b);
uint prod = 0, prod_shift = b << bit_diff, shift = 1 << bit_diff;

result.quot = 0;
do {
if (prod + prod_shift <= a) {
prod += prod_shift;
result.quot += shift;
}
shift >>= 1;
prod_shift >>= 1;
} while (shift);
result.rem = a - prod;

return result;
}

// utility function to find number of bits
// required to represent the value of a.
uint bit_count(uint a) {
unsigned int i = 0;
do {
++i;
} while ((a = a >> 1));
return i;
}

The above algorithm was tested on $4294967293$ pairs of $(a,b), a,b \in [1, INT\_MAX]$ and compared with std::div -- running all these combinations just take around $2$ minutes on a mac book pro. It makes sense because that $\log(a) \leq 32$ (word size on the machine) hence we were just running around four billion constant operations. Typically in all algorithm analysis the word size is considered a $O(1)$, so as long as the inputs bit-widths are within the word size of the machine the runtime of the algorithm can be considered $O(1)$.

Notice that the algorithm is greedy in nature and will provide a correctness next.

Let $c=\sum_{i=0}^{\log(a)-\log(b)}t_i\times 2^{i}$ and $t_{\log(a)-\log(b)}t_{\log(a)-\log(b)-1}\ldots t_0$ the bit representation of $c$. The algorithm tries to pick the powers of $2$ in a greedy manner (i.e. $2^{\log(a)-\log(b)}\geq 2^{\log(a)-\log(b)-1} \geq \ldots 2^{0}$). Consider the choice of $i^{th}$ bit of $c$, if $2^{i}\times a \leq b$ then any $i$-bit integer with a $0$ in the $i^{th}$ bit is strictly less than $2^{i}$. Since $2^{i-1} + 2^{i-2} + \ldots 2^{0} = \frac{2^{i-1+1} -1}{2-1} = 2^{i}-1 < 2^{i}$. We can use induction of $i$ to complete the proof. Hence we can claim that if we greedily pick the powers of $2$ (i.e ${2^i\,\,| \log(a)-\log(b)\leq i \leq 0 }$) it will yield the largest integer $c$ s.t $a\times c \leq b$. The integer division algorithm is just corollary of this basic fact.

## Saturday, June 11, 2016

### What kind of strings have Kolmogorov complexity equal to their length ?

I have been doing some reading about Information Theory and Kolmogorov complexity lately. I just thought about the following problem. Consider a language of all C programs (strings) which print themselves, formally Let $$L_{self} = \{ w\,\, | w\,\, \mbox{is a C program which prints itself} \}$$.
Following is an example of a C program which prints itself (on a side-note the trick behind the construction is to use the ASCII values for quotes and newlines). Clearly the Kolmogorov complexity of each w is less than or equal to the length of w -- given that we have constructed a program of length same as the string w.
I'm wondering if this is indeed optimal ? -- that is given a w, if its possible to produce a program to generate w and size strictly less than |w|.

char qt=34;
char nl=10;
char *str="char qt=34;%cchar nl=10;%cchar *str=%c%s%c;%cint main() {%c printf(str,nl,nl,qt,str,qt,nl,nl,nl,nl);%c}%c";
int main() {
printf(str,nl,nl,qt,str,qt,nl,nl,nl,nl);
}

BTW, compile the program with
gcc -w
to suppress warnings about header inclusion.

## Saturday, March 07, 2015

### Some Special Cases Of the Sequence Assembly Problem.

Sequence Assembly (SA) is a well studied problem in Bioinformatics. To quickly summarize the problem: let $\Sigma$ be some alphabet and $S =\{ s_1 , s_2 \ldots s_n\}$ be a set of strings from $\Sigma$. The goal of the SA problem is re-construct a string $s_a$ such that every $s_i \in S$ is a sub-string in $s_a$. Although its very tempting to reduce the SA problem to the Shortest Super String Problem (which happens to be $\mathcal{NP}$-hard), the repeats in genome makes the SA problem complex and ambiguous.

Any way one of my friends showed me this problem yesterday, which is a much simpler special case of the SA problem. Let $s_{\pi}$ be a some permutation (unknown) of $\Sigma$. Now given a set $S' = \{s'_1,s'_2\ldots s'_m\}$ of sub-sequences of $s_{\pi}$, we want to re-construct this un-known $s_{\pi}$. Just like the SA problem for a given $S'$ there could be several answers for $s_{\pi}$ -- in which case we want to find out the lexicographically smallest permutation. If $\Sigma_i |s'_i| = N$ the problem could be solved efficiently by reducing it to a topological sorting algorithm on a DAG. Following is a quick solution to this problem I wrote while watching the RSA vs PAK world cup cricket match -- I'm little disappointed that RSA lost despite a ferocious knock by AB. The construction of the DAG from $S'$ takes $O(Nlog(N)$ operations, however the cost of topological sorting to report the answer is still $O(N)$ operations.

// 03/06/2015: vamsik (at) engr.uconn.edu//
#include <cmath>
#include <cstdio>
#include <map>
#include <list>
#include <unordered_map>
#include <iostream>
#include <algorithm>
#include <cassert>
using namespace std;

struct CDirectedGraph{
map<unsigned int, map<unsigned int, bool> > m_aEdges;
unordered_map<unsigned int, unsigned int> m_aInDegree;

inline CDirectedGraph() : m_aEdges(), m_aInDegree() {}

inline void add_edge(unsigned int i, unsigned int j){
if((m_aEdges[i]).find(j) == (m_aEdges[i]).end()){
if(m_aInDegree.find(j) == m_aInDegree.end()){
m_aInDegree[j] = 0;
}
m_aInDegree[j]++;
(m_aEdges[i])[j] = true;
}
}

inline void top_sort(){
map<unsigned int, bool> snodes;
unsigned int v;
for(auto vItr=m_aEdges.begin(); vItr!=m_aEdges.end(); ++vItr){
if(!m_aInDegree[vItr->first]){
snodes[vItr->first] = true;
}
}

while(!snodes.empty()){
auto vItr = snodes.begin();
v = vItr->first;
auto eEnd = m_aEdges[v].end();
auto eBeg = m_aEdges[v].begin();

for(auto eItr=eBeg;  eItr!=eEnd; ++eItr){
m_aInDegree[eItr->first]--;
if(!m_aInDegree[eItr->first]){
snodes[eItr->first] = true;
}
}
printf("%u ", v);
snodes.erase(vItr);
}
}
};

int main() {
unsigned int N;
unsigned int K, pnode, cnode;
scanf("%u", &N);
CDirectedGraph diGraph;

for(unsigned int i=0; i<N; i++){
scanf("%u", &K);
scanf("%u", &pnode);
for(unsigned int j=0; j<K-1; j++){
scanf("%u", &cnode);
diGraph.add_edge(pnode, cnode);
pnode = cnode;
}
}
diGraph.top_sort();
return 0;
}

## Sunday, November 16, 2014

### Pseudo Polynomial Dynamic Programming

Sub Set Sum and Integer Knapsack are few examples of combinatorial optimization problems with pseudo polynomial running time algorithms. Often we the implementation of these pseudo polynomial dynamic programming algorithms is done using a huge array to hold intermediate sub problems (e.g. $OPT[K][n]$ would indicate a presence of a subset in the first $n$ integer elements ($x_1,x_2\ldots x_n$) which sum to $K$). If $N$ is the bound on the subset value you were looking for then the space complexity would $\Theta(Nn)$ -- which is clearly exponential on the input size which is $O(n\log(N))$ bits. To build robust algorithms -- reduce the dependence on the value of $N$ -- for these pseudo polynomial algorithms, it would be efficient if we could exploit the underlying DAG (Directed Acyclic Graph) nature of any dynamic programming formulation. But that might be a little bit more code to traverse the sub-problems in a topological order. However we can make our life a little easy by using a hash table to keep track of the sub-problems, which might add to the worst case asymptotic runtime by a factor of $O(log(nN))$. On the other hand the average case of sub-problem lookup would be a constant and would run fairly efficiently in practice.

The following code is a quick illustration of using a hash tables to keep track (and also lookup) of sub-problems in pseudo polynomial dynamic programming algorithms. You can read the problem statement here . Notice the obvious generalization of the problem (scheduling with constant number of resources) from the way I have formulated the dynamic program.

// 11/15/2014: vamsik@engr.uconn.edu//
#include
#include
#include
#include
#include
#include
using namespace std;
typedef std::pair KeyType;
struct KeyTypeHasher{
inline size_t operator()(const KeyType& a) const {
return (a.first)*1729 + a.second;
}
};
typedef std::unordered_map SubProbType;
typedef typename SubProbType::iterator SubProbItr;

bool IsFeasible(const std::vector& A, unsigned int G){
SubProbType sub_problems_even, sub_problems_odd;
SubProbItr itr;

//init//
if(A[0] <= G){
sub_problems_even[KeyType(A[0], 0)] = true;
sub_problems_even[KeyType(0, A[0])] = true;
}

for(size_t i=1; i < A.size(); i++){
SubProbType& curr_problems =
(i%2) ? sub_problems_even : sub_problems_odd;
SubProbType& next_problems =
(i%2) ? sub_problems_odd : sub_problems_even;

if(curr_problems.empty()){
return false;
}
next_problems.clear();
//create new sub-problems in the next level//
for(itr = curr_problems.begin();
itr != curr_problems.end(); ++itr){
const KeyType &a = itr->first;

if( (a.first + A[i]) <= G){
next_problems[ KeyType((a.first+A[i]), a.second)] = true;
}

if( (a.second + A[i]) <= G){
next_problems[ KeyType(a.first, (a.second + A[i]))] = true;
}
}
curr_problems.clear();
}

return ((A.size())%2) ? !(sub_problems_even.empty()) : !(sub_problems_odd.empty());
}

int main() {
/* Enter your code here. Read input from STDIN. Print output to STDOUT */
size_t T,N;
unsigned int G, a;
scanf("%lu", &T);
for(size_t i=0; i < T; i++){
std::vector A;
scanf("%lu %u",&N, &G);
for(size_t i=0; i < N; i++){
scanf("%u", &a);
A.push_back(a);
}
printf("%s\n", IsFeasible(A, G) ? "YES" : "NO");
}
return 0;
}


text:
# Adventures in the Python Visualisation Landscape

Calling the current python visualisation landscape fragmented would probably be an understatement, since itself requires a visualisation to even begin to comprehend. For various reasons I have been unhappy with the tool I was using for visualisation at that time and I have been searching for the one visualisation package to rule them all (spoiler: it doesn't exist)

About a year ago I was using Matplotlib for all my figures, which enabled me to create anything I wanted, usually with a stackoverflow answer giving me a working example to adapt. Problem was, nearly everything required a stackoverflow post or the documentation, making the process of building these figures slow and tedious. Around the same time I was getting fed up with Matplotlib, I found Jake VanderPlas' PyCon talk on The Python Visualisation Landscape which has since acted as both a map and a list of achievements to unlock.

### Altair 1.2

Somewhat naturally the first package I looked at for improving my visualisation workflow was Altair which was at the time a 1.2 release. At the time I thought it was fine for simple figures, however my data was not formatted to make the most of what Altair could offer, which as a relative newbie to pandas was frustrating. Another sticking point with Altair was working out how to customise figures, After Matplotlib, I was expecting to search what I wanted to in Google and have appropriate stack overflow answer on the first page of results. With Altair being a new package, this how-to style documentation didn't exist, which meant I never worked out how to customise the figure. The most frustrating of the default settings was the axis labels which use the SI Prefix for large numbers i.e. 1Mm for $1 \times 10^6$ and 1no for $1 \times 10^{-9}$, something which I never worked out how to change at the time. Another problem I had with Altair (and Matplotlib) was the lack of interactivity. I was in the early stages of data investigation which required both a high level overview of the trends, while also being able to look at smaller regions in more detail. Having to constantly change axis ranges was making it slow and frustrating to create a figure with the appropriate view.

### Bokeh

With an emphasis on interactivity, Bokeh was the tool I was looking for to investigate data on different scales. With the ability to generate figures in a notebook, for quick visualisations of datasets to understand the data, and compare it with previous studies, and also as a Bokeh server application, to perform standard analyses on large volumes of data, like ensuring a simulation is running properly. Despite having these interactive visualisations, bokeh is lacking in the same way as Matplotlib, it takes a long time to specify everything required to create the figures. While there was a start on a high level plotting interface in the form of bkcharts, the code is now unmaintained and directs users to Holoviews.

### Holoviews

Like Altair, Holoviews is a declarative plotting interface, providing a quick and simple interface for constructing figures. Holoviews is structures around the idea of describing data on creation of the dataset, rather than the construction of the figure, allowing for simple figure definitions. Rather than being a library that actually creates a figure, Holoviews performs the reasoning about the dataset and passing the rendering off to Bokeh, Matplotlib, or Plotly. I initially saw this as a big strength of Holoviews, that I could generate interactive visualisations in bokeh, change the output to Matplotlib and have a configurable figure for publication. In practice this is not so simple, with some modifications to plot style parameters changing between the different output formats and having a limited range of customisations.

The drawback of having this special annotated data object is that you lose all the flexibility of having a pandas DataFrame and the vast array of operations that allows. This is important to me because my field of science has a long history of researchers taking some quantities and combining them to give a describable temperature dependence. This lack of flexibility led me back to Altair, this time version 2.0.

### Altair 2.0

Why Altair again? When I first tried Altair I was approaching it from Matplotlib, with its extensive documentation in the form of the technical reference and the numerous how to guides. This time I was approaching it from Bokeh and Holoviews, both of which are relatively new libraries with some teething problems. As a consequence I had become much better at problem solving issues and navigating technical reference materials. Another major moment of discovery was making the connection that since Altair implements the Vega-Lite specification, I should have a look at the Vega-Lite documentation. This turned out to be particularly helpful because the Vega-Lite documentation is currently more extensive than for Altair.

It was in reading the Vega-Lite documentation that I finally understood how to use transform functions in Altair. These are a set of functions that perform computations on the input dataset to generate the resulting figure. This allows me to have a single canonical dataset, with data transformations like ratios of two quantities tied to the figure, rather than following awkwardly named variables around. One example of these functions is to make the data in the cars dataset useful for the majority of the world's population by converting the units as part of the figure definition.

import altair as alt
from vega_datasets import data

cars = data.cars()

chart = alt.Chart(cars).mark_circle().transform_filter(
alt.expr.datum.Miles_per_Gallon > 0
).transform_calculate(
'Fuel Economy (L/100 km)', '235.2 / datum.Miles_per_Gallon'
).transform_calculate(
'Weight (kg)', 'datum.Weight_in_lbs * 0.45'
)
chart.encode(
x='Fuel Economy (L/100 km):Q',
y='Weight (kg):Q',
size='Acceleration:Q'
)

This allows for only storing the fundamental values in the dataset and being able to compute derived values as part of the figure, something that is particularly useful in my workflow.

This computing of values, also extends to the computation of histograms, complete with shortened notation. Using the chart object from above it is possible to easily create a histogram

chart.mark_bar().encode(
x=alt.X('Fuel Economy (L/100 km):Q', bin=True),
y='count():Q',
)

Where setting bin=True will create bins with the default parameters, and the count():Q on the y axis counts the elements in each bin. Instead of count() it is also possible to perform other aggregations, like computing the mean of a column.

chart.mark_bar().encode(
x=alt.X('Fuel Economy (L/100 km):Q', bin=True),
y='mean(Weight (kg)):Q',
)

For a more comprehensive view of using Altiar, have a look at either the Example Gallery, or a case study.

Each of the visualisation libraries in python have their own strengths and weaknesses, types of visualisations they excel at, and others you wouldn't want to try. For me, while Altair does still have a some quirks, most notably in the handling of large datasets, and a somewhat complicated method of setting titles, it provides a simple and intuitive interface to data which at the time of writing makes it the first tool I will reach for to understand a dataset.

format: plain


text:
# A derivative word search

This square is composed of elements, all in the same set, written either horizontally or vertically. There is no overlap, each letter belongs to exactly one word. However, the words are anagrammed. Please find the words, and the specific smallest recognised set which includes these.

STREVIES
CEDRRSLM
ALRNELEE
LOAHWANS
AJFEBTMN
PUAYEEUI
LOTVHOMS
MOBOULCE

• I'm not going to overdo these, just the one before was quickly solved, and want to see if this lasts longer. – Tom Apr 18 at 20:27
• no it doesn't :P but thanks for composing it! – Glorfindel Apr 18 at 20:55
• Over two legs, Glorfindel 2 : 0 Tom, both matches won before 20 mins. – Tom Apr 18 at 20:59
• Tom, it's your turn – Glorfindel Apr 19 at 10:22
• Glorfindel 2 : 1 Tom. An honourable final result, in which the last match was a scrappy affair on foreign soil, which extended over 50 mins. – Tom Apr 19 at 11:52

They all are

SI units, more specifically derived SI units which explains the title, as @oALT notes in the comments

namely

Coulomb (row 8, starting from column 1)
Sievert (row 1, starting from column 2)
Joule (column 2, starting from column 2)
Lumen (column 7, starting from column 2)
Henry (column 4, starting from column 2)
Farad (column 3, starting from column 2)
Weber (column 5, starting from column 2)
Tesla (column 6, starting from column 2)
Ohm (row 7, starting from column 1)
Volt (row 7, starting from column 4)
Pascal (column 1, starting from column 2)
Siemens (column 8, starting from row 2)

• Sorry, hadn't seen your edit, will delete. – hexomino Apr 18 at 21:07
• @hexomino with so many edits, that's completely understandable:) – Glorfindel Apr 18 at 21:08
• doh, I read the 'elements, all in the same set' and was looking for the halogens and so on :P – Beastly Gerbil Apr 18 at 22:39
• (those SI units are also "derived" ones, which explains the title) – oAlt Apr 19 at 3:00