Experiments with S(x) = sin(x) + sin(2x)/2 + sin(4x)/4 + ...

sumsin_mpmath.py

"""Experiments with the 'sumsin' function: maximization and plotting

   S(x) = sin(x) + sin(2x)/2 + sin(4x)/4 + ...

it appears that S(x) reaches maximum at
   xmax = 0.8905302010175791857059461558*7027*
whereas
    36pi/127=0.8905302010175791857059461558*9025*

Plots of the function:
* One period of S(x): http://imgur.com/yDcTy9E
* Behavior of S(x) near 36pi/127, coordinates relative to (36pi/127,S(36pi/127)): http://imgur.com/Btkd1kH
* Same plot, zoomed to show that left peak is higher: http://imgur.com/hbpjufM

Language: Python 3
Requires: mpmath
"""
from mpmath import mp


def sumsin(x):
    """Calculate S(x) = sin(x) + sin(2x)/2 + ...."""
    pi2 = mp.pi*2
    if x < 0 or x > pi2:
        raise ValueError("x must be in 0...2pi range, for convenience")
    
    s = mp.mpf(0.0)
    k = mp.mpf(1.0)

    while k > mp.eps:
        s += mp.sin(x) * k
        x *= 2
        if x > pi2:
            x = x-pi2
        k /= 2
    return s

def _intervals_above_treshold( xs, ys, treshold ):
    """Returns list of (x0,x1) pairs, such that corresponding y is above treshold"""
    a = None
    xprev = None
    for xi, yi in zip(xs, ys):
        if a is None: #not yet entered range
            if yi >= treshold:
                #start range
                a = xprev if xprev is not None else xi
        else:
            if yi <= treshold:
                #end range
                yield (a, xi)
                a = None
        xprev = xi
    #final: end range if it has started
    if a is not None:
        yield (a, xi)
        
def search_intervals( func, a, b, eps, subdivisions=50, percent = 0.5, verbose=True, parallel=True ):
    """Extremely redundant maximizer.
    Search for maximum of func(x) when x in [a,b]
    With X tolerance eps.

    On each step, determine set of intervals, where maximum is searched,
    subdivide them, then take only sub-intervals above treshold.

    Can use parallel processing.
    """
    
    a,b,eps, percent = map(mp.mpf, (a,b,eps,percent))
    if subdivisions < 2: raise ValueError("Subdivisions must be > 2")
    if percent <= 0 or percent >=1: raise ValueError("percent must be in (0,1)")

    print_ = print if verbose else lambda *x:None
    
    print_ ("Start search fom {a} to {b}, subdivisions: {subdivisions}".format(**locals()))
    if parallel:
        from multiprocessing import Pool
        pool = Pool()
        map_ = pool.map
        print_ ("Using  multiprocessing")
    else:
        def map_(f,arg): return list(map(f,arg))

    #list of intervals, where maximum is searched for.
    intervals = [ (a,b) ]
    step = 0
    while True:
        print_ ("Step: {step}, intervals: {nintervals}".format(nintervals=len(intervals), **locals()))
        span = intervals[-1][1] - intervals[0][0]
        if span <= eps:
            break
        print_ ("    X span:", float(span))
        
        xvalues = sum( (mp.linspace(ai,bi,subdivisions) for (ai, bi) in intervals),
                       [] )
        
        yvalues = map_(func, xvalues)
        
        ymin = min(yvalues)
        ymax = max(yvalues)
        print_("    Y span: {ymin} ... {ymax}".format(**locals()))
        ytreshold = ymax - (ymax-ymin)*percent
        print_ ("    Y treshold:", ytreshold)
        intervals = list(_intervals_above_treshold( xvalues, yvalues, ytreshold ))
            
        step += 1
    
    return intervals[-1][1], intervals[0][0]

def plot_near_extremum():
    """Make the plot of S(x), deeply zoomed around x=36pi/127"""
    from matplotlib import pyplot as pp
    N = 1500                  #number of points
    xcenter = mp.pi*36/127    #central point
    dx = mp.mpf('-2e-29')*1.5 #horizontal span to plot
    
    
    xs = mp.linspace(xcenter - dx, xcenter + dx, N)    
    ys = [sumsin(xi) for xi in xs]
    
    ymin = min(ys)
    ymax = max(ys)
    yspan = ymax - ymin
    
    dxs_float = list(map(float, mp.linspace(- dx, dx, N)))
    dys_float = [float(yi-ymax) for yi in ys]

    pp.plot( dxs_float, dys_float )
    pp.grid()
    pp.show()

def plot_overview():
    """Plot whole period of S(x)"""
    from matplotlib import pyplot as pp
    N = 1000
    
    xs = [float(xi) for xi in mp.linspace(0, 2*mp.pi, N)]
    ys = [float(sumsin(xi)) for xi in xs]
    
    pp.plot( xs, ys )
    pp.grid()
    pp.show()

def search_maximum():
    """Search for maximum of S(x)"""
    pi36_127 = mp.pi*36/127
    
    print ("S(36pi/127={s0}".format(s0 = sumsin(pi36_127)))

    #set disired tolerance 2 digits less than 
    xtol = mp.eps*100
    
    print ("Searching argmax(S(x)) with tolerance {xtol}".format(xtol=float(xtol)))
    a, b = search_intervals( sumsin, 0, mp.pi*2, xtol,
                             subdivisions=50,  #on each step, divide each interval in this number of points.
                                               #more - slower & reliable, less - faster & error-prone
                             percent = 0.2,    #on each step, percent of Y span taken to the next step
                             verbose = True,
                             parallel = True)
    c=(a+b)/2
    print ("Finished search, xmax={c}".format(c=c))
    print ("X tolerance = {xtol}".format(xtol=float(b-a)))
    print ("xmax - 36pi/237 = {diff}".format(diff=float(c-pi36_127)))
    print ("S(xmax) - S(36pi/127) = {df}".format( df = float(sumsin(c) - sumsin(pi36_127))))

if __name__ == "__main__":
    #more then 31 to get interesting result
    mp.dps=40
    print ("Set precision to {d} digits".format(d=mp.dps))
    
    #plot_all()
    #plot_near_extremum()
    search_maximum()