Example/bs (#755)

* EXAMPLE update example_bs.cpp

Author: densmirn

dpnp/backend/examples/example_bs.cpp

@@ -42,107 +42,161 @@
 void black_scholes(double* price,
                    double* strike,
                    double* t,
-                   double* mrs,
-                   double* vol_vol_twos,
-                   double* quarters,
-                   double* ones,
-                   double* halfs,
+                   const double rate,
+                   const double vol,
                    double* call,
                    double* put,
                    const size_t size)
 {
+    const size_t ndim = 1;
+    const size_t scalar_size = 1;
+
+    double* mr = (double*)dpnp_memory_alloc_c(1 * sizeof(double));
+    mr[0] = -rate;
+
+    double* vol_vol_two = (double*)dpnp_memory_alloc_c(1 * sizeof(double));
+    vol_vol_two[0] = vol * vol * 2;
+
+    double* quarter = (double*)dpnp_memory_alloc_c(1 * sizeof(double));
+    quarter[0] = 0.25;
+
+    double* one = (double*)dpnp_memory_alloc_c(1 * sizeof(double));
+    one[0] = 1.;
+
+    double* half = (double*)dpnp_memory_alloc_c(1 * sizeof(double));
+    half[0] = 0.5;
+
     double* P = price;
     double* S = strike;
     double* T = t;
 
-    double* PdivS = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_divide_c<double, double, double>(P, S, PdivS, size); // P / S
+    double* p_div_s = (double*)dpnp_memory_alloc_c(size * sizeof(double));
+    // p_div_s = P / S
+    dpnp_divide_c<double, double, double>(p_div_s, P, size, &size, ndim, S, size, &size, ndim, NULL);
     double* a = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_log_c<double, double>(PdivS, a, size); // np.log(P / S)
-    dpnp_memory_free_c(PdivS);
+    dpnp_log_c<double, double>(p_div_s, a, size); // a = np.log(p_div_s)
+    dpnp_memory_free_c(p_div_s);
 
     double* b = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_multiply_c<double, double, double>(T, mrs, b, size); // T * mrs
-
+    // b = T * mr
+    dpnp_multiply_c<double, double, double>(b, T, size, &size, ndim, mr, scalar_size, &scalar_size, ndim, NULL);
+    dpnp_memory_free_c(mr);
     double* z = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_multiply_c<double, double, double>(T, vol_vol_twos, z, size); // T * vol_vol_twos
+    // z = T * vol_vol_twos
+    dpnp_multiply_c<double, double, double>(z,
+                                            T, size, &size, ndim,
+                                            vol_vol_two, scalar_size, &scalar_size, ndim,
+                                            NULL);
+    dpnp_memory_free_c(vol_vol_two);
 
     double* c = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_multiply_c<double, double, double>(quarters, z, c, size); // quarters * z
+    // c = quarters * z
+    dpnp_multiply_c<double, double, double>(c, quarter, scalar_size, &scalar_size, ndim, z, size, &size, ndim, NULL);
+    dpnp_memory_free_c(quarter);
 
     double* sqrt_z = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_sqrt_c<double, double>(z, sqrt_z, size); // np.sqrt(z)
+    dpnp_sqrt_c<double, double>(z, sqrt_z, size); // sqrt_z = np.sqrt(z)
     dpnp_memory_free_c(z);
     double* y = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_divide_c<double, double, double>(ones, sqrt_z, y, size); // ones/np.sqrt(z)
+    // y = ones / np.sqrt(z)
+    dpnp_divide_c<double, double, double>(y, one, scalar_size, &scalar_size, ndim, sqrt_z, size, &size, ndim, NULL);
     dpnp_memory_free_c(sqrt_z);
+    dpnp_memory_free_c(one);
 
     double* a_sub_b = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_subtract_c<double, double, double>(a, b, a_sub_b, size); // a - b
+    // a_sub_b = a - b
+    dpnp_subtract_c<double, double, double>(a_sub_b, a, size, &size, ndim, b, size, &size, ndim, NULL);
     dpnp_memory_free_c(a);
     double* a_sub_b_add_c = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_add_c<double, double, double>(a_sub_b, c, a_sub_b_add_c, size); // a - b + c
+    // a_sub_b_add_c = a_sub_b + c
+    dpnp_add_c<double, double, double>(a_sub_b_add_c, a_sub_b, size, &size, ndim, c, size, &size, ndim, NULL);
     double* w1 = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_multiply_c<double, double, double>(a_sub_b_add_c, y, w1, size); // (a - b + c) * y
+    // w1 = a_sub_b_add_c * y
+    dpnp_multiply_c<double, double, double>(w1, a_sub_b_add_c, size, &size, ndim, y, size, &size, ndim, NULL);
     dpnp_memory_free_c(a_sub_b_add_c);
 
     double* a_sub_b_sub_c = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_subtract_c<double, double, double>(a_sub_b, c, a_sub_b_sub_c, size); // a - b - c
+    // a_sub_b_sub_c = a_sub_b - c
+    dpnp_subtract_c<double, double, double>(a_sub_b_sub_c, a_sub_b, size, &size, ndim, c, size, &size, ndim, NULL);
     dpnp_memory_free_c(a_sub_b);
     dpnp_memory_free_c(c);
     double* w2 = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_multiply_c<double, double, double>(a_sub_b_sub_c, y, w2, size); // (a - b - c) * y
+    // w2 = a_sub_b_sub_c * y
+    dpnp_multiply_c<double, double, double>(w2, a_sub_b_sub_c, size, &size, ndim, y, size, &size, ndim, NULL);
     dpnp_memory_free_c(a_sub_b_sub_c);
     dpnp_memory_free_c(y);
 
     double* erf_w1 = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_erf_c<double>(w1, erf_w1, size); // np.erf(w1)
+    dpnp_erf_c<double>(w1, erf_w1, size); // erf_w1 = np.erf(w1)
     dpnp_memory_free_c(w1);
     double* halfs_mul_erf_w1 = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_multiply_c<double, double, double>(halfs, erf_w1, halfs_mul_erf_w1, size); // halfs * np.erf(w1)
+    // halfs_mul_erf_w1 = halfs * erf_w1
+    dpnp_multiply_c<double, double, double>(halfs_mul_erf_w1,
+                                            half, scalar_size, &scalar_size, ndim,
+                                            erf_w1, size, &size, ndim,
+                                            NULL);
     dpnp_memory_free_c(erf_w1);
     double* d1 = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_add_c<double, double, double>(halfs, halfs_mul_erf_w1, d1, size); // halfs + halfs * np.erf(w1)
+    // d1 = halfs + halfs_mul_erf_w1
+    dpnp_add_c<double, double, double>(d1,
+                                       half, scalar_size, &scalar_size, ndim,
+                                       halfs_mul_erf_w1, size, &size, ndim,
+                                       NULL);
     dpnp_memory_free_c(halfs_mul_erf_w1);
 
     double* erf_w2 = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_erf_c<double>(w2, erf_w2, size); // np.erf(w2)
+    dpnp_erf_c<double>(w2, erf_w2, size); // erf_w2 = np.erf(w2)
     dpnp_memory_free_c(w2);
     double* halfs_mul_erf_w2 = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_multiply_c<double, double, double>(halfs, erf_w2, halfs_mul_erf_w2, size); // halfs * np.erf(w2)
+    // halfs_mul_erf_w2 = halfs * erf_w2
+    dpnp_multiply_c<double, double, double>(halfs_mul_erf_w2, 
+                                            half, scalar_size, &scalar_size, ndim,
+                                            erf_w2, size, &size, ndim,
+                                            NULL);
     dpnp_memory_free_c(erf_w2);
     double* d2 = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_add_c<double, double, double>(halfs, halfs_mul_erf_w2, d2, size); // halfs + halfs * np.erf(w2)
+    // d2 = halfs + halfs_mul_erf_w2
+    dpnp_add_c<double, double, double>(d2,
+                                       half, scalar_size, &scalar_size, ndim,
+                                       halfs_mul_erf_w2, size, &size, ndim,
+                                       NULL);
     dpnp_memory_free_c(halfs_mul_erf_w2);
+    dpnp_memory_free_c(half);
 
     double* exp_b = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_exp_c<double, double>(b, exp_b, size); // np.exp(b)
+    dpnp_exp_c<double, double>(b, exp_b, size); // exp_b = np.exp(b)
     double* Se = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_multiply_c<double, double, double>(exp_b, S, Se, size); // np.exp(b) * S
+    // Se = exp_b * S
+    dpnp_multiply_c<double, double, double>(Se, exp_b, size, &size, ndim, S, size, &size, ndim, NULL);
     dpnp_memory_free_c(exp_b);
     dpnp_memory_free_c(b);
 
-    double* Pmul_d1 = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_multiply_c<double, double, double>(P, d1, Pmul_d1, size); // P * d1
+    double* P_mul_d1 = (double*)dpnp_memory_alloc_c(size * sizeof(double));
+    // P_mul_d1 = P * d1
+    dpnp_multiply_c<double, double, double>(P_mul_d1, P, size, &size, ndim, d1, size, &size, ndim, NULL);
     dpnp_memory_free_c(d1);
     double* Se_mul_d2 = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_multiply_c<double, double, double>(Se, d2, Se_mul_d2, size); // Se * d2
+    // Se_mul_d2 = Se * d2
+    dpnp_multiply_c<double, double, double>(Se_mul_d2, Se, size, &size, ndim, d2, size, &size, ndim, NULL);
     dpnp_memory_free_c(d2);
     double* r = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_subtract_c<double, double, double>(Pmul_d1, Se_mul_d2, r, size); // P * d1 - Se * d2
+    // r = P_mul_d1 - Se_mul_d2
+    dpnp_subtract_c<double, double, double>(r, P_mul_d1, size, &size, ndim, Se_mul_d2, size, &size, ndim, NULL);
     dpnp_memory_free_c(Se_mul_d2);
-    dpnp_memory_free_c(Pmul_d1);
+    dpnp_memory_free_c(P_mul_d1);
 
     dpnp_copyto_c<double, double>(call, r, size); // call[:] = r
-    double* r_subP = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_subtract_c<double, double, double>(r, P, r_subP, size); // r - P
+    double* r_sub_P = (double*)dpnp_memory_alloc_c(size * sizeof(double));
+    // r_sub_P = r - P
+    dpnp_subtract_c<double, double, double>(r_sub_P, r, size, &size, ndim, P, size, &size, ndim, NULL);
     dpnp_memory_free_c(r);
-    double* r_subPaddSe = (double*)dpnp_memory_alloc_c(size * sizeof(double));
-    dpnp_add_c<double, double, double>(r_subP, Se, r_subPaddSe, size); // r - P + Se
-    dpnp_memory_free_c(r_subP);
+    double* r_sub_P_add_Se = (double*)dpnp_memory_alloc_c(size * sizeof(double));
+    // r_sub_P_add_Se = r_sub_P + Se
+    dpnp_add_c<double, double, double>(r_sub_P_add_Se, r_sub_P, size, &size, ndim, Se, size, &size, ndim, NULL);
+    dpnp_memory_free_c(r_sub_P);
     dpnp_memory_free_c(Se);
-    dpnp_copyto_c<double, double>(put, r_subPaddSe, size); // put[:] = r - P + Se
-    dpnp_memory_free_c(r_subPaddSe);
+    dpnp_copyto_c<double, double>(put, r_sub_P_add_Se, size); // put[:] = r_sub_P_add_Se
+    dpnp_memory_free_c(r_sub_P_add_Se);
 }
 
 int main(int, char**)
@@ -174,64 +228,30 @@ int main(int, char**)
     double* mone = (double*)dpnp_memory_alloc_c(1 * sizeof(double));
     mone[0] = -1.;
 
-    double* mr = (double*)dpnp_memory_alloc_c(1 * sizeof(double));
-    mr[0] = -RISK_FREE;
-
-    double* vol_vol_two = (double*)dpnp_memory_alloc_c(1 * sizeof(double));
-    vol_vol_two[0] = VOLATILITY * VOLATILITY * 2.;
-
-    double* quarter = (double*)dpnp_memory_alloc_c(1 * sizeof(double));
-    quarter[0] = 0.25;
-
-    double* one = (double*)dpnp_memory_alloc_c(1 * sizeof(double));
-    one[0] = 1.;
-
-    double* half = (double*)dpnp_memory_alloc_c(1 * sizeof(double));
-    half[0] = 0.5;
-
     double* call = (double*)dpnp_memory_alloc_c(SIZE * sizeof(double));
     double* put = (double*)dpnp_memory_alloc_c(SIZE * sizeof(double));
-    double* mrs = (double*)dpnp_memory_alloc_c(SIZE * sizeof(double));
-    double* vol_vol_twos = (double*)dpnp_memory_alloc_c(SIZE * sizeof(double));
-    double* quarters = (double*)dpnp_memory_alloc_c(SIZE * sizeof(double));
-    double* ones = (double*)dpnp_memory_alloc_c(SIZE * sizeof(double));
-    double* halfs = (double*)dpnp_memory_alloc_c(SIZE * sizeof(double));
-
-    dpnp_full_c<double>(zero, call, SIZE);                // np.full(SIZE, 0., dtype=DTYPE)
-    dpnp_full_c<double>(mone, put, SIZE);                 // np.full(SIZE, -1., dtype=DTYPE)
-    dpnp_full_c<double>(mr, mrs, SIZE);                   // np.full(SIZE, -RISK_FREE, dtype=DTYPE)
-    dpnp_full_c<double>(vol_vol_two, vol_vol_twos, SIZE); // np.full(SIZE, VOLATILITY * VOLATILITY * 2., dtype=DTYPE)
-    dpnp_full_c<double>(quarter, quarters, SIZE);         // np.full(SIZE, 0.25, dtype=DTYPE)
-    dpnp_full_c<double>(one, ones, SIZE);                 // np.full(SIZE, 1., dtype=DTYPE)
-    dpnp_full_c<double>(half, halfs, SIZE);               // np.full(SIZE, 0.5, dtype=DTYPE)
 
-    dpnp_memory_free_c(half);
-    dpnp_memory_free_c(one);
-    dpnp_memory_free_c(quarter);
-    dpnp_memory_free_c(vol_vol_two);
-    dpnp_memory_free_c(mr);
+    dpnp_full_c<double>(zero, call, SIZE);            // np.full(SIZE, 0., dtype=DTYPE)
+    dpnp_full_c<double>(mone, put, SIZE);             // np.full(SIZE, -1., dtype=DTYPE)
+
     dpnp_memory_free_c(mone);
     dpnp_memory_free_c(zero);
 
-    black_scholes(price, strike, t, mrs, vol_vol_twos, quarters, ones, halfs, call, put, SIZE);
+    black_scholes(price, strike, t, RISK_FREE, VOLATILITY, call, put, SIZE);
 
-    std::cout << std::endl;
+    std::cout << "call: ";
     for (size_t i = 0; i < 10; ++i)
     {
         std::cout << call[i] << ", ";
     }
-    std::cout << std::endl;
+    std::cout << "..." << std::endl;
+    std::cout << "put: ";
     for (size_t i = 0; i < 10; ++i)
     {
         std::cout << put[i] << ", ";
     }
-    std::cout << std::endl;
+    std::cout << "..." << std::endl;
 
-    dpnp_memory_free_c(halfs);
-    dpnp_memory_free_c(ones);
-    dpnp_memory_free_c(quarters);
-    dpnp_memory_free_c(vol_vol_twos);
-    dpnp_memory_free_c(mrs);
     dpnp_memory_free_c(put);
     dpnp_memory_free_c(call);
 

examples/example_bs.py

@@ -55,27 +55,30 @@
 VOLATILITY = 0.2
 
 
-def black_scholes(price, strike, t, mrs, vol_vol_twos, quarters, ones, halfs, call, put):
+def black_scholes(price, strike, t, rate, vol, call, put):
+    mr = -rate
+    sig_sig_two = vol * vol * 2
+
     P = price
     S = strike
     T = t
 
     a = np.log(P / S)
-    b = T * mrs
+    b = T * mr
 
-    z = T * vol_vol_twos
-    c = quarters * z
-    y = ones / np.sqrt(z)
+    z = T * sig_sig_two
+    c = 0.25 * z
+    y = 1./np.sqrt(z)
 
     w1 = (a - b + c) * y
     w2 = (a - b - c) * y
 
-    d1 = halfs + halfs * np.erf(w1)
-    d2 = halfs + halfs * np.erf(w2)
+    d1 = 0.5 + 0.5 * np.erf(w1)
+    d2 = 0.5 + 0.5 * np.erf(w2)
 
     Se = np.exp(b) * S
 
-    r = P * d1 - Se * d2
+    r =  P * d1 - Se * d2
     call[:] = r  # temporary `r` is necessary for faster `put` computation
     put[:] = r - P + Se
 
@@ -85,18 +88,9 @@ def black_scholes(price, strike, t, mrs, vol_vol_twos, quarters, ones, halfs, ca
 strike = np.random.uniform(SL, SH, SIZE)
 t = np.random.uniform(TL, TH, SIZE)
 
-call = np.full(SIZE, 0, dtype=DTYPE)
-put = np.full(SIZE, -1, dtype=DTYPE)
-
-mrs = np.full(SIZE, -RISK_FREE, dtype=DTYPE)
-
-vol = VOLATILITY
-vol_vol_twos = np.full(SIZE, vol * vol * 2, dtype=DTYPE)
-
-quarters = np.full(SIZE, 0.25, dtype=DTYPE)
-ones = np.full(SIZE, 1, dtype=DTYPE)
-halfs = np.full(SIZE, 0.5, dtype=DTYPE)
+call = np.zeros(SIZE, dtype=DTYPE)
+put  = -np.ones(SIZE, dtype=DTYPE)
 
-black_scholes(price, strike, t, mrs, vol_vol_twos, quarters, ones, halfs, call, put)
+black_scholes(price, strike, t, RISK_FREE, VOLATILITY, call, put)
 print(call[:10])
 print(put[:10])