Initial guess selection with `multicube`: a multiple component case

This example picks up where example.ipynb notebook ended. This notebook illustrates an example usage the SubCube routines to select an optimal starting guess for the case of double velocity component model.



In [1]:

    
import numpy as np
%matplotlib inline
import matplotlib.pylab as plt
import pyspeckit
from multicube.subcube import SubCube
from multicube.astro_toolbox import make_test_cube, get_ncores
from IPython.utils import io
import warnings
warnings.filterwarnings('ignore')

Let's first make a test FITS file, containing a mixture of synthetic signal with some noise put on top of it. The created spectral cube will be 10x10 pixels wide in the plain of sky and 300 pixels "long" along its spectral axis.



In [2]:

    
make_test_cube((300,10,10), outfile='foo.fits', sigma=(10,5))
sc = SubCube('foo.fits')

To make things interesting, let's introduce a radial velocity gradient in our cube along with a second, weaker component.



In [3]:

    
# TODO: move this to astro_toolbox.py
#       as a general synthetic cube generator routine
def tinker_ppv(arr):
    scale_roll = 15
    rel_shift  = 30
    rel_str    = 5
    shifted_component = np.roll(arr, rel_shift) / rel_str
    for y,x in np.ndindex(arr.shape[1:]):
        roll  = np.sqrt((x-5)**2 + (y-5)**2) * scale_roll
        arr[:,y,x] = np.roll(arr[:,y,x], int(roll))
    return arr + shifted_component
sc.cube = tinker_ppv(sc.cube)

This is how a sample spectrum looks like at x,y = (3,7):



In [4]:

    
sc.plot_spectrum(3,7)

For this case, all we have to do is to pass an extra triple of gaussian parameters to multicube. Similarly to a the other example, the routine then generates a model grid in the parameter space with given spacing set by finesse



In [5]:

    
npeaks = 2
sc.update_model('gaussian')
sc.specfit.fitter.npeaks = npeaks
minpars = [0.1, sc.xarr.min().value, 0.1] + [0.1, -13, 0.5]
maxpars = [2.0, sc.xarr.max().value, 2.0] + [0.5,  -5, 1.5]
finesse = [5, 10, 5] + [5, 20, 5]
sc.make_guess_grid(minpars, maxpars, finesse)









    



INFO: Selected gaussian model [multicube.subcube]
INFO: Binning the 6-dimensional parameter space into a (5, 10, 5, 5, 20, 5)-shaped grid [multicube.subcube]






    Out[5]:





array([[  0.41666667, -21.78181818,   0.41666667,   0.16666667,
        -12.61904762,   0.66666667],
       [  0.41666667, -21.78181818,   0.41666667,   0.16666667,
        -12.61904762,   0.83333333],
       [  0.41666667, -21.78181818,   0.41666667,   0.16666667,
        -12.61904762,   1.        ],
       ..., 
       [  1.68333333,   2.68181818,   1.68333333,   0.43333333,
         -5.38095238,   1.        ],
       [  1.68333333,   2.68181818,   1.68333333,   0.43333333,
         -5.38095238,   1.16666667],
       [  1.68333333,   2.68181818,   1.68333333,   0.43333333,
         -5.38095238,   1.33333333]])

This grid, stored under sc.guess_grid, can be used to generate a number of spectral models with pyspeckit, and the guesses that have the least residual rms can be selected for the whole cube:

sc.best_map stored the map between x,y pixel numbers and the numbers of corresponding best models
sc.best_model is the number of the model suited best for the pixel with the highest S/N ratio.



In [6]:

    
sc.generate_model()
sc.get_snr_map()
sc.best_guess()









    



INFO: Generating spectral models from the guess grid . . . [multicube.subcube]
INFO: Calculating residuals for generated models . . . [multicube.subcube]






    



WARNING: The available free memory might not be enough for broadcasting model grid to the spectral cube. Will iterate over all the XY pairs instead. Coffee time! [multicube.subcube]






    



INFO: Overall best model: selected #13540 [  0.42 -19.06   1.05   0.17  -9.57   0.67] [multicube.subcube]
INFO: Best model @ highest SNR: #64938 [ 0.42 -8.19  1.68  0.43 -9.95  1.17] [multicube.subcube]






    



WARNING: Selected model is best only for less than %5 of the cube, consider using the map of guesses. [multicube.subcube]

