def assert_almost_equal_inf(x, y, decimal=6, msg=None):
    x = np.atleast_1d(x)
    y = np.atleast_1d(y)
    assert_equal(np.isposinf(x), np.isposinf(y))
    assert_equal(np.isneginf(x), np.isneginf(y))
    assert_equal(np.isnan(x), np.isnan(y))
    assert_almost_equal(x[np.isfinite(x)], y[np.isfinite(y)])

def non_matches(arr, val):
    '''
    Given a ndarray and an arbitrary 
    value, including np.nan, np.inf, etc.,
    return an ndarray that contains 
    only elements that are *not* equal 
    to val.  
    
    :param arr: n-dimensional numpy array
    :type arr: numpy.ndarray
    :param val: value, including special values numpy.nan, numpy.inf, numpy.neginf, etc.
    :type val: ANY.
    '''
    
    # Special value?
    if np.isfinite(val):
        # No, just normal value:
        return arr[arr != val]
    # Is special value, such as numpy.nan.
    # Create ndarray with True/False entries
    # that reflect which entries are not equal
    # to val:
    elif np.isnan(val):
        cond = np.logical_not(np.isnan(arr))
    elif np.isinf(val):
        cond = np.logical_not(np.isinf(arr))
    elif np.isneginf(val):
        cond = np.logical_not(np.isneginf(arr))
    elif np.isposinf(val):
        cond = np.logical_not(np.isposinf(arr))
        
    # Use the True/False ndarray as a mask
    # over arr:
    return arr[cond]

def compute_frequencydomaincoef(x, data_length_sec, sampling_frequency, nfreq_bands, win_length_sec, stride_sec):
    n_channels = x.shape[1]
    n_timesteps = (data_length_sec - win_length_sec) / stride_sec + 1
    n_fbins = nfreq_bands
    xfreq = np.zeros((n_timesteps, 136))
    x2 = np.zeros((n_channels, n_fbins, n_timesteps))
    for i in range(n_channels):
        xc = np.zeros((n_fbins, n_timesteps))
        for frame_num, w in enumerate(range(0, data_length_sec - win_length_sec + 1, stride_sec)):
            xw = x[w * sampling_frequency: (w + win_length_sec) * sampling_frequency,i]#window
            fft = np.log10(np.absolute(np.fft.rfft(xw)))
            fft_freq = np.fft.rfftfreq(n=xw.shape[-1], d=1.0 / sampling_frequency) ## return the FFT sample  frequency
            xc[:nfreq_bands, frame_num] = group_into_bands(fft, fft_freq, nfreq_bands)
        x2[i, :, :] = xc
    for j in range(n_timesteps):

        x2[:, :, j][np.isneginf(x2[:, :, j])] = 0
        scaled = preprocessing.scale(x2[:, :, j], axis=0)

        matrix = CorrelationMatrix().apply(scaled)
        matrix[np.isneginf(matrix)] = 0
        matrix[np.isnan(matrix)] = 0

        eigenvalue = Eigenvalues().apply(matrix)

        freqdomaincor = upper_right_triangle(matrix)
        xfreq[j, :] = np.concatenate((freqdomaincor, eigenvalue))
    xfreq[np.isneginf(xfreq)] = 0
    xfreq[np.isnan(xfreq)] = 0
    return xfreq

def min_sum_diffs(filename, args):
    """Sum of the differences (in dB) between responses and a reference response.
        
    Args:
        filename (str): Name of output file
        args (dict): 'refresp' key with path & filename of reference response; 'outputs' key with a list of names (IDs) of outputs (rxs) from input file
        
    Returns:
        diffdB (float): Sum of the differences (in dB) between responses and a reference response
    """

    # Load (from gprMax output file) the reference response
    f = h5py.File(args['refresp'], 'r')
    tmp = f['/rxs/rx1/']
    fieldname = list(tmp.keys())[0]
    refresp = np.array(tmp[fieldname])

    # Load (from gprMax output file) the response
    f = h5py.File(filename, 'r')
    nrx = f.attrs['nrx']
    
    diffdB = 0
    outputs = 0
    for rx in range(1, nrx + 1):
        output = f['/rxs/rx' + str(rx) + '/']
        if output.attrs['Name'] in args['outputs']:
            outputname = list(output.keys())[0]
            modelresp = np.array(output[outputname])
            # Calculate sum of differences
            tmp = 20 * np.log10(np.abs(modelresp - refresp) / np.amax(np.abs(refresp)))
            tmp = np.abs(np.sum(tmp[-np.isneginf(tmp)])) / len(tmp[-np.isneginf(tmp)])
            diffdB += tmp
            outputs += 1

    return diffdB / outputs

def imagesDiffer(imageArr1, imageArr2, skipMaskArr=None, rtol=1.0e-05, atol=1e-08):
    """Compare the pixels of two image arrays; return True if close, False otherwise
    
    Inputs:
    - image1: first image to compare
    - image2: second image to compare
    - skipMaskArr: pixels to ignore; nonzero values are skipped
    - rtol: relative tolerance (see below)
    - atol: absolute tolerance (see below)
    
    rtol and atol are positive, typically very small numbers.
    The relative difference (rtol * abs(b)) and the absolute difference "atol" are added together
    to compare against the absolute difference between "a" and "b".
    
    Return a string describing the error if the images differ significantly, an empty string otherwise
    """
    retStrs = []
    if skipMaskArr != None:
        maskedArr1 = numpy.ma.array(imageArr1, copy=False, mask = skipMaskArr)
        maskedArr2 = numpy.ma.array(imageArr2, copy=False, mask = skipMaskArr)
        filledArr1 = maskedArr1.filled(0.0)
        filledArr2 = maskedArr2.filled(0.0)
    else:
        filledArr1 = imageArr1
        filledArr2 = imageArr2

    nan1 = numpy.isnan(filledArr1)
    nan2 = numpy.isnan(filledArr2)
    if numpy.any(nan1 != nan2):
        retStrs.append("NaNs differ")

    posinf1 = numpy.isposinf(filledArr1)
    posinf2 = numpy.isposinf(filledArr2)
    if numpy.any(posinf1 != posinf2):
        retStrs.append("+infs differ")

    neginf1 = numpy.isneginf(filledArr1)
    neginf2 = numpy.isneginf(filledArr2)
    if numpy.any(neginf1 != neginf2):
        retStrs.append("-infs differ")

    # compare values that should be comparable (are neither infinite, nan nor masked)
    valSkipMaskArr = nan1 | nan2 | posinf1 | posinf2 | neginf1 | neginf2
    if skipMaskArr != None:
        valSkipMaskArr |= skipMaskArr
    valMaskedArr1 = numpy.ma.array(imageArr1, copy=False, mask = valSkipMaskArr)
    valMaskedArr2 = numpy.ma.array(imageArr2, copy=False, mask = valSkipMaskArr)
    valFilledArr1 = valMaskedArr1.filled(0.0)
    valFilledArr2 = valMaskedArr2.filled(0.0)
    
    if not numpy.allclose(valFilledArr1, valFilledArr2, rtol=rtol, atol=atol):
        errArr = numpy.abs(valFilledArr1 - valFilledArr2)
        maxErr = errArr.max()
        maxPosInd = numpy.where(errArr==maxErr)
        maxPosTuple = (maxPosInd[1][0], maxPosInd[0][0])
        errStr = "maxDiff=%s at position %s; value=%s vs. %s" % \
            (maxErr, maxPosTuple, valFilledArr1[maxPosInd][0], valFilledArr2[maxPosInd][0])
        retStrs.insert(0, errStr)
    return "; ".join(retStrs)

