diff --git a/README.md b/README.md index c659c44..cd7cded 100644 --- a/README.md +++ b/README.md @@ -1,16 +1,6 @@ ---- TO BE REMOVED --- +# MIDSTREAM: Model Testing and Forecasting — Official Application Repository -This is a template for an official application repository. After creating a new repository, fill in the placeholder fields and remove this text. - -✅ Follow [Section 1: Migrating an Application to Git repository](https://docs.cyfronet.pl/x/W4MdDg) when migrating an existing application from `apps_to_migrate`. - -✅ Follow [Section 2: Managing Code Changes in Official Application Repositories](https://docs.cyfronet.pl/x/W4MdDg) to develop or apply updates. - ---- END TO BE REMOVED --- - -# — Official Application Repository - -This repository contains the source code and configuration files for the `` application used within the [EPISODES Platform](https://EpisodesPlatform.eu/). +This repository contains the source code and configuration files for the `MIDSTREAM: Model Testing and Forecasting` application used within the [EPISODES Platform](https://EpisodesPlatform.eu/). 📦 To test or modify this application in the EPISODES Platform environment, follow the guide: https://docs.cyfronet.pl/display/ISDOC/Editing+application+codes+guide diff --git a/src/calcrate_exponential.m b/src/calcrate_exponential.m new file mode 100644 index 0000000..c7d13ee --- /dev/null +++ b/src/calcrate_exponential.m @@ -0,0 +1,66 @@ +function[TEST_seism_time, TEST_seism_Nev, TEST_IR_time_resampled, ... + TEST_IR_data_resampled, Forecast_cumulated_Nev, CumSeis, Nsmpls] = calcrate_exponential(... + Parameter_Samples, TEST_IR_time, TEST_IR_data, TEST_Tev_i, HighPrctile, LowPrctile, ... + Dt_fcst) + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Resample the IR data for considering the Dt for forecasts: + +Ttmp = TEST_IR_time(1); +j=1; +while Ttmp + Dt_fcst <= TEST_IR_time(end) + ind_ir_sel = find(and(TEST_IR_time>=Ttmp, TEST_IR_time=Ttmp, TEST_Tev_i0 + TEST_IR_data_resampled(j,1) = max(TEST_IR_data(ind_ir_sel)); + else + TEST_IR_data_resampled(j,1) = 0; + end + + TEST_seism_time(j,1) = TEST_IR_time_resampled(j,1); %mean(TEST_Tev_i(ind_seism)); + TEST_seism_Nev(j,1) = length(ind_seism); + j=j+1; + Ttmp = Ttmp + Dt_fcst; + clear ind_ir_sel ind_seism +end + +%Find IR values = 0 and set them as NaN: +clear ind +ind = find(TEST_IR_data_resampled<1e-6); +TEST_IR_data_resampled(ind) = NaN; %min(TEST_IR_data_resampled); + +[Nsmpls, Npars] = size(Parameter_Samples); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Calculate seism. rtaes from IR data using the Power Law model +clear SeismRate +for i=1:Nsmpls + %seism rate, Nev/day, rescaled to the forecasting window Dt_fcst: + %SeismRate(:,i) = Dt_fcst.* ( (1./(10.^((Parameter_Samples(i,2).*(TEST_IR_data_resampled)) + Parameter_Samples(i,1)))) ); % Nev/Dt_fcst + %indtmp = isnan(SeismRate(:,i)); + %SeismRate(indtmp,i) = 0; + %clear indtmp + %CumSeis(:,i) = cumsum(SeismRate(:,i)); + MU(:,i) = Parameter_Samples(i,1).*(10.^((TEST_IR_data_resampled).*(Parameter_Samples(i,2)))); % days + indtmp = isnan(MU(:,i)); + MU(indtmp,i) = median(MU(~isnan(MU(:,i)),i)); + %SeismRate(:,i) = Dt_fcst.*(1./MU(:,i)); + for j=1:numel(TEST_IR_data_resampled) + SeismRate(j,i) = Dt_fcst.*(1/median(exprnd(MU(j,i), 100,1))); + end + clear indtmp mu_tmp + CumSeis(:,i) = cumsum(SeismRate(:,i)); + +end + +for i=1:numel(TEST_IR_time_resampled) + Forecast_cumulated_Nev(i,:) = [prctile(CumSeis(i,:), LowPrctile) median(CumSeis(i,:)) prctile(CumSeis(i,:), HighPrctile)]; +end + + + + + + diff --git a/src/calcrate_powerlaw.m b/src/calcrate_powerlaw.m new file mode 100644 index 0000000..a87e5b2 --- /dev/null +++ b/src/calcrate_powerlaw.m @@ -0,0 +1,71 @@ +function[TEST_seism_time, TEST_seism_Nev, TEST_IR_time_resampled, ... + TEST_IR_data_resampled, Forecast_cumulated_Nev, CumSeis, Nsmpls] = calcrate_powerlaw(... + Parameter_Samples, TEST_IR_time, TEST_IR_data, TEST_Tev_i, HighPrctile, LowPrctile, ... + Dt_fcst) + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Resample the IR data for considering the Dt for forecasts: + +Ttmp = TEST_IR_time(1); +j=1; +while Ttmp + Dt_fcst <= TEST_IR_time(end) + ind_ir_sel = find(and(TEST_IR_time>=Ttmp, TEST_IR_time=Ttmp, TEST_Tev_i0 + TEST_IR_data_resampled(j,1) = max(TEST_IR_data(ind_ir_sel)); + else + TEST_IR_data_resampled(j,1) = 0; + end + + TEST_seism_time(j,1) = TEST_IR_time_resampled(j,1); %mean(TEST_Tev_i(ind_seism)); + TEST_seism_Nev(j,1) = length(ind_seism); + j=j+1; + Ttmp = Ttmp + Dt_fcst; + clear ind_ir_sel ind_seism +end + +%Find IR values = 0 and set them as NaN: +clear ind +ind = find(TEST_IR_data_resampled<1e-6); +TEST_IR_data_resampled(ind) = NaN; %min(TEST_IR_data_resampled); + +[Nsmpls, Npars] = size(Parameter_Samples); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Calculate seism. rtaes from IR data using the Power Law model +clear SeismRate +% for i=1:Nsmpls +% %seism rate, Nev/day, rescaled to the forecasting window Dt_fcst: +% SeismRate(:,i) = Dt_fcst.* ( (1./(10.^((Parameter_Samples(i,2).*log10(TEST_IR_data_resampled)) + Parameter_Samples(i,1)))) ); % Nev/Dt_fcst +% indtmp = isnan(SeismRate(:,i)); +% SeismRate(indtmp,i) = 0; +% clear indtmp +% CumSeis(:,i) = cumsum(SeismRate(:,i)); +% end + +for i=1:Nsmpls + %seism rate, Nev/day, rescaled to the forecasting window Dt_fcst: + %SeismRate(:,i) = Dt_fcst.* ( (1./(10.^((Parameter_Samples(i,2).*log10(TEST_IR_data_resampled)) + Parameter_Samples(i,1)))) ); % Nev/Dt_fcst + %MU(:,i) = 10.^((Parameter_Samples(i,2).*log10(TEST_IR_data_resampled)) + Parameter_Samples(i,1)); % days + MU(:,i) = Parameter_Samples(i,1).*((TEST_IR_data_resampled).^(Parameter_Samples(i,2))); % days + indtmp = isnan(MU(:,i)); + MU(indtmp,i) = median(MU(~isnan(MU(:,i)),i)); + %SeismRate(:,i) = Dt_fcst.*(1./MU(:,i)); + for j=1:numel(TEST_IR_data_resampled) + SeismRate(j,i) = Dt_fcst.*(1/median(exprnd(MU(j,i), 100,1))); + end + clear indtmp mu_tmp + CumSeis(:,i) = cumsum(SeismRate(:,i)); +end + +for i=1:numel(TEST_IR_time_resampled) + Forecast_cumulated_Nev(i,:) = [prctile(CumSeis(i,:), LowPrctile) median(CumSeis(i,:)) prctile(CumSeis(i,:), HighPrctile)]; +end + + + + + + diff --git a/src/calcrate_stationary.m b/src/calcrate_stationary.m new file mode 100644 index 0000000..d020533 --- /dev/null +++ b/src/calcrate_stationary.m @@ -0,0 +1,55 @@ +function[TEST_seism_time, TEST_seism_Nev, TEST_IR_time_resampled, ... + TEST_IR_data_resampled, Forecast_cumulated_Nev, CumSeis, Nsmpls] = calcrate_stationary(... + Parameter_Samples, TEST_IR_time, TEST_IR_data, TEST_Tev_i, HighPrctile, LowPrctile, ... + Dt_fcst) + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Resample the IR data for considering the Dt for forecasts: + +Ttmp = TEST_IR_time(1); +j=1; +while Ttmp + Dt_fcst <= TEST_IR_time(end) + ind_ir_sel = find(and(TEST_IR_time>=Ttmp, TEST_IR_time=Ttmp, TEST_Tev_i0 + TEST_IR_data_resampled(j,1) = max(TEST_IR_data(ind_ir_sel)); + else + TEST_IR_data_resampled(j,1) = 0; + end + TEST_seism_time(j,1) = TEST_IR_time_resampled(j,1); %mean(TEST_Tev_i(ind_seism)); + TEST_seism_Nev(j,1) = length(ind_seism); + j=j+1; + Ttmp = Ttmp + Dt_fcst; + clear ind_ir_sel ind_seism +end + +%Find IR values = 0 and set them as NaN: +clear ind +ind = find(TEST_IR_data_resampled<1e-6); +TEST_IR_data_resampled(ind) = NaN; %min(TEST_IR_data_resampled); +[Nsmpls, Npars] = size(Parameter_Samples); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Calculate seism. rtaes from IR data using the Power Law model +clear SeismRate +for i=1:Nsmpls + %seism rate, Nev/day, rescaled to the forecasting window Dt_fcst: + SeismRate(:,i) = Dt_fcst.* ( (1./(Parameter_Samples(i,1).*(ones(numel(TEST_IR_data_resampled),1))))); + indtmp = isnan(SeismRate(:,i)); + SeismRate(indtmp,i) = median(SeismRate(~isnan(SeismRate(:,i)),i)); + %SeismRate(indtmp,i) = 0; + clear indtmp + CumSeis(:,i) = cumsum(SeismRate(:,i)); +end + +for i=1:numel(TEST_IR_time_resampled) + Forecast_cumulated_Nev(i,:) = [prctile(CumSeis(i,:), LowPrctile) median(CumSeis(i,:)) prctile(CumSeis(i,:), HighPrctile)]; +end + + + + + + diff --git a/src/midstream_model_testing.m b/src/midstream_model_testing.m new file mode 100644 index 0000000..5fa15fd --- /dev/null +++ b/src/midstream_model_testing.m @@ -0,0 +1,244 @@ +function[WarningMessages, Exit_ID] = midstream_model_testing(Training_Model_Samples, ... + Testing_dataset_for_testing, model, Model_Selection, ir_units, ... + iet_units, ir_threshold, iet_threshold, Dt_fcst, FcastTime_window, ... + Path_output_results, HighPrctile, LowPrctile, ScaleNumEventForecast, ... + Script_Mode) + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Load data +disp(['Reading model training parameter samples (model: ' model.folder{Model_Selection} ')']) +Parameter_Samples = csvread(Training_Model_Samples, 1); + +disp(['Reading input data saved for testing: ' Testing_dataset_for_testing]) +% +if Script_Mode.ID == 1 + WarningMessages = {'Using the tool for testing the forecasting performance with the testing dataset',... + 'List of warning messages:'}; + disp(cell2mat(WarningMessages(1))) + load(Testing_dataset_for_testing); +else + WarningMessages = {'Using the tool for forecasting seismicity using a custom timeseries of Injection rate data', ... + 'Error:', ... + 'This functionality is not yet available! Program Exit without results'}; + disp(cell2mat(WarningMessages(1))) + return +end +% It loads: +% CatMagnitude 1x2 4 char +% DelayTime 1x1 8 double +% Mc 1x1 8 double +% TEST_IETij 93x1 744 double +% TEST_IR_data_raw 37077x1 296616 double +% TEST_IR_time_raw 37077x1 296616 double +% TEST_IRij 93x1 744 double +% TEST_Magev_i 94x1 752 double +% TEST_Tev_i 94x1 752 double + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% check if events in the test dataset are above Mc: +%WarningMessages = {'List of warning messages:'}; + +ind_mc = find(TEST_Magev_i0 + WarningMessages(numel(WarningMessages)+1) = {'ERROR: events below Mc have been found. Inconsistency in input data. Check'}; +end +clear ind_mc + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% select data above min. IR & IET defined thresholds +ind = find(and(TEST_IRij>ir_threshold, TEST_IETij>iet_threshold)); +%Keep source data in these variables (useful for plots) +IR = TEST_IRij(ind); +IET = TEST_IETij(ind); +clear TEST_IRij TEST_IETij +TEST_IRij = IR; +TEST_IETij = IET; +clear IR IET + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% % Check if the sampling rate of IR data and observed seimsicity rates are +% coherent with the Dt time window selected for Seism. rate forecasts. +% Eventually write a warning message + +Mean_iet = median(TEST_IETij); +Mean_dt_ir = median(diff(TEST_IR_time_raw)); + +disp(['Dt - seism. rate forecast (' num2str(Dt_fcst) ' days; ' num2str(Dt_fcst*24) ' hours; ' num2str(Dt_fcst*24*60) ' seconds)']) +disp(['Median IET seismicity for testing (' num2str(Mean_iet) ' days; ' num2str(Mean_iet*24) ' hours; ' num2str(Mean_iet*24*60) ' seconds)']) +disp(['Median sampling rate of IR data (' num2str(Mean_dt_ir) ' days; ' num2str(Mean_dt_ir*24) ' hours; ' num2str(Mean_dt_ir*24*60) ' seconds)']) + +if Dt_fcst < Mean_iet + WarningMessages(numel(WarningMessages)+1) = {'Warning: Dt for forecast is lower than the median IET. Increase Dt.'}; +elseif Dt_fcst < Mean_dt_ir + WarningMessages(numel(WarningMessages)+1) = {'Warning: Dt for forecast is lower than the median Sampling rate if IR data. Increase Dt.'}; +end + +%start forecast from the time of the firat recorded event in the testing +%period +To = min(TEST_Tev_i); +ind_ir = find(TEST_IR_time_raw >= To); + +TEST_IR_time = TEST_IR_time_raw(ind_ir); +TEST_IR_data = TEST_IR_data_raw(ind_ir); + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Seism rate calculation for differnet models: + + +if Model_Selection== 1 %stationary + disp('Testing the stationary model...') + [TEST_seism_time, TEST_seism_Nev, TEST_IR_time_resampled, ... + TEST_IR_data_resampled, Forecast_cumulated_Nev, ... + CumSeis, Nsmpls] = calcrate_stationary(Parameter_Samples, TEST_IR_time, ... + TEST_IR_data, TEST_Tev_i, HighPrctile, LowPrctile, Dt_fcst); + +elseif Model_Selection== 2 % Power Law + disp('Testing the Power-law model...') + [TEST_seism_time, TEST_seism_Nev, TEST_IR_time_resampled, ... + TEST_IR_data_resampled, Forecast_cumulated_Nev, ... + CumSeis, Nsmpls] = calcrate_powerlaw(Parameter_Samples, TEST_IR_time, ... + TEST_IR_data, TEST_Tev_i, HighPrctile, LowPrctile, Dt_fcst); + +elseif Model_Selection== 3 % Exponential + disp('Testing the Exponential model...') + [TEST_seism_time, TEST_seism_Nev, TEST_IR_time_resampled, ... + TEST_IR_data_resampled, Forecast_cumulated_Nev, ... + CumSeis, Nsmpls] = calcrate_exponential(Parameter_Samples, TEST_IR_time, ... + TEST_IR_data, TEST_Tev_i, HighPrctile, LowPrctile, Dt_fcst); +end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Plot results, save output data and plots + +%% +%Create folder +path2SaveOutData = [Path_output_results model.folder{Model_Selection} '/Forecast_Nev_per_' cell2mat(FcastTime_window)]; +try + mkdir(path2SaveOutData) +catch + disp('Forder for output data already exists') +end + +%% +figure(1) +set(1, 'Position', [943 264 964 670]) + +% Times normalized respect to the injection start: +To=TEST_IR_time(1); + +for i=1:Nsmpls + if i==1 + plt1 = plot(TEST_IR_time_resampled - To, 1+CumSeis(:,i), '.', 'Color', [0.8 0.8 0.8]); + else + plot(TEST_IR_time_resampled - To, 1+CumSeis(:,i), '.', 'Color', [0.8 0.8 0.8]); + end + hold on + %pause +end + +plt2 = plot(TEST_seism_time - To, cumsum(TEST_seism_Nev), '.r', 'LineWidth', 1.5); +plt3 = plot(TEST_IR_time_resampled - To, 1+Forecast_cumulated_Nev(:,2), '.-', 'Color', 'k', 'LineWidth', 1.5); +plt4 = plot(TEST_IR_time_resampled - To, 1+Forecast_cumulated_Nev(:,1), '.-', 'Color', 'b', 'LineWidth', 1.); +plot(TEST_IR_time_resampled - To, 1+Forecast_cumulated_Nev(:,3), '.-', 'Color', 'b', 'LineWidth', 1.); +hold off +grid on +legend([plt2 plt3 plt4 plt1], 'Cumulated number of events (observed in the testing period)', ... + 'Cumulated n. events (median forecast)', ... + ['Cumulated n. events (' num2str(LowPrctile) 'th - ' num2str(HighPrctile) 'th prctiles)'], ... + 'Cumulated number of events (samples of model parameters)', 'Location', 'southeast'); +try + set(gca, 'YScale', ScaleNumEventForecast) +catch + %default: log plot + set(gca, 'YScale', 'log') +end + +xlim([min(TEST_IR_time_resampled)-To max(TEST_IR_time_resampled)-To]) +xlabel(['Time (days) after ' datestr(To)]) +ylabel('Cumulated number of events') +title({'Comparison between observed and forecast seimsicity', ... + ['(cumulated number of events considering event rates / ' cell2mat(FcastTime_window) ')']}) + + +%% +figure(2) +set(2, 'Position', [1321 107 588 846]) +subplot(3,1,1) +plot(TEST_Tev_i-To, TEST_Magev_i, 'o', 'MarkerSize',6 , 'MarkerEdgeColor','k', 'MarkerFaceColor', [0.8 0.8 0.8]) +hold on +plot([min(TEST_Tev_i)-To max(TEST_Tev_i)-To], [Mc Mc], '--r', 'LineWidth', 1.5) +xlim([min(TEST_Tev_i)-To max(TEST_Tev_i)-To]) +grid on +xlabel(['Time (days) after ' datestr(To)]) +ylabel('Magnitude') +legend('Seismic event', 'Mc', 'Location', 'northeast') +CumNum_observed_events = cumsum(ones(length(TEST_Tev_i),1)); +% disp(['Time injection start: ' datestr(min(TEST_IR_time_raw))]) +% disp(['Time first event: ' datestr(min(TEST_Tev_i))]) +% Tau = min(TEST_Tev_i) -min(TEST_IR_time_raw); +% disp(['\tau (Tev-Tinj) = ' num2str(Tau) ' days']) + +subplot(3,1,2) +colororder({'k','k'}) +yyaxis left +plot(TEST_Tev_i-To, CumNum_observed_events, 'k', 'LineWidth', 1.5) +ylabel({'Cumulated number of recorded', ['events (above Mc= ' num2str(Mc) ')']}) +grid on +yyaxis right +plot(TEST_IR_time_raw(ind_ir)-To, TEST_IR_data_raw(ind_ir), 'r', 'LineWidth', 1.5) +xlim([min(TEST_Tev_i)-To max(TEST_Tev_i)-To]) +ylabel(['Injection rate, ' ir_units]) +xlabel(['Time (days) after ' datestr(To)]) +legend('cumulated number of events', 'injection rate', 'Location', 'northwest') +title('Original data') + +subplot(3,1,3) +colororder({'k','k'}) +yyaxis left +plot(TEST_seism_time, cumsum(TEST_seism_Nev), 'k', 'LineWidth', 1.5) +xlabel('time') +ylabel({'Cumulated number of recorded', ['events (above Mc= ' num2str(Mc) ')']}) +grid on +yyaxis right +plot(TEST_IR_time_resampled, TEST_IR_data_resampled, 'r', 'LineWidth', 1.5) +ylabel(['Injection rate, ' ir_units]) +legend('cumulated number of events', 'injection rate') +title(['Data averaged in time windows of 1 ' cell2mat(FcastTime_window)]) + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% save plot +fout1 = [path2SaveOutData '/' model.folder{Model_Selection} '_SeismicityRateForecast_comparison_Nev_per_' cell2mat(FcastTime_window) '.jpg']; +saveas(1,fout1,'jpg') +fout2 = [path2SaveOutData '/' model.folder{Model_Selection} '_ProcessedData_Nev_per_' cell2mat(FcastTime_window) '.jpg']; +saveas(2,fout2,'jpg') + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Save output parameter samples +NomOutSolution_model = [path2SaveOutData '/' model.folder{Model_Selection} '_Seismicity_rate_forecast_Nev_per_' cell2mat(FcastTime_window) '.txt']; +cHeader = {['year'] ['month'] ['day'] ['hour'] ['minute'] ['second'] ... + ['lenght_window(days)'] ['ObservedSeismicity(cumulated)'] ... + ['SeismicityForecast(cumulated)-' num2str(LowPrctile) 'prctile'] ... + ['SeismicityForecast(cumulated)-median'] ... + ['SeismicityForecast(cumulated)-' num2str(HighPrctile) 'prctile']}; +commaHeader = [cHeader;repmat({','},1,numel(cHeader))]; +commaHeader = commaHeader(:)'; +textHeader = cell2mat(commaHeader); %cHeader in text with commas +%write header to file +fid = fopen(NomOutSolution_model,'w'); +fprintf(fid,'%s\n',textHeader); +fclose(fid); +%write data to end of file +%Data2Write = [datevec(TEST_IR_time_resampled) repmat(Dt_fcst,numel(TEST_IR_time_resampled),1) cumsum(TEST_seism_Nev) Forecast_cumulated_Nev]; +dlmwrite(NomOutSolution_model,[datevec(TEST_IR_time_resampled) repmat(Dt_fcst,numel(TEST_IR_time_resampled),1) cumsum(TEST_seism_Nev) Forecast_cumulated_Nev],'-append'); + + + + +