Wednesday, July 20, 2016

google interview: find number in multi disk list better than O(logn)


// http://www.impactinterview.com/2009/10/140-google-interview-questions/#software_engineer
//
// Find or determine non existence of a number in a sorted list
// of N numbers where the numbers range over M, M>> N and N large
// enough to span multiple disks. Algorithm to beat O(log n) bonus
// points for constant time algorithm.

// the const time algo is called a "bloom filter" if you work with large
// databases you (hopefully) learn this..  what a bloom filter does on insert
// is hash your key in 3 ways and then light up flags in 3 indexes..  on lookup
// it hashs again and checks the 3 places..  if the any one the flags is off
// then you have a item that doesnt exist..  otherwise commence your normal
// logN search

#include <iostream>
#include <vector>
#include <stdint.h>
#include <memory>
#include <cstring>
#include <cstdlib>
#include <chrono>

#define BLOOMFILTER

class PagedVector
{
public:
    typedef uint64_t                Key;         // ideally this would be a template param
    typedef std::vector<Key>        Keys;
    typedef std::shared_ptr<Keys>   KeysHandle;
    typedef std::vector<KeysHandle> PagedKeys;

    const uint32_t filterSize = 256;             // ideally this would be tunable to reduce misses/mem use
    const uint32_t filterMask = filterSize - 1;

private:
    std::size_t threashold_;

    PagedKeys pages_;

    std::vector<bool> presence_[3]; // the bloom filter

    // the divide part of a divide and conqurer search
    template <typename Iterator,
              typename Compare,
              typename Splitter>
    struct Divider
    {
        Compare&  comparer_;
        Splitter& splitter_;

        Divider(Compare&  comparer,
                Splitter& splitter) :
            comparer_(comparer),
            splitter_(splitter)
        {}

        void operator()(Iterator& a,
                        Iterator& b)
        {
            // should never happen unless we are empty..
            if (b==a) return;

            Iterator n = splitter_(a,b);

            if (comparer_(n))
            {
                b = n;
                return;
            }
            a = n + 1;
        }
    };

    // All divide and conqurer searches do 1 thing..  they split the data in
    // some way and repeat until the desired location is found..  In this
    // system we have mutliple uses of the search.. they are used to search and find insertion points

    template <typename Iterator,
              typename Compare,
              typename Splitter>
    Iterator search(Iterator begin,
                    Iterator end,
                    Compare  comparer,
                    Splitter splitter)
    {
        // point at the match or the point to insert in if it doesnt match
        Divider<Iterator, Compare, Splitter> Divider(comparer, splitter);

        // narrows the sub range in which to keep searching when the range is 1 item we are done..
        while (1)
        {
            Divider(begin,end);
            if (end == begin)
                return begin;
        }
    }

    struct FindInPage
    {
        Key target_;

        FindInPage(Key target) :
            target_(target)
        {}

        bool operator()(const PagedKeys::iterator& a)
        {
            // ok slight strangeness here.. pages are a range of numers
            //  t < [b,e] then true
            //  t > [b,e] then false
            // but
            // t > b, t < e its equal..  so it has to be absolete less than!

            const KeysHandle& pageHandle = *a;
            if (pageHandle->size() > 0)
            {
                return target_ <= (*pageHandle)[pageHandle->size()-1];

            }
            return true; // empty... should be the last if it exists..
        }
    };

    struct SplitPage
    {
        Key target_;

        SplitPage()
        {}

        PagedKeys::iterator operator()(const PagedKeys::iterator& a,
                                       const PagedKeys::iterator& b)
        {
            std::size_t offset = (b-a)/2;
            return a + offset;
        }
    };

    struct FindInVector
    {
        Key target_;

        FindInVector(Key target) :
            target_(target)
        {}

        bool operator()(const Keys::iterator& a)
        {
            return target_ <= (*a);
        }
    };

    struct SplitVector
    {
        Key target_;

        SplitVector(Key target) :
            target_(target)
        {}

        Keys::iterator operator()(const Keys::iterator& a,
                                  const Keys::iterator& b)
        {
            // std::size_t offset = (b-a)/2;
            // return a + offset;

            if ((b-a) <= 2)
            {
                std::size_t offset = (b-a)/2;
                return a + offset;
            }

            Key high = *(b-1);
            Key low  = *a;

            // since we have the values...
            if (target_ <= low)  return a;
            if (target_ >= high) return b-1;

            double splitPercent;
            splitPercent = static_cast<double>(target_ - low) / static_cast<double>(high-low);
            splitPercent = (splitPercent < 0.0) ? 0.0 : splitPercent;

            std::size_t range  = (b-a);
            std::size_t offset = range * splitPercent;

            return  a + offset;
        }
    };

    // retun a fax iterator to the page and page offset where the item is or should be
    std::pair<PagedKeys::iterator, Keys::iterator> find(Key target)
    {
        std::pair<PagedKeys::iterator, Keys::iterator> iter;

        iter.first = search(pages_.begin(),
                            pages_.end(),
                            FindInPage(target),
                            SplitPage());

        if (iter.first == pages_.end())
            iter.first = pages_.end() - 1;

        KeysHandle& page = *(iter.first);

        // find the location to add into the list
        iter.second = search(page->begin(),
                             page->end(),
                             FindInVector(target),
                             SplitVector(target));
        return iter;
    }

public:
    PagedVector(std::size_t threashold) :
        threashold_(threashold),
        pages_(),
        presence_()
    {
        for (int i = 0; i < 3; ++i)
            presence_[i].resize(filterSize);

        // make it a bit more easy on the boarder conditions
        KeysHandle newPage(new Keys);
        pages_.push_back(newPage);
    }

    // add the target to the paged data struture
    void insert(Key target)
    {
#ifdef BLOOMFILTER
        // not the best way but assume we are doing a item by item insert
        uint64_t hash = std::hash<Key>()(target);

        // light the bloomfilters presence bits according to the hash
        presence_[0][ hash        & filterMask] = true;
        presence_[1][(hash >> 8)  & filterMask] = true;
        presence_[2][(hash >> 16) & filterMask] = true;
#endif

        // find the relevent part of the list
        std::pair<PagedKeys::iterator, Keys::iterator> iter = find(target);

        KeysHandle& page = *(iter.first);

        page->insert(iter.second, target);

        // then threshold check.
        if (page->size() > threashold_)
        {
            // ok over the threshold.. so divide the list in half
            KeysHandle newPage(new Keys);
            uint32_t spliceAt = page->size() / 2;
            uint32_t sizeAfter = page->size() - spliceAt;
            newPage->resize(sizeAfter);
            std::memcpy(&((*newPage)[0]), &((*page)[spliceAt]), sizeAfter * sizeof(Key));
            page->resize(spliceAt);

            pages_.insert(iter.first + 1, newPage);
        }
    }

    bool exists(Key target)
    {
#ifdef BLOOMFILTER
        uint64_t hash = std::hash<Key>()(target);

        // a quick O(1) pre check before the more costly real check
        if (not (presence_[0][ hash        & filterMask] and
                 presence_[1][(hash >> 8)  & filterMask] and
                 presence_[2][(hash >> 16) & filterMask]))
            return false;
#endif

        // ok technically the question was about this.. so no std:: here..
        // basically have to implement a better then log(N) search
        std::pair<PagedKeys::iterator, Keys::iterator> iter = find(target);

        return *(iter.second) == target;
    }

    const PagedKeys& pages() const
    {
        return pages_;
    }

    void print(std::ostream& os) const
    {
        for ( KeysHandle page : pages_)
        {
            os << "\n ########################################## \n";
            for ( Key item : *page)
            {
                os << " " << item;
            }
        }
    }
};

std::ostream& operator<<(std::ostream& os, const PagedVector& list)
{
    list.print(os);
    return os;
}

void validate_sorted(PagedVector& theList)
{
    bool failed = false;
    bool first = true;
    PagedVector::Key prev;

    for ( PagedVector::KeysHandle page : theList.pages() )
    {
        for ( PagedVector::Key item : *page)
        {
            if (first)
            {
                first = false;
                prev = item;
            }
            else
            {
                if (prev > item)
                    failed = true;
                prev = item;
            }
        }
    }
    std::cout << (failed ? "FAILED" : "passed")
              << "\n";
}

void test_random()
{
    PagedVector theList(5000);

    auto istart = std::chrono::system_clock::now();
    for (int i = 0; i < 1000000;  ++i)
    {
        PagedVector::Key item = (std::rand() % 100000)*2;

        // std::cout << "adding:" << item << "\n";
        theList.insert(item);
    }
    auto iend = std::chrono::system_clock::now();

    std::cout << "inserted time:"
              << std::chrono::duration_cast<std::chrono::milliseconds>(iend - istart).count()
              << "mSec \n";
    // std::cout << theList << "\n";

    int count = 0;
    int exist = 0;

    auto cstart = std::chrono::system_clock::now();
    for (int i = 0; i < 1000000;  ++i)
    {
        PagedVector::Key item = std::rand() % 200000;
        if (theList.exists(item))
            ++exist;
        ++count;

        // std::cout << " " << item << " "
        //           << (theList.exists(item) ? "exists" : "missing")
        //           << "\n";
    }
    auto cend = std::chrono::system_clock::now();

    std::cout << "checking time:"
              << std::chrono::duration_cast<std::chrono::milliseconds>(cend - cstart).count()
              << "mSec\n\n";

    // expecting a 50% hit rate
    std::cout << "checked:"
              << count
              << " exist:" << exist
              << "\n";

    validate_sorted(theList);
}

int main()
{
    test_random();
}

Output with a bloomfilter is:
inserted time:2613mSec 
checking time:735mSec

checked:1000000 exist:499965
passed

The output without the filter is
inserted time:2584mSec 
checking time:1292mSec

checked:1000000 exist:499965
passed

Google interview: implement sort X and discuss its space/time complexity

The problem with sort algorithms is there are many tricks to making them quicker. These tricks also make it more complex to do in an interview.

I have failed a google interview in the past because of this, at that time my interviewer asked me to code a *stable* quick sort on the white board.. The moment it occurred to just how tricky the code will be to do on a white board i simply froze up and promptly forgot everything.. I went down the rabbit hole so to speak, worst interview ever..

So i have tried to keep the algorithms as simple as possible. My tactic here is to break them into an to easy to remember set of sub operations. They are not optimal versions.

#include <iostream>
#include <vector>
#include <functional>

// ***********************************************
// ************** ITERATIVE SORTS ****************
// ***********************************************

void insertionSort(std::vector<int>& list)
{
    // ok insertion sort takes current object a inserts it in the correct place
    // inplace (excluding swap memory)
    // stable ordering
    //  -- because the item only moves if less than.. equals stop the inner loop
    // average case: O(n^2)
    // -- because it has to do some percent of the inner loop lets say 0.5..
    //     hence 0.5*N*N is still O(N^2)
    // best case: O(n)
    //  -- sorted and nothing moves
    // worst case: O(n^2)
    //  -- reversed and everything moves

    for (int i = 1 ; i < list.size(); ++i)
    {
        // now swap item down until its in place
        int j = i;
        while (j > 0 and list[j] < list[j-1])
        {
            int temp = list[j];
            list[j] = list[j-1];
            list[j-1] = temp;
            --j;
        }
    }
}

void selectionSort(std::vector<int>& list)
{
    // selection sort finds the correct object and for the current place
    //
    // inplace
    // stable ordering
    //  -- because the first of of equal items is selected
    // average case: O(n^2)
    // -- same reason as best
    // best case: O(n^2)
    //  -- even sorted u have to check all items left over
    // worst case: O(n^2)
    //  -- everything moves

    for (int i = 0 ; i < list.size(); ++i)
    {
        int selection = i;
        for (int j = i; j < list.size(); ++j)
        {
            if (list[selection] > list[j])
                selection = j;
        }

        if (selection != i)
        {
            int temp = list[i];
            list[i] = list[selection];
            list[selection] = temp;
        }
    }
}

// ***********************************************
// *********** DIVIDE AND CONQUER SORTS **********
// ***********************************************

void swap(std::vector<int>::iterator a,
          std::vector<int>::iterator b)
{
    int temp = *a;
    *a = *b;
    *b = temp;
}

void quickSort(std::vector<int>::iterator begin,
               std::vector<int>::iterator end)
{
    // ok divide and conquer by partition on a pivot
    // inplace
    // unstable
    //  -- because pivot moves out of order and then we swap left end woith right end
    // average case O(nlogn)
    // -- we linear scan the array and swap, with hopefully leaves a 50/50
    //    split and recurse for log(n) deep
    // best case O(nlogn)
    // -- still have to scan for the swaps..
    // worst case O(n^2)
    // -- if the pivot value sucks it will split all 2 one side

    // the range is [begin, end)..
    if ((end - begin) < 2) return;

    // choose pivot.. lets say midway...
    std::vector<int>::iterator p = begin + (end - begin)/2;
    int pivot = *(p);

    std::vector<int>::iterator head = begin;
    std::vector<int>::iterator tail = end - 1;

    // swap partition out of the way - If you dont remove a value then *when*
    // partitioning fails you end up stuck in an infinite loop
    swap(p, tail);
    --tail;

    // partition
    while (head != tail)
    {
        if (*head <= pivot)
            // head is in right place move head up
            ++head;
        else
        {
            // swivel head to tail and move tail down
            swap(head, tail);
            --tail;
        }
    }

    // check last item againest the pivot
    if (*head <= pivot) ++head;

    // swap the pivot back to the middle.
    swap(head, end-1);

    // note carefully this is infinite loop avoidance...  the head now points
    // at the pivot so exclude that from recursion
    quickSort(begin, head);   // Note  [begin, head)
    quickSort(head+1, end);   // Note         (head, end)
}

std::vector<int>::iterator rotate(std::vector<int>::iterator begin,
                                  std::vector<int>::iterator mid,
                                  std::vector<int>::iterator end)
{
    std::vector<int>::iterator p1 = begin;
    std::vector<int>::iterator p2 = mid;

    while (1)
    {
        swap(p1,p2);

        ++p1;
        ++p2;

        if (p1 >= mid and p2 == end) break;
        if (p2 == end) p2 = mid;
    }

    return begin + (end - mid);
}

std::vector<int>::iterator stablePartition(std::vector<int>::iterator cursor,
                                           std::vector<int>::iterator end,
                                           std::function<bool (int)>  op)
{
    // Note "op" will be true if its in the left side

    // move through the data until we find the end of the *first* left
    // section
    while (cursor < end and op(*cursor)) ++cursor;

    // record the start of the right section
    std::vector<int>::iterator rightStart = cursor;

    // find the end of the right section and the start of out of order
    // left section
    while (cursor < end and not op(*cursor)) ++cursor;

    while (cursor < end)
    {
        // if we are not at the end of the entire range theb the opFlase secton
        // has to be *out of place*

        // So record the end of the out of place right section
        std::vector<int>::iterator rightEnd = cursor;

        // find end of the out of place left section
        while (cursor < end and op(*cursor)) ++cursor;

        // then rotate the out of order left section in front of the right
        // section, and update the rightStart to the new divide point
        rightStart = rotate(rightStart, rightEnd, cursor);

        // then move the cursor to the new end of the right section, starting
        // from the end of the prior searched region
        while (cursor < end and not op(*cursor)) ++cursor;
    }

    // return the divide point between the left and rightsections
    return rightStart;
}


void stableQuickSort(std::vector<int>::iterator begin,
                     std::vector<int>::iterator end)
{
    // divide and conqure using *stable* partitioning of data
    // the trick with stable quick sort is that you rotate the sequances of
    // data to instead of swap the start and end to keep the data in order.
    //
    // It is a real pain..  and tricky to get right...  and would be very
    // difficult to do in an interview because partition failure avoidance is
    // much harder..  u need to break the sort range into 3 parts so u can leave
    // out the middle an avoid the infinite loop when partitioning fails
    //
    // AND NOTE BEST i have been asked this in a interview with google before.. (and i failed)
    //
    // inplace
    // stable
    //  -- the use of rotates moves the entire sequence keeping order
    //  -- however care must be taken so that pivot is also stable.  The way I
    //     see it there is only 1 real option you have to rotate and form 3
    //     sections the lessThan, equalTo and greaterThan sections. You can
    //     then avoid infinite loops on partition failures by excluding the
    //     equalTo section from the recurse
    // average case O(log(n) log(n) n)
    // -- http://www.codeproject.com/Articles/26048/Fastest-In-Place-Stable-Sort
    // worst case O(n^3) (i think)
    // -- if the pivoit value sucks it will split all 2 one side and do N
    //    recursions, if that partition is N and a worst will repeatedly call
    //    rotate which is worst case N then we have N*N*N or N^3

    // a list of less than 2 items is sorted by default
    if ((end - begin) < 2) return;

    // choose pivot.. mid point...
    std::vector<int>::iterator p = begin + (end - begin)/2;

    // record the pivot value..
    int pivot = *(p);

    // form the lessThan and greaterThanEqualTo sections
    std::vector<int>::iterator greaterThanEqualToStart
        = stablePartition(begin,
                          end,
                          [pivot](int v) { return v < pivot; });

    // form the equalTo and greaterThan sections
    std::vector<int>::iterator greaterThanStart
        = stablePartition(greaterThanEqualToStart,
                          end,
                          [pivot](int v) { return not (pivot < v); });

    // now recurse on the unhanded lessThan and greaterThan partitions *only*
    stableQuickSort(begin,            greaterThanEqualToStart);
    stableQuickSort(greaterThanStart, end);
}

void mergeSort(std::vector<int>::iterator begin,
               std::vector<int>::iterator end,
               std::vector<int>& tmp)
{
    // divide and conqure sort using *merging* of data
    //
    // stable
    //  -- because the left side merges in before right if the values are the same
    // average case O(nlogn)
    // -- we always divide in 2 and recurse for log(n) deep then do 2n in merge actions
    // best case O(nlogn)
    // -- same as average..
    // worst case O(nlogn)
    // -- same as best and average
    // worst space complexity O(n)
    // -- in theory u can grow tmp as u need to min the space..but practically
    //    this may result in issues because resizing a large vector temporally
    //    2x the space as it has to copy the old buffer to the new one, you can
    //    avoid this by allocating new pages.. but most of time that not what i
    //    see in practice

    if ((end - begin) < 2) return;

    // split
    std::vector<int>::iterator n = begin + (end - begin)/2;

    // recurse
    mergeSort(begin,n,tmp);
    mergeSort(n,end,tmp);

    //merge
    std::vector<int>::iterator left   = begin;
    std::vector<int>::iterator right  = n;
    std::vector<int>::iterator wcursor = tmp.begin();

    // stop the moment left runs out there is no need to copy all of right to tmp
    while ( left < n )
    {
        if      (right == end)   *(wcursor++) = *(left++ );
        else if (*right < *left) *(wcursor++) = *(right++);
        else                     *(wcursor++) = *(left++ );
    }

    // copy the section we merged into tmp back out (note this may not be the full length)
    std::vector<int>::iterator ccursor = tmp.begin();
    while(ccursor != wcursor)
    {
        *(begin++) = *(ccursor++);
    }
}


void inplaceMergeSort(std::vector<int>::iterator begin,
                      std::vector<int>::iterator end)
{
    // divide and conquer sort using *inplace* merging of data
    // the trick with inplace merge sort is that you rotate the sequences of
    // data to instead of copy then into the new buffer
    //
    // inplace
    //  -- the use of rotates swaps data multiple times to keep the in the same locations
    // stable
    // average case O(log(n) log(n) n)
    // -- http://www.codeproject.com/Articles/26048/Fastest-In-Place-Stable-Sort

    if ((end - begin) < 2) return;

    // split
    std::vector<int>::iterator n = begin + (end - begin)/2;

    // recurse
    inplaceMergeSort(begin,n);
    inplaceMergeSort(n,end);

    //merge
    std::vector<int>::iterator left  = begin;
    std::vector<int>::iterator right = n;
    while ( left < right)
    {
        // find rotate start ..
        while (not (*right < *left) and left < right) ++left;
        if (not (left < right)) break;

        // note rotate mid
        std::vector<int>::iterator mid = right;

        // find rotate end
        while (*right < *left and right < end) ++right;
        if (mid == right) break;

        // swap the sections of the array
        left = rotate(left,mid,right);
    }
}

// ***********************************************
// ****************** HEAP SORT ******************
// ***********************************************

std::vector<int>::iterator parent(std::vector<int>::iterator begin,
                                  std::vector<int>::iterator i)
{
    return begin + (i - begin - 1)/2;
}

std::vector<int>::iterator child(std::vector<int>::iterator begin,
                                 std::vector<int>::iterator i)
{
    return begin + 2 * (i - begin) + 1;
}

void upShift(std::vector<int>::iterator begin,
             std::vector<int>::iterator i)
{
    while(begin != i)
    {
        std::vector<int>::iterator p = parent(begin, i);
        if (*p > *i) return;
        swap(p,i);

        i = p;
    }
}

void makeHeap(std::vector<int>::iterator begin,
              std::vector<int>::iterator end)
{
    // version 1: top down O(nlogn)..
    if((end - begin) < 2) return;

    for(std::vector<int>::iterator i = begin + 1;
        i != end;
        ++i)
    {
        upShift(begin, i);
    }
}

void reheap(std::vector<int>::iterator begin,
            std::vector<int>::iterator end)
{
    std::vector<int>::iterator i = begin;
    while(1)
    {
        std::vector<int>::iterator kid1 = child(begin, i);
        std::vector<int>::iterator kid2 = kid1 + 1;

        if (not (kid1 < end)) return;

        std::vector<int>::iterator kidmax
            = ((kid2 < end) and (*kid1 < *kid2)) ? kid2 : kid1;

        if (not (*kidmax > *i)) return;

        swap(kidmax, i);
        i = kidmax;
    }
}

void heapSort(std::vector<int>::iterator begin,
              std::vector<int>::iterator end)
{
    // Heap sort, is the same as selection. Ie:  Find the correct item for the
    // current location but there is a critical difference, heap sort prepares
    // the unsorted data as a tree that forms a priority index with the
    // larger numbers the parents of children nodes
    //
    // inplace
    // unstable
    //   -- the heap operators dont work to keep stable order
    // best case: O(nlogn)  ...
    //  -- building a heap can be done as O(n) but this implementation doesn't use it

    if ((end - begin) < 2) return;

    // make the heap
    makeHeap(begin,end);

    //then starting at the back
    std::vector<int>::iterator i = end - 1;
    while (i != begin) // dont bother with the last it has to be the largest
    {
        // move the largest to the end
        swap(i, begin);

        // and repair the heap
        reheap(begin,i);
        --i;
    }
}

// ***********************************************
// ************** OVERLOAD WRAPPERS **************
// ***********************************************

// should be a real overload but std::functional cant seem to match the
// overload requested

void quickSortWrap(std::vector<int>& list)
{
    quickSort(list.begin(), list.end());
}

void stableQuickSortWrap(std::vector<int>& list)
{
    stableQuickSort(list.begin(), list.end());
}

void mergeSortWrap(std::vector<int>& list)
{
    std::vector<int> tmp;
    tmp.resize(list.size());

    mergeSort(list.begin(), list.end(), tmp);
}

void inplaceMergeSortWrap(std::vector<int>& list)
{
    inplaceMergeSort(list.begin(), list.end());
}

void heapSortWrap(std::vector<int>& list)
{
    heapSort(list.begin(), list.end());
}

// ***********************************************
// ******************** TESTING ******************
// ***********************************************

bool check(std::vector<int>& list)
{
    bool failed = false;
    bool first = true;
    int prev;

    for ( int i : list)
    {
        if (first)
        {
            prev = i;
            first = false;
        }
        else if (i < prev)
        {
            std::cout << "*";
            failed = true;
        }
        prev = i;
        std::cout << i << " ";
    }
    std::cout << " :" << (failed ? "FAILED" : "passed") << "\n";
}

