.packageName <- "fpc"
tdecomp <- function(m){
  wm <- eigen(m, symmetric=TRUE)
  p <- ncol(m)
  wmd <- wm$values
  for (i in 1:p){
    if (abs(wmd[i])<1e-6)
      wmd[i] <- 1e-6
  }
  out <- t(wm$vectors %*% diag(sqrt(wmd)))
  out
}

discrcoord <- function(xd, clvecd, pool="n", ...) {
  x <- as.matrix(xd)
  clvec <- as.integer(clvecd)
  n <- nrow(x)
  p <- ncol(x)
  clf <- factor(clvec)
  cll <- as.integer(levels(clf)) 
  clnum <- length(cll)
  cln <- rep(0, times=clnum)
  for (i in 1:clnum){
    cln[i] <- sum(clvec==cll[i])
  }
  W <- rep(0, times=p*p)
  dim(W) <- c(p,p)
  for (i in 1:clnum){
    clx <- rep(0, times=p*cln[i])
    dim(clx) <- c(cln[i],p)
    for (j in 1:p){
      clx[,j] <- x[,j][clvec==cll[i]]
    }
    if (pool=="n")
      W <- W + ((cln[i]-1)*cov(clx))
    else
      W <- W + (n-1)*cov(clx)/clnum
  }
  Tm <- tdecomp(W)
  Tinv <- solve(Tm)
  S <- (n-1)*cov(x)
  B <- S-W
  Z <- t(Tinv) %*% B %*% Tinv
  dc <- eigen(Z, symmetric=TRUE)
  units <- Tinv %*% dc$vectors * sqrt(n-clnum)
  proj <- x %*% units    
  list(ev=dc$values, units=units, proj=proj, W=W)
}

batvarcoord <- function(xd, clvecd, clnum=1){
  x <- as.matrix(xd)
  clvec <- as.integer(as.integer(clvecd)==as.integer(clnum))
  n <- nrow(x)
  p <- ncol(x)
  cll <- c(0,1)
  cln <- rep(0, times=2)
  for (i in 1:2)
    cln[i] <- sum(clvec==cll[i])
  clx <- list()
  for (i in 1:2){
    clx[[i]] <- rep(0, times=p*cln[i])
    dim(clx[[i]]) <- c(cln[i],p)
    for (j in 1:p)
      clx[[i]][,j] <- x[,j][clvec==cll[i]]
  }
  S1 <- cov(clx[[1]])
  S2 <- cov(clx[[2]])
  W <- solve(S2) %*% S1
  Weigen <- eigen(W)
  rev <- Weigen$values + 1/Weigen$values + 2
  dw <- diag(t(Weigen$vectors) %*% S2 %*% Weigen$vectors)
  svw <- matrix(ncol=p, nrow=p)
  for (i in 1:p)
    svw[,i] <- Weigen$vectors[,i]/sqrt(dw[i])
  units <- svw[,(order(-rev))]
  proj <- x %*% units    
  list(ev=Weigen$values[order(-rev)], rev=rev[(order(-rev))], 
       units=units, proj=proj, W=W, S1=S1, S2=S2)
}

batcoord <- function(xd, clvecd, clnum=1, dom="mean"){
  x <- as.matrix(xd)
  clvec <- as.integer(clvecd)
  clf <- factor(clvec)
  cll <- as.integer(levels(clf)) 
  clz <- length(cll)
  p <- ncol(x)
  if (clz!=2){
    clvec <- as.integer(clvecd==clnum)
    cll <- c(0,1)
    print("Cluster indicator has more than 2 values")
  }
  if (dom=="mean"){
    dcx <- discrcoord(x, clvec, pool="equal")
    x2 <- dcx$proj[,2:p]
    batx <- batvarcoord(x2, as.integer(clvecd==clnum))
    ev <- c(dcx$ev[1],batx$ev)
    rev <- c(max(batx$rev)+1, batx$rev)
    units <- cbind(dcx$units[,1],dcx$units[,2:p] %*% batx$units)
    proj <- cbind(dcx$proj[,1],batx$proj)
  }
  else{
    batx <- batvarcoord(x, as.integer(clvecd==clnum))
    cln <- rep(0, times=2)
    mx <- matrix(nrow=nrow(x),ncol=p)
    for (i in 1:2)
      cln[i] <- sum(clvec==cll[i])
    for (i in 1:2){
      clx <- rep(0, times=p*cln[i])
      dim(clx) <- c(cln[i],p)
      for (j in 1:p){
        clx[,j] <- x[,j][clvec==cll[i]]
      }
      mx[i,] <- colMeans(clx)
    }
    mdiff <- mx[2,]-mx[1,]
    ev <- batx$ev
    rev <- rep(0, times=p)
    for (i in 1:p)
      rev[i] <- (batx$units[,i] %*% mdiff)^2/(1+ev[i])+log(ev[i]+1/ev[i]+2)
    units <- batx$units[,(order(-rev))]
    proj <- batx$proj[,(order(-rev))]
  }    
  list(ev=ev[order(-rev)], rev=rev[order(-rev)], 
       units=units, proj=proj)
}

# discriminant plot
# bw: black/white
plotcluster <- function(x, clvecd, clnum=1,
                        method=ifelse(identical(range(as.integer(clvecd)),
                          as.integer(c(0,1))),"awc","dc"),bw=FALSE, xlab=NULL,
                        ylab=NULL, pch=NULL, col=NULL, ...){
  asym <- any(method==c("bc","vbc","adc","awc","arc","anc"))
  if (asym)
    clvec <- as.integer(as.integer(clvecd)==as.integer(clnum))
  else
    clvec <- as.integer(clvecd)
  cx <- discrproj(x, clvecd, method, clnum, ...)$proj
  if (is.null(xlab))
    xlab <- paste(method,"1")
  if (is.null(ylab))
    ylab <- paste(method,"2")
  if (is.null(pch))
    pch <- if (bw){
             1+clvec 
           }
           else 1
  if (is.null(col))
    col <- if (bw) 1
           else{
             1+clvec
           }
  plot(cx, xlab=xlab, ylab=ylab, pch=pch, col=col, ...)
}

discrproj <- function(x, clvecd, method="awc", clnum=1, ...){
  result <- switch(method,
                   dc=discrcoord(x, clvecd, ...),
                   bc=batcoord(x, clvecd, clnum),
                   vbc=batcoord(x, clvecd, clnum, dom="var"),
                   mvdc=mvdcoord(x, clvecd, clnum, ...),
                   adc=adcoord(x, clvecd, clnum),
                   awc=awcoord(x, clvecd, clnum, ...),
                   arc=awcoord(x, clvecd, clnum, method="mcd", ...),
                   nc=ncoord(x, clvecd, ...),
                   wnc=ncoord(x, clvecd, weighted=TRUE, ...),
                   anc=ancoord(x, clvecd, clnum, ...))
  result
}
                   
                   







# nearest neighbor pooled linear dimension reduction according to
# Hastie and Tibshirani,
# IEEE Trans. Pattern Analysis and Machine Intelligence 18 (1996), 607-616
ncoord <- function(xd, clvecd, nn=50, weighted=FALSE,
                    sphere="mcd", orderall=TRUE, countmode=1000, ...){
  z <- x <- as.matrix(xd)
  if (is.matrix(sphere)){
    cv <- sphere
    sphere <- "matrix"
  }
  if (sphere!="none"){
    if (as.numeric(R.version$major)<=1 & as.numeric(R.version$minor)<9)
      require(lqs)
    else require(MASS)
    if (sphere=="matrix")
      Sig <- cv
    else
      Sig <- cov.rob(x, method=sphere, nsamp=500)$cov
    Tds <- tdecomp(Sig)
    Y <- solve(Tds)
    z <- x %*% Y
  }
  clvec <- as.integer(clvecd)
  clf <- factor(clvec)
  cll <- as.integer(levels(clf)) 
  cln <- length(cll)
  p <- ncol(z)
  n <- nrow(z)
  B <- matrix(0,ncol=p,nrow=p)
  for (i in 1:n){
    if(countmode*round(i/countmode)==i)
      cat("Processing point ",i," of ",n,"\n")
    Bi <- matrix(0,ncol=p,nrow=p)
    za <- sweep(z,2,z[i,])
    mds <- rowSums(za*za)
    if (orderall)
      omah <- order(mds)[1:nn]
    else{
      omah <- c()
      maxmds <- max(mds)
      for (j in 1:nn){
        argmin <- which.min(mds)
        omah <- c(omah,argmin)
        mds[argmin] <- maxmds+1
      }
    }
    mi <- colMeans(z[omah,])
    for (j in 1:cln){
      nj <- sum(clvec[omah]==cll[j])
      pij <- nj/nn
      v <- if (pij==0) rep(0,p)
           else colMeans(z[omah,][clvec[omah]==cll[j],,drop=FALSE])-mi
      Bi <- Bi+ pij*(v %*% rbind(v))
    }
    if (weighted){
      sb <- sum(diag(Bi))
      if (sb>0)
        Bi <- Bi/sb
    }
    B <- B+Bi
  }
  B <- B/n
  em <- eigen(B, symmetric=TRUE)
  units <- em$vectors
  if (sphere!="none")
    units <- Y %*% units
  proj <- x %*% units 
  list(ev=em$values, units=units, proj=proj)
}

# asymmetric robust
# nearest neighbor pooled linear dimension reduction 
ancoord <- function(xd, clvecd, clnum=1, nn=50, method="mcd",
                    countmode=1000, ...){
  if (as.numeric(R.version$major)<=1 & as.numeric(R.version$minor)<9)
      require(lqs)
  else require(MASS)
  x <- as.matrix(xd)
  p <- ncol(x)
  n <- nrow(x)
  clvec <- as.integer(clvecd)
  dcl <- as.integer(clnum)
  ci <- clvec==dcl
  clxf <- x[ci,]
  nc <- sum(ci)
  quant <- min(floor(3*(nrow(clxf) + ncol(clxf) + 1)/4),nrow(clxf)-2)
  cv <- cov.rob(clxf,quantile.used=quant,method=method,nsamp=500)
  S1 <- cv$cov
  cinv <- solvecov(S1)$inv
  B <- matrix(0,ncol=p,nrow=p)
  w <- 0
  repeat{
    for (i in 1:nc){
      if(countmode*round(i/countmode)==i)
        cat("Processing point ",i," of ",nc,"\n")
      Bi <- matrix(0,ncol=p,nrow=p)
      mds <- mahalanobis(x,center=clxf[i,],cov=cinv,inverted=TRUE)
      omah <- order(mds)[1:nn]
      wi <- 1    
      mi <- colMeans(x[omah,])
      ni <- sum(clvec[omah]==dcl)
      nr <- nn-ni
      wi <- ni*nr
      if (wi>0){
        vi <- colMeans(x[omah,][ci[omah],,drop=FALSE])-mi
        vr <- colMeans(x[omah,][!ci[omah],,drop=FALSE])-mi
        Bi <- Bi+ ni*(vi %*% rbind(vi))+ nr*(vr %*% rbind(vr))
      }
      sb <- sum(diag(Bi))
      if (sb>0)
        Bi <- Bi/(nn*sb)
      B <- B+Bi
    }
    if (!identical(B,matrix(0,ncol=p,nrow=p)))
      break
    else{
      if (nn<nc+1)
        nn <- nc+1
      else{
        warning("Estimated between groups matrix is zero!")
        break
      }
    }
  }
  Tm <- tdecomp(S1)
  Tinv <- solve(Tm)
  Z <- t(Tinv) %*% B %*% Tinv
  dc <- eigen(Z, symmetric=TRUE)
  units <- Tinv %*% dc$vectors
  proj <- x %*% units    
  list(ev=dc$values, units=units, proj=proj, nn=nn)
}

# quadratic dimension reduction according to Young, Marco and Odell,
# Journal Stat. Plann. Inf. 17 (1986), 307-319; computation according to
# Roehl and Weihs, in Gaul & Locarek-Junge (1999), 253.
mvdcoord <- function(xd, clvecd, clnum=1, sphere="mcd", ...){
  x <- as.matrix(xd)
  if (is.matrix(sphere)){
    cv <- sphere
    sphere="matrix"
  }
  if (as.numeric(R.version$major)<=1 & as.numeric(R.version$minor)<9)
      require(lqs)
  else require(MASS)
  if (sphere!="none"){
    if (sphere=="matrix")
      Sig <- cv
    else
      Sig <- cov.rob(x, method=sphere, nsamp=500)$cov
    Tds <- tdecomp(Sig)
    Y <- solve(Tds)
    z <- x %*% Y
  }
  else
    z <- x
  clvec <- as.integer(clvecd)
  clf <- factor(clvec)
  cll <- as.integer(levels(clf)) 
  clnum <- length(cll)
  p <- ncol(z)
  mx <- vx <- list()
  for (i in 1:clnum){
    mx[[i]] <- colMeans(z[clvecd==cll[i],])
    vx[[i]] <- cov(z[clvecd==cll[i],])
  }
  meandiff <- vardiff <- c()
  for (i in 2:clnum){
    meandiff <- cbind(meandiff,mx[[i]]-mx[[1]])
    vardiff <- cbind(vardiff,vx[[i]]-vx[[1]])
  }
  M <- cbind(meandiff,vardiff)
  em <- eigen(M %*% t(M), symmetric=TRUE)
  units <- em$vectors
  if (sphere!="none")
    units <- Y %*% units
  proj <- x %*% units 
  list(ev=em$values, units=units, proj=proj)
}

# mahalanodisc=vector of mahalanobis distances from n1 points of x1
# to n-n1 points from x2
# modus: see mahal in robcoord
mahalanodisc <- function (x2, mg, covg, modus="square") {
  covinv <- solvecov(covg)$inv
  md <- switch(modus,
         md=sqrt(mahalanobis(x2,mg,covinv,inverted=TRUE)),
         mahalanobis(x2,mg,covinv,inverted=TRUE))
  md
}
# dist:n-n1 Mahalanobis distances, mg: mean(x1), covg: Covariance(x1)

# weight function for robcoord
c.weight <- function(x,ca){
  out <- 1
  if (x > ca)
    out <- ca/x
  out
}

adcoord <- function(xd, clvecd, clnum=1) {
  x <- as.matrix(xd)
  clvec <- as.integer(clvecd)
  dcl <- as.integer(clnum)
  ci <- clvec==dcl
  n <- nrow(x)
  p <- ncol(x)
  cln <- sum(ci)
  clx <- rep(0, times=p*cln)
  clxc <- rep(0, times=p*(n-cln))
  dim(clx) <- c(cln,p)
  dim(clxc) <- c((n-cln),p)
  for (j in 1:p){
    clx[,j] <- x[,j][ci]
    clxc[,j] <- x[,j][!ci]
  }
  S1 <- cov(clx)
  S2 <- cov(clxc)
  W1 <- (cln-1)*S1
  S <- cov(x)
  B <- (n*(n-1)*S - cln*W1 - (n-cln)*(n-cln-1)*S2)/(n-cln)
  Tm <- tdecomp(S1)
  Tinv <- solve(Tm)
  Z <- t(Tinv) %*% B %*% Tinv
  dc <- eigen(Z, symmetric=TRUE)
  units <- Tinv %*% dc$vectors
  proj <- x %*% units    
  list(ev=dc$values, units=units, proj=proj)
}

# "robustifizierte" 1-Cluster-Diskriminanzkoordinaten (durchschnittlicher
# Innerhalb-Abstand vs. Abstand nach ausserhalb, letzterer gewichtet
# mit c(Mahal gross)/Mahal(x_j-mean1).
# Projektionen, Eigenwerte
# x: Daten, clvec: Clusterindikatorvektor, clnum: Nummer des
# zu trennenden Clusters ,
# mahal="square": Squared Mahalanobis distance is used
# mahal="md": Mahalanobis distance is used
# subsample: size of subsample of cluster to use (0=all)
# countmode=output of current point number
awcoord <- function(xd, clvecd, clnum=1, mahal="square", method="classical",
                     clweight=switch(method,classical=FALSE,TRUE), alpha=0.99,
                     subsample=0, countmode=1000, ...) {
  x <- as.matrix(xd)
  if (as.numeric(R.version$major)<=1 & as.numeric(R.version$minor)<9)
      require(lqs)
  else require(MASS)
  n <- nrow(x)
  p <- ncol(x)
  dcl <- as.integer(clnum)
  clfull <- as.integer(clvecd)
  cln <- sum(clfull==dcl)
  if (subsample==0){
    clvec <- as.integer(clfull==dcl)
    cn <- cln
    subs <- NULL
    clxf <- clx <- x[clvec==1,]
  }
  else{
    subs <- sample((1:n)[clfull==dcl],subsample)
    clvec <- 2*(clfull==dcl)
    clvec[subs] <- 1
    clxf <- x[clvec==2,]
    clx <- x[clvec==1,]
    cn <- subsample
  }
  clxc <- x[clvec==0,]
  clxa <- if (clweight) x else clxc
  quant <- min(floor(3*(nrow(clxf) + ncol(clxf) + 1)/4),nrow(clxf)-2)
  cv <- cov.rob(clxf, quantile.used=quant,method=method,nsamp=500)
  S1 <- cv$cov
  mg <- cv$center
  mah <- mahalanodisc(clxa, mg, S1, modus=mahal)
  wg <- switch(mahal,
     md=sapply(mah,c.weight,sqrt(qchisq(alpha,p))),
     sapply(mah,c.weight,qchisq(alpha,p))) 
   if (clweight){
     wg0 <- wg[clvec==0]
     wg1 <- wg[clvec==1]
     wsum0 <- sum(wg0)
     wsum1 <- sum(wg1)
   }
  else{
    wg0 <- wg
    wsum0 <- sum(wg)
  }
   d1 <- d2 <- rep(0,p)
   D12 <- D22 <- matrix(0,ncol=p, nrow=p)
   for(i in 1:cn){
     if (clweight)
       D12 <- D12 + wg1[i]*clx[i,] %*% t(clx[i,])
     else
       D12 <- D12 + clx[i,] %*% t(clx[i,])
     if (clweight)
       d1 <- d1+wg1[i]*clx[i,]
     else
       d1 <- d1+clx[i,]
   }
   for (j in 1:(n-cln)){
     D22 <- D22 + wg0[j]*clxc[j,] %*% t(clxc[j,])
     d2 <- d2+wg0[j]*clxc[j,]
   }
   D21 <- d1 %*% t(d2)
   if (clweight)
     B <- wsum0*D12-D21-t(D21)+wsum1*D22
   else
     B <- wsum0*D12-D21-t(D21)+cn*D22
  if (clweight)
    wsum <- sum(outer(wg0,wg1))
  else              
     wsum <- wsum0*cn
  B <- B/wsum
  Tm <- tdecomp(S1)
  Tinv <- solve(Tm)
  Z <- t(Tinv) %*% B %*% Tinv
  dc <- eigen(Z, symmetric=TRUE)
  units <- Tinv %*% dc$vectors
  proj <- x %*% units    
  list(ev=dc$values, units=units, proj=proj, wg=wg)
}

#
# fixreg utilities
#
# Generation of ca by formula
can <- function (n,p){
    ca <- 3+33/(n*2^(-(p-1)/2))^(1/3)+2900000/(n*2^(-(p-1)/2))^3
  ca
}  
      
# clusexpect= Expectation of times found for a cluster of given size under
# assumption that it will be always found with only its points and never else.
clusexpect <- function (n, p, cn, ir) {
  prob <- choose (cn, p+2)/choose (n, p+2)
  result <- ir*prob
  result
}
 
# Calculates iteration number such that theoretical probability that cluster
# with size cn is found at least mtf times exceeds prob.   
itnumber <- function (n, p, cn, mtf, prob=0.95, maxir=20000){
  minir <- 1
  repeat
  {
     if (qbinom(1-prob,minir,clusexpect(n,p,cn,1))>=mtf){
       ir <- minir
       break
     }
     if (qbinom(1-prob,maxir,clusexpect(n,p,cn,1))<mtf){
       ir <- maxir
       warning("Iteration number too large. maxir chosen.")
       break
     }
     ir <- round((minir+maxir)/2)
     if (qbinom(1-prob,ir,clusexpect(n,p,cn,1))<mtf)
       minir <- ir
     else
       maxir <- ir
     if (maxir <= minir+1)
       break
  }     # repeat
  ir
}

# Minimum size of FPCs that will be found at least mtf times with probability
# prob under ir iteration runs.
minsize <- function(n, p, ir, mtf, prob=0.5){
  mincn <- p+2
  maxcn <- n-1
  repeat
  {
     if (qbinom(1-prob,ir,clusexpect(n,p,mincn,1))>=mtf){
       cn <- mincn
       break
     }
     if (qbinom(1-prob,ir,clusexpect(n,p,maxcn,1))<mtf){
       cn <- maxcn
       warning("Too few iteration runs. Informative FPCs are not probable.")
       break
     }
     cn <- round((mincn+maxcn)/2)
     if (qbinom(1-prob,ir,clusexpect(n,p,cn,1))<mtf)
       mincn <- cn
     else
       maxcn <- cn
     if (maxcn <= mincn+1)
       break
  }     # repeat
  cn
}
     

