numpy.where

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

172 Examples 7

Example 1

Project: finmarketpy Source File: techindicator.py
    def create_tech_ind(self, data_frame_non_nan, name, tech_params, data_frame_non_nan_early = None):
        self._signal = None
        self._techind = None

        data_frame = data_frame_non_nan.fillna(method="ffill")

        if data_frame_non_nan_early is not None:
            data_frame_early = data_frame_non_nan_early.fillna(method="ffill")

        if name == "SMA":

            if (data_frame_non_nan_early is not None):
                # calculate the lagged sum of the n-1 point
                rolling_sum = data_frame.shift(1).rolling(center=False, window=tech_params.sma_period - 1).sum()

                # add non-nan one for today
                rolling_sum = rolling_sum + data_frame_early

                # calculate average = sum / n
                self._techind = rolling_sum / tech_params.sma_period

                narray = numpy.where(data_frame_early > self._techind, 1, -1)
            else:
                # self._techind = pandas.rolling_mean(data_frame, tech_params.sma_period)
                self._techind = data_frame.rolling(window=tech_params.sma_period, center=False).mean()

                narray = numpy.where(data_frame > self._techind, 1, -1)

            self._signal = pandas.DataFrame(index = data_frame.index, data = narray)
            self._signal.loc[0:tech_params.sma_period] = numpy.nan
            self._signal.columns = [x + " SMA Signal" for x in data_frame.columns.values]

            self._techind.columns = [x + " SMA" for x in data_frame.columns.values]

        elif name == "EMA":

            # self._techind = pandas.ewma(data_frame, span = tech_params.ema_period)
            self._techind = data_frame.ewm(ignore_na=False,span=tech_params.ema_period,min_periods=0,adjust=True).mean()

            narray = numpy.where(data_frame > self._techind, 1, -1)

            self._signal = pandas.DataFrame(index = data_frame.index, data = narray)
            self._signal.loc[0:tech_params.ema_period] = numpy.nan
            self._signal.columns = [x + " EMA Signal" for x in data_frame.columns.values]

            self._techind.columns = [x + " EMA" for x in data_frame.columns.values]

        elif name == "ROC":

            if (data_frame_non_nan_early is not None):
                self._techind = data_frame_early / data_frame.shift(tech_params.roc_period) - 1
            else:
                self._techind = data_frame / data_frame.shift(tech_params.roc_period) - 1

            narray = numpy.where(self._techind > 0, 1, -1)

            self._signal = pandas.DataFrame(index = data_frame.index, data = narray)
            self._signal.loc[0:tech_params.roc_period] = numpy.nan
            self._signal.columns = [x + " ROC Signal" for x in data_frame.columns.values]

            self._techind.columns = [x + " ROC" for x in data_frame.columns.values]

        elif name == "polarity":
            self._techind = data_frame

            narray = numpy.where(self._techind > 0, 1, -1)

            self._signal = pandas.DataFrame(index = data_frame.index, data = narray)
            self._signal.columns = [x + " Polarity Signal" for x in data_frame.columns.values]

            self._techind.columns = [x + " Polarity" for x in data_frame.columns.values]

        elif name == "SMA2":
            sma = data_frame.rolling(window=tech_params.sma_period, center=False).mean()
            sma2 = data_frame.rolling(window=tech_params.sma2_period, center=False).mean()

            narray = numpy.where(sma > sma2, 1, -1)

            self._signal = pandas.DataFrame(index = data_frame.index, data = narray)
            self._signal.columns = [x + " SMA2 Signal" for x in data_frame.columns.values]

            sma.columns = [x + " SMA" for x in data_frame.columns.values]
            sma2.columns = [x + " SMA2" for x in data_frame.columns.values]
            most = max(tech_params.sma_period, tech_params.sma2_period)
            self._signal.loc[0:most] = numpy.nan
            self._techind = pandas.concat([sma, sma2], axis = 1)

        elif name in ['RSI']:
            # delta = data_frame.diff()
            #
            # dUp, dDown = delta.copy(), delta.copy()
            # dUp[dUp < 0] = 0
            # dDown[dDown > 0] = 0
            #
            # rolUp = pandas.rolling_mean(dUp, tech_params.rsi_period)
            # rolDown = pandas.rolling_mean(dDown, tech_params.rsi_period).abs()
            #
            # rsi = rolUp / rolDown

            # Get the difference in price from previous step
            delta = data_frame.diff()
            # Get rid of the first row, which is NaN since it did not have a previous
            # row to calculate the differences
            delta = delta[1:]

            # Make the positive gains (up) and negative gains (down) Series
            up, down = delta.copy(), delta.copy()
            up[up < 0] = 0
            down[down > 0] = 0

            # Calculate the EWMA
            roll_up1 = pandas.stats.moments.ewma(up, tech_params.rsi_period)
            roll_down1 = pandas.stats.moments.ewma(down.abs(), tech_params.rsi_period)

            # Calculate the RSI based on EWMA
            RS1 = roll_up1 / roll_down1
            RSI1 = 100.0 - (100.0 / (1.0 + RS1))

            # Calculate the SMA
            roll_up2 = up.rolling(window=tech_params.rsi_period, center=False).mean()
            roll_down2 = down.abs().rolling(window=tech_params.rsi_period, center=False).mean()

            # Calculate the RSI based on SMA
            RS2 = roll_up2 / roll_down2
            RSI2 = 100.0 - (100.0 / (1.0 + RS2))

            self._techind = RSI2
            self._techind.columns = [x + " RSI" for x in data_frame.columns.values]

            signal = data_frame.copy()

            sells = (signal.shift(-1) < tech_params.rsi_lower) & (signal > tech_params.rsi_lower)
            buys = (signal.shift(-1) > tech_params.rsi_upper) & (signal < tech_params.rsi_upper)

            # print (buys[buys == True])

            # buys
            signal[buys] =  1
            signal[sells] = -1
            signal[~(buys | sells)] = numpy.nan
            signal = signal.fillna(method = 'ffill')

            self._signal = signal

            self._signal.loc[0:tech_params.rsi_period] = numpy.nan
            self._signal.columns = [x + " RSI Signal" for x in data_frame.columns.values]

        elif name in ["BB"]:
            ## calcuate Bollinger bands
            mid = data_frame.rolling(center=False, window=tech_params.bb_period).mean(); mid.columns = [x + " BB Mid" for x in data_frame.columns.values]
            std_dev = data_frame.rolling(center=False, window=tech_params.bb_period).std()
            BB_std = tech_params.bb_mult * std_dev

            lower = pandas.DataFrame(data = mid.values - BB_std.values, index = mid.index,
                            columns = data_frame.columns)

            upper = pandas.DataFrame(data = mid.values + BB_std.values, index = mid.index,
                            columns = data_frame.columns)

            ## calculate signals
            signal = data_frame.copy()

            buys = signal > upper
            sells = signal < lower

            signal[buys] = 1
            signal[sells] = -1
            signal[~(buys | sells)] = numpy.nan
            signal = signal.fillna(method = 'ffill')

            self._signal = signal
            self._signal.loc[0:tech_params.bb_period] = numpy.nan
            self._signal.columns = [x + " " + name + " Signal" for x in data_frame.columns.values]

            lower.columns = [x + " BB Lower" for x in data_frame.columns.values]
            upper.columns = [x + " BB Mid" for x in data_frame.columns.values]
            upper.columns = [x + " BB Lower" for x in data_frame.columns.values]

            self._techind = pandas.concat([lower, mid, upper], axis = 1)
        elif name == "long-only":
            ## have +1 signals only
            self._techind = data_frame  # the technical indicator is just "prices"

            narray = numpy.ones((len(data_frame.index), len(data_frame.columns)))

            self._signal = pandas.DataFrame(index = data_frame.index, data = narray)
            self._signal.columns = [x + " Long Only Signal" for x in data_frame.columns.values]

            self._techind.columns = [x + " Long Only" for x in data_frame.columns.values]

        self.create_custom_tech_ind(data_frame_non_nan, name, tech_params, data_frame_non_nan_early)

        # TODO create other indicators
        if hasattr(tech_params, 'only_allow_longs'):
            self._signal[self._signal < 0] = 0

        # TODO create other indicators
        if hasattr(tech_params, 'only_allow_shorts'):
           self._signal[self._signal > 0] = 0

        # apply signal multiplier (typically to flip signals)
        if hasattr(tech_params, 'signal_mult'):
            self._signal = self._signal * tech_params.signal_mult

        if hasattr(tech_params, 'strip_signal_name'):
            if tech_params.strip_signal_name:
                self._signal.columns = data_frame.columns

        return self._techind

Example 2

Project: pyscf Source File: geom.py
    def search_possible_rotations(self, zaxis=None):
        '''If zaxis is given, the rotation axis is parallel to zaxis'''
        maybe_cn = []
        for lst in self.group_atoms_by_distance:
            natm = len(lst)
            if natm > 1:
                coords = self.atoms[lst,1:]
# possible C2 axis
                for i in range(1, natm):
                    if abs(coords[0]+coords[i]).sum() > TOLERANCE:
                        maybe_cn.append((coords[0]+coords[i], 2))
                    else: # abs(coords[0]-coords[i]).sum() > TOLERANCE:
                        maybe_cn.append((coords[0]-coords[i], 2))

# atoms of equal distances may be associated with rotation axis > C2.
                r0 = coords - coords[0]
                distance = pyscf.lib.norm(r0, axis=1)
                eq_distance = abs(distance[:,None] - distance) < TOLERANCE
                for i in range(2, natm):
                    for j in numpy.where(eq_distance[i,:i])[0]:
                        cos = numpy.dot(r0[i],r0[j]) / (distance[i]*distance[j])
                        ang = numpy.arccos(cos)
                        nfrac = numpy.pi*2 / (numpy.pi-ang)
                        n = int(numpy.around(nfrac))
                        if abs(nfrac-n) < TOLERANCE:
                            maybe_cn.append((numpy.cross(r0[i],r0[j]),n))

        # remove zero-vectors and duplicated vectors
        vecs = numpy.vstack([x[0] for x in maybe_cn])
        idx = pyscf.lib.norm(vecs, axis=1) > TOLERANCE
        ns = numpy.hstack([x[1] for x in maybe_cn])
        vecs = _normalize(vecs[idx])
        ns = ns[idx]

        if zaxis is not None:  # Keep parallel rotation axes
            cos = numpy.dot(vecs, _normalize(zaxis))
            vecs = vecs[(abs(cos-1) < TOLERANCE) | (abs(cos+1) < TOLERANCE)]
            ns = ns[(abs(cos-1) < TOLERANCE) | (abs(cos+1) < TOLERANCE)]

        possible_cn = []
        seen = numpy.zeros(len(vecs), dtype=bool)
        for k, v in enumerate(vecs):
            if not seen[k]:
                where1 = numpy.einsum('ix->i', abs(vecs[k:] - v)) < TOLERANCE
                where1 = numpy.where(where1)[0] + k
                where2 = numpy.einsum('ix->i', abs(vecs[k:] + v)) < TOLERANCE
                where2 = numpy.where(where2)[0] + k
                seen[where1] = True
                seen[where2] = True

                vk = _normalize((numpy.einsum('ix->x', vecs[where1]) -
                                 numpy.einsum('ix->x', vecs[where2])))
                for n in (set(ns[where1]) | set(ns[where2])):
                    possible_cn.append((vk,n))
        return possible_cn

Example 3

Project: PyGazeAnalyser Source File: traces.py
def interpolate_blink(signal, mode='auto', velthresh=5, maxdur=500, margin=10, invalid=-1, edfonly=False):
	
	"""Returns signal with interpolated results, based on a cubic or linear
	interpolation of all blinks detected in the signal; based on:
	https://github.com/smathot/exparser/blob/master/exparser/TraceKit.py
	
	arguments
	signal	--	a vector (i.e. a NumPy array) containing a single
				trace of your signal; alternatively a trial gaze data
				dict as is returned by edfreader can be passed; in this
				case the blink ending events will be used to find blinks
				before the pupil size velocity algorithm will be used
				(NOTE: this means both will be used successively!)
	
	keyword arguments
	mode		--	string indicating what kind of interpolation to use:
				'linear' for a linear interpolation
				'cubic' for a cubic interpolation
				'auto' for a cubic interpolation is possible (i.e.
					when more than four data points are available)
					and linear when this is not the case
				(default = 'auto')
	velthresh	--	pupil size change velocity threshold in arbitrary
				units per sample (default = 5)
	maxdur	--	maximal duration of the blink in samples
				(default = 500)
	margin	--	margin (in samples) to compensate for blink duration
				underestimatiom; blink is extended for detected start
				minus margin, and detected end plus margin
				(default = 10)
	edfonly	--	Boolean indicating whether blinks should ONLY be
				detected using the EDF logs and NOT algorithmically
	
	returns
	signal	--	a NumPy array containing the interpolated signal
	"""
	
	# # # # #
	# input errors
	
	# wrong interpolation method	
	if mode not in ['auto','linear','cubic']:
		raise Exception("Error in pyenalysis.interpolate_missing: mode '%s' is not supported, please use one of the following: 'auto','linear','cubic'" % mode)
	# wrong signal dimension
	if type(signal) != dict:
		if signal.ndim != 1:
			raise Exception("Error in pyenalysis.interpolate_missing: input is not a single signal trace, but has %d dimensions; please provide a 1-dimension array" % signal.ndim)

	# # # # #
	# find blinks
	
	# empty lists, to store blink starts and endings
	starts = []
	ends = []
	
	# edfreader data
	if type(signal) == dict:

		# loop through blinks
		for st, et, dur in signal['events']['Eblk']: # Eblk - list of lists, each containing [starttime, endtime, duration]
			
			# edf time to sample number
			st = numpy.where(signal['edftime']==st)[0]
			et = numpy.where(signal['edftime']==et)[0]
			# if the starting or ending time did not appear in the trial,
			# correct the blink starting or ending point to the first or
			# last sample, respectively
			if len(st) == 0:
				st = 0
			else:
				st = st[0]
			if len(et) == 0:
				et = len(signal['edftime'])
			else:
				et = et[0]
			# compensate for underestimation of blink duration
			if st-margin >= 0:
				st -= margin
			if et+margin < len(signal['size']):
				et += margin
			# do not except blinks that exceed maximal blink duration
			if et-st <= maxdur:
				# append start time and ending time
				starts.append(st)
				ends.append(et)
		# extract pupil size data from signal
		signal = signal['size']
	
	if not edfonly:
		# signal in NumPy array
		# create a velocity profile of the signal
		vprof = signal[1:]-signal[:-1]
		
		# start detection
		ifrom = 0
		while True:
			# blink onset is detected when pupil size change velocity exceeds
			# threshold
			l = numpy.where(vprof[ifrom:] < -velthresh)[0]
			# break when no blink start is detected
			if len(l) == 0:
				break
			# blink start index
			istart = l[0]+ifrom
			if ifrom == istart:
				break
			# reversal (opening of the eye) is detected when pupil size
			# starts to increase with a super-threshold velocity
			l = numpy.where(vprof[istart:] > velthresh)[0]
			# if no reversal is detected, start detection process at istart
			# next run
			if len(l) == 0:
				ifrom = istart
				# reloop
				continue
			# index number of somewhat halfway blink process
			imid = l[0] + istart
			# a blink ending is detected when pupil size increase velocity
			# falls back to zero
			l = numpy.where(vprof[imid:] < 0)[0]
			# if no ending is detected, start detection process from imid
			# next run
			if len(l) == 0:
				ifrom = imid
				# reloop
				continue
			# blink end index
			iend = l[0]+imid
			# start detection process from current blink ending next run
			ifrom = iend
			# compensate for underestimation of blink duration
			if istart-margin >= 0:
				istart -= margin
			if iend+margin < len(signal):
				iend += margin
			# do not except blinks that exceed maximal blink duration
			if iend-istart > maxdur:
				# reloop
				continue
			# if all is well, we append start and ending to their respective
			# lists
			starts.append(istart)
			ends.append(iend)

#	# DEBUG #
#	pyplot.figure()
#	pyplot.title("" % ())
#	pyplot.plot(signal,'ko')
#	pyplot.plot(vprof,'b')
#	# # # # #

	# # # # #
	# interpolate

	# loop through all starting and ending positions
	for i in range(len(starts)):
		
		# empty list to store data points for interpolation
		pl = []
		
		# duration in samples
		duration = ends[i]-starts[i]

		# starting point
		if starts[i] - duration >= 0:
			pl.extend([starts[i]-duration])
		# central points (data between these points will be replaced)
		pl.extend([starts[i],ends[i]])
		# ending point
		if ends[i] + duration < len(signal):
			pl.extend([ends[i]+duration])
		
		# choose interpolation type
		if mode == 'auto':
			# if our range is wide enough, we can interpolate cubicly
			if len(pl) >= 4:
				kind = 'cubic'
			# if not, we use a linear interpolation
			else:
				kind = 'linear'
		else:
			kind = mode[:]
		
		# select values for interpolation function
		x = numpy.array(pl)
		y = signal[x]
		
		# replace any invalid values with trial average
		y[y==invalid] = numpy.mean(signal[signal!=invalid])
		
		# create interpolation function
		intfunc = interp1d(x,y,kind=kind)
		
		# do interpolation
		xint = numpy.arange(starts[i],ends[i])
		yint = intfunc(xint)
		
		# insert interpolated values into signal
		signal[xint] = yint
	
#		# DEBUG #
#		y = numpy.zeros(len(pl)) + max(signal)
#		pyplot.plot(pl,y,'ro')
#	pyplot.plot(signal,'r')
#	# # # # #
		
	return signal

Example 4

Project: biggles Source File: func.py
def get_log_plot_range(x, err=None, input_range=None, get_good=False):
    """
    Get a plot range in the case of log axes
    """
    if input_range is not None:
        if len(input_range) < 2:
            raise ValueError("expected [xmin,xmax] for input range")
        if input_range[0] <= 0. or input_range[1] <= 0.:
            raise ValueError("cannot use plot range < 0 for log plots, got [%s,%s]" % tuple(input_range))
        if get_good:
            w,=numpy.where((x >= input_range[0]) & (x <= input_range[1]))
            return input_range, w
        else:
            return input_range

    w,=numpy.where(x > 0.)
    if w.size == 0:
        raise ValueError("No values are greater than zero in log plot")

    minval = min(x[w])
    if err is not None:
        w2, = numpy.where( (x[w] - err[w]) > 0 )
        if w2.size > 0:
            minval2 =  min(x[w[w2]] - err[w[w2]])
            minval = min(minval,minval2)

        maxval = max(x+err)
    else:
        maxval = max(x)

    minval *= 0.5
    maxval *= 2

    if get_good:
        return [minval,maxval], w
    else:
        return [minval,maxval]

Example 5

Project: pyscf Source File: newton_ah.py
def gen_g_hop_uhf(mf, mo_coeff, mo_occ, fock_ao=None):
    mol = mf.mol
    occidxa = numpy.where(mo_occ[0]>0)[0]
    occidxb = numpy.where(mo_occ[1]>0)[0]
    viridxa = numpy.where(mo_occ[0]==0)[0]
    viridxb = numpy.where(mo_occ[1]==0)[0]
    nocca = len(occidxa)
    noccb = len(occidxb)
    nvira = len(viridxa)
    nvirb = len(viridxb)

    if fock_ao is None:
        dm0 = mf.make_rdm1(mo_coeff, mo_occ)
        fock_ao = mf.get_hcore() + mf.get_veff(mol, dm0)
    focka = reduce(numpy.dot, (mo_coeff[0].T, fock_ao[0], mo_coeff[0]))
    fockb = reduce(numpy.dot, (mo_coeff[1].T, fock_ao[1], mo_coeff[1]))

    g = numpy.hstack((focka[viridxa[:,None],occidxa].reshape(-1),
                      fockb[viridxb[:,None],occidxb].reshape(-1)))

    h_diaga =(focka[viridxa,viridxa].reshape(-1,1) - focka[occidxa,occidxa])
    h_diagb =(fockb[viridxb,viridxb].reshape(-1,1) - fockb[occidxb,occidxb])
    h_diag = numpy.hstack((h_diaga.reshape(-1), h_diagb.reshape(-1)))

    if hasattr(mf, 'xc') and hasattr(mf, '_numint'):
        if mf.grids.coords is None:
            mf.grids.build()
        hyb = mf._numint.hybrid_coeff(mf.xc, spin=(mol.spin>0)+1)
        rho0, vxc, fxc = mf._numint.cache_xc_kernel(mol, mf.grids, mf.xc,
                                                    mo_coeff, mo_occ, 1)
        #dm0 =(numpy.dot(mo_coeff[0]*mo_occ[0], mo_coeff[0].T),
        #      numpy.dot(mo_coeff[1]*mo_occ[1], mo_coeff[1].T))
        dm0 = None
    else:
        hyb = None

    def h_op(x):
        x1a = x[:nvira*nocca].reshape(nvira,nocca)
        x1b = x[nvira*nocca:].reshape(nvirb,noccb)
        x2a  = numpy.einsum('rp,rq->pq', focka[viridxa[:,None],viridxa], x1a)
        x2a -= numpy.einsum('sq,ps->pq', focka[occidxa[:,None],occidxa], x1a)
        x2b  = numpy.einsum('rp,rq->pq', fockb[viridxb[:,None],viridxb], x1b)
        x2b -= numpy.einsum('sq,ps->pq', fockb[occidxb[:,None],occidxb], x1b)

        d1a = reduce(numpy.dot, (mo_coeff[0][:,viridxa], x1a,
                                 mo_coeff[0][:,occidxa].T))
        d1b = reduce(numpy.dot, (mo_coeff[1][:,viridxb], x1b,
                                 mo_coeff[1][:,occidxb].T))
        dm1 = numpy.array((d1a+d1a.T,d1b+d1b.T))
        if hyb is None:
            v1 = mf.get_veff(mol, dm1)
        else:
            v1 = mf._numint.nr_uks_fxc(mol, mf.grids, mf.xc, dm0, dm1,
                                       0, 1, rho0, vxc, fxc)
            if abs(hyb) < 1e-10:
                vj = mf.get_j(mol, dm1)
                v1 += vj[0] + vj[1]
            else:
                vj, vk = mf.get_jk(mol, dm1)
                v1 += vj[0]+vj[1] - vk * hyb * .5
        x2a += reduce(numpy.dot, (mo_coeff[0][:,viridxa].T, v1[0],
                                  mo_coeff[0][:,occidxa]))
        x2b += reduce(numpy.dot, (mo_coeff[1][:,viridxb].T, v1[1],
                                  mo_coeff[1][:,occidxb]))
        return numpy.hstack((x2a.ravel(), x2b.ravel()))

    return g, h_op, h_diag

Example 6

Project: biggles Source File: func.py
    def _set_range_and_subpts(self):
        """
        keys['xerr'] and 'yerr' could be modified to be 1-d arrays in the range,
        so make sure keys was already a copy of original keys
        """

        x=self.x
        y=self.y
        xerr=self.xerr
        yerr=self.yerr

        xrng = self.get('xrange',None)
        yrng = self.get('yrange',None)

        # For log, Don't plot points less than zero
        w=None
        if self['xlog'] and self['ylog']:
            xrng = get_log_plot_range(x, err=xerr, input_range=xrng)
            yrng = get_log_plot_range(y, err=yerr, input_range=yrng)
            w,=numpy.where( (x > xrng[0]) & (y > yrng[0]) )
        elif self['xlog']:
            xrng = get_log_plot_range(x, err=xerr, input_range=xrng)
            w,=numpy.where( x > xrng[0])
        elif self['ylog']:
            yrng = get_log_plot_range(y, err=yerr, input_range=yrng)
            w,=numpy.where( y > yrng[0])

        if w is not None:
            if w.size == 0:
                raise ValueError("no points are in range")
        self.indices=w
        self.xrng=xrng
        self.yrng=yrng

Example 7

Project: nupic.research Source File: temporal_pooler.py
  def _adaptSynapses(self, inputVector, activeColumns, predictedActiveCells):
    """
    This is the primary learning method. It updates synapses' permanence based
    on the bottom-up input to the TP and the TP's active cells.
    For each active cell, its synapses' permanences are updated as follows:

    1. if pre-synaptic input is ON due to a correctly predicted cell,
       increase permanence by _synPredictedInc
    2. else if input is ON due to an active cell, increase permanence by
       _synPermActiveInc
    3. else input is OFF, decrease permanence by _synPermInactiveDec

    Parameters:
    ----------------------------
    inputVector:    a numpy array whose ON bits represent the active cells from
                    temporal memory
    activeColumns:  an array containing the indices of the columns that 
                    survived the inhibition step
    predictedActiveCells: a numpy array with numInputs elements. A 1 indicates
                          that this cell switched from predicted state in
                          the previous time step to active state in the current
                          timestep
    """
    inputIndices = numpy.where(inputVector > 0)[0]
    predictedIndices = numpy.where(predictedActiveCells > 0)[0]
    permChanges = numpy.zeros(self._numInputs)

    # Decrement inactive TM cell -> active TP cell connections
    permChanges.fill(-1 * self._synPermInactiveDec)

    # Increment active TM cell -> active TP cell connections
    permChanges[inputIndices] = self._synPermActiveInc

    # Increment correctly predicted TM cell -> active TP cell connections
    permChanges[predictedIndices] = self._synPredictedInc

    if self._spVerbosity > 4:
      print "\n============== _adaptSynapses ======"
      print "Active input indices:",inputIndices
      print "predicted input indices:",predictedIndices
      print "\n============== _adaptSynapses ======\n"

    for i in activeColumns:
      # Get the permanences of the synapses of TP cell i
      perm = self._permanences.getRow(i)

      # Only consider connections in column's potential pool (receptive field)
      maskPotential = numpy.where(self._potentialPools.getRow(i) > 0)[0]
      perm[maskPotential] += permChanges[maskPotential]
      self._updatePermanencesForColumn(perm, i, raisePerm=False)

Example 8

Project: academicmarkdown Source File: tools.py
def addLineNumbersToPDF(inFile, outFile, color='#d3d7cf'):

	"""
	desc:
		Adds line numbers to a PDF file.

	arguments:
		inFile:
			desc:	The name of the input PDF.
			type:	str, unicode
		outFile:
			desc:	The name of the output PDF.
			type:	str, unicode

	keywords:
		color:
			desc:	An HTML-style color name.
			type:	str, unicode
	"""

	import os
	import shutil
	import subprocess
	from scipy import ndimage
	import numpy as np
	from PIL import Image, ImageDraw, ImageFont

	#fontFile = '/usr/share/fonts/truetype/msttcorefonts/Times_New_Roman.ttf'
	fontFile = '/usr/share/fonts/truetype/freefont/FreeSans.ttf'
	fontSize = 20
	tmpFolder = u'line-numbers-tmp'
	pageFolder = u'%s/page' % tmpFolder
	watermarkFolder = u'%s/watermark' % tmpFolder

	try:
		shutil.rmtree(tmpFolder)
	except:
		pass
	os.makedirs(pageFolder)
	os.makedirs(watermarkFolder)

	print(u'Adding line numbers to PDF')
	print(u'Converting ...')
	cmd = u'convert -density 150 %s %s' % (inFile, os.path.join(pageFolder,
		u'%03d.png'))
	subprocess.call(cmd.split())
	print(u'Done!')
	# Create watermarks for all pages
	for path in os.listdir(pageFolder):
		im = ndimage.imread(os.path.join(pageFolder, path), flatten=True)
		# Create a list of indices that have text on them
		nonEmptyRows = np.where(im.mean(axis=1) != 255)[0]
		# Store the rows (i.e.) y coordinates of all to-be-numbered-rows
		numberRows =[]
		firstRow = None
		for row in nonEmptyRows:
			if im[row-1].mean() == 255:
				numberRows.append(row)
		print(u'Found %d lines!' % len(numberRows))
		# Create watermark image
		print(u'Creating watermark ...')
		font = ImageFont.truetype(fontFile, fontSize)
		wm = Image.new('RGBA', (im.shape[1], im.shape[0]))
		dr = ImageDraw.Draw(wm)
		i = 1
		for row in numberRows:
			dr.text((32, row), '%s' % i, font=font, fill=color)
			i += 1
		wm.save(os.path.join(watermarkFolder, path))
		print(u'Done!')

	print(u'Creating watermark pdf ...')
	cmd = 'convert %s/*.png watermark.pdf' % watermarkFolder
	subprocess.call(cmd.split())
	print(u'Done!')

	print(u'Merging watermark and source docuement ...')
	cmd = u'pdftk %s multibackground watermark.pdf output %s' \
		% (inFile, outFile)
	subprocess.call(cmd.split())
	print(u'Done!')

	print(u'Cleaning up ...')
	shutil.rmtree(tmpFolder)
	print(u'Done')

Example 9

Project: summpy Source File: lexrank.py
def lexrank(sentences, continuous=False, sim_threshold=0.1, alpha=0.9,
            use_divrank=False, divrank_alpha=0.25):
    '''
    compute centrality score of sentences.

    Args:
      sentences: [u'こんにちは.', u'私の名前は飯沼です.', ... ]
      continuous: if True, apply continuous LexRank. (see reference)
      sim_threshold: if continuous is False and smilarity is greater or
        equal to sim_threshold, link the sentences.
      alpha: the damping factor of PageRank and DivRank
      divrank: if True, apply DivRank instead of PageRank
      divrank_alpha: strength of self-link [0.0-1.0]
        (it's not the damping factor, see divrank.py)

    Returns: tuple
      (
        {
          # sentence index -> score
          0: 0.003,
          1: 0.002,
          ...
        },
        similarity_matrix
      )
    
    Reference:
      Günes Erkan and Dragomir R. Radev.
      LexRank: graph-based lexical centrality as salience in text
      summarization. (section 3)
      http://www.cs.cmu.edu/afs/cs/project/jair/pub/volume22/erkan04a-html/erkan04a.html
    '''
    # configure ranker
    ranker_params = {'max_iter': 1000}
    if use_divrank:
        ranker = divrank_scipy
        ranker_params['alpha'] = divrank_alpha
        ranker_params['d'] = alpha
    else:
        ranker = networkx.pagerank_scipy
        ranker_params['alpha'] = alpha

    graph = networkx.DiGraph()

    # sentence -> tf
    sent_tf_list = []
    for sent in sentences:
        words = tools.word_segmenter_ja(sent)
        tf = collections.Counter(words)
        sent_tf_list.append(tf)

    sent_vectorizer = DictVectorizer(sparse=True)
    sent_vecs = sent_vectorizer.fit_transform(sent_tf_list)

    # compute similarities between senteces
    sim_mat = 1 - pairwise_distances(sent_vecs, sent_vecs, metric='cosine')

    if continuous:
        linked_rows, linked_cols = numpy.where(sim_mat > 0)
    else:
        linked_rows, linked_cols = numpy.where(sim_mat >= sim_threshold)

    # create similarity graph
    graph.add_nodes_from(range(sent_vecs.shape[0]))
    for i, j in zip(linked_rows, linked_cols):
        if i == j:
            continue
        weight = sim_mat[i,j] if continuous else 1.0
        graph.add_edge(i, j, {'weight': weight})

    scores = ranker(graph, **ranker_params)
    return scores, sim_mat

Example 10

Project: opticspy Source File: axis3d.py
    def draw(self, renderer):
        self.label._transform = self.axes.transData
        renderer.open_group('axis3d')

        # code from XAxis
        majorTicks = self.get_major_ticks()
        majorLocs = self.major.locator()

        info = self._axinfo
        index = info['i']

        # filter locations here so that no extra grid lines are drawn
        locmin, locmax = self.get_view_interval()
        if locmin > locmax:
            locmin, locmax = locmax, locmin

        # Rudimentary clipping
        majorLocs = [loc for loc in majorLocs if
                     locmin <= loc <= locmax]
        self.major.formatter.set_locs(majorLocs)
        majorLabels = [self.major.formatter(val, i)
                       for i, val in enumerate(majorLocs)]

        mins, maxs, centers, deltas, tc, highs = self._get_coord_info(renderer)

        # Determine grid lines
        minmax = np.where(highs, maxs, mins)

        # Draw main axis line
        juggled = info['juggled']
        edgep1 = minmax.copy()
        edgep1[juggled[0]] = get_flip_min_max(edgep1, juggled[0], mins, maxs)

        edgep2 = edgep1.copy()
        edgep2[juggled[1]] = get_flip_min_max(edgep2, juggled[1], mins, maxs)
        pep = proj3d.proj_trans_points([edgep1, edgep2], renderer.M)
        centpt = proj3d.proj_transform(centers[0], centers[1], centers[2], renderer.M)
        self.line.set_data((pep[0][0], pep[0][1]), (pep[1][0], pep[1][1]))
        self.line.draw(renderer)

        # Grid points where the planes meet
        xyz0 = []
        for val in majorLocs:
            coord = minmax.copy()
            coord[index] = val
            xyz0.append(coord)

        # Draw labels
        peparray = np.asanyarray(pep)
        # The transAxes transform is used because the Text object
        # rotates the text relative to the display coordinate system.
        # Therefore, if we want the labels to remain parallel to the
        # axis regardless of the aspect ratio, we need to convert the
        # edge points of the plane to display coordinates and calculate
        # an angle from that.
        # TODO: Maybe Text objects should handle this themselves?
        dx, dy = (self.axes.transAxes.transform(peparray[0:2, 1]) -
                  self.axes.transAxes.transform(peparray[0:2, 0]))

        lxyz = 0.5*(edgep1 + edgep2)

        labeldeltas = info['label']['space_factor'] * deltas
        axmask = [True, True, True]
        axmask[index] = False
        lxyz = move_from_center(lxyz, centers, labeldeltas, axmask)
        tlx, tly, tlz = proj3d.proj_transform(lxyz[0], lxyz[1], lxyz[2], \
                renderer.M)
        self.label.set_position((tlx, tly))
        if self.get_rotate_label(self.label.get_text()):
            angle = art3d.norm_text_angle(math.degrees(math.atan2(dy, dx)))
            self.label.set_rotation(angle)
        self.label.set_va(info['label']['va'])
        self.label.set_ha(info['label']['ha'])
        self.label.draw(renderer)


        # Draw Offset text

        # Which of the two edge points do we want to
        # use for locating the offset text?
        if juggled[2] == 2 :
            outeredgep = edgep1
            outerindex = 0
        else :
            outeredgep = edgep2
            outerindex = 1

        pos = copy.copy(outeredgep)
        pos = move_from_center(pos, centers, labeldeltas, axmask)
        olx, oly, olz = proj3d.proj_transform(pos[0], pos[1], pos[2], renderer.M)
        self.offsetText.set_text( self.major.formatter.get_offset() )
        self.offsetText.set_position( (olx, oly) )
        angle = art3d.norm_text_angle(math.degrees(math.atan2(dy, dx)))
        self.offsetText.set_rotation(angle)
        # Must set rotation mode to "anchor" so that
        # the alignment point is used as the "fulcrum" for rotation.
        self.offsetText.set_rotation_mode('anchor')

        #-----------------------------------------------------------------------
        # Note: the following statement for determining the proper alignment of
        #       the offset text. This was determined entirely by trial-and-error
        #       and should not be in any way considered as "the way".  There are
        #       still some edge cases where alignment is not quite right, but
        #       this seems to be more of a geometry issue (in other words, I
        #       might be using the wrong reference points).
        #
        #   (TT, FF, TF, FT) are the shorthand for the tuple of
        #     (centpt[info['tickdir']] <= peparray[info['tickdir'], outerindex],
        #      centpt[index] <= peparray[index, outerindex])
        #
        #   Three-letters (e.g., TFT, FTT) are short-hand for the array
        #    of bools from the variable 'highs'.
        # ---------------------------------------------------------------------
        if centpt[info['tickdir']] > peparray[info['tickdir'], outerindex] :
            # if FT and if highs has an even number of Trues
            if (centpt[index] <= peparray[index, outerindex]
                and ((len(highs.nonzero()[0]) % 2) == 0)) :
                # Usually, this means align right, except for the FTT case,
                # in which offset for axis 1 and 2 are aligned left.
                if highs.tolist() == [False, True, True] and index in (1, 2) :
                    align = 'left'
                else :
                    align = 'right'
            else :
                # The FF case
                align = 'left'
        else :
            # if TF and if highs has an even number of Trues
            if (centpt[index] > peparray[index, outerindex]
                and ((len(highs.nonzero()[0]) % 2) == 0)) :
                # Usually mean align left, except if it is axis 2
                if index == 2 :
                    align = 'right'
                else :
                    align = 'left'
            else :
                # The TT case
                align = 'right'

        self.offsetText.set_va('center')
        self.offsetText.set_ha(align)
        self.offsetText.draw(renderer)

        # Draw grid lines
        if len(xyz0) > 0:
            # Grid points at end of one plane
            xyz1 = copy.deepcopy(xyz0)
            newindex = (index + 1) % 3
            newval = get_flip_min_max(xyz1[0], newindex, mins, maxs)
            for i in range(len(majorLocs)):
                xyz1[i][newindex] = newval

            # Grid points at end of the other plane
            xyz2 = copy.deepcopy(xyz0)
            newindex = (index + 2) %  3
            newval = get_flip_min_max(xyz2[0], newindex, mins, maxs)
            for i in range(len(majorLocs)):
                xyz2[i][newindex] = newval

            lines = list(zip(xyz1, xyz0, xyz2))
            if self.axes._draw_grid:
                self.gridlines.set_segments(lines)
                self.gridlines.set_color([info['grid']['color']] * len(lines))
                self.gridlines.draw(renderer, project=True)

        # Draw ticks
        tickdir = info['tickdir']
        tickdelta = deltas[tickdir]
        if highs[tickdir]:
            ticksign = 1
        else:
            ticksign = -1

        for tick, loc, label in zip(majorTicks, majorLocs, majorLabels):
            if tick is None:
                continue

            # Get tick line positions
            pos = copy.copy(edgep1)
            pos[index] = loc
            pos[tickdir] = edgep1[tickdir] + info['tick']['outward_factor'] * \
                                             ticksign * tickdelta
            x1, y1, z1 = proj3d.proj_transform(pos[0], pos[1], pos[2], \
                    renderer.M)
            pos[tickdir] = edgep1[tickdir] - info['tick']['inward_factor'] * \
                                             ticksign * tickdelta
            x2, y2, z2 = proj3d.proj_transform(pos[0], pos[1], pos[2], \
                    renderer.M)

            # Get position of label
            labeldeltas = [info['ticklabel']['space_factor'] * x for
                           x in deltas]
            axmask = [True, True, True]
            axmask[index] = False
            pos[tickdir] = edgep1[tickdir]
            pos = move_from_center(pos, centers, labeldeltas, axmask)
            lx, ly, lz = proj3d.proj_transform(pos[0], pos[1], pos[2], \
                    renderer.M)

            tick_update_position(tick, (x1, x2), (y1, y2), (lx, ly))
            tick.set_label1(label)
            tick.set_label2(label)
            tick.draw(renderer)

        renderer.close_group('axis3d')

Example 11

Project: pyNastran Source File: cart3d_io.py
    def _create_box(self, cart3d_filename, ID, form, cases, icase, regions):
        stack = True
        dirname = os.path.dirname(os.path.abspath(cart3d_filename))
        input_c3d_filename = os.path.join(dirname, 'input.c3d')
        input_cntl_filename = os.path.join(dirname, 'input.cntl')
        mach = None
        alpha = None
        beta = None
        if os.path.exists(input_cntl_filename):
            cntl = read_input_cntl(input_cntl_filename, log=self.log, debug=self.debug)
            mach, alpha, beta = cntl.get_flow_conditions()
            bcs = cntl.get_boundary_conditions()
            bc_xmin, bc_xmax, bc_ymin, bc_ymax, bc_xmin, bc_xmax, surfbcs = bcs
            #stack = False

            if surfbcs:
                bc_form = [
                    ('Rho', icase, []),
                    ('xVelocity', icase + 1, []),
                    ('yVelocity', icase + 2, []),
                    ('zVelocity', icase + 3, []),
                    ('Mach', icase + 4, []),
                    ('Pressure', icase + 5, []),
                ]
                icase += 5
                nelements = self.nElements
                rho = zeros(nelements, dtype='float32')
                xvel = zeros(nelements, dtype='float32')
                yvel = zeros(nelements, dtype='float32')
                zvel = zeros(nelements, dtype='float32')
                vel = zeros(nelements, dtype='float32')
                pressure = zeros(nelements, dtype='float32')

                uregions = set(unique(regions))
                surf_bc_regions = set(surfbcs.keys())
                invalid_regions = surf_bc_regions - uregions
                if len(invalid_regions) != 0:
                    assert len(invalid_regions) == 0, invalid_regions

                for bc_id, bc_values in sorted(iteritems(surfbcs)):
                    rhoi, xveli, yveli, zveli, pressi = bc_values
                    i = where(regions == bc_id)[0]
                    rho[i] = rhoi
                    xvel[i] = xveli
                    yvel[i] = yveli
                    zvel[i] = zveli
                    pressure[i] = pressi

                mach = sqrt(xvel ** 2 + yvel ** 2 + zvel ** 2)

                rho_res = GuiResult(ID, header='Rho', title='Rho',
                                    location='centroid', scalar=rho)
                xvel_res = GuiResult(ID, header='xVelocity', title='xVelocity',
                                     location='centroid', scalar=xvel)
                yvel_res = GuiResult(ID, header='yVelocity', title='yVelocity',
                                     location='centroid', scalar=yvel)
                zvel_res = GuiResult(ID, header='zVelocity', title='zVelocity',
                                     location='centroid', scalar=zvel)
                mach_res = GuiResult(ID, header='Mach', title='Mach',
                                     location='centroid', scalar=mach)
                pressure_res = GuiResult(ID, header='Pressure', title='Pressure',
                                         location='centroid', scalar=pressure)

                cases[icase] = (rho_res, (ID, 'Rho'))
                cases[icase + 1] = (xvel_res, (ID, 'xVelocity'))
                cases[icase + 2] = (yvel_res, (ID, 'yVelocity'))
                cases[icase + 3] = (zvel_res, (ID, 'zVelocity'))
                cases[icase + 4] = (mach_res, (ID, 'Mach'))
                cases[icase + 5] = (pressure_res, (ID, 'Pressure'))
                form.append(('Boundary Conditions', None, bc_form))


        if os.path.exists(input_c3d_filename):
            nodes, elements = read_input_c3d(input_c3d_filename, stack=stack, log=self.log, debug=self.debug)

            # Planes
            # ----------
            # xmin, xmax
            # ymin, ymax
            # zmin, zmax

            if stack:
                red = (1., 0., 0.)
                color = red
                self.set_quad_grid('box', nodes, elements, color, line_width=1, opacity=1.)
            else:
                red = (1., 0., 0.)
                inflow_nodes = []
                inflow_elements = []

                green = (0., 1., 0.)
                symmetry_nodes = []
                symmetry_elements = []

                colori = (1., 1., 0.)
                outflow_nodes = []
                outflow_elements = []

                blue = (0., 0., 1.)
                farfield_nodes = []
                farfield_elements = []

                ifarfield = 0
                isymmetry = 0
                iinflow = 0
                ioutflow = 0

                nfarfield_nodes = 0
                nsymmetry_nodes = 0
                ninflow_nodes = 0
                noutflow_nodes = 0
                for bcsi, nodesi, elementsi in zip(bcs, nodes, elements):
                    # 0 = FAR FIELD
                    # 1 = SYMMETRY
                    # 2 = INFLOW  (specify all)
                    # 3 = OUTFLOW (simple extrap)
                    self.log.info('bcsi = %s' % bcsi)
                    nnodes = nodesi.shape[0]
                    bc = bcsi
                    if isinstance(bc, int):
                        if bc == 0:
                            farfield_nodes.append(nodesi)
                            farfield_elements.append(elementsi + nfarfield_nodes)
                            nfarfield_nodes += nnodes
                            ifarfield += 1
                        elif bc == 1:
                            symmetry_nodes.append(nodesi)
                            symmetry_elements.append(elementsi + nsymmetry_nodes)
                            nsymmetry_nodes += nnodes
                            isymmetry += 1
                        elif bc == 2:
                            inflow_nodes.append(nodesi)
                            inflow_elements.append(elementsi + ninflow_nodes)
                            ninflow_nodes += nnodes
                            iinflow += 1
                        elif bc == 3:
                            outflow_nodes.append(nodesi)
                            outflow_elements.append(elementsi + noutflow_nodes)
                            noutflow_nodes += nnodes
                            ioutflow += 1
                        else:
                            msg = 'bc=%s' % str(bc)
                            raise NotImplementedError(msg)
                    elif isinstance(bc, dict):
                        continue
                    else:
                        msg = 'bc=%s' % str(bc)
                        raise NotImplementedError(msg)

                if ifarfield:
                    color = blue
                    nodes = vstack(farfield_nodes)
                    elements = vstack(farfield_elements)
                    self.set_quad_grid('farfield', nodes, elements, color, line_width=1, opacity=1.)

                if isymmetry:
                    color = green
                    nodes = vstack(symmetry_nodes)
                    elements = vstack(symmetry_elements)
                    self.set_quad_grid('symmetry', nodes, elements, color, line_width=1, opacity=1.)

                if iinflow:
                    color = red
                    nodes = vstack(inflow_nodes)
                    elements = vstack(inflow_elements)
                    self.set_quad_grid('inflow', nodes, elements, color, line_width=1, opacity=1.)

                if ioutflow:
                    color = colori
                    nodes = vstack(outflow_nodes)
                    elements = vstack(outflow_elements)
                    self.set_quad_grid('outflow', nodes, elements, color, line_width=1, opacity=1.)

                #i = 0
                #for nodesi, elementsi in zip(nodes, elements):
                    #self.set_quad_grid('box_%i' % i, nodesi, elementsi, color, line_width=1, opacity=1.)
                    #i += 1
        return mach, alpha, beta

Example 12

Project: pyhawkes Source File: preprocess.py
Function: process_dataset
def process_dataset(K=100,
                    suffix="_iNet1_Size100_CC01inh.txt",
                    dir="data/chalearn/small",
                    outfile="network1c.pkl"):

    # Get the full filenames
    fluor_file = os.path.join(dir, "fluorescence" + suffix)
    net_file   = os.path.join(dir, "network" + suffix)
    pos_file   = os.path.join(dir, "networkPositions" + suffix)

    # Parse the files
    F       = parse_fluorescence_file(fluor_file, K)

    # # Load the oopsi processed fluorescence
    # data_path = os.path.join("data", "chalearn", "small", "network1c.pkl.gz")
    # with gzip.open(data_path, 'r') as f:
    #     _, F, C, network, pos = cPickle.load(f)

    network = parse_network_file(net_file, K)
    pos     = parse_position_file(pos_file, K)

    # Discretize the fluorescence signal
    # Hardcode the bins
    # bins = np.array([-10, 0.17, 10]).reshape((1,3)).repeat(K, axis=0)
    # S, _ = discretize_fluorescence(F, edges=bins, binsui=False)
    # # S, bins = discretize_fluorescence(F, nbins=2, binsui=True)
    # S, bins = discretize_fluorescence(C, nbins=2, binsui=True)
    # # S = remove_double_spikes(S)

    # Get the spike times with oopsi
    # fast-oopsi,
    S,C = extract_spike_oopsi(F, dt=0.02)


    # Plot a segment of fluorescence traces and spikes
    start = 0
    end   = 10000
    k     = 0
    spks  = np.where(S[start:end, k])
    plt.figure()
    plt.plot(F[start:end, k], '-k')
    plt.plot(spks, F[spks,k], 'ro')
    plt.show()

    # Scatter plot the positions
    plt.figure()

    pres,posts = network.nonzero()
    for i,j in zip(pres,posts):
        if np.random.rand() < 0.25:
            plt.plot([pos[i,0], pos[j,0]],
                     [pos[i,1], pos[j,1]],
                     '-k', lw=0.5)
    plt.scatter(pos[:,0], pos[:,1], s=10, c='r', marker='o', facecolor='k')
    plt.xlabel("X")
    plt.ylabel("Y")
    plt.show()

    # Plot the network as a function of X position
    perm = np.argsort(pos[:,0])
    plt.figure()
    plt.imshow(network[np.ix_(perm, perm)], cmap="Greys", interpolation="none")
    plt.xlabel("Postsynaptic")
    plt.ylabel("Presynaptic")
    plt.show()

    with gzip.open(os.path.join(dir, outfile + ".gz"), 'w') as f:
        cPickle.dump((S, F, C, network, pos), f, protocol=-1)

Example 13

Project: mlxtend Source File: adaline.py
Function: fit
    def _fit(self, X, y, init_params=True):
        self._check_target_array(y, allowed={(0, 1)})
        y_data = np.where(y == 0, -1., 1.)

        if init_params:
            self.b_, self.w_ = self._init_params(
                weights_shape=(X.shape[1], 1),
                bias_shape=(1,),
                random_seed=self.random_seed)
            self.cost_ = []

        if self.minibatches is None:
            self.b_, self.w_ = self._normal_equation(X, y_data)

        # Gradient descent or stochastic gradient descent learning
        else:
            self.init_time_ = time()
            rgen = np.random.RandomState(self.random_seed)
            for i in range(self.epochs):

                for idx in self._yield_minibatches_idx(
                        rgen=rgen,
                        n_batches=self.minibatches,
                        data_ary=y_data,
                        shuffle=True):

                    y_val = self._net_input(X[idx])
                    errors = (y_data[idx] - y_val)
                    self.w_ += (self.eta *
                                X[idx].T.dot(errors).reshape(self.w_.shape))
                    self.b_ += self.eta * errors.sum()

                cost = self._sum_squared_error_cost(y_data, self._net_input(X))
                self.cost_.append(cost)
                if self.print_progress:
                    self._print_progress(iteration=(i + 1),
                                         n_iter=self.epochs,
                                         cost=cost)

        return self

Example 14

Project: spectralDNS Source File: MKM.py
def initialize(solver, context):
    # Initialize with pertubation ala perturbU (https://github.com/wyldckat/perturbU) for openfoam
    U = context.U
    X = context.X
    params = config.params
    
    Y = where(X[0]<0, 1+X[0], 1-X[0])
    utau = params.nu * params.Re_tau
    #Um = 46.9091*utau
    Um = 56.*utau 
    Xplus = Y*params.Re_tau
    Yplus = X[1]*params.Re_tau
    Zplus = X[2]*params.Re_tau
    duplus = Um*0.2/utau  #Um*0.25/utau 
    alfaplus = params.L[1]/500.
    betaplus = params.L[2]/200.
    sigma = 0.00055 # 0.00055
    epsilon = Um/200.   #Um/200.
    U[:] = 0
    U[1] = Um*(Y-0.5*Y**2)
    dev = 1+0.005*random.randn(Y.shape[0], Y.shape[1], Y.shape[2])
    dd = utau*duplus/2.0*Xplus/40.*exp(-sigma*Xplus**2+0.5)*cos(betaplus*Zplus)*dev
    U[1] += dd
    U[2] += epsilon*sin(alfaplus*Yplus)*Xplus*exp(-sigma*Xplus**2)*dev    
    
    U_hat = solver.set_velocity(**context)
    U = solver.get_velocity(**context)
    U_hat = solver.set_velocity(**context)
    
    if "KMM" in params.solver:
        context.g[:] = 1j*context.K[1]*U_hat[2] - 1j*context.K[2]*U_hat[1]

    # Set the flux
    
    #flux[0] = context.FST.dx(U[1], context.ST.quad)
    #solver.comm.Bcast(flux)
    
    if solver.rank == 0:
        print("Flux {}".format(flux[0]))
    
    if not params.solver in ("KMM", "KMMRK3"):
        P_hat = solver.compute_pressure(**context)
        P = context.FST.ifst(P_hat, context.P, context.SN)
        
    context.U_hat0[:] = context.U_hat[:]
    context.H_hat1[:] = solver.get_convection(**context)

Example 15

Project: zipline Source File: test_factor.py
    @parameter_space(
        seed_value=range(1, 2),
        normalizer_name_and_func=[
            ('demean', lambda row: row - nanmean(row)),
            ('zscore', lambda row: (row - nanmean(row)) / nanstd(row)),
        ],
        add_nulls_to_factor=(False, True,),
    )
    def test_normalizations_randomized(self,
                                       seed_value,
                                       normalizer_name_and_func,
                                       add_nulls_to_factor):

        name, func = normalizer_name_and_func

        shape = (7, 7)

        # All Trues.
        nomask = self.ones_mask(shape=shape)
        # Falses on main diagonal.
        eyemask = self.eye_mask(shape=shape)
        # Falses on other diagonal.
        eyemask90 = rot90(eyemask)
        # Falses on both diagonals.
        xmask = eyemask & eyemask90

        # Block of random data.
        factor_data = self.randn_data(seed=seed_value, shape=shape)
        if add_nulls_to_factor:
            factor_data = where(eyemask, factor_data, nan)

        # Cycles of 0, 1, 2, 0, 1, 2, ...
        classifier_data = (
            (self.arange_data(shape=shape, dtype=int64_dtype) + seed_value) % 3
        )
        # With -1s on main diagonal.
        classifier_data_eyenulls = where(eyemask, classifier_data, -1)
        # With -1s on opposite diagonal.
        classifier_data_eyenulls90 = where(eyemask90, classifier_data, -1)
        # With -1s on both diagonals.
        classifier_data_xnulls = where(xmask, classifier_data, -1)

        f = self.f
        c = C()
        c_with_nulls = OtherC()
        m = Mask()
        method = getattr(f, name)
        terms = {
            'vanilla': method(),
            'masked': method(mask=m),
            'grouped': method(groupby=c),
            'grouped_with_nulls': method(groupby=c_with_nulls),
            'both': method(mask=m, groupby=c),
            'both_with_nulls': method(mask=m, groupby=c_with_nulls),
        }

        expected = {
            'vanilla': apply_along_axis(func, 1, factor_data,),
            'masked': where(
                eyemask,
                grouped_apply(factor_data, eyemask, func),
                nan,
            ),
            'grouped': grouped_apply(
                factor_data,
                classifier_data,
                func,
            ),
            # If the classifier has nulls, we should get NaNs in the
            # corresponding locations in the output.
            'grouped_with_nulls': where(
                eyemask90,
                grouped_apply(factor_data, classifier_data_eyenulls90, func),
                nan,
            ),
            # Passing a mask with a classifier should behave as though the
            # classifier had nulls where the mask was False.
            'both': where(
                eyemask,
                grouped_apply(
                    factor_data,
                    classifier_data_eyenulls,
                    func,
                ),
                nan,
            ),
            'both_with_nulls': where(
                xmask,
                grouped_apply(
                    factor_data,
                    classifier_data_xnulls,
                    func,
                ),
                nan,
            )
        }

        self.check_terms(
            terms=terms,
            expected=expected,
            initial_workspace={
                f: factor_data,
                c: classifier_data,
                c_with_nulls: classifier_data_eyenulls90,
                Mask(): eyemask,
            },
            mask=self.build_mask(nomask),
        )

Example 16

Project: pyNastran Source File: write_mesh.py
    def _write_elements_properties(self, bdf_file, size):
        """
        Writes the elements and properties in and interspersed order
        """
        #return self._write_elements_properties2(f, size)
        msg = []
        missing_properties = []
        if self.properties:
            bdf_file.write('$ELEMENTS_WITH_PROPERTIES\n')

        eids_written = []
        #pids = sorted(self.properties.keys())

        ptypes = [
            self.properties_shell.pshell,
            self.properties_shell.pcomp,
            self.pshear,
            self.prod,
            self.properties_solid.psolid,

            #self.properties_bar.pbar,
            #self.properties_bar.pbarl,
            #self.properties_beam.pbeam,
            #self.properties_beam.pbeaml,
        ]

        n = 0
        pids_all = None  # the actual properties
        for t in ptypes:
            if t.n and n == 0:
                pids_all = t.property_id
                n = 1
            elif t.n:
                self.log.debug(pids_all)
                self.log.debug(t.property_id)
                try:
                    pids_all = concatenate(pids_all, t.property_id)
                except ValueError:
                    pids_all = array(list(pids_all) + list(t.property_id))

        etypes = (self.elements_shell._get_types() +
                  self.elements_solid._get_types() +
                  [self.crod, self.cshear])

        #pids_set = None
        if pids_all is None:
            bdf_file.write('$MISSING_ELEMENTS because there are no properties\n')
            for t in etypes:
                #print("t.type =", t.type)
                t.write_card(bdf_file, size=size)
            return

        # there are properties
        pids_set = set(list(pids_all))

        n = 0
        pids = None
        for t in etypes:
            #print("t.type =", t.type)
            if t.n and n == 0:
                eids = t.element_id
                pids = t.property_id
                n = 1
            elif t.n:
                try:
                    eids = concatenate(eids, t.element_id)
                #except AttributeError:
                    #eids = array(list(eids) + list(t.element_id))
                except TypeError:
                    #print eids
                    #print t.element_id
                    eids = array(list(eids) + list(t.element_id))
                except ValueError:
                    #print eids
                    #print t.element_id
                    eids = array(list(eids) + list(t.element_id))

                try:
                    pids = concatenate(pids, t.property_id)
                except AttributeError:
                    pids = array(list(pids) + list(t.property_id))
                except TypeError:
                    pids = array(list(pids) + list(t.property_id))
                except ValueError:
                    pids = array(list(pids) + list(t.property_id))
            #else:
                #print t.type

        elements_by_pid = {}
        if pids is not None:
            pids_unique = unique(pids)
            self.log.debug("pids_unique = %s" % pids_unique)
            pids_unique.sort()
            if len(pids_unique) > 0:
                bdf_file.write('$ELEMENTS_WITH_PROPERTIES\n')

            for pid in pids_all:
                i = where(pid==pids)[0]
                eids2 = eids[i]

                for t in ptypes:
                    if t.n and pid in t.property_id:
                        self.log.debug("prop.type = %s" % t.type)
                        t.write_card(bdf_file, size=size, property_ids=[pid])
                        pids_set.remove(pid)
                n = 0
                for t in etypes:
                    if not t.n:
                        continue
                    eids3 = intersect1d(t.element_id, eids2, assume_unique=False)
                    #print("eids3[pid=%s]" %(pid), eids3)
                    if n == 0 and len(eids3):
                        elements_by_pid[pid] = eids3
                        n = 1
                    elif len(eids3):
                        try:
                            c = concatenate(elements_by_pid[pid], eids3)
                        except TypeError:
                            c = array(list(elements_by_pid[pid]) + list(eids3))
                        except ValueError:
                            c = array(list(elements_by_pid[pid]) + list(eids3))
                        elements_by_pid[pid] = c
                    else:
                        continue
                    try:
                        t.write_card(bdf_file, size=size, element_ids=eids3)
                    except TypeError:
                        print("t.type = %s" % t.type)
                        raise
                    del eids3
            #for pid, elements in elements_by_pid.items():
                #print("pid=%s n=%s" % (pid, len(elements)))
            #print elements_by_pid

        # missing properties
        if pids_set:
            pids_list = list(pids_set)
            bdf_file.write('$UNASSOCIATED_PROPERTIES\n')
            for pid in pids_list:
                for t in ptypes:
                    if t.n and pid in t.property_id:
                        t.write_card(bdf_file, size=size, property_ids=[pid])

        #.. todo:: finish...
        bdf_file.write('$UNASSOCIATED_ELEMENTS\n')

Example 17

Project: scikit-beam Source File: roi.py
def circular_average(image, calibrated_center, threshold=0, nx=100,
                     pixel_size=(1, 1), min_x=None, max_x=None, mask=None):
    """Circular average of the the image data
    The circular average is also known as the radial integration

    Parameters
    ----------
    image : array
        Image to compute the average as a function of radius
    calibrated_center : tuple
        The center of the image in pixel units
        argument order should be (row, col)
    threshold : int, optional
        Ignore counts below `threshold`
        default is zero
    nx : int, optional
        number of bins in x
        defaults is 100 bins
    pixel_size : tuple, optional
        The size of a pixel (in a real unit, like mm).
        argument order should be (pixel_height, pixel_width)
        default is (1, 1)
    min_x : float, optional number of pixels
        Left edge of first bin defaults to minimum value of x
    max_x : float, optional number of pixels
        Right edge of last bin defaults to maximum value of x
    mask : mask for 2D data. Assumes 1 is non masked and 0 masked.
        None defaults to no mask.

    Returns
    -------
    bin_centers : array
        The center of each bin in R. shape is (nx, )
    ring_averages : array
        Radial average of the image. shape is (nx, ).

    See Also
    --------
    bad_to_nan_gen : Create a mask with np.nan entries
    bin_grid : Bin and integrate an image, given the radial array of pixels
        Useful for nonlinear spacing (Ewald curvature)
    """
    radial_val = utils.radial_grid(calibrated_center, image.shape, pixel_size)

    if mask is not None:
        w = np.where(mask == 1)
        radial_val = radial_val[w]
        image = image[w]

    bin_edges, sums, counts = utils.bin_1D(np.ravel(radial_val),
                                           np.ravel(image), nx,
                                           min_x=min_x,
                                           max_x=max_x)
    th_mask = counts > threshold
    ring_averages = sums[th_mask] / counts[th_mask]

    bin_centers = utils.bin_edges_to_centers(bin_edges)[th_mask]

    return bin_centers, ring_averages

Example 18

Project: ComputationalNeurodynamics Source File: IzNetwork.py
Function: neuronupdate
  def NeuronUpdate(self, i, t):
    """
    Izhikevich neuron update function. Update one layer for 1 millisecond
    using the Euler method.

    Inputs:
    i -- Number of layer to update
    t -- Current timestep. Necessary to sort out the synaptic delays.
    """

    # Euler method step size in ms
    dt = 0.2

    # Calculate current from incoming spikes
    for j in xrange(self.Nlayers):

      # If layer[i].S[j] exists then layer[i].factor[j] and
      # layer[i].delay[j] have to exist
      if j in self.layer[i].S:
        S = self.layer[i].S[j]  # target neuron->rows, source neuron->columns

        # Firings contains time and neuron idx of each spike.
        # [t, index of the neuron in the layer j]
        firings = self.layer[j].firings

        # Find incoming spikes taking delays into account
        delay = self.layer[i].delay[j]
        F = self.layer[i].factor[j]

        # Sum current from incoming spikes
        k = len(firings)
        while k > 0 and (firings[k-1, 0] > (t - self.Dmax)):
          idx = delay[:, firings[k-1, 1]] == (t-firings[k-1, 0])
          self.layer[i].I[idx] += F * S[idx, firings[k-1, 1]]
          k = k-1

    # Update v and u using the Izhikevich model and Euler method
    for k in xrange(int(1/dt)):
      v = self.layer[i].v
      u = self.layer[i].u

      self.layer[i].v += dt*(0.04*v*v + 5*v + 140 - u + self.layer[i].I)
      self.layer[i].u += dt*(self.layer[i].a*(self.layer[i].b*v - u))

      # Find index of neurons that have fired this millisecond
      fired = np.where(self.layer[i].v >= 30)[0]

      if len(fired) > 0:
        for f in fired:
          # Add spikes into spike train
          if len(self.layer[i].firings) != 0:
            self.layer[i].firings = np.vstack([self.layer[i].firings, [t, f]])
          else:
            self.layer[i].firings = np.array([[t, f]])

          # Reset the membrane potential after spikes
          self.layer[i].v[f]  = self.layer[i].c[f]
          self.layer[i].u[f] += self.layer[i].d[f]

    return

Example 19

Project: scikit-learn Source File: svmlight_format.py
def _dump_svmlight(X, y, f, multilabel, one_based, comment, query_id):
    X_is_sp = int(hasattr(X, "tocsr"))
    y_is_sp = int(hasattr(y, "tocsr"))
    if X.dtype.kind == 'i':
        value_pattern = u("%d:%d")
    else:
        value_pattern = u("%d:%.16g")

    if y.dtype.kind == 'i':
        label_pattern = u("%d")
    else:
        label_pattern = u("%.16g")

    line_pattern = u("%s")
    if query_id is not None:
        line_pattern += u(" qid:%d")
    line_pattern += u(" %s\n")

    if comment:
        f.write(b("# Generated by dump_svmlight_file from scikit-learn %s\n"
                % __version__))
        f.write(b("# Column indices are %s-based\n"
                  % ["zero", "one"][one_based]))

        f.write(b("#\n"))
        f.writelines(b("# %s\n" % line) for line in comment.splitlines())

    for i in range(X.shape[0]):
        if X_is_sp:
            span = slice(X.indptr[i], X.indptr[i + 1])
            row = zip(X.indices[span], X.data[span])
        else:
            nz = X[i] != 0
            row = zip(np.where(nz)[0], X[i, nz])

        s = " ".join(value_pattern % (j + one_based, x) for j, x in row)

        if multilabel:
            if y_is_sp:
                nz_labels = y[i].nonzero()[1]
            else:
                nz_labels = np.where(y[i] != 0)[0]
            labels_str = ",".join(label_pattern % j for j in nz_labels)
        else:
            if y_is_sp:
                labels_str = label_pattern % y.data[i]
            else:
                labels_str = label_pattern % y[i]

        if query_id is not None:
            feat = (labels_str, query_id[i], s)
        else:
            feat = (labels_str, s)

        f.write((line_pattern % feat).encode('ascii'))

Example 20

