Python numpy.linspace() Examples

The following are 30 code examples for showing how to use numpy.linspace(). These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example.

You may check out the related API usage on the sidebar.

You may also want to check out all available functions/classes of the module numpy , or try the search function .

Example 1
Project: xrft   Author: xgcm   File: xrft.py    License: MIT License 6 votes vote down vote up
def _radial_wvnum(k, l, N, nfactor):
    """ Creates a radial wavenumber based on two horizontal wavenumbers
    along with the appropriate index map
    """

    # compute target wavenumbers
    k = k.values
    l = l.values
    K = np.sqrt(k[np.newaxis,:]**2 + l[:,np.newaxis]**2)
    nbins = int(N/nfactor)
    if k.max() > l.max():
        ki = np.linspace(0., l.max(), nbins)
    else:
        ki = np.linspace(0., k.max(), nbins)

    # compute bin index
    kidx = np.digitize(np.ravel(K), ki)
    # compute number of points for each wavenumber
    area = np.bincount(kidx)
    # compute the average radial wavenumber for each bin
    kr = (np.bincount(kidx, weights=K.ravel())
          / np.ma.masked_where(area==0, area))

    return ki, kr[1:-1] 
Example 2
Project: xrft   Author: xgcm   File: test_xrft.py    License: MIT License 6 votes vote down vote up
def test_cross_phase_2d(self, dask):
        Ny, Nx = (32, 16)
        x = np.linspace(0, 1, num=Nx, endpoint=False)
        y = np.ones(Ny)
        f = 6
        phase_offset = np.pi/2
        signal1 = np.cos(2*np.pi*f*x)  # frequency = 1/(2*pi)
        signal2 = np.cos(2*np.pi*f*x - phase_offset)
        da1 = xr.DataArray(data=signal1*y[:,np.newaxis], name='a',
                          dims=['y','x'], coords={'y':y, 'x':x})
        da2 = xr.DataArray(data=signal2*y[:,np.newaxis], name='b',
                          dims=['y','x'], coords={'y':y, 'x':x})
        with pytest.raises(ValueError):
            xrft.cross_phase(da1, da2, dim=['y','x'])

        if dask:
            da1 = da1.chunk({'x': 16})
            da2 = da2.chunk({'x': 16})
        cp = xrft.cross_phase(da1, da2, dim=['x'])
        actual_phase_offset = cp.sel(freq_x=f).values
        npt.assert_almost_equal(actual_phase_offset, phase_offset) 
Example 3
Project: FRIDA   Author: LCAV   File: doa.py    License: MIT License 6 votes vote down vote up
def compute_mode(self):
        """
        Pre-compute mode vectors from candidate locations (in spherical 
        coordinates).
        """
        if self.num_loc is None:
            raise ValueError('Lookup table appears to be empty. \
                Run build_lookup().')
        self.mode_vec = np.zeros((self.max_bin,self.M,self.num_loc), 
            dtype='complex64')
        if (self.nfft % 2 == 1):
            raise ValueError('Signal length must be even.')
        f = 1.0 / self.nfft * np.linspace(0, self.nfft / 2, self.max_bin) \
            * 1j * 2 * np.pi
        for i in range(self.num_loc):
            p_s = self.loc[:, i]
            for m in range(self.M):
                p_m = self.L[:, m]
                if (self.mode == 'near'):
                    dist = np.linalg.norm(p_m - p_s, axis=1)
                if (self.mode == 'far'):
                    dist = np.dot(p_s, p_m)
                # tau = np.round(self.fs*dist/self.c) # discrete - jagged
                tau = self.fs * dist / self.c  # "continuous" - smoother
                self.mode_vec[:, m, i] = np.exp(f * tau) 
Example 4
Project: MPContribs   Author: materialsproject   File: pre_submission.py    License: MIT License 6 votes vote down vote up
def load_RSM(filename):
    om, tt, psd = xu.io.getxrdml_map(filename)
    om = np.deg2rad(om)
    tt = np.deg2rad(tt)
    wavelength = 1.54056

    q_y = (1 / wavelength) * (np.cos(tt) - np.cos(2 * om - tt))
    q_x = (1 / wavelength) * (np.sin(tt) - np.sin(2 * om - tt))

    xi = np.linspace(np.min(q_x), np.max(q_x), 100)
    yi = np.linspace(np.min(q_y), np.max(q_y), 100)
    psd[psd < 1] = 1
    data_grid = griddata(
        (q_x, q_y), psd, (xi[None, :], yi[:, None]), fill_value=1, method="cubic"
    )
    nx, ny = data_grid.shape

    range_values = [np.min(q_x), np.max(q_x), np.min(q_y), np.max(q_y)]
    output_data = (
        Panel(np.log(data_grid).reshape(nx, ny, 1), minor_axis=["RSM"])
        .transpose(2, 0, 1)
        .to_frame()
    )

    return range_values, output_data 
Example 5
Project: dustmaps   Author: gregreen   File: test_bayestar.py    License: GNU General Public License v2.0 6 votes vote down vote up
def atest_plot_samples(self):
        dm = np.linspace(4., 19., 1001)
        samples = []

        for dm_k in dm:
            d = 10.**(dm_k/5.-2.)
            samples.append(self._interp_ebv(self._test_data[0], d))

        samples = np.array(samples).T
        # print samples

        import matplotlib.pyplot as plt
        fig = plt.figure()
        ax = fig.add_subplot(1,1,1)
        for s in samples:
            ax.plot(dm, s, lw=2., alpha=0.5)

        plt.show() 
Example 6
Project: StructEngPy   Author: zhuoju36   File: spectrum.py    License: MIT License 6 votes vote down vote up
def __init__(self,alpha_max,Tg,xi):
        gamma=0.9+(0.05-xi)/(0.3+6*xi)
        eta1=0.02+(0.05-xi)/(4+32*xi)
        eta1=eta1 if eta1>0 else 0
        eta2=1+(0.05-xi)/(0.08+1.6*xi)
        eta2=eta2 if eta2>0.55 else 0.55
        T=np.linspace(0,6,601)
        alpha=[]
        for t in T:
            if t<0.1:
                alpha.append(np.interp(t,[0,0.1],[0.45*alpha_max,eta2*alpha_max]))
            elif t<Tg:
                alpha.append(eta2*alpha_max)
            elif t<5*Tg:
                alpha.append((Tg/t)**gamma*eta2*alpha_max)
            else:
                alpha.append((eta2*0.2**gamma-eta1*(t-5*Tg))*alpha_max)
        self.__spectrum={'T':T,'alpha':alpha} 
Example 7
Project: neuropythy   Author: noahbenson   File: core.py    License: GNU Affero General Public License v3.0 6 votes vote down vote up
def curve_length(self, start=None, end=None, precision=0.01):
        '''
        Calculates the length of the curve by dividing the curve up
        into pieces of parameterized-length <precision>.
        '''
        if start is None: start = self.t[0]
        if end is None: end = self.t[-1]
        from scipy import interpolate
        if self.order == 1:
            # we just want to add up along the steps...
            ii = [ii for (ii,t) in enumerate(self.t) if start < t and t < end]
            ts = np.concatenate([[start], self.t[ii], [end]])
            xy = np.vstack([[self(start)], self.coordinates[:,ii].T, [self(end)]])
            return np.sum(np.sqrt(np.sum((xy[1:] - xy[:-1])**2, axis=1)))
        else:
            t = np.linspace(start, end, int(np.ceil((end-start)/precision)))
            dt = t[1] - t[0]
            dx = interpolate.splev(t, self.splrep[0], der=1)
            dy = interpolate.splev(t, self.splrep[1], der=1)
            return np.sum(np.sqrt(dx**2 + dy**2)) * dt 
Example 8
Project: TOPFARM   Author: DTUWindEnergy   File: test_topfarm.py    License: GNU Affero General Public License v3.0 6 votes vote down vote up
def test_2x3(self):
        # Loading the water depth map
        dat = loadtxt('data/WaterDepth1.dat')
        X, Y = meshgrid(linspace(0., 1000., 50), linspace(0., 1000., 50))
        depth = array(zip(X.flatten(), Y.flatten(), dat.flatten()))
        borders = array([[200, 200], [150, 500], [200, 800], [600, 900], [700, 700], [900, 500], [800, 200], [500, 100], [200, 200]])
        baseline = array([[587.5, 223.07692308], [525., 346.15384615], [837.5, 530.76923077], [525., 530.76923077], [525., 838.46153846], [837.5, 469.23076923]])

        wt_desc = WTDescFromWTG('data/V80-2MW-offshore.wtg').wt_desc
        wt_layout = GenericWindFarmTurbineLayout([WTPC(wt_desc=wt_desc, position=pos) for pos in baseline])

        t = Topfarm(
            baseline_layout = wt_layout,
            borders = borders,
            depth_map = depth,
            dist_WT_D = 5.0,
            distribution='spiral',
            wind_speeds=[4., 8., 20.],
            wind_directions=linspace(0., 360., 36)[:-1]
        )

        t.run()

        self.fail('make save function')
        t.save() 
