from __future__ import division

import numpy as np
import scipy as sp
import matplotlib as mpl
import matplotlib.pyplot as plt

#IPython magic command for inline plotting
%matplotlib inline
#a better plot shape for IPython
mpl.rcParams['figure.figsize']=[15,3]

#Example to find the area of a square
def areaSquare(a):
    return a ** 2

print areaSquare(3)



x=np.linspace(1,10,10)
print x
print sumArray(x)

#Pen.py 
class Pen(object):
    def __init__ (self):
        self.state = False

    def click(self):
        self.state = not self.state

    def write(self, in_str):
        if (self.state):
            print "Writing message '%s'.\n"%in_str;
        else:
            print "Oops!  The pen is retracted.\n";


bic = Pen();

bic.click();
bic.write("hello world");
print bic.state

my_string = 'sunday sunday sunday'
my_string.upper()

'philadelphia'.capitalize()

# The following use of a member method `join` is astonishingly useful.
elements = ['fire', 'earth', 'water', 'air']
', '.join(elements)

# As is the use of `split`, particularly in data processing.
pmu_data = '03-Apr-2014 00:55:07.500,60.0347,Good,0,Good,0,Good,121.275,Good,102.8688,Good,60.0353,Good,0,Good,0,Good,60.035,Good,0,Good,0,Good,123.285,Good,14.69637,Good,123.783,Good,50.32861,Good,0,Good,60.035,Good,0.01,Good,0,Good,112.985,Good,162.6308,Good,0.1530308,Good,132.8789,Good'
pmu_data.split(',')

u'己所不欲，勿施於人。 (孔夫子)'
#print u'己所不欲，勿施於人。 (孔夫子)'

u'São Paolo, Brasil'
#print u'São Paolo, Brasil'

print  "\tTab me over"
print r"\tDon't tab me over"

from string import ascii_uppercase
ascii_uppercase

a = [66.25, 333, 333, 1, 1234.5]
a.insert(2, -1)
a.append(333)

dogs = ['newfoundland', 'labrador', 'chow-chow', 'mutt']
print 'I have %d dogs.'%len(dogs)
for dog in dogs:
    print 'One kind of dog is a %s.'%dog

squares = []
for x in range(10):
     squares.append(x**2)

squares

#Same as
squares = [x**2 for x in range(10)]

 matrix = [ [1, 2, 3, 4],
    [5, 6, 7, 8],
     [9, 10, 11, 12],
]
    
#List comprehension to evaluate transpose
[[row[i] for row in matrix] for i in range(4)]

myfirstdict = {'a':1, 'b':2, 'c':3}

myfirstdict['a']

bic.__dict__

class Dog:
    breed='abc'
    size=0
    name='xyz'

# Now let's create an object or an instance of the class Dog
d=Dog()

# The dot operator (.) gives you access to an object’s state and behavior (instance variables and methods).
d.name= 'Jack'
print d.name

# Now let us add a method to the class
class Dog:
    breed='abc'
    size=0
    name='xyz'
    
    def bark(self):
        print "Ruff Ruff!"

d=Dog()
d.bark()

x=Dog()
x.bark()

x.name = 'Wolverine'
x.breed = 'hound'
size = 24 # inches tall at the shoulder

class Dog:
    ''' Attributes include breed, size and name'''
    
    # constructor
    def __init__(self, name='xyz', breed='abc', size=0):
        self.name= name
        self.breed= breed
        self.size= size
    
    def bark(self):
        print "Ruff Ruff!"

y=Dog('Jack','lab',10)
print y.name
print y.breed

y=Dog()
print y.name
print y.breed



class Account:
    def __init__(self, holder, number, balance, credit_line=1500): 
        self.holder = holder 
        self.number = number 
        self.balance = balance
        self.credit_line = credit_line
    def deposit(self, amount): 
        self.balance = amount

acct = Account('Carrie Fisher', '492727ZED', 1300)
print acct.holder, acct.number, acct.balance, acct.credit_line

class Account:
    def __init__(self, holder, number, balance, credit_line=1500): 
        self._holder = holder 
        self.number = number 
        self.__balance = balance
        self.__credit_line = credit_line
    def deposit(self, amount): 
        self.balance = amount