Project: pyspeckit Source File: collapse_gaussfit.py
def wrap_collapse_gauss(filename,outprefix,redo='no'):
    """

    redo - if not equal to 'no', then...
    if collapse_gaussfit succeeded (to the extent that the .pysav files were written),
    but some part of the file writing or successive procedures failed, re-do those 
    procedures without redoing the whole collapse
    """
    fitsfile = pyfits.open(filename)
    dv,v0,p3 = fitsfile[0].header['CD3_3'],fitsfile[0].header['CRVAL3'],fitsfile[0].header['CRPIX3']

    cube = fitsfile[0].data
    cube = where(numpy.isnan(cube),0,cube)

    if redo=='no':
        doubleB = asarray(collapse_double_gaussfit(cube,axis=0))
        doubleB[numpy.isnan(doubleB)] = 0
        pickle.dump(doubleB,open('%s_doubleB.pysav' % outprefix,'w'))
    else:
        doubleB = pickle.load(open('%s_doubleB.pysav' % outprefix,'r'))
    db = doubleB
    gcd = double_gaussian(db[3],db[4],db[0],db[1],db[5],db[6])(indices(cube.shape)[0])
    fitsfile[0].data = gcd
    fitsfile.writeto('%s_doublegausscube.fits' % outprefix,clobber=True)

    gcd[numpy.isnan(gcd)] = 0
    doubleResids = cube-gcd
    fitsfile[0].data = doubleResids
    fitsfile.writeto('%s_doublegaussresids.fits' % outprefix,clobber=True)


    #doubleB[4] = (doubleB[4]-v0) / dv + p3-1
    #doubleB[3] = (doubleB[3]-v0) / dv + p3-1
    doubleB[4] = (doubleB[4]-p3+1) * dv + v0
    doubleB[3] = (doubleB[3]-p3+1) * dv + v0
    fitsfile[0].data = asarray(doubleB)
    fitsfile.writeto('%s_doublegausspars.fits' % outprefix,clobber=True)

    if redo=='no':
        singleB = asarray(collapse_gaussfit(cube,axis=0))
        pickle.dump(singleB,open('%s_singleB.pysav' % outprefix,'w'))
    else:
        singleB = pickle.load(open('%s_singleB.pysav' % outprefix,'r'))
    gc = gaussian(singleB[1],singleB[0],singleB[2])(indices(cube.shape)[0])
    singleB[1] = (singleB[1]-p3+1) * dv + v0
    fitsfile[0].data = gc
    fitsfile.writeto('%s_singlegausscube.fits' % outprefix,clobber=True)

    gc[numpy.isnan(gc)]=0
    singleResids = cube-gc
    fitsfile[0].data = singleResids
    fitsfile.writeto('%s_singlegaussresids.fits' % outprefix,clobber=True)
    fitsfile[0].data = asarray(singleB)
    fitsfile.writeto('%s_singlegausspars.fits' % outprefix,clobber=True)

    fitsfile[0].header.__delitem__('CD3_3')
    fitsfile[0].header.__delitem__('CRVAL3')
    fitsfile[0].header.__delitem__('CRPIX3')
    fitsfile[0].header.__delitem__('CUNIT3')
    fitsfile[0].header.__delitem__('CTYPE3')

    doubleResids[numpy.isnan(doubleResids)] = 0
    totalDResids = doubleResids.sum(axis=0)
    fitsfile[0].data = totalDResids
    fitsfile.writeto('%s_doublegauss_totalresids.fits' % outprefix,clobber=True)

    singleResids[numpy.isnan(singleResids)] = 0
    totalSResids = singleResids.sum(axis=0)
    fitsfile[0].data = totalSResids
    fitsfile.writeto('%s_singlegauss_totalresids.fits' % outprefix,clobber=True)

    return singleB,doubleB

Example 21

Project: zipline Source File: test_statistical.py
    @parameter_space(returns_length=[2, 3], correlation_length=[3, 4])
    def test_correlation_factors(self, returns_length, correlation_length):
        """
        Tests for the built-in factors `RollingPearsonOfReturns` and
        `RollingSpearmanOfReturns`.
        """
        assets = self.assets
        my_asset = self.my_asset
        my_asset_column = self.my_asset_column
        dates = self.dates
        start_date = self.pipeline_start_date
        end_date = self.pipeline_end_date
        start_date_index = self.start_date_index
        end_date_index = self.end_date_index
        num_days = self.num_days
        run_pipeline = self.run_pipeline

        returns = Returns(window_length=returns_length)
        masks = (self.cascading_mask, self.alternating_mask, NotSpecified)
        expected_mask_results = (
            self.expected_cascading_mask_result,
            self.expected_alternating_mask_result,
            self.expected_no_mask_result,
        )

        for mask, expected_mask in zip(masks, expected_mask_results):
            pearson_factor = RollingPearsonOfReturns(
                target=my_asset,
                returns_length=returns_length,
                correlation_length=correlation_length,
                mask=mask,
            )
            spearman_factor = RollingSpearmanOfReturns(
                target=my_asset,
                returns_length=returns_length,
                correlation_length=correlation_length,
                mask=mask,
            )

            columns = {
                'pearson_factor': pearson_factor,
                'spearman_factor': spearman_factor,
            }
            pipeline = Pipeline(columns=columns)
            if mask is not NotSpecified:
                pipeline.add(mask, 'mask')

            results = run_pipeline(pipeline, start_date, end_date)
            pearson_results = results['pearson_factor'].unstack()
            spearman_results = results['spearman_factor'].unstack()
            if mask is not NotSpecified:
                mask_results = results['mask'].unstack()
                check_arrays(mask_results.values, expected_mask)

            # Run a separate pipeline that calculates returns starting
            # (correlation_length - 1) days prior to our start date. This is
            # because we need (correlation_length - 1) extra days of returns to
            # compute our expected correlations.
            results = run_pipeline(
                Pipeline(columns={'returns': returns}),
                dates[start_date_index - (correlation_length - 1)],
                dates[end_date_index],
            )
            returns_results = results['returns'].unstack()

            # On each day, calculate the expected correlation coefficients
            # between the asset we are interested in and each other asset. Each
            # correlation is calculated over `correlation_length` days.
            expected_pearson_results = full_like(pearson_results, nan)
            expected_spearman_results = full_like(spearman_results, nan)
            for day in range(num_days):
                todays_returns = returns_results.iloc[
                    day:day + correlation_length
                ]
                my_asset_returns = todays_returns.iloc[:, my_asset_column]
                for asset, other_asset_returns in todays_returns.iteritems():
                    asset_column = int(asset) - 1
                    expected_pearson_results[day, asset_column] = pearsonr(
                        my_asset_returns, other_asset_returns,
                    )[0]
                    expected_spearman_results[day, asset_column] = spearmanr(
                        my_asset_returns, other_asset_returns,
                    )[0]

            expected_pearson_results = DataFrame(
                data=where(expected_mask, expected_pearson_results, nan),
                index=dates[start_date_index:end_date_index + 1],
                columns=assets,
            )
            assert_frame_equal(pearson_results, expected_pearson_results)

            expected_spearman_results = DataFrame(
                data=where(expected_mask, expected_spearman_results, nan),
                index=dates[start_date_index:end_date_index + 1],
                columns=assets,
            )
            assert_frame_equal(spearman_results, expected_spearman_results)

Example 22

Project: scikit-learn Source File: base.py
Function: radius_neighbors
    def radius_neighbors(self, X=None, radius=None, return_distance=True):
        """Finds the neighbors within a given radius of a point or points.

        Return the indices and distances of each point from the dataset
        lying in a ball with size ``radius`` around the points of the query
        array. Points lying on the boundary are included in the results.

        The result points are *not* necessarily sorted by distance to their
        query point.

        Parameters
        ----------
        X : array-like, (n_samples, n_features), optional
            The query point or points.
            If not provided, neighbors of each indexed point are returned.
            In this case, the query point is not considered its own neighbor.

        radius : float
            Limiting distance of neighbors to return.
            (default is the value passed to the constructor).

        return_distance : boolean, optional. Defaults to True.
            If False, distances will not be returned

        Returns
        -------
        dist : array, shape (n_samples,) of arrays
            Array representing the distances to each point, only present if
            return_distance=True. The distance values are computed according
            to the ``metric`` constructor parameter.

        ind : array, shape (n_samples,) of arrays
            An array of arrays of indices of the approximate nearest points
            from the population matrix that lie within a ball of size
            ``radius`` around the query points.

        Examples
        --------
        In the following example, we construct a NeighborsClassifier
        class from an array representing our data set and ask who's
        the closest point to [1, 1, 1]:

        >>> import numpy as np
        >>> samples = [[0., 0., 0.], [0., .5, 0.], [1., 1., .5]]
        >>> from sklearn.neighbors import NearestNeighbors
        >>> neigh = NearestNeighbors(radius=1.6)
        >>> neigh.fit(samples) # doctest: +ELLIPSIS
        NearestNeighbors(algorithm='auto', leaf_size=30, ...)
        >>> rng = neigh.radius_neighbors([[1., 1., 1.]])
        >>> print(np.asarray(rng[0][0])) # doctest: +ELLIPSIS
        [ 1.5  0.5]
        >>> print(np.asarray(rng[1][0])) # doctest: +ELLIPSIS
        [1 2]

        The first array returned contains the distances to all points which
        are closer than 1.6, while the second array returned contains their
        indices.  In general, multiple points can be queried at the same time.

        Notes
        -----
        Because the number of neighbors of each point is not necessarily
        equal, the results for multiple query points cannot be fit in a
        standard data array.
        For efficiency, `radius_neighbors` returns arrays of objects, where
        each object is a 1D array of indices or distances.
        """
        if self._fit_method is None:
            raise NotFittedError("Must fit neighbors before querying.")

        if X is not None:
            query_is_train = False
            X = check_array(X, accept_sparse='csr')
        else:
            query_is_train = True
            X = self._fit_X

        if radius is None:
            radius = self.radius

        n_samples = X.shape[0]
        if self._fit_method == 'brute':
            # for efficiency, use squared euclidean distances
            if self.effective_metric_ == 'euclidean':
                dist = pairwise_distances(X, self._fit_X, 'euclidean',
                                          n_jobs=self.n_jobs, squared=True)
                radius *= radius
            else:
                dist = pairwise_distances(X, self._fit_X,
                                          self.effective_metric_,
                                          n_jobs=self.n_jobs,
                                          **self.effective_metric_params_)

            neigh_ind_list = [np.where(d <= radius)[0] for d in dist]

            # See https://github.com/numpy/numpy/issues/5456
            # if you want to understand why this is initialized this way.
            neigh_ind = np.empty(n_samples, dtype='object')
            neigh_ind[:] = neigh_ind_list

            if return_distance:
                dist_array = np.empty(n_samples, dtype='object')
                if self.effective_metric_ == 'euclidean':
                    dist_list = [np.sqrt(d[neigh_ind[i]])
                                 for i, d in enumerate(dist)]
                else:
                    dist_list = [d[neigh_ind[i]]
                                 for i, d in enumerate(dist)]
                dist_array[:] = dist_list

                results = dist_array, neigh_ind
            else:
                results = neigh_ind

        elif self._fit_method in ['ball_tree', 'kd_tree']:
            if issparse(X):
                raise ValueError(
                    "%s does not work with sparse matrices. Densify the data, "
                    "or set algorithm='brute'" % self._fit_method)
            results = self._tree.query_radius(X, radius,
                                              return_distance=return_distance)
            if return_distance:
                results = results[::-1]
        else:
            raise ValueError("internal: _fit_method not recognized")

        if not query_is_train:
            return results
        else:
            # If the query data is the same as the indexed data, we would like
            # to ignore the first nearest neighbor of every sample, i.e
            # the sample itself.
            if return_distance:
                dist, neigh_ind = results
            else:
                neigh_ind = results

            for ind, ind_neighbor in enumerate(neigh_ind):
                mask = ind_neighbor != ind

                neigh_ind[ind] = ind_neighbor[mask]
                if return_distance:
                    dist[ind] = dist[ind][mask]

            if return_distance:
                return dist, neigh_ind
            return neigh_ind

Example 23

Project: scipy Source File: slsqp.py
def _minimize_slsqp(func, x0, args=(), jac=None, bounds=None,
                    constraints=(),
                    maxiter=100, ftol=1.0E-6, iprint=1, disp=False,
                    eps=_epsilon, callback=None,
                    **unknown_options):
    """
    Minimize a scalar function of one or more variables using Sequential
    Least SQuares Programming (SLSQP).

    Options
    -------
    ftol : float
        Precision goal for the value of f in the stopping criterion.
    eps : float
        Step size used for numerical approximation of the jacobian.
    disp : bool
        Set to True to print convergence messages. If False,
        `verbosity` is ignored and set to 0.
    maxiter : int
        Maximum number of iterations.

    """
    _check_unknown_options(unknown_options)
    fprime = jac
    iter = maxiter
    acc = ftol
    epsilon = eps

    if not disp:
        iprint = 0

    # Constraints are triaged per type into a dictionnary of tuples
    if isinstance(constraints, dict):
        constraints = (constraints, )

    cons = {'eq': (), 'ineq': ()}
    for ic, con in enumerate(constraints):
        # check type
        try:
            ctype = con['type'].lower()
        except KeyError:
            raise KeyError('Constraint %d has no type defined.' % ic)
        except TypeError:
            raise TypeError('Constraints must be defined using a '
                            'dictionary.')
        except AttributeError:
            raise TypeError("Constraint's type must be a string.")
        else:
            if ctype not in ['eq', 'ineq']:
                raise ValueError("Unknown constraint type '%s'." % con['type'])

        # check function
        if 'fun' not in con:
            raise ValueError('Constraint %d has no function defined.' % ic)

        # check jacobian
        cjac = con.get('jac')
        if cjac is None:
            # approximate jacobian function.  The factory function is needed
            # to keep a reference to `fun`, see gh-4240.
            def cjac_factory(fun):
                def cjac(x, *args):
                    return approx_jacobian(x, fun, epsilon, *args)
                return cjac
            cjac = cjac_factory(con['fun'])

        # update constraints' dictionary
        cons[ctype] += ({'fun': con['fun'],
                         'jac': cjac,
                         'args': con.get('args', ())}, )

    exit_modes = {-1: "Gradient evaluation required (g & a)",
                    0: "Optimization terminated successfully.",
                    1: "Function evaluation required (f & c)",
                    2: "More equality constraints than independent variables",
                    3: "More than 3*n iterations in LSQ subproblem",
                    4: "Inequality constraints incompatible",
                    5: "Singular matrix E in LSQ subproblem",
                    6: "Singular matrix C in LSQ subproblem",
                    7: "Rank-deficient equality constraint subproblem HFTI",
                    8: "Positive directional derivative for linesearch",
                    9: "Iteration limit exceeded"}

    # Wrap func
    feval, func = wrap_function(func, args)

    # Wrap fprime, if provided, or approx_jacobian if not
    if fprime:
        geval, fprime = wrap_function(fprime, args)
    else:
        geval, fprime = wrap_function(approx_jacobian, (func, epsilon))

    # Transform x0 into an array.
    x = asfarray(x0).flatten()

    # Set the parameters that SLSQP will need
    # meq, mieq: number of equality and inequality constraints
    meq = sum(map(len, [atleast_1d(c['fun'](x, *c['args'])) for c in cons['eq']]))
    mieq = sum(map(len, [atleast_1d(c['fun'](x, *c['args'])) for c in cons['ineq']]))
    # m = The total number of constraints
    m = meq + mieq
    # la = The number of constraints, or 1 if there are no constraints
    la = array([1, m]).max()
    # n = The number of independent variables
    n = len(x)

    # Define the workspaces for SLSQP
    n1 = n + 1
    mineq = m - meq + n1 + n1
    len_w = (3*n1+m)*(n1+1)+(n1-meq+1)*(mineq+2) + 2*mineq+(n1+mineq)*(n1-meq) \
            + 2*meq + n1 + ((n+1)*n)//2 + 2*m + 3*n + 3*n1 + 1
    len_jw = mineq
    w = zeros(len_w)
    jw = zeros(len_jw)

    # Decompose bounds into xl and xu
    if bounds is None or len(bounds) == 0:
        xl = np.empty(n, dtype=float)
        xu = np.empty(n, dtype=float)
        xl.fill(np.nan)
        xu.fill(np.nan)
    else:
        bnds = array(bounds, float)
        if bnds.shape[0] != n:
            raise IndexError('SLSQP Error: the length of bounds is not '
                             'compatible with that of x0.')

        bnderr = where(bnds[:, 0] > bnds[:, 1])[0]
        if bnderr.any():
            raise ValueError('SLSQP Error: lb > ub in bounds %s.' %
                             ', '.join(str(b) for b in bnderr))
        xl, xu = bnds[:, 0], bnds[:, 1]

        # Mark infinite bounds with nans; the Fortran code understands this
        infbnd = ~isfinite(bnds)
        xl[infbnd[:, 0]] = np.nan
        xu[infbnd[:, 1]] = np.nan

    # Initialize the iteration counter and the mode value
    mode = array(0,int)
    acc = array(acc,float)
    majiter = array(iter,int)
    majiter_prev = 0

    # Print the header if iprint >= 2
    if iprint >= 2:
        print("%5s %5s %16s %16s" % ("NIT","FC","OBJFUN","GNORM"))

    while 1:

        if mode == 0 or mode == 1:  # objective and constraint evaluation requird

            # Compute objective function
            try:
                fx = float(np.asarray(func(x)))
            except:
                raise ValueError("Objective function must return a scalar")
            # Compute the constraints
            if cons['eq']:
                c_eq = concatenate([atleast_1d(con['fun'](x, *con['args']))
                                     for con in cons['eq']])
            else:
                c_eq = zeros(0)
            if cons['ineq']:
                c_ieq = concatenate([atleast_1d(con['fun'](x, *con['args']))
                                     for con in cons['ineq']])
            else:
                c_ieq = zeros(0)

            # Now combine c_eq and c_ieq into a single matrix
            c = concatenate((c_eq, c_ieq))

        if mode == 0 or mode == -1:  # gradient evaluation required

            # Compute the derivatives of the objective function
            # For some reason SLSQP wants g dimensioned to n+1
            g = append(fprime(x),0.0)

            # Compute the normals of the constraints
            if cons['eq']:
                a_eq = vstack([con['jac'](x, *con['args'])
                               for con in cons['eq']])
            else:  # no equality constraint
                a_eq = zeros((meq, n))

            if cons['ineq']:
                a_ieq = vstack([con['jac'](x, *con['args'])
                                for con in cons['ineq']])
            else:  # no inequality constraint
                a_ieq = zeros((mieq, n))

            # Now combine a_eq and a_ieq into a single a matrix
            if m == 0:  # no constraints
                a = zeros((la, n))
            else:
                a = vstack((a_eq, a_ieq))
            a = concatenate((a,zeros([la,1])),1)

        # Call SLSQP
        slsqp(m, meq, x, xl, xu, fx, c, g, a, acc, majiter, mode, w, jw)

        # call callback if major iteration has incremented
        if callback is not None and majiter > majiter_prev:
            callback(x)

        # Print the status of the current iterate if iprint > 2 and the
        # major iteration has incremented
        if iprint >= 2 and majiter > majiter_prev:
            print("%5i %5i % 16.6E % 16.6E" % (majiter,feval[0],
                                               fx,linalg.norm(g)))

        # If exit mode is not -1 or 1, slsqp has completed
        if abs(mode) != 1:
            break

        majiter_prev = int(majiter)

    # Optimization loop complete.  Print status if requested
    if iprint >= 1:
        print(exit_modes[int(mode)] + "    (Exit mode " + str(mode) + ')')
        print("            Current function value:", fx)
        print("            Iterations:", majiter)
        print("            Function evaluations:", feval[0])
        print("            Gradient evaluations:", geval[0])

    return OptimizeResult(x=x, fun=fx, jac=g[:-1], nit=int(majiter),
                          nfev=feval[0], njev=geval[0], status=int(mode),
                          message=exit_modes[int(mode)], success=(mode == 0))

Example 24

Project: zipline Source File: test_statistical.py
    @parameter_space(returns_length=[2, 3], regression_length=[3, 4])
    def test_regression_of_returns_factor(self,
                                          returns_length,
                                          regression_length):
        """
        Tests for the built-in factor `RollingLinearRegressionOfReturns`.
        """
        assets = self.assets
        my_asset = self.my_asset
        my_asset_column = self.my_asset_column
        dates = self.dates
        start_date = self.pipeline_start_date
        end_date = self.pipeline_end_date
        start_date_index = self.start_date_index
        end_date_index = self.end_date_index
        num_days = self.num_days
        run_pipeline = self.run_pipeline

        # The order of these is meant to align with the output of `linregress`.
        outputs = ['beta', 'alpha', 'r_value', 'p_value', 'stderr']

        returns = Returns(window_length=returns_length)
        masks = self.cascading_mask, self.alternating_mask, NotSpecified
        expected_mask_results = (
            self.expected_cascading_mask_result,
            self.expected_alternating_mask_result,
            self.expected_no_mask_result,
        )

        for mask, expected_mask in zip(masks, expected_mask_results):
            regression_factor = RollingLinearRegressionOfReturns(
                target=my_asset,
                returns_length=returns_length,
                regression_length=regression_length,
                mask=mask,
            )

            columns = {
                output: getattr(regression_factor, output)
                for output in outputs
            }
            pipeline = Pipeline(columns=columns)
            if mask is not NotSpecified:
                pipeline.add(mask, 'mask')

            results = run_pipeline(pipeline, start_date, end_date)
            if mask is not NotSpecified:
                mask_results = results['mask'].unstack()
                check_arrays(mask_results.values, expected_mask)

            output_results = {}
            expected_output_results = {}
            for output in outputs:
                output_results[output] = results[output].unstack()
                expected_output_results[output] = full_like(
                    output_results[output], nan,
                )

            # Run a separate pipeline that calculates returns starting
            # (regression_length - 1) days prior to our start date. This is
            # because we need (regression_length - 1) extra days of returns to
            # compute our expected regressions.
            results = run_pipeline(
                Pipeline(columns={'returns': returns}),
                dates[start_date_index - (regression_length - 1)],
                dates[end_date_index],
            )
            returns_results = results['returns'].unstack()

            # On each day, calculate the expected regression results for Y ~ X
            # where Y is the asset we are interested in and X is each other
            # asset. Each regression is calculated over `regression_length`
            # days of data.
            for day in range(num_days):
                todays_returns = returns_results.iloc[
                    day:day + regression_length
                ]
                my_asset_returns = todays_returns.iloc[:, my_asset_column]
                for asset, other_asset_returns in todays_returns.iteritems():
                    asset_column = int(asset) - 1
                    expected_regression_results = linregress(
                        y=other_asset_returns, x=my_asset_returns,
                    )
                    for i, output in enumerate(outputs):
                        expected_output_results[output][day, asset_column] = \
                            expected_regression_results[i]

            for output in outputs:
                output_result = output_results[output]
                expected_output_result = DataFrame(
                    where(expected_mask, expected_output_results[output], nan),
                    index=dates[start_date_index:end_date_index + 1],
                    columns=assets,
                )
                assert_frame_equal(output_result, expected_output_result)

