ตัวอย่าง code สำหรับsampling แบบขนานใน cmdstan

Mac OS/Linux

for i in {1..4} 
do 
    ./my_model sample random seed=12345 id=$i data file=my_data output file=samples$i.csv & 
done

Windows

for /l %x in (1, 1, 4) do start /b model sample random seed=12345 id=%x data file=my_data output file=samples%x.csv

เอามาจาก cmdstan-guide-2.20.0.pdf

Makefile สำหรับใช้กับ แพ็คเกจที่ใช้ RcppArmadillo

พอดีว่าเขียนแพ็คเกจที่มีเรียกใช้งาน armadillo ใน Rcpp แล้วมีปัญหาเวลา compile ในเรื่องของ linking

คำแนะนำที่เจอบ่อยๆก็คือให้สร้างไฟล์ Makefile หรือ Makefile.win ไว้ที่เดียวกันกับ src โดยที่ใน Makefile ก็มีการประกาศตัวแปรตามข้างล่างนี้เลยครับ

CXX_STD = CXX11

PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) 

PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)

หรือแต่ถ้าไม่อยากให้มีปัญหานี้ตั้งแต่แรกก็ให้ใช้ Rcpp แล้วก็สร้างแพ็คเกจจากตัวอย่างที่เขาทำไว้แล้วจากคำสั่ง RcppArmadillo.package.skeleton ครับ

Kung package แพ็คเกจกุ้ง

แพ็คเกจกุ้งเป็นแพ็คเกจภาษา R ที่ผมเขียนขึ้นมาสำหรับใช้สร้าง Shiny App สำหรับโค้ดของโมเดลที่สร้างจากสมการอนุพันธ์หรือ ODE ที่ใช้ตัว solver จากแพ็คเกจที่ชื่อ deSolve ครับ

ไอเดียก็มีเพียงว่าจากสมการ ode ที่สร้างขึ้นมาด้วยถ้าตัวแปรในสมการมีการเปลี่ยนแปลงค่า ผลลัพท์ที่ออกมามันจะเปลี่ยนแปลงไปอย่างไร การเขียนใน R นั้เราก็อาจจะใช้คำสั่ง manipulate มาช่วยได้แต่มันก็ยังมีข้อจำกัดเรื่องจำนวนของตัวแปรที่สามารถเรียกมาใช้ได้ เพราะมันขึ้นอยู่กับขนาดหน้าจอของตัว plot ผมเลยคิดว่าถ้าใช้ shiny มันน่าจะดูดีกว่าและสามารถใส่ตัวแปรได้มากกว่า ขึ้นกับการออกแบบ ui แต่การเขียน shiny ก็ไม่ใช่ว่าจะง่ายสำหรับคนเริ่มเรียนรู้ R ผมเลยคิดว่าน่าจะเขียนแพ็คเกจที่สามารถช่วยเรื่องนี้ได้

แพ็คเกจที่เขียนก็พยายามทำให้มันใช้ง่ายมากที่สุดจาก code ที่มีอยู่แล้ว โดยผู้ใช้เพิ่มคำสั่งไม่มากก็สามารถที่จะสร้าง shiny application ได้แล้ว ตัวอย่างเช่น ถ้าผู้ใช้มี code อยู่แล้วตามนี้ ซึ่งมันจะคำนวณและ plot กราฟให้ตามค่าของตัวแปรพารามิเตอร์ที่กำหนดไว้ เช่น gamma = 0.14286 และ beta = 0.6

library(deSolve)

init <- c(S = 1-1e-6, I = 1e-6, R= 0.0)
times <- seq(0, 70, by = 1)

sir <-function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    dS <- -beta*S*I
    dI <- beta*S*I - gamma*I
    dR <- gamma*I
    
    return(list(c(dS, dI, dR)))
  })
}

parameters <- c(
  gamma = 0.14286, 
  beta = 0.6
)

out <- as.data.frame(ode(y = init, times = times, func = sir, parms = parameters))

matplot(times, out[,c("S","I","R")], type = "l", xlab = "Time", ylab = "Susceptibles and Recovereds", main = "SIR Model", lwd = 1, lty = 1, bty = "l", col = 2:4)
legend(40, 0.7, c("Susceptibles", "Infecteds", "Recovereds"), pch = 1, col = 2:4)

ถ้าผู้ใช้ต้องการดูว่ากราฟจะเปลี่ยนแปลงอย่างไรถ้าตัวแปร gamma หรือ beta มีการเปลี่ยน ก็เพียงเพิ่ม code ลงไปในส่วนต่างๆ ตามนี้

!Start

library(deSolve)

