diff --git a/examples/BayesNet.py b/examples/BayesNet.py
index e99b7ebcaf399d3b322899d15de0985e5ad317d0..710d6d5fc46e7bfe3299027efa049515ceaaf0ac 100755
--- a/examples/BayesNet.py
+++ b/examples/BayesNet.py
@@ -1,10 +1,10 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-from  primo.core import BayesNet
-from  primo.reasoning import DiscreteNode
+from  primo.networks import BayesianNetwork
+from  primo.nodes import DiscreteNode
 import numpy
 
-bn = BayesNet()
+bn = BayesianNetwork()
 burglary = DiscreteNode("Burglary", ["Intruder","Safe"])
 alarm = DiscreteNode("Alarm", ["Ringing", "Silent","Kaputt"])
 earthquake = DiscreteNode("Earthquake", ["Shaking", "Calm"])
diff --git a/examples/BaysianDecisionNetwork.py b/examples/BaysianDecisionNetwork.py
index e827f85d559d8372257cc5ce3259471e531c24a5..4d85b4df3a1240b0d903fac8528974b9a8fca5fd 100644
--- a/examples/BaysianDecisionNetwork.py
+++ b/examples/BaysianDecisionNetwork.py
@@ -1,8 +1,8 @@
-from primo.core import BayesianDecisionNetwork
-from primo.decision import DecisionNode
-from primo.decision import UtilityNode
-from primo.reasoning import DiscreteNode
-from primo.decision.make_decision import MakeDecision
+from primo.networks import BayesianDecisionNetwork
+from primo.nodes import DecisionNode
+from primo.nodes import UtilityNode
+from primo.nodes import DiscreteNode
+from primo.inference.decision import MakeDecision
 
 import numpy
 
diff --git a/examples/CalculateMarginals.py b/examples/CalculateMarginals.py
index b416cf49c94f52cf47db5cdf2d304d5a961a1831..23d834ffde554bc8a18d39973a779eb99f05d2f0 100644
--- a/examples/CalculateMarginals.py
+++ b/examples/CalculateMarginals.py
@@ -6,15 +6,15 @@ the probability of evidence.
 @author: djohn
 """
 
-from  primo.core import BayesNet
-from  primo.reasoning import DiscreteNode
-from primo.reasoning.factorelemination import EasiestFactorElimination
-from primo.reasoning.factorelemination import FactorTreeFactory
+from primo.networks import BayesianNetwork
+from primo.nodes import DiscreteNode
+from primo.inference.factor import EasiestFactorElimination
+from primo.inference.factor import FactorTreeFactory
 import numpy
 
 #==============================================================================
 #Builds the example BayesNet
-bn = BayesNet()
+bn = BayesianNetwork()
 burglary = DiscreteNode("Burglary", ["Intruder","Safe"])
 alarm = DiscreteNode("Alarm", ["Ringing", "Silent"])
 earthquake = DiscreteNode("Earthquake", ["Shaking", "Calm"])
diff --git a/examples/Compare.py b/examples/Compare.py
index 2ae438bf144b1d0b243855bc274a629488c76030..28a58c8bf867750cdba2be97ce7e0cb9d482c5c9 100644
--- a/examples/Compare.py
+++ b/examples/Compare.py
@@ -6,18 +6,17 @@ the probability of evidence.
 @author: djohn, mbaum
 """
 
-from  primo.core import BayesNet
-from  primo.reasoning import DiscreteNode
-from primo.reasoning.density import ProbabilityTable
-#from primo.reasoning.factorelemination import EasiestFactorElimination
-from primo.reasoning.factorelemination import FactorTreeFactory
-from primo.reasoning import MCMC
-from primo.reasoning import EvidenceEqual as EvEq
+from primo.networks import BayesianNetwork
+from primo.nodes import DiscreteNode
+from primo.densities import ProbabilityTable
+from primo.inference.factor import FactorTreeFactory
+from primo.inference.mcmc import MCMC
+from primo.evidence import EvidenceEqual as EvEq
 import numpy
 
 #==============================================================================
 #Builds the example BayesNet
-bn = BayesNet()
+bn = BayesianNetwork()
 burglary = DiscreteNode("Burglary", ["Intruder","Safe"])
 alarm = DiscreteNode("Alarm", ["Ringing", "Silent"])
 earthquake = DiscreteNode("Earthquake", ["Shaking", "Calm"])
@@ -119,7 +118,7 @@ print "EarthquakeFT: " + str(factorTree.calculate_marginal([earthquake]))
 
 
 
-mcmc_ask=MCMC(bn)
+mcmc_ask=MCMC(bn,1000)
 
 print "====MCMC===="
 
diff --git a/examples/DBN_Robot_Localization.py b/examples/DBN_Robot_Localization.py
index d423b3ed56c2c18f537299e381da5a6027d5c1c0..3c6d10d76ac3955c53526e5c0e5e8b43da4f20af 100755
--- a/examples/DBN_Robot_Localization.py
+++ b/examples/DBN_Robot_Localization.py
@@ -1,18 +1,21 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-from primo.core import BayesNet
-from primo.core import DynamicBayesNet
-from primo.core import TwoTBN
-from primo.reasoning import DiscreteNode
-import primo.reasoning.particlebased.ParticleFilterDBN as pf
-import numpy
-from primo.utils import XMLBIF
 import random
 
+import numpy
+
+from primo.networks import BayesianNetwork
+from primo.networks import DynamicBayesianNetwork
+from primo.networks import TwoTBN
+from primo.nodes import DiscreteNode
+import primo.inference.particlefilter as pf
+from primo.io import XMLBIF
+
+
 
 # Construct a DynmaicBayesianNetwork
-dbn = DynamicBayesNet()
-B0 = BayesNet()
+dbn = DynamicBayesianNetwork()
+B0 = BayesianNetwork()
 twoTBN = TwoTBN(XMLBIF.read("Robot_Localization.xmlbif"))
 
 # Configure TwoTBN
diff --git a/examples/EarthquakeNet.py b/examples/EarthquakeNet.py
index 8c866f2aad9242b832a3754ace29577538c2881b..18c851cf1f2ce66a5357c4a3851bb0094ebbe27f 100644
--- a/examples/EarthquakeNet.py
+++ b/examples/EarthquakeNet.py
@@ -6,12 +6,14 @@ This example shows how to create a BayesNet
 @author: djohn
 """
 
-from  primo.core import BayesNet
-from  primo.reasoning import DiscreteNode
 import numpy
 
+from  primo.networks import BayesianNetwork
+from  primo.nodes import DiscreteNode
+
+
 #initialize a new BayesNet
-bn = BayesNet()
+bn = BayesianNetwork()
 
 #create Nodes with Name and the possible values
 burglary = DiscreteNode("Burglary", ["Intruder","Safe"])
@@ -61,6 +63,6 @@ john_calls.set_probability(0.01,[(alarm,"Silent"),(john_calls,"Calling")])
 john_calls.set_probability(0.99,[(alarm,"Silent"),(john_calls,"Not Calling")])
 
 #draws the BayesNet
-#bn.draw()
+bn.draw()
 
 
diff --git a/examples/Gauss.py b/examples/Gauss.py
index 0a63c15ec99b6ae0f690e424ff67ed7a30ccc444..0fcece543cdee05a533f60bbea3e6bd1948b9bcd 100644
--- a/examples/Gauss.py
+++ b/examples/Gauss.py
@@ -1,9 +1,10 @@
 import numpy
 
-from primo.core import BayesNet
-from primo.reasoning import GaussNode
+from primo.networks import BayesianNetwork
+from primo.nodes import ContinuousNodeFactory
 
-age = GaussNode("Plant_age")
+cnf = ContinuousNodeFactory()
+age = cnf.createGaussNode("Plant_age")
 
 #parameterization
 age_b0=numpy.array([4])
diff --git a/examples/MCMC_Continuous.py b/examples/MCMC_Continuous.py
index c184da12e68f86e822d17cd00bfba44ed84b128d..63e5a8a5b716a6b5d70b14c3253505e28c2a98b5 100644
--- a/examples/MCMC_Continuous.py
+++ b/examples/MCMC_Continuous.py
@@ -1,14 +1,15 @@
-from primo.core import BayesNet
-from primo.reasoning import ContinuousNodeFactory
-from primo.reasoning.density import ExponentialParameters
-from primo.reasoning.density import BetaParameters
-from primo.reasoning.density import GaussParameters
-from primo.reasoning.density import NDGauss
-from primo.reasoning import MCMC
-from primo.reasoning.convergence_test import ConvergenceTestSimpleCounting
-
-from primo.reasoning import EvidenceEqual as EvEqual
-from primo.reasoning import EvidenceLower as EvLower
+from primo.networks import BayesianNetwork
+from primo.nodes import ContinuousNode
+from primo.nodes import ContinuousNodeFactory
+from primo.densities import ExponentialParameters
+from primo.densities import BetaParameters
+from primo.densities import GaussParameters
+from primo.densities import NDGauss
+from primo.inference.mcmc import MCMC
+from primo.inference.mcmc import ConvergenceTestSimpleCounting
+
+from primo.evidence import EvidenceEqual as EvEqual
+from primo.evidence import EvidenceLower as EvLower
 
 import numpy
 
@@ -26,7 +27,7 @@ import numpy
 #Construct some simple BayesianNetwork. 
 
 #topology
-bn = BayesNet()
+bn = BayesianNetwork()
 cnf=ContinuousNodeFactory()
 age = cnf.createExponentialNode("Plant_age")
 height = cnf.createGaussNode("Plant_height")
diff --git a/examples/MCMC_Continuous_Plant.py b/examples/MCMC_Continuous_Plant.py
index 1f678120e9daa8d4620fb16e08ad77b28564e8a7..21b6084c872957594ecb0bc741cc0ea0a8d132b3 100644
--- a/examples/MCMC_Continuous_Plant.py
+++ b/examples/MCMC_Continuous_Plant.py
@@ -1,13 +1,13 @@
-from primo.core import BayesNet
-from primo.reasoning import ContinuousNodeFactory
-from primo.reasoning.density import ExponentialParameters
-from primo.reasoning.density import BetaParameters
-from primo.reasoning.density import GaussParameters
-from primo.reasoning.density import NDGauss
-from primo.reasoning import MCMC
+from primo.networks import BayesianNetwork
+from primo.nodes import ContinuousNodeFactory
+from primo.densities import ExponentialParameters
+from primo.densities import BetaParameters
+from primo.densities import GaussParameters
+from primo.densities import NDGauss
+from primo.inference.mcmc import MCMC
 
-from primo.reasoning import EvidenceEqual as EvEqual
-from primo.reasoning import EvidenceInterval as EvInterval
+from primo.evidence import EvidenceEqual as EvEqual
+from primo.evidence import EvidenceInterval as EvInterval
 
 import numpy
 
@@ -16,7 +16,7 @@ import numpy
 #it has grown up to (+noise)
 
 
-bn = BayesNet()
+bn = BayesianNetwork()
 cnf=ContinuousNodeFactory()
 
 #create the nodes
diff --git a/examples/MCMC_Discrete.py b/examples/MCMC_Discrete.py
index 245aed52875efcd779a1bc2523cd86f332fe5c9a..539c08e542868e5f4c0f867937c10517093e679b 100644
--- a/examples/MCMC_Discrete.py
+++ b/examples/MCMC_Discrete.py
@@ -1,11 +1,11 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-from  primo.core import BayesNet
-from  primo.reasoning import DiscreteNode
-from primo.reasoning.density import ProbabilityTable
-from primo.reasoning import MCMC
-from primo.reasoning import EvidenceEqual as EvEq
-from primo.reasoning import GibbsTransitionModel
+from primo.networks import BayesianNetwork
+from primo.nodes import DiscreteNode
+from primo.densities import ProbabilityTable
+from primo.inference.mcmc import MCMC
+from primo.evidence import EvidenceEqual as EvEq
+from primo.inference.mcmc import GibbsTransitionModel
 import numpy
 import pdb
 
@@ -19,7 +19,7 @@ import pdb
 #-Maximum a-posteriori hypothesis
 
 #Construct some simple BayesianNetwork
-bn = BayesNet()
+bn = BayesianNetwork()
 burglary = DiscreteNode("Burglary", ["Intruder","Safe"])
 alarm = DiscreteNode("Alarm", ["Ringing", "Silent","Destroyed"])
 
diff --git a/examples/Robot_Localization.xmlbif.bak b/examples/Robot_Localization.xmlbif.bak
deleted file mode 100644
index 604be2f88fce097f3233884b7d8d316178cac1d5..0000000000000000000000000000000000000000
--- a/examples/Robot_Localization.xmlbif.bak
+++ /dev/null
@@ -1,79 +0,0 @@
-<?xml version="1.0" encoding="US-ASCII"?>
-
-<!--
-  Bayesian network in XMLBIF v0.3 (BayesNet Interchange Format)
-  Produced by SamIam http://reasoning.cs.ucla.edu/samiam
-  Output created May 13, 2013 12:09:30 PM
--->
-
-<BIF VERSION="0.3">
-<NETWORK>
-  <NAME>bayesiannetwork</NAME>
-
-  <VARIABLE TYPE="nature">
-    <NAME>x</NAME>
-    <OUTCOME>p0</OUTCOME>
-    <OUTCOME>p1</OUTCOME>
-    <OUTCOME>p2</OUTCOME>
-    <OUTCOME>p3</OUTCOME>
-    <OUTCOME>p4</OUTCOME>
-    <OUTCOME>p5</OUTCOME>
-    <OUTCOME>p6</OUTCOME>
-    <OUTCOME>p7</OUTCOME>
-    <OUTCOME>p8</OUTCOME>
-    <OUTCOME>p9</OUTCOME>
-    <PROPERTY>position = (477, -200)</PROPERTY>
-  </VARIABLE>
-
-  <VARIABLE TYPE="nature">
-    <NAME>door</NAME>
-    <OUTCOME>True</OUTCOME>
-    <OUTCOME>False</OUTCOME>
-    <PROPERTY>position = (479, -65)</PROPERTY>
-  </VARIABLE>
-
-  <VARIABLE TYPE="nature">
-    <NAME>x0</NAME>
-    <OUTCOME>p0</OUTCOME>
-    <OUTCOME>p1</OUTCOME>
-    <OUTCOME>p2</OUTCOME>
-    <OUTCOME>p3</OUTCOME>
-    <OUTCOME>p4</OUTCOME>
-    <OUTCOME>p5</OUTCOME>
-    <OUTCOME>p6</OUTCOME>
-    <OUTCOME>p7</OUTCOME>
-    <OUTCOME>p8</OUTCOME>
-    <OUTCOME>p9</OUTCOME>
-    <PROPERTY>position = (260, -200)</PROPERTY>
-  </VARIABLE>
-
-  <DEFINITION>
-    <FOR>x</FOR>
-    <GIVEN>x0</GIVEN>
-    <TABLE>
-    0.05 0.9 0.05 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-    0.0 0.05 0.9 0.05 0.0 0.0 0.0 0.0 0.0 0.0
-    0.0 0.0 0.05 0.9 0.05 0.0 0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 0.05 0.9 0.05 0.0 0.0 0.0 0.0
-    0.0 0.0 0.0 0.0 0.05 0.9 0.05 0.0 0.0 0.0
-    0.0 0.0 0.0 0.0 0.0 0.05 0.9 0.05 0.0 0.0
-    0.0 0.0 0.0 0.0 0.0 0.0 0.05 0.9 0.05 0.0
-    0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.05 0.9 0.05
-    0.05 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.05 0.9 
-    0.9 0.05 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.05
-    </TABLE>
-  </DEFINITION>
-
-  <DEFINITION>
-    <FOR>door</FOR>
-    <GIVEN>x</GIVEN>
-    <TABLE>0.01 0.99 0.99 0.01 0.01 0.99 0.99 0.01 0.01 0.99 0.01 0.99 0.01 0.99 0.99 0.01 0.01 0.99 0.01 0.99 </TABLE>
-  </DEFINITION>
-
-  <DEFINITION>
-    <FOR>x0</FOR>
-    <TABLE>0.10000000000000002 0.10000000000000002 0.10000000000000002 0.10000000000000002 0.10000000000000002 0.10000000000000002 0.10000000000000002 0.10000000000000002 0.10000000000000002 0.10000000000000002 </TABLE>
-  </DEFINITION>
-
-</NETWORK>
-</BIF>
diff --git a/examples/Robot_Localization.xmlbif.bak2 b/examples/Robot_Localization.xmlbif.bak2
deleted file mode 100644
index 5f5fd0c0ca105a4bd3b011b8101bc080e5352c9c..0000000000000000000000000000000000000000
--- a/examples/Robot_Localization.xmlbif.bak2
+++ /dev/null
@@ -1,79 +0,0 @@
-<?xml version="1.0" encoding="US-ASCII"?>
-
-<!--
-  Bayesian network in XMLBIF v0.3 (BayesNet Interchange Format)
-  Produced by SamIam http://reasoning.cs.ucla.edu/samiam
-  Output created May 13, 2013 12:09:30 PM
--->
-
-<BIF VERSION="0.3">
-<NETWORK>
-  <NAME>bayesiannetwork</NAME>
-
-  <VARIABLE TYPE="nature">
-    <NAME>x</NAME>
-    <OUTCOME>p0</OUTCOME>
-    <OUTCOME>p1</OUTCOME>
-    <OUTCOME>p2</OUTCOME>
-    <OUTCOME>p3</OUTCOME>
-    <OUTCOME>p4</OUTCOME>
-    <OUTCOME>p5</OUTCOME>
-    <OUTCOME>p6</OUTCOME>
-    <OUTCOME>p7</OUTCOME>
-    <OUTCOME>p8</OUTCOME>
-    <OUTCOME>p9</OUTCOME>
-    <PROPERTY>position = (477, -200)</PROPERTY>
-  </VARIABLE>
-
-  <VARIABLE TYPE="nature">
-    <NAME>door</NAME>
-    <OUTCOME>True</OUTCOME>
-    <OUTCOME>False</OUTCOME>
-    <PROPERTY>position = (479, -65)</PROPERTY>
-  </VARIABLE>
-
-  <VARIABLE TYPE="nature">
-    <NAME>x0</NAME>
-    <OUTCOME>p0</OUTCOME>
-    <OUTCOME>p1</OUTCOME>
-    <OUTCOME>p2</OUTCOME>
-    <OUTCOME>p3</OUTCOME>
-    <OUTCOME>p4</OUTCOME>
-    <OUTCOME>p5</OUTCOME>
-    <OUTCOME>p6</OUTCOME>
-    <OUTCOME>p7</OUTCOME>
-    <OUTCOME>p8</OUTCOME>
-    <OUTCOME>p9</OUTCOME>
-    <PROPERTY>position = (260, -200)</PROPERTY>
-  </VARIABLE>
-
-  <DEFINITION>
-    <FOR>x</FOR>
-    <GIVEN>x0</GIVEN>
-    <TABLE>
-    0.1 0.73 0.1 0.01 0.01 0.01 0.01 0.01 0.01 0.01
-    0.01 0.1 0.73 0.1 0.01 0.01 0.01 0.01 0.01 0.01
-    0.01 0.01 0.1 0.73 0.1 0.01 0.01 0.01 0.01 0.01
-    0.01 0.01 0.01 0.1 0.73 0.1 0.01 0.01 0.01 0.01
-    0.01 0.01 0.01 0.01 0.1 0.73 0.1 0.01 0.01 0.01
-    0.01 0.01 0.01 0.01 0.01 0.1 0.73 0.1 0.01 0.01
-    0.01 0.01 0.01 0.01 0.01 0.01 0.1 0.73 0.1 0.01
-    0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.1 0.73 0.1
-    0.1 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.1 0.73 
-    0.73 0.1 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.1
-    </TABLE>
-  </DEFINITION>
-
-  <DEFINITION>
-    <FOR>door</FOR>
-    <GIVEN>x</GIVEN>
-    <TABLE>0.1 0.9 0.9 0.1 0.1 0.9 0.9 0.1 0.1 0.9 0.1 0.9 0.1 0.9 0.9 0.1 0.1 0.9 0.1 0.9 </TABLE>
-  </DEFINITION>
-
-  <DEFINITION>
-    <FOR>x0</FOR>
-    <TABLE>0.10000000000000002 0.10000000000000002 0.10000000000000002 0.10000000000000002 0.10000000000000002 0.10000000000000002 0.10000000000000002 0.10000000000000002 0.10000000000000002 0.10000000000000002 </TABLE>
-  </DEFINITION>
-
-</NETWORK>
-</BIF>
diff --git a/examples/teilprojekt3_djohn_mbaum.py b/examples/teilprojekt3_djohn_mbaum.py
deleted file mode 100644
index 6a0be4da49458ec0b8669a3b5b0dde79b1f699da..0000000000000000000000000000000000000000
--- a/examples/teilprojekt3_djohn_mbaum.py
+++ /dev/null
@@ -1,752 +0,0 @@
-# -*- coding: utf-8 -*-
-# Copyright 2012 Manuel Baum and Denis John
-# Contact:
-# djohn@techfak.uni-bielefeld.de
-# mbaum@techfak.uni-bielefeld.de
-
-#    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 3 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, see <http://www.gnu.org/licenses/>.
-
-
-
-
-import random
-import itertools
-import copy
-import re
-import numpy
-import xml.dom.minidom as dom
-
-
-def weightedRandom(weights):
-    counter = random.random() * sum(weights)
-    for i,w in enumerate(weights):
-        counter -= w
-        if counter <=0:
-            return i
-    
-
-# Klasse Directed Graph
-# Diese Klasse stellt einen allgemeinen Graphen dar, der noch nicht unbedingt
-# ein Bayesnetz sein muss
-class DirectedGraph(dict):
-
-    def __init__(self):
-        self.nodeFactory = NodeFactory()
-
-    def insertNode(self,nodeName):
-        print "Inserting Node: "+nodeName
-        if nodeName in self:
-            raise Exception('Name ' + nodeName + ' already exists!')
-        self[nodeName] = self.nodeFactory.node()
-
-    def removeNode(self,nodeName):
-        if not nodeName in self:
-            raise Exception('Name ' + nodeName + ' does not exists!')
-
-        children=self[nodeName].children
-        parents=self[nodeName].parents
-        map(lambda x: self.removeEdge(nodeName,x),children)
-        map(lambda x: self.removeEdge(x,nodeName),parents)
-
-        del self[nodeName]
-
-    def calculatePositions(self):
-        q=[]
-        p = []
-        alreadySeen=[]
-        xStep = 150
-        yStep = 100
-        xPos = 0
-        yPos = 0
-        for k in self:
-            if not self.hasParent(k):
-                q.append(k)
-                alreadySeen.append(k)
-        while q:
-            p = q
-            q = []
-            yPos += yStep
-            xPos = xStep
-            while p:
-                n = p.pop()
-                self[n].position = (xPos,yPos)
-                xPos += xStep
-
-                for child in self[n].children:
-                    if not child in alreadySeen:
-                        q.append(child)
-                        alreadySeen.append(child)
-
-
-
-
-    def insertEdge(self,parentName,childName):
-        print "Inserting Edge: "+parentName+" -> "+childName
-        if self.existsEdge(parentName, childName):
-            raise Exception('An Edge already exists from '+parentName+' to '+childName)
-        self[parentName].children.append(childName)
-        self[childName].parents.append(parentName)
-        if not self.isAcyclic():
-            self.removeEdge(parentName,childName)
-            raise Exception('This BayesNet becomes cyclic by inserting Edge '+parentName+"->"+childName)
-
-    def existsEdge(self, parentName, childName):
-        return childName in self[parentName].children
-
-    def removeEdge(self,parentName,childName):
-        self[parentName].children.remove(childName)
-        self[childName].parents.remove(parentName)
-
-    def removeEdgeUnsafe(self,parentName,childName):
-        self[parentName].children.remove(childName)
-        self[childName].parents.remove(parentName)
-
-    def insertEdgeUnsafe(self,parentName,childName):
-        if self.existsEdge(parentName, childName):
-            raise Exception('An Edge already exists from '+parentName+' to '+childName)
-        self[parentName].children.append(childName)
-        self[childName].parents.append(parentName)
-
-    def isAcyclic(self):
-        graph = copy.deepcopy(self)
-        q=[]
-        for k in graph:
-            if not graph.hasParent(k):
-                q.append(k)
-        while q:
-            n = q.pop()
-            for child in graph[n].children[:]:
-                graph.removeEdge(n,child)
-                if not graph.hasParent(child):
-                    q.append(child)
-
-        for n in graph.values():
-            if n.children != []:
-                return False
-        return True
-
-    def hasParent(self, name):
-        node = self[name]
-        return node.parents
-
-    def __str__(self):
-        word='digraph G {\n'
-        for n in self.keys():
-            word+=str(n)+';\n'
-            for child in self[n].children:
-                word+=str(n)+'->'+str(child)+';\n'
-
-        word += '}'
-        return word
-
-    def writeToFile(self,filepath):
-        FILE = open(filepath,'w')
-        FILE.write(self.__str__())
-
-    def readFromFile(self, filepath):
-        with open(filepath,'r') as FILE:
-            string = FILE.read()
-        #whitespace entfernen
-        string = re.sub(r'\s', '', string)
-        #auf die geschweiften klammern zurechtschneiden
-        strings=string.split('{')
-        strings=strings[1].split('}')
-        #die einzelnen edges/nodes splitten
-        strings=strings[0].split(';')
-        #mindestens nach dem letzten semikolon bleibt ein '' übrig
-        while '' in strings:
-            strings.remove('')
-
-        for s in strings:
-            subs = s.split("->")
-            if not subs[0] in self.keys():
-                self.insertNode(subs[0])
-
-            if len(subs) > 1:
-                if not subs[1] in self.keys():
-                    self.insertNode(subs[1])
-                self.insertEdge(subs[0],subs[1])
-
-# Diese Klasse wird in einen DirectedGraph geladen
-class NodeFactory(object):
-    def node(self):
-        return Node()
-
-class RandomNodeFactory(object):
-    def node(self):
-        return RandomNode()
-
-class Node(object):
-    def __init__(self):
-        self.parents = []
-        self.children = []
-        self.position = (0,0)
-
-class RandomNode(Node):
-    def __init__(self):
-        Node.__init__(self)
-        #The CPT is a numpy.array. The first dimension (dim: 0) contains the
-        #values of the node itself. The following dimensions reflect the values
-        #of the parent nodes
-        self.cpt = numpy.array(0)
-        self.values = []
-
-    def normalize(self):
-        sumMat = numpy.sum(self.cpt,0)
-
-        if sumMat.shape == () and sumMat == 0:
-            self.cpt = numpy.ones(self.cpt.shape)
-            sumMat = numpy.sum(self.cpt,0)
-        elif not sumMat.shape == ():
-            addMat = numpy.zeros(sumMat.shape)
-            addMat[sumMat==0] += 1
-            self.cpt += addMat
-            sumMat = numpy.sum(self.cpt,0)
-
-
-        #finally normalize
-        self.cpt /= sumMat
-
-
-    def uniformize(self):
-        self.cpt = numpy.ones(self.cpt.shape)
-        self.normalize()
-
-    def setCpt(self, newCpt):
-        if newCpt.shape == self.cpt.shape:
-            self.cpt = newCpt
-            self.normalize()
-        else:
-            raise Exception('Can not set new values for this Cpt, as the new Cpt has wrong dimensionality.')
-
-    def appendCptZeroSubmatrix(self, dimension):
-        shape = self.cpt.shape
-        if shape == ():
-            newCpt = numpy.zeros((1))
-        else:
-            shapeOfAppendix=(shape[:dimension])+(1,)+(shape[dimension+1:])
-            newCpt = numpy.concatenate((self.cpt,numpy.zeros(shapeOfAppendix)),dimension)
-        self.cpt = newCpt
-
-    def deleteCptLine(self, dim, idx):
-        if idx==0:
-            B,C=numpy.split(self.cpt,[idx+1],axis=dim)
-            self.cpt=C
-        elif idx==self.cpt.shape[dim]-1:
-            A,B=numpy.split(self.cpt,[idx],axis=dim)
-            self.cpt=A
-        else:
-            A,B,C = numpy.split(self.cpt,[idx,idx+1],axis=dim)
-            self.cpt = numpy.concatenate((A,C),axis = dim)
-
-
-    def getXMLTag(self,nodeName):
-        tag_var = dom.Element("VARIABLE")
-        tag_own = dom.Element("NAME")
-        tag_pos = dom.Element("PROPERTY")
-
-        tag_var.setAttribute("TYPE","nature")
-
-        txt_name = dom.Text()
-        txt_name.data = str(nodeName)
-        tag_var.appendChild(tag_own)
-        tag_own.appendChild(txt_name)
-
-        for val in self.values:
-            tag_val = dom.Element("OUTCOME")
-            txt_val = dom.Text()
-            txt_val.data = str(val)
-            tag_val.appendChild(txt_val)
-            tag_var.appendChild(tag_val)
-
-        txt_pos = dom.Text()
-        x,y = self.position
-        txt_pos.data = "position = (" + str(x) + ", " + str(y) + ")"
-        tag_pos.appendChild(txt_pos)
-        tag_var.appendChild(tag_pos)
-
-        return tag_var
-
-
-    ### Functions for Inference
-
-    def setEvidence(self,evidenceName):        
-        idx = self.values.index(evidenceName)
-
-        newCpt = numpy.zeros(self.cpt.shape)
-
-        tmpCpt = self.cpt[idx]
-        newCpt[idx] = tmpCpt
-        self.cpt = newCpt
-
-    def rollBackParent(self,parentName):
-        if len(self.parents) == self.cpt.ndim:
-            idx = self.parents.index(parentName)
-        else:
-            idx = self.parents.index(parentName) +1
-        self.cpt = numpy.rollaxis(self.cpt,idx,self.cpt.ndim)
-        self.parents.remove(parentName)
-        self.parents.append(parentName)
-
-    def rollFrontParent(self,parentName):
-        if len(self.parents) == self.cpt.ndim:
-            idx = self.parents.index(parentName)
-        else:
-            idx = self.parents.index(parentName) +1
-        self.cpt = numpy.rollaxis(self.cpt,idx)
-        self.parents.remove(parentName)
-        self.parents.insert(0,parentName)
-
-class BayesNet(DirectedGraph):
-    def __init__(self):
-        self.nodeFactory=RandomNodeFactory()
-
-    def insertValue(self, nodeName, value):
-        print "Inserting Value: "+value+" in "+nodeName
-        node = self[nodeName]
-        if not value in node.values:
-            node.values.append(value)
-            node.appendCptZeroSubmatrix(0)
-            #node.normalize()
-
-            #node.normalize()
-            for child in node.children:
-                dim,idx = self.cptDimensionIndex(child,nodeName,value)
-                self[child].appendCptZeroSubmatrix(dim)
-                #self[child].normalize()
-
-        else:
-            raise Exception('Cannot insert the value '+value+" in the RandomNode "+name+" because a value with this name already exists in this Node")
-
-
-
-    def removeValue(self, nodeName, value):
-        node = self[nodeName]
-        if not value in node.values:
-            raise Exception('Cannot delete the value '+value+" in the RandomNode "+nodeName+" because it does not exists in this Node")
-        elif len(node.values) <= 1:
-            raise Exception('Cannot delete the value '+value+" in the RandomNode "+nodeName+" because it is the last value for this Node.")
-        else:
-            #Handle the CPTs of the Children
-            for child in node.children:
-                matDim, idx = self.cptDimensionIndex(child, nodeName, value)
-                self[child].deleteCptLine(matDim, idx)
-            #Handle the CPT of this Node
-            idx = node.values.index(value)
-            node.values.remove(value)
-            node.deleteCptLine(0,idx)
-            #node.normalize()
-
-    def insertEdge(self, fromName, toName):
-        super(BayesNet, self).insertEdge(fromName,toName)
-        node = self[toName]
-        ax = node.cpt.ndim
-        node.cpt=numpy.expand_dims(node.cpt,ax)
-        node.cpt=numpy.repeat(node.cpt,len(self[fromName].values),axis = ax)
-        #node.normalize()
-
-
-    def insertNode(self, name, values=["TRUE"]):
-        super(BayesNet, self).insertNode(name)
-        for value in values:
-            self.insertValue(name, value)
-
-    def cptDimension(self, nodeName):
-        node = self[nodeName]
-        dim=[len(node.values)]
-        for parent in node.parents:
-            dim.append(len(self[parent].values))
-        return dim
-
-#This Function returns the MatrixDimension and the Index in this
-#dimension at the which the specified valueName is located,
-#this is needed at least for removeValue
-    def cptDimensionIndex(self, nodeName, parentName, valueName):
-        matDim = 1 + self[nodeName].parents.index(parentName)
-        idx = self[parentName].values.index(valueName)
-        return (matDim, idx)
-
-    def generate_XML(self,netName):
-        print "Generate XML for: " + netName
-        self.calculatePositions()
-        mainNode = dom.Document()
-        tag_bif = mainNode.createElement("BIF")
-        tag_net = mainNode.createElement("NETWORK")
-        tag_bif.setAttribute("VERSION","0.3")
-        mainNode.appendChild(tag_bif)
-        tag_bif.appendChild(tag_net)
-
-        tag_name = dom.Element("NAME")
-        text = dom.Text()
-        text.data = str(netName)
-        tag_name.appendChild(text)
-        tag_net.appendChild(tag_name)
-
-        for nodeName in self.keys():
-            curNode = self[nodeName]
-            node_tag = curNode.getXMLTag(nodeName)
-            tag_net.appendChild(node_tag)
-
-        #Generate CPTs
-        for nodeName in self.keys():
-            curNode = self[nodeName]
-            tag_def = dom.Element("DEFINITION")
-            tag_for = dom.Element("FOR")
-            txt_for = dom.Text()
-            txt_for.data = str(nodeName)
-            tag_for.appendChild(txt_for)
-            tag_def.appendChild(tag_for)
-
-            for parent in reversed(curNode.parents):
-                tag_par = dom.Element("GIVEN")
-                txt_par = dom.Text()
-                txt_par.data = str(parent)
-                tag_par.appendChild(txt_par)
-                tag_def.appendChild(tag_par)
-
-            tag_cpt = dom.Element("TABLE")
-            txt_cpt = dom.Text()
-            txt = ""
-            for elem in curNode.cpt.T.flat:
-                txt += str(elem) + " "
-
-            txt_cpt.data = txt
-            tag_cpt.appendChild(txt_cpt)
-            tag_def.appendChild(tag_cpt)
-
-            tag_net.appendChild(tag_def)
-
-
-        return mainNode
-
-    def normalize(self):
-        for nodeName in self:
-            self[nodeName].normalize()
-
-    def updateCpts(self):
-        for nodeName in self:
-            self.updateCpt(nodeName)
-
-
-#### Inference Algorithms
-    #This function computes the prior marginal-distribution over a set of query-variables
-    def priorMarginal(self,queryVariables = []):
-        bn = copy.deepcopy(self)
-        #Assume all variables are of interest in case of [] as queryVariables
-        if(queryVariables == []):
-            queryVariables = bn.keys()
-
-        #Look at all Nodes in the BN
-        q = bn.keys()
-        while q:
-            curNodeName = q.pop()
-            curNode = bn[curNodeName]
-            #If this Node has no parents, it's children become potentially new
-            #unconditioned Variables when this node has been deleted
-            if len(curNode.parents) == 0:
-                q.extend(curNode.children)
-                bn.conditionElimination(curNodeName)
-
-        #Construct the Output as list of (nodename, node) pairs
-        returnList = []
-        for name in bn.keys():
-            if name in queryVariables:
-                returnList.append((name,bn[name].cpt))
-
-        return returnList
-                
-        
-    #This function removes edges to children of this node and eliminates the conditioning of the children-nodes on the specified one
-    def conditionElimination(self,nodeName):
-        node = self[nodeName]
-        #It is necessary, that this node has no parents, otherwise something has gone wrong
-        #Before this function has been called
-        if node.parents:
-            raise Exception("There shouldn't be a parent for this node at this time")
-
-        #In the CPT of each child, the values of this node are being marginalized out
-        for childName in node.children[:]:
-            childNode = self[childName]
-            childNode.rollBackParent(nodeName)
-            #Practically dot-product combines pointwise multiplication and summation over
-            #possible values of this node
-            childNode.cpt = numpy.dot(childNode.cpt,node.cpt)
-            self.removeEdgeUnsafe(nodeName,childName)
-        
-      
-
-    def probabilityOfEvidence(self,evidences = []):
-        bn = copy.deepcopy(self)
-
-        #First of all, we have to set the specified evidence in the bn
-        for evName,value in evidences:
-            if not evName in bn:
-                raise Exception("The node " + evName + " doesn't exist for calculating the probability of evidence")
-            elif not value in bn[evName].values:
-                raise Exception("The value " + value + " doesn't exist in node " + evName + " for calculating the probability of evidence")
-            bn[evName].setEvidence(value)
-
-        #Insert all nodes without parents in q
-        q = []
-        for k in bn:
-            if not bn.hasParent(k):
-                q.append(k)
-
-        #from these nodes go down recursively to the leaves of the DAG and eliminate upwards from there
-        for p in q:
-            bn.eliminateChildren(p)
-
-        #Finally multiply the partial results for divided graphs
-        result = 1
-        for k in bn:
-            result *= bn[k].cpt
-
-        return result
-
-    #This function recursively eliminates the children of this node
-    #and sums the different possible values of this node out too
-    def eliminateChildren(self,nodeName):
-        node = self[nodeName]
-        #Before we can eliminate this node, we will eliminate all of it's children
-        #recursively
-        for child in node.children:
-            self.eliminateChildren(child)
-
-        #Now that the values of the children have been summed out, we can
-        #combine their cpts and sum ourself out
-        self.uniteNodes(node.children[:],nodeName)
-        self.sumOutVariable(nodeName)
-
-    def sumOutVariable(self,nodeName):
-        node = self[nodeName]
-        #Check if summation is correct as we make expectations based on our bottom-up elimination of the DAG
-        if len(node.children) > 1:
-            raise Exception("sumOutVariable: The node has more than one child")
-        if len(node.parents) +1 < node.cpt.ndim:
-            raise Exception("sumOutVariable: The cpt of the node has the wrong number of dims")
-
-
-        if len(node.parents) == node.cpt.ndim:
-            #Already summed out, this happens if this node has multiple parents
-            return
-
-        #Trivial case: node has no children, just sum it's vector-shaped cpt
-        if len(node.children) == 0:
-            node.cpt = sum(node.cpt,0)
-        #Not so trivial case: node has children (necessarily only one as the childrens CPTs have already been combined into a single matrix)
-        else:
-            childName = node.children[0]
-            childNode = self[childName]
-
-            #combine the beliefs of the children nodes (already combined in one matrix) and our own cpt. This is done by dotwise multiplication and then summation over our own values.
-            node.cpt = numpy.rollaxis(node.cpt,0,node.cpt.ndim)
-            childNode.rollFrontParent(nodeName)           
-            node.cpt = numpy.dot(node.cpt,childNode.cpt)
-
-            #Finally remove this nodes child and connect this node to the other parents of that child
-            self.removeEdgeUnsafe(nodeName,childName)
-            for p in childNode.parents:
-                self.insertEdgeUnsafe(p,nodeName)
-            self.removeNode(childName)
-
-
-
-    def uniteNodes(self,nodeList,parentName):
-        if not nodeList:
-            return
-        #Check if it is sensible to unite the given nodes
-        for nodeName in nodeList:
-            node = self[nodeName]
-            if not node.cpt.ndim == len(node.parents):
-                raise Exception("Unite Nodes: CptDim and Size parents aren't the same")
-            if not parentName in node.parents:
-                raise Exception("Unite Nodes: At least one node doesn't have the given parent")
-
-        #Iteratively unite all nodes with the first one, so all information is contained in one matrix
-        firstNode = nodeList.pop()
-        while nodeList:
-            curNode = nodeList.pop()
-            firstNode = self.uniteTwoNodes(firstNode,curNode,parentName)
-
-    
-    
-    def uniteTwoNodes(self,nodeName1,nodeName2,parentName):
-        node1 = self[nodeName1]
-        node2 = self[nodeName2]
-        idx1 = node1.parents.index(parentName)
-        idx2 = node2.parents.index(parentName)
-
-        #Roll parent to the front
-        node1.parents.remove(parentName)
-        node1.parents.insert(0,parentName)
-        node2.parents.remove(parentName)
-        node2.parents.insert(0,parentName)
-        #Roll parent to the front in the cpt
-        cpt1 = numpy.rollaxis(node1.cpt,idx1)
-        cpt2 = numpy.rollaxis(node2.cpt,idx2)
-
-        #The CPTS are being combined
-        node1.cpt = self.uniteCpts(cpt1,cpt2)
-
-        #Remove the second Node and redirect the edges from it's parents to the first node
-        self.removeEdgeUnsafe(parentName,nodeName2)
-        for p in node2.parents:
-            self.insertEdgeUnsafe(p,nodeName1)
-        self.removeNode(nodeName2)
-
-        
-
-    def uniteCpts(self,cpt1,cpt2):
-        #Assuming the common dimension at zero
-        
-        #New shape for new cpt
-        newShape = cpt2.shape[1:] + cpt1.shape
-        newCpt = numpy.ones(newShape)
-        ndim1 = cpt1.ndim
-        ndim2 = cpt2.ndim
-
-        #Insert Values of the first cpt
-        newCpt *= cpt1
-
-        #Rotate Axis so insertion of the values for cpt2 becomes a pointwise matrixdimension
-        newCpt = numpy.rollaxis(newCpt,ndim2-1,newCpt.ndim)
-        for i in range(ndim2-1):
-            newCpt = numpy.rollaxis(newCpt,0,newCpt.ndim)
-        newCpt *= cpt2
-        
-        newCpt = numpy.rollaxis(newCpt,ndim1-1)
-        
-        return newCpt
-
-    def mcmc(self, evidences, queries, times):
-        bn = copy.deepcopy(self)
-
-        #Evidence and NonEvidence Variables
-        eVariables = []
-        neVariables = []
-
-        #First of all, we have to set the specified evidence in the bn
-        for evName,value in evidences:
-            if not evName in bn:
-                raise Exception("The node " + evName + " doesn't exist in this BayesNet")
-            elif not value in bn[evName].values:
-                raise Exception("The value " + value + " doesn't exist in node " + evName + " in this BayesNet")
-            elif evName in queries:
-                raise Exception("The node " + evName + " is being queried although evidence has been provided for it.")
-            eVariables += [evName]
-            bn[evName].state=bn[evName].values.index(value)
-
-        #Get the non-evidence-variables
-        neVariables = [x for x in self.keys() if x not in eVariables]
-        
-        #Insert the histograms into the non-evidence-nodes
-        for nodeName in neVariables:
-            node = bn[nodeName] 
-            node.histogram = numpy.zeros(len(node.values))
-            node.state = random.randint(0,len(node.values)-1)
-
-        #Sample times times
-        for t in range(0,times):
-            for nodeName in neVariables:
-                node=bn[nodeName]
-                node.histogram[node.state]+=1
-
-            for nodeName in neVariables:
-                node = bn[nodeName]
-
-                #Extract the vector of probabilities according to the parent-states
-                prob = node.cpt                
-                for parentDim,parentName in enumerate(node.parents):
-                    parentState = bn[parentName].state
-                    prob = prob.take([parentState],axis=parentDim+1)
-                prob = prob.flatten()
-                for childName in node.children:
-                    child = bn[childName]
-                    childProb = child.cpt                
-                    for parentDim,parentName in [(dim, name) for (dim,name) in enumerate(child.parents) if not name == nodeName]:
-                        parentState = bn[parentName].state
-                        childProb = childProb.take([parentState],axis=parentDim+1)
-                    childProb = childProb.take([child.state],axis=0)
-                    prob *= childProb.flatten() / numpy.sum(childProb.flatten())
-                prob /= numpy.sum(prob)
-                node.state = weightedRandom(prob)
-
-        ret = []
-        for nodeName,node in [(name, bn[name]) for name in queries]:
-            node.histogram /= numpy.sum(node.histogram)
-            ret += [(nodeName,node.histogram)]
-        return ret
-
-
-    def markovBlanket(self, nodeName):
-        node = self[nodeName]
-        blanket = node.parents + node.children
-        for child in node.children:
-            blanket += [p for p in child.parents if p not in blanket]
-        return blanket
-        
-bne = BayesNet()
-
-# Nodes and Values
-bne.insertNode("Earthquake")
-bne.insertValue("Earthquake","FALSE")
-
-bne.insertNode("Burglary")
-bne.insertValue("Burglary","FALSE")
-
-bne.insertNode("Alarm")
-bne.insertValue("Alarm","FALSE")
-
-bne.insertNode("JohnCalls")
-bne.insertValue("JohnCalls","FALSE")
-
-bne.insertNode("BaumCalls")
-bne.insertValue("BaumCalls","FALSE")
-
-#Edges
-bne.insertEdge("Burglary","Alarm")
-bne.insertEdge("Earthquake","Alarm")
-bne.insertEdge("Alarm","JohnCalls")
-bne.insertEdge("Alarm","BaumCalls")
-
-#CPTs
-bne["Burglary"].setCpt(numpy.array(\
-    [0.001 , 0.999]))
-
-bne["Earthquake"].setCpt(numpy.array(\
-    [0.002 , 0.998]))
-
-bne["Alarm"].setCpt(numpy.array(\
-    [[[ 0.95,0.94],\
-      [ 0.29,0.001]],\
-     [[ 0.05,0.06],\
-      [ 0.71,0.999]]]))
-
-bne["BaumCalls"].setCpt(numpy.array(\
-    [[0.9,0.05],\
-     [0.1,0.95]]))
-bne["JohnCalls"].setCpt(numpy.array(\
-    [[0.7,0.01],\
-     [0.3,0.99]]))
-
-###Inference Algorithms
-print "Probability of Evidence"
-#print bne.probRec([("JohnCalls","TRUE"),("BaumCalls","TRUE")])
-print bne.probabilityOfEvidence([("Earthquake","FALSE")])
-
-print "Prior Marginal"
-print  bne.priorMarginal()
-
-print "Posterior Marginal"
-print bne.mcmc([("Alarm","TRUE"),("Earthquake","TRUE")],["JohnCalls"],10000)
diff --git a/primo/densities.py b/primo/densities.py
index da0b009a87bc1107f17384e6b1357d70ae3b09bb..5d4d229d1ec51032f5b9e7f2d8cb5a4775bf5ea6 100644
--- a/primo/densities.py
+++ b/primo/densities.py
@@ -1,13 +1,13 @@
 # -*- coding: utf-8 -*-
-
 import copy
 import math
 import random
 
 import numpy
-from scipy.stats
+import scipy.stats
 
-from primo.nodes
+import primo.nodes
+from primo.util import weighted_random
 
 class Density(object):
     '''TODO: write doc'''
@@ -116,7 +116,7 @@ class Beta(Density):
         self.q=parameters.q
 
     def add_variable(self, variable):
-        if( not isinstance(variable,ContinuousNode.ContinuousNode)):
+        if( not isinstance(variable,primo.nodes.ContinuousNode)):
             raise Exception("Tried to add Variable as parent, but is not a ContinuousNode")
         self.p[variable]=0.0
         self.q[variable]=0.0
@@ -223,7 +223,7 @@ class Exponential(Density):
         to be parents because of their value range. Instead it should be evaluated
         if they can yield parameters for this distribution that are permitted. This
         can in any case happen under bad influence coefficients'''
-        if( not isinstance(variable,ContinuousNode.ContinuousNode)):
+        if( not isinstance(variable,primo.nodes.ContinuousNode)):
             raise Exception("Tried to add Variable as parent, but is not a ContinuousNode")
         self.b[variable]=0.0
 
diff --git a/primo/inference/factor.py b/primo/inference/factor.py
index fe31e053b3a8570c1f9d425103cacc97aab301b4..ae3c618e4b667b043bcd4b518d31cac485ea313b 100644
--- a/primo/inference/factor.py
+++ b/primo/inference/factor.py
@@ -1,3 +1,6 @@
+import copy
+import operator
+
 import networkx as nx
 
 import primo.densities
@@ -177,7 +180,7 @@ class FactorTreeFactory(object):
         for n in allNodes:
             sortNodeList.append((len(n.get_cpd().get_variables()),n))
         #sort node list
-        sortNodeList = sorted(sortNodeList,key=itemgetter(0),reverse=True)
+        sortNodeList = sorted(sortNodeList,key=operator.itemgetter(0),reverse=True)
         sortNodeList =  zip(*sortNodeList)
         sortNodeList = list(sortNodeList[1])
         #root node with the most variables
diff --git a/primo/inference/mcmc.py b/primo/inference/mcmc.py
index e8c25c13c25ba69820f26cbaf9612be9ff3057e4..e89d14609afdcce0ee498e2fd9b7c3e28eba7b65 100644
--- a/primo/inference/mcmc.py
+++ b/primo/inference/mcmc.py
@@ -1,9 +1,7 @@
-
 import copy
 import random
 
-
-from primo.reasoning.convergence_test import ConvergenceTestSimpleCounting
+from primo.util import weighted_random
 
 
 class MCMC(object):
@@ -123,19 +121,6 @@ class MCMC(object):
             state[var]=var.sample_global(state, evidence)
         return state
 
-
-def weighted_random(weights):
-    '''
-    Implements roulette-wheel-sampling.
-    @param weights: A List of float-values.
-    @returns: Index of the selected entity
-    '''
-    counter = random.random() * sum(weights)
-    for i,w in enumerate(weights):
-        counter -= w
-        if counter <=0:
-            return i
-
 class GibbsTransitionModel(object):
     '''
     Implements Gibbs-sampling. Can be used to constuct a Markov Chain and is
diff --git a/primo/inference/particlefilter.py b/primo/inference/particlefilter.py
index 72f5b1d7916d8da52b185dc3228db324b00124d4..2b8dc245e2a76bab73f5724f828e63f886fe76cd 100644
--- a/primo/inference/particlefilter.py
+++ b/primo/inference/particlefilter.py
@@ -5,6 +5,7 @@ import random
 import time
 
 import primo.networks
+from primo.util import weighted_random
 
 class Particle(object):
     '''
@@ -38,13 +39,6 @@ class Particle(object):
         '''
         pass
 
-def weighted_random(weights):
-    counter = random.random() * sum(weights)
-    for i, w in enumerate(weights):
-        counter -= w
-        if counter <= 0:
-            return i
-
 def wighted_sample_with_replacement(samples = [], weights = [], N = 0):
     '''
     The population is resampled to generate a new population of N samples.
@@ -90,7 +84,7 @@ def weighted_sample(network, evidence = {}):
     '''
     w = 1.0
     state = {}
-    if not isinstance(network, primo.network.BayesianNetwork):
+    if not isinstance(network, primo.networks.BayesianNetwork):
         raise Exception("The given network is not an instance of BayesNet.")
 
     nodes = network.get_nodes_in_topological_sort()
diff --git a/primo/io.py b/primo/io.py
index 8c97a77232eac5edc07f0efdf515e8b1095e3880..78e0854de2dc9b02a32da005ecf86d961f606675 100644
--- a/primo/io.py
+++ b/primo/io.py
@@ -1,9 +1,9 @@
 # -*- coding: utf-8 -*-
-import xml.dom.minidom as minidom
-from primo.core import BayesNet
-from primo.core import Node
-from primo.reasoning import DiscreteNode
 import re
+import xml.dom.minidom as minidom
+
+import primo.networks
+import primo.nodes
 
 
 class XMLBIF(object):
@@ -35,7 +35,7 @@ class XMLBIF(object):
         self.newl = newl
         self.addindent = addindent
         self.root = minidom.Document()
-        if isinstance(network, BayesNet):
+        if isinstance(network, primo.networks.BayesianNetwork):
             self.network = network
         else:
             raise Exception("Given network is not a BayesNet.")
@@ -82,7 +82,7 @@ class XMLBIF(object):
 
         for node_name in self.network.node_lookup:
             current_node = self.network.node_lookup[node_name]
-            if not isinstance(current_node, DiscreteNode):
+            if not isinstance(current_node, primo.nodes.DiscreteNode):
                 raise Exception("Node " + str(current_node) + " is not a DiscreteNode.")
             node_tag = self.create_node_tag(current_node)
             tag_net.appendChild(node_tag)
@@ -139,7 +139,7 @@ class XMLBIF(object):
 
         Returns a XMLBIF conform "variable" tag
         '''
-        if not isinstance(node, Node):
+        if not isinstance(node, primo.nodes.Node):
             raise Exception("Node " + str(node) + " is not a Node.")
         tag_var = minidom.Element("VARIABLE")
         tag_own = minidom.Element("NAME")
@@ -228,7 +228,7 @@ class XMLBIF(object):
 
         This method is used internally. Do not call it outside this class.
         '''
-        network = BayesNet()
+        network = primo.networks.BayesianNetwork()
         bif_nodes = root.getElementsByTagName("BIF")
         if len(bif_nodes) != 1:
             raise Exception("More than one or none <BIF>-tag in document.")
@@ -248,7 +248,7 @@ class XMLBIF(object):
             for position_node in variable_node.getElementsByTagName("PROPERTY"):
                 position = XMLBIF.get_node_position_from_text(position_node.childNodes)
                 break
-            new_node = DiscreteNode(name, value_range)
+            new_node = primo.nodes.DiscreteNode(name, value_range)
             new_node.position = position
             network.add_node(new_node)
         definition_nodes = network_nodes[0].getElementsByTagName("DEFINITION")
diff --git a/primo/networks.py b/primo/networks.py
index a86706db6eb99b29b3af2fe4a089253a55de7e8b..38c10b76f006db592aa58bf730d068230c1bc46e 100644
--- a/primo/networks.py
+++ b/primo/networks.py
@@ -1,5 +1,6 @@
 import networkx as nx
 
+import primo.densities
 import primo.nodes
 
 class BayesianNetwork(object):
@@ -152,20 +153,20 @@ class BayesianDecisionNetwork(BayesianNetwork):
                 return False
         decisionNodeList = []
         for node in self.get_all_nodes():
-            if isinstance(node, DecisionNode):
+            if isinstance(node, primo.nodes.DecisionNode):
                 decisionNodeList.append(node)
 
         return all([nx.has_path(self.graph, x, y) == True for x in decisionNodeList for y in decisionNodeList])
 
     def add_node(self, node):
-        if isinstance(node, Node):
+        if isinstance(node, primo.nodes.Node):
             if node.name in self.node_lookup.keys():
                 raise Exception("Node name already exists in Bayesnet: "+node.name)
-            if isinstance(node, DiscreteNode):
+            if isinstance(node, primo.nodes.DiscreteNode):
                 self.random_nodes.append(node)
-            elif isinstance(node, UtilityNode):
+            elif isinstance(node, primo.nodes.UtilityNode):
                 self.utility_nodes.append(node)
-            elif isinstance(node, DecisionNode):
+            elif isinstance(node, primo.nodes.DecisionNode):
                 self.decision_nodes.append(node)
             else:
                 raise Exception("Tried to add a node which the Bayesian Decision Network can not work with")
@@ -193,9 +194,9 @@ class BayesianDecisionNetwork(BayesianNetwork):
         node_from -- Node from where the edge shall begin
         node_to -- Node where the edge shall end
         """
-        if isinstance(node_from, DecisionNode) and isinstance(node_to, DecisionNode):
+        if isinstance(node_from, primo.nodes.DecisionNode) and isinstance(node_to, primo.nodes.DecisionNode):
             raise Exception("Tried to add an edge from a DecisionNode to a DecisionNode")
-        if isinstance(node_from, UtilityNode) and isinstance(node_to, UtilityNode):
+        if isinstance(node_from, primo.nodes.UtilityNode) and isinstance(node_to, primo.nodes.UtilityNode):
             raise Exception("Tried to add an edge from a UtilityNode to a UtilityNode")
         if node_from in self.graph.nodes() and node_to in self.graph.nodes():
             self.graph.add_edge(node_from, node_to)
@@ -315,7 +316,7 @@ class TwoTBN(BayesianNetwork):
         new value.
         '''
         for (node, node_t) in self.__initial_nodes:
-            cpd = ProbabilityTable()
+            cpd = primo.densities.ProbabilityTable()
             cpd.add_variable(node)
             node.set_cpd(cpd)
             if not initial:
diff --git a/primo/nodes.py b/primo/nodes.py
index 9abc383514829ccc5edc95fd0054693a9b34d56c..9697208eed1879282a60e5c403d108dd4a44a569 100644
--- a/primo/nodes.py
+++ b/primo/nodes.py
@@ -5,6 +5,7 @@ import re
 import scipy.stats
 
 import primo.densities
+import primo.inference.decision
 
 
 class Node(object):
@@ -376,7 +377,7 @@ class UtilityNode(Node):
         name -- The name of this node
         """
         super(UtilityNode, self).__init__(name)
-        self.ut = UtilityTable()
+        self.ut = primo.inference.decision.UtilityTable()
 
     def announce_parent(self, node):
         """
diff --git a/primo/tests/BayesNet_test.py b/primo/tests/BayesNet_test.py
index 24975a17dda00d79b7ed63e1b04b7da5488dac5b..ca85815402fae1f2be22fd5f4adaa85355c40140 100644
--- a/primo/tests/BayesNet_test.py
+++ b/primo/tests/BayesNet_test.py
@@ -1,12 +1,12 @@
 import unittest
 
-from primo.core import BayesNet
-from primo.reasoning import DiscreteNode
+from primo.networks import BayesianNetwork
+from primo.nodes import DiscreteNode
 
 
 class NodeAddAndRemoveTestCase(unittest.TestCase):
     def setUp(self):
-        self.bn = BayesNet()
+        self.bn = BayesianNetwork()
 
     def tearDown(self):
         self.bn = None
diff --git a/primo/tests/DynamicBayesNet_test.py b/primo/tests/DynamicBayesNet_test.py
index c8b4bfc9976118320171ddada8eab710d6d77ee8..ad384b20241607650ea85ffc4cdc3b91d6e78b64 100644
--- a/primo/tests/DynamicBayesNet_test.py
+++ b/primo/tests/DynamicBayesNet_test.py
@@ -1,23 +1,23 @@
 # -*- coding: utf-8 -*-
 import unittest
 
-from primo.core import DynamicBayesNet
-from primo.reasoning import DiscreteNode
+from primo.networks import DynamicBayesianNetwork
+from primo.nodes import DiscreteNode
 
 
 class DynamicBayesNetTest(unittest.TestCase):
     def setUp(self):
-        self.dbn = DynamicBayesNet()
+        self.dbn = DynamicBayesianNetwork()
 
     def tearDown(self):
         self.dbn = None
 
     def test_add_node(self):
         self.dbn.clear()
-        n = DiscreteNode("Some Node", [True, False])
+        n = DiscreteNode("Some_Node", [True, False])
         self.dbn.add_node(n)
-        self.assertEqual(n, self.dbn.get_node("Some Node"))
-        self.assertTrue(n in self.dbn.get_nodes(["Some Node"]))
+        self.assertEqual(n, self.dbn.get_node("Some_Node"))
+        self.assertTrue(n in self.dbn.get_nodes(["Some_Node"]))
 
     def test_temporal_edges(self):
         self.dbn.clear()
@@ -28,7 +28,8 @@ class DynamicBayesNetTest(unittest.TestCase):
         self.assertTrue(self.dbn.is_valid())
         self.dbn.add_edge(n1, n1)
         self.assertFalse(self.dbn.is_valid())
-        self.dbn.add_edge(n1, n1, True)
+        self.dbn.remove_edge(n1, n1)
+        self.dbn.add_edge(n1, n2)
         self.assertTrue(self.dbn.is_valid())
 
 
diff --git a/primo/tests/ProbabilityTable_test.py b/primo/tests/ProbabilityTable_test.py
index f187009f0f657aad3e21fce1e9d401c531370c4e..dff7b382799f077118d664bc44e647914de66e91 100644
--- a/primo/tests/ProbabilityTable_test.py
+++ b/primo/tests/ProbabilityTable_test.py
@@ -1,8 +1,8 @@
 import unittest
 import numpy
 
-from primo.reasoning.density import ProbabilityTable
-from primo.reasoning import DiscreteNode
+from primo.densities import ProbabilityTable
+from primo.nodes import DiscreteNode
 
 
 class MultiplicationTest(unittest.TestCase):
diff --git a/primo/tests/XMLBIF_test.py b/primo/tests/XMLBIF_test.py
index 9c10942d6de41b37ac5a2dcf71b4efb0c48bb3b7..81e739a5474544651918ec9ebd2e763f5ddfbe8c 100644
--- a/primo/tests/XMLBIF_test.py
+++ b/primo/tests/XMLBIF_test.py
@@ -1,15 +1,18 @@
+
+import os
 import unittest
 
-from primo.utils import XMLBIF
-from primo.core import BayesNet
-from primo.reasoning import DiscreteNode
 import numpy
-import os
+
+from primo.io import XMLBIF
+from primo.networks import BayesianNetwork
+from primo.nodes import DiscreteNode
+
 
 class ImportExportTest(unittest.TestCase):
     def setUp(self):
         # Create BayesNet
-        self.bn = BayesNet();
+        self.bn = BayesianNetwork();
         # Create Nodes
         weather0 = DiscreteNode("Weather0", ["Sun", "Rain"])
         weather = DiscreteNode("Weather", ["Sun", "Rain"])
@@ -39,17 +42,19 @@ class ImportExportTest(unittest.TestCase):
         # read BN
         bn2 = XMLBIF.read("test_out.xmlbif")
         for node1 in self.bn.get_nodes():
-            
             name_found = False
             cpd_equal = False
             value_range_equal = False
             str_equal = False
             pos_equal = False
             for node2 in bn2.get_nodes():
+                print(node2.name)
+                print(node2.get_cpd())
                 # Test node names
+                print(node2.name)
                 if node1.name == node2.name:
                     name_found = True
-                    cpd_equal = node1.get_cpd == node2.get_cpd
+                    cpd_equal = node1.get_cpd() == node2.get_cpd()
                     value_range_equal = node1.get_value_range() == node2.get_value_range()
                     str_equal = str(node1) == str(node2)
                     pos_equal = node1.pos == node2.pos
diff --git a/primo/util.py b/primo/util.py
new file mode 100644
index 0000000000000000000000000000000000000000..ce82f54342abc78e6efb8a7267103c84f6256c2c
--- /dev/null
+++ b/primo/util.py
@@ -0,0 +1,13 @@
+import random
+
+def weighted_random(weights):
+    '''
+    Implements roulette-wheel-sampling.
+    @param weights: A List of float-values.
+    @returns: Index of the selected entity
+    '''
+    counter = random.random() * sum(weights)
+    for i, w in enumerate(weights):
+        counter -= w
+        if counter <= 0:
+            return i
\ No newline at end of file