notebook.community

Edit and run

Module 1 Python





In [1]:



    


# Import Packages
from SimPEG import *
from IPython.html.widgets import interactive
import numpy as np
%pylab inline



    


















    






Efficiency Warning: Interpolation will be slow, use setup.py!

            python setup.py build_ext --inplace
    
Populating the interactive namespace from numpy and matplotlib










    






WARNING: pylab import has clobbered these variables: ['interactive']
`%matplotlib` prevents importing * from pylab and numpy

















In [2]:



    


# Begin by creating a ficticious set of model data

N=1000                      # Set number of model parameters 
x=np.linspace(0,1,N+1)      # Define 1D domain on nodes
xc = 0.5*(x[1:] + x[0:-1])  # Define 1D domain on cell centers
m = np.zeros(N)             # preallocate m array 

# Define Gaussian function:
gauss = lambda x, a, mean, SD: a*np.exp(-(x-mean)**2./(2.*SD**2.)) # create a Gaussian function

# Choose parameters for Gaussian, evaluate, and store in an array, f.
SD=6                           
mean=0
a=1/(SD*sqrt(2*pi))
x2=np.linspace(-20,20,0.2*N)  
f = gauss(x2,a, mean,SD)

# Define a boxcar function:
box = 0.4*np.ones(0.2*N)

# Insert the Gaussian and Boxcar into m:
m[0.2*N:0.4*N]=box
m[0.6*N:0.8*N]=10*f
   
# Plot      
pylab.plot(xc,m)
#pylab.plot(f)
print size(f)
pylab.xlabel('x')
pylab.ylabel('m(x)')
pylab.title('Model, $m(x)$')
print shape(xc)



    


















    






200
(1000L,)










    
























In [3]:



    


# Make the set of kernel functions

g = lambda x, i, p, q: np.exp(-p*i*x)*np.cos(2*np.pi*q*i*x) # create an anonymous function as immediately above
p = 0.01      # Set values for p, q
q = 0.15
M = 20        # specify number of output data
Gn = np.zeros((N+1,M))

for i in range(M):
    f = g(x,i,p,q)
    Gn[:,i] = f
    
# Plot  
pylab.plot(x,Gn)
pylab.xlabel('x')
pylab.ylabel('g(x)')
pylab.title('Kernel functions, $g(x)$')



    


















    
Out[3]:








<matplotlib.text.Text at 0xa58cbe0>










    
























In [4]:



    


# Make Averaging Matrix
n=N+1            # Define n as the N+1 dimension of the matrix   
w=n-1            # Define w and the n-1  dimentsion of the matix   
s = (N,n)        # Store matrix values 
Av = np.zeros(s) # Create a matrix of zeros of the correct dimensions 
                 # and fill in with elements usin the loop below (note the 1/2 is included in here). 
for i in range(N):
    j=i
    k=i+1
    Av[i,j] = 0.5  
    Av[i,k] = 0.5
print shape(Av)

# Plot  
pylab.title('Subset of averaging matrix')
pyplot.imshow(Av[:5,:6], cmap='Greys',  interpolation='nearest')
print(Av[:5,:6])



    


















    






(1000L, 1001L)
[[ 0.5  0.5  0.   0.   0.   0. ]
 [ 0.   0.5  0.5  0.   0.   0. ]
 [ 0.   0.   0.5  0.5  0.   0. ]
 [ 0.   0.   0.   0.5  0.5  0. ]
 [ 0.   0.   0.   0.   0.5  0.5]]










    
























In [5]:



    


# make the Volume, "delta x" array

Deltax = diff(x)
V = np.diag(Deltax)     # create diagonal matrix     
print shape(V)
print np.sum(Deltax)
print shape(x)



    


















    






(1000L, 1000L)
1.0
(1001L,)

















In [6]:



    


G=np.dot(np.transpose(np.dot(Av, Gn)), V)
d = np.dot(np.dot(np.transpose(np.dot(Av, Gn)), V),m)   # Create synthetic data
xd=np.linspace(0,1,len(d))                              # make x array for data 

# Plot

pylab.plot(xd,d)
pylab.xlabel('')
pylab.ylabel('d')
pylab.title('Synthetic Data $d$')



    


















    
Out[6]:








<matplotlib.text.Text at 0xa247898>

Module 2 Python





In [7]:



    


# Add noise to our synthetic data

noise = lambda fl, length, data, s: fl + randn(length)*data*s # introduce noise using a floor (f), length (l) and scaling factor (s)

r = noise(0, len(d), d, 0.04)

dobs= d + r

#pylab.plot(xd,r)
pylab.plot(xd,d)
pylab.plot(xd,dobs)
pylab.xlabel('')
pylab.ylabel('d')
pylab.title('Synthetic Data and noise, $d$')



    


















    
Out[7]:








<matplotlib.text.Text at 0xa537358>










    
























In [8]:



    


# Calculate the data misfit, phi_d #
####################################
invvar = lambda floor, percent, data: 1./(floor + percent*np.abs(data))

L2 = lambda M, w: dot(dot(w.T,M.T),dot(M,w))                                # Anonymous function for 2-Norm

# assign a floor
flr = 0.015*dot(d.T,d)**0.5/len(d)

# Make Wd
eps = invvar(flr,0.02,dobs) # define epsilon and Wd
Wd=np.diag(eps) 

# Take the 2-norm
phi_d = L2(Wd, np.dot(G,m)-dobs)
print phi_d
print flr



    


















    






25.2963685434
0.000200320925096

















In [9]:



    


# Reference Model#
##################
## A constant reference model ##
mref = np.mean(m)*np.ones(N)   

# Plot      
pylab.plot(xc,mref)
pylab.xlabel('x')
pylab.ylabel('m(x)')
pylab.title('Reference Model, $m(x)$')



    


















    
Out[9]:








<matplotlib.text.Text at 0xa80c550>










    
























In [10]:



    


# Domains #
###############################                       # Set up domains:
l = np.power(x[1:len(x)]-x[0:len(x)-1],0.5)           # vector of distances between nodes, l
midx=np.dot(Av,x)                                     # vector of midpoints, midx 
h = np.power(midx[1:len(midx)]-x[0:len(midx)-1], 0.5) # Calculate distances between midpoints & take sqr roots of each value in vector h



    











In [11]:



    


#  Calculate the model regularization, phi_m #
##############################################

# Create Ws, Wx matrices
Ws=np.diag(l)                                  # put length vector into a diagonal matrix, Ws
Wx = np.zeros((len(m), len(m)))                # preaollocate array and enter values into matrix, Wx
for i in range(len(h)):
    j=i
    k=i+1
    Wx[i,j] = h[i]  
    Wx[i,k] = h[i]  
    
alpha_s=0.1                                    # Set alpha_s and alpha_x values  
alpha_x=0.15

# build Wm #
Wm=np.concatenate((alpha_s*Ws, alpha_x*Wx), axis=0)

beta = 1000000
# Set beta value

# Get recovered model
mrecovered = lambda G,Wd,Wm,data,mref, beta: dot(np.linalg.inv(dot(dot(G.T,Wd.T),dot(Wd,G)) + beta*dot(Wm.T,Wm)),dot(dot(G.T,Wd),dot(Wd,dobs)) + beta*dot(dot(Wm.T,Wm),mref))
# mrec = dot(np.linalg.inv(dot(dot(G.T,Wd.T),dot(Wd,G)) + beta*dot(Wm.T,Wm)),dot(dot(G.T,Wd),dot(Wd,dobs)) + beta*dot(dot(Wm.T,Wm),mref))

mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)

# Get residual of mref and mrec
phi_d2 = L2(Wd,np.dot(G,mrec)-dobs)
phi_m2 = L2(Wm, mref-mrec)


print phi_d2
print phi_m2

# Plot      
pylab.plot(xc,m)
pylab.plot(xc,mrec)
pylab.xlabel('x')
pylab.ylabel('m(x)')
pylab.title('Model, $m(x)$')



    


















    






560.724989307
0.00266347462601










    
Out[11]:








<matplotlib.text.Text at 0xa831e10>










    
























In [58]:



    


# Get predicted data

dpred = np.dot(G,mrec)

#pylab.plot(xd,r)
pylab.plot(xd,dpred,'o')
pylab.plot(xd,dobs)
pylab.plot(xd,d,'x')
pylab.xlabel('')
pylab.ylabel('d')
pylab.title('Predicted and Observed Data')



    


















    
Out[58]:








<matplotlib.text.Text at 0x15eca9b0>

Beta cooling

The $\beta$ Cooling method

mar 25





In [12]:



    


print phi_d



    


















    






25.2963685434

















In [13]:



    


# Make the Tikhonov Curve
beta = 100000                                  # Set beta value
itns = 10                                      # set maximum number of iterations
phdTik=np.zeros(itns)                          # preallocate arrays
phmTik=np.zeros(itns)
foo = np.zeros(itns)
b=np.zeros(itns)


for i in range(itns):
    mrec = dot(np.linalg.inv(dot(dot(G.T,Wd.T),dot(Wd,G)) + beta*dot(Wm.T,Wm)),dot(dot(G.T,Wd),dot(Wd,dobs)) + beta*dot(dot(Wm.T,Wm),mref))
    b[i] = beta
    phdTik[i] = L2(Wd, np.dot(G,mrec)-dobs)
    phmTik[i] = L2(Wm, mref-mrec)
    foo[i] = np.absolute(phdTik[i]-phdstar)
    beta=beta*0.5


pylab.plot(phmTik,phdTik,'o')
pylab.xlabel('$\phi_m$')
pylab.ylabel('$\phi_d$')
pylab.title('Tikhonov Curve')



    


















    
Out[13]:








<matplotlib.text.Text at 0xadbc160>

This cell finds beta





In [14]:



    


beta_0=1000000
beta=beta_0
tol=0.00001
i=1
itns=10
mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)
phid = L2(Wd,np.dot(G,mrec)-dobs)
phi_dstar=20

while np.absolute(phid - phi_dstar)>tol:  
    j=i
    mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)
    phid = L2(Wd,np.dot(G,mrec)-dobs)
    print("Iteration = ",j, "beta=", beta, "phid at start=",phid)
    if phid>phi_dstar:
        beta = beta_0*0.5**i
        print("beta=", beta)
    if phid<phi_dstar:
        beta_abv = beta_0*0.5**(i-2)
        beta_bel = beta_0*0.5**(i-1)
        beta_mid = 0.5*(beta_abv+beta_bel)
        mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
        phidmid = L2(Wd,np.dot(G,mmid)-dobs)  
        #for k in range(itns):
        while np.absolute(phidmid - phi_dstar)>tol:
            if phidmid<phi_dstar:
                beta_abv = beta_0*0.5**(i-2)
                beta_bel = beta_mid
                beta_mid = 0.5*(beta_abv+beta_bel)
                beta=beta_mid
                mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
                phidmid = L2(Wd,np.dot(G,mmid)-dobs)           
            if phidmid>phi_dstar:
                beta_abv = beta_mid
                beta_bel = beta_bel
                beta_mid = 0.5*(beta_abv+beta_bel)
                beta=beta_mid
                mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
                phidmid = L2(Wd,np.dot(G,mmid)-dobs)                
                if phidmid<phi_dstar:
                    beta_abv = beta_abv
                    beta_bel = beta_mid
                    beta_mid = 0.5*(beta_abv+beta_bel)
                    beta=beta_mid
                mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
                phidmid = L2(Wd,np.dot(G,mmid)-dobs)              
            phid=phidmid
            mrec=mmid
    i=i+1   
