Mathematical-Optimization/ScaledDifferentiableFunction.py at main · laschoe/Mathematical-Optimization · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
import numpy as np
import math
from typing import Union, Callable, Optional
from Set import MultidimensionalInterval
from DifferentiableFunction import DifferentiableFunction
from BayesianOptimization import BO

class ScaledDifferentiableFunction(object):
    def __init__(self):
        super().__init__()

    @classmethod
    def getScalingParamsDim(cls, input_scalar: Union[int, float, np.array], output_scalar: Union[int, float, np.array], input_offset: Union[int, float, np.array], output_offset: Union[int, float, np.array], input_dim: int, output_dim: int) -> tuple:
        # Convert scalars to arrays if necessary
        if isinstance(input_scalar, (int, float)):
            input_scalar = np.array([input_scalar])
        if isinstance(input_offset, (int, float)):
            input_offset = np.array([input_offset])
        if isinstance(output_scalar, (int, float)):
            output_scalar = np.array([output_scalar])
        if isinstance(output_offset, (int, float)):
            output_offset = np.array([output_offset])
        # print shapes
        # print("Shapes: ", input_scalar.shape, input_offset.shape, output_scalar.shape, output_offset.shape)
        # Adjust shapes if necessary
        if input_scalar.shape == (1,):
            input_scalar = np.repeat(input_scalar, input_dim)
        if output_scalar.shape == (1,):
            output_scalar = np.repeat(output_scalar, output_dim)
        if input_offset.shape == (1,):
            input_offset = np.repeat(input_offset, input_dim)
        if output_offset.shape == (1,):
            output_offset = np.repeat(output_offset, output_dim)

        # print("Shapes: ", input_scalar.shape, input_offset.shape, output_scalar.shape, output_offset.shape)

        # Convert to diagonal matrices
        input_scalar_matrix = np.diag(input_scalar)
        output_scalar_matrix = np.diag(output_scalar)
        # print("Shapes: ", input_scalar_matrix.shape, input_offset.shape, output_scalar_matrix.shape, output_offset.shape)

        assert input_scalar_matrix.shape == (input_dim, input_dim)
        assert input_offset.shape == (input_dim,)
        assert output_scalar_matrix.shape == (output_dim, output_dim)
        assert output_offset.shape == (output_dim,)

        return input_scalar_matrix, output_scalar_matrix, input_offset, output_offset, input_scalar, output_scalar

    @classmethod
    def getScaledFunction(cls, f: DifferentiableFunction, input_scalar: Optional[Union[int, float, np.array]] = None, output_scalar: Optional[Union[int, float, np.array]] = None, input_offset: Optional[Union[int, float, np.array]] = None, output_offset: Optional[Union[int, float, np.array]] = None, auto_scale: Optional[bool] = False) -> DifferentiableFunction:
        """
        Returns a function that scales both input and output and can add offsets to input and output:
        x -> output_scalar * f(input_scalar * x + input_offset) + output_offset
        """
        output_dim = f.evaluate(f.domain.point()).shape[0]
        input_dim = f.domain._ambient_dimension

        if input_scalar is None:
            input_scalar = 1
        if output_scalar is None:
            output_scalar = 1
        if input_offset is None:
            input_offset = 0
        if output_offset is None:
            output_offset = 0

        input_scalar_matrix, output_scalar_matrix, input_offset, output_offset, input_scalar, output_scalar = cls.getScalingParamsDim(input_scalar, output_scalar, input_offset, output_offset, input_dim, output_dim)

        if auto_scale:
            # MulidimensionalInterval [-1,1] as domain
            domain = MultidimensionalInterval(np.array([-1]*input_dim), np.array([1]*input_dim))
        else:
            domain = f.domain

        return DifferentiableFunction(
            name=f"{output_scalar} * {f.name}({input_scalar} * x + {input_offset}) + {output_offset}",
            domain=domain,
            evaluate=lambda x: np.matmul(output_scalar_matrix, f.evaluate(np.matmul(input_scalar_matrix, x) + input_offset)) + output_offset if isinstance(x, np.ndarray) else output_scalar * f.evaluate(input_scalar * x + input_offset) + output_offset,
            jacobian=lambda x: np.matmul(output_scalar_matrix, np.matmul(f.jacobian(np.matmul(input_scalar_matrix, x) + input_offset), input_scalar_matrix)) if isinstance(x, np.ndarray) else output_scalar * f.jacobian(input_scalar * x + input_offset) * input_scalar
        )

    @classmethod
    def getAutoScaledFunction(cls, f: DifferentiableFunction, samples: Optional[int] = 1000, BO_option: Optional[bool] = False):
        """
        Returns the autoscaled function of the given function f, so that the output values are in [-1, 1] and fixed for inputs in [-1,1].
        Only implemented for MultidimensionalInterval domains.
        :param f: The function to be scaled.
        :param samples: The number of samples to be used for min max calculation if BO_option is False.
        :param BO_option: If True, the min max values are calculated via Bayesian Optimization, otherwise via sampling.
        :return: The autoscaled function.
        """
        # assert domain is MultiDimensionalInterval
        domain = f.domain
        assert isinstance(domain, MultidimensionalInterval), "The domain of the function must be a MultidimensionalInterval."
        # assert single value output for possible min max calculation
        assert (f.evaluate(f.domain.point())).shape == (1,), "The function should give a single value output"

        # BO option
        BO_option = BO_option
        if BO_option:
            # get Min and Max of the function via Bayesian Optimization
            bo = BO()
            # single value output
            min_point = round(float(f.evaluate(bo.Minimize(f))[0]),6)
            max_point = round(float(f.evaluate(bo.Minimize((-1)*f))[0]),6)
        else:
            # get Min and Max of the function via Sampling
            samples = samples
            x = np.array([domain.point() for i in range(samples)])
            y = np.array([f.evaluate(x[i]) for i in range(samples)])
            min_point = round(float(np.min(y)),6)
            max_point = round(float(np.max(y)),6)

        # print("Min: ", min_point, "Max: ", max_point)

        # calculate the output scalar and offset so that the output is scaled to [-1, 1]
        output_scalar = 2 / (max_point - min_point)
        output_offset = -1 - (2 * min_point / (max_point - min_point))

        # calculate the input scalar and offset so that the input is scaled to [-1, 1]
        input_scalar = 1/(2 / (domain._upper_bounds - domain._lower_bounds))
        input_offset = 1 * (domain._upper_bounds + domain._lower_bounds) / (domain._upper_bounds - domain._lower_bounds)

        # print("Params: ", output_scalar, output_offset, input_scalar, input_offset)

        return cls.getScaledFunction(f, input_scalar=input_scalar, input_offset=input_offset, output_scalar=output_scalar, output_offset=output_offset, auto_scale=True)