Use gsl brent solvergsl

slower than the jpm's one

5 files changed, 80 insertions, 51 deletions

diff --git a/c_layer/cdsbootstrap.c b/c_layer/cdsbootstrap.c
index f3ae975..6b0a4cd 100644
--- a/c_layer/cdsbootstrap.c
+++ b/c_layer/cdsbootstrap.c
@@ -1,12 +1,6 @@
 #include "cdsbootstrap.h"
 
-int cdsBootstrapPointFunction
-(double   hazardRate,
- void    *data,
- double  *pv)
-{
-    int         status    = FAILURE;
-
+double cdsBootstrapPointFunction(double   hazardRate, void    *data) {
     cds_bootstrap_ctx *context = (cds_bootstrap_ctx*)data;
 
     TCurve         *discountCurve = context->discountCurve;
@@ -27,7 +21,6 @@ int cdsBootstrapPointFunction
                                cdsCurve,
                                context->recoveryRate,
                                &pvC) != SUCCESS)
-        goto done;
 
     if (JpmcdsFeeLegPV(context->fl,
                        cdsBaseDate,
@@ -37,10 +30,8 @@ int cdsBootstrapPointFunction
                        cdsCurve,
                        1,
                        &pvF) != SUCCESS)
-        goto done;
     /* Note: price is discounted to cdsBaseDate */
-    *pv = pvC - context->spread * pvF;
-    status = SUCCESS;
+    return pvC - context->spread * pvF;
 
  done:
 
diff --git a/c_layer/cdsbootstrap.h b/c_layer/cdsbootstrap.h
index 2e19189..5a50962 100644
--- a/c_layer/cdsbootstrap.h
+++ b/c_layer/cdsbootstrap.h
@@ -21,6 +21,5 @@ typedef struct
 } cds_bootstrap_ctx;
 
 
-int cdsBootstrapPointFunction(double   hazardRate,
-                              void    *data,
-                              double  *pv);
+double cdsBootstrapPointFunction(double   hazardRate,
+                                 void    *data);
diff --git a/pyisda/flat_hazard.pyx b/pyisda/flat_hazard.pyx
index daeac41..2bcbda3 100644
--- a/pyisda/flat_hazard.pyx
+++ b/pyisda/flat_hazard.pyx
@@ -11,6 +11,8 @@ from date cimport pydate_to_TDate, ACT_360, ACT_365F
 from cdsone cimport JpmcdsStringToStubMethod, TStubMethod
 from libc.stdlib cimport free, malloc
 from libc.stdio cimport printf
+from libc.math cimport fabs
+cimport gsl
 cimport cython
 
 np.import_array()
@@ -79,7 +81,7 @@ def strike_vec(double[:] spreads, YieldCurve yc, trade_date, value_date,
         JpmcdsFreeTCurve(sc)
     return default_leg_pv, coupon_leg_pv
 
-cdef extern from "cdsbootstrap.h":
+cdef extern from "cdsbootstrap.h" nogil:
     ctypedef struct cds_bootstrap_ctx:
         TDate           stepinDate
         TDate           cashSettleDate
@@ -89,23 +91,23 @@ cdef extern from "cdsbootstrap.h":
         double          spread
         TContingentLeg *cl
         TFeeLeg        *fl
-    int cdsBootstrapPointFunction(...) nogil
+    double cdsBootstrapPointFunction(double, void*)
 
-cdef extern from "isda/rtbrent.h" nogil:
-    ctypedef int (*TObjectFunc) (double x, void * para, double *f)
+# cdef extern from "isda/rtbrent.h" nogil:
+#     ctypedef int (*TObjectFunc) (double x, void * para, double *f)
 
-    int JpmcdsRootFindBrent(
-        TObjectFunc funcd,              # (I) function to be solved */
-        void        *data,              # (I) data to pass into funcd */
-        double      boundLo,            # (I) lower bound on legal X */
-        double      boundHi,            # (I) upper bound on legal X */
-        int         numIterations,      # (I) Maximum number of iterations */
-        double      guess,              # (I) Initial guess */
-        double      initalXStep,        # (I) Size of step in x */
-        double      initialFDeriv,      # (I) Derivative, defaults to 0 */
-        double      xacc,               # (I) X accuracy tolerance */
-        double      facc,               # (I) function accuracy tolerance */
-        double      *solution)          # (O) root found */
+#     int JpmcdsRootFindBrent(
+#         TObjectFunc funcd,              # (I) function to be solved */
+#         void        *data,              # (I) data to pass into funcd */
+#         double      boundLo,            # (I) lower bound on legal X */
+#         double      boundHi,            # (I) upper bound on legal X */
+#         int         numIterations,      # (I) Maximum number of iterations */
+#         double      guess,              # (I) Initial guess */
+#         double      initalXStep,        # (I) Size of step in x */
+#         double      initialFDeriv,      # (I) Derivative, defaults to 0 */
+#         double      xacc,               # (I) X accuracy tolerance */
+#         double      facc,               # (I) function accuracy tolerance */
+#         double      *solution)          # (O) root found */
 
 @cython.boundscheck(False)
 @cython.cdivision(True)
@@ -119,17 +121,20 @@ def pv_vec(double[:] spreads, YieldCurve yc, trade_date, value_date,
     cdef TDate start_date_c = pydate_to_TDate(start_date)
     cdef TDate end_date_c = pydate_to_TDate(end_date)
     cdef np.npy_intp nspreads = spreads.shape[0]
-    cdef double h, guess, hi, lo, coupon_leg_pv
+    cdef double h, hi, lo, coupon_leg_pv
     cdef np.ndarray[np.float_t] pv = np.empty_like(spreads)
     cdef TStubMethod stub_type
     cdef TCurve* sc
     cdef TContingentLeg* default_leg
     cdef TFeeLeg* fee_leg
     cdef cds_bootstrap_ctx* params
-    cdef Py_ssize_t i
+    cdef Py_ssize_t i, j
     JpmcdsStringToStubMethod(b"f/s", &stub_type)
     cdef double eta = 365./360.
-
+    cdef gsl.root_fsolver* s = gsl.root_fsolver_alloc(gsl.gsl_root_fsolver_brent)
+    cdef gsl.function f
+    f.function = &cdsBootstrapPointFunction
+    cdef int status
     with nogil:
         h = spreads[0] / (1 - recovery_rate)
         sc = JpmcdsMakeTCurve(trade_date_c,
@@ -162,31 +167,34 @@ def pv_vec(double[:] spreads, YieldCurve yc, trade_date, value_date,
         params.recoveryRate = recovery_rate
         params.cl = default_leg
         params.fl = fee_leg
-
+        f.params = <void*>params
         for i in range(nspreads):
             params.spread = spreads[i]
-            lo = guess * 0.9
-            guess = spreads[i] / ( 1 - recovery_rate) * eta
-            hi = guess * 1.1
-            if JpmcdsRootFindBrent(<TObjectFunc>cdsBootstrapPointFunction,
-                                   <void*> params,
-                                   lo,   # boundLo */
-                                   hi,   # boundHi */
-                                   100,    # numIterations */
-                                   guess,
-                                   0.0005, # initialXstep */
-                                   0,      # initialFDeriv */
-                                   1e-10,  # xacc */
-                                   1e-10,  # facc */
-                                   &h) != 0:
-                printf("failed to find the root")
-            sc.fArray[0].fRate = h
+            lo = spreads[i] / ( 1 - recovery_rate) * eta * 0.99
+            hi = spreads[i] / ( 1 - recovery_rate) * eta * 1.01
+            if gsl.root_fsolver_set(s, &f, lo, hi) == gsl.GSL_EINVAL:
+                break
+            j = 0
+            while j < 100:
+                j += 1
+                gsl.root_fsolver_iterate(s)
+                status = gsl.root_test_interval(gsl.root_fsolver_x_lower(s),
+                                                gsl.root_fsolver_x_upper(s),
+                                                1e-8,
+                                                1e-8)
+                if status == gsl.GSL_SUCCESS:
+                    break
+            sc.fArray[0].fRate = gsl.root_fsolver_root(s)
             if JpmcdsFeeLegPV(fee_leg, trade_date_c, step_in_date_c, value_date_c,
                               yc._thisptr, sc, True, &coupon_leg_pv) != 0:
                 printf("Something went wrong")
             pv[i] = coupon_leg_pv * (spreads[i] - fixed_rate)
+
+
         free(params)
         JpmcdsFeeLegFree(fee_leg)
         free(default_leg)
         JpmcdsFreeTCurve(sc)
+        gsl.root_fsolver_free(s)
+
     return pv
diff --git a/pyisda/gsl.pxd b/pyisda/gsl.pxd
new file mode 100644
index 0000000..43ee9d5
--- /dev/null
+++ b/pyisda/gsl.pxd
@@ -0,0 +1,30 @@
+cdef extern from "gsl/gsl_roots.h" nogil:
+    ctypedef struct root_fsolver "gsl_root_fsolver":
+        pass
+    ctypedef struct gsl_root_fsolver_type:
+        pass
+    root_fsolver* root_fsolver_alloc "gsl_root_fsolver_alloc" (const gsl_root_fsolver_type * T)
+    void root_fsolver_free "gsl_root_fsolver_free" (root_fsolver * s)
+    int root_fsolver_set "gsl_root_fsolver_set"(root_fsolver * s,
+                                           function * f,
+                                           double x_lower, double x_upper)
+    int root_fsolver_iterate "gsl_root_fsolver_iterate" (root_fsolver * s)
+    const char * root_fsolver_name "gsl_root_fsolver_name" (const root_fsolver * s)
+    double root_fsolver_root "gsl_root_fsolver_root" (const root_fsolver * s)
+    double root_fsolver_x_lower "gsl_root_fsolver_x_lower" (const root_fsolver * s)
+    double root_fsolver_x_upper "gsl_root_fsolver_x_upper" (const root_fsolver * s)
+    const gsl_root_fsolver_type* gsl_root_fsolver_bisection
+    const gsl_root_fsolver_type* gsl_root_fsolver_brent
+    const gsl_root_fsolver_type* gsl_root_fsolver_falsepos
+    int root_test_interval "gsl_root_test_interval" (double x_lower, double x_upper, double epsabs, double epsrel)
+
+cdef extern from "gsl/gsl_math.h" nogil:
+    ctypedef struct function "gsl_function":
+        double (* function) (double x, void * params)
+        void *params
+
+cdef extern from "gsl/gsl_errno.h" nogil:
+    cdef enum:
+        GSL_EINVAL = 4
+        GSL_SUCCESS = 0
+        GSL_CONTINUE = -2
diff --git a/setup.py b/setup.py
index 5958dd8..54336ab 100644
--- a/setup.py
+++ b/setup.py
@@ -9,7 +9,8 @@ all_extensions = Extension("*", ["pyisda/*.pyx"],
 c_extension = Extension("pyisda.flat_hazard",
                         include_dirs = ['c_layer'],
                         sources = ['pyisda/flat_hazard.pyx', 'c_layer/cdsbootstrap.c'],
-                        libraries = ['cds'])
+                        libraries = ['cds', 'gsl', 'cblas'],
+                        extra_compile_args = ['-O3', '-march=native'])
 all_extensions = cythonize([c_extension, all_extensions], nthreads = 4,
                            compiler_directives={'embedsignature':True})