def _run_interface(self, runtime):
        nii1 = nb.load(self.inputs.volume1)
        nii2 = nb.load(self.inputs.volume2)

        origdata1 = np.logical_not(np.logical_or(nii1.get_data() == 0, np.isnan(nii1.get_data())))
        origdata2 = np.logical_not(np.logical_or(nii2.get_data() == 0, np.isnan(nii2.get_data())))

        if isdefined(self.inputs.mask_volume):
            maskdata = nb.load(self.inputs.mask_volume).get_data()
            maskdata = np.logical_not(np.logical_or(maskdata == 0, np.isnan(maskdata)))
            origdata1 = np.logical_and(maskdata, origdata1)
            origdata2 = np.logical_and(maskdata, origdata2)

        for method in ("dice", "jaccard"):
            setattr(self, "_" + method, self._bool_vec_dissimilarity(origdata1, origdata2, method=method))

        self._volume = int(origdata1.sum() - origdata2.sum())

        both_data = np.zeros(origdata1.shape)
        both_data[origdata1] = 1
        both_data[origdata2] += 2

        nb.save(nb.Nifti1Image(both_data, nii1.get_affine(), nii1.get_header()), self.inputs.out_file)

        return runtime

def evolution_of_votes_singleMP(dates, votes, wa_all, wa_party, name, asciiname):
    if not do_plots: return
    f = plt.figure(figsize=figsize_long)
    f.suptitle(u'Гласове и отсъствия на %s през годините.'%name)
    absences = f.add_subplot(3,1,3)
    with_all = f.add_subplot(3,1,1, sharex=absences)
    with_party = f.add_subplot(3,1,2, sharex=absences)

    all_votes_no_abs = np.sum(votes[:,:3], 1)
    all_votes = np.sum(votes, 1)
    mask_no_abs = np.logical_not(all_votes_no_abs)
    mask = np.logical_not(all_votes)
    with_all_array = np.ma.masked_array(100*wa_all[:,0], mask=mask_no_abs)/all_votes_no_abs
    with_party_array = np.ma.masked_array(100*wa_party[:,0], mask=mask_no_abs)/all_votes_no_abs
    absences_array = np.ma.masked_array(100*votes[:,3], mask=mask)/all_votes

    with_all.plot(dates, with_all_array, '.-', alpha=0.3, linewidth=0.1)
    with_all.legend([u'% съгласие с мнозинството (без отсъствия)'])
    with_party.plot(dates, with_party_array, '.-', alpha=0.3, linewidth=0.1)
    with_party.legend([u'% съгласие с партията (без отсъствия)'])
    absences.plot(dates, absences_array, '.-', alpha=0.3, linewidth=0.1)
    absences.legend([u'% отсъствия'])

    with_all.set_yticks([25, 50, 75])
    with_party.set_yticks([25, 50, 75])
    absences.set_yticks([25, 50, 75])
    with_all.set_ylim(0, 100)
    with_party.set_ylim(0, 100)
    absences.set_ylim(0, 100)
    absences.set_xlim(dates[0], dates[-1])
    f.autofmt_xdate()
    f.savefig('generated_html/vote_evol_%s.png'%asciiname)
    plt.close()

    def __dectree_train(self, X, Y, L, R, F, T, next, depth, minParent, maxDepth, minScore, nFeatures):
        """
        This is a recursive helper method that recusively trains the decision tree. Used in:
            train

        TODO:
            compare for numerical tolerance
        """
        n,d = mat(X).shape

        # check leaf conditions...
        if n < minParent or depth >= maxDepth or np.var(Y) < minScore:
            assert n != 0, ('TreeRegress.__dectree_train: tried to create size zero node')
            return self.__output_leaf(Y, n, L, R, F, T, next)

        best_val = np.inf
        best_feat = -1
        try_feat = np.random.permutation(d)

        # ...otherwise, search over (allowed) features
        for i_feat in try_feat[0:nFeatures]:
            dsorted = arr(np.sort(X[:,i_feat].T)).ravel()                       # sort data...
            pi = np.argsort(X[:,i_feat].T)                                      # ...get sorted indices...
            tsorted = Y[pi].ravel()                                             # ...and sort targets by feature ID
            can_split = np.append(arr(dsorted[:-1] != dsorted[1:]), 0)          # which indices are valid split points?

            if not np.any(can_split):          # no way to split on this feature?
                continue

            # find min weighted variance among split points
            val,idx = self.__min_weighted_var(tsorted, can_split, n)

            # save best feature and split point found so far
            if val < best_val:
                best_val = val
                best_feat = i_feat
                best_thresh = (dsorted[idx] + dsorted[idx + 1]) / 2

        # if no split possible, output leaf (prediction) node
        if best_feat == -1:         
            return self.__output_leaf(Y, n, L, R, F, T, next)

        # split data on feature i_feat, value (tsorted[idx] + tsorted[idx + 1]) / 2
        F[next] = best_feat
        T[next] = best_thresh
        go_left = X[:,F[next]] < T[next]
        my_idx = next
        next += 1

        # recur left
        L[my_idx] = next    
        L,R,F,T,next = self.__dectree_train(X[go_left,:], Y[go_left], L, R, F, T, 
            next, depth + 1, minParent, maxDepth, minScore, nFeatures)

        # recur right
        R[my_idx] = next    
        L,R,F,T,next = self.__dectree_train(X[np.logical_not(go_left),:], Y[np.logical_not(go_left)], L, R, F, T, 
            next, depth + 1, minParent, maxDepth, minScore, nFeatures)

        return (L,R,F,T,next)

def _extrapolate_out_mask(data, mask, iterations=1):
    """ Extrapolate values outside of the mask.
    """
    if iterations > 1:
        data, mask = _extrapolate_out_mask(data, mask,
                                          iterations=iterations - 1)
    new_mask = ndimage.binary_dilation(mask)
    larger_mask = np.zeros(np.array(mask.shape) + 2, dtype=np.bool)
    larger_mask[1:-1, 1:-1, 1:-1] = mask
    # Use nans as missing value: ugly
    masked_data = np.zeros(larger_mask.shape + data.shape[3:])
    masked_data[1:-1, 1:-1, 1:-1] = data.copy()
    masked_data[np.logical_not(larger_mask)] = np.nan
    outer_shell = larger_mask.copy()
    outer_shell[1:-1, 1:-1, 1:-1] = np.logical_xor(new_mask, mask)
    outer_shell_x, outer_shell_y, outer_shell_z = np.where(outer_shell)
    extrapolation = list()
    for i, j, k in [(1, 0, 0), (-1, 0, 0), 
                    (0, 1, 0), (0, -1, 0),
                    (0, 0, 1), (0, 0, -1)]:
        this_x = outer_shell_x + i
        this_y = outer_shell_y + j
        this_z = outer_shell_z + k
        extrapolation.append(masked_data[this_x, this_y, this_z])

    extrapolation = np.array(extrapolation)
    extrapolation = (np.nansum(extrapolation, axis=0)
                     / np.sum(np.isfinite(extrapolation), axis=0))
    extrapolation[np.logical_not(np.isfinite(extrapolation))] = 0
    new_data = np.zeros_like(masked_data)
    new_data[outer_shell] = extrapolation
    new_data[larger_mask] = masked_data[larger_mask]
    return new_data[1:-1, 1:-1, 1:-1], new_mask

def test_multilabel_accuracy_score_subset_accuracy():
    # Dense label indicator matrix format
    y1 = np.array([[0, 1, 1], [1, 0, 1]])
    y2 = np.array([[0, 0, 1], [1, 0, 1]])

    assert_equal(accuracy_score(y1, y2), 0.5)
    assert_equal(accuracy_score(y1, y1), 1)
    assert_equal(accuracy_score(y2, y2), 1)
    assert_equal(accuracy_score(y2, np.logical_not(y2)), 0)
    assert_equal(accuracy_score(y1, np.logical_not(y1)), 0)
    assert_equal(accuracy_score(y1, np.zeros(y1.shape)), 0)
    assert_equal(accuracy_score(y2, np.zeros(y1.shape)), 0)

    with ignore_warnings():  # sequence of sequences is deprecated
        # List of tuple of label
        y1 = [(1, 2,), (0, 2,)]
        y2 = [(2,), (0, 2,)]

        assert_equal(accuracy_score(y1, y2), 0.5)
        assert_equal(accuracy_score(y1, y1), 1)
        assert_equal(accuracy_score(y2, y2), 1)
        assert_equal(accuracy_score(y2, [(), ()]), 0)
        assert_equal(accuracy_score(y1, y2, normalize=False), 1)
        assert_equal(accuracy_score(y1, y1, normalize=False), 2)
        assert_equal(accuracy_score(y2, y2, normalize=False), 2)
        assert_equal(accuracy_score(y2, [(), ()], normalize=False), 0)

