Useful blog "getting genetics done"

http://gettinggeneticsdone.blogspot.com/

latexdiff - mark up differences between latex files

Determine and mark up significant differences between latex files.

http://www.ctan.org/tex-archive/support/latexdiff

This perl script will generate a latex file with changes marked up (similar to Word) from two latex files

latexdiff-so fileold.tex filenew.tex > dif.tex

ps2pdf to all ps files in current directory

for x in $(ls *.ps); do echo $x; ps2pdf $x $(echo $x|sed -e 's/ps$/pdf/');done

or more efficiently according to Yarko's comment

for x in *.ps; do
echo $x;
ps2pdf $x ${x/%ps/pdf};
done

mixed effects model in R

Book
http://lme4.r-forge.r-project.org/book/

Slides
http://www.stat.wisc.edu/~bates/PotsdamGLMM/LMMD.pdf

awk: verify same number of columns

## prints number of columns if it changed from previous line

> awk -F"\t" 'NF != oldnf { print NF; oldnf = NF }' filename |head

|head added just in case file is very long

## using a script

> awk -f script.awk filename | head

script.awk:

BEGIN {
 FS="\t"
 oldnf = 1
}
{
  if(NF != oldnf) 
  {
 print NF
 oldnf = NF
  }
}

Other examples in 
http://sparky.rice.edu/~hartigan/awk.html

merge eps files into one

convert *.eps all.pdf

running R64 in aquamacs

http://www.matthewckeller.com/html/aquamacs.html

USING 64-BIT R FROM WITHIN AQUAMACS

1) Aquamacs is terrific because it comes bundled with about everything you need to 
start using R with it. No need to separately download ESS (emacs speaks statistics) 
or mess with the .ess-lisp files. Download Aquamacs from http://aquamacs.org/ and 
install it.

2) Now any file with an “.R” extension that is opened with aquamacs will 
automatically have syntax highlighting AND will allow you to start 64-bit R. I 
highly recommend checking out this page, http://ess.r-project.org/, and downloading 
the ESS reference card therein if you are new. Emacs has a bit of a learning curve, 
but you’ll begin loving its capabilities as an aid to programming very soon.

3) To start 64-bit R from within aquamacs hit Control-u then Alt-x then “R” then 
. Aquamacs will ask you for “Starting Args”. Here’s where you tell it you 
want to run 64-bit R. Type “--arch=x86_64” then . A new window will open up 
that is the R session.

4) That’s almost it. Just one more issue: graphing will crash (something to do with 
the new R and X11). Aquamacs has X11 as the default (null) graphics device. To 
change this, type in R: “options(device=”quartz”)”. Now your null graphics device is 
quartz and it should all work out. If you really want to use X11 to write graphics, 
you can use x11(type="Xlib")' for each new graph. What I do is to set 
“options(device=”quartz”) within my Rprofile.site file.

trim white spaces

check ?gsub

pvalue plot vs uniform; FDR

qqunif  = function(p,BH=T,...)
{
  nn = length(p)
  xx =  -log10((1:nn)/(nn+1))
  plot( xx,  -sort(log10(p)),
     xlab=expression(Expected~~-log[10](italic(p))),
        ylab=expression(Observed~~-log[10](italic(p))),
       ... )
  abline(0,1,col='gray')
  if(BH)
    {
      lines( xx, -log10(0.05*(1:nn)/(nn+1)), col='red')
      lines( xx, -log10(0.10*(1:nn)/(nn+1)), col='orange')
      lines( xx, -log10(0.25*(1:nn)/(nn+1)), col='yellow')
      legend('bottomright', c("FDR = 0.05","FDR = 0.10","FDR = 0.25"),
             col=c('red','orange','yellow'),lty=1, cex=1.4)
    }
}

imagesc in CDAT

The imagesc function in MATLAB creates plots that have color bars scaled to a user-specified range. This allows you to simply create a whole series of plots that can be easily intercompared.

The following function imitates imagesc:

    import vcs
    def imagesc(x, mini, maxi):
        t = vcs.init()
        a = t.getboxfill('quick')
        a.level_1 = mini
        a.level_2 = maxi
        t.plot(x, a)

x is the input array and mini and maxi define the range in values to plot. The function establishes a VCS canvas, imports the "quick" boxfill graphics method, sets the level range to mini and maxi, and plots x using this customized graphics method. 


other tips for CDAT/Python check http://www.johnny-lin.com/cdat_tips/

likelihood ratio test

> logLik.test
function(fit0,fit)
{
  l0 = logLik(fit0)
  l1 = logLik(fit)
  lratio = abs(as.numeric(l0)-as.numeric(l1))
  df = abs( attr(l1,'df') - attr(l0,'df') )
  pval = 1-pchisq( 2*lratio, df)
  return(pval)
}

log likelihood from lm, glm, polr

> fit = lm(y ~ x)
> logLik(fit)
'log Lik.' -36.4 (df=4)
> attr(logLik(fit),'df')
[1] 4

R - How to save graphs in files

" R-wiki on how to save graphs in file "

how to embed a pdf or ppt on a blog

get the code from

https://docs.google.com/viewer

plots boxplot, stripchart and mean

## CAN BE REPLACED BY qplot in ggplot2 package

qplot(x,y,geom=c("boxplot","jitter"))


## plots boxplot, stripchart and mean

boxjitter = function(x,y=NULL,border='gray',jitter=.1,vertical=T,...)
{
  if(!is.null(y))
  {
    boxplot(y~x,border=border,ylim = range(y,na.rm=T),outline=F,...)
    stripchart(y~x,vertical=vertical,method="jitter",jitter=jitter,add=T)
    meanvec = tapply(y,x,mean.na)
    points(1:length(meanvec),meanvec,pch="M",col='blue')
  } else
  {
    x = x[!is.na(x)]
    boxplot(x,border=border,ylim = range(x),outline=F,...)
    stripchart(x,vertical=vertical,method="jitter",jitter=jitter,add=T)
    points(1,mean.na(x),pch="M",col='blue')
  }
}

calculate number of minor alleles from A G format

num.minal = function(a1,a2)
{
  aall = c(a1,a2)
  tab = table(aall)
  miname = names(tab)[tab==min(tab)]
  apply(cbind(a1==miname,a2==miname),1,sum)  
}

plink coding --recodeAD

plink --file data --recodeAD

which, assuming C is the minor allele, will recode genotypes as follows:

     SNP       SNP_A ,  SNP_HET
     ---       -----    -----
     A A   ->    0   ,   0
     A C   ->    1   ,   1
     C C   ->    2   ,   0
     0 0   ->   NA   ,  NA

R table, percentages and chi square test like in stata

table.all = function(...)
{
  args = list(...)
  dnn = names(args)
  x = args[[1]]
  y = args[[2]]
  print(table(x,y,dnn=dnn))
  print("percentage by row")
  print( round( prop.table(table(x,y,dnn=dnn),1) * 100 ) )
  print("percentage by col")
  print( round( prop.table(table(x,y,dnn=dnn),2) * 100 ) )
  print(chisq.test(table(x,y,dnn=dnn)) )
}

order data frames in R

say data has read and prog as colums

data[order(read, prog), ]

this will order first by read then by prog

data[order(-read), ]

- reverses the order

** there can be conflicts if read is a variable outside of the data frame. Use data[order(data[,1]),] assuming read is column 1

from: http://www.ats.ucla.edu/stat/R/faq/sort.htm

list names from data that contain charname

## list names from data frames that contain charname

grepnames = function(charname,data)
  {
    names(data)[grep(charname,names(data))]
  }

power calculation in SAS

proc power; 
  twosamplefreq test=pchi 
  groupproportions = (.5034  .6459) 
  power = .
  groupns = ( 1494 373)
  alpha = 0.00005;
run;

proc power; 
      twosamplesurvival test=logrank 
         gexphs = 0.6733 | 0.4155
         accrualtime = .00001
   follouptime = 1 
         groupns = 200 | 200 
   grouplossexphazards = (0 0)
         power = .; 
proc power; 
      twosamplesurvival test=logrank 
         curve("Control")   = (1):(.45) 
         curve("Treatment") = (1):(.66) 
         groupsurvival = "Control" | "Treatment" 
   GROUPLOSSEXPHAZARDS = (.2231 .2231)
   /*grouplossexphazards = (0,0)*/
         accrualtime = .00001
         follouptime = 1 
         groupns = 120 | 120 
         power = .; 