Example 9
Project: pywr   Author: pywr   File: thames.py    License: GNU General Public License v3.0 6 votes vote down vote up
def figures(ext, show):

    for name, df in TablesRecorder.generate_dataframes('thames_output.h5'):
        df.columns = ['Very low', 'Low', 'Central', 'High', 'Very high']

        fig, (ax1, ax2) = plt.subplots(figsize=(12, 4), ncols=2, sharey='row',
                                       gridspec_kw={'width_ratios': [3, 1]})
        df['2100':'2125'].plot(ax=ax1)
        df.quantile(np.linspace(0, 1)).plot(ax=ax2)

        if name.startswith('reservoir'):
            ax1.set_ylabel('Volume [$Mm^3$]')
        else:
            ax1.set_ylabel('Flow [$Mm^3/day$]')

        for ax in (ax1, ax2):
            ax.set_title(name)
            ax.grid(True)
        plt.tight_layout()

        if ext is not None:
            fig.savefig(f'{name}.{ext}', dpi=300)

    if show:
        plt.show() 
Example 10
Project: pywr   Author: pywr   File: figures.py    License: GNU General Public License v3.0 6 votes vote down vote up
def plot_percentiles(A, B, ax=None):
    if ax is None:
        ax = plt.gca()
    percentiles = np.linspace(0.001, 0.999, 1000) * 100
    A_pct = scipy.stats.scoreatpercentile(A.values, percentiles)
    B_pct = scipy.stats.scoreatpercentile(B.values, percentiles)
    percentiles = percentiles / 100.0
    ax.plot(percentiles, B_pct[::-1], color=c["Bfill"], clip_on=False, linewidth=2)
    ax.plot(percentiles, A_pct[::-1], color=c["Afill"], clip_on=False, linewidth=2)
    ax.set_xlabel("Cumulative frequency")
    ax.grid(True)
    ax.xaxis.grid(True, which="both")
    set_000formatter(ax.get_yaxis())
    ax.set_xscale("logit")
    xticks = ax.get_xticks()
    xticks_minr = ax.get_xticks(minor=True)
    ax.set_xticklabels([], minor=True)
    ax.set_xticks([0.01, 0.1, 0.5, 0.9, 0.99])
    ax.set_xticklabels(["1", "10", "50", "90", "99"])
    ax.set_xlim(0.001, 0.999)
    ax.legend([B.name, A.name], loc="best")
    return ax 
Example 11
Project: easy-faster-rcnn.pytorch   Author: potterhsu   File: region_proposal_network.py    License: MIT License 6 votes vote down vote up
def generate_anchors(self, image_width: int, image_height: int, num_x_anchors: int, num_y_anchors: int) -> Tensor:
        center_ys = np.linspace(start=0, stop=image_height, num=num_y_anchors + 2)[1:-1]
        center_xs = np.linspace(start=0, stop=image_width, num=num_x_anchors + 2)[1:-1]
        ratios = np.array(self._anchor_ratios)
        ratios = ratios[:, 0] / ratios[:, 1]
        sizes = np.array(self._anchor_sizes)

        # NOTE: it's important to let `center_ys` be the major index (i.e., move horizontally and then vertically) for consistency with 2D convolution
        # giving the string 'ij' returns a meshgrid with matrix indexing, i.e., with shape (#center_ys, #center_xs, #ratios)
        center_ys, center_xs, ratios, sizes = np.meshgrid(center_ys, center_xs, ratios, sizes, indexing='ij')

        center_ys = center_ys.reshape(-1)
        center_xs = center_xs.reshape(-1)
        ratios = ratios.reshape(-1)
        sizes = sizes.reshape(-1)

        widths = sizes * np.sqrt(1 / ratios)
        heights = sizes * np.sqrt(ratios)

        center_based_anchor_bboxes = np.stack((center_xs, center_ys, widths, heights), axis=1)
        center_based_anchor_bboxes = torch.from_numpy(center_based_anchor_bboxes).float()
        anchor_bboxes = BBox.from_center_base(center_based_anchor_bboxes)

        return anchor_bboxes 
Example 12
Project: EXOSIMS   Author: dsavransky   File: test_DulzPlavchan.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_dist_albedo(self):
        """
        Test that albedos outside of the range have zero probability

        """

        spec = copy.deepcopy(self.spec)
        spec['modules']['PlanetPhysicalModel'] = 'FortneyMarleyCahoyMix1'
        with RedirectStreams(stdout=self.dev_null):
            pp = DulzPlavchan(**spec)

        p = np.linspace(pp.prange[0]-1,pp.prange[1]+1,100)

        fp = pp.dist_albedo(p)
        self.assertTrue(np.all(fp[p < pp.prange[0]] == 0),'dist_albedo high bound failed for DulzPlavchan')
        self.assertTrue(np.all(fp[p > pp.prange[1]] == 0),'dist_albedo low bound failed for DulzPlavchan')
        self.assertTrue(np.all(fp[(p >= pp.prange[0]) & (p <= pp.prange[1])] > 0),'dist_albedo generates zero probabilities within range for DulzPlavchan') 
Example 13
Project: EXOSIMS   Author: dsavransky   File: test_PlanetPopulation.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_dist_eccen(self):
        """
        Test that eccentricities outside of the range have zero probability

        """
        for mod in self.allmods:
            if 'dist_eccen' in mod.__dict__:
                with RedirectStreams(stdout=self.dev_null):
                    pp = mod(**self.spec)

                e = np.linspace(pp.erange[0]-1,pp.erange[1]+1,100)

                fe = pp.dist_eccen(e)
                self.assertTrue(np.all(fe[e < pp.erange[0]] == 0),'dist_eccen high bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fe[e > pp.erange[1]] == 0),'dist_eccen low bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fe[(e >= pp.erange[0]) & (e <= pp.erange[1])] > 0),'dist_eccen generates zero probabilities within range for %s'%mod.__name__) 
Example 14
Project: EXOSIMS   Author: dsavransky   File: test_PlanetPopulation.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_dist_albedo(self):
        """
        Test that albedos outside of the range have zero probability

        """

        exclude_mods = ['KeplerLike1',  'AlbedoByRadiusDulzPlavchan',  'DulzPlavchan']
        for mod in self.allmods:
            if (mod.__name__ not in exclude_mods) and ('dist_albedo' in mod.__dict__):
                with RedirectStreams(stdout=self.dev_null):
                    pp = mod(**self.spec)

                p = np.linspace(pp.prange[0]-1,pp.prange[1]+1,100)

                fp = pp.dist_albedo(p)
                self.assertTrue(np.all(fp[p < pp.prange[0]] == 0),'dist_albedo high bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fp[p > pp.prange[1]] == 0),'dist_albedo low bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fp[(p >= pp.prange[0]) & (p <= pp.prange[1])] > 0),'dist_albedo generates zero probabilities within range for %s'%mod.__name__) 
Example 15
Project: EXOSIMS   Author: dsavransky   File: test_PostProcessing.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_ppFact_fits(self):
        # get fits file path for ppFact test
        classpath = os.path.split(inspect.getfile(self.__class__))[0]
        ppFactPath = os.path.join(classpath,'test_PostProcessing_ppFact.fits')

        # fits file has values for WA in [0.1,0.2]
        testWA = np.linspace(0.1,0.2,100)*u.arcsec

        for mod in self.allmods:
            with RedirectStreams(stdout=self.dev_null):
                obj = mod(ppFact=ppFactPath,**self.specs)

            vals = obj.ppFact(testWA)

            self.assertTrue(np.all(vals > 0),'negative value of ppFact for %s'%mod.__name__)
            self.assertTrue(np.all(vals <= 1),'ppFact > 1 for %s'%mod.__name__) 
Example 16
Project: EXOSIMS   Author: dsavransky   File: test_PostProcessing.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_FAdMag0_fits(self):
        # get fits file path for FAdMag0 test
        classpath = os.path.split(inspect.getfile(self.__class__))[0]
        FAdMag0Path = os.path.join(classpath,'test_PostProcessing_FAdMag0.fits')

        # fits file has values for WA in [0.1, 0.2] and FAdMag0 in [10, 20]
        testWA = np.linspace(0.1, 0.2, 100)*u.arcsec

        for mod in self.allmods:
            with RedirectStreams(stdout=self.dev_null):
                obj = mod(FAdMag0=FAdMag0Path,**self.specs)

            vals = obj.FAdMag0(testWA)

            self.assertTrue(np.all(vals >= 10),'value below range of FAdMag0 for %s'%mod.__name__)
            self.assertTrue(np.all(vals <= 20),'value above range of FAdMag0 for %s'%mod.__name__) 
Example 17
Project: EXOSIMS   Author: dsavransky   File: test_Observatory.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_cent(self):
        r"""Test cent method.

        Approach: Probes for a range of inputs.
        """

        print('cent()')
        obs = self.fixture
        # origin at 12:00 on 2000.Jan.01
        t_ref_string = '2000-01-01T12:00:00.0'
        t_ref = Time(t_ref_string, format='isot', scale='utc')
        self.assertEqual(obs.cent(t_ref), 0.0)

        # even-julian-year probes
        t_probe = np.linspace(1950.0, 2050.0, 101)
        for t_ref in t_probe:
            # get the Time object, required by the cent() method
            t_probe_2 = Time(t_ref, format='jyear')
            # exosims century (offset from 2000)
            t_exo = obs.cent(t_probe_2)
            # reference century
            t_ref_cent = (t_ref - 2000.0) / 100.0
            # they are not exactly equal
            self.assertAlmostEqual(t_exo, t_ref_cent, places=10) 
Example 18
Project: EXOSIMS   Author: dsavransky   File: FakeCatalog.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def inverse_method(self,N,d):
        
        t = np.linspace(1e-3,0.999,N)
        f = np.log( t / (1 - t) )
        f = f/f[0]
        
        psi= np.pi*f
        cosPsi = np.cos(psi)
        sinTheta = ( np.abs(cosPsi) + (1-np.abs(cosPsi))*np.random.rand(len(cosPsi)))
        
        theta = np.arcsin(sinTheta)
        theta = np.pi-theta + (2*theta - np.pi)*np.round(np.random.rand(len(t)))
        cosPhi = cosPsi/sinTheta
        phi = np.arccos(cosPhi)*(-1)**np.round(np.random.rand(len(t)))
        
        coords = SkyCoord(phi*u.rad,(np.pi/2-theta)*u.rad,d*np.ones(len(phi))*u.pc)

        return coords 
Example 19
Project: libTLDA   Author: wmkouw   File: viz.py    License: MIT License 5 votes vote down vote up
def plotc(parameters, ax=[], color='k', gridsize=(101, 101)):
    """
    Plot a linear classifier in a 2D scatterplot.

    INPUT   (1) tuple 'parameters': consists of a list of class proportions
                (1 by K classes), an array of class means (K classes by
                D features), an array of class-covariance matrices (D features
                by D features by K classes)
            (2) object 'ax': axes of a pyplot figure or subject (def: empty)
            (3) str 'colors': colors of the contours in the plot (def: 'k')
            (4) tuple 'gridsize': number of points in the grid
                (def: (101, 101))
    OUTPUT  None
    """
    # Check for figure object
    if fig:
        ax = fig.gca()
    else:
        fig, ax = plt.subplots()

    # Get axes limits
    xl = ax.get_xlim()
    yl = ax.get_ylim()

    # Define grid
    gx = np.linspace(xl[0], xl[1], gridsize[0])
    gy = np.linspace(yl[0], yl[1], gridsize[1])
    x, y = np.meshgrid(gx, gy)
    xy = np.vstack((x.ravel(), y.ravel())).T

    # Values of grid
    z = np.dot(xy, parameters[:-1, :]) + parameters[-1, :]
    z = np.reshape(z[:, 0] - z[:, 1], gridsize)

    # Plot grid
    ax.contour(x, y, z, levels=0, colors=colors) 
Example 20
Project: EDeN   Author: fabriziocosta   File: graph_layout.py    License: MIT License 5 votes vote down vote up
def _compute_initial_pos(self, graph):
        _radius = 1
        _offset = 0
        n = len(graph)
        pos = {id: np.array([_radius * math.cos(theta - math.pi / 2) + _offset,
                             _radius * math.sin(theta - math.pi / 2) + _offset]
                            )
               for id, theta in enumerate(
            np.linspace(0, 2 * math.pi * (1 - 1 / float(n)), num=n))}
        return pos 
