From a66362a111b07231703a096cfb97195a7e560695 Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Fri, 1 Apr 2022 15:04:52 +0200 Subject: [PATCH 01/21] Calculate the cdf of the dose for several distributions --- cProbOpt/apm_calcCdfDose.m | 77 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 77 insertions(+) create mode 100644 cProbOpt/apm_calcCdfDose.m diff --git a/cProbOpt/apm_calcCdfDose.m b/cProbOpt/apm_calcCdfDose.m new file mode 100644 index 0000000..e5ca97d --- /dev/null +++ b/cProbOpt/apm_calcCdfDose.m @@ -0,0 +1,77 @@ +function u = apm_calcCdfDose(d, mu, sigma, cdf_model, is_upper, dmin, dmax) +% Returns CDF(d; par1,par2), where CDF uses the model in cdf_model, +% of parameters par1, par2, set such that the pdf of the dose has the +% specified mean and standard deviation (mu, sigma). +% If is_upper == 1, returns 1-CDF(d; par1,par2) +if nargin<5 + is_upper = 0; +end +% Support of the shifted beta distribution +if nargin<6 + dmin = 0; +end +if nargin<7 + dmax = 1; +end + +switch true + case (strcmpi(cdf_model,'gauss')==1) + if is_upper == 0 + u = normcdf(d, mu, sigma); + else + u = normcdf(d, mu, sigma, 'upper'); + end + + case (strcmpi(cdf_model,'loggauss')==1) + [lnMu, lnVar] = apm_transformMeanCovarianceToLogNormalParameters(mu, sigma^2); + % Use that the lognormal cdf is the same as the normcdf + % evaluated at log(d) + if is_upper == 0 + u = normcdf(log(d), lnMu, lnVar); + else + u = normcdf(log(d), lnMu, lnVar, 'upper'); + end + + case (strcmpi(cdf_model,'beta')==1) + [alpha, beta] = apm_transformMeanVarianceToBetaParameters(mu, sigma); + if is_upper == 0 + u = betacdf(d, alpha, beta); + else + u = betacdf(d, alpha, beta, 'upper'); + end + + case (strcmpi(cdf_model,'shiftedbeta')==1) + [mu_prime,sigma_prime] = apm_transformMeanVarianceShiftedBetaToBetaMeanVariance(mu, sigma, dmin, dmax); + [alpha,beta] = apm_transformMeanVarianceToBetaParameters(mu_prime, sigma_prime); + d_prime = (d - dmin) ./ (dmax - dmin); + u = zeros(1,numel(d_prime)); + if is_upper == 0 + u(d_prime<0) = 0; + u(d_prime>1) = 1; + u((d_prime>=0) & (d_prime<=1)) = betacdf(d_prime((d_prime>=0) & (d_prime<=1)), alpha,beta); + else + u(d_prime<0) = 1; + u(d_prime>1) = 0; + u((d_prime>=0) & (d_prime<=1)) = betacdf(d_prime((d_prime>=0) & (d_prime<=1)), alpha,beta, 'upper'); + end + + + case (strcmpi(cdf_model,'gamma')==1) + [k,theta] = apm_transformMeanVarianceToGammaParameters(mu,sigma); + if is_upper == 0 + u = gamcdf(d, k, theta); + else + u = gamcdf(d, k, theta, 'upper'); + end + + case (strcmpi(cdf_model,'gumbel')==1) + [epsilon,alpha] = apm_transformMeanVarianceToGumbelParameters(mu,sigma); + if is_upper == 0 + u = evcdf(d, epsilon, alpha); + else + u = evcdf(d, epsilon, alpha, 'upper'); + end + otherwise + error(['Invalid copula model']); +end +end From 7b10e3c022f4a632f0d7cf9b5d5d469d8cfb541c Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Fri, 1 Apr 2022 15:08:13 +0200 Subject: [PATCH 02/21] Compute different Copulas --- cProbOpt/apm_calcCopula.m | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 cProbOpt/apm_calcCopula.m diff --git a/cProbOpt/apm_calcCopula.m b/cProbOpt/apm_calcCopula.m new file mode 100644 index 0000000..ce642e2 --- /dev/null +++ b/cProbOpt/apm_calcCopula.m @@ -0,0 +1,36 @@ +function c = apm_calcCopula(u1, u2, par, copula_model) +% Returns Copula(u1,u2; par), using the copula model in the input, +% where 'Copula' is the joint cummulative distribution function of u1 and u2 +% with parameter par. + +switch true + case (strcmpi(copula_model,'gauss') == 1) % Gaussian Copula + r = par; % Correlation coefficient + v1 = norminv(u1); % Inverse Gaussian CDF + v2 = norminv(u2); + c = arrayfun(@(v1,v2,r) bvn(-Inf,v1,-Inf,v2,r), v1, v2, r); + + case (strcmpi(copula_model,'amh') == 1) % Ali-Mikhail-Haq + c = u1 .* u2 ./ (1 - par .* (1-u1) .* (1-u2)); + + case (strcmpi(copula_model, 'gumbel') == 1) % Gumbel + c = copulacdf('Gumbel', [u1(:),u2(:)], par); + c = reshape(c,numel(u1(:,1)), numel(u2(:,2))); + + case (strcmpi(copula_model, 'clayton') == 1) % Clayton + c = copulacdf('Clayton', [u1(:),u2(:)], par); + c = reshape(c,numel(u1(:,1)), numel(u2(:,2))); + + case (strcmpi(copula_model, 'frank') == 1) % Frank + c = copulacdf('frank', [u1(:),u2(:)], par); + c = reshape(c,numel(u1(:,1)), numel(u2(:,2))); + %case (strcmpi(copula_model, 'uniform') == 1) % Uniform perfect correlation + % c = min(u1, u2); + % if par ~= 1 + % disp(['Unapropriate Copula model!']); + % end + otherwise + disp (['Invalid Copula model']); +end + +end From a5d61f3e4050fbd8833f6445c612a63102f52e24 Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 12:15:35 +0200 Subject: [PATCH 03/21] Add the computation of the moments of the probabilistic DVH using a loggauss distribution of the dose or a Copula of user-defined model and marginals --- cProbOpt/apm_DVHprob.m | 203 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 197 insertions(+), 6 deletions(-) diff --git a/cProbOpt/apm_DVHprob.m b/cProbOpt/apm_DVHprob.m index 73c6af9..df8e237 100644 --- a/cProbOpt/apm_DVHprob.m +++ b/cProbOpt/apm_DVHprob.m @@ -1,4 +1,4 @@ -function [expDvh,stdDvh,covDvh] = apm_DVHprob(expDose,covDose,nBins,dMax,method) +function [expDvh,stdDvh,covDvh] = apm_DVHprob(expDose,covDose,nBins,dMax,method,copula_model,copula_marginals,dmin_shiftedbeta,dmax_shiftedbeta) if nargin < 3 || isempty(nBins) nBins = 100; @@ -12,6 +12,25 @@ method = 'int_gauss'; end +% Define parameters of the Copula +if nargin < 6 + copula_model = 'gauss'; +end + +if nargin < 7 + copula_marginals = cell(1,numel(expDose)); + copula_marginals(:) = {'gauss'}; +end + +% Define support of the shifted beta +if nargin < 8 + dmin_shiftedbeta= zeros(1,numel(expDose)); +end + +if nargin < 9 + dmax_shiftedbeta = ones(1, numel(expDose)); +end + dAxis = linspace(0,dMax,nBins); expDvh(1,:) = dAxis; stdDvh(1,:) = dAxis; @@ -35,17 +54,28 @@ reducedVoxNum = numel(expDose); switch method + case 'copula' + [expDvh(2,:),stdDvh(2,:)] = probDVHCopulaExpStd(expDose,covDose,dAxis,copula_model,copula_marginals,dmin_shiftedbeta,dmax_shiftedbeta); + if nargout > 2 + covDvh = calcCovDvhCopula(expDose,covDose,dAxis,expDvh(2,:),copula_model,copula_marginals,dmin_shiftedbeta,dmax_shiftedbeta); + end case 'int_gauss' [expDvh(2,:),stdDvh(2,:)] = probDVHGaussExpStd(expDose,covDose,dAxis); if nargout > 2 covDvh = calcCovDvhGauss(expDose,covDose,dAxis,expDvh(2,:)); end + + case 'int_loggauss' + [expDvh(2,:),stdDvh(2,:)] = probDVHLnGaussExpStd(expDose,covDose,dAxis); + if nargout > 2 + covDvh = calcCovDvhLnGauss(expDose,covDose,dAxis,expDvh(2,:)); + end + case 'int_gauss_cuda' disp('Setting up kernel'); tic - - + %single precision %{ cuKernel = parallel.gpu.CUDAKernel('C:\dev\APM\CUDA\QI\DVH\x64\Release\dvh.ptx','float* , float*, const float*, const float*, const float*, int, int','dvhRawMoments'); @@ -145,6 +175,85 @@ end +function [expDvh,stdDvh] = probDVHCopulaExpStd(expDose,covDose,dAxis,copula_model,copula_marginals,dmin_shiftedbeta,dmax_shiftedbeta) + nBins = numel(dAxis); + biCum = zeros(1,nBins); + diagCum = zeros(1,nBins); + + for i = 1:numel(expDose) + mean_i = expDose(i); + sig_i = sqrt(covDose(i,i)); + if ~isreal(sig_i) || isnan(sig_i) + continue; + end + P1 = apm_calcCdfDose(dAxis, mean_i*ones(numel(1,nBins)), sig_i*ones(numel(1,nBins)), copula_marginals(i), 1,dmin_shiftedbeta(i),dmax_shiftedbeta(i)); + + + if any(isnan(P1)) + warning('CDF evaluates to NaN'); + P1(isnan(P1)) = 0; + end + diagCum = diagCum + P1; + + for l = i+1:numel(expDose) + mean_l = expDose(l); + sig_l = sqrt(covDose(l,l)); + if ~isreal(sig_l) || isnan(sig_l) + continue; + end + r = (covDose(i,l) / (sig_i*sig_l)) * ones(1,nBins); + u_i = apm_calcCdfDose(dAxis, mean_i, sig_i, copula_marginals(i),0,dmin_shiftedbeta(i),dmax_shiftedbeta(i)); + u_l = apm_calcCdfDose(dAxis, mean_l, sig_l, copula_marginals(l),0,dmin_shiftedbeta(l),dmax_shiftedbeta(l)); + P2D = apm_calcCopula(u_i, u_l, r, copula_model) + 1 - u_i - u_l; + + if any(isnan(P2D)) + warning('CDF (2D) evaluates to NaN'); + P2D(isnan(P2D)) = 0; + end + + biCum = biCum + P2D; + end + matRad_progress(i,numel(expDose)); + end + + % biCum has only the terms from l=i+1 to numel(expDose), + % then we need to add the terms l=i (diagCum) and l from 0 to i-1 (biCum). + biCum = diagCum + 2*biCum; + + expDvh = diagCum ./ numel(expDose); + + stdDvh = sqrt( 1/numel(expDose)^2 * (biCum) - expDvh.^2); +end + +function covDvh = calcCovDvhCopula(expDose,covDose,dAxis,expDvh,copula_model,copula_marginals,dmin_shiftedbeta,dmax_shiftedbeta) + nBins = numel(dAxis); + covDvh = zeros(nBins); + + for p = 1:nBins + for q = 1:nBins + s = 0; + for i = 1:numel(expDose) + mean_i = expDose(i); + sig_i = sqrt(covDose(i,i)); + for l = 1:numel(expDose) + mean_l = expDose(l); + sig_l = sqrt(covDose(l,l)); + r = covDose(i,l) / (sig_i*sig_l); + u_i = apm_calcCdfDose(dAxis(p), mean_i, sig_i, copula_marginals(i),0,dmin_shiftedbeta(i),dmax_shiftedbeta(i)); + u_l = apm_calcCdfDose(dAxis(q), mean_l, sig_l, copula_marginals(l),0,dmin_shiftedbeta(l),dmax_shiftedbeta(l)); + s = s + apm_calcCopula(u_i, u_l, r, copula_model) + 1 - u_i - u_l; + end + end + + covDvh(p,q) = s; + end + matRad_progress(p,nBins); + end + covDvh = 1/numel(expDose)^2 * covDvh; + covDvh = covDvh - transpose(expDvh)*expDvh; +end + + function [expDvh,stdDvh] = probDVHGaussExpStd(expDose,covDose,dAxis) nBins = numel(dAxis); biCum = zeros(1,nBins); @@ -175,7 +284,6 @@ r = covDose(i,l) / (sig_i*sig_l); P2D = arrayfun(@(dParam) bvn((dParam - expDose(i)) / sig_i,Inf,(dParam - expDose(l)) / sig_l,Inf,r),dAxis); - %P2D = arrayfun(@(dParam) mexBVNcdf(-[dParam dParam],-[expDose(i) expDose(l)],[covDose(i,i) covDose(i,l); covDose(l,i) covDose(l,l)]),dAxis); if any(isnan(P2D)) warn('CDF (2D) evaluates to NaN'); @@ -186,11 +294,14 @@ end matRad_progress(i,numel(expDose)); end + + % biCum has only the terms from l=i+1 to numel(expDose), + % then we need to add the terms l=i (diagCum) and l from 0 to i-1 (biCum). + biCum = diagCum + 2*biCum; expDvh = diagCum ./ numel(expDose); - stdDvh = 1/numel(expDose)^2 * (diagCum + 2*biCum); - stdDvh = sqrt(stdDvh - expDvh.^2); + stdDvh = sqrt( 1/numel(expDose)^2 * (biCum) - expDvh.^2); end function covDvh = calcCovDvhGauss(expDose,covDose,dAxis,expDvh) @@ -216,3 +327,83 @@ covDvh = 1/numel(expDose)^2 * covDvh; covDvh = covDvh - transpose(expDvh)*expDvh; end + + +function [expDvh,stdDvh] = probDVHLnGaussExpStd(expDose,covDose,dAxis) + nBins = numel(dAxis); + biCum = zeros(1,nBins); + diagCum = zeros(1,nBins); + + %Compute parameters of the lognormal + [lnExpDose, lnCovDose] = apm_transformMeanCovarianceToLogNormalParameters(expDose,covDose); + + for i = 1:numel(lnExpDose) %for each voxel + lnSig_i = sqrt(lnCovDose(i,i)); + + if ~isreal(lnSig_i) || isnan(lnSig_i) + continue; + end + + P1 = 1 - apm_normcdf(log(dAxis),lnExpDose(i)*ones(numel(1,nBins)),lnSig_i*ones(numel(1,nBins))); + + if any(isnan(P1)) + warning('CDF evaluates to NaN'); + P1(isnan(P1)) = 0; + end + diagCum = diagCum + P1; + + for l = i+1:numel(lnExpDose) + lnSig_l = sqrt(lnCovDose(l,l)); + + if ~isreal(lnSig_l) || isnan(lnSig_l) + continue; + end + + r = lnCovDose(i,l) / (lnSig_i*lnSig_l); + P2D = arrayfun(@(dParam) bvn((log(dParam) - lnExpDose(i)) / lnSig_i,Inf,(log(dParam) - lnExpDose(l)) / lnSig_l,Inf,r), dAxis); + + if any(isnan(P2D)) + warn('CDF (2D) evaluates to NaN'); + P2D(isnan(P2D)) = 0; + end + + biCum = biCum + P2D; + end + matRad_progress(i,numel(expDose)); + end + + % biCum has only the terms from l=i+1 to numel(expDose), + % then we need to add the terms l=i (diagCum) and l from 0 to i-1 (biCum). + biCum = diagCum + 2*biCum; + + expDvh = diagCum ./ numel(lnExpDose); + + stdDvh = sqrt( 1/numel(lnExpDose)^2 * (biCum) - expDvh.^2); +end + +function covDvh = calcCovDvhLnGauss(expDose,covDose,dAxis,expDvh) + nBins = numel(dAxis); + covDvh = zeros(nBins); + + %Compute parameters of the lognormal + [lnExpDose, lnCovDose] = apm_transformMeanCovarianceToLogNormalParameters(expDose,covDose); + + for p = 1:nBins + for q = 1:nBins + s = 0; + for i = 1:numel(lnExpDose) + lnSig_i = sqrt(lnCovDose(i,i)); + for l = 1:numel(lnExpDose) + lnSig_l = sqrt(lnCovDose(l,l)); + r = lnCovDose(i,l) / (lnSig_i*lnSig_l); + s = s + bvn((log(dAxis(p)) - lnExpDose(i)) / lnSig_i,Inf,(log(dAxis(q)) - lnExpDose(l)) / lnSig_l,inferiorto,r); + end + end + + covDvh(p,q) = s; + end + matRad_progress(p,nBins); + end + covDvh = 1/numel(lnExpDose)^2 * covDvh; + covDvh = covDvh - transpose(expDvh)*expDvh; +end From af310e196510f6297a26de76ec31b9b6f8365cf6 Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 12:22:51 +0200 Subject: [PATCH 04/21] Fix bug --- .../apm_transformMeanVarianceToBetaParameters.m | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/cProbOpt/transforms/apm_transformMeanVarianceToBetaParameters.m b/cProbOpt/transforms/apm_transformMeanVarianceToBetaParameters.m index 37cbbb1..de32425 100644 --- a/cProbOpt/transforms/apm_transformMeanVarianceToBetaParameters.m +++ b/cProbOpt/transforms/apm_transformMeanVarianceToBetaParameters.m @@ -1,11 +1,8 @@ function [alpha,beta] = apm_transformMeanVarianceToBetaParameters(mu,sigma) -%This function gives alpha and beta parameters for a beta-distributed -%variable, if their shape parameter is known. +% This function gives alpha and beta parameters for a beta-distributed +% variable, if their shape parameter is known. -%alpha = (mu.^2 - mu.^3 - mu.*(sigma.^2))./ sigma.^2; -%beta = (mu - 1) .* (mu.^2 - mu + sigma.^2) ./ sigma.^2; - -alpha = ((1 - mu.^2)./sigma.^2 - 1./mu) .* mu.^2; +alpha = ((mu .* (1 - mu))./ (sigma.^2) - 1) .* mu; beta = alpha.*(1./mu - 1); if any(sigma.^2 > mu .* (1 - mu)) From 1f129042aa3754f36d1a84e6673f8bec6d5ae956 Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 12:27:20 +0200 Subject: [PATCH 05/21] Return paremters k and theta --- .../apm_transformMeanVarianceToGammaParameters.m | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/cProbOpt/transforms/apm_transformMeanVarianceToGammaParameters.m b/cProbOpt/transforms/apm_transformMeanVarianceToGammaParameters.m index 0cac466..c6124a3 100644 --- a/cProbOpt/transforms/apm_transformMeanVarianceToGammaParameters.m +++ b/cProbOpt/transforms/apm_transformMeanVarianceToGammaParameters.m @@ -1,7 +1,7 @@ -function [p,b] = apm_transformMeanVarianceToGammaParameters(mu,sigma) -%UNTITLED3 Summary of this function goes here -% Detailed explanation goes here - p = mu.^2 ./ sigma.^2; - b = mu ./ sigma.^2; +function [k,theta] = apm_transformMeanVarianceToGammaParameters(mu,sigma) +% Method of the moments for the shape and scale parameters of the gamma distribution +% These are the parameters of gampdf + k = mu.^2 ./ sigma.^2; %shape + theta = sigma.^2 ./ mu; %scale end From 5b85e39382a15810eba65bf1bce426f321a859dc Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 12:41:07 +0200 Subject: [PATCH 06/21] Change name --- ...arameters.m => apm_transformMeanStdToBetaParameters.m} | 4 ++-- .../transforms/apm_transformMeanStdToGammaParameters.m | 7 +++++++ .../apm_transformMeanStdToLogNormalParameters.m | 8 ++++++++ .../apm_transformMeanVarianceToGammaParameters.m | 7 ------- .../apm_transformMeanVarianceToLogNormalParameters.m | 8 -------- 5 files changed, 17 insertions(+), 17 deletions(-) rename cProbOpt/transforms/{apm_transformMeanVarianceToBetaParameters.m => apm_transformMeanStdToBetaParameters.m} (62%) create mode 100644 cProbOpt/transforms/apm_transformMeanStdToGammaParameters.m create mode 100644 cProbOpt/transforms/apm_transformMeanStdToLogNormalParameters.m delete mode 100644 cProbOpt/transforms/apm_transformMeanVarianceToGammaParameters.m delete mode 100644 cProbOpt/transforms/apm_transformMeanVarianceToLogNormalParameters.m diff --git a/cProbOpt/transforms/apm_transformMeanVarianceToBetaParameters.m b/cProbOpt/transforms/apm_transformMeanStdToBetaParameters.m similarity index 62% rename from cProbOpt/transforms/apm_transformMeanVarianceToBetaParameters.m rename to cProbOpt/transforms/apm_transformMeanStdToBetaParameters.m index de32425..98d471d 100644 --- a/cProbOpt/transforms/apm_transformMeanVarianceToBetaParameters.m +++ b/cProbOpt/transforms/apm_transformMeanStdToBetaParameters.m @@ -1,5 +1,5 @@ -function [alpha,beta] = apm_transformMeanVarianceToBetaParameters(mu,sigma) -% This function gives alpha and beta parameters for a beta-distributed +function [alpha,beta] = apm_transformMeanStdToBetaParameters(mu,sigma) +% Return the parameters alpha and beta for a beta-distributed % variable, if their shape parameter is known. alpha = ((mu .* (1 - mu))./ (sigma.^2) - 1) .* mu; diff --git a/cProbOpt/transforms/apm_transformMeanStdToGammaParameters.m b/cProbOpt/transforms/apm_transformMeanStdToGammaParameters.m new file mode 100644 index 0000000..b3b47f4 --- /dev/null +++ b/cProbOpt/transforms/apm_transformMeanStdToGammaParameters.m @@ -0,0 +1,7 @@ +function [k,theta] = apm_transformMeanStdToGammaParameters(mu,sigma) +% Return the shape and scale parameters of the gamma distribution +% (these are the input parameters of gampdf) + k = mu.^2 ./ sigma.^2; %shape + theta = sigma.^2 ./ mu; %scale +end + diff --git a/cProbOpt/transforms/apm_transformMeanStdToLogNormalParameters.m b/cProbOpt/transforms/apm_transformMeanStdToLogNormalParameters.m new file mode 100644 index 0000000..fad8660 --- /dev/null +++ b/cProbOpt/transforms/apm_transformMeanStdToLogNormalParameters.m @@ -0,0 +1,8 @@ +function [lnMu,lnSigma] = apm_transformMeanStdToLogNormalParameters(mu, sigma) +% Return the parameters of the lognormal distribution + +lnSigma = sqrt(log(1 + (sigma.^2 ./ mu.^2))); +lnMu = log(mu) - lnSigma.^2 ./ 2; + +end + diff --git a/cProbOpt/transforms/apm_transformMeanVarianceToGammaParameters.m b/cProbOpt/transforms/apm_transformMeanVarianceToGammaParameters.m deleted file mode 100644 index c6124a3..0000000 --- a/cProbOpt/transforms/apm_transformMeanVarianceToGammaParameters.m +++ /dev/null @@ -1,7 +0,0 @@ -function [k,theta] = apm_transformMeanVarianceToGammaParameters(mu,sigma) -% Method of the moments for the shape and scale parameters of the gamma distribution -% These are the parameters of gampdf - k = mu.^2 ./ sigma.^2; %shape - theta = sigma.^2 ./ mu; %scale -end - diff --git a/cProbOpt/transforms/apm_transformMeanVarianceToLogNormalParameters.m b/cProbOpt/transforms/apm_transformMeanVarianceToLogNormalParameters.m deleted file mode 100644 index b079bc1..0000000 --- a/cProbOpt/transforms/apm_transformMeanVarianceToLogNormalParameters.m +++ /dev/null @@ -1,8 +0,0 @@ -function [lnMu,lnSigma] = apm_transformMeanVarianceToLogNormalParameters(mu, sigma) -%UNTITLED4 Summary of this function goes here -% Detailed explanation goes here -lnSigma = sqrt(log(1 + (sigma.^2 ./ mu.^2))); -lnMu = log(mu) - lnSigma.^2 ./ 2;%log(mu.^2 ./ sqrt(sigma.^2 +mu.^2)); - -end - From 7bcf0bdb7f8f9d9b5c8dccfcef8e0cd51a24c5c5 Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 12:42:21 +0200 Subject: [PATCH 07/21] Add function explanation --- ...transformMeanCovarianceToLogNormalParameters.m | 4 ++-- .../apm_transformMeanStdToGumbelParameters.m | 12 ++++++++++++ .../apm_transformMeanStdToShiftedBetaParameters.m | 15 +++++++++++++++ 3 files changed, 29 insertions(+), 2 deletions(-) create mode 100644 cProbOpt/transforms/apm_transformMeanStdToGumbelParameters.m create mode 100644 cProbOpt/transforms/apm_transformMeanStdToShiftedBetaParameters.m diff --git a/cProbOpt/transforms/apm_transformMeanCovarianceToLogNormalParameters.m b/cProbOpt/transforms/apm_transformMeanCovarianceToLogNormalParameters.m index fe101f4..b89e930 100644 --- a/cProbOpt/transforms/apm_transformMeanCovarianceToLogNormalParameters.m +++ b/cProbOpt/transforms/apm_transformMeanCovarianceToLogNormalParameters.m @@ -1,6 +1,6 @@ function [lnMu,lnSigma] = apm_transformMeanCovarianceToLogNormalParameters(mu, Sigma) -%UNTITLED4 Summary of this function goes here -% Detailed explanation goes here +% Return the mean and covariance matrix of a multivariate lognormal distribution + lnSigma = real(log(1 + (Sigma ./ (mu*mu')))); lnSigma(isnan(lnSigma)) = 0; lnMu = log(mu) - diag(lnSigma) ./ 2;%log(mu.^2 ./ sqrt(sigma.^2 +mu.^2)); diff --git a/cProbOpt/transforms/apm_transformMeanStdToGumbelParameters.m b/cProbOpt/transforms/apm_transformMeanStdToGumbelParameters.m new file mode 100644 index 0000000..856a42c --- /dev/null +++ b/cProbOpt/transforms/apm_transformMeanStdToGumbelParameters.m @@ -0,0 +1,12 @@ +function [epsilon,alpha] = apm_transformMeanStdToGumbelParameters(mu, sigma) +% Return the parameters of a Gumbel distribution + +% See Mahdi, Smail & Cenac, Myrtene. (2012). "Estimating Parameters of Gumbel Distribution using the Methods of Moments, probability weighted Moments and maximum likelihood". Revista de Matemática: Teoría y Aplicaciones. 12. 151. 10.15517/rmta.v12i1-2.259. + +J = 1.978; +gamma = 0.577215664901532; +alpha = sqrt(sigma.^2 ./ (J-gamma.^2)); +epsilon = mu - gamma .* alpha; + +end + diff --git a/cProbOpt/transforms/apm_transformMeanStdToShiftedBetaParameters.m b/cProbOpt/transforms/apm_transformMeanStdToShiftedBetaParameters.m new file mode 100644 index 0000000..4b2939a --- /dev/null +++ b/cProbOpt/transforms/apm_transformMeanStdToShiftedBetaParameters.m @@ -0,0 +1,15 @@ +function [alpha,beta] = apm_transformMeanStdToShiftedBetaParameters(mu,sigma,a,b) +% Return the parameters alpha and beta for a shifted beta-distributed variable +% with support in [a,b], if their location and shape parameters are known. + +aux = (a .* b - (a+b) .* mu + mu.^2 + sigma.^2) / (sigma.^2 .* (b-a)); + +alpha = (a - mu) .* aux; +beta = - (b - mu) .* aux; + +if (any(alpha<0) | any(beta<0)) + warning('Moment to shape transformation actually not defined!'); +end + +end + From 4ca8fac76e1acd580192581a40efb28fae3e3f5e Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 13:54:10 +0200 Subject: [PATCH 08/21] No change --- cProbOpt/apm_calcExpDoseInfluenceLateral.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cProbOpt/apm_calcExpDoseInfluenceLateral.m b/cProbOpt/apm_calcExpDoseInfluenceLateral.m index e708f1e..e57c011 100644 --- a/cProbOpt/apm_calcExpDoseInfluenceLateral.m +++ b/cProbOpt/apm_calcExpDoseInfluenceLateral.m @@ -4,7 +4,7 @@ sigma = [spots(:).sigma]; mu = [spots(:).mu]; -sigmaAddSqr = sigma.^2 + diag(ucm.covSys)' + diag(ucm.covRand)'; +sigmaAddSqr = sigma.^2 + diag(ucm.covSys)' + diag(ucm.covRand)'; edij = 1./sqrt( 2*pi*ones(nVox,1)*sigmaAddSqr) .* exp( - ( bsxfun(@minus,ones(nVox,1)*mu,x).^2 ) ./ ( 2*ones(nVox,1)*sigmaAddSqr) ); From a815484440f568c3e8d26ff734a76a2590f442c0 Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 13:56:28 +0200 Subject: [PATCH 09/21] Write as a module to be reachable from any other module --- cProbOpt/apm_calcDose.m | 4 ++++ cProbOpt/apm_calcExpDose.m | 4 ++++ 2 files changed, 8 insertions(+) create mode 100644 cProbOpt/apm_calcDose.m create mode 100644 cProbOpt/apm_calcExpDose.m diff --git a/cProbOpt/apm_calcDose.m b/cProbOpt/apm_calcDose.m new file mode 100644 index 0000000..9fe79cc --- /dev/null +++ b/cProbOpt/apm_calcDose.m @@ -0,0 +1,4 @@ +function d = apm_calcDose(dij, w) +% Compute nominal dose +d = dij*w; +end \ No newline at end of file diff --git a/cProbOpt/apm_calcExpDose.m b/cProbOpt/apm_calcExpDose.m new file mode 100644 index 0000000..339e09e --- /dev/null +++ b/cProbOpt/apm_calcExpDose.m @@ -0,0 +1,4 @@ +function d = apm_calcExpDose(edij, w) +% Compute expected dose +d = edij*w; +end \ No newline at end of file From eb1009bdf94e3ffb4d356f3e0a3b8ddeacf373f9 Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 14:26:03 +0200 Subject: [PATCH 10/21] Change names of transform calls --- cProbOpt/apm_calcCdfDose.m | 10 +++++----- cProbOpt/apm_expPowerX.m | 2 +- cProbOpt/apm_getPDVHfromProbDVH.m | 4 ++-- cProbOpt/apm_plotProbDVH.m | 2 +- 4 files changed, 9 insertions(+), 9 deletions(-) diff --git a/cProbOpt/apm_calcCdfDose.m b/cProbOpt/apm_calcCdfDose.m index e5ca97d..fce45e8 100644 --- a/cProbOpt/apm_calcCdfDose.m +++ b/cProbOpt/apm_calcCdfDose.m @@ -33,7 +33,7 @@ end case (strcmpi(cdf_model,'beta')==1) - [alpha, beta] = apm_transformMeanVarianceToBetaParameters(mu, sigma); + [alpha, beta] = apm_transformMeanStdToBetaParameters(mu, sigma); if is_upper == 0 u = betacdf(d, alpha, beta); else @@ -41,8 +41,8 @@ end case (strcmpi(cdf_model,'shiftedbeta')==1) - [mu_prime,sigma_prime] = apm_transformMeanVarianceShiftedBetaToBetaMeanVariance(mu, sigma, dmin, dmax); - [alpha,beta] = apm_transformMeanVarianceToBetaParameters(mu_prime, sigma_prime); + [mu_prime,sigma_prime] = apm_transformMeanStdShiftedBetaToBetaMeanStd(mu, sigma, dmin, dmax); + [alpha,beta] = apm_transformMeanStdToBetaParameters(mu_prime, sigma_prime); d_prime = (d - dmin) ./ (dmax - dmin); u = zeros(1,numel(d_prime)); if is_upper == 0 @@ -57,7 +57,7 @@ case (strcmpi(cdf_model,'gamma')==1) - [k,theta] = apm_transformMeanVarianceToGammaParameters(mu,sigma); + [k,theta] = apm_transformMeanStdToGammaParameters(mu,sigma); if is_upper == 0 u = gamcdf(d, k, theta); else @@ -65,7 +65,7 @@ end case (strcmpi(cdf_model,'gumbel')==1) - [epsilon,alpha] = apm_transformMeanVarianceToGumbelParameters(mu,sigma); + [epsilon,alpha] = apm_transformMeanStdToGumbelParameters(mu,sigma); if is_upper == 0 u = evcdf(d, epsilon, alpha); else diff --git a/cProbOpt/apm_expPowerX.m b/cProbOpt/apm_expPowerX.m index 1c6ea15..e269edf 100644 --- a/cProbOpt/apm_expPowerX.m +++ b/cProbOpt/apm_expPowerX.m @@ -23,7 +23,7 @@ %Transform mu and sig to lognormal %mu_log = log(mu./sqrt(1 + sig.^2 ./ mu.^2)); %sig_log = sqrt(log(1 + sig.^2 ./ mu.^2)); - [mu_log,sig_log] = apm_transformMeanVarianceToLogNormalParameters(mu,sig); + [mu_log,sig_log] = apm_transformMeanStdToLogNormalParameters(mu,sig); %Enjoy the easy lognormal form expPowerX = exp(k * mu_log + k^2 * sig_log.^2 / 2); diff --git a/cProbOpt/apm_getPDVHfromProbDVH.m b/cProbOpt/apm_getPDVHfromProbDVH.m index 1dc7ff1..0329cd4 100644 --- a/cProbOpt/apm_getPDVHfromProbDVH.m +++ b/cProbOpt/apm_getPDVHfromProbDVH.m @@ -24,7 +24,7 @@ case 'normal' pdvh = arrayfun(@(mu,sigma) getNormalQuantile(mu,sigma,p,trunc),expDvh,stdDvh); case 'beta' - [a,b] = apm_transformMeanVarianceToBetaParameters(expDvh,stdDvh); + [a,b] = apm_transformMeanStdToBetaParameters(expDvh,stdDvh); pdvh = arrayfun(@(a,b) getBetaQuantile(a,b,p),a,b); otherwise error(['Function not implemented for Distributions of type ' distName ]); @@ -44,4 +44,4 @@ function pVal = getBetaQuantile(a,b,p) dist = makedist('beta','a',a,'b',b); pVal = dist.icdf(p); -end \ No newline at end of file +end diff --git a/cProbOpt/apm_plotProbDVH.m b/cProbOpt/apm_plotProbDVH.m index 1fac4a7..de6d8f0 100644 --- a/cProbOpt/apm_plotProbDVH.m +++ b/cProbOpt/apm_plotProbDVH.m @@ -53,7 +53,7 @@ function apm_plotProbDVH(h,expDvh,stdDvh,color,styles,mode) %parametric transform - [a,b] = apm_transformMeanVarianceToBetaParameters(mu,sig); + [a,b] = apm_transformMeanStdToBetaParameters(mu,sig); aMat = repmat(a,numel(dvSpace),1); bMat = repmat(b,numel(dvSpace),1); From a48d85e5b3a5046bc80065b708ca2b74d48e94e8 Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 14:27:31 +0200 Subject: [PATCH 11/21] Get support of a shifted beta distribution of the dose in a voxel --- cProbOpt/transforms/apm_getShiftedBetaRange.m | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) create mode 100644 cProbOpt/transforms/apm_getShiftedBetaRange.m diff --git a/cProbOpt/transforms/apm_getShiftedBetaRange.m b/cProbOpt/transforms/apm_getShiftedBetaRange.m new file mode 100644 index 0000000..ec61093 --- /dev/null +++ b/cProbOpt/transforms/apm_getShiftedBetaRange.m @@ -0,0 +1,17 @@ +function [dmin, dmax] = apm_getShiftedBetaRange(mu, sigma, x, vois) +% Get the support of the shifted beta function (heuristic) + +% For data +dmin = max(mu - 3.*sigma, 0); % No negative dose +dmax = mu + 3.*sigma; + +if nargin>2 % For simulations with perfect correlation + dmin = max(mu - 1.25.*sigma, 0); % No negative dose + dmax = mu + 1.25.*sigma; + % Set the range to [0, 1+] outside the ill area + a buffer zone + dmin(find(xvois(1).xU*0.8)) = 0; + dmax(find(xvois(1).xU*0.8)) = max(1, mu(find(x>vois(1).xU*0.8)) + 1.25 .* sigma(find(x>vois(1).xU*0.8))); +end +end From b7fdfafd68ec0f1a2aa6810fce73beef4d5c2f9f Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 14:28:34 +0200 Subject: [PATCH 12/21] Get the mean and std of a beta distribution from that of a shifted beta --- .../apm_transformMeanStdShiftedBetaToBetaMeanStd.m | 9 +++++++++ 1 file changed, 9 insertions(+) create mode 100644 cProbOpt/transforms/apm_transformMeanStdShiftedBetaToBetaMeanStd.m diff --git a/cProbOpt/transforms/apm_transformMeanStdShiftedBetaToBetaMeanStd.m b/cProbOpt/transforms/apm_transformMeanStdShiftedBetaToBetaMeanStd.m new file mode 100644 index 0000000..6f0e5a1 --- /dev/null +++ b/cProbOpt/transforms/apm_transformMeanStdShiftedBetaToBetaMeanStd.m @@ -0,0 +1,9 @@ +function [mu,sigma] = apm_transformMeanStdToShiftedBetaMeanStd(mu_shifted,sigma_shifted,a,b) +% Return the mean and sigma of the beta distribution, which has support in [0,1], +% that is equivalent to a shifted beta which has support in [a,b]. + +mu = (mu_shifted - a) ./ (b - a); +sigma = (sigma_shifted) ./ (b - a); + +end + From eae6abca646ba3c6c3d1daf8ab533a5013c61bbd Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 14:30:12 +0200 Subject: [PATCH 13/21] Get random multivariate gaussian samples --- cProbOpt/mb_mgd.m | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) create mode 100644 cProbOpt/mb_mgd.m diff --git a/cProbOpt/mb_mgd.m b/cProbOpt/mb_mgd.m new file mode 100644 index 0000000..8d12b3e --- /dev/null +++ b/cProbOpt/mb_mgd.m @@ -0,0 +1,32 @@ +function r = mb_mgd(cases,mu,sigma) + +[m1 n1] = size(mu); +c = max([m1 n1]); +if m1 .* n1 ~= c + error('Mu must be a vector.'); +end + +[m n] = size(sigma); +if m ~= n + error('Sigma must be square'); +end + +if m ~= c + error('The length of mu must equal the number of rows in sigma.'); +end + +[T p] = chol(sigma); +if p ~= 0 + error('Sigma must be a positive definite matrix.'); +end + + +if m1 == c + mu = mu'; +end + +mu = mu(ones(cases,1),:); + +r = randn(cases,c) * T + mu; + +end From 30c204447acbf5100f736a682c97788b6fee66f4 Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 14:31:48 +0200 Subject: [PATCH 14/21] Change names of transform calls --- cProbOpt/apm_probDvhMaxQuantileConstraint.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cProbOpt/apm_probDvhMaxQuantileConstraint.m b/cProbOpt/apm_probDvhMaxQuantileConstraint.m index 6f8cf88..a7917d9 100644 --- a/cProbOpt/apm_probDvhMaxQuantileConstraint.m +++ b/cProbOpt/apm_probDvhMaxQuantileConstraint.m @@ -40,7 +40,7 @@ case 'normal' finv = dvExp + dvStd*sqrt(2)*erfinv(2*obj.p - 1); %case 'gamma' - % [p,b] = apm_transformMeanVarianceToGammaParameters(dvExp,dvStd); + % [p,b] = apm_transformMeanStdToGammaParameters(dvExp,dvStd); end c = finv-obj.maxV; From b95090b15ecbd29edb62218002eee25fc0de81d3 Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 16:05:03 +0200 Subject: [PATCH 15/21] Modularize --- cProbOpt/apm_sampleTreatments.m | 55 +++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) create mode 100644 cProbOpt/apm_sampleTreatments.m diff --git a/cProbOpt/apm_sampleTreatments.m b/cProbOpt/apm_sampleTreatments.m new file mode 100644 index 0000000..950e5cd --- /dev/null +++ b/cProbOpt/apm_sampleTreatments.m @@ -0,0 +1,55 @@ +function doseSampleS = apm_sampleTreatments(nS_S, nFrac, mu, sigma, ucm, xSamp, w, showWaitbar, showScatterPlots) +% Generates nS_S treatments of nFrac fractions each, with systematic and +% random uncertainties in the positions of the spots of weight w, +% evaluated at the voxels in positions xSamp + +if nargin < 8 + showWaitbar = true; +end + +if nargin < 9 + showScatterPlots = false; +end + +% Sample the systematics in the positions of the spots (same for all fractions) +S_S = mb_mgd(nS_S,mu,ucm.covSys)'; + +% Create the matrix that will store the sampled treatments +doseSampleS = zeros(numel(xSamp),nS_S); + +for i = 1:nS_S + % Sample the random fluctuations in the positions of the spots (that are + % already shifted from the systematics) + S_R = mb_mgd(nFrac,S_S(:,i)',ucm.covRand)'; + + doseSampleF = zeros(numel(xSamp), nFrac); + + for f = 1:nFrac + shifted_Spots = struct('mu',num2cell(S_R(:,f)),'sigma', num2cell(sigma)); + + % Calculate dose at voxels in xSamp + doseSampleF_ij = apm_calcDoseInfluenceLateral(xSamp,shifted_Spots); + doseSampleF(:,f) = apm_calcDose(doseSampleF_ij,w); + end + doseSampleS(:,i) = sum(doseSampleF,2)/nFrac; + + % Plot the accumulated treatment dose sample + if (showScatterPlots) + plot(xSamp,doseSampleS(:,i),'.','Color',[0.6 0.6 1],'MarkerSize',2); + drawnow; + end + + %plot(xSamp,std(yR),'.','Color',[1 0.6 1],'MarkerSize',2) + if (showWaitbar) + if (i == 1) + h = waitbar(i/nS_S,['Sample Treatment ' num2str(i) ' of ' num2str(nS_S)]); + else + waitbar(i/nS_S,h,['Sample Treatment ' num2str(i) ' of ' num2str(nS_S)]); + end + else + matRad_progress(i,nS_S); + end +end +if (showWaitbar) + delete(h); +end From 7c2f528c84c9f25b5de85b68e43b092fe749b45c Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 16:47:47 +0200 Subject: [PATCH 16/21] Modularize the sample validation and add optional plots --- cProbOpt/apm_runSampleValidation.m | 135 +++++++++++++++++++++++++++++ 1 file changed, 135 insertions(+) create mode 100644 cProbOpt/apm_runSampleValidation.m diff --git a/cProbOpt/apm_runSampleValidation.m b/cProbOpt/apm_runSampleValidation.m new file mode 100644 index 0000000..7fb8021 --- /dev/null +++ b/cProbOpt/apm_runSampleValidation.m @@ -0,0 +1,135 @@ +function [doseSampleS, nS_S] = apm_runSampleValidation(w, axDVH, expDose, stdDose, nSamplesTotal, nFrac, mu, sigma, ucm, xLowRes, ixLowRes, xStar, axProfile, t, probDvhMethod, nDvhBins, vois, showWaitbar, showScatterPlots, showAllSampDVHs, showAvgSampDVHs) + % Simulate a number of treatments with fractionation. + % In the profile plot, plot the mean dose +- sdev of the mean, + % and stdev +- stdev of the stdev, and for a specific xStar, + % plot the distribution of the dose. + % Separately, plot the nominal DVHs of each sample. + + if nargin<18 + showWaitbar=true; + end + if nargin<19 + showScatterPlots=false; + end + % Compute number of treatments + nS_S = floor(nSamplesTotal/nFrac); + + % Sample the treatments + doseSampleS = apm_sampleTreatments(nS_S, nFrac, mu, sigma, ucm, xLowRes, w, showWaitbar, showScatterPlots); + doseSampleS_xStar = apm_sampleTreatments(nS_S, nFrac, mu, sigma, ucm, xStar, w, showWaitbar, showScatterPlots); + + % Compute statistics + doseSampleS_mean = mean(doseSampleS,2); + doseSampleS_std = std(doseSampleS,[],2); + %doseSampleS_cov = cov(doseSampleS'); + doseSampleS_mean_stdErr = doseSampleS_std / sqrt(nS_S); + doseSampleS_std_stdErr = doseSampleS_std * apm_calcSampledStdAccRel(nS_S); + + doseSampleS_xStar_mean = mean(doseSampleS_xStar,2); + doseSampleS_xStar_std = std(doseSampleS_xStar,[],2); + + % Plot mean and stdev dose as a function of x + errorbar(axProfile,xLowRes,doseSampleS_mean,doseSampleS_mean_stdErr,'rs','MarkerSize',2); + errorbar(axProfile,xLowRes,doseSampleS_std,doseSampleS_std_stdErr,'ms','MarkerSize',2); + %plot(xLowRes,kurtosis(yS),'y*','MarkerSize',10); + %plot(xLowRes,skewness(yS),'k*','MarkerSize',10); + + + % Plot normalized histogram of dose at xStar + axProfile2 = axes(t); + [counts, centers] = hist(doseSampleS_xStar,floor(sqrt(nS_S))); + binWidths = diff(centers); + counts = (counts/binWidths(1))/nS_S; + % Adjust bar heights to be contained within plot limits + %if ((max(counts)+xStar)>max(vois(1).xU, vois(2).xU)) + % counts = counts*(max(vois(1).xU, vois(2).xU))*0.75/(max(counts)+xStar); + %end + %counts = counts*((max(vois(1).xU, vois(2).xU) - min(vois(1).xL, vois(2).xL)))*0.15/(max(counts)); + barh(axProfile2,centers, counts,'EdgeColor','none','BarWidth',1,'FaceColor',.8*[1 1 1]);%counts+xStar, 'BaseValue', xStar); + plot(axProfile,[xStar xStar],[0 1.2],'--','Color','k','MarkerFaceColor','k'); + + % Plot modelled distribution of dose at xStar + axProfile3 = axes(t); + maxPdf = apm_plotModelledDosePdf(axProfile3, probDvhMethod, doseSampleS_xStar_mean, doseSampleS_xStar_std); + + % Make the histogram and modelled distribution occupy the left 15% of the plot + xlim(axProfile3, [0, max(real(max(counts)/0.15), real(maxPdf/0.15))]); + axProfile2.XAxisLocation = 'top'; + axProfile3.XAxisLocation = 'top'; + xlabel(axProfile3, 'rel. counts'); + axProfile2.Color = 'none'; + axProfile3.Color = 'none'; + axProfile.Box = 'off'; + axProfile2.Box = 'off'; + axProfile3.Box = 'off'; + linkaxes([axProfile axProfile2 axProfile3], 'y'); + linkaxes([axProfile2 axProfile3], 'x'); + + % Compute nominal DVHs for the samples + disp(['Computing nominal DVHs for samples']); + sampDVH_x = zeros(1, nDvhBins); + sampDVH_y_v1 = zeros(1, nDvhBins); + sampDVH_y_v2 = zeros(1, nDvhBins); + for i = 1:nS_S + for v = 1:numel(vois) + % Get the indices of the voxels in ixLowRes that correspond to + % those of vois(v) + isxLowResInVoi = ismember(ixLowRes,find(vois(v).ix==1)); + idxs = find(isxLowResInVoi==1); + % Compute the nominal DVHs + sampDVH_iv = apm_DVH(doseSampleS(idxs,i),nDvhBins,1.1); + % Get the values of the dose (same for all vois and samples) + sampDVH_x = sampDVH_iv(1,:); + if v==1 + sampDVH_y_v1 = vertcat(sampDVH_y_v1, sampDVH_iv(2,:)); + end + if v==2 + sampDVH_y_v2 = vertcat(sampDVH_y_v2, sampDVH_iv(2,:)); + end + %plot(axSamples,sampDVH_iv(1,:),sampDVH_iv(2,:),'LineWidth',1); + end + end + sampDVH_y_v1 = sampDVH_y_v1(2:(nS_S+1),:); + sampDVH_y_v2 = sampDVH_y_v2(2:(nS_S+1),:); + + % Compute statistics on the samples + sampDVH_v1_mean = mean(sampDVH_y_v1,1); + sampDVH_v1_std = std(sampDVH_y_v1,[],1); + sampDVH_v2_mean = mean(sampDVH_y_v2,1); + sampDVH_v2_std = std(sampDVH_y_v2,[],1); + + % Plot the individual sampled DVHs + if showAllSampDVHs + figure; + axSamples = axes(); + hold(axSamples, 'on'); + for i = 1:nS_S + plot(axSamples,sampDVH_y_v1(i,:),sampDVH_y_v2(i,:),'LineWidth',1); + end + box(axSamples,'on'); + grid(axSamples,'on'); + ylim(axSamples,[0 1]); + xlim(axSamples,[0 1.1]); + xlabel(axSamples,'rel. dose'); + ylabel(axSamples,'rel. volume'); + title('Sample DHVs'); + end + + % Plot the average sampled DVHs + if showAvgSampDVHs + plot(axDVH, sampDVH_x, sampDVH_v1_mean, 'ks', 'MarkerSize',0.5); + plot(axDVH, sampDVH_x, sampDVH_v2_mean, 'ks', 'MarkerSize',0.5); + fill(axDVH, [sampDVH_x'; flipud(sampDVH_x')], [(sampDVH_v1_mean - sampDVH_v1_std)'; flipud((sampDVH_v1_mean + sampDVH_v1_std)')],'k','FaceAlpha',0.1,'LineStyle','none'); + fill(axDVH, [sampDVH_x'; flipud(sampDVH_x')], [(sampDVH_v2_mean - sampDVH_v2_std)'; flipud((sampDVH_v2_mean + sampDVH_v2_std)')],'k','FaceAlpha',0.1,'LineStyle','none'); + end + + %Evaluate accuracy of sampled vs analytical dose(x) + dose_mean_atSamples = expDose(ixLowRes); + dose_std_atSamples = stdDose(ixLowRes); + + rmse_mean = sqrt(mean((dose_mean_atSamples - doseSampleS_mean).^2)); + rmse_std = sqrt(mean((dose_std_atSamples - doseSampleS_std).^2)); + + disp(['Mean RMSE: ' num2str(rmse_mean)]); + disp(['Std RMSE: ' num2str(rmse_std)]); +end \ No newline at end of file From d6d848724db2c998f8b091f6eedf89a7b4808539 Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 17:04:33 +0200 Subject: [PATCH 17/21] Modularize and extend functionality for other dose distribution assumptions --- cProbOpt/apm_Lateral.m | 560 ++++++++++++++++++----------------------- 1 file changed, 248 insertions(+), 312 deletions(-) diff --git a/cProbOpt/apm_Lateral.m b/cProbOpt/apm_Lateral.m index b81bd0b..8faa602 100644 --- a/cProbOpt/apm_Lateral.m +++ b/cProbOpt/apm_Lateral.m @@ -1,56 +1,64 @@ -% clear environment +% Clear environment clear %clf %close all; set(groot,'defaultlegendinterpreter','latex'); set(groot,'defaulttextinterpreter','latex'); set(groot,'defaultaxesfontsize',14); -% set state of random engine +% Set state of random engine stream = RandStream('mt19937ar','Seed',datenum(datetime)); %stream = RandStream('mt19937ar','Seed',26689); RandStream.setGlobalStream(stream); plotColors = apm_plotColors; -% some parameters we need -fullCovariance = true; -sampleValidation = false; +% Anatomy +nVox = 100; %100 +res = 1; %[mm] +relativeTargetSize = 0.4; +relativeOARSize = 0.2; +% Steering parameters +fullCovariance = true; +sampleValidation = true; +probabilisticOptimization = true; separateDVHplots = false; plotBodyDVH = false; dvhPlotMode = 'tripleband'; %dvhPlotMode = 'beta'; - -%Anatomy -nVox = 100; -res = 1; %[mm] -relativeTargetSize = 0.4; -relativeOARSize = 0.2; - -%Irradiation geometry -spotDistance = 3; %Spot Spacing +probDvhMethod = 'int_gauss'; %'int_gauss', 'int_gauss_cuda', 'int_loggauss' or 'copula' +% Define copula model and marginals if 'copula' is used as probDvhMethod +copula_model = 'amh'; % 'amh' or 'gauss' +copula_marginals = cell(1,nVox); +%copula_marginals(1,1:floor(0.1*nVox)) = {'loggauss'}; +%copula_marginals(1,(floor(0.1*nVox)+1):floor(0.8*nVox)) = {'shiftedbeta'}; +%copula_marginals(1,(floor(0.8*nVox)+1):nVox) = {'loggauss'}; +copula_marginals(:) = {'gauss'}; %'gauss', 'loggauss', 'beta', 'shiftedbeta', 'gamma', or 'gumbel' + +% Irradiation geometry +spotDistance = 3; % Spot Spacing spotWidth = 5; nFrac = 1; % Number of Fractions +% DVHs nDvhBins = 500; - - -%Uncertainty +% Uncertainty sigmaS = 1; sigmaR = 2; correlationModel = 'perfect'; +% Miscallaneous %latCutOff = 3; - tikzFolder = 'tikz/'; %Local tikz folder -%Create Anatomy and Irradiation geometry +% Create Anatomy and Irradiation geometry [x,vois] = apm_createAnatomy1D(nVox,res,relativeTargetSize,relativeOARSize); -figure; -axProfile = axes();%axes('Position',get(axHist,'Position'),'ActivePositionProperty','Position'); +t=tiledlayout(1,1); +axProfile = axes(t); +%axProfile = axes('Position',get(axHist,'Position'),'ActivePositionProperty','Position'); apm_anatomyPlot(axProfile,x,vois); spots = apm_setLateralBeamletGrid(x,vois,2.5,spotDistance,3*sqrt(sigmaR^2 + sigmaS^2)); @@ -63,34 +71,34 @@ plotLim1 = x(1); plotLim2 = x(end); -%sample location +% Parameters for sample validation +lowResFac = 1; % Use one every lowResFac voxels +nLowResVox = nVox/lowResFac; +if (floor(nLowResVox) ~= nLowResVox) + error(['ERROR: Choose a lowResFac that is an exact divider of nVox']); +end +ixLowRes = linspace(1, nLowResVox, nLowResVox) * lowResFac; % Indices of xLowRes in x +xLowRes = x(ixLowRes); +nSamplesTotal = 100; +samplingMethod = 'fractions'; %'independent' +showWaitbar = false; +showScatterPlots = false; +showAllSampDVHs = true; +showAvgSampDVHs = true; +saveSampToFile = true; + +% Sample location xStar = mean([xLim1 xLim2])+.3*rand*(xLim2-xLim1); -%Maybe find multiple ways of setting sigmas, e.g. all equal or varying +% Maybe find multiple ways of setting sigmas, e.g. all equal or varying %sigma = 2.5*ones(1,nSpots); - mu = [spots(:).mu]'; sigma = [spots(:).sigma]'; - wStart = ones(nSpots,1);%4*rand(n,1)+2; - - -nSamplesTotal = 100*nFrac; - -samplingMethod = 'fractions'; %'fractions'; %'independent' -showWaitbar = false; -showScatterPlots = false; - -Kn = @(n) sqrt((n-1)/2).*exp(gammaln((n-1)/2) - gammaln(n/2)); -Vn = @(n) 2*((n-1)/2 - exp(2*gammaln(n/2) - 2*gammaln((n-1)/2))); -sampledStdAccRel = @(n) Kn(n).*sqrt(Vn(n))./sqrt(n-1); - multigauss = @(x,mu,Sigma) 1/sqrt((2*pi)^numel(x) * det(Sigma)) * exp(-0.5*(x - mu)' / Sigma * (x-mu)); - - modeStr = ''; if ~strcmp(correlationModel,'perfect') modeStr = ['_' correlationModel]; @@ -101,19 +109,16 @@ %% Compute probabilistic dose influence -calcDose = @(dij,w) dij*w; -calcExpDose = @(edij,w) edij*w; - -%calculate nominal dose influence +% Calculate nominal dose influence dose_ij = apm_calcDoseInfluenceLateral(x,spots); -%calculate expected dose influence +% Calculate expected dose influence expDose_ij = apm_calcExpDoseInfluenceLateral(x,spots,ucm); -%Calculate covariance Influence +% Calculate covariance influence [covInfluenceSys,covInfluenceRand] = apm_calcCovarianceInfluenceLateral(x,spots,ucm,expDose_ij,true); -%Deduct variance influence +% Deduce variance influence varInfluenceSys = zeros(nVox,nSpots,nSpots); varInfluenceRand = zeros(nVox,nSpots,nSpots); for i=1:nVox @@ -125,14 +130,12 @@ covInfluence = 1/nFrac * covInfluenceRand + covInfluenceSys; varInfluence = 1/nFrac * varInfluenceRand + varInfluenceSys; - - - %dij.physicalDose = dose_ij; %dij.physicalExpDose = expDose_ij; %dij.physicalCovDose = covInfluence; %dij.physicalVarDose = varInfluence; + %% Nominal optimization %Create cell arrays for objectives @@ -149,7 +152,6 @@ %constr = 'DVHmax'; %constr = ''; - %Set a default constraint set via helper function vois = apm_setDefaultObjectivesAndConstraints(vois,obj,constr); @@ -168,21 +170,32 @@ [w,fVal] = fmincon(fungrad,wStart,[],[],[],[],zeros(nSpots,1),Inf*ones(nSpots,1),nonlcon,options); - %Map influence -dose = calcDose(dose_ij,w); - -expDose = calcExpDose(expDose_ij,w); - +dose = apm_calcDose(dose_ij,w); +expDose = apm_calcExpDose(expDose_ij,w); covDose = apm_calcCovDose(covInfluence,w); varDose = apm_calcVarDose(varInfluence,w); stdDose = sqrt(varDose); -%calculate DVHs +% Get support of shifted beta distrib. in case they are used as copula marginals +[dmin_shiftedbeta,dmax_shiftedbeta] = apm_getShiftedBetaRange(expDose,stdDose,x,vois); + +% Calculate DVHs for v = 1:numel(vois) disp(['Computing nominal and probabilistic DVHs for VOI ' vois(v).name '...']); vois(v).nomDVH = apm_DVH(dose(vois(v).ix),nDvhBins,1.1); - [vois(v).expDVH,vois(v).stdDVH] = apm_DVHprob(expDose(vois(v).ix),covDose(vois(v).ix,vois(v).ix),nDvhBins,1.1,'int_gauss_cuda'); + expDose(vois(v).ix) + covDose(vois(v).ix) + [vois(v).expDVH,vois(v).stdDVH] = apm_DVHprob(expDose(vois(v).ix),covDose(vois(v).ix,vois(v).ix),nDvhBins,1.1,probDvhMethod,copula_model,copula_marginals(vois(v).ix),dmin_shiftedbeta(vois(v).ix),dmax_shiftedbeta(vois(v).ix)); + %if (v==1) + % figure; + % ax=axes(); + % hold(ax, 'on'); + % plot(ax,linspace(0,1.1,nDvhBins),vois(v).nomDVH(2,:),'Color',plotColors.nomDose,'LineWidth',2); + % plot(ax,linspace(0,1.1,nDvhBins),vois(v).expDVH(2,:),'Color',plotColors.expDose,'LineWidth',2); + % plot(ax,linspace(0,1.1,nDvhBins),vois(v).stdDVH(2,:),'Color',plotColors.stdDose,'LineWidth',2); + % legend({'Nominal', 'Expected', 'Stdev'},'Location','southwest'); + %end end %% Plot result @@ -190,11 +203,12 @@ figure; axProfile = axes();%axes('Position',get(axHist,'Position'),'ActivePositionProperty','Position'); apm_anatomyPlot(axProfile,x,vois); + title('Dose Profile'); end apm_profilePlot(axProfile,x,dose,expDose,stdDose,'--'); -%DVH +% DVH figure; axDVH = axes; hold(axDVH,'on'); @@ -202,7 +216,6 @@ for v = 1:numel(vois) %dvhs{v} = apm_DVH(dose(vois(v).ix),100,1.1); - % %dvhsProb{v} = apm_DVH(doseProb(vois(v).ix),100,1.1); %[expDvhs{v},stdDvhs{v}] = apm_DVHprob(expDose(vois(v).ix),covDose(vois(v).ix,vois(v).ix),100,1.1,'int_gauss_cuda'); if ~strcmp(vois(v).type,'BODY') || plotBodyDVH @@ -216,6 +229,7 @@ xlim(axDVH,[0 1.1]); xlabel(axDVH,'rel. dose'); ylabel(axDVH,'rel. volume'); +title('Conventional optimization'); apm_plotObjConstrInDVH(axDVH,vois,false,plotBodyDVH); @@ -283,279 +297,201 @@ % plot([vois(2).probObjFunc{1}.dMax vois(2).probObjFunc{1}.dMax],[0 1],'--k<','MarkerFaceColor','k'); % end - - - - %% Simulate fractions for each treatment scenario if sampleValidation - % 1st plot histogram around one point in backgorund - % look at one special point and show the distribution within this graph - % by simulating a number of treatments with fractionation - nS_S = floor(nSamplesTotal/nFrac); - muS = mb_mgd(nS_S,mu,ucm.covSys)'; - %yStarR = ones(nS_R,1); - stdMeanR = ones(nS_S,1); - doseStarS = ones(nS_S,1); - for s=1:nS_S - muR = mb_mgd(nFrac,muS(:,s)',ucm.covRand)'; - tmpFracDoses = (1./ sqrt( 2*pi*ones(nFrac,1)*sigma.^2) .* exp( - ( bsxfun(@minus,xStar,muR').^2 ) ./ ( 2*ones(nFrac,1)*sigma.^2) ) )*w; - doseStarS(s) = sum(tmpFracDoses)/nFrac; + [doseSampleS_convOpt, nS_S_convOpt] = apm_runSampleValidation(w, axDVH, expDose, stdDose, nSamplesTotal, nFrac, mu, sigma, ucm, xLowRes, ixLowRes, xStar, axProfile, t, probDvhMethod, nDvhBins, vois, showWaitbar, showScatterPlots, showAllSampDVHs, showAvgSampDVHs); + if saveSampToFile + filename = 'SampledDosesConvOpt.mat'; + save(filename, 'doseSampleS_convOpt', 'nS_S_convOpt', 'xLowRes', 'x', 'dose', 'expDose', 'stdDose', 'vois'); end - [histX,histY]=hist(doseStarS,50); - % normalize histogram - binWidths = diff(histY); - histX = (histX/binWidths(1))/nS_S; - histAxisFac = 5; - axis([0 max(histX)*histAxisFac 0 yLimAx1]) - barh(axHist,histY,histX,'EdgeColor','none','BarWidth',1,'FaceColor',.8*[1 1 1]) +end + + +%% Optimize Probabilistic +if (probabilisticOptimization == true) + numProbConstraints = sum(cellfun(@numel,{vois(:).probCFunc})); + % fmincon + options = optimoptions('fmincon',... + 'algorithm','interior-point',... interior-point, sqp + 'Display','iter-detailed',... + 'SpecifyObjectiveGradient',true,... + 'SpecifyConstraintGradient',true,... + 'AlwaysHonorConstraints', 'bounds',... + 'ScaleProblem',true);%,... + %'PlotFcn',{@optimplotx,@optimplotfunccount,@optimplotfval,@optimplotfirstorderopt,@optimplotconstrviolation,@optimplotstepsize,@optimplotfirstorderopt}); + %@(x,optimvals,state) apm_fminconPlotGradientFcn(x,optimvals,state,probFunGrad,probNonlcon)}); + %,'OutputFcn',@apm_fminconOutputFcn); + + probObjFunc = @(x) apm_fminconObjFuncWrapper(x,@(x) apm_probObjFunc(expDose_ij,covInfluence,x,vois),@(x) apm_probObjGrad(expDose_ij,covInfluence,x,vois)); + probNonlcon = @(x) apm_fminconConstrFuncWrapper(x,@(x) apm_probCFunc(expDose_ij,covInfluence,x,vois),[],@(x) apm_probCJacob(expDose_ij,covInfluence,x,vois),[]); + [wProb,fVal] = fmincon(probObjFunc,wStart,[],[],[],[],zeros(nSpots,1),Inf*ones(nSpots,1),probNonlcon,options); - nS_S = floor(nSamplesTotal/nFrac); + %minimize + %[wProb,fVal] = minimize(probObjFunc,wStart,[],[],[],[],zeros(nSpots,1),Inf*ones(nSpots,1),probNonlcon); - %Sample the systematic Part - S_S = mb_mgd(nS_S,mu,ucm.covSys)'; - doseSampleS = NaN*(ones(lowResFac*nSpots,1)*ones(1,nS_S)); + %ipopt + % funcs.objective = @(x) apm_probObjFunc(expDose_ij,covInfluence,x,vois); + % funcs.gradient = @(x) apm_probObjGrad(expDose_ij,covInfluence,x,vois); + % funcs.constraints = @(x) apm_probCFunc(expDose_ij,covInfluence,x,vois); + % funcs.jacobian = @(x) sparse(transpose(apm_probCJacob(expDose_ij,covInfluence,x,vois))); + % funcs.iterfunc = @(n,f,auxdata) drawnow('update'); + % + % startJacob = funcs.jacobian(wStart); + % funcs.jacobianstructure = @() sparse(ones(size(startJacob))); + % + % ipoptOptions.lb = zeros(nSpots,1); + % ipoptOptions.ub = Inf*ones(nSpots,1); + % ipoptOptions.cu = zeros(numProbConstraints,1); + % ipoptOptions.cl = -Inf*ones(numProbConstraints,1); + % + % ipoptOptions.ipopt.hessian_approximation = 'limited-memory'; + % %options.ipopt.print_level = 0; + % %options.ipopt.hessian_approximation = 'limited-memory'; + % ipoptOptions.ipopt.derivative_test = 'first-order'; + % [wProb,info] = ipopt(wStart,funcs,ipoptOptions); - %Now we obtain mean, error of mean and standard deviation for each sampled - %treatment - for i = 1:nS_S - %Simulate nFrac fractions - S_R = mb_mgd(nFrac,S_S(:,i)',ucm.covRand)'; - - doseSampleF = zeros(nSpots*lowResFac,nFrac); - for f = 1:nFrac - doseSampleF(:,f) = (1./ sqrt( 2*pi*ones(nSpots*lowResFac,1)*sigma.^2) .* exp( - ( bsxfun(@minus,S_R(:,f)',xLowRes).^2 ) ./ ( 2*ones(nSpots*lowResFac,1)*sigma.^2) ))*w; - end - doseSampleS(:,i) = sum(doseSampleF,2)/nFrac; - - %Plot the accumulated treatment dose sample - if (showScatterPlots) - plot(xLowRes,doseSampleS(:,i),'.','Color',[0.6 0.6 1],'MarkerSize',2); - drawnow; - end - %plot(xLowRes,std(yR),'.','Color',[1 0.6 1],'MarkerSize',2) - if (showWaitbar) - if (i == 1) - h = waitbar(i/nS_S,['Sample Treatment ' num2str(i) ' of ' num2str(nS_S)]); - else - waitbar(i/nS_S,h,['Sample Treatment ' num2str(i) ' of ' num2str(nS_S)]); - end - else - matRad_progress(i,nS_S); - end - end - if (showWaitbar) - delete(h); - end + %disp(['Final Constraint Function value: ' num2str(probConstrFunc(wProb))]); - doseSampleS_mean = mean(doseSampleS,2); - doseSampleS_std = std(doseSampleS,[],2); - doseSampleS_cov = cov(doseSampleS'); - doseSampleS_mean_stdErr = doseSampleS_std / sqrt(nS_S); - doseSampleS_std_stdErr = doseSampleS_std * sampledStdAccRel(nS_S); + %Map influence + doseProb = apm_calcDose(dose_ij,wProb); - %QI - errorbar(axProfile,xLowRes,doseSampleS_mean,doseSampleS_mean_stdErr,'rs','MarkerSize',2) - errorbar(axProfile,xLowRes,doseSampleS_std,doseSampleS_std_stdErr,'ms','MarkerSize',2) - %plot(xLowRes,kurtosis(yS),'y*','MarkerSize',10); - %plot(xLowRes,skewness(yS),'k*','MarkerSize',10); -end - - - -%Evaluate accuracy sample/analytical -if sampleValidation - resolutionFactor = highResFac / lowResFac; - ix = 1:numel(doseSampleS_mean); + expDoseProb = apm_calcExpDose(expDose_ij,wProb); - dose_mean_atSamples = expDose(ix*resolutionFactor); - dose_std_atSamples = stdDose(ix*resolutionFactor); + covDoseProb = apm_calcCovDose(covInfluence,wProb); + varDoseProb = apm_calcVarDose(varInfluence,wProb); + stdDoseProb = sqrt(varDoseProb); - rmse_mean = sqrt(mean((dose_mean_atSamples - doseSampleS_mean).^2)); - rmse_std = sqrt(mean((dose_std_atSamples - doseSampleS_std).^2)); + % For shifted beta marginals + [dmin_shiftedbeta,dmax_shiftedbeta] = apm_getShiftedBetaRange(expDoseProb,stdDoseProb,x,vois); + + % Calculate DVHs + for v = 1:numel(vois) + disp(['Computing nominal and probabilistic DVHs for VOI ' vois(v).name '...']); + vois(v).nomDVHprob = apm_DVH(doseProb(vois(v).ix),nDvhBins,1.1); + [vois(v).expDVHprob,vois(v).stdDVHprob] = apm_DVHprob(expDoseProb(vois(v).ix),covDoseProb(vois(v).ix,vois(v).ix),nDvhBins,1.1,probDvhMethod,copula_model,copula_marginals(vois(v).ix),dmin_shiftedbeta(vois(v).ix),dmax_shiftedbeta(vois(v).ix)); + end - disp(['Mean RMSE: ' num2str(rmse_mean)]); - disp(['Std RMSE: ' num2str(rmse_std)]); -end - -%% Optimize Probabilistic - -numProbConstraints = sum(cellfun(@numel,{vois(:).probCFunc})); - -% fmincon -options = optimoptions('fmincon',... - 'algorithm','interior-point',... interior-point, sqp - 'Display','iter-detailed',... - 'SpecifyObjectiveGradient',true,... - 'SpecifyConstraintGradient',true,... - 'AlwaysHonorConstraints', 'bounds',... - 'ScaleProblem',true);%,... - %'PlotFcn',{@optimplotx,@optimplotfunccount,@optimplotfval,@optimplotfirstorderopt,@optimplotconstrviolation,@optimplotstepsize,@optimplotfirstorderopt}); - %@(x,optimvals,state) apm_fminconPlotGradientFcn(x,optimvals,state,probFunGrad,probNonlcon)}); - %,'OutputFcn',@apm_fminconOutputFcn); - + %% Plot results + if ~isvalid(axProfile) + figure; + axProfile = axes(); + %axProfile = axes('Position',get(axHist,'Position'),'ActivePositionProperty','Position'); + apm_anatomyPlot(axProfile,x,vois); + end + apm_profilePlot(axProfile,x,doseProb,expDoseProb,stdDoseProb,'-'); -probObjFunc = @(x) apm_fminconObjFuncWrapper(x,@(x) apm_probObjFunc(expDose_ij,covInfluence,x,vois),@(x) apm_probObjGrad(expDose_ij,covInfluence,x,vois)); -probNonlcon = @(x) apm_fminconConstrFuncWrapper(x,@(x) apm_probCFunc(expDose_ij,covInfluence,x,vois),[],@(x) apm_probCJacob(expDose_ij,covInfluence,x,vois),[]); -[wProb,fVal] = fmincon(probObjFunc,wStart,[],[],[],[],zeros(nSpots,1),Inf*ones(nSpots,1),probNonlcon,options); - -%minimize -%[wProb,fVal] = minimize(probObjFunc,wStart,[],[],[],[],zeros(nSpots,1),Inf*ones(nSpots,1),probNonlcon); - -%ipopt -% funcs.objective = @(x) apm_probObjFunc(expDose_ij,covInfluence,x,vois); -% funcs.gradient = @(x) apm_probObjGrad(expDose_ij,covInfluence,x,vois); -% funcs.constraints = @(x) apm_probCFunc(expDose_ij,covInfluence,x,vois); -% funcs.jacobian = @(x) sparse(transpose(apm_probCJacob(expDose_ij,covInfluence,x,vois))); -% funcs.iterfunc = @(n,f,auxdata) drawnow('update'); -% -% startJacob = funcs.jacobian(wStart); -% funcs.jacobianstructure = @() sparse(ones(size(startJacob))); -% -% ipoptOptions.lb = zeros(nSpots,1); -% ipoptOptions.ub = Inf*ones(nSpots,1); -% ipoptOptions.cu = zeros(numProbConstraints,1); -% ipoptOptions.cl = -Inf*ones(numProbConstraints,1); -% -% ipoptOptions.ipopt.hessian_approximation = 'limited-memory'; -% %options.ipopt.print_level = 0; -% %options.ipopt.hessian_approximation = 'limited-memory'; -% ipoptOptions.ipopt.derivative_test = 'first-order'; -% [wProb,info] = ipopt(wStart,funcs,ipoptOptions); - - -%disp(['Final Constraint Function value: ' num2str(probConstrFunc(wProb))]); - - -%Map influence -doseProb = calcDose(dose_ij,wProb); - -expDoseProb = calcExpDose(expDose_ij,wProb); - -covDoseProb = apm_calcCovDose(covInfluence,wProb); -varDoseProb = apm_calcVarDose(varInfluence,wProb); -stdDoseProb = sqrt(varDoseProb); - -%calculate DVHs -for v = 1:numel(vois) - disp(['Computing nominal and probabilistic DVHs for VOI ' vois(v).name '...']); - vois(v).nomDVHprob = apm_DVH(doseProb(vois(v).ix),nDvhBins,1.1); - [vois(v).expDVHprob,vois(v).stdDVHprob] = apm_DVHprob(expDoseProb(vois(v).ix),covDoseProb(vois(v).ix,vois(v).ix),nDvhBins,1.1,'int_gauss_cuda'); -end - -%% Plot results -if ~isvalid(axProfile) - figure; - axProfile = axes();%axes('Position',get(axHist,'Position'),'ActivePositionProperty','Position'); - apm_anatomyPlot(axProfile,x,vois); -end -apm_profilePlot(axProfile,x,doseProb,expDoseProb,stdDoseProb,'-'); - -if ~exist('hAxProbDvh','var') || ~isvalid(hAxProbDvh) - hFigProbDvh = figure; - hAxProbDvh = axes(hFigProbDvh); -else - clear hAxProbDvh; -end - -hold(hAxProbDvh,'on'); -for v = 1:numel(vois) - if ~strcmp(vois(v).type,'BODY') || plotBodyDVH - plot(hAxProbDvh,vois(v).nomDVHprob(1,:),vois(v).nomDVHprob(2,:),'LineWidth',0.5,'LineStyle','--','Color',vois(v).dvhColor); - %dvhsProb{v} = apm_DVH(doseProb(vois(v).ix),100,1.1); + if ~exist('hAxProbDvh','var') || ~isvalid(hAxProbDvh) + hFigProbDvh = figure; + hAxProbDvh = axes(hFigProbDvh); + else + clear hAxProbDvh; + end - apm_plotProbDVH(hAxProbDvh,vois(v).expDVHprob,vois(v).stdDVHprob,vois(v).dvhColor,{'-','-.'},dvhPlotMode); + hold(hAxProbDvh,'on'); + for v = 1:numel(vois) + if ~strcmp(vois(v).type,'BODY') || plotBodyDVH + plot(hAxProbDvh,vois(v).nomDVHprob(1,:),vois(v).nomDVHprob(2,:),'LineWidth',0.5,'LineStyle','--','Color',vois(v).dvhColor); + %dvhsProb{v} = apm_DVH(doseProb(vois(v).ix),100,1.1); + + apm_plotProbDVH(hAxProbDvh,vois(v).expDVHprob,vois(v).stdDVHprob,vois(v).dvhColor,{'-','-.'},dvhPlotMode); + end + %plot(hAxProbDvh,vois(v).expDVHprob(1,:),vois(v).expDVHprob(2,:),'LineWidth',2,'LineStyle','-','Color',vois(v).dvhColor); + %hAxProbDvh.ColorOrderIndex = hAxProbDvh.ColorOrderIndex - 1; + %plot(hAxProbDvh,vois(v).expDVHprob(1,:),vois(v).expDVHprob(2,:)+vois(v).stdDVHprob(2,:),'LineWidth',1,'LineStyle','-.','Color',vois(v).dvhColor); + %hAxProbDvh.ColorOrderIndex = hAxProbDvh.ColorOrderIndex - 1; + %plot(hAxProbDvh,vois(v).expDVHprob(1,:),vois(v).expDVHprob(2,:)-vois(v).stdDVHprob(2,:),'LineWidth',1,'LineStyle','-.','Color',vois(v).dvhColor); end - %plot(hAxProbDvh,vois(v).expDVHprob(1,:),vois(v).expDVHprob(2,:),'LineWidth',2,'LineStyle','-','Color',vois(v).dvhColor); - %hAxProbDvh.ColorOrderIndex = hAxProbDvh.ColorOrderIndex - 1; - %plot(hAxProbDvh,vois(v).expDVHprob(1,:),vois(v).expDVHprob(2,:)+vois(v).stdDVHprob(2,:),'LineWidth',1,'LineStyle','-.','Color',vois(v).dvhColor); - %hAxProbDvh.ColorOrderIndex = hAxProbDvh.ColorOrderIndex - 1; - %plot(hAxProbDvh,vois(v).expDVHprob(1,:),vois(v).expDVHprob(2,:)-vois(v).stdDVHprob(2,:),'LineWidth',1,'LineStyle','-.','Color',vois(v).dvhColor); + box(hAxProbDvh,'on'); + grid(hAxProbDvh,'on'); + ylim(hAxProbDvh,[0 1]); + xlim(hAxProbDvh,[0 1.1]); + xlabel(hAxProbDvh,'rel. dose'); + ylabel(hAxProbDvh,'rel. volume'); + title('Probabilistic optimization'); + + apm_plotObjConstrInDVH(hAxProbDvh,vois,true,plotBodyDVH); + + % switch constr + % case 'DVHmin' + % plot(dvhDparam,dvhMinVol,'^k','MarkerFaceColor','k'); + % case 'DVHmax' + % plot(dvhDparam,dvhMaxVol,'vk','MarkerFaceColor','k'); + % case 'minDose' + % plot([minDose minDose],[0 1],'--k>','MarkerFaceColor','k'); + % [minDoseMu,minDoseStd] = apm_eudProb(expDoseProb(vois(1).ix),covDoseProb(vois(1).ix,vois(1).ix),-100); + % tmp_d = linspace(1e-6,vois(1).dPres*1.2,100); + % tmp_gauss = 1/sqrt(2*pi*minDoseStd^2) * exp(-0.5*(tmp_d-minDoseMu).^2 ./ minDoseStd^2); + % tmp_gauss = 0.5 * tmp_gauss / max(tmp_gauss); + % plot(tmp_d,tmp_gauss,'Color',0.5*[1 1 1]); + % case 'maxDose' + % plot([maxDose maxDose],[0 max(tmp_gauss)],'--k<','MarkerFaceColor','k'); + % [maxDoseMu,maxDoseStd] = apm_eudProb(expDoseProb(vois(2).ix),covDoseProb(vois(2).ix,vois(2).ix),100); + % tmp_d = linspace(0,1.1,100); + % tmp_gauss = 1/sqrt(2*pi*maxDoseStd^2) * exp(-0.5*(tmp_d-maxDoseMu).^2 ./ maxDoseStd^2); + % tmp_gauss = 0.5 * tmp_gauss / max(tmp_gauss); + % plot(tmp_d,tmp_gauss,'Color',0.5*[1 1 1]); + % %plot([maxDose maxDose],[0 max(tmp_gauss)],'--k<','MarkerFaceColor','k'); + % case 'EUDmin' + % [minEudMu,minEudStd] = apm_eudProb(expDoseProb(vois(1).ix),covDoseProb(vois(1).ix,vois(1).ix),eudK); + % %plot([minDose minDose],[0 1],'--k<','MarkerFaceColor','k'); + % tmp_d = linspace(0,1.1,100); + % tmp_gauss = 1/sqrt(2*pi*minEudStd^2) * exp(-0.5*(tmp_d-minEudMu).^2 ./ minEudStd^2); + % tmp_gauss = 0.5 * tmp_gauss / max(tmp_gauss); + % plot(tmp_d,tmp_gauss,'Color',0.5*[1 1 1]); + % plot([eudMin eudMin],[0 max(tmp_gauss)],'-k>','MarkerFaceColor','k'); + % case 'EUDmax' + % [maxEudMu,maxEudStd] = apm_eudProb(expDoseProb(vois(2).ix),covDoseProb(vois(2).ix,vois(2).ix),eudK); + % %plot([minDose minDose],[0 1],'--k<','MarkerFaceColor','k'); + % tmp_d = linspace(0,1.1,100); + % tmp_gauss = 1/sqrt(2*pi*maxEudStd^2) * exp(-0.5*(tmp_d-maxEudMu).^2 ./ maxEudStd^2); + % tmp_gauss = 0.5 * tmp_gauss / max(tmp_gauss); + % + % plot(tmp_d,tmp_gauss,'Color',0.5*[1 1 1]); + % + % tmp_pVal = maxEudMu + maxEudStd*sqrt(2)*erfinv(2*eudMaxProbability - 1); + % [~, tmp_ix ] = min( abs( tmp_d-tmp_pVal ) ); + % + % plot([eudMax eudMax],[0 tmp_gauss(tmp_ix)],'-k<','MarkerFaceColor','k'); + % + % nomEud = apm_eud(doseProb(vois(2).ix),vois(2).eudK); + % plot(nomEud,0,'k*','MarkerFaceColor','k'); + % + % case 'meanMax' + % [maxEudMu,maxEudStd] = apm_eudProb(expDoseProb(vois(2).ix),covDoseProb(vois(2).ix,vois(2).ix),1); + % %plot([minDose minDose],[0 1],'--k<','MarkerFaceColor','k'); + % tmp_d = linspace(0,1.1,100); + % tmp_gauss = 1/sqrt(2*pi*maxEudStd^2) * exp(-0.5*(tmp_d-maxEudMu).^2 ./ maxEudStd^2); + % tmp_gauss = 0.5 * tmp_gauss / max(tmp_gauss); + % + % plot(tmp_d,tmp_gauss,'Color',0.5*[1 1 1]); + % + % tmp_pVal = maxEudMu + maxEudStd*sqrt(2)*erfinv(2*eudMaxProbability - 1); + % [~, tmp_ix ] = min( abs( tmp_d-tmp_pVal ) ); + % + % plot([eudMax eudMax],[0 tmp_gauss(tmp_ix)],'-k<','MarkerFaceColor','k'); + % + % nomEud = apm_eud(doseProb(vois(2).ix),1); + % plot(nomEud,0,'k*','MarkerFaceColor','k'); + % otherwise + % end + % + % if contains(obj,'pwSqDev') + % plot([vois(2).probObjFunc{1}.dMax vois(2).probObjFunc{1}.dMax],[0 1],'--k<','MarkerFaceColor','k'); + % end + if sampleValidation + [doseSampleS_probOpt, nS_S_probOpt] = apm_runSampleValidation(wProb, hAxProbDvh, expDoseProb, stdDoseProb, nSamplesTotal, nFrac, mu, sigma, ucm, xLowRes, ixLowRes, xStar, axProfile, t, probDvhMethod, nDvhBins, vois, showWaitbar, showScatterPlots, showAllSampDVHs, showAvgSampDVHs); + if saveSampToFile + filename = 'SampledDosesProbOpt.mat'; + save(filename, 'doseSampleS_probOpt', 'nS_S_probOpt', 'xLowRes', 'x', 'doseProb', 'expDoseProb', 'stdDoseProb', 'vois'); + end + end end - -box(hAxProbDvh,'on'); -grid(hAxProbDvh,'on'); -ylim(hAxProbDvh,[0 1]); -xlim(hAxProbDvh,[0 1.1]); -xlabel(hAxProbDvh,'rel. dose'); -ylabel(hAxProbDvh,'rel. volume'); - -apm_plotObjConstrInDVH(hAxProbDvh,vois,true,plotBodyDVH); - -% switch constr -% case 'DVHmin' -% plot(dvhDparam,dvhMinVol,'^k','MarkerFaceColor','k'); -% case 'DVHmax' -% plot(dvhDparam,dvhMaxVol,'vk','MarkerFaceColor','k'); -% case 'minDose' -% plot([minDose minDose],[0 1],'--k>','MarkerFaceColor','k'); -% [minDoseMu,minDoseStd] = apm_eudProb(expDoseProb(vois(1).ix),covDoseProb(vois(1).ix,vois(1).ix),-100); -% tmp_d = linspace(1e-6,vois(1).dPres*1.2,100); -% tmp_gauss = 1/sqrt(2*pi*minDoseStd^2) * exp(-0.5*(tmp_d-minDoseMu).^2 ./ minDoseStd^2); -% tmp_gauss = 0.5 * tmp_gauss / max(tmp_gauss); -% plot(tmp_d,tmp_gauss,'Color',0.5*[1 1 1]); -% case 'maxDose' -% plot([maxDose maxDose],[0 max(tmp_gauss)],'--k<','MarkerFaceColor','k'); -% [maxDoseMu,maxDoseStd] = apm_eudProb(expDoseProb(vois(2).ix),covDoseProb(vois(2).ix,vois(2).ix),100); -% tmp_d = linspace(0,1.1,100); -% tmp_gauss = 1/sqrt(2*pi*maxDoseStd^2) * exp(-0.5*(tmp_d-maxDoseMu).^2 ./ maxDoseStd^2); -% tmp_gauss = 0.5 * tmp_gauss / max(tmp_gauss); -% plot(tmp_d,tmp_gauss,'Color',0.5*[1 1 1]); -% %plot([maxDose maxDose],[0 max(tmp_gauss)],'--k<','MarkerFaceColor','k'); -% case 'EUDmin' -% [minEudMu,minEudStd] = apm_eudProb(expDoseProb(vois(1).ix),covDoseProb(vois(1).ix,vois(1).ix),eudK); -% %plot([minDose minDose],[0 1],'--k<','MarkerFaceColor','k'); -% tmp_d = linspace(0,1.1,100); -% tmp_gauss = 1/sqrt(2*pi*minEudStd^2) * exp(-0.5*(tmp_d-minEudMu).^2 ./ minEudStd^2); -% tmp_gauss = 0.5 * tmp_gauss / max(tmp_gauss); -% plot(tmp_d,tmp_gauss,'Color',0.5*[1 1 1]); -% plot([eudMin eudMin],[0 max(tmp_gauss)],'-k>','MarkerFaceColor','k'); -% case 'EUDmax' -% [maxEudMu,maxEudStd] = apm_eudProb(expDoseProb(vois(2).ix),covDoseProb(vois(2).ix,vois(2).ix),eudK); -% %plot([minDose minDose],[0 1],'--k<','MarkerFaceColor','k'); -% tmp_d = linspace(0,1.1,100); -% tmp_gauss = 1/sqrt(2*pi*maxEudStd^2) * exp(-0.5*(tmp_d-maxEudMu).^2 ./ maxEudStd^2); -% tmp_gauss = 0.5 * tmp_gauss / max(tmp_gauss); -% -% plot(tmp_d,tmp_gauss,'Color',0.5*[1 1 1]); -% -% tmp_pVal = maxEudMu + maxEudStd*sqrt(2)*erfinv(2*eudMaxProbability - 1); -% [~, tmp_ix ] = min( abs( tmp_d-tmp_pVal ) ); -% -% plot([eudMax eudMax],[0 tmp_gauss(tmp_ix)],'-k<','MarkerFaceColor','k'); -% -% nomEud = apm_eud(doseProb(vois(2).ix),vois(2).eudK); -% plot(nomEud,0,'k*','MarkerFaceColor','k'); -% -% case 'meanMax' -% [maxEudMu,maxEudStd] = apm_eudProb(expDoseProb(vois(2).ix),covDoseProb(vois(2).ix,vois(2).ix),1); -% %plot([minDose minDose],[0 1],'--k<','MarkerFaceColor','k'); -% tmp_d = linspace(0,1.1,100); -% tmp_gauss = 1/sqrt(2*pi*maxEudStd^2) * exp(-0.5*(tmp_d-maxEudMu).^2 ./ maxEudStd^2); -% tmp_gauss = 0.5 * tmp_gauss / max(tmp_gauss); -% -% plot(tmp_d,tmp_gauss,'Color',0.5*[1 1 1]); -% -% tmp_pVal = maxEudMu + maxEudStd*sqrt(2)*erfinv(2*eudMaxProbability - 1); -% [~, tmp_ix ] = min( abs( tmp_d-tmp_pVal ) ); -% -% plot([eudMax eudMax],[0 tmp_gauss(tmp_ix)],'-k<','MarkerFaceColor','k'); -% -% nomEud = apm_eud(doseProb(vois(2).ix),1); -% plot(nomEud,0,'k*','MarkerFaceColor','k'); -% otherwise -% end -% -% if contains(obj,'pwSqDev') -% plot([vois(2).probObjFunc{1}.dMax vois(2).probObjFunc{1}.dMax],[0 1],'--k<','MarkerFaceColor','k'); -% end - -%% export +%% Export %cleanfigure('handle',axProfile.Parent); %cleanfigure('handle',hAxProbDvh.Parent); %cleanfigure('handle',axDVH.Parent); From 5b4daea0e1e008361a106317b26ae118f5412a69 Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 17:05:16 +0200 Subject: [PATCH 18/21] Modularize --- cProbOpt/apm_calcSampledStdAccRel.m | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 cProbOpt/apm_calcSampledStdAccRel.m diff --git a/cProbOpt/apm_calcSampledStdAccRel.m b/cProbOpt/apm_calcSampledStdAccRel.m new file mode 100644 index 0000000..40687e4 --- /dev/null +++ b/cProbOpt/apm_calcSampledStdAccRel.m @@ -0,0 +1,8 @@ +function sampledStdAccRel = apm_calcSampledStdAccRel(n) + +Kn = @(n) sqrt((n-1)/2).*exp(gammaln((n-1)/2) - gammaln(n/2)); +Vn = @(n) 2*((n-1)/2 - exp(2*gammaln(n/2) - 2*gammaln((n-1)/2))); +sampledStdAccRel = @(n) Kn(n).*sqrt(Vn(n))./sqrt(n-1); + +end + From c2cdb9b213a6f261ee5db1ff1afed0a6661af25c Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 17:19:59 +0200 Subject: [PATCH 19/21] Modularize and extend functionality for other dose distribution assumptions --- cProbOpt/apm_Lateral.m | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/cProbOpt/apm_Lateral.m b/cProbOpt/apm_Lateral.m index 8faa602..fb72dc5 100644 --- a/cProbOpt/apm_Lateral.m +++ b/cProbOpt/apm_Lateral.m @@ -83,7 +83,7 @@ samplingMethod = 'fractions'; %'independent' showWaitbar = false; showScatterPlots = false; -showAllSampDVHs = true; +showAllSampDVHs = false; showAvgSampDVHs = true; saveSampToFile = true; @@ -184,8 +184,6 @@ for v = 1:numel(vois) disp(['Computing nominal and probabilistic DVHs for VOI ' vois(v).name '...']); vois(v).nomDVH = apm_DVH(dose(vois(v).ix),nDvhBins,1.1); - expDose(vois(v).ix) - covDose(vois(v).ix) [vois(v).expDVH,vois(v).stdDVH] = apm_DVHprob(expDose(vois(v).ix),covDose(vois(v).ix,vois(v).ix),nDvhBins,1.1,probDvhMethod,copula_model,copula_marginals(vois(v).ix),dmin_shiftedbeta(vois(v).ix),dmax_shiftedbeta(vois(v).ix)); %if (v==1) % figure; From 21f14f4606b2aaae05dafa217cb0bd6a69fc906a Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 17:20:12 +0200 Subject: [PATCH 20/21] Correct plotting --- cProbOpt/apm_runSampleValidation.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cProbOpt/apm_runSampleValidation.m b/cProbOpt/apm_runSampleValidation.m index 7fb8021..9566e4c 100644 --- a/cProbOpt/apm_runSampleValidation.m +++ b/cProbOpt/apm_runSampleValidation.m @@ -104,7 +104,8 @@ axSamples = axes(); hold(axSamples, 'on'); for i = 1:nS_S - plot(axSamples,sampDVH_y_v1(i,:),sampDVH_y_v2(i,:),'LineWidth',1); + plot(axSamples,sampDVH_x(1,:),sampDVH_y_v1(i,:),'LineWidth',1); + plot(axSamples,sampDVH_x(1,:),sampDVH_y_v2(i,:),'LineWidth',1); end box(axSamples,'on'); grid(axSamples,'on'); From a6f3696bf9679fa11f3a24cec9b40489aaaa8f48 Mon Sep 17 00:00:00 2001 From: fgesualdi Date: Mon, 4 Apr 2022 17:21:28 +0200 Subject: [PATCH 21/21] Fix bug --- cProbOpt/apm_calcSampledStdAccRel.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cProbOpt/apm_calcSampledStdAccRel.m b/cProbOpt/apm_calcSampledStdAccRel.m index 40687e4..a4c5294 100644 --- a/cProbOpt/apm_calcSampledStdAccRel.m +++ b/cProbOpt/apm_calcSampledStdAccRel.m @@ -2,7 +2,7 @@ Kn = @(n) sqrt((n-1)/2).*exp(gammaln((n-1)/2) - gammaln(n/2)); Vn = @(n) 2*((n-1)/2 - exp(2*gammaln(n/2) - 2*gammaln((n-1)/2))); -sampledStdAccRel = @(n) Kn(n).*sqrt(Vn(n))./sqrt(n-1); +sampledStdAccRel = Kn(n).*sqrt(Vn(n))./sqrt(n-1); end