Tag: Mathematica

Categories
Mathematica

ประมาณค่า Pi

มีคนถามมาว่าช่วยทำตัวอย่างให้ดูหน่อยในการใช้ OpenCL ใน Mathematica ซึ่งผมก็บอกไปแล้วว่ามันตัวอย่างเยอะเลยในhelp Mathematica เอง 🙂 แต่ก็ยังอยากให้ผมทำให้ดู ตัวอย่างที่ทำให้ดูก็เป็นการประมาณค่า Pi โดยอาศัยที่เรียกว่า Monte carlo โดยประมาณจากการสุ่มจุดลงในสี่เหลี่ยมจัตุรัสที่มีวงกลมรัศมีครึ่งหนึ่งของความยาวด้านสี่เหลี่ยมจัตุรัส แล้วดูว่ามีกี่จุดที่อยู่ในวงกลม ซึ่งค่า Pi จะสามรถประมาณได้จาก

Pi = 4 * จำนวนจุดในวงกลม / จำนวนจุดที่สุ่มมาทั้งหมด

อันนี้เป็น code ที่ผมเขียนโดยใช้ OpenCL ใน Mathematica เพื่อหาดูว่าจุดไหนบ้างที่อยู่ในวงกลมซึ่งมีจุดศูนย์กลาง (0.5,0.5) ในสี่เหลี่ยนจัตุรัส

