merge

examples/multi_instrument/OPTIONS

@@ -2,12 +2,12 @@
 # Put comments at the top, or at the end of the line.
 2	# Number of particles
 5000	# new level interval
-2000	# save interval
-50	# threadSteps: number of steps each thread does independently before communication
+3000	# save interval
+100	# threadSteps: number of steps each thread does independently before communication
 0	# maximum number of levels
 10	# Backtracking scale length (lambda)
 100	# Strength of effect to force histogram to equal push (beta)
-10000	# Maximum number of saves (0 = infinite)
+20000	# Maximum number of saves (0 = infinite)
     # (optional) samples file
     # (optional) sample_info file
     # (optional) levels file

examples/multi_instrument/kima_setup.cpp

@@ -24,6 +24,9 @@ RVmodel::RVmodel():fix(true),npmax(1)
     // set the prior for the between-instrument offsets
     offsets_prior = new Uniform(-RVspan, RVspan);
 
+    // note: there is one extra white noise (jitter) parameter per instrument
+    // they all share the same prior, defined as Jprior
+
     save_setup();
 }
 

pykima/display.py

@@ -616,21 +616,26 @@ def make_plot5(self, show=True, save=False):
 
         # available_etas = [v for v in dir(self) if v.startswith('eta')]
         available_etas = ['eta1', 'eta2', 'eta3', 'eta4']
-        labels = [r'$s$'] + [r'$\eta_%d$' % (i+1) for i,_ in enumerate(available_etas)]
-        units = ['m/s', 'm/s', 'days', 'days', None]
+        labels = [r'$s$']*self.n_jitters
+        labels += [r'$\eta_%d$' % (i+1) for i,_ in enumerate(available_etas)]
+        units = ['m/s']*self.n_jitters + ['m/s', 'days', 'days', None]
         xlabels = []
         for label, unit in zip(labels, units):
             xlabels.append(label + ' (%s)' % unit 
                                 if unit is not None else label)
 
         ### all Np together
-        variables = [self.extra_sigma]
+        if self.multi:
+            variables = list(self.extra_sigma.T)
+        else:
+            variables = [self.extra_sigma]
+
         for eta in available_etas:
             variables.append(getattr(self, eta))
 
         self.post_samples = np.vstack(variables).T
 
-        ranges = [1.]*(len(available_etas)+1)
+        ranges = [1.]*(len(available_etas) + self.extra_sigma.shape[1])
         # ranges[3] = (self.pmin, self.pmax)
 
         c = corner.corner
@@ -640,7 +645,7 @@ def make_plot5(self, show=True, save=False):
                             # fill_contours=True, smooth=True,
                             # contourf_kwargs={'cmap':plt.get_cmap('afmhot'), 'colors':None},
                             hexbin_kwargs={'cmap':plt.get_cmap('afmhot_r'), 'bins':'log'},
-                            hist_kwargs={'normed':True}, 
+                            hist_kwargs={'density':True},
                             range=ranges, data_kwargs={'alpha':1},
                             )
         except AssertionError as exc:

src/RVmodel.cpp

@@ -81,11 +81,16 @@ void RVmodel::from_prior(RNG& rng)
 void RVmodel::calculate_C()
 {
     // Get the data
-    const vector<double>& t = Data::get_instance().get_t();
-    const vector<double>& sig = Data::get_instance().get_sig();
-    int N = Data::get_instance().get_t().size();
+    auto data = Data::get_instance();
+    const vector<double>& t = data.get_t();
+    const vector<double>& sig = data.get_sig();
+    const vector<int>& obsi = data.get_obsi();
+    int N = data.N();
+    double jit;
 
-    // auto begin = std::chrono::high_resolution_clock::now();  // start timing
+    #if TIMING
+    auto begin = std::chrono::high_resolution_clock::now();  // start timing
+    #endif
 
     for(size_t i=0; i<N; i++)
     {
@@ -95,14 +100,30 @@ void RVmodel::calculate_C()
                         -2.0*pow(sin(M_PI*(t[i] - t[j])/eta3)/eta4, 2) );
 
             if(i==j)
-                C(i, j) += sig[i]*sig[i] + extra_sigma*extra_sigma;
+            {
+                if (multi_instrument)
+                {
+                    jit = jitters[obsi[i]-1];
+                    C(i, j) += sig[i]*sig[i] + jit*jit;
+                }
+                else
+                {
+                    C(i, j) += sig[i]*sig[i] + extra_sigma*extra_sigma;
+                }
+            }
             else
+            {
                 C(j, i) = C(i, j);
+            }
         }
     }
 
-    // auto end = std::chrono::high_resolution_clock::now();
-    // cout << "GP: " << std::chrono::duration_cast<std::chrono::nanoseconds>(end-begin).count() << " ns" << "\t"; // << std::endl;
+    #if TIMING
+    auto end = std::chrono::high_resolution_clock::now();
+    cout << "GP build matrix: ";
+    cout << std::chrono::duration_cast<std::chrono::nanoseconds>(end-begin).count();
+    cout << " ns" << "\t"; // << std::endl;
+    #endif
 }
 
 void RVmodel::calculate_mu()
@@ -238,7 +259,15 @@ double RVmodel::perturb(RNG& rng)
         }
         else if(rng.rand() <= 0.5)
         {
-            Jprior->perturb(extra_sigma, rng);
+            if(multi_instrument)
+            {
+                for(int i=0; i<jitters.size(); i++)
+                    Jprior->perturb(jitters[i], rng);
+            }
+            else
+            {
+                Jprior->perturb(extra_sigma, rng);
+            }
             calculate_C();
         }
         else