Category: Mathematica

เรื่องทั่วไปเกี่ยวกับการใช้งานโปรแกรม Mathematica

Categories
Mathematica

สมการชูสามนิ้ว

เห็นน้องๆสร้างสมการคณิตศาสตร์กันเป็นรูปมือกันผมเลยลองเล่นบ้าง เป็นรูปมือชู 3 นิ้ว (พระพุทธ พระธรรม พระสงฆ์ หรือเปล่า ?)

สมการที่ผมหามาได้มีหน้าตาตามนี้ครับ *มีคนถามมาว่าหามายังไง เดี๋ยว​มาอธิบายให้ครับหลักๆก็อาศัย Fourier series และภาพต้นฉบับก็เอามาจากเวบ https://www.123rf.com/stock-photo/three_finger.html?sti=lei4tdvr1lskwqekk0| ครับ

h3Function[
  t_] := {((30359/31 - 1/17 Sin[23/15 - 16 t] - 
        4/43 Sin[55/36 - 12 t] - 6/23 Sin[29/19 - 9 t] - 
        13/36 Sin[36/23 - 8 t] - 20/31 Sin[39/25 - 6 t] - 
        33/43 Sin[31/20 - 5 t] - 3091/309 Sin[36/23 - 2 t] - 
        92/29 Sin[36/23 - t] + 27/40 Sin[69/44 + 3 t] + 
        9/44 Sin[29/18 + 4 t] + 8/19 Sin[68/43 + 7 t] + 
        2/35 Sin[164/35 + 10 t] + 4/21 Sin[59/37 + 11 t] + 
        1/78 Sin[151/101 + 13 t] + 4/47 Sin[70/43 + 14 t] + 
        2/29 Sin[55/34 + 15 t] + 2/19 Sin[53/33 + 17 t] + 
        2/27 Sin[18/11 + 18 t] + 1/25 Sin[43/27 + 19 t] + 
        1/31 Sin[41/25 + 20 t] + 2/23 Sin[73/45 + 21 t] + 
        1/52 Sin[45/26 + 22 t] + 1/24 Sin[115/72 + 23 t] + 
        1/29 Sin[49/29 + 24 t] + 2/39 Sin[57/35 + 25 t] + 
        1/81 Sin[37/22 + 26 t] + 1/27 Sin[28/17 + 27 t] + 
        1/49 Sin[77/46 + 28 t] + 1/36 Sin[18/11 + 29 t] + 
        1/131 Sin[33/19 + 30 t] + 1/36 Sin[118/71 + 31 t] + 
        1/95 Sin[37/22 + 32 t] + 1/58 Sin[18/11 + 33 t] + 
        1/148 Sin[85/48 + 34 t] + 1/51 Sin[29/18 + 35 t] + 
        1/223 Sin[50/27 + 36 t] + 1/73 Sin[44/27 + 37 t] + 
        1/183 Sin[66/37 + 38 t] + 1/68 Sin[49/30 + 39 t] + 
        1/373 Sin[44/23 + 40 t] + 1/82 Sin[50/31 + 41 t] + 
        1/222 Sin[19/11 + 42 t]) UnitStep[
       15 \[Pi] - t] UnitStep[-11 \[Pi] + t] + (82913/46 - 
        6/31 Sin[12/25 - 34 t] - 2/17 Sin[45/31 - 28 t] - 
        2/7 Sin[21/16 - 26 t] - 15/22 Sin[26/35 - 17 t] - 
        97/39 Sin[69/44 - 9 t] - 116/31 Sin[19/39 - 8 t] - 
        281/108 Sin[57/41 - 7 t] - 138/31 Sin[20/53 - 6 t] - 
        864/23 Sin[31/21 - 2 t] + 13015/73 Sin[39/17 + t] + 
        2183/38 Sin[231/68 + 3 t] + 473/59 Sin[118/47 + 4 t] + 
        253/35 Sin[53/15 + 5 t] + 88/63 Sin[142/41 + 10 t] + 
        39/77 Sin[103/38 + 11 t] + 61/33 Sin[107/33 + 12 t] + 
        32/23 Sin[17/9 + 13 t] + 2/23 Sin[54/25 + 14 t] + 
        5/16 Sin[451/450 + 15 t] + 5/11 Sin[6/53 + 16 t] + 
        9/29 Sin[37/9 + 18 t] + 3/16 Sin[173/38 + 19 t] + 
        3/7 Sin[198/61 + 20 t] + 8/21 Sin[37/15 + 21 t] + 
        7/26 Sin[59/33 + 22 t] + 7/31 Sin[57/65 + 23 t] + 
        1/6 Sin[14/29 + 24 t] + 9/34 Sin[5/16 + 25 t] + 
        3/32 Sin[445/99 + 27 t] + 13/66 Sin[131/34 + 29 t] + 
        7/29 Sin[103/40 + 30 t] + 3/31 Sin[100/99 + 31 t] + 
        2/41 Sin[37/16 + 32 t] + 3/20 Sin[13/28 + 33 t] + 
        9/50 Sin[23/5 + 35 t] + 2/23 Sin[105/23 + 36 t] + 
        3/28 Sin[310/67 + 37 t] + 2/15 Sin[143/45 + 38 t] + 
        4/33 Sin[48/29 + 39 t] + 1/20 Sin[34/57 + 40 t] + 
        1/32 Sin[125/43 + 41 t] + 3/37 Sin[49/44 + 42 t]) UnitStep[
       11 \[Pi] - t] UnitStep[-7 \[Pi] + t] + (28702/23 - 
        4/39 Sin[33/23 - 42 t] - 53/37 Sin[51/76 - 40 t] - 
        4/21 Sin[26/43 - 38 t] - 77/37 Sin[7/25 - 33 t] - 
        5/41 Sin[77/58 - 32 t] - 106/77 Sin[20/19 - 29 t] - 
        71/19 Sin[23/24 - 27 t] - 97/42 Sin[3/28 - 26 t] - 
        127/17 Sin[11/26 - 20 t] - 42/19 Sin[10/71 - 19 t] - 
        142/31 Sin[49/45 - 17 t] - 465/29 Sin[23/20 - 14 t] - 
        418/37 Sin[1/14 - 13 t] - 51/25 Sin[37/29 - 12 t] - 
        365/53 Sin[34/27 - 9 t] - 6795/137 Sin[19/36 - 7 t] - 
        241/16 Sin[11/28 - 6 t] - 5584/29 Sin[19/39 - 2 t] + 
        21149/52 Sin[49/11 + t] + 1701/17 Sin[29/35 + 3 t] + 
        862/35 Sin[26/9 + 4 t] + 667/28 Sin[92/25 + 5 t] + 
        577/78 Sin[79/44 + 8 t] + 193/18 Sin[1/27 + 10 t] + 
        178/27 Sin[71/19 + 11 t] + 359/57 Sin[34/33 + 15 t] + 
        135/38 Sin[91/24 + 16 t] + 91/44 Sin[111/31 + 18 t] + 
        160/29 Sin[118/27 + 21 t] + 84/29 Sin[9/40 + 22 t] + 
        52/19 Sin[114/37 + 23 t] + 113/46 Sin[161/36 + 24 t] + 
        10/33 Sin[2977/744 + 25 t] + 27/16 Sin[37/11 + 28 t] + 
        39/20 Sin[77/31 + 30 t] + 7/9 Sin[51/13 + 31 t] + 
        47/34 Sin[145/31 + 34 t] + 14/29 Sin[86/39 + 35 t] + 
        29/28 Sin[532/145 + 36 t] + 37/28 Sin[80/43 + 37 t] + 
        3/16 Sin[10/13 + 39 t] + 5/26 Sin[35/17 + 41 t]) UnitStep[
       7 \[Pi] - t] UnitStep[-3 \[Pi] + t] + (40485/29 - 
        255/23 Sin[2/21 - 23 t] - 1023/40 Sin[25/16 - 16 t] - 
        148/11 Sin[47/31 - 15 t] - 1781/43 Sin[27/34 - 10 t] - 
        9508/15 Sin[48/37 - t] + 3391/22 Sin[17/25 + 2 t] + 
        2003/35 Sin[22/19 + 3 t] + 3390/61 Sin[89/21 + 4 t] + 
        3313/25 Sin[82/27 + 5 t] + 671/8 Sin[70/41 + 6 t] + 
        574/13 Sin[58/35 + 7 t] + 1293/16 Sin[78/67 + 8 t] + 
        1462/37 Sin[17/21 + 9 t] + 727/32 Sin[96/23 + 11 t] + 
        1747/37 Sin[367/100 + 12 t] + 68/3 Sin[45/91 + 13 t] + 
        1063/31 Sin[41/35 + 14 t] + 54/11 Sin[55/16 + 17 t] + 
        165/8 Sin[224/225 + 18 t] + 135/38 Sin[1/18 + 19 t] + 
        583/63 Sin[113/24 + 20 t] + 205/23 Sin[220/73 + 21 t] + 
        461/33 Sin[7/18 + 22 t] + 53/5 Sin[99/23 + 24 t] + 
        102/25 Sin[82/47 + 25 t] + 63/11 Sin[31/16 + 26 t] + 
        224/27 Sin[117/25 + 27 t] + 121/16 Sin[58/33 + 28 t] + 
        49/9 Sin[21/17 + 29 t] + 107/50 Sin[3/7 + 30 t] + 
        38/9 Sin[139/34 + 31 t] + 55/23 Sin[46/19 + 32 t] + 
        136/39 Sin[44/29 + 33 t] + 13/7 Sin[14/39 + 34 t] + 
        11/3 Sin[26/15 + 35 t] + 285/89 Sin[118/29 + 36 t] + 
        28/19 Sin[13/17 + 37 t] + 53/35 Sin[89/34 + 38 t] + 
        18/19 Sin[114/47 + 39 t] + 39/44 Sin[9/13 + 40 t] + 
        16/29 Sin[185/52 + 41 t] + 69/26 Sin[56/33 + 42 t]) UnitStep[
       3 \[Pi] - t] UnitStep[\[Pi] + t]) UnitStep[Sqrt[
    Sign[Sin[t/
      2]]]], ((-(248369/91) - 1/102 Sin[127/85 - 28 t] - 
        1/31 Sin[59/38 - 18 t] - 1/110 Sin[71/47 - 12 t] - 
        11/28 Sin[57/37 - 5 t] - 3/5 Sin[69/44 - 4 t] - 
        903/139 Sin[36/23 - 2 t] - 773/119 Sin[47/30 - t] + 
        55/23 Sin[41/26 + 3 t] + 9/64 Sin[57/35 + 6 t] + 
        13/23 Sin[65/41 + 7 t] + 1/21 Sin[71/42 + 8 t] + 
        17/37 Sin[27/17 + 9 t] + 18/47 Sin[29/18 + 10 t] + 
        11/40 Sin[51/32 + 11 t] + 5/19 Sin[53/33 + 13 t] + 
        3/16 Sin[21/13 + 14 t] + 3/35 Sin[71/44 + 15 t] + 
        6/37 Sin[47/29 + 16 t] + 2/9 Sin[21/13 + 17 t] + 
        1/19 Sin[139/87 + 19 t] + 1/13 Sin[41/25 + 20 t] + 
        1/113 Sin[37/24 + 21 t] + 1/231 Sin[55/29 + 22 t] + 
        1/12 Sin[34/21 + 23 t] + 1/213 Sin[119/27 + 24 t] + 
        1/199 Sin[41/29 + 25 t] + 1/131 Sin[71/41 + 26 t] + 
        1/57 Sin[28/17 + 27 t] + 1/45 Sin[57/35 + 29 t] + 
        1/98 Sin[22/13 + 30 t] + 1/69 Sin[28/17 + 31 t] + 
        1/836 Sin[91/20 + 32 t] + 1/53 Sin[18/11 + 33 t] + 
        1/246 Sin[64/35 + 34 t] + 1/79 Sin[51/32 + 35 t] + 
        1/126 Sin[44/25 + 36 t] + 1/54 Sin[18/11 + 37 t] + 
        1/306 Sin[65/34 + 38 t] + 1/66 Sin[18/11 + 39 t] + 
        1/171 Sin[97/54 + 40 t] + 1/80 Sin[77/47 + 41 t] + 
        1/212 Sin[79/45 + 42 t]) UnitStep[
       15 \[Pi] - t] UnitStep[-11 \[Pi] + t] + (-(27062/15) - 
        1/21 Sin[49/39 - 42 t] - 1/12 Sin[29/68 - 36 t] - 
        3/22 Sin[12/19 - 33 t] - 2/33 Sin[16/21 - 30 t] - 
        1/15 Sin[2/19 - 27 t] - 3/35 Sin[1/325 - 24 t] - 
        9/20 Sin[58/37 - 19 t] - 30/89 Sin[59/42 - 16 t] - 
        13/24 Sin[22/19 - 13 t] - 16/21 Sin[32/25 - 10 t] + 
        7171/19 Sin[48/29 + t] + 118/5 Sin[39/14 + 2 t] + 
        8283/202 Sin[49/32 + 3 t] + 276/31 Sin[149/36 + 4 t] + 
        245/19 Sin[71/32 + 5 t] + 116/55 Sin[39/22 + 6 t] + 
        8/13 Sin[68/41 + 7 t] + 51/35 Sin[61/19 + 8 t] + 
        13/14 Sin[29/37 + 9 t] + 35/61 Sin[83/25 + 11 t] + 
        19/22 Sin[37/35 + 12 t] + 14/23 Sin[41/18 + 14 t] + 
        13/30 Sin[10/19 + 15 t] + 7/34 Sin[20/7 + 17 t] + 
        13/38 Sin[4/47 + 18 t] + 13/42 Sin[115/41 + 20 t] + 
        13/40 Sin[25/63 + 21 t] + 5/31 Sin[108/25 + 22 t] + 
        17/67 Sin[103/69 + 23 t] + 2/37 Sin[142/33 + 25 t] + 
        2/23 Sin[72/25 + 26 t] + 5/32 Sin[203/53 + 28 t] + 
        3/22 Sin[29/14 + 29 t] + 3/25 Sin[116/37 + 31 t] + 
        6/41 Sin[20/19 + 32 t] + 1/12 Sin[105/31 + 34 t] + 
        3/29 Sin[28/23 + 35 t] + 3/34 Sin[239/60 + 37 t] + 
        4/33 Sin[67/36 + 38 t] + 2/27 Sin[4/31 + 39 t] + 
        2/41 Sin[35/12 + 40 t] + 2/19 Sin[11/17 + 41 t]) UnitStep[
       11 \[Pi] - t] UnitStep[-7 \[Pi] + t] + (-(256299/275) - 
        7/16 Sin[4/35 - 40 t] - 31/29 Sin[7/22 - 39 t] - 
        16/31 Sin[7/22 - 34 t] - 28/23 Sin[29/32 - 33 t] - 
        34/19 Sin[27/26 - 31 t] - 47/25 Sin[21/37 - 27 t] - 
        183/64 Sin[13/15 - 26 t] - 697/209 Sin[13/16 - 20 t] - 
        103/27 Sin[22/51 - 19 t] - 41/42 Sin[4/23 - 15 t] - 
        271/53 Sin[3/43 - 14 t] - 207/38 Sin[59/60 - 13 t] - 
        286/23 Sin[22/27 - 11 t] - 517/54 Sin[33/65 - 8 t] - 
        171/22 Sin[13/30 - 7 t] - 7885/34 Sin[22/59 - t] + 
        8146/43 Sin[28/27 + 2 t] + 10347/22 Sin[28/13 + 3 t] + 
        7433/38 Sin[11/38 + 4 t] + 3353/72 Sin[37/20 + 5 t] + 
        331/46 Sin[170/37 + 6 t] + 1573/131 Sin[37/27 + 9 t] + 
        1999/41 Sin[63/40 + 10 t] + 173/31 Sin[35/32 + 12 t] + 
        656/101 Sin[49/29 + 16 t] + 437/32 Sin[23/22 + 17 t] + 
        114/23 Sin[4 + 18 t] + 71/20 Sin[3/22 + 21 t] + 
        84/37 Sin[69/32 + 22 t] + 143/43 Sin[61/35 + 23 t] + 
        231/68 Sin[5/22 + 24 t] + 92/55 Sin[112/33 + 25 t] + 
        3/2 Sin[1/7 + 28 t] + 41/18 Sin[55/27 + 29 t] + 
        59/46 Sin[153/76 + 30 t] + 16/35 Sin[1/10 + 32 t] + 
        26/77 Sin[8/33 + 35 t] + 50/41 Sin[45/23 + 36 t] + 
        13/25 Sin[8/3 + 37 t] + 40/37 Sin[346/77 + 38 t] + 
        7/34 Sin[37/30 + 41 t] + 1/3 Sin[46/15 + 42 t]) UnitStep[
       7 \[Pi] - t] UnitStep[-3 \[Pi] + t] + (-(43691/25) - 
        108/35 Sin[30/29 - 40 t] - 118/43 Sin[29/22 - 39 t] - 
        89/25 Sin[6/17 - 35 t] - 278/101 Sin[11/37 - 32 t] - 
        89/24 Sin[11/36 - 30 t] - 338/49 Sin[93/94 - 26 t] - 
        155/14 Sin[1/19 - 23 t] - 520/27 Sin[25/34 - 18 t] - 
        3729/233 Sin[15/13 - 17 t] - 1935/22 Sin[59/38 - 13 t] - 
        1624/33 Sin[31/23 - 9 t] - 299/12 Sin[26/43 - 6 t] - 
        8789/23 Sin[36/37 - 4 t] - 794/7 Sin[15/41 - 3 t] - 
        19163/23 Sin[43/36 - 2 t] + 24528/29 Sin[97/30 + t] + 
        2255/24 Sin[144/37 + 5 t] + 29915/277 Sin[233/58 + 7 t] + 
        1494/35 Sin[71/95 + 8 t] + 991/22 Sin[65/14 + 10 t] + 
        3579/61 Sin[31/27 + 11 t] + 537/32 Sin[106/57 + 12 t] + 
        463/42 Sin[18/55 + 14 t] + 1027/29 Sin[31/35 + 15 t] + 
        362/17 Sin[46/13 + 16 t] + 88/27 Sin[131/28 + 19 t] + 
        91/37 Sin[173/39 + 20 t] + 374/29 Sin[48/25 + 21 t] + 
        253/24 Sin[155/33 + 22 t] + 227/33 Sin[1/608 + 24 t] + 
        75/37 Sin[49/30 + 25 t] + 242/23 Sin[249/56 + 27 t] + 
        133/18 Sin[7/5 + 28 t] + 54/25 Sin[1/358 + 29 t] + 
        49/31 Sin[109/28 + 31 t] + 50/33 Sin[1/22 + 33 t] + 
        7/5 Sin[40/43 + 34 t] + 214/55 Sin[89/24 + 36 t] + 
        241/43 Sin[12/25 + 37 t] + 20/13 Sin[280/279 + 38 t] + 
        105/79 Sin[53/19 + 41 t] + 87/26 Sin[65/37 + 42 t]) UnitStep[
       3 \[Pi] - t] UnitStep[\[Pi] + t]) UnitStep[Sqrt[
    Sign[Sin[t/2]]]]}