init <- c(S = 1-1e-6, I = 1e-6, R= 0.0)
times <- seq(0, 70, by = 1)

sir <-function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    dS <- -beta*S*I
    dI <- beta*S*I - gamma*I
    dR <- gamma*I

    return(list(c(dS, dI, dR)))
  })
}

!Parameters
parameters <- c(
  gamma = 0.14286,
  beta = 0.6
)

!ODECMD
out <- as.data.frame(ode(y = init, times = times, func = sir, parms = parameters))

!PostProcess

!Plots
matplot(times, out[,c("S","I","R")], type = "l", xlab = "Time", ylab = "Susceptibles and Recovereds", main = "SIR Model", lwd = 1, lty = 1, bty = "l", col = 2:4)
legend(40, 0.7, c("Susceptibles", "Infecteds", "Recovereds"), pch = 1, col = 2:4)

!Controls
sliderInput("beta","beta", min = 0,max = 10,step = 0.001,value = 1.4),
sliderInput("gamma","gamma", min = 0,max = 1,step = 0.001,value = 0.14)

!End

จะเห็นว่ามี keywords ที่เพิ่มเข้าไปดังนี้

!Start – เป็นตัวบอกเริ่มต้น เราสามารถบอกรายละเอียดของตัวแปรอย่าง เวลา และค่าเริ่มต้นของ state หรือ compartment ต่างๆและ function ของระบบสมการหลักที่จะแก้ด้วย deSolve::ode ใน keywords นี้ได้เลยครับ

!Parameters – เป็นตัวแปร parameters ของโมเดล ode ครับ ซึ่งต้องกำหนดเป็น vector ตามตัวอย่างนี้เลยครับ

!ODECMD – สำหรับคำสั่งที่ต้องใช้คำสั่ง ode จาก deSolve ครับ สำหรับตอนนี้ต้องให้ตัวแปร output ชื่อ out นะครับ

!PostProcess – ถ้ามีการคำนวณค่าอะไรบางอย่างจาก out ก็ให้ใส่ใน keyword นี้ได้เลยครับ

!Controls – สำหรับใส่ sliders ของตัวแปรที่สนใจครับ

!End – สำหรับบอกว่าจบแล้ว 🙂

โดยหลังจากที่ใส่ keywords และsave เป็นไฟล์ใหม่แล้ว เช่นชื่อ mysystem.R เราสามารถสร้าง shiny app ได้ด้วยการรันคำสั่งนี้ครับ

runSystem(‘mysystem.R’)

ผลที่ได้ก็จะประมาณนี้ครับ

ถ้าสนใจอยากใช้งานก็ทำการติดตั้งได้โดยพิมพ์

devtools::install_github('slphyx/Kung')

ดูเพิ่มเติมที่ https://github.com/slphyx/Kung

เรื่อง plane geometry ใน Mathematica 12

Wolfram Mathematica 12 ที่เพิ่งจะเปิดตัวไปมีความสามารถเพิ่มขึ้นมาเยอะพอสมควรเลย และหนึ่งจากหลายความสามารถที่เพิ่มเข้ามานี้ที่ผมชอบมากก็คือความสามารถในเรื่องเรขาคณิตที่เรียกว่า Synthetic Geometry

เจ้า Synthetic Geometry ก็เป็นสาขาที่ศึกษาเรขาคณิตที่หาข้อสรุปจากการใช้เรื่องสัจพจน์หรือข้อความคาดการณ์ต่างๆเป็นเครื่องมือหลักในการแก้ปัญหา ลองดูความหมายเพิ่มเติมที่ https://en.wikipedia.org/wiki/Synthetic_geometry  นะครับ

ใน Mathematica 12 นี้เราสามารถที่จะสร้างไดอะแกรมของปัญหาได้จากคำสั่ง GeometricScene เช่น

ถ้ามีสามเหลี่ยมที่มีมุมทั้งสามมุมอยู่ที่ จุด (0,0),(1,0),(0,1) หรือสามเหลี่ยมมุมฉากที่มีฐานและความสูงคือ 1 หน่วย เราก็สามารถสร้างได้แบบนี้ครับ

จากคำสั่งด้านบน GeometricScene เราเพียงแต่บอกว่าเรามีจุดอะไรบ้าง ซึ่งในที่นี้ก็คือจุด a,b,c ที่อยู่ที่ {0,0},{1,0},{0,1} ตามลำดับ และก็บอกว่าให้สร้าง สามเหลี่ยมโดยใช้จุดที่กำหนดไว้

จากนั้นผมทดลองหาพื้นที่ของสามเหลี่ยมนี้ ซึ่งก็ผลลัพธ์ที่ได้ก็คือ 1/2

