In [3]:
using HDF5, JLD
In [4]:
# add additional processes to serve as workers
n_procs = 6
addprocs(n_procs);
In [5]:
# add present working directory to path from which to load modules
push!(LOAD_PATH, ".")
# load module DynamicProgramming on each of the processes
using DynamicProgramming
using Interpolations
In [6]:
# nominal parameter value
@everywhere T1 = 0.680 # seconds in white matter
@everywhere T2 = 0.090 # seconds in white matter
@everywhere dT = 0.2 # seconds between acquisitions
@everywhere sigma = 0.1 # noise strength
@everywhere theta0 = [exp(-dT/T1); exp(-dT/T2)]
# time horizon
@everywhere N = 30
# state dimension
@everywhere n = 6
# initial state
@everywhere x0 = [0.; 0.; 1.; 0; 0; 0]
# define rotation matrices
@everywhere function Rx(alpha::Float64)
return [1 0 0 ;
0 cos(alpha) -sin(alpha);
0 sin(alpha) cos(alpha)]
end
@everywhere function Ry(alpha::Float64)
return [cos(alpha) 0 -sin(alpha);
0 1 0 ;
sin(alpha) 0 cos(alpha)]
end
@everywhere function Rz(alpha::Float64)
return [cos(alpha) -sin(alpha) 0;
sin(alpha) cos(alpha) 0;
0 0 1]
end
# dynamics
@everywhere function f(t::Int64, x::Array{Float64, 1}, u::Array{Float64, 1}, theta::Array{Float64, 1})
x_state = x[1:3]
x_sensitivity = x[4:6]
U = Rz(u[1])*Ry(u[2])*Rx(u[3])
D = [theta[2], theta[2], theta[1]]
x_state_plus = U*(D.*x_state) + [0; 0; 1-theta[1]]
x_sensitivity_plus = U*(D.*x_sensitivity) + U*([0; 0; 1].*x_state) + [0; 0; -1]
return [x_state_plus; x_sensitivity_plus]
end
# test dynamics
f(1, x0, [0.; 0.; 0.], theta0)
# reward function
@everywhere function g(k::Int64, x::Array{Float64, 1}, u::Array{Float64, 1}, theta::Array{Float64, 1})
return x[4]^2 + x[5]^2
end
# terminal reward
@everywhere function g_terminal(x::Array{Float64, 1})
return x[4]^2 + x[5]^2
end
# test reward function
g(0, x0, [0.; 0.; 0.], theta0)
# define state grid
@everywhere nx1 = 6 # number of grid points in each dimension
@everywhere nx2 = 6
@everywhere nx3 = 6
@everywhere nx4 = 7
@everywhere nx5 = 7
@everywhere nx6 = 7
@everywhere xgrid = (
linspace(-1., 1., nx1),
linspace(-1., 1., nx2),
linspace(-1., 1., nx3),
linspace(-50., 50., nx4),
linspace(-50., 50., nx5),
linspace(-100., 100., nx6)
)
# define input grid
# make nu's even so that there is a grid knot at zero
@everywhere nu1 = 1
@everywhere nu2 = 16
@everywhere nu3 = 8
@everywhere ugrid = (
linspace(0, 0, nu1),
linspace(-pi, pi-2*pi/nu2, nu2),
linspace(0, pi-pi/nu3, nu3)
)
In [7]:
# define invariants
@everywhere function rho(x::Array{Float64, 1})
r = sqrt(x[1]*x[1] + x[2]*x[2])
return [r;
x[3];
1/r*(x[2]*x[4] - x[1]*x[5]);
1/r*(x[1]*x[4] + x[2]*x[5]);
x[6]]
end
@everywhere function rho_bar_inverse(x_bar::Array{Float64, 1})
return [0; x_bar]
end
In [10]:
# @everywhere nx1_reduced = 6
# @everywhere nx2_reduced = 15
# @everywhere nx3_reduced = 20
# @everywhere nx4_reduced = 20
# @everywhere nx5_reduced = 20
@everywhere nx1_reduced = 6
@everywhere nx2_reduced = 10
@everywhere nx3_reduced = 15
@everywhere nx4_reduced = 15
@everywhere nx5_reduced = 15
@everywhere xgrid_reduced = (
linspace(0.01, 0.6, nx1_reduced),
linspace(1-2*theta0[1], 1., nx2_reduced),
linspace(-5., 5., nx3_reduced),
linspace(-5., 5., nx4_reduced),
linspace(-5., 5., nx5_reduced)
)
In [6]:
@time J = dp_reduced(f, g, g_terminal, rho, rho_bar_inverse, ugrid, xgrid_reduced, theta0, N);
In [8]:
# save value function computed at grid points
# save("J_MR_fingerprinting.jld", "data", J)
In [8]:
# reload value function
J = load("J_MR_fingerprinting.jld")["data"];
In [11]:
# define a finer input grid for policy rollout
@everywhere nu1 = 32
@everywhere nu2 = 32
@everywhere nu3 = 16
@everywhere ugrid = (
linspace(0, 0, 1),
linspace(-pi, pi-2*pi/nu2, nu2),
linspace(0, pi-pi/nu3, nu3)
)
# compute optimal trajectory from initial condition x0
x0 = [0.; 0.; 1.; 0.; 0.; 0.]
x_opt, u_opt = dp_rollout_reduced(J, x0, f, g, rho, ugrid, xgrid_reduced, theta0, N);
In [12]:
using Colors
using Gadfly
berkeley_blue = RGB((1/256*[ 45., 99., 127.])...)
berkeley_gold = RGB((1/256*[224., 158., 25.])...);
gray = RGB((1/256*[200., 200., 200.])...);
In [14]:
# plot optimal input
p = plot(layer(x=0:N-1, y=180/pi*u_opt[1, :], Geom.point, Geom.line, Theme(default_color=berkeley_blue)),
layer(x=0:N-1, y=180/pi*u_opt[2, :], Geom.point, Geom.line, Theme(default_color=berkeley_gold)),
layer(x=0:N-1, y=180/pi*u_opt[3, :], Geom.point, Geom.line, Theme(default_color=gray)),
Guide.XLabel("time"),
Guide.YLabel("input angle (degrees)"),
Guide.manual_color_key("Angle", ["α", "β", "δ"], [berkeley_blue, berkeley_gold, gray])
)
draw(PNG("MRI_optimal_input.png", 15cm, 10cm), p)
p
Out[14]:
In [13]:
# plot optimal trajectory in transverse plane
p = plot(x=x_opt[1, :], y=x_opt[2, :],
Guide.XLabel("x_1 magnetization"),
Guide.YLabel("x_2 magnetization"),
Geom.path, Geom.point, Coord.cartesian(fixed=true),
Theme(default_color=berkeley_gold)
)
draw(SVG("MRI_full_trajectory_transverse_plane.svg", 10cm, 20cm), p)
p
Out[13]:
In [14]:
# plot optimal trajectory in transverse vs. longitudinal
p = plot(x=sqrt(x_opt[1, :].^2 + x_opt[2, :].^2), y=x_opt[3, :],
Guide.XLabel("transverse magnetization"),
Guide.YLabel("longitudinal magnetization"),
Geom.path, Geom.point, Coord.cartesian(fixed=true),
Theme(default_color=berkeley_gold)
)
draw(SVG("MRI_full_trajectory_transverse_vs_longitudinal.svg", 10cm, 20cm), p)
p
Out[14]:
In [15]:
p = plot(layer(x=0:N, y=x_opt[1, :], Geom.point, Geom.line, Theme(default_color=berkeley_blue)),
layer(x=0:N, y=x_opt[2, :], Geom.point, Geom.line, Theme(default_color=berkeley_gold)),
layer(x=0:N, y=x_opt[3, :], Geom.point, Geom.line, Theme(default_color=gray)),
Guide.XLabel("time"),
Guide.YLabel("state"),
# Guide.manual_color_key("Angle", ["α", "β", "δ"], [berkeley_blue, berkeley_gold, gray])
)
Out[15]:
In [15]:
p = plot(layer(x=0:N, y=sqrt(x_opt[1, :].^2 + x_opt[2, :].^2), Geom.point, Geom.line, Theme(default_color=berkeley_blue)),
layer(x=0:N, y=x_opt[3, :], Geom.point, Geom.line, Theme(default_color=berkeley_gold)),
# layer(x=0:N, y=x_opt[6, :], Geom.point, Geom.line, Theme(default_color=gray)),
Guide.XLabel("time"),
Guide.YLabel("magnetization"),
Guide.manual_color_key("", ["transverse |(x_1, x_2)|", "longitudinal (x_3)"], [berkeley_blue, berkeley_gold])
)
draw(PNG("MRI_optimal_magnetization.png", 15cm, 10cm), p)
p
Out[15]:
In [17]:
p = plot(layer(x=0:N, y=x_opt[4, :], Geom.point, Geom.line, Theme(default_color=berkeley_blue)),
layer(x=0:N, y=x_opt[5, :], Geom.point, Geom.line, Theme(default_color=berkeley_gold)),
layer(x=0:N, y=x_opt[6, :], Geom.point, Geom.line, Theme(default_color=gray)),
Guide.XLabel("time"),
Guide.YLabel("sensitivity"),
# Guide.manual_color_key("Angle", ["α", "β", "δ"], [berkeley_blue, berkeley_gold, gray])
)
Out[17]:
In [16]:
p = plot(layer(x=0:N, y=sqrt(x_opt[4, :].^2 + x_opt[5, :].^2), Geom.point, Geom.line, Theme(default_color=berkeley_blue)),
layer(x=0:N, y=x_opt[6, :], Geom.point, Geom.line, Theme(default_color=berkeley_gold)),
# layer(x=0:N, y=x_opt[6, :], Geom.point, Geom.line, Theme(default_color=gray)),
Guide.XLabel("time"),
Guide.YLabel("sensitivity"),
Guide.manual_color_key("", ["transverse |(x_4, x_5)|", "longitudinal (x_6)"], [berkeley_blue, berkeley_gold])
)
draw(PNG("MRI_optimal_sensitivity.png", 15cm, 10cm), p)
p
Out[16]:
In [19]:
Pkg.status()
In [ ]: