Find subsequence with max product

max_prod_subsequence.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Find subsequence with maximum product in array\n",
    "\n",
    "## Description\n",
    "\n",
    "Given array of **N** numbers, find a subsecuence of **_k_** numbers with maximum product (multiple of all elements).\n",
    "\n",
    "## Solution\n",
    "\n",
    "We can solve in O(n) time. Simply calculate product of first _k_ elements, then move through the array, multiplying current (_i_-th) element with current product and dividing to (_i - k_)-th element. At each _i_-th position this will give us product of last _k_ elements. Since multiplication on 0 gives 0, there is no need to calculate sequences that include 0. Each time we encounter 0 in the array we start a new subsequence from the next element. If array does not contain _k_ subsequent non-zero elements, its maximum subsequent product is 0 and starts anywhere.\n",
    "\n",
    "Depending on the size of array, its elements and _k_ this algorithm may overflow quite quickly. Additionally using floating point calculations will introduce accumulating error that, however, should not affect result significantly. We will use power of Julia type system to create two generic functions for integer and floating point arrays.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First - integer array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = rand(Int8, 1000); k=13;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check number of zeroes in the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "count(i->i==0, A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "find_max_prod (generic function with 3 methods)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Julia dynamic dispatch magic: Arrays of integers will trigger method\n",
    "# with BigInt, arrays of floats will be processed as Float64, while\n",
    "# algorithm is written in generic form\n",
    "\n",
    "find_max_prod(A::AbstractArray{T}, k) where {T <:Integer} = find_max_prod(A,k,BigInt(0))\n",
    "find_max_prod(A::AbstractArray{T}, k) where {T <:AbstractFloat} = find_max_prod(A,k,Float64(0))\n",
    "function find_max_prod(A::AbstractArray, k, max_so_far)\n",
    "    # Using BigInt we can't overflow\n",
    "    len = 1; prod::typeof(max_so_far) = 1; max_ending_pos = -1\n",
    "    for (i,a) in enumerate(A)\n",
    "        if a != 0\n",
    "            prod *= a\n",
    "            if len < k\n",
    "                len += 1\n",
    "            else\n",
    "                if prod > max_so_far\n",
    "                    max_so_far = prod\n",
    "                    max_ending_pos = i\n",
    "                end\n",
    "                prod /= A[i-k+1]\n",
    "            end\n",
    "        else\n",
    "            len = 1; prod = 1\n",
    "        end\n",
    "    end\n",
    "    return (max_so_far, max_ending_pos)\n",
    "end                  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum product subsequence of 13 elements is 4375218897328327276061400 and starts at index 587\n"
     ]
    }
   ],
   "source": [
    "(product, i) = find_max_prod(A, k)\n",
    "println(\"Maximum product subsequence of $k elements is $product and starts at index $(i-k+1)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The subsequence is:\n",
      "Int8[-41, -79, 113, -127, -119, 113, -118, 83, -37, 31, 67, -93, -100]\n",
      "And its product is: 4375218897328327276061400\n"
     ]
    }
   ],
   "source": [
    "println(\"The subsequence is:\")\n",
    "println(A[(i-k+1):i])\n",
    "prod_2 = Base.prod(convert(Array{BigInt},A[(i-k+1):i]))\n",
    "println(\"And its product is: $prod_2\")\n",
    "@assert(product == prod_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And little benchmark:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  467.277 μs (11198 allocations: 316.66 KiB)\n"
     ]
    }
   ],
   "source": [
    "using BenchmarkTools\n",
    "@btime find_max_prod(A, k);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now for floats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normal distribution -1:1\n",
    "Af = randn(Float32, 1000); k=13;\n",
    "\n",
    "# Scale up a little bit to get positive numbers\n",
    "Af *= 10;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia's dynamic dispatch will choose appropriate method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum product subsequence of 13 elements is 2.544447215909553e13 and starts at index 416\n"
     ]
    }
   ],
   "source": [
    "(product, i) = find_max_prod(Af, k)\n",
    "println(\"Maximum product subsequence of $k elements is $product and starts at index $(i-k+1)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The subsequence is:\n",
      "Float32[-24.5258, 11.9536, 14.3918, 12.4983, -14.8106, -2.193, 12.2479, 10.3853, -5.10423, -15.2527, 9.37286, 17.1246, -9.34643]\n",
      "And its product is: 2.5444473e13\n"
     ]
    }
   ],
   "source": [
    "println(\"The subsequence is:\")\n",
    "println(Af[(i-k+1):i])\n",
    "prod_2 = Base.prod(Af[(i-k+1):i])\n",
    "println(\"And its product is: $prod_2\")\n",
    "@assert(product ≈ prod_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  4.773 μs (1 allocation: 32 bytes)\n"
     ]
    }
   ],
   "source": [
    "@btime find_max_prod(Af, k);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see using floats is approximately 100 times more efficient than using BigInts."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.0.0",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}