4448-3. Optimization.sagews

I only have a few minutes to spend on this, but here is "Sage-a-love-story" based on http://www.stiglerdiet.com/2012/02/10/mathematica-a-love-story/

 

1. Powerful Symbolic Computations

f(x) = 1/((x+2)^2*(x^2+1))
g = f.partial_fraction().expand()
g

\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ -\frac{4 \, x}{25 \, {\left(x^{2} + 1\right)}} + \frac{4}{25 \, {\left(x + 2\right)}} + \frac{1}{5 \, {\left(x + 2\right)}^{2}} + \frac{3}{25 \, {\left(x^{2} + 1\right)}}

var('p')
solve(20 * (3*(1-p)^4*p^2 - 3*(1-p)^2*p^4) == 0, p)

\newcommand{\Bold}[1]{\mathbf{#1}}\left[p = \left(\frac{1}{2}\right), p = 0, p = 1\right]

var('q,r,theta,a,R,phi')

assume(R-a>0)
assume(R>0)
assume(a>0)

integrate(integrate(integrate(
    q^2/r^4*r^2*sin(theta), (r, R, a)), (theta, 0, pi)), (phi, 0, 2*pi))

\newcommand{\Bold}[1]{\mathbf{#1}}4 \, {\left(\frac{1}{R} - \frac{1}{a}\right)} \pi q^{2}

2. Functional Programming

Sage (via Python) supports procedural and functional programming.  And like R... Sage includes R!

x = 0
while x<=1:
    print x
    x += 1

0 1

reset('r')

%r
v <- c(1,2,7)
sd(v)

[1] 3.21455

p = [2, 4, 6, 8]
d = [4, 14, 10, 16]

def accumulate(v):
    if not len(v): return []
    w = [v[0]]
    for a in v[1:]:
        w.append(w[-1]+a)
    return w

x = [1,4,3,2]    
accumulate([p[i-1] for i in x])

\newcommand{\Bold}[1]{\mathbf{#1}}\left[2, 10, 16, 20\right]

Or just use numpy:

import numpy as np
v = np.array([p[i-1] for i in x]).cumsum(); v

\newcommand{\Bold}[1]{\mathbf{#1}}\verb|[|\phantom{x}\verb|2|\phantom{x}\verb|10|\phantom{x}\verb|16|\phantom{x}\verb|20]|

It seems like the original article is wrong... ?

import numpy as np
w = v - np.array([d[i-1] for i in x]); w

\newcommand{\Bold}[1]{\mathbf{#1}}\verb|[-2|\phantom{x}\verb|-6|\phantom{xx}\verb|6|\phantom{xx}\verb|6]|

[max(int(a), 0) for a in w]

\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, 0, 6, 6\right]

sum(_)

\newcommand{\Bold}[1]{\mathbf{#1}}12

3. Optimization

cvxopt can do that minimize problem fine, but it is currently very difficult to use from Sage, to put it mildly.   The docs to illustrate how to do it are here: 

http://abel.ee.ucla.edu/cvxopt/userguide/coneprog.html#linear-programming

model = MixedIntegerLinearProgram(maximization=False) 
# don't miss the [0] ...
c = model.new_variable()[0]
o = model.new_variable()[0]
model.set_objective(.5 * c + .5 * o)
model.add_constraint(172 * c + 1742 * o, min= 2000)
model.add_constraint(.019 * c + .015 * o, min = .06)
model.add_constraint(c, min=0)
model.add_constraint(o, min=0)
model.show()

Minimization: 0.5 x_0 +0.5 x_1 Constraints: 2000.0 <= 172.0 x_0 +1742.0 x_1 0.06 <= 0.019 x_0 +0.015 x_1 0.0 <= x_0 0.0 <= x_1 Variables: x_0 is a continuous variable (min=0.0, max=+oo) x_1 is a continuous variable (min=0.0, max=+oo)

print 'Objective Value:', model.solve()

Objective Value: 1.6744216528

print 'c = %s' % model.get_values(c)
print 'o = %s' % model.get_values(o)

c = 2.44183760404 o = 0.907005701553

pcost dcost gap pres dres k/t 0: -5.9153e-01 -6.7231e+00 9e+00 1e+00 2e+03 1e+00 1: 7.9188e-02 -5.9908e-01 1e+00 1e-01 2e+02 1e-02 2: 5.3296e-04 -7.6202e-03 1e-02 1e-03 2e+00 8e-04 3: 3.6898e-04 -6.5570e-05 5e-04 6e-05 1e-01 1e-05 4: 1.3107e-06 -3.0559e-05 4e-05 5e-06 1e-02 6e-06 5: 1.3482e-08 -3.0563e-07 4e-07 5e-08 1e-04 6e-08 6: 1.3482e-10 -3.0563e-09 4e-09 5e-10 1e-06 6e-10 7: 1.3482e-12 -3.0563e-11 4e-11 5e-12 1e-08 6e-12 Optimal solution found. [ 4.70e-13] [ 2.23e-12]

Memoization is dead easy in Sage/Python (and you can even make it cache to a file for longterm use if you want with disk_cached_function...)

@cached_function
def f(n):
    sleep(1)
    return n*n

time f(2)

\newcommand{\Bold}[1]{\mathbf{#1}}4

Time: CPU 0.00 s, Wall: 1.00 s

time f(2)

\newcommand{\Bold}[1]{\mathbf{#1}}4

Time: CPU 0.00 s, Wall: 0.00 s

4. Graphics

Graphics ready to include in Latex, etc.

var('x')
plot(sin(x), (x,0,10)).save('plot.pdf')

5. Documentation

10,000+ pages of Sage documentation is available online. 

6. Naming Conventions

In Sage, names of functions are usually complete English words spelled out, with underscores.

7. Interactivity

@interact
def f(theta=(pi/100, pi)):
    x = var('x')
    show(plot(sin(1/x), (x, -theta, theta)), figsize=4)

var('x')
show(animate([plot(sin(1/x), (x, -theta, theta)) for theta in [pi/100, pi/10, .., pi]], figsize=4))