function(                                                      # PLOT ACTIVITY/STATE OVER TIME FOR INDIVIDUAL UNITS
         df=factor(levels=dataframes(eN))[[1]],                  # Dataframe                      # data frame to plot
         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
         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
         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
         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
         panels=factor(NULL,levels=colnames(df)),                # Multi-panel               # Plot multiple panels
         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
         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
         ythresh=eN_model_parameter(df,"id_conn::threshold")             # Cosmetic                  # Where to draw horizontal threshold
) {
  tryCatch(
    {
      
      #   if (is.null(ythresh) |  is.na(ythresh)) {
      #     lazyNut_exec("((default_model) parameters) copy tmp.params")
      #     param.df<-eN["tmp.params"]
      #     parameters<-param.df["parameter"]
      #     threshold_param<-parameters[grep("threshold",parameters[,1]),]
      #     ythresh=eN_model_parameter(df,threshold_param[1])
      #   }
        
        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)>0)
       {
         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") 
      
         if(length(colnames(df))>2)
           inter<-data.frame(unit=cols$levels,
                            x=unlist(lapply(df,
                                     function(x) if(is.numeric(x)) max(x,na.rm=T) else NA 
                                    ))
                            [cols$levelcolb])
         else
           inter<-data.frame(unit=cols$levels,x=max(df[,2]))
      
         
      
      
      #  inter=with(pdat,aggregate(value,list(unit=unit),FUN=max))
        exactly=exactly[exactly%in%inter$unit]
      
        if(length(exactly)==0) {
      
          always=always[always%in%inter$unit]
          never=never[never%in%inter$unit]
      
      
          prefer=prefer[prefer%in%inter$unit&!prefer%in%c(always,never)]
          disprefer=disprefer[disprefer%in%inter$unit&!disprefer%in%c(always,never)]
      
      
          if(length(always)>max_n_units) max_n_units<-length(always)
      
          if((length(inter$x)-length(never))<min_n_units) min_n_units<-length(inter$x)-length(never)
      
      
          rang<-max(inter$x)-min(inter$x)+1
      
          inter$x[inter$unit%in%always]=inter$x[inter$unit%in%always]+10*rang
          inter$x[inter$unit%in%never]=inter$x[inter$unit%in%never]-10*rang
      
          inter$x[inter$unit%in%prefer]=inter$x[inter$unit%in%prefer]+rang
          inter$x[inter$unit%in%disprefer]=inter$x[inter$unit%in%disprefer]-rang
        } else {
          threshold=0
          min_n_units=max_n_units=length(exactly)
          inter$x[inter$unit%in%exactly]=1
          inter$x[!inter$unit%in%exactly]=-1
        }
      
        interest<-subset(inter,x>threshold)$unit
      #  if(length(lev>0)) 
      
        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<-subset(inter,rank(-x,ties.method="min")<=min_n_units)$unit # try and pick at least min units, keeping ties together 
            if(length(interest)>max_n_units) interest<-subset(inter,rank(-x,ties.method="max")<=max_n_units)$unit # if that's too many, drop the tie that forces it over
            if(length(interest)<min_n_units) interest<-subset(inter,rank(-x,ties.method="first")<=max_n_units)$unit # 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<-subset(inter,rank(-x,ties.method="max")<=max_n_units)$unit # try and pick at most max units, dropping ties together
            if(length(interest)<min_n_units) interest<-subset(inter,rank(-x,ties.method="min")<=min_n_units)$unit # if that's too few, add the next block of tied units
            if(length(interest)>max_n_units) interest<-subset(inter,rank(-x,ties.method="first")<=max_n_units)$unit # if that's too many, drop dome from that tie down to max
          }
        }
      
      #  rwid=.2
      
        interest<-droplevels(interest)
      #  crp.rg <- colorRampPalette(c("black", "darkgray","lightgray","blue","purple","red","cyan","green"))
      #  colors=crp.rg(length(interest))
         colors=rep(c("black","blue","darkgreen","red","purple","cyan","orange","lightgray"),times=8)
         ltys=rep(each=8,c("solid","dashed","dotted","dotdash"))
      
      #  leveling<-1:length(interest)
      #  pdat<-droplevels(subset(pdat,unit%in%interest))
        relcols=paste(cols$name,interest,sep="")
        pdat<-reshape(df[c(relcols,"time",panels)],direction="long",
          varying=list(relcols),
          times=interest,
          v.names="value",
          timevar="unit",
          idvar="time",
          new.row.names=1:((length(relcols))*(length(rownames(df))))
        ) 
          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
        leveling<-order(lev,decreasing=T)
      #  print(reo)
      #  print(lev)
      #  print(leveling)
      
        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"
                             ), 
                     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"
                             ),   
                     ylim=c(min.y,max.y),
                     panel=function(...)
      {
        panel.xyplot(...)
        panel.refline(h=ythresh,col="black",lty="dotted",lwd=5)
      }
                    )
             )
      }
    }, warning = function(w) { print("warning"); print(w); print(plot(0,0,"n",xaxt="n",yaxt="n",xlab="",ylab="")); return
      
    }, error = function(e) { print("error"); print(e); print(plot(0,0,"n",xaxt="n",yaxt="n",xlab="",ylab="")); return
      
    }, finally = { print("done"); 
      
    }
  )
  
  
  
}