class ChildAcc(Account):
    def getHolder(self):
        return self._holder # with two __ will fail, with one _ succeeds

acct = Account('Carrie Fisher', '492727ZED', 1300)
print acct.number

#print acct.holder, acct.number, acct.balance, acct.credit_line # will fail

chld = ChildAcc('Harrison Ford', 'C3POR2D2', 1250)
print chld.getHolder()

# How to access a 
class ChildAcc(Account):
    def __getHolder(self):
        return self._holder # with two __ will fail, with one _ succeeds

chld = ChildAcc('Frank Oz', 'Y0D4', 900)
print chld._ChildAcc__getHolder()

class Point:
    def __init( self, x=0, y=0):
        self.x = x
        self.y = y
        
    def __del__(self):
        class_name = self.__class__.__name__
        print class_name, "deleted"

pt1 = Point()

# The destructor `__del__` is called by `del`.
del pt1

class Rectangle(Canvas):
    def printRect(self):
        print "Inside Rectangle Class"

r = Rectangle(10, 5)
print r.area()

r.printRect()

class Square(Rectangle):
    def __init__(self,a):
        self.side=a
    def area(self):
        return self.side**2

s = Square(4)
print s.area()

s.printRect()



# Test your code as follows.
g = Graduate('John', '609248986', 'M', 'CS', 'PhD')
g.getDegree()

items = [1, 2, 3, 4, 5]
squared = []
for x in items:
	squared.append(x ** 2)

squared
[1, 4, 9, 16, 25]

def sqr(x): return x ** 2
list(map(sqr, items))

list(range(-5,5))
[-5, -4, -3, -2, -1, 0, 1, 2, 3, 4]
list( filter((lambda x: x < 0), range(-5,5)))
[-5, -4, -3, -2, -1]

#Sieve of Eratosthenes
nums = range(2, 50) 
for i in range(2, 8): 
    nums = filter(lambda x: x == i or x % i, nums)
print nums

reduce(lambda x, y: x+y, [1, 2, 3, 4, 5]) 

x = 5              # x <- 5
#x = eval('%d + 6' % x)     # x <- 11
#x = eval('abs(%d)' % -100) # x <- 100
#x = eval('x = 5')          # INVALID; assignment is not an expression.
#x = eval('if 1: x = 4')    # INVALID; if is a statement, not an expression.

x

exec 'print(5)'          # prints 5.
exec 'print(5)\nprint(6)'  # prints 5{newline}6.
exec 'if True: print(6)'  # prints 6.
exec '5'                  # does nothing and returns nothing.

def greet_me(**kwargs):
    if kwargs is not None:
        for key, value in kwargs.iteritems():
            print "%s == %s" %(key,value)
 
greet_me(name="abcd")

import site
print(site.getsitepackages())  # Return the path of the user base directory
print(site.getusersitepackages()) # Return the path of the user-specific site-packages directory, USER_SITE

from __future__ import division
import numpy as np
import numpy.linalg as la
import scipy.special as sp


class CompositeLegendreDiscretization:
    """A discrete function space on a 1D domain consisting of multiple
    subintervals, each of which is discretized as a polynomial space of
    maximum degree *order*. (see constructor arguments)

    There are :attr:`nintervals` * :attr:`npoints` degrees of freedom
    representing a function on the domain. On each subinterval, the
    function is represented by its function values at Gauss-Legendre
    nodes (see :func:`scipy.speical.legendre`) mapped into that
    subinterval.

    .. note::

        While the external interface of this discretization
        is exclusively in terms of vectors, it may be practical
        to internally reshape (see :func:`numpy.reshape`) these
        vectors into 2D arrays of shape *(nintervals, npoints)*.

    The object has the following attributes:

    .. attribute:: intervals

        The *intervals* constructor argument.

    .. attribute:: nintervals

        Number of subintervals in the discretization.

    .. attribute:: npoints

        Number of points on each interval. Equals *order+1*.

    .. attributes:: nodes

        A vector of ``(nintervals*npoints)`` node locations, consisting of
        Gauss-Legendre nodes that are linearly (or, to be technically correct,
        affinely) mapped into each subinterval.

    """


    def __init__(self, intervals, order):
        """
        :arg intervals: determines the boundaries of the subintervals to
            be used. If this is ``[a,b,c]``, then there are two subintervals
            :math:`(a,b)` and :math:`(b,c)`.
            (and the overall domain is :math:`(a,c)`)

        :arg order: highest polynomial degree being used
        """
        self.intervals=intervals

        self.nintervals=len(intervals)-1

        self.npoints=order + 1

        #Initializing shifted_nodes
        shifted_nodes=np.zeros(self.npoints*self.nintervals) 

        #Calling the scipy function to obtain the unshifted nodes
        unshifted_nodes=sp.legendre(self.npoints).weights[:,0]

        #Initializing shifted weights
        shifted_weights=np.zeros(self.nintervals*self.npoints)

        #Calling the scipy function to obtain the unshifted weights
        unshifted_weights=sp.legendre(self.npoints).weights[:,1]

        #Linearly mapping the unshifted nodes and weights to get the shifted
        #nodes and weights
        for i in range(self.nintervals):
                shifted_nodes[i*self.npoints:(i+1)*self.npoints]=(self.intervals[i]+ self.intervals[i+1])/2 + (self.intervals[i+1]-self.intervals[i])*(unshifted_nodes[0:self.npoints])/2
                shifted_weights[i*self.npoints:(i+1)*self.npoints]=(self.intervals[i+1]-self.intervals[i])*(unshifted_weights[0:self.npoints])/2
       
        #Setting nodes and weights attributes
        self.nodes=np.reshape(shifted_nodes,(self.nintervals,self.npoints))
        self.weights=np.reshape(shifted_weights,(self.nintervals,self.npoints))
    
        #Obtaining Vandermonde and RHS matrices to get A
        def vandermonde_rhs(m,arr):
            X=np.zeros((m,m))
            RHS=np.zeros((m,m))
            for i in range(m):
                for j in range(m):
                    X[i][j]=arr[i]**j
                    RHS[i][j]=((arr[i]**(j+1))-((-1)**(j+1)))/(j+1)
            return X,RHS
        
        A=np.zeros((self.npoints,self.npoints))
        X,RHS=vandermonde_rhs(self.npoints,unshifted_nodes)

        #Solving for spectral integration matrix
        A=np.dot(RHS,la.inv(X))
        
        self.A=A

        

    def integral(self, f):
        r"""Use Gauss-Legendre quadrature on each subinterval to approximate
        and return the value of

        .. math::

            \int_a^b f(x) dx

        where :math:`a` and :math:`b` are the left-hand and right-hand edges of
        the computational domain.

        :arg f: the function to be integrated, given as function values
            at :attr:`nodes`
        """
        int_val=0
        #Using Composite Gauss Quadrature to evaluate the whole integral

        for i in range(self.nintervals):
            f_val_arr=f[i,:]
            dot_val=np.dot(self.weights[i,:],f_val_arr.T)
            int_val+=dot_val

        return int_val




    def left_indefinite_integral(self, f):
        r"""Use a spectral integration matrix on each subinterval to
        approximate and return the value of

        .. math::

            g(x) = \int_a^x f(x) dx

        at :attr:`nodes`, where :math:`a` is the left-hand edge of the
        computational domain.

        The computational cost of this routine is linear in
        the number of degrees of freedom.

        :arg f: the function to be integrated, given as function values
            at :attr:`nodes`
        """
        #Initializing the left_indefinite_integral
        I=np.zeros((self.nintervals,self.npoints))
        current_quad_val=0

        for i in range(self.nintervals):
            
            if i>0:
                    f_val_prev=f[i-1:i,:]
                    weights_prev=self.weights[i-1:i,:]
                    current_quad_val+=np.dot(weights_prev,f_val_prev.T)
                
            current_f_val=f[i,:]

            #Scaling the Spectral integration matrix by the interval length
            A_scale=np.dot((self.intervals[i+1]-self.intervals[i])/2,self.A)

            spectral_factor=np.dot(A_scale,current_f_val.T)

            #Adding the  spectral factor to the evaluated definite integral
            indefinite_int=current_quad_val + spectral_factor
            I[i,:]=indefinite_int
            J=np.reshape(I,(self.nintervals*self.npoints,1))
        

        return I

    def right_indefinite_integral(self, f):
        r"""Use a spectral integration matrix on each subinterval to
        approximate and return the value of

        .. math::

            g(x) = \int_x^b f(x) dx

        at :attr:`nodes`, where :math:`b` is the left-hand edge of the
        computational domain.

        The computational cost of this routine is linear in
        the number of degrees of freedom.

        :arg f: the function to be integrated, given as function values
            at :attr:`nodes`
        """
        #Initializing the right integral
        I1=np.zeros(self.nintervals*self.npoints)

        #Unravelling the weights and function values at nodes
        f_right=f.ravel()    
        weights_right=self.weights.ravel()  
        

        for i in range(self.nintervals):
            
            #Evaluating Gauss Quadrature for all intervals from current to end
            f_next=f_right[i*self.npoints:self.npoints*self.nintervals]
            weights_next=weights_right[i*self.npoints:self.npoints*self.nintervals]
            current_quad_val=np.dot(weights_next,f_next.T)
            current_f_val=f_right[i*self.npoints:(i+1)*self.npoints]
            #Scaling the spectral integration matrix
            A_scale=np.dot((self.intervals[i+1]-self.intervals[i])/2,self.A)
            spectral_factor=np.dot(A_scale,current_f_val.T)
            
            #Subtracting spectral factor from Gauss Quadrature value
            I1[(i)*self.npoints:(i+1)*self.npoints]=np.dot(weights_next,f_next.T) - spectral_factor
            
            J=np.reshape(I1,(self.nintervals,self.npoints))

    
        return J
        


from __future__ import division
#from legendre_discr import CompositeLegendreDiscretization
import numpy as np
import matplotlib.pyplot as pt


def get_left_int_error(n, order):
    a = 2
    b = 30
    intervals = np.linspace(0, 1, n, endpoint=True) ** 2 * (b-a) + a
    discr = CompositeLegendreDiscretization(intervals, order)

    x = discr.nodes
 
    assert abs(discr.integral(1+0*x) - (b-a)) < 1e-13


    alpha = 4
    from scipy.special import jv, jvp
    f = jvp(alpha, x)

    num_int_f = jv(alpha, a) + discr.left_indefinite_integral(f)
    int_f = jv(alpha, x)

    if 0:
        pt.plot(x.ravel(), num_int_f.ravel())
        pt.plot(x.ravel(), int_f.ravel())
        pt.show()

    L2_err = np.sqrt(discr.integral((num_int_f - int_f)**2))
    return 1/n, L2_err


def get_right_int_error(n, order):
    a = 2
    b = 30
    intervals = np.linspace(0, 1, n, endpoint=True) ** 2 * (b-a) + a
    discr = CompositeLegendreDiscretization(intervals, order)

    x = discr.nodes

    assert abs(discr.integral(1+0*x) - (b-a)) < 1e-13

    alpha = 4
    from scipy.special import jv, jvp
    f = jvp(alpha, x)

    num_int_f = jv(alpha, b) - discr.right_indefinite_integral(f)
    int_f = jv(alpha, x)

    if 0:
        pt.plot(x.ravel(), num_int_f.ravel())
        pt.plot(x.ravel(), int_f.ravel())
        pt.show()
    L2_err = np.sqrt(discr.integral((num_int_f - int_f)**2))
    return 1/n, L2_err


def estimate_order(f, point_counts):
    n1, n2 = point_counts
    h1, err1 = f(n1)
    h2, err2 = f(n2)

    print "h=%g err=%g" % (h1, err1)
    print "h=%g err=%g" % (h2, err2)

    from math import log
    est_order = (log(err2/err1) / log(h2/h1))
    print "%s: EOC: %g" % (f.__name__, est_order)
    print

    return est_order


if __name__ == "__main__":
    for order in [2, 3, 5, 7]:
        print "---------------------------------"
        print "ORDER", order
        print "---------------------------------"
        assert (estimate_order(lambda n: get_left_int_error(n, order), [10, 30])
                >= order-0.5)
        assert (estimate_order(lambda n: get_right_int_error(n, order), [10, 30])
                >= order-0.5)