src/common/muc_randvar.c

#ifndef lint
static const char RCSid[] = "$Id$";
#endif
#include <stdio.h>
#include <stdlib.h>

#include "fvect.h"
#include "muc_randvar.h"

int             muc_rvar_init(struct muc_rvar* rv)
{
        rv->n = 0;
        rv->w = 0.0;
        if (!(rv->samples = g3fl_create(1)))
                return 0;
        rv->sum = NULL;
        rv->sum_sqr = NULL;
        rv->min = NULL;
        rv->max = NULL;
        rv->save_samp = 1;
        return muc_rvar_set_dim(rv,1);
}

struct  muc_rvar*       muc_rvar_create()
{
        struct muc_rvar* res;

        if (!(res = (struct  muc_rvar*) malloc(sizeof(struct  muc_rvar))))
                return NULL;
        if (!(muc_rvar_init(res))) {
                free(res);
                return NULL;
        }
        return res;
}

void    muc_rvar_free(struct muc_rvar* rv)
{
        free(rv->sum);
        free(rv->sum_sqr);
        free(rv->min);
        free(rv->max);
        g3fl_free(rv->samples);
        free(rv);
}

int             muc_rvar_set_dim(struct muc_rvar* rv,int dim)
{
        if (rv->n != 0) {
                fprintf(stderr,"muc_rvar_set_dim: variable not cleared\n");
                return 0;
        }
        if (!(rv->sum = realloc(rv->sum,dim*sizeof(double))))
                return 0;
        if (!(rv->sum_sqr = realloc(rv->sum_sqr,dim*sizeof(double))))
                return 0;
        if (!(rv->min = realloc(rv->min,dim*sizeof(double))))
                return 0;
        if (!(rv->max = realloc(rv->max,dim*sizeof(double))))
                return 0;
        rv->samples->comp_size = dim;
        muc_rvar_clear(rv);
        return 1;
}

void    muc_rvar_clear(struct muc_rvar* rv)
{
        int i;

        rv->samples->size=0;
        rv->w = 0.0;
        rv->n = 0;
        for(i=0;i<muc_rvar_get_dim(rv);i++) {
                rv->sum[i] = 0.0;
                rv->sum_sqr[i] = 0.0;
                rv->min[i] = FHUGE;
                rv->max[i] = -FHUGE;
        }
}

void    muc_rvar_no_samples(struct muc_rvar* rv)
{
        rv->samples->size = 0;
        rv->save_samp = 0;
}

int             muc_rvar_store_sample(struct muc_rvar* rv, int on)
{
        if (rv->n != 0) {
                fprintf(stderr,"muc_rvar_store_sample: clear samples first\n");
                return 0;
        }
        rv->save_samp = on;
        return 1;
}

int             muc_rvar_add_sample(struct muc_rvar* rv,double* s)
{
        return muc_rvar_add_weighted_sample(rv,1.0,s);
}

int     muc_rvar_add_weighted_sample(struct muc_rvar* rv,double w,double* s)
{
        int i;
        double val;

        if (rv->save_samp)
                if (!(g3fl_append(rv->samples,s)))
                        return 0;
        for(i=0;i<muc_rvar_get_dim(rv);i++) {
                val = s[i]*w;
                rv->sum[i] += val;
                rv->sum_sqr[i] += s[i]*s[i]*w;
                if (rv->min[i] > val)
                        rv->min[i] = val;
                if (rv->max[i] < val)
                        rv->max[i] = val;
        }
        rv->w += w;
        rv->n++;
        return 1;
}
        
int             muc_rvar_get_vx(struct muc_rvar* rv,double* vx)
{
        int i;
        double ex,ex2;

        if (rv->w == 0.0) {
                return 0;
        }
        for(i=0;i<muc_rvar_get_dim(rv);i++) {
                ex = rv->sum[i]/rv->w;
                ex2 = rv->sum_sqr[i]/rv->w;
                vx[i] = ex2 - ex*ex;
        }
        return 1;
}

int             muc_rvar_get_ex(struct muc_rvar* rv,double* ex)
{
        int i;

        if (rv->w == 0.0) {
                return 0;
        }
        for(i=0;i<muc_rvar_get_dim(rv);i++) {
                ex[i] = rv->sum[i]/rv->w;
        }
        return 1;
}
int             muc_rvar_get_sum(struct muc_rvar* rv,double* sum2)
{
        int i;

        for(i=0;i<muc_rvar_get_dim(rv);i++) {
                sum2[i] = rv->sum[i];
        }
        return 1;
}

int             muc_rvar_get_median(struct muc_rvar* rv, double* median)
{
        double val;
        int i;

        if (!(rv->save_samp)) {
                fprintf(stderr,"muc_rvar_get_median: samples are not stored\n");
                return 0;
        }
        for(i=0;i<muc_rvar_get_dim(rv);i++) {
                g3fl_sort(rv->samples, i);
                val = g3fl_get(rv->samples, rv->n/2)[i];
                
                if (rv->n % 2 == 0) {
                        val += g3fl_get(rv->samples,rv->n/2 - 1)[i];
                        val /= 2.0;
                }
                median[i] = val;
        }
        return 1;
}
int             muc_rvar_get_percentile(struct muc_rvar* rv, double* median, double percentile)
{
        double val;
        int i;

        if (!(rv->save_samp)) {
                fprintf(stderr,"muc_rvar_get_median: samples are not stored\n");
                return 0;
        }
        for(i=0;i<muc_rvar_get_dim(rv);i++) {
                g3fl_sort(rv->samples, i);
                val = g3fl_get(rv->samples, rv->n*percentile)[i];
                
                if (rv->n % 1/percentile == 0) {
                        val += g3fl_get(rv->samples,rv->n*percentile - 1)[i];
                        val /= 2.0;
                }
                median[i] = val;
        }
        return 1;
}

int             muc_rvar_get_bounding_box(struct muc_rvar* rv, double* bounds)
{
        int i;

        if (rv->n == 0) 
                return 0;
        for(i=0;i<muc_rvar_get_dim(rv);i++) {
                bounds[2*i] = rv->min[i];
                bounds[2*i + 1] = rv->max[i];
        }
        return 1;
}

Revision:	2.1
Committed:	Thu Jul 14 17:32:12 2016 UTC (7 years, 10 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R4, rad5R2, rad5R3, rad5R1, HEAD
Log Message:	Changes to evalglare by Jan Wienold
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2			static const char RCSid[] = "$Id$";
3			#endif
4			#include <stdio.h>
5			#include <stdlib.h>
6
7			#include "fvect.h"
8			#include "muc_randvar.h"
9
10			int muc_rvar_init(struct muc_rvar* rv)
11			{
12			rv->n = 0;
13			rv->w = 0.0;
14			if (!(rv->samples = g3fl_create(1)))
15			return 0;
16			rv->sum = NULL;
17			rv->sum_sqr = NULL;
18			rv->min = NULL;
19			rv->max = NULL;
20			rv->save_samp = 1;
21			return muc_rvar_set_dim(rv,1);
22			}
23
24			struct muc_rvar* muc_rvar_create()
25			{
26			struct muc_rvar* res;
27
28			if (!(res = (struct muc_rvar*) malloc(sizeof(struct muc_rvar))))
29			return NULL;
30			if (!(muc_rvar_init(res))) {
31			free(res);
32			return NULL;
33			}
34			return res;
35			}
36
37			void muc_rvar_free(struct muc_rvar* rv)
38			{
39			free(rv->sum);
40			free(rv->sum_sqr);
41			free(rv->min);
42			free(rv->max);
43			g3fl_free(rv->samples);
44			free(rv);
45			}
46
47			int muc_rvar_set_dim(struct muc_rvar* rv,int dim)
48			{
49			if (rv->n != 0) {
50			fprintf(stderr,"muc_rvar_set_dim: variable not cleared\n");
51			return 0;
52			}
53			if (!(rv->sum = realloc(rv->sum,dim*sizeof(double))))
54			return 0;
55			if (!(rv->sum_sqr = realloc(rv->sum_sqr,dim*sizeof(double))))
56			return 0;
57			if (!(rv->min = realloc(rv->min,dim*sizeof(double))))
58			return 0;
59			if (!(rv->max = realloc(rv->max,dim*sizeof(double))))
60			return 0;
61			rv->samples->comp_size = dim;
62			muc_rvar_clear(rv);
63			return 1;
64			}
65
66			void muc_rvar_clear(struct muc_rvar* rv)
67			{
68			int i;
69
70			rv->samples->size=0;
71			rv->w = 0.0;
72			rv->n = 0;
73			for(i=0;i<muc_rvar_get_dim(rv);i++) {
74			rv->sum[i] = 0.0;
75			rv->sum_sqr[i] = 0.0;
76			rv->min[i] = FHUGE;
77			rv->max[i] = -FHUGE;
78			}
79			}
80
81			void muc_rvar_no_samples(struct muc_rvar* rv)
82			{
83			rv->samples->size = 0;
84			rv->save_samp = 0;
85			}
86
87			int muc_rvar_store_sample(struct muc_rvar* rv, int on)
88			{
89			if (rv->n != 0) {
90			fprintf(stderr,"muc_rvar_store_sample: clear samples first\n");
91			return 0;
92			}
93			rv->save_samp = on;
94			return 1;
95			}
96
97			int muc_rvar_add_sample(struct muc_rvar* rv,double* s)
98			{
99			return muc_rvar_add_weighted_sample(rv,1.0,s);
100			}
101
102			int muc_rvar_add_weighted_sample(struct muc_rvar* rv,double w,double* s)
103			{
104			int i;
105			double val;
106
107			if (rv->save_samp)
108			if (!(g3fl_append(rv->samples,s)))
109			return 0;
110			for(i=0;i<muc_rvar_get_dim(rv);i++) {
111			val = s[i]*w;
112			rv->sum[i] += val;
113			rv->sum_sqr[i] += s[i]s[i]w;
114			if (rv->min[i] > val)
115			rv->min[i] = val;
116			if (rv->max[i] < val)
117			rv->max[i] = val;
118			}
119			rv->w += w;
120			rv->n++;
121			return 1;
122			}
123
124			int muc_rvar_get_vx(struct muc_rvar* rv,double* vx)
125			{
126			int i;
127			double ex,ex2;
128
129			if (rv->w == 0.0) {
130			return 0;
131			}
132			for(i=0;i<muc_rvar_get_dim(rv);i++) {
133			ex = rv->sum[i]/rv->w;
134			ex2 = rv->sum_sqr[i]/rv->w;
135			vx[i] = ex2 - ex*ex;
136			}
137			return 1;
138			}
139
140			int muc_rvar_get_ex(struct muc_rvar* rv,double* ex)
141			{
142			int i;
143
144			if (rv->w == 0.0) {
145			return 0;
146			}
147			for(i=0;i<muc_rvar_get_dim(rv);i++) {
148			ex[i] = rv->sum[i]/rv->w;
149			}
150			return 1;
151			}
152			int muc_rvar_get_sum(struct muc_rvar* rv,double* sum2)
153			{
154			int i;
155
156			for(i=0;i<muc_rvar_get_dim(rv);i++) {
157			sum2[i] = rv->sum[i];
158			}
159			return 1;
160			}
161
162			int muc_rvar_get_median(struct muc_rvar* rv, double* median)
163			{
164			double val;
165			int i;
166
167			if (!(rv->save_samp)) {
168			fprintf(stderr,"muc_rvar_get_median: samples are not stored\n");
169			return 0;
170			}
171			for(i=0;i<muc_rvar_get_dim(rv);i++) {
172			g3fl_sort(rv->samples, i);
173			val = g3fl_get(rv->samples, rv->n/2)[i];
174
175			if (rv->n % 2 == 0) {
176			val += g3fl_get(rv->samples,rv->n/2 - 1)[i];
177			val /= 2.0;
178			}
179			median[i] = val;
180			}
181			return 1;
182			}
183			int muc_rvar_get_percentile(struct muc_rvar* rv, double* median, double percentile)
184			{
185			double val;
186			int i;
187
188			if (!(rv->save_samp)) {
189			fprintf(stderr,"muc_rvar_get_median: samples are not stored\n");
190			return 0;
191			}
192			for(i=0;i<muc_rvar_get_dim(rv);i++) {
193			g3fl_sort(rv->samples, i);
194			val = g3fl_get(rv->samples, rv->n*percentile)[i];
195
196			if (rv->n % 1/percentile == 0) {
197			val += g3fl_get(rv->samples,rv->n*percentile - 1)[i];
198			val /= 2.0;
199			}
200			median[i] = val;
201			}
202			return 1;
203			}
204
205			int muc_rvar_get_bounding_box(struct muc_rvar* rv, double* bounds)
206			{
207			int i;
208
209			if (rv->n == 0)
210			return 0;
211			for(i=0;i<muc_rvar_get_dim(rv);i++) {
212			bounds[2*i] = rv->min[i];
213			bounds[2*i + 1] = rv->max[i];
214			}
215			return 1;
216			}
217