numpy.polyfit

Here are the examples of the python api numpy.polyfit taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.

50 Examples 7

Example 1

Project: chaco Source File: regression_lasso.py
    def _selection_changed_fired(self, event):
        indices = self.selection_datasource.metadata["selection"]
        if any(indices):
            x = compress(indices, self.component.index.get_data())
            y = compress(indices, self.component.value.get_data())
            if len(x) < 2 or len(y) < 2:
                self.fit_params = None
                self.centroid = None
            else:
                self.fit_params = tuple(polyfit(x,y,1))
                self.centroid = (sum(x)/len(x)), (sum(y)/len(y))
        else:
            self.fit_params = None
            self.centroid = None
        return

Example 2

Project: wordgraph Source File: analysers.py
    def get_validity(self):
        x_values = [point.x for point in self.points]
        y_values = [point.y for point in self.points]
        par = np.polyfit(x_values, y_values, 1, full=True)

        self.gradient = par[0][0]
        self.constant = par[0][1]

        residuals = np.var([(self.gradient * point.x + self.constant - point.y)
            for point in self.points])
        # variance = np.var(y_values)
        # Rsqr = np.round(1 - residuals / variance, decimals=2)
        return 1 - residuals

Example 3

Project: spinalcordtoolbox Source File: sct_smooth_spinal_cord_shifting_centerline.py
Function: polynome_centerline
def polynome_centerline(x_centerline,y_centerline,z_centerline):
    """Fit polynomial function through centerline"""

    # Fit centerline in the Z-X plane using polynomial function
    print '\nFit centerline in the Z-X plane using polynomial function...'
    coeffsx = np.polyfit(z_centerline, x_centerline, deg=5)
    polyx = np.poly1d(coeffsx)
    x_centerline_fit = np.polyval(polyx, z_centerline)

    #Fit centerline in the Z-Y plane using polynomial function
    print '\nFit centerline in the Z-Y plane using polynomial function...'
    coeffsy = np.polyfit(z_centerline, y_centerline, deg=5)
    polyy = np.poly1d(coeffsy)
    y_centerline_fit = np.polyval(polyy, z_centerline)


    return x_centerline_fit,y_centerline_fit

Example 4

Project: spinalcordtoolbox Source File: straighten.py
Function: polynome_centerline
def polynome_centerline(x_centerline,y_centerline,z_centerline):
    """Fit polynomial function through centerline"""
    
    # Fit centerline in the Z-X plane using polynomial function
    print '\nFit centerline in the Z-X plane using polynomial function...'
    coeffsx = numpy.polyfit(z_centerline, x_centerline, deg=param.deg_poly)
    polyx = numpy.poly1d(coeffsx)
    x_centerline_fit = numpy.polyval(polyx, z_centerline)
    
    #Fit centerline in the Z-Y plane using polynomial function
    print '\nFit centerline in the Z-Y plane using polynomial function...'
    coeffsy = numpy.polyfit(z_centerline, y_centerline, deg=param.deg_poly)
    polyy = numpy.poly1d(coeffsy)
    y_centerline_fit = numpy.polyval(polyy, z_centerline)
    
    return x_centerline_fit,y_centerline_fit,polyx,polyy

Example 5

Project: stanza Source File: trigger.py
Function: slope
    def slope(self):
        """
        :return: the esitmated slope for points in the current window
        """
        x = range(self.window_size)
        y = self.vals
        slope, bias = np.polyfit(x, y, 1)
        return slope

Example 6

