new functions added

NonLinearSolver/modelReduction/Tree.cpp

@@ -33,7 +33,7 @@ TreeNode::~TreeNode()
 
 void TreeNode::printData() const
 {
-  printf("Node: depth %d- loc %d  -mat index %d- child order %d - weight = %e - nb of direction vars %d",depth,location,matIndex,childOrder,weight,direction.size());
+  printf("Node: depth %d- loc %d  -mat index %d- child order %d - weight = %e - nb of direction vars %ld",depth,location,matIndex,childOrder,weight,direction.size());
   for (int i=0; i< direction.size(); i++)
   {
     printf(" %e",direction[i]);
@@ -49,7 +49,7 @@ void TreeNode::saveToFile(FILE* f) const
   {
     parentloc = parent->location;
   }
-  fprintf(f,"%d %d %d %d %d %d %d",depth,location, matIndex, childOrder,parentloc,childs.size(), direction.size());
+  fprintf(f,"%d %d %d %d %d %ld %ld",depth,location, matIndex, childOrder,parentloc,childs.size(), direction.size());
   //
   if (childs.size() > 0)
   {
@@ -163,7 +163,6 @@ void Tree::createSpecialBinaryTree(int numNodes, int numDirVars, const std::stri
 {
   clear();
   _root = new TreeNode(0,0,NULL,0);
-  _nodeMap[std::pair<int,int>(0,0)]=_root;
   
   TreeNode* np = _root;
   int depth = 1;
@@ -176,7 +175,6 @@ void Tree::createSpecialBinaryTree(int numNodes, int numDirVars, const std::stri
     np->childs.resize(2);
     np->direction.resize(numDirVars);
     np->childs[0] = new TreeNode(depth,0,np,0,affunc);
-    _nodeMap[std::pair<int,int>(depth,0)]=np->childs[0];
     np->childs[0]->matIndex = matStart;
     
     if (numAllocatedNodes >= numNodes)
@@ -195,7 +193,6 @@ void Tree::createSpecialBinaryTree(int numNodes, int numDirVars, const std::stri
     {
       np->childs[1] = new TreeNode(depth,1,np,1);
     }
-    _nodeMap[std::pair<int,int>(depth,1)]=np->childs[1];
     if (matStart == 0) matStart = 1;
     else matStart= 0;
     depth++;
@@ -207,8 +204,6 @@ void Tree::createSpecialTripleTree(int numNodes, int numDirVars, const std::stri
 {
   clear();
   _root = new TreeNode(0,0,NULL,0);  
-  _nodeMap[std::pair<int,int>(0,0)]=_root;
-  
   TreeNode* np = _root;
   int depth = 1;
   int numAllocatedNodes = 1;
@@ -227,16 +222,13 @@ void Tree::createSpecialTripleTree(int numNodes, int numDirVars, const std::stri
       np->direction.resize(numDirVars);
     }
     np->childs[0] = new TreeNode(depth,0,np,0,affunc);
-    _nodeMap[std::pair<int,int>(depth,0)]=np->childs[0];
     np->childs[0]->matIndex = 0;
     np->childs[1] = new TreeNode(depth,1,np,1,affunc);
-    _nodeMap[std::pair<int,int>(depth,1)]=np->childs[1];
     np->childs[1]->matIndex = 1;
     
     if (numAllocatedNodes < numNodes)
     {
       np->childs[2] = new TreeNode(depth,2,np,2);
-      _nodeMap[std::pair<int,int>(depth,2)]=np->childs[2];
       np = np->childs[2];
     }
     
@@ -328,7 +320,7 @@ void Tree::createRandomTreeSameDepthForNodes(int nbLeaves, int numChildMax, int
   
   auto printVector = [](const std::vector<int>& vv)
   {
-    printf("%d:",vv.size());
+    printf("%ld:",vv.size());
     for (int i=0; i< vv.size(); i++)
     {
       printf("%d ",vv[i]);
@@ -370,7 +362,6 @@ void Tree::createRandomTreeSameDepthForNodes(int nbLeaves, int numChildMax, int
   
   std::vector<TreeNode*> listPrev, listCur;  
   _root = new TreeNode(0,0,NULL,0);
-  _nodeMap[std::pair<int,int>(0,0)]=_root;
   
   listPrev.push_back(_root);
   int depth = 1;
@@ -400,7 +391,6 @@ void Tree::createRandomTreeSameDepthForNodes(int nbLeaves, int numChildMax, int
           np->childs[j] = new TreeNode(depth,inode,np,j);
           listCur.push_back(np->childs[j]);
         }
-        _nodeMap[std::pair<int,int>(depth,inode)]=np->childs[j];
         inode++;
       }
     }
@@ -542,10 +532,7 @@ void Tree::createRandomTree(int nbLeaves, int numChildMax, int leafPerNode, int
   std::vector<TreeNode*> listPrev, listCur;
   std::vector<int> listCurLeaves, listPrevLeaves;
 
-  
   _root = new TreeNode(0,0,NULL,0);
-  _nodeMap[std::pair<int,int>(0,0)]=_root;
-  
   listPrev.push_back(_root);
   listPrevLeaves.push_back(nbLeaves);
   
@@ -598,7 +585,6 @@ void Tree::createRandomTree(int nbLeaves, int numChildMax, int leafPerNode, int
           listCurLeaves.push_back(partIndexes[j]);
           listCur.push_back(np->childs[j]);
         }
-        _nodeMap[std::pair<int,int>(depth,inode)]=np->childs[j];
         inode++;
       }
     }
@@ -667,7 +653,6 @@ void Tree::createMixedTree(int maxDepth, int numChildMax, int leafPerNode, int n
     if (depth == 0)
     {
       _root = new TreeNode(0,0,NULL,0);
-      _nodeMap[std::pair<int,int>(0,0)]=_root;
       listPrev.push_back(_root);
     }
     else
@@ -699,7 +684,6 @@ void Tree::createMixedTree(int maxDepth, int numChildMax, int leafPerNode, int n
           {
             np->childs[j] = new TreeNode(depth,inode,np,j);
           }
-          _nodeMap[std::pair<int,int>(depth,inode)]=np->childs[j];
           listCur.push_back(np->childs[j]);
           inode++;
         }
@@ -729,7 +713,6 @@ void Tree::createPerfectTree(int maxDepth, int numChildRegular, int numChildLeaf
     if (depth == 0)
     {
       _root = new TreeNode(0,0,NULL,0);
-      _nodeMap[std::pair<int,int>(0,0)]=_root;
       listPrev.push_back(_root);
     }
     else
@@ -756,7 +739,6 @@ void Tree::createPerfectTree(int maxDepth, int numChildRegular, int numChildLeaf
           {
             np->childs[j] = new TreeNode(depth,inode,np,j);
           }
-          _nodeMap[std::pair<int,int>(depth,inode)]=np->childs[j];
           listCur.push_back(np->childs[j]);
           inode++;
         }
@@ -787,7 +769,6 @@ void Tree::createPerfectTree(int maxDepth, int numChilds, int numDirVars,const s
     if (depth == 0)
     {
       _root = new TreeNode(0,0,NULL,0);
-      _nodeMap[std::pair<int,int>(0,0)]=_root;
       listPrev.push_back(_root);
     }
     else
@@ -809,7 +790,6 @@ void Tree::createPerfectTree(int maxDepth, int numChilds, int numDirVars,const s
           {
             np->childs[j] = new TreeNode(depth,inode,np,j);
           }
-          _nodeMap[std::pair<int,int>(depth,inode)]=np->childs[j];
           listCur.push_back(np->childs[j]);
           inode++;
         }
@@ -865,7 +845,7 @@ void Tree::loadDataFromFile(std::string filename)
 {
   // clear tree
   clear();
-  _nodeMap.clear();
+  std::map<std::pair<int,int>, TreeNode* > nodeMap;
   printf("loading data from file %s\n",filename.c_str());
   //
   FILE* fp = fopen(filename.c_str(),"r");
@@ -889,20 +869,20 @@ void Tree::loadDataFromFile(std::string filename)
       std::vector<int> allChilds(numChild);
       for (int i=0; i< numChild; i++)
       {
-        fscanf(fp, "%d",&allChilds[i]);
+        int fl=fscanf(fp, "%d",&allChilds[i]);
         printf("child %d location %d\n",i,allChilds[i]);
       }
       double weight;
-      fscanf(fp, "%lf",&weight);
+      int fl=fscanf(fp, "%lf",&weight);
       printf("weight = %f\n",weight);
       std::vector<double> direction(numDir);
       for (int i=0; i< numDir; i++)
       {
-        fscanf(fp, "%lf",&direction[i]);
+        fl=fscanf(fp, "%lf",&direction[i]);
         printf("direction %d value %f\n",i,direction[i]);
       }
       char what[256];
-      fscanf(fp, "%s", what);
+      fl=fscanf(fp, "%s", what);
       printf("aftype = %s\n",what);
       
       // create node
@@ -926,7 +906,7 @@ void Tree::loadDataFromFile(std::string filename)
           node->childs.clear();
         }
         _root = node;
-        _nodeMap[std::pair<int,int>(depth,location)]=node;
+        nodeMap[std::pair<int,int>(depth,location)]=node;
       }
       else
       {
@@ -937,7 +917,7 @@ void Tree::loadDataFromFile(std::string filename)
         node->depth = depth;
         node->location = location;
         node->childOrder = childOrder;
-        node->parent = _nodeMap[std::pair<int,int>(depth-1,parentIdLoc)];
+        node->parent = nodeMap[std::pair<int,int>(depth-1,parentIdLoc)];
         node->parent->childs[childOrder] = node;
         node->matIndex = matIndex;
         if (numChild > 0)
@@ -950,7 +930,7 @@ void Tree::loadDataFromFile(std::string filename)
         }
         printf("node %d %d:\n",depth,location);
         node->parent->printData();
-        _nodeMap[std::pair<int,int>(depth,location)]=node;
+        nodeMap[std::pair<int,int>(depth,location)]=node;
       }
     }
     fclose(fp);
@@ -1027,7 +1007,7 @@ bool Tree::removeZeroLeaves(double tol, bool iteration)
         otherChild->parent = parent->parent;
         otherChild->childOrder = parent->childOrder;
         grandParent->childs[parent->childOrder] = otherChild;
-        Msg::Info("otherChild weight = %e parent weight = %e",otherChild->weight,parent->weight);
+        Msg::Info("child weight = %e otherChild weight = %e parent weight = %e",leaf->af->getVal(leaf->weight),otherChild->af->getVal(otherChild->weight),parent->weight);
       }
       else if (parent->childs.size() > 2)
       {
@@ -1263,6 +1243,15 @@ bool Tree::treeMerging()
                 }
               }
             };
+            
+            Msg::Info("\nSTART merging ");
+            parent->printData();
+            printf("merge loc (number of phases = %ld): %d (mat %d) ",parent->childs.size(),sameMat[0].second, (sameMat[0].first)->matIndex);
+            for (std::set<int>::iterator itpos = position.begin(); itpos != position.end(); itpos++)
+            {
+              printf(" %d (mat %d)",*itpos,parent->childs[*itpos]->matIndex);
+            }
+            printf("\n----------------------\n");
             parent->childs = newChilds;
             for (int j=0; j <parent->childs.size(); j++)
             {
@@ -1270,6 +1259,7 @@ bool Tree::treeMerging()
             }
             parent->direction = newNormals;
             parent->printData();
+            Msg::Info("DONE merging\n");
             //
             break;
           }
@@ -1281,6 +1271,110 @@ bool Tree::treeMerging()
   return false;
 };
 
+bool Tree::removeZeroBranches(double tol)
+{
+  int depth = 0;
+  bool ok = false;
+  while (true)
+  {
+    std::map<int,std::vector<TreeNode*> > allNodesByDepth;
+    getAllNodesByDepth(allNodesByDepth);
+    std::map<int,std::vector<TreeNode*> >::iterator itF= allNodesByDepth.find(depth);
+    if ( itF!= allNodesByDepth.end())
+    {
+      std::vector<TreeNode*>& depthNodes = itF->second;
+      for (int i=0; i<depthNodes.size(); i++)
+      {
+        TreeNode* node = depthNodes[i];
+        double Wt = getNonNegativeWeight(node);
+        //
+        bool removed = false;
+        std::vector<TreeNode*> newChilds;
+        std::vector<double> newDirection;
+        std::vector<int> removeLoc;
+        //
+        int numberNormal = node->childs.size()-1;
+        int totalNumberDirVars = node->direction.size();
+        int numberVarsPerNormal = totalNumberDirVars/numberNormal;
+        
+        for (int j=0; j< node->childs.size(); j++)
+        {
+          double Wchild = getNonNegativeWeight(node->childs[j]);
+          if (Wchild < tol* Wt)
+          {
+            removed = true;
+            removeLoc.push_back(j);
+          }
+          else
+          {
+            newChilds.push_back(node->childs[j]);
+            if (newChilds.size() > 1)
+            {
+              if (numberVarsPerNormal == 1)
+              {
+                newDirection.push_back(node->direction[j-1]);
+              }
+              else if (numberVarsPerNormal == 2)
+              {
+                newDirection.push_back(node->direction[(j-1)*numberVarsPerNormal]);
+                newDirection.push_back(node->direction[(j-1)*numberVarsPerNormal+1]);
+              }
+              else
+              {
+                Msg::Error("numberVarsPerNormal = %d is not correct !!!",numberVarsPerNormal);
+                Msg::Exit(0);
+              }
+            }
+          }
+        }
+        if (removed)
+        {
+          Msg::Info("\nSTART removing ");
+          node->printData();
+          ok = true;
+          printf("pos (number of phases= %ld): ",node->childs.size());
+          for (int j=0; j< removeLoc.size(); j++)
+          {
+            printf("%d (f=%e) ",removeLoc[j],getNonNegativeWeight(node->childs[removeLoc[j]])/Wt);
+          }
+          printf("\n");
+          if (newChilds.size() >1)
+          {
+            node->childs = newChilds;
+            for (int j=0; j <node->childs.size(); j++)
+            {
+              node->childs[j]->childOrder = j;
+            }
+            node->direction = newDirection;   
+            Msg::Info("------------------");
+            node->printData();     
+          }
+          else
+          {
+            TreeNode* grandParent = node->parent;
+            TreeNode* child = newChilds[0];
+            //
+            child->depth = node->depth;
+            child->location = node->location;
+            child->parent = node->parent;
+            child->childOrder = node->childOrder;
+            grandParent->childs[node->childOrder] = child;
+            Msg::Info("------------------");
+            Msg::Info("node is replaced by child");
+          }
+          Msg::Info("DONE removing\n");
+        }
+      }
+    }
+    else
+    {
+      break;
+    }
+    depth++;
+  };
+  return ok;
+};
+
 bool Tree::removeLeavesWithZeroContribution(double tol)
 {
   int iterRemove = 0;
@@ -1294,7 +1388,7 @@ bool Tree::removeLeavesWithZeroContribution(double tol)
       Msg::Error("removeLeavesWithZeroContribution: maximal number of iterations (%d) reaches", maxIterRemove);
       Msg::Exit(0);
     }
-    bool okLeaf = removeZeroLeaves(tol);
+    bool okLeaf = removeZeroBranches(tol);
     bool okMerge = treeMerging();
     // check
     if (okLeaf || okLeaf)
@@ -1396,7 +1490,7 @@ void Tree::printTree() const
   
   nodeContainer allLeaves;
   getAllLeaves(allLeaves);
-  printf("all leaves = %d\n",allLeaves.size());
+  printf("all leaves = %ld\n",allLeaves.size());
   for (int i=0; i< allLeaves.size(); i++)
   {
     allLeaves[i]->printData();
@@ -1412,7 +1506,7 @@ void Tree::printTreeInteraction(const std::string fname, bool colorMat, int dir)
   listPrev.push_back(_root);
   nodeContainer allLeaves;
   getAssociatedLeavesForNode(_root,allLeaves);
-  printf("*|%d|*\n",allLeaves.size());
+  printf("*|%ld|*\n",allLeaves.size());
   while (true)
   {
     listCur.clear();
@@ -1429,22 +1523,22 @@ void Tree::printTreeInteraction(const std::string fname, bool colorMat, int dir)
         {
           if (j==0)
           {
-            printf("%d(%d)",allLeaves.size(),nj->matIndex);
+            printf("%ld(%d)",allLeaves.size(),nj->matIndex);
           }
           else
           {
-            printf("-%d(%d)",allLeaves.size(),nj->matIndex);
+            printf("-%ld(%d)",allLeaves.size(),nj->matIndex);
           }
         }
         else
         {
           if (j==0)
           {
-            printf("%d",allLeaves.size());
+            printf("%ld",allLeaves.size());
           }
           else
           {
-            printf("-%d",allLeaves.size());
+            printf("-%ld",allLeaves.size());
           }
           listCur.push_back(nj);
         }
@@ -1596,11 +1690,11 @@ void Tree::getRefToAllLeaves(std::vector<TreeNode*>& allLeaves)
   };
 };
 
-void Tree::getAllNodesByDepth(std::map<int,nodeContainer >& allNodes) const
+void Tree::getAllNodesByDepth(std::map<int,std::vector<TreeNode*> >& allNodes)
 {
   allNodes.clear();
   if (_root == NULL) return;
-  nodeContainer listPrev, listCur;
+  std::vector<TreeNode*> listPrev, listCur;
   listPrev.push_back(_root);
   //
   allNodes[_root->depth].push_back(_root);
@@ -1609,7 +1703,7 @@ void Tree::getAllNodesByDepth(std::map<int,nodeContainer >& allNodes) const
     listCur.clear();    
     for (int i =0; i< listPrev.size(); i++)
     {
-      const TreeNode* np = listPrev[i];
+      TreeNode* np = listPrev[i];
       for (int j=0; j< np->childs.size(); j++)
       {
         listCur.push_back(np->childs[j]);
@@ -1624,24 +1718,42 @@ void Tree::getAllNodesByDepth(std::map<int,nodeContainer >& allNodes) const
   };
 };
 
-const TreeNode* Tree::getNode(int depth, int pos) const
+void Tree::getAllNodesByDepth(std::map<int,nodeContainer >& allNodes) const
+{
+  allNodes.clear();
+  if (_root == NULL) return;
+  nodeContainer listPrev, listCur;
+  listPrev.push_back(_root);
+  //
+  allNodes[_root->depth].push_back(_root);
+  while (true)
+  {
+    listCur.clear();    
+    for (int i =0; i< listPrev.size(); i++)
+    {
+      const TreeNode* np = listPrev[i];
+      for (int j=0; j< np->childs.size(); j++)
       {
-  std::map<std::pair<int,int>, TreeNode* >::const_iterator it = _nodeMap.find(std::pair<int,int>(depth,pos));
-  if (it == _nodeMap.end())
+        listCur.push_back(np->childs[j]);
+        allNodes[np->childs[j]->depth].push_back(np->childs[j]);
+      }
+    }
+    listPrev = listCur;
+    if (listCur.size() == 0)
     {
-    printf("node in depth %d pos %d is not found \n",depth,pos);
-    return NULL;
+      break;
     }
-  return it->second;
+  };
 };
 
+
 void Tree::printBackTrackingPath(const TreeNode* node) const
 {
   nodeContainer allNodes;
   getBackTrackingPath(node,allNodes);
   printf("-----------------------------\n");
   node->printData();
-  printf("backtracking laves to node: %d \n",allNodes.size());
+  printf("backtracking laves to node: %ld \n",allNodes.size());
   for (int i=0; i< allNodes.size(); i++)
   {
     allNodes[i]->printData();
@@ -1656,7 +1768,7 @@ void Tree::printAssociatedLeavesForNode(const TreeNode* node) const
   
   printf("-----------------------------\n");
   node->printData();
-  printf("associated laves to node: %d \n",allLeaves.size());
+  printf("associated laves to node: %ld \n",allLeaves.size());
   for (int i=0; i< allLeaves.size(); i++)
   {
     allLeaves[i]->printData();
@@ -1676,6 +1788,54 @@ void Tree::getBackTrackingPath(const TreeNode* node, nodeContainer& allNodes) co
   };
 };
 
+void Tree::getAssociatedLeavesForNode(TreeNode* node, std::vector<TreeNode*>& allLeaves)
+{
+  allLeaves.clear();
+  std::vector<TreeNode*> listPrev, listCur;
+  listPrev.push_back(node);
+  while (true)
+  {
+    listCur.clear();    
+    for (int i =0; i< listPrev.size(); i++)
+    {
+      TreeNode* np = listPrev[i];
+      if (np->childs.size() == 0)
+      {
+        allLeaves.push_back(np);
+      }
+      for (int j=0; j< np->childs.size(); j++)
+      {
+        listCur.push_back(np->childs[j]);
+      }
+    }
+    listPrev = listCur;
+    if (listCur.size() == 0)
+    {
+      break;
+    }
+  };
+};
+
+double Tree::getNonNegativeWeight(const TreeNode* node) const
+{
+  if (node->childs.size() == 0)
+  {
+    return node->af->getVal(node->weight);
+  }
+  else
+  {
+    nodeContainer allLeaves;
+    getAssociatedLeavesForNode(node,allLeaves);
+    double W = 0;
+    for (int i=0; i< allLeaves.size(); i++)
+    {
+      const TreeNode* leaf = allLeaves[i];
+      W += leaf->af->getVal(leaf->weight);
+    };
+    return W;
+  }
+}
+
 void Tree::getAssociatedLeavesForNode(const TreeNode* node, nodeContainer& allLeaves) const
 {
   allLeaves.clear();
@@ -1763,7 +1923,7 @@ void Tree::printAllRegularNodes() const
   nodeContainer allNodes;
   getAllRegularNodes(allNodes);
   printf("-----------------------------\n");
-  printf("all regular nodes: %d \n",allNodes.size());
+  printf("all regular nodes: %ld \n",allNodes.size());
   for (int i=0; i< allNodes.size(); i++)
   {
     allNodes[i]->printData();
@@ -1926,7 +2086,7 @@ void Tree::initialize(bool rand)
   //
   std::vector<TreeNode*> allNodes;
   getRefToAllNodes(allNodes);
-  printf("num of nodes = %d\n",allNodes.size());
+  printf("num of nodes = %ld\n",allNodes.size());
   for (int i=0; i< allNodes.size(); i++)
   {
     TreeNode& node = *(allNodes[i]);

NonLinearSolver/modelReduction/Tree.h

@@ -61,7 +61,6 @@ class Tree
     void clear();
     void initialize(bool rand=true);
     const TreeNode* getRootNode() const {return _root;};
-    const TreeNode* getNode(int depth, int pos) const;
     
     void assignMaterialLaws(int numPhases);
     
@@ -78,6 +77,7 @@ class Tree
     void printNodeFraction() const;
     
     bool removeLeavesWithZeroContribution(double tol=1e-6);
+    bool removeZeroBranches(double tol);
     bool removeZeroLeaves(double tol, bool iteration=true);
     bool treeMerging();
     
@@ -86,9 +86,11 @@ class Tree
     void getAllLeaves(nodeContainer& allLeaves) const;
     void getAllNodes(nodeContainer& allNodes) const;
     void getAllNodesByDepth(std::map<int,nodeContainer>& allNodes) const;
+    void getAllNodesByDepth(std::map<int,std::vector<TreeNode*> >& allNodes);
     void getAssociatedLeavesForNode(const TreeNode* node, nodeContainer& allLeaves) const;
+    void getAssociatedLeavesForNode(TreeNode* node, std::vector<TreeNode*>& allLeaves);
     void getBackTrackingPath(const TreeNode* node, nodeContainer& path) const;
-    
+    double getNonNegativeWeight(const TreeNode* node) const;
     void getRefToAllLeaves(std::vector<TreeNode*>& allLeaves);
     void getRefToAllNodes(std::vector<TreeNode*>& allNodes);  
     #endif //SWIG
@@ -96,7 +98,6 @@ class Tree
   protected:
     #ifndef SWIG
     TreeNode* _root;
-    std::map<std::pair<int,int>, TreeNode* > _nodeMap;
     #endif //SWIG
 };
 #endif //_TREE_H_