Categories
Mathematica

DateListPlot

มีคนถามมาว่ามีข้อมูลตามมด้านล่างนี้จะเอาไป plot ได้ยังไง

20110508000128    1
20110508000315    5
…………. .
20110508235824    x

ซึ่งใน column แรกข้อมูลคือวันที่และเวลาที่อยู่ติดกัน และ column ที่สองก็คือข้อมูลจากการวัดอะไรสักอย่าง

ผมก็บอกคนที่ถามไปว่าไม่ยาก:)  ใน Mathematica ก็น่าจะทำได้

ก่อนอื่นเลยก็ต้องทำการจัดรูปแบบของข้อมูลให้ตรงตามที่ Mathematica ต้องการก่อน โดยที่ผมคิดว่าจะเอาไปใช้กับ

function ชื่อ DateListPlot ครับ ดูรายละเอียดได้ที่ http://reference.wolfram.com/mathematica/ref/DateListPlot.html

อันนี้เป็น function ที่ผมเขียนขึ้นมาใช้ในการเปลี่ยนรูปแบบของข้อมูลครับ

convertdata[yourdatapoint_] := Module[
{output, ls, dat1, dat2, year, month, date, hour, minute, second},
dat1 = yourdatapoint[[1]];
dat2 = yourdatapoint[[2]];
ls = Characters@ToString@dat1;
year = ToExpression@StringJoin@ls[[1 ;; 4]];
month = ToExpression@StringJoin@ls[[5 ;; 6]];
date = ToExpression@StringJoin@ls[[7 ;; 8]];
hour = ToExpression@StringJoin@ls[[9 ;; 10]];
minute = ToExpression@StringJoin@ls[[11 ;; 12]];
second = ToExpression@StringJoin@ls[[13 ;; 14]];
output = {{year, month, date, hour, minute, second}, dat2};
output
];

หรือจะเขียนแบบนี้โดยใช้ StringTake ก็ได้เหมือนกันครับ

convertdata[yourdatapoint_] :=
Module[{output, ls, dat1, dat2, year, month, date, hour, minute,
second},
dat1 = yourdatapoint[[1]];
dat2 = yourdatapoint[[2]];
ls = ToString@dat1;
{year, month, date, hour, minute, second} =
ToExpression /@
StringTake[
ls, {{1, 4}, {5, 6}, {7, 8}, {9, 10}, {11, 12}, {13, 14}}];
output = {{year, month, date, hour, minute, second}, dat2};
output
];

เอาไปใช้กับข้อมูลก็สามารถทำได้แบบนี้ครับ

data = {{20110508000128, 1}, {20110508000315, 5}, {20110508000411,
10}, {20110508001130, 20}};

DateListPlot[convertdata2 /@dat a, Joined -> True]

 

Categories
Mathematica Uncategorized

วงกลมแห่งความตาย

จากที่มีคนถามปัญหาที่ http://www.pantip.com/cafe/wahkor/topic/X10570620/X10570620.html

มีชายหนึ่งล้านคนยืนกันเป็นวงกลม

ชายเหล่านี้ตั้งใจ จะฆ่าตัวตายหมู่ แต่เนื่องจากไม่มีใคร
กล้าฆ่าตัวตายจึงตกลงทำดังนี้ ให้ทุกคนล้อมเป็นวงกลม
จากนั้นให้เริ่มต้นที่คนแรก นำดาบฆ่าคนที่ 2
แล้วส่งดาบไปให้คนที่ 3 แล้วคนที่ 3 นำดาบนั้นฆ่าคนที่ 4
แล้วส่งดาบไปให้คนถัดไป ทำเช่นนี้ จนเหลือผู้รอดชีวิตอยู่คนเดียว
จึงค่อยฆ่าตัวตาย
ถ้าคุณบังเอิญอยู่ในกลุ่มคนพวกนี้ และเกิดไม่อยากถูกใครฆ่าขึ้นมา
คุณจึงพยายามไปยืนตำแหน่งที่จะเหลือรอดเป็นคนสุดท้าย
ถามว่าคุณจะไปยืนเป็นคนที่เท่าไหร่

 