Example 21
Project: EDeN   Author: fabriziocosta   File: estimator.py    License: MIT License 5 votes vote down vote up
def learning_curve(self, graphs, targets,
                       cv=5, n_steps=10, start_fraction=0.1):
        """learning_curve."""
        graphs, targets = paired_shuffle(graphs, targets)
        x = self.transform(graphs)
        train_sizes = np.linspace(start_fraction, 1.0, n_steps)
        scoring = 'roc_auc'
        train_sizes, train_scores, test_scores = learning_curve(
            self.model, x, targets,
            cv=cv, train_sizes=train_sizes,
            scoring=scoring)
        return train_sizes, train_scores, test_scores 
Example 22
Project: xrft   Author: xgcm   File: test_xrft.py    License: MIT License 5 votes vote down vote up
def test_data_1d(request):
    """Create one dimensional test DataArray."""
    Nx = 16
    Lx = 1.0
    x = np.linspace(0, Lx, Nx)
    dx = x[1] - x[0]
    coords = None if request.param == 'nocoords' else [x]
    da = xr.DataArray(np.random.rand(Nx), coords=coords, dims=['x'])
    if request.param == 'dask':
        da = da.chunk()
    return da 
Example 23
Project: xrft   Author: xgcm   File: test_xrft.py    License: MIT License 5 votes vote down vote up
def test_cross_phase_1d(self, dask):
        N = 32
        x = np.linspace(0, 1, num=N, endpoint=False)
        f = 6
        phase_offset = np.pi/2
        signal1 = np.cos(2*np.pi*f*x)  # frequency = 1/(2*pi)
        signal2 = np.cos(2*np.pi*f*x - phase_offset)
        da1 = xr.DataArray(data=signal1, name='a', dims=['x'], coords={'x': x})
        da2 = xr.DataArray(data=signal2, name='b', dims=['x'], coords={'x': x})

        if dask:
            da1 = da1.chunk({'x': 32})
            da2 = da2.chunk({'x': 32})
        cp = xrft.cross_phase(da1, da2, dim=['x'])

        actual_phase_offset = cp.sel(freq_x=f).values
        npt.assert_almost_equal(actual_phase_offset, phase_offset)
        assert cp.name == 'a_b_phase'

        xrt.assert_equal(xrft.cross_phase(da1, da2), cp)

        with pytest.raises(ValueError):
            xrft.cross_phase(da1, da2.isel(x=0).drop('x'))

        with pytest.raises(ValueError):
            xrft.cross_phase(da1, da2.rename({'x':'y'})) 
Example 24
Project: FRIDA   Author: LCAV   File: generators.py    License: MIT License 5 votes vote down vote up
def gen_sig_at_mic(sigmak2_k, phi_k, pos_mic_x,
                   pos_mic_y, omega_band, sound_speed,
                   SNR, Ns=256):
    """
    generate complex base-band signal received at microphones
    :param sigmak2_k: the variance of the circulant complex Gaussian signal
                emitted by the K sources
    :param phi_k: source locations (azimuths)
    :param pos_mic_x: a vector that contains microphones' x coordinates
    :param pos_mic_y: a vector that contains microphones' y coordinates
    :param omega_band: mid-band (ANGULAR) frequency [radian/sec]
    :param sound_speed: speed of sound
    :param SNR: SNR for the received signal at microphones
    :param Ns: number of snapshots used to estimate the covariance matrix
    :return: y_mic: received (complex) signal at microphones
    """
    num_mic = pos_mic_x.size
    xk, yk = polar2cart(1, phi_k)  # source locations in cartesian coordinates
    # reshape to use broadcasting
    xk = np.reshape(xk, (1, -1), order='F')
    yk = np.reshape(yk, (1, -1), order='F')
    pos_mic_x = np.reshape(pos_mic_x, (-1, 1), order='F')
    pos_mic_y = np.reshape(pos_mic_y, (-1, 1), order='F')

    t = np.reshape(np.linspace(0, 10 * np.pi, num=Ns), (1, -1), order='F')
    K = sigmak2_k.size
    sigmak2_k = np.reshape(sigmak2_k, (-1, 1), order='F')

    # x_tilde_k size: K x length_of_t
    # circular complex Gaussian process
    x_tilde_k = np.sqrt(sigmak2_k / 2.) * (np.random.randn(K, Ns) + 1j *
                                           np.random.randn(K, Ns))
    y_mic = np.dot(np.exp(-1j * (xk * pos_mic_x + yk * pos_mic_y) / (sound_speed / omega_band)),
                   x_tilde_k * np.exp(1j * omega_band * t))
    signal_energy = linalg.norm(y_mic, 'fro') ** 2
    noise_energy = signal_energy / 10 ** (SNR * 0.1)
    sigma2_noise = noise_energy / (Ns * num_mic)
    noise = np.sqrt(sigma2_noise / 2.) * (np.random.randn(*y_mic.shape) + 1j *
                                          np.random.randn(*y_mic.shape))
    y_mic_noisy = y_mic + noise
    return y_mic_noisy, y_mic 
