.packageName <- "tlemix"
#line 1 "/private/tmp/RtmpAzwPqW/R.INSTALL10d63af1/tlemix/R/FLXglm.Density.R"
### #############################################################################################
### #############################################################################################
### ####	  TLE.Flexmix binomial,poisson & gaussian model driver			       ##############
### ####										                                   ##############
### ####	  Density function according current parameter estimate.		       ##############
### ####										                                   ############## 
### ####	  data	- model data, expected to be model.frame with family attribute;##############
### ####	  pars	- model estimate.						                       ##############
### ####		  Current parameters extracted via coef function;		           ##############
### ####                                                                           ##############
### #############################################################################################

flexmix.Density	=function(data,pars,model,family) {
  ##
  ## family <- attr(data,"family") ## family is now user supplied..
  ## this should be written in a nicer way - TODO for the next version)
  est	= coefglm(pars,family) # changed to coefglm # get current parameter estimate
  par	= est$coefficients
  sigma	= est$sigma
  n	= dim(data)[1]
  ## nc	= pars@k
  nc	= length(est$prior)             # number of components
  ##
  ##    family	= attr(data,"family")						# model family
  if(is.null(family)) family="guassian"
  ##
  ## model	= attr(data,"terms")						# model formula
  ##	  attributes(model)=NULL
  ## model	= as.formula(model) # now supplied as argument
  mat	= model.matrix(model,data)      # design matrix
                                        #
  ## this is a work-around to get the response.
  ## is there a better way to do this?
  yresp   = model.response(model.frame(model,data))        # responce variable(s)
  linp    = mat%*%par                   # linear predictor
  if(family=="poisson") lambda = exp(linp) # poisson model parameter
  if(family=="binomial") {prob = exp(linp); prob = prob/(1+prob) # binomial model probability
                          success = yresp[,1]   # number of successes
                          size = rowSums(yresp) # size of cell
                        }
                                        #
                                        #.......................................................................................
                                        #
  ll	= matrix(data =	0, nrow	= n, ncol = nc)
                                        #
  for (i in 1:nc)	{
    if(family=="gaussian") tmp = try(dnorm(yresp,linp[,i],sigma[i]))
    if(family=="poisson")  tmp = try(dpois(yresp,lambda[,i]))
    if(family=="binomial") tmp = try(dbinom(success,size,prob[,i]))
    if(inherits(tmp,"try-error")) return(NULL)
    ll[,i]	= tmp
  }
  if(nc>1)na  = apply(ll,1,function(x)any(is.na(x)))
  ##
  ##       Cluster affiliations and Loglikelihoods by cases:
  ##
  c	= apply(ll,1,function(x)sort(x,index.return=TRUE,decreasing=TRUE)$ix[1]) # cluster
  cc  	= NULL
  for(i in 1:nc) {q=c; q[q!=i]=0; cc=cbind(cc,q)} # clusters as matrix
  cc[cc!=0]=1
  ##
  lik	= apply(ll,1,function(x)sum(est$prior*x)) # Loglikelihood
  if(nc>1)lik[na] = NA
  ##
  return(list(ll=ll,c=c,cc=cc,lik=lik))
}
#line 1 "/private/tmp/RtmpAzwPqW/R.INSTALL10d63af1/tlemix/R/FLXglm.Estimate.R"
#########################################################################################
#######	  TLE-flexmix gaussian, poisson & binomial model driver			#########
#######										#########
#######	  Returns flexmix class estimate.					#########
#######	  data	- model data, expected to be model.frame;			#########
#######	  ind	- data subset;							#########
#######										#########
#######	  ...	- other parameters of current estimate (GLM) procedure.         #########
#######										#########
#######   ntry = maximum number of trials
##########################################################################################

flexmix.Estimate=function(data,ind=NULL,
  nc,class="hard",cluster=NULL,niter=200,minprior=0.1,model=NULL,family=NULL,ntry=9) {
  ##	classification=c("auto", "weighted", "hard", "random")
  ##require(flexmix)                      
                                        
  ## removed getting family as attribute because it is supplied as argument
  ## set default family
  if(is.null(family)) family = "gaussian" # gaussian, poisson and binomial family

  ## model needs to be supplied
  ## if it's not supplied, stop
  if(is.null(model)) {
    stop("No formula supplied.")
  }
                                        
  if(!is.null(ind)) {
    data=data[ind,]                     # substr attribute is not appropriate
    if(!is.null(cluster)) cluster=cluster[ind,]
  } else data=data[,]
  if(!is.null(cluster)) nc=dim(cluster)[2]
  nobs=dim(data)[1]
                                        #
  control	= list(iter=niter,tol=1.e-9,class=class,verb=0,minprior=minprior)

  ## maximum number of trials, now user supplied, default 9
  ## ntry = 9
  itry = 0								# current try number
                                        
  while(itry<=ntry){
    nmix <-	try(flexmix(model,                 
                        data=data.frame(data), 
                        k=nc,                  
                        cluster=cluster, # =NULL	or cluster=d$c,
                        control=control, # control=list(verb=verb,iter=iter,tol=tol,class=class)
                        model=FLXglm(model, family=family)),silent=TRUE)
    if(!inherits(nmix, "try-error")) break
    if(itry>ntry) break                 # no more than ntry tries
    itry<-itry+1
  }
  if(itry>ntry) return(NULL)            # no solution
  stop=attr(nmix,"stop")                # ready
  if(is.null(stop)) stop=F

  if(T&&!stop){ ### 
	if(class!="weighted" && class!="auto")
      {	control$class="weighted"
		control$iter=nmix@iter      # no more iteration than in previous step
                                        #control$iter=3
		tmp_nmix<-nmix

		nk = nmix@k                     # number of components
		c  = nmix@cluster               # cluster
		cc = NULL                       # Problem: it needs to be matrix
		for(i in 1:nk) {q=c; q[q!=i]=0; cc=cbind(cc,q)} # clusters as matrix
		cc[cc!=0]=1

		old.logLik = nmix@logLik
		itry =0                         # current try number
		while(itry<=ntry)
          {	#Sys.sleep(0.01)					# needed for synchronisation of	the other windows
			nmix <-	try(flexmix(model,  # as.matrix(d)~1,
                                data=data.frame(data),
                                k=nk,       #
                                cluster=cc, # =Previous estimate
                                control=control, # control=list(verb=verb,iter=iter,tol=tol,class=class)
                                model=FLXglm(model, family=family)),silent=TRUE)
			if(!inherits(nmix, "try-error")) break
			if(itry>ntry) break         # no more than ntry tries
			itry<-itry+1
          }
		if(itry>ntry) return(NULL)      # no solution
		new.logLik = nmix@logLik
        
        ## if(new.logLik>old.logLik) cat(">")
        ## else if(new.logLik<old.logLik) cat("\<")
        ## else cat("=")
        
      }
  } #..............................................................................
  nc    =nmix@k                   # number of the	components in the estimate
  if(nc==0) return(NULL)          # no solution

  return(nmix)
}
###
### ######################################################################################
### Returns list of coefficents for glm-class models
coefglm=function(nmix,family)
{
	nc = length(nmix@prior)							# number of the	components in the estimate
	if(nc==0) return(NULL)							# no solution
	coef = rep(list(NULL),nc)
	param =	attr(nmix@components[[1]][[1]],"parameters")
	coef = param$coef
	if(family=="gaussian") sigma = param$sigma
	prior =	nmix@prior
	if(nc>1) for( i	in 2:nc)	{
		param=attr(nmix@components[[i]][[1]],"parameters")
		coef = cbind(coef,param$coef)
		if(family=="gaussian") sigma = c(sigma,param$sigma)
	}
	coef = as.matrix(coef)
	colnames(coef) = paste("C",1:nc,sep="")
	names(prior) = colnames(coef)
	if(family=="gaussian") names(sigma) = colnames(coef)

	if(family=="binomial") return(list(coefficients=coef,prior=prior))
	if(family=="poisson")  return(list(coefficients=coef,prior=prior))
	if(family=="gaussian") return(list(coefficients=coef,sigma=sigma,prior=prior))
	return(NULL)
}

#line 1 "/private/tmp/RtmpAzwPqW/R.INSTALL10d63af1/tlemix/R/FLXmclust.Density.R"
##########################################################################################
#######   TLE-FLXmclust - MvtNormal model driver                                ##########
#######                                                                         ##########
#######   Density function according current parameter estimate.                ##########
#######                                                                         ##########
#######   data  - model data, expected to be matrix     ;                       ##########
#######   estim - model estimate.                                               ##########
#######           Current parameters extracted via coef function;               ##########
#######                                                                         ##########
##########################################################################################

### ... is used for user-supplied stuff and to enable TLE to call Density with 4 paramters
FLXmclust.Density =function(data,estim,model,...) {
  mat <- data[[1]]
  ## require(mvtnorm)
  pars    = coefmclust(estim)
  n       = dim(mat)[1]                # number of cases, changed from data to mat
  nc      = length(pars)                # number of components/clusters
                                        #
  lc=ll   = matrix(data = 0, nrow = n, ncol = nc) # loglikelihood by clusters
  prior   = rep(0,nc)

  ## remove old stuff that exptected model.frame
  ## model   = attr(data,"terms")          # model formula
  ## attributes(model)=NULL
  ## model   = as.formula(model) 
  ## changed because if different input parameters...
  ##mat = model.matrix(model,data)        # design matrix
  ## mat     = as.matrix(data) #model.matrix(model,data)                    # design matrix
                                        #
  for (i in 1:nc)
    {       tmp      = try(dmvnorm(mat, pars[[i]]$center, pars[[i]]$cov)) # density
            if(inherits(tmp,"try-error")) return(NULL)
            prior[i] = pars[[i]]$prior
            ll[,i]   = tmp*prior[i]
            lc[,i]  = tmp
          }
  if(nc>1)na  = apply(ll,1,function(x)any(is.na(x)))
  ##
  ##       Cluster affiliations and Loglikelihoods by cases:
  ##
  c	= apply(ll,1,function(x)sort(x,index.return=TRUE,decreasing=TRUE)$ix[1]) # cluster
  cc  	= NULL
  for(i in 1:nc) {q=c; q[q!=i]=0; cc=cbind(cc,q)} # clusters as matrix
  cc[cc!=0]=1
  ##
  lik     = rowSums(ll)
  if(nc>1)lik[na] = NA
  ##
  return(list(ll=ll,lc=lc,c=c,cc=cc,lik=lik))
}
#line 1 "/private/tmp/RtmpAzwPqW/R.INSTALL10d63af1/tlemix/R/FLXmclust.Estimate.R"
##########################################################################################
#######   TLE-FLXmclust - MvtNormal model driver                                 #########
#######                                                                         #########
#######   Returns flexmix class estimate.                                       #########
#######   data  - model data, expected to be model.frame;                       #########
#######   ind   - data subset;                                                  #########
#######                                                                         #########
#######   ...   - other parameters of current estimate (GLM) procedure.         #########
#######                                                                         #########
#########################################################################################

FLXmclust.Estimate=function(data,ind=NULL, 
  nc,class="hard",cluster=NULL,niter=200,minprior=0.1,model=NULL,ntry=9,...) {
  ##       classification=c("auto", "weighted", "hard", "random")
  ## require(flexmix)
  ##
  if(is.null(model)) {                  # extract model formula
    ## model = attr(data,"terms")
    ## attributes(model)=NULL
    ## model = as.formula(model)
    stop("Model formula must be supplied.")
  }
  ##
  ## create model matrix, but data must be a matrix, not a data.frame
  ## TODO: ensure that this is the case!
  ## mat=data
  ## mat = as.matrix(data)                 # user-supplied data frame
  var_name <- names(data) # extract variable name
  mat <- data[[1]] # get data from list
  if(!is.null(ind)) {
    mat=mat[ind,]                       # substr attribute is not appropriate
    if(!is.null(cluster)) cluster=cluster[ind,]
  } 

  else mat=mat[,]

  if(!is.null(cluster)) nc=dim(cluster)[2]
  nobs=dim(mat)[1]
                                        
  if(is.null(nobs)) return(NULL)
                                        
  control	= list(iter=niter,tol=1.e-9,class=class,verb=0,minprior=minprior)

  ## build data list again
  data <- list(mat)
  names(data) <- var_name
  
  ## moved this to function call, should be user supplied
  ##  ntry = 9
  ## maximal number of trials
  itry = 0                              # current try number
                                        
  while(itry<=ntry) {
    nmix <- try(flexmix(model, #mat~1,                                      # as.matrix(d)~1,
                        data=data,
                        k=nc,            #
                        cluster=cluster, # =NULL or cluster=d$c,
                        control=control, # control=list(verb=verb,iter=iter,tol=tol,class=class)
                        model=FLXmclust(diag=FALSE)),silent=TRUE)
    if(!inherits(nmix, "try-error")) break
    if(itry>ntry) break                 # no more than ntry tries
    itry<-itry+1
  }
	
			
			
	
  if(itry>ntry) return(NULL)            # no solution
  stop=attr(nmix,"stop")                # ready
  if(is.null(stop)) stop=F

  ##................................................................................
  if(T&&!stop) { # the last iterations must be carried out with posterior weights!
	if(class!="weighted" && class!="auto") 
      {	control$class="weighted"
		control$iter=nmix@iter      # no more iteration than in previous step
                                        #control$iter=3	
		tmp_nmix<-nmix
		
		nk = nmix@k                     # number of components
		c  = nmix@cluster               # cluster 
		cc = NULL                       # Problem: it needs to be matrix
		for(i in 1:nk) {q=c; q[q!=i]=0; cc=cbind(cc,q)} # clusters as matrix
		cc[cc!=0]=1
		
		old.logLik = nmix@logLik
		itry =0                         # current try number
		while(itry<=ntry)
          ##nmix <-	try(flexmix(mat~1,				# as.matrix(d)~1,
          nmix <- try(flexmix(model,
                              data,
                              k=nk,       # 
                              cluster=cc, # =Previous estimate
                              control=control, # control=list(verb=verb,iter=iter,tol=tol,class=class)
                              model=FLXmclust(diag=FALSE)),silent=TRUE)
        if(!inherits(nmix, "try-error")) break
        if(itry>ntry) break             # no more than ntry tries
        itry<-itry+1
      }
    if(itry>ntry) return(NULL)          # no solution
    new.logLik = nmix@logLik
    ## if(new.logLik>old.logLik) cat(">")
    ## else if(new.logLik<old.logLik) cat("\<")
    ## else cat("=")
  }
                                        #..............................................................................
  nc    =nmix@k                 # number of the	components in the estimate
  if(nc==0) return(NULL)        # no solution

  return(nmix)
}

###
#########################################################################################
### # FLXmclust coef function

coefmclust <- function(nmix) {
  nc = length(nmix@prior)          # number of the components in the estimate
  if(nc==0) 
    return(NULL)                        # no solution
  coef    = rep(list(NULL),nc)
  names(coef) = paste("C",1:nc,sep="")
  for( i in 1:nc) {       
    par     = parameters(nmix,component=i,model=1,simplify=FALSE)
    size    = nmix@size[i]
    names(size) = NULL
    coef[[i]] = list(center=par$center,cov=par$cov,prior=nmix@prior[i],size=size)
  }
  return(coef)
}




