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

แสดง output จาก ggplot ใน Mathematica

ลองดูตัวอย่างที่ผมเขียนนี้ครับ  ไอเดียก็คือว่าให้มันเขียนภาพเป็นไฟล์ pdf ด้วยคำสั่ง ggsave แล้วเรียกกลับมาใน Mathematica ด้วย Import  โค้ดที่เขียนนี้ดัดแปลงมาจาก https://mathematica.stackexchange.com/questions/16726/how-to-get-a-plot-generated-by-r-returned-in-an-output-cell

 

 

 

RDotNet 1.7 ใช้งานได้กับ R.3.3-4 เท่านั้น

ผมเขียนโปรแกรมเล็กๆที่เรียกใช้งาน RDotNet ไว้หลายตัว ซึ่งมันก็ทำให้สะดวกมากๆในการที่จะเรียก R มาใช้งานใน C# ซึ่งก่อนหน้านี้จะใช้งานกับ R ที่ค่อนข้างเก่าหน่อยพวก R 3.3-4 แต่พอมี library บางตัวที่บังคับให้ต้อง upgrade ตัว R ไป version ที่สูงขึ้นอย่าง 3.5 นี้ทำให้โปรแกรมที่เขียนหลายตัวมันมีปัญหาทันที เหตุผลก็เพราะมัน GetInstance() ไม่ได้ มีการพูดถึงปัญหานี้ที่ https://github.com/jmp75/rdotnet/issues/70 ตอนนี้ที่ทำได้ก็คือถ้า library ที่อยากใช้นั้นถ้าไม่ได้ใช้ features ที่มันมีมาใหม่นั้นจริงๆก็ต้องใช้ R ตัวเก่าไปก่อน

 

RcppGSL บน Windows 10

ใครที่จะใช้ RcppGSL บน Windows 10 อาจจะมีเรื่องที่ต้องทำความเข้าใจหน่อยนะครับ มีคนถามผมมาว่ามันเซ็ตอย่างไร เพราะลองเอา code ตัวอย่าง มาใช้แล้วก็ยังไม่ผ่านสักที และดูเหมือนว่าคนถามผมเองยังสับสนและไม่เข้าอยู่เล็กน้อย  เล่าคร่าวๆ ก่อนนะครับ เจ้า GSL หรือ GNU Scientific library เนี้ยเป็นlibraryที่รวบรวมคำสั่งหรือ functions ที่ใช้กันมากในวงการคณิตศาสตร์และฟิสิกส์ มันถูกเขียนขึ้นด้วยภาษา C โดยกลุ่มนักฟิสิกส์และแจกจ่ายให้ใช้กันได้ฟรีครับ

ที่นี่มันก็มีคนอยากใช้งานจาก GSL นี้ใน R ก็เลยมีคนเขียน package ขึ้นมากันหลายตัว และ  RcppGSL ก็เป็นหนึ่งในนั้น เจ้า RcppGSL นี้ถ้าเราจะใช้มันเนี้ย เราต้องมี GSL ติดตั้งให้เรียบร้อยก่อน (เดี๋ยวบอกอีกทีว่าทำอย่างไรครับ) โดยที่เราจะต้องเซ็ตพารามิเตอร์สำหรับ Windows ที่ชื่อว่า LIB_GSL ไปยัง path ของ GSL ที่ติดตั้งไว้ครับ

อย่างเช่นผมติดตั้งไว้ที่ D:\Data-Work\Programs\gsl\gsl_2.5 ผมก็เพียงเพิ่มเข้าใน environment variables ตามนี้ครับ สังเกตุว่ามันมี “” ด้วยนะครับ

จากนั้นก็เปิด R หรือ RStudio แล้วลองใช้คำสั่ง Rcpp::sournceCpp กับไฟล์ที่ต้องการอีกที่นะครับ

อีกเรื่องที่อาจต้องเพิ่มเข้าไปใช้ code คือบรรทัดนี้ครับ

// [[Rcpp::depends(RcppGSL)]]

 

Rcpp::sourceCpp กับปัญหา -Wunused-variable

ใครที่ใช้ sourceCpp สำหรับ compile โค้ดของ Rcpp แล้วเจอปัญหา unsed-variable แล้วอยากที่จะ  ignore มัน ก็สามารถทำได้ด้วยการสร้างไฟล์ Makevars ที่มีบรรทัดนี้ครับ

CXXFLAGS += -O3 -Wall -pipe -Wno-unused

หรือไม่ก็แก้ไขไฟล์ Makeconf ที่ R_HOME/bin/R_ARCH (ไม่ค่อยแนะนำ) โดยเพิ่ม -O3 -Wall -pipe -Wno-unused  ในส่วนของ CXXFLAGS  ครับ

 

 

ปัญหา unlock_solver จาก แพ็คเกจ deSolve

ถ้าใครที่เจอปัญหา unlock solver ตอนที่ใช้คำสั่งพวก ode solver จาก แพ็คเกจ deSolve เช่น

Error in .C(“unlock_solver”) :
“unlock_solver” not resolved from current namespace (deSolve)

 

