Skip to content
Snippets Groups Projects
Commit fa97eb90 authored by Saad Jbabdi's avatar Saad Jbabdi
Browse files

added rbf

parent e74965d4
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
talks/matlab_vs_python/rbf/.ipynb_checkpoints/fit_model-checkpoint.ipynb 0 → 100644
+ 157
− 0
View file @ fa97eb90
%% Cell type:code id: tags:
``` python
# Fit a model to some data
# Model is:
# prediction = M0 * exp(-TE/T2)*(1-exp(-TR/T1))
# where M0,T1,T2 are unknown parameters and TE/TR are experimental parameters
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
TEs = np.array([10,40,60,80]) # TE values in ms
TRs = np.array([.5,1,1.5,2]) # TR in seconds
# All combinations of TEs/TRs
combinations = np.array([(x,y) for x in TEs for y in TRs])
TEs,TRs = combinations[:,0],combinations[:,1]
# function for our model
def forward(p):
M0,T1,T2 = p
return M0*np.exp(-TEs/T2)*(1-np.exp(-TRs/T1))
# simulate data using model
true_p = [100,.8,50]
data = forward(true_p)
data = data + np.random.randn(data.size)
```
%% Output
[<matplotlib.lines.Line2D at 0x121ca4278>]
%% Cell type:code id: tags:
``` python
# Now for the fitting
# we need a cost function:
def cf(p):
pred = forward(p)
return np.mean((pred-data)**2)/2
# always a good idea to calculate gradient
def forward_deriv(p):
M0,T1,T2 = p
E1,E2 = np.exp(-TEs/T2),np.exp(-TRs/T1)
dfdM0 = E2*(1-E1)
dfdT1 = M0*E2*(-E1/T1**2)
dfdT2 = M0*(E2/T2**2)*(1-E1)
return np.array([dfdM0,dfdT1,dfdT2])
def gradient(p):
pred = forward(p)
deriv = forward_deriv(p)
return np.mean( deriv * (pred-data)[None,:],axis=1)
# get ready to minimize
p0 = [200,70,1000] # some random guess
method = 'TNC'
arguments = {'x0':p0,'method':method,'jac':gradient}
result = minimize(cf,**arguments)
```
%% Cell type:code id: tags:
``` python
plt.plot(data)
plt.plot(forward(result.x))
```
%% Output
[<matplotlib.lines.Line2D at 0x121f34eb8>]
%% Cell type:code id: tags:
``` python
```
talks/matlab_vs_python/rbf/.ipynb_checkpoints/play_rbf-checkpoint.ipynb 0 → 100644
+ 97
− 0
View file @ fa97eb90
%% Cell type:code id: tags:
``` python
#!/usr/env/python
# Python version of RBF fitting
import numpy as np
import matplotlib.pyplot as plt
# Generate noisy sine wave
x = np.linspace(0,10,100)
y = np.sin(3*x) + np.random.randn(x.size)*.5
# Define RBF atom
# two options:
# inline function definition (not available in matlab)
# lambda : like @ in matlab
sig = 2
# option 1
# def rbf(x,c):
# return np.exp(-(x-c)**2/sig**2)
# option 2
rbf = lambda x,c : np.exp(-(x-c)**2/sig**2)
# create design matrix
# (use list comprehension to show off)
xi = np.linspace(0,10,15)
desmat = [rbf(x,c) for c in xi]
desmat = np.asarray(desmat).T
# invert model
beta = np.linalg.pinv(desmat)@y.T
# plot fit
plt.figure()
plt.plot(x,y,'.')
plt.plot(x,desmat,'k')
plt.plot(x,desmat@beta)
# make it pretty
plt.grid()
plt.xlabel('x')
plt.ylabel('y')
plt.title('RBF fitting')
plt.savefig('/Users/saad/Desktop/RBF.pdf')
```
%% Output
%% Cell type:code id: tags:
``` python
```
talks/matlab_vs_python/rbf/fit_model.ipynb 0 → 100644
+ 232
− 0
View file @ fa97eb90
%% Cell type:code id: tags:
``` python
# Fit a model to some data
# Model is:
# prediction = M0 * exp(-TE/T2)*(1-exp(-TR/T1))
# where M0,T1,T2 are unknown parameters and TE/TR are experimental parameters
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
TEs = np.array([10,40,60,80]) # TE values in ms
TRs = np.array([.5,1,1.5,2]) # TR in seconds
# All combinations of TEs/TRs
combinations = np.array([(x,y) for x in TEs for y in TRs])
TEs,TRs = combinations[:,0],combinations[:,1]
# function for our model
def forward(p):
M0,T1,T2 = p
return M0*np.exp(-TEs/T2)*(1-np.exp(-TRs/T1))
# simulate data using model
true_p = [100,.8,50]
data = forward(true_p)
data = data + np.random.randn(data.size)
```
%% Output
[<matplotlib.lines.Line2D at 0x121ca4278>]
%% Cell type:code id: tags:
``` python
# Now for the fitting
# we need a cost function:
def cf(p):
pred = forward(p)
return np.mean((pred-data)**2)/2
# always a good idea to calculate gradient
def forward_deriv(p):
M0,T1,T2 = p
E1,E2 = np.exp(-TEs/T2),np.exp(-TRs/T1)
dfdM0 = E2*(1-E1)
dfdT1 = M0*E2*(-E1/T1**2)
dfdT2 = M0*(E2/T2**2)*(1-E1)
return np.array([dfdM0,dfdT1,dfdT2])
def forward_deriv2(p):
M0,T1,T2 = p
E1,E2 = np.exp(-TEs/T2),np.exp(-TRs/T1)
dfdM0dM0 = 0
dfdM0dT1 = -E2
dfdM0dT2 = (1-E1)
dfdT1dM0 = E2*(-E1/T1**2)
dfdT1dT1 = M0*E2*(-E1/T1**4)
dfdT1dT2 = M0*(E2/T2**4)*(1-E1)
dfdT2dM0 = (E2/T2**2)*(1-E1)
dfdT2dT1 = M0*(E2/T2**2)*(1-E1/T1**2)
dfdT2dT2 = M0*(E2/T2**4)*(1-E1)
return np.array([dfdM0dM0,dfdM0dT1,dfdM0dT2],[dfdT1dM0,dfdT1dT1,dfdT1dT2],[dfdT2dM0,dfdT2dT1,dfdT2dT2])
def gradient(p):
pred = forward(p)
deriv = forward_deriv(p)
return np.mean( deriv * (pred-data)[None,:],axis=1)
def hess(p):
pred = forward(p)
deriv = forward_deriv(p)
d2Fdp2 = forward_deriv2(p)
deriv*deriv
# get ready to minimize
p0 = [200,1,70] # some random guess
method = 'Nelder-Mead'
arguments = {'x0':p0,'method':method,'jac':gradient}
result = minimize(cf,**arguments)
```
%% Output
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/scipy/optimize/_minimize.py:506: RuntimeWarning: Method Nelder-Mead does not use gradient information (jac).
RuntimeWarning)
%% Cell type:code id: tags:
``` python
plt.plot(data)
plt.plot(forward(result.x))
```
%% Output
[<matplotlib.lines.Line2D at 0x123074ef0>]
%% Cell type:code id: tags:
``` python
np.mean(forward_deriv(p0),axis=1).shape
```
%% Output
(3,)
%% Cell type:code id: tags:
``` python
true_p
```
%% Output
[100, 0.8, 50]
%% Cell type:code id: tags:
``` python
```
talks/matlab_vs_python/rbf/play_rbf.ipynb 0 → 100644
+ 97
− 0
View file @ fa97eb90
%% Cell type:code id: tags:
``` python
#!/usr/env/python
# Python version of RBF fitting
import numpy as np
import matplotlib.pyplot as plt
# Generate noisy sine wave
x = np.linspace(0,10,100)
y = np.sin(3*x) + np.random.randn(x.size)*.5
# Define RBF atom
# two options:
# inline function definition (not available in matlab)
# lambda : like @ in matlab
sig = 2
# option 1
# def rbf(x,c):
# return np.exp(-(x-c)**2/sig**2)
# option 2
rbf = lambda x,c : np.exp(-(x-c)**2/sig**2)
# create design matrix
# (use list comprehension to show off)
xi = np.linspace(0,10,15)
desmat = [rbf(x,c) for c in xi]
desmat = np.asarray(desmat).T
# invert model
beta = np.linalg.pinv(desmat)@y.T
# plot fit
plt.figure()
plt.plot(x,y,'.')
plt.plot(x,desmat,'k')
plt.plot(x,desmat@beta)
# make it pretty
plt.grid()
plt.xlabel('x')
plt.ylabel('y')
plt.title('RBF fitting')
plt.savefig('/Users/saad/Desktop/RBF.pdf')
```
%% Output
%% Cell type:code id: tags:
``` python
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment