Git Product home page Git Product logo

starry_process's Introduction


Documentation Status DOI badge

Interpretable Gaussian processes for stellar light curves using starry.

A Gaussian Process for Stellar Variability

The starry_process code implements an interpretable Gaussian process (GP) for modeling stellar light curves. Whether your goal is to marginalize over the stellar variability signal (if you think of it as noise) or to understand the surface features that generated it (if you think of it as data), this code is for you. The GP implemented here works just like any other GP you might already use in your analysis, except that its hyperparameters are physically interpretable. These are (among others) the radius of the spots, the mean and variance of the latitude distribution, the spot contrast, and the number of spots. Users can also specify things like the rotational period of the star, the limb darkening parameters, and the inclination (or marginalize over the inclination if it is not known).

The code is written in Python and relies on the Theano package, so a little familiarity with that is recommended. Check out the crash course here. If you would like to report an issue or contribute to the project, please check out CONTRIBUTING.md.

Installation

The quickest way is via pip:

pip install starry-process

Note that the starry_process package requires Python 3.6 or later.

Quickstart

Import the main interface:

from starry_process import StarryProcess

Draw samples from a Gaussian process with small mid-latitude spots:

import numpy as np
import matplotlib.pyplot as plt

# Instantiate the GP
sp = StarryProcess(
  r=10,               # spot radius in degrees
  mu=30,              # central spot latitude in degrees
  sigma=5,            # latitude std. dev. in degrees
  c=0.1,              # fractional spot contrast
  n=10                # number of spots
)

# Draw & visualize a spherical harmonic sample
y = sp.sample_ylm().eval()
sp.visualize(y)

# Compute & plot the flux at some inclination
t = np.linspace(0, 4, 1000)
flux = sp.flux(y, t, i=60).eval()[0]
plt.plot(t, flux)

Same as above, but for high-latitude spots:

sp = StarryProcess(r=10, mu=0, sigma=10, c=0.1, n=10)

Large equatorial spots:

sp = StarryProcess(r=30, mu=0, sigma=10, c=0.1, n=10)

Small, approximately isotropic spots:

sp = StarryProcess(r=10, mu=0, sigma=40, c=0.1, n=10)

For more information check out the full Quickstart tutorial and the complete documentation.

References & Attribution

The code is described in this JOSS paper. It is the backbone of the Mapping Stellar Surfaces paper series, including:

If you make use of this code in your research, please cite