<< OpenCLLink`

scr = "
  __kernel void myPie( __global double *X, __global double *Y, 
__global mint *inC, mint np) {
      int indx = get_global_id(0);
      double d; //distance
  	if (indx < np) {
  		d = sqrt((X[indx]-0.5)*(X[indx]-0.5)+(Y[indx]-0.5)*(Y[indx]-0.5));
          if(d <= 0.5)
inC[indx] = 1;   
  	}
  }
  ";

fn = OpenCLFunctionLoad[scr, 
  "myPie", {{_Real, _, "Input"}, {_Real, _, "Input"}, {_Integer, _, 
    "Output"}, _Integer}, 16 ];

np = 10^6;
x = OpenCLMemoryLoad[RandomReal[1, np]];
y = OpenCLMemoryLoad[RandomReal[1, np]];
inc = OpenCLMemoryLoad[ConstantArray[0, np]];
Print["my Pi : ", (OpenCLMemoryGet@(fn[x, y, inc, np] // First) // 
    Total)*4.0/np]

OpenCLMemoryUnload[x, y, inc]

ทดลองรันดูสัก 100 ครั้ง โดยที่แต่ล่ะครั้งสุ่มมา 10^7 จุด แล้วดูว่ามันให้ค่าประมาณแตกต่างเท่าไหร่เมื่อเทียบกับค่าจริง (จาก Mathematica) จากรูปด้านล่าง เส้นประสีดำคือที่ได้จากการประมาณ สีแดงคือเส้นค่า Pi จาก Mathematica

Categories
Mathematica Uncategorized

ว่าด้วยเรื่อง Paclets ของภาษา Wolfram

เพิ่งอัพเดท Mathematica ไปใช้ 12.1 แล้วเจอเรื่อง paclet ที่เปลี่ยนไป (ในทางที่ดีขึ้น) เลยเดี๋ยวจะมาเล่าให้ฟัง ขอแปะแหล่งข้อมูลไว้ก่อน

1. https://www.wolframcloud.com/obj/tgayley/Published/PacletDevelopment.nb

2. https://www.wolfram.com/broadcast/video.php?sx=paclet&v=2833

3. https://reference.wolfram.com/language/guide/Paclets.html

Categories
Mathematica

แบบจำลองการเพิ่มขึ้นของเชื้อมาลาเรียในคน

กำลังเตรียมdraft งานเก่าที่ทำมาหลายปีแล้วไม่ได้มีโอกาสส่งตีพิมพ์สักที เลยเอามาปัดฝุ่นใหม่โดยใช้ Mathematica รู้สึกได้เลยว่า code มันเขียนไม่กี่บรรทัดเอง แถมแชร์ก็ง่าย https://www.wolframcloud.com/obj/sompob/Published/NJW_Im4a.nb

Categories
Mathematica

เรียก compiled function จาก compiled function

โดยปกติ Mathematica หรือภาษา Wolfram จะมีปัญหากับการเรียก function ที่สร้างจาก Compile ใน function ที่จะสร้างจาก Compile อีกทีได้ (ส่วนใหญ่เป็นเรื่องที่เกี่ยวกับว่า Funcions ไหนที่มัน compile ได้หรือไม่ได้) แต่ใน version 12 นี้เราสามารถที่จะทำแบบนั้นได้ ถ้า function เหล่านั้นถูก compile ด้วยคำสั่ง FunctionCompile ซึ่งผมมองว่ามันสะดวกอย่างมาก และที่สำคัญเราสามารถที่จะ export สร้างเป็น library ไฟล์ได้ด้วย โดยคำสั่ง FunctionCompileExport สกุลไฟล์ที่สร้างได้ก็ได้แก่

ตัวอย่างการเรียกใช้ compiled function ใน compiled function อีกที ก็สมมุติให้ผมมี function สำหรับทำ bisection เพื่อหาค่าที่ทำให้ function ที่สนใจเป็น 0 หรือใกล้ 0 มากๆ 

BisectionMethod = FunctionCompile[
    Function[
        {Typed[f, {"Real64"} -> "Real64"], Typed[lim0, "Real64"], Typed[lim1, "Real64"]},
      Module[{Maxiter = 100, tol = 10.^-8, iter = 0, x0 = lim0, x1 = lim1, mid = 0.0, f0, f1, fmid},
      f0 = f[x0]; 
      f1 = f[x1];
            While[Abs[x1 - x0] >= tol && iter++ < Maxiter,
                mid = (x1 + x0)/2.;
                fmid = f[mid];
                If[Xor[fmid > 0, f0 > 0],
                    x1 = mid; f1 = fmid,
                    x0 = mid; f0 = fmid]];
               mid]]]

โดยfunction ที่ผมจะหา root คือ f(x) = sin(x)+exp(x)

f = FunctionCompile[Function[Typed[arg, "Real64"], Sin[arg] + Exp[arg]]]

จะเห็นได้ว่าทั้ง BisectionMethod และ f ถูกสร้างจาก FunctionCompile แต่ f จะถูกเรียกเข้าไปใช้ใน BisectionMethod อีกที

BisectionMethod[f, -1.0, 1.0]-FindRoot[Sin[x] + Exp[x], {x, -1}][[1, 2]]
(* -5.18667*10^-11 *)
Categories
Mathematica

การแพร่กระจายของโรคไข้เลือดออกในจีน

พอดีว่าต้องทำ slides ไปโชว์ว่างานที่รับอยู่ไปถึงไหนแล้ว มันมีส่วนหนึ่งที่ต้องทำให้เห็นว่าผลที่ได้จากโมเดลการระบาดของไข้เลือดออกในจีนกับข้อมูลจริงนั้นมันใกล้เคียงกัน พอดีผมไปเจองานอยู่งานหนึ่งที่น่าสนใจคือ https://bmcmedicine.biomedcentral.com/articles/10.1186/s12916-015-0336-1

Figure 5

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

โดยผมเขียน code ตามนี้ครับ

Categories
Mathematica

เรื่องการเรียกใช้ package ใน package อีกทีของภาษา Wolfram

พอดีกำลังรับงานพัฒนาโมเดลตัวหนึ่งด้วยภาษา Wolfram แล้วเจอปัญหาเรื่องการเรียกใช้ Package ในอีก Package หนึ่ง code ที่เขาส่งมาให้ช่วยดูให้มันมีการเรียกใช้ package ตั้งแต่คำสั่ง BeginPackage เช่น

BeginPackage[ "Package`", {"Package1`", "Package2`"}]
  ...

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

ผมก็เลยเปลี่ยนวิธีการเรียก packages ด้วย Needs ไปใช้วิธีเรียกทีล่ะตัวหลังคำสั่ง BeginPackage แทน เช่น

BeginPackage[ "Package`"]
  Needs[ "Package1`"]
  Needs[ "Package2`"]
  ...

