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