Skip to content
Snippets Groups Projects
Commit b7ceb11e authored by Jacobs Matthieu's avatar Jacobs Matthieu
Browse files

Merge branch 'main' of gitlab.mpcdf.mpg.de:mbj/python-scripts

parents 0637ee43 d42083bd
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Matlab Noise Analysis/CorrelatedNoiseSphericalGrid.asv 0 → 100644
+ 119
− 0
View file @ b7ceb11e
% Generate noise with spatial correlation in poloidal, radial and toroidal
% directions
% Sampling on spherical grid with spatial correlation
close all
clear all
clc
Rmax = 500;
rmin = 100;
rstart = 10;
L = 2*pi*Rmax;
% Parameters of the noise distribution
mu = 0; % mean of the noise (spatially constant)
sigma = 1; % standard deviation for the noise (spatially constant)
Lr = rmin/5; % correlation length in radial direction
Lp = 2*pi*rmin/10; % Correlation "length" in poloidal direction
Lt = 2*pi*Rmax/2; % Correlation "length" in toroidal direction
nr = 10; % number of grid points in x-direction
np = 5; % number of grid points in y-direction
nt = 10;
r = (rstart/rmin:0.9/(nr-1):1)*rmin;
p = (0:1/(np-1):1)*2*pi;
[P,R] = meshgrid(r,p);
% Vector of mean at each grid point
Mu = mu*ones(nr*np,1);
% Construct covariance matrices
Sigma_corr = zeros(nr*np);
Sigma_no_corr = zeros(nr*np);
for j = 1:np
for i = 1:nr
k = i + (j-1)*nr;
% Main diagonal: variance
Sigma_no_corr(k,k) = sigma^2;
Sigma_corr(k,k) = sigma^2;
% Off-diagonal: spatial correlation
for n = 1:np
for m = 1:nr
l = m + (n-1)*nr;
if l ~= k
Sigma_corr(k,l) = sigma^2*exp(-sqrt((R(m,n)-R(i,j))^2)/Lr)*exp(-(abs(P(m,n)-P(i,j))*(R(m,n)+R(i,j))/2)/Lp);
end
end
end
end
end
b = issymmetric(Sigma_corr)
d = eig(Sigma_corr)
% Generate the noise from multivariate normal distributions
noise_no_corr = mvnrnd(Mu,Sigma_no_corr);
noise_corr = mvnrnd(Mu,Sigma_corr);
% Reshape onto mesh
noise_no_corr = reshape(noise_no_corr,nr,np);
noise_corr = reshape(noise_corr,nr,np);
% Plot
figure;surf(R,P,noise_no_corr);title('Independent');
figure;surf(R,P,noise_corr); title('Spatial correlation');
% Also check toroidal
t = (0:1/(nt-1):1)*2*pi;
p = (0:1/(np-1):1)*2*pi;
[T,P] = meshgrid(t,p);
% Vector of mean at each grid point
Mu = mu*ones(nt*np,1);
% Construct covariance matrices
Sigma_corr = zeros(nt*np);
Sigma_no_corr = zeros(nt*np);
for j = 1:np
for i = 1:nt
k = i + (j-1)*nt;
% Main diagonal: variance
Sigma_no_corr(k,k) = sigma^2;
Sigma_corr(k,k) = sigma^2;
% Off-diagonal: spatial correlation
for n = 1:np
for m = 1:nt
l = m + (n-1)*np;
if l ~= k
Sigma_corr(k,l) = sigma^2*exp(-(abs(T(m,n)-T(i,j))*(Rmax+R(m,n)*cos(P(m,n))+Rmax+R(i,j)*cos(P(i,j)))/2)/Lt) ...
*exp(-(abs(P(m,n)-P(i,j))*(R(m,n)+R(i,j))/2)/Lp);
end
end
end
end
end
b = issymmetric(Sigma_corr)
d = eig(Sigma_corr)
% Generate the noise from multivariate normal distributions
noise_no_corr = mvnrnd(Mu,Sigma_no_corr);
noise_corr = mvnrnd(Mu,Sigma_corr);
% Reshape onto mesh
noise_no_corr = reshape(noise_no_corr,nt,np);
noise_corr = reshape(noise_corr,nt,np);
% Plot
figure;surf(T,P,noise_no_corr);title('Independent');
figure;surf(T,P,noise_corr); title('Spatial correlation');
Matlab Noise Analysis/CorrelatedNoiseSphericalGrid3D.asv 0 → 100644
+ 77
− 0
View file @ b7ceb11e
% Generate noise with spatial correlation in poloidal, radial and toroidal
% directions
% Sampling on spherical grid with spatial correlation
% 3D case
close all
clear all
clc
Rmax = 500;
rmin = 100;
rstart = 10;
L = 2*pi*Rmax;
% Parameters of the noise distribution
mu = 0; % mean of the noise (spatially constant)
sigma = 1; % standard deviation for the noise (spatially constant)
Lr = rmin/10; % correlation length in radial direction
Lp = 2*pi*rmin/10; % Correlation "length" in poloidal direction
Lt = 2*pi*Rmax/10; % Correlation "length" in toroidal direction
nr = 5;
np = 10;
nt = 2;
r = (rstart/rmin:0.9/(nr-1):1)*rmin;
p = (0:1/(np-1):1)*2*pi;
t = (0:1/(nt-1):1)*2*pi;
[R,P,T] = meshgrid(r,p,t);
% Vector of mean at each grid point
Mu = mu*ones(nr*np*nt,1);
% Construct covariance matrices
Sigma_corr = zeros(nr*np*nt);
Sigma_no_corr = zeros(nr*np*nt);
for h = 1:nt
for j = 1:np
for i = 1:nr
k = i + (j-1)*nr + (h-1)*nr*np;
% Main diagonal: variance
Sigma_no_corr(k,k) = sigma^2;
Sigma_corr(k,k) = sigma^2;
% Off-diagonal: spatial correlation
for p = 1:nt
for n = 1:np
for m = 1:nr
l = m + (n-1)*np + (p-1)*nr*np;
if l ~= k
Sigma_corr(k,l) = sigma^2*exp(-sqrt((R(m,n,p)-R(i,j,h))^2)/Lr)...
*exp(-(abs(P(m,n,p)-P(i,j,h))*(R(m,n,p)+R(i,j,h))/2)/Lp) ...
*exp(-(abs(T(m,n,p)-T(i,j,h))*(Rmax+R(m,n,p)*cos(P(m,n,p))+Rmax+R(i,j,h)*cos(P(i,j,h)))/2)/Lt);
end
end
end
end
end
end
end
b = issymmetric(Sigma_corr)
d = eig(Sigma_corr)
% Generate the noise from multivariate normal distributions
noise_no_corr = mvnrnd(Mu,Sigma_no_corr);
noise_corr = mvnrnd(Mu,Sigma_corr);
% Reshape onto mesh
noise_no_corr = reshape(noise_no_corr,nr,np,nt);
noise_corr = reshape(noise_corr,nr,np,nt);
% Plot
%figure;surf(R,P,T);title('Independent');
%figure;surf(R,P,T); title('Spatial correlation');
\ No newline at end of file
Matlab Noise Analysis/Sparse2DWouter.m 0 → 100644
+ 65
− 0
View file @ b7ceb11e
% Sampling on cartesian grid with spatial correlation
close all
clear all
clc
% Parameters of the noise distribution
mu = 0; % mean of the noise (spatially constant)
sigma = 1; % standard deviation for the noise (spatially constant)
L = 0.25; % correlation length (spatially constant)
Ln =0.65;
% Grid to sample on
Lx = 1; % length of domain in x-direction
Ly = 1; % length of domain in y-direction
nx = 10; % number of grid points in x-direction
ny = 10; % number of grid points in y-direction
x = (0:1/(nx-1):1)*Lx;
y = (0:1/(ny-1):1)*Ly;
[X,Y] = meshgrid(x,y);
% Vector of mean at each grid point
Mu = mu*ones(nx*ny,1);
% Construct covariance matrices
Sigma_corr = zeros(nx*ny);
Sigma_no_corr = zeros(nx*ny);
for j = 1:ny
for i = 1:nx
k = i + (j-1)*ny;
% Main diagonal: variance
Sigma_no_corr(k,k) = sigma^2;
Sigma_corr(k,k) = sigma^2;
% Off-diagonal: spatial correlation
for n = 1:ny
for m = 1:nx
l = m + (n-1)*ny;
if l ~= k
if ((X(m,n)-X(i,j))^2 + (Y(m,n)-Y(i,j))^2)^0.5 < Ln
Sigma_corr(k,l) = sigma^2*exp(-((X(m,n)-X(i,j))^2 + (Y(m,n)-Y(i,j))^2)^0.5/L);
end
end
end
end
end
end
% Generate the noise from multivariate normal distributions
noise_no_corr = mvnrnd(Mu,Sigma_no_corr);
noise_corr = mvnrnd(Mu,Sigma_corr);
% Reshape onto mesh
noise_no_corr = reshape(noise_no_corr,nx,ny);
noise_corr = reshape(noise_corr,nx,ny);
% Plot
figure;surf(X,Y,noise_no_corr);title('Independent');
figure;surf(X,Y,noise_corr); title('Spatial correlation');
\ No newline at end of file
Matlab Noise Analysis/SparseCorr2D.m 0 → 100644
+ 146
− 0
View file @ b7ceb11e
% Generate noise with spatial correlation in poloidal, radial and toroidal
% directions
% Sampling on spherical grid with spatial correlation
close all
clear all
clc
Rmax = 500;
rmin = 100;
rstart = 10;
L = 2*pi*Rmax;
% Parameters of the noise distribution
mu = 0; % mean of the noise (spatially constant)
sigma = 1; % standard deviation for the noise (spatially constant)
Lr = rmin/10; % correlation length in radial direction
Lp = 2*pi*rmin/20; % Correlation "length" in poloidal direction
Lt = 2*pi*Rmax; % Correlation "length" in toroidal direction
%Lp = Lr;
Lcorr = Lr/2;
nr = 100; % number of grid points in x-direction
np = 100; % number of grid points in y-direction
nt = 10;
r = (rstart/rmin:0.9/(nr-1):1)*rmin;
p = (0:1/(np-1):1)*2*pi;
[P,R] = meshgrid(p,r);
% Vector of mean at each grid point
Mu = mu*ones(nr*np,1);
is = (1:np*nr);
js = (1:np*nr);
vs = sigma^2*ones(np*nr,1);
% Construct covariance matrices
Sigma_corr = sparse(is,js,vs,nr*np,nr*np);
%Sigma_corr = zeros(nr*np);
Sigma_no_corr = sparse(is,js,vs,nr*np,nr*np);
for j = 1:np
for i = 1:nr
k = i + (j-1)*nr;
% Main diagonal: variance
% Sigma_no_corr(k,k) = sigma^2;
% Sigma_corr(k,k) = sigma^2;
% Off-diagonal: spatial correlation
nnb = 4;
% for n = j-nnb:j+nnb
% for m = i-nnb:i+nnb
% if (m ==i || n==j) && not(m==i && n==j) %&& not(m==0 || m==nr+1 ||n==0 ||n == np+1)
% mc = m+nr*(m<=0) - nr*(m>=nr+1);
% nc = n+np*(n<=0) - np*(n>=np+1);
% l = mc + (nc-1)*nr;
% Sigma_corr(k,l) = sigma^2*exp(-sqrt((R(mc,nc)-R(i,j))^2)/Lr)*exp(-(abs(P(mc,nc)-P(i,j))*(R(mc,nc)+R(i,j))/2)/Lp);
% end
% end
%end
% Alternatief
for n=1:np
for m=1:nr
l = m+(n-1)*nr;
if l ~= k
distance = ((R(m,n)-R(i,j))^2+ (abs(P(m,n)-P(i,j))*(R(m,n)+R(i,j))/2))^0.5;
if distance < Lcorr
Sigma_corr(k,l) = sigma^2*exp(-sqrt((R(m,n)-R(i,j))^2)/Lr)*exp(-(abs(P(m,n)-P(i,j))*(R(m,n)+R(i,j))/2)/Lp);
%Sigma_corr(k,l) = sigma^2*exp(-distance/Lcorr);
end
end
end
end
end
end
b = issymmetric(Sigma_corr)
%d = eig(Sigma_corr)
[U,S,V] = svd(Sigma_corr);
z = 0+randn(nr*np,1)*sigma;
noise_corr = U *sqrt(S)*z;
% Generate the noise from multivariate normal distributions
noise_no_corr = mvnrnd(Mu,Sigma_no_corr);
%noise_corr = mvnrnd(Mu,Sigma_corr);
% Reshape onto mesh
noise_no_corr = reshape(noise_no_corr,nr,np);
noise_corr = reshape(noise_corr,nr,np);
% Plot
figure;surf(R,P,noise_no_corr);title('Independent');
figure;surf(R,P,noise_corr); title('Spatial correlation');
% Also check toroidal
t = (0:1/(nt-1):1)*2*pi;
p = (0:1/(np-1):1)*2*pi;
[P,T] = meshgrid(p,t);
% Vector of mean at each grid point
Mu = mu*ones(nt*np,1);
% Construct covariance matrices
Sigma_corr = zeros(nt*np);
Sigma_no_corr = zeros(nt*np);
for j = 1:np
for i = 1:nt
k = i + (j-1)*nt;
% Main diagonal: variance
Sigma_no_corr(k,k) = sigma^2;
Sigma_corr(k,k) = sigma^2;
% Off-diagonal: spatial correlation
for n = 1:np
for m = 1:nt
l = m + (n-1)*nt;
if l ~= k
Sigma_corr(k,l) = sigma^2*exp(-(abs(T(m,n)-T(i,j))*(Rmax+R(m,n)*cos(P(m,n))+Rmax+R(i,j)*cos(P(i,j)))/2)/Lt) ...
*exp(-(abs(P(m,n)-P(i,j))*(R(m,n)+R(i,j))/2)/Lp);
end
end
end
end
end
b = issymmetric(Sigma_corr)
d = eig(Sigma_corr)
% Generate the noise from multivariate normal distributions
noise_no_corr = mvnrnd(Mu,Sigma_no_corr);
noise_corr = mvnrnd(Mu,Sigma_corr);
% Reshape onto mesh
noise_no_corr = reshape(noise_no_corr,nt,np);
noise_corr = reshape(noise_corr,nt,np);
% Plot
figure;surf(T,P,noise_no_corr);title('Independent');
figure;surf(T,P,noise_corr); title('Spatial correlation');
Matlab Noise Analysis/SparseCorr3D.m 0 → 100644
+ 105
− 0
View file @ b7ceb11e
% Generate noise with spatial correlation in poloidal, radial and toroidal
% directions
% Sampling on spherical grid with spatial correlation
close all
clear all
clc
Rmax = 500;
rmin = 100;
rstart = 10;
L = 2*pi*Rmax;
% Parameters of the noise distribution
mu = 0; % mean of the noise (spatially constant)
sigma = 1; % standard deviation for the noise (spatially constant)
Lr = rmin/10; % correlation length in radial direction
Lp = 2*pi*rmin/20; % Correlation "length" in poloidal direction
Lt = 2*pi*Rmax; % Correlation "length" in toroidal direction
%Lp = Lr;
Lcorr = Lr/2;
nr = 10; % number of grid points in x-direction
np = 10; % number of grid points in y-direction
nt = 10;
tol = 10^-15;
r = (rstart/rmin:0.9/(nr-1):1)*rmin;
p = (0:1/(np-1):1)*2*pi;
t = (0:1/(nt-1):1)*2*pi;
[P,R,T] = meshgrid(p,r,t);
% Vector of mean at each grid point
Mu = mu*ones(nr*np*nt,1);
is = (1:np*nr*nt);
js = (1:np*nr*nt);
vs = sigma^2*ones(np*nr*nt,1);
% Construct covariance matrices
Sigma_corr = sparse(is,js,vs,nr*np*nt,nr*np*nt);
%Sigma_corr = zeros(nr*np);
Sigma_no_corr = sparse(is,js,vs,nr*np*nt,nr*np*nt);
for h = 1:nt
for j = 1:np
for i = 1:nr
k = i + (j-1)*nr + (h-1)*nr*np;
% Main diagonal: variance
% Sigma_no_corr(k,k) = sigma^2;
% Sigma_corr(k,k) = sigma^2;
% Off-diagonal: spatial correlation
nnb = 4;
% for n = j-nnb:j+nnb
% for m = i-nnb:i+nnb
% if (m ==i || n==j) && not(m==i && n==j) %&& not(m==0 || m==nr+1 ||n==0 ||n == np+1)
% mc = m+nr*(m<=0) - nr*(m>=nr+1);
% nc = n+np*(n<=0) - np*(n>=np+1);
% l = mc + (nc-1)*nr;
% Sigma_corr(k,l) = sigma^2*exp(-sqrt((R(mc,nc)-R(i,j))^2)/Lr)*exp(-(abs(P(mc,nc)-P(i,j))*(R(mc,nc)+R(i,j))/2)/Lp);
% end
% end
%end
% Alternatief
for p = 1:nt
for n=1:np
for m=1:nr
l = m+(n-1)*nr+(p-1)*nr*np;
if l ~= k
distance = ((R(m,n,p)-R(i,j,h))^2+ (abs(P(m,n,p)-P(i,j,h))*(R(m,n,p)+R(i,j,h))/2))^0.5;
if distance < Lcorr
Sigma_corr(k,l) = sigma^2*exp(-sqrt((R(m,n,p)-R(i,j,h))^2)/Lr)...
*exp(-(abs(P(m,n,p)-P(i,j,h))*(R(m,n,p)+R(i,j,h))/2)/Lp)...
*exp(-(abs(T(m,n,p)-T(i,j,h))*(Rmax+R(m,n,p)*cos(P(m,n,p))+Rmax+R(i,j,h)*cos(P(i,j,h)))/2)/Lt);
b=1; % test bp
%Sigma_corr(k,l) = sigma^2*exp(-distance/Lcorr);
end
end
end
end
end
end
end
end
b = issymmetric(Sigma_corr)
%d = eig(Sigma_corr)
[U,S,V] = svds(Sigma_corr,nr*np*nt);
z = 0+randn(nr*np*nt,1)*sigma;
noise_corr = U *sqrt(S)*z;
% Generate the noise from multivariate normal distributions
noise_no_corr = mvnrnd(Mu,Sigma_no_corr);
%noise_corr = mvnrnd(Mu,Sigma_corr);
% Reshape onto mesh
noise_no_corr = reshape(noise_no_corr,nr,np,nt);
noise_corr = reshape(noise_corr,nr,np,nt);
% Plot
%figure;surf(R,P,noise_no_corr);title('Independent');
%figure;surf(R,P,noise_corr); title('Spatial correlation');
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment