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)

 

 

maemod (แม่มด)

ผมเขียน R package ชื่อ maemod สำหรับช่วยให้คนเพิ่งเริ่มเรียนรู้ใช้งาน R และอยากจะแก้สมการเชิงอนุพันธ์เบื้องต้น (ordinary differential equation: ode) ใน R ให้ทำได้ง่ายขึ้นครับ โดยที่ผู้ใช้งานก็เพียงเขียนสมการ ODE (ในรูปแบบคำสั่งของ R), ค่าของตัวแปรพารามิเตอร์ และค่าเริ่มต้นต่าง ๆ ของสมการ ตามตัวอย่างนี้ครับ
1. สร้างไฟล์สมมุติชื่อ test.txt ที่มีสมการตามนี้ครับ

!Equations
 dY1<-Y2
 dY2<-a*sin(Y1)+sin(t)
 
!Parameters
 a= -1.0
 
!Inits
 Y1=0,Y2=0
 
!Outputs
 c(dY1,dY2)
 
!ExtraFunctions

!Plots
!MAEMOD_End

โดยที่บรรทัดที่อยู่ถัดจาก !Equations คือสมการ ode, !Parameters คือค่าของตัวแปรต่าง ๆ, !Inits ค่าเริ่มต้นของแต่ล่ะสมการ ode, !Outputs คือlistของoutputที่ต้องการโดยที่คอลัมน์แรกซึ่งคือค่าของเวลาจะถูกเพิ่มเข้าไปอัตโนมัติเสมอ, !ExtraFunctions ในกรณีที่มี function หรือตัวแปรพิเศษก็สามารถใส่ในส่วนนี้ได้ครับ และสุดท้าย !MAEMOD_End เป็นตัวบอกว่าสมการหรือค่าต่างๆสิ้นสุดที่บรรทัดนี้ครับ พูดง่ายๆก็คือมีคีย์เวิร์ด อยู่ 6 คำคือ !Equations, !Parameters, !Inits, !Outputs, !ExtraFunctions และ !MAEMOD_End ครับ ไม่จำเป็นต้องเรียงตามลำดับนี้นะครับ อันจะขึ้นก่อนหลังก็ได้
2. วิธีรันก็เพียงใช้คำสั่ง maemod.ode เช่น
ต้องการจะคำนวณสมการข้างบนโดยที่ ค่า t เริ่มต้นที่ 0 เพิ่มขึ้นที่ล่ะ 0.1 จนถึง 20*pi

out<-maemod.ode(input.filename = "test.txt",timegrid = seq(0,20*pi,0.1))
head(out)  time Y1 Y2 
[1,] 0.0 0.0000000000 0.000000000 
[2,] 0.1 0.0001663662 0.004991383 
[3,] 0.2 0.0013277223 0.019866019 
[4,] 0.3 0.0044590899 0.044327146 
[5,] 0.4 0.0104962890 0.077883467 
[6,] 0.5 0.0203162441 0.119856470
plot(out[,c(2,3)],col="red")

heart

อีกตัวอย่างครับ สำหรับแก้สมการลอเร้นซ์   แต่แทนที่จะเซฟสมการและพารามิเตอร์ใน text file เราสามารถสร้างเป็นตัวแปรเก็บข้อความของระบบที่มีตัวแปรต่างๆได้ตามนี้เลยครับ (สังเกตุตัวแปรที่ชื่อ Lorenz ครับ)

library(scatterplot3d)
library(maemod)
lorenz<-"
 !Equations
 dx<-sigma*(y-x)
 dy<-x*(rho-z)-y
 dz<-x*y-beta*z
 
 !Parameters
 sigma=10, beta=2.66666, rho=28
 
 !Inits
 x=-10, y=-12, z=30.05
 
 !Outputs
 c(dx,dy,dz)

 !Plots
 !ExtraFunctions
 
 !MAEMOD_End
 "
lorenzout<-maemod.ode(input.text = lorenz, timegrid = seq(0,30,0.01))

scatterplot3d(lorenzout[,c(2,3,4)],type = 'l')

lorenzถ้าสนใจอยากติดตั้ง package นี้ก็สามารถทำได้ตามนี้เลยครับ

 library(devtools)
 install_github("slphyx/maemod")

code ทั้งหมดอยู่ที่ https://github.com/slphyx/maemod ครับ

จริงๆแล้วใน maemod นี้เราสามารถเลือกใช้ methods ต่างๆ เช่น lsoda, rk4 และอื่นๆอีกมากมาย สำหรับแก้ระบบสมการอนุพันธ์นี้ได้ครับ โดยวิธีการเลือกใช้ methods เหล่านี้สามารถดูได้ที่ package deSolve (http://desolve.r-forge.r-project.org/index.html) เลยครับเพราะค่าพารามิเตอร์สำหรับmethods ต่างๆนี้จะถูกส่งไปเรียกจาก deSolve อีกทีครับ

หากอยากช่วยเขียนปรับปรุงก็ Pull requests หรือ Fork ได้เลยครับ