@article{Luger2021a,
  author        = {{Luger}, Rodrigo and {Foreman-Mackey}, Daniel and {Hedges}, Christina and {Hogg}, David W.},
  title         = {{Mapping stellar surfaces I: Degeneracies in the rotational light curve problem}},
  journal       = {arXiv e-prints},
  keywords      = {Astrophysics - Solar and Stellar Astrophysics, Astrophysics - Instrumentation and Methods for Astrophysics},
  year          = 2021,
  month         = jan,
  eid           = {arXiv:2102.00007},
  pages         = {arXiv:2102.00007},
  archiveprefix = {arXiv},
  eprint        = {2102.00007},
  primaryclass  = {astro-ph.SR},
  adsurl        = {https://ui.adsabs.harvard.edu/abs/2021arXiv210200007L},
  adsnote       = {Provided by the SAO/NASA Astrophysics Data System}
}
@article{Luger2021b,
  author        = {{Luger}, Rodrigo and {Foreman-Mackey}, Daniel and {Hedges}, Christina},
  title         = {{Mapping stellar surfaces II: An interpretable Gaussian process model for light curves}},
  journal       = {arXiv e-prints},
  keywords      = {Astrophysics - Solar and Stellar Astrophysics, Astrophysics - Earth and Planetary Astrophysics, Astrophysics - Instrumentation and Methods for Astrophysics},
  year          = 2021,
  month         = feb,
  eid           = {arXiv:2102.01697},
  pages         = {arXiv:2102.01697},
  archiveprefix = {arXiv},
  eprint        = {2102.01697},
  primaryclass  = {astro-ph.SR},
  adsurl        = {https://ui.adsabs.harvard.edu/abs/2021arXiv210201697L},
  adsnote       = {Provided by the SAO/NASA Astrophysics Data System}
}
@article{Luger2021c,
  author        = {{Luger}, Rodrigo and {Foreman-Mackey}, Daniel and {Hedges}, Christina},
  title         = {{starry\_process: Interpretable Gaussian processes for stellar light curves}},
  journal       = {arXiv e-prints},
  keywords      = {Astrophysics - Solar and Stellar Astrophysics, Astrophysics - Earth and Planetary Astrophysics, Astrophysics - Instrumentation and Methods for Astrophysics},
  year          = 2021,
  month         = feb,
  eid           = {arXiv:2102.01774},
  pages         = {arXiv:2102.01774},
  archiveprefix = {arXiv},
  eprint        = {2102.01774},
  primaryclass  = {astro-ph.SR},
  adsurl        = {https://ui.adsabs.harvard.edu/abs/2021arXiv210201774L},
  adsnote       = {Provided by the SAO/NASA Astrophysics Data System}
}

starry_process's People

Contributors

arfon avatar dfm avatar rodluger avatar toihr avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar

starry_process's Issues

Fix theano dependencies

The current code requires theano<=1.0.5 and fails if theano-pymc3 is installed. We need to make it forward-compatible.

To make things a bit easier, I removed the explicit dependency of the code on starry, which installs pymc3>=3.10.0 (and thus theano-pymc3) if we're not careful. Note that starry is required for the tests but it won't be installed unless the [tests] extra flag is passed to pip.

We'll have to fix all these issues for starry as well, but let's take it one step at a time.

Following quickstart something goes wrong with Theano

Hi, after I installed dependencies and tried to run the quickstart, appears this error.

In [2]: import numpy as np

In [3]: import matplotlib.pyplot as plt

In [4]: from tqdm.auto import tqdm

In [5]: sp = StarryProcess()

In [6]: sp.r
Out[6]: Elemwise{mul,no_inplace}.0

In [7]: sp.r.eval()
Out[7]: array(20.)

In [8]: y = sp.sample_ylm().eval()

You can find the C code in this temporary file: /tmp/theano_compilation_error_iie8tb6w
ERROR (theano.graph.opt): Optimization failure due to: constant_folding
ERROR (theano.graph.opt): node: LatitudeIntegralOp{ydeg=15, udeg=2, compile_args=()}(TensorConstant{54.598150033144236}, TensorConstant{8.971048493895104})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
File "/home/jserna/Astro/TOOLS/anaconda3/lib/python3.7/site-packages/theano/graph/opt.py", line 2017, in process_node
replacements = lopt.transform(fgraph, node)
File "/home/jserna/Astro/TOOLS/anaconda3/lib/python3.7/site-packages/theano/graph/opt.py", line 1209, in transform
return self.fn(*args, **kwargs)
File "/home/jserna/Astro/TOOLS/anaconda3/lib/python3.7/site-packages/theano/tensor/opt.py", line 7007, in constant_folding
node, storage_map, compute_map, no_recycling=[], impl=impl
File "/home/jserna/Astro/TOOLS/anaconda3/lib/python3.7/site-packages/theano/graph/op.py", line 634, in make_thunk
return self.make_c_thunk(node, storage_map, compute_map, no_recycling)
File "/home/jserna/Astro/TOOLS/anaconda3/lib/python3.7/site-packages/theano/graph/op.py", line 601, in make_c_thunk
input_storage=node_input_storage, output_storage=node_output_storage
File "/home/jserna/Astro/TOOLS/anaconda3/lib/python3.7/site-packages/theano/link/c/basic.py", line 1204, in make_thunk
input_storage, output_storage, storage_map
File "/home/jserna/Astro/TOOLS/anaconda3/lib/python3.7/site-packages/theano/link/c/basic.py", line 1142, in compile
storage_map,
File "/home/jserna/Astro/TOOLS/anaconda3/lib/python3.7/site-packages/theano/link/c/basic.py", line 1634, in cthunk_factory
module = get_module_cache().module_from_key(key=key, lnk=self)
File "/home/jserna/Astro/TOOLS/anaconda3/lib/python3.7/site-packages/theano/link/c/cmodule.py", line 1191, in module_from_key
module = lnk.compile_cmodule(location)
File "/home/jserna/Astro/TOOLS/anaconda3/lib/python3.7/site-packages/theano/link/c/basic.py", line 1550, in compile_cmodule
preargs=preargs,
File "/home/jserna/Astro/TOOLS/anaconda3/lib/python3.7/site-packages/theano/link/c/cmodule.py", line 2547, in compile_str
f"Compilation failed (return status={status}): {compile_stderr}"
Exception: ('Compilation failed (return status=1): g++: error: unrecognized command line option ‘-std=c++14’. ', 'FunctionGraph(*1 -> LatitudeIntegralOp{ydeg=15, udeg=2, compile_args=()}(TensorConstant{54.598150033144236}, TensorConstant{8.971048493895104}), *1::1, *1::2, *1::3, *1::4, *1::5)')

Any idea of how to solve it?
Appreciate if you find a solution.

My comments on the manuscript

I will continue to edit this as I go through the paper.

Abstract:

  • Most people wouldn't agree that GPs are used for their "computational efficiency" :D
  • I don't agree with the statement that physically motivated GPs haven't been used to good effect. There are plenty of solid examples of measuring rotation periods and asteroseismic parameters using GPs. I think the better statement is that the specific things that you want to measure are hard, regardless of method. I now explicitly and exclusively refer to "star spot" variability. Let me know if that's better.
  • "ridden by degeneracies"? Is "ridden" the right word?
  • "of the Gaussian process that best describes" -> "of an effective GP model to describe"
  • "when employed in an inference setting"

Introduction

  • I think this section should be re-structured to lead with the astronomical application rather than an introduction to GPs. In fact, I think that some of the pedagogy could be removed or moved to another section/appendix. Instead, start with why it's interesting to measure these spot properties and why it's hard.
  • I think the word "extremely" is being overused: it is neither "extremely" easy or efficient to work with GPs. Perhaps
    "well-defined"? Lol. Ctrl+F yields 27 instances of "extremely"! I should look at a word cloud for this paper...
  • Eq 3-6, I'm not sure that it's necessary to include all these example kernels here. I want to keep them because I reference them later (in the temporal GP discussion of the Extensions section) and I need to motivate how many of their hyperparameters are not interpretable.
  • There has been some work on empirically interpreting the parameters of these standard GP kernels and their relationship to physical quantities by fitting simulated data: e.g. https://arxiv.org/abs/2012.01862 and references therein. These efforts should be cited!

[JOSS review] List dependencies on README.md and/or documentation

Hey @rodluger!

A very minor item on my review (openjournals/joss-reviews#3071) is to please list the dependencies of starry_process clearly, either on the README.md and/or the documentation. I saw these are actually in the setup.py (and actually handled automatically as written in the installation instructions!), but listing them might help users figure out what they need (e.g., only Python 3 apparently as a scipy version that doesn't support Python 2.7 at least is listed!).

Thanks in advance!

N.

[JOSS] README example doesn't execute

README example doesn't execute with the pip installable version (0.9.4), there's no draw_ylm function (sample_ylm?), and no visualize function that I can execute on the object either.

You also need to import numpy at the top. <3

name 'p0' is not defined

when i run the following code(after optimizing in Ensemble analysis), it shows the following error


# Number of parameters
ndim = p0.shape[1]

# Instantiate the sampler
sampler = emcee.EnsembleSampler(nwalkers, ndim, mci.logp)

# Run the chains
np.random.seed(0)
nsteps = 500
_ = sampler.run_mcmc(p0, nsteps, progress=True)

Numerical issues at high SNR

When the error on the flux is < 1 ppt, inference breaks down. This can be seen in this test, which often fails (by a lot!) if ferr < 1e-3. Need to investigate this in more detail.

Comments on manuscript

Awesome work, here are some comments

  • "Dark spots arise due to the suppression of convection in regions of intense magnetic field, resulting in a locally cooler (and hence darker) photosphere." Can you add a citation or link to nice review? Ignore if it's Berdyugina 2005 It is!

  • "In discrete spot models like the ones discussed above, the degeneracy-breaking prior is (typically) the assumption that the spots must be circular, have uniform contrast, and sit atop an otherwise uniform photosphere." Do the gridded models also assume an active latitude? If yes, can you add those words here? Not usually, I think, but it would be easy to incorporate such a prior.

  • Section 4.1, first 3 paragraphs: You did already introduce the null space in enough detail in the introduction, I think you can cut these paragraphs down fairly significantly, but I don't feel very strongly. I think this is now fine given the splitting into two papers.

  • "This has the effect of reversing the null space: under only linear limb darkening, it is the even modes that would like in the null space." Do you mean be in the null space?

  • "In reality, the true limb darkening profile of a stellar surface is more complicated than a two parameter quadratic model can capture" Do you have a nice citation for a work we could read about limb darkening laws, and why they are more complicated?

  • "In practice, of course, we will never know the limb darkening coefficients exactly, so we must either model them or marginalize over them (see §5.3). This will dilute our ability to infer surface properties exactly." Planets and occultations would greatly help you constrain limb darkening, might be nice to nod to that here.

  • "As a rule of thumb, if the amplitude of variability exceeds 10%, we recommend not normalizing the light curve in this way," I think in the normalization problem section you mentioned that 2-4% was the cut off?

  • "This implies that there is no preferred time (or phase) and that the spatial covariance is the same at all points in time. " + "(i.e. spots to do not preferentially form at the same longitude on the surface)"

  • At the end you list future work, but don't say you're going to implement the algorithm.......maybe add a bullet point saying you'll apply it to satisfy an observation focused reviewer?

  • Time evolution stuff -could- be a separate paper

  • Section 5.4, by using in conjunction with starry, will this greatly improve the constraint on the spot properties, and allow to constrain full map (inc null space)?

[JOSS review] List Community Guidelines

Hey @rodluger!

One check we need to do as reviewers of a JOSS submission (openjournals/joss-reviews#3071) are Community Guidelines; basically, a text snippet outlining how to: 1) Contribute to the software 2) Report issues or problems with the software 3) Seek support. While to me the most likely path you are suggesting is to perform these through issue handling/pull requests (as the issue-reporting is implicitly defined in the documentation), I think these should be described explicitly (either in the README and/or documentation) so the workflow is defined without ambiguities for users.

Thanks again for this fantastic piece of software!
Néstor

ValueError from Quickstart

Hi @rodluger,

I tried running the first example in the Quickstart and got a ValueError when plotting the flux (which boils down to the shape of the output)

ValueError: x and y must have same first dimension, but have shapes (1000,) and (1, 1000)

(This is part of joss-reviews/issues/3071)

error: 'M_PI' was not declared in this scope

I am trying to use starry_process on Win10. And i keep getting this kind of error: Even if i just want to do this simple line:
y = sp.sample_ylm().eval()

The Error Logs are something about M_PI was not declared in this scope:

Problem occurred during compilation with the command line below:
"D:\Software\Programming\mingw-64\mingw64\bin\g++.exe" -shared -g -O3 -fno-math-errno -Wno-unused-label -Wno-unused-variable -Wno-write-strings -std=c++14 -O2 -DNDEBUG -DSP__LMAX=15 -DSP__UMAX=2 -Wno-c++11-narrowing -march=haswell -mmmx -mno-3dnow -msse -msse2 -msse3 -mssse3 -mno-sse4a -mcx16 -msahf -mmovbe -maes -mno-sha -mpclmul -mpopcnt -mabm -mno-lwp -mfma -mno-fma4 -mno-xop -mbmi -mno-sgx -mbmi2 -mno-pconfig -mno-wbnoinvd -mno-tbm -mavx -mavx2 -msse4.2 -msse4.1 -mlzcnt -mno-rtm -mno-hle -mrdrnd -mf16c -mfsgsbase -mno-rdseed -mno-prfchw -mno-adx -mfxsr -mxsave -mxsaveopt -mno-avx512f -mno-avx512er -mno-avx512cd -mno-avx512pf -mno-prefetchwt1 -mno-clflushopt -mno-xsavec -mno-xsaves -mno-avx512dq -mno-avx512bw -mno-avx512vl -mno-avx512ifma -mno-avx512vbmi -mno-avx5124fmaps -mno-avx5124vnniw -mno-clwb -mno-mwaitx -mno-clzero -mno-pku -mno-rdpid -mno-gfni -mno-shstk -mno-avx512vbmi2 -mno-avx512vnni -mno-vaes -mno-vpclmulqdq -mno-avx512bitalg -mno-movdiri -mno-movdir64b --param l1-cache-size=32 --param l1-cache-line-size=64 --param l2-cache-size=8192 -mtune=haswell -DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION -m64 -DMS_WIN64 -I"D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include" -I"D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\vendor\eigen_3.3.5" -I"D:\Software\Programming\Anaconda\lib\site-packages\numpy\core\include" -I"D:\Software\Programming\Anaconda\include" -I"D:\Software\Programming\Anaconda\lib\site-packages\theano\link\c\c_code" -L"D:\Software\Programming\Anaconda\libs" -L"D:\Software\Programming\Anaconda" -o "C:\Users\diese\AppData\Local\Theano\compiledir_Windows-10-10.0.19041-SP0-Intel64_Family_6_Model_60_Stepping_3_GenuineIntel-3.8.8-64\tmpom5gfujq\m49317b78e657a6e50b4a70b6f485f9b198c08562c7599af663e1f62ca091d93b.pyd" "C:\Users\diese\AppData\Local\Theano\compiledir_Windows-10-10.0.19041-SP0-Intel64_Family_6_Model_60_Stepping_3_GenuineIntel-3.8.8-64\tmpom5gfujq\mod.cpp" -lpython38In file included from C:\Users\diese\AppData\Local\Theano\compiledir_Windows-10-10.0.19041-SP0-Intel64_Family_6_Model_60_Stepping_3_GenuineIntel-3.8.8-64\tmpom5gfujq\mod.cpp:8:
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/special.h: In function 'double sp::special::digamma::psi(const double&)':
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/special.h:134:10: error: 'M_PI' was not declared in this scope
     y = -M_PI / tan(M_PI * r);
          ^~~~
In file included from C:\Users\diese\AppData\Local\Theano\compiledir_Windows-10-10.0.19041-SP0-Intel64_Family_6_Model_60_Stepping_3_GenuineIntel-3.8.8-64\tmpom5gfujq\mod.cpp:12:
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/flux.h: In function 'void sp::flux::computerT(int, sp::utils::RowVector<Scalar, -1>&)':
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/flux.h:28:10: error: 'M_PI' was not declared in this scope
   amp0 = M_PI;
          ^~~~
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/flux.h: In function 'void sp::flux::amp(int, Eigen::MatrixBase<Derived>&)':
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/flux.h:202:18: error: 'M_PI' was not declared in this scope
   M /= (2 * sqrt(M_PI));
                  ^~~~
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/flux.h: In function 'void sp::flux::computeA1(int, Eigen::SparseMatrix<Scalar>&)':
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/flux.h:248:28: error: 'M_PI' was not declared in this scope
   Scalar norm = 2.0 / sqrt(M_PI);
                            ^~~~
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/flux.h: In member function 'void sp::flux::LimbDark<Scalar>::computeU1()':
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/flux.h:332:30: error: 'M_PI' was not declared in this scope
     Scalar norm = 2.0 / sqrt(M_PI);
                              ^~~~
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/flux.h: In member function 'void sp::flux::LimbDark<Scalar>::computerTA1L(sp::utils::Vector<Scalar, 2>&, ROWVECTOR&)':
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/flux.h:512:17: error: 'M_PI' was not declared in this scope
     p *= norm * M_PI;
                 ^~~~
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/flux.h: In member function 'void sp::flux::LimbDark<Scalar>::computerTA1L(sp::utils::Vector<Scalar, 2>&, ROWVECTOR&, VECTOR&)':
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/flux.h:540:17: error: 'M_PI' was not declared in this scope
     p *= norm * M_PI;
                 ^~~~
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/flux.h: In member function 'void sp::flux::LimbDark<Scalar>::computeL(sp::utils::Vector<Scalar, 2>&, ROWMATRIX&)':
D:\Software\Programming\Anaconda\lib\site-packages\starry_process\ops\include/flux.h:574:17: error: 'M_PI' was not declared in this scope
     p *= norm * M_PI;
                 ^~~~
At global scope:
cc1plus.exe: warning: unrecognized command line option '-Wno-c++11-narrowing

What am i doing wrong?

Sample from the conditional Ylm distribution

We need to figure out how to sample from the distribution over spherical harmonic coefficients conditioned on a measurement of the flux. If we are not (1) normalizing the light curve and (2) marginalizing over inclination, this can be done by solving the least squares problem, in which case the covariance of the desired distribution is

(A^T . C^-1 . A + L^-1)^-1

where A is the design matrix that transforms spherical harmonics to flux, C is the flux covariance, and L is the prior on the spherical harmonics (equal to the covariance of the starry-process GP).

But if we normalize the light curve, the transformation from sph. harms. to flux is nonlinear, so this won't work. If we marginalize over inclination it also won't work because the transformation is stochastic!

@dfm Any ideas if this is even possible to do with linear algebra? Or would it require numerical sampling?

# TODO
if self._marginalize_over_inclination:
raise NotImplementedError(
"Not yet implemented when marginalizing over inclination."
)
# TODO
if self._normalized:
raise NotImplementedError(
"Not yet implemented when the flux is normalized."
)
# TODO
if self._time_variable:
raise NotImplementedError(
"Not yet implemented for time-variable maps."
)

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.