diff options
| -rw-r--r-- | python/optimization.py | 39 |
1 files changed, 39 insertions, 0 deletions
diff --git a/python/optimization.py b/python/optimization.py index 2651760e..15fd0ec8 100644 --- a/python/optimization.py +++ b/python/optimization.py @@ -1,5 +1,6 @@ import numpy as np import numpy.linalg as la +from scipy.interpolate import PchipInterpolator def decr_fun(a, b): return np.sqrt(np.square(a).sum() + np.square(b).sum()) @@ -56,6 +57,44 @@ def KLfit(P, q, b): return {"obj": sum(x*(np.log(x) - np.log(q))), "weight":x, "status": decr} +def interpweights(w, v1, v2): + """ Given L=(w, v1), computes neww such that newL=(neww, v2) = L in distribution""" + cumw = np.cumsum(w) + neww = PchipInterpolator(v1, cumw)(v2, der = 1) + interpweights = new/np.sum(neww) + return interpweights + +def interpvalues(w, v, neww): + """ Given a distribution D=(w,v), compute new weights + such that Dnew=(neww, newv) equals D in distribution """ + cumw = np.cumsum(w) + cdf = PchipInterpolator(v, cumw) + eps = 1e-4 + neww = np.zeros(neww.size) + cumneww = np.cumsum(neww) + mid = 0 + for i in range(neww.size): + iter = 0 + hi = cdf(1) + lo = mid + if hi < cumneww[i]: + neww[i] = hi + continue + if cdf(lo) > cumneww[i]: + neww[i] = lo + continue + mid = (lo+hi)/2 + iter = 0 + while( abs(cdf(mid)-cumneww[i]) > eps ): + if cdf(mid) > cumneww[i]: + hi = mid + else: + lo = mid + mid = (lo + hi)/2 + newv[i] = mid + return newv + + if __name__=="__main__": #write some small test P = np.array([[5, 4, 3, 2, 1],[3, 3, 4, 3, 3]]) |