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
2
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
2
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
2
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
2
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]
2
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
2
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
2
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)
0
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')
0
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
0
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')
0
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
0
Example 12
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)
0
Example 13
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
0
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)
0
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),
)
0
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')
0
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
0
Example 18
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
0
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'))
0
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
0
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)
0
Example 22
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
0
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))
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)
0
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
0
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
0
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[:]
0
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
0
Example 29
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()
0
Example 30
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
0
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
0
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
0
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
0
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
0
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
0
Example 36
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
0
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')
0
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)),
)
0
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
0
Example 40
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
0
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
0
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
0
Example 43
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.')
0
Example 44
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
0
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
0
Example 46
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)
0
Example 47
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
0
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
0
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
0
Example 50
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