int testHeap(std::vector<int> data)
{
    // check indexing
    if (data.size() > 1)
    {
        // some hard checks
        // index 1 -> parent 0
        if (data.size() > 1)
            std::cout << " parent1: "
                      << ((parent(data.begin(), data.begin()+1) == data.begin()) ?
                          "passed" : "FAILED")
                      << "\n";

        // index 2 -> parent 0
        if (data.size() > 2)
            std::cout << " parent2: "
                      << ((parent(data.begin(), data.begin()+2) == data.begin()) ?
                          "passed" : "FAILED")
                      << "\n";

        // index 3 -> parent 1
        if (data.size() > 3)
            std::cout << " parent3: "
                      << ((parent(data.begin(), data.begin()+3) == (data.begin()+1)) ?
                          "passed" : "FAILED")
                      << "\n";

        // index 4 -> parent 1
        if (data.size() > 4)
            std::cout << " parent4: "
                      << ((parent(data.begin(), data.begin()+4) == (data.begin()+1)) ?
                          "passed" : "FAILED")
                      << "\n";

        // the soft general check
        bool failed = false;
        for (std::vector<int>::iterator i = data.begin() + 1;
             i != data.end();
             ++i)
        {
            std::vector<int>::iterator p    = parent(data.begin(),i);
            std::vector<int>::iterator kid1 = child(data.begin(),p);
            std::vector<int>::iterator kid2 = kid1 + 1;

            bool bad = (kid1 != i) and (kid2 != i);
            if (bad)
                std::cout << " indexing bad at"
                          << " i:" << (i-data.begin())
                          << " p:" << (p-data.begin())
                          << " k1:" << (kid1-data.begin())
                          << " k2:" << (kid2-data.begin())
                          << "\n";

            failed |= bad;
        }
        std::cout << " general heap indexing:" << (failed ? "FAILED" : "passed") << "\n";
    }

    // check heap making
    makeHeap(data.begin(),
             data.end());

    bool failed = false;

    if (data.size() > 0)
        std::cout << " # " << data[0];

    for (int i = 1; i < data.size(); ++i)
    {
        int p = (i-1)/2;
        bool bad = (data[i] > data[p]);
        failed |= bad;
        std::cout << (bad ? "|" : " ") << data[i];
    }

    std::cout <<" " << (failed ? "FAILED" : "passed")
              << "\n";
}

int testRotate()
{
    std::vector<int> data({1,2,3,4,5,6,7,8});

    std::vector<int> data1(data);
    rotate(data1.begin(),
           data1.begin()+2,
           data1.end());

    for (std::vector<int>::iterator a = data1.begin();
         a != data1.end();
         ++a)
        std::cout << " " << *a;
    std::cout << "\n";

    std::cout << (data1 == std::vector<int>({3,4,5,6,7,8,1,2}) ? "passed" : "FAILED")
              << "\n";

    std::vector<int> data2(data);
    rotate(data2.begin(),
           data2.begin()+6,
           data2.end());

    for (std::vector<int>::iterator a = data2.begin();
         a != data2.end();
         ++a)
        std::cout << " " << *a;
    std::cout << "\n";

    std::cout << (data2 == std::vector<int>({7,8,1,2,3,4,5,6}) ? "passed" : "FAILED")
              << "\n";
}

int test_one(std::vector<int> data,
             std::function<void (std::vector<int>&)> sort)
{
    std::vector<int> list(data);
    sort(list);
    check(list);
}

void test(std::initializer_list<int> data)
{
    std::cout << "testing heap parts\n";
    testHeap(data);

    std::cout << "testing sorts\n";
    test_one(data,insertionSort);
    test_one(data,selectionSort);
    test_one(data,mergeSortWrap);
    test_one(data,inplaceMergeSortWrap);
    test_one(data,quickSortWrap);
    test_one(data,stableQuickSortWrap);
    test_one(data,heapSortWrap);
}

int main()
{
    std::cout << "testing rotate\n";
    testRotate();

    std::cout << "test checker\n";
    std::vector<int> list1({6,7,8,8,6,5});
    check(list1);

    std::cout << "testing sets\n";
    test({});
    test({1});
    test({1,1}); // sort breaker

    test({1,2});
    test({2,1});

    test({6,7,8,8,6,5});

    test({2,2,2,2,2,2,2,2,2,2,2,2,2}); // sort breaker...

    test({3,5,2,6,8,3,2,6,5,3,4,8,7});
}

Saturday, July 2, 2016

google interview: uniform sampling from an infinite stream

//http://www.impactinterview.com/2009/10/140-google-interview-questions/#software_engineer
// You have a stream of infinite queries (ie: real time Google search queries
// that people are entering).  Describe how you would go about finding a good
// estimate of 1000 samples from this never ending set of data and then write
// code for it.

#include <iostream>
#include <iomanip>
#include <vector>
#include <ctime>

// Reservoir Sampling
// did it here http://code-slim-jim.blogspot.jp/2010/06/reservoir-sampling.html
// but wow 2010 is so long ago..

// ok so the explanation:
// you have stream and u need to sample from N "good" items it..  assumably
// "good" means uniform here so until you have the first N items u just take
// every thing
//
// now the tricky part when u get the N+1 item each item is therefore supposed
// to have a N/(N+1) probability of being in the output so the new item has a
// N/(N+1) chance of selection and all the existing items have a 1/N chance of
// begin rejected..
//
// then we can generalize to the (N+x)-th sample.  when u get the N+x item each
// item is therefore supposed to have a N/(N+x) probability of being in the
// output so the new item has a N/(N+x) chance of selection and all the
// existing items have a 1/N chance of begin rejected


template <typename Type>
class ReservoirSample
{
 public:
    typedef std::vector<Type> Samples;

 private:
    std::size_t       size_;
    std::size_t       count_;
    std::vector<Type> samples_;

 public:
    ReservoirSample(unsigned int size) :
        size_(size),
        count_(0),
        samples_()
    {}

    float rand()
    {
        return static_cast<float>(std::rand() % 1000) / 1000.0;
    }

    void sample(Type& item)
    {
        count_ += 1;
        if (samples_.size() < size_)
        {
            samples_.push_back(item);
        }
        else
        {
            float chance = static_cast<float>(size_)/static_cast<float>(count_);

            if (rand() < chance)
            {
                std::size_t idx = size_ * rand();
                samples_[idx] = item;
            }
        }
    }

    std::vector<Type>& samples()
    {
        return samples_;
    }
};

void test(int selection,
          int streamSize)
{
    // testing random stuff can be a pain.. you either de-random it
    // or you repeat and confirm the expected distrubution (uniform...)
    const int cycles = 100000;
    std::vector<int> freq(streamSize,0);

    for (int j = 0; j < cycles; ++j)
    {
        ReservoirSample<int> sampler(selection);

        for (int i = 0; i < streamSize; ++i)
        {
            sampler.sample(i);
        }

        const std::vector<int>& samples = sampler.samples();
        for ( std::vector<int>::const_iterator sit = samples.begin();
              sit != samples.end();
              ++sit)
        {
            ++(freq[*sit]);
        }
    }

    // each number will come up once in a cycle and "selection" are choosen out of "streamSize"
    float expected
        = (static_cast<float>(cycles)
           * static_cast<float>(selection))
        / static_cast<float>(streamSize);

    // and lets give it +/- 5% .. (which will fail now and then..)
    float expectMin = 0.95 * expected;
    float expectMax = 1.05 * expected;

    for (int i = 0; i < streamSize; ++i)
    {
        bool passed = freq[i] > expectMin and freq[i] < expectMax;
        std::cout << std::setw(3) << i << ":" << freq[i]
                  << " " << (passed ? "passed" : "FAILED")
                  << "\n";
    }
    std::cout << "\n";
}

int main()
{
    std::srand(std::time(0));

    test(1,30);
    test(3,30);
    test(7,30);
    test(15,30);
}

0:3307 passed
  1:3302 passed
  2:3237 passed
  3:3344 passed
  4:3308 passed
  5:3233 passed
  6:3317 passed
  7:3257 passed
  8:3230 passed
  9:3407 passed
 10:3310 passed
 11:3267 passed
 12:3303 passed
 13:3414 passed
 14:3267 passed
 15:3327 passed
 16:3333 passed
 17:3371 passed
 18:3351 passed
 19:3364 passed
 20:3306 passed
 21:3364 passed
 22:3402 passed
 23:3440 passed
 24:3326 passed
 25:3390 passed
 26:3351 passed
 27:3451 passed
 28:3289 passed
 29:3432 passed

  0:10056 passed
  1:9997 passed
  2:10060 passed
  3:10012 passed
  4:9919 passed
  5:10021 passed
  6:9984 passed
  7:9945 passed
  8:9987 passed
  9:10040 passed
 10:9919 passed
 11:9965 passed
 12:9918 passed
 13:9904 passed
 14:9864 passed
 15:9876 passed
 16:10043 passed
 17:9945 passed
 18:10064 passed
 19:10025 passed
 20:9878 passed
 21:10021 passed
 22:10108 passed
 23:9796 passed
 24:10100 passed
 25:10025 passed
 26:10206 passed
 27:10003 passed
 28:10172 passed
 29:10147 passed

  0:23484 passed
  1:23126 passed
  2:23259 passed
  3:23345 passed
  4:23290 passed
  5:23339 passed
  6:23497 passed
  7:23348 passed
  8:23290 passed
  9:23179 passed
 10:23603 passed
 11:23154 passed
 12:23223 passed
 13:23333 passed
 14:23394 passed
 15:23193 passed
 16:23546 passed
 17:23219 passed
 18:23401 passed
 19:23201 passed
 20:23246 passed
 21:23240 passed
 22:23621 passed
 23:23224 passed
 24:23410 passed
 25:23376 passed
 26:23518 passed
 27:23268 passed
 28:23348 passed
 29:23325 passed

  0:49663 passed
  1:49911 passed
  2:50401 passed
  3:49617 passed
  4:49691 passed
  5:50278 passed
  6:49804 passed
  7:49940 passed
  8:50497 passed
  9:49687 passed
 10:49582 passed
 11:50552 passed
 12:49871 passed
 13:49670 passed
 14:50304 passed
 15:50184 passed
 16:49988 passed
 17:49965 passed
 18:50077 passed
 19:49945 passed
 20:50290 passed
 21:50173 passed
 22:49972 passed
 23:50051 passed
 24:49943 passed
 25:50053 passed
 26:50288 passed
 27:49852 passed
 28:50009 passed
 29:49742 passed

Google interview question: write an execl col label to integer converter

// http://www.impactinterview.com/2009/10/140-google-interview-questions/#software_engineer
// Write a function (with helper functions if needed) called to Excel that takes an excel column value
// (A,B,C,D…AA,AB,AC,… AAA..) and returns a corresponding integer value (A=1,B=2,… AA=26..).

#include <iostream>
#include <stdexcept>

// the function "toExecl" will convert *from* excel format to integers... ignoring that confusion...
int fromExcel(std::string val)
{
    // there is an oddity here
    // 0A -> A but A0 -> doesnt exist! note the +1 in the conversion line as a result

    int out = 0;
    for (std::string::iterator vit = val.begin();
         vit != val.end();
         ++vit)
    {
        if (*vit < 'A' or *vit > 'Z')
            throw std::runtime_error("doesnt look like an execl col");
        out = (out*26) + (*vit - 'A' + 1);
    }

    return out;
}

int test(std::string in, int exp)
{
    int out = fromExcel(in);

    std::cout << (out == exp ? "passed" : "FAILED")
              << " in:" << in
              << " out:" << out
              << " exp:" << exp
              << "\n";
}

int main()
{
    try
    {
        test("012", 0);
        std::cout << "FAILED: did not detect bad inputs\n";
    }
    catch (std::exception& e)
    {
        std::cout << "pass: throw check worked\n";
    }

    test(  "A",                     1);
    test(  "Z",                    26);
    test( "AA",             1*26 +  1);
    test( "AZ",             1*26 + 26);
    test( "BA",             2*26 +  1);
    test( "ZZ",            26*26 + 26);
    test("AAA", 1*26*26 +   1*26 +  1);
    test("AZZ", 1*26*26 +  26*26 + 26);
    test("ZAA",26*26*26 +   1*26 +  1);
    test("ZZZ",26*26*26 +  26*26 + 26);
}

Output

pass: throw check worked
passed in:A out:1 exp:1
passed in:Z out:26 exp:26
passed in:AA out:27 exp:27
passed in:AZ out:52 exp:52
passed in:BA out:53 exp:53
passed in:ZZ out:702 exp:702
passed in:AAA out:703 exp:703
passed in:AZZ out:1378 exp:1378
passed in:ZAA out:17603 exp:17603
passed in:ZZZ out:18278 exp:18278

Google interview questions: print matching chars in order of first string


// http://www.impactinterview.com/2009/10/140-google-interview-questions/#software_engineer
// Write a function f(a, b) which takes two character string arguments and returns a string containing
// only the characters found in both strings in the order of a. Write a version which is order N-squared
// and one which is order N.

#include <iostream>
#include <string>
#include <functional>
#include <vector>
#include <cstring>

std::string match_the_stupid_way(std::string a, std::string b)
{
    // the silly way is just for each char in string a O(n)
    // look through string b O(n) double loop hence O(n^2)

    std::string res;

    for (std::string::iterator ait = a.begin();
         ait != a.end();
         ++ait)
    {
        // ok match the first each char
        std::string::iterator bit = b.begin();
        while (bit != b.end() and
               *bit != *ait)
        {
            ++bit;
        }

        if (bit != b.end())
            res += *ait;
    }
    return res;
}

std::string match_the_fast_way(std::string a, std::string b)
{
    // O(N) .. this imples reading over a string (which takes O(N)) and
    // digesting into an O(1) data struture and then check over the other...

    // well an O(1) data structure is a hash or table..
    std::vector<bool> found_char(256,false);

    // the O(N) read is
    for (std::string::iterator bit = b.begin();
         bit != b.end();
         ++bit)
    {
        found_char[*bit] = true;
    }

    // and the second O(N) check is
    std::string res;
    for (std::string::iterator ait = a.begin();
         ait != a.end();
         ++ait)
    {
        if (found_char[*ait])
            res += *ait;
    }

    return res;
}

void test(std::function<std::string (std::string,std::string)> func,
          std::string a,
          std::string b,
          std::string exp)
{
    std::string out = func(a,b);

    std::cout << (out == exp ? "passed" : "FAILED")
              << " a:\""     << a
              << "\" b:\""   << b
              << "\" out:\"" << out
              << "\" exp:\"" << exp
              << "\"\n";
}

void test_set(std::function<std::string (std::string,std::string)> func)
{
    test(func, "a", "cba", "a");
    test(func, "aaa", "cba", "aaa");  // hmm question isnt clear on repeats.. assume it does repeats
    test(func, "abcdefg", "gfedcba", "abcdefg");
    test(func, "the quick brown fox", "peter piper picked a pepper", "te ick r ");
    std::cout << "*** SET DONE ***\n";
}

int main()
{
    test_set(match_the_stupid_way);
    test_set(match_the_fast_way);
}

Output is

passed a:"a" b:"cba" out:"a" exp:"a"
passed a:"aaa" b:"cba" out:"aaa" exp:"aaa"
passed a:"abcdefg" b:"gfedcba" out:"abcdefg" exp:"abcdefg"
passed a:"the quick brown fox" b:"peter piper picked a pepper" out:"te ick r " exp:"te ick r "
*** SET DONE ***
passed a:"a" b:"cba" out:"a" exp:"a"
passed a:"aaa" b:"cba" out:"aaa" exp:"aaa"
passed a:"abcdefg" b:"gfedcba" out:"abcdefg" exp:"abcdefg"
passed a:"the quick brown fox" b:"peter piper picked a pepper" out:"te ick r " exp:"te ick r "
*** SET DONE ***

google interview question: build a random generator

Google interview questions

// http://www.impactinterview.com/2009/10/140-google-interview-questions/
// Given a function which produces a random integer in the
// range 1 to 5, write a function which produces a random
// integer in the range 1 to 7.

#include <cstring>
#include <cstdlib>
#include <iostream>
#include <functional>
#include <chrono>

// Note stuff the [1,5] thats just +1. -1 everywhere and a waste of time lets go [0,4]

int rand5()
{
    return std::rand() % 5;
}

// ok so the problem is that you need to combine the rand5 to make rand7
// the question doesnt say anything about *uniform* distribution so in theroy
// u just call rand5 2 times and mod 7.. lets call this the weasel solution
// cause the question was badly written

int rand7_weasel()
{
    return (rand5() + rand5()) % 7;

}

// ok so now the questioner has realized his mistake and corrected.
// a "uniform" rand7() generator
//
// The above is not uniform because if use two 6 sided dice and added the rolls
// the most common number is 7. Ie 7 has the most number of combos that add to it.
//
// So lets simplify and assume for the moment we have a more natural problem
// we can generate digits [0,9] uniformly and need to create a uniform number
// from [0,999].. well the answer is much more clear now you generate 3 numbers
// and use them as the digits of the final number
//
// The same idea works for the [0,4], we have 5 numbers so we simply use base 5
// instead of base 10. This way we can make a very large uniform distributed number
//
// once we have a very large range of uniformly distributed numbers we can module
// it by a much smaller number the binning error will be negligible

// lets minimize the binning error by making the max number possible (overflow it once)
// 
// 2^32 = 5^N
// 32 log2 = N log5
// N = 13.78

int rand7_overkill()
{
    unsigned int sum = 0;
    for (int i = 0; i < 14; i++)
        sum = (sum*5) + rand5();

    return sum % 7;
}

// of course if you know your maths the central limit theorem says that for most distributions
// when you take many independent samples and summate them the resulting distribution will 
// approach a another stable one.. for uniform distributions thats the normal distribution.
//  refer to: http://demonstrations.wolfram.com/CentralLimitTheoremForTheContinuousUniformDistribution/
//
// Then add to this the fact that a small enough modulus on a normal distribution (much
// much less that its variance) will make the resulting distribution look somewhat uniform
//  ie the wrapped normal distribution where sigma^2 approaches infinite
//  refer to: https://en.wikipedia.org/wiki/Wrapped_normal_distribution
//
// So technically this works.. but a smart interviewer is going to ask you why it works.. and
// your likely to choke on it.. (i know the maths would kill me during an interview)

int rand7_central_limit_weasel()
{
    unsigned int sum = 0;
    for (int i = 0; i < 14; i++)
        sum += rand5();

    return sum % 7;
}

// ok so now we have a working solution.. but guess what your interviewing for google..
// and they want to curve ball it.. so they ask the following.. can you improve the
// performance?.. looping 14 times is too much.. a very "innocent" question
//
// now you might think ok well just reduce the number of iterations it that does it..
// but that is not really going to work well
//
// So here is the trick.. the real question is about testing your knowledge of pseudo
// random generators.. so just build one.. now a pseudo random generator isn't anything
// special but you have to hack one out in an interview and i don't have primes. pi and
// what not memorized by heart...
//
// What i do know is how to generate a uniform distribution that ranges from [0,infinite]..
// so all u have to do is make one big continuous rolling uniform number from the above code
// and then sample from it as needed..
//
// now this does raise some questions about the independence of the samples... but its faster...

int rand7_pseudo()
{
    // inject a the random sample to seeded the hash
    static unsigned int sum = rand5();

    // the uniform sum on with the prior...
    sum = (sum*5) + rand5();

    return sum % 7;
}

void test(std::function<int()> randy)
{
    int freq[7] = {0};
    std::memset(freq, 0, sizeof(freq));

    auto start = std::chrono::steady_clock::now();

    for (int i = 0; i < 100000; ++i)
    {
        ++freq[randy()];
    }

    auto finish = std::chrono::steady_clock::now();
    double elapsed_seconds = std::chrono::duration_cast<
        std::chrono::duration<double> >(finish - start).count();

    for (int i=0;i<7;++i)
    {
        std::cout << "\n " << i << ":" << freq[i];
    }
    std::cout << "\n speed:" << elapsed_seconds
              << "\n";
}


int main()
{
    test(rand7_weasel);
    test(rand7_overkill);

    test(rand7_central_limit_weasel);
    test(rand7_pseudo);
}

and output looks like this:

 0:11949
 1:12073
 2:11885
 3:16060
 4:19993
 5:16024
 6:12016
 speed:0.0190048

 0:14130
 1:14541
 2:14177
 3:14296
 4:14157
 5:14400
 6:14299
 speed:0.056561

 0:14264
 1:14538
 2:14425
 3:14064
 4:14476
 5:14076
 6:14157
 speed:0.0409873

 0:14272
 1:14230
 2:14308
 3:14110
 4:14214
 5:14422
 6:14444
 speed:0.0131068

Tuesday, June 7, 2016

c++ like template AutoGrad - How to compute the Derivative Of a function using templates

You may have heard of a system called AutoGrad. It is a system that when given an equation computes the derivate of it automatically. Here is an example of how to build the same system using c++ templates.

I was inspired to write this after i realized i messed up the differentiation of the for the neural networks back prop path in my last post.

Clearly you see where im headed with this.. a few more mods to make it work with tensors and ill have full autograd system for c++ that can automatically compute the backdrop of a neural network.. and will also compile in c4driod.. and then toss in some thrust code and it will be CUDA capable able to use PC with GPUs

from smartphone to GPU in a few short steps...

Anyway here is the proof of concept code:

output looks like
x:2 m:3 c:4 p:6 y:10 z:16 q:100
 dm/dx:0 dx/dx:1 dp/dx:1 dy/dx:3 dz/dx:12 dq/dx:60

// compile with "g++ std=c++11 ..."
// or in c4driod on an andriod phone..

#include <tuple>
#include <iostream>
#include <cmath>

template <int VALUE>
struct Const
{
    template< typename... Types >
    static double op(std::tuple<Types...>& params )
    {
        return VALUE;
    }
};

template <int ID>
struct Var
{
    template< typename... Types >
    static double op(std::tuple<Types...>& params )
    {
        return std::get<ID>(params);
    }
};

template <typename L, typename R>
struct OpAdd
{
    template< typename... Types >
    static double op(std::tuple<Types...>& params)
    {
        return L::op(params) + R::op(params);
    }
};

template <typename L, typename R>
struct OpMult
{
    template< typename... Types >
    static double op(std::tuple<Types...>& params)
    {
        return L::op(params) * R::op(params);
    }
};

template <typename L, int POW>
struct OpPower
{
    template< typename... Types >
    static double op(std::tuple<Types...>& params)
    {
        return std::pow(L::op(params), POW);
    }
};

template <typename Num, typename Denum>
struct DerivativeOf {};

template <int V, int ID_D>
struct DerivativeOf<Const<V>, Var<ID_D> >
{
    typedef Const<0> Type;
};

template <int ID_N, int ID_D>
struct DerivativeOf<Var<ID_N>, Var<ID_D> >
{
    typedef Const<0> Type;
};

template <int ID_N>
struct DerivativeOf<Var<ID_N>, Var<ID_N> >
{
    typedef Const<1> Type;
};

template <typename L, typename R, typename D>
struct DerivativeOf<OpAdd<L,R>, D >
{
    typedef OpAdd<typename DerivativeOf<L,D>::Type,
                  typename DerivativeOf<R,D>::Type> Type;
};

template <typename L, typename R, typename D>
struct DerivativeOf<OpMult<L,R>, D >
{
    typedef OpAdd<OpMult<typename DerivativeOf<L, D>::Type, R>,
                  OpMult<L, typename DerivativeOf<R, D>::Type> > Type;
};

template <typename L, int POW>
struct DerivativeOf<OpPower<L,POW>, L >
{
    typedef OpMult<Const<POW>, OpPower<L, POW-1> > Type;
};

template <typename L, int POW, typename D>
struct DerivativeOf<OpPower<L,POW>, D >
{
    typedef OpMult<OpMult<Const<POW>, OpPower<L, POW-1> >,
                   typename DerivativeOf<L, D>::Type>  Type;
};

int main()
{
    typedef Var<0> X;
    typedef Var<1> M;
    typedef Var<2> C;
    typedef OpAdd<X,C> P;                      // p = x + c
    typedef OpAdd<OpMult<M,X>,C> Y;            // y = m*x + c
    typedef OpAdd<OpMult<M,OpPower<X,2>>,C> Z; // z = m*x^2 + c
    typedef OpPower<OpAdd<OpMult<M,X>,C>,2> Q; // q = = y^2 = (m*x + c)^2

    typedef DerivativeOf<M, X>::Type  dM_dX;   // dm/dx = 0
    typedef DerivativeOf<X, X>::Type  dX_dX;   // dx/dx = 1
    typedef DerivativeOf<P, X>::Type  dP_dX;   // dP/dx = 1
    typedef DerivativeOf<Y, X>::Type  dY_dX;   // dY/Dx = m
    typedef DerivativeOf<Z, X>::Type  dZ_dX;   // dZ/dx = 2x
    typedef DerivativeOf<Q, X>::Type  dQ_dX;   // dQ/dx = dq/dy + dy/dx = 2*(m*x + c)*m

    X x;
    M m;
    C c;
    P p;
    Y y;
    Z z;
    Q q;

    dX_dX  dx_dx;
    dM_dX  dm_dx;
    dP_dX  dp_dx;
    dY_dX  dy_dx;
    dZ_dX  dz_dx;
    dQ_dX  dq_dx;

    std::tuple<double,double,double> params = std::make_tuple(2.0,3.0,4.0);

    std::cout << " x:" << x.op(params)
              << " m:" << m.op(params)
              << " c:" << c.op(params)
              << " p:" << p.op(params)
              << " y:" << y.op(params)
              << " z:" << z.op(params)
              << " q:" << q.op(params)
              << "\n";

    std::cout << " dm/dx:"  << dm_dx.op(params)
              << " dx/dx:"  << dx_dx.op(params)
              << " dp/dx:"  << dp_dx.op(params)
              << " dy/dx:"  << dy_dx.op(params)
              << " dz/dx:"  << dz_dx.op(params)
              << " dq/dx:"  << dq_dx.op(params)
              << "\n";
}