randconf <- function (n,p){
  gv <- rep(FALSE, times=n)
  m <- sample(n,p)
  gv[m] <- TRUE
  gv
}

#
# Regression FPC iteration;
# output: coef=regression coefficients, var=residual variance, g=indicator vector
rfpi <- function (indep, dep, p, gv, ca, maxit, plot) {
  cachange <- (is.na(ca))
  n <- length(gv)
  ir <- 0
  gva <- rep(FALSE, times=n)
#  print("gva") 
  change <- TRUE
  while (change) {
    ir <- ir+1
    change <- !identical(gv,gva)
    if (ir>maxit){
      change <- FALSE
      warning("Maximum iteration number exceeded.")
    }
#    oon <- on
#    on <- sum(gva)
#    orv <- rv
#    if(sum(is.na(gv))>0)
#      cat("coefs :",rc,"\nresiduals :",res,"\ngv :",gv,"\ngva :",gva,
#          "\nrv : ",rv," ca= ",ca,"\noon= ",oon," on= ",on," orv= ",orv,
#          "\nogv :",ogv,"\nrv1= ",rv1,"\nres1= ",res1,"\n")
    gva <- gv
#    print("Now reg:")
#    print(length(gv))
#    print(length(dep))
    if (p>0)
      reg <- lsfit(indep, dep, wt=gv)
    else
      reg <- lsfit(indep, dep, intercept=FALSE, wt=gv)
    rc <- coef(reg)
    res <- resid(reg)
    rv <- sum(res[as.logical(gv)]^2)/(sum(gv)-p-1)
    if (cachange){
      ca <- can(sum(gv),p,"lowresid")
#      cat("n= ",sum(gv)," p= ",p," ca= ",ca," change= ",change,"\n")
    }  
    gv <- res^2<=ca*rv
    coll <- FALSE
    if (p>0)
      coll <- (reg$qr$rank<ncol(as.matrix(indep))+1)
    if (coll)
      break
#    print("gv")
#    print(change)
    if (plot)
      plot(cbind(indep,dep)[,c(1,p+1)],col=1+gv,pch=1+gv)
  }
  out <- list(coef=rc, var=rv, g=gv, coll=coll, ca=ca)
  out
}

#
# fixreg main and output
#
fixreg <- function (indep=rep(1,n), dep, n=length(dep),
                    p=ncol(as.matrix(indep)),
                    ca=NA, mnc=NA, mtf=3, ir=NA, irnc=NA,
                    irprob=0.95, mncprob=0.5, maxir=20000, maxit=5*n,
                    distcut=0.85, init.group=list(), 
                    ind.storage=FALSE, countmode=100, 
                    plot=FALSE){

# Initializations, WDC

# print("fixreg")
  if (is.na(ca))
    ca <- can(n,p)
  if(is.na(ir) & is.na(irnc))
    irnc <- round(n/5)
  if(is.na(ir))
      ir <- itnumber(n,p,irnc,mtf,irprob,maxir)
  if(is.na(mnc) & ir>0)
      mnc <- minsize(n,p,ir,mtf,mncprob)
  if(is.na(mnc) & ir==0)
      mnc <- 1
  if (ir<mtf)
    mtf <- 1
  tsc <- 0          # number of too small FPCs
  ncoll <- 0        # number of iterations with collinear regressors
  gv <- rep (TRUE, times=n)
  fpc <- rfpi(indep, dep, p, gv, ca, maxit, plot)   # WDC-iteration
  if (fpc$coll){
    print("Independent variables of whole data collinear. Break forced.")
    break
  }
  imatrix <- c(sum(fpc$g)) # matrix of FPC size and intersections
  smatrix <- c(sum(fpc$g)) # matrix of FPC size and dissimilarities
  nc <- 1                  # number of FPCs
  clist <- list(fpc$coef)  # list of FPC regression coefficients
  vlist <- list(fpc$var)   # list of FPC residual variances
  nfound <- c(1)           # times found per FPC
  expectratio <- c()       # nfound/nfound expected
  if (ind.storage)
    glist <- list(fpc$g)   # list of FPC indicator vectors

# Standard iterations

  i <- 1
# print("standard it.")
  while(i<=ir){
    if(countmode*round(i/countmode)==i)
      cat("Iteration run no. ", i, " of ",ir, "\n")
    gv <- randconf(n,p+2)
# cat(length(gv)," ",sum(gv)," ")
# print("iteration")
    fpc <- rfpi(indep, dep, p, gv, ca, maxit, plot)
# print("neu?")
    neu <- TRUE
    if(!fpc$coll){
      j <- 1
      while(j<=nc) {
        if (identical(clist[[j]],fpc$coef) & identical(vlist[[j]],fpc$var)) {
          neu <- FALSE
          cnum <- j
        }          # if j found
        j <- j+1
      }            # for j
      ng <- sum(fpc$g)
      if (neu & (ng>=mnc)){
# print("neu! smatrix")
        nc <- nc +1
        nfound[nc] <- 1
        clist[[nc]] <- fpc$coef
        vlist[[nc]] <- fpc$var
# cat("ind.storage= ",ind.storage,  "nc= ", nc, "\n")
        if (ind.storage){
          glist[[nc]] <- fpc$g
          j <- 1
          while( j<=(nc-1)) {
            imatrix[sseg(nc,j)] <- sum(fpc$g & glist[[j]])
# cat("   (storage T) smatrix: ",nc,j,smatrix[sseg(nc,j)],"\n")
            j<- j+1
          }        # for j
        }          # if ind.storage
        else{
          for(j in 1:(nc-1)){
            glistj <- ((dep - cbind(1,indep) %*% clist[[j]])^2 <= ca*vlist[[j]])
            imatrix[sseg(nc,j)] <- sum(fpc$g & glistj)
# cat("   smatrix: ",nc,j,smatrix[sseg(nc,j)],"\n")
          }       # for j
        }          # else (!ind.storage)
        imatrix[sseg(nc,nc)] <- sum(fpc$g)
      }            # if neu & ng>=mnc
      if (ng<mnc)
        tsc <- tsc+1
      if (!neu)
        nfound[cnum] <- nfound[cnum]+1
    } # if !coll
    else
      ncoll <- ncoll+1
    i <- i+1
# print("schleifenende")
  }              # while i
        
# Iterations with init.group

#  print(length(init.group))
  grfpc <- FALSE
  if(length(init.group)>0){ 
    i <- 1
    grfpc <- rep(0, times=length(init.group))
    while(i<=length(init.group)){
      gv <- init.group[[i]]  
#      print(sum(gv))
#      print(length(gv))
      fpc <- rfpi(indep, dep, p, gv, ca, maxit, plot)
#      print("fpc")
      neu <- TRUE
      if(!fpc$coll){
        cnum <- nc +1
        for (j in 1:nc) {
          if (identical(clist[[j]],fpc$coef) & identical(vlist[[j]],fpc$var)) {
            neu <- FALSE
            cnum <- j
          }         # if j found
#        print(j)
        }           # for j
        ng <- sum(fpc$g)
        if (neu & (ng>=mnc)){
          nc <- nc +1
          nfound[nc] <- 1
          clist[[nc]] <- fpc$coef
          vlist[[nc]] <- fpc$var
          if (ind.storage){
            glist[[nc]] <- fpc$g
            for(j in 1:(nc-1))
              imatrix[sseg(nc,j)] <- sum(fpc$g & glist[[j]])
          }         # if ind.storage
          else{
            for(j in 1:(nc-1)){
              glistj <- ((dep - cbind(1,indep) %*% clist[[j]])^2 <= ca*vlist[[j]])
              imatrix[sseg(nc,j)] <- sum(fpc$g & glistj)
            }       # for j
          }         # else (!ind.storage)
        imatrix[sseg(nc,nc)] <- sum(fpc$g)
        }           # if neu & ng>=mnc
        if (ng<mnc)
          tsc <- tsc+1
        if(!neu)
          nfound[cnum] <- nfound[cnum]+1
      } # if !coll
      else{
        ncoll <- ncoll+1
        cnum <- 0
      } # else (coll)
      grfpc[i] <- cnum
      i <- i+1
    }             # for i
  }               # if init.group

# Find structures

# print(nc)
 if(nc>1){
    for (i in 1:(nc-1)){
      smatrix[sseg(i,i)] <- imatrix[sseg(i,i)]
      for (j in (i+1):nc){
        smatrix[sseg(i,j)] <- 2*imatrix[sseg(i,j)]/(imatrix[sseg(i,i)]+
                              imatrix[sseg(j,j)])
#        cat("smatrix: ",i,j,smatrix[sseg(i,j)],"\n")
      } # for j
    }
    smatrix[sseg(nc,nc)] <- imatrix[sseg(nc,nc)]
    comat <- matrix(nrow=nc,ncol=nc)
    for (i in 1:(nc-1))
      for(j in (i+1):nc)
        comat[i,j] <- comat[j,i] <- (smatrix[sseg(i,j)]>=distcut)
  } #if nc>1
  else
    comat <- matrix(1,1)
#  print(comat)
  struc <- con.comp(comat)
#  print(struc[1])
  stn <- max(struc)
  rm(comat)

# print("description of structures")

  tf <- rep(0, times=stn)   # structure: times found
  maxf <- rep(0, times=stn) # structure, rep. fpc: expectratio
  sfpc <- rep(0, times=stn) # structure: representative fpc
  for(i in 1:nc){
    expect <- clusexpect(n, p, imatrix[sseg(i,i)], ir)
    expectratio[i] <- nfound[i] / expect
    tf[struc[i]] <- tf[struc[i]] + nfound[i]
    if(expectratio[i]==maxf[struc[i]])
      if(imatrix[sseg(i,i)]<imatrix[sseg(sfpc[struc[i]],sfpc[struc[i]])])
        sfpc[struc[i]] <- i
    if(expectratio[i]>maxf[struc[i]]){
        sfpc[struc[i]] <- i
        maxf[struc[i]] <- expectratio[i]
    }
  } # for i
  stn <- sum(tf>=mtf)
#  print(stn)
  tsc <- tsc+sum(tf)-sum(tf[tf>=mtf])
  
# Output
  if (!ind.storage)
    glist <- FALSE
  out <- list(nc=nc, g=glist, coefs=clist, vars=vlist,
                nfound=nfound, er=expectratio, tsc=tsc, ncoll=ncoll,
                grto=grfpc, imatrix=imatrix,
                smatrix=smatrix, 
                stn=stn, stfound=tf,
                sfpc=sfpc, ssig=sfpc[tf>=mtf], sto=order(-tf), 
                struc=struc, n=n, p=p, ca=ca, ir=ir,
                mnc=mnc, mtf=mtf, distcut=distcut)  
  class(out) <- "rfpc"
# print ("fixreg ende")
  out   
}

summary.rfpc <- function(object, ...){
#  print("Beginn reducereg")
  clist <- list(0)
  vlist <- list(0)
  expectratio <- c()
  tf <- c()
  sn <- c()
  strn <- object$stn
#  print(strn)
  if (strn>0)
    for(i in 1:strn){
  #  cat(i, object$coefs[[object$sfpc[object$sto[i]]]], "\n")
      clist[[i]] <- object$coefs[[object$sfpc[object$sto[i]]]]
  #  cat(i, object$vars[[object$sfpc[object$sto[i]]]], "\n")
      vlist[[i]] <- object$vars[[object$sfpc[object$sto[i]]]]
  #  cat(i, object$stfound[object$sto[i]], "\n")
      tf[i] <- object$stfound[object$sto[i]]
      sn[i] <- object$imatrix[sseg(object$sfpc[object$sto[i]],object$sfpc[object$sto[i]])]
      expect <- clusexpect(object$n, object$p, sn[i], object$ir)
      expectratio[i] <- tf[i] / expect
    }
#  print(expectratio)
  sim <- simmatrix(object)
  out <- list(coefs=clist, vars=vlist, stn=strn, stfound=tf, sn=sn, 
              ser=expectratio, tsc=object$tsc, sim=sim,
              ca=object$ca, ir=object$ir, mnc=object$mnc, mtf=object$mtf)
  class(out) <- "summary.rfpc"
  out
}

fpclusters.rfpc <- function(object, indep=NA, dep=NA, ca=object$ca, ...){
  glist <- list()
  if (identical(indep,NA) & !identical(dep,NA))
    indep <- rep(1,n)
# cat("stn= ",object$stn,"\n")
  if(object$stn>0)
    for(i in 1:object$stn){
      if (object$g==FALSE){
        rc <- object$coefs[[ object$sfpc[object$sto[i]] ]]
        if (object$p==0)
          glist[[i]] <- (dep - rc)^2 <=
                         ca*object$vars[[ object$sfpc[object$sto[i]] ]]
        else
          glist[[i]] <- (dep - cbind(1, indep) %*% rc)^2 <=
                         ca*object$vars[[ object$sfpc[object$sto[i]] ]]
        if(sum(glist[[i]])<length(object$coefs[[1]])+1){
          cat("Warning! FPC ",i," too small, presumably because of rounding error\n")
          cat("To get the correct indicator vector, run fixreg again with 'ind.storage=T'\n") 
        } # if sum glist too small
      }   # if g==F
      else
        glist[[i]] <- object$g[[object$sfpc[object$sto[i]]]]
  # print(i)
    } # for i
  else
    warning("No FPCs were found often enough.")
  glist
}

# Visualization of the representative FPC no. no; x-axis: regression 
# linera combination => theoretical line= identity.
plot.rfpc <- function(x, indep=rep(1,n), dep, no, bw=TRUE,
                      main=c("Representative FPC No. ",no),
                      xlab="Linear combination of independents",
                      ylab=deparse(substitute(indep)),
                      xlim=NULL, ylim=range(dep), 
                      pch=NULL, col=NULL,...){
  n <- x$n
  p <- x$p
  sumobj <- summary(x)
  rc <- sumobj$coefs[[no]]
  lc <- cbind(1, indep) %*% rc
  ind <- ((dep-lc)^2 <= sumobj$ca*sumobj$vars[[no]])
  if (is.null(pch))
    pch <- ifelse(bw, 2, 1)
  if (is.null(col))
    col <- ifelse(bw, 1, 2)
  if (is.null(xlim))
    xlim <- range(lc)
  if (is.null(ylim))
    ylim <- range(dep)
  plot(lc[ind], dep[ind], main=main, xlab=xlab, ylab=ylab,
       xlim=xlim, ylim=ylim, pch=pch, col=col, ...)
  points(lc[!ind], dep[!ind], pch=1, col=1, ...)
  abline(c(0,1))
  abline(c(-sqrt(sumobj$ca*sumobj$vars[[no]]),1), lty="dotted")
  abline(c(sqrt(sumobj$ca*sumobj$vars[[no]]),1), lty="dotted")
  invisible()
}

print.rfpc <- function(x, ...){
  cat("Linear Regression Fixed Point Cluster object\n")  
  cat(x$stn," representative stable fixed point clusters\n") 
  cat(" of totally ",x$nc," found fixed point clusters.\n")
  invisible(x)
}

# Summary output of rfpc objects; choose fpcobj <- summary(fpcobj) to
# replace fpcobj by its reduction    
print.summary.rfpc <- function(x, maxnc=30, ...){
#  print("Beginn Summary")
  minnc <- min(x$stn,maxnc)
  cat("  *  Fixed Point Clusters  *\n\n")
  cat("Often a clear cluster in the data leads to several similar FPCs.\n")
  cat("The summary shows the representative FPCs of groups of similar FPCs,\n")
  cat("which were found at least ",x$mtf," times.\n\n")
  cat("Constant ca= ",x$ca,"\n")
  cat("Number of representative FPCs: ", x$stn, "\n\n")
  cat("FPCs with less than ",x$mnc," points were skipped.\n")
  cat(x$tsc," iterations led to skipped FPCs.\n\n")
  if (x$stn>maxnc)
    cat("Warning! Only ",maxnc," clusters are displayed. \nSpecify maxnc in call of summary.rfpc if you want enlarged display or \nrun fixreg with larger ca, mnc or mtf to get fewer clusters. \n\n")
  if (x$stn==0)
    cat("No FPCs were found often enough.\n")
  else{
    for(i in 1:minnc){
      cat(" FPC ",i, "\n")
      cat("  Times found (group members): ",x$stfound[i], "\n")
      if (x$ir>0)
        cat("  Ratio to estimated expectation: ",x$ser[i], "\n")
      cat("  Regression parameters:\n")
      print(x$coefs[[i]])
      cat("  Error variance: ", x$vars[[i]], "\n")
      cat("  Number of points: ",x$sn[[i]], "\n\n")
    }
    cat("Number of points in intersection of  representative FPCs\n")
    sm <- rep(0,times=minnc^2)
    dim(sm) <- c(minnc,minnc)
    for(i in 1:minnc)
      for(j in 1:minnc)
        sm[i,j] <- x$sim[sseg(i,j)]
    print.matrix(sm)
  }
  invisible(x)
}

.packageName <- "fpc"
.First.lib <- function(lib,pkg){
  require(cluster)
  if (as.numeric(R.version$major)<=1 & as.numeric(R.version$minor)<9){
    require(lqs)
    require(mva)
  }
  else require(MASS)
}
  

#
# generic functions
#
fpclusters <- function(object, ...)
  UseMethod("fpclusters")

#
# general/set distance utilities
#
#
# sseg= position of distance between i and j in smatrix/imatrix
sseg <- function (i,j){
  if (i>=j)
    out <- (i-1)*i/2+j
  else
    out <- (j-1)*j/2+i
  out
}

# Connectivity components of a graph by depth-first search
# Input: Boolean coincidence matrix, output: Vector of component numbers
con.comp <- function(comat){
  nc <- ncol(comat)
  ccn <- rep(0, times=nc) # con.comp number for each point
  fhist <- rep(FALSE, times=nc) # indicator if point had been under consideration  
  stn <- 0                  # current cc number
  pn <- 1                  # point no. to which similar objects are looked for
  while(pn>0){
    stn <- stn+1
    repeat{
      sm <- 0              # smallest new point no.
      ccn[pn] <- stn
      fhist[pn] <- TRUE
      if(nc>1)
      {
  	for(i in 2:nc)
        {
  #          cat(i, pn, ccn[i], comat[i,pn],"\n")
  	  if((ccn[i]==0) & (comat[i,pn]))
            ccn[i] <- stn
  	  if ((sm==0) & (ccn[i]==stn) & (fhist[i]==FALSE))
            sm <- i
  	} # for i
  #      cat("stn=", stn, "sm=", sm, "\n")
      } # if nc>1
      if (sm>0)
        pn <- sm
      else
        break
    } # repeat
#    print("repeat terminated")
    pn <- 0
    i <- 2
    while(i<=nc){
      if(ccn[i]==0){
        pn <- i
        i <- nc
      } # if
      i <- i+1
    } # while i
  } # while pn>0 (stn-loop)
  ccn
}    

# Intersection matrix between significant fpc groups
simmatrix <- function(fpcobj){
  sim <- c()
#  print(stn)
  for(i in 1:fpcobj$stn)
    for(j in i:fpcobj$stn){
#      cat(i,j,fpcobj$sfpc[fpcobj$sto[i]],fpcobj$sfpc[fpcobj$sto[j]])
      sim[sseg(i,j)] <- fpcobj$imatrix[sseg(fpcobj$sfpc[fpcobj$sto[i]],fpcobj$sfpc[fpcobj$sto[j]])]
    }
  sim
}

#
# fixmahal utilities
#
# Weighted ML-covariance matrix
cov.wml <- function (x, wt = rep(1/nrow(x), nrow(x)),
                     cor = FALSE, center = TRUE) 
{
    if (is.data.frame(x)) 
        x <- as.matrix(x)
    else if (!is.matrix(x)) 
        stop("x must be a matrix or a data frame")
    if (!all(is.finite(x))) 
        stop("x must contain finite values only")
    n <- nrow(x)
    if (with.wt <- !missing(wt)) {
        if (length(wt) != n) 
            stop("length of wt must equal the number of rows in x")
        if (any(wt < 0) || (s <- sum(wt)) == 0) 
            stop("weights must be non-negative and not all zero")
        wt <- wt/s
    }
    if (is.logical(center)) {
        center <- if (center) 
            colSums(wt * x)
        else 0
    }
    else {
        if (length(center) != ncol(x)) 
            stop("length of center must equal the number of columns in x")
    }
    x <- sqrt(wt) * sweep(x, 2, center)
    cov <- (t(x) %*% x)
    y <- list(cov = cov, center = center, n.obs = n)
    if (with.wt) 
        y$wt <- wt
    if (cor) {
        sdinv <- diag(1/sqrt(diag(cov)), nrow(cov))
        y$cor <- sdinv %*% cov %*% sdinv
    }
    y
}

