00001 #ifndef __enumerate
00002 #define __enumerate
00003
00004 #include<algorithm>
00005 #include<cmath>
00006 #include<vector>
00007
00008 #include"weight.hpp"
00009
00010 #ifdef DEBUG_VERBOSE
00011 #include<iostream>
00012 #endif
00013
00014 namespace graphgen{
00015
00016
00017 class EnumA{
00018 int Nv;
00019 int Nl;
00020
00021 int digits;
00022 int b;
00023 bool valid;
00024
00025
00026 std::vector<int> A;
00027
00028
00029 int aij( int i, int j ) const {
00030
00031 if(i==j) return 2*A[i*(i+3)/2];
00032
00033 return A[(std::max(i,j)*(std::max(i,j)+1))/2 + std::min(i,j)];
00034 }
00035
00036 public:
00037
00038 EnumA(int NV, int NL, int base): Nv(NV), Nl(NL), digits(Nv*(Nv+1)/2), b(base), A(digits,0){
00039
00040 int n=Nl;
00041 for( int i=0; (i<digits) && (n>0); i++ ) n -= (A[i]=std::min(b-1,n));
00042 valid = (n==0);
00043 }
00044
00045 int operator()( int i, int j ) const { return aij(i,j); }
00046
00047
00048 EnumA& operator++(){
00049
00050 int sum = Nl;
00051
00052 int carry;
00053 do{
00054 carry = 1;
00055
00056 for( int i=0; (carry>0)&&(i<digits); i++ ){
00057 if((A[i]+=carry)>=b){
00058 sum-=(A[i]-1);
00059 A[i]=0;
00060 carry=1;
00061 } else {
00062 carry = 0;
00063 sum++;
00064 }
00065 }
00066 }while( sum != Nl && carry==0 );
00067
00068 valid = (carry == 0);
00069 return *this;
00070 }
00071
00072 operator bool() const { return valid; }
00073
00074 int NV() const {return Nv;}
00075 int NL() const {return Nl;}
00076
00077 int Q( int i ) const {
00078 int Qi=0;
00079 for( int j=0; j<Nv; j++ ) Qi+=aij(i,j);
00080 return Qi;
00081 }
00082
00083 void printdigits( std::ostream& os ) const {
00084 for( int i=0; i<digits; i++ ) os<<A[i];
00085 os << std::endl;
00086 }
00087
00088 void printmatrix( std::ostream& os ) const{
00089 for( int i=0; i<Nv; i++ ){
00090 for( int j=0; j<Nv; j++ ) os << aij(i,j) <<"\t";
00091 os<<std::endl;
00092 }
00093 }
00094 };
00095
00096
00097
00098
00099
00100
00101 namespace conservation{
00102
00103
00104
00105 struct microcanonic{
00106 static bool is_microcan(){ return true; }
00107 static bool is_canonic(){ return false; }
00108 static bool is_grandcan(){ return false; }
00109
00110 bool operator()( const EnumA& A, const std::vector<int>& Q ){
00111
00112
00113
00114 for( int i=0; i<A.NV(); i++ ) if(A.Q(i)!=Q[i]) return false;
00115 return true;
00116 }
00117 };
00118
00119 struct canonic{
00120 static bool is_microcan(){ return false; }
00121 static bool is_canonic(){ return true; }
00122 static bool is_grandcan(){ return false; }
00123
00124 bool operator()( const EnumA& ){
00125
00126
00127 return true;
00128 }
00129 };
00130
00131 struct grandcanonic{
00132 static bool is_microcan(){ return false; }
00133 static bool is_canonic(){ return false; }
00134 static bool is_grandcan(){ return true; }
00135
00136 bool operator()( const EnumA& ){
00137 return true;
00138 }
00139 };
00140
00141
00142 struct graph{
00143 static const bool allow_multi = true;
00144 static const bool allow_self = true;
00145
00146 bool operator()( const EnumA& ){
00147 return true;
00148 }
00149 };
00150
00151 struct multi{
00152 static const bool allow_multi = true;
00153 static const bool allow_self = false;
00154
00155 bool operator()( const EnumA& A ){
00156 for( int i=0; i<A.NV(); i++ ) if(A(i,i)!=0) return false;
00157 return true;
00158 }
00159 };
00160
00161 struct simple{
00162 static const bool allow_multi = false;
00163 static const bool allow_self = false;
00164
00165 bool operator()( const EnumA& A ){
00166 for( int i=0; i<A.NV(); i++ ) for( int j=0; j<i; j++ ) if(A(i,j)>1) return false;
00167 for( int i=0; i<A.NV(); i++ ) if(A(i,i)!=0) return false;
00168 return true;
00169 }
00170 };
00171
00172
00173 struct tree{
00174 };
00175
00176 }
00177
00178
00179
00180
00181
00182 template<value_type::weight_value_type WV, template<value_type::weight_value_type WV>class W> inline
00183 double graph_weight( const EnumA& A, const vertex_weight<WV>& w );
00184
00185 template<value_type::weight_value_type WV> inline
00186 double graph_weight( const EnumA& A, const vertex_weight<WV>& w ){
00187 if(w.logarithmic){
00188 double weight = 0;
00189 for( int i=0; i<A.NV(); i++ ) weight += w.p( A.Q(i) );
00190
00191 return weight;
00192 } else {
00193 double weight=1;
00194 for( int i=0; i<A.NV(); i++ ) weight *= w.p( A.Q(i) );
00195 return weight;
00196 }
00197 }
00198
00199 template<value_type::weight_value_type WV> inline
00200 double graph_weight( const EnumA& A, const link_weight<WV>& w ){
00201 if(w.logarithmic()){
00202 double weight = 0;
00203 for( int i=0; i<A.NV(); i++ ) for( int j=0; j<i; j++ )
00204 weight += A(i,j)*w.weight( A.Q(i), A.Q(j));
00205 return weight;
00206 } else {
00207 double weight=1;
00208 for( int i=0; i<A.NV(); i++ ) for( int j=0; j<i; j++ ) if(A(i,j))
00209 weight *= pow(w.weight( A.Q(i), A.Q(j) ), A(i,j));
00210
00211 return weight;
00212 }
00213 }
00214
00215
00216 template<class R>
00217 class enumerate{
00218 int nv;
00219 int nl;
00220
00221 double addln(double l1, double l2)
00222 {
00223 return std::max(l1,l2)+log1p(exp(-abs(l2-l1)));
00224 }
00225
00226 public:
00227 enumerate( int NV, int NL ): nv(NV),nl(NL) {}
00228
00229 template<class E, class W>
00230 void eval( const W&w, std::vector<double>° )
00231 {
00232 E ensemble_conserved;
00233 R restrict_conserved;
00234
00235 double SumW = 0;
00236
00237 if(E::is_canonic()){
00238
00239 int base=2;
00240 if(R::allow_multi) base=nl+1;
00241
00242 for( EnumA A(nv,nl,base); A; ++A ) {
00243
00244 if( ensemble_conserved(A) && restrict_conserved(A) ){
00245 double WA=graph_weight(A,w);
00246 for( int i=0; i<nv; i++ ) deg[A.Q(i)] += WA;
00247 SumW += WA;
00248 #ifdef DEBUG_VERBOSE
00249 A.printmatrix(std::cerr);
00250 std::cerr<<"(W="<<WA<<std::endl<<std::endl;
00251 #endif
00252 }
00253 }
00254 } else if(E::is_microcan()) {
00255
00256
00257
00258 } else if(E::is_grandcan()) {
00259
00260
00261 }
00262
00263 for( unsigned int q=0; q<deg.size(); q++ ) deg[q] /= (nv*SumW);
00264 }
00265 };
00266
00267 }
00268
00269 #endif