boundary_transformation.c 5.84 KB
Newer Older
1
2
3
4
5
6
7
8
#include <stddef.h>
#include <stdlib.h>
#include <math.h>
#include <float.h> /* DBL_MAX */
#include <stdio.h>
#include "boundary_transformation.h"

static int do_assertions = 1;
9
10
static unsigned long _index(boundary_transformation_t* t, unsigned long i);
static void _FatalError(const char* s);
11
12
13
static double default_lower[1];
static double default_upper[1];

14
15
16
17
void boundary_transformation_init(boundary_transformation_t* t,
    double const* lower_bounds,
    double const* upper_bounds,
    unsigned long len_of_bounds)
18
{
19
20
21
22
23
24
25
26
27
28
29
30
    unsigned i;
    double const* ub, *lb;
    default_lower[0] = DBL_MAX / -1e2;
    default_upper[0] = DBL_MAX / 1e2;
    t->lower_bounds = lower_bounds;
    t->upper_bounds = upper_bounds;
    t->len_of_bounds = len_of_bounds;

    if (lower_bounds == NULL && len_of_bounds <= 1) /* convenience default */
        t->lower_bounds = default_lower;
    if (upper_bounds == NULL && len_of_bounds <= 1)
        t->upper_bounds = default_upper;
31
    if (len_of_bounds == 0) {
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
        t->lower_bounds = default_lower;
        t->upper_bounds = default_upper;
        t->len_of_bounds = 1;
    }

    if (t->lower_bounds == NULL || t->upper_bounds == NULL)
        _FatalError("init: input upper_bounds or lower_bounds was NULL and "
                    "len_of_bounds > 1");

    /* compute boundaries in pre-image space, al and au */
    t->al = calloc(t->len_of_bounds, sizeof(double));
    t->au = calloc(t->len_of_bounds, sizeof(double));
    if (!t->al || !t->au)
        _FatalError(" in _init(): could not allocate memory");

    lb = t->lower_bounds;
    ub = t->upper_bounds;
49
    for (i = 0; i < t->len_of_bounds; ++i) {
50
51
52
53
54
55
56
        if (lb[i] == ub[i])
            _FatalError("in _init: lower and upper bounds must be different in all "
                        "variables");
        /* between lb+al and ub-au transformation is the identity */
        t->al[i] = fmin((ub[i] - lb[i]) / 2., (1. + fabs(lb[i])) / 20.);
        t->au[i] = fmin((ub[i] - lb[i]) / 2., (1. + fabs(ub[i])) / 20.);
    }
57
}
58
void boundary_transformation_exit(boundary_transformation_t* t)
59
{
60
61
62
63
    if (t->al)
        free(t->al);
    if (t->au)
        free(t->au);
64
65
}

66
67
void boundary_transformation(boundary_transformation_t* t, double const* x,
    double* y, unsigned long len)
68
{
69
70
71
    double lb, ub, al, au;
    unsigned long i;
    boundary_transformation_shift_into_feasible_preimage(t, x, y, len);
72
    for (i = 0; i < len; ++i) {
73
74
75
76
        lb = t->lower_bounds[_index(t, i)];
        ub = t->upper_bounds[_index(t, i)];
        al = t->al[_index(t, i)];
        au = t->au[_index(t, i)];
77
78
79
80
        if (y[i] < lb + al)
            y[i] = lb + (y[i] - (lb - al)) * (y[i] - (lb - al)) / 4. / al;
        else if (y[i] > ub - au)
            y[i] = ub - (y[i] - (ub + au)) * (y[i] - (ub + au)) / 4. / au;
81
    }
82
}
83
84
85
void boundary_transformation_discrete(boundary_transformation_t* t,
    double const* x, double* y,
    unsigned long len)
86
87
{

88
89
    double lb, ub, al, au;
    unsigned long i;
90

91
    // boundary_transformation_shift_into_feasible_preimage(t, x, y, len);
92
    for (i = 0; i < len; ++i) {
93
94
95
96
        lb = t->lower_bounds[_index(t, i)];
        ub = t->upper_bounds[_index(t, i)];
        al = t->al[_index(t, i)];
        au = t->au[_index(t, i)];
97

98
        y[i] = x[i];
99

100
101
102
103
104
        if (y[i] < lb)
            y[i] = lb;
        else if (y[i] > ub)
            y[i] = ub;
    }
105
106
}
void boundary_transformation_shift_into_feasible_preimage(
107
108
    boundary_transformation_t* t, double const* x, double* y,
    unsigned long len)
109
{
110
111
112
    double lb, ub, al, au, r, xlow, xup;
    unsigned long i;

113
    for (i = 0; i < len; ++i) {
114
115
116
117
118
119
        lb = t->lower_bounds[_index(t, i)];
        ub = t->upper_bounds[_index(t, i)];
        al = t->al[_index(t, i)];
        au = t->al[_index(t, i)];
        xlow = lb - 2 * al - (ub - lb) / 2.0;
        xup = ub + 2 * au + (ub - lb) / 2.0;
120
121
122
123
        r = 2 * (ub - lb + al + au); /* == xup - xlow == period of the transformation */

        y[i] = x[i];

124
        if (y[i] < xlow) { /* shift up */
125
126
            y[i] += r * (1 + (int)((xlow - y[i]) / r));
        }
127
        if (y[i] > xup) { /* shift down */
128
129
130
131
132
            y[i] -= r * (1 + (int)((y[i] - xup) / r));
            /* printf(" \n%f\n", fmod(y[i] - ub - au, r)); */
        }
        if (y[i] < lb - al) /* mirror */
            y[i] += 2 * (lb - al - y[i]);
133
        if (y[i] > ub + au)
134
135
            y[i] -= 2 * (y[i] - ub - au);

136
        if ((y[i] < lb - al - 1e-15) || (y[i] > ub + au + 1e-15)) {
137
138
139
140
141
142
            printf("BUG in boundary_transformation_shift_into_feasible_preimage: "
                   "lb=%f, ub=%f, al=%f au=%f, y=%f\n",
                lb, ub, al, au, y[i]);
            _FatalError("BUG");
        }
    }
143
}
144
145
146
void boundary_transformation_inverse(boundary_transformation_t* t,
    double const* x, double* y,
    unsigned long len)
147
{
148
149
150
    double lb, ub, al, au;
    unsigned long i;

151
    for (i = 0; i < len; ++i) {
152
153
154
155
156
157
158
159
160
161
        lb = t->lower_bounds[_index(t, i)];
        ub = t->upper_bounds[_index(t, i)];
        al = t->al[_index(t, i)];
        au = t->al[_index(t, i)];
        y[i] = x[i];
        if (y[i] < lb + al)
            y[i] = (lb - al) + 2 * sqrt(al * (y[i] - lb));
        else if (y[i] > ub - au)
            y[i] = (ub + au) - 2 * sqrt(au * (ub - y[i]));
    }
162
    if (11 < 3 || do_assertions) {
163
164
165
166
167
168
169
170
171
172
173
174
        double* z = calloc(len, sizeof(double));
        for (i = 0; i < len; ++i)
            z[i] = y[i];
        boundary_transformation(t, z, y, len);
        for (i = 0; i < len; ++i)
            if (fabs(y[i] - x[i]) > 1e-14)
                printf("  difference for index %ld should be zero, is %f ", i,
                    y[i] - x[i]);
        for (i = 0; i < len; ++i)
            y[i] = z[i];
        free(z);
    }
175
}
176
static unsigned long _index(boundary_transformation_t* t, unsigned long i)
177
{
178
    return i < t->len_of_bounds ? i : t->len_of_bounds - 1;
179
}
180
static void _FatalError(const char* s)
181
{
182
183
    printf("Fatal error in boundary_transformation: %s\n", s);
    exit(1);
184
}