(Dies zu spät für Ihre Hausaufgaben, aber hoffe, es hilft trotzdem.)
Zunächst ist die dynamische Programmierung in Python/numpy für k = 4
nur, Sie verstehen, wie die dynamische Programmierung arbeitet zu helfen; Sobald Sie das verstehen, sollte das Schreiben einer Schleife für jedes k einfach sein.
Auch Cost[]
ist eine 2D-Matrix, Raum O (n^2); die Notizen am Ende sehen für den Weltraum O immer nach unten (n k)
#!/usr/bin/env python
""" split4.py: min-cost split into 4 pieces, dynamic programming k=4 """
from __future__ import division
import numpy as np
__version__ = "2014-03-09 mar denis"
#...............................................................................
def split4(Cost, verbose=1):
""" split4.py: min-cost split into 4 pieces, dynamic programming k=4
min Cost[0:a] + Cost[a:b] + Cost[b:c] + Cost[c:n]
Cost[a,b] = error in least-squares line fit to xy[a] .. xy[b] *including b*
or error in lsq horizontal lines, sum (y_j - av y) ^2 for each piece --
o--
o-
o---
o----
| | | |
0 2 5 9
(Why 4 ? to walk through step by step, then put in a loop)
"""
# speedup: maxlen 2 n/k or so
Cost = np.asanyarray(Cost)
n = Cost.shape[1]
# C2 C3 ... costs, J2 J3 ... indices of best splits
J2 = - np.ones(n, dtype=int) # -1, NaN mark undefined/bug
C2 = np.ones(n) * np.NaN
J3 = - np.ones(n, dtype=int)
C3 = np.ones(n) * np.NaN
# best 2-splits of the left 2 3 4 ...
for nleft in range(1, n):
J2[nleft] = j = np.argmin([ Cost[0,j-1] + Cost[j,nleft] for j in range(1, nleft+1)]) + 1
C2[nleft] = Cost[0,j-1] + Cost[j,nleft]
# an idiom for argmin j, min value c together
# best 3-splits of the left 3 4 5 ...
for nleft in range(2, n):
J3[nleft] = j = np.argmin([ C2[j-1] + Cost[j,nleft] for j in range(2, nleft+1)]) + 2
C3[nleft] = C2[j-1] + Cost[j,nleft]
# best 4-split of all n --
j4 = np.argmin([ C3[j-1] + Cost[j,n-1] for j in range(3, n)]) + 3
c4 = C3[j4-1] + Cost[j4,n-1]
j3 = J3[j4]
j2 = J2[j3]
jsplit = np.array([ 0, j2, j3, j4, n ])
if verbose:
print "split4: len %s pos %s cost %.3g" % (np.diff(jsplit), jsplit, c4)
print "split4: J2 %s C2 %s" %(J2, C2)
print "split4: J3 %s C3 %s" %(J3, C3)
return jsplit
#...............................................................................
if __name__ == "__main__":
import random
import sys
import spread
n = 10
ncycle = 2
plot = 0
seed = 0
# run this.py a=1 b=None c=[3] 'd = expr' ... in sh or ipython
for arg in sys.argv[1:]:
exec(arg)
np.set_printoptions(1, threshold=100, edgeitems=10, linewidth=100, suppress=True)
np.random.seed(seed)
random.seed(seed)
print "\n", 80 * "-"
title = "Dynamic programming least-square horizontal lines %s n %d seed %d" % (
__file__, n, seed)
print title
x = np.arange(n + 0.)
y = np.sin(2*np.pi * x * ncycle/n)
# synthetic time series ?
print "y: %s av %.3g variance %.3g" % (y, y.mean(), np.var(y))
print "Cost[j,k] = sum (y - av y)^2 --" # len * var y[j:k+1]
Cost = spread.spreads_allij(y)
print Cost # .round().astype(int)
jsplit = split4(Cost)
# split4: len [3 2 3 2] pos [ 0 3 5 8 10]
if plot:
import matplotlib.pyplot as pl
title += "\n lengths: %s" % np.diff(jsplit)
pl.title(title)
pl.plot(y)
for js, js1 in zip(jsplit[:-1], jsplit[1:]):
if js1 <= js: continue
yav = y[js:js1].mean() * np.ones(js1 - js + 1)
pl.plot(np.arange(js, js1 + 1), yav)
# pl.legend()
pl.show()
Dann wird der folgende Code macht Cost[]
für horizontale Linien lediglich, Neigung 0; wird es als Übung ausgegeben, es auf Linienabschnitte mit beliebiger Steigung in der Zeit O (n) zu verlängern.
""" spreads(all y[:j]) in time O(n)
define spread(y[]) = sum (y - average y)^2
e.g. spread of 24 hourly temperatures y[0:24] i.e. y[0] .. y[23]
around a horizontal line at the average temperature
(spread = 0 for constant temperature,
24 c^2 for constant + [c -c c -c ...],
24 * variance(y))
How fast can one compute all 24 spreads
1 hour (midnight to 1 am), 2 hours ... all 24 ?
A simpler problem: compute all 24 averages in time O(n):
N = np.arange(1, len(y)+1)
allav = np.cumsum(y)/N
= [ y0, (y0 + y1)/2, (y0 + y1 + y2)/3 ...]
An identity:
spread(y) = sum(y^2) - n * (av y)^2
Voila: the code below, all spreads() in time O(n).
Exercise: extend this to spreads around least-squares lines
fit to [ y0, [y0 y1], [y0 y1 y2] ... ], not just horizontal lines.
"""
from __future__ import division
import sys
import numpy as np
#...............................................................................
def spreads(y):
""" [ spread y[:1], spread y[:2] ... spread y ] in time O(n)
where spread(y[]) = sum (y - average y)^2
= n * variance(y)
"""
N = np.arange(1, len(y)+1)
return np.cumsum(y**2) - np.cumsum(y)**2/N
def spreads_allij(y):
""" -> A[i,j] = sum (y - av y)^2, spread of y around its average
for all y[i:j+1]
time, space O(n^2)
"""
y = np.asanyarray(y, dtype=float)
n = len(y)
A = np.zeros((n,n))
for i in range(n):
A[i,i:] = spreads(y[i:])
return A
Bisher haben wir eine n x n Kostenmatrix, Raum O (n^2). Um nach unten zu Raum O (nk) zu erhalten, genau hinsehen, um das Muster der Cost[i,j]
greift im dyn-prog-Code:
for nleft .. to n:
Cost_nleft = Cost[j,nleft ] -- time nleft or nleft^2
for k in 3 4 5 ...:
min [ C[k-1, j-1] + Cost_nleft[j] for j .. to nleft ]
Hier Cost_nleft
eine Zeile der vollständigen nxn-Kosten-Matrix ist, ~ n Segmente, generiert nach Bedarf. Dies kann in der Zeit O (n) für Liniensegmente erfolgen. Aber wenn "Fehler für ein Segment durch einen Satz von n Punkten dauert O (n) Zeit", scheint es, wir sind bis zu der Zeit O (n^3). Kommentare jemand?
'O (n^(2k))' oder 'O (n^2 * k)'? –
Muss die Funktion auch kontinuierlich sein? –
Es ist O (n^2 * k) und die Funktion muss nicht kontinuierlich sein. – user3019792