BBMCSP’s overarching ideas are described in the pseudocode below:

Input: A simple graph G(V, E)

Output: A maximum clique in G in Smax

- Compute k-core decomposition of G: K(G)

- Compute an initial clique greedily: Smax

- Remove vertices from G such that K(v)<|Smax|: Gr (Vr, Er)

- Sort Vr according to degeneracy

- Iteratively branch on vertices in Vr in reverse order: v

- | Compute subproblem derived from v: W ⊆ Vr

- | Compute a greedy coloring of W: C(W). if |C(W)|< |Smax| return

- | Compute k-core decomposition of W: K(W). if |K(W)|< |Smax| return

- | Remove vertices w from W such that K(w)<|Smax|: W’

- | Sort W’ according to degeneracy

- | SEARCH for a maximum clique in W’ with BBMC sparse solver: Smax

- |_ Remove v from Vr

- return Smax

All steps except 11 constitute preprocessing tailored for massive graphs. Core hierarchical decomposition in step 1 is a must for any decent graph software package. Each k-core determines a maximal subgraph with the property that all its vertices have degree at least k (cf. previous post). Thus, the core number of a vertex is the order the highest core which contains it and the core number of the graph is the order of its highest core. Note that, in the case of real graphs, the core number is usually much lower than the maximum degree of the graph—e.g. in networks with a power law distribution only a few individuals are massively connected—. Step 3 removes vertices that cannot possibly improve the initial greedy clique in Smax because of a low core number. After that vertices are sorted in step 4 according to degeneracy, i.e. non increasing core number.

Step 5 is the main loop which may be regarded as a first level unrolling of the search tree of a typical exact solver. For each vertex v picked in reverse order, the child subproblem—its neighbor set— is also preprocessed trying to produce an early cut. Step 7 is concerned with approximate color cuts, i.e. the current best solution in Smax cannot be improved by any subproblem in which a coloring of size less than | Smax|. Moreover, Step 8 prunes the subproblem if its core number is less than| Smax|. Finally, if the subproblem cannot be pruned directly in the previous two steps, concrete vertices are removed when their core number is less than | Smax| in step 9. Remaining vertices are sorted by degeneracy in step 10.

The actual NP-hard search occurs in step 11 over the filtered subproblem W'. BBMCSP uses the exact bitstring solver BBMC tailored for large or massive graphs by a specifically designed sparse encoding, as mentioned in the alluded previous post.

Figure 1 contains a simple graph G of order 10 --notation for node labels x:y:x is x for index, y for kcore number and z for degree--. The graph core number K(G) is 4 and its maximum degree is 6 (given by vertex 6). Step 2 of the algorithm finds an initial greedy clique | Smax| of 3 so vertices 1 and 5 are explicitly removed from G in step 3. The remaining vertices are ordered according to degeneracy in step 4 (see figure 1, left).

Figure 2 below contains the resulting behaviour of BBMCSP in the main loop (steps 5-12), which includes the first level unrolling and the actual NP-hard search. Vertices are picked in reverse order (header v) and, after each iteration, removed from Gr. The first iteration in the loop calls vertex 8 and analyses subproblem W={2, 3, 5, 7}. The vertex coloring C(W) has color size 3, same as | Smax| so the subproblem cannot be pruned at step 7. The subproblem is pruned, however, in step 8 by core analysis. Subsequently vertex 8 is removed and vertex 7 startsa new iteration (second row).

The critical step occurs when vertex 6 is chosen, which contains a maximum clique {1, 2, 4, 6}. In this case the entire subproblem cannot obviously be pruned so vertex 9 is reached and vertex 5 is removed since K(5) is lower than | Smax|. The search routine is finally called in step 11 and the clique si found. Remaining vertices in Gr are pruned more or less trivially thereafter.

BBMCSP steps 5-12 taken in the exact maximum clique problem solving

K-core analysis is critical during preprocessing. If the reader wants to implement it by himself we warn him that the typical solution runs in O(|V|2) which is not adequate for large or massive sparse graphs. We recommend the O(|E|) algorithm of Bagatelij—the implementation used by BBMCSP can be found in the pablodev/graph block in the Biicode repo—. A side result is that vertices are also sorted according to degeneracy, which is necessary for steps 4 and 10 as explained.

A sequential greedy coloring in step 7 assigns, to each vertex in turn, the lowest possible color label consistent with the current partial coloring. The bitstring implementation employed by BBMCSP may be found in the pablodev/copt block (InitColor class).

Additionally, the sparse encoding used by BBMC is available as part of the BITSCAN library—pablodev/bitscan block— and the GRAPH library—pablodev/graph— block. BBMCSP uses the sparse_graph available in GRAPH and currently most of the source code is also available in the pablodev/copt block.

In reply to the comments of some of our readers, maximum clique is still NP-hard, yes, and uniform graphs of say a thousand vertices and 0.8 density remain a very difficult challenge for today’s best exact solvers. The reason why BBMCSP is so successful lies mainly on the structure of real networks, which is typically much simplified after the heavy tailored preprocessing described.

That this is so can be seen in the following example: the California road network roadNet-CA with 1.9 million nodes, 2.7 million edges and a clique number of 4 was trivially solved in less than a second of preprocessing by BBMCSP—step 11, the NP-hard search step, is actually never called—. The same graph was proposed as benchmark in a recent Big Data Conference ^{[1]}. There it took 153 seconds to solve using 75 processor nodes!

Raw results of BMCSP against more than 200 real networks, as well as results of state-of-the-art algorithms are available here. We are currently working on a heavy refactoring and a more sophisticated command line parameter interface both for Windows and Linux binaries. As always, we await comments and suggestions from readers.

[1] Hagan, R. D. et al.; Toward an Efficient, Highly Scalable Maximum Clique Solver for Massive Graphs, IEEE Conf. on Big Data, 2014.

]]>I always wanted to learn more about complex numbers. Not just the fact that the root of -1 was equal to **i**. I was particularly interested in figures with fractal structure. I wanted to understand what this meant and how to make such visualization. Having some time from work, I found a book on the theory of functions of a complex variable on the shelf and decided that it was time to start.

Read more]]>

Queues can be different. They can differ in the number of producers and consumers (single/multi producer — single/multi consumer, 4 variants). They can be bounded, on the basis of a pre-allocated buffer, and unbounded, on the basis of a list. They can either support priorities, or not. They can be lock-free, wait-free or lock-based, with the strict (fair) and not strict (unfair) conformance to FIFO, etc. These types of queues are described in detail in the article by Dmitry Vyukov. As a rule, the more specialized requirements to a queue are, the more efficient its algorithm is. In this article, we will consider the most common version of queues – multi-producer/multi-consumer unbounded concurrent queues without the support of priorities.

I guess the queue is the favorite data structure for researchers. On the one hand, it’s really simple; on the other hand, it is not as easy as the stack, as it has two ends, not just one. Since there are two ends, interesting problems arise, such as: how to manage them in a multithreaded environment? The number of publications with different variations of the queue algorithm is off the charts, so it impossible to cover all of them. I’ll dwell briefly on the most popular of them, and start with the classic queue.

Read more]]>

To elucidate my iteration cycle, I am sharing the story of building Rocket Romeo.

*Here’s how the game now looks.*

Like I have said, I have been appreciated to be really creative while building this game. But on the other hand in my opinion — I was not creative, but iterative.

Why I started building an Android game is not the matter of discussion here. So, I will start the story from the point when I felt determined to develop a gaming app. The ideation phase started around 12 Janurary ‘15.

*‘Cursed Bee’. The concept went like this:*

“A Bee has been cursed and can fly clockwise in dark and anticlockwise in light. The game will have hurdles which the Bee has to dodge. The player will TAP to switch between light/dark mode and adjust the flight of bee. ”

Once the basic mechanics was formulated, I focused on graphics building. So I started preparing the Main Screen of the game.

*You can see a bee that traces the dotted path and the honey sticks were the hurdles that it had to dodge.*

*After 3–4 hours of working with graphics, I felt frustrated with the outcome. It seemed repulsive to me.*

I decided to scrape the whole work and start afresh.

But before starting again, I thought why not build a small MVP that would help me in understanding the essence of the game.

So instead of working on the illustration, I jumped on to Android studio and built a quick prototype which looked like this:

*The small square served as a Bee and the red hurdles served as the sticks.*

*I played the MVP for about 1 hr and realized that the mechanics of the game was bad. Not at all interesting and challenging.*

I dropped the idea and decided to move to some other idea.

I really liked the game play of Road Fighter and thought why not build a Multi-player version of the game in which races will be organised at regular interval of time. Players around the global will compete. The race will last for 2 min and the player who travels farther will win. Looked interesting and challenging from development perspective as well.

This time MVP was not important because all I wanted was the exact experience of playing Road Fighter. Seemingly, I started to design the graphics.

The above image shows various iterations of the look ‘n’ feel of the game. After working on this concept for 2 days, I realized that the game would not provide a unique experience. Moreover, there are thousands of better racing games available on Play Store. Another idea dropped.

Now I was out of ideas. I thought about puzzle games, light saber based games, number-summing games and various other forms. But none seemed worth working on.

I have been a fan of Ketchapp’s games. All of them are simple, innovative and fun to play. I wanted to create something simple & addictive which could be built easily by an Indie Developer like me.

Things were not falling in place. All ideas seemed complicated and brazen. But just when I was hopeless, my mind came up with a new idea.

The idea was simple. A falling object with some mechanics that can cause propulsion of the object and slow down its speed. The falling object will cross rectangular openings and points will be awarded for the same.

In no time I was sitting on my laptop with Android Studio and started building a quick MVP. Everything build with rectangular rendering.

*Just to get an understanding of how the game’s experience would be, I built this prototype. The square served as the character that was falling & rectangular blocks were hurdles that were supposed to be crossed.*

*Like the previous MVPs, I played this one for more than an hour and it seemed really interesting with time.*

Once a game developer is ready with the core mechanics, the next job is to put layers of illusion. I needed a story that could cover the reason why the character was falling.

Now, somebody was falling and the tap was resulting in propulsion of the falling person: The story should give the player an impression that the character was on a mission.

It took me around 4–5 hours to write a small story which fit the concept. After multiple iterations, I finalized the following version:

Help Romeo of the Chicken-land to rescue the love of his life, Juliet, from the holds of the Dark Dragon. The Dark Dragon has held Juliet captive inside a volcano. Romeo has got a rocket at his back which will help him descend into the volcano, pass through the various mazes and get to Juliet. Be a superhero and help Romeo find Juliet !

Bingo !

So the background was supposed to be a volcano and the character should have a rocket on his back. The path should have some bridges/hurdles/guards.

I was on a roll at that time. Everything was clear then, all I had to do was build a beautiful Main scene of the game.

After working with Adobe Illustrator for 2–3 days I constructed a few scenes. Following shows the stages of the background & character construction in order of iterations:

__Initial Versions of Romeo and Scene Background.__

Above designs gave me hope and I thought I coulddo better. I started browsing Dribbble for more inspirations. For the next 1–2 days I was browsing Behance and Dribbble, trying to figure the ideal design that can create an amazing gaming experience.

__I finally started to feel better with my design and ended up with these textures for the game. Looked really clean.__

The background was fancy with the exact feel of being inside a cave, transitioning into the end of volcano. Looked amazing to me, so I decided to go ahead.

Next 10 days, I was working on the code, changing atlas to suit the needs of various scene and reading Mario Zechner’s book ‘Android Game Development’. The books is simpley fantabulous. I have used the framework mentioned in it. I wasn’t much aware about gaming engines like Unity3D, AndEngine, Libgdx back then.

In 10 days I had a working copy of the game on my Gingerbread 2.3.2 However, the game was not that interesting as the propulsion was very artificial and the experience was not charming. At this point of time the logic of propulsion was straight forward with normal rocket propulsion theory. The real world effect with rocket was not working here, so I molded the logic and used propulsion theory only for a certain range of velocity. This resulted in a logic which was slightly away from the real life physics, but it created a good experience. [No further iterations were done on the movement of character after this.]

The basic game play was now ready, and the game was in a stage where I could hand it to my friend for feedback. Majority of them liked it but they said that the background was interfering a lot with the game play, which was absolutely true. I never realized that the background was so dominating.

So, I was back on square 1, starting again, designing the background. Re-designing the background meant a redesign of characters & blocks as well.*(By this time, I felt that the green Romeo looked ugly and artificial too)*

