cndaqiang Web Linux DFT

R语言快速上手

2019-05-09
cndaqiang
R
RSS

R语言快速上手

参考

.Fortran()返回任何结果
卡巴, 科弗, 高涛, 等. R 语 言 实 战[J]. 2013.
dnorm, pnorm, qnorm与rnorm的区别
R语言的流程控制
R语言入门
R 语言–接收命令行参数
R里面有没有函数可以输出到已有文件末尾?
快速傅里叶变换 - MATLAB fft- MathWorks 中国
【R语言数据导出txt】 write.table 函数用法

说明

  • 运行脚本方法
    1) R CMD [option] [xxx.R] [输出]
    2) 与python和bash一样,在脚本里指明程序,添加可执行权限
    #!/usr/bin/env Rscript
    age <- c(1,2,3,4,5) 
    weight <- c(2,3,4,5,6) 
    plot(age,weight)
    
  • 注释用#
  • 同行书写,用;代表换行
  • 变量名区分大小写,字母数字下划线
  • 变量无需定义,随用随赋值,新创建,赋值为新数据=删除旧变量并重新定义赋值
  • R的数组和矩阵的元素类型不用相同
  • 变量赋值,教程中是b <- 5,在R version 3.3.1可以直接b=5;b=c(1,2,3)
  • 删除变量rm(变量名)
  • 使用[” “]提取元素,有,双引号
  • 在脚本或程序中多次执行plot,每次的图会变成一个pdf的一页,程序执行结束后,有个总图的pdf保存在当前目录
  • 使用source像模块一样调用
  • 当前目录若存在.RData .Rhistory文件,好像会影响脚本运行结果

安装程序包

安装R语言

  • 下载/编译 安装清华-mirrors
    可以下载安装包,也可以源码编译
  • homebrew安装
    #MacOS也可以
    brew install R
    
  • anaconda也可以安装
    cndaqiang@girl:~/work/dfft-yimin/out$ conda install R
    

CRAN 进入R后安装,无需root

install.packages("inline")

指定清华镜像

options(repos=structure(c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))) 

Bioconductor

Bioconductor是大部分的生信相关的包存放的地方,网址https://www.bioconductor.org/ 安装使用代码如下

## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("package.name")#将package.name替换为自己的要安装的包的名字即可。

GitHub

install.packages("devtools")
library(devtools)
install_github("user.name/package.name")#其中user.name是github作者的用户名。

查看程序包帮助

help(package = "sxtTools")
#如果想看某个具体的函数的信息,可以使用下面的代码:
?installBioc

数据类型

R的基本数据元素有数值,字符,逻辑型

基本数据元素

直接赋值,和数组就个c的区别

> a=5.0     #数值型
> a
[1] 5
> a="hello" #字符类型
> a
[1] "hello"
> a=T       #逻辑型
> a
[1] TRUE

数据类型转换

  • as.numeric(a) 转数值
  • as.character(b) 将数字b转化成字符格式
  • format(1.2345,digits = 3) 将数字转成字符格式,并设置有效位数为3
  • as.matrix(a)转矩阵
  • data.frame(Ave) 转数据框
  • as.numeric(as.vector(unlist(strsplit(Dipole,",")))) 组合操作,拆分字符串
  • options(digits=N),设置数据精度,屏幕显示的小数点位数

数组

数组元素可以是数值,字符,逻辑型混合在一起,不必相同类型

定义c

> a=4
> a
[1] 4
> a <- c(2,3)
> a
[1] 2 3

其他定义方式

> 1:4
[1] 1 2 3 4

调用[ ]

a[n]

矩阵

矩阵是由向量构成,因此元素可以是数值,字符,逻辑型混合在一起,不必相同类型

定义matrix

matrix(向量,[nrow=行数,ncol=列数][,byrow=T/F][,dimname=list(行名矩阵,列名矩阵)]  )

提取向量中的前nrow*ncol个元素填充矩阵
若不设置行数,默认为列向量
byrow=T向量按行填充,F按列填充, 示例

> matrix(c(1,2,"hello",3),nrow=2,ncol=1)
     [,1]
[1,] "1"
[2,] "2"
> matrix(c(1,2,"hello",3),nrow=4,ncol=1)
     [,1]
[1,] "1"
[2,] "2"
[3,] "hello"
[4,] "3"
> matrix(1:4,nrow=2,ncol=2,byrow=F,dimname=list(c("hang1","hang2"),c("lie1","lie2")))
      lie1 lie2
hang1    1    3
hang2    2    4

调用[ ]

> b=matrix(1:4,nrow=2,ncol=2,byrow=F,dimname=list(c("hang1","hang2"),c("lie1","lie2")))
> b[1,1]
[1] 1
> b["hang1","lie2"]
[1] 3
> aa
     [,1] [,2] [,3] [,4]
[1,]    1    3    5    7
[2,]    2    4    6    8
> aa[1,2:3]
[1] 3 5

也支持像matlab一样,条件调用

> x=1:10
> x[x>5]
[1]  6  7  8  9 10

注意,a:b/2其实是c(a:b)/2, 应该写 a:(b/2),在提取矩阵时用, mydata[1:(n/2),2]

数组array

矩阵只能是二维的,数组与矩阵一样,区别是可以是高维的

定义array

示例

 b=array(1:8,c(2,2,2),dimname=list(c("hang1","hang2"),c("lie1","lie2"),c("matrix1","matrix2")))
 > b
, , matrix1

      lie1 lie2
hang1    1    3
hang2    2    4

, , matrix2

      lie1 lie2
hang1    5    7
hang2    6    8

特定数组

x=1:10 !等差数组
> seq(1,5,2)
[1] 1 3 5

提取[ ]

> b["hang2","lie1","matrix2"]
[1] 6

数据框

是将多个变量保存到一个数据框中
可以根据变量名从数据框中提取数据

定义data.frame(col1,col2,col3...)

每列元素要一样多

> col1=1:5
> col2=2:6
> col3=c("q","w","e","r","t")
> data.frame(col1,col2,col3)
  col1 col2 col3
1    1    2    q
2    2    3    w
3    3    4    e
4    4    5    r
5    5    6    t
#指定行名称
> data.frame(col1,col2,col3,row.names=col3)
  col1 col2 col3
q    1    2    q
w    2    3    w
e    3    4    e
r    4    5    r
t    5    6    t
#后期修改列名
dipoledate=data.frame(result$getDipolePython())
names(dipoledate)=c("Time","x","y","z")

调用[ ]

个人觉得$挺好用的

> b=data.frame(col1,col2,col3)
> b["col2"]
  col2
1    2
2    3
3    4
4    5
5    6
> b[1:3,2:3] #和matlab一样呀,[  [行],列]内只有一个数组则从列提取
  col2 col3
1    2    q
2    3    w
3    4    e
> b$col1
[1] 1 2 3 4 5
> summary(b$col1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
      1       2       3       3       4       5

也可以通过attach将数据框加到R的搜索路径中,就可以直接用

> rm(col1)
> col1
Error: object 'col1' not found
> attach(b)     #加入搜索路径
The following objects are masked _by_ .GlobalEnv:

    col2, col3

> col1
[1] 1 2 3 4 5
> detach(b)     #移除搜索路径
> col1
Error: object 'col1' not found

也可以使用with(数据框,{命令})

> with(b,{
+ col1
+ })
[1] 1 2 3 4 5

因子

因子可以理解为对已有数据的分类
一般含有字符型的数组,默认就是factor因子

> sex=factor(c("男","男","女","男"))
> sex
[1]    
Levels:  

指定设置分类顺序levels=c(...)

> sex=factor(c("男","男","女","男"),levels=c("男","女"))
> sex
[1]    
Levels:  

分类后,并不影响原数据,可以更好的用于数据分析

> sex[1]
[1] 
Levels:  

常将分类后的向量组成数据框,用于进一步分析

> sex
[1]                Unkown
Levels: Unkown  
> num
[1] 12 14 12 15
Levels: 12 14 15
> name
[1] "Li" "Ya" "CQ" "LL"
> data=data.frame(name,sex,num)
> data
  name    sex num
1   Li       12
2   Ya       14
3   CQ       12
4   LL Unkown  15
> str(data)
'data.frame':   4 obs. of  3 variables:
 $ name: Factor w/ 4 levels "CQ","Li","LL",..: 2 4 1 3
 $ sex : Factor w/ 3 levels "Unkown","女",..: 3 3 2 1
 $ num : Factor w/ 3 levels "12","14","15": 1 2 1 3
> summary(data)
 name       sex    num
 CQ:1   Unkown:1   12:2
 Li:1       :1   14:1
 LL:1       :2   15:1
 Ya:1

列表

就像一个垃圾站,什么都能往里丢

定义

list(keywords1=...,keywords2=a,b,c...,keywords3....)
> list(one=num,two=name,sex)
$one
[1] 12 14 12 15
Levels: 12 14 15

$two
[1] "Li" "Ya" "CQ" "LL"

[[3]]
[1]                Unkown
Levels: Unkown  

调用listname[" keywords"]

> lista["two"]
$two
[1] "Li" "Ya" "CQ" "LL"
> lista$two
[1] "Li" "Ya" "CQ" "LL"

其他

对变量注释names

> a=c("xiaoming",13)
> a
[1] "xiaoming" "13"
> names(a)[1]="Name"
> names(a)[2]="age"
> a
      Name        age
"xiaoming"       "13"

变量(数据对象)操作

class(变量名)查看变量类型

# 在拿到数据之前第一步要做的就是概览数据
head(df,n=10)#查看数据集的前十条记录
str(df)#查看数据集结构
tail(df)#查看数据集最后部分
sapply(df)#查看数据集类型及名称
summary(df)#查看数据集各变量描述统计

函 数 用 途

length(object) 显示对象中元素/成分的数量
dim(object) 显示某个对象的维度
str(object) 显示某个对象的结构
class(object) 显示某个对象的类或类型
mode(object) 显示某个对象的模式
names(object) 显示某对象中各成分的名称
c(object, object,) 将对象合并入一个向量 

cbind(object, object, ) 按列合并对象
rbind(object, object, ) 按行合并对象
Object 输出某个对象
head(object) 列出某个对象的开始部分
tail(object) 列出某个对象的最后部分
ls() 显示当前的对象列表
rm(object, object, ) 删除一个或更多个对象语句rm(list = ls())
将删除当前工作环境中的几乎所有对象*
newobject <- edit(object) 编辑对象并另存为newobject
fix(object) 直接编辑对象调用vi编辑

变量(数据对象)处理函数

字符

函数grep()、 sub()和strsplit() 搜索某个文本字符串( fixed=TRUE)
搜索正则表达式( fixed=FALSE,默认值为FALSE)
正则表达式可参考[gitbook文章重置] 正则表达式 grep sed awk

转义字符串

要使用特殊字符时,要在前面加\\,如

name=unlist(strsplit(as.character(mydata[id,1]),"\\|"))[3]

其他处理函数

矩阵操作

转置t(a)

> a=matrix(1:10,nrow=2)
> t(a)
     [,1] [,2]
[1,]    1    2
[2,]    3    4
[3,]    5    6
[4,]    7    8
[5,]    9   10

整合重构数据,暂略

数学函数

和matlab差不多啊

x=1:10
plot(x,sin(x))

也可以画x图

数学运算

abs, 取整, cos, sin, acos,asin,三角函数,ln=log(x,base=exp(1))

统计函数

mean,sum

概率函数

在R中,概率函数形如:

[dpqr][分布名称](参数表)

[dpqr]表示其所指分布的某一方面:

  • d = 密度函数( density)!就是统计函数值
  • p = 分布函数( distribution function) !统计函数的积分
  • q = 分位数函数( quantile function)
  • r = 生成随机数(随机偏差)!生成符合分布的一组数x,x符合这一分布函数

以正态分布norm为例
其他分布名称

随机数种子

R每次计算伪随机数都会使用不同的随机数种子,我们可以指定随机数种子
这样可以重复出(自己/别人)的计算结果,如下
**注意:种子只生效一次,若仍使用该种子,要继续set.seed() **

> rnorm(10)
 [1]  1.2240818  0.3598138  0.4007715  0.1106827 -0.5558411  1.7869131
 [7]  0.4978505 -1.9666172  0.7013559 -0.4727914
> rnorm(10)
 [1] -1.0678237 -0.2179749 -1.0260044 -0.7288912 -0.6250393 -1.6866933
 [7]  0.8377870  0.1533731 -1.1381369  1.2538149
> set.seed(123)
> rnorm(10)
 [1] -0.56047565 -0.23017749  1.55870831  0.07050839  0.12928774  1.71506499
 [7]  0.46091621 -1.26506123 -0.68685285 -0.44566197
> set.seed(123)
> rnorm(10)
 [1] -0.56047565 -0.23017749  1.55870831  0.07050839  0.12928774  1.71506499
 [7]  0.46091621 -1.26506123 -0.68685285 -0.44566197

多元正态数据mvrnorm

将变量(矩阵,数组,数据框…)作为函数的输入参数

直接置入

> x=1:10
> sin(x)
 [1]  0.8414710  0.9092974  0.1411200 -0.7568025 -0.9589243 -0.2794155
 [7]  0.6569866  0.9893582  0.4121185 -0.5440211
> xx=matrix(1:10,nrow=2)
> sin(xx)
          [,1]       [,2]       [,3]      [,4]       [,5]
[1,] 0.8414710  0.1411200 -0.9589243 0.6569866  0.4121185
[2,] 0.9092974 -0.7568025 -0.2794155 0.9893582 -0.5440211
> mean(xx)
[1] 5.5

apply (数组,矩阵)某一维度

apply(变量名,维度,函数名)

> apply(xx,1,sin)
           [,1]       [,2]
[1,]  0.8414710  0.9092974
[2,]  0.1411200 -0.7568025
[3,] -0.9589243 -0.2794155
[4,]  0.6569866  0.9893582
[5,]  0.4121185 -0.5440211
> apply(xx,1,mean)
[1] 5 6

lapply,sapply (列表list)的某一维度,暂略

自编函数

判断

if

}后跟新的else

if (条件)
{ 
执行
}else if
{
执行    
}else{
执行 
}

示例

> if(a==10)
+ {
+ a=5
+ }
> if(a==5)
+ {
+ a=5
+ }else if(a==6)
+ {
+ a=6
+ }else
+ {a=10
+ }
> a
[1] 5

swith

swith(表达式,case1,case2,cse3)
#表达式为整数时,返回表达式对应的case值
#表达式为字符串时,case必须有名称,见下例,返回case值
#非整数非字符串,系统会转换成整数

好像就只能这么用

> a=2
> switch(a,"ming","hong","li")
[1] "hong"
> switch("a",a="ming",b="hong",c="li")
[1] "ming"

与或非

&、|、!、&&、||、xor

循环

for ( i in seq)

若省略{ }只循环for下面的一句,类似c

> for (i in 1:10)
+ {
+ y[i]=i
+ }

while (条件)

与for一样

> > while (i>5)
+ {
+ i=i-2
+ y=i
+ }

循环控制break和next

自定义函数

先定义后使用 示例

myfun=function(x,y=2,z=3)
{
a=x*y+z
return(a)
}
myfun(x=0)

调用其他文件中的函数,先

source("myfun.R")
y=myfun(x=0)

把source当模块用

source("myfun.R")

绘图

函数

plot 散点点线
hist 直方图
boxplot
par 
...等等

线性,比例尺等也可以调控,都可以查表,或者单独做一个教程,此文不多讲,(估计没)有时间重开一个文章写绘图

pretty生成特定坐标,用于画图

pretty(x,n) 创建美观的分割点通过选取n+1个等间距的取整值将一个连续型变量x分割为n个区间
pretty(1:15)          # 0  2  4  6  8 10 12 14 16
pretty(1:15, h = 2)   # 0  5 10 15
pretty(1:15, n = 4)   # 0  5 10 15
pretty(1:15 * 2)      # 0  5 10 15 20 25 30
pretty(1:20)          # 0  5 10 15 20
pretty(1:20, n = 2)   # 0 10 20
pretty(1:20, n = 10)  # 0  2  4 ... 20

格式化输入输出

cat输出信息到屏幕/文件

#源码
#!/usr/bin/env Rscript
cat("the total number is:",3+5,"\n")
#sep指定各个间隔,file="输出到文件"
cat(1:4,"\n",sep="\t",file="catout") #默认清除源文件输入
cat(1:4,sep=c(",","!","@","\n"),file="catout", append = TRUE) # append = TRUE追加到末尾
#sep也可以指定每一个字符之间的间隔
#也可使用   sep=c(rep("\t",3),"\n") 
cat(1:4,sep=c(rep("\t",3),"\n"),file="catout", append = TRUE)
x=3+4
print(x)
#运行
cndaqiang@win10:~/code/R/Rcode> ./test.R
the total number is: 8
[1] 7
cndaqiang@win10:~/code/R/Rcode> cat catout
1       2       3       4
1,2!3@4
1       2       3       4

还可以这样

com=system("ls")
cat(com)
#结果
cndaqiang@win10:~/code/R/Rcode> ./test.R
catout    Cv.dat    myf.f90  myf.so           RCv.R       test      windowscode
cnq.R     ggplot.r  myf.mod  myfun.R          Rgraph.R    test.dat
cnq.Rout  input.R   myf.o    PErrorBar_CNQ.R  Rplots.pdf  test.R

print的输出会带很多信息,如上

文件读写

手输入edit()

new=edit()

然后填入c(1,2,3)即等价于new=c(1,2,3)
也可以在edit()括号中加入其他变量,进行修改,不影响原数据

> new=edit(c(2,3,4))

从文本文件读入

read.table("文件名"[,header=T/F,sep="间隔标识符",row.name="变量名"]) 
Cv=read.table("Cv.dat",skip=1) #skip跳过行号
  • header是否将数据文件的第一行识别为变量名,如
    > read.table("my.dat")
     V1  V2  V3
    1 num one two
    2   1  11  21
    3   2  12  22
    4   3  13  23
    > read.table("my.dat",header=T)
    num one two
    1   1  11  21
    2   2  12  22
    3   3  13  23
    
  • sep数据间隔,默认为一个或多个空格、制表符、换行符或回车符,默认就行
  • row.name=”变量名”,将变量名所在列作为行名
    > read.table("my.dat",header=T,row.name="num")
    one two
    1  11  21
    2  12  22
    3  13  23
    

保存数据

write.csv(Dipole,"RDipole.csv",row.names =F)
#======保存图片
#
pdf("RDipole.pdf")
#png("filename.png")
#jpeg("filename.jpeg")
#bmp("filename.bmp")
#
#
#
#
for ( direct in 2:4)
{
plot(Dipole[,1],Dipole[,direct],type="l",main = "Dipole", xlab = names(Dipole)[1],  ylab = names(Dipole)[direct])
}
dev.off()

#控制画板
png(file="my.png",width=800,height=30000)

write.csv默认sep是,,其他具体参数同write.table

write.table (x, file ="", sep ="", row.names =TRUE, col.names =TRUE, quote =TRUE)
  • x:需要导出的数据
  • file:导出的文件路径
  • sep:分隔符,默认为空格(” “),也就是以空格为分割列
  • row.names:是否导出行序号,默认为TRUE,也就是导出行序号
  • col.names:是否导出列名,默认为TRUE,也就是导出列名
  • quote:字符串是否使用引号表示,默认为TRUE,也就是使用引号表示

貌似write.csv在输出data.frame时还是比不过write.table

> write.table(TDEFIELD_all,file="test",row.names =F, col.names =FALSE)
> write.csv(TDEFIELD_all,file="test",row.names =F, col.names =FALSE)
Warning message:
In write.csv(TDEFIELD_all, file = "test", row.names = F, col.names = FALSE) :
  不能修改'col.names'

其他

R调用Fortran

R会自动编译文件
myf.90文件

subroutine myfun(a,a2)
    REAL(8) :: a,a2 #Fortran必须使用双精度
    open(unit=20,file="test")
    write(20,*) a,a2
    close(20)
end subroutine myfun

subroutine ofun(a,b)
    REAL(8) :: a,b
    b=a**10
end subroutine ofun

cnq.R文件

#!/usr/bin/env Rscript
system("R CMD SHLIB myf.f90")  #编译动态库
#好像会自动根据时间戳,判断是否需要重新编译
dyn.load("myf.so") #导入动态库
#.Fortran("函数名",[索引1=]as.类型(输入数据),[索引2=]as.类型(输入数据).....)
#索引任意也可以没有
#计算结果为list型,使用[ ]或$提取计算结果
input=0
output=0
c=.Fortran("ofun",inp=as.double(2),outp=as.double(4))["outp"]
.Fortran("myfun",a=as.double(2),a2=as.double(c))

查看计算结果

cndaqiang@win10:/mnt/e/work/CODE/R/Rcode> cat test
   2.0000000000000000        1024.0000000000000

R调用Python

# 安装reticulate包
install.packages("reticulate")
# 加载reticulate包
library(reticulate)
> os=import("os")
> os$system("python --version")
Python 2.7.16 :: Anaconda, Inc.
[1] 0
result=import("pyramids.io.result")
dipo=result$getDipolePython()

 

有趣

assign将字符串转为变量

> x
 [1]  1  2  3  4  5  6  7  8  9 10
> assign('ww',x)#转为变量并赋值x
> ww
 [1]  1  2  3  4  5  6  7  8  9 10
 
#get():返回与字符串同名的变量的值,并转成字符串
#assign():为字符串变量的字符串赋值
> get("xx")
 [1] "1" "1" "1" "1" "1" "1" "1" "1" "0" "0"
> 

system("ls")

调用系统命令system("系统命令")

com=system(c("cd ..;ls"))
#查看保存变量
cat(com)
#执行系统命令,并将结果保存到com
或者
> com=system(c("cd ..;ls"),intern=T)
> com
[1] "Ag55"        "Ag55.tar.gz" "graphenTDAP" "H2O"         "test-O3"

存储一些数据,函数,直接加载调用

> a='hello'
> myfun=function()
+ {
+ cat("hello R")
+ }
> save.image("cndaqiang.RData")
> q()

然后导入

# 导入数据
> load("cndaqiang.Rdata")
# 已经有了
> a
[1] "hello"
> myfun
function()
{
cat("hello R")
}
#查看变量空间有哪些内容
> ls.str()
a :  chr "hello"
myfun : function ()

将R语言,变成Linux的小脚本使用

在超算上编译安装R到个人目录

编译时候发现,R语言编译时有lapack,不错不错

