diff --git a/ComparisonWeakSINDy/Main_NoiseCompare01.m b/ComparisonWeakSINDy/Main_NoiseCompare01.m new file mode 100644 index 0000000..60d6e9e --- /dev/null +++ b/ComparisonWeakSINDy/Main_NoiseCompare01.m @@ -0,0 +1,106 @@ +%% This file will swipe the noise level +clc;clear all;close all; +%% First load the noise and simulation data +load('GeneratedNoise.mat') +addpath("Result") +addpath("Functions") +%% Define the parameters for the WeakSINDy_AG +% Library order +polyorder=2; +usesine=0; + +% Optimization parameter +gamma=0; +tau = 1; + +% Number of test function +K = 200; + +% Test function order +poly = 14; + +% Test function support size +s = 31; +r_whm = 8; + +% Parameters for the WeakSINDy +wsindy_params = {s, K, poly, tau}; + +% Define the thresholding parameters +lam=linspace(0.01,0.95,30); +lamNum=length(lam); + +% Define a matrix to store the values +Theta_test=poolData(x_test,stateVar,polyorder,usesine); +Xi0=zeros(lamNum,size(Theta_test,2),stateVar); +TestingScore=zeros(lamNum,1); +SuccessOrNot_WS=zeros(N_run,NoiseNum); +Xi_WS=zeros(lamNum,size(Theta_test,2),stateVar); + +% Store the parameter error +Eparm_WS=zeros(N_run,NoiseNum); + +Xi_base=[-10,28,0; + 10,-1,0; + 0,0,-8/3; + 0,0,0; + 0,0,1; + 0,-1,0; + 0,0,0; + 0,0,0; + 0,0,0]; +%% Next run a for loop to test each noise level +pin=0; + +for i=1:N_run + fprintf(strcat("This is the run:",string(i),"\n")) + for j=1:NoiseNum + fprintf(strcat("\t","Using noise level as: ",string(NoisePercentageToSwipe(j)),"\n")) + pin=pin+1; + % Get the observation data + xes=squeeze(xes_List(pin,:,:)); + Theta=poolData(xes,stateVar,polyorder,usesine); + + % Swipe the thresholding parameter + parfor k=1:lamNum + fprintf(strcat("\t\t","Swiping lambda, using lambda as: ",string(lam(k)),"\n")) + + % Get the result for current lam + Xi0(k,:,:)=WeakSINDy_AG(t,xes,Theta,stateVar,lam(k),gamma,r_whm,wsindy_params); + + % Calculate the error on the test data + TestingScore(k,1)=norm(dx_test-Theta_test*squeeze(Xi0(k,:,:)))/norm(dx_test); + + end + + % Test whether the structure of the model is correct + for k=1:lamNum + if norm((Xi_base~=0)-(squeeze(Xi0(k,:,:))~=0))==0 + SuccessOrNot_WS(i,j)=1; + Xi_dummy=squeeze(Xi0(k,:,:)); + Xi_WS(pin,:,:)=Xi_dummy; + end + end + + % If no model has the correct structure, choose the model that has + % lowest testing error as final model + if SuccessOrNot_WS(i,j)==0 + % Now select the best Xi + [minVal,index]=min(TestingScore); + + % Store the current result + Xi_dummy=squeeze(Xi0(index,:,:)); + Xi_WS(pin,:,:)=Xi_dummy; + end + + % Store the parameter error + Eparm_WS(i,j)=norm(Xi_base-Xi_dummy)/norm(Xi_base); + + end +end + +%% Save the results +Ep_WS_01=Eparm_WS; +SuccessOrNot_WS_01=SuccessOrNot_WS; + +%save ("Result\W_SINDy_dt_01.mat","Ep_WS_01","SuccessOrNot_WS_01")