*I started experiment with the background again. These are some of the iterations that I tried, but all in vain. All this time I was trying to stick to the story of Romeo falling inside the cave.*

I was not willing to give away the story, it seemed so convincing. But lately I gave up the thought and decided to go for a sky background — simpler and flat.

*The story was put on hold. I had no clue why he was falling down, but I took a leap of faith that I will figure out what he was doing in the story.*

__Final iterations on background & characters with sky background.__

The blue sky seemed clean and promising. So, I decided to go with it. The game was now back on track. I build the main scene, included the Blue guards in it. Added dragons to increase the difficulty as you descend down. To make life easier for player, introduced the concept of collecting eggs. When the player collects 50 eggs, an extra life is awarded.

Apart from the iterations on the Illustration front, I ran a lot of experimentation with character and effects in the game. The game was planned to have teleporting pipes which had random behavior of taking you either forward/backward in the game. The recharge fuels were initially designed in a different way.

*The players were supposed to land on blocks that will contain rocket fuel. However, I dropped it because the experience was not pleasing.*

By the end of the whole development process, I began to think that the game had something missing for people to engage. For simple games like these the most important factor question at the end is ‘How much have other people scored’. Hence, the in-app leader board was constructed. It gives your global rank and tell you where you stand. It drives you to play again and beat others out there.

It took me around 45 days to release the app after bug fixes, testing and all the Hit-n-trials. All these days I was experimenting, testing and iterating. I formed theories, started working with assumptions and tried to put things in place.

The result of all the sleepless is the final **Rocket Romeo**, which is gaining traction at heavy pace. You can find the Application on Play store here.

The whole crux of the story and process taught me a simple lesson that Creativity is all about trying things out and if it doesn’t work try something else. Try to solve problems in different ways.

If you have reached so far, I am grateful to you. I hope my small story will inspire you to work hard and iterate.

If you like this article, please recommend it for the community. If you loved this story, you will love Rocket Romeo too. Download the app from Play Store and don’t forget to Rate & Review it. Here is the link.]]>

(a quote from a chat in Skype):

I read an article about how a dude in the subway fished out a USB flash drive from the outer pocket of some guy’s bag. The USB drive had “128” written on it. He came home, inserted it into his laptop and burnt half of it down. He wrote “129” on the USB drive and now has it in the outer pocket of his bag…

A Picture to Attract Attention:

Read more]]>

I want to show that the dark corners of C are much closer than it seems, and even trivial code lines may contain *undefined behavior*.

The article is organized as a set of questions and answers. All the examples are separate files of the source code.

Read more]]>

These are the main points:

- If you are aware that the market does not provide a decent solution of your real problem, you should think about creating one;

- Graphical user interface is the most important part of the application;

- Even the most simple application can become really popular;

- Even the most simple application can be profitable;

- Do not spend money on reviews in the App Store;

- Advertising in applications is not so bad.

I guess that’s it. You are most welcome under the cut.

Read more]]>

Transducers have become quite popular recently. It’s a new feature of the non-released Clojure 1.7. As of writing this, Сlojure 1.7-alpha5 version is the latest development release, but there also appeared a fair number of ports of transducers for such languages as Python, Ruby, JavaScript, PHP, Java, C++, Lua, Erlang. To tell the truth, it’s a bit disconcerting, the library of reducers was added long ago (in Clojure 1.5). But no one said anything about them; nothing was ported, though they seem to be performing similar things. Or not?

Let’s find out what these reducers & transducers are made for (do we really need them?), how they work and how we can use them… Finally, we will decide whether it’s time to throw reducers away.

It would be wrong to describe concepts emerged in Clojure outside the context of the given language. Therefore, below you will find many listings in Clojure. Plus, there will be no mathematical analysis. Anyway, it’s good to have the basic knowledge of Clojure (especially the concept of **sequences**), but it’s not necessary to know Haskell.

Warning you that all the provided here listings of standard functions are changed a lot, sometimes even “slightly” broken. All for the sake of simplicity. Oh, and on the picture is that very __burrito__.

Read more]]>

Tech blogging is on the rise though most of startups see it as a secondary thing. They’re totally right but there’s one condition to it — it requires skill as well as time. A lot of companies are comfortable with solving this problem by outsourcing content creation and hiring agencies or freelancers.

That’s debatable whether you should really do it or not but let’s focus on one detail that may save you from getting zero engagement, because nobody gives a rat’s ass about your blogs. It’s how you pay for the content you get and we’ll try to explain why exactly pay-per-word just doesn’t work for tech blogging.

Read more]]>

I have not dealt with time series in practice, but I definitely read about them (mostly at school) and had some idea about the way the analysis is carried out. But it is well known that what told in textbooks on statistics and machine learning does not always reflect the real situation.

I guess a lot of people follow the pirouettes made by the curve of oil prices. The chart looks either chaotic, or too regular, so making any predictions on it is quite a thankless job. Of course, we can unleash the full power of statistical, economic and mathematical, and expert methods on time series, but let’s try to deal with the technical analysis – of course, on the basis of R.

When working with regular time series, we can use a standard approach:

Read more]]>

Major tech companies have their own marketing teams, startups tend to hire growth hackers and move towards spending most of their time working on core products. Both are still viable ways of doing marketing but there's a couple of moves that can minimize your co's burn rate. Whether it's a small to mid-sized tech company or a startup, this is going to work for you.

Marketing teams: forget about efficiency here. Low output, high burn rate, poor content quality and tons of press releases – no need to explain the rest. **Reason to this**: dev folks, designers and real workers are off the mic.

There's no need to reinvent the wheel and think up any stories but this is what most of marketing teams do. Do a Skype-call with your tech lead, transcribe it, let him/her edit and finalise your draft. This is it.

We often consult tech co's and startups and see people who have no idea on how to communicate with the real world and talk about their products. They believe that there's nothing interesting to offer, no big things going on at their startup etc.

In-house folks tend to feel good talking about their skills. We rarely hear marketers expose their weak spots, fight for communication with the rest of the company and plan their education. They are too scared to offer any changes and this is when companies lose.

And we definitely love what they have to say about it:

«It's going fine. I'll manage it on my own.»

One of the most staggering things we see is content planning. It's cool to have all your content creation tasks on the plate and we actually help tech co's with their content strategies it but here's what happens in the real world:

**Nothing**.

Everybody is too busy with plans, ideas and disputes. Nobody is into actual content production. Outcome: stories with no expertise, posts smelling like ads, audience skipping through and business standing still (best you can get here).

There's no way you're going to impress people if you're constantly trying to make up something viral, ideal or unique. Put it out there, get feedback and improve. Simple as cake.

Compare the situation with everything «cloud». Think about it as a leverage. Companies and startups utilising Dropbox, Google Apps etc. improve their processes and move in the fast lane while you're still communicating via pigeon post.

Let your team breather and work on core features that hold your business together. Outsourcing social media management, blogging, content strategy and PR-activities is here to stay and be used as the leverage by co's willing to move faster.

Whether you're a small to mid-sized tech company with or without in-house marketing, you should look for hustler's mentality. Hustling is putting every minute and all your effort into achieving the goal at hand. This is what your co's content really needs, so get it right now.

]]>Let’s start from the very beginning. If one event occurrence either increases, or decreases the probability of other event occurrence, such events are called dependent. Probability theory does not study cause-and-effect relations. Therefore, dependent events are not necessarily effects of each other, the relation is not always obvious. For example, “*a person has blue eyes*” and “*a person knows Arabic*” are dependent events, as Arabian people rarely have blue eyes.

Read more]]>

Hustling for attention is hard and doing it for your company’s blog seems even harder but there’s a simple set of tips to be used in case you don’t want to waste time and resources on fruitless blogging.

Read more]]>

Find the Nth Fibonacci Number in O(N) time of arithmetic operations.

Thinking about it, I realized that the only solutions coming to my mind were those operating in O(n) time. But I found a better solution later.

I am going to use the following denotation of sets:

- non-negative integers

— positive integers

According to various mathematic traditions, a set of natural numbers can either contain 0 or not. Therefore, it’s preferred to indicate it explicitly in international mathematic texts.

Donald E. Knuth provides a matrix identity of the following form:

The identity is provided without any proof, but it’s quite easy to prove it by induction.

The matrix on the right is called Q-matrix.

Let’s indicate:

According to the identity, . Thus, to calculate , we should calculate matrix and take the first element of the first line (the enumeration starts with 1).

Since calculation comes to raising the matrix to power, let’s take a look at this process in details.

Let's say there is a matrix to be raised to power. Suppose, is the power of 2. Thus, .

We can represent in the form of a tree:

Meaning that:

So, to calculate matrix, we should calculate matrix and multiply it by itself. To calculate we would have to do the same with and so on.

Obviously, the tree height is .

Let’s estimate calculation time. matrix is of the same size in any power. Therefore, we can perform the multiplication of two matrices in any power in . We should perform of such multiplications. So, complexity is .

There is a question now: what should we do if is not a power of 2? We can factor down any natural number as a sum of numbers that are the power of 2, without any repetitions (we do it all the time when converting a number from a binary to decimal numeral system). I.e.:

where is a set of powers, with the help of which we express . Remembering that is the power of a matrix, we get:

.

Matrix multiplication is not commutative though, i.e. the order of operands during the multiplication is important. But the commutative property is kept for the so-called permutation matrices. Matrix is permutational for , . Thus, we will not have to consider the order of operands during the multiplication task, which simplifies the process a bit.

So, we can represent the algorithm of calculation if the form of the following two steps:

- Factor down into the sum of the powers of two in the form of .
- Calculate all the elements of set.
- Calculate .

Let’s estimate the execution time of the given algorithm.

The first step performs in , where is the number of binary digits in .

Since we have to perform no more than matrix exponentiation, the second step takes time.

The third step performs in , as we should multiply matrices no more than times.

Is it possible to improve the algorithm execution time? Yes, it is. You should note that the second step does not describe the method of calculating matrices that belong to set. We know that is the power of 2 for each of them. Getting back to the picture with the tree, it means that all of them lie in some or other tree tiers and to calculate the greater ones, we should calculate the inferior ones. If we apply Memoization technique and cache all the calculated matrices of every tier, this will allow us to shorten the execution time of the second step up to , as well as the total execution time of a whole algorithm. That’s exactly what we need according to the task.

Let’s get down to coding. I implemented the algorithm in Python due to the following reasons:

- It’s expressive enough to replace the pseudo code

- Its arbitrary-precision arithmetic is clear.

It’s going to look like this:

```
class MatrixFibonacci:
Q = [[1, 1],
[1, 0]]
def __init__(self):
self.__memo = {}
def __multiply_matrices(self, M1, M2):
"""Matrices miltiplication
(the matrices are expected in the form of a list of 2x2 size)."""
a11 = M1[0][0]*M2[0][0] + M1[0][1]*M2[1][0]
a12 = M1[0][0]*M2[0][1] + M1[0][1]*M2[1][1]
a21 = M1[1][0]*M2[0][0] + M1[1][1]*M2[1][0]
a22 = M1[1][0]*M2[0][1] + M1[1][1]*M2[1][1]
r = [[a11, a12], [a21, a22]]
return r
def __get_matrix_power(self, M, p):
"""Matrix exponentiation (it is expected that p that is equal to the power of 2)."""
if p == 1:
return M
if p in self.__memo:
return self.__memo[p]
K = self.__get_matrix_power(M, int(p/2))
R = self.__multiply_matrices(K, K)
self.__memo[p] = R
return R
def get_number(self, n):
"""Getting the nth Fibonacci number
(a non-negative integer number is expected as n)."""
if n == 0:
return 0
if n == 1:
return 1
# Factoring down the passed power into the powers that are equal to the power of 2),
# i.e. 62 = 2^5 + 2^4 + 2^3 + 2^2 + 2^0 = 32 + 16 + 8 + 4 + 1.
powers = [int(pow(2, b))
for (b, d) in enumerate(reversed(bin(n-1)[2:])) if d == '1']
# The same, but less pythonic: http://pastebin.com/h8cKDkHX
matrices = [self.__get_matrix_power(MatrixFibonacci.Q, p)
for p in powers]
while len(matrices) > 1:
M1 = matrices.pop()
M2 = matrices.pop()
R = self.__multiply_matrices(M1, M2)
matrices.append®
return matrices[0][0][0]
mfib = MatrixFibonacci()
for n in range(0, 128):
num = mfib.get_number(n)
print(num)
```

I wanted to compare the execution time of my algorithm with a classical iteration algorithm that is based on the consecutive calculation saving the previous two numbers. I implemented it like this:

```
def get_number(self, n):
if n == 0:
return 0
a = 0
b = 1
c = 1
for i in range(n-1):
c = a + b
a = b
b = c
return c
```

I used the following technique of performance testing. number accepts values. Objects of **MatrixFibonacci** and **IterationFibonacci** (it’s the class that calculates Fibonacci numbers iteratively) classes are created for each . Then, an arbitrary integer is generated 10 000 times.

The test outputs a lot of strings like: n <tab> T1 <tab> T2. These data is used for building a chart.

As you can see from the chart, the matrix algorithm leaves the iterative algorithm behind starting from (its worth noting that will not fit into unsigned 64 bits). I guess, it's possible to say that the algorithm is impractical, though its theoretical speed is quite good.

You will find the complete source code with the benchmark here.

The code for gnuplot to build charts:

```
set terminal png
set output 'fibgr.png'
set grid
set style data linespoints
set title "Calculation of 10000 Fibonacci numbers (Fn)"
set xlabel "Approximate number n"
set ylabel "Time, s"
plot "benchmark.txt" using 1:2 with lines title 'iterative', \
"benchmark.txt" using 1:3 with lines title 'matrix'
```

- Donald E. Knuth. The Art of Computer Programming. Volume 1. Fundamental Algorithms.

C++ programs have been compiled with gcc-4.7.2 in C++11 mode, using online compiler. Programs in Rust have been built with Rust 0.11, using the rust play app.

I know that there are a lot of changes in C++14 as well as in latest Rust release. But I would still be glad to read your feedback. You’re also welcome to share any information about D.

C++ author has not been happy as for templates implementation in the language. He called them «compile-time duck typing» in one of his Lang-NEXT talks.

The thing is that it’s not always clear, what we should use to instantiate a template looking at its declaration. The situation is worsened by monster-like error messages. Let’s try to build the following program:

```
#include <vector>
#include <algorithm>
int main()
{
int a;
std::vector< std::vector <int> > v;
std::vector< std::vector <int> >::const_iterator it = std::find( v.begin(), v.end(), a );
}
```

Just imagine the happiness of a person reading the longest error message.

As for the Rust, templates and their correctness are verified before instantiation. Therefore, there’s a sharp distinction between errors in the template itself (there should be no errors if you’re using a third-party template) and in the place of instantiation. All you should do there is meet requirements to the type that are described in the template:

```
trait Sortable {}
fn sort<T: Sortable>(array: &mut [T]) {}
fn main() {
sort(&mut [1,2,3]);
}
```

This code does not compile due to an obvious reason:

```
demo:5:5: 5:9 error: failed to find an implementation of trait Sortable for int
demo:5 sort(&mut [1,2,3]);
```

There’s a series of problems in C++ that are represented by an undefined behavior, which occurs if you try to use the already freed memory.

An example:

```
int main() {
int *x = new int(1);
delete x;
*x = 0;
}
```

Since there are no commands to deallocate memory, such problems are impossible in Rust. Memory on the stack is «alive» as long as it’s in the scope. Rust compiler automatically frees the memory for you when a pointer goes out of scope. If the memory has been allocated in a heap, a pointer to it (Box<T>) behaves as a regular variable on the stack (it’s deleted when leaving the scope). There’s a reference counter (std::rc::Rc<T>) and a garbage collector (std::gc::Gc<T>) for the data sharing. Both of them are implemented as third-party classes (you can write your own ones).

С++ variant:

```
#include <stdio.h>
int *bar(int *p) {
return p;
}
int* foo(int n) {
return bar(&n);
}
int main() {
int *p1 = foo(1);
int *p2 = foo(2);
printf("%d, %d\n", *p1, *p2);
}
```

At the output:

`2, 2`

Rust variant:

```
fn bar<'a>(p: &'a int) -> &'a int {
return p;
}
fn foo(n: int) -> &int {
bar(&n)
}
fn main() {
let p1 = foo(1);
let p2 = foo(2);
println!("{}, {}", *p1, *p2);
}
```

The compiler returns the following:

```
demo:5:10: 5:11 error: `n` does not live long enough
demo:5 bar(&n)
^
demo:4:24: 6:2 note: reference must be valid for the anonymous lifetime #1 defined on the block at 4:23…
demo:4 fn foo(n: int) -> &int {
demo:5 bar(&n)
demo:6 }
demo:4:24: 6:2 note: ...but borrowed value is only valid for the block at 4:23
demo:4 fn foo(n: int) -> &int {
demo:5 bar(&n)
demo:6 }
```

```
#include <stdio.h>
int minval(int *A, int n) {
int currmin;
for (int i=0; i<n; i++)
if (A[i] < currmin)
currmin = A[i];
return currmin;
}
int main() {
int A[] = {1,2,3};
int min = minval(A,3);
printf("%d\n", min);
}
```

Prints 0. The result is actually ambiguous, though. Here’s the same in Rust:

```
fn minval(A: &[int]) -> int {
let mut currmin;
for a in A.iter() {
if *a < currmin {
currmin = *a;
}
}
currmin
}
fn main() {
let A = [1i,2i,3i];
let min = minval(A.as_slice());
println!("{}", min);
}
```

The code does not compile and returns the following error:

`use of possibly uninitialized variable: `currmin``

More idiomatic (and operable) variant of this function would look like the following:

```
fn minval(A: &[int]) -> int {
A.iter().fold(A[0], |u,&a| {
if a<u {a} else {u}
})
}
```

```
struct A{
int *x;
A(int v): x(new int(v)) {}
~A() {delete x;}
};
int main() {
A a(1), b=a;
}
```

The code is built, but crashes during the runtime:

`*** glibc detected *** demo: double free or corruption (fasttop): 0x0000000000601010 ***`

The same in Rust:

```
struct A{
x: Box<int>
}
impl A {
pub fn new(v: int) -> A {
A{ x: box v }
}
}
impl Drop for A {
fn drop(&mut self) {} //there’s no point in it. It’s provided for the exact copy of C++ code
}
fn main() {
let a = A::new(1);
let _b = a;
}
```

The code builds and runs with no error. There’s no copying as the object does not implement *trait Copy*.

Rust won’t do anything behind your back. Do you want an automatic implementation of *Eq* or *Clone*? Simply add *deriving* to your structure:

```
#[deriving(Clone, Eq, Hash, PartialEq, PartialOrd, Ord, Show)]
struct A{
x: Box<int>
}
```

```
#include <stdio.h>
struct X { int a, b; };
void swap_from(X& x, const X& y) {
x.a = y.b; x.b = y.a;
}
int main() {
X x = {1,2};
swap_from(x,x);
printf("%d,%d\n", x.a, x.b);
}
```

Output:

`2,2`

The function doesn’t explicitly expect to receive references to the same object. There’s restrict in C99 to persuade the compiler that references are unique. It serves as a hint for the optimizer and does not guarantee that there will be no overlaps: the program will be compiled and executed as before.

Let’s try to do the same in Rust:

```
struct X { pub a: int, pub b: int }
fn swap_from(x: &mut X, y: &X) {
x.a = y.b; x.b = y.a;
}
fn main() {
let mut x = X{a:1, b:2};
swap_from(&mut x, &x);
}
```

Returns the following:

```
demo:7:24: 7:25 error: cannot borrow `x` as immutable because it is also borrowed as mutable
demo:7 swap_from(&mut x, &x);
^
demo:7:20: 7:21 note: previous borrow of `x` occurs here; the mutable borrow prevents subsequent moves, borrows, or modification of `x` until the borrow ends
demo:7 swap_from(&mut x, &x);
^
demo:7:26: 7:26 note: previous borrow ends here
demo:7 swap_from(&mut x, &x);
```

As you can see, the compiler does not allow to reference to the same variable via "&mut" and "&" at once. Therefore, it guarantees that no one will be able to read or modify the mutable variable as long as &mut reference is active. These guarantees are calculated during the building process and do not slow down the execution of the program itself. Moreover, this code compiles as is we used **restrict** pointers in C99 (Rust provides LLVM information as for references uniqueness). It gives the optimizer a green light.

```
#include <vector>
int main() {
std::vector<int> v;
v.push_back(1);
v.push_back(2);
for(std::vector<int>::const_iterator it=v.begin(); it!=v.end(); ++it) {
if (*it < 5)
v.push_back(5-*it);
}
}
```

The code is built without errors, but crashes during the runtime.

`Segmentation fault (core dumped)`

Let’s try to convert it to Rust:

```
fn main() {
let mut v: Vec<int> = Vec::new();
v.push(1);
v.push(2);
for x in v.iter() {
if *x < 5 {
v.push(5-*x);
}
}
}
```

The compiler doesn’t allow to run it by pointing out that we can’t change the vector whilst traversing it

```
demo:7:13: 7:14 error: cannot borrow `v` as mutable because it is also borrowed as immutable
demo:7 v.push(5-*x);
^
demo:5:14: 5:15 note: previous borrow of `v` occurs here; the immutable borrow prevents subsequent moves or mutable borrows of `v` until the borrow ends
demo:5 for x in v.iter() {
^
demo:10:2: 10:2 note: previous borrow ends here
demo:5 for x in v.iter() {
demo:6 if *x < 5 {
demo:7 v.push(5-*x);
demo:8 }
demo:9 }
demo:10 }
```

```
#include <stdio.h>
enum {RED, BLUE, GRAY, UNKNOWN} color = GRAY;
int main() {
int x;
switch(color) {
case GRAY: x=1;
case RED:
case BLUE: x=2;
}
printf("%d", x);
}
```

Returns 2. In Rust you’re obliged to enumerate all variants when pattern matching. Besides, the code doesn’t automatically jump to the next case, unless it meets break. A proper reaction in Rust will look like the following:

```
enum Color {RED, BLUE, GRAY, UNKNOWN}
fn main() {
let color = GRAY;
let x = match color {
GRAY => 1,
RED | BLUE => 2,
_ => 3,
};
println!("{}", x);
}
```

```
int main() {
int pixels = 1;
for (int j=0; j<5; j++);
pixels++;
}
```

In Rust, bodies of loop structures must be wrapped in braces. It's a simple thing, but helps to reduce a set of typos.

```
#include <stdio.h>
#include <pthread.h>
#include <unistd.h>
class Resource {
int *value;
public:
Resource(): value(NULL) {}
~Resource() {delete value;}
int *acquire() {
if (!value) {
value = new int(0);
}
return value;
}
};
void* function(void *param) {
int *value = ((Resource*)param)->acquire();
printf("resource: %p\n", (void*)value);
return value;
}
int main() {
Resource res;
for (int i=0; i<5; ++i) {
pthread_t pt;
pthread_create(&pt, NULL, function, &res);
}
//sleep(10);
printf("done\n");
}
```

Spawns multiple resources instead of just one:

```
done
resource: 0x7f229c0008c0
resource: 0x7f22840008c0
resource: 0x7f228c0008c0
resource: 0x7f22940008c0
resource: 0x7f227c0008c0
```

It is a common problem of thread synchronization. It occurs when the object is changed concurrently by several threads. Let’s try to write the same in Rust:

```
struct Resource {
value: Option<int>,
}
impl Resource {
pub fn new() -> Resource {
Resource{ value: None }
}
pub fn acquire<'a>(&'a mut self) -> &'a int {
if self.value.is_none() {
self.value = Some(1);
}
self.value.get_ref()
}
}
fn main() {
let mut res = Resource::new();
for _ in range(0,5) {
spawn(proc() {
let ptr = res.acquire();
println!("resource {}", ptr)
})
}
}
```

Since we can’t modify the object that is the same for all threads, it throws an error.

```
demo:20:23: 20:26 error: cannot borrow immutable captured outer variable in a proc `res` as mutable
demo:20 let ptr = res.acquire();
```

That’s how the fixed code meeting the compiler requirements can look like:

```
extern crate sync;
use sync::{Arc, RWLock};
struct Resource {
value: Option<Box<int>>,
}
impl Resource {
pub fn new() -> Resource {
Resource{ value: None }
}
pub fn acquire(&mut self) -> *int {
if self.value.is_none() {
self.value = Some(box 1)
}
&**self.value.get_ref() as *int
}
}
fn main() {
let arc_res = Arc::new(RWLock::new(Resource::new()));
for _ in range(0,5) {
let child_res = arc_res.clone();
spawn(proc() {
let ptr = child_res.write().acquire();
println!("resource: {}", ptr)
})
}
}
```