wget https://mirrors.tuna.tsinghua.edu.cn/CRAN/src/base/R-3/R-3.6.0.tar.gz
tar xzvf R-3.6.0.tar.gz 
cd R-3.6.0/
module load  mpi/mvapich2/gnu/2.3b
./configure --prefix=/public/home/cndaqiang/soft/gcc-MVAPICH/R-3.6.0
make #不支持-j20
make install
PATH=/public/home/chendq/soft/gcc-MVAPICH/R-3.6.0/bin:$PATH
#--------------------------------------------
#Error :configure: error: libcurl >= 7.22.0 library and headers are required with support for https
#参考https://www.chengzi520.com/?p=1912
wget https://curl.haxx.se/download/curl-7.57.0.tar.gz
tar zxvf curl-7.57.0.tar.gz
cd curl-7.57.0
./configure --prefix=/home/cndaqiang/software/gcc-4.8.5-openmpi-1.10.7/curl-7.57.0
make && make install 
# 在~/.bashrc添加:
#R
CURLDIR=/home/cndaqiang/software/gcc-4.8.5-openmpi-1.10.7/curl-7.57.0
export PATH=$CURLDIR/bin:$PATH
export LD_LIBRARY_PATH=$CURLDIR/lib:$LD_LIBRARY_PATH
#使环境变量生效:
source ~/.bashrc
#--------------------------------------------------
#WARNING:configure: WARNING: neither inconsolata.sty nor zi4.sty found: PDF vignettes and package manuals will not be rendered optimally 警告
#我们忽视了这个错误,有人解决如下
#https://my.oschina.net/u/2429108/blog/552977
3) 下载安装 inconsolata.sty
下载 :wget http://mirrors.ctan.org/fonts/inconsolata.zip
解压 :Unzip inconsolata.zip 
将文件拷贝到目录下:
        cp -Rfp inconsolata/* /usr/share/texmf/
 刷新sty :
        mktexlsr

脚本

文件开头

#!/usr/bin/env Rscript
margs=commandArgs()
#参数
#print(margs)
#cndaqiang@win10:/mnt/e/work/CODE/R/Rcode> ./input.R 1 2 3
#[1] "/usr/lib64/R/bin/exec/R" "--slave"
#[3] "--no-restore"            "--file=./input.R"
#[5] "--args"                  "1"
#[7] "2"                       "3"
#linux中,第6个开始才是用户输入的参数
inpara=margs[6:length(margs)]
#常用包
#ggplot2
library(ggplot2)

快速画图

#!/usr/bin/env Rscript
##################################
# Author : cndaqiang             #
# Update : 2019-05-11            #
# Build  : 2019-05-11            #
# What   : 画x,y图               #
##################################
#头文件
#函数库
source("~/scripts/R-module/head.R")
#source(paste(CNQMODULEPATH,"xxx.R"))
datafile=inpara[1] #第一个输入参数为数据文件名
skipline=inpara[2] #第二个输入参数为是否跳过第一行
rplot=function(datafile,xlable="x",ylable="y",qtitle="",skipline=0)
{
#datafile="Cv.dat" #输入文件
#xlable="x" #x标签
#ylable="y" #y标签
#skipline=1 #跳过开头的几行
#读入数据
    mydata=read.table(datafile,skip=skipline)
    datacol=length(mydata[1,])
    datarow=length(mydata[,1])
    if (datacol == 1)
    {
        plot(mydata,type="b",main = qtitle, sub = "", xlab = xlable,  ylab = ylable)
    }else
    {

        ydata=mydata[,2:datacol] #
        xdata=mydata[,1] #第一列
        for ( i in 2:datacol )
        {
            y=mydata[,i]
            x=mydata[,1]
            plot(x,y,type="b",main = qtitle, sub = "", xlab = xlable,  ylab = ylable)
        }
    }
}

rplot(datafile,xlable=inpara[3],ylable=inpara[4],qtitle=inpara[5],skipline=skipline)

集成到makefile里,编译计算,画图,简直不要太爽

test:$(TEST)
	$(FC)  $(FFLAGS) -o test $(OBJS) $(LIBS)
	./test > result.dat
	./plot.r

调用shell的grep组合操作

提取矩阵

Dipole=paste("grep 'Electric dipole (a.u.)' ",result,"|grep 'TDAP' |awk '{printf \"%f,%f,%f,\",  $(NF-2),$(NF-1),$NF }' ")
Dipole=system(Dipole,intern=T)
Dipole=as.numeric(as.vector(unlist(strsplit(Dipole,",")))) #分离后,默认是list,转为向量和数值型
Dipole=c(Dipole)
Dipole=matrix(Dipole,ncol=3,byrow=T)

读入fdf格式

#======读入时间
MD.LengthTimeStep=paste("grep -i '^MD.LengthTimeStep'",result,"|awk '{ print $2 }'")
#此处默认单位是fs了,不对fs进行判断
MD.LengthTimeStep=as.numeric(system(MD.LengthTimeStep,intern=T))
MD.FinalTimeStep=paste("grep -i '^MD.FinalTimeStep'",result,"|awk '{ print $2 }'")
MD.FinalTimeStep=as.numeric(system(MD.FinalTimeStep,intern=T))

读入fdf block

#行号
phyHang=paste("grep -ni '^%block myblock'",input.fdf,"|awk -F: '{ print $1 }'")
phyHang=paste("hang=$(",phyHang,")")
phyHang=paste(phyHang,"hang=$(echo ${hang}+1|bc)",sep=";")
#读入
phy=paste(phyHang,";sed -n ${hang}p",input.fdf)
phy=system(phy,intern=T)
#
phy=as.vector(unlist(strsplit(phy," +"))) #以空格分开,strsplit生成的是list,要转为向量
phy=c(as.numeric(phy))

phydata=phy[1]       #有物理含义
phydata=phy[2]   #有物理含义
phydata=phy[3]       #有物理含义
phydata=phy[4]           #有物理含义
phydata=phy[5]           #有物理含义

plot

使用plot进行画图时,数据框的一列,一定要写成mydata[,lie],把一列提成数组,那个,不能少 示例,

pdf(file="RCurrentFFT.pdf")
for ( direct in 2:colnum)
{
plot(myfft[1:myn/2,1],myfft[1:myn/2,direct],type="l",main = "Current FFT", xlab = paste(names(myfft)[1],"FFT"), ylab = paste(names(myfft)[direct],"FFT") )
#----下面为加text到图中
minfft=max(myfft[1:myn/2,direct])/10
peak=(myfft[2:(myn/2),direct] > myfft[1:(myn/2-1),direct]) & (myfft[2:(myn/2),direct] > myfft[3:(myn/2+1),direct]) & ( myfft[2:(myn/2),direct] > minfft ) 
peak=c(T,peak,(1:(myn-2) < 0))
text(myfft[peak,1],myfft[peak,direct],format(myfft[peak,1],digits = 3),cex=1,pos=1,col="red")
#x,y,显示的内容
}
dev.off()

ggplot2

代码的发展,很多函数都变了,注意,有的教程给的函数不能用

qplot

qplot(数据框索引,数据框索引,data=数据框)
> price=1:10
> num=sin(price)
> mydata=data.frame(price,num)
qplot(price,num,data=mydata)
#也可以使用数学函数处理
qplot(price,log(x),data=mydata)
#也可以是数据框内多组数据的计算结果
qplot(price,log(x)*num,data=mydata)

根据分类,确定颜色和尺寸

> qplot(price,num,data=data,colour=good,shape=xx,size=num,alpha=num/10)
#shape=字符型,不然报错
#color字符数子都行,图形颜色
#size指定大小
#alpha透明度[0-1]之间
> mydata
   price        num xx good
1      1  0.8414710  1    0
2      2  0.9092974  1    0
3      3  0.1411200  1    0
4      4 -0.7568025  1    0
5      5 -0.9589243  1    0
6      6 -0.2794155  1    1
7      7  0.6569866  1    1
8      8  0.9893582  1    1
9      9  0.4121185  0    1
10    10 -0.5440211  0    1

geom指定形状

qplot(price,num,data=data,geom="line")
#geom="path","line"连线
#geom="point" 默认散点
#"smooth"平滑,使用span,method等继续控制,
#qplot(num,price,data=data,geom=c("smooth"),span=0.5,method="gam")
#也可以组合多种图形
qplot(price,num,data=data,geom=c("point","line"))

保存

ggsave("my.png")

貌似只保存最近的图

ggplot

#画图
png(outfile)
mypicture=ggplot(mydata,aes(x=xdata,y=ydata,color=orderdata))+geom_point()
#加水平垂直参考线(阈值线),由输入参数Qylimit和Qxlimit控制
mypicture=mypicture + geom_hline(yintercept=Qylimit,linetype=4)
mypicture=mypicture + geom_vline(xintercept=c(-1*Qxlimit,Qxlimit),linetype=4)
#设置x,y上下限
mypicture=mypicture + xlim(Lxlimit,Hxlimit)
mypicture=mypicture + ylim(Lylimit,Hylimit)
#设置点的颜色
mypicture=mypicture + scale_color_manual(values =c("Up" = "red", "Down" = "blue", "Normal" = "grey"))
#x,y轴label
#expression是为了输出上下标
mypicture=mypicture + labs(x=expression(log[2](FC)),y=expression(-log[10](FDR)),color="significant")

mypicture   
dev.off()

3d画图

zypper install Mesa-libGL-devel

R语言用于科学计算

FFT

FFT

数学课本f(t)-c(w)
代码t=t0+n*dt, f&n->FFT->C(n)
自己将n转为f(n),T(n)

#!/usr/bin/env Rscript
##################################
# Author : cndaqiang             #
# Update : 2019-05-20            #
# Build  : 2019-05-20            #
# What   : FFT分析               #
##################################

fn1=5;fn2=10;fn3=100;fn4=233 #频率
dt=0.002 #采样间隔
t=seq(0,1,dt) #样本点
n=length(t)  #采样数
y=sin(2*pi*fn1*t)+3*cos(2*pi*fn2*t)+4*sin(2*pi*fn3*t)+5*sin(2*pi*fn4*t) #采样值
#plot(x,y)
z=abs(fft(y)/n)  #频谱 
f=(1:n-1)/n/dt  #频率
png("fft.png")
plot(f[1:n/2],z[1:n/2],type="l")
dev.off()

结果图
图中的峰值比例不对,加密N点减少dt后就准了,就不补图了
由于输入是实数,这里f=(1:n-1)/n/dt #频率应该变为f=(1:n-1)/(n/2)/dt #频率

IFFT示例

FFT的结果需要除以N才能得到正确的强度

> a=1000
> x=1:a;x=x*2*pi/a;y=exp(1i*x);z=fft(y);  z1=fft(fft(y),inverse = TRUE);  z2=fft(fft(y)/a,inverse = TRUE)
> max(abs(z))
[1] 1000
> max(abs(z1))
[1] 1000
> max(abs(zz))
[1] 10000
> max(abs(z2))
[1] 1
> 

本文首发于我的博客@cndaqiang.
本博客所有文章除特别声明外,均采用 CC BY-SA 4.0 协议 ,转载请注明出处!



评论


广告

目录

广告
访客数据