text:
function planC = translateStruct(structNum,xyzT,structName,planC)
%
% function planC = translateStruct(structNum,xyzT,structName,planC)
%
% This function creates a new structure by translating structNum by amounts
% xT, yT and zT in x,y,z directions respectively. Note that xyzT = [xT yT zT].
% structName must be a string to name the new structure.
%
% Example:
%   planC = translateStruct(12,[-1 0 0],'moveStruct12');
%
% APA, 9/12/06
%
% LM DK compatible with CERR 3.0
%
% See also GETUNIFORMSTR MASKTOCERRSTRUCTURE
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('planC')

    global planC

end

indexS = planC{end};

%obtain scanNum associated to the structNum for CERR3

%scanNum = getAssociatedScan(planC{indexS.structures}(structNum).assocScanUID);

%scanNum is always 1 for CERR2

scanNum =  getStructureAssociatedScan(structNum, planC);

%obtain r,c,s coordinates of scan x,y,z vals

scanS = planC{indexS.scan}(scanNum);

[xUnifScanValsV, yUnifScanValsV, zUnifScanValsV] = getUniformScanXYZVals(scanS);

[rScanValsV, cScanValsV, sScanValsV] = xyztom(xUnifScanValsV, yUnifScanValsV, zUnifScanValsV,scanNum, planC, 1);

%obtain r,c,s coordinates of structure based on its associated scan

rcsStructValsV = getUniformStr(structNum);

scanUnifSiz = getUniformScanSize(scanS);

rcsStructValsV = find(rcsStructValsV);

[rStructValsV, cStructValsV, sStructValsV] = ind2sub(scanUnifSiz,rcsStructValsV);

%obtain x,y,z coordinates of voxels included within the structure

xStructValsV = xUnifScanValsV(cStructValsV);

yStructValsV = yUnifScanValsV(rStructValsV);

zStructValsV = zUnifScanValsV(sStructValsV);

%translate the x,y,z pints included within the structure by specified amount

xStructValsV = xStructValsV + xyzT(1);

yStructValsV = yStructValsV + xyzT(2);

zStructValsV = zStructValsV + xyzT(3);

%convert the translated x,y,z vals of structure to r,c,s of scanNum

[rStructValsV, cStructValsV, sStructValsV] = xyztom(xStructValsV, yStructValsV, zStructValsV,scanNum, planC, 1);

rStructValsV = round(rStructValsV);

rStructValsV = clip(rStructValsV,min(rScanValsV),max(rScanValsV),'limits');

cStructValsV = round(cStructValsV);

cStructValsV = clip(cStructValsV,min(cScanValsV),max(cScanValsV),'limits');

sStructValsV = round(sStructValsV);

sStructValsV = clip(sStructValsV,min(sScanValsV),max(sScanValsV),'limits');

%generate uniformized mask for this new structure

maskM = zeros(scanUnifSiz);

indicesWithinSkinV = sub2ind(scanUnifSiz,rStructValsV,cStructValsV,sStructValsV);

maskM(indicesWithinSkinV) = 1;

%generate contours on slices out of uniform mask and add to planC

planC = maskToCERRStructure(maskM, 1, scanNum, structName, planC);


text:
% eegplot helper function to read key strokes

function eegplot_readkey(src,evnt)
    if strcmp(evnt.Key, 'rightarrow')==1
        eegplot('drawp',4);
    elseif strcmp(evnt.Key, 'leftarrow')==1
        eegplot('drawp',1);
    end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%-----------------------------------------------------------------------
% Job configuration created by cfg_util (rev $Rev: 3355 $)
%-----------------------------------------------------------------------
matlabbatch{1}.cfg_basicio.cfg_named_input.name = 'a';
matlabbatch{1}.cfg_basicio.cfg_named_input.input = '<UNDEFINED>';
matlabbatch{2}.cfg_basicio.cfg_named_input.name = 'b';
matlabbatch{2}.cfg_basicio.cfg_named_input.input = '<UNDEFINED>';
matlabbatch{3}.cfg_toy{1}.add2{1}.cfg_example_add1.a(1) = cfg_dep;
matlabbatch{3}.cfg_toy{1}.add2{1}.cfg_example_add1.a(1).tname = 'Input a';
matlabbatch{3}.cfg_toy{1}.add2{1}.cfg_example_add1.a(1).tgt_spec{1}.name = 'class';
matlabbatch{3}.cfg_toy{1}.add2{1}.cfg_example_add1.a(1).tgt_spec{1}.value = 'cfg_entry';
matlabbatch{3}.cfg_toy{1}.add2{1}.cfg_example_add1.a(1).sname = 'Named Input: a';
matlabbatch{3}.cfg_toy{1}.add2{1}.cfg_example_add1.a(1).src_exbranch = substruct('.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{3}.cfg_toy{1}.add2{1}.cfg_example_add1.a(1).src_output = substruct('.','input');
matlabbatch{3}.cfg_toy{1}.add2{1}.cfg_example_add1.b(1) = cfg_dep;
matlabbatch{3}.cfg_toy{1}.add2{1}.cfg_example_add1.b(1).tname = 'Input b';
matlabbatch{3}.cfg_toy{1}.add2{1}.cfg_example_add1.b(1).tgt_spec{1}.name = 'class';
matlabbatch{3}.cfg_toy{1}.add2{1}.cfg_example_add1.b(1).tgt_spec{1}.value = 'cfg_entry';
matlabbatch{3}.cfg_toy{1}.add2{1}.cfg_example_add1.b(1).sname = 'Named Input: b';
matlabbatch{3}.cfg_toy{1}.add2{1}.cfg_example_add1.b(1).src_exbranch = substruct('.','val', '{}',{2}, '.','val', '{}',{1});
matlabbatch{3}.cfg_toy{1}.add2{1}.cfg_example_add1.b(1).src_output = substruct('.','input');
matlabbatch{4}.cfg_toy{1}.cfg_example_div.a(1) = cfg_dep;
matlabbatch{4}.cfg_toy{1}.cfg_example_div.a(1).tname = 'Input a';
matlabbatch{4}.cfg_toy{1}.cfg_example_div.a(1).tgt_spec = {};
matlabbatch{4}.cfg_toy{1}.cfg_example_div.a(1).sname = 'Add1: Add1: a + b';
matlabbatch{4}.cfg_toy{1}.cfg_example_div.a(1).src_exbranch = substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{4}.cfg_toy{1}.cfg_example_div.a(1).src_output = substruct('()',{1});
matlabbatch{4}.cfg_toy{1}.cfg_example_div.b(1) = cfg_dep;
matlabbatch{4}.cfg_toy{1}.cfg_example_div.b(1).tname = 'Input b';
matlabbatch{4}.cfg_toy{1}.cfg_example_div.b(1).tgt_spec{1}.name = 'class';
matlabbatch{4}.cfg_toy{1}.cfg_example_div.b(1).tgt_spec{1}.value = 'cfg_entry';
matlabbatch{4}.cfg_toy{1}.cfg_example_div.b(1).sname = 'Named Input: b';
matlabbatch{4}.cfg_toy{1}.cfg_example_div.b(1).src_exbranch = substruct('.','val', '{}',{2}, '.','val', '{}',{1});
matlabbatch{4}.cfg_toy{1}.cfg_example_div.b(1).src_output = substruct('.','input');
matlabbatch{5}.cfg_basicio.cfg_assignin.name = '<UNDEFINED>';
matlabbatch{5}.cfg_basicio.cfg_assignin.output(1) = cfg_dep;
matlabbatch{5}.cfg_basicio.cfg_assignin.output(1).tname = 'Output Item';
matlabbatch{5}.cfg_basicio.cfg_assignin.output(1).tgt_spec = {};
matlabbatch{5}.cfg_basicio.cfg_assignin.output(1).sname = 'div: a div b: mod';
matlabbatch{5}.cfg_basicio.cfg_assignin.output(1).src_exbranch = substruct('.','val', '{}',{4}, '.','val', '{}',{1});
matlabbatch{5}.cfg_basicio.cfg_assignin.output(1).src_output = substruct('.','mod');

author: spm


text:
function varargout = AssessmentClassificationResultsViewerUI(varargin)
% ASSESSMENTCLASSIFICATIONRESULTSVIEWERUI MATLAB code for AssessmentClassificationResultsViewerUI.fig
%      ASSESSMENTCLASSIFICATIONRESULTSVIEWERUI, by itself, creates a new ASSESSMENTCLASSIFICATIONRESULTSVIEWERUI or raises the existing
%      singleton*.
%
%      H = ASSESSMENTCLASSIFICATIONRESULTSVIEWERUI returns the handle to a new ASSESSMENTCLASSIFICATIONRESULTSVIEWERUI or the handle to
%      the existing singleton*.
%
%      ASSESSMENTCLASSIFICATIONRESULTSVIEWERUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in ASSESSMENTCLASSIFICATIONRESULTSVIEWERUI.M with the given input arguments.
%
%      ASSESSMENTCLASSIFICATIONRESULTSVIEWERUI('Property','Value',...) creates a new ASSESSMENTCLASSIFICATIONRESULTSVIEWERUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before AssessmentClassificationResultsViewerUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to AssessmentClassificationResultsViewerUI_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 AssessmentClassificationResultsViewerUI

% Last Modified by GUIDE v2.5 11-Nov-2019 17:39:39

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @AssessmentClassificationResultsViewerUI_OpeningFcn, ...
                   'gui_OutputFcn',  @AssessmentClassificationResultsViewerUI_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 AssessmentClassificationResultsViewerUI is made visible.
function AssessmentClassificationResultsViewerUI_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 AssessmentClassificationResultsViewerUI (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes AssessmentClassificationResultsViewerUI wait for user response (see UIRESUME)
% uiwait(handles.mainFigure);

% --- Outputs from this function are returned to the command line.
function varargout = AssessmentClassificationResultsViewerUI_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 selection change in signalsList.
function signalsList_Callback(hObject, eventdata, handles)
% hObject    handle to signalsList (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns signalsList contents as cell array
%        contents{get(hObject,'Value')} returns selected item from signalsList

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

% Hint: listbox 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 selection change in algorithmsList.
function algorithmsList_Callback(hObject, eventdata, handles)
% hObject    handle to algorithmsList (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns algorithmsList contents as cell array
%        contents{get(hObject,'Value')} returns selected item from algorithmsList

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

% Hint: listbox 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 visualizeButton.
function visualizeButton_Callback(hObject, eventdata, handles)
% hObject    handle to visualizeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on selection change in filesList.
function filesList_Callback(hObject, eventdata, handles)
% hObject    handle to filesList (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns filesList contents as cell array
%        contents{get(hObject,'Value')} returns selected item from filesList

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

% Hint: listbox 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 loadDataButton.
function loadDataButton_Callback(hObject, eventdata, handles)
% hObject    handle to loadDataButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function [scanArraySuperior, isCompress, isRemote] = getScanArraySuperior(scanIndex, planC)
%"getScanArray"
%   Returns the scanArray stored in a planC, in slot scanIndex.  The last
%   returned scanArray is cached in a persistent variable to avoid delays
%   in loading remotely stored or compressed scans.
%
%   If a scan array is compressed it is decompressed before being returned.
%
%   If planC is not passed in, the global planC is used.
%
%   If scanIndex is not valid, or if no valid planC can be found an error is
%   returned.
%
%   isCompress and isRemote are booleans to let the calling function know
%   if the scanArray was compressed and/or remote.
%
%   JRA 10/14/04
%
%Usage:
%   function scanArray = getScanArray(scanIndex, planC)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

persistent lastScanStruct;
persistent lastScanArray;
persistent isLastScanCompressed;
persistent isLastScanRemote;

%Check if a scanStruct, not a scanIndex was passed.
if isstruct(scanIndex) && isfield(scanIndex, 'scanArray')

    scanStruct = scanIndex(1);

else
    %An index was passed, extract the scanStruct.
	if ~exist('planC','var')
        global planC
	end

	if ~iscell(planC)
        error('Cannot get scan. PlanC is not a valid cell array.');
	end

	indexS = planC{end};

	if scanIndex > length(planC{indexS.scan}) || scanIndex < 1
        error('Cannot get scan.  Requested scanIndex scan not exist.')
	end

    scanStruct = planC{indexS.scan}(scanIndex);

end

scanArraySuperior = scanStruct.scanArraySuperior;

%If remote or compressed, and struct is same as last time, return cached array.
if (isCompressed(scanArraySuperior) || ...
        ~isLocal(scanArraySuperior)) && isequal(scanStruct, lastScanStruct)
    scanArraySuperior = lastScanArray;
    isCompress = isLastScanCompressed;
    isRemote = isLastScanRemote;
    return;
%If remote or compressed and NOT the same as last time, clear cache.
elseif (isCompressed(scanArraySuperior) || ...
        ~isLocal(scanArraySuperior)) && ~isequal(scanStruct, lastScanStruct)
	lastScanArray           = [];
	lastScanStruct          = [];
	isLastScanCompressed    = [];
	isLastScanRemote        = [];
end

%To be output for calling function's use.
isCompress  = 0;
isRemote    = 0;

%Decompress and follow all file pointers until get to an array.
while isCompressed(scanArraySuperior) || ~isLocal(scanArraySuperior)
    if isCompressed(scanArraySuperior)
        scanArraySuperior   = decompress(scanArraySuperior);
        isCompress  = 1;
    elseif ~isLocal(scanArraySuperior)
        scanArraySuperior   = getRemoteVariable(scanArraySuperior);
        isRemote    = 1;
    end
end

lastScanArray           = scanArraySuperior;
lastScanStruct          = scanStruct;
isLastScanCompressed    = isCompress;
isLastScanRemote        = isRemote;

%temp

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function cdf = stdnormal_cdf (x)
% For each component of x, compute the CDF of the standard normal
% distribution at x.

% Description: CDF of the standard normal distribution

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
% Copyright (C) 1995-2011 Kurt Hornik
% Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       KH <Kurt.Hornik@wu-wien.ac.at>
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

if (nargin ~= 1)
    error('Requires one input arguments.');
end

if (numel(x) == 0)
    error ('stdnormal_cdf: X must not be empty');
end

cdf = erfc (x / (-sqrt(2))) / 2;

author: goGPS-Project


text:
%% Machine Learning Online Class - Exercise 2: Logistic Regression
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the second part
%  of the exercise which covers regularization with logistic regression.
%
%  You will need to complete the following functions in this exericse:
%
%     sigmoid.m
%     costFunction.m
%     predict.m
%     costFunctionReg.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear all; close all; clc

%% Load Data
%  The first two columns contains the exam scores and the third column
%  contains the label.

data = csvread('ex2data2.txt');
X = data(:, [1, 2]); y = data(:, 3);

plotData(X, y);

% Put some labels
hold on;

% Labels and Legend
xlabel('Microchip Test 1')
ylabel('Microchip Test 2')

% Specified in plot order
legend('y = 1', 'y = 0')
hold off;

%% =========== Part 1: Regularized Logistic Regression ============
%  In this part, you are given a dataset with data points that are not
%  linearly separable. However, you would still like to use logistic
%  regression to classify the data points.
%
%  To do so, you introduce more features to use -- in particular, you add
%  polynomial features to our data matrix (similar to polynomial
%  regression).
%

% Add Polynomial Features

% Note that mapFeature also adds a column of ones for us, so the intercept
% term is handled
X = mapFeature(X(:,1), X(:,2));

% Initialize fitting parameters
initial_theta = zeros(size(X, 2), 1);

% Set regularization parameter lambda to 1
lambda = 1;

% Compute and display initial cost and gradient for regularized logistic
% regression
[cost, grad] = costFunctionReg(initial_theta, X, y, lambda);

fprintf('Cost at initial theta (zeros): %f\n', cost);

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============= Part 2: Regularization and Accuracies =============
%  Optional Exercise:
%  In this part, you will get to try different values of lambda and
%  see how regularization affects the decision coundart
%
%  Try the following values of lambda (0, 1, 10, 100).
%
%  How does the decision boundary change when you vary lambda? How does
%  the training set accuracy vary?
%

% Initialize fitting parameters
initial_theta = zeros(size(X, 2), 1);

% Set regularization parameter lambda to 1 (you should vary this)
lambda = 1;

% Set Options
options = optimset('GradObj', 'on', 'MaxIter', 400);

% Optimize
[theta, J, exit_flag] = ...
	fminunc(@(t)(costFunctionReg(t, X, y, lambda)), initial_theta, options);

% Plot Boundary
plotDecisionBoundary(theta, X, y);
hold on;
title(sprintf('lambda = %g', lambda))

% Labels and Legend
xlabel('Microchip Test 1')
ylabel('Microchip Test 2')

legend('y = 1', 'y = 0', 'Decision boundary')
hold off;

% Compute accuracy on our training set
p = predict(theta, X);

fprintf('Train Accuracy: %f\n', mean(double(p == y)) * 100);

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
% The COBRAToolbox: testGeneMCS.m
%
% Purpose:
%     - Check the function of geneMCS using different solvers
%
% Authors:
%     - Luis V. Valcarcel 2017-11-17
%

global CBTDIR
requiredSolvers = {'ibm_cplex'};
prepareTest('requiredSolvers',requiredSolvers);

% save the current path
currentDir = pwd;

% initialize the test
testDir = fileparts(which('testGeneMCS'));
cd(testDir);

% define the solver packages to be used to run this test
solverPkgs = {'ibm_cplex', 'glpk', 'gurobi'};

% Load Toy Example
fileName = [CBTDIR filesep 'tutorials' filesep 'analysis' filesep 'gMCS' filesep 'gMCStoyExample.mat'];
if 0
    %TODO gMCStoyExample.mat needs to be made compatible with
%     testGeneMCS.m:
%     Error using readCbModel (line 232)
%     There were no valid models in the mat file.
%     Please load the model manually via ' load /home/rfleming/work/sbgCloud/code/fork-cobratoolbox/tutorials/analysis/gMCS/gMCStoyExample.mat' and check it with verifyModel() to validate it.
%     Try using convertOldStyleModel() before verifyModel() to bring the model structure up to date.
%     Error in testGeneMCS (line 25)
%     model = readCbModel([CBTDIR filesep 'tutorials' filesep 'analysis' filesep 'gMCS' filesep 'gMCStoyExample.mat']);
    model = readCbModel(fileName);
else
    load(fileName)
end
% expected solution
true_gmcs = cell(3,1);
true_gmcs{1} = {'g5'}';
true_gmcs{2} = {'g1' 'g4'}';
true_gmcs{3} = {'g2' 'g3' 'g4'}';

for k = 1:length(solverPkgs)
    fprintf(' -- Running testGeneMCS using the solver interface: %s ... ', solverPkgs{k});

    solverLPOK = changeCobraSolver(solverPkgs{k}, 'MILP', 0);

    if solverLPOK
        % Eliminate G-matrix if it exist
        if exist([currentDir filesep 'G_toy_example_gMCS.mat'], 'file')
            delete G_toy_example_gMCS.mat
        end

        % Check errors when missing argument
        assert(verifyCobraFunctionError('calculateGeneMCS', 'inputs', {model, 20, 5}));
        assert(verifyCobraFunctionError('calculateGeneMCS', 'inputs', {'toy_example_gMCS', [], 20,5}));
        assert(verifyCobraFunctionError('calculateGeneMCS', 'inputs', {'toy_example_gMCS', model, [],5}));
        assert(verifyCobraFunctionError('calculateGeneMCS', 'inputs', {'toy_example_gMCS', model, 20}));

        % Calculate GMCS
        [gmcs, gmcs_time] = calculateGeneMCS('toy_example_gMCS', model, 20,5);

        % Check if the solution obtained is the same as the expected
        % solution
        gmcsIsTrue = zeros(size(true_gmcs));
        for i = 1:numel(gmcs)
            for j=1:numel(true_gmcs)
                aux1 = gmcs{i};
                aux2 = true_gmcs{j};
                if isequal(aux1,aux2)
                    gmcsIsTrue(j) = gmcsIsTrue(j)+1;
                    break
                end
            end
        end
        assert(sum(~logical(gmcsIsTrue))==0);
        %Now, test with a gene_set
        options = struct();
%         options.gene_set = model.genes([1 2 4 5 6]);
        [gmcs, gmcs_time] = calculateGeneMCS('toy_example_gMCS', model, 20, 5, 'gene_set', model.genes([1 2 4 5 6]));
            % Check the gMCS
        [IsCutSet, IsMinimal, geneNotMinimal] = checkGeneMCS(model, gmcs);
        assert(all(IsMinimal));
        assert(all(IsCutSet));
        %Assert that all correct solutions are there
        assert(all(cellfun(@(x) any(cellfun(@(y) isempty(setxor(x,y)),gmcs)), {{'g5'},{'g1','g4'}})))
        %and, that there are no surplus solutions
        assert(all(cellfun(@(x) any(cellfun(@(y) isempty(setxor(x,y)),{{'g5'},{'g1','g4'}})), gmcs)))
        %Finally test this for gMCS containing a specific knockout.
        options.KO = 'g5';
        [gmcs, gmcs_time] = calculateGeneMCS('toy_example_gMCS', model, 20, 5, options);
        assert(isequal(gmcs,{{'g5'}}));
        % Check the gMCS
        [IsCutSet, IsMinimal, geneNotMinimal] = checkGeneMCS(model, gmcs);
        assert(IsMinimal);
        assert(IsCutSet);
        %assert using one worker
        options = struct();
        options.numWorkers = 1;
        assert(~verifyCobraFunctionError('calculateGeneMCS', 'inputs', {'toy_example_gMCS', model, 20, 5, options}));
    else
        warning('The test testGeneMCS cannot run using the solver interface: %s. The solver interface is not installed or not configured properly.\n', solverPkgs{k});
    end

    % Eliminate generated files
    if exist([testDir filesep 'G_toy_example_gMCS.mat'], 'file')
        delete G_toy_example_gMCS.mat
    end
    if exist([testDir filesep 'CobraMILPSolver.log'], 'file')
        delete CobraMILPSolver.log
    end
    if exist([testDir filesep 'MILPProblem.mat'], 'file')
        delete MILPProblem.mat
    end
    if exist([testDir filesep 'tmp.mat'], 'file')
        delete tmp.mat
    end

    % output a success message
    fprintf('Done.\n');
end

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function show_map(calc_type, in2, compare_window_dur)
    % ZMAP script show_map.m. Creates Dialog boxes for Z-map calculation
    % does the calculation and makes displays the map
    % stefan wiemer 11/94
    %
    % make dialog interface and call maxzlta
    %
    % turned into function by Celso G Reyes 2017

    %FIXME (maybe) changing compare_window_dur doesn't change the global version.

    % Input Rubberband
    %
    report_this_filefun();

    % compare_window_dur is shared among the parts of this function
    it=[]

    if in2 ~= 'calma'
        create_the_figure(calc_type)
    else
        do_calma(calc_type)
    end

    function create_the_figure(calc_type)
        %initial values
        ZG=ZmapGlobal.Data; % used by get_zmap_globals
        ZG.compare_window_dur_v3 = years(1.5);
        it = t0b +1;
        figure(mess);
        clf
        set(gca,'visible','off')
        set(gcf,'Units','pixel','NumberTitle','off','Name','Input Parameters');

        set(gcf,'pos',[ ZG.welcome_pos, ZG.welcome_len +[200, -50]]);

        % creates a dialog box to input some parameters
        %
        add_timecut_controls(it)

        if calc_type == "rub" || calc_type == "lta"
            add_compareduration_controls(compare_window_dur)
        end

        uicontrol('Style','Pushbutton',...
            'Position', [.60 .05 .15 .15 ],...
            'Units','normalized',...
            'Callback',@(~,~)close(),...
            'String','Cancel');

        uicontrol('Style','Pushbutton',...
            'Position',[.25 .05 .15 .15 ],...
            'Units','normalized',...
            'callback',@cb_go,...
            'String','Go');

        set(gcf,'visible','on');
        watchoff
    end

    function do_calma(calc_type)
        % check if time are with limits
        %
        if ~exist('it', 'var')
            it = t0b + (teb-t0b)/2;
        end
        if it + compare_window_dur > teb || it < t0b
            errordlg('Time out of limits')
            in2 = 'nocal';
            show_map
        end

        % initial parameter
        winlen_days = floor(compare_window_dur/ZG.bin_dur);
        ti = floor((it -t0b)/days(ZG.bin_dur));
        [len, ncu] = size(cumuall);
        len = len-2;
        var1 = zeros(1,ncu);
        var2 = zeros(1,ncu);
        mean1 = zeros(1,ncu);
        mean2 = zeros(1,ncu);
        as = zeros(1,ncu);

        % loop over all grid points for percent
        %
        %
        switch calc_type
            case 'per'
                [as, strib, stri2] = calc_percent_change();
            case 'rub'
                [as, var1, var2] = calc_rubberband(ncu, ti, winlen_days);
            case 'ast'
                [as, var1, var2] = calc_ast(ncu, ti, len);
            case 'lta'
                [as, var1, var2] = calc_lta(ncu, ti, winlen_days, len, cumuall);
            case 'maz'
                [maxlta, maxlta2, mean1, mean2] = calc_maz(ncu, cumuall, len);
            otherwise
                error('ZMAP:show_map:unknownOperation','not sure which operation to do')
        end

        %normalisation of lap1
        normlap1=nan(length(tmpgri(:,1)),1)
        normlap2=nan(length(tmpgri(:,1)),1)

        lv = logical(ll) ;  % having to do with the grid(?)
        normlap2(lv)= as(:);
        %construct a matrix for the color plot
        valueMap=reshape(normlap2,length(yvect),length(xvect));

        [n1, n2] = size(cumuall);
        s = cumuall(n1,:);
        normlap2(lv)= s(:);
        %construct a matrix for the color plot
        r = reshape(normlap2, length(yvect), length(xvect));
        ZG.tresh_km = max(r(:));
        % find max and min of data for automatic scaling
        %
        ZG.maxc = max(valueMap(:));
        ZG.maxc = fix(max(valueMap(:)))+1;
        ZG.minc = min(valueMap(:));
        ZG.minc = fix(min(valueMap(:)))-1;
        %plot imge
        %
        det = 'nop';
        old = valueMap;
        view_max(valueMap,gx,gy,stri,'nop')

    end

    function add_timecut_controls(default_value)
        uicontrol('Style','edit',...  # was inp2_field
            'Position',[.80 .775 .18 .15],...
            'Units','normalized',...
            'String',num2str(default_value),...
            'callback',@cb_timecut_year);

        text(...  #was txt2
            'Position',[0. 0.9 0 ],...
            'FontWeight','bold',...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'String','Please input time of cut in years (e.g. 86.5):');
    end

    function add_compareduration_controls(compare_window_dur)
        text(...  # was txt3
            'Position',     [0. 0.65 0 ],...
            'FontWeight',   'bold',...
            'FontSize',     ZmapGlobal.Data.fontsz.m ,...
            'String',       'Please input window length in years (e.g. 1.5):');
        uicontrol('Style',  'edit',...
            'Position',     [.80 .575 .18 .15],...
            'Units',        'normalized',...
            'String',       num2str(years(compare_window_dur)),...
            'callback',     @cb_window_duration_edit);
    end

    %% callbacks

    function cb_timecut_year(mysrc,~)
        it = str2double(mysrc.String);
    end

    function cb_window_duration_edit(mysrc,~)
        compare_window_dur = years(str2double(mysrc.String));
    end

    function cb_go(mysrc, myevt)
        watchon;
        drawnow;
        do_calma(calc_type);
    end
end

%% calculation functions

function [as, strib, stri2] = calc_percent_change(cumuall, ti, len)
    for i = ncu:-1:1
        mean1(i) = mean(cumuall(1:ti,i));
        mean2(i) = mean(cumuall(ti:len,i));
    end    %for i

    as = -((mean1-mean2)./mean1)*100;

    strib = 'Change in Percent';
    stri2 = ['ti=' num2str(ti*days(ZG.bin_dur) + t0b)  ];
end

function [as, var1, var2] = calc_rubberband(ncu, ti, winlen_days)
    for i = ncu:-1:1
        mean1(i) = mean(cumuall(1:ti,i));
        mean2(i) = mean(cumuall(ti+1:ti+winlen_days,i));
        var1(i) = cov(cumuall(1:ti,i));
        var2(i) = cov(cumuall(ti+1:ti+winlen_days,i));
    end %  for i ;
    as = (mean1 - mean2)./(sqrt(var1/ti+var2/winlen_days));
end

function [as, var1, var2] = calc_ast(ncu, ti, len)
    for i = ncu:-1:1
        mean1(i) = mean(cumuall(1:ti,i));
        var1(i) = cov(cumuall(1:ti,i));
        mean2(i) = mean(cumuall(ti+1:len,i));
        var2(i) = cov(cumuall(ti+1:len,i));
    end
    as = (mean1 - mean2)./(sqrt(var1/ti+var2/(len-ti)));
end

function [as, var1, var2] = calc_lta(ncu, ti, winlen_days, len, cumuall)
    mean1 = mean([cumuall(1:ti-1,:) ; cumuall(ti+winlen_days+1:len,:)]);
    mean2 = mean(cumuall(ti:ti+winlen_days,:));
    for i = ncu : -1 : 1
        var1(i) = cov([cumuall(1:ti-1,i) ; cumuall(ti+winlen_days+1:len,i)]);
        var2(i) = cov(cumuall(ti:ti+winlen_days,i));
    end     % for i
    as = (mean1 - mean2)./(sqrt(var1/(len-winlen_days)+var2/winlen_days));
end

function [maxlta, maxlta2, mean1, mean2] = calc_maz(ncu, cumuall, len)
            maxlta = zeros(1,ncu);
            maxlta = maxlta -5;
            mean1 = mean(cumuall(1:len,:));
            wai = waitbar(0,'Please wait...');
            set(wai,'Color',[0.8 0.8 0.8], 'NumberTitle','off', 'Name','Percent done');

            for i = ncu:-1:1
                var1(i) = cov(cumuall(1:len,i));
            end

            for ti = 1:step: len - winlen_days
                waitbar(ti/len)
                mean1 = mean(cumuall(1:len,:));
                mean2 = mean(cumuall(ti:ti+winlen_days,:));
                for i = 1:ncu
                    var1(i) = cov(cumuall(1:len,i));
                    var2(i) = cov(cumuall(ti:ti+winlen_days,i));
                end
                as = (mean1 - mean2)./(sqrt(var1/len+var2/winlen_days));
                maxlta2 = [maxlta ;  as ];
                maxlta = max(maxlta2);
            end
            %as = reshape(maxlta,length(gy),length(gx));
            close(wai)
end

author: CelsoReyes


text:
%CREATECLIQUETREE Takes in a list of factors F, Evidence and returns a
%clique tree after calling ComputeInitialPotentials at the end.
%
%   C = CREATECLIQUETREE(F) Takes a list of factors and creates a clique
%   tree . The value of the cliques should be initialized to
%   the initial potential.
%   It returns a clique tree that has the following fields:
%   - .edges: Contains indices of the nodes that have edges between them.
%   - .factorList: Contains the list of factors used to build the Clique
%   tree.
%

% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function P = CreateCliqueTree(F)

C.nodes = {};

V = unique([F(:).var]);

% Setting up the cardinality for the variables since we only get a list
% of factors.
C.card = zeros(1, length(V));
for i = 1 : length(V),

    for j = 1 : length(F)
        if (~isempty(find(F(j).var == i)))
            C.card(i) = F(j).card(find(F(j).var == i));
            break;
        end
    end
end

C.factorList = F;

% Setting up the adjaceny matrix.
edges = zeros(length(V));

for i = 1:length(F)
    for j = 1:length(F(i).var)
        for k = 1:length(F(i).var)
            edges(F(i).var(j), F(i).var(k)) = 1;
        end
    end
end

cliquesConsidered = 0;

while cliquesConsidered < length(V)

    % Using Min-Neighbors where you prefer to eliminate the variable that has
    % the smallest number of edges connected to it.
    % Everytime you enter the loop, you look at the state of the graph and
    % pick the variable to be eliminated.

    bestClique = 0;
    bestScore = inf;
    for i=1:size(edges,1)
        score = sum(edges(i,:));
        if score > 0 && score < bestScore
            bestScore = score;
            bestClique = i;
        end
    end

    cliquesConsidered = cliquesConsidered + 1;
    [F, C, edges] = EliminateVar(F, C, edges, bestClique);

end

% Pruning the tree.
C = PruneTree(C);

% Assume that C now has correct cardinality, variables, nodes and edges.
% Here we make the function call to assign factors to cliques and compute the
% initial potentials for clusters.

P = ComputeInitialPotentials(C);

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function handle = vectorModify(handle, values)

% VECTORMODIFY Helper code for visualisation of vectorial data.

% MLTOOLS

set(handle, 'YData', values);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [clpot, seppot] = init_pot(engine, clqs, pots, pot_type, onodes, ndx)
% INIT_POT Initialise potentials with evidence (jtree_inf)
% function [clpot, seppot] = init_pot(engine, clqs, pots, pot_type, onodes)

cliques = engine.cliques;
bnet = bnet_from_engine(engine);
% Set the clique potentials to all 1s
C = length(cliques);
clpot = cell(1,C);
for i=1:C
  clpot{i} = mk_initial_pot(pot_type, cliques{i}, bnet.node_sizes(:), bnet.cnodes(:), onodes);
end

% Multiply on specified potentials
for i=1:length(clqs)
  c = clqs(i);
  clpot{c} = multiply_by_pot(clpot{c}, pots{i});
end

seppot = cell(C,C); % implicitely initialized to 1

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function plotData(x, y)
%PLOTDATA Plots the data points x and y into a new figure
%   PLOTDATA(x,y) plots the data points and gives the figure axes labels of
%   population and profit.

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the training data into a figure using the
%               "figure" and "plot" commands. Set the axes labels using
%               the "xlabel" and "ylabel" commands. Assume the
%               population and revenue data have been passed in
%               as the x and y arguments of this function.
%
% Hint: You can use the 'rx' option with plot to have the markers
%       appear as red crosses. Furthermore, you can make the
%       markers larger by using plot(..., 'rx', 'MarkerSize', 10);

figure; % open a new figure window

plot(x,y,'rx','MarkerSize',10);%rx means color is red and show symbol is x,and x's size is 10
ylabel=('Profit in $10,000s');
xlabel=('Population of City in 10,000s');

% ============================================================

end

repo: MachineLearningEx


text:
function bvec_sub_test ( )

%*****************************************************************************80
%
%% BVEC_SUB_TEST tests BVEC_SUB;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2015
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  seed = 123456789;
  test_num = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BVEC_SUB_TEST\n' );
  fprintf ( 1, '  BVEC_SUB subtracts binary vectors representing integers;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        I        J        I - J   BVEC_SUB\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    [ i, seed ] = i4_uniform_ab ( -100, 100, seed );
    [ j, seed ] = i4_uniform_ab ( -100, 100, seed );

    k = i - j;

    bvec1 = i4_to_bvec ( i, n );
    bvec2 = i4_to_bvec ( j, n );
    bvec4 = bvec_sub ( n, bvec1, bvec2 );
    l = bvec_to_i4 ( n, bvec4 );

    fprintf ( 1, '  %8d  %8d  %8d  %8d\n', i, j, k, l );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function A = import_data(fname)
%IMPORT_DATA Import tensor-related data to a file.
%
%   A = IMPORT_DATA(FNAME) imports an object A from the file named FNAME.
%   The supported data types and formatting of the file are explained in
%   EXPORT_DATA.
%
%   See also TENSOR, EXPORT_DATA
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Open file
fid = fopen(fname,'r');
if (fid == -1)
    error('Cannot open file %s',fname);
end

%% Get the type of object
line = fgets(fid);
type = sscanf(line, '%s');

%% Import the object

if strcmpi(type,'tensor')

    sz = import_size(fid);
    data = import_array(fid, prod(sz));
    A = tensor(data, sz);

elseif strcmpi(type,'matrix')

    sz = import_size(fid);
    data = import_array(fid, prod(sz));
    A = reshape(data, sz);

else

    error('Invalid data type for export');

end

%% Close file
fclose(fid);

function sz = import_size(fid)
% Import the size of something from a file
line = fgets(fid);
n = sscanf(line, '%d');
line = fgets(fid);
sz = sscanf(line, '%d');
sz = sz';
if (size(sz,2) ~= n)
    error('Imported dimensions are not of expected size');
end

function data = import_array(fid, n)
% Export dense data that supports numel and linear indexing
data = fscanf(fid, '%e', n);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
classdef DehomogenizingSingularitiesTest < handle

    properties (Access = private)
        mesh
        orientation
        levelSet
    end

    properties (Access = private)
        testName
        meshSize
        singularitiesData
        xmin
        xmax
        ymin
        ymax
        widthH
        widthW
        nCells
    end

    methods (Access = public)

        function obj = DehomogenizingSingularitiesTest(cParams)
            obj.init(cParams);
            obj.createMesh();
            obj.createOrientation();
            obj.dehomogenize();
        end

        function passed = hasPassed(obj)
            d = load(obj.testName);
            ls = obj.levelSet;
            errI = zeros(numel(ls),1);
            for iCell = 1:numel(ls)
                lSI = d.levelSet{iCell};
                errI(iCell) = norm(ls{iCell} - lSI)/norm(lSI);
            end
            itIs = sum(errI) < 1e-12;
            passed = itIs;
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.testName = cParams.testName;
            %obj.meshSize = 0.00521;
            obj.meshSize = 0.09;%0.0221;%0.09;%0.0221;%0.0521 %0.0221;0.0921
            obj.nCells   = [60 62];%linspace(60,62,40);%45;   %45
            obj.xmin = 0.5;
            obj.xmax = 2.0;
            obj.ymin = 0.25;
            obj.ymax = 1.75;
            obj.singularitiesData = [0.32,-0.8];
            obj.widthH = 0.87;
            obj.widthW = 0.87;
        end

        function createMesh(obj)
            h = obj.meshSize;
            xv = obj.xmin:h:obj.xmax;
            yv = obj.ymin:h:obj.ymax;
            [X,Y] = meshgrid(xv,yv);
            s.coord(:,1) = X(:);
            s.coord(:,2) = Y(:);
            s.connec = delaunay(s.coord);
            m = Mesh(s);
            obj.mesh = m;
        end

        function createOrientation(obj)
            m = obj.mesh;
            s1 = obj.singularitiesData(:,1);
            s2 = obj.singularitiesData(:,2);
            x1 = m.coord(:,1);
            x2 = m.coord(:,2);
            v(:,1) = cos(pi*(x1 + s1*x2));
            v(:,2) = cos(pi*(x2 + s2*x1));
            beta = atan2(v(:,2),v(:,1));
            alpha = beta/2;
            obj.orientation(:,1) = cos(alpha);
            obj.orientation(:,2) = sin(alpha);
        end

        function plotOrientation(obj)
            figure()
            x = obj.mesh.coord(:,1);
            y = obj.mesh.coord(:,2);
            t  = obj.orientation;
            ct = cos(t(:,1));
            st = sin(t(:,1));
            quiver(x,y,ct,st)
        end

        function s = createLevelSetCellParams(obj)
            I        = ones(size(obj.mesh.coord,1),1);
            s.type   = 'rectangleInclusion';
            s.widthH = obj.widthH*I;
            s.widthV = obj.widthW*I;
            s.ndim   = 2;
        end

        function dehomogenize(obj)
            s.nCells             = obj.nCells;
            s.cellLevelSetParams = obj.createLevelSetCellParams();
            s.mesh               = obj.mesh;
            s.theta              = atan2(obj.orientation(:,2),obj.orientation(:,1));
            d = Dehomogenizer(s);
            ls = d.compute();
            d.plot();
            obj.levelSet = ls;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function S = MSAFinit(w0,mu,N,L,H,F)

% MSAFinit      Initialize Parameter Structure for the MSAF Algorithm
%                 Psuedo-QMF CMFB is Used by Default, H and F are Used Otherwise
%
% Arguments:
% w0            Coefficients of FIR filter at start
% mu            Step size
% N             Number of subbands
% H             Analysis filter bank (optional), each column represents a filter
% F             Synthesis filter bank (optional), each column represents a filter

if nargin > 4
    if (size(H,2)~=N)|(size(F,2)~=N)
        error('Columns of H (%d) or F (%d) not match with N = %d',size(H,2),size(F,2),N);
    end
else

% Defualt filter bank: Pseudo-QMF CMFB

    [hopt,passedge] = opt_filter(L-1,N); % Generate a prototype lowpass filter
    [H,F] = make_bank(hopt,N);           % Generate filter banks using cosine modulation
    H = sqrt(N)*H';                      % Analysis section
    F = sqrt(N)*F';                      % Synthesis section
end

% Assign structure fields

S.coeffs        = w0(:);                 % Convert to column vector of length M
S.step          = mu;
S.analysis      = H;
S.synthesis     = F;
S.iter          = 0;                     % Iteration count
S.alpha         = ones(1,N)*1e-4;        % Small positive constant
S.AdaptStart    = L + length(w0);        % Running effect of analysis and adaptive filter,
                                         %   minimum L + M

author: CharlesThaCat

repo: acoustic-interference-cancellation

path: github-repos/MATLAB/CharlesThaCat-acoustic-interference-cancellation/acoustic-interference-cancellation-edb394499ea6f9c96445a3e9613bd64a854c289e/Subband processing/Book/After reading Chapter 6 -- Multiband-structured subband adaptive filters/MSAFinit.m


text:
function h = polarScatter(az, decl, point_size, color, flag, plot_bg)
% Equivalent of polar scatter but with support for older Matlab versions
% SYNTAX:
%    polarScatter(az, decl, point_size, color, <flag>)
%
% INPUT
%   az      azimuth      [rad]
%   decl    declination  [rad]
%   size    size of the scattered point
%   color   data field for scatter
%   flag    at the moment can only be 'filled'
%
% OUTPUT
%   h       handle to the scattered points

%--------------------------------------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Mirko Reguzzoni
%  Contributors:     Giulio Tagliaferro, Andrea Gatti ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%    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/>.
%----------------------------------------------------------------------------------------------

%%% INTERNAL PARAMETER
    scale = 1;
    %%%

    decl_n = decl/(pi/2)*scale;
    x = sin(az) .* decl_n;
    y = cos(az) .* decl_n;
    if nargin > 4 && strcmp(flag,'filled')
        h = scatter(x,y,point_size,color,'filled');
    else
        h = scatter(x,y,point_size,color);
    end
    is_hold = ishold();
    if nargin < 6
        plot_bg = ~is_hold;
    end
    if plot_bg
        hold on
        %plot parallel
        az_l = [0:pi/200:2*pi];
        d_step = 15/180*pi;
        decl_s = ([0:d_step:pi/2]/(pi/2))*scale;
        for d = decl_s
            x = cos(az_l).*d;
            y = sin(az_l).*d;
            plot(x,y,'color',[0.6 0.6 0.6]);
            text(cos(80/180*pi)*d,sin(80/180*pi)*d,sprintf('%d',round(d*90)),'HorizontalAlignment','center', 'FontWeight', 'bold', 'FontSize', 13);
        end
        %plot meridian
        az_step = 30/180 *pi;
        az_s = [0:az_step:2*pi];
        decl_l = ([0 1])*scale;
        for a = az_s
            x = cos(a).*decl_l;
            y = sin(a).*decl_l;
            plot(x,y,'color',[0.6 0.6 0.6]);
            if abs(a-2*pi) > 0.0001
                text(cos(a)*1.1,sin(a)*1.1,sprintf('%d', mod(round((2*pi - a + pi/2) / pi * 180), 360)), 'HorizontalAlignment','center', 'FontWeight', 'bold', 'FontSize', 13);
            end
        end
        axis equal
        % xlim([-2 2])
        % ylim([-2 2])
        axis off
        set(gcf,'color','w');
        if ~is_hold
            hold off
        end
        xlim([-1.15 1.15]); ylim([-1.15 1.15]);
    end
end


text:
function  [gx] = g_Kalman4bandits (x, P, ~, in)
% // VBA toolbox //////////////////////////////////////////////////////////
%
% [gx] = g_Kalman4bandits (x, P, u, in)
% softmax decision rule for N-armed bandit task as in Daw et al. 2006
%
% /////////////////////////////////////////////////////////////////////////

% Get parameter values
% -------------------------------------------------------------------------
% inverse temperature
beta = exp (P(1)); % exp: [-Inf,Inf] -> [0 Inf]
% bonus to exploration
phi = P(2);

% Extract expectation and variance of each bandit
% -------------------------------------------------------------------------
nBandits = in.nBandits;
mu = x(1 : nBandits);
sigma2 = x(nBandits + (1 : nBandits)) .^ 2;

% apply softmax
% -------------------------------------------------------------------------
gx = VBA_softmax(beta * (mu + phi * sigma2));

author: MBB-team

repo: VBA-toolbox


text:
classdef S2FunHarmonicSym < S2FunHarmonic
% a class represeneting a symmetrised function on the sphere

properties
  s = []; % symmetry
end

properties (Dependent = true)
  CS
  SS
end

methods
  function sF = S2FunHarmonicSym(fhat, s)
    if nargin == 0, return; end
    sF.fhat = fhat;
    sF.s = s;
  end

  function CS = get.CS(sF)
    CS = sF.s;
  end

  function SS = get.SS(sF)
    SS = sF.s;
  end

end

methods (Static = true)
  sF = quadrature(f, varargin);
end

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [elem,bdFlag] = sortelem3(elem,bdFlag)
%% SORTELEM3 sort elem in ascend ordering
%
% [elem,bdFlag] = sortelem3(elem,bdFlag) sorts the elem such that
% elem(t,1)< elem(t,2)< elem(t,3)<elem(t,4). A simple sort(elem,2) cannot
% switch bdFlag.
%
% See also  sortelem
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Step 1: elem(:,4) is the largest one
[tempvar,idx] = max(elem,[],2);  %#ok<*ASGLU>
elem(idx==1,1:4) = elem(idx==1,[2 4 3 1]);
elem(idx==2,1:4) = elem(idx==2,[3 4 1 2]);
elem(idx==3,1:4) = elem(idx==3,[4 2 1 3]);
if exist('bdFlag','var')
    bdFlag(idx==1,1:4) = bdFlag(idx==1,[2 4 3 1]);
    bdFlag(idx==2,1:4) = bdFlag(idx==2,[3 4 1 2]);
    bdFlag(idx==3,1:4) = bdFlag(idx==3,[4 2 1 3]);
end

%% Step 2: elem(:,1) is the smallest one
[tempvar,idx] = min(elem(:,1:3),[],2);
% elem(idx==1,1:3) = elem(idx==1,[1 2 3]);
elem(idx==2,1:3) = elem(idx==2,[2 3 1]);
elem(idx==3,1:3) = elem(idx==3,[3 1 2]);
if exist('bdFlag','var')
    bdFlag(idx==2,1:3) = bdFlag(idx==2,[2 3 1]);
    bdFlag(idx==3,1:3) = bdFlag(idx==3,[3 1 2]);
end

%% Step 3: sort elem(:,2)<elem(:3)
idx = (elem(:,3) < elem(:,2));
elem(idx,[2 3]) = elem(idx,[3 2]);
if exist('bdFlag','var')
    bdFlag(idx,[2 3]) = bdFlag(idx,[3 2]);
end

%% Output
if ~exist('bdFlag','var')
    bdFlag = [];
end

repo: ifem


text:
% the following program calculates the
% Colebrool-White friction factors at various
% roughess heights for Reynolds numbers
% between 5000 and 500000
% ---------------------------------------------------------------
%  The MATLAB function was created by Tibor Balint, December 1998
% TBoreal Research Corporation, Toronto, Ont. Canada
% (tibor@netcom.ca) and also, University of Warwick, UK
% ---------------------------------------------------------------

clear
Re=linspace(5000,500000,200); % set the Reynolds numbers
Dh=0.008;  % set the hydraulic diameter in meters

rough=0.00000075;  % set the 1st EQUIVALENT ROUGHNESS HEIGHT in m

for i=1:200
   miller(i)=0.25/(log10(rough/(3.7*0.009)+5.74/Re(i)^0.9))^2;
   fcw(i)=ffcw(Re(i), Dh, rough);
end

rough=0.000001;  % EQUIVALENT ROUGHNESS HEIGHT
for i=1:200
%   miller1(i)=0.25/(log10(rough/(3.7*0.009)+5.74/Re(i)^0.9))^2;
   fcw1(i)=ffcw(Re(i), Dh, rough);
end

rough=0.00000125;  % EQUIVALENT ROUGHNESS HEIGHT
for i=1:200
%   miller2(i)=0.25/(log10(rough/(3.7*0.009)+5.74/Re(i)^0.9))^2;
   fcw2(i)=ffcw(Re(i), Dh, rough);
end

rough=0.0000015;  % EQUIVALENT ROUGHNESS HEIGHT
for i=1:200
%   miller3(i)=0.25/(log10(rough/(3.7*0.009)+5.74/Re(i)^0.9))^2;
   fcw3(i)=ffcw(Re(i), Dh, rough);
end

rough=0.00000175;  % EQUIVALENT ROUGHNESS HEIGHT
for i=1:200
%   miller4(i)=0.25/(log10(rough/(3.7*0.009)+5.74/Re(i)^0.9))^2;
   fcw4(i)=ffcw(Re(i), Dh, rough);
end

%plot the figures
figure(1)
orient landscape;
subplot(2,1,1);
plot(Re,fcw,'b',Re,miller,'r--')
legend('f_{cw}','f_{miller}');
grid on;
axis tight;
xlabel('Reynolds number');
ylabel('Friction factor');
title ('Hydraulic Diameter (D_h=8 mm), Eq. Roughness Height (\epsilon = 0.75 \mum)');

subplot(2,1,2);
plot(Re,fcw,'b',Re,fcw1,'r',Re,fcw2,'g',Re,fcw3,'m',Re,fcw4,'k')
legend('\epsilon=0.75 \mum','\epsilon=1.00 \mum','\epsilon=1.25 \mum','\epsilon=1.50 \mum','\epsilon=1.75 \mum');
grid on;
axis tight;
xlabel('Reynolds number');
ylabel('Friction factor');
title ('Colebrook-White Friction Factor (f_{cw})');

author: Sable

repo: mcbench-benchmarks


text:
function par_est
% parameter estimation with derivatives              Holzbecher September 2005
% for exponential fit of exponent lambda for c(0)=c0
global tfit cfit c0

% specify fitting data
tfit =  [0.25 1 2 4 8];
cfit = [0.7716 0.5791 0.4002 0.1860 0.1019];
c0 = .8;

lambda = fzero(@myfun,0.11);
normc = norm(cfit - c0*exp(-lambda*tfit));
display (['Best fit for lambda = ' num2str(lambda)]);
display (['Norm of residuals =' num2str(normc)]);
tmax = tfit(size(tfit,2));
t = [0:0.01*tmax:tmax];
figure; plot (tfit,cfit,'or',t,c0*exp(-lambda*t),'-');
legend ('given','modelled');
text(0.5*tmax,c0*0.7,['\lambda: ' num2str(lambda)]);
text(0.5*tmax,c0*0.6,['norm of residuals: ' num2str(normc)]);
xlabel ('time'); ylabel ('concentration')

function f = myfun(lambda);
global tfit cfit c0

t = tfit;
c = c0*exp(-lambda*t);    % solve linear decay equation for c with c(0)=c0
clambda = -c.*t;          % equation for dc/dlambda
f = (cfit-c)*clambda';    % specify function f to vanish


text:
function success = verifyCobraFunctionError(functionCall, varargin)
% Tests whether the provided call throws an error and if the error message
% fits to the message provided
%
% USAGE:
%
%    success = verifyCobraFunctionError(functionCall, message)
%
% INPUTS:
%    functionCall:       The name of a function.
%
% OPTIONAL INPUTS:
%
%    varargin:           Parameter Value pairs for specific options.
%                         * outputArgCount -     The number of output arguments requested, as this can influence the functionality (Default: 0).
%                         * testMessage -        The message that should be thrown in this error. If no message is supplied, any thrown error will be accepted
%                         * inputs -             If the function requires Any inputs, provide them in a cell array.
%
% OUTPUTS:
%    success:            Whether an error was thrown, and the type and
%                        message match (if provided)
% EXAMPLES:
%    %Test whether pow throws an error when using chars
%    verifyCobraFunctionError('pow','inputs',{'abc','dce'});
%    %Test, whether initCobraToolbox throws an error
%    verifyCobraFunctionError('initCobraToolbox');
%    %Test whether targetFunctions error message is 'blubb'
%    verifyCobraFunctionError('targetFunctions','message','blubb');
%    %Test whether fluxvariability throws an error when 5 outputs are defined
%    verifyCobraFunctionError('fluxVariability','inputs',{model,0},'outputArgCount',10);

testMessage = false;
parser = inputParser();
parser.addRequired('functionCall',@ischar)
parser.addParamValue('inputs',{},@iscell)
parser.addParamValue('testMessage','',@ischar)
parser.addParamValue('outputArgCount',0,@(x) isnumeric(x) && mod(x,1) == 0);

parser.parse(functionCall,varargin{:});

outputArgcount = parser.Results.outputArgCount;
testMessage = ~any(ismember('testMessage',parser.UsingDefaults));
message = parser.Results.testMessage;
inputs = parser.Results.inputs;

functionCall = str2func(functionCall);
success = true;
try
    if outputArgcount > 0
        outArgs = cell(outputArgcount,1);
        [outArgs{:}] = functionCall(inputs{:});
    else
        functionCall(inputs{:});
    end
catch ME
    if testMessage
        success = strcmp(ME.message,message);
    end
    %Now, we checked the message if necessary, and we are obvously in the
    %catch block, so an error was thrown, i.e. this is successful.
    return
end
%If we reach this point, the function did not throw an error. i.e. the
%verification is unsuccessful.
success = false;


text:
function [im] = blend(im1, im2, c1, c2)

[height, width, ~] = size(im1);

alpha1 = zeros(size(im1));
alpha2 = zeros(size(im2));

im1 = im2double(im1);
im2 = im2double(im2);

im = (alpha2 .* im2 + alpha1 .* im1) ./ (alpah1 + alpha2);

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
function plot_catalog(mCatalog, vMain,bFaults)

if bFaults
    load ~/data/faults/mCA-faults.mat
end
figure
plot(mCatalog(:,1),mCatalog(:,2),'k.');
hold on;plot(mCatalog(~vMain,1),mCatalog(~vMain,2),'r.');
xlim([-118.2 -115]);
ylim([32 36]);
set(gca,'FontSize',16);
if bFaults
    hold on; plot(mFaults(:,1),mFaults(:,2),'k-','LineWidth',1);
end

figure
plot(mCatalog(:,3),cumsum(ones(size(mCatalog,1),1)),'r',...
    'LineWidth',2);
hold on;plot(mCatalog(:,3),cumsum(vMain),'k',...
    'LineWidth',2);
set(gca,'FontSize',16);


text:
% DEMGPDISIMMEF2 Run experiments on Mef2 data. The raw data is pre-processed by the PUMA package.

% SHEFFIELDML

%/~
% path(path, '/local/Matlab/underDevelopment/GPSIM016');
%~/

clear; close all;
expNo = 2;
type = 'Dros';
warning('off');
saveFigures = 0;

tf = 'mef2';

%/~
if exist('./data/mef2Data.mat') == 2
%~/
load('./data/mef2Data.mat');
%/~
else

  drosLoadMef2Data;
  targetsFull = drosFindTargets(drosmef2chip);
  % targets = ([targetsFull(1:6); targetsFull(8:10)])';
  % selection = [8 34 19 26 2];
  selection = [25 34 19 37 21 40];
  targets = targetsFull(selection)';

  % targets = drosTargets(tf);
  tflabel = drosTF.labels(strcmp(tf, drosTF.names));

  genes = [tflabel, targets];

  [y, yvar, gene, times, scale, rawExp, rawVar] = gpdisimGetDrosData(drosexp, ...
                                                    genes);

  save('./data/mef2Data.mat', 'y', 'yvar', 'gene', 'times', 'scale', 'rawVar', ...
       'rawExp', 'genes', 'targets');

end
%~/
genenames = genes;

% Get the default options structure.
options = gpsimOptions;
options.includeNoise = 1;
% Fix one decay (from the fourth gene --- p21) to 0.8 hr^-1, and
% the corresponding sensitivity (see just after eqn 2 in the
% mathematical methods of Barenco et al.)
options.fix(1).index = 2;
options.fix(1).value = expTransform(1, 'xtoa');
options.fix(2).index = 6;
options.fix(2).value = expTransform(1, 'xtoa');
%options.fix(2).index = 1;
%options.fix(2).value = expTransform(.3, 'xtoa');;

% initialise the model.
model.type = 'cgpdisim';  % This new model type is a hack to run
                        % the model in a hierarchical manner.
                        % need to do this more elegantly later.
for i =1:3          %% 3 original
  model.comp{i} = gpdisimCreate(length(targets), 1, times, y{i}, yvar{i}, options);
end

% Learn the model.
model = modelOptimise(model, [], [], 1, 3000);

tfName = tf;
tfName(1) = upper(tfName(1));
fileName = ['dem' tfName type num2str(expNo)];
save(fileName);

% Plot
% figure(1); clf; plot_tf_exps(exp_struct, expse_struct);
% figure(2); clf; plot_expros(genes, exp_struct, expse_struct, genenames);
% figure(3); clf;

numGenes = model.comp{1}.numGenes;
genenames{2} = 'Rya-r44F';
genenames{4} = 'ttk';

for j = 1:length(model.comp)

  % Generate predictions of the functions.
  % to do this we need to compute the K_xf portions of the kernel
  % (simXrbfKernCompute does this for us).
  predt = [1:0.1:12 model.comp{j}.t']';
  proteinKern = kernCreate(model.comp{j}.t, 'sim');
  proteinKern.inverseWidth = model.comp{j}.kern.comp{1}.inverseWidth;
  proteinKern.decay = model.comp{j}.delta;
  proteinKern.variance = model.comp{j}.kern.comp{2}.di_variance;
  K = simXrbfKernCompute(proteinKern, model.comp{j}.kern.comp{1}, ...
                               predt, model.comp{j}.t);
  for i=2:model.comp{j}.kern.numBlocks
    blockK =  disimXsimKernCompute(model.comp{j}.kern.comp{i}, proteinKern, ...
                               model.comp{j}.t, predt);
    K = [K blockK'];
  end

  ymean = reshape(ones(length(times),1)*y{j}(1,:), length(model.comp{j}.y), ...
                                                     1);
%  ymean = mean(model.comp{j}.y);

  predF = K*model.comp{j}.invK*(model.comp{j}.y-ymean);
  varF = kernDiagCompute(proteinKern, predt) - sum(K'.* ...
                                                   (model.comp{j}.invK* ...
                                                    K'))';
%  varF = proteinKern.variance - sum(K'.*(model.comp{j}.invK*K'))';

%   % predict gene via the model
%   predF(end-length(model.comp{j}.t)+1:end) = [];
%   varF(end-length(model.comp{j}.t)+1:end) = [];
%   predt(end-length(model.comp{j}.t)+1:end) = [];
%   model.comp{j}.mapt = predt;
%   model.comp{j}.g = predF;
%   model.comp{j}.numMapPts = length(model.comp{j}.mapt);
%   model.comp{j}.step = 0.1;
%
%   model.comp{j}.times_index = [];
%   for i = 1:length(times)
%     model.comp{j}.times_index(i) = find((times(i) - model.comp{j}.mapt)==0);
%   end
%
%   model.comp{j} = gpsimMapUpdateYpred(model.comp{j});
%
%   predExprs = zeros(length(predt), model.comp{j}.numGenes+1);
%   varExprs = zeros(length(predt), model.comp{j}.numGenes+1);
%
%   predExprs(:,2:end) = model.comp{j}.ypred;
%   predFdelay = zeros(size(predF));
%   predFdelay(1) = predF(1);
%   predFdelay(2:end) = predF(1:end-1);
%   df = (predF - predFdelay)/0.1;
%
%   predDi = (df + model.comp{j}.delta*predF)/ ...
%       model.comp{j}.sigma;
%   scaleExprs = sqrt(var(predExprs))./sqrt(var(y{j}));
%   meanExprs = mean(predExprs);
%
%   model.comp{j}.sigma = model.comp{j}.sigma*scaleExprs(1);
%   model.comp{j}.B = (model.comp{j}.B-model.comp{j}.D.*meanExprs(2:end))./ ...
%       scaleExprs(2:end) + model.comp{j}.D.*mean(y{j}(2:end));
%   model.comp{j}.B = model.comp{j}.B./scaleExprs(2:end);
%   model.comp{j}.S = model.comp{j}.S./scaleExprs(2:end);
%
%   model.comp{j}.kern.comp{2}.di_variance = model.comp{j}.sigma^2;
%   for i = 2:model.comp{j}.kern.numBlocks
%     model.comp{j}.kern.comp{i}.variance = model.comp{j}.S(i-1)^2;
%   end
%
%   model.comp{j} = gpsimMapUpdateYpred(model.comp{j});
%   predExprs(:,2:end) = model.comp{j}.ypred;
%   predExprs(:,1) = (df + model.comp{j}.delta*model.comp{j}.g)/ ...
%       model.comp{j}.sigma;

  % Predicted Gene Expressions
  Kxx = multiKernCompute(model.comp{j}.kern, predt, model.comp{j}.t);
  meanPredX = reshape(ones(length(predt),1)*([0 model.comp{j}.B./ ...
                      model.comp{j}.D]), length(predt)*(numGenes+1), 1);
  predX = meanPredX + real(Kxx*model.comp{j}.invK*(model.comp{j}.y-ymean));
  varX = real(kernDiagCompute(model.comp{j}.kern, predt) - sum(Kxx'.* ...
                           (model.comp{j}.invK*Kxx'), 1)');

  % Take out predictions at data points.
  % Use them to get the scale for the other data.
  numData = length(predX)/(numGenes+1);
  predExprs = reshape(predX, numData, numGenes+1);
  meanExprs = ones(numData, 1)*mean(predExprs);
  scaleExprs = ones(numData, 1)*(sqrt(var(predExprs))./sqrt(var(y{j})));
  % Driving input can only be adjusted by the scale constant.
%  predExprs(:,1) = predExprs(:,1)./scaleExprs(:,1);
  % predictions of other genes can be generally normalised.
%   predExprs(:,2:end) =  ones(numData, 1)*mean(y{j}(:,2:end)) + ...
%       (predExprs(:,2:end) - meanExprs(:,2:end))./scaleExprs(:,2:end);
  predExprs(end-length(model.comp{j}.t)+1:end,:) = [];
  varExprs = reshape(varX, numData, numGenes+1);
  varExprs = varExprs./scaleExprs./scaleExprs;
  varExprs(end-length(model.comp{j}.t)+1:end,:) = [];
  predF(end-length(model.comp{j}.t)+1:end) = [];
  varF(end-length(model.comp{j}.t)+1:end) = [];
  predt(end-length(model.comp{j}.t)+1:end) = [];

  yscale = ones(length(times), 1)*sqrt(var(y{j}));
  ymean = ones(length(times),1)*mean(y{j});
  ynormal = ymean + (y{j}-ymean)./yscale;
  yvarNormal = yvar{j}./yscale./yscale;

  scalePred = sqrt(var(predExprs));

  figure;
%  subplot(length(model.comp), 1, j);
  lin = plot(predt, predF, '-');
  hold on,
  bh = plot(predt, predF + 2*sqrt(varF), '--');
  bh = [bh plot(predt, predF - 2*sqrt(varF), '--')];
  hold off
%  ylabel(['Replica' num2str(j)]);
  title('Inferred Mef2 Protein', 'fontsize', 20);
  set(bh, 'lineWidth', 2);
  set(lin, 'lineWidth', 3);
  %set(lin, 'markersize', 20);
  set(gca, 'fontname', 'arial', 'fontsize', 24, 'xlim', [min(model.comp{j}.t) ...
                                                    max(model.comp{j} ...
                                                    .t)]);
  set(gca, 'ylim', [-0.1 0.4]);
  if saveFigures==1
    saveFileName = [fileName 'TF_profile_Rep' num2str(j)];
    print('-dpng', ['./results/' saveFileName]);
    pos = get(gcf, 'paperposition');
    origpos = pos;
    pos(3) = pos(3);
    pos(4) = pos(4);
    set(gcf, 'paperposition', pos);
    lineWidth = get(gca, 'lineWidth');
    set(gca, 'lineWidth', lineWidth*2);
    print('-deps', ['./results/' saveFileName]);
    set(gca, 'lineWidth', lineWidth);
    set(gcf, 'paperposition', origpos);
  end

  for index = 1:numGenes+1
    figure;
    lin = plot(predt, predExprs(:,index), '-');
    hold on,
    bh = plot(predt, predExprs(:,index)+2*real(sqrt(varExprs(:,index))), '--');
    bh = [bh plot(predt, predExprs(:,index)-2*real(sqrt(varExprs(:,index))), '--')];
    lin = [lin plot(times, y{j}(:,index), 'rx')];
    lin1 = errorbar(times, y{j}(:,index), 2*sqrt(yvar{j}(:,index)), ...
                    'rx');
    if index == 1
      titleText = ['Driving Input mRNA'];
%      lin = [lin plot(predt, predDi, 'm-')];
    else
      titleText = ['Gene ' genenames{index} ' mRNA'];
    end
    title(titleText, 'fontsize', 20);
    set(bh, 'lineWidth', 3);
    set(lin, 'lineWidth', 4);
    set(lin, 'markersize', 20);
    set(lin1, 'lineWidth', 2);
    set(gca, 'fontname', 'arial', 'fontsize', 24, 'xlim', [min(predt) ...
                        max(predt)]);

    if saveFigures
      saveFileName = [fileName '_ExprsProfile_Rep' num2str(j) '_Gene' num2str(index)];
      print('-deps', ['./results/' saveFileName]);
      pos = get(gcf, 'paperposition');
      origpos = pos;
      pos(3) = pos(3);
      pos(4) = pos(4);
      set(gcf, 'paperposition', pos);
      lineWidth = get(gca, 'lineWidth');
      set(gca, 'lineWidth', lineWidth*2);
      print('-dpng', ['./results/' saveFileName]);
      set(gca, 'lineWidth', lineWidth);
      set(gcf, 'paperposition', origpos);
    end
  end

end


text:
dirs = {'potentials/Tables', ...
	'CPDs/@discrete_CPD', ...
      'inference/static/@jtree_sparse_inf_engine', ...
      'inference/static/@gibbs_sampling_inf_engine/private'};

BNT_HOME = '/home/ai2/murphyk/matlab/FullBNT'; % edit this
%global BNT_HOME

for d=1:length(dirs)
  f = fullfile(BNT_HOME, 'BNT', dirs{d});
  fprintf('removing Cmex files from %s\n', f);
  cd(f)
  delete *.mex*
  delete *.dll
  delete *.obj
  delete *.o
end


text:
function options = drawCbMap(map,options,varargin)
%Draws a map with the specified optional parameters
%
% drawCbMap(map,options)
%
%INPUTS
% map                       COBRA map structure
%
%OPTIONAL INPUTS
% options                  Structure containing optional parameters
%   nodeWeight              Size of primary metabolite nodes
%   nodeWeightSecondary     Size of secondary metabolite nodes
%   nodeColor               Color of metabolite nodes
%   edgeColor               Color of reaction arrows2
%   edgeArrowColor          Color of reaction arrowheads
%   edgeWeight              Width of reaction arrows
%   textSize                Font size of metabolite text
%   textColor               Text color for metaboltes
%   rxnTextSize             Font size of reaction text
%   rxnTextColor            Text color for reactions
%   otherTextColor          Color of other text
%   fileName                Name of output file
%
% varargin                  optional parameter name / parameter value pairs
%
%OUTPUT
% Displays map in a matlab figure ('matlab') or target.svg ('svg') file
% depending on value of CB_MAP_OUTPUT.
% options                  Structure containing optional parameters
%
%
% Turner Conrad     6/12/12     Added rxnTextSize default (8) to fix error
%                               in writing reaction texts to .svg

global mapHandle
global CB_MAP_OUTPUT

%check for render output
if ~exist('CB_MAP_OUTPUT', 'var') || isempty(CB_MAP_OUTPUT)
    error('No render target specified.  Call changeCbMapOutput(outputFormat)');
end

if (nargin < 2)
    options = [];
end

if mod(length(varargin),2)==0
    for i=1:2:length(varargin)-1
        switch lower(varargin{i})
            case 'nodeweight', options.nodeWeight = cell2mat(varargin(i+1));
            case 'nodecolor', options.nodeColor = cell2mat(varargin(i+1));
            case 'edgeweight', options.edgeWeight = cell2mat(varargin(i+1));
            case 'edgecolor', options.edgeColor = cell2mat(varargin(i+1));
            case 'edgearrowcolor', options.edgeArrowColor = cell2mat(varargin(i+1));
            case 'textsize', options.textSize = cell2mat(varargin(i+1));
            case 'textcolor', options.textColor = cell2mat(varargin(i+1));
            case 'othertextsize', options.otherTextSize = cell2mat(varargin(i+1));
            case 'othertextcolor', options.otherTextColor = cell2mat(varargin(i+1));
            case 'filename', options.fileName = varargin{i+1};
            otherwise, options.(varargin{i}) = varargin{i+1};
        end
    end
else
    error('Invalid number of parameters/values');
end

%%%% Compelete the missing parts of the option
nNodes = size(map.molName,1);
nEdges = size(map.connection,1);
%Node size
if ~isfield(options,'nodeWeight')
    options.nodeWeight = ones(nNodes,1)*15;
    if strcmp(CB_MAP_OUTPUT,'svg')
        options.nodeWeight = ones(nNodes,1)*25;
    end
end

if ~isfield(options,'nodeWeightSecondary')
    options.nodeWeightSecondary = ones(nNodes,1)*10;
    if strcmp(CB_MAP_OUTPUT,'svg')
        options.nodeWeightSecondary = ones(nNodes,1)*15;
    end
end
%Node color
if ~isfield(options,'nodeColor')
    options.nodeColor = repmat([255,160,128],nNodes,1);
end
%Edge color
if ~isfield(options,'edgeColor')
    options.edgeColor = repmat([0,191,255],nEdges,1);
end
%Arrowhead color
if ~isfield(options,'edgeArrowColor')
    options.edgeArrowColor = repmat([0,0,255],nEdges,1);
end
%Edge thickness
if ~isfield(options,'edgeWeight')
    options.edgeWeight = ones(nEdges,1)*2;
    if strcmp(CB_MAP_OUTPUT,'svg')
        options.edgeWeight = ones(nEdges,1)*4;
    end
end
%Font Size
if ~isfield(options,'textSize')
    options.textSize = ones(max(nNodes,nEdges),1)*12;
    if strcmp(CB_MAP_OUTPUT,'svg')
        options.textSize = ones(max(nNodes,nEdges),1)*10;
    end
end
%Font Color
if ~isfield(options,'textColor')
    options.textColor = zeros(nNodes,3);
end

% if ~isfield(options,'otherTextSize')
%     options.otherTextSize = ones(size(map.text,1),1)*12;
%     if strcmp(CB_MAP_OUTPUT,'svg')
%         options.otherTextSize = ones(size(map.text,1),1)*135;
%     end
% end

if ~isfield(options,'otherTextColor')
    options.otherTextColor= zeros(size(map.text,1),3);
end

nodeWeight = options.nodeWeightSecondary;
nodeWeight(strcmp(map.molPrime,'Y')) = options.nodeWeight(strcmp(map.molPrime,'Y'));

if ~isfield(options,'fileName')
    options.fileName = 'target.svg';
end

if ~isfield(options,'rxnDir')
    options.rxnDir = zeros(size(map.connectionAbb,1),1);
end

if ~isfield(options,'rxnDirMultiplier')
    options.rxnDirMultiplier = 2;
end

%%%%%%%% initialization
if strcmp(CB_MAP_OUTPUT,'matlab')    % use matlab to draw the map
    clf; % this was in line 41 before
    % setting the color bar
    figure(1);colormap(cool(100))
    colorbar('location','southoutside');
    axis equal;
    hold on
elseif strcmp(CB_MAP_OUTPUT, 'java')
    % use Java/OpenGL to draw the map
    plotcbmap;
    % send the transformation coordinates
    R=map.molPosition';
    R=sort(R);
    a=max(R);
    b=min(R);
    xmax=a(1,1);
    xmin=b(1,1);
    ymax=a(1,2);
    ymin=b(1,2);
    settrans(mapHandle,xmax,xmin,ymax,ymin);
elseif strcmp(CB_MAP_OUTPUT, 'svg')
    %check fileName extension
    if isempty(regexp(lower(options.fileName),'.svg$'))
        options.fileName = strcat(options.fileName,'.svg');
    end
    textPos = (map.textPos);
    x1 = min(map.molPosition(1,:));
    if min(textPos(:,1))<x1
        x1 = min(textPos(:,1));
    end
    y1 = min(map.molPosition(2,:));
    if min(textPos(:,2))<y1
        y1 = min(textPos(:,2));
    end
    x2 = max(map.molPosition(1,:));
    if max(textPos(:,1))>x2
        x2 = max(textPos(:,1));
    end
    y2 = max(map.molPosition(2,:));
    if max(textPos(:,2))>y2
        y2 = max(textPos(:,2));
    end
    if isfield(options,'colorScale')
        numColorBins = size(options.colorScale,1);
        colorScaleWidth = 0.25*(x2-x1);
        binWidth = colorScaleWidth/numColorBins;
        colorScaleHeight = 0.05*colorScaleWidth;
        colorScaleBuffer = colorScaleHeight*1.5;
        y2 = y2+colorScaleBuffer; %add buffer for scale
    end
    if isfield(options,'fluxVarColor')
        colorScaleHeight = 0.0125*(x2-x1);
        colorScaleBuffer = colorScaleHeight*1.5;
        colorWidth = 0.025*(x2-x1);
        y2 = y2+colorScaleBuffer;
    end
    [x1,y1,x2,y2] = deal(x1-200, y1-200, x2+200, y2+200); % add buffer
    SF = .25;
    mapHandle = fopen(options.fileName, 'w');
    fprintf(mapHandle, '<?xml version="1.0" encoding="UTF-8" standalone="yes"?>\n');
    fprintf(mapHandle,'<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.0//EN" "http://www.w3.org/TR/2001/REC-SVG-20010904/DTD/svg10.dtd">\n');
    fprintf(mapHandle,'<svg height="%+.2f" width="%+.2f" viewBox="%+.2f %+.2f %+.2f %+.2f" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">\n',(y2-y1+20)*SF,(x2-x1+20)*SF,(x1-10)*SF,(y1-10)*SF,(x2-x1+0)*SF,(y2-y1+20)*SF);
    fprintf(mapHandle,'<g transform="scale(%5.3f)">\n',SF);

    %Add Scale
    if isfield(options,'overlayType')
        for i=1:numColorBins
            color = strcat('rgb(',num2str(options.colorScale(i,1)),',',num2str(options.colorScale(i,2)),',',num2str(options.colorScale(i,3)),')');
            fprintf(mapHandle,'<g style="stroke-linecap: round; stroke-linejoin: round; stroke-miterlimit: 20; fill: %s; ">\n',color);
            fprintf(mapHandle,'<rect x="%.2f" y="%.2f" width="%.2f" height="%.2f" />\n</g>\n',x2-colorScaleWidth-200*2+(binWidth*(i-1)),y2-colorScaleHeight*1.25,binWidth,colorScaleHeight);
        end
        fprintf(mapHandle, '<g style="font-family: sans; stroke: none; text-anchor:end">\n');
        fprintf(mapHandle,'<text style="fill: rgb(0,0,0); text-rendering: optimizeLegibility;" x="%.2f" y="%.2f" font-size="%dpx">%s</text>\n</g>\n',x2-colorScaleWidth-200*2-0.25*colorScaleHeight,y2-colorScaleHeight*.4,colorScaleHeight,['Scale: ' options.scaleTypeLabel '    ' options.overlayType ': ' options.overlayLB]);
        fprintf(mapHandle, '<g style="font-family: sans; stroke: none; text-anchor:start">\n');
        fprintf(mapHandle,'<text style="fill: rgb(0,0,0); text-rendering: optimizeLegibility;" x="%.2f" y="%.2f" font-size="%dpx">%s</text>\n</g>\n',x2-200*2+0.25*colorScaleHeight,y2-colorScaleHeight*.4,colorScaleHeight,[options.overlayUB]);
    end
    if isfield(options,'fluxVarColor')
        colorTextLabel = {'Bidirectional / reversible:', 'Unidirectional / reversible forward:', 'Unidirectional / reversible reverse:', 'Unidirectional / irreversible:'};
        color{1} = strcat('rgb(',num2str(options.fluxVarColor.biDirColor(1)),',',num2str(options.fluxVarColor.biDirColor(2)),',',num2str(options.fluxVarColor.biDirColor(3)),')');
        color{2} = strcat('rgb(',num2str(options.fluxVarColor.uniDirFwdColor(1)),',',num2str(options.fluxVarColor.uniDirFwdColor(2)),',',num2str(options.fluxVarColor.uniDirFwdColor(3)),')');
        color{3} = strcat('rgb(',num2str(options.fluxVarColor.uniDirRevColor(1)),',',num2str(options.fluxVarColor.uniDirRevColor(2)),',',num2str(options.fluxVarColor.uniDirRevColor(3)),')');
        color{4} = strcat('rgb(',num2str(options.fluxVarColor.uniDirIrrColor(1)),',',num2str(options.fluxVarColor.uniDirIrrColor(2)),',',num2str(options.fluxVarColor.uniDirIrrColor(3)),')');
        for i=2:3:11
            fprintf(mapHandle, '<g style="font-family: sans; stroke: none; text-anchor:end">\n');
            fprintf(mapHandle,'<text style="fill: rgb(0,0,0); text-rendering: optimizeLegibility;" x="%.2f" y="%.2f" font-size="%dpx">%s</text>\n</g>\n',i*(x2-x1)/12,y2-colorScaleHeight*.4,colorScaleHeight,colorTextLabel{ceil(i/3)});
            fprintf(mapHandle,'<g style="stroke-linecap: round; stroke-linejoin: round; stroke-miterlimit: 20; fill: %s; ">\n',color{ceil(i/3)});
            fprintf(mapHandle,'<rect x="%.2f" y="%.2f" width="%.2f" height="%.2f" />\n</g>\n',i*(x2-x1)/12,y2-colorScaleHeight*1.25,colorWidth,colorScaleHeight);
        end
    end
end

%%%%% actual map drawing code
% draw other shapes
if isfield(map,'shapeThickness')
    for i = 1:size((map.shapeThickness),1)
        drawShape(map.shapeType(i,1),map.shapePos(i,1:2),map.shapeSize(i,1:2),map.shapeColor(i,1:3),map.shapeThickness(i,1),map.shapeStyle(i,1));
    end
end
% draw the connection segments traversing through the connection matrix
for i = 1:(size((map.connection),1))
    drawLine(map.connection(i,1),map.connection(i,2),map,options.edgeColor(i,:),options.edgeArrowColor(i,:),options.edgeWeight(i),nodeWeight,options.rxnDir(i),options.rxnDirMultiplier);
end
% draw the circles representing molecules
for i = 1:size((map.molPosition),2)
    drawCircle(map.molPosition(:,i),nodeWeight(i),options.nodeColor(i,:));
end
% draw texts
for i = 1:length(map.text)
    textFont =map.textFont{i};
    if regexp(textFont,'@')
        [textFont, textSize] = strtok(textFont,'@');
        textSize = str2num(regexprep(textSize,'@',''));
    elseif(map.textSize(i) >= 60)
        textSize = 60;
    else
        textSize = map.textSize(i);
    end
    if find(regexp(textFont,'Italic'))
        textStyle = 'italic;';
    else
        textStyle = '';
    end
    %textFont
    % if find(regexp(textFont,' B')) || find(regexp(textFont(end),'B'))
    if (find(regexp(textFont,'B')))
        textWeight = 'bold';
        textFont = regexprep(textFont,' B','');
    else
        textWeight = '';
    end
    if isfield(options,'otherTextSize'), textSize = options.otherTextSize(i); end
    drawText(map.textPos(i,1),map.textPos(i,2),map.text{i,1},textSize,textStyle,options.otherTextColor(i,:),lower(textFont),textWeight,true);
end
% Write Metabolite Label
for i = 1:size((map.molPosition),2)
    % write the labels for molecules
    if(options.textSize(i) ~= 0)
        drawText(map.molLabelPos(i,1),map.molLabelPos(i,2),map.molAbbreviation{i},options.textSize(i),'',options.textColor(i,:));
    end
end
% Write Reaction Label
for i = 1:size(map.rxnLabelPosition,2)
    if ~any(isnan(map.rxnLabelPosition(:,i)))
      if isfield(options, 'rxnTextSize')
        drawText(map.rxnLabelPosition(1,i),map.rxnLabelPosition(2,i),map.connectionAbb{find(map.rxnIndex(i)==map.connection,1)},options.rxnTextSize(i),'italic');
      else
        drawText(map.rxnLabelPosition(1,i),map.rxnLabelPosition(2,i),map.connectionAbb{find(map.rxnIndex(i)==map.connection,1)},10,'italic');
      end
    end
end

if strcmp(CB_MAP_OUTPUT,'matlab')
    hold off;
elseif strcmp(CB_MAP_OUTPUT,'java')

elseif strcmp(CB_MAP_OUTPUT,'svg')

    fprintf(mapHandle,'</g>\n');
    fprintf(mapHandle,'</svg>\n');
    fclose(mapHandle);
    display('Document Written')
end


text:
%Script to run Dose-recomputation using QIB
%
%APA, 03/27/2009
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC
indexS = planC{end};

% --------------------------------- Option to Downsample CT
ButtonName = questdlg('Do you wish to downsample CT by factor of 2?');
if strcmpi(ButtonName,'yes')
    planC = getplanCDownSample(planC, planC{indexS.CERROptions}, 2);
end

% --------------------------------- Target Struture used to Scale dose (D98)
ptvStructNum = 13;

% --------------------------------- Clinical Dose used to Scale dose
clinicDoseNum = 1;

% --------------------------------- Name IMRTP
imrtpName = 'Leaf Seq 2';

%% --------------------------------- Inputs for QIB
MCsolver = 3;
scatterThreshold = 0.1;
structNumsV = [5 6 7]; %Structure indices where dose needs to be calculated.
sampleRateV = [4 2 2];   %NOTE: Skin can be sampled at rate of 8, whereas targets and critical
%                          structures can be sampled at rate of 2. This helps speed-up calculation
%                          while maintaining dose accuracy for important structures. Same length as structNumsV

%% --------------------------------- Inputs for DPM
% MCsolver = 1; %1:DPM, 2:VMC++
% scatterThreshold = 11.2; %DPM
% structNumsV = []; %leave structure number empty since DPM computes dose over entire body
% sampleRateV = []; %leave empty for DPM calculation

%Change nhist and batch
nhist = 1e5; %per cm^2, use 1000000 for testing; %1M is 2% uncertainty, 4M is 1%
batch = 101; %must be unique for each calculation
planC_File = ''; %Leave empty since it is global.
leak = 0.032;
spectrum_File = '6MV10x10MDA.spectrum';
OutputError = 0; %not working for 1
whichBeam = 1;  %redefined later
PBMaxWidth = 10; %
gradsense = 25;
% MCsolver = 1; %1:DPM, 2:VMC++
saveIM = 0;
sourceModel = 0; %1 to use source model
doseToWater = 0;
fillWater = 0;
useWedge = 0;
inputPB = 0;
inputIM = 1;
Softening = 1;
UseFlatFilter = 1;
MLC = 0; %1?
TongueGroove = 0;
interactiveMode = 1;
LS_flag = 0;
K = 1;
% scatterThreshold = 11.2; %DPM
% structNumsV = []; %leave structure number empty since DPM computes dose over entire body
% sampleRateV = []; %leave empty for DPM calculation

% % --------------------------------- Inputs for QIB
% MCsolver = 3;
% scatterThreshold = 0.1;
% structNumsV = [2 1 5 6 11]; %Structure indices where dose needs to be calculated.
% sampleRateV = [8 2 2 2 2]; %NOTE: Skin can be sampled at rate of 8, whereas targets and critical
% %                           structures can be sampled at rate of 2. This helps speed-up calculation
% %                           while maintaining dose accuracy for important structures
%
% structNumsV = [11 12 9 6];
% sampleRateV = [2 2 2 2]; %same length as structNumsV

%% --------------------------------- Call dose calculation engine
try
    numBeams = planC{indexS.beams}(1).FractionGroupSequence.Item_1.NumberofBeams;
catch
    numBeams = planC{indexS.beams}(1).FractionGroupSequence.Item_1.NumberOfBeams;
end
beamletWeightC = {};
for whichBeam = 1:numBeams
    %IMwDPM = beam2MCdose(leak, spectrum_File, planC_File, nhist, OutputError, whichBeam, PBMaxWidth, gradsense, MCsolver, saveIM, sourceModel, doseToWater, fillWater, useWedge, inputPB, inputIM, Softening, UseFlatFilter, MLC, TongueGroove, batch, interactiveMode, LS_flag, K, scatterThreshold);
    [IMwQIB(whichBeam), beamletWeightC{whichBeam}] = beam2MCdose_with_QIB(leak, spectrum_File, planC_File, nhist, OutputError, whichBeam, PBMaxWidth, gradsense, MCsolver, saveIM, sourceModel, doseToWater, fillWater, useWedge, inputPB, inputIM, Softening, UseFlatFilter, MLC, TongueGroove, batch, interactiveMode, LS_flag, K, scatterThreshold, structNumsV, sampleRateV);
end

if MCsolver == 1 %DPM
    % --------------------------------- Read dose from disk and apply Meter-Set
    [dose3Dsum] = calcDoseByBeamMeterset(planC, nhist, batch);
    doseName = 'DPM reCalc';

elseif MCsolver == 3 %QIB

    planC{indexS.IM}(1).IMDosimetry.beams = [IMwQIB.beams];
    planC{indexS.IM}.IMDosimetry.goals = IMwQIB(1).goals;
    planC{indexS.IM}.IMDosimetry.params = IMwQIB(1).params;
    clear IMwQIB
    planC{indexS.IM}.IMDosimetry.isFresh = 1;
    planC{indexS.IM}.IMDosimetry.name = imrtpName;
    planC{indexS.IM}.IMDosimetry.solutions = [];
    planC{indexS.IM}.IMDosimetry.assocScanUID = planC{indexS.scan}(1).scanUID;
    planC{indexS.IM}.IMDosimetry.solutions = [beamletWeightC{:}];
    [dose3D] = getIMDose(planC{indexS.IM}(end).IMDosimetry, [beamletWeightC{:}], structNumsV);
    doseName = 'QIB reCalc';
end

%% --------------------------------- Scale to clinical dose
doseNum = length(planC{indexS.dose});

%Scale based on PTV D98
reCalc_metric   = Dx(planC, ptvStructNum, doseNum, 98);
Clinical_metric = Dx(planC, ptvStructNum, clinicDoseNum, 98);

%Scale based on PTV meanDose
reCalc_metric   = meanDose(planC, ptvStructNum, doseNum, 'Absolute');
Clinical_metric = meanDose(planC, ptvStructNum, clinicDoseNum, 'Absolute');

planC{indexS.dose}(doseNum).doseArray = planC{indexS.dose}(doseNum).doseArray*Clinical_metric/reCalc_metric;

%Show dose in CERR
showIMDose(dose3D,doseName)


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

%Fourier Transfrom and  inverse Fourier Transfrom  of various functions
% the order of the operations is a little different
%from the order that they appear on the book

%F{exp(-t^2)}
syms t w
x=exp(-t^2);
fourier(x)

int(x*exp(-j*w*t),t,-inf,inf)

%F^-1{1/(1+j*w)}
X=1/(1+j*w);
ifourier(X)

X=1/(1+j*w);
ifourier(X,t)

syms n
X=1/(1+j*w);
ifourier(X,n)

%F{exp(-t)*u(t)}
syms t w
x=exp(-t)*heaviside(t);
X=fourier(x,w)

%F{1}
x=1;
fourier(x,w)
syms s
fourier(x,s)
fourier(x)


text:
function t = ktensor(varargin)
%KTENSOR Tensor stored as a Kruskal operator (decomposed).
%
%   K = KTENSOR(lambda,U1,U2,...,UM) creates a Kruskal tensor from its
%   constituent parts. Here lambda is a k-vector and each Um is a
%   matrix with k columns.
%
%   K = KTENSOR(lambda, U) is the same as above except that U is a
%   cell array containing matrix Um in cell m.
%
%   K = KTENSOR(U) assumes U is a cell array containing matrix Um in
%   cell m and assigns the weight of each factor to be one.
%
%   K = KTENSOR(T) creates a ktensor by copying an existing ktensor.
%
%   Examples
%   K = ktensor([3; 2], rand(4,2), rand(5,2), rand(3,2))
%
%   See also TENSOR, TTENSOR, KTENSOR/FULL.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

% EMPTY CONSTRUCTOR
if nargin == 0
    t.lambda = [];
    t.u = {};
    t = class(t,'ktensor');
    return;
end

% Copy CONSTRUCTOR
if (nargin == 1) && isa(varargin{1}, 'ktensor')
    t.lambda = varargin{1}.lambda;
    t.u = varargin{1}.u;
    t = class(t, 'ktensor');
    return;
end

if isa(varargin{1},'cell')

    u = varargin{1};
    t.lambda = ones(size(u{1},2),1);
    t.u = u;

else

    t.lambda = varargin{1};
    if ~isa(t.lambda,'numeric') || ndims(t.lambda) ~=2 || size(t.lambda,2) ~= 1
	error('LAMBDA must be a column vector.');
    end

    if isa(varargin{2},'cell')
	t.u = varargin{2};
    else
	for i = 2 : nargin
	    t.u{i-1} = varargin{i};
	end
    end

end

% Check that each Um is indeed a matrix
for i = 1 : length(t.u)
    if ndims(t.u{i}) ~= 2
	error(['Matrix U' int2str(i) ' is not a matrix!']);
    end
end

% Size error checking
k = length(t.lambda);
for i = 1 : length(t.u)
    if  size(t.u{i},2) ~= k
       error(['Matrix U' int2str(i) ' does not have ' int2str(k) ' columns.']);
    end
end

t = class(t, 'ktensor');
return;