Project

General

Profile

Statistics
| Branch: | Revision:

root / src / layout / forcedirectedembedding.cc @ 8aeaaccc

History | View | Annotate | Download (13.1 KB)

1
//=========================================================================
2
//  FORCEDIRECTEDEMBEDDING.H - part of
3
//                  OMNeT++/OMNEST
4
//           Discrete System Simulation in C++
5
//
6
//  Author: Levente Meszaros
7
//
8
//=========================================================================
9

    
10
/*--------------------------------------------------------------*
11
  Copyright (C) 2006-2008 OpenSim Ltd.
12

13
  This file is distributed WITHOUT ANY WARRANTY. See the file
14
  `license' for details on this and other legal matters.
15
*--------------------------------------------------------------*/
16

    
17
#include <float.h>
18
#include "forcedirectedembedding.h"
19
#include "forcedirectedparameters.h"
20
#include "lcgrandom.h"
21

    
22
USING_NAMESPACE
23

    
24
ForceDirectedEmbedding::ForceDirectedEmbedding() {
25
    debugLevel = 0;
26
    inspected = false;
27
    initialized = false;
28
    finished = false;
29

    
30
    parameters = getParameters();
31
}
32

    
33
ForceDirectedEmbedding::~ForceDirectedEmbedding() {
34
    for(std::vector<Variable *>::iterator it = variables.begin(); it != variables.end(); it++)
35
        delete *it;
36

    
37
    for(std::vector<IForceProvider *>::iterator it = forceProviders.begin(); it != forceProviders.end(); it++)
38
        delete *it;
39

    
40
    for(std::vector<IBody *>::iterator it = bodies.begin(); it != bodies.end(); it++)
41
        delete *it;
42
}
43

    
44
ForceDirectedParameters ForceDirectedEmbedding::getParameters(int32 seed) {
45
    LCGRandom lcgRandom(seed);
46

    
47
    ForceDirectedParameters parameters;
48

    
49
    parameters.defaultBodySize = Rs(10, 10);
50
    parameters.defaultBodyMass = 10;
51
    parameters.defaultBodyCharge = 1;
52

    
53
    parameters.defaultSpringCoefficient = 0.1 + 0.9 * lcgRandom.next01();
54
    parameters.defaultSpringReposeLength = 50;
55

    
56
    parameters.electricRepulsionCoefficient = 10000 + 90000 * lcgRandom.next01();
57
    parameters.defaultElectricRepulsionLinearityDistance = -1;
58
    parameters.defaultElectricRepulsionMaxDistance = -1;
59

    
60
    parameters.frictionCoefficient = 1 + 4 * lcgRandom.next01();
61

    
62
    parameters.defaultSlippery = false;
63
    parameters.defaultPointLikeDistance = false;
64

    
65
    parameters.timeStep = 1;
66
    parameters.minTimeStep = 0;
67
    parameters.maxTimeStep = DBL_MAX;
68
    parameters.timeStepMultiplier = 2;
69

    
70
    parameters.minAccelerationError = 0.2;
71
    parameters.maxAccelerationError = 0.5;
72

    
73
    parameters.defaultMaxForce = 1000;
74
    parameters.maxVelocity = 100;
75

    
76
    parameters.velocityRelaxLimit = 0.2;
77
    parameters.accelerationRelaxLimit = 1;
78

    
79
    parameters.maxCycle = 1000;
80
    parameters.maxCalculationTime = 1000 + 9000 * lcgRandom.next01();
81

    
82
    return parameters;
83
}
84

    
85
void ForceDirectedEmbedding::reinitialize() {
86
    initialized = true;
87
    finished = false;
88

    
89
    relaxFactor = 0;
90
    cycle = 0;
91
    probCycle = 0;
92
    elapsedTime = 0;
93
    elapsedTicks = 0;
94
    elapsedCalculationTime = 0;
95
    kineticEnergySum = 0;
96
    totalMass = 0;
97

    
98
    // internal state arrays
99
    pn = createPtArray();
100
    vn = createPtArray();
101
    an = createPtArray();
102
    a1 = createPtArray();
103
    a2 = createPtArray();
104
    a3 = createPtArray();
105
    a4 = createPtArray();
106
    dpn = createPtArray();
107
    tpn = createPtArray();
108
    dvn = createPtArray();
109
    tvn = createPtArray();
110

    
111
    // reinitialize parts
112
    for(std::vector<Variable *>::iterator it = variables.begin(); it != variables.end(); it++)
113
        (*it)->reinitialize();
114
    for(std::vector<IForceProvider *>::iterator it = forceProviders.begin(); it != forceProviders.end(); it++)
115
        (*it)->reinitialize();
116
    for(std::vector<IBody *>::iterator it = bodies.begin(); it != bodies.end(); it++)
117
        (*it)->reinitialize();
118

    
119
    // reinitialize positions and velocities
120
    for (int i = 0; i < (int)variables.size(); i++) {
121
        Variable *variable = variables[i];
122
        pn[i].assign(variable->getPosition());
123
        vn[i].assign(variable->getVelocity());
124
        double variableMass = 0;
125
        for(std::vector<IBody *>::iterator it = bodies.begin(); it != bodies.end(); it++) {
126
            IBody *body = *it;
127
            if (body->getVariable() == variable) {
128
                variableMass += body->getMass();
129
            }
130
        }
131
        variable->setMass(variableMass);
132
        totalMass += variable->getMass();
133
    }
134

    
135
    // some values needs to be reinitialized
136
    updatedTimeStep = parameters.timeStep;
137

    
138
    if (parameters.frictionCoefficient == -1)
139
        parameters.frictionCoefficient = getEnergyBasedFrictionCoefficient(parameters.maxCalculationTime);
140
}
141

    
142
/**
143
 * Find the solution for the differential equation using a
144
 * modified Runge-Kutta 4th order algorithm.
145
 *
146
 * a1 = a[pn, vn]
147
 * a2 = a[pn + h / 2 * vn + h * h / 8 * a1, vn + h / 2 * a1]
148
 * a3 = a[pn + h / 2 * vn + h * h / 8 * a1, vn + h / 2 * a2]
149
 * a4 = a[pn + h * vn + h * h / 2 * a3, vn + h * a3]
150
 *
151
 * pn+1 = pn + h * vn + h * h / 6 * [a1 + a2 + a3]
152
 * vn+1 = vn + h / 6 * (a1 + 2 * a2 + 2 * a3 + a4)
153
 *
154
 * This algorithm adaptively modifies timeStep and friction.
155
 */
156
void ForceDirectedEmbedding::embed() {
157
    if (finished)
158
        return;
159

    
160
    if (!initialized)
161
        reinitialize();
162

    
163
    long begin = clock();
164

    
165
    // modified Runge Kutta 4th order
166
    do {
167
        double hMultiplier = 0;
168
        cycle++;
169
        double nextUpdatedTimeStep = POSITIVE_INFINITY;
170

    
171
        while (true) {
172
            probCycle++;
173

    
174
            // update acceleration errors
175
            elapsedCalculationTime = ticksToMilliseconds(elapsedTicks + clock() - begin);
176
            if (elapsedCalculationTime > parameters.maxCalculationTime)
177
                break;
178

    
179
            // a1 = a[pn, vn]
180
            a(a1, pn, vn);
181

    
182
            // a2 = a[pn + h / 2 * vn + h * h / 8 * a1, vn + h / 2 * a1]
183
            addMultiplied(tpn, pn, updatedTimeStep / 2, vn);
184
            incrementWithMultiplied(tpn, updatedTimeStep * updatedTimeStep / 8, a1);
185
            addMultiplied(tvn, vn, updatedTimeStep / 2, a1);
186
            a(a2, tpn, tvn);
187

    
188
            // a3 = a[pn + h / 2 * vn + h * h / 8 * a1, vn + h / 2 * a2]
189
            addMultiplied(tpn, pn, updatedTimeStep / 2, vn);
190
            incrementWithMultiplied(tpn, updatedTimeStep * updatedTimeStep / 8, a1);
191
            addMultiplied(tvn, vn, updatedTimeStep / 2, a2);
192
            a(a3, tpn, tvn);
193

    
194
            // a4 = a[pn + h * vn + h * h / 2 * a3, vn + h * a3]
195
            addMultiplied(tpn, pn, updatedTimeStep, vn);
196
            incrementWithMultiplied(tpn, updatedTimeStep * updatedTimeStep / 2, a3);
197
            addMultiplied(tvn, vn, updatedTimeStep, a3);
198
            a(a4, tpn, tvn);
199

    
200
            // Adjust time step (h)
201
            //lastAccelerationError = maximumDifference(a1, a2, a3, a4);
202
            lastAccelerationError = averageRelativeError(a1, a2, a3, a4);
203

    
204
            if (debugLevel >= 3) {
205
                std::cout << "Prob cycle ";
206
                writeDebugInformation(std::cout);
207
            }
208

    
209
            if (lastAccelerationError == 0)
210
                break;
211

    
212
            // stop if acceleration error and time step are within range
213
            if (parameters.minAccelerationError < lastAccelerationError && lastAccelerationError < parameters.maxAccelerationError &&
214
                parameters.minTimeStep < updatedTimeStep && updatedTimeStep < parameters.maxTimeStep)
215
                break;
216

    
217
            // find now the update time step
218
            if (lastAccelerationError < parameters.maxAccelerationError) {
219
                if (hMultiplier == 0)
220
                    hMultiplier = parameters.timeStepMultiplier;
221
                else if (hMultiplier == 1.0 / parameters.timeStepMultiplier)
222
                    break;
223
            }
224
            else {
225
                if (hMultiplier == 0 || hMultiplier == parameters.timeStepMultiplier)
226
                    hMultiplier = 1.0 / parameters.timeStepMultiplier;
227
            }
228

    
229
            // stop if time step would be out of valid time step range
230
            nextUpdatedTimeStep = updatedTimeStep * hMultiplier;
231
            if (nextUpdatedTimeStep < parameters.minTimeStep || parameters.maxTimeStep < nextUpdatedTimeStep)
232
                break;
233

    
234
            updatedTimeStep = nextUpdatedTimeStep;
235
        }
236

    
237
        Assert(elapsedCalculationTime > parameters.maxCalculationTime ||
238
               lastAccelerationError <= parameters.maxAccelerationError ||
239
               nextUpdatedTimeStep <= parameters.minTimeStep);
240

    
241
        // pn+1 = pn + h * vn + h * h / 6 * [a1 + a2 + a3]
242
        add(dpn, a1, a2);
243
        increment(dpn, a3);
244
        multiply(dpn, updatedTimeStep * updatedTimeStep / 6);
245
        incrementWithMultiplied(dpn, updatedTimeStep, vn);
246
        increment(pn, dpn);
247

    
248
        // vn+1 = vn + h / 6 * (a1 + 2 * a2 + 2 * a3 + a4)
249
        addMultiplied(dvn, a1, 2, a2);
250
        incrementWithMultiplied(dvn, 2, a3);
251
        increment(dvn, a4);
252
        multiply(dvn, updatedTimeStep / 6);
253
        increment(vn, dvn);
254

    
255
        // Maximize velocity
256
        for (int i = 0; i < (int)vn.size(); i++) {
257
            if (vn[i].getLength() > parameters.maxVelocity) {
258
                vn[i].normalize();
259
                vn[i].multiply(parameters.maxVelocity);
260
            }
261
        }
262

    
263
        elapsedTime += updatedTimeStep;
264
        elapsedCalculationTime = ticksToMilliseconds(elapsedTicks + clock() - begin);
265
        if (debugLevel >= 2) {
266
            std::cout << "Cycle ";
267
            writeDebugInformation(std::cout);
268
        }
269

    
270
        // check if relaxed
271
        lastMaxVelocity = 0;
272
        lastMaxAcceleration = 0;
273
        for (int i = 0; i < (int)vn.size(); i++) {
274
            double velocity = vn[i].getLength();
275
            if (velocity > lastMaxVelocity)
276
                lastMaxVelocity = velocity;
277

    
278
            an[i].assign(a1[i]).add(a2[i]).add(a3[i]).add(a4[i]).divide(4);
279
            double acceleration = an[i].getLength();
280
            if (acceleration > lastMaxAcceleration)
281
                lastMaxAcceleration = acceleration;
282
        }
283
        double velocityRelaxFactor = std::min(1.0, parameters.velocityRelaxLimit / lastMaxVelocity);
284
        double accelerationRelaxFactor = std::min(1.0, parameters.accelerationRelaxLimit / lastMaxAcceleration);
285
        relaxFactor = std::max(relaxFactor, std::min(velocityRelaxFactor, accelerationRelaxFactor));
286

    
287
        // stop condition
288
        if (cycle > parameters.maxCycle ||
289
            elapsedCalculationTime > parameters.maxCalculationTime ||
290
            (lastMaxVelocity <= parameters.velocityRelaxLimit && lastMaxAcceleration <= parameters.accelerationRelaxLimit))
291
            finished = true;
292
    }
293
    while (!inspected && !finished);
294

    
295
    elapsedTicks += clock() - begin;
296
    elapsedCalculationTime = ticksToMilliseconds(elapsedTicks);
297

    
298
    if (debugLevel >= 1) {
299
        if (finished) {
300
            std::cout << "Runge-Kutta-4 ended ";
301
            writeDebugInformation(std::cout);
302
        }
303

    
304
        std::cout.flush();
305
    }
306
}
307

    
308
/**
309
 * an = b[pn, vn]
310
 */
311
void ForceDirectedEmbedding::a(std::vector<Pt>& an, const std::vector<Pt>& pn, const std::vector<Pt>& vn)
312
{
313
    if (debugLevel >= 4) {
314
        std::cout << "\nCalling a() with:\n";
315
        for (int i = 0; i < (int)pn.size(); i++) {
316
            std::cout << "p[" << i << "] = " << pn[i].x << ", " << pn[i].y << ", " << pn[i].z << " ";
317
            std::cout << "v[" << i << "] = " << vn[i].x << ", " << vn[i].y << ", " << vn[i].z << "\n";
318
        }
319
    }
320

    
321
    // set positions and velocities and reset forces
322
    for (int i = 0; i < (int)variables.size(); i++) {
323
        Variable *variable = variables[i];
324
        variable->assignPosition(pn[i]);
325
        variable->assignVelocity(vn[i]);
326
        variable->resetForce();
327

    
328
        if (inspected)
329
            variable->resetForces();
330
    }
331

    
332
    // calculate forces
333
    for (std::vector<IForceProvider *>::iterator it = forceProviders.begin(); it != forceProviders.end(); it++)
334
        (*it)->applyForces();
335

    
336
    // return accelerations
337
    for (int i = 0; i < (int)variables.size(); i++) {
338
        Variable *variable = variables[i];
339
        an[i].assign(variable->getAcceleration());
340
    }
341

    
342
    if (debugLevel >= 4) {
343
        std::cout << "Returning from a() with:\n";
344
        for (int i = 0; i < (int)pn.size(); i++)
345
            std::cout << "a[" << i << "] = " << an[i].x << ", " << an[i].y << ", " << an[i].z << "\n";
346
    }
347
}
348

    
349
Rc ForceDirectedEmbedding::getBoundingRectangle()
350
{
351
    double top = DBL_MAX, bottom = DBL_MIN;
352
    double left = DBL_MAX, right = DBL_MIN;
353

    
354
    const std::vector<IBody *>& bodies = getBodies();
355
    for (int i = 0; i < (int)bodies.size(); i++) {
356
        IBody *body = bodies[i];
357

    
358
        if (!dynamic_cast<WallBody *>(body)) {
359
            top = std::min(top, body->getTop());
360
            bottom = std::max(bottom, body->getBottom());
361
            left = std::min(left, body->getLeft());
362
            right = std::max(right, body->getRight());
363
        }
364
    }
365

    
366
    return Rc(left, top, 0, right - left, bottom - top);
367
}
368

    
369
void ForceDirectedEmbedding::writeDebugInformation(std::ostream& ostream)
370
{
371
    if (initialized)
372
        ostream << "at real time: " << elapsedCalculationTime << " time: " << elapsedTime << " relaxFactor: " << relaxFactor << " h: " << updatedTimeStep << " friction: " << parameters.frictionCoefficient << " min acceleration error: " << parameters.minAccelerationError << " max acceleration error: " << parameters.maxAccelerationError << " last acceleration error: " << lastAccelerationError << " cycle: " << cycle << " prob cycle: " << probCycle << " last max velocity: " << lastMaxVelocity << " last max acceleration: " << lastMaxAcceleration << "\n";
373
}