rtrwalker's profile - activity

The sage tutorial Loading and Attaching Sage files states that if I create a file such as walk.sage containing:

print "Hello world"
print 2^4

and then load or attach it:

sage: attach 'walk.sage'
Hello world
16

I should get a walk.py file containing python code in the same directory as walk.sage. I can find no such file (I've even searched my the rest of my computer). Does anyone know where the file is?

How can I get a list/tuple of the variables in a symbolic function with the same ordering as when the function was defined? e.g. for the function below I would want (z,t) not the alphabetically ordered (t, z) I get with .variables() of .arguments(). The ordering has to be stored/used somewhere in sage because I can differentiate with respect to z and get D0(z,t) as an answer where the '0' corresponds to 'z'.

sage: var('z,t')
(z, t)
sage: f = function('u',z,t)
sage: print (f)
u(z, t)
sage: f.variables()
(t, z)
sage: f.arguments()
(t, z)
sage: f.diff(z)
D[0](u)(z, t)
sage: f.diff(t)
D[1](u)(z, t)

How can I transform, step by step, the partial differential equation (pde) u.diff(z,z)k/gw == mvu.diff(t) into u.diff(Z,Z)=u.diff(Tv) where Z=z/H, Tv = cvt/H^2 and cv = k/(gwmv)? As you can see from my attempt below I can isolate cv; after that I have no idea how to do the change of variables. I can do this example by hand: let Z=z/H --> du/dz = du/dZ * dZ/dz --> du/dz = du/dZ *(1/H) ...etc. But Other examples may not be so easy (and I don't want to do it by hand!).

var('u,z,t,Z,H,Tv,k,gw,mv,cv')
u = function('u',z,t)
eq1 = u.diff(z,z)*k/gw == mv * u.diff(t)
print(eq1)
mycoeff=eq1.left_hand_side().coeff(u.diff(z,z),1)
print(mycoeff)
eq2 = eq1.divide_both_sides(mycoeff)
print(eq2)
mycoeff = eq2.right_hand_side().coeff(u.diff(t),1)
eq3 = cv == 1/mycoeff
print(eq3) 
eq4 = solve(eq3,eq3.right_hand_side().default_variable(), solution_dict=true)
print(eq4)
eq5 = eq2.subs(eq4[0])
print(eq5)

which gives me the following output:

k*D[0, 0](u)(z, t)/gw == mv*D[1](u)(z, t)
k/gw
D[0, 0](u)(z, t) == gw*mv*D[1](u)(z, t)/k
cv == k/(gw*mv)
[{gw: k/(cv*mv)}]
D[0, 0](u)(z, t) == D[1](u)(z, t)/cv

To get your 'answer code to work I need: 'from sympy.utilities.codegen import codegen' and 'from sympy import Eq'

Thanks @DSM. I've tried to generalize your solution with the following:

def get3(fn, constraint_nums, x0, iobj=0):
    """
    Constrained mimimization of a function that is output, along with constraint equations, from function fn.

    Arguments:
        fn - function that returns a tuple containing objective function and constraint equations

        constraint_nums - tuple or list containing index positions of constraint equations in fn output

        x0 - initial guess

        iobj - index position of object function in fn output.  default = 0
    """  

    a = minimize_constrained(lambda x,j=iobj: fn(x)[j], 
                             [lambda x,j=i: fn(x)[j] for i in constraint_nums],
                             x0)
    return a

running get3(func_wrap,range(1,5),[2,3]) results in the correct (45.0, 6.25). I had some trouble with the list comprehension of the lambda functions for the constraints. If I used [lambda x: fn(x)[i] for i in constraint_nums] then the [i] wasn't hardcoded in and each of my lambda equations had a parameter i in them; the list comphehension looped through the i's meaning all the lambda functions were the same with i equal to the final i such that only my last constraint equation was used. [lambda x,j=i: fn(x)[j] for i in constraint_nums] seems to work.

Two further questions: 1.Is there a way to determine the number of values returned by a function without actually evaluating the function? If there was I could further generalise get3 by giving constraint_nums a default value. 2. What actually gets cached when I use FunctionCached? Is everything in the function stored or just the output values? When should I be using fn.clear_cache()?

On a different note, within get2 I did try to use a list construction to specify my constraints: [lambda x: fn(x)[i] for i in range(1,5)]. This gave the answer (739.437607818, -608.579643471) for the non-cached version and still gave errors for the cached version; not sure how/why that approach stuffed up.

@DSM I can successfully run your example but if I adapt your solution to the example in the minimize_constrained documentation then I get an error for the cached version. For the non-cached version I get the correct result (45.0, 6.25).

def func(p):
    f = -p[0]-p[1]+50
    c_1 = p[0]-45
    c_2 = p[1]-5
    c_3 = -50*p[0]-24*p[1]+2400
    c_4 = -30*p[0]-33*p[1]+2100
    return f, c_1, c_2, c_3, c_4

@cached_function    
def func_cached(p):
    f = -p[0]-p[1]+50
    c_1 = p[0]-45
    c_2 =p[1]-5
    c_3 =-50*p[0]-24*p[1]+2400
    c_4 =-30*p[0]-33*p[1]+2100
    return f, c_1, c_2, c_3, c_4

def get2(fn):
    a = minimize_constrained(lambda x: fn(x)[0], 
                             [lambda x: fn(x)[1],lambda x: fn(x)[2],lambda x: fn(x)[3],lambda x: fn(x)[4]],
                             [2,3])
    return a

If I run get2(func) the result is (45.0, 6.25). If I run get2(func_cached) I get errors. I did try rewritting func using x, y arguments (and fn(*x) in get2) but I got the same results.

I'm trying to think of a way to do constrained optimization where both the objective function that I want to minimise and my constraint function are calculated within the same overarching function. Lets say I have a function like the following:

def complicated_function(x,y,z):
    #...some lengthy and complicated actions
    return f, g

Is there a way to minimize f by changing x,y,z subject to g>=0?

Looking at the documentation for sage.numerical.optimize.minimize_constrained it looks as if I have to define all my constraint equations individually with the varibales as a tuple. I could maybe wrap my complicated_function up like this:

def funcObjective(x):
    f, g = complicated_function(x[0],x[1],x[2])
    return f
def funcConstraint(x):
    f, g = complicated_function(x[0],x[1],x[2])
    return g

But that would mean complicated_function would be called multiple times for each optimization iteration which seems highly inefficient. The problem would only get worse if there were more contraints. Any ideas how to define the constraint without reevaluating complicated_function?

On my internet wanderings I found the recently (June 6, 2011) released pyOpt 1.0 (journal article) which at first glance looks well suited to the problem. I see OpenOpt is an experimental package for sage. I'm not sure if openOpt is suitable; the pyOpt documentation is, at first glance, clearer. Any chance pyOpt could be made an optional package with sage, it's published under the GNU Lesser General Public License?

In addition to my first comment to @kcrisman for the case where g = 1/4 if I run solve() on the incomplete answer that solve() originally gave me then I get the correct answer: solve(1/2*sqrt(4*x + 1) == h,x) gives: x == h^2 - 1/4 Why doesn't solve() give me that answer in the first place?

Also, in my original question if I replace sol2 = solve(eqs[1].subs(sol1[0]),x_0) with sol2 = solve(eqs[1].subs(sol1[0]).expand(),x_0), i.e. expand before solving, I get x0 = some function of x0 which is wrong.

@kcrisman here's a smaller example that starts misbehaving for reals: reset() forget() var('h,x,g') assume (h>0) print(solve(sqrt(x+g)==h,x)[0]) print(solve(sqrt(x+1)==h,x)[0]) print(solve(sqrt(x+2.0)==h,x)[0]) print(solve(sqrt(x+1/4)==h,x)[0]) print(solve(sqrt(x+2.5)==h,x)[0]) myrr.<rr> = PolynomialRing(RR) print(solve(sqrt(x+rr)==h,x)[0]) which gives: x == h^2 - g x == h^2 - 1 x == h^2 - 2 1/2*sqrt(4*x + 1) == h 1/2*sqrt(2*x + 5)*sqrt(2) == h and a traceback error on the last example with rr (not enough space in comment to show). Other polynomial rings (CC,ZZ,QQ) gave the correct answer.