print("At end beta=", beta,"phid=", phid)



    


















    






('Iteration = ', 1, 'beta=', 1000000, 'phid at start=', 560.72498930749271)
('beta=', 500000.0)
('Iteration = ', 2, 'beta=', 500000.0, 'phid at start=', 190.13392440741501)
('beta=', 250000.0)
('Iteration = ', 3, 'beta=', 250000.0, 'phid at start=', 69.672931123447967)
('beta=', 125000.0)
('Iteration = ', 4, 'beta=', 125000.0, 'phid at start=', 33.528159668743008)
('beta=', 62500.0)
('Iteration = ', 5, 'beta=', 62500.0, 'phid at start=', 22.443001524284853)
('beta=', 31250.0)
('Iteration = ', 6, 'beta=', 31250.0, 'phid at start=', 18.622133840371461)
('At end beta=', 43639.77606408298, 'phid=', 20.000008521815857)

This is the recovered model for the given beta





In [1]:



    


mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)
phid = L2(Wd,np.dot(G,mrec)-dobs)
print (beta)
print phid
# Plot      
pylab.plot(xc,m)
pylab.plot(xc,mrec)
pylab.xlabel('x')
pylab.ylabel('m(x)')
pylab.title('Model, $m(x)$')



    


















    






---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-ac0c8d279f2d> in <module>()
----> 1 mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)
      2 phid = L2(Wd,np.dot(G,mrec)-dobs)
      3 print (beta)
      4 print phid
      5 # Plot

NameError: name 'mrecovered' is not defined

This is a version that shows what is going on inside every iteration





In [65]:



    


beta_0=100000 
beta=beta_0
tol=0.1
i=1
itns=10
mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)
phid = L2(Wd,np.dot(G,mrec)-dobs)
phi_dstar=20

while np.absolute(phid - phi_dstar)>tol:  
    j=i
    mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)
    phid = L2(Wd,np.dot(G,mrec)-dobs)
    print("Iteration = ",j, "beta=", beta, "phid at start=",phid)
    if phid>phi_dstar:
        beta = beta_0*0.5**i
        print("beta=", beta)
    if phid<phi_dstar:
        beta_abv = beta_0*0.5**(i-2)
        beta_bel = beta_0*0.5**(i-1)
        beta_mid = 0.5*(beta_abv+beta_bel)
        mabv = mrecovered(G,Wd,Wm,dobs,mref,beta_abv)
        mbel = mrecovered(G,Wd,Wm,dobs,mref,beta_bel)
        mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
        phidabv = L2(Wd,np.dot(G,mabv)-dobs)
        phidbel = L2(Wd,np.dot(G,mbel)-dobs)
        phidmid = L2(Wd,np.dot(G,mmid)-dobs)  
        print("Itn # is:",j, "below threshold", "beta above=", beta_abv, "beta below=", beta_bel, "beta mid=",beta_mid)
        print("Itn # is:",j, "below threshold", "phi above=",  phidabv, "phi below=",  phidbel, "phi mid=", phidmid) 
        for k in range(itns):
        #while np.absolute(phidmid - phi_dstar)>tol:
            if phidmid<phi_dstar:
                beta_abv = beta_0*0.5**(i-2)
                beta_bel = beta_mid
                beta_mid = 0.5*(beta_abv+beta_bel)
                beta=beta_mid
                mabv = mrecovered(G,Wd,Wm,dobs,mref,beta_abv)
                mbel = mrecovered(G,Wd,Wm,dobs,mref,beta_bel)
                mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
                phidabv = L2(Wd,np.dot(G,mabv)-dobs)
                phidbel = L2(Wd,np.dot(G,mbel)-dobs)
                phidmid = L2(Wd,np.dot(G,mmid)-dobs) 
                print("Itn # is:",j, "in R1", "beta above=", beta_abv, "beta below=", beta_bel, "beta mid=",beta_mid)
                print("Itn # is:",j, "in R1", "phi above=",  phidabv, "phi below=",  phidbel, "phi mid=", phidmid)                  
            if phidmid>phi_dstar:
                beta_abv = beta_mid
                beta_bel = beta_bel
                beta_mid = 0.5*(beta_abv+beta_bel)
                beta=beta_mid
                mabv = mrecovered(G,Wd,Wm,dobs,mref,beta_abv)
                mbel = mrecovered(G,Wd,Wm,dobs,mref,beta_bel)
                mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
                phidabv = L2(Wd,np.dot(G,mabv)-dobs)
                phidbel = L2(Wd,np.dot(G,mbel)-dobs)
                phidmid = L2(Wd,np.dot(G,mmid)-dobs)  
                print("Itn # is:",j, "in R2", "beta above=", beta_abv, "beta below=", beta_bel, "beta mid=",beta_mid)
                print("Itn # is:",j, "in R2", "phi above=",  phidabv, "phi below=",  phidbel, "phi mid=", phidmid)                
                if phidmid<phi_dstar:
                    beta_abv = beta_abv
                    beta_bel = beta_mid
                    beta_mid = 0.5*(beta_abv+beta_bel)
                    beta=beta_mid
                mabv = mrecovered(G,Wd,Wm,dobs,mref,beta_abv)
                mbel = mrecovered(G,Wd,Wm,dobs,mref,beta_bel)
                mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
                phidabv = L2(Wd,np.dot(G,mabv)-dobs)
                phidbel = L2(Wd,np.dot(G,mbel)-dobs)
                phidmid = L2(Wd,np.dot(G,mmid)-dobs)
                print("Itn # is:",j, "in R3", "beta above=", beta_abv, "beta below=", beta_bel, "beta mid=",beta_mid)
                print("Itn # is:",j, "in R3", "phi above=",  phidabv, "phi below=",  phidbel, "phi mid=", phidmid)                  
            phid=phidmid
            mrec=mmid
    i=i+1   