Example 25

Project: cartopy Source File: patch.py
def path_to_geos(path, force_ccw=False):
    """
    Creates a list of Shapely geometric objects from a
    :class:`matplotlib.path.Path`.

    Args:

    * path
        A :class:`matplotlib.path.Path` instance.

    Kwargs:

    * force_ccw
        Boolean flag determining whether the path can be inverted to enforce
        ccw.

    Returns:
        A list of :class:`shapely.geometry.polygon.Polygon`,
        :class:`shapely.geometry.linestring.LineString` and/or
        :class:`shapely.geometry.multilinestring.MultiLineString` instances.

    """
    # Convert path into numpy array of vertices (and associated codes)
    path_verts, path_codes = path_segments(path, curves=False)

    # Split into subarrays such that each subarray consists of connected
    # line segments based on the start of each one being marked by a
    # matplotlib MOVETO code.
    verts_split_inds = np.where(path_codes == Path.MOVETO)[0]
    verts_split = np.split(path_verts, verts_split_inds)
    codes_split = np.split(path_codes, verts_split_inds)

    # Iterate through the vertices generating a list of
    # (external_geom, [internal_polygons]) tuples.
    other_result_geoms = []
    collection = []
    for path_verts, path_codes in zip(verts_split, codes_split):
        if len(path_verts) == 0:
            continue

        # XXX A path can be given which does not end with close poly, in that
        # situation, we have to guess?
        verts_same_as_first = np.all(path_verts[0, :] == path_verts[1:, :],
                                     axis=1)
        if all(verts_same_as_first):
            geom = sgeom.Point(path_verts[0, :])
        elif (path_verts.shape[0] > 2 and
                (path_codes[-1] == Path.CLOSEPOLY or
                 verts_same_as_first[-1])):
            if path_codes[-1] == Path.CLOSEPOLY:
                geom = sgeom.Polygon(path_verts[:-1, :])
            else:
                geom = sgeom.Polygon(path_verts)
        else:
            geom = sgeom.LineString(path_verts)

        # If geom is a Polygon and is contained within the last geom in
        # collection, add it to its list of internal polygons, otherwise
        # simply append it as a new external geom.
        if geom.is_empty:
            pass
        elif (len(collection) > 0 and
                isinstance(collection[-1][0], sgeom.Polygon) and
                isinstance(geom, sgeom.Polygon) and
                collection[-1][0].contains(geom.exterior)):
            collection[-1][1].append(geom.exterior)
        elif isinstance(geom, sgeom.Point):
            other_result_geoms.append(geom)
        else:
            collection.append((geom, []))

    # Convert each (external_geom, [internal_polygons]) pair into a
    # a shapely Polygon that encapsulates the internal polygons, if the
    # external geom is a LineString leave it alone.
    geom_collection = []
    for external_geom, internal_polys in collection:
        if internal_polys:
            # XXX worry about islands within lakes
            geom = sgeom.Polygon(external_geom.exterior, internal_polys)
        else:
            geom = external_geom

        # Correctly orientate the polygon (ccw)
        if isinstance(geom, sgeom.Polygon):
            if force_ccw and not geom.exterior.is_ccw:
                geom = sgeom.polygon.orient(geom)

        geom_collection.append(geom)

    # If the geom_collection only contains LineStrings combine them
    # into a single MultiLinestring.
    if geom_collection and all(isinstance(geom, sgeom.LineString) for
                               geom in geom_collection):
        geom_collection = [sgeom.MultiLineString(geom_collection)]

    # Remove any zero area Polygons
    def not_zero_poly(geom):
        return ((isinstance(geom, sgeom.Polygon) and not geom._is_empty and
                 geom.area != 0) or
                not isinstance(geom, sgeom.Polygon))

    result = list(filter(not_zero_poly, geom_collection))

    return result + other_result_geoms

Example 26

Project: iris Source File: geometry.py
def _extract_relevant_cube_slice(cube, geometry):
    """
    Given a shapely geometry object, this helper method returns
    the tuple
    (subcube, x_coord_of_subcube, y_coord_of_subcube,
     (min_x_index, min_y_index, max_x_index, max_y_index))

    If cube and geometry don't overlap, returns None.

    """

    # Validate the input parameters
    if not cube.coords(axis='x') or not cube.coords(axis='y'):
        raise ValueError('The cube must contain x and y axes.')

    x_coords = cube.coords(axis='x')
    y_coords = cube.coords(axis='y')
    if len(x_coords) != 1 or len(y_coords) != 1:
        raise ValueError('The cube must contain one, and only one, coordinate '
                         'for each of the x and y axes.')

    x_coord = x_coords[0]
    y_coord = y_coords[0]
    if not (x_coord.has_bounds() and y_coord.has_bounds()):
        raise ValueError('Both horizontal coordinates must have bounds.')

    if x_coord.ndim != 1:
        raise iris.exceptions.CoordinateMultiDimError(x_coord)
    if y_coord.ndim != 1:
        raise iris.exceptions.CoordinateMultiDimError(y_coord)

    # bounds of cube dimensions
    x_bounds = x_coord.bounds
    y_bounds = y_coord.bounds

    # identify ascending/descending coordinate dimensions
    x_ascending = x_coord.points[1] - x_coord.points[0] > 0.
    y_ascending = y_coord.points[1] - y_coord.points[0] > 0.

    # identify upper/lower bounds of coordinate dimensions
    x_bounds_lower = x_bounds[:, 0] if x_ascending else x_bounds[:, 1]
    y_bounds_lower = y_bounds[:, 0] if y_ascending else y_bounds[:, 1]
    x_bounds_upper = x_bounds[:, 1] if x_ascending else x_bounds[:, 0]
    y_bounds_upper = y_bounds[:, 1] if y_ascending else y_bounds[:, 0]

    # find indices of coordinate bounds to fully cover geometry
    x_min_geom, y_min_geom, x_max_geom, y_max_geom = geometry.bounds
    try:
        x_min_ix = np.where(x_bounds_lower <= x_min_geom)[0]
        x_min_ix = x_min_ix[np.argmax(x_bounds_lower[x_min_ix])]
    except ValueError:
        warnings.warn("The geometry exceeds the cube's x dimension at the "
                      "lower end.", UserWarning)
        x_min_ix = 0 if x_ascending else x_coord.points.size - 1

    try:
        x_max_ix = np.where(x_bounds_upper >= x_max_geom)[0]
        x_max_ix = x_max_ix[np.argmin(x_bounds_upper[x_max_ix])]
    except ValueError:
        warnings.warn("The geometry exceeds the cube's x dimension at the "
                      "upper end.", UserWarning)
        x_max_ix = x_coord.points.size - 1 if x_ascending else 0

    try:
        y_min_ix = np.where(y_bounds_lower <= y_min_geom)[0]
        y_min_ix = y_min_ix[np.argmax(y_bounds_lower[y_min_ix])]
    except ValueError:
        warnings.warn("The geometry exceeds the cube's y dimension at the "
                      "lower end.", UserWarning)
        y_min_ix = 0 if y_ascending else y_coord.points.size - 1

    try:
        y_max_ix = np.where(y_bounds_upper >= y_max_geom)[0]
        y_max_ix = y_max_ix[np.argmin(y_bounds_upper[y_max_ix])]
    except ValueError:
        warnings.warn("The geometry exceeds the cube's y dimension at the "
                      "upper end.", UserWarning)
        y_max_ix = y_coord.points.size - 1 if y_ascending else 0

    # extract coordinate values at these indices
    x_min = x_bounds_lower[x_min_ix]
    x_max = x_bounds_upper[x_max_ix]
    y_min = y_bounds_lower[y_min_ix]
    y_max = y_bounds_upper[y_max_ix]

    # switch min and max if necessary, to create slice objects later on
    if x_min_ix > x_max_ix:
        x_min_ix, x_max_ix = x_max_ix, x_min_ix
    if y_min_ix > y_max_ix:
        y_min_ix, y_max_ix = y_max_ix, y_min_ix
    bnds_ix = x_min_ix, y_min_ix, x_max_ix, y_max_ix

    # cut the relevant part from the original cube
    coord_constr = {x_coord.name(): lambda x: x_min <= x.point <= x_max,
                    y_coord.name(): lambda y: y_min <= y.point <= y_max}
    constraint = iris.Constraint(coord_values=coord_constr)
    subcube = cube.extract(constraint)

    if subcube is None:
        return None
    else:
        x_coord = subcube.coord(axis='x')
        y_coord = subcube.coord(axis='y')
        return subcube, x_coord, y_coord, bnds_ix

Example 27

Project: python-skyfield Source File: relativity.py
def add_deflection(position, observer, ephemeris, t,
                   include_earth_deflection, count=3):
    """Update `position` for how solar system masses will deflect its light.

    Given the ICRS `position` [x,y,z] of an object (au) that is being
    viewed from the `observer` also expressed as [x,y,z], and given an
    ephemeris that can be used to determine solar system body positions,
    and given the time `t` and Boolean `apply_earth` indicating whether
    to worry about the effect of Earth's mass, and a `count` of how many
    major solar system bodies to worry about, this function updates
    `position` in-place to show how the masses in the solar system will
    deflect its image.

    """
    # Compute light-time to observed object.

    tlt = length_of(position) / C_AUDAY

    # Cycle through gravitating bodies.

    jd_tdb = t.tdb
    ts = t.ts
    for name in deflectors[:count]:
        try:
            deflector = ephemeris[name]
        except KeyError:
            deflector = ephemeris[name + ' barycenter']

        # Get position of gravitating body wrt ss barycenter at time 't_tdb'.

        bposition = deflector.at(ts.tdb(jd=jd_tdb)).position.au  # TODO

        # Get position of gravitating body wrt observer at time 'jd_tdb'.

        gpv = bposition - observer

        # Compute light-time from point on incoming light ray that is closest
        # to gravitating body.

        dlt = light_time_difference(position, gpv)

        # Get position of gravitating body wrt ss barycenter at time when
        # incoming photons were closest to it.

        tclose = jd_tdb

        # if dlt > 0.0:
        #     tclose = jd - dlt

        tclose = where(dlt > 0.0, jd_tdb - dlt, tclose)
        tclose = where(tlt < dlt, jd_tdb - tlt, tclose)

        # if tlt < dlt:
        #     tclose = jd - tlt

        bposition = deflector.at(ts.tdb(jd=tclose)).position.au  # TODO
        rmass = rmasses[name]
        _add_deflection(position, observer, bposition, rmass)

    # If observer is not at geocenter, add in deflection due to Earth.

    if include_earth_deflection.any():
        deflector = ephemeris['earth']
        bposition = deflector.at(ts.tdb(jd=tclose)).position.au  # TODO
        rmass = rmasses['earth']
        # TODO: Make the following code less messy, maybe by having
        # _add_deflection() return a new vector instead of modifying the
        # old one in-place.
        deflected_position = position.copy()
        _add_deflection(deflected_position, observer, bposition, rmass)
        if include_earth_deflection.shape:
            position[:,include_earth_deflection] = (
                deflected_position[:,include_earth_deflection])
        else:
            position[:] = deflected_position[:]

Example 28

Project: gensim Source File: ldamodel.py
    def top_topics(self, corpus, num_words=20):
        """
        Calculate the Umass topic coherence for each topic. Algorithm from
        **Mimno, Wallach, Talley, Leenders, McCallum: Optimizing Semantic Coherence in Topic Models, CEMNLP 2011.**
        """
        is_corpus, corpus = utils.is_corpus(corpus)
        if not is_corpus:
            logger.warning("LdaModel.top_topics() called with an empty corpus")
            return

        topics = []
        str_topics = []
        for topic in self.state.get_lambda():
            topic = topic / topic.sum()  # normalize to probability distribution
            bestn = matutils.argsort(topic, topn=num_words, reverse=True)
            topics.append(bestn)
            beststr = [(topic[id], self.id2word[id]) for id in bestn]
            str_topics.append(beststr)

        # top_ids are limited to every topics top words. should not exceed the
        # vocabulary size.
        top_ids = set(chain.from_iterable(topics))

        # create a docuement occurence sparse matrix for each word
        doc_word_list = {}
        for id in top_ids:
            id_list = set()
            for n, docuement in enumerate(corpus):
                if id in frozenset(x[0] for x in docuement):
                    id_list.add(n)

            doc_word_list[id] = id_list

        coherence_scores = []
        for t, top_words in enumerate(topics):
            # Calculate each coherence score C(t, top_words)
            coherence = 0.0
            # Sum of top words m=2..M
            for m in top_words[1:]:
                # m_docs is v_m^(t)
                m_docs = doc_word_list[m]
                m_index = numpy.where(top_words == m)[0]

                # Sum of top words l=1..m-1
                # i.e., all words ranked higher than the current word m
                for l in top_words[:m_index - 1]:
                    # l_docs is v_l^(t)
                    l_docs = doc_word_list[l]

                    # make sure this word appears in some docuements.
                    if len(l_docs) > 0:
                        # co_doc_frequency is D(v_m^(t), v_l^(t))
                        co_doc_frequency = len(m_docs.intersection(l_docs))

                        # add to the coherence sum for these two words m, l
                        coherence += numpy.log((co_doc_frequency + 1.0) / len(l_docs))

            coherence_scores.append((str_topics[t], coherence))

        top_topics = sorted(coherence_scores, key=lambda t: t[1], reverse=True)
        return top_topics

Example 29

Project: pyhawkes Source File: parents.py
Function: init
    def __init__(self, model, T, S, F):
        """
        Initialize a parent array Z of size TxKxKxB to model the
        event parents for data matrix S (TxK) which has been filtered
        to create filtered data array F (TxKxB).

        Also create a background parent array of size TxK to specify
        how many events are attributed to the background process.

        :param T: Number of time bins
        :param K: Number of processes
        :param B: Number of basis functions
        :param S: Data matrix (TxK)
        :param F: Filtered data matrix (TxKxB)
        """
        self.model = model
        self.dt = model.dt
        self.K = model.K
        self.B = model.B

        # TODO: Remove dependencies on S and F
        self.T = T
        self.S = S
        self.F = F

        # Save sparse versions of S and F
        self.ts = []
        self.Ts = []
        self.Ns = []
        self.Ss = []
        self.Fs = []
        for k in xrange(self.K):
            # Find the rows where S[:,k] is nonzero
                tk = np.where(S[:,k])[0]
                self.ts.append(tk)
                self.Ts.append(len(tk))
                self.Ss.append(S[tk,k].astype(np.uint32))
                self.Ns.append(S[tk,k].sum())
                self.Fs.append(F[tk])
        self.Ns = np.array(self.Ns)


        # The base class handles the parent variables
        # We use a sparse representation that only considers times (rows)
        # where there is a spike
        self._Z = None
        self._EZ = None

        # Initialize GSL RNGs for resampling Z
        self.pyrngs = initialize_pyrngs()

Example 30

Project: mpop Source File: tools.py
Function: estimate_cth
def estimate_cth(IR_108, cth_atm="standard"):

    '''
    Estimation of the cloud top height using the 10.8 micron channel
    limitations: this is the most simple approach
                 a simple fit of the ir108 to the temperature profile
                 * no correction for water vapour or any other trace gas
                 * no viewing angle dependency
                 * no correction for semi-transparent clouds

    optional input:
      cth_atm    * "standard", "tropics", "midlatitude summer", "midlatitude winter", "subarctic summer", "subarctic winter"
                  Matching the 10.8 micron temperature with atmosphere profile
                  (s)  AFGL atmospheric constituent profile. U.S. standard atmosphere 1976. (AFGL-TR-86-0110) 
                  (t)  AFGL atmospheric constituent profile. tropical.                      (AFGL-TR-86-0110)
                  (mw) AFGL atmospheric constituent profile. midlatitude summer.            (AFGL-TR-86-0110) 
                  (ms) AFGL atmospheric constituent profile. midlatitude winter.            (AFGL-TR-86-0110)
                  (ss) AFGL atmospheric constituent profile. subarctic summer.              (AFGL-TR-86-0110) 
                  (sw) AFGL atmospheric constituent profile. subarctic winter.              (AFGL-TR-86-0110)
                  Ulrich Hamann (MeteoSwiss)
                * "tropopause"
                  Assuming a fixed tropopause height and a fixed temperature gradient
                  Richard Mueller (DWD)
    output: 
      parallax corrected channel
                 the content of the channel will be parallax corrected.
                 The name of the new channel will be
                 *original_chan.name+'_PC'*, eg. "IR_108_PC". This name is
                 also stored to the info dictionary of the originating channel.

    Versions: 05.07.2016 initial version
              Ulrich Hamann (MeteoSwiss), Richard Mueller (DWD)
    '''

    print "*** estimating CTH using the 10.8 micro meter brightness temperature "

    if cth_atm.lower() != "tropopause":

        # define atmospheric temperature profile    
        import os
        from numpy import loadtxt, zeros, where, logical_and
        import mpop 

        mpop_dir = os.path.dirname(mpop.__file__)
        afgl_file = mpop_dir+"/afgl.dat"
        print "... assume ", cth_atm, " atmosphere for temperature profile"

        if cth_atm.lower()=="standard" or cth_atm.lower()=="s":
            z, T = loadtxt(afgl_file, usecols=(0, 1), unpack=True, comments="#")
        elif cth_atm.lower()=="tropics" or cth_atm.lower()=="t":
            z, T = loadtxt(afgl_file, usecols=(0, 2), unpack=True, comments="#")
        elif cth_atm.lower()=="midlatitude summer" or cth_atm.lower()=="ms":
            z, T = loadtxt(afgl_file, usecols=(0, 3), unpack=True, comments="#")
        elif cth_atm.lower()=="midlatitude winter" or cth_atm.lower()=="ws":
            z, T = loadtxt(afgl_file, usecols=(0, 4), unpack=True, comments="#")
        elif cth_atm.lower()=="subarctic summer" or cth_atm.lower()=="ss":
            z, T = loadtxt(afgl_file, usecols=(0, 5), unpack=True, comments="#")
        elif cth_atm.lower()=="subarctic winter" or cth_atm.lower()=="ss":
            z, T = loadtxt(afgl_file, usecols=(0, 6), unpack=True, comments="#")
        else:
            print "*** Error in estimate_cth (mpop/tools.py)"
            print "unknown temperature profiel for CTH estimation: cth_atm = ", cth_atm
            quit()

        height = zeros(IR_108.shape)
        # warmer than lowest level -> clear sky 
        height[where(IR_108 > T[-1])] = -1.
        print "     z0(km)   z1(km)   T0(K)   T1(K)  number of pixels"
        print "------------------------------------------------------"
        for i in range(z.size)[::-1]:

            # search for temperatures between layer i-1 and i
            ind =  np.where( logical_and( T[i-1]< IR_108, IR_108 < T[i]) )
            # interpolate CTH according to ir108 temperature
            height[ind] = z[i] + (IR_108[ind]-T[i])/(T[i-1]-T[i]) * (z[i-1]-z[i])
            # verbose output
            print " {0:8.1f} {1:8.1f} {2:8.1f} {3:8.1f} {4:8d}".format(z[i], z[i-1], T[i], T[i-1], len(ind[0]))

            # if temperature increases above 8km -> tropopause detected
            if z[i]>=8. and T[i] <= T[i-1]:
                # no cloud above tropopose
                break
            # no cloud heights above 20km
            if z[i]>=20.:
                break

        # if height is still 0 -> cloud colder than tropopause -> cth == tropopause height
        height[np.where( height == 0 )] = z[i]
        
    else:

        Htropo=11.0 # km
        # this is an assumption it should be optimized 
        # by making it dependent on region and season. 
        # It might be good to include the ITC in the  
        # region of interest, that would make a fixed Htropo 
        # value more reliable. 
        Tmin = np.amin(IR_108) 
        # for Tmin it might be better to use the 5th or 10th percentile 
        # else overshoting tops induces further uncertainties  
        # in the calculation of the cloud height. 
        # However numpy provides weird results for 5th percentile. 
        # Hence, for the working version the minima is used 

        print "... assume tropopause height ", Htropo, ", tropopause temperature ", Tmin, "K (", Tmin-273.16, "deg C)"
        print "    and constant temperature gradient 6.5 K/km"

        height = -(IR_108 - Tmin)/6.5 + Htropo 
        # calculation of the height, the temperature gradient 
        # 6.5 K/km is an assumption  
        # derived from USS and MPI standard profiles. It 
        # has to be improved as well 

    # convert to masked array
    # convert form km to meter
    height = np.ma.masked_where(height <= 0, height, copy=False) * 1000.

    if False:
        from trollimage.image import Image as trollimage
        from trollimage.colormap import rainbow
        from copy import deepcopy 
        # cloud top height
        prop = height
        min_data = prop.min()
        max_data = prop.max()
        print " estimated CTH(meter) (min/max): ", min_data, max_data
        min_data =     0
        max_data = 12000    
        colormap = deepcopy(rainbow)
        colormap.set_range(min_data, max_data)
        img = trollimage(prop, mode="L") #, fill_value=[0,0,0]
        img.colorize(colormap)
        img.show()

    # return cloud top height in meter
    return height

Example 31

Project: PyFITS Source File: image.py
    def _get_scaled_image_data(self, offset, shape):
        """
        Internal function for reading image data from a file and apply scale
        factors to it.  Normally this is used for the entire image, but it
        supports alternate offset/shape for Section support.
        """

        code = BITPIX2DTYPE[self._orig_bitpix]

        raw_data = self._get_raw_data(shape, code, offset)
        raw_data.dtype = raw_data.dtype.newbyteorder('>')

        if self._do_not_scale_image_data or (
                self._orig_bzero == 0 and self._orig_bscale == 1 and
                self._blank is None):
            # No further conversion of the data is necessary
            return raw_data

        try:
            if self._file.strict_memmap:
                raise ValueError("Cannot load a memory-mapped image: "
                                 "BZERO/BSCALE/BLANK header keywords present. "
                                 "Set memmap=False.")
        except AttributeError:  # strict_memmap not set
            pass

        data = None
        if not (self._orig_bzero == 0 and self._orig_bscale == 1):
            data = self._convert_pseudo_unsigned(raw_data)

        if data is None:
            # In these cases, we end up with floating-point arrays and have to
            # apply bscale and bzero. We may have to handle BLANK and convert
            # to NaN in the resulting floating-point arrays.
            # The BLANK keyword should only be applied for integer data (this
            # is checked in __init__ but it can't hurt to double check here)
            blanks = None

            if self._blank is not None and self._bitpix > 0:
                blanks = raw_data.flat == self._blank
                # The size of blanks in bytes is the number of elements in
                # raw_data.flat.  However, if we use np.where instead we will
                # only use 8 bytes for each index where the condition is true.
                # So if the number of blank items is fewer than
                # len(raw_data.flat) / 8, using np.where will use less memory
                if blanks.sum() < len(blanks) / 8:
                    blanks = np.where(blanks)

            new_dtype = self._dtype_for_bitpix()
            if new_dtype is not None:
                data = np.array(raw_data, dtype=new_dtype)
            else:  # floating point cases
                if self._file is not None and self._file.memmap:
                    data = raw_data.copy()
                elif not raw_data.flags.writeable:
                    # create a writeable copy if needed
                    data = raw_data.copy()
                # if not memmap, use the space already in memory
                else:
                    data = raw_data

            del raw_data

            if self._orig_bscale != 1:
                np.multiply(data, self._orig_bscale, data)
            if self._orig_bzero != 0:
                data += self._orig_bzero

            if self._blank:
                data.flat[blanks] = np.nan

        return data

Example 32

Project: fast-rcnn Source File: minibatch.py
def _sample_rois(roidb, fg_rois_per_image, rois_per_image, num_classes):
    """Generate a random sample of RoIs comprising foreground and background
    examples.
    """
    # label = class RoI has max overlap with
    labels = roidb['max_classes']
    overlaps = roidb['max_overlaps']
    rois = roidb['boxes']

    # Select foreground RoIs as those with >= FG_THRESH overlap
    fg_inds = np.where(overlaps >= cfg.TRAIN.FG_THRESH)[0]
    # Guard against the case when an image has fewer than fg_rois_per_image
    # foreground RoIs
    fg_rois_per_this_image = np.minimum(fg_rois_per_image, fg_inds.size)
    # Sample foreground regions without replacement
    if fg_inds.size > 0:
        fg_inds = npr.choice(fg_inds, size=fg_rois_per_this_image,
                             replace=False)

    # Select background RoIs as those within [BG_THRESH_LO, BG_THRESH_HI)
    bg_inds = np.where((overlaps < cfg.TRAIN.BG_THRESH_HI) &
                       (overlaps >= cfg.TRAIN.BG_THRESH_LO))[0]
    # Compute number of background RoIs to take from this image (guarding
    # against there being fewer than desired)
    bg_rois_per_this_image = rois_per_image - fg_rois_per_this_image
    bg_rois_per_this_image = np.minimum(bg_rois_per_this_image,
                                        bg_inds.size)
    # Sample foreground regions without replacement
    if bg_inds.size > 0:
        bg_inds = npr.choice(bg_inds, size=bg_rois_per_this_image,
                             replace=False)

    # The indices that we're selecting (both fg and bg)
    keep_inds = np.append(fg_inds, bg_inds)
    # Select sampled values from various arrays:
    labels = labels[keep_inds]
    # Clamp labels for the background RoIs to 0
    labels[fg_rois_per_this_image:] = 0
    overlaps = overlaps[keep_inds]
    rois = rois[keep_inds]

    bbox_targets, bbox_loss_weights = \
            _get_bbox_regression_labels(roidb['bbox_targets'][keep_inds, :],
                                        num_classes)

    return labels, overlaps, rois, bbox_targets, bbox_loss_weights

Example 33

Project: statsmodels Source File: data.py
    def _handle_constant(self, hasconst):
        if hasconst is not None:
            if hasconst:
                self.k_constant = 1
                self.const_idx = None
            else:
                self.k_constant = 0
                self.const_idx = None
        elif self.exog is None:
            self.const_idx = None
            self.k_constant = 0
        else:
            # detect where the constant is
            check_implicit = False
            const_idx = np.where(self.exog.ptp(axis=0) == 0)[0].squeeze()
            self.k_constant = const_idx.size

            if self.k_constant == 1:
                if self.exog[:, const_idx].mean() != 0:
                    self.const_idx = const_idx
                else:
                    # we only have a zero column and no other constant
                    check_implicit = True
            elif self.k_constant > 1:
                # we have more than one constant column
                # look for ones
                values = []  # keep values if we need != 0
                for idx in const_idx:
                    value = self.exog[:, idx].mean()
                    if value == 1:
                        self.k_constant = 1
                        self.const_idx = idx
                        break
                    values.append(value)
                else:
                    # we didn't break, no column of ones
                    pos = (np.array(values) != 0)
                    if pos.any():
                        # take the first nonzero column
                        self.k_constant = 1
                        self.const_idx = const_idx[pos.argmax()]
                    else:
                        # only zero columns
                        check_implicit = True
            elif self.k_constant == 0:
                check_implicit = True
            else:
                # shouldn't be here
                pass

            if check_implicit:
                # look for implicit constant
                # Compute rank of augmented matrix
                augmented_exog = np.column_stack(
                            (np.ones(self.exog.shape[0]), self.exog))
                rank_augm = np_matrix_rank(augmented_exog)
                rank_orig = np_matrix_rank(self.exog)
                self.k_constant = int(rank_orig == rank_augm)
                self.const_idx = None

Example 34