def figS4(data_dir=mydir, figname = 'FigS4', saveAs = 'eps'):
    models = ['lognorm', 'mete', 'zipf']
    fig = plt.figure()
    count = 0
    gs = gridspec.GridSpec(4, 4)
    #gs.update(wspace=0.1, hspace=0.1)
    for i in range(0, 4, 2):
        for j in range(0, 4, 2):
            if count < 2:
                ax = plt.subplot(gs[i:i+2, j:j+2], adjustable='box-forced')
                count += 1
            else:
                ax = plt.subplot(gs[i:i+2, 1:3], adjustable='box-forced')
            if i == 0 and j == 0:
                NSR2 = importData.import_NSR2_data(data_dir + \
                'data/NSR2/Stratified/lognorm_pln_NSR2_stratify.txt')
                ax.set_title("Lognormal", fontsize = 18)
                ll = np.asarray(list(((NSR2["ll"]))))
                ll = ll[np.isneginf(ll) == False]
                print 'Lognorm: mean = ' + str(np.mean(ll)) + ' std = ' + str(np.std(ll))
            elif i == 0 and j == 2:
                NSR2 = importData.import_NSR2_data(data_dir + \
                'data/NSR2/Stratified/zipf_mle_NSR2_stratify.txt')
                ax.set_title("Zipf", fontsize = 18)
                ll = np.asarray(list(((NSR2["ll"]))))
                ll = ll[np.isneginf(ll) == False]
                print 'Zipf: mean = ' + str(np.mean(ll)) + ' std = ' + str(np.std(ll))
            elif i == 2 and j == 0:
                NSR2 = importData.import_NSR2_data(data_dir + \
                'data/NSR2/Stratified/mete_NSR2_stratify.txt')
                ax.set_title("Log-series", fontsize = 18)
                ll = np.asarray(list(((NSR2["ll"]))))
                ll = ll[np.isneginf(ll) == False]
                print 'Log-series: mean = ' + str(np.mean(ll)) + ' std = ' + str(np.std(ll))
            else:
                continue

            ax.set( adjustable='box-forced')
            KDE = mo.CV_KDE(ll)
            #ax.hist(ll, 30, fc='gray', histtype='stepfilled', alpha=0.5, normed=True)
            ax.plot(KDE[0], KDE[1], linewidth=3, alpha=0.8 , color = 'blue')
            ax.yaxis.set_major_formatter(mticker.FormatStrFormatter('%.0E'))
            ax.xaxis.set_major_formatter(mticker.FormatStrFormatter('%.0E'))

            ax.set_xlim([min(KDE[0]), 0])
            plt.xticks(fontsize = 7)
            plt.yticks(fontsize = 7)
            ax.set_xlabel('Log-likelihood', fontsize = 16)
            ax.set_ylabel('Probability density', fontsize = 14)
            plt.setp(ax.get_xticklabels()[::2], visible=False)
            plt.setp(ax.get_yticklabels()[::2], visible=False)

    fig_name = str(mydir + 'figures/' + figname + '_RGB.' + saveAs)
    fig.subplots_adjust(left=0.1, bottom = 0.1,hspace=0.1)
    fig.tight_layout()#pad=1.2, w_pad=0.8, h_pad=0.8
    #fig.text(0.50, 0.017, 'Log-likelihood', ha='center', va='center', fontsize=15)
    #fig.text(0.04, 0.5, 'Probability', ha='center', va='center', rotation='vertical', fontsize=20)
    plt.savefig(fig_name, dpi=600, format = saveAs)
    plt.close()

 def test_infinity_neg(self):
     x = -numpy.inf
     y = self.sendAndReceive(x)
     self.assertEqual(y, x)
     self.assert_(numpy.isneginf(x))
     self.assert_(numpy.isneginf(y))
     self.assertEqual(numpy.array(x).shape, y.shape)
     self.assertEqual(numpy.array(x).dtype, y.dtype)

