def grid_search(self):
C_range = np.logspace(-5, 15, 21, base=2)
param_grid = dict(C=C_range)
cv = StratifiedShuffleSplit(self.y_ex, n_iter=5, test_size=0.2, random_state=42)
grid = GridSearchCV(SVC(kernel='poly', max_iter=10000), param_grid=param_grid, cv=cv, n_jobs=1, verbose=0)
logger.info('start grid search for Linear')
grid.fit(self.X_ex, self.y_ex)
logger.info('end grid search for Linear')
scores = [x[1] for x in grid.grid_scores_]
# final train
clf = grid.best_estimator_
pred_train = clf.predict(self.X_ex)
pred_val = clf.predict(self.val_x)
pred_test = clf.predict(self.test_x)
r = Result(self.name + ' (X)', 'Poly', len(self.X_ex),
sm.accuracy_score(self.y_ex, pred_train),
sm.accuracy_score(self.val_y, pred_val),
sm.accuracy_score(self.test_y, pred_test))
return r
python类logspace()的实例源码
def wage_data_linear():
X, y = wage()
gam = LinearGAM(n_splines=10)
gam.gridsearch(X, y, lam=np.logspace(-5,3,50))
XX = generate_X_grid(gam)
plt.figure()
fig, axs = plt.subplots(1,3)
titles = ['year', 'age', 'education']
for i, ax in enumerate(axs):
ax.plot(XX[:, i], gam.partial_dependence(XX, feature=i+1))
ax.plot(XX[:, i], *gam.partial_dependence(XX, feature=i+1, width=.95)[1],
c='r', ls='--')
if i == 0:
ax.set_ylim(-30,30);
ax.set_title(titles[i])
fig.tight_layout()
plt.savefig('imgs/pygam_wage_data_linear.png', dpi=300)
def plot_pdf_log2(x, nbins=10, **kwargs):
'''
Adds a log-log PDF plot to the current axes. The PDF is binned with
logarithmic binning of base 2.
Arguments
---------
x : array_like
The data to plot
nbins : integer
The number of bins to take
Additional keyword arguments are passed to `matplotlib.pyplot.loglog`.
'''
x = np.asarray(x)
exp_max = np.ceil(np.log2(x.max()))
bins = np.logspace(0, exp_max, exp_max + 1, base=2)
ax = plt.gca()
hist, _ = np.histogram(x, bins=bins)
binsize = np.diff(np.asfarray(bins))
hist = hist / binsize
ax.loglog(bins[1:], hist, 'ow', **kwargs)
return ax
def get_model():
if FLAGS.model == 'logistic':
return linear_model.LogisticRegressionCV(class_weight='balanced',
scoring='roc_auc',
n_jobs=FLAGS.n_jobs,
max_iter=10000, verbose=1)
elif FLAGS.model == 'random_forest':
return ensemble.RandomForestClassifier(n_estimators=100,
n_jobs=FLAGS.n_jobs,
class_weight='balanced',
verbose=1)
elif FLAGS.model == 'svm':
return grid_search.GridSearchCV(
estimator=svm.SVC(kernel='rbf', gamma='auto',
class_weight='balanced'),
param_grid={'C': np.logspace(-4, 4, 10)}, scoring='roc_auc',
n_jobs=FLAGS.n_jobs, verbose=1)
else:
raise ValueError('Unrecognized model %s' % FLAGS.model)
def __init__(self, model, ax=None, alphas=None,
cv=None, scoring=None, **kwargs):
# Check to make sure this is not a "RegressorCV"
name = model.__class__.__name__
if name.endswith("CV"):
raise YellowbrickTypeError((
"'{}' is a CV regularization model;"
" try AlphaSelection instead."
).format(name))
# Call super to initialize the class
super(ManualAlphaSelection, self).__init__(model, ax=ax, **kwargs)
# Set manual alpha selection parameters
self.alphas = alphas or np.logspace(-10, -2, 200)
self.errors = None
self.score_method = partial(cross_val_score, cv=cv, scoring=scoring)
def create_grid(self,abins=None,zbins=None):
if abins is None and zbins is None:
filenames = glob.glob(self.get_dirname()+'/%s_*.dat'%(self._prefix))
data = np.array([self.filename2params(f) for f in filenames])
if not len(data):
msg = "No isochrone files found in: %s"%self.get_dirname()
raise Exception(msg)
arange = np.unique(data[:,0])
zrange = np.unique(data[:,1])
elif abins is not None and zbins is not None:
# Age in units of Gyr
arange = np.linspace(abins[0],abins[1],abins[2]+1)
# Metallicity sampled logarithmically
zrange = np.logspace(np.log10(zbins[0]),np.log10(zbins[1]),zbins[2]+1)
else:
msg = "Must specify both `abins` and `zbins` or neither"
raise Exception(msg)
aa,zz = np.meshgrid(arange,zrange)
return aa.flatten(),zz.flatten()
def search_queue_params():
df = []
data_batch_sizes = np.logspace(0, 8, num=9, base=2, dtype=int)
capacities = np.logspace(0, 12, num=13, base=2, dtype=int)
nthreads = np.logspace(0, 5, num=6, base=2, dtype=int)
for nth in nthreads:
for data_batch_size in data_batch_sizes:
for capacity in capacities:
cap = nth * capacity
tf.reset_default_graph()
d = DataHDF5(batch_size=data_batch_size)
queue = data.Queue(d.node, d,
queue_type='fifo',
batch_size=BATCH_SIZE,
capacity=cap,
n_threads=nth)
queue.kind = '{} / {} / {}'.format(nth, data_batch_size, capacity)
durs = time_tf(queue)
durs['data batch size'] = data_batch_size
durs['queue capacity'] = cap
durs['nthreads'] = nth
df.append(durs)
d.cleanup()
df = pandas.concat(df, ignore_index=True)
df.kind = df.kind.astype('category', ordered=True, categories=df.kind.unique())
df.to_pickle('/home/qbilius/mh17/computed/search_queue_params.pkl')
print(df.groupby(['nthreads', 'data batch size', 'queue capacity']).dur.mean())
def scaled_histogram(data, num_points, scale):
if scale == 'linear':
hist, edges = np.histogram(data, bins=max([5, int(num_points / 50)]))
else:
# Conditional catches an empty data input.
h1, h2 = ((np.log10(min(data)), np.log10(max(data))) if len(data) > 0 else (0, 1))
hist, edges = np.histogram(data, bins=np.logspace(h1, h2, 1 + max([5, int(num_points / 50)])))
hist_max = max(hist) * 1.1
return hist, edges, hist_max
def bernoulli_gaussian_trial(M=250,N=500,L=1000,pnz=.1,kappa=None,SNR=40):
A = np.random.normal(size=(M, N), scale=1.0 / math.sqrt(M)).astype(np.float32)
if kappa >= 1:
# create a random operator with a specific condition number
U,_,V = la.svd(A,full_matrices=False)
s = np.logspace( 0, np.log10( 1/kappa),M)
A = np.dot( U*(s*np.sqrt(N)/la.norm(s)),V).astype(np.float32)
A_ = tf.constant(A,name='A')
prob = TFGenerator(A=A,A_=A_,pnz=pnz,kappa=kappa,SNR=SNR)
prob.name = 'Bernoulli-Gaussian, random A'
bernoulli_ = tf.to_float( tf.random_uniform( (N,L) ) < pnz)
xgen_ = bernoulli_ * tf.random_normal( (N,L) )
noise_var = pnz*N/M * math.pow(10., -SNR / 10.)
ygen_ = tf.matmul( A_,xgen_) + tf.random_normal( (M,L),stddev=math.sqrt( noise_var ) )
prob.xval = ((np.random.uniform( 0,1,(N,L))<pnz) * np.random.normal(0,1,(N,L))).astype(np.float32)
prob.yval = np.matmul(A,prob.xval) + np.random.normal(0,math.sqrt( noise_var ),(M,L))
prob.xinit = ((np.random.uniform( 0,1,(N,L))<pnz) * np.random.normal(0,1,(N,L))).astype(np.float32)
prob.yinit = np.matmul(A,prob.xinit) + np.random.normal(0,math.sqrt( noise_var ),(M,L))
prob.xgen_ = xgen_
prob.ygen_ = ygen_
prob.noise_var = noise_var
return prob
def test_gradE(adjcube):
"""Tests the gradient of `E` using finite difference methods.
"""
from pydft.bases.fourier import gradE, E
from numpy.matlib import randn
cell = adjcube
V = QHO(cell)
Ns=4
#He set the random seed; we could do the same, but the
#implementation is probably different between numpy and matlab:
#randn('seed', 0.2004)
W = np.array(randn(np.prod(cell.S), Ns) + 1j*randn(np.prod(cell.S), Ns))
# Compute intial energy and gradient
E0 = E(V, W, cell)
g0 = gradE(V, W, cell)
# Choose a random direction to explore
dW = np.array(randn(W.shape) + 1j*randn(W.shape))
# Explore a range of step sizes decreasing by powers of ten
steps = np.logspace(np.log10(1e-3), np.log10(1e-7), 8)
for delta in steps:
# Directional derivative formula
dE = 2*np.real(np.trace(np.dot(g0.conjugate().T, delta*dW)))
# Print ratio of actual change to expected change, along with estimate
# of the error in this quantity due to rounding
ratio = abs(1.-(E(V, W+delta*dW, cell)-E0)/dE)
print(int(np.log10(ratio)), int(np.log10(delta)), ratio)
assert abs(int(np.log10(ratio)) - int(np.log10(delta))) <= 2
def fcn_FDEM_InductionSpherePlaneWidget(xtx,ytx,ztx,m,orient,x0,y0,z0,a,sig,mur,xrx,yrx,zrx,logf,Comp,Phase):
sig = 10**sig
f = 10**logf
fvec = np.logspace(0,8,41)
xmin, xmax, dx, ymin, ymax, dy = -30., 30., 0.3, -30., 30., 0.4
X,Y = np.mgrid[xmin:xmax+dx:dx, ymin:ymax+dy:dy]
X = np.transpose(X)
Y = np.transpose(Y)
Obj = SphereFEM(m,orient,xtx,ytx,ztx)
Hx,Hy,Hz,Habs = Obj.fcn_ComputeFrequencyResponse(f,sig,mur,a,x0,y0,z0,X,Y,zrx)
Hxi,Hyi,Hzi,Habsi = Obj.fcn_ComputeFrequencyResponse(fvec,sig,mur,a,x0,y0,z0,xrx,yrx,zrx)
fig1 = plt.figure(figsize=(17,6))
Ax1 = fig1.add_axes([0.04,0,0.43,1])
Ax2 = fig1.add_axes([0.6,0,0.4,1])
if Comp == 'x':
Ax1 = plotAnomalyXYplane(Ax1,f,X,Y,ztx,Hx,Comp,Phase)
Ax1 = plotPlaceTxRxSphereXY(Ax1,xtx,ytx,xrx,yrx,x0,y0,a)
Ax2 = plotResponseFEM(Ax2,f,fvec,Hxi,Comp)
elif Comp == 'y':
Ax1 = plotAnomalyXYplane(Ax1,f,X,Y,ztx,Hy,Comp,Phase)
Ax1 = plotPlaceTxRxSphereXY(Ax1,xtx,ytx,xrx,yrx,x0,y0,a)
Ax2 = plotResponseFEM(Ax2,f,fvec,Hyi,Comp)
elif Comp == 'z':
Ax1 = plotAnomalyXYplane(Ax1,f,X,Y,ztx,Hz,Comp,Phase)
Ax1 = plotPlaceTxRxSphereXY(Ax1,xtx,ytx,xrx,yrx,x0,y0,a)
Ax2 = plotResponseFEM(Ax2,f,fvec,Hzi,Comp)
elif Comp == 'abs':
Ax1 = plotAnomalyXYplane(Ax1,f,X,Y,ztx,Habs,Comp,Phase)
Ax1 = plotPlaceTxRxSphereXY(Ax1,xtx,ytx,xrx,yrx,x0,y0,a)
Ax2 = plotResponseFEM(Ax2,f,fvec,Habsi,Comp)
plt.show(fig1)
def calc_IndCurrent_TD_offtime(self):
"""Gives FD induced current spectrum"""
#INITIALIZE ATTRIBUTES
Bpx = self.Bpx
Bpz = self.Bpz
a2 = self.a2
azm = np.pi*self.azm/180.
R = self.R
L = self.L
t = np.logspace(-6,0,101)
Ax = np.pi*a2**2*np.sin(azm)
Az = np.pi*a2**2*np.cos(azm)
Phi = (Ax*Bpx + Az*Bpz)
Is = (Phi/L)*np.exp(-(R/L)*t)
V = (Phi*R/L)*np.exp(-(R/L)*t) - (Phi*R/L**2)*np.exp(-(R/L)*t)
return V,Is
###########################################
# PLOTTING FUNCTIONS
###########################################
def plot_InducedCurrent_FD(self,Ax,Is,fi):
FS = 20
R = self.R
L = self.L
Imax = np.max(-np.real(Is))
f = np.logspace(0,8,101)
Ax.grid('both', linestyle='-', linewidth=0.8, color=[0.8, 0.8, 0.8])
Ax.semilogx(f,-np.real(Is),color='k',linewidth=4,label="$I_{Re}$")
Ax.semilogx(f,-np.imag(Is),color='k',ls='--',linewidth=4,label="$I_{Im}$")
Ax.semilogx(fi*np.array([1.,1.]),np.array([0,1.1*Imax]),color='r',ls='-',linewidth=3)
handles, labels = Ax.get_legend_handles_labels()
Ax.legend(handles, labels, loc='upper left', fontsize=FS)
Ax.set_xlabel('Frequency [Hz]',fontsize=FS+2)
Ax.set_ylabel('$\mathbf{- \, I_s (\omega)}$ [A]',fontsize=FS+2,labelpad=-10)
Ax.set_title('Frequency Response',fontsize=FS)
Ax.set_ybound(0,1.1*Imax)
Ax.tick_params(labelsize=FS-2)
Ax.yaxis.set_major_formatter(FormatStrFormatter('%.1e'))
#R_str = '{:.3e}'.format(R)
#L_str = '{:.3e}'.format(L)
#f_str = '{:.3e}'.format(fi)
#EMF_str = '{:.2e}j'.format(EMFi.imag)
#I_str = '{:.2e} - {:.2e}j'.format(float(np.real(Isi)),np.abs(float(np.imag(Isi))))
#Ax.text(1.4,1.01*Imax,'$R$ = '+R_str+' $\Omega$',fontsize=FS)
#Ax.text(1.4,0.94*Imax,'$L$ = '+L_str+' H',fontsize=FS)
#Ax.text(1.4,0.87*Imax,'$f$ = '+f_str+' Hz',fontsize=FS,color='r')
#Ax.text(1.4,0.8*Imax,'$V$ = '+EMF_str+' V',fontsize=FS,color='r')
#Ax.text(1.4,0.73*Imax,'$I_s$ = '+I_str+' A',fontsize=FS,color='r')
return Ax
def plot_InducedCurrent_TD(self,Ax,Is,ti,Vi,Isi):
FS = 20
R = self.R
L = self.L
Imax = np.max(Is)
t = np.logspace(-6,0,101)
Ax.grid('both', linestyle='-', linewidth=0.8, color=[0.8, 0.8, 0.8])
Ax.semilogx(t,Is,color='k',linewidth=4)
Ax.semilogx(ti*np.array([1.,1.]),np.array([0,1.3*Imax]),color='r',ls='-',linewidth=3)
Ax.set_xlabel('Time [s]',fontsize=FS+2)
Ax.set_ylabel('$\mathbf{I_s (\omega)}$ [A]',fontsize=FS+2,labelpad=-10)
Ax.set_title('Transient Induced Current',fontsize=FS)
Ax.set_ybound(0,1.2*Imax)
Ax.tick_params(labelsize=FS-2)
Ax.yaxis.set_major_formatter(FormatStrFormatter('%.1e'))
#R_str = '{:.3e}'.format(R)
#L_str = '{:.3e}'.format(L)
#t_str = '{:.3e}'.format(ti)
#V_str = '{:.3e}'.format(Vi)
#I_str = '{:.3e}'.format(Isi)
#Ax.text(1.4e-6,1.12*Imax,'$R$ = '+R_str+' $\Omega$',fontsize=FS)
#Ax.text(1.4e-6,1.04*Imax,'$L$ = '+L_str+' H',fontsize=FS)
#Ax.text(4e-2,1.12*Imax,'$t$ = '+t_str+' s',fontsize=FS,color='r')
#Ax.text(4e-2,1.04*Imax,'$V$ = '+V_str+' V',fontsize=FS,color='r')
#Ax.text(4e-2,0.96*Imax,'$I_s$ = '+I_str+' A',fontsize=FS,color='r')
return Ax
def __init__(self):
self.genMesh()
self.getCoreDomain()
# url = "http://em.geosci.xyz/_images/disc_dipole.png"
# response = requests.get(url)
# self.im = Image.open(StringIO(response.content))
self.time = np.logspace(-5, -2, 41)
def test_from_float_hex(self):
# IEEE doubles and floats only, otherwise the float32
# conversion may fail.
tgt = np.logspace(-10, 10, 5).astype(np.float32)
tgt = np.hstack((tgt, -tgt)).astype(np.float)
inp = '\n'.join(map(float.hex, tgt))
c = TextIO()
c.write(inp)
for dt in [np.float, np.float32]:
c.seek(0)
res = np.loadtxt(c, dtype=dt)
assert_equal(res, tgt, err_msg="%s" % dt)
def FrequencyAxis(Fmin, Fmax, Fnum, Log=True):
"""
Compute a lin/log spaced frequency axis.
"""
# Computing frequency axis
if Log:
Freq = _np.logspace(_np.log10(Fmin), _np.log10(Fmax), Fnum)
else:
Freq = _np.linspace(Fmin, Fmax, Fnum)
return Freq
def test_call(self):
"""Test that the bump EOS can be called
"""
bump_eos = EOSBump()
vol = np.logspace(np.log10(.1), np.log10(1), 50)
pressure = bump_eos(vol)
self.assertIsInstance(pressure, np.ndarray)
self.assertEqual(50, len(pressure))
def test_derivative(self):
"""Test the derivative function
"""
bump_eos = EOSBump()
vol = np.logspace(np.log10(.1), np.log10(1), 50)
pressure = bump_eos.derivative()(vol)
self.assertIsInstance(pressure, np.ndarray)
self.assertEqual(50, len(pressure))
def test_bad_derivative(self):
"""Tests that derivative errors are caught
"""
bump_eos = EOSBump()
vol = np.logspace(np.log10(.1), np.log10(1), 50)
with self.assertRaises(IOError):
pressure = bump_eos.derivative(order=2)(vol)
# end
p_fun = lambda v: 2.56e9 / v**3
