List of packages to deploy:

library(tidyverse)
library(rio)
library(ggfortify)
library(ggforce)
library(GGally)
library(ggsci)
library(RColorBrewer)
library(cowplot)
library(cluster)
library(factoextra)
library(randomForest)
library(randomForestExplainer)

Functions

binwidth_FD <- function(x){
  2*IQR(x, na.rm = T)*sum(!is.na(x))^(-1/3)}
getmode <- function(v){
  d <- density(v, na.rm = T)
  d$x[which.max(d$y)]
}
quantInv <- function(distr, value) ecdf(distr)(value)
df <- import("dataset_preprocessed.csv") %>% 
  mutate_at(vars(c(kg,kb)), funs(log(.,10))) %>% select(-contains("V"))
raw_data <- import("dataset_raw.csv")
moe <- import("moe_descriptors.csv")

1. Cluster analysis

1.1. Model \(G = k_g \Delta C^g\) (mod_1)

  • Standardization of kg and g, and selection of cluster by silhouette.
mod_1 <- left_join(df, moe, by = c("solute"="compound")) %>% filter(growth.rate %in% c("G=kg?Cg"))
kinetic_parameters_1 <- mod_1 %>% select(c(g,kg)) %>% mutate_all(funs(scale(.)))

fviz_nbclust(kinetic_parameters_1, hcut, method = "silhouette")

  • Based on silhouette criterion, the optimal number of clusters is 3 as it maximizes the index. Thus, AHC was carried out using k = 3 and method = “ward.D” with N = 92, corresponding to 35 solutes.
