diff --git a/README.md b/README.md
index 2e2e4850437d8d1b7ae2cc3a3c457c9822a76bd5..d3e4b7bbd5bc899229ddab688722b643be75bc7d 100644
--- a/README.md
+++ b/README.md
@@ -1,92 +1,34 @@
-# Cone Dictionary Learning
-
-
-
-## Getting started
-
-To make it easy for you to get started with GitLab, here's a list of recommended next steps.
-
-Already a pro? Just edit this README.md and make it your own. Want to make it easy? [Use the template at the bottom](#editing-this-readme)!
-
-## Add your files
-
-- [ ] [Create](https://docs.gitlab.com/ee/user/project/repository/web_editor.html#create-a-file) or [upload](https://docs.gitlab.com/ee/user/project/repository/web_editor.html#upload-a-file) files
-- [ ] [Add files using the command line](https://docs.gitlab.com/ee/gitlab-basics/add-file.html#add-a-file-using-the-command-line) or push an existing Git repository with the following command:
-
-```
-cd existing_repo
-git remote add origin https://gitlab.cs.pub.ro/asydil/cone-dictionary-learning.git
-git branch -M main
-git push -uf origin main
-```
-
-## Integrate with your tools
-
-- [ ] [Set up project integrations](https://gitlab.cs.pub.ro/asydil/cone-dictionary-learning/-/settings/integrations)
-
-## Collaborate with your team
-
-- [ ] [Invite team members and collaborators](https://docs.gitlab.com/ee/user/project/members/)
-- [ ] [Create a new merge request](https://docs.gitlab.com/ee/user/project/merge_requests/creating_merge_requests.html)
-- [ ] [Automatically close issues from merge requests](https://docs.gitlab.com/ee/user/project/issues/managing_issues.html#closing-issues-automatically)
-- [ ] [Enable merge request approvals](https://docs.gitlab.com/ee/user/project/merge_requests/approvals/)
-- [ ] [Automatically merge when pipeline succeeds](https://docs.gitlab.com/ee/user/project/merge_requests/merge_when_pipeline_succeeds.html)
-
-## Test and Deploy
-
-Use the built-in continuous integration in GitLab.
-
-- [ ] [Get started with GitLab CI/CD](https://docs.gitlab.com/ee/ci/quick_start/index.html)
-- [ ] [Analyze your code for known vulnerabilities with Static Application Security Testing(SAST)](https://docs.gitlab.com/ee/user/application_security/sast/)
-- [ ] [Deploy to Kubernetes, Amazon EC2, or Amazon ECS using Auto Deploy](https://docs.gitlab.com/ee/topics/autodevops/requirements.html)
-- [ ] [Use pull-based deployments for improved Kubernetes management](https://docs.gitlab.com/ee/user/clusters/agent/)
-- [ ] [Set up protected environments](https://docs.gitlab.com/ee/ci/environments/protected_environments.html)
-
-***
-
-# Editing this README
-
-When you're ready to make this README your own, just edit this file and use the handy template below (or feel free to structure it however you want - this is just a starting point!). Thank you to [makeareadme.com](https://www.makeareadme.com/) for this template.
-
-## Suggestions for a good README
-Every project is different, so consider which of these sections apply to yours. The sections used in the template are suggestions for most open source projects. Also keep in mind that while a README can be too long and detailed, too long is better than too short. If you think your README is too long, consider utilizing another form of documentation rather than cutting out information.
-
-## Name
-Choose a self-explaining name for your project.
-
-## Description
-Let people know what your project can do specifically. Provide context and add a link to any reference visitors might be unfamiliar with. A list of Features or a Background subsection can also be added here. If there are alternatives to your project, this is a good place to list differentiating factors.
-
-## Badges
-On some READMEs, you may see small images that convey metadata, such as whether or not all the tests are passing for the project. You can use Shields to add some to your README. Many services also have instructions for adding a badge.
-
-## Visuals
-Depending on what you are making, it can be a good idea to include screenshots or even a video (you'll frequently see GIFs rather than actual videos). Tools like ttygif can help, but check out Asciinema for a more sophisticated method.
-
-## Installation
-Within a particular ecosystem, there may be a common way of installing things, such as using Yarn, NuGet, or Homebrew. However, consider the possibility that whoever is reading your README is a novice and would like more guidance. Listing specific steps helps remove ambiguity and gets people to using your project as quickly as possible. If it only runs in a specific context like a particular programming language version or operating system or has dependencies that have to be installed manually, also add a Requirements subsection.
-
-## Usage
-Use examples liberally, and show the expected output if you can. It's helpful to have inline the smallest example of usage that you can demonstrate, while providing links to more sophisticated examples if they are too long to reasonably include in the README.
-
-## Support
-Tell people where they can go to for help. It can be any combination of an issue tracker, a chat room, an email address, etc.
-
-## Roadmap
-If you have ideas for releases in the future, it is a good idea to list them in the README.
-
-## Contributing
-State if you are open to contributions and what your requirements are for accepting them.
-
-For people who want to make changes to your project, it's helpful to have some documentation on how to get started. Perhaps there is a script that they should run or some environment variables that they need to set. Make these steps explicit. These instructions could also be useful to your future self.
-
-You can also document commands to lint the code or run tests. These steps help to ensure high code quality and reduce the likelihood that the changes inadvertently break something. Having instructions for running tests is especially helpful if it requires external setup, such as starting a Selenium server for testing in a browser.
-
-## Authors and acknowledgment
-Show your appreciation to those who have contributed to the project.
-
-## License
-For open source projects, say how it is licensed.
-
-## Project status
-If you have run out of energy or time for your project, put a note at the top of the README saying that development has slowed down or stopped completely. Someone may choose to fork your project or volunteer to step in as a maintainer or owner, allowing your project to keep going. You can also make an explicit request for maintainers.
+# Dictionary learning with cone atoms and application to anomaly detection
+
+These programs implement the algorithms from the paper
+"Dictionary learning with cone atoms and application
+to anomaly detection", by Andra BÄltoiu, Denis C. Ilie-Ablachim,
+and Bogdan Dumitrescu, August 2023.
+
+The programs use data that should be downloaded from
+[http://asydil.upb.ro/wp-content/uploads/2023/08/dependency_30.zip](http://asydil.upb.ro/wp-content/uploads/2023/08/dependency_30.zip)
+and installed in a folder '../DATA_DEPENDENCY/'
+(or in any other folder that you like, but then you need to modify
+the test programs).
+
+The programs use a dictionary learning function called `DL`.
+Such a function can be found in the toolbox
+[https://gitlab.com/pirofti/dl-box](https://gitlab.com/pirofti/dl-box).
+The `DL` function is used only for comparison with standard dictionary learning algorithms.
+
+**Run first:**
+- `draw_err_conedl.m`: generates a plot like those in Figure 8 from the paper.
+ Change the variable `nf` to select another dataset.
+ The first figure (musk) was already generated and it will be loaded instantly.
+ For another figure, you'll have to wait a minute or so.
+
+- `test_all_dl_ad`: gives anomaly detection results (ROC AUC) like in the paper.
+ There are only two datasets, and the result is averaged on only `n_tests=2` tests,
+ for only 2 combinations of parameters, in order to have a short execution time.
+ The radii are uniformly distributed between `rad_min` and `rad_max`.
+
+**Important algorithms:**
+- `omp_cone_single.m`: sparse representation of a signal using Cone-OMP.
+- `dl_cone_paksvd_opt.m`: Cone-DL.
+
+Send any queries to [bogdan.dumitrescu@upb.ro](mailto:bogdan.dumitrescu@upb.ro).
diff --git a/adjust_radii.m b/adjust_radii.m
new file mode 100644
index 0000000000000000000000000000000000000000..2cfca2340664cd5ac51358e68ba880813b10e949
--- /dev/null
+++ b/adjust_radii.m
@@ -0,0 +1,54 @@
+function rad = adjust_radii(D, rad, min_dist)
+% Adjust radii such that the minimum distance between atoms is satisfied.
+%
+% Input:
+% D - dictionary (normalized atoms)
+% rad - cone radii (if scalar, then all atoms have the same radius)
+% min_dist - minimum distance between cones
+%
+% Output:
+% rad - adjusted radii
+
+% BD 22.05.2023
+
+tol = 1e-12;
+s_list = check_cone_superposition(D, rad, min_dist);
+
+i = size(s_list,1);
+%options = optimset('TolX',1e-4);
+while ~isempty(s_list) && i>0 % check list one by one
+ i1 = s_list(1,1);
+ i2 = s_list(1,2); % indices of atoms that are too close
+ r = sqrt(2-2*abs(D(:,i1)'*D(:,i2))) - min_dist+tol; % desired distance between the two atoms
+ t = rad(i2)/rad(i1);
+% rad(i1) = fzero(@(x) x*sqrt(1-x*x*t*t/4) + x*t*sqrt(1-x*x/4) - r, rad(i1)*r/(rad(i1)+rad(i2)), options);
+ rad(i1) = my_bisection(t, r, 0, rad(i1));
+ rad(i2) = rad(i1)*t;
+ s_list = check_cone_superposition(D, rad, min_dist);
+ i=i-1;
+end
+
+if ~isempty(s_list)
+ s_list
+end
+
+end
+
+function a = my_bisection(t, r, a, b)
+ f = @(x) x*sqrt(1-x*x*t*t/4) + x*t*sqrt(1-x*x/4) - r;
+ tol = 1e-6;
+ if 2/t < b % adjust upper bound to keep the function real
+ b = 2/t;
+ end
+ c = (a+b)/2;
+ fc = f(c);
+ while abs(fc) > tol
+ c = (a+b)/2;
+ fc = f(c);
+ if f(c)<0
+ a=c;
+ else
+ b=c;
+ end
+ end
+end
diff --git a/check_cone_superposition.m b/check_cone_superposition.m
new file mode 100644
index 0000000000000000000000000000000000000000..ba2b0884c7e3843beb8eaa35abf8a3d44c27b044
--- /dev/null
+++ b/check_cone_superposition.m
@@ -0,0 +1,43 @@
+function s_list = check_cone_superposition(D, rad, min_dist)
+% Check whether cone atoms have intersections
+%
+% Input:
+% D - dictionary (normalized atoms)
+% rad - cone radii (if scalar, then all atoms have the same radius)
+% min_dist - minimum distance between cones (default = 0)
+%
+% Output:
+% s_list - matrix with 2 columns; on each row, the indices of superposing atoms
+
+% BD 31.07.2022
+
+if nargin < 3
+ min_dist = 0;
+end
+
+tol = 1e-12;
+
+[~,n] = size(D);
+if length (rad) == 1
+ rad = rad * ones(n,1);
+end
+if length (rad) ~= n
+ error('Radii vector has wrong length')
+end
+rad = rad(:);
+
+G = eye(n);
+for i = 1 : n
+ for j = 1 : i
+% if D(:,i)'*D(:,j) < 0
+% G(i,j) = norm(D(:,i) + D(:,j));
+% else
+% G(i,j) = norm(D(:,i) - D(:,j));
+% end
+ G(i,j) = sqrt(2 - 2*abs(D(:,i)'*D(:,j)));
+ end
+end
+R = rad * sqrt(1 - rad'.*rad'/4) + sqrt(1 - rad.*rad/4) * rad' + min_dist;
+
+[rl, cl] = find(tril(G,-1) + tol < tril(R,-1));
+s_list = [rl, cl];
\ No newline at end of file
diff --git a/compute_rep_cone.m b/compute_rep_cone.m
new file mode 100644
index 0000000000000000000000000000000000000000..fc6e2acec016577f76c3525504003e1916167b23
--- /dev/null
+++ b/compute_rep_cone.m
@@ -0,0 +1,13 @@
+function [err, X] = compute_rep_cone(D, Y, s, rad, iter_cd, err_thresh)
+% Compute representations of many signals by calling omp_cone_single on
+% each of them. The main purpose is to get the representation error.
+
+[m,N] = size(Y);
+err = zeros(N,1);
+X = zeros(size(D,2),N);
+for ii = 1 : N
+ y = Y(:,ii);
+ [x, Da, supp] = omp_cone_single(y, D, s, rad, iter_cd, err_thresh);
+ err(ii) = sqrt(norm(y-Da*x)^2)/sqrt(m);
+ X(supp,ii) = x;
+end
diff --git a/data/dependency_Hepatitis.mat b/data/dependency_Hepatitis.mat
new file mode 100644
index 0000000000000000000000000000000000000000..4193ea5fa7db9e01eecdebf637d8434bba43f2ea
Binary files /dev/null and b/data/dependency_Hepatitis.mat differ
diff --git a/data/dependency_Ionosphere.mat b/data/dependency_Ionosphere.mat
new file mode 100644
index 0000000000000000000000000000000000000000..7fc50987297786369a1a6d93b7e201e199e689ba
Binary files /dev/null and b/data/dependency_Ionosphere.mat differ
diff --git a/data/dependency_Lymphography.mat b/data/dependency_Lymphography.mat
new file mode 100644
index 0000000000000000000000000000000000000000..e759540c3b58600d3310584c7003615f882a156c
Binary files /dev/null and b/data/dependency_Lymphography.mat differ
diff --git a/data/dependency_Pima.mat b/data/dependency_Pima.mat
new file mode 100644
index 0000000000000000000000000000000000000000..1027c5b3132dd97dd64dfbb21d08a87ecf4daa4d
Binary files /dev/null and b/data/dependency_Pima.mat differ
diff --git a/data/dependency_SpamBase.mat b/data/dependency_SpamBase.mat
new file mode 100644
index 0000000000000000000000000000000000000000..90d4965e64068ecf460ef06ae0dbf6b4c7b92fa1
Binary files /dev/null and b/data/dependency_SpamBase.mat differ
diff --git a/data/dependency_Stamps.mat b/data/dependency_Stamps.mat
new file mode 100644
index 0000000000000000000000000000000000000000..19077f3098aa3c26ba36f63bd072743a55332f9e
Binary files /dev/null and b/data/dependency_Stamps.mat differ
diff --git a/data/dependency_WBC.mat b/data/dependency_WBC.mat
new file mode 100644
index 0000000000000000000000000000000000000000..daa16674433dab3863d37588bd1dc4cfc579d036
Binary files /dev/null and b/data/dependency_WBC.mat differ
diff --git a/data/dependency_WDBC.mat b/data/dependency_WDBC.mat
new file mode 100644
index 0000000000000000000000000000000000000000..4ac7010014ccecafc641d7d40a2c65cd078030e6
Binary files /dev/null and b/data/dependency_WDBC.mat differ
diff --git a/data/dependency_WPBC.mat b/data/dependency_WPBC.mat
new file mode 100644
index 0000000000000000000000000000000000000000..d053a2ed0d6fa8e079c638b492f016635da93073
Binary files /dev/null and b/data/dependency_WPBC.mat differ
diff --git a/data/dependency_Waveform.mat b/data/dependency_Waveform.mat
new file mode 100644
index 0000000000000000000000000000000000000000..5d6a8300154697ef35ae86846050b8ca56b3b898
Binary files /dev/null and b/data/dependency_Waveform.mat differ
diff --git a/data/dependency_Wilt.mat b/data/dependency_Wilt.mat
new file mode 100644
index 0000000000000000000000000000000000000000..3cdfdd079cd8eafb5de8b2c6db9967f119904042
Binary files /dev/null and b/data/dependency_Wilt.mat differ
diff --git a/data/dependency_annthyroid.mat b/data/dependency_annthyroid.mat
new file mode 100644
index 0000000000000000000000000000000000000000..5c56318aa5aad7e647ce25b16b886a3b47d2204c
Binary files /dev/null and b/data/dependency_annthyroid.mat differ
diff --git a/data/dependency_breastw.mat b/data/dependency_breastw.mat
new file mode 100644
index 0000000000000000000000000000000000000000..89e07292d17cdfd235ec624398b36d0342c793d5
Binary files /dev/null and b/data/dependency_breastw.mat differ
diff --git a/data/dependency_cardio.mat b/data/dependency_cardio.mat
new file mode 100644
index 0000000000000000000000000000000000000000..1d1aec7742eacce1cf81c0bde3dcc929a48bde8f
Binary files /dev/null and b/data/dependency_cardio.mat differ
diff --git a/data/dependency_cover.mat b/data/dependency_cover.mat
new file mode 100644
index 0000000000000000000000000000000000000000..49fc65c24e37a78a5b14ccae2aede49a54af97b0
Binary files /dev/null and b/data/dependency_cover.mat differ
diff --git a/data/dependency_fault.mat b/data/dependency_fault.mat
new file mode 100644
index 0000000000000000000000000000000000000000..f938d22d27f44b43af65dc7389013a18a012706b
Binary files /dev/null and b/data/dependency_fault.mat differ
diff --git a/data/dependency_glass.mat b/data/dependency_glass.mat
new file mode 100644
index 0000000000000000000000000000000000000000..88b11f3138716b42c3011a751e3d7eb787459e66
Binary files /dev/null and b/data/dependency_glass.mat differ
diff --git a/data/dependency_landsat.mat b/data/dependency_landsat.mat
new file mode 100644
index 0000000000000000000000000000000000000000..83adc97004f9eff570f5bb8cceeb79019de5a17f
Binary files /dev/null and b/data/dependency_landsat.mat differ
diff --git a/data/dependency_letter.mat b/data/dependency_letter.mat
new file mode 100644
index 0000000000000000000000000000000000000000..0d64aa7c6608290dd3beae06e5c8bbb689c05302
Binary files /dev/null and b/data/dependency_letter.mat differ
diff --git a/data/dependency_magic.gamma.mat b/data/dependency_magic.gamma.mat
new file mode 100644
index 0000000000000000000000000000000000000000..af2a40648371a5fcb61beb2646e09a65d393a5de
Binary files /dev/null and b/data/dependency_magic.gamma.mat differ
diff --git a/data/dependency_mammography.mat b/data/dependency_mammography.mat
new file mode 100644
index 0000000000000000000000000000000000000000..7a0977857d24893610dd07df8fcda727141660e7
Binary files /dev/null and b/data/dependency_mammography.mat differ
diff --git a/data/dependency_musk.mat b/data/dependency_musk.mat
new file mode 100644
index 0000000000000000000000000000000000000000..f3e907ec48c944de920c4f10f1197a0203e05aec
Binary files /dev/null and b/data/dependency_musk.mat differ
diff --git a/data/dependency_pendigits.mat b/data/dependency_pendigits.mat
new file mode 100644
index 0000000000000000000000000000000000000000..a34d85d5c08dc29d80e3d3f876413d01322bcb84
Binary files /dev/null and b/data/dependency_pendigits.mat differ
diff --git a/data/dependency_satellite.mat b/data/dependency_satellite.mat
new file mode 100644
index 0000000000000000000000000000000000000000..fab7b5990b7be217624373a85dd6e106765f8300
Binary files /dev/null and b/data/dependency_satellite.mat differ
diff --git a/data/dependency_satimage.mat b/data/dependency_satimage.mat
new file mode 100644
index 0000000000000000000000000000000000000000..c37df27405d7dedd21dc27a68b1a0f18d316ed99
Binary files /dev/null and b/data/dependency_satimage.mat differ
diff --git a/data/dependency_thyroid.mat b/data/dependency_thyroid.mat
new file mode 100644
index 0000000000000000000000000000000000000000..09ede8a7aa2742578af40222d54d14c70ae8e6bf
Binary files /dev/null and b/data/dependency_thyroid.mat differ
diff --git a/data/dependency_vertebral.mat b/data/dependency_vertebral.mat
new file mode 100644
index 0000000000000000000000000000000000000000..0006540b26749d2a0fcd18a725be5ef88df10fb5
Binary files /dev/null and b/data/dependency_vertebral.mat differ
diff --git a/data/dependency_vowels.mat b/data/dependency_vowels.mat
new file mode 100644
index 0000000000000000000000000000000000000000..a44b4ae7d2a6ce6f24adbe06ee972afd58907c8d
Binary files /dev/null and b/data/dependency_vowels.mat differ
diff --git a/data/dependency_wine.mat b/data/dependency_wine.mat
new file mode 100644
index 0000000000000000000000000000000000000000..e84f7b8f35c7695d7596fd489dc8778ce9263ed3
Binary files /dev/null and b/data/dependency_wine.mat differ
diff --git a/data/dependency_yeast.mat b/data/dependency_yeast.mat
new file mode 100644
index 0000000000000000000000000000000000000000..cf9d602720331cfd56a193774a7f5314d85fb9ef
Binary files /dev/null and b/data/dependency_yeast.mat differ
diff --git a/disjoin_cones.m b/disjoin_cones.m
new file mode 100644
index 0000000000000000000000000000000000000000..78279cbbc87905856df0a23619e01ea0c13a9989
--- /dev/null
+++ b/disjoin_cones.m
@@ -0,0 +1,128 @@
+function Dup = disjoin_cones(D, Dup, rad, min_dist)
+% Ensure that cones of new dictionary are disjoint, by moving them towards
+% the old dictionary when necessary.
+% If some cones of the old dictionary are too close, the corresponding new
+% cones keep their values or take the old ones, if the distance between them
+% is larger or smaller, respectively, than the old distance.
+% Input:
+% D - old dictionary (central atoms); it is assumed that the cones
+% generated by D are disjoint
+% Dup - updated dictionary (central atoms)
+% rad - cone radii (if scalar, then all atoms have the same radius)
+% min_dist - minimum distance between cones ("zero" = the sum of radii)
+%
+% Output:
+% Dup - updated dictionary with disjoint atoms
+
+% BD 31.07.2022
+
+if nargin < 4
+ min_dist = 0;
+end
+
+tol = 1e-12;
+
+[~,n] = size(D);
+
+if length (rad) == 1
+ rad = rad * ones(n,1);
+end
+
+% align atoms (in case they are opposed)
+for i = 1:n
+ if D(:,i)'*Dup(:,i) < 0
+ Dup(:,i) = -Dup(:,i);
+ end
+end
+
+% check if there are cones too close
+s_list = check_cone_superposition(Dup, rad, min_dist);
+it_nr = 0;
+old_too_close_list = []; % list of pairs in the old dictionary that are too close
+while ~isempty(s_list)
+ for k = 1 : size(s_list,1)
+ i1 = s_list(k,1);
+ i2 = s_list(k,2); % indices of atoms that are too close
+ if D(:,i1)'*D(:,i2) < 0
+ D(:,i2) = - D(:,i2);
+ end
+ if Dup(:,i1)'*Dup(:,i2) < 0
+ D(:,i2) = - D(:,i2);
+ end
+
+ %d = norm(D(:,i1) - D(:,i2)); % distance between old atoms
+ %dup = norm(Dup(:,i1) - Dup(:,i2)); % distance between new atoms
+ d = sqrt(2 - 2*abs(D(:,i1)'*D(:,i2))); % distance between old atoms
+ dup = sqrt(2 - 2*abs(Dup(:,i1)'*Dup(:,i2))); % distance between new atoms
+ dok = rad(i1)*sqrt(1 - rad(i2)*rad(i2)/4) + rad(i2)*sqrt(1 - rad(i1)*rad(i1)/4) + min_dist; % desired distance
+ if d + tol < dok % check whether D respected the distance rule
+ if it_nr == 0 % action is needed only at the first iteration
+ if dup < d % the distance has decreased, keep the old dictionary
+ Dup(:,i1) = D(:,i1);
+ Dup(:,i2) = D(:,i2);
+ end % otherwise, take the new atoms as they are
+ % update the list of atoms that are too close and cannot be moved
+ old_too_close_list = [old_too_close_list; s_list(k,:)];
+ end
+ else % move towards old atoms as much as necessary
+ if 0 % successive approximation
+ for ii = 1 : 5 % approximate, using proportionality on straight lines (not on sphere)
+ f = (d - dok) / (d - dup); % fraction of distance with respect to old atom
+ Dup(:,i1) = normc((1-f)*D(:,i1) + f*Dup(:,i1));
+ Dup(:,i2) = normc((1-f)*D(:,i2) + f*Dup(:,i2));
+ dup = norm(Dup(:,i1) - Dup(:,i2));
+ end
+ else
+ % bisection algorithm
+ fi = 0;
+ ff = 1;
+ d1i = D(:,i1);
+ d1f = Dup(:,i1);
+ d2i = D(:,i2);
+ d2f = Dup(:,i2);
+ for ii = 1:12
+ f = (fi+ff)/2;
+ d1 = normc((1-f)*d1i + f*d1f);
+ d2 = normc((1-f)*d2i + f*d2f);
+ % if norm(d1-d2) < dok
+ if sqrt(2-2*abs(d1'*d2)) < dok
+ ff=f;
+ d1f = d1;
+ d2f = d2;
+ else
+ fi = f;
+ d1i = d1;
+ d2i = d2;
+ end
+ end
+ Dup(:,i1) = d1i;
+ Dup(:,i2) = d2i;
+ end
+ end
+ end
+ s_list = check_cone_superposition(Dup, rad, min_dist); % check again after solving, maybe new conflicts appeared
+
+ if it_nr > 10
+ s_list
+ end
+
+ if isempty(s_list) % cut it short if dictionary is ok
+ break
+ end
+
+ % check whether the list contains only old atoms that are too close
+ if ~isempty(old_too_close_list)
+ only_old = 1;
+ for k = 1 : size(s_list,1)
+ if ~ismember(s_list(k,:), old_too_close_list, 'rows') % not old
+ only_old = 0;
+ break
+ end
+ end
+ if only_old % no more solvable conflicts
+ break
+ end
+ end
+
+ it_nr = it_nr+1;
+end
diff --git a/dl_cone_paksvd_opt.m b/dl_cone_paksvd_opt.m
new file mode 100644
index 0000000000000000000000000000000000000000..fc5042a9c7c012cae455771cf76c4072df8e100c
--- /dev/null
+++ b/dl_cone_paksvd_opt.m
@@ -0,0 +1,86 @@
+function [D, X, err] = dl_cone_paksvd_opt(Y, D, s, iternum, rad, iter_cd, err_thresh, ...
+ min_d_cones, update_quota, q_dec_rate)
+% Cone dictionary learning. Dictionary update is made in Parallel AK-SVD style,
+% using an optimal update of the atoms.
+% The radii of the cones are not changed.
+%
+% Input:
+% Y - set of signals
+% D - dictionary (normalized atoms)
+% s - sparsity level
+% iternum - number of learning iterations
+% rad - cone radii (if scalar, then all atoms have the same radius)
+% iter_cd - number of iterations in the coordinate descent refinement
+% (default = 1)
+% err_thresh - representation error for which OMP is stopped
+% min_d_cones - minimum distance between cones
+% update_quota- fraction of atoms that are updated in one pass over data (default = 1)
+% q_dec_rate - quota decrease rate (update_quota *= q_dec_rate in each iteration)
+%
+% Output:
+% D - trained dictionary
+% X - representation matrix
+% err - vector of errors at each iteration
+
+% BD 24.07.2022
+
+if nargin < 8
+ min_d_cones = 0;
+end
+if nargin < 9
+ update_quota = 1;
+end
+if nargin < 10
+ q_dec_rate = 1;
+end
+
+[m,n] = size(D);
+N = size(Y,2);
+
+err = zeros(1,iternum);
+for k = 1 : iternum
+ n_update = round(update_quota*n);
+ X = zeros(n,N);
+ Dav = zeros(m,n);
+ e = 0;
+ if n_update < n
+ crt_atoms = randperm(n);
+ crt_atoms = crt_atoms(1:n_update);
+ end
+ for i = 1 : N
+ % compute representation of a signal
+ y = Y(:,i);
+ [x, Da, support] = omp_cone_single(y, D, s, rad, iter_cd, err_thresh);
+ X(support,i) = x;
+ r = y - Da*x; % representation residual
+ if n_update == n % fully parallel case (don't waste time with column computation)
+ for j = 1 : length(support)
+ f = r + Da(:,j)*x(j);
+ Dav(:,support(j)) = Dav(:,support(j)) + x(j)*(f - x(j)*(Da(:,j) - D(:,support(j))));
+ end
+ else
+ for j = 1 : length(support)
+ if any(support(j)==crt_atoms)
+ f = r + Da(:,j)*x(j);
+ Dav(:,support(j)) = Dav(:,support(j)) + x(j)*(f - x(j)*(Da(:,j) - D(:,support(j))));
+ end
+ end
+ end
+ e = e + norm(r)^2;
+ end
+ e = sqrt(e)/sqrt(m*N);
+ err(k) = e;
+
+ % if new atoms are zero, replace with old atoms
+ for j = 1:n
+ if norm(Dav(:,j)) == 0
+ Dav(:,j) = D(:,j);
+ end
+ end
+
+ % new atoms
+ D = disjoin_cones(D, normc(Dav), rad, min_d_cones);
+ %min_dist_dict(D, rad)
+
+ update_quota = update_quota * q_dec_rate;
+end
diff --git a/draw_err_conedl.asv b/draw_err_conedl.asv
new file mode 100644
index 0000000000000000000000000000000000000000..1f3bf08daae47514ffab36abb5b75f056adc72f8
--- /dev/null
+++ b/draw_err_conedl.asv
@@ -0,0 +1,171 @@
+% comparison between Cone-DL and AK-SVD+Cone-OMP
+% for data from ADBench (but not only)
+
+% BD 26.06.2023
+
+datasets_sort_m = {'musk', 'SpamBase', 'landsat', 'satellite', 'satimage', ...
+ 'ionosphere', 'WPBC', 'letter', 'WDBC', 'fault', ...
+ 'cardio', 'Waveform', 'Hepatitis', 'Lymphography', 'pendigits', ...
+ 'wine', 'vowels', 'cover', 'magic.gamma.mat', 'breastw', ...
+ 'Stamps', 'WBC', 'yeast', 'glass', 'annthyroid', ...
+ 'mammography', 'Pima', 'thyroid', 'vertebral', 'Wilt', ...
+ 'skin', 'smtp' };
+
+nf = 1; % number of dataset in th
+data_file_name = ['../DATA_DEPENDENCY/dependency_', datasets_sort_m{nf}];
+
+n_tests = 10;
+
+% parameters and sizes
+overcompleteness = 3;
+s = 3;
+rad_max = 0.2; %0.25;
+rad_min = 0.01; %0.1;
+iternum = 200;
+iter_cd = 10; %10;
+err_thresh = 1e-7;
+min_dist_cones = 0.01; %0.1;
+init_dl = 0;
+iter_init = 10; % iterations for DL initialization
+train_perc = 0.8; % for Selective Cone-DL
+train_drop_perc = 0.2;
+
+alg_legend = {'Cone-DL', 'Cone-DL + swap', 'AK-SVD + Cone-OMP', 'AK-SVD + Cone-OMP + swap'};
+
+% check if data file already exists
+fig_data_name = ['./figs/fig_data_', datasets_sort_m{nf}, '_s_', num2str(s), ...
+ '_c_', num2str(overcompleteness), '_r_', num2str(rad_min), '_', num2str(rad_max), '.mat'];
+if isfile(fig_data_name)
+ load(fig_data_name)
+else % compute if data do not exist
+
+ % data loading and normalization
+ fprintf("%s\n", data_file_name)
+ load(data_file_name)
+ Y=Y';
+ [m,N] = size(Y);
+ vm = mean(Y,2); % Z-score normalization
+ Y = double(Y) - repmat(vm,1,N); % subtract mean
+ sm = std(Y,0,2);
+ Y = Y ./ repmat(sm,1,N);
+
+ n = overcompleteness * m;
+
+ err1_av = 0;
+ err2_av = 0;
+ err3_av = 0;
+ err4_av = 0;
+ for it = 1 : n_tests
+ dist_ok = 0;
+ i0 = 0;
+ while ~dist_ok
+ % generate radii
+ rad = rad_min + (rad_max-rad_min) * rand(n,1);
+ D0 = randn(m,n);
+ D0 = normc(D0); % initial dictionary
+ if init_dl
+ D0 = DL(Y, D0, s, iter_init);
+ end
+ if isempty(check_cone_superposition(D0, rad))
+ dist_ok = 1;
+ end
+ i0 = i0+1;
+ if i0 == 10
+ break
+ end
+ end
+ if i0==10
+ error('D0 cones superposition')
+ end
+
+ % dictionary learning
+ % Cone-DL
+ [D1, ~, err1] = dl_cone_paksvd_opt(Y, D0, s, iternum, rad, iter_cd, err_thresh, min_dist_cones);
+ err1_av = err1_av + err1;
+ % swap radii by use
+ [~, Xs] = compute_rep_cone(D1, Y, s, rad, iter_cd, err_thresh);
+ atom_use = sum(squeeze(Xs)~=0,2);
+ [~,Ip] = sort(atom_use,'descend');
+ [~,Ir] = sort(rad,'descend');
+ rad2 = zeros(n,1);
+ rad2(Ip) = rad(Ir); % swap original radii!!!
+ rad2 = adjust_radii(D1, rad2, min_dist_cones);
+
+ err2 = 0;
+ for i = 1 : N
+ y = Y(:,i);
+ [x, Da, support] = omp_cone_single(y, D1, s, rad2, iter_cd, err_thresh);
+ err2 = err2 + norm(y - Da*x)^2;
+ end
+ err2 = sqrt(err2)/sqrt(m*N);
+ err2_av = err2_av + err2;
+
+ % Selective Cone-DL
+ %[D2, err2] = dl_cone_paksvd_opt_selective(Y, D0, s, iternum, rad, iter_cd, err_thresh, ...
+ % train_perc, train_drop_perc, min_dist_cones);
+
+ % AK-SVD
+ update = 'aksvd';
+ replatom = 'random';
+ D3 = DL(Y, D0, s, iternum, str2func(update), {}, 'replatom', replatom); % standard DL
+ % run Cone-OMP with initial radii
+ rad3 = adjust_radii(D3, rad, min_dist_cones); % ensure that cones are disjoint for AK-SVD
+
+ err3 = 0;
+ for i = 1 : N
+ y = Y(:,i);
+ [x, Da, support] = omp_cone_single(y, D3, s, rad3, iter_cd, err_thresh);
+ err3 = err3 + norm(y - Da*x)^2;
+ end
+ err3 = sqrt(err3)/sqrt(m*N);
+ err3_av = err3_av + err3;
+
+ % swap radii by use
+ [~, Xs] = compute_rep_cone(D3, Y, s, rad3, iter_cd, err_thresh);
+ atom_use = sum(squeeze(Xs)~=0,2);
+ [~,Ip] = sort(atom_use,'descend');
+ [~,Ir] = sort(rad3,'descend');
+ rad_swap = zeros(n,1);
+ rad_swap(Ip) = rad(Ir); % swap original radii!!!
+ rad_swap = adjust_radii(D3, rad_swap, min_dist_cones);
+
+ err4 = 0;
+ for i = 1 : N
+ y = Y(:,i);
+ [x, Da, support] = omp_cone_single(y, D3, s, rad_swap, iter_cd, err_thresh);
+ err4 = err4 + norm(y - Da*x)^2;
+ end
+ err4 = sqrt(err4)/sqrt(m*N);
+ err4_av = err4_av + err4;
+ end
+
+ best_err = [min(err1_av) err2_av err3_av err4_av] / n_tests;
+ % show sorted best errors
+ for i = 1 : length(alg_legend)
+ fprintf("%24s: %8.6f\n", alg_legend{i}, best_err(i));
+ end
+
+ % save data
+ save(fig_data_name, "err1_av", "err2_av", "err3_av", "err4_av", "n_tests")
+end
+
+% plot
+font_size = 12;
+line_width = 2;
+figure, hold on
+plot(log10(err1_av/n_tests), 'LineWidth', line_width)
+%plot(log10(err2))
+plot(log10(err2_av/n_tests*ones(1,iternum)), 'LineWidth', line_width)
+plot(log10(err3_av/n_tests*ones(1,iternum)), 'LineWidth', line_width)
+plot(log10(err4_av/n_tests*ones(1,iternum)), 'LineWidth', line_width)
+legend(alg_legend, 'FontSize', font_size)
+xlabel('iteration', 'FontSize', font_size)
+ylabel('log10(error)', 'FontSize', font_size)
+title(datasets_sort_m{nf}, 'FontSize', font_size)
+grid
+hold off
+
+fig_name = ['./figs/fig_', datasets_sort_m{nf}, '_s_', num2str(s), '_c_', num2str(overcompleteness), ...
+ '_r_', num2str(rad_min), '_', num2str(rad_max), '.eps'];
+print(fig_name, '-depsc')
+
diff --git a/draw_err_conedl.m b/draw_err_conedl.m
new file mode 100644
index 0000000000000000000000000000000000000000..b98226a5adad63430da88dbf2d44306eb4f30083
--- /dev/null
+++ b/draw_err_conedl.m
@@ -0,0 +1,163 @@
+% comparison between Cone-DL and AK-SVD+Cone-OMP
+% for data from ADBench (but not only)
+
+% BD 26.06.2023
+
+datasets_sort_m = {'musk', 'SpamBase', 'landsat', 'satellite', 'satimage', ...
+ 'ionosphere', 'WPBC', 'letter', 'WDBC', 'fault', ...
+ 'cardio', 'Waveform', 'Hepatitis', 'Lymphography', 'pendigits', ...
+ 'wine', 'vowels', 'cover', 'magic.gamma.mat', 'breastw', ...
+ 'Stamps', 'WBC', 'yeast', 'glass', 'annthyroid', ...
+ 'mammography', 'Pima', 'thyroid', 'vertebral', 'Wilt'};
+
+nf = 1; % number of dataset in the list above
+data_file_name = ['../DATA_DEPENDENCY/dependency_', datasets_sort_m{nf}];
+
+n_tests = 10; % number of tests on which the average is made
+
+% parameters and sizes
+overcompleteness = 3;
+s = 3;
+rad_max = 0.2;
+rad_min = 0.01;
+iternum = 200;
+iter_cd = 10;
+err_thresh = 1e-7;
+min_dist_cones = 0.01;
+init_dl = 0;
+iter_init = 10; % iterations for DL initialization
+
+alg_legend = {'Cone-DL', 'Cone-DL + swap', 'AK-SVD + Cone-OMP', 'AK-SVD + Cone-OMP + swap'};
+
+% check if data file already exists
+fig_data_name = ['./figs/fig_data_', datasets_sort_m{nf}, '_s_', num2str(s), ...
+ '_c_', num2str(overcompleteness), '_r_', num2str(rad_min), '_', num2str(rad_max), '.mat'];
+if isfile(fig_data_name)
+ load(fig_data_name)
+else % compute if data do not exist
+
+ % data loading and normalization
+ fprintf("%s\n", data_file_name)
+ load(data_file_name)
+ Y=Y';
+ [m,N] = size(Y);
+ vm = mean(Y,2); % Z-score normalization
+ Y = double(Y) - repmat(vm,1,N); % subtract mean
+ sm = std(Y,0,2);
+ Y = Y ./ repmat(sm,1,N);
+
+ n = overcompleteness * m;
+
+ err1_av = 0;
+ err2_av = 0;
+ err3_av = 0;
+ err4_av = 0;
+ for it = 1 : n_tests
+ dist_ok = 0;
+ i0 = 0;
+ while ~dist_ok
+ % generate radii
+ rad = rad_min + (rad_max-rad_min) * rand(n,1);
+ D0 = randn(m,n);
+ D0 = normc(D0); % initial dictionary
+ if init_dl
+ D0 = DL(Y, D0, s, iter_init);
+ end
+ if isempty(check_cone_superposition(D0, rad))
+ dist_ok = 1;
+ end
+ i0 = i0+1;
+ if i0 == 10
+ break
+ end
+ end
+ if i0==10
+ error('D0 cones superposition')
+ end
+
+ % dictionary learning
+ % Cone-DL
+ [D1, ~, err1] = dl_cone_paksvd_opt(Y, D0, s, iternum, rad, iter_cd, err_thresh, min_dist_cones);
+ err1_av = err1_av + err1;
+ % swap radii by use
+ [~, Xs] = compute_rep_cone(D1, Y, s, rad, iter_cd, err_thresh);
+ atom_use = sum(squeeze(Xs)~=0,2);
+ [~,Ip] = sort(atom_use,'descend');
+ [~,Ir] = sort(rad,'descend');
+ rad2 = zeros(n,1);
+ rad2(Ip) = rad(Ir); % swap original radii!!!
+ rad2 = adjust_radii(D1, rad2, min_dist_cones);
+
+ err2 = 0;
+ for i = 1 : N
+ y = Y(:,i);
+ [x, Da, support] = omp_cone_single(y, D1, s, rad2, iter_cd, err_thresh);
+ err2 = err2 + norm(y - Da*x)^2;
+ end
+ err2 = sqrt(err2)/sqrt(m*N);
+ err2_av = err2_av + err2;
+
+ % AK-SVD
+ update = 'aksvd';
+ replatom = 'random';
+ D3 = DL(Y, D0, s, iternum, str2func(update), {}, 'replatom', replatom); % standard DL
+ % run Cone-OMP with initial radii
+ rad3 = adjust_radii(D3, rad, min_dist_cones); % ensure that cones are disjoint for AK-SVD
+
+ err3 = 0;
+ for i = 1 : N
+ y = Y(:,i);
+ [x, Da, support] = omp_cone_single(y, D3, s, rad3, iter_cd, err_thresh);
+ err3 = err3 + norm(y - Da*x)^2;
+ end
+ err3 = sqrt(err3)/sqrt(m*N);
+ err3_av = err3_av + err3;
+
+ % swap radii by use
+ [~, Xs] = compute_rep_cone(D3, Y, s, rad3, iter_cd, err_thresh);
+ atom_use = sum(squeeze(Xs)~=0,2);
+ [~,Ip] = sort(atom_use,'descend');
+ [~,Ir] = sort(rad3,'descend');
+ rad_swap = zeros(n,1);
+ rad_swap(Ip) = rad(Ir); % swap original radii!!!
+ rad_swap = adjust_radii(D3, rad_swap, min_dist_cones);
+
+ err4 = 0;
+ for i = 1 : N
+ y = Y(:,i);
+ [x, Da, support] = omp_cone_single(y, D3, s, rad_swap, iter_cd, err_thresh);
+ err4 = err4 + norm(y - Da*x)^2;
+ end
+ err4 = sqrt(err4)/sqrt(m*N);
+ err4_av = err4_av + err4;
+ end
+
+ best_err = [min(err1_av) err2_av err3_av err4_av] / n_tests;
+ % show sorted best errors
+ for i = 1 : length(alg_legend)
+ fprintf("%24s: %8.6f\n", alg_legend{i}, best_err(i));
+ end
+
+ % save data
+ save(fig_data_name, "err1_av", "err2_av", "err3_av", "err4_av", "n_tests")
+end
+
+% plot
+font_size = 12;
+line_width = 2;
+figure, hold on
+plot(log10(err1_av/n_tests), 'LineWidth', line_width)
+plot(log10(err2_av/n_tests*ones(1,iternum)), 'LineWidth', line_width)
+plot(log10(err3_av/n_tests*ones(1,iternum)), 'LineWidth', line_width)
+plot(log10(err4_av/n_tests*ones(1,iternum)), 'LineWidth', line_width)
+legend(alg_legend, 'FontSize', font_size)
+xlabel('iteration', 'FontSize', font_size)
+ylabel('log10(error)', 'FontSize', font_size)
+title(datasets_sort_m{nf}, 'FontSize', font_size)
+grid
+hold off
+
+fig_name = ['./figs/fig_', datasets_sort_m{nf}, '_s_', num2str(s), '_c_', num2str(overcompleteness), ...
+ '_r_', num2str(rad_min), '_', num2str(rad_max), '.eps'];
+print(fig_name, '-depsc')
+
diff --git a/figs/fig_data_musk_s_3_c_3_r_0.01_0.2.mat b/figs/fig_data_musk_s_3_c_3_r_0.01_0.2.mat
new file mode 100644
index 0000000000000000000000000000000000000000..6ecfe5f6cc6e2deb6e90c8e3487a29e798b3a5b9
Binary files /dev/null and b/figs/fig_data_musk_s_3_c_3_r_0.01_0.2.mat differ
diff --git a/figs/fig_musk_s_3_c_3_r_0.01_0.2.eps b/figs/fig_musk_s_3_c_3_r_0.01_0.2.eps
new file mode 100644
index 0000000000000000000000000000000000000000..5ba6d45ed39da2b5d2cb94ce39459152c25ad857
--- /dev/null
+++ b/figs/fig_musk_s_3_c_3_r_0.01_0.2.eps
@@ -0,0 +1,1995 @@
+%!PS-Adobe-3.0 EPSF-3.0
+%%Creator: (MATLAB, The Mathworks, Inc. Version 9.14.0.2254940 \(R2023a\) Update 2. Operating System: Windows 10)
+%%Title: ./figs/fig_musk_s_3_c_3_r_0.01_0.2.eps
+%%CreationDate: 2023-08-12T11:34:42
+%%Pages: (atend)
+%%BoundingBox: 18 4 388 304
+%%LanguageLevel: 3
+%%EndComments
+%%BeginProlog
+%%BeginResource: procset (Apache XML Graphics Std ProcSet) 1.2 0
+%%Version: 1.2 0
+%%Copyright: (Copyright 2001-2003,2010 The Apache Software Foundation. License terms: http://www.apache.org/licenses/LICENSE-2.0)
+/bd{bind def}bind def
+/ld{load def}bd
+/GR/grestore ld
+/GS/gsave ld
+/RM/rmoveto ld
+/C/curveto ld
+/t/show ld
+/L/lineto ld
+/ML/setmiterlimit ld
+/CT/concat ld
+/f/fill ld
+/N/newpath ld
+/S/stroke ld
+/CC/setcmykcolor ld
+/A/ashow ld
+/cp/closepath ld
+/RC/setrgbcolor ld
+/LJ/setlinejoin ld
+/GC/setgray ld
+/LW/setlinewidth ld
+/M/moveto ld
+/re {4 2 roll M
+1 index 0 rlineto
+0 exch rlineto
+neg 0 rlineto
+cp } bd
+/_ctm matrix def
+/_tm matrix def
+/BT { _ctm currentmatrix pop matrix _tm copy pop 0 0 moveto } bd
+/ET { _ctm setmatrix } bd
+/iTm { _ctm setmatrix _tm concat } bd
+/Tm { _tm astore pop iTm 0 0 moveto } bd
+/ux 0.0 def
+/uy 0.0 def
+/F {
+ /Tp exch def
+ /Tf exch def
+ Tf findfont Tp scalefont setfont
+ /cf Tf def /cs Tp def
+} bd
+/ULS {currentpoint /uy exch def /ux exch def} bd
+/ULE {
+ /Tcx currentpoint pop def
+ gsave
+ newpath
+ cf findfont cs scalefont dup
+ /FontMatrix get 0 get /Ts exch def /FontInfo get dup
+ /UnderlinePosition get Ts mul /To exch def
+ /UnderlineThickness get Ts mul /Tt exch def
+ ux uy To add moveto Tcx uy To add lineto
+ Tt setlinewidth stroke
+ grestore
+} bd
+/OLE {
+ /Tcx currentpoint pop def
+ gsave
+ newpath
+ cf findfont cs scalefont dup
+ /FontMatrix get 0 get /Ts exch def /FontInfo get dup
+ /UnderlinePosition get Ts mul /To exch def
+ /UnderlineThickness get Ts mul /Tt exch def
+ ux uy To add cs add moveto Tcx uy To add cs add lineto
+ Tt setlinewidth stroke
+ grestore
+} bd
+/SOE {
+ /Tcx currentpoint pop def
+ gsave
+ newpath
+ cf findfont cs scalefont dup
+ /FontMatrix get 0 get /Ts exch def /FontInfo get dup
+ /UnderlinePosition get Ts mul /To exch def
+ /UnderlineThickness get Ts mul /Tt exch def
+ ux uy To add cs 10 mul 26 idiv add moveto Tcx uy To add cs 10 mul 26 idiv add lineto
+ Tt setlinewidth stroke
+ grestore
+} bd
+/QT {
+/Y22 exch store
+/X22 exch store
+/Y21 exch store
+/X21 exch store
+currentpoint
+/Y21 load 2 mul add 3 div exch
+/X21 load 2 mul add 3 div exch
+/X21 load 2 mul /X22 load add 3 div
+/Y21 load 2 mul /Y22 load add 3 div
+/X22 load /Y22 load curveto
+} bd
+/SSPD {
+dup length /d exch dict def
+{
+/v exch def
+/k exch def
+currentpagedevice k known {
+/cpdv currentpagedevice k get def
+v cpdv ne {
+/upd false def
+/nullv v type /nulltype eq def
+/nullcpdv cpdv type /nulltype eq def
+nullv nullcpdv or
+{
+/upd true def
+} {
+/sametype v type cpdv type eq def
+sametype {
+v type /arraytype eq {
+/vlen v length def
+/cpdvlen cpdv length def
+vlen cpdvlen eq {
+0 1 vlen 1 sub {
+/i exch def
+/obj v i get def
+/cpdobj cpdv i get def
+obj cpdobj ne {
+/upd true def
+exit
+} if
+} for
+} {
+/upd true def
+} ifelse
+} {
+v type /dicttype eq {
+v {
+/dv exch def
+/dk exch def
+/cpddv cpdv dk get def
+dv cpddv ne {
+/upd true def
+exit
+} if
+} forall
+} {
+/upd true def
+} ifelse
+} ifelse
+} if
+} ifelse
+upd true eq {
+d k v put
+} if
+} if
+} if
+} forall
+d length 0 gt {
+d setpagedevice
+} if
+} bd
+/RE { % /NewFontName [NewEncodingArray] /FontName RE -
+ findfont dup length dict begin
+ {
+ 1 index /FID ne
+ {def} {pop pop} ifelse
+ } forall
+ /Encoding exch def
+ /FontName 1 index def
+ currentdict definefont pop
+ end
+} bind def
+%%EndResource
+%%BeginResource: procset (Apache XML Graphics EPS ProcSet) 1.0 0
+%%Version: 1.0 0
+%%Copyright: (Copyright 2002-2003 The Apache Software Foundation. License terms: http://www.apache.org/licenses/LICENSE-2.0)
+/BeginEPSF { %def
+/b4_Inc_state save def % Save state for cleanup
+/dict_count countdictstack def % Count objects on dict stack
+/op_count count 1 sub def % Count objects on operand stack
+userdict begin % Push userdict on dict stack
+/showpage { } def % Redefine showpage, { } = null proc
+0 setgray 0 setlinecap % Prepare graphics state
+1 setlinewidth 0 setlinejoin
+10 setmiterlimit [ ] 0 setdash newpath
+/languagelevel where % If level not equal to 1 then
+{pop languagelevel % set strokeadjust and
+1 ne % overprint to their defaults.
+{false setstrokeadjust false setoverprint
+} if
+} if
+} bd
+/EndEPSF { %def
+count op_count sub {pop} repeat % Clean up stacks
+countdictstack dict_count sub {end} repeat
+b4_Inc_state restore
+} bd
+%%EndResource
+%FOPBeginFontDict
+%%IncludeResource: font Courier-Oblique
+%%IncludeResource: font Courier-BoldOblique
+%%IncludeResource: font Courier-Bold
+%%IncludeResource: font ZapfDingbats
+%%IncludeResource: font Symbol
+%%IncludeResource: font Helvetica
+%%IncludeResource: font Helvetica-Oblique
+%%IncludeResource: font Helvetica-Bold
+%%IncludeResource: font Helvetica-BoldOblique
+%%IncludeResource: font Times-Roman
+%%IncludeResource: font Times-Italic
+%%IncludeResource: font Times-Bold
+%%IncludeResource: font Times-BoldItalic
+%%IncludeResource: font Courier
+%FOPEndFontDict
+%%BeginResource: encoding WinAnsiEncoding
+/WinAnsiEncoding [
+/.notdef /.notdef /.notdef /.notdef /.notdef
+/.notdef /.notdef /.notdef /.notdef /.notdef
+/.notdef /.notdef /.notdef /.notdef /.notdef
+/.notdef /.notdef /.notdef /.notdef /.notdef
+/.notdef /.notdef /.notdef /.notdef /.notdef
+/.notdef /.notdef /.notdef /.notdef /.notdef
+/.notdef /.notdef /space /exclam /quotedbl
+/numbersign /dollar /percent /ampersand /quotesingle
+/parenleft /parenright /asterisk /plus /comma
+/hyphen /period /slash /zero /one
+/two /three /four /five /six
+/seven /eight /nine /colon /semicolon
+/less /equal /greater /question /at
+/A /B /C /D /E
+/F /G /H /I /J
+/K /L /M /N /O
+/P /Q /R /S /T
+/U /V /W /X /Y
+/Z /bracketleft /backslash /bracketright /asciicircum
+/underscore /quoteleft /a /b /c
+/d /e /f /g /h
+/i /j /k /l /m
+/n /o /p /q /r
+/s /t /u /v /w
+/x /y /z /braceleft /bar
+/braceright /asciitilde /bullet /Euro /bullet
+/quotesinglbase /florin /quotedblbase /ellipsis /dagger
+/daggerdbl /circumflex /perthousand /Scaron /guilsinglleft
+/OE /bullet /Zcaron /bullet /bullet
+/quoteleft /quoteright /quotedblleft /quotedblright /bullet
+/endash /emdash /asciitilde /trademark /scaron
+/guilsinglright /oe /bullet /zcaron /Ydieresis
+/space /exclamdown /cent /sterling /currency
+/yen /brokenbar /section /dieresis /copyright
+/ordfeminine /guillemotleft /logicalnot /sfthyphen /registered
+/macron /degree /plusminus /twosuperior /threesuperior
+/acute /mu /paragraph /middot /cedilla
+/onesuperior /ordmasculine /guillemotright /onequarter /onehalf
+/threequarters /questiondown /Agrave /Aacute /Acircumflex
+/Atilde /Adieresis /Aring /AE /Ccedilla
+/Egrave /Eacute /Ecircumflex /Edieresis /Igrave
+/Iacute /Icircumflex /Idieresis /Eth /Ntilde
+/Ograve /Oacute /Ocircumflex /Otilde /Odieresis
+/multiply /Oslash /Ugrave /Uacute /Ucircumflex
+/Udieresis /Yacute /Thorn /germandbls /agrave
+/aacute /acircumflex /atilde /adieresis /aring
+/ae /ccedilla /egrave /eacute /ecircumflex
+/edieresis /igrave /iacute /icircumflex /idieresis
+/eth /ntilde /ograve /oacute /ocircumflex
+/otilde /odieresis /divide /oslash /ugrave
+/uacute /ucircumflex /udieresis /yacute /thorn
+/ydieresis
+] def
+%%EndResource
+%FOPBeginFontReencode
+/Courier-Oblique findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding WinAnsiEncoding def
+ currentdict
+end
+/Courier-Oblique exch definefont pop
+/Courier-BoldOblique findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding WinAnsiEncoding def
+ currentdict
+end
+/Courier-BoldOblique exch definefont pop
+/Courier-Bold findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding WinAnsiEncoding def
+ currentdict
+end
+/Courier-Bold exch definefont pop
+/Helvetica findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding WinAnsiEncoding def
+ currentdict
+end
+/Helvetica exch definefont pop
+/Helvetica-Oblique findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding WinAnsiEncoding def
+ currentdict
+end
+/Helvetica-Oblique exch definefont pop
+/Helvetica-Bold findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding WinAnsiEncoding def
+ currentdict
+end
+/Helvetica-Bold exch definefont pop
+/Helvetica-BoldOblique findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding WinAnsiEncoding def
+ currentdict
+end
+/Helvetica-BoldOblique exch definefont pop
+/Times-Roman findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding WinAnsiEncoding def
+ currentdict
+end
+/Times-Roman exch definefont pop
+/Times-Italic findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding WinAnsiEncoding def
+ currentdict
+end
+/Times-Italic exch definefont pop
+/Times-Bold findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding WinAnsiEncoding def
+ currentdict
+end
+/Times-Bold exch definefont pop
+/Times-BoldItalic findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding WinAnsiEncoding def
+ currentdict
+end
+/Times-BoldItalic exch definefont pop
+/Courier findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding WinAnsiEncoding def
+ currentdict
+end
+/Courier exch definefont pop
+%FOPEndFontReencode
+%%EndProlog
+%%Page: 1 1
+%%PageBoundingBox: 0 0 420 315
+%%BeginPageSetup
+N
+ 18 4 M
+ 406 4 L
+ 406 308 L
+ 18 308 L
+cp
+clip
+[1 0 0 -1 0 315] CT
+%%EndPageSetup
+GS
+[0.3 0 0 0.3 0 0] CT
+1 GC
+N
+0 0 1400 1050 re
+f
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+1 GC
+N
+0 0 1400 1050 re
+f
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+1 GC
+N
+182 934 M
+1267 934 L
+1267 79 L
+182 79 L
+cp
+f
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+182 934 M
+182 79 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+290.5 934 M
+290.5 79 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+399 934 M
+399 79 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+507.5 934 M
+507.5 79 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+616 934 M
+616 79 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+724.5 934 M
+724.5 79 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+833 934 M
+833 79 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+941.5 934 M
+941.5 79 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+1050 934 M
+1050 79 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+1158.5 934 M
+1158.5 79 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+1267 934 M
+1267 79 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+1267 934 M
+182 934 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+1267 811.857 M
+182 811.857 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+1267 689.714 M
+182 689.714 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+1267 567.571 M
+182 567.571 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+1267 445.429 M
+182 445.429 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+1267 323.286 M
+182 323.286 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+1267 201.143 M
+182 201.143 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.873 GC
+1 LJ
+1.667 LW
+N
+1267 79 M
+182 79 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+182 934 M
+1267 934 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+182 934 M
+182 923.15 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+290.5 934 M
+290.5 923.15 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+399 934 M
+399 923.15 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+507.5 934 M
+507.5 923.15 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+616 934 M
+616 923.15 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+724.5 934 M
+724.5 923.15 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+833 934 M
+833 923.15 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+941.5 934 M
+941.5 923.15 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+1050 934 M
+1050 923.15 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+1158.5 934 M
+1158.5 923.15 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+1267 934 M
+1267 923.15 L
+S
+GR
+GS
+[0.3 0 0 0.3 54.6 284.19999] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-9.5 34 moveto
+1 -1 scale
+(0) t
+GR
+GR
+GS
+[0.3 0 0 0.3 87.15 284.19999] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-19 34 moveto
+1 -1 scale
+(20) t
+GR
+GR
+GS
+[0.3 0 0 0.3 119.7 284.19999] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-19 34 moveto
+1 -1 scale
+(40) t
+GR
+GR
+GS
+[0.3 0 0 0.3 152.25 284.19999] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-19 34 moveto
+1 -1 scale
+(60) t
+GR
+GR
+GS
+[0.3 0 0 0.3 184.8 284.19999] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-19 34 moveto
+1 -1 scale
+(80) t
+GR
+GR
+GS
+[0.3 0 0 0.3 217.35 284.19999] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-28 34 moveto
+1 -1 scale
+(100) t
+GR
+GR
+GS
+[0.3 0 0 0.3 249.9 284.19999] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-28 34 moveto
+1 -1 scale
+(120) t
+GR
+GR
+GS
+[0.3 0 0 0.3 282.45 284.19999] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-28 34 moveto
+1 -1 scale
+(140) t
+GR
+GR
+GS
+[0.3 0 0 0.3 315 284.19999] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-28 34 moveto
+1 -1 scale
+(160) t
+GR
+GR
+GS
+[0.3 0 0 0.3 347.55 284.19999] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-28 34 moveto
+1 -1 scale
+(180) t
+GR
+GR
+GS
+[0.3 0 0 0.3 380.1 284.19999] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-28 34 moveto
+1 -1 scale
+(200) t
+GR
+GR
+GS
+[0.3 0 0 0.3 217.35015 297.80001] CT
+0.149 GC
+/Helvetica 40 F
+GS
+[1 0 0 1 0 0] CT
+-71.5 41 moveto
+1 -1 scale
+(iteration) t
+GR
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+182 934 M
+182 79 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+182 934 M
+192.85 934 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+182 811.857 M
+192.85 811.857 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+182 689.714 M
+192.85 689.714 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+182 567.571 M
+192.85 567.571 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+182 445.429 M
+192.85 445.429 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+182 323.286 M
+192.85 323.286 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+182 201.143 M
+192.85 201.143 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+2 setlinecap
+1 LJ
+1.667 LW
+N
+182 79 M
+192.85 79 L
+S
+GR
+GS
+[0.3 0 0 0.3 50.6 280.2] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-58 13 moveto
+1 -1 scale
+(-1.6) t
+GR
+GR
+GS
+[0.3 0 0 0.3 50.6 243.55714] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-58 13 moveto
+1 -1 scale
+(-1.4) t
+GR
+GR
+GS
+[0.3 0 0 0.3 50.6 206.91429] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-58 13 moveto
+1 -1 scale
+(-1.2) t
+GR
+GR
+GS
+[0.3 0 0 0.3 50.6 170.27142] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-30 13 moveto
+1 -1 scale
+(-1) t
+GR
+GR
+GS
+[0.3 0 0 0.3 50.6 133.62858] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-58 13 moveto
+1 -1 scale
+(-0.8) t
+GR
+GR
+GS
+[0.3 0 0 0.3 50.6 96.98572] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-58 13 moveto
+1 -1 scale
+(-0.6) t
+GR
+GR
+GS
+[0.3 0 0 0.3 50.6 60.34286] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-58 13 moveto
+1 -1 scale
+(-0.4) t
+GR
+GR
+GS
+[0.3 0 0 0.3 50.6 23.7] CT
+0.149 GC
+/Helvetica 33.333 F
+GS
+[1 0 0 1 0 0] CT
+-58 13 moveto
+1 -1 scale
+(-0.2) t
+GR
+GR
+GS
+[0 -0.3 0.3 0 30.2 151.94987] CT
+0.149 GC
+/Helvetica 40 F
+GS
+[1 0 0 1 0 0] CT
+-104.5 -9 moveto
+1 -1 scale
+(log10\(error\)) t
+GR
+GR
+GS
+[0.3 0 0 0.3 217.35022 22.8] CT
+/Helvetica-Bold 40 F
+GS
+[1 0 0 1 0 0] CT
+-52.5 -9 moveto
+1 -1 scale
+(musk) t
+GR
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0 0.447 0.741 RC
+1 LJ
+6.667 LW
+N
+187.425 81.373 M
+192.85 376.397 L
+198.275 480.538 L
+203.7 516.821 L
+209.125 538.059 L
+214.55 554.589 L
+219.975 566.648 L
+225.4 572.919 L
+230.825 584.201 L
+236.25 595.283 L
+241.675 607.289 L
+247.1 615.136 L
+252.525 620.785 L
+257.95 627.775 L
+263.375 633.38 L
+268.8 637.722 L
+274.225 640.328 L
+279.65 648.322 L
+285.075 649.552 L
+290.5 656.415 L
+295.925 663.018 L
+301.35 666.936 L
+306.775 674.196 L
+312.2 678.009 L
+317.625 681.848 L
+323.05 684.637 L
+328.475 685.549 L
+333.9 687.504 L
+339.325 690.576 L
+344.75 694.518 L
+350.175 694.242 L
+355.6 697.497 L
+361.025 702.392 L
+366.45 703.41 L
+371.875 705.767 L
+377.3 704.705 L
+382.725 709.031 L
+388.15 711.608 L
+393.575 711.562 L
+399 712.844 L
+404.425 717.301 L
+409.85 723.875 L
+415.275 728.184 L
+420.7 729.812 L
+426.125 732.763 L
+431.55 731.143 L
+436.975 731.683 L
+442.4 734.248 L
+447.825 734.662 L
+453.25 734.869 L
+458.675 737.828 L
+464.1 737.112 L
+469.525 739.814 L
+474.95 742.185 L
+480.375 739.191 L
+485.8 736.793 L
+491.225 736.998 L
+496.65 736.684 L
+502.075 737.869 L
+507.5 740.904 L
+512.925 737.527 L
+518.35 741.534 L
+523.775 748.768 L
+529.2 749.75 L
+534.625 750.64 L
+540.05 751.231 L
+545.475 749.948 L
+550.9 750.481 L
+556.325 753.048 L
+561.75 755.451 L
+567.175 757.761 L
+572.6 760.23 L
+578.025 762.263 L
+583.45 765.641 L
+588.875 763.47 L
+594.3 762.574 L
+599.725 761.147 L
+605.15 760.901 L
+610.575 767.052 L
+616 769.503 L
+621.425 772.101 L
+626.85 774.146 L
+632.275 771.363 L
+637.7 773.883 L
+643.125 771.471 L
+648.55 773.904 L
+653.975 779.732 L
+659.4 779.904 L
+664.825 776.973 L
+670.25 779.246 L
+675.675 778.792 L
+681.1 780.018 L
+686.525 783.556 L
+691.95 781.595 L
+697.375 779.428 L
+702.8 782.161 L
+708.225 785.196 L
+713.65 786.415 L
+719.075 782.723 L
+724.5 784.279 L
+729.925 781.294 L
+735.35 782.117 L
+740.775 785.948 L
+746.2 783.902 L
+751.625 787.814 L
+757.05 788.342 L
+762.475 791.375 L
+767.9 791.931 L
+773.325 790.413 L
+778.75 792.001 L
+784.175 788.758 L
+789.6 789.969 L
+795.025 793.365 L
+800.45 794.205 L
+805.875 793.334 L
+811.3 797.607 L
+816.725 799.775 L
+822.15 800.973 L
+827.575 800.322 L
+833 798.875 L
+838.425 804.079 L
+843.85 800.927 L
+849.275 800.752 L
+854.7 800.266 L
+860.125 800.234 L
+865.55 806.741 L
+870.975 809.502 L
+876.4 808.699 L
+881.825 810.421 L
+887.25 812.044 L
+892.675 815.177 L
+898.1 812.741 L
+903.525 816.025 L
+908.95 815.67 L
+914.375 816.935 L
+919.8 815.677 L
+925.225 814.181 L
+930.65 818.054 L
+936.075 820.586 L
+941.5 818.657 L
+946.925 823.746 L
+952.35 823.531 L
+957.775 829.665 L
+963.2 825.796 L
+968.625 826.736 L
+974.05 827.832 L
+979.475 830.722 L
+984.9 830.251 L
+990.325 829.368 L
+995.75 831.22 L
+1001.175 829.669 L
+1006.6 833.057 L
+1012.025 836.294 L
+1017.45 833.687 L
+1022.875 834.771 L
+1028.3 834.352 L
+1033.725 835.85 L
+1039.15 836.733 L
+1044.575 838.564 L
+1050 839.807 L
+1055.425 837.524 L
+1060.85 839.064 L
+1066.275 841.039 L
+1071.7 844.409 L
+1077.125 846.42 L
+1082.55 844.041 L
+1087.975 846.785 L
+1093.4 849.453 L
+1098.825 848.551 L
+1104.25 846.61 L
+1109.675 847.625 L
+1115.1 849.693 L
+1120.525 851.489 L
+1125.95 851.85 L
+1131.375 852.596 L
+1136.8 850.245 L
+1142.225 853.117 L
+1147.65 851.102 L
+1153.075 853.116 L
+1158.5 856.584 L
+1163.925 859.125 L
+1169.35 857.581 L
+1174.775 859.201 L
+1180.2 859.119 L
+1185.625 860.576 L
+1191.05 862.678 L
+1196.475 865.282 L
+1201.9 862.701 L
+1207.325 860.253 L
+1212.75 860.906 L
+1218.175 859.882 L
+1223.6 862.518 L
+1229.025 861.444 L
+1234.45 863.464 L
+1239.875 862.883 L
+1245.3 864.336 L
+1250.725 862.253 L
+1256.15 857.542 L
+1261.575 860.828 L
+1267 861.737 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.851 0.325 0.098 RC
+1 LJ
+6.667 LW
+N
+187.425 635.547 M
+192.85 635.547 L
+198.275 635.547 L
+203.7 635.547 L
+209.125 635.547 L
+214.55 635.547 L
+219.975 635.547 L
+225.4 635.547 L
+230.825 635.547 L
+236.25 635.547 L
+241.675 635.547 L
+247.1 635.547 L
+252.525 635.547 L
+257.95 635.547 L
+263.375 635.547 L
+268.8 635.547 L
+274.225 635.547 L
+279.65 635.547 L
+285.075 635.547 L
+290.5 635.547 L
+295.925 635.547 L
+301.35 635.547 L
+306.775 635.547 L
+312.2 635.547 L
+317.625 635.547 L
+323.05 635.547 L
+328.475 635.547 L
+333.9 635.547 L
+339.325 635.547 L
+344.75 635.547 L
+350.175 635.547 L
+355.6 635.547 L
+361.025 635.547 L
+366.45 635.547 L
+371.875 635.547 L
+377.3 635.547 L
+382.725 635.547 L
+388.15 635.547 L
+393.575 635.547 L
+399 635.547 L
+404.425 635.547 L
+409.85 635.547 L
+415.275 635.547 L
+420.7 635.547 L
+426.125 635.547 L
+431.55 635.547 L
+436.975 635.547 L
+442.4 635.547 L
+447.825 635.547 L
+453.25 635.547 L
+458.675 635.547 L
+464.1 635.547 L
+469.525 635.547 L
+474.95 635.547 L
+480.375 635.547 L
+485.8 635.547 L
+491.225 635.547 L
+496.65 635.547 L
+502.075 635.547 L
+507.5 635.547 L
+512.925 635.547 L
+518.35 635.547 L
+523.775 635.547 L
+529.2 635.547 L
+534.625 635.547 L
+540.05 635.547 L
+545.475 635.547 L
+550.9 635.547 L
+556.325 635.547 L
+561.75 635.547 L
+567.175 635.547 L
+572.6 635.547 L
+578.025 635.547 L
+583.45 635.547 L
+588.875 635.547 L
+594.3 635.547 L
+599.725 635.547 L
+605.15 635.547 L
+610.575 635.547 L
+616 635.547 L
+621.425 635.547 L
+626.85 635.547 L
+632.275 635.547 L
+637.7 635.547 L
+643.125 635.547 L
+648.55 635.547 L
+653.975 635.547 L
+659.4 635.547 L
+664.825 635.547 L
+670.25 635.547 L
+675.675 635.547 L
+681.1 635.547 L
+686.525 635.547 L
+691.95 635.547 L
+697.375 635.547 L
+702.8 635.547 L
+708.225 635.547 L
+713.65 635.547 L
+719.075 635.547 L
+724.5 635.547 L
+729.925 635.547 L
+735.35 635.547 L
+740.775 635.547 L
+746.2 635.547 L
+751.625 635.547 L
+757.05 635.547 L
+762.475 635.547 L
+767.9 635.547 L
+773.325 635.547 L
+778.75 635.547 L
+784.175 635.547 L
+789.6 635.547 L
+795.025 635.547 L
+800.45 635.547 L
+805.875 635.547 L
+811.3 635.547 L
+816.725 635.547 L
+822.15 635.547 L
+827.575 635.547 L
+833 635.547 L
+838.425 635.547 L
+843.85 635.547 L
+849.275 635.547 L
+854.7 635.547 L
+860.125 635.547 L
+865.55 635.547 L
+870.975 635.547 L
+876.4 635.547 L
+881.825 635.547 L
+887.25 635.547 L
+892.675 635.547 L
+898.1 635.547 L
+903.525 635.547 L
+908.95 635.547 L
+914.375 635.547 L
+919.8 635.547 L
+925.225 635.547 L
+930.65 635.547 L
+936.075 635.547 L
+941.5 635.547 L
+946.925 635.547 L
+952.35 635.547 L
+957.775 635.547 L
+963.2 635.547 L
+968.625 635.547 L
+974.05 635.547 L
+979.475 635.547 L
+984.9 635.547 L
+990.325 635.547 L
+995.75 635.547 L
+1001.175 635.547 L
+1006.6 635.547 L
+1012.025 635.547 L
+1017.45 635.547 L
+1022.875 635.547 L
+1028.3 635.547 L
+1033.725 635.547 L
+1039.15 635.547 L
+1044.575 635.547 L
+1050 635.547 L
+1055.425 635.547 L
+1060.85 635.547 L
+1066.275 635.547 L
+1071.7 635.547 L
+1077.125 635.547 L
+1082.55 635.547 L
+1087.975 635.547 L
+1093.4 635.547 L
+1098.825 635.547 L
+1104.25 635.547 L
+1109.675 635.547 L
+1115.1 635.547 L
+1120.525 635.547 L
+1125.95 635.547 L
+1131.375 635.547 L
+1136.8 635.547 L
+1142.225 635.547 L
+1147.65 635.547 L
+1153.075 635.547 L
+1158.5 635.547 L
+1163.925 635.547 L
+1169.35 635.547 L
+1174.775 635.547 L
+1180.2 635.547 L
+1185.625 635.547 L
+1191.05 635.547 L
+1196.475 635.547 L
+1201.9 635.547 L
+1207.325 635.547 L
+1212.75 635.547 L
+1218.175 635.547 L
+1223.6 635.547 L
+1229.025 635.547 L
+1234.45 635.547 L
+1239.875 635.547 L
+1245.3 635.547 L
+1250.725 635.547 L
+1256.15 635.547 L
+1261.575 635.547 L
+1267 635.547 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.929 0.694 0.125 RC
+1 LJ
+6.667 LW
+N
+187.425 700.862 M
+192.85 700.862 L
+198.275 700.862 L
+203.7 700.862 L
+209.125 700.862 L
+214.55 700.862 L
+219.975 700.862 L
+225.4 700.862 L
+230.825 700.862 L
+236.25 700.862 L
+241.675 700.862 L
+247.1 700.862 L
+252.525 700.862 L
+257.95 700.862 L
+263.375 700.862 L
+268.8 700.862 L
+274.225 700.862 L
+279.65 700.862 L
+285.075 700.862 L
+290.5 700.862 L
+295.925 700.862 L
+301.35 700.862 L
+306.775 700.862 L
+312.2 700.862 L
+317.625 700.862 L
+323.05 700.862 L
+328.475 700.862 L
+333.9 700.862 L
+339.325 700.862 L
+344.75 700.862 L
+350.175 700.862 L
+355.6 700.862 L
+361.025 700.862 L
+366.45 700.862 L
+371.875 700.862 L
+377.3 700.862 L
+382.725 700.862 L
+388.15 700.862 L
+393.575 700.862 L
+399 700.862 L
+404.425 700.862 L
+409.85 700.862 L
+415.275 700.862 L
+420.7 700.862 L
+426.125 700.862 L
+431.55 700.862 L
+436.975 700.862 L
+442.4 700.862 L
+447.825 700.862 L
+453.25 700.862 L
+458.675 700.862 L
+464.1 700.862 L
+469.525 700.862 L
+474.95 700.862 L
+480.375 700.862 L
+485.8 700.862 L
+491.225 700.862 L
+496.65 700.862 L
+502.075 700.862 L
+507.5 700.862 L
+512.925 700.862 L
+518.35 700.862 L
+523.775 700.862 L
+529.2 700.862 L
+534.625 700.862 L
+540.05 700.862 L
+545.475 700.862 L
+550.9 700.862 L
+556.325 700.862 L
+561.75 700.862 L
+567.175 700.862 L
+572.6 700.862 L
+578.025 700.862 L
+583.45 700.862 L
+588.875 700.862 L
+594.3 700.862 L
+599.725 700.862 L
+605.15 700.862 L
+610.575 700.862 L
+616 700.862 L
+621.425 700.862 L
+626.85 700.862 L
+632.275 700.862 L
+637.7 700.862 L
+643.125 700.862 L
+648.55 700.862 L
+653.975 700.862 L
+659.4 700.862 L
+664.825 700.862 L
+670.25 700.862 L
+675.675 700.862 L
+681.1 700.862 L
+686.525 700.862 L
+691.95 700.862 L
+697.375 700.862 L
+702.8 700.862 L
+708.225 700.862 L
+713.65 700.862 L
+719.075 700.862 L
+724.5 700.862 L
+729.925 700.862 L
+735.35 700.862 L
+740.775 700.862 L
+746.2 700.862 L
+751.625 700.862 L
+757.05 700.862 L
+762.475 700.862 L
+767.9 700.862 L
+773.325 700.862 L
+778.75 700.862 L
+784.175 700.862 L
+789.6 700.862 L
+795.025 700.862 L
+800.45 700.862 L
+805.875 700.862 L
+811.3 700.862 L
+816.725 700.862 L
+822.15 700.862 L
+827.575 700.862 L
+833 700.862 L
+838.425 700.862 L
+843.85 700.862 L
+849.275 700.862 L
+854.7 700.862 L
+860.125 700.862 L
+865.55 700.862 L
+870.975 700.862 L
+876.4 700.862 L
+881.825 700.862 L
+887.25 700.862 L
+892.675 700.862 L
+898.1 700.862 L
+903.525 700.862 L
+908.95 700.862 L
+914.375 700.862 L
+919.8 700.862 L
+925.225 700.862 L
+930.65 700.862 L
+936.075 700.862 L
+941.5 700.862 L
+946.925 700.862 L
+952.35 700.862 L
+957.775 700.862 L
+963.2 700.862 L
+968.625 700.862 L
+974.05 700.862 L
+979.475 700.862 L
+984.9 700.862 L
+990.325 700.862 L
+995.75 700.862 L
+1001.175 700.862 L
+1006.6 700.862 L
+1012.025 700.862 L
+1017.45 700.862 L
+1022.875 700.862 L
+1028.3 700.862 L
+1033.725 700.862 L
+1039.15 700.862 L
+1044.575 700.862 L
+1050 700.862 L
+1055.425 700.862 L
+1060.85 700.862 L
+1066.275 700.862 L
+1071.7 700.862 L
+1077.125 700.862 L
+1082.55 700.862 L
+1087.975 700.862 L
+1093.4 700.862 L
+1098.825 700.862 L
+1104.25 700.862 L
+1109.675 700.862 L
+1115.1 700.862 L
+1120.525 700.862 L
+1125.95 700.862 L
+1131.375 700.862 L
+1136.8 700.862 L
+1142.225 700.862 L
+1147.65 700.862 L
+1153.075 700.862 L
+1158.5 700.862 L
+1163.925 700.862 L
+1169.35 700.862 L
+1174.775 700.862 L
+1180.2 700.862 L
+1185.625 700.862 L
+1191.05 700.862 L
+1196.475 700.862 L
+1201.9 700.862 L
+1207.325 700.862 L
+1212.75 700.862 L
+1218.175 700.862 L
+1223.6 700.862 L
+1229.025 700.862 L
+1234.45 700.862 L
+1239.875 700.862 L
+1245.3 700.862 L
+1250.725 700.862 L
+1256.15 700.862 L
+1261.575 700.862 L
+1267 700.862 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.494 0.184 0.557 RC
+1 LJ
+6.667 LW
+N
+187.425 689.944 M
+192.85 689.944 L
+198.275 689.944 L
+203.7 689.944 L
+209.125 689.944 L
+214.55 689.944 L
+219.975 689.944 L
+225.4 689.944 L
+230.825 689.944 L
+236.25 689.944 L
+241.675 689.944 L
+247.1 689.944 L
+252.525 689.944 L
+257.95 689.944 L
+263.375 689.944 L
+268.8 689.944 L
+274.225 689.944 L
+279.65 689.944 L
+285.075 689.944 L
+290.5 689.944 L
+295.925 689.944 L
+301.35 689.944 L
+306.775 689.944 L
+312.2 689.944 L
+317.625 689.944 L
+323.05 689.944 L
+328.475 689.944 L
+333.9 689.944 L
+339.325 689.944 L
+344.75 689.944 L
+350.175 689.944 L
+355.6 689.944 L
+361.025 689.944 L
+366.45 689.944 L
+371.875 689.944 L
+377.3 689.944 L
+382.725 689.944 L
+388.15 689.944 L
+393.575 689.944 L
+399 689.944 L
+404.425 689.944 L
+409.85 689.944 L
+415.275 689.944 L
+420.7 689.944 L
+426.125 689.944 L
+431.55 689.944 L
+436.975 689.944 L
+442.4 689.944 L
+447.825 689.944 L
+453.25 689.944 L
+458.675 689.944 L
+464.1 689.944 L
+469.525 689.944 L
+474.95 689.944 L
+480.375 689.944 L
+485.8 689.944 L
+491.225 689.944 L
+496.65 689.944 L
+502.075 689.944 L
+507.5 689.944 L
+512.925 689.944 L
+518.35 689.944 L
+523.775 689.944 L
+529.2 689.944 L
+534.625 689.944 L
+540.05 689.944 L
+545.475 689.944 L
+550.9 689.944 L
+556.325 689.944 L
+561.75 689.944 L
+567.175 689.944 L
+572.6 689.944 L
+578.025 689.944 L
+583.45 689.944 L
+588.875 689.944 L
+594.3 689.944 L
+599.725 689.944 L
+605.15 689.944 L
+610.575 689.944 L
+616 689.944 L
+621.425 689.944 L
+626.85 689.944 L
+632.275 689.944 L
+637.7 689.944 L
+643.125 689.944 L
+648.55 689.944 L
+653.975 689.944 L
+659.4 689.944 L
+664.825 689.944 L
+670.25 689.944 L
+675.675 689.944 L
+681.1 689.944 L
+686.525 689.944 L
+691.95 689.944 L
+697.375 689.944 L
+702.8 689.944 L
+708.225 689.944 L
+713.65 689.944 L
+719.075 689.944 L
+724.5 689.944 L
+729.925 689.944 L
+735.35 689.944 L
+740.775 689.944 L
+746.2 689.944 L
+751.625 689.944 L
+757.05 689.944 L
+762.475 689.944 L
+767.9 689.944 L
+773.325 689.944 L
+778.75 689.944 L
+784.175 689.944 L
+789.6 689.944 L
+795.025 689.944 L
+800.45 689.944 L
+805.875 689.944 L
+811.3 689.944 L
+816.725 689.944 L
+822.15 689.944 L
+827.575 689.944 L
+833 689.944 L
+838.425 689.944 L
+843.85 689.944 L
+849.275 689.944 L
+854.7 689.944 L
+860.125 689.944 L
+865.55 689.944 L
+870.975 689.944 L
+876.4 689.944 L
+881.825 689.944 L
+887.25 689.944 L
+892.675 689.944 L
+898.1 689.944 L
+903.525 689.944 L
+908.95 689.944 L
+914.375 689.944 L
+919.8 689.944 L
+925.225 689.944 L
+930.65 689.944 L
+936.075 689.944 L
+941.5 689.944 L
+946.925 689.944 L
+952.35 689.944 L
+957.775 689.944 L
+963.2 689.944 L
+968.625 689.944 L
+974.05 689.944 L
+979.475 689.944 L
+984.9 689.944 L
+990.325 689.944 L
+995.75 689.944 L
+1001.175 689.944 L
+1006.6 689.944 L
+1012.025 689.944 L
+1017.45 689.944 L
+1022.875 689.944 L
+1028.3 689.944 L
+1033.725 689.944 L
+1039.15 689.944 L
+1044.575 689.944 L
+1050 689.944 L
+1055.425 689.944 L
+1060.85 689.944 L
+1066.275 689.944 L
+1071.7 689.944 L
+1077.125 689.944 L
+1082.55 689.944 L
+1087.975 689.944 L
+1093.4 689.944 L
+1098.825 689.944 L
+1104.25 689.944 L
+1109.675 689.944 L
+1115.1 689.944 L
+1120.525 689.944 L
+1125.95 689.944 L
+1131.375 689.944 L
+1136.8 689.944 L
+1142.225 689.944 L
+1147.65 689.944 L
+1153.075 689.944 L
+1158.5 689.944 L
+1163.925 689.944 L
+1169.35 689.944 L
+1174.775 689.944 L
+1180.2 689.944 L
+1185.625 689.944 L
+1191.05 689.944 L
+1196.475 689.944 L
+1201.9 689.944 L
+1207.325 689.944 L
+1212.75 689.944 L
+1218.175 689.944 L
+1223.6 689.944 L
+1229.025 689.944 L
+1234.45 689.944 L
+1239.875 689.944 L
+1245.3 689.944 L
+1250.725 689.944 L
+1256.15 689.944 L
+1261.575 689.944 L
+1267 689.944 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+1 GC
+N
+1240 319 M
+1240 106 L
+579 106 L
+579 319 L
+cp
+f
+GR
+GS
+[0.3 0 0 0.3 208.92336 40.76842] CT
+/Helvetica 40 F
+GS
+[1 0 0 1 0 0] CT
+0 16 moveto
+1 -1 scale
+(Cone-DL) t
+GR
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0 0.447 0.741 RC
+1 LJ
+6.667 LW
+N
+588.992 135.895 M
+688.917 135.895 L
+S
+GR
+GS
+[0.3 0 0 0.3 208.92336 56.08947] CT
+/Helvetica 40 F
+GS
+[1 0 0 1 0 0] CT
+0 16 moveto
+1 -1 scale
+(Cone-DL + swap) t
+GR
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.851 0.325 0.098 RC
+1 LJ
+6.667 LW
+N
+588.992 186.965 M
+688.917 186.965 L
+S
+GR
+GS
+[0.3 0 0 0.3 208.92336 71.41052] CT
+/Helvetica 40 F
+GS
+[1 0 0 1 0 0] CT
+0 16 moveto
+1 -1 scale
+(AK-SVD + Cone-OMP) t
+GR
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.929 0.694 0.125 RC
+1 LJ
+6.667 LW
+N
+588.992 238.035 M
+688.917 238.035 L
+S
+GR
+GS
+[0.3 0 0 0.3 208.92336 86.73158] CT
+/Helvetica 40 F
+GS
+[1 0 0 1 0 0] CT
+0 16 moveto
+1 -1 scale
+(AK-SVD + Cone-OMP + swap) t
+GR
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.494 0.184 0.557 RC
+1 LJ
+6.667 LW
+N
+588.992 289.105 M
+688.917 289.105 L
+S
+GR
+GS
+[0.3 0 0 0.3 0 0] CT
+0.149 GC
+10.0 ML
+1.667 LW
+N
+579 319 M
+579 106 L
+1240 106 L
+1240 319 L
+cp
+S
+GR
+%%Trailer
+%%Pages: 1
+%%EOF
diff --git a/omp_cone_single.m b/omp_cone_single.m
new file mode 100644
index 0000000000000000000000000000000000000000..ea69e26af7f984a1a02274c6dd6b8f1ecd0ee143
--- /dev/null
+++ b/omp_cone_single.m
@@ -0,0 +1,116 @@
+function [x, Da, support] = omp_cone_single(y, D, s, rad, iter_cd, err_thresh)
+
+% Optimal Matching Pursuit with cone atoms
+%
+% Input:
+% y - signal
+% D - dictionary (normalized atoms)
+% s - sparsity level
+% rad - cone radii (if scalar, then all atoms have the same radius)
+% iter_cd - number of iterations in the coordinate descent refinement
+% (default = 1)
+% err_thresh- representation error for which OMP is stopped
+%
+% Output:
+% X - sparse representations matrix
+% Da - actual dictionary for last representation (meaningful when N=1)
+% support - indices of atoms used in representation
+
+% BD 14.06.2022
+
+if nargin < 5
+ iter_cd = 1;
+end
+
+if nargin < 6
+ err_thresh = 1e-14;
+end
+
+% sizes and constants
+[m,n] = size(D);
+if length(rad) == 1
+ rad = rad * ones(n,1);
+end
+rad2 = rad.*rad/2;
+lam = 1 - rad2;
+
+x = zeros(s,1);
+Da = zeros(m,s);
+
+r = y; % init residual
+support = [];
+for k = 1 : s
+ rn = r / norm(r); % normalized residual, to be used in convex combination with best atom
+ % find best current atom
+ proj = D'*rn; % projections, like in OMP
+ proj(support) = 0; % avoid atom duplication
+ [~,jmax] = max(lam.*abs(proj) + sqrt((1-lam.*lam).*(1-proj.*proj))); % find nearest cone
+ if proj(jmax) < 0
+ rn = -rn;
+ end
+ p = abs(proj(jmax));
+ if p + rad2(jmax) >= 1 % residual is inside cone
+ dc = rn; % hence new atom is the residual
+ res_inside_cone = 1;
+ else % residual outside cone
+ b = sqrt(1/(1-p*p));
+ q = -b*p*D(:,jmax) + b*rn; % vector orthogonal on atom, in the plane atom-residual
+ a = 1 - rad2(jmax);
+ dc = a*D(:,jmax) + sqrt(1-a*a)*q; % the new atom
+ dc = dc/norm(dc); % normalize, to be sure (although it should have norm 1)
+ res_inside_cone = 0;
+ end
+ support = [support jmax]; % update support with new atom
+
+ % compute (nearly) optimal representation
+ Da(:,k) = dc; % current dictionary
+ if res_inside_cone % with new atom, residual becomes zero
+ x(k) = dc'*r;
+ x = x(1:k); % cut representation to useful part
+ Da = Da(:,1:k);
+ return
+ else
+ Dk = Da(:,1:k);
+ xk = Dk \ y; % least square solution, with current atoms
+ if k > 1
+ for icd = 1 : iter_cd
+ r = y - Dk*xk;
+ %norm(r)/norm(y)
+ for i = 1 : k % a round of coordinate descent
+ r = r + Dk(:,i)*xk(i);
+ rn = r / norm(r);
+ p = D(:,support(i))'*rn;
+ if p<0
+ rn = -rn;
+ end
+ p = abs(p);
+ if p + rad2(support(i)) >= 1 % residual in the cone => ready
+ Dk(:,i) = rn;
+ xk(i) = rn'*r;
+ x = xk;
+ Da = Dk;
+ return
+ else
+ b = sqrt(1/(1-p*p));
+ q = -b*p*D(:,support(i)) + b*rn;
+ a = 1 - rad2(support(i));
+ dc = a*D(:,support(i)) + sqrt(1-a*a)*q; % the new atom
+ dc = dc/norm(dc); % ???
+ Dk(:,i) = dc;
+ xk(i) = dc'*r;
+ r = r - dc*xk(i);
+ end
+ end
+ end
+ end
+ r = y - Dk*xk;
+ err = norm(r)/norm(y);
+ x(1:k) = xk;
+ Da(:,1:k) = Dk;
+ if err < err_thresh
+ x = x(1:k);
+ Da = Da(:,1:k);
+ return
+ end
+ end
+end
diff --git a/test_all_dl_ad.m b/test_all_dl_ad.m
new file mode 100644
index 0000000000000000000000000000000000000000..b76b05373a43b255d2315a5e76d530e63504c641
--- /dev/null
+++ b/test_all_dl_ad.m
@@ -0,0 +1,223 @@
+% test several DL algorithms on PYOD benchmark datasets
+
+% BD 16.05.2023
+
+%% ----------------------------------------
+data_file_path = '../DATA_DEPENDENCY/dependency_';
+% data_file_names = {'musk', 'SpamBase', 'landsat', 'satellite', 'satimage', ...
+% 'ionosphere', 'WPBC', 'letter', 'WDBC', 'fault', ...
+% 'cardio', 'Waveform', 'Hepatitis', 'Lymphography', 'pendigits', ...
+% 'wine', 'vowels', 'cover', 'magic.gamma.mat', 'breastw', ...
+% 'Stamps', 'WBC', 'yeast', 'glass', 'annthyroid', ...
+% 'mammography', 'Pima', 'thyroid', 'vertebral', 'Wilt'};
+data_file_names = {'landsat', 'cardio'}; % uncomment full list above for more datasets
+results_file_name = ['res_all_odds_', int2str(randi(9999)), '.csv'];
+
+n_tests = 2; % 10
+train_data_ratio = 0.7;
+
+overcompleteness_list = [2 3]; % overcompleteness factor of the dictionary
+s_list = [2];
+
+rad_max = 0.05; %0.1;
+rad_min = 0.05; %0.01;
+iternum = 100;
+iter_cd = 10;
+err_thresh = 1e-7;
+min_dist_cones = 0.01;
+
+if rad_min ~= rad_max
+ swap_rad = 1;
+ method_names = {'AK-SVD + OMP', 'Cone-DL', 'AK-SVD + Cone-OMP', 'Cone-DL + swap', 'AK-SVD + Cone-OMP + swap'};
+else
+ swap_rad = 0;
+ method_names = {'AK-SVD + OMP', 'Cone-DL', 'AK-SVD + Cone-OMP'};
+end
+
+% ----------------------------------------
+n_datasets = length(data_file_names);
+n_methods = length(method_names);
+n_params = length(overcompleteness_list)*length(s_list);
+
+% prepare file for results; write header
+fid_res = fopen(results_file_name,'w');
+fprintf(fid_res, "%s,%s,", 'Dataset', 'Method');
+for i_over = 1:length(overcompleteness_list)
+ for i_s = 1 : length(s_list)
+ fprintf(fid_res, "%s,", ['o=',num2str(overcompleteness_list(i_over)), ' s=', int2str(s_list(i_s))]);
+ end
+end
+fprintf(fid_res, "\n");
+
+for i=1:n_methods
+ res_methods{i} = zeros(n_datasets,n_params);
+end
+
+for id = 1 : n_datasets
+ fprintf("%s\n", data_file_names{id})
+ filename = [data_file_path, data_file_names{id}];
+ load(filename)
+ X = double(Y); % data are in X, labels are in y
+ y = labels;
+ clear Y, clear labels
+ X = X'; % signals are on columns
+ [m,N] = size(X);
+
+ % preprocessing: Z-score normalization
+ vm = mean(X,2);
+ X = double(X) - repmat(vm,1,N); % subtract mean
+ sm = std(X,0,2);
+ X = X ./ repmat(sm,1,N);
+
+ col_normal = find(y==0);
+ col_anomaly = find(y~=0);
+ Nanomaly = length(col_anomaly);
+ Nnormal = N - Nanomaly;
+ anomaly_ratio = Nanomaly/N;
+
+ Ntrain = round(train_data_ratio*N);
+ Ntest = N - Ntrain;
+ N_anomaly_train = round(anomaly_ratio*Ntrain);
+ N_normal_train = Ntrain - N_anomaly_train;
+ N_anomaly_test = Nanomaly - N_anomaly_train;
+ N_normal_test = Nnormal - N_normal_train;
+
+ for im = 1:n_methods
+ res_data{im} = zeros(n_tests, n_params);
+ res_err{im} = zeros(n_tests, n_params);
+ end
+
+ for it = 1 : n_tests
+ % split train and test data, respecting proportion of anomaly
+ p1 = randperm(Nnormal);
+ p2 = randperm(Nanomaly);
+ Ytrain = X(:,[col_normal(p1(1:N_normal_train)); col_anomaly(p2(1:N_anomaly_train))]);
+ Ytest = X(:,[col_normal(p1(N_normal_train+1:end)); col_anomaly(p2(N_anomaly_train+1:end))]);
+ ytest = [zeros(N_normal_test,1); y(col_anomaly(p2(N_anomaly_train+1:end)))];
+ N_anomaly_test_w = sum(ytest);
+
+ i_param = 0;
+
+ for i_over = 1:length(overcompleteness_list)
+ n = round(m*overcompleteness_list(i_over));
+ rad_unif = rad_min + (rad_max-rad_min) * rand(n,1);
+ dist_ok = 0;
+ i0 = 0;
+ while ~dist_ok
+ D0 = normc(randn(m,n));
+ i1 = check_cone_superposition(D0, rad_unif, min_dist_cones);
+ dist_ok = isempty(i1);
+ i0 = i0+1;
+ if i0 == 10 && dist_ok == 0
+ %error('Cannot generate dictionary with imposed distance')
+ rad_unif = 0.95*rad_unif;
+ i0=0;
+ end
+ end
+ %save('D0last', 'D0', 'rad_unif');
+ %max(rad_unif)
+
+ for i_s = 1 : length(s_list)
+ i_param = i_param+1;
+ s = s_list(i_s);
+ rad = rad_unif;
+
+ % training phase
+ % standard DL (AK-SVD)
+ D1 = DL(Ytrain, D0, s, iternum);
+
+ % Cone-DL
+ D2 = dl_cone_paksvd_opt(Ytrain, D0, s, iternum, rad, iter_cd, err_thresh, min_dist_cones);
+
+ % testing phase
+ % AK-SVD + OMP
+ Xomp = omp(Ytest, D1, s);
+ err = Ytest - D1*Xomp;
+ err = sqrt(sum(err.*err));
+ res_err{1}(it,i_param) = mean(err);
+ [~,~,~,roc_auc] = perfcurve(ytest,err,1);
+ res_data{1}(it,i_param) = roc_auc;
+
+ % AK-SVD + Cone-OMP
+ rad1 = adjust_radii(D1, rad, min_dist_cones); % ensure that cones are disjoint
+ err = zeros(Ntest,1);
+ for k = 1 : Ntest
+ [x, Da, support] = omp_cone_single(Ytest(:,k), D1, s, rad1, iter_cd, err_thresh);
+ err(k) = norm(Ytest(:,k) - Da*x);
+ end
+ res_err{3}(it,i_param) = mean(err);
+ [~,~,~,roc_auc] = perfcurve(ytest,err,1);
+ res_data{3}(it,i_param) = roc_auc;
+
+ % AK-SVD + Cone-OMP + swap
+ if swap_rad
+ % swap radii
+ [~, Xs] = compute_rep_cone(D1, Ytrain, s, rad1, iter_cd, err_thresh);
+ atom_use = sum(squeeze(Xs)~=0,2);
+ [~,Ip] = sort(atom_use,'descend');
+ [~,Ir] = sort(rad,'descend');
+ rad_swap = zeros(n,1);
+ rad_swap(Ip) = rad(Ir);
+ rad_swap = adjust_radii(D1, rad_swap, min_dist_cones);
+
+ % compute representation error
+ err = zeros(Ntest,1);
+ for k = 1 : Ntest
+ [x, Da, support] = omp_cone_single(Ytest(:,k), D1, s, rad_swap, iter_cd, err_thresh);
+ err(k) = norm(Ytest(:,k) - Da*x);
+ end
+ res_err{5}(it,i_param) = mean(err);
+ [~,~,~,roc_auc] = perfcurve(ytest,err,1);
+ res_data{5}(it,i_param) = roc_auc;
+ end
+
+ % Cone-DL
+ err = zeros(Ntest,1);
+ for k = 1 : Ntest
+ [x, Da, support] = omp_cone_single(Ytest(:,k), D2, s, rad, iter_cd, err_thresh);
+ err(k) = norm(Ytest(:,k) - Da*x);
+ end
+ res_err{2}(it,i_param) = mean(err);
+ [~,~,~,roc_auc] = perfcurve(ytest,err,1);
+ res_data{2}(it,i_param) = roc_auc;
+
+ % Cone-DL + swap
+ if swap_rad
+ % swap radii
+ [~, Xs] = compute_rep_cone(D2, Ytrain, s, rad, iter_cd, err_thresh);
+ atom_use = sum(squeeze(Xs)~=0,2);
+ [~,Ip] = sort(atom_use,'descend');
+ [~,Ir] = sort(rad,'descend');
+ rad_swap = zeros(n,1);
+ rad_swap(Ip) = rad(Ir);
+ rad_swap = adjust_radii(D2, rad_swap, min_dist_cones);
+
+ % compute representation error
+ err = zeros(Ntest,1);
+ for k = 1 : Ntest
+ [x, Da, support] = omp_cone_single(Ytest(:,k), D2, s, rad_swap, iter_cd, err_thresh);
+ err(k) = norm(Ytest(:,k) - Da*x);
+ end
+ res_err{4}(it,i_param) = mean(err);
+ [~,~,~,roc_auc] = perfcurve(ytest,err,1);
+ res_data{4}(it,i_param) = roc_auc;
+ end
+ end
+ end
+ end
+ %fprintf("ROC AUC: %f\n", mean(res_roc_auc(:,id)))
+ for im = 1 : n_methods
+ res_methods{im}(id,:) = mean(res_data{im});
+ fprintf("%40s: %4f\n", method_names{im}, max(res_methods{im}(id,:)))
+ end
+ % write in file
+ for im = 1 : n_methods
+ fprintf(fid_res, "%s,%s,", data_file_names{id}, method_names{im});
+ for ip = 1 : n_params
+ fprintf(fid_res, "%f,", res_methods{im}(id,ip));
+ end
+ fprintf(fid_res, "\n");
+ end
+end
+
+fclose(fid_res);