Depth First Development in NYC

By: Ilya Kavalerov

[ilyakava]

Apple Magic Mouse on Ubuntu 18

Used this driver, needed to do the troubleshhot, and every reconnect:

~/Downloads/Linux-Magic-Trackpad-2-Driver/linux/drivers/hid sudo rmmod hid_magicmouse sudo insmod ./hid-magicmouse.ko

To change speed to be slower:

xinput --set-prop "My Mouse" "Coordinate Transformation Matrix" 0.6 0 0 0 0.6 0 0 0 2 xinput --set-prop "My Mouse" "libinput Accel Speed" -0.7

Can get name with xinput --list --short

[ilyakava]

Terminal Web Browser with Images macOS

Sometimes it's better to browse in simplicity, but without images many websites become hard to navigate. Here's what web browsing in the terminal with images looks like, which even works on a normally ad filled website like business insider:

First, you'll have to use iTerm 3.3 or later. Install the image converter for the terminal:

brew install libsixel

`

I had already installed w3m earlier, but there was no image support by default. To enable it:

brew edit w3m

Change --disable-image to --enable-image and save. Then:

brew reinstall --build-from-source w3m

Launch with:

W3M_IMG2SIXEL=/usr/local/bin/img2sixel w3m -sixel -o display_image=1 google.com

Scrolling with the mouse lags and it is better to page up/down by holding the fn key with up and down arrows. Press Shift+h for the help menu.

The above works on OSX 10.14.6 Source Helpful iterm2 cursor movement settings iTerm color theme like the Bloomberg terminal

[ilyakava]

Record System Audio macOS Mojave

Install soundflower, then follow the directions on that website:

start playing audio

Open Audio MIDI Setup: (found in /Applications/Utilities)

hit the '+' button in the bottom left corner and select "Create Multi Output Device"

in the the panel that appears on the right, select "Built-in Output" AND "Soundflower (2ch)"

then hit the button with the gear icon bottom left and select "use this device for sound output" (you should still be hearing your audio it is now going to both Soundflower and the built-in speakers)

Then either in Quicktime, or in Audacity, select "Soundflower (2ch)" as the input. In quicktime Do "Edit > New Audio recording" and hit the arrow near the record button. In audacity use the dropdown immediately to the right of the microphone symbol, right of the "Core Audio" dropdown.

[ilyakava]

Some notes on MBR and GPT

I had a mysterious problem with a 10TB HD. In the end I reformatted it after zapping the gpt with gdisk, but perhaps it was salvagable... It had 1 200MB EFI partition and 1 NTFS partition originally. Then plugging into windows gave a "Disk corrupted" popup. ckdisk couldn't run. And right clicking and getting properties on the disk and running the check under the "Tools" tab also wasn't permitted because it was recognized as corrupted.

What started it all:

$ sudo gpt -vv -r show /dev/disk2s2 gpt show: /dev/disk2s2: mediasize=10000619536384; sectorsize=512; blocks=19532460032 gpt show: /dev/disk2s2: MBR at sector 0 gpt show: error: bogus map gpt show: unable to open device '/dev/disk2s2': Undefined error: 0

$ diskutil list /dev/disk2 (external, physical): #: TYPE NAME SIZE IDENTIFIER 0: GUID_partition_scheme *10.0 TB disk2 1: EFI EFI 209.7 MB disk2s1 2: Microsoft Basic Data 10.0 TB disk2s2

$ sudo gpt -r -vvv show /dev/disk2s2 gpt show: /dev/disk2s2: mediasize=10000619536384; sectorsize=512; blocks=19532460032 gpt show: /dev/disk2s2: MBR at sector 0 gpt show: /dev/disk2s2: MBR part: type=114, start=1920221984, size=1816210284 gpt show: /dev/disk2s2: MBR part: type=108, start=1936028192, size=1953653108 gpt show: error: bogus map gpt show: unable to open device '/dev/disk2s2': Undefined error: 0

$ sudo gpt -vv -r show /dev/disk2s1 gpt show: /dev/disk2s1: mediasize=209715200; sectorsize=512; blocks=409600 gpt show: /dev/disk2s1: MBR at sector 0 start size index contents 0 1 MBR 1 409599

$ sudo gpt -vv -r show /dev/disk2 Password: gpt show: /dev/disk2: mediasize=10000831348736; sectorsize=512; blocks=19532873728 gpt show: /dev/disk2: PMBR at sector 0 gpt show: /dev/disk2: Pri GPT at sector 1 gpt show: /dev/disk2: Sec GPT at sector 19532873727 start size index contents 0 1 PMBR 1 1 Pri GPT header 2 32 Pri GPT table 34 6 40 409600 1 GPT part - C12A7328-F81F-11D2-BA4B-00A0C93EC93B 409640 2008 411648 19532460032 2 GPT part - EBD0A0A2-B9E5-4433-87C0-68B6B72699C7 19532871680 2015 19532873695 32 Sec GPT table 19532873727 1 Sec GPT header

$ sudo fdisk /dev/disk2 Password: Disk: /dev/disk2 geometry: 146467/255/63 [2353004544 sectors] Signature: 0xAA55 Starting Ending #: id cyl hd sec - cyl hd sec [ start - size] ------------------------------------------------------------------------ 1: EE 1023 254 63 - 1023 254 63 [ 1 - 4294967294] 2: 00 0 0 0 - 0 0 0 [ 0 - 0] unused 3: 00 0 0 0 - 0 0 0 [ 0 - 0] unused 4: 00 0 0 0 - 0 0 0 [ 0 - 0] unused

According to wikipedia MBR limits the maximum addressable storage space of a partitioned disk to 2**32 * 512 bytes. Note that for this 10TB HD the sector size is 512. 2**32 - 4294967294 = 2, so this MBR is phony.

$ sudo fdisk /dev/disk2s2 Disk: /dev/disk2s2 geometry: 146442/255/63 [2352590848 sectors] Signature: 0xAA55 Starting Ending #: id cyl hd sec - cyl hd sec [ start - size] ------------------------------------------------------------------------ 1: 72 80 13 10 - 371 101 51 [1920221984 - 1816210284] 2: 6C 357 43 4 - 367 32 52 [1936028192 - 1953653108] 3: 00 0 0 0 - 0 0 0 [ 0 - 0] unused 4: 00 0 0 0 - 0 0 0 [ 0 - 0] unused

Both those entries above are 0.85 and 0.9 TB. The disk had 4TB of data on it before things soured.

$ sudo fdisk /dev/disk2s1 Disk: /dev/disk2s1 geometry: 812/8/63 [409600 sectors] Signature: 0xAA55 Starting Ending #: id cyl hd sec - cyl hd sec [ start - size] ------------------------------------------------------------------------ 1: 00 0 0 0 - 0 0 0 [ 0 - 0] unused 2: 00 0 0 0 - 0 0 0 [ 0 - 0] unused 3: 00 0 0 0 - 0 0 0 [ 0 - 0] unused 4: 00 0 0 0 - 0 0 0 [ 0 - 0] unused

$ sudo gdisk /dev/disk2 GPT fdisk (gdisk) version 1.0.5 Partition table scan: MBR: protective BSD: not present APM: not present GPT: present

One of the most suspicious things was:

$ diskutil info /dev/disk2s2 Device Identifier: disk2s2 Device Node: /dev/disk2s2 Whole: No Part of Whole: disk2 Volume Name: Mounted: No Partition Type: Microsoft Basic Data File System Personality: NTFS Type (Bundle): ntfs Name (User Visible): Windows NT File System (NTFS) OS Can Be Installed: No Media Type: Generic Protocol: USB SMART Status: Not Supported Disk / Partition UUID: 8B773B6E-F8E4-44BF-99AA-A228114F3377 Partition Offset: 210763776 Bytes (411648 512-Byte-Device-Blocks) Disk Size: 10.0 TB (10000619536384 Bytes) (exactly 19532460032 512-Byte-Units) Device Block Size: 512 Bytes Volume Total Space: 0 B (0 Bytes) (exactly 0 512-Byte-Units) Volume Free Space: 0 B (0 Bytes) (exactly 0 512-Byte-Units) Read-Only Media: No Read-Only Volume: Not applicable (not mounted) Device Location: External Removable Media: Fixed

0 Bytes for total space and free space.

$ diskutil info /dev/disk2s1 Device Identifier: disk2s1 Device Node: /dev/disk2s1 Whole: No Part of Whole: disk2 Volume Name: EFI Mounted: No Partition Type: EFI File System Personality: MS-DOS FAT32 Type (Bundle): msdos Name (User Visible): MS-DOS (FAT32) OS Can Be Installed: No Media Type: Generic Protocol: USB SMART Status: Not Supported Volume UUID: 0E239BC6-F960-3107-89CF-1C97F78BB46B Disk / Partition UUID: AB192D63-96F2-49A1-932D-3709E7CAFA36 Partition Offset: 20480 Bytes (40 512-Byte-Device-Blocks) Disk Size: 209.7 MB (209715200 Bytes) (exactly 409600 512-Byte-Units) Device Block Size: 512 Bytes Volume Total Space: 0 B (0 Bytes) (exactly 0 512-Byte-Units) Volume Free Space: 0 B (0 Bytes) (exactly 0 512-Byte-Units) Read-Only Media: No Read-Only Volume: Not applicable (not mounted) Device Location: External Removable Media: Fixed

What I tried to do

sudo gpt destroy /dev/disk2 sudo gpt create -f /dev/disk2 sudo gpt add -i 1 -b 40 -s 409600 -t C12A7328-F81F-11D2-BA4B-00A0C93EC93B /dev/disk2 sudo gpt add -i 2 -b 409640 -s XXX -t EBD0A0A2-B9E5-4433-87C0-68B6B72699C7 /dev/disk2

Where XXX was one of 19532460032, 19532464052. I got the idea for the start -b so that there was no gap in the gpt table. And I got the idea to make the gap 3 in the gpt table between the NTFS and the second copy of the GPT with this -s 19532464052. Windows changed its behaviour towards the drive. At this point I should have tried ckdisk on windows again. Perhaps I got the start of the partition right, but not the end. It's worth noting that once, I deleted everything in the GPT table except the EFI partition. This led to windows only recognizing the 200MB EFI and having no alerts/problems. But this was the first time ever I had seen this 200MB partition pop up in windows.

[ilyakava]

Erase all data and all partitions

Removing the EFI parition on a backup drive for windows proved tricky with normal means. The following solved this.

brew install gdrive

diskutil list to find out which hard drive.

sudo gpt -r -vvv show /dev/disk2 to read gpt table on hard drive.

If you want to wipe everything, better than gpt destroy is gdisk's zap. sudo gdisk /dev/disk2, then 3 to use gpt (in my case was already after gpt destroy but everything was fine). Then the options are:

x z Y Y

You can also write 1MB of zeros to the front of the drive with:

sudo dd if=/dev/zero bs=512 count=2048 of=/dev/disk2

And eject. Then plugged into windows and formated.

source

[ilyakava]

My 1st Kaggle ConvNet: Getting to 3rd Percentile in 3 months

The Diabetic Retinopathy challenge on Kaggle has just finished. The goal of the competition was to predict the presence and severity of the disease Diabetic Retinopathy from photographs of eyes. I finished in 20th place using a Convolutional Neural Network (ConvNet). In this post I'll explain my learning process and progress as I implemented my first ConvNet over the last 3 months. Throughout, I'll link to the implementations in my code, which is available on github for anyone who wishes to replicate my score.

My progress over all my 170+ experiments. See the Misc section at the end for a list of the improvements that each point represents here as written along the x-axis. Each point represents an experiment that set a personal record high (on a validation set) of the Kappa score that the competition is judged by. Each point contains a description of the change that caused the improvement (each improvement is accumulated over the experiments, i.e. later runs include all the past improvements).

Introduction

The Problem

Diabetic Retinopathy (DR) is one of the most significant complications of diabetes and is a leading cause of blindness. Early detection and treatment is essential for preventing blindness. Ophthalmologists can use a lens to look through the dilated pupils of a patient and see the retina at the back of the eyeball, looking for symptoms that indicate changes in blood vessels (NIH). At worst this means that new blood vessels are growing (proliferative DR or PDR) and disturbing the retina, otherwise the patient has non-proliferative DR (NPDR). For the challenge, there are 5 stages of DR severity that have symptoms (from here):

non-pathological, no NPDR

Mild NPDR, microaneurysms (red blotches) which are the source of hard exudate (high contrast yellow spots) sometimes in circulate patterns

Moderate NPDR. "More than just microaneurysms," perhaps cotton wool spots (fuzzy light blotches)

Severe NPDR: IRMA (shunt vessels), venous bleeding in 2+ quadrants, 20+ intra-retinal hemorrhages, no signs PDR

Neovascularization (often vessels with loops or very squiggly vessels), vitreous/preretinal hemorrhage, PDR

Examples of pathological classes 1,2,3,4 from left to right

The goal of the competition is to build a classifier that takes in these images, and outputs an integer diagnosis 0-4.

My Approach: Convnets

A ConvNet correctly predicts class 2 This plot shows the original 128x128 training image (41188_right) with four heatmaps. Each heatmap corresponds to one of four pathological classes. Each pixel in the heatmap represents the probability that that class is true given that pixel is corrupted in the original image. The heatmaps were created by moving around a 11x11 pixel block which hid a region of the image at testing time. The low probability areas (blue) mean that without those pixels that class would not be predicted. Here the network has learned to find the subtle hard exudate and (less so) cotton wool spots indicate classes 1 and 2 (since those classes would not be likely if those areas were obfuscated).

A ConvNet predicts 0 for a true class 2 This time we see, in the probability heatmap for class 1, the cluster of hard exudates at the bottom middle and single exudate top right are picked up as important features. It is also interesting to see that when the bottom letterbox and the middle lens blur defects are covered up the probability of the pathological class rises. However, the probability of each of the pathological cases was just below the threshold for prediction (0.5) when no part of the image was covered up, so we have a false negative here.

Since my last Machine Learning class, I've been looking forward to using ConvNets because of the promise of end-to-end learning: learning a feature extractor and a classifier simultaneously. This property allows accurate classifiers to be created without much domain knowledge. So in May, I read through a Stanford tutorial and Columbia reading list and toured theano while implementing the best performing ConvNet on MNIST. After that, I moved onto the more challenging DR dataset, with guidance from Sander Dieleman's posts on his two Kaggle ConvNet wins: classifying galaxies and plankton.

Software and Hardware

I used same software setup as Sander Dieleman in his Galaxy post. For details see my repo. For hardware, I had access to a Nvidia M2090 with 6Gb RAM, and at the last moment a K40. Because of these assumptions, I use the lasagne.layers.cuda_convnet module which requires a GPU to run.

Of the networks I include, 128x128 runs (kappa ~0.68) took me 3.5 mins/epoch, 192x192 runs (kappa ~0.72) took 7.8 mins/epoch, and 256x256 runs (kappa ~0.74) took 14 mins/epoch.

Adding resolution provided huge gains As seen in these validation errors for runs 120 through 162 (run descriptions in Misc section at end). A whole new cliff of lower errors is discovered simply by using larger images!

Preprocessing

Resizing with Graphicsmagick

The supplied data consisted of JPEGs which were often 16 megapixels. For every size that I experimented with (128,152,192,236,256,292) my downsampling was the same strategy: crop out the surrounding black, size down the width, then vertically crop/letterbox until the image is square. I used PNGs because they are my favorite lossless format. I did not play around with different ways to downsample with graphicsmagick, though I am curious if it would have made a difference.

I played around with the idea of removing black pixels by using a log-polar transformation on the images, but was unhappy with the distortion effect since the images are not all aligned in the same way (macula is not always centered).

Input = Raw Pixels

Other than subtracting the mean image, and dividing by the standard deviation image, I made no modifications to the images entering the network. Early runs used grayscale images to lessen runtime.

Normalization had no gains

For input into the network, histogram normalization with graphicsmagick gave marginal improvements using grayscale images, and slightly worse performance when using color images. I never tried local contrast normalization simply because it has fallen out of favor.

Exploiting Invariances by adding Noise

Adding random flipping and color casting noise was essential to delay overfitting and continue decaying the error (as seen in the validation error of runs 40 through 120 plotted here (run descriptions in Misc section at end)). Notice how the earlier runs are bowl shaped while later runs are a cliff: adding color deepens the bowl, a custom error function widens it, but adding noise to the dataset turns it into a cliff.

Flipping

I randomly flipped/didn't flip an image along the horizontal axis, and then again for the vertical axis every time it entered the network. I liked that this transformation didn't change the image quality or the appearance of the letterbox which was generally symmetric. I didn't try rotating the eyeball a random degree because of the potential extra delay when loading the image into the network (and also because I feared the rotated letterbox might add harmful noise to the dataside), but given more time I would have tried this.

Color Casting

I first heard about this from a paper out of Baidu. I randomly decided on each channel whether or not to add/subtract a constant, and then drew this constant from a centered gaussian distribution of a standard deviation 10. This means that 99% of the values I added were between [-30,30]. I experimented with other ranges, ±30 worked best for me.

Reducing Noise

Alignment

Early on, I had the belief that any reduction in noise in the training set would improve my convergence. For this reason, I decided to align all the training set images. Some images display a tab jutting out of the eyeball in the right half of the image. When this is the case, this means the image is inverted. I wrote a pattern recognition style tab detector that with 90% accuracy would detect this tab, and with the additional information of right/left eye from the image name, would output the flip that would put the optic nerve in the right of the image.

A sample of 115 images after alignment, most of the time the optic nerve is on the right.

This pre-alignment gave a 1% improvement over the same run with the images oriented as they are found in the training set. However random flipping provided a 10% benefit, so this work did not prove useful.

Noise decay

After some experiments I started to believe that for ConvNets, adding the right kind of noise prevents overfitting by roughening the error surface of the network, and lowering the energy barriers to hop out of bad local minima. I thought that by having high noise in the beginning of training, and less noise at the end, I would have the advantage of getting out of bad local minima early on, and being able to stay in good minima later on in training.

So I tried decaying the amount of noise: randomly flipping pre-aligned images less often as training progressed. I was surprised to find that shortly after every time I reduced the noise, my network began to overfit. I rationalized this as me reducing the amount of data (total unique images) that the net had access to over time.

Network Architecture

The basis for my network is the popular VGGNet from Oxford that Andrej Karpathy recommends.

I tried several variations, all of which are in my network_specs.json, but the most successful one looked like:

9 Weight Layers Output Shape (batch size 128) input 256x256x3 (3, 256, 256, 128) conv3-32 (32, 254, 254, 128) maxpool size 3 stride 2 (32, 127, 127, 128) conv3-64 (64, 125, 125, 128) maxpool size 3 stride 2 (64, 62, 62, 128) conv3-128 (128, 60, 60, 128) conv3-128 (128, 58, 58, 128) maxpool size 3 stride 2 (128, 29, 29, 128) conv3-128 (128, 27, 27, 128) maxpool size 3 stride 2 (128, 13, 13, 128) conv3-128 (256, 11, 11, 128) maxpool size 3 stride 2 (256, 6, 6, 128) FC-2048 (128, 2048) maxpool size 2 (128, 1024) FC-2048 (128, 2048) maxpool size 2 (128, 1024) FC-4 (128, 4) sigmoid

This network had 21.66 million parameters and used up almost all of my available 12Gb of GPU RAM. This same network with an input image of size 192x192 would fit in 6Gb of GPU RAM. Each convolutional layer (except the first) has dropout with p=0.1, and every convolutional layer has an LReLu non-linearity, and each FC layer has dropout with p=0.5.

The other successful network in network_specs.json is similar, but would take smaller input images, and has 1 fewer convolution layer.

The initial network decisions were made with the help of the Stanford tutorial. Given more time, I would have experimented with larger filters on larger images (3x3 filters were best on 128x128 images), as well as larger overlaps in the pooling to preserve spatial information.

Parameter Sharing Attempts

I tried the same solution Sander Dieleman did in his plankton challenge: splitting the image into quarters, running each quarter through the convolutional layers independently, and then connecting the features from the four quarters to the same series of 3 FC layers. I called this Fold4xBatchesLayer in the code because the pixels are being folded across the batch dimension in the 4D input data tensor. This marginally worsened results. Given more time, I would have combined this folding pixels across batches with the pre-alignment strategy.

I also tried a similar strategy but with folding pixels across channels, which led to a considerable runtime speedup but 10% worse performance.

I would have also liked to try sharing parameters between pairs of images (left right). One problem is that not all image pairs share the same diagnosis.

Error Function and Number of Output Nodes

I started by using the same error function I did for MNIST: categorical cross entropy after a softmax (all class probabilities sum to 1) non-linearity. This has the downside of not encoding any information that the classes are ordinal (4>3>2>1>0 in the severity of DR) and not differentiating between errors of different magnitudes (unlike the metric that the competition is judged on).

For this reason I followed the advice on the Kaggle forums and used an nn-rank target matrix with relative entropy as described in this paper after a sigmoid layer (each class probability is between [0,1]). This resulted in the following target matrix for a four node output network:

[ [0, 0, 0, 0], [1, 0, 0, 0], [1, 1, 0, 0], [1, 1, 1, 0], [1, 1, 1, 1] ]

instead of the standard "one-hot" target matrix (which requires an additional output node):

[ [1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 1] ]

I also tried an nn-rank target matrix with a 5 node output, but did not get superior results.

Variants

To emphasize gross errors, I modified the target matrix to be more punitive on errors of larger magnitude:

[ [ 0, 0, 0, 0 ], [ 1, 0, 0, 0 ], [ 1.3, 0.6, 0, 0 ], [ 1.5, 1, 0.5, 0 ], [ 1.6, 1.2, 0.8, 0.4] ]

When I tried matrices that were more steep (concentrated most of the weight on more discrepant label-prediction pairs), I noticed that the network made fewer underestimates, and more overestimates. I saw the opposite effect as I tried more flat matrices (somewhere between the last one shown and the binary nn-rank target matrix). The best performing target matrix I tried is shown here and comes from the code here.

If I had more time, my next step would have been to experiment directly with exaggerating the underestimate penalty and under-emphasizing the overestimate penalty to exert more control over the false-negative <-> false-positive tradeoff.

Prediction

The actual prediction was pretty trivial. Since each output node after the sigmoid non-linearity take a value between 0 and 1, I count the number of nodes greater than. There were very few discontinuities here (98% of examples always had a lower probability for a more pathological case). But I also feel like maybe I could have done something more elaborate here to boost performance.

Training/Optimization

Batch Selection

Each minibatch had the same proportions of labels as the entire training set. This led to a 2 to 3% performance improvement, and reduced the noise in the training error.

Validation

I laid away 15% of the training set for validation, and held this set constant for my last 100 runs. Luckily, my validation set score was always ±0.5% of my Kaggle Public leaderboard submissions.

SGD

Training was done with SGD and Nesterov momentum (always 0.9). The majority of runs were with minibatches of size 128, so there were ~240 gradient steps per epoch.

Init

Initialization proved important. Normal initialization prevented training altogether, and GlorotUniform worked best.

Learning Rate Decay

Decaying the learning rate after the validation error stalled for 3 consecutive epochs led to minor improvements of 1 to 2%.

Regularization

Very few experiments were run with L1 and L2 regularization. L1 always did considerable damage, and L2 didn't seem to help. With more time, I would have explore this.

Dropout

Dropout was essential for the network to even train for long. This plot shows the validation error on runs 14 through 40 in the summary (run descriptions in Misc section at end). Before GlorotUniform, the error was a totally flat line at 1. Adding dropout and pooling prevented the overfitting seen in GlorotUniform, and overlapping pooling led to a quicker error drop.

Ensemble

I got about 1 to 2% improvements when I averaged the raw outputs of several networks (last 4 nodes after the nonlinearity for each example) and then made my prediction. I wonder if there would have been a difference if I averaged before the non-linearity.

Misc

List of improvements (Left to Right on Summary graph)

Numbers correspond to the experiment number in the first summary graph that each cumulative change corresponds to.

run 13: vgg_mini7

run 14: GlorotUniform Init

run 16: All Conv dropout

run 18: Both FC pooling

run 19: 1 more FC dropout

run 21: Overlap Pooling

run 40: LReLu

run 55: controlled batch distributions

run 80: nnrank-re

run 83: color

run 91: 4 outputs

run 93: random flips

run 120: ±20 color cast

run 122: ±30 color cast

run 136: kappa weighted error func

run 149: 152px

run 151: 192px + extra Pool

run 152: 192px + extra ConvPool

run 162: 256px

Runs 14 through 40

Runs 40 through 120

Runs 120 through 162

Dealing with Class Imbalance

Since this is a real world medical dataset, it can be expected that pathological cases are greatly outnumbered by healthy ones. In this case, the class proportions are: [0.73, 0.07, 0.15, 0.02, 0.02] for classes [0,1,2,3,4]. I tried changing the class proportions in each minibatch to be more uniform. As I leveled out the populations, performance worsened as a result of the majority of the 0 labels being classified as pathological. This underperformance on class 0 was likely a result of undersampling that class, while overfitting from oversampling the rarer pathological cases set in.

I also tried training the network to discriminate the hardest classes first (only training on 0,1 examples, then only training on 0,1,2 examples, and so on), but this led to no performance increase either.

Batch Size

I didn't experiment much with the batch size (stuck with 128 throughout), but would have liked to effectively increase it over the duration of the experiment (by accumulating gradients across minibatches before updating) to have less noise in the gradient steps.

Thinner Conv Layers

Thinning out the CONV layers led to 2 to 3% worse performance, but could reduce memory consumption and runtime by 1.5-2x. Compensating with additional input image resolution would have been nice.

Conclusion

I had a lot of fun, and am thankful to Kaggle and the Sponsors for making such an exciting and challenging dataset available, as well as the open source authors behind theano, Lasagne, and pylearn2. Next steps for me will be experimenting with even larger images (I had trouble above 256px with my current setup), writing my own Cuda code for fun and potential performance gains, learning how to parallelize across multiple GPUs, and maybe experimenting with other frameworks.

Links

HN

r/MachineLearning

Kaggle

DataTau

[ilyakava]

Facebook React Meetup Summary

Last Wednesday I went to a meetup at Facebook, the creator of React (a frontend JS framework), to hear about React’s inception from Martin B on the photo sharing team.

Martin’s team developed React, while also internally testing immutable-js before it was open-sourced.

The main motivation behind React is to ensure that your code grows in complexity as slowly as possible. Making new feature work easier for existing developers, and especially for new developers to on-board quickly. Martin focused on the benefit of transitioning from imperative code (based on modeling transitions between application states), to declarative code (based on modeling instances of state).

Typically, when writing an application with global synthetic events, its tempting to think of the states that the application can take and write transitions between them that are triggered by some event. Facebook’s example is a bell icon to represent user notifications that can be either blank, have a badge with a number, or have a “on fire” decoration in addition to a badge with a large number. After the page is loaded from the server and the user interacts with it, you can imagine that the page needs to have code that can transition between these three states, that is, there must be 6 transitions hard-coded (A->B, A->C, B->C, B->A, C->A, C->B).

Think about how much more code we have to write when we model a transition… Since one transition corresponds to each unique combination of two states, the number of transitions as a function of n states is a n-choose-2 problem, and thus grows order O(n(n-1)). In other words, transitions grow quadratically with respect to your n states. So why not forgo the imperative style of modeling transitions, and instead model states declaratively?

To do this, we move rendering to the client, and instead of having a number of event based functions to have access to transition the state of a element on the DOM, we update a single source of state in memory, which would lead to the parts of the DOM that rely on that portion of mutated state to self re-render.

To do this while keeping the web-page performant, the React needs to reconcile changes in state very quickly. Since the DOM is slow, we can’t keep the state in the page’s DOM. And since tracking object changes between two deeply nested and complicated objects is a hard problem, we need a fast identity based comparison between instances of data.

Both these things are achieved by immutable data-structures that share data and self-identify by reference. As shown in the picture below, when running something like: ys = xs.insertNode(8) on a tree xs, we return a reference to a brand new tree ys that shares its mutual structure with xs, saving memory with shared redundancy. Comparing the two trees is easy, since they have totally different references (1 in blue and 1 in orange).

Introducing immutable data-structures thus led to reduced memory consumption. It also led to huge perf(ormance) gains (in terms of frame-rate) which follow from functional lazy data-structures that avoid creating intermediate data-representations during sequential execution.

[ilyakava]

Saving Bootcamp after repartitioning

Recently I decided to upgrade my secondary MacbookPro5,4 (6 years old!) from double-boot (Windows 7 and OSX 10.8.5) triple boot, adding the latest Debian Jesse.

Somewhere along the way, I ended up doing this to my windows Bootcamp partition:

The Cause

Specifically how I got here

I went from 2 partitions (OSX + Bootcamp), to adding extra NTFS and HFS partitions (without altering the Bootcamp partition) with Mac’s Disk Utility, to adding reFind, and then partitioning during my Debian install to add a swap-space partition from the space I had created for Debian.

More generally, what’s wrong?

A macbook’s hard-drive has a couple of ways that it keeps information about its partitions. Primarily there is GPT which OSX uses, and MBR which Windows needs to boot from (read more here).

On my other primary macbook, (single-boot and factory 1 partition hard-drive), this is what GPT looks like:

$ sudo gpt -vv -r show /dev/disk0 gpt show: /dev/disk0: mediasize=500277790720; sectorsize=512; blocks=977105060 gpt show: /dev/disk0: PMBR at sector 0 gpt show: /dev/disk0: Pri GPT at sector 1 gpt show: /dev/disk0: Sec GPT at sector 977105059 start size index contents 0 1 PMBR 1 1 Pri GPT header 2 32 Pri GPT table 34 6 40 409600 1 GPT part - C12A7328-F81F-11D2-BA4B-00A0C93EC93B 409640 975425848 2 GPT part - 53746F72-6167-11AA-AA11-00306543ECAC 975835488 1269536 3 GPT part - 426F6F74-0000-11AA-AA11-00306543ECAC 977105024 3 977105027 32 Sec GPT table 977105059 1 Sec GPT header

And this is what the unused MBR looks like:

$ sudo fdisk /dev/disk0 Disk: /dev/disk0 geometry: 60821/255/63 [977105060 sectors] Signature: 0xAA55 Starting Ending #: id cyl hd sec - cyl hd sec [ start - size] ------------------------------------------------------------------------ 1: EE 1023 254 63 - 1023 254 63 [ 1 - 977105059] <Unknown ID> 2: 00 0 0 0 - 0 0 0 [ 0 - 0] unused 3: 00 0 0 0 - 0 0 0 [ 0 - 0] unused 4: 00 0 0 0 - 0 0 0 [ 0 - 0] unused

So, as you can guess already from the neglected state of the MBR, what most likely killed Windows in my case, is using Disk utility to move around the partitions. Since it likely updated my GPT which the OSX side needs, but not the MBR which the Windows side needs.

Fixing the Incongruence

The solution is to sync up the important bits of the changed GPT with the MBR. For me this is my old Windows partition (which we’ll see was moved to a different place on my hard-drive, even though a different partition was shrunken by Disk Utility). In my case this is my Windows partition, but this could just as easily be a linux partition (example).

The Specific Problem

On the now triple boot macbook, I take a look at the MBR, which I copied into an Excel table to sort by partition size. I did this since all the computer names were indecipherable (contents column), but I remembered the relative sizes of the partitions I had created and labeled them (my name column).

After running sudo gpt -vv -r show /dev/disk1 and labeling the output:

start size index contents my name in MBR slot 409640 1340525080 2 GPT part - 48465300-0000-11AA-AA11-00306543ECAC Mac OSX 2 1341196864 370117188 3 GPT part - C12A7328-F81F-11D2-BA4B-00A0C93EC93B Debian 3 1742800896 126998528 4 GPT part - EBD0A0A2-B9E5-4433-87C0-68B6B72699C7 Extra NTFS space 1869813760 83709952 5 GPT part - EBD0A0A2-B9E5-4433-87C0-68B6B72699C7 Bootcamp 1711314944 31485952 6 GPT part - 0657FD6D-A4AB-43C4-84E5-0933C84B4F4F Debian Swap space 40 409600 1 GPT part - C12A7328-F81F-11D2-BA4B-00A0C93EC93B 1 1340934720 262144 1869799424 14336 1953523712 1423 1711314052 892 2 32 Pri GPT table 1953525135 32 Sec GPT table 34 6 0 1 MBR 1 1 Pri GPT header 1953525167 1 Sec GPT header

Compare this to my MBR table:

21:57:40:☯ sudo fdisk /dev/disk0 dyld: DYLD_ environment variables being ignored because main executable (/usr/bin/sudo) is setuid or setgid Disk: /dev/disk0 geometry: 121601/255/63 [1953525168 sectors] Signature: 0xAA55 Starting Ending #: id cyl hd sec - cyl hd sec [ start - size] ------------------------------------------------------------------------ 1: EE 0 0 1 - 1023 254 63 [ 1 - 409639] <Unknown ID> 2: AF 1023 254 63 - 1023 254 63 [ 409640 - 1340525080] HFS+ 3: 83 1023 254 63 - 1023 254 63 [1341196864 - 370117188] Linux files* *4: 07 1023 254 63 - 1023 254 63 [1711314052 - 83711408] HPFS/QNX/AUX

Slot 4 in my MBR table clearly has a bad start position and size compared to the above GPT table. This is what we need to fix. So, with reference to the above GPT table’s values (start/offset 1869813760 and size 83709952) I did:

21:57:47:☯ sudo fdisk -e /dev/disk0 dyld: DYLD_ environment variables being ignored because main executable (/usr/bin/sudo) is setuid or setgid Password: fdisk: could not open MBR file /usr/standalone/i386/boot0: No such file or directory Enter 'help' for information fdisk: 1> edit 4 Starting Ending #: id cyl hd sec - cyl hd sec [ start - size] ------------------------------------------------------------------------ *4: 07 1023 254 63 - 1023 254 63 [1869813760 - 83709952] HPFS/QNX/AUX Partition id ('0' to disable) [0 - FF]: [7] (? for help) 07 Do you wish to edit in CHS mode? [n] n Partition offset [0 - 1953525168]: [1711314052] 1869813760 Partition size [1 - 83711408]: [83711408] 83709952 fdisk:*1> flag 4 Partition 4 marked active. fdisk:*1> write

For editing a partition slot for linux you would follow an analogous procedure, but use a Partition id of 83.

And lo and behold, I restarted, held down option during boot, and successfully entered my unaltered Windows 7 installation!

References

Difference between GPT and MBR Same problem after installing Yosemite Similar problem

[ilyakava]

3 makefile musts

I recently had a burst in makefiles required in some of my projects, and struggled with a few make no brainers:

1) setting a variable dynamically

Although it looks harder than it should be:

$(eval CURR_USER="$(shell whoami;echo $)")

will set CURR_USER to be a string in double quotes to the result of your normal whoami call in bash.

This can come in handy for iterating through files:

$(eval EPS=$(shell ls ./plots/*.eps;echo $)) allEpstopdf: for f in $(EPS); do epstopdf $$f; done

2) if checks

The tricky thing here is that the condition of the ifeq has to be indented the same level as the task definition, and the body of the conditions has to be indented 2 tabs in:

download_geolite: # normally, everything is indented 1 tab $(info 'Checking for GeoLiteCity.dat file...') ifeq (GeoLiteCity.dat, $(shell ls ./3-enrich/scala-kinesis-enrich/ | grep GeoLiteCity.dat)) $(info 'Found GeoLiteCity.dat file...') else $(info 'No GeoLiteCity.dat file found. Downloading...') wget http://geolite.maxmind.com/download/geoip/database/GeoLiteCity.dat.gz -O ./3-enrich/scala-kinesis-enrich/GeoLiteCity.dat.gz gunzip ./3-enrich/scala-kinesis-enrich/GeoLiteCity.dat.gz endif

ifneq is also available. I had trouble using variables within the ifeq‘s parentheses, perhaps some sort of interpolation issue.

3) you can change directories

This one was news to me, but to have your commands truly chained, they should be on the “same line,” separated by semicolons and backslashes if they fall to multiple lines:

cd ./deploy; git init

will work differently than:

cd .deploy git init

Bonus.

$(info 'Checking for GeoLiteCity.dat file...')

is much nicer than

echo 'Checking for GeoLiteCity.dat file...'

since it only prints the message once.

[ilyakava]

N Languages N Months Meetup - Haskell

Last night I went to a meetup where Eric Jones introduced Haskell while live coding some circuits, code here.

The theme of Eric’s presentation was the idea that the choice of a programming language for a project is a choice of the culture that will come out of it, and that Haskell is one of many languages that breeds a fast-paced development culture.

What we need in a language

First, Eric introduced the abstraction that when people speak with each other, we can reduce their conversation to a recurring cycle with two nodes: Communication <-> Verification. One person talks, the other tries to understand, and talks about what they understood or did not understand. He stressed the point of communication as the action that starts this infinite feedback loop.

Eric then rephrased this social interaction loop for software developers, which instead starts at a more passive point than initiated “communication.” For an individual software engineer on a team, the social existence loop is more like: Observe -> Orient -> Decide -> Act -> …and cycle back… This is General John Boyd’s, a military strategist’s, OODA loop.

Eric then used this loop to phrase the central task of a developer: to always be oriented well in the codebase. And as a result, code should always be written to best communicate to the developer reading it. The code that developers write should be a domain specific language (DSL). It should “[build] up a common, rigorous language between developers and users,” what Eric Evans called “Ubiquitous Language” (Martin Fowler). At the same time, we should be aware that the readability of a language depends on what the particular reader brings to the text, and we should perhaps augment our writing to account for an already existent culture we may want to integrate into. Assumptions about syntax (ex: the ' (prime) symbol indicating that x and x' are related yet different) become a form of redundancy and familiar comfort that welcomes certain readers.

Why Haskell

Eric believes that Haskell is a great language for a group to work in because it is well suited to build a DSL, and thus it is a tool (one of many other suitable languages) to control the intellectual complexity of the programs that come out of a domain. Eric then live coded some circuits (code here), a half adder and a full adder, and while doing so he emphasized some advantages of Haskell:

It is easy to create types and to build your DSL on the spot. This gives a ‘pure’ feeling to the code we write, we can work directly in our high domain terms.

Logic through pattern matching builds referential transparency. It was easy to combine to half adders into a full adder and have it work as expected.

Conclusions

After this detailed live coding session, Eric returned to the question he opened with: Do the languages that we speak shape the way we think? to paraphrase the Wall Street Journal. Ultimately, the Meetup’s consensus was that they do, and for programming languages specifically, the things that a language makes easy and makes hard will determine what people do with it, and thus the culture that emerges out of it.

Further reading

Structure and Interpretation of Computer Programs, online text

Eric Jones at StrangeLoop

QuickCheck for testing the correctness of your programs with Property based testing, and data generation

Haskell in production environments

[ilyakava]

Introduction The Indianapolis Museum of Art (IMA) recently shared thousands of high-resolution images from its permanent collection with Artsy, …

Discuss on hn

Download as a Pdf here

[ilyakava]

Demystifying PCA

I've read some magical things about PCA, such as: it identifies the most important dimensions in your data, and therefore tells you which kinds of measurements are redundant. I've seen many references discuss it like a black box, always saying something appealing like: find out which parts of your data account for the most variance, and thereby find a better way to re-express your dataset. PCA sounded so alluring to me, so I couldn't help but be frustrated by my inability to interpret the output of an inscrutable Rlang method that did all the work for me.

I'm collecting my newfound insight here because it was difficult to find the corners of the Internet that documented PCA without assuming the reader could derive it themselves. In my search for enlightenment (which I primarily found in Jon Shlens' great tutorial), I realized that PCA isn't nearly as magical as I thought it was.

What you need to enjoy PCA

PCA is applied statistics and linear algebra, so an understanding of variance and covariance from Ross' book helped me greatly. Before learning about eigenvectors and eigenvalues from Bretscher's incredible book, PCA didn't really mean anything to me. Learning the ins and outs of SVD also helped, since the two are similar techniques, but SVD is easier to understand in my opinion.

The Motivation: A perfect covariance matrix

The juice of PCA is in the Linear Algebra, but its motivation is a piece of statistics:

A typical problem with data is figuring out what measurements are "useful." Often we are interested in some kind of change in data, and knowing which of our measurements are related to that change. An example: we were video-taping a falling ball from street level a few blocks away from the tower of Pisa with a sort-of-sideways stationary camera. When we watch that tape now, we record the x,y coordinates of the ball on the screen, as it moves diagonally on our television. In this case we know that there is only 1 dimension important for a ball flying "straight down," but our shoddy camera work has led to two intertwingled variables.

We want to cut out the redundancy in our dataset: we wish that our measurements were unrelated to each other, that they did not vary together in any organized way, so that we had 1 and only 1 measurement per independent factor in the process we were recording. A bad dataset would have non-zero covariances between variables, a perfect dataset would have zero covariance between any two measurements. Covariance between two random variables X and Y (of which each single measurement x or y is an instance of), is defined like so:

\(Cov(X,Y) =E[(X-E[X])(Y-E[Y])]\)

In the case of datasets, the measurements of a single variable X (say the x axis on the video tape in our example) across all trials (individual frames of the video tape) would be held in a vector. The expectation E[X] of that variable would simply be the mean of that vector. The covariance between the X and Y vectors would simply be their dot product (after their mean was subtracted). If we had many variables, we would create a square (and symmetric) covariance matrix where the *ij*th entry would be the covariance between the *i*th and *j*th variable. The diagonal entries would be the covariance between a given variable and itself - which is equal to that same variable's variance.

So in the perfect case, where no measurement is intertwingled with another's, we have a covariance matrix that has zeroes everywhere, but along the diagonal (each measurement is independent of every other, but has some variance). In the bad case, our covariance matrix is not such a diagonal matrix, since each variable is dependent on another in some way, and has non-zero covariances.

It's Basic Linear Algebra

Once we've stated our goal: that we wished our data had a diagonal covariance matrix, we only need to rearrange a few equations too see how we need to adjust our dataset.

Our initial state is some data X with the covariance matrix: \(S_X = \frac{1}{1-n}XX^{T}=cXX^{T}\)

But we imagine a better version of the data Y with the (diagonal) covariance matrix: \(S_Y = \frac{1}{1-n}YY^{T}=cYY^{T}\) (eq. 1)

So lets just assume that there exists some matrix P that can linearly transform our data X into the perfect data Y: \(Y=PX\) (eq. 2)

All we have to do is substitute eq. 2 into eq. 1: \(S_Y = c(PX)(PX)^{T}=cP(XX^{T})P^{T}=PS_XP^{T}\) (eq. 3)

Lets introduce a new variable, \(P=E^{T}\), that will rephrase eq. 3 (into something very familiar) in terms of \(S_X =cXX^{T}\). Since \(S_Y = PS_XP^{T}\) we can rearrange into: \(S_X = P^{T}S_YP\) and then substitute in our new variable to get: \(S_X = ES_YE^{T}\) (eq. 4).

This is very informative! Eq. 4 is the typical change of basis formula from linear algebra (when E is orthonormal as it is for necessarily symmetric covariance matrices). E diagonalizes \(S_X\), and is an eigenbasis (see Caveat section) of \(S_X\). We also see that \(S_Y\), the perfect covariance matrix, has the eigenvalues of \(S_X\) on its diagonal (the variances of Y). I said "an" eigenbasis because we can always rearrange the rows and columns of E and \(S_Y\) for the equation to still be valid (for PCA we want to arrange them so that the eigenvalues in \(S_Y\)'s diagonal are descending).

If we wanted to get to "perfect data" Y from X, since we know E, and \(P=E^{T}\), we just multiply the transpose of the eigenbasis by our data: \(Y=E^{T}X\). Since E is orthonormal, \(P=E^{T}=E^{-1}\), this transformation is a rotation. The rows of \(E^{T}\) are called loading vectors, and they define the directions of greatest variance. The rows of Y, which are the data X projected onto the loading vectors, are called the Principle Components of X (or Score Vectors).

Big Picture

Wow. PCA is just a change of coordinates problem, a kind of problem found in the beginning chapters of Linear Algebra texts. We have our data, with one kind of measurement dependent on another, and it feels messy. All we are really doing is looking for a new coordinate system, which is a linear combination of the old one, where each axis is independent of all others (like in the only coordinate system that we are used to: the Cartesian coordinate system).

Caveat

I've been referring to E as an eigenbasis, but most likely it will not be a real basis since the covariance matrix will probably not have full rank. If the geometric multiplicity of the covariance matrix (number of repeat eigenvalues) is greater than zero, then you technically don't have enough eigenvectors to fill all the dimensions of a space, and thus you don't have a basis, but just a collection of eigenvectors.

It's a piece of SVD

SVD is a matrix decomposition technique that I talk about exclusively in another post of mine. SVD factorizes a matrix X into sums of pieces of X (U and V are orthonormal bases of \(XX^{T}\) and \(X^{T}X\) respectively).

\[X=UΣV^{T}=u_1σ_1v_1^{T}+u_2σ_2v_2^{T}+\cdot\cdot\cdot+u_rσ_rv_r^{T}\]

We can actually get the result we want for PCA from SVD. SVD finds a way to break up a matrix A into 3 matrices: \(X = UΣV^{T}\). Notice that if we multiply each side by its transpose:

\[ XX^{T}=(UΣV^{T})(UΣV^{T})^{T}=(UΣV^{T})V(UΣ)^{T}=UΣ(V^{T}V)Σ^{T}U^{T}=UDU^{T} \]

This just tells us that we can perform SVD on our data matrix, the U matrix will be an eigenbasis of the covariance of our data.

This is useful to know because it can be dangerous to multiply a matrix by its transpose - in some cases you can lose a lot of information. Computing the SVD of a matrix can forgo this step, and is preferred, say in the way R computes SVD when using its prcomp function.

What do we do with E?

From \(E^{T}\)

\(E^{T}\)'s rows are the principle component loading vectors, and the eigenvectors of X's covariance matrix.

Relate features to one another

An example can be found in James' Book (p328) where, in a State Crime dataset, features Murder, Assault, and Rape occur prominently in the 1st loading vector, but not the 2nd. This prompts an understanding as these 3 correlated features belonging to a single "Serious Crime" theme.

From Y

Y contains the data X projected onto the principle component loading vectors, and its rows are the principle components.

Generalized insights

In the above crime example, we could generalize that states with a large 1st principle component value have high crime rates: reading 1 dimension summarizes 3 features (the 3 kinds of crimes measured).

The same procedure could generate a single scalar performance metrics for businesses to focus on, which can summarize multiple related goals

Find groups

Still in the same above crime example, while the theme of the 1st principle component was "serious crime," the 2nd was "urban population." Plotting states (trials) in a plot with the 1st vs 2nd principle component could separate groups of states and classify them somehow - some fitting into the high crime category - and another in high population low crime category.

Be forewarned however, that PCA will not do any "separation work" on the data. It will only capture the major pattern in fewer dimensions.

As a preprocessor

As discussed in the linear algebra section, the eigenvalues along the diagonal of \(S_Y\) correspond to the variances in \(S_X\) captured in each eigenvector in \(E^{T}\). We can use a metric for each eigenvector's "proportion of variance explained" that is its eigenvalue divided by the sum of eigenvalues. We might find that, cumulatively, the "proportion of variance explained" among the first few eigenvectors is 80% or 90%. This would allow us to throw away the last few rows of Y in confidence that those principle components explain little variance in the data and are just noise. After reducing the size of our dataset in Y, we could then perform other techniques on this now purified and distilled data.

principle components regression

clustering

When will this not work, and how can we tell?

In PCA, we are optimizing a transformation of our data that gives us the directions of the largest variance: these might not always correspond to the axes in the underlying basis of the system. This is because we are optimizing a covariance matrix. We make the assumption that the information that this matrix contains, the mean and variance of our random variables, is enough to describe our random variables. This is actually only the case for variables with a Gaussian distribution. Thankfully, most natural data is Gaussian distributed.

John Shlens shows how PCA fails with exponentially distributed data, and introduces another technique, ICA, for dealing with non-Gaussian data.

Sources and Additional reading

John Shlen's tutorial on PCA is great for a more detailed look.

The ISL book is great for describing PCA with R, and has labs/example problems, but assumes an understanding of PCA's derivation.

Johnson's book gives a good look at PCA's description/derivation and has some good interpretations on performing PCA.

This hw assignment has a good couple of exercises that helped my understanding.

For a quick review of the background concepts, this seems like a good source.

[ilyakava]

Thinking with Streams - NY Open Stat Meetup

Last night I went to a meetup where Mike Dewar, a data scientist with the NY Times, talked about a tool his team created to aid exploratory data analysis of streaming data. An introduction of the visual programming tool, streamtools (written in go and JavaScript), was followed by a demo. The content of the talk can also be found on the streamtools homepage.

Motivation behind Streamtools

Dewar introduced his technical team at the Times as a forward thinking team. Streamtools was the answer they created to the question: in 3-5 years, what kind of questions do we need to be able to answer effectively and quickly?

His team speculated that in the future:

1) All data comes as a stream (a continuous infinite source) 2) Online (rather than batch) data analysis techniques will be the standard 3) Responsive tools (like if this then that) will be the norm for making sense of things (more via Bret Victor)

Some consequences of this are: - data analysis will be closer to the sensors: instead of querying a structured database, measurements straight off of a sensor will be examined - the responsibility of data storage will be abstracted from the data source

In addition, Dewar’s team hoped that more people in the future would be data-curious. They wanted to create a tool that could be used by writing code for it, but also have a gui component geared towards excel users.

Demo

The tool itself is reminiscent of MaxMSP (a visual programming tool for artists). The idea is that you build a sort of circuit out of discrete modules called blocks. Each block may have an input, a transformation of that input, and an output.

The demo included a board where a module received POSTs from mobile phones reporting their orientation angle, and then sent this to a sparkline plot of the average phone angle that updated several times a second. This was used as a voting mechanism pitting userRs and pythonistas against each other (the pythonistas had a slight lead in this crowd).

Screenshots of the rotate-to-vote page on a phone

Another demo recreated a Poisson process by using a ticker module (periodic event emitter), to sample from an exponential distribution module, to feed back into the ticker module setting a new event emitting rate, to record the event in a log module.

Use cases

Streamtools is billed primarily as a thinking tool, but is also used in permanent production roles at the NYTimes

NYTimes employees (Newsroom analysts and a/b testers) make stat boards (like geckoboards) with it

Non-technical users who have a clearly defined problem, and a well formatted data stream to start analyzing, have no problem picking it up

Go lang supports ARM architecture so using streamtools with a Raspberry Pi or a Beaglebone (with its i/o pins) is easy

Streams are the main target data, but csv parsers also exist

Streamtools can back other data visualizations - like this globe showing pageviews around the world of the NYTimes site

[ilyakava]

Creating my first epub - pandoc rules!

Today I found an excuse to make my first epub. Hearing some buzz about the next great startup-related book, I ran into the blog that originated it, which I had only partially read a long time ago. There is a pdf version available of it, but I prefer epubs, so here is how I made one.

I first jQueried out the urls of each blogpost into a file hrefs.txt (19 urls for 19 chapters). Then I piped each to pandoc (brew install pandoc) to generate a markdown file for each:

for NUM in $(seq 1 19) do pandoc -s -r html $(sed -n "${NUM}p" hrefs.txt) -o "ch${NUM}.blake.masters.markdown" done

Above, we iterate through the chapter numbers, and each is used with sed to reference a line in the file with the hrefs. I did some monkeying around and deleted the comments sections and some twitter headers (manually :(). Though perhaps the files were regular enough that I could have regexed for some lines indicative of content.

The second and last step with pandoc was just to pass the markdown files as arguments to pandoc: pandoc -o cs183_notes.epub $(ls | grep markdown)

So the lesson was that it's super easy to make an epub of your own from web content! For an example of what the final result looks like, compare the first blogpost to the whole epub.

[ilyakava]

Play w/ Scala and Precompiled Assets: a match made in hell

The story goes that I have a play (2.2) app that uses the css preprocessor stylus. There is a play plugin that compiles the css in the target directory (excluded from git). Everything works really well in development, so I decided that I should replicate my local setup in production, on heroku.

To do this, I would need to install both the node and scala buildpacks on the same instance, using the multi buildpack to install both. Although the installation of both went fine, I was stuck with a java.io.IOException: Cannot run program "stylus" Exception. Who knows what path stylus is available in? Going off of this SO question I was inspired to export additional directories into the PATH on heroku. Through this, I discovered that tampering with the BUILD_DIR env var on the heroku cli (via something like heroku config:add BUILD_DIR=/app) would cause a preposterous sbt.version error recorded here (deleting the app and recreating was quick, so that’s how I reverted). I realized that instead, this additional PATH export had to be made in the buildpack. So I forked the buildpack, exported all sorts of new paths where I thought my stylus executable might be installed, like:

$HOME/node_modules/.bin $HOME/node_modules/stylus/bin /app/node_modules/stylus/bin /app/node_modules/.bin $BUILD_DIR/node_modules/.bin $BUILD_DIR/node_modules/stylus/bin

But all to no avail, the same java.io.IOException: Cannot run program "stylus" exception continued.

At this point, I spoke to another developer, and understood that compiling the styl to css on the server is pretty wasteful anyway, and ditched this idea. Instead, I’m making a task on Jenkins CI that will compile and upload these assets to S3, and then have a config setting in my Play app to change the asset location lookup between production and development.

Would be nice if serving assets in play was as easy and mindless as in rails :(.

[ilyakava]

OSS maven hosting with scala/sbt

Today I found an excuse for learning how to get a scala library on a public maven repository.

Thankfully, the instructions for publishing a scala library to an open source repo are available online (though it was difficult to surface), and pretty easy to follow (much easier than the sbt instructions). My only emendation would be on registering the project with sonatype. I was confused about the naming of the groupId for the project, for which I found the answer, for libraries hosted on github, here. The answer being to use com.github.username, and that’s it (at first I submitted com.github.username.projectname but it was amended by my reviewer). I was also amazed at the speed of review, it only took an hour after I submitted! The groupId also needs to be set to organization in build.sbt, something I overlooked until I got errors.

On getting an identity for the pgp signing (which seems similar to the workings of ssh, and also seems to be popular for sending ‘secret’ email messages): I generated the keypair from the gpg cmd line tool (brew install gpg). A key with a pass-phrase does need to be generated (hitting enter to skip the pass-phrase is okay for the cmd line tool, but hitting enter while publishing aborts the process).

Right after uploading, I was able to resolve the package on sonatype at this address.

[ilyakava]

Gerhard Richter, SVD, and me

Inspired by Gilbert Strang's mention of using a Linear Algebra technique to compress images, I tried it out for fun, and discovered that this form of lossy compression has an interesting form of distortion resembling Gerhard Richter's squeegee blur.

I would have liked to choose some cheerier photos that Gerhard Richter painted off of and then blurred, but examples of these were hard to come by. These images are from the aftermath of the Baader-Meinhof group and are still available today since they were printed in german newspapers (source). Although these images are depressing, Richter's paintings of them (The October 18 1977 cycle) are characteristic of his artistic interests: memory, trauma, and the memory of trauma.

Reproducibility

I've packaged my code into a gem, and the code is on github.

gem install image_svd image_svd --help

A prerequisite is to have imagemagick, and a ruby version 1.9.2 or greater installed.

image_svd -i ~/Downloads/svd_photos/in_RichterfuneralM.jpg -n 10 will create a new image in the current directory using 10 singular values.

Be warned! Compression is by no means quick! A 200x300 image takes 15 seconds on my computer. Finding the eigenvalues of a matrix is no easy computational task (O(n3) in practice, O(n2.3) by super fancy approximation).

The gem also allows archiving an image as the set of singular vectors that make it up (the option -a true). After reading the section below, you might think that this would be the most space saving method to save images (this would be true if the precision of each number in the singular vectors would be the same as that of the original matrix - hint, this is not the case).

The Implementation

Intro

An image, at least a grayscale image, might as well be a matrix where each pixel's brightness is represented by a number. As Gilbert Strang mentions (in Linear Algebra and its Applications on p332), if we need to send the information contained in the image/matrix, it would be better to only send the "essential" information. Luckily, we can break down any matrix into a product of three matricies, which can be broken down into a sum of the columns and rows of these matricies:

\[A=UΣV^{T}=u_1σ_1v_1^{T}+u_2σ_2v_2^{T}+\cdot\cdot\cdot+u_rσ_rv_r^{T}\]

Let's say A is a 1000x1000 image/matrix. A contains 1,0002 = 1,000,000 numbers, while each chunk of the sum contains (2 * 1,000) + 1 = 2,001 numbers. We will see that the most important parts of the image are contained within the first tens of chunks. If we saved the numbers in the first 50 chunks of the SVD, we would save 1/10th the numbers for the same image and decrease disk usage 10 fold!*

SVD by hand

First, lets see what we are going to do on an image by hand so that our strategy is set before we write any code.

Imagine that we have the following 2x3 image/matrix (blown up here, but really is 2x3 pixels):

This is the matrix:

\(A= \begin{bmatrix} 0 & 1 & 1 \\\ 1 & 1 & 0 \end{bmatrix}\)

A quick digression before we start - what is our motivation?

Motivation behind SVD

We are often interested in expressing a matrix A with another diagonal matrix (diagonalizing it) because it makes A easier to work with (easier because it makes an abstract vector space look like Rn). Usually, this means finding an eigenbasis for A (a matrix whose columns are the axes along which A is stable as a transformation). However an eigenbasis has the constraint that its matrix be square, an inconvenient limitation for image processing, and even then, it may not be orthogonal. Instead of finding a single matrix that diagonalizes A:

\(S^{-1}AS=D\) or \(A = SDS^{-1}\)

We will find two: \(A = UΣV^{T}\)

We have two orthogonal bases here. And we see that one is the eigenbasis of ATA and the other of AAT:

\[ A^{T}A=(UΣV^{T})^{T}(UΣV^{T})=(ΣV^{T})^{T}U^{T}(UΣV^{T})=VΣ^{T}(U^{T}U)ΣV^{T}=VDV^{-1} \\\ AA^{T}=(UΣV^{T})(UΣV^{T})^{T}=(UΣV^{T})V(UΣ)^{T}=UΣ(V^{T}V)Σ^{T}U^{T}=UDU^{-1} \]

Together they contain enough information on how to diagonalize A. We could say more, but that is enough. A good proof is in another of Strang's books.

Back to doing it by hand

Right, so what are we going to code up anyway?

First, we find the eigenbasis of ATA:

\(A^{T}A= \begin{bmatrix}0 & 1 \\\ 1 & 1 \\\ 1 & 0 \end{bmatrix} \begin{bmatrix}0 & 1 & 1 \\\ 1 & 1 & 0 \end{bmatrix} =\begin{bmatrix}1 & 1 & 0 \\\ 1 & 2 & 1 \\\ 0 & 1 & 1 \end{bmatrix}\)

We do this by finding its characteristic roots, i.e. eigenvalues:

\[ \begin{align} f_{A^{T}A}(λ)&= det\begin{bmatrix}1-λ & 1 & 0 \\\ 1 & 2-λ & 1 \\\ 0 & 1 & 1-λ \end{bmatrix} \\\ f_{A^{T}A}(λ)&=(1-λ)((2-λ)(1-λ)-1)-1(1-λ) \\\ &=(2+λ^{2}-3λ-1)(1-λ) \\\ &=(λ^{2}-3λ)(1-λ) \\\ &=λ(λ-3)(1-λ) \end{align} \\\ λ_1,λ_2,λ_3=3,1,0 \]

The square roots of these will be our singular values, and make up the diagonal of Σ. Let's find their eigenvectors that span its eigenspaces:

\[ \begin{align} E_3&= ker\begin{bmatrix}-2 & 1 & 0 \\\ 1 & -1 & 1 \\\ 0 & 1 & -2 \end{bmatrix}=\begin{bmatrix}1 \\\ 2 \\\ 1 \end{bmatrix} \\\ E_1&= ker\begin{bmatrix}0 & 1 & 0 \\\ 1 & 1 & 1 \\\ 0 & 1 & 0 \end{bmatrix}=\begin{bmatrix}1 \\\ 0 \\\ -1 \end{bmatrix} \\\ E_0&= ker\begin{bmatrix}1 & 1 & 0 \\\ 1 & 2 & 1 \\\ 0 & 1 & 1 \end{bmatrix}=\begin{bmatrix}1 \\\ -1 \\\ 1 \end{bmatrix} \end{align} \]

Let's skip straight to the good part, the sum of the chunks of the SVD (a sum of rank 1 matrices).

We have the eigenvectors sorted by their eigenvalues, and we normalize them (make their lengths 1) to form v1 and v2.

For us, u1 and u2 can be found by transforming v1 and v2 by A, and dividing by the singular values:

\( v_1= \frac{1}{\sqrt{6}}\begin{bmatrix}1 \\\ 2 \\\ 1 \end{bmatrix} ,v_2= \frac{1}{\sqrt{2}}\begin{bmatrix}1 \\\ 0 \\\ -1 \end{bmatrix} \)

\[ \begin{align} u_1&=\frac{1}{\sqrt{λ_1}}Av_1 \\\ &= \frac{1}{\sqrt{3}} \frac{1}{\sqrt{6}} \begin{bmatrix}0 & 1 & 1 \\\ 1 & 1 & 0 \end{bmatrix} \begin{bmatrix}1 \\\ 2 \\\ 1 \end{bmatrix} \\\ &= \frac{1}{\sqrt{2}} \begin{bmatrix}1 \\\ 1 \end{bmatrix} \end{align} \]

Similarly:

\( u_2= \frac{1}{\sqrt{2}} \begin{bmatrix}-1 \\\ 1 \end{bmatrix} \)

And boom, their combination yields the original matrix:

\[ \begin{align} A&=UΣV^{T}=u_1σ_1v_1^{T}+u_2σ_2v_2^{T}=u_1{\sqrt{λ_1}}v_1^{T}+u_2{\sqrt{λ_2}}v_2^{T} \\\ &= \frac{1}{\sqrt{2}}\begin{bmatrix}1 \\\ 1 \end{bmatrix} {\sqrt{3}} \frac{1}{\sqrt{6}}\begin{bmatrix}1 \\\ 2 \\\ 1 \end{bmatrix} + \frac{1}{\sqrt{2}}\begin{bmatrix}-1 \\\ 1 \end{bmatrix} \frac{1}{\sqrt{2}}\begin{bmatrix}1 \\\ 0 \\\ -1 \end{bmatrix} \\\ &= \frac{1}{2}\begin{bmatrix}1&2&1 \\\ 1&2&1\end{bmatrix} + \frac{1}{2}\begin{bmatrix}-1&0&1 \\\ 1&0&-1\end{bmatrix} \\\ &= \begin{bmatrix}0&1&1 \\\ 1&1&0\end{bmatrix} \end{align} \]

Or in images:

(caveat - negative valued pixels don't exist in my file format, so I took some artistic liberty with representing those colors)

Implementation

Very cool! In the section above we verified that a potentially big matrix can be decomposed into sums of vector-matrix products (where each element in the sum is a rank 1 matrix). Now comes the step of translating the work by hand into so code.

I chose to write the code to automate this in ruby. The reason being that the tediousness of finding the characteristic roots of a matrix is bundled into the language, but the manipulation of these elements into A's SVD still has to be done. I don't want to write 10,000 lines of code, but writing a hundred or so lines of code is more fun than writing one (like with R: svd(matrix)).

The excerpt repeated here captures everything done by hand above:

def decompose(m_A) m_AT = m_A.transpose @m, @n = m_A.to_a.length, m_A.to_a.first.length m_ATA = m_AT * m_A dcmp = Matrix::EigenvalueDecomposition.new(m_ATA) evs = dcmp.eigenvalues # eigenvectors are already normalized and in same order as eigenvalues sorted_eigenvectors = dcmp.eigenvectors.each_with_index .sort_by { |_v, i| -evs[i] } .map { |v, _idx| v } both = (0...@num_singular_values).map do |idx| u = sorted_eigenvectors[idx] sigma_vT = (m_A * u).covector [u, sigma_vT] end @sigma_vTs = both.map { |p| p.last } @us = both.map { |p| p.first } end # rubocop:enable MethodLength def reconstruct_matrix zero_matrix = Matrix[*Array.new(@n) { Array.new(@m) { 0 } }] (0...@num_singular_values).reduce(zero_matrix) do |acc, idx| acc + (@us[idx] * @sigma_vTs[idx]) end.transpose end