ไปเจอมา ขำดี แต่มันก็จริงนะครับ 555
Computational scientists solve tomorrow’s problems with yesterday’s computers; computer scientists seem to do it the other way around.
– anonymous
Categories
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)
Categories
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];
เพิ่งอัพเดท 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
Categories
เว้นวรรค
มีคนส่งมาถามว่ามีข้อความเป็นภาษาอังกฤษแต่อยากให้เว้นวรรคหรือใส่space ให้หน่อยถ้ามันขึ้นต้นด้วยตัวพิมพ์ใหญ่ เช่น “AasdfasfdAsdfasd” อยากให้มันเป็น “Aasdfasfd Asdfasd” ผมเขียน script ง่ายให้ดังนี้ครับ
AddSpace[str_] := (StringReplace[str, x_?UpperCaseQ :> ” ” <> x ] // StringTrim)
ซึ่งเป็นภาษา Wolfram ครับ
อันนี้เอามาลองกับชื่อจังหวัดที่พิมพ์ติดกันมา
ที่อยากจะเน้นคือเรื่อง pattern ในภาษา Wolfram นี่มีประโยชน์มาก ทำให้เราทำงานได้เร็วขึ้นครับ
กำลังเตรียมdraft งานเก่าที่ทำมาหลายปีแล้วไม่ได้มีโอกาสส่งตีพิมพ์สักที เลยเอามาปัดฝุ่นใหม่โดยใช้ Mathematica รู้สึกได้เลยว่า code มันเขียนไม่กี่บรรทัดเอง แถมแชร์ก็ง่าย https://www.wolframcloud.com/obj/sompob/Published/NJW_Im4a.nb
Categories
มันจะกระจายยังไงถ้าไม่ทำอะไรเลย
ทำเล่นๆขำ ว่าถ้าไวรัสที่ว่ามันกระจายไปตามจังหวัดต่างๆที่ว่า จำนวนคนติดเชื้อมันพุ่งเยอะขนาดไหน แบบว่าถ้าคนมันคิดว่าเป็นหวัดธรรมดา
โห … ตัวเลขพุ่งขึ้นเร็วมากกก ถ้าเห็นโมเดลทำนายเรื่องอะไรก็ตามอยากให้คิดกันเสมอว่า
“All models are wrong, but some are useful”
ฉะนั้นอย่าเพิ่งไปเชื่ออะไรง่ายๆครับ ทุกโมเดลมันมีสมมุติฐานในการสร้างของมัน ระบบในธรรมชาติมันซับซ้อนมากครับ ไม่มีใครที่จะโมเดลมันได้หมดหรอกครับ ทุกโมเดลมันจะตัดความยุ่งยาก ซับซ้อนออกเพื่อให้โมเดลได้ง่ายครับ
Categories
มาลองใช้ Stan แก้ปัญหา
พอดีว่าง ระหว่างรอผลคำนวณอะไรบางอย่าง ผมเห็นมีคนโพสท์ถามที่เวบพันทิปตามนี้ครับ
จากคำถามนี้เราสามารถที่จะใช้ Stan แก้ปัญหาได้ถ้าผมมองว่าการวัดเปรียบเทียบเครื่องมือมาตรฐาน(จริง)กับที่ปรับปรุงขึ้นมานั้นในแต่ล่ะครั้งนั้นไม่ได้เกี่ยวกันเลย และค่าจากวัดของเครื่องที่ปรับปรุงนั้นมีการกระจายรอบค่าจากเครื่องมาตราฐานแบบ normal distribution โดยมีค่า standard deviation หรือ SD อยู่ค่าหนึ่งที่เป็นตัวกำหนดประสิทธิภาพของเครื่อง อย่างเช่น จากข้อมูลที่ให้มาเมื่อวัดเทียบกับเครื่องมือจริงที่วัดได้ 2 แต่ค่าจากเครื่องปรับปรุงวัดมา 3 ครั้งได้ (1.8,2.1,1.9) เราจะสมมุติให้ทั้ง 3 ค่านี้กระจายรอบค่าใดค่าหนึ่ง โดยที่ ถ้าเราปรับปรุงได้เจ๋ง ค่าที่มันกระจายรอบนี้มันก็ควรจะได้เท่ากับค่าจากเครื่องจริงโดยที่มีค่า SD น้อยที่สุดเท่าที่จะเป็นไปได้ แต่ในความเป็นจริงก็จะมีผลผันผวนหรือerrorต่างๆเข้ามาเกี่ยวข้อง ที่ผมจะทำให้ดูนี้เราจะลองใช้โปรแกรมอย่าง Stan มาช่วยหาดูว่าค่ากลางที่เครื่องปรับปรุงหรือสร้างขึ้นมานั้นมันกระจายรอบในแต่ล่ะครั้งของการวัดมันคือค่าอะไรและมี SD เท่าใด
ผมใช้โปรแกรมที่ชื่อ Stan ในภาษา R ผ่าน library ที่เรียกว่า rstan ครับ ด้านล่างนี้ก็เป็น code ที่ผมใช้กับปัญหานี้ครับ
#โหลด library ที่จะใช้งานครับ
#ในที่นี้มี 2 ตัวคือ rstan กับ bayesplot (ใช้วาดกราฟสรุปผล)
library(rstan)
library(bayesplot)
#เตรียมข้อมูล โดยผมแยกตามค่าของเครื่องจริง
data <- list(
ob2 = c(1.8,2.1,1.9),
ob3 = c(3.2,3.0,3.2),
ob4 = c(3.9,4.3,4.4),
ob5 = c(5.0, 5.2, 5.5),
ob6 = c(6.1,5.9,6.5),
n=3
)
# code ของโมเดล
model <- "
data{
//จำแหนกประเภทของตัวแปรของข้อมูล ว่าเป็นเลขจำนวนเต็มหรือทศนิยม พร้อมกำหนดขนาด
int<lower=1> n;
real ob2[n];
real ob3[n];
real ob4[n];
real ob5[n];
real ob6[n];
}
parameters{
// กำหนดประเภทของตัวแปรที่จะใช้ในโมเดล ซึ่งในที่นี้คือค่ากลางที่ข้อมูลมันกระจายรอบ ซึ่งแบ่งตามค่าจริง
// sig เป็นค่า SD
real ob2mu;
real ob3mu;
real ob4mu;
real ob5mu;
real ob6mu;
real<lower=0> sig;
}
model{
//กำหนดว่าค่ากลางมันและSD อยู่ในช่วงไหน จากการกระจายแบบไหน
//ซึ่งในที่นี้ผมให้มันมาจาก uniform distribution โดยใส่ช่วงที่คิดว่าค่ามันจะอยู่ในนั้น
sig ~ uniform(0,3);
ob2mu ~ uniform(0,7);
ob3mu ~ uniform(0,7);
ob4mu ~ uniform(0,7);
ob5mu ~ uniform(0,7);
ob6mu ~ uniform(0,7);
// กำหนดให้ค่าที่วัดมาในแต่ล่ะค่าจริงกระจายแบบ normal รอบค่ากลางอันหนึ่ง
ob2 ~ normal(ob2mu, sig);
ob3 ~ normal(ob3mu, sig);
ob4 ~ normal(ob4mu, sig);
ob5 ~ normal(ob5mu, sig);
ob6 ~ normal(ob6mu, sig);
}
"
fit <- stan(model_code=model, data = data, iter = 10000)ผลลัพท์ที่ได้ก็จะประมาณนี้ครับ
เห็นได้ว่าค่ากลางหรือ median ที่ 50% นั้นค่อนข้างจะใกล้กับค่าจริงคือ เช่นที่ค่าจริง = 2 เครื่องปรับปรุงทำได้ 1.93(95%CI: (1.63,2.24)) หรือที่ 6 เครื่องปรับปรุงทำได้ที่ 6.15 (CI: (5.86,6.48)) และ SD = 0.25 CI: (0.17-0.44)
posterior <- as.matrix(fit)
plot_title <- ggtitle("Posterior distributions",
"with medians and 95% intervals")
mcmc_areas(posterior,
pars = c("ob2mu","ob3mu","ob4mu","ob5mu","ob6mu", "sig" ),
prob = 0.95) + plot_title
posterior2 <- extract(fit, inc_warmup = T, permuted = FALSE)
color_scheme_set("mix-blue-pink")
p <- mcmc_trace(posterior2, pars = c("ob2mu","ob3mu","ob4mu","ob5mu","ob6mu", "sig" ),
n_warmup = 5000,
facet_args = list(nrow = 6, labeller = label_parsed))
p + facet_text(size = 15)จากกราฟการกระจายที่ได้จะเห็นได้ว่าผมใส่เป็น uniform distribution ที่ช่วงระหว่าง 0-7 เลยแต่ผลลัพท์ที่ได้เป็น normal distribution ครับ
โดยปกติ Mathematica หรือภาษา Wolfram จะมีปัญหากับการเรียก function ที่สร้างจาก Compile ใน function ที่จะสร้างจาก Compile อีกทีได้ (ส่วนใหญ่เป็นเรื่องที่เกี่ยวกับว่า Funcions ไหนที่มัน compile ได้หรือไม่ได้) แต่ใน version 12 นี้เราสามารถที่จะทำแบบนั้นได้ ถ้า function เหล่านั้นถูก compile ด้วยคำสั่ง FunctionCompile ซึ่งผมมองว่ามันสะดวกอย่างมาก และที่สำคัญเราสามารถที่จะ export สร้างเป็น library ไฟล์ได้ด้วย โดยคำสั่ง FunctionCompileExport สกุลไฟล์ที่สร้างได้ก็ได้แก่
ตัวอย่างการเรียกใช้ compiled function ใน compiled function อีกที ก็สมมุติให้ผมมี function สำหรับทำ bisection เพื่อหาค่าที่ทำให้ function ที่สนใจเป็น 0 หรือใกล้ 0 มากๆ
BisectionMethod = FunctionCompile[
Function[
{Typed[f, {"Real64"} -> "Real64"], Typed[lim0, "Real64"], Typed[lim1, "Real64"]},
Module[{Maxiter = 100, tol = 10.^-8, iter = 0, x0 = lim0, x1 = lim1, mid = 0.0, f0, f1, fmid},
f0 = f[x0];
f1 = f[x1];
While[Abs[x1 - x0] >= tol && iter++ < Maxiter,
mid = (x1 + x0)/2.;
fmid = f[mid];
If[Xor[fmid > 0, f0 > 0],
x1 = mid; f1 = fmid,
x0 = mid; f0 = fmid]];
mid]]]โดยfunction ที่ผมจะหา root คือ f(x) = sin(x)+exp(x)
f = FunctionCompile[Function[Typed[arg, "Real64"], Sin[arg] + Exp[arg]]]จะเห็นได้ว่าทั้ง BisectionMethod และ f ถูกสร้างจาก FunctionCompile แต่ f จะถูกเรียกเข้าไปใช้ใน BisectionMethod อีกที
BisectionMethod[f, -1.0, 1.0]-FindRoot[Sin[x] + Exp[x], {x, -1}][[1, 2]]
(* -5.18667*10^-11 *)