ssp_ranst_gaussian.c (6103B)
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/mem_allocator.h> 19 #include <rsys/ref_count.h> 20 21 using normal_dist = RAN_NAMESPACE::normal_distribution<double>; 22 23 using normal_dist_float = RAN_NAMESPACE::normal_distribution<float>; 24 25 struct ssp_ranst_gaussian { 26 ref_T ref; 27 struct mem_allocator* allocator; 28 normal_dist* distrib; 29 double mu; 30 double K1; /* 1.0 / sigma */ 31 double K2; /* 1.0 / (sigma * SQRT_2_PI) */ 32 }; 33 34 struct ssp_ranst_gaussian_float { 35 ref_T ref; 36 struct mem_allocator* allocator; 37 normal_dist_float* distrib; 38 float mu; 39 float K1; /* 1.0 / sigma */ 40 float K2; /* 1.0 / (sigma * SQRT_2_PI) */ 41 }; 42 43 /******************************************************************************* 44 * Helper function 45 ******************************************************************************/ 46 static void 47 gaussian_release(ref_T* ref) 48 { 49 ssp_ranst_gaussian* ran; 50 ASSERT(ref); 51 ran = CONTAINER_OF(ref, ssp_ranst_gaussian, ref); 52 if(ran->distrib) { 53 ran->distrib->~normal_dist(); 54 MEM_RM(ran->allocator, ran->distrib); 55 } 56 MEM_RM(ran->allocator, ran); 57 } 58 59 static void 60 gaussian_float_release(ref_T* ref) 61 { 62 ssp_ranst_gaussian_float* ran; 63 ASSERT(ref); 64 ran = CONTAINER_OF(ref, ssp_ranst_gaussian_float, ref); 65 if(ran->distrib) { 66 ran->distrib->~normal_dist_float(); 67 MEM_RM(ran->allocator, ran->distrib); 68 } 69 MEM_RM(ran->allocator, ran); 70 } 71 72 /******************************************************************************* 73 * Exported functions 74 ******************************************************************************/ 75 res_T 76 ssp_ranst_gaussian_create 77 (struct mem_allocator* allocator, 78 struct ssp_ranst_gaussian** out_ran) 79 { 80 ssp_ranst_gaussian* ran = nullptr; 81 void* mem = nullptr; 82 res_T res = RES_OK; 83 84 if(!out_ran) 85 return RES_BAD_ARG; 86 87 allocator = allocator ? allocator : &mem_default_allocator; 88 89 ran = static_cast<ssp_ranst_gaussian*> 90 (MEM_CALLOC(allocator, 1, sizeof(struct ssp_ranst_gaussian))); 91 if(!ran) { 92 res = RES_MEM_ERR; 93 goto error; 94 } 95 ref_init(&ran->ref); 96 97 mem = MEM_ALLOC(allocator, sizeof(normal_dist)); 98 if(!mem) { 99 res = RES_MEM_ERR; 100 goto error; 101 } 102 ran->allocator = allocator; 103 ran->distrib = static_cast<normal_dist*>(new(mem) normal_dist); 104 ran->K1 = -1; /* invalid */ 105 if(out_ran) *out_ran = ran; 106 107 exit: 108 return res; 109 error: 110 if(ran) { 111 SSP(ranst_gaussian_ref_put(ran)); 112 ran = nullptr; 113 } 114 goto exit; 115 } 116 117 res_T 118 ssp_ranst_gaussian_float_create 119 (struct mem_allocator* allocator, 120 struct ssp_ranst_gaussian_float** out_ran) 121 { 122 ssp_ranst_gaussian_float* ran = nullptr; 123 void* mem = nullptr; 124 res_T res = RES_OK; 125 126 if(!out_ran) 127 return RES_BAD_ARG; 128 129 allocator = allocator ? allocator : &mem_default_allocator; 130 131 ran = static_cast<ssp_ranst_gaussian_float*> 132 (MEM_CALLOC(allocator, 1, sizeof(struct ssp_ranst_gaussian_float))); 133 if(!ran) { 134 res = RES_MEM_ERR; 135 goto error; 136 } 137 ref_init(&ran->ref); 138 139 mem = MEM_ALLOC(allocator, sizeof(normal_dist_float)); 140 if(!mem) { 141 res = RES_MEM_ERR; 142 goto error; 143 } 144 ran->allocator = allocator; 145 ran->distrib = static_cast<normal_dist_float*>(new(mem) normal_dist_float); 146 ran->K1 = -1; /* invalid */ 147 if (out_ran) *out_ran = ran; 148 149 exit: 150 return res; 151 error: 152 if(ran) { 153 SSP(ranst_gaussian_float_ref_put(ran)); 154 ran = nullptr; 155 } 156 goto exit; 157 } 158 159 res_T 160 ssp_ranst_gaussian_setup 161 (struct ssp_ranst_gaussian* ran, 162 const double mu, 163 const double sigma) 164 { 165 if(!ran || sigma < 0) 166 return RES_BAD_ARG; 167 168 normal_dist::param_type p{mu, sigma}; 169 ran->distrib->param(p); 170 ran->mu = mu; 171 ran->K1 = 1 / sigma; 172 ran->K2 = 1 / (sigma * SQRT_2_PI); 173 174 return RES_OK; 175 } 176 177 res_T 178 ssp_ranst_gaussian_float_setup 179 (struct ssp_ranst_gaussian_float* ran, 180 const float mu, 181 const float sigma) 182 { 183 if(!ran || sigma < 0) 184 return RES_BAD_ARG; 185 186 normal_dist_float::param_type p{ mu, sigma }; 187 ran->distrib->param(p); 188 ran->mu = mu; 189 ran->K1 = 1 / sigma; 190 ran->K2 = 1 / (sigma * (float)SQRT_2_PI); 191 192 return RES_OK; 193 } 194 195 res_T 196 ssp_ranst_gaussian_ref_get 197 (struct ssp_ranst_gaussian* ran) 198 { 199 if(!ran) return RES_BAD_ARG; 200 ref_get(&ran->ref); 201 return RES_OK; 202 } 203 204 res_T 205 ssp_ranst_gaussian_float_ref_get 206 (struct ssp_ranst_gaussian_float* ran) 207 { 208 if(!ran) return RES_BAD_ARG; 209 ref_get(&ran->ref); 210 return RES_OK; 211 } 212 213 res_T 214 ssp_ranst_gaussian_ref_put 215 (struct ssp_ranst_gaussian* ran) 216 { 217 if(!ran) return RES_BAD_ARG; 218 ref_put(&ran->ref, gaussian_release); 219 return RES_OK; 220 } 221 222 res_T 223 ssp_ranst_gaussian_float_ref_put 224 (struct ssp_ranst_gaussian_float* ran) 225 { 226 if(!ran) return RES_BAD_ARG; 227 ref_put(&ran->ref, gaussian_float_release); 228 return RES_OK; 229 } 230 231 double 232 ssp_ranst_gaussian_get 233 (const struct ssp_ranst_gaussian* ran, struct ssp_rng* rng) 234 { 235 ASSERT(ran->K1 > 0); 236 return wrap_ran(*rng, *ran->distrib); 237 } 238 239 float 240 ssp_ranst_gaussian_float_get 241 (const struct ssp_ranst_gaussian_float* ran, struct ssp_rng* rng) 242 { 243 ASSERT(ran->K1 > 0); 244 return wrap_ran(*rng, *ran->distrib); 245 } 246 247 double 248 ssp_ranst_gaussian_pdf 249 (const struct ssp_ranst_gaussian* ran, const double x) 250 { 251 const double tmp = (x - ran->mu) * ran->K1; 252 ASSERT(ran->K1 > 0); 253 return ran->K2 * exp(-0.5 * tmp * tmp); 254 } 255 256 float 257 ssp_ranst_gaussian_float_pdf 258 (const struct ssp_ranst_gaussian_float* ran, const float x) 259 { 260 const float tmp = (x - ran->mu) * ran->K1; 261 ASSERT(ran->K1 > 0); 262 return ran->K2 * expf(-0.5f * tmp * tmp); 263 } 264