Example 25
Project: dustmaps   Author: gregreen   File: plot_bh.py    License: GNU General Public License v2.0 5 votes vote down vote up
def main():
    w,h = (2056,1024)
    l_0 = 0.

    # Create a grid of coordinates
    print('Creating grid of coordinates...')
    l = np.linspace(-180.+l_0, 180.+l_0, 2*w)
    b = np.linspace(-90., 90., 2*h+2)
    b = b[1:-1]
    l,b = np.meshgrid(l, b)

    l += (np.random.random(l.shape) - 0.5) * 360./(2.*w)
    b += (np.random.random(l.shape) - 0.5) * 180./(2.*h)

    coords = SkyCoord(l*u.deg, b*u.deg, frame='galactic')

    # Set up BH query object
    print('Loading BH map...')
    bh = BHQuery()

    print('Querying map...')
    ebv = bh.query(coords)

    # Convert the output array to a PIL image and save
    print('Saving image...')
    img = numpy2pil(ebv[::-1,::-1], 0., 1.5)
    img = img.resize((w,h), resample=PIL.Image.LANCZOS)
    fname = 'bh.png'
    img.save(fname)

    return 0 
Example 26
Project: dustmaps   Author: gregreen   File: plot_marshall.py    License: GNU General Public License v2.0 5 votes vote down vote up
def main():
    w,h = (2*2056, 2*int(2056*(20./200.)))
    l_0 = 0.

    # Set up MarshallQuery object
    print('Loading Marshall map...')
    query = MarshallQuery()

    # Create a grid of coordinates
    print('Creating grid of coordinates...')
    l = np.linspace(-100.+l_0, 100.+l_0, 2*w)
    b = np.linspace(-10., 10., 2*h)
    dl = l[1] - l[0]
    db = b[1] - b[0]
    l,b = np.meshgrid(l, b)

    l += (np.random.random(l.shape) - 0.5) * dl
    b += (np.random.random(l.shape) - 0.5) * db

    A = np.empty(l.shape+(3,), dtype='f8')

    for k,d in enumerate([1., 2.5, 5.]):
        coords = SkyCoord(l*u.deg, b*u.deg, d*u.kpc, frame='galactic')

        # Get the mean dust extinction at each coordinate
        print('Querying map...')
        A[:,:,k] = query(coords, return_sigma=False)

    A[:,:,2] -= A[:,:,1]
    A[:,:,1] -= A[:,:,0]

    # Convert the output array to a PIL image and save
    print('Saving image...')
    img = numpy2pil(A[::-1,::-1,:], 0., 1., fill=255)
    img = img.resize((w,h), resample=PIL.Image.LANCZOS)
    fname = 'marshall.png'
    img.save(fname)

    return 0 
Example 27
Project: dustmaps   Author: gregreen   File: plot_bayestar.py    License: GNU General Public License v2.0 5 votes vote down vote up
def main():
    w,h = (2056,1024)
    l_0 = 130.

    # Set up Bayestar query object
    print('Loading bayestar map...')
    bayestar = BayestarQuery(max_samples=1)

    # Create a grid of coordinates
    print('Creating grid of coordinates...')
    l = np.linspace(-180.+l_0, 180.+l_0, 2*w)
    b = np.linspace(-90., 90., 2*h+2)
    b = b[1:-1]
    l,b = np.meshgrid(l, b)

    l += (np.random.random(l.shape) - 0.5) * 360./(2.*w)
    b += (np.random.random(l.shape) - 0.5) * 180./(2.*h)

    ebv = np.empty(l.shape+(3,), dtype='f8')

    for k,d in enumerate([0.5, 1.5, 5.]):
        # d = 5.    # We'll query integrated reddening to a distance of 5 kpc
        coords = SkyCoord(l*u.deg, b*u.deg, d*u.kpc, frame='galactic')

        # Get the dust median reddening at each coordinate
        print('Querying map...')
        ebv[:,:,k] = bayestar.query(coords, mode='median')

    ebv[:,:,2] -= ebv[:,:,1]
    ebv[:,:,1] -= ebv[:,:,0]

    # Convert the output array to a PIL image and save
    print('Saving image...')
    img = numpy2pil(ebv[::-1,::-1,:], 0., 1.5)
    img = img.resize((w,h), resample=PIL.Image.LANCZOS)
    fname = 'bayestar.png'
    img.save(fname)

    return 0 
Example 28
Project: dustmaps   Author: gregreen   File: plot_lenz2017.py    License: GNU General Public License v2.0 5 votes vote down vote up
def main():
    w,h = (2056,1024)
    l_0 = 0.

    # Create a grid of coordinates
    print('Creating grid of coordinates...')
    l = np.linspace(-180.+l_0, 180.+l_0, 2*w)
    b = np.linspace(-90., 90., 2*h+2)
    b = b[1:-1]
    l,b = np.meshgrid(l, b)

    l += (np.random.random(l.shape) - 0.5) * 360./(2.*w)
    b += (np.random.random(l.shape) - 0.5) * 180./(2.*h)

    coords = SkyCoord(l*u.deg, b*u.deg, frame='galactic')

    # Set up Lenz+(2017) query object
    print('Loading Lenz+(2017) map...')
    q = Lenz2017Query()

    print('Querying map...')
    ebv = q.query(coords)

    # Convert the output array to a PIL image and save
    print('Saving image...')
    img = numpy2pil(ebv[::-1,::-1], 0., 0.05)
    img = img.resize((w,h), resample=PIL.Image.LANCZOS)
    fname = 'lenz2017.png'
    img.save(fname)

    return 0 
Example 29
Project: dustmaps   Author: gregreen   File: plot_planck.py    License: GNU General Public License v2.0 5 votes vote down vote up
def main():
    w,h = (2056,1024)
    l_0 = 0.

    # Create a grid of coordinates
    print('Creating grid of coordinates...')
    l = np.linspace(-180.+l_0, 180.+l_0, 2*w)
    b = np.linspace(-90., 90., 2*h+2)
    b = b[1:-1]
    l,b = np.meshgrid(l, b)

    l += (np.random.random(l.shape) - 0.5) * 360./(2.*w)
    b += (np.random.random(l.shape) - 0.5) * 180./(2.*h)

    coords = SkyCoord(l*u.deg, b*u.deg, frame='galactic')

    planck_components = [
        ('ebv', 0., 1.5),
        ('radiance', 0., 1.5),
        ('tau', 0., 1.5),
        ('temp', 15.*u.K, 25.*u.K),
        ('err_temp', 0.*u.K, 4.*u.K),
        ('beta', 1., 3.),
        ('err_beta', 0., 0.2)]

    for component,vmin,vmax in planck_components:
        # Set up Planck query object
        print('Loading Planck map...')
        planck = PlanckQuery(component=component)

        print('Querying map...')
        res = planck.query(coords)

        # Convert the output array to a PIL image and save
        print('Saving image...')
        img = numpy2pil(res[::-1,::-1], vmin, vmax)
        img = img.resize((w,h), resample=PIL.Image.LANCZOS)
        fname = 'planck_{}.png'.format(component)
        img.save(fname)

    return 0 
Example 30
Project: dustmaps   Author: gregreen   File: plot_sfd.py    License: GNU General Public License v2.0 5 votes vote down vote up
def main():
    w,h = (2056,1024)
    l_0 = 0.

    # Create a grid of coordinates
    print('Creating grid of coordinates...')
    l = np.linspace(-180.+l_0, 180.+l_0, 2*w)
    b = np.linspace(-90., 90., 2*h+2)
    b = b[1:-1]
    l,b = np.meshgrid(l, b)

    l += (np.random.random(l.shape) - 0.5) * 360./(2.*w)
    b += (np.random.random(l.shape) - 0.5) * 180./(2.*h)

    coords = SkyCoord(l*u.deg, b*u.deg, frame='galactic')

    # Set up SFD query object
    print('Loading SFD map...')
    sfd = SFDQuery()

    print('Querying map...')
    ebv = sfd.query(coords)

    # Convert the output array to a PIL image and save
    print('Saving image...')
    img = numpy2pil(ebv[::-1,::-1], 0., 1.5)
    img = img.resize((w,h), resample=PIL.Image.LANCZOS)
    fname = 'sfd.png'
    img.save(fname)

    return 0