Spotting Dodgy Domain Names

These are various approaches employing machine learning to differentiate between good domain names and bad ones. By bad, I mean domains that are used to trick people into thinking they're clicking on a legitimate address (www.goog1e.com, for instance).

Data

The data is the top 1 million domain names as recorded by Cisco. You can find it here.

The data was stripped of the top level domains to remove elements that were not useful.

Then, what is left was one-hot encoded converted to bigrams of characters leading to vectors of length 1444 (that is, 38 x 38 possible ASCII characters). The code for this lives here.

This data set was then split down a 95/5 ratio of training to holdout.

We created test data from this holdout data when we deliberately corrupted it. This was done by changing either a single 'o' to become an '0' or a single 'l' to become a '1'. If there were no such characters to corrupt, the data point was discarded.

Kullback-Leibler Results

We actually use a variant of KL divergence in these results that handles zeros in the data - the Jensen Shannon metric.

In the following histograms, red indicates bad domains and green good ones.

The tables represent the entire holdout ("Good") and the entire test ("Bad") data sets with their Jensen-Shannon metric calculated against the training data.

Note that these metrics are calculated by summing the columns of the data sets. This leads to something of an unrepresentative description of the data since the original is one-hot encoded. Therefore, in any subset of 38 elements of a real vector, only one can be 1 and all the rest must be 0. That is, the elements in a vector for a given domain are not independent of each other.

No normalisation

KL Score Histogram with no normalisation

Note that "+4.914e6" in the bottom right hand corner. Indeed the KL scores are close:

Class	KL Score
Good	4 319 810.24
Bad	4 169 380.40

There's a hair's breadth between them so this is probably going to be hard to differentiate the classes.

L1-Normalise everything

Here, the control group's KL was 6370.17 and the bad domains scored 6370.92 - very close. The histograms unsurprisingly look similar:

KL Score Histogram with all vectors L1-Normalised

Hmm, still not much to work with, so let's try combinations of the two. First:

Normalised Baselines, No Normalisation for Others

In this trial, the baseline is L1-normalised but the other vectors are not.

Class	KL Score
Good	97 285.45
Bad	139 889.62

The histogram for the holdout and bad domains now looks like:

KL Score Histogram with the baseline L1 normalised; all other vectors unnormalised

This is good. There are now two distinct distributions with separate peaks.

L2-normalisation gave very similar KL scores and a graph that looked like:

KL Score Histogram with the baseline L2-normalised; all other vectors unnormalised

Let's try:

Unnormalised Baseline, L1-Normalisation for Others

... and it looks like we're back to square one. The KL scores are amazingly close:

Class	KL Score
Good	4 914 404.73
Bad	4 914 407.38

so, not surprisingly are the distributions of the holdout and test data:

KL Score Histogram with the baseline unnormalized; all other vectors L1-normalized

Again, note that: +4.9144e6 in the bottom right hand corner.

So, we seem to be going backwards.

Aside: Normalise then sum baseline, all others unnormalised

I tried a few other variations like normalising then summing, first with L1:

Class	KL Score
Good	773 262.74	KL Score Histogram with the baseline L2-normalised then summed; all other vectors unnormalised
Bad	704 254.41

then L2:

Class	KL Score
Good	94 506.53	KL Score Histogram with the baseline L1-normalised then summed; all other vectors unnormalised
Bad	83 559.17

But their summed distributions didn't give me a difference in KL scores as good as the "Normalised Baselines, No Normalisation for Others" results, so I stuck with those and I discuss just that distribution in the next section.

The ROC

Running our model against the test data, the ROC looks likes somewhat underwhelming:

ROC for "Normalised Baselines, No Normalisation for Others" KL

Or, in 3d where we can see the threshold value:

where a threshold value of about 10 is the closest the curve comes to the top, left hand corner.

It seems clear that out tool cannot with great confidence determine if a domain name is suspect or not. But then could a human? Which of these domain names would you say are bogus and which are genuine?

mxa-00133b02.gslb.pphosted.com
m7.mwlzwwr.biz
live800plus.jp
lf.wangsu.cmcdn.cdn.10086.cn
x10.mhtjwmxf.com
mailex1.palomar.edu
3gppnetwork.org
mimicromaxfinallb-1513904418.ap-south-1.elb.amazonaws.com
mkt4137.com
modt1thf4yr7dff-yes29yy7h9.stream
jj40.com