Fitting the cube best models for each x,y pixel:



In [7]:

    
sc.fiteach(fittype   = sc.fittype,
           guesses   = sc.best_guesses, 
           multicore = get_ncores(),
           errmap    = sc._rms_map,
           **sc.fiteach_args);









    



INFO: Finished fit     53 of    100 at (   3,   7) s/n=  9.6. Elapsed time is 0.3 seconds.  %53 [multicube.subcube]
INFO: Finished fit      1 of    100 at (   0,   0) s/n=  4.0. Elapsed time is 0.3 seconds.  %1 [multicube.subcube]
INFO: Finished fit     27 of    100 at (   5,   1) s/n=  4.0. Elapsed time is 0.3 seconds.  %27 [multicube.subcube]
INFO: Finished fit     89 of    100 at (   9,   5) s/n=  3.9. Elapsed time is 0.3 seconds.  %89 [multicube.subcube]
INFO: Finished fit     65 of    100 at (   6,   6) s/n= 18.6. Elapsed time is 0.3 seconds.  %65 [multicube.subcube]
INFO: Finished fit     14 of    100 at (   2,   3) s/n=  4.5. Elapsed time is 0.5 seconds.  %14 [multicube.subcube]
INFO: Finished fit      2 of    100 at (   1,   0) s/n=  3.7. Elapsed time is 0.5 seconds.  %2 [multicube.subcube]
INFO: Finished fit     54 of    100 at (   7,   3) s/n=  8.9. Elapsed time is 0.5 seconds.  %54 [multicube.subcube]
INFO: Finished fit     28 of    100 at (   1,   5) s/n=  4.0. Elapsed time is 0.6 seconds.  %28 [multicube.subcube]
INFO: Finished fit     90 of    100 at (   5,   9) s/n=  4.0. Elapsed time is 0.6 seconds.  %90 [multicube.subcube]
INFO: Finished fit     66 of    100 at (   8,   3) s/n=  4.4. Elapsed time is 0.6 seconds.  %66 [multicube.subcube]
INFO: Finished fit     40 of    100 at (   4,   5) s/n= 21.2. Elapsed time is 0.7 seconds.  %40 [multicube.subcube]
INFO: Finished fit     77 of    100 at (   7,   6) s/n= 16.5. Elapsed time is 0.7 seconds.  %77 [multicube.subcube]
INFO: Finished fit     29 of    100 at (   2,   5) s/n=  7.2. Elapsed time is 0.9 seconds.  %29 [multicube.subcube]
INFO: Finished fit     15 of    100 at (   3,   2) s/n=  4.5. Elapsed time is 0.9 seconds.  %15 [multicube.subcube]
INFO: Finished fit      3 of    100 at (   0,   1) s/n=  3.8. Elapsed time is 0.9 seconds.  %3 [multicube.subcube]
INFO: Finished fit     67 of    100 at (   3,   8) s/n=  4.4. Elapsed time is 1.0 seconds.  %67 [multicube.subcube]
INFO: Finished fit     91 of    100 at (   8,   7) s/n=  4.5. Elapsed time is 1.0 seconds.  %91 [multicube.subcube]
INFO: Finished fit     78 of    100 at (   2,   9) s/n=  4.1. Elapsed time is 1.0 seconds.  %78 [multicube.subcube]
INFO: Finished fit     41 of    100 at (   5,   4) s/n= 23.2. Elapsed time is 1.1 seconds.  %41 [multicube.subcube]
INFO: Finished fit     30 of    100 at (   5,   2) s/n=  7.2. Elapsed time is 1.2 seconds.  %30 [multicube.subcube]
INFO: Finished fit     16 of    100 at (   0,   4) s/n=  3.6. Elapsed time is 1.2 seconds.  %16 [multicube.subcube]
INFO: Finished fit     92 of    100 at (   7,   8) s/n=  4.7. Elapsed time is 1.2 seconds.  %92 [multicube.subcube]
INFO: Finished fit      4 of    100 at (   1,   1) s/n=  3.8. Elapsed time is 1.3 seconds.  %4 [multicube.subcube]
INFO: Finished fit     79 of    100 at (   9,   2) s/n=  3.9. Elapsed time is 1.3 seconds.  %79 [multicube.subcube]
INFO: Finished fit     55 of    100 at (   5,   6) s/n= 21.8. Elapsed time is 1.4 seconds.  %55 [multicube.subcube]
INFO: Finished fit     17 of    100 at (   4,   0) s/n=  4.0. Elapsed time is 1.5 seconds.  %17 [multicube.subcube]
INFO: Finished fit      5 of    100 at (   0,   2) s/n=  3.8. Elapsed time is 1.5 seconds.  %5 [multicube.subcube]
INFO: Finished fit     31 of    100 at (   4,   4) s/n= 22.1. Elapsed time is 1.6 seconds.  %31 [multicube.subcube]
INFO: Finished fit     93 of    100 at (   9,   6) s/n=  3.9. Elapsed time is 1.6 seconds.  %93 [multicube.subcube]
INFO: Finished fit     80 of    100 at (   5,   8) s/n=  7.2. Elapsed time is 1.6 seconds.  %80 [multicube.subcube]
INFO: Finished fit     68 of    100 at (   5,   7) s/n= 17.4. Elapsed time is 1.6 seconds.  %68 [multicube.subcube]
INFO: Finished fit     42 of    100 at (   6,   3) s/n= 19.1. Elapsed time is 1.7 seconds.  %42 [multicube.subcube]
INFO: Finished fit      6 of    100 at (   2,   0) s/n=  3.9. Elapsed time is 1.8 seconds.  %6 [multicube.subcube]
INFO: Finished fit     18 of    100 at (   1,   4) s/n=  4.0. Elapsed time is 1.8 seconds.  %18 [multicube.subcube]
INFO: Finished fit     56 of    100 at (   6,   5) s/n= 21.0. Elapsed time is 1.8 seconds.  %56 [multicube.subcube]
INFO: Finished fit     94 of    100 at (   6,   9) s/n=  4.0. Elapsed time is 1.9 seconds.  %94 [multicube.subcube]
INFO: Finished fit     81 of    100 at (   8,   5) s/n=  7.0. Elapsed time is 1.9 seconds.  %81 [multicube.subcube]
INFO: Finished fit     19 of    100 at (   4,   1) s/n=  4.0. Elapsed time is 2.0 seconds.  %19 [multicube.subcube]
INFO: Finished fit      7 of    100 at (   1,   2) s/n=  3.9. Elapsed time is 2.0 seconds.  %7 [multicube.subcube]
INFO: Finished fit     95 of    100 at (   8,   8) s/n=  4.0. Elapsed time is 2.1 seconds.  %95 [multicube.subcube]
INFO: Finished fit     82 of    100 at (   3,   9) s/n=  3.8. Elapsed time is 2.2 seconds.  %82 [multicube.subcube]
INFO: Finished fit     69 of    100 at (   7,   5) s/n= 18.3. Elapsed time is 2.2 seconds.  %69 [multicube.subcube]
INFO: Finished fit     43 of    100 at (   3,   6) s/n= 19.0. Elapsed time is 2.3 seconds.  %43 [multicube.subcube]
INFO: Finished fit     57 of    100 at (   0,   8) s/n=  3.7. Elapsed time is 2.3 seconds.  %57 [multicube.subcube]
INFO: Finished fit      8 of    100 at (   2,   1) s/n=  3.7. Elapsed time is 2.4 seconds.  %8 [multicube.subcube]
INFO: Finished fit     20 of    100 at (   3,   3) s/n= 10.2. Elapsed time is 2.4 seconds.  %20 [multicube.subcube]
INFO: Finished fit     32 of    100 at (   3,   5) s/n= 19.9. Elapsed time is 2.4 seconds.  %32 [multicube.subcube]
INFO: Finished fit     96 of    100 at (   7,   9) s/n=  3.9. Elapsed time is 2.4 seconds.  %96 [multicube.subcube]
INFO: Finished fit     70 of    100 at (   4,   8) s/n=  6.1. Elapsed time is 2.4 seconds.  %70 [multicube.subcube]
INFO: Finished fit     83 of    100 at (   9,   3) s/n=  3.9. Elapsed time is 2.5 seconds.  %83 [multicube.subcube]
INFO: Finished fit     44 of    100 at (   7,   0) s/n=  3.7. Elapsed time is 2.5 seconds.  %44 [multicube.subcube]
INFO: Finished fit     58 of    100 at (   8,   0) s/n=  3.5. Elapsed time is 2.5 seconds.  %58 [multicube.subcube]
INFO: Finished fit     97 of    100 at (   9,   7) s/n=  3.9. Elapsed time is 2.6 seconds.  %97 [multicube.subcube]
INFO: Finished fit     21 of    100 at (   4,   2) s/n=  6.3. Elapsed time is 2.7 seconds.  %21 [multicube.subcube]
INFO: Finished fit     84 of    100 at (   4,   9) s/n=  3.8. Elapsed time is 2.8 seconds.  %84 [multicube.subcube]
INFO: Finished fit      9 of    100 at (   2,   2) s/n=  3.8. Elapsed time is 2.8 seconds.  %9 [multicube.subcube]
INFO: Finished fit     45 of    100 at (   0,   7) s/n=  3.8. Elapsed time is 2.8 seconds.  %45 [multicube.subcube]
INFO: Finished fit     71 of    100 at (   8,   4) s/n=  6.3. Elapsed time is 2.8 seconds.  %71 [multicube.subcube]
INFO: Finished fit     98 of    100 at (   9,   8) s/n=  4.0. Elapsed time is 2.8 seconds.  %98 [multicube.subcube]
INFO: Finished fit     22 of    100 at (   2,   4) s/n=  6.8. Elapsed time is 3.0 seconds.  %22 [multicube.subcube]
INFO: Finished fit     33 of    100 at (   5,   3) s/n= 20.0. Elapsed time is 3.0 seconds.  %33 [multicube.subcube]
INFO: Finished fit     10 of    100 at (   3,   0) s/n=  3.9. Elapsed time is 3.1 seconds.  %10 [multicube.subcube]
INFO: Finished fit     85 of    100 at (   9,   4) s/n=  3.8. Elapsed time is 3.1 seconds.  %85 [multicube.subcube]
INFO: Finished fit     99 of    100 at (   8,   9) s/n=  3.9. Elapsed time is 3.1 seconds.  %99 [multicube.subcube]
INFO: Finished fit     72 of    100 at (   9,   0) s/n=  3.9. Elapsed time is 3.1 seconds.  %72 [multicube.subcube]
INFO: Finished fit     46 of    100 at (   5,   5) s/n= 25.1. Elapsed time is 3.2 seconds.  %46 [multicube.subcube]
INFO: Finished fit     23 of    100 at (   0,   5) s/n=  3.9. Elapsed time is 3.2 seconds.  %23 [multicube.subcube]
INFO: Finished fit     59 of    100 at (   4,   7) s/n= 18.0. Elapsed time is 3.2 seconds.  %59 [multicube.subcube]
INFO: Finished fit     86 of    100 at (   7,   7) s/n= 10.2. Elapsed time is 3.3 seconds.  %86 [multicube.subcube]
INFO: Finished fit     34 of    100 at (   0,   6) s/n=  3.8. Elapsed time is 3.3 seconds.  %34 [multicube.subcube]
INFO: Finished fit     11 of    100 at (   0,   3) s/n=  3.9. Elapsed time is 3.4 seconds.  %11 [multicube.subcube]
INFO: Finished fit     73 of    100 at (   0,   9) s/n=  3.7. Elapsed time is 3.4 seconds.  %73 [multicube.subcube]
INFO: Finished fit    100 of    100 at (   9,   9) s/n=  3.8. Elapsed time is 3.4 seconds.  %100 [multicube.subcube]
INFO: Finished fit     87 of    100 at (   6,   8) s/n=  7.0. Elapsed time is 3.5 seconds.  %87 [multicube.subcube]
INFO: Finished fit     47 of    100 at (   1,   7) s/n=  3.8. Elapsed time is 3.5 seconds.  %47 [multicube.subcube]
INFO: Finished fit     35 of    100 at (   6,   0) s/n=  3.9. Elapsed time is 3.6 seconds.  %35 [multicube.subcube]
INFO: Finished fit     60 of    100 at (   1,   8) s/n=  3.8. Elapsed time is 3.6 seconds.  %60 [multicube.subcube]
INFO: Finished fit     12 of    100 at (   1,   3) s/n=  3.9. Elapsed time is 3.7 seconds.  %12 [multicube.subcube]
INFO: Finished fit     24 of    100 at (   3,   4) s/n= 18.6. Elapsed time is 3.7 seconds.  %24 [multicube.subcube]
INFO: Finished fit     88 of    100 at (   8,   6) s/n=  6.7. Elapsed time is 3.8 seconds.  %88 [multicube.subcube]
INFO: Finished fit     48 of    100 at (   7,   1) s/n=  3.9. Elapsed time is 3.8 seconds.  %48 [multicube.subcube]
INFO: Finished fit     74 of    100 at (   9,   1) s/n=  4.0. Elapsed time is 3.8 seconds.  %74 [multicube.subcube]
INFO: Finished fit     61 of    100 at (   8,   1) s/n=  3.9. Elapsed time is 3.8 seconds.  %61 [multicube.subcube]
INFO: Finished fit     13 of    100 at (   3,   1) s/n=  3.9. Elapsed time is 3.9 seconds.  %13 [multicube.subcube]
INFO: Finished fit     36 of    100 at (   1,   6) s/n=  3.9. Elapsed time is 3.9 seconds.  %36 [multicube.subcube]
INFO: Finished fit     49 of    100 at (   6,   4) s/n= 20.7. Elapsed time is 4.0 seconds.  %49 [multicube.subcube]
INFO: Finished fit     25 of    100 at (   5,   0) s/n=  3.9. Elapsed time is 4.0 seconds.  %25 [multicube.subcube]
INFO: Finished fit     75 of    100 at (   1,   9) s/n=  3.8. Elapsed time is 4.0 seconds.  %75 [multicube.subcube]
INFO: Finished fit     62 of    100 at (   7,   4) s/n= 17.2. Elapsed time is 4.2 seconds.  %62 [multicube.subcube]
INFO: Finished fit     37 of    100 at (   6,   1) s/n=  4.0. Elapsed time is 4.3 seconds.  %37 [multicube.subcube]
INFO: Finished fit     50 of    100 at (   4,   6) s/n= 20.7. Elapsed time is 4.3 seconds.  %50 [multicube.subcube]
INFO: Finished fit     63 of    100 at (   8,   2) s/n=  3.9. Elapsed time is 4.4 seconds.  %63 [multicube.subcube]
INFO: Finished fit     76 of    100 at (   6,   7) s/n= 17.7. Elapsed time is 4.5 seconds.  %76 [multicube.subcube]
INFO: Finished fit     38 of    100 at (   2,   6) s/n=  6.5. Elapsed time is 4.5 seconds.  %38 [multicube.subcube]
INFO: Finished fit     26 of    100 at (   4,   3) s/n= 18.8. Elapsed time is 4.5 seconds.  %26 [multicube.subcube]
INFO: Finished fit     51 of    100 at (   7,   2) s/n=  4.6. Elapsed time is 4.5 seconds.  %51 [multicube.subcube]
INFO: Finished fit     64 of    100 at (   2,   8) s/n=  3.9. Elapsed time is 4.6 seconds.  %64 [multicube.subcube]
INFO: Finished fit     39 of    100 at (   6,   2) s/n=  6.3. Elapsed time is 4.7 seconds.  %39 [multicube.subcube]
INFO: Finished fit     52 of    100 at (   2,   7) s/n=  4.6. Elapsed time is 4.7 seconds.  %52 [multicube.subcube]
INFO: Finished final fit 100.  Elapsed time was 4.8 seconds [multicube.subcube]



In [8]:

    
sc.get_modelcube()
import matplotlib.gridspec as gridspec
plt.figure(figsize=(10,10))
gs = gridspec.GridSpec(*sc.cube.shape[1:])
gs.update(wspace=0., hspace=0.)
for ii in range(np.prod(sc.cube.shape[1:])):
    xy = np.unravel_index(ii, dims=sc.cube.shape[1:])
    ax = plt.subplot(gs[ii])
    sc.plot_spectrum(*xy, axis=ax)
    ax.set_xticklabels([])
    ax.set_yticklabels([])
    ax.set_xlabel('')
    ax.set_ylabel('')
    ax.tick_params(axis=u'both', which=u'both',length=0)
    ax.set_ylim(-0.15, 1.3)
    ax.plot(sc.xarr.value, sc.model_grid[sc._best_map[xy]], label='Initial guess\nat (x,y)={}'.format(xy))
    ax.plot(sc.xarr.value, sc._modelcube[:, xy[1], xy[0]], label='Final fit\nat (x,y)={}'.format(xy))

Initial guess selection with multicube: a multiple component case

Fitting the cube best models for each x,y pixel:

Initial guess selection with `multicube`: a multiple component case