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