https://github.com/cran/kappalab
Tip revision: c35c3217a80758cae846a0c1336934ca8aee9b2b authored by Ivan Kojadinovic on 23 September 2006, 00:00:00 UTC
version 0.3-0
version 0.3-0
Tip revision: c35c321
core.c
/*##############################################################################
#
# Copyright © 2005 Michel Grabisch and Ivan Kojadinovic
#
# Ivan.Kojadinovic@polytech.univ-nantes.fr
#
# This software is a package for the statistical system GNU R:
# http://www.r-project.org
#
# This software is governed by the CeCILL license under French law and
# abiding by the rules of distribution of free software. You can use,
# modify and/ or redistribute the software under the terms of the CeCILL
# license as circulated by CEA, CNRS and INRIA at the following URL
# "http://www.cecill.info".
#
# As a counterpart to the access to the source code and rights to copy,
# modify and redistribute granted by the license, users are provided only
# with a limited warranty and the software's author, the holder of the
# economic rights, and the successive licensors have only limited
# liability.
#
# In this respect, the user's attention is drawn to the risks associated
# with loading, using, modifying and/or developing or reproducing the
# software by the user in light of its specific status of free software,
# that may mean that it is complicated to manipulate, and that also
# therefore means that it is reserved for developers and experienced
# professionals having in-depth computer knowledge. Users are therefore
# encouraged to load and test the software's suitability as regards their
# requirements in conditions enabling the security of their systems and/or
# data to be ensured and, more generally, to use and operate it in the
# same conditions as regards security.
#
# The fact that you are presently reading this means that you have had
# knowledge of the CeCILL license and that you accept its terms.
#
##############################################################################
*/
/*****************************************************************************
Non-additive measure and integral manipulation routines.
C functions source file. Binary coding of sets.
Univers X, X={0,1,2,3,...,n-1} with n=32 maximum.
Initial work: Michel Grabisch, 21/12/94
Contributions: Ivan Kojadinovic, 01/03/2005
*****************************************************************************/
#include <string.h>
#include <R.h>
#include "core.h"
/*****************************************************************************
Conversion list of elements -> binary
x : set to be converted
nx : number of elements of x
*****************************************************************************/
int subset2binary(int *x, int nx) {
int i;
int b;
b=0;
for(i=0; i<nx; i++)
b += 1 << x[i];
return(b);
}
/*****************************************************************************
Conversion binary -> list of elements
b : binary code to be converted
x : result in an array (ordered liste of elements)
*****************************************************************************/
void binary2subset(int n, int b, int *x) {
int i;
for(i=0; i<n; i++)
if(b & 1 << i) {
*x = i;
x++;
}
}
/*****************************************************************************
Pointer version of the binary2subset function, interfaced in kappalab.R
l: the real length of x
*****************************************************************************/
void binary2subsetR(int *n, int *b, int *x, int *l) {
int i;
*l=0;
for(i=0; i<*n; i++)
if(*b & 1 << i) {
*x = i+1; /* R arrays are 1-based */
x++;
(*l)++;
}
}
/*****************************************************************************
Sorting by tournament (M. Hebert)
n: dimension of vec and tourn
vec: vector to be sorted
iv: indices of the sorted elements (out)
*****************************************************************************/
void tri(register int n, register int *tourn, register double *vec,
register int *iv) {
int i;
register int k, kk;
double x;
int nb, nb1, nnb1, id;
tourn--;
vec--;
iv--;
/* adjonction dans le tournoi */
for(i=1; i<=n; i++) {
x = vec[i];
k = i;
while((k > 1) && (x < vec[tourn[(kk = k >> 1)]])) {
tourn[k] = tourn[kk];
k = kk;
}
tourn[k] = i;
}
for(i=1; i<=n; i++) {
nb = n - i + 1;
iv[i] = tourn[1] - 1;
/* extraction de la racine (min) */
id = tourn[nb];
k = 1;
nb1 = nb + 1;
nnb1 = nb1 >> 1;
while(k < nnb1)
if(vec[tourn[k<<1]] > vec[tourn[(k<<1)|1]]) {
tourn[k] = tourn[(k << 1) | 1];
k = (k << 1) | 1;
}
else {
tourn[k] = tourn[k << 1];
k = k << 1;
}
if(nb1 == ((k << 1) | 1)) {
tourn[k] = tourn[k << 1];
k = k << 1;
}
while((k > 1) && (vec[id] <= vec[tourn[(kk = k >> 1)]])) {
tourn[k] = tourn[kk];
k = kk;
}
tourn[k] = tourn[nb];
}
}
/*****************************************************************************
Counts the number of bits equal to 1
*****************************************************************************/
int cardinal(int n) {
int i;
for(i=0; n; n >>= 1)
i += n & 1;
return(i);
}
/**************************************************************************
Conjugate transform of a set function
mu: set function
The result is written in mu_out
***************************************************************************/
void setfunction2conjugate(double *mu, int *n, double *mu_out) {
int i;
int pow = 1<<*n;
for (i=0;i<pow;i++)
mu_out[i] = mu[pow-1] - mu[pow-1-i];
}
/**************************************************************************
Quick computation of the lower cardinality transform
c: factor used in the transformation
mu: set function
The result is written in mu
***************************************************************************/
void fast_lower_cardinality_transform(double *mu, double c, int n) {
int i;
int j,k;
for (i=1;i<=n;i++)
for (j=1;j<(1<<i);j+=2)
for (k=0;k<(1<<(n-i));k++)
mu[j*(1<<(n-i))+k] += c*mu[(j-1)*(1<<(n-i))+k];
}
/****************************************************************************
Computation of the Mobius transform using Robert Kennes algorithm (1990)
The result is written in mu.
***************************************************************************/
void setfunction2Mobius(double *mu, int *n) {
fast_lower_cardinality_transform(mu,-1.0,*n);
}
/****************************************************************************
Computation of the inverse Mobius transform.
Old version: Robert Kennes algorithm (1990)
Pre: v is given in binary order.
void Mobius2setfunction(double *v, int *n)
{
fast_lower_cardinality_transform(v,1.0,*n);
}
Current version: less efficient
Pre: a is given in "natural" order
***************************************************************************/
void Mobius2setfunction(int *n, int *k, double *a, int *subset, double *mu) {
int i,j;
int sb = (int)sum_binom(*n,*k);
for (i=0; i<(1<<*n); i++) {
mu[i] = a[0];
for (j=1; j<sb; j++) {
if ((subset[j] & i) == subset[j])
mu[i] += a[j];
if (subset[j] == i)
break;
}
}
}
/****************************************************************************
Normalizes a Mobius transform ( mu(X) = 1 )
***************************************************************************/
void normalize_Mobius(int n, int k, double *a) {
int i;
int sb = (int)sum_binom(n,k);
double s = 0.0;
for (i=0; i<sb; i++)
s += a[i];
for (i=0; i<sb; i++)
a[i] /= s;
}
/*****************************************************************************
k-truncation of a set function mu.
Returns its truncated Mobius representation.
Result in a (dim sum_binom(n,k)).
*****************************************************************************/
void k_truncation(int *n, int *k, double *mu, int *subset, double *a) {
int i;
setfunction2Mobius(mu, n);
for (i=0;i<sum_binom(*n,*k);i++)
a[i] = mu[subset[i]];
}
/*****************************************************************************
Generates a cardinal set function from an array of size n+1 (csf)
Post: sf given in binary order
*****************************************************************************/
void cardinal2setfunction(int *n, double *csf, double *sf) {
int i;
for (i=0; i<(1<<*n); i++)
sf[i] = csf[cardinal(i)];
}
/*****************************************************************************
Generates an array of size n+1 (csf) from a cardinal set function
Pre: sf given in "natural" order and may be truncated
*****************************************************************************/
void setfunction2cardinal(int *n, int *k, double *sf, double *csf) {
int i;
int sb = 0;
for (i=0; i<=*k; i++) {
csf[i] = sf[sb];
sb += binom(*n,i);
}
for (i=*k+1;i<=*n;i++)
csf[i] = 0.0;
}
/*****************************************************************************
Calculation of the Choquet integral.
mu : game (2^n coefficients)
f : function (n coefficients)
*****************************************************************************/
void Choquet_integral_game(int *n, double *mu, double *f, double *resul) {
int i;
int tourn[NMAX];
int *index = (int *) R_alloc(*n, sizeof(int));
tri(*n, tourn, f, index);
*resul = f[index[0]] * mu[(1<<*n)-1];
for(i=1; i<*n; i++)
*resul += (f[index[i]] - f[index[i - 1]])
* mu[subset2binary(index + i, *n - i)];
index = NULL;
}
/*****************************************************************************
Calculation of the Choquet integral.
a : Mobius representation of a game
f : function (n coefficients)
*****************************************************************************/
void Choquet_integral_Mobius(int *n, int *k, double *a, int *subset,
double *f, double *resul) {
int i,j,l;
int sb = (int)sum_binom(*n,*k);
double min_f;
*resul = 0.0;
for (i=1;i<sb;i++)
{
for (j=0;j<*n;j++)
if (subset[i] & 1<<j)
{
min_f = f[j];
break;
}
for (l=j+1;l<*n;l++)
if (subset[i] & 1<<l)
min_f = INF(min_f, f[l]);
*resul += a[i] * min_f;
}
}
/*****************************************************************************
Calculation of the Sugeno integral.
mu : game (2^n coefficients)
f : function (n coefficients)
*****************************************************************************/
void Sugeno_integral_game(int *n, double *mu, double *f, double *resul) {
int i;
int tourn[NMAX];
int *index = (int *) R_alloc(*n, sizeof(int));
tri(*n, tourn, f, index);
*resul = INF(f[index[0]], mu[(1<<*n) -1]);
for(i=1; i<*n; i++)
*resul = SUP(*resul, INF(f[index[i]], mu[subset2binary(index+i, *n-i)]));
index = NULL;
}
/*****************************************************************************
Calculation of the Sugeno integral.
a : Mobius representation of a game
f : function (n coefficients)
VERY VERY VERY SUBOPTIMAL!!!
*****************************************************************************/
void Sugeno_integral_Mobius(int *n, int *k, double *a, int *subset,
double *f, double *resul) {
int i,j,s;
int sb = (int)sum_binom(*n,*k);
int tourn[NMAX];
int *index = (int *) R_alloc(*n, sizeof(int));
double mu;
tri(*n, tourn, f, index);
mu = 0.0;
for (i=1;i<sb;i++)
mu += a[i];
*resul = INF(f[index[0]], mu);
for(i=1; i<*n; i++) {
s = subset2binary(index+i, *n-i);
mu = 0.0;
for (j=1; j < sb; j++)
if ((subset[j] & s) == subset[j])
mu += a[j];
*resul = SUP(*resul, INF(f[index[i]], mu));
}
index = NULL;
}
/*****************************************************************************
Calculation of the Sipos integral.
mu : game (2^n coefficients)
f : function (n coefficients)
*****************************************************************************/
void Sipos_integral_game(int *n, double *mu, double *f, double *resul) {
int i, p;
int tourn[NMAX];
int *index = (int *) R_alloc(*n, sizeof(int));
/* sorts f */
tri(*n, tourn, f, index);
/* Computes p such that
f[index[0]] <= ... <= f[index[p-1]] < 0 <= f[index[p]] <= ...
<= f[index[n-1]]
*/
for(p=0; (p<*n) && (f[index[p]]<0); p++);
*resul=0.0;
/* Computes the terms of the integral for i<=p */
if(p>0) {
for(i=0; i<p-1; i++)
*resul += (f[index[i]] - f[index[i+1]])
* mu[subset2binary(index, i+1)];
*resul += f[index[p-1]] * mu[subset2binary(index, p)];
}
/* Computes the terms of the integral for i>p */
if(p<*n) {
*resul += f[index[p]] * mu[subset2binary(index+p, *n-p)];
for(i=p+1; i<*n; i++)
*resul += (f[index[i]] - f[index[i-1]])
* mu[subset2binary(index+i, *n-i)];
}
index = NULL;
}
/*****************************************************************************
Generation of the first k + 1 levels of the power set of X
in the "natural" order. Recursive function.
*****************************************************************************/
void k_power_set_rec(int n, int k, int last, int *power_set, int *b) {
int i, istart;
/* look for the leftmost 1 in b and start to fill blank cases
with 1 left from this position */
istart = n;
if(*b > 0)
while(!(*b & 1<<(istart-1)))
istart--;
else
istart = 0;
for(i=istart; i<n; i++) {
last++;
*(power_set+last) = *b + (1<<i);
}
if(last != (int)sum_binom(n,k) - 1)
k_power_set_rec(n, k, last, power_set, b+1);
}
/*****************************************************************************
Generation of the first k + 1 levels of the power set of X
in the "natural" order. Wrapps the previous function.
*****************************************************************************/
void k_power_set(int *n, int *k, int *power_set) {
power_set[0] = 0;
k_power_set_rec(*n, *k, 0, power_set, power_set);
}
/*****************************************************************************
Converts the k power set of X in the "natural" order to char**.
*****************************************************************************/
void k_power_set_char(int *n, int *k, int *k_power_set, char **subset) {
int i, j;
int x[NMAX];
char string[255];
sprintf(subset[0],"{}");
for(i=1; i<sum_binom(*n,*k); i++) {
for(j=0; j<*n; j++)
x[j]=0;
binary2subset(*n,k_power_set[i],x);
sprintf(subset[i],"{%d",x[0]+1);
for(j=1; j<cardinal(k_power_set[i]); j++) {
sprintf(string,",%d", x[j]+1);
strcat(subset[i],string);
}
strcat(subset[i],"}");
}
}
/*****************************************************************************
Converts the power set of X in the binary order to char**.
*****************************************************************************/
void power_set_binary_char(int *n, char **power_set) {
int i, j;
int x[NMAX];
char string[255];
sprintf(power_set[0],"{}");
for(i=1; i<(1<<*n); i++) {
for(j=0; j<*n; j++)
x[j]=0;
binary2subset(*n,i,x);
sprintf(power_set[i],"{%d",x[0]+1);
for(j=1; j<cardinal(i); j++) {
sprintf(string,",%d", x[j]+1);
strcat(power_set[i],string);
}
strcat(power_set[i],"}");
}
}
/*****************************************************************************
Prints the set function in "natural" order in the R terminal
*****************************************************************************/
void Rprint_setfunction(int *n, int *k, double *mu, int *subset, int *mobius)
{
int i,j;
int x[NMAX];
Rprintf("{}\t\t%lf\n",mu[0]);
for(i=1; i<sum_binom(*n,*k); i++)
{
for(j=0; j<*n; j++)
x[j]=0;
binary2subset(*n,subset[i],x);
Rprintf("{%d",x[0]+1);
for(j=1; j<cardinal(subset[i]); j++)
Rprintf(",%d",x[j]+1);
if (*mobius)
Rprintf("}\t\t%lf\n",mu[i]);
else
Rprintf("}\t\t%lf\n",mu[subset[i]]);
}
}
/*****************************************************************************
Writing a set function given in binary order in the "natural" order
Pre: power_set contains the power_set in "natural" order
*****************************************************************************/
void binary2natural(int *n, double *sf, int *power_set, double *sf_out) {
int i;
for(i=0; i<(1<<*n); i++)
sf_out[i] = sf[power_set[i]];
}
/*****************************************************************************
Writing a set function given in "natural" order in the binary order
Pre: power_set contains the power_set in "natural" order
*****************************************************************************/
void natural2binary(int *n, double *sf, int *power_set, double *sf_out) {
int i;
for(i=0; i<(1<<*n); i++)
sf_out[power_set[i]] = sf[i];
}
/*****************************************************************************
Is sf a k-cardinal set function? Used to test both the cardinality of
a set function or of its Mobius representation.
Pre: sf is given in "natural" order.
*****************************************************************************/
void is_kcardinal(int *n, int *k, double *sf, int *flag) {
int i, j, l;
int cni;
*flag = 0;
l = 1;
for(i=1; i<=INF(*k,*n-1); i++) {
cni = (int)binom(*n,i);
for (j=1; j < cni; j++) {
if (sf[l] != sf[l+1]) {
*flag = 1;
return;
}
l++;
}
l++;
}
}
/*****************************************************************************
Is mu a cardinal set function?
Pre: mu given in binary order
*****************************************************************************/
void is_cardinal_setfunction(int *n, double *mu, int *power_set, int *flag) {
double *mu_out = (double *) R_alloc((1<<*n), sizeof(double));
binary2natural(n, mu, power_set, mu_out);
is_kcardinal(n, n, mu_out, flag);
mu = NULL;
}
/*****************************************************************************
Is mu a k-additive set function?
Pre: 1 <= k <= n
kmax: the order of truncation of mu.
a: Mobius representation of mu, given in the "natural" order
*****************************************************************************/
void is_kadditive_Mobius(int *n, int *kmax, int *k, double *a, double *epsilon,
int *flag) {
int i;
int sb = (int)sum_binom(*n,*k-1);
int cnk = (int)binom(*n,*k);
*flag = 1;
for (i=0; i<cnk; i++)
if (fabs(a[sb + i]) > *epsilon) {
*flag = 0;
break;
}
if (*flag == 1)
return;
sb += cnk;
for (i=sb; i<(int)sum_binom(*n,*kmax); i++)
if (fabs(a[i]) > *epsilon) {
*flag = 1;
return;
}
}
/*****************************************************************************
Is mu a k-additive set function?
Pre: mu given in binary order
*****************************************************************************/
void is_kadditive_setfunction(int *n, int *k, double *mu, int *power_set,
double *epsilon, int *flag) {
double *mu_out = (double *) R_alloc((1<<*n), sizeof(double));
setfunction2Mobius(mu,n);
binary2natural(n, mu, power_set, mu_out);
is_kadditive_Mobius(n, n, k, mu_out, epsilon, flag);
mu = NULL;
}
/*****************************************************************************
Adds a strict veto criterion to mu_init (n criteria) at position n.
Result in mu (n+1 criteria).
*****************************************************************************/
void add_veto_setfunction(int *n, double *mu_init, double *mu) {
int i;
int power2n = 1<<*n;
for(i=0; i<power2n; i++)
mu[i+power2n] = mu_init[i];
}
/*****************************************************************************
Searchs for the upper neighbors
*****************************************************************************/
void search_upper_neighbors(int n, int node, int *neighbors) {
int i, count = 0;
for(i=0; i<n; i++)
if(!(1<<i & node)) {
neighbors[count] = node + (1<<i);
count++;
}
}
/*****************************************************************************
Searchs for the lower neighbors
*****************************************************************************/
void search_lower_neighbors(int n, int node, int *neighbors) {
int i, count = 0;
for(i=0; i<n; i++)
if(1<<i & node) {
neighbors[count] = node - (1<<i);
count++;
}
}
/*****************************************************************************
Is the set function mu monotone?
Pre: mu is given in binary order
*****************************************************************************/
void is_monotone_setfunction(int *n, double *mu, int *print, double *epsilon,
int *flag) {
int i, j, k, level, c;
int neighbors[NMAX];
int subset[NMAX];
*flag = 0;
for(i=0; i<(1<<*n)-1; i++) {
search_upper_neighbors(*n, i, neighbors);
level = cardinal(i);
for(j=0; j<*n-level; j++) {
if(mu[ i ] - mu[neighbors[j]] > *epsilon) {
*flag = 1;
if(*print) {
Rprintf("Violation of monotonicity constraint between {");
binary2subset(*n, i, subset);
c = cardinal(i);
for(k=0; k<cardinal(i); k++)
Rprintf(" %d", subset[k] + 1);
Rprintf(" } and {");
binary2subset(*n, neighbors[j], subset);
for(k=0; k<c+1; k++)
Rprintf(" %d", subset[k] + 1);
Rprintf(" }. \n");
}
else
return;
}
}
}
}
/*****************************************************************************
Does the Mobius representation a correspond to a monotone set function?
Pre: a given in "natural" order
*****************************************************************************/
void is_monotone_Mobius(int *n, int *k, double *a, int *subset, int *print,
double *epsilon, int *flag) {
int i,j,l,m,c,r;
int pow = 1<<*n;
int sb = (int)sum_binom(*n,*k);
double s;
int list[NMAX];
*flag = 0;
for (i=0; i<*n; i++)
for (j=1; j<pow; j++)
if(j & 1<<i) {
s = 0.0;
for(l=1;l<sb;l++)
if (((subset[l] & j) == subset[l]) && (subset[l] & 1<<i))
s += a[l];
if (s < - *epsilon) {
*flag = 1;
if (*print) {
r = j ^ 1<<i;
Rprintf("Violation of monotonicity constraint between {");
binary2subset(*n, r, list);
c = cardinal(r);
for(m=0; m<c; m++)
Rprintf(" %d", list[m] + 1);
Rprintf(" } and {");
binary2subset(*n, j, list);
for(m=0; m<c+1; m++)
Rprintf(" %d", list[m] + 1);
Rprintf(" }. \n");
}
else
return;
}
}
}
/*****************************************************************************
Factorial n
*****************************************************************************/
double fact(int n) {
int i;
double f = 1;
for(i=n; i>0; i--)
f *= i;
return(f);
}
/*****************************************************************************
The gamma function
*****************************************************************************/
double gamm(int a, int n) {
return(fact(n - a - 1) * fact(a) / fact(n));
}
/*****************************************************************************
The zeta function
*****************************************************************************/
double zeta(int a, int n) {
return(fact(n - a - 2) * fact(a) / fact(n - 1));
}
/*****************************************************************************
Choose(n,k)
*****************************************************************************/
double binom(int n, int k) {
/* return(fact(n) / (fact(k) * fact(n-k))); */
if ((k==0) || (k==n))
return 1.0;
else
return binom(n-1,k) + binom(n-1,k-1);
}
/*****************************************************************************
sum_{i=0}^k Choose(n,i)
*****************************************************************************/
double sum_binom(int n, int k) {
int i;
double s = 1.0;
for (i=1;i<=k;i++)
s += binom(n,i);
return s;
}
/*****************************************************************************
Calculation of the Shapley value of mu.
Result in phi (dimension n).
*****************************************************************************/
void Shapley_value_setfunction(int *n, double *mu, double *phi) {
int i, k;
int pow = 1<<*n;
for(k=0; k<*n; k++) {
phi[k] = 0.0;
for(i=0; i<pow; i++)
if(!(i & 1<<k))
phi[k] += gamm(cardinal(i), *n) * (mu[i + (1<<k)] - mu[i]);
}
}
/*****************************************************************************
Calculation of the Shapley value of a set function using its
Mobius representation a.
Result in phi (dimension n).
*****************************************************************************/
void Shapley_value_Mobius(int *n, int *k, double *a, int *subset,
double *phi) {
int i,j;
int sb = (int)sum_binom(*n,*k);
for (i=0;i<*n;i++) {
phi[i] = 0.0;
for (j=1; j<sb; j++)
if (subset[j] & 1<<i)
phi[i] += a[j] / (double)cardinal(subset[j]);
}
}
/*****************************************************************************
Calculation of the Shapley interaction indices of mu.
Result in phi (dimension n*n).
*****************************************************************************/
void interaction_indices_setfunction(int *n, double *mu, double *phi) {
int i, j, k;
int pow = 1<<*n;
for(k=0; k<*n; k++) {
for(j=k+1; j<*n; j++) {
phi[k*(*n)+j] = 0.0;
for(i=0; i<pow; i++)
if(!(i & 1<<k) && !(i & 1<<j))
phi[k*(*n)+j] += zeta(cardinal(i), *n) *
(mu[i+(1<<k)+(1<<j)] - mu[i+(1<<k)] - mu[i+(1<<j)]+ mu[i]);
}
}
for(k=0; k<*n; k++)
for(j=k+1; j<*n;j++)
phi[j*(*n)+k] = phi[k*(*n)+j];
for(k=0; k<*n; k++)
phi[k*(*n)+k] = 0.0;
}
/*****************************************************************************
Calculation of the Shapley interaction indices of a set function using its
Mobius representation a.
Result in phi (dimension n*n).
*****************************************************************************/
void interaction_indices_Mobius(int *n, int *k, double *a, int *subset,
double *phi) {
int i, j, l;
int sb = (int)sum_binom(*n,*k);
for(i=0; i<*n; i++)
for(j=i+1; j<*n; j++) {
phi[i*(*n)+j] = 0.0;
for(l=1; l<sb; l++)
if((subset[l] & 1<<i) && (subset[l] & 1<<j))
phi[i*(*n)+j] += a[l]/(double)(cardinal(subset[l])-1);
}
for(i=0; i<*n; i++)
for(j=i+1; j<*n;j++)
phi[j*(*n)+i] = phi[i*(*n)+j];
for(i=0; i<*n; i++)
phi[i*(*n)+i] = 0.0;
}
/*****************************************************************************
Calculation of veto indices of a capacity mu
Result in v (dimension n)
*****************************************************************************/
void veto_capacity(int *n, double *mu, double *v) {
int i, k;
int pow = 1<<*n;
for(k=0; k<*n; k++) {
v[k] = 0.0;
for(i=1; i<pow; i++)
if(!(i & 1<<k))
v[k] += mu[i] / binom(*n-1,cardinal(i));
v[k] /= (double)(*n - 1) * mu[pow-1];
v[k] = 1.0 - v[k];
}
}
/*****************************************************************************
Calculation of veto indices from the Mobius representation of a capacity
Result in v (dimension n)
*****************************************************************************/
void veto_Mobius(int *n, int *k, double *a, int *subset, double *v) {
int i,j;
int sb = (int)sum_binom(*n,*k);
normalize_Mobius(*n,*k,a);
for(i=0; i<*n; i++) {
v[i] = 0.0;
for(j=1; j<sb; j++)
if(!(subset[j] & 1<<i))
v[i] += a[j] / (double)(cardinal(subset[j]) + 1);
v[i] *= (double)(*n)/(double)(*n - 1);
v[i] = 1.0 - v[i];
}
}
/*****************************************************************************
Calculation of favor indices of a capacity mu
Result in favor (dimension n)
*****************************************************************************/
void favor_capacity(int *n, double *mu, double *f) {
int i, k;
int pow = 1<<*n;
for(k=0; k<*n; k++) {
f[k] = 0.0;
for(i=0; i<pow; i++)
if(!(i & 1<<k))
f[k] += mu[i + (1<<k)] / binom(*n-1,cardinal(i));
f[k] /= (double)(*n - 1) * mu[pow-1];
f[k] -= 1.0/(double)(*n - 1);
}
}
/*****************************************************************************
Calculation of favor indices from the Mobius representation of a capacity
Result in f (dimension n)
*****************************************************************************/
void favor_Mobius(int *n, int *k, double *a, int *subset, double *f) {
int i,j,l;
int sb = (int)sum_binom(*n,*k);
normalize_Mobius(*n,*k,a);
for(i=0; i<*n; i++) {
f[i] = 0.0;
for(j=0; j<sb; j++)
if(!(subset[j] & 1<<i)) {
for (l=j+1;l<sb;l++)
if ((subset[j] | 1<<i) == subset[l])
break;
if (l < sb)
f[i] += (a[j] + a[l])
/ (double)(cardinal(subset[j]) + 1);
else
f[i] += a[j] / (double)(cardinal(subset[j]) + 1);
}
f[i] *= (double)(*n)/(double)(*n - 1);
f[i] -= 1.0/(double)(*n - 1);
}
}
/*****************************************************************************
Calculation of the orness of a capacity mu
*****************************************************************************/
void orness_capacity(int *n, double *mu, double *resul) {
int i;
int pow = 1<<*n;
*resul = 0.0;
for(i=1; i<pow-1; i++) /* could be optimized */
*resul += mu[i] / binom(*n, cardinal(i));
*resul /= (double)(*n - 1) * mu[pow-1];
}
/*****************************************************************************
Calculation of the orness from the Mobius representation of a capacity
*****************************************************************************/
void orness_Mobius(int *n, int *k, double *a, int *subset, double *resul) {
int i,c;
int sb = (int)sum_binom(*n,*k);
normalize_Mobius(*n,*k,a);
*resul = 0.0;
for(i=1; i<sb; i++) {
c = cardinal(subset[i]);
*resul += (double)(*n - c) * a[i] / (double)(c + 1);
}
*resul /= (double)(*n - 1);
}
/*****************************************************************************
Calculation of the variance of a capacity mu
*****************************************************************************/
void variance_capacity(int *n, double *mu, double *resul) {
int i, k;
int pow = 1<<*n;
*resul = 0.0;
for(k=0; k<*n; k++)
for(i=0; i<pow; i++)
if(!(i & 1<<k))
*resul += gamm(cardinal(i), *n)
* SQR((mu[i + (1<<k)] - mu[i])
/ mu[pow-1]);
*resul -= 1.0/(double)*n;
*resul /= 1 - 1.0/(double)*n;
}
/*****************************************************************************
Calculation of the variance from the Mobius representation of a capacity
*****************************************************************************/
void variance_Mobius(int *n, int *k, double *a, int *subset, double *resul) {
int i,j,l;
int pow = 1<<*n;
int sb = (int)sum_binom(*n,*k);
double s;
normalize_Mobius(*n,*k,a);
for (i=0; i<*n; i++)
for (j=1; j<pow; j++)
if(j & 1<<i) {
s = 0.0;
for(l=1;l<sb;l++)
if (((subset[l] & j) == subset[l])
&& (subset[l] & 1<<i))
s += a[l];
*resul += gamm(cardinal(j)-1, *n) * SQR(s);
}
*resul -= 1.0/(double)*n;
*resul /= 1 - 1.0/(double)*n;
}
/*****************************************************************************
Calculation of the normalized Marichal entropy of a capacity mu
*****************************************************************************/
void entropy_capacity(int *n, double *mu, double *resul) {
int i, k;
int pow = 1<<*n;
*resul = 0.0;
for(k=0; k<*n; k++)
for(i=0; i<pow; i++)
if(!(i & 1<<k))
*resul += gamm(cardinal(i), *n)
* H((mu[i + (1<<k)] - mu[i])
/mu[pow-1]);
*resul /= log(*n);
}
/*****************************************************************************
Calculation of the normalized Marichal
from the Mobius representation of a capacity
*****************************************************************************/
void entropy_Mobius(int *n, int *k, double *a, int *subset, double *resul) {
int i,j,l;
int pow = 1<<*n;
int sb = (int)sum_binom(*n,*k);
double s;
normalize_Mobius(*n,*k,a);
for (i=0; i<*n; i++)
for (j=1; j<pow; j++)
if(j & 1<<i) {
s = 0.0;
for(l=1;l<sb;l++)
if (((subset[l] & j) == subset[l])
&& (subset[l] & 1<<i))
s += a[l];
*resul += gamm(cardinal(j)-1, *n) * H(s);
}
*resul /= log(*n);
}
/****************************************************************************/