maemod package

จากที่เคยเขียนไว้เกี่ยวกับ r package อันหนึ่งที่ผมเขียนชื่อ maemod สำหรับช่วยให้คนที่สนใจอยากคำนวณพวก ode ได้ง่ายขึ้น (ดูข้างล่าง)

maemod (แม่มด)

มีคนสนใจว่าถ้าพวกตัวแปร state ต่างๆนั้นเป็นแบบ array จะทำอย่างไร ผมเลยเขียนตัวอย่างพร้อมกับเพิ่มความสามารถด้านarrayนี้เข้าไป พอใช้ได้ไปคร่าวๆก่อน ดูตัวอย่างข้างล่างนี้ครับ

# Example from Berkeley Madonna
 # for using Array
 # METHOD RK4
 #
 # STARTTIME = 0
 # STOPTIME = 20
 # DT = 0.02
 #
 # d/dt (A[1]) = -k[1]*A[1]
 # d/dt (A[2..n-1]) = k[i-1]*A[i-1]-k[i]*A[i]
 # d/dt (A[n]) = k[n-1]*A[n-1]
 #
 # init A[1] = 100
 # init A[2..n]=0
 #
 # k[1..n] =1
 # n = 14
 #
 #

arrayEx <-'
 !ExtraFunctions
 n<-14

!Equations
 dA <-rep(0,n)
 dA[1] <- -k[1]*A[1]
 #ArrayB
 "dA[i.indx] <- k[i.indx-1] *
 A[i.indx-1]-k[i.indx]*A[i.indx]"%@@%list("i.indx"%=>%2:(n-1))
 #ArrayE
 dA[n] <- k[n-1]*A[n-1]

!Parameters
 1,1,1,1,1,1,1,1,1,1,1

!Inits
 100,0,0,0,0,0,0,0,0,0,0,0,0,0

!Outputs
 c(dA)

!Plots

!MAEMOD_End
 '
 out<-maemod.ode(input.text = arrayEx ,timegrid = seq(0,20,0.02),sys.template = Maemod_Array)
 plot(out[,c(1,2)],type = 'l')
 cols<-rainbow(14)
 for(i in 3:n) lines(out[,c(1,i)],col=cols[i])

ในเบื้องต้นนี้ ตัวstate กับ parameter จะต้องเป็น A กับ k เท่านั้น จากตัวอย่างเราสามารถที่จะdefine ตัวแปรที่ index ต่างได้ เช่น ต้องการ

d/dt (A[2..n-1]) = k[i-1]*A[i-1]-k[i]*A[i]

โดยที่ i มีค่าตั้งแต่ 2 ถึง  n

เราก็สามารถพิมพ์ตามนี้ได้

#ArrayB

"dA[i] <- k[i-1]*A[i-1]-k[i]*A[i] "%@@%list('i'%=>%2:n)

#ArrayE

เจ้าเครื่องหมาย %@@% เป็นตัวบอกว่าเราจะใช้ array กับ สมการด้านซ้ายมือของเคื่องหมายและindexจะอยู่ด้านขวาของเครื่องหมาย ซึ่งจะต้องเป็นตัวแปรแบบ list

ส่วนเครื่องหมาย %=>% บอกว่าสัญลักษณ์หรือรูปแแบบindexด้านขวาจะมีค่าเปลี่ยนไปตาม vector ด้านซ้าย เช่น ‘i’%=>%2:n จะหมายถึง i มีค่าตั้งแต่ 2 ถึง n

ส่วน #ArrayB และ #ArrayE เป็น block ที่บอกว่ามีการใช้ ตัวแปรแบบ array ในข้อความที่อยู่ใน block นี้ครับ

codeที่มีarray blockนี้เวลาเอาไปใช้กับ maemod.ode จะต้องใช้กับtemplate เฉพาะที่ชื่อ Maemod_Array นั่นก็คือต้องเซ็ต option ชื่อ sys.template = Maemod_Array ครับ

สนใจอยากลองใช้งานก็ดูได้ที่

https://github.com/slphyx/maemod

 

ใช้ประโยชน์จาก rstan

ถ้าใครได้ใช้ deSolve สำหรับแก้ปัญหาเชิงตัวเลขของระบบสมการเชิงอนุพันธ์(ode) แล้วจะเห็นว่ามันมีข้อจำกัดเรื่องความเร็วอยู่ เพราะถ้าเขียน function ของสมการด้วย R แล้วจะช้ามาก ยิ่งถ้ามีการเอาfunction ไปใช้ fitting ด้วยเทคนิคอย่าง Markov chain monte carlo แล้วนี่เลิกคิดได้เลยสำหรับระบบใหญ่ ๆ

deSolve เองก็พยายามแก้ปัญหานี้ โดยผู้ใช้สามารถเรียกใช้ function ที่เขียนใน C/C++ ที่คอมไพล์แล้วมาใช้ได้ แต่ก็ยังดูยากพอสมควรสำหรับคนที่ไม่เคยเขียน C/C++ เลย แต่มันก็ยังพอมีวิธีที่ดีกว่าหน่อยที่อยากจะแนะนำสำหรับคนที่อยากจะเพิ่มความเร็วสำหรับใช้กับ deSolve  นั่นก็คือเขียน ในภาษาของ rstan แทนครับ

โดยหลักการก็คือเขียนสมการใน rstan แล้วใช้คำสั่งของ rstan ที่เรียกว่า expose_stan_functions เพื่อที่จะคอมไพล์ function แล้วเราก็จะยังเรียกใช้ได้ใน .GlobalEnv นั่นหมายความว่าเราสามารถเรียกใช้ function ที่คอมไพล์แล้วนี้กับ deSolve ได้ครับ

ตัวอย่างเช่น ผมเขียนสมการ lorenz ในไฟล์ชื่อ lorenz.stan

##lorenz.stan
functions{
  vector lorenz(real t,vector y,vector parms){
    vector[3] ydot;
    // parms[1]=beta ,[2]=sigma,[3]=rho
    
    ydot[1]=parms[2]*(y[2]-y[1]);
    ydot[2]=y[1]*(parms[3]-y[3])-y[2];
    ydot[3]=y[1]*y[2]-parms[1]*y[3];

    return ydot;
  }
}
model{
}

จากนั้นก็ใช้ expose_stan_functions เรียกจากใน R ได้เลยตามนี้ครับ

library(rstan)

expose_stan_functions("lorenz.stan")

# parms[1]=beta ,[2]=sigma,[3]=rho
parameters1<-c(2.66666,10,28)

# x=-10, y=-12, z=30.05
init.state<-c(-10,-12,30.05)

nderivs<-function(t,state,parameters){
  tmp<-lorenz(t,state,parameters)
  return(list(tmp))
}

run1<-ode(y=init.state,times=seq(0,15,0.01), 
   func=nderivs, parms=parameters1)

plot(run1)