ด้วยวิธีนี้ตัว คำสั่งจาก package1 หรือ package2 จะถูกเรียกใช้ได้ภายใน package แต่ด้านนอก package จะเรียกไม่ได้ พูดง่ายๆคือไม่มีผลอะไรกับตัวแปรนอก package

Categories
Uncategorized

เรียกใช้ MathNet ใน Mathematica

ที่ต้องระวังคือถ้า methods ใน C# มีการใช้ตัวแปรเป็น Matrix<> หรือ Vector<> เวลาจะส่งผ่านค่า Matrix หรือ Vector จาก Mathematica เข้าไปซึ่งมันเป็น List จะต้องทำการแปลงเป็น Matrix หรือ Vector โดยอาศัย constructors จากตัว MathNet ก่อน

ตัวอย่างเช่น มี Method Neighbors ที่ต้องการ Matrix<double> ตามนี้

เวลาเรียกจาก Mathematica ก็ทำตามนี้ครับ

สังเกตตัวแปร mat ถูกสร้างขึ้นมาจาก คำสั่งนี้ครับ

NETNew["MathNet.Numerics.LinearAlgebra.Double.DenseMatrix", 10, 10, 
RandomReal[1, 100]]

Categories
Mathematica

อ่านตัวเลขจากภาพ

มีคนถามคำถามนี้ใน Facebook ของสมาคมโปรแกรมเมอร์ไทย

คำถามนี้สามารถทำได้ง่ายๆด้วยภาษา Wolfram ครับ โดยใช้คำสั่ง TextRecognize ตามนี้ครับ 

https://www.wolframcloud.com/obj/sompob/Published/demo_3Oct19.nb
Categories
Mathematica

แสดงผลลัพธ์จาก 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 ตามวิธีที่เขาแนะนำไว้ครับ

Categories
Uncategorized

เรื่อง plane geometry ใน Mathematica 12

Wolfram Mathematica 12 ที่เพิ่งจะเปิดตัวไปมีความสามารถเพิ่มขึ้นมาเยอะพอสมควรเลย และหนึ่งจากหลายความสามารถที่เพิ่มเข้ามานี้ที่ผมชอบมากก็คือความสามารถในเรื่องเรขาคณิตที่เรียกว่า Synthetic Geometry

เจ้า Synthetic Geometry ก็เป็นสาขาที่ศึกษาเรขาคณิตที่หาข้อสรุปจากการใช้เรื่องสัจพจน์หรือข้อความคาดการณ์ต่างๆเป็นเครื่องมือหลักในการแก้ปัญหา ลองดูความหมายเพิ่มเติมที่ https://en.wikipedia.org/wiki/Synthetic_geometry  นะครับ

ใน Mathematica 12 นี้เราสามารถที่จะสร้างไดอะแกรมของปัญหาได้จากคำสั่ง GeometricScene เช่น

ถ้ามีสามเหลี่ยมที่มีมุมทั้งสามมุมอยู่ที่ จุด (0,0),(1,0),(0,1) หรือสามเหลี่ยมมุมฉากที่มีฐานและความสูงคือ 1 หน่วย เราก็สามารถสร้างได้แบบนี้ครับ

จากคำสั่งด้านบน GeometricScene เราเพียงแต่บอกว่าเรามีจุดอะไรบ้าง ซึ่งในที่นี้ก็คือจุด a,b,c ที่อยู่ที่ {0,0},{1,0},{0,1} ตามลำดับ และก็บอกว่าให้สร้าง สามเหลี่ยมโดยใช้จุดที่กำหนดไว้

จากนั้นผมทดลองหาพื้นที่ของสามเหลี่ยมนี้ ซึ่งก็ผลลัพธ์ที่ได้ก็คือ 1/2

มาลองทำอะไรที่มันดูซับซ้อนขึ้นมาอีกหน่อยครับ  สมมุติว่าผมต้องการหาพื้นที่ของสามเหลี่ยมที่มีเหลี่ยมทั้งสามที่จุด k,l,n

และมีวงกลมที่ผ่านจุดทั้งสามนี้

