SC-7-1-FinalVersion/pfit.py at master · orebas/SC-7-1-FinalVersion · 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
import numpy
import ptest

def pFit(xsubk, fsubk):
    n =len(xsubk)
    vandermonde = numpy.zeros((n,n))
    for i in range(n):
        for k  in range(n):
            vandermonde[i][k] = xsubk[i] ** k
    asubk = numpy.linalg.solve(vandermonde,fsubk)
    cond = numpy.linalg.cond(vandermonde)
    return (asubk,cond)


#below fits a cubic spline to input data.  We assume the xs are sorted in increasing order.
# we return them as part of our data as we expect them to be passed around as part of the
#spline definition.
#some notes to make it readable:
# the polynomials are sum of a(i,k)(x-x_k)^i for i=0 to 3
# a(i,k) is stored in coeffs[4*k+i]


def splinefit(xsubk,fsubk):
    A = []
    B = []
    n = len(xsubk)
    coeffs = [0.0] * (4 * (n-1))  #there's a cubic poly for each interval, i.e. for all but the last x
    for k in range(0, n-1):
        row = [0.0] * (4 * (n-1))
        row[4*k+0] = 1
        A.append(row)
        B.append(fsubk[k])
    for k in range(0, n-1):
        t = xsubk[k+1] - xsubk[k]
        row = [0.0] * (4 * (n - 1))
        row[4 * k + 0] = 1
        row[4 * k + 1] = t
        row[4 * k + 2] = t**2
        row[4 * k + 3] = t**3
        A.append(row)
        B.append(fsubk[k+1])
    for k in range(0,n-2):
        t = xsubk[k + 1] - xsubk[k]
        row = [0.0] * (4 * (n - 1))
        row[4*k+1] = 1.0
        row[4*k+2] = 2.0*t
        row[4 * k + 3] = 3.0 * t * t
        row[4 * (k+1) + 1] = -1
        A.append(row)
        B.append(0)
    for k in range(0,n-2):
        t = xsubk[k + 1] - xsubk[k]
        row = [0.0] * (4 * (n - 1))
        row[4*k+2] = 2.0
        row[4*k+3] = 6.0*t
        row[4 * (k+1) + 2] = -2
        A.append(row)
        B.append(0)

    t = xsubk[k + 1] - xsubk[k]
    row = [0.0] * (4 * (n - 1))
    row[4 * 0 + 2] = 2.0
    A.append(row)
    B.append(0)

    t = xsubk[-1] - xsubk[-2]
    row = [0.0] * (4 * (n - 1))
    row[4 * (n-2) + 2] = 2.0
    row[4 * (n - 2) + 3] = 6.0 * t
    A.append(row)
    B.append(0)
    A = numpy.array(A)
    B = numpy.array(B)
    #print(A)
    #print(B)

    coeffs = numpy.linalg.solve(A,B)
    cond = numpy.linalg.cond(A)
    #print (coeffs)
    return (coeffs,cond)

def rfit(xs,fs,l):
    n = len(xs)
    A = []
    B = []
    for i in range(n):
        row = [0] * n
        for k in range(n):
            row[k] = ptest.phi(xs[i] - xs[k], l)
        A.append(row)
        B.append(fs[i])
    w = numpy.linalg.solve(A,B)
    return(w, numpy.linalg.cond(A))

#dummy comment