print("At end beta=", beta,"phid=", phid)



    


















    






('Iteration = ', 1, 'beta=', 100000, 'phid at start=', 30.796129097570208)
('beta=', 50000.0)
('Iteration = ', 2, 'beta=', 50000.0, 'phid at start=', 22.240411728510566)
('beta=', 25000.0)
('Iteration = ', 3, 'beta=', 25000.0, 'phid at start=', 19.201693409779164)
('Itn # is:', 3, 'below threshold', 'beta above=', 50000.0, 'beta below=', 25000.0, 'beta mid=', 37500.0)
('Itn # is:', 3, 'below threshold', 'phi above=', 22.240411728510566, 'phi below=', 19.201693409779164, 'phi mid=', 20.610353891506911)
('Itn # is:', 3, 'in R2', 'beta above=', 37500.0, 'beta below=', 25000.0, 'beta mid=', 31250.0)
('Itn # is:', 3, 'in R2', 'phi above=', 20.610353891506911, 'phi below=', 19.201693409779164, 'phi mid=', 19.875865209693295)
('Itn # is:', 3, 'in R3', 'beta above=', 37500.0, 'beta below=', 31250.0, 'beta mid=', 34375.0)
('Itn # is:', 3, 'in R3', 'phi above=', 20.610353891506911, 'phi below=', 19.875865209693295, 'phi mid=', 20.236017406995565)
('Itn # is:', 3, 'in R2', 'beta above=', 34375.0, 'beta below=', 31250.0, 'beta mid=', 32812.5)
('Itn # is:', 3, 'in R2', 'phi above=', 20.236017406995565, 'phi below=', 19.875865209693295, 'phi mid=', 20.05412390871237)
('Itn # is:', 3, 'in R3', 'beta above=', 34375.0, 'beta below=', 31250.0, 'beta mid=', 32812.5)
('Itn # is:', 3, 'in R3', 'phi above=', 20.236017406995565, 'phi below=', 19.875865209693295, 'phi mid=', 20.05412390871237)
('Itn # is:', 3, 'in R2', 'beta above=', 32812.5, 'beta below=', 31250.0, 'beta mid=', 32031.25)
('Itn # is:', 3, 'in R2', 'phi above=', 20.05412390871237, 'phi below=', 19.875865209693295, 'phi mid=', 19.964533792945403)
('Itn # is:', 3, 'in R3', 'beta above=', 32812.5, 'beta below=', 32031.25, 'beta mid=', 32421.875)
('Itn # is:', 3, 'in R3', 'phi above=', 20.05412390871237, 'phi below=', 19.964533792945403, 'phi mid=', 20.009214504197004)
('Itn # is:', 3, 'in R2', 'beta above=', 32421.875, 'beta below=', 32031.25, 'beta mid=', 32226.5625)
('Itn # is:', 3, 'in R2', 'phi above=', 20.009214504197004, 'phi below=', 19.964533792945403, 'phi mid=', 19.986845459027418)
('Itn # is:', 3, 'in R3', 'beta above=', 32421.875, 'beta below=', 32226.5625, 'beta mid=', 32324.21875)
('Itn # is:', 3, 'in R3', 'phi above=', 20.009214504197004, 'phi below=', 19.986845459027418, 'phi mid=', 19.998022822238564)
('Itn # is:', 3, 'in R1', 'beta above=', 50000.0, 'beta below=', 32324.21875, 'beta mid=', 41162.109375)
('Itn # is:', 3, 'in R1', 'phi above=', 22.240411728510566, 'phi below=', 19.998022822238564, 'phi mid=', 21.066317451636024)
('Itn # is:', 3, 'in R2', 'beta above=', 41162.109375, 'beta below=', 32324.21875, 'beta mid=', 36743.1640625)
('Itn # is:', 3, 'in R2', 'phi above=', 21.066317451636024, 'phi below=', 19.998022822238564, 'phi mid=', 20.518426873550297)
('Itn # is:', 3, 'in R3', 'beta above=', 41162.109375, 'beta below=', 32324.21875, 'beta mid=', 36743.1640625)
('Itn # is:', 3, 'in R3', 'phi above=', 21.066317451636024, 'phi below=', 19.998022822238564, 'phi mid=', 20.518426873550297)
('Itn # is:', 3, 'in R2', 'beta above=', 36743.1640625, 'beta below=', 32324.21875, 'beta mid=', 34533.69140625)
('Itn # is:', 3, 'in R2', 'phi above=', 20.518426873550297, 'phi below=', 19.998022822238564, 'phi mid=', 20.254690532818699)
('Itn # is:', 3, 'in R3', 'beta above=', 36743.1640625, 'beta below=', 32324.21875, 'beta mid=', 34533.69140625)
('Itn # is:', 3, 'in R3', 'phi above=', 20.518426873550297, 'phi below=', 19.998022822238564, 'phi mid=', 20.254690532818699)
('Itn # is:', 3, 'in R2', 'beta above=', 34533.69140625, 'beta below=', 32324.21875, 'beta mid=', 33428.955078125)
('Itn # is:', 3, 'in R2', 'phi above=', 20.254690532818699, 'phi below=', 19.998022822238564, 'phi mid=', 20.125457930242412)
('Itn # is:', 3, 'in R3', 'beta above=', 34533.69140625, 'beta below=', 32324.21875, 'beta mid=', 33428.955078125)
('Itn # is:', 3, 'in R3', 'phi above=', 20.254690532818699, 'phi below=', 19.998022822238564, 'phi mid=', 20.125457930242412)
('Itn # is:', 3, 'in R2', 'beta above=', 33428.955078125, 'beta below=', 32324.21875, 'beta mid=', 32876.5869140625)
('Itn # is:', 3, 'in R2', 'phi above=', 20.125457930242412, 'phi below=', 19.998022822238564, 'phi mid=', 20.061513584347065)
('Itn # is:', 3, 'in R3', 'beta above=', 33428.955078125, 'beta below=', 32324.21875, 'beta mid=', 32876.5869140625)
('Itn # is:', 3, 'in R3', 'phi above=', 20.125457930242412, 'phi below=', 19.998022822238564, 'phi mid=', 20.061513584347065)
('Itn # is:', 3, 'in R2', 'beta above=', 32876.5869140625, 'beta below=', 32324.21875, 'beta mid=', 32600.40283203125)
('Itn # is:', 3, 'in R2', 'phi above=', 20.061513584347065, 'phi below=', 19.998022822238564, 'phi mid=', 20.029711228024652)
('Itn # is:', 3, 'in R3', 'beta above=', 32876.5869140625, 'beta below=', 32324.21875, 'beta mid=', 32600.40283203125)
('Itn # is:', 3, 'in R3', 'phi above=', 20.061513584347065, 'phi below=', 19.998022822238564, 'phi mid=', 20.029711228024652)
('Itn # is:', 3, 'in R2', 'beta above=', 32600.40283203125, 'beta below=', 32324.21875, 'beta mid=', 32462.310791015625)
('Itn # is:', 3, 'in R2', 'phi above=', 20.029711228024652, 'phi below=', 19.998022822238564, 'phi mid=', 20.013852745745062)
('Itn # is:', 3, 'in R3', 'beta above=', 32600.40283203125, 'beta below=', 32324.21875, 'beta mid=', 32462.310791015625)
('Itn # is:', 3, 'in R3', 'phi above=', 20.029711228024652, 'phi below=', 19.998022822238564, 'phi mid=', 20.013852745745062)
('At end beta=', 32462.310791015625, 'phid=', 20.013852745745062)

This makes the bisection algorithm into a function





In [76]:



    


def bisect_beta(tol, phid, phi_dstar, G,Wd,Wm,dobs,mref,beta):
    i=1
    while np.absolute(phid - phi_dstar)>tol:          
        mrec = mrecovered(G,Wd,Wm,dobs,mref,beta)
        phid = L2(Wd,np.dot(G,mrec)-dobs)
        print("Iteration = ",j, "beta=", beta, "phid at start=",phid)
        if phid>phi_dstar:
            beta = beta_0*0.5**i
            print("beta=", beta)
        if phid<phi_dstar:
            beta_abv = beta_0*0.5**(i-2)
            beta_bel = beta_0*0.5**(i-1)
            beta_mid = 0.5*(beta_abv+beta_bel)
            mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
            phidmid = L2(Wd,np.dot(G,mmid)-dobs)  
            #for k in range(itns):
            while np.absolute(phidmid - phi_dstar)>tol:
                if phidmid<phi_dstar:
                    beta_abv = beta_0*0.5**(i-2)
                    beta_bel = beta_mid
                    beta_mid = 0.5*(beta_abv+beta_bel)
                    beta=beta_mid
                    mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
                    phidmid = L2(Wd,np.dot(G,mmid)-dobs)           
                if phidmid>phi_dstar:
                    beta_abv = beta_mid
                    beta_bel = beta_bel
                    beta_mid = 0.5*(beta_abv+beta_bel)
                    beta=beta_mid
                    mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
                    phidmid = L2(Wd,np.dot(G,mmid)-dobs)                
                    if phidmid<phi_dstar:
                        beta_abv = beta_abv
                        beta_bel = beta_mid
                        beta_mid = 0.5*(beta_abv+beta_bel)
                        beta=beta_mid
                    mmid = mrecovered(G,Wd,Wm,dobs,mref,beta_mid)
                    phidmid = L2(Wd,np.dot(G,mmid)-dobs)              
                phid=phidmid
                mrec=mmid
                return beta, phid, mrec 
        i=i+1

Content source: jokulhaup/directed_studies

Similar notebooks:

notebook.community | gallery | about