def _generate_colorbar_ticks_label(
    data_transform=False, colorbarlabel=None, trans_base_list=None, forcelabel=None, plotlev=None, plotlab=None
):
    """
    Return (colorbar_ticks,colorbar_labels)
    """
    # data_transform==True and levels!=None
    if data_transform == True:
        if colorbarlabel != None:
            colorbarlabel = pb.iteflat(colorbarlabel)
            transformed_colorbarlabel_ticks, x, y, trans_base_list = mathex.plot_array_transg(
                colorbarlabel, trans_base_list, copy=True
            )

        # Note if/else blocks are organized in 1st tire by check if the two
        # ends are -inf/inf and 2nd tire by check if colorbarlabel is None
        if np.isneginf(plotlab[0]) and np.isposinf(plotlab[-1]):
            if colorbarlabel != None:
                ftuple = (transformed_colorbarlabel_ticks, colorbarlabel)
            else:
                ftuple = (plotlev, plotlab[1:-1])
        elif np.isneginf(plotlab[0]) or np.isposinf(plotlab[-1]):
            raise ValueError("It's strange to set only side as infitive")
        else:
            if colorbarlabel != None:
                ftuple = (transformed_colorbarlabel_ticks, colorbarlabel)
            else:
                ftuple = (plotlev, plotlab)

    # data_transform==False
    else:
        if np.isneginf(plotlab[0]) and np.isposinf(plotlab[-1]):
            # if colorbarlabel is forced, then ticks and ticklabels will be forced.
            if colorbarlabel != None:
                ftuple = (colorbarlabel, colorbarlabel)
            # This by default will be done, it's maintained here only for clarity.
            else:
                ftuple = (plotlab[1:-1], plotlab[1:-1])
        elif np.isneginf(plotlab[0]) or np.isposinf(plotlab[-1]):
            raise ValueError("It's strange to set only side as infitive")
        else:
            if colorbarlabel != None:
                ftuple = (colorbarlabel, colorbarlabel)
            else:
                ftuple = (plotlab, plotlab)

    ftuple = list(ftuple)
    if forcelabel != None:
        if len(forcelabel) != len(ftuple[1]):
            raise ValueError(
                """the length of the forcelabel and the
                length of labeled ticks is not equal!"""
            )
        else:
            ftuple[1] = forcelabel

    return ftuple

def _transform_data(pdata, levels, data_transform):
    """
    Return [pdata,plotlev,plotlab,extend,trans_base_list];
    if data_transform == False, trans_base_list = None.

    Notes:
    ------
    pdata: data used for contourf plotting.
    plotlev: the levels used in contourf plotting.
    extend: the value for parameter extand in contourf.
    trans_base_list: cf. mathex.plot_array_transg
    """
    if levels == None:
        ftuple = (pdata, None, None, "neither")
        if data_transform == True:
            raise Warning("Strange levels is None but data_transform is True")
    else:
        if data_transform == True:
            # make the data transform before plotting.
            pdata_trans, plotlev, plotlab, trans_base_list = mathex.plot_array_transg(pdata, levels, copy=True)
            if np.isneginf(plotlab[0]) and np.isposinf(plotlab[-1]):
                ftuple = (pdata_trans, plotlev[1:-1], plotlab, "both")
            elif np.isneginf(plotlab[0]) or np.isposinf(plotlab[-1]):
                raise ValueError(
                    """only one extreme set as infinitive, please
                    set both as infinitive if arrow colorbar is wanted."""
                )
            else:
                ftuple = (pdata_trans, plotlev, plotlab, "neither")
        # data_transform==False
        else:
            plotlev = pb.iteflat(levels)
            plotlab = pb.iteflat(levels)
            if np.isneginf(plotlab[0]) and np.isposinf(plotlab[-1]):
                # here the levels would be like [np.NINF,1,2,3,np.PINF]
                # in following contourf, all values <1 and all values>3 will be
                # automatically plotted in the color of two arrows.
                # easy to see in this example:
                # a=np.tile(np.arange(10),10).reshape(10,10);
                # fig,ax=g.Create_1Axes();
                # cs=ax.contourf(a,levels=np.arange(2,7),extend='both');
                # plt.colorbar(cs)
                ftuple = (pdata, plotlev[1:-1], plotlab, "both")
            elif np.isneginf(plotlab[0]) or np.isposinf(plotlab[-1]):
                raise ValueError(
                    """only one extreme set as infinitive, please
                    set both as infinitive if arrow colorbar is wanted."""
                )
            else:
                ftuple = (pdata, plotlev, plotlab, "neither")
    datalist = list(ftuple)

    if data_transform == True:
        datalist.append(trans_base_list)
    else:
        datalist.append(None)
    return datalist

def _diagnose(self):

    # Update log.
    self.logger.debug("diagnose: data: shape: " + str(self.data.shape))
    self.logger.debug("diagnose: data: dtype: " + str(self.data.dtype))
    self.logger.debug("diagnose: data: size: %.2fMB", self.data.nbytes * 9.53674e-7)
    self.logger.debug("diagnose: data: nans: " + str(np.sum(np.isnan(self.data))))
    self.logger.debug("diagnose: data: -inf: " + str(np.sum(np.isneginf(self.data))))
    self.logger.debug("diagnose: data: +inf: " + str(np.sum(np.isposinf(self.data))))
    self.logger.debug("diagnose: data: positives: " + str(np.sum(self.data > 0)))
    self.logger.debug("diagnose: data: negatives: " + str(np.sum(self.data < 0)))
    self.logger.debug("diagnose: data: mean: " + str(np.mean(self.data)))
    self.logger.debug("diagnose: data: min: " + str(np.min(self.data)))
    self.logger.debug("diagnose: data: max: " + str(np.max(self.data)))

    self.logger.debug("diagnose: data_white: shape: " + str(self.data_white.shape))
    self.logger.debug("diagnose: data_white: dtype: " + str(self.data_white.dtype))
    self.logger.debug("diagnose: data_white: size: %.2fMB", self.data_white.nbytes * 9.53674e-7)
    self.logger.debug("diagnose: data_white: nans: " + str(np.sum(np.isnan(self.data_white))))
    self.logger.debug("diagnose: data_white: -inf: " + str(np.sum(np.isneginf(self.data_white))))
    self.logger.debug("diagnose: data_white: +inf: " + str(np.sum(np.isposinf(self.data_white))))
    self.logger.debug("diagnose: data_white: positives: " + str(np.sum(self.data_white > 0)))
    self.logger.debug("diagnose: data_white: negatives: " + str(np.sum(self.data_white < 0)))
    self.logger.debug("diagnose: data_white: mean: " + str(np.mean(self.data_white)))
    self.logger.debug("diagnose: data_white: min: " + str(np.min(self.data_white)))
    self.logger.debug("diagnose: data_white: max: " + str(np.max(self.data_white)))

    self.logger.debug("diagnose: data_dark: shape: " + str(self.data_dark.shape))
    self.logger.debug("diagnose: data_dark: dtype: " + str(self.data_dark.dtype))
    self.logger.debug("diagnose: data_dark: size: %.2fMB", self.data_dark.nbytes * 9.53674e-7)
    self.logger.debug("diagnose: data_dark: nans: " + str(np.sum(np.isnan(self.data_dark))))
    self.logger.debug("diagnose: data_dark: -inf: " + str(np.sum(np.isneginf(self.data_dark))))
    self.logger.debug("diagnose: data_dark: +inf: " + str(np.sum(np.isposinf(self.data_dark))))
    self.logger.debug("diagnose: data_dark: positives: " + str(np.sum(self.data_dark > 0)))
    self.logger.debug("diagnose: data_dark: negatives: " + str(np.sum(self.data_dark < 0)))
    self.logger.debug("diagnose: data_dark: mean: " + str(np.mean(self.data_dark)))
    self.logger.debug("diagnose: data_dark: min: " + str(np.min(self.data_dark)))
    self.logger.debug("diagnose: data_dark: max: " + str(np.max(self.data_dark)))

    self.logger.debug("diagnose: theta: shape: " + str(self.theta.shape))
    self.logger.debug("diagnose: theta: dtype: " + str(self.theta.dtype))
    self.logger.debug("diagnose: theta: size: %.2fMB", self.theta.nbytes * 9.53674e-7)
    self.logger.debug("diagnose: theta: nans: " + str(np.sum(np.isnan(self.theta))))
    self.logger.debug("diagnose: theta: -inf: " + str(np.sum(np.isneginf(self.theta))))
    self.logger.debug("diagnose: theta: +inf: " + str(np.sum(np.isposinf(self.theta))))
    self.logger.debug("diagnose: theta: positives: " + str(np.sum(self.theta > 0)))
    self.logger.debug("diagnose: theta: negatives: " + str(np.sum(self.theta < 0)))
    self.logger.debug("diagnose: theta: mean: " + str(np.mean(self.theta)))
    self.logger.debug("diagnose: theta: min: " + str(np.min(self.theta)))
    self.logger.debug("diagnose: theta: max: " + str(np.max(self.theta)))

    self.logger.info("diagnose [ok]")

    def test_neginf(self):
        arr =np.empty(100)
        arr[:] = -np.inf
        for np_func, acml_func in self.vector_funcs:
            np_out = np_func(arr)
            acml_out = acml_func(arr)

            equal_nan = np.isnan(np_out) == np.isnan(acml_out)
            equal_posinf = np.isposinf(np_out) == np.isposinf(acml_out)
            equal_neginf = np.isneginf(np_out) == np.isneginf(acml_out)
            self.assertTrue( np.alltrue(equal_nan), msg="NaN-test failed for %s" % acml_func)
            self.assertTrue( np.alltrue(equal_posinf), msg="posinf-test failed for %s" % acml_func)
            self.assertTrue( np.alltrue(equal_neginf), msg="neginf-test failed for %s" % acml_func)

 def filt_butter(data, samp_freq, butter_freq, axis=-1):
     '''
     Filter data with a 2nd order butterworth filter.
     
     Parameters
     ==========
       data: ndarray
       samp_freq: sampling period (s)
       butter_freq: [cutoff_low, cutoff_high] (Hz), can be infinite
       axis (optional): axis along which to filter, default = -1
     Returns
     =======
       filtNs: filtered version of data
     '''
     order = 2
     ny = 0.5 / samp_freq # Nyquist frequency
     cof = butter_freq / ny # normalized cutoff freq
     if np.isneginf(cof[0]) and np.isfinite(cof[1]):
         # lowpass
         cof1 = cof[1]
         b, a = scipy.signal.butter(order, cof1, btype='low')
         filtNs = scipy.signal.filtfilt(b, a, data, axis=axis)
     elif np.isfinite(cof[0]) and np.isinf(cof[1]):
         # highpass
         cof1 = cof[0]
         b, a = scipy.signal.butter(order, cof1, btype='high')
         filtNs = scipy.signal.filtfilt(b, a, data, axis=axis)
     elif np.isfinite(cof[0]) and np.isfinite(cof[1]):
         # bandpass
         b, a = scipy.signal.butter(order, cof, btype='band')
         filtNs = scipy.signal.filtfilt(b, a, data, axis=axis)
     else:
         raise Exception('filt_butter called with bad cutoff frequency')
     filtNs /= samp_freq # normalize to rate
     return filtNs

    def optimize_A(self, A):
        """Find optimal transformation matrix A by minimization.

        Parameters
        ----------
        A : ndarray
        The transformation matrix A.

        Returns
        -------
        A : ndarray
            The transformation matrix.
        """
        flat_map, square_map = get_maps(A)
        alpha = to_flat(1.0 * A, flat_map)

        obj = lambda x: -1 * self.objective_function(x, self.T, self.right_eigenvectors, square_map, self.populations)
        self.obj = obj
        self.alpha = alpha.copy()

        logger.info("Initial value of objective function: f = %f", obj(alpha))

        alpha = scipy.optimize.anneal(obj, alpha, lower=0.0, maxiter=1, schedule="boltzmann", dwell=1000, feps=1E-3, boltzmann=2.0, T0=1.0)[0]

        alpha = scipy.optimize.fmin(obj, alpha, full_output=True, xtol=1E-4, ftol=1E-4, maxfun=5000, maxiter=100000)[0]

        logger.info("Final value: f = %f" % (obj(alpha)))

        if np.isneginf(obj(alpha)):
            raise(ValueError("Error: minimization has not located a feasible point."))

        A = to_square(alpha, square_map)

        return A

