welltestdeconvolution/testFunctions.py at master · legrabow/welltestdeconvolution · 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
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
### run tests on the code

## import build-in
import numpy as np
import pandas as pd
from scipy.optimize import check_grad, approx_fprime
from datetime import datetime

## import other functions
from functions import *
from data_preparation import *
from variableProjection import *
from weightFunctions import *


def check_example():
    ## test 1: Uniform spacing of nodes in the time domain
    a = 20.4
    timeseries = np.arange(0, 20,dtype=float)
    nodes = np.log(timeseries[2:])
    nodes = np.insert(nodes, 0, np.log(0.5))
    rates = np.arange(100,119)
    z = np.full(len(nodes), np.log(a))

    convMat = generate_convMatrix(z, len(rates), nodes, timeseries)
    print(pd.DataFrame(convMat))

    v1Truncated = list()
    for t in xrange(2,len(timeseries)):
        res = (np.log(timeseries[t]) - np.log(timeseries[t - 1])) * a
        v1Truncated.append(res)

    print(v1Truncated)

    jacobian = generate_jacobianConvolution(nodes, z, rates, len(rates), timeseries)
    print(pd.DataFrame(jacobian))


def check_gradients_convMat():
    ## test 2: Check jacobian by evaluating each row of the error function based on how
    ## good its corresponding gradient approximates the numerical gradient
    timeseries = np.arange(0, 31)
    timeRange = [datetime.strptime("01/01/2007", "%d/%m/%Y"), datetime.strptime("31/01/2007","%d/%m/%Y")]
    z = np.random.normal(loc = -12, size = len(timeseries))
    rates = np.random.normal(loc = 900,scale = 100, size = len(timeseries) - 1)
    waterlevel = np.random.normal(loc = 0.6, size = len(timeseries) - 1)
    nodes = get_nodes(startNode = 0.4, timeseries = timeseries, amountNodes=None, interpolation="linear")
    wlNatIn = max(waterlevel)
    weights = dict()
    weights["rew"] = get_rateError_weight(wlNat = wlNatIn, waterlevel = waterlevel, rates = rates)
    weights["dw"] = get_derivate_weight(wlNat = wlNatIn, waterlevel = waterlevel)
    weighingFunctions = {
        "tidals":{"Consider":False, "Envelope":True},
        "holidays":{"Consider":False, "Variance":1},
        "beginning":{"Consider":True, "Sharpness":np.random.randint(4,50),"Shift":np.random.randint(0,10)},
        "pumping":{"Consider":False},
        "derivativeContrain":{"Consider":False,"unconstrainedTime":30,"relWeight":2000}
    }
    weights["totalWeightMatrix"] = get_total_weight(weighingFunctions, nodes, rates, len(timeseries), timeRange)
    yIn = rates + np.random.normal(loc = 0,scale = 10, size = len(timeseries) - 1)
    total = np.concatenate((yIn, z))
    total = np.insert(total, 0 , wlNatIn)
    total = list(total)
    rates = rates.reshape((len(rates),1))
    waterlevel = waterlevel.reshape((len(waterlevel), 1))
    for row in xrange(0, 62):
        resid = check_grad(error_measure, grad_measure, total, rates, weights, waterlevel, nodes, timeseries, row)
        eps = np.sqrt(np.finfo(float).eps)
        jacPython = approx_fprime(total, error_measure, eps, rates, weights, waterlevel, nodes, timeseries, row)
        jacMe = grad_measure(total, rates, weights, waterlevel, nodes, timeseries, row)
        print("################ row: " + str(row))
        print("Residuum: " + str(resid))
        print("Numeric gradient:")
        print(jacPython)
        print("Program's gradient:")
        print(jacMe)


def grad_measure(total, rates, weights, waterlevel, nodes, timeseries, row):
    ## Gets the jacobian gradient of a given row
    x = total[:(len(rates) + 1)]
    x = np.array(x)
    x = x.reshape((len(x), 1))
    z = total[(len(rates) + 1):]
    z = np.array(z)
    z = z.reshape((len(z), 1))
    fMat = generate_fMatrix(weights["rew"], z, len(rates), len(waterlevel), nodes, timeseries, weights["totalWeightMatrix"])
    zJacobian = generate_jacobian(nodes, z, x[1:], weights["dw"], len(rates), timeseries, weights["totalWeightMatrix"])
    totalJacobian = np.hstack((fMat, zJacobian))
    return list(totalJacobian[row,:])