It uses Arc (*Atomically Reference Counted*) and RWLock (to lock shared modifications) the synchronization primitives. We get the following at the output:

```
resource: 0x7ff4b0010378
resource: 0x7ff4b0010378
resource: 0x7ff4b0010378
resource: 0x7ff4b0010378
resource: 0x7ff4b0010378
```

Of course, we can write it properly in C++, as well as in the assembler. Rust doesn’t let you shoot in the foot by protecting you from your own mistakes. As a rule, if the program is built, it operates. It’s better to spend half an hour for adjusting the code to the form that is acceptable by the compiler. Otherwise you’ll have to fix synchronization errors for months (the cost of bugs fixing).

Rust allows playing with raw pointers as much as you like, but not inside unsafe{} block. It’s the case when you ask your compiler not to interfere as you know what you’re doing. For instance, all “foreign” functions (from the library written in C you’re merging with) are automatically marked as unsafe. The language philosophy lies in isolating small parts of the unsafe code from the main part (of the normal code) utilizing secured interfaces. For example, you can find unsafe parts in implementations of Cell and Mutex classes. Isolating the unsafe code allows not only to narrow the search of the occurred problem, but to cover it with tests (we’re on friendly terms with TDD!).

Please let us know what you think and how we can improve. Really! Your opinion is super important. There is a comments section below the article. Otherwise, simply shoot us an email at support@kukuruku.co.

In addition to redesign, we rolled out a new Companies feature.

As you may or may not know, Kukuruku is based on hubs. A hub is a type of weblog in which posts are written and published by more than one author. Every hub is dedicated to a specific topic. It helps us to keep the information categorized and structured.

The main goal of «Companies» functionality is to provide IT companies with an opportunity to have their own corporate hub on Kukuruku.

Whilst regular hubs (such as C++, Haskell, Web Development, etc.) are accessible to any user, corporate hub is a collaborative weblog maintained by company employees only. Our goal is to let companies/startups a simple way to spread a word about their business and stay in touch with professionals from all over the world.

This is **all for free** and we do not charge *ANY* fees.

We are still working on making Kukuruku project more useful, more interesting and more unique. Our team is very small, but we are constantly learning and do our best.

Stay healthy, stay hungry, don't drink too much alcohol and have a great day!

]]>Deep Blue beat Kasparov at chess in 1997.

Watson beat the brightest trivia minds at Jeopardy in 2011.

Can you tell Fido from Mittens in 2013?

The picture and quote are taken from the challenge at Kaggle (contest to benchmark the latest computer vision and deep learning approaches) that took place last year in Autumn.

Looking ahead, we can fairly enough reply “Yes” to the question. The ten leaders coped with the task by 98.8%, which can not but impress.

But where does such question come from? Why have classification tasks been (and still are) beyond programs’ capacity, though a four year old child can cope with them? Why is it more difficult to identify objects of the world than to play chess? What is deep learning and why are cats mentioned so often in related posts? Let’s talk about it.

Suppose there are two categories and a great number of pictures to be apportioned to the two corresponding piles. According to what principle are we going to do that?

Nobody knows, actually, but a common approach is as follows: we will look for some “interesting” pieces of data that can be met in one of the categories only. Such data pieces are called **features**, and the approach itself is called **feature detection**. There are plenty of arguments for the fact that biological brain operates in a similar manner. The first of them is the famous experiment of Hubel & Wiesel on the cells of cat (again) visual cortex.

We never know beforehand, which parts of our picture can be used as good features. Anything can be associated with them: image parts, form, size or color. A feature itself can even be missing at the picture. But it can be expressed in a parameter that has been obtained from the initial data, for example, after the edge detection use. Okay, let’s take a look at some examples with increasing complication:

Let’s assume that we want to make a Google car that can distinguish right from the left turns and twist the steering wheel accordingly. We can compound the rule of identifying a good feature in simple words. Cut off the upper half of the picture and highlight the area of a definite shade (asphalt). Then apply some logarithmic curve to it on the left. If the whole asphalt is located under the curve, it’s the right turn, if otherwise – the left turn. We can have several curves in case of turns of different curvature, and, of course, a set of asphalt shades including dry and wet states. Our feature will be useless on dirt roads, though.

An example from the data set of hand-written MNIST numbers. I guess everyone being familiar with machine learning has seen this picture. Each number has typical geometric entities that define the kind of number it is – a squiggle at the bottom of two, an oblique stroke through the entire field of one, two joined circles of eight, etc. We can create a set of filters that will choose these significant elements, and then apply these filters by turns to our image. The one showing the best result is most likely the right answer.

These filters will look something like this:

It’s the picture from the course by Geoffrey Hinton “Neural networks for machine learning”. Pay attention to 7 and 9 numbers. They miss the lower part. The thing is that it’s the same for seven and nine and contains no useful information for identifying. Therefore, a neuronet that had chosen these features ignored the element. We usually use regular single-layer neuronetworks or something similar to get such feature filters.

Let’s stick to the point. What about this?

There are plenty of differences between these two pictures, such as brightness and color levels. White color prevails in the left part of the left picture. As for the right picture, white color prevails in its right part. But what we need to choose is the differences that will unambiguously detect cats or dogs. For example, the following two pictures can be identified as belonging to one category:

If we look at them for a long time and try to define things they have in common, the only thing coming to mind is the form of their ears. They are more or less the same but a bit bent on the right. But it’s also a coincidence. It’s easy to imagine (and find examples from the data set) a picture, in which the cat looks in a different direction, tilts the head, or it’s captured from the back. All other things differ. The scale, hair length and color, the pose, the background… They have nothing in common. However, a small device in your head can accurately assign these two pictures to one category. As for the upper two pictures, it will identify them as belonging to different categories. As for me, I cannot but admire the fact that such a powerful device is so close to us. Nevertheless, we can not comprehend the way it operates.

Okay. But what if we try and ask an innocent question? How do cats and dogs differ? We can easily begin this list: hair, fluffiness, the form of paws, characteristic poses. Or, for example, the fact that cats have no eyebrows. The problem is that all of these distinctive features are expressed in terms of pixels. We can not use them in an algorithm until we explain to it, what eyebrows are and where they should be located, or what paws are and where they grow from. Moreover, we execute all of detection algorithms in order to understand that it is a cat – a creature we can apply the notions of “whiskers”, “paws” and “tail” to. Without it, we can not even say with reasonable confidence, where the wallpapers or the sofa ends and the cat begins. We’ve come full circle.

And still, we can draw some conclusion from it. When formulating the features in the previous examples, we judged by possible object variability. The road turn can only be either left, or right and there are no other variants (besides going straight, but there’s nothing to do there anyway). Plus, the road building standards guarantee that the turn will be quite smooth, and not at a straight angle. Therefore, we build our feature in such a way so that it would accept various curvature of the turn, a specified set of hues of the road surface and that would be it for the possible variability. The next example: 1 number can be written in various handwritings and all variants will differ from each other. But it should not necessarily contain a straight (or an inclined) vertical line, otherwise it will be stop being one. When preparing the feature filter, we leave some space for the classifier variability.

If we look at the picture one more time, we will see that the active part of the filter for *one* represents a thick line that allows to draw a line with different tilts and a possible acute angle in the upper part.

In case with cats “the space for maneuver” for our objects becomes immeasurably huge. The cats in the picture can be of various breeds, big or small, against any background. They can be partially obscured by some object and they can definitely pose in thousands of different ways. We have not even mentioned rotation and cropping – the bane of all classifiers. It seems impossible to create a flat filter that similar to the previous one, and will be able to take all of these changes into account. Let’s imagine the combination of thousands of different poses in one picture. We will get a shapeless spot that will respond positively to everything. Thus, the sought features should represent a more complex structure. It’s not clear yet, what kind of structure it should be, but it should definitely have a feasibility to consider all these possible changes in it.

The “not clear yet” has lasted for quite a while, — during the bigger part of the deep learning history. Suddenly, people realized the following idea about the world. It sounds something like:

**All things consist of other smaller and more elementary things. **

Saying “all things”, I literally mean everything we are able to learn. Since this post is about the vision, I am talking about depicted outward objects. We can represent any visible object in the form of a composition of some stable elements. The latter consist of geometric figures being a combination of lines and angles located in a certain order. Something like the following:

(Since I haven’t found a good informative picture, I cut this one from Andre Ng (Coursera founder) speech about deep learning.

By, the way, within the limits of blue sky thinking, we could say that our speech and a natural language (that have been considered to be problems of the artificial intelligence for quite a while) represent a structural hierarchy, in which letters form words, and words form word-combinations that form sentences and texts. When coming across a new word, we do not have to learn its letters all over again. As for unknown texts, we do not consider them as something requiring specific memorizing and learning. Checking with the history, there are plenty of approaches (for the most part, they were more scientific) that have already stated it:

1. In 1959 the already mentioned Hubel and Wiesel discovered cells in the visual cortex that reacted to certain symbols on screen. They also discovered existence of other cells “at a higher level” that reacted to certain stable combination of signals from cells of the first level. In this basis, they supposed that there existed an entire hierarchy of similar cells-detectors.

A great video fragment of the experiment that demonstrates the way they discovered the necessary feature by chance. It made the neuron fire after they moved the sample a bit aside and a piece of glass slipped into the projector.

*Warning: If you’re a sensitive person, you’d better not watch the video as they carry out the experiment on a cat. *

2. Somewhere in 2000s there appeared the **deep learning** term with regard to neural networks that have not just one layer of neurons, but plenty of them. Thus, neuronetworks can be taught several feature levels. This architecture has strongly valid advantages. The more levels there are in a network, the more complex structures it can express. But there arises a problem of the way such networks can be taught. The backpropagation algorithm that was commonly used, can not operate well with a big number of layers. There appear several different models for these purposes: autoencoders, restricted Boltzmann machines, etc.

3. In 2004 Jeff Hawkins wrote in his “On Intelligence” book that the hierarchic approach is great and future belongs to it. He was a bit late, but I still want to mention him as he derived this thought in plain language from absolutely ordinary things. He is the person being far from machine learning and even used to say that neural networks were a bad idea. You should definitely read this book as it’s really inspiring.

Thus, we have a hypothesis. Instead of stuffing the learning algorithm with 1024x768 pixels and watch it missing memory and being unable to understand what pixels are important to identify, we will extract some hierarchic structure from the picture. It will consist of different levels. The first one will contain some basic, structurally simple picture elements — its structural bricks: boundaries, strokes and segments. The higher level will contain stable feature combinations of the first level (for example, angles). The next level will contain features built from the previous ones (geometric figures, etc.). But where can we get such a structure for a single picture from?

Now, let’s talk about code.

To represent an object of the real world in the computer, we use some set of rules to convert this object piece by piece to the digital form. For example, one letter can be represented as an integer from ASCII table. As for the picture, it will be divided into plenty of small pixels. Each of them will be expressed with a set of numbers rendering brightness and color information. There are great many of color representing models. It does matter, which one of them we are going to use. But for simplicity, let’s imagine a black-and-white world, in which one pixel is represented with a number from 0 to 1 that will express its brightness – from black to white.

What’s wrong with this representing? Each pixel here is independent. It conveys just a small part of the information about the final picture. On the one hand, it is good when we need to save a picture somewhere or pass it through the network, as it occupies less space. On the other hand, it is uncomfortable for detection. In our case we can see an oblique stroke (with a bit of curve) in the lower part of the picture. It’s hard to tell, but it is a part of the nose contour from the picture of the face. In the given case, important pixels are the ones leaving this strike. The boundary between black and white is important. As for the subtle play of light in light-gray hues of the upper part of the square, it is not important at all. Thus, there’s no need to spend computing resources on it. This view requires dealing with all pixels at once. Each of them is not better than another one.

Now, let’s imagine another code. Separate this square into a linear sum of other squares. Each of them is multiplied by a coefficient. Let’s imagine a lot of dark glass plates with various transparency degrees. There are various strokes on each of the plates: vertical, horizontal, etc. Put the plates on top of each other and form a pile. Adjust transparency in a way to get something similar to our picture. It should not be perfect, but enough for identifying purposes.

The new code consists of functional elements. Each of them says something about the presence of some separate component in the initial square. A component with a vertical stroke has 0.01 coefficient. So, we understand that there’s little “verticality” in the sample. But there’re plenty of oblique strokes in this component (see the first coefficient). If we select components of the new code, its dictionary, in an independent way, we can expect that there will be few non-zero coefficients. Such a code is called **sparse**.

We can see useful features of such representation by the example of one of applications named **denoising autoencoder**. If we divide this picture into small squares of 10x10 size and select appropriate code, we will be able to perform this image denoising with impressive efficiency and convert the noised image to code and restore it (here you will find an example). It shows that code is insensitive to random noise and saves those parts of the image being necessary for object identifying. That’s why we think that there became “less” noise after restoration.

On the other hand, it turns out that the new code is more heavyweight, depending on the number of components. For example, a square of 10x10 pixels can become much heavier.

There is evidence that the visual cortex encodes 14x14 pixels (of 196 size) with the help of around 100,000 neurons.

We have suddenly got hierarchy of the first level that consists of elements of the code dictionary. The elements represent strokes and boundaries. The only thing we have to get is the dictionary.

Let’s use scikit-learn package — the library for machine learning for SciPy (Python). In particular, we will utilize MiniBatchDictionaryLearning class. We chose MiniBatch, as the algorithm will not be for the entire data set at once. It will be applied for small, randomly selected data packs. The process is quite simple and can be written in ten lines of code:

```
from sklearn.decomposition import MiniBatchDictionaryLearning
from sklearn.feature_extraction.image import extract_patches_2d
from sklearn import preprocessing
from scipy.misc import lena
lena = lena() / 256.0 # test image
data = extract_patches_2d(lena, (10, 10), max_patches=1000) # extract a thousand of 10x10 pieces – they are the teaching selection
data = preprocessing.scale(data.reshape(data.shape[0], -1)) # rescaling – move values symmetrically to null; the standard deviation should be equal to 1
learning = MiniBatchDictionaryLearning(n_components=49)
features = learning.fit(data).components_
```

If we drew those parts being located in features, we would get something of the kind:

The output via pylab:

```
import pylab as pl
for i, feature in enumerate(features):
pl.subplot(7, 7, i + 1)
pl.imshow(feature.reshape(10, 10),
cmap=pl.cm.gray_r, interpolation='nearest')
pl.xticks(())
pl.yticks(())
pl.show()
```

That’s where we can stop for a little while and recall why we have started all of this. We wanted to get a set of quite independent “structural bricks” forming the depicted object. In order to gain it, we have cut a great many of small square pieces and applied an algorithm. Then we noticed that it is possible to represent with adequate authenticity degree all of the small square pieces in the form of such elements. Since we face only edges and boundaries on the level of 10x10 pixels, we will get them as a result and all of them will be different.

We can use the encoded representation as a detector. To understand, whether the randomly chosen piece of the picture is an edge or a boundary, we ask scikit to choose the equivalent code. Like this:

```
patch = lena[0:10, 0:10]
code = learning.transform(patch)
```

If any of the code components has quite a big coefficient in contrast to other ones, we know that this signals about the presence of the appropriate vertical, horizontal or any other stroke. If all components are about the same, it means that a single-layer noise or a background is depicted here. And we know that it is of no interest.

But let’s move on. We are going to need a few more transformations.

So, we can express any 10x10 fragment with the help of a chain consisting of 49 numbers. Each of them will denote the coefficient of transparency for the corresponding component of the picture above. Now, let’s write these 49 numbers in the form of a square 7x7 matrix and draw what we’ve got.

Here are two examples for clarity:

On the left is the fragment of the original image. On the right is its encoded representation, in which each pixel is the level of presence of the appropriate component (the lighter, the more). We can notice that the first (upper) fragment has no distinct stroke. It code looks like a merge of everything in a weak pale-gray intensity. The second fragment has one component only. All other are equal to null.

In order to improve the second hierarchy level, let’s take a bigger fragment of the picture (so that it could contain several small ones, say, 30x30). Now let’s cut it into small garments and represent each of them encoded. Then, join all of them and teach another DictionaryLearning on the basis of the data. The logic is simple: if the primary idea is right, the adjoined edges and boundaries should also form stable and repeating combinations.

At first glance, it does not look like something meaningful, but it only seems so. That’s what we get at the second hierarchy level, which we practiced on people faces.

However, the size of a fragment here is bigger — 25x25 instead of 10x10. One of the unpleasant features of this approach is the need to adjust the size of the most "smallest unit of meaning".

There occur some difficulties with drawing the obtained “dictionary”, as the second level is taught on the first level code. Its components will look like a mixed hash of points from the picture above. That’s why we should take a step down and split these components into parts, then “decode” them with the help of the first level. We are not going to review this process in details.

Then we grow the levels on the same principle. For example, here’s the third level. We are going to see something interesting in it:

Each face here is a 160x160 feature. These are the most common viewings: frontal, half-turn, right and left, plus various skin colors. Each feature consists of two layers. First of all, they allow performing a quick check of test images for validity. Secondly, they provide additional “freedom”: contours and edges can deviate from perfect lines. But until they are within the limits of features of their level, they can signal upward about their presence.

Not bad, hah?

Obviously, we haven’t. If we run the script I draw all these sets with, we will see quite a dispiriting situation with the sought data set about cats and dogs. Level after level, it will return us the same features representing a bit curved boundaries.

We fetched one of dogs, but it is a pure coincidence — because similar silhouette met in the sample, for example, twice. If we run the script again, it may not appear.

Our approach suffers due to the same reason we criticized regular feed-forward neuron networks. During the process of learning, DictionaryLearning tries to find some common places and structural components of the chosen picture fragments. In case with faces, it turned out well as they are more or less alike: elongate oval forms with some number of deviations. As for several hierarchy levels, they provide more freedom in this regard. In case with cats, it did not turn out well as it was hard to find two alike silhouettes. The algorithm wouldn’t find anything common between the pictures in the test selection, excluding the first levels, in which we dealt with strokes and boundaries. Fail. It’s another dead end.

Actually, come to think of it, selection with a great number of different cats is good in a way that it covers varieties of breeds, poses, sizes and colors. But it may be not so successful for teaching our brain. After all, we teach by means of frequent repetitions and object study and not by looking through all possible variations. To learn playing piano, we should play scales all the time. It would be great if we would have to just listen to a thousand of classics.

Thus, the first idea is to step aside from variety in a selection and focus on one object in the same scene, but, say, in different positions.

The second idea follows the first one. It has already been mentioned by Jeff Hawkins and also by other people. We should try to benefit from time. After all, variety of forms and poses belonging to one object is temporal. To begin with, we can group successively coming pictures knowing that they depict one and the same cat that poses a bit differently in each of them. It means that we will have to change the teaching selection cardinally and use Youtube videos found by “kitty wakes up” query. We are going to talk about it in the following post.

… at github. Launch it with python train.py myimage.jpg (you can also specify the folder with pictures). There are also additional setting parameters: number of levels, size of fragments, etc. It requires scipy, scikit-learn and matplotlib.

- A Primer on Deep Learning. An informative post with the history of the problem, a short introduction and much more beautiful pictures.

- UFLDL Tutorial. A tutorial by Andrew Ng from Stanford. To get your hands dirty. It contains literally everything you need to get familiar with the way all of it works. The introduction, the process mathematics, parallels with feed-forward networks, illustrated examples and exercises in Matlab/Octave

- Neural Networks and Deep Learning. A free online book. Unfortunately, it’s not finished yet. It describes the basis in simple phrases, starting with perceptrons, neuron models, etc.

- Geoffrey Hinton tells about new generations of neuron networks

- The last talk with Hawkins, in which he covers briefly the ideas of his book, but more specifically. He tells about the things an intelligent algorithm should be able to do, about the fact that well-known features of a human brain tell us about it. He also dwells on neuron networks drawbacks and sparse coding advantages.

Read more]]>

A trie is an ordered tree data structure that is used to store a dynamic set or associative array where the keys are usually strings. We are going to take a look at the implementation of prefix trees for storing ASCII strings. Each string ends with a terminal symbol that will never repeat in a string again. To indicate terminal symbols in pictures and examples, we will use the dollar sign. To indicate them in the code, we’re going to use the null symbol. Thus, in our implementation the strings will be standard C ones, which allows to use the standard string.h library. All the mentioned above should be easily carried over to other alphabets.

In a regular radix tree each edge is assigned with some symbol. As for root nodes, they store some service(user) information. Thus, any route from a tree root to one of its leaves defines precisely only one string. This string is considered to be stored in the specified tree. For instance, a trie stores the following set of strings: {abab$, aba$, bc$, b$, bac$, baca$} (see the picture below). The root is on the left, leaves (their number coincides with the number of strings) are on the right. The root of baca$ string is highlighted in red.

A cautious reader will notice that this way a trie represents a regular deterministic finite-state machine, the states of which corresponds to leaves of a trie. Trie implementation in libdatrie library is based exactly on this concept. The automaton is represented with the help of arrays.

To make an automaton smaller, represent a trie in a partially compressed form. Such form means that all “tails” of the tree (areas with no branches) are packed in strings (see the picture A below). We can also compress in a similar way all the chains with no branches (picture B). In this case the tree will no longer be a finite-state automaton, or rather it will still be an automaton, but for an alphabet of strings, not symbols. Working with a tree like with an automaton relies on the fact that the alphabet size is finite and we can define its serial number by any symbol in the alphabet in O(1) of time (the same will not work for strings). Thus, its implementation with the help of arrays becomes problematic. But when implementing with the help of pointers, there should be no problems with strings storage.

Let’s see how we can use pointers to represent trie edges. First of all, we should realize that in this case, it will be really uncomfortable to store information in trie edges. Therefore, we should relocate symbol chains from edges to root nodes. Use an obvious feature of oriented trees: just one edge enters any node, except for the root one.

We are going to get the following structure:

There’s another problem, as an outdegree of a tree node can be random (even of the current alphabet size). To solve this problem, let’s apply a standard method: store a list of children nodes in each node. The list will be singly linked, so we will store only a head (older sister) in a node. This will allow us to reject an empty root. Now, a trie will be represented with a pointer to the list head of children nodes of the old root (i.e. we replace a tree with a forest). Thus, each node will store precisely two pointers: *link* – to the oldest node, *next* –to its younger sister. The following picture shows the process of such transformation; blue arrows correspond to *link* pointers, red arrows – to *next* pointers.

In order to understand the principle of working with a tree, you should keep in mind the scheme that is presented in the picture on the left. Meanwhile, the actual representation of a tree will be like in the picture on the right.

So, we have moved on from the tree with a variable outdegree (a regular tree) to the *binary tree*, in which *link* and *next* are pointers to the right and left subtrees, like in the picture below:

Now we can get down to implementation. A tree node is a chain of *key* symbols, its *len* length and *link* and *next* pointers. The chain is not required to end with a terminal symbol, so we should know its length. There’s also a minimal constructor that creates a trivial tree consisting of one node with a given key (to copy symbols, use a standard strncpy function) and an even more minimal destructor.

```
struct node // a structure to represent tree nodes
{
char* key;
int len;
node* link;
node* next;
node(char* x, int n) : len(n), link(0), next(0)
{
key = new char[n];
strncpy(key,x,n);
}
~node() { delete[] key; }
};
```

The first operation we will take a look at is inserting a new string in a prefix tree. The search idea is standard. We will move from the tree root. If the root is empty, the search is unsuccessful. If otherwise, we will compare the key in the node with a current x string. Use the following function that will calculate the length of the biggest common prefix of two strings of the given length.

```
int prefix(char* x, int n, char* key, int m) // length of the biggest common prefix of x and key strings
{
for( int k=0; k<n; k++ )
if( k==m || x[k]!=key[k] )
return k;
return n;
}
```

When searching, we are interested in the following cases:

- A common prefix can be empty. In this case we will continue to search recursively at the child node of the given one. i.e. we will pass by
*next*reference;

- If a common prefix is equal to the sought
*x*string, the search is successful, we’ve found the node. That’s when we use the fact that we can find the end of a string, which has a terminal symbol, in the tree leaf only.

- A common prefix is congruent with the
*key*, but is not congruent with*x*. In this case we will recursively go to the parent child node by using a*link*reference and passing it*x*string without the found prefix.

If we have found the common prefix, but it is not congruent with the *key*, the search is unsuccessful as well.