#line 1 "/private/tmp/RtmpAzwPqW/R.INSTALL10d63af1/tlemix/R/TLE.R"
###########################################################################################
###   version info removed, hae a look at the svn info
###	TLE procedure
###
###		kStar	= k*- size of the initial random subsample;
###		kTrim	= trimming parameter: size of the C-steps random subsample;
###		data	= data	frame containing the x and y variables;
###	.................................................................................
###		nit	= number of iterations;
###		msglvl	= level of messages;
###		result	= restart/continuation information;
###		cit	= number of iteration in refinement steps;
###		test	= expected true loglikelihood of the model. stop procedure if reached;
###	.................................................................................
###		Density	= density function of type - function(data,solution);
###		Estimate = specific estimation procedure interface: function(data,ind,...);
###		...	= specific parameters for Estimate procedure.
###
#########################################################################################
###
###
### Changelog:
### 2008-12-16:  - added sensible default parameters
###              - changed number of iterations from 100 to 10
###
##########################################################################################
##### setClass TLE
##########################################################################################
setClass("TLE",
		representation(estimate="ANY",
						iterbest="numeric",
						it="numeric",
						maxloglik="numeric",
						indbest="integer",
						indout="integer",
						tleweights="numeric",
						tlelogliks="numeric",
						tleclusters="ANY",
						kTrim="numeric",
						kStar="numeric",
						mcomp="ANY",
						nobs="integer",
						stop="logical",
						call="call"))

##########################################################################################
##### setMethod show for TLE object
##########################################################################################
if(!isGeneric("show")) setGeneric("show")
setMethod("show", "TLE",
          function(object){
            cat("\nCall:", deparse(object@call,0.75*getOption("width")),
                sep="\n")
            cat("\nEstimation:\n")
            print(object@estimate)

          })

##########################################################################################
##### setClass summary.TLE - class for summary object
##########################################################################################
setClass("summary.TLE",
         representation(call="call",
                		estimate="ANY",
                        kTrim="numeric",
                        kStar="numeric",
                        maxloglik="numeric",
                        nobs="integer",
                        nout="integer"))

##########################################################################################
##### setMethod for show() of summary.TLE objects
##########################################################################################
setMethod("show", "summary.TLE",
          function(object){
            cat("\nCall:", deparse(object@call,0.75*getOption("width")),
                sep="\n")
            cat("\n")
            print(object@estimate)
            cat("\n")
            cat("kTrim:", object@kTrim, "   Number of Observations:", object@nobs, "   Number of Outliers:",object@nout,"\n")
            cat("\n")       
          })

#######################################
### plot
#######################################
### In vielen Fällen möchte man auch plot-Methoden für eine S4-Klasse erstellen.
### Dies funktioniert z.B. folgendermaßen:
if(!isGeneric("plot")) setGeneric("plot")
setMethod("plot", signature("TLE","missing"), 
	#function(x, y = NULL, ...){ 
	function(x,y=NULL,...) {
	plot(x@estimate)
})



##########################################################################################
##### setGeneric function TLE
##########################################################################################
setGeneric("TLE",
        function(formula,family,data, kStar=NULL,kTrim=NULL, nit=10, msglvl=0, result=NULL, cit=9, test=NULL, nc, Density, Estimate, ...)
        standardGeneric("TLE"))

### just to make sure that some S3 generics are available in S4
setGeneric("summary")