Project: waveform-analyzer Source File: common.py
def parabolic_polyfit(f, x, n):
    """
    Use the built-in polyfit() function to find the peak of a parabola

    f is a vector and x is an index for that vector.

    n is the number of samples of the curve used to fit the parabola.
    """
    a, b, c = polyfit(arange(x-n//2, x+n//2+1), f[x-n//2:x+n//2+1], 2)
    xv = -0.5 * b/a
    yv = a * xv**2 + b * xv + c
    return (xv, yv)

Example 7

Project: spinalcordtoolbox Source File: msct_smooth.py
def polynomial_fit(x,y,degree):

    import numpy as np

    coeffs = np.polyfit(x, y, degree)
    poly = np.poly1d(coeffs)
    y_fit = np.polyval(poly, x)
 
    return y_fit, poly

Example 8

Project: mtpy Source File: pek1dplotting.py
Function: get_profile
    def get_profile(self):
        """
        get location of profile by linear regression
        """
        if self.Model_suite.x is None:
            print "can't get profile, no x y locations"
            return
        
        x = self.Model_suite.x
        y = self.Model_suite.y
        
        self.profile = np.polyfit(x,y,1)

Example 9

Project: EyeTab Source File: ransac_eyelids.py
Function: fit_parabola
def fit_parabola(xs, ys):
    
    a, b, c = np.polyfit(xs, ys, 2)

    if a < 0: raise BadFitShape()
    return a, b, c

Example 10

Project: seaborn Source File: test_linearmodels.py
Function: test_residplot
    def test_residplot(self):

        x, y = self.df.x, self.df.y
        ax = lm.residplot(x, y)

        resid = y - np.polyval(np.polyfit(x, y, 1), x)
        x_plot, y_plot = ax.collections[0].get_offsets().T

        npt.assert_array_equal(x, x_plot)
        npt.assert_array_almost_equal(resid, y_plot)

Example 11

Project: quant_at Source File: hurst.py
def hurst(ts):
	"""Returns the Hurst Exponent of the time series vector ts"""
	# Create the range of lag values
	lags = range(2, 100)

	# Calculate the array of the variances of the lagged differences
	tau = [sqrt(std(subtract(ts[lag:], ts[:-lag]))) for lag in lags]

	# Use a linear fit to estimate the Hurst Exponent
	poly = polyfit(log(lags), log(tau), 1)

	# Return the Hurst exponent from the polyfit output
	return poly[0]*2.0

Example 12

Project: spinalcordtoolbox Source File: sct_flatten_sagittal.py
Function: polynome_centerline
def polynome_centerline(x_centerline,y_centerline,z_centerline):
    """Fit polynomial function through centerline"""
    
    # Fit centerline in the Z-X plane using polynomial function
    print '\nFit centerline in the Z-X plane using polynomial function...'
    coeffsx = numpy.polyfit(z_centerline, x_centerline, deg=5)
    polyx = numpy.poly1d(coeffsx)
    x_centerline_fit = numpy.polyval(polyx, z_centerline)
    
    # Fit centerline in the Z-Y plane using polynomial function
    print '\nFit centerline in the Z-Y plane using polynomial function...'
    coeffsy = numpy.polyfit(z_centerline, y_centerline, deg=5)
    polyy = numpy.poly1d(coeffsy)
    y_centerline_fit = numpy.polyval(polyy, z_centerline)
    
    
    return x_centerline_fit, y_centerline_fit

Example 13

Project: battery-status Source File: battery-status-graph.py
def guess_expiry(x, y, zero = 0):
    fit = np.polyfit(data['energy_full'], data['timestamp'], 1)
    #print "fit: %s" % fit

    fit_fn = np.poly1d(fit)
    #print "fit_fn: %s" % fit_fn
    return datetime.datetime.fromtimestamp(fit_fn(zero))

Example 14

Project: seisflows Source File: math.py
def lsq2(x, f):
    # parabolic least squares fit
    p = np.polyfit(x, f, 2)
    if p[0] > 0:
        return -p[1]/(2*p[0])
    else:
        print -1
        raise Exception()

Example 15

Project: TheCannon Source File: snr_test.py
def quad_fit(x, y, yerr):
    """ Fit a qudratic to the SNR to make a lookup error table """
    qfit = np.polyfit(x, y, deg=2, w = 1 / yerr)
    plt.figure()
    plt.scatter(x, y)
    plt.errorbar(x, y, yerr=yerr)
    xvals = np.linspace(min(x), max(x), 100)
    yvals = qfit[2] + qfit[1]*xvals + qfit[0]*xvals**2
    plt.plot(xvals, yvals, color='r', lw=2)
    plt.show()
    print(qfit)

Example 16

Project: pymatgen Source File: elastic.py
    @classmethod
    def from_stress_dict(cls, stress_dict, tol=0.1, vasp=True, symmetry=False):
        """
        Constructs the elastic tensor from IndependentStrain-Stress dictionary
        corresponding to legacy behavior of elasticity package.

        Args:
            stress_dict (dict): dictionary of stresses indexed by corresponding
                IndependentStrain objects.
            tol (float): tolerance for zeroing small values of the tensor
            vasp (boolean): flag for whether the stress tensor should be
                converted based on vasp units/convention for stress
            symmetry (boolean): flag for whether or not the elastic tensor
                should fit from data based on symmetry
        """
        inds = [(0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)]
        c_ij = np.array([[np.polyfit([strain[ind1] for strain in list(stress_dict.keys())
                                      if (strain.i, strain.j) == ind1],
                                     [stress_dict[strain][ind2] for strain
                                      in list(stress_dict.keys())
                                      if (strain.i, strain.j) == ind1], 1)[0]
                          for ind1 in inds] for ind2 in inds])
        if vasp:
            c_ij *= -0.1  # Convert units/sign convention of vasp stress tensor
        c_ij[0:, 3:] = 0.5 * c_ij[0:, 3:]  # account for voigt doubling of e4,e5,e6
        c = cls.from_voigt(c_ij)
        c = c.zeroed()
        return c

Example 17

Project: distalgo Source File: plot.py
def load_graphline(graphline, ax=plt):
    """Plot a line on the graph, return its handle."""
    datafile = graphline.dataset.data_file
    results = load_data(datafile)
    xset = []
    yset = []
    for xaxis, datapoint in graphline.dataset:
        datas = [(data[graphline.key] if not graphline.avg_over_procs
                  else data[graphline.key] / data['Total_processes'])
                 / graphline.avg_factor
                 for config, data, ts in results if config == datapoint]
        if len(datas) == 0:
            print("No data for ", graphline.key, datapoint)
        xset.append(xaxis)
        yset.append(graphline.aggregate(datas))
    if graphline.fit_degree is not None:
        pol = np.polyfit(xset, yset, graphline.fit_degree)
        sample_xs = np.linspace(xset[0], xset[-1], len(xset) * 10)
        fitline_args = graphline.line_properties
        datapoint_args = graphline.line_properties
        fitline_args.pop('marker')
        datapoint_args.pop('linestyle')
        datapoint_args.pop('label')
        ax.plot(sample_xs, np.polyval(pol, sample_xs), **fitline_args)
        ax.plot(xset, yset, label='_nolegend_', linestyle="", **datapoint_args)
        return Line2D([0, 1], [0, 1], **graphline.line_properties)
    else:
        return ax.plot(xset, yset, **graphline.line_properties)

Example 18

Project: pymatgen Source File: eos_2.py
Function: deltafactor_polyfit
def deltafactor_polyfit(volumes, energies):
    """
    This is the routine used to compute V0, B0, B1 in the deltafactor code.
    """
    fitdata = np.polyfit(volumes**(-2./3.), energies, 3, full=True)

    deriv0 = np.poly1d(fitdata[0])
    deriv1 = np.polyder(deriv0, 1)
    deriv2 = np.polyder(deriv1, 1)
    deriv3 = np.polyder(deriv2, 1)

    for x in np.roots(deriv1):
        if x > 0 and deriv2(x) > 0:
            v0 = x**(-3./2.)
            break
    else:
        raise EOSError("No minimum could be found")

    derivV2 = 4./9. * x**5. * deriv2(x)
    derivV3 = (-20./9. * x**(13./2.) * deriv2(x) - 8./27. * x**(15./2.) * deriv3(x))
    b0 = derivV2 / x**(3./2.)
    b1 = -1 - x**(-3./2.) * derivV3 / derivV2

    # e0, b0, b1, v0, eos_params
    return np.poly1d(fitdata[0])(v0**(-2./3.)), b0, b1, v0, fitdata[0]

Example 19

Project: FaST-LMM Source File: feature_selection_cv.py
Function: fit_parabola
    def fit_parabola(self, deltas, perf, output_prefix=None, create_pdf=True):
        """
        for best k, fit parabola to 3 points in delta dimension and determine optimum accordingly,
        assume convex function
        """

        #self.run_once() - don't call this because we don't need snp data loaded

        assert len(perf) == len(deltas)
        assert perf[1] <= perf[0]
        assert perf[1] <= perf[2]
        assert deltas[0] < deltas[1] < deltas[2]

        coef = np.polyfit(deltas, perf, 2)

        xfit = np.linspace(min(deltas), max(deltas), num=100)
        yfit =  np.polyval(coef, xfit)

        best_idx = np.argmin(yfit)
        best_delta = xfit[best_idx]
        best_obj = yfit[best_idx]

        #if create_pdf and (output_prefix!=None):
        #    pylab.figure()
        #    pylab.plot(deltas, perf, '.')
        #    pylab.plot(xfit, yfit, '-')
        #    pylab.grid(True)
        #    plot_fn = output_prefix + "_parabola.pdf"
        #    util.create_directory_if_necessary(plot_fn)
        #    pylab.savefig(plot_fn)

        return best_delta, best_obj

Example 20

Project: jcvi Source File: synteny.py
Function: simple
def simple(args):
    """
    %prog simple anchorfile --qbed=qbedfile --sbed=sbedfile [options]

    Write the block ends for each block in the anchorfile.
    GeneA1    GeneA2    GeneB1    GeneB2   +/-      score

    Optional additional columns:
    orderA1   orderA2   orderB1   orderB2  sizeA    sizeB   size    block_id

    With base coordinates (--coords):
    block_id  seqidA    startA    endA     bpSpanA  GeneA1   GeneA2  geneSpanA
    block_id  seqidB    startB    endB     bpSpanB  GeneB1   GeneB2  geneSpanB
    """
    p = OptionParser(simple.__doc__)
    p.add_option("--rich", default=False, action="store_true", \
                help="Output additional columns [default: %default]")
    p.add_option("--coords", default=False, action="store_true",
                help="Output columns with base coordinates [default: %default]")
    p.add_option("--bed", default=False, action="store_true",
                help="Generate BED file for the blocks")
    p.add_option("--noheader", default=False, action="store_true",
                help="Don't output header [default: %default]")
    p.set_beds()
    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    anchorfile, = args
    additional = opts.rich
    coords = opts.coords
    header = not opts.noheader
    bed = opts.bed
    if bed:
        coords = True
        bbed = Bed()

    ac = AnchorFile(anchorfile)
    simplefile = anchorfile.rsplit(".", 1)[0] + ".simple"

    qbed, sbed, qorder, sorder, is_self = check_beds(anchorfile, p, opts)
    pf = "-".join(anchorfile.split(".", 2)[:2])
    blocks = ac.blocks

    if coords:
        h = "Block|Chr|Start|End|Span|StartGene|EndGene|GeneSpan|Orientation"
    else:
        h = "StartGeneA|EndGeneA|StartGeneB|EndGeneB|Orientation|Score"
        if additional:
            h += "|StartOrderA|EndOrderA|StartOrderB|EndOrderB|"\
                  "SizeA|SizeB|Size|Block"

    fws = open(simplefile, "w")
    if header:
        print >> fws, "\t".join(h.split("|"))

    atotalbase = btotalbase = 0
    for i, block in enumerate(blocks):

        a, b, scores = zip(*block)
        a = [qorder[x] for x in a]
        b = [sorder[x] for x in b]
        ia, oa = zip(*a)
        ib, ob = zip(*b)

        astarti, aendi = min(ia), max(ia)
        bstarti, bendi = min(ib), max(ib)
        astart, aend = min(a)[1].accn, max(a)[1].accn
        bstart, bend = min(b)[1].accn, max(b)[1].accn

        sizeA = len(set(ia))
        sizeB = len(set(ib))
        size = len(block)

        slope, intercept = np.polyfit(ia, ib, 1)
        orientation = "+" if slope >= 0 else '-'
        aspan = aendi - astarti + 1
        bspan = bendi - bstarti + 1
        score = int((aspan * bspan) ** .5)
        score = str(score)
        block_id = pf + "-block-{0}".format(i)

        if coords:

            aseqid, astartbase, aendbase = \
                    get_boundary_bases(astart, aend, qorder)
            bseqid, bstartbase, bendbase = \
                    get_boundary_bases(bstart, bend, sorder)
            abase = aendbase - astartbase + 1
            bbase = bendbase - bstartbase + 1
            atotalbase += abase
            btotalbase += bbase

            # Write dual lines
            aargs = [block_id, aseqid, astartbase, aendbase,
                     abase, astart, aend, aspan, "+"]
            bargs = [block_id, bseqid, bstartbase, bendbase,
                     bbase, bstart, bend, bspan, orientation]

            if bed:
                bbed.append(BedLine("\t".join(str(x) for x in \
                           (bseqid, bstartbase - 1, bendbase,
                           "{}:{}-{}".format(aseqid, astartbase, aendbase),
                           size, orientation))))

            for args in (aargs, bargs):
                print >> fws, "\t".join(str(x) for x in args)
            continue

        args = [astart, aend, bstart, bend, score, orientation]
        if additional:
            args += [astarti, aendi, bstarti, bendi,
                     sizeA, sizeB, size, block_id]
        print >> fws, "\t".join(str(x) for x in args)

    fws.close()
    logging.debug("A total of {0} blocks written to `{1}`.".format(i + 1, simplefile))

    if coords:
        print >> sys.stderr, "Total block span in {0}: {1}".format(qbed.filename, \
                        human_size(atotalbase, precision=2))
        print >> sys.stderr, "Total block span in {0}: {1}".format(sbed.filename, \
                        human_size(btotalbase, precision=2))
        print >> sys.stderr, "Ratio: {0:.1f}x".format(\
                        max(atotalbase, btotalbase) * 1. / min(atotalbase, btotalbase))

    if bed:
        bedfile = simplefile + ".bed"
        bbed.print_to_file(filename=bedfile, sorted=True)
        logging.debug("Bed file written to `{}`".format(bedfile))

Example 21

Project: alchemical-analysis Source File: parser_amber.py
def _extrapol(x, y, scheme):
    """Simple extrapolation schemes."""

    y0 = None
    y1 = None

    if scheme == 'linear':
        if 0.0 not in x:
            y0 = (x[0] * y[1] - x[1] * y[0]) / (x[0] - x[1])

        if 1.0 not in x:
            y1 = (((x[-2] - 1.0) * y[-1] + ((1.0 - x[-1]) * y[-2])) /
                  (x[-2] - x[-1]))
    elif scheme == 'polyfit':
        nval = len(x)

        if nval < 6:
            deg = nval - 1
        else:
            deg = 6

        coeffs = np.polyfit(x, y, deg)

        if 0.0 not in x:
            y0 = coeffs[-1]

        if 1.0 not in x:
            y1 = sum(coeffs)
    else:
        raise SystemExit('ERROR: Unsupported extrapolation scheme: %s' %
                         scheme)

    return y0, y1

Example 22

Project: AWS-Lambda-ML-Microservice-Skeleton Source File: test_regression.py
Function: test_polyfit_build
    def test_polyfit_build(self):
        # Ticket #628
        ref = [-1.06123820e-06, 5.70886914e-04, -1.13822012e-01,
               9.95368241e+00, -3.14526520e+02]
        x = [90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
             104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115,
             116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 129,
             130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
             146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157,
             158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
             170, 171, 172, 173, 174, 175, 176]
        y = [9.0, 3.0, 7.0, 4.0, 4.0, 8.0, 6.0, 11.0, 9.0, 8.0, 11.0, 5.0,
             6.0, 5.0, 9.0, 8.0, 6.0, 10.0, 6.0, 10.0, 7.0, 6.0, 6.0, 6.0,
             13.0, 4.0, 9.0, 11.0, 4.0, 5.0, 8.0, 5.0, 7.0, 7.0, 6.0, 12.0,
             7.0, 7.0, 9.0, 4.0, 12.0, 6.0, 6.0, 4.0, 3.0, 9.0, 8.0, 8.0,
             6.0, 7.0, 9.0, 10.0, 6.0, 8.0, 4.0, 7.0, 7.0, 10.0, 8.0, 8.0,
             6.0, 3.0, 8.0, 4.0, 5.0, 7.0, 8.0, 6.0, 6.0, 4.0, 12.0, 9.0,
             8.0, 8.0, 8.0, 6.0, 7.0, 4.0, 4.0, 5.0, 7.0]
        tested = np.polyfit(x, y, 4)
        assert_array_almost_equal(ref, tested)

Example 23

Project: workload-automation Source File: __init__.py
def fit_polynomial(s, n):
    # pylint: disable=no-member
    coeffs = np.polyfit(s.index, s.values, n)
    poly = np.poly1d(coeffs)
    return poly(s.index)

Example 24

Project: thunder Source File: series.py
Function: detrend
    def detrend(self, method='linear', order=5):
        """
        Detrend series data with linear or nonlinear detrending.

        Preserve intercept so that subsequent operations can adjust the baseline.

        Parameters
        ----------
        method : str, optional, default = 'linear'
            Detrending method

        order : int, optional, default = 5
            Order of polynomial, for non-linear detrending only
        """
        check_options(method, ['linear', 'nonlinear'])

        if method == 'linear':
            order = 1

        def func(y):
            x = arange(len(y))
            p = polyfit(x, y, order)
            p[-1] = 0
            yy = polyval(p, x)
            return y - yy

        return self.map(func)

Example 25

Project: spinalcordtoolbox Source File: sct_get_centerline_automatic.py
def main():

    # Initialization
    fname_anat = ''
    fname_point = ''
    slice_gap = param.gap
    remove_tmp_files = param.remove_tmp_files
    gaussian_kernel = param.gaussian_kernel
    start_time = time.time()
    verbose = 1

    # get path of the toolbox
    status, path_sct = commands.getstatusoutput('echo $SCT_DIR')
    path_sct = sct.slash_at_the_end(path_sct, 1)

    # Parameters for debug mode
    if param.debug == 1:
        sct.printv('\n*** WARNING: DEBUG MODE ON ***\n\t\t\tCurrent working directory: '+os.getcwd(), 'warning')
        status, path_sct_testing_data = commands.getstatusoutput('echo $SCT_TESTING_DATA_DIR')
        fname_anat = path_sct_testing_data+'/t2/t2.nii.gz'
        fname_point = path_sct_testing_data+'/t2/t2_centerline_init.nii.gz'
        slice_gap = 5

    else:
        # Check input param
        try:
            opts, args = getopt.getopt(sys.argv[1:],'hi:p:g:r:k:')
        except getopt.GetoptError as err:
            print str(err)
            usage()
        if not opts:
            usage()
        for opt, arg in opts:
            if opt == '-h':
                usage()
            elif opt in ('-i'):
                fname_anat = arg
            elif opt in ('-p'):
                fname_point = arg
            elif opt in ('-g'):
                slice_gap = int(arg)
            elif opt in ('-r'):
                remove_tmp_files = int(arg)
            elif opt in ('-k'):
                gaussian_kernel = int(arg)

    # display usage if a mandatory argument is not provided
    if fname_anat == '' or fname_point == '':
        usage()

    # check existence of input files
    sct.check_file_exist(fname_anat)
    sct.check_file_exist(fname_point)

    # extract path/file/extension
    path_anat, file_anat, ext_anat = sct.extract_fname(fname_anat)
    path_point, file_point, ext_point = sct.extract_fname(fname_point)

    # extract path of schedule file
    # TODO: include schedule file in sct
    # TODO: check existence of schedule file
    file_schedule = path_sct + param.schedule_file

    # Get input image orientation
    input_image_orientation = get_orientation(fname_anat)

    # Display arguments
    print '\nCheck input arguments...'
    print '  Anatomical image:     '+fname_anat
    print '  Orientation:          '+input_image_orientation
    print '  Point in spinal cord: '+fname_point
    print '  Slice gap:            '+str(slice_gap)
    print '  Gaussian kernel:      '+str(gaussian_kernel)
    print '  Degree of polynomial: '+str(param.deg_poly)

    # create temporary folder
    print('\nCreate temporary folder...')
    path_tmp = 'tmp.'+time.strftime("%y%m%d%H%M%S")
    sct.create_folder(path_tmp)
    print '\nCopy input data...'
    sct.run('cp '+fname_anat+ ' '+path_tmp+'/tmp.anat'+ext_anat)
    sct.run('cp '+fname_point+ ' '+path_tmp+'/tmp.point'+ext_point)

    # go to temporary folder
    os.chdir(path_tmp)

    # convert to nii
    convert('tmp.anat'+ext_anat, 'tmp.anat.nii')
    convert('tmp.point'+ext_point, 'tmp.point.nii')

    # Reorient input anatomical volume into RL PA IS orientation
    print '\nReorient input volume to RL PA IS orientation...'
    #sct.run(sct.fsloutput + 'fslswapdim tmp.anat RL PA IS tmp.anat_orient')
    set_orientation('tmp.anat.nii', 'RPI', 'tmp.anat_orient.nii')
    # Reorient binary point into RL PA IS orientation
    print '\nReorient binary point into RL PA IS orientation...'
    # sct.run(sct.fsloutput + 'fslswapdim tmp.point RL PA IS tmp.point_orient')
    set_orientation('tmp.point.nii', 'RPI', 'tmp.point_orient.nii')

    # Get image dimensions
    print '\nGet image dimensions...'
    nx, ny, nz, nt, px, py, pz, pt = Image('tmp.anat_orient.nii').dim
    print '.. matrix size: '+str(nx)+' x '+str(ny)+' x '+str(nz)
    print '.. voxel size:  '+str(px)+'mm x '+str(py)+'mm x '+str(pz)+'mm'

    # Split input volume
    print '\nSplit input volume...'
    split_data('tmp.anat_orient.nii', 2, '_z')
    file_anat_split = ['tmp.anat_orient_z'+str(z).zfill(4) for z in range(0, nz, 1)]
    split_data('tmp.point_orient.nii', 2, '_z')
    file_point_split = ['tmp.point_orient_z'+str(z).zfill(4) for z in range(0, nz, 1)]

    # Extract coordinates of input point
    # sct.printv('\nExtract the slice corresponding to z='+str(z_init)+'...', verbose)
    #
    data_point = Image('tmp.point_orient.nii').data
    x_init, y_init, z_init = unravel_index(data_point.argmax(), data_point.shape)
    sct.printv('Coordinates of input point: ('+str(x_init)+', '+str(y_init)+', '+str(z_init)+')', verbose)

    # Create 2D gaussian mask
    sct.printv('\nCreate gaussian mask from point...', verbose)
    xx, yy = mgrid[:nx, :ny]
    mask2d = zeros((nx, ny))
    radius = round(float(gaussian_kernel+1)/2)  # add 1 because the radius includes the center.
    sigma = float(radius)
    mask2d = exp(-(((xx-x_init)**2)/(2*(sigma**2)) + ((yy-y_init)**2)/(2*(sigma**2))))

    # Save mask to 2d file
    file_mask_split = ['tmp.mask_orient_z'+str(z).zfill(4) for z in range(0,nz,1)]
    nii_mask2d = Image('tmp.anat_orient_z0000.nii')
    nii_mask2d.data = mask2d
    nii_mask2d.setFileName(file_mask_split[z_init]+'.nii')
    nii_mask2d.save()
    #
    # # Get the coordinates of the input point
    # print '\nGet the coordinates of the input point...'
    # data_point = Image('tmp.point_orient.nii').data
    # x_init, y_init, z_init = unravel_index(data_point.argmax(), data_point.shape)
    # print '('+str(x_init)+', '+str(y_init)+', '+str(z_init)+')'

    # x_init, y_init, z_init = (data > 0).nonzero()
    # x_init = x_init[0]
    # y_init = y_init[0]
    # z_init = z_init[0]
    # print '('+str(x_init)+', '+str(y_init)+', '+str(z_init)+')'
    #
    # numpy.unravel_index(a.argmax(), a.shape)
    #
    # file = nibabel.load('tmp.point_orient.nii')
    # data = file.get_data()
    # x_init, y_init, z_init = (data > 0).nonzero()
    # x_init = x_init[0]
    # y_init = y_init[0]
    # z_init = z_init[0]
    # print '('+str(x_init)+', '+str(y_init)+', '+str(z_init)+')'
    #
    # # Extract the slice corresponding to z=z_init
    # print '\nExtract the slice corresponding to z='+str(z_init)+'...'
    # file_point_split = ['tmp.point_orient_z'+str(z).zfill(4) for z in range(0,nz,1)]
    # nii = Image('tmp.point_orient.nii')
    # data_crop = nii.data[:, :, z_init:z_init+1]
    # nii.data = data_crop
    # nii.setFileName(file_point_split[z_init]+'.nii')
    # nii.save()
    #
    # # Create gaussian mask from point
    # print '\nCreate gaussian mask from point...'
    # file_mask_split = ['tmp.mask_orient_z'+str(z).zfill(4) for z in range(0,nz,1)]
    # sct.run(sct.fsloutput+'fslmaths '+file_point_split[z_init]+' -s '+str(gaussian_kernel)+' '+file_mask_split[z_init])
    #
    # # Obtain max value from mask
    # print '\nFind maximum value from mask...'
    # file = nibabel.load(file_mask_split[z_init]+'.nii')
    # data = file.get_data()
    # max_value_mask = numpy.max(data)
    # print '..'+str(max_value_mask)
    #
    # # Normalize mask beween 0 and 1
    # print '\nNormalize mask beween 0 and 1...'
    # sct.run(sct.fsloutput+'fslmaths '+file_mask_split[z_init]+' -div '+str(max_value_mask)+' '+file_mask_split[z_init])

    ## Take the square of the mask
    #print '\nCalculate the square of the mask...'
    #sct.run(sct.fsloutput+'fslmaths '+file_mask_split[z_init]+' -mul '+file_mask_split[z_init]+' '+file_mask_split[z_init])

    # initialize variables
    file_mat = ['tmp.mat_z'+str(z).zfill(4) for z in range(0,nz,1)]
    file_mat_inv = ['tmp.mat_inv_z'+str(z).zfill(4) for z in range(0,nz,1)]
    file_mat_inv_cuemul = ['tmp.mat_inv_cuemul_z'+str(z).zfill(4) for z in range(0,nz,1)]

    # create identity matrix for initial transformation matrix
    fid = open(file_mat_inv_cuemul[z_init], 'w')
    fid.write('%i %i %i %i\n' %(1, 0, 0, 0) )
    fid.write('%i %i %i %i\n' %(0, 1, 0, 0) )
    fid.write('%i %i %i %i\n' %(0, 0, 1, 0) )
    fid.write('%i %i %i %i\n' %(0, 0, 0, 1) )
    fid.close()

    # initialize centerline: give value corresponding to initial point
    x_centerline = [x_init]
    y_centerline = [y_init]
    z_centerline = [z_init]
    warning_count = 0

    # go up (1), then down (2) in reference to the binary point
    for iUpDown in range(1, 3):

        if iUpDown == 1:
            # z increases
            slice_gap_signed = slice_gap
        elif iUpDown == 2:
            # z decreases
            slice_gap_signed = -slice_gap
            # reverse centerline (because values will be appended at the end)
            x_centerline.reverse()
            y_centerline.reverse()
            z_centerline.reverse()

        # initialization before looping
        z_dest = z_init # point given by user
        z_src = z_dest + slice_gap_signed

        # continue looping if 0 < z < nz
        while 0 <= z_src and z_src <= nz-1:

            # print current z:
            print 'z='+str(z_src)+':'

            # estimate transformation
            sct.run(fsloutput+'flirt -in '+file_anat_split[z_src]+' -ref '+file_anat_split[z_dest]+' -schedule '+file_schedule+ ' -verbose 0 -omat '+file_mat[z_src]+' -cost normcorr -forcescaling -inweight '+file_mask_split[z_dest]+' -refweight '+file_mask_split[z_dest])

            # display transfo
            status, output = sct.run('cat '+file_mat[z_src])
            print output

            # check if transformation is bigger than 1.5x slice_gap
            tx = float(output.split()[3])
            ty = float(output.split()[7])
            norm_txy = linalg.norm([tx, ty],ord=2)
            if norm_txy > 1.5*slice_gap:
                print 'WARNING: Transformation is too large --> using previous one.'
                warning_count = warning_count + 1
                # if previous transformation exists, replace current one with previous one
                if os.path.isfile(file_mat[z_dest]):
                    sct.run('cp '+file_mat[z_dest]+' '+file_mat[z_src])

            # estimate inverse transformation matrix
            sct.run('convert_xfm -omat '+file_mat_inv[z_src]+' -inverse '+file_mat[z_src])

            # compute cuemulative transformation
            sct.run('convert_xfm -omat '+file_mat_inv_cuemul[z_src]+' -concat '+file_mat_inv[z_src]+' '+file_mat_inv_cuemul[z_dest])

            # apply inverse cuemulative transformation to initial gaussian mask (to put it in src space)
            sct.run(fsloutput+'flirt -in '+file_mask_split[z_init]+' -ref '+file_mask_split[z_init]+' -applyxfm -init '+file_mat_inv_cuemul[z_src]+' -out '+file_mask_split[z_src])

            # open inverse cuemulative transformation file and generate centerline
            fid = open(file_mat_inv_cuemul[z_src])
            mat = fid.read().split()
            x_centerline.append(x_init + float(mat[3]))
            y_centerline.append(y_init + float(mat[7]))
            z_centerline.append(z_src)
            #z_index = z_index+1

            # define new z_dest (target slice) and new z_src (moving slice)
            z_dest = z_dest + slice_gap_signed
            z_src = z_src + slice_gap_signed


    # Reconstruct centerline
    # ====================================================================================================

    # reverse back centerline (because it's been reversed once, so now all values are in the right order)
    x_centerline.reverse()
    y_centerline.reverse()
    z_centerline.reverse()

    # fit centerline in the Z-X plane using polynomial function
    print '\nFit centerline in the Z-X plane using polynomial function...'
    coeffsx = polyfit(z_centerline, x_centerline, deg=param.deg_poly)
    polyx = poly1d(coeffsx)
    x_centerline_fit = polyval(polyx, z_centerline)
    # calculate RMSE
    rmse = linalg.norm(x_centerline_fit-x_centerline)/sqrt( len(x_centerline) )
    # calculate max absolute error
    max_abs = max( abs(x_centerline_fit-x_centerline) )
    print '.. RMSE (in mm): '+str(rmse*px)
    print '.. Maximum absolute error (in mm): '+str(max_abs*px)

    # fit centerline in the Z-Y plane using polynomial function
    print '\nFit centerline in the Z-Y plane using polynomial function...'
    coeffsy = polyfit(z_centerline, y_centerline, deg=param.deg_poly)
    polyy = poly1d(coeffsy)
    y_centerline_fit = polyval(polyy, z_centerline)
    # calculate RMSE
    rmse = linalg.norm(y_centerline_fit-y_centerline)/sqrt( len(y_centerline) )
    # calculate max absolute error
    max_abs = max( abs(y_centerline_fit-y_centerline) )
    print '.. RMSE (in mm): '+str(rmse*py)
    print '.. Maximum absolute error (in mm): '+str(max_abs*py)

    # display
    if param.debug == 1:
        import matplotlib.pyplot as plt
        plt.figure()
        plt.plot(z_centerline,x_centerline,'.',z_centerline,x_centerline_fit,'r')
        plt.legend(['Data','Polynomial Fit'])
        plt.title('Z-X plane polynomial interpolation')
        plt.show()

        plt.figure()
        plt.plot(z_centerline,y_centerline,'.',z_centerline,y_centerline_fit,'r')
        plt.legend(['Data','Polynomial Fit'])
        plt.title('Z-Y plane polynomial interpolation')
        plt.show()

    # generate full range z-values for centerline
    z_centerline_full = [iz for iz in range(0, nz, 1)]

    # calculate X and Y values for the full centerline
    x_centerline_fit_full = polyval(polyx, z_centerline_full)
    y_centerline_fit_full = polyval(polyy, z_centerline_full)

    # Generate fitted transformation matrices and write centerline coordinates in text file
    print '\nGenerate fitted transformation matrices and write centerline coordinates in text file...'
    file_mat_inv_cuemul_fit = ['tmp.mat_inv_cuemul_fit_z'+str(z).zfill(4) for z in range(0,nz,1)]
    file_mat_cuemul_fit = ['tmp.mat_cuemul_fit_z'+str(z).zfill(4) for z in range(0,nz,1)]
    fid_centerline = open('tmp.centerline_coordinates.txt', 'w')
    for iz in range(0, nz, 1):
        # compute inverse cuemulative fitted transformation matrix
        fid = open(file_mat_inv_cuemul_fit[iz], 'w')
        fid.write('%i %i %i %f\n' %(1, 0, 0, x_centerline_fit_full[iz]-x_init) )
        fid.write('%i %i %i %f\n' %(0, 1, 0, y_centerline_fit_full[iz]-y_init) )
        fid.write('%i %i %i %i\n' %(0, 0, 1, 0) )
        fid.write('%i %i %i %i\n' %(0, 0, 0, 1) )
        fid.close()
        # compute forward cuemulative fitted transformation matrix
        sct.run('convert_xfm -omat '+file_mat_cuemul_fit[iz]+' -inverse '+file_mat_inv_cuemul_fit[iz])
        # write centerline coordinates in x, y, z format
        fid_centerline.write('%f %f %f\n' %(x_centerline_fit_full[iz], y_centerline_fit_full[iz], z_centerline_full[iz]) )
    fid_centerline.close()


    # Prepare output data
    # ====================================================================================================

    # write centerline as text file
    for iz in range(0, nz, 1):
        # compute inverse cuemulative fitted transformation matrix
        fid = open(file_mat_inv_cuemul_fit[iz], 'w')
        fid.write('%i %i %i %f\n' %(1, 0, 0, x_centerline_fit_full[iz]-x_init) )
        fid.write('%i %i %i %f\n' %(0, 1, 0, y_centerline_fit_full[iz]-y_init) )
        fid.write('%i %i %i %i\n' %(0, 0, 1, 0) )
        fid.write('%i %i %i %i\n' %(0, 0, 0, 1) )
        fid.close()

    # write polynomial coefficients
    savetxt('tmp.centerline_polycoeffs_x.txt',coeffsx)
    savetxt('tmp.centerline_polycoeffs_y.txt',coeffsy)

    # apply transformations to data
    print '\nApply fitted transformation matrices...'
    file_anat_split_fit = ['tmp.anat_orient_fit_z'+str(z).zfill(4) for z in range(0,nz,1)]
    file_mask_split_fit = ['tmp.mask_orient_fit_z'+str(z).zfill(4) for z in range(0,nz,1)]
    file_point_split_fit = ['tmp.point_orient_fit_z'+str(z).zfill(4) for z in range(0,nz,1)]
    for iz in range(0, nz, 1):
        # forward cuemulative transformation to data
        sct.run(fsloutput+'flirt -in '+file_anat_split[iz]+' -ref '+file_anat_split[iz]+' -applyxfm -init '+file_mat_cuemul_fit[iz]+' -out '+file_anat_split_fit[iz])
        # inverse cuemulative transformation to mask
        sct.run(fsloutput+'flirt -in '+file_mask_split[z_init]+' -ref '+file_mask_split[z_init]+' -applyxfm -init '+file_mat_inv_cuemul_fit[iz]+' -out '+file_mask_split_fit[iz])
        # inverse cuemulative transformation to point
        sct.run(fsloutput+'flirt -in '+file_point_split[z_init]+' -ref '+file_point_split[z_init]+' -applyxfm -init '+file_mat_inv_cuemul_fit[iz]+' -out '+file_point_split_fit[iz]+' -interp nearestneighbour')

    # Merge into 4D volume
    print '\nMerge into 4D volume...'
    # sct.run(fsloutput+'fslmerge -z tmp.anat_orient_fit tmp.anat_orient_fit_z*')
    # sct.run(fsloutput+'fslmerge -z tmp.mask_orient_fit tmp.mask_orient_fit_z*')
    # sct.run(fsloutput+'fslmerge -z tmp.point_orient_fit tmp.point_orient_fit_z*')
    concat_data(glob.glob('tmp.anat_orient_fit_z*.nii'), 'tmp.anat_orient_fit.nii', dim=2)
    concat_data(glob.glob('tmp.mask_orient_fit_z*.nii'), 'tmp.mask_orient_fit.nii', dim=2)
    concat_data(glob.glob('tmp.point_orient_fit_z*.nii'), 'tmp.point_orient_fit.nii', dim=2)

    # Copy header geometry from input data
    print '\nCopy header geometry from input data...'
    copy_header('tmp.anat_orient.nii', 'tmp.anat_orient_fit.nii')
    copy_header('tmp.anat_orient.nii', 'tmp.mask_orient_fit.nii')
    copy_header('tmp.anat_orient.nii', 'tmp.point_orient_fit.nii')

    # Reorient outputs into the initial orientation of the input image
    print '\nReorient the centerline into the initial orientation of the input image...'
    set_orientation('tmp.point_orient_fit.nii', input_image_orientation, 'tmp.point_orient_fit.nii')
    set_orientation('tmp.mask_orient_fit.nii', input_image_orientation, 'tmp.mask_orient_fit.nii')

    # Generate output file (in current folder)
    print '\nGenerate output file (in current folder)...'
    os.chdir('..')  # come back to parent folder
    #sct.generate_output_file('tmp.centerline_polycoeffs_x.txt','./','centerline_polycoeffs_x','.txt')
    #sct.generate_output_file('tmp.centerline_polycoeffs_y.txt','./','centerline_polycoeffs_y','.txt')
    #sct.generate_output_file('tmp.centerline_coordinates.txt','./','centerline_coordinates','.txt')
    #sct.generate_output_file('tmp.anat_orient.nii','./',file_anat+'_rpi',ext_anat)
    #sct.generate_output_file('tmp.anat_orient_fit.nii', file_anat+'_rpi_align'+ext_anat)
    #sct.generate_output_file('tmp.mask_orient_fit.nii', file_anat+'_mask'+ext_anat)
    fname_output_centerline = sct.generate_output_file(path_tmp+'/tmp.point_orient_fit.nii', file_anat+'_centerline'+ext_anat)

    # Delete temporary files
    if remove_tmp_files == 1:
        print '\nRemove temporary files...'
        sct.run('rm -rf '+path_tmp)

    # print number of warnings
    print '\nNumber of warnings: '+str(warning_count)+' (if >10, you should probably reduce the gap and/or increase the kernel size'

    # display elapsed time
    elapsed_time = time.time() - start_time
    print '\nFinished! \n\tGenerated file: '+fname_output_centerline+'\n\tElapsed time: '+str(int(round(elapsed_time)))+'s\n'

Example 26

Project: compare-codecs Source File: graph_metrics.py
Function: bdrate
def BdRate(group1, group2):
  """Compute the BD-rate between two score groups.

  The returned object also contains the range of PSNR values used
  to compute the result.

  Bjontegaard's metric allows to compute the average % saving in bitrate
  between two rate-distortion curves [1].

  rate1,psnr1 - RD points for curve 1
  rate2,psnr2 - RD points for curve 2

  adapted from code from: (c) 2010 Giuseppe Valenzise
  copied from code by [email protected], [email protected]

  """
  # pylint: disable=too-many-locals
  metric_set1 = group1.dataPoints()
  metric_set2 = group2.dataPoints()

  # numpy plays games with its exported functions.
  # pylint: disable=no-member
  # pylint: disable=bad-builtin
  psnr1 = [x[1] for x in metric_set1]
  psnr2 = [x[1] for x in metric_set2]

  log_rate1 = map(math.log, [x[0] for x in metric_set1])
  log_rate2 = map(math.log, [x[0] for x in metric_set2])

  # Best cubic poly fit for graph represented by log_ratex, psrn_x.
  poly1 = numpy.polyfit(psnr1, log_rate1, 3)
  poly2 = numpy.polyfit(psnr2, log_rate2, 3)

  # Integration interval.
  min_int = max([min(psnr1), min(psnr2)])
  max_int = min([max(psnr1), max(psnr2)])

  # find integral
  p_int1 = numpy.polyint(poly1)
  p_int2 = numpy.polyint(poly2)

  # Calculate the integrated value over the interval we care about.
  int1 = numpy.polyval(p_int1, max_int) - numpy.polyval(p_int1, min_int)
  int2 = numpy.polyval(p_int2, max_int) - numpy.polyval(p_int2, min_int)

  # Calculate the average improvement.
  avg_exp_diff = (int2 - int1) / (max_int - min_int)

  # In really bad formed data the exponent can grow too large.
  # clamp it.
  if avg_exp_diff > 200:
    avg_exp_diff = 200

  # Convert to a percentage.
  avg_diff = (math.exp(avg_exp_diff) - 1) * 100

  return {'difference': avg_diff, 'psnr':[min_int, max_int]}

Example 27

Project: iris Source File: polynomial_fit.py
def main():
    # Enable a future option, to ensure that the netcdf load works the same way
    # as in future Iris versions.
    iris.FUTURE.netcdf_promote = True

    # Load some test data.
    fname = iris.sample_data_path('A1B_north_america.nc')
    cube = iris.load_cube(fname)

    # Extract a single time series at a latitude and longitude point.
    location = next(cube.slices(['time']))

    # Calculate a polynomial fit to the data at this time series.
    x_points = location.coord('time').points
    y_points = location.data
    degree = 2

    p = np.polyfit(x_points, y_points, degree)
    y_fitted = np.polyval(p, x_points)

    # Add the polynomial fit values to the time series to take
    # full advantage of Iris plotting functionality.
    long_name = 'degree_{}_polynomial_fit_of_{}'.format(degree, cube.name())
    fit = iris.coords.AuxCoord(y_fitted, long_name=long_name,
                               units=location.units)
    location.add_aux_coord(fit, 0)

    qplt.plot(location.coord('time'), location, label='data')
    qplt.plot(location.coord('time'),
              location.coord(long_name),
              'g-', label='polynomial fit')
    plt.legend(loc='best')
    plt.title('Trend of US air temperature over time')

    qplt.show()

Example 28

Project: compare-codecs Source File: visual_metrics.py
def bdsnr(metric_set1, metric_set2):
  """
  BJONTEGAARD    Bjontegaard metric calculation
  Bjontegaard's metric allows to compute the average gain in psnr between two
  rate-distortion curves [1].
  rate1,psnr1 - RD points for curve 1
  rate2,psnr2 - RD points for curve 2

  returns the calculated Bjontegaard metric 'dsnr'

  code adapted from code written by : (c) 2010 Giuseppe Valenzise
  http://www.mathworks.com/matlabcentral/fileexchange/27798-bjontegaard-metric/content/bjontegaard.m
  """
  # pylint: disable=too-many-locals
  # numpy seems to do tricks with its exports.
  # pylint: disable=no-member
  # map() is recommended against.
  # pylint: disable=bad-builtin
  rate1 = [x[0] for x in metric_set1]
  psnr1 = [x[1] for x in metric_set1]
  rate2 = [x[0] for x in metric_set2]
  psnr2 = [x[1] for x in metric_set2]

  log_rate1 = map(math.log, rate1)
  log_rate2 = map(math.log, rate2)

  # Best cubic poly fit for graph represented by log_ratex, psrn_x.
  poly1 = numpy.polyfit(log_rate1, psnr1, 3)
  poly2 = numpy.polyfit(log_rate2, psnr2, 3)

  # Integration interval.
  min_int = max([min(log_rate1), min(log_rate2)])
  max_int = min([max(log_rate1), max(log_rate2)])

  # Integrate poly1, and poly2.
  p_int1 = numpy.polyint(poly1)
  p_int2 = numpy.polyint(poly2)

  # Calculate the integrated value over the interval we care about.
  int1 = numpy.polyval(p_int1, max_int) - numpy.polyval(p_int1, min_int)
  int2 = numpy.polyval(p_int2, max_int) - numpy.polyval(p_int2, min_int)

  # Calculate the average improvement.
  if max_int != min_int:
    avg_diff = (int2 - int1) / (max_int - min_int)
  else:
    avg_diff = 0.0
  return avg_diff

Example 29

Project: spinalcordtoolbox Source File: sct_get_centerline.py
def get_centerline_from_point(input_image, point_file, gap=4, gaussian_kernel=4, remove_tmp_files=1):

    # Initialization
    fname_anat = input_image
    fname_point = point_file
    slice_gap = gap
    remove_tmp_files = remove_tmp_files
    gaussian_kernel = gaussian_kernel
    start_time = time()
    verbose = 1

    # get path of the toolbox
    status, path_sct = commands.getstatusoutput('echo $SCT_DIR')
    path_sct = sct.slash_at_the_end(path_sct, 1)

    # Parameters for debug mode
    if param.debug == 1:
        sct.printv('\n*** WARNING: DEBUG MODE ON ***\n\t\t\tCurrent working directory: '+os.getcwd(), 'warning')
        status, path_sct_testing_data = commands.getstatusoutput('echo $SCT_TESTING_DATA_DIR')
        fname_anat = path_sct_testing_data+'/t2/t2.nii.gz'
        fname_point = path_sct_testing_data+'/t2/t2_centerline_init.nii.gz'
        slice_gap = 5

    # check existence of input files
    sct.check_file_exist(fname_anat)
    sct.check_file_exist(fname_point)

    # extract path/file/extension
    path_anat, file_anat, ext_anat = sct.extract_fname(fname_anat)
    path_point, file_point, ext_point = sct.extract_fname(fname_point)

    # extract path of schedule file
    # TODO: include schedule file in sct
    # TODO: check existence of schedule file
    file_schedule = path_sct + param.schedule_file

    # Get input image orientation
    input_image_orientation = get_orientation_3d(fname_anat, filename=True)

    # Display arguments
    print '\nCheck input arguments...'
    print '  Anatomical image:     '+fname_anat
    print '  Orientation:          '+input_image_orientation
    print '  Point in spinal cord: '+fname_point
    print '  Slice gap:            '+str(slice_gap)
    print '  Gaussian kernel:      '+str(gaussian_kernel)
    print '  Degree of polynomial: '+str(param.deg_poly)

    # create temporary folder
    print('\nCreate temporary folder...')
    path_tmp = 'tmp.'+strftime('%y%m%d%H%M%S')
    sct.create_folder(path_tmp)
    print '\nCopy input data...'
    sct.run('cp '+fname_anat+ ' '+path_tmp+'/tmp.anat'+ext_anat)
    sct.run('cp '+fname_point+ ' '+path_tmp+'/tmp.point'+ext_point)

    # go to temporary folder
    os.chdir(path_tmp)

    # convert to nii
    im_anat = convert('tmp.anat'+ext_anat, 'tmp.anat.nii')
    im_point = convert('tmp.point'+ext_point, 'tmp.point.nii')

    # Reorient input anatomical volume into RL PA IS orientation
    print '\nReorient input volume to RL PA IS orientation...'
    set_orientation(im_anat, 'RPI')
    im_anat.setFileName('tmp.anat_orient.nii')
    # Reorient binary point into RL PA IS orientation
    print '\nReorient binary point into RL PA IS orientation...'
    # sct.run(sct.fsloutput + 'fslswapdim tmp.point RL PA IS tmp.point_orient')
    set_orientation(im_point, 'RPI')
    im_point.setFileName('tmp.point_orient.nii')

    # Get image dimensions
    print '\nGet image dimensions...'
    nx, ny, nz, nt, px, py, pz, pt = Image('tmp.anat_orient.nii').dim
    print '.. matrix size: '+str(nx)+' x '+str(ny)+' x '+str(nz)
    print '.. voxel size:  '+str(px)+'mm x '+str(py)+'mm x '+str(pz)+'mm'

    # Split input volume
    print '\nSplit input volume...'
    im_anat_split_list = split_data(im_anat, 2)
    file_anat_split = []
    for im in im_anat_split_list:
        file_anat_split.append(im.absolutepath)
        im.save()

    im_point_split_list = split_data(im_point, 2)
    file_point_split = []
    for im in im_point_split_list:
        file_point_split.append(im.absolutepath)
        im.save()

    # Extract coordinates of input point
    data_point = Image('tmp.point_orient.nii').data
    x_init, y_init, z_init = unravel_index(data_point.argmax(), data_point.shape)
    sct.printv('Coordinates of input point: ('+str(x_init)+', '+str(y_init)+', '+str(z_init)+')', verbose)

    # Create 2D gaussian mask
    sct.printv('\nCreate gaussian mask from point...', verbose)
    xx, yy = mgrid[:nx, :ny]
    mask2d = zeros((nx, ny))
    radius = round(float(gaussian_kernel+1)/2)  # add 1 because the radius includes the center.
    sigma = float(radius)
    mask2d = exp(-(((xx-x_init)**2)/(2*(sigma**2)) + ((yy-y_init)**2)/(2*(sigma**2))))

    # Save mask to 2d file
    file_mask_split = ['tmp.mask_orient_Z'+str(z).zfill(4) for z in range(0, nz, 1)]
    nii_mask2d = Image('tmp.anat_orient_Z0000.nii')
    nii_mask2d.data = mask2d
    nii_mask2d.setFileName(file_mask_split[z_init]+'.nii')
    nii_mask2d.save()

    # initialize variables
    file_mat = ['tmp.mat_Z'+str(z).zfill(4) for z in range(0, nz, 1)]
    file_mat_inv = ['tmp.mat_inv_Z'+str(z).zfill(4) for z in range(0, nz, 1)]
    file_mat_inv_cuemul = ['tmp.mat_inv_cuemul_Z'+str(z).zfill(4) for z in range(0, nz, 1)]

    # create identity matrix for initial transformation matrix
    fid = open(file_mat_inv_cuemul[z_init], 'w')
    fid.write('%i %i %i %i\n' % (1, 0, 0, 0))
    fid.write('%i %i %i %i\n' % (0, 1, 0, 0))
    fid.write('%i %i %i %i\n' % (0, 0, 1, 0))
    fid.write('%i %i %i %i\n' % (0, 0, 0, 1))
    fid.close()

    # initialize centerline: give value corresponding to initial point
    x_centerline = [x_init]
    y_centerline = [y_init]
    z_centerline = [z_init]
    warning_count = 0

    # go up (1), then down (2) in reference to the binary point
    for iUpDown in range(1, 3):

        if iUpDown == 1:
            # z increases
            slice_gap_signed = slice_gap
        elif iUpDown == 2:
            # z decreases
            slice_gap_signed = -slice_gap
            # reverse centerline (because values will be appended at the end)
            x_centerline.reverse()
            y_centerline.reverse()
            z_centerline.reverse()

        # initialization before looping
        z_dest = z_init  # point given by user
        z_src = z_dest + slice_gap_signed

        # continue looping if 0 <= z < nz
        while 0 <= z_src < nz:

            # print current z:
            print 'z='+str(z_src)+':'

            # estimate transformation
            sct.run(fsloutput+'flirt -in '+file_anat_split[z_src]+' -ref '+file_anat_split[z_dest]+' -schedule ' +
                    file_schedule + ' -verbose 0 -omat ' + file_mat[z_src] +
                    ' -cost normcorr -forcescaling -inweight ' + file_mask_split[z_dest] +
                    ' -refweight '+file_mask_split[z_dest])

            # display transfo
            status, output = sct.run('cat '+file_mat[z_src])
            print output

            # check if transformation is bigger than 1.5x slice_gap
            tx = float(output.split()[3])
            ty = float(output.split()[7])
            norm_txy = linalg.norm([tx, ty], ord=2)
            if norm_txy > 1.5*slice_gap:
                print 'WARNING: Transformation is too large --> using previous one.'
                warning_count = warning_count + 1
                # if previous transformation exists, replace current one with previous one
                if os.path.isfile(file_mat[z_dest]):
                    sct.run('cp '+file_mat[z_dest]+' '+file_mat[z_src])

            # estimate inverse transformation matrix
            sct.run('convert_xfm -omat '+file_mat_inv[z_src]+' -inverse '+file_mat[z_src])

            # compute cuemulative transformation
            sct.run('convert_xfm -omat '+file_mat_inv_cuemul[z_src]+' -concat '+file_mat_inv[z_src]+' '+file_mat_inv_cuemul[z_dest])

            # apply inverse cuemulative transformation to initial gaussian mask (to put it in src space)
            sct.run(fsloutput+'flirt -in '+file_mask_split[z_init]+' -ref '+file_mask_split[z_init]+' -applyxfm -init '+file_mat_inv_cuemul[z_src]+' -out '+file_mask_split[z_src])

            # open inverse cuemulative transformation file and generate centerline
            fid = open(file_mat_inv_cuemul[z_src])
            mat = fid.read().split()
            x_centerline.append(x_init + float(mat[3]))
            y_centerline.append(y_init + float(mat[7]))
            z_centerline.append(z_src)
            #z_index = z_index+1

            # define new z_dest (target slice) and new z_src (moving slice)
            z_dest = z_dest + slice_gap_signed
            z_src = z_src + slice_gap_signed


    # Reconstruct centerline
    # ====================================================================================================

    # reverse back centerline (because it's been reversed once, so now all values are in the right order)
    x_centerline.reverse()
    y_centerline.reverse()
    z_centerline.reverse()

    # fit centerline in the Z-X plane using polynomial function
    print '\nFit centerline in the Z-X plane using polynomial function...'
    coeffsx = polyfit(z_centerline, x_centerline, deg=param.deg_poly)
    polyx = poly1d(coeffsx)
    x_centerline_fit = polyval(polyx, z_centerline)
    # calculate RMSE
    rmse = linalg.norm(x_centerline_fit-x_centerline)/sqrt( len(x_centerline) )
    # calculate max absolute error
    max_abs = max(abs(x_centerline_fit-x_centerline))
    print '.. RMSE (in mm): '+str(rmse*px)
    print '.. Maximum absolute error (in mm): '+str(max_abs*px)

    # fit centerline in the Z-Y plane using polynomial function
    print '\nFit centerline in the Z-Y plane using polynomial function...'
    coeffsy = polyfit(z_centerline, y_centerline, deg=param.deg_poly)
    polyy = poly1d(coeffsy)
    y_centerline_fit = polyval(polyy, z_centerline)
    # calculate RMSE
    rmse = linalg.norm(y_centerline_fit-y_centerline)/sqrt( len(y_centerline) )
    # calculate max absolute error
    max_abs = max( abs(y_centerline_fit-y_centerline) )
    print '.. RMSE (in mm): '+str(rmse*py)
    print '.. Maximum absolute error (in mm): '+str(max_abs*py)

    # display
    if param.debug == 1:
        import matplotlib.pyplot as plt
        plt.figure()
        plt.plot(z_centerline,x_centerline,'.',z_centerline,x_centerline_fit,'r')
        plt.legend(['Data','Polynomial Fit'])
        plt.title('Z-X plane polynomial interpolation')
        plt.show()

        plt.figure()
        plt.plot(z_centerline,y_centerline,'.',z_centerline,y_centerline_fit,'r')
        plt.legend(['Data','Polynomial Fit'])
        plt.title('Z-Y plane polynomial interpolation')
        plt.show()

    # generate full range z-values for centerline
    z_centerline_full = [iz for iz in range(0, nz, 1)]

    # calculate X and Y values for the full centerline
    x_centerline_fit_full = polyval(polyx, z_centerline_full)
    y_centerline_fit_full = polyval(polyy, z_centerline_full)

    # Generate fitted transformation matrices and write centerline coordinates in text file
    print '\nGenerate fitted transformation matrices and write centerline coordinates in text file...'
    file_mat_inv_cuemul_fit = ['tmp.mat_inv_cuemul_fit_z'+str(z).zfill(4) for z in range(0,nz,1)]
    file_mat_cuemul_fit = ['tmp.mat_cuemul_fit_z'+str(z).zfill(4) for z in range(0,nz,1)]
    fid_centerline = open('tmp.centerline_coordinates.txt', 'w')
    for iz in range(0, nz, 1):
        # compute inverse cuemulative fitted transformation matrix
        fid = open(file_mat_inv_cuemul_fit[iz], 'w')
        fid.write('%i %i %i %f\n' % (1, 0, 0, x_centerline_fit_full[iz]-x_init))
        fid.write('%i %i %i %f\n' % (0, 1, 0, y_centerline_fit_full[iz]-y_init))
        fid.write('%i %i %i %i\n' % (0, 0, 1, 0))
        fid.write('%i %i %i %i\n' % (0, 0, 0, 1))
        fid.close()
        # compute forward cuemulative fitted transformation matrix
        sct.run('convert_xfm -omat '+file_mat_cuemul_fit[iz]+' -inverse '+file_mat_inv_cuemul_fit[iz])
        # write centerline coordinates in x, y, z format
        fid_centerline.write('%f %f %f\n' %(x_centerline_fit_full[iz], y_centerline_fit_full[iz], z_centerline_full[iz]) )
    fid_centerline.close()


    # Prepare output data
    # ====================================================================================================

    # write centerline as text file
    for iz in range(0, nz, 1):
        # compute inverse cuemulative fitted transformation matrix
        fid = open(file_mat_inv_cuemul_fit[iz], 'w')
        fid.write('%i %i %i %f\n' % (1, 0, 0, x_centerline_fit_full[iz]-x_init))
        fid.write('%i %i %i %f\n' % (0, 1, 0, y_centerline_fit_full[iz]-y_init))
        fid.write('%i %i %i %i\n' % (0, 0, 1, 0))
        fid.write('%i %i %i %i\n' % (0, 0, 0, 1))
        fid.close()

    # write polynomial coefficients
    savetxt('tmp.centerline_polycoeffs_x.txt',coeffsx)
    savetxt('tmp.centerline_polycoeffs_y.txt',coeffsy)

    # apply transformations to data
    print '\nApply fitted transformation matrices...'
    file_anat_split_fit = ['tmp.anat_orient_fit_z'+str(z).zfill(4) for z in range(0, nz, 1)]
    file_mask_split_fit = ['tmp.mask_orient_fit_z'+str(z).zfill(4) for z in range(0, nz, 1)]
    file_point_split_fit = ['tmp.point_orient_fit_z'+str(z).zfill(4) for z in range(0, nz, 1)]
    for iz in range(0, nz, 1):
        # forward cuemulative transformation to data
        sct.run(fsloutput+'flirt -in '+file_anat_split[iz]+' -ref '+file_anat_split[iz]+' -applyxfm -init '+file_mat_cuemul_fit[iz]+' -out '+file_anat_split_fit[iz])
        # inverse cuemulative transformation to mask
        sct.run(fsloutput+'flirt -in '+file_mask_split[z_init]+' -ref '+file_mask_split[z_init]+' -applyxfm -init '+file_mat_inv_cuemul_fit[iz]+' -out '+file_mask_split_fit[iz])
        # inverse cuemulative transformation to point
        sct.run(fsloutput+'flirt -in '+file_point_split[z_init]+' -ref '+file_point_split[z_init]+' -applyxfm -init '+file_mat_inv_cuemul_fit[iz]+' -out '+file_point_split_fit[iz]+' -interp nearestneighbour')

    # Merge into 4D volume
    print '\nMerge into 4D volume...'
    # im_anat_list = [Image(fname) for fname in glob.glob('tmp.anat_orient_fit_z*.nii')]
    fname_anat_list = glob.glob('tmp.anat_orient_fit_z*.nii')
    im_anat_concat = concat_data(fname_anat_list, 2)
    im_anat_concat.setFileName('tmp.anat_orient_fit.nii')
    im_anat_concat.save()

    # im_mask_list = [Image(fname) for fname in glob.glob('tmp.mask_orient_fit_z*.nii')]
    fname_mask_list = glob.glob('tmp.mask_orient_fit_z*.nii')
    im_mask_concat = concat_data(fname_mask_list, 2)
    im_mask_concat.setFileName('tmp.mask_orient_fit.nii')
    im_mask_concat.save()

    # im_point_list = [Image(fname) for fname in 	glob.glob('tmp.point_orient_fit_z*.nii')]
    fname_point_list = glob.glob('tmp.point_orient_fit_z*.nii')
    im_point_concat = concat_data(fname_point_list, 2)
    im_point_concat.setFileName('tmp.point_orient_fit.nii')
    im_point_concat.save()

    # Copy header geometry from input data
    print '\nCopy header geometry from input data...'
    im_anat = Image('tmp.anat_orient.nii')
    im_anat_orient_fit = Image('tmp.anat_orient_fit.nii')
    im_mask_orient_fit = Image('tmp.mask_orient_fit.nii')
    im_point_orient_fit = Image('tmp.point_orient_fit.nii')
    im_anat_orient_fit = copy_header(im_anat, im_anat_orient_fit)
    im_mask_orient_fit = copy_header(im_anat, im_mask_orient_fit)
    im_point_orient_fit = copy_header(im_anat, im_point_orient_fit)
    for im in [im_anat_orient_fit, im_mask_orient_fit, im_point_orient_fit]:
        im.save()

    # Reorient outputs into the initial orientation of the input image
    print '\nReorient the centerline into the initial orientation of the input image...'
    set_orientation('tmp.point_orient_fit.nii', input_image_orientation, 'tmp.point_orient_fit.nii')
    set_orientation('tmp.mask_orient_fit.nii', input_image_orientation, 'tmp.mask_orient_fit.nii')

    # Generate output file (in current folder)
    print '\nGenerate output file (in current folder)...'
    os.chdir('..')  # come back to parent folder
    fname_output_centerline = sct.generate_output_file(path_tmp+'/tmp.point_orient_fit.nii', file_anat+'_centerline'+ext_anat)

    # Delete temporary files
    if remove_tmp_files == 1:
        print '\nRemove temporary files...'
        sct.run('rm -rf '+path_tmp, error_exit='warning')

    # print number of warnings
    print '\nNumber of warnings: '+str(warning_count)+' (if >10, you should probably reduce the gap and/or increase the kernel size'

    # display elapsed time
    elapsed_time = time() - start_time
    print '\nFinished! \n\tGenerated file: '+fname_output_centerline+'\n\tElapsed time: '+str(int(round(elapsed_time)))+'s\n'

Example 30

Project: librosa Source File: spectral.py
def poly_features(y=None, sr=22050, S=None, n_fft=2048, hop_length=512,
                  order=1, freq=None):
    '''Get coefficients of fitting an nth-order polynomial to the columns
    of a spectrogram.

    Parameters
    ----------
    y : np.ndarray [shape=(n,)] or None
        audio time series

    sr : number > 0 [scalar]
        audio sampling rate of `y`

    S : np.ndarray [shape=(d, t)] or None
        (optional) spectrogram magnitude

    n_fft : int > 0 [scalar]
        FFT window size

    hop_length : int > 0 [scalar]
        hop length for STFT. See `librosa.core.stft` for details.

    order : int > 0
        order of the polynomial to fit

    freq : None or np.ndarray [shape=(d,) or shape=(d, t)]
        Center frequencies for spectrogram bins.
        If `None`, then FFT bin center frequencies are used.
        Otherwise, it can be a single array of `d` center frequencies,
        or a matrix of center frequencies as constructed by
        `librosa.core.ifgram`

    Returns
    -------
    coefficients : np.ndarray [shape=(order+1, t)]
        polynomial coefficients for each frame

    Examples
    --------
    >>> y, sr = librosa.load(librosa.util.example_audio_file())
    >>> S = np.abs(librosa.stft(y))

    Line features

    >>> line = librosa.feature.poly_features(S=S, sr=sr)
    >>> line
    array([[ -2.406e-08,  -5.051e-06, ...,  -1.103e-08,  -5.651e-09],
           [  3.445e-04,   3.834e-02, ...,   2.661e-04,   2.239e-04]])

    Quadratic features

    >>> quad = librosa.feature.poly_features(S=S, order=2)
    >>> quad
    array([[  6.276e-12,   2.010e-09, ...,   1.493e-12,   1.000e-13],
           [ -9.325e-08,  -2.721e-05, ...,  -2.749e-08,  -6.754e-09],
           [  4.715e-04,   7.902e-02, ...,   2.963e-04,   2.259e-04]])

    Plot the results for comparison

    >>> import matplotlib.pyplot as plt
    >>> plt.figure()
    >>> plt.subplot(3, 1, 1)
    >>> librosa.display.specshow(line)
    >>> plt.colorbar()
    >>> plt.title('Line coefficients')
    >>> plt.subplot(3, 1, 2)
    >>> librosa.display.specshow(quad)
    >>> plt.colorbar()
    >>> plt.title('Quadratic coefficients')
    >>> plt.subplot(3, 1, 3)
    >>> librosa.display.specshow(librosa.logamplitude(S**2, ref_power=np.max),
    ...                          y_axis='log', x_axis='time')
    >>> plt.title('log Power spectrogram')
    >>> plt.colorbar(format='%+2.0f dB')
    >>> plt.tight_layout()

    '''

    S, n_fft = _spectrogram(y=y, S=S, n_fft=n_fft, hop_length=hop_length)

    # Compute the center frequencies of each bin
    if freq is None:
        freq = fft_frequencies(sr=sr, n_fft=n_fft)

    # If frequencies are constant over frames, then we only need to fit once
    if freq.ndim == 1:
        coefficients = np.polyfit(freq, S, order)
    else:
        # Else, fit each frame independently and stack the results
        coefficients = np.concatenate([[np.polyfit(freq[:, i], S[:, i], order)]
                                       for i in range(S.shape[1])], axis=0).T

    return coefficients

Example 31

Project: pymatgen Source File: eos.py
Function: deltafactor_polyfit
def deltafactor_polyfit(volumes, energies):
    """
    This is the routine used to compute V0, B0, B1 in the deltafactor code.

    Taken from deltafactor/eosfit.py
    """
    fitdata = np.polyfit(volumes**(-2./3.), energies, 3, full=True)
    ssr = fitdata[1]
    sst = np.sum((energies - np.average(energies))**2.)
    residuals0 = ssr/sst
    deriv0 = np.poly1d(fitdata[0])
    deriv1 = np.polyder(deriv0, 1)
    deriv2 = np.polyder(deriv1, 1)
    deriv3 = np.polyder(deriv2, 1)

    v0 = 0
    x = 0
    for x in np.roots(deriv1):
        if x > 0 and deriv2(x) > 0:
            v0 = x**(-3./2.)
            break
    else:
        raise EOSError("No minimum could be found")

    derivV2 = 4./9. * x**5. * deriv2(x)
    derivV3 = (-20./9. * x**(13./2.) * deriv2(x) - 8./27. * x**(15./2.) * deriv3(x))
    b0 = derivV2 / x**(3./2.)
    b1 = -1 - x**(-3./2.) * derivV3 / derivV2

    #print('deltafactor polyfit:')
    #print('e0, b0, b1, v0')
    #print(fitdata[0], b0, b1, v0)

    n = collections.namedtuple("DeltaFitResults", "v0 b0 b1 poly1d")
    return n(v0, b0, b1, fitdata[0])

Example 32

Project: jcvi Source File: synteny.py
Function: get_extent
    def get_extent(self, i, order, debug=True):
        # Some blocks file, such as ones manually edited, will have garbled
        # order, which prompts the hack below
        acol = [order[x][0] for x in self.columns[0] if x in order]
        bcol = [order[x][0] for x in self.columns[i] if x in order]
        elen = min(len(acol), len(bcol))
        ia, ib = acol[:elen], bcol[:elen]
        slope, intercept = np.polyfit(ia, ib, 1)
        orientation = '+' if slope >= 0 else '-'

        ocol = [order[x] for x in self.columns[i] if x in order]
        #orientation = '+' if ocol[0][0] <= ocol[-1][0] else '-'
        si, start = min(ocol)
        ei, end = max(ocol)
        same_chr = (start.seqid == end.seqid)
        chr = start.seqid if same_chr else None
        ngenes = ei - si + 1
        if debug:
            r = "{0}:{1}-{2}".format(chr, start.start, end.end)
            print >> sys.stderr, "Column {0}: {1} - {2} ({3})".\
                    format(i, start.accn, end.accn, r)
            print >> sys.stderr, "  {0} .. {1} ({2}) features .. {3}".\
                    format(chr, ngenes, len(ocol), orientation)

        span = abs(start.start - end.end)

        return start, end, si, ei, chr, orientation, span

Example 33

Project: pymatgen Source File: eos_2.py
Function: initial_guess
    def initial_guess(self):
        """
        Quadratic fit to get an initial guess for the parameters.

        Returns:
            tuple: (minimum energy(e0), buk modulus(b0),
                derivative of bulk modulus wrt pressure(b1),
                minimum volume(v0))
        """
        a, b, c = np.polyfit(self.volumes, self.energies, 2)
        self.eos_params = [a, b, c]

        v0 = -b / (2 * a)
        e0 = a * v0 ** 2 + b * v0 + c
        b0 = 2 * a * v0
        b1 = 4  # b1 is usually a small number like 4

        vmin, vmax = self.volumes.min(), self.volumes.max()

        if not vmin < v0 and v0 < vmax:
            raise EOSError('The minimum volume of a fitted parabola is '
                           'not in the input volumes\n.')

        return e0, b0, b1, v0

Example 34

Project: pandashells Source File: p_regplot.py
def main():
    msg = textwrap.dedent(
        """
        Create a single variable regression plot of specified order.

        -----------------------------------------------------------------------
        Examples:
            * Fit a line to synthetic data with boostrap errors.
                p.linspace 0 10 20 \\
                | p.df 'df["y_true"] = .2 * df.x' \\
                       'df["noise"] = np.random.randn(20)' \\
                        'df["y"] = df.y_true + df.noise' --names x \\
                | p.regplot -x x -y y

            * Fit a quadratic to synthetic data with boostrap errors.
                p.linspace 0 10 40 \\
                | p.df 'df["y_true"] = .5 * df.x  + .3 * df.x ** 2'\\
                       'df["noise"] = np.random.randn(40)' \\
                        'df["y"] = df.y_true + df.noise' --names x \\
                | p.regplot -x x -y y --order 2

            * Fit sealevel data with no bootstrap
                p.example_data -d sealevel\\
                | p.regplot -x year -y sealevel_mm --n_boot 1


        -----------------------------------------------------------------------
        """
    )

    #  read command line arguments
    parser = argparse.ArgumentParser(
        formatter_class=argparse.RawDescriptionHelpFormatter, description=msg)

    arg_lib.add_args(parser, 'io_in', 'io_out', 'decorating')

    msg = 'Column for dependent variable'
    parser.add_argument('-x', nargs=1, type=str, dest='x', metavar='col',
                        help=msg, required=True)

    msg = 'Column for independent variable'
    parser.add_argument('-y', nargs=1, type=str, dest='y',
                        metavar='col', help=msg, required=True)

    msg = 'The order of the polynomial to fit (default = 1)'
    parser.add_argument('--order', help=msg, nargs=1, default=[1], type=int)

    msg = 'Number of bootstrap samples for uncertainty region (default=1000)'
    parser.add_argument(
        '--n_boot', help=msg, nargs=1, default=[1000], type=int)

    parser.add_argument('-a', '--alpha', help='Set opacity',
                        nargs=1, default=[0.5], type=float)

    # parse arguments
    args = parser.parse_args()

    # get the input dataframe
    df = io_lib.df_from_input(args)

    # extract command line params
    x = df[args.x[0]].values
    y = df[args.y[0]].values

    # do a polyfit with the specified order
    coeffs = np.polyfit(x, y, args.order[0])

    label = make_label(coeffs, args.savefig)

    sns.regplot(
        x, y, order=args.order[0], n_boot=args.n_boot[0],
        line_kws={'label': label, 'color': CC[2], 'alpha': .5},
        scatter_kws={'alpha': args.alpha[0], 'color': CC[0]})

    pl.legend(loc='best')
    pl.xlabel(args.x[0])
    pl.ylabel(args.y[0])
    plot_lib.refine_plot(args)
    plot_lib.show(args)

Example 35

Project: odl Source File: operator.py
    def _adjoint_definition(self):
        """Verify ``<Ax, y> == <x, A^* y>``."""
        left_inner_vals = []
        right_inner_vals = []

        with FailCounter(
                test_name='Verifying the identity <Ax, y> = <x, A^T y>',
                err_msg='error = |<Ax, y< - <x, A^* y>| / ||A|| ||x|| ||y||',
                logger=self.log) as counter:

            for [name_x, x], [name_y, y] in samples(self.operator.domain,
                                                    self.operator.range):
                x_norm = x.norm()
                y_norm = y.norm()

                l_inner = self.operator(x).inner(y)
                r_inner = x.inner(self.operator.adjoint(y))

                denom = self.operator_norm * x_norm * y_norm
                error = 0 if denom == 0 else abs(l_inner - r_inner) / denom

                if error > self.tol:
                    counter.fail('x={:25s} y={:25s} : error={:6.5f}'
                                 ''.format(name_x, name_y, error))

                left_inner_vals.append(l_inner)
                right_inner_vals.append(r_inner)

        scale = np.polyfit(left_inner_vals, right_inner_vals, 1)[0]
        self.log('\nThe adjoint seems to be scaled according to:')
        self.log('(x, A^T y) / (Ax, y) = {}. Should be 1.0'.format(scale))

Example 36

Project: FreeCAD_assembly2 Source File: lineSearches.py
def quadraticLineSearch( f, x1, f1, intialStep, it, debugPrintLevel, printF, tol_stag=3, tol_x=10**-6): 
    if norm(intialStep) == 0:
        printF('    quadraticLineSearch: norm search direction is 0, aborting!')
        return x1
    def LSEval(lam, fv=None):
        return LineSearchEvaluation( f,  x1, intialStep, lam, fv )
    Y = [ LSEval(     0.0,  f1), LSEval(       1 ), LSEval(       2 )]
    y_min_prev = min(Y)
    count_stagnation = 0
    tol_lambda = tol_x / norm(intialStep)
    for k in range(it):
        Y.sort()
        if debugPrintLevel > 0:
            printF('    quadratic line search   it %i, fmin %1.2e, lam  %1.6f  %1.6f  %1.6f, f(lam) %1.2e  %1.2e  %1.2e'%( k+1, Y[0].fv, Y[0].lam,Y[1].lam,Y[2].lam,Y[0].fv,Y[1].fv,Y[2].fv ))
        #``p[0]*x**(N-1) + p[1]*x**(N-2) + ... + p[N-2]*x + p[N-1]``
        quadraticCoefs, residuals, rank, singular_values, rcond = numpy.polyfit( [y.lam for y in Y], [y.fv for y in Y], 2, full=True)
        if quadraticCoefs[0] > 0 and rank == 3:
            lam_c = -quadraticCoefs[1] / (2*quadraticCoefs[0]) #diff poly a*x**2 + b*x + c -> grad_poly = 2*a*x + b
            lam_c = min( max( [y.lam for y in Y])*4, lam_c)
            if lam_c < 0:
                if debugPrintLevel > 1:  printF('    quadratic line search lam_c < 0')
                lam_c = 1.0 / (k + 1) ** 2
        else:
            if debugPrintLevel > 1: printF('    quadratic fit invalid, using interval halving instead')
            lam_c = ( Y[0].lam + Y[1].lam )/2
        del Y[2] # Y sorted at start of each iteration
        Y.append(  LSEval( lam_c ))
        y_min = min(Y)
        if y_min == y_min_prev:
            count_stagnation = count_stagnation + 1
            if count_stagnation > tol_stag:
                if debugPrintLevel > 0:  printF('    terminating quadratic line search as count_stagnation > tol_stag')
                break
        else:
            y_min_prev = y_min
            count_stagnation = 0
        Lam = [y.lam for y in Y] 
        if max(Lam) - min(Lam) < tol_lambda:
            if debugPrintLevel > 0:  printF('    terminating quadratic max(Lam)-min(Lam) < tol_lambda (%e < %e)' % (max(Lam) - min(Lam), tol_lambda))
            break

    return min(Y).xv

Example 37

Project: chimera Source File: autofocus.py
Function: fit
    @staticmethod
    def fit(position, fwhm, temperature=None, minmax=None):

        if minmax and len(minmax) >= 2:
            idxs = (fwhm >= minmax[0]) & (fwhm <= minmax[1])
            position = position[idxs]
            fwhm = fwhm[idxs]

        A, B, C = N.polyfit(position, fwhm, 2)

        fwhm_fit = N.polyval([A, B, C], position)

        err = sqrt(sum((fwhm_fit - fwhm)**2) / len(position))

        fit = FocusFit()
        fit.position = position
        fit.fwhm = fwhm
        fit.temperature = temperature
        fit.minmax = minmax

        fit.A, fit.B, fit.C = A, B, C
        fit.err = err
        fit.fwhm_fit = fwhm_fit

        return fit

Example 38

Project: robothon Source File: test_regression.py
Function: test_polyfit_build
    def test_polyfit_build(self,level=rlevel):
        """Ticket #628"""
        ref = [-1.06123820e-06, 5.70886914e-04, -1.13822012e-01,
                9.95368241e+00, -3.14526520e+02]
        x = [90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
             104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115,
             116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 129,
             130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
             146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157,
             158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
             170, 171, 172, 173, 174, 175, 176]
        y = [9.0, 3.0, 7.0, 4.0, 4.0, 8.0, 6.0, 11.0, 9.0, 8.0, 11.0, 5.0,
             6.0, 5.0, 9.0, 8.0, 6.0, 10.0, 6.0, 10.0, 7.0, 6.0, 6.0, 6.0,
             13.0, 4.0, 9.0, 11.0, 4.0, 5.0, 8.0, 5.0, 7.0, 7.0, 6.0, 12.0,
             7.0, 7.0, 9.0, 4.0, 12.0, 6.0, 6.0, 4.0, 3.0, 9.0, 8.0, 8.0,
             6.0, 7.0, 9.0, 10.0, 6.0, 8.0, 4.0, 7.0, 7.0, 10.0, 8.0, 8.0,
             6.0, 3.0, 8.0, 4.0, 5.0, 7.0, 8.0, 6.0, 6.0, 4.0, 12.0, 9.0,
             8.0, 8.0, 8.0, 6.0, 7.0, 4.0, 4.0, 5.0, 7.0]
        tested = np.polyfit(x, y, 4)
        assert_array_almost_equal(ref, tested)

Example 39

Project: xraylarch Source File: fitpeak.py
Function: guess_starting_values
    def guess_starting_values(self, y, x):
        sval, oval = np.polyfit(x, y, 1)
        self.set_initval(self.params.offset, oval)
        self.set_initval(self.params.slope, sval)

Example 40

Project: hyperspy Source File: polynomial.py
Function: estimate_parameters
    def estimate_parameters(self, signal, x1, x2, only_current=False):
        """Estimate the parameters by the two area method

        Parameters
        ----------
        signal : Signal1D instance
        x1 : float
            Defines the left limit of the spectral range to use for the
            estimation.
        x2 : float
            Defines the right limit of the spectral range to use for the
            estimation.

        only_current : bool
            If False estimates the parameters for the full dataset.

        Returns
        -------
        bool

        """
        super(Polynomial, self)._estimate_parameters(signal)
        axis = signal.axes_manager.signal_axes[0]
        binned = signal.metadata.Signal.binned
        i1, i2 = axis.value_range_to_indices(x1, x2)
        if only_current is True:
            estimation = np.polyfit(axis.axis[i1:i2],
                                    signal()[i1:i2],
                                    self.get_polynomial_order())
            if binned is True:
                self.coefficients.value = estimation / axis.scale
            else:
                self.coefficients.value = estimation
            return True
        else:
            if self.coefficients.map is None:
                self._create_arrays()
            nav_shape = signal.axes_manager._navigation_shape_in_array
            with signal.unfolded():
                dc = signal.data
                # For polyfit the spectrum goes in the first axis
                if axis.index_in_array > 0:
                    dc = dc.T             # Unfolded, so simply transpose
                cmaps = np.polyfit(axis.axis[i1:i2], dc[i1:i2, :],
                                   self.get_polynomial_order())
                if axis.index_in_array > 0:
                    cmaps = cmaps.T       # Transpose back if needed
                # Shape needed to fit coefficients.map:
                cmap_shape = nav_shape + (self.get_polynomial_order() + 1, )
                self.coefficients.map['values'][:] = cmaps.reshape(cmap_shape)
                if binned is True:
                    self.coefficients.map["values"] /= axis.scale
                self.coefficients.map['is_set'][:] = True
            self.fetch_stored_values()
            return True

Example 41

Project: emailinsight Source File: kerasPlottingUtils.py
    def plot_regression_line(self, plot_ax, data, batches, future_batches,
                             n_backwards, linestyle, label):
        """Calculates and plots a regression line.
        Args:
            plot_ax: The ax on which to plot the line.
            data: The data used to perform the regression
                (e.g. training loss values).
            batchs: List of all batchs (0, 1, 2, ...).
            future_batchs: List of the future batchs for which values are
                ought to be predicted.
            n_backwards: How far back to go in time (in batchs) in order
                to compute the regression. (E.g. 10 = calculate it on the
                last 10 values max.)
            linestyle: Linestyle of the regression line.
            label: Label of the regression line.
        """
        # dont try to draw anything if the data list is empty or it's the
        # first batch
        if len(data) > 1:
            poly = np.poly1d(np.polyfit(batches[-n_backwards:],
                                        data[-n_backwards:], self.poly_degree))
            future_values = [poly(i) for i in future_batches]
            plot_ax.plot(future_batches, future_values, linestyle, label=label)

Example 42

Project: xraylarch Source File: fitpeak.py
Function: guess_starting_values
    def guess_starting_values(self, y, x):
        qval, sval, oval = np.polyfit(x, y, 2)
        self.set_initval(self.params.offset, oval)
        self.set_initval(self.params.slope, sval)
        self.set_initval(self.params.quad, qval)

Example 43

Project: landlab Source File: fault_facet_finder.py
    def analyse_fault_trace(self, fault_trace_node_ids):
        """
        This method takes the grid and an array listing the (contiguous) node
        ids of cells that contain a single fault segment trace of interest.

        It sets and returns the azimuth of the fault trace, az,
        -pi/2 < az <= pi/2.
        (i.e., smallest angle between north and the trace).
        """
        self.ft_trace_node_ids = fault_trace_node_ids
        x = self.grid.node_x[fault_trace_node_ids]
        y = self.grid.node_y[fault_trace_node_ids]
        (grad, offset) = np.polyfit(x, y, 1)
        angle = np.arctan(grad)
        if grad >= 0.:
            self.az = np.pi / 2. - angle
        else:
            self.az = -np.pi / 2. - angle

        return self.az

Example 44

Project: odl Source File: operator.py
    def self_adjoint(self):
        """Verify ``<Ax, y> == <x, Ay>``."""
        left_inner_vals = []
        right_inner_vals = []

        with FailCounter(
                test_name='Verifying the identity <Ax, y> = <x, Ay>',
                err_msg='error = |<Ax, y> - <x, Ay>| / ||A|| ||x|| ||y||',
                logger=self.log) as counter:

            for [name_x, x], [name_y, y] in samples(self.operator.domain,
                                                    self.operator.range):
                x_norm = x.norm()
                y_norm = y.norm()

                l_inner = self.operator(x).inner(y)
                r_inner = x.inner(self.operator(y))

                denom = self.operator_norm * x_norm * y_norm
                error = 0 if denom == 0 else abs(l_inner - r_inner) / denom

                if error > self.tol:
                    counter.fail('x={:25s} y={:25s} : error={:6.5f}'
                                 ''.format(name_x, name_y, error))

                left_inner_vals.append(l_inner)
                right_inner_vals.append(r_inner)

        scale = np.polyfit(left_inner_vals, right_inner_vals, 1)[0]
        self.log('\nThe adjoint seems to be scaled according to:')
        self.log('(x, Ay) / (Ax, y) = {}. Should be 1.0'.format(scale))

Example 45

Project: compare-codecs Source File: visual_metrics.py
Function: bdrate
def bdrate(metric_set1, metric_set2):
  """
  BJONTEGAARD    Bjontegaard metric calculation
  Bjontegaard's metric allows to compute the average % saving in bitrate
  between two rate-distortion curves [1].

  rate1,psnr1 - RD points for curve 1
  rate2,psnr2 - RD points for curve 2

  adapted from code from: (c) 2010 Giuseppe Valenzise

  """
  # numpy plays games with its exported functions.
  # pylint: disable=no-member
  # pylint: disable=too-many-locals
  # pylint: disable=bad-builtin
  rate1 = [x[0] for x in metric_set1]
  psnr1 = [x[1] for x in metric_set1]
  rate2 = [x[0] for x in metric_set2]
  psnr2 = [x[1] for x in metric_set2]

  log_rate1 = map(math.log, rate1)
  log_rate2 = map(math.log, rate2)

  # Best cubic poly fit for graph represented by log_ratex, psrn_x.
  poly1 = numpy.polyfit(psnr1, log_rate1, 3)
  poly2 = numpy.polyfit(psnr2, log_rate2, 3)

  # Integration interval.
  min_int = max([min(psnr1), min(psnr2)])
  max_int = min([max(psnr1), max(psnr2)])

  # find integral
  p_int1 = numpy.polyint(poly1)
  p_int2 = numpy.polyint(poly2)

  # Calculate the integrated value over the interval we care about.
  int1 = numpy.polyval(p_int1, max_int) - numpy.polyval(p_int1, min_int)
  int2 = numpy.polyval(p_int2, max_int) - numpy.polyval(p_int2, min_int)

  # Calculate the average improvement.
  avg_exp_diff = (int2 - int1) / (max_int - min_int)

  # In really bad formed data the exponent can grow too large.
  # clamp it.
  if avg_exp_diff > 200:
    avg_exp_diff = 200

  # Convert to a percentage.
  avg_diff = (math.exp(avg_exp_diff) - 1) * 100

  return avg_diff

Example 46

Project: histwords Source File: plothelper.py
def trendline(xd, yd, order=1, c='r', alpha=1, plot_r=False, text_pos=None):
    """Make a line of best fit"""

    #Calculate trendline
    coeffs = np.polyfit(xd, yd, order)

    intercept = coeffs[-1]
    slope = coeffs[-2]
    if order == 2: power = coeffs[0]
    else: power = 0

    minxd = np.min(xd)
    maxxd = np.max(xd)

    xl = np.array([minxd, maxxd])
    yl = power * xl ** 2 + slope * xl + intercept

    #Plot trendline
    plt.plot(xl, yl, color=c, alpha=alpha)

    #Calculate R Squared
    r = sp.stats.pearsonr(xd, yd)[0]

    if plot_r == False:
        #Plot R^2 value
        if text_pos == None:
            text_pos = (0.9 * maxxd + 0.1 * minxd, 0.9 * np.max(yd) + 0.1 * np.min(yd),)
        plt.text(text_pos[0], text_pos[1], '$R = %0.2f$' % r)
    else:
        #Return the R^2 value:
        return r

Example 47

Project: homer Source File: fingerprint.py
def fingerprint_page(page, grid_size=16):
    "Extract scale-invariant fingerprint from the systems on the page"
    # Find the average angle of each system to undo rotation.
    # This may be a useful step to undo rotation on the actual image.
    # Also, fit a quadratic for deskewing?
    angles = []
    for sys in page.systems:
        xs = sys['barlines'][:, :, 0].mean(1)
        ys = sys['barlines'][:, :, 1].mean(1)
        coefs = np.polyfit(xs, ys, 1)
        theta = np.arctan(coefs[0])
        angles.append(theta)
    theta = np.mean(angles)
    # Rotate all x values by -theta
    rotation = np.array([[np.cos(-theta), -np.sin(-theta)],
                         [np.sin(-theta),  np.cos(-theta)]])
    systems_xs = []
    for sys in page.systems:
        def barline_rotated_x(barline):
            center = barline.mean(0)
            rotated_center = rotation.dot(center[:, None])
            return rotated_center[0, 0]
        systems_xs.append(map(barline_rotated_x, sys['barlines']))
    # Make margin-invariant
    xmin = min(min(xs) for xs in systems_xs)
    xmax = max(max(xs) for xs in systems_xs)
    width = xmax - xmin
    systems_xs = [[x - xmin for x in xs] for xs in systems_xs]
    systems_str = StringIO()
    for xs in systems_xs:
        grid_pos = ((np.asarray(xs) * grid_size) / width).astype(int)
        grid_pos[grid_pos >= grid_size] = grid_size - 1
        systems_str.write(','.join(map(str, grid_pos)))
        systems_str.write('\n')
    systems_str.seek(0)
    systems_str = systems_str.read()
    sha = hashlib.sha1(systems_str).digest()
    fingerprint, = struct.unpack("<L", sha[:4])
    return fingerprint

Example 48

Project: memopol2 Source File: views.py
def mep(request, mep_id):
    mep_ = MEP.get(mep_id)
    mep_['achievements']=autoTrophies(mep_)
    positions = Position.objects.filter(mep_id=mep_id)
    score_list = mep_.scores
    for score in score_list:
        score['color'] = score_to_color(int(score['value']))
    score_list.sort(key = lambda k : datetime.strptime(k['date'], "%d/%m/%Y"))
    scores = [s['value'] for s in mep_.scores]

    if score_list:
        try:
            import numpy
            import matplotlib
            matplotlib.use("Agg")
            from matplotlib import pyplot

            pyplot.plot(scores, 'bo')
            a, b = numpy.polyfit(range(len(scores)), [int(x) for x in scores], 1)
            pyplot.plot([a*int(x) + b for x in range(len(scores))])
            pyplot.legend(('Scores', 'Mediane'), 'best', shadow=True)
            pyplot.plot(scores)
            pyplot.axis([0, len(scores) - 1, 0, 102])
            pyplot.title("%s - Votes notes evolution over time" % (mep_.infos['name']['full']))
            pyplot.xticks(range(len(scores)), [k['date'] for k in score_list])
            pyplot.xlabel("Votes dates")
            pyplot.ylabel("Scores on votes")
            pyplot.savefig(realpath("./memopol2/%simg/trends/meps/%s-scores.png" % (settings.MEDIA_URL, mep_id)), format="png")
            pyplot.clf()
        except ImportError:
            pass

    context = {
        'mep_id': mep_id,
        'mep': mep_,
        'positions': positions,
        'visible_count': len([x for x in positions if x.visible]),
        'average': sum(scores)/len(scores) if len(scores) > 0 else "",
        'score_list' : score_list,
    }
    return direct_to_template(request, 'meps/mep.html', context)

Example 49

Project: xraylarch Source File: xrd_calc.py
Function: fit_background
def fit_background(q,I):
    
    ## Working on background calculation
    ## mkak 2016.09.28
    
    x = q
    y = I
    pfit = np.polyfit(x,y,4)
    yfit = np.polyval(pfit,x)
    #panel.plot(xrd_spectra[0], xrd_spectra[1]-yfit, label='no bkg')
    #panel.plot(xrd_spectra[0], yfit, color='blue', label='bkg')
    
    ### calculation works, but plotting here wipes previous plots - only shows last
    import matplotlib as plt
    plt.figure()
    plt.plot(x,y,label='raw data')
    plt.plot(x,yfit,label='background')
    plt.plot(x,y-yfit,label='background subtracted')
    plt.legend()
    plt.show()

Example 50

Project: landlab Source File: fault_facet_finder.py
    def fit_slopes_to_facet_lines(self, polynomial_degree=4, curvature_threshold=0.0004):
        """
        Fits (linear) lines of best fit to extracted profiles, already stored as
        class properties.
        """
        avg_slopes_linear = []
        avg_slopes_poly = []
        curv_of_flattest_part_list = []
        slope_min_curv = []
        rsqd_list = []
        big_slope_small_curv = []
        elev_at_bssc = []
        for i in range(len(self.profile_x_facet_pts)):
            x = self.profile_x_facet_pts[i]
            z = self.profile_z_facet_pts[i]
            (grad, offset) = np.polyfit(x, z, 1)
            coeffs, residuals = np.polyfit(
                x, z, polynomial_degree, full=True)[:2]
            rsqd = 1. - residuals / (z.size * z.var())
            # differentiate the coeffs to get slope:
            diff_multiplier = np.arange(polynomial_degree + 1)[::-1]
            curv_multiplier = np.arange(polynomial_degree)[::-1]
            z_equ = np.poly1d(coeffs)
            S_equ = np.poly1d((coeffs * diff_multiplier)[:-1])
            curv_equ = np.poly1d(
                ((coeffs * diff_multiplier)[:-1] * curv_multiplier)[:-1])
            S_at_each_pt = S_equ(x)
            curv_at_each_pt = curv_equ(x)
            avg_slopes_linear.append(abs(grad))
            avg_slopes_poly.append(np.amax(np.fabs(S_at_each_pt)))
            loc_of_flattest_part = np.argmin(
                np.fabs(curv_at_each_pt[2:-2])) + 2
            curv_of_flattest_part = curv_at_each_pt[loc_of_flattest_part]
            S_at_min_curve_untested = abs(S_at_each_pt[loc_of_flattest_part])
            small_curves = np.less(
                np.fabs(curv_at_each_pt[2:-2]), curvature_threshold)
            try:
                big_slope_small_curv.append(
                    np.amax(S_at_each_pt[small_curves]))
                elev_at_bssc.append(z[np.argmax(S_at_each_pt[small_curves])])
            except ValueError:
                big_slope_small_curv.append(np.nan)
                elev_at_bssc.append(np.nan)
            slope_min_curv.append(S_at_min_curve_untested)
            curv_of_flattest_part_list.append(curv_of_flattest_part)
            rsqd_list.append(rsqd)
            # figure(8)
            #synthetic_z = grad*x + offset
            synthetic_z = z_equ(x)
            plot(x, z, 'x')
            plot(x, synthetic_z, '-')
        self.avg_slopes_linear = np.array(avg_slopes_linear)
        self.avg_slopes_poly = np.array(avg_slopes_poly)
        self.curv_of_flattest_part = np.array(curv_of_flattest_part_list)
        self.slope_min_curv = np.array(slope_min_curv)
        self.big_slope_small_curv = np.array(big_slope_small_curv)
        self.elev_at_bssc = np.array(elev_at_bssc)
        self.rsqd = np.array(rsqd_list)