Project: pyspace Source File: frequency_features.py
    def _execute(self, x):
        """ Extract the Fourier features from the given data x """
        # Lazily set NFFT and noverlap
        if not hasattr(self, "NFFT"):
            self.NFFT = (x.shape[0] // 4) * 2
            self.noverlap = self.NFFT * 9 // 10
            
        # For each channel, we compute the band power in the specified
        # frequency bands
        features = []
        feature_names = []
        for channel_name in x.channel_names:
            channel_index = x.channel_names.index(channel_name)
            # Compute spectrogram
            (Pxx, freqs, bins) = mlab.specgram(x[:, channel_index],
                                               Fs = x.sampling_frequency,
                                               NFFT = self.NFFT,
                                               noverlap = self.noverlap)
            
            for min_frequency, max_frequency in self.frequency_bands:
                # Just to make sure....
                min_frequency = float(min_frequency)
                max_frequency = float(max_frequency)
    
                # Compute band powers in the specified frequency range for all
                # bins
                selected_freqs = Pxx[numpy.where((freqs >= min_frequency)
                                                & (freqs <= max_frequency))]
                band_powers = numpy.sum(selected_freqs, 0)
                                                          
                # Compute the corresponding feature based on the criterion
                # that was specified
                if self.criterion == "max":
                    channel_features = [10 * numpy.log10(max(band_powers))]
                elif self.criterion == "min":
                    channel_features = [10 * numpy.log10(min(band_powers))]
                elif self.criterion == "median":
                    channel_features = \
                           [10 * numpy.log10(numpy.median(band_powers))]
                elif self.criterion == "max_change":
                    channel_features = [10*numpy.log10((max(band_powers)- \
                                                           min(band_powers)))]
                elif self.criterion == "all":
                    channel_features = 10 * numpy.log10(band_powers)
    
                # Append the feature
                features.extend(channel_features)
                # Determine feature names
                if self.criterion != "all":
                    feature_names.append("BP_%s_%.2fHz_%.2fHz_%s" % 
                                         (channel_name, min_frequency,
                                          max_frequency, self.criterion))
                else:
                    feature_names.extend("BP_%s_%.2fHz_%.2fHz_%s_%.2fsec" % 
                                         (channel_name, min_frequency,
                                          max_frequency, self.criterion,
                                          bin) for bin in bins)
        
        feature_vector = \
         FeatureVector(numpy.atleast_2d(features).astype(numpy.float64),
         feature_names)
        return feature_vector

Example 35

Project: pyNastran Source File: structured_grid_maker.py
def Grid(object):
    def __init__(self, bdf_name, debug=True):
        self.bdf = BDF(debug=debug)
        self.bdf.read_bdf(bdf_name, xref=True)

        self.xyz = {}

    def build_2d_mesh(self):
        eids = arange(4719, 9878)
        starting_eid = 9876
        starting_element = self.bdf.elements[starting_eid]

        (edges_dict, edge_to_eids_dict) = self.build_edge_dicts()

        # x
        edges = edges_dict[starting_eid]

        # di = [eid1, eid2, eid3, eid4]
        d = [edge_to_eids_dict[edge] for edge in edges]

        # Li = [1, 1, 2, 2]
        L = array([len(di) for di in d])
        print("L = %s" % L)
        i = where(L==2)[0]
        assert len(i) == 2, 'eid=%s is not a corner' % starting_eid
        x = i[0]
        y = i[1]

        edge_x = edges[x]
        dx = d[x]
        next_eid_x = dx.pop(starting_eid)


        edge_y = edges[y]
        dy = d[y]
        next_eid_y = dy.pop(starting_eid)

        # x direction
        nx, eids_x = self.walk_edge(starting_eid, edge_x)
        ny, eids_y = self.walk_edge(starting_eid, edge_y)

        starting_nids = starting_element.node_ids
        set_starting_nids = set(starting_nids)
        non_corner_nids = set(edge_x) + set(edge_y)
        corner_nid = set_starting_nids - non_corner_nids
        assert len(corner_nid) == 1, corner_nid

        positions = {}
        for nid, node in self.bdf.nodes:
            positions[nid] = node.get_position()

        # size the array
        xyz = array( (nx, ny, 3), 'float32')
        xyz[0, 0, :] = positions[corner_nid]

        # 7----8
        # | C  |
        # 4----3----6
        # | A  |  B |
        # 1----2----5  ---> x
        #
        # A - starting element
        # 2-3 is edge_x
        # 4-3 is edge_y
        # 1 - corner node and is on the boundary
        # 2 - part of edge_x that's not in edge_y

        # 5-6 is edge_x2
        # A*B = |A||B|cos(theta)
        # 5 - the node that has the smaller angle
        node_2x = list(set(edge_x) - set(edge_y))
        assert len(node_2x) == 1
        node_2x = node_2x[0]
        v12 = positions[node_2x] - positions[corner_nid]
        n12 = norm(v12)

        node_4y = list(set(edge_y) - set(edge_x))  # node 4
        assert len(node_4y) == 1
        node_4y = node_4y[0]
        v14 = positions[node_4y] - positions[corner_nid]
        n14 = norm(v14)

        # find the 0th row
        for icol, eid in enumerate(eids_x):
            # walk_coords(corner_nid)
            element = self.bdf.elements[eid]
            nids = element.node_ids
            nids.pop(corner_nid)

            (n2, n3, n4) = nids
            edge1 = [corner_nid, n2]
            edge2 = [corner_nid, n3]
            edge3 = [corner_nid, n4]
            #edge1.sort()
            #edge2.sort()
            #edge3.sort()
            edges = [tuple(edge1), tuple(edge2), tuple(edge3)]

            A = array([ (positions[nA] - positions[nB]) for (nA, nB) in edges])
            cosTs = dot(v12, A)/(n12*norm(A))
            thetas = acos(cosTs)

            itheta_min = where(theta==theta.min())[0][0]
            edge_nid = edge_nid_choices[itheta_min]
            xyzi = positions[edge_nid]
            xyz[0, icol, :] = xyzi

        for ix in range(nx):
            for iy in range(ny):
                xyz[0, icol, :] = xyzi


    def walk_edge(self, starting_eid, edge_x):
        nx = 1
        eids_x = [starting_eid]
        while 1:
            eid = edge_to_eids_dict[edge_x][0]
            edges = edges_dict[eid]
            iedge0 = edges.index(edge_x)

            # go to the opposite edge of the element
            iedge_new = iedge0 + 2

            # correct the edge number
            iedge = iedge_new if iedge_new < 4 else iedge_new - 4
            new_edge = edges[iedge]
            eids_next = edge_to_eids_dict[edge]
            eids_next.pop(eid)
            if len(eids_next) == 0:
                break
            eid = eids_next[0]

            nx += 1
            eids_x.append(eid)
        return nx, eids_x

    def build_edge_dicts(self, eids):
        edges_dict = {}
        edge_to_eids_dict = {}

        for eid, element in self.bdf.elements:
            if isinstance(element, CQUAD4):
                (n1, n2, n3, n4) = element.node_ids
                edge1 = [n1, n2]
                edge2 = [n2, n3]
                edge3 = [n3, n4]
                edge4 = [n4, n1]
                edge1.sort()
                edge2.sort()
                edge3.sort()
                edge4.sort()
                edges = [tuple(edge1), tuple(edge2), tuple(edge3), tuple(edge4)]
            else:
                raise RuntimeError('only CQUAD4s are supported')

            edges_dict[eid] = edges
            for edge in edges:
                if edge in edge_to_eids_dict:
                    edge_to_eids_dict[edge].append(eid)
                else:
                    edge_to_eids_dict[edge] = [eid]
        return edges_dict, edge_to_eids_dict

Example 36

Project: rlpy Source File: KernelizediFDD.py
Function: post_discover
    def post_discover(self, s, terminal, a, td_error, phi_s=None):
        if phi_s is None:
            phi_s = self.phi(s, terminal)
        phi_s_unnorm = self.phi_raw(s, terminal)
        discovered = 0
        Q = self.Qs(s, terminal, phi_s=phi_s).reshape(-1, 1)
        # indices of active features
        active_indices = list(
            np.where(phi_s_unnorm > self.active_threshold)[0])
        # "active indices", active_indices
        # gather all dimensions regarded by active features
        active_dimensions = np.zeros((len(s)), dtype="int")
        closest_neighbor = np.zeros((len(s)))
        for i in active_indices:
            for j in self.features[i].dim:
                active_dimensions[j] += 1
                closest_neighbor[j] = max(closest_neighbor[j], phi_s_unnorm[i])

        # add new base features for all dimension not regarded
        for j in xrange(len(s)):
            if active_dimensions[j] < self.max_active_base_feat and (closest_neighbor[j] < self.max_base_feat_sim or active_dimensions[j] < 1):
                active_indices.append(self.add_base_feature(s, j, Q=Q))
                discovered += 1

        # update relevance statistics of all feature candidates
        if discovered:
            phi_s = self.phi(s, terminal)
        la = len(active_indices)
        if la * (la - 1) < len(self.candidates):
            for ind, cand in self.candidates.items():
                g, h = ind
                rel = self.update_relevance_stat(
                    cand,
                    g,
                    h,
                    td_error,
                    s,
                    a,
                    phi_s)
                self.max_relevance = max(rel, self.max_relevance)
                # add if relevance is high enough
                if rel > self.discover_threshold:
                    self.add_refined_feature(g, h, Q=Q)
                    discovered += 1

        else:
            # the result of both branches can be very different as this one
            # updates only combinations which are considered active.
            for g, h in combinations(active_indices, 2):
                # note: g, h are ordered as active_indices are ordered
                cand = self.candidates.get((g, h))
                if cand is None:
                    continue
                rel = self.update_relevance_stat(
                    cand,
                    g,
                    h,
                    td_error,
                    s,
                    a,
                    phi_s)
                self.max_relevance = max(rel, self.max_relevance)
                # add if relevance is high enough
                if rel > self.discover_threshold:
                    self.add_refined_feature(g, h, Q=Q)
                    discovered += 1

        if discovered:
            self.max_relevance = 0.
        return discovered

Example 37

Project: pyNastran Source File: write_mesh.py
    def _write_elements_properties(self, bdf_file, size):
        """
        Writes the elements and properties in and interspersed order
        """
        #return self._write_elements_properties2(f, size)
        msg = []
        missing_properties = []
        if self.properties:
            bdf_file.write('$ELEMENTS_WITH_PROPERTIES\n')

        eids_written = []
        #pids = sorted(self.properties.keys())

        ptypes = [
            self.properties_shell.pshell,
            self.properties_shell.pcomp,
            self.pshear,
            self.prod,
            self.properties_solid.psolid,

            #self.properties_bar.pbar,
            #self.properties_bar.pbarl,
            #self.properties_beam.pbeam,
            #self.properties_beam.pbeaml,
        ]

        n = 0
        pids_all = None  # the actual properties
        for t in ptypes:
            if t.n and n == 0:
                pids_all = t.property_id
                n = 1
            elif t.n:
                self.log.debug(pids_all)
                self.log.debug(t.property_id)
                try:
                    pids_all = concatenate(pids_all, t.property_id)
                except ValueError:
                    pids_all = array(list(pids_all) + list(t.property_id))

        etypes = (self.elements_shell._get_types() +
                  self.elements_solid._get_types() +
                  [self.crod, self.cshear])

        #pids_set = None
        if pids_all is None:
            bdf_file.write('$MISSING_ELEMENTS because there are no properties\n')
            for t in etypes:
                #print("t.type =", t.type)
                t.write_card(bdf_file, size=size)
            return

        # there are properties
        pids_set = set(list(pids_all))

        n = 0
        pids = None
        for t in etypes:
            #print("t.type =", t.type)
            if t.n and n == 0:
                eids = t.element_id
                pids = t.property_id
                n = 1
            elif t.n:
                try:
                    eids = concatenate(eids, t.element_id)
                #except AttributeError:
                    #eids = array(list(eids) + list(t.element_id))
                except TypeError:
                    #print eids
                    #print t.element_id
                    eids = array(list(eids) + list(t.element_id))
                except ValueError:
                    #print eids
                    #print t.element_id
                    eids = array(list(eids) + list(t.element_id))

                try:
                    pids = concatenate(pids, t.property_id)
                except AttributeError:
                    pids = array(list(pids) + list(t.property_id))
                except TypeError:
                    pids = array(list(pids) + list(t.property_id))
                except ValueError:
                    pids = array(list(pids) + list(t.property_id))
            #else:
                #print t.type

        elements_by_pid = {}
        if pids is not None:
            pids_unique = unique(pids)
            self.log.debug("pids_unique = %s" % pids_unique)
            pids_unique.sort()
            if len(pids_unique) > 0:
                bdf_file.write('$ELEMENTS_WITH_PROPERTIES\n')

            for pid in pids_all:
                i = where(pid == pids)[0]
                eids2 = eids[i]

                for t in ptypes:
                    if t.n and pid in t.property_id:
                        self.log.debug("prop.type = %s" % t.type)
                        t.write_card(bdf_file, size=size, property_ids=[pid])
                        pids_set.remove(pid)
                n = 0
                for t in etypes:
                    if not t.n:
                        continue
                    eids3 = intersect1d(t.element_id, eids2, assume_unique=False)
                    #print("eids3[pid=%s]" %(pid), eids3)
                    if n == 0 and len(eids3):
                        elements_by_pid[pid] = eids3
                        n = 1
                    elif len(eids3):
                        try:
                            c = concatenate(elements_by_pid[pid], eids3)
                        except TypeError:
                            c = array(list(elements_by_pid[pid]) + list(eids3))
                        except ValueError:
                            c = array(list(elements_by_pid[pid]) + list(eids3))
                        elements_by_pid[pid] = c
                    else:
                        continue
                    try:
                        t.write_card(bdf_file, size=size, element_ids=eids3)
                    except TypeError:
                        print("t.type = %s" % t.type)
                        raise
                    del eids3
            #for pid, elements in elements_by_pid.items():
                #print("pid=%s n=%s" % (pid, len(elements)))
            #print elements_by_pid

        # missing properties
        if pids_set:
            pids_list = list(pids_set)
            bdf_file.write('$UNASSOCIATED_PROPERTIES\n')
            for pid in pids_list:
                for t in ptypes:
                    if t.n and pid in t.property_id:
                        t.write_card(bdf_file, size=size, property_ids=[pid])

        #.. todo:: finish...
        bdf_file.write('$UNASSOCIATED_ELEMENTS\n')

Example 38

Project: zipline Source File: test_factor.py
    @parameter_space(seed=[1, 2, 3])
    def test_quantiles_masked(self, seed):
        permute = partial(permute_rows, seed)

        # 7 x 7 so that we divide evenly into 2/3/6-tiles after including the
        # nan value in each row.
        shape = (7, 7)

        # Shuffle the input rows to verify that we don't depend on the order.
        # Take the log to ensure that we don't depend on linear scaling or
        # integrality of inputs
        factor_data = permute(log1p(arange(49, dtype=float).reshape(shape)))
        factor_data_w_nans = where(
            permute(rot90(self.eye_mask(shape=shape))),
            factor_data,
            nan,
        )
        mask_data = permute(self.eye_mask(shape=shape))

        f = F()
        f_nans = OtherF()
        m = Mask()

        # Apply the same shuffle we applied to the input rows to our
        # expectations. Doing it this way makes it obvious that our
        # expectation corresponds to our input, while still testing against
        # a range of input orderings.
        permuted_array = compose(permute, partial(array, dtype=int64_dtype))

        self.check_terms(
            terms={
                '2_masked': f.quantiles(bins=2, mask=m),
                '3_masked': f.quantiles(bins=3, mask=m),
                '6_masked': f.quantiles(bins=6, mask=m),
                '2_nans': f_nans.quantiles(bins=2),
                '3_nans': f_nans.quantiles(bins=3),
                '6_nans': f_nans.quantiles(bins=6),
            },
            initial_workspace={
                f: factor_data,
                f_nans: factor_data_w_nans,
                m: mask_data,
            },
            expected={
                # Expected results here are the same as in
                # test_quantiles_unmasked, except with diagonals of -1s
                # interpolated to match the effects of masking and/or input
                # nans.
                '2_masked': permuted_array([[-1, 0,  0,  0,  1,  1,  1],
                                            [0, -1,  0,  0,  1,  1,  1],
                                            [0,  0, -1,  0,  1,  1,  1],
                                            [0,  0,  0, -1,  1,  1,  1],
                                            [0,  0,  0,  1, -1,  1,  1],
                                            [0,  0,  0,  1,  1, -1,  1],
                                            [0,  0,  0,  1,  1,  1, -1]]),
                '3_masked': permuted_array([[-1, 0,  0,  1,  1,  2,  2],
                                            [0, -1,  0,  1,  1,  2,  2],
                                            [0,  0, -1,  1,  1,  2,  2],
                                            [0,  0,  1, -1,  1,  2,  2],
                                            [0,  0,  1,  1, -1,  2,  2],
                                            [0,  0,  1,  1,  2, -1,  2],
                                            [0,  0,  1,  1,  2,  2, -1]]),
                '6_masked': permuted_array([[-1, 0,  1,  2,  3,  4,  5],
                                            [0, -1,  1,  2,  3,  4,  5],
                                            [0,  1, -1,  2,  3,  4,  5],
                                            [0,  1,  2, -1,  3,  4,  5],
                                            [0,  1,  2,  3, -1,  4,  5],
                                            [0,  1,  2,  3,  4, -1,  5],
                                            [0,  1,  2,  3,  4,  5, -1]]),
                '2_nans': permuted_array([[0,  0,  0,  1,  1,  1, -1],
                                          [0,  0,  0,  1,  1, -1,  1],
                                          [0,  0,  0,  1, -1,  1,  1],
                                          [0,  0,  0, -1,  1,  1,  1],
                                          [0,  0, -1,  0,  1,  1,  1],
                                          [0, -1,  0,  0,  1,  1,  1],
                                          [-1, 0,  0,  0,  1,  1,  1]]),
                '3_nans': permuted_array([[0,  0,  1,  1,  2,  2, -1],
                                          [0,  0,  1,  1,  2, -1,  2],
                                          [0,  0,  1,  1, -1,  2,  2],
                                          [0,  0,  1, -1,  1,  2,  2],
                                          [0,  0, -1,  1,  1,  2,  2],
                                          [0, -1,  0,  1,  1,  2,  2],
                                          [-1, 0,  0,  1,  1,  2,  2]]),
                '6_nans': permuted_array([[0,  1,  2,  3,  4,  5, -1],
                                          [0,  1,  2,  3,  4, -1,  5],
                                          [0,  1,  2,  3, -1,  4,  5],
                                          [0,  1,  2, -1,  3,  4,  5],
                                          [0,  1, -1,  2,  3,  4,  5],
                                          [0, -1,  1,  2,  3,  4,  5],
                                          [-1, 0,  1,  2,  3,  4,  5]]),
            },
            mask=self.build_mask(self.ones_mask(shape=shape)),
        )

Example 39

Project: pyspeckit Source File: collapse_gaussfit.py
def wrap_collapse_adaptive(filename,outprefix,redo='no',nsig=5,nrsig=2,doplot=True):
    """
    redo - if not equal to 'no', then...
    if collapse_gaussfit succeeded (to the extent that the .pysav files were written),
    but some part of the file writing or successive procedures failed, re-do those 
    procedures without redoing the whole collapse
    """
    fitsfile = pyfits.open(filename)
    dv,v0,p3 = fitsfile[0].header['CD3_3'],fitsfile[0].header['CRVAL3'],fitsfile[0].header['CRPIX3']
    dr,r0,p1 = fitsfile[0].header['CD1_1'],fitsfile[0].header['CRVAL1'],fitsfile[0].header['CRPIX1']
    dd,d0,p2 = fitsfile[0].header['CD2_2'],fitsfile[0].header['CRVAL2'],fitsfile[0].header['CRPIX2']
    xtora = lambda x: (x-p1+1)*dr+r0    # convert pixel coordinates to RA/Dec/Velocity
    ytodec = lambda y: (y-p2+1)*dd+d0
    vconv = lambda v: (v-p3+1)*dv+v0

    cube = fitsfile[0].data
    cube = where(numpy.isnan(cube),0,cube)

    if redo=='no':
        adaptB = asarray(adaptive_collapse_gaussfit(cube,axis=0,prefix=outprefix+'_triple',
            nsig=nsig,nrsig=nrsig,vconv=vconv,xtora=xtora,ytodec=ytodec,doplot=doplot))
        adaptB[numpy.isnan(adaptB)] = 0
        pickle.dump(adaptB,open('%s_adaptB.pysav' % outprefix,'w'))
    else:
        adaptB = pickle.load(open('%s_adaptB.pysav' % outprefix,'r'))
    db = adaptB
    gcd = double_gaussian(db[3],db[4],db[0],db[1],db[5],db[6])(indices(cube.shape)[0])
    fitsfile[0].data = gcd
    fitsfile.writeto('%s_adaptgausscube.fits' % outprefix,clobber=True)

    gcd[numpy.isnan(gcd)] = 0
    adaptResids = cube-gcd
    fitsfile[0].data = adaptResids
    fitsfile.writeto('%s_adaptgaussresids.fits' % outprefix,clobber=True)


    #adaptB[4] = (adaptB[4]-v0) / dv + p3-1
    #adaptB[3] = (adaptB[3]-v0) / dv + p3-1
    adaptB[4] = (adaptB[4]-p3+1) * dv + v0
    adaptB[3] = (adaptB[3]-p3+1) * dv + v0
    fitsfile[0].data = asarray(adaptB)
    fitsfile.writeto('%s_adaptgausspars.fits' % outprefix,clobber=True)

    fitsfile[0].header.__delitem__('CD3_3')
    fitsfile[0].header.__delitem__('CRVAL3')
    fitsfile[0].header.__delitem__('CRPIX3')
    fitsfile[0].header.__delitem__('CUNIT3')
    fitsfile[0].header.__delitem__('CTYPE3')

    adaptResids[numpy.isnan(adaptResids)] = 0
    totalDResids = adaptResids.sum(axis=0)
    fitsfile[0].data = totalDResids
    fitsfile.writeto('%s_adaptgauss_totalresids.fits' % outprefix,clobber=True)

    return adaptB

Example 40

Project: zipline Source File: data_portal.py
Function: get_history_window
    def get_history_window(self, assets, end_dt, bar_count, frequency, field,
                           ffill=True):
        """
        Public API method that returns a dataframe containing the requested
        history window.  Data is fully adjusted.

        Parameters
        ----------
        assets : list of zipline.data.Asset objects
            The assets whose data is desired.

        bar_count: int
            The number of bars desired.

        frequency: string
            "1d" or "1m"

        field: string
            The desired field of the asset.

        ffill: boolean
            Forward-fill missing values. Only has effect if field
            is 'price'.

        Returns
        -------
        A dataframe containing the requested data.
        """
        if field not in OHLCVP_FIELDS and field != 'sid':
            raise ValueError("Invalid field: {0}".format(field))

        if frequency == "1d":
            if field == "price":
                df = self._get_history_daily_window(assets, end_dt, bar_count,
                                                    "close")
            else:
                df = self._get_history_daily_window(assets, end_dt, bar_count,
                                                    field)
        elif frequency == "1m":
            if field == "price":
                df = self._get_history_minute_window(assets, end_dt, bar_count,
                                                     "close")
            else:
                df = self._get_history_minute_window(assets, end_dt, bar_count,
                                                     field)
        else:
            raise ValueError("Invalid frequency: {0}".format(frequency))

        # forward-fill price
        if field == "price":
            if frequency == "1m":
                data_frequency = 'minute'
            elif frequency == "1d":
                data_frequency = 'daily'
            else:
                raise Exception(
                    "Only 1d and 1m are supported for forward-filling.")

            assets_with_leading_nan = np.where(isnull(df.iloc[0]))[0]

            history_start, history_end = df.index[[0, -1]]
            initial_values = []
            for asset in df.columns[assets_with_leading_nan]:
                last_traded = self.get_last_traded_dt(
                    asset,
                    history_start,
                    data_frequency,
                )
                if isnull(last_traded):
                    initial_values.append(nan)
                else:
                    initial_values.append(
                        self.get_adjusted_value(
                            asset,
                            field,
                            dt=last_traded,
                            perspective_dt=history_end,
                            data_frequency=data_frequency,
                        )
                    )

            # Set leading values for assets that were missing data, then ffill.
            df.ix[0, assets_with_leading_nan] = np.array(
                initial_values,
                dtype=np.float64
            )
            df.fillna(method='ffill', inplace=True)

            # forward-filling will incorrectly produce values after the end of
            # an asset's lifetime, so write NaNs back over the asset's
            # end_date.
            normed_index = df.index.normalize()
            for asset in df.columns:
                if history_end >= asset.end_date:
                    # if the window extends past the asset's end date, set
                    # all post-end-date values to NaN in that asset's series
                    df.loc[normed_index > asset.end_date, asset] = nan
        return df

Example 41

Project: sfepy Source File: equations.py
    def get_graph_conns(self, any_dof_conn=False, rdcs=None, cdcs=None,
                        active_only=True):
        """
        Get DOF connectivities needed for creating tangent matrix graph.

        Parameters
        ----------
        any_dof_conn : bool
            By default, only volume DOF connectivities are used, with
            the exception of trace surface DOF connectivities. If True,
            any kind of DOF connectivities is allowed.
        rdcs, cdcs : arrays, optional
            Additional row and column DOF connectivities, corresponding
            to the variables used in the equations.
        active_only : bool
            If True, the active DOF connectivities have reduced size and are
            created with the reduced (active DOFs only) numbering.

        Returns
        -------
        rdcs, cdcs : arrays
            The row and column DOF connectivities defining the matrix
            graph blocks.
        """
        if rdcs is None:
            rdcs = []
            cdcs = []

        elif cdcs is None:
            cdcs = copy(rdcs)

        else:
            assert_(len(rdcs) == len(cdcs))
            if rdcs is cdcs: # Make sure the lists are not the same object.
                rdcs = copy(rdcs)

        adcs = self.variables.adof_conns

        # Only volume dof connectivities are used, with the exception of trace
        # surface dof connectivities.
        shared = set()
        for key, ii, info in iter_dict_of_lists(self.conn_info,
                                                return_keys=True):
            rvar, cvar = info.virtual, info.state
            if (rvar is None) or (cvar is None):
                continue

            is_surface = rvar.is_surface or cvar.is_surface

            dct = info.dc_type.type
            if not (dct in ('volume', 'scalar', 'custom') or is_surface
                    or info.is_trace or any_dof_conn):
                continue

            rreg_name = info.get_region_name(can_trace=False)
            creg_name = info.get_region_name()

            rname = rvar.get_primary_name()
            rkey = (rname, rreg_name, dct, False)
            ckey = (cvar.name, creg_name, dct, info.is_trace)

            dc_key = (rkey, ckey)

            if not dc_key in shared:
                rdc = adcs[rkey]
                cdc = adcs[ckey]
                if not active_only:
                    ii = nm.where(rdc < 0)
                    rdc = rdc.copy()
                    rdc[ii] = -1 - rdc[ii]

                    ii = nm.where(cdc < 0)
                    cdc = cdc.copy()
                    cdc[ii] = -1 - cdc[ii]

                rdcs.append(rdc)
                cdcs.append(cdc)

                shared.add(dc_key)

        return rdcs, cdcs

Example 42

Project: sfepy Source File: lcbc_operators.py
    def __init__(self, name, regions, dof_names, dof_map_fun, shift_fun,
                 variables, ts, functions):
        LCBCOperator.__init__(self, name, regions, dof_names, dof_map_fun,
                              variables, functions=functions)

        self.shift_fun = get_condition_value(shift_fun, functions,
                                             'LCBC', 'shift')

        mvar = variables[self.var_names[0]]
        svar = variables[self.var_names[1]]

        mfield = mvar.field
        sfield = svar.field

        nmaster = mfield.get_dofs_in_region(regions[0], merge=True)
        nslave = sfield.get_dofs_in_region(regions[1], merge=True)

        if nmaster.shape != nslave.shape:
            msg = 'shifted EPBC node list lengths do not match!\n(%s,\n %s)' %\
                  (nmaster, nslave)
            raise ValueError(msg)

        mcoor = mfield.get_coor(nmaster)
        scoor = sfield.get_coor(nslave)

        i1, i2 = self.dof_map_fun(mcoor, scoor)
        self.mdofs = expand_nodes_to_equations(nmaster[i1], dof_names[0],
                                               self.all_dof_names[0])
        self.sdofs = expand_nodes_to_equations(nslave[i2], dof_names[1],
                                               self.all_dof_names[1])

        self.shift = self.shift_fun(ts, scoor[i2], regions[1])

        meq = mvar.eq_map.eq[self.mdofs]
        seq = svar.eq_map.eq[self.sdofs]

        # Ignore DOFs with EBCs or EPBCs.
        mia = nm.where(meq >= 0)[0]
        sia = nm.where(seq >= 0)[0]

        ia = nm.intersect1d(mia, sia)

        meq = meq[ia]
        seq = seq[ia]

        num = len(ia)

        ones = nm.ones(num, dtype=nm.float64)
        n_dofs = [variables.adi.n_dof[ii] for ii in self.var_names]
        mtx = sp.coo_matrix((ones, (meq, seq)), shape=n_dofs)

        self.mtx = mtx.tocsr()

        self.rhs = self.shift.ravel()[ia]

        self.ameq = meq
        self.aseq = seq

        self.n_mdof = len(nm.unique(meq))
        self.n_sdof = len(nm.unique(seq))
        self.n_new_dof = 0

Example 43

Project: qutip Source File: permute.py
Function: permute
def _permute(Q, order):
    if Q.isket:
        dims, perm = _perm_inds(Q.dims[0], order)
        nzs = Q.data.nonzero()[0]
        wh = np.where(perm == nzs)[0]
        data = np.ones(len(wh), dtype=int)
        cols = perm[wh].T[0]
        perm_matrix = sp.coo_matrix((data, (wh, cols)),
                                    shape=(Q.shape[0], Q.shape[0]), dtype=int)
        perm_matrix = perm_matrix.tocsr()
        return perm_matrix * Q.data, Q.dims

    elif Q.isbra:
        dims, perm = _perm_inds(Q.dims[1], order)
        nzs = Q.data.nonzero()[1]
        wh = np.where(perm == nzs)[0]
        data = np.ones(len(wh), dtype=int)
        rows = perm[wh].T[0]
        perm_matrix = sp.coo_matrix((data, (rows, wh)),
                                    shape=(Q.shape[1], Q.shape[1]), dtype=int)
        perm_matrix = perm_matrix.tocsr()
        return Q.data * perm_matrix, Q.dims

    elif Q.isoper:
        dims, perm = _perm_inds(Q.dims[0], order)
        data = np.ones(Q.shape[0], dtype=int)
        rows = np.arange(Q.shape[0], dtype=int)
        perm_matrix = sp.coo_matrix((data, (rows, perm.T[0])),
                                    shape=(Q.shape[1], Q.shape[1]), dtype=int)
        perm_matrix = perm_matrix.tocsr()
        dims_part = list(dims[order])
        dims = [dims_part, dims_part]
        return (perm_matrix * Q.data) * perm_matrix.T, dims

    elif Q.issuper or Q.isoperket:
        # For superoperators, we expect order to be something like
        # [[0, 2], [1, 3]], which tells us to permute according to
        # [0, 2, 1 ,3], and then group indices according to the length
        # of each sublist.
        # As another example,
        # permuting [[[1, 2, 3], [1, 2, 3]], [[1, 2, 3], [1, 2, 3]]] by
        # [[0, 3], [1, 4], [2, 5]] should give
        # [[[1, 1], [2, 2], [3, 3]], [[1, 1], [2, 2], [3, 3]]].
        #
        # Get the breakout of the left index into dims.
        # Since this is a super, the left index itself breaks into left
        # and right indices, each of which breaks down further.
        # The best way to deal with that here is to flatten dims.
        flat_order = sum(order, [])
        q_dims_left = sum(Q.dims[0], [])
        dims, perm = _perm_inds(q_dims_left, flat_order)
        dims = dims.flatten()

        data = np.ones(Q.shape[0], dtype=int)
        rows = np.arange(Q.shape[0], dtype=int)

        perm_matrix = sp.coo_matrix((data, (rows, perm.T[0])),
                                    shape=(Q.shape[0], Q.shape[0]), dtype=int)

        perm_matrix = perm_matrix.tocsr()
        dims_part = list(dims[flat_order])

        # Finally, we need to restructure the now-decomposed left index
        # into left and right subindices, so that the overall dims we return
        # are of the form specified by order.
        dims_part = list(_chunk_dims(dims_part, order))
        perm_left = (perm_matrix * Q.data)
        if Q.type == 'operator-ket':
            return perm_left, [dims_part, [1]]
        elif Q.type == 'super':
            return perm_left * perm_matrix.T, [dims_part, dims_part]
    else:
        raise TypeError('Invalid quantum object for permutation.')

Example 44

Project: SHARPpy Source File: map.py
Function: load_dat
    def _loadDat(self, name, res):
        """
        Code shamelessly lifted from Basemap's data file parser by Jeff Whitaker.
        http://matplotlib.org/basemap/
        """
        def segmentPath(b, lb_lat, ub_lat):
            paths = []
            if b[:, 1].max() <= lb_lat or b[:, 1].min() >= ub_lat:
                return paths

            idxs = np.where((b[:, 1] >= lb_lat) & (b[:, 1] <= ub_lat))[0]
            if len(idxs) < 2:
                return paths

            segs = (np.diff(idxs) == 1)
            try:
                breaks = np.where(segs == 0)[0] + 1
            except IndexError:
                breaks = []
            breaks = [ 0 ] + list(breaks) + [ -1 ]
 
            for idx in xrange(len(breaks) - 1):
                if breaks[idx + 1] == -1:
                    seg_idxs = idxs[breaks[idx]:]
                else:
                    seg_idxs = idxs[breaks[idx]:breaks[idx + 1]]

                if len(seg_idxs) >= 2:
                    paths.append(b[seg_idxs, ::-1])

            return paths

        bdatfile = open(os.path.join(Mapper.data_dir, name + '_' + res + '.dat'), 'rb')
        bdatmetafile = open(os.path.join(Mapper.data_dir, name + 'meta_' + res + '.dat'), 'r')

        projs = ['npstere', 'merc', 'spstere']
        paths = dict( (p, []) for p in projs )

#       old_proj = self.proj
        for line in bdatmetafile:
            lats, lons = [], []
            linesplit = line.split()
            area = float(linesplit[1])
            south = float(linesplit[3])
            north = float(linesplit[4])
            if area < 0:
                area = 1e30

            if area > 1500.:
                typ = int(linesplit[0])
                npts = int(linesplit[2])
                offsetbytes = int(linesplit[5])
                bytecount = int(linesplit[6])
                bdatfile.seek(offsetbytes,0)
                # read in binary string convert into an npts by 2
                # numpy array (first column is lons, second is lats).
                polystring = bdatfile.read(bytecount)
                # binary data is little endian.
                b = np.array(np.fromstring(polystring,dtype='<f4'),'f8')
                b.shape = (npts, 2)

                if np.any(b[:, 0] > 180):
                    b[:, 0] -= 360

                for proj in projs:
                    lb_lat, ub_lat = Mapper.min_lat[proj], Mapper.max_lat[proj]
                    path = segmentPath(b, lb_lat, ub_lat)
                    paths[proj].extend(path)

        return paths

Example 45

Project: GPy Source File: gradient_checker.py
    def checkgrad_block(self, analytic_hess, numeric_hess, verbose=False, step=1e-6, tolerance=1e-3, block_indices=None, plot=False):
        """
        Checkgrad a block matrix
        """
        if analytic_hess.dtype is np.dtype('object'):
            #Make numeric hessian also into a block matrix
            real_size = get_block_shapes(analytic_hess)
            num_elements = np.sum(real_size)
            if (num_elements, num_elements) == numeric_hess.shape:
                #If the sizes are the same we assume they are the same
                #(we have not fixed any values so the numeric is the whole hessian)
                numeric_hess = get_blocks(numeric_hess, real_size)
            else:
                #Make a fake empty matrix and fill out the correct block
                tmp_numeric_hess = get_blocks(np.zeros((num_elements, num_elements)), real_size)
                tmp_numeric_hess[block_indices] = numeric_hess.copy()
                numeric_hess = tmp_numeric_hess

        if block_indices is not None:
            #Extract the right block
            analytic_hess = analytic_hess[block_indices]
            numeric_hess = numeric_hess[block_indices]
        else:
            #Unblock them if they are in blocks and you aren't checking a single block (checking whole hessian)
            if analytic_hess.dtype is np.dtype('object'):
                analytic_hess = unblock(analytic_hess)
                numeric_hess = unblock(numeric_hess)

        ratio = numeric_hess / (numpy.where(analytic_hess==0, 1e-10, analytic_hess))
        difference = numpy.abs(analytic_hess - numeric_hess)

        check_passed = numpy.all((numpy.abs(1 - ratio)) < tolerance) or numpy.allclose(numeric_hess, analytic_hess, atol = tolerance)

        if verbose:
            if block_indices:
                print("\nBlock {}".format(block_indices))
            else:
                print("\nAll blocks")

            header = ['Checked', 'Max-Ratio', 'Min-Ratio', 'Min-Difference', 'Max-Difference']
            header_string = map(lambda x: ' | '.join(header), [header])
            separator = '-' * len(header_string[0])
            print('\n'.join([header_string[0], separator]))
            min_r = '%.6f' % float(numpy.min(ratio))
            max_r = '%.6f' % float(numpy.max(ratio))
            max_d = '%.6f' % float(numpy.max(difference))
            min_d = '%.6f' % float(numpy.min(difference))
            cols = [max_r, min_r, min_d, max_d]

            if check_passed:
                checked = "\033[92m  True \033[0m"
            else:
                checked = "\033[91m  False \033[0m"

            grad_string = "{} | {}  | {} |    {}    |   {} ".format(checked, cols[0], cols[1], cols[2], cols[3])
            print(grad_string)

            if plot:
                from matplotlib import pyplot as pb
                fig, axes = pb.subplots(2, 2)
                max_lim = numpy.max(numpy.vstack((analytic_hess, numeric_hess)))
                min_lim = numpy.min(numpy.vstack((analytic_hess, numeric_hess)))
                msa = axes[0,0].matshow(analytic_hess, vmin=min_lim, vmax=max_lim)
                axes[0,0].set_title('Analytic hessian')
                axes[0,0].xaxis.set_ticklabels([None])
                axes[0,0].yaxis.set_ticklabels([None])
                axes[0,0].xaxis.set_ticks([None])
                axes[0,0].yaxis.set_ticks([None])
                msn = axes[0,1].matshow(numeric_hess, vmin=min_lim, vmax=max_lim)
                pb.colorbar(msn, ax=axes[0,1])
                axes[0,1].set_title('Numeric hessian')
                axes[0,1].xaxis.set_ticklabels([None])
                axes[0,1].yaxis.set_ticklabels([None])
                axes[0,1].xaxis.set_ticks([None])
                axes[0,1].yaxis.set_ticks([None])
                msr = axes[1,0].matshow(ratio)
                pb.colorbar(msr, ax=axes[1,0])
                axes[1,0].set_title('Ratio')
                axes[1,0].xaxis.set_ticklabels([None])
                axes[1,0].yaxis.set_ticklabels([None])
                axes[1,0].xaxis.set_ticks([None])
                axes[1,0].yaxis.set_ticks([None])
                msd = axes[1,1].matshow(difference)
                pb.colorbar(msd, ax=axes[1,1])
                axes[1,1].set_title('difference')
                axes[1,1].xaxis.set_ticklabels([None])
                axes[1,1].yaxis.set_ticklabels([None])
                axes[1,1].xaxis.set_ticks([None])
                axes[1,1].yaxis.set_ticks([None])
                if block_indices:
                    fig.suptitle("Block: {}".format(block_indices))
                pb.show()

        return check_passed

Example 46

Project: pyunicorn Source File: recurrence_network.py
Function: init
    def __init__(self, time_series, metric="supremum", normalize=False,
                 missing_values=False, silence_level=0, **kwds):
        """
        Initialize an instance of RecurrenceNetwork.

        Creates an embedding of the given time series, calculates a recurrence
        plot from the embedding and then creates a Network object from the
        recurrence plot, interpreting the recurrence matrix as the adjacency
        matrix of a complex network.

        Either recurrence threshold ``threshold``/``threshold_std``, recurrence
        rate ``recurrence_rate`` or local recurrence rate
        ``local_recurrence_rate`` have to be given as keyword arguments.

        Embedding is only supported for scalar time series. If embedding
        dimension ``dim`` and delay ``tau`` are **both** given as keyword
        arguments, embedding is applied. Multidimensional time series are
        processed as is by default.

        :type time_series: 2D array (time, dimension)
        :arg time_series: The time series to be analyzed, can be scalar or
            multi-dimensional.
        :arg str metric: The metric for measuring distances in phase space
            ("manhattan", "euclidean", "supremum").
        :arg bool normalize: Decide whether to normalize the time series to
            zero mean and unit standard deviation.
        :arg bool missing_values: Toggle special treatment of missing values in
            :attr:`.RecurrencePlot.time_series`.
        :arg number silence_level: Inverse level of verbosity of the object.
        :arg number threshold: The recurrence threshold keyword for generating
            the recurrence network using a fixed threshold.
        :arg number threshold_std: The recurrence threshold keyword for
            generating the recurrence plot using a fixed threshold in units of
            the time series' STD.
        :arg number recurrence_rate: The recurrence rate keyword for generating
            the recurrence network using a fixed recurrence rate.
        :arg number local_recurrence_rate: The local recurrence rate keyword
            for generating the recurrence plot using a fixed local recurrence
            rate (same number of recurrences for each state vector).
        :arg number adaptive_neighborhood_size: The adaptive neighborhood size
            parameter for generating recurrence plots based on the algorithm in
            [Xu2008]_.
        :arg number dim: The embedding dimension.
        :arg number tau: The embedding delay.
        :type node_weights: 1D array (time)
        :arg node_weights: The sequence of weights associated with each
            node for calculating n.s.i. network measures.
        """
        #  Initialize the underlying RecurrencePlot object
        RecurrencePlot.__init__(self, time_series=time_series, metric=metric,
                                normalize=normalize,
                                missing_values=missing_values,
                                silence_level=silence_level, **kwds)

        #  Set diagonal of R to zero to avoid self-loops in the recurrence
        #  network
        A = self.R.copy()
        A.flat[::self.N+1] = 0

        #  Get keywords
        local_recurrence_rate = kwds.get("local_recurrence_rate")
        node_weights = kwds.get("node_weights")

        #  Assign to each embedded state vector the weight of the earliest
        #  sample contained in it (too simple!)
        if node_weights is not None:
            node_weights = node_weights[:self.N]

        #  Remove state vectors containing missing values
        if missing_values:
            A = np.delete(A, np.where(self.missing_value_indices), axis=0)
            A = np.delete(A, np.where(self.missing_value_indices), axis=1)

        #  Create a Network object interpreting the recurrence matrix as the
        #  graph adjacency matrix.
        if local_recurrence_rate is not None:
            #  A recurrence network with fixed local recurrence rate (Eckmann
            #  definition of a recurrence plot) is directed by definition.
            Network.__init__(self, A, directed=True, node_weights=node_weights,
                             silence_level=silence_level)
        else:
            #  Create a Network object interpreting the recurrence matrix as
            #  the graph adjacency matrix. Recurrence networks are undirected
            #  by definition.
            Network.__init__(self, A, directed=False,
                             node_weights=node_weights,
                             silence_level=silence_level)

Example 47

Project: ComputationalNeurodynamics Source File: QIFNetwork.py
Function: neuronupdate
  def NeuronUpdate(self, i, t):
    """
    QIF neuron update function. Update one layer for 1 millisecond
    using the Euler method.

    Inputs:
    i -- Number of layer to update
    t -- Current timestep. Necessary to sort out the synaptic delays.
    """

    # Euler method step size in ms
    dt = 0.2

    # Calculate current from incoming spikes
    for j in xrange(self.Nlayers):

      # If layer[i].S[j] exists then layer[i].factor[j] and
      # layer[i].delay[j] have to exist
      if j in self.layer[i].S:
        S = self.layer[i].S[j]  # target neuron->rows, source neuron->columns

        # Firings contains time and neuron idx of each spike.
        # [t, index of the neuron in the layer j]
        firings = self.layer[j].firings

        # Find incoming spikes taking delays into account
        delay = self.layer[i].delay[j]
        F = self.layer[i].factor[j]

        # Sum current from incoming spikes
        k = len(firings)
        while k > 0 and (firings[k-1, 0] > (t - self.Dmax)):
          idx = delay[:, firings[k-1, 1]] == (t-firings[k-1, 0])
          self.layer[i].I[idx] += F * S[idx, firings[k-1, 1]]
          k = k-1

    # Update v using the QIF equations and Euler method
    for k in xrange(int(1/dt)):
      v = self.layer[i].v

      self.layer[i].v += dt*(
          self.layer[i].a*(self.layer[i].vr - v)*(self.layer[i].vc - v) +
          self.layer[i].R*self.layer[i].I) / self.layer[i].tau

      # Find index of neurons that have fired this millisecond
      fired = np.where(self.layer[i].v >= 30)[0]

      if len(fired) > 0:
        for f in fired:
          # Add spikes into spike train
          if len(self.layer[i].firings) != 0:
            self.layer[i].firings = np.vstack([self.layer[i].firings, [t, f]])
          else:
            self.layer[i].firings = np.array([[t, f]])

          # Reset the membrane potential after spikes
          # Here's a little hack to see if vr is array or scalar
          if hasattr(self.layer[i].vr, "__len__"):
            self.layer[i].v[f]  = self.layer[i].vr[f]
          else:
            self.layer[i].v[f]  = self.layer[i].vr

    return

Example 48

Project: mpop Source File: HRWimage.py
def HRW_2dfield( HRW_data, obj_area, interpol_method=None, hrw_channels=None, min_correlation=None, level=''):

    print "... calculate 2d wind field (HRW_2dfield)"

    if min_correlation != None:
        print "    filter for min_correlation = ", min_correlation
        inds = where(HRW_data.correlation > min_correlation)
        HRW_data.subset(inds)

    xx, yy = obj_area.get_xy_from_lonlat( HRW_data.lon, HRW_data.lat, outside_error=False, return_int=False) #, return_int=True

    yy = obj_area.y_size - yy  

    uu = - HRW_data.wind_speed * sin(radians(HRW_data.wind_direction))
    vv = - HRW_data.wind_speed * cos(radians(HRW_data.wind_direction))

    # get rid of all vectors outside of the field 
    index = nonzero(xx)
    xx = xx[index]
    yy = yy[index]
    uu = uu[index]
    vv = vv[index]

    points = transpose(append([xx], [yy], axis=0))
    #print type(uu), uu.shape
    #print type(points), points.shape
    #print points[0], yy[0], xx[0]
    #print uu[0]

    nx = obj_area.x_size
    ny = obj_area.y_size

    x2 = arange(nx)
    y2 = (ny-1) - arange(ny)

    grid_x, grid_y = meshgrid(x2, y2)
    
    if interpol_method == None:
        # we need at least 2 winds to interpolate 
        if uu.size < 4:  
            print "*** Warning, not wnough wind data available, n_winds = ", uu.size
            fake = empty(grid_x.shape)
            fake[:,:] = nan
            HRW_data.interpol_method = None
            return fake, fake

        elif uu.size < 50:
            interpol_method = "RBF"

        else:
            interpol_method = "linear + nearest"
     
    #interpol_method = "nearest"
    #interpol_method = "cubic + nearest" # might cause unrealistic overshoots
    #interpol_method = "kriging"
    #interpol_method = "..."

    print "... min windspeed (org data): ", HRW_data.wind_speed.min()
    print "... max windspeed (org data): ", HRW_data.wind_speed.max()

    for i_iteration in [0,1]:

        if interpol_method == "nearest":

            print '... fill with nearest neighbour'
            # griddata, see http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.griddata.html
            grid_u1x = griddata(points, uu, (grid_x, grid_y), method='nearest')
            grid_v1x = griddata(points, vv, (grid_x, grid_y), method='nearest')

        elif interpol_method == "RBF":

            print '... inter- and extrapolation using radial basis functions'
            # https://www.youtube.com/watch?v=_cJLVhdj0j4
            print "... start Rbf"
            from scipy.interpolate import Rbf
            # rbfu = Rbf(xx, yy, uu, epsilon=0.1) #
            rbfu = Rbf(xx, yy, uu, epsilon=0.2)
            grid_u1x = rbfu(grid_x, grid_y)
            rbfv = Rbf(xx, yy, vv, epsilon=0.1) #
            grid_v1x = rbfv(grid_x, grid_y)
            print "... finish Rbf"
            # !very! slow for a large number of observations 

        elif interpol_method == "linear + nearest" or interpol_method == "cubic + nearest":

            if interpol_method == "linear + nearest":
                print '... calculate linear interpolation'
                # griddata, see http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.griddata.html
                grid_u1 = griddata(points, uu, (grid_x, grid_y), method='linear')
                grid_v1 = griddata(points, vv, (grid_x, grid_y), method='linear')
            elif interpol_method == "cubic + nearest":
                # smoother, but can cause unrealistic overshoots
                print '... calculate cubic interpolation'
                grid_u1 = griddata(points, uu, (grid_x, grid_y), method='cubic')
                grid_v1 = griddata(points, vv, (grid_x, grid_y), method='cubic')
            else:
                print "*** Error in mpop/imageo/HRWimage.py"
                print "    unknown interpolation method: ", interpol_method
                quit()

            if 1==1:
                # use faster function to extrapolate with closest neighbour
                print "... fill outside area with closest value"
                grid_u1x = fill_with_closest_pixel(grid_u1, invalid=None) 
                grid_v1x = fill_with_closest_pixel(grid_v1, invalid=None) 
            else:
                # use griddata to extrapolate with closest neighbour
                points2 = transpose(append([grid_x.flatten()], [grid_y.flatten()], axis=0))
                print type(grid_x.flatten()), grid_x.flatten().shape
                print type(points2), points2.shape
                mask = ~isnan(grid_v1.flatten())
                inds = where(mask)[0]
                grid_u1x = griddata(points2[inds], grid_u1.flatten()[inds], (grid_x, grid_y), method='nearest')
                grid_v1x = griddata(points2[inds], grid_v1.flatten()[inds], (grid_x, grid_y), method='nearest')

            if 1==0:
                # add othermost points as additional data
                y_add = [0,    0, ny-1, ny-1]
                x_add = [0, nx-1,    0, nx-1]
                for (i,j) in zip(x_add,y_add):
                    uu = append(uu, grid_u0[i,j])
                    vv = append(vv, grid_v0[i,j])
                xx = append(xx, x_add)
                yy = append(yy, y_add)
                points = transpose(append([yy], [xx], axis=0))

                print 'calc extent1'
                grid_u1e = griddata(points, uu, (grid_x, grid_y), method='linear')
                grid_v1e = griddata(points, vv, (grid_x, grid_y), method='linear')

        else:
            print "*** Error in mpop/imageo/HRWimage.py"
            print "    unknown interpol_method", interpol_method
            quit()

        ##http://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html
        ##http://stackoverflow.com/questions/3526514/problem-with-2d-interpolation-in-scipy-non-rectangular-grid
        #print "SmoothBivariateSpline:"
        #from scipy.interpolate import SmoothBivariateSpline
        #fitu = SmoothBivariateSpline( xx, yy, uu, s=1000) # , kx=3, ky=3, s = smooth * z2sum m 
        #from numpy import empty 
        #grid_u_SBP = empty(grid_x.shape)
        #for i in range(0,nx-1):       # starting upper right going down
        #    for j in range(0,ny-1):   # starting lower right going right
        #        #print i,j
        #        grid_u_SBP[j,i] = fitu(j,i)

        #grid_u_SBP = np.array([k.predict([x,y]) for x,y in zip(np.ravel(grid_x), np.ravel(grid_y))])
        #grid_u_SBP = grid_u_SBP.reshape(grid_x.shape)

        ##print x2
        ##print y2
        #grid_u_SBP = fitu(x2,y2)
        ##print "grid_u_SBP.shape", grid_u_SBP.shape
        ###print grid_u_SBP
        #print "End SmoothBivariateSpline:"

        #print "bisplrep:"
        #from scipy import interpolate
        #tck = interpolate.bisplrep(xx, yy, uu)
        #grid_u_BSR = interpolate.bisplev(grid_x[:,0], grid_y[0,:], tck)
        #print grid_u_BSR.shape
        #print "bisplrep"
        #print "grid_v1x.shape", grid_v1x.shape

        extent=(0,nx,0,ny)
        origin='lower'
        origin='upper'
        origin=None 

        # show different stages of 2d inter- and extra-polation 
        if 1==0:
            print 'make matplotlib.pyplot'
            import matplotlib.pyplot as plt
            vmin=-10
            vmax=10
            fig = plt.figure()
            plt.subplot(221)
            plt.title('u '+interpol_method)
            plt.plot(points[:,0], ny-1-points[:,1], 'k.', ms=1)
            plt.imshow(grid_u1x, vmin=vmin, vmax=vmax) #, extent=extent
            #plt.colorbar()
            plt.subplot(222)
            plt.title('v '+interpol_method)
            plt.plot(points[:,0], ny-1-points[:,1], 'k.', ms=1)
            plt.imshow(grid_v1x, origin=origin, vmin=vmin, vmax=vmax) #, extent=extent
            #plt.colorbar()

            # standard calculation for comparison 
            print '... calculate linear interpolation'
            grid_u1 = griddata(points, uu, (grid_x, grid_y), method='linear')
            grid_v1 = griddata(points, vv, (grid_x, grid_y), method='linear')
            grid_u1xx = fill_with_closest_pixel(grid_u1, invalid=None) 
            grid_v1xx = fill_with_closest_pixel(grid_v1, invalid=None) 

            plt.subplot(223)
            plt.title('U Linear+Nearest')
            plt.plot(points[:,0], ny-1-points[:,1], 'k.', ms=1)
            plt.imshow(grid_u1xx, origin=origin, vmin=vmin, vmax=vmax) #, extent=extent
            #plt.colorbar()
            plt.subplot(224)
            plt.title('V Linear+Nearest') 
            plt.plot(points[:,0], ny-1-points[:,1], 'k.', ms=1)
            #plt.title('Cubic')
            plt.imshow(grid_v1xx, origin=origin, vmin=vmin, vmax=vmax) #, extent=extent
            #plt.colorbar()
            plt.gcf().set_size_inches(6, 6)
            #plt.show()  # does not work with AGG
            tmpfile="test_hrw"+level+".png"
            fig.savefig(tmpfile)
            print "display "+tmpfile+" &"


        if grid_u1x.min() < -150 or grid_v1x.min() < -150 or grid_u1x.max() > 150 or grid_v1x.max() > 150:
            print "*** Warning, numerical instability detected, interpolation method: ", interpol_method
            print "    min u windspeed (u 2dimensional): ", grid_u1x.min()
            print "    min v windspeed (v 2dimensional): ", grid_v1x.min()
            print "    max u windspeed (u 2dimensional): ", grid_u1x.max()
            print "    max v windspeed (v 2dimensional): ", grid_v1x.max()
            interpol_method = "glinear + nearest"
            print "... try another interpolation method: ", interpol_method
        else:
            # (hopefully) numerical stable interpolation, exit the interpolation loop
            break 

    HRW_data.interpol_method = interpol_method

    return grid_u1x, grid_v1x

Example 49

Project: pyoptools Source File: misc.py
def interpolate_g(xi,yi,zi,xx,yy,knots=10, error=False,mask=None):
    """Create a grid of zi values interpolating the values from xi,yi,zi
    
    xi,yi,zi  1D Lists  or arrays containing the values to use as base for the interpolation
    xx,yy     1D vectors or lists containing the output coordinates
    samples  tuple containing the shape of the output array.
    knots     number of knots to be used in each direction
    error     if set to true, half of the points (x, y, z) are used to create 
              the interpolation, and half are used to evaluate the interpolation error
    
    
    """
    xi=array(xi)
    yi=array(yi)
    zi=array(zi)
    
    #print xi
    #print yi
    #print zi
    assert xi.ndim==1 ,"xi must ba a 1D array or list" 
    assert yi.ndim==1 ,"yi must ba a 1D array or list" 
    assert zi.ndim==1 ,"zi must ba a 1D array or list" 
    
    assert xx.ndim==1 ,"xx must ba a 1D array or list" 
    assert yy.ndim==1 ,"yy must ba a 1D array or list" 
    
    assert len(xi)==len(yi) and len(xi)==len(zi), "xi, yi, zi must have the same number of items"
    
    
    
    if error==True:
        # Create a list of indexes to be able to select the points that are going 
        # to be used as spline generators, and as control points
        idx=where(arange(len(xi)) %2 ==0, False, True)

   # Use only half of the samples to create the Spline, 
    if error == True:
        isp=argwhere(idx==True)
        ich=argwhere(idx==False)
        
        xsp=xi[isp]
        ysp=yi[isp]
        zsp=zi[isp]
        
        xch=xi[ich]
        ych=yi[ich]
        zch=zi[ich]
    else:
        xsp=xi
        ysp=yi
        zsp=zi
    
    #Distribute humogeneously the knots
    xk=linspace(xsp.min(), xsp.max(),knots)
    yk=linspace(ysp.min(), ysp.max(),knots)
    
    # LSQBivariateSpline using some knots gives smaller error than 
    # SmoothBivariateSpline
    di=interpolate.LSQBivariateSpline(xsp, ysp, zsp, xk[1:-1],  yk[1:-1])
    #print xsp,ysp,zsp
    #di=interpolate.SmoothBivariateSpline(xsp, ysp, zsp)
        
    # Evaluate error
    if error==True:
        zch1=di.ev(xch, ych)
        er=(zch.flatten()-zch1).std()

    if mask==None:
        #d=griddata(xi,  yi,  zi,  xx, yy) #
        d=di(xx,yy).transpose()
    else:
        d=ma_array(di(xx,yy).transpose(), mask=mask)
    
    if error==True: return d, er
    else: return d

Example 50

Project: Chimp Source File: replay_memory.py
Function: minibatch
    def minibatch(self):
        ''' Uniformly sample (o,a,r,o') experiences from the replay dataset.

        Args:
            batch_size: size of mini-batch

        Returns:
            Five numpy arrays that corresponds to o, a, r, o', and the terminal
            state indicator.
        '''
        batch_size = self.batch_size
        if batch_size >= self.valid:
            raise ValueError("Can't draw sample of size %d from replay dataset of size %d"
                       % (batch_size, self.valid))

        ohist_size, ahist_size, rhist_size = self.ohist_size, self.ahist_size, self.rhist_size
        max_hist = self.max_size

        indices = self.get_indices(batch_size)

        self.clear_history()

        # TODO: can we get rid of this loop by sorting inidces and then reshaping? 
        for i in xrange(batch_size):
            # all end on the same index
            endi = indices[i]
            starti = endi - max_hist
            # starting indecies if no terminal states
            starto, starta, startr = endi-ohist_size, endi-ahist_size, endi-rhist_size

            # look backwards and find first terminal state
            termarr = np.where(self.terminals[starti:endi-1]==True)[0]
            termidx = starti
            if termarr.size is not 0:
                termidx = endi - (endi-starti - termarr[-1]) + 1

            # if starts before terminal, change start index
            starto = termidx if starto < termidx else starto
            starta = termidx if starta < termidx else starta
            startr = termidx if startr < termidx else startr

            ohl, ahl, rhl = (endi - starto), (endi - starta), (endi - startr)
        
            # load from memory
            self.ohist[i, ohist_size-ohl:] = self.observations[xrange(starto, endi)]
            self.ophist[i, ohist_size-ohl:] = self.next_observations[xrange(starto, endi)]
            self.ahist[i, ahist_size-ahl:] = self.actions[xrange(starta, endi)]
            self.rhist[i, rhist_size-rhl:] = self.rewards[xrange(startr, endi)]
            self.thist[i, ohist_size-ohl:] = self.terminals[xrange(starto, endi)]

        return self.ohist, self.ahist, self.rhist, self.ophist, self.thist
See More Examples - Go to Next Page
Page 1 Selected Page 2 Page 3 Page 4