Merge branch 'no-edges2'
[libdai.git] / src / x2x.cpp
1 /*
2 Copyright (C) 2005 Martijn Leisink m.leisink at science.ru.nl
3
4 This program is free software; you can redistribute it and/or
5 modify it under the terms of the GNU General Public License
6 as published by the Free Software Foundation; either version 2
7 of the License, or (at your option) any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17
18 CHANGES:
19 2006-11-20 Joris Mooij
20 * removed MATLAB interface
21 * put into namespace "x2x"
22 */
23
24
25 #include <cmath>
26 #include <cstring>
27
28
29
30 namespace x2x {
31
32 // helper functions to compute the sum over all partitions
33 double psum (double *x, long s, int n=0);
34 double psumx (double *x, long a, long s, int n) {
35 // recursively process all one bits remaining after psum
36 // put some of them into the same subset and call psum again for the rest
37 if((s>>n)>1) {
38 while(!((s>>++n)&1))
39 {};
40 return(psumx(x,a,s,n)+psumx(x,a^(1l<<n),s,n));
41 } else {
42 return(x[a]*psum(x,s^a,0));
43 };
44 }
45 double psum (double *x, long s, int n) {
46 // take the first one bit and put it in the first subset, then call psumx
47 if(s>>n) {
48 while(!((s>>n)&1)) ++n;
49 return(psumx(x,1l<<n,s,n));
50 } else return(1);
51 }
52
53 // convert moments to cumulants upto order k
54 void m2c (int N, double *x, int k) {
55 int *c=new int[k+1];
56 long s;
57 int z;
58 x[0]=log(x[0]); // to get the correct answer if not normalized
59 for(int b=1;b<=k;++b) {
60 // start with marginals, then two-point correlations and so on
61 c[z=b]=-1;
62 s=0; // index into x
63 do {
64 while(--z>=0) s^=(1l<<(c[z]=c[z+1]+1));
65 // c_ijk = m_ijk - c_i*c_jk - c_j*c_ik - c_k*c_ij - c_i*c_j*c_k
66 x[s]=2*x[s]-psum(x,s);
67 // increment b indices
68 for(z=0;z<b&&(s^=3l<<c[z],++c[z]==N-z);++z) s^=1l<<c[z];
69 } while(z<b);
70 };
71 delete[] c;
72 }
73
74 // convert cumulants to moments upto order k
75 void c2m (int N, double *x, int k) {
76 int *c=new int[k+1];
77 long s;
78 int z;
79 for(int b=k;b>=1;--b) {
80 // start with k-point cumulants, then k-1-point cumulants and so on
81 c[z=b]=-1;
82 s=0; // index into x
83 do {
84 while(--z>=0) s^=(1l<<(c[z]=c[z+1]+1));
85 // m_ijk = c_ijk + c_i*c_jk + c_j*c_ik + c_k*c_ij + c_i*c_j*c_k
86 x[s]=psum(x,s);
87 // increment b indices
88 for(z=0;z<b&&(s^=3l<<c[z],++c[z]==N-z);++z) s^=1l<<c[z];
89 } while(z<b);
90 };
91 x[0]=exp(x[0]); // to get the correct answer if not normalized
92 delete[] c;
93 }
94
95 // convert (generalized) weights to log probability or energy
96 void w2logp (int N, double *x) {
97 for(long s=1l<<N;s>>=1;) for(long j=1l<<N;(j-=(s<<1))>=0;)
98 for(long k=j+s;--k>=j;) { x[k+s]+=x[k]; x[k]=2*x[k]-x[k+s]; };
99 }
100
101 // convert log probability or energy to (generalized) weights
102 void logp2w (int N, double *x) {
103 for(long s=1l<<N;s>>=1;) for(long j=1l<<N;(j-=(s<<1))>=0;)
104 for(long k=j+s;--k>=j;) { x[k]=(x[k]+x[k+s])/2; x[k+s]-=x[k]; };
105 }
106
107 // convert probability to moments
108 void p2m (int N, double *x) {
109 for(long s=1l<<N;s>>=1;) for(long j=1l<<N;(j-=(s<<1))>=0;)
110 for(long k=j+s;--k>=j;) { x[k+s]-=x[k]; x[k]=2*x[k]+x[k+s]; };
111 }
112
113 // convert moments to probability
114 void m2p (int N, double *x) {
115 for(long s=1l<<N;s>>=1;) for(long j=1l<<N;(j-=(s<<1))>=0;)
116 for(long k=j+s;--k>=j;) { x[k]=(x[k]-x[k+s])/2; x[k+s]+=x[k]; };
117 }
118
119 // convert log probability to probability
120 void logp2p (int N, double *x) {
121 for(long s=1l<<N;s--;) x[s]=exp(x[s]);
122 }
123
124 // convert probability to log probability
125 void p2logp (int N, double *x) {
126 for(long s=1l<<N;s--;) x[s]=log(x[s]);
127 }
128
129 // normalize a log probability table
130 void logpnorm (int N, double *x) {
131 double f=x[0];
132 for(long s=1l<<N;--s;)
133 if(f>x[s]) f+=log1p(exp(x[s]-f));
134 else f=x[s]+log1p(exp(f-x[s]));
135 for(long s=1l<<N;s--;) x[s]-=f;
136 }
137
138 // normalize a probability table, use logpnorm whenever possible
139 void pnorm (int N, double *x) {
140 double z=0;
141 for(long s=1l<<N;s--;) z+=x[s];
142 for(long s=1l<<N;s--;) x[s]/=z;
143 }
144
145 // fills table with v for all entries with more than k indices
146 // used for example when cumulants or moments are converted upto some order
147 void fill (int N, double *x, int k, double v) {
148 if(k>N) return;
149 long ss=0,s;
150 int n=0;
151 for(long i=1l<<N;i--;) { // make use of gray code to count number of bits
152 if(n>k) x[ss]=v;
153 s=i^(i>>1);
154 if(s>ss) ++n; else --n;
155 ss=s;
156 };
157 }
158
159 } // end of namespace x2x