มาลองทำอะไรที่มันดูซับซ้อนขึ้นมาอีกหน่อยครับ  สมมุติว่าผมต้องการหาพื้นที่ของสามเหลี่ยมที่มีเหลี่ยมทั้งสามที่จุด k,l,n

และมีวงกลมที่ผ่านจุดทั้งสามนี้

โดยผมกำหนดให้ว่าระยะระหว่างจุด k และ n  คือ 3 หน่วย และมีจุด m อีก หนึ่งจุดที่อยู่ด้านนอกวงกลม โดยกำหนดให้ มุมที่ทำระหว่าง จุด l,m,n เท่ากับ 120 องศา และแขนของมุมทั้งสอง (เส้นตรงจากจุด m ไป n กับ เส้นตรงจาก จุด m ไป l) นี้ต้องเป็นเส้นสัมผัสวงกลม และก็ให้ เส้นตรง จากจุด m ไป l กับเส้นตรง จุด k n นี้ขนานกัน  555 ดูจะสับสนซับซ้อนอย่างค่อยเป็นค่อยไปนะครับ แต่ถ้าเห็น code และจะเข้าใจดีขึ้นครับ

สังเกตุนะครับว่าเราไม่ได้ระบุพิกัดของจุดใดๆเลย หลังจากที่อธิบายไดอะแกรมของเราเรียบร้อยแล้วผ่านคำสั่ง GeometricScene แล้วเราก็ให้ Mathematica สร้างภาพออกมาโดยใช้คำสั่ง RandomInstance ซึ่งมันก็จะสุ่มจุดขึ้นมาและวาดตามเงื่อนไขที่เรากำหนดไป ซึ่งถ้าเรารัน คำสั่ง RandomInstance และไม่กำหนด RandomSeeding ภาพวาดที่ออกมามันจะไม่ซ้ำกันเลยครับ ดังนั้นเราสามารถกำหนดหรือ fix ภาพ ได้จากการกำหนดค่า RandomSeeding ครับ

จากภาพที่วาดออกมาเราสามารถหาพื้นที่ของสามเหลี่ยน kln ได้ตามนี้ครับ

จากตัวอย่างข้างบนจะเห็นว่าผมใช้คำสั่ง GeometricAssertion เป็นตัวบอกว่า object ที่วาดขึ้นนั้น เช่น เส้นตรง วงกลม หรือจุดต่าง มันมีปฏิสัมพันธ์กันอย่างไร อย่างเช่น ผมต้องการ เส้นตรง kn กับเส้นตรง lm ขนานกัน ผมก็เพียงพิมพ์

GeometricAssertion[{Line[{k,n}],Line[{l,m}]},”Parallel”]

แต่ความเจ๋งของ Mathematica 12 ในเรื่อง Synthetic Geometry มันอยู่ตรง คำสั่ง FindGeometricConjecture นี้ครับ เจ้าคำสั่ง FindGemetricConjecture เป็นคำสั่งที่จะช่วยหาข้อสรุปต่างๆ จากตัวภาพที่วาดจากคำสั่ง GemetricScene ครับ เพื่อให้เห็นความสามารถของคำสั่งนี้ ผมจะลองให้มันหาดูว่ามุมภายในครึ่งวงกลมจะเป็นมุมฉากหรือเปล่าตามทฤษฎีของเธลิส เริ่มจากกำหนดจุด {a,b,c} อยู่บนวงกลมที่มีจุดศูนย์กลางที่จุด o และจุด o เป็นจุดกึ่งกลางระหว่างจุด a กับ c พอวาดเสร็จก็ให้หาข้อสรุปดู ตาม code ด้านล่างนี้

จะเห็นได้ว่าจาก คำสั่ง FindGeometricConjectures มันบอกว่าเส้นตรง ab กับ เส้นตรง cb มันตั้งฉากกัน ซึ่งมันก็โชว์อีกว่ามุม b ที่เกิดจากจุด abc มันทำมุม 90 องศา

ลองดูอีกสักตัวอย่าง เพื่อดูว่า จุด X,Y,Z ที่เกิดจากการตัดกันตามภาพล่างนี้อยู่บนระนาบเดียวกันหรือเปล่า

 ซึ่งจากผลที่ได้ X,Y,Z อยู่บนเส้นตรงเดียวกันครับ

ใครที่สนใจอยากดูตัวอย่างเพิ่มเติมลองไปดูได้ที่ https://www.wolfram.com/language/12/plane-geometry/

ซูมกราฟใน Mathematica

กราฟที่ได้จากคำสั่ง Plot ใน Mathematica นี่จะ click ซูมเข้าไปดูในบริเวณที่สนใจไม่ได้ ต้องคอยมากำหนด Range ที่จะ plot เอาเองเพื่อดูบริเวณที่สนใจ แต่ล่าสุดก็มีคนเขียนเป็น Function สำหรับทำให้เราสามารถที่จะเลือกซูมดูบริเวณที่สนใจในกราฟได้ครับ

ลองไปดูได้ที่ https://resources.wolframcloud.com/FunctionRepository/resources/DragZoomPlot

ผมลองโหลดมาเล่นแล้วนี่สะดวกมากๆครับ

ตัวอย่างใช้งาน notebook ของ Mathematica


ผมเขียน code เพื่อสร้างกราฟเพื่อดูว่าผลลัพธ์จากโมเดลที่ผมสร้างขึ้นมานั้นมันเหมือนหรือต่างข้อมูลจริงจากคนไข้อย่างไร โดยเขียนให้มันสร้าง notebook ขึ้นมาพร้อมกับรัน code สำหรับสร้างกราฟเลย โดยที่ในแต่ล่ะ notebook จะมีข้อมูลของคนไข้ประมาณ 300 คน ตามตัวอย่างนี้ครับ

ทีนี่ผมอยากจะ save กราฟที่สร้างขึ้นในแต่ล่ะ cell เป็น pdf จากไฟล์ notebook แต่ล่ะไฟล์ ซึ่งในแต่ล่ะ notebook จะมีประมาณ 5-10 cells ได้ ถ้าไปมัวเลือก save ในแต่ล่ะ cells คงต้องเสียเวลามาก เพราะมันมีหลายไฟล์ และผมก็มีหลายโมเดล

วิธีที่ผมทำก็คือเขียน code มาเปิด notebook แต่ล่ะไฟล์แล้วก็ให้มันเลือกเซฟเฉพาะ กราฟของแต่ล่ะ cell ที่ต้องการโดย code ที่เขียนก็จะประมาณนี้ครับ

ExportPdf[nbfilename_, outnames_List] := Module[{nb, nbr, i},
   nb = NotebookOpen[NotebookDirectory[] <> nbfilename];
   SelectionMove[nb, Before, Notebook];
   For[i = 1, i <= Length@outnames, i++,
    SelectionMove[nb, Next, CellContents, 2];
    nbr = NotebookRead[nb];
    If[nbr =!= {},
     Export[outnames[[i]] <> ".pdf", nbr];
     Print[outnames[[i]] <> ".pdf has been created."];
     ];
    ];
   NotebookClose[nb];
   ];

โดยคำสั่งหลักที่ใช้ก็คือ NotebookOpen, SelectionMove และ NotebookRead ครับ

NotebookOpen ก็เปิดไฟล์ notebook แล้วสร้างเป็น object จากนั้นก็เคลื่อนที่ภายใน Notebook ผ่านคำสั่ง SelectionMove พร้อมกับอ่านค่าตรงตำแหน่งที่เคลื่อนที่ไปด้วยคำสั่ง NotebookRead แล้วก็เขียนค่าที่ได้ออกมาเป็นไฟล์ใน format ที่ต้องการ

ลองดูเพิ่มเติมการใช้งาน notebook ของ Mathematica ได้ที่ https://reference.wolfram.com/language/tutorial/ManipulatingNotebooksOverview.html

คำแนะนำการใช้งาน Mathematica notebook

Notebooks จาก Mathematica นี่ปฏิวัติเรื่อง workflow ของงานด้านการคำนวณเลยที่เดียวครับ ไอเดียที่แสดงผลลัพธ์ต่อท้ายคำสั่ง สามารถที่เซฟและส่งต่อผลลัพธ์ในลักษณะเหมือนสมุดรายงานมีทั้งภาพ เสียง การเคลื่อนไหว แถมมีความสามารถเหมือน ms word กับ ms power point รวมกัน บวกกับสามารถที่จะ export เป็นได้หลายformat มาก และความสามารถอื่นๆอีกมากมาย ไม่แปลกใจที่ไอเดียจะถูกลอกไปในภาษาอื่นๆ โดยเฉพาะ R กับ Python แต่ความสามารถของ notebooks จากทั้งสองภาษายังห่างไกลจากของ Mathematica มาก

Video ที่ผมแชร์มานี้ก็เป็นเพียงตัวอย่างการใช้งาน notebook ของ Mathematica ให้เกิดประสิทธิภาพสูงสุด แต่ก็ยังมีอีกหลายความสามารถที่ไม่ได้พูดถึงครับ