ทดลอง plot

Categories
Mathematica

Protected: แยกไข่พยาธิด้วย Neural Network

This content is password protected. To view it please enter your password below:

Categories
Mathematica

OpenCL ในภาษา Wolfram

มีคนถามมาว่าช่วยทำตัวอย่างให้ดูหน่อยในการใช้ 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

ส่วนอีกอันก็เป็นการหมุนภาพ โดยที่มีจุดหมุนอยู่ที่กลางภาพครับ 

scr = "
  __kernel void test1( __global float * mat1, __global float * mat2, 
float theta, mint width, mint height) {
      int xIndex = get_global_id(0);
  	int yIndex = get_global_id(1);
  	int indx = xIndex + yIndex*width;
  
  	 if (xIndex >= width || yIndex >= height)
  		return ;
  
  	//image centre
  	float x0 = width/2;
  	float y0 = height/2;
  
  	//relative position
  	int xprime = xIndex-x0;
  	int yprime = yIndex-y0;
  
  	float sinTheta = sin(theta);
  	float cosTheta = cos(theta);
  
  	int nx =(xprime*cosTheta - yprime*sinTheta + x0);
  	int ny =(xprime*sinTheta + yprime*cosTheta + y0);
  	mat2[indx] = mat1[nx + width*ny];
  }
  
  ";