```
node* find(node* t, char* x, int n=0) // x key search in t tree
{
if( !n ) n = strlen(x)+1;
if( !t ) return 0;
int k = prefix(x,n,t->key,t->len);
if( k==0 ) return find(t->next,x,n); // let’s look for the child’s node
if( k==n ) return t;
if( k==t->len ) return find(t->link,x+k,n-k); // let’s look for at the child’s node
return 0;
}
```

By default, *n* length of *x* string is zero so that we could call the search function by specifying a tree root and a search string.

```
...
node* p = find(t,"baca");
...
```

The following picture shows the process of searching three strings (one of them is successful, the other ones – not so much) in the provided above tree.

We should note that in case of success, the search function returns a pointer to the tree leaf, where the search was finished. This node should contain all the information connected with the sought string.

A new key insertion (as well as in binary search tress) is similar to key search. Of course, there are some differences. First of all, in case of an empty tree, we should create a node (a trivial tree) with the defined key and return a pointer to this node. Secondly, if the length of a common prefix of the current *key* and the current *x* string is more than zero, but less than the *key* length (the second case of an unsuccessful search), we should split the current node into two, leave the found prefix in the parent node and place the remaining part of the *key* into *p* child node. To perform this operation, use *split* function. After splitting the node, go on inserting *x* string in a *p* node without a found prefix.

The node splitting code:

```
void split(node* t, int k) // dividing t node according to k key symbol
{
node* p = new node(t->key+k,t->len-k);
p->link = t->link;
t->link = p;
char* a = new char[k];
strncpy(a,t->key,k);
delete[] t->key;
t->key = a;
t->len = k;
}
```

Insertion code:

```
node* insert(node* t, char* x, int n=0) // inserting x key in t tree
{
if( !n ) n = strlen(x)+1;
if( !t ) return new node(x,n);
int k = prefix(x,n,t->key,t->len);
if( k==0 ) t->next = insert(t->next,x,n);
else if( k<n )
{
if( k<t->len ) // cut or not to cut?
split(t,k);
t->link = insert(t->link,x+k,n-k);
}
return t;
}
```

The following picture represents abaca$ and abcd$ keys insertion in the mentioned above tree:

You should note, that if the given *x* string is already in the tree, there will be made no insertion. In this case a prefix tree behaves as a proper set.

As always, key deletion is the most difficult operation. Though in case of tries, it is not that scary. The thing is that when we delete a key, we delete just one leaf node corresponding to some suffix of the key. At first, we should find this node. If the search is successful, we delete it ands return the pointer to the younger sister of the given node (since it is a leaf, it has no children, but can have sisters).

At that, we could finish the deleting process, but there’s one problem. After the node deletion, there can form a chain of *t* and *p* nodes, in which *t* node has just one *p* child node. Therefore, if we want to keep the tree compressed, we should join these two nodes into one by performing merge operation.

The code of merge operation is quite trite. We create a new key, then move the subtree of a *p* node to *t* and delete *p* node.

```
void join(node* t) // t and t->link nodes merge
{
node* p = t->link;
char* a = new char[t->len+p->len];
strncpy(a,t->key,t->len);
strncpy(a+t->len,p->key,p->len);
delete[] t->key;
t->key = a;
t->len += p->len;
t->link = p->link;
delete p;
}
```

The older node is in charge of merging as the younger one does not have the information about its parent. Merge execution has the following criteria:

- Key deletion by
*link*reference, not*next*;

- After deletion, a new
*link*reference does not have*next*reference. There’s just one child node, so we can merge it with the current one.

```
node* remove(node* t, char* x, int n=0) // deleting x key from t tree
{
if( !n ) n = strlen(x)+1;
if( !t ) return 0;
int k = prefix(x,n,t->key,t->len);
if( k==n ) // deleting a leaf
{
znode* p = t->next;
delete t;
return p;
}
if( k==0 ) t->next = remove(t->next, x, n);
else if( k==t->len )
{
t->link = remove(t->link, x+k, n-k);
if( t->link && !t->link->next ) // does t node have just one sister node?
join(t);
}
return t;
}
```

Examples of deleting keys without merging:

with merging:

There have been carried out a computational investigation on comparing AVL and prefix trees as for operation time and memory consumption. As for me, the result was a bit bewildering. Anyway, let’s go step by step. There have been formed 8 test string sets. The following table represents their features:

First of all, they measured the build time of a tree according to the given set of strings and time for searching in it all the keys of the same set (i.e. successful search only). The below graph presents the comparison of AVL tree and a prefix tree (radix tree) build time. You can see that a prefix tree builds up a bit faster.

The following chart compares the time spent on searching all keys for the same two trees. A balanced binary tree spends two times less time on searching, than a prefix one.

Finally, let’s take a look at memory requirements for one symbol. The results are presented in the following graph:

On average, there are about 2 bytes per one symbol, which is not bad. But in case of flags a prefix tree spends less than 1 byte on 1 symbol (there are plenty of common long prefixes).

Thus, there have been found no winner. It’s worth to compare the operation time of these two trees as for the number of strings. But according to the provided above charts, there would be nothing revolutionary. Of course, it would also be interesting to compare these two approaches with hash tables…

Thanks for your attention!

You’re most welcome to share your thoughts and comments.

P.S. Does anyone know where to get a large datasets open to the public?

Read more]]>

According to classical mechanics, time is absolute and invariable in the sense that when moving from one reference system to another one, time frames do not change:

(x, y, z, t) are coordinates in the old frame, while (x', y', z', t') are coordinates in the new one. It is expected that one system moves in-parallel along x axis at v speed.

It’s the so-called Galilean transformation, which is what happens to coordinates when reference systems are changed. In the Galilean sense, there’s just one “time flow” in the Universe and time coordinates are the same for all objects. At that, classical mechanics does not treat the time arrow originality in any way. Moreover, the concept of the flow of time does not work with the Newtonian formulas.

By the way, we are the ones to introduce the movement from past to future in classical mechanics. Let’s assume that there is a given set of material points (coordinates and speed) and operating forces. Then we provide *dt* interval and watch the system evolve in time. We can also move in the opposite direction and see what happened to the system in the past.

But “time travels”, meaning moving from one particular object to the past along *t* axis, is forbidden by Newtonian mechanics (see above – there’s just one time flow in the Universe).

The situation had changed, when Maxwell introduced electrodynamics and then Einstein created the Special relativity (SR), when trying to solve the contradictions between electrodynamics and classical mechanics.

Within the limits of SR when moving from one inertial reference system to another one, we are going to use the following transformations:

where

As you can see, each reference system has “its own” time scale. Within the limits of SR, neither space, nor time intervals between points are saved. Only the difference of their squares is.

ds^{2} = c^{2}dt^{2} — dx^{2} — dy^{2} — dz^{2}

What is written here? It’s a description of the special four-dimensional pseudo-Euclidian space. It’s the so-called Minkowski space, in which the distance between the points is equal to the difference of squares of the coordinate differences (not the sum, as in the regular Euclidian space).

From the point of view of SR, each body movement represents a trajectory in the four-dimensional Minkowski space, the so-called “world line”.

Cones at the picture are the so-called “light cones”. They are trajectories of objects moving at the speed of light.

Since each object has its own time that can change arbitrary, we can theoretically assume that there can be “unfolding” of the world line, meaning time traveling.

In the special theory of relativity (SR), it’s quite simple to travel forward in time. It’s enough to move at the speed that is close to the speed of light. This effect has been confirmed successfully. Therefore, we are going to talk about traveling back in time.

The first one lies in the fact that a material object can never reach the boundary of the light cone, as it would be necessary to apply an infinite amount of energy for that purpose.

The second one lies in peculiarities of Lorentz transformations. Suppose, we have managed to move an object from A point to B faster than the speed of light. Then there will surely be a reference system, in which the object will show up in B point at first and then will leave from A point. Put simply, movement faster than the speed of light is actually traveling in time in SR, but in an extremely nontrivial sense.

But how should we understand the phrase “there will be a system…”? It appears that for a part of observers there are travels back in time, but there aren’t any for another part. It depends on the speed of the observer movements.

Secondly, there’s a complete mess with cause and effect. It is obvious that the object appearing in B point is the effect of it leaving A point. But for some observers the effect takes place later than the cause. But for others it happens earlier. We are going to get back to this question a bit later.

As a result, moving at the speed of light is forbidden in Special Relativity. Or rather, it’s forbidden to transfer information between the objects that are outside the light cone of each other. Otherwise SR faces unsolvable contradictions.

Moving faster than the speed of light is initially permitted, as long as the systems moving faster than light do not exchange any information. For example, according to the theory of Eternal inflation, at the very beginning of the Big Bang the family of newly-born Universes moves away from each other at the speed that exceeds the speed of light. There’s no contradiction here, as all Universes move away from each other faster than light. The light of any of them will never reach another one and the information will never be transferred.

But let’s get back to our talk about the nature of time. To tell the truth, SR did not bring much understanding to this problem. What is “internal time” of an object and how do we measure it? Why is the temporal constant of Minkowski space opposed to spatial ones? Finally, what does the “inertial reference system” mean and where can we get it in reality?

Einstein was the one to understand all of the mentioned difficulties. Therefore, he had invented General relativity that confused everyone even more.

To begin with, all coordinates in General relativity (including temporal one) are just numbers of points on the axis. Since they have no other meaning, you can set any numbers.

It is necessary to determine the so-called fundamental metric tensor — g_{ab} in the coordinates for each point of the space-time. What is it? It’s a 4х4 matrix that defines the distance between neighboring points. Suppose you are in some (x^{0}, x^{1}, x^{2}, x^{3}) point and move towards (dx^{0}, dx^{1}, dx^{2}, dx^{3}). You are going to cover the distance like this:

Hereinafter, repeated indices mean summation, i.e. the sum of 16 components is written here.

For the regular Euclidean space with the orthonormal basis the matrix is equal to one. In Minkowski space g_{00} = 1, g_{11} = g_{22} = g_{33} = -1.

In the General relativity space there can be any components of the metric tensor. The distance between neighboring points is (almost) an arbitrary function of the difference of coordinates. Keep in mind that coordinates are absolutely relative.

If we want to move from one reference system to another one, all tensors, including the fundamental metric one, will undergo the following changes:

Guess what, this part turns out to be the simplest one :). To cut a long story short, you can assign any coordinates, as well as any reference point in General relativity. The only thing you will have to calculate is the fundamental metric tensor. You can number the moments of time as you like.

We seem to have clarified some things about the space-time geometry. Now, let’s move on to the objects in this space. We should form the so-called Stress-energy tensor (Energy-momentum tensor, T_{ab}) from the distribution of impulses and energies of the bodies. There’s a separate tensor for each point of the space-time. It defines “density” of masses and energies.

Finally, we are going to need derivatives of the fundamental tensor for all directions, as well as derivatives (speeds) of all the objects.

Now, we can state the way geometry of space and moves of bodies influence each other.

The geometry of space is measured the following way:

This equation is known as Einstein’s equation. R_{ab} is the so-called Ricci curvature tensor that is expressed via the second partial derivatives of the fundamental metric tensor. (By the way, Grigori Perelman has proved Puankare theorem in his works dedicated to Ricci flow).

What does all of this mean? Knowing the energy-momentum tensor, the metric tensor and its derivatives, we can calculate the second derivatives of the metric tensor. Thus, we are defining the way energy and impulses influence the space geometry.

In order to define the way geometry of space influences the objects moving in it, we direct ds^{2} → min covered distance for each material point. Expanding this condition with the help of, say, Lagrange–d'Alembert principle, we will get expressions for the second object derivatives (“acceleration”). In particular, for non-mass particles (of light) the motion occurs along the geodesic line that is expressed by the following equation:

Г^{a}_{bc} are Christoffel symbols that are expressed via the metric tensor.

There is quite a smart analogy describing the meaning of such a motion in clear terms. Let’s imagine that we have stretched a sheet and put heavy balls in some parts of it. The sheet has weighed down in these places, just like space distorts when there’s some impulse or energy.

Now, let’s roll a light ball on the sheet. It will move in a straight line on the level parts. But when coming close to massive points, it will crook its trajectory, as if heavy balls have “drawn” it. That is the gravity effect caused by the crook of “space” the ball is moving in.

Okay, let’s sum up the way objects move with regard to General relativity:

- we should know the current geometry of space (the metric tensor) and locations of energies and impulses in it (the energy-momentum tensor);

- we should also know the “speeds”, meaning the first partial derivatives of both of them;

- then we will be able to calculate the second derivatives of the metric tensor, as well as the world lines of objects, and make a modeling step.

To tell the truth, it is not that simple to set some initial state of the system. It can turn out that this location is unacceptable. There’s no algorithm allowing to construct the state. The sum of two acceptable states can turn out to be the unacceptable state. There are four accurate analytic solutions of Einstein’s equations: Schwarzschild, Kerr, Reissner–Nordström and Kerr–Newman metrics describing space with one material spinning/non-spinning charged/non-charged point.

Numerical solution of equations is not either possible due to their complexity. We would have to store four 16-digit matrices for each point of the four-dimensional space and solve 32 equations during each step. Considering modern capacities, it would be impossible.

Anyway, let’s get down to more particular questions. Does General relativity forbid the motion with time decrease? It does not. Directing ds^{2} → min, we can get dt < 0 in some conditions for the optimal motion trajectory. To tell the truth, it’s not simple to formalize “some conditions”. For instance, we can suppose that there exist such a phenomenon as “wormholes”. They are peculiar distortions of the space geometry that allow to travel (perhaps, safely even for humans) at speeds faster than the light speed (which is, in some way, traveling back in time).

Has General relativity clarified the nature of time? On the one hand, space, time, energy and impulse have acquired new meanings within the limits of the theory. On the other hand, questions to General relativity are not clearly formed, just like in Newtonian mechanics. For instance, all the General relativity is covariant. It is written with the help of tensors with subscripts (tensors with superscripts are called contravariant). By the way, I forgot to mention that we cannot change the location of an index for no reason as it will change the expression for transfer from one coordinate system to another one. But the reason General relativity is contravariant is one of philosophical questions that are not simple to realize, let alone answering them.

Taking a look at such modern theories as String ones, we can say that they are more or less about the same. But they operate not on the four-dimensional, but on the 10- or even 11- dimensional space, in which all dimensions, except four, “folded” so that their size is null for us. Multidimensional surfaces are considered in this space (branes are multidimensional membranes). String-objects move on the surfaces. An object moves in a way, so that its trajectory would not “sweep” the minimal “area” on the “brane” surface. The sense is approximately the same, but there are more dimensions.

We could also talk about Quantum mechanics, but, to tell the truth, it’s beyond me. There’s almost nothing interesting about it there. Quantum mechanics is more classical than General relativity in this regard. Its laws are obviously written like derivatives with time. Besides, the temporal coordinate of an object, just like the spatial one, subordinates to the indeterminacy principle and we cannot determine it accurately together with energy.

Now, let’s get back to the causality principle. We can see that General relativity does not provide any insight here. As if there’s no causality. But we do know that it exists and can see it in the world provided in sensations.

General Theory of Relativity (and String theories) have complicated everything even more. While Newtonian mechanics did not allow any time travels, General relativity states that the causality principle is the only thing preventing from traveling in time. Nothing forbids time travels mathematically.

To get some scientific reply to this question, we would have to understand the meaning of the causality principle.

There’s a theory deriving from causality as from the basis of the Universe. It’s the so-called theory of Causal dynamical triangulation. You can google it. Alas, there isn’t anything special in the theory, string theories seem more interesting. But some physicists think otherwise. For example, Lee Smolin in his “The trouble with physics: the rise of string theory, the fall of a science, and what comes next”.

So, that’s it. After all, I have disappointed you. I would be really surprised to find out that anything became clearer for you. For instance, the way time is built and whether time travels are possible. To tell the truth, I became more confused when writing this article. But I do hope that now you realize that it’s impossible to solve this question easily by dreaming up some nonsense :).

Read more]]>

First of all, an AVL Tree is a Binary Search Tree (BST), the keys of which meet standard requirements: a key of any tree node is not less than the key in the left subtree of the given node and not more than any key in the right subtree of this node. This means that in order to find the necessary key in an AVL tree, we can use a standard algorithm. For simplicity, we will consider that all keys in a tree are integers and not repeated.

AVL trees are peculiar as for their balance in a way that *for any tree node the height of its right subtree differs from the left subtree height for not more than one*. It has been proved that this property is enough for the tree height to depend logarithmically in the number of its nodes. **h** height of an AVL tree with **n** keys lies in the range from log2(n + 1) to 1.44 log2(n + 2) − 0.328. Since all major operations on a BST (search, insertion and deletion) depend linearly on its height, we get a guaranteed logarithmic dependence of these algorithms operation time from the number of keys that are stored in a tree. Reminding you, that randomized search trees provide the balance in probability sense only. Though the probability of getting a badly imbalanced tree having high **n** values, is negligible, it is still not equal to zero.

We will represent AVL tree nodes with the help of the following structure:

```
struct node // the structure for representing tree nodes
{
int key;
unsigned char height;
node* left;
node* right;
node(int k) { key = k; left = right = 0; height = 1; }
};
```

**key** field stores the node key, **height** field stores the height of the subtree with the root in the given node, **left** and **right** fields are pointers to the left and the right subtrees. A simple constructor creates a new node (of 1 height) with the specified **k** key.

Traditionally, AVL tree nodes do not store the height, but the difference between the height of the left and the right subtrees (the so-called *balance factor*) that can accept three values only: -1, 0 and 1. But we should note that this difference is still stored in a variable, the size of which is equal to at least 1 byte (if we do not think out any schemes of an “efficient” package of such sizes). Let’s recall that h < 1.44 log2(n + 2). It can mean that when n=10^9 (one milliard of keys, more than 10 GB of memory allocated for nodes storage) the tree height will not exceed h=44 size that can be easily placed in the same 1 memory byte, just as the *balance factor*. Thus, height storage does not increase the memory amount that is allocated for the tree nodes. On the other hand, it efficiently simplifies implementation of some operations.

Let’s define three helper functions related to height. The first one is a wrapper for **height** field. It can operate with NULL pointers (empty trees) as well:

```
unsigned char height(node* p)
{
return p ? p->height : 0;
}
```

The second one calculates the *balance factor* of the given node. It operates with nonzero pointers only:

```
int bfactor(node* p)
{
return height(p->right) - height(p->left);
}
```

The third function retrieves the correct value of **height** field of the given node (provided that this field value in the left and the right child nodes are correct)

```
void fixheight(node* p)
{
unsigned char hl = height(p->left);
unsigned char hr = height(p->right);
p->height = (hl>hr ? hl : hr) + 1;
}
```

We should note that all of the three functions are nonrecursive, i.e. their operation time is О(1).

When adding or deleting nodes in an AVL tree, the balance factor of some nodes can be equal either to 2, or -2. Thus, the subtree is imbalanced. To solve this problem, we should apply the well known rotations round some tree nodes. Reminding you that a simple right (left) rotation causes the following transformations of the tree:

Let’s take a look at the code implementing a right rotation (as usual, each function which modifies a tree, returns a root of a new tree):

```
node* rotateright(node* p) // the right rotation round p
{
node* q = p->left;
p->left = q->right;
q->right = p;
fixheight(p);
fixheight(q);
return q;
}
```

The left rotation is the symmetric copy of the right one:

```
node* rotateleft(node* q) // the left rotation round q
{
node* p = q->right;
q->right = p->left;
p->left = q;
fixheight(q);
fixheight(p);
return p;
}
```

Let’s review the situation of imbalance, when the height of the right subtree of **p** node is greater by 2, than of the left subtree (the opposite case is symmetric and is implemented in much the same way). Let’s assume that **q** is the right child node of a **p** node, while **s** is the left child node of a **q** node.

The analysis of possible cases within the limits of the given situation shows that in order to get rid of imbalance in **p** node, it’s enough to perform either a simple rotation to the left round **p**, or the so-called big rotation to the left round the same **p**. A simple rotation is performed when the height of the left subtree of **g** node is more than the height of its right subtree: h(s)≤h(D).

We can apply a big rotation when h(s)>h(D) and perform two simple ones – at first the right rotation round **q**, then the left one round **p**.

The balance executing code comes to checking conditions and performing rotations:

```
node* balance(node* p) // p node balance
{
fixheight(p);
if( bfactor(p)==2 )
{
if( bfactor(p->right) < 0 )
p->right = rotateright(p->right);
return rotateleft(p);
}
if( bfactor(p)==-2 )
{
if( bfactor(p->left) > 0 )
p->left = rotateleft(p->left);
return rotateright(p);
}
return p; // no balance needed
}
```

The described rotation and balance functions contain neither loops, nor recursion. Thus, they’re executed in constant-time that does not depend on the AVL tree size.

A new key insertion in an AVL tree is executed in much the same way as in simple search trees. We go down the tree choosing either the right, or the left move direction depending on the result of comparing the key in the current node with the inserted one.

The only difference is that we balance current node when returning from recursion (i.e. after the key is inserted either in the left or the right subtree and this tree is balanced). It is proven that during insertion, the imbalance occurring in any node along the path is not more than 2. Thus, the use of the mentioned above balance function is correct.

```
node* insert(node* p, int k) // k key insertion in the tree with p root
{
if( !p ) return new node(k);
if( k<p->key )
p->left = insert(p->left,k);
else
p->right = insert(p->right,k);
return balance(p);
}
```

To check the correspondence of the implemented insertion algorithm to theoretical estimates for AVL trees height, a simple calculating experiment has been carried out. There was an array with randomly generated numbers from 1 to 1000. Then these numbers were subsequently inserted in an initially empty AVL tree and the tree height was measured after each insertion. The obtained results were averaged according to 1000 calculations. The following chart represents the dependency of **n** on average height (the red line), the minimal height (the green line) and the maximum height (the blue line). It also depicts the upper and the lower theoretical estimates.

We can see that the experimentally found heights for random key sequences are within the theoretical limits even with some spare. We can reach the lower bound (at least in some points) if the initial key sequence is arranged in ascending order.

As for nodes deletion from an AVL tree, it is not as simple as with Randomized Search Trees. I have not been able to either find, or create a method that would be based on joining two trees together. Therefore, as a basis I used the variant that is described almost everywhere. It is usually applied for deleting nodes from a standard Binary Search Tree as well.

Its main concept is that we find **p** node with the specified **k** key (if it is not there, we do not have to perform any actions). In the right subtree find min node with the least key and replace the deleted **p** node with the found min.

There are some nuances may occur during the implementation. First of all, if the found **p** node does not have a right subtree, according to an AVL tree property, this node either should have just one child node (the tree of 1 height) on the left, or it is a list. In both cases we should just delete **p** node and return a pointer to the left child node of **p**.

Assume that now **p** node has the right subtree. Find the minimal key in it. According to Binary Trees, this key is found in the end of the left branch, starting from the tree root. Use the recursive function:

```
node* findmin(node* p) // searching the node with the minimal key in p tree
{
return p->left?findmin(p->left):p;
}
```

Another function will be in charge of deleting the minimal element from the given tree. And again, according to the AVL tree property, the minimal element either has a single node, or it is empty. In both cases we should simply return the pointer to the right node and on our way back (when returning from recursion) perform the balancing. We will not delete the minimal node as it will come in handy later.

```
node* removemin(node* p) // deleting the node with the minimal key from p tree
{
if( p->left==0 )
return p->right;
p->left = removemin(p->left);
return balance(p);
}
```

Now all is ready to implement key deletion from the AVL tree. First, find the necessary node by performing the same actions as when inserting a key.

```
node* remove(node* p, int k) // k key deletion from p tree
{
if( !p ) return 0;
if( k < p->key )
p->left = remove(p->left,k);
else if( k > p->key )
p->right = remove(p->right,k);
```

As soon as we have found the **k** key, turn to plan B: Memorize **q** and **r** roots of the left and the right subtrees of **p** node, and then delete **p**. If the right subtree is empty, return the pointer to the left subtree. If the right subtree is not empty, we should find the minimal **min** element and extract it from there. Hook the **q** to **min** on the left. On the right we will hook what was obtained from **r**. Return **min** after its balancing.

```
else // k == p->key
{
node* q = p->left;
node* r = p->right;
delete p;
if( !r ) return q;
node* min = findmin®;
min->right = removemin®;
min->left = q;
return balance(min);
}
```

When quitting recursion, do not forget to execute balancing:

```
return balance(p);
}
```

That’s pretty much it! The search of the minimal node and its extract can be implemented in one function. At that, we will have to solve a problem (not really difficult one) of returning a pair of pointers from the function. But we can gain some time during one passage along the right subtree.

It is obvious that insertion and deletion functions (and also a simpler one search operation) are performed in time that is proportional to the tree height. During the time of these operations performance we execute the descent from the node to the given node. On each level some fixed number of actions are performed. Since an AVL tree is balanced, its height depends logarithmically on the number of nodes. Thus, the operation time of all three base operations is guaranteed to depend logarithmically on the number of the tree nodes.

Tanks for your time!

Complete code:

```
struct node
{
int key;
unsigned char height;
node* left;
node* right;
node(int k) { key = k; left = right = 0; height = 1; }
};
unsigned char height(node* p)
{
return p?p->height:0;
}
int bfactor(node* p)
{
return height(p->right)-height(p->left);
}
void fixheight(node* p)
{
unsigned char hl = height(p->left);
unsigned char hr = height(p->right);
p->height = (hl>hr?hl:hr)+1;
}
node* rotateright(node* p)
{
node* q = p->left;
p->left = q->right;
q->right = p;
fixheight(p);
fixheight(q);
return q;
}
node* rotateleft(node* q)
{
node* p = q->right;
q->right = p->left;
p->left = q;
fixheight(q);
fixheight(p);
return p;
}
node* balance(node* p) // balancing the p node
{
fixheight(p);
if( bfactor(p)==2 )
{
if( bfactor(p->right) < 0 )
p->right = rotateright(p->right);
return rotateleft(p);
}
if( bfactor(p)==-2 )
{
if( bfactor(p->left) > 0 )
p->left = rotateleft(p->left);
return rotateright(p);
}
return p; // balancing is not required
}
node* insert(node* p, int k) // insert k key in a tree with p root
{
if( !p ) return new node(k);
if( k<p->key )
p->left = insert(p->left,k);
else
p->right = insert(p->right,k);
return balance(p);
}
node* findmin(node* p) // find a node with minimal key in a p tree
{
return p->left?findmin(p->left):p;
}
node* removemin(node* p) // deleting a node with minimal key from a p tree
{
if( p->left==0 )
return p->right;
p->left = removemin(p->left);
return balance(p);
}
node* remove(node* p, int k) // deleting k key from p tree
{
if( !p ) return 0;
if( k < p->key )
p->left = remove(p->left,k);
else if( k > p->key )
p->right = remove(p->right,k);
else // k == p->key
{
node* q = p->left;
node* r = p->right;
delete p;
if( !r ) return q;
node* min = findmin®;
min->right = removemin®;
min->left = q;
return balance(min);
}
return balance(p);
}
```

- B. Pfaff. An Introduction to Binary Search Trees and Balanced Trees —
*libavl library description*

- N. Wirth. Algorithms and Data Structures —
*according to Wirth, balanced trees are AVL trees*

- Thomas H. Cormen and Others. Introduction to Algorithms —
*AVL trees are mentioned in exercises to the chapter about Red-black trees*

- D. Knuth. The Art of Computer Programming —
*6.2.3 section is dedicated to theoretical analysis of AVL trees*

Read more]]>

Let’s take a look at *Universe* data type that is defined the following way:

`data Universe a = Universe [a] a [a]`

It’s a doubly infinite list focusing on some element that we can shift using the following functions:

```
left, right :: Universe a -> Universe a
left (Universe (a:as) x bs) = Universe as a (x:bs)
right (Universe as x (b:bs)) = Universe (x:as) b bs
```

It’s basically a zipper, but we can consider it as a constant C-pointer to the infinite memory area as all increment and decrement operations are applicable to it. But how do we dereference it? Let’s define the function that will get a focused value:

```
extract :: Universe a -> a
extract (Universe _ x _) = x
```

For example, Universe [-1, -2..] 0 [1, 2..] represents all integers. Nevertheless, Universe [0, -1..] 1 [2, 3..] are also integers but with slightly changed context (we point to another element).

To get all powers of 2, apply (2**) function to Universe of integers. It’s quite simple to determine the instance of Fucntor class that follows all the laws:

```
instance Functor Universe where
fmap f (Universe as x bs) = Universe (fmap f as) (f x) (fmap f bs)
-- accordingly
powersOf2 = fmap (2**) (Universe [-1, -2..] 0 [1, 2..])
-- ..0.25, 0.5, 1, 2, 4..
```

In a cellular automaton cell values depend on the values of all other cells of the previous step. Therefore, we can create Universe of all shifts and a rule for their convolution.

```
duplicate :: Universe a -> Universe (Universe a)
duplicate u = Universe (tail $ iterate left u) u (tail $ iterate right u)
```

Convolution rule should be of Universe a -> a type. Thus, a rule example for Universe Bool can be the following:

```
rule :: Universe Bool -> Bool
rule u = lx /= cx
where lx = extract $ left u
cx = extract u
```

Having applied the rule to Universe of all shifts, we get the following state of the automaton:

```
next :: Universe a -> (Universe a -> a) -> Universe a
next u r = fmap r (duplicate u)
-- accordingly
un = Universe (repeat False) True (repeat False) `next` rule
```

We can see that our functions follow the following rules:

```
extract . duplicate = id
fmap extract . duplicate = id
duplicate . duplicate = fmap duplicate . duplicate
```

Therefore, Universe forms a comonad and next function corresponds to (=>>) operator. A comonad is a monad dual. Thus, we can see the following analogies between their operations. For instance, *join* superposes embedded scopes, while *duplicate*, on the contrary, doubles the scope; *return* locates into the scope and *extract* extracts from it, etc.

Now, we can just as well implement a two-dimensional cellular automaton. To begin with, let’s declare a type of the two-dimensional Universe:

`newtype Universe2 a = Universe2 { getUniverse2 :: Universe (Universe a) }`

In Haskell, it’s really simple to apply a function to embedded containers with the help of *fmap* composition. Thus, it’s no problem to write an instance of Functor class for Universe2.

```
instance Functor Universe2 where
fmap f = Universe2 . (fmap . fmap) f . getUniverse2
```

We can make a comonad instance by analogy with a regular Universe. Since Universe2 is just a wrapper, we can define the methods using the current terms.

For example, it’s quite simple to execute *extract* twice. As for *duplicate*, in order to get shifts of embedded scopes, we should define a helper function:

```
instance Comonad Universe2 where
extract = extract . extract . getUniverse2
duplicate = fmap Universe2 . Universe2 . shifted . shifted . getUniverse2
where shifted :: Universe (Universe a) -> Universe (Universe (Universe a))
shifted u = Universe (tail $ iterate (fmap left) u) u (tail $ iterate (fmap right) u)
```

That’s almost it! We just should define the rule and apply it with the help of (=>>). In The Game of Life, a new state of a cell depends on the state of neighboring cells. Thus, let’s define the function of their location:

```
nearest3 :: Universe a -> [a]
nearest3 u = fmap extract [left u, u, right u]
neighbours :: (Universe2 a) -> [a]
neighbours u =
[ nearest3 . extract . left
, pure . extract . left . extract
, pure . extract . right . extract
, nearest3 . extract . right
] >>= ($ getUniverse2 u)
Here’s the rule itself:
data Cell = Dead | Alive
deriving (Eq, Show)
rule :: Universe2 Cell -> Cell
rule u
| nc == 2 = extract u
| nc == 3 = Alive
| otherwise = Dead
where nc = length $ filter (==Alive) (neighbours u)
```

Thus, we can implement any cellular automaton by simply defining *rule* function. Thanks to lazy calculations, we get the infinite field as a present, though it leads to linear memory consumption.

Since we apply the rule to each element of the infinite list, to calculate the cells that have not been referred to, we will have to go through all the previous steps. Therefore, we should keep them in memory.

The source code of both files:

Universe.hs:

```
module Universe where
import Control.Comonad
data Universe a = Universe [a] a [a]
newtype Universe2 a = Universe2 { getUniverse2 :: Universe (Universe a) }
left :: Universe a -> Universe a
left (Universe (a:as) x bs) = Universe as a (x:bs)
right :: Universe a -> Universe a
right (Universe as x (b:bs)) = Universe (x:as) b bs
makeUniverse fl fr x = Universe (tail $ iterate fl x) x (tail $ iterate fr x)
instance Functor Universe where
fmap f (Universe as x bs) = Universe (fmap f as) (f x) (fmap f bs)
instance Comonad Universe where
duplicate = makeUniverse left right
extract (Universe _ x _) = x
takeRange :: (Int, Int) -> Universe a -> [a]
takeRange (a, b) u = take (b-a+1) x
where Universe _ _ x
| a < 0 = iterate left u !! (-a+1)
| otherwise = iterate right u !! (a-1)
instance Functor Universe2 where
fmap f = Universe2 . (fmap . fmap) f . getUniverse2
instance Comonad Universe2 where
extract = extract . extract . getUniverse2
duplicate = fmap Universe2 . Universe2 . shifted . shifted . getUniverse2
where shifted :: Universe (Universe a) -> Universe (Universe (Universe a))
shifted = makeUniverse (fmap left) (fmap right)
takeRange2 :: (Int, Int) -> (Int, Int) -> Universe2 a -> [[a]]
takeRange2 (x0, y0) (x1, y1)
= takeRange (y0, y1)
. fmap (takeRange (x0, x1))
. getUniverse2
```

__Life.hs:__

```
import Control.Comonad
import Control.Applicative
import System.Process (rawSystem)
import Universe
data Cell = Dead | Alive
deriving (Eq, Show)
nearest3 :: Universe a -> [a]
nearest3 u = fmap extract [left u, u, right u]
neighbours :: (Universe2 a) -> [a]
neighbours u =
[ nearest3 . extract . left
, pure . extract . left . extract
, pure . extract . right . extract
, nearest3 . extract . right
] >>= ($ getUniverse2 u)
rule :: Universe2 Cell -> Cell
rule u
| nc == 2 = extract u
| nc == 3 = Alive
| otherwise = Dead
where nc = length $ filter (==Alive) (neighbours u)
renderLife :: Universe2 Cell -> String
renderLife = unlines . map concat . map (map renderCell) . takeRange2 (-7, -7) (20, 20)
where renderCell Alive = "██"
renderCell Dead = " "
fromList :: a -> [a] -> Universe a
fromList d (x:xs) = Universe (repeat d) x (xs ++ repeat d)
fromList2 :: a -> [[a]] -> Universe2 a
fromList2 d = Universe2 . fromList ud . fmap (fromList d)
where ud = Universe (repeat d) d (repeat d)
cells = [ [ Dead, Alive, Dead]
, [Alive, Dead, Dead]
, [Alive, Alive, Alive] ]
main = do
gameLoop $ fromList2 Dead cells
gameLoop :: Universe2 Cell -> IO a
gameLoop u = do
getLine
rawSystem "clear" []
putStr $ renderLife u
gameLoop (u =>> rule)
```

Read more]]>

We can’t control the order in which operations are performed. Therefore, we also can’t control the result of the program execution. Even if we get an error, it won’t be easy to step in the same river twice. I’d like to suggest a recipe of testing a multithread application.

We’ll need the following ingredients: haskell, QuickCheck, some monads, salt/pepper to your taste.

Read more]]>

Read more]]>

Parallel programming is difficult. When using systems with common memory, you can’t go without synchronization of parallel processes/threads access to the common resource (memory). The following are used for it:

- locks (mutex);
- lock-free algorithms (lockless);
- transactional memory.

Transactional memory is a technology of concurrent threads synchronization. It simplifies the parallel programming by extracting instruction groups to atomic transactions. Concurrent threads operate paralleled till they start to modify the same memory chunk. For example, operations of nodes adding to the red/black tree (animation in the heading) can operate in parallel in several threads.

Read more]]>

This post implies a good intro into N2O.

In order to learn about Erlang/OTP Web Framework N2O main features you should visit the github page or the official SynRC website. You will find there all the charts and presentations which you like so much.

In this article I will cover the principles of framework operation.

Version under consideration: N2O: 1.1.0.

Read more]]>

Dear

I would like to share with you my first-hand experience in creating a Website on CppCMS (library-template engine on C++). It can also be named as “help for beginners on CppCMS”.

Read more]]>

Read more]]>

I hope that this article will give a good start for a series of notes about lock-free data structures. I would like to share my experience with community, monitoring and thoughts about what lock-free structures are, how to implement them and whether the concepts of Standard Template Library (STL) containers are suitable for the lock-free containers, and when its worth to apply lock-free data structures.

Read more]]>