run;

proc power; 
      twosamplesurvival test=logrank 
         curve("Control")   = (1):(.51) 
         curve("Treatment") = (1):(.66) 
         groupsurvival = "Control" | "Treatment" 
   GROUPLOSSEXPHAZARDS = (.2231 .2231)
   /*grouplossexphazards = (0,0)*/
         accrualtime = .00001
         follouptime = 1 
         groupns = 100 | 100 
         power = .; 
run;

proc power; 
      twosamplesurvival test=logrank 
         curve("Control")   = (1):(.45) 
         curve("Treatment") = (1):(.70) 
         groupsurvival = "Control" | "Treatment" 
   GROUPLOSSEXPHAZARDS = (.2231 .2231)
   /*grouplossexphazards = (0,0)*/
         accrualtime = .00001
         follouptime = 1 
         groupns = 200 | 200 
         power = .; 
run;
 

how to get rid of emac's automatic switching of _ into ->

Solution by Paul R. (thanks!)

M-x
set-variable 
ess-S-assign 
"-"  

That's it.

File Management in Python 
( Page 1 of 5 )

File management is a basic function, and an important part of many applications. Python makes file management surprisingly easy, especially when compared to other languages. Peyton McCullough explains the basics.
 
Introduction
 
The game you played yesterday uses files to store game saves. The order you placed yesterday was saved in a file. That report you typed up this morning was, obviously, stored in a file as well.
 
File management is an important part of many applications written in nearly every language. Python is no exception to this. In this article, we will explore the task of manipulating files using several modules. We'll read, write to, append and do other strange things to files. Let's get started.
 
Reading and Writing
 
The most basic tasks involved in file manipulation are reading data from files and writing data to files. This is a very simple task to learn. Let's open a file for writing: 

fileHandle = open ( 'test.txt', 'w' ) 

The "w" indicates that we will be writing to the file, and the rest is pretty simple to understand. The next step is to write data to the file: 

fileHandle.write ( 'This is a test.\nReally, it is.' ) 

This will write the string "This is a test." to the file's first line and "Really, it is." to the file's second line. Finally, we need to clean up after ourselves and close the file: 

fileHandle.close() 

As you can see, it's very easy, especially with Python's object orientation. Note that when you use the "w" mode to write to the file again, all its contents will be deleted. To get past this, use the "a" mode to append data to a file, adding data to the bottom: 

fileHandle = open ( 'test.txt', 'a' )
fileHandle.write ( '\n\n\nBottom line.' )
fileHandle.close() 

Now let's read our file and display the contents: 

fileHandle = open ( 'test.txt' )
print fileHandle.read()
fileHandle.close() 

This will read the entire file and print the data within it. We can also read a single line in the file: 

fileHandle = open ( 'test.txt' )
print fileHandle.readline() # "This is a test."
fileHandle.close() 

It is also possible to store the lines of a file into a list: 

fileHandle = open ( 'test.txt' )
fileList = fileHandle.readlines()
for fileLine in fileList:
   print '>>', fileLine
fileHandle.close() 

When reading a file, Python's place in the file will be remembered, illustrated in this example: 

fileHandle = open ( 'test.txt' )
garbage = fileHandle.readline()
fileHandle.readline() # "Really, it is."
fileHandle.close() 

Only the second line is displayed. We can, however, get past this by telling Python to resume reading from a different position: 

fileHandle = open ( 'test.txt' )
garbage = fileHandle.readline()
fileHandle.seek ( 0 )
print fileHandle.readline() # "This is a test."
fileHandle.close() 

In the above example, we tell Python to continue reading from the first byte in the file. Thus, the first line is printed. We can also request Python's place within the file: 

fileHandle = open ( 'test.txt' )
print fileHandle.readline() # "This is a test."
print fileHandle.tell() # "17"
print fileHandle.readline() # "Really, it is." 

It is also possible to read the file a few bytes at a time: 

fileHandle = open ( 'test.txt' )
print fileHandle.read ( 1 ) # "T"
fileHandle.seek ( 4 )
print FileHandle.read ( 1 ) # "T" 

When working with Windows and Macintosh, sometimes you are required to read and write files in binary mode, such as images or executional files. To do this, simply append "b" to the file mode: 

fileHandle = open ( 'testBinary.txt', 'wb' )
fileHandle.write ( 'There is no spoon.' )
fileHandle.close()
 
fileHandle = open ( 'testBinary.txt', 'rb' )
print fileHandle.read()
fileHandle.close() 


#write to file
f = open('c:/temp/workfile','w')
f.write('hallo')
f.close
 
#open from file
f = open('file.txt')
 
#read linewise
for i in f:
    print i
 
#read characterwise
for i in f.read():
    print i

myfunctions.r

## logit/expit
logit = function(x,b=1) log(x/(b-x))
expit = function(x,b=1) b/(1+exp(-x))

## matlab equivalent of repmat in R
repmat = function(X,m,n){
  ## tiles matrix X, m by n times
  mx = dim(X)[1]
  nx = dim(X)[2]
  matrix(t(matrix(X,mx,nx*n)),mx*m,nx*n,byrow=T)
}

acfi = function(x, lag.max = NULL,
         type = c("correlation", "covariance", "partial"),
         plot = TRUE,  demean = TRUE,...)
{
acf(x, lag.max = lag.max,
         type = type,
         plot = plot, na.action = na.pass, demean = TRUE, ...)
}

"%&%" <- function(a, b) paste(a, b, sep="")


rangi = function(x){
 range(x, na.rm = T)
}

cori = function(x,y=NULL){
cor(x, y = y, use = "complete.obs")
}


## myoptim
## minimizes fun using optim with nelder mead
## algorithm. 

myoptim = function(param,fun,epsi=1e-5,hess=F){
  dife = 10
  eini = fun(param)
  while(dife>epsi){
    res=optim(param,fun,hessian=hess)
    dife = abs(eini-res$value)
    param=res$par
    eini=res$value
  }
  res
}

##======================
## mean squared error sum((x-y)^2)/n
mse <- function(x,y) {
  n <- length(x)
  sum((x-y)^2)/n}
                                                                                
rmse <- function(x,y) sqrt(mse(x,y))
##======================
regress <- function(form,plot=F){
  fit=lm(form)
  fitsum=summary(fit);
  print(fitsum$coeff);
  print(c("R squared","residual se"))
  print(round(c(fitsum$r.squared,rmse(fit$res,0)),3))
  }

## Great circle distance from library fields, default is km.
gcdist <- function (loc1, loc2, miles = FALSE, R = NULL) 
{
    if (is.null(R)) {
        if (miles) 
            R <- 3963.34
        else R <- 6378.388
    }
    if (missing(loc2)) 
        loc2 <- loc1
    coslat1 <- cos((loc1[, 2] * pi)/180)
    sinlat1 <- sin((loc1[, 2] * pi)/180)
    coslon1 <- cos((loc1[, 1] * pi)/180)
    sinlon1 <- sin((loc1[, 1] * pi)/180)
    coslat2 <- cos((loc2[, 2] * pi)/180)
    sinlat2 <- sin((loc2[, 2] * pi)/180)
    coslon2 <- cos((loc2[, 1] * pi)/180)
    sinlon2 <- sin((loc2[, 1] * pi)/180)
    pp <- cbind(coslat1 * coslon1, coslat1 * sinlon1, sinlat1) %*% 
        t(cbind(coslat2 * coslon2, coslat2 * sinlon2, sinlat2))
    R * acos(ifelse(pp > 1, 1, pp))
}
## cord distance, default is km.
## loc1 = [lon, lat ]
cordist <- function (loc1, loc2, miles = FALSE, R = NULL) 
{
    if (is.null(R)) {
        if (miles) 
            R <- 3963.34
        else R <- 6378.388
    }
    if (missing(loc2)) 
        loc2 <- loc1
    coslat1 <- cos((loc1[, 2] * pi)/180)
    sinlat1 <- sin((loc1[, 2] * pi)/180)
    coslon1 <- cos((loc1[, 1] * pi)/180)
    sinlon1 <- sin((loc1[, 1] * pi)/180)
    coslat2 <- cos((loc2[, 2] * pi)/180)
    sinlat2 <- sin((loc2[, 2] * pi)/180)
    coslon2 <- cos((loc2[, 1] * pi)/180)
    sinlon2 <- sin((loc2[, 1] * pi)/180)
    pp <- cbind(coslat1 * coslon1, coslat1 * sinlon1, sinlat1) %*% 
        t(cbind(coslat2 * coslon2, coslat2 * sinlon2, sinlat2))
    R * sqrt( 2* (1-ifelse(pp > 1, 1, pp) ) )
}



## matern function with Handcock's parameterization
## hmatern
hmatern = function(u,range,smoo){
  if (is.vector(u)){
    names(u) <- NULL
    y = u
    y[]=1
  }
  if (is.matrix(u)) {
    dimnames(u) <- list(NULL, NULL)
    y = u
    y[,]=1
  }
  unonzero= u>0
  if(smoo<100) {
    urange = u*2*sqrt(smoo)/range
    y[unonzero] = (urange[unonzero]^smoo) * besselK(urange[unonzero],smoo)/
      (2^(smoo - 1))/ gamma(smoo)}
  else
    {y[unonzero] = exp(-u[unonzero]^2)}
  y
  }
##===========================================================
## simple useful functions
##===========================================================
suma <- function(x) {
  indnna = !is.na(x)
  if(sum(indnna)>0) res = sum(x[indnna]) else res=NA
  res
}
maxi <- function(x){
  indnna = !is.na(x)
  if(sum(indnna)>0) res=max(x[indnna]) else res=NA
  res}

mini <- function(x){
  indnna = !is.na(x)
  if(sum(indnna)>0) res=min(x[indnna]) else res=NA
  res}

meana <- function(x){
  indnna = !is.na(x)
  if(sum(indnna)>0) res=mean(x[indnna]) else res=NA
  res}

vari <- function(x){
  indnna = !is.na(x)
  if(sum(indnna)>0) res=var(x[indnna]) else res=NA
  res
}

longi <- function(data){length(data[,1])}

plotfit <- function(fit)
         {
           par( mfrow=c(2,2) )
           plot(fit)
           par( mfrow=c(1,1) )
         }

algunos <- function(x,cant=10) x[1:cant,]

## returns the row and column of an index matrix == 1
indmat <- function(x)
{
  n <- length(x[,1])
  tm1 <- seq(1,n)
  tm2 <- matrix(rep(1,n*n),n)
  rowmat <- diag(tm1)%*%tm2
  colmat <- tm2%*%diag(tm1)
  ro <- rowmat[x==1]
  co <- colmat[x==1]
  cbind(ro,co)
}

reducing margins in R plots

reduce margins
par(mar=c(2,2,1,1)) # default in R is c(5, 4, 4, 2) + 0.1
c(bottom, left, top, right)

confidence intervals ACF

under independence assumption: clim0 = qnorm((1+.95)/2)/sqrt(N)
ARIMA: at k-th lag: clim0 * sqrt( 1+2 \sum_1^k( acf(i)^2 ) ) 

rules of thumbs for ARIMA model selection

 Prof. R Nau's lecture notes

theoretical ACF of ARMA

ARMAacf(ar = numeric(0), ma = numeric(0), lag.max = r, pacf = FALSE)

plot in log scale with axis

plot(1:100, log = "y", yaxt = "n") # do not show y axis
axis(2, c(1,10,100)) # draw y axis with required labels 

Other useful R functions blogs and pages

Paul Johnson's R Page

Li Chen's R page

James Forester R Functions' Blog 

R blog

calculating angle from xy coordinates

## from north (-pi,pi) positive toward east
angle.north = function(xx,yy){sign(xx)*(yy<0)*pi + atan(xx/yy)}

## from east (0,2.pi) positive toward north
angle = function(xx,yy){atan(yy/xx) - sign(xx)*sign(yy)*pi*(xx<0)}

empirical distributions CDF

"empirical CDF's are not good at showing differences in the tails of a distribution"
Owen, A. B., Empirical Likelihood, 2001, Chapman & Hall/CRC, pg 10.

rsync excluding Rdata files

rsync -avz --exclude 'temper.*.Rdata' asthma aitken:~/donnpc/

plot with confidence intervals

## modified version of the code I found on the web.
## I forgot who wrote the origial version

plotCI <- function (x, y = NULL, uiw, liw = uiw, ylim=NULL,..., sfrac = 0.01)  {
 col2 = 'gray'
 if (is.list(x)) {
   y <- x$y
   x <- x$x
 }
 if (is.null(y)) {
   if (is.null(x))
     stop("both x and y NULL")
   y <- as.numeric(x)
   x <- seq(along = x)
 }
 ui <- y + uiw
 li <- y - liw
 if(is.null(ylim)) {
   plot(x, y, ylim = range(c(y, ui, li),na.rm=T), ...)} else
   {plot(x, y, ylim = ylim, ...)}
 smidge <- diff(par("usr")[1:2]) * sfrac
 segments(x, li, x, ui,col=col2)
 x2 <- c(x, x)
 ul <- c(li, ui)
 segments(x2 - smidge, ul, x2 + smidge, ul,col=col2)
 if(is.null(ylim))  ylim=range(c(y, ui, li),na.rm=T)
 points(x, y, ...)
 invisible(list(x = x, y = y))
}

image with gray scale

imagebw = function(x,y,z,...) image(x,y,z,col=gray((33:64)/64),...)

Mac Eye for the Linux Guy

Mac Eye for the Linux Guy

rsync excluding subdirectories contents

rsync -dvlptgoDz source.server:~/source.dir/ target.dir

R documentation with graphics

R documentation with graphics

specifying printer lpr

    lpr -#2 -sP dj thesis.txt
This command will create a symbolic link to the file thesis.txt in the spool 
directory for the printer named dj, where it would be processed by lpd. 
It would then print a second copy of thesis.txt.

bold greek letters latex

% use bm package
\usepackage{bm}

% inside math mode
\[  ... {\bm beta} ... \]

better than \boldmath because it works within math mode directly, 
no need to use in text mode.

logit - expit

%% maps (0,b) into (-inf,inf)
%% inverse of expit
%% usage: logit(x,b) or logit(x) 
%% b defaults to 1
function y = logit(x,b) 
  if nargin<2
    b = 1;
  end
  y = log( x ./ (b - x) );

%% maps (-inf, inf) into (0,b)
%% inverse of logit
%% usage: expit(x,b) or expit(x) 
%% b defaults to 1
function y = expit(x,b) 
  if nargin<2
    b = 1;
  end
  y =  b./(1+exp(-x)) ;

xcdroast

to copy files to CD in linux
> xcdroast

psnup

psnup -d -8 oldfile.ps newfile.ps
many postscript pages in one

R equivalent of repmat (matlab)

repmat = function(X,m,n){
##R equivalent of repmat (matlab)
mx = dim(X)[1]
nx = dim(X)[2]
matrix(t(matrix(X,mx,nx*n)),mx*m,nx*n,byrow=T)

ssh without password howto

 
1. Run "ssh-keygen -t rsa".  

2. Save the generated key in the default location.

3. Leave the passphrase blank.  

4. Change to the .ssh directory, and append the contents of id_rsa.pub
to the file "authorized_keys." If no authorized_keys file exists, you
can just run "cp id_rsa.pub authorized_keys".

Now you should be able to ssh to department machines without typing
your password. Note that if you're doing this from a non-department
Linux machine, you should do steps 1, 2, and 3 from the command prompt
of your machine, and step 4 from the command prompt of any department
Linux machine. If you're doing this from a department Linux machine,
do all steps on that machine.

- Elliot

R legacy help galton

## calling C functions from R
gcc -c kendall.c
R CMD SHLIB -o libkendall.so kendall.o
R
dyn.load('libkendall.so')
xd <- c(44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 45.2, 60.1)
yd <- c( 2.6,  3.1,  2.5,  5.0,  3.6,  4.0,  5.2,  2.8,  3.8)
.C("kendall",as.double(xd),as.double(yd),as.integer(9),
result=as.integer(1))$result

## header
rm(list=ls())
options(contrasts=c("contr.treatment","contr.poly"))
setwd("/mnt/win_c/my documents/stat/ammonia")

## plot function
plot(fun,-5, 5, ylim = c(-.2, 1))

eval.parent(expr) 
## example eval.parent(media <- 1), used inside
## functions to change global variables
media <<- 1 ## does the same

rep(x,n) replicates x n-times

mtext("texto",side=3,line=.5)

install.packages("fields","/aa/haky/R/Library")
library("fields",lib="/aa/haky/R/Library")
detach(package:fields)

setenv R_LIBS /aa/haky/R/Library (then no need to specify lib=...)

xmat <- matrix(NA,n,n)

#print
dev.print()
dev.copy2eps(file="filename.ps")
#
postscript(file="plot.ps",horizontal=F,heigth=9,width=7)
dev.off()

as.numeric(x<1)

# stop watch
tic<-Sys.time()
....do something here
toc<-Sys.time()
print(toc-tic)

paste(s1,s2,sep="")

system("mail -s something haky@galton.uchicago.edu < output")<

## strip an object's attributes:
attributes(x) <- NULL
x # now just a vector of length 6

kk[names(kk)=="site"] is same as kk$site

unlist(x)

debug(myfunction)
undebug(myfunction)

call R from emacs (ESS): alt-x R

par(mfrow=c(3,2))

legacy galton help

# synchronize file system
rsync -avz --exclude '*.eps' Penalty aitken:~/donnpc/ .
# dry run - shows only what would have been transferred
rsync -avzn --exclude '*.eps' Penalty aitken:~/donnpc/ . | less

# indent region mode dependent
M - C - \

# rigid tab region
C - x - TAB

# restart pcmcia, network
/etc/rc.d/init.d/pcmcia restart
/etc/init.d/network restart

# create symbolic links
ln -s /home/azubrow/CMAQ_v4.3/data/emis/M_36_EUSA ./emisfiles

# run identical command in all nodes
sshall uptime
sshall w

To run batch jobs in ravana
chmod u+x batch
./batch &

batch:
#!/bin/sh
/usr/bin/nohup nice -5 R BATCH mplefft.r output

Ctrl-z (suspend)
bg (sends process to background)
fg(bring process to foreground)

rm
rmdir
mkdir
mv
chmod 700 mydir (no access to others)

ftp
$ scp * ravana:~/

gunzip filename: unzip
gzip filename: zip
tar xvf filename.tar: extract tar
tar cvf filename.tar filename or directory: merges files

get rid of the ^M:  in vi type
:%s/(ctrl-v)(ctrl-m)//g
in emacs
M-x replace-regexp RET C-q C-m $ RET RET

chsh #changes shell in unix

mail -s "subject" email@galton <> mkdir work_dir
> cd work_dir
> Splus5 CHAPTER
> touch .Data/.Audit
> chmod 0 .Data/.Audit
emacs kk.s   alt-x S+6

cdms quick tutorial

> cdat This starts the CDAT shell, identical to the Python shell.
> import cdms Imports the CDMS module.
> f=cdms.open('wind_comps.nc') Opens the NetCDF file and assigns the file object to variable 'f'.
> f.listvariables() Shows the variable names in the file.
['u', 'v']
> u_wind=f('u') Reads the variable 'u' and assigns it to the Python variable 'u_wind'.
> u_wind.attributes Displays all the attributes pertaining to the variable 'u_wind'.
> v_wind=f('v') Reads the variable 'v' and assigns it to the Python variable 'v_wind'.
> wspd=(u_wind**2+v_wind**2)**0.5
Calculates the combined wind speed and assigns it to the Python variable 'wspd'.
> wspd.id='wspd' Sets the 'id' attribute on the 'wspd' variable.
> wspd.long_name='Wind speed' Sets the 'long_name' attribute on the 'wspd' variable.
> wspd.units='m s**-1' Sets the 'units' attribute on the 'wspd' variable.
> import vcs Imports the VCS (Visualisation Control System) module.
> p=vcs.init() Initialise a VCS canvas and assign Python variable 'p' to it.
> p.plot(wspd) Plot the wind speed variable on canvas 'p'.
> fout=cdms.open('output.nc', 'w') Open a file to write out the wind speed variable to.
> fout.write(wspd) Write the wind speed variable to the output file.

> fout.close() Close the output file.
> CTRL^D Press Control and D to close CDAT.