import plainEnglishCoding
import warnings
warnings.filterwarnings("ignore")
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import ipywidgets as widgets

M = np.array(
    [
        [3,6,2],
        [10,3,8],
        [1,7,5],
    ]
)
costs = np.array([2,3,5])
totals = np.dot(M,costs)
eigvals,eigvec = np.linalg.eig(M)
print(eigvals)
print(totals)

[15.16204133  2.41626789 -6.57830922]
[34 69 48]

M = np.array(
    [
        [3,6,2],
        [10,3,8],
        [1,7,5],
    ]
)
p = np.array([2,3,5])
np.dot(M,p)

array([34, 69, 48])

M = np.array(
    [
        [3,6,2],
        [10,3,8],
        [1,7,5],
    ]
)
M1 = np.array(
    [
        [41/241,16/241,-42/241],
        [42/241,-13/241,4/241],
        [-67/241,15/241,51/241],
    ]
)
c = np.array([34,69,48])
np.dot(M1,c)

array([2., 3., 5.])

# np.dot(M1,M)
np.round(np.dot(M1,M),decimals=10)

array([[ 1., -0.,  0.],
       [-0.,  1., -0.],
       [-0., -0.,  1.]])

# using symbolic programming to check my work throughout
from sympy import *
M = Matrix(
    [
        [3,6,2,1,0,0],
        [10,3,8,0,1,0],
        [1,7,5,0,0,1]
    ])
M[0,:] /= 3
M[1,:] = M[1,:]-10*M[0,:]
M[2,:] = M[2,:]-1*M[0,:]
M[1,:] /= -17
M[2,:] = M[2,:]-5*M[1,:]
M[2,:] /= M[2,2]
M[1,:] -= M[1,2]*M[2,:]
M[0,:] -= M[0,2]*M[2,:]
M[0,:] -= M[0,1]*M[1,:]
# M
print(latex(M))
print(M[:,3:]*Matrix([34,69,48]))
M[:,3:]*Matrix([[3,6,2],[10,3,8],[1,7,5]])

\left[\begin{matrix}1 & 0 & 0 & \frac{41}{241} & \frac{16}{241} & - \frac{42}{241}\\0 & 1 & 0 & \frac{42}{241} & - \frac{13}{241} & \frac{4}{241}\\0 & 0 & 1 & - \frac{67}{241} & \frac{15}{241} & \frac{51}{241}\end{matrix}\right]
Matrix([[2], [3], [5]])

# using symbolic programming to check check my work throughout
from sympy import *
M = Matrix(
    [
        [3,6,2,1,0,0],
        [10,3,8,0,1,0],
        [20,6,16,0,0,1]
    ])
M[0,:] /= M[0,0]
M[1,:] = M[1,:]-M[1,0]*M[0,:]
M[2,:] = M[2,:]-M[2,0]*M[0,:]
# M[3,:] = M[3,:]-M[3,0]*M[0,:]
M[1,:] /= M[1,1]
M[2,:] = M[2,:]-M[2,1]*M[1,:]
print(latex(M))
M
# print(M[:,3:]*Matrix([34,69,48]))
# M[:,3:]*Matrix([[3,6,2],[10,3,8],[1,7,5]])

\left[\begin{matrix}1 & 2 & \frac{2}{3} & \frac{1}{3} & 0 & 0\\0 & 1 & - \frac{4}{51} & \frac{10}{51} & - \frac{1}{17} & 0\\0 & 0 & 0 & 0 & -2 & 1\end{matrix}\right]

M = Matrix(
    [
        [3,6,2,1,0,0],
        [10,3,8,0,1,0],
        [20,6,16,0,0,1]
    ])
M = M[:,:3]
M.det()

M = Matrix(
    [
        [3,6,2,1,0,0],
        [10,3,8,0,1,0],
        [20,6,16,0,0,1],
        [1,7,5,0,0,1],
    ])
print(latex(M))
M[0,:] /= M[0,0]
M[1,:] = M[1,:]-M[1,0]*M[0,:]
M[2,:] = M[2,:]-M[2,0]*M[0,:]
M[3,:] = M[3,:]-M[3,0]*M[0,:]
M[1,:] /= M[1,1]
M[2,:] = M[2,:]-M[2,1]*M[1,:]
M[3,:] -= M[3,1]*M[1,:]
M[3,:] /= M[3,2]
# M[2,:] /= M[2,2]
# M[1,:] -= M[1,2]*M[2,:]
# M[0,:] -= M[0,2]*M[2,:]
# M[0,:] -= M[0,1]*M[1,:]
print(latex(M))
M

\left[\begin{matrix}3 & 6 & 2 & 1 & 0 & 0\\10 & 3 & 8 & 0 & 1 & 0\\20 & 6 & 16 & 0 & 0 & 1\\1 & 7 & 5 & 0 & 0 & 1\end{matrix}\right]
\left[\begin{matrix}1 & 2 & \frac{2}{3} & \frac{1}{3} & 0 & 0\\0 & 1 & - \frac{4}{51} & \frac{10}{51} & - \frac{1}{17} & 0\\0 & 0 & 0 & 0 & -2 & 1\\0 & 0 & 1 & - \frac{67}{241} & \frac{15}{241} & \frac{51}{241}\end{matrix}\right]

# problem setup
# matrix of individuals' purchasing quantities
M = np.array(
    [
        [3,6,2],
        [10,3,8],
        [1,7,5],
    ]
)
# prices per item
prices = np.array([2,3,5])
# calculate total costs per individual
totalCosts = np.dot(M,prices)
totalCosts

array([34, 69, 48])

# using Numpy to invert a matrix
Minv = np.linalg.inv(M)
print(Minv)

# test that Minv*M yields identity matrix
# note '@' is an alias for np.dot
# so A@B is equivalent to np.dot(A,B)
print(Minv@M)

# same thing but round resulting matrix components
print(np.round(Minv@M,decimals=10)) 

# use Minv to solve for grocery prices
prices2 = np.dot(Minv,totalCosts)
prices2

[[ 0.17012448  0.06639004 -0.17427386]
 [ 0.17427386 -0.05394191  0.01659751]
 [-0.2780083   0.06224066  0.21161826]]
[[ 1.00000000e+00  1.66533454e-16  5.55111512e-17]
 [-2.39391840e-16  1.00000000e+00 -1.14491749e-16]
 [ 1.66533454e-16 -2.77555756e-16  1.00000000e+00]]
[[ 1.  0.  0.]
 [-0.  1. -0.]
 [ 0. -0.  1.]]

array([2., 3., 5.])

# another option is the Numpy solve() function to 
# treat the same problem in it's original form.
# Behind this scenes, this approach is still inverting M.
prices3 = np.linalg.solve(M,totalCosts)
prices3

array([2., 3., 5.])

# Numpy does calculations on binary objects
a,b = 10,2
print("a+b = %d" % int(a+b))

# equivalent calculation for simple unsigned bit strings
import bitstring
bit_a = bitstring.Bits(uint=a,length=16)
bit_b = bitstring.Bits(uint=b,length=16)
print("bit_a: %s" % bit_a.b)
print("bit_b: %s" % bit_b.b)

def addBitStrings(b1,b2):
    # function for adding two binary strings
    # using bitwise operations.
    # This mimics how CPUs perform calculations.
    carry = 0
    Z = bitstring.Bits(uint=0,length=16)
    while b2 != Z:
        carry = b1 & b2
        b1 = b1 ^ b2
        b2 = carry << 1
    return b1

bit_c = addBitStrings(bit_a,bit_b)
print("bit_c: %s" % bit_c.b)
print("bit_c as an integer: %d" % bit_c.i)

a+b = 12
bit_a: 0000000000001010
bit_b: 0000000000000010
bit_c: 0000000000001100
bit_c as an integer: 12

# normal calculations of fractions results in decimals
a = 1/2
b = 2/3
print(a+b)

# import SymPy
import sympy,fractions
# define symbols instead of variables
x = fractions.Fraction(1,2)
y = fractions.Fraction(2,3)
e = x+y
print(e)
print(7/6)

1.1666666666666665
7/6
1.1666666666666667

import pprint
def symbolicGaussianElimination(M):
    n = M.shape[0]
    aM = M.row_join(sympy.eye(n))
    # From top to bottom
    print("Going down the rows!")
    for i in range(n):
        # normalize the ith row to create 1 in the (i,i) position
        aM[i,:] /= aM[i,i]
        for j in range(i+1,n):
            # remove multiple of row i from all subsequent rows
            aM[j,:] -= aM[j,i]*aM[i,:]
        #print current step
        pprint.pp(aM)
    print("Going back up the rows!")
    for i in range(n,0,-1):
        i -= 1
        for j in range(0,i):
            # for row above the ith row, subtract multiple of row i
            aM[j,:] -= aM[j,i]*aM[i,:]
        #print current step
        pprint.pp(aM)
    return aM
        
        
M = sympy.Matrix(
    [
        [3,6,2],
        [10,3,8],
        [1,7,5],
    ]
)
symbolicGaussianElimination(M);

Going down the rows!
Matrix([
[1,   2,  2/3,   1/3, 0, 0],
[0, -17,  4/3, -10/3, 1, 0],
[0,   5, 13/3,  -1/3, 0, 1]])
Matrix([
[1, 2,    2/3,    1/3,     0, 0],
[0, 1,  -4/51,  10/51, -1/17, 0],
[0, 0, 241/51, -67/51,  5/17, 1]])
Matrix([
[1, 2,   2/3,     1/3,      0,      0],
[0, 1, -4/51,   10/51,  -1/17,      0],
[0, 0,     1, -67/241, 15/241, 51/241]])
Going back up the rows!
Matrix([
[1, 2, 0, 125/241, -10/241, -34/241],
[0, 1, 0,  42/241, -13/241,   4/241],
[0, 0, 1, -67/241,  15/241,  51/241]])
Matrix([
[1, 0, 0,  41/241,  16/241, -42/241],
[0, 1, 0,  42/241, -13/241,   4/241],
[0, 0, 1, -67/241,  15/241,  51/241]])
Matrix([
[1, 0, 0,  41/241,  16/241, -42/241],
[0, 1, 0,  42/241, -13/241,   4/241],
[0, 0, 1, -67/241,  15/241,  51/241]])

# Use Sympy to represent a symbolic matrix
m11,m12,m21,m22 = sympy.symbols("m11 m12 m21 m22")
M = sympy.Matrix(
    [
        [m11,m12],
        [m21,m22]
    ]
)
M

# use Sympy to symbolically invert 2x2 matrix
M1 = M.inv()
M1

f = sympy.lambdify([m11,m12,m21,m22],M1)
f(5,7,2,3)

array([[ 3., -7.],
       [-2.,  5.]])

# the symbolic solution is nice enough for 2x2 matrices, but what about 3x3?
m11,m12,m13,m21,m22,m23,m31,m32,m33 = sympy.symbols("m11 m12 m13 m21 m22 m23 m31 m32 m33")
M = sympy.Matrix(
    [
        [m11,m12,m13],
        [m21,m22,m23],
        [m31,m32,m33]
    ]
)
M

# craziness
M.inv()

# using code from https://github.com/adamdempsey90/gravitational-lensing/tree/master
from matplotlib import cm
import matplotlib.pyplot as plt
# from scipy.misc import imread
from ipywidgets import interact,fixed,Checkbox,Text
import IPython.display as disp
import dm_funcs as dm