# inversion of cov-matrices; if singular, Eigenvalues below 1/cmax are set
# to 1/cmax.
solvecov <- function(m, cmax=1e10){
  options(show.error.messages = FALSE)
  covinv <- try(solve(m))
  if(class(covinv)!="try-error")
     coll=FALSE 
  else{ 
    p <- nrow(m)
    cove <- eigen(m, symmetric=TRUE)
    coll <- TRUE
    if (min(cove$values)<1/cmax){
      covewi <- diag(p)
      for (i in 1:p)
        if (cove$values[i]<1/cmax) covewi[i,i] <- cmax
        else covewi[i,i] <- 1/cove$values[i]
    }
    else covewi <- diag(1/cove$values,nrow=length(cove$values))
#  cat("covewi ",covewi)
    covinv <- cove$vectors %*% covewi %*% t(cove$vectors)
  }
  options(show.error.messages = TRUE)
  out <- list(inv=covinv,coll=coll)
}

# Generation of ca by "recursion for convergence"
cmahal <- function(n,p,nmin,cmin,nc1,c1=cmin,q=1){
  cs <- rep(c1,max(n,nc1))
  for (i in (nc1-1):nmin)
    cs[i] <- cs[i+1]+q*(cs[i+1]-p)/(i+1-q)
  cs[1:(nmin-1)] <- rep(cs[nmin],nmin-1)
  if (n>nc1)
    for (i in (nc1+1):n){
      cs[i] <- cs[i-1]-q*(cs[i-1]-p)/i
      if (cs[i]<=cmin){
        cs[i] <- cmin
        i <- n
      }
    }
  cs
}

wfu <- function(md,ca,ca2,a1=1/(ca-ca2),a0=-a1*ca2){
  v <- rep(0,length(md))
  v[md<=ca] <- 1
  mind <- (md>ca) & (md<=ca2)
  v[mind] <- a1*md[mind]+a0
  v
}

# mahalanofix=vector of mahalanobis distances from all points to
# center of points indexed by gv (default=all)
mahalanofix <- function (x, n=nrow(as.matrix(x)), p=ncol(as.matrix(x)),
                         gv=rep(1, times=n), 
                         cmax=1e+10, method="ml") {
#  print("Mahalanofix")
  gv <- as.logical(gv)
  ng <- sum(gv)
  xg <- x[gv,1:p,drop=FALSE]
  if (method=="ml"){
    mg <- colMeans(xg)
    if (method=="ml")
      covg <- (ng-1)*cov(xg)/ng
    else
      covg <- cov(xg)
  }
  else{
    if (as.numeric(R.version$major)<=1 & as.numeric(R.version$minor)<9)
      require(lqs)
    else require(MASS)
    grob <- cov.rob(xg, method=method)
    mg <- grob$center
    covg <- grob$cov
  }
  cm <- solvecov(covg,cmax=cmax)
#  print("robmahal")
  md <- mahalanobis(x,mg,cm$inv,inverted=TRUE)
#  print("end robmahal")
  out <- list(md=md, mg=mg, covg=covg, covinv=cm$inv, coll=cm$coll)
}

# vector of mahalanobis distances from all points to
# center of points indexed by "fuzzy" gv (default=all)
mahalanofuz <- function (x, n=nrow(as.matrix(x)), p=ncol(as.matrix(x)),
                         gv=rep(1, times=n), 
                         cmax=1e+10) {
#  print("Mahalanofix")
  mg <- cov.wml(x,gv)
  covg <- mg$cov
  mg <- mg$center
  cm <- solvecov(covg, cmax=cmax)
#  print("Mahalanobis")
  md <- mahalanobis(x,mg,cm$inv,inverted=TRUE)
#  print("End Mahalanofix")
  out <- list(md=md, mg=mg, covg=covg, covinv=cm$inv, coll=cm$coll)
}

# Start configuration for fixmahal, point no. no
# startn points, cov/center of whole data
mahalconf <- function(x, no, startn, covall, plot){
  p <- ncol(x)
  n <- nrow(x)
  gv <- rep(FALSE, n)
  om <- order(mahalanobis(x,center=x[no,],cov=solvecov(covall)$inv, inverted=TRUE))
  gv[om[1:(p+1)]] <- TRUE
#  cat("conf ",gv," ",om,"\n")
  if (plot=="start" || plot=="both")
    plot(x,col=1+gv+((1:n)==no),pch=1+gv+((1:n)==no))
  if (startn>p+1)
    for (pn in (p+2):startn){
        om <- order(mahalanofix(x,n,p,gv)$md)
      j <- 0
      while (sum(gv)<pn){
        j <- j+1
        if (!gv[om[j]])
          gv[om[j]] <- TRUE
      }
    }
#   print(sum(gv))
  gv
}                

#
# fixmahal iteration
#
fpmi <- function (dat, n=nrow(as.matrix(dat)), p=ncol(as.matrix(dat)),
                  gv, ca, ca2, method="ml", plot,
                  maxit=5*n, iter=n*1e-6) {
  gva <- rep(0, times=n)
  change <- TRUE
  coll <- FALSE
  ir <- 0
  while (change) {
    ir <- ir+1
    if (method=="fuzzy")
      change <- sum((gv-gva)^2>iter)
    else
      change <- !identical(gv,gva)
    if (ir>maxit){
      change <- FALSE
      warning("Maximum number of iterations exceeded.")
    }
    gva <- gv
#  print("mahalanobis")
    if (method=="fuzzy")
      mfix <- mahalanofuz(x=dat, n=n, p=p, gv=gv)
    else
      mfix <- mahalanofix(x=dat, n=n, p=p, gv=gv, method=method)
#  print(ca)
    md <- mfix$md
    mg <- mfix$mg
    if (mfix$coll) coll=TRUE
    covg <- mfix$covg
#  print(md)
    if (method=="fuzzy"){
      gv <- wfu(md, ca, ca2)
#      print(gv)
      cn <- ca
#     plot(x,col=1+gv)
    }
    else{
      ng <- sum(gv)
      if (length(ca)==1)
        cn <- ca
      else
        cn <- ca[ng]
#    print(ng)
#    cat("Target function = ",
#      det(covg*(ng-1)/ng)*prod(exp((2-ca)/(1:ng))), "\n")
      gv <- md<=cn
    }
#  print(gv)
     if (plot=="iteration" || plot=="both")
       plot(dat,col=1+gv,pch=1+gv)
  } # while change
  out <- list(mg=mg, covg=covg, md=md, gv=gv, coll=coll, method=method, ca=cn)
  out
}

