בכל זאת Cython

מזמן לא כתבתי על שיפור ביצועים עם פייתון; אבל מכיוון שכולם עכשיו מסכמים שנה, והשנה שלי עמדה בסימן מחשוב מדעי עם פייתון, אולי הגיע הזמן לעוד פרק בסדרה.

הפעם אני חוזר ל-Cython, פתרון הכלאיים לייצור תוכנות מקומפלות לשפת מכונה ע"י תחביר דמוי פייתון ככל האפשר. זה נכון שבעבר הערתי שהפתרון הזה פוגע בניידות הקוד וכו', אבל אי אפשר להתווכח עם כך שהוא משפר את הביצועים במשהו. מכיוון שבזמן האחרון יש לי קצת זמן פנוי בשעה שאני מחכה שנושא התזה שלי יאושר או יידחה החלטתי לחזור ולבחון את העניין.

הסיבה לבחירה ב-Cython (מעכשיו סייתון) דווקא, ולא באביה הרוחני Pyrex, היא שבגרסאות חדשות יותר של סייתון (שאותן משתמשי אובונטו יכולים להתקין מהמאגרים של Jaunty אפילו באינטרפיד ללא קושי) נכנס קוד טרי טרי מה-Google Summer of Code האחרון. קוד זה שמספק תמיכה ספציפית במערכים של numpy מקל מאוד על העבודה. בפעם האחרונה שבדקתי את הנושא, למעשה כל שימוש במערך כזה או גישה לאלמנט שלו נעשו ע"י קריאה למתודות המתאימות על אובייקט פייתון שמייצג את המערך, מה שחרב את הביצועים לגמרי (ובכל זאת הצלחתי להשיג קוד שרץ במהירות כפולה).

עכשיו יש אפשרות להתייחס למערך numpy כאל חוצץ (buffer) בזיכרון, וזה הופך את הגישה אליו למהירה בהרבה. למעשה, חישובים על תאים של מערך כזה מתורגמים פחות או יותר לקוד C שהייתי כותב עבור ערכים בודדים.

ניתן וכדאי לוודא זאת ע"י בדיקת קוד ה-C שמיוצר ע"י סייתון. כדאי - כי סייתון מכיל בתוכו כמה המרות מובלעות, שאותן הוא מחליט אם לבצע או לא בלי להודיע לך, בהתבסס על השימוש שלך במשתני פייתון רגילים. הנה לדוגמא קוד שיומר לקוד איטי:
import numpy as N
cimport numpy as N

ctypedef N.float64_t DTYPE_t

...

cdef DTYPE_t coeff = N.pi*L*k # L, k are defined as DTYPE_t

במקרה זה, אפילו ש-k,L הם משתני C פשוטים, הנוכחות של אובייקט פייתון pi מתוך numpy מספיקה כדי להפוך את כל ההכפלות לפעולות על אובייקטי פייתון, שהן איטיות להחריד; ואם זה לא מספיק, עכשיו כל קוד שנוגע ב-coeff יסבול מאותה בעיה.

הקוד הבא דווקא רץ מהר:
DEF     M_PI = 3.14159265358979323846
...
cdef DTYPE_t coeff = M_PI*L*k
כאן השתמשתי בפקודה DEF כדי להגדיר קבוע (מועתק מ-math.h). עכשיו שרשרת ההכפלות הזו יוצרת קוד C פשוט ללא זכר לפייתון, והכל מתחיל לטוס.

שיטת הפעולה שלי היתה לקחת קוד שממילא נסמך רבות על numpy ולממש אותו בסייתון, ע"י פתיחתו ללולאה (כפולה) שמיישמת כל חישוב על מערך דו-מימדי כסדרת חישובים על ערכים בודדים, כל אחד בתורו. מעניין לשאול למה זה בכלל עובד, שהרי כל הפונקציות של numpy מיושמות ב-C - היית מצפה שלא תזכה לרווח גדול כל כך. תשובה שמצאתי באינטרנט היא שכאשר עובדים על מערכים גדולים, ומכיוון שבפייתון כל תוצאת ביניים תחושב במלואה לפני שהחישוב הבא ייקרה, יוצא שהערכים שמשתתפים בחישוב תמיד נדחפים החוצה מהמטמון של המעבד דווקא כשצריך אותם. כשמיישמים את כל חישובי הביניים כחלק מאיטרציה אחת של הלולאה, אפשר לשמור על הקרבה-בזמן של הנתונים לחישוב שלהם. לדוגמא:
    for ii in range(iimax):
        for jj in range(jjmax):
            bdiv = b[ii,jj] / a
            cdiv = c[ii,jj] / a
            Q = -(bdiv**2)/8.
            U = (0.5*Q + sqrt(Q**2/4. - cdiv**3/27.))**(1./3.)
            y = -cdiv/3./U - U
            W = sqrt(2*y)
            ret[ii,jj] = (sqrt(2*(bdiv/W - y)) - W) /2.

הקוד הזה, שפותר משוואה מסויימת ממעלה רביעית, רצוף חישובי ביניים, וכולם נעשים בקרבת-זמן אחד לשני. בגרסת הפייתון, החלוקה השניה היתה דוחפת מהמטמון את החלוקה הראשונה, שדווקא תוצאותיה דרושות בשורה הבאה.

אגב, גם הקוד הזה רץ באותה מהירות כמו קוד הפייתון המקורי בהתחלה, משום ששכחתי להגדיר את משתני הביניים מראש ע"י שימוש ב-cdef, ולכן הם טופלו כמשתני פייתון רגילים. לאחר הגדרה כמו שצריך, הקוד המלא (שאינו מובא כאן) רץ במהירות כפולה בערך. אולי עוד כדאי להעיר שהסיבה לכך שהשיפור לא גדול אף יותר היא שרוב הזמן מוקדש לחישוב של פונקציות שורש, שמיושמות ע"י libc כקריאה ישירה לפונקציית שורש מתוך סט פקודות MMX של המעבד (אצלי לפחות, בדקתי בחבילת קוד המקור המתאימה מאובונטו), ואת זה קשה להביס.

מכיוון ששתי הדוגמאות שהבאתי כאן לקוחות מלולאות פנימיות חשובות בתכנה שלי, בסך הכל הצלחתי להוריד את זמן הביצוע של מקרה המבחן שלי מ-48 שניות ל-34 שניות (כולל כל הקוד שלא קשור לקטעים האלה), שזה ממש לא רע כשזוכרים שסבבי הרצה מלאים כוללים מאות או אלפי מקרים דומים למקרה זה.

עכשיו רק נשאר לי לקוות שנושא התזה שלי יאושר כדי שיהיה שימוש לכל השיפור הזה :)

אה, ודבר אחרון: מי שבאמת הולך להכנס לזה, יימצא את המדריך למשתמשי numpy שימושי ביותר.