source: trunk/ADOL-C/examples/additional_examples/hessmat/hessmat.cpp

Last change on this file was 387, checked in by kulshres, 7 years ago

switch extra cast depending on advanced branching or not , else the example won't compile without the advanced branching

  • Property svn:keywords set to Author Date Id Revision
File size: 8.0 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     hessmat.cpp
4 Revision: $Id: hessmat.cpp 387 2013-01-25 15:47:43Z stefan $
5 Contents: example for testing the routines:
6           hov_wk_forward    ( = Higher Order Vector forward With Keep )
7           hos_ov_reverse    ( = Higher Order Scalar reverse over vectors)
8
9 Copyright (c) Andrea Walther, Andreas Kowarz, Olaf Vogel
10 
11 This file is part of ADOL-C. This software is provided as open source.
12 Any use, reproduction, or distribution of the software constitutes
13 recipient's acceptance of the terms of the accompanying license file.
14 
15---------------------------------------------------------------------------*/
16
17/****************************************************************************/
18/*                                                                 INCLUDES */
19#include <adolc/adolc.h>
20
21#include <stdlib.h>
22#include <iostream>
23using namespace std;
24
25/****************************************************************************/
26/*                                                                     MAIN */
27int main() {
28    int i,j,l,m,n,d,q,bd, keep;
29
30    /*--------------------------------------------------------------------------*/
31    /* inputs */
32    cout << "vector x Hessian x matrix for the function \n\n";
33    cout << " y[0] = cos(x[0])* ...*cos(x[n]) \n";
34    cout << " y[1] = x[0]^n \n";
35    cout << " y[2] = condassign(y[i],y[0]>y[1],y[1],y[0]) \n";
36    cout << " y[3] = sin(x[0])+ ...+sin(x[n]) \n";
37    cout << " y[4] = exp(x[0])- ...-exp(x[n]) \n";
38    cout << " y[5] = pow(y[1],3) \n";
39    cout << " y[6] += y[5]*y[4] \n";
40    cout << " y[7] -= y[6]*y[5] \n";
41    cout << " y[j] = 1/x[0]/ .../x[n], j > 3 \n\n";
42
43    cout << " Number of independents = ?\n ";
44    cin >> n;
45    cout << " Number of dependents =  ?\n ";
46    cin >> m;
47    cout << " Degree d (for forward) =  ?\n";
48    cin >> d;
49    cout << " keep (degree of corresponding reverse = keep-1) =  ?\n";
50    cout << "       keep <= d+1 must be valid \n";
51    cin >> keep;
52    cout << " Number of directions =  ?\n ";
53    cin >> q;
54
55    /*--------------------------------------------------------------------------*/
56    /* allocations and inits */
57
58    double* xp = new double[n];                      /* passive indeps        */
59    double* yp = new double[m];                      /* passive depends       */
60
61    /* vector x Hessian x matrix = Upp x H x XPPP */
62
63    double* Up = myalloc(m);        /* vector on left-hand side                */
64    double** Upp = myalloc(m,d+1);   /* vector on left-hand side                */
65    double*** Xppp = myalloc(n,q,d); /* matrix on right-hand side             */
66    double*** Zppp = myalloc(q,n,d+1); /* result of Up x H x XPPP             */
67
68    double*** Yppp = myalloc(m,q,d); /* results of needed hos_wk_forward      */
69
70    /* check results with usual lagra-Hess-vec  */
71
72    double** Xpp = myalloc(n,d);
73    double** V   = myalloc(n,q);
74    double** W   = myalloc(q,n);
75    double** H   = myalloc(n,n);
76    double** Ypp = myalloc(m,d);
77    double** Zpp = myalloc(n,d+1);
78    /* inits */
79
80    for (l=0; l<d; l++)                    /* first everything is set to zero */
81        for (i=0; i<n; i++)
82            for (j=0;j<q;j++)
83                Xppp[i][j][l] = 0;
84
85    /* now carthesian directions as choosen as    */
86    /* matrix on right-hand side of Up x H x XPPP */
87    bd = (n<q)?n:q;
88    for (j=0;j<bd;j++)
89        Xppp[j][j][0] = 1;
90
91    for (i=0; i<m; i++)         /* vector on left-hand side of Up x H x XPPP  */
92    {Up[i] = 1;                /* is initialised with 1's                    */
93        Upp[i][0] = 1;
94        for (j=1;j<=d;j++)
95            Upp[i][j] = 0;
96    }
97
98    for (i=0; i<n; i++)                    /* first everything is set to zero */
99        for (j=0;j<d;j++)
100            Xpp[i][j] = 0;
101    Xpp[0][0] = 1;                 /* now one carthesian direction as choosen */
102    /* as vector for lagra-Hess-vec            */
103
104    for (i=0; i<n; i++)                            /* inits of passive indeps */
105        xp[i] = (i+1.0)/(2.0+i);
106
107    for (i=0; i<n; i++) {
108        for (j=0;j<q;j++)
109            V[i][j] = 0;
110        if (i < q)
111            V[i][i] = 1;
112    }
113
114    /*--------------------------------------------------------------------------*/
115    trace_on(1);                                      /* tracing the function */
116
117    adouble* x = new adouble[n];                     /* active indeps         */
118    adouble* y = new adouble[m];                     /* active depends        */
119    for(i=0;i<m;i++)
120        y[i] = 1;
121    for (i=0; i<n; i++) {
122        x[i] <<= xp[i];
123        y[0] *= cos(x[i]);
124    }
125
126    for(i=1;i<m;i++)
127        for(j=0;j<n;j++) {
128            switch (i) {
129                case 1 :
130                    y[i] *= x[0];
131                    break;
132                case 2 :
133#ifndef ADOLC_ADVANCED_BRANCHING
134                    condassign(y[i],adouble(y[0]>y[1]),y[1],y[0]);
135#else
136                    condassign(y[i],(y[0]>y[1]),y[1],y[0]);
137#endif
138                    break;
139                case 3 :
140                    y[i] -= sin(x[j]);
141                    break;
142                case 4 :
143                    y[i] -= exp(x[j]);
144                    break;
145                case 5 :
146                    y[5] = pow(y[1],3);
147                case 6 :
148                    y[6] += y[5]*y[4];
149                case 7 :
150                    y[7] -= y[6]*y[5];
151                default :
152                    y[i] /= x[j];
153            }
154        }
155    for (i=0; i<m; i++)
156        y[i] >>= yp[i] ;
157    trace_off();
158
159
160    /*--------------------------------------------------------------------------*/
161    /* work on the tape */
162
163
164    /* compute results of lagra_hess_vec */
165    /* the following is equal to calls inside of lagra_hess_vec(..) */
166    /* direct calls to the basic routines hos_forward and hos_reverse */
167    /* seem to be faster than call of lagra_hess_vec(..) */
168    /* at least in some of our test cases */
169
170    hos_forward(1,m,n,d,keep,xp,Xpp,yp,Ypp);
171    hos_reverse(1,m,n,keep-1,Up,Zpp);
172
173    printf("\n Results of hos_reverse:\n\n");
174
175    for (i=0; i<=d; i++) {
176        printf(" d = %d \n",i);
177        for (j=0;j<n;j++)
178            printf(" %6.3f ",Zpp[j][i]);
179        printf("\n");
180    }
181
182    /* The new drivers. First, hov_wk_forward(..) is called.
183       So far, it was impossible to store the results of
184       a higher-order-vector (=hov) forward in order to perform
185       a corresponding reverse sweep (for no particular reason.
186       Now we have hov with keep (=wk) and the results needed on
187       the way back are stored in a specific tape */
188
189    hov_wk_forward(1,m,n,d,keep,q,xp,Xppp,yp,Yppp);
190
191    /* The corresponding reverse sweep
192       So far we had only a higher-order-scalar (=hos, scalar because
193       only one vector on the left-hand-side) for a scalar forward
194       call.
195       Now, we use the stored vector information (= hos vector)
196       to compute multiple lagra_hess_vec at once */
197
198    hos_ov_reverse(1,m,n,keep-1,q,Upp,Zppp);
199
200    printf("\n Results of hosv_reverse:\n");
201
202    for (l=0; l<q; l++) {
203        for (i=0; i<=d; i++) {
204            printf(" d = %d \n",i);
205            for (j=0;j<n;j++)
206                printf(" %6.3f ",Zppp[l][j][i]);
207            printf("\n");
208        }
209        printf("\n\n");
210    }
211
212    if (m==1) {
213        printf("hess_mat:\n");
214        hess_mat(1,n,q,xp,V,W);
215        for (i=0; i<q; i++) {
216            for (j=0;j<n;j++)
217                printf(" %6.3f ",W[i][j]);
218            printf("\n");
219        }
220        printf("hessian2:\n");
221        hessian2(1,n,xp,H);
222        for (i=0; i<n; i++) {
223            for (j=0;j<n;j++)
224                printf(" %6.3f ",H[i][j]);
225            printf("\n");
226        }
227    }
228
229    myfree(Zpp);
230    myfree(Ypp);
231    myfree(H);
232    myfree(W);
233    myfree(V);
234    myfree(Xpp);
235    myfree(Zppp);
236    myfree(Yppp);
237    myfree(Xppp);
238    myfree(yp);
239    myfree(xp);
240    myfree(Up);
241    return 1;
242}
243
244
245/****************************************************************************/
246/*                                                               THAT'S ALL */
247
248
249
250
Note: See TracBrowser for help on using the repository browser.