### #######################################################################################
### #######################################################################################
### #######################################################################################
### ## setMethod TLE
### #######################################################################################
### if called with existing TLE object (as result=), then restart estimate...
### #######################################################################################
setMethod("TLE",
          signature(formula="formula",family="character",data="ANY",kStar="ANY",
                    kTrim="ANY", nit="ANY", msglvl="ANY", cit="ANY",
                    result="ANY",test="ANY",nc="numeric", Density="function", Estimate="function"),
          function(formula,family,data, kStar=NULL,kTrim=NULL, nit=100, msglvl=0, result=NULL,
                   cit=9, test=NULL,nc, Density, Estimate, ...)
          {
			
            mycall = match.call()
            nobs = dim(data)[1]         # number of observations
            dimension = dim(data)[2]    # number of dimensions

            ## proposed default parameter
            k = (nobs + nc * (dimension + 1) + 1)/2
            if(is.null(kTrim)) {
              kTrim = k;
            }
            if(kTrim<k) {
              warning("kTrim should be greater than (nobs+nc*(dimension+1)+1)/2")
            }

          
            if(is.null(kStar)) {
              ## set kStar to the average of kTrim and g(p+1).
              ## round to nearest integer value.
              kStar = round((kTrim + (nc*(dimension+1)))/2,0)
            }
            ## cat(paste("k* and ktrim chosen: ",kStar,kTrim))
        


            ## check restart object	
            ## scan if current run is restart run			
            if( 	is.null(result) ||
               is.null(result@kTrim) || 
               result@kTrim!=kTrim ||
               is.null(result@nobs) || 
               result@nobs!=nobs ||
               is.null(result@maxloglik) ||
               is.null(result@estimate)   ) {	
              if(msglvl>0) print("New estimate start")
              estimate	= NULL          # current model estimate
              iterbest	= 0				# iteration of the best estimate
              it        = 0         # current iteration number
              maxloglik	= NULL          # maximal loglikelihood
              mcomp		= 1				# initial number of clusters if any
            } else {                    # no, it's a restart...
              if(msglvl>0) print("Restart estimate!")
              ## save previous estimate elements
              estimate	= result@estimate # iteration of the best estimate
              iterbest	= result@iterbest
              it			= result@it          # current iteration number
              maxloglik	= result@maxloglik       # maximal loglikelihood
              indbest		= result@indbest     # selected best subsample
              indout		= result@indout      # outliers
              tleweights	= result@tleweights  # tle-weights
              tlelogliks	= result@tlelogliks  # tle-logliks
              tleclusters	= result@tleclusters # tle-clusters
              mcomp		= result@mcomp           # number of components
            }
            Qnew	= NULL              # loglikelihood of selected data

            ## if model based clustering, we have to create a list
            ## this is some kind of workaround
            if(family=="mvtnormal") {
              var_name <- as.character(substitute(data))
              data <- list(data)
              names(data) <- var_name
            }
            ## #######################################################################################
            ## #	main TLE loop
            ## #######################################################################################
            oldit=it    # set oldit to it (0 if new estimate, else result@it)
            for(it in (1:nit+oldit)) # 'nit' new iteratins, (1+oldit):(nit+oldit)
              {
                ##	set the	initial	random subsample of rank k*:
                ind <-sort(sample(nobs,kStar))
                                        
                ##	control	prints:
                if(msglvl==1) {	if((it-1)%%80==0) cat("\n.") else cat(".") }
		        if(msglvl>1) cat("\n iteration:",it,"\n")
                                        

                ## Initial Step. Estimate using random subsample of size k*
                est	= Estimate(data,ind,nc=nc,model=formula,family=family, ...)
                if(is.null(est)) {
                  next() # if we don't get an estimate, go to the next iteration
                }

                ## get estimated densities for all data; update: added model and family
                wrk	= Density(data,est,model=formula,family=family)		
                if(is.null(wrk)) {
                  next() # if we don't get an estimate, go to the next iteration
                }
                loglik	= log(wrk$lik) # estimated Loglikelihood vector
                ## the difference between $c and $cc is that $cc is in matrix form whereas $c is just a vector
                ## ie $cc has one column for every component;
                cluster	= wrk$c			# estimated clusters
                ccomp	= wrk$cc		# estimated clusters
                ncomp	= dim(wrk$cc)[2] # number of components
                                        #
                qtmp	= loglik[ind] ## loglikelihood of the data points in the subsample
                if(any(is.na(qtmp))) next()	# incorrect estimate
                Qold	= sum(qtmp)         # LogLikelihood of selected data
                                        
                if(msglvl>1) cat("\n Initial sample of size   ",kStar,
                                 "\n Estimated LogLik=	        ",Qold,
                                 "\n sample observations:       ",ind,
                                 "\n")
                ## ###################################################################################
                ## C-steps:	##########################################################################
                for(cstep in 1:cit) {
                  ## if(msglvl>1)progress.bar((it-oldit-1)*9+cstep-1,nit*9+1,title="TLE progress")
                  ## Improvement:
                  ## Get kTrim best cases.
                  ## Replace all NA values with machine minimal value,
                  ## in order to obtain correct kTrim indeces:
                  ##	if(any(is.na(loglik))) loglik[is.na(loglik)]=.Machine$double.xmin
				
                  kSort	= sort(loglik,index.return=T,decreasing=TRUE)$ix[1:kTrim]
				
                  ## Use appropriate estimation function interface
                  est	= Estimate(data,kSort,cluster=ccomp, nc=nc,model=formula,family=family,...) # added model
                  if(is.null(est)) break()
                  ##
		  
                  wrk	= Density(data,est,model=formula,family=family) # get estimate densities for all data, added model and family
                  if(is.null(wrk)) break()
                  loglik	= log(wrk$lik) # estimated Loglikelihood vector
                  cluster	= wrk$c        # estimated clusters
                  ccomp	= wrk$cc           # estimated clusters
                  ncomp	= dim(wrk$cc)[2]   # number of components
                                        
                  qtmp	= loglik[kSort]
                  if(any(is.na(qtmp))) break() # incorrect estimate
                  Qnew	= sum(qtmp)     # LogLikelihood of selected data
                                        
                  if(is.na(Qnew)) break() # If not all logliks are incorrect (Binary model!)
                                        
                  totlik	= sum(loglik,na.rm=TRUE) # LogLikelihood of all data
                  weights	= rep(0,nobs)            # set all weights to 0
                  weights[kSort] = 1 # set weights of objects in the subsample to 1; TLE weights
                                        #
                  if(msglvl>1) cat("\n C-step",cstep," of size ",length(kSort),
                                   "\n Estimates LogLik=	        ",Qnew,
                                   "\n sample observations:       ",sort(kSort),
                                   "\n")
                                        
                  if( is.null(maxloglik) || Qnew > maxloglik ) { # save new "best" attributes
					maxloglik	= Qnew   # LogLikelihood of the best estimate
					iterbest	= it     # iteration of the best subsample
					totloglik	= totlik # LogLikelihood of all data
					tleweights	= weights     # TLE weights
					tlelogliks	= loglik      # TLE logliks
					tleclusters	= cluster     # TLE clusters
					indbest		= sort(kSort) # ordered indexes of best subsample
					indout		= (1:nobs)[-kSort] # indices of outliers
					indout		= sort(indout)     # sorted indices outliers
					estimate	= est              # estimate & coefficients
					if(is.null(ncomp)) ncomp=1
					mcomp		= ncomp # number of clusters/components
                                        
					if(msglvl>0)   cat("\n =========================================",
                                       "\n loglik updated!",maxloglik,iterbest,mcomp, 
                                       if(is.null(test)) "" else paste("(",format(maxloglik-test,digits=3),")",sep=""),
                                       "\n =========================================",
                                       "\n")
                  }
                                        
                  if(msglvl>1 && length(ind)==kTrim && all(ind==sort(kSort))) {
					print("TLE	C-step converged"); break () 
                  }
                  ind<-sort(kSort)
                  Qold = Qnew
                } 
                ## ########################
                ## End of the Improvement loop (C-steps)
                ## ########################
                stop=!is.null(test)&&!is.null(Qnew)&&test<=Qnew
                if(stop) {
                  print("target OK")
                  break()
                }
              }
            ## End of the main	iteration loop
            ##
            ##	Final print:
            if(msglvl>1) {	
              if(!is.null(estimate))
                {	
                  li=length(indbest)
                  if(!is.null(maxloglik))
                    {     
                      cat("\n best subsample found:",
                          "\n ==========================================================================",
                          "\n loglik	          = ",maxloglik,
                          "\n num of components = ",mcomp,
                          "\n iter of best           =	",iterbest,
                          "\n total num.iter.     =	",it,
                          "\n sorted subsample   =	",indbest[1:(li%/%4)],
                          "\n		          =	",indbest[(   li%/%4+1):(2*li%/%4)],
                          "\n		          =	",indbest[(2*li%/%4+1):(3*li%/%4)],
                          "\n		          =	",indbest[(3*li%/%4+1):li],
                          "\n outliers	          =	",indout,
                          "\n ==========================================================================",
                          "\n")
                    }
                } else
              {	cat(	"\n No improvement \n")
              }
            }
            cat("\n")
            ## if(msglvl>1)progress.bar(1,0)
            if(!is.null(estimate))      # there is an estimate
              {	
                result = new("TLE", estimate=estimate,
                  iterbest=iterbest,it=it, maxloglik=maxloglik, indbest=indbest, indout=indout,
                  tleweights=tleweights, tlelogliks=tlelogliks, tleclusters=tleclusters,
                  kTrim=kTrim, kStar=kStar, mcomp=mcomp, nobs=nobs, stop=stop)
                ## result=list(estimate=estimate,
                ## iterbest=iterbest,it=it,maxloglik=maxloglik,indbest=indbest,indout=indout,tleweights=tleweights,
                ## kTrim=kTrim,kStar=kStar,mcomp=mcomp,nobs=nobs)
                ## attr(result,"stop")=stop
                result@call <- mycall
                result
              }
            else return(new("TLE"))
          })