This is of course a trick question. They're all genuine URLs that Cisco have logged.

However, if our tool is used as part of a suite of metrics it might identify nefarious activity.

Conclusion

Our tool is definitely better than the monkey score but can we improve it? I have a neural net that looks promising but is computationally very expensive. The KL calculations (and variants of them) are very fast and cheap. I'll compare them to a neural net solution in another post.

Everything you needed to know about Spark Structured Streaming

Back Pressure

Backpressure is defined at Wikipedia in the context of routing "as an algorithm for dynamically routing traffic over a multi-hop network by using congestion gradients."

Note that back pressure within Spark was once an option (see the Spark property spark.streaming.backpressure.enabled). However, it appears that back pressure is not necessary in Spark Structured Streaming from Kafka (StackOverflow):

"Structured Streaming cannot do real backpressure, because, such as, Spark cannot tell other applications to slow down the speed of pushing data into Kafka. 

"If you mean dynamically changing the size of each internal batch in Structured Streaming, then NO. ...Generally, Structured Streaming will try to process data as fast as possible by default. There are options in each source to allow to control the processing rate, such as maxFilesPerTrigger in File source, and maxOffsetsPerTrigger in Kafka source."

In general, Kafka consumers don't need back pressure. Note what the creators of the .NET Kafka client write: "The Kafka consumer will only pull from Kafka as fast as you are handling the messages. If you are forwarding the messages consumed from Kafka onto another queue, simply pause before adding more to that queue if it is full... If you question was just to not poll too fast on consumer side (to avoid taking too much memory), then yes, pause will be ok when available. You can also simply not call Poll when your buffer is full"

Partitions and Parallelism

What does the notion of a DataFrame's mean in the world of streams? "When you retrieve the data at first, the number of partitions will be equal to the number of Kafka partition... If your Kafka topic has only 1 partition, that means that prior to groupByKey, your internal stream will contain a single partition, which won't be parallalized but executed on a single executor. As long as your Kafka partition count is greater than 1, your processing will be parallel. After the shuffle boundary, Spark will re-partition the data to contain the amount of partitions specified by the spark.sql.shuffle.partitions." (StackOverflow)

Many, Small Files

Note that the nature of streaming means many files may be created (at least one per mapper per trigger duration - "interval between checkpoints"). Too many files can swamp the Hadoop Name Node so you may want to curb this. Evo Eftimov (in his blog) talks of ways you can do this. Amongst his ideas, you can increase the trigger time or coalesce them in a batch job. We do the latter but have not got to the tuning stage yet. The option, maxFilesPerTrigger, appeared to make no difference to us.

"Note that when writing DataFrame to Parquet even in “Append Mode”, Spark Streaming does NOT append to already existing parquet files – it simply adds new small parquet files to the same output directory... The columns in parquet are stored sequentially one after another and the next such table data segment has to be ready before beginning to write to a parquet file." [ibid]

So, although you might have set the org.apache.spark.sql.streaming.OutputMode to be Append, the files are not actually appended. A new file is created each trigger time.  If the OutputMode for Parquet is Complete, you'll get "Data source parquet does not support Complete output mode" because the structure of a Parquet file is such that its columns are uninterrupted and lie back-to-back.

"In terms of purging old parquet files you may want to partition the data and then periodically delete old partitions as needed. Otherwise you can't just delete rows if all the data is being written to the same output path." (StackOverflow). Consequently, we partition the incoming streams on time and then delete old directories after they have been collated by a Spark batch job.


Unions and watermarks

Streams can be unioned as they are persisted to the same sink. However, at first I was getting the error "There are [1] sources in the checkpoint offsets and now there are [2] sources requested by the query. Cannot continue." This appeared to be just remarking that what it thought was a single source is actually now 2. Adding withWatermark appears to have fixed it

"A streaming query can have multiple input streams that are unioned or joined together... the global watermark will safely move at the pace of the slowest stream and the query output will be delayed accordingly" (the Spark docs).


Caveat

To summarise our architecture:

1. Reads from a union multiple streams and writes partitioned on a timestamp.
2. Coalesces each partition one-by-one sorting the data on field X as it goes. 

It is currently working fine in QA but has yet to meet production levels of data.

Hands-on with Variational Autoencoders

I've hidden 25 samples in 10 000 that are different to the rest then used a Variational Auto-encoder (VAE) to find them. Here is my account of trying to find the anomalies.

Data

Each sample has 50 points in time (represented by a Long) that may be either bunched together in a few hours or scattered randomly across a calendar year.

By far the single biggest improvement came from normalising this data. Without normalising, the neural net was pretty useless.

So, before normalization, a single sample looks like:

1371877007, 1386677403, 1371954061, 1361641428, 1366151894, 1366819029, 1380334620, 1379574699, 1359865022, 1377141715, 1370407230, 1358989583, 1373813009, 1364038087, 1361247093, 1367920808, 1379825490, 1379755109, 1363559641, 1373945939, ...

and after normalization, it may look something like this:

0.2737,   -0.0451,   -0.6842,    1.6797,   -1.3887,   -0.0844,   -0.6952,    0.9683,    0.7747,    1.6273,   -1.0817,   -0.0380,    1.3321,    0.2864,    0.9135,   -1.3018,    1.0786,    0.0830,   -0.3311,   -1.6751,    1.6270,    1.4007,    0.8983,    ...

Note that the normalized data is (roughly) zero-centred and (very roughly) in the region of -1 to 1. See below for why this is relevant.

Aside: it's really, really important for the data to be reproducible through runs. That is, although the data is random, it must be reproducibly random. I wasted a lot of time being fooled by randomness in the data.


What are VAEs?

"It is an autoencoder that learns a latent variable model for its input data. So instead of letting your neural network learn an arbitrary function, you are learning the parameters of a probability distribution modeling your data. If you sample points from this distribution, you can generate new input data samples: a VAE is a generative model." (Keras Blog)

With vanilla encoders, "If the [latent space] has discontinuities (eg. gaps between clusters) and you sample/generate a variation from there, the decoder will simply generate an unrealistic output, because the decoder has no idea how to deal with that region of the latent space.

"Variational Autoencoders (VAEs) have one fundamentally unique property that separates them from vanilla autoencoders, and it is this property that makes them so useful for generative modeling: their latent spaces are, by design, continuous, allowing easy random sampling and interpolation." (TowardsDataScience)


Tuning

According to Andrej Karpathy:

"The most common hyperparameters in context of Neural Networks include:

• the initial learning rate
• learning rate decay schedule (such as the decay constant)
• regularization strength (L2 penalty, dropout strength)"

But first, let's look at:


Activation Functions

"If you know the outputs have certain bounds, it makes sense to use an activation function that constrains you to those bounds." (StackOverflow)

Given our data, one might think that HARDSIGMOID, SIGMOID, SWISH etc or even TANH would yield the best results (SWISH is just x*sigmoid(x)). Wheras RELU, ELU etc don't model it at all.

From Shruti Jadon on Medium.com

(Graphic from here - EDIT see a more comprehensive list here).

But there is an interesting opinion at TowardsDataScience:

"The question was which one is better to use? 

"Answer to this question is that nowadays we should use ReLu which should only be applied to the hidden layers. And if your model suffers form dead neurons during training we should use leaky ReLu or Maxout function. 

"It’s just that Sigmoid and Tanh should not be used nowadays due to the vanishing Gradient Problem which causes a lots of problems to train,degrades the accuracy and performance of a deep Neural Network Model."

However, I found no difference in accuracy with my VAE using  ELU, LEAKYRELU nor RELU. In fact, playing with the 21 activation functions that came with DL4J, I did not see any variety when applying them to the hidden layers.

I only saw a big difference when using it at the bottleneck layer and in the reconstruction distribution (see below).

Regularization

Setting the L2 regularization parameter gave me the following results

L2	Mean	Accuracy(%)	Standard deviation
10-5	15.2	60.8	0.837
10-4	15.2	60.8	0.837
10-3	15.2	60.8	0.837
10-2	15.2	60.8	0.837
10-1	16	64	0.707
100	16.2	64.8	0.447
101	16	64	0
102	16	64	0

All using the SWISH activation function.

Batch Sizes

Accuracy hovered around 16 or 17 up to and including a batch size of 64. After that, it dropped off quickly to an accuracy of 13 (53%) and 6 (24%) for batch sizes of 128 and 256.


