Software & Service

Algorithmic pricing with Redis

2017-09-11T23:23:00.000-05:00

The algorithmic pricing is an exciting new area, and it combines engineering and mathematics. Chen’s paper has introduced the algorithmic pricing on Amazon Marketplace. This post is to discuss the implementation of an algorithmic pricing based on Redis from the perspective of the sellers.

The background

Each of Amazon’s ASINs will have many sellers that compete each other.
Amazon has a ranking mechanism for Buy Box, say, to punish a new seller. But for the same ASIN, the seller who has the lowest price usually wins the Buy Box.
For each ASIN, each of many sellers will have an optimal price (the price they want to sell) and a lowest acceptable price (they cannot sell if the price is below it).
Amazon allows a seller to change price for an ASIN every 15 minutes.
Why 15 minutes? Because Amazon said that it takes up to 15 minutes for all systems to converge to the new price. But it is not exactly true. Amazon MWS uses 4 data centers in North America and data synchronization is not that fast for each of the databases or each of the data centers.
Amazon MWS has APIs. But don’t rely on it, and build your own crawlers and price adjustor. The reasons are explained later.
The sellers’ target variables are mainly
- Buy Box take-over percentage (how they beat the other competitors)
- Profit margin (the higher prices they sell goods at, the better)

How an algorithmic pricing engine works

The overall infrastructure may include three parts.

1. The distributed crawler

Amazon most time allows a crawler running 1 QPS/IP. But all the data centers and all ASINs from the seller and his competitors have to be closely watched. So a distributed approach will be safer. The response time will be crucial to decide which one is the winner within the game. In a common chasing diagram below, clearly Seller 1 has the upper hand and Seller 2 is losing the group, since the crawler from Seller 1 is faster.

2. The algorithm core

The embedded algorithm will have two purposes

Analyze the data from the crawler and detect other competitors’ optimal price/lowest acceptable price and strategy.
Adjust price according to an algorithm. The strategy can be either as simple as "minus one cent from the current competitor's price", or as complicated as machine learning or deep learning.

3. The price adjustor

The computed price will enter Amazon and become effective. Amazon will rank lower for the sellers whom he thinks doing algorithmic pricing. Sometimes Amazon even bans the seller. So don’t let Amazon to find that a computer is manipulating its APIs. To emulate human’s behavior on a browser, the two options are phantom.js and headless Chrome.

The role of Redis

Redis is a in-memory data store, which supports persistence and sharding. The first usage for this algorithmic pricing engine is that it can be used as a centralized task scheduler for the crawlers. Besides that, some other interesting fields in algorithmic pricing could be explored and utilized with Redis.

1. Predicative pricing

Redis as a cache has support for TTL(time to live). With the accumulation of the data, the competitors’ price changing time could be predicted. In the publisher-subscriber model, each time the predicted duration for next price changing can be inputted as an expiring key with TTL. Once the key expires, the publisher dispatches a crawling task and makes the price adjustment. The good thing for this approach is that the crawlers don’t need to tap a single web page from Amazon every second that brings the risk of being banned.

+--------------------+
|  Subscriber        |
+--------------------+
|                    |
| + psubscribe():void|  +--------------+
|                    |  |      PubSub  |
+--------------------+  +--------+-----+
                                 ^
                                 |
+--------------------------------+-----+
|       Publisher                      |
+--------------------------------------+
|                                      |
|  + calculateNextTTL():int            |
|  + onPMessage():void                 |
|  + dispatchCrawler():boolean         |
|                                      |
+--------------------------------------+

2. Cross-data-center hedging

The synchronization mechanism across the data centers costs time sometimes hours. That is the reason why we see different prices at different IPs of Amazon at the same time. People will also have different purchasing behavior pattern at different time. Since a seller has the option to change price at a specified data center with IP instead of domain name, it will be an interesting topic to utilize the cool down time for the price to spread to the overall network and make a hedging. Redis’ capacity to keep all prices as hashes in memory will be helpful to spot those valuable occasions.

Good math, bad engineering

2017-02-17T10:38:00.001-06:00

As a formal statistician and a current engineer, I feel that a successful engineering project may require both the mathematician’s ability to find the abstraction and the engineer’s ability to find the implementation.

For a typical engineering problem, the steps are usually -
- 1. Abstract the problem with a formula or some pseudocodes
- 2. Solve the problem with the formula
- 3. Iterate the initial solution until it achieves the optimal time complexity and space complexity

I feel that a mathematician would like dynamic programming or DP questions most, because they are too similar to the typical deduction question in math. An engineer will feel it challenging, since it needs the imagination and some sense of math.

The formula is the most important: without it, try-and-error or debugging does not help. Once the the formula is figured out, the rest becomes a piece of the cake. However, sometimes things are not that straightforward. Good mathematics does not always lead to good engineering.

Let’s see one question from Leetcode.

You are given a list of non-negative integers, a1, a2, ..., an, and a target, S. Now you have 2 symbols + and -. For each integer, you should choose one from + and - as its new symbol.

Find out how many ways to assign symbols to make sum of integers equal to target S.

Example 1:
Input: nums is [1, 1, 1, 1, 1], S is 3. 
Output: 5
Explanation: 

-1+1+1+1+1 = 3
+1-1+1+1+1 = 3
+1+1-1+1+1 = 3
+1+1+1-1+1 = 3
+1+1+1+1-1 = 3

There are 5 ways to assign symbols to make the sum of nums be target 3.

1. The quick solution

For each of the element of a list, it has two options: plus or minus. So the question asks how many ways to get a special number by all possible paths. Of course, if the sum of numbers is unrealistic, we just need to return 0.

Sounds exactly like a DP question. If we have a pencil and a paper, we can start to explore the relationship between dp(n) and dp(n-1). For example, our goal is to get a sum of 5, and we are given a list of [1, 1, 1, 1, 1]. If th a smaller tuple/list is (1, 1, 1, 1) and some paths get 4, that is exactly what we want since it adds 1 and becomes 5. Similarly, if they could get 6, that is fine as well. We add simply both paths together, since there are only two paths.

The formula is is dp(n, s) = dp(n-1, s-x) + dp(n-1, s+x), where n is the size of the list, s is the sum of the numbers and x is the one that adds to the previous list. OK, the second step is easy.

def findTargetSumWays_1(nums, S):
    """
    :type nums: Tuple[int]
    :type S: int
    :rtype: int
    """
    if not nums:
        if S == 0:
            return 1
        else:
            return 0
    return findTargetSumWays_1(nums[1:], S+nums[0]) + findTargetSumWays_1(nums[1:], S-nums[0]) 

small_test_nums = (1, 1, 1, 1, 1)
small_test_S = 3

%time findTargetSumWays_1(small_test_nums, small_test_S)

It is theoretically correct and works perfectly with small test cases. But we know that it is going to a nightmare for an engineering application, because it has a hefty time complexity of O(2^N). So math part is done, and We have to move to the third step.

2. The third step that is hard

So we need to find a data structure to record all the paths. If it is the Fibonacci number problem, a simple linear data structure like a list will slash O(2^N) to O(N).

But the hard part is: what data structure is going to be used here. Since the get operation in Hashtable is O(1), a rolling dictionary will help record the previous states. However, Python’s dictionary does not support change/add ops while it is in a loop, then we have to manually replace it. The overall path will be like a tree structure. So the ideal solution will be like -

def findTargetSumWays_2(nums, S):
    if not nums:
        return 0
    dic = {nums[0]: 1, -nums[0]: 1} if nums[0] != 0 else {0: 2}
    for i in range(1, len(nums)):
        tdic = {}
        for d in dic:
            tdic[d + nums[i]] = tdic.get(d + nums[i], 0) + dic.get(d, 0)
            tdic[d - nums[i]] = tdic.get(d - nums[i], 0) + dic.get(d, 0)
        dic = tdic
    return dic.get(S, 0)

big_test_nums = tuple(range(100))
big_test_S = sum(range(88))
%time findTargetSumWays_2(big_test_nums, big_test_S)

The time is exactly what we need. However, the codes are not elegant and hard to understand.

CPU times: user 189 ms, sys: 4.77 ms, total: 194 ms
Wall time: 192 ms

3. Finally the easy solution

If we don’t want things to get complicated. Here we just want a cache and Python 3 provides a lru_cache decorator. Then adding one line to the first solution will quickly solve the problem.

@lru_cache(10000000)
def findTargetSumWays_3(nums, S):
    if not nums:
        if S == 0:
            return 1
        else:
            return 0
    return findTargetSumWays_3(nums[1:], S+nums[0]) + findTargetSumWays_3(nums[1:], S-nums[0]) 

%time findTargetSumWays_3(big_test_nums, big_test_S)

CPU times: user 658 ms, sys: 19.7 ms, total: 677 ms
Wall time: 680 ms

Conclusion

Good math cannot solve all the engineering problems. It has to combine with the details of the languange, the application and the system to avoid bad engineering implementation.

The Jupyter notebook is at Github. If you have any comment, please email me wm@sasanalysis.com.

Use Elasticsearch and Kibana for large BI system

2017-02-16T13:20:00.001-06:00

Nowadays Elasticsearch is more and more popular. Besides it original search functionalities, I found Elasticsearch can be

used as a logging container. That is what the ELK stack is created for.
utilized as a JSON server with richful APIs, which can be combined with its Kibana as BI servers.

That is the data store I see everyday

10PB stocking data
average 30TB incoming data everyday
various data sources including binary files such PDF
including very complicated SQL queries (fortunately no stored procedures)
millions of JSON creations daily

People want to know what is going on with such data. So a business intelligence or an OLAP system is needed to visualize/aggregate the data and its flow. Since Elasticsearch is so easy to scale out, it beats other solutions for big data on the market.

1. Batch Worker

There are many options to implement a batch worker. Finally the decision falls to either Spring Data Batch or writing a library from the scratch in Python.

1.1 Spring Data Batch v.s. Python

Spring Data Batch

Pros:
- is a full framework that includes rollback, notification and scheduler features
- provides great design pattern for Dependency Injection, such as factory and singleton, which help multiple persons work together. For example -

@Bean
public Step IndexMySQLJob01() {
   return stepBuilderFactory.get("IndexMySQLJob01")
           .<Data, Data> chunk(10)
           .reader(reader())
           .processor(processor())
           .writer(writer())
           .build();
}

Cons:
- is difficult to handle nested JSONs
- Elasticsearch Java client is not as good as its Python client

Python

Pros:
- less codes; flexible
- super easy from dictionary to JSON
- has official/3rd party libraries for everything
Cons:
- you create your own library/framework
- if the pattern like Spring Data Batch is deployed, has to inject dependencies manually, such as -

class IndexMySQLJob01(object):
     def __init__(self, reader, processor, writer, listener):
          self.reader = reader
          self.processor = processor
          self.writer = writer
          self.listener = listener
     ...

Eventually Python is picked, because the overall scenario is more algorithm-bound instead of language-bound.

1.2 The algorithms

Since the data size is pretty big, time and space are always considered. The direct way to decrease the time complexity is using the hash tables, as long as the memory can hold the data. For example, a join between an N rows table and an M rows table can be optimized from O(M*N) to O(M).

To save the space, a generator chain is used to stream data from the start to the end, instead of materializing sizable objects.

class JsonTask01(object):
    ...
    def get_json(self, generator1, hashtable1):
        for each_dict in generator1:
            key = each_dict.get('key')
            each_dict.update(hashtable1.get(key))
            yield each_dict

1.3 The scheduler

A scheduler is a must: cron is enough for simple tasking, while a bigger system requires a work flow. Airflow is the one that helps organize and schedule. It has a web UI and is written in Python, which is easy to be integrated with the batch worker.

1.4 High availability with zero downtime

Indexing of large quantity of data will impose significant impact. For mission-critical indexes that need 100% up time, the zero down algorithm is implemented and we keep two copies of an index for maximum safety. The alias will switch between the two copies once the indexing is finished.

    def add_alias(self, idx):
        LOGGER.warn("The alias {} will point to {}.".format(self.index, idx))
        self.es.indices.put_alias(idx, self.index)

    def delete_alias(self, idx):
        LOGGER.warn("The alias {} will be removed from {}.".format(self.index, idx))
        self.es.indices.delete_alias(idx, self.index)

2. Elasticsearch Cluster

2.1 Three kinds of nodes

An Elasticsearch node can choose one of three roles: master node, data node and ingest node(previously called client node). It is commonly seen to dedicate a node as master and ingest and data all together. For a large system, it is always helpful to assign the three roles to different machines/VMs. Therefore, once a node is down/up, it will be quicker to failover or recover.
elasticsearch-head can clearly visualize the data transfer process of the shards once an accident occurs.

With the increased number of cluster nodes, the deployment becomes painful. I feel the best tool so far is ansible-elasticsearch. With ansible-playbook -i hosts ./your-playbook.yml -c paramiko, the cluster is on the fly.

2.2 Memory is cheap but heap is expensive

The rules of thumb for Elasticsearch are -

Give (less than) Half Your Memory to Lucene

Don’t Cross 32 GB!

The result causes an awkward situation: if you have a machine that has more than 64GB memory, then the additional memory will mean nothing to Elasticsearch. Actually it is meaningful to run two or more Elasticsearch instances side by side to save the hardware. For example, there is a machine with 96GB memory. We can allocate 31GB for an ingest node, 31 GB for a data node and the rest for the OS. However, two data nodes in a single machine will compete for the disk IO that damages the performance, while a master node and a data node will increase the risk of downtime.

The great thing for Elasticsearch is that it provides richful REST APIs, such as http://localhost:9200/_nodes/stats?pretty. We could use Xpack(paid) or other customized tools to monitor them. I feel that the three most important statistics for the heap and therefore the performance are -

The heap usage and the old GC duration