def max_diff(vec1, vec2, tol):
    mask = np.logical_and(np.logical_not(np.isnan(vec1)), np.logical_not(np.isnan(vec2)))
    vec1 = vec1[mask]
    vec2 = vec2[mask]
    err = np.max(np.abs((vec1 - vec2)))
    print "Max Diff: ", err, "(> ", tol, ")"
    return err < tol

def binary_hit_or_miss(input, structure1 = None, structure2 = None,
                       output = None, origin1 = 0, origin2 = None):
    """Multi-dimensional binary hit-or-miss transform.

    An output array can optionally be provided. The origin parameters
    controls the placement of the structuring elements. If the first
    structuring element is not given one is generated with a squared
    connectivity equal to one. If the second structuring element is
    not provided, it set equal to the inverse of the first structuring
    element. If the origin for the second structure is equal to None
    it is set equal to the origin of the first.
    """
    input = numpy.asarray(input)
    if structure1 is None:
        structure1 = generate_binary_structure(input.ndim, 1)
    if structure2 is None:
        structure2 = numpy.logical_not(structure1)
    origin1 = _ni_support._normalize_sequence(origin1, input.ndim)
    if origin2 is None:
        origin2 = origin1
    else:
        origin2 = _ni_support._normalize_sequence(origin2, input.ndim)

    tmp1 = _binary_erosion(input, structure1, 1, None, None, 0, origin1,
                           0, False)
    inplace = isinstance(output, numpy.ndarray)
    result = _binary_erosion(input, structure2, 1, None, output, 0,
                             origin2, 1, False)
    if inplace:
        numpy.logical_not(output, output)
        numpy.logical_and(tmp1, output, output)
    else:
        numpy.logical_not(result, result)
        return numpy.logical_and(tmp1, result)

def compare_fixlens(samp_fixlen, fixlendist, eps=.000000001):
    nonan_samp_fixlen = samp_fixlen[np.logical_not(np.isnan(samp_fixlen))]
    nonan_fixlendist = fixlendist[np.logical_not(np.isnan(fixlendist))]
    print nonan_samp_fixlen, nonan_fixlendist
    ks, p = sts.ks_2samp(nonan_samp_fixlen, nonan_fixlendist)
    print ks, p
    return np.log(p + eps)

    def _eucl_max(self, nii1, nii2):
        origdata1 = nii1.get_data()
        origdata1 = np.logical_not(
            np.logical_or(origdata1 == 0, np.isnan(origdata1)))
        origdata2 = nii2.get_data()
        origdata2 = np.logical_not(
            np.logical_or(origdata2 == 0, np.isnan(origdata2)))

        if isdefined(self.inputs.mask_volume):
            maskdata = nb.load(self.inputs.mask_volume).get_data()
            maskdata = np.logical_not(
                np.logical_or(maskdata == 0, np.isnan(maskdata)))
            origdata1 = np.logical_and(maskdata, origdata1)
            origdata2 = np.logical_and(maskdata, origdata2)

        if origdata1.max() == 0 or origdata2.max() == 0:
            return np.NaN

        border1 = self._find_border(origdata1)
        border2 = self._find_border(origdata2)

        set1_coordinates = self._get_coordinates(border1, nii1.affine)
        set2_coordinates = self._get_coordinates(border2, nii2.affine)
        distances = cdist(set1_coordinates.T, set2_coordinates.T)
        mins = np.concatenate(
            (np.amin(distances, axis=0), np.amin(distances, axis=1)))

        return np.max(mins)

 def updateSingle(self,event=None):
     """Update the 2D data for the region of interest and intensity"""
     (vMin,vMax) = self.opPanel.getIntensityRange()
     (xMin,xMax,yMin,yMax) = self.opPanel.getRoi()
     data = self.flatdata[:,:]
   
     #Mask to zero during the summing parts
     data[np.logical_not(self.mask)] = 0
     self.posPanel.data = data
     self.posPanel.setRange(xMin,yMin,xMax,yMax)
     x=np.arange(128,0,-1)
     y=np.sum(data[:,xMin:xMax],axis=1)
     self.yPanel.SetPlot(x,y)
     #handle the x-plot
     x=np.arange(0,16,1)
     y=np.sum(data[yMin:yMax,:],axis=0)
     self.xPanel.SetPlot(x,y)
     if vMin is None:
         vMin = np.min(data)
     if vMax is None:
         vMax = np.max(data)
     self.colorbar.setRange(vMin,vMax)
     self.colorbar.update()
     #mask to vmin for the plotting
     data[np.logical_not(self.mask)] = vMin
     self.imPanel.update(self.flatdata,vMin,vMax)

 def build_tree_vector(points_r,points_c,levels_left,local_out_array):
     
     tile_rs = tile[points_r,points_c].reshape( -1,fs);
     local_out_array[0,:] = ma.mean(tile_rs,axis=0)
     
         #plt.plot(points_r,points_c,'o')
     if levels_left > 1:
         remaining_out_array = local_out_array[1:,:]
         mean_r = np.mean(points_r);
         mean_c = np.mean(points_c)
         
         offset_size = remaining_out_array.shape[0]/4
 
         top = points_r < mean_r
         bottom = np.logical_not(top)
         left = points_c < mean_c
         right = np.logical_not(left)
         
         quadrents = [ (top,right),(top,left),(bottom,left),(bottom,right)  ]
         
         #Fill the solution for all 4 quadrents 
         for idx,quadrent in enumerate(quadrents):
             q = np.logical_and(quadrent[0],quadrent[1])
             q_out = remaining_out_array[ idx*offset_size : (idx+1)*offset_size, : ]
             build_tree_vector(points_r[q],points_c[q],levels_left - 1,q_out)
         #renormilize 
         remaining_out_array *= .25

    def data(self, t=None, extrapolate=np.nan, return_indices=False):
        if t is None:
            d = self.D
            ix = np.arange(len(d))
        else:
            t = np.array(t)
            t0 = self.starttime()
            t1 = self.endtime()
            ix = np.array(np.round((t-t0)/self.dT)).astype(int)
            in_range = np.logical_and(t>=t0, t<= t1)

            if extrapolate is None:
                ix = ix[in_range]
            elif extrapolate is False:
                ix[t<t0] = 0
                ix[t>t1] = self.nD
            else:
                if any(np.logical_not(in_range)):
                    ix[np.logical_not(in_range)] = extrapolate

            d = selectalonglastdimension(self.D,ix)
        if return_indices:
            return (d,ix)
        else:
            return d

    def resample(self):
        """
        :return:
            Return the data with majority samples that form a Tomek link
            removed.
        """

        from sklearn.neighbors import NearestNeighbors

        # Find the nearest neighbour of every point
        nn = NearestNeighbors(n_neighbors=2)
        nn.fit(self.x)
        nns = nn.kneighbors(self.x, return_distance=False)[:, 1]

        # Send the information to is_tomek function to get boolean vector back
        if self.verbose:
            print("Looking for majority Tomek links...")
        links = self.is_tomek(self.y, nns, self.minc, self.verbose)

        if self.verbose:
            print("Under-sampling "
                  "performed: " + str(Counter(self.y[logical_not(links)])))

        # Return data set without majority Tomek links.
        return self.x[logical_not(links)], self.y[logical_not(links)]

def fix_nonfinite(data):
    bad_indexes = np.logical_not(np.isfinite(data))
    good_indexes = np.logical_not(bad_indexes)
    good_data = data[good_indexes]
    interpolated = np.interp(bad_indexes.nonzero()[0], good_indexes.nonzero()[0], good_data)
    data[bad_indexes] = interpolated
    return data

 def run(self, outputs_requested, **kwargs):
     # TODO find some interface that doesn't involve string parsing
     # modeled after pandas.Dataframe.query:
     #     http://pandas.pydata.org/pandas-docs/dev/generated/pandas.DataFrame.query.html
     # which implements its own computation engine:
     #     http://pandas.pydata.org/pandas-docs/dev/generated/pandas.eval.html
     # supports numpy arithmetic comparison operators:
     #     http://docs.scipy.org/doc/numpy/reference/arrays.ndarray.html#arithmetic-and-comparison-operations
     in_table = kwargs['input'].to_np()
     col_names = in_table.dtype.names
     query = self.__get_ast(col_names)
     mask = eval(compile(query, '<string>', 'eval'))
     ret = {}
     if 'output' in outputs_requested:
         uo_out = UObject(UObjectPhase.Write)
         uo_out.from_np(in_table[mask])
         ret['output'] = uo_out
     if 'complement' in outputs_requested:
         uo_comp = UObject(UObjectPhase.Write)
         uo_comp.from_np(in_table[np.logical_not(mask)])
         ret['complement'] = uo_comp
     if 'output_inds' in outputs_requested:
         uo_out_inds = UObject(UObjectPhase.Write)
         uo_out_inds.from_np(np.where(mask)[0])
         ret['output_inds'] = uo_out_inds
     if 'complement_inds' in outputs_requested:
         uo_comp_inds = UObject(UObjectPhase.Write)
         uo_comp_inds.from_np(np.where(np.logical_not(mask))[0])
         ret['complement_inds'] = uo_comp_inds
     return ret

def rms_diff(vec1, vec2, tol):
    mask = np.logical_and(np.logical_not(np.isnan(vec1)), np.logical_not(np.isnan(vec2)))
    vec1 = vec1[mask]
    vec2 = vec2[mask]
    err = np.sqrt(np.mean(((vec1 - vec2) / vec2) ** 2))
    print "RMS Error: ", err, "(> ", tol, ")"
    return err < tol

def find_large_empty_regions(grayscale_vol, min_background_voxel_count=100):
    """
    Returns mask that excludes large background (0-valued) regions, if any exist.
    """    
    if not (grayscale_vol == 0).any():
        # No background pixels.
        # We could return all ones, but we are also allowed
        # by convention to return 'None', which is faster.
        return None

    # Produce a mask that excludes 'background' pixels
    # (typically zeros around the volume edges)
    background_mask = numpy.zeros(grayscale_vol.shape, dtype=numpy.uint8)
    background_mask[grayscale_vol == 0] = 1

    # Compute connected components (cc) and toss out the small components
    cc = scipy.ndimage.label(background_mask)[0]
    cc_sizes = numpy.bincount(cc.ravel())
    small_cc_selections = cc_sizes < min_background_voxel_count
    small_cc_locations = small_cc_selections[cc]
    background_mask[small_cc_locations] = 0

    if not background_mask.any():
        # No background pixels.
        # We could return all ones, but we are also allowed
        # by convention to return 'None', which is faster.
        return None
    
    # Now background_mask == 1 for background and 0 elsewhere, so invert.
    numpy.logical_not(background_mask, out=background_mask)
    return background_mask.view(numpy.bool_)

def is_leap_year(year, gregorian=True):
    """Return True if this is a leap year in the Julian or Gregorian calendars

    Arguments:
      - `year` : (int) year

    Keywords:
      - `gregorian` : (bool, default=True) If True, use Gregorian calendar,
        else use Julian calendar

    Returns:
      - (bool) True is this is a leap year, else False.

    """
    year = np.atleast_1d(year).astype(np.int64)
    x = np.fmod(year, 4)
    if gregorian:
        x = np.fmod(year, 4)
        y = np.fmod(year, 100)
        z = np.fmod(year, 400)
        return _scalar_if_one(
            np.logical_and(np.logical_not(x),
                           np.logical_or(y, np.logical_not(z))))
    else:
        return _scalar_if_one(x == 0)

def fitAndPredict(X, y, train_mask):
    #  partition into train/test set
    Xtrain = X[train_mask]
    ytrain = y[train_mask]
    Xtest = X[np.logical_not(train_mask)]
    ytest = y[np.logical_not(train_mask)]
    
    # Fit model
    linreg = lm.LinearRegression(fit_intercept = False)
    linreg.fit(Xtrain, ytrain)
    
    # Extract parameters
    coef = linreg.coef_
    names = X.columns.tolist()
    print([name + ':' + str(round(w,3)) for name, w in zip(names, coef)])
    
    # Measure train error
    yhat = linreg.predict(Xtrain)
    (mse, stderr) = L2loss(yhat, ytrain) 
    print 'Train mse {0:0.3f}, stderr {1:0.3f}'.format(mse, stderr) 
    
    # Measure test error
    yhat = linreg.predict(Xtest)
    (mse, stderr) = L2loss(yhat, ytest) 
    print 'Test mse {0:0.3f}, stderr {1:0.3f}'.format(mse, stderr) 

 def test_insert_nan(self):
     """Test fetching of null values"""
     b = self.rel.fetch('value', order_by='id')
     assert_true((np.isnan(self.a) == np.isnan(b)).all(),
                 'incorrect handling of Nans')
     assert_true(np.allclose(self.a[np.logical_not(np.isnan(self.a))], b[np.logical_not(np.isnan(b))]),
                 'incorrect storage of floats')