Add files via upload

Author: javedali99

Fit_Surrogate Models_MLR_SI_RBF_ and Generate_Results_River WSE.m

@@ -0,0 +1,97 @@
+
+
+%%% This script fit Surrogate Models and Generate results (River Surface water lavels, SWL) from the simulated boundary conditions
+
+%%   1. read the HEC-RAS simulated data 
+%%   2. read the simulated boundary conditions
+%%   3. Fit the models and predict the Surface Water Level
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%% Load HEC-RAS simulated data (SWL, Q, and WL) %% This data will be used for model fitting
+
+%%Load MDA selected dischrage and sea water level and HEC-RAS simulated river water level 
+load('C:\Potomac_conditions_all_ft.mat');
+Inp = [wl_surrogate,Q_surrogate];
+
+%%% HEC-RAS simulations
+load('C:\Data\PotomacRiver\Potomac_sim.mat');
+
+
+cd('C:\Data\PotomacRiver\synthetic_simulations\Simulation data');
+
+FF=dir('*.txt');
+
+%%
+for r=1:length(FF);
+AA=importdata(strcat(FF(r).folder,'\',FF(r).name));
+A=AA.data;
+
+WLs = A(:,6);
+Diss = A(:,5);
+
+Inp2 = [WLs,Diss];
+
+%%
+for i = 1:length(STAT) 
+    disp(strcat('Analysis performing for Simultion'   ,num2str(r),'  XXXXXXXXX', '  Transect No:',   num2str(i)));
+    WL = WS_ELEV(i,:);
+    Res = WL';
+    
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
+    %  MLRM
+    mdl = fitlm(Inp,Res);
+    [m] = predict(mdl,Inp2);
+    preds_tmp(:,i) = m;
+%%	
+%%%%%%%%%% Radial Basis Function %%%%%
+	ZI = rbfinterp([Inp2(:,1)'; Inp2(:,2)'], rbfcreate([Inp(:,1)'; Inp(:,2)'], Res','RBFFunction', 'linear', 'RBFConstant', 0.01,'RBFSmooth', 10));
+    % reshape
+    ZI = reshape(ZI, size(Inp2(:,1)));
+    preds_RBF_tmp(:,i) = ZI;
+	
+%%%%%%%%%%%%%RADIAN BASIS FUNCT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%
+%%%%%%%%%%%%%%%%%  Scatter Interpolant %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    F = scatteredInterpolant(Inp(:,1),Inp(:,2),Res);
+	F.Method = 'natural';% methods
+	WSE_I = F(Inp2(:,1),Inp2(:,2));
+	preds_SI_tmp(:,i) = WSE_I;
+%%%%%%%%%%%%%%%%%  Scatter Interpolant %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+end;
+
+preds_MLR(:,:,r)=preds_tmp;
+preds_RBF(:,:,r)=preds_RBF_tmp;
+preds_SI(:,:,r)=preds_SI_tmp;
+
+Q_simu(:,r)=Inp2(:,2);
+SWL_simu(:,r)=Inp2(:,1);
+
+end;
+
+%%%Save the data
+cd('C:\Surroate Model_Results');
+
+save('surrogate_model_results_potomacRiver.mat','preds_MLR','preds_RBF','preds_SI','Q_simu','SWL_simu','-v7.3,'-nocompression');
+
+
+%%%%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX%%%%%%%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
+         %%XXXXXXXXXXXXX    END   END XXXXXXXXXXXXXXXXXXXXXXXXXXXXX
+%%%%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX%%%%%%%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

rbfcheck.m

@@ -0,0 +1,16 @@
+function maxdiff = rbfcheck(options)
+
+tic;
+
+nodes     = options.('x');
+    y     = options.('y');
+
+s = rbfinterp(nodes, options);
+
+fprintf('RBF Check\n');
+fprintf('max|y - yi| = %e \n', max(abs(s-y)) );
+
+if (strcmp(options.('Stats'),'on'))
+    fprintf('%d points were checked in %e sec\n', length(y), toc);    
+end;
+fprintf('\n');

rbfcreate.m

@@ -0,0 +1,194 @@
+function options = rbfcreate(x, y, varargin)
+%RBFCREATE Creates an RBF interpolation
+%   OPTIONS = RBFSET(X, Y, 'NAME1',VALUE1,'NAME2',VALUE2,...) creates an   
+%   radial base function interpolation 
+%   
+%   RBFCREATE with no input arguments displays all property names and their
+%   possible values.
+%   
+%RBFCREATE PROPERTIES
+% 
+
+%
+% Alex Chirokov, alex.chirokov@gmail.com
+% 16 Feb 2006
+tic;
+% Print out possible values of properties.
+if (nargin == 0) & (nargout == 0)
+  fprintf('               x: [ dim by n matrix of coordinates for the nodes ]\n');
+  fprintf('               y: [   1 by n vector of values at nodes ]\n');
+  fprintf('     RBFFunction: [ gaussian  | thinplate | cubic | multiquadrics | {linear} ]\n');
+  fprintf('     RBFConstant: [ positive scalar     ]\n');
+  fprintf('       RBFSmooth: [ positive scalar {0} ]\n');
+  fprintf('           Stats: [ on | {off} ]\n');
+  fprintf('\n');
+  return;
+end
+Names = [
+    'RBFFunction      '
+    'RBFConstant      '
+    'RBFSmooth        '
+    'Stats            '
+];
+[m,n] = size(Names);
+names = lower(Names);
+
+options = [];
+for j = 1:m
+  options.(deblank(Names(j,:))) = [];
+end
+
+%**************************************************************************
+%Check input arrays 
+%**************************************************************************
+[nXDim nXCount]=size(x);
+[nYDim nYCount]=size(y);
+
+if (nXCount~=nYCount)
+  error(sprintf('x and y should have the same number of rows'));
+end;
+
+if (nYDim~=1)
+  error(sprintf('y should be n by 1 vector'));
+end;
+
+options.('x')           = x;
+options.('y')           = y;
+%**************************************************************************
+%Default values 
+%**************************************************************************
+options.('RBFFunction') = 'linear';
+options.('RBFConstant') = (prod(max(x')-min(x'))/nXCount)^(1/nXDim); %approx. average distance between the nodes 
+options.('RBFSmooth')   = 0;
+options.('Stats')       = 'off';
+
+%**************************************************************************
+% Argument parsing code: similar to ODESET.m
+%**************************************************************************
+
+i = 1;
+% A finite state machine to parse name-value pairs.
+if rem(nargin-2,2) ~= 0
+  error('Arguments must occur in name-value pairs.');
+end
+expectval = 0;                          % start expecting a name, not a value
+while i <= nargin-2
+  arg = varargin{i};
+    
+  if ~expectval
+    if ~isstr(arg)
+      error(sprintf('Expected argument %d to be a string property name.', i));
+    end
+    
+    lowArg = lower(arg);
+    j = strmatch(lowArg,names);
+    if isempty(j)                       % if no matches
+      error(sprintf('Unrecognized property name ''%s''.', arg));
+    elseif length(j) > 1                % if more than one match
+      % Check for any exact matches (in case any names are subsets of others)
+      k = strmatch(lowArg,names,'exact');
+      if length(k) == 1
+        j = k;
+      else
+        msg = sprintf('Ambiguous property name ''%s'' ', arg);
+        msg = [msg '(' deblank(Names(j(1),:))];
+        for k = j(2:length(j))'
+          msg = [msg ', ' deblank(Names(k,:))];
+        end
+        msg = sprintf('%s).', msg);
+        error(msg);
+      end
+    end
+    expectval = 1;                      % we expect a value next
+    
+  else
+    options.(deblank(Names(j,:))) = arg;
+    expectval = 0;      
+  end
+  i = i + 1;
+end
+
+if expectval
+  error(sprintf('Expected value for property ''%s''.', arg));
+end
+
+    
+%**************************************************************************
+% Creating RBF Interpolatin
+%**************************************************************************
+
+switch lower(options.('RBFFunction'))
+      case 'linear'          
+        options.('rbfphi')   = @rbfphi_linear;
+      case 'cubic'
+        options.('rbfphi')   = @rbfphi_cubic;
+      case 'multiquadric'
+        options.('rbfphi')   = @rbfphi_multiquadrics;
+      case 'thinplate'
+        options.('rbfphi')   = @rbfphi_thinplate;
+      case 'gaussian'
+        options.('rbfphi')   = @rbfphi_gaussian;
+    otherwise
+        options.('rbfphi')   = @rbfphi_linear;
+end
+
+phi       = options.('rbfphi');
+
+A=rbfAssemble(x, phi, options.('RBFConstant'), options.('RBFSmooth'));
+
+b=[y'; zeros(nXDim+1, 1)];                       
+
+%inverse
+rbfcoeff=A\b;
+
+%SVD
+% [U,S,V] = svd(A);
+% 
+% for i=1:1:nXCount+1
+%     if (S(i,i)>0) S(i,i)=1/S(i,i); end;   
+% end;    
+% rbfcoeff = V*S'*U*b;
+
+
+options.('rbfcoeff') = rbfcoeff;
+
+
+if (strcmp(options.('Stats'),'on'))
+    fprintf('%d point RBF interpolation was created in %e sec\n', length(y), toc);  
+    fprintf('\n');
+end;
+
+function [A]=rbfAssemble(x, phi, const, smooth)
+[dim n]=size(x);
+A=zeros(n,n);
+for i=1:n
+    for j=1:i
+        r=norm(x(:,i)-x(:,j));
+        temp=feval(phi,r, const);
+        A(i,j)=temp;
+        A(j,i)=temp;
+    end
+    A(i,i) = A(i,i) - smooth;
+end
+% Polynomial part
+P=[ones(n,1) x'];
+A = [ A      P
+      P' zeros(dim+1,dim+1)];
+
+%**************************************************************************
+% Radial Base Functions
+%************************************************************************** 
+function u=rbfphi_linear(r, const)
+u=r;
+
+function u=rbfphi_cubic(r, const)
+u=r.*r.*r;
+
+function u=rbfphi_gaussian(r, const)
+u=exp(-0.5*r.*r/(const*const));
+
+function u=rbfphi_multiquadrics(r, const)
+u=sqrt(1+r.*r/(const*const));
+
+function u=rbfphi_thinplate(r, const)
+u=r.*r.*log(r+1);

rbfinterp.m

@@ -0,0 +1,37 @@
+function [f] = rbfinterp(x, options)
+tic;
+phi       = options.('rbfphi');
+rbfconst  = options.('RBFConstant');
+nodes     = options.('x');
+rbfcoeff  = (options.('rbfcoeff'))';
+
+
+[dim              n] = size(nodes);
+[dimPoints  nPoints] = size(x);
+
+if (dim~=dimPoints)
+  error(sprintf('x should have the same number of rows as an array used to create RBF interpolation'));
+end;
+
+f = zeros(1, nPoints);
+r = zeros(1, n);
+
+for i=1:1:nPoints
+	s=0;
+    r =  (x(:,i)*ones(1,n)) - nodes;
+    r = sqrt(sum(r.*r, 1));
+%     for j=1:n
+%          r(j) =  norm(x(:,i) - nodes(:,j));
+%     end
+    
+     s = rbfcoeff(n+1) + sum(rbfcoeff(1:n).*feval(phi, r, rbfconst));
+ 
+	for k=1:dim
+       s=s+rbfcoeff(k+n+1)*x(k,i);     % linear part
+	end
+	f(i) = s;
+end;
+
+if (strcmp(options.('Stats'),'on'))
+    fprintf('Interpolation at %d points was computed in %e sec\n', length(f), toc);    
+end;