Complex vectors and shape for arrays

2020-05-03 10:14发布

I am slightly new to python and I am trying to convert some code.This is an approximation method. Which isn't important. In my oddev function I get returned

        c2[1:modes+1] = v* 1j
ValueError: could not broadcast input array from shape (25) into shape (25,1) 

When I do this Matlab I believe it automatically casts it, and will store the complex array. The function is a getting the coefficient from a partial sine transform to do this. At first I tried storing the random matrix which just an array using np.matlib method and this had the same shape but I believe I will lose the real values of the filter when I cast it. How do I store this?

import math
import numpy as np
def quickcontmin(datain):

n = np.shape(datain)[0]
m = math.floor(n / 2)
modes = math.floor(m / 2)
addl = 20
nn = 20 * n
chi = 10 ** -13

def evenhp(xv):
    "Even high pass"
    n1 = np.shape(xv)[0]
    vx = np.array(xv[:-1])
    vx = vx[::-1]
    c1 = np.append(xv,vx)
    c1 = np.fft.fft(c1)       
    c1[0:modes-1] = 0.0
    c1[-1 - modes + 2:-1] = 0.0
    evenl = np.real(np.fft.ifft(c1))
    even = evenl[0:n1-1]
    return even

def evenhpt(xv):
    " Transpose of EvenHP"
    n1 = np.shape(xv)[0]
    xy = np.zeros((n1- 2, 1))
    c1 = np.append(xv,xy)
    c1 = np.fft.fft(c1)
    c1[0:modes-1] = 0.0
    c1[-1 - modes + 1:-1] = 0.0
    evenl = np.real(np.fft.ifft(c1))
    even = evenl[0:n1-1]
    even[1:-2] = even[1:-2] + evenl[-1:-1:n1+1]
    return even``

def evenlp(xv):
    " Low pass cosine filter"
    n1 = np.shape(xv)[0]
    vx = np.array(xv[:-1])
    vx = vx[::-1]
    c1 = np.append(xv,vx)
    c1 = np.fft.fft(c1)
    c1[modes + 1:-1 - modes + 1] = 0.0
    evenl = np.real(np.fft.ifft(c1))
    even = evenl[0:n1-1]
    return even

def oddev(xv):
    "Evaluate the sine modes on the grid"
    c2 = np.zeros((2 *n - 2, 1))*1j
    v = np.array(xv[:])
    v1 = v[:-1]
    v1 = v[::-1]
    c2[1:modes+1] = v* 1j
    c2[-1 - modes + 1:-1] = -v1* 1j
    evall = np.fft.ifft(c2) * math.sqrt(2 * n - 2)
    eva = evall[0:n-1]
    return eva

def oddevt(xv):
    " Transpose the sine modes on the function OddEv"
    c1 = np.array(xv[1:-2])
    c1 = np.insert(c1,0.0,0)
    c1 = np.append(c1,0.0)
    c1 = np.append(c1,xv[-2:-1:2])
    c1a = np.divide(np.fft.fft(c1),math.sqrt(2 * n - 2))
    fcoef = np.imag(c1a[1:modes])
    return fcoef

def eextnd(xv):
    "Obtain cosine coefficients and evalue on the refined grid"  
    vx = np.array(xv[:-1])
    vx = vx[::-1]
    c1  = np.append(xv,vx)
    c1 = np.fft.fft(c1)
    cL = np.zeros((2*nn-2,1))
    cL[0:modes-1] = c1[0:modes-1]
    cL[-1 - modes + 1:-1] = c1[-1 - modes + 1:-1]
    evenexL = np.multiply(np.fft.ifft(cL) , (nn - 1) / (n - 1))
    evenex = evenexL[0:nn-1]
    return evenex

def oextnd(xv):
    "Evaluate sine coefficients on the refined grid"
    c2 = np.zeros((2 * nn - 2, 1))
    c2[0] = 0.0
    c2[1:modes + 1] = np.multiply(xv[0:-1],1j)
    c2[-1 - modes + 1:-1] = np.multiply(-xv[-1:-1:1],1j)
    evall = np.real(np.multiply(np.fft.ifft(c2), math.sqrt(2 * n - 2) * (2 *nn - 2) / (2 * n - 2)))
    oox = evall[0:nn-1]
    return oox

dc = evenlp(datain)
#L in paper, number of vectors used to sample the columnspace
lll = round(4 * math.log(m )/ math.log(2)) + addl
lll = int(lll)
#The following should be straightforward from the psuedo-code
w=2 * np.random.rand(modes , lll) - 1 
p=np.matlib.zeros(shape=(n,lll))
for j in range(lll):
    p[:,j] = evenhp(oddev(w[:,j]))
q,r = np.linalg.qr(p , mode='reduced')
z = np.zeros(shape=(modes,lll))
for j in range(lll):
    z[:,j]= oddevt(evenhpt(q[:,j]))
un,s,v = np.linalg.svd(z,full_matrices='False')
ds=np.diag(s)
aa=np.extract(np.diag(s)>(chi))
aa[-1] = aa
aa = int(aa)
s = 0 * s
for j in range(aa):
    s[j,j] = 1.0 / ds(j)
#find the sine coefficents
b=un*s* v.T* q.T* evenhp(datain)
#Constructing the continuation
exs=oddev(b)
pexs = evenlp(exs)

dataCont=exs-pexs+dc
dataCont[n+1:2*n-2]=-exs[-2:-1:1]-pexs[-2:-1:1]+dc[-2:-1:1]
#Evaluate the continuation on the refined grid
dataRefined=eextnd(dc-exs)+oextnd(b)

return dataRefined, dataCont   



n1 = 100
t = np.linspace(0,2*math.pi,n1)
y = np.sin(t)
data = quickcontmin(y)    
dc1 = data[1]
dc1 = dc1[0:n1-1]`

1条回答
【Aperson】
2楼-- · 2020-05-03 11:14

Replacing c2[1:modes+1] = v* 1j by c2[1:modes+1, 0] = v* 1j should fix that specific error. More consistent would be to replace:

v = np.array(xv[:])
v1 = v[:-1]
v1 = v[::-1]

by

v = xv
v1 = v[:-1]

v is already a column vector so you don't need to transform it into a 1d vector when you later need a column vector.

查看更多
登录 后发表回答