function(                                                      # PLOT ACTIVITY/STATE OVER TIME FOR GROUPS OF SIMILAR UNITS
         df=factor(levels=dataframes(eN))[[1]],                  # Dataframe                      # data frame to plot
         threshold=0.0,                                           # Pick units by threshold   # units exceeding this threshold will normally be included
         min_n_units=1,                                          # Pick units by threshold   # additional units will be shown to reach this minimum number of units
         max_n_units=32,                                         # Pick units by threshold   # units will be dropped so as to not exceed this maximum number of units
         plot.type=factor("l",levels=c("b","l","p"))[[1]],         # Cosmetic               # l = lines, p = points, b= both
         line.width=3,                                             # Cosmetic                # width of lines in plot
         panels=factor(NULL,levels=colnames(df)),                # Multi-panel               # Plot multiple panels
         min.y=eN_model_parameter(df,"min_act"),                                              # Cosmetic                # lower limit on y-axis
         max.y=eN_model_parameter(df,"max_act")                                               # Cosmetic                # upper limit on y-axis
) {
  tryCatch(
    {    
#         exactly=factor(NULL,levels=eN_colnamesplit(df)$levels), # Pick units by name               # if set with non-zero length, IGNORE OTHER SETTINGS and plot only these units
#         disprefer=factor(NULL,levels=eN_colnamesplit(df)$levels),# Exceptions to threshold   # units that will be dispreferred to all other units, and only included to meet min_n_units
#         always=factor(NULL,levels=eN_colnamesplit(df)$levels),  # Exceptions to threshold     # units that will be included, regardless of max_n_units
#         prefer=factor(NULL,levels=eN_colnamesplit(df)$levels),  # Exceptions to threshold     # units that will be preferred to all other units, but least active will be excluded to meet max_n_units
#         never=factor(NULL,levels=eN_colnamesplit(df)$levels),   # Exceptions to threshold   # units that will be excluded, regardless of min_n_units
  exactly=disprefer=always=prefer=never=NULL
  require(lattice)
  cols<-eN_colnamesplit(df)

  panels<-as.character(panels)
 #print(panels)
 if(length(panels)==0)
 {
   panels<-setdiff(colnames(df)[!cols$levelcolb],"time")
 }
 if(length(panels)==1)
 {
   df[panels]<-ordered(df[[panels]],levels=unique(df[[panels]]))
 if(panels[1]=="") panels<-NULL
 } 
  if(min_n_units<1) min_n_units=1;



#  pdat<-reshape(df,direction="long",varying=list(names(df)),times=cols$levels,v.names="value",timevar="unit",idvar="time") 


   cnt<-df[cols$levelcolb]

   gpc<-data.frame()
   gpn<-list()
   for(n in 1:ncol(cnt))
   {
     cmp<-unlist(lapply(gpc,function(j,i) all(j==i),cnt[n]))
     if(!any(cmp)) {
       if(ncol(gpc)==0) { gpc<-data.frame(cnt[n])
       } else {      gpc<-cbind(gpc,cnt[n])}
       gpn<-c(gpn,cols$levels[n])
     } else {
       gpn[cmp]<-paste(gpn[cmp],cols$levels[n])       
     }
   }

   gpn<-paste(
    lapply(gpn,function(x) substr(x,1,40)),sep="",
    ifelse(nchar(gpn)>40,".... etc.","")
   )

   gpn<-gsub("([^ ]*) ([^ ]*) ","\\1 \\2\n",gpn)


   gpn<-unlist(gpn)
   colnames(gpc)<-gpn

   if(length(colnames(df))>2)
     inter<-data.frame(unit=colnames(gpc),
                      x=unlist(lapply(gpc,
                               function(x) if(is.numeric(x)) max(x,na.rm=T) else NA 
                              ))
                      )
   else
     inter<-data.frame(unit=cols$levels,x=max(gpc[,2]))


  interest<-which(inter$x>threshold)

#  print("int")
#  print(inter)

  if(min_n_units<=max_n_units) { # only use limits if sensible
    if(length(interest)<min_n_units) { # if threshold yields too few units
      interest<-which(rank(-inter$x,ties.method="min")<=min_n_units) # try and pick at least min units, keeping ties together 
      if(length(interest)>max_n_units) interest<-which(rank(-inter$x,ties.method="max")<=max_n_units) # if that's too many, drop the tie that forces it over
      if(length(interest)<min_n_units) interest<-which(rank(-inter$x,ties.method="first")<=max_n_units) # if that's too few, go up to max, selecting from that tie
    } else if(length(interest)>max_n_units) { # if threshold yields too many units
      interest<-which(rank(-inter$x,ties.method="max")<=max_n_units) # try and pick at most max units, dropping ties together
      if(length(interest)<min_n_units) interest<-which(rank(-inter$x,ties.method="min")<=min_n_units) # if that's too few, add the next block of tied units
      if(length(interest)>max_n_units) interest<-which(rank(-inter$x,ties.method="first")<=max_n_units) # if that's too many, drop dome from that tie down to max
    }
  }

   relcols=colnames(gpc)[interest]

   gpc<-gpc[,interest]



#  rwid=.2

#  crp.rg <- colorRampPalette(c("black", "darkgray","lightgray","blue","purple","red","cyan","green"))
#  colors=crp.rg(length(interest))
   colors=rep(c("black","cadetblue1","chartreuse2","red","purple","darkcyan","orange","lightgray"),times=length(interest)/8+1)
   ltys=rep(rep(each=8,c("solid","dashed","dotted","dotdash")),times=length(interest)/32+1)


#  leveling<-1:length(interest)
#  pdat<-droplevels(subset(pdat,unit%in%interest))
#  relcols=paste(cols$name,interest,sep="")


  pdat<-data.frame()
  if(length(interest)>1)
  {
    pdat<-reshape(cbind(gpc,df[c("time",panels)]),direction="long",
      varying=list(relcols),
      times=gpn[interest],
      v.names="value",
      timevar="unit",
      idvar="time",
      new.row.names=1:((length(interest))*(length(rownames(df))))
    )
  } else {
    pdat<-data.frame(unit=relcols,value=gpc,df[c("time",panels)])
  }
   print("ooo")
    eminy=min(pdat$value)
    emaxy=max(pdat$value)
  if(is.na(min.y))
    min.y=1.05*eminy-.05*emaxy
  if(is.na(max.y))
    max.y=1.05*emaxy-.05*eminy
  if(min.y==max.y)
  {
    min.y=1.05*eminy-.05*emaxy
    max.y=1.05*emaxy-.05*eminy
  }
  if(min.y==max.y)
  {
    min.y=0;
    max.y=1;
  }

#  reo<-order(inter$unit)
#  fs<-order(levels(factor(interest)))
#  lev<-subset(inter[reo,],unit%in%interest)$x

#  print(inter)
  leveling<-order(inter$x[interest],decreasing=T)
#  print(reo)
#  print(lev)
#  print(leveling)
#  print(xyplot(c(1,2)~c(3,4))); return(NULL);

  if(length(panels)>0) {
pnl<-paste(panels,collapse=",")
  fmla<-paste("value~time|paste(",pnl,")")
  print(fmla)
  print(xyplot(
               as.formula(fmla),groups=unit,type=plot.type,lwd=line.width,data=pdat,main=cols$name,
               sub=eN_list_format(attr(df,"settings")),
               col=colors, 
               lty=ltys,
               key=list(text=list(levels(factor(pdat$unit))[leveling]),
                        lines=list(col=colors[leveling],lwd=line.width,lty=ltys[leveling]),
                        space="right",padding.text=5
                       ), 
               ylim=c(min.y,max.y),as.table=T
              )
       )
} else {

  print(xyplot(
               value~time,groups=unit,type=plot.type,lwd=line.width,data=pdat,main=cols$name,
               sub=eN_list_format(attr(df,"settings")),
               col=colors, 
               lty=ltys,
               key=list(text=list(levels(factor(pdat$unit))[leveling]),
                        lines=list(col=colors[leveling],lwd=line.width,lty=ltys[leveling]),
                        space="right",padding.text=5
                       ),   
               ylim=c(min.y,max.y)
              )
       )
}
    }, warning = function(w) { print(w); plt <- print(plot(0,0,"n",xaxt="n",yaxt="n",xlab="",ylab="")); return
    
    }, error = function(e) { print(e); plt <- print(plot(0,0,"n",xaxt="n",yaxt="n",xlab="",ylab="")); return
    
    }, finally = { print ("done");
      
    }
  )       
}
