Randomness....! Is random the same as random?

Often I need a random number (as a seed for a split, or to randomly order records). Now you might think that,  if we have a randomnumber generator, you could make things MORE random if you multiply these (random * random is randomer). Let's see if this holds true by just trying! 

Say we have 20,000 random numbers between 0 and 1 we can plot the number against their frequency using a histogram (100 bins):

Now this looks pretty random, 100 bins means that the space is divided in 100 sections. So 1 bin (section) with all values between 0.00 - 0.01 , 1 bin with all values between 0.01-0.02. 

We observe that most bins have a value between 0.009 & 0.011 and all between 0.008 & 0.012. which means that there are about 0.010 * 20,000 = +/- 200 points generated with a value between 0 and 0.01. 

There are also 200 point generated with a value between 0.01 and 0.02 and so on. This indeed looks rather random.  We can also plot the index (1^st number generated to 20,000^th number generated) versus the value:

Again this looks al pretty random, and so it should. The whole space is filled evenly.

Now here comes the kicker: If we multiply these random numbers (the same in the same order) with another random number things start to look less random:

Again we can plot the index versus the number:

What we see is that the numbers are actually getting less random. So a given value occurs more frequently. Now let's look at the bin plot again:

Remember just now, there were approximately 200 points with a value between 0.00 and 0.01. Now we see 0.06 * 20,000 = 1200, conversely we see only 0.00005 * 20,000 = 1 datapoint with a value between 0.99 and 1.00! This is looking much less random. 

And this effects gets worse when we use 3 times random numbers multiplied:

0.165 * 20,000 = 3310 values between 0.00 - 0.01  and 0 values between 0.99 - 1.00. The highest bin is between 0.89 and 0.90 with 0.0001 * 20,000 = 2 data points. 

And 4 times random numbers multiplied…

0.33030 * 20,000 = 6606 values between 0.00 and 0.01 & again 0 in the highest bin. There is 0.00005 * 20,000 = 1 datapoint with a value between 0.86 and 0.87. (perhaps not very visible in the second plot)