def map_to_range(v, oi, oa, ni, na):
    if numpy.isinf(v):
        return na
    elif numpy.isneginf(v):
        return ni
    else:
        return(((v - oi) * (na - ni)) / (oa - oi)) + ni

def traverse_data(datum, is_numpy=is_numpy, use_numpy=True):
    """recursively dig until a flat list is found
    if numpy is available convert the flat list to a numpy array
    and send off to transform_array() to handle nan, inf, -inf
    otherwise iterate through items in array converting non-json items

    Args:
        datum (list) : a list of values or lists
        is_numpy: True if numpy is present (see imports)
        use_numpy: toggle numpy as a dependency for testing purposes
    """
    is_numpy = is_numpy and use_numpy
    if is_numpy and not any(isinstance(el, (list, tuple)) for el in datum):
        return transform_array(np.asarray(datum))
    datum_copy = []
    for item in datum:
        if isinstance(item, (list, tuple)):
            datum_copy.append(traverse_data(item))
        elif isinstance(item, float):
            if np.isnan(item):
                item = 'NaN'
            elif np.isposinf(item):
                item = 'Infinity'
            elif np.isneginf(item):
                item = '-Infinity'
            datum_copy.append(item)
        else:
            datum_copy.append(item)
    return datum_copy

 def likelihood_function ( self, theta ):
     """For example! This function ought to be overridden by the user, and maybe extended with whatever extra parameters you need to get hold of your observations, or model driver parameters.
     This function method calculates the likelihood function for a vector M{\theta}. Usually, you have a model you run with these parameters as inputs (+ some driver data), and some observations that go with the output of the forward model output. These two sets of values are combined in some sort of cost function/likelihood function. A common criterion is to assume that the model is able to perfectly replicate the observations (given a proper parametrisation). The only mismatch between model output and observations is then given by the uncertainty with which the measurement is performed, and we can encode this as a zero-mean Normal distribution. The variance of this distribution is then related to the observational error. If different measurements are used, a multivariate normal is useful, and correlation between observations can also be included, if needs be.
     """
     means = numpy.matrix([-3.0, 2.8])
     means = numpy.matrix([-5.0, 5])
     sigma1 = 1.0
     sigma2 = 2#0.5
     rho = -0.5#-0.1
     covar = numpy.matrix([[sigma1*sigma1,rho*sigma1*sigma2],[rho*sigma1*sigma2,sigma2*sigma2]])
     inv_covar = numpy.linalg.inv ( covar ) # numpy.matrix([[  5.26315789,   9.47368421],\
                     #[  9.47368421,  21.05263158]])
     det_covar = numpy.linalg.det( covar ) #0.047499999999999987
     N = means.shape[0]
     X = numpy.matrix(means- theta)
     #X = theta
     #p = full_gauss_den(X, means, covar, True)
     #This is just lazy... Using libraries to invert a 2x2 matrix & calc. its determinant....
     #Also, the log calculations could be done more efficiently and stored, but...
     p = pow(1.0/(2*numpy.pi), N/2.)
     p = p / numpy.sqrt ( numpy.linalg.det (covar))
     #p = 0.73025296137109341 # Precalc'ed
     #p = p *    numpy.exp (-0.5*X*inv_covar*X.transpose())
     a = X*inv_covar*X.T
     p = p*numpy.exp(-0.5*a)
     #pdb.set_trace()
     p = numpy.log(p)
     if numpy.isneginf(p):
         p = numpy.log(1.0E-300)
     return p

