set.seed(1)
library(RColorBrewer)
rf<-colorRampPalette(rev(brewer.pal(11,'Spectral')))
r<-rf(512)
x<-c();y<-c()
x[1]<-0;y[1]<-0
n<-100
for(i in 1:(n*5)) {
  sdx<-rbinom(1,1,.02)*50+1
  
  mx<-rnorm(1,x[length(x)],sdx); sx<-runif(1)*.5+.5
  my<-rnorm(1,y[length(y)],sdx); sy<-runif(1)*.5+.5
  x<-c(x,rnorm(mean=mx,sd=sx,n))
  y<-c(y,rnorm(mean=my,sd=sy,n))
  }
df<-data.frame(x,y)

#plot(df,pch=16,col='black',cex=.5)

library(gplots)
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
par(bg='black')
h2<-hist2d(df,col=c('black', r),axes=F)

cls<-kmeans(df,5, nstart=11)
cl<-cls$cluster
prob=sapply(cl, function(x) c(1,1,1000,1,3000)[x])
pr=sum(cl == 3 | cl==5)/length(cl)
ind=sample(1:nrow(df),floor(nrow(df)*(pr)),prob=prob)
df=df[ind,]

adddata<-function(cen) {
  tdf=lapply(cen, function (c) {
    x<-c();y<-c()
    x[1]<-cls$centers[c,1];y[1]<-cls$centers[c,2]
    for(i in 1:20) {
      sdx<-rbinom(1,1,.03)*20+1
      mx<-rnorm(1,x[length(x)],sdx); sx<-runif(1)*.5+.5
      my<-rnorm(1,y[length(y)],sdx); sy<-runif(1)*.5+.5
      x<-c(x,rnorm(mean=mx,sd=sx,100))
      y<-c(y,rnorm(mean=my,sd=sy,100))
    }
  cdf<-data.frame(x,y)
  })
  
  ttdf=tdf[[1]]
  for(i in 2:(length(tdf)))
    ttdf=rbind(ttdf,tdf[[i]])
  ttdf
}

df=rbind(df,adddata(c(3,5)))

h2<-hist2d(df,col=c('black', r),axes=F)

cls<-kmeans(df,9, nstart=11)
cl<-cls$cluster
prob=sapply(cl, function(x) c(10,1,1,1,3000,1,1,1,1)[x])
pr=sum(cl == 1 | cl==3)/length(cl)
ind=sample(1:nrow(df),floor(nrow(df)*(pr)),prob=prob)
df=df[ind,]

df=rbind(df,adddata(c(1,5)))
h2<-hist2d(df,col=c('black', r),axes=F)

cls<-kmeans(df,5, nstart=11)
cl<-cls$cluster
prob=sapply(cl, function(x) c(1,2000,3000,1,1)[x])
pr=sum(cl == 3 | cl==2)/length(cl)
ind=sample(1:nrow(df),floor(nrow(df)*(pr)),prob=prob)
df=df[ind,]

df=rbind(df,adddata(c(2,2)))
h2<-hist2d(df,col=c('black', r))

cls<-kmeans(df,7, nstart=11)
cl<-cls$cluster
prob=sapply(cl, function(x) c(3000,3000,3000,1,1,1,1)[x])
pr=sum(cl==1 | cl==2 | cl==3)/length(cl)
ind=sample(1:nrow(df),floor(nrow(df)*(pr)),prob=prob)
df=df[ind,]

adddata1<-function(cen) {
  tdf=lapply(cen, function (c) {
    x<-c();y<-c()
    x[1]<-cls$centers[c,1];y[1]<-cls$centers[c,2]
    for(i in 1:20) {
      sdx<-rbinom(1,1,.03)*1+1
      mx<-rnorm(1,x[length(x)],sdx); sx<-runif(1)*.5+.5
      my<-rnorm(1,y[length(y)],sdx); sy<-runif(1)*.5+.5
      x<-c(x,rnorm(mean=mx,sd=sx,100))
      y<-c(y,rnorm(mean=my,sd=sy,100))
    }
    cdf<-data.frame(x,y)
  })
  
  ttdf=tdf[[1]]
  for(i in 2:(length(tdf)))
    ttdf=rbind(ttdf,tdf[[i]])
  ttdf
}

df=rbind(df,adddata1(c(3,2,1,3,3,2,1,2,3)))
h2<-hist2d(df,col=c('black', r))