ahc1 <- hcut(kinetic_parameters_1, 3, hc_method = "ward.D", graph = T) 
fviz_cluster(ahc1, geom = "point") +
  geom_text(aes(label = abbreviate(mod_1$solute, 3, strict = TRUE))) + 
  scale_fill_jco() + scale_color_jco() + theme_classic() + theme(legend.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(y = expression(log~k[g]))

  • Samples were assigned a group and the groups were characterized in term of their values of \(k_g\) and g
mod_1 <- mod_1 %>% mutate(cluster = ahc1$cluster)
mod_1 %>% group_by(cluster) %>% summarise(n(),mean(kg), median(kg), min(kg), max(kg), sd(kg), IQR(kg))
mod_1 %>% group_by(cluster) %>% summarise(n(),mean(g), median(g), min(g), max(g), sd(g),IQR(g))

1.2. Model \(G = k_g (S-1)^g\) (mod_2)

  • Standardization of kg and g, and selection of cluster by silhouette.
mod_2 <- left_join(df, moe, by = c("solute"="compound")) %>% filter(growth.rate %in% c("G=kg(S-1)g"))
kinetic_parameters_2 <- mod_2  %>%  select(c(g,kg)) %>% mutate_all(funs(scale(.)))
fviz_nbclust(kinetic_parameters_2, hcut, method = "silhouette")

  • Based on silhoutte criterion, the optimal number of clusters is 2 as it maximizes the index. However, AHC was carried out using k = 3, given it provides better discrimination and silhoutte index is similar to k = 2. Method = “ward.D2” with N = 68, corresponding to 26 solutes
ahc2 <- hcut(kinetic_parameters_2, 3, hc_method = "ward.D", graph = T) 
fviz_cluster(ahc2, geom = "point") +
  geom_text(aes(label = abbreviate(mod_2$solute, 3, strict = TRUE))) + 
  scale_fill_jco() + scale_color_jco() + theme_classic() + theme(legend.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(y = expression(log~k[g]))

mod_2 <- mod_2 %>% mutate(cluster = ahc2$cluster)
mod_2 %>% group_by(cluster) %>% summarise(n(),mean(kg), median(kg), min(kg), max(kg), sd(kg), IQR(kg))
mod_2 %>% group_by(cluster) %>% summarise(n(),mean(g), median(g), min(g), max(g), sd(g),IQR(g))

1.3. Model \(B = k_b \Delta C^b\) (mod_3)

mod_3 <- left_join(df, moe, by = c("solute"="compound")) %>% filter(nucleation.rate %in% c("B=kb?Cb"))
kinetic_parameters_3 <- mod_3  %>%  select(c(b,kb)) %>% mutate_all(funs(scale(.)))
fviz_nbclust(kinetic_parameters_3, hcut, method = "silhouette")

ahc3 <- hcut(kinetic_parameters_3, 5, hc_method = "ward.D", graph = T) 
fviz_cluster(ahc3, geom = "point") +
  geom_text(aes(label = abbreviate(mod_3$solute, 3, strict = TRUE))) + 
  scale_fill_jco() + scale_color_jco() + theme_classic() + theme(legend.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(y = expression(log~k[b]))

mod_3 <- mod_3 %>% mutate(cluster = ahc3$cluster)
mod_3 %>% group_by(cluster) %>% summarise(n(),mean(kb), median(kb), min(kb), max(kb), sd(kb), IQR(kb))
mod_3 %>% group_by(cluster) %>% summarise(n(),mean(b), median(b), min(kb), max(b), sd(b), IQR(b))

2. Random Forest

2.1. Model \(G = k_g \Delta C^g\)

Co-crystallization datum is removed since descriptors are limited to pure substances (agomelatine-citric acid and caffeine-glutaric acid) (6 data point). RF is run using mtry = 10, ntree = 15000, set.seed = 50

system <- raw_data %>% select(solvent, id) %>% 
  mutate_at(vars(solvent), funs(if_else(str_detect(., "water:"), "aqueous", .)))

excluded_vars <- c("id", "solute", "solute.mw", "density", "b", "kb",   "g", "kg",  "En",   "Eg",   "growth.rate",  "nucleation.rate",  "driving.force", "units.driving.force", "comments", "bidimensional")
mod_1.rf <- mod_1 %>% 
  filter(!(id %in% c(18:22,69) | seeding == "both")) %>% 
  left_join(., system, by = "id") %>% 
  select(-any_of(excluded_vars)) %>% select(c(cluster, solvent, everything())) %>%  
  mutate_at(vars(cluster, solvent, class, method, seeding), funs(as.factor(.)))

set.seed(50)
renamed_vars <- quo(paste0("x", 1:length(mod_1.rf)))
rf1 <- randomForest(x1~., data = mod_1.rf %>% rename_all(renamed_vars), ntree = 15000, mtry = 10, importance = T, localimp = T)
rf1

Call:
 randomForest(formula = x1 ~ ., data = mod_1.rf %>% rename_all(renamed_vars),      ntree = 15000, mtry = 10, importance = T, localimp = T) 
               Type of random forest: classification
                     Number of trees: 15000
No. of variables tried at each split: 10

        OOB estimate of  error rate: 25.88%
Confusion matrix:
  1  2  3 class.error
1 7  0  4   0.3636364
2 1 16  6   0.3043478
3 5  6 40   0.2156863
  • Cross-validation
rf1.cv <- data.frame(actual = rep(NA, nrow(mod_1.rf)), predicted = rep(NA, nrow(mod_1.rf)))

for (i in 1:nrow(mod_1.rf)) {
  set.seed(50)
  rf <- randomForest(x1~., data = (mod_1.rf %>% rename_all(renamed_vars))[-i,], ntree = 15000, mtry = 10)

  rf1.cv[i,1] <- mod_1.rf[i,1]
  rf1.cv[i,2] <- predict(rf, (mod_1.rf %>% rename_all(renamed_vars))[i,])

}
sum(diag(table(rf1.cv)))/nrow(rf1.cv)
[1] 0.7411765
  • Importance by MDA
imp_est_1 <- data.frame(variable = colnames(mod_1.rf)[-1] , MDA = rf1$importance[,4])
imp_est_1 %>% 
  arrange(desc(MDA)) %>% filter(row_number() < 16) %>% ggplot(aes(reorder(variable, MDA), MDA)) + 
  geom_col() + coord_flip() + labs(y = "MDA", x = "Variable")

2.2. Model \(G = k_g (S-1)^g\)

  • Bidimentional growth data were removed
mod_2.rf <- mod_2 %>% 
  filter(bidimensional != "yes", seeding != "both") %>% 
  left_join(., system, by = "id") %>% 
  select(-any_of(excluded_vars)) %>% select(c(cluster, solvent, everything())) %>%  
  mutate_at(vars(cluster, solvent, class, method, seeding), funs(as.factor(.))) %>% 
  na.omit()

set.seed(50)
renamed_vars <- quo(paste0("x", 1:length(mod_2.rf)))
rf2 <- randomForest(x1~., data = mod_2.rf %>% rename_all(renamed_vars), ntree = 15000, mtry = 10, importance = T, localimp = T)
rf2

Call:
 randomForest(formula = x1 ~ ., data = mod_2.rf %>% rename_all(renamed_vars),      ntree = 15000, mtry = 10, importance = T, localimp = T) 
               Type of random forest: classification
                     Number of trees: 15000
No. of variables tried at each split: 10

        OOB estimate of  error rate: 14.55%
Confusion matrix:
   1 2 3 class.error
1 35 4 0   0.1025641
2  3 7 0   0.3000000
3  1 0 5   0.1666667
  • Cross-validation
rf2.cv <- data.frame(actual = rep(NA, nrow(mod_2.rf)), predicted = rep(NA, nrow(mod_2.rf)))

for (i in 1:nrow(mod_2.rf)) {
  set.seed(50)
  rf <- randomForest(x1~., data = (mod_2.rf %>% rename_all(renamed_vars))[-i,], ntree = 15000, mtry = 10)

  rf2.cv[i,1] <- mod_2.rf[i,1]
  rf2.cv[i,2] <- predict(rf, (mod_2.rf %>% rename_all(renamed_vars))[i,])

}
sum(diag(table(rf2.cv)))/nrow(rf2.cv)
[1] 0.8545455
  • Importance by MDA
imp_est_2 <- data.frame(variable = colnames(mod_2.rf)[-1] , MDA = rf2$importance[,4])
imp_est_2 %>% 
  arrange(desc(MDA)) %>% filter(row_number() < 16) %>% ggplot(aes(reorder(variable, MDA), MDA)) + 
  geom_col() + coord_flip() + labs(y = "MDA", x = "Variable")

2.3. Model \(B = k_b \Delta C^g\)

mod_3.rf <- mod_3 %>% 
  filter(!(id %in% c(18,19))) %>% 
  left_join(., system, by = "id") %>% 
  select(-any_of(excluded_vars)) %>% select(c(cluster, solvent, everything())) %>%  
  mutate_at(vars(cluster, solvent, class, method, seeding), funs(as.factor(.)))

set.seed(50)
renamed_vars <- quo(paste0("x", 1:length(mod_3.rf)))
rf3 <- randomForest(x1~., data = mod_3.rf %>% rename_all(renamed_vars), ntree = 15000, mtry = 10, importance = T, localimp = T)
rf3

Call:
 randomForest(formula = x1 ~ ., data = mod_3.rf %>% rename_all(renamed_vars),      ntree = 15000, mtry = 10, importance = T, localimp = T) 
               Type of random forest: classification
                     Number of trees: 15000
No. of variables tried at each split: 10

        OOB estimate of  error rate: 16.95%
Confusion matrix:
   1 2 4 5 class.error
1 33 1 0 0  0.02941176
2  0 8 0 0  0.00000000
4  7 0 2 0  0.77777778
5  2 0 0 6  0.25000000
  • Cross-validation
rf3.cv <- data.frame(actual = rep(NA, nrow(mod_3.rf)), predicted = rep(NA, nrow(mod_3.rf)))

for (i in 1:nrow(mod_3.rf)) {
  set.seed(50)
  rf <- randomForest(x1~., data = (mod_3.rf %>% rename_all(renamed_vars))[-i,], ntree = 15000, mtry = 10)

  rf3.cv[i,1] <- mod_3.rf[i,1]
  rf3.cv[i,2] <- predict(rf, (mod_3.rf %>% rename_all(renamed_vars))[i,])

}
sum(diag(table(rf3.cv)))/nrow(rf3.cv)
[1] 0.8305085
  • Importance by MDA
imp_est_3 <- data.frame(variable = colnames(mod_3.rf)[-1] , MDA = rf3$importance[,4])
imp_est_3 %>% 
  arrange(desc(MDA)) %>% filter(row_number() < 16) %>% ggplot(aes(reorder(variable, MDA), MDA)) + 
  geom_col() + coord_flip() + labs(y = "MDA", x = "Variable")

LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQpMaXN0IG9mIHBhY2thZ2VzIHRvIGRlcGxveToNCg0KYGBge3IsIHdhcm5pbmc9RkFMU0V9DQpsaWJyYXJ5KHRpZHl2ZXJzZSkNCmxpYnJhcnkocmlvKQ0KbGlicmFyeShnZ2ZvcnRpZnkpDQpsaWJyYXJ5KGdnZm9yY2UpDQpsaWJyYXJ5KEdHYWxseSkNCmxpYnJhcnkoZ2dzY2kpDQpsaWJyYXJ5KFJDb2xvckJyZXdlcikNCmxpYnJhcnkoY293cGxvdCkNCmxpYnJhcnkoY2x1c3RlcikNCmxpYnJhcnkoZmFjdG9leHRyYSkNCmxpYnJhcnkocmFuZG9tRm9yZXN0KQ0KbGlicmFyeShyYW5kb21Gb3Jlc3RFeHBsYWluZXIpDQoNCmBgYA0KRnVuY3Rpb25zDQpgYGB7cn0NCmJpbndpZHRoX0ZEIDwtIGZ1bmN0aW9uKHgpew0KICAyKklRUih4LCBuYS5ybSA9IFQpKnN1bSghaXMubmEoeCkpXigtMS8zKX0NCmdldG1vZGUgPC0gZnVuY3Rpb24odil7DQogIGQgPC0gZGVuc2l0eSh2LCBuYS5ybSA9IFQpDQogIGQkeFt3aGljaC5tYXgoZCR5KV0NCn0NCnF1YW50SW52IDwtIGZ1bmN0aW9uKGRpc3RyLCB2YWx1ZSkgZWNkZihkaXN0cikodmFsdWUpDQpgYGANCg0KKiBMb2FkaW5nIGRhdGEgYW5kIG1vbGVjdWxhciBkZXNjcmlwdG9ycw0KKiBDaGFuZ2Ugb2Ygc2NhbGUgb2Yga2IgYW5kIGtnDQoNCmBgYHtyLCB3YXJuaW5nPUZBTFNFfQ0KZGYgPC0gaW1wb3J0KCJkYXRhc2V0X3ByZXByb2Nlc3NlZC5jc3YiKSAlPiUgDQogIG11dGF0ZV9hdCh2YXJzKGMoa2csa2IpKSwgZnVucyhsb2coLiwxMCkpKSAlPiUgc2VsZWN0KC1jb250YWlucygiViIpKQ0KcmF3X2RhdGEgPC0gaW1wb3J0KCJkYXRhc2V0X3Jhdy5jc3YiKQ0KbW9lIDwtIGltcG9ydCgibW9lX2Rlc2NyaXB0b3JzLmNzdiIpDQpgYGANCg0KIyAxLiBDbHVzdGVyIGFuYWx5c2lzDQoNCiMjIDEuMS4gTW9kZWwgJEcgPSBrX2cgXERlbHRhIENeZyQgKG1vZF8xKQ0KDQoqIFN0YW5kYXJkaXphdGlvbiBvZiBrZyBhbmQgZywgYW5kIHNlbGVjdGlvbiBvZiBjbHVzdGVyIGJ5IHNpbGhvdWV0dGUuDQoNCmBgYHtyfQ0KbW9kXzEgPC0gbGVmdF9qb2luKGRmLCBtb2UsIGJ5ID0gYygic29sdXRlIj0iY29tcG91bmQiKSkgJT4lIGZpbHRlcihncm93dGgucmF0ZSAlaW4lIGMoIkc9a2c/Q2ciKSkNCmtpbmV0aWNfcGFyYW1ldGVyc18xIDwtIG1vZF8xICU+JSBzZWxlY3QoYyhnLGtnKSkgJT4lIG11dGF0ZV9hbGwoZnVucyhzY2FsZSguKSkpDQoNCmZ2aXpfbmJjbHVzdChraW5ldGljX3BhcmFtZXRlcnNfMSwgaGN1dCwgbWV0aG9kID0gInNpbGhvdWV0dGUiKQ0KYGBgDQoqIEJhc2VkIG9uIHNpbGhvdWV0dGUgY3JpdGVyaW9uLCB0aGUgb3B0aW1hbCBudW1iZXIgb2YgY2x1c3RlcnMgaXMgMyBhcyBpdCBtYXhpbWl6ZXMgdGhlIGluZGV4LiBUaHVzLCBBSEMgd2FzIGNhcnJpZWQgb3V0IHVzaW5nIGsgPSAzIGFuZCBtZXRob2QgPSDigJx3YXJkLkTigJ0gd2l0aCBOID0gOTIsIGNvcnJlc3BvbmRpbmcgdG8gMzUgc29sdXRlcy4NCg0KYGBge3IsIHdhcm5pbmc9RkFMU0V9DQphaGMxIDwtIGhjdXQoa2luZXRpY19wYXJhbWV0ZXJzXzEsIDMsIGhjX21ldGhvZCA9ICJ3YXJkLkQiLCBncmFwaCA9IFQpIA0KZnZpel9jbHVzdGVyKGFoYzEsIGdlb20gPSAicG9pbnQiKSArDQogIGdlb21fdGV4dChhZXMobGFiZWwgPSBhYmJyZXZpYXRlKG1vZF8xJHNvbHV0ZSwgMywgc3RyaWN0ID0gVFJVRSkpKSArIA0KICBzY2FsZV9maWxsX2pjbygpICsgc2NhbGVfY29sb3JfamNvKCkgKyB0aGVtZV9jbGFzc2ljKCkgKyB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAidG9wIikgKw0KICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSAwLCBsaW5ldHlwZSA9ICJkYXNoZWQiKSArDQogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IDAsIGxpbmV0eXBlID0gImRhc2hlZCIpICsNCiAgbGFicyh5ID0gZXhwcmVzc2lvbihsb2d+a1tnXSkpDQoNCmBgYA0KKiBTYW1wbGVzIHdlcmUgYXNzaWduZWQgYSBncm91cCBhbmQgdGhlIGdyb3VwcyB3ZXJlIGNoYXJhY3Rlcml6ZWQgaW4gdGVybSBvZiB0aGVpciB2YWx1ZXMgb2YgJGtfZyQgYW5kIGcNCg0KYGBge3J9DQptb2RfMSA8LSBtb2RfMSAlPiUgbXV0YXRlKGNsdXN0ZXIgPSBhaGMxJGNsdXN0ZXIpDQptb2RfMSAlPiUgZ3JvdXBfYnkoY2x1c3RlcikgJT4lIHN1bW1hcmlzZShuKCksbWVhbihrZyksIG1lZGlhbihrZyksIG1pbihrZyksIG1heChrZyksIHNkKGtnKSwgSVFSKGtnKSkNCm1vZF8xICU+JSBncm91cF9ieShjbHVzdGVyKSAlPiUgc3VtbWFyaXNlKG4oKSxtZWFuKGcpLCBtZWRpYW4oZyksIG1pbihnKSwgbWF4KGcpLCBzZChnKSxJUVIoZykpDQpgYGANCg0KIyMgMS4yLiBNb2RlbCAkRyA9IGtfZyAoUy0xKV5nJCAobW9kXzIpDQoNCiogU3RhbmRhcmRpemF0aW9uIG9mIGtnIGFuZCBnLCBhbmQgc2VsZWN0aW9uIG9mIGNsdXN0ZXIgYnkgc2lsaG91ZXR0ZS4NCg0KYGBge3J9DQptb2RfMiA8LSBsZWZ0X2pvaW4oZGYsIG1vZSwgYnkgPSBjKCJzb2x1dGUiPSJjb21wb3VuZCIpKSAlPiUgZmlsdGVyKGdyb3d0aC5yYXRlICVpbiUgYygiRz1rZyhTLTEpZyIpKQ0Ka2luZXRpY19wYXJhbWV0ZXJzXzIgPC0gbW9kXzIgICU+JSAgc2VsZWN0KGMoZyxrZykpICU+JSBtdXRhdGVfYWxsKGZ1bnMoc2NhbGUoLikpKQ0KZnZpel9uYmNsdXN0KGtpbmV0aWNfcGFyYW1ldGVyc18yLCBoY3V0LCBtZXRob2QgPSAic2lsaG91ZXR0ZSIpDQpgYGANCiogQmFzZWQgb24gc2lsaG91dHRlIGNyaXRlcmlvbiwgdGhlIG9wdGltYWwgbnVtYmVyIG9mIGNsdXN0ZXJzIGlzIDIgYXMgaXQgbWF4aW1pemVzIHRoZSBpbmRleC4gSG93ZXZlciwgQUhDIHdhcyBjYXJyaWVkIG91dCB1c2luZyBrID0gMywgZ2l2ZW4gaXQgcHJvdmlkZXMgYmV0dGVyIGRpc2NyaW1pbmF0aW9uIGFuZCBzaWxob3V0dGUgaW5kZXggaXMgc2ltaWxhciB0byBrID0gMi4gTWV0aG9kID0g4oCcd2FyZC5EMuKAnSB3aXRoIE4gPSA2OCwgY29ycmVzcG9uZGluZyB0byAyNiBzb2x1dGVzDQoNCmBgYHtyLCB3YXJuaW5nPUZBTFNFfQ0KYWhjMiA8LSBoY3V0KGtpbmV0aWNfcGFyYW1ldGVyc18yLCAzLCBoY19tZXRob2QgPSAid2FyZC5EIiwgZ3JhcGggPSBUKSANCmZ2aXpfY2x1c3RlcihhaGMyLCBnZW9tID0gInBvaW50IikgKw0KICBnZW9tX3RleHQoYWVzKGxhYmVsID0gYWJicmV2aWF0ZShtb2RfMiRzb2x1dGUsIDMsIHN0cmljdCA9IFRSVUUpKSkgKyANCiAgc2NhbGVfZmlsbF9qY28oKSArIHNjYWxlX2NvbG9yX2pjbygpICsgdGhlbWVfY2xhc3NpYygpICsgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gInRvcCIpICsNCiAgZ2VvbV92bGluZSh4aW50ZXJjZXB0ID0gMCwgbGluZXR5cGUgPSAiZGFzaGVkIikgKw0KICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSAwLCBsaW5ldHlwZSA9ICJkYXNoZWQiKSArDQogIGxhYnMoeSA9IGV4cHJlc3Npb24obG9nfmtbZ10pKQ0KYGBgDQpgYGB7cn0NCm1vZF8yIDwtIG1vZF8yICU+JSBtdXRhdGUoY2x1c3RlciA9IGFoYzIkY2x1c3RlcikNCm1vZF8yICU+JSBncm91cF9ieShjbHVzdGVyKSAlPiUgc3VtbWFyaXNlKG4oKSxtZWFuKGtnKSwgbWVkaWFuKGtnKSwgbWluKGtnKSwgbWF4KGtnKSwgc2Qoa2cpLCBJUVIoa2cpKQ0KbW9kXzIgJT4lIGdyb3VwX2J5KGNsdXN0ZXIpICU+JSBzdW1tYXJpc2UobigpLG1lYW4oZyksIG1lZGlhbihnKSwgbWluKGcpLCBtYXgoZyksIHNkKGcpLElRUihnKSkNCmBgYA0KIyAxLjMuIE1vZGVsICRCID0ga19iIFxEZWx0YSBDXmIkIChtb2RfMykNCg0KKiBTdGFuZGFyZGl6YXRpb24gb2Yga2IgYW5kIGIsIGFuZCBzZWxlY3Rpb24gb2YgY2x1c3RlciBieSBzaWxob3VldHRlLg0KDQpgYGB7cn0NCm1vZF8zIDwtIGxlZnRfam9pbihkZiwgbW9lLCBieSA9IGMoInNvbHV0ZSI9ImNvbXBvdW5kIikpICU+JSBmaWx0ZXIobnVjbGVhdGlvbi5yYXRlICVpbiUgYygiQj1rYj9DYiIpKQ0Ka2luZXRpY19wYXJhbWV0ZXJzXzMgPC0gbW9kXzMgICU+JSAgc2VsZWN0KGMoYixrYikpICU+JSBtdXRhdGVfYWxsKGZ1bnMoc2NhbGUoLikpKQ0KZnZpel9uYmNsdXN0KGtpbmV0aWNfcGFyYW1ldGVyc18zLCBoY3V0LCBtZXRob2QgPSAic2lsaG91ZXR0ZSIpDQpgYGANCg0KKiBCYXNlZCBvbiBzaWxob3VldHRlIGNyaXRlcmlvbiwgdGhlIG9wdGltYWwgbnVtYmVyIG9mIGNsdXN0ZXJzIGlzIDIgYXMgaXQgbWF4aW1pemVzIHRoZSBpbmRleDsgYnV0LCBrID0gNSBpcyBhbHNvIGdvb2Qgc2luY2UgdGhlIGluZGV4IGlzIGp1c3Qgc2xpZ2h0bHkgbG93ZXIgdGhhbiBrID0gMiBhbmQgcHJvdmlkZSBhIGJldHRlciBzZXBhcmF0aW9uLiBUaHVzLCBBSEMgd2FzIGNhcnJpZWQgb3V0IHVzaW5nIGsgPSA1IGFuZCBtZXRob2QgPSDigJx3YXJkLkTigJ0gd2l0aCBOID0gNjEsIGNvcnJlc3BvbmRpbmcgdG8gMjEgc29sdXRlcy4NCg0KYGBge3IsIHdhcm5pbmc9RkFMU0V9DQphaGMzIDwtIGhjdXQoa2luZXRpY19wYXJhbWV0ZXJzXzMsIDUsIGhjX21ldGhvZCA9ICJ3YXJkLkQiLCBncmFwaCA9IFQpIA0KZnZpel9jbHVzdGVyKGFoYzMsIGdlb20gPSAicG9pbnQiKSArDQogIGdlb21fdGV4dChhZXMobGFiZWwgPSBhYmJyZXZpYXRlKG1vZF8zJHNvbHV0ZSwgMywgc3RyaWN0ID0gVFJVRSkpKSArIA0KICBzY2FsZV9maWxsX2pjbygpICsgc2NhbGVfY29sb3JfamNvKCkgKyB0aGVtZV9jbGFzc2ljKCkgKyB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAidG9wIikgKw0KICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSAwLCBsaW5ldHlwZSA9ICJkYXNoZWQiKSArDQogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IDAsIGxpbmV0eXBlID0gImRhc2hlZCIpICsNCiAgbGFicyh5ID0gZXhwcmVzc2lvbihsb2d+a1tiXSkpDQpgYGANCg0KKiBTYW1wbGVzIHdlcmUgYXNzaWduZWQgYSBncm91cCBhbmQgdGhlIGdyb3VwcyB3ZXJlIGNoYXJhY3Rlcml6ZWQgaW4gdGVybSBvZiB0aGVpciB2YWx1ZXMgb2YgJGtfYiQgYW5kIGINCg0KYGBge3IsIHdhcm5pbmc9RkFMU0V9DQptb2RfMyA8LSBtb2RfMyAlPiUgbXV0YXRlKGNsdXN0ZXIgPSBhaGMzJGNsdXN0ZXIpDQptb2RfMyAlPiUgZ3JvdXBfYnkoY2x1c3RlcikgJT4lIHN1bW1hcmlzZShuKCksbWVhbihrYiksIG1lZGlhbihrYiksIG1pbihrYiksIG1heChrYiksIHNkKGtiKSwgSVFSKGtiKSkNCm1vZF8zICU+JSBncm91cF9ieShjbHVzdGVyKSAlPiUgc3VtbWFyaXNlKG4oKSxtZWFuKGIpLCBtZWRpYW4oYiksIG1pbihrYiksIG1heChiKSwgc2QoYiksIElRUihiKSkNCmBgYA0KIyAyLiBSYW5kb20gRm9yZXN0DQoNCiMjIDIuMS4gTW9kZWwgJEcgPSBrX2cgXERlbHRhIENeZyQNCg0KQ28tY3J5c3RhbGxpemF0aW9uIGRhdHVtIGlzIHJlbW92ZWQgc2luY2UgZGVzY3JpcHRvcnMgYXJlIGxpbWl0ZWQgdG8gcHVyZSBzdWJzdGFuY2VzIChhZ29tZWxhdGluZS1jaXRyaWMgYWNpZCBhbmQgY2FmZmVpbmUtZ2x1dGFyaWMgYWNpZCkgKDYgZGF0YSBwb2ludCkuIFJGIGlzIHJ1biB1c2luZyBtdHJ5ID0gMTAsIG50cmVlID0gMTUwMDAsIHNldC5zZWVkID0gNTANCg0KYGBge3J9DQpzeXN0ZW0gPC0gcmF3X2RhdGEgJT4lIHNlbGVjdChzb2x2ZW50LCBpZCkgJT4lIA0KICBtdXRhdGVfYXQodmFycyhzb2x2ZW50KSwgZnVucyhpZl9lbHNlKHN0cl9kZXRlY3QoLiwgIndhdGVyOiIpLCAiYXF1ZW91cyIsIC4pKSkNCg0KZXhjbHVkZWRfdmFycyA8LSBjKCJpZCIsICJzb2x1dGUiLCAic29sdXRlLm13IiwgImRlbnNpdHkiLCAiYiIsCSJrYiIsCSJnIiwgImtnIiwJIkVuIiwJIkVnIiwJImdyb3d0aC5yYXRlIiwJIm51Y2xlYXRpb24ucmF0ZSIsCSJkcml2aW5nLmZvcmNlIiwgInVuaXRzLmRyaXZpbmcuZm9yY2UiLAkiY29tbWVudHMiLAkiYmlkaW1lbnNpb25hbCIpDQpgYGANCg0KYGBge3IsIHdhcm5pbmc9RkFMU0V9DQptb2RfMS5yZiA8LSBtb2RfMSAlPiUgDQogIGZpbHRlcighKGlkICVpbiUgYygxODoyMiw2OSkgfCBzZWVkaW5nID09ICJib3RoIikpICU+JSANCiAgbGVmdF9qb2luKC4sIHN5c3RlbSwgYnkgPSAiaWQiKSAlPiUgDQogIHNlbGVjdCgtYW55X29mKGV4Y2x1ZGVkX3ZhcnMpKSAlPiUgc2VsZWN0KGMoY2x1c3Rlciwgc29sdmVudCwgZXZlcnl0aGluZygpKSkgJT4lICANCiAgbXV0YXRlX2F0KHZhcnMoY2x1c3Rlciwgc29sdmVudCwgY2xhc3MsIG1ldGhvZCwgc2VlZGluZyksIGZ1bnMoYXMuZmFjdG9yKC4pKSkNCg0Kc2V0LnNlZWQoNTApDQpyZW5hbWVkX3ZhcnMgPC0gcXVvKHBhc3RlMCgieCIsIDE6bGVuZ3RoKG1vZF8xLnJmKSkpDQpyZjEgPC0gcmFuZG9tRm9yZXN0KHgxfi4sIGRhdGEgPSBtb2RfMS5yZiAlPiUgcmVuYW1lX2FsbChyZW5hbWVkX3ZhcnMpLCBudHJlZSA9IDE1MDAwLCBtdHJ5ID0gMTAsIGltcG9ydGFuY2UgPSBULCBsb2NhbGltcCA9IFQpDQpyZjENCg0KYGBgDQoqIENyb3NzLXZhbGlkYXRpb24NCg0KYGBge3J9DQpyZjEuY3YgPC0gZGF0YS5mcmFtZShhY3R1YWwgPSByZXAoTkEsIG5yb3cobW9kXzEucmYpKSwgcHJlZGljdGVkID0gcmVwKE5BLCBucm93KG1vZF8xLnJmKSkpDQoNCmZvciAoaSBpbiAxOm5yb3cobW9kXzEucmYpKSB7DQogIHNldC5zZWVkKDUwKQ0KICByZiA8LSByYW5kb21Gb3Jlc3QoeDF+LiwgZGF0YSA9IChtb2RfMS5yZiAlPiUgcmVuYW1lX2FsbChyZW5hbWVkX3ZhcnMpKVstaSxdLCBudHJlZSA9IDE1MDAwLCBtdHJ5ID0gMTApDQoNCiAgcmYxLmN2W2ksMV0gPC0gbW9kXzEucmZbaSwxXQ0KICByZjEuY3ZbaSwyXSA8LSBwcmVkaWN0KHJmLCAobW9kXzEucmYgJT4lIHJlbmFtZV9hbGwocmVuYW1lZF92YXJzKSlbaSxdKQ0KDQp9DQpzdW0oZGlhZyh0YWJsZShyZjEuY3YpKSkvbnJvdyhyZjEuY3YpDQpgYGANCiogSW1wb3J0YW5jZSBieSBNREENCg0KYGBge3J9DQppbXBfZXN0XzEgPC0gZGF0YS5mcmFtZSh2YXJpYWJsZSA9IGNvbG5hbWVzKG1vZF8xLnJmKVstMV0gLCBNREEgPSByZjEkaW1wb3J0YW5jZVssNF0pDQppbXBfZXN0XzEgJT4lIA0KICBhcnJhbmdlKGRlc2MoTURBKSkgJT4lIGZpbHRlcihyb3dfbnVtYmVyKCkgPCAxNikgJT4lIGdncGxvdChhZXMocmVvcmRlcih2YXJpYWJsZSwgTURBKSwgTURBKSkgKyANCiAgZ2VvbV9jb2woKSArIGNvb3JkX2ZsaXAoKSArIGxhYnMoeSA9ICJNREEiLCB4ID0gIlZhcmlhYmxlIikNCg0KYGBgDQoNCiMjIDIuMi4gTW9kZWwgJEcgPSBrX2cgKFMtMSleZyQNCg0KKiBCaWRpbWVudGlvbmFsIGdyb3d0aCBkYXRhIHdlcmUgcmVtb3ZlZA0KDQpgYGB7cn0NCm1vZF8yLnJmIDwtIG1vZF8yICU+JSANCiAgZmlsdGVyKGJpZGltZW5zaW9uYWwgIT0gInllcyIsIHNlZWRpbmcgIT0gImJvdGgiKSAlPiUgDQogIGxlZnRfam9pbiguLCBzeXN0ZW0sIGJ5ID0gImlkIikgJT4lIA0KICBzZWxlY3QoLWFueV9vZihleGNsdWRlZF92YXJzKSkgJT4lIHNlbGVjdChjKGNsdXN0ZXIsIHNvbHZlbnQsIGV2ZXJ5dGhpbmcoKSkpICU+JSAgDQogIG11dGF0ZV9hdCh2YXJzKGNsdXN0ZXIsIHNvbHZlbnQsIGNsYXNzLCBtZXRob2QsIHNlZWRpbmcpLCBmdW5zKGFzLmZhY3RvciguKSkpICU+JSANCiAgbmEub21pdCgpDQoNCnNldC5zZWVkKDUwKQ0KcmVuYW1lZF92YXJzIDwtIHF1byhwYXN0ZTAoIngiLCAxOmxlbmd0aChtb2RfMi5yZikpKQ0KcmYyIDwtIHJhbmRvbUZvcmVzdCh4MX4uLCBkYXRhID0gbW9kXzIucmYgJT4lIHJlbmFtZV9hbGwocmVuYW1lZF92YXJzKSwgbnRyZWUgPSAxNTAwMCwgbXRyeSA9IDEwLCBpbXBvcnRhbmNlID0gVCwgbG9jYWxpbXAgPSBUKQ0KcmYyDQoNCmBgYA0KKiBDcm9zcy12YWxpZGF0aW9uDQoNCmBgYHtyfQ0KcmYyLmN2IDwtIGRhdGEuZnJhbWUoYWN0dWFsID0gcmVwKE5BLCBucm93KG1vZF8yLnJmKSksIHByZWRpY3RlZCA9IHJlcChOQSwgbnJvdyhtb2RfMi5yZikpKQ0KDQpmb3IgKGkgaW4gMTpucm93KG1vZF8yLnJmKSkgew0KICBzZXQuc2VlZCg1MCkNCiAgcmYgPC0gcmFuZG9tRm9yZXN0KHgxfi4sIGRhdGEgPSAobW9kXzIucmYgJT4lIHJlbmFtZV9hbGwocmVuYW1lZF92YXJzKSlbLWksXSwgbnRyZWUgPSAxNTAwMCwgbXRyeSA9IDEwKQ0KDQogIHJmMi5jdltpLDFdIDwtIG1vZF8yLnJmW2ksMV0NCiAgcmYyLmN2W2ksMl0gPC0gcHJlZGljdChyZiwgKG1vZF8yLnJmICU+JSByZW5hbWVfYWxsKHJlbmFtZWRfdmFycykpW2ksXSkNCg0KfQ0Kc3VtKGRpYWcodGFibGUocmYyLmN2KSkpL25yb3cocmYyLmN2KQ0KYGBgDQoqIEltcG9ydGFuY2UgYnkgTURBDQoNCmBgYHtyfQ0KaW1wX2VzdF8yIDwtIGRhdGEuZnJhbWUodmFyaWFibGUgPSBjb2xuYW1lcyhtb2RfMi5yZilbLTFdICwgTURBID0gcmYyJGltcG9ydGFuY2VbLDRdKQ0KaW1wX2VzdF8yICU+JSANCiAgYXJyYW5nZShkZXNjKE1EQSkpICU+JSBmaWx0ZXIocm93X251bWJlcigpIDwgMTYpICU+JSBnZ3Bsb3QoYWVzKHJlb3JkZXIodmFyaWFibGUsIE1EQSksIE1EQSkpICsgDQogIGdlb21fY29sKCkgKyBjb29yZF9mbGlwKCkgKyBsYWJzKHkgPSAiTURBIiwgeCA9ICJWYXJpYWJsZSIpDQoNCmBgYA0KDQojIyAyLjMuIE1vZGVsICRCID0ga19iIFxEZWx0YSBDXmckDQoNCmBgYHtyfQ0KbW9kXzMucmYgPC0gbW9kXzMgJT4lIA0KICBmaWx0ZXIoIShpZCAlaW4lIGMoMTgsMTkpKSkgJT4lIA0KICBsZWZ0X2pvaW4oLiwgc3lzdGVtLCBieSA9ICJpZCIpICU+JSANCiAgc2VsZWN0KC1hbnlfb2YoZXhjbHVkZWRfdmFycykpICU+JSBzZWxlY3QoYyhjbHVzdGVyLCBzb2x2ZW50LCBldmVyeXRoaW5nKCkpKSAlPiUgIA0KICBtdXRhdGVfYXQodmFycyhjbHVzdGVyLCBzb2x2ZW50LCBjbGFzcywgbWV0aG9kLCBzZWVkaW5nKSwgZnVucyhhcy5mYWN0b3IoLikpKQ0KDQpzZXQuc2VlZCg1MCkNCnJlbmFtZWRfdmFycyA8LSBxdW8ocGFzdGUwKCJ4IiwgMTpsZW5ndGgobW9kXzMucmYpKSkNCnJmMyA8LSByYW5kb21Gb3Jlc3QoeDF+LiwgZGF0YSA9IG1vZF8zLnJmICU+JSByZW5hbWVfYWxsKHJlbmFtZWRfdmFycyksIG50cmVlID0gMTUwMDAsIG10cnkgPSAxMCwgaW1wb3J0YW5jZSA9IFQsIGxvY2FsaW1wID0gVCkNCnJmMw0KYGBgDQoqIENyb3NzLXZhbGlkYXRpb24NCmBgYHtyfQ0KcmYzLmN2IDwtIGRhdGEuZnJhbWUoYWN0dWFsID0gcmVwKE5BLCBucm93KG1vZF8zLnJmKSksIHByZWRpY3RlZCA9IHJlcChOQSwgbnJvdyhtb2RfMy5yZikpKQ0KDQpmb3IgKGkgaW4gMTpucm93KG1vZF8zLnJmKSkgew0KICBzZXQuc2VlZCg1MCkNCiAgcmYgPC0gcmFuZG9tRm9yZXN0KHgxfi4sIGRhdGEgPSAobW9kXzMucmYgJT4lIHJlbmFtZV9hbGwocmVuYW1lZF92YXJzKSlbLWksXSwgbnRyZWUgPSAxNTAwMCwgbXRyeSA9IDEwKQ0KDQogIHJmMy5jdltpLDFdIDwtIG1vZF8zLnJmW2ksMV0NCiAgcmYzLmN2W2ksMl0gPC0gcHJlZGljdChyZiwgKG1vZF8zLnJmICU+JSByZW5hbWVfYWxsKHJlbmFtZWRfdmFycykpW2ksXSkNCg0KfQ0Kc3VtKGRpYWcodGFibGUocmYzLmN2KSkpL25yb3cocmYzLmN2KQ0KYGBgDQoNCiogSW1wb3J0YW5jZSBieSBNREENCmBgYHtyfQ0KaW1wX2VzdF8zIDwtIGRhdGEuZnJhbWUodmFyaWFibGUgPSBjb2xuYW1lcyhtb2RfMy5yZilbLTFdICwgTURBID0gcmYzJGltcG9ydGFuY2VbLDRdKQ0KaW1wX2VzdF8zICU+JSANCiAgYXJyYW5nZShkZXNjKE1EQSkpICU+JSBmaWx0ZXIocm93X251bWJlcigpIDwgMTYpICU+JSBnZ3Bsb3QoYWVzKHJlb3JkZXIodmFyaWFibGUsIE1EQSksIE1EQSkpICsgDQogIGdlb21fY29sKCkgKyBjb29yZF9mbGlwKCkgKyBsYWJzKHkgPSAiTURBIiwgeCA9ICJWYXJpYWJsZSIpDQpgYGANCg0KDQo=