global x,y,galaxy_properties,gal
@interact(n=(10,400,50),sigma=(.01,.5,.01),xcen=(-1,1,.1),ycen=(-1,1,.1),axratio=(.2,4,.1),posangle=(0,180,5))
def set_number_of_points(n=100,sigma=0.2,xcen=0,ycen=0,axratio=2.0,posangle=30):
    """Set up a mock galaxy with the input parameters.
        n:              sets the number of points in each direction.
        sigma:          sets the 'width' of the galaxy.
        xcen & ycen:    set the center of the galaxy.
        axratio:        sets the ellipticity of the galaxy.
        posangle:       sets the orientation of the galaxy. """
    
    global x,y,galaxy_properties,gal
    (x,y) = dm.load_grid(n)
    galaxy_properties={'peak':1,'sigma':sigma,'xcen':xcen,'ycen':ycen,'axratio':axratio,'posangle':posangle}
    gal = dm.create_galaxy(x,y,**galaxy_properties)
    plt.imshow(gal);

interactive(children=(IntSlider(value=100, description='n', max=400, min=10, step=50), FloatSlider(value=0.2, …

def compute_deflection2(x,y,**kwargs):
        """ Compute the deflection and shear from a singular isothermal ellipsoid potential """
        mass = kwargs['Mass']
        xcen = kwargs['xcen']
        ycen = kwargs['ycen']
        axratio = np.abs(kwargs['axratio'])
        posangle = kwargs['posangle']

        tol = .001

        if axratio>1:
                axratio = 1.0/axratio
                posangle += 90
        phi = np.deg2rad(posangle)
        

        def rotate(x,y,phi,clockwise=False):
            # rotation matrix 
            if clockwise:
                R = np.array(
                    [
                        [np.cos(phi),np.sin(phi)],
                        [-np.sin(phi),np.cos(phi)]
                    ]
                )
            else:
                R = np.array(
                    [
                        [np.cos(phi),np.sin(phi)],
                        [-np.sin(phi),np.cos(phi)]
                    ]
                )
            P = np.stack((x.ravel(),y.ravel()))
            P = R@P
            x_new = P[0,:].reshape(x.shape)
            y_new = P[1,:].reshape(y.shape)
            return x_new,y_new
            
        x_new,y_new = rotate(x-xcen,y-ycen,phi,clockwise=True)

        r = np.sqrt(axratio*x_new*x_new + y_new * y_new /axratio)

        q = np.sqrt(1./axratio - axratio)
        softening = (r == 0).astype(float)
        fac_x = x_new/(r + softening)
        fac_y = y_new/(r + softening)



        if q>=tol:
                x_def = (mass/q)*np.arctan(q*fac_x)
                y_def = (mass/q)*np.arctanh(q*fac_y)
        else:
                x_def = mass*fac_x
                y_def = mass*fac_y

        x_out,y_out = rotate(x_def,y_def,phi,clockwise=False)

        return x_out,y_out

@interact(
    DarkMatterSize=(0,3,.1), 
    DarkMatterXCenter=(-1,1,.01),
    DarkMatterYCenter=(-1,1,.01),
    DarkMatterAxisRatio=(.01,4,.03),
    DarkMatterPositionAngle=(0,180,5)
)
def plot_galaxy(DarkMatterSize=.1,DarkMatterXCenter=0,DarkMatterYCenter=0,DarkMatterAxisRatio=1,DarkMatterPositionAngle=0):
        global x,y,galaxy_properties,gal
       
        plt.close();
#         xl,yl = dm.compute_deflection(x,y,Mass=DarkMatterSize,xcen=DarkMatterXCenter, 
#                                       ycen=DarkMatterYCenter,axratio=DarkMatterAxisRatio,
#                                       posangle=DarkMatterPositionAngle)
        xl,yl = compute_deflection2(x,y,Mass=DarkMatterSize,xcen=DarkMatterXCenter, 
                                      ycen=DarkMatterYCenter,axratio=DarkMatterAxisRatio,
                                      posangle=DarkMatterPositionAngle)
        gal = dm.create_galaxy(x-xl,y-yl,**galaxy_properties)
 
        plt.imshow(gal)

interactive(children=(FloatSlider(value=0.1, description='DarkMatterSize', max=3.0), FloatSlider(value=0.0, de…

Systems of Linear Equations¶

Backsolving and Inverting Matrices¶

Perils in Backsolving¶

Python Programming and Inverting Matrices¶

Backsolving Example: Gravitational Lensing¶