def error_measure(total, rates, weights, waterlevel, nodes, timeseries, row):
    ## Gets the error measure of a given row
    x = total[:(len(rates) + 1)]
    x = np.array(x)
    x = x.reshape((len(x), 1))
    z = total[(len(rates) + 1):]
    z = np.array(z)
    z = z.reshape((len(z), 1))
    fMat = generate_fMatrix(weights["rew"], z, len(rates), len(waterlevel), nodes, timeseries, weights["totalWeightMatrix"])
    vVec = generate_vVector(weights["rew"], weights["dw"], z, waterlevel, rates, nodes, weights["totalWeightMatrix"])
    error = fMat.dot(x) - vVec
    return error[row][0]


def check_gradients_convMat_NORATES():
    ## test 3: Same as test2, but check for the case of no rate-errors
    timeseries = np.arange(0, 31)
    timeRange = [datetime.strptime("01/01/2007", "%d/%m/%Y"), datetime.strptime("31/01/2007","%d/%m/%Y")]
    z = np.random.normal(loc = -12, size = len(timeseries))
    rates = np.random.normal(loc = 900,scale = 100, size = len(timeseries) - 1)
    waterlevel = np.random.normal(loc = 0.6, size = len(timeseries) - 1)
    nodes = get_nodes(startNode = 0.4, timeseries = timeseries, amountNodes=None, interpolation="linear")
    wlNatIn = max(waterlevel)
    weights = dict()
    weights["rew"] = 0
    weights["dw"] = get_derivate_weight(wlNat = wlNatIn, waterlevel = waterlevel)
    weighingFunctions = {
        "tidals":{"Consider":False, "Envelope":True},
        "holidays":{"Consider":False, "Variance":1},
        "beginning":{"Consider":True, "Sharpness":np.random.randint(4,50),"Shift":np.random.randint(0,10)},
        "pumping":{"Consider":False},
        "derivativeContrain":{"Consider":False,"unconstrainedTime":30,"relWeight":2000}
    }
    weights["totalWeightMatrix"] = get_total_weight(weighingFunctions, nodes, rates, len(timeseries), timeRange)
    total = np.insert(z, 0 , wlNatIn)
    total = list(total)
    rates = rates.reshape((len(rates),1))
    waterlevel = waterlevel.reshape((len(waterlevel), 1))
    for row in xrange(0, 62):
        resid = check_grad(error_measure_NORATES, grad_measure_NORATES, total, rates, weights, waterlevel, nodes, timeseries, row)
        eps = np.sqrt(np.finfo(float).eps)
        jacPython = approx_fprime(total, error_measure_NORATES, eps, rates, weights, waterlevel, nodes, timeseries, row)
        jacMe = grad_measure_NORATES(total, rates, weights, waterlevel, nodes, timeseries, row)
        print("################ row: " + str(row))
        print("Residuum: " + str(resid))
        print("Numeric gradient:")
        print(jacPython)
        print("Program's gradient:")
        print(jacMe)


def grad_measure_NORATES(total, rates, weights, waterlevel, nodes, timeseries, row):
    z = total[1:]
    z = np.array(z)
    z = z.reshape((len(z), 1))
    fMat = generate_fMatrix(weights["rew"], z, len(rates), len(waterlevel), nodes, timeseries, weights["totalWeightMatrix"])
    A = fMat[:,0]
    A = A.reshape((len(A), 1))
    zJacobian = generate_jacobian(nodes, z, rates, weights["dw"], len(rates), timeseries, weights["totalWeightMatrix"])
    totalJacobian = np.hstack((A, zJacobian))
    return list(totalJacobian[row,:])

def error_measure_NORATES(total, rates, weights, waterlevel, nodes, timeseries, row):
    x = np.vstack([total[0], rates])
    z = total[1:]
    z = np.array(z)
    z = z.reshape((len(z), 1))
    fMat = generate_fMatrix(weights["rew"], z, len(rates), len(waterlevel), nodes, timeseries, weights["totalWeightMatrix"])
    vVec = generate_vVector(weights["rew"], weights["dw"], z, waterlevel, rates, nodes, weights["totalWeightMatrix"])
    error = fMat.dot(x) - vVec
    return error[row][0]