Friday, December 5, 2014

Use a vector to print Pascal's triangle in SAS

Yesterday Rick Wicklin showed a cool SAS/IML function to use a matrix and print a Pascal’s triangle. I come up with an alternative solution by using a vector in SAS/IML.

Method

Two functions are used, including a main function PascalRule and a helper function _PascalRule. The helper function recycles the vector every time and fills the updated values; the main function increases the length of the vector from 1 to n.

Pro

Get the nth row directly, for example, return the 10th row by PascalRule(10); no need to use a matrix or matrix related operators; use less memory to fit a possibly bigger n.

Con

More lines of codes; slowlier to print the triangle, since there is no data structure such as matrix to remember the transient numbers.
proc iml;
/* The main function */
start PascalRule(n);
if n <= 1 then
return({1})
;
answer = {1, 1};
do i = 1 to n - 2 ;
answer = _PascalRule(answer);
end;
return(answer);
finish;
/* The helper function */
start _PascalRule(vector);
previous = 1;
do i = 2 to nrow(vector);
current = vector[i];
vector[i] = previous + current;
previous = current;
end;
vector = vector // {1};
return(vector);
finish;
/* Print the pascal's triangle */
do i = 1 to 10;
x = PascalRule(i);
x = x`;
print x;
end;
quit;
Theoretically, Rick’s solution has a time complexity of O(N^2) and a space complexity of O(N^2), while my solution has a time complexity of O(N^3) (unfortunately have to use three times of do loop in IML) and a space complexity of O(N). Actually it's a trade-off between speed and memory cost.
​

Friday, November 28, 2014

Minimize complexity by Spark

There is always a trade-off between time complexity and space complexity for computer programs. Deceasing the time cost will increase space cost, and vice versa, The ideal solution to parallelize the program to multiple cores if there is a multiple-core computer, or even scale it out to multiple machines across a cluster, which would eventually reduce both time complexity and space complexity.
Spark is currently the hottest platform for cluster computing on top of Hadoop, and its Python interface provides map, reduce and many other methods, which allow a mapRecdue job in a straightforward way, and therefore easily migrate an algorithm from a single machine to a cluster of many machines.
  • Minimize space complexity
There is a question to look for the only single number from a mostly paired-number array.
Single Number
  Given an array of integers, every element appears twice except for one.
Find that single one.
Note:
Your algorithm should have a linear runtime complexity.
Could you implement it without using extra memory?
The optimal space complexity for this question is O(1) by using the bit manipulator xor. For a cluster, since Spark aggregates memory acrosss machines, the space complexity may become O(1/k), where k is the number of the machines in the cluster.
# Space.py
import pyspark
from random import shuffle
from operator import xor
sc = pyspark.Spark.Context()

# Create the test case and the target is 99
testCase = range(0, 99) + range(0, 100)
shuffle(testCase)
# Run the testing with Spark
result = sc.parallelize(testCase).reduce(xor)
# Show the result
print result
sc.stop()
  • Minimize time complexity
There is a question to implement the function (or a method) that returns the square root of an integer.
Sqrt(x)
Implement int sqrt(int x).
Compute and return the square root of x.
The optimal solution could achieve the time complexity of O(lgN) by using binary search. If we pass the sqrt function to Spark, then the time complexity will decreased to O(lgN/k), where k is the number of the machines in the cluster.
# Time.py
import pyspark
sc = pyspark.Spark.Context()
# Implement binary search for square root function
def sqrt(x):
if x < 0 or not isinstance(x, int):
raise ValueError
hi, lo = x/2 + 1, 0
while hi >= lo:
mid = (hi + lo) / 2
if mid * mid > x:
hi = mid - 1
else:
lo = mid + 1
return int(hi)

# Test the square root algorithm
testCase = sc.parallelize(xrange(1, 100))
result = testCase.map(lambda x: sqrt(x))
# Show the result
for x in result.collect():
print x
sc.stop()
  • Find the worst rating by accounts
There is a question to find the worst one among a few rating letters for each of the account numbers.
Want to find the worst rating for each account number.
sample.txt is below
Account_number    Rating
1 A
1 B
2 A
2 B
2 C
3 A
3 C
the desired result should be like
1            B
2 C
3 C
The question is essentially one of the grouping questons. Spark’s pair RDD, which reflects the key-value relationship for groups, supplies a one-line solution for it.
import pyspark
sc = pyspark.SparkContext()

# Squeeze the letters by keys
rdd = sc.textFile('sample.txt')
result = rdd.map(lambda x: x.split()).filter(x: x[0].isdigit()).reduceByKey(max)
# Show the result
for x in result.collect():
print x
sc.stop()
In a conclusion, Spark significantly changes the way we think about data analysis.
​

Monday, October 20, 2014

Automated testing by pytest

The most hard part in testing is to write test cases, which is time-consuming and error-prone. Fortunately, besides Python built-in modules such as doctest, unittest, there are quite a few third-party packages that could help with automated testing. My favorite one is pytest, which enjoys proven record and syntax sugar.

Step 1: test-driven development

For example, there is a coding challenge on Leetcode:
Find Minimum in Rotated Sorted Array
Suppose a sorted array is rotated at some pivot unknown to you beforehand.
(i.e., 0 1 2 4 5 6 7 might become 4 5 6 7 0 1 2).
Find the minimum element.
You may assume no duplicate exists in the array.
The straightforward way to find a minimal element in an array(or list in Python) is sequential searching, which goes through every element and has a time complexity of O(N). If the array is sorted, then the minimal one is the first element that only costs O(1).
However, this question provides a rotated sorted array, which suggests a binary search and reduces the complexity from O(N) to O(logN).
As usual, write the test cases first. The great thing for pytest is that it significantly simplies the effort to code the test cases: in this example, I only use 3 lines to generate 101 test cases to cover all conditions from 0 to 99 and also include an null test.
Next step is to code the function. It is easy to transplant the iterative approach of binary search to this question. If the pointer is between a sorted segment, then return the most left element as minimal. Otherwise, adjust the right boundary and the left boundary.
# test1.py
import pytest

# Prepare 101 test cases
array = list(range(100))
_testdata = [[array[i: ] + array[ :i], 0] for i in range(100)]
_testdata += [pytest.mark.empty(([], None))]

# Code the initial binary search function
def findMinPrev(num):
lo, hi = 0, len(num) - 1
while lo <= hi:
if num[lo] <= num[hi]:
return num[lo]
mid = (lo + hi) / 2
if num[mid] < num[hi]:
hi = mid - 1
else:
lo = mid + 1

@pytest.mark.parametrize('input, expected', _testdata)
def test_findMinPrev(input, expected):
assert findMinPrev(input) == expected
After running the py.test -v test1.py command, part of the results shows below. 65 tests passed and 36 failed; the failed cases return the much bigger values that suggests out of boundary during loops, and the selection of the boudaries may be too aggresive.
test1.py:20: AssertionError
_________________________ test_findMinPrev[input98-0] _________________________

input = [98, 99, 0, 1, 2, 3, ...], expected = 0

@pytest.mark.parametrize('input, expected', _testdata)
def test_findMinPrev(input, expected):
> assert findMinPrev(input) == expected
E assert 98 == 0
E + where 98 = findMinPrev([98, 99, 0, 1, 2, 3, ...])

test1.py:20: AssertionError
==================== 36 failed, 65 passed in 0.72 seconds =====================
Now I adjust the right boundary slightly and finally come up with a solution that passes all the tests.
def findMin(num):
lo, hi = 0, len(num) - 1
while lo <= hi:
if num[lo] <= num[hi]:
return num[lo]
mid = (lo + hi) / 2
if num[mid] < num[hi]:
hi = mid
else:
lo = mid + 1

Step 2: performance profiling

Besides the right solution, I am also interested in if the binary search method has indeed improved the performance. This step I choose line_profiler given its line-by-line ability of profiling. I take the most basic one (the sequential search) as benchmark, and also include the method that applies the min function since a few functions similar to it in Pyhton implement vectorizaiton to speed up. The test case is a rotated sorted array with 10 million elements.
# test2.py
from line_profiler import LineProfiler
from sys import maxint

@profile
def findMinRaw(num):
"""Sequential searching"""
if not num:
return
min_val = maxint
for x in num:
if x < min_val:
min_val = x
return min_val

@profile
def findMinLst(num):
"""Searching by list comprehension"""
if not num:
return
return min(num)

@profile
def findMin(num):
""""Binary search"""
lo, hi = 0, len(num) - 1
while lo <= hi:
if num[lo] <= num[hi]:
return num[lo]
mid = (lo + hi) / 2
if num[mid] < num[hi]:
hi = mid
else:
lo = mid + 1

# Prepare a rotated array
array = list(range(10000000))
_testdata = array[56780: ] + array[ :56780]
# Test the three functions
findMinRaw(_testdata)
findMinLst(_testdata)
findMin(_testdata)
After running kernprof -l -v test2.py, I have the output as below. The sequential search has hit the loops 10000001 times and costs almost 14 seconds. The min function encapsulate all details inside and uses 0.5 seconds which is 28 times faster. On the contrary, the binary search method only takes 20 loops to find the minimal value and spends just 0.0001 seconds. As a result, while dealing with large number, an improved algorithm can really save time.
Total time: 13.8512 s
File: test2.py
Function: findMinRaw at line 4

Line # Hits Time Per Hit % Time Line Contents
==============================================================
4 @profile
5 def findMinRaw(num):
6 1 13 13.0 0.0 if not num:
7 return
8 1 3 3.0 0.0 min_val = maxint
9 10000001 16031900 1.6 47.5 for x in num:
10 10000000 17707821 1.8 52.5 if x < min_val:
11 2 5 2.5 0.0 min_val = x
12 1 3 3.0 0.0 return min_val

Total time: 0.510298 s
File: test2.py
Function: findMinLst at line 15

Line # Hits Time Per Hit % Time Line Contents
==============================================================
15 @profile
16 def findMinLst(num):
17 1 4 4.0 0.0 if not num:
18 return
19 1 1243016 1243016.0 100.0 return min(num)

Total time: 0.000101812 s
File: test2.py
Function: findMin at line 22

Line # Hits Time Per Hit % Time Line Contents
==============================================================
22 @profile
23 def findMin(num):
24 1 15 15.0 6.0 lo, hi = 0, len(num) - 1
25 20 40 2.0 16.1 while lo <= hi:
26 20 48 2.4 19.4 if num[lo] <= num[hi]:
27 1 2 2.0 0.8 return num[lo]
28 19 54 2.8 21.8 mid = (lo + hi) / 2
29 19 50 2.6 20.2 if num[mid] < num[hi]:
30 5 10 2.0 4.0 hi = mid
31 else:
32 14 29 2.1 11.7 lo = mid + 1
​

Wednesday, September 24, 2014

One example of test-driven development in Python

Function or method is the most basic unit in Python programming. Test-driven development is a key for a developer to assure the code quality of those units. In his book, Harry Percival illustrated a few great examples about how to use TDD with Python and Django. It seems that for web development, TDD including unit testing and integration testing is the cornerstone for every success. For data analysis, coding mostly relies on built-in packages instead large framework like Django, which makes TDD easier. In my opnion, TDD in data analysis could have three steps.
  • Step 1: requirement analysis
    Before writing any code for data analysis, the programmer should seriously ask the customer or himself about the requirements.
    • What are the input parameter?
    • what if the input data doesn’t fit the assumptions?
    • What is the purpose of this funtction or method? what are the desired outputs?
  • For example, there is a recent coding challenge called Maximum Product Subarray.
      > Find the contiguous subarray within an array (containing at least one number) which has the largest product.
    For example, given the array [2,3,-2,4],
    the contiguous subarray [2,3] has the largest product = 6.
  • OK, understanding this question is quite straight-forward. Given a array(or a list in Python), you return the integer that is the maximum product from a continuous subarry out of the input array.
      def maxProduct(A):
    """ A function to find the maximum product value for a continuous subarray.
    :param A: an array or list
    :type A: list
    """

    if A is None or not isinstance(A, list):
    return None
    if len(A) == 1:
    return A[0]
    pass
    • A production version of the codes above should be more like:
      class FunctionInputError(Exception):
      pass

      def maxProduct(A):
      """ A function to find the maximum product value for a continuous subarray.
      :param A: an array or list
      :type A: list
      """

      if A is None or not isinstance(A, list):
      raise FunctionInputError('must give a list as input')
      if len(A) == 1:
      return A[0]
      pass
  • Step 2: write test cases
    Given not a single line of logic codes has been writen yet, I call the current step as black-box testing, which means that I want this funtion to fail every test cases. Python has a built-in module doctest, which allows embedding the test cases within the docstring. I write six test cases, run the psedu-function below and arbitrarily specify the result to be -1. As expected, it fails all the six test cases with horrible red warnings under the Python shell. That is a good thing: it proves that the testing works.
      import doctest
    def maxProduct(A):
    """ A function to find the maximum product value for a continuous subarray.
    :param A: an array or list
    :type A: list

    - testcase1
    >>> maxProduct([0, 2, 3,-2,4, 1, -1])
    48

    - testcase2
    >>> maxProduct([2, 3, 0, 2, 4, 0, 3, 4])
    12

    - testcase3
    >>> maxProduct([0, 5, 3, -1, -2, 0, -2, 4, 0, 3, 4])
    30

    - testcase4
    >>> maxProduct([0, 1, -3, -4, -5, -1, 1, 2, 1, -1, -100, 0, -100000])
    12000

    - testcase5
    >>> maxProduct([0, -3000, -2, 0, -100, -100, 0, -9, -8, 1, 1, 2])
    10000

    - testcase6
    >>> maxProduct([0, -2, 0])
    0
    """

    if A is None or not isinstance(A, list):
    return None
    if len(A) == 1:
    return A[0]
    return -1

    doctest.testmod()
  • Step 3: implement the logic
    It’s time to tackle the most difficult part: write the real function. Think about time complexity (it is best to use only one iteration around the input array which means O(n)), and space complexity (it is best not to use extra space). Run testmod() again and again to find mistakes and modify the codes accordingly. Finally I come with a solution with a helper shadow function _maxProduct. And it passes the six test cases. Althoug I am not sure that this function does not have any bug, at least it works now.
      import doctest
    from sys import maxint

    def maxProduct(A):
    """ A function to find the maximum product value for a continuous subarray.
    :param A: an array or list
    :type A: list

    - testcase1
    >>> maxProduct([0, 2, 3,-2,4, 1, -1])
    48

    - testcase2
    >>> maxProduct([2, 3, 0, 2, 4, 0, 3, 4])
    12

    - testcase3
    >>> maxProduct([0, 5, 3, -1, -2, 0, -2, 4, 0, 3, 4])
    30

    - testcase4
    >>> maxProduct([0, 1, -3, -4, -5, -1, 1, 2, 1, -1, -100, 0, -100000])
    12000

    - testcase5
    >>> maxProduct([0, -3000, -2, 0, -100, -100, 0, -9, -8, 1, 1, 2])
    10000

    - testcase6
    >>> maxProduct([0, -2, 0])
    0
    """

    if A is None or not isinstance(A, list):
    return None
    if len(A) == 1:
    return A[0]
    return max(_maxProduct(A), _maxProduct([a for a in reversed(A)]))

    def _maxProduct(A):
    max_val_forward = 1
    rst = -maxint
    for a in A:
    if a != 0:
    max_val_forward *= a
    rst = max(rst, max_val_forward)
    else:
    rst = max(0, rst)
    max_val_forward = 1
    return rst

    if __name__ == "__main__":
    doctest.testmod()
In conclusion, the most important thing about TDD in data analysis is writing test cases, which really needs a lot of training and exercises.
​

Sunday, August 10, 2014

Translate SAS's sas7bdat format to SQLite and Pandas

Thanks Jared Hobbs’ sas7bdat package, Python can read SAS’s data sets quickly and precisely. And it will be great to have a few extension functions to enhance this package with SQLite and Pandas.
The good things to transfer SAS libraries to SQLite:
  1. Size reduction:
    SAS’s sas7bdat format is verbose. So far successfully loaded 40GB SAS data to SQLite with 85% reduction of disk usage.
  2. Save the cost to buy SAS/ACCESS
    SAS/ACCESS costs around $8,000 a year for a server, while SQLite is accessible for most common softwares.
The good things to transfer SAS data set to Pandas:
  1. Pandas’ powerful Excel interface:
    Write very large Excel file quickly as long as memory can hold data.
  2. Validation of statistics
    Pandas works well with statsmodels and scikit-learn. Easy to validate  SAS’s outputs. 

Tuesday, June 10, 2014

Remove tabs from SAS code files

By default, SAS records the indent by pressing the tab key by tab, which causes many problem to use the code files under a different environment. There are actually two ways to eliminate the tab character in SAS and replace with empty spaces.
  • Regular expression
    Press Ctrl + H → Replace window pops out → Choose Regular expression search → At the box of Find text input \t→ At the box of Replace input multiple\s, say four

  • Editor option
    Click Tools → Options → Enhanced Editors… → Choose Insert spaces for tabs → Choose Replace tabs with spaces on file open
​