### #################################################################################################
### ####### summary method for TLE object
### #################################################################################################
setMethod("summary", "TLE",
          function(object, ...) {       
            z <- new("summary.TLE",
                     call = object@call,
                     estimate = object@estimate,
                     kTrim = object@kTrim,
                     kStar = object@kStar,
                     ## glik = object@maxloglik,
                     nobs = object@nobs,
                     nout = length(object@indout))
            z
          } )


if(!isGeneric("estimate"))
	setGeneric("estimate", function(object) standardGeneric("estimate"))

### # Set method for estimate() function used to return the estimate
setMethod("estimate", signature="TLE", 
          function(object) 
          {
			object@estimate
          } )
#line 1 "/private/tmp/RtmpAzwPqW/R.INSTALL10d63af1/tlemix/R/tleplot.R"
setGeneric("tleplot",
		function(object,data, ...)
			standardGeneric("tleplot"))


setMethod("tleplot",
		signature(object="TLE",data="data.frame"),	
		function(object, data,...) {
			nobs <- dim(data)[1] # number of observations
			outliers <- object@indout # index of outliers
			normalobs <- object@indbest # "normal" observations
			col <- rep(1,nobs) # create initial color vector
			
			
			# set color to clusternumber of the observations
			# that where selected by flexmix
			#col[normalobs] <- object@estimate@cluster+1
			col[normalobs] <- object@tleclusters[object@indbest]+1
			
			## assert: col is a vector of length nobs containing
			# ones (=black) for outliers, and > 1 for normal observations
			dottype <- rep(16,nobs)
			dottype[outliers] <- rep(2,length(outliers))
			# draw a scatterplot for twodimensional
			#### CAVE, TODO: x,y selected by hand, needs to be generalized
			plot(data$x,data$y,col=col,pch=dottype,...)
			
			## Legend
			cf <- as.factor(object@estimate@cluster)
			n <- length(levels(cf))
			string <- character(length=n+1)
			for(i in 1:n) {
				string[i] <- paste("Cluster", i)
			}
			string[n+1] <- "Outliers"
			legend("bottomright", string, col=(n+1):1, pch=cbind(rep(1,n),2))
		})
