ext/OGDF/src/energybased/LinearQuadtreeExpansion.cpp

   1 /*
   2  * $Revision: 2552 $
   3  *
   4  * last checkin:
   5  *   $Author: gutwenger $
   6  *   $Date: 2012-07-05 16:45:20 +0200 (Do, 05. Jul 2012) $
   7  ***************************************************************/
   8
   9 /** \file
  10  * \brief Implementation of class LinearQuadtreeExpansion.
  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 #include "LinearQuadtreeExpansion.h"
  44 #include "ComplexDouble.h"
  45 #include "WSPD.h"
  46 #include <complex>
  47
  48 using namespace ogdf::sse;
  49
  50 namespace ogdf {
  51
  52 LinearQuadtreeExpansion::LinearQuadtreeExpansion(__uint32 precision, const LinearQuadtree& tree) : m_tree(tree), m_numCoeff(precision), binCoef(2*m_numCoeff)
  53 {
  54         m_numExp = m_tree.maxNumberOfNodes();
  55         allocate();
  56 }
  57
  58
  59 LinearQuadtreeExpansion::~LinearQuadtreeExpansion(void)
  60 {
  61         deallocate();
  62 }
  63
  64
  65 void LinearQuadtreeExpansion::allocate()
  66 {
  67         m_multiExp = (double*)MALLOC_16(m_numCoeff*sizeof(double)*2*m_numExp);
  68         m_localExp = (double*)MALLOC_16(m_numCoeff*sizeof(double)*2*m_numExp);
  69 }
  70
  71
  72 void LinearQuadtreeExpansion::deallocate()
  73 {
  74         FREE_16(m_multiExp);
  75         FREE_16(m_localExp);
  76 }
  77
  78
  79 void LinearQuadtreeExpansion::P2M(__uint32 point, __uint32 receiver)
  80 {
  81         double* receiv_coeff = m_multiExp + receiver*(m_numCoeff<<1);
  82         const double q = (double)m_tree.pointSize(point);
  83         const double x = (double)m_tree.pointX(point);
  84         const double y = (double)m_tree.pointY(point);
  85         const double centerX = (double)m_tree.nodeX(receiver);
  86         const double centerY = (double)m_tree.nodeY(receiver);
  87         // a0 += q_i
  88         receiv_coeff[0] += q;
  89         // a_1..m
  90         ComplexDouble ak;
  91         // p - z0
  92         ComplexDouble delta(ComplexDouble(x, y) - ComplexDouble(centerX, centerY));
  93         // (p - z0)^k
  94         ComplexDouble delta_k(delta);
  95         for (__uint32 k=1; k<m_numCoeff; k++)
  96         {
  97                 ak.load(receiv_coeff+(k<<1));
  98                 ak -= delta_k*(q/(double)k);
  99                 ak.store(receiv_coeff+(k<<1));
 100                 delta_k *= delta;
 101         }
 102 }
 103
 104
 105 void LinearQuadtreeExpansion::L2P(__uint32 source, __uint32 point, float& fx, float& fy)
 106 {
 107         const double* source_coeff = m_localExp + source*(m_numCoeff << 1);
 108         const double x = (double)m_tree.pointX(point);
 109         const double y = (double)m_tree.pointY(point);
 110         const double centerX = (double)m_tree.nodeX(source);
 111         const double centerY = (double)m_tree.nodeY(source);
 112         ComplexDouble ak;
 113         ComplexDouble res;
 114         ComplexDouble delta(ComplexDouble(x,y) - ComplexDouble(centerX, centerY));
 115         ComplexDouble delta_k(1);
 116         for (__uint32 k=1; k<m_numCoeff; k++)
 117         {
 118                 ak.load(source_coeff+(k<<1));
 119                 res += ak*delta_k*(double)k;
 120                 delta_k *= delta;
 121         }
 122         res = res.conj();
 123         double resTemp[2];
 124         res.store_unaligned(resTemp);
 125
 126         fx -= ((float)resTemp[0]);
 127         fy -= ((float)resTemp[1]);
 128 }
 129
 130
 131 void LinearQuadtreeExpansion::M2M(__uint32 source, __uint32 receiver)
 132 {
 133         double* receiv_coeff = m_multiExp + receiver*(m_numCoeff<<1);
 134         double* source_coeff = m_multiExp + source*(m_numCoeff<<1);
 135
 136         const double center_x_source   = (double)m_tree.nodeX(source);
 137         const double center_y_source   = (double)m_tree.nodeY(source);
 138         const double center_x_receiver = (double)m_tree.nodeX(receiver);
 139         const double center_y_receiver = (double)m_tree.nodeY(receiver);
 140
 141         ComplexDouble delta(ComplexDouble(center_x_source, center_y_source) - ComplexDouble(center_x_receiver, center_y_receiver));
 142
 143         ComplexDouble a(source_coeff);
 144         ComplexDouble b(receiv_coeff);
 145         b += a;
 146         b.store(receiv_coeff);
 147         for (__uint32 l=1;l<m_numCoeff;l++)
 148         {
 149                 b.load(receiv_coeff+(l<<1));
 150                 ComplexDouble delta_k(1.0,0.0);
 151                 for (__uint32 k=0;k<l;k++)
 152                 {
 153                         a.load(source_coeff+((l-k)<<1));
 154                         b += a*delta_k*binCoef.value(l-1, k);
 155                         delta_k *= delta;
 156                 }
 157                 a.load(source_coeff);
 158                 b -= a*delta_k*(1/(double)l);
 159                 b.store(receiv_coeff+(l<<1));
 160         }
 161 }
 162
 163
 164 void LinearQuadtreeExpansion::L2L(__uint32 source, __uint32 receiver)
 165 {
 166         double* receiv_coeff = m_localExp + receiver*(m_numCoeff<<1);
 167         double* source_coeff = m_localExp + source*(m_numCoeff<<1);
 168
 169         const double center_x_source   = (double)m_tree.nodeX(source);
 170         const double center_y_source   = (double)m_tree.nodeY(source);
 171         const double center_x_receiver = (double)m_tree.nodeX(receiver);
 172         const double center_y_receiver = (double)m_tree.nodeY(receiver);
 173
 174         ComplexDouble center_receiver(center_x_receiver, center_y_receiver);
 175         ComplexDouble center_source(center_x_source, center_y_source);
 176         ComplexDouble delta(center_source - center_receiver);
 177
 178         for (__uint32 l=0;l<m_numCoeff;l++)
 179         {
 180                 ComplexDouble b(receiv_coeff+(l<<1));
 181                 ComplexDouble delta_k(1.0,0.0);
 182                 for (__uint32 k=l;k<m_numCoeff;k++)
 183                 {
 184                         ComplexDouble a(source_coeff+(k<<1));
 185                         b += a*delta_k*binCoef.value(k, l);
 186                         delta_k *= delta;
 187                 }
 188                 b.store(receiv_coeff+(l<<1));
 189         }
 190 }
 191
 192
 193 void LinearQuadtreeExpansion::M2L(__uint32 source, __uint32 receiver)
 194 {
 195         double* receiv_coeff = m_localExp + receiver*(m_numCoeff<<1);
 196         double* source_coeff = m_multiExp + source*(m_numCoeff<<1);
 197
 198         const float center_x_source   = (float)m_tree.nodeX(source);
 199         const float center_y_source   = (float)m_tree.nodeY(source);
 200         const float center_x_receiver = (float)m_tree.nodeX(receiver);
 201         const float center_y_receiver = (float)m_tree.nodeY(receiver);
 202
 203         ComplexDouble center_receiver(center_x_receiver, center_y_receiver);
 204         ComplexDouble center_source(center_x_source, center_y_source);
 205         ComplexDouble delta0(center_source - center_receiver);
 206
 207         ComplexDouble delta1 = -delta0;
 208         ComplexDouble delta1_l(delta1);
 209         ComplexDouble a;
 210         ComplexDouble a0(source_coeff);
 211         ComplexDouble b;
 212         ComplexDouble sum;
 213         for (__uint32 l=1;l<m_numCoeff;l++)
 214         {
 215                 b.load(receiv_coeff+(l<<1));
 216                 sum = a0*(-1/(double)l);
 217                 ComplexDouble delta0_k(delta0);
 218                 for (__uint32 k=1;k<m_numCoeff;k++)
 219                 {
 220                         a.load(source_coeff+(k<<1));
 221                         sum += (a*binCoef.value(l+k-1, k-1)) / delta0_k;
 222                         delta0_k *= delta0;
 223                 }
 224
 225                 b += sum/delta1_l;
 226                 b.store(receiv_coeff+(l<<1));
 227                 delta1_l *= delta1;
 228         }
 229
 230         // b0
 231         b.load(receiv_coeff);
 232         //std::complex<double> sdelta(m_center[(receiver<<1)] - m_center[(source<<1)],  m_center[(receiver<<1)+1] - m_center[(source<<1) +1]);//, m_x[receiver] - m_x[source], m_y[receiver] - m_y[source]);
 233         /*std::complex<double> sdelta(center_x_receiver - center_x_source,
 234                                                                 center_y_receiver - center_y_source);//, m_x[receiver] - m_x[source], m_y[receiver] - m_y[source]);
 235
 236         if ((sdelta.real() <=0) && (sdelta.imag() == 0)) //no cont. compl. log fct exists !!!
 237                 sdelta = log(sdelta + 0.00000001);
 238         else
 239                 sdelta = log(sdelta);
 240
 241         */
 242         double r = delta1.length();//= sqrt(sdelta.real()*sdelta.real()+sdelta.imag()*sdelta.imag());
 243         double phi = atan((center_x_receiver - center_x_source)/(center_y_receiver - center_y_source));
 244         // sum = a0*log(z1 - z0)
 245         b += a0*ComplexDouble(log(r), phi);
 246         // (z1 - z0)^1
 247         //b += a0*ComplexDouble(sdelta.real(), sdelta.imag());
 248         delta1_l = delta1;
 249         for (__uint32 k=1;k<m_numCoeff;k++)
 250         {
 251                 a.load(source_coeff+(k<<1));
 252                 b += a/delta1_l;
 253                 //(z1 - z0)^k
 254                 delta1_l *= delta1;
 255         }
 256         b.store(receiv_coeff);
 257 }
 258
 259 } // end of namespace ogdf