#
# fixmahal main and output
#
fixmahal <- function (dat, n=nrow(as.matrix(dat)), p=ncol(as.matrix(dat)), 
                      method="fuzzy", cgen="fixed", 
                      ca=NA, ca2=NA,
                      calpha=ifelse(method=="fuzzy",0.9,0.99),
                      calpha2=0.99,
                      pointit=TRUE, subset=n,
                      mnc=min(floor(n/2),10+2*p), nc1=100+20*p,
                      startn=min(floor(n/2),ifelse(cgen=="auto",mnc,2*mnc)),
                      mer=ifelse(pointit,0.1,0), distcut=0.85,  
                      maxit=5*n, iter=n*1e-5, 
                      init.group=list(), 
                      ind.storage=TRUE, countmode=100, 
                      plot="none"){

# ca2, calpha2 are ignored unless method="fuzzy"
# Initializations, WDC

# print("fixmahal")
  dat <- as.matrix(dat)
  if (method=="fuzzy"){
    if (is.na(ca)){
      ca <- qchisq(calpha,p)
      ca2 <- qchisq(calpha2,p)
    }
  }
  else{
    if (cgen=="fixed"){
      if (is.na(ca))
        ca <- qchisq(calpha,p)
      else
        calpha <- pchisq(ca,p)
    }
    else{
      if (is.na(ca))
        ca <- qchisq(calpha,p)
      ca <- cmahal(n,p,nmin=startn,cmin=ca,nc1=nc1,q=1)
    }
    ca2 <- 0
  }  
  dat <- dat[1:n,1:p,drop=FALSE]
  tsc <- 0          # number of too small FPCs
  ncoll <- 0        # number of iterations leading to collinear data
  gv <- rep (1, times=n)
# print("fpmi")
  fpc <- fpmi(dat, n, p, gv, ca, ca2, method, plot, maxit, iter)   # WDC-iteration
# print("1Ende")
  if (cgen=="auto")
    cv <- c(fpc$ca)
  imatrix <- c(sum(fpc$g)) # matrix of FPC size and intersections
  smatrix <- c(sum(fpc$g)) # matrix of FPC size and similarities
  nc <- 1                  # number of FPCs
  clist <- list(fpc$mg)  # list of FPC locations
  vlist <- list(fpc$covg)   # list of FPC covariance matrices
  nfound <- c(1)           # times found per FPC
  expectratio <- c()       # nfound/size
  if (ind.storage)
    glist <- list(fpc$g)   # list of FPC indicator vectors

# Standard iterations

  if (pointit){
    if(subset < n)
      itpoints <- sample(1:n,subset)
    else{
      itpoints <- 1:n
      subset <- n
    }
    if (method=="fuzzy" || method=="classical" || method=="ml")
      covall <- cov(dat)
    else
      covall <- cov.rob(dat, method=method)$cov
    for (ni in 1:subset){
      i <- itpoints[ni]
      if(countmode*round(ni/countmode)==i)
        cat("Iteration run no. ", ni, " of ",subset, "\n")
# print(covall)
      gv <- mahalconf(dat,i,startn,covall=covall, plot)
#  cat("i= ",i," gv=  ",(1:366)[gv]," \n")
# print("iteration")
  #      cat(sum(gv[1:50]),sum(gv[51:100]))
      fpc <- fpmi(dat, n, p, gv, ca, ca2, method, plot, maxit, iter)
#   print(fpc$mg)
  # print("neu?")
      neu <- TRUE
      j <- 1
      while(j<=nc) {
        if (sum((clist[[j]]-fpc$mg)^2)<1e-6 && sum((vlist[[j]]-fpc$covg)^2)<1e-5){
          neu <- FALSE
          cnum <- j
        }          # if j found
        j <- j+1
      }            # while j
      ng <- sum(fpc$g)
      if (neu & (ng>=mnc)){
        nc <- nc +1
        nfound[nc] <- 1
        clist[[nc]] <- fpc$mg
        vlist[[nc]] <- fpc$covg
        if (cgen=="auto")
          cv <- c(cv,fpc$ca)
# cat("ind.storage= ",ind.storage,  "nc= ", nc, "\n")
        if (ind.storage){
          glist[[nc]] <- fpc$g
          j <- 1
          while( j<=(nc-1)) {
            imatrix[sseg(nc,j)] <- sum(pmin(fpc$g,glist[[j]]))
# cat("   (storage T) smatrix: ",nc,j,smatrix[sseg(nc,j)],"\n")
            j<- j+1
          }        # for j
        }          # if ind.storage
        else{
          for(j in 1:(nc-1)){
            if (method=="fuzzy"){
              mah <- mahalanobis(dat,center=clist[[j]],
                           cov=solvecov(vlist[[j]])$inv,inverted=TRUE)
              glistj <- wfu(mah, ca, ca2)
            }
            else{
              if (cgen=="fixed"){
                  glistj <- (mahalanobis(dat,center=clist[[j]],
                           cov=solvecov(vlist[[j]])$inv,inverted=TRUE) <= ca)
              }
              else{
                  glistj <- (mahalanobis(dat,center=clist[[j]],
                           cov=solvecov(vlist[[j]])$inv,inverted=TRUE) <= cv[j])
              }
            }
            imatrix[sseg(nc,j)] <- sum(pmin(fpc$g,glistj))
# cat("   smatrix: ",nc,j,smatrix[sseg(nc,j)],"\n")
          }       # for j
        }          # else (!ind.storage)
        imatrix[sseg(nc,nc)] <- sum(fpc$g)
#        cat("Neu: Punkt ",i,", Cluster mit ",sum(fpc$g)," Punkten. \n")
      }            # if neu & ng>=mnc
      if (ng<mnc)
        tsc <- tsc+1
      if (!neu){
        nfound[cnum] <- nfound[cnum]+1
      }
      if (neu & fpc$coll)
        ncoll <- ncoll+1
    }              # for i
  }                # if pointit
# Iterations with init.group

#  print(length(init.group))
  grfpc <- FALSE
  if(length(init.group)>0){ 
    i <- 1
    grfpc <- rep(0, times=length(init.group))
    while(i<=length(init.group)){
      gv <- init.group[[i]]  
#      print(sum(gv))
#      print(length(gv))
      fpc <- fpmi(dat, n, p, gv, ca, ca2, method, plot, maxit, iter)
#      print("fpc")
      neu <- TRUE
      j <- 1
      while(j<=nc) {
        if (sum((clist[[j]]-fpc$mg)^2)<1e-6 && sum((vlist[[j]]-fpc$covg)^2)<1e-5){
          neu <- FALSE
          cnum <- j
        }          # if j found
        j <- j+1
      }            # while j
      ng <- sum(fpc$g)
      if (neu & (ng>=mnc)){
# print("neu! smatrix")
        nc <- nc +1
        nfound[nc] <- 1
        clist[[nc]] <- fpc$mg
        vlist[[nc]] <- fpc$covg
        if (cgen=="auto")
          cv <- c(cv,fpc$ca)
# cat("ind.storage= ",ind.storage,  "nc= ", nc, "\n")
        if (ind.storage){
          glist[[nc]] <- fpc$g
          j <- 1
          while( j<=(nc-1)) {
            imatrix[sseg(nc,j)] <- sum(pmin(fpc$g,glist[[j]]))
# cat("   (storage T) smatrix: ",nc,j,smatrix[sseg(nc,j)],"\n")
            j<- j+1
          }        # for j
        }          # if ind.storage
        else{
          for(j in 1:(nc-1)){
            if (method=="fuzzy"){
              mah <- mahalanobis(dat,center=clist[[j]],
                           cov=solvecov(vlist[[j]])$inv,inverted=TRUE)
              glistj <- wfu(mah, ca, ca2)
            }
            else{
              if (cgen=="fixed"){
                  glistj <- (mahalanobis(dat,center=clist[[j]],
                           cov=solvecov(vlist[[j]])$inv,inverted=TRUE) <= ca)
              }
              else{
                  glistj <- (mahalanobis(dat,center=clist[[j]],
                           cov=solvecov(vlist[[j]])$inv,inverted=TRUE) <= cv[j])
              }
            }
            imatrix[sseg(nc,j)] <- sum(pmin(fpc$g,glistj))
# cat("   smatrix: ",nc,j,smatrix[sseg(nc,j)],"\n")
          }       # for j
        }          # else (!ind.storage)
        imatrix[sseg(nc,nc)] <- sum(fpc$g)
      }            # if neu & ng>=mnc
      if (ng<mnc)
        tsc <- tsc+1
      if (!neu)
        nfound[cnum] <- nfound[cnum]+1
      else
        cnum <- nc
      if (neu & fpc$coll)
        ncoll <- ncoll+1
      grfpc[i] <- cnum
      i <- i+1
    }             # while i
  }               # if init.group

# Find structures

# print(nc)
 if(nc>1){
    for (i in 1:(nc-1)){
      smatrix[sseg(i,i)] <- imatrix[sseg(i,i)]
      for (j in (i+1):nc){
        smatrix[sseg(i,j)] <- 2*imatrix[sseg(i,j)]/(imatrix[sseg(i,i)]+
                              imatrix[sseg(j,j)])
#        cat("smatrix: ",i,j,smatrix[sseg(i,j)],"\n")
      } # for j
    }
    smatrix[sseg(nc,nc)] <- imatrix[sseg(nc,nc)]
    comat <- matrix(nrow=nc,ncol=nc)
    for (i in 1:(nc-1))
      for(j in (i+1):nc)
        comat[i,j] <- comat[j,i] <- (smatrix[sseg(i,j)]>=distcut)
  } #if nc>1
  else
    comat <- matrix(1,1)
#  print(comat)
  struc <- con.comp(comat)
#  print(struc)
  stn <- max(struc)
  rm(comat)

# print("description of structures")

  tf <- rep(0, times=stn)   # structure: times found
  maxf <- rep(0, times=stn) # structure, rep. fpc: expectratio
  sfpc <- rep(0, times=stn) # structure: representative fpc
  ser <- rep(0, times=stn) # structure: expectation ratio
  for(i in 1:nc){
    expectratio[i] <- nfound[i] / imatrix[sseg(i,i)]
    tf[struc[i]] <- tf[struc[i]] + nfound[i]
    if(expectratio[i]==maxf[struc[i]])
      if(imatrix[sseg(i,i)]<imatrix[sseg(sfpc[struc[i]],sfpc[struc[i]])])
        sfpc[struc[i]] <- i
    if(expectratio[i]>maxf[struc[i]]){
        sfpc[struc[i]] <- i
        maxf[struc[i]] <- expectratio[i]
    }
  } # for i
  for (i in 1:stn)
    ser[i] <- (tf[i]*n)/(imatrix[sseg(sfpc[i],sfpc[i])]*subset)
  skc <- stn
  stn <- sum(ser>=mer)
  skc <- skc-stn
#  print(stn)
  
# Output
  cvec <- ca
  if (cgen=="auto")
    ca <- cv
  if (!ind.storage)
    glist=FALSE
  out <- list(nc=nc, g=glist, means=clist, covs=vlist,
                nfound=nfound, er=expectratio,
                tsc=tsc, ncoll=ncoll, skc=skc, grto=grfpc,
                imatrix=imatrix, smatrix=smatrix,
                stn=stn, stfound=tf, ser=ser, 
                sfpc=sfpc, ssig=sfpc[ser>=mer], sto=order(-ser), 
                struc=struc, 
                n=n, p=p, method=method, cgen=cgen, 
                ca=ca, ca2=ca2, cvec=cvec, calpha=calpha,
                pointit=pointit, subset=subset, 
                mnc=mnc, startn=startn, mer=mer, distcut=distcut)
  class(out) <- "mfpc"
  out   
}

summary.mfpc <- function(object, ...){
#  print("Beginn reducemahal")
  clist <- list(0)
  vlist <- list(0)
  expectratio <- c()
  tf <- c()
  sn <- c()
  ca <- object$ca
  method <- object$method
  strn <- object$stn
  if (strn>0)
    for(i in 1:strn){
#    cat(i, object$coefs[[object$sfpc[object$sto[i]]]], "\n")
      clist[[i]] <- object$means[[object$sfpc[object$sto[i]]]]
#    cat(i, object$vars[[object$sfpc[object$sto[i]]]], "\n")
      vlist[[i]] <- object$covs[[object$sfpc[object$sto[i]]]]
#    cat(i, object$stfound[object$sto[i]], "\n")
      tf[i] <- object$stfound[object$sto[i]]
      sn[i] <- object$imatrix[sseg(object$sfpc[object$sto[i]],object$sfpc[object$sto[i]])]
      expectratio[i] <- object$ser[object$sto[i]]
      if (object$cgen=="auto")
        ca[i] <- object$ca[object$sfpc[object$sto[i]]]
    }
#    print(expectratio)
  tskip <- sum(object$nfound)-sum(tf)
  sim <- simmatrix(object)
  out <- list(means=clist, covs=vlist, stn=strn, stfound=tf, sn=sn, 
              ser=expectratio, 
              tskip=tskip, skc=object$skc, tsc=object$tsc, sim=sim, 
              ca=ca, ca2=object$ca2, calpha=object$calpha,
              mer=object$mer, mnc=object$mnc,
              method=method, cgen=object$cgen,
              pointit=object$pointit)
  class(out) <- "summary.mfpc"
  out
}

# clusters: List of indicator vectors of representative 
# fpcs of significant structures
# in order of stfound.
# Specify x if ind.storage in fixreg was F
fpclusters.mfpc <- function(object, dat=NA, ca=object$ca, p=object$p, ...){
  glist <- list()
  if (object$method=="fuzzy"){
    ca <- qchisq(0.9,object$p)
    ca2 <- 10*qchisq(0.9,object$p)
  }
# cat("stn= ",object$stn,"\n")
  if (object$stn>0)
    for(i in 1:object$stn){
      if (identical(object$g,FALSE)){
        mc <- object$means[[ object$sfpc[object$sto[i]] ]]
        if (p>1){
          cm <- solvecov(object$covs[[ object$sfpc[object$sto[i]] ]])
          if (object$method=="fuzzy")
            glist[[i]] <- wfu(mahalanobis(dat,mc,cm$inv,inverted=TRUE), ca, ca2)
          if (object$cgen=="fixed")
            glist[[i]] <- (mahalanobis(dat,mc,cm$inv,inverted=TRUE)<=ca)
          else
            glist[[i]] <- (mahalanobis(dat,mc,cm$inv,inverted=TRUE)<=ca[object$sfpc[object$sto[i]]] )
        }
        else{
          cm <- object$covs[[ object$sfpc[object$sto[i]] ]]
          if (object$cgen=="fixed")
            glist[[i]] <- as.vector((dat-mc)^2)<=as.real(ca*cm)
          else
            glist[[i]] <- as.vector((dat-mc)^2)<=as.real(ca[object$sfpc[object$sto[i]]]*cm)
        }
      }   # if g==F
      else
        glist[[i]] <- object$g[[object$sfpc[object$sto[i]]]]
  # print(i)
    } # for i
  else
    warning("No FPC was found often enough!")
  glist
}

# Bhattacharyya coordinates of FPC no. no
# bw: black/white
plot.mfpc <- function(x, dat, no, bw=FALSE,
                      main=c("Representative FPC No. ",no),
                      xlab=NULL, ylab=NULL,
                      pch=NULL, col=NULL, ...){
  n <- x$n
  p <- x$p
  ca <- x$ca
  sumobj <- summary(x)
  if (x$cgen=="auto")
    ca <- sumobj$ca[no]
  mg <- sumobj$means[[no]]
  covg <- sumobj$covs[[no]]
#  print("clusters")
  if (x$method=="fuzzy"){
    g <- fpclusters(x,dat)[[no]] 
    gv <- as.integer(g>0.5)
    gvp <- as.integer(g>0.99)+
      3*as.integer(g>0.01 & g<=0.99)
  }
  else
    gvp <- gv <- as.integer(fpclusters(x,dat)[[no]])
  if (n==sum(gv)){
    cat("Cluster equals whole dataset => No cluster plot is done.\n")
  }
  else{
    if (is.null(pch))
      pch <- if (bw) gvp+1 else 1
    if (is.null(col))
      col <- if (bw) 1 else gvp+1
#  print("plot")
    if (p==1){
        if (is.null(xlab))
          xlab <- "Index"
        if (is.null(ylab))
          ylab <- deparse(substitute(dat))
        plot(dat,col=col, pch=pch, main=main, xlab=xlab, ylab=ylab, ...)
        abline(c(mg,0))
        abline(c(mg-sqrt(ca*covg),0), lty="dotted")       
        abline(c(mg+sqrt(ca*covg),0), lty="dotted")
    }
    if (p==2){
        xy <- xy.coords(dat,NULL)
        if (is.null(xlab))
          xlab <- xy$xlab
        if (is.null(ylab))
          ylab <- xy$ylab
        plot(dat,col=col, pch=pch, main=main, xlab=xlab, ylab=ylab, ...)
        points(as.vector(mg[1]),as.vector(mg[2]),pch="M",
                 col=ifelse(bw, col, "blue"), ...)
        cge <- eigen(covg, symmetric=TRUE)
        poly <- rep(0, times=200)
        dim(poly) <- c(100,2)
        for(i in 1:100)
          poly[i,] <- (mg + c(sin(i*2*pi/100), cos(i*2*pi/100)) %*% 
            (t(cge$vectors) * sqrt(cge$values))* sqrt(ca))
        polygon(poly[,1], poly[,2], border="black")  
    }
    if (p>2){
      gc1 <- batcoord(dat,gv)
        if (is.null(xlab))
          xlab <- "Discriminant projection 1"
        if (is.null(ylab))
          ylab <- "Discriminant projection 2"
        plot(gc1$proj[,1:2],col=col, pch=pch,
             main=main, xlab=xlab, ylab=ylab, ...)
        points(as.vector((mg %*% gc1$units)[1]),
               as.vector((mg %*% gc1$units)[2]), pch="M",
                 col=ifelse(bw, col, "blue"), ...)
    }
  }
  invisible()
}

print.mfpc <- function(x, ...){
  cat("Mahalanobis Fixed Point Cluster object\n")  
  cat(x$stn," representative stable fixed point clusters\n") 
  cat(" of totally ",x$nc," found fixed point clusters.\n")
  invisible(x)
}

# Summary output of mfpc objects; choose fpcobj <- summary(fpcobj) to
# replace fpcobj by its reduction    
print.summary.mfpc <- function(x, maxnc=30, ...){
#  print("Beginn Summary")
#  print(class(x))
  minnc <- min(x$stn,maxnc)
  cat("  *  Mahalanobis Fixed Point Clusters  *\n\n")
  cat("Often a clear cluster in the data leads to several similar FPCs.\n")
  cat("The summary shows the representative FPCs of groups of similar FPCs.\n\n")
  cat("Method ",x$method," was used.\n")
  cat("Number of representative FPCs: ", x$stn, "\n\n")
  cat("FPCs with less than ",x$mnc," points were skipped.\n")
  if (x$mer>0) 
    cat("FPCs with ratio of times found to number of points less than ",
        x$mer," were skipped.\n")
  cat(x$tsc+x$tskip, " iteration runs led to ",x$skc+x$tsc," skipped clusters.\n")
  if (x$stn>maxnc)
    cat("Warning! Only ",maxnc," clusters are displayed. \nSpecify maxnc in call of summary.mfpc if you want enlarged display or \nrun fixmahal with larger ca, calpha, mnc or mer to get fewer clusters. \n")
  if (x$method=="fuzzy")
    cat("  Weight 1 for r^2<= ",x$ca," weight 0 for r^2> ",
    x$ca2," \n")   
  if (x$cgen=="fixed")
    cat("  Constant ca= ", x$ca,
      " corresponding to alpha= ",x$calpha,"\n")
  cat("\n")
  if (x$stn==0)
    cat("No FPC was found often enough.\n")
  else{
    for(i in 1:minnc){
      cat(" FPC ",i, "\n")
      cat("  Times found (group members): ",x$stfound[i], "\n")
      if (x$pointit)
        cat("  Ratio to size: ",x$ser[i], "\n")
      cat("  Mean:\n")
      print(x$means[[i]])
      cat("  Covariance matrix:\n")
      print(x$covs[[i]])
      if (x$cgen=="auto" && x$method!="fuzzy")
        cat("  Constant ca= ",x$ca[[i]],"\n")
      cat("  Number of points (sum of weights): ",x$sn[[i]], "\n\n")
    }
    cat("Number of points (rounded weights) in intersection of representative FPCs\n")
    sm <- rep(0,times=minnc^2)
    dim(sm) <- c(minnc,minnc)
    for(i in 1:minnc)
      for(j in 1:minnc)
        sm[i,j] <- round(x$sim[sseg(i,j)])
    print.matrix(sm)
  }
  invisible(x)
}











# d is a distance object or matrix, clustering is assumed to be numbered
# 1 to clusternumber
# If alt.clustering is another clustering, the corrected rand index will
# be computed.
# silhouette, G2 and G3 indicate if the corresponding statistics are
# computed. Silhouette requires library cluster, G2 and G3 may take very long
# for large n.
cluster.stats <- function(d,clustering,alt.clustering=NULL,
                          silhouette=TRUE,G2=FALSE,G3=FALSE){
  if (silhouette) require(cluster)
  cn <- max(clustering)
  n <- length(clustering)
  dmat <- as.matrix(d)
  diameter <- average.distance <- median.distance <- separation <-
    average.toother <- 
    cluster.size <- within.dist <- between.dist <- numeric(0)
  separation.matrix <- matrix(0,ncol=cn,nrow=cn)
  di <- list()
  for (i in 1:cn){
    cluster.size[i] <- sum(clustering==i)
    di <- as.dist(dmat[clustering==i,clustering==i])
    within.dist <- c(within.dist,di)
    diameter[i] <- max(di)
    average.distance[i] <- mean(di)
    median.distance[i] <- median(di)
    bv <- numeric(0)
    for (j in 1:cn){
      if (j!=i){
        sij <- dmat[clustering==i,clustering==j]
        bv <- c(bv,sij)
        if (i<j){
          separation.matrix[i,j] <- separation.matrix[j,i] <- min(sij)
          between.dist <- c(between.dist,sij)
        }
      }
    }
    separation[i] <- min(bv)
    average.toother[i] <- mean(bv)
  }
  average.between <- mean(between.dist)
  average.within <- mean(within.dist)
  nwithin <- length(within.dist)
  nbetween <- length(between.dist)
  clus.avg.widths <- avg.width <- NULL
  if (silhouette){
    sc <- summary(silhouette(clustering,dmatrix=dmat))
    clus.avg.widths <- sc$clus.avg.widths
    avg.width <- sc$avg.width
  }
  g2 <- g3 <- corrected.rand <- cn2 <- NULL
  if (G2){
    splus <- sminus <- 0
    for (i in 1:nwithin)
      for (j in 1:nbetween){
        if (within.dist[i]<between.dist[j]) splus <- splus+1
        if (within.dist[i]>between.dist[j]) sminus <- sminus+1
      }
    g2 <- (splus-sminus)/(splus+sminus)
  }
  if (G3){
    sdist <- sort(c(within.dist,between.dist))
    sr <- nwithin+nbetween
    dmin <- sum(sdist[1:nwithin])
    dmax <- sum(sdist[(sr-nwithin+1):sr])
    g3 <- (sum(within.dist)-dmin)/(dmax-dmin)
  }
  if (!is.null(alt.clustering)){
    choose2 <- function(v){
      out <- numeric(0)
      for (i in 1:length(v))
        out[i] <- ifelse(v[i]>=2,choose(v[i],2),0)
      out
    }
    cn2 <- max(alt.clustering)
    nij <- table(clustering,alt.clustering)
    dsum <- sum(choose2(nij))
    cs2 <- numeric(0)
    for (i in 1:cn2)
      cs2[i] <- sum(alt.clustering==i)
    sum1 <- sum(choose2(cluster.size))
    sum2 <- sum(choose2(cs2))
    corrected.rand <- (dsum-sum1*sum2/choose2(n))/
      ((sum1+sum2)/2-sum1*sum2/choose2(n))
  }
  hubertgamma <- cor(c(within.dist,between.dist),c(rep(0,nwithin),
                                                   rep(1,nbetween)))
  dunn <- min(as.dist(separation))/max(diameter)
  out <- list(n=n,
              cluster.number=cn,
              cluster.size=cluster.size, # vector of cluster sizes
              diameter=diameter, # vector of cluster diameters
              average.distance=average.distance,
                                        # vector of within cl. av. dist.
              median.distance=median.distance,
                                        # vector of within cl. median dist.
              separation=separation, # vector of min. clusterwise between dist.
              average.toother=average.toother, 
                                     # vector of mean clusterwise between dist.
              separation.matrix=separation.matrix,
                                     # clusterwise matrix of min. between dist.
              average.between=average.between, # mean between cl. distance
              average.within=average.within, # mean within cl. distance
              n.between=nbetween, # number of between cl. distances
              n.within=nwithin, # number of within cl. distances
              clus.avg.silwidths=clus.avg.widths,
                                # vector of cluster avg. silhouette widths
              avg.silwidth=avg.width, # average silhouette width
              g2=g2, # Goodman and Kruskal coefficient, see Gordon p. 62
              g3=g3, # G3 index, see Gordon p. 62
              hubertgamma=hubertgamma, # Correlation between distances and
                                        # 0-1-vector same/different cluster
              dunn=dunn, # Dunn index, see Halkidi et al. (2002)
                         # Min. sepatation / max. diameter
              wb.ratio=average.within/average.between,
              corrected.rand=corrected.rand) # Corrected rand index between
                                        # clustering and alt.clustering
#  class(out) <- "cluster.stats"
  out
}


#
# face benchmark dataset by M. Maechler and C. Hennig
#
## MM: -  function(n, p)  {where `n' was a bit tricky}
##     -  return grouping as well

rFace <- function(n, p = 6, nrep.top = 2, smile.coef = 0.6,
                  dMoNo = 1.2, dNoEy = 1)
{
    ## Purpose: Generate random "Face" data set -- to be "Hard for Clustering"
    ## -------------------------------------------------------------------------
    ## Arguments: (n,p)    : dimension of result
    ##            nrep.top : #{repetitions} of the top point / "hair tip"
    ##            dMoNo    : distance{Mouth, Nose}
    ##            dNoEy    : distance{Nose,  Eyes} {vertically only}
    ##
    ## MM- TODO's:   o  provide more arguments (face shape)
    ##			    --> separation of clusters
    ## -------------------------------------------------------------------------
    ## Author: Christian Hennig & Martin Maechler, 26 Jun 2002

    if((p <- as.integer(p)) < 2) stop("number of variables p must be at least 2")
    if((n <- as.integer(n)) < 10) stop("number of points  n  must be at least 10")
    if((nrep.top <- as.integer(nrep.top)) < 1) stop("`nrep.top' must be positive")

    ntips <- nrep.top + 2
    n0 <- n - ntips
    m <- n0 %/% 5
    n.5 <- n0 %/% 2
    n.7 <- n.5+m
    n.9 <- n.7+m
    ## shouldn't happen:
    if(m < 1 || n.9 >= n0) stop("number of points n is too small")

    ## Indices of the different groups :
    Gr <- list(chin = 1:m,
               mouth= (m+1):n.5,
               nose = (n.5+1):n.7,
               rEye = (n.7+1):n.9,
               lEye = (n.9+1):n0,
               tips = (n0+1):n)

    face <- matrix(nrow = n, ncol = p)
    ## chin :
    face[Gr$chin, 1] <- U <- runif(m, -3, 3)
    face[Gr$chin, 2] <- rnorm(m, mean = U^2, sd=0.1)
    ## mouth:
    m0m <- 3 # lower mouth mean
    face[Gr$mouth, 1] <- Z <- rnorm(n.5-m, sd= 0.5)
    face[Gr$mouth, 2] <- rnorm(n.5-m, mean = m0m + smile.coef * Z^2, sd= 0.2)
    ## nose :
    face[Gr$nose, 1] <-  rnorm(m, mean=0, sd=0.2)
    yEye <- 17
    ##face[Gr$nose, 2] <- rnorm(m, mean=9, sd= 2.5)
    face[Gr$nose, 2] <- pmin(yEye - dNoEy,
                             (m0m + dMoNo) + rgamma(m, shape= 0.8, scale= 2.5))
    ## right eye  U[circle] :
    rangle <- runif(m, 0, 2*pi)
    rpos   <- runif(m)
    face[Gr$rEye, 1] <- rpos*cos(rangle) +  2
    face[Gr$rEye, 2] <- rpos*sin(rangle) + yEye
    ## left eye:
    face[Gr$lEye, 1] <- rnorm(n0-n.9, mean= -2,   sd= 0.5)
    face[Gr$lEye, 2] <- rnorm(n0-n.9, mean= yEye, sd= 0.5)
    ## `ntips'  ``hair tips'', the last one `nrep.top' times:
    face[n0+1,1:2] <- c(-4.5, 25)
    face[n0+2,1:2] <- c( 4.5, 25)
    for(k in 1:nrep.top)
        face[n-k+1, 1:2] <- c(0,32)

    ##-- Extra coordinates with noise ---

    if(p >= 3) {
        face[, p] <- rexp(n)
        if(p >= 4) {
            face[, p-1] <- rt(n,df=1)
            if(p >= 5)
                for(k in 3:(p-2))
                    face[,k] <- rnorm(n)
        }
    }
    gr <- character(n)
    ng <- names(Gr)
    for(i in seq(ng)) gr[Gr[[i]]] <- ng[i]
    structure(face, grouping = as.factor(gr), indexlist = Gr)
}


# randcmatrix=random partition matrix for n observations to cln clusters
randcmatrix <- function (n,cln,p){
  ct <- 0
  while(ct<p+3){
    m <- rep(0, times=n*cln)
    summ <- rep(0, times=cln)
    dim(m) <- c(n,cln)
    for (i in 1:n){
      nummer <- round(0.5+cln*runif(1))
#      print(m[i])
      m[i,nummer] <- 1
      summ[nummer] <- summ[nummer] + 1
    }
    ct <- min(summ)
  }
  m
}


# Regression EM iteration;
# m = random partition matrix, cln = number of clusters, icrit = iteration
# stopping criterion, minsig = minimal error variance
# output: coef=regression coefficients, var=residual variances, eps=cluster 
# proportions, z=posterior probabilities, loglik= loglikelihood, warn= T
# if too small or collinear cluster 
regem <- function (indep, dep, m, cln, icrit=1.e-5, minsig=1.e-6,
                         warnings=FALSE) {
  n <- length(dep)
  p <- ncol(as.matrix(indep))
  loglik <- (-1.e8)
  eps <- rep(0,cln)
  fv <- rep(0,n*cln)
  dim(fv) <- c(n,cln)
  rc <- rep(0,(p+1)*cln)
  dim(rc) <- c(p+1,cln)
  rv <- rep(0,cln)
  stm <- rep(0,n)  
  change <- TRUE
  smallcluster <- FALSE
  while (change & !smallcluster) {
#    plot(indep,dep)
    for(i in 1:cln){
      eps[i] <- sum(m[,i])/n
      if (sum(m[,i]>0.01) < p+2){
        if (warnings) warning("Too small cluster")
        smallcluster <- TRUE
      } # if too small cluster
      else{
        reg <- lm(dep~indep, weights=m[,i])
        fv[,i] <- fitted.values(reg)
        rc[,i] <- coefficients(reg)
#        abline(rc[,i],col=i)
        for (j in 2:(p+1))
          if (is.na(rc[j,i])){
            smallcluster <- TRUE
            if (warnings) warning("Collinear regressors")
          } # if collinearity
        res <- residuals(reg)
        rv[i] <- weighted.mean(res^2,m[,i])
        if (rv[i]<minsig){
          rv[i] <- minsig
          if (warnings) warning("Error variance smaller than minimum.")
        } # if error variance below minimum
      } # else (cluster large enough)
    } # for i
    if(!smallcluster){
        for (i in 1:cln)
          for (j in 1:n){
# cat("i= ",i," j= ",j, "dep[j]= ",dep[j]," mean= ",fv[j,i]," sd= ",sqrt(rv[i]),    "\n")
            m[j,i] <- eps[i]*dnorm(dep[j],mean=fv[j,i], sd=sqrt(rv[i]))
          }
        for (j in 1:n){
          stm[j] <- sum(m[j,])
          for (i in 1:cln)
            m[j,i] <- m[j,i]/stm[j]
        } # for j
        oldlog <- loglik        
        loglik <- sum(log(stm))
        change <- (loglik - oldlog > icrit)
    } # if no collinearity & clusters large enough
  } # while change
  g <- c()
  for (i in 1:n)
    g[i] <- which.max(m[i,])
  out <- list(coef=rc, vars=rv, z=m, g=g, eps=eps, loglik=loglik,
              warn=smallcluster)
  out
} # regem     



# Regression mixture analysis (DeSarbo and Cron), 
# ir=iteration runs, nclust= cluster numbers vector, icrit=iteration stopping 
# criterion, minsig = minimum error variance
regmix <- function (indep, dep,
                    ir=1, nclust=1:7, icrit=1.e-5, minsig=1.e-6,
                    warnings=FALSE){
  n <- length(dep)
  p <- ncol(as.matrix(indep))
  clnopt <- min(nclust)
  czmax <- max(nclust)  
  bic <- loglik <- (-1.e9)
  clbic <- rep((-1.e9), czmax)
  eps <- rep(0, czmax)
  rc <- rep(0,(p+1)*czmax)
  dim(rc) <- c(p+1,czmax)
  rv <- rep(0,czmax)
  z <- rep(0, n*czmax)
  dim(z) <- c(n,czmax)
  for (cln in nclust){
    for (i in 1:ir){
      cat("Iteration ",i," for ",cln," clusters.\n")
      emi <- regem(indep, dep, m=randcmatrix(n,cln,p), cln=cln,
                         icrit=icrit, minsig=minsig, warnings=warnings)
      if (emi$warn)   
        emi <- regem(indep, dep, m=randcmatrix(n,cln,p), cln=cln,
                           icrit=icrit, minsig=minsig, warnings=warnings)
      if (!emi$warn){
        bicval <- 2*emi$loglik - log(n)*((p+3)*cln-1)
        if (bicval > clbic[cln])
          clbic[cln] <- bicval
        if (bicval > bic){
          clnopt <- cln
          bic <- bicval
          loglik <- emi$loglik
          eps[1:cln] <- emi$eps
          rc[,1:cln] <- emi$coef
          rv[1:cln] <- emi$var
          z[,1:cln] <- emi$z
        } # if bicval>bic
      }   # if no warning
    }     # for i
  }       # for cln
  g <- c()
  for (i in 1:n)
    g[i] <- which.max(z[i,1:clnopt])
  out <- list(clnopt=clnopt, loglik=loglik, bic=clbic,
              coef=rc[,1:clnopt], var=rv[1:clnopt], eps=eps[1:clnopt], 
              z=z[,1:clnopt], g=g)
  out
# clnopt: Optimal number of clusters, loglik: Loglikelihood, bic: Vector of
# BIC values, coef: Regression coefficients, var: Error variances: 
# eps: cluster proportions, z:a posteriori probabilities, g:optimal
# classification
}          
        





