The two statistics intertwined together. The high heap usage, such as 75%, will lead to a GC, while GC with high heap usage will take longer time. We have to keep both numbers as low as possible.

     "jvm" : {
        "mem" : {
          "heap_used_percent" : 89,
          ... 
        "gc" : {
          "collectors" : {
            ...
            "old" : {
              "collection_count" : 225835,
              "collection_time_in_millis" : 22624857
            }
          }
        }

Thread pools

There are three kinds of thread pools: active, queue, and reject. It is useful to visualize the real time change. Once there are a lot of queued threads or rejected threads, it is good time to think about scale up or scale out.

"thread_pool" : {
        "bulk" : {
          "threads" : 4,
          "queue" : 0,
          "active" : 0,
          "rejected" : 0,
          "largest" : 4,
          "completed" : 53680
        },

The segments’ number/size

The segments are the in-memory inverted indexes corresponding to the indexes on the hard disk, which are persistent in the physical memory and GC will have no effect on them. The segment will have the footage on every search thread. The size of the segments are important because they will be multiplied by a factor of the number of threads.

        "segments" : {
          "count" : 215,
          "memory_in_bytes" : 15084680,
        },

The number of shards actually controls the number of the segments. The shards increase, then the size of the segments decreases and the number of the segments increases. So we cannot increase the number of shards as many as we want. If there are many small segments, the heap usage will turn much higher. The solution is Force merge, which is time-consuming but effecitve.

3. Kibana as BI dashboard

3.1 Plugins

Kibana integreated DevTools(previously The Sense UI) for free. DevTools has code assistance and is a powerful tool for debugging. If the budget is not an issue, Xpack is also highly recommended. As for Elasticsearch, since 5.0, ingest-geoip is now a plugin. We will have to write it to Ansible YAML such as -

es_plugins:
    - plugin: ingest-geoip

3.2 Instant Aggregation

There are quite a few KPIs that need system-wide term aggregations. From 5.0 the request cache will be enabled by default for all requests with size:0.
For example -

POST /big_data_index/data/_search
{   "size": 0,
    "query": {
        "bool": {
            "must_not": {
                "exists": {
                    "field": "interesting_field"
                }
            }
        }
    }
}

The Fore merge as mentioned above, such asPOST /_forcemerge?max_num_segments=1, will combine the segments and dramatically increase the aggregation speed.

3.3 Proxy

Nginx is possibly the best proxy as the frontend toward Kibana. There are two advantages: first the proxy can cache the static resources of Kibana; second we can always check the Nginx logs to figure out what causes problem for Kibana.

Conclusion

Elasticsearch and Kibana together provide high availability and high scalability for large BI system.

if you have any comment, please email me wm@sasanalysis.com

Use Slack bot to monitor the server

2016-09-07T14:19:00.003-05:00

I used to install Datadog or other SaaS to monitor my Linux boxes on the cloud. Most times they are just overkill for my tiny servers with only 1GB or 2GB memory. Actually what I am most interested is the up-and-running processes, or/and the exact memory usage. And I need a mobile solution to monitor on-the-go.

Now with the coming of Slack bot, and its real time Python client, I can just use a simple Python script to realize the purposes.

from slackclient import SlackClient
from subprocess import  getoutput
import logging
import time

message_channel = '#my-server-001'
api_key = 'xoxb-slack-token'
client = SlackClient(api_key)

if client.rtm_connect():
    while True:
        last_read = client.rtm_read()
        if last_read:
            try:
                parsed = last_read[0]['text']
                if parsed and 'status' in parsed:
                    result = getoutput('pstree')
                    result += '\n\n' + getoutput('free -h')
                    client.rtm_send_message(message_channel, str(result))
            except Exception as e:
                logging.error(e)
        time.sleep(1)

Then I use systemd or other tools to daemonize it. No matter where and when I am, I enter status at the #my-server-001 channel on my phone, I will instantly get the result like -

systemd-+-accounts-daemon-+-{gdbus}
       |                 `-{gmain}
       |-agetty
       |-cron
       |-dbus-daemon
       |-fail2ban-server---2*[{fail2ban-server}]
       |-login---bash
       |-nginx---nginx
       |-postgres---5*[postgres]
       |-python---sh---pstree
       |-redis-server---2*[{redis-server}]
       |-rsyslogd-+-{in:imklog}
       |          |-{in:imuxsock}
       |          `-{rs:main Q:Reg}
       |-sshd-+-3*[sshd---sshd---bash]
       |      `-sshd---sshd
       |-2*[systemd---(sd-pam)]
       |-systemd-journal
       |-systemd-logind
       |-systemd-timesyn---{sd-resolve}
       |-systemd-udevd
       `-uwsgi---uwsgi---5*[uwsgi]

             total        used        free      shared  buff/cache   available
Mem:           2.0G        207M        527M         26M        1.2G        1.7G
Swap:          255M          0B        255M

Deploy edx spark environment to DigitalOcean

2015-07-29T08:24:00.000-05:00

This summer I took the Spark courses at edx CS100 and CS190, and had wonderful experience.

The two classes apply a Vagrant virtual machine containing Spark and all teaching materials. There are two challenges with the virtual machine —

The labs usually take long time to finish, say 8-10 hours. If the host machine is closed, the RDDs will be lost and the pipeline has to be run again.
Some RDD operations take a lot computation/communication powers, such as groupByKey and distinct. Many of my 50k classmates complained about the waiting time. And my most used laptop is a Chromebook and doesn’t even have options to install Virtual Box.

To deploy the learning environment to a cloud may be an alternative. DigitalOcean is a good choice because it uses mirrors for most packages, and the network speed is amazingly fast that is almost 100MB/s (thanks to the SSD infrastructure DigitalOcean implements for the cloud, otherwise the hard disk may not stand this rapid IO; see my deployment records GitHub).

I found that a Linux box with 1 GB memory and 1 CPU at DigitalOcean that costs 10 dollars a month will handle most labs fairly easy with IPython and Spark. A 2 GB memory and 2 CPU droplet will be ideal since it is the minimal requirement for a simulated cluster. It costs 20 dollars a month, but is still much cheaper than the cost to earn the big data certificate that is $100 (50 for each). I just need to write Python scripts to install IPython notebook with SSL, and download Spark and the course materials.

The DevOps tool is Fabric and the fabfile is at GitHub.
The deployment pipeline is also at GitHub

Transform SAS files to Parquet through Spark

2015-07-17T11:25:00.000-05:00

The demo pipeline is at GitHub.

Since the version 1.3, Spark has introduced the new data structure DataFrame. A data analyst now could easily scale out the exsiting codes based on the DataFrame from Python or R to a cluster hosting Hadoop and Spark.

There are quite a few practical scenarios that DataFrame fits well. For example, a lot of data files including the hardly read SAS files want to merge into a single data store. Apache Parquet is a popular column store in a distributed environment, and especially friendly to structured or semi-strucutred data. It is an ideal candidate for a univeral data destination.

I copy three SAS files called prdsale, prdsal2 and prdsal3, which are about a simulated sales record, from the SASHELP library to a Linux directory. And then I launch the SQL context from Spark 1.4.

The three SAS files now have the size of 4.2MB. My overall strategy is to build a pipeline to realize my purpose such as SAS --> Python --> Spark --> Parquet.

import os
try:
    import sas7bdat
    import pandas
except ImportError:
    print('try to install the packags first')

print('Spark verion is {}'.format(sc.version))

if type(sqlContext) != pyspark.sql.context.HiveContext:
    print('reset the Spark SQL context')

os.chdir('/root/playground')

def print_bytes(filename):
    print('{} has {:,} bytes'.format(filename, os.path.getsize(filename)))

print_bytes('prdsale.sas7bdat')
print_bytes('prdsal2.sas7bdat')
print_bytes('prdsal3.sas7bdat')

!du -ch --exclude=test_parquet

Spark verion is 1.4.0
prdsale.sas7bdat has 148,480 bytes
prdsal2.sas7bdat has 2,790,400 bytes
prdsal3.sas7bdat has 1,401,856 bytes
4.2M    .
4.2M    total

1. Test DataFrame in Python and Spark

First I transform a SAS sas7bdat file to a pandas DataFrame. The great thing in Spark is that a Python/pandas DataFrame could be translated to Spark DataFrame by the createDataFrame method. Now I have two DataFrames: one is a pandas DataFrame and the other is a Spark DataFrame.

with sas7bdat.SAS7BDAT('prdsale.sas7bdat') as f:
     pandas_df = f.to_data_frame()
print('-----Data in Pandas dataframe-----')
print(pandas_df.head())

print('-----Data in Spark dataframe-----')
spark_df = sqlContext.createDataFrame(pandas_df)
spark_df.show(5)

-----Data in Pandas dataframe-----
   ACTUAL COUNTRY   DIVISION  MONTH  PREDICT   PRODTYPE PRODUCT  QUARTER  \
0     925  CANADA  EDUCATION  12054      850  FURNITURE    SOFA        1   
1     999  CANADA  EDUCATION  12085      297  FURNITURE    SOFA        1   
2     608  CANADA  EDUCATION  12113      846  FURNITURE    SOFA        1   
3     642  CANADA  EDUCATION  12144      533  FURNITURE    SOFA        2   
4     656  CANADA  EDUCATION  12174      646  FURNITURE    SOFA        2   

  REGION  YEAR  
0   EAST  1993  
1   EAST  1993  
2   EAST  1993  
3   EAST  1993  
4   EAST  1993  
-----Data in Spark dataframe-----
+------+-------+---------+-------+-------+---------+-------+-------+------+------+
|ACTUAL|COUNTRY| DIVISION|  MONTH|PREDICT| PRODTYPE|PRODUCT|QUARTER|REGION|  YEAR|
+------+-------+---------+-------+-------+---------+-------+-------+------+------+
| 925.0| CANADA|EDUCATION|12054.0|  850.0|FURNITURE|   SOFA|    1.0|  EAST|1993.0|
| 999.0| CANADA|EDUCATION|12085.0|  297.0|FURNITURE|   SOFA|    1.0|  EAST|1993.0|
| 608.0| CANADA|EDUCATION|12113.0|  846.0|FURNITURE|   SOFA|    1.0|  EAST|1993.0|
| 642.0| CANADA|EDUCATION|12144.0|  533.0|FURNITURE|   SOFA|    2.0|  EAST|1993.0|
| 656.0| CANADA|EDUCATION|12174.0|  646.0|FURNITURE|   SOFA|    2.0|  EAST|1993.0|
+------+-------+---------+-------+-------+---------+-------+-------+------+------+

The two should be the identical length. Here both show 1,440 rows.

print(len(pandas_df))
print(spark_df.count())

1440
1440

2. Automate the transformation

I write a pipeline function to automate the transformation. As the result, the all three SAS files are saved to the same directory as Parquet format.

def sas_to_parquet(filelist, destination):
    """Save SAS file to parquet
    Args:
        filelist (list): the list of sas file names
        destination (str): the path for parquet
    Returns:
        None
    """
    rows = 0
    for i, filename in enumerate(filelist):
        with sas7bdat.SAS7BDAT(filename) as f:
            pandas_df = f.to_data_frame()
            rows += len(pandas_df)
        spark_df = sqlContext.createDataFrame(pandas_df)
        spark_df.save("{0}/key={1}".format(destination, i), "parquet")
    print('{0} rows have been transformed'.format(rows))

sasfiles = [x for x in os.listdir('.') if x[-9:] == '.sas7bdat']
print(sasfiles)

sas_to_parquet(sasfiles, '/root/playground/test_parquet')

['prdsale.sas7bdat', 'prdsal2.sas7bdat', 'prdsal3.sas7bdat']
36000 rows has been transformed

Then I read from the newly created Parquet data store. The query shows that the data has been successfully saved.

df = sqlContext.load("/root/playground/test_parquet", "parquet")
print(df.count())
df.filter(df.key == 0).show(5)

36000
+------+-------+------+----+-------+-------+---------+-------+-------+-----+------+-----+---------+------+---+
|ACTUAL|COUNTRY|COUNTY|DATE|  MONTH|PREDICT| PRODTYPE|PRODUCT|QUARTER|STATE|  YEAR|MONYR| DIVISION|REGION|key|
+------+-------+------+----+-------+-------+---------+-------+-------+-----+------+-----+---------+------+---+
| 925.0| CANADA|  null|null|12054.0|  850.0|FURNITURE|   SOFA|    1.0| null|1993.0| null|EDUCATION|  EAST|  0|
| 999.0| CANADA|  null|null|12085.0|  297.0|FURNITURE|   SOFA|    1.0| null|1993.0| null|EDUCATION|  EAST|  0|
| 608.0| CANADA|  null|null|12113.0|  846.0|FURNITURE|   SOFA|    1.0| null|1993.0| null|EDUCATION|  EAST|  0|
| 642.0| CANADA|  null|null|12144.0|  533.0|FURNITURE|   SOFA|    2.0| null|1993.0| null|EDUCATION|  EAST|  0|
| 656.0| CANADA|  null|null|12174.0|  646.0|FURNITURE|   SOFA|    2.0| null|1993.0| null|EDUCATION|  EAST|  0|
+------+-------+------+----+-------+-------+---------+-------+-------+-----+------+-----+---------+------+---+

3. Conclusion

There are multiple advantages to tranform data from various sources to Parquet.

It is an open format that could be read and written by major softwares.
It could be well distributed to HDFS.
It compresses data.

For example, the original SAS files add up to 4.2 megabyte. Now as Parquet, it only weighs 292KB and achieves 14X compression ratio.

os.chdir('/root/playground/test_parquet/')
!du -ahc 

4.0K    ./key=2/._metadata.crc
4.0K    ./key=2/._SUCCESS.crc
0    ./key=2/_SUCCESS
4.0K    ./key=2/_common_metadata
4.0K    ./key=2/.part-r-00001.gz.parquet.crc
4.0K    ./key=2/._common_metadata.crc
4.0K    ./key=2/_metadata
60K    ./key=2/part-r-00001.gz.parquet
88K    ./key=2
4.0K    ./key=0/._metadata.crc
4.0K    ./key=0/._SUCCESS.crc
0    ./key=0/_SUCCESS
4.0K    ./key=0/_common_metadata
4.0K    ./key=0/.part-r-00001.gz.parquet.crc
4.0K    ./key=0/._common_metadata.crc
4.0K    ./key=0/_metadata
12K    ./key=0/part-r-00001.gz.parquet
40K    ./key=0
4.0K    ./key=1/._metadata.crc
4.0K    ./key=1/._SUCCESS.crc
0    ./key=1/_SUCCESS
4.0K    ./key=1/_common_metadata
4.0K    ./key=1/.part-r-00001.gz.parquet.crc
4.0K    ./key=1/._common_metadata.crc
4.0K    ./key=1/_metadata
132K    ./key=1/part-r-00001.gz.parquet
160K    ./key=1
292K    .
292K    total

A bar plot visualizes the signifcant size difference between the two formats. It shows an order of magnitude space deduction.

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
index = np.arange(2)
bar_width = 0.35
data = [4200, 292]
header = ['SAS files', 'Parquet']
plt.bar(index, data)
plt.grid(b=True, which='major', axis='y')
plt.ylabel('File Size by KB')
plt.xticks(index + bar_width, header)
plt.show()

Two alternative ways to query large dataset in SAS

2015-06-19T10:43:00.000-05:00

I really appreciate those wonderful comments on my SAS posts by the readers (1, 2, 3). They gave me a lot of inspirations. Due to SAS or SQL’s inherent limitation, recently I feel difficult in deal with some extremely large SAS datasets (it means that I exhausted all possible traditional ways). Here I conclude two alternative solutions in these extreme cases as a follow-up to the comments.

Read Directly
- Use a scripting language such as Python to Reading SAS datasets directly
Code Generator
- Use SAS or other scripting languages to generate SAS/SQL codes

The examples still use sashelp.class, which has 19 rows. The target variable is weight.

*In SAS
data class;
    set sashelp.class;
    row = _n_;
run;

Example 1: Find the median

SQL Query

In the comment, Anders SköllermoFebruary wrote

Hi! About 1. Calculate the median of a variable:
If you look at the details in the SQL code for calculation the median, then you find that the intermediate file is of size N*N obs, where N is the number of obs in the SAS data set.
So this is OK for very small files. But for a file with 10000 obs, you have an intermediate file of size 100 million obs. / Br Anders Anders Sköllermo Ph.D., Reuma and Neuro Data Analyst

The SQL query below is simple and pure, so that it can be ported to any other SQL platform. However, just like what Anders said, it is just way too expensive.

*In SAS
proc sql;
   select avg(weight) as Median
   from (select e.weight
       from class e, class d
       group by e.weight
       having sum(case when e.weight = d.weight then 1 else 0 end)
          >= abs(sum(sign(e.weight - d.weight)))
    );
quit;

PROC UNIVARIATE

In the comment, Anonymous wrote:
I noticed the same thing - we tried this on one of our 'smaller' datasets (~2.9 million records), and it took forever.
Excellent solution, but maybe PROC UNIVARIATE will get you there faster on a large dataset.
Indeed PROC UNIVARIATE is the best solution in SAS to find the median, which utilizes SAS's built-in powers.

Read Directly

When the extreme cases come, say SAS cannot even open the entire dataset, we may have to use the streaming method to Reading the sas7bdat file line by line. The sas7bdat format has been decoded by Java, R and Python. Theoretically we don't need to have SAS to query a SAS dataset.

Heap is an interesting data structure, which easily finds a min or a max. ream the values, we could build a max heap and a min heap to cut the incoming stream into half in Python. The algorithm looks like a heap sorting. The good news is that it only Reading one variable each time and thus saves a lot of space.

#In Python
import heapq
from sas7bdat import SAS7BDAT
class MedianStream(object):
    def __init__(self):
        self.first_half = [] # will be a max heap
        self.second_half = [] # will be a min heap, 1/2 chance has one more element
        self.N = 0

    def insert(self, x):
        heapq.heappush(self.first_half, -x)
        self.N += 1
        if len(self.second_half) == len(self.first_half):
            to_second, to_first = map(heapq.heappop, [self.first_half, self.second_half])
            heapq.heappush(self.second_half, -to_second)
            heapq.heappush(self.first_half, -to_first)
        else:
            to_second = heapq.heappop(self.first_half)
            heapq.heappush(self.second_half, -to_second)

    def show_median(self):
        if self.N == 0:
            raise IOError('please use the insert method first')
        elif self.N % 2 == 0:
            return (-self.first_half[0] + self.second_half[0]) / 2.0
        return -self.first_half[0]

if __name__ == "__main__": 
    stream = MedianStream()
    with SAS7BDAT('class.sas7bdat') as infile:
        for i, line in enumerate(infile):
            if i == 0:
                continue
            stream.insert(float(line[-1]))
    print stream.show_median()

99.5

Example 2: Find top K by groups

SQL Query

This query below is very expensive. We have a self-joining O(N^2) and a sorting O(NlogN), and the total time complexity is a terrible O(N^2 + Nlog(N)).

* In SAS
proc sql; 
    select a.sex, a.name, a.weight, (select count(distinct b.weight) 
            from class as b where b.weight >= a.weight and a.sex = b.sex ) as rank 
    from class as a
    where calculated rank <= 3
    order by sex, rank
;quit;

Code Generator

The overall thought is break-and-conquer. If we synthesize SAS codes from a scripting tool such as Python, we essentially get many small SAS codes segments. For example, the SQL code below is just about sorting. So the time comlexity is largely decreased to O(Nlog(N)).

# In Python
def create_sql(k, candidates):
    template = """
    proc sql outobs = {0};
    select *
    from {1}
    where sex = '{2}' 
    order by weight desc
    ;
    quit;"""
    for x in candidates:
        current = template.format(k, 'class', x)
        print current
if __name__ == "__main__":
    create_sql(3, ['M', 'F'])


    proc sql outobs = 3;
    select *
    from class
    where sex = 'M' 
    order by weight desc
    ;
    quit;

    proc sql outobs = 3;
    select *
    from class
    where sex = 'F' 
    order by weight desc
    ;
    quit;

Read Directly

This time we use the data structure of heap again in Python. To find the k top rows for each group, we just need to prepare the min heaps with the k size for each group. With the smaller values popped out everytime, we finally get the top k values for each group. The optimized time complexity is O(Nlog(k))

#In Python
from sas7bdat import SAS7BDAT
from heapq import heappush, heappop

def get_top(k, sasfile):
    minheaps = [[], []]
    sexes = ['M', 'F']
    with SAS7BDAT(sasfile) as infile:
        for i, row in enumerate(infile):
            if i == 0:
                continue
            sex, weight = row[1], row[-1]
            i = sexes.index(sex)
            current = minheaps[i]
            heappush(current, (weight, row))
            if len(current) > k:
                heappop(current)
    for x in minheaps:
        for _, y in x:
            print y

if __name__ == "__main__":
    get_top(3, 'class.sas7bdat')


[u'Robert', u'M', 12.0, 64.8, 128.0]
[u'Ronald', u'M', 15.0, 67.0, 133.0]
[u'Philip', u'M', 16.0, 72.0, 150.0]
[u'Carol', u'F', 14.0, 62.8, 102.5]
[u'Mary', u'F', 15.0, 66.5, 112.0]
[u'Janet', u'F', 15.0, 62.5, 112.5]

Example 3: Find Moving Window Maxium

At the daily work, I always want to find three statistics for a moving window: mean, max, and min. The sheer data size poses challenges.
In his blog post, Liang Xie showed three advanced approaches to calculated the moving averages, including `PROC EXPAND`, `DATA STEP` and `PROC SQL`. Apparently `PROC EXPAND` is the winner throughout the comparison. As conclusion, self-joining is very expensive and always O(N^2) and we should avoid it as much as possible.
The question to find the max or the min is somewhat different other than to find the mean, since for the mean only the mean is memorized, while for the max/min the locations of the past min/max should also be memorized.

Code Generator

The strategy is very straightforward: we choose three rows from the table sequentially and calculate the means. The time complexity is O(k*N). The generated SAS code is very lengthy, but the machine should feel comfortable to Reading it.

In addition, if we want to save the results, we could insert those maximums to an empty table.

# In Python
def create_sql(k, N):
    template = """
    select max(weight)
    from class
    where row in ({0})
    ;"""
    SQL = ""
    for x in range(1, N - k + 2):
        current = map(str, range(x, x + 3))
        SQL += template.format(','.join(current))
    print "proc sql;" + SQL + "quit;"
if __name__ == "__main__":
    create_sql(3, 19)


proc sql;
    select max(weight)
    from class
    where row in (1,2,3)
    ;
    select max(weight)
    from class
    where row in (2,3,4)
    ;
    select max(weight)
    from class
    where row in (3,4,5)
    ;
    select max(weight)
    from class
    where row in (4,5,6)
    ;
    select max(weight)
    from class
    where row in (5,6,7)
    ;
    select max(weight)
    from class
    where row in (6,7,8)
    ;
    select max(weight)
    from class
    where row in (7,8,9)
    ;
    select max(weight)
    from class
    where row in (8,9,10)
    ;
    select max(weight)
    from class
    where row in (9,10,11)
    ;
    select max(weight)
    from class
    where row in (10,11,12)
    ;
    select max(weight)
    from class
    where row in (11,12,13)
    ;
    select max(weight)
    from class
    where row in (12,13,14)
    ;
    select max(weight)
    from class
    where row in (13,14,15)
    ;
    select max(weight)
    from class
    where row in (14,15,16)
    ;
    select max(weight)
    from class
    where row in (15,16,17)
    ;
    select max(weight)
    from class
    where row in (16,17,18)
    ;
    select max(weight)
    from class
    where row in (17,18,19)
    ;quit;

Read Directly

Again, if we want to further decrease the time complexity, say O(N), we have to use better data structure, such as queue. SAS doesn't have queue, so we may switch to Python. Actually it has two loops which adds up to O(2N). However, it is still better than any other methods.

# In Python
from sas7bdat import SAS7BDAT
from collections import deque

def maxSlidingWindow(A, w):
    N = len(A)
    ans =[0] * (N - w + 1)
    myqueue = deque()
    for i in range(w):
        while myqueue and A[i] >= A[myqueue[-1]]:
            myqueue.pop()
        myqueue.append(i)
    for i in range(w, N):
        ans[i - w] = A[myqueue[0]]
        while myqueue and A[i] >= A[myqueue[-1]]:
            myqueue.pop()
        while myqueue and myqueue[0] <= i-w:
            myqueue.popleft()
        myqueue.append(i)
    ans[-1] = A[myqueue[0]]
    return ans

if __name__ == "__main__":
    weights = []
    with SAS7BDAT('class.sas7bdat') as infile:
        for i, row in enumerate(infile):
            if i == 0:
                continue
            weights.append(float(row[-1]))

    print maxSlidingWindow(weights, 3)

[112.5, 102.5, 102.5, 102.5, 102.5, 112.5, 112.5, 112.5, 99.5, 99.5, 90.0, 112.0, 150.0, 150.0, 150.0, 133.0, 133.0]

Conclusion

While data is expanding, we should more and more consider three things -
Time complexity: we don't want run data for weeks.
Space complexity: we don't want the memory overflow.
Clean codes: the colleagues should easily Reading and modify the codes.

saslib: a simple Python tool to lookup SAS metadata

2015-06-03T10:59:00.000-05:00

saslib is an HTML report generator to lookup the metadata (or the head information) like PROC CONTENTS in SAS.

It reads the sas7bdat files directly and quickly, and does not need SAS installed.
Emulate PROC CONTENTS by jQuery and DataTables.
Extract the meta data from all SAS7bdat files under the specified directory.
Support IE(>=10), firefox, chrome and any other modern browser.

Installation

pip install saslib

saslib requires sas7bdat and jinjia2.

Usage

The module is very simple to use. For example, the SAS data sets under the SASHELP library could be viewed —

from saslib import PROCcontents

sasdata = PROCcontents('c:/Program Files/SASHome/SASFoundation/9.3/core/sashelp')
sasdata.show()

The resulting HTML file from the codes above will be like here.

Deploy a minimal Spark cluster

2015-03-20T13:52:00.001-05:00

Requirements

Since Spark is rapidly evolving, I need to deploy and maintain a minimal Spark cluster for the purpose of testing and prototyping. A public cloud is the best fit for my current demand.

Intranet speed
The cluster should easily copy the data from one server to another. MapReduce always shuffles a large chunk of data throughout the HDFS. It’s best that the hard disk is SSD.
Elasticity and scalability
Before scaling the cluster out to more machines, the cloud should have some elasticity to size up or size down.
Locality of Hadoop
Most importantly, the Hadoop cluster and the Spark cluster should have one-to-one mapping relationship like below. The computation and the storage should always be on the same machines.

Hadoop	Cluster Manager	Spark	MapReduce
Name Node	Master	Driver	Job Tracker
Data Node	Slave	Executor	Task Tracker

Choice of public cloud:

I simply compare two cloud service provider: AWS and DigitalOcean. Both have nice Python-based monitoring tools(Boto for AWS and python-digitalocean for DigitalOcean).

From storage to computation
Hadoop’s S3 is a great storage to keep data and load it into the Spark/EC2 cluster. Or the Spark cluster on EC2 can directly read S3 bucket such as s3n://file (the speed is still acceptable). On DigitalOcean, I have to upload data from local to the cluster’s HDFS.
DevOps tools:
- AWS: spark-ec2.py
  - With default setting after running it, you will get
    - 2 HDFSs: one persistent and one ephemeral
    - Spark 1.3 or any earlier version
    - Spark’s stand-alone cluster manager
  - A minimal cluster with 1 master and 3 slaves will be consist of 4 m1.xlarge EC2 instances
    - Pros: large memory with each node having 15 GB memory
    - Cons: not SSD; expensive (cost $0.35 * 6 = $2.1 per hour)
- DigitalOcean: https://digitalocean.mesosphere.com/
  - With default setting after running it, you will get
    - HDFS
    - no Spark
    - Mesos
    - OpenVPN
  - A minimal cluster with 1 master and 3 slaves will be consist of 4 2GB/2CPUs droplets
    - Pros: as low as $0.12 per hour; Mesos provide fine-grained control over the cluster(down to 0.1 CPU and 16MB memory); nice to have VPN to guarantee the security
    - Cons: small memory(each has 2GB memory); have to install Spark manually

Add Spark to DigitalOcean cluster

Tom Faulhaber has a quick bash script for deployment. To install Spark 1.3.0, I write it into a fabfile for Python’s Fabric.
Then all the deployment onto the DigitOcean is just one command line.

# 10.1.2.3 is the internal IP address of the master
fab -H 10.1.2.3 deploy_spark

The source codes above are available at my Github

Solve the Top N questions in SAS/SQL

2015-02-03T08:33:00.000-06:00

This is a following post after my previous post about SAS/SQL.

SAS’s SQL procedure has a basic SQL syntax. I found that the most challenging work is to use PROC SQL to solve the TOP N (or TOP N by Group) questions. Comparing with other modern database systems, PROC SQL is lack of -

The ranking functions such as RANK() or the SELECT TOP clause such as
```
select TOP 3 * 
from class
;
```

The partition by clause such as

select sex, name, weight
from (select sex, name, max(weight) over(partition by sex) max_weight
    from class)
where  weight = max_weight
;

However, there are always some alternative solutions in SAS. I list a few question from an ascending difficulty below to explore the possibilities.

Prepare the data

First a SASHELP.CLASS dataset is used as a demo (availabe for every SAS copy). It is a small weight and height dataset from a faked class of 19 children. Now I only keep the weight variable as target column.

data class;
    set sashelp.class;
    keep name sex weight;
run;

proc sort;
    by descending weight;
run;

Name	Sex	Age	Weight
Philip	M	16	150
Ronald	M	15	133
Robert	M	12	128
Alfred	M	14	112.5
Janet	F	15	112.5
Mary	F	15	112
William	M	15	112
Carol	F	14	102.5
Henry	M	14	102.5
John	M	12	99.5
Barbara	F	13	98
Judy	F	14	90
Thomas	M	11	85
Jane	F	12	84.5
Alice	F	13	84
Jeffrey	M	13	84
James	M	12	83
Louise	F	12	77
Joyce	F	11	50.5

1. Select highest value

It is straightforward to use the outobs option at the begining to single out the highest weight.

title "Select highest weight overall";
proc sql outobs = 1;
    select name, weight
    from class
    order by weight desc
;quit;

Name	Weight
Philip	150

2. Select second highest value

How about the second highest weight? The logic is simple — if we remove the highest weight first, then the second highest weight will take the first row.

title "Select second highest weight overall";
proc sql outobs = 1;
    select name, weight 
    from class
    where weight not in (select max(weight) from class)
    order by weight desc
;quit;

Name	Weight
Ronald	133

3. Select Nth highest value

Now it comes to the hard part. How about the Nth highest value, say, the fourth highest weight? Now we have to do a self-joining to let the distinct value point to 3. Since there are two children with the weight 112.5, the query returns the two tied names.

title "Select Nth highest weight";
%let n = 4;

proc sql;
    select distinct a.name, a.weight
    from class as a
    where (select count(distinct b.weight)
        from class as b
        where b.weight > a.weight
        ) = &n - 1;
quit;

Name	Weight
Alfred	112.5
Janet	112.5

4. Select highest values by group

There are two groups Male and Female in the class, and the easiest way to find the highest weight for each category is select max for female union select max for male. However, a more scalable solution is to use the group by clause that fits more than two groups.

title "Select highest weights by group";
proc sql;
    select sex, name, weight
    from class
    group by sex
    having weight = max(weight)
;quit;

Sex	Name	Weight
F	Janet	112.5
M	Philip	150

5. Rank all values

The ultimate solution to solve all the question above is to derive a rank column for the target. There are two solutions: the first one use a subquery in the select clause, while the second one utilizes a subquery in the where clause.

The subquery in the first solution is independent to the main query, which uses less codes and is easier to recall in practice. The second one is actually a self-joining that is faster than the first solution.

/* Solution I */
proc sql; 
    select name, weight, (select count(distinct b.weight) 
            from class as b where b.weight >= a.weight) as Rank
    from class as a
    order by rank
;quit;

/* Solution II */
proc sql;
    select a.name, a.weight, count(b.weight) as rank
    from class as a, (select distinct weight
           from class
           ) as b
    where a.weight <= b.weight
    group by a.name, a.weight
    order by a.weight desc
;quit;

Name	Weight	Rank
Philip	150	1
Ronald	133	2
Robert	128	3
Alfred	112.5	4
Janet	112.5	4
Mary	112	5
William	112	5
Henry	102.5	6
Carol	102.5	6
John	99.5	7
Barbara	98	8
Judy	90	9
Thomas	85	10
Jane	84.5	11
Alice	84	12
Jeffrey	84	12
James	83	13
Louise	77	14
Joyce	50.5	15

6. Select top N values by group

Once with the rank column at hand, many perplexing problems could be easily solved. For example, we can use it to find the top 3 heaviest people for each category of male and female. And it is also scalable to more than two groups.

title "Select Top N weights by group";
proc sql; 
    select a.sex, a.name, a.weight, (select count(distinct b.weight) 
            from class as b where b.weight >= a.weight and a.sex = b.sex ) as rank 
    from class as a
    where calculated rank <= 3
    order by sex, rank
;quit;

Sex	Name	Weight	rank
F	Janet	112.5	1
F	Mary	112	2
F	Carol	102.5	3
M	Philip	150	1
M	Ronald	133	2
M	Robert	128	3

Deploy a MongoDB powered Flask app in 5 minutes

2015-02-01T18:34:00.000-06:00

This is a quick tutorial to deploy a web service (a social network) by the LNMP (Linux, Nginx, MongoDB, Python) infrastructure on any IaaS cloud. The repo at Github is at https://github.com/dapangmao/minitwit-mongo-ubuntu.

Stack

The stack is built on the tools in the ecosystem of Python below.

Tool	Name	Advantage
Cloud	DigitalOcean	Cheap but fast
Server distro	Ubuntu 14.10 x64	Everything is latest
WSGI proxy	Gunicorn	Manage workers automatically
Web proxy	Nginx	Fast and easy to configure
Framework	Flask	Single file approach for MVC
Data store	MongoDB	No scheme needed and scalable
DevOps	Fabric	Agentless and Pythonic

In addition, a Supervisor running on the server provides a daemon to protect the Gunicorn-Flask process.

The MiniTwit app

The MiniTwit application is an example provided by Flask, which is a prototype of Twitter like multiple-user social network. The original application depends on SQLite. However, the data store could be modified to fit the category of NoSQL such as Google Data Store or MongoDB. A live MintiTwit demo is hosted at http://minitwit-123.appspot.com/public

Deployment

1. Install Fabric and clone the Github repo

The DevOps tool is fabric that is simply based on SSH. The fabfile.py and the staging flask files are stored on Github. We should install fabric and download the fabfile.py on the local machine before the deployment.

sudo pip install fabric 
wget https://raw.githubusercontent.com/dapangmao/minitwit-mongo-ubuntu/master/fabfile.py
fab -l

2. Enter IP from the virtual machine

A new VM from ausually emails IP address and the root password. Then we could modify the head part of the fabfile.py accordingly. There are quite a less expensive cloud providers for prototyping other than the costly Amazon EC2. For example, a minimal instance from DigitalOcean only costs five dollars a month. If SSH key has been uploaded, the password could be ignored.

env.hosts = ['YOUR IP ADDRESS'] #  Enter IP
env.user = 'root'
env.password = 'YOUR PASSWORD'  #  Enter password

3. Fire up Fabric

Now it is time to formally deploy the application. With the command below, the fabric will first install pip, git, nginx, gunicorn, supervisor and the latest MongodB, and configure them sequentially. In less than 5 minutes, a Flask and MongoDB application will be ready for use. Since DigitalOcean has its own software repository for Ubuntu, and its VMs are on SSD, the deployment is even faster, which is usually finished in one minute.

fab deploy_minitwit

Spark practice(4): malicious web attack

2015-01-08T17:45:00.001-06:00

Suppose there is a website tracking user activities to prevent robotic attack on the Internet. Please design an algorithm to identify user IDs that have more than 500 clicks within any given 10 minutes.
Sample.txt: anonymousUserID timeStamp clickCount
123    9:45am    10
234    9:46am    12
234    9:50am    20
456    9:53am    100
123    9:55am    33
456    9:56am    312
123    10:03am    110
123    10:16am    312
234    10:20am    201
456    10:23am    180
123    10:25am    393
456    10:27am    112
999    12:21pm    888

Thought

This is a typical example of stream processing. The key is to build a fixed-length window to slide through all data, count data within and return the possible malicious IDs.

Single machine solution

Two data structures are used: a queue and a hash table. The queue is scanning the data and only keeps the data within a 10-minute window. Once a new data entry is filled, the old ones out of the window are popped out. The hash table counts the data in the queue and will be updated with the changing queue. Any ID with more than 500 clicks will be added to a set.

from datetime import datetime
import time
from collections import deque

def get_minute(s, fmt = '%I:%M%p'):
    return time.mktime(datetime.strptime(s, fmt).timetuple())

def get_diff(s1, s2):
    return int(get_minute(s2) - get_minute(s1)) / 60

def find_ids(infile, duration, maxcnt):
    queue, htable, ans = deque(), {}, set()
    with open(infile, 'rt') as _infile:
        for l in _infile:
            line = l.split()
            line[2] = int(line[2])
            current_id, current_time, current_clk = line
            if current_id not in htable:
                htable[current_id] = current_clk
            else:
                htable[current_id] += current_clk
            queue.append(line)
            while queue and get_diff(queue[0][1], current_time) > duration:
                past_id, _, past_clk = queue.popleft()
                htable[past_id] -= past_clk
            if htable[current_id] > maxcnt:
                ans.add(current_id)
    return ans

if __name__ == "__main__":
    print find_ids('sample.txt', 10, 500)

Cluster solution

The newest Spark (version 1.2.0) starts to support Python streaming. However, the document is still scarce — wait to see if this problem can be done by the new API.

To be continued

Spark practice (3): clean and sort Social Security numbers

2014-12-23T21:52:00.000-06:00

Sample.txt

Requirements:
1. separate valid SSN and invalid SSN
2. count the number of valid SSN

402-94-7709 
283-90-3049 
124-01-2425 
1231232
088-57-9593 
905-60-3585 
44-82-8341
257581087
327-84-0220
402-94-7709

Thoughts

SSN indexed data is commonly seen and stored in many file systems. The trick to accelerate the speed on Spark is to build a numerical key and use the sortByKey operator. Besides, the accumulator provides a global variable existing across machines in a cluster, which is especially useful for counting data.

Single machine solution

#!/usr/bin/env python
# coding=utf-8
htable = {}
valid_cnt = 0
with open('sample.txt', 'rb') as infile, open('sample_bad.txt', 'wb') as outfile:
    for l in infile:
        l = l.strip()
        nums = l.split('-')
        key = -1
        if l.isdigit() and len(l) == 9:
            key = int(l)
        if len(nums) == 3 and map(len, nums) == [3, 2, 4]:
            key = 1000000*int(nums[0]) + 10000*int(nums[1]) + int(nums[2])
        if key == -1:
            outfile.write(l + '\n')
        else:
            if key not in htable:
                htable[key] = l
                valid_cnt += 1

with open('sample_sorted.txt', 'wb') as outfile:
    for x in sorted(htable):
        outfile.write(htable[x] + '\n')

print valid_cnt

Cluster solution

#!/usr/bin/env python
# coding=utf-8
import pyspark
sc = pyspark.SparkContext()
valid_cnt = sc.accumulator(0)

def is_validSSN(l):
    l = l.strip()
    nums = l.split('-')
    cdn1 = (l.isdigit() and len(l) == 9)
    cdn2 = (len(nums) == 3 and map(len, nums) == [3, 2, 4])
    if cdn1 or cdn2:
        return True
    return False

def set_key(l):
    global valid_cnt
    valid_cnt += 1
    l = l.strip()
    if len(l) == 9:
        return (int(l), l)
    nums = l.split('-')
    return (1000000*int(nums[0]) + 10000*int(nums[1]) + int(nums[2]), l)

rdd = sc.textFile('sample.txt')
rdd1 = rdd.filter(lambda x: not is_validSSN(x))

rdd2 = rdd.filter(is_validSSN).distinct() \
    .map(lambda x: set_key(x)) \
    .sortByKey().map(lambda x: x[1])

for x in rdd1.collect():
    print 'Invalid SSN\t', x

for x in rdd2.collect():
    print 'valid SSN\t', x

print '\nNumber of valid SSN is {}'.format(valid_cnt)

# Save RDD to file system
rdd1.saveAsTextFile('sample_bad')
rdd2.saveAsTextFile('sample_sorted')
sc.stop()

Spark practice (2): query text using SQL

2014-12-12T12:36:00.001-06:00

In a class of a few children, use SQL to find those who are male and weight over 100.

class.txt (including Name Sex Age Height Weight)

Alfred M 14 69.0 112.5 
Alice F 13 56.5 84.0 
Barbara F 13 65.3 98.0 
Carol F 14 62.8 102.5 
Henry M 14 63.5 102.5 
James M 12 57.3 83.0 
Jane F 12 59.8 84.5 
Janet F 15 62.5 112.5 
Jeffrey M 13 62.5 84.0 
John M 12 59.0 99.5 
Joyce F 11 51.3 50.5 
Judy F 14 64.3 90.0 
Louise F 12 56.3 77.0 
Mary F 15 66.5 112.0 
Philip M 16 72.0 150.0 
Robert M 12 64.8 128.0 
Ronald M 15 67.0 133.0 
Thomas M 11 57.5 85.0 
William M 15 66.5 112.0

Thoughts

The challenge is to transform unstructured data to structured data. In this question, a schema has to be applied including column name and type, so that the syntax of SQL is able to query the pure text.

Single machine solution

Straight-forward and simple if with Python’s built-in module sqlite3.

import sqlite3

conn = sqlite3.connect(':memory:')
c = conn.cursor()
c.execute("""CREATE TABLE class
             (Name text, Sex text, Age real, Height real, Weight real)""")

with open('class.txt') as infile:
    for l in infile:
        line = l.split()
        c.execute('INSERT INTO class VALUES (?,?,?,?,?)', line)
conn.commit()

for x in c.execute("SELECT * FROM class WHERE Sex = 'M' AND Weight > 100"):
    print x
conn.close()

Cluster solution

Spark SQL is built on Hive, and seamlessly queries the JSON formatted data that is semi-structured. To dump the JSON file on the file system will be the first step.

import os
import subprocess
import json
from pyspark import SparkContext
from pyspark.sql import HiveContext
sc = SparkContext()
hiveCtx = HiveContext(sc)
def trans(x):
    return {'Name': x[0], 'Sex': x[1], 'Age': int(x[2]), \
           'Height': float(x[3]), 'Weight': float(x[4])}
# Remove the output directory for JSON if it exists
if 'class-output' in os.listdir('.'):
    subprocess.call(['rm', '-rf', 'class-output'])

rdd = sc.textFile('class.txt')
rdd1 = rdd.map(lambda x: x.split()).map(lambda x: trans(x))
rdd1.map(lambda x: json.dumps(x)).saveAsTextFile('class-output')

infile = hiveCtx.jsonFile("class-output/part-00000")
infile.registerTempTable("class")

query = hiveCtx.sql("""SELECT * FROM class WHERE Sex = 'M' AND Weight > 100
      """)
for x in query.collect():
    print x

sc.stop()

In a conclusion, JSON should be considered if SQL is desired on Spark.

Spark practice (1): find the stranger that shares the most friends with me

2014-12-07T17:05:00.002-06:00

Given the friend pairs in the sample text below (each line contains two people who are friends), find the stranger that shares the most friends with me.
sample.txt
me Alice
Henry me
Henry Alice
me Jane
Alice John
Jane John
Judy Alice
me Mary
Mary Joyce
Joyce Henry
Judy me
Judy Jane
John Carol 
Carol me
Mary Henry
Louise Ronald
Ronald Thomas
William Thomas

Thoughts

The scenario is commonly seen for a social network user. Spark has three methods to query such data:

MapReduce
GraphX
Spark SQL

If I start with the simplest MapReduce approach, then I would like to use two hash tables in Python. First I scan all friend pairs and store the friends for each person in a hash table. Second I use another hash table to count my friends’ friends and pick out the strangers to me.

Single machine solution

#!/usr/bin/env python
# coding=utf-8
htable1 = {}
with open('sample.txt', 'rb') as infile:
    for l in infile:
        line = l.split()
        if line[0] not in htable1:
            htable1[line[0]] = [line[1]]
        else:
            htable1[line[0]] += [line[1]]
        if line[1] not in htable1:
            htable1[line[1]] = [line[0]]
        else:
            htable1[line[1]] += [line[0]]

lst = htable1['me']
htable2 = {}
for key, value in htable1.iteritems():
    if key in lst:
        for x in value:
            if x not in lst and x != 'me': # should only limit to strangers
                if x not in htable2:
                    htable2[x] = 1
                else:
                    htable2[x] += 1

for x in sorted(htable2, key = htable2.get, reverse = True):
    print "The stranger {} has {} common friends with me".format(x, \
        htable2[x])

The result shows that John has three common friends like I do, followed by Joyce who has two. Therefore, John will be the one who is most likely to be recommended by the social network.

Cluster solution

If the log file for the friend pairs is quite big, say, like several GB size, the single machine solution is not able to load the data into the memory and we have to seek help from a cluster.

Spark provides the pair RDD that is similar to a hash table and essentially a key-value structure. To translate the single machine solution to a cluster one, I use the operators from Spark’s Python API including map, reduceByKey, filter, union and sortByKey.

#!/usr/bin/env python
# coding=utf-8
import pyspark
sc = pyspark.SparkContext()
# Load data from hdfs
rdd = sc.textFile('hdfs://sample.txt') 
# Build the first pair RDD
rdd1 = rdd.map(lambda x: x.split()).union(rdd.map(lambda x: x.split()[::-1]))
# Bring my friend list to local
lst = rdd1.filter(lambda x: x[0] == 'me').map(lambda x: x[1]).collect()
# Build the second pari RDD
rdd2 = rdd1.filter(lambda x: x[0] in lst).map(lambda x: x[1]) \
    .filter(lambda x: x != 'me' and x not in lst) \
    .map(lambda x: (x, 1)).reduceByKey(lambda a, b: a + b) \
    .map(lambda (x, y): (y, x)).sortByKey(ascending = False)
# Save the result to hdfs
rdd2.saveAsTextFile("hdfs://sample_output")
# Bring the result to local since the sample is small
for x, y in rdd2.collect():
    print "The stranger {} has {} common friends with me".format(y, x)

sc.stop()

The result is the same. In this experiment, most time is spent on the data loading process from HDFS to the memory. The following MapReduce operations actually costs just a small fraction of overall time. In conclusion, Spark fits well on an iterative data analysis against existing RDD.

Use a vector to print Pascal's triangle in SAS

2014-12-05T13:57:00.001-06:00

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.

Minimize complexity by Spark

2014-11-28T13:59:00.001-06:00

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.

Automated testing by pytest

2014-10-20T21:18:00.000-05:00

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

One example of test-driven development in Python

2014-09-24T14:08:00.000-05:00

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.

Translate SAS's sas7bdat format to SQLite and Pandas

2014-08-10T08:30:00.000-05:00

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:

Size reduction:
SAS’s sas7bdat format is verbose. So far successfully loaded 40GB SAS data to SQLite with 85% reduction of disk usage.
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:

Pandas’ powerful Excel interface:
Write very large Excel file quickly as long as memory can hold data.
Validation of statistics
Pandas works well with statsmodels and scikit-learn. Easy to validate SAS’s outputs.

Remove tabs from SAS code files

2014-06-10T09:03:00.002-05:00

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

Steps to deploy Flask's minitwit on Google App Enginee

2014-05-26T09:58:00.000-05:00

Flask is a light-weight web framework for Python, which is well documented and clearly written. Its Github depository provides a few examples, which includes minitwit. The minittwit website enjoys a few basic features of social network such as following, login/logout. The demo site on GAE is http://minitwit-123.appspot.com. The Github repo is https://github.com/dapangmao/minitwit.

Google App Engine or GAE is a major public clouder service besides Amazon EC2. Among the four languages(Java/Python/Go/PHP) it supports, GAE is friendly to Python users, possibly because Guido van Rossum worked there and personally created Python datastore interface. As for me, it is a good choice for a Flask app.

Step1: download GAE SDK and GAE Flask skeleton

GAE’s Python SDK tests the staging app and eventuall pushes the app to the cloud.

A Flask skeleton can be dowloaded from Google Developer Console. It contains three files:

app.yaml: specify the entrance of run-time
appengine_config.py: add the external libraries such as Flask to system path
main.py: the root Python program

Step2: schema design

The dabase used for the original minitwit is SQLite. The schema consists of three tables: user, follower and message, which makes a normalized database together. GAE has two Datastore APIs: DB and NDB. Since neither of them supports joining (in this case one-to-many joining for user to follower), I move the follwer table as an nested text propery into the user table, which eliminatse the need for joining.

As the result, the main.py has two data models: User and Message. They will create and maintain two kinds (or we call them as tables) with the same names in Datastore.

class User(ndb.Model):
  username = ndb.StringProperty(required=True)
  email = ndb.StringProperty(required=True)
  pw_hash = ndb.StringProperty(required=True)
  following = ndb.IntegerProperty(repeated=True)
  start_date = ndb.DateTimeProperty(auto_now_add=True)

class Message(ndb.Model):
  author = ndb.IntegerProperty(required=True)
  text = ndb.TextProperty(required=True)
  pub_date = ndb.DateTimeProperty(auto_now_add=True)
  email = ndb.StringProperty(required=True)
  username = ndb.StringProperty(required=True)

Step3: replace SQL statements

The next step is to replace SQL operations in each of the routing functions with NDB’s methods. NDB’s two fundamental methods are get() that retrieves data from Datastore as a list, and put() that pushes list to Datastore as a row. In short, data is created and manipulated as individual object.

For example, if a follower needs to add to a user, I first retrieve the user by its ID that returns a list like [username, email, pw_hash, following, start_date], where following itself is a list. Then I insert the new follower into the following element and save it back again.

u = User.get_by_id(cid)
if u.following is None:
  u.following = [whom_id]
  u.put()
else:
  u.following.append(whom_id)
  u.put()

People with experience in ORM such as SQLAlchemy will be comfortable to implement the changes.

Setp4: testing and deployment

Without the schema file, now the minitwit is a real single file web app. It’s time to use GAE SDK to test it locally, or eventually push it to the cloud. On GAE, We can check any error or warning through the Logs tab to find bugs, or view the raw data through the Datastore Viewer tab.

In conclusion, GAE has a few advantages and disadvantages to work with Flask as a web app.

Pro:
- It allows up to 25 free apps (great for exercises)
- Use of database is free
- Automatical memoryCached for high IO
Con:
- Database is No-SQL, which makes data hard to port
- More expensive for production than EC2

Use recursion and gradient ascent to solve logistic regression in Python

2014-05-21T16:22:00.003-05:00

In his book Machine Learning in Action, Peter Harrington provides a solution for parameter estimation of logistic regression . I use pandas and ggplot to realize a recursive alternative. Comparing with the iterative method, the recursion costs more space but may bring the improvement of performance.

# -*- coding: utf-8 -*-
"""
Use recursion and gradient ascent to solve logistic regression in Python
"""

import pandas as pd
from ggplot import *

def sigmoid(inX):
    return 1.0/(1+exp(-inX))

def grad_ascent(dataMatrix, labelMat, cycle):
    """
    A function to use gradient ascent to calculate the coefficients
    """
    if isinstance(cycle, int) == False or cycle < 0:
        raise ValueError("Must be a valid value for the number of iterations")
    m, n = shape(dataMatrix)
    alpha = 0.001
    if cycle == 0:
        return ones((n, 1))
    else:
        weights = grad_ascent(dataMatrix, labelMat, cycle-1)
        h = sigmoid(dataMatrix * weights)
        errors = (labelMat - h)
        return weights + alpha * dataMatrix.transpose()* errors

def plot(vector):
    """
    A funtion to use ggplot to visualize the result
    """
    x = arange(-3, 3, 0.1)
    y = (-vector[0]-vector[1]*x) / vector[2]
    new = pd.DataFrame()
    new['x'] = x
    new['y'] = array(y).flatten()
    infile.classlab = infile.classlab.astype(str)
    p = ggplot(aes(x='x', y='y', colour='classlab'), data=infile) + geom_point()
    return p + geom_line

# Use pandas to manipulate data
if __name__ == '__main__':
    infile = pd.read_csv("https://raw.githubusercontent.com/pbharrin/machinelearninginaction/master/Ch05/testSet.txt", sep='\t', header=None, names=['x', 'y', 'classlab'])
    infile['one'] = 1
    mat1 = mat(infile[['one', 'x', 'y']])
    mat2 = mat(infile['classlab']).transpose()
    result1 = grad_ascent(mat1, mat2, 500)
    print plot(result1)

Count large chunk of data in Python

2014-04-30T16:09:00.001-05:00

The line-by-line feature in Python allows it to count hard disk-bound data. The most frequently used data structures in Python are list and dictionary. Many cases the dictionary has advantages since it is a basically a hash table that many realizes O(1) operations.

However, for the tasks of counting values, the two options make no much difference and we can choose any of them for convenience. I listed two examples below.

Use a dictionary as a counter

There is a question to count the strings in Excel.

Count the unique values in one column in EXCEL 2010. The worksheet has 1 million rows and 10 columns.
or numbers.
For example,
A5389579_10
A1543848_6
A5389579_8
Need to cut off the part after (including) underscore such as from A5389579_10 to A5389579

Commonly Excel on a desktop can’t handle this size of data, while Python would easily handle the job.

# Load the Excel file by the xlrd package
import xlrd
book = xlrd.open_workbook("test.xlsx")
sh = book.sheet_by_index(0)
print sh.name, sh.nrows, sh.ncols
print "Cell D30 is", sh.cell_value(rowx=29, colx=3)

# Count the unique values in a dictionary
c = {} 
for rx in range(sh.nrows):
    word = str(sh.row(rx)[1].value)[:-3]
    try:
        c[word] += 1 
    except:
        c[word] = 1

print c

Use a list as a counter

There is a question to count emails.

A 3-column data set includes sender, receiver and timestamp. How to calculate the time between the sender sends the email
and the receiver sends the reply email?

The challenge is to scale up the small sample data to larger size. The solution I have has the complexity of O(nlogn), which is only limited by the sorting step.

raw_data = """
    SENDER|RECEIVER|TIMESTAMP
    A B 56
    A A 7
    A C 5
    C D 9
    B B 12
    B A 8
    F G 12
    B A 18
    G F 2
    A B 20
    """

# Transform the raw data to a nested list
data = raw_data.split()
data.pop(0) # Remove the Head
data = zip(data[0::3], data[1::3], map(lambda x: int(x), data[2::3]))

# Sort the nested list by the timestamp 
from operator import itemgetter
data.sort(key=itemgetter(2))
for r in data:
    print r

# Count the time difference in a list
c = []
while len(data) != 1:
    y = data.pop(0)
    for x in data:
        if x[0] == y[1] and x[1] == y[0]:
            diff = x[2] - y[2] 
            print y, x, '---->', diff
            c.append(diff)
            break # Only find the quickest time to respond
print c

P.S.

I come up with the O(n) solution below, which utilizes two hash tables to decrease the complexity.

__author__ = 'dapangmao'

def find_duration(data):
    # Construct two hash tables
    h1 = {}
    h2 = {}
    # Find the starting time for each ID pair
    for x in data:
        if x[0] != x[1]:
            key = x[0] + x[1]
            try:
                h1[key] = x[2]
            except:
                h1[key] = min(h1[key], x[2])
    # Find the minimum duration for each ID pair
    for x in data:
        key = x[1] + x[0]
        if h1.has_key(key):
            duration = x[2] - h1[key]
            try:
                h2[key] = duration
            except:
                h2[key] = min(h2[key], duration)
    return h2

if __name__ == "__main__":
    raw_data = """
        SENDER|RECEIVER|TIMESTAMP
        A B 56
        A A 7
        A C 5
        C D 9
        B B 12
        B A 8
        F G 12
        B A 18
        G F 2
        A B 20
        """

    # Transform the raw data to a nested list
    data = raw_data.split()
    data.pop(0) # Remove the Head
    data = zip(data[0::3], data[1::3], map(lambda x: int(x), data[2::3]))
    # Verify the result
    print find_duration(data)

10 popular Linux commands for Hadoop

2014-04-06T08:41:00.002-05:00

The Hadoop system has its unique shell language, which is called FS. Comparing with the common Bash shell within the Linux ecosystem, the FS shell has much fewer commands. To deal with the humongous size of data distributively stored at the Hadoop nodes, in my practice, I have 10 popular Linux command to facilitate my daily work.

1. sort
A good conduct of running Hadoop is to always test the map/reduce programs at the local machine before releasing the time-consuming map/reduce codes to the cluster environment. The sort command simulates the sort and shuffle step necessary for the map/redcue process. For example, I can run the piped commands below to verify whether the Python codes have any bugs.
./mapper.py | sort | ./reducer.py

2. tail
Interestingly, the FS shell at Hadoop only supports the tail command instead of the head command. Then I can only grab the bottom lines of the data stored at Hadoop.
hadoop fs -tail 5 data/web.log.9

3. sed
Sine the FS shell doesn’t provide the head command, the alternative solution is to use the sed command that actually has more flexible options.
hadoop fs -cat data/web.log.9 | sed '1,+5!d'

4. stat
The stat command allows me to know the time when the file has been touched.
hadoop fs -stat data/web.log.9

5. awk
The commands that the FS shell supports usually have very few options. For example the du command under the FS shell does not support -sh option to aggregate the disk usage of the sub-directories. In this case, I have to look for help from the awk command to satisfy my need.
hadoop fs -du data | awk '{sum+=$1} END {print sum}'

6. wc
One of the most important things to understand a file located at the Hadoop is to find the number of its total lines.
hadoop fs -cat data/web.log.9 | wc -l

7. cut
The cut command is convenient to select the specified columns at the file. For example, I am able to count the lines for each of the unique groups from the column between the position of #5 and #14.
hadoop fs -cat data/web.log.9 | cut -c 5-14 | uniq -c

8. getmerge
The great thing for the getmerge command is that I am able to fetch all the result after map/reduce to the local file system as a single file.
hadoop fs -getmerge result result_merged.txt

9. grep
I can start a mapper-only job only with the grep command form the Bash shell to search the lines which contain the key words I am interested in. And this is a map-only task.

hadoop jar $STREAMING -D mapred.reduce.tasks=0 -input data -output result -mapper "bash -c 'grep -e Texas'"

10. at and crontab
The at and crontab commnands are my favorite to schedule a job at Hadoop. For example, I would like to use the order below to clean the map/reduce results at midnight.

at 0212
at > hadoop fs -rmr result

Software & Service

Algorithmic pricing with Redis

The background

How an algorithmic pricing engine works

1. The distributed crawler

2. The algorithm core

3. The price adjustor

The role of Redis

1. Predicative pricing

2. Cross-data-center hedging

Good math, bad engineering

1. The quick solution

2. The third step that is hard

3. Finally the easy solution

Conclusion

Use Elasticsearch and Kibana for large BI system

1. Batch Worker

1.1 Spring Data Batch v.s. Python

Spring Data Batch

Python

1.2 The algorithms

1.3 The scheduler

1.4 High availability with zero downtime

2. Elasticsearch Cluster

2.1 Three kinds of nodes

2.2 Memory is cheap but heap is expensive

The heap usage and the old GC duration

Thread pools

The segments’ number/size

3. Kibana as BI dashboard

3.1 Plugins

3.2 Instant Aggregation

3.3 Proxy

Conclusion

Use Slack bot to monitor the server

Deploy edx spark environment to DigitalOcean

Transform SAS files to Parquet through Spark

1. Test DataFrame in Python and Spark

2. Automate the transformation

3. Conclusion

Two alternative ways to query large dataset in SAS

Example 1: Find the median

SQL Query

PROC UNIVARIATE

Read Directly

Example 2: Find top K by groups

SQL Query

Code Generator

Read Directly

Example 3: Find Moving Window Maxium

Code Generator

Read Directly

Conclusion

While data is expanding, we should more and more consider three things -Time complexity: we don't want run data for weeks.Space complexity: we don't want the memory overflow.Clean codes: the colleagues should easily Reading and modify the codes.

saslib: a simple Python tool to lookup SAS metadata

Installation

Usage

Deploy a minimal Spark cluster

Requirements

Choice of public cloud:

Add Spark to DigitalOcean cluster

Solve the Top N questions in SAS/SQL

Prepare the data

1. Select highest value

2. Select second highest value

3. Select Nth highest value

4. Select highest values by group

5. Rank all values

6. Select top N values by group

Deploy a MongoDB powered Flask app in 5 minutes

Stack

The MiniTwit app

Deployment

1. Install Fabric and clone the Github repo

2. Enter IP from the virtual machine

3. Fire up Fabric

Spark practice(4): malicious web attack

Thought

Single machine solution

Cluster solution

While data is expanding, we should more and more consider three things -
Time complexity: we don't want run data for weeks.
Space complexity: we don't want the memory overflow.
Clean codes: the colleagues should easily Reading and modify the codes.