def merciless_print(i, node, fn):
    """Debugging theano. Prints inputs and outputs at every point.
    In case NaN, Inf or -Inf is detected, fires up the pdb debugger."""
    print ''
    print '-------------------------------------------------------'
    print 'Node %s' % str(i)
    theano.printing.debugprint(node)
    print 'Inputs : %s' % [input for input in fn.inputs]
    print 'Outputs: %s' % [output for output in fn.outputs]
    print 'Node:'
    for output in fn.outputs:
        try:
            if numpy.isnan(output[0]).any():
                print '*** NaN detected ***'
                theano.printing.debugprint(node)
                print 'Inputs : %s' % [input[0] for input in fn.inputs]
                print 'Outputs: %s' % [output[0] for output in fn.outputs]
                pdb.set_trace()
                raise ValueError('Found NaN in computation!')
            if numpy.isposinf(output[0]).any() or numpy.isneginf(output[0]).any():
                print '*** Inf detected ***'
                theano.printing.debugprint(node)
                print 'Inputs : %s' % [input[0] for input in fn.inputs]
                print 'Outputs: %s' % [output[0] for output in fn.outputs]
                pdb.set_trace()
                raise ValueError('Found Inf in computation!')
        except TypeError:
            logging.debug('Couldn\'t check node for NaN/Inf: {0}'.format(node))

def takeRatio(num,den):
    toReturn = num.copy()
    toReturn['data'] = np.log(num['data']/den['data'])
    whereBad = np.isnan(toReturn['data']) | np.isinf(toReturn['data']) | np.isneginf(toReturn['data'])
    toReturn['data'][whereBad] = 0.0
    
    return toReturn

    def Draw(self, args=None):
        """Draw the various functions"""

        if not args or "SAME" not in args:
            # make a 'blank' function to occupy the complete range of x values:
            lower_lim = min([lim[0] for lim in self.functions_dict.keys()])
            if np.isneginf(lower_lim):
                lower_lim = -999
            upper_lim = max([lim[1] for lim in self.functions_dict.keys()])
            if np.isposinf(upper_lim):
                upper_lim = 999
            blank = ROOT.TF1("blank" + str(np.random.randint(0, 10000)), "1.5", lower_lim, upper_lim)
            blank.Draw()
            max_value = max([func.GetMaximum(lim[0], lim[1])
                             for lim, func in self.functions_dict.iteritems()]) * 1.1
            blank.SetMaximum(max_value)
            min_value = min([func.GetMinimum(lim[0], lim[1])
                             for lim, func in self.functions_dict.iteritems()]) * 0.9
            blank.SetMinimum(min_value)
            ROOT.SetOwnership(blank, False)  # NEED THIS SO IT ACTUALLY GETS DRAWN. SERIOUSLY, WTF?!
            blank.SetLineColor(ROOT.kWhite)

        # now draw the rest of the functions
        args = "" if not args else args
        for func in self.functions_dict.values():
            func.Draw("SAME" + args)

def gradient_desc_ridge(X, Y, W, alpha, lambd, num_iter=1000, conv_tol=0.01, check_interval=500):
    c = float("inf")
    log("Learn Rate", alpha)
    for i in range(num_iter):
        #
        # delta =  2/N SIGMA[(XW - Y)*x] + 2 * \lambd * W
        diff = predict(X, W) - Y
        delta = np.sum(np.multiply(X, diff), axis=0)  # sum top to bottom for each attribute
        delta = delta * 2.0 / len(Y)
        delta = np.array([delta]).transpose()  # restore vector shape of (n_attr x 1)
        delta = delta + (2 * lambd * W)  # Vectors addition

        W = W - alpha * delta

        if i % check_interval == 0:
            predY = predict(X, W)
            newcost = MSECost(predY, Y)

            log("#%d, cost = %.8g" % (i, newcost))
            if np.isnan(newcost) or np.isinf(newcost) or np.isneginf(newcost):
                raise Exception("ERROR: number overflow, please adjust learning rate")
            diff = abs(newcost - c)
            c = newcost
            if diff < conv_tol:
                log("Converged with tolerance %f " % conv_tol)
                break
        if not quiet and i % (check_interval * 10) == 0:
            print(W.flatten())
    return W