ext/OGDF/src/energybased/FMEKernel.h

   1 /*
   2  * $Revision: 2565 $
   3  *
   4  * last checkin:
   5  *   $Author: gutwenger $
   6  *   $Date: 2012-07-07 17:14:54 +0200 (Sa, 07. Jul 2012) $
   7  ***************************************************************/
   8
   9 /** \file
  10  * \brief Declaration of FME kernel.
  11  *
  12  * \author Martin Gronemann
  13  *
  14  * \par License:
  15  * This file is part of the Open Graph Drawing Framework (OGDF).
  16  *
  17  * \par
  18  * Copyright (C)<br>
  19  * See README.txt in the root directory of the OGDF installation for details.
  20  *
  21  * \par
  22  * This program is free software; you can redistribute it and/or
  23  * modify it under the terms of the GNU General Public License
  24  * Version 2 or 3 as published by the Free Software Foundation;
  25  * see the file LICENSE.txt included in the packaging of this file
  26  * for details.
  27  *
  28  * \par
  29  * This program is distributed in the hope that it will be useful,
  30  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  31  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  32  * GNU General Public License for more details.
  33  *
  34  * \par
  35  * You should have received a copy of the GNU General Public
  36  * License along with this program; if not, write to the Free
  37  * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  38  * Boston, MA 02110-1301, USA.
  39  *
  40  * \see  http://www.gnu.org/copyleft/gpl.html
  41  ***************************************************************/
  42
  43 #ifdef _MSC_VER
  44 #pragma once
  45 #endif
  46
  47 #ifndef OGDF_FME_KERNEL_H
  48 #define OGDF_FME_KERNEL_H
  49
  50 #include <ogdf/basic/basic.h>
  51 #include "FastUtils.h"
  52 #include "ArrayGraph.h"
  53 #include "FMEThread.h"
  54 #include <list>
  55
  56 namespace ogdf {
  57
  58 class FMEKernel
  59 {
  60 public:
  61         FMEKernel(FMEThread* pThread) : m_pThread(pThread) { }
  62
  63         inline void sync() { m_pThread->sync(); }
  64
  65         //! returns the index of the thread ( 0.. numThreads()-1 )
  66         inline __uint32 threadNr() const { return m_pThread->threadNr(); }
  67
  68         //! returns the total number of threads in the pool
  69         inline __uint32 numThreads() const { return m_pThread->numThreads(); }
  70
  71         //! returns true if this is the main thread ( the main thread is always the first thread )
  72         inline bool isMainThread() const { return m_pThread->isMainThread(); }
  73
  74         //! returns true if this run only uses one thread )
  75         inline bool isSingleThreaded() const { return (m_pThread->numThreads() == 1); };
  76
  77 private:
  78         FMEThread* m_pThread;
  79 };
  80
  81
  82 #define FME_KERNEL_USE_OLD
  83
  84 #define COMPUTE_FORCE_PROTECTION_FACTOR 0.25f
  85 // makro for force computation via SSE   s / max(s*0.5, (dx*dx + dy*dy))
  86 #define _MM_COMPUTE_FORCE(dx,dy,s) _mm_div_ps((s),_mm_max_ps(_mm_mul_ps((s),_mm_set1_ps(COMPUTE_FORCE_PROTECTION_FACTOR)), _mm_add_ps(_mm_mul_ps((dx),(dx)), _mm_mul_ps((dy),(dy)))))
  87 #define COMPUTE_FORCE(dx,dy,s) (s/(max<float>(s*COMPUTE_FORCE_PROTECTION_FACTOR, (dx)*(dx) + (dy)*(dy))))
  88
  89 inline double move_nodes(float* x, float* y, const __uint32 begin, const __uint32 end, const float* fx, const float* fy, const float t)
  90 {
  91         double dsq_max = 0.0;
  92         for (__uint32 i=begin; i <= end; i++)
  93         {
  94                 double dsq = fx[i]*fx[i] + fy[i]*fy[i];
  95                 x[i] += fx[i]*t;
  96                 y[i] += fy[i]*t;
  97                 dsq_max = max(dsq_max, dsq);
  98         }
  99         return dsq_max;
 100 }
 101
 102
 103 inline void eval_edges(const ArrayGraph& graph, const __uint32 begin, const __uint32 end, float* fx, float* fy)
 104 {
 105         const float* x = graph.nodeXPos();
 106         const float* y = graph.nodeYPos();
 107         const float* e = graph.desiredEdgeLength();
 108         for (__uint32 i=begin; i <= end; i++)
 109         {
 110                 const EdgeAdjInfo& e_info = graph.edgeInfo(i);
 111                 const __uint32 a = e_info.a;
 112                 const __uint32 b = e_info.b;
 113                 const NodeAdjInfo& a_info = graph.nodeInfo(a);
 114                 const NodeAdjInfo& b_info = graph.nodeInfo(b);
 115
 116                 float dx = x[a] - x[b];
 117                 float dy = y[a] - y[b];
 118                 float dsq = dx*dx + dy*dy;
 119                 float f = (logf(dsq)*0.5f-logf(e[i])) * 0.25f;
 120                 float fa = (float)(f/((float)a_info.degree));
 121                 float fb = (float)(f/((float)b_info.degree));
 122                 fx[a] -= dx*fa;
 123                 fy[a] -= dy*fa;
 124                 fx[b] += dx*fb;
 125                 fy[b] += dy*fb;
 126         }
 127 }
 128
 129
 130 //! kernel function to evaluate forces between n points with coords x, y directly. result is stored in fx, fy
 131 inline void eval_direct(float* x, float* y, float* s, float* fx, float* fy, size_t n)
 132 {
 133         for (__uint32 i=0; i < n; i++)
 134         {
 135                 for (__uint32 j=i+1; j < n; j++)
 136                 {
 137                         float dx = x[i] - x[j];
 138                         float dy = y[i] - y[j];
 139 #ifdef FME_KERNEL_USE_OLD
 140                         float s_sum = s[i]+s[j];
 141 #else
 142                         float s_sum = s[i]*s[j];
 143 #endif
 144                         float f = COMPUTE_FORCE(dx, dy, s_sum);
 145                         fx[i] += dx*f;
 146                         fy[i] += dy*f;
 147                         fx[j] -= dx*f;
 148                         fy[j] -= dy*f;
 149                 }
 150         }
 151 }
 152
 153
 154 //! kernel function to evaluate forces between two sets of points with coords x1, y1 (x2, y2) directly. result is stored in fx1, fy1 (fx2, fy2
 155 inline void eval_direct(float* x1, float* y1, float* s1, float* fx1, float* fy1, size_t n1,
 156                                                 float* x2, float* y2, float* s2, float* fx2, float* fy2, size_t n2)
 157 {
 158         for (__uint32 i=0; i < n1; i++)
 159         {
 160                 for (__uint32 j=0; j < n2; j++)
 161                 {
 162                         float dx = x1[i] - x2[j];
 163                         float dy = y1[i] - y2[j];
 164 #ifdef FME_KERNEL_USE_OLD
 165                         float s_sum = s1[i]+s2[j];
 166 #else
 167                         float s_sum = s1[i]*s2[j];
 168 #endif
 169                         float f = COMPUTE_FORCE(dx, dy, s_sum);
 170                         fx1[i] += dx*f;
 171                         fy1[i] += dy*f;
 172                         fx2[j] -= dx*f;
 173                         fy2[j] -= dy*f;
 174                 }
 175         }
 176 }
 177
 178
 179 #ifndef OGDF_FME_KERNEL_USE_SSE_DIRECT
 180 //! kernel function to evaluate forces between n points with coords x, y directly. result is stored in fx, fy
 181 inline void eval_direct_fast(float* x, float* y, float* s, float* fx, float* fy, size_t n)
 182 {
 183         eval_direct(x, y, s, fx, fy, n);
 184 }
 185
 186 //! kernel function to evaluate forces between two sets of points with coords x1, y1 (x2, y2) directly. result is stored in fx1, fy1 (fx2, fy2
 187 inline void eval_direct_fast(float* x1, float* y1, float* s1, float* fx1, float* fy1, size_t n1,
 188                                                          float* x2, float* y2, float* s2, float* fx2, float* fy2, size_t n2)
 189 {
 190         eval_direct(x1, y1, s1, fx1, fy1, n1,
 191                                 x2, y2, s2, fx2, fy2, n2);
 192 }
 193
 194 #else
 195 //! kernel function to evaluate forces between n points with coords x, y directly. result is stored in fx, fy
 196 void eval_direct_fast(float* x, float* y, float* s, float* fx, float* fy, size_t n);
 197 //! kernel function to evaluate forces between two sets of points with coords x1, y1 (x2, y2) directly. result is stored in fx1, fy1 (fx2, fy2
 198 void eval_direct_fast(
 199         float* x1, float* y1, float* s1, float* fx1, float* fy1, size_t n1,
 200         float* x2, float* y2, float* s2, float* fx2, float* fy2, size_t n2);
 201 #endif
 202
 203 //! kernel function to evalute a local expansion at point x,y result is added to fx, fy
 204 void fast_multipole_l2p(double* localCoeffiecients, __uint32 numCoeffiecients, double centerX, double centerY,
 205                                                                         float x, float y, float q, float& fx, float &fy);
 206
 207 void fast_multipole_p2m(double* mulitCoeffiecients, __uint32 numCoeffiecients, double centerX, double centerY,
 208                                                                         float x, float y, float q);
 209
 210 class FMEBasicKernel
 211 {
 212 public:
 213         inline void edgeForces(const ArrayGraph& graph, float* fx, float* fy)
 214         {
 215                 eval_edges(graph, 0, graph.numEdges()-1, fx, fy);
 216         }
 217
 218         inline void repForces(ArrayGraph& graph, float* fx, float* fy)
 219         {
 220                 eval_direct_fast(graph.nodeXPos(), graph.nodeYPos(), graph.nodeSize(), fx, fy,  graph.numNodes());
 221         }
 222
 223         inline double moveNodes(ArrayGraph& graph, float* fx, float* fy, float timeStep)
 224         {
 225                 return move_nodes(graph.nodeXPos(), graph.nodeYPos(), 0, graph.numNodes()-1, fx, fy, timeStep);
 226         }
 227
 228         inline double simpleIteration(ArrayGraph& graph, float* fx, float* fy, float timeStep)
 229         {
 230                 repForces(graph, fx, fy);
 231                 edgeForces(graph, fx, fy);
 232                 return moveNodes(graph, fx, fy, timeStep);
 233         }
 234
 235
 236         inline double simpleEdgeIteration(ArrayGraph& graph, float* fx, float* fy, float timeStep)
 237         {
 238                 edgeForces(graph, fx, fy);
 239                 return moveNodes(graph, fx, fy, timeStep);
 240         }
 241
 242
 243         inline void simpleForceDirected(ArrayGraph& graph, float timeStep, __uint32 minIt, __uint32 maxIt, __uint32 preProcIt, double threshold)
 244         {
 245                 bool earlyExit = false;
 246                 float* fx = (float*)MALLOC_16(sizeof(float)*graph.numNodes());
 247                 float* fy = (float*)MALLOC_16(sizeof(float)*graph.numNodes());
 248
 249                 for (__uint32 i = 0; i<preProcIt; i++)
 250                 {
 251                         for (__uint32 j = 0; j<graph.numNodes(); j++)
 252                         {
 253                                 fx[j] = 0.0f;
 254                                 fy[j] = 0.0f;
 255                         }
 256                         simpleEdgeIteration(graph, fx, fy, timeStep);
 257                 }
 258                 for(__uint32 i = 0; (i<maxIt) && (!earlyExit); i++)
 259                 {
 260                         for (__uint32 j = 0; j<graph.numNodes(); j++)
 261                         {
 262                                 fx[j] = 0.0f;
 263                                 fy[j] = 0.0f;
 264                         }
 265                         double dsq = simpleIteration(graph, fx, fy, timeStep);
 266                         if (dsq < threshold && i>(minIt))
 267                                 earlyExit = true;
 268                 }
 269
 270                 FREE_16(fx);
 271                 FREE_16(fy);
 272         }
 273
 274 private:
 275 };
 276
 277
 278 class FMESingleKernel : FMEBasicKernel
 279 {
 280 public:
 281         //FMESingleKernel(FMEThread* pThread) : FMEKernel(pThread) {};
 282         void operator()(ArrayGraph& graph, float timeStep, __uint32 minIt, __uint32 maxIt, double threshold)
 283         {
 284                 simpleForceDirected(graph, timeStep, minIt, maxIt, 20, threshold);
 285         }
 286 };
 287
 288 } // end of namespace ogdf
 289
 290 #endif