EconForge / interpolation.py

BSD 2-Clause "Simplified" License
123 stars 35 forks source link

Fix derivative of complete polynomial and add additional test #17

Closed cc7768 closed 7 years ago

cc7768 commented 7 years ago

Fixed the algorithm for computing the derivative of the complete polynomial when using matrices. The c1, t1, c2, t2 weren't being updated as we changed over different points. Key difference summarized by the change below:

Before

    if d == 3:
        for i1 in range(nvar):
            for i2 in range(i1, nvar):
                ix += 1
                for k in range(nobs):
                    c1, t1 = (1, 1.0) if i1==der else (0, z[i1, k])
                    c2, t2 = (c1+1, 1.0) if i2==der else (c1, z[i2, k])
                    out[ix, k] = c2*t1*t2*z[der, k]**(c2-1) if c2>0 else 0.0

                for i3 in range(i2, nvar):
                    ix += 1
                    for k in range(nobs):
                        c3, t3 = (c2+1, 1.0) if i3==der else (c2, z[i3, k])
                        out[ix, k] = c3*t1*t2*t3*z[der, k]**(c3-1) if c3>0 else 0.0

        return

After

    if d == 3:
        for i1 in range(nvar):
            for i2 in range(i1, nvar):
                ix += 1
                for k in range(nobs):
                    c1, t1 = (1, 1.0) if i1==der else (0, z[i1, k])
                    c2, t2 = (c1+1, 1.0) if i2==der else (c1, z[i2, k])
                    out[ix, k] = c2*t1*t2*z[der, k]**(c2-1) if c2>0 else 0.0

                for i3 in range(i2, nvar):
                    ix += 1
                    for k in range(nobs):
                        c1, t1 = (1, 1.0) if i1==der else (0, z[i1, k])
                        c2, t2 = (c1+1, 1.0) if i2==der else (c1, z[i2, k])
                        c3, t3 = (c2+1, 1.0) if i3==der else (c2, z[i3, k])
                        out[ix, k] = c3*t1*t2*t3*z[der, k]**(c3-1) if c3>0 else 0.0

        return