เพื่อที่จะแก้ปัญหานี้แบบคนขี้เกียจผมก็ลองเขียนโปรแกรมโดยใช้ Mathematica ดูปรากฏว่าคำตอบที่ได้คือ

ต้องยืนที่ตำแหน่ง 951425

อันนี้โปรแกรมที่เขียนครับ

drp = Compile[{{ls, _Integer, 1}},
Select[ls, MemberQ[Drop[ls, {1, Length@ls, 2}], #] == False &]]

fn[ls_] :=
Module[{tmp}, tmp = drp[ls];(*Select[ls,MemberQ[Drop[ls,{1,Length@ls,
2}],#]==False&];*)
If[EvenQ[Length@ls], Developer`ToPackedArray@tmp,
Developer`ToPackedArray@RotateRight@tmp]
]

ls = Developer`ToPackedArray[Table[i, {i, 1, 10^6}]];

Nest[fn2, ls, 30]

{951425}

 

ซึ่งจริงๆแล้วปัญหานี้ก็คือรูปแบบหนึ่งของปัญหาที่เรียกว่า Josephus problem ครับ เราสามารถเขียนด้วย Mathematica

สั้นๆได้ตามนี้เลยครับ

Needs[“Combinatorica`”]
Last@InversePermutation[Josephus[10^6, 2]]
951425

 

ปล. อาจจะคำนวณนานหน่อยนะครับ 🙂

Categories
Molecular Dynamics

มาลองทำ Molecular Dynamics Simulation ของโปรตีนกันครับ 2

ต่อจากตอนที่แล้ว (http://www.sakngoi.com/?p=189)ครับ

Molecular dynamics (MD) ก็เป็นเทคนิคที่ฮิตฮอตสำคัญอันหนึ่งที่นักวิทยาศาสตร์ใช้ในการศึกษาเกี่ยวกับการเคลื่อนไหวหรือสั่นไหวของโมเลกุลต่างๆ โดยเค้าจะมองว่าพันธะที่คอยเชื่อมต่ออะตอมเข้าด้วยกันเป็นโมเลกุลนี้เป็นเหมือนสปริงที่เชื่อมลูกบอลเข้าไว้ด้วยกัน  ซึ่งพลังงานศักย์ทั้งหมดของระบบก็จะมาจากผลรวมทั้งหมดของ interaction ระหว่างลูกบอลที่เชื่อมกันผ่านสปริง(bond) และระหว่างสปริงด้วยกัน(non-bonded) ถ้าใครที่เรียนมาก็คงจะจำได้ว่าแรงของปริงนี้เป็นแรงอนุรักษ์ซึ่งเราจะเขียนได้ว่าแรงทั้งหมดที่กระทำกับลูกบอล(อะตอม) i นี้เขียนได้ว่า
F_i=-\nabla E

โดย E ก็คือพลังงานศักย์ ซึ่งเราเราสามารถเขียนได้ว่า E_{total}=E_{bond}+E_{non-bonded} โดยเทอมที่ไม่เป็นพันธะสามารถเขียนได้ว่าเป็นผลรวมศักย์จากแรงวันเดอร์วาลส์กับแรงจากประจุไฟฟ้า E_{non-bond}=E_{vdW}+E_{elec} โดย E_{elec}=\sum_{i=1}^{N} \sum_{j=i+1}^{N}\frac{q_iq_j}{4\pi\varepsilon_0 r_{ij}} และ E_{vdW}= \sum_{i=1}^{N} \sum_{j=i+1}^{N}4\varepsilon\left(\frac{\sigma_{ij}^{12}}{r_{ij}^{12}}-\frac{\sigma_{ij}^{6}}{r_{ij}^{6}}\right)

ส่วนพลังงานศักย์ของแรงจากพันธะก็สามารถเขียนได้ว่าเป็นรวมจากรูปแบบของพันธะที่เป็นไปได้ทั้งหมด (ดูรูป เดี๋ยวมาใส่ให้ครับ)

E_{bond}=E_{stretch}+E_{angle}+E_{out-of-plane}+E_{str-str}+...
หรือ
E_{bond}=\frac{1}{2}\sum_\text{angles}k_a\left(\theta-\theta_0\right)^2+\frac{1}{2}\sum_\text{bonds}k_b\left(l-l_0\right)^2+

\sum_\text{bonds}k_d\left(1 + \cos\left(n\phi-\phi_0\right)\right)+\sum_{\text{improper}}k_i\left(\omega-\omega_0\right)^2
สมการที่เขียนมานี่ก็เป็นเพียงตัวอย่างเท่านั้นครับซึ่งมันจะแตกต่างกันไป ขึ้นกับว่าเราเป็นสาวกของสำนักไหน 🙂 แต่ล่ะสำนักก็จะมีรูปแบบของฟังก์ชั่นและค่า constants ต่างๆที่ใช้แตกต่างกัน (เราเรียก พวกฟังก์ชั่นนี้ว่า force field)
force field ที่เป็นที่นิยมกันก็ได้แก่

AMBER (Assisted Model Building with Energy Refinement, http://ambermd.org/)

OPLS, OPLS-AA (Optimized Potentials for Liquid Simulations,http://zarbi.chem.yale.edu/)

CHARMM (Chemistry at HARvard Macromolecular Mechanic, http://mackerell.umaryland.edu/MacKerell_Lab.html)

GROMOS (GROningen Molecular Simulation, http://www.gromos.net/)

โดยปกติเค้าจะไม่เอาแต่ล่ะ force field มาใช้ปนกันครับ

มาต่อที่เรื่องของสมการนิวตันที่เราได้ครับ F_i=m_i\frac{d^2r_i}{dt^2}=-\nabla E  ซึ่งเป็นสมการที่แสดงความสัมพันธ์ของพลังงานศักย์กับระยะทางหรือตำแหน่งของอะตอมที่เราสนใจ เราสามารถแก้สมการนี้โดยอาศัยวิธีที่เรียกว่า Verlet หรือ Leap-fr0g ใครที่เรียนเรื่อง numerical analysis มาคงเข้าใจนะครับ เพราะมันก็เป็นวิธีการอินทีเกรตทำธรรมดา แตกกันก็เพียงว่าจะเอากี่เทอม เริ่มที่ตรงไหนก่อนเท่านั้นเอง

ใครที่สนใจเรื่องการเขียนโปรแกรมเกี่ยวกับ MD นี้ลองดูที่นี่ครับ http://www.sakngoi.com/?p=64

ต่อตอนที่ 3

Categories
Molecular Dynamics

มาลองทำ Molecular Dynamics Simulation ของโปรตีนกันครับ 1

มาลองทำ Molecular Dynamics Simulation ของโปรตีนกันครับ

โปรแกรมที่ผมใช้ในการคำนวณคือ NAMD ส่วนที่เอาไว้ใช้ดูโครงสร้างสวยๆของโปรตีนคือ VMD

ทั้งสองโปรแกรมนี้ฟรีครับ

โปรตีนที่ผมพูดถึงนี้ก็คือ HIV 1-RT (reverse transcriptase) ครับ เราจะมาลองดูครับว่าเจ้าสาย DNA หรือโครงสร้างของโปรตีนแถวๆ Active site ของมันจะมีการเปลี่ยนแปลงไปยังไงบ้างถ้าเราใส่ตัว Inhibitor หรือยาอย่าง adefovir หรือ tenofovir เข้าไป โครงสร้างของโปรตีนที่ผมใช้คือ 1T05 ซึ่งสามารถdownload ได้จากเวบ www.pdb.org

รูปด้านล่างนี้แสดงลำดับของนิวคลีโอไทด์กับตำแหน่งของกรดอะมิโนจากโปรตีนที่ผมว่าจะเอาไว้ดูว่ามันเจ้าสายDNA กับ HIV-RT มันมีอะไรกันหรือเปล่า จริงๆแล้วมันก็คือสะพานเกลือครับ 🙂 เรียกตรงๆเลย ผมสนใจว่าเจ้าฟอสเฟตจากสาย DNA กับพวกขั้วบวกจากกรดอะมิโนอย่างไลซีนกับอาร์จินินมันจะแนบแน่นขนาดไหนตลอดเวลาหรือเปล่า ซึ่งผมก็หาอยู่ตั้งนานก็เจอตัวที่ใกล้ที่กันมากที่สุดก็ตามรูปเลยครับ

ส่วนเจ้าสาย DNA เราก็จะดูว่า furanose ring มันมีการผับตัวหรือบิดยังไง (sugar puckering) ในตำแหน่งต่างๆ

แต่ก่อนจะไปดูว่าเราจะทำ Molecular dynamics ได้อย่างไรเราก็ควรจะเข้าใจก่อนสักนิดหนึ่งว่า Molecular dynamics คืออะไรซึ่งเดี๋ยวผมจะเล่าให้ฟังคร่าวๆนะครับส่วนใครที่สนรายละเอียดทางด้านทฤษฎีของ Molecular dynamics ผมแนะนำให้อ่านหนังสือตามนี้เลยครับ

The Art of Molecular Dynamics Simulation by D. C. Rapaport

Understanding Molecular Simulation, Second Edition: From Algorithms to Applications (Computational Science) by Daan Frenkel and Berend Smit

Molecular Modeling and Simulation by Tamar Schlick

Molecular Modelling: Principles and Applications (2nd Edition) by Andrew R. Leach

Computer Simulation of Liquids by M. P. Allen and D. J. Tildesley

ทำไมถึงต้องเป็นห้าเล่มนี้ ง่ายๆเลยครับเพราะผมมีแค่ห้าเล่มนี้เอง 🙂

ต่อตอน 2

Categories
Mathematica

switching between notebooks in Mathematica

ถ้าจะสวิชระหว่าง notebooks ใน Mathematica กด Ctrl+F6

shortcut อื่นๆดูได้ที่
http://reference.wolfram.com/mathematica/tutorial/KeyboardShortcutListing.html

Categories
Mathematica

การใส่ไฟล์ CDF ใน HTML

เครื่องที่จะเปิดดูไฟล์ต้องลง CDF player ก่อนนะครับ

ถ้าจะ embed ไฟล์ cdf ลงใน html สามารถทำได้ตามนี้ครับ

<embed src="ไฟล์.cdf" width="588" height="380">

หรือ

<object classid="clsid:612AB921-E294-41AA-8E98-87E7E057EF33"
width="500" height="300"
type="application/vnd.wolfram.cdf.text">
<param name="src" value="ไฟล์.cdf">
<embed width="500" height="300" src="ไฟล์.cdf"
type="application/vnd.wolfram.cdf.text">
</object>

หรือจะใช้ script

<script type="text/javascript" src="http://www.wolfram.com
/cdf-player/plugin/v2.1/cdfplugin.js"></script>
<script type="text/javascript">
	var cdf = new cdfplugin();
	cdf.embed('/path/to/ไฟล์.cdf', width, height);
</script>

ที่มา

http://wolfram.com/cdf/adopting-cdf/deploying-cdf/web-delivery.html

 

####

เสียดายที่ CDF ใช้ไม่ได้กับ browser ตัวใหม่ๆแล้ว ทาง Wolfram เองก็แนะนำให้ใช้กับ Wolframcloud แทน ดูเพิ่มเติมได้ที่นี่ครับ http://www.wolfram.com/cdf/adopting-cdf/deploying-cdf/web-delivery-cloud.html

Categories
IT Mathematica Uncategorized

MathLink แบบง่ายๆ

เดี๋ยวมาเล่าให้ฟังครับ เขียนหัวข้อไว้ก่อน 😛

MathLink

Categories
IT Mathematica

มาใช้ Mathematica ในการเรียนการสอนกัน ด้วยคำสั่ง Manipulate 1

ผมจะมาเล่าให้ฟังครับว่าเราจะใช้โปรแกรมที่ชื่อ Mathematica ในการเรียนการสอนได้อย่างไร ด้วยคำสั่ง Manipulate แต่ก่อนจะไปถึงตอนนั้น ผมขอแนะนำว่าคนที่ไม่รู้ว่า Mathematica คืออะไรแล้วมันใช้ทำอะไรได้บ้างลองดูจากวีดีโอแนะนำนี้ดูครับ  QuickTour ส่วนใครที่กำลังเริ่มต้นเรียนรู้ใช้งานสามารถอ่านคำแนะนำการใช้งานเบื้องต้นของท่านอ.พงศกร สายเพ็ชร์ ได้ครับหรือจะไปพูดคุยสอบถามการใช้งาน Mathematica ได้ที่ http://mpec.sc.mahidol.ac.th/forums/ ครับ

ตัวอย่างจากการใช้คำสั่ง Manipulate ในการสร้างapplicationแบบต่างๆมีตัวอย่างมากมายที่ http://demonstrations.wolfram.com/ ต้องติดตั้ง Wolfram CDF Player ก่อนนะครับถึงจะเปิดดูพวกตัวอย่างนั้นได้  เจ้าตัว CDF Player นี้เองที่เราจะเปิดอ่านไฟล์พวก .nb หรือ .cdf ได้แต่จะไม่สามารถแก้ไขได้ครับ ต้องใช้ Mathematica อย่างเดียวครับ ถ้าจะสร้างหรือแก้ไฟล์ .nb  ส่วนไฟล์ .cdf นั้น Mathematica เปิดดูได้อย่างเดียวครับจะไม่สามารถแก้ไขตัว application ได้ครับ

ต่อตอนสองครับ(http://www.sakngoi.com/?p=290)

Categories
IT

External Hard Disk กับ Linux

ผมซื้อ External Hard Disk(HD) ตัวหนึ่งยี่ห้อ Western Digital (WD) ขนาด 1 TB มาเอาไว้ใช้สำรองข้อมูลจากเครื่อง laptop ของผมซึ่งเป็นระบบปฏิบัติการ Linux Debian (Lenny) ผมมีปัญหาอยู่ว่าพอต่อ HD ผ่าน USB port กับ laptop แล้วมันมองเห็นข้อมูลใน HD แต่ไม่สามารถเขียนลงไปได้  ซึ่งผมก็ได้ลองใช้คำสั่ง

fdisk -l

เพื่อดูว่า HD ใช้ filesystem แบบไหนก็พบว่าเป็น NTFS  แล้วใช้ คำสั่ง

mount -t ntfs /dev/sdb1 /media/XHD -o rw

ก็พบว่ายังไม่สามารถทำให้สามารถเขียนข้อมูลลงในHDได้  ก็เลยหาข้อมูลดูพบว่า Debian ตัวที่ผมใช้อยู่นี้ยังไม่ได้ลงpackage ที่จะช่วยให้สามารถเขียนข้อมูลบนfilesystem ที่เป็น ntfs ไว้ ซึ่ง package ที่ว่านี้ก็คือ ntfs-3g ผมจึงได้ทำการติดตั้ง package นี้ผ่าน Synaptic package manager ซึ่งหลังจากติดตั้งเสร็จผมก็ทำการ mount HD ด้วยคำสั่ง

mount -t ntfs-3g  /dev/sdb1 /media/XHD -o rw

ซึ่งก็สำเร็จครับ สามารถเขียน ลบ อ่าน ข้อมูลบน HD ตัวนี้ได้ตามที่ตั้งใจ

 

Categories
Mathematica

การใส่สีพื้นหลังของกราฟจากคำสั่ง Plot

การใส่สีพื้นหลังของกราฟจากคำสั่ง Plot มีหลายวิธีครับ ด้านล่างนี้เป็นตัวอย่างครับ

Plot[Sin[x], {x, 0, 4 \[Pi]}, PlotStyle -> Red, Background -> Gray]

 

Graphics[Inset[
Plot[Sin[x], {x, 0, 4 \[Pi]}, PlotStyle -> White,
Frame -> {True, True, False, False},
FrameStyle -> Directive[White, 12]]], Background -> Gray]

 

bg = Graphics[
Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}},
VertexColors -> {Black, Black, White, White}],
AspectRatio -> Full];
Plot[Sin[x], {x, 0, 4 \[Pi]}, PlotStyle -> Red,
Prolog -> Inset[bg, Scaled[{0, 0}], Scaled[{0, 0}], Scaled[{1, 1}]]]