So there you have it, while you might think that random * random is more random than random, it is NOT! This is stated by the central limit theorem. In probability theory, the central limit theorem (CLT) states that, given certain conditions, the arithmetic mean of a sufficiently large number of iterates of independent random variables, each with a well-defined expected value and well-defined variance (all will be between 0 and 1), will be approximately normally distributed, regardless of the underlying distribution. (see [1] http://www.math.uah.edu/stat/sample/CLT.html, [2] Rice, John (1995), Mathematical Statistics and Data Analysis (Second ed.), Duxbury Press, ISBN 0-534-20934-3), also [3] wikipedia, https://en.wikipedia.org/wiki/Central_limit_theorem )

Part V: ZFS - Poor man's deduplication - Howto: Ubuntu Home Server (ZFS + virtual IpFire)

Series : Ubuntu Home Server (ZFS + virtual IpFire) 

Part I   : Ubuntu install and setting up IpFire 

Part II  : Setting up ZFS

Part III : Setting up hardware monitoring and webmin

Part IV: How much space do you lose with a RaidZ1/RaidZ2/RaidZ3?
-->Part V : Poor man's deduplication  (this post)<--

Poor man's deduplication..

Table of Contents:

Why deduplication?

How to check for duplicate files

Processing and sorting your list
Consolidating directories
Removing duplicates

Why deduplication?

I use our ZFS server / NAS mostly for backup and network storage. As I am working in cheminformatics my job produces a lot of data, supervising students add to this.. Much of this data and old files need to be stored (in particular if you have published on this data, it needs to be reproducible). I like to think I am pretty organised in the way I store data, but turns out I am not.... I work from home sometimes, do work om my workstation at work, sometimes on a laptop (used to make backups on external drives in the pre-zfs era). 

When I installed my NAS I consciously turned off deduplication for two reasons, firstly the server hardware was likely not powerful enough (X3 with 16GB RAM). Secondly because I thought I did not need it (convinced that my file system was well organised...)

However ZFS showed that I had a number of duplicate blocks.. In fact the output from zdb -S showed that I could in theory obtain a 1.1 reduction (or 9%) by deduplication. The output is listed here:

root@KarelDoorman:~# zdb -S zfspool
Simulated DDT histogram:

bucket              allocated                       referenced          
______   ______________________________   ______________________________
refcnt   blocks   LSIZE   PSIZE   DSIZE   blocks   LSIZE   PSIZE   DSIZE
------   ------   -----   -----   -----   ------   -----   -----   -----
     1    35.7M   4.43T   4.13T   4.15T    35.7M   4.43T   4.13T   4.15T
     2    2.70M    333G    274G    278G    5.80M    717G    589G    597G
     4     275K   30.7G   23.6G   24.3G    1.29M    146G    113G    116G
     8    32.4K   3.39G   2.66G   2.74G     317K   33.0G   26.2G   27.0G
    16    4.02K    254M    207M    224M    81.7K   4.86G   3.94G   4.28G
    32      913   16.9M   9.36M   15.1M    40.1K    685M    365M    625M
    64       59   1.49M    590K    975K    4.83K    124M   46.0M   77.5M
   128       18    574K     15K    144K    3.28K    117M   2.72M   26.2M
   256        6    390K   4.50K   48.0K    2.01K    151M   1.58M   16.0M
   512        5    258K      4K   40.0K    4.10K    151M   3.12M   32.8M
    1K        4    257K      3K   32.0K    5.52K    379M   4.24M   44.1M
    2K        2      1K      1K   16.0K    5.18K   2.59M   2.59M   41.4M
    8K        1    128K      1K   7.99K    10.1K   1.27G   10.1M   81.0M
   16K        1    128K      1K   7.99K    25.8K   3.22G   25.8M    206M
 Total    38.7M   4.79T   4.42T   4.44T    43.3M   5.32T   4.85T   4.87T

dedup = 1.10, compress = 1.10, copies = 1.01, dedup * compress / copies = 1.20

However, this also shows that there are in total 38.7 million unique blocks in my file system, which would require 38.7 million * 320 bytes ~ 11.5 GB of RAM for just the dedup table (requiring in total 4 times that as ARC should be limited to 25 %, so 48 GB).

Now at the current prices (Tweaker.net) this costs about 800 EUR! Corresponding to 5.7 disks of 2TB. So it seems to me that deduplication in ZFS is out of my budget. Hence, poor man's (manual) deduplication.

If you're interested, this is what it looks like graphically (note the logarithmic y scale):

This shows that there are 35,700,000 blocks that are unique, 2,700,000 that exist in duplicate, 27,5000 that exist in triplicate, etc. It even shows that there is a single block present in 16,000 fold. So hence I thought I'd try to do this on the file level, maybe there was something to gain. 

How to check for duplicate files? 

Checking for duplicates can be rather tedious. Files can be named differently, have different timestamps etc. For this a brilliant program has been written : Fdupes. Github page : https://github.com/adrianlopezroche/fdupes. It even has it's own wikipedia page: http://en.wikipedia.org/wiki/Fdupes

"The program first compares file size and MD5 signatures and then performs a byte-by-byte check for verification." 

So you can be rather sure that a duplicate is an actual duplicate. I ran this on my /zfspool folder (which is the root folder of all zfs datasets. In total about 1.7 million files are stored on my zfs pool (taking up approximately 5 TB of space). The results were as follows (note that I processed and grouped the output using pipelining tools (e.g. KNIME or Pipeline Pilot):

So it turns out that of the 1.7 million files about 665692 were duplicates taking up 893.51 Gigabytes!

so much for an organised file system as the majority of these duplicates were actually in my work folder....

Processing and sorting your list

By default fdupes outputs the duplicates with their full location, this can be ported to a text file and then the files are sorted with blank lines separating the groups. 

Now I wanted to easily gain a lot os space and hence I wanted to start with the largest files, moreover I wanted a uniqie identifier per files to browse through. So using two simple bash scripts I did the following:

1 removed the empty lines 

#!/bin/sh

files="/media/dupes/duplicates.txt"

for i in $files

do

  sed '/^$/d' $i > duplicates_out.txt

done

2 for each file calculate the MD5 hash 9also because I am paranoid) and add the size in bytes (for later sorting)

#!/bin/sh

while read name

do

hash=`md5sum "$name" | awk '{print $1}'`

size=`ls -all "$name" | awk '{print $5}'`

echo -e  "$name\t$hash\t$size"

done

Afterward you get a file that lists 3 columns seperated by tabs: name, hash, size (bytes).

Using Pipeline pilot I created a unique ID (hash_size), calculated the size in MB / GB, and calculated first occurences for each ID. (protocol is available HERE, but this can easily be done with KNIME). 

Consolidating directories

From this I observed several directories to be duplicate from other directories, with Rsync you can easily merge them (I used timestamps and kept a log).

rsync -avhP /directory1/ /directory2/ > /loggind_dir/transferlogs/1.txt

Removing duplicates

After this I repeated fdupes and the processing and deleted the remaining duplicates. In total I freed about 700 GB!

Component to calculate Matthews Correlation, Sensitivity, Specificity, PPV and NPV with PP

(I am seeing about 5-10 views a day on the Pipeline Pilot pages, please be so kind to acknowledge / cite my blog when you use these tools and guides)

Why would we want such a thing?
Like with the regression validation parameters, I found that PP lacked a component to calculate correlation coefficients between two properties in the data stream in classification. 

Therefore I have written a component to do just that. One of the features I find useful is the option to include a bar chart that displays the values of the calculated properties on a scale between 0 and 1. This allows a quick visual inspection of your model reliability. When applied to the 'KNN classification of Estrogen Antagonists' from the example protocols, it looks like this:

In addition it outputs the parameters in a shaded table: 

The component calculates these parameters between two properties. Therefore, when using in external validation of a model, the modeled values have to be pre-calculated. 

The component requires that you give the names of the properties carying the measured value, the modeled value and the classes that were modeled. Currently it can only be used in a two class classification. In addition you can choose to also output the original unmodified records through the fail port, while the correlation plot and table are output through the pass port. I have not written a 'how to use' in the help, but will do so next week. In the meantime, the component can be found 

here.

Component for calculating the R^2 zero with PP

(I am seeing about 5-10 views a day on the Pipeline Pilot pages, please be so kind to acknowledge / cite my blog when you use these tools and guides)

Why would we want such a thing?
During the time I have been using PP, I found it inconvenient that there was no component to calculate the correlation coëfficient between two properties present in the data stream (for instance when performing external validation of a model).

Therefore I have written a component to do just that. One of the features I find useful is the option to include both an upper and lower error margin line. Allowing a quick visual inspection of your model reliability.

While in the latest version (8.5) there is a component called "Regression Model Evaluation Viewer" which calculates an RMSE and R², this component has some downsides.

1. The component calculates the modeled values internally, so it cannot be used to calculate the correlation between two sets of values obtained from external sources.
2. The component only calculates the R² and RMSE, while for a proper evaluation R₀² and k-slope are also required.

My component is on my website and compatible with PP 8.5 and up, it can be found 

here.

It has been tested up to a maximum of approx. 20,000 records and works fine. In addition the parameters that are also calculated in the 'Regression Model Evaluation Viewer' and 'R-statistics fit plots' are identical. 


So what does it do?
The component calculates correlation parameters according to Tropsha (2010) ¹ between two properties present in the stream. These properties are defined as 'Activity' (Y-values) and 'Model' (X-values). These have to be present in the stream and therefore need to be pre-calculated in the case of a model. In addition, a scatter plot containing all values is output. Both the parameters and the plot are output as reporting items.

The following values are calculated:

1. RMS Error (RMSE)
2. R² (R2)
3. R₀² (R2_zero)
4. R₀²' (R2_zero_acc)
5. k-Slope (Slope_K)
6. k-Slope ' (Slope_K_acc)
7. % Difference between R² and R₀² (Perc_Diff_R2_with_R2_zero)
8. % Difference between R² and R₀²' (Perc_Diff_R2_with_R2_zero_acc)
9. Absolute difference between R₀² and R₀²' ( Absolute_diff_R2_zero_and_R2_zero_acc) 

Additional Settings:

• Under 'Plot Parameters' variables for the x-y scatter plot can be defined. Furthermore the range of the upper and lower error lines can be set (default 0.5 from the line of unity).
  • 'Auto_range'; when set to 'True' the scale of the axis is automatically defined to the scale of the data. Alternatively; when set to 'False' (default), a range can be entered manually for 'Activity' (y-value) and 'Model' (x-value)(Default is 2.0 - 12.0).
  • 'Uncertainty' defines the margin between the line of unity and the uncertainty lines (default 0.5 units away from line of unity).
  • If 'Uncertainty_in_plot' is set to 'True' (default) then two lines indicating a lower and upper error line are drawn in the plot.
• If 'Output_Records' is set to 'True' all values are output unchanged to the 'Fail' port while the plot and correlation parameters are output to the 'Pass' port.

The examples are made in the example protocol "08 Calculate logP using the R_logP_SVM Model" , listed under Examples/R Statistics/Learning and Clustering/R Learn Models...

RMSE	R2_zero	R2	R2_zero_acc	Slope_K acc	Slope_K	Perc_Diff_R2 with_R2_zero	Perc_Diff_R2 with_R2_zero_acc	Absolute_diff_R2_zero and_R2_zero_acc
0.679	0.839	0.839	0.827	0.997	0.928	0.000	0.015	0.012

If for some reason you are heaving trouble with the component, please contact me!

1. Tropsha, A. (2010). Predictive Quantitative Structure-Activity Relationships Modeling. Handbook of Chemoinformatics Algorithms. J. Faulon and A. Bender.

R-Statistics Error messages in Pipeline Pilot

Updated!
(I am seeing about 5-10 views a day on the Pipeline Pilot pages, please be so kind to acknowledge / cite my blog when you use these tools and guides).

Over the last years I have been using R to create my models. However the interface running on top of R (doing the data shaping and fingerprint folding) was pipeline pilot. This works quite nice and efficient (although one could think of better solutions, but for my works this set up suffices). 

When there are errors in your data though, things go wrong. Not all error messages are as intuitive as you would like. The pipeline pilot help can't really help here either, so over the last years I have kept a list of error codes and what they mean in practice. I have listed it here so that anyone else struggling with an unknown error might find it. however this is also convenient for myself as online these things are retrieved quicker than on network share xxx :). 

The organisation is as follows, the closed dot with italic characters is the actual error message received (trimmed), the white dot with regular text contains a possible cause, the closed square a solution.

Related to SVM as performed in the “e1071” package:

• Error in svm.default(x, y, scale = scale, ..., na.action = na.action) : 
• dependent variable has to be of factor or integer type for classification mode. 
 Calls: doCV -> modelfunc -> svm -> svm.formula -> svm.default
• Fingerprint properties are not recognized as fingerprints
• Set property type of properties to learn from to “fingerprint” (like 'SciTegic.value.IntegerFingerprintValue')
• Set option convert fingerprints to “Fixed-Length array of bits”
• Possibly due to merge there are array properties present (multiple values for one property)

• Error in …. Subscript out of bounds
• The property to learn is incorrect
• Two values  present in one property where there should be one
• Only actives are present
• No properties present to learn from
• Possibly through ignore properties)

• Empty beginning of file
• The property to learn is incorrect.
• Either not present in the stream
• The name is misspelled

• Missing properties in file
• Problem with the fingerprints that are being input in a learned model.
• The ‘change fingerprints to fixed length bit size’ is executed wrongly,
• This specific property is missing
• Set property type to fingerprint has not been performed ('SciTegic.value.IntegerFingerprintValue')

• "Error in svd(x, nu = 0) : 0 extent dimensions"
• When performing a PCA, (multiple) properties are not considered to be numeric.
• Decimal comma instead of dot

• “Error in svm.default(x,y,scale,…..): C <= 0!”

• The allocation of a cost value is incorrect.
• Decimal comma instead of dot

•  “Error in matrix(ret$dec, nrow = nrow(newdata), byrow = TRUE, dimnames = list(rowns,  :   matrix: invalid 'ncol' value (< 0)Execution halted”
• Properties to learn from defined incorrect
• “allpropertiesonfirstdata” instead of “user set” when not all properties are present in all records

• Error in svm.default(x, y, scale = scale, ..., na.action = na.action) :           Need numeric dependent variable for regression. In addition: Warning message:data length exceeds size of matrix
• Property to learn from contains non-numeric characters
• Continuous model selected for classification data

• Error in cor(preds[[1]], preds[[2]], method = "pearson") : missing observations in cov/cor. In addition: Warning messages: 1-5: data length exceeds size of matrix
• Non numeric properties are used to learn from in regression.
• Use ‘IgnoreProperties’ to exclude non numeric properties
• Possibly, property should be changed to ('SciTegic.value.IntegerFingerprintValue') while using regression.

• Error in c(1e-05/nx, 0.001/nx, 1/nx, ) : argument 4 is empty
• Gamma values to be sampled ends with comm rather than value
• Remove comma at the end or add value

• Error in svm.default(x, y, scale = scale, ..., na.action = na.action) :  NA/NaN/Inf in foreign function call (arg 4) Calls: doCV ... modelfunc -> svm -> svm.formula -> svm.default -> .C
• Property to learn from non-numeric
• ‘Inf’ rather than numeric

• Error in withCallingHandlers(expr, warning = function(w) invokeRestart("muffleWarning")) : invalid multibyte string at '<b2>II' (or at '<a0>hydra) Calls: readxy -> cleandata -> FactorOrNumber Execution halted
• Array property present formatted as blabla[1], blabla[2], etc.
• Flatten to single properties (eg turn binary flag on for proeprties present named by value in array property.

Related to decision tree forests as performed in the “randomForest” package:

• Error in randomForest.default(xy[-1], y, ntree = 500, mtry = mtry, importance = imp) :   NA/NaN/Inf in foreign function call (arg 2) Calls: randomForest -> randomForest.default -> .C
• Property to learn from non-numeric
• ‘Inf’ rather than numeric

• Error in randomForest.default(xy[-1], y, ntree = 70, mtry = mtry, importance = imp,  :  
• NA not permitted in predictors
• Property to learn from numeric when classifying

• Error in comps[c1, c2] <- round(roc12, digits = 4) : replacement has length zero Calls: print -> genroc
• One of the classes might be present once, making ou-of-bag validation impossible

• Error in read.table(file = file, header = header, sep = sep, quote = quote,  :
    empty beginning of file Calls: readxy -> read.csv -> read.table
• Property to learn from is missing from the data
• Possibly removed using keep / remove properties



• Error in `rownames<-`(`*tmp*`, value = row.names(x)) :  attempt to set rownames on object with no dimensions  Calls: randomForest ... randomForest.default -> is.na -> is.na.data.frame -> rownames<-
• One of observations has an incomplete set of variables, one or more descriptors are missing (n/a) 

• Error in predict.randomForest(model, x, type = "response") :  New factor levels not present in the training data  Calls: predict -> predict.randomForest
• One of observations has a level for  set one of the variables that was not observed in the training set (e.g. present in the training set : 0,1,2,3 ; value in the test set 6)
• Make sure each observation is seen in the training set
• alternatively use continuous variables to describe the datapoints rather than categorical

• Error in withCallingHandlers(expr, warning = function(w) invokeRestart("muffleWarning")) : invalid multibyte string at '<b2>II' (or at '<a0>hydra) Calls: readxy -> cleandata -> FactorOrNumber Execution halted
• Array property present formatted as blabla[1], blabla[2], etc.
• Flatten to single properties (eg turn binary flag on for proeprties present named by value in array property.

Hope this helps anyone when stuck (and that this page is indexed by Google, well probably isn't the case)