'IS there something im doing wrong with my python code?
I am trying to convert a matlab script into python. I dont seem to see what im doing wrong but i run into some errors.
import numpy as np
nargin = 10.
if nargin == 0:
rseed = 56
dx = 25
dz = dx #dx=dz
ax = 750 #initial 750....a = 1;
az = 200 #initial 200
nu = 0.6 #initial 0.6
ix = 6000 #6000
iz = 6000 #6000
pop = 2
med = 3
vel = np.array([6000, 6600, 7000, 2000, 5000]) #6000 6600 7000 2000 5000
frac = np.array([0.5, 0.5, 0.5, 0.3, 0.1]) #.5 .5 .5 .3 .1
if nargin == 10:
frac = [1,2,3] #the same as [.16 .33 .5]
vel = [2000,3000,4000]
if pop == 1:
D = 3 - nu #POP=GAUSSIAN
elif pop == 2:
D = 3 - nu/2 #POP=PDF
#ceil runs off to the next integer.
nx = np.ceil(ix/dx)
nz = np.ceil(iz/dz)
if (nx/2)!= round(nx/2):
ddx=1
nx=nx+ddx
else:
ddx=0
if (nz/2) != round(nz/2):
ddz = 1
nz = nz+ddz
else:
ddz=0
x=np.arange(1,nx+1,1)*dx
z=np.arange(1,nz+1,1)*dz
dkx = 1/(nx*dx)
dkz = 1/(nz*dz)
a1= -nx/2
a2= (nx/2)-1
#aa= np.arange(0, nx/2, -nx/2)
[kx, kz] = np.meshgrid(2*np.pi*dkx*np.arange(1,nx/2,-nx/2), 2*np.pi*dkz*np.arange(1,nz/2,-nz/2))
k = np.sqrt(kx**2 *ax**2 + kz**2 * az**2)
if med==1:
# Gaussian Chh
#SQRT of the FOURIER TRANSROMF OF Gaussian CORR fkt.
#disp([' Calculating filter : Chh=Gaussian'])
expcorr=((ax*az)/2) * np.exp(-(k**2/4)) # (Exact F(C_hh(x))
expcorr=expcorr/np.max(np.max(expcorr)) # normalizing sqr(F(C_hh))
expcorr=np.sqrt(expcorr)
if med==2:
#Exponential Chh
#SQRT of the FOURIER TRANSROMF OF exp CORR fkt.
#disp([' Calculating filter : Chh=exponential'])
#expcorr=(ax^2+az^2)./((1+k.^2).^(1.5)); % (Exact F(C_hh(x))
expcorr=1/((1+k**2)**(1.5)) # (Exact F(C_hh(x))
expcorr=expcorr/np.max(np.max(expcorr)) # normalizing sqr(F(C_hh))
expcorr=np.sqrt(expcorr)
if med==3:
# von Karman
# SQRT of the FOURIER TRANSROMF OF vonk CORR fkt.
# disp([' Calculating filter : Chh=vonKarman'])
#disp([' Fractal Dimension : ',num2str(D)])
#expcorr=((ax*az)/2)./((1+k.^2).^(nu+1)); # (Exact F(C_hh(x))
expcorr=1/((1+k**2)**(nu+1)) # (Exact F(C_hh(x))
expcorr=expcorr/np.max(np.max(expcorr)) # normalizing sqr(F(C_hh))
expcorr=np.sqrt(expcorr)
Below is the matlab code:
function [rdata,x,z,data,expcorr]=vonk2d(rseed,dx,dz,ax,az,ix,iz,pop,med,nu,vel,frac)
if nargin==0,
help vonk2d
rseed=56;
dx=25;
dz=dx;%dx=dz
ax=750;%initial 750....a = 1;
az=200; %initial 200
nu=0.6; %initial 0.6
ix=6000; %6000
iz=6000; %6000
pop=2;
med=3;
vel=[6000 6600 7000 2000 5000]; %6000 6600 7000 2000 5000
frac = [0.5 0.5 0.5 0.3 0.1]; %.5 .5 .5 .3 .1
end
if nargin==10,
frac=[1,2,3]; % the same as [.16 .33 .5]
vel=[2000,3000,4000];
end
if pop==1, % POP=GAUSSIAN
D=3-nu;
elseif pop==2 % POP=PDF
D=3-nu/2;
end
%ceil runs off to the next integer.
nx=ceil(ix/dx);
nz=ceil(iz/dz);
%disp([' VONK2D : 2D RANDOM MEDIA GENERATOR'])
%disp([' Using (nx,nz)=(',num2str(nx),',',num2str(nz),')'])
if (nx/2)~=round(nx/2), ddx=1; nx=nx+ddx; else ddx=0; end
if (nz/2)~=round(nz/2), ddz=1; nz=nz+ddz; else ddz=0; end
% SPACEDOMAIN GRID
%[x,z]=meshgrid(dx*[-nx/2:1:nx/2-1],dz*[-nz/2:1:nz/2-1]);
x=[1:1:nx]*dx;
z=[1:1:nz]*dz;
% WAVENUMBER DOMAIN GRID
dkx=1/(nx*dx);
dkz=1/(nz*dz);
[kx,kz]=meshgrid(2*pi*dkx*[-nx/2:1:nx/2-1],2*pi*dkz*[-nz/2:1:nz/2-1]);
k=sqrt(kx.^2.*ax.^2+kz.^2.*az.^2);
if med==1,
%
% Gaussian Chh
%
% SQRT of the FOURIER TRANSROMF OF Gaussian CORR fkt.
%disp([' Calculating filter : Chh=Gaussian'])
expcorr=((ax*az)/2).*exp(-(k.^2./4)); % (Exact F(C_hh(x))
expcorr=expcorr./max(max(expcorr)); % normalizing sqr(F(C_hh))
expcorr=sqrt(expcorr); %
end
if med==2,
%
% Exponential Chh
%
% SQRT of the FOURIER TRANSROMF OF exp CORR fkt.
%disp([' Calculating filter : Chh=exponential'])
%expcorr=(ax^2+az^2)./((1+k.^2).^(1.5)); % (Exact F(C_hh(x))
expcorr=1./((1+k.^2).^(1.5)); % (Exact F(C_hh(x))
expcorr=expcorr./max(max(expcorr)); % normalizing sqr(F(C_hh))
expcorr=sqrt(expcorr); %
end
if med==3,
%
% von Karman
%
% SQRT of the FOURIER TRANSROMF OF vonk CORR fkt.
%disp([' Calculating filter : Chh=vonKarman'])
%disp([' Fractal Dimension : ',num2str(D)])
%expcorr=((ax*az)/2)./((1+k.^2).^(nu+1)); % (Exact F(C_hh(x))
expcorr=1./((1+k.^2).^(nu+1)); % (Exact F(C_hh(x))
expcorr=expcorr./max(max(expcorr)); % normalizing sqr(F(C_hh))
expcorr=sqrt(expcorr); %
end
% DATA
%a = 500; b = 5000;
rand('seed',rseed);
data=(rand(nz,nx));%-0.5); %random distribution are between 0 & 1.
%data = (b-a).*data+a; %to spread values between 3000 & 2000.
% GOING TO FOURIER DOMAIN
%disp([' Going to Fourier Domain(fft)'])
fdata=fftshift(fft2(data));
% MULTIPLYING fdata by sqrt(C_hh)
%disp([' Applying filter'])
newfdata=fdata.*expcorr;
%FROM FOURIER TO SPACE DOMAIN
%disp([' Going to spacedomain (ifft)'])
randdata=real(ifft2(fftshift(newfdata)));
rdata=randdata;
%return
if pop==1,
% scaling filed according to vel
rdata=randdata.*2 *vel(1);
%data=data.*2 *vel(1);
end
if pop==2
%disp([' Using pdf population'])
frac=cumsum(frac./sum(frac)); % normalize and cumsum
%disp(' Sorting data')
sdata=sort(randdata(:));
nn=nx*nz;
% Calculate critical values in randdata to resemble fraction
fraclimit=zeros(length(frac+1),1);
fraclimit(1)=sdata(1);
for n=1:length(frac); fraclimit(n+1)=sdata(round(nn*frac(n))); end
fraclimit(1) = fraclimit(1)- 0.1*abs(fraclimit(1));
fraclimit(length(fraclimit)) = fraclimit(length(fraclimit)) + 0.1*abs(fraclimit(length(fraclimit)));
rdata=zeros(size(randdata));
for n=1:length(frac),
%disp([' Assigning velocities to fraction',num2str(n)])
% OLD SLOW SOLUTION
% [x1,z1]=find(randdata>fraclimit(n) & randdata<fraclimit(n+1) ;
% ixz=find(randdata>fraclimit(n) & randdata<fraclimit(n+1) );
mask = randdata>fraclimit(n) & randdata<=fraclimit(n+1);
rdata = rdata + mask*vel(n);
end
end
rdata=rdata(1:1:nz-ddz,1:1:nx-ddx);
if nargin==0,
figure
imagesc(x,z,rdata)
title(['rseed=',num2str(rseed),', dx=dz=',num2str(dx),', ax=',num2str(ax),', az=',num2str(az),', med=',num2str(med),', pop=',num2str(pop),', Hurst =',num2str(nu),', D=',num2str(D)])
colorbar; colormap(copper)
% rseed,dx,dz,ax,az,ix,iz,pop,med,nu,vel,frac)
end
Below is the error generated:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Input In [124], in <cell line: 16>()
16 if med==3:
17 # von Karman
18 # SQRT of the FOURIER TRANSROMF OF vonk CORR fkt.
19 # disp([' Calculating filter : Chh=vonKarman'])
20 #disp([' Fractal Dimension : ',num2str(D)])
21 #expcorr=((ax*az)/2)./((1+k.^2).^(nu+1)); # (Exact F(C_hh(x))
22 expcorr=1/((1+k**2)**(nu+1)) # (Exact F(C_hh(x))
---> 23 expcorr=expcorr/np.max(np.max(expcorr)) # normalizing sqr(F(C_hh))
24 expcorr=np.sqrt(expcorr)
File <__array_function__ internals>:5, in amax(*args, **kwargs)
File ~\Anaconda3\envs\pyml\lib\site-packages\numpy\core\fromnumeric.py:2754, in amax(a, axis, out, keepdims, initial, where)
2638 @array_function_dispatch(_amax_dispatcher)
2639 def amax(a, axis=None, out=None, keepdims=np._NoValue, initial=np._NoValue,
2640 where=np._NoValue):
2641 """
2642 Return the maximum of an array or maximum along an axis.
2643
(...)
2752 5
2753 """
-> 2754 return _wrapreduction(a, np.maximum, 'max', axis, None, out,
2755 keepdims=keepdims, initial=initial, where=where)
File ~\Anaconda3\envs\pyml\lib\site-packages\numpy\core\fromnumeric.py:86, in _wrapreduction(obj, ufunc, method, axis, dtype, out, **kwargs)
83 else:
84 return reduction(axis=axis, out=out, **passkwargs)
---> 86 return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
ValueError: zero-size array to reduction operation maximum which has no identity
Sources
This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.
Source: Stack Overflow
| Solution | Source |
|---|