--- /dev/null
+ GNU GENERAL PUBLIC LICENSE
+ Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.
+ 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The licenses for most software are designed to take away your
+freedom to share and change it. By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users. This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it. (Some other Free Software Foundation software is covered by
+the GNU Library General Public License instead.) You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+ To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have. You must make sure that they, too, receive or can get the
+source code. And you must show them these terms so they know their
+rights.
+
+ We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+ Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software. If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+ Finally, any free program is threatened constantly by software
+patents. We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary. To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+\f
+ GNU GENERAL PUBLIC LICENSE
+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+ 0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License. The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language. (Hereinafter, translation is included without limitation in
+the term "modification".) Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope. The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+ 1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+ 2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+ a) You must cause the modified files to carry prominent notices
+ stating that you changed the files and the date of any change.
+
+ b) You must cause any work that you distribute or publish, that in
+ whole or in part contains or is derived from the Program or any
+ part thereof, to be licensed as a whole at no charge to all third
+ parties under the terms of this License.
+
+ c) If the modified program normally reads commands interactively
+ when run, you must cause it, when started running for such
+ interactive use in the most ordinary way, to print or display an
+ announcement including an appropriate copyright notice and a
+ notice that there is no warranty (or else, saying that you provide
+ a warranty) and that users may redistribute the program under
+ these conditions, and telling the user how to view a copy of this
+ License. (Exception: if the Program itself is interactive but
+ does not normally print such an announcement, your work based on
+ the Program is not required to print an announcement.)
+\f
+These requirements apply to the modified work as a whole. If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works. But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+ 3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+ a) Accompany it with the complete corresponding machine-readable
+ source code, which must be distributed under the terms of Sections
+ 1 and 2 above on a medium customarily used for software interchange; or,
+
+ b) Accompany it with a written offer, valid for at least three
+ years, to give any third party, for a charge no more than your
+ cost of physically performing source distribution, a complete
+ machine-readable copy of the corresponding source code, to be
+ distributed under the terms of Sections 1 and 2 above on a medium
+ customarily used for software interchange; or,
+
+ c) Accompany it with the information you received as to the offer
+ to distribute corresponding source code. (This alternative is
+ allowed only for noncommercial distribution and only if you
+ received the program in object code or executable form with such
+ an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it. For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable. However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+\f
+ 4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License. Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+ 5. You are not required to accept this License, since you have not
+signed it. However, nothing else grants you permission to modify or
+distribute the Program or its derivative works. These actions are
+prohibited by law if you do not accept this License. Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+ 6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions. You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+ 7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all. For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices. Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+\f
+ 8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded. In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+ 9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number. If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation. If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+ 10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission. For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this. Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+ NO WARRANTY
+
+ 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+ 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+ END OF TERMS AND CONDITIONS
+\f
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+ <one line to give the program's name and a brief idea of what it does.>
+ Copyright (C) <year> <name of author>
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+ Gnomovision version 69, Copyright (C) year name of author
+ Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary. Here is a sample; alter the names:
+
+ Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+ `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+ <signature of Ty Coon>, 1 April 1989
+ Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs. If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library. If this is what you want to do, use the GNU Library General
+Public License instead of this License.
--- /dev/null
+libDAI-0.2.1 (2008-05-26)
+-------------------------
+
+Bugfix release.
+* added missing cstdio header in util.h
+* fixed Properties in MR_CLAMPING_* and MR_EXACT_*
+* added description of the factor graph fileformat
+* improved Makefile
+
+
+libDAI-0.2.0 (2006-11-30)
+-------------------------
+
+First public release.
--- /dev/null
+This file describes the .fg fileformat used in libDAI to store factor graphs.
+
+Markov Random Fields are special cases of factor graphs, as are Bayesian
+networks. A factor graph can be specified as follows: for each factor, one has
+to specify which variables occur in the factor, what their respective
+cardinalities (i.e., number of possible values) are, and a table listing all
+the values of that factor for all possible configurations of these variables.
+A .fg file is not much more than that. It starts with a line containing the
+number of factors in that graph, followed by an empty line. Then all factors
+are specified, one block for each factor, where the blocks are seperated by
+empty lines. Each variable occurring in the factor graph has a unique
+identifier, its index (which should be a nonnegative integer). Comment lines
+start with #.
+
+Each block starts with a line containing the number of variables in that
+factor. The second line contains the indices of these variables, seperated by
+spaces (indices are nonnegative integers and to avoid confusion, it is
+suggested to start counting at 0). The third line contains the number of
+possible values of each of these variables, also seperated by spaces. Note that
+there is some redundancy here, since if a variable appears in more than one
+factor, the cardinality of that variable appears several times in the .fg file.
+The fourth line contains the number of nonzero entries in the factor table.
+The rest of the lines contain these nonzero entries; each entry consists of a
+table index, followed by white-space, followed by the value corresponding to
+that table index. The most difficult part is getting the indexing right. The
+convention that is used is that the left-most variables cycle through their
+values the fastest (similar to MATLAB indexing of multidimensional arrays). An
+example block describing one factor is:
+
+3
+4 8 7
+3 2 2
+11
+0 0.1
+1 3.5
+2 2.8
+3 6.3
+4 8.4
+6 7.4
+7 2.4
+8 8.9
+9 1.3
+10 1.6
+12 6.4
+11 2.6
+
+which corresponds to the following factor:
+
+x_4 x_8 x_7 | value
+ 0 0 0 | 0.1
+ 1 0 0 | 3.5
+ 2 0 0 | 2.8
+ 0 1 0 | 6.3
+ 1 1 0 | 8.4
+ 2 1 0 | 0.0
+ 0 0 1 | 7.4
+ 1 0 1 | 2.4
+ 2 0 1 | 8.9
+ 0 1 1 | 1.3
+ 1 1 1 | 1.6
+ 2 1 1 | 2.6
+
+Note that the value of x_4 changes fastest, followed by that of x_8, and x_7
+varies the slowest, corresponding to the second line of the block ("4 8 7").
+Further, x_4 can take on three values, and x_8 and x_7 each have two possible
+values, as described in the third line of the block ("3 2 2"). The table
+contains 11 non-zero entries (all except for the fifth entry). Note that the
+eleventh and twelveth entries are interchanged.
+
+A final note: the internal representation in libDAI of the factor above is
+different, because the variables are ordered according to their indices
+(i.e., the ordering would be x_4 x_7 x_8) and the values of the table are
+stored accordingly, with the variable having the smallest index changing
+fastest:
+
+x_4 x_7 x_8 | value
+ 0 0 0 | 0.1
+ 1 0 0 | 3.5
+ 2 0 0 | 2.8
+ 0 1 0 | 7.4
+ 1 1 0 | 2.4
+ 2 1 0 | 8.9
+ 0 0 1 | 6.3
+ 1 0 1 | 8.4
+ 2 0 1 | 0.0
+ 0 1 1 | 1.3
+ 1 1 1 | 1.6
+ 2 1 1 | 2.6
--- /dev/null
+# Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+# Radboud University Nijmegen, The Netherlands
+#
+# This file is part of libDAI.
+#
+# libDAI is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# libDAI is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with libDAI; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+
+# We use the BOOST Program Options library
+BOOSTFLAGS = -lboost_program_options
+
+# Compile using GNU C++ Compiler
+CC = g++
+
+# Flags for the C++ compiler
+CCFLAGS = -Wall -fpic -g -DDEBUG -I. -I$(HOME)/include -O3 #-static #-pg #-DVERBOSE
+
+# To enable the Matlab interface, define WITH_MATLAB = yes
+WITH_MATLAB =
+ifdef WITH_MATLAB
+# Replace the following by the directory where Matlab has been installed
+MATLABDIR = /opt/matlab/bin
+MEX = $(MATLABDIR)/mex
+MEXFLAGS = -g -I. -DDEBUG -largeArrayDims #-g means debugging
+endif
+
+# Replace the following with the extension of compiled MEX files on this platform, e.g. .mexglx for x86
+MEXEXT = .mexglx
+
+HEADERS = bipgraph.h diffs.h index.h var.h factor.h varset.h prob.h daialg.h properties.h alldai.h enum.h x2x.h
+
+# target matlabs is disabled by default since it only compiles with a very recent MatLab version
+TARGETS = tests utils libdai.a example testregression
+ifdef WITH_MATLAB
+TARGETS := $(TARGETS) matlabs
+endif
+all : $(TARGETS)
+ echo -e "\a"
+
+matlabs : matlab/dai.$(MEXEXT) matlab/dai_readfg.$(MEXEXT) matlab/dai_writefg.$(MEXEXT) matlab/dai_removeshortloops.$(MEXEXT) matlab/dai_potstrength.$(MEXEXT)
+
+libdai.a : daialg.o alldai.o bp.o clustergraph.o factorgraph.o hak.o jtree.o lc.o mf.o mr.o properties.o regiongraph.o util.o treeep.o weightedgraph.o x2x.o
+ ar rcs libdai.a daialg.o alldai.o bp.o clustergraph.o factorgraph.o hak.o jtree.o lc.o mf.o mr.o properties.o regiongraph.o util.o treeep.o weightedgraph.o x2x.o
+
+tests : tests/test
+
+utils : utils/createfg utils/fg2dot utils/remove_short_loops utils/fginfo
+
+testregression : tests/test
+ echo Testing...this can take a while...
+ cd tests; ./testregression; cd ..
+
+doc : *.h *.cpp doxygen.conf
+ doxygen doxygen.conf
+
+clean :
+ rm *.o *.$(MEXEXT) example matlab/*.$(MEXEXT) matlab/*.o tests/test utils/fg2dot utils/createfg utils/remove_short_loops utils/fginfo libdai.a; echo
+ rm -R doc; echo
+
+
+daialg.o : daialg.cpp $(HEADERS)
+ $(CC) $(CCFLAGS) -c daialg.cpp
+
+bp.o : bp.cpp bp.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c bp.cpp
+
+lc.o : lc.cpp lc.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c lc.cpp
+
+mf.o : mf.cpp mf.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c mf.cpp
+
+factorgraph.o : factorgraph.cpp factorgraph.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c factorgraph.cpp
+
+util.o : util.cpp util.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c util.cpp
+
+regiongraph.o : regiongraph.cpp regiongraph.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c regiongraph.cpp
+
+hak.o : hak.cpp hak.h $(HEADERS) regiongraph.h
+ $(CC) $(CCFLAGS) -c hak.cpp
+
+clustergraph.o : clustergraph.cpp clustergraph.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c clustergraph.cpp
+
+jtree.o : jtree.cpp jtree.h $(HEADERS) weightedgraph.h clustergraph.h regiongraph.h
+ $(CC) $(CCFLAGS) -c jtree.cpp
+
+treeep.o : treeep.cpp treeep.h $(HEADERS) weightedgraph.h clustergraph.h regiongraph.h jtree.h
+ $(CC) $(CCFLAGS) -c treeep.cpp
+
+weightedgraph.o : weightedgraph.cpp weightedgraph.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c weightedgraph.cpp
+
+mr.o : mr.cpp mr.h $(HEADERS)
+ $(CC) $(CCFLAGS) -c mr.cpp
+
+properties.o : properties.cpp $(HEADERS)
+ $(CC) $(CCFLAGS) -c properties.cpp
+
+alldai.o : alldai.cpp $(HEADERS)
+ $(CC) $(CCFLAGS) -c alldai.cpp
+
+x2x.o : x2x.cpp $(HEADERS)
+ $(CC) $(CCFLAGS) -c x2x.cpp
+
+# EXAMPLE
+##########
+
+example : example.cpp $(HEADERS) libdai.a
+ $(CC) $(CCFLAGS) -o example example.cpp -L. libdai.a
+
+# TESTS
+########
+
+tests/test : tests/test.cpp $(HEADERS) libdai.a
+ $(CC) $(CCFLAGS) -o tests/test tests/test.cpp -L. -ldai $(BOOSTFLAGS)
+
+
+# MATLAB INTERFACE
+###################
+
+matlab/dai.$(MEXEXT) : matlab/dai.cpp $(HEADERS) matlab/matlab.o daialg.o alldai.o bp.o clustergraph.o factorgraph.o hak.o jtree.o lc.o mf.o mr.o properties.o regiongraph.o util.o treeep.o weightedgraph.o x2x.o
+ $(MEX) $(MEXFLAGS) -o matlab/dai matlab/dai.cpp matlab/matlab.o daialg.o alldai.o bp.o clustergraph.o factorgraph.o hak.o jtree.o lc.o mf.o mr.o properties.o regiongraph.o util.o treeep.o weightedgraph.o x2x.o
+
+matlab/dai_readfg.$(MEXEXT) : matlab/dai_readfg.cpp $(HEADERS) factorgraph.o matlab/matlab.o
+ $(MEX) $(MEXFLAGS) -o matlab/dai_readfg matlab/dai_readfg.cpp factorgraph.o matlab/matlab.o
+
+matlab/dai_writefg.$(MEXEXT) : matlab/dai_writefg.cpp $(HEADERS) factorgraph.o matlab/matlab.o
+ $(MEX) $(MEXFLAGS) -o matlab/dai_writefg matlab/dai_writefg.cpp factorgraph.o matlab/matlab.o
+
+matlab/dai_removeshortloops.$(MEXEXT) : matlab/dai_removeshortloops.cpp $(HEADERS) factorgraph.o matlab/matlab.o
+ $(MEX) $(MEXFLAGS) -o matlab/dai_removeshortloops matlab/dai_removeshortloops.cpp factorgraph.o matlab/matlab.o
+
+matlab/dai_potstrength.$(MEXEXT) : matlab/dai_potstrength.cpp $(HEADERS) matlab/matlab.o
+ $(MEX) $(MEXFLAGS) -o matlab/dai_potstrength matlab/dai_potstrength.cpp matlab/matlab.o
+
+matlab/matlab.o : matlab/matlab.cpp matlab/matlab.h $(HEADERS)
+ $(MEX) $(MEXFLAGS) -outdir matlab -c matlab/matlab.cpp
+
+
+# UTILS
+########
+
+utils/createfg : utils/createfg.cpp $(HEADERS) factorgraph.o weightedgraph.o util.o
+ $(CC) $(CCFLAGS) -o utils/createfg utils/createfg.cpp factorgraph.o weightedgraph.o util.o $(BOOSTFLAGS)
+
+utils/fg2dot : utils/fg2dot.cpp $(HEADERS) factorgraph.o
+ $(CC) $(CCFLAGS) -o utils/fg2dot utils/fg2dot.cpp factorgraph.o
+
+utils/remove_short_loops : utils/remove_short_loops.cpp $(HEADERS) factorgraph.o
+ $(CC) $(CCFLAGS) -o utils/remove_short_loops utils/remove_short_loops.cpp factorgraph.o
+
+utils/fginfo : utils/fginfo.cpp $(HEADERS) factorgraph.o
+ $(CC) $(CCFLAGS) -o utils/fginfo utils/fginfo.cpp factorgraph.o
-Test. First commit...
-Test continues. Second commit.
-Test continues again. Third commit.
+libDAI - A free/open source C++ library for Discrete Approximate Inference methods
+==================================================================================
+
+v 0.2.1 - May 26, 2008
+
+
+Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+Radboud University Nijmegen, The Netherlands
+
+
+----------------------------------------------------------------------------------
+This file is part of libDAI.
+
+libDAI is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+libDAI is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with libDAI; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+----------------------------------------------------------------------------------
+
+
+
+What is libDAI?
+---------------
+libDAI is a free/open source C++ library (licensed under GPL, see the file
+COPYING for more details) that provides implementations of various
+(deterministic) approximate inference methods for discrete graphical models.
+libDAI supports arbitrary factor graphs with discrete variables (this includes
+discrete Markov Random Fields and Bayesian Networks).
+
+libDAI is not intended to be a complete package for approximate inference.
+Instead, it should be considered as an "inference engine", providing various
+inference methods. In particular, it contains no GUI, currently only supports
+its own file format for input and output (although support for standard file
+formats may be added), and provides no visualization.
+
+Because libDAI is implemented in C++, it is very fast compared with e.g. MatLab
+implementations. libDAI does provide a MatLab interface for easy integration
+with MatLab. Currently, libDAI supports the following deterministic approximate
+inference methods:
+
+ * Mean Field
+ * (Loopy) Belief Propagation
+ * Tree Expectation Propagation
+ * Generalized Belief Propagation
+ * Double-loop GBP
+ * Loop Corrected Approximate Inference
+
+Exact inference by JunctionTree is also provided.
+
+Many of these algorithms are not yet available in similar open source software,
+to the best of the author's knowledge (open source packages supporting both
+directed and undirected graphical models are Murphy's BNT, Intel's PNL and gR).
+
+The library is targeted at researchers; to be able to use the library, a good
+understanding of graphical models is needed. However, the code will hopefully
+find its way into real-world applications as well.
+
+
+Rationale
+---------
+In my opinion, the lack of open source reference implementations hampers
+progress in research on approximate inference. Methods differ widely in terms
+of quality and performance characteristics, which also depend in different ways
+on various properties of the graphical models. Finding the best approximate
+inference method for a particular application therefore often requires
+empirical comparisons. However, implementing and debugging these methods takes
+a lot of time which could otherwise be spent on research. I hope that this code
+will aid researchers to be able to easily compare various (existing as well as
+new) approximate inference methods, in this way accelerating research and
+stimulating real-world applications of approximate inference.
+
+
+Releases
+--------
+Releases can be obtained from http://www.mbfys.ru.nl/~jorism/libDAI.
+License: GNU Public License v2 (or higher).
+
+libDAI-0.2 December 1, 2006
+libDAI-0.2.1 May 26, 2008
+
+
+Acknowledgments
+---------------
+The development reported here is part of the Interactive Collaborative
+Information Systems (ICIS) project, supported by the Dutch Ministry of Economic
+Affairs, grant BSIK03024. I would like to thank Martijn Leisink for providing
+the basis on which libDAI has been built.
+
+
+Known issues
+------------
+Due to a bug in GCC 3.3.x and earlier (http://gcc.gnu.org/bugzilla/show_bug.cgi?id=20358)
+it doesn't compile with these versions (it does compile with GCC version 3.4 and higher).
+Workaround: replace the two NAN's in factor.h causing the error messages by -1.
+
+
+Documentation
+-------------
+Almost nonexistant. But I'm working on it. In the meantime, I'll provide limited support
+by email. The following gives an overview of different methods and their properties
+(can be slightly obsolete):
+
+BP
+ updates UpdateType SEQFIX,SEQRND,SEQMAX,PARALL
+ tol double
+ maxiter size_t
+ verbose size_t
+MF
+ tol double
+ maxiter size_t
+ verbose size_t
+HAK
+ clusters MIN,DELTA,LOOP
+ loopdepth
+ doubleloop bool
+ tol double
+ maxiter size_t
+ verbose size_t
+JTREE
+ updates UpdateType HUGIN,SHSH
+ verbose size_t
+MR
+ updates UpdateType FULL,LINEAR
+ inits InitType RESPPROP,CLAMPING,EXACT
+ verbose size_t
+TREEEP
+ type TypeType ORG,ALT
+ tol double
+ maxiter size_t
+ verbose size_t
+LC
+ cavity CavityType FULL,PAIR,PAIR2,UNIFORM
+ updates UpdateType SEQFIX,SEQRND(,NONE)
+ reinit bool
+ cavainame string
+ cavaiopts Properties
+ tol double
+ maxiter size_t
+ verbose size_t
+
+
+
+Quick start
+-----------
+You need:
+- a recent version of gcc (version 3.4 at least)
+- the Boost C++ libraries (under Debian/Ubuntu you can install them using
+ "apt-get install libboost-dev libboost-program-options-dev")
+- GNU make
+
+To build the source, edit the Makefile and then run
+
+ make
+
+If the build was successful, you can test the example program:
+
+ ./example tests/alarm.fg
+
+or the more elaborate test program:
+
+ tests/test --aliases tests/aliases.conf --filename tests/alarm.fg --methods JTREE_HUGIN BP_SEQMAX
+
+A description of the factor graph (.fg) file format can be found in the file FILEFORMAT.
--- /dev/null
+IMPORTANT
+
+- Find out which methods need connected factor graphs (at least TreeEP and JunctionTree).
+
+- Forwardport the Bart patch
+
+- Get rid of multiple inheritance
+
+- Introduce a namespace libDAI
+
+- Investigate the performance penalty from the use of Properties map
+
+- Another design question that needs to be answered: should a BP object own a
+ FactorGraph, or only store a reference to a FactorGraph (which can optimize
+ memory), or should it be a FactorGraph with additional variables and functions?
+ Probably, the first or second option would be the best. Since FactorGraphs are
+ assumed to be rather static, it would make sense to store *references* to
+ FactorGraphs inside AIAlgs. Or maybe some smart_ptr? It would make sense to
+ prevent unnecessary copying of FactorGraphs. Also, general marginalization
+ functions now copy a complete object. Instead, it would make more sense that
+ they construct a new object without copying the FactorGraph. Or they can be made
+ simply methods of the general InfAlg class.
+
+
+DIFFICULT
+
+- Bug: TreeEP sometimes returns NANs.
+
+- Bug: TreeEP::logZ() seems to be wrong (?).
+
+- Bug: strange things happen for at least JTree and TreeEP if the factor graph is not connected.
+ Should be fixed by finding out which methods demand connected factor graphs and these
+ methods should check whether the supplied factor graph is connected.
+
+- Kees' trick for preventing NANs in GBP updates: if quantity smaller than epsilon, replace by epsilon.
+
+- Bug: MF_SEQRANDOM does not work on alarm2.fg (normalization error: Z = 0)
+
+
+OPTIONAL
+
+- Restore brute force EXACT method.
+
+- Alternative implementation of undo factor changes: the only things that have to be
+undone are setting a factor to all ones and setting a factor to a Kronecker delta. This
+could also be implemented in the factor itself, which could maintain its state
+(one/delta/full) and act accordingly.
+
+- Regression test could be enhanced by
+ * reporting number of iterations
+ * reporting all single node marginal errors
+
+- Think about the _fac2ind variable: a better implementation of this feature is
+probably possible.
+
+- Changed() methods?
+
+- It might be an idea to forbid modification of the structure of a factor graph
+ (except in constructors). Then we could do away with Regenerate and similar methods.
+
+- Think about whether using the L-1 norm as a distance measure between new and
+ old messages/beliefs would be better, since || P(x_A, x_B) - Q(x_A, xB) ||_1 < epsilon
+ implies that || P(x_A) - Q(x_A) ||_1 < epsilon. However, for a product of messages
+ we have the weaker || P1(x_A) P2(x_A) - Q1(x_A) Q2(x_A) ||_1 < 2 epsilon if
+ || P1(x_A) - Q1(x_A) ||_1 < epsilon AND || P2(x_A) - Q1(x_A) ||_1 < epsilon...
+
+- Implement damping.
+
+- Optimize GBP with precalculated indices.
+
+- Variable elimination should be implemented generically using a function
+ object that tells you which variable to delete.
+
+- Improve general support for graphs and trees.
+
+- What to do in case of divide-by-zero? For HUGIN, the correct way to handle
+ this is to define x/0 = 0. For other methods, this may not be the case...
+ Maybe a good solution would be to make different flavours of Prob which differ
+ in the way that they handle zeroes (two issues here: normalization of an all-zero
+ factor, which appears to be a sign that something is wrong, and division-by-zero
+ in operator/, operator/= and inverse(), which is probably harmless and can be fixed
+ by defining x/0 = 0. I think it is harmless because usually it occurs by calculating
+ some marginal _in the absence of some factor_ as an intermediate step. Afterwards,
+ this factor is always multiplied in again, which sets any value on these positions
+ to zero. Another difference would be how to deal with underflows: GBP often has
+ numbers that converge to zero, but instead this gives an underflow. A solution is
+ to set everything smaller than some epsilon to zero. Finally, one could decide to
+ add support for sparse potentials. Oh, and whether it can handle negative entries.
+
+- Optimize BP_SEQMAX (it should use a better data structure than a vector for the residuals).
+ Also it does strange things in case of negative potentials.
+
+- Investigate the possibility to calculate logZ *during* the inference algorithm
+ by keeping track of all normalizations.
+
+- Write documentation.
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#include "alldai.h"
+
+
+InfAlg *newInfAlg( const string &name, const FactorGraph &fg, const Properties &opts ) {
+ if( name == BP::Name )
+ return new BP (fg, opts);
+ else if( name == MF::Name )
+ return new MF (fg, opts);
+ else if( name == HAK::Name )
+ return new HAK (fg, opts);
+ else if( name == LC::Name )
+ return new LC (fg, opts);
+ else if( name == TreeEP::Name )
+ return new TreeEP (fg, opts);
+ else if( name == MR::Name )
+ return new MR (fg, opts);
+ else if( name == JTree::Name )
+ return new JTree (fg, opts);
+ else
+ throw "Unknown inference algorithm";
+}
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __ALLDAI_H__
+#define __ALLDAI_H__
+
+
+#include "daialg.h"
+#include "bp.h"
+#include "lc.h"
+#include "hak.h"
+#include "mf.h"
+#include "jtree.h"
+#include "treeep.h"
+#include "mr.h"
+
+
+/// newInfAlg constructs a new approximate inference algorithm named name for the
+/// FactorGraph fg with optionts opts and returns a pointer to the new object.
+/// The caller needs to delete it (maybe some sort of smart_ptr might be useful here).
+InfAlg *newInfAlg( const string &name, const FactorGraph &fg, const Properties &opts );
+
+
+/// AINames contains the names of all approximate inference algorithms
+static const char* DAINames[] = {BP::Name, MF::Name, HAK::Name, LC::Name, TreeEP::Name, MR::Name, JTree::Name};
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __BIPGRAPH_H__
+#define __BIPGRAPH_H__
+
+
+#include <vector>
+#include <algorithm>
+
+
+/// A BipartiteGraph represents a graph with two types of nodes
+template <class T1, class T2> class BipartiteGraph {
+ public:
+ typedef std::pair<size_t,size_t> _edge_t;
+ typedef std::vector<size_t> _nb_t;
+ typedef _nb_t::const_iterator _nb_cit;
+
+
+ private:
+ /// Vertices of first type
+ std::vector<T1> _V1;
+
+ /// Vertices of second type
+ std::vector<T2> _V2;
+
+ /// Edges, which are pairs (v1,v2) with v1 in _V1 and v2 in _V2
+ std::vector<_edge_t> _E12;
+
+ /// Conversion matrix that computes which index of _E12 corresponds to a vertex index pair (v1,v2)
+ std::vector<std::vector<size_t> > _E12ind;
+
+ /// Neighbour indices of vertices of first type
+ std::vector<_nb_t> _nb1;
+
+ /// Neighbour indices of vertices of second type
+ std::vector<_nb_t> _nb2;
+
+
+ public:
+ /// Default constructor
+ BipartiteGraph<T1,T2> () {};
+
+ /// Copy constructor
+ BipartiteGraph<T1,T2> ( const BipartiteGraph<T1,T2> & x ) : _V1(x._V1), _V2(x._V2), _E12(x._E12), _E12ind(x._E12ind), _nb1(x._nb1), _nb2(x._nb2) {};
+
+ /// Assignment operator
+ BipartiteGraph<T1,T2> & operator=(const BipartiteGraph<T1,T2> & x) {
+ if( this != &x ) {
+ _V1 = x._V1;
+ _V2 = x._V2;
+ _E12 = x._E12;
+ _E12ind = x._E12ind;
+ _nb1 = x._nb1;
+ _nb2 = x._nb2;
+ }
+ return *this;
+ }
+
+ /// Provides read access to node of first type
+ const T1 & V1( size_t i1 ) const { return _V1[i1]; }
+ /// Provides full access to node of first type
+ T1 & V1( size_t i1 ) { return _V1[i1]; }
+ /// Provides read access to all nodes of first type
+ const std::vector<T1> & V1s() const { return _V1; }
+ /// Provides full access to all nodes of first type
+ std::vector<T1> & V1s() { return _V1; }
+
+ /// Provides read access to node of second type
+ const T2 & V2( size_t i2 ) const { return _V2[i2]; }
+ /// Provides full access to node of second type
+ T2 & V2( size_t i2 ) { return _V2[i2]; }
+ /// Provides read access to all nodes of second type
+ const std::vector<T2> & V2s() const { return _V2; }
+ /// Provides full access to all nodes of second type
+ std::vector<T2> & V2s() { return _V2; }
+
+ /// Provides read access to edge
+ const _edge_t & edge(size_t ind) const { return _E12[ind]; }
+ /// Provides full access to edge
+ _edge_t & edge(size_t ind) { return _E12[ind]; }
+ /// Provides read access to all edges
+ const std::vector<_edge_t> & edges() const { return _E12; }
+ /// Provides full access to all edges
+ std::vector<_edge_t> & edges() { return _E12; }
+ /// Returns number of edges
+ size_t nr_edges() const { return _E12.size(); }
+
+ /// Provides read access to neighbours of node of first type
+ const _nb_t & nb1( size_t i1 ) const { return _nb1[i1]; }
+ /// Provides full access to neighbours of node of first type
+ _nb_t & nb1( size_t i1 ) { return _nb1[i1]; }
+
+ /// Provides read access to neighbours of node of second type
+ const _nb_t & nb2( size_t i2 ) const { return _nb2[i2]; }
+ /// Provides full access to neighbours of node of second type
+ _nb_t & nb2( size_t i2 ) { return _nb2[i2]; }
+
+ /// Converts the pair of indices (i1,i2) to the corresponding edge index
+ size_t VV2E( const size_t i1, const size_t i2 ) const { return _E12ind[i1][i2]; }
+
+
+ /// Regenerates internal structures (_E12ind, _nb1 and _nb2) based upon _V1, _V2 and _E12
+ void Regenerate() {
+ // Calculate _nb1 and _nb2
+
+ // Start with empty vectors
+ _nb1.clear();
+ _nb1.resize(_V1.size());
+ // Start with empty vectors
+ _nb2.clear();
+ _nb2.resize(_V2.size());
+ // Each edge yields a neighbour pair
+ for( std::vector<_edge_t>::const_iterator e = _E12.begin(); e != _E12.end(); e++ ) {
+ _nb1[e->first].push_back(e->second);
+ _nb2[e->second].push_back(e->first);
+ }
+ // Remove duplicates from _nb1
+ for( size_t i1 = 0; i1 < _V1.size(); i1++ ) {
+ _nb_t::iterator new_end = unique(_nb1[i1].begin(), _nb1[i1].end());
+ _nb1[i1].erase( new_end, _nb1[i1].end() );
+ }
+ // Remove duplicates from _nb2
+ for( size_t i2 = 0; i2 < _V2.size(); i2++ ) {
+ _nb_t::iterator new_end = unique(_nb2[i2].begin(), _nb2[i2].end());
+ _nb2[i2].erase( new_end, _nb2[i2].end() );
+ }
+
+ // Calculate _E12ind
+
+ // Allocate data structures
+ _E12ind.clear();
+ std::vector<size_t> col(_V2.size(),-1);
+ _E12ind.assign(_V1.size(), col);
+ // Assign elements
+ for( size_t k = 0; k < _E12.size(); k++ )
+ _E12ind[_E12[k].first][_E12[k].second] = k;
+ }
+};
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#include <iostream>
+#include <sstream>
+#include <map>
+#include <set>
+#include <algorithm>
+#include "bp.h"
+#include "diffs.h"
+#include "util.h"
+#include "properties.h"
+
+
+using namespace std;
+
+
+const char *BP::Name = "BP";
+
+
+bool BP::checkProperties() {
+ if( !HasProperty("updates") )
+ return false;
+ if( !HasProperty("tol") )
+ return false;
+ if (!HasProperty("maxiter") )
+ return false;
+ if (!HasProperty("verbose") )
+ return false;
+
+ ConvertPropertyTo<double>("tol");
+ ConvertPropertyTo<size_t>("maxiter");
+ ConvertPropertyTo<size_t>("verbose");
+ ConvertPropertyTo<UpdateType>("updates");
+
+ return true;
+}
+
+
+void BP::Regenerate() {
+ DAIAlgFG::Regenerate();
+
+ // clear messages
+ _messages.clear();
+ _messages.reserve(nr_edges());
+
+ // clear indices
+ _indices.clear();
+ _indices.reserve(nr_edges());
+
+ // create messages and indices
+ for( vector<_edge_t>::const_iterator iI=edges().begin(); iI!=edges().end(); iI++ ) {
+ _messages.push_back( Prob( var(iI->first).states() ) );
+
+ vector<size_t> ind( factor(iI->second).stateSpace(), 0 );
+ Index i (var(iI->first), factor(iI->second).vars() );
+ for( size_t j = 0; i >= 0; ++i,++j )
+ ind[j] = i;
+ _indices.push_back( ind );
+ }
+
+ // create new_messages
+ _newmessages = _messages;
+}
+
+
+void BP::init() {
+ assert( checkProperties() );
+ for( vector<Prob>::iterator mij = _messages.begin(); mij != _messages.end(); mij++ )
+ mij->fill(1.0 / mij->size());
+ _newmessages = _messages;
+}
+
+
+void BP::calcNewMessage (size_t iI) {
+ // calculate updated message I->i
+ size_t i = edge(iI).first;
+ size_t I = edge(iI).second;
+
+/* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION
+
+ Factor prod( factor( I ) );
+ for( _nb_cit j = nb2(I).begin(); j != nb2(I).end(); j++ )
+ if( *j != i ) { // for all j in I \ i
+ for( _nb_cit J = nb1(*j).begin(); J != nb1(*j).end(); J++ )
+ if( *J != I ) { // for all J in nb(j) \ I
+ prod *= Factor( *j, message(*j,*J) );
+ Factor marg = prod.marginal(var(i));
+*/
+
+ Prob prod( factor(I).p() );
+
+ // Calculate product of incoming messages and factor I
+ for( _nb_cit j = nb2(I).begin(); j != nb2(I).end(); j++ )
+ if( *j != i ) { // for all j in I \ i
+ // ind is the precalculated Index(j,I) i.e. to x_I == k corresponds x_j == ind[k]
+ _ind_t* ind = &(index(*j,I));
+
+ // prod_j will be the product of messages coming into j
+ Prob prod_j( var(*j).states() );
+ for( _nb_cit J = nb1(*j).begin(); J != nb1(*j).end(); J++ )
+ if( *J != I ) // for all J in nb(j) \ I
+ prod_j *= message(*j,*J);
+
+ // multiply prod with prod_j
+ for( size_t r = 0; r < prod.size(); r++ )
+ prod[r] *= prod_j[(*ind)[r]];
+ }
+
+ // Marginalize onto i
+ Prob marg( var(i).states(), 0.0 );
+ // ind is the precalculated Index(i,I) i.e. to x_I == k corresponds x_i == ind[k]
+ _ind_t* ind = &(index(i,I));
+ for( size_t r = 0; r < prod.size(); r++ )
+ marg[(*ind)[r]] += prod[r];
+ marg.normalize( _normtype );
+
+ // Store result
+ _newmessages[iI] = marg;
+}
+
+
+// BP::run does not check for NANs for performance reasons
+// Somehow NaNs do not often occur in BP...
+double BP::run() {
+ if( Verbose() >= 1 )
+ cout << "Starting " << identify() << "...";
+ if( Verbose() >= 3)
+ cout << endl;
+
+ clock_t tic = toc();
+ Diffs diffs(nrVars(), 1.0);
+
+ vector<size_t> edge_seq;
+ vector<double> residuals;
+
+ vector<Factor> old_beliefs;
+ old_beliefs.reserve( nrVars() );
+ for( size_t i = 0; i < nrVars(); i++ )
+ old_beliefs.push_back(belief1(i));
+
+ size_t iter = 0;
+
+ if( Updates() == UpdateType::SEQMAX ) {
+ // do the first pass
+ for(size_t iI = 0; iI < nr_edges(); iI++ )
+ calcNewMessage(iI);
+
+ // calculate initial residuals
+ residuals.reserve(nr_edges());
+ for( size_t iI = 0; iI < nr_edges(); iI++ )
+ residuals.push_back( dist( _newmessages[iI], _messages[iI], Prob::DISTLINF ) );
+ } else {
+ edge_seq.reserve( nr_edges() );
+ for( size_t i = 0; i < nr_edges(); i++ )
+ edge_seq.push_back( i );
+ }
+
+ // do several passes over the network until maximum number of iterations has
+ // been reached or until the maximum belief difference is smaller than tolerance
+ for( iter=0; iter < MaxIter() && diffs.max() > Tol(); iter++ ) {
+ if( Updates() == UpdateType::SEQMAX ) {
+ // Residuals-BP by Koller et al.
+ for( size_t t = 0; t < nr_edges(); t++ ) {
+ // update the message with the largest residual
+ size_t iI = max_element(residuals.begin(), residuals.end()) - residuals.begin();
+ _messages[iI] = _newmessages[iI];
+ residuals[iI] = 0;
+
+ // I->i has been updated, which means that residuals for all
+ // J->j with J in nb[i]\I and j in nb[J]\i have to be updated
+ size_t i = edge(iI).first;
+ size_t I = edge(iI).second;
+ for( _nb_cit J = nb1(i).begin(); J != nb1(i).end(); J++ )
+ if( *J != I )
+ for( _nb_cit j = nb2(*J).begin(); j != nb2(*J).end(); j++ )
+ if( *j != i ) {
+ size_t jJ = VV2E(*j,*J);
+ calcNewMessage(jJ);
+ residuals[jJ] = dist( _newmessages[jJ], _messages[jJ], Prob::DISTLINF );
+ }
+ }
+ } else if( Updates() == UpdateType::PARALL ) {
+ // Parallel updates
+ for( size_t t = 0; t < nr_edges(); t++ )
+ calcNewMessage(t);
+
+ for( size_t t = 0; t < nr_edges(); t++ )
+ _messages[t] = _newmessages[t];
+ } else {
+ // Sequential updates
+ if( Updates() == UpdateType::SEQRND )
+ random_shuffle( edge_seq.begin(), edge_seq.end() );
+
+ for( size_t t = 0; t < nr_edges(); t++ ) {
+ size_t k = edge_seq[t];
+ calcNewMessage(k);
+ _messages[k] = _newmessages[k];
+ }
+ }
+
+ // calculate new beliefs and compare with old ones
+ for( size_t i = 0; i < nrVars(); i++ ) {
+ Factor nb( belief1(i) );
+ diffs.push( dist( nb, old_beliefs[i], Prob::DISTLINF ) );
+ old_beliefs[i] = nb;
+ }
+
+ if( Verbose() >= 3 )
+ cout << "BP::run: maxdiff " << diffs.max() << " after " << iter+1 << " passes" << endl;
+ }
+
+ updateMaxDiff( diffs.max() );
+
+ if( Verbose() >= 1 ) {
+ if( diffs.max() > Tol() ) {
+ if( Verbose() == 1 )
+ cout << endl;
+ cout << "BP::run: WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.max() << endl;
+ } else {
+ if( Verbose() >= 3 )
+ cout << "BP::run: ";
+ cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
+ }
+ }
+
+ return diffs.max();
+}
+
+
+Factor BP::belief1( size_t i ) const {
+ Prob prod( var(i).states() );
+ for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); I++ )
+ prod *= newMessage(i,*I);
+
+ prod.normalize( Prob::NORMPROB );
+ return( Factor( var(i), prod ) );
+}
+
+
+Factor BP::belief (const Var &n) const {
+ return( belief1( findVar( n ) ) );
+}
+
+
+vector<Factor> BP::beliefs() const {
+ vector<Factor> result;
+ for( size_t i = 0; i < nrVars(); i++ )
+ result.push_back( belief1(i) );
+ for( size_t I = 0; I < nrFactors(); I++ )
+ result.push_back( belief2(I) );
+ return result;
+}
+
+
+Factor BP::belief( const VarSet &ns ) const {
+ if( ns.size() == 1 )
+ return belief( *(ns.begin()) );
+ else {
+ size_t I;
+ for( I = 0; I < nrFactors(); I++ )
+ if( factor(I).vars() >> ns )
+ break;
+ assert( I != nrFactors() );
+ return belief2(I).marginal(ns);
+ }
+}
+
+
+Factor BP::belief2 (size_t I) const {
+ Prob prod( factor(I).p() );
+
+ for( _nb_cit j = nb2(I).begin(); j != nb2(I).end(); j++ ) {
+ // ind is the precalculated Index(j,I) i.e. to x_I == k corresponds x_j == ind[k]
+ const _ind_t *ind = &(index(*j, I));
+
+ // prod_j will be the product of messages coming into j
+ Prob prod_j( var(*j).states() );
+ for( _nb_cit J = nb1(*j).begin(); J != nb1(*j).end(); J++ )
+ if( *J != I ) // for all J in nb(j) \ I
+ prod_j *= newMessage(*j,*J);
+
+ // multiply prod with prod_j
+ for( size_t r = 0; r < prod.size(); r++ )
+ prod[r] *= prod_j[(*ind)[r]];
+ }
+
+ Factor result( factor(I).vars(), prod );
+ result.normalize( Prob::NORMPROB );
+
+ return( result );
+
+/* UNOPTIMIZED VERSION
+
+ Factor prod( factor(I) );
+ for( _nb_cit i = nb2(I).begin(); i != nb2(I).end(); i++ ) {
+ for( _nb_cit J = nb1(*i).begin(); J != nb1(*i).end(); J++ )
+ if( *J != I )
+ prod *= Factor( var(*i), newMessage(*i,*J)) );
+ }
+ return prod.normalize( Prob::NORMPROB );*/
+}
+
+
+Complex BP::logZ() const {
+ Complex sum = 0.0;
+ for(size_t i = 0; i < nrVars(); i++ )
+ sum += Complex(1.0 - nb1(i).size()) * belief1(i).entropy();
+ for( size_t I = 0; I < nrFactors(); I++ )
+ sum -= KL_dist( belief2(I), factor(I) );
+ return sum;
+}
+
+
+string BP::identify() const {
+ stringstream result (stringstream::out);
+ result << Name << GetProperties();
+ return result.str();
+}
+
+
+void BP::init( const VarSet &ns ) {
+ for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ ) {
+ size_t ni = findVar( *n );
+ for( _nb_cit I = nb1(ni).begin(); I != nb1(ni).end(); I++ )
+ message(ni,*I).fill( 1.0 );
+ }
+}
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __BP_H__
+#define __BP_H__
+
+
+#include "daialg.h"
+#include "factorgraph.h"
+#include "enum.h"
+
+
+using namespace std;
+
+
+class BP : public DAIAlgFG {
+ protected:
+ typedef vector<size_t> _ind_t;
+
+ vector<_ind_t> _indices;
+ vector<Prob> _messages, _newmessages;
+
+ public:
+ ENUM4(UpdateType,SEQFIX,SEQRND,SEQMAX,PARALL)
+ UpdateType Updates() const { return GetPropertyAs<UpdateType>("updates"); }
+
+ // default constructor
+ BP() : DAIAlgFG() {};
+ // copy constructor
+ BP(const BP & x) : DAIAlgFG(x), _indices(x._indices), _messages(x._messages), _newmessages(x._newmessages) {};
+ BP* clone() const { return new BP(*this); }
+ // construct BP object from FactorGraph
+ BP(const FactorGraph & fg, const Properties &opts) : DAIAlgFG(fg, opts) {
+ assert( checkProperties() );
+ Regenerate();
+ }
+ // assignment operator
+ BP & operator=(const BP & x) {
+ if(this!=&x) {
+ DAIAlgFG::operator=(x);
+ _messages = x._messages;
+ _newmessages = x._newmessages;
+ _indices = x._indices;
+ }
+ return *this;
+ }
+
+ static const char *Name;
+
+ Prob & message(size_t i1, size_t i2) { return( _messages[VV2E(i1,i2)] ); }
+ const Prob & message(size_t i1, size_t i2) const { return( _messages[VV2E(i1,i2)] ); }
+ Prob & newMessage(size_t i1, size_t i2) { return( _newmessages[VV2E(i1,i2)] ); }
+ const Prob & newMessage(size_t i1, size_t i2) const { return( _newmessages[VV2E(i1,i2)] ); }
+ _ind_t & index(size_t i1, size_t i2) { return( _indices[VV2E(i1,i2)] ); }
+ const _ind_t & index(size_t i1, size_t i2) const { return( _indices[VV2E(i1,i2)] ); }
+
+ string identify() const;
+ void Regenerate();
+ void init();
+ void calcNewMessage(size_t iI);
+ double run();
+ Factor belief1 (size_t i) const;
+ Factor belief2 (size_t I) const;
+ Factor belief (const Var &n) const;
+ Factor belief (const VarSet &n) const;
+ vector<Factor> beliefs() const;
+ Complex logZ() const;
+
+ void init( const VarSet &ns );
+ void undoProbs( const VarSet &ns ) { FactorGraph::undoProbs(ns); init(ns); }
+ bool checkProperties();
+};
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#include <set>
+#include <vector>
+#include <iostream>
+#include "varset.h"
+#include "clustergraph.h"
+
+
+using namespace std;
+
+
+ClusterGraph ClusterGraph::VarElim( const vector<Var> & ElimSeq ) const {
+ const long verbose = 0;
+
+ // Make a copy
+ ClusterGraph _Cl(*this);
+
+ ClusterGraph result;
+ _Cl.eraseNonMaximal();
+
+ // Do variable elimination
+ for( vector<Var>::const_iterator n = ElimSeq.begin(); n != ElimSeq.end(); n++ ) {
+ assert( _Cl.vars() && *n );
+
+ if( verbose >= 1 )
+ cout << "Cost of eliminating " << *n << ": " << _Cl.eliminationCost( *n ) << " new edges" << endl;
+
+ result.insert( _Cl.Delta(*n) );
+
+ if( verbose >= 1 )
+ cout << "_Cl = " << _Cl << endl;
+
+ if( verbose >= 1 )
+ cout << "After inserting " << _Cl.delta(*n) << ", _Cl = ";
+ _Cl.insert( _Cl.delta(*n) );
+ if( verbose >= 1 )
+ cout << _Cl << endl;
+
+ if( verbose >= 1 )
+ cout << "After erasing clusters that contain " << *n << ", _Cl = ";
+ _Cl.eraseSubsuming( *n );
+ if( verbose >= 1 )
+ cout << _Cl << endl;
+
+ if( verbose >= 1 )
+ cout << "After erasing nonmaximal clusters, _Cl = ";
+ _Cl.eraseNonMaximal();
+ if( verbose >= 1 )
+ cout << _Cl << endl;
+ }
+
+ return result;
+}
+
+
+ClusterGraph ClusterGraph::VarElim_MinFill() const {
+ const long verbose = 0;
+
+ // Make a copy
+ ClusterGraph _Cl(*this);
+ VarSet _vars( vars() );
+
+ ClusterGraph result;
+ _Cl.eraseNonMaximal();
+
+ // Do variable elimination
+ while( !_vars.empty() ) {
+ if( verbose >= 1 )
+ cout << "Var Eliminiation cost" << endl;
+ VarSet::iterator lowest = _vars.end();
+ size_t lowest_cost = -1UL;
+ for( VarSet::const_iterator n = _vars.begin(); n != _vars.end(); n++ ) {
+ size_t cost = _Cl.eliminationCost( *n );
+ if( verbose >= 1 )
+ cout << *n << " " << cost << endl;
+ if( lowest == _vars.end() || lowest_cost > cost ) {
+ lowest = n;
+ lowest_cost = cost;
+ }
+ }
+ Var n = *lowest;
+
+ if( verbose >= 1 )
+ cout << "Lowest: " << n << " (" << lowest_cost << ")" << endl;
+
+ result.insert( _Cl.Delta(n) );
+
+ _Cl.insert( _Cl.delta(n) );
+ _Cl.eraseSubsuming( n );
+ _Cl.eraseNonMaximal();
+ _vars /= n;
+
+ }
+
+ return result;
+}
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __CLUSTERGRAPH_H__
+#define __CLUSTERGRAPH_H__
+
+
+#include <set>
+#include <vector>
+#include "varset.h"
+
+
+using namespace std;
+
+
+/// A ClusterGraph is a hypergraph with VarSets as nodes.
+/// It is implemented as a set<VarSet> in which the adjacency
+/// relationship is computed realtime. It may be better to
+/// implement it as a bipartitegraph, though. The additional
+/// functionality compared to a simple set<VarSet> is
+/// finding maximal clusters, finding cliques, etc...
+class ClusterGraph : public set<VarSet> {
+ public:
+ /// Default constructor
+ ClusterGraph() : set<VarSet>() {}
+
+ /// Construct from vector<VarSet>
+ ClusterGraph( const vector<VarSet> & cls ) {
+ insert( cls.begin(), cls.end() );
+ }
+
+ /// Copy constructor
+ ClusterGraph(const ClusterGraph &x) : set<VarSet>(x) {}
+
+ /// Assignment operator
+ ClusterGraph& operator=( const ClusterGraph &x ) {
+ if( this != &x ) {
+ set<VarSet>::operator=( x );
+ }
+ return *this;
+ }
+
+ /// Returns true if ns is a maximal member of *this under inclusion (VarSet::operator<<)
+ bool isMaximal( const VarSet &ns ) const {
+ if( count( ns ) ) {
+ // ns is a member
+ bool maximal = true;
+ for( const_iterator x = begin(); x != end() && maximal; x++ )
+ if( (ns << *x) && (ns != *x) )
+ maximal = false;
+ return maximal;
+ } else
+ return false;
+ }
+
+ /// Erase all VarSets that are not maximal
+ ClusterGraph& eraseNonMaximal() {
+ for( iterator x = begin(); x != end(); )
+ if( !isMaximal(*x) )
+ erase(x++);
+ else
+ x++;
+ return *this;
+ }
+
+ /// Return union of all members
+ VarSet vars() const {
+ VarSet result;
+ for( iterator x = begin(); x != end(); x++ )
+ result |= *x;
+ return result;
+ }
+
+ /// Returns true if n1 and n2 are adjacent, i.e.\ by
+ /// definition, are both contained in some cluster in *this
+ bool adj( const Var& n1, const Var& n2 ) {
+ bool result = false;
+ for( iterator x = begin(); (x != end()) && (!result); x++ )
+ if( (*x && n1) && (*x && n2) )
+ result = true;
+ return result;
+ }
+
+ /// Returns union of clusters that contain n, minus n
+ VarSet Delta( const Var& n ) const {
+ VarSet result;
+ for( iterator x = begin(); x != end(); x++ )
+ if( (*x && n) )
+ result |= *x;
+ return result;
+ }
+
+ /// Returns union of clusters that contain n, minus n
+ VarSet delta( const Var& n ) const {
+ return Delta( n ) / n;
+ }
+
+ /// Erases all members that contain n
+ ClusterGraph& eraseSubsuming( const Var& n ) {
+ for( iterator x = begin(); x != end(); )
+ if( (*x && n) )
+ erase(x++);
+ else
+ x++;
+ return *this;
+ }
+
+ /// Send to output stream
+ friend std::ostream & operator << ( std::ostream & os, const ClusterGraph & cl ) {
+ os << "{";
+ ClusterGraph::const_iterator x = cl.begin();
+ if( x != cl.end() )
+ os << *(x++);
+ for( ; x != cl.end(); x++ )
+ os << ", " << *x;
+ os << "}";
+ return os;
+ }
+
+ /// Convert to vector<VarSet>
+ vector<VarSet> toVector() const {
+ vector<VarSet> result;
+ result.reserve( size() );
+ for( const_iterator x = begin(); x != end(); x++ )
+ result.push_back( *x );
+ return result;
+ }
+
+ /// Calculate cost of eliminating variable n
+ /// using as a measure "number of added edges in the adjacency graph"
+ /// where the adjacency graph has the variables as its nodes and
+ /// connects nodes i1 and i2 iff i1 and i2 occur in some common factor I
+ size_t eliminationCost( const Var& n ) {
+ VarSet d_n = delta( n );
+ size_t cost = 0;
+
+ // for each unordered pair {i1,i2} adjacent to n
+ for( VarSet::const_iterator i1 = d_n.begin(); i1 != d_n.end(); i1++ ) {
+ VarSet d_i1 = delta( *i1 );
+ for( VarSet::const_iterator i2 = i1; (++i2) != d_n.end(); ) {
+ // if i1 and i2 are not adjacent, eliminating n would make them adjacent
+ if( !adj(*i1, *i2) )
+ cost++;
+ }
+ }
+ return cost;
+ }
+
+ /// Perform Variable Elimination without Probs, i.e. only keeping track of
+ /// the interactions that are created along the way.
+ /// Input: a set of outer clusters and an elimination sequence
+ /// Output: a set of elimination "cliques"
+ ClusterGraph VarElim( const vector<Var> &ElimSeq ) const;
+
+ /// As Taylan does it
+ ClusterGraph VarElim_MinFill() const;
+};
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#include "daialg.h"
+
+
+/// Calculate the marginal of obj on ns by clamping
+/// all variables in ns and calculating logZ for each joined state
+Factor calcMarginal( const InfAlg & obj, const VarSet & ns, bool reInit ) {
+ Factor Pns (ns);
+
+ multind mi( ns );
+
+ InfAlg *clamped = obj.clone();
+ if( !reInit )
+ clamped->init();
+
+ Complex logZ0;
+ for( size_t j = 0; j < mi.max(); j++ ) {
+ // save unclamped factors connected to ns
+ clamped->saveProbs( ns );
+
+ // set clamping Factors to delta functions
+ vector<size_t> vi = mi.vi( j );
+ size_t k = 0;
+ for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++, k++ )
+ clamped->clamp( *n, vi[k] );
+
+ // run DAIAlg, calc logZ, store in Pns
+ if( clamped->Verbose() >= 2 )
+ cout << j << ": ";
+ if( reInit )
+ clamped->init();
+ clamped->run();
+
+ Complex Z;
+ if( j == 0 ) {
+ logZ0 = clamped->logZ();
+ Z = 1.0;
+ } else {
+ // subtract logZ0 to avoid very large numbers
+ Z = exp(clamped->logZ() - logZ0);
+ if( fabs(imag(Z)) > 1e-5 )
+ cout << "Marginal:: WARNING: complex Z (" << Z << ")" << endl;
+ }
+
+ Pns[j] = real(Z);
+
+ // restore clamped factors
+ clamped->undoProbs( ns );
+ }
+
+ delete clamped;
+
+ return( Pns.normalized(Prob::NORMPROB) );
+}
+
+
+vector<Factor> calcPairBeliefs( const InfAlg & obj, const VarSet& ns, bool reInit ) {
+ // convert ns to vector<VarSet>
+ size_t N = ns.size();
+ vector<Var> vns;
+ vns.reserve( N );
+ for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ )
+ vns.push_back( *n );
+
+ vector<Factor> pairbeliefs;
+ pairbeliefs.reserve( N * N );
+ for( size_t j = 0; j < N; j++ )
+ for( size_t k = 0; k < N; k++ )
+ if( j == k )
+ pairbeliefs.push_back(Factor());
+ else
+ pairbeliefs.push_back(Factor(vns[j] | vns[k]));
+
+ InfAlg *clamped = obj.clone();
+ if( !reInit )
+ clamped->init();
+
+ Complex logZ0;
+ for( size_t j = 0; j < N; j++ ) {
+ // clamp Var j to its possible values
+ for( size_t j_val = 0; j_val < vns[j].states(); j_val++ ) {
+ if( obj.Verbose() >= 2 )
+ cout << j << "/" << N-1 << " (" << j_val << "/" << vns[j].states() << "): ";
+
+ // save unclamped factors connected to ns
+ clamped->saveProbs( ns );
+
+ clamped->clamp( vns[j], j_val );
+ if( reInit )
+ clamped->init();
+ clamped->run();
+
+ //if( j == 0 )
+ // logZ0 = obj.logZ();
+ double Z_xj = 1.0;
+ if( j == 0 && j_val == 0 ) {
+ logZ0 = clamped->logZ();
+ } else {
+ // subtract logZ0 to avoid very large numbers
+ Complex Z = exp(clamped->logZ() - logZ0);
+ if( fabs(imag(Z)) > 1e-5 )
+ cout << "calcPairBelief:: Warning: complex Z: " << Z << endl;
+ Z_xj = real(Z);
+ }
+
+ for( size_t k = 0; k < N; k++ )
+ if( k != j ) {
+ Factor b_k = clamped->belief(vns[k]);
+ for( size_t k_val = 0; k_val < vns[k].states(); k_val++ )
+ if( vns[j].label() < vns[k].label() )
+ pairbeliefs[j * N + k][j_val + (k_val * vns[j].states())] = Z_xj * b_k[k_val];
+ else
+ pairbeliefs[j * N + k][k_val + (j_val * vns[k].states())] = Z_xj * b_k[k_val];
+ }
+
+ // restore clamped factors
+ clamped->undoProbs( ns );
+ }
+ }
+
+ delete clamped;
+
+ // Calculate result by taking the geometric average
+ vector<Factor> result;
+ result.reserve( N * (N - 1) / 2 );
+ for( size_t j = 0; j < N; j++ )
+ for( size_t k = j+1; k < N; k++ )
+ result.push_back( (pairbeliefs[j * N + k] * pairbeliefs[k * N + j]) ^ 0.5 );
+
+ return result;
+}
+
+
+Factor calcMarginal2ndO( const InfAlg & obj, const VarSet& ns, bool reInit ) {
+ // returns a a probability distribution whose 1st order interactions
+ // are unspecified, whose 2nd order interactions approximate those of
+ // the marginal on ns, and whose higher order interactions are absent.
+
+ vector<Factor> pairbeliefs = calcPairBeliefs( obj, ns, reInit );
+
+ Factor Pns (ns);
+ for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
+ Pns *= pairbeliefs[ij];
+
+ return( Pns.normalized(Prob::NORMPROB) );
+}
+
+
+vector<Factor> calcPairBeliefsNew( const InfAlg & obj, const VarSet& ns, bool reInit ) {
+ vector<Factor> result;
+ result.reserve( ns.size() * (ns.size() - 1) / 2 );
+
+ InfAlg *clamped = obj.clone();
+ if( !reInit )
+ clamped->init();
+
+ Complex logZ0;
+ VarSet::const_iterator nj = ns.begin();
+ for( long j = 0; j < (long)ns.size() - 1; j++, nj++ ) {
+ size_t k = 0;
+ for( VarSet::const_iterator nk = nj; (++nk) != ns.end(); k++ ) {
+ Factor pairbelief( *nj | *nk );
+
+ // clamp Vars j and k to their possible values
+ for( size_t j_val = 0; j_val < nj->states(); j_val++ )
+ for( size_t k_val = 0; k_val < nk->states(); k_val++ ) {
+ // save unclamped factors connected to ns
+ clamped->saveProbs( ns );
+
+ clamped->clamp( *nj, j_val );
+ clamped->clamp( *nk, k_val );
+ if( reInit )
+ clamped->init();
+ clamped->run();
+
+ double Z_xj = 1.0;
+ if( j_val == 0 && k_val == 0 ) {
+ logZ0 = clamped->logZ();
+ } else {
+ // subtract logZ0 to avoid very large numbers
+ Complex Z = exp(clamped->logZ() - logZ0);
+ if( fabs(imag(Z)) > 1e-5 )
+ cout << "calcPairBelief:: Warning: complex Z: " << Z << endl;
+ Z_xj = real(Z);
+ }
+
+ // we assume that j.label() < k.label()
+ // i.e. we make an assumption here about the indexing
+ pairbelief[j_val + (k_val * nj->states())] = Z_xj;
+
+ // restore clamped factors
+ clamped->undoProbs( ns );
+ }
+
+ result.push_back( pairbelief );
+ }
+ }
+
+ delete clamped;
+
+ assert( result.size() == (ns.size() * (ns.size() - 1) / 2) );
+
+ return result;
+}
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __DAIALG_H__
+#define __DAIALG_H__
+
+
+#include <string>
+#include <iostream>
+#include <vector>
+#include "factorgraph.h"
+#include "regiongraph.h"
+#include "properties.h"
+
+
+/// The InfAlg class is the common denominator of the various approximate inference algorithms.
+/// A InfAlg object represents a discrete factorized probability distribution over multiple variables
+/// together with an inference algorithm.
+class InfAlg {
+ private:
+ /// Properties of the algorithm (replaces _tol, _maxiter, _verbose)
+ Properties _properties;
+
+ /// Maximum difference encountered so far
+ double _maxdiff;
+
+
+ public:
+ /// Default constructor
+ InfAlg() : _properties(), _maxdiff(0.0) {}
+
+ /// Constructor with options
+ InfAlg( const Properties &opts ) : _properties(opts), _maxdiff(0.0) {}
+
+ /// Copy constructor
+ InfAlg( const InfAlg & x ) : _properties(x._properties), _maxdiff(x._maxdiff) {}
+
+ /// Clone (virtual copy constructor)
+ virtual InfAlg* clone() const = 0;
+
+ /// Assignment operator
+ InfAlg & operator=( const InfAlg & x ) {
+ if( this != &x ) {
+ _properties = x._properties;
+ _maxdiff = x._maxdiff;
+ }
+ return *this;
+ }
+
+ /// Virtual desctructor
+ // (this is needed because this class contains virtual functions)
+ virtual ~InfAlg() {}
+
+ /// Returns true if a property with the given key is present
+ bool HasProperty(const PropertyKey &key) const { return _properties.hasKey(key); }
+
+ /// Gets a property
+ const PropertyValue & GetProperty(const PropertyKey &key) const { return _properties.Get(key); }
+
+ /// Gets a property, casted as ValueType
+ template<typename ValueType>
+ ValueType GetPropertyAs(const PropertyKey &key) const { return _properties.GetAs<ValueType>(key); }
+
+ /// Sets a property
+ void SetProperty(const PropertyKey &key, const PropertyValue &val) { _properties[key] = val; }
+
+ /// Converts a property from string to ValueType, if necessary
+ template<typename ValueType>
+ void ConvertPropertyTo(const PropertyKey &key) { _properties.ConvertTo<ValueType>(key); }
+
+ /// Gets all properties
+ const Properties & GetProperties() const { return _properties; }
+
+ /// Sets properties
+ void SetProperties(const Properties &p) { _properties = p; }
+
+ /// Sets tolerance
+ void Tol( double tol ) { SetProperty("tol", tol); }
+ /// Gets tolerance
+ double Tol() const { return GetPropertyAs<double>("tol"); }
+
+ /// Sets maximum number of iterations
+ void MaxIter( size_t maxiter ) { SetProperty("maxiter", maxiter); }
+ /// Gets maximum number of iterations
+ size_t MaxIter() const { return GetPropertyAs<size_t>("maxiter"); }
+
+ /// Sets verbosity
+ void Verbose( size_t verbose ) { SetProperty("verbose", verbose); }
+ /// Gets verbosity
+ size_t Verbose() const { return GetPropertyAs<size_t>("verbose"); }
+
+ /// Sets maximum difference encountered so far
+ void MaxDiff( double maxdiff ) { _maxdiff = maxdiff; }
+ /// Gets maximum difference encountered so far
+ double MaxDiff() const { return _maxdiff; }
+ /// Updates maximum difference encountered so far
+ void updateMaxDiff( double maxdiff ) { if( maxdiff > _maxdiff ) _maxdiff = maxdiff; }
+ /// Sets maximum difference encountered so far to zero
+ void clearMaxDiff() { _maxdiff = 0.0; }
+
+ /// Identifies itself for logging purposes
+ virtual std::string identify() const = 0;
+
+ /// Get single node belief
+ virtual Factor belief( const Var &n ) const = 0;
+
+ /// Get general belief
+ virtual Factor belief( const VarSet &n ) const = 0;
+
+ /// Get all beliefs
+ virtual vector<Factor> beliefs() const = 0;
+
+ /// Get log partition sum
+ virtual Complex logZ() const = 0;
+
+ /// Clear messages and beliefs
+ virtual void init() = 0;
+
+ /// The actual approximate inference algorithm
+ virtual double run() = 0;
+
+ /// Save factor I
+ virtual void saveProb( size_t I ) = 0;
+ /// Save Factors involving ns
+ virtual void saveProbs( const VarSet &ns ) = 0;
+
+ /// Restore factor I
+ virtual void undoProb( size_t I ) = 0;
+ /// Restore Factors involving ns
+ virtual void undoProbs( const VarSet &ns ) = 0;
+
+ /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$)
+ virtual void clamp( const Var & n, size_t i ) = 0;
+
+ /// Return all variables that interact with n
+ virtual VarSet delta( const Var & n ) const = 0;
+
+ /// Set all factors interacting with n to 1
+ virtual void makeCavity( const Var & n ) = 0;
+
+ /// Set factor I to 1
+ virtual void makeFactorCavity( size_t I ) = 0;
+
+ /// Get index of variable n
+ virtual size_t findVar( const Var & n ) const = 0;
+
+ /// Get index of first factor involving ns
+ virtual size_t findFactor( const VarSet &ns ) const = 0;
+
+ /// Get number of variables
+ virtual size_t nrVars() const = 0;
+
+ /// Get number of factors
+ virtual size_t nrFactors() const = 0;
+
+ /// Get const reference to variable i
+ virtual const Var & var(size_t i) const = 0;
+
+ /// Get reference to variable i
+ virtual Var & var(size_t i) = 0;
+
+ /// Get const reference to factor I
+ virtual const Factor & factor( size_t I ) const = 0;
+
+ /// Get reference to factor I
+ virtual Factor & factor( size_t I ) = 0;
+
+ /// Factor I has been updated
+ virtual void updatedFactor( size_t I ) = 0;
+
+ /// Checks whether all necessary properties have been set
+ /// and casts string-valued properties to other values if necessary
+ virtual bool checkProperties() = 0;
+};
+
+
+template <class T>
+class DAIAlg : public InfAlg, public T {
+ public:
+ /// Default constructor
+ DAIAlg() : InfAlg(), T() {}
+
+ /// Construct DAIAlg with empty T but using the specified properties
+ DAIAlg( const Properties &opts ) : InfAlg( opts ), T() {}
+
+ /// Construct DAIAlg using the specified properties
+ DAIAlg( const T & t, const Properties &opts ) : InfAlg( opts ), T(t) {}
+
+ /// Copy constructor
+ DAIAlg( const DAIAlg & x ) : InfAlg(x), T(x) {}
+
+ /// Save factor I (using T::saveProb)
+ void saveProb( size_t I ) { T::saveProb( I ); }
+ /// Save Factors involving ns (using T::saveProbs)
+ void saveProbs( const VarSet &ns ) { T::saveProbs( ns ); }
+
+ /// Restore factor I (using T::undoProb)
+ void undoProb( size_t I ) { T::undoProb( I ); }
+ /// Restore Factors involving ns (using T::undoProbs)
+ void undoProbs( const VarSet &ns ) { T::undoProbs( ns ); }
+
+ /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$) (using T::clamp)
+ void clamp( const Var & n, size_t i ) { T::clamp( n, i ); }
+
+ /// Return all variables that interact with n (using T::delta)
+ VarSet delta( const Var & n ) const { return T::delta(n); }
+
+ /// Set all factors interacting with n to 1 (using T::makeCavity)
+ void makeCavity( const Var & n ) { T::makeCavity(n); }
+
+ /// Set factor I to 1 (using T::makeFactorCavity)
+ void makeFactorCavity( size_t I ) { T::makeFactorCavity(I); }
+
+ /// Get index of variable n (using T::findVar)
+ size_t findVar( const Var & n ) const { return T::findVar(n); }
+
+ /// Get index of first factor involving ns (using T::findFactor)
+ size_t findFactor( const VarSet &ns ) const { return T::findFactor(ns); }
+
+ /// Get number of variables (using T::nrFactors)
+ size_t nrVars() const { return T::nrVars(); }
+
+ /// Get number of factors (using T::nrFactors)
+ size_t nrFactors() const { return T::nrFactors(); }
+
+ /// Get const reference to variable i (using T::var)
+ const Var & var(size_t i) const { return T::var(i); }
+
+ /// Get reference to variable i (using T::var)
+ Var & var(size_t i) { return T::var(i); }
+
+ /// Get const reference to factor I (using T::factor)
+ const Factor & factor( size_t I ) const { return T::factor(I); }
+
+ /// Get reference to factor I (using T::factor)
+ Factor & factor( size_t I ) { return T::factor(I); }
+
+ /// Factor I has been updated (using T::updatedFactor)
+ void updatedFactor( size_t I ) { T::updatedFactor(I); }
+};
+
+
+typedef DAIAlg<FactorGraph> DAIAlgFG;
+typedef DAIAlg<RegionGraph> DAIAlgRG;
+
+
+/// Calculate the marginal of obj on ns by clamping
+/// all variables in ns and calculating logZ for each joined state
+Factor calcMarginal( const InfAlg & obj, const VarSet & ns, bool reInit );
+
+
+/// Calculate beliefs of all pairs in ns (by clamping
+/// nodes in ns and calculating logZ and the beliefs for each state)
+vector<Factor> calcPairBeliefs( const InfAlg & obj, const VarSet& ns, bool reInit );
+
+
+/// Calculate beliefs of all pairs in ns (by clamping
+/// pairs in ns and calculating logZ for each joined state)
+vector<Factor> calcPairBeliefsNew( const InfAlg & obj, const VarSet& ns, bool reInit );
+
+
+/// Calculate 2nd order interactions of the marginal of obj on ns
+Factor calcMarginal2ndO( const InfAlg & obj, const VarSet& ns, bool reInit );
+
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __DIFFS_H__
+#define __DIFFS_H__
+
+
+#include <vector>
+
+using namespace std;
+
+
+class Diffs : public vector<double> {
+ private:
+ size_t _maxsize;
+ double _def;
+ vector<double>::iterator _pos;
+ vector<double>::iterator _maxpos;
+ public:
+ Diffs(long maxsize, double def) : vector<double>(), _maxsize(maxsize), _def(def) {
+ this->reserve(_maxsize);
+ _pos = begin();
+ _maxpos = begin();
+ };
+ double max() {
+ if( size() < _maxsize )
+ return _def;
+ else
+ return( *_maxpos );
+ }
+ void push(double x) {
+ if( size() < _maxsize ) {
+ push_back(x);
+ _pos = end();
+ _maxpos = max_element(begin(),end());
+ }
+ else {
+ if( _pos == end() )
+ _pos = begin();
+ if( _maxpos == _pos ) {
+ *_pos++ = x;
+ _maxpos = max_element(begin(),end());
+ } else {
+ if( x > *_maxpos )
+ _maxpos = _pos;
+ *_pos++ = x;
+ }
+ }
+ }
+};
+
+#endif
--- /dev/null
+# Doxyfile 1.4.2
+
+# This file describes the settings to be used by the documentation system
+# doxygen (www.doxygen.org) for a project
+#
+# All text after a hash (#) is considered a comment and will be ignored
+# The format is:
+# TAG = value [value, ...]
+# For lists items can also be appended using:
+# TAG += value [value, ...]
+# Values that contain spaces should be placed between quotes (" ")
+
+#---------------------------------------------------------------------------
+# Project related configuration options
+#---------------------------------------------------------------------------
+
+# The PROJECT_NAME tag is a single word (or a sequence of words surrounded
+# by quotes) that should identify the project.
+
+PROJECT_NAME = libDAI
+
+# The PROJECT_NUMBER tag can be used to enter a project or revision number.
+# This could be handy for archiving the generated documentation or
+# if some version control system is used.
+
+PROJECT_NUMBER =
+
+# The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute)
+# base path where the generated documentation will be put.
+# If a relative path is entered, it will be relative to the location
+# where doxygen was started. If left blank the current directory will be used.
+
+OUTPUT_DIRECTORY = doc
+
+# If the CREATE_SUBDIRS tag is set to YES, then doxygen will create
+# 4096 sub-directories (in 2 levels) under the output directory of each output
+# format and will distribute the generated files over these directories.
+# Enabling this option can be useful when feeding doxygen a huge amount of
+# source files, where putting all generated files in the same directory would
+# otherwise cause performance problems for the file system.
+
+CREATE_SUBDIRS = NO
+
+# The OUTPUT_LANGUAGE tag is used to specify the language in which all
+# documentation generated by doxygen is written. Doxygen will use this
+# information to generate all constant output in the proper language.
+# The default language is English, other supported languages are:
+# Brazilian, Catalan, Chinese, Chinese-Traditional, Croatian, Czech, Danish,
+# Dutch, Finnish, French, German, Greek, Hungarian, Italian, Japanese,
+# Japanese-en (Japanese with English messages), Korean, Korean-en, Norwegian,
+# Polish, Portuguese, Romanian, Russian, Serbian, Slovak, Slovene, Spanish,
+# Swedish, and Ukrainian.
+
+OUTPUT_LANGUAGE = English
+
+# This tag can be used to specify the encoding used in the generated output.
+# The encoding is not always determined by the language that is chosen,
+# but also whether or not the output is meant for Windows or non-Windows users.
+# In case there is a difference, setting the USE_WINDOWS_ENCODING tag to YES
+# forces the Windows encoding (this is the default for the Windows binary),
+# whereas setting the tag to NO uses a Unix-style encoding (the default for
+# all platforms other than Windows).
+
+USE_WINDOWS_ENCODING = NO
+
+# If the BRIEF_MEMBER_DESC tag is set to YES (the default) Doxygen will
+# include brief member descriptions after the members that are listed in
+# the file and class documentation (similar to JavaDoc).
+# Set to NO to disable this.
+
+BRIEF_MEMBER_DESC = YES
+
+# If the REPEAT_BRIEF tag is set to YES (the default) Doxygen will prepend
+# the brief description of a member or function before the detailed description.
+# Note: if both HIDE_UNDOC_MEMBERS and BRIEF_MEMBER_DESC are set to NO, the
+# brief descriptions will be completely suppressed.
+
+REPEAT_BRIEF = YES
+
+# This tag implements a quasi-intelligent brief description abbreviator
+# that is used to form the text in various listings. Each string
+# in this list, if found as the leading text of the brief description, will be
+# stripped from the text and the result after processing the whole list, is
+# used as the annotated text. Otherwise, the brief description is used as-is.
+# If left blank, the following values are used ("$name" is automatically
+# replaced with the name of the entity): "The $name class" "The $name widget"
+# "The $name file" "is" "provides" "specifies" "contains"
+# "represents" "a" "an" "the"
+
+ABBREVIATE_BRIEF =
+
+# If the ALWAYS_DETAILED_SEC and REPEAT_BRIEF tags are both set to YES then
+# Doxygen will generate a detailed section even if there is only a brief
+# description.
+
+ALWAYS_DETAILED_SEC = NO
+
+# If the INLINE_INHERITED_MEMB tag is set to YES, doxygen will show all
+# inherited members of a class in the documentation of that class as if those
+# members were ordinary class members. Constructors, destructors and assignment
+# operators of the base classes will not be shown.
+
+INLINE_INHERITED_MEMB = NO
+
+# If the FULL_PATH_NAMES tag is set to YES then Doxygen will prepend the full
+# path before files name in the file list and in the header files. If set
+# to NO the shortest path that makes the file name unique will be used.
+
+FULL_PATH_NAMES = YES
+
+# If the FULL_PATH_NAMES tag is set to YES then the STRIP_FROM_PATH tag
+# can be used to strip a user-defined part of the path. Stripping is
+# only done if one of the specified strings matches the left-hand part of
+# the path. The tag can be used to show relative paths in the file list.
+# If left blank the directory from which doxygen is run is used as the
+# path to strip.
+
+STRIP_FROM_PATH =
+
+# The STRIP_FROM_INC_PATH tag can be used to strip a user-defined part of
+# the path mentioned in the documentation of a class, which tells
+# the reader which header file to include in order to use a class.
+# If left blank only the name of the header file containing the class
+# definition is used. Otherwise one should specify the include paths that
+# are normally passed to the compiler using the -I flag.
+
+STRIP_FROM_INC_PATH =
+
+# If the SHORT_NAMES tag is set to YES, doxygen will generate much shorter
+# (but less readable) file names. This can be useful is your file systems
+# doesn't support long names like on DOS, Mac, or CD-ROM.
+
+SHORT_NAMES = NO
+
+# If the JAVADOC_AUTOBRIEF tag is set to YES then Doxygen
+# will interpret the first line (until the first dot) of a JavaDoc-style
+# comment as the brief description. If set to NO, the JavaDoc
+# comments will behave just like the Qt-style comments (thus requiring an
+# explicit @brief command for a brief description.
+
+JAVADOC_AUTOBRIEF = NO
+
+# The MULTILINE_CPP_IS_BRIEF tag can be set to YES to make Doxygen
+# treat a multi-line C++ special comment block (i.e. a block of //! or ///
+# comments) as a brief description. This used to be the default behaviour.
+# The new default is to treat a multi-line C++ comment block as a detailed
+# description. Set this tag to YES if you prefer the old behaviour instead.
+
+MULTILINE_CPP_IS_BRIEF = NO
+
+# If the DETAILS_AT_TOP tag is set to YES then Doxygen
+# will output the detailed description near the top, like JavaDoc.
+# If set to NO, the detailed description appears after the member
+# documentation.
+
+DETAILS_AT_TOP = NO
+
+# If the INHERIT_DOCS tag is set to YES (the default) then an undocumented
+# member inherits the documentation from any documented member that it
+# re-implements.
+
+INHERIT_DOCS = YES
+
+# If member grouping is used in the documentation and the DISTRIBUTE_GROUP_DOC
+# tag is set to YES, then doxygen will reuse the documentation of the first
+# member in the group (if any) for the other members of the group. By default
+# all members of a group must be documented explicitly.
+
+DISTRIBUTE_GROUP_DOC = NO
+
+# If the SEPARATE_MEMBER_PAGES tag is set to YES, then doxygen will produce
+# a new page for each member. If set to NO, the documentation of a member will
+# be part of the file/class/namespace that contains it.
+
+SEPARATE_MEMBER_PAGES = NO
+
+# The TAB_SIZE tag can be used to set the number of spaces in a tab.
+# Doxygen uses this value to replace tabs by spaces in code fragments.
+
+TAB_SIZE = 8
+
+# This tag can be used to specify a number of aliases that acts
+# as commands in the documentation. An alias has the form "name=value".
+# For example adding "sideeffect=\par Side Effects:\n" will allow you to
+# put the command \sideeffect (or @sideeffect) in the documentation, which
+# will result in a user-defined paragraph with heading "Side Effects:".
+# You can put \n's in the value part of an alias to insert newlines.
+
+ALIASES =
+
+# Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C
+# sources only. Doxygen will then generate output that is more tailored for C.
+# For instance, some of the names that are used will be different. The list
+# of all members will be omitted, etc.
+
+OPTIMIZE_OUTPUT_FOR_C = NO
+
+# Set the OPTIMIZE_OUTPUT_JAVA tag to YES if your project consists of Java sources
+# only. Doxygen will then generate output that is more tailored for Java.
+# For instance, namespaces will be presented as packages, qualified scopes
+# will look different, etc.
+
+OPTIMIZE_OUTPUT_JAVA = NO
+
+# Set the SUBGROUPING tag to YES (the default) to allow class member groups of
+# the same type (for instance a group of public functions) to be put as a
+# subgroup of that type (e.g. under the Public Functions section). Set it to
+# NO to prevent subgrouping. Alternatively, this can be done per class using
+# the \nosubgrouping command.
+
+SUBGROUPING = YES
+
+#---------------------------------------------------------------------------
+# Build related configuration options
+#---------------------------------------------------------------------------
+
+# If the EXTRACT_ALL tag is set to YES doxygen will assume all entities in
+# documentation are documented, even if no documentation was available.
+# Private class members and static file members will be hidden unless
+# the EXTRACT_PRIVATE and EXTRACT_STATIC tags are set to YES
+
+EXTRACT_ALL = YES
+
+# If the EXTRACT_PRIVATE tag is set to YES all private members of a class
+# will be included in the documentation.
+
+EXTRACT_PRIVATE = NO
+
+# If the EXTRACT_STATIC tag is set to YES all static members of a file
+# will be included in the documentation.
+
+EXTRACT_STATIC = NO
+
+# If the EXTRACT_LOCAL_CLASSES tag is set to YES classes (and structs)
+# defined locally in source files will be included in the documentation.
+# If set to NO only classes defined in header files are included.
+
+EXTRACT_LOCAL_CLASSES = YES
+
+# This flag is only useful for Objective-C code. When set to YES local
+# methods, which are defined in the implementation section but not in
+# the interface are included in the documentation.
+# If set to NO (the default) only methods in the interface are included.
+
+EXTRACT_LOCAL_METHODS = NO
+
+# If the HIDE_UNDOC_MEMBERS tag is set to YES, Doxygen will hide all
+# undocumented members of documented classes, files or namespaces.
+# If set to NO (the default) these members will be included in the
+# various overviews, but no documentation section is generated.
+# This option has no effect if EXTRACT_ALL is enabled.
+
+HIDE_UNDOC_MEMBERS = NO
+
+# If the HIDE_UNDOC_CLASSES tag is set to YES, Doxygen will hide all
+# undocumented classes that are normally visible in the class hierarchy.
+# If set to NO (the default) these classes will be included in the various
+# overviews. This option has no effect if EXTRACT_ALL is enabled.
+
+HIDE_UNDOC_CLASSES = NO
+
+# If the HIDE_FRIEND_COMPOUNDS tag is set to YES, Doxygen will hide all
+# friend (class|struct|union) declarations.
+# If set to NO (the default) these declarations will be included in the
+# documentation.
+
+HIDE_FRIEND_COMPOUNDS = NO
+
+# If the HIDE_IN_BODY_DOCS tag is set to YES, Doxygen will hide any
+# documentation blocks found inside the body of a function.
+# If set to NO (the default) these blocks will be appended to the
+# function's detailed documentation block.
+
+HIDE_IN_BODY_DOCS = NO
+
+# The INTERNAL_DOCS tag determines if documentation
+# that is typed after a \internal command is included. If the tag is set
+# to NO (the default) then the documentation will be excluded.
+# Set it to YES to include the internal documentation.
+
+INTERNAL_DOCS = NO
+
+# If the CASE_SENSE_NAMES tag is set to NO then Doxygen will only generate
+# file names in lower-case letters. If set to YES upper-case letters are also
+# allowed. This is useful if you have classes or files whose names only differ
+# in case and if your file system supports case sensitive file names. Windows
+# and Mac users are advised to set this option to NO.
+
+CASE_SENSE_NAMES = YES
+
+# If the HIDE_SCOPE_NAMES tag is set to NO (the default) then Doxygen
+# will show members with their full class and namespace scopes in the
+# documentation. If set to YES the scope will be hidden.
+
+HIDE_SCOPE_NAMES = NO
+
+# If the SHOW_INCLUDE_FILES tag is set to YES (the default) then Doxygen
+# will put a list of the files that are included by a file in the documentation
+# of that file.
+
+SHOW_INCLUDE_FILES = YES
+
+# If the INLINE_INFO tag is set to YES (the default) then a tag [inline]
+# is inserted in the documentation for inline members.
+
+INLINE_INFO = YES
+
+# If the SORT_MEMBER_DOCS tag is set to YES (the default) then doxygen
+# will sort the (detailed) documentation of file and class members
+# alphabetically by member name. If set to NO the members will appear in
+# declaration order.
+
+SORT_MEMBER_DOCS = YES
+
+# If the SORT_BRIEF_DOCS tag is set to YES then doxygen will sort the
+# brief documentation of file, namespace and class members alphabetically
+# by member name. If set to NO (the default) the members will appear in
+# declaration order.
+
+SORT_BRIEF_DOCS = NO
+
+# If the SORT_BY_SCOPE_NAME tag is set to YES, the class list will be
+# sorted by fully-qualified names, including namespaces. If set to
+# NO (the default), the class list will be sorted only by class name,
+# not including the namespace part.
+# Note: This option is not very useful if HIDE_SCOPE_NAMES is set to YES.
+# Note: This option applies only to the class list, not to the
+# alphabetical list.
+
+SORT_BY_SCOPE_NAME = NO
+
+# The GENERATE_TODOLIST tag can be used to enable (YES) or
+# disable (NO) the todo list. This list is created by putting \todo
+# commands in the documentation.
+
+GENERATE_TODOLIST = YES
+
+# The GENERATE_TESTLIST tag can be used to enable (YES) or
+# disable (NO) the test list. This list is created by putting \test
+# commands in the documentation.
+
+GENERATE_TESTLIST = YES
+
+# The GENERATE_BUGLIST tag can be used to enable (YES) or
+# disable (NO) the bug list. This list is created by putting \bug
+# commands in the documentation.
+
+GENERATE_BUGLIST = YES
+
+# The GENERATE_DEPRECATEDLIST tag can be used to enable (YES) or
+# disable (NO) the deprecated list. This list is created by putting
+# \deprecated commands in the documentation.
+
+GENERATE_DEPRECATEDLIST= YES
+
+# The ENABLED_SECTIONS tag can be used to enable conditional
+# documentation sections, marked by \if sectionname ... \endif.
+
+ENABLED_SECTIONS =
+
+# The MAX_INITIALIZER_LINES tag determines the maximum number of lines
+# the initial value of a variable or define consists of for it to appear in
+# the documentation. If the initializer consists of more lines than specified
+# here it will be hidden. Use a value of 0 to hide initializers completely.
+# The appearance of the initializer of individual variables and defines in the
+# documentation can be controlled using \showinitializer or \hideinitializer
+# command in the documentation regardless of this setting.
+
+MAX_INITIALIZER_LINES = 30
+
+# Set the SHOW_USED_FILES tag to NO to disable the list of files generated
+# at the bottom of the documentation of classes and structs. If set to YES the
+# list will mention the files that were used to generate the documentation.
+
+SHOW_USED_FILES = YES
+
+# If the sources in your project are distributed over multiple directories
+# then setting the SHOW_DIRECTORIES tag to YES will show the directory hierarchy
+# in the documentation.
+
+SHOW_DIRECTORIES = YES
+
+# The FILE_VERSION_FILTER tag can be used to specify a program or script that
+# doxygen should invoke to get the current version for each file (typically from the
+# version control system). Doxygen will invoke the program by executing (via
+# popen()) the command <command> <input-file>, where <command> is the value of
+# the FILE_VERSION_FILTER tag, and <input-file> is the name of an input file
+# provided by doxygen. Whatever the progam writes to standard output
+# is used as the file version. See the manual for examples.
+
+FILE_VERSION_FILTER =
+
+#---------------------------------------------------------------------------
+# configuration options related to warning and progress messages
+#---------------------------------------------------------------------------
+
+# The QUIET tag can be used to turn on/off the messages that are generated
+# by doxygen. Possible values are YES and NO. If left blank NO is used.
+
+QUIET = NO
+
+# The WARNINGS tag can be used to turn on/off the warning messages that are
+# generated by doxygen. Possible values are YES and NO. If left blank
+# NO is used.
+
+WARNINGS = YES
+
+# If WARN_IF_UNDOCUMENTED is set to YES, then doxygen will generate warnings
+# for undocumented members. If EXTRACT_ALL is set to YES then this flag will
+# automatically be disabled.
+
+WARN_IF_UNDOCUMENTED = YES
+
+# If WARN_IF_DOC_ERROR is set to YES, doxygen will generate warnings for
+# potential errors in the documentation, such as not documenting some
+# parameters in a documented function, or documenting parameters that
+# don't exist or using markup commands wrongly.
+
+WARN_IF_DOC_ERROR = YES
+
+# This WARN_NO_PARAMDOC option can be abled to get warnings for
+# functions that are documented, but have no documentation for their parameters
+# or return value. If set to NO (the default) doxygen will only warn about
+# wrong or incomplete parameter documentation, but not about the absence of
+# documentation.
+
+WARN_NO_PARAMDOC = NO
+
+# The WARN_FORMAT tag determines the format of the warning messages that
+# doxygen can produce. The string should contain the $file, $line, and $text
+# tags, which will be replaced by the file and line number from which the
+# warning originated and the warning text. Optionally the format may contain
+# $version, which will be replaced by the version of the file (if it could
+# be obtained via FILE_VERSION_FILTER)
+
+WARN_FORMAT = "$file:$line: $text"
+
+# The WARN_LOGFILE tag can be used to specify a file to which warning
+# and error messages should be written. If left blank the output is written
+# to stderr.
+
+WARN_LOGFILE =
+
+#---------------------------------------------------------------------------
+# configuration options related to the input files
+#---------------------------------------------------------------------------
+
+# The INPUT tag can be used to specify the files and/or directories that contain
+# documented source files. You may enter file names like "myfile.cpp" or
+# directories like "/usr/src/myproject". Separate the files or directories
+# with spaces.
+
+INPUT =
+
+# If the value of the INPUT tag contains directories, you can use the
+# FILE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp
+# and *.h) to filter out the source-files in the directories. If left
+# blank the following patterns are tested:
+# *.c *.cc *.cxx *.cpp *.c++ *.java *.ii *.ixx *.ipp *.i++ *.inl *.h *.hh *.hxx
+# *.hpp *.h++ *.idl *.odl *.cs *.php *.php3 *.inc *.m *.mm
+
+FILE_PATTERNS =
+
+# The RECURSIVE tag can be used to turn specify whether or not subdirectories
+# should be searched for input files as well. Possible values are YES and NO.
+# If left blank NO is used.
+
+RECURSIVE = NO
+
+# The EXCLUDE tag can be used to specify files and/or directories that should
+# excluded from the INPUT source files. This way you can easily exclude a
+# subdirectory from a directory tree whose root is specified with the INPUT tag.
+
+EXCLUDE =
+
+# The EXCLUDE_SYMLINKS tag can be used select whether or not files or
+# directories that are symbolic links (a Unix filesystem feature) are excluded
+# from the input.
+
+EXCLUDE_SYMLINKS = NO
+
+# If the value of the INPUT tag contains directories, you can use the
+# EXCLUDE_PATTERNS tag to specify one or more wildcard patterns to exclude
+# certain files from those directories.
+
+EXCLUDE_PATTERNS =
+
+# The EXAMPLE_PATH tag can be used to specify one or more files or
+# directories that contain example code fragments that are included (see
+# the \include command).
+
+EXAMPLE_PATH =
+
+# If the value of the EXAMPLE_PATH tag contains directories, you can use the
+# EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp
+# and *.h) to filter out the source-files in the directories. If left
+# blank all files are included.
+
+EXAMPLE_PATTERNS =
+
+# If the EXAMPLE_RECURSIVE tag is set to YES then subdirectories will be
+# searched for input files to be used with the \include or \dontinclude
+# commands irrespective of the value of the RECURSIVE tag.
+# Possible values are YES and NO. If left blank NO is used.
+
+EXAMPLE_RECURSIVE = NO
+
+# The IMAGE_PATH tag can be used to specify one or more files or
+# directories that contain image that are included in the documentation (see
+# the \image command).
+
+IMAGE_PATH =
+
+# The INPUT_FILTER tag can be used to specify a program that doxygen should
+# invoke to filter for each input file. Doxygen will invoke the filter program
+# by executing (via popen()) the command <filter> <input-file>, where <filter>
+# is the value of the INPUT_FILTER tag, and <input-file> is the name of an
+# input file. Doxygen will then use the output that the filter program writes
+# to standard output. If FILTER_PATTERNS is specified, this tag will be
+# ignored.
+
+INPUT_FILTER =
+
+# The FILTER_PATTERNS tag can be used to specify filters on a per file pattern
+# basis. Doxygen will compare the file name with each pattern and apply the
+# filter if there is a match. The filters are a list of the form:
+# pattern=filter (like *.cpp=my_cpp_filter). See INPUT_FILTER for further
+# info on how filters are used. If FILTER_PATTERNS is empty, INPUT_FILTER
+# is applied to all files.
+
+FILTER_PATTERNS =
+
+# If the FILTER_SOURCE_FILES tag is set to YES, the input filter (if set using
+# INPUT_FILTER) will be used to filter the input files when producing source
+# files to browse (i.e. when SOURCE_BROWSER is set to YES).
+
+FILTER_SOURCE_FILES = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to source browsing
+#---------------------------------------------------------------------------
+
+# If the SOURCE_BROWSER tag is set to YES then a list of source files will
+# be generated. Documented entities will be cross-referenced with these sources.
+# Note: To get rid of all source code in the generated output, make sure also
+# VERBATIM_HEADERS is set to NO.
+
+SOURCE_BROWSER = NO
+
+# Setting the INLINE_SOURCES tag to YES will include the body
+# of functions and classes directly in the documentation.
+
+INLINE_SOURCES = NO
+
+# Setting the STRIP_CODE_COMMENTS tag to YES (the default) will instruct
+# doxygen to hide any special comment blocks from generated source code
+# fragments. Normal C and C++ comments will always remain visible.
+
+STRIP_CODE_COMMENTS = YES
+
+# If the REFERENCED_BY_RELATION tag is set to YES (the default)
+# then for each documented function all documented
+# functions referencing it will be listed.
+
+REFERENCED_BY_RELATION = YES
+
+# If the REFERENCES_RELATION tag is set to YES (the default)
+# then for each documented function all documented entities
+# called/used by that function will be listed.
+
+REFERENCES_RELATION = YES
+
+# If the VERBATIM_HEADERS tag is set to YES (the default) then Doxygen
+# will generate a verbatim copy of the header file for each class for
+# which an include is specified. Set to NO to disable this.
+
+VERBATIM_HEADERS = YES
+
+#---------------------------------------------------------------------------
+# configuration options related to the alphabetical class index
+#---------------------------------------------------------------------------
+
+# If the ALPHABETICAL_INDEX tag is set to YES, an alphabetical index
+# of all compounds will be generated. Enable this if the project
+# contains a lot of classes, structs, unions or interfaces.
+
+ALPHABETICAL_INDEX = NO
+
+# If the alphabetical index is enabled (see ALPHABETICAL_INDEX) then
+# the COLS_IN_ALPHA_INDEX tag can be used to specify the number of columns
+# in which this list will be split (can be a number in the range [1..20])
+
+COLS_IN_ALPHA_INDEX = 5
+
+# In case all classes in a project start with a common prefix, all
+# classes will be put under the same header in the alphabetical index.
+# The IGNORE_PREFIX tag can be used to specify one or more prefixes that
+# should be ignored while generating the index headers.
+
+IGNORE_PREFIX =
+
+#---------------------------------------------------------------------------
+# configuration options related to the HTML output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_HTML tag is set to YES (the default) Doxygen will
+# generate HTML output.
+
+GENERATE_HTML = YES
+
+# The HTML_OUTPUT tag is used to specify where the HTML docs will be put.
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be
+# put in front of it. If left blank `html' will be used as the default path.
+
+HTML_OUTPUT = html
+
+# The HTML_FILE_EXTENSION tag can be used to specify the file extension for
+# each generated HTML page (for example: .htm,.php,.asp). If it is left blank
+# doxygen will generate files with .html extension.
+
+HTML_FILE_EXTENSION = .html
+
+# The HTML_HEADER tag can be used to specify a personal HTML header for
+# each generated HTML page. If it is left blank doxygen will generate a
+# standard header.
+
+HTML_HEADER =
+
+# The HTML_FOOTER tag can be used to specify a personal HTML footer for
+# each generated HTML page. If it is left blank doxygen will generate a
+# standard footer.
+
+HTML_FOOTER =
+
+# The HTML_STYLESHEET tag can be used to specify a user-defined cascading
+# style sheet that is used by each HTML page. It can be used to
+# fine-tune the look of the HTML output. If the tag is left blank doxygen
+# will generate a default style sheet. Note that doxygen will try to copy
+# the style sheet file to the HTML output directory, so don't put your own
+# stylesheet in the HTML output directory as well, or it will be erased!
+
+HTML_STYLESHEET =
+
+# If the HTML_ALIGN_MEMBERS tag is set to YES, the members of classes,
+# files or namespaces will be aligned in HTML using tables. If set to
+# NO a bullet list will be used.
+
+HTML_ALIGN_MEMBERS = YES
+
+# If the GENERATE_HTMLHELP tag is set to YES, additional index files
+# will be generated that can be used as input for tools like the
+# Microsoft HTML help workshop to generate a compressed HTML help file (.chm)
+# of the generated HTML documentation.
+
+GENERATE_HTMLHELP = NO
+
+# If the GENERATE_HTMLHELP tag is set to YES, the CHM_FILE tag can
+# be used to specify the file name of the resulting .chm file. You
+# can add a path in front of the file if the result should not be
+# written to the html output directory.
+
+CHM_FILE =
+
+# If the GENERATE_HTMLHELP tag is set to YES, the HHC_LOCATION tag can
+# be used to specify the location (absolute path including file name) of
+# the HTML help compiler (hhc.exe). If non-empty doxygen will try to run
+# the HTML help compiler on the generated index.hhp.
+
+HHC_LOCATION =
+
+# If the GENERATE_HTMLHELP tag is set to YES, the GENERATE_CHI flag
+# controls if a separate .chi index file is generated (YES) or that
+# it should be included in the master .chm file (NO).
+
+GENERATE_CHI = NO
+
+# If the GENERATE_HTMLHELP tag is set to YES, the BINARY_TOC flag
+# controls whether a binary table of contents is generated (YES) or a
+# normal table of contents (NO) in the .chm file.
+
+BINARY_TOC = NO
+
+# The TOC_EXPAND flag can be set to YES to add extra items for group members
+# to the contents of the HTML help documentation and to the tree view.
+
+TOC_EXPAND = NO
+
+# The DISABLE_INDEX tag can be used to turn on/off the condensed index at
+# top of each HTML page. The value NO (the default) enables the index and
+# the value YES disables it.
+
+DISABLE_INDEX = NO
+
+# This tag can be used to set the number of enum values (range [1..20])
+# that doxygen will group on one line in the generated HTML documentation.
+
+ENUM_VALUES_PER_LINE = 4
+
+# If the GENERATE_TREEVIEW tag is set to YES, a side panel will be
+# generated containing a tree-like index structure (just like the one that
+# is generated for HTML Help). For this to work a browser that supports
+# JavaScript, DHTML, CSS and frames is required (for instance Mozilla 1.0+,
+# Netscape 6.0+, Internet explorer 5.0+, or Konqueror). Windows users are
+# probably better off using the HTML help feature.
+
+GENERATE_TREEVIEW = NO
+
+# If the treeview is enabled (see GENERATE_TREEVIEW) then this tag can be
+# used to set the initial width (in pixels) of the frame in which the tree
+# is shown.
+
+TREEVIEW_WIDTH = 250
+
+#---------------------------------------------------------------------------
+# configuration options related to the LaTeX output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_LATEX tag is set to YES (the default) Doxygen will
+# generate Latex output.
+
+GENERATE_LATEX = YES
+
+# The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put.
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be
+# put in front of it. If left blank `latex' will be used as the default path.
+
+LATEX_OUTPUT = latex
+
+# The LATEX_CMD_NAME tag can be used to specify the LaTeX command name to be
+# invoked. If left blank `latex' will be used as the default command name.
+
+LATEX_CMD_NAME = latex
+
+# The MAKEINDEX_CMD_NAME tag can be used to specify the command name to
+# generate index for LaTeX. If left blank `makeindex' will be used as the
+# default command name.
+
+MAKEINDEX_CMD_NAME = makeindex
+
+# If the COMPACT_LATEX tag is set to YES Doxygen generates more compact
+# LaTeX documents. This may be useful for small projects and may help to
+# save some trees in general.
+
+COMPACT_LATEX = NO
+
+# The PAPER_TYPE tag can be used to set the paper type that is used
+# by the printer. Possible values are: a4, a4wide, letter, legal and
+# executive. If left blank a4wide will be used.
+
+PAPER_TYPE = a4wide
+
+# The EXTRA_PACKAGES tag can be to specify one or more names of LaTeX
+# packages that should be included in the LaTeX output.
+
+EXTRA_PACKAGES =
+
+# The LATEX_HEADER tag can be used to specify a personal LaTeX header for
+# the generated latex document. The header should contain everything until
+# the first chapter. If it is left blank doxygen will generate a
+# standard header. Notice: only use this tag if you know what you are doing!
+
+LATEX_HEADER =
+
+# If the PDF_HYPERLINKS tag is set to YES, the LaTeX that is generated
+# is prepared for conversion to pdf (using ps2pdf). The pdf file will
+# contain links (just like the HTML output) instead of page references
+# This makes the output suitable for online browsing using a pdf viewer.
+
+PDF_HYPERLINKS = NO
+
+# If the USE_PDFLATEX tag is set to YES, pdflatex will be used instead of
+# plain latex in the generated Makefile. Set this option to YES to get a
+# higher quality PDF documentation.
+
+USE_PDFLATEX = NO
+
+# If the LATEX_BATCHMODE tag is set to YES, doxygen will add the \\batchmode.
+# command to the generated LaTeX files. This will instruct LaTeX to keep
+# running if errors occur, instead of asking the user for help.
+# This option is also used when generating formulas in HTML.
+
+LATEX_BATCHMODE = NO
+
+# If LATEX_HIDE_INDICES is set to YES then doxygen will not
+# include the index chapters (such as File Index, Compound Index, etc.)
+# in the output.
+
+LATEX_HIDE_INDICES = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to the RTF output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_RTF tag is set to YES Doxygen will generate RTF output
+# The RTF output is optimized for Word 97 and may not look very pretty with
+# other RTF readers or editors.
+
+GENERATE_RTF = NO
+
+# The RTF_OUTPUT tag is used to specify where the RTF docs will be put.
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be
+# put in front of it. If left blank `rtf' will be used as the default path.
+
+RTF_OUTPUT = rtf
+
+# If the COMPACT_RTF tag is set to YES Doxygen generates more compact
+# RTF documents. This may be useful for small projects and may help to
+# save some trees in general.
+
+COMPACT_RTF = NO
+
+# If the RTF_HYPERLINKS tag is set to YES, the RTF that is generated
+# will contain hyperlink fields. The RTF file will
+# contain links (just like the HTML output) instead of page references.
+# This makes the output suitable for online browsing using WORD or other
+# programs which support those fields.
+# Note: wordpad (write) and others do not support links.
+
+RTF_HYPERLINKS = NO
+
+# Load stylesheet definitions from file. Syntax is similar to doxygen's
+# config file, i.e. a series of assignments. You only have to provide
+# replacements, missing definitions are set to their default value.
+
+RTF_STYLESHEET_FILE =
+
+# Set optional variables used in the generation of an rtf document.
+# Syntax is similar to doxygen's config file.
+
+RTF_EXTENSIONS_FILE =
+
+#---------------------------------------------------------------------------
+# configuration options related to the man page output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_MAN tag is set to YES (the default) Doxygen will
+# generate man pages
+
+GENERATE_MAN = NO
+
+# The MAN_OUTPUT tag is used to specify where the man pages will be put.
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be
+# put in front of it. If left blank `man' will be used as the default path.
+
+MAN_OUTPUT = man
+
+# The MAN_EXTENSION tag determines the extension that is added to
+# the generated man pages (default is the subroutine's section .3)
+
+MAN_EXTENSION = .3
+
+# If the MAN_LINKS tag is set to YES and Doxygen generates man output,
+# then it will generate one additional man file for each entity
+# documented in the real man page(s). These additional files
+# only source the real man page, but without them the man command
+# would be unable to find the correct page. The default is NO.
+
+MAN_LINKS = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to the XML output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_XML tag is set to YES Doxygen will
+# generate an XML file that captures the structure of
+# the code including all documentation.
+
+GENERATE_XML = NO
+
+# The XML_OUTPUT tag is used to specify where the XML pages will be put.
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be
+# put in front of it. If left blank `xml' will be used as the default path.
+
+XML_OUTPUT = xml
+
+# The XML_SCHEMA tag can be used to specify an XML schema,
+# which can be used by a validating XML parser to check the
+# syntax of the XML files.
+
+XML_SCHEMA =
+
+# The XML_DTD tag can be used to specify an XML DTD,
+# which can be used by a validating XML parser to check the
+# syntax of the XML files.
+
+XML_DTD =
+
+# If the XML_PROGRAMLISTING tag is set to YES Doxygen will
+# dump the program listings (including syntax highlighting
+# and cross-referencing information) to the XML output. Note that
+# enabling this will significantly increase the size of the XML output.
+
+XML_PROGRAMLISTING = YES
+
+#---------------------------------------------------------------------------
+# configuration options for the AutoGen Definitions output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_AUTOGEN_DEF tag is set to YES Doxygen will
+# generate an AutoGen Definitions (see autogen.sf.net) file
+# that captures the structure of the code including all
+# documentation. Note that this feature is still experimental
+# and incomplete at the moment.
+
+GENERATE_AUTOGEN_DEF = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to the Perl module output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_PERLMOD tag is set to YES Doxygen will
+# generate a Perl module file that captures the structure of
+# the code including all documentation. Note that this
+# feature is still experimental and incomplete at the
+# moment.
+
+GENERATE_PERLMOD = NO
+
+# If the PERLMOD_LATEX tag is set to YES Doxygen will generate
+# the necessary Makefile rules, Perl scripts and LaTeX code to be able
+# to generate PDF and DVI output from the Perl module output.
+
+PERLMOD_LATEX = NO
+
+# If the PERLMOD_PRETTY tag is set to YES the Perl module output will be
+# nicely formatted so it can be parsed by a human reader. This is useful
+# if you want to understand what is going on. On the other hand, if this
+# tag is set to NO the size of the Perl module output will be much smaller
+# and Perl will parse it just the same.
+
+PERLMOD_PRETTY = YES
+
+# The names of the make variables in the generated doxyrules.make file
+# are prefixed with the string contained in PERLMOD_MAKEVAR_PREFIX.
+# This is useful so different doxyrules.make files included by the same
+# Makefile don't overwrite each other's variables.
+
+PERLMOD_MAKEVAR_PREFIX =
+
+#---------------------------------------------------------------------------
+# Configuration options related to the preprocessor
+#---------------------------------------------------------------------------
+
+# If the ENABLE_PREPROCESSING tag is set to YES (the default) Doxygen will
+# evaluate all C-preprocessor directives found in the sources and include
+# files.
+
+ENABLE_PREPROCESSING = YES
+
+# If the MACRO_EXPANSION tag is set to YES Doxygen will expand all macro
+# names in the source code. If set to NO (the default) only conditional
+# compilation will be performed. Macro expansion can be done in a controlled
+# way by setting EXPAND_ONLY_PREDEF to YES.
+
+MACRO_EXPANSION = NO
+
+# If the EXPAND_ONLY_PREDEF and MACRO_EXPANSION tags are both set to YES
+# then the macro expansion is limited to the macros specified with the
+# PREDEFINED and EXPAND_AS_PREDEFINED tags.
+
+EXPAND_ONLY_PREDEF = NO
+
+# If the SEARCH_INCLUDES tag is set to YES (the default) the includes files
+# in the INCLUDE_PATH (see below) will be search if a #include is found.
+
+SEARCH_INCLUDES = YES
+
+# The INCLUDE_PATH tag can be used to specify one or more directories that
+# contain include files that are not input files but should be processed by
+# the preprocessor.
+
+INCLUDE_PATH =
+
+# You can use the INCLUDE_FILE_PATTERNS tag to specify one or more wildcard
+# patterns (like *.h and *.hpp) to filter out the header-files in the
+# directories. If left blank, the patterns specified with FILE_PATTERNS will
+# be used.
+
+INCLUDE_FILE_PATTERNS =
+
+# The PREDEFINED tag can be used to specify one or more macro names that
+# are defined before the preprocessor is started (similar to the -D option of
+# gcc). The argument of the tag is a list of macros of the form: name
+# or name=definition (no spaces). If the definition and the = are
+# omitted =1 is assumed. To prevent a macro definition from being
+# undefined via #undef or recursively expanded use the := operator
+# instead of the = operator.
+
+PREDEFINED =
+
+# If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then
+# this tag can be used to specify a list of macro names that should be expanded.
+# The macro definition that is found in the sources will be used.
+# Use the PREDEFINED tag if you want to use a different macro definition.
+
+EXPAND_AS_DEFINED =
+
+# If the SKIP_FUNCTION_MACROS tag is set to YES (the default) then
+# doxygen's preprocessor will remove all function-like macros that are alone
+# on a line, have an all uppercase name, and do not end with a semicolon. Such
+# function macros are typically used for boiler-plate code, and will confuse
+# the parser if not removed.
+
+SKIP_FUNCTION_MACROS = YES
+
+#---------------------------------------------------------------------------
+# Configuration::additions related to external references
+#---------------------------------------------------------------------------
+
+# The TAGFILES option can be used to specify one or more tagfiles.
+# Optionally an initial location of the external documentation
+# can be added for each tagfile. The format of a tag file without
+# this location is as follows:
+# TAGFILES = file1 file2 ...
+# Adding location for the tag files is done as follows:
+# TAGFILES = file1=loc1 "file2 = loc2" ...
+# where "loc1" and "loc2" can be relative or absolute paths or
+# URLs. If a location is present for each tag, the installdox tool
+# does not have to be run to correct the links.
+# Note that each tag file must have a unique name
+# (where the name does NOT include the path)
+# If a tag file is not located in the directory in which doxygen
+# is run, you must also specify the path to the tagfile here.
+
+TAGFILES =
+
+# When a file name is specified after GENERATE_TAGFILE, doxygen will create
+# a tag file that is based on the input files it reads.
+
+GENERATE_TAGFILE =
+
+# If the ALLEXTERNALS tag is set to YES all external classes will be listed
+# in the class index. If set to NO only the inherited external classes
+# will be listed.
+
+ALLEXTERNALS = NO
+
+# If the EXTERNAL_GROUPS tag is set to YES all external groups will be listed
+# in the modules index. If set to NO, only the current project's groups will
+# be listed.
+
+EXTERNAL_GROUPS = YES
+
+# The PERL_PATH should be the absolute path and name of the perl script
+# interpreter (i.e. the result of `which perl').
+
+PERL_PATH = /usr/bin/perl
+
+#---------------------------------------------------------------------------
+# Configuration options related to the dot tool
+#---------------------------------------------------------------------------
+
+# If the CLASS_DIAGRAMS tag is set to YES (the default) Doxygen will
+# generate a inheritance diagram (in HTML, RTF and LaTeX) for classes with base
+# or super classes. Setting the tag to NO turns the diagrams off. Note that
+# this option is superseded by the HAVE_DOT option below. This is only a
+# fallback. It is recommended to install and use dot, since it yields more
+# powerful graphs.
+
+CLASS_DIAGRAMS = YES
+
+# If set to YES, the inheritance and collaboration graphs will hide
+# inheritance and usage relations if the target is undocumented
+# or is not a class.
+
+HIDE_UNDOC_RELATIONS = YES
+
+# If you set the HAVE_DOT tag to YES then doxygen will assume the dot tool is
+# available from the path. This tool is part of Graphviz, a graph visualization
+# toolkit from AT&T and Lucent Bell Labs. The other options in this section
+# have no effect if this option is set to NO (the default)
+
+HAVE_DOT = NO
+
+# If the CLASS_GRAPH and HAVE_DOT tags are set to YES then doxygen
+# will generate a graph for each documented class showing the direct and
+# indirect inheritance relations. Setting this tag to YES will force the
+# the CLASS_DIAGRAMS tag to NO.
+
+CLASS_GRAPH = YES
+
+# If the COLLABORATION_GRAPH and HAVE_DOT tags are set to YES then doxygen
+# will generate a graph for each documented class showing the direct and
+# indirect implementation dependencies (inheritance, containment, and
+# class references variables) of the class with other documented classes.
+
+COLLABORATION_GRAPH = YES
+
+# If the GROUP_GRAPHS and HAVE_DOT tags are set to YES then doxygen
+# will generate a graph for groups, showing the direct groups dependencies
+
+GROUP_GRAPHS = YES
+
+# If the UML_LOOK tag is set to YES doxygen will generate inheritance and
+# collaboration diagrams in a style similar to the OMG's Unified Modeling
+# Language.
+
+UML_LOOK = NO
+
+# If set to YES, the inheritance and collaboration graphs will show the
+# relations between templates and their instances.
+
+TEMPLATE_RELATIONS = NO
+
+# If the ENABLE_PREPROCESSING, SEARCH_INCLUDES, INCLUDE_GRAPH, and HAVE_DOT
+# tags are set to YES then doxygen will generate a graph for each documented
+# file showing the direct and indirect include dependencies of the file with
+# other documented files.
+
+INCLUDE_GRAPH = YES
+
+# If the ENABLE_PREPROCESSING, SEARCH_INCLUDES, INCLUDED_BY_GRAPH, and
+# HAVE_DOT tags are set to YES then doxygen will generate a graph for each
+# documented header file showing the documented files that directly or
+# indirectly include this file.
+
+INCLUDED_BY_GRAPH = YES
+
+# If the CALL_GRAPH and HAVE_DOT tags are set to YES then doxygen will
+# generate a call dependency graph for every global function or class method.
+# Note that enabling this option will significantly increase the time of a run.
+# So in most cases it will be better to enable call graphs for selected
+# functions only using the \callgraph command.
+
+CALL_GRAPH = NO
+
+# If the GRAPHICAL_HIERARCHY and HAVE_DOT tags are set to YES then doxygen
+# will graphical hierarchy of all classes instead of a textual one.
+
+GRAPHICAL_HIERARCHY = YES
+
+# If the DIRECTORY_GRAPH, SHOW_DIRECTORIES and HAVE_DOT tags are set to YES
+# then doxygen will show the dependencies a directory has on other directories
+# in a graphical way. The dependency relations are determined by the #include
+# relations between the files in the directories.
+
+DIRECTORY_GRAPH = YES
+
+# The DOT_IMAGE_FORMAT tag can be used to set the image format of the images
+# generated by dot. Possible values are png, jpg, or gif
+# If left blank png will be used.
+
+DOT_IMAGE_FORMAT = png
+
+# The tag DOT_PATH can be used to specify the path where the dot tool can be
+# found. If left blank, it is assumed the dot tool can be found in the path.
+
+DOT_PATH =
+
+# The DOTFILE_DIRS tag can be used to specify one or more directories that
+# contain dot files that are included in the documentation (see the
+# \dotfile command).
+
+DOTFILE_DIRS =
+
+# The MAX_DOT_GRAPH_WIDTH tag can be used to set the maximum allowed width
+# (in pixels) of the graphs generated by dot. If a graph becomes larger than
+# this value, doxygen will try to truncate the graph, so that it fits within
+# the specified constraint. Beware that most browsers cannot cope with very
+# large images.
+
+MAX_DOT_GRAPH_WIDTH = 1024
+
+# The MAX_DOT_GRAPH_HEIGHT tag can be used to set the maximum allows height
+# (in pixels) of the graphs generated by dot. If a graph becomes larger than
+# this value, doxygen will try to truncate the graph, so that it fits within
+# the specified constraint. Beware that most browsers cannot cope with very
+# large images.
+
+MAX_DOT_GRAPH_HEIGHT = 1024
+
+# The MAX_DOT_GRAPH_DEPTH tag can be used to set the maximum depth of the
+# graphs generated by dot. A depth value of 3 means that only nodes reachable
+# from the root by following a path via at most 3 edges will be shown. Nodes
+# that lay further from the root node will be omitted. Note that setting this
+# option to 1 or 2 may greatly reduce the computation time needed for large
+# code bases. Also note that a graph may be further truncated if the graph's
+# image dimensions are not sufficient to fit the graph (see MAX_DOT_GRAPH_WIDTH
+# and MAX_DOT_GRAPH_HEIGHT). If 0 is used for the depth value (the default),
+# the graph is not depth-constrained.
+
+MAX_DOT_GRAPH_DEPTH = 0
+
+# Set the DOT_TRANSPARENT tag to YES to generate images with a transparent
+# background. This is disabled by default, which results in a white background.
+# Warning: Depending on the platform used, enabling this option may lead to
+# badly anti-aliased labels on the edges of a graph (i.e. they become hard to
+# read).
+
+DOT_TRANSPARENT = NO
+
+# Set the DOT_MULTI_TARGETS tag to YES allow dot to generate multiple output
+# files in one run (i.e. multiple -o and -T options on the command line). This
+# makes dot run faster, but since only newer versions of dot (>1.8.10)
+# support this, this feature is disabled by default.
+
+DOT_MULTI_TARGETS = NO
+
+# If the GENERATE_LEGEND tag is set to YES (the default) Doxygen will
+# generate a legend page explaining the meaning of the various boxes and
+# arrows in the dot generated graphs.
+
+GENERATE_LEGEND = YES
+
+# If the DOT_CLEANUP tag is set to YES (the default) Doxygen will
+# remove the intermediate dot files that are used to generate
+# the various graphs.
+
+DOT_CLEANUP = YES
+
+#---------------------------------------------------------------------------
+# Configuration::additions related to the search engine
+#---------------------------------------------------------------------------
+
+# The SEARCHENGINE tag specifies whether or not a search engine should be
+# used. If set to NO the values of all tags below this one will be ignored.
+
+SEARCHENGINE = NO
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __ENUM_H
+#define __ENUM_H
+
+
+#include <cstring>
+#include <iostream>
+
+
+// C++ enums are too limited for my purposes. This defines wrapper classes
+// that provide much more functionality than a simple enum. The only
+// disadvantage is that one wrapper class needs to be written for each
+// number of values an enum can take... a better solution is needed.
+
+
+#define ENUM2(x,a,b) class x {\
+ public:\
+ enum value {a, b};\
+\
+ x() : v(a) {}\
+\
+ x(value w) : v(w) {}\
+\
+ x(char const *w) {\
+ static char const* labels[] = {#a, #b};\
+ size_t i = 0;\
+ for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
+ if( strcmp( w, labels[i] ) == 0 ) {\
+ v = (value)i;\
+ break;\
+ }\
+ if( i == sizeof(labels) / sizeof(char const *) )\
+ throw "Unknown " #x " value";\
+ }\
+\
+ operator value () const { return v; }\
+\
+ operator size_t () const { return (size_t)v; }\
+\
+ operator char const* () const {\
+ static char const* labels[] = {#a, #b};\
+ return labels[v];\
+ }\
+\
+ friend std::istream& operator >> (std::istream& is, x& y) {\
+ std::string s;\
+ is >> s;\
+ y = x(s.c_str());\
+ return is;\
+ }\
+\
+ friend std::ostream& operator << (std::ostream& os, const x& y) {\
+ os << (const char *)y;\
+ return os;\
+ }\
+\
+ private:\
+ value v;\
+};
+
+
+#define ENUM3(x,a,b,c) class x {\
+ public:\
+ enum value {a, b, c};\
+\
+ x() : v(a) {}\
+\
+ x(value w) : v(w) {}\
+\
+ x(char const *w) {\
+ static char const* labels[] = {#a, #b, #c};\
+ size_t i = 0;\
+ for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
+ if( strcmp( w, labels[i] ) == 0 ) {\
+ v = (value)i;\
+ break;\
+ }\
+ if( i == sizeof(labels) / sizeof(char const *) )\
+ throw "Unknown " #x " value";\
+ }\
+\
+ operator value () const { return v; }\
+\
+ operator size_t () const { return (size_t)v; }\
+\
+ operator char const* () const {\
+ static char const* labels[] = {#a, #b, #c};\
+ return labels[v];\
+ }\
+\
+ friend std::istream& operator >> (std::istream& is, x& y) {\
+ std::string s;\
+ is >> s;\
+ y = x(s.c_str());\
+ return is;\
+ }\
+\
+ friend std::ostream& operator << (std::ostream& os, const x& y) {\
+ os << (const char *)y;\
+ return os;\
+ }\
+\
+ private:\
+ value v;\
+};
+
+
+#define ENUM4(x,a,b,c,d) class x {\
+ public:\
+ enum value {a, b, c, d};\
+\
+ x() : v(a) {}\
+\
+ x(value w) : v(w) {}\
+\
+ x(char const *w) {\
+ static char const* labels[] = {#a, #b, #c, #d};\
+ size_t i = 0;\
+ for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
+ if( strcmp( w, labels[i] ) == 0 ) {\
+ v = (value)i;\
+ break;\
+ }\
+ if( i == sizeof(labels) / sizeof(char const *) )\
+ throw "Unknown " #x " value";\
+ }\
+\
+ operator value () const { return v; }\
+\
+ operator size_t () const { return (size_t)v; }\
+\
+ operator char const* () const {\
+ static char const* labels[] = {#a, #b, #c, #d};\
+ return labels[v];\
+ }\
+\
+ friend std::istream& operator >> (std::istream& is, x& y) {\
+ std::string s;\
+ is >> s;\
+ y = x(s.c_str());\
+ return is;\
+ }\
+\
+ friend std::ostream& operator << (std::ostream& os, const x& y) {\
+ os << (const char *)y;\
+ return os;\
+ }\
+\
+ private:\
+ value v;\
+};
+
+
+#define ENUM5(x,a,b,c,d,e) class x {\
+ public:\
+ enum value {a, b, c, d, e};\
+\
+ x() : v(a) {}\
+\
+ x(value w) : v(w) {}\
+\
+ x(char const *w) {\
+ static char const* labels[] = {#a, #b, #c, #d, #e};\
+ size_t i = 0;\
+ for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
+ if( strcmp( w, labels[i] ) == 0 ) {\
+ v = (value)i;\
+ break;\
+ }\
+ if( i == sizeof(labels) / sizeof(char const *) )\
+ throw "Unknown " #x " value";\
+ }\
+\
+ operator value () const { return v; }\
+\
+ operator size_t () const { return (size_t)v; }\
+\
+ operator char const* () const {\
+ static char const* labels[] = {#a, #b, #c, #d, #e};\
+ return labels[v];\
+ }\
+\
+ friend std::istream& operator >> (std::istream& is, x& y) {\
+ std::string s;\
+ is >> s;\
+ y = x(s.c_str());\
+ return is;\
+ }\
+\
+ friend std::ostream& operator << (std::ostream& os, const x& y) {\
+ os << (const char *)y;\
+ return os;\
+ }\
+\
+ private:\
+ value v;\
+};
+
+
+#define ENUM6(x,a,b,c,d,e,f) class x {\
+ public:\
+ enum value {a, b, c, d, e, f};\
+\
+ x() : v(a) {}\
+\
+ x(value w) : v(w) {}\
+\
+ x(char const *w) {\
+ static char const* labels[] = {#a, #b, #c, #d, #e, #f};\
+ size_t i = 0;\
+ for( ; i < sizeof(labels) / sizeof(char const *); i++ )\
+ if( strcmp( w, labels[i] ) == 0 ) {\
+ v = (value)i;\
+ break;\
+ }\
+ if( i == sizeof(labels) / sizeof(char const *) )\
+ throw "Unknown " #x " value";\
+ }\
+\
+ operator value () const { return v; }\
+\
+ operator size_t () const { return (size_t)v; }\
+\
+ operator char const* () const {\
+ static char const* labels[] = {#a, #b, #c, #d, #e, #f};\
+ return labels[v];\
+ }\
+\
+ friend std::istream& operator >> (std::istream& is, x& y) {\
+ std::string s;\
+ is >> s;\
+ y = x(s.c_str());\
+ return is;\
+ }\
+\
+ friend std::ostream& operator << (std::ostream& os, const x& y) {\
+ os << (const char *)y;\
+ return os;\
+ }\
+\
+ private:\
+ value v;\
+};
+
+
+#endif
--- /dev/null
+/* Copyright 2006-2008 Joris Mooij
+ jorism@jorismooij.nl
+
+ This file is part of AI.
+
+ AI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ AI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with AI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#include <iostream>
+#include "alldai.h"
+
+
+using namespace std;
+
+
+int main( int argc, char *argv[] ) {
+ if ( argc != 2 ) {
+ cout << "Usage: " << argv[0] << " <filename.fg>" << endl << endl;
+ cout << "Reads factor graph <filename.fg> and runs" << endl;
+ cout << "Belief Propagation and JunctionTree on it." << endl << endl;
+ return 1;
+ } else {
+ FactorGraph fg;
+
+ if( fg.ReadFromFile(argv[1]) ) {
+ cout << "Error reading " << argv[1] << endl;
+ return 2;
+ } else {
+ size_t maxiter = 1000;
+ double tol = 1e-9;
+ size_t verb = 1;
+
+ Properties opts;
+ opts.Set("maxiter",maxiter);
+ opts.Set("tol",tol);
+ opts.Set("verbose",verb);
+
+ JTree jt( fg, opts("updates",string("HUGIN")) );
+ jt.init();
+ jt.run();
+
+ cout << "Exact single node marginals:" << endl;
+ for( size_t i = 0; i < fg.nrVars(); i++ )
+ cout << jt.belief(fg.var(i)) << endl;
+
+ BP bp(fg, opts("updates",string("SEQMAX")));
+ bp.init();
+ bp.run();
+
+ cout << "Exact single node marginals:" << endl;
+ for( size_t i = 0; i < fg.nrVars(); i++ )
+ cout << bp.belief(fg.var(i)) << endl;
+ }
+ }
+
+ return 0;
+}
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __FACTOR_H__
+#define __FACTOR_H__
+
+
+#include <iostream>
+#include <cmath>
+#include "prob.h"
+#include "varset.h"
+#include "index.h"
+
+
+template<typename T> class TFactor;
+typedef TFactor<Real> Factor;
+typedef TFactor<Complex> CFactor;
+
+
+// predefine friends
+template<typename T> Real dist( const TFactor<T> & x, const TFactor<T> & y, Prob::DistType dt );
+template<typename T> Complex KL_dist( const TFactor<T> & p, const TFactor<T> & q );
+template<typename T> std::ostream& operator<< (std::ostream& os, const TFactor<T>& P);
+
+
+// T should be castable from and to double and to complex
+template <typename T> class TFactor {
+ protected:
+ VarSet _vs;
+ TProb<T> _p;
+
+ public:
+ // Default constructor
+ TFactor () : _vs(), _p(1,1.0) {}
+
+ // Construct Factor from VarSet
+ TFactor( const VarSet& ns ) : _vs(ns), _p(_vs.stateSpace()) {}
+
+ // Construct Factor from VarSet and initial value
+ TFactor( const VarSet& ns, Real p ) : _vs(ns), _p(_vs.stateSpace(),p) {}
+
+ // Construct Factor from VarSet and initial array
+ TFactor( const VarSet& ns, const Real* p ) : _vs(ns), _p(_vs.stateSpace(),p) {}
+
+ // Construct Factor from VarSet and TProb<T>
+ TFactor( const VarSet& ns, const TProb<T> p ) : _vs(ns), _p(p) {
+#ifdef DEBUG
+ assert( _vs.stateSpace() == _p.size() );
+#endif
+ }
+
+ // Construct Factor from Var
+ TFactor( const Var& n ) : _vs(n), _p(n.states()) {}
+
+ // Copy constructor
+ TFactor( const TFactor<T> &x ) : _vs(x._vs), _p(x._p) {}
+
+ // Assignment operator
+ TFactor<T> & operator= (const TFactor<T> &x) {
+ if( this != &x ) {
+ _vs = x._vs;
+ _p = x._p;
+ }
+ return *this;
+ }
+
+ const TProb<T> & p() const { return _p; }
+ TProb<T> & p() { return _p; }
+ const VarSet & vars() const { return _vs; }
+ size_t stateSpace() const {
+#ifdef DEBUG
+ assert( _vs.stateSpace() == _p.size() );
+#endif
+ return _p.size();
+ }
+
+ T operator[] (size_t i) const { return _p[i]; }
+ T& operator[] (size_t i) { return _p[i]; }
+ TFactor<T> & fill (T p)
+ { _p.fill( p ); return(*this); }
+ TFactor<T> & randomize ()
+ { _p.randomize(); return(*this); }
+ TFactor<T> operator* (T x) const {
+ Factor result = *this;
+ result.p() *= x;
+ return result;
+ }
+ TFactor<T>& operator*= (T x) {
+ _p *= x;
+ return *this;
+ }
+ TFactor<T> operator/ (T x) const {
+ Factor result = *this;
+ result.p() /= x;
+ return result;
+ }
+ TFactor<T>& operator/= (T x) {
+ _p /= x;
+ return *this;
+ }
+ TFactor<T> operator* (const TFactor<T>& Q) const;
+ TFactor<T>& operator*= (const TFactor<T>& Q) { return( *this = (*this * Q) ); }
+ TFactor<T> operator+ (const TFactor<T>& Q) const {
+#ifdef DEBUG
+ assert( Q._vs == _vs );
+#endif
+ TFactor<T> sum(*this);
+ sum._p += Q._p;
+ return sum;
+ }
+ TFactor<T> operator- (const TFactor<T>& Q) const {
+#ifdef DEBUG
+ assert( Q._vs == _vs );
+#endif
+ TFactor<T> sum(*this);
+ sum._p -= Q._p;
+ return sum;
+ }
+ TFactor<T>& operator+= (const TFactor<T>& Q) {
+#ifdef DEBUG
+ assert( Q._vs == _vs );
+#endif
+ _p += Q._p;
+ return *this;
+ }
+ TFactor<T>& operator-= (const TFactor<T>& Q) {
+#ifdef DEBUG
+ assert( Q._vs == _vs );
+#endif
+ _p -= Q._p;
+ return *this;
+ }
+
+ TFactor<T> operator^ (Real a) const { TFactor<T> x; x._vs = _vs; x._p = _p^a; return x; }
+ TFactor<T>& operator^= (Real a) { _p ^= a; return *this; }
+
+ TFactor<T>& makeZero( Real epsilon ) {
+ _p.makeZero( epsilon );
+ return *this;
+ }
+
+ TFactor<T> inverse() const {
+ TFactor<T> inv;
+ inv._vs = _vs;
+ inv._p = _p.inverse(true); // FIXME
+ return inv;
+ }
+
+ TFactor<T> divided_by( const TFactor<T>& denom ) const {
+#ifdef DEBUG
+ assert( denom._vs == _vs );
+#endif
+ TFactor<T> quot(*this);
+ quot._p /= denom._p;
+ return quot;
+ }
+
+ TFactor<T>& divide( const TFactor<T>& denom ) {
+#ifdef DEBUG
+ assert( denom._vs == _vs );
+#endif
+ _p /= denom._p;
+ return *this;
+ }
+
+ TFactor<T> exp() const {
+ TFactor<T> e;
+ e._vs = _vs;
+ e._p = _p.exp();
+ return e;
+ }
+
+ TFactor<T> log() const {
+ TFactor<T> l;
+ l._vs = _vs;
+ l._p = _p.log();
+ return l;
+ }
+
+ TFactor<T> log0() const {
+ TFactor<T> l0;
+ l0._vs = _vs;
+ l0._p = _p.log0();
+ return l0;
+ }
+
+ CFactor clog0() const {
+ CFactor l0;
+ l0._vs = _vs;
+ l0._p = _p.clog0();
+ return l0;
+ }
+
+ T normalize( typename Prob::NormType norm ) { return _p.normalize( norm ); }
+ TFactor<T> normalized( typename Prob::NormType norm ) const {
+ TFactor<T> result;
+ result._vs = _vs;
+ result._p = _p.normalized( norm );
+ return result;
+ }
+
+ // returns slice of this factor where the subset ns is in state ns_state
+ Factor slice( const VarSet & ns, size_t ns_state ) const {
+ assert( ns << _vs );
+ VarSet nsrem = _vs / ns;
+ Factor result( nsrem, 0.0 );
+
+ // OPTIMIZE ME
+ Index i_ns (ns, _vs);
+ Index i_nsrem (nsrem, _vs);
+ for( size_t i = 0; i < stateSpace(); i++, ++i_ns, ++i_nsrem )
+ if( (size_t)i_ns == ns_state )
+ result._p[i_nsrem] = _p[i];
+
+ return result;
+ }
+
+ // returns unnormalized marginal
+ TFactor<T> part_sum(const VarSet & ns) const;
+ // returns normalized marginal
+ TFactor<T> marginal(const VarSet & ns) const { return part_sum(ns).normalized( Prob::NORMPROB ); }
+
+ bool hasNaNs() const { return _p.hasNaNs(); }
+ bool hasNegatives() const { return _p.hasNegatives(); }
+ T totalSum() const { return _p.totalSum(); }
+ T maxAbs() const { return _p.maxAbs(); }
+ T max() const { return _p.max(); }
+ Complex entropy() const { return _p.entropy(); }
+ T strength( const Var &i, const Var &j ) const;
+
+ friend Real dist( const TFactor<T> & x, const TFactor<T> & y, Prob::DistType dt ) {
+ if( x._vs.empty() || y._vs.empty() )
+ return -1;
+ else {
+#ifdef DEBUG
+ assert( x._vs == y._vs );
+#endif
+ return dist( x._p, y._p, dt );
+ }
+ }
+ friend Complex KL_dist <> (const TFactor<T> & p, const TFactor<T> & q);
+ template<class U> friend std::ostream& operator<< (std::ostream& os, const TFactor<U>& P);
+};
+
+
+template<typename T> TFactor<T> TFactor<T>::part_sum(const VarSet & ns) const {
+#ifdef DEBUG
+ assert( ns << _vs );
+#endif
+
+ TFactor<T> res( ns, 0.0 );
+
+ Index i_res( ns, _vs );
+ for( size_t i = 0; i < _p.size(); i++, ++i_res )
+ res._p[i_res] += _p[i];
+
+ return res;
+}
+
+
+template<typename T> std::ostream& operator<< (std::ostream& os, const TFactor<T>& P) {
+ os << "(" << P.vars() << " <";
+ for( size_t i = 0; i < P._p.size(); i++ )
+ os << P._p[i] << " ";
+ os << ">)";
+ return os;
+}
+
+
+template<typename T> TFactor<T> TFactor<T>::operator* (const TFactor<T>& Q) const {
+ TFactor<T> prod( _vs | Q._vs, 0.0 );
+
+ Index i1(_vs, prod._vs);
+ Index i2(Q._vs, prod._vs);
+
+ for( size_t i = 0; i < prod._p.size(); i++, ++i1, ++i2 )
+ prod._p[i] += _p[i1] * Q._p[i2];
+
+ return prod;
+}
+
+
+template<typename T> Complex KL_dist(const TFactor<T> & P, const TFactor<T> & Q) {
+ if( P._vs.empty() || Q._vs.empty() )
+ return -1;
+ else {
+#ifdef DEBUG
+ assert( P._vs == Q._vs );
+#endif
+ return KL_dist( P._p, Q._p );
+ }
+}
+
+
+// calculate N(psi, i, j)
+template<typename T> T TFactor<T>::strength( const Var &i, const Var &j ) const {
+#ifdef DEBUG
+ assert( _vs && i );
+ assert( _vs && j );
+ assert( i != j );
+#endif
+ VarSet ij = i | j;
+
+ T max = 0.0;
+ for( size_t alpha1 = 0; alpha1 < i.states(); alpha1++ )
+ for( size_t alpha2 = 0; alpha2 < i.states(); alpha2++ )
+ if( alpha2 != alpha1 )
+ for( size_t beta1 = 0; beta1 < j.states(); beta1++ )
+ for( size_t beta2 = 0; beta2 < j.states(); beta2++ )
+ if( beta2 != beta1 ) {
+ size_t as = 1, bs = 1;
+ if( i < j )
+ bs = i.states();
+ else
+ as = j.states();
+ T f1 = slice( ij, alpha1 * as + beta1 * bs ).p().divide( slice( ij, alpha2 * as + beta1 * bs ).p() ).max();
+ T f2 = slice( ij, alpha2 * as + beta2 * bs ).p().divide( slice( ij, alpha1 * as + beta2 * bs ).p() ).max();
+ T f = f1 * f2;
+ if( f > max )
+ max = f;
+ }
+
+ return std::tanh( 0.25 * std::log( max ) );
+}
+
+
+template<typename T> TFactor<T> RemoveFirstOrderInteractions( const TFactor<T> & psi ) {
+ TFactor<T> result = psi;
+
+ VarSet vars = psi.vars();
+ for( size_t iter = 0; iter < 100; iter++ ) {
+ for( VarSet::const_iterator n = vars.begin(); n != vars.end(); n++ )
+ result = result * result.part_sum(*n).inverse();
+ result.normalize( Prob::NORMPROB );
+ }
+
+ return result;
+}
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#include <iostream>
+#include <iterator>
+#include <map>
+#include <set>
+#include <fstream>
+#include <string>
+#include <algorithm>
+#include <functional>
+#include "factorgraph.h"
+
+
+using namespace std;
+
+
+FactorGraph::FactorGraph( const vector<Factor> &P ) : BipartiteGraph<Var,Factor>(), _undoProbs(), _hasNegatives(false), _normtype(Prob::NORMPROB) {
+ // add Factors
+ set<Var> _vars;
+ for(vector<Factor>::const_iterator p2 = P.begin(); p2 != P.end(); p2++ ) {
+ V2s().push_back(*p2);
+ if( !_hasNegatives && p2->hasNegatives() )
+ _hasNegatives = true;
+ for( set<Var>::const_iterator i = p2->vars().begin(); i != p2->vars().end(); i++ )
+ _vars.insert(*i);
+ }
+
+ // if negative factors are present, use LINF norm
+ if( _hasNegatives )
+ _normtype = Prob::NORMLINF;
+
+ // add _vars
+ for(VarSet::const_iterator p1 = _vars.begin(); p1 != _vars.end(); p1++ )
+ vars().push_back(*p1);
+
+ // create edges
+ for(size_t i2 = 0; i2 < nrFactors(); i2++ ) {
+ VarSet ns = factor(i2).vars();
+ for(VarSet::const_iterator q = ns.begin(); q != ns.end(); q++ ) {
+ for(size_t i1=0; i1 < nrVars(); i1++ ) {
+ if (var(i1) == *q) {
+ edges().push_back(_edge_t(i1,i2));
+ break;
+ }
+ }
+ }
+ }
+
+ // calc neighbours and adjacency matrix
+ Regenerate();
+
+ // Check for short loops
+ if( hasShortLoops(P) )
+ cerr << "FactorGraph::FactorGraph(): WARNING: short loops are present" << endl;
+}
+
+
+/*FactorGraph& FactorGraph::addFactor( const Factor &I ) {
+ // add Factor
+ _V2.push_back( I );
+ if( !_hasNegatives && I.hasNegatives() )
+ _hasNegatives = true;
+
+ // if negative factors are present, use LINF norm
+ if( _hasNegatives )
+ _normtype = Prob::NORMLINF;
+
+ // add new vars in Factor
+ for( VarSet::const_iterator i = I.vars().begin(); i != I.vars().end(); i++ ) {
+ size_t i_ind = find(vars().begin(), vars().end(), *i) - vars().begin();
+ if( i_ind == vars().size() )
+ _V1.push_back( *i );
+ _E12.push_back( _edge_t( i_ind, nrFactors() - 1 ) );
+ }
+
+ Regenerate();
+ return(*this);
+}*/
+
+
+ostream& operator << (ostream& os, const FactorGraph& fg) {
+ os << fg.nrFactors() << endl;
+
+ for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+ os << endl;
+ os << fg.factor(I).vars().size() << endl;
+ for( VarSet::const_iterator i = fg.factor(I).vars().begin(); i != fg.factor(I).vars().end(); i++ )
+ os << i->label() << " ";
+ os << endl;
+ for( VarSet::const_iterator i = fg.factor(I).vars().begin(); i != fg.factor(I).vars().end(); i++ )
+ os << i->states() << " ";
+ os << endl;
+ size_t nr_nonzeros = 0;
+ for( size_t k = 0; k < fg.factor(I).stateSpace(); k++ )
+ if( fg.factor(I)[k] != 0.0 )
+ nr_nonzeros++;
+ os << nr_nonzeros << endl;
+ for( size_t k = 0; k < fg.factor(I).stateSpace(); k++ )
+ if( fg.factor(I)[k] != 0.0 ) {
+ char buf[20];
+ sprintf(buf,"%18.14g", fg.factor(I)[k]);
+ os << k << " " << buf << endl;
+ }
+ }
+
+ return(os);
+}
+
+
+istream& operator >> (istream& is, FactorGraph& fg) {
+ long verbose = 0;
+
+ try {
+ vector<Factor> factors;
+ size_t nr_f;
+ string line;
+
+ while( (is.peek()) == '#' )
+ getline(is,line);
+ is >> nr_f;
+ if( is.fail() )
+ throw "ReadFromFile: unable to read number of Factors";
+ if( verbose >= 2 )
+ cout << "Reading " << nr_f << " factors..." << endl;
+
+ getline (is,line);
+ if( is.fail() )
+ throw "ReadFromFile: empty line expected";
+
+ for( size_t I = 0; I < nr_f; I++ ) {
+ if( verbose >= 3 )
+ cout << "Reading factor " << I << "..." << endl;
+ size_t nr_members;
+ while( (is.peek()) == '#' )
+ getline(is,line);
+ is >> nr_members;
+ if( verbose >= 3 )
+ cout << " nr_members: " << nr_members << endl;
+
+ vector<long> labels;
+ for( size_t mi = 0; mi < nr_members; mi++ ) {
+ long mi_label;
+ while( (is.peek()) == '#' )
+ getline(is,line);
+ is >> mi_label;
+ labels.push_back(mi_label);
+ }
+ if( verbose >= 3 ) {
+ cout << " labels: ";
+ copy (labels.begin(), labels.end(), ostream_iterator<int>(cout, " "));
+ cout << endl;
+ }
+
+ vector<size_t> dims;
+ for( size_t mi = 0; mi < nr_members; mi++ ) {
+ size_t mi_dim;
+ while( (is.peek()) == '#' )
+ getline(is,line);
+ is >> mi_dim;
+ dims.push_back(mi_dim);
+ }
+ if( verbose >= 3 ) {
+ cout << " dimensions: ";
+ copy (dims.begin(), dims.end(), ostream_iterator<int>(cout, " "));
+ cout << endl;
+ }
+
+ // add the Factor
+ VarSet I_vars;
+ for( size_t mi = 0; mi < nr_members; mi++ )
+ I_vars.insert( Var(labels[mi], dims[mi]) );
+ factors.push_back(Factor(I_vars,0.0));
+
+ // calculate permutation sigma (internally, members are sorted)
+ vector<long> sigma(nr_members,0);
+ VarSet::iterator j = I_vars.begin();
+ for( size_t mi = 0; mi < nr_members; mi++,j++ ) {
+ long search_for = j->label();
+ vector<long>::iterator j_loc = find(labels.begin(),labels.end(),search_for);
+ sigma[mi] = j_loc - labels.begin();
+ }
+ if( verbose >= 3 ) {
+ cout << " sigma: ";
+ copy( sigma.begin(), sigma.end(), ostream_iterator<int>(cout," "));
+ cout << endl;
+ }
+
+ // calculate multindices
+ vector<size_t> sdims(nr_members,0);
+ for( size_t k = 0; k < nr_members; k++ ) {
+ sdims[k] = dims[sigma[k]];
+ }
+ multind mi(dims);
+ multind smi(sdims);
+ if( verbose >= 3 ) {
+ cout << " mi.max(): " << mi.max() << endl;
+ cout << " ";
+ for( size_t k=0; k < nr_members; k++ )
+ cout << labels[k] << " ";
+ cout << " ";
+ for( size_t k=0; k < nr_members; k++ )
+ cout << labels[sigma[k]] << " ";
+ cout << endl;
+ }
+
+ // read values
+ size_t nr_nonzeros;
+ while( (is.peek()) == '#' )
+ getline(is,line);
+ is >> nr_nonzeros;
+ if( verbose >= 3 )
+ cout << " nonzeroes: " << nr_nonzeros << endl;
+ for( size_t k = 0; k < nr_nonzeros; k++ ) {
+ size_t li;
+ double val;
+ while( (is.peek()) == '#' )
+ getline(is,line);
+ is >> li;
+ while( (is.peek()) == '#' )
+ getline(is,line);
+ is >> val;
+
+ vector<size_t> vi = mi.vi(li);
+ vector<size_t> svi(vi.size(),0);
+ for( size_t k = 0; k < vi.size(); k++ )
+ svi[k] = vi[sigma[k]];
+ size_t sli = smi.li(svi);
+ if( verbose >= 3 ) {
+ cout << " " << li << ": ";
+ copy( vi.begin(), vi.end(), ostream_iterator<size_t>(cout," "));
+ cout << "-> ";
+ copy( svi.begin(), svi.end(), ostream_iterator<size_t>(cout," "));
+ cout << ": " << sli << endl;
+ }
+ factors.back()[sli] = val;
+ }
+ }
+
+ if( verbose >= 3 ) {
+ cout << "factors:" << endl;
+ copy(factors.begin(), factors.end(), ostream_iterator<Factor>(cout,"\n"));
+ }
+
+ fg = FactorGraph(factors);
+ } catch (char *e) {
+ cout << e << endl;
+ }
+
+ return is;
+}
+
+
+VarSet FactorGraph::delta(const Var & n) const {
+ // calculate Markov Blanket
+ size_t i = findVar( n );
+
+ VarSet del;
+ for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); I++ )
+ for( _nb_cit j = nb2(*I).begin(); j != nb2(*I).end(); j++ )
+ if( *j != i )
+ del |= var(*j);
+
+ return del;
+}
+
+
+VarSet FactorGraph::Delta(const Var & n) const {
+ return( delta(n) | n );
+}
+
+
+void FactorGraph::makeFactorCavity(size_t I) {
+ // fill Factor I with ones
+ factor(I).fill(1.0);
+}
+
+
+void FactorGraph::makeCavity(const Var & n) {
+ // fills all Factors that include Var n with ones
+ size_t i = findVar( n );
+
+ for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); I++ )
+ factor(*I).fill(1.0);
+}
+
+
+/*FactorGraph & FactorGraph::DeleteFactor(size_t I) {
+ // Go through all edges
+ for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
+ if( edge->second >= I ) {
+ if( edge->second == I )
+ edge->second = -1UL;
+ else
+ (edge->second)--;
+ }
+ // Remove all edges containing I
+ for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
+ if( edge->second == -1UL )
+ edge = _E12.erase( edge );
+// vector<_edge_t>::iterator new_end = _E12.remove_if( _E12.begin(), _E12.end(), compose1( bind2nd(equal_to<size_t>(), -1), select2nd<_edge_t>() ) );
+// _E12.erase( new_end, _E12.end() );
+
+ // Erase the factor
+ _V2.erase( _V2.begin() + I );
+
+ Regenerate();
+
+ return *this;
+}
+
+
+FactorGraph & FactorGraph::DeleteVar(size_t i) {
+ // Go through all edges
+ for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
+ if( edge->first >= i ) {
+ if( edge->first == i )
+ edge->first = -1UL;
+ else
+ (edge->first)--;
+ }
+ // Remove all edges containing i
+ for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
+ if( edge->first == -1UL )
+ edge = _E12.erase( edge );
+
+// vector<_edge_t>::iterator new_end = _E12.remove_if( _E12.begin(), _E12.end(), compose1( bind2nd(equal_to<size_t>(), -1), select1st<_edge_t>() ) );
+// _E12.erase( new_end, _E12.end() );
+
+ // Erase the variable
+ _V1.erase( _V1.begin() + i );
+
+ Regenerate();
+
+ return *this;
+}*/
+
+
+long FactorGraph::ReadFromFile(const char *filename) {
+ ifstream infile;
+ infile.open (filename);
+ if (infile.is_open()) {
+ infile >> *this;
+ infile.close();
+ return 0;
+ } else {
+ cout << "ERROR OPENING FILE" << endl;
+ return 1;
+ }
+}
+
+
+long FactorGraph::WriteToFile(const char *filename) const {
+ ofstream outfile;
+ outfile.open (filename);
+ if (outfile.is_open()) {
+ try {
+ outfile << *this;
+ } catch (char *e) {
+ cout << e << endl;
+ return 1;
+ }
+ outfile.close();
+ return 0;
+ } else {
+ cout << "ERROR OPENING FILE" << endl;
+ return 1;
+ }
+}
+
+
+long FactorGraph::WriteToDotFile(const char *filename) const {
+ ofstream outfile;
+ outfile.open (filename);
+ if (outfile.is_open()) {
+ try {
+ outfile << "graph G {" << endl;
+ outfile << "graph[size=\"9,9\"];" << endl;
+ outfile << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
+ for( size_t i = 0; i < nrVars(); i++ )
+ outfile << "\tx" << var(i).label() << ";" << endl;
+ outfile << "node[shape=box,style=filled,color=lightgrey,width=0.3,height=0.3,fixedsize=true];" << endl;
+ for( size_t I = 0; I < nrFactors(); I++ )
+ outfile << "\tp" << I << ";" << endl;
+ for( size_t iI = 0; iI < nr_edges(); iI++ )
+ outfile << "\tx" << var(edge(iI).first).label() << " -- p" << edge(iI).second << ";" << endl;
+ outfile << "}" << endl;
+ } catch (char *e) {
+ cout << e << endl;
+ return 1;
+ }
+ outfile.close();
+ return 0;
+ } else {
+ cout << "ERROR OPENING FILE" << endl;
+ return 1;
+ }
+}
+
+
+bool hasShortLoops( const vector<Factor> &P ) {
+ bool found = false;
+ vector<Factor>::const_iterator I, J;
+ for( I = P.begin(); I != P.end(); I++ ) {
+ J = I;
+ J++;
+ for( ; J != P.end(); J++ )
+ if( (I->vars() & J->vars()).size() >= 2 ) {
+ found = true;
+ break;
+ }
+ if( found )
+ break;
+ }
+ return found;
+}
+
+
+void RemoveShortLoops(vector<Factor> &P) {
+ bool found = true;
+ while( found ) {
+ found = false;
+ vector<Factor>::iterator I, J;
+ for( I = P.begin(); I != P.end(); I++ ) {
+ J = I;
+ J++;
+ for( ; J != P.end(); J++ )
+ if( (I->vars() & J->vars()).size() >= 2 ) {
+ found = true;
+ break;
+ }
+ if( found )
+ break;
+ }
+ if( found ) {
+ cout << "Merging factors " << I->vars() << " and " << J->vars() << endl;
+ *I *= *J;
+ P.erase(J);
+ }
+ }
+}
+
+
+Factor FactorGraph::ExactMarginal(const VarSet & x) const {
+ Factor P;
+ for( size_t I = 0; I < nrFactors(); I++ )
+ P *= factor(I);
+ return P.marginal(x);
+}
+
+
+Real FactorGraph::ExactlogZ() const {
+ Factor P;
+ for( size_t I = 0; I < nrFactors(); I++ )
+ P *= factor(I);
+ return std::log(P.totalSum());
+}
+
+
+vector<VarSet> FactorGraph::Cliques() const {
+ vector<VarSet> result;
+
+ for( size_t I = 0; I < nrFactors(); I++ ) {
+ bool maximal = true;
+ for( size_t J = 0; (J < nrFactors()) && maximal; J++ )
+ if( (factor(J).vars() >> factor(I).vars()) && !(factor(J).vars() == factor(I).vars()) )
+ maximal = false;
+
+ if( maximal )
+ result.push_back( factor(I).vars() );
+ }
+
+ return result;
+}
+
+
+void FactorGraph::clamp( const Var & n, size_t i ) {
+ assert( i <= n.states() );
+
+/* if( do_surgery ) {
+ size_t ni = find( vars().begin(), vars().end(), n) - vars().begin();
+
+ if( ni != nrVars() ) {
+ for( _nb_cit I = nb1(ni).begin(); I != nb1(ni).end(); I++ ) {
+ if( factor(*I).size() == 1 )
+ // Remove this single-variable factor
+ // I = (_V2.erase(I))--;
+ _E12.erase( _E12.begin() + VV2E(ni, *I) );
+ else {
+ // Replace it by the slice
+ Index ind_I_min_n( factor(*I), factor(*I) / n );
+ Index ind_n( factor(*I), n );
+ Factor slice_I( factor(*I) / n );
+ for( size_t ind_I = 0; ind_I < factor(*I).stateSpace(); ++ind_I, ++ind_I_min_n, ++ind_n )
+ if( ind_n == i )
+ slice_I[ind_I_min_n] = factor(*I)[ind_I];
+ factor(*I) = slice_I;
+
+ // Remove the edge between n and I
+ _E12.erase( _E12.begin() + VV2E(ni, *I) );
+ }
+ }
+
+ Regenerate();
+
+ // remove all unconnected factors
+ for( size_t I = 0; I < nrFactors(); I++ )
+ if( nb2(I).size() == 0 )
+ DeleteFactor(I--);
+
+ DeleteVar( ni );
+
+ // FIXME
+ }
+ } */
+
+ // The cheap solution (at least in terms of coding time) is to multiply every factor
+ // that contains the variable with a delta function
+
+ Factor delta_n_i(n,0.0);
+ delta_n_i[i] = 1.0;
+
+ // For all factors that contain n
+ for( size_t I = 0; I < nrFactors(); I++ )
+ if( factor(I).vars() && n )
+ // Multiply it with a delta function
+ factor(I) *= delta_n_i;
+
+ return;
+}
+
+
+void FactorGraph::saveProb( size_t I ) {
+ map<size_t,Prob>::iterator it = _undoProbs.find( I );
+ if( it != _undoProbs.end() )
+ cout << "FactorGraph::saveProb: WARNING: _undoProbs[I] already defined!" << endl;
+ _undoProbs[I] = factor(I).p();
+}
+
+
+void FactorGraph::undoProb( size_t I ) {
+ map<size_t,Prob>::iterator it = _undoProbs.find( I );
+ if( it != _undoProbs.end() ) {
+ factor(I).p() = (*it).second;
+ _undoProbs.erase(it);
+ }
+}
+
+
+void FactorGraph::saveProbs( const VarSet &ns ) {
+ if( !_undoProbs.empty() )
+ cout << "FactorGraph::saveProbs: WARNING: _undoProbs not empy!" << endl;
+ for( size_t I = 0; I < nrFactors(); I++ )
+ if( factor(I).vars() && ns )
+ _undoProbs[I] = factor(I).p();
+}
+
+
+void FactorGraph::undoProbs( const VarSet &ns ) {
+ for( map<size_t,Prob>::iterator uI = _undoProbs.begin(); uI != _undoProbs.end(); ) {
+ if( factor((*uI).first).vars() && ns ) {
+// cout << "undoing " << factor((*uI).first).vars() << endl;
+// cout << "from " << factor((*uI).first).p() << " to " << (*uI).second << endl;
+ factor((*uI).first).p() = (*uI).second;
+ _undoProbs.erase(uI++);
+ } else
+ uI++;
+ }
+}
+
+
+bool FactorGraph::isConnected() const {
+ if( nrVars() == 0 )
+ return false;
+ else {
+ Var n = var( 0 );
+
+ VarSet component = n;
+
+ VarSet remaining;
+ for( size_t i = 1; i < nrVars(); i++ )
+ remaining |= var(i);
+
+ bool found_new_vars = true;
+ while( found_new_vars ) {
+ VarSet new_vars;
+ for( VarSet::const_iterator m = remaining.begin(); m != remaining.end(); m++ )
+ if( delta(*m) && component )
+ new_vars |= *m;
+
+ if( new_vars.empty() )
+ found_new_vars = false;
+ else
+ found_new_vars = true;
+
+ component |= new_vars;
+ remaining /= new_vars;
+ };
+ return remaining.empty();
+ }
+}
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __FACTORGRAPH_H__
+#define __FACTORGRAPH_H__
+
+
+#include <iostream>
+#include <map>
+#include "bipgraph.h"
+#include "factor.h"
+
+
+using namespace std;
+
+
+class FactorGraph : public BipartiteGraph<Var,Factor> {
+ protected:
+ map<size_t,Prob> _undoProbs;
+ bool _hasNegatives;
+ Prob::NormType _normtype;
+
+ public:
+ /// Default constructor
+ FactorGraph() : BipartiteGraph<Var,Factor>(), _undoProbs(), _hasNegatives(false), _normtype(Prob::NORMPROB) {};
+ /// Copy constructor
+ FactorGraph(const FactorGraph & x) : BipartiteGraph<Var,Factor>(x), _undoProbs(), _hasNegatives(x._hasNegatives), _normtype(x._normtype) {};
+ /// Construct FactorGraph from vector of Factors
+ FactorGraph(const vector<Factor> &P);
+ /// Assignment operator
+ FactorGraph & operator=(const FactorGraph & x) {
+ if(this!=&x) {
+ BipartiteGraph<Var,Factor>::operator=(x);
+ _undoProbs = x._undoProbs;
+ _hasNegatives = x._hasNegatives;
+ _normtype = x._normtype;
+ }
+ return *this;
+ }
+ virtual ~FactorGraph() {}
+
+ // aliases
+ Var & var(size_t i) { return V1(i); }
+ const Var & var(size_t i) const { return V1(i); }
+ const vector<Var> & vars() const { return V1s(); }
+ vector<Var> & vars() { return V1s(); }
+ size_t nrVars() const { return V1s().size(); }
+ Factor & factor(size_t I) { return V2(I); }
+ const Factor & factor(size_t I) const { return V2(I); }
+ const vector<Factor> & factors() const { return V2s(); }
+ vector<Factor> & factors() { return V2s(); }
+ size_t nrFactors() const { return V2s().size(); }
+
+ /// Provides read access to neighbours of variable
+ const _nb_t & nbV( size_t i1 ) const { return nb1(i1); }
+ /// Provides full access to neighbours of variable
+ _nb_t & nbV( size_t i1 ) { return nb1(i1); }
+ /// Provides read access to neighbours of factor
+ const _nb_t & nbF( size_t i2 ) const { return nb2(i2); }
+ /// Provides full access to neighbours of factor
+ _nb_t & nbF( size_t i2 ) { return nb2(i2); }
+
+ size_t findVar(const Var & n) const {
+ size_t i = find( vars().begin(), vars().end(), n ) - vars().begin();
+ assert( i != nrVars() );
+ return i;
+ }
+ size_t findFactor(const VarSet &ns) const {
+ size_t I;
+ for( I = 0; I < nrFactors(); I++ )
+ if( factor(I).vars() == ns )
+ break;
+ assert( I != nrFactors() );
+ return I;
+ }
+
+ friend ostream& operator << (ostream& os, const FactorGraph& fg);
+ friend istream& operator >> (istream& is, FactorGraph& fg);
+
+ VarSet delta(const Var & n) const;
+ VarSet Delta(const Var & n) const;
+ virtual void makeFactorCavity(size_t I);
+ virtual void makeCavity(const Var & n);
+
+ long ReadFromFile(const char *filename);
+ long WriteToFile(const char *filename) const;
+ long WriteToDotFile(const char *filename) const;
+
+ Factor ExactMarginal(const VarSet & x) const;
+ Real ExactlogZ() const;
+
+ virtual void clamp( const Var & n, size_t i );
+
+ bool hasNegatives() const { return _hasNegatives; }
+ Prob::NormType NormType() const { return _normtype; }
+
+ vector<VarSet> Cliques() const;
+
+ virtual void undoProbs( const VarSet &ns );
+ void saveProbs( const VarSet &ns );
+ virtual void undoProb( size_t I );
+ void saveProb( size_t I );
+
+ bool isConnected() const;
+
+ virtual void updatedFactor( size_t I ) {};
+};
+
+
+bool hasShortLoops(const vector<Factor> &P);
+void RemoveShortLoops(vector<Factor> &P);
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#include <map>
+#include "hak.h"
+#include "util.h"
+#include "diffs.h"
+
+
+const char *HAK::Name = "HAK";
+
+
+bool HAK::checkProperties() {
+ if( !HasProperty("tol") )
+ return false;
+ if (!HasProperty("maxiter") )
+ return false;
+ if (!HasProperty("verbose") )
+ return false;
+ if( !HasProperty("doubleloop") )
+ return false;
+ if( !HasProperty("clusters") )
+ return false;
+
+ ConvertPropertyTo<double>("tol");
+ ConvertPropertyTo<size_t>("maxiter");
+ ConvertPropertyTo<size_t>("verbose");
+ ConvertPropertyTo<bool>("doubleloop");
+ ConvertPropertyTo<ClustersType>("clusters");
+
+ if( HasProperty("loopdepth") )
+ ConvertPropertyTo<size_t>("loopdepth");
+ else if( Clusters() == ClustersType::LOOP )
+ return false;
+
+ return true;
+}
+
+
+void HAK::constructMessages() {
+ // Create outer beliefs
+ _Qa.clear();
+ _Qa.reserve(nr_ORs());
+ for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+ _Qa.push_back( Factor( OR(alpha).vars() ) );
+
+ // Create inner beliefs
+ _Qb.clear();
+ _Qb.reserve(nr_IRs());
+ for( size_t beta = 0; beta < nr_IRs(); beta++ )
+ _Qb.push_back( Factor( IR(beta) ) );
+
+ // Create messages
+ _muab.clear();
+ _muab.reserve(nr_Redges());
+ _muba.clear();
+ _muba.reserve(nr_Redges());
+ for( vector<R_edge_t>::const_iterator ab = Redges().begin(); ab != Redges().end(); ab++ ) {
+ _muab.push_back( Factor( IR(ab->second) ) );
+ _muba.push_back( Factor( IR(ab->second) ) );
+ }
+}
+
+
+HAK::HAK(const RegionGraph & rg, const Properties &opts) : DAIAlgRG(rg, opts) {
+ assert( checkProperties() );
+
+ constructMessages();
+}
+
+
+void HAK::findLoopClusters( const FactorGraph & fg, set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
+ for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
+ VarSet ind = fg.delta( *in );
+ if( (newcl.size()) >= 2 && (ind >> root) ) {
+ allcl.insert( newcl | *in );
+ }
+ else if( length > 1 )
+ findLoopClusters( fg, allcl, newcl | *in, root, length - 1, ind / newcl );
+ }
+}
+
+
+HAK::HAK(const FactorGraph & fg, const Properties &opts) : DAIAlgRG(opts) {
+ assert( checkProperties() );
+
+ vector<VarSet> cl;
+ if( Clusters() == ClustersType::MIN ) {
+ cl = fg.Cliques();
+ } else if( Clusters() == ClustersType::DELTA ) {
+ for( size_t i = 0; i < fg.nrVars(); i++ )
+ cl.push_back(fg.Delta(fg.var(i)));
+ } else if( Clusters() == ClustersType::LOOP ) {
+ cl = fg.Cliques();
+ set<VarSet> scl;
+ for( vector<Var>::const_iterator i0 = fg.vars().begin(); i0 != fg.vars().end(); i0++ ) {
+ VarSet i0d = fg.delta(*i0);
+ if( LoopDepth() > 1 )
+ findLoopClusters( fg, scl, *i0, *i0, LoopDepth() - 1, fg.delta(*i0) );
+ }
+ for( set<VarSet>::const_iterator c = scl.begin(); c != scl.end(); c++ )
+ cl.push_back(*c);
+ if( Verbose() >= 3 ) {
+ cout << "HAK uses the following clusters: " << endl;
+ for( vector<VarSet>::const_iterator cli = cl.begin(); cli != cl.end(); cli++ )
+ cout << *cli << endl;
+ }
+ } else
+ throw "Invalid Clusters type";
+
+ RegionGraph rg(fg,cl);
+ RegionGraph::operator=(rg);
+ constructMessages();
+
+ if( Verbose() >= 3 )
+ cout << "HAK regiongraph: " << *this << endl;
+}
+
+
+string HAK::identify() const {
+ stringstream result (stringstream::out);
+ result << Name << GetProperties();
+ return result.str();
+}
+
+
+void HAK::init( const VarSet &ns ) {
+ for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
+ if( alpha->vars() && ns )
+ alpha->fill( 1.0 / alpha->stateSpace() );
+
+ for( size_t beta = 0; beta < nr_IRs(); beta++ )
+ if( IR(beta) && ns ) {
+ _Qb[beta].fill( 1.0 / IR(beta).stateSpace() );
+ for( R_nb_cit alpha = nbIR(beta).begin(); alpha != nbIR(beta).end(); alpha++ ) {
+ muab(*alpha,beta).fill( 1.0 / IR(beta).stateSpace() );
+ muba(beta,*alpha).fill( 1.0 / IR(beta).stateSpace() );
+ }
+ }
+}
+
+
+void HAK::init() {
+ assert( checkProperties() );
+
+ for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
+ alpha->fill( 1.0 / alpha->stateSpace() );
+
+ for( vector<Factor>::iterator beta = _Qb.begin(); beta != _Qb.end(); beta++ )
+ beta->fill( 1.0 / beta->stateSpace() );
+
+ for( size_t ab = 0; ab < nr_Redges(); ab++ ) {
+ _muab[ab].fill( 1.0 / _muab[ab].stateSpace() );
+ _muba[ab].fill( 1.0 / _muba[ab].stateSpace() );
+ }
+}
+
+
+double HAK::doGBP() {
+ if( Verbose() >= 1 )
+ cout << "Starting " << identify() << "...";
+ if( Verbose() >= 3)
+ cout << endl;
+
+ clock_t tic = toc();
+
+ // Check whether counting numbers won't lead to problems
+ for( size_t beta = 0; beta < nr_IRs(); beta++ )
+ assert( nbIR(beta).size() + IR(beta).c() != 0.0 );
+
+ // Keep old beliefs to check convergence
+ vector<Factor> old_beliefs;
+ old_beliefs.reserve( nrVars() );
+ for( size_t i = 0; i < nrVars(); i++ )
+ old_beliefs.push_back( belief( var(i) ) );
+
+ // Differences in single node beliefs
+ Diffs diffs(nrVars(), 1.0);
+
+ size_t iter = 0;
+ // do several passes over the network until maximum number of iterations has
+ // been reached or until the maximum belief difference is smaller than tolerance
+ for( iter = 0; iter < MaxIter() && diffs.max() > Tol(); iter++ ) {
+ for( size_t beta = 0; beta < nr_IRs(); beta++ ) {
+ for( R_nb_cit alpha = nbIR(beta).begin(); alpha != nbIR(beta).end(); alpha++ )
+ muab(*alpha,beta) = _Qa[*alpha].marginal(IR(beta)).divided_by( muba(beta,*alpha) );
+
+ Factor Qb_new;
+ for( R_nb_cit alpha = nbIR(beta).begin(); alpha != nbIR(beta).end(); alpha++ )
+ Qb_new *= muab(*alpha,beta) ^ (1 / (nbIR(beta).size() + IR(beta).c()));
+ Qb_new.normalize( _normtype );
+ if( Qb_new.hasNaNs() ) {
+ cout << "HAK::doGBP: Qb_new has NaNs!" << endl;
+ return NAN;
+ }
+// _Qb[beta] = Qb_new.makeZero(1e-100); // damping?
+ _Qb[beta] = Qb_new;
+
+ for( R_nb_cit alpha = nbIR(beta).begin(); alpha != nbIR(beta).end(); alpha++ ) {
+ muba(beta,*alpha) = _Qb[beta].divided_by( muab(*alpha,beta) );
+
+ Factor Qa_new = OR(*alpha);
+ for( R_nb_cit gamma = nbOR(*alpha).begin(); gamma != nbOR(*alpha).end(); gamma++ )
+ Qa_new *= muba(*gamma,*alpha);
+ Qa_new ^= (1.0 / OR(*alpha).c());
+ Qa_new.normalize( _normtype );
+ if( Qa_new.hasNaNs() ) {
+ cout << "HAK::doGBP: Qa_new has NaNs!" << endl;
+ return NAN;
+ }
+// _Qa[*alpha] = Qa_new.makeZero(1e-100); // damping?
+ _Qa[*alpha] = Qa_new;
+ }
+ }
+
+ // Calculate new single variable beliefs and compare with old ones
+ for( size_t i = 0; i < nrVars(); i++ ) {
+ Factor new_belief = belief( var( i ) );
+ diffs.push( dist( new_belief, old_beliefs[i], Prob::DISTLINF ) );
+ old_beliefs[i] = new_belief;
+ }
+
+ if( Verbose() >= 3 )
+ cout << "HAK::doGBP: maxdiff " << diffs.max() << " after " << iter+1 << " passes" << endl;
+ }
+
+ updateMaxDiff( diffs.max() );
+
+ if( Verbose() >= 1 ) {
+ if( diffs.max() > Tol() ) {
+ if( Verbose() == 1 )
+ cout << endl;
+ cout << "HAK::doGBP: WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.max() << endl;
+ } else {
+ if( Verbose() >= 2 )
+ cout << "HAK::doGBP: ";
+ cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
+ }
+ }
+
+ return diffs.max();
+}
+
+
+double HAK::doDoubleLoop() {
+ if( Verbose() >= 1 )
+ cout << "Starting " << identify() << "...";
+ if( Verbose() >= 3)
+ cout << endl;
+
+ clock_t tic = toc();
+
+ // Save original outer regions
+ vector<FRegion> org_ORs = ORs();
+
+ // Save original inner counting numbers and set negative counting numbers to zero
+ vector<double> org_IR_cs( nr_IRs(), 0.0 );
+ for( size_t beta = 0; beta < nr_IRs(); beta++ ) {
+ org_IR_cs[beta] = IR(beta).c();
+ if( IR(beta).c() < 0.0 )
+ IR(beta).c() = 0.0;
+ }
+
+ // Keep old beliefs to check convergence
+ vector<Factor> old_beliefs;
+ old_beliefs.reserve( nrVars() );
+ for( size_t i = 0; i < nrVars(); i++ )
+ old_beliefs.push_back( belief( var(i) ) );
+
+ // Differences in single node beliefs
+ Diffs diffs(nrVars(), 1.0);
+
+ size_t outer_maxiter = MaxIter();
+ double outer_tol = Tol();
+ size_t outer_verbose = Verbose();
+ double org_maxdiff = MaxDiff();
+
+ // Set parameters for inner loop
+ MaxIter( 5 );
+ Verbose( outer_verbose ? outer_verbose - 1 : 0 );
+
+ size_t outer_iter = 0;
+ for( outer_iter = 0; outer_iter < outer_maxiter && diffs.max() > outer_tol; outer_iter++ ) {
+ // Calculate new outer regions
+ for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
+ OR(alpha) = org_ORs[alpha];
+ for( R_nb_cit beta = nbOR(alpha).begin(); beta != nbOR(alpha).end(); beta++ )
+ OR(alpha) *= _Qb[*beta] ^ ((IR(*beta).c() - org_IR_cs[*beta]) / nbIR(*beta).size());
+ }
+
+ // Inner loop
+ if( isnan( doGBP() ) )
+ return NAN;
+
+ // Calculate new single variable beliefs and compare with old ones
+ for( size_t i = 0; i < nrVars(); i++ ) {
+ Factor new_belief = belief( var( i ) );
+ diffs.push( dist( new_belief, old_beliefs[i], Prob::DISTLINF ) );
+ old_beliefs[i] = new_belief;
+ }
+
+ if( Verbose() >= 3 )
+ cout << "HAK::doDoubleLoop: maxdiff " << diffs.max() << " after " << outer_iter+1 << " passes" << endl;
+ }
+
+ // restore _maxiter, _verbose and _maxdiff
+ MaxIter( outer_maxiter );
+ Verbose( outer_verbose );
+ MaxDiff( org_maxdiff );
+
+ updateMaxDiff( diffs.max() );
+
+ // Restore original outer regions
+ ORs() = org_ORs;
+
+ // Restore original inner counting numbers
+ for( size_t beta = 0; beta < nr_IRs(); beta++ )
+ IR(beta).c() = org_IR_cs[beta];
+
+ if( Verbose() >= 1 ) {
+ if( diffs.max() > Tol() ) {
+ if( Verbose() == 1 )
+ cout << endl;
+ cout << "HAK::doDoubleLoop: WARNING: not converged within " << outer_maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.max() << endl;
+ } else {
+ if( Verbose() >= 3 )
+ cout << "HAK::doDoubleLoop: ";
+ cout << "converged in " << outer_iter << " passes (" << toc() - tic << " clocks)." << endl;
+ }
+ }
+
+ return diffs.max();
+}
+
+
+double HAK::run() {
+ if( DoubleLoop() )
+ return doDoubleLoop();
+ else
+ return doGBP();
+}
+
+
+Factor HAK::belief( const VarSet &ns ) const {
+ vector<Factor>::const_iterator beta;
+ for( beta = _Qb.begin(); beta != _Qb.end(); beta++ )
+ if( beta->vars() >> ns )
+ break;
+ if( beta != _Qb.end() )
+ return( beta->marginal(ns) );
+ else {
+ vector<Factor>::const_iterator alpha;
+ for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
+ if( alpha->vars() >> ns )
+ break;
+ assert( alpha != _Qa.end() );
+ return( alpha->marginal(ns) );
+ }
+}
+
+
+Factor HAK::belief( const Var &n ) const {
+ return belief( (VarSet)n );
+}
+
+
+vector<Factor> HAK::beliefs() const {
+ vector<Factor> result;
+ for( size_t beta = 0; beta < nr_IRs(); beta++ )
+ result.push_back( Qb(beta) );
+ for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+ result.push_back( Qa(alpha) );
+ return result;
+}
+
+
+Complex HAK::logZ() const {
+ Complex sum = 0.0;
+ for( size_t beta = 0; beta < nr_IRs(); beta++ )
+ sum += Complex(IR(beta).c()) * Qb(beta).entropy();
+ for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
+ sum += Complex(OR(alpha).c()) * Qa(alpha).entropy();
+ sum += (OR(alpha).log0() * Qa(alpha)).totalSum();
+ }
+ return sum;
+}
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __HAK_H__
+#define __HAK_H__
+
+
+#include "daialg.h"
+#include "regiongraph.h"
+#include "enum.h"
+
+
+using namespace std;
+
+
+/// HAK provides an implementation of the single and double-loop algorithms by Heskes, Albers and Kappen
+class HAK : public DAIAlgRG {
+ protected:
+ vector<Factor> _Qa;
+ vector<Factor> _Qb;
+ vector<Factor> _muab;
+ vector<Factor> _muba;
+
+ public:
+ /// Default constructor
+ HAK() : DAIAlgRG() {};
+
+ /// Copy constructor
+ HAK(const HAK & x) : DAIAlgRG(x), _Qa(x._Qa), _Qb(x._Qb), _muab(x._muab), _muba(x._muba) {};
+
+ /// Clone function
+ HAK* clone() const { return new HAK(*this); }
+
+ /// Construct from RegionGraph
+ HAK(const RegionGraph & rg, const Properties &opts);
+
+ /// Construct from RactorGraph using "clusters" option
+ HAK(const FactorGraph & fg, const Properties &opts);
+
+ /// Assignment operator
+ HAK & operator=(const HAK & x) {
+ if( this != &x ) {
+ DAIAlgRG::operator=(x);
+ _Qa = x._Qa;
+ _Qb = x._Qb;
+ _muab = x._muab;
+ _muba = x._muba;
+ }
+ return *this;
+ }
+
+ static const char *Name;
+
+ ENUM3(ClustersType,MIN,DELTA,LOOP)
+ ClustersType Clusters() const { return GetPropertyAs<ClustersType>("clusters"); }
+ bool DoubleLoop() { return GetPropertyAs<bool>("doubleloop"); }
+ size_t LoopDepth() { return GetPropertyAs<size_t>("loopdepth"); }
+
+ Factor & muab( size_t alpha, size_t beta ) { return _muab[ORIR2E(alpha,beta)]; }
+ Factor & muba( size_t beta, size_t alpha ) { return _muba[ORIR2E(alpha,beta)]; }
+ const Factor& Qa( size_t alpha ) const { return _Qa[alpha]; };
+ const Factor& Qb( size_t beta ) const { return _Qb[beta]; };
+
+// void Regenerate();
+ double doGBP();
+ double doDoubleLoop();
+ double run();
+ void init();
+ string identify() const;
+ Factor belief( const Var &n ) const;
+ Factor belief( const VarSet &ns ) const;
+ vector<Factor> beliefs() const;
+ Complex logZ () const;
+
+ void init( const VarSet &ns );
+ void undoProbs( const VarSet &ns ) { RegionGraph::undoProbs( ns ); init( ns ); }
+ bool checkProperties();
+
+ private:
+ void constructMessages();
+ void findLoopClusters( const FactorGraph &fg, set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars );
+};
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __INDEX_H__
+#define __INDEX_H__
+
+#include <vector>
+#include "varset.h"
+
+/* Example:
+ *
+ * Index i ({s_j_1,s_j_2,...,s_j_m}, {s_1,...,s_N}); // j_k in {1,...,N}
+ * for( ; i>=0; ++i ) {
+ * // loops over all states of (s_1,...,s_N)
+ * // i is linear index of corresponding state of (s_j_1, ..., s_j_m)
+ * }
+ */
+
+class Index
+{
+private:
+ long _index;
+ std::vector<int> _count,_max,_sum;
+public:
+ Index () { _index=-1; };
+ Index (const VarSet& P, const VarSet& ns)
+ {
+ long sum=1;
+ VarSet::const_iterator j=ns.begin();
+ for(VarSet::const_iterator i=P.begin();i!=P.end();++i)
+ {
+ for(;j!=ns.end()&&j->label()<=i->label();++j)
+ {
+ _count.push_back(0);
+ _max.push_back(j->states());
+ _sum.push_back((i->label()==j->label())?sum:0);
+ };
+ sum*=i->states();
+ };
+ for(;j!=ns.end();++j)
+ {
+ _count.push_back(0);
+ _max.push_back(j->states());
+ _sum.push_back(0);
+ };
+ _index=0;
+ };
+ Index (const Index & ind) : _index(ind._index), _count(ind._count), _max(ind._max), _sum(ind._sum) {};
+ Index & operator=(const Index & ind) {
+ if(this!=&ind) {
+ _index = ind._index;
+ _count = ind._count;
+ _max = ind._max;
+ _sum = ind._sum;
+ }
+ return *this;
+ }
+ Index& clear ()
+ {
+ for(unsigned i=0;i!=_count.size();++i) _count[i]=0;
+ _index=0;
+ return(*this);
+ };
+ operator long () const { return(_index); };
+ Index& operator ++ ()
+ {
+ if(_index>=0)
+ {
+ unsigned i;
+ for(i=0;(i<_count.size())
+ &&(_index+=_sum[i],++_count[i]==_max[i]);++i)
+ {
+ _index-=_sum[i]*_max[i];
+ _count[i]=0;
+ };
+ if(i==_count.size()) _index=-1;
+ };
+ return(*this);
+ };
+};
+
+class multind {
+ private:
+ std::vector<size_t> _dims; // dimensions
+ std::vector<size_t> _pdims; // products of dimensions
+
+ public:
+ multind(const std::vector<size_t> di) {
+ _dims = di;
+ size_t prod = 1;
+ for( std::vector<size_t>::const_iterator i=di.begin(); i!=di.end(); i++ ) {
+ _pdims.push_back(prod);
+ prod = prod * (*i);
+ }
+ _pdims.push_back(prod);
+ }
+ multind(const VarSet& ns) {
+ _dims.reserve( ns.size() );
+ _pdims.reserve( ns.size() + 1 );
+ size_t prod = 1;
+ for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ ) {
+ _pdims.push_back( prod );
+ prod *= n->states();
+ _dims.push_back( n->states() );
+ }
+ _pdims.push_back( prod );
+ }
+ std::vector<size_t> vi(size_t li) const { // linear index to vector index
+ std::vector<size_t> v(_dims.size(),0);
+ assert(li >= 0 && li < _pdims.back());
+ for( long j = v.size()-1; j >= 0; j-- ) {
+ size_t q = li / _pdims[j];
+ v[j] = q;
+ li = li - q * _pdims[j];
+ }
+ return v;
+ }
+ size_t li(const std::vector<size_t> vi) const { // linear index
+ size_t s = 0;
+ assert(vi.size() == _dims.size());
+ for( size_t j = 0; j < vi.size(); j++ )
+ s += vi[j] * _pdims[j];
+ return s;
+ }
+ size_t max() const { return( _pdims.back() ); };
+
+ // FIXME add an iterator, which increases a vector index just using addition
+};
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#include <iostream>
+#include "jtree.h"
+
+
+using namespace std;
+
+
+const char *JTree::Name = "JTREE";
+
+
+bool JTree::checkProperties() {
+ if (!HasProperty("verbose") )
+ return false;
+ if( !HasProperty("updates") )
+ return false;
+
+ ConvertPropertyTo<size_t>("verbose");
+ ConvertPropertyTo<UpdateType>("updates");
+
+ return true;
+}
+
+
+JTree::JTree( const FactorGraph &fg, const Properties &opts, bool automatic) : DAIAlgRG(fg, opts), _RTree(), _Qa(), _Qb(), _mes(), _logZ() {
+ assert( checkProperties() );
+
+ if( automatic ) {
+ ClusterGraph _cg;
+
+ // Copy factors
+ for( size_t I = 0; I < nrFactors(); I++ )
+ _cg.insert( factor(I).vars() );
+ if( Verbose() >= 3 )
+ cout << "Initial clusters: " << _cg << endl;
+
+ // Retain only maximal clusters
+ _cg.eraseNonMaximal();
+ if( Verbose() >= 3 )
+ cout << "Maximal clusters: " << _cg << endl;
+
+ vector<VarSet> ElimVec = _cg.VarElim_MinFill().eraseNonMaximal().toVector();
+ if( Verbose() >= 3 ) {
+ cout << "VarElim_MinFill result: {" << endl;
+ for( size_t i = 0; i < ElimVec.size(); i++ ) {
+ if( i != 0 )
+ cout << ", ";
+ cout << ElimVec[i];
+ }
+ cout << "}" << endl;
+ }
+
+ GenerateJT( ElimVec );
+ }
+}
+
+
+void JTree::GenerateJT( const vector<VarSet> &Cliques ) {
+ // Construct a weighted graph (each edge is weighted with the cardinality
+ // of the intersection of the nodes, where the nodes are the elements of
+ // Cliques).
+ WeightedGraph<int> JuncGraph;
+ for( size_t i = 0; i < Cliques.size(); i++ )
+ for( size_t j = i+1; j < Cliques.size(); j++ ) {
+ size_t w = (Cliques[i] & Cliques[j]).size();
+ JuncGraph[UEdge(i,j)] = w;
+ }
+
+ // Construct maximal spanning tree using Prim's algorithm
+ _RTree = MaxSpanningTreePrim( JuncGraph );
+
+ // Construct corresponding region graph
+
+ // Create outer regions
+ ORs().reserve( Cliques.size() );
+ for( size_t i = 0; i < Cliques.size(); i++ )
+ ORs().push_back( FRegion( Factor(Cliques[i], 1.0), 1.0 ) );
+
+ // For each factor, find an outer region that subsumes that factor.
+ // Then, multiply the outer region with that factor.
+ for( size_t I = 0; I < nrFactors(); I++ ) {
+ size_t alpha;
+ for( alpha = 0; alpha < nr_ORs(); alpha++ )
+ if( OR(alpha).vars() >> factor(I).vars() ) {
+// OR(alpha) *= factor(I);
+ _fac2OR[I] = alpha;
+ break;
+ }
+ assert( alpha != nr_ORs() );
+ }
+ RecomputeORs();
+
+ // Create inner regions and edges
+ IRs().reserve( _RTree.size() );
+ Redges().reserve( 2 * _RTree.size() );
+ for( size_t i = 0; i < _RTree.size(); i++ ) {
+ Redges().push_back( R_edge_t( _RTree[i].n1, IRs().size() ) );
+ Redges().push_back( R_edge_t( _RTree[i].n2, IRs().size() ) );
+ // inner clusters have counting number -1
+ IRs().push_back( Region( Cliques[_RTree[i].n1] & Cliques[_RTree[i].n2], -1.0 ) );
+ }
+
+ // Regenerate BipartiteGraph internals
+ Regenerate();
+
+ // Create messages and beliefs
+ _Qa.clear();
+ _Qa.reserve( nr_ORs() );
+ for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+ _Qa.push_back( OR(alpha) );
+
+ _Qb.clear();
+ _Qb.reserve( nr_IRs() );
+ for( size_t beta = 0; beta < nr_IRs(); beta++ )
+ _Qb.push_back( Factor( IR(beta), 1.0 ) );
+
+ _mes.clear();
+ _mes.reserve( nr_Redges() );
+ for( size_t e = 0; e < nr_Redges(); e++ )
+ _mes.push_back( Factor( IR(Redge(e).second), 1.0 ) );
+
+ // Check counting numbers
+ Check_Counting_Numbers();
+
+ if( Verbose() >= 3 ) {
+ cout << "Resulting regiongraph: " << *this << endl;
+ }
+}
+
+
+string JTree::identify() const {
+ stringstream result (stringstream::out);
+ result << Name << GetProperties();
+ return result.str();
+}
+
+
+Factor JTree::belief( const VarSet &ns ) const {
+ vector<Factor>::const_iterator beta;
+ for( beta = _Qb.begin(); beta != _Qb.end(); beta++ )
+ if( beta->vars() >> ns )
+ break;
+ if( beta != _Qb.end() )
+ return( beta->marginal(ns) );
+ else {
+ vector<Factor>::const_iterator alpha;
+ for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
+ if( alpha->vars() >> ns )
+ break;
+ assert( alpha != _Qa.end() );
+ return( alpha->marginal(ns) );
+ }
+}
+
+
+vector<Factor> JTree::beliefs() const {
+ vector<Factor> result;
+ for( size_t beta = 0; beta < nr_IRs(); beta++ )
+ result.push_back( _Qb[beta] );
+ for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+ result.push_back( _Qa[alpha] );
+ return result;
+}
+
+
+Factor JTree::belief( const Var &n ) const {
+ return belief( (VarSet)n );
+}
+
+
+// Needs no init
+void JTree::runHUGIN() {
+ for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+ _Qa[alpha] = OR(alpha);
+
+ for( size_t beta = 0; beta < nr_IRs(); beta++ )
+ _Qb[beta].fill( 1.0 );
+
+ // CollectEvidence
+ _logZ = 0.0;
+ for( size_t i = _RTree.size(); (i--) != 0; ) {
+// Make outer region _RTree[i].n1 consistent with outer region _RTree[i].n2
+// IR(i) = seperator OR(_RTree[i].n1) && OR(_RTree[i].n2)
+ Factor new_Qb = _Qa[_RTree[i].n2].part_sum( IR( i ) );
+ _logZ += log(new_Qb.normalize( Prob::NORMPROB ));
+ _Qa[_RTree[i].n1] *= new_Qb.divided_by( _Qb[i] );
+ _Qb[i] = new_Qb;
+ }
+ if( _RTree.empty() )
+ _logZ += log(_Qa[0].normalize( Prob::NORMPROB ) );
+ else
+ _logZ += log(_Qa[_RTree[0].n1].normalize( Prob::NORMPROB ));
+
+ // DistributeEvidence
+ for( size_t i = 0; i < _RTree.size(); i++ ) {
+// Make outer region _RTree[i].n2 consistent with outer region _RTree[i].n1
+// IR(i) = seperator OR(_RTree[i].n1) && OR(_RTree[i].n2)
+ Factor new_Qb = _Qa[_RTree[i].n1].marginal( IR( i ) );
+ _Qa[_RTree[i].n2] *= new_Qb.divided_by( _Qb[i] );
+ _Qb[i] = new_Qb;
+ }
+
+ // Normalize
+ for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+ _Qa[alpha].normalize( Prob::NORMPROB );
+}
+
+
+// Really needs no init! Initial messages can be anything.
+void JTree::runShaferShenoy() {
+ // First pass
+ _logZ = 0.0;
+ for( size_t e = _RTree.size(); (e--) != 0; ) {
+ // send a message from _RTree[e].n2 to _RTree[e].n1
+ // or, actually, from the seperator IR(e) to _RTree[e].n1
+
+ size_t i = _RTree[e].n2;
+ size_t j = _RTree[e].n1;
+
+ Factor piet = OR(i);
+ for( R_nb_cit k = nbOR(i).begin(); k != nbOR(i).end(); k++ )
+ if( *k != e )
+ piet *= message( i, *k );
+ message( j, e ) = piet.part_sum( IR(e) );
+ _logZ += log( message(j,e).normalize( Prob::NORMPROB ) );
+ }
+
+ // Second pass
+ for( size_t e = 0; e < _RTree.size(); e++ ) {
+ size_t i = _RTree[e].n1;
+ size_t j = _RTree[e].n2;
+
+ Factor piet = OR(i);
+ for( R_nb_cit k = nbOR(i).begin(); k != nbOR(i).end(); k++ )
+ if( *k != e )
+ piet *= message( i, *k );
+ message( j, e ) = piet.marginal( IR(e) );
+ }
+
+ // Calculate beliefs
+ for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
+ Factor piet = OR(alpha);
+ for( R_nb_cit k = nbOR(alpha).begin(); k != nbOR(alpha).end(); k++ )
+ piet *= message( alpha, *k );
+ if( _RTree.empty() ) {
+ _logZ += log( piet.normalize( Prob::NORMPROB ) );
+ _Qa[alpha] = piet;
+ } else if( alpha == _RTree[0].n1 ) {
+ _logZ += log( piet.normalize( Prob::NORMPROB ) );
+ _Qa[alpha] = piet;
+ } else
+ _Qa[alpha] = piet.normalized( Prob::NORMPROB );
+ }
+
+ // Only for logZ (and for belief)...
+ for( size_t beta = 0; beta < nr_IRs(); beta++ )
+ _Qb[beta] = _Qa[nbIR(beta)[0]].marginal( IR(beta) );
+}
+
+
+double JTree::run() {
+ if( Updates() == UpdateType::HUGIN )
+ runHUGIN();
+ else if( Updates() == UpdateType::SHSH )
+ runShaferShenoy();
+ return 0.0;
+}
+
+
+Complex JTree::logZ() const {
+ Complex sum = 0.0;
+ for( size_t beta = 0; beta < nr_IRs(); beta++ )
+ sum += Complex(IR(beta).c()) * _Qb[beta].entropy();
+ for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
+ sum += Complex(OR(alpha).c()) * _Qa[alpha].entropy();
+ sum += (OR(alpha).log0() * _Qa[alpha]).totalSum();
+ }
+ return sum;
+}
+
+
+
+size_t JTree::findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t PreviousRoot ) const {
+ // find new root clique (the one with maximal statespace overlap with ns)
+ size_t maxval = 0, maxalpha = 0;
+ for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
+ size_t val = (ns & OR(alpha).vars()).stateSpace();
+ if( val > maxval ) {
+ maxval = val;
+ maxalpha = alpha;
+ }
+ }
+
+// for( size_t e = 0; e < _RTree.size(); e++ )
+// cout << OR(_RTree[e].n1).vars() << "->" << OR(_RTree[e].n2).vars() << ", ";
+// cout << endl;
+ // grow new tree
+ Graph oldTree;
+ for( DEdgeVec::const_iterator e = _RTree.begin(); e != _RTree.end(); e++ )
+ oldTree.insert( UEdge(e->n1, e->n2) );
+ DEdgeVec newTree = GrowRootedTree( oldTree, maxalpha );
+// cout << ns << ": ";
+// for( size_t e = 0; e < newTree.size(); e++ )
+// cout << OR(newTree[e].n1).vars() << "->" << OR(newTree[e].n2).vars() << ", ";
+// cout << endl;
+
+ // identify subtree that contains variables of ns which are not in the new root
+ VarSet nsrem = ns / OR(maxalpha).vars();
+// cout << "nsrem:" << nsrem << endl;
+ set<DEdge> subTree;
+ // for each variable in ns that is not in the root clique
+ for( VarSet::const_iterator n = nsrem.begin(); n != nsrem.end(); n++ ) {
+ // find first occurence of *n in the tree, which is closest to the root
+ size_t e = 0;
+ for( ; e != newTree.size(); e++ ) {
+ if( OR(newTree[e].n2).vars() && *n )
+ break;
+ }
+ assert( e != newTree.size() );
+
+ // track-back path to root and add edges to subTree
+ subTree.insert( newTree[e] );
+ size_t pos = newTree[e].n1;
+ for( ; e > 0; e-- )
+ if( newTree[e-1].n2 == pos ) {
+ subTree.insert( newTree[e-1] );
+ pos = newTree[e-1].n1;
+ }
+ }
+ if( PreviousRoot != (size_t)-1 && PreviousRoot != maxalpha) {
+ // find first occurence of PreviousRoot in the tree, which is closest to the new root
+ size_t e = 0;
+ for( ; e != newTree.size(); e++ ) {
+ if( newTree[e].n2 == PreviousRoot )
+ break;
+ }
+ assert( e != newTree.size() );
+
+ // track-back path to root and add edges to subTree
+ subTree.insert( newTree[e] );
+ size_t pos = newTree[e].n1;
+ for( ; e > 0; e-- )
+ if( newTree[e-1].n2 == pos ) {
+ subTree.insert( newTree[e-1] );
+ pos = newTree[e-1].n1;
+ }
+ }
+// cout << "subTree: " << endl;
+// for( set<DEdge>::const_iterator sTi = subTree.begin(); sTi != subTree.end(); sTi++ )
+// cout << OR(sTi->n1).vars() << "->" << OR(sTi->n2).vars() << ", ";
+// cout << endl;
+
+ // Resulting Tree is a reordered copy of newTree
+ // First add edges in subTree to Tree
+ Tree.clear();
+ for( DEdgeVec::const_iterator e = newTree.begin(); e != newTree.end(); e++ )
+ if( subTree.count( *e ) ) {
+ Tree.push_back( *e );
+// cout << OR(e->n1).vars() << "->" << OR(e->n2).vars() << ", ";
+ }
+// cout << endl;
+ // Then add edges pointing away from nsrem
+ // FIXME
+/* for( DEdgeVec::const_iterator e = newTree.begin(); e != newTree.end(); e++ )
+ for( set<DEdge>::const_iterator sTi = subTree.begin(); sTi != subTree.end(); sTi++ )
+ if( *e != *sTi ) {
+ if( e->n1 == sTi->n1 || e->n1 == sTi->n2 ||
+ e->n2 == sTi->n1 || e->n2 == sTi->n2 ) {
+ Tree.push_back( *e );
+// cout << OR(e->n1).vars() << "->" << OR(e->n2).vars() << ", ";
+ }
+ }*/
+ // FIXME
+/* for( DEdgeVec::const_iterator e = newTree.begin(); e != newTree.end(); e++ )
+ if( find( Tree.begin(), Tree.end(), *e) == Tree.end() ) {
+ bool found = false;
+ for( VarSet::const_iterator n = nsrem.begin(); n != nsrem.end(); n++ )
+ if( (OR(e->n1).vars() && *n) ) {
+ found = true;
+ break;
+ }
+ if( found ) {
+ Tree.push_back( *e );
+ cout << OR(e->n1).vars() << "->" << OR(e->n2).vars() << ", ";
+ }
+ }
+ cout << endl;*/
+ size_t subTreeSize = Tree.size();
+ // Then add remaining edges
+ for( DEdgeVec::const_iterator e = newTree.begin(); e != newTree.end(); e++ )
+ if( find( Tree.begin(), Tree.end(), *e ) == Tree.end() )
+ Tree.push_back( *e );
+
+ return subTreeSize;
+}
+
+
+// Cutset conditioning
+// assumes that run() has been called already
+Factor JTree::calcMarginal( const VarSet& ns ) {
+ vector<Factor>::const_iterator beta;
+ for( beta = _Qb.begin(); beta != _Qb.end(); beta++ )
+ if( beta->vars() >> ns )
+ break;
+ if( beta != _Qb.end() )
+ return( beta->marginal(ns) );
+ else {
+ vector<Factor>::const_iterator alpha;
+ for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
+ if( alpha->vars() >> ns )
+ break;
+ if( alpha != _Qa.end() )
+ return( alpha->marginal(ns) );
+ else {
+ // Find subtree to do efficient inference
+ DEdgeVec T;
+ size_t Tsize = findEfficientTree( ns, T );
+
+ // Find remaining variables (which are not in the new root)
+ VarSet nsrem = ns / OR(T.front().n1).vars();
+ Factor Pns (ns, 0.0);
+
+ multind mi( nsrem );
+
+ // Save _Qa and _Qb on the subtree
+ map<size_t,Factor> _Qa_old;
+ map<size_t,Factor> _Qb_old;
+ vector<size_t> b(Tsize, 0);
+ for( size_t i = Tsize; (i--) != 0; ) {
+ size_t alpha1 = T[i].n1;
+ size_t alpha2 = T[i].n2;
+ size_t beta;
+ for( beta = 0; beta < nr_IRs(); beta++ )
+ if( UEdge( _RTree[beta].n1, _RTree[beta].n2 ) == UEdge( alpha1, alpha2 ) )
+ break;
+ assert( beta != nr_IRs() );
+ b[i] = beta;
+
+ if( !_Qa_old.count( alpha1 ) )
+ _Qa_old[alpha1] = _Qa[alpha1];
+ if( !_Qa_old.count( alpha2 ) )
+ _Qa_old[alpha2] = _Qa[alpha2];
+ if( !_Qb_old.count( beta ) )
+ _Qb_old[beta] = _Qb[beta];
+ }
+
+ // For all states of nsrem
+ for( size_t j = 0; j < mi.max(); j++ ) {
+ vector<size_t> vi = mi.vi( j );
+
+ // CollectEvidence
+ double logZ = 0.0;
+ for( size_t i = Tsize; (i--) != 0; ) {
+ // Make outer region T[i].n1 consistent with outer region T[i].n2
+ // IR(i) = seperator OR(T[i].n1) && OR(T[i].n2)
+
+ size_t k = 0;
+ for( VarSet::const_iterator n = nsrem.begin(); n != nsrem.end(); n++, k++ )
+ if( _Qa[T[i].n2].vars() >> *n ) {
+ Factor piet( *n, 0.0 );
+ piet[vi[k]] = 1.0;
+ _Qa[T[i].n2] *= piet;
+ }
+
+ Factor new_Qb = _Qa[T[i].n2].part_sum( IR( b[i] ) );
+ logZ += log(new_Qb.normalize( Prob::NORMPROB ));
+ _Qa[T[i].n1] *= new_Qb.divided_by( _Qb[b[i]] );
+ _Qb[b[i]] = new_Qb;
+ }
+ logZ += log(_Qa[T[0].n1].normalize( Prob::NORMPROB ));
+
+ Factor piet( nsrem, 0.0 );
+ piet[j] = exp(logZ);
+ Pns += piet * _Qa[T[0].n1].part_sum( ns / nsrem ); // OPTIMIZE ME
+
+ // Restore clamped beliefs
+ for( map<size_t,Factor>::const_iterator alpha = _Qa_old.begin(); alpha != _Qa_old.end(); alpha++ )
+ _Qa[alpha->first] = alpha->second;
+ for( map<size_t,Factor>::const_iterator beta = _Qb_old.begin(); beta != _Qb_old.end(); beta++ )
+ _Qb[beta->first] = beta->second;
+ }
+
+ return( Pns.normalized(Prob::NORMPROB) );
+ }
+ }
+}
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __JTREE_H__
+#define __JTREE_H__
+
+
+#include <vector>
+#include "daialg.h"
+#include "varset.h"
+#include "regiongraph.h"
+#include "factorgraph.h"
+#include "clustergraph.h"
+#include "weightedgraph.h"
+#include "enum.h"
+
+
+using namespace std;
+
+
+class JTree : public DAIAlgRG {
+ protected:
+ DEdgeVec _RTree; // rooted tree
+ vector<Factor> _Qa;
+ vector<Factor> _Qb;
+ vector<Factor> _mes;
+ double _logZ;
+
+
+ public:
+ ENUM2(UpdateType,HUGIN,SHSH)
+ UpdateType Updates() const { return GetPropertyAs<UpdateType>("updates"); }
+
+ JTree() : DAIAlgRG(), _RTree(), _Qa(), _Qb(), _mes(), _logZ() {};
+ JTree( const JTree& x ) : DAIAlgRG(x), _RTree(x._RTree), _Qa(x._Qa), _Qb(x._Qb), _mes(x._mes), _logZ(x._logZ) {};
+ JTree* clone() const { return new JTree(*this); }
+ JTree & operator=( const JTree& x ) {
+ if( this != &x ) {
+ DAIAlgRG::operator=(x);
+ _RTree = x._RTree;
+ _Qa = x._Qa;
+ _Qb = x._Qb;
+ _mes = x._mes;
+ _logZ = x._logZ;
+ }
+ return *this;
+ }
+ JTree( const FactorGraph &fg, const Properties &opts, bool automatic=true );
+ void GenerateJT( const vector<VarSet> &Cliques );
+
+ Factor & message(size_t i1, size_t i2) { return( _mes[ORIR2E(i1,i2)] ); }
+ const Factor & message(size_t i1, size_t i2) const { return( _mes[ORIR2E(i1,i2)] ); }
+
+ static const char *Name;
+ string identify() const;
+// void Regenerate();
+ void init() {
+ assert( checkProperties() );
+ }
+ void runHUGIN();
+ void runShaferShenoy();
+ double run();
+ Factor belief( const Var &n ) const;
+ Factor belief( const VarSet &ns ) const;
+ vector<Factor> beliefs() const;
+ Complex logZ() const;
+
+ void init( const VarSet &ns ) {}
+ void undoProbs( const VarSet &ns ) { RegionGraph::undoProbs( ns ); init( ns ); }
+
+ size_t findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t PreviousRoot=(size_t)-1 ) const;
+ Factor calcMarginal( const VarSet& ns );
+ bool checkProperties();
+};
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#include <iostream>
+#include <algorithm>
+#include <map>
+#include <set>
+#include "lc.h"
+#include "diffs.h"
+#include "util.h"
+#include "alldai.h"
+#include "x2x.h"
+
+
+using namespace std;
+
+
+const char *LC::Name = "LC";
+
+
+bool LC::checkProperties() {
+ if( !HasProperty("cavity") )
+ return false;
+ if( !HasProperty("updates") )
+ return false;
+ if( !HasProperty("tol") )
+ return false;
+ if (!HasProperty("maxiter") )
+ return false;
+ if (!HasProperty("verbose") )
+ return false;
+
+ ConvertPropertyTo<CavityType>("cavity");
+ ConvertPropertyTo<UpdateType>("updates");
+ ConvertPropertyTo<double>("tol");
+ ConvertPropertyTo<size_t>("maxiter");
+ ConvertPropertyTo<size_t>("verbose");
+
+ if (HasProperty("cavainame") )
+ ConvertPropertyTo<string>("cavainame");
+ if (HasProperty("cavaiopts") )
+ ConvertPropertyTo<Properties>("cavaiopts");
+ if( HasProperty("reinit") )
+ ConvertPropertyTo<bool>("reinit");
+
+ return true;
+}
+
+
+LC::LC(const FactorGraph & fg, const Properties &opts) : DAIAlgFG(fg, opts) {
+ assert( checkProperties() );
+
+ // calc iI
+ for( size_t i=0; i < nrVars(); i++ ) {
+ for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); I++ ) {
+ _iI_type _iI_entry;
+ _iI_entry.i = i;
+ _iI_entry.I = *I;
+
+ _iI.push_back(_iI_entry);
+ }
+ }
+
+ // create pancakes
+ _pancakes.resize(nrVars());
+
+ // create cavitydists
+ for( size_t i=0; i < nrVars(); i++ )
+ _cavitydists.push_back(Factor(delta(var(i))));
+
+ // create phis
+ _phis.reserve(nr_edges());
+ for( size_t iI = 0; iI < nr_edges(); iI++ ) {
+ size_t i = edge(iI).first;
+ size_t I = edge(iI).second;
+ _phis.push_back( Factor( factor(I).vars() / var(i) ) );
+ }
+
+ // create beliefs
+ for( size_t i=0; i < nrVars(); i++ )
+ _beliefs.push_back(Factor(var(i)));
+}
+
+
+string LC::identify() const {
+ stringstream result (stringstream::out);
+ result << Name << GetProperties();
+ return result.str();
+}
+
+
+void LC::CalcBelief (size_t i) {
+ _beliefs[i] = _pancakes[i].marginal(var(i));
+}
+
+
+double LC::CalcCavityDist (size_t i, const string &name, const Properties &opts) {
+ Factor Bi;
+ double maxdiff = 0;
+
+ if( Verbose() >= 2 )
+ cout << "Initing cavity " << var(i) << "(" << delta(var(i)).size() << " vars, " << delta(var(i)).stateSpace() << " states)" << endl;
+
+ if( Cavity() == CavityType::UNIFORM )
+ Bi = Factor(delta(var(i)));
+ else {
+ InfAlg *cav = newInfAlg( name, *this, opts );
+ cav->makeCavity( var(i) );
+
+ if( Cavity() == CavityType::FULL )
+ Bi = calcMarginal( *cav, cav->delta(var(i)), reInit() );
+ else if( Cavity() == CavityType::PAIR )
+ Bi = calcMarginal2ndO( *cav, cav->delta(var(i)), reInit() );
+ else if( Cavity() == CavityType::PAIR2 ) {
+ vector<Factor> pairbeliefs = calcPairBeliefsNew( *cav, cav->delta(var(i)), reInit() );
+ for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
+ Bi *= pairbeliefs[ij];
+ } else if( Cavity() == CavityType::PAIRINT ) {
+ Bi = calcMarginal( *cav, cav->delta(var(i)), reInit() );
+
+ // Set interactions of order > 2 to zero
+ size_t N = delta(var(i)).size();
+ double *p = &(*Bi.p().p().begin());
+ x2x::p2logp (N, p);
+ x2x::logp2w (N, p);
+ x2x::fill (N, p, 2, 0.0);
+ x2x::w2logp (N, p);
+// x2x::logpnorm (N, p);
+ x2x::logp2p (N, p);
+ } else if( Cavity() == CavityType::PAIRCUM ) {
+ Bi = calcMarginal( *cav, cav->delta(var(i)), reInit() );
+
+ // Set cumulants of order > 2 to zero
+ size_t N = delta(var(i)).size();
+ double *p = &(*Bi.p().p().begin());
+ x2x::p2m (N, p);
+ x2x::m2c (N, p, N);
+ x2x::fill (N, p, 2, 0.0);
+ x2x::c2m (N, p, N);
+ x2x::m2p (N, p);
+ }
+ maxdiff = cav->MaxDiff();
+ delete cav;
+ }
+ Bi.normalize( _normtype );
+ _cavitydists[i] = Bi;
+
+ return maxdiff;
+}
+
+
+double LC::InitCavityDists (const string &name, const Properties &opts) {
+ clock_t tic = toc();
+
+ if( Verbose() >= 1 ) {
+ cout << "LC::InitCavityDists: ";
+ if( Cavity() == CavityType::UNIFORM )
+ cout << "Using uniform initial cavity distributions" << endl;
+ else if( Cavity() == CavityType::FULL )
+ cout << "Using full " << name << opts << "...";
+ else if( Cavity() == CavityType::PAIR )
+ cout << "Using pairwise " << name << opts << "...";
+ else if( Cavity() == CavityType::PAIR2 )
+ cout << "Using pairwise(new) " << name << opts << "...";
+ }
+
+ double maxdiff = 0.0;
+ for( size_t i = 0; i < nrVars(); i++ ) {
+ double md = CalcCavityDist(i, name, opts);
+ if( md > maxdiff )
+ maxdiff = md;
+ }
+ init();
+
+ if( Verbose() >= 1 ) {
+ cout << "used " << toc() - tic << " clocks." << endl;
+ }
+
+ return maxdiff;
+}
+
+
+long LC::SetCavityDists (vector<Factor> &Q) {
+ if( Verbose() >= 1 )
+ cout << "LC::SetCavityDists: Setting initial cavity distributions" << endl;
+ if( Q.size() != nrVars() )
+ return -1;
+ for( size_t i = 0; i < nrVars(); i++ ) {
+ if( _cavitydists[i].vars() != Q[i].vars() ) {
+ return i+1;
+ } else
+ _cavitydists[i] = Q[i];
+ }
+ init();
+ return 0;
+}
+
+
+void LC::init() {
+ for( size_t iI = 0; iI < nr_edges(); iI++ ) {
+ if( Updates() == UpdateType::SEQRND )
+ _phis[iI].randomize();
+ else
+ _phis[iI].fill(1.0);
+ }
+ for( size_t i = 0; i < nrVars(); i++ ) {
+ _pancakes[i] = _cavitydists[i];
+
+ for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); I++ ) {
+ _pancakes[i] *= factor(*I);
+ if( Updates() == UpdateType::SEQRND )
+ _pancakes[i] *= _phis[VV2E(i,*I)];
+ }
+
+ _pancakes[i].normalize( _normtype );
+
+ CalcBelief(i);
+ }
+}
+
+
+Factor LC::NewPancake (size_t iI, bool & hasNaNs) {
+ size_t i = _iI[iI].i;
+ size_t I = _iI[iI].I;
+ iI = VV2E(i, I);
+
+ Factor piet = _pancakes[i];
+
+ // recalculate _pancake[i]
+ VarSet Ivars = factor(I).vars();
+ Factor A_I;
+ for( VarSet::const_iterator k = Ivars.begin(); k != Ivars.end(); k++ )
+ if( var(i) != *k )
+ A_I *= (_pancakes[findVar(*k)] * factor(I).inverse()).part_sum( Ivars / var(i) );
+ if( Ivars.size() > 1 )
+ A_I ^= (1.0 / (Ivars.size() - 1));
+ Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[iI].inverse()).part_sum( Ivars / var(i) );
+ Factor quot = A_I.divided_by(A_Ii);
+
+ piet *= quot.divided_by( _phis[iI] ).normalized( _normtype );
+ _phis[iI] = quot.normalized( _normtype );
+
+ piet.normalize( _normtype );
+
+ if( piet.hasNaNs() ) {
+ cout << "LC::NewPancake(" << iI << "): has NaNs!" << endl;
+ hasNaNs = true;
+ }
+
+ return piet;
+}
+
+
+double LC::run() {
+ if( Verbose() >= 1 )
+ cout << "Starting " << identify() << "...";
+ if( Verbose() >= 2 )
+ cout << endl;
+
+ clock_t tic = toc();
+ Diffs diffs(nrVars(), 1.0);
+
+ updateMaxDiff( InitCavityDists(GetPropertyAs<string>("cavainame"), GetPropertyAs<Properties>("cavaiopts")) );
+
+ vector<Factor> old_beliefs;
+ for(size_t i=0; i < nrVars(); i++ )
+ old_beliefs.push_back(belief(i));
+
+ bool hasNaNs = false;
+ for( size_t i=0; i < nrVars(); i++ )
+ if( _pancakes[i].hasNaNs() ) {
+ hasNaNs = true;
+ break;
+ }
+ if( hasNaNs ) {
+ cout << "LC::run: initial _pancakes has NaNs!" << endl;
+ return NAN;
+ }
+
+ vector<long> update_seq(nr_iI(),0);
+ for( size_t k=0; k < nr_iI(); k++ )
+ update_seq[k] = k;
+
+ size_t iter=0;
+
+ // do several passes over the network until maximum number of iterations has
+ // been reached or until the maximum belief difference is smaller than tolerance
+ for( iter=0; iter < MaxIter() && diffs.max() > Tol(); iter++ ) {
+ // Sequential updates
+ if( Updates() == UpdateType::SEQRND )
+ random_shuffle( update_seq.begin(), update_seq.end() );
+
+ for( size_t t=0; t < nr_iI(); t++ ) {
+ long iI = update_seq[t];
+ long i = _iI[iI].i;
+ _pancakes[i] = NewPancake(iI, hasNaNs);
+ if( hasNaNs )
+ return NAN;
+ CalcBelief(i);
+ }
+
+ // compare new beliefs with old ones
+ for(size_t i=0; i < nrVars(); i++ ) {
+ diffs.push( dist( belief(i), old_beliefs[i], Prob::DISTLINF ) );
+ old_beliefs[i] = belief(i);
+ }
+
+ if( Verbose() >= 3 )
+ cout << "LC::run: maxdiff " << diffs.max() << " after " << iter+1 << " passes" << endl;
+ }
+
+ updateMaxDiff( diffs.max() );
+
+ if( Verbose() >= 1 ) {
+ if( diffs.max() > Tol() ) {
+ if( Verbose() == 1 )
+ cout << endl;
+ cout << "LC::run: WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.max() << endl;
+ } else {
+ if( Verbose() >= 2 )
+ cout << "LC::run: ";
+ cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
+ }
+ }
+
+ return diffs.max();
+}
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __LC_H__
+#define __LC_H__
+
+
+#include "daialg.h"
+#include "enum.h"
+#include "factorgraph.h"
+
+
+using namespace std;
+
+
+class LC : public DAIAlgFG {
+ protected:
+ typedef struct { size_t i; size_t I; } _iI_type;
+
+ vector<Factor> _pancakes; // used by all LC types (psi_I is stored in the pancake)
+ vector<Factor> _cavitydists; // used by all LC types to store the approximate cavity distribution
+ /// _phis[VV2E(i,I)] corresponds to \f$ \phi^{\setminus i}_I(x_{I \setminus i}) \f$
+ vector<Factor> _phis;
+
+ /// Single variable beliefs
+ vector<Factor> _beliefs;
+
+ /// For each pair (i,j) with j in delta(i), store i and the common factor I
+ vector<_iI_type> _iI;
+
+ public:
+ ENUM6(CavityType,FULL,PAIR,PAIR2,PAIRINT,PAIRCUM,UNIFORM)
+ ENUM3(UpdateType,SEQFIX,SEQRND,NONE)
+
+ CavityType Cavity() const { return GetPropertyAs<CavityType>("cavity"); }
+ UpdateType Updates() const { return GetPropertyAs<UpdateType>("updates"); }
+ bool reInit() const { return GetPropertyAs<bool>("reinit"); }
+
+ /// Default constructor
+ LC() : DAIAlgFG() {};
+ /// Copy constructor
+ LC(const LC & x) : DAIAlgFG(x), _pancakes(x._pancakes), _cavitydists(x._cavitydists), _phis(x._phis), _beliefs(x._beliefs), _iI(x._iI) {};
+ /// Clone function
+ LC* clone() const { return new LC(*this); }
+ /// Construct LC object from a FactorGraph and parameters
+ LC(const FactorGraph & fg, const Properties &opts);
+ /// Assignment operator
+ LC& operator=(const LC & x) {
+ if( this != &x ) {
+ DAIAlgFG::operator=(x);
+ _pancakes = x._pancakes;
+ _cavitydists = x._cavitydists;
+ _phis = x._phis;
+ _beliefs = x._beliefs;
+ _iI = x._iI;
+ }
+ return *this;
+ }
+
+ static const char *Name;
+ double CalcCavityDist (size_t i, const string &name, const Properties &opts);
+ double InitCavityDists (const string &name, const Properties &opts);
+ long SetCavityDists (vector<Factor> &Q);
+
+ void init();
+ Factor NewPancake (size_t iI, bool & hasNaNs);
+ double run();
+
+ string identify() const;
+ Factor belief (const Var &n) const { return( _beliefs[findVar(n)] ); }
+ Factor belief (const VarSet &ns) const { assert( 0 == 1 ); }
+ vector<Factor> beliefs() const { return _beliefs; }
+ Complex logZ() const { return NAN; }
+ void CalcBelief (size_t i);
+ const Factor & belief (size_t i) const { return _beliefs[i]; };
+ const Factor & pancake (size_t i) const { return _pancakes[i]; };
+ const Factor & cavitydist (size_t i) const { return _cavitydists[i]; };
+ size_t nr_iI() const { return _iI.size(); };
+
+ void clamp( const Var &n, size_t i ) { assert( 0 == 1 ); }
+ void undoProbs( const VarSet &ns ) { assert( 0 == 1 ); }
+ void saveProbs( const VarSet &ns ) { assert( 0 == 1 ); }
+ void makeFactorCavity(size_t I) { assert( 0 == 1 ); }
+ virtual void makeCavity(const Var & n) { assert( 0 == 1 ); }
+ bool checkProperties();
+};
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+/*=================================================================*
+ * *
+ * This is a MEX-file for MATLAB. *
+ * *
+ * [logZ,q,md] = dai(psi,method,opts); *
+ * *
+ *=================================================================*/
+
+
+#include <iostream>
+#include "matlab.h"
+#include "mex.h"
+#include "alldai.h"
+
+
+using namespace std;
+
+
+/* Input Arguments */
+
+#define PSI_IN prhs[0]
+#define METHOD_IN prhs[1]
+#define OPTS_IN prhs[2]
+#define NR_IN 3
+
+
+/* Output Arguments */
+
+#define LOGZ_OUT plhs[0]
+#define Q_OUT plhs[1]
+#define MD_OUT plhs[2]
+#define NR_OUT 3
+
+
+void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] )
+{
+ size_t buflen;
+
+ /* Check for proper number of arguments */
+ if( (nrhs != NR_IN) || (nlhs != NR_OUT) ) {
+ mexErrMsgTxt("Usage: [logZ,q,md] = dai(psi,method,opts)\n\n"
+ "\n"
+ "INPUT: psi = linear cell array containing the factors \n"
+ " psi{i} should be a structure with a Member field\n"
+ " and a P field, like a CPTAB).\n"
+ " method = name of the method (see README)\n"
+ " opts = string of options (see README)\n"
+ "\n"
+ "OUTPUT: logZ = approximation of the logarithm of the partition sum.\n"
+ " q = linear cell array containing all final beliefs.\n"
+ " md = maxdiff (final linf-dist between new and old single node beliefs).\n");
+ }
+
+ char *method;
+ char *opts;
+
+
+ // Get psi and construct factorgraph
+ vector<Factor> factors = mx2Factors(PSI_IN, 0);
+ FactorGraph fg(factors);
+ long nr_v = fg.nrVars();
+
+ // Get method
+ buflen = mxGetN( METHOD_IN ) + 1;
+ method = (char *)mxCalloc( buflen, sizeof(char) );
+ mxGetString( METHOD_IN, method, buflen );
+
+ // Get options string
+ buflen = mxGetN( OPTS_IN ) + 1;
+ opts = (char *)mxCalloc( buflen, sizeof(char) );
+ mxGetString( OPTS_IN, opts, buflen );
+ // Convert to options object props
+ stringstream ss;
+ ss << opts;
+ Properties props;
+ ss >> props;
+
+ // Construct InfAlg object, init and run
+ InfAlg *obj = newInfAlg( method, fg, props );
+ obj->init();
+ obj->run();
+
+
+ // Save logZ
+ double logZ = real( obj->logZ() );
+
+ // Save maxdiff
+ double maxdiff = obj->MaxDiff();
+
+
+ // Hand over results to MATLAB
+ LOGZ_OUT = mxCreateDoubleMatrix(1,1,mxREAL);
+ *(mxGetPr(LOGZ_OUT)) = logZ;
+
+ Q_OUT = Factors2mx(obj->beliefs());
+
+ MD_OUT = mxCreateDoubleMatrix(1,1,mxREAL);
+ *(mxGetPr(MD_OUT)) = maxdiff;
+
+
+ return;
+}
--- /dev/null
+% [logZ,q,md] = dai (psi, method, opts)
+%
+% INPUT: psi = linear cell array containing the factors
+% psi{i} should be a structure with a Member field
+% and a P field, like a CPTAB).
+% method = name of the method (see README)
+% opts = string of options (see README)
+%
+% OUTPUT: logZ = approximation of the logarithm of the partition sum.
+% q = linear cell array containing all final beliefs.
+% md = maxdiff (final linf-dist between new and old single node beliefs).
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+/*=================================================================*
+ * *
+ * This is a MEX-file for MATLAB. *
+ * *
+ * N = dai_potstrength(psi,i,j); *
+ * *
+ *=================================================================*/
+
+
+#include <iostream>
+#include "mex.h"
+#include "matlab.h"
+#include "../factor.h"
+
+
+using namespace std;
+
+
+/* Input Arguments */
+
+#define PSI_IN prhs[0]
+#define I_IN prhs[1]
+#define J_IN prhs[2]
+#define NR_IN 3
+
+
+/* Output Arguments */
+
+#define N_OUT plhs[0]
+#define NR_OUT 1
+
+
+void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] )
+{
+ long ilabel, jlabel;
+
+ // Check for proper number of arguments
+ if ((nrhs != NR_IN) || (nlhs != NR_OUT)) {
+ mexErrMsgTxt("Usage: N = dai_potstrength(psi,i,j);\n\n"
+ "\n"
+ "INPUT: psi = structure with a Member field and a P field, like a CPTAB.\n"
+ " i = label of a variable in psi.\n"
+ " j = label of another variable in psi.\n"
+ "\n"
+ "OUTPUT: N = strength of psi in direction i->j.\n");
+ }
+
+ // Get input parameters
+ Factor psi = mx2Factor(PSI_IN);
+ ilabel = (long)*mxGetPr(I_IN);
+ jlabel = (long)*mxGetPr(J_IN);
+
+ // Find variable in psi with label ilabel
+ Var i;
+ for( VarSet::const_iterator n = psi.vars().begin(); n != psi.vars().end(); n++ )
+ if( n->label() == ilabel ) {
+ i = *n;
+ break;
+ }
+ assert( i.label() == ilabel );
+
+ // Find variable in psi with label jlabel
+ Var j;
+ for( VarSet::const_iterator n = psi.vars().begin(); n != psi.vars().end(); n++ )
+ if( n->label() == jlabel ) {
+ j = *n;
+ break;
+ }
+ assert( j.label() == jlabel );
+
+ // Calculate N(psi,i,j);
+ double N = psi.strength( i, j );
+
+ // Hand over result to MATLAB
+ N_OUT = mxCreateDoubleMatrix(1,1,mxREAL);
+ *(mxGetPr(N_OUT)) = N;
+
+
+ return;
+}
--- /dev/null
+% N = dai_potstrength (psi, i, j)
+%
+% INPUT: psi = structure with a Member field and a P field, like a CPTAB.
+% i = label of a variable in psi.
+% j = label of another variable in psi.
+%
+% OUTPUT: N = strength of psi in direction i->j.
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+/*=================================================================*
+ * *
+ * This is a MEX-file for MATLAB. *
+ * *
+ * [psi] = dai_readfg(filename); *
+ * *
+ *=================================================================*/
+
+
+#include <iostream>
+#include "mex.h"
+#include "matlab.h"
+#include "factorgraph.h"
+
+
+using namespace std;
+
+
+/* Input Arguments */
+
+#define FILENAME_IN prhs[0]
+#define NR_IN 1
+
+
+/* Output Arguments */
+
+#define PSI_OUT plhs[0]
+#define NR_OUT 1
+
+
+void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] )
+{
+ char *filename;
+
+
+ // Check for proper number of arguments
+ if ((nrhs != NR_IN) || (nlhs != NR_OUT)) {
+ mexErrMsgTxt("Usage: [psi] = dai_readfg(filename);\n\n"
+ "\n"
+ "INPUT: filename = filename of a .fg file\n"
+ "\n"
+ "OUTPUT: psi = linear cell array containing the factors\n"
+ " (psi{i} is a structure with a Member field\n"
+ " and a P field, like a CPTAB).\n");
+ }
+
+ // Get input parameters
+ size_t buflen;
+ buflen = mxGetN( FILENAME_IN ) + 1;
+ filename = (char *)mxCalloc( buflen, sizeof(char) );
+ mxGetString( FILENAME_IN, filename, buflen );
+
+
+ // Read factorgraph
+ FactorGraph fg;
+ if( fg.ReadFromFile( filename ) ) {
+ mexErrMsgTxt("dai_readfg: error reading file\n");
+ }
+
+
+ // Save factors
+ vector<Factor> psi;
+ for( size_t I = 0; I < fg.nrFactors(); I++ )
+ psi.push_back(fg.factor(I));
+
+
+ // Hand over results to MATLAB
+ PSI_OUT = Factors2mx(psi);
+
+
+ return;
+}
--- /dev/null
+% [psi] = dai_readfg (filename)
+%
+% INPUT: filename = filename of a .fg file
+%
+% OUTPUT: psi = linear cell array containing the factors
+% (psi{i} is a structure with a Member field
+% and a P field, like a CPTAB).
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+/*=================================================================*
+ * *
+ * This is a MEX-file for MATLAB. *
+ * *
+ * [psi_out] = dai_removeshortloops(psi_in); *
+ * *
+ *=================================================================*/
+
+
+#include <iostream>
+#include "mex.h"
+#include "matlab.h"
+#include "factorgraph.h"
+
+
+using namespace std;
+
+
+/* Input Arguments */
+
+#define PSI_IN prhs[0]
+#define NR_IN 1
+
+
+/* Output Arguments */
+
+#define PSI_OUT plhs[0]
+#define NR_OUT 1
+
+
+void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] ) {
+ // Check for proper number of arguments
+ if ((nrhs != NR_IN) || (nlhs != NR_OUT)) {
+ mexErrMsgTxt("Usage: [psi_out] = dai_removeshortloops(psi_in);\n\n"
+ "\n"
+ "INPUT: psi_in = linear cell array containing the factors\n"
+ " (psi{i} is a structure with a Member field\n"
+ " and a P field, like a CPTAB).\n"
+ "\n"
+ "OUTPUT: psi_out = linear cell array containing the factors of psi_in,\n"
+ " where factors have been merged to prevent short loops\n"
+ " of length 4 in the factor graph (i.e. loops of type\n"
+ " var1-factor1-var2-factor2-var1),\n");
+ }
+
+ // Get the factors from PSI_IN
+ vector<Factor> psi = mx2Factors(PSI_IN, 0);
+
+ // Remove the short loops
+ RemoveShortLoops(psi);
+
+ // Hand over results to MATLAB
+ PSI_OUT = Factors2mx(psi);
+
+
+ return;
+}
--- /dev/null
+% [psi_out] = dai_removeshortloops (psi_in)
+%
+% INPUT: psi_in = linear cell array containing the factors
+% (psi{i} is a structure with a Member field
+% and a P field, like a CPTAB).
+%
+% OUTPUT: psi_out = linear cell array containing the factors of psi_in,
+% where factors have been merged to prevent short loops
+% of length 4 in the factor graph (i.e. loops of type
+% var1-factor1-var2-factor2-var1).
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+/*=================================================================*
+ * *
+ * This is a MEX-file for MATLAB. *
+ * *
+ * dai_writefg(psi, filename); *
+ * *
+ *=================================================================*/
+
+
+#include <iostream>
+#include "mex.h"
+#include "matlab.h"
+#include "factorgraph.h"
+
+
+using namespace std;
+
+
+/* Input Arguments */
+
+#define PSI_IN prhs[0]
+#define FILENAME_IN prhs[1]
+#define NR_IN 2
+
+
+/* Output Arguments */
+
+#define NR_OUT 0
+
+
+void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] )
+{
+ char *filename;
+
+
+ // Check for proper number of arguments
+ if ((nrhs != NR_IN) || (nlhs != NR_OUT)) {
+ mexErrMsgTxt("Usage: dai_writefg(psi,filename);\n\n"
+ "\n"
+ "INPUT: psi = linear cell array containing the factors\n"
+ " (psi{i} should be a structure with a Member field\n"
+ " and a P field, like a CPTAB).\n"
+ " filename = filename of a .fg file\n");
+ }
+
+ // Get input parameters
+ vector<Factor> factors = mx2Factors(PSI_IN,0);
+
+ size_t buflen;
+ buflen = mxGetN( FILENAME_IN ) + 1;
+ filename = (char *)mxCalloc( buflen, sizeof(char) );
+ mxGetString( FILENAME_IN, filename, buflen );
+
+ // Construct factorgraph
+ FactorGraph fg(factors);
+ long nr_v = fg.nrVars();
+ long nr_f = fg.nrFactors();
+
+ if( fg.WriteToFile( filename ) ) {
+ mexErrMsgTxt("dai_writefg: error reading file\n");
+ }
+
+ return;
+}
--- /dev/null
+% dai_writefg(psi,filename)
+%
+% INPUT: psi = linear cell array containing the factors
+% (psi{i} should be a structure with a Member field
+% and a P field, like a CPTAB).
+% filename = filename of a .fg file
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#include <iostream>
+#include "matlab.h"
+
+
+using namespace std;
+
+
+/* Convert vector<Factor> structure to a cell vector of CPTAB-like structs */
+mxArray *Factors2mx(const vector<Factor> &Ps) {
+ size_t nr = Ps.size();
+
+ mxArray *psi = mxCreateCellMatrix(nr,1);
+
+ const char *fieldnames[2];
+ fieldnames[0] = "Member";
+ fieldnames[1] = "P";
+
+ size_t I_ind = 0;
+ for( vector<Factor>::const_iterator I = Ps.begin(); I != Ps.end(); I++, I_ind++ ) {
+ mxArray *Bi = mxCreateStructMatrix(1,1,2,fieldnames);
+
+ mxArray *BiMember = mxCreateDoubleMatrix(1,I->vars().size(),mxREAL);
+ double *BiMember_data = mxGetPr(BiMember);
+ size_t i = 0;
+ vector<mwSize> dims;
+ for( VarSet::iterator j = I->vars().begin(); j != I->vars().end(); j++,i++ ) {
+ BiMember_data[i] = j->label();
+ dims.push_back( j->states() );
+ }
+
+ // mxArray *BiP = mxCreateDoubleMatrix(I->stateSpace(),1,mxREAL);
+ mxArray *BiP = mxCreateNumericArray(I->vars().size(), &(*(dims.begin())), mxDOUBLE_CLASS, mxREAL);
+ double *BiP_data = mxGetPr(BiP);
+ for( size_t j = 0; j < I->stateSpace(); j++ )
+ BiP_data[j] = (*I)[j];
+
+ mxSetField(Bi,0,"Member",BiMember);
+ mxSetField(Bi,0,"P",BiP);
+
+ mxSetCell(psi, I_ind, Bi);
+ }
+ return( psi );
+}
+
+
+/* Convert cell vector of CPTAB-like structs to vector<Factor> */
+vector<Factor> mx2Factors(const mxArray *psi, long verbose) {
+ set<Var> vars;
+ vector<Factor> factors;
+
+ int n1 = mxGetM(psi);
+ int n2 = mxGetN(psi);
+ if( n2 != 1 && n1 != 1 )
+ mexErrMsgTxt("psi should be a Nx1 or 1xN cell matrix.");
+ size_t nr_f = n1;
+ if( n1 == 1 )
+ nr_f = n2;
+
+ // interpret psi, linear cell array of cptabs
+ for( size_t cellind = 0; cellind < nr_f; cellind++ ) {
+ if( verbose >= 3 )
+ cout << "reading factor " << cellind << ": " << endl;
+ mxArray *cell = mxGetCell(psi, cellind);
+ mxArray *mx_member = mxGetField(cell, 0, "Member");
+ size_t nr_mem = mxGetN(mx_member);
+ double *members = mxGetPr(mx_member);
+ const mwSize *dims = mxGetDimensions(mxGetField(cell,0,"P"));
+ double *factordata = mxGetPr(mxGetField(cell, 0, "P"));
+
+ // add variables
+ VarSet factorvars;
+ vector<long> labels(nr_mem,0);
+ if( verbose >= 3 )
+ cout << " vars: ";
+ for( size_t mi = 0; mi < nr_mem; mi++ ) {
+ labels[mi] = (long)members[mi];
+ if( verbose >= 3 )
+ cout << labels[mi] << "(" << dims[mi] << ") ";
+ vars.insert( Var(labels[mi], dims[mi]) );
+ factorvars.insert( Var(labels[mi], dims[mi]) );
+ }
+ factors.push_back(Factor(factorvars));
+
+ // calculate permutation matrix
+ vector<size_t> perm(nr_mem,0);
+ VarSet::iterator j = factorvars.begin();
+ for( size_t mi = 0; mi < nr_mem; mi++,j++ ) {
+ long gezocht = j->label();
+ vector<long>::iterator piet = find(labels.begin(),labels.end(),gezocht);
+ perm[mi] = piet - labels.begin();
+ }
+
+ if( verbose >= 3 ) {
+ cout << endl << " perm: ";
+ for( vector<size_t>::iterator r=perm.begin(); r!=perm.end(); r++ )
+ cout << *r << " ";
+ cout << endl;
+ }
+
+ // read Factor
+ vector<size_t> di(nr_mem,0);
+ vector<size_t> pdi(nr_mem,0);
+ for( size_t k = 0; k < nr_mem; k++ ) {
+ di[k] = dims[k];
+ pdi[k] = dims[perm[k]];
+ }
+ multind mi(di);
+ multind pmi(pdi);
+ if( verbose >= 3 ) {
+ cout << " mi.max(): " << mi.max() << endl;
+ cout << " ";
+ for( size_t k = 0; k < nr_mem; k++ )
+ cout << labels[k] << " ";
+ cout << " ";
+ for( size_t k = 0; k < nr_mem; k++ )
+ cout << labels[perm[k]] << " ";
+ cout << endl;
+ }
+ for( size_t li = 0; li < mi.max(); li++ ) {
+ vector<size_t> vi = mi.vi(li);
+ vector<size_t> pvi(vi.size(),0);
+ for( size_t k = 0; k < vi.size(); k++ )
+ pvi[k] = vi[perm[k]];
+ size_t pli = pmi.li(pvi);
+ if( verbose >= 3 ) {
+ cout << " " << li << ": ";
+ for( vector<size_t>::iterator r=vi.begin(); r!=vi.end(); r++)
+ cout << *r << " ";
+ cout << "-> ";
+ for( vector<size_t>::iterator r=pvi.begin(); r!=pvi.end(); r++)
+ cout << *r << " ";
+ cout << ": " << pli << endl;
+ }
+ factors.back()[pli] = factordata[li];
+ }
+ }
+
+ if( verbose >= 3 ) {
+ for(vector<Factor>::const_iterator I=factors.begin(); I!=factors.end(); I++ )
+ cout << *I << endl;
+ }
+
+ return( factors );
+}
+
+
+/* Convert CPTAB-like struct to Factor */
+Factor mx2Factor(const mxArray *psi) {
+ mxArray *mx_member = mxGetField(psi, 0, "Member");
+ size_t nr_mem = mxGetN(mx_member);
+ double *members = mxGetPr(mx_member);
+ const mwSize *dims = mxGetDimensions(mxGetField(psi,0,"P"));
+ double *factordata = mxGetPr(mxGetField(psi, 0, "P"));
+
+ // add variables
+ VarSet vars;
+ vector<long> labels(nr_mem,0);
+ for( size_t mi = 0; mi < nr_mem; mi++ ) {
+ labels[mi] = (long)members[mi];
+ vars.insert( Var(labels[mi], dims[mi]) );
+ }
+ Factor factor(vars);
+
+ // calculate permutation matrix
+ vector<size_t> perm(nr_mem,0);
+ VarSet::iterator j = vars.begin();
+ for( size_t mi = 0; mi < nr_mem; mi++,j++ ) {
+ long gezocht = j->label();
+ vector<long>::iterator piet = find(labels.begin(),labels.end(),gezocht);
+ perm[mi] = piet - labels.begin();
+ }
+
+ // read Factor
+ vector<size_t> di(nr_mem,0);
+ vector<size_t> pdi(nr_mem,0);
+ for( size_t k = 0; k < nr_mem; k++ ) {
+ di[k] = dims[k];
+ pdi[k] = dims[perm[k]];
+ }
+ multind mi(di);
+ multind pmi(pdi);
+ for( size_t li = 0; li < mi.max(); li++ ) {
+ vector<size_t> vi = mi.vi(li);
+ vector<size_t> pvi(vi.size(),0);
+ for( size_t k = 0; k < vi.size(); k++ )
+ pvi[k] = vi[perm[k]];
+ size_t pli = pmi.li(pvi);
+ factor[pli] = factordata[li];
+ }
+
+ return( factor );
+}
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#include "mex.h"
+#include "../factor.h"
+
+
+/* Convert vector<Factor> structure to a cell vector of CPTAB-like structs */
+mxArray *Factors2mx(const std::vector<Factor> &Ps);
+
+/* Convert cell vector of CPTAB-like structs to vector<Factor> */
+std::vector<Factor> mx2Factors(const mxArray *psi, long verbose);
+
+/* Convert CPTAB-like struct to Factor */
+Factor mx2Factor(const mxArray *psi);
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#include <iostream>
+#include <sstream>
+#include <map>
+#include <set>
+#include "mf.h"
+#include "diffs.h"
+#include "util.h"
+
+
+using namespace std;
+
+
+const char *MF::Name = "MF";
+
+
+bool MF::checkProperties() {
+ if( !HasProperty("tol") )
+ return false;
+ if (!HasProperty("maxiter") )
+ return false;
+ if (!HasProperty("verbose") )
+ return false;
+
+ ConvertPropertyTo<double>("tol");
+ ConvertPropertyTo<size_t>("maxiter");
+ ConvertPropertyTo<size_t>("verbose");
+
+ return true;
+}
+
+
+void MF::Regenerate() {
+ DAIAlgFG::Regenerate();
+
+ // clear beliefs
+ _beliefs.clear();
+ _beliefs.reserve( nrVars() );
+
+ // create beliefs
+ for( vector<Var>::const_iterator i = vars().begin(); i != vars().end(); i++ )
+ _beliefs.push_back(Factor(*i));
+}
+
+
+string MF::identify() const {
+ stringstream result (stringstream::out);
+ result << Name << GetProperties();
+ return result.str();
+}
+
+
+void MF::init() {
+ assert( checkProperties() );
+
+ for( vector<Factor>::iterator qi = _beliefs.begin(); qi != _beliefs.end(); qi++ )
+ qi->fill(1.0);
+}
+
+
+double MF::run() {
+ clock_t tic = toc();
+
+ if( Verbose() >= 1 )
+ cout << "Starting " << identify() << "...";
+
+ size_t pass_size = _beliefs.size();
+ Diffs diffs(pass_size * 3, 1.0);
+
+ size_t t=0;
+ for( t=0; t < (MaxIter()*pass_size) && diffs.max() > Tol(); t++ ) {
+ // choose random Var i
+ size_t i = (size_t) (nrVars() * rnd_uniform());
+
+ Factor jan;
+ Factor piet;
+ for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); I++ ) {
+
+ Factor henk;
+ for( _nb_cit j = nb2(*I).begin(); j != nb2(*I).end(); j++ ) // for all j in I \ i
+ if( *j != i )
+ henk *= _beliefs[*j];
+ piet = factor(*I).log0();
+ piet *= henk;
+ piet = piet.part_sum(var(i));
+ piet = piet.exp();
+ jan *= piet;
+ }
+
+ jan.normalize( _normtype );
+
+ if( jan.hasNaNs() ) {
+ cout << "MF::run(): ERROR: jan has NaNs!" << endl;
+ return NAN;
+ }
+
+ diffs.push( dist( jan, _beliefs[i], Prob::DISTLINF ) );
+
+ _beliefs[i] = jan;
+ }
+
+ updateMaxDiff( diffs.max() );
+
+ if( Verbose() >= 1 ) {
+ if( diffs.max() > Tol() ) {
+ if( Verbose() == 1 )
+ cout << endl;
+ cout << "MF::run: WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.max() << endl;
+ } else {
+ if( Verbose() >= 2 )
+ cout << "MF::run: ";
+ cout << "converged in " << t / pass_size << " passes (" << toc() - tic << " clocks)." << endl;
+ }
+ }
+
+ return diffs.max();
+}
+
+
+Factor MF::belief1 (size_t i) const {
+ Factor piet;
+ piet = _beliefs[i];
+ piet.normalize( Prob::NORMPROB );
+ return(piet);
+}
+
+
+Factor MF::belief (const VarSet &ns) const {
+ if( ns.size() == 1 )
+ return belief( *(ns.begin()) );
+ else {
+ assert( ns.size() == 1 );
+ return Factor();
+ }
+}
+
+
+Factor MF::belief (const Var &n) const {
+ return( belief1( findVar( n) ) );
+}
+
+
+vector<Factor> MF::beliefs() const {
+ vector<Factor> result;
+ for( size_t i = 0; i < nrVars(); i++ )
+ result.push_back( belief1(i) );
+ return result;
+}
+
+
+Complex MF::logZ() const {
+ Complex sum = 0.0;
+
+ for(size_t i=0; i < nrVars(); i++ )
+ sum -= belief1(i).entropy();
+ for(size_t I=0; I < nrFactors(); I++ ) {
+ Factor henk;
+ for( _nb_cit j = nb2(I).begin(); j != nb2(I).end(); j++ ) // for all j in I
+ henk *= _beliefs[*j];
+ henk.normalize( Prob::NORMPROB );
+ Factor piet;
+ piet = factor(I).log0();
+ piet *= henk;
+ sum -= Complex( piet.totalSum() );
+ }
+
+ return -sum;
+}
+
+
+void MF::init( const VarSet &ns ) {
+ for( size_t i = 0; i < nrVars(); i++ ) {
+ if( ns && var(i) )
+ _beliefs[i].fill( 1.0 );
+ }
+}
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __MF_H__
+#define __MF_H__
+
+
+#include "daialg.h"
+#include "factorgraph.h"
+
+using namespace std;
+
+
+class MF : public DAIAlgFG {
+ protected:
+ vector<Factor> _beliefs;
+
+ public:
+ // default constructor
+ MF() : DAIAlgFG(), _beliefs() {};
+ // copy constructor
+ MF(const MF & x) : DAIAlgFG(x), _beliefs(x._beliefs) {};
+ MF* clone() const { return new MF(*this); }
+ // construct MF object from FactorGraph
+ MF(const FactorGraph & fg, const Properties &opts) : DAIAlgFG(fg, opts) {
+ assert( checkProperties() );
+ Regenerate();
+ }
+ // assignment operator
+ MF & operator=(const MF & x) {
+ if(this!=&x) {
+ DAIAlgFG::operator=(x);
+ _beliefs = x._beliefs;
+ }
+ return *this;
+ }
+
+ static const char *Name;
+ string identify() const;
+ void Regenerate();
+ void init();
+ double run();
+ Factor belief1 (size_t i) const;
+ Factor belief (const Var &n) const;
+ Factor belief (const VarSet &ns) const;
+ vector<Factor> beliefs() const;
+ Complex logZ() const;
+
+ void init( const VarSet &ns );
+ void undoProbs( const VarSet &ns ) { FactorGraph::undoProbs(ns); init(ns); }
+ bool checkProperties();
+};
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#include "mr.h"
+#include "bp.h"
+#include "jtree.h"
+#include <stdio.h>
+#include <time.h>
+#include <math.h>
+#include <stdlib.h>
+#include "util.h"
+#include "diffs.h"
+
+
+using namespace std;
+
+
+const char *MR::Name = "MR";
+
+
+bool MR::checkProperties() {
+ if( !HasProperty("updates") )
+ return false;
+ if( !HasProperty("inits") )
+ return false;
+ if( !HasProperty("verbose") )
+ return false;
+ if( !HasProperty("tol") )
+ return false;
+
+ ConvertPropertyTo<UpdateType>("updates");
+ ConvertPropertyTo<InitType>("inits");
+ ConvertPropertyTo<size_t>("verbose");
+ ConvertPropertyTo<double>("tol");
+
+ return true;
+}
+
+
+// init N, con, nb, tJ, theta
+void MR::init(size_t _N, double *_w, double *_th) {
+ size_t i,j;
+
+ N = _N;
+
+ con.resize(N);
+ nb.resize(N);
+ tJ.resize(N);
+ for(i=0; i<N; i++ ) {
+ nb[i].resize(kmax);
+ tJ[i].resize(kmax);
+ con[i]=0;
+ for(j=0; j<N; j++ )
+ if( _w[i*N+j] != 0.0 ) {
+ nb[i][con[i]] = j;
+ tJ[i][con[i]] = tanh(_w[i*N+j]);
+ con[i]++;
+ }
+ }
+
+ theta.resize(N);
+ for(i=0; i<N; i++)
+ theta[i] = _th[i];
+}
+
+
+// calculate cors
+double MR::init_cor_resp() {
+ size_t j,k,l, runx,i2;
+ double variab1, variab2;
+ double md, maxdev;
+ double thbJsite[kmax];
+ double xinter;
+ double rinter;
+ double res[kmax];
+ size_t s2;
+ size_t flag;
+ size_t concav;
+ size_t runs = 3000;
+ double eps = 0.2;
+ size_t cavity;
+
+ vector<vector<double> > tJ_org;
+ vector<vector<size_t> > nb_org;
+ vector<size_t> con_org;
+ vector<double> theta_org;
+
+ vector<double> xfield(N*kmax,0.0);
+ vector<double> rfield(N*kmax,0.0);
+ vector<double> Hfield(N,0.0);
+ vector<double> devs(N*kmax,0.0);
+ vector<double> devs2(N*kmax,0.0);
+ vector<double> dev(N,0.0);
+ vector<double> avmag(N,0.0);
+
+ // save original tJ, nb
+ nb_org = nb;
+ tJ_org = tJ;
+ con_org = con;
+ theta_org = theta;
+
+ maxdev = 0.0;
+ for(cavity=0; cavity<N; cavity++){ // for each spin to be removed
+ con = con_org;
+ concav=con[cavity];
+
+ nb = nb_org;
+ tJ = tJ_org;
+
+ // Adapt the graph variables nb[], tJ[] and con[]
+ for(size_t i=0; i<con[cavity]; i++) {
+ size_t ij = nb[cavity][i];
+ flag=0;
+ j=0;
+ do{
+ if(nb[ij][j]==cavity){
+ while(j<(con[ij]-1)){
+ nb[ij][j]=nb[ij][j+1];
+ tJ[ij][j] = tJ[ij][j+1];
+ j++;
+ }
+ flag=1;
+ }
+ j++;
+ } while(flag==0);
+ }
+ for(size_t i=0; i<con[cavity]; i++)
+ con[nb[cavity][i]]--;
+ con[cavity] = 0;
+
+
+ // Do everything starting from the new graph********
+
+ makekindex();
+ theta = theta_org;
+
+ for(size_t i=0; i<kmax*N; i++)
+ xfield[i] = 3.0*(2*rnd_uniform()-1.);
+
+ for(i2=0; i2<concav; i2++){ // Subsequently apply a field to each cavity spin ****************
+
+ s2 = nb[cavity][i2]; // identify the index of the cavity spin
+ for(size_t i=0; i<con[s2]; i++)
+ rfield[kmax*s2+i] = 1.;
+
+ runx=0;
+ do { // From here start the response and belief propagation
+ runx++;
+ md=0.0;
+ for(k=0; k<N; k++){
+ if(k!=cavity) {
+ for(size_t i=0; i<con[k]; i++)
+ thbJsite[i] = tJ[k][i];
+ for(l=0; l<con[k]; l++){
+ xinter = 1.;
+ rinter = 0.;
+ if(k==s2) rinter += 1.;
+ for(j=0; j<con[k]; j++)
+ if(j!=l){
+ variab2 = tanh(xfield[kmax*nb[k][j]+kindex[k][j]]);
+ variab1 = thbJsite[j]*variab2;
+ xinter *= (1.+variab1)/(1.-variab1);
+
+ rinter += thbJsite[j]*rfield[kmax*nb[k][j]+kindex[k][j]]*(1-variab2*variab2)/(1-variab1*variab1);
+ }
+
+ variab1 = 0.5*log(xinter);
+ xinter = variab1 + theta[k];
+ devs[kmax*k+l] = xinter-xfield[kmax*k+l];
+ xfield[kmax*k+l] = xfield[kmax*k+l]+devs[kmax*k+l]*eps;
+ if( fabs(devs[kmax*k+l]) > md )
+ md = fabs(devs[kmax*k+l]);
+
+ devs2[kmax*k+l] = rinter-rfield[kmax*k+l];
+ rfield[kmax*k+l] = rfield[kmax*k+l]+devs2[kmax*k+l]*eps;
+ if( fabs(devs2[kmax*k+l]) > md )
+ md = fabs(devs2[kmax*k+l]);
+ }
+ }
+ }
+ } while((md > Tol())&&(runx<runs)); // Precision condition reached -> BP and RP finished
+ if(runx==runs)
+ if( Verbose() >= 2 )
+ cout << "init_cor_resp: Convergence not reached (md=" << md << ")..." << endl;
+ if(md > maxdev)
+ maxdev = md;
+
+ // compute the observables (i.e. magnetizations and responses)******
+
+ for(size_t i=0; i<concav; i++){
+ rinter = 0.;
+ xinter = 1.;
+ if(i!=i2)
+ for(j=0; j<con[nb[cavity][i]]; j++){
+ variab2 = tanh(xfield[kmax*nb[nb[cavity][i]][j]+kindex[nb[cavity][i]][j]]);
+ variab1 = tJ[nb[cavity][i]][j]*variab2;
+ rinter += tJ[nb[cavity][i]][j]*rfield[kmax*nb[nb[cavity][i]][j]+kindex[nb[cavity][i]][j]]*(1-variab2*variab2)/(1-variab1*variab1);
+ xinter *= (1.+variab1)/(1.-variab1);
+ }
+ xinter = tanh(0.5*log(xinter)+theta[nb[cavity][i]]);
+ res[i] = rinter*(1-xinter*xinter);
+ }
+
+ // *******************
+
+ for(size_t i=0; i<concav; i++)
+ if(nb[cavity][i]!=s2)
+ // if(i!=i2)
+ cors[cavity][i2][i] = res[i];
+ else
+ cors[cavity][i2][i] = 0;
+ } // close for i2 = 0...concav
+ }
+
+ // restore nb, tJ, con
+ tJ = tJ_org;
+ nb = nb_org;
+ con = con_org;
+ theta = theta_org;
+
+ return maxdev;
+}
+
+
+double MR::T(size_t i, sub_nb A) {
+ // i is a variable index
+ // A is a subset of nb[i]
+ //
+ // calculate T{(i)}_A as defined in Rizzo&Montanari e-print (2.17)
+
+ sub_nb _nbi_min_A(con[i]);
+ for( size_t __j = 0; __j < A.size(); __j++ )
+ _nbi_min_A -= A[__j];
+
+ double res = theta[i];
+ for( size_t __j = 0; __j < _nbi_min_A.size(); __j++ ) {
+ size_t _j = _nbi_min_A[__j];
+ res += atanh(tJ[i][_j] * M[i][_j]);
+ }
+ return tanh(res);
+}
+
+
+double MR::T(size_t i, size_t _j) {
+ sub_nb j(con[i]);
+ j.clear();
+ j += _j;
+ return T(i,j);
+}
+
+
+double MR::Omega(size_t i, size_t _j, size_t _l) {
+ sub_nb jl(con[i]);
+ jl.clear();
+ jl += _j;
+ jl += _l;
+ double Tijl = T(i,jl);
+ return Tijl / (1.0 + tJ[i][_l] * M[i][_l] * Tijl);
+}
+
+
+double MR::Gamma(size_t i, size_t _j, size_t _l1, size_t _l2) {
+ sub_nb jll(con[i]);
+ jll.clear();
+ jll += _j;
+ double Tij = T(i,jll);
+ jll += _l1;
+ jll += _l2;
+ double Tijll = T(i,jll);
+
+ return (Tijll - Tij) / (1.0 + tJ[i][_l1] * tJ[i][_l2] * M[i][_l1] * M[i][_l2] + tJ[i][_l1] * M[i][_l1] * Tijll + tJ[i][_l2] * M[i][_l2] * Tijll);
+}
+
+
+double MR::Gamma(size_t i, size_t _l1, size_t _l2) {
+ sub_nb ll(con[i]);
+ ll.clear();
+ double Ti = T(i,ll);
+ ll += _l1;
+ ll += _l2;
+ double Till = T(i,ll);
+
+ return (Till - Ti) / (1.0 + tJ[i][_l1] * tJ[i][_l2] * M[i][_l1] * M[i][_l2] + tJ[i][_l1] * M[i][_l1] * Till + tJ[i][_l2] * M[i][_l2] * Till);
+}
+
+
+double MR::_tJ(size_t i, sub_nb A) {
+ // i is a variable index
+ // A is a subset of nb[i]
+ //
+ // calculate the product of all tJ[i][_j] for _j in A
+
+ size_t Asize = A.size();
+ switch( Asize ) {
+ case 0: return 1.0;
+// case 1: return tJ[i][A[0]];
+// case 2: return tJ[i][A[0]] * tJ[i][A[1]];
+// case 3: return tJ[i][A[0]] * tJ[i][A[1]] * tJ[i][A[2]];
+ default:
+ size_t __j = Asize - 1;
+ size_t _j = A[__j];
+ sub_nb A_j = A - _j;
+ return tJ[i][_j] * _tJ(i, A_j);
+ }
+
+}
+
+
+double MR::appM(size_t i, sub_nb A) {
+ // i is a variable index
+ // A is a subset of nb[i]
+ //
+ // calculate the moment of variables in A from M and cors, neglecting higher order cumulants,
+ // defined as the sum over all partitions of A into subsets of cardinality two at most of the
+ // product of the cumulants (either first order, i.e. M, or second order, i.e. cors) of the
+ // entries of the partitions
+
+ size_t Asize = A.size();
+ switch( Asize ) {
+ case 0: return 1.0;
+// case 1: return M[i][A[0]];
+// case 2: return M[i][A[0]] * M[i][A[1]] + cors[i][A[0]][A[1]];
+// case 3: return M[i][A[0]] * M[i][A[1]] * M[i][A[2]] + M[i][A[0]] * cors[i][A[1]][A[2]] + M[i][A[1]] * cors[i][A[0]][A[2]] + M[i][A[2]] * cors[i][A[0]][A[1]];
+ default:
+ size_t __j = Asize - 1;
+ size_t _j = A[__j];
+ sub_nb A_j = A - _j;
+
+ double result = M[i][_j] * appM(i, A_j);
+ for( size_t __k = 0; __k < __j; __k++ ) {
+ size_t _k = A[__k];
+ result += cors[i][_j][_k] * appM(i,A_j - _k);
+ }
+
+ return result;
+ }
+}
+
+
+void MR::sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd) {
+ // j is a variable index
+ // A is a subset of nb[j]
+
+ // calculate sum over all even/odd subsets B of A of _tJ(j,B) appM(j,B)
+
+ *sum_even = 0.0;
+ *sum_odd = 0.0;
+
+ sub_nb S(A.size());
+ S.clear();
+ do { // for all subsets of A
+
+ // construct subset B of A corresponding to S
+ sub_nb B = A;
+ B.clear();
+ size_t Ssize = S.size();
+ for( size_t bit = 0; bit < Ssize; bit++ )
+ B += A[S[bit]];
+
+ if( S.size() % 2 )
+ *sum_odd += _tJ(j,B) * appM(j,B);
+ else
+ *sum_even += _tJ(j,B) * appM(j,B);
+
+ ++S;
+ } while( !S.empty() );
+}
+
+
+void MR::solvemcav() {
+ double sum_even, sum_odd;
+ double maxdev;
+ size_t maxruns = 1000;
+
+ makekindex();
+ for(size_t i=0; i<N; i++)
+ for(size_t _j=0; _j<con[i]; _j++)
+ M[i][_j]=0.1;
+
+ size_t run=0;
+ do {
+ maxdev=0.0;
+ run++;
+ for(size_t i=0; i<N; i++){ // for all i
+ for(size_t _j=0; _j<con[i]; _j++){ // for all j in N_i
+ size_t _i = kindex[i][_j];
+ size_t j = nb[i][_j];
+ assert( nb[j][_i] == i );
+
+ double newM = 0.0;
+ if( Updates() == UpdateType::FULL ) {
+ // find indices in nb[j] that do not correspond with i
+ sub_nb _nbj_min_i(con[j]);
+ _nbj_min_i -= kindex[i][_j];
+
+ // find indices in nb[i] that do not correspond with j
+ sub_nb _nbi_min_j(con[i]);
+ _nbi_min_j -= _j;
+
+ sum_subs(j, _nbj_min_i, &sum_even, &sum_odd);
+ newM = (tanh(theta[j]) * sum_even + sum_odd) / (sum_even + tanh(theta[j]) * sum_odd);
+
+ sum_subs(i, _nbi_min_j, &sum_even, &sum_odd);
+ double denom = sum_even + tanh(theta[i]) * sum_odd;
+ double numer = 0.0;
+ for(size_t _k=0; _k<con[i]; _k++) if(_k != _j) {
+ sum_subs(i, _nbi_min_j - _k, &sum_even, &sum_odd);
+ numer += tJ[i][_k] * cors[i][_j][_k] * (tanh(theta[i]) * sum_even + sum_odd);
+ }
+ newM -= numer / denom;
+ } else if( Updates() == UpdateType::LINEAR ) {
+ newM = T(j,_i);
+ for(size_t _l=0; _l<con[i]; _l++) if( _l != _j )
+ newM -= Omega(i,_j,_l) * tJ[i][_l] * cors[i][_j][_l];
+ for(size_t _l1=0; _l1<con[j]; _l1++) if( _l1 != _i )
+ for( size_t _l2=_l1+1; _l2<con[j]; _l2++) if( _l2 != _i)
+ newM += Gamma(j,_i,_l1,_l2) * tJ[j][_l1] * tJ[j][_l2] * cors[j][_l1][_l2];
+ }
+
+ double dev = newM - M[i][_j];
+// dev *= 0.02;
+ if( fabs(dev) >= maxdev )
+ maxdev = fabs(dev);
+
+ newM = M[i][_j] + dev;
+ if( fabs(newM) > 1.0 )
+ newM = sign(newM);
+ M[i][_j] = newM;
+ }
+ }
+ } while((maxdev>Tol())&&(run<maxruns));
+
+ updateMaxDiff( maxdev );
+
+ if(run==maxruns){
+ if( Verbose() >= 1 )
+ cout << "solve_mcav: Convergence not reached (maxdev=" << maxdev << ")..." << endl;
+ }
+}
+
+
+void MR::solveM() {
+ for(size_t i=0; i<N; i++) {
+ if( Updates() == UpdateType::FULL ) {
+ // find indices in nb[i]
+ sub_nb _nbi(con[i]);
+
+ // calc numerator1 and denominator1
+ double sum_even, sum_odd;
+ sum_subs(i, _nbi, &sum_even, &sum_odd);
+
+ Mag[i] = (tanh(theta[i]) * sum_even + sum_odd) / (sum_even + tanh(theta[i]) * sum_odd);
+
+ } else if( Updates() == UpdateType::LINEAR ) {
+ sub_nb empty(con[i]);
+ empty.clear();
+ Mag[i] = T(i,empty);
+
+ for(size_t _l1=0; _l1<con[i]; _l1++)
+ for( size_t _l2=_l1+1; _l2<con[i]; _l2++)
+ Mag[i] += Gamma(i,_l1,_l2) * tJ[i][_l1] * tJ[i][_l2] * cors[i][_l1][_l2];
+ }
+ if(fabs(Mag[i])>1.)
+ Mag[i] = sign(Mag[i]);
+ }
+}
+
+
+double MR::init_cor() {
+ for( size_t i = 0; i < nrVars(); i++ ) {
+ vector<Factor> pairq;
+ if( Inits() == InitType::CLAMPING ) {
+ BP bpcav(*this, Properties()("updates",string("SEQMAX"))("tol", string("1e-9"))("maxiter", string("1000UL"))("verbose", string("0UL")));
+ bpcav.makeCavity( var(i) );
+ pairq = calcPairBeliefs( bpcav, delta(var(i)), false );
+ } else if( Inits() == InitType::EXACT ) {
+ JTree jtcav(*this, Properties()("updates",string("HUGIN"))("verbose", string("0UL")) );
+ jtcav.makeCavity( var(i) );
+ pairq = calcPairBeliefs( jtcav, delta(var(i)), false );
+ }
+ for( size_t jk = 0; jk < pairq.size(); jk++ ) {
+ VarSet::const_iterator kit = pairq[jk].vars().begin();
+ size_t j = findVar( *(kit) );
+ size_t k = findVar( *(++kit) );
+ pairq[jk].normalize(Prob::NORMPROB);
+ double cor = (pairq[jk][3] - pairq[jk][2] - pairq[jk][1] + pairq[jk][0]) - (pairq[jk][3] + pairq[jk][2] - pairq[jk][1] - pairq[jk][0]) * (pairq[jk][3] - pairq[jk][2] + pairq[jk][1] - pairq[jk][0]);
+ for( size_t _j = 0; _j < con[i]; _j++ ) if( nb[i][_j] == j )
+ for( size_t _k = 0; _k < con[i]; _k++ ) if( nb[i][_k] == k ) {
+ cors[i][_j][_k] = cor;
+ cors[i][_k][_j] = cor;
+ }
+ }
+ }
+}
+
+
+string MR::identify() const {
+ stringstream result (stringstream::out);
+ result << Name << GetProperties();
+ return result.str();
+}
+
+
+double MR::run() {
+ if( supported ) {
+ if( Verbose() >= 1 )
+ cout << "Starting " << identify() << "...";
+
+ clock_t tic = toc();
+// Diffs diffs(nrVars(), 1.0);
+
+ M.resize(N);
+ for(size_t i=0; i<N; i++)
+ M[i].resize(kmax);
+
+ cors.resize(N);
+ for(size_t i=0; i<N; i++)
+ cors[i].resize(kmax);
+ for(size_t i=0; i<N; i++)
+ for(size_t j=0; j<kmax; j++)
+ cors[i][j].resize(kmax);
+
+ kindex.resize(N);
+ for(size_t i=0; i<N; i++)
+ kindex[i].resize(kmax);
+
+ if( Inits() == InitType::RESPPROP )
+ updateMaxDiff( init_cor_resp() );
+ else if( Inits() == InitType::EXACT )
+ init_cor(); // FIXME no MaxDiff() calculation
+ else if( Inits() == InitType::CLAMPING )
+ init_cor(); // FIXME no MaxDiff() calculation
+
+ solvemcav();
+
+ Mag.resize(N);
+ solveM();
+
+ if( Verbose() >= 1 )
+ cout << "MR needed " << toc() - tic << " clocks." << endl;
+
+ return 0.0;
+ } else
+ return NAN;
+}
+
+
+void MR::makekindex() {
+ for(size_t i=0; i<N; i++)
+ for(size_t j=0; j<con[i]; j++) {
+ size_t ij = nb[i][j]; // ij is the j'th neighbour of spin i
+ size_t k=0;
+ while( nb[ij][k] != i )
+ k++;
+ kindex[i][j] = k; // the j'th neighbour of spin i has spin i as its k'th neighbour
+ }
+}
+
+
+Factor MR::belief( const Var &n ) const {
+ if( supported ) {
+ size_t i = findVar( n );
+
+ double x[2];
+ x[0] = 0.5 - Mag[i] / 2.0;
+ x[1] = 0.5 + Mag[i] / 2.0;
+
+ return Factor( n, x );
+ } else
+ return Factor();
+}
+
+
+vector<Factor> MR::beliefs() const {
+ vector<Factor> result;
+ for( size_t i = 0; i < nrVars(); i++ )
+ result.push_back( belief( var(i) ) );
+ return result;
+}
+
+
+
+MR::MR( const FactorGraph &fg, const Properties &opts ) : DAIAlgFG(fg, opts), supported(true) {
+ // check whether all vars in fg are binary
+ // check whether connectivity is <= kmax
+ for( size_t i = 0; i < fg.nrVars(); i++ )
+ if( (fg.var(i).states() > 2) || (fg.delta(fg.var(i)).size() > kmax) ) {
+ supported = false;
+ break;
+ }
+
+ if( !supported )
+ return;
+
+ // check whether all interactions are pairwise or single
+ for( size_t I = 0; I < fg.nrFactors(); I++ )
+ if( fg.factor(I).vars().size() > 2 ) {
+ supported = false;
+ break;
+ }
+
+ if( !supported )
+ return;
+
+ // create w and th
+ size_t _N = fg.nrVars();
+
+ double *w = new double[_N*_N];
+ double *th = new double[_N];
+
+ for( size_t i = 0; i < _N; i++ ) {
+ th[i] = 0.0;
+ for( size_t j = 0; j < _N; j++ )
+ w[i*_N+j] = 0.0;
+ }
+
+ for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+ const Factor &psi = fg.factor(I);
+ if( psi.vars().size() == 1 ) {
+ size_t i = fg.findVar( *(psi.vars().begin()) );
+ th[i] += 0.5 * log(psi[1] / psi[0]);
+ } else if( psi.vars().size() == 2 ) {
+ size_t i = fg.findVar( *(psi.vars().begin()) );
+ VarSet::const_iterator jit = psi.vars().begin();
+ size_t j = fg.findVar( *(++jit) );
+
+ w[i*_N+j] += 0.25 * log(psi[3] * psi[0] / (psi[2] * psi[1]));
+ w[j*_N+i] += 0.25 * log(psi[3] * psi[0] / (psi[2] * psi[1]));
+
+ th[i] += 0.25 * log(psi[3] / psi[2] * psi[1] / psi[0]);
+ th[j] += 0.25 * log(psi[3] / psi[1] * psi[2] / psi[0]);
+ }
+ }
+
+ init(_N, w, th);
+
+ delete th;
+ delete w;
+}
--- /dev/null
+/* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
+ Radboud University Nijmegen, The Netherlands
+
+ This file is part of libDAI.
+
+ libDAI is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ libDAI is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with libDAI; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+
+#ifndef __MR_H__
+#define __MR_H__
+
+
+#include <vector>
+#include "factorgraph.h"
+#include "daialg.h"
+#include "enum.h"
+
+
+using namespace std;
+
+
+class sub_nb;
+
+
+class MR : public DAIAlgFG {
+ private:
+ bool supported; // is the underlying factor graph supported?
+
+ vector<size_t> con; // con[i] = connectivity of spin i
+ vector<vector<size_t> > nb; // nb[i] are the neighbours of spin i
+ vector<vector<double> > tJ; // tJ[i][_j] is the tanh of the interaction between spin i and its neighbour nb[i][_j]
+ vector<double> theta; // theta[i] is the local field on spin i
+ vector<vector<double> > M; // M[i][_j] is M^{(i)}_j
+ vector<vector<size_t> > kindex; // the _j'th neighbour of spin i has spin i as its kindex[i][_j]'th neighbour
+ vector<vector<vector<double> > > cors;
+
+ static const size_t kmax = 31;
+
+ size_t N;
+
+ vector<double> Mag;
+
+ public:
+ ENUM2(UpdateType,FULL,LINEAR)
+ ENUM3(InitType,RESPPROP,CLAMPING,EXACT)
+
+ UpdateType Updates() const { return GetPropertyAs<UpdateType>("updates"); }
+ InitType Inits() const { return GetPropertyAs<InitType>("inits"); }
+
+ MR( const FactorGraph & fg, const Properties &opts );
+ void init(size_t _N, double *_w, double *_th);
+ void makekindex();
+ void read_files();
+ double init_cor();
+ double init_cor_resp();
+ void solvemcav();
+ void solveM();
+ double run();
+ Factor belief( const Var &n ) const;
+ Factor belief( const VarSet &ns ) const { assert( 0 == 1 ); }
+ vector<Factor> beliefs() const;
+ Complex logZ() const { return NAN; }
+ void init() { assert( checkProperties() );
projects / libdai.git / commitdiff