Updater

Adam with an initial value of 10^-4 seemed to give better accuracy at 17.8 / 71.2% (sd. 0.422) than RmsProp(10^-3) and AdaDelta (which both yielded an accuracy of 16 (64%), standard deviation of 0).


Reconstruction Distribution

Now fiddling with this knob did make quite a difference.

All the results so far were using a BernoulliReconstructionDistribution with a SIGMOID. This was because I had cribbed the code from somewhere else where the Bernoulli distribution was more appropriate as it represents "binary or 0 to 1 data only".

My data was not best approximated by a Bernoulli but a Gaussian. So, using a GaussianReconstructionDistribution with a TANH gave better results.

The DL4J JavaDocs state: "For activation functions, identity and perhaps tanh are typical - though tanh (unlike identity) implies a minimum/maximum possible value for mean and log variance. Asymmetric activation functions such as sigmoid or relu should be avoided". However, I didn't find SIGMOID or RELU made much difference to my data/ANN combination (although using CUBE led to zero anomalies being found).

This is similar to what I blogged last year that when modelling the features: features should (in a very loose sense) model your output.

Anyway, using a Gaussian reconstruction distribution, accuracy jumped to 18.6 (74.4%) albeit with a large standard deviation of 3.438.

Then, through brute force, I discovered that using SOFTPLUS in both pzxActivationFunction and GaussianReconstructionDistribution gave me an average accuracy 19.1 (sd. 3.542). This was the high-water marker of my investigation.


Architecture

All the results so far were using just a single encoder and a single decoder layer that was half the size of the input vector. Let's call this value x.

Using hidden layers of size [x, x, 2, 2, x, x] did not change the best accuracy. Neither did [x, x/2, 2, 2, x/2, x] nor [x, x/2, x/4, 2, 2, x/4, x/2, x] nor even [x, x/2, x/4, 1, 1, x/4, x/2, x].

So, this avenue proved fruitless.

Conclusion

I am still something of a neophyte to neural nets but although I can improve the accuracy it still seems more like guesswork than following a process. There was no a priori way I know of that would have indicated that SOFTPLUS was the best activation function to use in the reconstruction, for instance.

It's clear that there are some rules-of-thumb but I wish somebody would publish a full list. Even then, it seems very data-dependent. "For most data sets only a few of the hyper-parameters really matter, but [...] different hyper-parameters are important on different data sets" (Random Search for Hyper-Parameter Optimization).

Statistical Covariance

The meaning of statistical covariance is completely different to that used in programming languages but it's very important to data scientists. Here are some miscellaneous notes I made playing around.

The covariance matrix

... can be calculated for matrix X with something like this:

    x = X - np.mean(X, axis = 0)
    C = np.dot(x, x.T) / (n - 1)

For x, we're just taking the means of each column and subtracting each mean from each element in its respective column.

Note that in the diagram above, the three vectors from the covariance all live in the same plane. This is not a coincidence. The covariance matrix “is ALWAYS singular for ANY square matrix because you subtract off the column means. This guarantees that you reduces the rank by one (unless it is already singular) before multiplying the matrix with its transpose.” (from MathWorks).

Take this R code:

NROW = 10 NCOL = 10
res <- rep(NA, 1000)
for (i in 1:1000) {
   x <- matrix(runif(NROW*NCOL), ncol = NCOL, nrow = NROW)
   C <- cov(x)
   res[i] <- det(C)
}
hist(res, breaks = 100, main="det(C)")

and look at the histogram:

The reason it's not zero all the time is simply rounding errors. You can see from the spread that it should really be zero.

Anyway, there are many properties of an NxN matrix whose determinant is zero that are all equivalent (here are some) but the one we are interested in is that "the columns of the matrix are dependent vectors in ℝ^N ".

The proof in 2-D looks like this. Take the matrix:

a	c
b	d

The means of the columns are:

μ₁ = (a+b)/2

μ₂ = (c+d)/2

So, centering this matrix results in:

a-μ1	c-μ2
b-μ1	d-μ2

Substituting in the values for μ₁ and μ₂ gives:

(a-b)/2	(c-d)/2
(b-a)/2	(d-c)/2

Multiplying this matrix by its transpose gives:

1/4	(a-b)2+(c-d)2 -(a-b)2(c-d)2 -(a-b)2-(c-d)2 (a-b)2+(c-d)2

and the determinant is therefore:

(1/4²) [ ((a-b)²+(c-d)²)² - ((a-b)²+(c-d)²)² ] = 0

The matrix is always singular meaning necessarily that there is linear dependency among the vectors.

This generalizes to higher dimensions.

Note that we might choose the correlation matrix rather than the covariance matrix (see StackOverflow).

Note also that the “units [of covariance] are the product of the units of X and Y. So the covariance of weight and height might be in units of kilogram-meters, which doesn’t mean much.” [Think Stats p108].

Relationship with Cosine Similarities and Pearson Correlation

The triangle inequality is a defining property of norms and measures of distance. It simply says that for vectors x, y and z that make a triange, x + y <= z. It is a consequence of the law of cosines and a defining property of distances, that is, functions that satisfy:

non-negativity:      f(x, y)   > 0
identity:            f(x,y)    = 0 means x =y
symmetry:            f(x,y)    = f(y,x)
triangle-inequality: f(x,z)    <= f(x,y) + f(y,z)

Only symmetry is true for cosine similarities therefore it is not a true distance  although they can be converted to one (see Wikipedia). Despite this, we can still use it for comparing data.

"Cosine similarity is not invariant to shifts. If x was shifted to x+1, the cosine similarity would change. What is invariant, though, is the Pearson correlation" [1] which is very similar. The relationship between the two looks like this:

Pearson correlation = cosine similarity (x - x̄, y - ȳ)

"People usually talk about cosine similarity in terms of vector angles, but it can be loosely thought of as a correlation, if you think of the vectors as paired samples." [1]

Pearson correlation is the normalized form of the covariance. Or, to put it in Dirac notation:

covariance = <x - x̄, y - ȳ> / n

where n=2 as we're finding the mean which the same as:

cov(X, Y) = E[(X - E[X])(Y - E[Y])]

And so all these concepts are related.

[1] Brendan O'Connor's blog

Setting up a Hadoop/Spark/Hive/Kafka cluster

We wanted full control of what was running in our cluster so we installed our own rather than use Cloudera or HortonWorks. These are some things I learned along the way.

HDFS

Getting a simple distributed filesystem up and running was relatively straight forward. This StackOverflow answer and this gives the minimum work required.

Spark

Spark doesn't need Hadoop to run. So, simplest thing that works, I started a Spark cluster by following the instructions at DataWookie. This made running a Spark shell in a cluster as simple as adding the switch --master spark:SPARK_MASTER_HOST:7077.

Although simple, such configuration brought down my cluster when I tried to do something that was computationally expensive (my Linux boxes died with oom_reaper messages in /var/log/messages).

So, I did the sensible thing and ran it under YARN with the --master yarn switch.

YARN

We didn't find YARN needed much additional configuration. The defaults were sufficient with the exception of:

1. yarn.nodemanager.vmem-check-enabled should be set to false as it kept saying my Spark jobs did not have enough virtual memory (see StackOverflow).
2. yarn.resourcemanager.hostname should be set to the machine where YARN's ResourceManager runs. Failing to do this will likely lead to HDFS running on the other nodes in the cluster but not YARN jobs (see StackOverflow).
3. yarn-site.xml yarn.nodemanager.resource.memory-mb should be set to be a high percentage of the total memory of the cluster (StackOverflow). Note that the default is a measly 8gb so your Spark jobs will silently use few resources if you don't change this.

These are to be set in yarn-site.xml. The only other file to configure was capacity-scheduler.xml as I was the only person using the cluster and wanted to hog the resources (in particular, yarn.scheduler.capacity.maximum-am-resource-percent).

Please do remember to copy the $HADOOP_HOME/etc directory onto all nodes in the cluster when you make changes and restart the cluster for good measure.

Check that all your nodes are running by executing:

yarn node -list

Hive

Hive was the most annoying to get running so beware. It appears that in addition to copying the JARs from this article (Apache) on how to integrate it with Spark, you also must copy the spark-yarn and scala-reflect JARs into Hive's lib directory and also remove the Hive 1.2.1 JARs from the HDFS directory.