Clear[test];
test = OpenCLFunctionLoad[scr, 
  "test1", {{"Float", _, "Input"}, {"Float", _, "Output"}, 
   "Float", _Integer, _Integer}, {16, 16}, 
  "ShellOutputFunction" -> Print];
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

เว้นวรรค

มีคนส่งมาถามว่ามีข้อความเป็นภาษาอังกฤษแต่อยากให้เว้นวรรคหรือใส่space ให้หน่อยถ้ามันขึ้นต้นด้วยตัวพิมพ์ใหญ่ เช่น “AasdfasfdAsdfasd” อยากให้มันเป็น “Aasdfasfd Asdfasd” ผมเขียน script ง่ายให้ดังนี้ครับ

AddSpace[str_] := (StringReplace[str, x_?UpperCaseQ :> ” ” <> x ] // StringTrim)

ซึ่งเป็นภาษา Wolfram ครับ

อันนี้เอามาลองกับชื่อจังหวัดที่พิมพ์ติดกันมา

ที่อยากจะเน้นคือเรื่อง pattern ในภาษา Wolfram นี่มีประโยชน์มาก ทำให้เราทำงานได้เร็วขึ้นครับ

Categories
Mathematica

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

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

Categories
Mathematica

มันจะกระจายยังไงถ้าไม่ทำอะไรเลย

ทำเล่นๆขำ ว่าถ้าไวรัสที่ว่ามันกระจายไปตามจังหวัดต่างๆที่ว่า จำนวนคนติดเชื้อมันพุ่งเยอะขนาดไหน แบบว่าถ้าคนมันคิดว่าเป็นหวัดธรรมดา

โห … ตัวเลขพุ่งขึ้นเร็วมากกก ถ้าเห็นโมเดลทำนายเรื่องอะไรก็ตามอยากให้คิดกันเสมอว่า

“All models are wrong, but some are useful”

ฉะนั้นอย่าเพิ่งไปเชื่ออะไรง่ายๆครับ ทุกโมเดลมันมีสมมุติฐานในการสร้างของมัน ระบบในธรรมชาติมันซับซ้อนมากครับ ไม่มีใครที่จะโมเดลมันได้หมดหรอกครับ ทุกโมเดลมันจะตัดความยุ่งยาก ซับซ้อนออกเพื่อให้โมเดลได้ง่ายครับ

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