qemegs.anovaplot <- function(datafr, grouping, svgfile, plottype = "bar", cmap = "heat", errorbar=TRUE, errorbartype="Standard error", errorbarfactor=1) 
{
    library(reshape2);
    if (R.Version()$os=="linux-gnu")
    {
        library(RSvgDevice);
    }

    
    debug = FALSE
    
    idvars = colnames(datafr);
    idvars = idvars[-length(idvars)];
    datafrmelt = melt(datafr, id.vars=idvars, value.name="dv");
    castform = grouping[1];
    plottitle = grouping[1];
    if (length(grouping)==1)
    {
        castform = paste(castform, "~ .");
    }
    else
    {
        for (i in 2:length(grouping))
        {
            castform = paste(castform, "~");
            castform = paste(castform, grouping[i]);
            plottitle = paste(  plottitle, " x ",sep="");
            plottitle = paste(plottitle, grouping[i],sep="");
        }
    }
    datafrmeltcastmean = acast(datafr, castform, fun.aggregate=mean, value.var="dv");
    datafrmeltcastsd = acast(datafr, castform, fun.aggregate=sd, value.var="dv");
    datafrmeltcastlength = acast(datafr, castform, fun.aggregate=length, value.var="dv");
    datafrmeltcastsqrtn= sqrt(datafrmeltcastlength);
    
    
    
    if (errorbartype=="Standard error")
    {
        datafrmeltcasterror = errorbarfactor * (datafrmeltcastsd / datafrmeltcastsqrtn);
    }
    else if (errorbartype=="Confidence interval")
    {
        datafrmeltcasterror = errorbarfactor * qt(0.975, datafrmeltcastlength-1) * (datafrmeltcastsd / datafrmeltcastsqrtn);
    }
    else if (errorbartype=="Standard deviation")
    {
        datafrmeltcasterror = errorbarfactor * datafrmeltcastsd;
    }
    else
    {
        print("error bar type not found!")
        return()
    }
    datafrmeltcasterrorbackup = datafrmeltcasterror;
    
    
    
    maxval = max(datafrmeltcastmean + datafrmeltcasterror);
    minval = min(datafrmeltcastmean - datafrmeltcasterror);
    if (plottype=="bar")
    {
        if ((maxval>0 && minval>0) || (maxval<0 && minval<0))
        {
            minval = 0;
        }
    }
    
    range = abs(maxval - minval);
    YLimVec = c(minval- range*0.1, maxval+ range*0.1);
    
    colInds = which(colnames(DataFr)==grouping[1]);
    for (i in 2:length(grouping))
    {
        colInds = c(colInds, which(colnames(DataFr)==grouping[i]));
    }
    YLabel = "Regional amplitude [MicroV]";
    
    # 1 grouping factor
    if (length(grouping)==1)
    {
        if (!debug)
        {
            if (R.Version()$os=="linux-gnu")
            {
                devSVG(svgfile, width=7,  height=7)
            }
            else
            {
                svg(svgfile, width=7,  height=7, onefile=TRUE);
            }
        }
        else
        {
            dev.new(width=7,  height=7)
        }
        
        figmatrixdata = c(0, 0.8, 0.05, 0.95,     0.8, 1.0, 0.05, 0.95);
        figmatrix = matrix(figmatrixdata, nrow=2, ncol=4, byrow=TRUE);
        split.screen(figmatrix);
        
        screen(1)
        if (PlotType=="bar")
        {
            CategoryNames = colnames(datafrmeltcastmean);
        }
        else
        {
            CategoryNames = rownames(datafrmeltcastmean);
        }
        ColorTmp = qemegs.oneplot(datafrmeltcastmean, cmap, errorbar, datafrmeltcasterror, plottype, "", YLabel, YLimVec, CategoryNames);
        
        title(outer=TRUE, adj=0.5, main = list(plottitle, cex=1.1, col="black", font=2), line = -2) 
        
        screen(2)
        if (PlotType=="bar")
        {
            add_legend("topright", rownames(datafrmeltcastmean),  fill=ColorTmp,  cex=0.8, inset=c(0.02, 0.3))
        }
        else
        {
            add_legend("topright", colnames(datafrmeltcastmean),  fill=ColorTmp,  cex=0.8, inset=c(0.02, 0.3))
        }
        close.screen(all =TRUE);
        
        if (!debug)
        {
            dev.off();
        }
    }
    # 2 grouping factors
    else if (length(grouping)==2)
    { 
        
        if (!debug)
        {
            if (R.Version()$os=="linux-gnu")
            {
                devSVG(svgfile, width=10,  height=7)
            }
            else
            {
                svg(svgfile, width=10,  height=7, onefile=TRUE);
            }
        }
        else
        {
            dev.new(width=10,  height=7)
        }
        
        
        
        figmatrixdata = c(0, 0.8, 0.05, 0.95,     0.8, 1.0, 0.05, 0.95);
        figmatrix = matrix(figmatrixdata, nrow=2, ncol=4, byrow=TRUE);
        split.screen(figmatrix);
        
        screen(1)

        if (PlotType=="bar")
        {
            CategoryNames = colnames(datafrmeltcastmean);
        }
        else
        {
            CategoryNames = rownames(datafrmeltcastmean);
        }

        ColorTmp = qemegs.oneplot(datafrmeltcastmean, cmap, errorbar, datafrmeltcasterror, plottype, "", YLabel, YLimVec, CategoryNames);
        
        title(outer=TRUE, adj=0.5, main = list(plottitle, cex=1.1, col="black", font=2), line = -2) 
        
        screen(2)
        if (PlotType=="bar")
        {
            add_legend("topright", rownames(datafrmeltcastmean),  fill=ColorTmp,  cex=0.8, inset=c(0.02, 0.3))
        }
        else
        {
            add_legend("topright", colnames(datafrmeltcastmean),  fill=ColorTmp,  cex=0.8, inset=c(0.02, 0.3))
        }
        close.screen(all =TRUE);
        
        
        if (!debug)
        {
            dev.off();
        }
    }
    # 3 grouping factors
    else if (length(grouping)==3)
    {
        
        
        NCols = dim(datafrmeltcastmean)[3];
        
        if (!debug)
        {
            if (R.Version()$os=="linux-gnu")
            {
                devSVG(svgfile, width=4*NCols,  height=7);
            }
            else
            {
                svg(svgfile, width=4*NCols,  height=7, onefile=TRUE);
            }
        }
        else
        {
            dev.new(svgfile, width=4*NCols,  height=7);
        }
        #dev.new( width=4*NCols,  height=7);  #sizes for new windows don't work in RStudio... test with command line!
        
        OneScreenWidth = 0.9/NCols;
        figmatrixdata = c(0, OneScreenWidth, 0.05, 0.95);
        for (i in 2:NCols)
        {
            figmatrixdata = c(figmatrixdata, c((i-1)*OneScreenWidth, i*OneScreenWidth, 0.05, 0.95));
        }
        figmatrixdata = c(figmatrixdata, c(0.9, 1.0, 0.05, 0.95));
        figmatrix = matrix(figmatrixdata, nrow=NCols+1, ncol=4, byrow=TRUE);
        split.screen(figmatrix);
        
        LevelNames = levels(DataFr[,colInds[3]])
        for (i in 1:NCols)
        {
            screen(n=i, new=TRUE);
            datatmp = drop(datafrmeltcastmean[,,i])
            errtmp = drop(datafrmeltcasterror[,,i])
            
            if (PlotType=="bar")
            {
                CategoryNames = colnames(datatmp);
            }
            else
            {
                CategoryNames = rownames(datatmp);
            }
            if (i>1)
            {
                ColorTmp = qemegs.oneplot(datatmp, cmap, errorbar, errtmp, plottype, LevelNames[i], "", YLimVec, CategoryNames); 
            }
            else
            {
                ColorTmp = qemegs.oneplot(datatmp, cmap, errorbar, errtmp, plottype, LevelNames[i], YLabel, YLimVec, CategoryNames);
            }
        }
        title(outer=TRUE, adj=0.5, main = list(plottitle, cex=1.1, col="black", font=2), line = -2) 
        
        screen(n=NCols+1, new=TRUE);
        if (PlotType=="bar")
        {
            add_legend("topright", rownames(datatmp),  fill=ColorTmp,  cex=0.8, inset=c(0.02, 0.3))
        }
        else
        {
            add_legend("topright", colnames(datatmp),  fill=ColorTmp,  cex=0.8, inset=c(0.02, 0.3))
        }
        
        close.screen(all =TRUE);
        
        if (!debug)
        {
            dev.off();
        }
        #white <polygon> elements hide second axis and later below them - remove in QEMEGS
    }
    # 4 grouping factors
    else if (length(grouping)==4)
    {
        NCols = dim(datafrmeltcastmean)[3];
        NRows = dim(datafrmeltcastmean)[4];
        
        if (!debug)
        {
            if (R.Version()$os=="linux-gnu")
            {
                devSVG(svgfile, width=4*NCols,  height=4*NRows);
            }
            else
            {
                svg(svgfile, width=4*NCols,  height=4*NRows, onefile=TRUE);
            }
        }
        else
        {
            dev.new( width=6*NCols,  height=NRows*6); 
        }
         #sizes for new windows don't work in RStudio... test with command line!
        
        OneScreenWidth = 0.9/NCols;
        OneScreenHeight = 0.9/NRows;
        figmatrixdata = NULL
        for (j in 1:NRows)
        {
            for (i in 1:NCols)
            {
                figmatrixdata = c(figmatrixdata, c((i-1)*OneScreenWidth, i*OneScreenWidth, 0.1 + (NRows-j)*OneScreenHeight, 0.1 + (NRows-j+1)*OneScreenHeight));
            }
        }
        figmatrixdata = c(figmatrixdata, c(0.9, 1.0, 0.05, 0.95));  #left right bottom top!!!!!
        figmatrix = matrix(figmatrixdata, nrow=NCols*NCols+1, ncol=4, byrow=TRUE);
        screenvec = split.screen(figmatrix);
        
        LevelNamesCol = levels(DataFr[,colInds[3]])
        LevelNamesRow = levels(DataFr[,colInds[4]])
        countscreen = 1
        for (j in 1:NRows)
        {
            for (i in 1:NCols)
            {
                screen(screenvec[countscreen], new=TRUE);
                countscreen = countscreen + 1
                datatmp = drop(datafrmeltcastmean[,,i,j])
                errtmp = drop(datafrmeltcasterror[,,i,j])
                
                
                if (PlotType=="bar")
                {
                    CategoryNames = colnames(datatmp);
                }
                else
                {
                    CategoryNames = rownames(datatmp);
                }
                
                #skip category lables in upper rows
                if (j<NRows)
                {
                    for (k in 1:length(CategoryNames))
                    {
                        CategoryNames[k] = "";
                    }
                    
                    #leave a little margin for line plot, not for bar lot
                    if (PlotType=="bar")
                    {
                        par(mar=c(0,6,4,1));#bottom left top right
                    }
                    else
                    {
                        par(mar=c(1,6,4,1));#bottom left top right
                    }
                    
                    #skip axis lable in right columns
                    if (i>1)
                    {
                        ColorTmp = qemegs.oneplot(datatmp, cmap, errorbar, errtmp, plottype, ""              , ""              , YLimVec, CategoryNames); 
                    }
                    else
                    {
                        ColorTmp = qemegs.oneplot(datatmp, cmap, errorbar, errtmp, plottype, ""              , LevelNamesRow[j], YLimVec, CategoryNames);
                    }
                }
                #add category lables and subtitle in lower row
                else
                {
                    #leave a little margin for line plot, not for bar lot
                    if (PlotType=="bar")
                    {
                        par(mar=c(6,6,0,1)); #bottom left top right
                    }
                    else
                    {
                        par(mar=c(6,6,1,1));
                    }
                    #skip axis lable in right columns
                    if (i>1)
                    {
                        ColorTmp = qemegs.oneplot(datatmp, cmap, errorbar, errtmp, plottype, LevelNamesCol[i], ""              , YLimVec, CategoryNames); 
                    }
                    else
                    {
                        ColorTmp = qemegs.oneplot(datatmp, cmap, errorbar, errtmp, plottype, LevelNamesCol[i], LevelNamesRow[j], YLimVec, CategoryNames);
                    }
                }
            }
        }
        title(outer=TRUE, adj=0.5, main = list(plottitle, cex=1.1, col="black", font=2), line = -2) 
        
        screen(n=screenvec[length(screenvec)], new=TRUE);
        if (PlotType=="bar")
        {
            add_legend("topright", rownames(datatmp),  fill=ColorTmp,  cex=0.8, inset=c(0.02, 0.3))
        }
        else
        {
            add_legend("topright", colnames(datatmp),  fill=ColorTmp,  cex=0.8, inset=c(0.02, 0.3))
        }
        
        close.screen(all =TRUE);
        
        if (!debug)
        {
            dev.off();
        }
    }
    # more grouping factors
    else 
    {
      print("too many grouping factors!")
      return()
    }
    
    
    
}


qemegs.oneplot <- function(Data, ColorMap, ErrorBarStatus, DataErr, PlotType, SubTitle, YLabel, YLimVec, CategoryNames, ...)
{
    if (PlotType=="bar")
    {
        if (ColorMap=="Terrain")
        {
            Colors =terrain.colors(dim(Data)[1]);
        }
        else if(ColorMap=="Gray")
        {
            Colors =gray.colors(dim(Data)[1]);
        }
        else if(ColorMap=="Cm")
        {
            Colors =cm.colors(dim(Data)[1]);
        }
        else if(ColorMap=="Topo")
        {
            Colors =topo.colors(dim(Data)[1]);
        }
        else if(ColorMap=="Heat")
        {
            Colors =heat.colors(dim(Data)[1]);
        }
        else if(ColorMap=="Rainbow")
        {
            Colors =rainbow(dim(Data)[1]);
        }
    }
    else  if (PlotType=="line")
    {
        if (ColorMap=="Terrain")
        {
            Colors =terrain.colors(dim(Data)[2]);
        }
        else if(ColorMap=="Gray")
        {
            Colors =gray.colors(dim(Data)[2]);
        }
        else if(ColorMap=="Cm")
        {
            Colors =cm.colors(dim(Data)[2]);
        }
        else if(ColorMap=="Topo")
        {
            Colors =topo.colors(dim(Data)[2]);
        }
        else if(ColorMap=="Heat")
        {
            Colors =heat.colors(dim(Data)[2]);
        }
        else if(ColorMap=="Rainbow")
        {
            Colors =rainbow(dim(Data)[2]);
        }
    }
        
        
    if (PlotType=="bar")
    {
            xpos = barplot(Data, main="", sub=SubTitle, beside=TRUE, col=Colors, legend.text=FALSE,
                           names.arg=CategoryNames, ylab=YLabel, ylim=YLimVec,  ... );       
            
        if (ErrorBarStatus)
        {
            segments(xpos, Data+ DataErr, xpos, Data - DataErr);
            errorbartickwidth = (max(xpos) - min(xpos))*0.03;
            segments(xpos- errorbartickwidth, Data+ DataErr, xpos+errorbartickwidth, Data+ DataErr);
            segments(xpos- errorbartickwidth, Data- DataErr, xpos+errorbartickwidth, Data- DataErr);
        }
    }
    else  if (PlotType=="line")
    {
        XEnlarge = (dim(Data)[1]) * 0.05;
        XLimVec = c(1-XEnlarge , dim(Data)[1] + XEnlarge);
        plot(0,0, main="", sub=SubTitle,  type="b", xlab="", ylab=YLabel, xlim=XLimVec, ylim=YLimVec, xaxt="n", cex.sub = 1.2);
    
        
        xpos = 1:dim(Data)[1];
        jitterwidth = (max(xpos) - min(xpos))*0.04;
        xjitter = seq(-(dim(Data)[2]-1)/2*jitterwidth, (dim(Data)[2]-1)/2*jitterwidth, length=dim(Data)[2])  
        for (i in 1:dim(Data)[2])
        { 
          lines( xpos+ xjitter[i] , Data[,i], col = Colors[i],type = 'b'); 
        }
        axis(side=1,xpos, labels=CategoryNames, cex.axis=1.2);
        if (ErrorBarStatus)
        {
            errorbartickwidth = (max(xpos) - min(xpos))*0.03;
            for (i in 1:dim(Data)[2])
            { 
                segments(xpos + xjitter[i], Data[,i] + DataErr[,i], xpos + xjitter[i] , Data[,i] - DataErr[,i], col = Colors[i]);
                segments(xpos + xjitter[i] -errorbartickwidth, Data[,i] + DataErr[,i], xpos + xjitter[i] +errorbartickwidth, Data[,i] + DataErr[,i], col = Colors[i]);
                segments(xpos + xjitter[i] -errorbartickwidth, Data[,i] - DataErr[,i], xpos + xjitter[i] +errorbartickwidth, Data[,i] - DataErr[,i], col = Colors[i]);
            }
        }
    }
    
    return(Colors);
}


add_legend <- function(...) {
    opar <- par(mar=c(0, 0, 0, 0), new=TRUE)
    on.exit(par(opar))
    plot(0, 0, type='n', bty='n', xaxt='n', yaxt='n')
    legend(...)
}
