From c1c53ba7dbffab14cb59bc8dd4625cb87be658f9 Mon Sep 17 00:00:00 2001 From: Silas Date: Mon, 7 Oct 2024 09:37:55 +0200 Subject: [PATCH] Add Datastorage --- Classes/Warehouse_class/classes/DataStorage.m | 224 ++++++++++ Classes/Warehouse_class/classes/Parameter.m | 44 ++ .../functions/fwm_plots/CompleteRoutine.m | 251 +++++++++++ .../CompleteRoutine_DifferentChannels_Fig3.m | 250 +++++++++++ .../functions/fwm_plots/automate_JLT_plots.m | 394 ++++++++++++++++++ .../functions/fwm_plots/dispersion_only.m | 83 ++++ .../dispersion_validation_miniskript.m | 52 +++ .../functions/fwm_plots/generatePlots.m | 61 +++ .../functions/fwm_plots/plot3dCurve.m | 248 +++++++++++ .../functions/fwm_plots/plotBerVsZDW.m | 300 +++++++++++++ .../fwm_plots/plotBerVsZdwFailureRate.m | 157 +++++++ .../fwm_plots/plotChannelSpacingAna.m | 51 +++ .../functions/fwm_plots/plotCurve.m | 278 ++++++++++++ .../functions/fwm_plots/plotHistogram.m | 151 +++++++ .../functions/fwm_plots/plotViolin.m | 206 +++++++++ .../fwm_plots/plot_ber_distribution.m | 58 +++ .../Warehouse_class/functions/interpCurve.m | 9 + .../functions/phase_predist_plots/build_wh.m | 129 ++++++ Classes/Warehouse_class/plotjob_example.mat | Bin 0 -> 672 bytes 19 files changed, 2946 insertions(+) create mode 100644 Classes/Warehouse_class/classes/DataStorage.m create mode 100644 Classes/Warehouse_class/classes/Parameter.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/CompleteRoutine.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/CompleteRoutine_DifferentChannels_Fig3.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/automate_JLT_plots.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/dispersion_only.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/dispersion_validation_miniskript.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/generatePlots.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/plot3dCurve.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/plotBerVsZDW.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/plotBerVsZdwFailureRate.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/plotChannelSpacingAna.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/plotCurve.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/plotHistogram.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/plotViolin.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/plot_ber_distribution.m create mode 100644 Classes/Warehouse_class/functions/interpCurve.m create mode 100644 Classes/Warehouse_class/functions/phase_predist_plots/build_wh.m create mode 100644 Classes/Warehouse_class/plotjob_example.mat diff --git a/Classes/Warehouse_class/classes/DataStorage.m b/Classes/Warehouse_class/classes/DataStorage.m new file mode 100644 index 0000000..5111ade --- /dev/null +++ b/Classes/Warehouse_class/classes/DataStorage.m @@ -0,0 +1,224 @@ +classdef DataStorage < handle + %DATASTORAGE Summary of this class goes here + % Detailed explanation goes here + + properties + inputParams = struct; + parameter = struct; + fn = []; + dim = []; + sto = {}; + % getPhysForIndex = struct; + % getIndexForPhys = struct; + end + + methods + + function obj = DataStorage(inputParams) + %DATASTORAGE Construct an instance of this class + % Detailed explanation goes here + + % scheiß Variablenname + obj.inputParams = inputParams; + + % Field Names + obj.fn = string(fieldnames(inputParams)); + + % _______________ + % Two dicts that map between physical and array index :-) + % Dictionary: 10 km -> 3 + % obj.getIndexForPhys = obj.buildIndexDict(); + + % Dictionary: 3 -> 10Km + % obj.getPhysForIndex = obj.buildPhysDict(); + + % _______________ + % This is the 2nd Idea -> ecery given Param will be a class + % instance of "Parameter", therin user can access the dicts and + % informations... Have not decided which way is best.. + obj = obj.buildParameter(); + + % get dimension of dataStorage + % e.g. if we have L = [1,2,10,80] and D=[8, 17, 21], we would + % need an array with dimesion [4,3]. + obj.dim = obj.getDimension(); + + % finally, create the main storage as cell array + obj.sto = struct; + + end + + function showInfo(obj) + disp("Data Structure with fields:"); + fprintf('%-12s', 'Name'); fprintf('%1s', '| '); fprintf('%0s ', 'Dimension'); fprintf('%4s', '| '); fprintf('%0s ', 'Physical Values'); fprintf('\n'); + disp('----------------------------------------------------------------'); + for i = 1:numel(obj.fn) + + fprintf('%-12s', char(obj.fn(i))); fprintf('%1s', '| '); fprintf('%8.1i ', obj.dim(i)); fprintf('%5s', '| '); fprintf('%-7s ', string(obj.parameter.(obj.fn(i)).values) ); fprintf('\n'); + + end + disp('----------------------------------------------------------------'); + + stofn = string(fieldnames(obj.sto)); + for s = 1:numel(stofn) + nonempty = numel(find(~cellfun(@isempty,obj.sto.(stofn(s))))); + overall = numel(obj.sto.(stofn(s))); + fprintf('%-8s', 'Storage '); fprintf('%-10s', char(stofn(s))); fprintf('%4s', 'filled with '); fprintf('%-5s', num2str(nonempty)); fprintf('%-1s', ' entries -> '); fprintf('%-5s', num2str(nonempty/overall*100)); fprintf('%-1s', '% filled'); fprintf('\n'); + + end + end + + function dim = getDimension(obj) + dim = zeros(1,numel(obj.fn)); + for p = 1:numel(obj.fn) + %loop all parameter names and add their length :-) + dim(p)=obj.parameter.(obj.fn(p)).length; + end + end + + function obj = buildParameter(obj) + + for p = 1:numel(obj.fn) + name = obj.fn(p); + values = obj.inputParams.(name); + obj.parameter.(name) = Parameter(name,values); + end + + end + + function addStorage(obj,varName) + % add a storage + + storage = cell(obj.dim); + + obj.sto.(string(varName)) = storage; + + end + + function addValueToStorage(obj, valueToStore ,storageVarName, varargin) + + if nargin-3 == numel(obj.fn) + lin_idx = obj.getIndicesByPhys(varargin); + obj.sto.(storageVarName){lin_idx} = valueToStore; + else + error('Specify all the indices to access the right place in storage!') + end + + end + + % Access Value(s) + function value = getStoValue(obj,storageVarName, varargin) + + if nargin-2 == numel(obj.fn) + %es wurden aausreichend argumente übergeben :-) + %es gibt jedoch erstmal keinen Check ob die Reihenfolge + %richtig ist + value = []; + lin_idx = obj.getIndicesByPhys(varargin); + errcnt = 0; + for i=1:numel(lin_idx) + try + tmp = obj.sto.(storageVarName){lin_idx(i)}; + if ~isempty(tmp) + value(i,:) = tmp ; + else + errcnt = errcnt+1; + + if errcnt < 3 + %get back the n-dimensional subiondices... + [sub{1:length(size(obj.sto.(storageVarName)))}] = ind2sub(size(obj.sto.(storageVarName)),lin_idx(i)); + + %get back the physical representaion + word = []; + for phys_idx = 1:numel(obj.fn) + parametername = obj.fn(phys_idx); + word = [word,char(parametername),': ', num2str(obj.parameter.(parametername).getPhysForIndex(sub{phys_idx})),' ;']; + end + % warning(['Requested Data is not in Warehouse ', word]); + elseif errcnt == 3 + % warning(['... ', word]); + end + + end + if errcnt > 2 + % warning([num2str(errcnt),' requested datapoint(s) not in warehouse.']); + end + + catch + error('Error in Datastorage: Something happened while looking up in warehouse.') + end + end + else + error('Wrong Request using ExampleWarehouse.getStoValue(*parameter set*). Give me all the Parameters! Please!') + end + + end + + % Mapping for several Indices (calls the mapping for single index) + function lin_idx = getIndicesByPhys(obj,varargin) + %map _all_ phys. to indices - several calls of single mapping + + inputsz = cellfun(@size,varargin{1},'UniformOutput',false); + inputmax = cell2mat(cellfun(@max,inputsz,'UniformOutput',false)); + + % I wrote this method for a single query, then refined it for vectorial + % inputs (which work fine), but lastly recognized that I + % destroyed single query... however, this if / else will fix it! + if sum(inputmax)==length(inputmax) + vecQuery = 1; %set any value to one + else + vecQuery = find(inputmax~=1); %position of vectorial queries + end + + q={}; + for i = 1:numel(vecQuery) + q{i} = varargin{1,1}{vecQuery(i)}; %the two vectors + end + + combination = combvec(q{:}); %combine all possible combinations + + for c = 1:length(combination) + %loop over all possible combinations + + indices = {}; + str = []; + + for r = 1:numel(vecQuery) + %replace varargin with current query (could have been renamed... however it works) + varargin{1}{vecQuery(r)} = combination(r,c); + end + + for p = 1:numel(obj.fn) + + curPhysQuery = varargin{1}{p}; %can be: a) single value // b) range + + curParameterName = obj.fn(p); + + indices{end+1} = obj.getIndexByPhys(curParameterName,curPhysQuery); + + % str = [str, ',indices{', num2str(p), '}']; + + end + + %append to index list :-) + fn_=fieldnames(obj.sto); + n_ = fn_{1}; + lin_idx(c,:) = sub2ind(size(obj.sto.(n_)),indices{:}); + % lin_idx(c,:) = eval(['sub2ind(size(obj.sto.',n_,')',str,');']); + + end + + + end + + % Mapping for single Index + function idx = getIndexByPhys(obj,fieldname,phys) + %map single phys to index + idx = obj.parameter.(fieldname).getIndexForPhys(phys); + end + + + + end +end + diff --git a/Classes/Warehouse_class/classes/Parameter.m b/Classes/Warehouse_class/classes/Parameter.m new file mode 100644 index 0000000..c5243f0 --- /dev/null +++ b/Classes/Warehouse_class/classes/Parameter.m @@ -0,0 +1,44 @@ +classdef Parameter < handle + %PARAMETER Summary of this class goes here + % Detailed explanation goes here + + properties + name + values + indices + length + getPhysForIndex = dictionary; + getIndexForPhys = dictionary; + end + + methods + function obj = Parameter(name, values) + %PARAMETER Construct + obj.name = name; + obj.values = values; + obj.indices = [1:numel(obj.values)]; + obj.length = length(obj.values); + obj = obj.buildIndexForPhysDict(); + obj = obj.buildPhysForIndexDict(); + end + + function obj = buildPhysForIndexDict(obj) + %take Index + %return Phys + for idx = 1:numel(obj.values) + obj.getPhysForIndex(idx) = obj.values(idx); + end + end + + + function obj = buildIndexForPhysDict(obj) + %take Phys + %return Index + for idx = 1:numel(obj.values) + obj.getIndexForPhys(obj.values(idx)) = idx; + end + end + + end +end + diff --git a/Classes/Warehouse_class/functions/fwm_plots/CompleteRoutine.m b/Classes/Warehouse_class/functions/fwm_plots/CompleteRoutine.m new file mode 100644 index 0000000..86cbe93 --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/CompleteRoutine.m @@ -0,0 +1,251 @@ +% Script, that shows the data management routine :-) + +loadExistingWareHouse = 0; + +if loadExistingWareHouse + + [file, path] = uigetfile(); + wh = load([path filesep file]); + wh = wh.wh; + wh.showInfo; + +else + + % 1) Define all your parameters, best practice directly constructs a + % structure + + params = struct; + + params.l = [2,10]; + + params.dispersion = [0]; + + params.sgm = [0]; + +% params.pol = ["YXYXYXYX","YXXYYXXY","YYYYYYYY"]; + params.pol = ["alternated","paired","copolarized"]; + + params.p_in = [3]; + + params.p_out = [-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2]; + + params.pmd = [0.1]; + + params.gamma = [0.0023]; + + params.realization = [1:20]; + + params.numchannels = [16]; + + params.center_wavelength = floor([getSweepWavelengths(35, 50e9, 1310)] .* 1000) ./ 1000 ; + params.center_wavelength = [1285 1287 1290 1292 1295]; + params.center_wavelength = 1310; + + params.channelspacing = [400e9]; + + params.random_zdw = [0]; + + %wh = warehouse :-) + wh = DataStorage(params); + + wh.showInfo; + + wh.addStorage("ber"); + + wh.addStorage("totalBer"); + +end + + + +%2) Simulate a bunch of data - TO BE IMPLEMENTED HERE - for now use scripts +%from Sebastian + +%3) Once the simulation folder is around, specifiy path and analyze dirs + +path = uigetdir('C:\Users\Silas\Documents\MATLAB\Raw_Cluster_Simulations'); + +allMat = getAllFilesInFolder(path,'.mat'); +allErr = getAllFilesInFolder(path,'.err'); + +%allMat = dir([path filesep '*.mat']); + +%allErr = dir([path filesep '*.err']); + +if numel(allMat) == 0 + warning('You defined an empty folder. Could not locate any .mat file.') +else + fprintf('%-20s', 'Err Files:'); fprintf('%-12s', num2str(numel(allErr))); fprintf('\n'); + fprintf('%-20s', 'Mat Files:'); fprintf('%-12s', num2str(numel(allMat))); fprintf('\n'); + fprintf('%-20s', 'Missing Mat Files:'); fprintf('%-12s', num2str(numel(allErr)-numel(allMat))); fprintf('\n'); +end + +%4) Now load that data + +f = waitbar(0,'Please wait...'); +cnt = 0; + +for num = 1:numel(allMat) + + fileName = allMat(num).name; + fileFolder = allMat(num).path; + fileExt = allMat(num).ext; +% + matFile = load([fileFolder filesep fileName fileExt]); + matFile = matFile.loop_data; + + % ____________________________________ + % FIND THE DATAPOINT CURRENTLY LOADED + zdw = 1310; + + channelplan = "symmetric"; + + channelspacing = str2double(strrep(regexp(fileName,'(_chsp)+([\d]*)','match'),'_chsp','')).*1e9; + + numchannels = str2double(strrep(regexp(fileName,'(ch)+(_)+([\d]*)','match'),'ch_','')); + + center_wavelength = str2double(insertAfter(strrep(regexp(fileName,'(lambda)+([\d]*)','match'),'lambda',''),4,'.')); + + center_wavelength = floor(center_wavelength * 1000) / 1000; + + if center_wavelength == 2192 + continue + end + + center_wavelength = 1310; + + random_zdw = str2double(strrep(regexp(fileName,'(rzwd)+([\d])','match'),'rzwd','')); + + l = str2double(strrep(regexp(fileName,'([L])+(_)+([\d]*)','match'),'L_','')); + + d = str2double(strrep(regexp(fileName,'([D])+(_)+([\d]*)','match'),'D_','')); + + if d == 0 + sgm = false; + else + sgm = true; + end + + + if numel(regexp(fileName,'(YYYY)','match')) > 1 + pol = "copolarized"; + elseif numel(regexp(fileName,'(YXXY)','match')) > 1 + pol = "paired"; + elseif numel(regexp(fileName,'(YXYX)','match')) > 1 + pol = "alternated"; + else + pol = "copolarized"; + end + + p_in = str2double(strrep(regexp(fileName,'(pow_)+([-,\d]{1})','match'),'pow_','')); + + pmd = 0.1; + + gamma = 0.0023; + + realiz = str2double(strrep(regexp(fileName,'(r)+([-,\d]{1,3})','match'),'r','')); + + + + % ____________________________________ + % Get the information you want from current file + rop=[]; + ber = []; + for pow = 2:12 + + module_number = ''; + for p = 1:11 %11 because there are 11 ROP branches in model + + % get ROP + if p == 1 + p_out = matFile.dp_optatten_para.atten; + else + p_out = matFile.("dp_optatten__"+(p)+"_para").atten; + end + + p_out = round(p_out-10*log10(numel(matFile.config.parameters.common.wavelengthPlan))); + + for c = 1:numel(matFile.config.parameters.common.wavelengthPlan) + + ber(c) = matFile.("prms_compare_wdm__"+(p+1)+"_out"){1,c}.ber; + + end + + totalBer = matFile.("prms_compare_wdm__"+(p+1)+"_out"){1,end}.totalBer; + + if totalBer > 0.2 && channelspacing == 400e9 && pol == "alternated" + disp("stopping here"); + pause; + end + + + % ____________________________________ + % Add value to warehouse at the correct position + + + + wh.addValueToStorage(ber,'ber',l,d,sgm,pol,p_in,p_out,pmd,gamma,realiz,numchannels, center_wavelength,channelspacing,random_zdw); + wh.getStoValue('ber',l,d,sgm,pol,p_in,p_out,pmd,gamma,realiz,numchannels, center_wavelength,channelspacing,random_zdw); + wh.addValueToStorage(totalBer,'totalBer',l,d,sgm,pol,p_in,p_out,pmd,gamma,realiz,numchannels,center_wavelength,channelspacing,random_zdw); + + end + + end + waitbar(num/numel(allMat),f,'Loading your data'); +end + +close(f) + + + + + +% 4) Hey! the warehouse is here and (hopefully) filled with data :-) + +% Create a save dialog +defaultDir = 'C:\Users\Silas\Documents\MATLAB\Raw_Cluster_Simulations\'; +defaultExt = '*.mat'; +[filename, pathname] = uiputfile(fullfile(defaultDir, defaultExt),'', 'wh.mat'); + +% Check if the user pressed Cancel +if isequal(filename, 0) || isequal(pathname, 0) + disp('Save operation canceled.'); +else + % Save the variable to the selected file + save(fullfile(pathname, filename), 'wh'); + disp(['Variable "wh" saved to: ', fullfile(pathname, filename)]); +end + + + + +function matFileStructArray = getAllFilesInFolder(folderPath,extension) + % Get a list of all files in the current folder + currentFolderFiles = dir(fullfile(folderPath, '*')); + + % Exclude '.' and '..' directories + currentFolderFiles = currentFolderFiles(~ismember({currentFolderFiles.name}, {'.', '..'})); + + % Initialize the structure array for .mat files + matFileStructArray = struct('path', {}, 'name', {}, 'ext', {}); + + % Loop over each file in the current folder + for i = 1:length(currentFolderFiles) + currentFile = currentFolderFiles(i); + + % Check if the current item is a file and has a .mat extension + if ~currentFile.isdir && endsWith(currentFile.name, extension, 'IgnoreCase', true) + % If it's a .mat file, add it to the structure array + [matFileStructArray(end + 1).path,matFileStructArray(end+1).name, matFileStructArray(end+1).ext] = fileparts(fullfile(folderPath, currentFile.name)); + elseif currentFile.isdir + % If it's a directory, recursively call the function + subfolderPath = fullfile(folderPath, currentFile.name); + subfolderMatFiles = getAllFilesInFolder(subfolderPath,extension); + + % Add .mat files from the subfolder to the structure array + matFileStructArray = [matFileStructArray, subfolderMatFiles]; + end + end +end + + diff --git a/Classes/Warehouse_class/functions/fwm_plots/CompleteRoutine_DifferentChannels_Fig3.m b/Classes/Warehouse_class/functions/fwm_plots/CompleteRoutine_DifferentChannels_Fig3.m new file mode 100644 index 0000000..aaec117 --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/CompleteRoutine_DifferentChannels_Fig3.m @@ -0,0 +1,250 @@ +% Script, that shows the data management routine :-) + +loadExistingWareHouse = 0; + +if loadExistingWareHouse + + [file, path] = uigetfile(); + wh = load([path filesep file]); + wh = wh.wh; + wh.showInfo; + +else + + % 1) Define all your parameters, best practice directly constructs a + % structure + + params = struct; + + params.l = [2, 10]; + + params.dispersion = [0, 3]; + + params.sgm = [0, 1]; + +% params.pol = ["YXYXYXYX","YXXYYXXY","YYYYYYYY"]; + params.pol = ["alternated","paired","copolarized"]; + + params.p_in = [3]; + + params.p_out = [-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2]; + + params.pmd = [0.1]; + + params.gamma = [0.0023]; + + params.realization = [1:20]; + + params.numchannels = [1,2,4,8,16]; + + params.center_wavelength = floor([getSweepWavelengths(35, 50e9, 1310)] .* 1000) ./ 1000 ; + params.center_wavelength = [1285 1287 1290 1292 1295]; + params.center_wavelength = 1310; + + params.channelspacing = [400e9]; + + params.random_zdw = [0,1]; + + %wh = warehouse :-) + wh = DataStorage(params); + + wh.showInfo; + + wh.addStorage("ber"); + + wh.addStorage("totalBer"); + +end + + + +%2) Simulate a bunch of data - TO BE IMPLEMENTED HERE - for now use scripts +%from Sebastian + +%3) Once the simulation folder is around, specifiy path and analyze dirs + +path = uigetdir('C:\Users\Silas\Documents\MATLAB\Raw_Cluster_Simulations'); + +allMat = getAllFilesInFolder(path,'.mat'); +allErr = getAllFilesInFolder(path,'.err'); + +%allMat = dir([path filesep '*.mat']); + +%allErr = dir([path filesep '*.err']); + +if numel(allMat) == 0 + warning('You defined an empty folder. Could not locate any .mat file.') +else + fprintf('%-20s', 'Err Files:'); fprintf('%-12s', num2str(numel(allErr))); fprintf('\n'); + fprintf('%-20s', 'Mat Files:'); fprintf('%-12s', num2str(numel(allMat))); fprintf('\n'); + fprintf('%-20s', 'Missing Mat Files:'); fprintf('%-12s', num2str(numel(allErr)-numel(allMat))); fprintf('\n'); +end + +%4) Now load that data + +f = waitbar(0,'Please wait...'); +cnt = 0; + +for num = 1:numel(allMat) + + fileName = allMat(num).name; + fileFolder = allMat(num).path; + fileExt = allMat(num).ext; +% + matFile = load([fileFolder filesep fileName fileExt]); + matFile = matFile.loop_data; + + % ____________________________________ + % FIND THE DATAPOINT CURRENTLY LOADED + zdw = 1310; + + channelplan = "symmetric"; + + channelspacing = str2double(strrep(regexp(fileName,'(_chsp)+([\d]*)','match'),'_chsp','')).*1e9; + + numchannels = str2double(strrep(regexp(fileName,'(ch)+(_)+([\d]*)','match'),'ch_','')); + + center_wavelength = str2double(insertAfter(strrep(regexp(fileName,'(lambda)+([\d]*)','match'),'lambda',''),4,'.')); + + center_wavelength = floor(center_wavelength * 1000) / 1000; + + if center_wavelength == 2192 + continue + end + + center_wavelength = 1310; + + random_zdw = str2double(strrep(regexp(fileName,'(rzwd)+([\d])','match'),'rzwd','')); + + l = str2double(strrep(regexp(fileName,'([L])+(_)+([\d]*)','match'),'L_','')); + + d = str2double(strrep(regexp(fileName,'([D])+(_)+([\d]*)','match'),'D_','')); + + if d == 0 + sgm = false; + else + sgm = true; + end + + + if numel(regexp(fileName,'(YYYY)','match')) > 1 + pol = "copolarized"; + elseif numel(regexp(fileName,'(YXXY)','match')) > 1 + pol = "paired"; + elseif numel(regexp(fileName,'(YXYX)','match')) > 1 + pol = "alternated"; + else + pol = "copolarized"; + end + + p_in = str2double(strrep(regexp(fileName,'(pow_)+([-,\d]{1})','match'),'pow_','')); + + pmd = 0.1; + + gamma = 0.0023; + + realiz = str2double(strrep(regexp(fileName,'(r)+([-,\d]{1,3})','match'),'r','')); + + + + % ____________________________________ + % Get the information you want from current file + rop=[]; + ber = []; + for pow = 2:12 + + module_number = ''; + for p = 1:11 %11 because there are 11 ROP branches in model + + % get ROP + if p == 1 + p_out = matFile.dp_optatten_para.atten; + else + p_out = matFile.("dp_optatten__"+(p)+"_para").atten; + end + + p_out = round(p_out-10*log10(numel(matFile.config.parameters.common.wavelengthPlan))); + + for c = 1:numel(matFile.config.parameters.common.wavelengthPlan) + + ber(c) = matFile.("prms_compare_wdm__"+(p+1)+"_out"){1,c}.ber; + + end + + totalBer = matFile.("prms_compare_wdm__"+(p+1)+"_out"){1,end}.totalBer; + + if totalBer > 0.2 && channelspacing == 400e9 && pol == "alternated" + disp("stopping here"); + pause; + end + + % ____________________________________ + % Add value to warehouse at the correct position + + + + wh.addValueToStorage(ber,'ber',l,d,sgm,pol,p_in,p_out,pmd,gamma,realiz,numchannels, center_wavelength,channelspacing,random_zdw); + wh.getStoValue('ber',l,d,sgm,pol,p_in,p_out,pmd,gamma,realiz,numchannels, center_wavelength,channelspacing,random_zdw); + wh.addValueToStorage(totalBer,'totalBer',l,d,sgm,pol,p_in,p_out,pmd,gamma,realiz,numchannels,center_wavelength,channelspacing,random_zdw); + + end + + end + waitbar(num/numel(allMat),f,'Loading your data'); +end + +close(f) + + + + + +% 4) Hey! the warehouse is here and (hopefully) filled with data :-) + +% Create a save dialog +defaultDir = 'C:\Users\Silas\Documents\MATLAB\Raw_Cluster_Simulations\'; +defaultExt = '*.mat'; +[filename, pathname] = uiputfile(fullfile(defaultDir, defaultExt),'', 'wh.mat'); + +% Check if the user pressed Cancel +if isequal(filename, 0) || isequal(pathname, 0) + disp('Save operation canceled.'); +else + % Save the variable to the selected file + save(fullfile(pathname, filename), 'wh'); + disp(['Variable "wh" saved to: ', fullfile(pathname, filename)]); +end + + + + +function matFileStructArray = getAllFilesInFolder(folderPath,extension) + % Get a list of all files in the current folder + currentFolderFiles = dir(fullfile(folderPath, '*')); + + % Exclude '.' and '..' directories + currentFolderFiles = currentFolderFiles(~ismember({currentFolderFiles.name}, {'.', '..'})); + + % Initialize the structure array for .mat files + matFileStructArray = struct('path', {}, 'name', {}, 'ext', {}); + + % Loop over each file in the current folder + for i = 1:length(currentFolderFiles) + currentFile = currentFolderFiles(i); + + % Check if the current item is a file and has a .mat extension + if ~currentFile.isdir && endsWith(currentFile.name, extension, 'IgnoreCase', true) + % If it's a .mat file, add it to the structure array + [matFileStructArray(end + 1).path,matFileStructArray(end+1).name, matFileStructArray(end+1).ext] = fileparts(fullfile(folderPath, currentFile.name)); + elseif currentFile.isdir + % If it's a directory, recursively call the function + subfolderPath = fullfile(folderPath, currentFile.name); + subfolderMatFiles = getAllFilesInFolder(subfolderPath,extension); + + % Add .mat files from the subfolder to the structure array + matFileStructArray = [matFileStructArray, subfolderMatFiles]; + end + end +end + + diff --git a/Classes/Warehouse_class/functions/fwm_plots/automate_JLT_plots.m b/Classes/Warehouse_class/functions/fwm_plots/automate_JLT_plots.m new file mode 100644 index 0000000..68b77b6 --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/automate_JLT_plots.m @@ -0,0 +1,394 @@ + + +%automate plots +[file, path] = uigetfile("C:\Users\Silas\Documents\MATLAB\Datensätze\Raw_Cluster_Simulations\session_februar_24\wh_mi_nacht.mat"); +wh = load([path filesep file]); +wh = wh.wh; + +plotJob = struct(); +width = 350; +height = 200; +plotJob.Position = [100 100 width 100+height]; +cols = cbrewer2("paired",12); +plotJob.color = cols(1,:); +plotJob.l = 2; +plotJob.ch = 16; +plotJob.d = 0; +plotJob.sgm = 0; +plotJob.pol = "copolarized"; +plotJob.p_in = 3; +plotJob.gamma = 0; +plotJob.pmd = 0; +plotJob.channelspacing = 400e9; +plotJob.randzdw = 0; + + +plotJob.plot_ber_curve = 1; +plotJob.plot_3dber_curve = 0; +plotJob.plot_violin = 0; +plotJob.plot_wavelength_sweep = 0; +plotJob.plot_wavelength_sweep_failure_rate = 0; + + +plotJob.dataStatArg = 'Lineplot with quartiles'; +plotJob.plotTypeArg = 'Lines'; +plotJob.displayname = 'a'; +plotJob.title = 'title'; +plotJob.figName = '16 Chann__'; +plotJob.xAxisLabel = 'ROP per Channel in dBm'; +plotJob.yAxisLabel = 'BER'; + +% createbercurves(wh,plotJob) +createviolinplots(wh,plotJob); +% createsweepplots(wh,plotJob); + + +%% 1 +function createbercurves(wh,plotJob) +width = 1650; +height = 400; +s = 100; +e = 100; + +cols = cbrewer2("paired",12); +numRows = 2; +numCols = 4; + +plotJob.figName = '16 Chann_200G'; +plotJob.channelspacing = 400e9; +plotJob.ch = 16; +Len = [2,2,2,2,10,10,10,10]; +Pol = ["copolarized","alternated","paired","copolarized","copolarized","alternated","paired","copolarized"]; +Title = ["Co Polarized","Alternating Pol. Interl.","Paired Pol. Interl.","Link Segmentation",]; +D = [0,0,0,3,0,0,0,3]; +Sgm = [0,0,0,1,0,0,0,1]; + +colidx = [4,8,6]; +P_launch = [0,3,6]; + +fig = figure('Name',plotJob.figName); +fig.Position = plotJob.Position; +fig.Units = "centimeters"; +fig.Position = [0 0 18 7]; +t = tiledlayout(numRows,numCols,'TileSpacing','compact','Padding','compact'); +for idx = 1:(numRows * numCols) + % Create subplot + % sp = subplot(numRows, numCols, idx); + nexttile; + plotJob.l = Len(idx); + plotJob.pol = Pol(idx); + plotJob.d = D(idx); + plotJob.sgm = Sgm(idx); + plotJob.randzdw = 1; + + for i = 1:3 + + plotJob.p_in = P_launch(i); + plotJob.color = cols(colidx(i),:); + hold on + plotCurve(wh, plotJob); + + end + % + if idx ~= 1 && idx ~= 5 % For example, hide y-axis for subplot 1 + set(gca, 'YTickLabel',[]); % Hide y-axis ticks and labels + set(gca,'YGrid','on'); + set(gca, 'YLabel', []); + end + if idx ~= 5 && idx ~= 6 && idx ~= 7 && idx ~= 8 + set(gca, 'XLabel', []); + set(gca, 'XTickLabel', []); + + end + grid on + + g = gca; + pos = g.Position; + if idx <= 4 + title(Title(idx),'FontSize',8); + % a = annotation('textbox', pos-[0.0020 -0.1434 0.0947 0.3121], 'String', "FEC: 3.8e-3","LineStyle","none","FontSize",8,'LineWidth',0.1,'Interpreter','latex',"FontUnits","points",'Vert','middle','FitBoxToText','on'); + % a = annotation('textbox', pos, 'String', "2 km","LineStyle","none","FontSize",8,'LineWidth',0.1,'Interpreter','latex',"FontUnits","points",'Vert','bottom','FitBoxToText','on'); + else + % a = annotation('textbox', pos, 'String', "FEC: 3.8e-3","LineStyle","none","FontSize",8,'LineWidth',0.1,'Interpreter','latex',"FontUnits","points",'Vert','middle','FitBoxToText','on'); + % a = annotation('textbox', pos, 'String', "10 km","LineStyle","none","FontSize",8,'LineWidth',0.1,'Interpreter','latex',"FontUnits","points",'Vert','bottom','FitBoxToText','on'); + end + +end + +% Create textbox +annotation(fig,'textbox',... + [0.0696078431372547 0.246851385390432 0.0656862745098043 0.0453400503778337],... + 'String','10 km',... + 'LineStyle','none',... + 'Interpreter','latex',... + 'FontSize',8,... + 'FitBoxToText','off'); + +% Create textbox +annotation(fig,'textbox',... + [0.288235294117647 0.239294710327459 0.0656862745098043 0.0453400503778338],... + 'String','10 km',... + 'LineStyle','none',... + 'Interpreter','latex',... + 'FontSize',8,... + 'FitBoxToText','off'); + +% Create textbox +annotation(fig,'textbox',... + [0.516666666666666 0.241813602015117 0.0656862745098042 0.0453400503778338],... + 'String','10 km',... + 'LineStyle','none',... + 'Interpreter','latex',... + 'FontSize',8,... + 'FitBoxToText','off'); + +% Create textbox +annotation(fig,'textbox',... + [0.742156862745097 0.236775818639802 0.0656862745098042 0.0453400503778339],... + 'String','10 km',... + 'LineStyle','none',... + 'Interpreter','latex',... + 'FontSize',8,... + 'FitBoxToText','off'); + +% Create textbox +annotation(fig,'textbox',... + [0.071996471804854 0.578899159967702 0.0656862745098039 0.0453400503778341],... + 'String',{'2 km'},... + 'LineStyle','none',... + 'Interpreter','latex',... + 'FontSize',8,... + 'FitBoxToText','off'); + +% Create textbox +annotation(fig,'textbox',... + [0.29596893566113 0.576253721089996 0.0656862745098041 0.0453400503778341],... + 'String',{'2 km'},... + 'LineStyle','none',... + 'Interpreter','latex',... + 'FontSize',8,... + 'FitBoxToText','off'); + +% Create textbox +annotation(fig,'textbox',... + [0.524883914268912 0.580785315705112 0.0656862745098041 0.0453400503778341],... + 'String',{'2 km'},... + 'LineStyle','none',... + 'Interpreter','latex',... + 'FontSize',8,... + 'FitBoxToText','off'); + +% Create textbox +annotation(fig,'textbox',... + [0.753758338909501 0.581291504465311 0.0656862745098037 0.0453400503778341],... + 'String',{'2 km'},... + 'LineStyle','none',... + 'Interpreter','latex',... + 'FontSize',8,... + 'FitBoxToText','off'); + +a=sgtitle(['N=',num2str(plotJob.ch),'; $\Delta f_{\mathrm{ch}}$= ',num2str(plotJob.channelspacing*1e-9),' GHz'],'FontSIze',10); +a.Interpreter = "latex"; + +lgd = legend('$P_{\mathrm{in}}=0$ dBm','$P_{\mathrm{in}}=3$ dBm','$P_{\mathrm{in}}=6$ dBm','Interpreter','latex'); +lgd.NumColumns = 3; +lgd.Layout.Tile = 'south'; + +copygraphics(t,'BackgroundColor','none'); +end + +%% 2 +function createviolinplots(wh,plotJob) +width = 350; +height = 200; +s = 100; +e = 100; + +cols = cbrewer2("paired",12); +numRows = 1; +numCols = 4; + + +plotJob.ch = 16; +plotJob.p_in = 3; +plotJob.randzdw = 0; + +Pol = ["copolarized","copolarized","alternated","paired",]; +Title = ["Co Pol.","Link Segmentation","Paired Pol. Interl.","Alternating Pol. Interl."]; +D = [0,3,0,0]; +Sgm = [0,1,0,0]; + +colidx = [2]; +Len = [2]; + +plotJob.figName = ['_Violin',num2str(plotJob.ch),' Channels; ',num2str(plotJob.channelspacing*1e-9),' GHz; ',num2str(plotJob.p_in),' dBm; randomized: ', num2str(plotJob.randzdw)]; +fig = figure('Name',plotJob.figName); +fig.Position = plotJob.Position; +fig.Units = "centimeters"; +fig.Position = [0 0 18 7]; + +t = tiledlayout(numRows,numCols,'TileSpacing','compact','Padding','compact'); +for idx = 1:(numRows * numCols) + % Create subplot + %subplot(numRows, numCols, idx); + nexttile; + + plotJob.pol = Pol(idx); + plotJob.d = D(idx); + plotJob.sgm = Sgm(idx); + + for i = 1 + + plotJob.color = cols(colidx(i),:); + plotJob.l = Len(i); + hold on + plotViolin(wh, plotJob); + + end + + if idx ~= 1 % For example, hide y-axis for subplot 1 + %set(gca, 'YTickLabel',[]); % Hide y-axis ticks and labels + set(gca, 'YGrid','on'); + set(gca, 'YLabel', []); + end + + if idx <= 4 + title(Title(idx)); + end + + +end + +% Create textbox +annotation(fig,'textbox',... + [0.300019607843137 0.816120906801009 0.108803921568628 0.0906801007556676],... + 'String','$P_{\mathrm{in}}=3$ dBm',... + 'LineStyle','none',... + 'Interpreter','latex',... + 'FontSize',8,... + 'FitBoxToText','off'); + +% Create textbox +annotation(fig,'textbox',... + [0.530411764705882 0.81612090680101 0.108803921568628 0.0906801007556676],... + 'String','$P_{\mathrm{in}}=3$ dBm',... + 'LineStyle','none',... + 'Interpreter','latex',... + 'FontSize',8,... + 'FitBoxToText','off'); + +% Create textbox +annotation(fig,'textbox',... + [0.755901960784313 0.816120906801011 0.108803921568628 0.0906801007556676],... + 'String','$P_{\mathrm{in}}=3$ dBm',... + 'LineStyle','none',... + 'Interpreter','latex',... + 'FontSize',8,... + 'FitBoxToText','off'); + +% Create textbox +annotation(fig,'textbox',... + [0.0696274509803918 0.584382871536529 0.108803921568628 0.0906801007556676],... + 'String','$P_{\mathrm{in}}=3$ dBm',... + 'LineStyle','none',... + 'Interpreter','latex',... + 'FontSize',8,... + 'FitBoxToText','off'); + +copygraphics(t,'BackgroundColor','none'); +% lgd = legend('$P_{\mathrm{in}}=0$ dBm','$P_{\mathrm{in}}=3$ dBm','$P_{\mathrm{in}}=6$ dBm','Interpreter','latex'); +% lgd.NumColumns = 3; +% lgd.Layout.Tile = 'south'; + +end + +%% 3 +function createsweepplots(wh,plotJob) + +width = 650; +height = 200; +s = 100; +e = 100; +plotJob.Position = [0 0 width e+height]; + +cols = cbrewer2("paired",12); +numRows = 1; +numCols = 4; + + +plotJob.channelspacing = 200e9; +plotJob.ch = 16; +plotJob.randzdw = 1; +plotJob.l = 10; + +plotJob.p_in = 3; + +Pol = ["copolarized","alternated","paired","copolarized"]; +Title = ["Co Polarized","Alternating Pol. Interl.","Paired Pol. Interl.","Link Segmentation",]; +D = [0,0,0,3]; +Sgm = [0,0,0,1]; +Channelspacing = [200e9, 200e9]; +PlotTypeArg = ["--","-"]; +colidx = [6,8,2,4]; +Len = [2,10]; + +plotJob.figName = [num2str(plotJob.ch),num2str(plotJob.channelspacing*1e-9),num2str(plotJob.p_in),'...']; +plotJob.figName = "10km 400ghz"; + + + + +fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName); + +if isvalid(fig) + figure(fig) + fig = get(fig); + AxesMain = fig.CurrentAxes; + hold on +else + fig = figure('name',char(plotJob.figName)); + AxesMain = gca; + hold on +end + +fig.Position = plotJob.Position; +fig.Units = "centimeters"; +fig.Position = [0 0 18 7]; + +for j = 1 + + plotJob.channelspacing = Channelspacing(j); + plotJob.plotTypeArg = PlotTypeArg(j); + + for idx = 1:4 + + plotJob.color = cols(colidx(idx),:); + plotJob.pol = Pol(idx); + plotJob.d = D(idx); + plotJob.sgm = Sgm(idx); + plotJob.displayname = [char(plotJob.pol)]; + hold on + plotBerVsZdwFailureRate(wh, plotJob); + + end +end +legend('Location', 'southoutside', 'Orientation', 'horizontal'); + + +%plot channel positions +hold on +chpos = calcWavelengthPlan(plotJob.ch, plotJob.channelspacing, 1310); +xline(chpos,'LineWidth',2,'Alpha',0.4,'HandleVisibility','off'); + +chpos = calcWavelengthPlan(plotJob.ch, plotJob.channelspacing, chpos(4)); +xline(chpos,'LineWidth',2,'LineStyle','--','Alpha',0.1,'HandleVisibility','off'); + +title(['N=',num2str(plotJob.ch),' $\Delta f_{\mathrm{ch}}$= ',num2str(plotJob.channelspacing*1e-9),' GHz'],'FontSize',10,'Interpreter','latex'); + + + + +end + + + diff --git a/Classes/Warehouse_class/functions/fwm_plots/dispersion_only.m b/Classes/Warehouse_class/functions/fwm_plots/dispersion_only.m new file mode 100644 index 0000000..a0e8a98 --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/dispersion_only.m @@ -0,0 +1,83 @@ +%automate plots +[file, path] = uigetfile("C:\Users\Silas\Documents\MATLAB\Datensätze\Raw_Cluster_Simulations\"); +wh = load([path filesep file]); +wh = wh.wh; + +plotJob = struct(); +width = 350; +height = 200; +plotJob.Position = [100 100 width 100+height]; +cols = cbrewer2("paired",12); +plotJob.color = cols(1,:); +plotJob.l = 1; +plotJob.ch = 1; + +plotJob.sgm = 1; +plotJob.pol = "copolarized"; +plotJob.p_in = 3; +plotJob.gamma = 0.0023; +plotJob.pmd = 0.1; +plotJob.channelspacing = 400e9; +plotJob.randzdw = 0; + + +plotJob.plot_ber_curve = 1; +plotJob.plot_3dber_curve = 0; +plotJob.plot_violin = 0; +plotJob.plot_wavelength_sweep = 0; +plotJob.plot_wavelength_sweep_failure_rate = 0; + + +plotJob.dataStatArg = 'Lineplot with quartiles'; +plotJob.plotTypeArg = 'Lines'; +plotJob.displayname = 'a'; +plotJob.title = 'title'; +plotJob.figName = '1 Chann__'; +plotJob.xAxisLabel = 'ROP per Channel in dBm'; +plotJob.yAxisLabel = 'BER'; + + +plotJob.d = 0; + +xAxis = wh.parameter.p_out.values; +D = wh.parameter.dispersion.values; + +figure() +ber_ = []; +for d_ = 0:39 + if d_ == 0 + plotJob.sgm = 0; + ber_(d_+1,:) = wh.getStoValue('ber',plotJob.l,d_,plotJob.sgm,string(plotJob.pol),plotJob.p_in,xAxis,plotJob.pmd,plotJob.gamma,1,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw)'; + else + plotJob.sgm = 1; + ber_(d_+1,:) = wh.getStoValue('ber',plotJob.l,d_,plotJob.sgm,string(plotJob.pol),plotJob.p_in,xAxis,plotJob.pmd,plotJob.gamma,1,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw)'; + end + hold on + plot(xAxis,ber_(d_+1,:)) + set(gca,'yscale','log'); +end +yline(3.8e-3); + + +hdfec = 3.8e-3.*ones(size(xAxis)); +for i = 1:size(ber_,1) + ber_series = ber_(i,:); + a = InterX([hdfec(:)';xAxis],[ber_series;xAxis]); + cross(i) = a(2); +end + +col = cbrewer2('Paired',8); +figure() +plot(D,cross,'Marker','o','MarkerSize',5,'MarkerEdgeColor',[1,1,1],'MarkerFaceColor',col(2,:),'Color',col(1,:),'LineWidth',1); +grid minor +xlabel('Accumulated Dispersion') +ylabel('Required ROP to reach FEC limit in dB') +line([D(16),D(16)],[-10,cross(16)],'linestyle','--') +line([0,D(16)],[cross(16),cross(16)],'linestyle','--') + +line([D(29),D(29)],[-10,cross(29)],'linestyle','--') +line([0,D(29)],[cross(29),cross(29)],'linestyle','--') + + + + diff --git a/Classes/Warehouse_class/functions/fwm_plots/dispersion_validation_miniskript.m b/Classes/Warehouse_class/functions/fwm_plots/dispersion_validation_miniskript.m new file mode 100644 index 0000000..802db0a --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/dispersion_validation_miniskript.m @@ -0,0 +1,52 @@ + +wh = load("C:\Users\Silas\Documents\MATLAB\Datensätze\Raw_Cluster_Simulations\session_dispersion_validation\wh_with_variation.mat"); +wh = wh.wh; + +lambda = 1295; + +figure(3) +plot(xAxis,getber(wh,lambda),'DisplayName',['w:', num2str(lambda), 'variation']); +yline(3.8e-3,'HandleVisibility','off'); +legend +set(gca,'yscale','log'); +grid(gca,'on'); +grid(gca,'minor'); +grid minor +fontsize(gca,8,"points") +fig.Units = "centimeters"; +fig.Position = [2 2 8.5 7]; +set(gca,'TickLabelInterpreter','latex') +ylim([1e-5,0.5]); +xlim([min(xAxis),-3]); + +wh = load("C:\Users\Silas\Documents\MATLAB\Datensätze\Raw_Cluster_Simulations\session_dispersion_validation\wh_no_variation.mat"); +wh = wh.wh; + + +figure(3) +hold on +plot(xAxis,getber(wh,lambda),'DisplayName',['w:', num2str(lambda), 'no variation']); +yline(3.8e-3,'HandleVisibility','off'); +legend +set(gca,'yscale','log'); +grid(gca,'on'); +grid(gca,'minor'); +grid minor +fontsize(gca,8,"points") +fig.Units = "centimeters"; +fig.Position = [2 2 8.5 7]; +set(gca,'TickLabelInterpreter','latex') +ylim([1e-5,0.5]); +xlim([min(xAxis),-3]); + + +function ber = getber(wh,lambda) + realization = wh.parameter.realization.values(1:end); + xAxis = wh.parameter.p_out.values; + ber = []; + for xl = 1:numel(xAxis) + p_out = xAxis(xl); + temp = wh.getStoValue('ber',10,0,0,"copolarized",3,p_out,0.1,0.0023,realization,1,lambda,400e9,1); + ber(xl) = mean(temp,'all'); + end +end \ No newline at end of file diff --git a/Classes/Warehouse_class/functions/fwm_plots/generatePlots.m b/Classes/Warehouse_class/functions/fwm_plots/generatePlots.m new file mode 100644 index 0000000..2e720ff --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/generatePlots.m @@ -0,0 +1,61 @@ +function generatePlots(wh,plotJob) + +% 0) Test for valid query: +p_out = wh.parameter.p_out.values(1); +realization = 9; + + +if 1 %~isempty(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310)) + % test violin + + baseName = plotJob.figName; + + width = 350; + height = 200; + s = 100; + e = 100; + + if plotJob.plot_ber_curve + plotJob.Position = [100 100 width e+height]; + plotJob.figName = [baseName, ' zdwvsber']; + plotCurve(wh, plotJob); + end + + if plotJob.plot_3dber_curve + plotJob.Position = [100 100 width e+height]; + plotJob.figName = [baseName, ' zdwvsber']; + plot3dCurve(wh, plotJob); + end + + if plotJob.plot_wavelength_sweep + plotJob.Position = [100 100 width e+height]; + plotJob.figName = [baseName, ' zdwvsber']; + plotBerVsZDW(wh, plotJob); + end + + if plotJob.plot_wavelength_sweep_failure_rate + plotJob.Position = [100 100 width e+height]; + plotJob.figName = [baseName, ' zdwvsber']; + plotBerVsZdwFailureRate(wh, plotJob); + end + + if plotJob.plot_violin + plotJob.Position = [s+width 100 width e+height]; + plotJob.figName = [baseName, ' violin']; + plotViolin(wh, plotJob); + end + + + if 0 + %2) plotHistogram + plotJob.Position = [s+2*width 100 width e+height]; + plotJob.figName = [baseName, ' FEC crossing']; + plotHistogram(wh,plotJob) + end + + +else + warndlg('The requested Datapoint is not available... This can occur for some edgecase constellations... ') +end + +end \ No newline at end of file diff --git a/Classes/Warehouse_class/functions/fwm_plots/plot3dCurve.m b/Classes/Warehouse_class/functions/fwm_plots/plot3dCurve.m new file mode 100644 index 0000000..37d8fcf --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/plot3dCurve.m @@ -0,0 +1,248 @@ +function plotCurve(wh,plotJob) +%PLOTCURVE Summary of this function goes here +% Detailed explanation goes here + +fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName); + +if isvalid(fig) + figure(fig) + fig = get(fig); + AxesMain = fig.CurrentAxes; + hold on +else + fig = figure('name',char(plotJob.figName)); + AxesMain = gca; + hold on +end + +col = plotJob.color; + +% we want to fetch all realizations +if plotJob.pmd == 0 + realization = 499; + % realization = 0:7; +else + realization = wh.parameter.realization.values(1:end); +end +realization = wh.parameter.realization.values(1:end); +% get all xAxis values +xAxis = wh.parameter.p_out.values; + +% Fetch Data from Warehouse +for xl = 1:numel(xAxis) + p_out = xAxis(xl); + if string(plotJob.dataStatArg) == "Worst" + temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw); + temp = removeZeros(temp); + ber(xl) = max(temp,[],'all'); + linew = 1.0; + markersz = 3; + linestyle = '-'; + elseif string(plotJob.dataStatArg) == "AVG" + temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw); + temp = removeZeros(temp); + ber(xl) = mean(temp,'all'); + linew = 1.0; + markersz = 3; + linestyle = '-'; + elseif string(plotJob.dataStatArg) == "Best" + temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw); + temp = removeZeros(temp); + ber(xl) = min(temp,[],'all'); + linew = 1.0; + markersz = 3; + linestyle = '-'; + + elseif string(plotJob.dataStatArg) == "All Channels; mean(PMD Realizations)" + temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw); + temp = removeZeros(temp); + + ber(:,xl) = mean(temp,1,"omitnan").'; + + linew = 1; + markersz = 3; + linestyle = '--'; + + elseif string(plotJob.dataStatArg) == "Lineplot with quartiles" + + temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw); + temp = removeZeros(temp); + ber(xl) = mean(temp,"all","omitnan").'; + + + upperq(xl) = quantile(temp,0.9,"all"); + lowerq(xl) = quantile(temp,0.1,"all"); + + if lowerq(xl) == 0 + lowerq(xl) = lowerq(xl-1); + end + % upperq(xl) = 0.5*std(tmp,1,'all','omitnan'); + % lowerq(xl) = 0.5*std(tmp,1,'all','omitnan'); + + % upperq(xl) = max(dataNoNans); + % lowerq(xl) = min(dataNoNans); + + % upperq(xl) = mean(tmp,"all","omitnan") + 1.96 * (std(tmp,1,'all','omitnan')/sqrt(numel(tmp))); + % lowerq(xl) = mean(tmp,"all","omitnan") - 1.96 * (std(tmp,1,'all','omitnan')/sqrt(numel(tmp))); + + linew = 1; + markersz = 1; + linestyle = '-'; + + elseif string(plotJob.dataStatArg) == "All Channels ;All PMD Realizations" + + tmp = reshape(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw).',[],1); + ber(1:size(tmp,1),xl) = tmp; + + linew = 0.3; + markersz = 2; + linestyle = ':'; + end +end + +xAxis = xAxis; + +% Plot Data +if string(plotJob.plotTypeArg) == "Scatter" + for rlz = 1:size(ber,1) + + if rlz < size(ber,1) + scatter(xAxis,ber(rlz,:),markersz,'MarkerEdgeColor',col,'MarkerFaceColor',col,'Parent', AxesMain(1),'HandleVisibility','off'); + else + scatter(xAxis,ber(rlz,:),markersz,'MarkerEdgeColor',col,'MarkerFaceColor',col,'Parent', AxesMain(1),'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname]); + end + + end + +elseif string(plotJob.plotTypeArg) == "Lines" + + % if ~anynan(ber) + % [xAxis,ber] = interpCurve(xAxis, ber); + % end + + for rlz = 1:size(ber,1) +% + if (string(plotJob.dataStatArg) == "All Channels ;All PMD Realizations")||(string(plotJob.dataStatArg) == "All Channels; mean(PMD Realizations)") + ch = mod(rlz,plotJob.ch); + if ch == 0; ch = plotJob.ch; end + else + ch = plotJob.dataStatArg; + end + + if rlz <= size(ber,1) + + + s = plot3(xAxis,repmat(ch,1,numel(xAxis)),ber(rlz,:),linestyle,'Marker',"o",'MarkerSize',markersz,'MarkerFaceColor',plotJob.color,'LineWidth',linew,'Color',col,'Parent', AxesMain,'HandleVisibility','off'); + + s.DataTipTemplate.Interpreter = "latex"; + s.DataTipTemplate.DataTipRows(1).Label = "Ch: "; + s.DataTipTemplate.DataTipRows(1).Value = repmat(ch,size(ber)); + s.DataTipTemplate.DataTipRows(1); + s.DataTipTemplate.DataTipRows(2) = []; + + + else + + if string(plotJob.dataStatArg) == "Lineplot with quartiles" + [hl,hp] = boundedline(xAxis,ber(rlz,:),([(ber(rlz,:)-lowerq(rlz,:));upperq(rlz,:)-ber(rlz,:)]'),'-*', 'alpha','Color',col,'transparency', 0.2); + hp.LineWidth = 1.2; + + ho = outlinebounds(hl,hp); + set(ho, 'linestyle', ':', 'color', col); + else + + s = plot(xAxis,ber(rlz,:),linestyle,'Marker',"o",'MarkerFaceColor',col,'MarkerSize',markersz,'LineWidth',linew,'Color',col,'Parent', AxesMain,'DisplayName',[plotJob.displayname]); + + s.DataTipTemplate.Interpreter = "latex"; + s.DataTipTemplate.DataTipRows(1).Label = "Ch: "; + s.DataTipTemplate.DataTipRows(1).Value = repmat(string(ch),size(ber)); + s.DataTipTemplate.DataTipRows(1) + s.DataTipTemplate.DataTipRows(2) = []; + end + + end + + end +end + +% Draw FEC Threshold Line +%get x data of first children: +%get all linear Values +if 0 + linear = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,wh.parameter.p_out.values,0,0,499,"symmetric"); + linear = mean(linear,2); + lincurve = findall(AxesMain, 'Type', 'line','DisplayName','Linear Baseline'); + % + if isempty(lincurve) + xdata = AxesMain.Children(1).XData; + hdfec = 3.8e-3.*ones(size(xdata)); + plot(xdata,linear,'-','Marker',"o",'MarkerSize',3,'LineWidth',4,'Color',[0.6400 0.6400 0.6400],'MarkerFaceColor',[0.6400 0.6400 0.6400],'Parent', AxesMain,'DisplayName','Linear Baseline'); + %h = get(gca,'Children'); + %set(gca,'Children',[h(2) h(1)]) + end +end + +feccurve = findall(AxesMain, 'Type', 'line','DisplayName','FEC $3.8*10^{-3}$'); +% +if isempty(feccurve) + xdata = xAxis; + hdfec = 3.8e-3.*ones(size(xdata)); + for ch = 1:plotJob.ch + plot3(xdata,repmat(ch,1,numel(xAxis)),hdfec,':','MarkerSize',4,'Color','black','MarkerFaceColor','black','LineWidth',0.5,'Parent', AxesMain,'DisplayName','FEC $3.8*10^{-3}$','HandleVisibility','off'); + end + %h = get(gca,'Children'); + %set(gca,'Children',[h(2) h(1)]) +end + +% Figure Settings +%title(AxesMain,plotJob.title,"Interpreter","none"); + +xlabel(AxesMain,plotJob.xAxisLabel,"Interpreter","none"); + +ylabel(AxesMain,plotJob.yAxisLabel,"Interpreter","none"); + +set(AxesMain,'zscale','log'); + +grid(AxesMain,'on'); + +grid(AxesMain,'minor'); + +grid minor + +view(AxesMain,[42.0619302949062 23.4176470588235]); +%legend(AxesMain); + +fontsize(AxesMain,8,"points") +fontname(AxesMain,"Arial") + +fig.Position = plotJob.Position; +fig.Units = "centimeters"; +fig.Position = [2 2 8.5 7]; + +set(AxesMain,'TickLabelInterpreter','none') + +set(AxesMain.Legend,'Interpreter','none') +% set(gcf,'Units','centimeters') +% set(gcf,'Position',[2 2 9 4.5]) + +zlim([1e-4,0.3]); + +xlim([min(xAxis),-3]); + + +annotation('textbox', [0.125, 0.32, 0.1, 0.1], 'String', "FEC 3.8e-3","LineStyle","none","FontSize",8,"FontUnits","points") + + +hold off + + +end + + +function vec = removeZeros(vec) + % Find rows that contain only zeros + rows_to_remove = all(vec == 0, 2); + + % Remove rows with only zeros + vec(rows_to_remove, :) = []; +end diff --git a/Classes/Warehouse_class/functions/fwm_plots/plotBerVsZDW.m b/Classes/Warehouse_class/functions/fwm_plots/plotBerVsZDW.m new file mode 100644 index 0000000..bebc5b5 --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/plotBerVsZDW.m @@ -0,0 +1,300 @@ +function plotBerVsZDW(wh,plotJob) + + fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName); + + if isvalid(fig) + figure(fig) + fig = get(fig); + AxesMain = fig.CurrentAxes; + hold on + else + fig = figure('name',char(plotJob.figName)); + AxesMain = gca; + hold on + end + + + col = plotJob.color; + + % we want to fetch all realizations + realization = wh.parameter.realization.values(1:end); + + % get all xAxis values + xAxis = wh.parameter.p_out.values; + + % get all center wavelengths + wavelengths = wh.parameter.center_wavelength.values; + +% totber = NaN(numel(realization),1,numel(xAxis),numel(wavelengths)); +% ber = zeros(numel(realization),plotJob.ch,numel(xAxis),numel(wavelengths)); + + %get BER values for query + for w = 2:numel(wavelengths) + + for xl = 1:numel(xAxis) + + c_wavelen = wavelengths(w); + p_out = xAxis(xl); + + % dim1 : realiz; dim2: channels, dim3: rop, dim4, c_wavelength + temp = removeZeros(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,c_wavelen,plotJob.channelspacing,plotJob.randzdw)); + ber(1:size(temp,1),:,xl,w) = temp; + + end + end + + hdfec = 3.8e-3.*ones(size(xAxis)); + zdw_ = []; + zdw_chann = []; + zdw_tot = []; + cf_tot = []; + cf_ = []; + S = []; + Stot = []; + S_chann = []; + cf_chann = []; + + cnt = 0; + %get fec thresholds + %linear = squeeze(linear); + for c_wavelen = 1:size(ber,4) + for realiz = 1:size(ber,1) + for chann = 1:size(ber,2) + + %finde Schnittpunkt zwischen FEC und BER Kurve + temp_ber = squeeze(ber(realiz,chann,:,c_wavelen)).'; + if ~all(temp_ber == 0) + %nur wenn nicht alles nullen sind + crossing_ch = InterX([hdfec;xAxis],[temp_ber;xAxis]); + else + continue + end + + %Req. FEC Ergebnis einsortieren + if ~isempty(crossing_ch) + if crossing_ch(2) == 0 + print("d") + end + S(realiz,chann,c_wavelen) = crossing_ch(2); + else + S(realiz,chann,c_wavelen) = -1; + cnt = cnt +1; + end + + + end + end + end + + temp_max = -inf; + for i = 1:plotJob.ch + hold on + %S:: 1.dim: realiz; 2.dim: channel; 3.dim: center wavelength + %squeeze a channel: + temp_data = squeeze(S(:,i,:)); + + %remove realizations that have no entry (only zero) + temp_data = removeZeros(temp_data); + + %replace zeros with NAN (e.g. for the wavelengths that have missing realizations) + temp_data(temp_data==0) = NaN; + + %plot required ROP for channel and all realizations that cross the + %FEC limit + scatter(wavelengths,temp_data ,5,plotJob.color,'Marker','.'); + +% %plot mean per channel +% temp_mean = mean(temp_data,'omitnan'); +% hold on +% plot(wavelengths,temp_mean,'Marker','*'); + + %get max overall value + temp_max = max(temp_max,max(temp_data)); + end + + + + %plot mean overall + mean_overall = squeeze(mean(S,2)); + mean_overall(mean_overall==0) = NaN; + %mean_overall(mean_overall==-1) = NaN; + mean_overall=mean(mean_overall,1,'omitnan'); + plot(wavelengths,mean_overall,'Color',plotJob.color); + + %plot max overall + scatter(wavelengths,temp_max ,35,plotJob.color,'Marker','v'); + + + %plot channel positions + hold on + chpos = calcWavelengthPlan(plotJob.ch, 400e9, 1310); + xline(chpos,'LineWidth',2,'Alpha',0.2); + chpos = calcWavelengthPlan(plotJob.ch, 400e9, chpos(4)); + xline(chpos,'LineWidth',2,'Alpha',0.2); + + fig.Position = plotJob.Position; + + ylabel('Penalty in dB'); + xlabel('Wavelength in nm'); + + xlim([min(wavelengths),max(wavelengths) ]); + + grid minor; + set(gca, 'color', 'none'); + legend = []; + + fontsize(AxesMain,8,"points") + + fig.Position = plotJob.Position; + fig.Units = "centimeters"; + fig.Position = [2 2 8.5 7]; + + set(AxesMain,'TickLabelInterpreter','latex') + + set(AxesMain.Legend,'Interpreter','latex') + + + + + + + % + % + % + % + % + % distinct_cf = unique(cf_chann); + % + % for i = 1:length(distinct_cf) + % indices = find(cf_chann==distinct_cf(i)); + % cf(i) = distinct_cf(i); + % worst_fec_cross(i) = max(S_chann(indices)); + % avg_fec_cross(i) = mean(S_chann(indices)); + % end + % + % avg_fec_cross = smooth(avg_fec_cross,5); + % + % figure(224) + % hold on + % scatter(cf_,S,10.*abs(S-mean(S)).*ones(size(S)),'DisplayName',['AVG'],'MarkerEdgeColor',col,'MarkerFaceColor',col,'Marker','.'); + % hold on + % scatter(cf(2:end),worst_fec_cross(2:end),15,'DisplayName',['Worst'],'MarkerEdgeColor',col,'MarkerFaceColor',col,'Marker','.','HandleVisibility','off'); + % plot(cf(2:end),avg_fec_cross(2:end),'DisplayName',['AVG'],'LineWidth',1,'LineStyle','-','Color',col,'Marker','none','MarkerFaceColor',col,'MarkerSize',2); + % plot(cf(2:end),worst_fec_cross(2:end),'DisplayName',['AVG'],'LineWidth',0.5,'LineStyle','-','Color',col,'Marker','none','MarkerFaceColor',col,'MarkerSize',2); + % ghzgrid = hz2nm(nm2hz(1310)+[0:1:20].*400e9); + % set(gca,'xtick',sort(ghzgrid)) + % xlim([min(cf(cf~=0)), 1310.1]); + % + % + % % With matlab internal errorbar function... + % figure(221) + % %plot(cf,avg_fec_cross,'DisplayName',['AVG'],'LineWidth',1,'Color',col,'Marker','none','MarkerFaceColor',col,'MarkerSize',2); + % hold on + % %plot(cf_tot,min(S_chann(1:length(cf_tot),:),[],2),'DisplayName',['AVG'],'LineWidth',1,'LineStyle',':','Color',col,'Marker','^','MarkerFaceColor',col,'MarkerSize',2); + % plot(cf,worst_fec_cross,'DisplayName',['AVG'],'LineWidth',2,'LineStyle','-','Color',col,'Marker','none','MarkerFaceColor',col,'MarkerSize',2); + % ghzgrid = hz2nm(nm2hz(1310)+[0:1:20].*400e9); + % set(gca,'xtick',sort(ghzgrid)) + % xlim([1302, 1310.1]); + % xline(ghzgrid,'LineStyle',':','Color',[.7 .7 .7]); + + + %with + % Stot = movmean(Stot,5); + % figure(222) + % [hl,hp] = boundedline(cf_tot,Stot,[(Stot'-min(S_chann(1:length(cf_tot),:),[],2)),(max(S_chann(1:length(cf_tot),:),[],2)-Stot')], 'alpha','Color',col,'transparency', 0.05); + % ho = outlinebounds(hl,hp); + % set(ho, 'linestyle', ':', 'color', col, 'marker', '.','linewidth',0.5); + % hold on + % + % ghzgrid = hz2nm(nm2hz(1310)+[0:1:12].*400e9); + % xline(ghzgrid); + + + + + %plot the total ber + + % %figure(22); + % hold on; + % b= movmean(Stot,3); + % plot(AxesMain,cf_tot,b,'LineWidth',2,'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname],'Color',col,'Marker','o'); + % hold on + + + + + + %scatter(AxesMain,zdw_,S,'Marker','+','MarkerEdgeColor',col,'MarkerFaceAlpha',0.4,'MarkerEdgeAlpha',0.4,'LineWidth',0.5,'Parent', AxesMain(1),'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname]); + %ylim(AxesMain,[-9.3 -7]); + + % for i = 1:numel(chp) + % hold on + % xline(AxesMain,chp(i),'Color',colr(i,:),'DisplayName',['CH: ', num2str(i)],'LineWidth',1.5); + % hold off + % end + + + + + % + % a = movmean(sortrows([zdw_; S]'),10,'Endpoints','discard'); + % + % sorted = sortrows([zdw_; S]'); + % figure(2) + % scatter(sorted(:,1),sorted(:,2)) + % + % ber_sorted = sort(S); + % mean(ber_sorted); + % std(ber_sorted); + % z1 = []; + % penalty = mean(ber_sorted):0.01:mean(ber_sorted)+2; + % for i = 1:length(penalty) + % l = penalty(i); + % if i == 1 + % z1 = [z1 sum(ber_sorted(1,:)l-0.01) / length(ber_sorted) ]; + % end + % + % end + % + % penalty_higherthan = 0.5; + % probability = sum(z1(find(penalty>mean(ber_sorted)+penalty_higherthan))); + % disp(['A penalty of more than 1dB has a probability of: ', num2str(probability)]); + % % + % stem(AxesMain,penalty,z1,"filled",'Marker','o','MarkerSize',2,'Color',col); + % + % + % [f1,x1]=ecdf(S(end,:)); + % %figure(23);plot(AxesMain,x1,f1,'r','LineWidth',3, 'Color',col); + % + % %plot(AxesMain,a(:,1),a(:,2),'Color',col+1,'Parent', AxesMain(1)); + % + % % histogram(AxesMain,S,1000,'EdgeColor','none','FaceAlpha',0.4); + % + % + % + % % + % %scatter(AxesMain,zdw_,S,'Marker','+','MarkerEdgeColor',col,'MarkerFaceAlpha',0.4,'MarkerEdgeAlpha',0.4,'LineWidth',0.5,'Parent', AxesMain(1),'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname]); + % hold on + % %scatter(AxesMain,zdw_chann(:,1),mean(S_chann,2),'Marker','diamond','MarkerEdgeColor',col,'MarkerFaceAlpha',0.6,'LineWidth',7,'MarkerFaceColor',col,'Parent', AxesMain(1),'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname]); + % hold off + % % + % % for rlz = 1:size(S,2) + % % zdwval = zdw_(rlz); + % % feccrossing = S(rlz); + % % scatter(AxesMain,zdwval,feccrossing,10,'MarkerEdgeColor',col,'MarkerFaceColor',col,'Parent', AxesMain(1),'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname]); + % % end + % + % xline([1302.03471732576,1304.30061112288,1306.57440520904,1308.85614097425,1311.14586009814,1313.44360455251,1315.74941660391,1318.06333881618]); + + +end + +function vec = removeZeros(vec) + % Find rows that contain only zeros + rows_to_remove = all(vec == 0, 2); + + % Remove rows with only zeros + vec(rows_to_remove, :) = []; +end \ No newline at end of file diff --git a/Classes/Warehouse_class/functions/fwm_plots/plotBerVsZdwFailureRate.m b/Classes/Warehouse_class/functions/fwm_plots/plotBerVsZdwFailureRate.m new file mode 100644 index 0000000..480641f --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/plotBerVsZdwFailureRate.m @@ -0,0 +1,157 @@ +function plotBerVsZdwFailureRate(wh,plotJob) + + fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName); + + if isvalid(fig) + figure(fig) + fig = get(fig); + AxesMain = fig.CurrentAxes; + hold on + else + fig = figure('name',char(plotJob.figName)); + AxesMain = gca; + hold on + end + + + col = plotJob.color; + + % we want to fetch all realizations + realization = wh.parameter.realization.values(1:end); + + % get all xAxis values + xAxis = wh.parameter.p_out.values; + + % get all center wavelengths + wavelengths = wh.parameter.center_wavelength.values; + %wavelengths = wavelengths(2:end); +% totber = NaN(numel(realization),1,numel(xAxis),numel(wavelengths)); +% ber = zeros(numel(realization),plotJob.ch,numel(xAxis),numel(wavelengths)); + + %get BER values for query + for w = 1:numel(wavelengths) + + for xl = 1:numel(xAxis) + + c_wavelen = wavelengths(w); + p_out = xAxis(xl); + + % dim1 : realiz; dim2: channels, dim3: rop, dim4, c_wavelength + temp = removeZeros(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,c_wavelen,plotJob.channelspacing,plotJob.randzdw)); + ber(1:size(temp,1),:,xl,w) = temp; + + end + end + + hdfec = 3.8e-3.*ones(size(xAxis)); + zdw_ = []; + zdw_chann = []; + zdw_tot = []; + cf_tot = []; + cf_ = []; + S = []; + Stot = []; + S_chann = []; + cf_chann = []; + + cnt = 0; + %get fec thresholds + %linear = squeeze(linear); + for c_wavelen = 1:size(ber,4) + for realiz = 1:size(ber,1) + for chann = 1:size(ber,2) + + %finde Schnittpunkt zwischen FEC und BER Kurve + temp_ber = squeeze(ber(realiz,chann,:,c_wavelen)).'; + if ~all(temp_ber == 0) + %nur wenn nicht alles nullen sind + crossing_ch = InterX([hdfec;xAxis],[temp_ber;xAxis]); + else + continue + end + + %Req. FEC Ergebnis einsortieren + if ~isempty(crossing_ch) + if crossing_ch(2) == 0 + print("d") + end + S(realiz,chann,c_wavelen) = crossing_ch(2); + else + S(realiz,chann,c_wavelen) = -1; + cnt = cnt +1; + end + + + end + end + end + + temp_max = -inf; + sum_FEC_not_crossed=[]; + sum_FEC_crossed=[]; + + threshold = plotJob.p_in - 10; + + for i = 1:plotJob.ch + hold on + %S:: 1.dim: realiz; 2.dim: channel; 3.dim: center wavelength + %squeeze a channel: + temp_data = squeeze(S(:,i,:)); + + %remove realizations that have no entry (only zero) + temp_data = removeZeros(temp_data); + + %replace zeros with NAN (e.g. for the wavelengths that have missing realizations) + temp_data(temp_data==0) = NaN; + + %for current channel + FEC_crossed = temp_data < threshold & ~isnan(temp_data); + FEC_not_crossed = temp_data >= threshold & ~isnan(temp_data); + + %sum over channels for overall picture + sum_FEC_not_crossed(i,:) = sum(FEC_not_crossed); + sum_FEC_crossed(i,:) = sum(FEC_crossed); + + failure_rate_channelwise(i,:) = sum_FEC_not_crossed(i,:)./ ( sum_FEC_crossed(i,:) + sum_FEC_not_crossed(i,:)); + + end + + failure_rate_total = sum(sum_FEC_not_crossed,1) ./ ( sum(sum_FEC_crossed,1) + sum(sum_FEC_not_crossed,1) ); + % plot failure rate (nbetween 0 and 1) + + plot(wavelengths,failure_rate_total,'Color',plotJob.color,'LineWidth',1,'LineStyle',plotJob.plotTypeArg,'Marker','x','MarkerSize',5,'MarkerFaceColor',plotJob.color,'DisplayName',plotJob.displayname); + + %plot max overall + %scatter(wavelengths,failure_rate_channelwise ,35,plotJob.color,'Marker','.'); + + ylabel('Failure Rate of Link'); + xlabel('Wavelength in nm'); + + xlim([min(wavelengths),max(wavelengths) ]); + ylim([0,1]); + + grid minor; + set(gca, 'color', 'none'); + legend = []; + +% fontsize(AxesMain,8,"points") + + fig.Position = plotJob.Position; + fig.Units = "centimeters"; + fig.Position = [0 0 12 5 7]; + + try + set(AxesMain,'TickLabelInterpreter','latex') + + set(AxesMain.Legend,'Interpreter','latex') + end + +end + +function vec = removeZeros(vec) + % Find rows that contain only zeros + rows_to_remove = all(vec == 0, 2); + + % Remove rows with only zeros + vec(rows_to_remove, :) = []; +end \ No newline at end of file diff --git a/Classes/Warehouse_class/functions/fwm_plots/plotChannelSpacingAna.m b/Classes/Warehouse_class/functions/fwm_plots/plotChannelSpacingAna.m new file mode 100644 index 0000000..44febc6 --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/plotChannelSpacingAna.m @@ -0,0 +1,51 @@ +function plotChannelSpacingAna(wh,plotJob) + +xAxis = wh.parameter.p_out.values; + + +realization = wh.parameter.realization.values(1:end); + +channelsp = wh.parameter.channelspacing.values(1:end); +channelsp = [200 400].*1e9; +for ch = 1:2 + channspacing = channelsp(ch); + for xl = 1:numel(xAxis) + p_out = xAxis(xl); + + curber = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.numchannels,channspacing); + curzdw = wh.getStoValue('zdw',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.numchannels,channspacing); + + ber(1:size(curber,1),1:size(curber,2),xl) = curber; + zdw(1:size(curber,1),1,xl) = curzdw; + + end + + ber = squeeze(mean(ber,1)); + zdw = squeeze(mean(zdw,1)); + + hdfec = 3.8e-3.*ones(size(xAxis)); + S = []; + wavelength={}; + + zdw_ = []; + zdw_chann = []; + zdw_tot = []; + S = []; + S_chann = []; + wl = round([1302.03471732576,1304.30061112288,1306.57440520904,1308.85614097425,1311.14586009814,1313.44360455251,1315.74941660391,1318.06333881618],2); + wl = 1:16; + wl = [1.2930 1.2953 1.2975 1.2998 1.3020 1.3043 1.3066 1.3089 1.3111 1.3134 1.3157 1.3181 1.3204 1.3227 1.3251 1.3274]; + + + a = InterX([hdfec(:)';xAxis],[mean(ber,1);xAxis]); + if ~isempty(a) + s(ch) = a(2); + else + s(ch) = NaN; + end +end + +figure(2224) +hold on +plot(channelsp,s,'LineWidth',1,'Color',plotJob.color,'Marker','o'); + diff --git a/Classes/Warehouse_class/functions/fwm_plots/plotCurve.m b/Classes/Warehouse_class/functions/fwm_plots/plotCurve.m new file mode 100644 index 0000000..a7134c9 --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/plotCurve.m @@ -0,0 +1,278 @@ +function plotCurve(wh,plotJob) +%PLOTCURVE Summary of this function goes here +% Detailed explanation goes here + +fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName); + +if isvalid(fig) + figure(fig) + fig = get(fig); + AxesMain = fig.CurrentAxes; + hold on +else + fig = figure('name',char(plotJob.figName)); + AxesMain = gca; + hold on +end + +col = plotJob.color; + +% we want to fetch all realizations +if plotJob.pmd == 0 + realization = 499; + % realization = 0:7; +else + realization = wh.parameter.realization.values(1:end); +end +realization = wh.parameter.realization.values(1:end); +% get all xAxis values +xAxis = wh.parameter.p_out.values; + +% Fetch Data from Warehouse +for xl = 1:numel(xAxis) + p_out = xAxis(xl); + if string(plotJob.dataStatArg) == "Worst" + temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw); + temp = removeZeros(temp); + ber(xl) = max(temp,[],'all'); + linew = 1.0; + markersz = 3; + linestyle = '-'; + + elseif string(plotJob.dataStatArg) == "AVG" + temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw); + temp = removeZeros(temp); + ber(xl) = mean(temp,'all'); + linew = 1.0; + markersz = 3; + linestyle = '--'; + + elseif string(plotJob.dataStatArg) == "Best" + temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw); + temp = removeZeros(temp); + ber(xl) = min(temp,[],'all'); + linew = 1.0; + markersz = 3; + linestyle = '-'; + + elseif string(plotJob.dataStatArg) == "All Channels; mean(PMD Realizations)" + temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw); + temp = removeZeros(temp); + + ber(:,xl) = mean(temp,1,"omitnan").'; + + linew = 1; + markersz = 1; + linestyle = '-'; + + elseif string(plotJob.dataStatArg) == "Lineplot with quartiles" + + temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw); + temp = removeZeros(temp); + ber(xl) = mean(temp,"all","omitnan").'; + + + upperq(xl) = quantile(temp,0.99,"all"); + lowerq(xl) = quantile(temp,0.04,"all"); + + if lowerq(xl) == 0 + lowerq(xl) = lowerq(xl-1); + end + % upperq(xl) = 0.5*std(tmp,1,'all','omitnan'); + % lowerq(xl) = 0.5*std(tmp,1,'all','omitnan'); + + % upperq(xl) = max(dataNoNans); + % lowerq(xl) = min(dataNoNans); + + % upperq(xl) = mean(tmp,"all","omitnan") + 1.96 * (std(tmp,1,'all','omitnan')/sqrt(numel(tmp))); + % lowerq(xl) = mean(tmp,"all","omitnan") - 1.96 * (std(tmp,1,'all','omitnan')/sqrt(numel(tmp))); + + linew = 1; + markersz = 1; + linestyle = '-'; + + elseif string(plotJob.dataStatArg) == "All Channels ;All PMD Realizations" + + raw_fetch = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw).'; + tmp = reshape(raw_fetch,[],1); + ber(1:size(tmp,1),xl) = tmp; + + ber_per_chann(:,:,xl) = raw_fetch; + + linew = 0.3; + markersz = 2; + linestyle = ':'; + + end +end + + +%routine to remove total outliers (here those wehere the rop curve has a mean BER greater than 0.1) +% ber_per_chann_clean = NaN(size(ber_per_chann)); +% for ch = 1:size(ber_per_chann,1) +% bla = squeeze(ber_per_chann(ch,:,:)); +% ber_per_chann(ch,find(mean(bla,2)>0.25),:) = NaN; +% cleaned = rmoutliers(bla,"mean",'ThresholdFactor',2); +% +% ber_per_chann_clean(ch,1:size(cleaned,1),1:size(cleaned,2)) = cleaned; +% +% end +% +% ber = []; +% for rop = 1:size(ber_per_chann,3) +% temp = squeeze(ber_per_chann_clean(:,:,rop)); +% ber(rop) = mean(temp,"all","omitnan").'; +% upperq(rop) = quantile(temp,0.9,"all"); +% lowerq(rop) = quantile(temp,0.1,"all"); +% end + + + + +xAxis = xAxis; + +% Plot Data +if string(plotJob.plotTypeArg) == "Scatter" + for rlz = 1:size(ber,1) + + if rlz < size(ber,1) + scatter(xAxis,ber(rlz,:),markersz,'MarkerEdgeColor',col,'MarkerFaceColor',col,'Parent', AxesMain(1),'HandleVisibility','off'); + else + scatter(xAxis,ber(rlz,:),markersz,'MarkerEdgeColor',col,'MarkerFaceColor',col,'Parent', AxesMain(1),'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname]); + end + + end + +elseif string(plotJob.plotTypeArg) == "Lines" + + % if ~anynan(ber) + % [xAxis,ber] = interpCurve(xAxis, ber); + % end + + for rlz = 1:size(ber,1) +% + if (string(plotJob.dataStatArg) == "All Channels ;All PMD Realizations")||(string(plotJob.dataStatArg) == "All Channels; mean(PMD Realizations)") + ch = mod(rlz,plotJob.ch); + if ch == 0; ch = plotJob.ch; end + else + ch = plotJob.dataStatArg; + end + + if rlz < size(ber,1) + + + s = plot(xAxis,ber(rlz,:),linestyle,'Marker',"none",'MarkerSize',markersz,'LineWidth',linew,'Color',col,'Parent', AxesMain,'HandleVisibility','off'); + + s.DataTipTemplate.Interpreter = "latex"; + s.DataTipTemplate.DataTipRows(1).Label = "Ch: "; + s.DataTipTemplate.DataTipRows(1).Value = repmat(ch,size(ber)); + s.DataTipTemplate.DataTipRows(1); + s.DataTipTemplate.DataTipRows(2) = []; + + + else + + if string(plotJob.dataStatArg) == "Lineplot with quartiles" + [hl,hp] = boundedline(xAxis,ber(rlz,:),([(ber(rlz,:)-lowerq(rlz,:));upperq(rlz,:)-ber(rlz,:)]'),'-o','alpha','Color',col,'transparency', 0.1,'linewidth',0.7); + hl.MarkerFaceColor = col; + hl.MarkerSize = 3; + set(hp,'HandleVisibility','off'); + %hp.LineWidth = 1.2; + + ho = outlinebounds(hl,hp); + set(ho, 'linestyle', ':', 'color', col,'Linewidth',0.5); + set(ho,'HandleVisibility','off'); + else + + s = plot(xAxis,ber(rlz,:),linestyle,'Marker',"o",'MarkerFaceColor',col,'MarkerSize',markersz,'LineWidth',linew,'Color',col,'Parent', AxesMain,'DisplayName',[plotJob.displayname]); + + s.DataTipTemplate.Interpreter = "latex"; + s.DataTipTemplate.DataTipRows(1).Label = "Ch: "; + s.DataTipTemplate.DataTipRows(1).Value = repmat(string(ch),size(ber)); + s.DataTipTemplate.DataTipRows(1) + s.DataTipTemplate.DataTipRows(2) = []; + end + + end + + end +end + +% Draw FEC Threshold Line +%get x data of first children: +%get all linear Values +if 0 + linear = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,wh.parameter.p_out.values,0,0,499,"symmetric"); + linear = mean(linear,2); + lincurve = findall(AxesMain, 'Type', 'line','DisplayName','Linear Baseline'); + % + if isempty(lincurve) + xdata = AxesMain.Children(1).XData; + hdfec = 3.8e-3.*ones(size(xdata)); + plot(xdata,linear,'-','Marker',"o",'MarkerSize',3,'LineWidth',4,'Color',[0.6400 0.6400 0.6400],'MarkerFaceColor',[0.6400 0.6400 0.6400],'Parent', AxesMain,'DisplayName','Linear Baseline'); + %h = get(gca,'Children'); + %set(gca,'Children',[h(2) h(1)]) + end +end + +feccurve = findall(AxesMain, 'Type', 'line','DisplayName','FEC $3.8*10^{-3}$'); +% +if isempty(feccurve) + xdata = AxesMain.Children(1).XData; + hdfec = 3.8e-3.*ones(size(xdata)); + plot(xdata,hdfec,'--','MarkerSize',4,'Color','black','MarkerFaceColor','black','LineWidth',1,'Parent', AxesMain,'DisplayName','FEC $3.8*10^{-3}$','HandleVisibility','off'); + %h = get(gca,'Children'); + %set(gca,'Children',[h(2) h(1)]) +end + +% Figure Settings +%title(AxesMain,plotJob.title,"Interpreter","none"); + +xlabel(AxesMain,plotJob.xAxisLabel,"Interpreter","latex"); + +ylabel(AxesMain,plotJob.yAxisLabel,"Interpreter","latex"); + +set(AxesMain,'yscale','log'); + +grid(AxesMain,'on'); + +grid(AxesMain,'minor'); + +grid minor + +%legend(AxesMain); + +fontsize(AxesMain,8,"points") +% fontname(AxesMain,"Arial") + +fig.Position = plotJob.Position; +fig.Units = "centimeters"; +fig.Position = [2 2 8.5 7]; + +set(AxesMain,'TickLabelInterpreter','latex') + +set(AxesMain.Legend,'Interpreter','latex') +% set(gcf,'Units','centimeters') +% set(gcf,'Position',[2 2 9 4.5]) + +ylim([1e-5,0.3]); + +xlim([min(xAxis),-3]); + + +% annotation('textbox', [0.125, 0.32, 0.1, 0.1], 'String', "FEC 3.8e-3","LineStyle","none","FontSize",8,"FontUnits","points") + + +hold off + + +end + + +function vec = removeZeros(vec) + % Find rows that contain only zeros + rows_to_remove = all(vec == 0, 2); + + % Remove rows with only zeros + vec(rows_to_remove, :) = []; +end diff --git a/Classes/Warehouse_class/functions/fwm_plots/plotHistogram.m b/Classes/Warehouse_class/functions/fwm_plots/plotHistogram.m new file mode 100644 index 0000000..023a3a1 --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/plotHistogram.m @@ -0,0 +1,151 @@ +function plotHistogram(wh,plotJob) + + +fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName); + +if isvalid(fig) + figure(fig) + fig = get(fig); + AxesMain = fig.CurrentAxes; + hold on +else + fig = figure('name',char(plotJob.figName)); + AxesMain = gca; + hold on +end + +col = plotJob.color; + +% we want to fetch all realizations +realization = wh.parameter.realization.values(1:end); + +% get all xAxis values +xAxis = wh.parameter.p_out.values; + + +% Fetch Data from Warehouse +for xl = 1:numel(xAxis) + p_out = xAxis(xl); + if string(plotJob.dataStatArg) == "Worst" + ber(xl) = max(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,string(plotJob.channelplan),string(plotJob.dsp)),[],'all'); + linew = 2.0; + markersz = 3; + linestyle = '-'; + elseif string(plotJob.dataStatArg) == "AVG" + ber(xl) = mean(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,string(plotJob.channelplan),string(plotJob.dsp)),'all'); + linew = 2.0; + markersz = 3; + linestyle = ':'; + elseif string(plotJob.dataStatArg) == "Best" + ber(xl) = min(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,string(plotJob.channelplan),string(plotJob.dsp)),[],'all'); + linew = 2.0; + markersz = 3; + linestyle = ':'; + + elseif string(plotJob.dataStatArg) == "All Channels; mean(PMD Realizations)" + tmp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,string(plotJob.channelplan),string(plotJob.dsp)); + + if numel(tmp(tmp==0)) ~= 0 + disp('Removed all zero values!'); + tmp(tmp==0) = NaN; + end + ber(:,xl) = mean(tmp,1,"omitnan").'; + + linew = 0.7; + markersz = 2; + linestyle = '-'; + + elseif string(plotJob.dataStatArg) == "All Channels ;All PMD Realizations" + + + + tmp = reshape(wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,string(plotJob.channelplan),string(plotJob.dsp)).',[],1); + ber(1:size(tmp,1),xl) = tmp; + + linew = 0.3; + markersz = 2; + linestyle = ':'; + end +end + +disp('Removed all zero values!'); +ber(ber==0) = NaN; +if ~anynan(ber) + [xAxis,ber] = interpCurve(xAxis, ber); +end + +% plot FEC Crossing as histogram + +hdfec = 3.8e-3.*ones(size(xAxis)); +S = []; +for i = 1:size(ber,1) + a = InterX([hdfec;xAxis],[ber(i,:);xAxis]); + if ~isempty(a) + S(:,i) = a; + end +end + + +%% SUB 1 +AxesMain = subplot(2,1,1); + +hold on + +if ~isempty(S) + histogram(S(end,:),300,'Normalization','probability','FaceColor',col,'EdgeColor',col,'Parent',AxesMain,'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname],'FaceAlpha',0.4,'EdgeAlpha',0.4); +end + +xlim([-3 ,9 ]); +ylim([0 .10]); + +% Figure Settings +title(AxesMain,['$P_{in}:$ ',num2str(plotJob.p_in-9), 'dBm; L : ', num2str(plotJob.l),'Km'],'Interpreter','latex'); + +xlabel(AxesMain,plotJob.xAxisLabel,'Interpreter','latex'); + +ylabel('PDF') + +grid(AxesMain,'on'); + +grid(AxesMain,'minor'); + +%legend(AxesMain,'Interpreter','latex'); + +fontsize(AxesMain,24,"pixels") + +hold off + + +%% SUB 2 +AxesMain = subplot(2,1,2); + +if ~isempty(S) +hold on + +[f1,x1]=ecdf(S(end,:)); + plot(x1,f1,'r','LineWidth',3, 'Color',col,'DisplayName',[plotJob.dataStatArg,' Ch.: ', plotJob.displayname]); +end +xlim([-3 ,9 ]); +ylim([0 1]); +% Figure Settings +title(AxesMain,['$P_{in}:$ ',num2str(plotJob.p_in-9), 'dBm; L : ', num2str(plotJob.l),'Km'],'Interpreter','latex'); + +xlabel(AxesMain,plotJob.xAxisLabel,'Interpreter','latex'); + +ylabel('CDF') + +grid(AxesMain,'on'); + +grid(AxesMain,'minor'); + +fontsize(AxesMain,24,"pixels") + +%legend(AxesMain,'Interpreter','latex'); + +fig.Position = plotJob.Position; + +hold off + + + +end \ No newline at end of file diff --git a/Classes/Warehouse_class/functions/fwm_plots/plotViolin.m b/Classes/Warehouse_class/functions/fwm_plots/plotViolin.m new file mode 100644 index 0000000..2453744 --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/plotViolin.m @@ -0,0 +1,206 @@ +function plotViolin(wh,plotJob) + +fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName); + +if isvalid(fig) + figure(fig) + fig = get(fig); + AxesMain = fig.CurrentAxes; + hold on +else + fig = figure('name',char(plotJob.figName)); + AxesMain = gca; + hold on +end + +%% Violin +col = plotJob.color; + +% we want to fetch all realizations +if plotJob.pmd == 0 + realization = 1; +else + realization = wh.parameter.realization.values(1:end); +end + +%realization = 0:8; + +% get all xAxis values +xAxis = wh.parameter.p_out.values; + +% ber = NaN(500,16,10); +% zdw = NaN(500,1,10); + +%get BER values for query +for xl = 1:numel(xAxis) + p_out = xAxis(xl); + + curber = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw); + curber= removeZeros(curber); + ber(1:size(curber,1),1:size(curber,2),xl) = curber; + +end + +% to remove outliers set the percentile range +% a = squeeze(mean(ber,2)); +% out = isoutlier(mean(a,2),"percentiles",[0 100]); +% ber = ber(~out,:,:); +% disp(sum(out)); + +hdfec = 3.8e-3.*ones(size(xAxis)); +S = []; +wavelength={}; + + +S = []; +S_chann = []; + +wl = calcWavelengthPlan(plotJob.ch,plotJob.channelspacing,1310); +%get fec thresholds +% [C,ia,ib] =intersect(linx,xAxis); +% linear = squeeze(linear); + +%ber(ber==0) = NaN; +S_chann_no_crossing = zeros(1,plotJob.ch); +for chann = 1:size(ber,2) + + for realiz = 1:size(ber,1) + + ber_series = squeeze(ber(realiz,chann,:)).'; + if mean(ber_series) > 0.1 + continue + end + + a = InterX([hdfec(:)';xAxis],[ber_series;xAxis]); + + if ~isempty(a) + S_chann(realiz,chann) = a(2); +% if a(2) > -7 && string(plotJob.pol) == "copolarized" +% continue +% end + S(end+1) = a(2); + wavelength{end+1} = num2str(wl(chann)); + else + S(end+1) = 0; + wavelength{end+1} = num2str(wl(chann)); + S_chann_no_crossing(realiz,chann) = 1; + S_chann(realiz,chann) = -1; + end + + end +end + +threshold = -6; +FEC_crossed = sum(S_chann < threshold & ~isnan(S_chann),1); +FEC_not_crossed = sum(S_chann >= threshold & ~isnan(S_chann),1); +failure_rate = FEC_not_crossed ./ (FEC_crossed + FEC_not_crossed) ; + + +S_chann(S_chann==0) = NaN; + +total_avg = mean(S_chann,"all","omitnan"); + +%figure(2024) +%C = flip(cbrewer2('Spectral',8)); +if numel(S) <= numel(wl) + vs = scatter(1:numel(S),S,50,'o','MarkerEdgeColor','black','MarkerFaceColor',plotJob.color,'LineWidth',1,'HandleVisibility','off'); + %vs = scatter(1,mean(S),50,'o','MarkerEdgeColor','black','MarkerFaceColor',plotJob.color,'LineWidth',1); + +else + + vs = violinplot(S,wavelength,... + 'ViolinColor',plotJob.color,... + 'ViolinAlpha',0.1,... + 'MarkerSize',1,... + 'ShowMedian',false,... + 'EdgeColor',plotJob.color,... + 'ShowWhiskers',false,... + 'ShowData',false,... + 'ShowBox',false,... + 'Bandwidth',0.051 ... + ); + + + hold on + + partly_failed = boolean(ceil(failure_rate)); + + avg = mean(S_chann,1,"omitnan"); + + notfailed = ~partly_failed .* avg; + notfailed(notfailed==0) = NaN; + scatter(1:size(S_chann,2),notfailed,10,'Marker','x','MarkerEdgeColor','black','MarkerFaceColor',plotJob.color,'LineWidth',0.5,'HandleVisibility','off'); + + + hold on + partly_failed = partly_failed.*avg; + partly_failed(partly_failed==0) = NaN; + s=scatter(1:numel(failure_rate),partly_failed,10,'Marker','x','LineWidth',0.5,'HandleVisibility','off','MarkerEdgeColor','red'); + s.DataTipTemplate.Interpreter = "latex"; + s.DataTipTemplate.DataTipRows(1).Label = "Fail Rate: "; + s.DataTipTemplate.DataTipRows(1).Value = failure_rate; + s.DataTipTemplate.DataTipRows(2) = []; + hold off +end + +% ax = gca; +% ax.XTicks + +hold on +yline(total_avg,'LineWidth',1,'LineStyle','--','DisplayName','System Avg.') +fig.Position = plotJob.Position; + +xticklabels(1:16); +ylabel('Penalty in dB'); +xlabel('Channel Number'); +ylim([-9.3,-3]); +xlim([0,plotJob.ch+1]); +grid minor; +set(gca, 'color', 'none'); +legend = []; + +fontsize(AxesMain,8,"points") + +fig.Position = plotJob.Position; +fig.Units = "centimeters"; +fig.Position = [2 2 8.5 7]; + +set(AxesMain,'TickLabelInterpreter','latex') + +set(AxesMain.Legend,'Interpreter','latex') + + + +if 0 +ber_sorted = sort(S); + +z1 = []; +penalty = mean(ber_sorted):0.01:mean(ber_sorted)+2; + +for i = 1:length(penalty) + l = penalty(i); + if i == 1 + z1 = [z1 sum(ber_sorted(1,:)l-0.01) / length(ber_sorted) ]; + end + +end + +penalty_higherthan = 0.5; +probability = sum(z1(find(penalty>mean(ber_sorted)+penalty_higherthan))); +disp(['A penalty of more than 0.5 dB has a probability of: ', num2str(probability)]); +end + +end + + +function vec = removeZeros(vec) + % Find rows that contain only zeros + rows_to_remove = all(vec == 0, 2); + + % Remove rows with only zeros + vec(rows_to_remove, :) = []; + + +end \ No newline at end of file diff --git a/Classes/Warehouse_class/functions/fwm_plots/plot_ber_distribution.m b/Classes/Warehouse_class/functions/fwm_plots/plot_ber_distribution.m new file mode 100644 index 0000000..1d43674 --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/plot_ber_distribution.m @@ -0,0 +1,58 @@ + +%automate plots +% [file, path] = uigetfile("C:\Users\Silas\Documents\MATLAB\Raw_Cluster_Simulations\session_januar24\wh_complete_at_1310.mat"); +% wh = load([path filesep file]); +% wh = wh.wh; + +plotJob = struct(); + +plotJob.l = 10; +plotJob.ch = 16; +plotJob.d = 3; +plotJob.sgm = 1; +plotJob.pol = "copolarized"; +plotJob.p_in = 3; +plotJob.gamma = 0.0023; +plotJob.pmd = 0.1; +plotJob.channelspacing = 400e9; +plotJob.randzdw = 0; + +ber_per_chann = []; +% get all xAxis values +xAxis = wh.parameter.p_out.values; +realization = wh.parameter.realization.values(1:end); +% Fetch Data from Warehouse +for xl = 1:numel(xAxis) + p_out = xAxis(xl); + raw_fetch = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw).'; + + ber_per_chann(:,:,xl) = raw_fetch; +end +%% + +figure(2023); + +for ch = [1,floor(plotJob.ch/2),ceil(plotJob.ch/2)+1,plotJob.ch] %1:15:size(ber_per_chann,1) + + for p = 5%1:size(ber_per_chann,3) + + % Extract data for the current row + row_data = squeeze(ber_per_chann(ch,:,p)); + [f, xi] = ksdensity(row_data); + % Identify the peak point + [max_density, max_index] = max(f); + peak_x = xi(max_index); + + end + % Create a histogram plot for the current row with a unique color + plot(xi, f, 'LineWidth', 2, 'DisplayName', ['Ch. ', num2str(ch)],'LineStyle','--'); + %histogram(row_data,100, 'DisplayName', ['Channel ', num2str(ch)], 'EdgeColor', 'none'); + + hold on; % Hold the plot for the next iteration + text(peak_x, max_density, ['Ch ', num2str(ch)], 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'left'); +end + + +legend show + +%% \ No newline at end of file diff --git a/Classes/Warehouse_class/functions/interpCurve.m b/Classes/Warehouse_class/functions/interpCurve.m new file mode 100644 index 0000000..5378f82 --- /dev/null +++ b/Classes/Warehouse_class/functions/interpCurve.m @@ -0,0 +1,9 @@ +function [fitx,fity] = interpCurve(inpx, inpy) + + fitx = linspace(inpx(1),inpx(end),1000); + %x = repmat(inpx, size(inpy) ./ size(inpx)); + fity = zeros(size(inpy,1),1000); + for line = 1:size(inpy,1) + fity(line,:) = interp1( inpx , inpy(line,:) , fitx ,"makima"); + end +end \ No newline at end of file diff --git a/Classes/Warehouse_class/functions/phase_predist_plots/build_wh.m b/Classes/Warehouse_class/functions/phase_predist_plots/build_wh.m new file mode 100644 index 0000000..b8c3981 --- /dev/null +++ b/Classes/Warehouse_class/functions/phase_predist_plots/build_wh.m @@ -0,0 +1,129 @@ +% Script, that shows the data management routine :-) + + + +% 1) Define all your parameters, best practice directly constructs a +% structure + +params = struct; + +params.D =[0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 ]; +params.a =[-2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 ]; +params.chirp =[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 ]; +params.rop =[-15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 ]; + +%wh = warehouse :-) +wh = DataStorage(params); + +wh.showInfo; + +wh.addStorage("ber"); + + +%2) Simulate a bunch of data - TO BE IMPLEMENTED HERE - for now use scripts +%from Sebastian + +%3) Once the simulation folder is around, specifiy path and analyze dirs + +path = uigetdir('C:\Users\Silas\Documents\MATLAB\move_it_model_collection\phase_predist_imdd\save2km'); + +allMat = getAllFilesInFolder(path,'.mat'); +allErr = getAllFilesInFolder(path,'.err'); + +%allMat = dir([path filesep '*.mat']); + +%allErr = dir([path filesep '*.err']); + +if numel(allMat) == 0 + warning('You defined an empty folder. Could not locate any .mat file.') +else + fprintf('%-20s', 'Err Files:'); fprintf('%-12s', num2str(numel(allErr))); fprintf('\n'); + fprintf('%-20s', 'Mat Files:'); fprintf('%-12s', num2str(numel(allMat))); fprintf('\n'); + fprintf('%-20s', 'Missing Mat Files:'); fprintf('%-12s', num2str(numel(allErr)-numel(allMat))); fprintf('\n'); +end + +%4) Now load that data + +f = waitbar(0,'Please wait...'); +cnt = 0; + +for num = 1:numel(allMat) + + fileName = allMat(num).name; + fileFolder = allMat(num).path; + fileExt = allMat(num).ext; + + matFile = load([fileFolder filesep fileName fileExt]); + matFile = matFile.saveStructTemp; + + % ____________________________________ + % FIND THE DATAPOINT CURRENTLY LOADED + + D =matFile.common.dispersion; + a1 =matFile.common.a1; + chirp =matFile.common.chirp_alpha; + rop =matFile.common.rec_opt_pow; + + ber = mean(matFile.prms_compare_state.BER); + + wh.addValueToStorage(ber,'ber',D,a1,chirp,rop); + + + waitbar(num/numel(allMat),f,'Loading your data'); +end + +close(f) + + + + + +% 4) Hey! the warehouse is here and (hopefully) filled with data :-) + +% Create a save dialog +defaultDir = 'C:\Users\Silas\Documents\MATLAB\move_it_model_collection\'; +defaultExt = '*.mat'; +[filename, pathname] = uiputfile(fullfile(defaultDir, defaultExt),'', 'wh.mat'); + +% Check if the user pressed Cancel +if isequal(filename, 0) || isequal(pathname, 0) + disp('Save operation canceled.'); +else + % Save the variable to the selected file + save(fullfile(pathname, filename), 'wh'); + disp(['Variable "wh" saved to: ', fullfile(pathname, filename)]); +end + + + + +function matFileStructArray = getAllFilesInFolder(folderPath,extension) +% Get a list of all files in the current folder +currentFolderFiles = dir(fullfile(folderPath, '*')); + +% Exclude '.' and '..' directories +currentFolderFiles = currentFolderFiles(~ismember({currentFolderFiles.name}, {'.', '..'})); + +% Initialize the structure array for .mat files +matFileStructArray = struct('path', {}, 'name', {}, 'ext', {}); + +% Loop over each file in the current folder +for i = 1:length(currentFolderFiles) + currentFile = currentFolderFiles(i); + + % Check if the current item is a file and has a .mat extension + if ~currentFile.isdir && endsWith(currentFile.name, extension, 'IgnoreCase', true) + % If it's a .mat file, add it to the structure array + [matFileStructArray(end + 1).path,matFileStructArray(end+1).name, matFileStructArray(end+1).ext] = fileparts(fullfile(folderPath, currentFile.name)); + elseif currentFile.isdir + % If it's a directory, recursively call the function + subfolderPath = fullfile(folderPath, currentFile.name); + subfolderMatFiles = getAllFilesInFolder(subfolderPath,extension); + + % Add .mat files from the subfolder to the structure array + matFileStructArray = [matFileStructArray, subfolderMatFiles]; + end +end +end + + diff --git a/Classes/Warehouse_class/plotjob_example.mat b/Classes/Warehouse_class/plotjob_example.mat new file mode 100644 index 0000000000000000000000000000000000000000..0ee6062af6109417337cd3ee03aa5eda9d2f7cbf GIT binary patch literal 672 zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2cQV4Jk_w+L}(NSU?zf5VRA z2Z|O6`MMsv^X5D|6tcMF&HLYP3ez6dGE`h=f3Ygvzmhxbn&HOqPqt4jPe19-S@rC8 z!S$a0ng)Yp$#$>gtCw#bnkg)_Gt_0y4u#5RcYNh8Tw&}ioRrY+_2^E4#ocE@hNt{}lVWDBaPq}|sFTXo&&z=5#m07`Sx!X_K$7TFB`|+gq<()lW zi$5r9^u3p*5{oG|F&9UucX<@+huWor_FwD`ouK$ zef`H%f3sH~zr7)@-%9TClo$MRQ300D^3P^?pX%5V6mPcoTwC1n7o6U5X8X1u`_6Fq z>aJ}ct2XRg*9R1=YPhFf`s>pzw)o@y6Lvl4TKG<5MeRIxX)7x`%Xtgt6nx1k-KjKd mQO1lb*8KIp=i48Co*Y$sPw(TcG@GFP+wa-F=hqdSD*^yp@G3U| literal 0 HcmV?d00001