ให้ลองแก้ปัญหาโดยโหลด deSolve ด้วยคำสั่งนี้แทนครับ

library.dynam.unload(“deSolve”, libpath=paste(.libPaths()[1], “//deSolve”, sep=””))
library.dynam(“deSolve”, package=”deSolve”, lib.loc=.libPaths()[1])

 

หรือไม่ก็ลองโหลด deSolve ก่อนที่จะเรียกใช้งานพวก shared library (.dll, .so) ครับ

 

 

 

a simple SIR simulation

ผมเขียน shiny app ง่ายๆเอาไว้สอนเรื่องโมเดลโรคระบาดแบบง่ายๆที่เรียกว่า SIR ครับ

code สามารถดูได้ที่ https://github.com/slphyx/HatGame

ใช้งาน gcc ที่มากับ Rtools ใน Mathematica

ใน Rtools จะมี compiler ของ gcc (mingw) มาด้วยแล้วทั้งที่เป็นแบบ 32 bits และ 64 bits หากอยากจะเอาไปใช้ใน Mathematica ก็ต้องเรียกผ่าน CCompilerDriver`GenericCCompiler` ครับ แล้วเพียงเซ็ต Path ของ gcc จาก Rtools นี้ให้ถูก เช่น

Needs["CCompilerDriver`GenericCCompiler`"]

ทดลองเรียกใช้งาน

greeter = CreateExecutable[StringJoin[
 "#include <stdio.h>\n",
 "int main(){\n",
 " printf(\"Hello MinGW-w64 world.\\n\");\n",
 "}\n"],
 "helloworld", "Compiler" -> GenericCCompiler, 
 "CompilerInstallation" -> "C:/Rtools/mingw_64", 
 "CompilerName" -> "x86_64-w64-mingw32-gcc.exe"]

Import["!\""<>greeter<>"\"","Text"]
Hello MinGW-w64 world.

ที่นี้ถ้าอยากจะให้ Mathematica มันเรียกใช้เจ้า mingw 64 จาก Rtools นี้ตลอด อย่างเช่นคำสั่ง Compile ก็สามารถทำได้โดยเซ็ตค่าที่มันเกี่ยวข้องอย่างเช่น Path กับตัวแปร $CCompiler ได้เลยครับ

$CCompiler = {"Compiler" -> GenericCCompiler, 
 "CompilerInstallation" -> "C:/Rtools/mingw_64", 
 "CompilerName" -> "x86_64-w64-mingw32-gcc.exe"};

f = Compile[{x, y}, Sqrt[x^2 + y^2], CompilationTarget -> "C"]

Table[{x, f[x, 5/4 x]}, {x, 0, 6, 0.5}]

 

Inkscape กับ ggplot

เสียเวลาพอสมควรที่ต้องมาคอยจัดระยะห่างในเอกสารโดยเฉพาะกราฟต่างๆที่ต้องนำมาวางให้อยู่ด้วยกันเป็น grid หรือตาราง โดยใน R นั้นคำสั่งที่ผมใช้เยอะก็พวก คำสั่ง plot_grid หรือแม้แต่ option ของ ggplot อย่าง plot.margin นั้นบางที่มันก็ทำให้กะขนาดเอายากพอสมควร นี่ยังไม่รวมเรื่อง scale หรือขนาดของ fonts อีก …. โอย….

วิธีที่พอช่วยเรื่องพวกนี้ได้บ้างก็คือ หลังจากที่ export เป็น pdf ของกราฟที่ต้องการจาก plot_grid แล้วก็อาศัย Inkscape มาช่วยในการขยับกราฟให้อยู่ในตำแหน่งที่ต้องการ วิธีการก็แสนจะง่าย เพียงแค่เปิดไฟล์ pdf ใน Inkscape จากนั้นก็ click หรือลากกราฟที่จะย้าย หรือจะเพิ่มเติมอะไรลงไปเลยยังได้

แนะนำให้ติดตั้ง Inkscape ไว้ที่เครื่องครับ อีกตัวก็เป็น sumatra pdf ครับมันจะเร็วดีครับในการเปิดเอกสาร pdf แต่จะแก้ไขเพิ่มเติมอะไรไม่ได้ครับ

 

เปลี่ยนรูปแบบของตัวเลขเป็นแบบ m*10^n บนแกน ใน ggplot

มีตัวอย่างที่ code ที่นำมาใช้ได้ที่ https://groups.google.com/forum/#!topic/ggplot2/a_xhMoQyxZ4 (โดย Brian Diggs)

fancy_scientific <- function(l) {
 # turn in to character string in scientific notation
 l <- format(l, scientific = TRUE)
 l <- gsub("0e\\+00","0",l)
 # quote the part before the exponent to keep all the digits
 l <- gsub("^(.*)e", "'\\1'e", l)
 # turn the 'e+' into plotmath format
 l <- gsub("e", "%*%10^", l)
 # return this as an expression
 parse(text=l)
}

เอาไปใช้งาน

library(ggplot2) 
library(scale)
 ggplot(data=df, aes(x=x, y=y)) +
 geom_point() +
 scale_y_continuous(labels=fancy_scientific)