You'll need to configure Hive to talk to the database of your choice to allow it to store its metadata. This StackOverflow answer is for MySQL but the principle remains the same for me when I was using Derby.

Derby

Start Derby with something like $DERBY_HOME/bin/startNetworkServer -h HOST where HOST is the IP address you've told Hive to bind to.

This Apache resource is good for checking your Derby server is up and running.

Kafka 

Kafka and Zookeeper were surprisingly easy to set up as the Apache Quickstart page was very well written. It just remains to say that Kafka "keys are used to determine the partition within a log" (StackOverflow) and that if you want to clean the logs this StackOverflow answer may be useful

Check your disk space

Note that the data for old jobs hang around in $SPARK_HOME/work so you regularly need to clean this up.

If you don't, jobs just won't start and Yarn will complain there are not enough resources but not tell you what those deficient resources are. That your disk space is running low is not obvious.

YARN barfs because yarn.nodemanager.disk-health-checker.max-disk-utilization-per-disk-percentage is set to 90%. So, it is not sufficient to have "enough" disk space in absolute terms. See Cloudera for more information.

In addition, you may want to set hadoop.tmp.dir in core-site.xml to a partition with lots of free space or you might see Spark jobs failing with "Check the YARN application logs for more details".

If in doubt...

... check the logs.

Keep checking $HADOOP_HOME/logs, particularly the yarn-*.log files. They'll give an indication of why things are not working although the messages can sometimes be misleading.

The disk space issue above manifested itself when Hive was running even though it was ultimately due to Yarn. Hive would throw timeout exceptions in its log. YARN would say there were not enough resources even though nothing else was running (yarn application -list). Looking at YARN's pending applications, I saw just my job sitting there but not running. The hint was that the equivalent job would run on Spark when running outside of Yarn and on the Spark master.

Hive's logs default to /tmp/HIVE_USER/hive.log and this is the first place to turn when things go wrong.

Chaining Monads

What is exactly going on when we chain monads? Here is some Scalaz code to demonstrate.

First, we create the monad:

import scalaz.Monad

sealed trait MonadX[+A] {
  def run(ctx: Context): A
}

object MonadX {

  def apply[A](f: Context => A): MonadX[A] = new MonadX[A] {
    override def run(ctx: Context): A = f(ctx)
  }

  implicit val monad = new Monad[MonadX] {
    override def bind[A, B](fa: MonadX[A])(f: A ⇒ MonadX[B]): MonadX[B] = 
      MonadX(ctx ⇒ f(fa.run(ctx)).run(ctx))

    override def point[A](a: ⇒ A): MonadX[A] = MonadX(_ ⇒ a)
  }

}

We created two such monads that we will chain:

    case class Context(aString: String, aLong: Long)

    val hello: MonadX[String] = MonadX { ctx: Context =>
      ctx.aString
    }
    val meaningOfLife: MonadX[Long] = MonadX { ctx: Context =>
      ctx.aLong
    }

Not very useful, are they? But you get the idea. Now, all we want is a for-comprehension, so:

  val concatLength: MonadX[Int] = for {
    x <- hello
    y <- meaningOfLife
  } yield (x + y).length

You can think of monads as programs, so let's run it:

  val ctx         = Context("hello", 42)

  val length: Int = concatLength.run(ctx)

Using this highly de-sugared and non-FP code to demonstrate, the flow of control can be given as this:

About to run for-comprehension
==============================
bind: Creating boundHello with fa=Hello, f=<function1>

This is just the first part of our for-comprehension (x <- greeting). Note that nothing further is executed as monads are lazy. All the bind operation did was create a new MonadX containing a function. We never applied that function.

Only when we run the outer monad (concatLength.run(ctx)) does the 'program' execute:

About to run boundHello
=======================
boundHello.run
    boundHello.f(ctx) = 
        Hello.run
            helloFn(ctx) = 
                'hello'
        Hello.run Finished
        <function1>(hello) = 
            bind: Creating boundMeaningOfLife with fa=MeaningOfLife, f=<function1>
            'boundMeaningOfLife'
        boundMeaningOfLife.run
            boundMeaningOfLife.f(ctx) = 
                MeaningOfLife.run
                    meaningOfLifeFn(ctx) = 
                        '42'
                MeaningOfLife.run Finished
                <function1>(42) = 
                    Creating point (7) [Integer]
                    'point'
                point.run
                    point.f(ctx) = 
                        '7'
                point.run Finished
                '7'
        boundMeaningOfLife.run Finished
        '7'
boundHello.run Finished

So, what's happened? Well, first our monads hello and meaningOfLife  have had bind called with them (bind is another word for flatMap in some languages). The reason is that anything in a for-comprehension will have to be flatMapped as that's what we're doing under the covers. Yes,  de-sugared for-comprehension invokes map but a map can be substituted for a flatMap and a point (sometimes called unit, see the monad laws here). And this is where the point comes from in the above flow.

Leveraging this substitution, the Scalaz map, Monad.map is defined as

map[A,B](fa: F[A])(f: A => B): F[B] = bind(fa)(a => point(f(a)))

Since you can't see map and flatMap functions needed by the Scala compiler, where do they come from? They're provided by Scalaz in scalaz.syntax.FunctorOps.map and scalaz.syntax.BindOps.flatMap.

In my heavily de-sugared version of this code, I have given my functions names. But the Scala compiler is also giving me anonymous functions (<function1>) . These appear to be y <- meaningOfLife for the first one and the yield function, (x + y).length, for the second.

So, in a brief, hand-wavey summary: 

1. the outermost monad is not hello but boundHello which wraps it. 
2. boundHello calls run on its hello.
3. It feeds the results from this into its function, f. This happens to be the block of code that is a result of y <- meaningOfLife. Since we're now using the bind/point substitution while mapping, we're given a boundMeaningOfLife.
4. boundHello runs this boundMeaningOfLife, which, being a recursive structure, runs the same steps as 1 and 2 but on its wrapped MeaningOfLife monad.
5. Again, like boundHello in step #3, boundMeaningOfLife calls its f function but this time the result is a point.
6. Again, since it's a recursive structure, run is called upon point which returns the result of the yield function.
7. Then this program's "stack" is popped all the way to the top with our result.

FP and Big Data

Distributed computing relies surprisingly heavily on mathematics.


Associativity

This just says that

A . (B . C) = (A . B) . C

for any operator '.'

Note that addition is obviously associative when subtraction isn't. However, there is a clever hack around this.

"Subtraction is not an associative operation. But it is associative under addition, where subtraction is defined to be the inverse operation, ⁻¹​ :N→N:n↦−n together with addition. This way, I can do subtraction, via addition of inverses, where 3 - 5 is desugared to 3 + 5⁻¹.

"Subtraction as a standalone function is not associative, because 3−(−5−2)≠(3−(−5))−2. But addition is, in fact, principled, associative, commutative, preserves identity, respects inverses, and is semantically what you want.

"We care about associativity ... in the context of distributed programming ... because it means we can arbitrarily partition data. We care about identity because it means we know that recursion will terminate when we reach the identity element (it is a fixed point modulo our operation), and we know that if a function is also commutative, then we can arbitrarily combine our results.

"This is the model that Spark uses underneath the hood - you'll notice it only takes commutative monoidal functions for its jobs (i.e. sets of data with an operation that has an identity within the dataset, and is associative + commutative). The algebra of distributed programming is quite simple."

[Emily Pillmore, Gitter]

"In this context this means that if we have pipelines A, B and C, and we wish to join them to form a single pipeline ABC, we can either join A and B and then join the result AB and C, or we could join A to pipeline BC. " [Stephen Zoio's blog]


Commutivity

"Suppose you wrote this:

someNumbers.reduce(_ - _)

What would you expect the result to be? The short answer is: we don’t know. Since the operation is not associative and is not commutative, we have broken both constraints we need to have this operation work well. " [Nicholos Tietz-Sokolosky's blog]

Domain

"Unless the set A [to which binary operations are applied] is chosen carefully, they [binary operations] may not always be defined. For example, if one restricts one's attention to the positive integers, then the expression 3-5 has no meaning. There are two conventions one could imagine adopting to this. One might decide not to insist that a binary operations should be defined for every pair of elements of A, and to regard it as a desirable extra property of an operation if it is defined everywhere. But the convention actually in force is that binary operations do have to be defined everywhere, so that "minus," though a perfectly good binary operation on the set of all integers, is not a binary operation on the set of all positive integers." [Princeton Companion to Mathematics]

Similarly, a bad choice of domain can blow up in something like Spark. If you're adding Ints, their total will not exceed the maximum value, right?


Composability

"Whenever one is looking for a general solution to composability, monads are normally not too far away. By composability, we mean the output of one process is the input into another. And that is precisely what we are trying to do. A data processing pipeline consists of several operations, each joined together to form the pipeline. And we can join two or more pipelines to form larger pipelines."
[Stephen Zoio]

Zoio gives an example where wrapping calls to Spark in Monads makes things more re-usable but although the article is otherwise excellent, I found this unconvincing as each Spark call was tied to another. Instead, I refactored his Scalaz/Monad code to make it a little more divisible here.


Equality

Rob Norris in his "Gazing at the Stars" lecture makes some interesting points. He describes how a simple type alias of type Angle = Double caused thousands of dollars of damage as it was not sufficient for all comparisons. (The lecture is largely about mapping between different units).

At 11'34" into the video, Norris tells us Scala 2 lets you compare anything for equality, for example a URL and a Classloader. That's a bug that should be fixed in Scala 3. FP programmers prefer type classes as a solution for example Eq (pronounced "eek").

Note that Norris uses CatsSuite for testing. His code for the lecture can be found here, here and here. Norris' use of CatsSuite throws a lot of random numbers at the unit tests to try and flush out errors.


Monoids

The mathematical definition of Monoids and some Scalaz examples can be found in a previous post of mine here.

Note that even Haskell doesn't (can't?) enforce the Monoid rules. "There must be a value that acts as the identity with respect to the binary function and that the binary function must be associative. It's possible to make instances of Monoid that don't follow these rules... Haskell doesn't enforce these laws so we need to be careful that our instances do indeed obey them." [LYAHFGG]

The best way I've seen so far to address this issue is by the types giving a hint. The function "reduceLeft is not parallelizable, as it explicitly assumes a form that is not associative: (B,A) => B. As long as you use an associative operator, reduce is parallelizable." [SO] Whereas, (B,B) => B would imply associativity.

You need a monoid for fold operations to work properly and a commutative semigroup for reduce. "Commutative semigroups consists of all those semigroups in which the binary operation satisfies the commutativity property that ab = ba for all elements a and b in the semigroup." [Wikipedia]. Monoids are a specialization of semigroups that have an identity element.


Optimization

Obeying the FP laws allows Spark etc to optimize the code on our behalf.

"We can do this because we know our programs are invariant under certain substitutions. Referential transparency gives us the most basic rule: we can inline anything, and factor anything out. Lawful abstractions give us more substitution rules for our toolbox. Functor tells us that fa.map(f).map(g) = fa.map(f andThen g) for instance, so we get an optimization opportunity that we can apply freely in our programs." [tpolecat Jan 09 17:33, Gitter]

However, as with all performance testing, it's an empirical science. Using the almost canonical example of an AST for basic arithmetic, Noel Walsh ("Unite Church and State") says: "each constructor in the FP implementation becomes a method in the OO implementation. This removes the need for the pattern matching. The transformation is known as Church encoding.

"We can go the opposite way as well: convert every method in the OO representation into a case class in a sealed trait (an algebraic data type), and then use a pattern match to chose the action to take (a structural recursion). This transformation is known as reification.

"We can think of OO and FP not as entirely different programming paradigms but choices we make to encode a solution to the problem according to the tradeoffs we want to make.  Let’s see an example where where this transformation is useful. One aspect ... is performance. In the FP representation we must allocate memory to hold the data structure that represents the operations we want to perform."

Frameless

The Frameless project brings a more FP approach to using Dataframes. Now, one can already use the more strongly-typed Datasets to give you compile-time error checking but "unfortunately, this syntax does not allow Spark to optimize the code." Frameless apparently gets around that problem will still giving type safety.

However, I did change a field name but foolishly did not change it elsewhere in the codebase where I called partitionBy. I asked on their Gitter channel how Frameless could save me from this but was told:

Cody Allen @ceedubs Jan 12 19:19
@PhillHenry as far as I can tell, Frameless doesn't specifically have support for that partitionBy method. But in general that is the sort of problem that frameless solves.

So, maybe a future pull request...?