โดยผมกำหนดให้ว่าระยะระหว่างจุด k และ n  คือ 3 หน่วย และมีจุด m อีก หนึ่งจุดที่อยู่ด้านนอกวงกลม โดยกำหนดให้ มุมที่ทำระหว่าง จุด l,m,n เท่ากับ 120 องศา และแขนของมุมทั้งสอง (เส้นตรงจากจุด m ไป n กับ เส้นตรงจาก จุด m ไป l) นี้ต้องเป็นเส้นสัมผัสวงกลม และก็ให้ เส้นตรง จากจุด m ไป l กับเส้นตรง จุด k n นี้ขนานกัน  555 ดูจะสับสนซับซ้อนอย่างค่อยเป็นค่อยไปนะครับ แต่ถ้าเห็น code และจะเข้าใจดีขึ้นครับ

สังเกตุนะครับว่าเราไม่ได้ระบุพิกัดของจุดใดๆเลย หลังจากที่อธิบายไดอะแกรมของเราเรียบร้อยแล้วผ่านคำสั่ง GeometricScene แล้วเราก็ให้ Mathematica สร้างภาพออกมาโดยใช้คำสั่ง RandomInstance ซึ่งมันก็จะสุ่มจุดขึ้นมาและวาดตามเงื่อนไขที่เรากำหนดไป ซึ่งถ้าเรารัน คำสั่ง RandomInstance และไม่กำหนด RandomSeeding ภาพวาดที่ออกมามันจะไม่ซ้ำกันเลยครับ ดังนั้นเราสามารถกำหนดหรือ fix ภาพ ได้จากการกำหนดค่า RandomSeeding ครับ

จากภาพที่วาดออกมาเราสามารถหาพื้นที่ของสามเหลี่ยน kln ได้ตามนี้ครับ

จากตัวอย่างข้างบนจะเห็นว่าผมใช้คำสั่ง GeometricAssertion เป็นตัวบอกว่า object ที่วาดขึ้นนั้น เช่น เส้นตรง วงกลม หรือจุดต่าง มันมีปฏิสัมพันธ์กันอย่างไร อย่างเช่น ผมต้องการ เส้นตรง kn กับเส้นตรง lm ขนานกัน ผมก็เพียงพิมพ์

GeometricAssertion[{Line[{k,n}],Line[{l,m}]},”Parallel”]

แต่ความเจ๋งของ Mathematica 12 ในเรื่อง Synthetic Geometry มันอยู่ตรง คำสั่ง FindGeometricConjecture นี้ครับ เจ้าคำสั่ง FindGemetricConjecture เป็นคำสั่งที่จะช่วยหาข้อสรุปต่างๆ จากตัวภาพที่วาดจากคำสั่ง GemetricScene ครับ เพื่อให้เห็นความสามารถของคำสั่งนี้ ผมจะลองให้มันหาดูว่ามุมภายในครึ่งวงกลมจะเป็นมุมฉากหรือเปล่าตามทฤษฎีของเธลิส เริ่มจากกำหนดจุด {a,b,c} อยู่บนวงกลมที่มีจุดศูนย์กลางที่จุด o และจุด o เป็นจุดกึ่งกลางระหว่างจุด a กับ c พอวาดเสร็จก็ให้หาข้อสรุปดู ตาม code ด้านล่างนี้

จะเห็นได้ว่าจาก คำสั่ง FindGeometricConjectures มันบอกว่าเส้นตรง ab กับ เส้นตรง cb มันตั้งฉากกัน ซึ่งมันก็โชว์อีกว่ามุม b ที่เกิดจากจุด abc มันทำมุม 90 องศา

ลองดูอีกสักตัวอย่าง เพื่อดูว่า จุด X,Y,Z ที่เกิดจากการตัดกันตามภาพล่างนี้อยู่บนระนาบเดียวกันหรือเปล่า

 ซึ่งจากผลที่ได้ X,Y,Z อยู่บนเส้นตรงเดียวกันครับ

ใครที่สนใจอยากดูตัวอย่างเพิ่มเติมลองไปดูได้ที่ https://www.wolfram.com/language/12/plane-geometry/