In [1]:
%matplotlib inline
import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import sklearn.linear_model as sl
import scipy.stats as st
%load_ext cython
ncaa = pd.read_csv("http://www4.stat.ncsu.edu/~boos/var.select/ncaa.data2.txt",
delim_whitespace = True)
x = ncaa.ix[:,:-1]
y = ncaa.ix[:,-1]
x.head()
Out[1]:
x1
x2
x3
x4
x5
x6
x7
x8
x9
x10
x11
x12
x13
x14
x15
x16
x17
x18
x19
0
13
17
9
15
28.0
0
-1.14045
3.660
4.490
3409
65.8
18
81
42.2
660000
77
100
59
1
1
28
20
32
18
18.4
18
-0.13719
2.594
3.610
7258
66.3
17
82
40.5
150555
88
94
41
25
2
32
20
20
20
34.8
18
1.55358
2.060
4.930
6405
75.0
19
71
46.5
415400
94
81
25
36
3
32
21
24
21
14.5
20
2.05712
2.887
3.876
18294
66.0
16
84
42.2
211000
93
88
26
13
4
24
20
16
20
21.8
13
-0.77082
2.565
4.960
8259
63.5
16
91
41.2
44000
90
92
32
31
In [2]:
pval = np.array([ 1.11022302e-16, 6.95109478e-05, 1.15845889e-02,
5.29907450e-03, 2.48342482e-03, 4.33151925e-02,
5.27341418e-02, 1.05631519e-01, 8.26270613e-02,
5.36356789e-02, 2.34958569e-01, 2.86440303e-01,
3.16318194e-01, 2.69726331e-01, 4.95325787e-01,
6.32648734e-01, 7.05641795e-01, 8.60540370e-01,
9.03197593e-01])
In [3]:
%%cython -a
import cython
import numpy as np
@cython.boundscheck(False)
@cython.wraparound(False)
def helper(double[:] pv_orig, int m1):
cdef int i,j
cdef double[:] pvm = np.zeros(m1)
cdef double[:] alpha = np.zeros(m1+1)
cdef double[:] alpha2 = np.zeros(m1)
cdef double[:] ghat2 = np.zeros(m1)
cdef double[:] S = np.zeros(m1+1)
cdef double[:] ghat = np.zeros(m1+1)
alpha[0] = 0
for i in range(0, m1):
pvm[i] = max(pv_orig[0:(i+1)])
# calculate model size
with cython.nogil:
for i in range(1, m1+1):
alpha[i] = pvm[i-1]
alpha2[i-1] = pvm[i-1] - 0.0000001
return (np.array(pvm), np.array(alpha), np.array(alpha2))
Out[3]:
Cython: _cython_magic_748f6478bc09a921a3b2c45d450b6d07.pyx
Generated by Cython 0.24
Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.
+01: import cython
__pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 1, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_1) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+02: import numpy as np
__pyx_t_1 = __Pyx_Import(__pyx_n_s_numpy, 0, -1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
03:
04: @cython.boundscheck(False)
05: @cython.wraparound(False)
+06: def helper(double[:] pv_orig, int m1):
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_748f6478bc09a921a3b2c45d450b6d07_1helper(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_748f6478bc09a921a3b2c45d450b6d07_1helper = {"helper", (PyCFunction)__pyx_pw_46_cython_magic_748f6478bc09a921a3b2c45d450b6d07_1helper, METH_VARARGS|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_46_cython_magic_748f6478bc09a921a3b2c45d450b6d07_1helper(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
__Pyx_memviewslice __pyx_v_pv_orig = { 0, 0, { 0 }, { 0 }, { 0 } };
int __pyx_v_m1;
PyObject *__pyx_r = 0;
__Pyx_RefNannyDeclarations
__Pyx_RefNannySetupContext("helper (wrapper)", 0);
{
static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_pv_orig,&__pyx_n_s_m1,0};
PyObject* values[2] = {0,0};
if (unlikely(__pyx_kwds)) {
Py_ssize_t kw_args;
const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
switch (pos_args) {
case 2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
case 1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
case 0: break;
default: goto __pyx_L5_argtuple_error;
}
kw_args = PyDict_Size(__pyx_kwds);
switch (pos_args) {
case 0:
if (likely((values[0] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_pv_orig)) != 0)) kw_args--;
else goto __pyx_L5_argtuple_error;
case 1:
if (likely((values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_m1)) != 0)) kw_args--;
else {
__Pyx_RaiseArgtupleInvalid("helper", 1, 2, 2, 1); __PYX_ERR(0, 6, __pyx_L3_error)
}
}
if (unlikely(kw_args > 0)) {
if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "helper") < 0)) __PYX_ERR(0, 6, __pyx_L3_error)
}
} else if (PyTuple_GET_SIZE(__pyx_args) != 2) {
goto __pyx_L5_argtuple_error;
} else {
values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
}
__pyx_v_pv_orig = __Pyx_PyObject_to_MemoryviewSlice_ds_double(values[0]); if (unlikely(!__pyx_v_pv_orig.memview)) __PYX_ERR(0, 6, __pyx_L3_error)
__pyx_v_m1 = __Pyx_PyInt_As_int(values[1]); if (unlikely((__pyx_v_m1 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 6, __pyx_L3_error)
}
goto __pyx_L4_argument_unpacking_done;
__pyx_L5_argtuple_error:;
__Pyx_RaiseArgtupleInvalid("helper", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 6, __pyx_L3_error)
__pyx_L3_error:;
__Pyx_AddTraceback("_cython_magic_748f6478bc09a921a3b2c45d450b6d07.helper", __pyx_clineno, __pyx_lineno, __pyx_filename);
__Pyx_RefNannyFinishContext();
return NULL;
__pyx_L4_argument_unpacking_done:;
__pyx_r = __pyx_pf_46_cython_magic_748f6478bc09a921a3b2c45d450b6d07_helper(__pyx_self, __pyx_v_pv_orig, __pyx_v_m1);
/* function exit code */
__Pyx_RefNannyFinishContext();
return __pyx_r;
}
static PyObject *__pyx_pf_46_cython_magic_748f6478bc09a921a3b2c45d450b6d07_helper(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_pv_orig, int __pyx_v_m1) {
int __pyx_v_i;
__Pyx_memviewslice __pyx_v_pvm = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_alpha = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_alpha2 = { 0, 0, { 0 }, { 0 }, { 0 } };
CYTHON_UNUSED __Pyx_memviewslice __pyx_v_ghat2 = { 0, 0, { 0 }, { 0 }, { 0 } };
CYTHON_UNUSED __Pyx_memviewslice __pyx_v_S = { 0, 0, { 0 }, { 0 }, { 0 } };
CYTHON_UNUSED __Pyx_memviewslice __pyx_v_ghat = { 0, 0, { 0 }, { 0 }, { 0 } };
PyObject *__pyx_r = NULL;
__Pyx_RefNannyDeclarations
__Pyx_RefNannySetupContext("helper", 0);
/* … */
/* function exit code */
__pyx_L1_error:;
__Pyx_XDECREF(__pyx_t_1);
__Pyx_XDECREF(__pyx_t_2);
__Pyx_XDECREF(__pyx_t_3);
__Pyx_XDECREF(__pyx_t_4);
__Pyx_XDECREF(__pyx_t_5);
__PYX_XDEC_MEMVIEW(&__pyx_t_6, 1);
__Pyx_XDECREF(__pyx_t_18);
__Pyx_XDECREF(__pyx_t_19);
__Pyx_AddTraceback("_cython_magic_748f6478bc09a921a3b2c45d450b6d07.helper", __pyx_clineno, __pyx_lineno, __pyx_filename);
__pyx_r = NULL;
__pyx_L0:;
__PYX_XDEC_MEMVIEW(&__pyx_v_pvm, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_alpha, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_alpha2, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_ghat2, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_S, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_ghat, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_pv_orig, 1);
__Pyx_XGIVEREF(__pyx_r);
__Pyx_RefNannyFinishContext();
return __pyx_r;
}
/* … */
__pyx_tuple__14 = PyTuple_Pack(10, __pyx_n_s_pv_orig, __pyx_n_s_m1, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_pvm, __pyx_n_s_alpha, __pyx_n_s_alpha2, __pyx_n_s_ghat2, __pyx_n_s_S, __pyx_n_s_ghat); if (unlikely(!__pyx_tuple__14)) __PYX_ERR(0, 6, __pyx_L1_error)
__Pyx_GOTREF(__pyx_tuple__14);
__Pyx_GIVEREF(__pyx_tuple__14);
/* … */
__pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_748f6478bc09a921a3b2c45d450b6d07_1helper, NULL, __pyx_n_s_cython_magic_748f6478bc09a921a3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 6, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_helper, __pyx_t_1) < 0) __PYX_ERR(0, 6, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_codeobj__15 = (PyObject*)__Pyx_PyCode_New(2, 0, 10, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__14, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_Rita_ipython_cython__cyth, __pyx_n_s_helper, 6, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__15)) __PYX_ERR(0, 6, __pyx_L1_error)
07: cdef int i,j
+08: cdef double[:] pvm = np.zeros(m1)
__pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 8, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_zeros); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 8, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_t_2 = __Pyx_PyInt_From_int(__pyx_v_m1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 8, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_4 = NULL;
if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_3))) {
__pyx_t_4 = PyMethod_GET_SELF(__pyx_t_3);
if (likely(__pyx_t_4)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
__Pyx_INCREF(__pyx_t_4);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_3, function);
}
}
if (!__pyx_t_4) {
__pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 8, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else {
__pyx_t_5 = PyTuple_New(1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 8, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_GIVEREF(__pyx_t_4); PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_4); __pyx_t_4 = NULL;
__Pyx_GIVEREF(__pyx_t_2);
PyTuple_SET_ITEM(__pyx_t_5, 0+1, __pyx_t_2);
__pyx_t_2 = 0;
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_5, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 8, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
}
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_ds_double(__pyx_t_1);
if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 8, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_v_pvm = __pyx_t_6;
__pyx_t_6.memview = NULL;
__pyx_t_6.data = NULL;
+09: cdef double[:] alpha = np.zeros(m1+1)
__pyx_t_3 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 9, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_zeros); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 9, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_3 = __Pyx_PyInt_From_long((__pyx_v_m1 + 1)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 9, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_2 = NULL;
if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_5))) {
__pyx_t_2 = PyMethod_GET_SELF(__pyx_t_5);
if (likely(__pyx_t_2)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_5);
__Pyx_INCREF(__pyx_t_2);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_5, function);
}
}
if (!__pyx_t_2) {
__pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_5, __pyx_t_3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 9, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else {
__pyx_t_4 = PyTuple_New(1+1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 9, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_GIVEREF(__pyx_t_2); PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_2); __pyx_t_2 = NULL;
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_4, 0+1, __pyx_t_3);
__pyx_t_3 = 0;
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_5, __pyx_t_4, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 9, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
}
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_ds_double(__pyx_t_1);
if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 9, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_v_alpha = __pyx_t_6;
__pyx_t_6.memview = NULL;
__pyx_t_6.data = NULL;
+10: cdef double[:] alpha2 = np.zeros(m1)
__pyx_t_5 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 10, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_n_s_zeros); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 10, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__pyx_t_5 = __Pyx_PyInt_From_int(__pyx_v_m1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 10, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__pyx_t_3 = NULL;
if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_4))) {
__pyx_t_3 = PyMethod_GET_SELF(__pyx_t_4);
if (likely(__pyx_t_3)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);
__Pyx_INCREF(__pyx_t_3);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_4, function);
}
}
if (!__pyx_t_3) {
__pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_5); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 10, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else {
__pyx_t_2 = PyTuple_New(1+1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 10, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_GIVEREF(__pyx_t_3); PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_t_3); __pyx_t_3 = NULL;
__Pyx_GIVEREF(__pyx_t_5);
PyTuple_SET_ITEM(__pyx_t_2, 0+1, __pyx_t_5);
__pyx_t_5 = 0;
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_2, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 10, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
}
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_ds_double(__pyx_t_1);
if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 10, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_v_alpha2 = __pyx_t_6;
__pyx_t_6.memview = NULL;
__pyx_t_6.data = NULL;
+11: cdef double[:] ghat2 = np.zeros(m1)
__pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 11, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_zeros); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 11, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__pyx_t_4 = __Pyx_PyInt_From_int(__pyx_v_m1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 11, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__pyx_t_5 = NULL;
if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_2))) {
__pyx_t_5 = PyMethod_GET_SELF(__pyx_t_2);
if (likely(__pyx_t_5)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_2);
__Pyx_INCREF(__pyx_t_5);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_2, function);
}
}
if (!__pyx_t_5) {
__pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_2, __pyx_t_4); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 11, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else {
__pyx_t_3 = PyTuple_New(1+1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 11, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_GIVEREF(__pyx_t_5); PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_5); __pyx_t_5 = NULL;
__Pyx_GIVEREF(__pyx_t_4);
PyTuple_SET_ITEM(__pyx_t_3, 0+1, __pyx_t_4);
__pyx_t_4 = 0;
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_3, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 11, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
}
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_ds_double(__pyx_t_1);
if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 11, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_v_ghat2 = __pyx_t_6;
__pyx_t_6.memview = NULL;
__pyx_t_6.data = NULL;
+12: cdef double[:] S = np.zeros(m1+1)
__pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 12, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_zeros); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 12, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_t_2 = __Pyx_PyInt_From_long((__pyx_v_m1 + 1)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 12, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_4 = NULL;
if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_3))) {
__pyx_t_4 = PyMethod_GET_SELF(__pyx_t_3);
if (likely(__pyx_t_4)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
__Pyx_INCREF(__pyx_t_4);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_3, function);
}
}
if (!__pyx_t_4) {
__pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 12, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else {
__pyx_t_5 = PyTuple_New(1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 12, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_GIVEREF(__pyx_t_4); PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_4); __pyx_t_4 = NULL;
__Pyx_GIVEREF(__pyx_t_2);
PyTuple_SET_ITEM(__pyx_t_5, 0+1, __pyx_t_2);
__pyx_t_2 = 0;
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_5, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 12, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
}
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_ds_double(__pyx_t_1);
if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 12, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_v_S = __pyx_t_6;
__pyx_t_6.memview = NULL;
__pyx_t_6.data = NULL;
+13: cdef double[:] ghat = np.zeros(m1+1)
__pyx_t_3 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_zeros); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_3 = __Pyx_PyInt_From_long((__pyx_v_m1 + 1)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_2 = NULL;
if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_5))) {
__pyx_t_2 = PyMethod_GET_SELF(__pyx_t_5);
if (likely(__pyx_t_2)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_5);
__Pyx_INCREF(__pyx_t_2);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_5, function);
}
}
if (!__pyx_t_2) {
__pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_5, __pyx_t_3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else {
__pyx_t_4 = PyTuple_New(1+1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_GIVEREF(__pyx_t_2); PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_2); __pyx_t_2 = NULL;
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_4, 0+1, __pyx_t_3);
__pyx_t_3 = 0;
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_5, __pyx_t_4, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
}
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_ds_double(__pyx_t_1);
if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_v_ghat = __pyx_t_6;
__pyx_t_6.memview = NULL;
__pyx_t_6.data = NULL;
14:
+15: alpha[0] = 0
__pyx_t_7 = 0;
*((double *) ( /* dim=0 */ (__pyx_v_alpha.data + __pyx_t_7 * __pyx_v_alpha.strides[0]) )) = 0.0;
+16: for i in range(0, m1):
__pyx_t_8 = __pyx_v_m1;
for (__pyx_t_9 = 0; __pyx_t_9 < __pyx_t_8; __pyx_t_9+=1) {
__pyx_v_i = __pyx_t_9;
+17: pvm[i] = max(pv_orig[0:(i+1)])
__pyx_t_6.data = __pyx_v_pv_orig.data;
__pyx_t_6.memview = __pyx_v_pv_orig.memview;
__PYX_INC_MEMVIEW(&__pyx_t_6, 0);
__pyx_t_10 = -1;
if (unlikely(__pyx_memoryview_slice_memviewslice(
&__pyx_t_6,
__pyx_v_pv_orig.shape[0], __pyx_v_pv_orig.strides[0], __pyx_v_pv_orig.suboffsets[0],
0,
0,
&__pyx_t_10,
0,
(__pyx_v_i + 1),
0,
1,
1,
0,
1) < 0))
{
__PYX_ERR(0, 17, __pyx_L1_error)
}
__pyx_t_1 = __pyx_memoryview_fromslice(__pyx_t_6, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__PYX_XDEC_MEMVIEW(&__pyx_t_6, 1);
__pyx_t_5 = PyTuple_New(1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_GIVEREF(__pyx_t_1);
PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_1);
__pyx_t_1 = 0;
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_builtin_max, __pyx_t_5, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__pyx_t_11 = __pyx_PyFloat_AsDouble(__pyx_t_1); if (unlikely((__pyx_t_11 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_t_12 = __pyx_v_i;
*((double *) ( /* dim=0 */ (__pyx_v_pvm.data + __pyx_t_12 * __pyx_v_pvm.strides[0]) )) = __pyx_t_11;
}
18:
19: # calculate model size
+20: with cython.nogil:
{
#ifdef WITH_THREAD
PyThreadState *_save;
Py_UNBLOCK_THREADS
#endif
/*try:*/ {
/* … */
/*finally:*/ {
/*normal exit:*/{
#ifdef WITH_THREAD
Py_BLOCK_THREADS
#endif
goto __pyx_L7;
}
__pyx_L7:;
}
}
+21: for i in range(1, m1+1):
__pyx_t_13 = (__pyx_v_m1 + 1);
for (__pyx_t_8 = 1; __pyx_t_8 < __pyx_t_13; __pyx_t_8+=1) {
__pyx_v_i = __pyx_t_8;
+22: alpha[i] = pvm[i-1]
__pyx_t_14 = (__pyx_v_i - 1);
__pyx_t_15 = __pyx_v_i;
*((double *) ( /* dim=0 */ (__pyx_v_alpha.data + __pyx_t_15 * __pyx_v_alpha.strides[0]) )) = (*((double *) ( /* dim=0 */ (__pyx_v_pvm.data + __pyx_t_14 * __pyx_v_pvm.strides[0]) )));
+23: alpha2[i-1] = pvm[i-1] - 0.0000001
__pyx_t_16 = (__pyx_v_i - 1);
__pyx_t_17 = (__pyx_v_i - 1);
*((double *) ( /* dim=0 */ (__pyx_v_alpha2.data + __pyx_t_17 * __pyx_v_alpha2.strides[0]) )) = ((*((double *) ( /* dim=0 */ (__pyx_v_pvm.data + __pyx_t_16 * __pyx_v_pvm.strides[0]) ))) - 0.0000001);
}
}
+24: return (np.array(pvm), np.array(alpha), np.array(alpha2))
__Pyx_XDECREF(__pyx_r);
__pyx_t_5 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_n_s_array); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__pyx_t_5 = __pyx_memoryview_fromslice(__pyx_v_pvm, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__pyx_t_3 = NULL;
if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_4))) {
__pyx_t_3 = PyMethod_GET_SELF(__pyx_t_4);
if (likely(__pyx_t_3)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);
__Pyx_INCREF(__pyx_t_3);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_4, function);
}
}
if (!__pyx_t_3) {
__pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_5); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else {
__pyx_t_2 = PyTuple_New(1+1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_GIVEREF(__pyx_t_3); PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_t_3); __pyx_t_3 = NULL;
__Pyx_GIVEREF(__pyx_t_5);
PyTuple_SET_ITEM(__pyx_t_2, 0+1, __pyx_t_5);
__pyx_t_5 = 0;
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_2, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
}
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_array); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_t_2 = __pyx_memoryview_fromslice(__pyx_v_alpha, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_3 = NULL;
if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_5))) {
__pyx_t_3 = PyMethod_GET_SELF(__pyx_t_5);
if (likely(__pyx_t_3)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_5);
__Pyx_INCREF(__pyx_t_3);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_5, function);
}
}
if (!__pyx_t_3) {
__pyx_t_4 = __Pyx_PyObject_CallOneArg(__pyx_t_5, __pyx_t_2); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_4);
} else {
__pyx_t_18 = PyTuple_New(1+1); if (unlikely(!__pyx_t_18)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_18);
__Pyx_GIVEREF(__pyx_t_3); PyTuple_SET_ITEM(__pyx_t_18, 0, __pyx_t_3); __pyx_t_3 = NULL;
__Pyx_GIVEREF(__pyx_t_2);
PyTuple_SET_ITEM(__pyx_t_18, 0+1, __pyx_t_2);
__pyx_t_2 = 0;
__pyx_t_4 = __Pyx_PyObject_Call(__pyx_t_5, __pyx_t_18, NULL); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_18); __pyx_t_18 = 0;
}
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__pyx_t_18 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_18)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_18);
__pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_18, __pyx_n_s_array); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_18); __pyx_t_18 = 0;
__pyx_t_18 = __pyx_memoryview_fromslice(__pyx_v_alpha2, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_18)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_18);
__pyx_t_3 = NULL;
if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_2))) {
__pyx_t_3 = PyMethod_GET_SELF(__pyx_t_2);
if (likely(__pyx_t_3)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_2);
__Pyx_INCREF(__pyx_t_3);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_2, function);
}
}
if (!__pyx_t_3) {
__pyx_t_5 = __Pyx_PyObject_CallOneArg(__pyx_t_2, __pyx_t_18); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_18); __pyx_t_18 = 0;
__Pyx_GOTREF(__pyx_t_5);
} else {
__pyx_t_19 = PyTuple_New(1+1); if (unlikely(!__pyx_t_19)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_19);
__Pyx_GIVEREF(__pyx_t_3); PyTuple_SET_ITEM(__pyx_t_19, 0, __pyx_t_3); __pyx_t_3 = NULL;
__Pyx_GIVEREF(__pyx_t_18);
PyTuple_SET_ITEM(__pyx_t_19, 0+1, __pyx_t_18);
__pyx_t_18 = 0;
__pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_19, NULL); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_DECREF(__pyx_t_19); __pyx_t_19 = 0;
}
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_t_2 = PyTuple_New(3); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_GIVEREF(__pyx_t_1);
PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_t_1);
__Pyx_GIVEREF(__pyx_t_4);
PyTuple_SET_ITEM(__pyx_t_2, 1, __pyx_t_4);
__Pyx_GIVEREF(__pyx_t_5);
PyTuple_SET_ITEM(__pyx_t_2, 2, __pyx_t_5);
__pyx_t_1 = 0;
__pyx_t_4 = 0;
__pyx_t_5 = 0;
__pyx_r = __pyx_t_2;
__pyx_t_2 = 0;
goto __pyx_L0;
In [4]:
%%cython -a
import cython
import numpy as np
from cython.parallel import parallel, prange
@cython.boundscheck(False)
@cython.wraparound(False)
def helper2(double[:] S, double[:] alpha, double[:] alpha2, int m1, int m, int ng):
cdef int i
cdef double[:] ghat2 = np.zeros(m1)
cdef double[:] ghat = np.zeros(m1+1)
with cython.nogil:
for i in range(1, m1+1):
ghat[i] = (m - S[i])*alpha[i]/(1 + S[i])
ghat2[i-1] = (m-S[i-1])*alpha2[i-1]/(1+S[i-1])
return (np.array(ghat2), np.array(ghat))
Out[4]:
Cython: _cython_magic_7ba0232d9e0da78a0a34c3307e217fdd.pyx
Generated by Cython 0.24
Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.
+01: import cython
__pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 1, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_1) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+02: import numpy as np
__pyx_t_1 = __Pyx_Import(__pyx_n_s_numpy, 0, -1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
03: from cython.parallel import parallel, prange
04:
05: @cython.boundscheck(False)
06: @cython.wraparound(False)
+07: def helper2(double[:] S, double[:] alpha, double[:] alpha2, int m1, int m, int ng):
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_7ba0232d9e0da78a0a34c3307e217fdd_1helper2(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_7ba0232d9e0da78a0a34c3307e217fdd_1helper2 = {"helper2", (PyCFunction)__pyx_pw_46_cython_magic_7ba0232d9e0da78a0a34c3307e217fdd_1helper2, METH_VARARGS|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_46_cython_magic_7ba0232d9e0da78a0a34c3307e217fdd_1helper2(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
__Pyx_memviewslice __pyx_v_S = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_alpha = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_alpha2 = { 0, 0, { 0 }, { 0 }, { 0 } };
int __pyx_v_m1;
int __pyx_v_m;
CYTHON_UNUSED int __pyx_v_ng;
PyObject *__pyx_r = 0;
__Pyx_RefNannyDeclarations
__Pyx_RefNannySetupContext("helper2 (wrapper)", 0);
{
static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_S,&__pyx_n_s_alpha,&__pyx_n_s_alpha2,&__pyx_n_s_m1,&__pyx_n_s_m,&__pyx_n_s_ng,0};
PyObject* values[6] = {0,0,0,0,0,0};
if (unlikely(__pyx_kwds)) {
Py_ssize_t kw_args;
const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
switch (pos_args) {
case 6: values[5] = PyTuple_GET_ITEM(__pyx_args, 5);
case 5: values[4] = PyTuple_GET_ITEM(__pyx_args, 4);
case 4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
case 3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
case 2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
case 1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
case 0: break;
default: goto __pyx_L5_argtuple_error;
}
kw_args = PyDict_Size(__pyx_kwds);
switch (pos_args) {
case 0:
if (likely((values[0] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_S)) != 0)) kw_args--;
else goto __pyx_L5_argtuple_error;
case 1:
if (likely((values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_alpha)) != 0)) kw_args--;
else {
__Pyx_RaiseArgtupleInvalid("helper2", 1, 6, 6, 1); __PYX_ERR(0, 7, __pyx_L3_error)
}
case 2:
if (likely((values[2] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_alpha2)) != 0)) kw_args--;
else {
__Pyx_RaiseArgtupleInvalid("helper2", 1, 6, 6, 2); __PYX_ERR(0, 7, __pyx_L3_error)
}
case 3:
if (likely((values[3] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_m1)) != 0)) kw_args--;
else {
__Pyx_RaiseArgtupleInvalid("helper2", 1, 6, 6, 3); __PYX_ERR(0, 7, __pyx_L3_error)
}
case 4:
if (likely((values[4] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_m)) != 0)) kw_args--;
else {
__Pyx_RaiseArgtupleInvalid("helper2", 1, 6, 6, 4); __PYX_ERR(0, 7, __pyx_L3_error)
}
case 5:
if (likely((values[5] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_ng)) != 0)) kw_args--;
else {
__Pyx_RaiseArgtupleInvalid("helper2", 1, 6, 6, 5); __PYX_ERR(0, 7, __pyx_L3_error)
}
}
if (unlikely(kw_args > 0)) {
if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "helper2") < 0)) __PYX_ERR(0, 7, __pyx_L3_error)
}
} else if (PyTuple_GET_SIZE(__pyx_args) != 6) {
goto __pyx_L5_argtuple_error;
} else {
values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
values[4] = PyTuple_GET_ITEM(__pyx_args, 4);
values[5] = PyTuple_GET_ITEM(__pyx_args, 5);
}
__pyx_v_S = __Pyx_PyObject_to_MemoryviewSlice_ds_double(values[0]); if (unlikely(!__pyx_v_S.memview)) __PYX_ERR(0, 7, __pyx_L3_error)
__pyx_v_alpha = __Pyx_PyObject_to_MemoryviewSlice_ds_double(values[1]); if (unlikely(!__pyx_v_alpha.memview)) __PYX_ERR(0, 7, __pyx_L3_error)
__pyx_v_alpha2 = __Pyx_PyObject_to_MemoryviewSlice_ds_double(values[2]); if (unlikely(!__pyx_v_alpha2.memview)) __PYX_ERR(0, 7, __pyx_L3_error)
__pyx_v_m1 = __Pyx_PyInt_As_int(values[3]); if (unlikely((__pyx_v_m1 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 7, __pyx_L3_error)
__pyx_v_m = __Pyx_PyInt_As_int(values[4]); if (unlikely((__pyx_v_m == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 7, __pyx_L3_error)
__pyx_v_ng = __Pyx_PyInt_As_int(values[5]); if (unlikely((__pyx_v_ng == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 7, __pyx_L3_error)
}
goto __pyx_L4_argument_unpacking_done;
__pyx_L5_argtuple_error:;
__Pyx_RaiseArgtupleInvalid("helper2", 1, 6, 6, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 7, __pyx_L3_error)
__pyx_L3_error:;
__Pyx_AddTraceback("_cython_magic_7ba0232d9e0da78a0a34c3307e217fdd.helper2", __pyx_clineno, __pyx_lineno, __pyx_filename);
__Pyx_RefNannyFinishContext();
return NULL;
__pyx_L4_argument_unpacking_done:;
__pyx_r = __pyx_pf_46_cython_magic_7ba0232d9e0da78a0a34c3307e217fdd_helper2(__pyx_self, __pyx_v_S, __pyx_v_alpha, __pyx_v_alpha2, __pyx_v_m1, __pyx_v_m, __pyx_v_ng);
/* function exit code */
__Pyx_RefNannyFinishContext();
return __pyx_r;
}
static PyObject *__pyx_pf_46_cython_magic_7ba0232d9e0da78a0a34c3307e217fdd_helper2(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_S, __Pyx_memviewslice __pyx_v_alpha, __Pyx_memviewslice __pyx_v_alpha2, int __pyx_v_m1, int __pyx_v_m, CYTHON_UNUSED int __pyx_v_ng) {
int __pyx_v_i;
__Pyx_memviewslice __pyx_v_ghat2 = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_ghat = { 0, 0, { 0 }, { 0 }, { 0 } };
PyObject *__pyx_r = NULL;
__Pyx_RefNannyDeclarations
__Pyx_RefNannySetupContext("helper2", 0);
/* … */
/* function exit code */
__pyx_L1_error:;
__Pyx_XDECREF(__pyx_t_1);
__Pyx_XDECREF(__pyx_t_2);
__Pyx_XDECREF(__pyx_t_3);
__Pyx_XDECREF(__pyx_t_4);
__Pyx_XDECREF(__pyx_t_5);
__PYX_XDEC_MEMVIEW(&__pyx_t_6, 1);
__Pyx_XDECREF(__pyx_t_19);
__Pyx_AddTraceback("_cython_magic_7ba0232d9e0da78a0a34c3307e217fdd.helper2", __pyx_clineno, __pyx_lineno, __pyx_filename);
__pyx_r = NULL;
__pyx_L0:;
__PYX_XDEC_MEMVIEW(&__pyx_v_ghat2, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_ghat, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_S, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_alpha, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_alpha2, 1);
__Pyx_XGIVEREF(__pyx_r);
__Pyx_RefNannyFinishContext();
return __pyx_r;
}
/* … */
__pyx_tuple__14 = PyTuple_Pack(9, __pyx_n_s_S, __pyx_n_s_alpha, __pyx_n_s_alpha2, __pyx_n_s_m1, __pyx_n_s_m, __pyx_n_s_ng, __pyx_n_s_i, __pyx_n_s_ghat2, __pyx_n_s_ghat); if (unlikely(!__pyx_tuple__14)) __PYX_ERR(0, 7, __pyx_L1_error)
__Pyx_GOTREF(__pyx_tuple__14);
__Pyx_GIVEREF(__pyx_tuple__14);
/* … */
__pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_7ba0232d9e0da78a0a34c3307e217fdd_1helper2, NULL, __pyx_n_s_cython_magic_7ba0232d9e0da78a0a); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 7, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_helper2, __pyx_t_1) < 0) __PYX_ERR(0, 7, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_codeobj__15 = (PyObject*)__Pyx_PyCode_New(6, 0, 9, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__14, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_Rita_ipython_cython__cyth, __pyx_n_s_helper2, 7, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__15)) __PYX_ERR(0, 7, __pyx_L1_error)
08: cdef int i
+09: cdef double[:] ghat2 = np.zeros(m1)
__pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 9, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_zeros); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 9, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_t_2 = __Pyx_PyInt_From_int(__pyx_v_m1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 9, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_4 = NULL;
if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_3))) {
__pyx_t_4 = PyMethod_GET_SELF(__pyx_t_3);
if (likely(__pyx_t_4)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
__Pyx_INCREF(__pyx_t_4);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_3, function);
}
}
if (!__pyx_t_4) {
__pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 9, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else {
__pyx_t_5 = PyTuple_New(1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 9, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_GIVEREF(__pyx_t_4); PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_4); __pyx_t_4 = NULL;
__Pyx_GIVEREF(__pyx_t_2);
PyTuple_SET_ITEM(__pyx_t_5, 0+1, __pyx_t_2);
__pyx_t_2 = 0;
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_5, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 9, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
}
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_ds_double(__pyx_t_1);
if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 9, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_v_ghat2 = __pyx_t_6;
__pyx_t_6.memview = NULL;
__pyx_t_6.data = NULL;
+10: cdef double[:] ghat = np.zeros(m1+1)
__pyx_t_3 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 10, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_zeros); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 10, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_3 = __Pyx_PyInt_From_long((__pyx_v_m1 + 1)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 10, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_2 = NULL;
if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_5))) {
__pyx_t_2 = PyMethod_GET_SELF(__pyx_t_5);
if (likely(__pyx_t_2)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_5);
__Pyx_INCREF(__pyx_t_2);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_5, function);
}
}
if (!__pyx_t_2) {
__pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_5, __pyx_t_3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 10, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else {
__pyx_t_4 = PyTuple_New(1+1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 10, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_GIVEREF(__pyx_t_2); PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_2); __pyx_t_2 = NULL;
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_4, 0+1, __pyx_t_3);
__pyx_t_3 = 0;
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_5, __pyx_t_4, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 10, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
}
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_ds_double(__pyx_t_1);
if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 10, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_v_ghat = __pyx_t_6;
__pyx_t_6.memview = NULL;
__pyx_t_6.data = NULL;
11:
+12: with cython.nogil:
{
#ifdef WITH_THREAD
PyThreadState *_save;
Py_UNBLOCK_THREADS
#endif
/*try:*/ {
/* … */
/*finally:*/ {
/*normal exit:*/{
#ifdef WITH_THREAD
Py_BLOCK_THREADS
#endif
goto __pyx_L5;
}
__pyx_L4_error: {
#ifdef WITH_THREAD
Py_BLOCK_THREADS
#endif
goto __pyx_L1_error;
}
__pyx_L5:;
}
}
+13: for i in range(1, m1+1):
__pyx_t_7 = (__pyx_v_m1 + 1);
for (__pyx_t_8 = 1; __pyx_t_8 < __pyx_t_7; __pyx_t_8+=1) {
__pyx_v_i = __pyx_t_8;
+14: ghat[i] = (m - S[i])*alpha[i]/(1 + S[i])
__pyx_t_9 = __pyx_v_i;
__pyx_t_10 = __pyx_v_i;
__pyx_t_11 = ((__pyx_v_m - (*((double *) ( /* dim=0 */ (__pyx_v_S.data + __pyx_t_9 * __pyx_v_S.strides[0]) )))) * (*((double *) ( /* dim=0 */ (__pyx_v_alpha.data + __pyx_t_10 * __pyx_v_alpha.strides[0]) ))));
__pyx_t_12 = __pyx_v_i;
__pyx_t_13 = (1.0 + (*((double *) ( /* dim=0 */ (__pyx_v_S.data + __pyx_t_12 * __pyx_v_S.strides[0]) ))));
if (unlikely(__pyx_t_13 == 0)) {
#ifdef WITH_THREAD
PyGILState_STATE __pyx_gilstate_save = PyGILState_Ensure();
#endif
PyErr_SetString(PyExc_ZeroDivisionError, "float division");
#ifdef WITH_THREAD
PyGILState_Release(__pyx_gilstate_save);
#endif
__PYX_ERR(0, 14, __pyx_L4_error)
}
__pyx_t_14 = __pyx_v_i;
*((double *) ( /* dim=0 */ (__pyx_v_ghat.data + __pyx_t_14 * __pyx_v_ghat.strides[0]) )) = (__pyx_t_11 / __pyx_t_13);
+15: ghat2[i-1] = (m-S[i-1])*alpha2[i-1]/(1+S[i-1])
__pyx_t_15 = (__pyx_v_i - 1);
__pyx_t_16 = (__pyx_v_i - 1);
__pyx_t_13 = ((__pyx_v_m - (*((double *) ( /* dim=0 */ (__pyx_v_S.data + __pyx_t_15 * __pyx_v_S.strides[0]) )))) * (*((double *) ( /* dim=0 */ (__pyx_v_alpha2.data + __pyx_t_16 * __pyx_v_alpha2.strides[0]) ))));
__pyx_t_17 = (__pyx_v_i - 1);
__pyx_t_11 = (1.0 + (*((double *) ( /* dim=0 */ (__pyx_v_S.data + __pyx_t_17 * __pyx_v_S.strides[0]) ))));
if (unlikely(__pyx_t_11 == 0)) {
#ifdef WITH_THREAD
PyGILState_STATE __pyx_gilstate_save = PyGILState_Ensure();
#endif
PyErr_SetString(PyExc_ZeroDivisionError, "float division");
#ifdef WITH_THREAD
PyGILState_Release(__pyx_gilstate_save);
#endif
__PYX_ERR(0, 15, __pyx_L4_error)
}
__pyx_t_18 = (__pyx_v_i - 1);
*((double *) ( /* dim=0 */ (__pyx_v_ghat2.data + __pyx_t_18 * __pyx_v_ghat2.strides[0]) )) = (__pyx_t_13 / __pyx_t_11);
}
}
16:
+17: return (np.array(ghat2), np.array(ghat))
__Pyx_XDECREF(__pyx_r);
__pyx_t_5 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_n_s_array); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__pyx_t_5 = __pyx_memoryview_fromslice(__pyx_v_ghat2, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__pyx_t_3 = NULL;
if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_4))) {
__pyx_t_3 = PyMethod_GET_SELF(__pyx_t_4);
if (likely(__pyx_t_3)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);
__Pyx_INCREF(__pyx_t_3);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_4, function);
}
}
if (!__pyx_t_3) {
__pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_5); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else {
__pyx_t_2 = PyTuple_New(1+1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_GIVEREF(__pyx_t_3); PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_t_3); __pyx_t_3 = NULL;
__Pyx_GIVEREF(__pyx_t_5);
PyTuple_SET_ITEM(__pyx_t_2, 0+1, __pyx_t_5);
__pyx_t_5 = 0;
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_2, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
}
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_array); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_t_2 = __pyx_memoryview_fromslice(__pyx_v_ghat, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_3 = NULL;
if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_5))) {
__pyx_t_3 = PyMethod_GET_SELF(__pyx_t_5);
if (likely(__pyx_t_3)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_5);
__Pyx_INCREF(__pyx_t_3);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_5, function);
}
}
if (!__pyx_t_3) {
__pyx_t_4 = __Pyx_PyObject_CallOneArg(__pyx_t_5, __pyx_t_2); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_4);
} else {
__pyx_t_19 = PyTuple_New(1+1); if (unlikely(!__pyx_t_19)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_19);
__Pyx_GIVEREF(__pyx_t_3); PyTuple_SET_ITEM(__pyx_t_19, 0, __pyx_t_3); __pyx_t_3 = NULL;
__Pyx_GIVEREF(__pyx_t_2);
PyTuple_SET_ITEM(__pyx_t_19, 0+1, __pyx_t_2);
__pyx_t_2 = 0;
__pyx_t_4 = __Pyx_PyObject_Call(__pyx_t_5, __pyx_t_19, NULL); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_19); __pyx_t_19 = 0;
}
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__pyx_t_5 = PyTuple_New(2); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_GIVEREF(__pyx_t_1);
PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_1);
__Pyx_GIVEREF(__pyx_t_4);
PyTuple_SET_ITEM(__pyx_t_5, 1, __pyx_t_4);
__pyx_t_1 = 0;
__pyx_t_4 = 0;
__pyx_r = __pyx_t_5;
__pyx_t_5 = 0;
goto __pyx_L0;
In [5]:
def reg_subset(x, y):
lm = sl.LinearRegression()
(n,m) = x.shape
in_ = []
out_ = list(range(m))
rss = np.zeros(m+1)
rss[0] = sum((y-np.mean(y))**2)
if(m>=n): ml = n-5
else: ml = m
for pi in range(m):
rss_find = []
for i in out_:
fit_X = pd.DataFrame(x.ix[:, in_ + [i]])
lm.fit(fit_X, y)
pred = lm.predict(fit_X)
rss_find.append(sum((pred-y)**2))
min_idx = np.argmin(rss_find)
min_var = out_[min_idx]
rss[pi+1] = np.min(rss_find)
in_.append(min_var)
del out_[min_idx]
in_ = np.array(in_)
in_var = x.columns[in_]
vm = np.array(range(ml))
pv_org = 1 - st.f.cdf((rss[vm] - rss[vm+1])*(n-(vm+2))/rss[vm+1],
1,n-(vm+2))
return (in_,np.array(in_var),rss, pv_org)
In [6]:
def fsr_fast_pv(pv_orig, m, gam0 = 0.05, digits = 4, printout = False, plot = False):
m1 = len(pv_orig)
pvm = np.zeros(m1)
# create monotone p-values
for i in range(0, m1):
pvm[i] = np.max(pv_orig[0:(i+1)])
alpha = np.concatenate([np.array([0]), pvm])
ng = len(alpha)
# calculate model size
S = np.zeros(ng)
for j in range(1, ng):
S[j] = sum(pvm <= alpha[j])
# calculate gamma hat
ghat = (m - S)*alpha/(1 + S)
# add additional points to make jumps
alpha2 = alpha[1:ng] - 0.0000001
ghat2 = (m - S[0:(ng - 1)])*alpha2/(1 + S[0:(ng - 1)])
zp = pd.DataFrame({'a': np.concatenate([alpha, alpha2]), 'g': np.concatenate([ghat, ghat2])})
zp.sort_values(by =['a', 'g'], ascending = [True, False], inplace = True)
# largest gamma hat and index
gamma_max = np.argmax(zp['g'])
alpha_max = zp['a'][gamma_max]
# model size with ghat just below gam0
ind = np.logical_and(ghat <= gam0, alpha <= alpha_max)*1
Sind = S[np.max(np.where(ind > 0))]
# calculate alpha_F
alpha_fast = (1 + Sind)*gam0/(m - Sind)
# size of model no intercept
size1 = sum(pvm <= alpha_fast)
# generate plot
if plot==True:
plt.plot(zp['a'], zp['g'], marker = 'o', markersize = 6)
plt.ylabel('Estimated Gamma')
plt.xlabel('Alpha')
pass
df1 = pd.DataFrame({'pval': pv_orig, 'pvmax': pvm, 'ghigh': ghat2, 'glow': ghat[1:ng]}, columns = ['pval', 'pvmax', 'ghigh', 'glow'])
df2 = pd.DataFrame({'m1': m1, 'm': m, 'gam0': gam0, 'size': size1, 'alphamax': alpha_max, 'alpha_fast': alpha_fast}, columns = ['m1', 'm', 'gam0', 'size', 'alphamax', 'alpha_fast'], index=[0])
if printout == True:
print(df1,df2)
return(np.round(df1, digits), np.round(df2, digits), ghat)
In [9]:
def fsr_fast_pv_cython(pv_orig, m, gam0 = 0.05, digits = 4, printout = False, plot = False):
m1 = len(pv_orig)
ng = m1+1
(pvm,alpha,alpha2) = helper(pv_orig, m1)
S = np.zeros(ng)
for j in range(1, ng):
S[j] = sum(pvm <= alpha[j])
# calculate gamma hat
(ghat2,ghat) = helper2(S,alpha,alpha2,m1,m,ng)
#ghat = (m - S)*alpha/(1 + S)
#ghat2 = (m - S[0:(ng - 1)])*alpha2/(1 + S[0:(ng - 1)])
zp = pd.DataFrame({'a': np.concatenate([alpha, alpha2]), 'g': np.concatenate([ghat, ghat2])})
zp.sort_values(by =['a', 'g'], ascending = [True, False], inplace = True)
# largest gamma hat and index
gamma_max = np.argmax(zp['g'])
alpha_max = zp['a'][gamma_max]
# model size with ghat just below gam0
ind = np.logical_and(ghat <= gam0, alpha <= alpha_max)*1
Sind = S[np.max(np.where(ind > 0))]
# calculate alpha_F
alpha_fast = (1 + Sind)*gam0/(m - Sind)
# size of model no intercept
size1 = sum(pvm <= alpha_fast)
# generate plot
if plot==True:
plt.plot(zp['a'], zp['g'], marker = 'o', markersize = 6)
plt.ylabel('Estimated Gamma')
plt.xlabel('Alpha')
pass
df1 = pd.DataFrame({'pval': pv_orig, 'pvmax': pvm, 'ghigh': ghat2, 'glow': ghat[1:ng]}, columns = ['pval', 'pvmax', 'ghigh', 'glow'])
df2 = pd.DataFrame({'m1': m1, 'm': m, 'gam0': gam0, 'size': size1, 'alphamax': alpha_max, 'alpha_fast': alpha_fast}, columns = ['m1', 'm', 'gam0', 'size', 'alphamax', 'alpha_fast'], index=[0])
if printout == True:
print(df1,df2)
return(np.round(df1, digits), np.round(df2, digits), ghat)
In [33]:
%timeit -r3 -n10 fsr_fast_pv(pval,x.shape[1])
10 loops, best of 3: 8.17 ms per loop
In [34]:
%timeit -r3 -n10 fsr_fast_pv_cython(pval,x.shape[1])
10 loops, best of 3: 7.17 ms per loop
Content source: CeciliaShi/STA-663-Final-Project
Similar notebooks: