Skip to content

Commit

Permalink
refactor of limbdark
Browse files Browse the repository at this point in the history
  • Loading branch information
@dfm
dfm committed Dec 5, 2018
1 parent d5651d6 commit df1bdf4
Show file tree
Hide file tree
Showing 6 changed files with 77 additions and 161 deletions.
2 changes: 1 addition & 1 deletion exoplanet/light_curves.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -166,4 +166,4 @@ def _compute_light_curve(self, b, r, los=None):
"""
if los is None:
los = -tt.ones_like(b)
return limbdark(self.c_norm, b, r, los)
return limbdark(self.c_norm, b, r, los)[0]
52 changes: 41 additions & 11 deletions exoplanet/theano_ops/starry/limbdark.cc
View file
Original file line number Diff line number Diff line change
@@ -1,11 +1,15 @@
#section support_code_apply

int APPLY_SPECIFIC(limbdark)(
PyArrayObject* input0, // Array of "cl"
PyArrayObject* input1, // Array of impact parameters "b"
PyArrayObject* input2, // Array of radius ratios "r"
PyArrayObject* input3, // Array of line-of-sight position "los"
PyArrayObject** output0)
PyArrayObject* input0, // Array of "cl"
PyArrayObject* input1, // Array of impact parameters "b"
PyArrayObject* input2, // Array of radius ratios "r"
PyArrayObject* input3, // Array of line-of-sight position "los"
PyArrayObject** output0, // Flux
PyArrayObject** output1, // dfdcl
PyArrayObject** output2, // dfdb
PyArrayObject** output3 // dfdr
)
{
npy_intp Nc, Nb, Nr, Nlos;
int success = get_size(input0, &Nc);
Expand All @@ -18,28 +22,54 @@ int APPLY_SPECIFIC(limbdark)(
return 1;
}

success += allocate_output(PyArray_NDIM(input1), PyArray_DIMS(input1), TYPENUM_OUTPUT_0, output0);
npy_intp ndim = PyArray_NDIM(input1);
npy_intp* dims = PyArray_DIMS(input1);
std::vector<npy_intp> shape(ndim + 1);
shape[0] = Nc;
for (npy_intp i = 0; i < ndim; ++i) shape[i+1] = dims[i];

success += allocate_output(ndim, dims, TYPENUM_OUTPUT_0, output0);
success += allocate_output(ndim+1, &(shape[0]), TYPENUM_OUTPUT_1, output1);
success += allocate_output(ndim, dims, TYPENUM_OUTPUT_2, output2);
success += allocate_output(ndim, dims, TYPENUM_OUTPUT_3, output3);
if (success) {
return 1;
}

DTYPE_INPUT_0* c = (DTYPE_INPUT_0*) PyArray_DATA(input0);
DTYPE_INPUT_1* b = (DTYPE_INPUT_1*) PyArray_DATA(input1);
DTYPE_INPUT_2* r = (DTYPE_INPUT_2*) PyArray_DATA(input2);
DTYPE_INPUT_0* c = (DTYPE_INPUT_0*) PyArray_DATA(input0);
DTYPE_INPUT_1* b = (DTYPE_INPUT_1*) PyArray_DATA(input1);
DTYPE_INPUT_2* r = (DTYPE_INPUT_2*) PyArray_DATA(input2);
DTYPE_INPUT_3* los = (DTYPE_INPUT_3*) PyArray_DATA(input3);
DTYPE_OUTPUT_0* f = (DTYPE_OUTPUT_0*)PyArray_DATA(*output0);

DTYPE_OUTPUT_0* f = (DTYPE_OUTPUT_0*)PyArray_DATA(*output0);
DTYPE_OUTPUT_1* dfdcl = (DTYPE_OUTPUT_1*)PyArray_DATA(*output1);
DTYPE_OUTPUT_2* dfdb = (DTYPE_OUTPUT_2*)PyArray_DATA(*output2);
DTYPE_OUTPUT_3* dfdr = (DTYPE_OUTPUT_3*)PyArray_DATA(*output3);

Eigen::Map<Eigen::Matrix<DTYPE_OUTPUT_1, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> dfdcl_mat(dfdcl, Nc, Nb);
dfdcl_mat.setZero();

Eigen::Map<Eigen::Matrix<DTYPE_INPUT_0, Eigen::Dynamic, 1>> cvec(c, Nc);
starry::limbdark::GreensLimbDark<DTYPE_OUTPUT_0> L(Nc-1);

for (npy_intp i = 0; i < Nb; ++i) {
f[i] = 0;
dfdb[i] = 0;
dfdr[i] = 0;

if (los[i] < 0) {
auto b_ = std::abs(b[i]);
auto r_ = std::abs(r[i]);
if (b_ < 1 + r_) {
L.compute(b_, r_, false);
L.compute(b_, r_, true);

// The value of the light curve
f[i] = L.S.dot(cvec) - 1;

// The gradients
dfdcl_mat.col(i) = L.S;
dfdb[i] = sgn(b[i]) * L.dSdb.dot(cvec);
dfdr[i] = sgn(r[i]) * L.dSdr.dot(cvec);
}
}
}
Expand Down
45 changes: 31 additions & 14 deletions exoplanet/theano_ops/starry/limbdark.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,11 +4,11 @@

__all__ = ["LimbDarkOp"]

import theano
from theano import gof
import theano.tensor as tt

from .base_op import StarryBaseOp
from .limbdark_rev import LimbDarkRevOp


class LimbDarkOp(StarryBaseOp):
Expand All @@ -17,27 +17,44 @@ class LimbDarkOp(StarryBaseOp):
func_file = "./limbdark.cc"
func_name = "APPLY_SPECIFIC(limbdark)"

def __init__(self):
self.grad_op = LimbDarkRevOp()
super(LimbDarkOp, self).__init__()

def make_node(self, c, b, r, los):
in_args = [
tt.as_tensor_variable(c),
tt.as_tensor_variable(b),
tt.as_tensor_variable(r),
tt.as_tensor_variable(los),
in_args = []
dtype = theano.config.floatX
for a in [c, b, r, los]:
try:
a = tt.as_tensor_variable(a)
except tt.AsTensorError:
pass
else:
dtype = theano.scalar.upcast(dtype, a.dtype)
in_args.append(a)

out_args = [
in_args[1].type(),
tt.TensorType(dtype=dtype,
broadcastable=[False] * (in_args[1].ndim + 1))(),
in_args[1].type(),
in_args[2].type(),
]
out_args = [in_args[1].type()]
return gof.Apply(self, in_args, out_args)

def infer_shape(self, node, shapes):
return shapes[1],
return (
shapes[1], list(shapes[0]) + list(shapes[1]),
shapes[1], shapes[2])

def grad(self, inputs, gradients):
c, b, r, los = inputs
bf, = gradients
bc, bb, br = self.grad_op(c, b, r, los, bf)
f, dfdcl, dfdb, dfdr = self(*inputs)
bf = gradients[0]
for i, g in enumerate(gradients[1:]):
if not isinstance(g.type, theano.gradient.DisconnectedType):
raise ValueError("can't propagate gradients wrt parameter {0}"
.format(i+1))
bc = tt.sum(tt.reshape(bf, (1, bf.size)) *
tt.reshape(dfdcl, (c.size, bf.size)), axis=-1)
bb = bf * dfdb
br = bf * dfdr
return bc, bb, br, tt.zeros_like(los)

def R_op(self, inputs, eval_points):
Expand Down
65 changes: 0 additions & 65 deletions exoplanet/theano_ops/starry/limbdark_rev.cc
View file

This file was deleted.

31 changes: 0 additions & 31 deletions exoplanet/theano_ops/starry/limbdark_rev.py
View file

This file was deleted.

43 changes: 4 additions & 39 deletions exoplanet/theano_ops/starry/limbdark_test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,6 @@
from theano.tests import unittest_tools as utt

from .limbdark import LimbDarkOp
from .limbdark_rev import LimbDarkRevOp


class TestLimbDark(utt.InferShapeTester):
Expand All @@ -24,7 +23,7 @@ def get_args(self):
b = tt.vector()
r = tt.vector()
los = tt.vector()
f = theano.function([c, b, r, los], self.op(c, b, r, los))
f = theano.function([c, b, r, los], self.op(c, b, r, los)[0])

c_val = np.array([-0.85, 2.5, -0.425, 0.1])
b_val = np.linspace(-1.5, 1.5, 100)
Expand All @@ -48,45 +47,11 @@ def test_los(self):
def test_infer_shape(self):
f, args, arg_vals = self.get_args()
self._compile_and_check(args,
[self.op(*args)],
self.op(*args),
arg_vals,
self.op_class)

def test_grad(self):
_, _, in_args = self.get_args()
utt.verify_grad(self.op, in_args)


class TestLimbDarkRev(utt.InferShapeTester):

def setUp(self):
super(TestLimbDarkRev, self).setUp()
self.op_class = LimbDarkRevOp
self.op = LimbDarkRevOp()

def get_args(self):
c = tt.vector()
b = tt.vector()
r = tt.vector()
los = tt.vector()
bf = tt.vector()
f = theano.function([c, b, r, los, bf], self.op(c, b, r, los, bf))

c_val = np.array([-0.85, 2.5, -0.425, 0.1])
b_val = np.linspace(-1.5, 1.5, 100)
r_val = 0.1 + np.zeros_like(b_val)
los_val = -np.ones_like(b_val)
bf_val = np.ones_like(b_val)

return f, [c, b, r, los, bf], [c_val, b_val, r_val, los_val, bf_val]

def test_basic(self):
f, _, in_args = self.get_args()
f(*in_args)

def test_infer_shape(self):
f, args, arg_vals = self.get_args()
self._compile_and_check(args,
self.op(*args),
arg_vals,
self.op_class)
func = lambda *args: self.op(*args)[0] # NOQA
utt.verify_grad(func, in_args)

0 comments on commit df1bdf4

Please sign in to comment.