ไปเจอมา ขำดี แต่มันก็จริงนะครับ 555
Computational scientists solve tomorrow’s problems with yesterday’s computers; computer scientists seem to do it the other way around.
– anonymous
สักงอย ดอท คอม
ไปเจอมา ขำดี แต่มันก็จริงนะครับ 555
Computational scientists solve tomorrow’s problems with yesterday’s computers; computer scientists seem to do it the other way around.
– anonymous
เป็นตัวอย่างที่ผมเขียนขึ้นตอนใช้สอนการเขียน 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 ใน 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];
เพิ่งอัพเดท 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