Author: Bright

computer scientists vs computational scientists

ไปเจอมา ขำดี แต่มันก็จริงนะครับ 555

Computational scientists solve tomorrow’s problems with yesterday’s computers; computer scientists seem to do it the other way around.

– anonymous

a simple pendulum in R

เป็นตัวอย่างที่ผมเขียนขึ้นตอนใช้สอนการเขียน function ในภาษา R ครับ

library(magick)
# the formula for the changing of angle over time 
# https://en.wikipedia.org/wiki/Pendulum
theta <- function(t,period,theta0){
  theta0*cos(2*pi*t/period)
}

# for drawing the pendulum at a given point (x,y), i.e. drawing a line between (0,0) and (x,y)
draw.pendulum <- function(pen.x, pen.y){
  
  pivot.x <- 0
  pivot.y <- 0

  plot(x=pivot.x,y=pivot.y,type = 'l',ylim = c(-1.5,0.1),asp = 1)
  lines(x=c(pivot.x,pen.x),y=c(pivot.y,pen.y))
  points(x = pen.x, y=pen.y,col='red')
}

# calculate the position of the pendulum over time
# l-length of the pendulum
# theta0 - the starting angle (in radian)
# runtime - the simulation time (in seconds)
# step - time step
pendulum <- function(theta0,l,runtime,step = 0.1){
  g <- 9.8
  period <- 2*pi*sqrt(l/g)
  
  times <- seq(0.0,runtime, by = step)

  pen.x <- l*sin(theta(times, period = period,theta0 = theta0)) 
  pen.y <- -l*cos(theta(times,period = period,theta0 = theta0))
  
  output <- matrix(c(times,pen.x,pen.y),ncol = 3)
  return(output)
}


pendulum.sim <- function(theta0, l, runtime, step = 0.1){
  
  sim.data <- pendulum(theta0 = theta0, l=l, runtime = runtime)
  
  ## export each frame as png
  for(i in 1:nrow(sim.data)){
    png(filename = paste0(tempdir(),"/",i,".png"))
    draw.pendulum(sim.data[i,2],sim.data[i,3])
    dev.off()
  }
  
  frames <- image_read(paste0(tempdir(),"/",1:nrow(sim.data),".png"))
  animate <- image_animate(image = frames, fps = 1/step)
  
  return(animate)
  
}

animate <- pendulum.sim(theta0 = pi/4, l = 1, runtime = 10)
print(animate)

OpenCL ในภาษา Wolfram

มีคนถามมาว่าช่วยทำตัวอย่างให้ดูหน่อยในการใช้ OpenCL ใน Mathematica ซึ่งผมก็บอกไปแล้วว่ามันตัวอย่างเยอะเลยในhelp Mathematica เอง 🙂 แต่ก็ยังอยากให้ผมทำให้ดู ตัวอย่างที่ทำให้ดูก็เป็นการประมาณค่า Pi โดยอาศัยที่เรียกว่า Monte carlo โดยประมาณจากการสุ่มจุดลงในสี่เหลี่ยมจัตุรัสที่มีวงกลมรัศมีครึ่งหนึ่งของความยาวด้านสี่เหลี่ยมจัตุรัส แล้วดูว่ามีกี่จุดที่อยู่ในวงกลม ซึ่งค่า Pi จะสามรถประมาณได้จาก

Pi = 4 * จำนวนจุดในวงกลม / จำนวนจุดที่สุ่มมาทั้งหมด

อันนี้เป็น code ที่ผมเขียนโดยใช้ OpenCL ใน Mathematica เพื่อหาดูว่าจุดไหนบ้างที่อยู่ในวงกลมซึ่งมีจุดศูนย์กลาง (0.5,0.5) ในสี่เหลี่ยนจัตุรัส

<< OpenCLLink`

scr = "
  __kernel void myPie( __global double *X, __global double *Y, 
__global mint *inC, mint np) {
      int indx = get_global_id(0);
      double d; //distance
  	if (indx < np) {
  		d = sqrt((X[indx]-0.5)*(X[indx]-0.5)+(Y[indx]-0.5)*(Y[indx]-0.5));
          if(d <= 0.5)
inC[indx] = 1;   
  	}
  }
  ";

fn = OpenCLFunctionLoad[scr, 
  "myPie", {{_Real, _, "Input"}, {_Real, _, "Input"}, {_Integer, _, 
    "Output"}, _Integer}, 16 ];

np = 10^6;
x = OpenCLMemoryLoad[RandomReal[1, np]];
y = OpenCLMemoryLoad[RandomReal[1, np]];
inc = OpenCLMemoryLoad[ConstantArray[0, np]];
Print["my Pi : ", (OpenCLMemoryGet@(fn[x, y, inc, np] // First) // 
    Total)*4.0/np]

OpenCLMemoryUnload[x, y, inc]

ทดลองรันดูสัก 100 ครั้ง โดยที่แต่ล่ะครั้งสุ่มมา 10^7 จุด แล้วดูว่ามันให้ค่าประมาณแตกต่างเท่าไหร่เมื่อเทียบกับค่าจริง (จาก Mathematica) จากรูปด้านล่าง เส้นประสีดำคือที่ได้จากการประมาณ สีแดงคือเส้นค่า Pi จาก Mathematica

ส่วนอีกอันก็เป็นการหมุนภาพ โดยที่มีจุดหมุนอยู่ที่กลางภาพครับ 

scr = "
  __kernel void test1( __global float * mat1, __global float * mat2, 
float theta, mint width, mint height) {
      int xIndex = get_global_id(0);
  	int yIndex = get_global_id(1);
  	int indx = xIndex + yIndex*width;
  
  	 if (xIndex >= width || yIndex >= height)
  		return ;
  
  	//image centre
  	float x0 = width/2;
  	float y0 = height/2;
  
  	//relative position
  	int xprime = xIndex-x0;
  	int yprime = yIndex-y0;
  
  	float sinTheta = sin(theta);
  	float cosTheta = cos(theta);
  
  	int nx =(xprime*cosTheta - yprime*sinTheta + x0);
  	int ny =(xprime*sinTheta + yprime*cosTheta + y0);
  	mat2[indx] = mat1[nx + width*ny];
  }
  
  ";

Clear[test];
test = OpenCLFunctionLoad[scr, 
  "test1", {{"Float", _, "Input"}, {"Float", _, "Output"}, 
   "Float", _Integer, _Integer}, {16, 16}, 
  "ShellOutputFunction" -> Print];

ว่าด้วยเรื่อง Paclets ของภาษา Wolfram

เพิ่งอัพเดท Mathematica ไปใช้ 12.1 แล้วเจอเรื่อง paclet ที่เปลี่ยนไป (ในทางที่ดีขึ้น) เลยเดี๋ยวจะมาเล่าให้ฟัง ขอแปะแหล่งข้อมูลไว้ก่อน

1. https://www.wolframcloud.com/obj/tgayley/Published/PacletDevelopment.nb

2. https://www.wolfram.com/broadcast/video.php?sx=paclet&v=2833

3. https://reference.wolfram.com/language/guide/Paclets.html