wip

src/core.py

@@ -471,7 +471,8 @@ class RandomResponse(trace.FrequencyResponse):
 
     def evaluate(self, freqs):
         n = freqs.size
-        return 1.0 + freqs*(self._rstate.normal(scale=self.scale, size=n) + 
+        return 1.0 + freqs*(
+            self._rstate.normal(scale=self.scale, size=n) +
             0.0J * self._rstate.normal(scale=self.scale, size=n))
 
 
@@ -562,13 +563,13 @@ class SyntheticTest(Object):
                 randomresponse.set_random_state(self._rstate)
 
                 tr = tr.transfer(tfade=tfade, freqlimits=freqlimits)
-                tr2 = tr2.transfer(tfade=tfade, freqlimits=freqlimits,
+                tr2 = tr2.transfer(
+                    tfade=tfade,
+                    freqlimits=freqlimits,
                     transfer_function=randomresponse)
 
                 tr.chop(tmin, tmax)
                 tr2.chop(tmin, tmax)
-
-                #trace.snuffle([tr, tr2])
                 return tr2
 
         return None
@@ -978,6 +979,7 @@ def solve(problem,
     isbad_mask = None
     accept_sum = num.zeros(1 + problem.nbootstrap, dtype=num.int)
     accept_hist = num.zeros(niter, dtype=num.int)
+    pnames = [p.name for p in problem.parameters]
 
     while iiter < niter:
 
@@ -1015,19 +1017,36 @@ def solve(problem,
                     if sampler_distribution == 'multivariate_normal':
                         ntries_sample = 0
 
+                        ntry = 0
+                        ok_mask_sum = num.zeros(npar, dtype=num.int)
                         while True:
                             ntries_sample += 1
                             vs = num.random.multivariate_normal(
                                 xb, factor*covs[jchoice])
 
-                            if (num.all(xbounds[:, 0] <= vs) and
-                                    num.all(vs <= xbounds[:, 1])):
+                            ok_mask = num.logical_and(
+                                xbounds[:, 0] <= vs, vs <= xbounds[:, 1])
+
+                            if num.all(ok_mask):
                                 break
 
+                            ok_mask_sum += ok_mask
+
+                            if ntry > 1000:
+                                raise GrondError(
+                                    'failed to produce a suitable candidate '
+                                    'sample from multivariate normal '
+                                    'distribution, (%s)' %
+                                    ', '.join('%s:%i' % xx for xx in
+                                              zip(pnames, ok_mask_sum)))
+
+                            ntry += 1
+
                         x = vs.tolist()
 
                     if sampler_distribution == 'normal':
                         for i in xrange(npar):
+                            ntry = 0
                             while True:
                                 v = num.random.normal(
                                     xb[i], math.sqrt(factor)*sbx[i])
@@ -1035,6 +1054,14 @@ def solve(problem,
                                 if xbounds[i, 0] <= v and v <= xbounds[i, 1]:
                                     break
 
+                                if ntry > 1000:
+                                    raise GrondError(
+                                        'failed to produce a suitable '
+                                        'candidate sample from normal '
+                                        'distribution')
+
+                                ntry += 1
+
                             x.append(v)
 
                 try:
@@ -1127,8 +1154,7 @@ def solve(problem,
                      'G best'))
 
                 for (pname, mbv, sbv, mgv, sgv, bgv) in zip(
-                        [p.name for p in problem.parameters],
-                        mbx, sbx, mgx, sgx, bgx):
+                        pnames, mbx, sbx, mgx, sgx, bgx):
 
                     lines.append(
                         '%-15s %15.4g %15.4g %15.4g %15.4g %15.4g' %
@@ -1190,9 +1216,13 @@ def bootstrap_outliers(problem, misfits, std_factor=1.0):
 def forward(rundir_or_config_path, event_names=None):
 
     if os.path.isdir(rundir_or_config_path):
-        config = guts.load(filename=op.join(rundir_or_config_path, 'config.yaml'))
+        rundir = rundir_or_config_path
+        config = guts.load(
+            filename=op.join(rundir, 'config.yaml'))
+
         config.set_basepath(rundir)
-        problem, xs, misfits = load_problem_info_and_data(rundir, subset='harvest')
+        problem, xs, misfits = load_problem_info_and_data(
+            rundir, subset='harvest')
 
         gms = problem.global_misfits(misfits)
         ibest = num.argmin(gms)

src/plot.py

@@ -380,7 +380,8 @@ def draw_jointpar_figures(
     iselected = 0
     for ipar in xrange(problem.ncombined):
         par = problem.combined[ipar]
-        if exclude and par.name in exclude or include and par.name not in include:
+        if exclude and par.name in exclude or \
+                include and par.name not in include:
             continue
 
         smap[iselected] = ipar
@@ -513,11 +514,6 @@ def draw_jointpar_figures(
                 xpar.scaled(fx), ypar.scaled(fy), 's',
                 mew=1.5, ms=5, color=ref_color_light, mec=ref_color)
 
-    #for jfig, figs_row in enumerate(figs):
-    #    for ifig, fig in enumerate(figs_row):
-    #        if fig is not None:
-    #            fig.savefig('jointpar-%i-%i.pdf' % (jfig, ifig))
-
 
 def draw_solution_figure(
         model, plt, misfit_cutoff=None, beachball_type='full'):
@@ -862,9 +858,8 @@ def draw_fits_figures(ds, model, plt):
 
     dtraces = []
     for target, result in zip(problem.targets, results):
+        print target.misfit_config.domain
         if result is None:
-            print 'xxx'
-            print target
             dtraces.append(None)
             continue
 
@@ -893,6 +888,10 @@ def draw_fits_figures(ds, model, plt):
         dtraces.append(dtrace)
         all_syn_trs.append(result.processed_syn)
 
+    if not all_syn_trs:
+        logger.warn('no traces to show')
+        return
+
     amin, amax = trace.minmax(all_syn_trs, lambda tr: None)[None]
 
     dmin, dmax = trace.minmax(
@@ -982,7 +981,6 @@ def draw_fits_figures(ds, model, plt):
                 if (iy, ix) not in frame_to_target:
                     continue
 
-
                 ixx = ix/nxmax
                 iyy = iy/nymax
                 if (iyy, ixx) not in figures: