import numpy as np

"""
Explanation: CSE 6040, Fall 2015 [16]: Numpy (and SciPy)
Numpy is a Python module that provides fast primitives for multidimensional arrays. It's well-suited to implementing numerical linear algebra algorithms, and for those can be much faster than Python's native list and dictionary types when you only need to store and operate on numerical data.
Some of the material from this lesson is lifted from the following more comprehensive tutorial: link
Quick demo. The recommended importing convention is:
End of explanation
"""


a_1d = np.array ([0, 1, 2, 3]) # a vector
print a_1d

b_1d = np.array ([4, 5, 6, 7]) # another vector
print b_1d

print a_1d + b_1d

print 5*a_1d

print a_1d**2

"""
Explanation: Numpy provides some natural types and operations on arrays. For instance:
End of explanation
"""


# Append '?' to get help on a specific routine
np.array?

# Search for key text
np.lookfor ("creating array")

# Wildcard search
np.con*?

"""
Explanation: Getting help. By the way, if you need help getting documentation from within this notebook, here are some handy shortcuts.
End of explanation
"""


n = 1000000

%timeit L = range (n)

L = range (n)
%timeit [i**2 for i in L]

np.arange (10) # Moral equivalent to `range`

%timeit A = np.arange (n)

A = np.arange (n)
%timeit A**2

"""
Explanation: Why bother with Numpy? A motivating example
We already have lists and dictionary types, which are pretty easy to use and very flexible. So why bother with this special type?
One reason to consider Numpy is that it "can be much faster," as noted above. But how much faster is that?
End of explanation
"""


from random import gauss # Generates random numbers from a Gaussian
from math import sqrt # Computes the square root of a number

n = 1000000
X_py = [gauss (0, 1) for i in range (n)]
X_np = np.array (X_py)

print ("==> Native Python lists:")
%timeit ...

print ("\n==> Numpy:")
%timeit ...

"""
Explanation: Exercise. Recall the definition of the 2-norm of a vector (or Euclidean length of a vector) from Da Kuang's notes on linear algebra. Compare its speed when using native Python lists versus Numpy arrays.

Hint: For Numpy, look for a routine that computes the norm.
End of explanation
"""


# Create a two-dimensional array of size 3 rows x 4 columns:
B = np.array([[0, 1, 2, 3],
              [4, 5, 6, 7],
              [8, 9, 10, 11]])

print B.ndim # What does this do?
print B.shape # What does this do?
print len (B) # What does this do?

C1 = [[0, 1, 2, 3],
      [4, 5, 6, 7],
      [8, 9, 10, 11]]

C2 = [[12, 13, 14, 15],
      [16, 17, 18, 19],
      [20, 21, 22, 23]]

C = np.array ([C1, C2])

print C.ndim
print C.shape
print len (C)

"""
Explanation: Creating multidimensional arrays
Beyond simple arrays, Numpy supports multidimensional arrays. To do more than one dimension, call numpy.array() but nest each new dimension within a list.
Huh? Let's look at some examples.
End of explanation
"""


print np.linspace (0, 1, 10)

print np.linspace (0, 1, 10, endpoint=False)

"""
Explanation: Besides arange(), you can also define an interval and a number of points. What does the following code do?
End of explanation
"""


print np.ones ((3, 4))

print np.zeros ((3, 4))

print np.eye (3)

print np.diag ([1, 2, 3])

"""
Explanation: There are routines for creating various kinds of structured matrices as well, which are similar to those found in MATLAB and Octave.
End of explanation
"""


n = 3
I = np.eye (n)

print ("I = eye(n):")
print (I)

u = [1] * n
I_u = np.diag (u)

print ("u:")
print u
print ("\ndiag (u):")
print I_u

assert np.all (D_u == E)

"""
Explanation: Exercise. The following code creates an identity matrix in two different ways, which are found to be equal according to the assertion. But in fact there is a subtle difference; can you spot it?
End of explanation
"""


A = np.empty ((3, 4)) # An empty 3 x 4 matrix
print A

"""
Explanation: You can also create empty (uninitialized) arrays. What does the following produce?
End of explanation
"""


# Recall: C
print C

"""
Explanation: Indexing and slicing
The usual 0-based slicing and indexing notation you know and love from lists is also supported for Numpy arrays. In the multidimensional case,  including their natural multidimensional analogues with index ranges separated by commas.
End of explanation
"""


print C[0, 2, :]

"""
Explanation: Exercise. What part of C will the following slice extract? Run the code to check your answer.
End of explanation
"""


print C[1, 0, ::-1]

"""
Explanation: Exercise. What will the following slice return? Run the code to check your answer.
End of explanation
"""


Z = ... # INSERT YOUR CODE HERE

print ("\n==> Orange slice:")
print Z[...] # Edit to print the orange-ish subset

print ("\n==> Red slice:")
print Z[...] # Edit to print the red-ish subset

print ("\n==> Blue slice:")
print Z[...] # Edit to print the blue-ish subset

print ("\n==> Green slice:")
print Z[...] # Edit to print the green-ish subset

"""
Explanation: Exercise. Consider the following $6 \times 6$ matrix, which has 4 different subsets highlighted.
<img src="https://raw.githubusercontent.com/rvuduc/cse6040-ipynbs/master/16--numpy-slicing-exercise.png" alt="Exercise: Extract these slices" width="240">
Write some code to generate this matrix, named Z. Then, for each subset illustrated above, write an indexing or slicing expression that extracts the subset.
End of explanation
"""


# Rich will demo this in class

"""
Explanation: Incidentally, there is a very cute way to create the above matrix. I would have never guessed this method, so I've put some print statements so you can see what it is doing.
End of explanation
"""


print ("==> Recall C: %s" % str (C.shape))
print C

C_view = C[1, 0:3:2, 1:4:2] # Question: What does this produce?
print ("==> C_view: %s" % str (C_view.shape))

C_view[:, :] = -C_view[::-1, ::-1] # Question: What does this do?

print C

"""
Explanation: Slices are views
To help save memory, when you slice a Numpy array, you are actually creating a view into that array. That means modifications through the view will modify the original array.
End of explanation
"""


C_copy = C[1, 0:3:2, 1:4:2].copy ()
C_copy[:, :] = -C_copy[::-1, ::-1]

print ("==> C_view:")
print (C_view)

print ("\n==> C_copy:")
print (C_copy)

"""
Explanation: You can force a copy using the .copy() method:
End of explanation
"""


print ("C and C_view share memory: %s" % np.may_share_memory (C, C_view))
print ("C and C_copy share memory: %s" % np.may_share_memory (C, C_copy))

"""
Explanation: And to check whether two Numpy array variables point to the same object, you can use the numpy.may_share_memory() function:
End of explanation
"""


from math import sqrt

n = 20
is_prime = np.ones (n+1, dtype=bool) # the "sieve"

# Initial values
is_prime[0:2] = False
is_prime[2:] = True

# Sieving loop
for i in range (2, int (sqrt (n))):
    # Fill in your code here
    is_prime[2*i::i] = False
    
# Prints your primes
print ("==> Primes through %d:" % n)
print np.nonzero (is_prime)[0]

"""
Explanation: Exercise. Complete the prime number sieve algorithm, which is illustrated below.
<img src="https://raw.githubusercontent.com/rvuduc/cse6040-ipynbs/master/16--numpy-prime-sieve.png" alt="Exercise: Extract these slices" width="480">
That is, given a positive integer $n$, the algorithm iterates from $i \in {2, 3, 4, \ldots, \left\lfloor\sqrt{n}\right\rfloor}$, repeatedly "crossing out" values that are strict multiples of $i$. "Crossing out" means maintaining an array of, say, booleans, and setting values that are multiples of $i$ to False.
End of explanation
"""


np.random.seed(3)
x = np.random.random_integers(0, 20, 15) # 15 random ints in [0, 20)
print x

# Find all positive multiples of 3
mask_mult_3 = (x % 3 == 0) & (x > 0)
print mask_mult_3
print x[mask_mult_3]

# Pull out an arbitrary subset of elements
inds = np.array ([3, 7, 8, 12])
print x[inds]

"""
Explanation: Indirect addressing
Another common indexing operation is using either a boolean mask or a set of integer indices.
End of explanation
"""


n = 5000
A_rowmaj = np.ones ((n, n), order='C') # row-major (C/C++ convention)
A_colmaj = np.ones ((n, n), order='F') # column-major (Fortran convention)

"""
Explanation: Matrix storage: Column-major versus row-major layouts
For linear algebra, we will be especially interested in 2-D arrays, which we will use to store matrices. There is a final subtle performance issue related to how matrices are stored in memory.
By way of background, physical storage---whether it be memory or disk---is basically one big array. And because of how physical storage is implemented, it turns out that it is much faster to access consecutive elements in memory than, say, to jump around randomly.
A matrix is a two-dimensional object. Thus, when it is stored in memory, it must be mapped in some way to the one-dimensional physical array. There are many possible mappings, but the two most common conventions are known as the row-major and column-major layouts:
<img src="https://raw.githubusercontent.com/rvuduc/cse6040-ipynbs/master/16--numpy-matrix-layout.png" alt="Exercise: Extract these slices" width="640">
In Numpy, you can ask for one layout or the other. The most common convention in numerical linear algebra libraries is column-major layout, which is also the Fortran language convention; by contrast, when you create a two-dimensional array in languages like C and C++, the convention is row-major layout.
The default in Numpy is C order (row-major), even though the linear algebra library convention is Fortran order (column-major). You can request either order using the order parameter. For linear algebra, we recommend using the column-major convention.
End of explanation
"""


def scale_colwise (A):
    """Given a matrix `A`, visits each column `A[:, j]` and scales it by `j`."""
    # Fill in this code
    pass

# Measure time to scale a row-major input column-wise
%timeit scale_colwise (A_rowmaj)

# Measure time to scale a column-major input column-wise
%timeit scale_colwise (A_colmaj)

"""
Explanation: Exercise. Write a function that iterates over the columns of a matrix and scales each column by its index. That is, column $A(:, j)$ should be replaced by $j*A(:, j)$.
Then use that routine to measure the difference in time when the input is row-major versus column-major.

Note: In Numpy, this operation can be easily done just by multiplying
End of explanation
"""


