CとRを動かす

Cの基本 R
  • 久しぶりにCの復習をしてみた
  • シェルでCとRを動かせるようにシェルスクリプトを書いた
    • Cを計算に使いたい
    • Rで図を描きたい
  • いままでマスで仕切ってマスごとに計算するというシミュレーションをやってきていたので、Cでも同じことをしてみる
  • マスの大きさをシェルスクリプトから与えられるようにしてみた
  • 下の3つを同じディレクトリに置いておいて、シェルスクリプトを実行すればよい
  • 全体の流れとしては
    • シェルで3次元空間の大きさを与える（Nx、Ny、Nz）
    • CがNx、Ny、Nzの大きさの空間内の各マスの中身を計算する（3次元の配列）
    • Cが3次元の配列の中身を書き出す（Numbers.txt）
    • RがNumbers.txtを読み込む
    • Rのrglパッケージで3dplotして図を描きだす（Array.png）
    • シェルでArray.pngを開き、Numbers.txtは消去する
#!/bin/sh

Nx=30
Ny=30
Nz=30

# Cをコンパイル、実行まで
gcc -o mat mat.c
./mat $Nx $Ny $Nz

# Rを実行
R --vanilla --slave --args < mat.R $Nx $Ny $Nz

# ファイルの実行と消去
open Array.png
rm Numbers.txt

/* mat.c */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define PI 3.1415926535

int main (int argc, char *argv[]){
	
	int Nx, Ny, Nz;
	int i,j,k;
	
	FILE *fpA;
	
	Nx = atoi(argv[1]);
	Ny = atoi(argv[2]);
	Nz = atoi(argv[3]);
		
	double A[Nx][Ny][Nz];
	
	for (k = 0; k < Nz; k++) {
	for (j = 0; j < Ny; j++) {
	for (i = 0; i < Nx; i++) {
		
		A[i][j][k] = (i + 1)*(j + 1)*(k + 1)*1.0;
		
		fpA = fopen("Numbers.txt","a");
		if(fpA == NULL)return 1;
		fprintf(fpA,"%f\n",A[i][j][k]);
		fclose(fpA);
	}}}
	
	return 0;
}

# mat.R
library(rgl)

Nx <- as.numeric(commandArgs(trailingOnly = TRUE)[1])
Ny <- as.numeric(commandArgs(trailingOnly = TRUE)[2])
Nz <- as.numeric(commandArgs(trailingOnly = TRUE)[3])

m <- data.matrix(read.table("Numbers.txt",header=FALSE,sep=""))
A<-array(m,c(Nx,Ny,Nz))

plot3d(slice.index(A,1),slice.index(A,2),slice.index(A,3),col=rainbow(256)[A/(max(A))*256+1],alpha=ifelse(A!=0,1,0))
rgl.snapshot("Array.png")