In [ ]:

    

import numpy as np


    

In [ ]:

    

from scipy.special import comb as nchoosek # n Choose k (n ueber k)

# As seen in A.4 jasniwiecz
# Intervalltransformation (Gleichung (4.36))
# p => [p_l, p_u] (p => [p_min, 1]) wird zu p^ => [-1,1] return neue Matrizen(!) mit diesen p^
#
# Q_in ist Liste aus Faktormatrizen für Polynom in p
# Q_out ist Liste aus Faktormatrizen für Polynom in p^
# a = (p_u - p_l)
# b = (p_u + p_l)
def transQ(Q_in, pl=0.01, pu=1):
    a = pu - pl
    b = pu + pl
    m = len(Q_in) # Anzahl der Koeffizienten in Q_in
    n = Q_in[0].shape

    #print Q_in
    Q_out = [np.zeros(Q_in[0].shape)] * m
    #print Q_out
    for j in range(0, m):
        #print j
        #Q_out[j] = Q_in[j]
        for i in range(j, m):
            #print Q_out[j]
            Q_out[j] = Q_out[j] + 2**(-i) * b**(i-j) * nchoosek(i, j) * Q_in[i]
        Q_out[j] = a**j*Q_out[j]
    return Q_out


    

In [ ]:

    

%load_ext oct2py.ipython


    

In [ ]:

    

p_min = np.random.rand(1)
n = int(100*np.random.rand(1))
print "n:", n
print "p_min:", p_min
R = [np.random.rand(n,n),
     np.random.rand(n,n),
     np.random.rand(n,n)]

%octave addpath('~/jupyter_nb/optimization')
%octave_push p_min
%octave_push R
%octave a = (1-p_min)/2;
%octave b = (1+p_min)/2;
%octave X = transQ(R, a, b);
%octave_pull X

X2 = transQ(R, pl=p_min, pu=1)

for i in range(3):
    print "Test:", i
    np.testing.assert_array_equal(X[0,i], X2[i])


    

In [ ]:

    

import math

def calcQSum(Q_in):
    from numpy import concatenate as con
    m = len(Q_in)
    n = Q_in[0].shape
    if m is 1:
        #([[2*Q_in[0],   np.zeros(n)], 
        #  [np.zeros(n), np.zeros(n)]])
        Q_sum = con((con((2*Q_in[0], np.zeros(n)), axis=1),
                     con((np.zeros(n), np.zeros(n)), axis=1)), axis=0)
    elif m is 2:
        # [[2*Q_in[0], Q_in[1]],
        #  [Q_in[1],   np.zeros(n)]]
        Q_sum = con((con((2*Q_in[0], Q_in[1]), axis=1),
                     con((Q_in[1], np.zeros(n)), axis=1)), axis=0)

    else:
        # [[2*Q_in[0], Q_in[1]],
        #  [Q_in[1], 2*Q_in[2]]]
        # m here is m+1 from Diss (starts at 0 there), so even and odd is interchanged
        Q_sum = con((con((2*Q_in[0], Q_in[1]), axis=1),
                     con((Q_in[1], 2*Q_in[2]), axis=1)), axis=0)
        
    for i1 in range(4, m+1, 2):
        if i1 is not m:
            #Q_sum = [[Q_sum], [[np.zeros((i1/2-1)*n, n)], [Q_in[i1]]], 
            #         [np.zeros(n, (i1/2-1)), Q_in[i1], 2*Q_in[i1+1]]]
            Q_sum = con((con((Q_sum, con((np.zeros((int(math.ceil(i1/2-1)*n[0]), n[1])), Q_in[i1-1]), axis=0)), axis=1),
                         con((np.zeros((n[0], int(math.ceil(i1/2-1)*n[1]))), Q_in[i1-1], 2*Q_in[i1]), axis=1)), axis=0)
        else:
            #Q_sum = [[Q_sum], [[np.zeros((i1/2-1)*n, n)], [Q_in[i1]]], 
            #         [np.zeros(n, (i1/2-1)), Q_in[i1], np.zeros(n)]]
            Q_sum = con((con((Q_sum, con((np.zeros((int(math.ceil(i1/2-1)*n[0]), n[1])), Q_in[i1-1]), axis=0)), axis=1),
                         con((np.zeros((n[0], int(math.ceil(i1/2-1)*n[1]))), Q_in[i1-1], np.zeros(n)), axis=1)), axis=0)
    Q_sum = 0.5*Q_sum
    return Q_sum


    

In [ ]:

    

n = int(100*np.random.rand(1))
m = int(10*np.random.rand(1))+1

#n=2
#m=2

print "n:", n
print "m:", m
R = [np.random.rand(n,n) for i in range(m)]

%octave addpath('~/jupyter_nb/optimization')
%octave_push R
%octave R;
%octave X = calcQsum(R);
%octave_pull X

X2 = calcQSum(R)

for i in range(m):
    print "Test:", i
    #print "X:", X
    #print "X2:", X2
    np.testing.assert_array_equal(X, X2)


    

In [ ]:

    




    

In [ ]:

    

import boris_tools as tools
n = int(10*np.random.rand(1))+1
print n
A = np.array([[0., 1., 0.],
              [1., 0., 1.],
              [0., 0., -0.005]])

A = np.random.rand(n,n)
#a = -A[-1][:].T #!!!!

b = np.array([[0], [1], [1.]])
b = np.random.rand(n,1)

d = 0
c = np.array([[1], [0], [0]])
c = np.random.rand(1,n)

roots_p1 = [-1, -1+1j, -1-1j]
roots_p1 = -10*np.random.rand(n)

roots_pmin = [-3, -3+2j, -3-2j]
roots_pmin = -20*np.random.rand(n)

p_min = 0.1

%octave addpath('~/jupyter_nb/optimization')
%octave_push A
%octave_push b
%octave_push roots_p1
%octave_push roots_pmin
%octave_push p_min

%octave [k0, k1] = k_explizit_Ab3(roots_p1, roots_pmin, p_min, A, b);
%octave_pull k0
%octave_pull k1

(k0, k1)
#print tools.k_explizit_Ab(roots_p1, roots_pmin, p_min, A, b)
k0x, k1x = tools.k_explizit_Ab2(roots_p1, roots_pmin, p_min, A, b)

np.testing.assert_allclose(k0, k0x)
np.testing.assert_allclose(k1, k1x)


    

In [ ]: