text:
function Lxx = nlp_hessfcn(om, x, lambda, cost_mult, Hs)
%NLP_HESSFCN  Evaluates Hessian of Lagrangian.
%   LXX = NLP_HESSFCN(OM, X, LAMBDA, COST_MULT)
%   LXX = NLP_HESSFCN(OM, X, LAMBDA, COST_MULT, HS)
%
%   Hessian evaluation function, suitable for use with MIPS, FMINCON, etc.
%
%   Inputs:
%     OM : Opt-Model object
%     X : optimization vector
%     LAMBDA (struct)
%       .eqnonlin : Lagrange multipliers on equality constraints
%       .ineqnonlin : Kuhn-Tucker multipliers on inequality constraints
%     COST_MULT : (optional) Scale factor to be applied to the cost
%          (default = 1).
%     HS : (optional) sparse matrix with tiny non-zero values specifying
%          the fixed sparsity structure that the resulting LXX should match
%
%   Outputs:
%     LXX : Hessian of the Lagrangian.
%
%   Examples:
%       Lxx = nlp_hessfcn(om, x, lambda, cost_mult);
%       Lxx = nlp_hessfcn(om, x, lambda, cost_mult, Hs);
%
%   See also NLP_COSTFCN, NLP_CONSFCN.

%   MP-Opt-Model
%   Copyright (c) 1996-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

%% ----- evaluate d2f -----
[f, df, d2f] = nlp_costfcn(om, x);
d2f = d2f * cost_mult;

%%----- evaluate Hessian of equality constraints -----
d2G = om.eval_nln_constraint_hess(x, lambda.eqnonlin, 1);

%%----- evaluate Hessian of inequality constraints -----
d2H = om.eval_nln_constraint_hess(x, lambda.ineqnonlin, 0);

%%-----  do numerical check using (central) finite differences  -----
if 0
    nx = length(x);
    step = 1e-5;
    num_d2f = sparse(nx, nx);
    num_d2G = sparse(nx, nx);
    num_d2H = sparse(nx, nx);
    for i = 1:nx
        xp = x;
        xm = x;
        xp(i) = x(i) + step/2;
        xm(i) = x(i) - step/2;
        % evaluate cost & gradients
        [fp, dfp] = nlp_costfcn(om, xp);
        [fm, dfm] = nlp_costfcn(om, xm);
        % evaluate constraints & gradients
        [Hp, Gp, dHp, dGp] = nlp_consfcn(om, xp);
        [Hm, Gm, dHm, dGm] = nlp_consfcn(om, xm);
        num_d2f(:, i) = cost_mult * (dfp - dfm) / step;
        num_d2G(:, i) = (dGp - dGm) * lambda.eqnonlin   / step;
        num_d2H(:, i) = (dHp - dHm) * lambda.ineqnonlin / step;
    end
    d2f_err = full(max(max(abs(d2f - num_d2f))));
    d2G_err = full(max(max(abs(d2G - num_d2G))));
    d2H_err = full(max(max(abs(d2H - num_d2H))));
    if d2f_err > 1e-6
        fprintf('Max difference in d2f: %g\n', d2f_err);
    end
    if d2G_err > 1e-5
        fprintf('Max difference in d2G: %g\n', d2G_err);
    end
    if d2H_err > 1e-6
        fprintf('Max difference in d2H: %g\n', d2H_err);
    end
end

Lxx = d2f + d2G + d2H;

%% force specified sparsity structure
if nargin > 4
    %% add sparse structure (with tiny values) to current matrices to
    %% ensure that sparsity structure matches that supplied
    Lxx = Lxx + Hs;

%     %% check sparsity structure against that supplied
%     if nnz(Lxx) ~= nnz(Hs)
%         fprintf('=====> nnz(Lxx) is %d, expected %d <=====\n', nnz(Lxx), nnz(Hs));
%     else
%         [iHs, jHs] = find(Hs);
%         [iH, jH] = find(Lxx);
%         if any(iH ~= iHs) || any(jH ~= jHs)
%             fprintf('=====> structure of Lxx is not as expected <=====\n');
%         end
%     end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function h=myarrow()

    %  [handle]=myarrow
    %
    %  Draws an arrow  by clicking between two points (tail, head)
    %  on the current figure, and returns the arrow's handle.
    %  Uses ARROW.
    %
    %  Richard G. Cobb 3/96

    axis(axis)
    [x,y]=ginput(2);
    start=[x(1),y(1)];
    stop=[x(2),y(2)];

    h=arrow(start,stop);
    %eof

author: CelsoReyes


text:
%--- help for dsge/accuracy ---
%
%  Solves dsge models
%
%  ::
%
%    Euler_errors=accuracy(m)
%    Euler_errors=accuracy(m,nxcuts)
%    Euler_errors=accuracy(m,nxcuts,nsims)
%    Euler_errors=accuracy(m,nxcuts,nsims,nburn)
%    Euler_errors=accuracy(m,nxcuts,nsims,nburn,reset_params)
%    [Euler_errors,outsims]=accuracy(...)
%
%  Args:
%
%     m (rise | dsge): scalar or vector of model objects.
%
%     nxcuts ({10} | integer): Number of quantiles per discretized
%       continuous shock
%
%     nsims ({1000} | integer| 1x2 cell): Number of simulations to generate
%       the state of the endogenous variables. Alternatively, if the input
%       is a cell, no simulations will be performed. In this case nsims is
%       of the form nsims={Y,regs}, where:
%       - Y : simulations of endogenous variables
%       - regs : corresponding simulations for regimes
%
%     nburn ({100} | integer): Number of burned simulations
%
%     reset_params ({true} | false): reset the parameters to their correct
%       values in evaluating the model equations.
%
%  Returns:
%     :
%
%     - **Euler_errors** [cell | struct]: if the number of models is greater
%       that one, the output is a cell array of structures. Each structure
%       has the following fields:
%       - EE [1 x nsols cell] : Euler equation errors for all the solutions
%       with each cell for one solution
%       - eqtn1, eqtn2,... : structures containing "mean" (mean of the euler
%       error across all simulations), "max" (maximum of the euler error
%       across all simulations) and "min" (minimum of the euler error across
%       all simulations)
%
%     - **outsims** [1 x 2 x nsols cell]: Simulations for the endogenous
%       variables (first element) and the regimes (second element) for each
%       solution (3rd dimension)
%
%  Note:
%
%     - The errors are log10(abs(error))
%
%     - The smaller the standard deviation of the shocks the higher the
%       accuracy even if the solution is not very accurate in the first
%       place.
%
%     - The accuracy of the solution may critically depend on the tolerance
%       level used for solving the model in the first place. For instance,
%       if an iterative algorithm (e.g. mfi) is used, the solution might not
%       be 100% accurate even in a linear model
%
%     - The errors are provided in absolute value for all equations, rather
%       than being normalized as is customarily done in the literature
%
%

author: jmaih


text:
function days = ymd_dif_common ( y1, m1, d1, y2, m2, d2 )

%*****************************************************************************80
%
%% YMD_DIF_COMMON gets the day difference between two Common YMD dates.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1, D1, the first YMD date.
%
%    Input, integer Y2, M2, D2, the second YMD date.
%
%    Output, integer DAYS, the number of days between the dates.
%
  days = 0;
