CVODE Constants

ช่วงนี้ปวดหัวกับการ fit model ด้วย MCMC ที่ระบบมัน stiff มาก ในงานที่ทำนั้นใช้ตัว ode solver ที่เรียกว่า cvode ของ sundials ซึ่งมันก็ใช้ได้ดี แต่ chain มัน converge ช้ามากกกก แถมบางครั้ง cvode ก็จะมี flag แปลกๆ กลับมา ข้างล่างนี้เป็นความหมายของ flag ต่างๆ ซึ่งก็เอามาจาก cvode manual

เขาว่าโจทย์เลขป.4

มีคนแชร์โจทย์นี้ในpantip โดยจากโจทย์ก็ให้หาว่า อักษรแต่ล่ะตัวแทนด้วยเลขอะไรบ้าง

ดูแล้วเด็กก็คงต้องเขียนโปรแกรมช่วยแล้วล่ะครับ 555 ผมลองเขียนโปรแกรมด้วยภาษา Wolfram สำหรับปัญหานี้ตามนี้ครับ

ls = Permutations[Range[0, 9], {8}];

SetSharedFunction[ParallelSow]
ParallelSow[expr_] := Sow[expr]

output = (ParallelTable[
      If[(FromDigits@{a, f, b, f} + FromDigits@{c, g, h, b} + 
               FromDigits@{d, a, f, g} + FromDigits@{a, e, a, b} == 
              FromDigits@{b, c, d, c}) /. {a -> #1, b -> #2, c -> #3, 
             d -> #4, e -> #5, f -> #6, g -> #7, h -> #8} & @@ # &@
        ls[[i]]
       ,
       ParallelSow[ls[[i]]]
       ]
      , {i, 1, Length@ls}
      ] // Reap)[[2, 1]];

คำตอบที่ได้ก็จะประมาณนี้ครับ

หรือให้แสดงผลสวย ๆ

Table[
  #[[1]] <> "+" <> #[[2]] <> "+" <> #[[3]] <> "+" <> #[[4]] <> 
     "=" <> #[[5]] &@(
    ToString /@ ({FromDigits@{a, f, b, f}, FromDigits@{c, g, h, b}, 
            FromDigits@{d, a, f, g}, FromDigits@{a, e, a, b}, 
            FromDigits@{b, c, d, c}} /. {a -> #1, b -> #2, c -> #3, 
            d -> #4, e -> #5, f -> #6, g -> #7, h -> #8} & @@ # &@
       output[[i]])), {i, 1, Length@output}] // TableForm

ก็จะได้

ตัวอย่าง 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

แสดงผลลัพธ์จาก wolfram script แบบกราฟฟิค

ตั้งแต่ Wolfram Engine สามารถใช้ได้ฟรีนี่ทำให้ผมมีอะไรได้ลองเล่นเยอะเลย ถึงแม้จะมี Wolfram Mathematica อยู่แล้วนะเนี้ย

ตัว Wolfram Engine ที่เราสามารถโหลดมาใช้ได้ฟรี(http://www.wolfram.com/engine/) นี้จะเป็นแบบ text mode ไม่มี UI มาให้ ก็อาจจะทำให้ลำบากหน่อยถ้าอยากจะดูผลลัพธ์ที่เป็นกราฟฟิค เพราะอาจต้อง Export เป็นไฟล์ไปก่อนแล้วค่อยเปิดดู แต่ถ้าอยากจะให้ตัว Wolfram Engine มันแสดงผลกราฟฟิคให้ดูเราก็สามารถที่จะใช้ตัวแพ็กเกจ JavaGraphics ช่วยได้ครับวิธีก็เพียงโหลด

<<JavaGraphics`

ก่อนที่จะให้มันแสดงผลลัพธ์กราฟฟิคของคำสั่งที่ต้องการ ตัวอย่างครับ

หรือถ้าใครสนใจจะใช้ผ่าน UI อย่าง jupyter ก็ลองดูที่ projects ตามนี้ดูนะครับ

สวนตัวแนะนำ WolframLanguageForJupyter ครับ

วิธการติดตั้งที่อยากแนะนำก็คือลง Anaconda พอลงเสร็จก็รัน Anaconda prompt แล้วรัน wolframscript

จากนั้นก็ติดตั้ง WolframLanguageForJupyter ตามวิธีที่เขาแนะนำไว้ครับ

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