star-wf

First-passage time of brownian motion
git clone git://git.meso-star.fr/star-wf.git
Log | Files | Refs | README | LICENSE

swf_HXd.h (3428B)

      1 /* Copyright (C) 2024 |Méso|Star> (contact@meso-star.com)
      2  *
      3  * This program is free software: you can redistribute it and/or modify
      4  * it under the terms of the GNU General Public License as published by
      5  * the Free Software Foundation, either version 3 of the License, or
      6  * (at your option) any later version.
      7  *
      8  * This program is distributed in the hope that it will be useful,
      9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     11  * GNU General Public License for more details.
     12  *
     13  * You should have received a copy of the GNU General Public License
     14  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     15 
     16 #include "swf.h"
     17 #include "swf_tabulation.h"
     18 
     19 #include <rsys/rsys.h>
     20 
     21 #ifndef SWF_DIMENSION
     22   #error "The SWF_DIMENSION macro must be defined"
     23 #endif
     24 
     25 #if SWF_DIMENSION != 2 && SWF_DIMENSION != 3
     26   #error "Invalid SWF_DIMEMSION"
     27 #endif
     28 
     29 #define Xd(Name) CONCAT(CONCAT(Name, SWF_DIMENSION), d)
     30 
     31 /*******************************************************************************
     32  * Helper generic dimension functions
     33  ******************************************************************************/
     34 static INLINE double
     35 Xd(H_inverse)
     36   (const double y,
     37    const double x_min,
     38    const double x_max,
     39    const double epsilon)
     40 {
     41   double x0 = x_min;
     42   double x1 = x_max;
     43   ASSERT(y >= 0 && y < 1 && x_min < x_max && epsilon > 0);
     44 
     45   if(y == 0) return 0;
     46 
     47   while(x1 - x0 > epsilon) {
     48     const double x = (x0 + x1) * 0.5;
     49     const double f_x = Xd(H)(x, NULL, NULL);
     50 
     51     if(f_x < y) {
     52       x0 = x;
     53     } else {
     54       x1 = x;
     55     }
     56   }
     57   return (x0 + x1) * 0.5;
     58 }
     59 
     60 static INLINE res_T
     61 Xd(H_tabulate)
     62   (const struct swf_H_tabulate_args* args,
     63    struct swf_tabulation** out_tab)
     64 {
     65   struct swf_tabulation* tab = NULL;
     66   size_t i = 0;
     67   double x = 0;
     68   res_T res = RES_OK;
     69 
     70   if(!out_tab) { res = RES_BAD_ARG; goto error; }
     71 
     72   if((res = check_swf_H_tabulate_args(args)) != RES_OK) goto error;
     73   if((res = tabulation_create(args->allocator, &tab)) != RES_OK) goto error;
     74 
     75   /* Setup the tabulation */
     76   i = 0;
     77   for(x = args->x_min; x <= args->x_max; x += x*args->step) {
     78     const struct item* pitem = darray_item_cdata_get(&tab->items) + i - 1;
     79     struct item item = ITEM_NULL;
     80 
     81     /* Evaluate H(x) and its derivatives up to order 2 */
     82     item.x = x;
     83     item.fx = Xd(H)(item.x, &item.dfx, &item.d2fx);
     84 
     85     if(i > 0) {
     86       /* Detect oscillations due to numerical problems.
     87        * These oscillations occur close to 0 */
     88       if(item.fx < pitem->fx) {
     89         res = RES_BAD_OP;
     90         goto error;
     91       }
     92       /* Do not tabulate argument whose value is (numerically) equal to the
     93        * previous one: it would artificially add staircase steps */
     94       if(item.fx == pitem->fx) {
     95         continue;
     96       }
     97     }
     98 
     99     ++i;
    100 
    101     res = darray_item_push_back(&tab->items, &item);
    102     if(res != RES_OK) goto error;
    103   }
    104   ASSERT(darray_item_size_get(&tab->items) == i);
    105 
    106   if(args->normalize) {
    107     struct item* items = darray_item_data_get(&tab->items);
    108     const size_t nitems = darray_item_size_get(&tab->items);
    109     const double norm = items[nitems - 1].fx;
    110 
    111     /* Normalize the tabulation */
    112     FOR_EACH(i, 0, nitems) items[i].fx /= norm;
    113   }
    114 
    115 exit:
    116   if(out_tab) *out_tab = tab;
    117   return res;
    118 error:
    119   if(tab) {
    120     SWF(tabulation_ref_put(tab));
    121     tab = NULL;
    122   }
    123   goto exit;
    124 }
    125 
    126 #undef SWF_DIMENSION
    127 #undef Xd