%
%  Check the dates.
%
  [ y1, m1, d1, ierror ] = ymd_check_common ( y1, m1, d1 );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'YMD_DIF_COMMON - Fatal error!\n' );
    fprintf ( 1, '  Y1/M1/D1 is illegal.\n' );
    error ( 'YMD_DIF_COMMON - Fatal error!' );
  end

  [ y2, m2, d2, ierror ] = ymd_check_common ( y2, m2, d2 );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'YMD_DIF_COMMON - Fatal error!\n' );
    fprintf ( 1, '  Y2/M2/D2 is illegal.\n' );
    error ( 'YMD_DIF_COMMON - Fatal error!' );
  end

  jed1 = ymd_to_jed_common ( y1, m1, d1 );

  jed2 = ymd_to_jed_common ( y2, m2, d2 );

  days = round ( jed2 - jed1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [data,units] = compute_min_wing_angle(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  % signed minimum: left wing will be negative, right wing positive
  data{i} = min(-trx(fly).wing_anglel,trx(fly).wing_angler);

end
units = parseunits('rad');

author: kristinbranson

repo: JAABA


text:
function [varargout]=ellipsoidFit_centered(varargin)

% function [M,s,R,MU]=ellipsoidFit_centered(X,MU)
% ------------------------------------------------------------------------
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 2014/11/11
%------------------------------------------------------------------------

%%

switch nargin
    case 1
        X=varargin{1};
        MU=mean(X,1); %Point set mean
    case 2
        X=varargin{1};
        MU=varargin{2};
end

if isempty(MU)
    MU=mean(X,1); %Point set mean
end

%%

%Centre on mean
X=X-MU(ones(size(X,1),1),:); %Centre points around mean

%Compute singular value decomposition to get principal directions
[~,D,V]=svd(X,0);

%Compute input radii
r = sqrt(sum(X.^2,2));

%Create sphere coordinates
Xr=(V\X')'; %Rotate back
Xc=(D\Xr')'; %Scale back to unit sphere
[theta,phi,~] = cart2sph(Xc(:,1),Xc(:,2),Xc(:,3));
[Xc(:,1),Xc(:,2),Xc(:,3)]=sph2cart(theta,phi,ones(size(theta)));

%Transforming sphere to ellipsoid according to fit
X_fit=Xc;
X_fit=(D*X_fit')'; %Scale
X_fit=(V*X_fit')'; %Rotate
[~,~,r_fit] = cart2sph(X_fit(:,1),X_fit(:,2),X_fit(:,3));

w=mean(r./r_fit); %Determine SVD weight factor

D=D*w; %Scale singular values using weight to convert to ellipsoid semi-axis constants

%Identity
I=eye(4,4);

%Rotation
R=I;
R(1:3,1:3)=V;

%Stretch factors

s=diag(D)';

S=I;
S(1:3,1:3)=D; %Stretch

T=I;
T(1:3,4)=MU(:); %Translation

M  = T * R * S; %The transformation matrix

%%
varargout{1}=M;
varargout{2}=s;
varargout{3}=R;
varargout{4}=MU;

%%

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%  Function arranges given network (net) according to given layers to use
%  (layers_to_use), and outputs the arranged nework.
%
%  @authors: ufukefe, kutalmisince
%  Created on March 23, 2021
%  @Middle East Technical University, Center for Image Analysis
%  Last Edited on July 1, 2021

function net = ArrangeNetwork(net,layers_to_use)

    % Loads the network into the wrapper
    net = dagnn.DagNN.loadobj(net);

    % Set network to test model.
    net.mode = 'test';

    % Remove unnecessary layers
    net.removeLayer({net.layers(getLayerIndex(net, layers_to_use{1})+1:size(net.layers,2)).name})

    % Remove unused activations from RAM
    net.conserveMemory = 1;

    % Keep activations that will be used
    for i=1:size(layers_to_use,1)
    net.vars(getVarIndex(net,net.layers(getLayerIndex(net, layers_to_use{i})).outputs)).precious = 1;
    end

    % Move network to gpu
    net.move('gpu');

end

author: ufukefe

repo: DFM

save_path: github-repos/MATLAB/ufukefe-DFM


text:
function [thickness, maxThick] = CalcGrayThickness(coneRadius);

% thickness = CalcGrayThickness(coneWidth);
%
% Calculate the thickness of the gray matter using mrGray nodes. Uses a
% search-cone and nearest-neighbor approach to associate voxels from layer
% 1 with ascending layers. The highest layer node that lies within the
% search cone is defined as the thickness at that node. This
% produces a voxel-wise integer measurement of the thickness that is
% returned as a vector.
%
% Ress, 08/05

mrGlobals

% Find anatomy path:
if exist(vANATOMYPATH)
  vPath = fileparts(vANATOMYPATH);
else
  vPath = uigetdir;
end
if vPath(1) == 0, return, end

if ~exist('coneWidth', 'var'), coneRadius = 1/sqrt(2); end

% Calculate gray graph for all nodes
waitH = mrvWaitbar(0, 'Getting gray graph...');
view = initHiddenGray;
view = loadAnat(view);
view.nodes = [view.allLeftNodes, view.allRightNodes];
view.edges = [view.allLeftEdges, view.allRightEdges];
view.coords = view.nodes([2 1 3], :);
view.grayConMat = makeGrayConMat(view.nodes, view.edges, 0);
layers = view.nodes(6, :);
nLayers = max(layers);

% Determine bounding box for gray matter:
bBox = zeros(3, 2);
for ii=1:3
  bBox(ii, 1) = min(view.coords(ii, :)');
  bBox(ii, 2) = max(view.coords(ii, :)');
end

mrvWaitbar(0, waitH, 'Build classification volume...')
% Create volume containing white-matter voxels set to 1, gray-matter voxels
% set to -1:
anat = double(permute(BuildWhiteVolume, [2 1 3]));
vDims = size(anat);
anat(anat == 0) = -1; % Set non-white matter to -1

% To reduce computations, restrict volume to bounding-box dimensions and
% adjust coordinates to match:
anat = anat(bBox(1, 1):bBox(1, 2), bBox(2, 1):bBox(2, 2), bBox(3, 1):bBox(3, 2));
vDims = size(anat);
for ii=1:3, view.coords(ii, :) = view.coords(ii, :) - bBox(ii, 1) + 1; end

% Form an isosurface at their boundary, that is, at an isodensity value of
% zero, and get surface normals
mrvWaitbar(0, waitH, 'Create gray-white interface surface...');
whiteSurface = isosurface(anat, 0);
whiteSurface.normals = isonormals(anat, whiteSurface.vertices);
whiteSurface.vertices = whiteSurface.vertices(:, [2 1 3])';
whiteSurface.faces = whiteSurface.faces(:, [2 1 3])';
whiteSurface.normals = whiteSurface.normals(:, [2 1 3])';

% Advance through the gray layers and associate layer 1 nodes with
% subsequent layers that within a search cone around the white-matter
% surface normal.

% Associate layer 1 with the white isosurface normals.
L1inds = find(layers == 1);
L1 = view.coords(:, L1inds);
L1toWhiteMap = nearpoints(L1, whiteSurface.vertices);
L1normals = whiteSurface.normals(:, L1toWhiteMap);
nLength = sqrt(sum(L1normals.^2));
L1normals = L1normals ./ repmat(nLength, [3 1]);

% Loop through all layer-1 vertices.
nL1 = length(L1inds);
thickness = ones(1, nL1);
maxThick = zeros(1, nL1);
mrvWaitbar(0, waitH, 'Calculating thickness...');
for iN=1:nL1
  mrvWaitbar(iN/nL1, waitH);
  cc = L1(:, iN);
  nn = L1normals(:, iN);
  nInds = iN;
  mThick = 0;
  % Dilate in the gray matter by number of layers:
  [neighbors, nInds] = DilateGrayCoords(view, nInds, nLayers);
  nNeighbors = size(neighbors, 2) - 1;
  % Calculate unit vectors pointing from the original point to each of the
  % neighbors:
  vec = neighbors - repmat(cc, [1, nNeighbors+1]);
  lVec = sqrt(sum(vec.^2));
  keep = find(lVec > 0);
  neighbors = neighbors(:, keep);
  nInds = nInds(keep);
  uVec = vec(:, keep) ./ repmat(lVec(keep), [3 1]);
  dotP = sum(uVec .* repmat(nn, [1 nNeighbors]));
  minDotP = cos(atan(coneRadius ./ layers(nInds)));
  % Thickness is defined as the maximum layer number within the search cone:
  thickVals = layers(nInds(dotP >= minDotP));
  if ~isempty(thickVals)
    thickness(iN) = max(thickVals);
    maxThick(iN) = max(layers(nInds));
  end
end

mrvWaitbar(1, waitH, 'Saving thickness results...')
fName = fullfile(vPath, 'thickness.mat');
save(fName, 'thickness');

close(waitH);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function script_rfcn_VOC0712_ResNet101_rpn()
% script_rfcn_VOC0712_ResNet101_rpn()
% RFCN training and testing with OHEM using ResNet101 model and RPN
% proposals
% --------------------------------------------------------
% R-FCN implementation
% Modified from MATLAB Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2016, Jifeng Dai
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

clc;
clear mex;
clear is_valid_handle; % to clear init_key
run(fullfile(fileparts(fileparts(mfilename('fullpath'))), 'startup'));
%% -------------------- CONFIG --------------------
opts.caffe_version          = 'caffe_rfcn';
opts.gpu_id                 = auto_select_gpu;
active_caffe_mex(opts.gpu_id, opts.caffe_version);

% model
model                       = Model.ResNet101_for_RFCN_VOC0712();
% cache name
opts.cache_name             = 'rfcn_VOC0712_ResNet101_rpn_resnet101';
% config
conf                        = rfcn_config_simple('image_means', model.mean_image);
% train/test data
fprintf('Loading dataset...')
dataset                     = [];
dataset                     = Dataset.voc0712_trainval_sp(dataset, 'train', conf.use_flipped, 'resnet101');
dataset                     = Dataset.voc2007_test_sp(dataset, 'test', false, 'resnet101');
fprintf('Done.\n');

% do validation, or not
opts.do_val                 = true;

%% -------------------- TRAINING --------------------

opts.rfcn_model        = rfcn_train(conf, dataset.imdb_train, dataset.roidb_train, ...
                                'do_val',           opts.do_val, ...
                                'imdb_val',         dataset.imdb_test, ...
                                'roidb_val',        dataset.roidb_test, ...
                                'solver_def_file',  model.solver_def_file, ...
                                'net_file',         model.net_file, ...
                                'cache_name',       opts.cache_name, ...
                                'caffe_version',    opts.caffe_version);
assert(exist(opts.rfcn_model, 'file') ~= 0, 'not found trained model');

%% -------------------- TESTING --------------------
                          rfcn_test(conf, dataset.imdb_test, dataset.roidb_test, ...
                                'net_def_file',     model.test_net_def_file, ...
                                'net_file',         opts.rfcn_model, ...
                                'cache_name',       opts.cache_name,...
                                'ignore_cache',     true);

end

repo: R-FCN


text:
function [X,E0,R0]=totalLeastSquares(A,B,d,t)
%%TOTALLEASTSQUARES Solve the total least squares problem. That is, given
%          an mXn matrix A and an mXk matrix or vector B, find the an nXk
%          matrix or vector X such that (A+E0)*X=(B+R0), where the error
%          terms E0 and R0 are chosen to minimize the Frobenius norm of
%          D*[E0,R0]*T where D=diag(d) and T=diag(t) are nonsingular
%          diagonal matrices. This problem is essentially just solving
%          A*X=B for X where both A and B can be noisy, whereas the normal
%          least squares problem solves for the case where only B is noisy.
%
%INPUTS: A A real mXn matrix.
%        B A real nXk matrix.
%        d Optional  mX1 vector of real diagonal elements for the matrix D,
%          which affects how the errors are weighted. If omitted, a vector
%          of ones is used (D is just an identity matrix).
%        t Optional (n+k)X1 vector of real diagonal elements for the matrix
%          T, which affects show errors are weighted. If omitted, a vector
%          of ones is used (T is just an identity matrix).
%
%OUTPUTS: X The nXk real solution of the total least squares problem.
%     E0,R0 The error matrices associated with the solution to the problem.
%           They are such that (A+E0)*X=(B+R0) and the Frobenius norm of
%           D*[E0,R0]*T, which is computed via norm(D*[E0,R0]*T,'fro'), is
%           minimized.
%
%The algorithm for solving the total least squares problem is taken from
%Chapter 6.3.1 of [1].
%
%REFERENCES:
%[1] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed.
%    Baltimore: Johns Hopkins University Press, 2013.
%
%April 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

m=size(A,1);
n=size(A,2);
k=size(B,2);

if(nargin<4)
    t=ones(n+k,1);
end

if(nargin<3)
    d=ones(m,1);
end

D=diag(d);
T=diag(t);

[U,Sigma,V]=svd(D*[A,B]*T);

V12=V(1:n,(n+1:end));
V22=V((n+1):end,(n+1):end);

T1=diag(t(1:n));
T2=diag(t((n+1):end));

X=-T1*V12/V22/T2;

if(nargout>1)
    U2=U(:,(n+1):end);
    Sigma2=Sigma((n+1):end,(n+1):end);

    ERMat=-D\(U2*Sigma2*[V12',V22'])/T;

    E0=ERMat(:,1:n);
    R0=ERMat(:,(n+1):end);
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% clc; clear all; close all;

%% hydro data
hydro = struct();
hydro = readWAMIT(hydro,'rm3.out',[]);
hydro = radiationIRF(hydro,60,[],[],[],[]);
hydro = radiationIRFSS(hydro,[],[]);
hydro = excitationIRF(hydro,157,[],[],[],[]);
writeBEMIOH5(hydro)

%% Plot hydro data
% plotBEMIO(hydro)

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function [ x, seed ] = cardioid_sample ( a, b, seed )

%*****************************************************************************80
%
%% CARDIOID_SAMPLE samples the Cardioid PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 <= B <= 0.5.
%
%    Input/output, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the PDF.
%    A - PI <= X <= A + PI.
%
  [ cdf, seed ] = r8_uniform_01 ( seed );

  x = cardioid_cdf_inv ( cdf, a, b );

  return
end


text:
classdef PTKLobesFromFissurePlane < PTKPlugin
    % PTKLobesFromFissurePlane. Plugin which is part of the lobar segmentation.
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     PTKLobesFromFissurePlane is an intermediate stage in segmenting the
    %     lobes.
    %
    %     For more information, see
    %     [Doel et al., Pulmonary lobe segmentation from CT images using
    %     fissureness, airways, vessels and multilevel B-splines, 2012]
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Lobes (fissure plane)'
        ToolTip = ''
        Category = 'Lobes'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = true
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
    end

    methods (Static)
        function results = RunPlugin(application, reporting)

            left_and_right_lungs = application.GetResult('PTKLeftAndRightLungs');
            fissure_plane = application.GetResult('PTKFissurePlane');
            lung_mask = application.GetResult('PTKLeftAndRightLungs');
            left_lung_template = application.GetTemplateImage(PTKContext.LeftLung).BlankCopy;
            right_lung_template = application.GetTemplateImage(PTKContext.RightLung).BlankCopy;
            results_left = PTKLobesFromFissurePlane.GetLeftLungResults(left_lung_template, lung_mask.Copy, fissure_plane.Copy, reporting);
            results_right = PTKLobesFromFissurePlane.GetRightLungResults(right_lung_template, lung_mask, fissure_plane, reporting);

            results = PTKCombineLeftAndRightImages(application.GetTemplateImage(PTKContext.LungROI), results_left, results_right, left_and_right_lungs);
            results.ImageType = PTKImageType.Colormap;
        end

        function results = GenerateImageFromResults(results, ~, ~)
        end
    end

    methods (Static, Access = private)
        function left_results = GetLeftLungResults(lung_template, lung_mask, fissure_plane, reporting)

            lung_mask.ChangeRawImage(uint8(lung_mask.RawImage == 2));

            lung_mask.ResizeToMatch(lung_template);

            fissure_plane.ResizeToMatch(lung_template);
            fissure_plane = find(fissure_plane.RawImage(:) == 4);

            left_results = PTKDivideVolumeUsingScatteredPoints(lung_mask, fissure_plane, 5, reporting);
            left_results.ChangeColourIndex(1, 5);
            left_results.ChangeColourIndex(2, 6);
        end

        function results_right = GetRightLungResults(lung_template, lung_mask, fissure_plane, reporting)

            lung_mask.ChangeRawImage(uint8(lung_mask.RawImage == 1));

            lung_mask.ResizeToMatch(lung_template);
            fissure_plane.ResizeToMatch(lung_template);
            fissure_plane_o = find(fissure_plane.RawImage(:) == 3);

            results_right = PTKDivideVolumeUsingScatteredPoints(lung_mask, fissure_plane_o, 5, reporting);
            results_right.ChangeColourIndex(2, 4);

            % Mid lobe
            fissure_plane_m = find(fissure_plane.RawImage(:) == 2);

            if ~isempty(fissure_plane_m)
                lung_mask_excluding_lower = lung_mask.Copy;
                lung_mask_excluding_lower.ChangeRawImage(results_right.RawImage == 1);

                results_mid_right = PTKDivideVolumeUsingScatteredPoints(lung_mask_excluding_lower, fissure_plane_m, 20, reporting);
                results_right.ChangeSubImageWithMask(results_mid_right, results_mid_right);
            else
                reporting.ShowWarning('PTKLobesFromFissurePlane:NoRightObliqueFissure', 'Unable to find the right horizontal fissure. No middle right lobe segmentation will be shown.', []);
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%MDL_COIL Create model of a coil manipulator
%
% MDL_COIL creates the workspace variable coil which describes the
% kinematic characteristics of a serial link manipulator with 50 joints
% that folds into a helix shape.
%
% MDL_BALL(N) as above but creates a manipulator with N joints.
%
% Also defines the workspace vectors:
%   q  joint angle vector for default helical configuration
%
% Reference::
% - "A divide and conquer articulated-body algorithm for parallel O(log(n))
%   calculation of rigid body dynamics, Part 2",
%   Int. J. Robotics Research, 18(9), pp 876-892.
%
% Notes::
% - Unlike most other mdl_xxx scripts this one is actually a function that
%   behaves like a script and writes to the global workspace.
%
% See also mdl_ball, SerialLink.

% MODEL: generic, coil, hyper redundant, 50DOF, standard_DH

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function r = mdl_coil(N)

    if nargin == 0
        N = 50;
    end

    % create the links
    for i=1:N
        links(i) = Link('d', 0, 'a', 1/N, 'alpha', 5*pi/N);
    end

    % and build a serial link manipulator
    robot = SerialLink(links, 'name', 'coil');

    % place the variables into the global workspace
    if nargin == 1
        r = robot;
    elseif nargin == 0
        assignin('caller', 'coil', robot);
        assignin('caller', 'q', 10*pi/N*ones(1,N));
    end

author: petercorke

repo: robotics-toolbox-matlab


text:
function f = ctranspose(f) %#ok<*INUSD>
%CTRANSPOSE   SINGFUN objects are not transposable.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

error('CHEBFUN:SINGFUN:ctranspose:notPossible', ...
      'SINGFUN objects are not transposable.')

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function route = perturb(route_old, method)
% PERTURB
% route = PERTURB(route_old, method) generate randomly a neighbouring route by
% perturb old route. perturb methods:
%                        ___________            ___________
%     1. reverse:   [1 2 3 4 5 6 7 8 9] -> [1 2 8 7 6 5 4 3 9]
%                        _         _            _         _
%     2. swap:      [1 2 3 4 5 6 7 8 9] -> [1 2 8 4 5 6 7 3 9]

route = route_old;
numbercities = length(route);
city1 = ceil(numbercities*rand);
city2 = ceil(numbercities*rand);
switch method
    case 'reverse'
        citymin = min(city1,city2);
        citymax = max(city1,city2);
        route(citymin:citymax) = route(citymax:-1:citymin);
    case 'swap'
        route([city1, city2]) = route([city2, city1]);
end

author: HuangCongQing


text:
function [L,E,err,iter] = trpca_snn(X,alpha,opts)

% Solve the Tensor Robust Principal Component Analysis (TRPCA) based on Sum of Nuclear Norm (SNN) problem by M-ADMM
%
% min_{L,E} \sum_i \alpha_i*||L_{i(i)}||_* + ||E||_1,
% s.t. X = L + E.
%
% ---------------------------------------------
% Input:
%       X       -    d1*d2*...dk tensor
%       alpha   -    k*1 vector, parameters
%       opts    -    Structure value in Matlab. The fields are
%           opts.tol        -   termination tolerance
%           opts.max_iter   -   maximum number of iterations
%           opts.mu         -   stepsize for dual variable updating in ADMM
%           opts.max_mu     -   maximum stepsize
%           opts.rho        -   rho>=1, ratio used to increase mu
%           opts.DEBUG      -   0 or 1
%
% Output:
%       L       -    d1*d2*...*dk tensor
%       E       -    d1*d2*...*dk tensor
%       err     -    residual
%       iter    -    number of iterations
%
% version 1.0 - 24/06/2016
%
% Written by Canyi Lu (canyilu@gmail.com)
%

tol = 1e-8;
max_iter = 500;
rho = 1.1;
mu = 1e-4;
max_mu = 1e10;
DEBUG = 0;

if ~exist('opts', 'var')
    opts = [];
end
if isfield(opts, 'tol');         tol = opts.tol;              end
if isfield(opts, 'max_iter');    max_iter = opts.max_iter;    end
if isfield(opts, 'rho');         rho = opts.rho;              end
if isfield(opts, 'mu');          mu = opts.mu;                end
if isfield(opts, 'max_mu');      max_mu = opts.max_mu;        end
if isfield(opts, 'DEBUG');       DEBUG = opts.DEBUG;          end

dim = size(X);
k = length(dim);

E = zeros(dim);
Y = cell(k,1);
L = Y;
for i = 1 : k
    Y{i} = E;
    L{i} = E;
end

iter = 0;
for iter = 1 : max_iter
    Lk = L;
    Ek = E;
    % first super block {L_i}
    sumtemp = zeros(dim);
    for i = 1 : k
        L{i} = Fold(prox_nuclear(Unfold(X-E-Y{i}/mu,dim,i), alpha(i)/mu),dim,i);
        sumtemp = sumtemp + L{i} + Y{i}/mu;
    end
    % second super block {E}
    E = prox_l1(X-sumtemp/k,1/(mu*k));

    chg = max(abs(Ek(:)-E(:)));
    err = 0;
    for i = 1 : k
        dY = L{i}+E-X;
        err = err+norm(dY(:))^2;
        Y{i} = Y{i}+mu*dY;
        chg = max([chg, max(abs(dY(:))), max(abs(Lk{i}(:)-L{i}(:)))]);
    end
    err = sqrt(err);

    if DEBUG
        if iter == 1 || mod(iter, 10) == 0
            disp(['iter ' num2str(iter) ', mu=' num2str(mu) ...
                    ', err=' num2str(err)]);
        end
    end
    if chg < tol
        break;
    end
    mu = min(rho*mu,max_mu);
end
L = L{1};

author: canyilu

repo: LibADMM-toolbox


text:
function opt= defopt_scalp_r(varargin)
%DEFOPT_SCALP_R - Properties for plot_scalp* functions, optimized for
%visualizing (signed, squared) r values

props= {'Extrapolation'          1
        'ExtrapolateToMean'      0
        'Resolution'             71
        'CLim'                   'sym'
        'ShrinkColorbar'         0.2
        'ColorOrder'             [0 0 0]
        'Colormap'               cmap_posneg(51)
        'LineWidth'              2
        'ChannelLineStyleOrder'  {'thick' 'thin'}
        'IvalColor'              [0.8 0.8 0.8; 0.6 0.6 0.6]
        'Contour'                5
        'ContourPolicy'          'withinrange'
        'ContourLineprop'        {'LineWidth' 0.3}
        'MarkMarkerProperties'   {'MarkerSize' 6, 'LineWidth' 1}
        'ChannelAtBottom'        1
        'GlobalCLim'             1
        'LegendPos'              'NorthWest'};

opt= opt_proplistToStruct(varargin{:});
opt= opt_setDefaults(opt, props);

author: bbci


text:
function opts = mkopts_cdp(K);

opts.algorithm = 'cdp';
opts.use_kd_tree = 0;
opts.initial_K = K;
opts.do_greedy = 0;
opts.do_split = 0;
opts.do_merge = 0;
opts.do_sort = 0;

% Local Variables: ***
% mode: matlab ***
% End: ***

repo: REST


text:
function [ a, b ] = p06_lim ( dim_num )

%*****************************************************************************80
%
%% P06_LIM returns the integration limits for problem 06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Output, real A(DIM_NUM), B(DIM_NUM), the lower and upper
%    limits of integration.
%
  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  return
end


text:
function h = hashmd5(inp)
% Computes the MD5 hash of input data.
%
% function h = hashmd5(inp)
%
% Returns a string containing the MD5 hash of the input variable. The input
% variable may be of any class that can be typecast to uint8 format, which
% is fairly non-restrictive.

% This file is part of Manopt: www.manopt.org.
% This code is a stripped version of more general hashing code by
% Michael Kleder, Nov 2005.
% Change log:
%
%   Aug. 8, 2013 (NB):
%       Made x a static (persistent) variable, in the hope it will speed
%       it up. Furthermore, the function is now Octave compatible.

    is_octave = exist('OCTAVE_VERSION', 'builtin');

    persistent x;
    if isempty(x) && ~is_octave
        x = java.security.MessageDigest.getInstance('MD5');
    end

    inp=inp(:);
    % Convert strings and logicals into uint8 format
    if ischar(inp) || islogical(inp)
        inp=uint8(inp);
    else % Convert everything else into uint8 format without loss of data
        inp=typecast(inp,'uint8');
    end

    % Create hash
    if ~is_octave
        x.update(inp);
        h = typecast(x.digest, 'uint8');
        h = dec2hex(h)';
        % Remote possibility: all hash bytes < 128, so pad:
        if(size(h,1))==1
            h = [repmat('0',[1 size(h,2)]);h];
        end
        h = lower(h(:)');
    else
        h = md5sum(char(inp'), true);
    end

end

author: NicolasBoumal

repo: manopt


text:
function s2de_test ( )

%*****************************************************************************80
%
%% S2DE_TEST tests the S2DE library.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/stokes_2d_exact/s2de_test.m
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 February 2015
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'S2DE_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the S2DE library.\n' );

  uvp_stokes1_test ( );
  resid_stokes1_test ( );
  gnuplot_stokes1_test ( );

  uvp_stokes2_test ( );
  resid_stokes2_test ( );
  gnuplot_stokes2_test ( );

  uvp_stokes3_test ( );
  resid_stokes3_test ( );
  gnuplot_stokes3_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'S2DE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  timestamp ( );

  return
end


text:
function [out] = interception_3(p1,In)
%interception_3

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Interception excess after a fraction is intercepted
% Constraints:  -
% @(Inputs):    p1   - fraction throughfall [-]
%               In   - incoming flux [mm/d]

out = p1.*In;

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function [ know, x ] = p24_sol ( n )

%*****************************************************************************80
%
%% P24_SOL returns the solution for problem 24.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the problem.  This value
%    is only needed for those problems with variable N.
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Output, real X(N), the solution, if known.
%
  know = 1;

  x(1:n) = -5.0;

  return
end


text:
function HaarTransformationMatrix=ConstructHaarWaveletTransformationMatrix(WidthOfSquareMatrix)
%--------------------------------------------------------------------------
%Create Haarwavelet transformation matrix H for the matrix vector
%mulplication implimentation of Haar wavelet transformation.
%This function uses the following nice formula to create the Haar
%transformation matrix:
%               H_n=1/sqrt(2)[H_(n/2) kron (1 1)
%                             I_(n/2) kron (1 -1)],
%                              where 'kron' denotes the kronecker product.
%The iteration starts with H_1=[1]. The normalization constant 1/sqrt(2)
%ensure that H_n^T*H_n=I, where I is identity matrix. Haar wavelets are the
%rows of H_n.
%--------------------------------------------------------------------------
%function HaarTransformationMatrix=ConstructHaarWaveletTransformationMatrix(WidthOfSquareMatrix)
% Input:
%       WidthOfSquareMatrix: the width of sqaure Haar wavelet
%                                  transformation matrix.
% Output:
%       HaarTransformationMatrix: Ceated Haar transformation matrix and
%       it's size is the power of 2,i.e., 2, 4,
%       8,16,32,64,128,256,512,1024,2048,4096,etc.
%--------------------------------------------------------------------------
%Author: Jin Qi
%Email: jqichina@hotmail.com
%Date:  11/4/2011
%--------------------------------------------------------------------------
%Test function for one dimentional signal
%
% n=8;
% S=rand(n,1); % for one dimentional signal
% H=ConstructHaarWaveletTransformationMatrix(n);
% C=H*S; %Haar wavelet transformation
% RS=H'*C;% Reconstruct signal
% fprintf('reconstruction error is %f',norm(S-RS));%plot reconstruction error
%--------------------------------------------------------------------------
%Test for two dimentional signal,ie. image
%
% n=8;
% S=rand(n,n); % for one dimentional signal
% H=ConstructHaarWaveletTransformationMatrix(n);
% C=H*S*H'; %Haar wavelet transformation
% RS=H'*C*H;% Reconstruct signal
% fprintf('reconstruction error is %f',norm(S-RS,'fro'));%plot reconstruction error

n=WidthOfSquareMatrix; % copy the parameter

% check input parameter and make sure it's the power of 2
Level=log2(n);
if 2^Level<n, error('please ensure the value of input parameter is the power of 2');end

%Initialization
H=[1];
NC=1/sqrt(2);%normalization constant
LP=[1 1];
HP=[1 -1];

% iteration from H=[1]
for i=1:Level
    H=NC*[kron(H,LP);kron(eye(size(H)),HP)];
end
HaarTransformationMatrix=H;

author: Sable

repo: mcbench-benchmarks


text:
function cvx_optval = det_rootn( X )

%DET_ROOTN   Internal cvx version.

error( nargchk( 1, 1, nargin ) );
n = size( X, 1 );
if ndims( X ) > 2,

    error( 'N-D arrays are not supported.' );

elseif size( X, 2 ) ~= n,

    error( 'Matrix must be square.' );

elseif nnz( X ) <= n && nnz( diag( X ) ) == nnz( X ),

    cvx_optval = geo_mean( diag( X ) );

elseif cvx_isconstant( X ),

    cvx_optval = cvx( det_rootn( cvx_constant( X ) ) );

elseif isreal( X ),

    cvx_begin
        variable Z(n,n) lower_triangular
        D = diag( Z );
        maximize( geo_mean( D ) );
        subject to
            [ diag( D ), Z' ; Z, X ] == semidefinite(2*n);
    cvx_end

else

    cvx_begin
        variable Z(n,n) lower_triangular complex
        D = diag( Z );
        maximize( geo_mean( real( D ) ) );
        subject to
            [ diag( D ), Z' ; Z, X ] == hermitian_semidefinite(2*n);
    cvx_end

end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.


text:
function update_plotTMPCPrediction(ax, handles, tPred, plot_tube, center_points, M_s, be_l, be_u, K_t, ...
    isConstrained, ub, lb, direction, isInput, TargetAvailable, target, LinearizationAvailable, lin)
% function plotStatePrediction
% Authors:      Martin Euler
%               Alexander Wischnewski
% Description:
%   helper function used to plot state predictions
% Inputs/parameters:
%   tPlot: vector containing the time window for plotting [tStart, tEnd]
%   tPred: vector with the time instants corresponding to the prediction
%   plot_tube: flag for tube plots
%   center_points: data for center points
%   M_s: tube shape matrices
%   be_l: lower bound of terminal set constraints
%   be_u: upper bound of terminal set constraints
%   K_t: tube controller
%   isConstrained: flag whether the signal is constrained or not
%   ub: upper limit for the signal
%   lb: lower limit for the signal
%   isInput: set to true in case this is an input (different calculation of tubes)
%   TargetAvailable: target trajectory available and should be plotted
%   target: actual target trajectory
%   LinearizationAvailable: linearization trajectory available and should be plotted
%   lin: actual linearization trajectory

% plot predicted signal
set(handles.pred, 'XData', tPred, 'YData', center_points);
% plot tube and terminal sets
if plot_tube
    % get tube bounds
    [y_low, y_up] = calcBoundsPlot(center_points, M_s, isInput, K_t, length(tPred), direction);
    % plot tube lower bounds
    set(handles.tube_ub, 'XData', tPred,'YData', y_low(1,:));
    % plot tube upper bounds
    set(handles.tube_lb, 'XData', tPred,'YData', y_up(1,:));
    % plot terminal set lower bounds and use a virtual time interval of one second (half a second
    % forth and back) to visualize it properly
    set(handles.terminal_lb, 'XData', [tPred(end)-0.5, tPred(end)+0.5],'YData',  be_l*ones(2, 1));
    % plot terminal set upper bounds
    set(handles.terminal_ub, 'XData', [tPred(end)-0.5, tPred(end)+0.5],'YData',  be_u*ones(2, 1));
end
if isConstrained
    set(handles.const_lb, 'XData', tPred, 'YData', ub);
    set(handles.const_ub, 'XData', tPred, 'YData', lb);
end
if TargetAvailable
    set(handles.target, 'XData', tPred, 'YData', target);
end
if LinearizationAvailable
    set(handles.lin, 'XData', tPred, 'YData', lin);
end

% adjust x limits appropriately such that it covers a certain range
xlim(ax, [tPred(1)-3, tPred(1)+5]);

end

author: TUMFTM


text:
function Offspring = OperatorGAhalf(Problem,Parent,Parameter)
%OperatorGAhalf - Crossover and mutation operators of genetic algorithm.
%
%   This function is the same to OperatorGA, while only the first half of
%   the offsprings are evaluated and returned.
%
%   See also OperatorGA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    if nargin > 2
        [proC,disC,proM,disM] = deal(Parameter{:});
    else
        [proC,disC,proM,disM] = deal(1,20,1,20);
    end
    if isa(Parent(1),'SOLUTION')
        evaluated = true;
        Parent    = Parent.decs;
    else
        evaluated = false;
    end
    Parent1   = Parent(1:floor(end/2),:);
    Parent2   = Parent(floor(end/2)+1:floor(end/2)*2,:);
    Offspring = zeros(size(Parent1,1),size(Parent1,2));
    Type      = arrayfun(@(i)find(Problem.encoding==i),1:5,'UniformOutput',false);
    if ~isempty([Type{1:2}])    % Real and integer variables
        Offspring(:,[Type{1:2}]) = GAreal(Parent1(:,[Type{1:2}]),Parent2(:,[Type{1:2}]),Problem.lower([Type{1:2}]),Problem.upper([Type{1:2}]),proC,disC,proM,disM);
    end
    if ~isempty(Type{3})        % Label variables
        Offspring(:,Type{3}) = GAlabel(Parent1(:,Type{3}),Parent2(:,Type{3}),proC,proM);
    end
    if ~isempty(Type{4})        % Binary variables
        Offspring(:,Type{4}) = GAbinary(Parent1(:,Type{4}),Parent2(:,Type{4}),proC,proM);
    end
    if ~isempty(Type{5})        % Permutation variables
        Offspring(:,Type{5}) = GApermutation(Parent1(:,Type{5}),Parent2(:,Type{5}),proC);
    end
    if evaluated
        Offspring = Problem.Evaluation(Offspring);
    end
end

function Offspring = GAreal(Parent1,Parent2,lower,upper,proC,disC,proM,disM)
% Genetic operators for real and integer variables

    %% Simulated binary crossover
    [N,D] = size(Parent1);
    beta  = zeros(N,D);
    mu    = rand(N,D);
    beta(mu<=0.5) = (2*mu(mu<=0.5)).^(1/(disC+1));
    beta(mu>0.5)  = (2-2*mu(mu>0.5)).^(-1/(disC+1));
    beta = beta.*(-1).^randi([0,1],N,D);
    beta(rand(N,D)<0.5) = 1;
    beta(repmat(rand(N,1)>proC,1,D)) = 1;
    Offspring = (Parent1+Parent2)/2+beta.*(Parent1-Parent2)/2;

    %% Polynomial mutation
    Lower = repmat(lower,N,1);
    Upper = repmat(upper,N,1);
    Site  = rand(N,D) < proM/D;
    mu    = rand(N,D);
    temp  = Site & mu<=0.5;
    Offspring       = min(max(Offspring,Lower),Upper);
    Offspring(temp) = Offspring(temp)+(Upper(temp)-Lower(temp)).*((2.*mu(temp)+(1-2.*mu(temp)).*...
                      (1-(Offspring(temp)-Lower(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1))-1);
    temp = Site & mu>0.5;
    Offspring(temp) = Offspring(temp)+(Upper(temp)-Lower(temp)).*(1-(2.*(1-mu(temp))+2.*(mu(temp)-0.5).*...
                      (1-(Upper(temp)-Offspring(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1)));
end

function Offspring = GAlabel(Parent1,Parent2,proC,proM)
% Genetic operators for label variables

    %% Uniform crossover
    [N,D] = size(Parent1);
    k     = rand(N,D) < 0.5;
    k(repmat(rand(N,1)>proC,1,D)) = false;
    Offspring    = Parent1;
    Offspring(k) = Parent2(k);

    %% Bitwise mutation
    Site = rand(N,D) < proM/D;
    Rand = randi(D,N,D);
    Offspring(Site) = Rand(Site);

    %% Repair
    for i = 1 : N
        Off = zeros(1,D);
        while ~all(Off)
            x = find(~Off,1);
            Off(Offspring(i,:)==Offspring(i,x)) = max(Off) + 1;
        end
        Offspring(i,:) = Off;
    end
end

function Offspring = GAbinary(Parent1,Parent2,proC,proM)
% Genetic operators for binary variables

    %% Uniform crossover
    [N,D] = size(Parent1);
    k     = rand(N,D) < 0.5;
    k(repmat(rand(N,1)>proC,1,D)) = false;
    Offspring    = Parent1;
    Offspring(k) = Parent2(k);

    %% Bit-flip mutation
    Site = rand(N,D) < proM/D;
    Offspring(Site) = ~Offspring(Site);
end

function Offspring = GApermutation(Parent1,Parent2,proC)
% Genetic operators for permutation variables

    %% Order crossover
    [N,D]     = size(Parent1);
    Offspring = Parent1;
    k         = randi(D,1,N);
    for i = 1 : N
        if rand < proC
            Offspring(i,k(i)+1:end) = setdiff(Parent2(i,:),Parent1(i,1:k(i)),'stable');
        end
    end

    %% Slight mutation
    k = randi(D,1,N);
    s = randi(D,1,N);
    for i = 1 : N
        if s(i) < k(i)
            Offspring(i,:) = Offspring(i,[1:s(i)-1,k(i),s(i):k(i)-1,k(i)+1:end]);
        elseif s(i) > k(i)
            Offspring(i,:) = Offspring(i,[1:k(i)-1,k(i)+1:s(i)-1,k(i),s(i):end]);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
close all;
clear all;
clc;
rng('default');

png_export = true;
pdf_export = false;

load ('bin/figure_1_spherical_dict_model_1_somp_success_with_k.mat');

mf = spx.graphics.Figures();

mf.new_figure('Recovery probability with K for SOMP');
hold all;
legends = cell(1, num_ss);
for ns=1:num_ss
    S = Ss(ns);
    plot(Ks, bp_success_with_k(ns, :));
    legends{ns} = sprintf('S=%d', S);
end
grid on;
xlabel('Sparsity Level');
ylabel('Empirical Recovery Rate');
legend(legends);

author: indigits

repo: sparse-plex


text:
classdef UnilateralConstraint < handle
    % UnilateralConstraint represents a scalar or vector inequality
    % constraints should be imposed on the continuous dynamical systems
    %
    % @author ayonga @date 2017-04-20
    %
    % Copyright (c) 2016-2017, AMBER Lab
    % All right reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted only in compliance with the BSD 3-Clause
    % license, see
    % http://www.opensource.org/licenses/bsd-license.php

    methods

        function obj = UnilateralConstraint(model, h, name, deps, varargin)
            % The class constructor function
            %
            % Parameters:
            % model: the dynamical system model in which the virtual
            % constraints are defined @type DynamicalSystem
            % h: the symbolic expression of the constraints
            % name: the name of the virtual constraints @type char
            % deps: the dependent variables (could be states or inputs) @type
            % cellstr
            % @type SymExpression
            % varargin: optional parameters. In details
            %  AuxData:
            %  ConstrLabel:

            if nargin == 0
                return;
            end

            % validate (model) argument
            validateattributes(model, {'DynamicalSystem'},...
                {'scalar'},...
                'UnilateralConstraint','model');

            % validates (name) argument
            validateName(obj, name);
            obj.Name = name;

            % validate (ya) argument
            validateattributes(h,{'SymExpression'},...
                {'nonempty','vector'},...
                'UnilateralConstraint','h');
            if isrow(h) % convert to column vector if it is a row vector
                h = vertcat(h(:));
            end
            dim = length(h);
            obj.Dimension = dim;

            % parse the input options
            args = struct(varargin{:});
            assert(isscalar(args),...
                'The values of optional properties are must be scalar data.');

            % validate and assign the label
            if isfield(args, 'ConstrLabel')
                obj.setConstrLabel(args.ConstrLabel);
            end

            if isfield(args, 'AuxData')
                auxdata = args.AuxData;
            else
                auxdata = [];
            end

            %% set the unilateral constraints
            if ~iscell(deps), deps = {deps}; end
            n_deps = length(deps);

            % get the variable group of each dependent variable
            var_group = cellfun(@(x)model.validateVarName(x), deps, 'UniformOutput', false);
            vars = cell(1,n_deps);

            % assume it is not input-dependent
            input_dep = false;
            for i=1:n_deps
                tmp = var_group{i};
                % check if it is input-dependent
                if strcmp(tmp{1},'Inputs')
                    input_dep = true;
                    vars{i} = model.(tmp{1}).(tmp{2}).(deps{i});
                else
                    vars{i} = model.(tmp{1}).(deps{i});
                end
            end
            obj.InputDependent = input_dep;

            if isa(h, 'SymFunction')
                % validate the dependent variables
                assert(length(h.Vars)==n_deps,...
                    ['The constraint SymFunction (h) must have the same',...
                    'number of dependent variables as specified in the argument (deps).']);
                for i=1:n_deps
                    assert(h.Vars{i} == vars{i},...
                        'The %d-th dependent variable must be the same as the variable %s.',i,deps{i});
                end
                obj.h_ = h;
                obj.DepLists = deps;

                if ~isempty(h.Params)
                    if isempty(auxdata)
                        error('The constraint function requires auxilary constant parameters.');
                    else
                        if ~iscell(auxdata), auxdata = {auxdata}; end
                        assert(length(h.Params) == length(auxdata),...
                            'The number of required auxilaray data (%d) and the provided auxilary data (%d) does not match.\n',...
                            length(h.Params),length(auxdata));
                        for i=1:numel(auxdata)
                            assert(numel(auxdata{i})== prod(size(h.Params{i})),...
                                'The length %d-th auxiliary parameters variable does not match the function definition.',i); %#ok<PSIZE>
                        end

                        obj.AuxData = auxdata;
                    end

                else
                    obj.AuxData = [];
                end

            elseif isa(h, 'SymExpression')
                % create a SymFunction object if the input is a SymExpression

                obj.h_ = SymFunction(['u_' name '_' model.Name], h, vars);
                obj.DepLists = deps;
                if ~isempty(auxdata)
                    error('If the constraint depends on auxilary constant parameters, please provide it as a SymFunction object directly.');
                end
                obj.AuxData = [];
            else
                error('The constraint expression must be given as an object of SymExpression or SymFunction.');
            end

        end

    end

    % properties determined internally
    properties (SetAccess=protected, GetAccess=public)
        % The dimension of the virtual constraints
        %
        % @type integer
        Dimension

        % An indicator whether the unilateral constraints depend on the
        % input variables of the dynamical system. If the constraints are
        % input-dependent, then it will requires to execute the
        % calcDynamics function to compute the required input signals.
        %
        % @type logical
        InputDependent
    end

    % properties must be determined by the users
    properties (SetAccess=protected, GetAccess=public)
        % The name of the virtual constraints
        %
        % @type char
        Name

        % The label of the holonomic constraint
        %
        % @type char
        ConstrLabel

        % The name list of dependent variables in the associated dynamical
        % system model
        %
        % @type cellstr
        DepLists

        % The list of auxiliary data to be used to call the function
        %
        % @type cell
        AuxData

    end

    properties (Dependent)
        % The holonomic constraint expression
        %
        % @type SymFunction
        ConstrExpr

    end

    methods
        function cstr = get.ConstrExpr(obj)
            cstr = obj.h_;
        end

    end

    properties (Access = protected)

        % The holonomic constraint expression
        %
        % @type SymFunction
        h_

    end

    methods
        function export(obj, export_path, varargin)
            % export the symbolic expressions of the constraints and
            % compile as MEX files.
            %
            % Parameters:
            %  export_path: the path to export the file @type char
            %  varargin: variable input parameters @type varargin
            %   Vars: a list of symbolic variables @type SymVariable
            %   File: the (full) file name of exported file @type char
            %   ForceExport: force the export @type logical
            %   BuildMex: flag whether to MEX the exported file @type logical
            %   Namespace: the namespace of the function @type char

            export(obj.h_,export_path, varargin{:});

        end

        function nlp = imposeNLPConstraint(obj, nlp)
            % impose holonomic objaints as NLP objaints in the trajectory
            % optimization problem 'nlp' of the dynamical system
            %
            %
            % Parameters:
            % nlp: the trajectory optimization NLP @type TrajectoryOptimization

            if ~isempty(obj.AuxData)
                nlp = addNodeConstraint(nlp, obj.h_, obj.DepLists, 'all',...
                    0, inf, 'Nonlinear', obj.AuxData);
            else
                nlp = addNodeConstraint(nlp, obj.h_, obj.DepLists, 'all',...
                    0, inf, 'Nonlinear');
            end

        end

        function f = calcConstraint(obj, varargin)
            % calculate the unilateral constraints
            %
            % Parameters:
            % varargin: input variable depends on the object
            %
            % Return values:
            % f: the value of the unilateral constraints
            if isempty(obj.AuxData)
                f = feval(obj.h_.Name, varargin{:});
            else
                f = feval(obj.h_.Name, varargin{:}, obj.AuxData{:});
            end
        end

    end

    % set functions
    methods

        % OutputLabel
        function obj = setConstrLabel(obj, label)
            % sets the naming labels of outputs
            %
            % Parameters:
            % label: the cell array of labels @type cellstr

            validateattributes(label,{'cell'},...
                {'nonempty','numel',obj.Dimension,'row'},...
                'UnilateralConstraint','ConstrLabel');
            cellfun( @(x) validateattributes(...
                x, {'char'},{}), label);

            obj.ConstrLabel = label;
        end

        % Name
        function name = validateName(~, name)
            validateattributes(name, {'char'},...
                {'nonempty','scalartext'},...
                'UnilateralConstraint','Name');

            assert(isempty(regexp(name, '\W', 'once')) || ~isempty(regexp(name, '\$', 'once')),...
                'UnilateralConstraint:invalidSymbol', ...
                'Invalid symbol string, can NOT contain special characters.');

            assert(isempty(regexp(name, '_', 'once')),...
                'UnilateralConstraint:invalidSymbol', ...
                'Invalid symbol string, can NOT contain ''_''.');

        end
    end
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [Xtrain, Ytrain, Xtest, Ytest, imageDim] = slaveLoadCIFAR(cifarFileDir, numSlaves, slaveIndex)

  % params
  samplesPerFile = 10000;
  numTrainFiles=5;
  numTestFiles=1;

  % compute numbers of samples
  trainSamplesPerSlave = samplesPerFile * numTrainFiles / numSlaves;
  trainIndex = 1:samplesPerFile*numTrainFiles;
  testSamplesPerSlave = samplesPerFile * numTestFiles / numSlaves;
  testIndex = 1:samplesPerFile*numTestFiles;

  % choose indices for this slave
  firstTrainSample = floor(trainSamplesPerSlave * (slaveIndex-1))+1;
  lastTrainSample = floor(trainSamplesPerSlave * (slaveIndex));
  myTrainIndices = trainIndex(firstTrainSample:lastTrainSample);

  firstTestSample = floor(testSamplesPerSlave * (slaveIndex-1))+1;
  lastTestSample = floor(testSamplesPerSlave * (slaveIndex));
  myTestIndices = testIndex(firstTestSample:lastTestSample);

  % load training data
  Xtrain = zeros(length(myTrainIndices), 3072);
  Ytrain = zeros(length(myTrainIndices), 1);
  for i=1:numTrainFiles
    firstFileSample = samplesPerFile*(i-1)+1;
    fileIndices = firstFileSample : samplesPerFile*i;
    myFileIndices = intersect(fileIndices, myTrainIndices);

    if (length(myFileIndices) > 0)
      fname = sprintf('%s/data_batch_%d.mat', cifarFileDir, i);
      f = load(fname);
      Xtrain(myFileIndices - firstTrainSample + 1, :) = f.data(myFileIndices - firstFileSample + 1, :);
      Ytrain(myFileIndices - firstTrainSample + 1) = double(f.labels(myFileIndices - firstFileSample + 1)) + 1;
    end
  end

  % load test data
  fname = sprintf('%s/test_batch.mat', cifarFileDir);
  f = load(fname);
  Xtest = double(f.data(myTestIndices, :));
  Ytest = double(f.labels(myTestIndices)) + 1;

  imageDim = [32 32 3];

  Xtrain=slaveRef(Xtrain);
  Ytrain=slaveRef(Ytrain);
  Xtest=slaveRef(Xtest);
  Ytest=slaveRef(Ytest);

author: huashiyiqike

repo: LSTM-MATLAB


text:
function varargout = odrfr(varargin)
% ODRFR M-file for odrfr.fig
%      ODRFR, by itself, creates a new ODRFR or raises the existing
%      singleton*.
%
%      H = ODRFR returns the handle to a new ODRFR or the handle to
%      the existing singleton*.
%
%      ODRFR('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in ODRFR.M with the given input arguments.
%
%      ODRFR('Property','Value',...) creates a new ODRFR or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before odrfr_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to odrfr_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help odrfr

% Last Modified by GUIDE v2.5 10-Oct-2004 00:15:21

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @odrfr_OpeningFcn, ...
                   'gui_OutputFcn',  @odrfr_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before odrfr is made visible.
function odrfr_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to odrfr (see VARARGIN)

% Choose default command line output for odrfr
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes odrfr wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = odrfr_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[filename, pathname] = uigetfile({'*.bmp';'*.jpg';'*.gif';'*.*'}, 'Pick an Image File');
S = imread([pathname,filename]);
axes(handles.axes1);
imshow(S);

handles.S = S;
guidata(hObject, handles);

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

S = handles.S;
axes(handles.axes1);
img_crop = imcrop(S);
axes(handles.axes2);
imshow(img_crop);
handles.img_crop = img_crop;

axes(handles.axes7);
I=S;
Igray = rgb2gray(I);
Ibw = im2bw(Igray,graythresh(Igray));
Iedge = edge(uint8(Ibw));
se = strel('square',2);
Iedge2 = imdilate(Iedge, se);
Ifill= imfill(Iedge2,'holes');
[Ilabel num] = bwlabel(Ifill);
Iprops = regionprops(Ilabel);
Ibox = [Iprops.BoundingBox];
[y,x]=size(Ibox);%
x=x/4;%
Ibox = reshape(Ibox,[4 x]);%

handles.Ibox=Ibox;
imshow(I)

selected_col = get(handles.edit5,'string');
selected_col = evalin('base',selected_col);

selected_ln = get(handles.edit6,'string');
selected_ln = evalin('base',selected_ln);

for cnt = 1:selected_ln * selected_col
    rectangle('position',Ibox(:,cnt),'edgecolor','r');
end
guidata(hObject, handles);
% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

img_crop = handles.img_crop;
imgGray = rgb2gray(img_crop);
bw = im2bw(img_crop,graythresh(imgGray));
axes(handles.axes3);
imshow(bw);
bw2 = charcrop(bw);
axes(handles.axes4);
imshow(bw2);
handles.bw2 = bw2;

guidata(hObject, handles);
% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

bw2  = handles.bw2;
charvec = figresize(bw2);
axes(handles.axes5);
plotchar(charvec);

handles.charvec = charvec;
guidata(hObject, handles);
% --- Executes on button press in pushbutton5.
function pushbutton5_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
net = handles.net ;
charvec = handles.charvec;
%selected_net = get(handles.edit1,'string');

result = sim(net,charvec);
[val, num] = max(result);
set(handles.edit2, 'string',num);
guidata(hObject, handles);

function edit1_Callback(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double
% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function edit2_Callback(hObject, eventdata, handles)
% hObject    handle to edit2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit2 as text
%        str2double(get(hObject,'String')) returns contents of edit2 as a double

% --- Executes during object creation, after setting all properties.
function edit2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in pushbutton6.
function pushbutton6_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[filename, pathname] = uigetfile({'*.bmp';'*.jpg';'*.gif';'*.*'}, 'Pick an Image File');
trimg = imread([pathname,filename]);

selected_col = get(handles.edit5,'string');
selected_col = evalin('base',selected_col);

selected_ln = get(handles.edit6,'string');
selected_ln = evalin('base',selected_ln);

img = preprocess(trimg,selected_col,selected_ln);

for cnt = 1:selected_ln * selected_col;
    bw2 = charcrop(img{cnt});
    charvec = figresize(bw2);
    out(:,cnt) = charvec;
end
P = out(:,1:40);
T = [eye(10) eye(10) eye(10) eye(10)];
%Ptest = out(:,(selected_ln -2)* selected_col :selected_ln * selected_col);

%% Creating and training of the Neural Network
net = createnn(P,T);
handles.net = net;

%% Testing the Neural Network
%[a,b]=max(sim(net,Ptest));
%disp(b);
assignin('base','net',net);
guidata(hObject, handles);
function edit5_Callback(hObject, eventdata, handles)
% hObject    handle to edit5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit5 as text
%        str2double(get(hObject,'String')) returns contents of edit5 as a double

% --- Executes during object creation, after setting all properties.
function edit5_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function edit6_Callback(hObject, eventdata, handles)
% hObject    handle to edit6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit6 as text
%        str2double(get(hObject,'String')) returns contents of edit6 as a double

% --- Executes during object creation, after setting all properties.
function edit6_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


text:
eeg_file= fullfile(BTB.DataDir, 'demoMat', 'VPiac_10_10_13', ...
                   'calibration_CenterSpellerMVEP_VPiac');

% Load data
try
  [cnt, mrk, mnt] = file_loadMatlab(eeg_file);
catch
  error('You need to run ''demo_convert_ERPSpeller'' first');
end

nBlocks= 6;

% define 'nBlocks' intervals ranging from the first to the last marker:
inter = linspace( mrk.time(1), mrk.time(end), nBlocks+1 );
blk = struct( 'ival', [inter(1:end-1)' inter(2:end)'] );

% define marker structure: one marker every second
mrkblk = mrk_evenlyInBlocks( blk, 1000 );
% assign class labels for each marker corresponding to the block
mrkblk.y = util_classind2labels( mrkblk.event.blkno );
mrkblk.className= str_cprintf('block %d', 1:6);

% calculate spectra (average within the defined blocks)
spec = proc_segmentation( cnt, mrkblk, [0 1000] );
spec = proc_spectrum( spec, [5 40], kaiser(cnt.fs,2) );

% visualize the block-wise spectra
fig_set(1);
cmap= cmap_hsvFade(nBlocks, [0 5/6], 1, 1);
grid_plot( spec, mnt , defopt_spec, 'ColorOrder',cmap );

% increasing noise level can be observed, e.g. in channels C6 and T8


text:
%% Projection onto the intersection of two sets
% Our goal is to find a point in the intersection of two convex sets C1 and C2.
% A similar goal is to find the projection of an arbitrary point onto the
% intersection (or if the intersection is empty, find the point(s) in C1
% that are as close as possible to C2). In this demo, the intersection
% of C1 and C2 is a singleton, so the concepts are the same.
%
% We assume that we know how to project on the sets C1 and C2 individually,
% but not onto C1 intersect C2.
%
% The basic algorithm to solve this problem is the alternating projection method,
% first studied by John von Neumann for the case where C1 and C2 are affine
% spaces. This algorithm can be extended to arbitrary convex sets, although
% you may not converge to the *projection* of the original point.
%
% This algorithm is very simple: starting at some point y, and with x <-- y,
% we update
%   x <-- Proj_{C1}( Proj_{C2}( x ) )
% where Proj_{C1} is the projector onto the set C1.
%
% Better methods (which actually give you the projection of y onto the intersection)
% are based on Dykstra's algorithm (not to be confused with the Dijkstra algorithm
% described here: http://en.wikipedia.org/wiki/Dijkstra%27s_algorithm ).
%
% Dykstra's algorithm also uses only Proj_{C1} and Proj_{C2}, but with a few
% intermediate steps. For an overview, see the 2011 book by Raydan:
% http://www.amazon.com/Alternating-Projection-Methods-Fundamentals-Algorithms/dp/1611971934
%
% Another overview of both alternating projection and Dykstra are in the book chapter
% "Proximal splitting methods in signal processing" by P. L. Combettes and J.-C. Pesquet,
% in the book 'Fixed-Point Algorithms for Inverse Problems in Science and Engineering',
% (ed.: H. H. Bauschke, R. S. Burachik, P. L. Combettes, V. Elser, D. R. Luke, and
% H. Wolkowicz, Editors), pp. 185-212. Springer, New York, 2011.
% The chapter can be downloaded at http://www.ann.jussieu.fr/~plc/prox.pdf
%
% The purpose of this demo is to show the above two methods, and also
% how to formulate this with TFOCS. The advantage of solving this with TFOCS
% is that we can use an accelerated solver (Nesterov-style) and use large
% step-sizes and line search (whereas Dykstra has no step-size parameter).
% In the first example, it is especially apparent that TFOCS avoids
% the slow convergence that both the alternating projection and Dykstra
% algorithms suffer from. In the second example, the performance
% of alternating projections and TFOCS are quite similar.

%% Mathematical formulation
% For a given point y, and two convex sets C1 and C2, we wish to find
% the closest point x to y, such that x is in both C1 and C2. We write this
% as:
%
% $$ \textrm{minimize}_{x\in C_1 \cap C_2} \,\, \mu/2||x\textrm{--}y||_2^2 $$
%
% The parameter $\mu > 0$ is arbitrary and does not affect the answer.
%
% This fits naturally into the TFOCS formulation since the primal
% problem is already strongly convex.

%% Demo with two polygon sets in 2D that intersect at (0,0)
% Set the dimension to 2
N = 2;  % 2D is best for visualization

%%
% Define some 2D polygons

mx = 1;
x1a = [1;.1];
x1b = [1;.2];
xx1 = [0,x1a(1),2*mx,x1b(1)];
yy1 = [0,x1a(2),2*mx,x1b(2)];

x2a = [1;.3];
x2b = [1;.35];
xx2 = [0,x2a(1),2*mx,x2b(1)];
yy2 = [0,x2a(2),2*mx,x2b(2)];

%%
% The solution is at x = (0,0), so our error function is simple

err = @(x) norm(x);
%%
% plot the regions

figure
red = [255,153,153]/255;
blue = [204,204,255]/255;
fill( xx1, yy1,red);
hold on
fill( xx2, yy2,blue);

xlim( [0,mx] ); ylim( [0,.5*mx] );
axis equal
text(.5,.07,'set C1');
text(.5,.22,'set C2');

%%
% Make some operators that will be used with TFOCS

addpath ~/Dropbox/TFOCS/  % modify this to wherever it is installed on your computer
op1     = project2DCone(x1a, x1b );
op2     = project2DCone(x2a, x2b );

offset1 = 0;
offset2 = 0;
op1_d   = prox_dualize(op1);
op2_d   = prox_dualize(op2);

% and some simpler operators that we will use with alternating projections
%   and with Dykstra
proj1   = @(x) callandmap( op1, 2, x, 1); % gamma is irrelevant
proj2   = @(x) callandmap( op2, 2, x, 1 ); % gamma is irrelevant

% for all methods, we need to specify a starting point

x0      = [1; 1/4];

%% solve in TFOCS
% We can pick any mu and should get the same result, although due to
% some scaling issues in stopping criteria, it does have a small effect.
mu = 1e-6;

opts    = struct('debug',false,'printEvery',20,'maxIts',200);
opts.errFcn{1}  = @(f,d,x) err(x);
opts.errFcn{2}  = @(f,d,x) recordPoints(x);
recordPoints(); % zero-out counter. We use this to plot the points later

opts.tol        = 1e-15;

affineOperator  = {eye(N),offset1;eye(N),offset2};
dualProx        = {op1_d,op2_d};
% Solve in TFOCS:
[x,out,optsOut] = tfocs_SCD( [], affineOperator, dualProx , mu, x0, [], opts );

% Record the path:
path    = recordPoints();
path    = [x0,path];

figure
subplot(1,2,1);
fill( xx1, yy1,red);
hold on
fill( xx2, yy2,blue);
xlim( [0,mx] ); ylim( [0,.5*mx] );
text(.5,.07,'set C1');
text(.5,.22,'set C2');
plot( path(1,:), path(2,:),'ko-','linewidth',2 )
title('Path for TFOCS solving the intersection of 2 polygons');

subplot(1,2,2);
semilogy( out.err(:,1) ,'o-'); xlabel('iterations'); ylabel('error');
title('Error for TFOCS method');
set(gcf, 'Position', [400 200 800 400]);
%% solve via alternating projection method
maxIts  = 500;
x       = x0;
path    = [x0];
errHist = [];
for k = 1:maxIts
    % basic alternating projection method:
    x   = proj1(x);
    path= [path, x];
    x   = proj2(x);
    path= [path, x];
    errHist     = [ errHist; err(x) ];
    if ~mod(k,50)
        fprintf('Iter %4d, error is %.2e\n', k, errHist(end) );
    end
end
figure
subplot(1,2,1);
fill( xx1, yy1,red);
hold on
fill( xx2, yy2,blue);
xlim( [0,mx] ); ylim( [0,.5*mx] );
text(.5,.07,'set C1');
text(.5,.22,'set C2');
plot( path(1,:), path(2,:),'ko-','linewidth',1 )
title('Path for alternating projection method solving the intersection of 2 polygons');

subplot(1,2,2);
semilogy( errHist, 'o-' ); xlabel('iterations'); ylabel('error');
errHist_1 = errHist;
title('Error for alternating projection method');
set(gcf, 'Position', [400 200 800 400]);
%% solve via Dykstra's algorithm

maxIts  = 500;
x       = x0;
[p,q]   = deal( 0*x0 );
p       = -.25*[1;1];
path    = [x0];
errHist = [];
for k = 1:maxIts
    % If x + p is feasible, the y = (x+p), so p=x+p-y = 0.
    y   = proj1( x + p );
    p   = x + p - y;
    x   = proj2( y + q );
    q   = y + q - x;
    path = [path,y,x];
    errHist     = [ errHist; err(x) ];
    if ~mod(k,50)
        fprintf('Iter %4d, error is %.2e\n', k, errHist(end) );
    end
end
figure
subplot(1,2,1);
fill( xx1, yy1,red);
hold on
fill( xx2, yy2,blue);
xlim( [0,mx] ); ylim( [0,.5*mx] );
text(.5,.07,'set C1');
text(.5,.22,'set C2');
plot( path(1,:), path(2,:),'ko-' )
title('Path for Dykstra''s algo solving the intersection of 2 polygons');

subplot(1,2,2);
semilogy( errHist, 'o-' ); xlabel('iterations'); ylabel('error');
title('Error for Dykstra''s method');
set(gcf, 'Position', [400 200 800 400]);

%%
% With Dykstra's algo (and alternating projections), we get
% very slow convergence, because of the very low angle between the
% two sets. Here is a zoom in on the graph of the iterates:
figure
fill( xx1, yy1,red);
hold on
fill( xx2, yy2,blue);
xlim( [0,mx] ); ylim( [0,.5*mx] );
text(.5,.07,'set C1');
text(.5,.22,'set C2');
plot( path(1,:), path(2,:),'ko-' )
title('Path for Dykstra''s algo solving the intersection of 2 polygons (zoom)');
xlim( [.25,.4] );
ylim( [.058,.1])
%% Plot all the errors together
figure
semilogy( out.err(:,1),'-' ,'linewidth',3);
hold all
semilogy( errHist_1,'-','linewidth',3);
semilogy( errHist,'--','linewidth',3);
legend('TFOCS','Alternating projection','Dykstra');
xlabel('iteration'); ylabel('error');
%% Demo with two circles that intersect at (0,0)
center1     = [0;.5];
radius1     = .5;
center2     = -center1;
radius2     = radius1;

figure
rectangle('Position',[-radius1,0,2*radius1,2*radius1],'Curvature',[1,1],'FaceColor',red)
rectangle('Position',[-radius1,-2*radius2,2*radius1,2*radius1],'Curvature',[1,1],'FaceColor',blue)
axis equal
hold on
text(0,.5,'set C1');
text(0,-.5,'set C2');

offset1     = center1;
offset2     = center2;

op1_d       = prox_l2(radius1);
op2_d       = prox_l2(radius2);
x0          = [1; .2];

% and some simpler operators that we will use with alternating projections
%   and with Dykstra
proj1   = @(x) radius1*(x-center1)/norm(x-center1) + center1;
proj2   = @(x) radius2*(x-center2)/norm(x-center2) + center2;
%% solve in TFOCS
opts.maxIts     = 300;
affineOperator  = {eye(N),offset1;eye(N),offset2};
dualProx        = {op1_d,op2_d};
% Solve in TFOCS:
[x,out,optsOut] = tfocs_SCD( [], affineOperator, dualProx , mu, x0, [], opts );

% Record the path:
path    = recordPoints();
path    = [x0,path];
% Plot:
figure
subplot(1,3,1);
rectangle('Position',[-radius1,0,2*radius1,2*radius1],'Curvature',[1,1],'FaceColor',red)
rectangle('Position',[-radius1,-2*radius2,2*radius1,2*radius1],'Curvature',[1,1],'FaceColor',blue)
hold on
text(0,.5,'set C1');
text(0,-.5,'set C2');
plot( path(1,:), path(2,:),'ko-','linewidth',2 )
title('Path for TFOCS solving the intersection of 2 circles');

subplot(1,3,2);
rectangle('Position',[-radius1,0,2*radius1,2*radius1],'Curvature',[1,1],'FaceColor',red)
rectangle('Position',[-radius1,-2*radius2,2*radius1,2*radius1],'Curvature',[1,1],'FaceColor',blue)
hold on
plot( path(1,:), path(2,:),'ko-','linewidth',2 )
xlim([0,.15]);
ylim([-.1,.1]);
title('(zoom)');

subplot(1,3,3);
semilogy( out.err,'o-')
title('Error of iterates for TFOCS method');
set(gcf, 'Position', [400 200 800 400]);
%% solve via alternating projection method
maxIts  = 300;
x       = x0;
path    = [x0];
errHist = [];
for k = 1:maxIts
    % basic alternating projection method:
    x   = proj1(x);
    path= [path, x];
    x   = proj2(x);
    path= [path, x];
    errHist     = [ errHist; err(x) ];
    if ~mod(k,50)
        fprintf('Iter %4d, error is %.2e\n', k, errHist(end) );
    end
end
figure
subplot(1,3,1);
rectangle('Position',[-radius1,0,2*radius1,2*radius1],'Curvature',[1,1],'FaceColor',red)
rectangle('Position',[-radius1,-2*radius2,2*radius1,2*radius1],'Curvature',[1,1],'FaceColor',blue)
hold on
text(0,.5,'set C1');
text(0,-.5,'set C2');
plot( path(1,:), path(2,:),'ko-' )
title('Path for alternating projection method solving the intersection of 2 circles');

subplot(1,3,2);
rectangle('Position',[-radius1,0,2*radius1,2*radius1],'Curvature',[1,1],'FaceColor',red)
rectangle('Position',[-radius1,-2*radius2,2*radius1,2*radius1],'Curvature',[1,1],'FaceColor',blue)
hold on
plot( path(1,:), path(2,:),'ko-','linewidth',2 )
xlim([0,.15]);
ylim([-.1,.1]);
title('(zoom)');

subplot(1,3,3);
semilogy( errHist,'o-')
errHist_1 = errHist;
title('Error of iterates for alternating projection method');
set(gcf, 'Position', [400 200 800 400]);
%% solve via Dykstra's algorithm

maxIts  = 300;
x       = x0;
[p,q]   = deal( 0*x0 );
path    = [x0];
errHist = [];
for k = 1:maxIts
    % If x + p is feasible, the y = (x+p), so p=x+p-y = 0.
    y   = proj1( x + p );
    p   = x + p - y;
    x   = proj2( y + q );
    q   = y + q - x;
    path = [path,y];
    path = [path,x];
    errHist     = [ errHist; err(x) ];
    if ~mod(k,50)
        fprintf('Iter %4d, error is %.2e\n', k, errHist(end) );
    end
end
figure
subplot(1,3,1);
rectangle('Position',[-radius1,0,2*radius1,2*radius1],'Curvature',[1,1],'FaceColor',red)
rectangle('Position',[-radius1,-2*radius2,2*radius1,2*radius1],'Curvature',[1,1],'FaceColor',blue)
hold on
text(0,.5,'set C1');
text(0,-.5,'set C2');
plot( path(1,:), path(2,:),'ko-' )
title('Path for Dykstra''s algo solving the intersection of 2 circles');

subplot(1,3,2);
rectangle('Position',[-radius1,0,2*radius1,2*radius1],'Curvature',[1,1],'FaceColor',red)
rectangle('Position',[-radius1,-2*radius2,2*radius1,2*radius1],'Curvature',[1,1],'FaceColor',blue)
hold on
plot( path(1,:), path(2,:),'ko-','linewidth',2 )
xlim([0,.15]);
ylim([-.1,.1]);
title('(zoom)');

subplot(1,3,3);
semilogy( errHist,'o-')
title('Error of iterates for Dykstra''s algo');
set(gcf, 'Position', [400 200 800 400]);

%% Plot all the errors together
figure
semilogy( out.err(:,1),'-' ,'linewidth',3);
hold all
semilogy( errHist_1,'-','linewidth',3);
semilogy( errHist,'--','linewidth',3);
legend('TFOCS','Alternating projection','Dykstra');
xlabel('iteration'); ylabel('error');

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% MRI data of a human knee, image courtesy by
% Thomas Netsch, Philips Medical Solutions, Hamburg
%==============================================================================

checkDataFile
if expfileExists, return; end;

infile = fullfile(FAIRpath,'kernel','data','knee3D.mat');
try
  load(infile);
catch
  fprintf('sorry, this data is not available');
  expfile = [];
  return;
end;

dataT = double(dataT);
dataR = double(dataR);
viewK = @(T) volView(T,omega,m,'isovalue',15,'view',[-20,-5],'colormap','bone(256)');

FAIRfigure(1,'color','w','figname',sprintf('%s/template',mfilename)); clf;
viewK(dataT); hold on; axis off; colormap(gray(128));
FAIRfigure(2,'color','w','figname',sprintf('%s/reference',mfilename)); clf;
viewK(dataR); hold on; axis off; colormap(gray(128));

[viewer,viewPara] = viewImage('reset','viewImage','imgmontage',...
  'direction','-zyx','colormap','bone(256)');
imgPara     = {'inter','linearInterMex'};
traPara     = {'trafo','affine3Dsparse'};
disPara     = {'distance','SSD'};
regPara     = {'regularizer','mfElastic','alpha',500,'mu',1,'lambda',0};

ML = getMultilevel({dataT,dataR},omega,m,'fig',2);

% save to outfile
save(expfile,'dataT','dataR','omega','m','ML',...
  'viewPara','imgPara','traPara','disPara','regPara');

checkDataFile
%==============================================================================


text:
function [fusedBox] = fusePETCT(ROIbox_PET,ROIbox_CT,maskBox_PET,CTinv,CTweight,wavelet)
% -------------------------------------------------------------------------
% function [fusedBox] = fusePETCT(ROIbox_PET,ROIbox_CT,maskBox_PET,CTinv,CTweight,wavelet)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function fuses the region of interest (ROI) of two registered PET and
% CT volumes using a technique based on the wavelet transform. See Ref. [2]
% for more details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% -------------------------------------------------------------------------
% INPUTS:
% - ROIbox_PET: 3D array of the smallest box containing the ROI of the PET
%               volume.
% - ROIbox_CT: 3D array of the smallest box containing the ROI of the CT
%               volume.
% - maskBox_PET: 3D array of the mask of the smallest box specifying the
%                ROI of the PET volume. Voxels within the ROI are assigned
%                value of 1, voxels outside a value of 0.
% - CTinv: String specifying if the intensities of the CT volume are
%          inverted prior to fusion with PET. Either 'Inv' for inversion,
%          or 'NoInv' for no inversion.
% - CTweight: Numerical value specifying the weight of the CT scan in the
%              fusion with PET.
% - wavelet: String specifying the name of the MATLAB wavelet basis used
%            wavelet basis used in the fusion process.
% -------------------------------------------------------------------------
% OUTPUTS:
% - fusedBox: 3D array of the smallest box containing the ROI of the fused
%             PET/CT volume.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: July 2015
% -------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

% PET PRE_PROCESSING
ROIbox_PET = sqrt(ROIbox_PET);

% RESAMPLING CT VOLUME TO PET IN-PLANE RESOLUTION
szPET = size(ROIbox_PET);
temp=zeros(szPET);
for k = 1:szPET(3)
    temp(:,:,k) = imresize(ROIbox_CT(:,:,k),[szPET(1),szPET(2)],'Method','cubic','Antialiasing',true);
end
ROIbox_CT = temp;

% NORMALIZATION (and inversion) OF VOLUMES
ROIOnly_CT = ROIbox_CT;
ROIOnly_PET = ROIbox_PET;
ROIOnly_CT(maskBox_PET==0) = NaN;
ROIOnly_PET(maskBox_PET==0) = NaN;
minCT = min(ROIOnly_CT(:));
ROIbox_CT = ROIbox_CT - minCT;
ROIOnly_CT = ROIOnly_CT - minCT;
ROIbox_CT = ROIbox_CT./max(ROIOnly_CT(:)).*255;
if strcmp(CTinv,'Inv')
    ROIbox_CT = 255 - ROIbox_CT;
end
minPET = min(ROIOnly_PET(:));
ROIbox_PET = ROIbox_PET - minPET;
ROIOnly_PET = ROIOnly_PET - minPET;
ROIbox_PET = ROIbox_PET./max(ROIOnly_PET(:)).*255;

% WAVELET DECOMPOSITION AND FUSION OF COEFFICIENTS
wdecCT = wavedec3(ROIbox_CT,1,wavelet);
wdecPET = wavedec3(ROIbox_PET,1,wavelet);
wdecFUSED = wdecPET;
szDEC = size(wdecFUSED.dec{1}); % All sub-bands are of the same size

% Fusion of the LLL sub-bands (dec{1})
b = 1;
wdecFUSED.dec{b} = (1-CTweight).*wdecPET.dec{b} + CTweight.*wdecCT.dec{b};

% Fusion of the HLL sub-bands (dec{2})
b = 2;
tempPET = zeros(szDEC);
tempCT = zeros(szDEC);
for k = 1:szDEC(3)
    [~,Gy] = imgradientxy(wdecPET.dec{b}(:,:,k)); tempPET(:,:,k) = sqrt(Gy.^2);
    [~,Gy] = imgradientxy(wdecCT.dec{b}(:,:,k)); tempCT(:,:,k) = sqrt(Gy.^2);
end
wdecFUSED.dec{b}(tempCT > tempPET) = wdecCT.dec{b}(tempCT > tempPET);

% Fusion of the LHL sub-bands (dec{3})
b = 3;
tempPET = zeros(szDEC);
tempCT = zeros(szDEC);
for k = 1:szDEC(3)
    [Gx,~] = imgradientxy(wdecPET.dec{b}(:,:,k)); tempPET(:,:,k) = sqrt(Gx.^2);
    [Gx,~] = imgradientxy(wdecCT.dec{b}(:,:,k)); tempCT(:,:,k) = sqrt(Gx.^2);
end
wdecFUSED.dec{b}(tempCT > tempPET) = wdecCT.dec{b}(tempCT > tempPET);

% Fusion of the HHL sub-bands (dec{4})
b = 4;
tempPET = zeros(szDEC);
tempCT = zeros(szDEC);
for k = 1:szDEC(3)
    [Gx,Gy] = imgradientxy(wdecPET.dec{b}(:,:,k)); tempPET(:,:,k) = sqrt(Gx.^2 + Gy.^2);
    [Gx,Gy] = imgradientxy(wdecCT.dec{b}(:,:,k)); tempCT(:,:,k) = sqrt(Gx.^2 + Gy.^2);
end
wdecFUSED.dec{b}(tempCT > tempPET) = wdecCT.dec{b}(tempCT > tempPET);

% Fusion of the LLH sub-bands (dec{5})
b = 5;
tempPET = zeros(szDEC);
tempCT = zeros(szDEC);
for i = 1:szDEC(1)
    [~,Gy_PET] = imgradientxy(squeeze(wdecPET.dec{b}(i,:,:))');
    [~,Gy_CT] = imgradientxy(squeeze(wdecCT.dec{b}(i,:,:))');
    for k = 1:szDEC(3)
        tempPET(i,:,k) = sqrt((Gy_PET(k,:)).^2);
        tempCT(i,:,k) = sqrt((Gy_CT(k,:)).^2);
    end
end
wdecFUSED.dec{b}(tempCT > tempPET) = wdecCT.dec{b}(tempCT > tempPET);

% Fusion of the HLH sub-bands (dec{6})
b = 6;
tempPET = zeros(szDEC);
tempCT = zeros(szDEC);
for k = 1:szDEC(3)
    [~,Gy] = imgradientxy(wdecPET.dec{b}(:,:,k)); tempPET(:,:,k) = Gy.^2;
    [~,Gy] = imgradientxy(wdecCT.dec{b}(:,:,k)); tempCT(:,:,k) = Gy.^2;
end
for i = 1:szDEC(1)
    [~,Gy_PET] = imgradientxy(squeeze(wdecPET.dec{b}(i,:,:))');
    [~,Gy_CT] = imgradientxy(squeeze(wdecCT.dec{b}(i,:,:))');
    for k = 1:szDEC(3)
        tempPET(i,:,k) = tempPET(i,:,k) + (Gy_PET(k,:)).^2;
        tempCT(i,:,k) = tempCT(i,:,k) + (Gy_CT(k,:)).^2;
    end
end
tempPET = sqrt(tempPET); tempCT = sqrt(tempCT);
wdecFUSED.dec{b}(tempCT > tempPET) = wdecCT.dec{b}(tempCT > tempPET);

% Fusion of the LHH sub-bands (dec{7})
b = 7;
tempPET = zeros(szDEC);
tempCT = zeros(szDEC);
for k = 1:szDEC(3)
    [Gx,~] = imgradientxy(wdecPET.dec{b}(:,:,k)); tempPET(:,:,k) = Gx.^2;
    [Gx,~] = imgradientxy(wdecCT.dec{b}(:,:,k)); tempCT(:,:,k) = Gx.^2;
end
for i = 1:szDEC(1)
    [~,Gy_PET] = imgradientxy(squeeze(wdecPET.dec{b}(i,:,:))');
    [~,Gy_CT] = imgradientxy(squeeze(wdecCT.dec{b}(i,:,:))');
    for k = 1:szDEC(3)
        tempPET(i,:,k) = tempPET(i,:,k) + (Gy_PET(k,:)).^2;
        tempCT(i,:,k) = tempCT(i,:,k) + (Gy_CT(k,:)).^2;
    end
end
tempPET = sqrt(tempPET); tempCT = sqrt(tempCT);
wdecFUSED.dec{b}(tempCT > tempPET) = wdecCT.dec{b}(tempCT > tempPET);

% Fusion of the HHH sub-bands (dec{8})
b = 8;
tempPET = zeros(szDEC);
tempCT = zeros(szDEC);
for k = 1:szDEC(3)
    [Gx,Gy] = imgradientxy(wdecPET.dec{b}(:,:,k)); tempPET(:,:,k) = Gx.^2 + Gy.^2;
    [Gx,Gy] = imgradientxy(wdecCT.dec{b}(:,:,k)); tempCT(:,:,k) = Gx.^2 + Gy.^2;
end
for i = 1:szDEC(1)
    [~,Gy_PET] = imgradientxy(squeeze(wdecPET.dec{b}(i,:,:))');
    [~,Gy_CT] = imgradientxy(squeeze(wdecCT.dec{b}(i,:,:))');
    for k = 1:szDEC(3)
        tempPET(i,:,k) = tempPET(i,:,k) + (Gy_PET(k,:)).^2;
        tempCT(i,:,k) = tempCT(i,:,k) + (Gy_CT(k,:)).^2;
    end
end
tempPET = sqrt(tempPET); tempCT = sqrt(tempCT);
wdecFUSED.dec{b}(tempCT > tempPET) = wdecCT.dec{b}(tempCT > tempPET);

% WAVELET RECONSTRUCTION
fusedBox = waverec3(wdecFUSED);
fusedOnly = fusedBox;
fusedOnly(maskBox_PET==0) = NaN;
minFus = min(fusedOnly(:));
fusedBox = fusedBox - minFus;
fusedOnly = fusedOnly - minFus;
fusedBox = fusedBox./max(fusedOnly(:)).*255;

end

author: mvallieres

repo: radiomics


text:
function [ node_num, node_dim ] = triangle_node_size_example ( )

%*****************************************************************************80
%
%% GMSH_MESH2D_NODE_SIZE_EXAMPLE: sizes of node information for the example.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NODE_NUM, the number of nodes.
%
%    Output, integer NODE_DIM, the spatial dimension.
%
  node_num = 21;
  node_dim = 2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function inpainted_depth = inpaint_depth_with_plane(pixel_coordinates, plane)
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here

number_of_pixels = size(pixel_coordinates, 1);

% At each pixel, return the **rectified** interpolation using the plane that is
% assumed to contain all input pixels, in order to always output positive depth
% values.
inpainted_depth = max([pixel_coordinates, ones(number_of_pixels, 1)] * plane,...
    0);

end

author: sakaridis