star-sp

Random number generators and distributions
git clone git://git.meso-star.fr/star-sp.git
Log | Files | Refs | README | LICENSE

ssp_ranst_discrete.c (4136B)


      1 /* Copyright (C) 2015-2025 |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 "ssp_rng_c.h"
     17 
     18 #include <rsys/dynamic_array_double.h>
     19 #include <rsys/mem_allocator.h>
     20 #include <rsys/ref_count.h>
     21 
     22 using discrete_dist=RAN_NAMESPACE::discrete_distribution<size_t>;
     23 
     24 struct ssp_ranst_discrete {
     25   struct darray_double pdf;
     26   ref_T ref;
     27   struct mem_allocator* allocator;
     28   discrete_dist* distrib;
     29 };
     30 
     31 /*******************************************************************************
     32  * Helper functions
     33  ******************************************************************************/
     34 static void
     35 discrete_release(ref_T* ref)
     36 {
     37   struct ssp_ranst_discrete* ran;
     38   ASSERT(ref);
     39   ran = CONTAINER_OF(ref, struct ssp_ranst_discrete, ref);
     40   if(ran->distrib) {
     41     ran->distrib->~discrete_dist();
     42     MEM_RM(ran->allocator, ran->distrib);
     43   }
     44   darray_double_release(&ran->pdf);
     45   MEM_RM(ran->allocator, ran);
     46 }
     47 
     48 /*******************************************************************************
     49  * Exported functions
     50  ******************************************************************************/
     51 res_T
     52 ssp_ranst_discrete_create
     53   (struct mem_allocator* allocator,
     54    struct ssp_ranst_discrete** out_ran)
     55 {
     56   struct ssp_ranst_discrete* ran = nullptr;
     57   void* mem = nullptr;
     58   res_T res = RES_OK;
     59 
     60   if(!out_ran) {
     61     res = RES_BAD_ARG;
     62     goto error;
     63   }
     64   allocator = allocator ? allocator : &mem_default_allocator;
     65 
     66   ran = static_cast<ssp_ranst_discrete*>
     67     (MEM_CALLOC(allocator, 1, sizeof(ssp_ranst_discrete)));
     68   if(!ran) {
     69     res = RES_MEM_ERR;
     70     goto error;
     71   }
     72   ref_init(&ran->ref);
     73 
     74   mem = MEM_ALLOC(allocator, sizeof(discrete_dist));
     75   if(!mem) {
     76     res = RES_MEM_ERR;
     77     goto error;
     78   }
     79   ran->allocator = allocator;
     80   ran->distrib = static_cast<discrete_dist*>(new(mem) discrete_dist);
     81   darray_double_init(allocator, &ran->pdf);
     82 
     83 exit:
     84   if(out_ran) *out_ran = ran;
     85   return res;
     86 error:
     87   if(ran) {
     88     SSP(ranst_discrete_ref_put(ran));
     89     ran = NULL;
     90   }
     91   goto exit;
     92 }
     93 
     94 res_T
     95 ssp_ranst_discrete_setup
     96   (struct ssp_ranst_discrete* ran,
     97    const double* weights,
     98    const size_t nweights)
     99 {
    100   double* pdf;
    101   double sum = 0;
    102   size_t i;
    103   res_T res = RES_OK;
    104 
    105   if(!ran || !weights || !nweights) {
    106     res = RES_BAD_ARG;
    107     goto error;
    108   }
    109 
    110   res = darray_double_resize(&ran->pdf, nweights);
    111   if(res != RES_OK) goto error;
    112 
    113   pdf = darray_double_data_get(&ran->pdf);
    114 
    115   FOR_EACH(i, 0, nweights) {
    116     if(weights[i] < 0.f) {
    117       res = RES_BAD_ARG;
    118       goto error;
    119     }
    120     sum += weights[i];
    121     pdf[i] = weights[i];
    122   }
    123   if(sum == 0) {
    124     res = RES_BAD_ARG;
    125     goto error;
    126   }
    127 
    128   sum = 1 / sum;
    129   FOR_EACH(i, 0, nweights) {
    130     pdf[i] *= sum;
    131   }
    132 
    133   {
    134     discrete_dist::param_type p{weights, weights + nweights};
    135     ran->distrib->param(p);
    136   }
    137 
    138 exit:
    139   return res;
    140 error:
    141   if(ran) {
    142     darray_double_clear(&ran->pdf);
    143   }
    144   goto exit;
    145 }
    146 
    147 res_T
    148 ssp_ranst_discrete_ref_get(struct ssp_ranst_discrete* ran)
    149 {
    150   if(!ran) return RES_BAD_ARG;
    151   ref_get(&ran->ref);
    152   return RES_OK;
    153 }
    154 
    155 res_T
    156 ssp_ranst_discrete_ref_put(struct ssp_ranst_discrete* ran)
    157 {
    158   if(!ran) return RES_BAD_ARG;
    159   ref_put(&ran->ref, discrete_release);
    160   return RES_OK;
    161 }
    162 
    163 size_t
    164 ssp_ranst_discrete_get(struct ssp_rng* rng, const struct ssp_ranst_discrete* ran)
    165 {
    166   return wrap_ran(*rng, *ran->distrib);
    167 }
    168 
    169 double
    170 ssp_ranst_discrete_pdf(const size_t i, const struct ssp_ranst_discrete* ran)
    171 {
    172   ASSERT(ran && i < darray_double_size_get(&ran->pdf));
    173   return darray_double_cdata_get(&ran->pdf)[i];
    174 }
    175