1 /**********************************************************
2 *
3 * Solve a quadratic function f(X) = X' A X + b X
4 * subject to a set of separation constraints cs
5 *
6 * Tim Dwyer, 2006
7 **********************************************************/
9 #include <math.h>
10 #include <stdlib.h>
11 #include <time.h>
12 #include <stdio.h>
13 #include <float.h>
14 #include <cassert>
15 #include <libvpsc/solve_VPSC.h>
16 #include <libvpsc/variable.h>
17 #include <libvpsc/constraint.h>
18 #include "gradient_projection.h"
19 #include <iostream>
20 #include "isnan.h"
21 #include "isinf.h"
24 using namespace std;
25 using namespace vpsc;
26 //#define CONMAJ_LOGGING 1
28 static void dumpVPSCException(char const *str, IncSolver* solver) {
29 cerr<<str<<endl;
30 unsigned m;
31 Constraint** cs = solver->getConstraints(m);
32 for(unsigned i=0;i<m;i++) {
33 cerr << *cs[i] << endl;
34 }
35 }
36 /*
37 * Use gradient-projection to solve an instance of
38 * the Variable Placement with Separation Constraints problem.
39 * Uses sparse matrix techniques to handle pairs of dummy
40 * vars.
41 */
42 unsigned GradientProjection::solve(double * b) {
43 unsigned i,j,counter;
44 if(max_iterations==0) return 0;
46 bool converged=false;
48 IncSolver* solver=NULL;
50 solver = setupVPSC();
51 //cerr << "in gradient projection: n=" << n << endl;
52 for (i=0;i<n;i++) {
53 assert(!isNaN(place[i]));
54 assert(!isinf(place[i]));
55 vars[i]->desiredPosition=place[i];
56 }
57 try {
58 solver->satisfy();
59 } catch (char const *str) {
60 dumpVPSCException(str,solver);
61 }
63 for (i=0;i<n;i++) {
64 place[i]=vars[i]->position();
65 }
66 for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
67 (*it)->updatePosition();
68 }
70 for (counter=0; counter<max_iterations&&!converged; counter++) {
71 converged=true;
72 // find steepest descent direction
73 // g = 2 ( b - Ax )
74 for (i=0; i<n; i++) {
75 old_place[i]=place[i];
76 g[i] = b[i];
77 for (j=0; j<n; j++) {
78 g[i] -= A[i][j]*place[j];
79 }
80 g[i] *= 2.0;
81 }
82 for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
83 (*it)->computeDescentVector();
84 }
85 // compute step size: alpha = ( g' g ) / ( 2 g' A g )
86 // g terms for dummy vars cancel out so don't consider
87 double numerator = 0, denominator = 0, r;
88 for (i=0; i<n; i++) {
89 numerator += g[i]*g[i];
90 r=0;
91 for (j=0; j<n; j++) {
92 r += A[i][j]*g[j];
93 }
94 denominator -= 2.0 * r*g[i];
95 }
96 double alpha = numerator/denominator;
98 // move to new unconstrained position
99 for (i=0; i<n; i++) {
100 place[i]-=alpha*g[i];
101 assert(!isNaN(place[i]));
102 assert(!isinf(place[i]));
103 vars[i]->desiredPosition=place[i];
104 }
105 for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
106 (*it)->steepestDescent(alpha);
107 }
109 //project to constraint boundary
110 try {
111 solver->satisfy();
112 } catch (char const *str) {
113 dumpVPSCException(str,solver);
114 }
115 for (i=0;i<n;i++) {
116 place[i]=vars[i]->position();
117 }
118 for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
119 (*it)->updatePosition();
120 }
121 // compute d, the vector from last pnt to projection pnt
122 for (i=0; i<n; i++) {
123 d[i]=place[i]-old_place[i];
124 }
125 // now compute beta, optimal step size from last pnt to projection pnt
126 // beta = ( g' d ) / ( 2 d' A d )
127 numerator = 0, denominator = 0;
128 for (i=0; i<n; i++) {
129 numerator += g[i] * d[i];
130 r=0;
131 for (j=0; j<n; j++) {
132 r += A[i][j] * d[j];
133 }
134 denominator += 2.0 * r * d[i];
135 }
136 for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
137 (*it)->betaCalc(numerator,denominator);
138 }
139 double beta = numerator/denominator;
141 // beta > 1.0 takes us back outside the feasible region
142 // beta < 0 clearly not useful and may happen due to numerical imp.
143 if(beta>0&&beta<1.0) {
144 for (i=0; i<n; i++) {
145 place[i]=old_place[i]+beta*d[i];
146 }
147 for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
148 (*it)->feasibleDescent(beta);
149 }
150 }
151 double test=0;
152 for (i=0; i<n; i++) {
153 test += fabs(place[i]-old_place[i]);
154 }
155 for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
156 test += (*it)->absoluteDisplacement();
157 }
158 if(test>tolerance) {
159 converged=false;
160 }
161 }
162 destroyVPSC(solver);
163 return counter;
164 }
165 // Setup an instance of the Variable Placement with Separation Constraints
166 // for one iteration.
167 // Generate transient local constraints --- such as non-overlap constraints
168 // --- that are only relevant to one iteration, and merge these with the
169 // global constraint list (including alignment constraints,
170 // dir-edge constraints, containment constraints, etc).
171 IncSolver* GradientProjection::setupVPSC() {
172 Constraint **cs;
173 //assert(lcs.size()==0);
175 for(DummyVars::iterator dit=dummy_vars.begin();
176 dit!=dummy_vars.end(); ++dit) {
177 (*dit)->setupVPSC(vars,lcs);
178 }
179 Variable** vs = new Variable*[vars.size()];
180 for(unsigned i=0;i<vars.size();i++) {
181 vs[i]=vars[i];
182 }
183 if(nonOverlapConstraints) {
184 Constraint** tmp_cs=NULL;
185 unsigned m=0;
186 if(k==HORIZONTAL) {
187 Rectangle::setXBorder(0.0001);
188 m=generateXConstraints(n,rs,vs,tmp_cs,true);
189 Rectangle::setXBorder(0);
190 } else {
191 m=generateYConstraints(n,rs,vs,tmp_cs);
192 }
193 for(unsigned i=0;i<m;i++) {
194 lcs.push_back(tmp_cs[i]);
195 }
196 }
197 cs = new Constraint*[lcs.size() + gcs.size()];
198 unsigned m = 0 ;
199 for(vector<Constraint*>::iterator ci = lcs.begin();ci!=lcs.end();++ci) {
200 cs[m++] = *ci;
201 }
202 for(vector<Constraint*>::iterator ci = gcs.begin();ci!=gcs.end();++ci) {
203 cs[m++] = *ci;
204 }
205 return new IncSolver(vars.size(),vs,m,cs);
206 }
207 void GradientProjection::clearDummyVars() {
208 for(DummyVars::iterator i=dummy_vars.begin();i!=dummy_vars.end();++i) {
209 delete *i;
210 }
211 dummy_vars.clear();
212 }
213 void GradientProjection::destroyVPSC(IncSolver *vpsc) {
214 if(acs) {
215 for(AlignmentConstraints::iterator ac=acs->begin(); ac!=acs->end();++ac) {
216 (*ac)->updatePosition();
217 }
218 }
219 unsigned m,n;
220 Constraint** cs = vpsc->getConstraints(m);
221 const Variable* const* vs = vpsc->getVariables(n);
222 delete vpsc;
223 delete [] cs;
224 delete [] vs;
225 for(vector<Constraint*>::iterator i=lcs.begin();i!=lcs.end();i++) {
226 delete *i;
227 }
228 lcs.clear();
229 //cout << " Vars count = " << vars.size() << " Dummy vars cnt=" << dummy_vars.size() << endl;
230 vars.resize(vars.size()-dummy_vars.size()*2);
231 for(DummyVars::iterator i=dummy_vars.begin();i!=dummy_vars.end();++i) {
232 DummyVarPair* p = *i;
233 delete p->left;
234 delete p->right;
235 }
236 }
238 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4 :